| blob | 30e9cd02d96bb70688e6a3770ee67558f59aebde |
1 !!
2 MODULE module_cu_sas
4 CONTAINS
6 !-----------------------------------------------------------------
7 SUBROUTINE CU_SAS(DT,ITIMESTEP,STEPCU, &
8 RTHCUTEN,RQVCUTEN,RQCCUTEN,RQICUTEN, &
9 RUCUTEN,RVCUTEN, & ! gopal's doing for SAS
10 RAINCV,PRATEC,HTOP,HBOT, &
11 U3D,V3D,W,T3D,QV3D,QC3D,QI3D,PI3D,RHO3D, &
12 DZ8W,PCPS,P8W,XLAND,CU_ACT_FLAG, &
13 P_QC, &
14 #if (NMM_CORE == 1)
15 STORE_RAND,MOMMIX, & ! gopal's doing
16 #endif
17 P_QI,P_FIRST_SCALAR, &
18 CUDT, CURR_SECS, ADAPT_STEP_FLAG, &
19 ids,ide, jds,jde, kds,kde, &
20 ims,ime, jms,jme, kms,kme, &
21 its,ite, jts,jte, kts,kte )
23 !-------------------------------------------------------------------
24 USE MODULE_GFS_MACHINE , ONLY : kind_phys
25 USE MODULE_GFS_FUNCPHYS , ONLY : gfuncphys
26 USE MODULE_GFS_PHYSCONS, grav => con_g, CP => con_CP, HVAP => con_HVAP &
27 &, RV => con_RV, FV => con_fvirt, T0C => con_T0C &
28 &, CVAP => con_CVAP, CLIQ => con_CLIQ &
29 &, EPS => con_eps, EPSM1 => con_epsm1 &
30 &, ROVCP => con_rocp, RD => con_rd
31 !-------------------------------------------------------------------
32 IMPLICIT NONE
33 !-------------------------------------------------------------------
34 !-- U3D 3D u-velocity interpolated to theta points (m/s)
35 !-- V3D 3D v-velocity interpolated to theta points (m/s)
36 !-- TH3D 3D potential temperature (K)
37 !-- T3D temperature (K)
38 !-- QV3D 3D water vapor mixing ratio (Kg/Kg)
39 !-- QC3D 3D cloud mixing ratio (Kg/Kg)
40 !-- QI3D 3D ice mixing ratio (Kg/Kg)
41 !-- P8w 3D pressure at full levels (Pa)
42 !-- Pcps 3D pressure (Pa)
43 !-- PI3D 3D exner function (dimensionless)
44 !-- rr3D 3D dry air density (kg/m^3)
45 !-- RUBLTEN U tendency due to
46 ! PBL parameterization (m/s^2)
47 !-- RVBLTEN V tendency due to
48 ! PBL parameterization (m/s^2)
49 !-- RTHBLTEN Theta tendency due to
50 ! PBL parameterization (K/s)
51 !-- RQVBLTEN Qv tendency due to
52 ! PBL parameterization (kg/kg/s)
53 !-- RQCBLTEN Qc tendency due to
54 ! PBL parameterization (kg/kg/s)
55 !-- RQIBLTEN Qi tendency due to
56 ! PBL parameterization (kg/kg/s)
57 !
58 !-- MOMMIX MOMENTUM MIXING COEFFICIENT (can be set in the namelist)
59 !-- RUCUTEN U tendency due to Cumulus Momentum Mixing (gopal's doing for SAS)
60 !-- RVCUTEN V tendency due to Cumulus Momentum Mixing (gopal's doing for SAS)
61 !
62 !-- CP heat capacity at constant pressure for dry air (J/kg/K)
63 !-- GRAV acceleration due to gravity (m/s^2)
64 !-- ROVCP R/CP
65 !-- RD gas constant for dry air (J/kg/K)
66 !-- ROVG R/G
67 !-- P_QI species index for cloud ice
68 !-- dz8w dz between full levels (m)
69 !-- z height above sea level (m)
70 !-- PSFC pressure at the surface (Pa)
71 !-- UST u* in similarity theory (m/s)
72 !-- PBL PBL height (m)
73 !-- PSIM similarity stability function for momentum
74 !-- PSIH similarity stability function for heat
75 !-- HFX upward heat flux at the surface (W/m^2)
76 !-- QFX upward moisture flux at the surface (kg/m^2/s)
77 !-- TSK surface temperature (K)
78 !-- GZ1OZ0 log(z/z0) where z0 is roughness length
79 !-- WSPD wind speed at lowest model level (m/s)
80 !-- BR bulk Richardson number in surface layer
81 !-- DT time step (s)
82 !-- rvovrd R_v divided by R_d (dimensionless)
83 !-- EP1 constant for virtual temperature (R_v/R_d - 1) (dimensionless)
84 !-- KARMAN Von Karman constant
85 !-- ids start index for i in domain
86 !-- ide end index for i in domain
87 !-- jds start index for j in domain
88 !-- jde end index for j in domain
89 !-- kds start index for k in domain
90 !-- kde end index for k in domain
91 !-- ims start index for i in memory
92 !-- ime end index for i in memory
93 !-- jms start index for j in memory
94 !-- jme end index for j in memory
95 !-- kms start index for k in memory
96 !-- kme end index for k in memory
97 !-- its start index for i in tile
98 !-- ite end index for i in tile
99 !-- jts start index for j in tile
100 !-- jte end index for j in tile
101 !-- kts start index for k in tile
102 !-- kte end index for k in tile
103 !-------------------------------------------------------------------
105 INTEGER :: ICLDCK
107 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
108 ims,ime, jms,jme, kms,kme, &
109 its,ite, jts,jte, kts,kte, &
110 ITIMESTEP, & !NSTD
111 P_FIRST_SCALAR, &
112 P_QC, &
113 P_QI, &
114 STEPCU
116 REAL, INTENT(IN) :: &
117 DT
120 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: &
121 RQCCUTEN, &
122 RQICUTEN, &
123 RQVCUTEN, &
124 RTHCUTEN
125 REAL, DIMENSION(ims:ime, jms:jme, kms:kme), INTENT(INOUT) :: &
126 RUCUTEN, & ! gopal's doing for SAS
127 RVCUTEN ! gopal's doing for SAS
128 #if (NMM_CORE == 1)
129 REAL, INTENT(IN) :: MOMMIX
130 REAL, DIMENSION( ims:ime , jms:jme ), &
131 INTENT(IN) :: STORE_RAND
133 #endif
135 REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: &
136 XLAND
138 REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: &
139 RAINCV, PRATEC
141 REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: &
142 HBOT, &
143 HTOP
145 LOGICAL, DIMENSION(IMS:IME,JMS:JME), INTENT(INOUT) :: &
146 CU_ACT_FLAG
149 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: &
150 DZ8W, &
151 P8w, &
152 Pcps, &
153 PI3D, &
154 QC3D, &
155 QI3D, &
156 QV3D, &
157 RHO3D, &
158 T3D, &
159 U3D, &
160 V3D, &
161 W
163 ! Adaptive time-step variables
164 REAL, INTENT(IN ) :: CUDT
165 REAL, INTENT(IN ) :: CURR_SECS
166 LOGICAL,INTENT(IN ) :: ADAPT_STEP_FLAG
168 !--------------------------- LOCAL VARS ------------------------------
170 REAL, DIMENSION(ims:ime, jms:jme) :: &
171 PSFC
174 REAL (kind=kind_phys) :: &
175 DELT, &
176 DPSHC, &
177 RDELT, &
178 RSEED
180 REAL (kind=kind_phys), DIMENSION(its:ite) :: &
181 CLDWRK, &
182 PS, &
183 RCS, &
184 RN, &
185 SLIMSK, &
186 XKT2
188 REAL (kind=kind_phys), DIMENSION(its:ite, kts:kte+1) :: &
189 PRSI
191 REAL (kind=kind_phys), DIMENSION(its:ite, kts:kte) :: &
192 DEL, &
193 DOT, &
194 PHIL, &
195 PRSL, &
196 PRSLK, &
197 Q1, &
198 T1, &
199 U1, &
200 V1, &
201 ZI, &
202 ZL
204 REAL (kind=kind_phys), DIMENSION(its:ite, kts:kte, 2) :: &
205 QL
207 INTEGER, DIMENSION(its:ite) :: &
208 KBOT, &
209 KTOP, &
210 KUO
212 INTEGER :: &
213 I, &
214 IGPVS, &
215 IM, &
216 J, &
217 JCAP, &
218 K, &
219 KM, &
220 KP, &
221 KX, &
222 NCLOUD
224 LOGICAL :: run_param
226 DATA IGPVS/0/
228 !-----------------------------------------------------------------------
229 !
230 !*** CHECK TO SEE IF THIS IS A CONVECTION TIMESTEP
231 !
232 if (adapt_step_flag) then
233 if ( (ITIMESTEP .eq. 0) .or. (cudt .eq. 0) .or. &
234 ( CURR_SECS + dt >= ( int( CURR_SECS / ( cudt * 60 ) ) + 1 ) * cudt * 60 ) ) then
235 run_param = .TRUE.
236 else
237 run_param = .FALSE.
238 endif
240 else
241 if (MOD(ITIMESTEP,STEPCU) .EQ. 0 .or. ITIMESTEP .eq. 0) then
242 run_param = .TRUE.
243 else
244 run_param = .FALSE.
245 endif
246 endif
248 !-----------------------------------------------------------------------
251 IF(run_param) THEN
253 DO J=JTS,JTE
254 DO I=ITS,ITE
255 CU_ACT_FLAG(I,J)=.TRUE.
256 ENDDO
257 ENDDO
259 IM=ITE-ITS+1
260 KX=KTE-KTS+1
261 JCAP=126
262 DPSHC=30_kind_phys
263 DELT=DT*STEPCU
264 RDELT=1./DELT
265 NCLOUD=1
268 DO J=jms,jme
269 DO I=ims,ime
270 PSFC(i,j)=P8w(i,kms,j)
271 ENDDO
272 ENDDO
274 if(igpvs.eq.0) CALL GFUNCPHYS
275 igpvs=1
277 !------------- J LOOP (OUTER) --------------------------------------------------
279 DO J=jts,jte
281 ! --------------- compute zi and zl -----------------------------------------
282 DO i=its,ite
283 ZI(I,KTS)=0.0
284 ENDDO
286 DO k=kts+1,kte
287 KM=K-1
288 DO i=its,ite
289 ZI(I,K)=ZI(I,KM)+dz8w(i,km,j)
290 ENDDO
291 ENDDO
293 DO k=kts+1,kte
294 KM=K-1
295 DO i=its,ite
296 ZL(I,KM)=(ZI(I,K)+ZI(I,KM))*0.5
297 ENDDO
298 ENDDO
300 DO i=its,ite
301 ZL(I,KTE)=2.*ZI(I,KTE)-ZL(I,KTE-1)
302 ENDDO
304 ! --------------- end compute zi and zl -------------------------------------
306 ! Based on some important findings from Morris Bender, XKT2 was defined in
307 ! terms of random number instead of random number based cloud tops
308 ! Also, these random numbers are stored and are changed only once in
309 ! approximately 5 minutes interval now. This is gopal's doing for HWRF.
311 ! call random_number(XKT2)
313 #if (EM_CORE == 1)
314 ! XKT2 was defined in terms of random number instead of random number based cloud tops
315 ! ZCX
316 call init_random_seed()
317 call random_number(XKT2)
318 #ifdef REGTEST
319 ! for regtest only
320 xkt2 = 0.1
321 #endif
322 #endif
323 !
324 #if (NMM_CORE == 1)
325 DO i=its,ite
326 XKT2(i) = STORE_RAND(i,j)
327 ENDDO
328 #endif
330 DO i=its,ite
331 PS(i)=PSFC(i,j)*.001
332 RCS(i)=1.
333 SLIMSK(i)=ABS(XLAND(i,j)-2.)
334 ENDDO
336 DO i=its,ite
337 PRSI(i,kts)=PS(i)
338 ENDDO
340 DO k=kts,kte
341 kp=k+1
342 DO i=its,ite
343 PRSL(I,K)=Pcps(i,k,j)*.001
344 PHIL(I,K)=ZL(I,K)*GRAV
345 DOT(i,k)=-5.0E-4*GRAV*rho3d(i,k,j)*(w(i,k,j)+w(i,kp,j))
346 ENDDO
347 ENDDO
349 DO k=kts,kte
350 DO i=its,ite
351 DEL(i,k)=PRSL(i,k)*GRAV/RD*dz8w(i,k,j)/T3D(i,k,j)
352 U1(i,k)=U3D(i,k,j)
353 V1(i,k)=V3D(i,k,j)
354 Q1(i,k)=QV3D(i,k,j)/(1.+QV3D(i,k,j))
355 T1(i,k)=T3D(i,k,j)
356 QL(i,k,1)=QI3D(i,k,j)/(1.+QI3D(i,k,j))
357 QL(i,k,2)=QC3D(i,k,j)/(1.+QC3D(i,k,j))
358 PRSLK(I,K)=(PRSL(i,k)*.01)**ROVCP
359 ENDDO
360 ENDDO
362 DO k=kts+1,kte+1
363 km=k-1
364 DO i=its,ite
365 PRSI(i,k)=PRSI(i,km)-del(i,km)
366 ENDDO
367 ENDDO
370 CALL SASCNV(IM,IM,KX,JCAP,DELT,DEL,PRSL,PS,PHIL, &
371 QL,Q1,T1,U1,V1,RCS,CLDWRK,RN,KBOT, &
372 KTOP,KUO,SLIMSK,DOT,XKT2,NCLOUD)
374 !!! make more like GFDL ... eliminate shallow convection.....
375 !!! CALL SHALCV(IM,IM,KX,DELT,DEL,PRSI,PRSL,PRSLK,KUO,Q1,T1,DPSHC)
376 #if (EM_CORE == 1)
377 CALL SHALCV(IM,IM,KX,DELT,DEL,PRSI,PRSL,PRSLK,KUO,Q1,T1,DPSHC)
378 #endif
380 DO I=ITS,ITE
381 RAINCV(I,J)=RN(I)*1000./STEPCU
382 PRATEC(I,J)=RN(I)*1000./(STEPCU * DT)
383 HBOT(I,J)=KBOT(I)
384 HTOP(I,J)=KTOP(I)
385 ENDDO
387 DO K=KTS,KTE
388 DO I=ITS,ITE
389 RTHCUTEN(I,K,J)=(T1(I,K)-T3D(I,K,J))/PI3D(I,K,J)*RDELT
390 RQVCUTEN(I,K,J)=(Q1(I,K)/(1.-q1(i,k))-QV3D(I,K,J))*RDELT
391 ENDDO
392 ENDDO
394 !===============================================================================
395 ! ADD MOMENTUM MIXING TERM AS TENDENCIES. This is gopal's doing for SAS
396 ! MOMMIX is the reduction factor set to 0.7 by default. Because NMM has
397 ! divergence damping term, a reducion factor for cumulum mixing may be
398 ! required otherwise storms were too weak.
399 !===============================================================================
400 !
401 #if (NMM_CORE == 1)
402 DO K=KTS,KTE
403 DO I=ITS,ITE
404 RUCUTEN(I,J,K)=MOMMIX*(U1(I,K)-U3D(I,K,J))*RDELT
405 RVCUTEN(I,J,K)=MOMMIX*(V1(I,K)-V3D(I,K,J))*RDELT
406 ENDDO
407 ENDDO
408 #endif
411 IF(P_QC .ge. P_FIRST_SCALAR)THEN
412 DO K=KTS,KTE
413 DO I=ITS,ITE
414 RQCCUTEN(I,K,J)=(QL(I,K,2)/(1.-ql(i,k,2))-QC3D(I,K,J))*RDELT
415 ENDDO
416 ENDDO
417 ENDIF
419 IF(P_QI .ge. P_FIRST_SCALAR)THEN
420 DO K=KTS,KTE
421 DO I=ITS,ITE
422 RQICUTEN(I,K,J)=(QL(I,K,1)/(1.-ql(i,k,1))-QI3D(I,K,J))*RDELT
423 ENDDO
424 ENDDO
425 ENDIF
427 ENDDO ! Outer most J loop
429 ENDIF
431 END SUBROUTINE CU_SAS
433 !====================================================================
434 SUBROUTINE sasinit(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQICUTEN, &
435 RUCUTEN,RVCUTEN, & ! gopal's doing for SAS
436 RESTART,P_QC,P_QI,P_FIRST_SCALAR, &
437 allowed_to_read, &
438 ids, ide, jds, jde, kds, kde, &
439 ims, ime, jms, jme, kms, kme, &
440 its, ite, jts, jte, kts, kte )
441 !--------------------------------------------------------------------
442 IMPLICIT NONE
443 !--------------------------------------------------------------------
444 LOGICAL , INTENT(IN) :: allowed_to_read,restart
445 INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
446 ims, ime, jms, jme, kms, kme, &
447 its, ite, jts, jte, kts, kte
448 INTEGER , INTENT(IN) :: P_FIRST_SCALAR, P_QI, P_QC
450 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
451 RTHCUTEN, &
452 RQVCUTEN, &
453 RQCCUTEN, &
454 RQICUTEN
455 REAL, DIMENSION( ims:ime , jms:jme , kms:kme ) , INTENT(OUT) :: &
456 RUCUTEN, & ! gopal's doing for SAS
457 RVCUTEN
459 INTEGER :: i, j, k, itf, jtf, ktf
461 jtf=min0(jte,jde-1)
462 ktf=min0(kte,kde-1)
463 itf=min0(ite,ide-1)
465 #ifdef HWRF
466 !zhang's doing
467 IF(.not.restart .or. .not.allowed_to_read)THEN
468 !end of zhang's doing
469 #else
470 IF(.not.restart)THEN
471 #endif
472 DO j=jts,jtf
473 DO k=kts,ktf
474 DO i=its,itf
475 RTHCUTEN(i,k,j)=0.
476 RQVCUTEN(i,k,j)=0.
477 RUCUTEN(i,j,k)=0. ! gopal's doing for SAS
478 RVCUTEN(i,j,k)=0. ! gopal's doing for SAS
479 ENDDO
480 ENDDO
481 ENDDO
483 IF (P_QC .ge. P_FIRST_SCALAR) THEN
484 DO j=jts,jtf
485 DO k=kts,ktf
486 DO i=its,itf
487 RQCCUTEN(i,k,j)=0.
488 ENDDO
489 ENDDO
490 ENDDO
491 ENDIF
493 IF (P_QI .ge. P_FIRST_SCALAR) THEN
494 DO j=jts,jtf
495 DO k=kts,ktf
496 DO i=its,itf
497 RQICUTEN(i,k,j)=0.
498 ENDDO
499 ENDDO
500 ENDDO
501 ENDIF
502 ENDIF
504 END SUBROUTINE sasinit
506 ! ------------------------------------------------------------------------
508 SUBROUTINE SASCNV(IM,IX,KM,JCAP,DELT,DEL,PRSL,PS,PHIL,QL, &
509 ! SUBROUTINE SASCNV(IM,IX,KM,JCAP,DLT,DEL,PRSL,PHIL,QL, &
510 & Q1,T1,U1,V1,RCS,CLDWRK,RN,KBOT,KTOP,KUO,SLIMSK, &
511 & DOT,XKT2,ncloud)
512 ! for cloud water version
513 ! parameter(ncloud=0)
514 ! SUBROUTINE SASCNV(KM,JCAP,DELT,DEL,SL,SLK,PS,QL,
515 ! & Q1,T1,U1,V1,RCS,CLDWRK,RN,KBOT,KTOP,KUO,SLIMSK,
516 ! & DOT,xkt2,ncloud)
517 !
518 USE MODULE_GFS_MACHINE , ONLY : kind_phys
519 USE MODULE_GFS_FUNCPHYS ,ONLY : fpvs
520 USE MODULE_GFS_PHYSCONS, grav => con_g, CP => con_CP, HVAP => con_HVAP &
521 &, RV => con_RV, FV => con_fvirt, T0C => con_T0C &
522 &, CVAP => con_CVAP, CLIQ => con_CLIQ &
523 &, EPS => con_eps, EPSM1 => con_epsm1
525 implicit none
526 !
527 ! include 'constant.h'
528 !
529 integer IM, IX, KM, JCAP, ncloud, &
530 & KBOT(IM), KTOP(IM), KUO(IM), J
531 real(kind=kind_phys) DELT
532 real(kind=kind_phys) PS(IM), DEL(IX,KM), PRSL(IX,KM), &
533 ! real(kind=kind_phys) DEL(IX,KM), PRSL(IX,KM),
534 & QL(IX,KM,2), Q1(IX,KM), T1(IX,KM), &
535 & U1(IX,KM), V1(IX,KM), RCS(IM), &
536 & CLDWRK(IM), RN(IM), SLIMSK(IM), &
537 & DOT(IX,KM), XKT2(IM), PHIL(IX,KM)
538 !
539 integer I, INDX, jmn, k, knumb, latd, lond, km1
540 !
541 real(kind=kind_phys) adw, alpha, alphal, alphas, &
542 & aup, beta, betal, betas, &
543 & c0, cpoel, dellat, delta, &
544 & desdt, deta, detad, dg, &
545 & dh, dhh, dlnsig, dp, &
546 & dq, dqsdp, dqsdt, dt, &
547 & dt2, dtmax, dtmin, dv1, &
548 & dv1q, dv2, dv2q, dv1u, &
549 & dv1v, dv2u, dv2v, dv3u, &
550 & dv3v, dv3, dv3q, dvq1, &
551 & dz, dz1, e1, edtmax, &
552 & edtmaxl, edtmaxs, el2orc, elocp, &
553 & es, etah, &
554 & evef, evfact, evfactl, fact1, &
555 & fact2, factor, fjcap, fkm, &
556 & fuv, g, gamma, onemf, &
557 & onemfu, pdetrn, pdpdwn, pprime, &
558 & qc, qlk, qrch, qs, &
559 & rain, rfact, shear, tem1, &
560 & tem2, terr, val, val1, &
561 & val2, w1, w1l, w1s, &
562 & w2, w2l, w2s, w3, &
563 & w3l, w3s, w4, w4l, &
564 & w4s, xdby, xpw, xpwd, &
565 & xqc, xqrch, xlambu, mbdt, &
566 & tem
567 !
568 !
569 integer JMIN(IM), KB(IM), KBCON(IM), KBDTR(IM), &
570 & KT2(IM), KTCON(IM), LMIN(IM), &
571 & kbm(IM), kbmax(IM), kmax(IM)
572 !
573 real(kind=kind_phys) AA1(IM), ACRT(IM), ACRTFCT(IM), &
574 & DELHBAR(IM), DELQ(IM), DELQ2(IM), &
575 & DELQBAR(IM), DELQEV(IM), DELTBAR(IM), &
576 & DELTV(IM), DTCONV(IM), EDT(IM), &
577 & EDTO(IM), EDTX(IM), FLD(IM), &
578 & HCDO(IM), HKBO(IM), HMAX(IM), &
579 & HMIN(IM), HSBAR(IM), UCDO(IM), &
580 & UKBO(IM), VCDO(IM), VKBO(IM), &
581 & PBCDIF(IM), PDOT(IM), PO(IM,KM), &
582 & PWAVO(IM), PWEVO(IM), &
583 ! & PSFC(IM), PWAVO(IM), PWEVO(IM), &
584 & QCDO(IM), QCOND(IM), QEVAP(IM), &
585 & QKBO(IM), RNTOT(IM), VSHEAR(IM), &
586 & XAA0(IM), XHCD(IM), XHKB(IM), &
587 & XK(IM), XLAMB(IM), XLAMD(IM), &
588 & XMB(IM), XMBMAX(IM), XPWAV(IM), &
589 & XPWEV(IM), XQCD(IM), XQKB(IM)
590 !
591 ! PHYSICAL PARAMETERS
592 PARAMETER(G=grav)
593 PARAMETER(CPOEL=CP/HVAP,ELOCP=HVAP/CP, &
594 & EL2ORC=HVAP*HVAP/(RV*CP))
595 PARAMETER(TERR=0.,C0=.002,DELTA=fv)
596 PARAMETER(FACT1=(CVAP-CLIQ)/RV,FACT2=HVAP/RV-FACT1*T0C)
597 ! LOCAL VARIABLES AND ARRAYS
598 real(kind=kind_phys) PFLD(IM,KM), TO(IM,KM), QO(IM,KM), &
599 & UO(IM,KM), VO(IM,KM), QESO(IM,KM)
600 ! cloud water
601 real(kind=kind_phys) QLKO_KTCON(IM), DELLAL(IM), TVO(IM,KM), &
602 & DBYO(IM,KM), ZO(IM,KM), SUMZ(IM,KM), &
603 & SUMH(IM,KM), HEO(IM,KM), HESO(IM,KM), &
604 & QRCD(IM,KM), DELLAH(IM,KM), DELLAQ(IM,KM),&
605 & DELLAU(IM,KM), DELLAV(IM,KM), HCKO(IM,KM), &
606 & UCKO(IM,KM), VCKO(IM,KM), QCKO(IM,KM), &
607 & ETA(IM,KM), ETAU(IM,KM), ETAD(IM,KM), &
608 & QRCDO(IM,KM), PWO(IM,KM), PWDO(IM,KM), &
609 & RHBAR(IM), TX1(IM)
610 !
611 LOGICAL TOTFLG, CNVFLG(IM), DWNFLG(IM), DWNFLG2(IM), FLG(IM)
612 !
613 real(kind=kind_phys) PCRIT(15), ACRITT(15), ACRIT(15)
614 ! SAVE PCRIT, ACRITT
615 DATA PCRIT/850.,800.,750.,700.,650.,600.,550.,500.,450.,400., &
616 & 350.,300.,250.,200.,150./
617 DATA ACRITT/.0633,.0445,.0553,.0664,.075,.1082,.1521,.2216, &
618 & .3151,.3677,.41,.5255,.7663,1.1686,1.6851/
619 ! GDAS DERIVED ACRIT
620 ! DATA ACRITT/.203,.515,.521,.566,.625,.665,.659,.688, &
621 ! & .743,.813,.886,.947,1.138,1.377,1.896/
622 !
623 real(kind=kind_phys) TF, TCR, TCRF, RZERO, RONE
624 parameter (TF=233.16, TCR=263.16, TCRF=1.0/(TCR-TF))
625 parameter (RZERO=0.0,RONE=1.0)
626 !-----------------------------------------------------------------------
627 !
628 km1 = km - 1
629 ! INITIALIZE ARRAYS
630 !
631 DO I=1,IM
632 RN(I)=0.
633 KBOT(I)=KM+1
634 KTOP(I)=0
635 KUO(I)=0
636 CNVFLG(I) = .TRUE.
637 DTCONV(I) = 3600.
638 CLDWRK(I) = 0.
639 PDOT(I) = 0.
640 KT2(I) = 0
641 QLKO_KTCON(I) = 0.
642 DELLAL(I) = 0.
643 ENDDO
644 !!
645 DO K = 1, 15
646 ACRIT(K) = ACRITT(K) * (975. - PCRIT(K))
647 ENDDO
648 DT2 = DELT
649 !cmr dtmin = max(dt2,1200.)
650 val = 1200.
651 dtmin = max(dt2, val )
652 !cmr dtmax = max(dt2,3600.)
653 val = 3600.
654 dtmax = max(dt2, val )
655 ! MODEL TUNABLE PARAMETERS ARE ALL HERE
656 MBDT = 10.
657 EDTMAXl = .3
658 EDTMAXs = .3
659 ALPHAl = .5
660 ALPHAs = .5
661 BETAl = .15
662 betas = .15
663 BETAl = .05
664 betas = .05
665 ! change for hurricane model
666 BETAl = .5
667 betas = .5
668 ! EVEF = 0.07
669 evfact = 0.3
670 evfactl = 0.3
671 ! change for hurricane model
672 evfact = 0.6
673 evfactl = .6
674 #if ( EM_CORE == 1 )
675 ! HAWAII TEST - ZCX
676 ALPHAl = .5
677 ALPHAs = .75
678 BETAl = .05
679 betas = .05
680 evfact = 0.5
681 evfactl = 0.5
682 #endif
683 PDPDWN = 0.
684 PDETRN = 200.
685 xlambu = 1.e-4
686 fjcap = (float(jcap) / 126.) ** 2
687 !cmr fjcap = max(fjcap,1.)
688 val = 1.
689 fjcap = max(fjcap,val)
690 fkm = (float(km) / 28.) ** 2
691 !cmr fkm = max(fkm,1.)
692 fkm = max(fkm,val)
693 W1l = -8.E-3
694 W2l = -4.E-2
695 W3l = -5.E-3
696 W4l = -5.E-4
697 W1s = -2.E-4
698 W2s = -2.E-3
699 W3s = -1.E-3
700 W4s = -2.E-5
701 !CCCC IF(IM.EQ.384) THEN
702 LATD = 92
703 lond = 189
704 !CCCC ELSEIF(IM.EQ.768) THEN
705 !CCCC LATD = 80
706 !CCCC ELSE
707 !CCCC LATD = 0
708 !CCCC ENDIF
709 !
710 ! DEFINE TOP LAYER FOR SEARCH OF THE DOWNDRAFT ORIGINATING LAYER
711 ! AND THE MAXIMUM THETAE FOR UPDRAFT
712 !
713 DO I=1,IM
714 KBMAX(I) = KM
715 KBM(I) = KM
716 KMAX(I) = KM
717 TX1(I) = 1.0 / PS(I)
718 ENDDO
719 !
720 DO K = 1, KM
721 DO I=1,IM
722 IF (prSL(I,K)*tx1(I) .GT. 0.45) KBMAX(I) = K + 1
723 IF (prSL(I,K)*tx1(I) .GT. 0.70) KBM(I) = K + 1
724 IF (prSL(I,K)*tx1(I) .GT. 0.04) KMAX(I) = MIN(KM,K + 1)
725 ENDDO
726 ENDDO
727 DO I=1,IM
728 KBMAX(I) = MIN(KBMAX(I),KMAX(I))
729 KBM(I) = MIN(KBM(I),KMAX(I))
730 ENDDO
731 !
732 ! CONVERT SURFACE PRESSURE TO MB FROM CB
733 !
734 !!
735 DO K = 1, KM
736 DO I=1,IM
737 if (K .le. kmax(i)) then
738 PFLD(I,k) = PRSL(I,K) * 10.0
739 PWO(I,k) = 0.
740 PWDO(I,k) = 0.
741 TO(I,k) = T1(I,k)
742 QO(I,k) = Q1(I,k)
743 UO(I,k) = U1(I,k)
744 VO(I,k) = V1(I,k)
745 DBYO(I,k) = 0.
746 SUMZ(I,k) = 0.
747 SUMH(I,k) = 0.
748 endif
749 ENDDO
750 ENDDO
752 !
753 ! COLUMN VARIABLES
754 ! P IS PRESSURE OF THE LAYER (MB)
755 ! T IS TEMPERATURE AT T-DT (K)..TN
756 ! Q IS MIXING RATIO AT T-DT (KG/KG)..QN
757 ! TO IS TEMPERATURE AT T+DT (K)... THIS IS AFTER ADVECTION AND TURBULAN
758 ! QO IS MIXING RATIO AT T+DT (KG/KG)..Q1
759 !
760 DO K = 1, KM
761 DO I=1,IM
762 if (k .le. kmax(i)) then
763 !jfe QESO(I,k) = 10. * FPVS(T1(I,k))
764 !
765 QESO(I,k) = 0.01 * fpvs(T1(I,K)) ! fpvs is in Pa
766 !
767 QESO(I,k) = EPS * QESO(I,k) / (PFLD(I,k) + EPSM1*QESO(I,k))
768 !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
769 val1 = 1.E-8
770 QESO(I,k) = MAX(QESO(I,k), val1)
771 !cmr QO(I,k) = max(QO(I,k),1.e-10)
772 val2 = 1.e-10
773 QO(I,k) = max(QO(I,k), val2 )
774 ! QO(I,k) = MIN(QO(I,k),QESO(I,k))
775 TVO(I,k) = TO(I,k) + DELTA * TO(I,k) * QO(I,k)
776 endif
777 ENDDO
778 ENDDO
780 !
781 ! HYDROSTATIC HEIGHT ASSUME ZERO TERR
782 !
783 DO K = 1, KM
784 DO I=1,IM
785 ZO(I,k) = PHIL(I,k) / G
786 ENDDO
787 ENDDO
788 ! COMPUTE MOIST STATIC ENERGY
789 DO K = 1, KM
790 DO I=1,IM
791 if (K .le. kmax(i)) then
792 ! tem = G * ZO(I,k) + CP * TO(I,k)
793 tem = PHIL(I,k) + CP * TO(I,k)
794 HEO(I,k) = tem + HVAP * QO(I,k)
795 HESO(I,k) = tem + HVAP * QESO(I,k)
796 ! HEO(I,k) = MIN(HEO(I,k),HESO(I,k))
797 endif
798 ENDDO
799 ENDDO
800 !
801 ! DETERMINE LEVEL WITH LARGEST MOIST STATIC ENERGY
802 ! THIS IS THE LEVEL WHERE UPDRAFT STARTS
803 !
804 DO I=1,IM
805 HMAX(I) = HEO(I,1)
806 KB(I) = 1
807 ENDDO
808 !!
809 DO K = 2, KM
810 DO I=1,IM
811 if (k .le. kbm(i)) then
812 IF(HEO(I,k).GT.HMAX(I).AND.CNVFLG(I)) THEN
813 KB(I) = K
814 HMAX(I) = HEO(I,k)
815 ENDIF
816 endif
817 ENDDO
818 ENDDO
819 ! DO K = 1, KMAX - 1
820 ! TOL(k) = .5 * (TO(I,k) + TO(I,k+1))
821 ! QOL(k) = .5 * (QO(I,k) + QO(I,k+1))
822 ! QESOL(I,k) = .5 * (QESO(I,k) + QESO(I,k+1))
823 ! HEOL(I,k) = .5 * (HEO(I,k) + HEO(I,k+1))
824 ! HESOL(I,k) = .5 * (HESO(I,k) + HESO(I,k+1))
825 ! ENDDO
826 DO K = 1, KM1
827 DO I=1,IM
828 if (k .le. kmax(i)-1) then
829 DZ = .5 * (ZO(I,k+1) - ZO(I,k))
830 DP = .5 * (PFLD(I,k+1) - PFLD(I,k))
831 !jfe ES = 10. * FPVS(TO(I,k+1))
832 !
833 ES = 0.01 * fpvs(TO(I,K+1)) ! fpvs is in Pa
834 !
835 PPRIME = PFLD(I,k+1) + EPSM1 * ES
836 QS = EPS * ES / PPRIME
837 DQSDP = - QS / PPRIME
838 DESDT = ES * (FACT1 / TO(I,k+1) + FACT2 / (TO(I,k+1)**2))
839 DQSDT = QS * PFLD(I,k+1) * DESDT / (ES * PPRIME)
840 GAMMA = EL2ORC * QESO(I,k+1) / (TO(I,k+1)**2)
841 DT = (G * DZ + HVAP * DQSDP * DP) / (CP * (1. + GAMMA))
842 DQ = DQSDT * DT + DQSDP * DP
843 TO(I,k) = TO(I,k+1) + DT
844 QO(I,k) = QO(I,k+1) + DQ
845 PO(I,k) = .5 * (PFLD(I,k) + PFLD(I,k+1))
846 endif
847 ENDDO
848 ENDDO
849 !
850 DO K = 1, KM1
851 DO I=1,IM
852 if (k .le. kmax(I)-1) then
853 !jfe QESO(I,k) = 10. * FPVS(TO(I,k))
854 !
855 QESO(I,k) = 0.01 * fpvs(TO(I,K)) ! fpvs is in Pa
856 !
857 QESO(I,k) = EPS * QESO(I,k) / (PO(I,k) + EPSM1*QESO(I,k))
858 !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
859 val1 = 1.E-8
860 QESO(I,k) = MAX(QESO(I,k), val1)
861 !cmr QO(I,k) = max(QO(I,k),1.e-10)
862 val2 = 1.e-10
863 QO(I,k) = max(QO(I,k), val2 )
864 ! QO(I,k) = MIN(QO(I,k),QESO(I,k))
865 HEO(I,k) = .5 * G * (ZO(I,k) + ZO(I,k+1)) + &
866 & CP * TO(I,k) + HVAP * QO(I,k)
867 HESO(I,k) = .5 * G * (ZO(I,k) + ZO(I,k+1)) + &
868 & CP * TO(I,k) + HVAP * QESO(I,k)
869 UO(I,k) = .5 * (UO(I,k) + UO(I,k+1))
870 VO(I,k) = .5 * (VO(I,k) + VO(I,k+1))
871 endif
872 ENDDO
873 ENDDO
874 ! k = kmax
875 ! HEO(I,k) = HEO(I,k)
876 ! hesol(k) = HESO(I,k)
877 ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
878 ! PRINT *, ' HEO ='
879 ! PRINT 6001, (HEO(I,K),K=1,KMAX)
880 ! PRINT *, ' HESO ='
881 ! PRINT 6001, (HESO(I,K),K=1,KMAX)
882 ! PRINT *, ' TO ='
883 ! PRINT 6002, (TO(I,K)-273.16,K=1,KMAX)
884 ! PRINT *, ' QO ='
885 ! PRINT 6003, (QO(I,K),K=1,KMAX)
886 ! PRINT *, ' QSO ='
887 ! PRINT 6003, (QESO(I,K),K=1,KMAX)
888 ! ENDIF
889 !
890 ! LOOK FOR CONVECTIVE CLOUD BASE AS THE LEVEL OF FREE CONVECTION
891 !
892 DO I=1,IM
893 IF(CNVFLG(I)) THEN
894 INDX = KB(I)
895 HKBO(I) = HEO(I,INDX)
896 QKBO(I) = QO(I,INDX)
897 UKBO(I) = UO(I,INDX)
898 VKBO(I) = VO(I,INDX)
899 ENDIF
900 FLG(I) = CNVFLG(I)
901 KBCON(I) = KMAX(I)
902 ENDDO
903 !!
904 DO K = 1, KM
905 DO I=1,IM
906 if (k .le. kbmax(i)) then
907 IF(FLG(I).AND.K.GT.KB(I)) THEN
908 HSBAR(I) = HESO(I,k)
909 IF(HKBO(I).GT.HSBAR(I)) THEN
910 FLG(I) = .FALSE.
911 KBCON(I) = K
912 ENDIF
913 ENDIF
914 endif
915 ENDDO
916 ENDDO
917 DO I=1,IM
918 IF(CNVFLG(I)) THEN
919 PBCDIF(I) = -PFLD(I,KBCON(I)) + PFLD(I,KB(I))
920 PDOT(I) = 10.* DOT(I,KBCON(I))
921 IF(PBCDIF(I).GT.150.) CNVFLG(I) = .FALSE.
922 IF(KBCON(I).EQ.KMAX(I)) CNVFLG(I) = .FALSE.
923 ENDIF
924 ENDDO
925 !!
926 TOTFLG = .TRUE.
927 DO I=1,IM
928 TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
929 ENDDO
930 IF(TOTFLG) RETURN
931 ! FOUND LFC, CAN DEFINE REST OF VARIABLES
932 6001 FORMAT(2X,-2P10F12.2)
933 6002 FORMAT(2X,10F12.2)
934 6003 FORMAT(2X,3P10F12.2)
936 !
937 ! DETERMINE ENTRAINMENT RATE BETWEEN KB AND KBCON
938 !
939 DO I = 1, IM
940 alpha = alphas
941 if(SLIMSK(I).eq.1.) alpha = alphal
942 IF(CNVFLG(I)) THEN
943 IF(KB(I).EQ.1) THEN
944 DZ = .5 * (ZO(I,KBCON(I)) + ZO(I,KBCON(I)-1)) - ZO(I,1)
945 ELSE
946 DZ = .5 * (ZO(I,KBCON(I)) + ZO(I,KBCON(I)-1)) &
947 & - .5 * (ZO(I,KB(I)) + ZO(I,KB(I)-1))
948 ENDIF
949 IF(KBCON(I).NE.KB(I)) THEN
950 !cmr XLAMB(I) = -ALOG(ALPHA) / DZ
951 XLAMB(I) = - LOG(ALPHA) / DZ
952 ELSE
953 XLAMB(I) = 0.
954 ENDIF
955 ENDIF
956 ENDDO
957 ! DETERMINE UPDRAFT MASS FLUX
958 DO K = 1, KM
959 DO I = 1, IM
960 if (k .le. kmax(i) .and. CNVFLG(I)) then
961 ETA(I,k) = 1.
962 ETAU(I,k) = 1.
963 ENDIF
964 ENDDO
965 ENDDO
966 DO K = KM1, 2, -1
967 DO I = 1, IM
968 if (k .le. kbmax(i)) then
969 IF(CNVFLG(I).AND.K.LT.KBCON(I).AND.K.GE.KB(I)) THEN
970 DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
971 ETA(I,k) = ETA(I,k+1) * EXP(-XLAMB(I) * DZ)
972 ETAU(I,k) = ETA(I,k)
973 ENDIF
974 endif
975 ENDDO
976 ENDDO
977 DO I = 1, IM
978 IF(CNVFLG(I).AND.KB(I).EQ.1.AND.KBCON(I).GT.1) THEN
979 DZ = .5 * (ZO(I,2) - ZO(I,1))
980 ETA(I,1) = ETA(I,2) * EXP(-XLAMB(I) * DZ)
981 ETAU(I,1) = ETA(I,1)
982 ENDIF
983 ENDDO
984 !
985 ! WORK UP UPDRAFT CLOUD PROPERTIES
986 !
987 DO I = 1, IM
988 IF(CNVFLG(I)) THEN
989 INDX = KB(I)
990 HCKO(I,INDX) = HKBO(I)
991 QCKO(I,INDX) = QKBO(I)
992 UCKO(I,INDX) = UKBO(I)
993 VCKO(I,INDX) = VKBO(I)
994 PWAVO(I) = 0.
995 ENDIF
996 ENDDO
997 !
998 ! CLOUD PROPERTY BELOW CLOUD BASE IS MODIFIED BY THE ENTRAINMENT PROCES
999 !
1000 DO K = 2, KM1
1001 DO I = 1, IM
1002 if (k .le. kmax(i)-1) then
1003 IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LE.KBCON(I)) THEN
1004 FACTOR = ETA(I,k-1) / ETA(I,k)
1005 ONEMF = 1. - FACTOR
1006 HCKO(I,k) = FACTOR * HCKO(I,k-1) + ONEMF * &
1007 & .5 * (HEO(I,k) + HEO(I,k+1))
1008 UCKO(I,k) = FACTOR * UCKO(I,k-1) + ONEMF * &
1009 & .5 * (UO(I,k) + UO(I,k+1))
1010 VCKO(I,k) = FACTOR * VCKO(I,k-1) + ONEMF * &
1011 & .5 * (VO(I,k) + VO(I,k+1))
1012 DBYO(I,k) = HCKO(I,k) - HESO(I,k)
1013 ENDIF
1014 IF(CNVFLG(I).AND.K.GT.KBCON(I)) THEN
1015 HCKO(I,k) = HCKO(I,k-1)
1016 UCKO(I,k) = UCKO(I,k-1)
1017 VCKO(I,k) = VCKO(I,k-1)
1018 DBYO(I,k) = HCKO(I,k) - HESO(I,k)
1019 ENDIF
1020 endif
1021 ENDDO
1022 ENDDO
1023 ! DETERMINE CLOUD TOP
1024 DO I = 1, IM
1025 FLG(I) = CNVFLG(I)
1026 KTCON(I) = 1
1027 ENDDO
1028 ! DO K = 2, KMAX
1029 ! KK = KMAX - K + 1
1030 ! IF(DBYO(I,kK).GE.0..AND.FLG(I).AND.KK.GT.KBCON(I)) THEN
1031 ! KTCON(I) = KK + 1
1032 ! FLG(I) = .FALSE.
1033 ! ENDIF
1034 ! ENDDO
1035 DO K = 2, KM
1036 DO I = 1, IM
1037 if (k .le. kmax(i)) then
1038 IF(DBYO(I,k).LT.0..AND.FLG(I).AND.K.GT.KBCON(I)) THEN
1039 KTCON(I) = K
1040 FLG(I) = .FALSE.
1041 ENDIF
1042 endif
1043 ENDDO
1044 ENDDO
1045 DO I = 1, IM
1046 IF(CNVFLG(I).AND.(PFLD(I,KBCON(I)) - PFLD(I,KTCON(I))).LT.150.) &
1047 & CNVFLG(I) = .FALSE.
1048 ENDDO
1049 TOTFLG = .TRUE.
1050 DO I = 1, IM
1051 TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
1052 ENDDO
1053 IF(TOTFLG) RETURN
1054 !
1055 ! SEARCH FOR DOWNDRAFT ORIGINATING LEVEL ABOVE THETA-E MINIMUM
1056 !
1057 DO I = 1, IM
1058 HMIN(I) = HEO(I,KBCON(I))
1059 LMIN(I) = KBMAX(I)
1060 JMIN(I) = KBMAX(I)
1061 ENDDO
1062 DO I = 1, IM
1063 DO K = KBCON(I), KBMAX(I)
1064 IF(HEO(I,k).LT.HMIN(I).AND.CNVFLG(I)) THEN
1065 LMIN(I) = K + 1
1066 HMIN(I) = HEO(I,k)
1067 ENDIF
1068 ENDDO
1069 ENDDO
1070 !
1071 ! Make sure that JMIN(I) is within the cloud
1072 !
1073 DO I = 1, IM
1074 IF(CNVFLG(I)) THEN
1075 JMIN(I) = MIN(LMIN(I),KTCON(I)-1)
1076 XMBMAX(I) = .1
1077 JMIN(I) = MAX(JMIN(I),KBCON(I)+1)
1078 ENDIF
1079 ENDDO
1080 !
1081 ! ENTRAINING CLOUD
1082 !
1083 do k = 2, km1
1084 DO I = 1, IM
1085 if (k .le. kmax(i)-1) then
1086 if(CNVFLG(I).and.k.gt.JMIN(I).and.k.le.KTCON(I)) THEN
1087 SUMZ(I,k) = SUMZ(I,k-1) + .5 * (ZO(I,k+1) - ZO(I,k-1))
1088 SUMH(I,k) = SUMH(I,k-1) + .5 * (ZO(I,k+1) - ZO(I,k-1)) &
1089 & * HEO(I,k)
1090 ENDIF
1091 endif
1092 enddo
1093 enddo
1094 !!
1095 DO I = 1, IM
1096 IF(CNVFLG(I)) THEN
1097 ! call random_number(XKT2)
1098 ! call srand(fhour)
1099 ! XKT2(I) = rand()
1100 KT2(I) = nint(XKT2(I)*float(KTCON(I)-JMIN(I))-.5)+JMIN(I)+1
1101 ! KT2(I) = nint(sqrt(XKT2(I))*float(KTCON(I)-JMIN(I))-.5) + JMIN(I) + 1
1102 ! KT2(I) = nint(ranf() *float(KTCON(I)-JMIN(I))-.5) + JMIN(I) + 1
1103 tem1 = (HCKO(I,JMIN(I)) - HESO(I,KT2(I)))
1104 tem2 = (SUMZ(I,KT2(I)) * HESO(I,KT2(I)) - SUMH(I,KT2(I)))
1105 if (abs(tem2) .gt. 0.000001) THEN
1106 XLAMB(I) = tem1 / tem2
1107 else
1108 CNVFLG(I) = .false.
1109 ENDIF
1110 ! XLAMB(I) = (HCKO(I,JMIN(I)) - HESO(I,KT2(I)))
1111 ! & / (SUMZ(I,KT2(I)) * HESO(I,KT2(I)) - SUMH(I,KT2(I)))
1112 XLAMB(I) = max(XLAMB(I),RZERO)
1113 XLAMB(I) = min(XLAMB(I),2.3/SUMZ(I,KT2(I)))
1114 ENDIF
1115 ENDDO
1116 !!
1117 DO I = 1, IM
1118 DWNFLG(I) = CNVFLG(I)
1119 DWNFLG2(I) = CNVFLG(I)
1120 IF(CNVFLG(I)) THEN
1121 if(KT2(I).ge.KTCON(I)) DWNFLG(I) = .false.
1122 if(XLAMB(I).le.1.e-30.or.HCKO(I,JMIN(I))-HESO(I,KT2(I)).le.1.e-30)&
1123 & DWNFLG(I) = .false.
1124 do k = JMIN(I), KT2(I)
1125 if(DWNFLG(I).and.HEO(I,k).gt.HESO(I,KT2(I))) DWNFLG(I)=.false.
1126 enddo
1127 ! IF(CNVFLG(I).AND.(PFLD(KBCON(I))-PFLD(KTCON(I))).GT.PDETRN)
1128 ! & DWNFLG(I)=.FALSE.
1129 IF(CNVFLG(I).AND.(PFLD(I,KBCON(I))-PFLD(I,KTCON(I))).LT.PDPDWN) &
1130 & DWNFLG2(I)=.FALSE.
1131 ENDIF
1132 ENDDO
1133 !!
1134 DO K = 2, KM1
1135 DO I = 1, IM
1136 if (k .le. kmax(i)-1) then
1137 IF(DWNFLG(I).AND.K.GT.JMIN(I).AND.K.LE.KT2(I)) THEN
1138 DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1139 ! ETA(I,k) = ETA(I,k-1) * EXP( XLAMB(I) * DZ)
1140 ! to simplify matter, we will take the linear approach here
1141 !
1142 ETA(I,k) = ETA(I,k-1) * (1. + XLAMB(I) * dz)
1143 ETAU(I,k) = ETAU(I,k-1) * (1. + (XLAMB(I)+xlambu) * dz)
1144 ENDIF
1145 endif
1146 ENDDO
1147 ENDDO
1148 !!
1149 DO K = 2, KM1
1150 DO I = 1, IM
1151 if (k .le. kmax(i)-1) then
1152 ! IF(.NOT.DWNFLG(I).AND.K.GT.JMIN(I).AND.K.LE.KT2(I)) THEN
1153 IF(.NOT.DWNFLG(I).AND.K.GT.JMIN(I).AND.K.LE.KTCON(I)) THEN
1154 DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1155 ETAU(I,k) = ETAU(I,k-1) * (1. + xlambu * dz)
1156 ENDIF
1157 endif
1158 ENDDO
1159 ENDDO
1160 ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
1161 ! PRINT *, ' LMIN(I), KT2(I)=', LMIN(I), KT2(I)
1162 ! PRINT *, ' KBOT, KTOP, JMIN(I) =', KBCON(I), KTCON(I), JMIN(I)
1163 ! ENDIF
1164 ! IF(LAT.EQ.LATD.AND.lon.eq.lond) THEN
1165 ! print *, ' xlamb =', xlamb
1166 ! print *, ' eta =', (eta(k),k=1,KT2(I))
1167 ! print *, ' ETAU =', (ETAU(I,k),k=1,KT2(I))
1168 ! print *, ' HCKO =', (HCKO(I,k),k=1,KT2(I))
1169 ! print *, ' SUMZ =', (SUMZ(I,k),k=1,KT2(I))
1170 ! print *, ' SUMH =', (SUMH(I,k),k=1,KT2(I))
1171 ! ENDIF
1172 DO I = 1, IM
1173 if(DWNFLG(I)) THEN
1174 KTCON(I) = KT2(I)
1175 ENDIF
1176 ENDDO
1177 !
1178 ! CLOUD PROPERTY ABOVE CLOUD Base IS MODIFIED BY THE DETRAINMENT PROCESS
1179 !
1180 DO K = 2, KM1
1181 DO I = 1, IM
1182 if (k .le. kmax(i)-1) then
1183 !jfe
1184 IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
1185 !jfe IF(K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
1186 FACTOR = ETA(I,k-1) / ETA(I,k)
1187 ONEMF = 1. - FACTOR
1188 fuv = ETAU(I,k-1) / ETAU(I,k)
1189 onemfu = 1. - fuv
1190 HCKO(I,k) = FACTOR * HCKO(I,k-1) + ONEMF * &
1191 & .5 * (HEO(I,k) + HEO(I,k+1))
1192 UCKO(I,k) = fuv * UCKO(I,k-1) + ONEMFu * &
1193 & .5 * (UO(I,k) + UO(I,k+1))
1194 VCKO(I,k) = fuv * VCKO(I,k-1) + ONEMFu * &
1195 & .5 * (VO(I,k) + VO(I,k+1))
1196 DBYO(I,k) = HCKO(I,k) - HESO(I,k)
1197 ENDIF
1198 endif
1199 ENDDO
1200 ENDDO
1201 ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
1202 ! PRINT *, ' UCKO=', (UCKO(I,k),k=KBCON(I)+1,KTCON(I))
1203 ! PRINT *, ' uenv=', (.5*(UO(I,k)+UO(I,k-1)),k=KBCON(I)+1,KTCON(I))
1204 ! ENDIF
1205 DO I = 1, IM
1206 if(CNVFLG(I).and.DWNFLG2(I).and.JMIN(I).le.KBCON(I)) &
1207 & THEN
1208 CNVFLG(I) = .false.
1209 DWNFLG(I) = .false.
1210 DWNFLG2(I) = .false.
1211 ENDIF
1212 ENDDO
1213 !!
1214 TOTFLG = .TRUE.
1215 DO I = 1, IM
1216 TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
1217 ENDDO
1218 IF(TOTFLG) RETURN
1219 !!
1220 !
1221 ! COMPUTE CLOUD MOISTURE PROPERTY AND PRECIPITATION
1222 !
1223 DO I = 1, IM
1224 AA1(I) = 0.
1225 RHBAR(I) = 0.
1226 ENDDO
1227 DO K = 1, KM
1228 DO I = 1, IM
1229 if (k .le. kmax(i)) then
1230 IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LT.KTCON(I)) THEN
1231 DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1232 DZ1 = (ZO(I,k) - ZO(I,k-1))
1233 GAMMA = EL2ORC * QESO(I,k) / (TO(I,k)**2)
1234 QRCH = QESO(I,k) &
1235 & + GAMMA * DBYO(I,k) / (HVAP * (1. + GAMMA))
1236 FACTOR = ETA(I,k-1) / ETA(I,k)
1237 ONEMF = 1. - FACTOR
1238 QCKO(I,k) = FACTOR * QCKO(I,k-1) + ONEMF * &
1239 & .5 * (QO(I,k) + QO(I,k+1))
1240 DQ = ETA(I,k) * QCKO(I,k) - ETA(I,k) * QRCH
1241 RHBAR(I) = RHBAR(I) + QO(I,k) / QESO(I,k)
1242 !
1243 ! BELOW LFC CHECK IF THERE IS EXCESS MOISTURE TO RELEASE LATENT HEAT
1244 !
1245 IF(DQ.GT.0.) THEN
1246 ETAH = .5 * (ETA(I,k) + ETA(I,k-1))
1247 QLK = DQ / (ETA(I,k) + ETAH * C0 * DZ)
1248 AA1(I) = AA1(I) - DZ1 * G * QLK
1249 QC = QLK + QRCH
1250 PWO(I,k) = ETAH * C0 * DZ * QLK
1251 QCKO(I,k) = QC
1252 PWAVO(I) = PWAVO(I) + PWO(I,k)
1253 ENDIF
1254 ENDIF
1255 endif
1256 ENDDO
1257 ENDDO
1258 DO I = 1, IM
1259 RHBAR(I) = RHBAR(I) / float(KTCON(I) - KB(I) - 1)
1260 ENDDO
1261 !
1262 ! this section is ready for cloud water
1263 !
1264 if(ncloud.gt.0) THEN
1265 !
1266 ! compute liquid and vapor separation at cloud top
1267 !
1268 DO I = 1, IM
1269 k = KTCON(I)
1270 IF(CNVFLG(I)) THEN
1271 GAMMA = EL2ORC * QESO(I,K) / (TO(I,K)**2)
1272 QRCH = QESO(I,K) &
1273 & + GAMMA * DBYO(I,K) / (HVAP * (1. + GAMMA))
1274 DQ = QCKO(I,K-1) - QRCH
1275 !
1276 ! CHECK IF THERE IS EXCESS MOISTURE TO RELEASE LATENT HEAT
1277 !
1278 IF(DQ.GT.0.) THEN
1279 QLKO_KTCON(I) = dq
1280 QCKO(I,K-1) = QRCH
1281 ENDIF
1282 ENDIF
1283 ENDDO
1284 ENDIF
1285 !
1286 ! CALCULATE CLOUD WORK FUNCTION AT T+DT
1287 !
1288 DO K = 1, KM
1289 DO I = 1, IM
1290 if (k .le. kmax(i)) then
1291 IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
1292 DZ1 = ZO(I,k) - ZO(I,k-1)
1293 GAMMA = EL2ORC * QESO(I,k-1) / (TO(I,k-1)**2)
1294 RFACT = 1. + DELTA * CP * GAMMA &
1295 & * TO(I,k-1) / HVAP
1296 AA1(I) = AA1(I) + &
1297 & DZ1 * (G / (CP * TO(I,k-1))) &
1298 & * DBYO(I,k-1) / (1. + GAMMA) &
1299 & * RFACT
1300 val = 0.
1301 AA1(I)=AA1(I)+ &
1302 & DZ1 * G * DELTA * &
1303 !cmr & MAX( 0.,(QESO(I,k-1) - QO(I,k-1))) &
1304 & MAX(val,(QESO(I,k-1) - QO(I,k-1)))
1305 ENDIF
1306 endif
1307 ENDDO
1308 ENDDO
1309 DO I = 1, IM
1310 IF(CNVFLG(I).AND.AA1(I).LE.0.) DWNFLG(I) = .FALSE.
1311 IF(CNVFLG(I).AND.AA1(I).LE.0.) DWNFLG2(I) = .FALSE.
1312 IF(CNVFLG(I).AND.AA1(I).LE.0.) CNVFLG(I) = .FALSE.
1313 ENDDO
1314 !!
1315 TOTFLG = .TRUE.
1316 DO I = 1, IM
1317 TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
1318 ENDDO
1319 IF(TOTFLG) RETURN
1320 !!
1321 !cccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
1322 !cccc PRINT *, ' AA1(I) BEFORE DWNDRFT =', AA1(I)
1323 !cccc ENDIF
1324 !
1325 !------- DOWNDRAFT CALCULATIONS
1326 !
1327 !
1328 !--- DETERMINE DOWNDRAFT STRENGTH IN TERMS OF WINDSHEAR
1329 !
1330 DO I = 1, IM
1331 IF(CNVFLG(I)) THEN
1332 VSHEAR(I) = 0.
1333 ENDIF
1334 ENDDO
1335 DO K = 1, KM
1336 DO I = 1, IM
1337 if (k .le. kmax(i)) then
1338 IF(K.GE.KB(I).AND.K.LE.KTCON(I).AND.CNVFLG(I)) THEN
1339 shear=rcs(I) * sqrt((UO(I,k+1)-UO(I,k)) ** 2 &
1340 & + (VO(I,k+1)-VO(I,k)) ** 2)
1341 VSHEAR(I) = VSHEAR(I) + SHEAR
1342 ENDIF
1343 endif
1344 ENDDO
1345 ENDDO
1346 DO I = 1, IM
1347 EDT(I) = 0.
1348 IF(CNVFLG(I)) THEN
1349 KNUMB = KTCON(I) - KB(I) + 1
1350 KNUMB = MAX(KNUMB,1)
1351 VSHEAR(I) = 1.E3 * VSHEAR(I) / (ZO(I,KTCON(I))-ZO(I,KB(I)))
1352 E1=1.591-.639*VSHEAR(I) &
1353 & +.0953*(VSHEAR(I)**2)-.00496*(VSHEAR(I)**3)
1354 EDT(I)=1.-E1
1355 !cmr EDT(I) = MIN(EDT(I),.9)
1356 val = .9
1357 EDT(I) = MIN(EDT(I),val)
1358 !cmr EDT(I) = MAX(EDT(I),.0)
1359 val = .0
1360 EDT(I) = MAX(EDT(I),val)
1361 EDTO(I)=EDT(I)
1362 EDTX(I)=EDT(I)
1363 ENDIF
1364 ENDDO
1365 ! DETERMINE DETRAINMENT RATE BETWEEN 1 AND KBDTR
1366 DO I = 1, IM
1367 KBDTR(I) = KBCON(I)
1368 beta = betas
1369 if(SLIMSK(I).eq.1.) beta = betal
1370 IF(CNVFLG(I)) THEN
1371 KBDTR(I) = KBCON(I)
1372 KBDTR(I) = MAX(KBDTR(I),1)
1373 XLAMD(I) = 0.
1374 IF(KBDTR(I).GT.1) THEN
1375 DZ = .5 * ZO(I,KBDTR(I)) + .5 * ZO(I,KBDTR(I)-1) &
1376 & - ZO(I,1)
1377 XLAMD(I) = LOG(BETA) / DZ
1378 ENDIF
1379 ENDIF
1380 ENDDO
1381 ! DETERMINE DOWNDRAFT MASS FLUX
1382 DO K = 1, KM
1383 DO I = 1, IM
1384 IF(k .le. kmax(i)) then
1385 IF(CNVFLG(I)) THEN
1386 ETAD(I,k) = 1.
1387 ENDIF
1388 QRCDO(I,k) = 0.
1389 endif
1390 ENDDO
1391 ENDDO
1392 DO K = KM1, 2, -1
1393 DO I = 1, IM
1394 if (k .le. kbmax(i)) then
1395 IF(CNVFLG(I).AND.K.LT.KBDTR(I)) THEN
1396 DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1397 ETAD(I,k) = ETAD(I,k+1) * EXP(XLAMD(I) * DZ)
1398 ENDIF
1399 endif
1400 ENDDO
1401 ENDDO
1402 K = 1
1403 DO I = 1, IM
1404 IF(CNVFLG(I).AND.KBDTR(I).GT.1) THEN
1405 DZ = .5 * (ZO(I,2) - ZO(I,1))
1406 ETAD(I,k) = ETAD(I,k+1) * EXP(XLAMD(I) * DZ)
1407 ENDIF
1408 ENDDO
1409 !
1410 !--- DOWNDRAFT MOISTURE PROPERTIES
1411 !
1412 DO I = 1, IM
1413 PWEVO(I) = 0.
1414 FLG(I) = CNVFLG(I)
1415 ENDDO
1416 DO I = 1, IM
1417 IF(CNVFLG(I)) THEN
1418 JMN = JMIN(I)
1419 HCDO(I) = HEO(I,JMN)
1420 QCDO(I) = QO(I,JMN)
1421 QRCDO(I,JMN) = QESO(I,JMN)
1422 UCDO(I) = UO(I,JMN)
1423 VCDO(I) = VO(I,JMN)
1424 ENDIF
1425 ENDDO
1426 DO K = KM1, 1, -1
1427 DO I = 1, IM
1428 if (k .le. kmax(i)-1) then
1429 IF(CNVFLG(I).AND.K.LT.JMIN(I)) THEN
1430 DQ = QESO(I,k)
1431 DT = TO(I,k)
1432 GAMMA = EL2ORC * DQ / DT**2
1433 DH = HCDO(I) - HESO(I,k)
1434 QRCDO(I,k) = DQ+(1./HVAP)*(GAMMA/(1.+GAMMA))*DH
1435 DETAD = ETAD(I,k+1) - ETAD(I,k)
1436 PWDO(I,k) = ETAD(I,k+1) * QCDO(I) - &
1437 & ETAD(I,k) * QRCDO(I,k)
1438 PWDO(I,k) = PWDO(I,k) - DETAD * &
1439 & .5 * (QRCDO(I,k) + QRCDO(I,k+1))
1440 QCDO(I) = QRCDO(I,k)
1441 PWEVO(I) = PWEVO(I) + PWDO(I,k)
1442 ENDIF
1443 endif
1444 ENDDO
1445 ENDDO
1446 ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.DWNFLG(I)) THEN
1447 ! PRINT *, ' PWAVO(I), PWEVO(I) =', PWAVO(I), PWEVO(I)
1448 ! ENDIF
1449 !
1450 !--- FINAL DOWNDRAFT STRENGTH DEPENDENT ON PRECIP
1451 !--- EFFICIENCY (EDT), NORMALIZED CONDENSATE (PWAV), AND
1452 !--- EVAPORATE (PWEV)
1453 !
1454 DO I = 1, IM
1455 edtmax = edtmaxl
1456 if(SLIMSK(I).eq.0.) edtmax = edtmaxs
1457 IF(DWNFLG2(I)) THEN
1458 IF(PWEVO(I).LT.0.) THEN
1459 EDTO(I) = -EDTO(I) * PWAVO(I) / PWEVO(I)
1460 EDTO(I) = MIN(EDTO(I),EDTMAX)
1461 ELSE
1462 EDTO(I) = 0.
1463 ENDIF
1464 ELSE
1465 EDTO(I) = 0.
1466 ENDIF
1467 ENDDO
1468 !
1469 !
1470 !--- DOWNDRAFT CLOUDWORK FUNCTIONS
1471 !
1472 !
1473 DO K = KM1, 1, -1
1474 DO I = 1, IM
1475 if (k .le. kmax(i)-1) then
1476 IF(DWNFLG2(I).AND.K.LT.JMIN(I)) THEN
1477 GAMMA = EL2ORC * QESO(I,k+1) / TO(I,k+1)**2
1478 DHH=HCDO(I)
1479 DT=TO(I,k+1)
1480 DG=GAMMA
1481 DH=HESO(I,k+1)
1482 DZ=-1.*(ZO(I,k+1)-ZO(I,k))
1483 AA1(I)=AA1(I)+EDTO(I)*DZ*(G/(CP*DT))*((DHH-DH)/(1.+DG)) &
1484 & *(1.+DELTA*CP*DG*DT/HVAP)
1485 val=0.
1486 AA1(I)=AA1(I)+EDTO(I)* &
1487 !cmr & DZ*G*DELTA*MAX( 0.,(QESO(I,k+1)-QO(I,k+1))) &
1488 & DZ*G*DELTA*MAX(val,(QESO(I,k+1)-QO(I,k+1)))
1489 ENDIF
1490 endif
1491 ENDDO
1492 ENDDO
1493 !cccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.DWNFLG2(I)) THEN
1494 !cccc PRINT *, ' AA1(I) AFTER DWNDRFT =', AA1(I)
1495 !cccc ENDIF
1496 DO I = 1, IM
1497 IF(AA1(I).LE.0.) CNVFLG(I) = .FALSE.
1498 IF(AA1(I).LE.0.) DWNFLG(I) = .FALSE.
1499 IF(AA1(I).LE.0.) DWNFLG2(I) = .FALSE.
1500 ENDDO
1501 !!
1502 TOTFLG = .TRUE.
1503 DO I = 1, IM
1504 TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
1505 ENDDO
1506 IF(TOTFLG) RETURN
1507 !!
1508 !
1509 !
1510 !--- WHAT WOULD THE CHANGE BE, THAT A CLOUD WITH UNIT MASS
1511 !--- WILL DO TO THE ENVIRONMENT?
1512 !
1513 DO K = 1, KM
1514 DO I = 1, IM
1515 IF(k .le. kmax(i) .and. CNVFLG(I)) THEN
1516 DELLAH(I,k) = 0.
1517 DELLAQ(I,k) = 0.
1518 DELLAU(I,k) = 0.
1519 DELLAV(I,k) = 0.
1520 ENDIF
1521 ENDDO
1522 ENDDO
1523 DO I = 1, IM
1524 IF(CNVFLG(I)) THEN
1525 DP = 1000. * DEL(I,1)
1526 DELLAH(I,1) = EDTO(I) * ETAD(I,1) * (HCDO(I) &
1527 & - HEO(I,1)) * G / DP
1528 DELLAQ(I,1) = EDTO(I) * ETAD(I,1) * (QCDO(I) &
1529 & - QO(I,1)) * G / DP
1530 DELLAU(I,1) = EDTO(I) * ETAD(I,1) * (UCDO(I) &
1531 & - UO(I,1)) * G / DP
1532 DELLAV(I,1) = EDTO(I) * ETAD(I,1) * (VCDO(I) &
1533 & - VO(I,1)) * G / DP
1534 ENDIF
1535 ENDDO
1536 !
1537 !--- CHANGED DUE TO SUBSIDENCE AND ENTRAINMENT
1538 !
1539 DO K = 2, KM1
1540 DO I = 1, IM
1541 if (k .le. kmax(i)-1) then
1542 IF(CNVFLG(I).AND.K.LT.KTCON(I)) THEN
1543 AUP = 1.
1544 IF(K.LE.KB(I)) AUP = 0.
1545 ADW = 1.
1546 IF(K.GT.JMIN(I)) ADW = 0.
1547 DV1= HEO(I,k)
1548 DV2 = .5 * (HEO(I,k) + HEO(I,k+1))
1549 DV3= HEO(I,k-1)
1550 DV1Q= QO(I,k)
1551 DV2Q = .5 * (QO(I,k) + QO(I,k+1))
1552 DV3Q= QO(I,k-1)
1553 DV1U= UO(I,k)
1554 DV2U = .5 * (UO(I,k) + UO(I,k+1))
1555 DV3U= UO(I,k-1)
1556 DV1V= VO(I,k)
1557 DV2V = .5 * (VO(I,k) + VO(I,k+1))
1558 DV3V= VO(I,k-1)
1559 DP = 1000. * DEL(I,K)
1560 DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1561 DETA = ETA(I,k) - ETA(I,k-1)
1562 DETAD = ETAD(I,k) - ETAD(I,k-1)
1563 DELLAH(I,k) = DELLAH(I,k) + &
1564 & ((AUP * ETA(I,k) - ADW * EDTO(I) * ETAD(I,k)) * DV1 &
1565 & - (AUP * ETA(I,k-1) - ADW * EDTO(I) * ETAD(I,k-1))* DV3 &
1566 & - AUP * DETA * DV2 &
1567 & + ADW * EDTO(I) * DETAD * HCDO(I)) * G / DP
1568 DELLAQ(I,k) = DELLAQ(I,k) + &
1569 & ((AUP * ETA(I,k) - ADW * EDTO(I) * ETAD(I,k)) * DV1Q &
1570 & - (AUP * ETA(I,k-1) - ADW * EDTO(I) * ETAD(I,k-1))* DV3Q &
1571 & - AUP * DETA * DV2Q &
1572 & +ADW*EDTO(I)*DETAD*.5*(QRCDO(I,k)+QRCDO(I,k-1))) * G / DP
1573 DELLAU(I,k) = DELLAU(I,k) + &
1574 & ((AUP * ETA(I,k) - ADW * EDTO(I) * ETAD(I,k)) * DV1U &
1575 & - (AUP * ETA(I,k-1) - ADW * EDTO(I) * ETAD(I,k-1))* DV3U &
1576 & - AUP * DETA * DV2U &
1577 & + ADW * EDTO(I) * DETAD * UCDO(I) &
1578 & ) * G / DP
1579 DELLAV(I,k) = DELLAV(I,k) + &
1580 & ((AUP * ETA(I,k) - ADW * EDTO(I) * ETAD(I,k)) * DV1V &
1581 & - (AUP * ETA(I,k-1) - ADW * EDTO(I) * ETAD(I,k-1))* DV3V &
1582 & - AUP * DETA * DV2V &
1583 & + ADW * EDTO(I) * DETAD * VCDO(I) &
1584 & ) * G / DP
1585 ENDIF
1586 endif
1587 ENDDO
1588 ENDDO
1589 !
1590 !------- CLOUD TOP
1591 !
1592 DO I = 1, IM
1593 IF(CNVFLG(I)) THEN
1594 INDX = KTCON(I)
1595 DP = 1000. * DEL(I,INDX)
1596 DV1 = HEO(I,INDX-1)
1597 DELLAH(I,INDX) = ETA(I,INDX-1) * &
1598 & (HCKO(I,INDX-1) - DV1) * G / DP
1599 DVQ1 = QO(I,INDX-1)
1600 DELLAQ(I,INDX) = ETA(I,INDX-1) * &
1601 & (QCKO(I,INDX-1) - DVQ1) * G / DP
1602 DV1U = UO(I,INDX-1)
1603 DELLAU(I,INDX) = ETA(I,INDX-1) * &
1604 & (UCKO(I,INDX-1) - DV1U) * G / DP
1605 DV1V = VO(I,INDX-1)
1606 DELLAV(I,INDX) = ETA(I,INDX-1) * &
1607 & (VCKO(I,INDX-1) - DV1V) * G / DP
1608 !
1609 ! cloud water
1610 !
1611 DELLAL(I) = ETA(I,INDX-1) * QLKO_KTCON(I) * g / dp
1612 ENDIF
1613 ENDDO
1614 !
1615 !------- FINAL CHANGED VARIABLE PER UNIT MASS FLUX
1616 !
1617 DO K = 1, KM
1618 DO I = 1, IM
1619 if (k .le. kmax(i)) then
1620 IF(CNVFLG(I).and.k.gt.KTCON(I)) THEN
1621 QO(I,k) = Q1(I,k)
1622 TO(I,k) = T1(I,k)
1623 UO(I,k) = U1(I,k)
1624 VO(I,k) = V1(I,k)
1625 ENDIF
1626 IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
1627 QO(I,k) = DELLAQ(I,k) * MBDT + Q1(I,k)
1628 DELLAT = (DELLAH(I,k) - HVAP * DELLAQ(I,k)) / CP
1629 TO(I,k) = DELLAT * MBDT + T1(I,k)
1630 !cmr QO(I,k) = max(QO(I,k),1.e-10)
1631 val = 1.e-10
1632 QO(I,k) = max(QO(I,k), val )
1633 ENDIF
1634 endif
1635 ENDDO
1636 ENDDO
1637 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1638 !
1639 !--- THE ABOVE CHANGED ENVIRONMENT IS NOW USED TO CALULATE THE
1640 !--- EFFECT THE ARBITRARY CLOUD (WITH UNIT MASS FLUX)
1641 !--- WOULD HAVE ON THE STABILITY,
1642 !--- WHICH THEN IS USED TO CALCULATE THE REAL MASS FLUX,
1643 !--- NECESSARY TO KEEP THIS CHANGE IN BALANCE WITH THE LARGE-SCALE
1644 !--- DESTABILIZATION.
1645 !
1646 !--- ENVIRONMENTAL CONDITIONS AGAIN, FIRST HEIGHTS
1647 !
1648 DO K = 1, KM
1649 DO I = 1, IM
1650 IF(k .le. kmax(i) .and. CNVFLG(I)) THEN
1651 !jfe QESO(I,k) = 10. * FPVS(TO(I,k))
1652 !
1653 QESO(I,k) = 0.01 * fpvs(TO(I,K)) ! fpvs is in Pa
1654 !
1655 QESO(I,k) = EPS * QESO(I,k) / (PFLD(I,k)+EPSM1*QESO(I,k))
1656 !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
1657 val = 1.E-8
1658 QESO(I,k) = MAX(QESO(I,k), val )
1659 TVO(I,k) = TO(I,k) + DELTA * TO(I,k) * QO(I,k)
1660 ENDIF
1661 ENDDO
1662 ENDDO
1663 DO I = 1, IM
1664 IF(CNVFLG(I)) THEN
1665 XAA0(I) = 0.
1666 XPWAV(I) = 0.
1667 ENDIF
1668 ENDDO
1669 !
1670 ! HYDROSTATIC HEIGHT ASSUME ZERO TERR
1671 !
1672 ! DO I = 1, IM
1673 ! IF(CNVFLG(I)) THEN
1674 ! DLNSIG = LOG(PRSL(I,1)/PS(I))
1675 ! ZO(I,1) = TERR - DLNSIG * RD / G * TVO(I,1)
1676 ! ENDIF
1677 ! ENDDO
1678 ! DO K = 2, KM
1679 ! DO I = 1, IM
1680 ! IF(k .le. kmax(i) .and. CNVFLG(I)) THEN
1681 ! DLNSIG = LOG(PRSL(I,K) / PRSL(I,K-1))
1682 ! ZO(I,k) = ZO(I,k-1) - DLNSIG * RD / G
1683 ! & * .5 * (TVO(I,k) + TVO(I,k-1))
1684 ! ENDIF
1685 ! ENDDO
1686 ! ENDDO
1687 !
1688 !--- MOIST STATIC ENERGY
1689 !
1690 DO K = 1, KM1
1691 DO I = 1, IM
1692 IF(k .le. kmax(i)-1 .and. CNVFLG(I)) THEN
1693 DZ = .5 * (ZO(I,k+1) - ZO(I,k))
1694 DP = .5 * (PFLD(I,k+1) - PFLD(I,k))
1695 !jfe ES = 10. * FPVS(TO(I,k+1))
1696 !
1697 ES = 0.01 * fpvs(TO(I,K+1)) ! fpvs is in Pa
1698 !
1699 PPRIME = PFLD(I,k+1) + EPSM1 * ES
1700 QS = EPS * ES / PPRIME
1701 DQSDP = - QS / PPRIME
1702 DESDT = ES * (FACT1 / TO(I,k+1) + FACT2 / (TO(I,k+1)**2))
1703 DQSDT = QS * PFLD(I,k+1) * DESDT / (ES * PPRIME)
1704 GAMMA = EL2ORC * QESO(I,k+1) / (TO(I,k+1)**2)
1705 DT = (G * DZ + HVAP * DQSDP * DP) / (CP * (1. + GAMMA))
1706 DQ = DQSDT * DT + DQSDP * DP
1707 TO(I,k) = TO(I,k+1) + DT
1708 QO(I,k) = QO(I,k+1) + DQ
1709 PO(I,k) = .5 * (PFLD(I,k) + PFLD(I,k+1))
1710 ENDIF
1711 ENDDO
1712 ENDDO
1713 DO K = 1, KM1
1714 DO I = 1, IM
1715 IF(k .le. kmax(i)-1 .and. CNVFLG(I)) THEN
1716 !jfe QESO(I,k) = 10. * FPVS(TO(I,k))
1717 !
1718 QESO(I,k) = 0.01 * fpvs(TO(I,K)) ! fpvs is in Pa
1719 !
1720 QESO(I,k) = EPS * QESO(I,k) / (PO(I,k) + EPSM1 * QESO(I,k))
1721 !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
1722 val1 = 1.E-8
1723 QESO(I,k) = MAX(QESO(I,k), val1)
1724 !cmr QO(I,k) = max(QO(I,k),1.e-10)
1725 val2 = 1.e-10
1726 QO(I,k) = max(QO(I,k), val2 )
1727 ! QO(I,k) = MIN(QO(I,k),QESO(I,k))
1728 HEO(I,k) = .5 * G * (ZO(I,k) + ZO(I,k+1)) + &
1729 & CP * TO(I,k) + HVAP * QO(I,k)
1730 HESO(I,k) = .5 * G * (ZO(I,k) + ZO(I,k+1)) + &
1731 & CP * TO(I,k) + HVAP * QESO(I,k)
1732 ENDIF
1733 ENDDO
1734 ENDDO
1735 DO I = 1, IM
1736 k = kmax(i)
1737 IF(CNVFLG(I)) THEN
1738 HEO(I,k) = G * ZO(I,k) + CP * TO(I,k) + HVAP * QO(I,k)
1739 HESO(I,k) = G * ZO(I,k) + CP * TO(I,k) + HVAP * QESO(I,k)
1740 ! HEO(I,k) = MIN(HEO(I,k),HESO(I,k))
1741 ENDIF
1742 ENDDO
1743 DO I = 1, IM
1744 IF(CNVFLG(I)) THEN
1745 INDX = KB(I)
1746 XHKB(I) = HEO(I,INDX)
1747 XQKB(I) = QO(I,INDX)
1748 HCKO(I,INDX) = XHKB(I)
1749 QCKO(I,INDX) = XQKB(I)
1750 ENDIF
1751 ENDDO
1752 !
1753 !
1754 !**************************** STATIC CONTROL
1755 !
1756 !
1757 !------- MOISTURE AND CLOUD WORK FUNCTIONS
1758 !
1759 DO K = 2, KM1
1760 DO I = 1, IM
1761 if (k .le. kmax(i)-1) then
1762 ! IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LE.KBCON(I)) THEN
1763 IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LE.KTCON(I)) THEN
1764 FACTOR = ETA(I,k-1) / ETA(I,k)
1765 ONEMF = 1. - FACTOR
1766 HCKO(I,k) = FACTOR * HCKO(I,k-1) + ONEMF * &
1767 & .5 * (HEO(I,k) + HEO(I,k+1))
1768 ENDIF
1769 ! IF(CNVFLG(I).AND.K.GT.KBCON(I)) THEN
1770 ! HEO(I,k) = HEO(I,k-1)
1771 ! ENDIF
1772 endif
1773 ENDDO
1774 ENDDO
1775 DO K = 2, KM1
1776 DO I = 1, IM
1777 if (k .le. kmax(i)-1) then
1778 IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LT.KTCON(I)) THEN
1779 DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1780 GAMMA = EL2ORC * QESO(I,k) / (TO(I,k)**2)
1781 XDBY = HCKO(I,k) - HESO(I,k)
1782 !cmr XDBY = MAX(XDBY,0.)
1783 val = 0.
1784 XDBY = MAX(XDBY,val)
1785 XQRCH = QESO(I,k) &
1786 & + GAMMA * XDBY / (HVAP * (1. + GAMMA))
1787 FACTOR = ETA(I,k-1) / ETA(I,k)
1788 ONEMF = 1. - FACTOR
1789 QCKO(I,k) = FACTOR * QCKO(I,k-1) + ONEMF * &
1790 & .5 * (QO(I,k) + QO(I,k+1))
1791 DQ = ETA(I,k) * QCKO(I,k) - ETA(I,k) * XQRCH
1792 IF(DQ.GT.0.) THEN
1793 ETAH = .5 * (ETA(I,k) + ETA(I,k-1))
1794 QLK = DQ / (ETA(I,k) + ETAH * C0 * DZ)
1795 XAA0(I) = XAA0(I) - (ZO(I,k) - ZO(I,k-1)) * G * QLK
1796 XQC = QLK + XQRCH
1797 XPW = ETAH * C0 * DZ * QLK
1798 QCKO(I,k) = XQC
1799 XPWAV(I) = XPWAV(I) + XPW
1800 ENDIF
1801 ENDIF
1802 ! IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LT.KTCON(I)) THEN
1803 IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
1804 DZ1 = ZO(I,k) - ZO(I,k-1)
1805 GAMMA = EL2ORC * QESO(I,k-1) / (TO(I,k-1)**2)
1806 RFACT = 1. + DELTA * CP * GAMMA &
1807 & * TO(I,k-1) / HVAP
1808 XDBY = HCKO(I,k-1) - HESO(I,k-1)
1809 XAA0(I) = XAA0(I) &
1810 & + DZ1 * (G / (CP * TO(I,k-1))) &
1811 & * XDBY / (1. + GAMMA) &
1812 & * RFACT
1813 val=0.
1814 XAA0(I)=XAA0(I)+ &
1815 & DZ1 * G * DELTA * &
1816 !cmr & MAX( 0.,(QESO(I,k-1) - QO(I,k-1))) &
1817 & MAX(val,(QESO(I,k-1) - QO(I,k-1)))
1818 ENDIF
1819 endif
1820 ENDDO
1821 ENDDO
1822 !cccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
1823 !cccc PRINT *, ' XAA BEFORE DWNDRFT =', XAA0(I)
1824 !cccc ENDIF
1825 !
1826 !------- DOWNDRAFT CALCULATIONS
1827 !
1828 !
1829 !--- DOWNDRAFT MOISTURE PROPERTIES
1830 !
1831 DO I = 1, IM
1832 XPWEV(I) = 0.
1833 ENDDO
1834 DO I = 1, IM
1835 IF(DWNFLG2(I)) THEN
1836 JMN = JMIN(I)
1837 XHCD(I) = HEO(I,JMN)
1838 XQCD(I) = QO(I,JMN)
1839 QRCD(I,JMN) = QESO(I,JMN)
1840 ENDIF
1841 ENDDO
1842 DO K = KM1, 1, -1
1843 DO I = 1, IM
1844 if (k .le. kmax(i)-1) then
1845 IF(DWNFLG2(I).AND.K.LT.JMIN(I)) THEN
1846 DQ = QESO(I,k)
1847 DT = TO(I,k)
1848 GAMMA = EL2ORC * DQ / DT**2
1849 DH = XHCD(I) - HESO(I,k)
1850 QRCD(I,k)=DQ+(1./HVAP)*(GAMMA/(1.+GAMMA))*DH
1851 DETAD = ETAD(I,k+1) - ETAD(I,k)
1852 XPWD = ETAD(I,k+1) * QRCD(I,k+1) - &
1853 & ETAD(I,k) * QRCD(I,k)
1854 XPWD = XPWD - DETAD * &
1855 & .5 * (QRCD(I,k) + QRCD(I,k+1))
1856 XPWEV(I) = XPWEV(I) + XPWD
1857 ENDIF
1858 endif
1859 ENDDO
1860 ENDDO
1861 !
1862 DO I = 1, IM
1863 edtmax = edtmaxl
1864 if(SLIMSK(I).eq.0.) edtmax = edtmaxs
1865 IF(DWNFLG2(I)) THEN
1866 IF(XPWEV(I).GE.0.) THEN
1867 EDTX(I) = 0.
1868 ELSE
1869 EDTX(I) = -EDTX(I) * XPWAV(I) / XPWEV(I)
1870 EDTX(I) = MIN(EDTX(I),EDTMAX)
1871 ENDIF
1872 ELSE
1873 EDTX(I) = 0.
1874 ENDIF
1875 ENDDO
1876 !
1877 !
1878 !
1879 !--- DOWNDRAFT CLOUDWORK FUNCTIONS
1880 !
1881 !
1882 DO K = KM1, 1, -1
1883 DO I = 1, IM
1884 if (k .le. kmax(i)-1) then
1885 IF(DWNFLG2(I).AND.K.LT.JMIN(I)) THEN
1886 GAMMA = EL2ORC * QESO(I,k+1) / TO(I,k+1)**2
1887 DHH=XHCD(I)
1888 DT= TO(I,k+1)
1889 DG= GAMMA
1890 DH= HESO(I,k+1)
1891 DZ=-1.*(ZO(I,k+1)-ZO(I,k))
1892 XAA0(I)=XAA0(I)+EDTX(I)*DZ*(G/(CP*DT))*((DHH-DH)/(1.+DG)) &
1893 & *(1.+DELTA*CP*DG*DT/HVAP)
1894 val=0.
1895 XAA0(I)=XAA0(I)+EDTX(I)* &
1896 !cmr & DZ*G*DELTA*MAX( 0.,(QESO(I,k+1)-QO(I,k+1))) &
1897 & DZ*G*DELTA*MAX(val,(QESO(I,k+1)-QO(I,k+1)))
1898 ENDIF
1899 endif
1900 ENDDO
1901 ENDDO
1902 !cccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.DWNFLG2(I)) THEN
1903 !cccc PRINT *, ' XAA AFTER DWNDRFT =', XAA0(I)
1904 !cccc ENDIF
1905 !
1906 ! CALCULATE CRITICAL CLOUD WORK FUNCTION
1907 !
1908 DO I = 1, IM
1909 ACRT(I) = 0.
1910 IF(CNVFLG(I)) THEN
1911 ! IF(CNVFLG(I).AND.SLIMSK(I).NE.1.) THEN
1912 IF(PFLD(I,KTCON(I)).LT.PCRIT(15))THEN
1913 ACRT(I)=ACRIT(15)*(975.-PFLD(I,KTCON(I))) &
1914 & /(975.-PCRIT(15))
1915 ELSE IF(PFLD(I,KTCON(I)).GT.PCRIT(1))THEN
1916 ACRT(I)=ACRIT(1)
1917 ELSE
1918 !cmr K = IFIX((850. - PFLD(I,KTCON(I)))/50.) + 2
1919 K = int((850. - PFLD(I,KTCON(I)))/50.) + 2
1920 K = MIN(K,15)
1921 K = MAX(K,2)
1922 ACRT(I)=ACRIT(K)+(ACRIT(K-1)-ACRIT(K))* &
1923 & (PFLD(I,KTCON(I))-PCRIT(K))/(PCRIT(K-1)-PCRIT(K))
1924 ENDIF
1925 ! ELSE
1926 ! ACRT(I) = .5 * (PFLD(I,KBCON(I)) - PFLD(I,KTCON(I)))
1927 ENDIF
1928 ENDDO
1929 DO I = 1, IM
1930 ACRTFCT(I) = 1.
1931 IF(CNVFLG(I)) THEN
1932 if(SLIMSK(I).eq.1.) THEN
1933 w1 = w1l
1934 w2 = w2l
1935 w3 = w3l
1936 w4 = w4l
1937 else
1938 w1 = w1s
1939 w2 = w2s
1940 w3 = w3s
1941 w4 = w4s
1942 ENDIF
1943 !C IF(CNVFLG(I).AND.SLIMSK(I).EQ.1.) THEN
1944 ! ACRTFCT(I) = PDOT(I) / W3
1945 !
1946 ! modify critical cloud workfunction by cloud base vertical velocity
1947 !
1948 IF(PDOT(I).LE.W4) THEN
1949 ACRTFCT(I) = (PDOT(I) - W4) / (W3 - W4)
1950 ELSEIF(PDOT(I).GE.-W4) THEN
1951 ACRTFCT(I) = - (PDOT(I) + W4) / (W4 - W3)
1952 ELSE
1953 ACRTFCT(I) = 0.
1954 ENDIF
1955 !cmr ACRTFCT(I) = MAX(ACRTFCT(I),-1.)
1956 val1 = -1.
1957 ACRTFCT(I) = MAX(ACRTFCT(I),val1)
1958 !cmr ACRTFCT(I) = MIN(ACRTFCT(I),1.)
1959 val2 = 1.
1960 ACRTFCT(I) = MIN(ACRTFCT(I),val2)
1961 ACRTFCT(I) = 1. - ACRTFCT(I)
1962 !
1963 ! modify ACRTFCT(I) by colume mean rh if RHBAR(I) is greater than 80 percent
1964 !
1965 ! if(RHBAR(I).ge..8) THEN
1966 ! ACRTFCT(I) = ACRTFCT(I) * (.9 - min(RHBAR(I),.9)) * 10.
1967 ! ENDIF
1968 !
1969 ! modify adjustment time scale by cloud base vertical velocity
1970 !
1971 DTCONV(I) = DT2 + max((1800. - DT2),RZERO) * &
1972 & (PDOT(I) - W2) / (W1 - W2)
1973 ! DTCONV(I) = MAX(DTCONV(I), DT2)
1974 ! DTCONV(I) = 1800. * (PDOT(I) - w2) / (w1 - w2)
1975 DTCONV(I) = max(DTCONV(I),dtmin)
1976 DTCONV(I) = min(DTCONV(I),dtmax)
1978 ENDIF
1979 ENDDO
1980 !
1981 !--- LARGE SCALE FORCING
1982 !
1983 DO I= 1, IM
1984 FLG(I) = CNVFLG(I)
1985 IF(CNVFLG(I)) THEN
1986 ! F = AA1(I) / DTCONV(I)
1987 FLD(I) = (AA1(I) - ACRT(I) * ACRTFCT(I)) / DTCONV(I)
1988 IF(FLD(I).LE.0.) FLG(I) = .FALSE.
1989 ENDIF
1990 CNVFLG(I) = FLG(I)
1991 IF(CNVFLG(I)) THEN
1992 ! XAA0(I) = MAX(XAA0(I),0.)
1993 XK(I) = (XAA0(I) - AA1(I)) / MBDT
1994 IF(XK(I).GE.0.) FLG(I) = .FALSE.
1995 ENDIF
1996 !
1997 !--- KERNEL, CLOUD BASE MASS FLUX
1998 !
1999 CNVFLG(I) = FLG(I)
2000 IF(CNVFLG(I)) THEN
2001 XMB(I) = -FLD(I) / XK(I)
2002 XMB(I) = MIN(XMB(I),XMBMAX(I))
2003 ENDIF
2004 ENDDO
2005 ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
2006 ! print *, ' RHBAR(I), ACRTFCT(I) =', RHBAR(I), ACRTFCT(I)
2007 ! PRINT *, ' A1, XA =', AA1(I), XAA0(I)
2008 ! PRINT *, ' XMB(I), ACRT =', XMB(I), ACRT
2009 ! ENDIF
2010 TOTFLG = .TRUE.
2011 DO I = 1, IM
2012 TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
2013 ENDDO
2014 IF(TOTFLG) RETURN
2015 !
2016 ! restore t0 and QO to t1 and q1 in case convection stops
2017 !
2018 do k = 1, km
2019 DO I = 1, IM
2020 if (k .le. kmax(i)) then
2021 TO(I,k) = T1(I,k)
2022 QO(I,k) = Q1(I,k)
2023 !jfe QESO(I,k) = 10. * FPVS(T1(I,k))
2024 !
2025 QESO(I,k) = 0.01 * fpvs(T1(I,K)) ! fpvs is in Pa
2026 !
2027 QESO(I,k) = EPS * QESO(I,k) / (PFLD(I,k) + EPSM1*QESO(I,k))
2028 !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
2029 val = 1.E-8
2030 QESO(I,k) = MAX(QESO(I,k), val )
2031 endif
2032 enddo
2033 enddo
2034 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2035 !
2036 !--- FEEDBACK: SIMPLY THE CHANGES FROM THE CLOUD WITH UNIT MASS FLUX
2037 !--- MULTIPLIED BY THE MASS FLUX NECESSARY TO KEEP THE
2038 !--- EQUILIBRIUM WITH THE LARGER-SCALE.
2039 !
2040 DO I = 1, IM
2041 DELHBAR(I) = 0.
2042 DELQBAR(I) = 0.
2043 DELTBAR(I) = 0.
2044 QCOND(I) = 0.
2045 ENDDO
2046 DO K = 1, KM
2047 DO I = 1, IM
2048 if (k .le. kmax(i)) then
2049 IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
2050 AUP = 1.
2051 IF(K.Le.KB(I)) AUP = 0.
2052 ADW = 1.
2053 IF(K.GT.JMIN(I)) ADW = 0.
2054 DELLAT = (DELLAH(I,k) - HVAP * DELLAQ(I,k)) / CP
2055 T1(I,k) = T1(I,k) + DELLAT * XMB(I) * DT2
2056 Q1(I,k) = Q1(I,k) + DELLAQ(I,k) * XMB(I) * DT2
2057 U1(I,k) = U1(I,k) + DELLAU(I,k) * XMB(I) * DT2
2058 V1(I,k) = V1(I,k) + DELLAV(I,k) * XMB(I) * DT2
2059 DP = 1000. * DEL(I,K)
2060 DELHBAR(I) = DELHBAR(I) + DELLAH(I,k)*XMB(I)*DP/G
2061 DELQBAR(I) = DELQBAR(I) + DELLAQ(I,k)*XMB(I)*DP/G
2062 DELTBAR(I) = DELTBAR(I) + DELLAT*XMB(I)*DP/G
2063 ENDIF
2064 endif
2065 ENDDO
2066 ENDDO
2067 DO K = 1, KM
2068 DO I = 1, IM
2069 if (k .le. kmax(i)) then
2070 IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
2071 !jfe QESO(I,k) = 10. * FPVS(T1(I,k))
2072 !
2073 QESO(I,k) = 0.01 * fpvs(T1(I,K)) ! fpvs is in Pa
2074 !
2075 QESO(I,k) = EPS * QESO(I,k)/(PFLD(I,k) + EPSM1*QESO(I,k))
2076 !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
2077 val = 1.E-8
2078 QESO(I,k) = MAX(QESO(I,k), val )
2079 !
2080 ! cloud water
2081 !
2082 if(ncloud.gt.0.and.CNVFLG(I).and.k.eq.KTCON(I)) THEN
2083 tem = DELLAL(I) * XMB(I) * dt2
2084 tem1 = MAX(RZERO, MIN(RONE, (TCR-t1(I,K))*TCRF))
2085 if (QL(I,k,2) .gt. -999.0) then
2086 QL(I,k,1) = QL(I,k,1) + tem * tem1 ! Ice
2087 QL(I,k,2) = QL(I,k,2) + tem *(1.0-tem1) ! Water
2088 else
2089 tem2 = QL(I,k,1) + tem
2090 QL(I,k,1) = tem2 * tem1 ! Ice
2091 QL(I,k,2) = tem2 - QL(I,k,1) ! Water
2092 endif
2093 ! QL(I,k) = QL(I,k) + DELLAL(I) * XMB(I) * dt2
2094 dp = 1000. * del(i,k)
2095 DELLAL(I) = DELLAL(I) * XMB(I) * dp / g
2096 ENDIF
2097 ENDIF
2098 endif
2099 ENDDO
2100 ENDDO
2101 ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I) ) THEN
2102 ! PRINT *, ' DELHBAR, DELQBAR, DELTBAR ='
2103 ! PRINT *, DELHBAR, HVAP*DELQBAR, CP*DELTBAR
2104 ! PRINT *, ' DELLBAR ='
2105 ! PRINT 6003, HVAP*DELLbar
2106 ! PRINT *, ' DELLAQ ='
2107 ! PRINT 6003, (HVAP*DELLAQ(I,k)*XMB(I),K=1,KMAX)
2108 ! PRINT *, ' DELLAT ='
2109 ! PRINT 6003, (DELLAH(i,k)*XMB(I)-HVAP*DELLAQ(I,k)*XMB(I), &
2110 ! & K=1,KMAX)
2111 ! ENDIF
2112 DO I = 1, IM
2113 RNTOT(I) = 0.
2114 DELQEV(I) = 0.
2115 DELQ2(I) = 0.
2116 FLG(I) = CNVFLG(I)
2117 ENDDO
2118 DO K = KM, 1, -1
2119 DO I = 1, IM
2120 if (k .le. kmax(i)) then
2121 IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
2122 AUP = 1.
2123 IF(K.Le.KB(I)) AUP = 0.
2124 ADW = 1.
2125 IF(K.GT.JMIN(I)) ADW = 0.
2126 rain = AUP * PWO(I,k) + ADW * EDTO(I) * PWDO(I,k)
2127 RNTOT(I) = RNTOT(I) + rain * XMB(I) * .001 * dt2
2128 ENDIF
2129 endif
2130 ENDDO
2131 ENDDO
2132 DO K = KM, 1, -1
2133 DO I = 1, IM
2134 if (k .le. kmax(i)) then
2135 DELTV(I) = 0.
2136 DELQ(I) = 0.
2137 QEVAP(I) = 0.
2138 IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
2139 AUP = 1.
2140 IF(K.Le.KB(I)) AUP = 0.
2141 ADW = 1.
2142 IF(K.GT.JMIN(I)) ADW = 0.
2143 rain = AUP * PWO(I,k) + ADW * EDTO(I) * PWDO(I,k)
2144 RN(I) = RN(I) + rain * XMB(I) * .001 * dt2
2145 ENDIF
2146 IF(FLG(I).AND.K.LE.KTCON(I)) THEN
2147 evef = EDT(I) * evfact
2148 if(SLIMSK(I).eq.1.) evef=EDT(I) * evfactl
2149 ! if(SLIMSK(I).eq.1.) evef=.07
2150 ! if(SLIMSK(I).ne.1.) evef = 0.
2151 QCOND(I) = EVEF * (Q1(I,k) - QESO(I,k)) &
2152 & / (1. + EL2ORC * QESO(I,k) / T1(I,k)**2)
2153 DP = 1000. * DEL(I,K)
2154 IF(RN(I).GT.0..AND.QCOND(I).LT.0.) THEN
2155 QEVAP(I) = -QCOND(I) * (1.-EXP(-.32*SQRT(DT2*RN(I))))
2156 QEVAP(I) = MIN(QEVAP(I), RN(I)*1000.*G/DP)
2157 DELQ2(I) = DELQEV(I) + .001 * QEVAP(I) * dp / g
2158 ENDIF
2159 if(RN(I).gt.0..and.QCOND(I).LT.0..and. &
2160 & DELQ2(I).gt.RNTOT(I)) THEN
2161 QEVAP(I) = 1000.* g * (RNTOT(I) - DELQEV(I)) / dp
2162 FLG(I) = .false.
2163 ENDIF
2164 IF(RN(I).GT.0..AND.QEVAP(I).gt.0.) THEN
2165 Q1(I,k) = Q1(I,k) + QEVAP(I)
2166 T1(I,k) = T1(I,k) - ELOCP * QEVAP(I)
2167 RN(I) = RN(I) - .001 * QEVAP(I) * DP / G
2168 DELTV(I) = - ELOCP*QEVAP(I)/DT2
2169 DELQ(I) = + QEVAP(I)/DT2
2170 DELQEV(I) = DELQEV(I) + .001*dp*QEVAP(I)/g
2171 ENDIF
2172 DELLAQ(I,k) = DELLAQ(I,k) + DELQ(I) / XMB(I)
2173 DELQBAR(I) = DELQBAR(I) + DELQ(I)*DP/G
2174 DELTBAR(I) = DELTBAR(I) + DELTV(I)*DP/G
2175 ENDIF
2176 endif
2177 ENDDO
2178 ENDDO
2179 ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I) ) THEN
2180 ! PRINT *, ' DELLAH ='
2181 ! PRINT 6003, (DELLAH(k)*XMB(I),K=1,KMAX)
2182 ! PRINT *, ' DELLAQ ='
2183 ! PRINT 6003, (HVAP*DELLAQ(I,k)*XMB(I),K=1,KMAX)
2184 ! PRINT *, ' DELHBAR, DELQBAR, DELTBAR ='
2185 ! PRINT *, DELHBAR, HVAP*DELQBAR, CP*DELTBAR
2186 ! PRINT *, ' PRECIP =', HVAP*RN(I)*1000./DT2
2187 !CCCC PRINT *, ' DELLBAR ='
2188 !CCCC PRINT *, HVAP*DELLbar
2189 ! ENDIF
2190 !
2191 ! PRECIPITATION RATE CONVERTED TO ACTUAL PRECIP
2192 ! IN UNIT OF M INSTEAD OF KG
2193 !
2194 DO I = 1, IM
2195 IF(CNVFLG(I)) THEN
2196 !
2197 ! IN THE EVENT OF UPPER LEVEL RAIN EVAPORATION AND LOWER LEVEL DOWNDRAF
2198 ! MOISTENING, RN CAN BECOME NEGATIVE, IN THIS CASE, WE BACK OUT OF TH
2199 ! HEATING AND THE MOISTENING
2200 !
2201 if(RN(I).lt.0..and..not.FLG(I)) RN(I) = 0.
2202 IF(RN(I).LE.0.) THEN
2203 RN(I) = 0.
2204 ELSE
2205 KTOP(I) = KTCON(I)
2206 KBOT(I) = KBCON(I)
2207 KUO(I) = 1
2208 CLDWRK(I) = AA1(I)
2209 ENDIF
2210 ENDIF
2211 ENDDO
2212 DO K = 1, KM
2213 DO I = 1, IM
2214 if (k .le. kmax(i)) then
2215 IF(CNVFLG(I).AND.RN(I).LE.0.) THEN
2216 T1(I,k) = TO(I,k)
2217 Q1(I,k) = QO(I,k)
2218 ENDIF
2219 endif
2220 ENDDO
2221 ENDDO
2222 !!
2223 RETURN
2224 END SUBROUTINE SASCNV
2226 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2228 SUBROUTINE SHALCV(IM,IX,KM,DT,DEL,PRSI,PRSL,PRSLK,KUO,Q,T,DPSHC)
2229 !
2230 USE MODULE_GFS_MACHINE , ONLY : kind_phys
2231 USE MODULE_GFS_PHYSCONS, grav => con_g, CP => con_CP, HVAP => con_HVAP &
2232 &, RD => con_RD
2234 implicit none
2235 !
2236 ! include 'constant.h'
2237 !
2238 integer IM, IX, KM, KUO(IM)
2239 real(kind=kind_phys) DEL(IX,KM), PRSI(IX,KM+1), PRSL(IX,KM), &
2240 & PRSLK(IX,KM), &
2241 & Q(IX,KM), T(IX,KM), DT, DPSHC
2242 !
2243 ! Locals
2244 !
2245 real(kind=kind_phys) ck, cpdt, dmse, dsdz1, dsdz2, &
2246 & dsig, dtodsl, dtodsu, eldq, g, &
2247 & gocp, rtdls
2248 !
2249 integer k,k1,k2,kliftl,kliftu,kt,N2,I,iku,ik1,ik,ii
2250 integer INDEX2(IM), KLCL(IM), KBOT(IM), KTOP(IM),kk &
2251 &, KTOPM(IM)
2252 !!
2253 ! PHYSICAL PARAMETERS
2254 PARAMETER(G=GRAV, GOCP=G/CP)
2255 ! BOUNDS OF PARCEL ORIGIN
2256 PARAMETER(KLIFTL=2,KLIFTU=2)
2257 LOGICAL LSHC(IM)
2258 real(kind=kind_phys) Q2(IM*KM), T2(IM*KM), &
2259 & PRSL2(IM*KM), PRSLK2(IM*KM), &
2260 & AL(IM*(KM-1)), AD(IM*KM), AU(IM*(KM-1))
2261 !-----------------------------------------------------------------------
2262 ! COMPRESS FIELDS TO POINTS WITH NO DEEP CONVECTION
2263 ! AND MOIST STATIC INSTABILITY.
2264 DO I=1,IM
2265 LSHC(I)=.FALSE.
2266 ENDDO
2267 DO K=1,KM-1
2268 DO I=1,IM
2269 IF(KUO(I).EQ.0) THEN
2270 ELDQ = HVAP*(Q(I,K)-Q(I,K+1))
2271 CPDT = CP*(T(I,K)-T(I,K+1))
2272 RTDLS = (PRSL(I,K)-PRSL(I,K+1)) / &
2273 & PRSI(I,K+1)*RD*0.5*(T(I,K)+T(I,K+1))
2274 DMSE = ELDQ+CPDT-RTDLS
2275 LSHC(I) = LSHC(I).OR.DMSE.GT.0.
2276 ENDIF
2277 ENDDO
2278 ENDDO
2279 N2 = 0
2280 DO I=1,IM
2281 IF(LSHC(I)) THEN
2282 N2 = N2 + 1
2283 INDEX2(N2) = I
2284 ENDIF
2285 ENDDO
2286 IF(N2.EQ.0) RETURN
2287 DO K=1,KM
2288 KK = (K-1)*N2
2289 DO I=1,N2
2290 IK = KK + I
2291 ii = index2(i)
2292 Q2(IK) = Q(II,K)
2293 T2(IK) = T(II,K)
2294 PRSL2(IK) = PRSL(II,K)
2295 PRSLK2(IK) = PRSLK(II,K)
2296 ENDDO
2297 ENDDO
2298 do i=1,N2
2299 ktopm(i) = KM
2300 enddo
2301 do k=2,KM
2302 do i=1,N2
2303 ii = index2(i)
2304 if (prsi(ii,1)-prsi(ii,k) .le. dpshc) ktopm(i) = k
2305 enddo
2306 enddo
2308 !-----------------------------------------------------------------------
2309 ! COMPUTE MOIST ADIABAT AND DETERMINE LIMITS OF SHALLOW CONVECTION.
2310 ! CHECK FOR MOIST STATIC INSTABILITY AGAIN WITHIN CLOUD.
2311 CALL MSTADBT3(N2,KM-1,KLIFTL,KLIFTU,PRSL2,PRSLK2,T2,Q2, &
2312 & KLCL,KBOT,KTOP,AL,AU)
2313 DO I=1,N2
2314 KBOT(I) = min(KLCL(I)-1, ktopm(i)-1)
2315 KTOP(I) = min(KTOP(I)+1, ktopm(i))
2316 LSHC(I) = .FALSE.
2317 ENDDO
2318 DO K=1,KM-1
2319 KK = (K-1)*N2
2320 DO I=1,N2
2321 IF(K.GE.KBOT(I).AND.K.LT.KTOP(I)) THEN
2322 IK = KK + I
2323 IKU = IK + N2
2324 ELDQ = HVAP * (Q2(IK)-Q2(IKU))
2325 CPDT = CP * (T2(IK)-T2(IKU))
2326 RTDLS = (PRSL2(IK)-PRSL2(IKU)) / &
2327 & PRSI(index2(i),K+1)*RD*0.5*(T2(IK)+T2(IKU))
2328 DMSE = ELDQ + CPDT - RTDLS
2329 LSHC(I) = LSHC(I).OR.DMSE.GT.0.
2330 AU(IK) = G/RTDLS
2331 ENDIF
2332 ENDDO
2333 ENDDO
2334 K1=KM+1
2335 K2=0
2336 DO I=1,N2
2337 IF(.NOT.LSHC(I)) THEN
2338 KBOT(I) = KM+1
2339 KTOP(I) = 0
2340 ENDIF
2341 K1 = MIN(K1,KBOT(I))
2342 K2 = MAX(K2,KTOP(I))
2343 ENDDO
2344 KT = K2-K1+1
2345 IF(KT.LT.2) RETURN
2346 !-----------------------------------------------------------------------
2347 ! SET EDDY VISCOSITY COEFFICIENT CKU AT SIGMA INTERFACES.
2348 ! COMPUTE DIAGONALS AND RHS FOR TRIDIAGONAL MATRIX SOLVER.
2349 ! EXPAND FINAL FIELDS.
2350 KK = (K1-1) * N2
2351 DO I=1,N2
2352 IK = KK + I
2353 AD(IK) = 1.
2354 ENDDO
2355 !
2356 ! DTODSU=DT/DEL(K1)
2357 DO K=K1,K2-1
2358 ! DTODSL=DTODSU
2359 ! DTODSU= DT/DEL(K+1)
2360 ! DSIG=SL(K)-SL(K+1)
2361 KK = (K-1) * N2
2362 DO I=1,N2
2363 ii = index2(i)
2364 DTODSL = DT/DEL(II,K)
2365 DTODSU = DT/DEL(II,K+1)
2366 DSIG = PRSL(II,K) - PRSL(II,K+1)
2367 IK = KK + I
2368 IKU = IK + N2
2369 IF(K.EQ.KBOT(I)) THEN
2370 CK=1.5
2371 ELSEIF(K.EQ.KTOP(I)-1) THEN
2372 CK=1.
2373 ELSEIF(K.EQ.KTOP(I)-2) THEN
2374 CK=3.
2375 ELSEIF(K.GT.KBOT(I).AND.K.LT.KTOP(I)-2) THEN
2376 CK=5.
2377 ELSE
2378 CK=0.
2379 ENDIF
2380 DSDZ1 = CK*DSIG*AU(IK)*GOCP
2381 DSDZ2 = CK*DSIG*AU(IK)*AU(IK)
2382 AU(IK) = -DTODSL*DSDZ2
2383 AL(IK) = -DTODSU*DSDZ2
2384 AD(IK) = AD(IK)-AU(IK)
2385 AD(IKU) = 1.-AL(IK)
2386 T2(IK) = T2(IK)+DTODSL*DSDZ1
2387 T2(IKU) = T2(IKU)-DTODSU*DSDZ1
2388 ENDDO
2389 ENDDO
2390 IK1=(K1-1)*N2+1
2391 CALL TRIDI2T3(N2,KT,AL(IK1),AD(IK1),AU(IK1),Q2(IK1),T2(IK1), &
2392 & AU(IK1),Q2(IK1),T2(IK1))
2393 DO K=K1,K2
2394 KK = (K-1)*N2
2395 DO I=1,N2
2396 IK = KK + I
2397 Q(INDEX2(I),K) = Q2(IK)
2398 T(INDEX2(I),K) = T2(IK)
2399 ENDDO
2400 ENDDO
2401 !-----------------------------------------------------------------------
2402 RETURN
2403 END SUBROUTINE SHALCV
2404 !-----------------------------------------------------------------------
2405 SUBROUTINE TRIDI2T3(L,N,CL,CM,CU,R1,R2,AU,A1,A2)
2406 !yt INCLUDE DBTRIDI2;
2407 !!
2408 USE MODULE_GFS_MACHINE , ONLY : kind_phys
2409 implicit none
2410 integer k,n,l,i
2411 real(kind=kind_phys) fk
2412 !!
2413 real(kind=kind_phys) &
2414 & CL(L,2:N),CM(L,N),CU(L,N-1),R1(L,N),R2(L,N), &
2415 & AU(L,N-1),A1(L,N),A2(L,N)
2416 !-----------------------------------------------------------------------
2417 DO I=1,L
2418 FK=1./CM(I,1)
2419 AU(I,1)=FK*CU(I,1)
2420 A1(I,1)=FK*R1(I,1)
2421 A2(I,1)=FK*R2(I,1)
2422 ENDDO
2423 DO K=2,N-1
2424 DO I=1,L
2425 FK=1./(CM(I,K)-CL(I,K)*AU(I,K-1))
2426 AU(I,K)=FK*CU(I,K)
2427 A1(I,K)=FK*(R1(I,K)-CL(I,K)*A1(I,K-1))
2428 A2(I,K)=FK*(R2(I,K)-CL(I,K)*A2(I,K-1))
2429 ENDDO
2430 ENDDO
2431 DO I=1,L
2432 FK=1./(CM(I,N)-CL(I,N)*AU(I,N-1))
2433 A1(I,N)=FK*(R1(I,N)-CL(I,N)*A1(I,N-1))
2434 A2(I,N)=FK*(R2(I,N)-CL(I,N)*A2(I,N-1))
2435 ENDDO
2436 DO K=N-1,1,-1
2437 DO I=1,L
2438 A1(I,K)=A1(I,K)-AU(I,K)*A1(I,K+1)
2439 A2(I,K)=A2(I,K)-AU(I,K)*A2(I,K+1)
2440 ENDDO
2441 ENDDO
2442 !-----------------------------------------------------------------------
2443 RETURN
2444 END SUBROUTINE TRIDI2T3
2445 !-----------------------------------------------------------------------
2447 SUBROUTINE MSTADBT3(IM,KM,K1,K2,PRSL,PRSLK,TENV,QENV, &
2448 & KLCL,KBOT,KTOP,TCLD,QCLD)
2449 !yt INCLUDE DBMSTADB;
2450 !!
2451 USE MODULE_GFS_MACHINE, ONLY : kind_phys
2452 USE MODULE_GFS_FUNCPHYS, ONLY : FTDP, FTHE, FTLCL, STMA
2453 USE MODULE_GFS_PHYSCONS, EPS => con_eps, EPSM1 => con_epsm1, FV => con_FVirt
2455 implicit none
2456 !!
2457 ! include 'constant.h'
2458 !!
2459 integer k,k1,k2,km,i,im
2460 real(kind=kind_phys) pv,qma,slklcl,tdpd,thelcl,tlcl
2461 real(kind=kind_phys) tma,tvcld,tvenv
2462 !!
2463 real(kind=kind_phys) PRSL(IM,KM), PRSLK(IM,KM), TENV(IM,KM), &
2464 & QENV(IM,KM), TCLD(IM,KM), QCLD(IM,KM)
2465 INTEGER KLCL(IM), KBOT(IM), KTOP(IM)
2466 ! LOCAL ARRAYS
2467 real(kind=kind_phys) SLKMA(IM), THEMA(IM)
2468 !-----------------------------------------------------------------------
2469 ! DETERMINE WARMEST POTENTIAL WET-BULB TEMPERATURE BETWEEN K1 AND K2.
2470 ! COMPUTE ITS LIFTING CONDENSATION LEVEL.
2471 !
2472 DO I=1,IM
2473 SLKMA(I) = 0.
2474 THEMA(I) = 0.
2475 ENDDO
2476 DO K=K1,K2
2477 DO I=1,IM
2478 PV = 1000.0 * PRSL(I,K)*QENV(I,K)/(EPS-EPSM1*QENV(I,K))
2479 TDPD = TENV(I,K)-FTDP(PV)
2480 IF(TDPD.GT.0.) THEN
2481 TLCL = FTLCL(TENV(I,K),TDPD)
2482 SLKLCL = PRSLK(I,K)*TLCL/TENV(I,K)
2483 ELSE
2484 TLCL = TENV(I,K)
2485 SLKLCL = PRSLK(I,K)
2486 ENDIF
2487 THELCL=FTHE(TLCL,SLKLCL)
2488 IF(THELCL.GT.THEMA(I)) THEN
2489 SLKMA(I) = SLKLCL
2490 THEMA(I) = THELCL
2491 ENDIF
2492 ENDDO
2493 ENDDO
2494 !-----------------------------------------------------------------------
2495 ! SET CLOUD TEMPERATURES AND HUMIDITIES WHEREVER THE PARCEL LIFTED UP
2496 ! THE MOIST ADIABAT IS BUOYANT WITH RESPECT TO THE ENVIRONMENT.
2497 DO I=1,IM
2498 KLCL(I)=KM+1
2499 KBOT(I)=KM+1
2500 KTOP(I)=0
2501 ENDDO
2502 DO K=1,KM
2503 DO I=1,IM
2504 TCLD(I,K)=0.
2505 QCLD(I,K)=0.
2506 ENDDO
2507 ENDDO
2508 DO K=K1,KM
2509 DO I=1,IM
2510 IF(PRSLK(I,K).LE.SLKMA(I)) THEN
2511 KLCL(I)=MIN(KLCL(I),K)
2512 CALL STMA(THEMA(I),PRSLK(I,K),TMA,QMA)
2513 ! TMA=FTMA(THEMA(I),PRSLK(I,K),QMA)
2514 TVCLD=TMA*(1.+FV*QMA)
2515 TVENV=TENV(I,K)*(1.+FV*QENV(I,K))
2516 IF(TVCLD.GT.TVENV) THEN
2517 KBOT(I)=MIN(KBOT(I),K)
2518 KTOP(I)=MAX(KTOP(I),K)
2519 TCLD(I,K)=TMA-TENV(I,K)
2520 QCLD(I,K)=QMA-QENV(I,K)
2521 ENDIF
2522 ENDIF
2523 ENDDO
2524 ENDDO
2525 !-----------------------------------------------------------------------
2526 RETURN
2527 END SUBROUTINE MSTADBT3
2529 #if (EM_CORE == 1)
2530 ! random seeds - ZCX
2531 SUBROUTINE init_random_seed()
2532 INTEGER :: i, n, clock
2533 INTEGER, DIMENSION(:), ALLOCATABLE :: seed
2535 CALL RANDOM_SEED(size = n)
2536 ALLOCATE(seed(n))
2538 CALL SYSTEM_CLOCK(COUNT=clock)
2540 seed = clock + 37 * (/ (i - 1, i = 1, n) /)
2541 CALL RANDOM_SEED(PUT = seed)
2543 DEALLOCATE(seed)
2544 END SUBROUTINE
2545 #endif
2546 END MODULE module_cu_sas