MODULE module_gocart_drydep CONTAINS subroutine gocart_drydep_driver(dtstep, & config_flags,numgas, & t_phy,moist,p8w,t8w,rmol,aer_res, & p_phy,chem,rho_phy,dz8w,ddvel,xland,hfx, & ivgtyp,tsk,vegfra,pbl,ust,znt,xlat,xlong, & dustdrydep_1,dustdrydep_2,dustdrydep_3, & dustdrydep_4,dustdrydep_5, & depvelocity, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) USE module_model_constants USE module_configure USE module_state_description IMPLICIT NONE INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte,numgas INTEGER,DIMENSION( ims:ime , jms:jme ) , & INTENT(IN ) :: & ivgtyp REAL, INTENT(IN ) :: dtstep REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ), & INTENT(IN ) :: moist REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), & INTENT(INOUT ) :: chem REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , & INTENT(IN ) :: & t_phy, & p_phy, & dz8w, & t8w,p8w,rho_phy REAL, DIMENSION( its:ite, jts:jte, num_chem ), & INTENT(INOUT) :: ddvel REAL, DIMENSION( ims:ime , jms:jme ) , & INTENT(IN) :: & tsk, & vegfra, & pbl, & ust, & xlat, & xlong, & rmol,xland,znt,hfx REAL, DIMENSION( its:ite, jts:jte ), & INTENT(IN) :: aer_res REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: & dustdrydep_1, dustdrydep_2, dustdrydep_3, & dustdrydep_4, dustdrydep_5, depvelocity !! .. Local Scalars .. INTEGER :: iland, iprt, iseason, jce, jcs, & n, nr, ipr, jpr, nvr, & idrydep_onoff,imx,jmx,lmx integer :: ii,jj,kk,i,j,k,nv integer,dimension (1,1) :: ilwi,ireg REAL :: clwchem, dvfog, dvpart, & rad, rhchem, ta, vegfrac, z1,zntt ! real*8, DIMENSION (1,1,1,3) :: erodin ! real*8, DIMENSION (5) :: tc,bems ! real*8, dimension (1,1) :: z0,w10m,gwet,airden,airmas,delz_sfc,hflux,ts,pblz,ustar,ps real*8, dimension (1,1) :: z0,airden,delz_sfc,hflux,ts,pblz,ustar,ps REAL*8 :: dvel(1,1), drydf(1,1) LOGICAL :: highnh3, rainflag, vegflag, wetflag TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags do nv=numgas+1,num_chem do j=jts,jte do i=its,ite ddvel(i,j,nv) = 0. enddo enddo enddo imx = 1 jmx = 1 lmx = 1 do j = max(jds+1,jts),min(jde-1,jte) do i = max(ids+1,its),min(ide-1,ite) ipr = 0 dvel(1,1) = 0._8 if( xland(i,j) > 1.5 ) then ilwi(1,1) = 1 else ilwi(1,1) = 0 endif ! for aerosols, ii=1 or ii=2 ii = 1 ! if(ivgtyp(i,j).eq.19.or.ivgtyp(i,j).eq.23)ii=1 ireg(1,1) = 1 airden(1,1) = real( rho_phy(i,kts,j),kind=8 ) delz_sfc(1,1) = real( dz8w(i,kts,j),kind=8 ) ustar(1,1) = max( 1.e-1_8,real(ust(i,j),kind=8) ) hflux(1,1) = real( hfx(i,j),kind=8 ) pblz(1,1) = real( pbl(i,j),kind=8 ) ps(1,1) = real(p8w(i,kts,j),kind=8)*.01_8 z0(1,1) = real( znt(i,j),kind=8 ) ts(1,1) = real( tsk(i,j),kind=8 ) ! if(i.eq.23.and.j.eq.74)ipr=1 call depvel_gocart(config_flags,ipr,ii,imx,jmx,lmx,& airden, delz_sfc, pblz, ts, ustar, hflux, ilwi, & ps, z0, dvel, drydf,g,rmol(i,j),aer_res(i,j)) do nv = p_p25,num_chem ddvel(i,j,nv) = real( dvel(1,1),kind=4 ) enddo ddvel(i,j,p_sulf) = real( dvel(1,1),kind=4 ) ddvel(i,j,p_msa) = real( dvel(1,1),kind=4 ) depvelocity(i,j) = ddvel(i,j,p_dust_5) ! drydep [ug/m2/s] = dvel [m/s] * chem [ug/kg] * airden [kg/m3] dustdrydep_1(i,j)=-dvel(1,1)*chem(i,1,j,p_dust_1)*airden(1,1) dustdrydep_2(i,j)=-dvel(1,1)*chem(i,1,j,p_dust_2)*airden(1,1) dustdrydep_3(i,j)=-dvel(1,1)*chem(i,1,j,p_dust_3)*airden(1,1) dustdrydep_4(i,j)=-dvel(1,1)*chem(i,1,j,p_dust_4)*airden(1,1) dustdrydep_5(i,j)=-dvel(1,1)*chem(i,1,j,p_dust_5)*airden(1,1) enddo enddo end subroutine gocart_drydep_driver SUBROUTINE depvel_gocart( config_flags,ipr,ii,imx,jmx,lmx, & airden, delz_sfc, pblz, ts, ustar, hflux, ilwi, & ps, z0, dvel, drydf,g0,rmol,aer_res) ! **************************************************************************** ! * * ! * Calculate dry deposition velocity. * ! * * ! * Input variables: * ! * AEROSOL(k) - Logical, T = aerosol species, F = gas species * ! * IREG(i,j) - # of landtypes in grid square * ! * ILAND(i,j,ldt) - Land type ID for element ldt =1,IREG(i,j) * ! * IUSE(i,j,ldt) - Fraction of gridbox area occupied by land type * ! * element ldt * ! * USTAR(i,j) - Friction velocity (m s-1) * ! * DELZ_SFC(i,j) - Thickness of layer above surface * ! * PBLZ(i,j) - Mixing depth (m) * ! * Z0(i,j) - Roughness height (m) * ! * * ! * Determined in this subroutine (local): * ! * OBK - Monin-Obukhov length (m): set to 1.E5 m under * ! * neutral conditions * ! * Rs(ldt) - Bulk surface resistance(s m-1) for species k to * ! * surface ldt * ! * Ra - Aerodynamic resistance. * ! * Rb - Sublayer resistance. * ! * Rs - Surface resistance. * ! * Rttl - Total deposition resistance (s m-1) for species k * ! * Rttl(k) = Ra + Rb + Rs. * ! * * ! * Returned: * ! * DVEL(i,j,k) - Deposition velocity (m s-1) of species k * ! * DRYDf(i,j,k) - Deposition frequency (s-1) of species k, * ! * = DVEL / DELZ_SFC * ! * * ! **************************************************************************** USE module_configure IMPLICIT NONE REAL*8, INTENT(IN) :: airden(imx,jmx), delz_sfc(imx,jmx) REAL*8, INTENT(IN) :: hflux(imx,jmx), ts(imx,jmx) REAL*8, INTENT(IN) :: ustar(imx,jmx), pblz(imx,jmx) REAL*8, INTENT(IN) :: ps(imx,jmx) INTEGER, INTENT(IN) :: ilwi(imx,jmx) INTEGER, INTENT(IN) :: imx,jmx,lmx REAL*8, INTENT(IN) :: z0(imx,jmx) REAL, INTENT(IN) :: g0,rmol,aer_res REAL*8, INTENT(OUT) :: dvel(imx,jmx), drydf(imx,jmx) TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags REAL*8 :: obk, vds, czh, rttl, frac, logmfrac, psi_h, cz, eps REAL*8 :: vd, ra, rb, rs INTEGER :: ipr,i, j, k, ldt, iolson, ii CHARACTER(LEN=50) :: msg REAL*8 :: prss, tempk, tempc, xnu, ckustr, reyno, aird, diam, xm, z REAL*8 :: frpath, speed, dg, dw, rt REAL*8 :: rad0, rix, gfact, gfaci, rdc, rixx, rluxx, rgsx, rclx REAL*8 :: dtmp1, dtmp2, dtmp3, dtmp4 REAL*8 :: biofit,vk ! executable statements j_loop: DO j = 1,jmx i_loop: DO i = 1,imx vk = .4_8 vd = 0._8 ra = 0._8 rb = 0._8 ! only required for gases (SO2) rs = 0.0_8 ! **************************************************************************** ! * Compute the the Monin-Obhukov length. * ! * The direct computation of the Monin-Obhukov length is: * ! * * ! * - Air density * Cp * T(surface air) * Ustar^3 * ! * OBK = ---------------------------------------------- * ! * vK * g * Sensible Heat flux * ! * * ! * Cp = 1000 J/kg/K = specific heat at constant pressure * ! * vK = 0.4 = von Karman's constant * ! **************************************************************************** IF (abs(hflux(i,j)) <= 1.e-5_8) THEN obk = 1.0E5_8 ELSE ! MINVAL(hflux), MINVAL(airden), MINVAL(ustar) =?? obk = -airden(i,j) * 1000.0_8 * ts(i,j) * (ustar(i,j))**3 & / (vk * real(g0,kind=8) * hflux(i,j)) ! -- debug: IF ( obk == 0.0_8 ) WRITE(*,211) obk, i, j 211 FORMAT(1X,'OBK=', E11.2, 1X,' i,j = ', 2I4) END IF ! write(0,*)1./obk,rmol if(rmol.ne.0.)then obk=1._8/real(rmol,kind=8) else obk=1.e5_8 endif ! cz = delz_sfc(i,j) / 2.0_8 ! center of the grid box above surface cz = 2._8 ! **************************************************************************** ! * (1) Aerosodynamic resistance Ra and sublayer resistance Rb. * ! * * ! * The Reynolds number REYNO diagnoses whether a surface is * ! * aerodynamically rough (REYNO > 10) or smooth. Surface is * ! * rough in all cases except over water with low wind speeds. * ! * * ! * For gas species over land and ice (REYNO >= 10) and for aerosol * ! * species for all surfaces: * ! * * ! * Ra = 1./VT (VT from GEOS Kzz at L=1, m/s). * ! * * ! * The following equations are from Walcek et al, 1986: * ! * * ! * For gas species when REYNO < 10 (smooth), Ra and Rb are combined * ! * as Ra: * ! * * ! * Ra = { ln(ku* z1/Dg) - Sh } / ku* eq.(13) * ! * * ! * where z1 is the altitude at the center of the lowest model layer * ! * (CZ); * ! * Sh is a stability correction function; * ! * k is the von Karman constant (0.4, vK); * ! * u* is the friction velocity (USTAR). * ! * * ! * Sh is computed as a function of z1 and L eq ( 4) and (5)): * ! * * ! * 0 < z1/L <= 1: Sh = -5 * z1/L * ! * z1/L < 0: Sh = exp{ 0.598 + 0.39*ln(E) - 0.09(ln(E))^2 } * ! * where E = min(1,-z1/L) (Balkanski, thesis). * ! * * ! * For gas species when REYNO >= 10, * ! * * ! * Rb = 2/ku* (Dair/Dg)**(2/3) eq.(12) * ! * where Dg is the gas diffusivity, and * ! * Dair is the air diffusivity. * ! * * ! * For aerosol species, Rb is combined with surface resistance as Rs. * ! * * ! **************************************************************************** frac = cz / obk IF (frac > 1.0_8) frac = 1.0_8 IF (frac > 0.0_8 .AND. frac <= 1.0_8) THEN psi_h = -5.0_8*frac ELSE IF (frac < 0.0_8) THEN eps = MIN(1.0_8, -frac) logmfrac = LOG(eps) psi_h = EXP( 0.598_8 + 0.39_8 * logmfrac - 0.09_8 * (logmfrac)**2 ) END IF !-------------------------------------------------------------- ! Aerosol species, Rs here is the combination of Rb and Rs. ra = (LOG(cz/z0(i,j)) - psi_h) / (vk*ustar(i,j)) vds = 0.002_8*ustar(i,j) IF (obk < 0.0_8) vds = vds * (1.0_8+(-300.0_8/obk)**0.6667_8) czh = pblz(i,j)/obk IF (czh < -30.0_8) vds = 0.0009_8*ustar(i,j)*(-czh)**0.6667_8 if(ipr.eq.1) write(0,*)ra,aer_res,vds IF( config_flags%chem_opt /= CHEM_VASH .and. & config_flags%chem_opt /= DUST )THEN ra = real(aer_res,kind=8) ENDIF ! --Set Vds to be less than VDSMAX (entry in input file divided -- ! by 1.E4). VDSMAX is taken from Table 2 of Walcek et al. [1986]. ! Invert to get corresponding R ! if(ii.eq.1) then ! rs=1.0_8/MIN(vds,2.0e-2_8) ! else rs=1.0_8/MIN(vds,2.0e-3_8) ! endif ! ------ Set max and min values for bulk surface resistances ------ rs= MAX(1.0_8, MIN(rs, 9.9990e+3_8)) ! **************************************************************************** ! * * ! * Compute dry deposition velocity. * ! * * ! * IUSE is the fraction of the grid square occupied by surface ldt in * ! * units of per mil (IUSE=500 -> 50% of the grid square). Add the * ! * contribution of surface type ldt to the deposition velocity; this is * ! * a loop over all surface types in the gridbox. * ! * * ! * Total resistance = Ra + Rb + Rs. ! * * ! **************************************************************************** rttl = ra + rb + rs vd = vd + 1._8/rttl ! ------ Load array DVEL ------ if(ipr.eq.1) write(0,*)rs,ra,rb,vd dvel(i,j) = vd !* 1.2 ! -- Set a minimum value for DVEL ! MIN(VdSO2) = 2.0e-3 m/s over ice ! = 3.0e-3 m/s over land ! MIN(vd_aerosol) = 1.0e-4 m/s IF (dvel(i,j) < 1.0E-4_8) dvel(i,j) = 1.0E-4_8 drydf(i,j) = dvel(i,j) / delz_sfc(i,j) END DO i_loop END DO j_loop END SUBROUTINE depvel_gocart END MODULE module_gocart_drydep