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