!********************************************************************************** ! This computer software was prepared by Battelle Memorial Institute, hereinafter ! the Contractor, under Contract No. DE-AC05-76RL0 1830 with the Department of ! Energy (DOE). NEITHER THE GOVERNMENT NOR THE CONTRACTOR MAKES ANY WARRANTY, ! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE. ! ! MOSAIC module: see module_mosaic_driver.F for references and terms of use !********************************************************************************** module module_mosaic_drydep ! rce 2005-feb-18 - several fixes for indices of dlo/hi_sect, volumlo/hi_sect, ! which are now (isize,itype) ! rce 2004-dec-03 - many changes associated with the new aerosol "pointer" ! variables in module_data_mosaic_asect contains !----------------------------------------------------------------------- subroutine mosaic_drydep_driver( & id, curr_secs, ktau, dtstep, config_flags, & gmt, julday, & t_phy, rho_phy, p_phy, & ust, aer_res, & moist, chem, ddvel, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) use module_configure, only: grid_config_rec_type, num_moist, num_chem use module_state_description, only: param_first_scalar use module_data_mosaic_asect use module_data_mosaic_other use module_mosaic_driver, only: mapaer_tofrom_host use module_peg_util, only: peg_error_fatal implicit none ! subr arguments integer, intent(in) :: & id, ktau, julday integer, intent(in) :: & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte real(kind=8), intent(in) :: curr_secs real, intent(in) :: dtstep, gmt real, intent(in), & dimension( ims:ime, kms:kme, jms:jme ) :: & t_phy, rho_phy, p_phy real, intent(in), & dimension( ims:ime, jms:jme ) :: & ust real, intent(in), & dimension( its:ite, jts:jte ) :: & aer_res real, intent(in), & dimension( ims:ime, kms:kme, jms:jme, 1:num_moist ) :: & moist real, intent(inout), & dimension( ims:ime, kms:kme, jms:jme, 1:num_chem ) :: & chem real, intent(inout), & dimension( its:ite, jts:jte, 1:num_chem ) :: & ddvel type(grid_config_rec_type), intent(in) :: config_flags ! local variables integer idum, jdum integer it, jt, kt integer iphase, itype integer ktmaps, ktmape integer ll, l1, n integer levdbg_err, levdbg_info integer idiagaa_dum, ijcount_dum real dum real vdep_aer(maxd_asize,maxd_atype,maxd_aphase) character*100 msg !!rcetestdd diagnostics -------------------------------------------------- ! print 93010, ' ' ! print 93010, 'rcetestdd diagnostics from mosaic_drydep_driver' ! print 93010, 'id, chem_opt, ktau, julday ', & ! id, config_flags%chem_opt, ktau, julday ! print 93010, 'ims/e, j, k', ims, ime, jms, jme, kms, kme ! print 93010, 'its/e, j, k', its, ite, jts, jte, kts, kte ! ! do jdum = 0, 2 ! do idum = 0, 2 ! jt = jts + ((jte-jts)/2)*jdum ! it = its + ((ite-its)/2)*idum !if (idum .eq. 2) it = ite !if (jdum .eq. 2) jt = jte ! kt = kts ! print 93050, 'it, jt, t, p, ust, aer_res', it, jt, & ! t_phy(it,kt,jt), p_phy(it,kt,jt), ust(it,jt), aer_res(it,jt) ! end do ! end do ! !93010 format( a, 8(1x,i6) ) !93020 format( a, 8(1p,e14.6) ) !93050 format( a, 2(1x,i4), 8(1p,e14.6) ) !!rcetestdd diagnostics -------------------------------------------------- ! ktmaps,ktmape = first/last wrf kt for which depvel is calculated ! ktot = number of levels at which depvel is calculated ktmaps = kts ktmape = kts ktot = 1 ! set some variables to their wrf-chem "standard" values lunerr = -1 lunout = -1 levdbg_err = 0 levdbg_info = 15 iymdcur = 20030822 ihmscur = 0 ihmscur = nint( mod( real(gmt*3600.,8)+curr_secs, 86400.0_8 ) ) t = curr_secs ncorecnt = ktau - 1 ! reset some variables to "box test" values ! (*** aboxtest_get_extra_args is for "box testing" only ! and should be not be called from wrf-chem ***) ! call aboxtest_get_extra_args( 20, & ! iymdcur, ihmscur, & ! aboxtest_units_convert, aboxtest_rh_method, & ! aboxtest_map_method, aboxtest_gases_fixed, & ! lunerr, lunout, t, dum ) ! set "pegasus" grid size variables itot = ite jtot = jte nsubareas = 1 ijcount_dum = 0 do 2920 jt = jts, jte do 2910 it = its, ite ijcount_dum = ijcount_dum + 1 call mapaer_tofrom_host( 0, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, & it, jt, ktmaps,ktmape, & num_moist, num_chem, moist, chem, & t_phy, p_phy, rho_phy ) ! compute deposition velocities idiagaa_dum = 1 idiagaa_dum = 0 if ((jt.ne.jts) .and. (jt.ne.jte) .and. & (jt.ne.(jts+jte)/2)) idiagaa_dum = 0 if ((it.ne.its) .and. (it.ne.ite) .and. & (it.ne.(its+ite)/2)) idiagaa_dum = 0 call mosaic_drydep_1clm( idiagaa_dum, it, jt, & ust(it,jt), aer_res(it,jt), vdep_aer ) ! rce 23-nov-2004 - change to using the "..._aer" species pointers do iphase = 1, nphase_aer do itype = 1, ntype_aer do n = 1, nsize_aer(itype) do ll = -2, ncomp_plustracer_aer(itype) if (ll .eq. -2) then l1 = numptr_aer(n,itype,iphase) else if (ll .eq. -1) then l1 = -1 if (iphase .eq. ai_phase) l1 = waterptr_aer(n,itype) else if (ll .eq. 0) then l1 = -1 if (iphase .eq. ai_phase) l1 = hyswptr_aer(n,itype) else l1 = massptr_aer(ll,n,itype,iphase) end if if (l1 .ge. param_first_scalar) then ddvel(it,jt,l1) = vdep_aer(n,itype,iphase) end if end do end do end do end do 2910 continue 2920 continue ! print 93010, 'leaving mosaic_drydep_driver' return end subroutine mosaic_drydep_driver !----------------------------------------------------------------------- subroutine mosaic_drydep_1clm( idiagaa, it, jt, & ustar_in, depresist_a_in, vdep_aer ) use module_configure, only: grid_config_rec_type use module_data_mosaic_asect use module_data_mosaic_other use module_mosaic_driver, only: mapaer_tofrom_host use module_peg_util, only: peg_error_fatal implicit none ! subr arguments integer, intent(in) :: idiagaa, it, jt ! friction velocity (m/s) real, intent(in) :: ustar_in ! aerodynamic resistance (s/m) real, intent(in) :: depresist_a_in ! deposition velocities (m/s) real, intent(inout) :: vdep_aer(maxd_asize,maxd_atype,maxd_aphase) ! local variables real, parameter :: densdefault = 2.0 real, parameter :: smallmassaa = 1.0e-20 real, parameter :: smallmassbb = 1.0e-30 real, parameter :: piover6 = pi/6.0 real, parameter :: onethird = 1.0/3.0 integer iphase, itype, k, ll, l1, m, n real airdens, airkinvisc real depresist_a, depresist_unstabpblfact real depresist_d0, depresist_d3 real depvel_a0, depvel_a3 real drydens, drydp, drymass, dryvol real dum, dumalnsg, dumfact, dummass real freepath real rnum real temp real ustar real vsettl_0, vsettl_3 real wetdgnum, wetdens, wetdp, wetmass, wetvol ! if (idiagaa>0) print *, ' ' k = 1 m = 1 ! temperature (K) temp = rsub(ktemp,k,m) ! air density (g/cm^3) ! airdens = cairclm(1)*xmwair airdens = cairclm(1)*28.966 ! air kinematic viscosity (cm^2/s) airkinvisc = ( 1.8325e-4 * (416.16/(temp+120.0)) * & ((temp/296.16)**1.5) ) / airdens ! air molecular freepath (cm) freepath = 7.39758e-4 * airkinvisc / sqrt(temp) ! friction velocity (cm/s) ustar = ustar_in * 100.0 ! aerodynamic resistance (s/cm) depresist_a = depresist_a_in * 0.01 ! enhancement factor for unstable pbl depresist_unstabpblfact = 1.0 ! initialize vdep_aer vdep_aer(:,:,:) = 0.0 ! *** for now, just calc vdep_aer for iphase = ai_phase iphase = ai_phase ! calculate vdep_aer do 2900 itype = 1, ntype_aer do 2800 n = 1, nsize_aer(itype) dryvol = 0.0 drymass = 0.0 do ll = 1, ncomp_aer(itype) l1 = massptr_aer(ll,n,itype,iphase) dummass = rsub(l1,k,m)*mw_aer(ll,itype) drymass = drymass + dummass dryvol = dryvol + dummass/dens_aer(ll,itype) end do l1 = waterptr_aer(n,itype) dummass = rsub(l1,k,m)*mw_water_aer wetmass = drymass + dummass wetvol = dryvol + dummass/dens_water_aer l1 = numptr_aer(n,itype,iphase) rnum = rsub(l1,k,m) if (drymass .le. smallmassbb) then drydp = dcen_sect(n,itype) drydens = densdefault wetdp = drydp wetdens = drydens goto 1900 end if !jdf if (drymass .le. smallmassbb) then if (drymass .le. smallmassaa) then wetmass = drymass wetvol = dryvol end if drydens = drymass/dryvol wetdens = wetmass/wetvol if (rnum .ge. dryvol/volumlo_sect(n,itype)) then drydp = dlo_sect(n,itype) else if (rnum .le. dryvol/volumhi_sect(n,itype)) then drydp = dhi_sect(n,itype) else drydp = (dryvol/(piover6*rnum))**onethird end if !jdf dumfact = (wetvol/dryvol)**onethird !jdf dumfact = min( dumfact, 10.0 ) if(abs(wetvol).gt.(1000.*abs(dryvol))) then dumfact=10.0 else dumfact=abs(wetvol/dryvol)**onethird dumfact=max(1.0,min(dumfact,10.0)) endif !jdf wetdp = drydp*dumfact 1900 continue ! ! get surface resistance and settling velocity for mass (moment 3) ! and number (moment 0) ! dumalnsg = log( 1.0 ) wetdgnum = wetdp * exp( -1.5*dumalnsg*dumalnsg ) call aerosol_depvel_2( & wetdgnum, dumalnsg, wetdens, & temp, airdens, airkinvisc, freepath, & ustar, depresist_unstabpblfact, & depresist_d0, vsettl_0, & depresist_d3, vsettl_3 ) ! ! compute overall deposition velocity (binkowski/shankar eqn a33) ! dum = depresist_a + depresist_d3 + & depresist_a*depresist_d3*vsettl_3 depvel_a3 = vsettl_3 + (1./dum) dum = depresist_a + depresist_d0 + & depresist_a*depresist_d0*vsettl_0 depvel_a0 = vsettl_0 + (1./dum) ! cm/s --> m/s vdep_aer(n,itype,iphase) = depvel_a3 * 0.01 ! ! diagnostic output ! if (idiagaa>0) print 9310, it, jt, n, itype, iphase, & dcen_sect(n,itype), drydp, wetdp, & drydens, wetdens, vdep_aer(n,itype,iphase), & vsettl_3, depresist_d3, depresist_a 9310 format( 'aerdep', 2i4, 3i3, 1p, 3e10.2, & 2x, 0p, 2f5.2, 2x, 1p, 4e10.2 ) 2800 continue 2900 continue return end subroutine mosaic_drydep_1clm !------------------------------------------------------------------------ !------------------------------------------------------------------------ subroutine aerosol_depvel_2( & dgnum, alnsg, aerodens, & temp, airdens, airkinvisc, freepath, & ustar, depresist_unstabpblfact, & depresist_d0, vsettl_0, & depresist_d3, vsettl_3 ) ! ! computes the surface layer resistance term and the ! gravitational settling velocity term for the 3rd moment ! of a log-normal aerosol mode ! ! input parameters ! dgnum - geometric mean diameter for aerosol number (cm) ! alnsg - natural logarithm of the geometric standard deviation ! for aerosol number ! aerodens - aerosol density (dgnum and aerodens are for the ! actual wet distribution) ! temp - temperature (K) ! airdens - air density (g/cm^3) ! airkinvisc - air kinematic viscosity (cm^2/s) ! freepath - air molecular freepath (cm) ! ustar - friction velocity (cm/s) ! depresist_unstabpblfact = weseley et al. 1985 factor for increasing ! depvel in unstable pbl -- either ! 1. + (-0.3*zi/L)**0.667 OR ! 1. + 0.24*((wstar/ustar)**2) ! output parameters ! depresist_d3 - surface layer resistance for 3rd (mass) moment (s/cm) ! vsettl_3 - gravitational settling velocity for 3rd moment (cm/s) ! depresist_d0 - surface layer resistance for 0th (number) moment (s/cm) ! vsettl_0 - gravitational settling velocity for 0th moment (cm/s) ! implicit none real dgnum, alnsg, aerodens, & temp, airdens, airkinvisc, freepath, & ustar, depresist_unstabpblfact, & depresist_d0, vsettl_0, & depresist_d3, vsettl_3 real aerodiffus_0, schmidt_0, stokes_0, facdepresist_d0 real aerodiffus_3, schmidt_3, stokes_3, facdepresist_d3 common / aerosol_depvel_cmn01 / & aerodiffus_0, schmidt_0, stokes_0, facdepresist_d0, & aerodiffus_3, schmidt_3, stokes_3, facdepresist_d3 real xknudsen, xknudsenfact, alnsg2, duma, dumb, & vsettl_dgnum, aerodiffus_dgnum real pi parameter (pi = 3.1415926536) ! gravity = gravitational acceleration in cm/s^2 real gravity parameter (gravity = 980.616) ! boltzmann constant in erg/deg-K real boltzmann parameter (boltzmann = 1.3807e-16) xknudsen = 2.*freepath/dgnum xknudsenfact = xknudsen*1.246 alnsg2 = alnsg*alnsg vsettl_dgnum = (gravity*aerodens*dgnum*dgnum)/ & (18.*airkinvisc*airdens) vsettl_0 = vsettl_dgnum * & ( exp(2.*alnsg2) + xknudsenfact*exp(0.5*alnsg2) ) vsettl_3 = vsettl_dgnum * & ( exp(8.*alnsg2) + xknudsenfact*exp(3.5*alnsg2) ) aerodiffus_dgnum = (boltzmann*temp)/ & (3.*pi*airkinvisc*airdens*dgnum) aerodiffus_0 = aerodiffus_dgnum * & ( exp(+0.5*alnsg2) + xknudsenfact*exp(+2.*alnsg2) ) aerodiffus_3 = aerodiffus_dgnum * & ( exp(-2.5*alnsg2) + xknudsenfact*exp(-4.*alnsg2) ) schmidt_0 = airkinvisc/aerodiffus_0 schmidt_3 = airkinvisc/aerodiffus_3 stokes_0 = ustar*ustar*vsettl_0/(gravity*airkinvisc) stokes_3 = ustar*ustar*vsettl_3/(gravity*airkinvisc) duma = (schmidt_0**(-0.66666666)) + & (10.**(-3./max(0.03,stokes_0))) ! dumb = duma*ustar*(1. + 0.24*wstaroverustar*wstaroverustar) dumb = duma*ustar*depresist_unstabpblfact depresist_d0 = 1./max( dumb, 1.e-20 ) facdepresist_d0 = duma duma = (schmidt_3**(-0.66666666)) + & (10.**(-3./max(0.03,stokes_3))) ! dumb = duma*ustar*(1. + 0.24*wstaroverustar*wstaroverustar) dumb = duma*ustar*depresist_unstabpblfact depresist_d3 = 1./max( dumb, 1.e-20 ) facdepresist_d3 = duma return end subroutine aerosol_depvel_2 !----------------------------------------------------------------------- end module module_mosaic_drydep