!################################################################# ! ! Purpose: ! -------- ! ! This module contains all subroutines needed for dry deposition calculations. ! The mean purpose is to perform on a high resolution of 1x1 degree: ! ! 0. allocate the vd on the model resolution ! ! 1. read fields needed for further calculations in subroutines: ! ! 2. calculate the tracer independent friction velocity, aerodynamic ! and sub-laminar resistance ! ! 3. calculate fields needed for resistance calculations ! ! 4. the tracer dependent vd ! ! 5. coarsen the vd ! ! 6. apply the coarsened deposition velocities to concentration field rm ! ! 7. deallocate vd ! ! Reference: ! --------- ! general documentationon ECMWF fields can be found on: ! http://www.ecmwf.int/research/ifsdocs/PHYSICS/ ! --------- ! Authors: ! ------- ! original code by Laurens Ganzeveld and Ad Jeuken (1997) ! revised and adapted for use in TM5 and use of ECMWF data ! by F. Dentener and M. Krol (2003) ! ! 24 Jan 2012 - P. Le Sager - modified for mpi lat-lon domain decomposition ! !### macro's ##################################################### ! #define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr #define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if #define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if ! #include "tm5.inc" ! !################################################################# module dry_deposition use GO , only : gol, goPr, goErr use dims , only : nregions use chem_param , only : ntrace use global_types, only : emis_data implicit none ! --- in/out ---------------------- private public :: Dry_Deposition_Init, Dry_Deposition_Done public :: dd_surface_fields public :: vd public :: vda ! --- const ----------------------- character(len=*), parameter :: mname = 'dry_deposition' ! --- var ------------------------- ! vd : final deposition velocities for all species on model resolution. ! vd_temp : to store the vd temporarily for one specie type(emis_data),dimension(nregions,ntrace) :: vd type(emis_data),dimension(nregions) :: vd_temp, vd_temp_global ! vda : final deposition velocities for all aerosols on model resolution. type(emis_data), pointer :: vda(:,:) contains ! ! Allocate emission data structure "vd", ! final deposition velocities for all species on model resolution ! subroutine Dry_Deposition_Init( status ) use dims , only : nregions, iglbsfc use chem_param, only : ndep use chem_param, only : nrdep use meteodata , only : lsmask_dat use meteodata , only : sr_ecm_dat, sr_ols_dat, ci_dat, sd_dat, swvl1_dat use meteodata , only : nsss_dat, ewss_dat, wspd_dat, slhf_dat, sshf_dat use meteodata , only : src_dat, d2m_dat, t2m_dat, ssr_dat use meteodata , only : nveg use meteodata , only : tv_dat, cvl_dat, cvh_dat use meteodata , only : Set use tm5_distgrid, only : dgrid, Get_DistGrid use ParTools, only : isRoot ! --- in/out ------------------------------------ integer, intent(out) :: status ! --- const ------------------------------------- character(len=*), parameter :: rname = mname//'/Dry_Deposition_Init' ! --- local ------------------------------------- integer :: imr,jmr,region,n integer :: iveg, i1, i2, j1, j2 ! --- begin ------------------------------------- ! deposition velocities ; ! always allocated, filled with zeros if no tracers are to be deposited: do region = 1, nregions CALL GET_DISTGRID( DGRID(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) do n=1,ntrace allocate(vd(region,n)%surf(i1:i2, j1:j2)) vd(region,n)%surf = 0.0 end do end do ! deposited tracers ? if ( ndep > 0 ) then ! select meteo to be used call Set( lsmask_dat(iglbsfc), status, used=.true. ) call Set( sr_ecm_dat(iglbsfc), status, used=.true. ) call Set( sr_ols_dat(iglbsfc), status, used=.true. ) call Set( ci_dat(iglbsfc), status, used=.true. ) call Set( sd_dat(iglbsfc), status, used=.true. ) call Set( swvl1_dat(iglbsfc), status, used=.true. ) call Set( ewss_dat(iglbsfc), status, used=.true. ) call Set( nsss_dat(iglbsfc), status, used=.true. ) call Set( wspd_dat(iglbsfc), status, used=.true. ) call Set( slhf_dat(iglbsfc), status, used=.true. ) call Set( sshf_dat(iglbsfc), status, used=.true. ) call Set( src_dat(iglbsfc), status, used=.true. ) call Set( d2m_dat(iglbsfc), status, used=.true. ) call Set( t2m_dat(iglbsfc), status, used=.true. ) call Set( ssr_dat(iglbsfc), status, used=.true. ) call Set( cvl_dat(iglbsfc), status, used=.true. ) call Set( cvh_dat(iglbsfc), status, used=.true. ) do iveg = 1, nveg call Set( tv_dat(iglbsfc,iveg), status, used=.true. ) end do ! temporary storage: do region = 1, nregions CALL GET_DISTGRID( DGRID(region), I_STRT=i1, I_STOP=i2, & J_STRT=j1, J_STOP=j2, & I_WRLD=imr, J_WRLD=jmr ) allocate(vd_temp(region)%surf(i1:i2,j1:j2)) if(isRoot) then allocate(vd_temp_global(region)%surf(imr,jmr)) else allocate(vd_temp_global(region)%surf(1,1)) end if end do end if ! deposited tracers ! aerosols ? if ( nrdep > 0 ) then allocate( vda(nregions,nrdep) ) do region = 1, nregions CALL GET_DISTGRID( DGRID(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) do n=1,nrdep allocate(vda(region,n)%surf(i1:i2,j1:j2)) vda(region,n)%surf = 0.0 end do end do end if ! ok status = 0 end subroutine Dry_Deposition_Init ! *** ! ! De-allocate emission data structure "vd", ! final deposition velocities for all species on model resolution ! subroutine Dry_Deposition_Done( status ) use dims , only : nregions use chem_param, only : ndep use chem_param, only : nrdep ! --- in/out ------------------------------------ integer, intent(out) :: status ! --- const ------------------------------------- character(len=*), parameter :: rname = mname//'/Dry_Deposition_Done' ! --- local ------------------------------------- integer :: region, n ! --- begin ------------------------------------- ! clear deposition velocities: do region = 1, nregions do n=1,ntrace deallocate(vd(region,n)%surf) end do end do ! clear temporary storage: if ( ndep > 0 ) then do region = 1, nregions deallocate(vd_temp(region)%surf) deallocate(vd_temp_global(region)%surf) end do end if ! deposited tracers ! aerosols ? if ( nrdep > 0 ) then do region = 1, nregions do n = 1, nrdep deallocate( vda(region,n)%surf ) end do end do deallocate( vda ) end if ! ok status = 0 end subroutine Dry_Deposition_Done !------------------------------------ ! ! Purpose: ! ------- ! this subroutine reads and prepares all surface fields ! ! External ! -------- ! dd_get_3_hourly_surface_fields ! dd_calc_ustar_raero_rb ! dd_calc_rstom_rahcan ! dd_calc_inisurf ! dd_calc_rs ! dd_coarsen_vd ! ! Reference ! --------- ! Ganzeveld and Lelieveld (1996) and references therein. ! !------------------------------------ subroutine dd_surface_fields( status ) use binas , only : vkarman use dims , only : idate use dims , only : okdebug use chem_param , only : names use chem_param , only : ndep use chem_param , only : nrdep #ifndef without_photolysis use photolysis_data, only : phot_dat #endif use ParTools, only : isRoot use dims , only : iglbsfc, nlat180, nlon360 use meteodata , only : lsmask_dat use meteodata , only : sr_ecm_dat, sr_ols_dat, ci_dat, sd_dat, swvl1_dat use meteodata , only : ewss_dat, nsss_dat, wspd_dat, slhf_dat, sshf_dat use meteodata , only : src_dat, d2m_dat, t2m_dat, ssr_dat use tm5_distgrid, only : dgrid, Get_DistGrid, scatter, gather, dist_arr_stat ! --- in/out ----------------------------- integer, intent(out) :: status ! --- const ------------------------------- character(len=*), parameter :: rname = 'dd_surface_fields' ! --- local ------------------------------ character(len=10) :: c_time real, dimension(:,:,:), allocatable :: soilph real, dimension(:,:), allocatable :: lai, lai1 real, dimension(:,:), allocatable :: vgrat_low, vgrat_high real, dimension(:,:), allocatable :: sstr, wind10m real, dimension(:,:), allocatable :: ustar_loc,rb,raero real, dimension(:,:,:), allocatable :: vd11 real, dimension(:,:), allocatable :: vd11a real, dimension(:,:), allocatable :: rahcan real, dimension(:,:), allocatable :: rstom real, dimension(:,:), allocatable :: snow_cover real, dimension(:,:), allocatable :: wet_skin real, dimension(:,:), allocatable :: fws real, dimension(:,:), allocatable :: rh2m real, dimension(:,:), allocatable :: ddepvel_h2 #ifndef without_photolysis real, dimension(:,:), allocatable :: ags #endif real, dimension(:,:), allocatable :: alpha, bubble, alphae real, dimension(:,:), allocatable :: ustar_land, ustar_sea, sr_sea real, dimension(:,:), allocatable :: global_field integer :: i, nsend, region,itrace integer :: i1, i2, j1, j2, imr, jmr ! no deposited tracers or aerosols ? then leave immediately: if ( ndep+nrdep == 0 ) then status=0; return end if ! From here until the remapping (coarsen), we work on the extra iglbsfc ! grid, on which the data to read are expected to be. CALL GET_DISTGRID( DGRID(iglbsfc), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2, I_WRLD=imr, J_WRLD=jmr ) if (isRoot) then allocate(global_field(imr,jmr)) else allocate(global_field(1,1)) end if allocate(soilph (i1:i2, j1:j2,5)) allocate(lai (i1:i2, j1:j2) ) allocate(lai1 (i1:i2, j1:j2) ) allocate(vgrat_low (i1:i2, j1:j2) ) allocate(vgrat_high(i1:i2, j1:j2) ) allocate(sstr (i1:i2, j1:j2) ) allocate(wind10m (i1:i2, j1:j2) ) allocate(ustar_loc (i1:i2, j1:j2) ) allocate(raero (i1:i2, j1:j2) ) allocate(rb (i1:i2, j1:j2) ) allocate(rahcan (i1:i2, j1:j2) ) allocate(rstom (i1:i2, j1:j2) ) allocate(snow_cover(i1:i2, j1:j2) ) allocate(wet_skin (i1:i2, j1:j2) ) allocate(fws (i1:i2, j1:j2) ) allocate(rh2m (i1:i2, j1:j2) ) allocate(ddepvel_h2(i1:i2, j1:j2) ) #ifndef without_photolysis allocate(ags (i1:i2, j1:j2) ) #endif allocate(vd11 (i1:i2, j1:j2,ndep)) if ( nrdep > 0 ) allocate(vd11a(i1:i2, j1:j2)) allocate(ustar_land(i1:i2, j1:j2)) allocate(ustar_sea (i1:i2, j1:j2)) allocate(sr_sea (i1:i2, j1:j2)) allocate(alpha (i1:i2, j1:j2)) allocate(bubble (i1:i2, j1:j2)) allocate(alphae (i1:i2, j1:j2)) ! -- for output of time header on debug fields write(c_time,'(i4,3(i2))') idate(1:4) !CMK corrected do i=1,10 if (c_time(i:i)==' ') c_time(i:i)='0' end do ! ! surface data sets may have the following characteristics ! 1. instanteneous. ! The time written in the file is valid from t-1.5 to t+1.5 ! 2. Accumulated. ! The time written in the file is valid from t to t+3 ! 3. Daily averaged. ! The time on the file is always 00 hours, and valid until t+24 ! ! *** It is essential to understand this timing when ! reading and applying these sets. ! ! fd2mk it has to be decided to 'save' fields like ! vgrat_low and daily average surface fields ! for later use, or to read it every 3 hours time step ! ! Read data call dd_get_soilph_lai ! non-ECMWF provided data: soilph and lai IF_NOTOK_RETURN(status=1) call dd_get_surface_vegetation ! vgrat_low, vgrat_high and lsm IF_NOTOK_RETURN(status=1) wind10m = wspd_dat(iglbsfc)%data(:,:,1) sstr = sqrt( ewss_dat(iglbsfc)%data(:,:,1)**2 + nsss_dat(iglbsfc)%data(:,:,1)**2 ) ! AS: ensure being non-zero where( sstr .lt. 1.E-10 ) sstr = 1.E-10 call dd_calc_ustar_raero_rb ! raero, ustar_loc, rb IF_NOTOK_RETURN(status=1) call dd_calc_rstom_rahcan ! rstom, rahcan call dd_calc_inisurf ! snow_cover, wet_skin, fws, rh2m, ags call dd_calc_rs ! vd IF_NOTOK_RETURN(status=1) ! coarsen vd to model resolution #ifndef without_photolysis call dd_coarsen_vd( vd11, ags, i1, j1, status ) #else call dd_coarsen_vd( vd11, i1, j1, status ) #endif IF_NOTOK_RETURN(status=1) deallocate(soilph) deallocate(lai) deallocate(lai1) deallocate(vgrat_low) deallocate(vgrat_high) deallocate(sstr) deallocate(wind10m) deallocate(ustar_loc) deallocate(raero) deallocate(rb) deallocate(rahcan) deallocate(rstom) deallocate(snow_cover) deallocate(wet_skin) deallocate(fws) deallocate(rh2m) deallocate(ddepvel_h2) #ifndef without_photolysis deallocate(ags) #endif deallocate(vd11) if ( nrdep > 0 ) deallocate( vd11a ) deallocate(alpha) deallocate(bubble) deallocate(alphae) deallocate(ustar_land) deallocate(ustar_sea) deallocate(sr_sea) !FD23012004 deallocate(global_field) ! PLS #ifdef MPI ! PLS do region = 1, nregions ! PLS nsend = im(region)*jm(region) ! PLS do itrace = 1, ntrace ! PLS call mpi_bcast(vd(region,itrace)%surf, nsend, my_real, root_k, & ! PLS com_lev, ierr) ! PLS end do ! PLS ! aerosols ? ! PLS if ( nrdep > 0 ) then ! PLS do itrace = 1, nrdep ! PLS call mpi_bcast( vda(region,itrace)%surf, nsend, MY_REAL, root_k, com_lev, ierr ) ! PLS end do ! PLS end if ! PLS #ifndef without_photolysis ! PLS ! send albedo computed in dd_calc_inisurf to all other processors: ! PLS call mpi_bcast(phot_dat(region)%albedo ,nsend, my_real, root_k, com_lev, ierr) ! PLS #endif ! PLS end do ! PLS if ( myid == root_k .and. okdebug ) then ! PLS write (gol,'("dd_surface_fields: Deposition velocities broadcasted")') ; call goPr ! PLS end if ! PLS #endif ! ok status = 0 contains !--------------------------- ! ! Purpose ! ------- ! Reads two independent datasets on soils and vegetation not provided by ECMWF ! ! Reference: ! --------- ! Data provided by Laurens Ganzeveld ! !-------------------------- SUBROUTINE DD_GET_SOILPH_LAI use global_data, only : inputdir USE MDF, ONLY : MDF_Open, MDF_NETCDF, MDF_READ, MDF_Inq_VarID, MDF_Get_Var, MDF_Close character(len=5) :: vname character(len=8) :: vname8 character(len=9) :: vname9 character(len=5),dimension(5) :: name_soil = (/'spec1','spec2','spec3','spec4','spec5'/) integer :: i, mo, n, varid, fid real, dimension(:,:), allocatable :: pland ! work array if(isRoot)then allocate( pland(nlon360,nlat180) ) else allocate( pland(1,1) ) end if ! ! -- Reading of input file with LAI data which is derived ! from the Olson ecosystem database (0.5 X 0.5 degrees), which ! discerns 46 ecosystems and their characteristics, e.g. LAI. ! Fields are monthly averaged. ! mo = idate(2) - 1 if (isRoot)then CALL MDF_Open( TRIM(inputdir)//'/land/lsmlai.nc4', MDF_NETCDF, MDF_READ, fid, status ) IF_NOTOK_RETURN(status=1) endif if (isRoot) then if ( mo < 10 ) then write(vname8,'("spec1__",i1)')mo CALL MDF_Inq_VarID( fid, vname8, varid, status ) IF_NOTOK_RETURN(status=1) else write(vname9,'("spec1__",i2)')mo CALL MDF_Inq_VarID( fid, vname9, varid, status ) IF_NOTOK_RETURN(status=1) endif CALL MDF_Get_Var( fid, varid, pland, status ) IF_NOTOK_RETURN(status=1) endif call scatter( dgrid(iglbsfc), lai, pland, 0, status) IF_NOTOK_RETURN(status=1) if (isRoot)then CALL MDF_Close( fid, status ) IF_NOTOK_RETURN(status=1) endif if ( okdebug ) then call dist_arr_stat(dgrid(iglbsfc), 'lai', lai, 0, status) IF_NOTOK_RETURN(status=1) end if ! ! -- Reading of input file with soil pH data,which are derived from ! a 0.5 x 0.5 degree soil pH database, Batjes, 1995 ! over land soilph 1 to 5 should add up to 1. ! SOILPH(I,J,1) - soil pH <5.5 ! SOILPH(I,J,2) - soil 5.5 0. ) then !>>> TvN ! Over land, the friction velocity ustar ! is calculated from the 10 m wind speed, u10. ! In IFS u10 is a diagnostic variable, calculated to be ! compatible with "open-terrain" wind speed observations. ! Over rough or inhomogeneous terrain (z0M > 0.03 m), ! it is calculated using an aerodynamic roughness ! length typical for open terrain with grassland (0.03 m) ! (see IFS cy31r1 documentation, p. 46). ! In the expressions applied in TM5, ! the stability profile functions Psi_M have disappeared, ! and only the logarithmic terms are kept. ! Apart from this, the expression for ustar was incorrect: ! 10./sr_land should be 75./sr_land. ! This has now been corrected. ! Moreover, over islands and near coast lines, ! zeros can occur in sr_ols_dat, ! which have to be removed. ! However, a lower bound of 1e-2 m = 1 cm ! seems too high for smooth surfaces, ! like deserts and ice caps. ! The minimum positive value in sr_ols_dat ! is 2.5e-4 m = 0.025 cm, ! which is obtained in parts of Antarctica. ! This is likely the result of regridding of ! a field with minimum value of 0.1 cm ! from 0.5x0.5 to 1x1 degrees. ! Please note that with the introduction of CY31R1, ! the orographic contribution to the aerodynamic roughness length ! in IFS has been replaced by an explicit specification of stress ! on model levels due to turbulent orographic form drag, ! and the climatological variable previously used ! has been replaced by a prognostic variable. ! In TM5, the climatological variable is still used, ! but it would be better to use the prognostic variable instead. ! Note also that the same calculation ! is done in diffustion.F90. !sr_land=max(1e-2,sr_ols_dat(iglbsfc)%data(i,j,1)) ! occurs at Islands, etc. sr_land=max(1e-3,sr_ols_dat(iglbsfc)%data(i,j,1)) sr_help=min(sr_ecm_dat(iglbsfc)%data(i,j,1),0.03) !fd ustar_land=vKarman*wind10m(i,j)/alog(10./sr_help) ! !ustar consistent with 'clipped' large scale roughness !ustar_land(i,j)=vKarman*wind10m(i,j)/alog(10./sr_help)*& ! alog(75./sr_help)/alog(10./sr_land) ustar_land(i,j)=vKarman*wind10m(i,j)/alog(10./sr_help)*& alog(75./sr_help)/alog(75./sr_land) ! <<< TvN else ustar_land(i,j)=0.0 sr_land=0. endif ustar_loc(i,j)=xland*ustar_land(i,j)+(1-xland)*ustar_sea(i,j) sr_mix(i,j) =xland*sr_land +(1-xland)*sr_sea(i,j) end do !i end do !j ! ! Calculation of quasi laminar resistance of the layer above the surface. ! E.g. Walcek, Atmos. Environ, 20, pp. 949-964, 1986, and earlier refs. ! do j=j1,j2 do i=i1,i2 !>>> TvN ! Instead of the current approach, ! it is proposed to use the surface roughness length ! for heat directly as it is calculated prognostically ! in IFS (GRIB number 245). !<<< TvN rb(i,j)=(rz0/(ustar_loc(i,j)*vkarman)) end do end do ! ! calculate the aerodynamic resistance ! ! ra,the aerodynamic resistence in s/m over a layer with height z is ! the integral from z0 to z of 1/K(z) dz with K(z)=k.u*.z/phi(h) ! phi(h)=0.74(1-9z/l)**-0.5 for l<0 and 0.74+4.7(z/l) for l>0. ! by integration we get ra=0.74/(k.u*).(ln(z0/z)+ghi(z0)-ghi(z)) ! ghi(z)=-6.4(z/l) for l>0 and 2.ln(y(z)+1) for l<0 with ! y(z)=(1-9.z/l)**0.5. do j=j1,j2 do i=i1,i2 buoy = -sshf_dat(iglbsfc)%data(i,j,1) / rhoCp & -0.61 * t2m_dat(iglbsfc)%data(i,j,1) * slhf_dat(iglbsfc)%data(i,j,1)/rhoLv tstv=-buoy/ustar_loc(i,j) obuk=1e-6 if( abs(tstv) .gt. tiny(tstv) ) then obuk = ustar_loc(i,j)*ustar_loc(i,j)*t2m_dat(iglbsfc)%data(i,j,1)/(tstv*grav*vKarman) end if if( obuk > 0. ) then ! stable conditions ra = 0.74*( alog(Href/sr_mix(i,j)) + 6.4*(Href-sr_mix(i,j))/obuk )& / (vKarman*ustar_loc(i,j)) else ! unstable y0 = sqrt(1.-9.*sr_mix(i,j)/obuk)+1. yra = sqrt(1.-9.*Href/obuk)+1. ra = 0.74*(alog(Href/sr_mix(i,j))+2.*(alog(y0)-alog(yra)))/ & (vKarman*ustar_loc(i,j)) end if raero(i,j)=max(10.,min(ra,1e10)) ! end do !i end do !j if ( okdebug ) then call dist_arr_stat(dgrid(iglbsfc), 'ustar_loc', ustar_loc, 0, status ) IF_NOTOK_RETURN(status=1) call dist_arr_stat(dgrid(iglbsfc), 'raero', raero , 0, status ) IF_NOTOK_RETURN(status=1) call dist_arr_stat(dgrid(iglbsfc), 'rb', rb , 0, status ) IF_NOTOK_RETURN(status=1) write(*,*) 'end calc_aero_ustar' end if deallocate(sr_mix) end subroutine dd_calc_ustar_raero_rb ! ! subroutine dd_calc_rstom_rahcan ! ----------------------------- ! Purpose ! --------- ! Calculate water vapour stomatal resistance from the PAR ! (Photosythetically Active Radiation) which is 0.55 times ! the net short wave radiation (ssr) at the surface. ! Calculate rahcan from ustar and lai ! ! External ! -------- ! none ! ! Reference ! ---------- ! none !-------------------------------------------------- ! ! -- constants used for the calculation of the stomatal resistance ! (see ECHAM3 manual) implicit none real,parameter :: a=5000. real,parameter :: b=10. real,parameter :: c=100. real,parameter :: vk=0.9 real,parameter :: vlt=1. ! real, parameter :: foresth=20. integer :: j, i real :: vpar, d real :: canht,laihelp ! ! -- recalculated rstom for a LAI of 1 (see ECHAM3 manual) ! do j=j1,j2 do i=i1,i2 ! ! -- calculation of PAR from net short wave radiation ! -- radiation is corrected for daily cycle ! rstom(i,j)=1e10 !high resistance during the night if (ssr_dat(iglbsfc)%data(i,j,1) > 0.) then vpar=max(1.,0.55*ssr_dat(iglbsfc)%data(i,j,1)) d=(a+b*c)/(c*vpar) rstom(i,j)=(vk*c)/((b/(d*vpar))*alog((d*exp(vk*vlt)+1.)/(d+1.))-& alog((d+exp(-vk*vlt))/(d+1.))) end if !ssr > 0. end do !i end do !j ! ! -- computation of in-canopy aerodynamic resistance from canopy height ! and the friction velocity and the Leaf Area Index ! (see Erisman & Van Pul) ! do j=j1,j2 do i=i1,i2 ! ! ! -- calculation of canopy height CANHT from effective fraction ! of high vegetation assuming that forest has a canopy height ! of 20 m (other vegetation 0 m) ! canht= vgrat_high(i,j)*foresth ! laihelp: the maximum values derived from ECMWF are more realistic laihelp=min(lai1(i,j),lai(i,j)) rahcan(i,j)=max(1.e-5,14.*laihelp*canht/ustar_loc(i,j)) end do !j end do !i !cmk if ( okdebug ) & !cmk call dumpfield(0,'rahcan'//c_time//'.hdf',rahcan,'rahcan') !cmk if ( okdebug ) call dumpfield(0,'rstom'//c_time//'.hdf',rstom,'rstom') if ( okdebug ) write(*,*) 'dd_calc_rstom_rahcan: end' ! end subroutine dd_calc_rstom_rahcan ! ! subroutine dd_calc_inisurf ! ------------------------ ! Purpose !--------- ! Calculate some surface fields later needed in the calculation ! ! External ! -------- ! none ! ! Reference ! ---------- ! none !-------------------------------------------------- use Binas, only: pi implicit none real,parameter :: sncr = 0.015 ! critical snow cover depth real,parameter :: tstar = 273. real,parameter :: wlmax = 2.e-4 real,parameter :: ewsmx = 0.323 real,parameter :: ewsmx_sat= 0.472 real,parameter :: wpwp = 0.171 !FD23012004 real,parameter :: eff=0.5 ! parameters needed for sea/bubble formation real,parameter :: rdrop=0.005 ! cm real,parameter :: zdrop=10.0 ! cm real,parameter :: eps=0.6 ! parameter related to bubble formation real :: qdrop,s,phi,alpha1,vk1,vk2,alpharat ! auxiliary parameters for bubble formation real :: sr_help ! surface roughnes consistent with 10 m wind. !FD23012004 ! ! for albedo calculations real :: snow, freesea, bare, desert real :: e,es,wcr,wlmx,xland,ahelp,vgr,laihelp integer :: i,j ! auxiliary variables ! --- begin --------------------------------------- do j=j1,j2 do i=i1,i2 xland = lsmask_dat(iglbsfc)%data(i,j,1)/100.0 ! land-sea fraction directly from ECMWF ! the maximum values derived from ECMWF are more realistic laihelp=min(lai1(i,j),lai(i,j)) ! ! -- calculation of monthly snow cover fraction sd using ! constant critical snow depth ! alternatively this critical snow depth may be chosen as ! a linear function of ln(sr_ecm) ! B. v.d. Hurk (2002) suggests 0.015 m (=sncr) for z0<0.25 ! and 1 m for z0>5m and log-linear in between in between. ! factor 0.9 is introduced to account for the fact that ! it is unlikely ! that 'high' vegetation is completely snow covered ! ! FD25.03.2003 snow_cover(i,j)=min(xland,sd_dat(iglbsfc)%data(i,j,1)/sncr) ! ! -- calculation of wet skin fraction from the skin reservoir ! content src (prognostic variable in ECMWF) ! the vegetation fractions and their attributed LAI, ! Wlmax which is the maximum amount of water that can be held ! on one layer of leaves or bare ground, Wlmx is the ! maximum skin reservoir content ! if ( xland > 0.0 ) then ! bare soil fraction small discrepancy when lakes are present bare=max(0.,1.-vgrat_high(i,j)-vgrat_low(i,j)) wlmx=wlmax*(bare+laihelp) wet_skin(i,j)=min(1.,src_dat(iglbsfc)%data(i,j,1)/wlmx) else wet_skin(i,j)=1.0 !for sea CMK bug 01-07-2003 end if ! ! -- calculation of water stress factor FWS with Wsmx is the ! field capacity. ! ! Field capacity is defined as the maximum amount of water ! that an column of soil can hold ! against gravity 24-48 hours after wetting of the soil ! Wcr is a critical value, Wpwp is the permanent wilting point ! and Ws is the total amount of water available in the ! root zone (ECHAM3 manual) ! fws(i,j)=1e-5 if (xland > 0) then wcr=0.5*ewsmx_sat ! used ewsmx_sat instead of ewsmx, is more consistent ! with values of swvl1 if ( swvl1_dat(iglbsfc)%data(i,j,1) > wpwp & .and. swvl1_dat(iglbsfc)%data(i,j,1) < wcr ) & fws(i,j)=(swvl1_dat(iglbsfc)%data(i,j,1)-wpwp)/(wcr-wpwp) if ( swvl1_dat(iglbsfc)%data(i,j,1) > wcr ) fws(i,j)=1. end if ! ! -- Computation of relative humidity (2 m) out of the air and dew ! temperature at 2 m which is used ! es=exp(19.59*(t2m_dat(iglbsfc)%data(i,j,1)-tstar)/t2m_dat(iglbsfc)%data(i,j,1)) e=exp(19.59*(1.-tstar/d2m_dat(iglbsfc)%data(i,j,1))) rh2m(i,j)=e/es ! !-- Calculation of H2 deposition over arable field ! Sanderson et al., J. Atmos. Chem. 2000 ! Yonemura et al., JGR 2000 ! units: m/sec ! ddepvel_h2(i,j) = min(10.e-4, max(1.e-7,(-41.39 * swvl1_dat(iglbsfc)%data(i,j,1) + 16.85) *1e-4)) ! ! calculate albedo (needed for photolysis rates) from ! different fractions ! ! vgr: coverage by low and high vegetation vgr= vgrat_low(i,j)+vgrat_high(i,j) ! ! -- if high vegetation is present and snow this may alter ! the effective snow albedo: let high vegetation prevail snow = snow_cover(i,j)- max(0.0,snow_cover(i,j)+vgrat_high(i,j)-1.0) bare = max(0.,1.0 - vgr - snow) ! soils with pH values larger than 7.3 desert = soilph(i,j,3) + soilph(i,j,4) ! when more desert is present than bare land... desert = desert - max(0.0,desert-bare) bare = bare -desert freesea = max(0.,1.-xland-ci_dat(iglbsfc)%data(i,j,1)) ! #ifndef without_photolysis !AJS: uv-vis albedo computed from land use; !AJS: in future, ecmwf field should be used: !AJS: name : UV visible albedo for direct radiation !AJS: abbrev : ALUVP ; unit : 0-1 !AJS: code : 15 ; table : 128 ags(i,j) = freesea*0.05 + ci_dat(iglbsfc)%data(i,j,1)*0.70 + & xland*(desert*0.10+bare*0.05+vgr*0.01+snow*0.70) ! AS: make sure that it does not exceed 0.7 in cases where ! sea ice and lsmask are not fully consistent ags(i,j) = min(0.7,ags(i,j)) #endif if (freesea>0.01) then ! ! bubble bursting effect,see Hummelshoj, equation 10 ! relationship by Wu (1988), note that Hummelshoj has not ! considered the cunningham factor which yields a different ! vb curve, with smaller values for small particles ! according to LG indeed 10 m windspeed (instead of 1 m windspeed in use in ECHAM5 alpha(i,j)=MIN(1.0,MAX(1.e-10,1.7e-6*wind10m(i,j)**3.75)) ! set maximum! qdrop=5.*(100.*alpha(i,j)) ! 100 is the flux of bubbles per cm^2/s (see Monohan(1988)) bubble(i,j)=((100.*ustar_sea(i,j))**2.)/(100.*wind10m(i,j)) + & eff*(2.*pi*rdrop**2.)*(2.*zdrop)*(qdrop/alpha(i,j)) !--- Correction of particle radius for particle growth close to the ! surface according to Fitzgerald, 1975, the relative humidity over ! the ocean is restricted to 98% (0.98) due to the salinity s=MIN(0.98,rh2m(i,j)) !fd23012004 beta(i,j)=EXP((0.00077*s)/(1.009-s)) THIS term was present in ECHAM code ! !; max. value reached for this parameter is 1.04 and is ignored here. phi=1.058-((0.0155*(s-0.97))/(1.02-s**1.4)) alpha1=1.2*EXP((0.066*s)/(phi-s)) vk1=10.2-23.7*s+14.5*s**2. vk2=-6.7+15.5*s-9.2*s**2. alpharat=1.-vk1*(1.-eps)-vk2*(1.-eps**2) alphae(i,j)=alpharat*alpha1 !--- Over land no correction for large humidity close to the surface: else ! land surface alpha(i,j)=0. bubble(i,j)=0. alphae(i,j)=1. endif end do end do if ( okdebug ) write(*,*) 'after dd_cal_inisurf' end subroutine dd_calc_inisurf !-------------------- ! ! This subroutine calculates the surface resistance as ! a function of a series of resistances ! ! purpose ! ------- ! determine surface resistance "rsurf" for dry deposition ! calculations ! ! external ! -------- ! ! reference ! --------- ! Ganzeveld and Lelieveld (1996) and references therein ! !------------------------------------------------------------------ ! ! -- resistances and auxiliary variables ! !------------------------------------------------------------------ subroutine dd_calc_rs use binas , only : vKarman, pi use chem_param, only : density_ref use chem_param, only : ndep, idep use chem_param, only : ddep_diffrb, ddep_rsoil, ddep_rwat, ddep_rws use chem_param, only : ddep_rsnow , ddep_rmes , ddep_rcut, ddep_diffcf use chem_param, only : diffcf_o3 use chem_param, only : iso2, iso4, inh3, ihno3, ino, ino2, io3, ico use chem_param, only : nrdep, lur use toolbox , only : coarsen_emission use ParTools, only : isRoot integer,parameter :: avg_field = 1 real, parameter :: dynvisc=1.789e-4*2. ! g cm-1 s-1 CHECKED FD OK, there is temp. dependence Perry p. 3-248. ! unit is g cm-1 s+1 FD ! checked with Seinfeld---> factor 2. came out (diameter ? radius) real, parameter :: cl=0.066*1e-4 ! mean free path [cm] (particle size also in cm) real, parameter :: bc= 1.38e-16 ! boltzman constant [g cm-2 s-1 K-1] (1.38e-23 J deg-1) =>binas real, parameter :: kappa=1. ! shapefactor real, parameter :: visc=0.15 ! KINEMATIC molecular viscocity [cm2 s-1] ! this is also function of temperature FD real :: xland1,xland2,xland3,xland4 real :: xwat1,xwat2,xsum real :: rstomx,rsveg,ra1 real :: vdland1,vdland2,vdland3,vdland4 real :: rssoil,rsws,vdwat1,vdwat2 real :: rssnow,rveg,rleaf,rcanopy real :: w10,zrsnhno3,ustarh real :: cvsh,cvwh,evgrath,tsoilph real :: xland,ts1,laihelp real,dimension(:,:), allocatable :: rwatso4 real,dimension(:,:), allocatable :: rsoilnh3,rsoilso2,rsoilso4 real,dimension(:,:), allocatable :: rsnowhno3,rsnowso2 ! integer :: jt,j,i,jdep, irdep character(len=8) :: adum real :: zr, vd_land, vd_sea, um, dc, cunning, relax, ust_land, st_land real :: sc, vb_land, vi_land, ust_sea, st_sea, re real :: vbsea, visea, vkdaccsea, vkd_sea, vkd_land real :: densaer real, allocatable :: curr_diffrb(:), curr_rsoil(:), curr_rwat(:), curr_rws(:) real, allocatable :: curr_rsnow (:), curr_rmes(:) , curr_rcut(:), curr_diffcf(:) ! Openmp parameters integer :: js, je ! --- begin --------------------------- ! ! *** deposition ! if ( ndep > 0 ) then allocate(rwatso4 (i1:i2, j1:j2)) allocate(rsoilnh3 (i1:i2, j1:j2)) allocate(rsoilso2 (i1:i2, j1:j2)) allocate(rsoilso4 (i1:i2, j1:j2)) allocate(rsnowhno3(i1:i2, j1:j2)) allocate(rsnowso2 (i1:i2, j1:j2)) allocate( curr_diffrb(ndep) ) allocate( curr_rsoil (ndep) ) allocate( curr_rwat (ndep) ) allocate( curr_rws (ndep) ) allocate( curr_rsnow (ndep) ) allocate( curr_rmes (ndep) ) allocate( curr_rcut (ndep) ) allocate( curr_diffcf(ndep) ) ! ! -- extension with the Wesely scheme, 1989 ! (Atmospheric Environ.,1989) ! in which the resistances of trace gases are estimated from ! the reactivity relative to ! ozone and the dissolvation relativeto SO2. The deposition process of ! these two trace is used as a reference for the estimation. ! ! calculate some tracer specific resistances ! curr_diffrb = ddep_diffrb curr_rsoil = ddep_rsoil curr_rwat = ddep_rwat curr_rws = ddep_rws curr_rsnow = ddep_rsnow curr_rmes = ddep_rmes curr_rcut = ddep_rcut curr_diffcf = ddep_diffcf ! do j=j1,j2 do i=i1,i2 ustarh=ustar_loc(i,j) ! -- calculation of the integrated resistance from the mass size ! distribution for rural continental aerosol with an unimodal distr. ! with the mean mass radius at about 0.4 um. ! (observations by Mehlmann (1986) ! as referenced in Warneck, 1988) and the friction velocity. ! Polynominal fit! ! Over land the surface resistance of bare soil and snow ! covered surfaces is calculated ! ! -- following distributions are used: ! 1. deposition over land using a rural continental size distribution ! i.e. mostly accumulation range aerosol ! 2. deposition over water using a marine size distribution ! i.e. bimodal distribution with a 30% coarse fraction ! fjd this is most relevant when explicit chemistry in seasalt ! is considered ! 3. deposition over water using rural continental size distribution ! the deposition of 'seasalt' sulfate may be calculated using ! the relation vd2=0.7*vd1+0.3*vdsssulfate ! ! 1st distribution: rsoilso4(i,j) = max( 1., & 100./(0.08-0.57*ustarh+1.66*ustarh**2-0.36*ustarh**3) ) ! 2nd distribution: ! rwatso4(i,j) = ! max(1.,100./(0.12-0.24*ustarh+5.25*ustarh**2 -1.84*ustarh**3)) ! 3rd distribution: rwatso4(i,j) = max( 1., & 100./(0.07-0.1*ustarh+2.1*ustarh**2 -0.20*ustarh**3) ) ! ! -- parameterization of snow/ice SO2 resistance as a function of the ! snow/ice temperature, based on observations by Dasch and Cadle ! ts1=max(200.,t2m_dat(iglbsfc)%data(i,j,1)) ! ! FD25032003 the measurements are not completely consistent; ! put a lower ! limit of vd=0.10 cm/s for SO2 to avoid excessive accumulation ! if snow is melting high uptake ! rsnowso2(i,j)=amin1(max(10.,10.**(-0.09*(ts1-273.)+2.4)),1.e+3) ! snow is melting changed (advise Wouter Greull) FD072003 if ( ts1 > (273.+3.) .and. sshf_dat(iglbsfc)%data(i,j,1) > 10. ) rsnowso2(i,j)=50. ! ! arbitrary attribution of soil resistance to NH3 deposition ! acid soils have lowest deposition ! tsoilph=sum(soilph(i,j,1:5)) rsoilnh3(i,j)=1e10 if (tsoilph > 0.) & rsoilnh3(i,j)= max(25.,soilph(i,j,1)*25. +soilph(i,j,2)*25.& +soilph(i,j,3)*200.+soilph(i,j,4)*200.+soilph(i,j,5)*100.) ! ! -- Computation of rsoil for SO2 from soil pH three different ! rh classes ! The rsoil values for each class are determined out of ! regression and the average value of pH of the class. ! Concerning rh, a wet, dry and very dry class ! is discerned with the threshold value of 60% and 40% ! rsoilso2(i,j)=1e10 ! ! -- for rh > 60 % ! if ( tsoilph > 0. ) rsoilso2(i,j)=max(25.,(soilph(i,j,1)*115.+& soilph(i,j,2)*65+soilph(i,j,3)*25.+ & soilph(i,j,4)*25.+soilph(i,j,5)*70.)/tsoilph+& amin1(max(0.,1000.*exp(269.-ts1)),1.e+5)) ! -- for rh less than 60% and larger than 40 % if ( rh2m(i,j) < 0.6 .and. rh2m(i,j) > 0.4) & !cmk rsoilso2(i,j)=max(50.,(rsoilso2(i,j)*3.41-85.)+ !cmk amin1(max(0.,1000.*exp(269.-ts1)),1.e+5)) rsoilso2(i,j)=max(50.,rsoilso2(i,j)*3.41-85.)+ & amin1(max(0.,1000.*exp(269.-ts1)),1.e+5) ! -- for rh less than 40 % if (rh2m(i,j).le.0.4) & !cmk rsoilso2(i,j)=max(50.,amin1(1000.,(rsoilso2(i,j)*3.41-& !cmk 85.+((0.4-rh2m(i,j))/0.4)*1.e+5)+& rsoilso2(i,j)=max(50.,amin1(1000.,(rsoilso2(i,j)*3.41-85.+ & ((0.4-rh2m(i,j))/0.4)*1.e+3)+& ! 1e3, see paper Laurens, JGR amin1(max(0.,1000.*exp(269.-ts1)),1.e+5))) ! -- Temperature correction term for HNO3 above snow surface and ice ! (Wesely, 1989), recomputation from K to 0C rsnowhno3(i,j)=amin1(max(10.,1000.*exp(269.-ts1)),1.e+5) ! snow is melting changed (advise Wouter Greull) FD072003 if ( ts1 > (273.+3.) .and. sshf_dat(iglbsfc)%data(i,j,1) > 10. ) rsnowhno3(i,j)=50. ! ! end do !nlon end do !nlat ! ! *** ! vd11(:,:,:) = 1e-10 !CMK rsurf(:,:,:)=1e+5 ! ! -- Loop for tracers, all timesteps ! do jdep=1,ndep jt = idep(jdep) ! ! -- Only if a value for DIFFCF (diffusivity) is defined, the program ! calculates a surface resistance ! if ( curr_diffcf(jdep) > 1e-5 ) then ! do j=j1,j2 do i=i1,i2 xland = lsmask_dat(iglbsfc)%data(i,j,1)/100.0 ! laihelp: the max values derived from ECMWF are more realistic laihelp=min(lai(i,j),lai1(i,j)) ! -- Assigning of values calculated in previous routine ! rsnow/rsoil of other components in data statement if (jt == iso2) then curr_rsnow(jdep)=rsnowso2(i,j) curr_rsoil(jdep)=rsoilso2(i,j) end if if (jt == inh3) curr_rsoil(jdep)=rsoilnh3(i,j) if (jt == iso4) then ! set snow resistance of so4 equal to soil resistance curr_rsnow(jdep)=rsoilso4(i,j) curr_rsoil(jdep)=rsoilso4(i,j) curr_rwat(jdep)=rwatso4(i,j) end if ! if (jt == ihno3) curr_rsnow(jdep)=rsnowhno3(i,j) ! -- Computation of stomatal resistance for component X, ! correction based on differences in molecular ! diffusion of H2O and X and the water ! stress is also considered, FWS rstomx=curr_diffcf(jdep)*rstom(i,j)/fws(i,j) ! -- Calculation of mesophyll resistance of NO ! and NO2 as function of ! the stomatal resistance of ozone if (jt == ino ) curr_rmes(jdep)= 5.*diffcf_o3*rstom(i,j)/fws(i,j) if (jt == ino2) curr_rmes(jdep)=0.5*diffcf_o3*rstom(i,j)/fws(i,j) rleaf=(1./((1./curr_rcut(jdep))+(1./(rstomx+curr_rmes(jdep))))) ! linear upscaling from leaf scale to canopy scale ! applying the LAI derived from Olson rcanopy=rleaf/max(1e-5,laihelp) rveg=(1./((1./(rahcan(i,j)+rb(i,j)*curr_diffrb(jdep)+curr_rsoil(jdep)))+ & (1./rcanopy))) ! -- Exception for CO: Use Sanderson et al. parameterization ! CO dry dep. scaled from H2 dry dep., using soil wetness. ! Note: This soil resistance is critical one, so one may ! neglect other terms (cuticle, boundary layer, ...) if (jt == ico) rveg = 1./(ddepvel_h2(i,j) * 0.6 ) !-- It can be assumed that HNO3 deposition only depends on ra ! so that surface resistance = 0, however definition of ! minimal surface resistance in order to avoid ! extreme Vd values. if (jt == ihno3) rveg=1. !-- Computation of surface resistance Rs (s m-1) above land ! Vd for surface with vegetation and already incorporated ! is the vegetation boudanry layer resistance rsveg=rveg+rb(i,j)*curr_diffrb(jdep) !-- Sulfate deposition calculation over vegetation according ! to the observations by Wesely (1985), ! see the Dry Deposition Inferential Model if (jt == iso4) then ! ! fd removed the factor stheta ! (standard deviation of the wind direction) if ( ssr_dat(iglbsfc)%data(i,j,1) > 250 ) then rsveg=1./(0.01*ustar_loc(i,j)) else rsveg=1./(0.002*ustar_loc(i,j)) end if ! ! -- for particle deposition the rb is already ! implicitly included rssoil=curr_rsoil(jdep) ! ! assumed that the wet skin fraction consists of vegetation ! rsws=rsveg rssnow=curr_rsnow(jdep) ! else !jt neq iso4 ! !-- Rs for surface with bare soil ! rssoil=curr_rsoil(jdep)+rb(i,j)*curr_diffrb(jdep) ! !-- Wet skin fraction, ! It is assumed that the wet skin fraction ! mainly consists of vegetation thus laminar ! resistance for vegetation may be applied ! rsws=curr_rws(jdep)+rb(i,j)*curr_diffrb(jdep) ! !-- Snow coverage ! rssnow=curr_rsnow(jdep)+rb(i,j)*curr_diffrb(jdep) ! end if ! !-- land surface deposition velocity ! cvsh=snow_cover(i,j) - & max(0.0,snow_cover(i,j)+vgrat_high(i,j)-1.0) ! fd same assumption as photolysis !cmk: need to correct the rsoil resistance here. ! fraction will become bare soil, which in the ! vegetation deposition calculation ! will result in a too low resistance, ! since the 'soil' part will be snow covered! cvwh=min(xland,wet_skin(i,j) ) evgrath=min(xland,vgrat_low(i,j)+vgrat_high(i,j)) ra1=raero(i,j) ! snow covered land fraction xland1=cvsh vdland1=1./(rssnow+ra1) ! wet skin covered land fraction (not covered by snow) xland2=max(0.,cvwh-xland1) vdland2=1./(rsws+ra1) ! vegetation covered land fraction ! (not covered by snow and wet skin) xland3=max(0.,evgrath-xland1-xland2) vdland3=1./(rsveg+ra1) ! bare soil covered land fraction ! (landfraction not covered by snow, wet skin, or vegetation) xland4=max(0.,xland-xland1-xland2-xland3) vdland4=1./(rssoil+ra1) xwat1=min(1-xland,ci_dat(iglbsfc)%data(i,j,1)) vdwat1=1./(curr_rsnow(jdep)+rb(i,j)*curr_diffrb(jdep)+ra1) xwat2 = max(0.,1.-xland-ci_dat(iglbsfc)%data(i,j,1)) vdwat2=1./(curr_rwat(jdep)+rb(i,j)*curr_diffrb(jdep)+ra1) if (jt == iso4) then vdwat1=1./(curr_rsnow(jdep)+ra1) vdwat2=1./(curr_rwat(jdep)+ra1) end if xsum=xland1+xland2+xland3+xland4+xwat1+xwat2 !fd if (xsum.ne.1.) then !fd print *,'sum',i,j,xsum,xland,xland1,xland2,'xl3',& !fd xland3,xland4,xwat1,xwat2,cvsh,cvwh,evgrath,ci_dat(iglbsfc)%data(i,j,1) !fd end if ! !-- Computation of deposition velocity ! vd11(i,j,jdep)=xland1*vdland1+xland2*vdland2+xland3*vdland3+& xland4*vdland4+xwat1*vdwat1+xwat2*vdwat2 end do !End of loop longitude end do !End of loop latitude adum=names(jt) i=len_trim(adum) if ( okdebug ) then call dist_arr_stat(dgrid(iglbsfc), 'vd'//adum(1:i), vd11(:,:,jdep), 0, status ) IF_NOTOK_RETURN(status=1) end if end if ! (curr_diffcf(jdep) > 1e-5) ! end do !End of loop tracer deallocate(rwatso4) deallocate(rsoilnh3) deallocate(rsoilso2) deallocate(rsoilso4) deallocate(rsnowhno3) deallocate(rsnowso2) deallocate( curr_diffrb ) deallocate( curr_rsoil ) deallocate( curr_rwat ) deallocate( curr_rws ) deallocate( curr_rsnow ) deallocate( curr_rmes ) deallocate( curr_rcut ) deallocate( curr_diffcf ) end if ! ndep > 0 ! ! *** aerosols ! if ( nrdep > 0 ) then ! look up table for different aerosol radii: do irdep = 1, nrdep !PLS !$OMP PARALLEL & !PLS !$OMP default (none) & !PLS !$OMP shared (irdep, lur, wind10m, lsmask_dat, alphae, t2m_dat, & !PLS !$OMP ustar_land, raero, ustar_sea, sr_sea, alpha, bubble, & !PLS !$OMP vd11a) & !PLS !$OMP private (i, j, js, je, zr, vd_land, vd_sea, um, xland, cunning, & !PLS !$OMP dc, densaer, relax, sc, ust_land, st_land, vb_land, & !PLS !$OMP vkd_land, ust_sea, st_sea, re, vbsea, visea, & !PLS !$OMP vkdaccsea, vkd_sea, vi_land) ! call my_omp_domain (1, nlat180, js, je) do j=j1,j2 do i=i1,i2 zr = 2.0e-4 * lur(irdep) ! diameter in cm ! vd_land=0. vd_sea=0. um=wind10m(i,j) xland = lsmask_dat(iglbsfc)%data(i,j,1)/100.0 ![]fraction !--- Cunningham factor: cunning=1.+(cl/(alphae(i,j)*zr))*(2.514+0.800*EXP(-0.55*zr/cl)) !--- Diffusivity: dc=(bc* t2m_dat(iglbsfc)%data(i,j,1)*cunning)/(3.*pi*dynvisc*(alphae(i,j)*zr)) ! [cm2/s] FD ! Relaxation: ! relax represents characteristic relaxation time scale [Seinfeld, p. 319] ! [ kg m-3 => g cm-3 ] densaer = density_ref ! reference density (e.q. 1800 kg/m3) relax=cunning*densaer*1.E-3*((alphae(i,j)*zr)**2. ) & /(18.*dynvisc*kappa) ! Sedimentation is calculated operator split in the subroutine sedimentation: ! ! sedspeed=(((((alphae(i,j)*zr(i,j)))**2.)* & ! densaer(jl,klev,jmod,jrow)*1.E-3*grav*cunning)/(18.*dynvisc)) ! note grav should be in cm ! [ kg m-3 => g cm-3 ] ! Calculation of schmidt sc =visc/dc ![cm2 s-1]/[cm2 s-1] dimensionless if (xland.gt.0.01) then ! note that in ECHAM there is a difference between vegetaton and snow/bare soil ! in TM5 there we assume this is incorporated in the 1x1 fields ust_land=ustar_land(i,j) ! ! calculation of stokes numbers st_land =max(0.1,(relax*(100.*ust_land)**2.)/visc) !--- Calculation of the dry deposition velocity ! See paper slinn and slinn, 1980, vd is related to d**2/3 ! over land, whereas over sea there is accounted for slip ! vb_ represents the contribution in vd of the brownian diffusion [cm s-1] ! and vi_ represents the impaction [cm s-1] ! vb_land =(1./vkarman)*((ust_land/um)**2)*100.*um*(sc**(-2./3.)) ![cm s-1] vi_land =(1./vkarman)*((ust_land/um)**2)*100.*um*(10.**(-3./st_land)) ![cm s-1] vkd_land =(vb_land+vi_land)*1e-2 ![m s-1] vd_land=1./(raero(i,j)+1./vkd_land) ! if (vkd_land.gt.1.) write(*,999) & ! 'after vb',i,j,'jtype',jtype,'jmod',jmod,'vb',vb_land,'vi',vi_land,& ! 'um10',um,'ust_land',ust_land,'dc',dc,'relax',relax,'schmidt',sc,'stokes',st_land endif ! xland.gt.0 if (xland.lt.0.99) then !--- Over sea: ! Brownian diffusion for rough elements, see Hummelshoj ! re is the reynolds stress: ust_sea=ustar_sea(i,j) st_sea =max(0.1,(relax*(100.*ust_sea)**2.)/visc) re =(100.*ust_sea*100.*sr_sea(i,j))/visc ! [cm/s]*[cm]/[cm2/s] vbsea =(1./3.)*100.*ust_sea*((sc**(-0.5))*re**(-0.5)) visea =100.*ust_sea*10.**(-3./st_sea) vkdaccsea =vbsea+visea vkd_sea =((1.-alpha(i,j))*vkdaccsea+alpha(i,j)*bubble(i,j))*1e-2 ! [m s-1] vd_sea=1./(raero(i,j)+1./vkd_sea) ! [m s-1] ! if (vkd_sea.gt.1.) write(*,999) 'sea',i,j,'jtype',jtype,'jmod',jmod,'vb',vbsea,& ! 'vi',visea,'vkd_sea',vkd_sea,'alpha*bubble',alpha(i,j)*bubble(i,j),'alpha',alpha(i,j) endif ! xland.lt.0.99 vd11a(i,j)= min(0.1,(1.-xland)*vd_sea + xland*vd_land) ! [m s-1] limit to 10 cm/s enddo ! j=1,nlat180 enddo ! i=1,nlon360 !$OMP END PARALLEL ! gather before coarsening call gather( dgrid(iglbsfc), vd11a, global_field, 0, status) IF_NOTOK_RETURN(status=1) if (isRoot) then call coarsen_emission('vd', imr, jmr, global_field, vd_temp_global, avg_field, status) IF_NOTOK_RETURN(status=1) end if ! scatter and store do region = 1, nregions call scatter( dgrid(region), vd_temp(region)%surf, vd_temp_global(region)%surf, 0, status) IF_NOTOK_RETURN(status=1) vda(region,irdep)%surf = vd_temp(region)%surf end do enddo ! irdep end if ! nrdep > 0 ! ! *** ! if ( okdebug ) write(*,*) 'dd_calc_rs: end' end subroutine dd_calc_rs end subroutine dd_surface_fields ! *** #ifndef without_photolysis subroutine dd_coarsen_vd( vd11, ags, I1, J1, status ) #else subroutine dd_coarsen_vd( vd11, I1, J1, status ) #endif use dims , only : iglbsfc, okdebug use chem_param, only : ndep, idep, vd_ncopy, vd_copy_itarget, vd_copy_isource use chem_param, only : ihno4, ih2o2, ino2 use chem_param, only : names use toolbox , only : coarsen_emission use toolbox , only : escape_tm #ifndef without_photolysis use photolysis_data, only : phot_dat #endif use tm5_distgrid, only : dgrid, Get_DistGrid, gather, scatter use ParTools, only : isRoot ! --- in/out ------------------------------- real, intent(inout) :: vd11(I1:,J1:,:) #ifndef without_photolysis real, intent(inout) :: ags(I1:,J1:) #endif integer, intent(in) :: i1, j1 integer, intent(out) :: status ! --- const ------------------------------- character(len=*), parameter :: rname = 'dd_coarsen_vd' integer,parameter :: avg_field = 1 ! --- local ------------------------------- integer :: i,j, region integer :: jt,jdep integer :: idno2,idh2o2,idhno4 !character(len=8) :: adum !character(len=2) :: bdum real, allocatable :: vdhelp(:,:) integer :: i01, i02, j01, j02, imr, jmr ! --- begin ------------------------------ CALL GET_DISTGRID( DGRID(iglbsfc), I_STRT=i01, I_STOP=i02, J_STRT=j01, J_STOP=j02, I_WRLD=imr, J_WRLD=jmr ) if (isRoot) then allocate(vdhelp(imr,jmr)) else allocate(vdhelp(1,1)) end if ! ! -- ! Scale the surface resistance of a number of trace gases and aerosols ! NOx deposition for all NOx components separately and scaling ! of NOx afterwards ! ! *** vd(HNO4) = ( vd(NO2) + vd(H2O2) ) / 2.0 ! search deposited tracer indices for NO2, H2O2, and HNO4 idno2 = -1 idh2o2 = -1 idhno4 = -1 ! try to reset ... do jdep = 1,ndep select case(idep(jdep)) case(ihno4) idhno4 = jdep case(ih2o2) idh2o2 = jdep case(ino2) idno2 = jdep case default end select end do ! HNO4 present ? then compute its deposition velocity: if ( idhno4 > 0 ) then ! check ... if ( any( (/idno2,idh2o2/) < 0 ) ) then write (gol,'("deposition velocity for HNO4 computed from NO2 and H2HO2, but at least one not found:")'); call goErr write (gol,'(" idno2, idh2o2 : ",2i4)') idno2, idh2o2; call goErr TRACEBACK; status=1; return end if ! vd(HNO4) = ( vd(NO2) + vd(H2O2) ) / 2.0 vd11(:,:,idhno4) = ( vd11(:,:,idno2) + vd11(:,:,idh2o2) ) / 2.0 end if ! *** do jdep=1,ndep jt = idep(jdep) call gather( dgrid(iglbsfc), vd11(:,:,jdep), vdhelp, 0, status) IF_NOTOK_RETURN(status=1) if (isRoot) then call coarsen_emission('vd', imr, jmr, vdhelp, vd_temp_global, avg_field, status) IF_NOTOK_RETURN(status=1) end if do region = 1,nregions call scatter( dgrid(region), vd_temp(region)%surf, vd_temp_global(region)%surf, 0, status) IF_NOTOK_RETURN(status=1) vd(region,jt)%surf = vd_temp(region)%surf !if ( okdebug ) then ! adum=names(jt) ! write(bdum,'(i2.2)') region ! i=len_trim(adum) ! !cmk call dumpfield(0,'vd.hdf',& ! !cmk vd(region,jt)%surf(:,:),'vd_'//adum(1:i)//bdum) !end if end do end do !jdep ! copy fields: do jdep = 1, vd_ncopy do region = 1, nregions vd(region,vd_copy_itarget(jdep))%surf = vd(region,vd_copy_isource(jdep))%surf end do end do #ifndef without_photolysis call gather( dgrid(iglbsfc), ags, vdhelp, 0, status) IF_NOTOK_RETURN(status=1) if (isRoot) then ! coarsen ags for photolysis...(temp in vd_temp) call coarsen_emission('ags',imr, jmr, vdhelp, vd_temp_global, avg_field, status) IF_NOTOK_RETURN(status=1) end if do region = 1,nregions call scatter( dgrid(region), vd_temp(region)%surf, vd_temp_global(region)%surf, 0, status) IF_NOTOK_RETURN(status=1) ! store coarsend albedo in photolysis data: phot_dat(region)%albedo = vd_temp(region)%surf ! update albedo statistics: !JEW: phot_dat(region)%ags_av = phot_dat(region)%ags_av + phot_dat(region)%albedo !JEW: phot_dat(region)%nalb_av = phot_dat(region)%nalb_av + 1 end do #endif ! Done deallocate(vdhelp) status = 0 end subroutine dd_coarsen_vd END MODULE DRY_DEPOSITION