CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
samfaerosols.F
1
4
5 implicit none
6
7 private
8
9 public :: samfdeepcnv_aerosols, samfshalcnv_aerosols
10
11 contains
12
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,
16 & qtr, qaero)
17
18 use machine , only : kind_phys
19 use physcons, only : g => con_g, qamin
20
21 implicit none
22
23c -- input arguments
24 integer, intent(in) :: im,
25 & ix, km, itc, ntc, ntr
26 real(kind=kind_phys), intent(in) :: delt,
27 & xlamde, xlamdd
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,
32 & xlamd, xmb
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
38c -- output arguments
39 real(kind=kind_phys), dimension(im,km,ntc), intent(out) :: qaero
40
41c -- local variables
42c -- general variables
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,
49 & dellae2
50c -- additional variables for tracers for wet deposition,
51 real(kind=kind_phys), dimension(im,km,ntc) :: chem_c, chem_pw,
52 & wet_dep
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
59
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 ! prevent division by zero
65
66c -- begin
67
68c -- check if aerosols are present
69 if ( ntc <= 0 .or. itc <= 0 .or. ntr <= 0 ) return
70 if ( ntr < itc + ntc - 3 ) return
71
72c -- initialize work variables
73 km1 = km - 1
74
75 chem_c = zero
76 chem_pw = zero
77 ctro2 = zero
78 dellae2 = zero
79 ecdo2 = zero
80 ecko2 = zero
81 qaero = zero
82 wet_dep = zero
83
84c -- set work arrays
85
86 do n = 1, ntc
87 it = n + itc - 1
88 do k = 1, km
89 do i = 1, im
90 if (k <= kmax(i)) qaero(i,k,n) = max(qamin, qtr(i,k,it))
91 enddo
92 enddo
93 enddo
94
95 do k = 1, km
96 do i = 1, im
97 if (cnvflg(i)) xmbp(i,k) = g * xmb(i) / delp(i,k)
98 enddo
99 enddo
100
101 do n = 1, ntc
102c -- interface level
103 do k = 1, km1
104 kp1 = k + 1
105 do i = 1, im
106 if (kp1 <= kmax(i)) ctro2(i,k,n) =
107 & half * (qaero(i,k,n) + qaero(i,kp1,n))
108 enddo
109 enddo
110c -- top level
111 do i = 1, im
112 ctro2(i,kmax(i),n) = qaero(i,kmax(i),n)
113 enddo
114 enddo
115
116 do n = 1, ntc
117 do k = 1, km
118 do i = 1, im
119 if (cnvflg(i) .and. (k <= kb(i)))
120 & ecko2(i,k,n) = ctro2(i,k,n)
121 enddo
122 enddo
123 enddo
124
125 do n = 1, ntc
126 do i = 1, im
127 if (cnvflg(i)) ecdo2(i,jmin(i),n) = ctro2(i,jmin(i),n)
128 enddo
129 enddo
130
131c do chemical tracers, first need to know how much reevaporates
132
133c aerosol re-evaporation is set to zero for now
134c uncomment and edit the following code to enable re-evaporation
135c chem_pwd = zero
136c pwdper = zero
137c pwav = zero
138c do i = 1, im
139c do k=1,jmin(i)
140c pwdper(i,k)= -edto(i)*pwdo(i,k)/pwavo(i)
141c enddo
142c enddo
143c
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)
147
148 do n = 1, ntc
149 do k = 2, km1
150 kk = k - 1
151 do i = 1, im
152 if (cnvflg(i)) then
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
158
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
162
163c how much will be scavenged
164c
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
169
170c fraction fscav is going into liquid
171 chem_c(i,k,n)=fscav(n)*ecko2(i,k,n)
172
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) !etah
176 ecko2(i,k,n)=tem+ecko2(i,k,n)-chem_c(i,k,n)
177
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)
181 endif
182 endif
183 enddo
184 enddo
185 do k = 1, km1
186 do i = 1, im
187 if (k >= ktcon(i)) ecko2(i,k,n)=ctro2(i,k,n)
188 enddo
189 enddo
190 enddo
191
192c reevaporation of some, pw and pwd terms needed later for dellae2
193
194 do n = 1, ntc
195 do k = km1, 1, -1
196 kp1 = k + 1
197 do i = 1, im
198 if (cnvflg(i) .and. (k < jmin(i))) then
199 dz = zi(i,kp1) - zi(i,k)
200 if (k >= kd94(i)) then
201 tem = xlamde * dz
202 tem1 = half * xlamdd * dz
203 else
204 tem = xlamde * dz
205 tem1 = half * (xlamd(i)+xlamdd) * dz
206 endif
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))
213 endif
214 enddo
215 enddo
216 enddo
217
218 do n = 1, ntc
219 do i = 1, im
220 if (cnvflg(i)) then
221c subsidence term treated in fct routine
222 dellae2(i,1,n) = edto(i)*etad(i,1)*ecdo2(i,1,n)*xmbp(i,1)
223 endif
224 enddo
225 enddo
226
227 do n = 1, ntc
228 do i = 1, im
229 if (cnvflg(i)) then
230 k = ktcon(i)
231 kk = k - 1
232c for the subsidence term already is considered
233 dellae2(i,k,n) = eta(i,kk) * ecko2(i,kk,n) * xmbp(i,k)
234 endif
235 enddo
236 enddo
237
238c --- for updraft & downdraft vertical transport
239c
240c initialize maximum allowed timestep for upstream difference approach
241c
242 dtime_max=delt
243 do k=2,km1
244 kk = k - 1
245 do i = 1, im
246 if (kk < ktcon(i)) dtime_max = min(dtime_max,half*delp(i,kk))
247 enddo
248 enddo
249
250c now for every chemistry tracer
251 do n = 1, ntc
252 do k = 2, km1
253 kk = k - 1
254 do i = 1, im
255 if (cnvflg(i) .and. (k < ktcon(i))) then
256 dz = zi(i,k) - zi(i,kk)
257 aup = one
258 if (k <= kb(i)) aup = zero
259 adw = one
260 if (k > jmin(i)) adw = zero
261
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))
265
266 tem = half * (xlamue(i,k) + xlamue(i,kk))
267 tem1 = half * (xlamud(i,k) + xlamud(i,kk))
268
269 if (k <= kd94(i)) then
270 ptem = xlamde
271 ptem1 = xlamd(i) + xlamdd
272 else
273 ptem = xlamde
274 ptem1 = xlamdd
275 endif
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)
283
284 wet_dep(i,k,n)=chem_pw(i,k,n)*g/delp(i,k)
285
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)
290 endif
291 if (k == kb(i))then
292 dellae2(i,k,n) = dellae2(i,k,n)
293 & -eta(i,k)*ctro2(i,k,n)*xmbp(i,k)
294 endif
295 endif
296 enddo
297 enddo
298
299 do i = 1, im
300 if (cnvflg(i)) then
301 if (kb(i) == 1) then
302 k=kb(i)
303 dellae2(i,k,n) = dellae2(i,k,n)
304 & -eta(i,k)*ctro2(i,k,n)*xmbp(i,k)
305 endif
306 endif
307 enddo
308
309 enddo
310
311c for every tracer...
312
313 do n = 1, ntc
314 flx_lo = zero
315 totlout = zero
316 clipout = zero
317c compute low-order mass flux, upstream
318 do k = 2, km1
319 kk = k - 1
320 do i = 1, im
321 if (cnvflg(i) .and. (kk < ktcon(i))) then
322 tem = zero
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
326 if (tem > zero) then
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)
330 endif
331 endif
332 enddo
333 enddo
334
335c --- make sure low-ord fluxes don't violate positive-definiteness
336 do k=1,km1
337 kp1 = k + 1
338 do i=1,im
339.and. if (cnvflg(i) (k <= ktcon(i))) then
340c time step / grid spacing
341 dtovdz = g * dtime_max / abs(delp(i,k))
342c total flux out
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))
346 endif
347 enddo
348 enddo
349
350c recompute upstream mass fluxes
351 do k = 2, km1
352 kk = k - 1
353 do i = 1, im
354.and. if (cnvflg(i) (kk < ktcon(i))) then
355 tem = zero
356 if (kk >= kb(i) ) tem = eta(i,kk)
357 if (kk <= jmin(i)) tem = tem - edto(i)*etad(i,kk)
358 if (tem > zero) then
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)
362 endif
363 endif
364 enddo
365 enddo
366
367c --- a positive-definite low-order (diffusive) solution for the subsidnce fluxes
368 do k=1,km1
369 kp1 = k + 1
370 do i=1,im
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
375 endif
376 enddo
377 enddo
378
379 enddo ! ctr
380
381c convert wet deposition to total mass deposited over dt and dp
382
383 do n = 1, ntc
384 do k = 1, km
385 do i = 1, im
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)
388 enddo
389 enddo
390 enddo
391
392c compute final aerosol concentrations
393
394 do n = 1, ntc
395 do k = 1, km
396 do i = 1, im
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)
402 qaero(i,k,n) = qamin
403 endif
404 endif
405 enddo
406 enddo
407 enddo
408
409 return
410 end subroutine samfdeepcnv_aerosols
411
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,
415 & qtr, qaero)
416
417 use machine , only : kind_phys
418 use physcons, only : g => con_g, qamin
419
420 implicit none
421
422c -- input arguments
423 integer, intent(in) :: im,
424 & ix, km, itc, ntc, ntr
425 real(kind=kind_phys), intent(in) :: delt
426! & xlamde, xlamdd
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) ::
432 & xmb, xlamud
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
438c -- output arguments
439 real(kind=kind_phys), dimension(im,km,ntc), intent(out) :: qaero
440
441c -- local variables
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,
452 & wet_dep
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
459
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
466
467c -- begin
468
469c -- check if aerosols are present
470.or..or. if ( ntc <= 0 itc <= 0 ntr <= 0 ) return
471 if ( ntr < itc + ntc - 3 ) return
472
473c -- initialize work variables
474 km1 = km - 1
475
476 chem_c = zero
477 chem_pw = zero
478 ctro2 = zero
479 dellae2 = zero
480 !ecdo2 = zero
481 ecko2 = zero
482 qaero = zero
483 wet_dep = zero
484
485c -- set work arrays
486
487 do n = 1, ntc
488 it = n + itc - 1
489 do k = 1, km
490 do i = 1, im
491 if (k <= kmax(i)) qaero(i,k,n) = max(qamin, qtr(i,k,it))
492 enddo
493 enddo
494 enddo
495
496 do k = 1, km
497 do i = 1, im
498 if (cnvflg(i)) xmbp(i,k) = g * xmb(i) / delp(i,k)
499 enddo
500 enddo
501
502 do n = 1, ntc
503c -- interface level
504 do k = 1, km1
505 kp1 = k + 1
506 do i = 1, im
507 if (kp1 <= kmax(i)) ctro2(i,k,n) =
508 & half * (qaero(i,k,n) + qaero(i,kp1,n))
509 enddo
510 enddo
511c -- top level
512 do i = 1, im
513 ctro2(i,kmax(i),n) = qaero(i,kmax(i),n)
514 enddo
515 enddo
516
517 do n = 1, ntc
518 do k = 1, km
519 do i = 1, im
520.and. if (cnvflg(i) (k <= kb(i)))
521 & ecko2(i,k,n) = ctro2(i,k,n)
522 enddo
523 enddo
524 enddo
525
526 !do n = 1, ntc
527 ! do i = 1, im
528 ! if (cnvflg(i)) ecdo2(i,jmin(i),n) = ctro2(i,jmin(i),n)
529 ! enddo
530 !enddo
531
532c do chemical tracers, first need to know how much reevaporates
533
534c aerosol re-evaporation is set to zero for now
535c uncomment and edit the following code to enable re-evaporation
536c chem_pwd = zero
537c pwdper = zero
538c pwav = zero
539c do i = 1, im
540c do k=1,jmin(i)
541c pwdper(i,k)= -edto(i)*pwdo(i,k)/pwavo(i)
542c enddo
543c enddo
544c
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)
548
549 do n = 1, ntc
550 do k = 2, km1
551 kk = k - 1
552 do i = 1, im
553 if (cnvflg(i)) then
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
560
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
564
565c how much will be scavenged
566c
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
571
572c fraction fscav is going into liquid
573 chem_c(i,k,n)=escav*fscav(n)*ecko2(i,k,n)
574
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)
579
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)
583 endif
584 endif
585 enddo
586 enddo
587 do k = 1, km1
588 do i = 1, im
589 if (k >= ktcon(i)) ecko2(i,k,n)=ctro2(i,k,n)
590 enddo
591 enddo
592 enddo
593
594c reevaporation of some, pw and pwd terms needed later for dellae2
595
596! do n = 1, ntc
597! do k = km1, 1, -1
598! kp1 = k + 1
599! do i = 1, im
600.and.! if (cnvflg(i) (k < jmin(i))) then
601! dz = zi(i,kp1) - zi(i,k)
602! if (k >= kbcon(i)) then
603! tem = xlamde * dz
604! tem1 = half * xlamdd * dz
605! else
606! tem = xlamde * dz
607! tem1 = half * (xlamd(i)+xlamdd) * dz
608! endif
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))
615! endif
616! enddo
617! enddo
618! enddo
619
620! do n = 1, ntc
621! do i = 1, im
622! if (cnvflg(i)) then
623c subsidence term treated in fct routine
624! dellae2(i,1,n) = edto(i)*etad(i,1)*ecdo2(i,1,n)*xmbp(i,1)
625! endif
626! enddo
627! enddo
628
629 do n = 1, ntc
630 do i = 1, im
631 if (cnvflg(i)) then
632 k = ktcon(i)
633 kk = k - 1
634c for the subsidence term already is considered
635 dellae2(i,k,n) = eta(i,kk) * ecko2(i,kk,n) * xmbp(i,k)
636 endif
637 enddo
638 enddo
639
640c --- for updraft & downdraft vertical transport
641c
642c initialize maximum allowed timestep for upstream difference approach
643c
644 dtime_max=delt
645 do k=2,km1
646 kk = k - 1
647 do i = 1, im
648 if (kk < ktcon(i)) dtime_max = min(dtime_max,half*delp(i,kk))
649 enddo
650 enddo
651
652c now for every chemistry tracer
653 do n = 1, ntc
654 do k = 2, km1
655 kk = k - 1
656 do i = 1, im
657.and. if (cnvflg(i) (k < ktcon(i))) then
658 dz = zi(i,k) - zi(i,kk)
659 aup = one
660 if (k <= kb(i)) aup = zero
661! adw = one
662! if (k > jmin(i)) adw = zero
663
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))
667
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 ))
671
672! if (k <= kbcon(i)) then
673! ptem = xlamde
674! ptem1 = xlamd(i) + xlamdd
675! else
676! ptem = xlamde
677! ptem1 = xlamdd
678! endif
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)
687 & ) * dz * xmbp(i,k)
688
689 wet_dep(i,k,n)=chem_pw(i,k,n)*g/delp(i,k)
690
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)
695! endif
696 if (k == kb(i))then
697 dellae2(i,k,n) = dellae2(i,k,n)
698 & -eta(i,k)*ctro2(i,k,n)*xmbp(i,k)
699 endif
700 endif
701 enddo
702 enddo
703
704 do i = 1, im
705 if (cnvflg(i)) then
706 if (kb(i) == 1) then
707 k=kb(i)
708 dellae2(i,k,n) = dellae2(i,k,n)
709 & -eta(i,k)*ctro2(i,k,n)*xmbp(i,k)
710 endif
711 endif
712 enddo
713
714 enddo
715
716c for every tracer...
717
718 do n = 1, ntc
719 flx_lo = zero
720 totlout = zero
721 clipout = zero
722c compute low-order mass flux, upstream
723 do k = 2, km1
724 kk = k - 1
725 do i = 1, im
726.and. if (cnvflg(i) (kk < ktcon(i))) then
727 tem = zero
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
731 if (tem > zero) then
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)
735 endif
736 endif
737 enddo
738 enddo
739
740c --- make sure low-ord fluxes don't violate positive-definiteness
741 do k=1,km1
742 kp1 = k + 1
743 do i=1,im
744 if (cnvflg(i) .and. (k <= ktcon(i))) then
745c time step / grid spacing
746 dtovdz = g * dtime_max / abs(delp(i,k))
747c total flux out
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))
751 endif
752 enddo
753 enddo
754
755c recompute upstream mass fluxes
756 do k = 2, km1
757 kk = k - 1
758 do i = 1, im
759 if (cnvflg(i) .and. (kk < ktcon(i))) then
760 tem = zero
761 if (kk >= kb(i) ) tem = eta(i,kk)
762! if (kk <= jmin(i)) tem = tem - edto(i)*etad(i,kk)
763 if (tem > zero) then
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)
767 endif
768 endif
769 enddo
770 enddo
771
772c --- a positive-definite low-order(diffusive) solution for the subsidnce fluxes
773 do k=1,km1
774 kp1 = k + 1
775 do i=1,im
776 if (cnvflg(i) .and. (k <= ktcon(i))) then
777 dtovdz = g * dtime_max / abs(delp(i,k)) ! time step /grid spacing
778 dellae2(i,k,n) = dellae2(i,k,n)
779 & -(flx_lo(i,kp1)-flx_lo(i,k))*dtovdz/dtime_max
780 endif
781 enddo
782 enddo
783
784 enddo ! ctr
785
786c convert wet deposition to total mass deposited over dt and dp
787
788 do n = 1, ntc
789 do k = 1, km
790 do i = 1, im
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)
793 enddo
794 enddo
795 enddo
796
797c compute final aerosol concentrations
798
799 do n = 1, ntc
800 do k = 1, km
801 do i = 1, im
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)
807 qaero(i,k,n) = qamin
808 endif
809 endif
810 enddo
811 enddo
812 enddo
813
814 return
815 end subroutine samfshalcnv_aerosols
816
817 end module samfcnv_aerosols