| blob | 910f9f4a62f29a0004a03a4461db1d099306e5b0 |
1 !WRF:MODEL_LAYER:PHYSICS
2 !
4 MODULE module_ctrans_grell
5 !USE module_data_radm2
6 USE module_cu_gd
7 USE module_dep_simple
8 !USE module_ctrans_aqchem
9 ! USE module_configure
10 ! USE module_state_description
11 ! Mole weight
12 ! REAL :: WTM(NUMCHEM_radm)
13 ! DATA WTM / 64.,96.,46.,30.,48.,63.,34.,44.,30.,48.,62., &
14 ! 76.,46.,60.,17.,108.,62.,121.,44.,72.,114.,30., &
15 ! 28.,28.,42.,56.,92.,106.,75.,147.,47.,79.,72., &
16 ! 58.,72.,87.,119.,108.,68.,17.,33. /
18 CONTAINS
20 !-------------------------------------------------------------
21 SUBROUTINE GRELLDRVCT(DT,itimestep,DX, &
22 id,config_flags,rho_phy,RAINCV,chem, &
23 U,V,t_phy,moist,dz8w,p_phy, &
24 XLV,CP,G,r_v,z,cu_co_ten, &
25 numgas, &
26 ids,ide, jds,jde, kds,kde, &
27 ims,ime, jms,jme, kms,kme, &
28 its,ite, jts,jte, kts,kte )
29 USE module_configure
30 USE module_state_description
31 !-------------------------------------------------------------
32 IMPLICIT NONE
33 !-------------------------------------------------------------
34 INTEGER, INTENT(IN ) :: &
35 id, numgas, &
36 ids,ide, jds,jde, kds,kde, &
37 ims,ime, jms,jme, kms,kme, &
38 its,ite, jts,jte, kts,kte
41 INTEGER, INTENT(IN ) :: ITIMESTEP
43 REAL, INTENT(IN ) :: XLV, R_v
44 REAL, INTENT(IN ) :: CP,G
46 REAL, DIMENSION( ims:ime , kms:kme , jms:jme,num_moist ) , &
47 INTENT(IN ) :: moist
48 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
49 INTENT(IN ) :: &
50 U, &
51 V, &
52 t_phy, &
53 z, &
54 p_phy, &
55 dz8w, &
56 rho_phy
57 !
58 ! on output for control only, purely diagnostic
59 !
60 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
61 INTENT(INOUT ) :: &
62 cu_co_ten
65 !
66 REAL, INTENT(IN ) :: DT, DX
67 !
68 REAL, DIMENSION( ims:ime , kms:kme , jms:jme, num_chem ), &
69 INTENT(INOUT) :: &
70 chem
73 REAL, DIMENSION( ims:ime , jms:jme ), &
74 INTENT(IN) :: RAINCV
76 ! LOCAL VARS
77 real, dimension (its:ite,kts:kte) :: &
78 OUTT,OUTQ,OUTQC
79 real, dimension (its:ite) :: &
80 pret, ter11
82 !
83 ! basic environmental input includes moisture convergence (mconv)
84 ! omega (omeg), windspeed (us,vs), and a flag (aaeq) to turn off
85 ! convection for this call only and at that particular gridpoint
86 !
87 real, dimension (its:ite,kts:kte) :: &
88 T,TN,q,qo,PO,P,US,VS,hstary
89 real, dimension (its:ite,kts:kte,num_chem) :: &
90 tracer,tracert
91 real, dimension (its:ite) :: &
92 Z1,PSUR,AAEQ
93 integer, dimension (its:ite) :: &
94 ktop
95 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
97 INTEGER :: nv,i,j,k,ICLDCK,ipr,jpr,npr
98 REAL :: tcrit,dp,dq,epsilc
99 INTEGER :: itf,jtf,ktf,iopt
100 epsilc=1.e-30
101 ! return
102 ! ipr=111
103 ! jpr=40
104 ! if(itimestep.lt.34.or.itimestep.gt.36)ipr=0
105 ! if(itimestep.lt.34.or.itimestep.gt.36)jpr=0
106 ! ipr=61
107 ! jpr=60
108 ipr=0
109 jpr=0
110 npr=p_co
111 tcrit=258.
112 iopt=0
113 itf=MIN(ite,ide-1)
114 ktf=MIN(kte,kde-1)
115 jtf=MIN(jte,jde-1)
116 !
117 !
118 ! DO J = jts,jtf
119 ! DO I=ITS,ITF
120 ! if(raincv(i,j).gt.0.)then
121 ! ipr=i
122 ! jpr=j
123 ! go to 123
124 ! endif
125 ! ENDDO
126 ! ENDDO
127 123 continue
128 ! print *,ipr,jpr
129 DO 100 J = jts,jtf
130 if(j.eq.jpr)print *,'dt = ',dt
131 DO I=ITS,ITF
132 ktop(i)=0
133 PSUR(I)=p_phy(I,kts,J)*.01
134 TER11(I)=z(i,kts,j)
135 aaeq(i)=0.
136 !
137 ! rainrate is input for this transport/wet-deposition routine
138 !
139 pret(i)=raincv(i,j)/dt
140 if(pret(i).le.0.)aaeq(i)=20.
141 ENDDO
142 DO K=kts,ktf
143 DO I=ITS,ITF
144 po(i,k)=p_phy(i,k,j)*.01
145 P(I,K)=PO(i,k)
146 US(I,K) =u(i,k,j)
147 VS(I,K) =v(i,k,j)
148 T(I,K)=t_phy(i,k,j)
149 q(I,K)=moist(i,k,j,p_qv)
150 IF(Q(I,K).LT.1.E-08)Q(I,K)=1.E-08
151 ENDDO
152 ENDDO
153 do nv=1,num_chem
154 DO K=kts,ktf
155 DO I=ITS,ITF
156 tracer(i,k,nv)=max(epsilc,chem(i,k,j,nv))
157 tracert(i,k,nv)=0.
158 ENDDO
159 ENDDO
160 ENDDO
161 DO K=kts,ktf
162 DO I=ITS,ITF
163 cu_co_ten(i,k,j)=0.
164 ! hstary(i,k)=hstar4(nv)*exp(dhr(nv)*(1./t(i,k)-1./298.))
165 if(i.eq.ipr.and.j.eq.jpr)then
166 print *,k,pret(i),tracer(i,k,npr),p(i,k),z(i,k,j)
167 endif
168 ENDDO
169 ENDDO
170 ! ENDDO
171 !
172 !---- CALL CUMULUS PARAMETERIZATION
173 !
174 CALL CUP_ct(ktop,tracer,j,AAEQ,T,Q,TER11,PRET,P,tracert, &
175 hstary,DT,PSUR,US,VS,tcrit, &
176 xlv,r_v,cp,g,ipr,jpr,npr,num_chem, &
177 ids,ide, jds,jde, kds,kde, &
178 ims,ime, jms,jme, kms,kme, &
179 its,ite, jts,jte, kts,kte )
181 do nv=1,num_chem
182 DO I=its,itf
183 if(pret(i).le.0.)then
184 DO K=kts,ktf
185 tracert(i,k,nv)=0.
186 ENDDO
187 endif
188 enddo
189 enddo
190 CALL neg_check_ct(pret,ktop,epsilc,dt,tracer,tracert,iopt,num_chem, &
191 its,ite,kts,kte,itf,ktf,ipr,jpr,npr,j)
192 do nv=1,num_chem
193 DO I=its,itf
194 if(pret(i).gt.0.)then
195 DO K=kts,ktf
196 chem(i,k,j,nv)=max(epsilc,chem(i,k,j,nv)+tracert(i,k,nv)*dt)
197 if(nv.eq.npr)then
198 cu_co_ten(i,k,j)=tracert(i,k,npr)*dt
199 if(i.eq.ipr.and.j.eq.jpr)print *,k,chem(i,k,j,nv),cu_co_ten(i,k,j)
200 endif
201 ENDDO
202 else
203 DO K=kts,ktf
204 tracert(i,k,nv)=0.
205 if(nv.eq.npr)cu_co_ten(i,k,j)=0.
206 enddo
207 endif
208 ENDDO
209 ENDDO
212 100 continue
214 END SUBROUTINE GRELLDRVCT
217 SUBROUTINE CUP_ct(ktop,tracer,J,AAEQ,T,Q,Z1, &
218 PRE,P,tracert,hstary,DTIME,PSUR,US,VS,TCRIT, &
219 xl,rv,cp,g,ipr,jpr,npr,num_chem, &
220 ids,ide, jds,jde, kds,kde, &
221 ims,ime, jms,jme, kms,kme, &
222 its,ite, jts,jte, kts,kte )
224 IMPLICIT NONE
226 integer &
227 ,intent (in ) :: &
228 num_chem,ids,ide, jds,jde, kds,kde, &
229 ims,ime, jms,jme, kms,kme, &
230 its,ite, jts,jte, kts,kte,ipr,jpr,npr
231 integer, intent (in ) :: &
232 j
233 !
234 !
235 !
236 !tracert = output temp tendency (per s)
237 ! pre = input precip
238 real, dimension (its:ite,kts:kte,num_chem) &
239 ,intent (inout ) :: &
240 tracert,tracer
241 real, dimension (its:ite) &
242 ,intent (inout ) :: &
243 pre
244 integer, dimension (its:ite) &
245 ,intent (inout ) :: &
246 ktop
247 integer, dimension (its:ite) :: &
248 kbcon
249 !
250 ! basic environmental input includes moisture convergence (mconv)
251 ! omega (omeg), windspeed (us,vs), and a flag (aaeq) to turn off
252 ! convection for this call only and at that particular gridpoint
253 !
254 real, dimension (its:ite,kts:kte) &
255 ,intent (in ) :: &
256 T,P,US,VS,HSTARY
257 real, dimension (its:ite,kts:kte) &
258 ,intent (inout) :: &
259 Q
260 real, dimension (its:ite) &
261 ,intent (in ) :: &
262 Z1,PSUR,AAEQ
265 real &
266 ,intent (in ) :: &
267 dtime,tcrit,xl,cp,rv,g
270 real, dimension (its:ite,1:3) :: &
271 edtc
272 !
273 !
274 !
275 !***************** the following are your basic environmental
276 ! variables. They carry a "_cup" if they are
277 ! on model cloud levels (staggered). They carry
278 ! an "o"-ending (z becomes zo), if they are the forced
279 ! variables. They are preceded by x (z becomes xz)
280 ! to indicate modification by some typ of cloud
281 !
282 ! z = heights of model levels
283 ! q = environmental mixing ratio
284 ! qes = environmental saturation mixing ratio
285 ! t = environmental temp
286 ! p = environmental pressure
287 ! he = environmental moist static energy
288 ! hes = environmental saturation moist static energy
289 ! z_cup = heights of model cloud levels
290 ! q_cup = environmental q on model cloud levels
291 ! qes_cup = saturation q on model cloud levels
292 ! t_cup = temperature (Kelvin) on model cloud levels
293 ! p_cup = environmental pressure
294 ! he_cup = moist static energy on model cloud levels
295 ! hes_cup = saturation moist static energy on model cloud levels
296 ! gamma_cup = gamma on model cloud levels
297 !
298 !
299 ! hcd = moist static energy in downdraft
300 ! zd normalized downdraft mass flux
301 ! dby = buoancy term
302 ! entr = entrainment rate
303 ! zd = downdraft normalized mass flux
304 ! entr= entrainment rate
305 ! hcd = h in model cloud
306 ! bu = buoancy term
307 ! zd = normalized downdraft mass flux
308 ! gamma_cup = gamma on model cloud levels
309 ! mentr_rate = entrainment rate
310 ! qcd = cloud q (including liquid water) after entrainment
311 ! qrch = saturation q in cloud
312 ! pwd = evaporate at that level
313 ! pwev = total normalized integrated evaoprate (I2)
314 ! entr= entrainment rate
315 ! z1 = terrain elevation
316 ! entr = downdraft entrainment rate
317 ! jmin = downdraft originating level
318 ! kdet = level above ground where downdraft start detraining
319 ! psur = surface pressure
320 ! z1 = terrain elevation
321 ! zd = downdraft normalized mass flux
322 ! zu = updraft normalized mass flux
323 ! mbdt = arbitrary numerical parameter
324 ! dtime = dt over which forcing is applied
325 ! kbcon = LFC of parcel from k22
326 ! k22 = updraft originating level
327 ! dby = buoancy term
328 ! ktop = cloud top (output)
329 ! xmb = total base mass flux
330 ! hc = cloud moist static energy
331 ! hkb = moist static energy at originating level
332 ! mentr_rate = entrainment rate
334 real, dimension (its:ite,kts:kte) :: &
335 he,hes,qes,z,pwdper, &
337 qes_cup,q_cup,he_cup,hes_cup,z_cup,p_cup,gamma_cup,t_cup, &
339 dby,qc,qrcd,pwd,pw,hcd,qcd,dbyd,hc,qrc,zu,zd,clw_all, &
341 ! cd = detrainment function for updraft
342 ! cdd = detrainment function for downdraft
344 cd,cdd,scr1,DELLAH,DELLAQ,DELLAT,DELLAQC
346 ! edt = epsilon
347 ! edt = epsilon
348 real, dimension (its:ite) :: &
349 edt,HKB,QKB, &
350 XMB,PWAV,PWEV,BU,cap_max,cap_max_increment
351 real, dimension (its:ite,kts:kte,num_chem) :: &
352 tr_c,tr_up,tr_dd,tre_cup,tr_pw,tr_pwd
353 real, dimension (its:ite,num_chem) :: &
354 trkb
355 integer, dimension (its:ite) :: &
356 kzdown,KDET,K22,KB,JMIN,kstabi,kstabm, & !-lxz
357 ierr,KBMAX
359 integer :: &
360 ki,I,K,KK
361 real :: &
362 day,dz,mbdt,entr_rate,radius,entrd_rate,mentr_rate,mentrd_rate, &
363 zcutdown,edtmax,edtmin,depth_min,zkbmax,z_detr,zktop, &
364 dh,cap_maxs
366 integer :: itf,jtf,ktf
368 itf=MIN(ite,ide-1)
369 ktf=MIN(kte,kde-1)
370 jtf=MIN(jte,jde-1)
372 !sms$distribute end
373 day=86400.
374 !
375 !--- specify entrainmentrate and detrainmentrate
376 !
377 radius=12000.
378 !
379 !--- gross entrainment rate (these may be changed later on in the
380 !--- program, depending what your detrainment is!!)
381 !
382 entr_rate=.2/radius
383 !
384 !--- entrainment of mass
385 !
386 mentrd_rate=0.
387 mentr_rate=entr_rate
388 !
389 !--- initial detrainmentrates
390 !
391 do k=kts,ktf
392 do i=its,itf
393 cd(i,k)=0.1*entr_rate
394 cdd(i,k)=0.
395 clw_all(i,k)=0.
396 enddo
397 enddo
398 !
399 !--- max/min allowed value for epsilon (ratio downdraft base mass flux/updraft
400 ! base mass flux
401 !
402 edtmax=.8
403 edtmin=.2
404 !
405 !--- minimum depth (m), clouds must have
406 !
407 depth_min=500.
408 !
409 !--- maximum depth (mb) of capping
410 !--- inversion (larger cap = no convection)
411 !
412 cap_maxs=175.
413 !sms$to_local(grid_dh: <1, mix :size>, <2, mjx :size>) begin
414 DO 7 i=its,itf
415 kbmax(i)=1
416 cap_max_increment(i)=0.
417 edt(i)=0.
418 kstabm(i)=ktf-1
419 IERR(i)=0
420 if(aaeq(i).ne.0.)then
421 ierr(i)=20
422 endif
423 7 CONTINUE
424 do i=its,itf
425 cap_max(i)=cap_maxs
426 enddo
427 !
428 !--- max height(m) above ground where updraft air can originate
429 !
430 zkbmax=4000.
431 !
432 !--- height(m) above which no downdrafts are allowed to originate
433 !
434 zcutdown=3000.
435 !
436 !--- depth(m) over which downdraft detrains all its mass
437 !
438 z_detr=1250.
439 !
440 mbdt=dtime*4.E-03
441 !
442 !--- calculate moist static energy, heights, qes
443 !
444 call cup_env(z,qes,he,hes,t,q,p,z1, &
445 psur,ierr,tcrit,0,xl,cp, &
446 ids,ide, jds,jde, kds,kde, &
447 ims,ime, jms,jme, kms,kme, &
448 its,ite, jts,jte, kts,kte)
449 !
450 !--- environmental values on cloud levels
451 !
452 call cup_env_clev(t,qes,q,he,hes,z,p,qes_cup,q_cup,he_cup, &
453 hes_cup,z_cup,p_cup,gamma_cup,t_cup,psur, &
454 ierr,z1,xl,rv,cp, &
455 ids,ide, jds,jde, kds,kde, &
456 ims,ime, jms,jme, kms,kme, &
457 its,ite, jts,jte, kts,kte)
458 call cup_env_clev_tr(tracer,tre_cup,num_chem,ierr, &
459 ids,ide, jds,jde, kds,kde, &
460 ims,ime, jms,jme, kms,kme, &
461 its,ite, jts,jte, kts,kte)
462 do i=its,itf
463 if(ierr(i).eq.0)then
464 !
465 do k=kts,ktf-2
466 if(z_cup(i,k).gt.zkbmax+z1(i))then
467 kbmax(i)=k
468 go to 25
469 endif
470 enddo
471 25 continue
472 !
473 !
474 !--- level where detrainment for downdraft starts
475 !
476 do k=kts,ktf
477 if(z_cup(i,k).gt.z_detr+z1(i))then
478 kdet(i)=k
479 go to 26
480 endif
481 enddo
482 26 continue
483 !
484 endif
485 enddo
486 !
487 !
488 !
489 !------- DETERMINE LEVEL WITH HIGHEST MOIST STATIC ENERGY CONTENT - K22
490 !
491 CALL cup_MAXIMI(HE_CUP,3,KBMAX,K22,ierr, &
492 ids,ide, jds,jde, kds,kde, &
493 ims,ime, jms,jme, kms,kme, &
494 its,ite, jts,jte, kts,kte)
495 DO 36 i=its,itf
496 IF(ierr(I).eq.0.)THEN
497 IF(K22(I).GE.KBMAX(i))ierr(i)=2
498 endif
499 36 CONTINUE
500 !
501 !--- DETERMINE THE LEVEL OF CONVECTIVE CLOUD BASE - KBCON
502 !
503 call cup_kbcon(cap_max_increment,1,k22,kbcon,he_cup,hes_cup, &
504 ierr,kbmax,p_cup,cap_max, &
505 ids,ide, jds,jde, kds,kde, &
506 ims,ime, jms,jme, kms,kme, &
507 its,ite, jts,jte, kts,kte)
508 !
509 !--- increase detrainment in stable layers
510 !
511 CALL cup_minimi(HEs_cup,Kbcon,kstabm,kstabi,ierr, &
512 ids,ide, jds,jde, kds,kde, &
513 ims,ime, jms,jme, kms,kme, &
514 its,ite, jts,jte, kts,kte)
515 do i=its,itf
516 IF(ierr(I).eq.0.)THEN
517 if(kstabm(i)-1.gt.kstabi(i))then
518 do k=kstabi(i),kstabm(i)-1
519 cd(i,k)=cd(i,k-1)+1.5*entr_rate
520 if(cd(i,k).gt.10.0*entr_rate)cd(i,k)=10.0*entr_rate
521 enddo
522 ENDIF
523 ENDIF
524 ENDDO
525 !
526 !--- calculate incloud moist static energy
527 !
528 call cup_up_he(k22,hkb,z_cup,cd,mentr_rate,he_cup,hc, &
529 kbcon,ierr,dby,he,hes_cup, &
530 ids,ide, jds,jde, kds,kde, &
531 ims,ime, jms,jme, kms,kme, &
532 its,ite, jts,jte, kts,kte)
534 !--- DETERMINE CLOUD TOP - KTOP
535 !
536 call cup_ktop(1,dby,kbcon,ktop,ierr, &
537 ids,ide, jds,jde, kds,kde, &
538 ims,ime, jms,jme, kms,kme, &
539 its,ite, jts,jte, kts,kte)
540 DO 37 i=its,itf
541 kzdown(i)=0
542 if(ierr(i).eq.0)then
543 zktop=(z_cup(i,ktop(i))-z1(i))*.6
544 zktop=min(zktop+z1(i),zcutdown+z1(i))
545 do k=kts,ktf
546 if(z_cup(i,k).gt.zktop)then
547 kzdown(i)=k
548 go to 37
549 endif
550 enddo
551 endif
552 37 CONTINUE
553 !
554 !--- DOWNDRAFT ORIGINATING LEVEL - JMIN
555 !
556 call cup_minimi(HEs_cup,K22,kzdown,JMIN,ierr, &
557 ids,ide, jds,jde, kds,kde, &
558 ims,ime, jms,jme, kms,kme, &
559 its,ite, jts,jte, kts,kte)
560 DO 100 i=its,ite
561 IF(ierr(I).eq.0.)THEN
562 !
563 !--- check whether it would have buoyancy, if there where
564 !--- no entrainment/detrainment
565 !
566 101 continue
567 if(jmin(i)-1.lt.KDET(I))kdet(i)=jmin(i)-1
568 if(jmin(i).ge.Ktop(I)-1)jmin(i)=ktop(i)-2
569 ki=jmin(i)
570 hcd(i,ki)=hes_cup(i,ki)
571 DZ=Z_cup(i,Ki+1)-Z_cup(i,Ki)
572 dh=dz*(HCD(i,Ki)-hes_cup(i,ki))
573 dh=0.
574 !
575 do k=ki-1,1,-1
576 hcd(i,k)=hes_cup(i,jmin(i))
577 DZ=Z_cup(i,K+1)-Z_cup(i,K)
578 dh=dh+dz*(HCD(i,K)-hes_cup(i,k))
579 if(dh.gt.0.)then
580 jmin(i)=jmin(i)-1
581 if(jmin(i).gt.3)then
582 go to 101
583 else if(jmin(i).le.3)then
584 ierr(i)=9
585 go to 100
586 endif
587 endif
588 enddo
590 IF(JMIN(I).LE.3)then
591 ierr(i)=4
592 endif
594 ENDIF
595 100 continue
596 !
597 ! - Must have at least depth_min m between cloud convective base
598 ! and cloud top.
599 !
600 do i=its,itf
601 IF(ierr(I).eq.0.)THEN
602 IF(-z_cup(I,KBCON(I))+z_cup(I,KTOP(I)).LT.depth_min)then
603 ierr(i)=6
604 endif
605 endif
606 enddo
608 !
609 !c--- normalized updraft mass flux profile
610 !
611 call cup_up_nms(zu,z_cup,mentr_rate,cd,kbcon,ktop,ierr,k22, &
612 ids,ide, jds,jde, kds,kde, &
613 ims,ime, jms,jme, kms,kme, &
614 its,ite, jts,jte, kts,kte)
615 !
616 !c--- normalized downdraft mass flux profile,also work on bottom detrainment
617 !--- in this routine
618 !
619 call cup_dd_nms(zd,z_cup,cdd,mentrd_rate,jmin,ierr, &
620 0,kdet,z1, &
621 ids,ide, jds,jde, kds,kde, &
622 ims,ime, jms,jme, kms,kme, &
623 its,ite, jts,jte, kts,kte)
624 !
625 !--- downdraft moist static energy
626 !
627 call cup_dd_he(hes_cup,zd,hcd,z_cup,cdd,mentrd_rate, &
628 jmin,ierr,he,dbyd,he_cup, &
629 ids,ide, jds,jde, kds,kde, &
630 ims,ime, jms,jme, kms,kme, &
631 its,ite, jts,jte, kts,kte)
632 !
633 !--- calculate moisture properties of downdraft
634 !
635 call cup_dd_moisture(zd,hcd,hes_cup,qcd,qes_cup, &
636 pwd,q_cup,z_cup,cdd,mentrd_rate,jmin,ierr,gamma_cup, &
637 pwev,bu,qrcd,q,he,t_cup,2,xl, &
638 ids,ide, jds,jde, kds,kde, &
639 ims,ime, jms,jme, kms,kme, &
640 its,ite, jts,jte, kts,kte)
641 !
642 !--- calculate moisture properties of updraft
643 !
644 call cup_up_moisture(ierr,z_cup,qc,qrc,pw,pwav, &
645 kbcon,ktop,cd,dby,mentr_rate,clw_all, &
646 q,GAMMA_cup,zu,qes_cup,k22,q_cup,xl, &
647 ids,ide, jds,jde, kds,kde, &
648 ims,ime, jms,jme, kms,kme, &
649 its,ite, jts,jte, kts,kte)
650 !
651 !--- DETERMINE DOWNDRAFT STRENGTH IN TERMS OF WINDSHEAR
652 !
653 call cup_dd_edt(ierr,us,vs,z,ktop,kbcon,edt,p,pwav, &
654 pwev,edtmax,edtmin,3,edtc, &
655 ids,ide, jds,jde, kds,kde, &
656 ims,ime, jms,jme, kms,kme, &
657 its,ite, jts,jte, kts,kte)
658 do i=its,itf
659 if(ierr(i).eq.0)then
660 edt(i)=edtc(i,2)
661 endif
662 enddo
663 !
664 ! massflux from precip and normalized cloud properties
665 !
666 pwdper=0.
667 do i=its,itf
669 if(ierr(i).gt.0)pre(i)=0.
670 if(ierr(i).eq.0)then
671 xmb(i)=pre(i)/(pwav(i)+edt(i)*pwev(i))
672 !
673 !--- percent of that that is evaporated (pwd is negative)
674 !
675 if(i.eq.ipr.and.j.eq.jpr)then
676 print *,'xmb,edt,pwav = ',xmb(i),edt(i),pwav(i)
677 print *,'k,pwdper(i,k),pw,pwd(i,k)',z1(i)
678 endif
679 do k=1,ktop(i)
680 pwdper(i,k)=-edt(i)*pwd(i,k)/pwav(i)
681 if(i.eq.ipr.and.j.eq.jpr)then
682 print *,k,pwdper(i,k),pw(i,k),pwd(i,k)
683 endif
684 enddo
685 endif
686 enddo
687 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
688 !
689 !!!!! NOW WE HAVE EVREYTHING TO CALCULATE TRACER TRANSPORT AND WET DEPOSITION !!!
690 !
691 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
692 !
693 !--- calculate incloud tracer distribution
694 !
695 if(j.eq.jpr)print *,'calling up_tracer'
696 call cup_up_tracer(ierr,tcrit,t,pre,z_cup,p,tracer,tre_cup,tr_up,tr_pw, &
697 tr_c,hstary,pw,clw_all,kbcon,ktop,cd,mentr_rate,zu,k22,&
698 num_chem,ids,ide, jds,jde, kds,kde, &
699 ims,ime, jms,jme, kms,kme, &
700 its,ite, jts,jte, kts,kte,ipr,jpr,j,npr)
701 if(j.eq.jpr)print *,'called up_tracer'
702 call cup_dd_tracer(ierr,z_cup,qrcd,tracer,tre_cup,tr_up,tr_dd, &
703 tr_pw,tr_pwd,jmin,cdd,mentrd_rate,zd,pwdper,k22, &
704 num_chem,ids,ide, jds,jde, kds,kde, &
705 ims,ime, jms,jme, kms,kme, &
706 its,ite, jts,jte, kts,kte)
707 if(j.eq.jpr)print *,'called dd_tracer'
710 !
711 if(j.eq.jpr)then
712 i=ipr
713 print *,'in 250 loop ',edt(ipr),ierr(ipr)
714 ! if(ierr(i).eq.0.or.ierr(i).eq.3)then
715 print *,k22(I),kbcon(i),ktop(i),jmin(i)
716 print *,edt(i)
717 do k=kts,ktf
718 print *,k,z(i,k),he(i,k),hes(i,k)
719 enddo
720 do k=1,ktop(i)+1
721 print *,zu(i,k),zd(i,k),pw(i,k),pwd(i,k)
722 enddo
723 print *,'tr_up(i,k,6),tr_dd(i,k,6),tr_pw(i,k,6),tr_pwd(i,k,6)'
724 do k=1,ktop(i)+1
725 print *,tr_up(i,k,npr),tr_dd(i,k,npr),tr_pw(i,k,npr),tr_pwd(i,k,npr)
726 enddo
727 endif
728 ! endif
729 !
730 !--- calculate transport tendencies
731 !
732 !--- 1. in bottom layer
733 !
734 call cup_dellabot_tr(ipr,jpr,tre_cup,ierr,z_cup,p,tr_dd,edt, &
735 zd,cdd,tracer,tracert,j,mentrd_rate,z,g,xmb, &
736 num_chem,ids,ide, jds,jde, kds,kde, &
737 ims,ime, jms,jme, kms,kme, &
738 its,ite, jts,jte, kts,kte)
739 !
740 !--- 2. everywhere else
741 !
743 call cup_dellas_tr(ierr,z_cup,p_cup,tr_dd,edt,zd,cdd, &
744 tracer,tracert,j,mentrd_rate,zu,g,xmb, &
745 cd,tr_up,ktop,k22,kbcon,mentr_rate,jmin,tre_cup,kdet, &
746 k22,ipr,jpr,npr,'deep',num_chem, &
747 ids,ide, jds,jde, kds,kde, &
748 ims,ime, jms,jme, kms,kme, &
749 its,ite, jts,jte, kts,kte )
750 if(j.eq.jpr)then
751 i=ipr
752 do k=kts,ktf
753 print *,k,tracer(i,k,npr),tracert(i,k,npr)
754 enddo
755 endif
756 !
757 ! may need more below for wet deposition......
758 !
759 !
760 ! call cup_output_wd ( &
761 ! ids,ide, jds,jde, kds,kde, &
762 ! ims,ime, jms,jme, kms,kme, &
763 ! its,ite, jts,jte, kts,kte)
765 END SUBROUTINE CUP_CT
767 SUBROUTINE cup_dellabot_tr(ipr,jpr,tre_cup,ierr,z_cup,p_cup, &
768 tr_dd,edt,zd,cdd,tracer,tracert,j,mentrd_rate,z,g,xmb, &
769 num_chem,ids,ide, jds,jde, kds,kde, &
770 ims,ime, jms,jme, kms,kme, &
771 its,ite, jts,jte, kts,kte )
773 IMPLICIT NONE
775 integer &
776 ,intent (in ) :: &
777 num_chem,ids,ide, jds,jde, kds,kde, &
778 ims,ime, jms,jme, kms,kme, &
779 its,ite, jts,jte, kts,kte
780 integer, intent (in ) :: &
781 j,ipr,jpr
782 !
783 ! ierr error value, maybe modified in this routine
784 !
785 real, dimension (its:ite,kts:kte,1:num_chem) &
786 ,intent (out ) :: &
787 tracert
788 real, dimension (its:ite,kts:kte,1:num_chem) &
789 ,intent (in ) :: &
790 tre_cup,tracer,tr_dd
791 real, dimension (its:ite,kts:kte) &
792 ,intent (in ) :: &
793 z_cup,p_cup,zd,cdd,z
794 real, dimension (its:ite) &
795 ,intent (in ) :: &
796 edt,xmb
797 real &
798 ,intent (in ) :: &
799 g,mentrd_rate
800 integer, dimension (its:ite) &
801 ,intent (inout) :: &
802 ierr
803 !
804 ! local variables in this routine
805 !
807 integer i
808 real detdo,detdo1,detdo2,entdo,dp,dz,subin, &
809 totmas
810 !
811 integer :: itf, ktf, nv, npr
812 npr=24
813 itf=MIN(ite,ide-1)
814 ktf=MIN(kte,kde-1)
815 !
816 !
817 if(j.eq.jpr)print *,'in cup dellabot '
818 tracert=0.
819 do 100 i=its,itf
820 if(ierr(i).ne.0)go to 100
821 dz=z_cup(i,2)-z_cup(i,1)
822 DP=100.*(p_cup(i,1)-P_cup(i,2))
823 detdo1=edt(i)*zd(i,2)*CDD(i,1)*DZ
824 detdo2=edt(i)*zd(i,1)
825 entdo=edt(i)*zd(i,2)*mentrd_rate*dz
826 subin=-EDT(I)*zd(i,2)
827 detdo=detdo1+detdo2-entdo+subin
828 do nv=1,num_chem
829 tracert(I,1,nv)=(detdo1*.5*(tr_dd(i,1,nv)+tr_dd(i,2,nv)) &
830 +detdo2*tr_dd(i,1,nv) &
831 +subin*tre_cup(i,2,nv) &
832 -entdo*tracer(i,1,nv))*g/dp*xmb(i)
833 enddo
834 if(j.eq.jpr.and.i.eq.ipr)print *,'in cup dellabot ',tracert(I,1,npr), &
835 detdo1,detdo2,subin,entdo,tr_dd(i,1,npr),tr_dd(i,2,npr),tracer(i,1,npr)
836 100 CONTINUE
838 END SUBROUTINE cup_dellabot_tr
841 SUBROUTINE cup_dellas_tr(ierr,z_cup,p_cup,tr_dd,edt,zd,cdd, &
842 tracer,tracert,j,mentrd_rate,zu,g,xmb, &
843 cd,tr_up,ktop,k22,kbcon,mentr_rate,jmin,tre_cup,kdet,kpbl, &
844 ipr,jpr,npr,name,num_chem, &
845 ids,ide, jds,jde, kds,kde, &
846 ims,ime, jms,jme, kms,kme, &
847 its,ite, jts,jte, kts,kte )
849 IMPLICIT NONE
851 integer &
852 ,intent (in ) :: &
853 num_chem,ids,ide, jds,jde, kds,kde, &
854 ims,ime, jms,jme, kms,kme, &
855 its,ite, jts,jte, kts,kte
856 integer, intent (in ) :: &
857 j,ipr,jpr,npr
858 !
859 ! ierr error value, maybe modified in this routine
860 !
861 real, dimension (its:ite,kts:kte,1:num_chem) &
862 ,intent (inout ) :: &
863 tracert
864 real, dimension (its:ite,kts:kte,1:num_chem) &
865 ,intent (in ) :: &
866 tr_up,tr_dd,tre_cup,tracer
867 real, dimension (its:ite,kts:kte) &
868 ,intent (in ) :: &
869 z_cup,p_cup,zd,cdd,cd,zu
870 real, dimension (its:ite) &
871 ,intent (in ) :: &
872 edt,xmb
873 real &
874 ,intent (in ) :: &
875 g,mentrd_rate,mentr_rate
876 integer, dimension (its:ite) &
877 ,intent (in ) :: &
878 kbcon,ktop,k22,jmin,kdet,kpbl
879 integer, dimension (its:ite) &
880 ,intent (inout) :: &
881 ierr
882 character *(*), intent (in) :: &
883 name
884 !
885 ! local variables in this routine
886 !
888 integer i,k,nv
889 real detdo1,detdo2,entdo,dp,dz,subin,detdo,entup, &
890 detup,subdown,entdoj,entupk,detupk,totmas
891 !
892 integer :: itf, ktf
893 ! npr=24
894 itf=MIN(ite,ide-1)
895 ktf=MIN(kte,kde-1)
896 !
897 !
898 i=ipr
899 if(j.eq.jpr)then
900 print *,'in dellas kpbl(i),k22(i),kbcon(i),ktop(i),jmin(i)'
901 print *,kpbl(i),k22(i),kbcon(i),ktop(i),jmin(i)
902 endif
903 do nv=1,num_chem
904 DO K=kts+1,kte
905 do i=its,itf
906 tracert(i,k,nv)=0.
907 enddo
908 enddo
909 enddo
910 !
911 DO 100 k=kts+1,ktf-1
912 DO 100 i=its,ite
913 IF(ierr(i).ne.0)GO TO 100
914 IF(K.Gt.KTOP(I))GO TO 100
915 !
916 !--- SPECIFY DETRAINMENT OF DOWNDRAFT, HAS TO BE CONSISTENT
917 !--- WITH ZD CALCULATIONS IN SOUNDD.
918 !
919 DZ=Z_cup(I,K+1)-Z_cup(I,K)
920 detdo=edt(i)*CDD(i,K)*DZ*ZD(i,k+1)
921 entdo=edt(i)*mentrd_rate*dz*zd(i,k+1)
922 subin=zu(i,k+1)-zd(i,k+1)*edt(i)
923 entup=0.
924 detup=0.
925 if(k.ge.kbcon(i).and.k.lt.ktop(i))then
926 entup=mentr_rate*dz*zu(i,k)
927 detup=CD(i,K+1)*DZ*ZU(i,k)
928 endif
929 subdown=(zu(i,k)-zd(i,k)*edt(i))
930 entdoj=0.
931 entupk=0.
932 detupk=0.
933 !
934 if(k.eq.jmin(i))then
935 entdoj=edt(i)*zd(i,k)
936 endif
938 if(k.eq.k22(i)-1)then
939 entupk=zu(i,kpbl(i))
940 endif
942 if(k.gt.kdet(i))then
943 detdo=0.
944 endif
946 if(k.eq.ktop(i)-0)then
947 detupk=zu(i,ktop(i))
948 subin=0.
949 endif
950 if(k.lt.kbcon(i))then
951 detup=0.
952 endif
953 !C
954 !C--- CHANGED DUE TO SUBSIDENCE AND ENTRAINMENT
955 !C
956 totmas=subin-subdown+detup-entup-entdo+ &
957 detdo-entupk-entdoj+detupk
958 if(j.eq.jpr.and.i.eq.ipr)print *,'k,totmas,sui,sud = ',k, &
959 totmas,subin,subdown
960 ! if(j.eq.jpr.and.i.eq.ipr)print *,'updr stuff = ',detup,
961 ! 1 entup,entupk,detupk
962 ! if(j.eq.jpr.and.i.eq.ipr)print *,'dddr stuff = ',entdo,
963 ! 1 detdo,entdoj
964 if(abs(totmas).gt.1.e-6)then
965 print *,'*********************',i,j,k,totmas,name
966 print *,kpbl(i),k22(i),kbcon(i),ktop(i)
967 !c print *,'updr stuff = ',subin,
968 !c 1 subdown,detup,entup,entupk,detupk
969 !c print *,'dddr stuff = ',entdo,
970 !c 1 detdo,entdoj
971 CALL wrf_error_fatal ( 'cup_dellas_tr: TOTMAS > CRITICAL VALUE')
972 endif
973 dp=100.*(p_cup(i,k-1)-p_cup(i,k))
974 do nv=1,num_chem
975 ! tracert(i,k,nv)=(subin*tre_cup(i,k+1,nv) &
976 ! -subdown*tre_cup(i,k,nv) &
977 tracert(i,k,nv)=(subin*tracer(i,k+1,nv) &
978 -subdown*tracer(i,k,nv) &
979 +detup*.5*(tr_up(i,K+1,nv)+tr_up(i,K,nv)) &
980 +detdo*.5*(tr_dd(i,K+1,nv)+tr_dd(i,K,nv)) &
981 -entup*tracer(i,k,nv) &
982 -entdo*tracer(i,k,nv) &
983 -entupk*tre_cup(i,k22(i),nv) &
984 -entdoj*tre_cup(i,jmin(i),nv) &
985 +detupk*tr_up(i,ktop(i),nv) &
986 )*g/dp*xmb(i)
987 enddo
988 if(i.eq.ipr.and.j.eq.jpr)then
989 print *,k,tracert(i,k,npr),subin*tre_cup(i,k+1,npr),subdown*tre_cup(i,k,npr), &
990 detdo*.5*(tr_dd(i,K+1,npr)+tr_dd(i,K,npr))
991 print *,k,detup*.5*(tr_up(i,K+1,npr)+tr_up(i,K,npr)),detupk*tr_up(i,ktop(i),npr), &
992 entup*tracer(i,k,npr),entdo*tracer(i,k,npr)
993 print *,k,entupk*tre_cup(i,k,npr),detupk,tr_up(i,ktop(i),npr)
994 endif
996 100 CONTINUE
998 END SUBROUTINE cup_dellas_tr
999 SUBROUTINE cup_env_clev_tr(tracer,tre_cup,num_chem,ierr, &
1000 ids,ide, jds,jde, kds,kde, &
1001 ims,ime, jms,jme, kms,kme, &
1002 its,ite, jts,jte, kts,kte)
1003 implicit none
1004 integer &
1005 ,intent (in ) :: &
1006 num_chem,ids,ide, jds,jde, kds,kde, &
1007 ims,ime, jms,jme, kms,kme, &
1008 its,ite, jts,jte, kts,kte
1009 integer, dimension (its:ite) &
1010 ,intent (in) :: &
1011 ierr
1013 real, dimension (its:ite,kts:kte,1:num_chem) &
1014 ,intent (in ) :: &
1015 tracer
1016 real, dimension (its:ite,kts:kte,1:num_chem) &
1017 ,intent (out ) :: &
1018 tre_cup
1019 !
1020 ! local variables in this routine
1021 !
1023 integer :: &
1024 i,k,nv,itf,ktf
1025 itf=MIN(ite,ide-1)
1026 ktf=MIN(kte,kde-1)
1027 do nv=1,num_chem
1028 do k=kts+1,ktf
1029 do i=its,ite
1030 if(ierr(i).eq.0)then
1031 tre_cup(i,k,nv)=.5*(tracer(i,k-1,nv)+tracer(i,k,nv))
1032 endif
1033 enddo
1034 enddo
1035 enddo
1036 do nv=1,num_chem
1037 do i=its,ite
1038 if(ierr(i).eq.0)then
1039 tre_cup(i,kts,nv)=tracer(i,kts,nv)
1040 endif
1041 enddo
1042 enddo
1045 END subroutine cup_env_clev_tr
1048 SUBROUTINE cup_up_tracer(ierr,tcrit,t,pre,z_cup,p,tracer,tre_cup,tr_up, &
1049 tr_pw,tr_c,hstary,cupclw,clw_all,kbcon,ktop,cd,mentr_rate,zu,k22, &
1050 num_cc,ids,ide, jds,jde, kds,kde, &
1051 ims,ime, jms,jme, kms,kme, &
1052 its,ite, jts,jte, kts,kte,ipr,jpr,j,npr)
1053 USE module_configure
1054 USE module_state_description
1055 USE module_ctrans_aqchem
1056 implicit none
1057 ! Aqeuous species pointers INCLUDE File
1059 !...........PARAMETERS and their descriptions:
1061 INTEGER NGAS ! number of gas phase species for AQCHEM
1062 PARAMETER ( NGAS = 11 )
1064 INTEGER NAER ! number of aerosol species for AQCHEM
1065 PARAMETER ( NAER = 23 )
1067 !...pointers for the AQCHEM array GAS
1069 INTEGER LSO2 ! local pointer to SO2
1070 PARAMETER ( LSO2 = 1 )
1072 INTEGER LHNO3 ! local pointer to HNO3
1073 PARAMETER ( LHNO3 = 2 )
1075 INTEGER LN2O5 ! local pointer to N2O5
1076 PARAMETER ( LN2O5 = 3 )
1077 INTEGER LCO2 ! local pointer to CO2
1078 PARAMETER ( LCO2 = 4 )
1080 INTEGER LNH3 ! local pointer to NH3
1081 PARAMETER ( LNH3 = 5 )
1083 INTEGER LH2O2 ! local pointer to H2O2
1084 PARAMETER ( LH2O2 = 6 )
1086 INTEGER LO3 ! local pointer to O3
1087 PARAMETER ( LO3 = 7 )
1089 INTEGER LFOA ! local pointer to FOA
1090 PARAMETER ( LFOA = 8 )
1092 INTEGER LMHP ! local pointer to MHP
1093 PARAMETER ( LMHP = 9 )
1095 INTEGER LPAA ! local pointer to PAA
1096 PARAMETER ( LPAA = 10 )
1098 INTEGER LH2SO4 ! local pointer to H2SO4
1099 PARAMETER ( LH2SO4 = 11 )
1101 !...pointers for the AQCHEM array AEROSOL
1103 INTEGER LSO4AKN ! local pointer to SO4I aerosol
1104 PARAMETER ( LSO4AKN = 1 )
1106 INTEGER LSO4ACC ! local pointer to SO4 aerosol
1107 PARAMETER ( LSO4ACC = 2 )
1109 INTEGER LNH4AKN ! local pointer to NH4I aerosol
1110 PARAMETER ( LNH4AKN = 3 )
1112 INTEGER LNH4ACC ! local pointer to NH4 aerosol
1113 PARAMETER ( LNH4ACC = 4 )
1115 INTEGER LNO3AKN ! local pointer to NO3I aerosol
1116 PARAMETER ( LNO3AKN = 5 )
1118 INTEGER LNO3ACC ! local pointer to NO3 aerosol
1119 PARAMETER ( LNO3ACC = 6 )
1121 INTEGER LNO3COR ! local pointer to course aerosol nitrate
1122 PARAMETER ( LNO3COR = 7 )
1124 INTEGER LORGAKN ! local pointer to organic I aerosol
1125 PARAMETER ( LORGAKN = 8 )
1127 INTEGER LORGACC ! local pointer to organic aerosol
1128 PARAMETER ( LORGACC = 9 )
1130 INTEGER LPRIAKN ! local pointer to primary I aerosol
1131 PARAMETER ( LPRIAKN = 10 )
1133 INTEGER LPRIACC ! local pointer to primary aerosol
1134 PARAMETER ( LPRIACC = 11 )
1136 INTEGER LPRICOR ! local pointer to primary I aerosol
1137 PARAMETER ( LPRICOR = 12 )
1139 INTEGER LCACO3 ! local pointer to CaCO3 aerosol
1140 PARAMETER ( LCACO3 = 13 )
1142 INTEGER LMGCO3 ! local pointer to MgCO3 aerosol
1143 PARAMETER ( LMGCO3 = 14 )
1145 INTEGER LNACL ! local pointer to NaCl aerosol
1146 PARAMETER ( LNACL = 15 )
1148 INTEGER LA3FE ! local pointer to Fe+++ aerosol
1149 PARAMETER ( LA3FE = 16 )
1151 INTEGER LB2MN ! local pointer to Mn++ aerosol
1152 PARAMETER ( LB2MN = 17 )
1154 INTEGER LKCL ! local pointer to NaCl aerosol
1155 PARAMETER ( LKCL = 18 )
1157 INTEGER LNUMAKN ! local pointer to # Aitken aerosol
1158 PARAMETER ( LNUMAKN = 19 )
1160 INTEGER LNUMACC ! local pointer to # accumulation aerosol
1161 PARAMETER ( LNUMACC = 20 )
1163 INTEGER LNUMCOR ! local pointer to # coarse aerosol
1164 PARAMETER ( LNUMCOR = 21 )
1166 INTEGER LSRFAKN ! local pointer to sfc area Aitken aerosol
1167 PARAMETER ( LSRFAKN = 22 )
1169 INTEGER LSRFACC ! local pntr to sfc area accumulation aerosol
1170 PARAMETER ( LSRFACC = 23 )
1173 !
1174 ! on input
1175 !
1177 ! only local wrf dimensions are need as of now in this routine
1179 integer &
1180 ,intent (in ) :: &
1181 num_cc,ids,ide, jds,jde, kds,kde, &
1182 ims,ime, jms,jme, kms,kme, &
1183 its,ite, jts,jte, kts,kte,ipr,jpr,j,npr
1184 real, dimension (its:ite,kts:kte) &
1185 ,intent (in ) :: &
1186 z_cup,cd,zu,p,hstary,t
1187 real, dimension (its:ite,kts:kte) &
1188 ,intent (inout ) :: &
1189 cupclw,clw_all
1190 real, dimension (its:ite,kts:kte,1:num_chem) &
1191 ,intent (inout ) :: &
1192 tr_up,tr_c,tr_pw
1193 real, dimension (its:ite,kts:kte,1:num_chem) &
1194 ,intent (in ) :: &
1195 tre_cup,tracer
1196 real, dimension (its:ite) &
1197 ,intent (in ) :: &
1198 pre
1200 ! entr= entrainment rate
1201 real &
1202 ,intent (in ) :: &
1203 mentr_rate,tcrit
1204 integer, dimension (its:ite) &
1205 ,intent (in ) :: &
1206 kbcon,ktop,k22
1207 ! ierr error value, maybe modified in this routine
1209 integer, dimension (its:ite) &
1210 ,intent (inout) :: &
1211 ierr
1212 ! local variables in this routine
1213 !
1214 real :: conc_equi,conc_mxr,partialp,taucld
1216 integer :: &
1217 iall,i,k,iwd,nv
1218 real :: &
1219 dh,qrch,c0,dz,radius,airm,dens
1220 integer :: &
1221 itf,ktf,iaer,igas
1222 !
1223 ! aerosol scavenging coeffs for aitken mode
1224 !
1225 real alfa0,alfa2,alfa3
1226 ! output variables
1227 ! hpwdep h+ deposition
1228 real, dimension (ngas) :: gas,gaswdep
1229 real, dimension (naer) :: aerosol,aerwdep
1230 real hpwdep
1231 alfa0=0.
1232 alfa2=0.
1233 alfa3=0.
1234 gas(lco2)=340.
1235 taucld=1800.
1236 qrch=0.
1237 itf=MIN(ite,ide-1)
1238 ktf=MIN(kte,kde-1)
1240 !
1241 iall=0
1242 c0=.002
1243 iwd=0
1244 !
1245 !--- no precip for small clouds
1246 !
1247 if(mentr_rate.gt.0.)then
1248 radius=.2/mentr_rate
1249 if(radius.lt.900.)c0=0.
1250 ! if(radius.lt.900.)iall=0
1251 endif
1252 do nv=1,num_chem
1253 do k=kts,ktf
1254 do i=its,itf
1255 tr_pw(i,k,nv)=0.
1256 tr_up(i,k,nv)=tre_cup(i,k,nv)
1257 tr_c(i,k,nv)=0.
1258 enddo
1259 enddo
1260 enddo
1261 do nv=1,num_chem
1262 do i=its,itf
1263 if(ierr(i).eq.0.)then
1264 do k=k22(i),kbcon(i)-1
1265 tr_up(i,k,nv)=tre_cup(i,k22(i),nv)
1266 enddo
1267 endif
1268 enddo
1269 enddo
1270 if(j.eq.jpr)print *,'p_so2,o_o3 = ',p_so2,p_o3
1271 DO 100 k=kts+1,kte-1
1272 DO 100 i=its,itf
1273 AEROSOL=0.
1274 GAS=0.
1275 IF(ierr(i).ne.0)GO TO 100
1276 IF(K.Lt.KBCON(I))GO TO 100
1277 IF(K.Gt.KTOP(I)+1)GO TO 100
1278 DZ=Z_cup(i,K)-Z_cup(i,K-1)
1279 if(cupclw(i,k).le.0.)cupclw(i,k)=0.
1280 if(clw_all(i,k).le.0.)clw_all(i,k)=0.
1281 !
1282 !------ 1. steady state plume equation, for what could
1283 !------ be in cloud before anything happens (kg/kg)
1284 !------ tr_up would be the concentration if tr would be conserved
1285 !
1286 !
1287 do nv=1,num_chem
1288 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)
1289 tr_up(i,K,nv)=(tr_up(i,K-1,nv)*(1.-.5*CD(i,K)*DZ)+mentr_rate* &
1290 DZ*tracer(i,K-1,nv))/(1.+mentr_rate*DZ-.5*cd(i,k)*dz)
1291 enddo
1292 !
1293 ! sources or sinks due to aq chem
1294 !
1296 dens=1000.*p(i,k)*100./t(i,k)/287./28.9628
1297 airm=dens*dz
1298 !...gas concentrations (ppm)
1300 GAS( LCO2 ) = 370.0
1301 GAS( LFOA ) = 0.0 ! ???
1302 GAS( LMHP ) = 0.0 ! ???
1304 GAS( LSO2 ) = tr_up(i,k,p_so2)
1305 GAS( LH2SO4 ) = tr_up(i,k,p_sulf)
1306 GAS( LNH3 ) = tr_up(i,k,p_nh3)
1307 GAS( LH2O2 ) = tr_up(i,k,p_h2o2)
1309 GAS( LO3 ) = tr_up(i,k,p_o3)
1310 GAS( LPAA ) = tr_up(i,k,p_paa)
1311 GAS( LHNO3 ) = tr_up(i,k,p_hno3)
1312 GAS( LN2O5 ) = tr_up(i,k,p_n2o5)
1313 !...convert to mol/mol
1315 DO IGAS=1,NGAS
1316 GAS( IGAS ) = GAS( IGAS ) * 1.0E-6
1317 END DO
1319 !...aerosol concentrations (ug/m3)
1321 ! AEROSOL( LSO4ACC ) = 20.0
1322 ! AEROSOL( LNH4ACC ) = 6.65
1323 ! AEROSOL( LNO3ACC ) = 10.0
1324 ! AEROSOL( LNACL ) = 1.71
1325 !! AEROSOL( LA3FE ) = 0.5
1326 ! AEROSOL( LB2MN ) = 0.02
1327 ! AEROSOL( LNO3COR ) = 0.0
1328 AEROSOL( LORGACC ) = 0.0
1329 AEROSOL( LPRIACC ) = 0.0
1330 ! AEROSOL( LCACO3 ) = 3.05
1331 ! AEROSOL( LMGCO3 ) = 0.0
1333 AEROSOL( LSO4ACC ) = tr_up(i,k,p_so4aj)
1334 AEROSOL( LNH4ACC ) = tr_up(i,k,p_nh4aj)
1335 AEROSOL( LNO3ACC ) = tr_up(i,k,p_no3aj)
1336 AEROSOL( LNACL ) = 0.
1337 AEROSOL( LA3FE ) = .5
1338 AEROSOL( LB2MN ) = .02
1339 AEROSOL( LNO3COR ) = 0.
1340 ! AEROSOL( LORGACC ) = tr_up(i,k,) + tr_up(i,k,) + tr_up(i,k,)
1341 ! AEROSOL( LPRIACC ) = tr_up(i,k,) + tr_up(i,k,)
1342 AEROSOL( LCACO3 ) = 0.
1343 AEROSOL( LMGCO3 ) = 0.
1346 !...convert to mol/mol
1347 !
1349 ! DO IAER=1,NAER
1350 ! AEROSOL( IAER ) = AEROSOL( IAER ) * 1.0E-6 * CTHK1
1351 ! & / ( SGRAERMW( IAER ) * AIRM )
1352 ! END DO
1353 DO IAER=1,NAER
1354 AEROSOL( IAER ) = AEROSOL( IAER ) * 1.0E-6
1355 END DO
1356 ! first clw is water, second is total
1358 GASWDEP=0.
1359 AERWDEP=0.
1360 HPWDEP=0.
1361 ! if(clw_all(i,k).gt.1.e-12)then
1362 ! if(cupclw(i,k).gt.1.e-12)then
1363 ! CALL AQCHEM (t(i,k),p(i,k)*100.,taucld,cupclw(i,k)/3600., &
1364 ! clw_all(i,k)*dens,clw_all(i,k)*dens,airm,ALFA0,ALFA2,ALFA3,GAS, &
1365 ! AEROSOL, GASWDEP, AERWDEP, HPWDEP )
1366 ! endif
1367 ! endif
1373 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1374 !
1375 ! FOLLOWING FOR WET DEPOSITION !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1376 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1377 !
1378 do nv=1,num_chem
1379 tr_c(i,k,nv)=0.
1380 tr_pw(i,k,nv)=c0*dz*tr_C(I,K,nv)*zu(i,k)
1381 if(tr_c(i,k,nv).le.0.)then
1382 tr_c(i,k,nv)=0.
1383 endif
1384 enddo
1385 !
1386 !--- iall.eq.1, if all cloudwater goes to rain
1387 !
1388 if(iall.eq.1)then
1389 tr_c(i,k,nv)=0.
1390 tr_pw(i,k,nv)=(tr_c(I,K,nv)-QRCH)*zu(i,k)
1391 if(tr_pw(i,k,nv).lt.0.)tr_pw(i,k,nv)=0.
1392 endif
1394 !
1395 !----- set next level
1396 ! tr_up(I,K,nv)=tr_c(I,K,nv)+qrch
1397 tr_up(i,k,p_so2)=gas(lso2)*1.e6
1398 tr_up(i,k,p_sulf)=gas(lh2so4)*1.e6
1399 tr_up(i,k,p_nh3)=gas(lnh3)*1.e6
1400 tr_up(i,k,p_h2o2)=gas(lh2o2)*1.e6
1402 tr_up(i,k,p_o3)=gas(lo3)*1.e6
1403 tr_up(i,k,p_paa)=gas(lpaa)*1.e6
1404 tr_up(i,k,p_hno3)=gas(lhno3)*1.e6
1405 tr_up(i,k,p_n2o5)=gas(ln2o5)*1.e6
1406 tr_up(i,k,p_so4aj)=AEROSOL( LSO4ACC )*1.e6
1407 tr_up(i,k,p_nh4aj)=AEROSOL( LNH4ACC )*1.e6
1408 tr_up(i,k,p_no3aj)=AEROSOL( LNO3ACC ) *1.e6
1410 tr_pw(i,k,p_so2)=gaswdep(lso2)*1.e6
1411 tr_pw(i,k,p_sulf)=gaswdep(lh2so4)*1.e6
1412 tr_pw(i,k,p_nh3)=gaswdep(lnh3)*1.e6
1413 tr_pw(i,k,p_h2o2)=gaswdep(lh2o2)*1.e6
1415 tr_pw(i,k,p_o3)=gaswdep(lo3)*1.e6
1416 tr_pw(i,k,p_paa)=gaswdep(lpaa)*1.e6
1417 tr_pw(i,k,p_hno3)=gaswdep(lhno3)*1.e6
1418 tr_pw(i,k,p_n2o5)=gaswdep(ln2o5)*1.e6
1419 tr_pw(i,k,p_so4aj)=AERwdep( LSO4ACC )*1.e6
1420 tr_pw(i,k,p_nh4aj)=AERwdep( LNH4ACC )*1.e6
1421 tr_pw(i,k,p_no3aj)=AERwdep( LNO3ACC ) *1.e6
1422 if(i.eq.ipr.and.j.eq.jpr)then
1423 write(6,*)'a',tr_up(i,k,npr),tracer(i,K-1,npr),tr_pw(i,k,npr)
1424 endif
1426 100 CONTINUE
1429 END subroutine cup_up_tracer
1433 SUBROUTINE cup_dd_tracer(ierr,z_cup,qrcd,tracer,tre_cup,tr_up,tr_dd, &
1434 tr_pw,tr_pwd,jmin,cdd,entr,zd,pwdper,k22, &
1435 numch,ids,ide, jds,jde, kds,kde, &
1436 ims,ime, jms,jme, kms,kme, &
1437 its,ite, jts,jte, kts,kte)
1438 USE module_configure
1439 USE module_state_description
1440 implicit none
1441 !
1442 ! on input
1443 !
1445 ! only local wrf dimensions are need as of now in this routine
1447 integer &
1448 ,intent (in ) :: &
1449 numch,ids,ide, jds,jde, kds,kde, &
1450 ims,ime, jms,jme, kms,kme, &
1451 its,ite, jts,jte, kts,kte
1452 real, dimension (its:ite,kts:kte) &
1453 ,intent (in ) :: &
1454 pwdper,zd,cdd,qrcd,z_cup
1455 real, dimension (its:ite,kts:kte,1:num_chem) &
1456 ,intent (inout ) :: &
1457 tr_dd,tr_pwd,tr_up
1458 real, dimension (its:ite,kts:kte,1:num_chem) &
1459 ,intent (in ) :: &
1460 tre_cup,tracer,tr_pw
1461 real, dimension (its:ite,1:num_chem) :: pwav
1463 ! entr= entrainment rate
1464 real &
1465 ,intent (in ) :: &
1466 entr
1467 integer, dimension (its:ite) &
1468 ,intent (in ) :: &
1469 jmin
1470 ! ierr error value, maybe modified in this routine
1472 integer, dimension (its:ite) &
1473 ,intent (inout) :: &
1474 ierr,k22
1475 ! local variables in this routine
1476 !
1478 integer :: &
1479 iall,i,k,nv,ki
1480 real :: &
1481 dh,qrch,c0,dz,radius
1482 integer :: &
1483 itf,ktf
1484 logical iaer (num_chem)
1485 iaer = .false.
1487 iaer(p_so4aj) = .true.
1488 iaer(p_nh4aj) = .true.
1489 iaer(p_no3aj) = .true.
1491 itf=MIN(ite,ide-1)
1492 ktf=MIN(kte,kde-1)
1493 !
1494 qrch=0.
1495 do nv=1,num_chem
1496 do k=kts+1,kte
1497 do i=its,ite
1498 tr_dd(i,k,nv)=0.
1499 tr_pwd(i,k,nv)=0.
1500 enddo
1501 enddo
1502 do i=its,ite
1503 pwav(i,nv)=0.
1504 IF(ierr(I).eq.0)then
1505 do k=kts,ktf
1506 pwav(i,nv)=pwav(i,nv)+tr_pw(i,k,nv)
1507 enddo
1508 endif
1509 enddo
1510 enddo
1511 !
1512 !--- in downdraft, do only transport of tracers, other
1513 !--- than evaporation of part of the rainwater (see below)
1514 !
1515 !
1516 do 100 i=its,ite
1517 IF(ierr(I).eq.0)then
1518 !
1519 !--- assume no gas takeup by rain during falling
1520 !--- for now
1521 !
1522 !
1523 do nv=1,num_chem
1524 tr_dd(i,jmin(i),nv)=tre_cup(i,jmin(i),nv)
1525 enddo
1526 do ki=jmin(i)-1,1,-1
1527 DZ=Z_cup(i,Ki+1)-Z_cup(i,Ki)
1528 do nv=1,num_chem
1529 tr_pwd(i,jmin(i),nv)=0.
1530 tr_dd(i,Ki,nv)=(tr_dd(i,Ki+1,nv)*(1.-.5*CDD(i,Ki)*DZ) &
1531 +entr*DZ*tracer(i,Ki,nv) &
1532 )/(1.+entr*DZ-.5*CDD(i,Ki)*DZ)
1533 !
1534 !--- if tracer conserved
1535 !
1536 qrch=tr_dd(i,Ki,nv)
1537 !
1538 !--- part of dissolved liquid phase material that is being evaporated
1539 ! need percentage of rainwater that evaporates at level
1540 ! pwdper
1541 ! qcd=qcd+pwdper
1542 !
1543 ! tr_pwd(i,ki,nv)=pwdper(i,ki)*pwav(i,nv)
1544 if(iaer(nv))then
1545 tr_pwd(i,ki,nv)=0.
1546 else
1547 tr_pwd(i,ki,nv)=pwdper(i,ki)*pwav(i,nv)
1548 endif
1549 tr_dd(i,ki,nv)=qrch+tr_pwd(i,ki,nv)
1550 enddo
1551 !
1552 !--- end loop over nv
1553 enddo
1554 endif
1555 100 continue
1557 END subroutine cup_dd_tracer
1562 SUBROUTINE neg_check_ct(pret,ktop,epsilc,dt,q,outq,iopt,num_chem, &
1563 its,ite,kts,kte,itf,ktf,ipr,jpr,npr,j)
1565 INTEGER, INTENT(IN ) :: iopt,num_chem,its,ite,kts,kte,itf,ktf,ipr,jpr,npr,j
1567 real, dimension (its:ite,kts:kte,num_chem ) , &
1568 intent(inout ) :: &
1569 q,outq
1570 real, dimension (its:ite ) , &
1571 intent(in ) :: &
1572 pret
1573 integer, dimension (its:ite ) , &
1574 intent(in ) :: &
1575 ktop
1576 real &
1577 ,intent (in ) :: &
1578 dt,epsilc
1579 real :: tracermin,tracermax,thresh,qmem,qmemf,qmem2,qtest,qmem1
1580 !
1581 ! check whether routine produces negative q's. This can happen, since
1582 ! tendencies are calculated based on forced q's. This should have no
1583 ! influence on conservation properties, it scales linear through all
1584 ! tendencies. Use iopt=0 to test for each tracer seperately, iopt=1
1585 ! for a more severe limitation...
1586 !
1587 thresh=epsilc
1588 ! thresh=1.e-30
1589 if(iopt.eq.0)then
1590 do nv=1,num_chem
1591 do 100 i=its,itf
1592 if(pret(i).le.0.)go to 100
1593 tracermin=q(i,kts,nv)
1594 tracermax=q(i,kts,nv)
1595 do k=kts+1,kte-1
1596 tracermin=min(tracermin,q(i,k,nv))
1597 tracermax=max(tracermax,q(i,k,nv))
1598 enddo
1599 tracermin=max(tracermin,thresh)
1600 qmemf=1.
1601 !
1602 ! first check for minimum restriction
1603 !
1604 do k=kts,ktop(i)
1605 !
1606 ! tracer tendency
1607 !
1608 qmem=outq(i,k,nv)
1609 !
1610 ! only necessary if there is a tendency
1611 !
1612 if(qmem.lt.0.)then
1613 qtest=q(i,k,nv)+outq(i,k,nv)*dt
1614 if(qtest.lt.tracermin)then
1615 !
1616 ! qmem2 would be the maximum allowable tendency
1617 !
1618 qmem1=outq(i,k,nv)
1619 qmem2=(tracermin-q(i,k,nv))/dt
1620 qmemf=min(qmemf,qmem2/qmem1)
1621 if(qmemf.gt.1.)print *,'something wrong in negct_1',qmem2,qmem1
1622 if(i.eq.ipr.and.j.eq.jpr.and.nv.eq.npr)then
1623 print *,k,qtest,qmem2,qmem1,qmemf
1624 endif
1625 qmemf=max(qmemf,0.)
1626 endif
1627 endif
1628 enddo
1629 do k=kts,ktop(i)
1630 outq(i,k,nv)=outq(i,k,nv)*qmemf
1631 enddo
1632 !
1633 ! now check max
1634 !
1635 qmemf=1.
1636 do k=kts,ktop(i)
1637 !
1638 ! tracer tendency
1639 !
1640 qmem=outq(i,k,nv)
1641 !
1642 ! only necessary if there is a tendency
1643 !
1644 if(qmem.gt.0.)then
1645 qtest=q(i,k,nv)+outq(i,k,nv)*dt
1646 if(qtest.gt.tracermax)then
1647 !
1648 ! qmem2 would be the maximum allowable tendency
1649 !
1650 qmem1=outq(i,k,nv)
1651 qmem2=(tracermax-q(i,k,nv))/dt
1652 qmemf=min(qmemf,qmem2/qmem1)
1653 if(qmemf.gt.1.)print *,'something wrong in negct_2',qmem2,qmem1
1654 if(i.eq.ipr.and.j.eq.jpr.and.nv.eq.npr)then
1655 print *,'2',k,qtest,qmem2,qmem1,qmemf
1656 endif
1657 qmemf=max(qmemf,0.)
1658 endif
1659 endif
1660 enddo
1661 do k=kts,ktop(i)
1662 outq(i,k,nv)=outq(i,k,nv)*qmemf
1663 enddo
1664 100 continue
1665 enddo
1666 !
1667 ! ELSE
1668 !
1669 elseif(iopt.eq.1)then
1670 do i=its,itf
1671 qmemf=1.
1672 do k=kts,ktop(i)
1673 do nv=1,num_chem
1674 !
1675 ! tracer tendency
1676 !
1677 qmem=outq(i,k,nv)
1678 !
1679 ! only necessary if tendency is larger than zero
1680 !
1681 if(qmem.lt.0.)then
1682 qtest=q(i,k,nv)+outq(i,k,nv)*dt
1683 if(qtest.lt.thresh)then
1684 !
1685 ! qmem2 would be the maximum allowable tendency
1686 !
1687 qmem1=outq(i,k,nv)
1688 qmem2=(thresh-q(i,k,nv))/dt
1689 qmemf=min(qmemf,qmem2/qmem1)
1690 qmemf=max(0.,qmemf)
1691 endif
1692 endif
1693 enddo
1694 enddo
1695 do nv=1,num_chem
1696 do k=kts,ktop(i)
1697 outq(i,k,nv)=outq(i,k,nv)*qmemf
1698 enddo
1699 enddo
1700 enddo
1701 endif
1703 END SUBROUTINE neg_check_ct
1706 !-------------------------------------------------------
1707 END MODULE module_ctrans_grell