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