| blob | 1dd6bf164e6a0a8f6007a27e7778c73d454b56f9 |
1 !WRF:ODEL_LAYER:PHYSICS
2 !
4 MODULE module_ctrans_grell
5 USE module_cu_gd
6 USE module_dep_simple
7 USE module_state_description, only:p_co,p_qv,p_so2,p_hno3,p_hno4,p_n2o5,p_nh3,p_h2o2, &
8 p_o3,p_ora1,p_op1,p_paa,p_sulf,p_so4aj,p_nh4aj,p_no3aj, &
9 p_bc1,p_bc2,p_oc1,p_oc2,p_seas_1,p_seas_2, &
10 p_seas_3,p_seas_4,p_dms
12 REAL, PARAMETER :: qcldwtr_cutoff = 1.0e-6 ! kg/kg
14 REAL, PARAMETER :: mwdry = 28.966 ! Molecular mass of dry air (g/mol)
15 REAL, PARAMETER :: mwso4 = 96.00 ! Molecular mass of SO4-- (g/mol)
16 REAL, PARAMETER :: mwno3 = 62.0 ! Molecular mass of NO3- (g/mol)
17 REAL, PARAMETER :: mwnh4 = 18.0985 ! Molecular mass of NH4+ (g/mol)
19 CONTAINS
21 !-------------------------------------------------------------
22 SUBROUTINE GRELLDRVCT(DT,itimestep,DX, &
23 rho_phy,RAINCV,chem, &
24 U,V,t_phy,moist,dz8w,p_phy, &
25 XLV,CP,G,r_v,z,cu_co_ten,ktop_out, &
26 wd_no3,wd_so4, &
27 k22_shallow,kbcon_shallow,ktop_shallow,xmb_shallow, &
28 ishallow,num_moist,numgas,num_chem,chemopt,scalaropt, &
29 conv_tr_wetscav,conv_tr_aqchem, &
30 ids,ide, jds,jde, kds,kde, &
31 ims,ime, jms,jme, kms,kme, &
32 its,ite, jts,jte, kts,kte )
33 ! USE module_configure
34 ! USE module_state_description
35 !-------------------------------------------------------------
36 IMPLICIT NONE
37 !-------------------------------------------------------------
38 INTEGER, INTENT(IN ) :: &
39 numgas,chemopt,scalaropt, &
40 ids,ide, jds,jde, kds,kde, &
41 ims,ime, jms,jme, kms,kme, &
42 its,ite, jts,jte, kts,kte, &
43 ishallow,num_chem,num_moist, &
44 conv_tr_wetscav,conv_tr_aqchem
46 INTEGER, INTENT(IN ) :: ITIMESTEP
48 REAL, INTENT(IN ) :: XLV, R_v
49 REAL, INTENT(IN ) :: CP,G
50 REAL, DIMENSION( ims:ime , kms:kme , jms:jme,num_moist ) , &
51 INTENT(IN ) :: moist
52 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
53 INTENT(IN ) :: &
54 U, &
55 V, &
56 t_phy, &
57 z, &
58 p_phy, &
59 dz8w, &
60 rho_phy
61 !
62 ! on output for control only, purely diagnostic
63 !
64 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
65 INTENT(INOUT ) :: &
66 cu_co_ten
68 INTEGER, DIMENSION( ims:ime , jms:jme ) , INTENT( OUT) :: ktop_out
70 !
71 REAL, INTENT(IN ) :: DT, DX
72 !
73 REAL, DIMENSION( ims:ime , kms:kme , jms:jme, num_chem ), &
74 INTENT(INOUT) :: &
75 chem
77 REAL, DIMENSION( ims:ime , jms:jme ), &
78 INTENT(IN) :: RAINCV, xmb_shallow
79 INTEGER, DIMENSION( ims:ime , jms:jme ), &
80 INTENT(IN) :: k22_shallow,kbcon_shallow,ktop_shallow
81 !
82 ! Accumulated wet deposition
83 !
84 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT) :: wd_no3,wd_so4
85 ! LOCAL VARS
86 real, dimension (its:ite,kts:kte) :: &
87 OUTT,OUTQ,OUTQC
88 real, dimension (its:ite) :: &
89 pret, ter11
90 ! Auxiliary wet deposition variables
91 REAL, DIMENSION (its:ite,num_chem) :: wetdep_1d
92 REAL, DIMENSION (ims:ime,jms:jme,num_chem) :: wetdep_2d
93 ! Wet deposition over the current time step
94 REAL, DIMENSION( ims:ime , jms:jme ) :: wdi_no3,wdi_so4
95 !
96 ! basic environmental input includes moisture convergence (mconv)
97 ! omega (omeg), windspeed (us,vs), and a flag (aaeq) to turn off
98 ! convection for this call only and at that particular gridpoint
99 !
100 real, dimension (its:ite,kts:kte) :: &
101 T,TN,q,qo,PO,P,US,VS,hstary
102 real, dimension (its:ite,kts:kte,num_chem) :: &
103 tracer,tracert,tracert3
104 real, dimension (its:ite) :: &
105 Z1,PSUR,AAEQ,xmb3
106 integer, dimension (its:ite) :: &
107 ktop,k23,kbcon3,ktop3
109 INTEGER :: nv,i,j,k,ICLDCK,ipr,jpr,npr
110 REAL :: tcrit,dp,dq,epsilc
111 INTEGER :: itf,jtf,ktf,iopt
112 epsilc=1.e-30
113 wetdep_1d(:,:) = 0.0
114 wetdep_2d(:,:,:) = 0.0
115 ! return
116 ! ipr=111
117 ! jpr=40
118 ! if(itimestep.lt.34.or.itimestep.gt.36)ipr=0
119 ! if(itimestep.lt.34.or.itimestep.gt.36)jpr=0
120 ! ipr=61
121 ! jpr=60
122 ipr=0
123 jpr=0
124 npr=1
125 if(p_co.gt.1)npr=p_co
126 tcrit=258.
127 iopt=0
128 itf=MIN(ite,ide-1)
129 ktf=MIN(kte,kde-1)
130 jtf=MIN(jte,jde-1)
131 !
132 !
133 123 continue
134 DO 100 J = jts,jtf
135 if(j.eq.jpr)print *,'dt = ',dt
136 DO I=ITS,ITF
137 ktop(i)=0
138 PSUR(I)=p_phy(I,kts,J)*.01
139 TER11(I)=z(i,kts,j)
140 aaeq(i)=0.
142 !
143 ! rainrate is input for this transport/wet-deposition routine
144 !
145 pret(i)=raincv(i,j)/dt
146 if(pret(i).le.0.)aaeq(i)=20.
147 ENDDO
148 if(Ishallow.eq.1)then
149 DO I=ITS,ITF
150 xmb3(i)=xmb_shallow(i,j)
151 ktop3(i)=ktop_shallow(i,j)
152 k23(i)=k22_shallow(i,j)
153 kbcon3(i)=kbcon_shallow(i,j)
154 ENDDO
155 endif
156 DO K=kts,ktf
157 DO I=ITS,ITF
158 po(i,k)=p_phy(i,k,j)*.01
159 P(I,K)=PO(i,k)
160 US(I,K) =u(i,k,j)
161 VS(I,K) =v(i,k,j)
162 T(I,K)=t_phy(i,k,j)
163 q(I,K)=moist(i,k,j,p_qv)
164 IF(Q(I,K).LT.1.E-08)Q(I,K)=1.E-08
165 ENDDO
166 ENDDO
167 do nv=2,num_chem
168 DO K=kts,ktf
169 DO I=ITS,ITF
170 tracer(i,k,nv)=max(epsilc,chem(i,k,j,nv))
171 tracert(i,k,nv)=0.
172 tracert3(i,k,nv)=0.
173 ENDDO
174 ENDDO
175 ENDDO
176 DO K=kts,ktf
177 DO I=ITS,ITF
178 cu_co_ten(i,k,j)=0.
179 ! hstary(i,k)=hstar4(nv)*exp(dhr(nv)*(1./t(i,k)-1./298.))
180 if(i.eq.ipr.and.j.eq.jpr)then
181 print *,k,pret(i),tracer(i,k,npr),p(i,k),z(i,k,j)
182 endif
183 ENDDO
184 ENDDO
185 ! ENDDO
186 !
187 !---- CALL NON_RESOLVED CONVECTIVE TRANSPORT
188 !
189 CALL CUP_ct(ktop,k23,kbcon3,ktop3,xmb3, &
190 tracer,j,AAEQ,T,Q,TER11,PRET,P,tracert,tracert3, &
191 hstary,DT,PSUR,US,VS,tcrit,wetdep_1d, &
192 xlv,r_v,cp,g,ipr,jpr,npr,num_chem,chemopt,scalaropt,&
193 conv_tr_wetscav,conv_tr_aqchem, &
194 ishallow,numgas,ids,ide, jds,jde, kds,kde, &
195 ims,ime, jms,jme, kms,kme, &
196 its,ite, jts,jte, kts,kte )
198 do nv=2,num_chem
199 DO I=its,itf
200 if(pret(i).le.0.)then
201 DO K=kts,ktf
202 tracert(i,k,nv)=0.
203 ENDDO
204 endif
205 enddo
206 enddo
207 if(ishallow.eq.1)then
208 CALL neg_check_ct('shallow',pret,ktop3,epsilc,dt,tracer,tracert3,iopt,num_chem, &
209 its,ite,kts,kte,itf,ktf,ipr,jpr,npr,j)
210 do nv=2,num_chem
211 DO I=its,itf
212 DO K=kts,ktf
213 chem(i,k,j,nv)=max(epsilc,chem(i,k,j,nv)+tracert3(i,k,nv)*dt)
214 enddo
215 ENDDO
216 ENDDO
217 endif
218 !
219 ! now deep convection
220 !
221 CALL neg_check_ct('deep',pret,ktop,epsilc,dt,tracer,tracert,iopt,num_chem, &
222 its,ite,kts,kte,itf,ktf,ipr,jpr,npr,j)
223 do nv=2,num_chem
224 DO I=its,itf
225 if(pret(i).gt.0.)then
226 DO K=kts,ktf
227 chem(i,k,j,nv)=max(epsilc,chem(i,k,j,nv)+tracert(i,k,nv)*dt)
228 if(nv.eq.npr)then
229 cu_co_ten(i,k,j)=tracert(i,k,npr)*dt
230 ! if(i.eq.ipr.and.j.eq.jpr)print *,k,chem(i,k,j,nv),cu_co_ten(i,k,j)
231 endif
232 ENDDO
233 else
234 DO K=kts,ktf
235 tracert(i,k,nv)=0.
236 if(nv.eq.npr)cu_co_ten(i,k,j)=0.
237 enddo
238 endif
239 ENDDO
240 ENDDO
242 wetdep_2d(its:ite,j,:) = wetdep_1d(its:ite,:)
243 ktop_out(its:itf,J) = ktop(its:itf)
245 100 continue
247 ! Calculate the wet deposition over the time step:
249 wdi_no3(:,:) = 0.0
250 wdi_so4(:,:) = 0.0
252 ! We use the indices of the chem array that point to aerosol outside of cloud water,
253 ! because that's what the cumulus scheme operates with.
255 if (p_no3aj .gt.1) wdi_no3(its:ite,jts:jte) = wdi_no3(its:ite,jts:jte) + wetdep_2d(its:ite,jts:jte,p_no3aj)*dt*0.001/mwno3 ! mmol/m2
256 if (p_hno3 .gt.1) wdi_no3(its:ite,jts:jte) = wdi_no3(its:ite,jts:jte) + wetdep_2d(its:ite,jts:jte,p_hno3)*dt ! mmol/m2
257 if (p_hno4 .gt.1) wdi_no3(its:ite,jts:jte) = wdi_no3(its:ite,jts:jte) + wetdep_2d(its:ite,jts:jte,p_hno4)*dt ! mmol/m2
259 if (p_so4aj .gt.1) wdi_so4(its:ite,jts:jte) = wdi_so4(its:ite,jts:jte) + wetdep_2d(its:ite,jts:jte,p_so4aj)*dt*0.001/mwso4 ! mmol/m2
260 if (p_sulf .gt.1) wdi_so4(its:ite,jts:jte) = wdi_so4(its:ite,jts:jte) + wetdep_2d(its:ite,jts:jte,p_sulf)*dt ! mmol/m2
261 if (p_so2 .gt.1) wdi_so4(its:ite,jts:jte) = wdi_so4(its:ite,jts:jte) + wetdep_2d(its:ite,jts:jte,p_so2)*dt ! mmol/m2
263 ! Update the accumulated wet deposition:
265 wd_no3(its:ite,jts:jte) = wd_no3(its:ite,jts:jte) + wdi_no3(its:ite,jts:jte) ! mmol/m2
266 wd_so4(its:ite,jts:jte) = wd_so4(its:ite,jts:jte) + wdi_so4(its:ite,jts:jte) ! mmol/m2
268 END SUBROUTINE GRELLDRVCT
271 SUBROUTINE CUP_ct(ktop,k23,kbcon3,ktop3,xmb3,tracer,J,AAEQ,T,Q,Z1, &
272 PRE,P,tracert,tracert3,hstary,DTIME,PSUR,US,VS,TCRIT, &
273 wetdep,xl,rv,cp,g,ipr,jpr,npr,num_chem,chemopt,scalaropt,&
274 conv_tr_wetscav,conv_tr_aqchem, &
275 ishallow,numgas,ids,ide, jds,jde, kds,kde, &
276 ims,ime, jms,jme, kms,kme, &
277 its,ite, jts,jte, kts,kte )
279 IMPLICIT NONE
281 integer &
282 ,intent (in ) :: &
283 num_chem,ids,ide, jds,jde, kds,kde, &
284 ims,ime, jms,jme, kms,kme,scalaropt,conv_tr_wetscav, &
285 conv_tr_aqchem,ishallow, &
286 its,ite, jts,jte, kts,kte,ipr,jpr,npr,chemopt,numgas
287 integer, intent (in ) :: &
288 j
289 !
290 !
291 !
292 !tracert = output temp tendency (per s)
293 ! pre = input precip
294 real, dimension (its:ite,kts:kte,num_chem) &
295 ,intent (inout ) :: &
296 tracert,tracer,tracert3
297 real, dimension (its:ite) &
298 ,intent (inout ) :: &
299 pre
300 integer, dimension (its:ite) &
301 ,intent (inout ) :: &
302 ktop,k23,kbcon3,ktop3
303 integer, dimension (its:ite) :: &
304 kbcon
305 !
306 ! basic environmental input includes moisture convergence (mconv)
307 ! omega (omeg), windspeed (us,vs), and a flag (aaeq) to turn off
308 ! convection for this call only and at that particular gridpoint
309 !
310 real, dimension (its:ite,kts:kte) &
311 ,intent (in ) :: &
312 T,P,US,VS,HSTARY
313 real, dimension (its:ite,kts:kte) &
314 ,intent (inout) :: &
315 Q
316 real, dimension (its:ite) &
317 ,intent (in ) :: &
318 Z1,PSUR,AAEQ,xmb3
321 real &
322 ,intent (in ) :: &
323 dtime,tcrit,xl,cp,rv,g
326 real, dimension (its:ite,1:3) :: &
327 edtc
328 !
329 !
330 !
331 !***************** the following are your basic environmental
332 ! variables. They carry a "_cup" if they are
333 ! on model cloud levels (staggered). They carry
334 ! an "o"-ending (z becomes zo), if they are the forced
335 ! variables. They are preceded by x (z becomes xz)
336 ! to indicate modification by some typ of cloud
337 !
338 ! z = heights of model levels
339 ! q = environmental mixing ratio
340 ! qes = environmental saturation mixing ratio
341 ! t = environmental temp
342 ! p = environmental pressure
343 ! he = environmental moist static energy
344 ! hes = environmental saturation moist static energy
345 ! z_cup = heights of model cloud levels
346 ! q_cup = environmental q on model cloud levels
347 ! qes_cup = saturation q on model cloud levels
348 ! t_cup = temperature (Kelvin) on model cloud levels
349 ! p_cup = environmental pressure
350 ! he_cup = moist static energy on model cloud levels
351 ! hes_cup = saturation moist static energy on model cloud levels
352 ! gamma_cup = gamma on model cloud levels
353 !
354 !
355 ! hcd = moist static energy in downdraft
356 ! zd normalized downdraft mass flux
357 ! dby = buoancy term
358 ! entr = entrainment rate
359 ! zd = downdraft normalized mass flux
360 ! entr= entrainment rate
361 ! hcd = h in model cloud
362 ! bu = buoancy term
363 ! zd = normalized downdraft mass flux
364 ! gamma_cup = gamma on model cloud levels
365 ! mentr_rate = entrainment rate
366 ! qcd = cloud q (including liquid water) after entrainment
367 ! qrch = saturation q in cloud
368 ! pwd = evaporate at that level
369 ! pwev = total normalized integrated evaoprate (I2)
370 ! entr= entrainment rate
371 ! z1 = terrain elevation
372 ! entr = downdraft entrainment rate
373 ! jmin = downdraft originating level
374 ! kdet = level above ground where downdraft start detraining
375 ! psur = surface pressure
376 ! z1 = terrain elevation
377 ! zd = downdraft normalized mass flux
378 ! zu = updraft normalized mass flux
379 ! mbdt = arbitrary numerical parameter
380 ! dtime = dt over which forcing is applied
381 ! kbcon = LFC of parcel from k22
382 ! k22 = updraft originating level
383 ! dby = buoancy term
384 ! ktop = cloud top (output)
385 ! xmb = total base mass flux
386 ! hc = cloud moist static energy
387 ! hkb = moist static energy at originating level
388 ! mentr_rate = entrainment rate
390 real, dimension (its:ite,kts:kte) :: &
391 he,hes,qes,z,pwdper, &
393 qes_cup,q_cup,he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup, &
395 dby,qc,qrcd,pwd,pw,hcd,qcd,dbyd,hc,qrc,zu,zd,clw_all,zd3, &
396 clw_all3,cd3,pw3,zu3, &
398 ! cd = detrainment function for updraft
399 ! cdd = detrainment function for downdraft
401 cd,cdd,cdd3,scr1,DELLAH,DELLAQ,DELLAT,DELLAQC
403 ! edt = epsilon
404 ! edt = epsilon
405 real, dimension (its:ite) :: &
406 edt,HKB,QKB,edt3, &
407 XMB,PWAV,PWEV,BU,cap_max,cap_max_increment
408 real, dimension (its:ite,kts:kte,num_chem) :: &
409 tr3_c,tr3_up,tr3_pw
410 real, dimension (its:ite,num_chem) :: &
411 trkb3
412 real, dimension (its:ite,kts:kte,num_chem) :: &
413 tr_c,tr_up,tr_dd,tr_dd3,tre_cup,tr_pw,tr_pwd
414 real, dimension (its:ite,num_chem) :: &
415 trkb,wetdep
416 integer, dimension (its:ite) :: &
417 kzdown,KDET,K22,KB,JMIN,kstabi,kstabm, & !-lxz
418 ierr,KBMAX,ierr5
420 integer :: &
421 ki,I,K,KK
422 real :: &
423 day,dz,mbdt,entr_rate,radius,entrd_rate,mentr_rate,mentrd_rate, &
424 zcutdown,edtmax,edtmin,depth_min,zkbmax,z_detr,zktop, &
425 dh,cap_maxs,entr_rate3,mentr_rate3
427 integer :: itf,jtf,ktf
429 itf=MIN(ite,ide-1)
430 ktf=MIN(kte,kde-1)
431 jtf=MIN(jte,jde-1)
433 !sms$distribute end
434 day=86400.
435 !
436 !--- specify entrainmentrate and detrainmentrate
437 !
438 radius=12000.
439 !
440 !--- gross entrainment rate (these may be changed later on in the
441 !--- program, depending what your detrainment is!!)
442 !
443 entr_rate=.2/radius
444 entr_rate3=.2/200.
445 !
446 !--- entrainment of mass
447 !
448 mentrd_rate=0.
449 mentr_rate=entr_rate
450 mentr_rate3=entr_rate3
451 !
452 !--- initial detrainmentrates
453 !
454 do k=kts,ktf
455 do i=its,itf
456 cd(i,k)=0.1*entr_rate
457 cd3(i,k)=entr_rate3
458 cdd(i,k)=0.
459 cdd3(i,k)=0.
460 clw_all(i,k)=0.
461 enddo
462 enddo
463 do ki=1,num_chem
464 do k=kts,ktf
465 do i=its,itf
466 tr_dd3(i,k,ki)=0.
467 enddo
468 enddo
469 enddo
470 !
471 !--- max/min allowed value for epsilon (ratio downdraft base mass flux/updraft
472 ! base mass flux
473 !
474 edtmax=.8
475 edtmin=.2
476 !
477 !--- minimum depth (m), clouds must have
478 !
479 depth_min=500.
480 !
481 !--- maximum depth (mb) of capping
482 !--- inversion (larger cap = no convection)
483 !
484 cap_maxs=175.
485 !sms$to_local(grid_dh: <1, mix :size>, <2, mjx :size>) begin
486 DO 7 i=its,itf
487 kbmax(i)=1
488 cap_max_increment(i)=0.
489 edt(i)=0.
490 edt3(i)=0.
491 kstabm(i)=ktf-1
492 IERR(i)=0
493 IERR5(i)=0
494 if(ktop3(i).lt.2)ierr5(i)=20
495 if(xmb3(i).eq.0)ierr5(i)=21
496 7 CONTINUE
497 do i=its,itf
498 cap_max(i)=cap_maxs
499 do k=kts,kte
500 zd3(i,k)=0.
501 enddo
502 enddo
503 !
504 !--- max height(m) above ground where updraft air can originate
505 !
506 zkbmax=4000.
507 !
508 !--- height(m) above which no downdrafts are allowed to originate
509 !
510 zcutdown=3000.
511 !
512 !--- depth(m) over which downdraft detrains all its mass
513 !
514 z_detr=1250.
515 !
516 mbdt=dtime*4.E-03
517 !
518 !--- calculate moist static energy, heights, qes
519 !
520 call cup_env(z,qes,he,hes,t,q,p,z1, &
521 psur,ierr,tcrit,0,xl,cp, &
522 ids,ide, jds,jde, kds,kde, &
523 ims,ime, jms,jme, kms,kme, &
524 its,ite, jts,jte, kts,kte)
525 !
526 !--- environmental values on cloud levels
527 !
528 call cup_env_clev(t,qes,q,he,hes,z,p,qes_cup,q_cup,he_cup, &
529 hes_cup,z_cup,p_cup,gamma_cup,t_cup,psur, &
530 ierr,z1,xl,rv,cp, &
531 ids,ide, jds,jde, kds,kde, &
532 ims,ime, jms,jme, kms,kme, &
533 its,ite, jts,jte, kts,kte)
534 call cup_env_clev_tr(tracer,tre_cup,num_chem,ierr, &
535 ids,ide, jds,jde, kds,kde, &
536 ims,ime, jms,jme, kms,kme, &
537 its,ite, jts,jte, kts,kte)
538 do i=its,itf
539 if(ierr(i).eq.0)then
540 !
541 do k=kts,ktf-2
542 if(z_cup(i,k).gt.zkbmax+z1(i))then
543 kbmax(i)=k
544 go to 25
545 endif
546 enddo
547 25 continue
548 !
549 !
550 !--- level where detrainment for downdraft starts
551 !
552 do k=kts,ktf
553 if(z_cup(i,k).gt.z_detr+z1(i))then
554 kdet(i)=k
555 go to 26
556 endif
557 enddo
558 26 continue
559 !
560 endif
561 enddo
562 !
563 !
564 !
565 !------- DETERMINE LEVEL WITH HIGHEST MOIST STATIC ENERGY CONTENT - K22
566 !
567 CALL cup_MAXIMI(HE_CUP,3,KBMAX,K22,ierr, &
568 ids,ide, jds,jde, kds,kde, &
569 ims,ime, jms,jme, kms,kme, &
570 its,ite, jts,jte, kts,kte)
571 DO 36 i=its,itf
572 IF(ierr(I).eq.0.)THEN
573 IF(K22(I).GE.KBMAX(i))ierr(i)=2
574 endif
575 36 CONTINUE
576 !
577 ! done with basic stuff that is needed everywhere, from here on do not do
578 ! deep convection where aaeq .ne.0
579 !
580 DO i=its,itf
581 if(aaeq(i).ne.0.)then
582 ierr(i)=20
583 endif
584 enddo
585 !
586 !--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON
587 !
588 call cup_kbcon(cap_max_increment,1,k22,kbcon,he_cup,hes_cup, &
589 ierr,kbmax,p_cup,cap_max, &
590 ids,ide, jds,jde, kds,kde, &
591 ims,ime, jms,jme, kms,kme, &
592 its,ite, jts,jte, kts,kte)
593 !
594 !--- increase detrainment in stable layers
595 !
596 CALL cup_minimi(HEs_cup,Kbcon,kstabm,kstabi,ierr, &
597 ids,ide, jds,jde, kds,kde, &
598 ims,ime, jms,jme, kms,kme, &
599 its,ite, jts,jte, kts,kte)
600 do i=its,itf
601 IF(ierr(I).eq.0.)THEN
602 if(kstabm(i)-1.gt.kstabi(i))then
603 do k=kstabi(i),kstabm(i)-1
604 cd(i,k)=cd(i,k-1)+1.5*entr_rate
605 if(cd(i,k).gt.10.0*entr_rate)cd(i,k)=10.0*entr_rate
606 enddo
607 ENDIF
608 ENDIF
609 ENDDO
610 !
611 !--- calculate incloud moist static energy
612 !
613 call cup_up_he(k22,hkb,z_cup,cd,mentr_rate,he_cup,hc, &
614 kbcon,ierr,dby,he,hes_cup, &
615 ids,ide, jds,jde, kds,kde, &
616 ims,ime, jms,jme, kms,kme, &
617 its,ite, jts,jte, kts,kte)
619 !--- DETERMINE CLOUD TOP - KTOP
620 !
621 call cup_ktop(1,dby,kbcon,ktop,ierr, &
622 ids,ide, jds,jde, kds,kde, &
623 ims,ime, jms,jme, kms,kme, &
624 its,ite, jts,jte, kts,kte)
625 DO 37 i=its,itf
626 kzdown(i)=0
627 if(ierr(i).eq.0)then
628 zktop=(z_cup(i,ktop(i))-z1(i))*.6
629 zktop=min(zktop+z1(i),zcutdown+z1(i))
630 do k=kts,ktf
631 if(z_cup(i,k).gt.zktop)then
632 kzdown(i)=k
633 go to 37
634 endif
635 enddo
636 endif
637 37 CONTINUE
638 !
639 !--- DOWNDRAFT ORIGINATING LEVEL - JMIN
640 !
641 call cup_minimi(HEs_cup,K22,kzdown,JMIN,ierr, &
642 ids,ide, jds,jde, kds,kde, &
643 ims,ime, jms,jme, kms,kme, &
644 its,ite, jts,jte, kts,kte)
645 DO 100 i=its,ite
646 IF(ierr(I).eq.0.)THEN
647 if(jmin(i).le.3)then
648 ierr(i)=9
649 go to 100
650 endif
651 !
652 !--- check whether it would have buoyancy, if there where
653 !--- no entrainment/detrainment
654 !
655 101 continue
656 if(jmin(i)-1.lt.KDET(I))kdet(i)=jmin(i)-1
657 if(jmin(i).ge.Ktop(I)-1)jmin(i)=ktop(i)-2
658 ki=jmin(i)
659 hcd(i,ki)=hes_cup(i,ki)
660 DZ=Z_cup(i,Ki+1)-Z_cup(i,Ki)
661 dh=dz*(HCD(i,Ki)-hes_cup(i,ki))
662 dh=0.
663 !
664 do k=ki-1,1,-1
665 hcd(i,k)=hes_cup(i,jmin(i))
666 DZ=Z_cup(i,K+1)-Z_cup(i,K)
667 dh=dh+dz*(HCD(i,K)-hes_cup(i,k))
668 if(dh.gt.0.)then
669 jmin(i)=jmin(i)-1
670 if(jmin(i).gt.3)then
671 go to 101
672 else if(jmin(i).le.3)then
673 ierr(i)=9
674 go to 100
675 endif
676 endif
677 enddo
679 IF(JMIN(I).LE.3)then
680 ierr(i)=4
681 endif
683 ENDIF
684 100 continue
685 !
686 ! - Must have at least depth_min m between cloud convective base
687 ! and cloud top.
688 !
689 do i=its,itf
690 IF(ierr(I).eq.0.)THEN
691 if(jmin(i).le.3)then
692 ierr(i)=9
693 endif
694 IF(-z_cup(I,KBCON(I))+z_cup(I,KTOP(I)).LT.depth_min)then
695 ierr(i)=6
696 endif
697 endif
698 enddo
700 !
701 !c--- normalized updraft mass flux profile
702 !
703 call cup_up_nms(zu,z_cup,mentr_rate,cd,kbcon,ktop,ierr,k22, &
704 ids,ide, jds,jde, kds,kde, &
705 ims,ime, jms,jme, kms,kme, &
706 its,ite, jts,jte, kts,kte)
707 !
708 !c--- normalized downdraft mass flux profile,also work on bottom detrainment
709 !--- in this routine
710 !
711 call cup_dd_nms(zd,z_cup,cdd,mentrd_rate,jmin,ierr, &
712 0,kdet,z1, &
713 ids,ide, jds,jde, kds,kde, &
714 ims,ime, jms,jme, kms,kme, &
715 its,ite, jts,jte, kts,kte)
716 !
717 !--- downdraft moist static energy
718 !
719 call cup_dd_he(hes_cup,zd,hcd,z_cup,cdd,mentrd_rate, &
720 jmin,ierr,he,dbyd,he_cup, &
721 ids,ide, jds,jde, kds,kde, &
722 ims,ime, jms,jme, kms,kme, &
723 its,ite, jts,jte, kts,kte)
724 !
725 !--- calculate moisture properties of downdraft
726 !
727 call cup_dd_moisture(zd,hcd,hes_cup,qcd,qes_cup, &
728 pwd,q_cup,z_cup,cdd,mentrd_rate,jmin,ierr,gamma_cup, &
729 pwev,bu,qrcd,q,he,t_cup,1,xl, &
730 ids,ide, jds,jde, kds,kde, &
731 ims,ime, jms,jme, kms,kme, &
732 its,ite, jts,jte, kts,kte)
733 !
734 !--- calculate moisture properties of updraft
735 !
736 call cup_up_moisture(ierr,z_cup,qc,qrc,pw,pwav, &
737 kbcon,ktop,cd,dby,mentr_rate,clw_all, &
738 q,GAMMA_cup,zu,qes_cup,k22,q_cup,xl, &
739 ids,ide, jds,jde, kds,kde, &
740 ims,ime, jms,jme, kms,kme, &
741 its,ite, jts,jte, kts,kte)
742 !
743 !--- DETERMINE DOWNDRAFT STRENGTH IN TERMS OF WINDSHEAR
744 !
745 call cup_dd_edt(ierr,us,vs,z,ktop,kbcon,edt,p,pwav, &
746 pwev,edtmax,edtmin,3,edtc, &
747 ids,ide, jds,jde, kds,kde, &
748 ims,ime, jms,jme, kms,kme, &
749 its,ite, jts,jte, kts,kte)
750 do i=its,itf
751 if(ierr(i).eq.0)then
752 edt(i)=edtc(i,2)
753 endif
754 enddo
755 !
756 ! massflux from precip and normalized cloud properties
757 !
758 pwdper=0.
759 do i=its,itf
761 if(ierr(i).gt.0)pre(i)=0.
762 if(ierr(i).eq.0)then
763 xmb(i)=pre(i)/(pwav(i)+edt(i)*pwev(i))
764 !
765 !--- percent of that that is evaporated (pwd is negative)
766 !
767 if(i.eq.ipr.and.j.eq.jpr)then
768 print *,'xmb,edt,pwav = ',xmb(i),edt(i),pwav(i)
769 print *,'k,pwdper(i,k),pw,pwd(i,k)',z1(i)
770 endif
771 do k=1,ktop(i)
772 pwdper(i,k)=-edt(i)*pwd(i,k)/pwav(i)
773 if(i.eq.ipr.and.j.eq.jpr)then
774 print *,k,pwdper(i,k),pw(i,k),pwd(i,k)
775 endif
776 enddo
777 endif
778 enddo
779 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
780 !
781 !!!!! NOW WE HAVE EVREYTHING TO CALCULATE TRACER TRANSPORT AND WET DEPOSITION !!!
782 !
783 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
784 !
785 !--- calculate incloud tracer distribution
786 !
787 ! if(j.eq.jpr)print *,'calling up_tracer'
788 call cup_up_tracer(ierr,tcrit,t,pre,z_cup,p,tracer,tre_cup,tr_up,tr_pw, &
789 tr_c,hstary,pw,clw_all,kbcon,ktop,cd,mentr_rate,zu,k22,&
790 numgas,chemopt,scalaropt,conv_tr_wetscav,conv_tr_aqchem, &
791 ids,ide, jds,jde, kds,kde, &
792 ims,ime, jms,jme, kms,kme, &
793 its,ite, jts,jte, kts,kte, &
794 ipr,jpr,j,npr,num_chem,'deep')
795 ! if(j.eq.jpr)print *,'called up_tracer'
796 !
797 ! for shallow convection, only 3 subroutines required
798 !
799 if(ishallow.eq.1)then
800 call cup_up_nms(zu3,z_cup,mentr_rate3,cd3,kbcon3,ktop3,ierr5,k23, &
801 ids,ide, jds,jde, kds,kde, &
802 ims,ime, jms,jme, kms,kme, &
803 its,ite, jts,jte, kts,kte)
804 call cup_up_tracer(ierr5,tcrit,t,pre,z_cup,p,tracer,tre_cup,tr3_up,tr3_pw, &
805 tr3_c,hstary,pw3,clw_all3,kbcon3,ktop3,cd3,mentr_rate3,zu3,k23,&
806 numgas,chemopt,scalaropt,conv_tr_wetscav,conv_tr_aqchem, &
807 ids,ide, jds,jde, kds,kde, &
808 ims,ime, jms,jme, kms,kme, &
809 its,ite, jts,jte, kts,kte, &
810 ipr,jpr,j,npr,num_chem,'shallow')
811 call cup_dellas_tr(ierr5,z_cup,p_cup,tr_dd3,edt3,zd3,cdd3, &
812 tracer,tracert3,j,mentrd_rate,zu3,g,xmb3, &
813 cd3,tr3_up,ktop3,k23,kbcon3,mentr_rate3,jmin,tre_cup,kdet, &
814 k23,ipr,jpr,npr,'shallow',num_chem, &
815 ids,ide, jds,jde, kds,kde, &
816 ims,ime, jms,jme, kms,kme, &
817 its,ite, jts,jte, kts,kte )
818 endif
820 call cup_dd_tracer(ierr,z_cup,qrcd,tracer,tre_cup,tr_up,tr_dd, &
821 tr_pw,tr_pwd,jmin,cdd,mentrd_rate,zd,pwdper,wetdep,xmb,k22, &
822 numgas,ids,ide, jds,jde, kds,kde, &
823 ims,ime, jms,jme, kms,kme, &
824 its,ite, jts,jte, kts,kte)
825 if(j.eq.jpr)print *,'called dd_tracer'
828 !
829 if(j.eq.jpr)then
830 i=ipr
831 print *,'in 250 loop ',edt(ipr),ierr(ipr)
832 ! if(ierr(i).eq.0.or.ierr(i).eq.3)then
833 print *,k22(I),kbcon(i),ktop(i),jmin(i)
834 print *,edt(i)
835 do k=kts,ktf
836 print *,k,z(i,k),he(i,k),hes(i,k)
837 enddo
838 do k=1,ktop(i)+1
839 print *,zu(i,k),zd(i,k),pw(i,k),pwd(i,k)
840 enddo
841 print *,'tr_up(i,k,6),tr_dd(i,k,6),tr_pw(i,k,6),tr_pwd(i,k,6)'
842 do k=1,ktop(i)+1
843 print *,tr_up(i,k,npr),tr_dd(i,k,npr),tr_pw(i,k,npr),tr_pwd(i,k,npr)
844 enddo
845 endif
846 ! endif
847 !
848 !--- calculate transport tendencies
849 !
850 !--- 1. in bottom layer
851 !
852 call cup_dellabot_tr(ipr,jpr,tre_cup,ierr,z_cup,p,tr_dd,edt, &
853 zd,cdd,tracer,tracert,j,mentrd_rate,z,g,xmb, &
854 num_chem,ids,ide, jds,jde, kds,kde, &
855 ims,ime, jms,jme, kms,kme, &
856 its,ite, jts,jte, kts,kte)
857 !
858 !--- 2. everywhere else
859 !
861 call cup_dellas_tr(ierr,z_cup,p_cup,tr_dd,edt,zd,cdd, &
862 tracer,tracert,j,mentrd_rate,zu,g,xmb, &
863 cd,tr_up,ktop,k22,kbcon,mentr_rate,jmin,tre_cup,kdet, &
864 k22,ipr,jpr,npr,'deep',num_chem, &
865 ids,ide, jds,jde, kds,kde, &
866 ims,ime, jms,jme, kms,kme, &
867 its,ite, jts,jte, kts,kte )
868 if(j.eq.jpr)then
869 i=ipr
870 do k=kts,ktf
871 print *,k,tracer(i,k,npr),tracert(i,k,npr)
872 enddo
873 endif
874 !
875 ! may need more below for wet deposition......
876 !
877 !
878 ! call cup_output_wd ( &
879 ! ids,ide, jds,jde, kds,kde, &
880 ! ims,ime, jms,jme, kms,kme, &
881 ! its,ite, jts,jte, kts,kte)
883 END SUBROUTINE CUP_CT
885 SUBROUTINE cup_dellabot_tr(ipr,jpr,tre_cup,ierr,z_cup,p_cup, &
886 tr_dd,edt,zd,cdd,tracer,tracert,j,mentrd_rate,z,g,xmb, &
887 num_chem,ids,ide, jds,jde, kds,kde, &
888 ims,ime, jms,jme, kms,kme, &
889 its,ite, jts,jte, kts,kte )
891 IMPLICIT NONE
893 integer &
894 ,intent (in ) :: &
895 num_chem,ids,ide, jds,jde, kds,kde, &
896 ims,ime, jms,jme, kms,kme, &
897 its,ite, jts,jte, kts,kte
898 integer, intent (in ) :: &
899 j,ipr,jpr
900 !
901 ! ierr error value, maybe modified in this routine
902 !
903 real, dimension (its:ite,kts:kte,1:num_chem) &
904 ,intent (out ) :: &
905 tracert
906 real, dimension (its:ite,kts:kte,1:num_chem) &
907 ,intent (in ) :: &
908 tre_cup,tracer,tr_dd
909 real, dimension (its:ite,kts:kte) &
910 ,intent (in ) :: &
911 z_cup,p_cup,zd,cdd,z
912 real, dimension (its:ite) &
913 ,intent (in ) :: &
914 edt,xmb
915 real &
916 ,intent (in ) :: &
917 g,mentrd_rate
918 integer, dimension (its:ite) &
919 ,intent (inout) :: &
920 ierr
921 !
922 ! local variables in this routine
923 !
925 integer i
926 real detdo,detdo1,detdo2,entdo,dp,dz,subin, &
927 totmas
928 !
929 integer :: itf, ktf, nv, npr
930 npr=24
931 itf=MIN(ite,ide-1)
932 ktf=MIN(kte,kde-1)
933 !
934 !
935 if(j.eq.jpr)print *,'in cup dellabot '
936 tracert=0.
937 do 100 i=its,itf
938 if(ierr(i).ne.0)go to 100
939 dz=z_cup(i,2)-z_cup(i,1)
940 DP=100.*(p_cup(i,1)-P_cup(i,2))
941 detdo1=edt(i)*zd(i,2)*CDD(i,1)*DZ
942 detdo2=edt(i)*zd(i,1)
943 entdo=edt(i)*zd(i,2)*mentrd_rate*dz
944 subin=-EDT(I)*zd(i,2)
945 detdo=detdo1+detdo2-entdo+subin
946 do nv=2,num_chem
947 tracert(I,1,nv)=(detdo1*.5*(tr_dd(i,1,nv)+tr_dd(i,2,nv)) &
948 +detdo2*tr_dd(i,1,nv) &
949 +subin*tre_cup(i,2,nv) &
950 -entdo*tracer(i,1,nv))*g/dp*xmb(i)
951 enddo
952 if(j.eq.jpr.and.i.eq.ipr)print *,'in cup dellabot ',tracert(I,1,npr), &
953 detdo1,detdo2,subin,entdo,tr_dd(i,1,npr),tr_dd(i,2,npr),tracer(i,1,npr)
954 100 CONTINUE
956 END SUBROUTINE cup_dellabot_tr
959 SUBROUTINE cup_dellas_tr(ierr,z_cup,p_cup,tr_dd,edt,zd,cdd, &
960 tracer,tracert,j,mentrd_rate,zu,g,xmb, &
961 cd,tr_up,ktop,k22,kbcon,mentr_rate,jmin,tre_cup,kdet,kpbl, &
962 ipr,jpr,npr,name,num_chem, &
963 ids,ide, jds,jde, kds,kde, &
964 ims,ime, jms,jme, kms,kme, &
965 its,ite, jts,jte, kts,kte )
967 IMPLICIT NONE
969 integer &
970 ,intent (in ) :: &
971 num_chem,ids,ide, jds,jde, kds,kde, &
972 ims,ime, jms,jme, kms,kme, &
973 its,ite, jts,jte, kts,kte
974 integer, intent (in ) :: &
975 j,ipr,jpr,npr
976 !
977 ! ierr error value, maybe modified in this routine
978 !
979 real, dimension (its:ite,kts:kte,1:num_chem) &
980 ,intent (inout ) :: &
981 tracert
982 real, dimension (its:ite,kts:kte,1:num_chem) &
983 ,intent (in ) :: &
984 tr_up,tr_dd,tre_cup,tracer
985 real, dimension (its:ite,kts:kte) &
986 ,intent (in ) :: &
987 z_cup,p_cup,zd,cdd,cd,zu
988 real, dimension (its:ite) &
989 ,intent (in ) :: &
990 edt,xmb
991 real &
992 ,intent (in ) :: &
993 g,mentrd_rate,mentr_rate
994 integer, dimension (its:ite) &
995 ,intent (in ) :: &
996 kbcon,ktop,k22,jmin,kdet,kpbl
997 integer, dimension (its:ite) &
998 ,intent (inout) :: &
999 ierr
1000 character *(*), intent (in) :: &
1001 name
1002 !
1003 ! local variables in this routine
1004 !
1006 integer i,k,nv
1007 real detdo1,detdo2,entdo,dp,dz,subin,detdo,entup, &
1008 detup,subdown,entdoj,entupk,detupk,totmas
1009 !
1010 integer :: itf, ktf, kstart
1011 ! npr=24
1012 itf=MIN(ite,ide-1)
1013 ktf=MIN(kte,kde-1)
1014 kstart=kts+1
1015 if(name.eq.'shallow')kstart=kts
1016 !
1017 !
1018 i=ipr
1019 if(j.eq.jpr)then
1020 print *,'in dellas kpbl(i),k22(i),kbcon(i),ktop(i),jmin(i)'
1021 print *,kpbl(i),k22(i),kbcon(i),ktop(i),jmin(i)
1022 endif
1023 do nv=2,num_chem
1024 DO K=kstart,kte
1025 do i=its,itf
1026 tracert(i,k,nv)=0.
1027 enddo
1028 enddo
1029 enddo
1030 !
1031 DO 100 k=kts+1,ktf-1
1032 DO 100 i=its,ite
1033 IF(ierr(i).ne.0)GO TO 100
1034 IF(K.Gt.KTOP(I))GO TO 100
1035 !
1036 !--- SPECIFY DETRAINMENT OF DOWNDRAFT, HAS TO BE CONSISTENT
1037 !--- WITH ZD CALCULATIONS IN SOUNDD.
1038 !
1039 DZ=Z_cup(I,K+1)-Z_cup(I,K)
1040 detdo=edt(i)*CDD(i,K)*DZ*ZD(i,k+1)
1041 entdo=edt(i)*mentrd_rate*dz*zd(i,k+1)
1042 subin=zu(i,k+1)-zd(i,k+1)*edt(i)
1043 entup=0.
1044 detup=0.
1045 if(k.ge.kbcon(i).and.k.lt.ktop(i))then
1046 entup=mentr_rate*dz*zu(i,k)
1047 detup=CD(i,K+1)*DZ*ZU(i,k)
1048 endif
1049 subdown=(zu(i,k)-zd(i,k)*edt(i))
1050 entdoj=0.
1051 entupk=0.
1052 detupk=0.
1053 !
1054 if(k.eq.jmin(i))then
1055 entdoj=edt(i)*zd(i,k)
1056 endif
1058 if(k.eq.k22(i)-1)then
1059 entupk=zu(i,kpbl(i))
1060 endif
1062 if(k.gt.kdet(i))then
1063 detdo=0.
1064 endif
1066 if(k.eq.ktop(i)-0)then
1067 detupk=zu(i,ktop(i))
1068 subin=0.
1069 endif
1070 if(k.lt.kbcon(i))then
1071 detup=0.
1072 endif
1073 !C
1074 !C--- CHANGED DUE TO SUBSIDENCE AND ENTRAINMENT
1075 !C
1076 totmas=subin-subdown+detup-entup-entdo+ &
1077 detdo-entupk-entdoj+detupk
1078 if(j.eq.jpr.and.i.eq.ipr)print *,'k,totmas,sui,sud = ',k, &
1079 totmas,subin,subdown
1080 ! if(j.eq.jpr.and.i.eq.ipr)print *,'updr stuff = ',detup,
1081 ! 1 entup,entupk,detupk
1082 ! if(j.eq.jpr.and.i.eq.ipr)print *,'dddr stuff = ',entdo,
1083 ! 1 detdo,entdoj
1084 if(abs(totmas).gt.1.e-6)then
1085 print *,'*********************',i,j,k,totmas,name
1086 print *,kpbl(i),k22(i),kbcon(i),ktop(i)
1087 !c print *,'updr stuff = ',subin,
1088 !c 1 subdown,detup,entup,entupk,detupk
1089 !c print *,'dddr stuff = ',entdo,
1090 !c 1 detdo,entdoj
1091 CALL wrf_error_fatal ( 'cup_dellas_tr: TOTMAS > CRITICAL VALUE')
1092 endif
1093 dp=100.*(p_cup(i,k-1)-p_cup(i,k))
1094 do nv=2,num_chem
1095 ! tracert(i,k,nv)=(subin*tre_cup(i,k+1,nv) &
1096 ! -subdown*tre_cup(i,k,nv) &
1097 tracert(i,k,nv)=(subin*tracer(i,k+1,nv) &
1098 -subdown*tracer(i,k,nv) &
1099 +detup*.5*(tr_up(i,K+1,nv)+tr_up(i,K,nv)) &
1100 +detdo*.5*(tr_dd(i,K+1,nv)+tr_dd(i,K,nv)) &
1101 -entup*tracer(i,k,nv) &
1102 -entdo*tracer(i,k,nv) &
1103 -entupk*tre_cup(i,k22(i),nv) &
1104 -entdoj*tre_cup(i,jmin(i),nv) &
1105 +detupk*tr_up(i,ktop(i),nv) &
1106 )*g/dp*xmb(i)
1107 enddo
1108 if(i.eq.ipr.and.j.eq.jpr)then
1109 print *,k,tracert(i,k,npr),subin*tre_cup(i,k+1,npr),subdown*tre_cup(i,k,npr), &
1110 detdo*.5*(tr_dd(i,K+1,npr)+tr_dd(i,K,npr))
1111 print *,k,detup*.5*(tr_up(i,K+1,npr)+tr_up(i,K,npr)),detupk*tr_up(i,ktop(i),npr), &
1112 entup*tracer(i,k,npr),entdo*tracer(i,k,npr)
1113 print *,k,entupk*tre_cup(i,k,npr),detupk,tr_up(i,ktop(i),npr)
1114 endif
1116 100 CONTINUE
1118 END SUBROUTINE cup_dellas_tr
1119 SUBROUTINE cup_env_clev_tr(tracer,tre_cup,num_chem,ierr, &
1120 ids,ide, jds,jde, kds,kde, &
1121 ims,ime, jms,jme, kms,kme, &
1122 its,ite, jts,jte, kts,kte)
1123 implicit none
1124 integer &
1125 ,intent (in ) :: &
1126 num_chem,ids,ide, jds,jde, kds,kde, &
1127 ims,ime, jms,jme, kms,kme, &
1128 its,ite, jts,jte, kts,kte
1129 integer, dimension (its:ite) &
1130 ,intent (in) :: &
1131 ierr
1133 real, dimension (its:ite,kts:kte,1:num_chem) &
1134 ,intent (in ) :: &
1135 tracer
1136 real, dimension (its:ite,kts:kte,1:num_chem) &
1137 ,intent (out ) :: &
1138 tre_cup
1139 !
1140 ! local variables in this routine
1141 !
1143 integer :: &
1144 i,k,nv,itf,ktf
1145 itf=MIN(ite,ide-1)
1146 ktf=MIN(kte,kde-1)
1147 do nv=2,num_chem
1148 do k=kts+1,ktf
1149 do i=its,ite
1150 if(ierr(i).eq.0)then
1151 tre_cup(i,k,nv)=.5*(tracer(i,k-1,nv)+tracer(i,k,nv))
1152 endif
1153 enddo
1154 enddo
1155 enddo
1156 do nv=2,num_chem
1157 do i=its,ite
1158 if(ierr(i).eq.0)then
1159 tre_cup(i,kts,nv)=tracer(i,kts,nv)
1160 endif
1161 enddo
1162 enddo
1165 END subroutine cup_env_clev_tr
1168 SUBROUTINE cup_up_tracer(ierr,tcrit,t,pre,z_cup,p,tracer,tre_cup,tr_up, &
1169 tr_pw,tr_c,hstary,cupclw,clw_all,kbcon,ktop,cd,mentr_rate,zu,k22, &
1170 numgas,chemopt,scalaropt,conv_tr_wetscav,conv_tr_aqchem, &
1171 ids,ide, jds,jde, kds,kde, &
1172 ims,ime, jms,jme, kms,kme, &
1173 its,ite, jts,jte, kts,kte,ipr,jpr,j,npr,num_chem,name)
1174 ! USE module_configure
1175 USE module_state_description, only: RADM2SORG,RADM2SORG_AQ,RACMSORG_AQ,RACMSORG_KPP, &
1176 RADM2SORG_KPP,RACM_ESRLSORG_KPP,RACM_SOA_VBS_KPP, &
1177 RADM2SORG_AQCHEM,RACMSORG_AQCHEM
1178 USE module_ctrans_aqchem
1179 USE module_input_chem_data, only: get_last_gas
1180 implicit none
1181 ! Aqeuous species pointers INCLUDE File
1183 !...........PARAMETERS and their descriptions:
1185 INTEGER, PARAMETER :: NGAS = 12 ! number of gas-phase species for AQCHEM
1186 INTEGER, PARAMETER :: NAER = 36 ! number of aerosol species for AQCHEM
1187 INTEGER, PARAMETER :: NLIQS = 41 ! number of liquid-phase species in AQCHEM
1189 !...pointers for the AQCHEM array GAS
1191 INTEGER, PARAMETER :: LSO2 = 1 ! Sulfur Dioxide
1192 INTEGER, PARAMETER :: LHNO3 = 2 ! Nitric Acid
1193 INTEGER, PARAMETER :: LN2O5 = 3 ! Dinitrogen Pentoxide
1194 INTEGER, PARAMETER :: LCO2 = 4 ! Carbon Dioxide
1195 INTEGER, PARAMETER :: LNH3 = 5 ! Ammonia
1196 INTEGER, PARAMETER :: LH2O2 = 6 ! Hydrogen Perioxide
1197 INTEGER, PARAMETER :: LO3 = 7 ! Ozone
1198 INTEGER, PARAMETER :: LFOA = 8 ! Formic Acid
1199 INTEGER, PARAMETER :: LMHP = 9 ! Methyl Hydrogen Peroxide
1200 INTEGER, PARAMETER :: LPAA = 10 ! Peroxyacidic Acid
1201 INTEGER, PARAMETER :: LH2SO4 = 11 ! Sulfuric Acid
1202 INTEGER, PARAMETER :: LHCL = 12 ! Hydrogen Chloride
1204 !...pointers for the AQCHEM array AEROSOL
1206 INTEGER, PARAMETER :: LSO4AKN = 1 ! Aitken-mode Sulfate
1207 INTEGER, PARAMETER :: LSO4ACC = 2 ! Accumulation-mode Sulfate
1208 INTEGER, PARAMETER :: LSO4COR = 3 ! Coarse-mode Sulfate
1209 INTEGER, PARAMETER :: LNH4AKN = 4 ! Aitken-mode Ammonium
1210 INTEGER, PARAMETER :: LNH4ACC = 5 ! Accumulation-mode Ammonium
1211 INTEGER, PARAMETER :: LNO3AKN = 6 ! Aitken-mode Nitrate
1212 INTEGER, PARAMETER :: LNO3ACC = 7 ! Accumulation-mode Nitrate
1213 INTEGER, PARAMETER :: LNO3COR = 8 ! Coarse-mode Nitrate
1214 INTEGER, PARAMETER :: LORGAAKN = 9 ! Aitken-mode anthropogenic SOA
1215 INTEGER, PARAMETER :: LORGAACC = 10 ! Accumulation-mode anthropogenic SOA
1216 INTEGER, PARAMETER :: LORGPAKN = 11 ! Aitken-mode primary organic aerosol
1217 INTEGER, PARAMETER :: LORGPACC = 12 ! Accumulation-mode primary organic aerosol
1218 INTEGER, PARAMETER :: LORGBAKN = 13 ! Aitken-mode biogenic SOA
1219 INTEGER, PARAMETER :: LORGBACC = 14 ! Accumulation-mode biogenic SOA
1220 INTEGER, PARAMETER :: LECAKN = 15 ! Aitken-mode elemental carbon
1221 INTEGER, PARAMETER :: LECACC = 16 ! Accumulation-mode elemental carbon
1222 INTEGER, PARAMETER :: LPRIAKN = 17 ! Aitken-mode primary aerosol
1223 INTEGER, PARAMETER :: LPRIACC = 18 ! Accumulation-mode primary aerosol
1224 INTEGER, PARAMETER :: LPRICOR = 19 ! Coarse-mode primary aerosol
1225 INTEGER, PARAMETER :: LNAAKN = 20 ! Aitken-mode Sodium
1226 INTEGER, PARAMETER :: LNAACC = 21 ! Accumulation-mode Sodium
1227 INTEGER, PARAMETER :: LNACOR = 22 ! Coarse-mode Sodium
1228 INTEGER, PARAMETER :: LCLAKN = 23 ! Aitken-mode Chloride ion
1229 INTEGER, PARAMETER :: LCLACC = 24 ! Accumulation-mode Chloride ion
1230 INTEGER, PARAMETER :: LCLCOR = 25 ! Coarse-mode Chloride ion
1231 INTEGER, PARAMETER :: LNUMAKN = 26 ! Aitken-mode number
1232 INTEGER, PARAMETER :: LNUMACC = 27 ! Accumulation-mode number
1233 INTEGER, PARAMETER :: LNUMCOR = 28 ! Coarse-mode number
1234 INTEGER, PARAMETER :: LSRFAKN = 29 ! Aitken-mode surface area
1235 INTEGER, PARAMETER :: LSRFACC = 30 ! Accumulation-mode surface area
1236 INTEGER, PARAMETER :: LNACL = 31 ! Sodium Chloride aerosol for AE3 only {depreciated in AE4}
1237 INTEGER, PARAMETER :: LCACO3 = 32 ! Calcium Carbonate aerosol (place holder)
1238 INTEGER, PARAMETER :: LMGCO3 = 33 ! Magnesium Carbonate aerosol (place holder)
1239 INTEGER, PARAMETER :: LA3FE = 34 ! Iron aerosol (place holder)
1240 INTEGER, PARAMETER :: LB2MN = 35 ! Manganese aerosol (place holder)
1241 INTEGER, PARAMETER :: LK = 36 ! Potassium aerosol (Cl- tracked separately) (place holder)
1243 !
1244 ! on input
1245 !
1247 ! only local wrf dimensions are need as of now in this routine
1249 integer &
1250 ,intent (in ) :: &
1251 ids,ide, jds,jde, kds,kde,scalaropt, &
1252 numgas, conv_tr_wetscav,conv_tr_aqchem, &
1253 num_chem,ims,ime, jms,jme, kms,kme,chemopt, &
1254 its,ite, jts,jte, kts,kte,ipr,jpr,j,npr
1255 real, dimension (its:ite,kts:kte) &
1256 ,intent (in ) :: &
1257 z_cup,cd,zu,p,hstary,t
1258 real, dimension (its:ite,kts:kte) &
1259 ,intent (inout ) :: &
1260 cupclw,clw_all
1261 real, dimension (its:ite,kts:kte,1:num_chem) &
1262 ,intent (inout ) :: &
1263 tr_up,tr_c,tr_pw
1264 real, dimension (its:ite,kts:kte,1:num_chem) &
1265 ,intent (in ) :: &
1266 tre_cup,tracer
1267 real, dimension (its:ite) &
1268 ,intent (in ) :: &
1269 pre
1271 ! entr= entrainment rate
1272 real &
1273 ,intent (in ) :: &
1274 mentr_rate,tcrit
1275 integer, dimension (its:ite) &
1276 ,intent (in ) :: &
1277 kbcon,ktop,k22
1278 ! ierr error value, maybe modified in this routine
1280 integer, dimension (its:ite) &
1281 ,intent (inout) :: &
1282 ierr
1283 character *(*), intent (in) :: &
1284 name
1285 ! local variables in this routine
1286 !
1287 real :: conc_equi,conc_mxr,partialp,taucld
1289 integer :: &
1290 iall,i,k,iwd,nv
1291 real :: &
1292 trcc,trch,dh,c0,dz,radius,airm,dens
1293 integer :: &
1294 itf,ktf,iaer,igas
1296 ! Gas/aqueous phase partitioning for wet scavenging/deposition of gas
1297 ! phase and aerosol species:
1298 real aq_gas_ratio
1300 ! I/O for AQCHEM:
1302 real precip ! Precipitation rate (mm/h)
1303 real, dimension (ngas) :: gas,gaswdep
1304 real, dimension (naer) :: aerosol,aerwdep
1305 real, dimension (nliqs) :: liquid
1306 real hpwdep
1307 real alfa0,alfa2,alfa3 ! Aerosol scavenging coefficients for Aitken mode
1309 itf=MIN(ite,ide-1)
1310 ktf=MIN(kte,kde-1)
1312 !
1313 iall=0
1314 c0=.002
1315 iwd=0
1316 !
1317 !--- no precip for small clouds
1318 !
1319 ! if(mentr_rate.gt.0.)then
1320 ! radius=.2/mentr_rate
1321 ! if(radius.lt.900.)c0=0.
1322 if(name.eq.'shallow')c0=0.
1323 ! if(radius.lt.900.)iall=0
1324 ! endif
1325 do nv=2,num_chem
1326 do k=kts,ktf
1327 do i=its,itf
1328 ! Initialize species mass in rain water:
1329 tr_pw(i,k,nv)=0.
1330 ! Initialize species mass in updraft:
1331 if(ierr(i).eq.0)tr_up(i,k,nv)=tre_cup(i,k,nv)
1332 ! Initialize species mass in cloud + rain water:
1333 tr_c(i,k,nv)=0.
1334 enddo
1335 enddo
1336 enddo
1338 do nv=2,num_chem
1339 do i=its,itf
1340 if(ierr(i).eq.0.)then
1341 do k=k22(i),kbcon(i)-1
1342 tr_up(i,k,nv)=tre_cup(i,k22(i),nv)
1343 enddo
1344 endif
1345 enddo
1346 enddo
1348 DO 100 k=kts+1,ktf-1
1349 DO 100 i=its,itf
1351 IF(ierr(i).ne.0)GO TO 100
1352 IF(K.Lt.KBCON(I))GO TO 100
1353 IF(K.Gt.KTOP(I)+1)GO TO 100
1355 DZ=Z_cup(i,K)-Z_cup(i,K-1)
1357 if(cupclw(i,k).le.0.)cupclw(i,k)=0.
1358 if(clw_all(i,k).le.0.)clw_all(i,k)=0.
1360 !
1361 ! Steady state plume equation, for what could be in cloud before anything
1362 ! happens (kg/kg). tr_up would be the concentration if tracers were conserved.
1363 !
1365 do nv=2,num_chem
1366 if(i.eq.ipr.and.j.eq.jpr.and.nv.eq.npr)print *,k,tr_up(i,K-1,nv),tr_up(i,K,nv),tr_pw(i,k-1,nv),clw_all(i,k),cupclw(i,k)
1367 tr_up(i,K,nv)=(tr_up(i,K-1,nv)*(1.-.5*CD(i,K)*DZ)+mentr_rate* &
1368 DZ*tracer(i,K-1,nv))/(1.+mentr_rate*DZ-.5*cd(i,k)*dz)
1369 tr_up(i,k,nv)=max(1.e-16,tr_up(i,K,nv))
1370 enddo
1372 !
1373 ! Aqueous chemistry
1374 !
1376 if ((chemopt .EQ. RADM2SORG .OR. chemopt .EQ. RADM2SORG_AQ .OR. chemopt .EQ. RACMSORG_AQ .OR. chemopt .EQ. RACMSORG_KPP .OR. chemopt .EQ. RADM2SORG_KPP .OR. &
1377 chemopt .EQ. RACM_ESRLSORG_KPP .OR. chemopt .EQ. RACM_SOA_VBS_KPP .OR. chemopt .EQ. RADM2SORG_AQCHEM .OR. chemopt .EQ. RACMSORG_AQCHEM) &
1378 .and. (conv_tr_aqchem == 1)) then
1380 !
1381 ! For MADE/SORGAM derived schemes with aqueous chemistry
1382 !
1384 ! Air mass density
1385 dens = 0.1*p(i,k)/t(i,k)*mwdry/8.314472 ! kg/m3
1387 ! Column air number density:
1388 airm = 1000.0*dens*dz/mwdry ! mol/m2
1390 ! Wet scavenging initialization for AQCHEM
1392 GASWDEP = 0.0
1393 AERWDEP = 0.0
1394 HPWDEP = 0.0
1396 ! We provide a precipitation rate and aerosol scavenging rates of zero,
1397 ! in order to prevent wet scavenging in AQCHEM (it is treated later):
1399 precip = 0.0 ! mm/hr
1401 alfa0 = 0.0
1402 alfa2 = 0.0
1403 alfa3 = 0.0
1405 ! Gas phase concentrations before aqueous phase chemistry
1406 ! (with units conversion ppm -> mol/mol)
1408 gas(:) = 0.0
1410 gas(lco2) = 380.0e-6
1412 gas(lso2) = tr_up(i,k,p_so2)*1.0e-6
1413 gas(lhno3) = tr_up(i,k,p_hno3)*1.0e-6
1414 gas(ln2o5) = tr_up(i,k,p_n2o5)*1.0e-6
1415 gas(lnh3) = tr_up(i,k,p_nh3)*1.0e-6
1416 gas(lh2o2) = tr_up(i,k,p_h2o2)*1.0e-6
1417 gas(lo3) = tr_up(i,k,p_o3)*1.0e-6
1418 gas(lfoa) = tr_up(i,k,p_ora1)*1.0e-6
1419 gas(lmhp) = tr_up(i,k,p_op1)*1.0e-6
1420 gas(lpaa) = tr_up(i,k,p_paa)*1.0e-6
1421 gas(lh2so4) = tr_up(i,k,p_sulf)*1.0e-6
1423 ! Aerosol mass concentrations before aqueous phase chemistry
1424 ! (with units conversion ug/kg -> mol/mol). Although AQCHEM
1425 ! accounts for much of the aerosol compounds in MADE, they are
1426 ! not treated at the moment by AQCHEM, as the mapping between
1427 ! the organic compound groups in MADE and AQCHEM is not obvious.
1429 aerosol(:) = 0.0
1431 ! We assume all accumulation mode particles are activated in cumulus clouds:
1433 aerosol(lso4acc) = tr_up(i,k,p_so4aj)*1.0e-9*mwdry/mwso4
1434 aerosol(lnh4acc) = tr_up(i,k,p_nh4aj)*1.0e-9*mwdry/mwnh4
1435 aerosol(lno3acc) = tr_up(i,k,p_no3aj)*1.0e-9*mwdry/mwno3
1437 ! Cloud lifetime:
1438 taucld = 1800.0
1440 if (clw_all(i,k)*dens .gt. qcldwtr_cutoff) then ! Cloud water > threshold
1441 CALL AQCHEM( &
1442 t(i,k), &
1443 p(i,k)*100., &
1444 taucld, &
1445 precip, &
1446 clw_all(i,k)*dens, &
1447 clw_all(i,k)*dens, &
1448 airm, &
1449 ALFA0, &
1450 ALFA2, &
1451 ALFA3, &
1452 GAS, &
1453 AEROSOL, &
1454 LIQUID, &
1455 GASWDEP, &
1456 AERWDEP, &
1457 HPWDEP)
1458 endif
1460 ! Gas phase concentrations after aqueous phase chemistry
1461 ! (with units conversion mol/mol -> ppm)
1463 tr_up(i,k,p_so2) = gas(lso2)*1.0e6
1464 tr_up(i,k,p_hno3) = gas(lhno3)*1.0e6
1465 tr_up(i,k,p_n2o5) = gas(ln2o5)*1.0e6
1466 tr_up(i,k,p_nh3) = gas(lnh3)*1.0e6
1467 tr_up(i,k,p_h2o2) = gas(lh2o2)*1.0e6
1468 tr_up(i,k,p_o3) = gas(lo3)*1.0e6
1469 tr_up(i,k,p_ora1) = gas(lfoa)*1.0e6
1470 tr_up(i,k,p_op1) = gas(lmhp)*1.0e6
1471 tr_up(i,k,p_paa) = gas(lpaa)*1.0e6
1472 tr_up(i,k,p_sulf) = gas(lh2so4)*1.0e6
1474 ! Aerosol mass concentrations
1475 ! (with units conversion mol/mol -> ug/kg)
1477 tr_up(i,k,p_so4aj) = aerosol(lso4acc)*1.0e9*mwso4/mwdry
1478 tr_up(i,k,p_nh4aj) = aerosol(lnh4acc)*1.0e9*mwnh4/mwdry
1479 tr_up(i,k,p_no3aj) = aerosol(lno3acc)*1.0e9*mwno3/mwdry
1481 endif
1483 ! wet scavenging option (turn off by setting
1484 !
1485 if (conv_tr_wetscav == 1) then
1488 do nv=2,num_chem
1490 aq_gas_ratio = 0.0
1492 ! Fraction of gas phase species that partions into the liquid phase:
1494 if (nv.eq.p_so2) aq_gas_ratio = 1.0
1495 if (nv.eq.p_sulf) aq_gas_ratio = 1.0
1496 if (nv.eq.p_nh3) aq_gas_ratio = 1.0
1497 if (nv.eq.p_hno3) aq_gas_ratio = 1.0
1498 if (nv.gt.numgas) aq_gas_ratio = 0.5
1499 if (nv.eq.p_bc1 .or. nv.eq.p_oc1 .or. nv.eq.p_dms) aq_gas_ratio=0.
1500 if (nv.eq.p_sulf .or. nv.eq.p_seas_1 .or. nv.eq.p_seas_2) aq_gas_ratio=1.
1501 if (nv.eq.p_seas_3 .or. nv.eq.p_seas_4) aq_gas_ratio=1.
1502 if (nv.eq.p_bc2 .or. nv.eq.p_oc2) aq_gas_ratio=0.8
1504 ! Assume that all non-gas species (aerosol mass and number)
1505 ! partion completely into the liquid phase:
1508 if (aq_gas_ratio > 0.0) then
1509 tr_c(i,k,nv) = aq_gas_ratio*tr_up(i,k,nv) ! Amount of species cloud + rain water
1510 trch = tr_up(i,k,nv)-tr_c(i,k,nv) ! Amount of species remaining in gas phase
1511 trcc = (tr_up(i,k,nv)-trch)/(1.+c0*dz*zu(i,k)) ! Amount of species cloud + rain water
1512 tr_pw(i,k,nv) = c0*dz*trcc*zu(i,k) ! Amount of species in rain water
1513 tr_up(i,k,nv) = trcc+trch ! Total amount of species in updraft = conc in liq water (trcc) + conc in gas phase (trch)
1514 endif
1516 enddo
1518 endif ! parameterization wetscavenging option
1521 100 CONTINUE
1523 END subroutine cup_up_tracer
1526 SUBROUTINE cup_dd_tracer(ierr,z_cup,qrcd,tracer,tre_cup,tr_up,tr_dd, &
1527 tr_pw,tr_pwd,jmin,cdd,entr,zd,pwdper,wetdep,xmb,k22, &
1528 numgas,ids,ide, jds,jde, kds,kde, &
1529 ims,ime, jms,jme, kms,kme, &
1530 its,ite, jts,jte, kts,kte)
1531 USE module_configure
1532 USE module_state_description
1533 USE module_model_constants, only: mwdry
1535 implicit none
1536 !
1537 ! on input
1538 !
1540 ! only local wrf dimensions are need as of now in this routine
1542 integer &
1543 ,intent (in ) :: &
1544 numgas,ids,ide, jds,jde, kds,kde, &
1545 ims,ime, jms,jme, kms,kme, &
1546 its,ite, jts,jte, kts,kte
1547 real, dimension (its:ite,kts:kte) &
1548 ,intent (in ) :: &
1549 pwdper,zd,cdd,qrcd,z_cup
1550 real, dimension (its:ite) &
1551 ,intent (in ) :: &
1552 xmb
1553 real, dimension (its:ite,kts:kte,1:num_chem) &
1554 ,intent (inout ) :: &
1555 tr_dd,tr_pwd,tr_up
1556 real, dimension (its:ite,kts:kte,1:num_chem) &
1557 ,intent (in ) :: &
1558 tre_cup,tracer,tr_pw
1559 real, dimension (its:ite,1:num_chem) &
1560 ,intent (out ) :: &
1561 wetdep
1563 ! entr= entrainment rate
1564 real &
1565 ,intent (in ) :: &
1566 entr
1567 integer, dimension (its:ite) &
1568 ,intent (in ) :: &
1569 jmin
1570 ! ierr error value, maybe modified in this routine
1572 integer, dimension (its:ite) &
1573 ,intent (inout) :: &
1574 ierr,k22
1575 ! local variables in this routine
1576 !
1578 integer :: &
1579 iall,i,k,nv,ki,kj
1580 real :: &
1581 dh,c0,dz,radius
1582 integer :: &
1583 itf,ktf
1585 real :: &
1586 evaporate,condensate
1588 itf=MIN(ite,ide-1)
1589 ktf=MIN(kte,kde-1)
1591 ! Initialize the tracer amount that evaporated from rain water:
1592 tr_pwd(its:ite,kts:kte,1:num_chem) = 0.0
1594 ! Initialize wet deposition:
1595 wetdep(its:ite,1:num_chem) = 0.0
1597 ! Calculate wet deposition with re-evaporation (based on wet scavenging in cup_up_tracer);
1598 ! assume no gas takeup by rain during fall for now
1600 do i=its,ite
1602 IF(ierr(I).eq.0)then
1604 do nv=1,numgas ! Only gas phase species evaporate along with rain water
1606 do k=ktf,kts,-1 ! Descending loop over all levels
1608 ! Start with initializing the condensate with the tracer amount in rain water on this level:
1609 condensate = tr_pw(i,k,nv)
1611 do kj=k,kts,-1 ! Descending loop over this and lower levels
1612 ! Evaporated tracer amount at the current level:
1613 evaporate = condensate*pwdper(i,kj)
1614 ! Accumulate the evaporated tracer amount:
1615 tr_pwd(i,kj,nv) = tr_pwd(i,kj,nv) + evaporate
1616 ! Remove the evaporated tracer amount from the condensate:
1617 condensate = max(0.0,condensate - evaporate)
1618 enddo
1620 enddo
1622 ! Calculate the wet deposition as the column integral over the
1623 ! tracer amount in rain water - evaporated tracer amount:
1625 do k=kts,ktf
1626 wetdep(i,nv) = wetdep(i,nv) + tr_pw(i,k,nv) - tr_pwd(i,k,nv)
1627 enddo
1629 wetdep(i,nv) = max(0.0,wetdep(i,nv))
1631 enddo
1633 endif
1635 enddo
1637 ! In downdraft, do only transport of tracers
1639 do i=its,ite
1641 IF(ierr(I).eq.0)then
1643 do nv=1,num_chem
1644 tr_dd(i,jmin(i),nv)=tre_cup(i,jmin(i),nv) ! Tracer amount in downdraft
1645 enddo
1647 do ki=jmin(i)-1,1,-1
1648 DZ=Z_cup(i,ki+1)-Z_cup(i,ki)
1649 do nv=1,num_chem
1650 tr_dd(i,ki,nv)=(tr_dd(i,ki+1,nv)*(1.-.5*CDD(i,ki)*DZ) &
1651 +entr*DZ*tracer(i,ki,nv) &
1652 )/(1.+entr*DZ-.5*CDD(i,ki)*DZ)
1653 enddo
1655 enddo
1657 endif
1659 enddo
1661 ! Add evaporation from rain water:
1663 do i=its,ite
1664 IF(ierr(I).eq.0)then
1665 do nv=1,num_chem
1666 do k=kts,ktf
1667 tr_dd(i,k,nv) = tr_dd(i,k,nv) + tr_pwd(i,k,nv)
1668 enddo
1669 enddo
1670 endif
1671 enddo
1673 ! To obtain wet deposition, multiply with base mass flux; units are given
1674 ! assuming that the base mass flux units [xmb] = kg(air)/m2/s:
1676 do nv=1,num_chem
1677 if (nv <= numgas) then
1678 do i=its,itf
1679 wetdep(i,nv)=wetdep(i,nv)*xmb(i)/mwdry ! mmol/m2/s for gas phase species
1680 enddo
1681 else
1682 do i=its,itf
1683 wetdep(i,nv)=wetdep(i,nv)*xmb(i) ! ug/m2/s for aerosol mass, #/m2/s for aerosol number
1684 enddo
1685 endif
1686 enddo
1688 END subroutine cup_dd_tracer
1691 SUBROUTINE neg_check_ct(name,pret,ktop,epsilc,dt,q,outq,iopt,num_chem, &
1692 its,ite,kts,kte,itf,ktf,ipr,jpr,npr,j)
1694 INTEGER, INTENT(IN ) :: iopt,num_chem,its,ite,kts,kte,itf,ktf,ipr,jpr,npr,j
1696 real, dimension (its:ite,kts:kte,num_chem ) , &
1697 intent(inout ) :: &
1698 q,outq
1699 real, dimension (its:ite ) , &
1700 intent(in ) :: &
1701 pret
1702 integer, dimension (its:ite ) , &
1703 intent(in ) :: &
1704 ktop
1705 real &
1706 ,intent (in ) :: &
1707 dt,epsilc
1708 real :: tracermin,tracermax,thresh,qmem,qmemf,qmem2,qtest,qmem1
1709 character *(*) name
1710 !
1711 ! check whether routine produces negative q's. This can happen, since
1712 ! tendencies are calculated based on forced q's. This should have no
1713 ! influence on conservation properties, it scales linear through all
1714 ! tendencies. Use iopt=0 to test for each tracer seperately, iopt=1
1715 ! for a more severe limitation...
1716 !
1717 thresh=epsilc
1718 ! thresh=1.e-30
1719 if(iopt.eq.0)then
1720 do nv=2,num_chem
1721 do 100 i=its,itf
1722 if(pret(i).le.0.and.name.eq.'deep')go to 100
1723 tracermin=q(i,kts,nv)
1724 tracermax=q(i,kts,nv)
1725 do k=kts+1,kte-1
1726 tracermin=min(tracermin,q(i,k,nv))
1727 tracermax=max(tracermax,q(i,k,nv))
1728 enddo
1729 tracermin=max(tracermin,thresh)
1730 qmemf=1.
1731 !
1732 ! first check for minimum restriction
1733 !
1734 do k=kts,ktop(i)
1735 !
1736 ! tracer tendency
1737 !
1738 qmem=outq(i,k,nv)
1739 !
1740 ! only necessary if there is a tendency
1741 !
1742 if(qmem.lt.0.)then
1743 qtest=q(i,k,nv)+outq(i,k,nv)*dt
1744 if(qtest.lt.tracermin)then
1745 !
1746 ! qmem2 would be the maximum allowable tendency
1747 !
1748 qmem1=outq(i,k,nv)
1749 qmem2=(tracermin-q(i,k,nv))/dt
1750 qmemf=min(qmemf,qmem2/qmem1)
1751 if(qmemf.gt.1.)print *,'something wrong in negct_1',qmem2,qmem1
1752 if(i.eq.ipr.and.j.eq.jpr.and.nv.eq.npr)then
1753 print *,k,qtest,qmem2,qmem1,qmemf
1754 endif
1755 qmemf=max(qmemf,0.)
1756 endif
1757 endif
1758 enddo
1759 do k=kts,ktop(i)
1760 outq(i,k,nv)=outq(i,k,nv)*qmemf
1761 enddo
1762 !
1763 ! now check max
1764 !
1765 qmemf=1.
1766 do k=kts,ktop(i)
1767 !
1768 ! tracer tendency
1769 !
1770 qmem=outq(i,k,nv)
1771 !
1772 ! only necessary if there is a tendency
1773 !
1774 if(qmem.gt.0.)then
1775 qtest=q(i,k,nv)+outq(i,k,nv)*dt
1776 if(qtest.gt.tracermax)then
1777 !
1778 ! qmem2 would be the maximum allowable tendency
1779 !
1780 qmem1=outq(i,k,nv)
1781 qmem2=(tracermax-q(i,k,nv))/dt
1782 qmemf=min(qmemf,qmem2/qmem1)
1783 if(qmemf.gt.1.)print *,'something wrong in negct_2',qmem2,qmem1
1784 if(i.eq.ipr.and.j.eq.jpr.and.nv.eq.npr)then
1785 print *,'2',k,qtest,qmem2,qmem1,qmemf
1786 endif
1787 qmemf=max(qmemf,0.)
1788 endif
1789 endif
1790 enddo
1791 do k=kts,ktop(i)
1792 outq(i,k,nv)=outq(i,k,nv)*qmemf
1793 enddo
1794 100 continue
1795 enddo
1796 !
1797 ! ELSE
1798 !
1799 elseif(iopt.eq.1)then
1800 do i=its,itf
1801 qmemf=1.
1802 do k=kts,ktop(i)
1803 do nv=2,num_chem
1804 !
1805 ! tracer tendency
1806 !
1807 qmem=outq(i,k,nv)
1808 !
1809 ! only necessary if tendency is larger than zero
1810 !
1811 if(qmem.lt.0.)then
1812 qtest=q(i,k,nv)+outq(i,k,nv)*dt
1813 if(qtest.lt.thresh)then
1814 !
1815 ! qmem2 would be the maximum allowable tendency
1816 !
1817 qmem1=outq(i,k,nv)
1818 qmem2=(thresh-q(i,k,nv))/dt
1819 qmemf=min(qmemf,qmem2/qmem1)
1820 qmemf=max(0.,qmemf)
1821 endif
1822 endif
1823 enddo
1824 enddo
1825 do nv=2,num_chem
1826 do k=kts,ktop(i)
1827 outq(i,k,nv)=outq(i,k,nv)*qmemf
1828 enddo
1829 enddo
1830 enddo
1831 endif
1833 END SUBROUTINE neg_check_ct
1835 !-------------------------------------------------------
1836 END MODULE module_ctrans_grell