13 subroutine samfdeepcnv_aerosols(im, ix, km, itc, ntc, ntr, delt,
14 & xlamde, xlamdd, cnvflg, jmin, kb, kmax, kd94, ktcon, fscav,
15 & edto, xlamd, xmb, c0t, eta, etad, zi, xlamue, xlamud, delp,
18 use machine ,
only : kind_phys
19 use physcons,
only : g => con_g, qamin
24 integer,
intent(in) :: im,
25 & ix, km, itc, ntc, ntr
26 real(kind=kind_phys),
intent(in) :: delt,
28 logical,
dimension(im),
intent(in) :: cnvflg
29 integer,
dimension(im),
intent(in) :: jmin,
30 & kb, kmax, kd94, ktcon
31 real(kind=kind_phys),
dimension(im),
intent(in) :: edto,
33 real(kind=kind_phys),
dimension(ntc),
intent(in) :: fscav
34 real(kind=kind_phys),
dimension(im,km),
intent(in) :: c0t,
35 & eta, etad, zi, xlamue, xlamud
36 real(kind=kind_phys),
dimension(ix,km),
intent(in) :: delp
37 real(kind=kind_phys),
dimension(ix,km,ntr+2),
intent(in) :: qtr
39 real(kind=kind_phys),
dimension(im,km,ntc),
intent(out) :: qaero
43 integer :: i, indx, it, k, kk, km1, kp1, n
44 real(kind=kind_phys) :: adw, aup, dtime_max, dv1q, dv2q, dv3q,
45 & dtovdz, dz, factor, ptem, ptem1, qamax, tem, tem1
46 real(kind=kind_phys),
dimension(ix,km) :: xmbp
47c -- chemical transport variables
48 real(kind=kind_phys),
dimension(im,km,ntc) :: ctro2, ecko2, ecdo2,
50c -- additional variables for tracers for wet deposition,
51 real(kind=kind_phys),
dimension(im,km,ntc) :: chem_c, chem_pw,
53c --
if reevaporation is enabled, uncomment lines below
54c real(kind=kind_phys),
dimension(im,ntc) :: pwav
55c real(kind=kind_phys),
dimension(im,km) :: pwdper
56c real(kind=kind_phys),
dimension(im,km,ntr) :: chem_pwd
57c -- additional variables for fct
58 real(kind=kind_phys),
dimension(im,km) :: flx_lo, totlout, clipout
60 real(kind=kind_phys),
parameter :: one = 1.0_kind_phys
61 real(kind=kind_phys),
parameter :: half = 0.5_kind_phys
62 real(kind=kind_phys),
parameter :: quarter = 0.25_kind_phys
63 real(kind=kind_phys),
parameter :: zero = 0.0_kind_phys
64 real(kind=kind_phys),
parameter :: epsil = 1.e-22_kind_phys
68c -- check
if aerosols are
present
69 if ( ntc <= 0 .or. itc <= 0 .or. ntr <= 0 )
return
70 if ( ntr < itc + ntc - 3 )
return
72c -- initialize work variables
90 if (k <= kmax(i)) qaero(i,k,n) = max(qamin, qtr(i,k,it))
97 if (cnvflg(i)) xmbp(i,k) = g * xmb(i) / delp(i,k)
106 if (kp1 <= kmax(i)) ctro2(i,k,n) =
107 & half * (qaero(i,k,n) + qaero(i,kp1,n))
112 ctro2(i,kmax(i),n) = qaero(i,kmax(i),n)
119 if (cnvflg(i) .and. (k <= kb(i)))
120 & ecko2(i,k,n) = ctro2(i,k,n)
127 if (cnvflg(i)) ecdo2(i,jmin(i),n) = ctro2(i,jmin(i),n)
131c
do chemical tracers, first need to know how much reevaporates
133c aerosol re-evaporation is set to zero for now
134c uncomment and edit the following code to enable re-evaporation
140c pwdper(i,k)= -edto(i)*pwdo(i,k)/pwavo(i)
144c calculate include mixing ratio(ecko2), how much goes into
145c rainwater to be rained out(chem_pw), and total scavenged,
146c
if not reevaporated(pwav)
153 if ((k > kb(i)) .and. (k < ktcon(i)))
then
154 dz = zi(i,k) - zi(i,kk)
155 tem = half * (xlamue(i,k)+xlamue(i,kk)) * dz
156 tem1 = quarter * (xlamud(i,k)+xlamud(i,kk)) * dz
157 factor = one + tem - tem1
159c
if conserved(not scavenging)
then
160 ecko2(i,k,n) = ((one-tem1)*ecko2(i,kk,n)
161 & + half*tem*(ctro2(i,k,n)+ctro2(i,kk,n)))/factor
163c how much will be scavenged
165c this choice was used in gf, and is also described in a
166c successful implementation into cesm in grl(yu et al. 2019),
167c it uses dimesnsionless scavenging coefficients(fscav),
168c but includes henry coeffs with gas phase chemistry
170c fraction fscav is going into liquid
171 chem_c(i,k,n)=fscav(n)*ecko2(i,k,n)
173c of that part is going into rain out(chem_pw)
174 tem=chem_c(i,k,n)/(one+c0t(i,k)*dz)
175 chem_pw(i,k,n)=c0t(i,k)*dz*tem*eta(i,kk)
176 ecko2(i,k,n)=tem+ecko2(i,k,n)-chem_c(i,k,n)
178c pwav needed fo reevaporation in downdraft
179c
if including reevaporation, please uncomment code below
180c pwav(i,n)=pwav(i,n)+chem_pw(i,k,n)
187 if (k >= ktcon(i)) ecko2(i,k,n)=ctro2(i,k,n)
192c reevaporation of some, pw and pwd terms needed later for dellae2
198 if (cnvflg(i) .and. (k < jmin(i)))
then
199 dz = zi(i,kp1) - zi(i,k)
200 if (k >= kd94(i))
then
202 tem1 = half * xlamdd * dz
205 tem1 = half * (xlamd(i)+xlamdd) * dz
207 factor = one + tem - tem1
208 ecdo2(i,k,n) = ((one-tem1)*ecdo2(i,kp1,n)
209 & +half*tem*(ctro2(i,k,n)+ctro2(i,kp1,n)))/factor
210c
if including reevaporation, please uncomment code below
211c ecdo2(i,k,n)=ecdo2(i,k,n)+pwdper(i,kp1)*pwav(i,n)
212c chem_pwd(i,k,n)=max(zero,pwdper(i,kp1)*pwav(i,n))
221c subsidence term treated in fct routine
222 dellae2(i,1,n) = edto(i)*etad(i,1)*ecdo2(i,1,n)*xmbp(i,1)
232c for the subsidence term already is considered
233 dellae2(i,k,n) = eta(i,kk) * ecko2(i,kk,n) * xmbp(i,k)
238c --- for updraft & downdraft vertical transport
240c initialize maximum allowed timestep for upstream difference approach
246 if (kk < ktcon(i)) dtime_max = min(dtime_max,half*delp(i,kk))
250c now for every chemistry tracer
255 if (cnvflg(i) .and. (k < ktcon(i)))
then
256 dz = zi(i,k) - zi(i,kk)
258 if (k <= kb(i)) aup = zero
260 if (k > jmin(i)) adw = zero
262 dv1q = half * (ecko2(i,k,n) + ecko2(i,kk,n))
263 dv2q = half * (ctro2(i,k,n) + ctro2(i,kk,n))
264 dv3q = half * (ecdo2(i,k,n) + ecdo2(i,kk,n))
266 tem = half * (xlamue(i,k) + xlamue(i,kk))
267 tem1 = half * (xlamud(i,k) + xlamud(i,kk))
269 if (k <= kd94(i))
then
271 ptem1 = xlamd(i) + xlamdd
276 dellae2(i,k,n) = dellae2(i,k,n) +
277c detrainment from updraft
278 & ( aup*tem1*eta(i,kk)*dv1q
279c entrainement into up and downdraft
280 & - (aup*tem*eta(i,kk)+adw*edto(i)*ptem*etad(i,k))*dv2q
281c detrainment from downdraft
282 & + (adw*edto(i)*ptem1*etad(i,k)*dv3q) ) * dz * xmbp(i,k)
284 wet_dep(i,k,n)=chem_pw(i,k,n)*g/delp(i,k)
286c sinks from
where updraft and downdraft start
287 if (k == jmin(i)+1)
then
288 dellae2(i,k,n) = dellae2(i,k,n)
289 & -edto(i)*etad(i,kk)*ctro2(i,kk,n)*xmbp(i,k)
292 dellae2(i,k,n) = dellae2(i,k,n)
293 & -eta(i,k)*ctro2(i,k,n)*xmbp(i,k)
303 dellae2(i,k,n) = dellae2(i,k,n)
304 & -eta(i,k)*ctro2(i,k,n)*xmbp(i,k)
317c compute low-order mass flux, upstream
321 if (cnvflg(i) .and. (kk < ktcon(i)))
then
323 if (kk >= kb(i) ) tem = eta(i,kk)
324 if (kk <= jmin(i)) tem = tem - edto(i)*etad(i,kk)
325c low-order flux,upstream
327 flx_lo(i,k) = -xmb(i) * tem * qaero(i,k,n)
328 elseif (tem < zero)
then
329 flx_lo(i,k) = -xmb(i) * tem * qaero(i,kk,n)
335c --- make sure low-ord fluxes don
't violate positive-definiteness
339.and.
if (cnvflg(i) (k <= ktcon(i))) then
340c time step / grid spacing
341 dtovdz = g * dtime_max / abs(delp(i,k))
343 totlout(i,k)=max(zero,flx_lo(i,kp1))-min(zero,flx_lo(i,k))
344 clipout(i,k)=min(one ,qaero(i,k,n)/max(epsil,totlout(i,k))
345 & / (1.0001_kind_phys*dtovdz))
350c recompute upstream mass fluxes
354.and.
if (cnvflg(i) (kk < ktcon(i))) then
356 if (kk >= kb(i) ) tem = eta(i,kk)
357 if (kk <= jmin(i)) tem = tem - edto(i)*etad(i,kk)
359 flx_lo(i,k) = flx_lo(i,k) * clipout(i,k)
360 elseif (tem < zero) then
361 flx_lo(i,k) = flx_lo(i,k) * clipout(i,kk)
367c --- a positive-definite low-order (diffusive) solution for the subsidnce fluxes
371.and.
if (cnvflg(i) (k <= ktcon(i))) then
372 dtovdz = g * dtime_max / abs(delp(i,k)) ! time step /grid spacing
373 dellae2(i,k,n) = dellae2(i,k,n)
374 & -(flx_lo(i,kp1)-flx_lo(i,k))*dtovdz/dtime_max
381c convert wet deposition to total mass deposited over dt and dp
386.and.
if (cnvflg(i) (k < ktcon(i)))
387 & wet_dep(i,k,n) = wet_dep(i,k,n)*xmb(i)*delt*delp(i,k)
392c compute final aerosol concentrations
397.and.
if (cnvflg(i) (k <= min(kmax(i),ktcon(i)))) then
398 qaero(i,k,n) = qaero(i,k,n) + dellae2(i,k,n) * delt
399 if (qaero(i,k,n) < zero) then
400c add negative mass to wet deposition
401 wet_dep(i,k,n) = wet_dep(i,k,n)-qaero(i,k,n)*delp(i,k)
410 end subroutine samfdeepcnv_aerosols
412 subroutine samfshalcnv_aerosols(im, ix, km, itc, ntc, ntr, delt,
413 & cnvflg, kb, kmax, kbcon, ktcon, fscav,
414 & xmb, c0t, eta, zi, xlamue, xlamud, delp,
417 use machine , only : kind_phys
418 use physcons, only : g => con_g, qamin
423 integer, intent(in) :: im,
424 & ix, km, itc, ntc, ntr
425 real(kind=kind_phys), intent(in) :: delt
427 logical, dimension(im), intent(in) :: cnvflg
428! integer, dimension(im), intent(in) :: jmin,
429 integer, dimension(im), intent(in) ::
430 & kb, kmax, kbcon, ktcon
431 real(kind=kind_phys), dimension(im), intent(in) ::
433 real(kind=kind_phys), dimension(ntc), intent(in) :: fscav
434 real(kind=kind_phys), dimension(im,km), intent(in) :: c0t,
435 & eta, zi, xlamue !, xlamud
436 real(kind=kind_phys), dimension(ix,km), intent(in) :: delp
437 real(kind=kind_phys), dimension(ix,km,ntr+2), intent(in) :: qtr
439 real(kind=kind_phys), dimension(im,km,ntc), intent(out) :: qaero
442c -- general variables
443 integer :: i, indx, it, k, kk, km1, kp1, n
444! real(kind=kind_phys) :: adw, aup, dtime_max, dv1q, dv2q, dv3q,
445 real(kind=kind_phys) :: aup, dtime_max, dv1q, dv2q, dv3q,
446 & dtovdz, dz, factor, ptem, ptem1, qamax, tem, tem1
447 real(kind=kind_phys), dimension(ix,km) :: xmbp
448c -- chemical transport variables
449 real(kind=kind_phys), dimension(im,km,ntc) :: ctro2,ecko2,dellae2
450c -- additional variables for tracers for wet deposition,
451 real(kind=kind_phys), dimension(im,km,ntc) :: chem_c, chem_pw,
453c -- if reevaporation is enabled, uncomment lines below
454c real(kind=kind_phys), dimension(im,ntc) :: pwav
455c real(kind=kind_phys), dimension(im,km) :: pwdper
456c real(kind=kind_phys), dimension(im,km,ntr) :: chem_pwd
457c -- additional variables for fct
458 real(kind=kind_phys), dimension(im,km) :: flx_lo, totlout, clipout
460 real(kind=kind_phys), parameter :: one = 1.0_kind_phys
461 real(kind=kind_phys), parameter :: half = 0.5_kind_phys
462 real(kind=kind_phys), parameter :: quarter = 0.25_kind_phys
463 real(kind=kind_phys), parameter :: zero = 0.0_kind_phys
464 real(kind=kind_phys), parameter :: epsil = 1.e-22_kind_phys ! prevent division by zero
465 real(kind=kind_phys), parameter :: escav = 0.8_kind_phys ! wet scavenging efficiency
469c -- check if aerosols are present
470.or..or.
if ( ntc <= 0 itc <= 0 ntr <= 0 ) return
471 if ( ntr < itc + ntc - 3 ) return
473c -- initialize work variables
491 if (k <= kmax(i)) qaero(i,k,n) = max(qamin, qtr(i,k,it))
498 if (cnvflg(i)) xmbp(i,k) = g * xmb(i) / delp(i,k)
507 if (kp1 <= kmax(i)) ctro2(i,k,n) =
508 & half * (qaero(i,k,n) + qaero(i,kp1,n))
513 ctro2(i,kmax(i),n) = qaero(i,kmax(i),n)
520.and.
if (cnvflg(i) (k <= kb(i)))
521 & ecko2(i,k,n) = ctro2(i,k,n)
528 ! if (cnvflg(i)) ecdo2(i,jmin(i),n) = ctro2(i,jmin(i),n)
532c do chemical tracers, first need to know how much reevaporates
534c aerosol re-evaporation is set to zero for now
535c uncomment and edit the following code to enable re-evaporation
541c pwdper(i,k)= -edto(i)*pwdo(i,k)/pwavo(i)
545c calculate include mixing ratio (ecko2), how much goes into
546c rainwater to be rained out (chem_pw), and total scavenged,
547c if not reevaporated (pwav)
554.and.
if ((k > kb(i)) (k < ktcon(i))) then
555 dz = zi(i,k) - zi(i,kk)
556 tem = half * (xlamue(i,k)+xlamue(i,kk)) * dz
557! tem1 = quarter * (xlamud(i,k)+xlamud(i,kk)) * dz
558 tem1 = quarter * (xlamud(i )+xlamud(i )) * dz
559 factor = one + tem - tem1
561c if conserved (not scavenging) then
562 ecko2(i,k,n) = ((one-tem1)*ecko2(i,kk,n)
563 & + half*tem*(ctro2(i,k,n)+ctro2(i,kk,n)))/factor
565c how much will be scavenged
567c this choice was used in GF, and is also described in a
568c successful implementation into CESM in GRL (Yu et al. 2019),
569c it uses dimesnsionless scavenging coefficients (fscav),
570c but includes henry coeffs with gas phase chemistry
572c fraction fscav is going into liquid
573 chem_c(i,k,n)=escav*fscav(n)*ecko2(i,k,n)
575c of that part is going into rain out (chem_pw)
576 tem=chem_c(i,k,n)/(one+c0t(i,k)*dz)
577 chem_pw(i,k,n)=c0t(i,k)*dz*tem*eta(i,kk) !etah
578 ecko2(i,k,n)=tem+ecko2(i,k,n)-chem_c(i,k,n)
580c pwav needed fo reevaporation in downdraft
581c if including reevaporation, please uncomment code below
582c pwav(i,n)=pwav(i,n)+chem_pw(i,k,n)
589 if (k >= ktcon(i)) ecko2(i,k,n)=ctro2(i,k,n)
594c reevaporation of some, pw and pwd terms needed later for dellae2
600.and.
! if (cnvflg(i) (k < jmin(i))) then
601! dz = zi(i,kp1) - zi(i,k)
602! if (k >= kbcon(i)) then
604! tem1 = half * xlamdd * dz
607! tem1 = half * (xlamd(i)+xlamdd) * dz
609! factor = one + tem - tem1
610! ecdo2(i,k,n) = ((one-tem1)*ecdo2(i,kp1,n)
611! & +half*tem*(ctro2(i,k,n)+ctro2(i,kp1,n)))/factor
612c if including reevaporation, please uncomment code below
613c ecdo2(i,k,n)=ecdo2(i,k,n)+pwdper(i,kp1)*pwav(i,n)
614c chem_pwd(i,k,n)=max(zero,pwdper(i,kp1)*pwav(i,n))
623c subsidence term treated in fct routine
624! dellae2(i,1,n) = edto(i)*etad(i,1)*ecdo2(i,1,n)*xmbp(i,1)
634c for the subsidence term already is considered
635 dellae2(i,k,n) = eta(i,kk) * ecko2(i,kk,n) * xmbp(i,k)
640c --- for updraft & downdraft vertical transport
642c initialize maximum allowed timestep for upstream difference approach
648 if (kk < ktcon(i)) dtime_max = min(dtime_max,half*delp(i,kk))
652c now for every chemistry tracer
657.and.
if (cnvflg(i) (k < ktcon(i))) then
658 dz = zi(i,k) - zi(i,kk)
660 if (k <= kb(i)) aup = zero
662! if (k > jmin(i)) adw = zero
664 dv1q = half * (ecko2(i,k,n) + ecko2(i,kk,n))
665 dv2q = half * (ctro2(i,k,n) + ctro2(i,kk,n))
666c dv3q = half * (ecdo2(i,k,n) + ecdo2(i,kk,n))
668 tem = half * (xlamue(i,k) + xlamue(i,kk))
669 !tem1 = half * (xlamud(i,k) + xlamud(i,kk))
670 tem1 = half * (xlamud(i ) + xlamud(i ))
672! if (k <= kbcon(i)) then
674! ptem1 = xlamd(i) + xlamdd
679 dellae2(i,k,n) = dellae2(i,k,n) +
680c detrainment from updraft
681 & ( aup*tem1*eta(i,kk)*dv1q
682c entrainement into up and downdraft
683! & - (aup*tem*eta(i,kk)+adw*edto(i)*ptem*etad(i,k))*dv2q
684 & - (aup*tem*eta(i,kk))*dv2q
685c detrainment from downdraft
686! & + (adw*edto(i)*ptem1*etad(i,k)*dv3q)
689 wet_dep(i,k,n)=chem_pw(i,k,n)*g/delp(i,k)
691c sinks from where updraft and downdraft start
692! if (k == jmin(i)+1) then
693! dellae2(i,k,n) = dellae2(i,k,n)
694! & -edto(i)*etad(i,kk)*ctro2(i,kk,n)*xmbp(i,k)
697 dellae2(i,k,n) = dellae2(i,k,n)
698 & -eta(i,k)*ctro2(i,k,n)*xmbp(i,k)
708 dellae2(i,k,n) = dellae2(i,k,n)
709 & -eta(i,k)*ctro2(i,k,n)*xmbp(i,k)
722c compute low-order mass flux, upstream
726.and.
if (cnvflg(i) (kk < ktcon(i))) then
728 if (kk >= kb(i) ) tem = eta(i,kk)
729! if (kk <= jmin(i)) tem = tem - edto(i)*etad(i,kk)
730c low-order flux,upstream
732 flx_lo(i,k) = -xmb(i) * tem * qaero(i,k,n)
733 elseif (tem < zero) then
734 flx_lo(i,k) = -xmb(i) * tem * qaero(i,kk,n)
740c --- make sure low-ord fluxes don't violate positive-definiteness
744 if (cnvflg(i) .and. (k <= ktcon(i)))
then
745c time step / grid spacing
746 dtovdz = g * dtime_max / abs(delp(i,k))
748 totlout(i,k)=max(zero,flx_lo(i,kp1))-min(zero,flx_lo(i,k))
749 clipout(i,k)=min(one ,qaero(i,k,n)/max(epsil,totlout(i,k))
750 & / (1.0001_kind_phys*dtovdz))
755c recompute upstream mass fluxes
759 if (cnvflg(i) .and. (kk < ktcon(i)))
then
761 if (kk >= kb(i) ) tem = eta(i,kk)
764 flx_lo(i,k) = flx_lo(i,k) * clipout(i,k)
765 elseif (tem < zero)
then
766 flx_lo(i,k) = flx_lo(i,k) * clipout(i,kk)
772c --- a positive-definite low-order(diffusive) solution for the subsidnce fluxes
776 if (cnvflg(i) .and. (k <= ktcon(i)))
then
777 dtovdz = g * dtime_max / abs(delp(i,k))
778 dellae2(i,k,n) = dellae2(i,k,n)
779 & -(flx_lo(i,kp1)-flx_lo(i,k))*dtovdz/dtime_max
786c convert wet deposition to total mass deposited over dt and dp
791 if (cnvflg(i) .and. (k < ktcon(i)))
792 & wet_dep(i,k,n) = wet_dep(i,k,n)*xmb(i)*delt*delp(i,k)
797c compute final aerosol concentrations
802 if (cnvflg(i) .and. (k <= min(kmax(i),ktcon(i))))
then
803 qaero(i,k,n) = qaero(i,k,n) + dellae2(i,k,n) * delt
804 if (qaero(i,k,n) < zero)
then
805c add negative mass to wet deposition
806 wet_dep(i,k,n) = wet_dep(i,k,n)-qaero(i,k,n)*delp(i,k)