| blob | 89edb29a8fb51e9e5cdb9d2ed24160475fedb4d1 |
1 !!
2 MODULE module_cu_sas
4 CONTAINS
6 !-----------------------------------------------------------------
7 SUBROUTINE CU_SAS(DT,ITIMESTEP,STEPCU, &
8 RTHCUTEN,RQVCUTEN,RQCCUTEN,RQICUTEN, &
9 RUCUTEN,RVCUTEN, &
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 MOMMIX, & ! gopal's doing
15 PGCON,sas_mass_flux, &
16 shalconv,shal_pgcon, &
17 HPBL2D,EVAP2D,HEAT2D, & !Kwon for shallow convection
18 P_QI,P_FIRST_SCALAR, &
19 CUDT, CURR_SECS, ADAPT_STEP_FLAG, &
20 CUDTACTTIME, &
21 ids,ide, jds,jde, kds,kde, &
22 ims,ime, jms,jme, kms,kme, &
23 its,ite, jts,jte, kts,kte )
25 !-------------------------------------------------------------------
26 USE MODULE_GFS_MACHINE , ONLY : kind_phys
27 USE MODULE_GFS_FUNCPHYS , ONLY : gfuncphys
28 USE MODULE_GFS_PHYSCONS, grav => con_g, CP => con_CP, HVAP => con_HVAP &
29 &, RV => con_RV, FV => con_fvirt, T0C => con_T0C &
30 &, CVAP => con_CVAP, CLIQ => con_CLIQ &
31 &, EPS => con_eps, EPSM1 => con_epsm1 &
32 &, ROVCP => con_rocp, RD => con_rd
33 !-------------------------------------------------------------------
34 IMPLICIT NONE
35 !-------------------------------------------------------------------
36 !-- U3D 3D u-velocity interpolated to theta points (m/s)
37 !-- V3D 3D v-velocity interpolated to theta points (m/s)
38 !-- TH3D 3D potential temperature (K)
39 !-- T3D temperature (K)
40 !-- QV3D 3D water vapor mixing ratio (Kg/Kg)
41 !-- QC3D 3D cloud mixing ratio (Kg/Kg)
42 !-- QI3D 3D ice mixing ratio (Kg/Kg)
43 !-- P8w 3D pressure at full levels (Pa)
44 !-- Pcps 3D pressure (Pa)
45 !-- PI3D 3D exner function (dimensionless)
46 !-- rr3D 3D dry air density (kg/m^3)
47 !-- RUBLTEN U tendency due to
48 ! PBL parameterization (m/s^2)
49 !-- RVBLTEN V tendency due to
50 ! PBL parameterization (m/s^2)
51 !-- RTHBLTEN Theta tendency due to
52 ! PBL parameterization (K/s)
53 !-- RQVBLTEN Qv tendency due to
54 ! PBL parameterization (kg/kg/s)
55 !-- RQCBLTEN Qc tendency due to
56 ! PBL parameterization (kg/kg/s)
57 !-- RQIBLTEN Qi tendency due to
58 ! PBL parameterization (kg/kg/s)
59 !
60 !-- MOMMIX MOMENTUM MIXING COEFFICIENT (can be set in the namelist)
61 !-- RUCUTEN U tendency due to Cumulus Momentum Mixing (gopal's doing for SAS)
62 !-- RVCUTEN V tendency due to Cumulus Momentum Mixing (gopal's doing for SAS)
63 !
64 !-- CP heat capacity at constant pressure for dry air (J/kg/K)
65 !-- GRAV acceleration due to gravity (m/s^2)
66 !-- ROVCP R/CP
67 !-- RD gas constant for dry air (J/kg/K)
68 !-- ROVG R/G
69 !-- P_QI species index for cloud ice
70 !-- dz8w dz between full levels (m)
71 !-- z height above sea level (m)
72 !-- PSFC pressure at the surface (Pa)
73 !-- UST u* in similarity theory (m/s)
74 !-- PBL PBL height (m)
75 !-- PSIM similarity stability function for momentum
76 !-- PSIH similarity stability function for heat
77 !-- HFX upward heat flux at the surface (W/m^2)
78 !-- QFX upward moisture flux at the surface (kg/m^2/s)
79 !-- TSK surface temperature (K)
80 !-- GZ1OZ0 log(z/z0) where z0 is roughness length
81 !-- WSPD wind speed at lowest model level (m/s)
82 !-- BR bulk Richardson number in surface layer
83 !-- DT time step (s)
84 !-- rvovrd R_v divided by R_d (dimensionless)
85 !-- EP1 constant for virtual temperature (R_v/R_d - 1) (dimensionless)
86 !-- KARMAN Von Karman constant
87 !-- ids start index for i in domain
88 !-- ide end index for i in domain
89 !-- jds start index for j in domain
90 !-- jde end index for j in domain
91 !-- kds start index for k in domain
92 !-- kde end index for k in domain
93 !-- ims start index for i in memory
94 !-- ime end index for i in memory
95 !-- jms start index for j in memory
96 !-- jme end index for j in memory
97 !-- kms start index for k in memory
98 !-- kme end index for k in memory
99 !-- its start index for i in tile
100 !-- ite end index for i in tile
101 !-- jts start index for j in tile
102 !-- jte end index for j in tile
103 !-- kts start index for k in tile
104 !-- kte end index for k in tile
105 !-------------------------------------------------------------------
107 INTEGER :: ICLDCK
109 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
110 ims,ime, jms,jme, kms,kme, &
111 its,ite, jts,jte, kts,kte, &
112 ITIMESTEP, & !NSTD
113 P_FIRST_SCALAR, &
114 P_QC, &
115 P_QI, &
116 STEPCU
118 REAL, INTENT(IN) :: &
119 DT
121 REAL, OPTIONAL, INTENT(IN) :: PGCON,sas_mass_flux,shal_pgcon
122 INTEGER, OPTIONAL, INTENT(IN) :: shalconv
123 REAL(kind=kind_phys) :: PGCON_USE,SHAL_PGCON_USE,massf
124 INTEGER :: shalconv_use
125 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: &
126 RQCCUTEN, &
127 RQICUTEN, &
128 RQVCUTEN, &
129 RTHCUTEN
130 REAL, DIMENSION(ims:ime, jms:jme, kms:kme), INTENT(INOUT) :: &
131 RUCUTEN, &
132 RVCUTEN
133 REAL, OPTIONAL, INTENT(IN) :: MOMMIX
135 REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
136 INTENT(IN) :: HPBL2D,EVAP2D,HEAT2D !Kwon for sha
138 REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: &
139 XLAND
141 REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: &
142 RAINCV, PRATEC
144 REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: &
145 HBOT, &
146 HTOP
148 LOGICAL, DIMENSION(IMS:IME,JMS:JME), INTENT(INOUT) :: &
149 CU_ACT_FLAG
152 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: &
153 DZ8W, &
154 P8w, &
155 Pcps, &
156 PI3D, &
157 QC3D, &
158 QI3D, &
159 QV3D, &
160 RHO3D, &
161 T3D, &
162 U3D, &
163 V3D, &
164 W
166 ! Adaptive time-step variables
167 REAL, INTENT(IN ) :: CUDT
168 REAL, INTENT(IN ) :: CURR_SECS
169 LOGICAL,INTENT(IN ) , OPTIONAL :: ADAPT_STEP_FLAG
170 REAL, INTENT (INOUT) :: CUDTACTTIME
172 !--------------------------- LOCAL VARS ------------------------------
174 REAL, DIMENSION(ims:ime, jms:jme) :: &
175 PSFC
178 REAL (kind=kind_phys) :: &
179 DELT, &
180 DPSHC, &
181 RDELT, &
182 RSEED
184 REAL (kind=kind_phys), DIMENSION(its:ite) :: &
185 CLDWRK, &
186 PS, &
187 RCS, &
188 RN, &
189 SLIMSK, &
190 HPBL,EVAP,HEAT !Kwon for shallow convection
192 REAL (kind=kind_phys), DIMENSION(its:ite, kts:kte+1) :: &
193 PRSI
195 REAL (kind=kind_phys), DIMENSION(its:ite, kts:kte) :: &
196 DEL, &
197 DOT, &
198 PHIL, &
199 PRSL, &
200 PRSLK, &
201 Q1, &
202 T1, &
203 U1, &
204 V1, &
205 ZI, &
206 ZL
208 REAL (kind=kind_phys), DIMENSION(its:ite, kts:kte, 2) :: &
209 QL
211 INTEGER, DIMENSION(its:ite) :: &
212 KBOT, &
213 KTOP, &
214 KCNV
216 INTEGER :: &
217 I, &
218 IGPVS, &
219 IM, &
220 J, &
221 JCAP, &
222 K, &
223 KM, &
224 KP, &
225 KX, &
226 NCLOUD
228 LOGICAL :: run_param , doing_adapt_dt , decided
230 DATA IGPVS/0/
232 !-----------------------------------------------------------------------
233 !
234 !*** CHECK TO SEE IF THIS IS A CONVECTION TIMESTEP
235 !
237 ! Initialization for adaptive time step.
239 doing_adapt_dt = .FALSE.
240 IF ( PRESENT(adapt_step_flag) ) THEN
241 IF ( adapt_step_flag ) THEN
242 doing_adapt_dt = .TRUE.
243 IF ( cudtacttime .EQ. 0. ) THEN
244 cudtacttime = curr_secs + cudt*60.
245 END IF
246 END IF
247 END IF
249 ! Do we run through this scheme or not?
251 ! Test 1: If this is the initial model time, then yes.
252 ! ITIMESTEP=1
253 ! Test 2: If the user asked for the cumulus to be run every time step, then yes.
254 ! CUDT=0 or STEPCU=1
255 ! Test 3: If not adaptive dt, and this is on the requested cumulus frequency, then yes.
256 ! MOD(ITIMESTEP,STEPCU)=0
257 ! Test 4: If using adaptive dt and the current time is past the last requested activate cumulus time, then yes.
258 ! CURR_SECS >= CUDTACTTIME
260 ! If we do run through the scheme, we set the flag run_param to TRUE and we set the decided flag
261 ! to TRUE. The decided flag says that one of these tests was able to say "yes", run the scheme.
262 ! We only proceed to other tests if the previous tests all have left decided as FALSE.
264 ! If we set run_param to TRUE and this is adaptive time stepping, we set the time to the next
265 ! cumulus run.
267 decided = .FALSE.
268 run_param = .FALSE.
269 IF ( ( .NOT. decided ) .AND. &
270 ( itimestep .EQ. 1 ) ) THEN
271 run_param = .TRUE.
272 decided = .TRUE.
273 END IF
275 IF ( ( .NOT. decided ) .AND. &
276 ( ( cudt .EQ. 0. ) .OR. ( stepcu .EQ. 1 ) ) ) THEN
277 run_param = .TRUE.
278 decided = .TRUE.
279 END IF
281 IF ( ( .NOT. decided ) .AND. &
282 ( .NOT. doing_adapt_dt ) .AND. &
283 ( MOD(itimestep,stepcu) .EQ. 0 ) ) THEN
284 run_param = .TRUE.
285 decided = .TRUE.
286 END IF
288 IF ( ( .NOT. decided ) .AND. &
289 ( doing_adapt_dt ) .AND. &
290 ( curr_secs .GE. cudtacttime ) ) THEN
291 run_param = .TRUE.
292 decided = .TRUE.
293 cudtacttime = curr_secs + cudt*60
294 END IF
296 !-----------------------------------------------------------------------
298 if(present(shalconv)) then
299 shalconv_use=shalconv
300 else
301 #if (NMM_CORE==1)
302 shalconv_use=0
303 #else
304 #if (EM_CORE==1)
305 shalconv_use=1
306 #else
307 shalconv_use=0
308 #endif
309 #endif
310 endif
312 if(present(pgcon)) then
313 pgcon_use = pgcon
314 else
315 ! pgcon_use = 0.7 ! Gregory et al. (1997, QJRMS)
316 pgcon_use = 0.55 ! Zhang & Wu (2003,JAS), used in GFS (25km res spectral)
317 ! pgcon_use = 0.2 ! HWRF, for model tuning purposes
318 ! pgcon_use = 0.3 ! GFDL, or so I am told
320 ! For those attempting to tune pgcon:
322 ! The value of 0.55 comes from an observational study of
323 ! synoptic-scale deep convection and 0.7 came from an
324 ! incorrect fit to the same data. That value is likely
325 ! correct for deep convection at gridscales near that of GFS,
326 ! but is questionable in shallow convection, or for scales
327 ! much finer than synoptic scales.
329 ! Then again, the assumptions of SAS break down when the
330 ! gridscale is near the convection scale anyway. In a large
331 ! storm such as a hurricane, there is often no environment to
332 ! detrain into since adjancent gridsquares are also undergoing
333 ! active convection. Each gridsquare will no longer have many
334 ! updrafts and downdrafts. At sub-convective timescales, you
335 ! will find unstable columns for many (say, 5 second length)
336 ! timesteps in a real atmosphere during a convection cell's
337 ! lifetime, so forcing it to be neutrally stable is unphysical.
339 ! Hence, in scales near the convection scale (cells have
340 ! ~0.5-4km diameter in hurricanes), this parameter is more of a
341 ! tuning parameter to get a scheme that is inappropriate for
342 ! that resolution to do a reasonable job.
344 ! Your mileage might vary.
346 ! - Sam Trahan
347 endif
349 if(present(sas_mass_flux)) then
350 massf=sas_mass_flux
351 ! Use this to reduce the fluxes added by SAS to prevent
352 ! computational instability as a result of large fluxes.
353 else
354 massf=9e9 ! large number to disable check
355 endif
357 if(present(shal_pgcon)) then
358 if(shal_pgcon>=0) then
359 shal_pgcon_use = shal_pgcon
360 else
361 ! shal_pgcon<0 means use deep pgcon
362 shal_pgcon_use = pgcon_use
363 endif
364 else
365 ! Default: Same as deep convection pgcon
366 shal_pgcon_use = pgcon_use
367 ! Read the warning above though. It may be advisable for
368 ! these to be different.
369 endif
371 if_run_param: IF(run_param) THEN
373 DO J=JTS,JTE
374 DO I=ITS,ITE
375 CU_ACT_FLAG(I,J)=.TRUE.
376 ENDDO
377 ENDDO
379 IM=ITE-ITS+1
380 KX=KTE-KTS+1
381 JCAP=126
382 DPSHC=30_kind_phys
383 DELT=DT*STEPCU
384 RDELT=1./DELT
385 NCLOUD=1
388 DO J=jms,jme
389 DO I=ims,ime
390 PSFC(i,j)=P8w(i,kms,j)
391 ENDDO
392 ENDDO
394 if(igpvs.eq.0) CALL GFUNCPHYS
395 igpvs=1
397 !------------- J LOOP (OUTER) --------------------------------------------------
399 big_outer_j_loop: DO J=jts,jte
401 ! --------------- compute zi and zl -----------------------------------------
402 DO i=its,ite
403 ZI(I,KTS)=0.0
404 ENDDO
406 DO k=kts+1,kte
407 KM=K-1
408 DO i=its,ite
409 ZI(I,K)=ZI(I,KM)+dz8w(i,km,j)
410 ENDDO
411 ENDDO
413 DO k=kts+1,kte
414 KM=K-1
415 DO i=its,ite
416 ZL(I,KM)=(ZI(I,K)+ZI(I,KM))*0.5
417 ENDDO
418 ENDDO
420 DO i=its,ite
421 ZL(I,KTE)=2.*ZI(I,KTE)-ZL(I,KTE-1)
422 ENDDO
424 ! --------------- end compute zi and zl -------------------------------------
426 DO i=its,ite
427 PS(i)=PSFC(i,j)*.001
428 RCS(i)=1.
429 SLIMSK(i)=ABS(XLAND(i,j)-2.)
430 ENDDO
432 #if (NMM_CORE == 1)
433 if(shalconv_use==1) then
434 DO i=its,ite
435 HPBL(I) = HPBL2D(I,J) !kwon for shallow convection
436 EVAP(I) = EVAP2D(I,J) !kwon for shallow convection
437 HEAT(I) = HEAT2D(I,J) !kwon for shallow convection
438 ENDDO
439 endif
440 #endif
442 DO i=its,ite
443 PRSI(i,kts)=PS(i)
444 ENDDO
446 DO k=kts,kte
447 kp=k+1
448 DO i=its,ite
449 PRSL(I,K)=Pcps(i,k,j)*.001
450 PHIL(I,K)=ZL(I,K)*GRAV
451 DOT(i,k)=-5.0E-4*GRAV*rho3d(i,k,j)*(w(i,k,j)+w(i,kp,j))
452 ENDDO
453 ENDDO
455 DO k=kts,kte
456 DO i=its,ite
457 DEL(i,k)=PRSL(i,k)*GRAV/RD*dz8w(i,k,j)/T3D(i,k,j)
458 U1(i,k)=U3D(i,k,j)
459 V1(i,k)=V3D(i,k,j)
460 Q1(i,k)=QV3D(i,k,j)/(1.+QV3D(i,k,j))
461 T1(i,k)=T3D(i,k,j)
462 QL(i,k,1)=QI3D(i,k,j)/(1.+QI3D(i,k,j))
463 QL(i,k,2)=QC3D(i,k,j)/(1.+QC3D(i,k,j))
464 PRSLK(I,K)=(PRSL(i,k)*.01)**ROVCP
465 ENDDO
466 ENDDO
468 DO k=kts+1,kte+1
469 km=k-1
470 DO i=its,ite
471 PRSI(i,k)=PRSI(i,km)-del(i,km)
472 ENDDO
473 ENDDO
475 CALL SASCNVN(IM,IM,KX,JCAP,DELT,DEL,PRSL,PS,PHIL, &
476 QL,Q1,T1,U1,V1,RCS,CLDWRK,RN,KBOT, &
477 KTOP,KCNV,SLIMSK,DOT,NCLOUD,PGCON_USE,massf)
479 if_shallow_conv: if(shalconv_use==1) then
480 #if (NMM_CORE == 1)
481 ! NMM calls the new shallow convection developed by J Han
482 ! (Added to WRF by Y.Kwon)
483 call shalcnv(im,im,kx,jcap,delt,del,prsl,ps,phil,ql, &
484 & q1,t1,u1,v1,rcs,rn,kbot,ktop,kcnv,slimsk, &
485 & dot,ncloud,hpbl,heat,evap,shal_pgcon_use)
486 #else
487 #if (EM_CORE == 1)
488 ! NOTE: ARW should be able to call the new shalcnv here, but
489 ! they need to add the three new variables, so I'm leaving the
490 ! old shallow convection call here - Sam Trahan
491 CALL OLD_ARW_SHALCV(IM,IM,KX,DELT,DEL,PRSI,PRSL,PRSLK,KCNV,Q1,T1,DPSHC)
492 #else
493 ! Shallow convection is untested for other cores.
494 #endif
495 #endif
496 endif if_shallow_conv
498 DO I=ITS,ITE
499 RAINCV(I,J)=RN(I)*1000./STEPCU
500 PRATEC(I,J)=RN(I)*1000./(STEPCU * DT)
501 HBOT(I,J)=KBOT(I)
502 HTOP(I,J)=KTOP(I)
503 ENDDO
505 DO K=KTS,KTE
506 DO I=ITS,ITE
507 RTHCUTEN(I,K,J)=(T1(I,K)-T3D(I,K,J))/PI3D(I,K,J)*RDELT
508 RQVCUTEN(I,K,J)=(Q1(I,K)/(1.-q1(i,k))-QV3D(I,K,J))*RDELT
509 ENDDO
510 ENDDO
512 !===============================================================================
513 ! ADD MOMENTUM MIXING TERM AS TENDENCIES. This is gopal's doing for SAS
514 ! MOMMIX is the reduction factor set to 0.7 by default. Because NMM has
515 ! divergence damping term, a reducion factor for cumulum mixing may be
516 ! required otherwise storms were too weak.
517 !===============================================================================
518 !
519 #if (NMM_CORE == 1)
520 DO K=KTS,KTE
521 DO I=ITS,ITE
522 ! RUCUTEN(I,J,K)=MOMMIX*(U1(I,K)-U3D(I,K,J))*RDELT
523 ! RVCUTEN(I,J,K)=MOMMIX*(V1(I,K)-V3D(I,K,J))*RDELT
524 RUCUTEN(I,J,K)=(U1(I,K)-U3D(I,K,J))*RDELT
525 RVCUTEN(I,J,K)=(V1(I,K)-V3D(I,K,J))*RDELT
526 ENDDO
527 ENDDO
528 #endif
531 IF(P_QC .ge. P_FIRST_SCALAR)THEN
532 DO K=KTS,KTE
533 DO I=ITS,ITE
534 RQCCUTEN(I,K,J)=(QL(I,K,2)/(1.-ql(i,k,2))-QC3D(I,K,J))*RDELT
535 ENDDO
536 ENDDO
537 ENDIF
539 IF(P_QI .ge. P_FIRST_SCALAR)THEN
540 DO K=KTS,KTE
541 DO I=ITS,ITE
542 RQICUTEN(I,K,J)=(QL(I,K,1)/(1.-ql(i,k,1))-QI3D(I,K,J))*RDELT
543 ENDDO
544 ENDDO
545 ENDIF
547 ENDDO big_outer_j_loop ! Outer most J loop
549 ENDIF if_run_param
551 END SUBROUTINE CU_SAS
553 !====================================================================
554 SUBROUTINE sasinit(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQICUTEN, &
555 RUCUTEN,RVCUTEN, &
556 RESTART,P_QC,P_QI,P_FIRST_SCALAR, &
557 allowed_to_read, &
558 ids, ide, jds, jde, kds, kde, &
559 ims, ime, jms, jme, kms, kme, &
560 its, ite, jts, jte, kts, kte )
561 !--------------------------------------------------------------------
562 IMPLICIT NONE
563 !--------------------------------------------------------------------
564 LOGICAL , INTENT(IN) :: allowed_to_read,restart
565 INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
566 ims, ime, jms, jme, kms, kme, &
567 its, ite, jts, jte, kts, kte
568 INTEGER , INTENT(IN) :: P_FIRST_SCALAR, P_QI, P_QC
570 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: &
571 RTHCUTEN, &
572 RQVCUTEN, &
573 RQCCUTEN, &
574 RQICUTEN
575 REAL, DIMENSION( ims:ime , jms:jme , kms:kme ) , INTENT(OUT) :: &
576 RUCUTEN, & ! gopal's doing for SAS
577 RVCUTEN
579 INTEGER :: i, j, k, itf, jtf, ktf
581 jtf=min0(jte,jde-1)
582 ktf=min0(kte,kde-1)
583 itf=min0(ite,ide-1)
585 #ifdef HWRF
586 !zhang's doing
587 IF(.not.restart .or. .not.allowed_to_read)THEN
588 !end of zhang's doing
589 #else
590 IF(.not.restart)THEN
591 #endif
592 DO j=jts,jtf
593 DO k=kts,ktf
594 DO i=its,itf
595 RTHCUTEN(i,k,j)=0.
596 RQVCUTEN(i,k,j)=0.
597 RUCUTEN(i,j,k)=0.
598 RVCUTEN(i,j,k)=0.
599 ENDDO
600 ENDDO
601 ENDDO
603 IF (P_QC .ge. P_FIRST_SCALAR) THEN
604 DO j=jts,jtf
605 DO k=kts,ktf
606 DO i=its,itf
607 RQCCUTEN(i,k,j)=0.
608 ENDDO
609 ENDDO
610 ENDDO
611 ENDIF
613 IF (P_QI .ge. P_FIRST_SCALAR) THEN
614 DO j=jts,jtf
615 DO k=kts,ktf
616 DO i=its,itf
617 RQICUTEN(i,k,j)=0.
618 ENDDO
619 ENDDO
620 ENDDO
621 ENDIF
622 ENDIF
624 END SUBROUTINE sasinit
626 ! ------------------------------------------------------------------------
628 SUBROUTINE SASCNV(IM,IX,KM,JCAP,DELT,DEL,PRSL,PS,PHIL,QL, &
629 ! SUBROUTINE SASCNV(IM,IX,KM,JCAP,DLT,DEL,PRSL,PHIL,QL, &
630 & Q1,T1,U1,V1,RCS,CLDWRK,RN,KBOT,KTOP,KUO,SLIMSK, &
631 & DOT,XKT2,ncloud)
632 ! for cloud water version
633 ! parameter(ncloud=0)
634 ! SUBROUTINE SASCNV(KM,JCAP,DELT,DEL,SL,SLK,PS,QL,
635 ! & Q1,T1,U1,V1,RCS,CLDWRK,RN,KBOT,KTOP,KUO,SLIMSK,
636 ! & DOT,xkt2,ncloud)
637 !
638 USE MODULE_GFS_MACHINE , ONLY : kind_phys
639 USE MODULE_GFS_FUNCPHYS ,ONLY : fpvs
640 USE MODULE_GFS_PHYSCONS, grav => con_g, CP => con_CP, HVAP => con_HVAP &
641 &, RV => con_RV, FV => con_fvirt, T0C => con_T0C &
642 &, CVAP => con_CVAP, CLIQ => con_CLIQ &
643 &, EPS => con_eps, EPSM1 => con_epsm1
645 implicit none
646 !
647 ! include 'constant.h'
648 !
649 integer IM, IX, KM, JCAP, ncloud, &
650 & KBOT(IM), KTOP(IM), KUO(IM), J
651 real(kind=kind_phys) DELT
652 real(kind=kind_phys) PS(IM), DEL(IX,KM), PRSL(IX,KM), &
653 ! real(kind=kind_phys) DEL(IX,KM), PRSL(IX,KM),
654 & QL(IX,KM,2), Q1(IX,KM), T1(IX,KM), &
655 & U1(IX,KM), V1(IX,KM), RCS(IM), &
656 & CLDWRK(IM), RN(IM), SLIMSK(IM), &
657 & DOT(IX,KM), XKT2(IM), PHIL(IX,KM)
658 !
659 integer I, INDX, jmn, k, knumb, latd, lond, km1
660 !
661 real(kind=kind_phys) adw, alpha, alphal, alphas, &
662 & aup, beta, betal, betas, &
663 & c0, cpoel, dellat, delta, &
664 & desdt, deta, detad, dg, &
665 & dh, dhh, dlnsig, dp, &
666 & dq, dqsdp, dqsdt, dt, &
667 & dt2, dtmax, dtmin, dv1, &
668 & dv1q, dv2, dv2q, dv1u, &
669 & dv1v, dv2u, dv2v, dv3u, &
670 & dv3v, dv3, dv3q, dvq1, &
671 & dz, dz1, e1, edtmax, &
672 & edtmaxl, edtmaxs, el2orc, elocp, &
673 & es, etah, &
674 & evef, evfact, evfactl, fact1, &
675 & fact2, factor, fjcap, fkm, &
676 & fuv, g, gamma, onemf, &
677 & onemfu, pdetrn, pdpdwn, pprime, &
678 & qc, qlk, qrch, qs, &
679 & rain, rfact, shear, tem1, &
680 & tem2, terr, val, val1, &
681 & val2, w1, w1l, w1s, &
682 & w2, w2l, w2s, w3, &
683 & w3l, w3s, w4, w4l, &
684 & w4s, xdby, xpw, xpwd, &
685 & xqc, xqrch, xlambu, mbdt, &
686 & tem
687 !
688 !
689 integer JMIN(IM), KB(IM), KBCON(IM), KBDTR(IM), &
690 & KT2(IM), KTCON(IM), LMIN(IM), &
691 & kbm(IM), kbmax(IM), kmax(IM)
692 !
693 real(kind=kind_phys) AA1(IM), ACRT(IM), ACRTFCT(IM), &
694 & DELHBAR(IM), DELQ(IM), DELQ2(IM), &
695 & DELQBAR(IM), DELQEV(IM), DELTBAR(IM), &
696 & DELTV(IM), DTCONV(IM), EDT(IM), &
697 & EDTO(IM), EDTX(IM), FLD(IM), &
698 & HCDO(IM), HKBO(IM), HMAX(IM), &
699 & HMIN(IM), HSBAR(IM), UCDO(IM), &
700 & UKBO(IM), VCDO(IM), VKBO(IM), &
701 & PBCDIF(IM), PDOT(IM), PO(IM,KM), &
702 & PWAVO(IM), PWEVO(IM), &
703 ! & PSFC(IM), PWAVO(IM), PWEVO(IM), &
704 & QCDO(IM), QCOND(IM), QEVAP(IM), &
705 & QKBO(IM), RNTOT(IM), VSHEAR(IM), &
706 & XAA0(IM), XHCD(IM), XHKB(IM), &
707 & XK(IM), XLAMB(IM), XLAMD(IM), &
708 & XMB(IM), XMBMAX(IM), XPWAV(IM), &
709 & XPWEV(IM), XQCD(IM), XQKB(IM)
710 !
711 ! PHYSICAL PARAMETERS
712 PARAMETER(G=grav)
713 PARAMETER(CPOEL=CP/HVAP,ELOCP=HVAP/CP, &
714 & EL2ORC=HVAP*HVAP/(RV*CP))
715 PARAMETER(TERR=0.,C0=.002,DELTA=fv)
716 PARAMETER(FACT1=(CVAP-CLIQ)/RV,FACT2=HVAP/RV-FACT1*T0C)
717 ! LOCAL VARIABLES AND ARRAYS
718 real(kind=kind_phys) PFLD(IM,KM), TO(IM,KM), QO(IM,KM), &
719 & UO(IM,KM), VO(IM,KM), QESO(IM,KM)
720 ! cloud water
721 real(kind=kind_phys) QLKO_KTCON(IM), DELLAL(IM), TVO(IM,KM), &
722 & DBYO(IM,KM), ZO(IM,KM), SUMZ(IM,KM), &
723 & SUMH(IM,KM), HEO(IM,KM), HESO(IM,KM), &
724 & QRCD(IM,KM), DELLAH(IM,KM), DELLAQ(IM,KM),&
725 & DELLAU(IM,KM), DELLAV(IM,KM), HCKO(IM,KM), &
726 & UCKO(IM,KM), VCKO(IM,KM), QCKO(IM,KM), &
727 & ETA(IM,KM), ETAU(IM,KM), ETAD(IM,KM), &
728 & QRCDO(IM,KM), PWO(IM,KM), PWDO(IM,KM), &
729 & RHBAR(IM), TX1(IM)
730 !
731 LOGICAL TOTFLG, CNVFLG(IM), DWNFLG(IM), DWNFLG2(IM), FLG(IM)
732 !
733 real(kind=kind_phys) PCRIT(15), ACRITT(15), ACRIT(15)
734 ! SAVE PCRIT, ACRITT
735 DATA PCRIT/850.,800.,750.,700.,650.,600.,550.,500.,450.,400., &
736 & 350.,300.,250.,200.,150./
737 DATA ACRITT/.0633,.0445,.0553,.0664,.075,.1082,.1521,.2216, &
738 & .3151,.3677,.41,.5255,.7663,1.1686,1.6851/
739 ! GDAS DERIVED ACRIT
740 ! DATA ACRITT/.203,.515,.521,.566,.625,.665,.659,.688, &
741 ! & .743,.813,.886,.947,1.138,1.377,1.896/
742 !
743 real(kind=kind_phys) TF, TCR, TCRF, RZERO, RONE
744 parameter (TF=233.16, TCR=263.16, TCRF=1.0/(TCR-TF))
745 parameter (RZERO=0.0,RONE=1.0)
746 !-----------------------------------------------------------------------
747 !
748 km1 = km - 1
749 ! INITIALIZE ARRAYS
750 !
751 DO I=1,IM
752 RN(I)=0.
753 KBOT(I)=KM+1
754 KTOP(I)=0
755 KUO(I)=0
756 CNVFLG(I) = .TRUE.
757 DTCONV(I) = 3600.
758 CLDWRK(I) = 0.
759 PDOT(I) = 0.
760 KT2(I) = 0
761 QLKO_KTCON(I) = 0.
762 DELLAL(I) = 0.
763 ENDDO
764 !!
765 DO K = 1, 15
766 ACRIT(K) = ACRITT(K) * (975. - PCRIT(K))
767 ENDDO
768 DT2 = DELT
769 !cmr dtmin = max(dt2,1200.)
770 val = 1200.
771 dtmin = max(dt2, val )
772 !cmr dtmax = max(dt2,3600.)
773 val = 3600.
774 dtmax = max(dt2, val )
775 ! MODEL TUNABLE PARAMETERS ARE ALL HERE
776 MBDT = 10.
777 EDTMAXl = .3
778 EDTMAXs = .3
779 ALPHAl = .5
780 ALPHAs = .5
781 BETAl = .15
782 betas = .15
783 BETAl = .05
784 betas = .05
785 ! change for hurricane model
786 BETAl = .5
787 betas = .5
788 ! EVEF = 0.07
789 evfact = 0.3
790 evfactl = 0.3
791 ! change for hurricane model
792 evfact = 0.6
793 evfactl = .6
794 #if ( EM_CORE == 1 )
795 ! HAWAII TEST - ZCX
796 ALPHAl = .5
797 ALPHAs = .75
798 BETAl = .05
799 betas = .05
800 evfact = 0.5
801 evfactl = 0.5
802 #endif
803 PDPDWN = 0.
804 PDETRN = 200.
805 xlambu = 1.e-4
806 fjcap = (float(jcap) / 126.) ** 2
807 !cmr fjcap = max(fjcap,1.)
808 val = 1.
809 fjcap = max(fjcap,val)
810 fkm = (float(km) / 28.) ** 2
811 !cmr fkm = max(fkm,1.)
812 fkm = max(fkm,val)
813 W1l = -8.E-3
814 W2l = -4.E-2
815 W3l = -5.E-3
816 W4l = -5.E-4
817 W1s = -2.E-4
818 W2s = -2.E-3
819 W3s = -1.E-3
820 W4s = -2.E-5
821 !CCCC IF(IM.EQ.384) THEN
822 LATD = 92
823 lond = 189
824 !CCCC ELSEIF(IM.EQ.768) THEN
825 !CCCC LATD = 80
826 !CCCC ELSE
827 !CCCC LATD = 0
828 !CCCC ENDIF
829 !
830 ! DEFINE TOP LAYER FOR SEARCH OF THE DOWNDRAFT ORIGINATING LAYER
831 ! AND THE MAXIMUM THETAE FOR UPDRAFT
832 !
833 DO I=1,IM
834 KBMAX(I) = KM
835 KBM(I) = KM
836 KMAX(I) = KM
837 TX1(I) = 1.0 / PS(I)
838 ENDDO
839 !
840 DO K = 1, KM
841 DO I=1,IM
842 IF (prSL(I,K)*tx1(I) .GT. 0.45) KBMAX(I) = K + 1
843 IF (prSL(I,K)*tx1(I) .GT. 0.70) KBM(I) = K + 1
844 IF (prSL(I,K)*tx1(I) .GT. 0.04) KMAX(I) = MIN(KM,K + 1)
845 ENDDO
846 ENDDO
847 DO I=1,IM
848 KBMAX(I) = MIN(KBMAX(I),KMAX(I))
849 KBM(I) = MIN(KBM(I),KMAX(I))
850 ENDDO
851 !
852 ! CONVERT SURFACE PRESSURE TO MB FROM CB
853 !
854 !!
855 DO K = 1, KM
856 DO I=1,IM
857 if (K .le. kmax(i)) then
858 PFLD(I,k) = PRSL(I,K) * 10.0
859 PWO(I,k) = 0.
860 PWDO(I,k) = 0.
861 TO(I,k) = T1(I,k)
862 QO(I,k) = Q1(I,k)
863 UO(I,k) = U1(I,k)
864 VO(I,k) = V1(I,k)
865 DBYO(I,k) = 0.
866 SUMZ(I,k) = 0.
867 SUMH(I,k) = 0.
868 endif
869 ENDDO
870 ENDDO
872 !
873 ! COLUMN VARIABLES
874 ! P IS PRESSURE OF THE LAYER (MB)
875 ! T IS TEMPERATURE AT T-DT (K)..TN
876 ! Q IS MIXING RATIO AT T-DT (KG/KG)..QN
877 ! TO IS TEMPERATURE AT T+DT (K)... THIS IS AFTER ADVECTION AND TURBULAN
878 ! QO IS MIXING RATIO AT T+DT (KG/KG)..Q1
879 !
880 DO K = 1, KM
881 DO I=1,IM
882 if (k .le. kmax(i)) then
883 !jfe QESO(I,k) = 10. * FPVS(T1(I,k))
884 !
885 QESO(I,k) = 0.01 * fpvs(T1(I,K)) ! fpvs is in Pa
886 !
887 QESO(I,k) = EPS * QESO(I,k) / (PFLD(I,k) + EPSM1*QESO(I,k))
888 !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
889 val1 = 1.E-8
890 QESO(I,k) = MAX(QESO(I,k), val1)
891 !cmr QO(I,k) = max(QO(I,k),1.e-10)
892 val2 = 1.e-10
893 QO(I,k) = max(QO(I,k), val2 )
894 ! QO(I,k) = MIN(QO(I,k),QESO(I,k))
895 TVO(I,k) = TO(I,k) + DELTA * TO(I,k) * QO(I,k)
896 endif
897 ENDDO
898 ENDDO
900 !
901 ! HYDROSTATIC HEIGHT ASSUME ZERO TERR
902 !
903 DO K = 1, KM
904 DO I=1,IM
905 ZO(I,k) = PHIL(I,k) / G
906 ENDDO
907 ENDDO
908 ! COMPUTE MOIST STATIC ENERGY
909 DO K = 1, KM
910 DO I=1,IM
911 if (K .le. kmax(i)) then
912 ! tem = G * ZO(I,k) + CP * TO(I,k)
913 tem = PHIL(I,k) + CP * TO(I,k)
914 HEO(I,k) = tem + HVAP * QO(I,k)
915 HESO(I,k) = tem + HVAP * QESO(I,k)
916 ! HEO(I,k) = MIN(HEO(I,k),HESO(I,k))
917 endif
918 ENDDO
919 ENDDO
920 !
921 ! DETERMINE LEVEL WITH LARGEST MOIST STATIC ENERGY
922 ! THIS IS THE LEVEL WHERE UPDRAFT STARTS
923 !
924 DO I=1,IM
925 HMAX(I) = HEO(I,1)
926 KB(I) = 1
927 ENDDO
928 !!
929 DO K = 2, KM
930 DO I=1,IM
931 if (k .le. kbm(i)) then
932 IF(HEO(I,k).GT.HMAX(I).AND.CNVFLG(I)) THEN
933 KB(I) = K
934 HMAX(I) = HEO(I,k)
935 ENDIF
936 endif
937 ENDDO
938 ENDDO
939 ! DO K = 1, KMAX - 1
940 ! TOL(k) = .5 * (TO(I,k) + TO(I,k+1))
941 ! QOL(k) = .5 * (QO(I,k) + QO(I,k+1))
942 ! QESOL(I,k) = .5 * (QESO(I,k) + QESO(I,k+1))
943 ! HEOL(I,k) = .5 * (HEO(I,k) + HEO(I,k+1))
944 ! HESOL(I,k) = .5 * (HESO(I,k) + HESO(I,k+1))
945 ! ENDDO
946 DO K = 1, KM1
947 DO I=1,IM
948 if (k .le. kmax(i)-1) then
949 DZ = .5 * (ZO(I,k+1) - ZO(I,k))
950 DP = .5 * (PFLD(I,k+1) - PFLD(I,k))
951 !jfe ES = 10. * FPVS(TO(I,k+1))
952 !
953 ES = 0.01 * fpvs(TO(I,K+1)) ! fpvs is in Pa
954 !
955 PPRIME = PFLD(I,k+1) + EPSM1 * ES
956 QS = EPS * ES / PPRIME
957 DQSDP = - QS / PPRIME
958 DESDT = ES * (FACT1 / TO(I,k+1) + FACT2 / (TO(I,k+1)**2))
959 DQSDT = QS * PFLD(I,k+1) * DESDT / (ES * PPRIME)
960 GAMMA = EL2ORC * QESO(I,k+1) / (TO(I,k+1)**2)
961 DT = (G * DZ + HVAP * DQSDP * DP) / (CP * (1. + GAMMA))
962 DQ = DQSDT * DT + DQSDP * DP
963 TO(I,k) = TO(I,k+1) + DT
964 QO(I,k) = QO(I,k+1) + DQ
965 PO(I,k) = .5 * (PFLD(I,k) + PFLD(I,k+1))
966 endif
967 ENDDO
968 ENDDO
969 !
970 DO K = 1, KM1
971 DO I=1,IM
972 if (k .le. kmax(I)-1) then
973 !jfe QESO(I,k) = 10. * FPVS(TO(I,k))
974 !
975 QESO(I,k) = 0.01 * fpvs(TO(I,K)) ! fpvs is in Pa
976 !
977 QESO(I,k) = EPS * QESO(I,k) / (PO(I,k) + EPSM1*QESO(I,k))
978 !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
979 val1 = 1.E-8
980 QESO(I,k) = MAX(QESO(I,k), val1)
981 !cmr QO(I,k) = max(QO(I,k),1.e-10)
982 val2 = 1.e-10
983 QO(I,k) = max(QO(I,k), val2 )
984 ! QO(I,k) = MIN(QO(I,k),QESO(I,k))
985 HEO(I,k) = .5 * G * (ZO(I,k) + ZO(I,k+1)) + &
986 & CP * TO(I,k) + HVAP * QO(I,k)
987 HESO(I,k) = .5 * G * (ZO(I,k) + ZO(I,k+1)) + &
988 & CP * TO(I,k) + HVAP * QESO(I,k)
989 UO(I,k) = .5 * (UO(I,k) + UO(I,k+1))
990 VO(I,k) = .5 * (VO(I,k) + VO(I,k+1))
991 endif
992 ENDDO
993 ENDDO
994 ! k = kmax
995 ! HEO(I,k) = HEO(I,k)
996 ! hesol(k) = HESO(I,k)
997 ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
998 ! PRINT *, ' HEO ='
999 ! PRINT 6001, (HEO(I,K),K=1,KMAX)
1000 ! PRINT *, ' HESO ='
1001 ! PRINT 6001, (HESO(I,K),K=1,KMAX)
1002 ! PRINT *, ' TO ='
1003 ! PRINT 6002, (TO(I,K)-273.16,K=1,KMAX)
1004 ! PRINT *, ' QO ='
1005 ! PRINT 6003, (QO(I,K),K=1,KMAX)
1006 ! PRINT *, ' QSO ='
1007 ! PRINT 6003, (QESO(I,K),K=1,KMAX)
1008 ! ENDIF
1009 !
1010 ! LOOK FOR CONVECTIVE CLOUD BASE AS THE LEVEL OF FREE CONVECTION
1011 !
1012 DO I=1,IM
1013 IF(CNVFLG(I)) THEN
1014 INDX = KB(I)
1015 HKBO(I) = HEO(I,INDX)
1016 QKBO(I) = QO(I,INDX)
1017 UKBO(I) = UO(I,INDX)
1018 VKBO(I) = VO(I,INDX)
1019 ENDIF
1020 FLG(I) = CNVFLG(I)
1021 KBCON(I) = KMAX(I)
1022 ENDDO
1023 !!
1024 DO K = 1, KM
1025 DO I=1,IM
1026 if (k .le. kbmax(i)) then
1027 IF(FLG(I).AND.K.GT.KB(I)) THEN
1028 HSBAR(I) = HESO(I,k)
1029 IF(HKBO(I).GT.HSBAR(I)) THEN
1030 FLG(I) = .FALSE.
1031 KBCON(I) = K
1032 ENDIF
1033 ENDIF
1034 endif
1035 ENDDO
1036 ENDDO
1037 DO I=1,IM
1038 IF(CNVFLG(I)) THEN
1039 PBCDIF(I) = -PFLD(I,KBCON(I)) + PFLD(I,KB(I))
1040 PDOT(I) = 10.* DOT(I,KBCON(I))
1041 IF(PBCDIF(I).GT.150.) CNVFLG(I) = .FALSE.
1042 IF(KBCON(I).EQ.KMAX(I)) CNVFLG(I) = .FALSE.
1043 ENDIF
1044 ENDDO
1045 !!
1046 TOTFLG = .TRUE.
1047 DO I=1,IM
1048 TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
1049 ENDDO
1050 IF(TOTFLG) RETURN
1051 ! FOUND LFC, CAN DEFINE REST OF VARIABLES
1052 6001 FORMAT(2X,-2P10F12.2)
1053 6002 FORMAT(2X,10F12.2)
1054 6003 FORMAT(2X,3P10F12.2)
1056 !
1057 ! DETERMINE ENTRAINMENT RATE BETWEEN KB AND KBCON
1058 !
1059 DO I = 1, IM
1060 alpha = alphas
1061 if(SLIMSK(I).eq.1.) alpha = alphal
1062 IF(CNVFLG(I)) THEN
1063 IF(KB(I).EQ.1) THEN
1064 DZ = .5 * (ZO(I,KBCON(I)) + ZO(I,KBCON(I)-1)) - ZO(I,1)
1065 ELSE
1066 DZ = .5 * (ZO(I,KBCON(I)) + ZO(I,KBCON(I)-1)) &
1067 & - .5 * (ZO(I,KB(I)) + ZO(I,KB(I)-1))
1068 ENDIF
1069 IF(KBCON(I).NE.KB(I)) THEN
1070 !cmr XLAMB(I) = -ALOG(ALPHA) / DZ
1071 XLAMB(I) = - LOG(ALPHA) / DZ
1072 ELSE
1073 XLAMB(I) = 0.
1074 ENDIF
1075 ENDIF
1076 ENDDO
1077 ! DETERMINE UPDRAFT MASS FLUX
1078 DO K = 1, KM
1079 DO I = 1, IM
1080 if (k .le. kmax(i) .and. CNVFLG(I)) then
1081 ETA(I,k) = 1.
1082 ETAU(I,k) = 1.
1083 ENDIF
1084 ENDDO
1085 ENDDO
1086 DO K = KM1, 2, -1
1087 DO I = 1, IM
1088 if (k .le. kbmax(i)) then
1089 IF(CNVFLG(I).AND.K.LT.KBCON(I).AND.K.GE.KB(I)) THEN
1090 DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1091 ETA(I,k) = ETA(I,k+1) * EXP(-XLAMB(I) * DZ)
1092 ETAU(I,k) = ETA(I,k)
1093 ENDIF
1094 endif
1095 ENDDO
1096 ENDDO
1097 DO I = 1, IM
1098 IF(CNVFLG(I).AND.KB(I).EQ.1.AND.KBCON(I).GT.1) THEN
1099 DZ = .5 * (ZO(I,2) - ZO(I,1))
1100 ETA(I,1) = ETA(I,2) * EXP(-XLAMB(I) * DZ)
1101 ETAU(I,1) = ETA(I,1)
1102 ENDIF
1103 ENDDO
1104 !
1105 ! WORK UP UPDRAFT CLOUD PROPERTIES
1106 !
1107 DO I = 1, IM
1108 IF(CNVFLG(I)) THEN
1109 INDX = KB(I)
1110 HCKO(I,INDX) = HKBO(I)
1111 QCKO(I,INDX) = QKBO(I)
1112 UCKO(I,INDX) = UKBO(I)
1113 VCKO(I,INDX) = VKBO(I)
1114 PWAVO(I) = 0.
1115 ENDIF
1116 ENDDO
1117 !
1118 ! CLOUD PROPERTY BELOW CLOUD BASE IS MODIFIED BY THE ENTRAINMENT PROCES
1119 !
1120 DO K = 2, KM1
1121 DO I = 1, IM
1122 if (k .le. kmax(i)-1) then
1123 IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LE.KBCON(I)) THEN
1124 FACTOR = ETA(I,k-1) / ETA(I,k)
1125 ONEMF = 1. - FACTOR
1126 HCKO(I,k) = FACTOR * HCKO(I,k-1) + ONEMF * &
1127 & .5 * (HEO(I,k) + HEO(I,k+1))
1128 UCKO(I,k) = FACTOR * UCKO(I,k-1) + ONEMF * &
1129 & .5 * (UO(I,k) + UO(I,k+1))
1130 VCKO(I,k) = FACTOR * VCKO(I,k-1) + ONEMF * &
1131 & .5 * (VO(I,k) + VO(I,k+1))
1132 DBYO(I,k) = HCKO(I,k) - HESO(I,k)
1133 ENDIF
1134 IF(CNVFLG(I).AND.K.GT.KBCON(I)) THEN
1135 HCKO(I,k) = HCKO(I,k-1)
1136 UCKO(I,k) = UCKO(I,k-1)
1137 VCKO(I,k) = VCKO(I,k-1)
1138 DBYO(I,k) = HCKO(I,k) - HESO(I,k)
1139 ENDIF
1140 endif
1141 ENDDO
1142 ENDDO
1143 ! DETERMINE CLOUD TOP
1144 DO I = 1, IM
1145 FLG(I) = CNVFLG(I)
1146 KTCON(I) = 1
1147 ENDDO
1148 ! DO K = 2, KMAX
1149 ! KK = KMAX - K + 1
1150 ! IF(DBYO(I,kK).GE.0..AND.FLG(I).AND.KK.GT.KBCON(I)) THEN
1151 ! KTCON(I) = KK + 1
1152 ! FLG(I) = .FALSE.
1153 ! ENDIF
1154 ! ENDDO
1155 DO K = 2, KM
1156 DO I = 1, IM
1157 if (k .le. kmax(i)) then
1158 IF(DBYO(I,k).LT.0..AND.FLG(I).AND.K.GT.KBCON(I)) THEN
1159 KTCON(I) = K
1160 FLG(I) = .FALSE.
1161 ENDIF
1162 endif
1163 ENDDO
1164 ENDDO
1165 DO I = 1, IM
1166 IF(CNVFLG(I).AND.(PFLD(I,KBCON(I)) - PFLD(I,KTCON(I))).LT.150.) &
1167 & CNVFLG(I) = .FALSE.
1168 ENDDO
1169 TOTFLG = .TRUE.
1170 DO I = 1, IM
1171 TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
1172 ENDDO
1173 IF(TOTFLG) RETURN
1174 !
1175 ! SEARCH FOR DOWNDRAFT ORIGINATING LEVEL ABOVE THETA-E MINIMUM
1176 !
1177 DO I = 1, IM
1178 HMIN(I) = HEO(I,KBCON(I))
1179 LMIN(I) = KBMAX(I)
1180 JMIN(I) = KBMAX(I)
1181 ENDDO
1182 DO I = 1, IM
1183 DO K = KBCON(I), KBMAX(I)
1184 IF(HEO(I,k).LT.HMIN(I).AND.CNVFLG(I)) THEN
1185 LMIN(I) = K + 1
1186 HMIN(I) = HEO(I,k)
1187 ENDIF
1188 ENDDO
1189 ENDDO
1190 !
1191 ! Make sure that JMIN(I) is within the cloud
1192 !
1193 DO I = 1, IM
1194 IF(CNVFLG(I)) THEN
1195 JMIN(I) = MIN(LMIN(I),KTCON(I)-1)
1196 XMBMAX(I) = .1
1197 JMIN(I) = MAX(JMIN(I),KBCON(I)+1)
1198 ENDIF
1199 ENDDO
1200 !
1201 ! ENTRAINING CLOUD
1202 !
1203 do k = 2, km1
1204 DO I = 1, IM
1205 if (k .le. kmax(i)-1) then
1206 if(CNVFLG(I).and.k.gt.JMIN(I).and.k.le.KTCON(I)) THEN
1207 SUMZ(I,k) = SUMZ(I,k-1) + .5 * (ZO(I,k+1) - ZO(I,k-1))
1208 SUMH(I,k) = SUMH(I,k-1) + .5 * (ZO(I,k+1) - ZO(I,k-1)) &
1209 & * HEO(I,k)
1210 ENDIF
1211 endif
1212 enddo
1213 enddo
1214 !!
1215 DO I = 1, IM
1216 IF(CNVFLG(I)) THEN
1217 ! call random_number(XKT2)
1218 ! call srand(fhour)
1219 ! XKT2(I) = rand()
1220 KT2(I) = nint(XKT2(I)*float(KTCON(I)-JMIN(I))-.5)+JMIN(I)+1
1221 ! KT2(I) = nint(sqrt(XKT2(I))*float(KTCON(I)-JMIN(I))-.5) + JMIN(I) + 1
1222 ! KT2(I) = nint(ranf() *float(KTCON(I)-JMIN(I))-.5) + JMIN(I) + 1
1223 tem1 = (HCKO(I,JMIN(I)) - HESO(I,KT2(I)))
1224 tem2 = (SUMZ(I,KT2(I)) * HESO(I,KT2(I)) - SUMH(I,KT2(I)))
1225 if (abs(tem2) .gt. 0.000001) THEN
1226 XLAMB(I) = tem1 / tem2
1227 else
1228 CNVFLG(I) = .false.
1229 ENDIF
1230 ! XLAMB(I) = (HCKO(I,JMIN(I)) - HESO(I,KT2(I)))
1231 ! & / (SUMZ(I,KT2(I)) * HESO(I,KT2(I)) - SUMH(I,KT2(I)))
1232 XLAMB(I) = max(XLAMB(I),RZERO)
1233 XLAMB(I) = min(XLAMB(I),2.3/SUMZ(I,KT2(I)))
1234 ENDIF
1235 ENDDO
1236 !!
1237 DO I = 1, IM
1238 DWNFLG(I) = CNVFLG(I)
1239 DWNFLG2(I) = CNVFLG(I)
1240 IF(CNVFLG(I)) THEN
1241 if(KT2(I).ge.KTCON(I)) DWNFLG(I) = .false.
1242 if(XLAMB(I).le.1.e-30.or.HCKO(I,JMIN(I))-HESO(I,KT2(I)).le.1.e-30)&
1243 & DWNFLG(I) = .false.
1244 do k = JMIN(I), KT2(I)
1245 if(DWNFLG(I).and.HEO(I,k).gt.HESO(I,KT2(I))) DWNFLG(I)=.false.
1246 enddo
1247 ! IF(CNVFLG(I).AND.(PFLD(KBCON(I))-PFLD(KTCON(I))).GT.PDETRN)
1248 ! & DWNFLG(I)=.FALSE.
1249 IF(CNVFLG(I).AND.(PFLD(I,KBCON(I))-PFLD(I,KTCON(I))).LT.PDPDWN) &
1250 & DWNFLG2(I)=.FALSE.
1251 ENDIF
1252 ENDDO
1253 !!
1254 DO K = 2, KM1
1255 DO I = 1, IM
1256 if (k .le. kmax(i)-1) then
1257 IF(DWNFLG(I).AND.K.GT.JMIN(I).AND.K.LE.KT2(I)) THEN
1258 DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1259 ! ETA(I,k) = ETA(I,k-1) * EXP( XLAMB(I) * DZ)
1260 ! to simplify matter, we will take the linear approach here
1261 !
1262 ETA(I,k) = ETA(I,k-1) * (1. + XLAMB(I) * dz)
1263 ETAU(I,k) = ETAU(I,k-1) * (1. + (XLAMB(I)+xlambu) * dz)
1264 ENDIF
1265 endif
1266 ENDDO
1267 ENDDO
1268 !!
1269 DO K = 2, KM1
1270 DO I = 1, IM
1271 if (k .le. kmax(i)-1) then
1272 ! IF(.NOT.DWNFLG(I).AND.K.GT.JMIN(I).AND.K.LE.KT2(I)) THEN
1273 IF(.NOT.DWNFLG(I).AND.K.GT.JMIN(I).AND.K.LE.KTCON(I)) THEN
1274 DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1275 ETAU(I,k) = ETAU(I,k-1) * (1. + xlambu * dz)
1276 ENDIF
1277 endif
1278 ENDDO
1279 ENDDO
1280 ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
1281 ! PRINT *, ' LMIN(I), KT2(I)=', LMIN(I), KT2(I)
1282 ! PRINT *, ' KBOT, KTOP, JMIN(I) =', KBCON(I), KTCON(I), JMIN(I)
1283 ! ENDIF
1284 ! IF(LAT.EQ.LATD.AND.lon.eq.lond) THEN
1285 ! print *, ' xlamb =', xlamb
1286 ! print *, ' eta =', (eta(k),k=1,KT2(I))
1287 ! print *, ' ETAU =', (ETAU(I,k),k=1,KT2(I))
1288 ! print *, ' HCKO =', (HCKO(I,k),k=1,KT2(I))
1289 ! print *, ' SUMZ =', (SUMZ(I,k),k=1,KT2(I))
1290 ! print *, ' SUMH =', (SUMH(I,k),k=1,KT2(I))
1291 ! ENDIF
1292 DO I = 1, IM
1293 if(DWNFLG(I)) THEN
1294 KTCON(I) = KT2(I)
1295 ENDIF
1296 ENDDO
1297 !
1298 ! CLOUD PROPERTY ABOVE CLOUD Base IS MODIFIED BY THE DETRAINMENT PROCESS
1299 !
1300 DO K = 2, KM1
1301 DO I = 1, IM
1302 if (k .le. kmax(i)-1) then
1303 !jfe
1304 IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
1305 !jfe IF(K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
1306 FACTOR = ETA(I,k-1) / ETA(I,k)
1307 ONEMF = 1. - FACTOR
1308 fuv = ETAU(I,k-1) / ETAU(I,k)
1309 onemfu = 1. - fuv
1310 HCKO(I,k) = FACTOR * HCKO(I,k-1) + ONEMF * &
1311 & .5 * (HEO(I,k) + HEO(I,k+1))
1312 UCKO(I,k) = fuv * UCKO(I,k-1) + ONEMFu * &
1313 & .5 * (UO(I,k) + UO(I,k+1))
1314 VCKO(I,k) = fuv * VCKO(I,k-1) + ONEMFu * &
1315 & .5 * (VO(I,k) + VO(I,k+1))
1316 DBYO(I,k) = HCKO(I,k) - HESO(I,k)
1317 ENDIF
1318 endif
1319 ENDDO
1320 ENDDO
1321 ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
1322 ! PRINT *, ' UCKO=', (UCKO(I,k),k=KBCON(I)+1,KTCON(I))
1323 ! PRINT *, ' uenv=', (.5*(UO(I,k)+UO(I,k-1)),k=KBCON(I)+1,KTCON(I))
1324 ! ENDIF
1325 DO I = 1, IM
1326 if(CNVFLG(I).and.DWNFLG2(I).and.JMIN(I).le.KBCON(I)) &
1327 & THEN
1328 CNVFLG(I) = .false.
1329 DWNFLG(I) = .false.
1330 DWNFLG2(I) = .false.
1331 ENDIF
1332 ENDDO
1333 !!
1334 TOTFLG = .TRUE.
1335 DO I = 1, IM
1336 TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
1337 ENDDO
1338 IF(TOTFLG) RETURN
1339 !!
1340 !
1341 ! COMPUTE CLOUD MOISTURE PROPERTY AND PRECIPITATION
1342 !
1343 DO I = 1, IM
1344 AA1(I) = 0.
1345 RHBAR(I) = 0.
1346 ENDDO
1347 DO K = 1, KM
1348 DO I = 1, IM
1349 if (k .le. kmax(i)) then
1350 IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LT.KTCON(I)) THEN
1351 DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1352 DZ1 = (ZO(I,k) - ZO(I,k-1))
1353 GAMMA = EL2ORC * QESO(I,k) / (TO(I,k)**2)
1354 QRCH = QESO(I,k) &
1355 & + GAMMA * DBYO(I,k) / (HVAP * (1. + GAMMA))
1356 FACTOR = ETA(I,k-1) / ETA(I,k)
1357 ONEMF = 1. - FACTOR
1358 QCKO(I,k) = FACTOR * QCKO(I,k-1) + ONEMF * &
1359 & .5 * (QO(I,k) + QO(I,k+1))
1360 DQ = ETA(I,k) * QCKO(I,k) - ETA(I,k) * QRCH
1361 RHBAR(I) = RHBAR(I) + QO(I,k) / QESO(I,k)
1362 !
1363 ! BELOW LFC CHECK IF THERE IS EXCESS MOISTURE TO RELEASE LATENT HEAT
1364 !
1365 IF(DQ.GT.0.) THEN
1366 ETAH = .5 * (ETA(I,k) + ETA(I,k-1))
1367 QLK = DQ / (ETA(I,k) + ETAH * C0 * DZ)
1368 AA1(I) = AA1(I) - DZ1 * G * QLK
1369 QC = QLK + QRCH
1370 PWO(I,k) = ETAH * C0 * DZ * QLK
1371 QCKO(I,k) = QC
1372 PWAVO(I) = PWAVO(I) + PWO(I,k)
1373 ENDIF
1374 ENDIF
1375 endif
1376 ENDDO
1377 ENDDO
1378 DO I = 1, IM
1379 RHBAR(I) = RHBAR(I) / float(KTCON(I) - KB(I) - 1)
1380 ENDDO
1381 !
1382 ! this section is ready for cloud water
1383 !
1384 if(ncloud.gt.0) THEN
1385 !
1386 ! compute liquid and vapor separation at cloud top
1387 !
1388 DO I = 1, IM
1389 k = KTCON(I)
1390 IF(CNVFLG(I)) THEN
1391 GAMMA = EL2ORC * QESO(I,K) / (TO(I,K)**2)
1392 QRCH = QESO(I,K) &
1393 & + GAMMA * DBYO(I,K) / (HVAP * (1. + GAMMA))
1394 DQ = QCKO(I,K-1) - QRCH
1395 !
1396 ! CHECK IF THERE IS EXCESS MOISTURE TO RELEASE LATENT HEAT
1397 !
1398 IF(DQ.GT.0.) THEN
1399 QLKO_KTCON(I) = dq
1400 QCKO(I,K-1) = QRCH
1401 ENDIF
1402 ENDIF
1403 ENDDO
1404 ENDIF
1405 !
1406 ! CALCULATE CLOUD WORK FUNCTION AT T+DT
1407 !
1408 DO K = 1, KM
1409 DO I = 1, IM
1410 if (k .le. kmax(i)) then
1411 IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
1412 DZ1 = ZO(I,k) - ZO(I,k-1)
1413 GAMMA = EL2ORC * QESO(I,k-1) / (TO(I,k-1)**2)
1414 RFACT = 1. + DELTA * CP * GAMMA &
1415 & * TO(I,k-1) / HVAP
1416 AA1(I) = AA1(I) + &
1417 & DZ1 * (G / (CP * TO(I,k-1))) &
1418 & * DBYO(I,k-1) / (1. + GAMMA) &
1419 & * RFACT
1420 val = 0.
1421 AA1(I)=AA1(I)+ &
1422 & DZ1 * G * DELTA * &
1423 !cmr & MAX( 0.,(QESO(I,k-1) - QO(I,k-1))) &
1424 & MAX(val,(QESO(I,k-1) - QO(I,k-1)))
1425 ENDIF
1426 endif
1427 ENDDO
1428 ENDDO
1429 DO I = 1, IM
1430 IF(CNVFLG(I).AND.AA1(I).LE.0.) DWNFLG(I) = .FALSE.
1431 IF(CNVFLG(I).AND.AA1(I).LE.0.) DWNFLG2(I) = .FALSE.
1432 IF(CNVFLG(I).AND.AA1(I).LE.0.) CNVFLG(I) = .FALSE.
1433 ENDDO
1434 !!
1435 TOTFLG = .TRUE.
1436 DO I = 1, IM
1437 TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
1438 ENDDO
1439 IF(TOTFLG) RETURN
1440 !!
1441 !cccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
1442 !cccc PRINT *, ' AA1(I) BEFORE DWNDRFT =', AA1(I)
1443 !cccc ENDIF
1444 !
1445 !------- DOWNDRAFT CALCULATIONS
1446 !
1447 !
1448 !--- DETERMINE DOWNDRAFT STRENGTH IN TERMS OF WINDSHEAR
1449 !
1450 DO I = 1, IM
1451 IF(CNVFLG(I)) THEN
1452 VSHEAR(I) = 0.
1453 ENDIF
1454 ENDDO
1455 DO K = 1, KM
1456 DO I = 1, IM
1457 if (k .le. kmax(i)) then
1458 IF(K.GE.KB(I).AND.K.LE.KTCON(I).AND.CNVFLG(I)) THEN
1459 shear=rcs(I) * sqrt((UO(I,k+1)-UO(I,k)) ** 2 &
1460 & + (VO(I,k+1)-VO(I,k)) ** 2)
1461 VSHEAR(I) = VSHEAR(I) + SHEAR
1462 ENDIF
1463 endif
1464 ENDDO
1465 ENDDO
1466 DO I = 1, IM
1467 EDT(I) = 0.
1468 IF(CNVFLG(I)) THEN
1469 KNUMB = KTCON(I) - KB(I) + 1
1470 KNUMB = MAX(KNUMB,1)
1471 VSHEAR(I) = 1.E3 * VSHEAR(I) / (ZO(I,KTCON(I))-ZO(I,KB(I)))
1472 E1=1.591-.639*VSHEAR(I) &
1473 & +.0953*(VSHEAR(I)**2)-.00496*(VSHEAR(I)**3)
1474 EDT(I)=1.-E1
1475 !cmr EDT(I) = MIN(EDT(I),.9)
1476 val = .9
1477 EDT(I) = MIN(EDT(I),val)
1478 !cmr EDT(I) = MAX(EDT(I),.0)
1479 val = .0
1480 EDT(I) = MAX(EDT(I),val)
1481 EDTO(I)=EDT(I)
1482 EDTX(I)=EDT(I)
1483 ENDIF
1484 ENDDO
1485 ! DETERMINE DETRAINMENT RATE BETWEEN 1 AND KBDTR
1486 DO I = 1, IM
1487 KBDTR(I) = KBCON(I)
1488 beta = betas
1489 if(SLIMSK(I).eq.1.) beta = betal
1490 IF(CNVFLG(I)) THEN
1491 KBDTR(I) = KBCON(I)
1492 KBDTR(I) = MAX(KBDTR(I),1)
1493 XLAMD(I) = 0.
1494 IF(KBDTR(I).GT.1) THEN
1495 DZ = .5 * ZO(I,KBDTR(I)) + .5 * ZO(I,KBDTR(I)-1) &
1496 & - ZO(I,1)
1497 XLAMD(I) = LOG(BETA) / DZ
1498 ENDIF
1499 ENDIF
1500 ENDDO
1501 ! DETERMINE DOWNDRAFT MASS FLUX
1502 DO K = 1, KM
1503 DO I = 1, IM
1504 IF(k .le. kmax(i)) then
1505 IF(CNVFLG(I)) THEN
1506 ETAD(I,k) = 1.
1507 ENDIF
1508 QRCDO(I,k) = 0.
1509 endif
1510 ENDDO
1511 ENDDO
1512 DO K = KM1, 2, -1
1513 DO I = 1, IM
1514 if (k .le. kbmax(i)) then
1515 IF(CNVFLG(I).AND.K.LT.KBDTR(I)) THEN
1516 DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1517 ETAD(I,k) = ETAD(I,k+1) * EXP(XLAMD(I) * DZ)
1518 ENDIF
1519 endif
1520 ENDDO
1521 ENDDO
1522 K = 1
1523 DO I = 1, IM
1524 IF(CNVFLG(I).AND.KBDTR(I).GT.1) THEN
1525 DZ = .5 * (ZO(I,2) - ZO(I,1))
1526 ETAD(I,k) = ETAD(I,k+1) * EXP(XLAMD(I) * DZ)
1527 ENDIF
1528 ENDDO
1529 !
1530 !--- DOWNDRAFT MOISTURE PROPERTIES
1531 !
1532 DO I = 1, IM
1533 PWEVO(I) = 0.
1534 FLG(I) = CNVFLG(I)
1535 ENDDO
1536 DO I = 1, IM
1537 IF(CNVFLG(I)) THEN
1538 JMN = JMIN(I)
1539 HCDO(I) = HEO(I,JMN)
1540 QCDO(I) = QO(I,JMN)
1541 QRCDO(I,JMN) = QESO(I,JMN)
1542 UCDO(I) = UO(I,JMN)
1543 VCDO(I) = VO(I,JMN)
1544 ENDIF
1545 ENDDO
1546 DO K = KM1, 1, -1
1547 DO I = 1, IM
1548 if (k .le. kmax(i)-1) then
1549 IF(CNVFLG(I).AND.K.LT.JMIN(I)) THEN
1550 DQ = QESO(I,k)
1551 DT = TO(I,k)
1552 GAMMA = EL2ORC * DQ / DT**2
1553 DH = HCDO(I) - HESO(I,k)
1554 QRCDO(I,k) = DQ+(1./HVAP)*(GAMMA/(1.+GAMMA))*DH
1555 DETAD = ETAD(I,k+1) - ETAD(I,k)
1556 PWDO(I,k) = ETAD(I,k+1) * QCDO(I) - &
1557 & ETAD(I,k) * QRCDO(I,k)
1558 PWDO(I,k) = PWDO(I,k) - DETAD * &
1559 & .5 * (QRCDO(I,k) + QRCDO(I,k+1))
1560 QCDO(I) = QRCDO(I,k)
1561 PWEVO(I) = PWEVO(I) + PWDO(I,k)
1562 ENDIF
1563 endif
1564 ENDDO
1565 ENDDO
1566 ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.DWNFLG(I)) THEN
1567 ! PRINT *, ' PWAVO(I), PWEVO(I) =', PWAVO(I), PWEVO(I)
1568 ! ENDIF
1569 !
1570 !--- FINAL DOWNDRAFT STRENGTH DEPENDENT ON PRECIP
1571 !--- EFFICIENCY (EDT), NORMALIZED CONDENSATE (PWAV), AND
1572 !--- EVAPORATE (PWEV)
1573 !
1574 DO I = 1, IM
1575 edtmax = edtmaxl
1576 if(SLIMSK(I).eq.0.) edtmax = edtmaxs
1577 IF(DWNFLG2(I)) THEN
1578 IF(PWEVO(I).LT.0.) THEN
1579 EDTO(I) = -EDTO(I) * PWAVO(I) / PWEVO(I)
1580 EDTO(I) = MIN(EDTO(I),EDTMAX)
1581 ELSE
1582 EDTO(I) = 0.
1583 ENDIF
1584 ELSE
1585 EDTO(I) = 0.
1586 ENDIF
1587 ENDDO
1588 !
1589 !
1590 !--- DOWNDRAFT CLOUDWORK FUNCTIONS
1591 !
1592 !
1593 DO K = KM1, 1, -1
1594 DO I = 1, IM
1595 if (k .le. kmax(i)-1) then
1596 IF(DWNFLG2(I).AND.K.LT.JMIN(I)) THEN
1597 GAMMA = EL2ORC * QESO(I,k+1) / TO(I,k+1)**2
1598 DHH=HCDO(I)
1599 DT=TO(I,k+1)
1600 DG=GAMMA
1601 DH=HESO(I,k+1)
1602 DZ=-1.*(ZO(I,k+1)-ZO(I,k))
1603 AA1(I)=AA1(I)+EDTO(I)*DZ*(G/(CP*DT))*((DHH-DH)/(1.+DG)) &
1604 & *(1.+DELTA*CP*DG*DT/HVAP)
1605 val=0.
1606 AA1(I)=AA1(I)+EDTO(I)* &
1607 !cmr & DZ*G*DELTA*MAX( 0.,(QESO(I,k+1)-QO(I,k+1))) &
1608 & DZ*G*DELTA*MAX(val,(QESO(I,k+1)-QO(I,k+1)))
1609 ENDIF
1610 endif
1611 ENDDO
1612 ENDDO
1613 !cccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.DWNFLG2(I)) THEN
1614 !cccc PRINT *, ' AA1(I) AFTER DWNDRFT =', AA1(I)
1615 !cccc ENDIF
1616 DO I = 1, IM
1617 IF(AA1(I).LE.0.) CNVFLG(I) = .FALSE.
1618 IF(AA1(I).LE.0.) DWNFLG(I) = .FALSE.
1619 IF(AA1(I).LE.0.) DWNFLG2(I) = .FALSE.
1620 ENDDO
1621 !!
1622 TOTFLG = .TRUE.
1623 DO I = 1, IM
1624 TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
1625 ENDDO
1626 IF(TOTFLG) RETURN
1627 !!
1628 !
1629 !
1630 !--- WHAT WOULD THE CHANGE BE, THAT A CLOUD WITH UNIT MASS
1631 !--- WILL DO TO THE ENVIRONMENT?
1632 !
1633 DO K = 1, KM
1634 DO I = 1, IM
1635 IF(k .le. kmax(i) .and. CNVFLG(I)) THEN
1636 DELLAH(I,k) = 0.
1637 DELLAQ(I,k) = 0.
1638 DELLAU(I,k) = 0.
1639 DELLAV(I,k) = 0.
1640 ENDIF
1641 ENDDO
1642 ENDDO
1643 DO I = 1, IM
1644 IF(CNVFLG(I)) THEN
1645 DP = 1000. * DEL(I,1)
1646 DELLAH(I,1) = EDTO(I) * ETAD(I,1) * (HCDO(I) &
1647 & - HEO(I,1)) * G / DP
1648 DELLAQ(I,1) = EDTO(I) * ETAD(I,1) * (QCDO(I) &
1649 & - QO(I,1)) * G / DP
1650 DELLAU(I,1) = EDTO(I) * ETAD(I,1) * (UCDO(I) &
1651 & - UO(I,1)) * G / DP
1652 DELLAV(I,1) = EDTO(I) * ETAD(I,1) * (VCDO(I) &
1653 & - VO(I,1)) * G / DP
1654 ENDIF
1655 ENDDO
1656 !
1657 !--- CHANGED DUE TO SUBSIDENCE AND ENTRAINMENT
1658 !
1659 DO K = 2, KM1
1660 DO I = 1, IM
1661 if (k .le. kmax(i)-1) then
1662 IF(CNVFLG(I).AND.K.LT.KTCON(I)) THEN
1663 AUP = 1.
1664 IF(K.LE.KB(I)) AUP = 0.
1665 ADW = 1.
1666 IF(K.GT.JMIN(I)) ADW = 0.
1667 DV1= HEO(I,k)
1668 DV2 = .5 * (HEO(I,k) + HEO(I,k+1))
1669 DV3= HEO(I,k-1)
1670 DV1Q= QO(I,k)
1671 DV2Q = .5 * (QO(I,k) + QO(I,k+1))
1672 DV3Q= QO(I,k-1)
1673 DV1U= UO(I,k)
1674 DV2U = .5 * (UO(I,k) + UO(I,k+1))
1675 DV3U= UO(I,k-1)
1676 DV1V= VO(I,k)
1677 DV2V = .5 * (VO(I,k) + VO(I,k+1))
1678 DV3V= VO(I,k-1)
1679 DP = 1000. * DEL(I,K)
1680 DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1681 DETA = ETA(I,k) - ETA(I,k-1)
1682 DETAD = ETAD(I,k) - ETAD(I,k-1)
1683 DELLAH(I,k) = DELLAH(I,k) + &
1684 & ((AUP * ETA(I,k) - ADW * EDTO(I) * ETAD(I,k)) * DV1 &
1685 & - (AUP * ETA(I,k-1) - ADW * EDTO(I) * ETAD(I,k-1))* DV3 &
1686 & - AUP * DETA * DV2 &
1687 & + ADW * EDTO(I) * DETAD * HCDO(I)) * G / DP
1688 DELLAQ(I,k) = DELLAQ(I,k) + &
1689 & ((AUP * ETA(I,k) - ADW * EDTO(I) * ETAD(I,k)) * DV1Q &
1690 & - (AUP * ETA(I,k-1) - ADW * EDTO(I) * ETAD(I,k-1))* DV3Q &
1691 & - AUP * DETA * DV2Q &
1692 & +ADW*EDTO(I)*DETAD*.5*(QRCDO(I,k)+QRCDO(I,k-1))) * G / DP
1693 DELLAU(I,k) = DELLAU(I,k) + &
1694 & ((AUP * ETA(I,k) - ADW * EDTO(I) * ETAD(I,k)) * DV1U &
1695 & - (AUP * ETA(I,k-1) - ADW * EDTO(I) * ETAD(I,k-1))* DV3U &
1696 & - AUP * DETA * DV2U &
1697 & + ADW * EDTO(I) * DETAD * UCDO(I) &
1698 & ) * G / DP
1699 DELLAV(I,k) = DELLAV(I,k) + &
1700 & ((AUP * ETA(I,k) - ADW * EDTO(I) * ETAD(I,k)) * DV1V &
1701 & - (AUP * ETA(I,k-1) - ADW * EDTO(I) * ETAD(I,k-1))* DV3V &
1702 & - AUP * DETA * DV2V &
1703 & + ADW * EDTO(I) * DETAD * VCDO(I) &
1704 & ) * G / DP
1705 ENDIF
1706 endif
1707 ENDDO
1708 ENDDO
1709 !
1710 !------- CLOUD TOP
1711 !
1712 DO I = 1, IM
1713 IF(CNVFLG(I)) THEN
1714 INDX = KTCON(I)
1715 DP = 1000. * DEL(I,INDX)
1716 DV1 = HEO(I,INDX-1)
1717 DELLAH(I,INDX) = ETA(I,INDX-1) * &
1718 & (HCKO(I,INDX-1) - DV1) * G / DP
1719 DVQ1 = QO(I,INDX-1)
1720 DELLAQ(I,INDX) = ETA(I,INDX-1) * &
1721 & (QCKO(I,INDX-1) - DVQ1) * G / DP
1722 DV1U = UO(I,INDX-1)
1723 DELLAU(I,INDX) = ETA(I,INDX-1) * &
1724 & (UCKO(I,INDX-1) - DV1U) * G / DP
1725 DV1V = VO(I,INDX-1)
1726 DELLAV(I,INDX) = ETA(I,INDX-1) * &
1727 & (VCKO(I,INDX-1) - DV1V) * G / DP
1728 !
1729 ! cloud water
1730 !
1731 DELLAL(I) = ETA(I,INDX-1) * QLKO_KTCON(I) * g / dp
1732 ENDIF
1733 ENDDO
1734 !
1735 !------- FINAL CHANGED VARIABLE PER UNIT MASS FLUX
1736 !
1737 DO K = 1, KM
1738 DO I = 1, IM
1739 if (k .le. kmax(i)) then
1740 IF(CNVFLG(I).and.k.gt.KTCON(I)) THEN
1741 QO(I,k) = Q1(I,k)
1742 TO(I,k) = T1(I,k)
1743 UO(I,k) = U1(I,k)
1744 VO(I,k) = V1(I,k)
1745 ENDIF
1746 IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
1747 QO(I,k) = DELLAQ(I,k) * MBDT + Q1(I,k)
1748 DELLAT = (DELLAH(I,k) - HVAP * DELLAQ(I,k)) / CP
1749 TO(I,k) = DELLAT * MBDT + T1(I,k)
1750 !cmr QO(I,k) = max(QO(I,k),1.e-10)
1751 val = 1.e-10
1752 QO(I,k) = max(QO(I,k), val )
1753 ENDIF
1754 endif
1755 ENDDO
1756 ENDDO
1757 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1758 !
1759 !--- THE ABOVE CHANGED ENVIRONMENT IS NOW USED TO CALULATE THE
1760 !--- EFFECT THE ARBITRARY CLOUD (WITH UNIT MASS FLUX)
1761 !--- WOULD HAVE ON THE STABILITY,
1762 !--- WHICH THEN IS USED TO CALCULATE THE REAL MASS FLUX,
1763 !--- NECESSARY TO KEEP THIS CHANGE IN BALANCE WITH THE LARGE-SCALE
1764 !--- DESTABILIZATION.
1765 !
1766 !--- ENVIRONMENTAL CONDITIONS AGAIN, FIRST HEIGHTS
1767 !
1768 DO K = 1, KM
1769 DO I = 1, IM
1770 IF(k .le. kmax(i) .and. CNVFLG(I)) THEN
1771 !jfe QESO(I,k) = 10. * FPVS(TO(I,k))
1772 !
1773 QESO(I,k) = 0.01 * fpvs(TO(I,K)) ! fpvs is in Pa
1774 !
1775 QESO(I,k) = EPS * QESO(I,k) / (PFLD(I,k)+EPSM1*QESO(I,k))
1776 !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
1777 val = 1.E-8
1778 QESO(I,k) = MAX(QESO(I,k), val )
1779 TVO(I,k) = TO(I,k) + DELTA * TO(I,k) * QO(I,k)
1780 ENDIF
1781 ENDDO
1782 ENDDO
1783 DO I = 1, IM
1784 IF(CNVFLG(I)) THEN
1785 XAA0(I) = 0.
1786 XPWAV(I) = 0.
1787 ENDIF
1788 ENDDO
1789 !
1790 ! HYDROSTATIC HEIGHT ASSUME ZERO TERR
1791 !
1792 ! DO I = 1, IM
1793 ! IF(CNVFLG(I)) THEN
1794 ! DLNSIG = LOG(PRSL(I,1)/PS(I))
1795 ! ZO(I,1) = TERR - DLNSIG * RD / G * TVO(I,1)
1796 ! ENDIF
1797 ! ENDDO
1798 ! DO K = 2, KM
1799 ! DO I = 1, IM
1800 ! IF(k .le. kmax(i) .and. CNVFLG(I)) THEN
1801 ! DLNSIG = LOG(PRSL(I,K) / PRSL(I,K-1))
1802 ! ZO(I,k) = ZO(I,k-1) - DLNSIG * RD / G
1803 ! & * .5 * (TVO(I,k) + TVO(I,k-1))
1804 ! ENDIF
1805 ! ENDDO
1806 ! ENDDO
1807 !
1808 !--- MOIST STATIC ENERGY
1809 !
1810 DO K = 1, KM1
1811 DO I = 1, IM
1812 IF(k .le. kmax(i)-1 .and. CNVFLG(I)) THEN
1813 DZ = .5 * (ZO(I,k+1) - ZO(I,k))
1814 DP = .5 * (PFLD(I,k+1) - PFLD(I,k))
1815 !jfe ES = 10. * FPVS(TO(I,k+1))
1816 !
1817 ES = 0.01 * fpvs(TO(I,K+1)) ! fpvs is in Pa
1818 !
1819 PPRIME = PFLD(I,k+1) + EPSM1 * ES
1820 QS = EPS * ES / PPRIME
1821 DQSDP = - QS / PPRIME
1822 DESDT = ES * (FACT1 / TO(I,k+1) + FACT2 / (TO(I,k+1)**2))
1823 DQSDT = QS * PFLD(I,k+1) * DESDT / (ES * PPRIME)
1824 GAMMA = EL2ORC * QESO(I,k+1) / (TO(I,k+1)**2)
1825 DT = (G * DZ + HVAP * DQSDP * DP) / (CP * (1. + GAMMA))
1826 DQ = DQSDT * DT + DQSDP * DP
1827 TO(I,k) = TO(I,k+1) + DT
1828 QO(I,k) = QO(I,k+1) + DQ
1829 PO(I,k) = .5 * (PFLD(I,k) + PFLD(I,k+1))
1830 ENDIF
1831 ENDDO
1832 ENDDO
1833 DO K = 1, KM1
1834 DO I = 1, IM
1835 IF(k .le. kmax(i)-1 .and. CNVFLG(I)) THEN
1836 !jfe QESO(I,k) = 10. * FPVS(TO(I,k))
1837 !
1838 QESO(I,k) = 0.01 * fpvs(TO(I,K)) ! fpvs is in Pa
1839 !
1840 QESO(I,k) = EPS * QESO(I,k) / (PO(I,k) + EPSM1 * QESO(I,k))
1841 !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
1842 val1 = 1.E-8
1843 QESO(I,k) = MAX(QESO(I,k), val1)
1844 !cmr QO(I,k) = max(QO(I,k),1.e-10)
1845 val2 = 1.e-10
1846 QO(I,k) = max(QO(I,k), val2 )
1847 ! QO(I,k) = MIN(QO(I,k),QESO(I,k))
1848 HEO(I,k) = .5 * G * (ZO(I,k) + ZO(I,k+1)) + &
1849 & CP * TO(I,k) + HVAP * QO(I,k)
1850 HESO(I,k) = .5 * G * (ZO(I,k) + ZO(I,k+1)) + &
1851 & CP * TO(I,k) + HVAP * QESO(I,k)
1852 ENDIF
1853 ENDDO
1854 ENDDO
1855 DO I = 1, IM
1856 k = kmax(i)
1857 IF(CNVFLG(I)) THEN
1858 HEO(I,k) = G * ZO(I,k) + CP * TO(I,k) + HVAP * QO(I,k)
1859 HESO(I,k) = G * ZO(I,k) + CP * TO(I,k) + HVAP * QESO(I,k)
1860 ! HEO(I,k) = MIN(HEO(I,k),HESO(I,k))
1861 ENDIF
1862 ENDDO
1863 DO I = 1, IM
1864 IF(CNVFLG(I)) THEN
1865 INDX = KB(I)
1866 XHKB(I) = HEO(I,INDX)
1867 XQKB(I) = QO(I,INDX)
1868 HCKO(I,INDX) = XHKB(I)
1869 QCKO(I,INDX) = XQKB(I)
1870 ENDIF
1871 ENDDO
1872 !
1873 !
1874 !**************************** STATIC CONTROL
1875 !
1876 !
1877 !------- MOISTURE AND CLOUD WORK FUNCTIONS
1878 !
1879 DO K = 2, KM1
1880 DO I = 1, IM
1881 if (k .le. kmax(i)-1) then
1882 ! IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LE.KBCON(I)) THEN
1883 IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LE.KTCON(I)) THEN
1884 FACTOR = ETA(I,k-1) / ETA(I,k)
1885 ONEMF = 1. - FACTOR
1886 HCKO(I,k) = FACTOR * HCKO(I,k-1) + ONEMF * &
1887 & .5 * (HEO(I,k) + HEO(I,k+1))
1888 ENDIF
1889 ! IF(CNVFLG(I).AND.K.GT.KBCON(I)) THEN
1890 ! HEO(I,k) = HEO(I,k-1)
1891 ! ENDIF
1892 endif
1893 ENDDO
1894 ENDDO
1895 DO K = 2, KM1
1896 DO I = 1, IM
1897 if (k .le. kmax(i)-1) then
1898 IF(CNVFLG(I).AND.K.GT.KB(I).AND.K.LT.KTCON(I)) THEN
1899 DZ = .5 * (ZO(I,k+1) - ZO(I,k-1))
1900 GAMMA = EL2ORC * QESO(I,k) / (TO(I,k)**2)
1901 XDBY = HCKO(I,k) - HESO(I,k)
1902 !cmr XDBY = MAX(XDBY,0.)
1903 val = 0.
1904 XDBY = MAX(XDBY,val)
1905 XQRCH = QESO(I,k) &
1906 & + GAMMA * XDBY / (HVAP * (1. + GAMMA))
1907 FACTOR = ETA(I,k-1) / ETA(I,k)
1908 ONEMF = 1. - FACTOR
1909 QCKO(I,k) = FACTOR * QCKO(I,k-1) + ONEMF * &
1910 & .5 * (QO(I,k) + QO(I,k+1))
1911 DQ = ETA(I,k) * QCKO(I,k) - ETA(I,k) * XQRCH
1912 IF(DQ.GT.0.) THEN
1913 ETAH = .5 * (ETA(I,k) + ETA(I,k-1))
1914 QLK = DQ / (ETA(I,k) + ETAH * C0 * DZ)
1915 XAA0(I) = XAA0(I) - (ZO(I,k) - ZO(I,k-1)) * G * QLK
1916 XQC = QLK + XQRCH
1917 XPW = ETAH * C0 * DZ * QLK
1918 QCKO(I,k) = XQC
1919 XPWAV(I) = XPWAV(I) + XPW
1920 ENDIF
1921 ENDIF
1922 ! IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LT.KTCON(I)) THEN
1923 IF(CNVFLG(I).AND.K.GT.KBCON(I).AND.K.LE.KTCON(I)) THEN
1924 DZ1 = ZO(I,k) - ZO(I,k-1)
1925 GAMMA = EL2ORC * QESO(I,k-1) / (TO(I,k-1)**2)
1926 RFACT = 1. + DELTA * CP * GAMMA &
1927 & * TO(I,k-1) / HVAP
1928 XDBY = HCKO(I,k-1) - HESO(I,k-1)
1929 XAA0(I) = XAA0(I) &
1930 & + DZ1 * (G / (CP * TO(I,k-1))) &
1931 & * XDBY / (1. + GAMMA) &
1932 & * RFACT
1933 val=0.
1934 XAA0(I)=XAA0(I)+ &
1935 & DZ1 * G * DELTA * &
1936 !cmr & MAX( 0.,(QESO(I,k-1) - QO(I,k-1))) &
1937 & MAX(val,(QESO(I,k-1) - QO(I,k-1)))
1938 ENDIF
1939 endif
1940 ENDDO
1941 ENDDO
1942 !cccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
1943 !cccc PRINT *, ' XAA BEFORE DWNDRFT =', XAA0(I)
1944 !cccc ENDIF
1945 !
1946 !------- DOWNDRAFT CALCULATIONS
1947 !
1948 !
1949 !--- DOWNDRAFT MOISTURE PROPERTIES
1950 !
1951 DO I = 1, IM
1952 XPWEV(I) = 0.
1953 ENDDO
1954 DO I = 1, IM
1955 IF(DWNFLG2(I)) THEN
1956 JMN = JMIN(I)
1957 XHCD(I) = HEO(I,JMN)
1958 XQCD(I) = QO(I,JMN)
1959 QRCD(I,JMN) = QESO(I,JMN)
1960 ENDIF
1961 ENDDO
1962 DO K = KM1, 1, -1
1963 DO I = 1, IM
1964 if (k .le. kmax(i)-1) then
1965 IF(DWNFLG2(I).AND.K.LT.JMIN(I)) THEN
1966 DQ = QESO(I,k)
1967 DT = TO(I,k)
1968 GAMMA = EL2ORC * DQ / DT**2
1969 DH = XHCD(I) - HESO(I,k)
1970 QRCD(I,k)=DQ+(1./HVAP)*(GAMMA/(1.+GAMMA))*DH
1971 DETAD = ETAD(I,k+1) - ETAD(I,k)
1972 XPWD = ETAD(I,k+1) * QRCD(I,k+1) - &
1973 & ETAD(I,k) * QRCD(I,k)
1974 XPWD = XPWD - DETAD * &
1975 & .5 * (QRCD(I,k) + QRCD(I,k+1))
1976 XPWEV(I) = XPWEV(I) + XPWD
1977 ENDIF
1978 endif
1979 ENDDO
1980 ENDDO
1981 !
1982 DO I = 1, IM
1983 edtmax = edtmaxl
1984 if(SLIMSK(I).eq.0.) edtmax = edtmaxs
1985 IF(DWNFLG2(I)) THEN
1986 IF(XPWEV(I).GE.0.) THEN
1987 EDTX(I) = 0.
1988 ELSE
1989 EDTX(I) = -EDTX(I) * XPWAV(I) / XPWEV(I)
1990 EDTX(I) = MIN(EDTX(I),EDTMAX)
1991 ENDIF
1992 ELSE
1993 EDTX(I) = 0.
1994 ENDIF
1995 ENDDO
1996 !
1997 !
1998 !
1999 !--- DOWNDRAFT CLOUDWORK FUNCTIONS
2000 !
2001 !
2002 DO K = KM1, 1, -1
2003 DO I = 1, IM
2004 if (k .le. kmax(i)-1) then
2005 IF(DWNFLG2(I).AND.K.LT.JMIN(I)) THEN
2006 GAMMA = EL2ORC * QESO(I,k+1) / TO(I,k+1)**2
2007 DHH=XHCD(I)
2008 DT= TO(I,k+1)
2009 DG= GAMMA
2010 DH= HESO(I,k+1)
2011 DZ=-1.*(ZO(I,k+1)-ZO(I,k))
2012 XAA0(I)=XAA0(I)+EDTX(I)*DZ*(G/(CP*DT))*((DHH-DH)/(1.+DG)) &
2013 & *(1.+DELTA*CP*DG*DT/HVAP)
2014 val=0.
2015 XAA0(I)=XAA0(I)+EDTX(I)* &
2016 !cmr & DZ*G*DELTA*MAX( 0.,(QESO(I,k+1)-QO(I,k+1))) &
2017 & DZ*G*DELTA*MAX(val,(QESO(I,k+1)-QO(I,k+1)))
2018 ENDIF
2019 endif
2020 ENDDO
2021 ENDDO
2022 !cccc IF(LAT.EQ.LATD.AND.lon.eq.lond.and.DWNFLG2(I)) THEN
2023 !cccc PRINT *, ' XAA AFTER DWNDRFT =', XAA0(I)
2024 !cccc ENDIF
2025 !
2026 ! CALCULATE CRITICAL CLOUD WORK FUNCTION
2027 !
2028 DO I = 1, IM
2029 ACRT(I) = 0.
2030 IF(CNVFLG(I)) THEN
2031 ! IF(CNVFLG(I).AND.SLIMSK(I).NE.1.) THEN
2032 IF(PFLD(I,KTCON(I)).LT.PCRIT(15))THEN
2033 ACRT(I)=ACRIT(15)*(975.-PFLD(I,KTCON(I))) &
2034 & /(975.-PCRIT(15))
2035 ELSE IF(PFLD(I,KTCON(I)).GT.PCRIT(1))THEN
2036 ACRT(I)=ACRIT(1)
2037 ELSE
2038 !cmr K = IFIX((850. - PFLD(I,KTCON(I)))/50.) + 2
2039 K = int((850. - PFLD(I,KTCON(I)))/50.) + 2
2040 K = MIN(K,15)
2041 K = MAX(K,2)
2042 ACRT(I)=ACRIT(K)+(ACRIT(K-1)-ACRIT(K))* &
2043 & (PFLD(I,KTCON(I))-PCRIT(K))/(PCRIT(K-1)-PCRIT(K))
2044 ENDIF
2045 ! ELSE
2046 ! ACRT(I) = .5 * (PFLD(I,KBCON(I)) - PFLD(I,KTCON(I)))
2047 ENDIF
2048 ENDDO
2049 DO I = 1, IM
2050 ACRTFCT(I) = 1.
2051 IF(CNVFLG(I)) THEN
2052 if(SLIMSK(I).eq.1.) THEN
2053 w1 = w1l
2054 w2 = w2l
2055 w3 = w3l
2056 w4 = w4l
2057 else
2058 w1 = w1s
2059 w2 = w2s
2060 w3 = w3s
2061 w4 = w4s
2062 ENDIF
2063 !C IF(CNVFLG(I).AND.SLIMSK(I).EQ.1.) THEN
2064 ! ACRTFCT(I) = PDOT(I) / W3
2065 !
2066 ! modify critical cloud workfunction by cloud base vertical velocity
2067 !
2068 IF(PDOT(I).LE.W4) THEN
2069 ACRTFCT(I) = (PDOT(I) - W4) / (W3 - W4)
2070 ELSEIF(PDOT(I).GE.-W4) THEN
2071 ACRTFCT(I) = - (PDOT(I) + W4) / (W4 - W3)
2072 ELSE
2073 ACRTFCT(I) = 0.
2074 ENDIF
2075 !cmr ACRTFCT(I) = MAX(ACRTFCT(I),-1.)
2076 val1 = -1.
2077 ACRTFCT(I) = MAX(ACRTFCT(I),val1)
2078 !cmr ACRTFCT(I) = MIN(ACRTFCT(I),1.)
2079 val2 = 1.
2080 ACRTFCT(I) = MIN(ACRTFCT(I),val2)
2081 ACRTFCT(I) = 1. - ACRTFCT(I)
2082 !
2083 ! modify ACRTFCT(I) by colume mean rh if RHBAR(I) is greater than 80 percent
2084 !
2085 ! if(RHBAR(I).ge..8) THEN
2086 ! ACRTFCT(I) = ACRTFCT(I) * (.9 - min(RHBAR(I),.9)) * 10.
2087 ! ENDIF
2088 !
2089 ! modify adjustment time scale by cloud base vertical velocity
2090 !
2091 DTCONV(I) = DT2 + max((1800. - DT2),RZERO) * &
2092 & (PDOT(I) - W2) / (W1 - W2)
2093 ! DTCONV(I) = MAX(DTCONV(I), DT2)
2094 ! DTCONV(I) = 1800. * (PDOT(I) - w2) / (w1 - w2)
2095 DTCONV(I) = max(DTCONV(I),dtmin)
2096 DTCONV(I) = min(DTCONV(I),dtmax)
2098 ENDIF
2099 ENDDO
2100 !
2101 !--- LARGE SCALE FORCING
2102 !
2103 DO I= 1, IM
2104 FLG(I) = CNVFLG(I)
2105 IF(CNVFLG(I)) THEN
2106 ! F = AA1(I) / DTCONV(I)
2107 FLD(I) = (AA1(I) - ACRT(I) * ACRTFCT(I)) / DTCONV(I)
2108 IF(FLD(I).LE.0.) FLG(I) = .FALSE.
2109 ENDIF
2110 CNVFLG(I) = FLG(I)
2111 IF(CNVFLG(I)) THEN
2112 ! XAA0(I) = MAX(XAA0(I),0.)
2113 XK(I) = (XAA0(I) - AA1(I)) / MBDT
2114 IF(XK(I).GE.0.) FLG(I) = .FALSE.
2115 ENDIF
2116 !
2117 !--- KERNEL, CLOUD BASE MASS FLUX
2118 !
2119 CNVFLG(I) = FLG(I)
2120 IF(CNVFLG(I)) THEN
2121 XMB(I) = -FLD(I) / XK(I)
2122 XMB(I) = MIN(XMB(I),XMBMAX(I))
2123 ENDIF
2124 ENDDO
2125 ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I)) THEN
2126 ! print *, ' RHBAR(I), ACRTFCT(I) =', RHBAR(I), ACRTFCT(I)
2127 ! PRINT *, ' A1, XA =', AA1(I), XAA0(I)
2128 ! PRINT *, ' XMB(I), ACRT =', XMB(I), ACRT
2129 ! ENDIF
2130 TOTFLG = .TRUE.
2131 DO I = 1, IM
2132 TOTFLG = TOTFLG .AND. (.NOT. CNVFLG(I))
2133 ENDDO
2134 IF(TOTFLG) RETURN
2135 !
2136 ! restore t0 and QO to t1 and q1 in case convection stops
2137 !
2138 do k = 1, km
2139 DO I = 1, IM
2140 if (k .le. kmax(i)) then
2141 TO(I,k) = T1(I,k)
2142 QO(I,k) = Q1(I,k)
2143 !jfe QESO(I,k) = 10. * FPVS(T1(I,k))
2144 !
2145 QESO(I,k) = 0.01 * fpvs(T1(I,K)) ! fpvs is in Pa
2146 !
2147 QESO(I,k) = EPS * QESO(I,k) / (PFLD(I,k) + EPSM1*QESO(I,k))
2148 !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
2149 val = 1.E-8
2150 QESO(I,k) = MAX(QESO(I,k), val )
2151 endif
2152 enddo
2153 enddo
2154 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2155 !
2156 !--- FEEDBACK: SIMPLY THE CHANGES FROM THE CLOUD WITH UNIT MASS FLUX
2157 !--- MULTIPLIED BY THE MASS FLUX NECESSARY TO KEEP THE
2158 !--- EQUILIBRIUM WITH THE LARGER-SCALE.
2159 !
2160 DO I = 1, IM
2161 DELHBAR(I) = 0.
2162 DELQBAR(I) = 0.
2163 DELTBAR(I) = 0.
2164 QCOND(I) = 0.
2165 ENDDO
2166 DO K = 1, KM
2167 DO I = 1, IM
2168 if (k .le. kmax(i)) then
2169 IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
2170 AUP = 1.
2171 IF(K.Le.KB(I)) AUP = 0.
2172 ADW = 1.
2173 IF(K.GT.JMIN(I)) ADW = 0.
2174 DELLAT = (DELLAH(I,k) - HVAP * DELLAQ(I,k)) / CP
2175 T1(I,k) = T1(I,k) + DELLAT * XMB(I) * DT2
2176 Q1(I,k) = Q1(I,k) + DELLAQ(I,k) * XMB(I) * DT2
2177 U1(I,k) = U1(I,k) + DELLAU(I,k) * XMB(I) * DT2
2178 V1(I,k) = V1(I,k) + DELLAV(I,k) * XMB(I) * DT2
2179 DP = 1000. * DEL(I,K)
2180 DELHBAR(I) = DELHBAR(I) + DELLAH(I,k)*XMB(I)*DP/G
2181 DELQBAR(I) = DELQBAR(I) + DELLAQ(I,k)*XMB(I)*DP/G
2182 DELTBAR(I) = DELTBAR(I) + DELLAT*XMB(I)*DP/G
2183 ENDIF
2184 endif
2185 ENDDO
2186 ENDDO
2187 DO K = 1, KM
2188 DO I = 1, IM
2189 if (k .le. kmax(i)) then
2190 IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
2191 !jfe QESO(I,k) = 10. * FPVS(T1(I,k))
2192 !
2193 QESO(I,k) = 0.01 * fpvs(T1(I,K)) ! fpvs is in Pa
2194 !
2195 QESO(I,k) = EPS * QESO(I,k)/(PFLD(I,k) + EPSM1*QESO(I,k))
2196 !cmr QESO(I,k) = MAX(QESO(I,k),1.E-8)
2197 val = 1.E-8
2198 QESO(I,k) = MAX(QESO(I,k), val )
2199 !
2200 ! cloud water
2201 !
2202 if(ncloud.gt.0.and.CNVFLG(I).and.k.eq.KTCON(I)) THEN
2203 tem = DELLAL(I) * XMB(I) * dt2
2204 tem1 = MAX(RZERO, MIN(RONE, (TCR-t1(I,K))*TCRF))
2205 if (QL(I,k,2) .gt. -999.0) then
2206 QL(I,k,1) = QL(I,k,1) + tem * tem1 ! Ice
2207 QL(I,k,2) = QL(I,k,2) + tem *(1.0-tem1) ! Water
2208 else
2209 tem2 = QL(I,k,1) + tem
2210 QL(I,k,1) = tem2 * tem1 ! Ice
2211 QL(I,k,2) = tem2 - QL(I,k,1) ! Water
2212 endif
2213 ! QL(I,k) = QL(I,k) + DELLAL(I) * XMB(I) * dt2
2214 dp = 1000. * del(i,k)
2215 DELLAL(I) = DELLAL(I) * XMB(I) * dp / g
2216 ENDIF
2217 ENDIF
2218 endif
2219 ENDDO
2220 ENDDO
2221 ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I) ) THEN
2222 ! PRINT *, ' DELHBAR, DELQBAR, DELTBAR ='
2223 ! PRINT *, DELHBAR, HVAP*DELQBAR, CP*DELTBAR
2224 ! PRINT *, ' DELLBAR ='
2225 ! PRINT 6003, HVAP*DELLbar
2226 ! PRINT *, ' DELLAQ ='
2227 ! PRINT 6003, (HVAP*DELLAQ(I,k)*XMB(I),K=1,KMAX)
2228 ! PRINT *, ' DELLAT ='
2229 ! PRINT 6003, (DELLAH(i,k)*XMB(I)-HVAP*DELLAQ(I,k)*XMB(I), &
2230 ! & K=1,KMAX)
2231 ! ENDIF
2232 DO I = 1, IM
2233 RNTOT(I) = 0.
2234 DELQEV(I) = 0.
2235 DELQ2(I) = 0.
2236 FLG(I) = CNVFLG(I)
2237 ENDDO
2238 DO K = KM, 1, -1
2239 DO I = 1, IM
2240 if (k .le. kmax(i)) then
2241 IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
2242 AUP = 1.
2243 IF(K.Le.KB(I)) AUP = 0.
2244 ADW = 1.
2245 IF(K.GT.JMIN(I)) ADW = 0.
2246 rain = AUP * PWO(I,k) + ADW * EDTO(I) * PWDO(I,k)
2247 RNTOT(I) = RNTOT(I) + rain * XMB(I) * .001 * dt2
2248 ENDIF
2249 endif
2250 ENDDO
2251 ENDDO
2252 DO K = KM, 1, -1
2253 DO I = 1, IM
2254 if (k .le. kmax(i)) then
2255 DELTV(I) = 0.
2256 DELQ(I) = 0.
2257 QEVAP(I) = 0.
2258 IF(CNVFLG(I).AND.K.LE.KTCON(I)) THEN
2259 AUP = 1.
2260 IF(K.Le.KB(I)) AUP = 0.
2261 ADW = 1.
2262 IF(K.GT.JMIN(I)) ADW = 0.
2263 rain = AUP * PWO(I,k) + ADW * EDTO(I) * PWDO(I,k)
2264 RN(I) = RN(I) + rain * XMB(I) * .001 * dt2
2265 ENDIF
2266 IF(FLG(I).AND.K.LE.KTCON(I)) THEN
2267 evef = EDT(I) * evfact
2268 if(SLIMSK(I).eq.1.) evef=EDT(I) * evfactl
2269 ! if(SLIMSK(I).eq.1.) evef=.07
2270 ! if(SLIMSK(I).ne.1.) evef = 0.
2271 QCOND(I) = EVEF * (Q1(I,k) - QESO(I,k)) &
2272 & / (1. + EL2ORC * QESO(I,k) / T1(I,k)**2)
2273 DP = 1000. * DEL(I,K)
2274 IF(RN(I).GT.0..AND.QCOND(I).LT.0.) THEN
2275 QEVAP(I) = -QCOND(I) * (1.-EXP(-.32*SQRT(DT2*RN(I))))
2276 QEVAP(I) = MIN(QEVAP(I), RN(I)*1000.*G/DP)
2277 DELQ2(I) = DELQEV(I) + .001 * QEVAP(I) * dp / g
2278 ENDIF
2279 if(RN(I).gt.0..and.QCOND(I).LT.0..and. &
2280 & DELQ2(I).gt.RNTOT(I)) THEN
2281 QEVAP(I) = 1000.* g * (RNTOT(I) - DELQEV(I)) / dp
2282 FLG(I) = .false.
2283 ENDIF
2284 IF(RN(I).GT.0..AND.QEVAP(I).gt.0.) THEN
2285 Q1(I,k) = Q1(I,k) + QEVAP(I)
2286 T1(I,k) = T1(I,k) - ELOCP * QEVAP(I)
2287 RN(I) = RN(I) - .001 * QEVAP(I) * DP / G
2288 DELTV(I) = - ELOCP*QEVAP(I)/DT2
2289 DELQ(I) = + QEVAP(I)/DT2
2290 DELQEV(I) = DELQEV(I) + .001*dp*QEVAP(I)/g
2291 ENDIF
2292 DELLAQ(I,k) = DELLAQ(I,k) + DELQ(I) / XMB(I)
2293 DELQBAR(I) = DELQBAR(I) + DELQ(I)*DP/G
2294 DELTBAR(I) = DELTBAR(I) + DELTV(I)*DP/G
2295 ENDIF
2296 endif
2297 ENDDO
2298 ENDDO
2299 ! IF(LAT.EQ.LATD.AND.lon.eq.lond.and.CNVFLG(I) ) THEN
2300 ! PRINT *, ' DELLAH ='
2301 ! PRINT 6003, (DELLAH(k)*XMB(I),K=1,KMAX)
2302 ! PRINT *, ' DELLAQ ='
2303 ! PRINT 6003, (HVAP*DELLAQ(I,k)*XMB(I),K=1,KMAX)
2304 ! PRINT *, ' DELHBAR, DELQBAR, DELTBAR ='
2305 ! PRINT *, DELHBAR, HVAP*DELQBAR, CP*DELTBAR
2306 ! PRINT *, ' PRECIP =', HVAP*RN(I)*1000./DT2
2307 !CCCC PRINT *, ' DELLBAR ='
2308 !CCCC PRINT *, HVAP*DELLbar
2309 ! ENDIF
2310 !
2311 ! PRECIPITATION RATE CONVERTED TO ACTUAL PRECIP
2312 ! IN UNIT OF M INSTEAD OF KG
2313 !
2314 DO I = 1, IM
2315 IF(CNVFLG(I)) THEN
2316 !
2317 ! IN THE EVENT OF UPPER LEVEL RAIN EVAPORATION AND LOWER LEVEL DOWNDRAF
2318 ! MOISTENING, RN CAN BECOME NEGATIVE, IN THIS CASE, WE BACK OUT OF TH
2319 ! HEATING AND THE MOISTENING
2320 !
2321 if(RN(I).lt.0..and..not.FLG(I)) RN(I) = 0.
2322 IF(RN(I).LE.0.) THEN
2323 RN(I) = 0.
2324 ELSE
2325 KTOP(I) = KTCON(I)
2326 KBOT(I) = KBCON(I)
2327 KUO(I) = 1
2328 CLDWRK(I) = AA1(I)
2329 ENDIF
2330 ENDIF
2331 ENDDO
2332 DO K = 1, KM
2333 DO I = 1, IM
2334 if (k .le. kmax(i)) then
2335 IF(CNVFLG(I).AND.RN(I).LE.0.) THEN
2336 T1(I,k) = TO(I,k)
2337 Q1(I,k) = QO(I,k)
2338 ENDIF
2339 endif
2340 ENDDO
2341 ENDDO
2342 !!
2343 RETURN
2344 END SUBROUTINE SASCNV
2346 ! ------------------------------------------------------------------------
2348 SUBROUTINE OLD_ARW_SHALCV(IM,IX,KM,DT,DEL,PRSI,PRSL,PRSLK,KUO,Q,T,DPSHC)
2349 !
2350 USE MODULE_GFS_MACHINE , ONLY : kind_phys
2351 USE MODULE_GFS_PHYSCONS, grav => con_g, CP => con_CP, HVAP => con_HVAP &
2352 &, RD => con_RD
2354 implicit none
2355 !
2356 ! include 'constant.h'
2357 !
2358 integer IM, IX, KM, KUO(IM)
2359 real(kind=kind_phys) DEL(IX,KM), PRSI(IX,KM+1), PRSL(IX,KM), &
2360 & PRSLK(IX,KM), &
2361 & Q(IX,KM), T(IX,KM), DT, DPSHC
2362 !
2363 ! Locals
2364 !
2365 real(kind=kind_phys) ck, cpdt, dmse, dsdz1, dsdz2, &
2366 & dsig, dtodsl, dtodsu, eldq, g, &
2367 & gocp, rtdls
2368 !
2369 integer k,k1,k2,kliftl,kliftu,kt,N2,I,iku,ik1,ik,ii
2370 integer INDEX2(IM), KLCL(IM), KBOT(IM), KTOP(IM),kk &
2371 &, KTOPM(IM)
2372 !!
2373 ! PHYSICAL PARAMETERS
2374 PARAMETER(G=GRAV, GOCP=G/CP)
2375 ! BOUNDS OF PARCEL ORIGIN
2376 PARAMETER(KLIFTL=2,KLIFTU=2)
2377 LOGICAL LSHC(IM)
2378 real(kind=kind_phys) Q2(IM*KM), T2(IM*KM), &
2379 & PRSL2(IM*KM), PRSLK2(IM*KM), &
2380 & AL(IM*(KM-1)), AD(IM*KM), AU(IM*(KM-1))
2381 !-----------------------------------------------------------------------
2382 ! COMPRESS FIELDS TO POINTS WITH NO DEEP CONVECTION
2383 ! AND MOIST STATIC INSTABILITY.
2384 DO I=1,IM
2385 LSHC(I)=.FALSE.
2386 ENDDO
2387 DO K=1,KM-1
2388 DO I=1,IM
2389 IF(KUO(I).EQ.0) THEN
2390 ELDQ = HVAP*(Q(I,K)-Q(I,K+1))
2391 CPDT = CP*(T(I,K)-T(I,K+1))
2392 RTDLS = (PRSL(I,K)-PRSL(I,K+1)) / &
2393 & PRSI(I,K+1)*RD*0.5*(T(I,K)+T(I,K+1))
2394 DMSE = ELDQ+CPDT-RTDLS
2395 LSHC(I) = LSHC(I).OR.DMSE.GT.0.
2396 ENDIF
2397 ENDDO
2398 ENDDO
2399 N2 = 0
2400 DO I=1,IM
2401 IF(LSHC(I)) THEN
2402 N2 = N2 + 1
2403 INDEX2(N2) = I
2404 ENDIF
2405 ENDDO
2406 IF(N2.EQ.0) RETURN
2407 DO K=1,KM
2408 KK = (K-1)*N2
2409 DO I=1,N2
2410 IK = KK + I
2411 ii = index2(i)
2412 Q2(IK) = Q(II,K)
2413 T2(IK) = T(II,K)
2414 PRSL2(IK) = PRSL(II,K)
2415 PRSLK2(IK) = PRSLK(II,K)
2416 ENDDO
2417 ENDDO
2418 do i=1,N2
2419 ktopm(i) = KM
2420 enddo
2421 do k=2,KM
2422 do i=1,N2
2423 ii = index2(i)
2424 if (prsi(ii,1)-prsi(ii,k) .le. dpshc) ktopm(i) = k
2425 enddo
2426 enddo
2428 !-----------------------------------------------------------------------
2429 ! COMPUTE MOIST ADIABAT AND DETERMINE LIMITS OF SHALLOW CONVECTION.
2430 ! CHECK FOR MOIST STATIC INSTABILITY AGAIN WITHIN CLOUD.
2431 CALL MSTADBT3(N2,KM-1,KLIFTL,KLIFTU,PRSL2,PRSLK2,T2,Q2, &
2432 & KLCL,KBOT,KTOP,AL,AU)
2433 DO I=1,N2
2434 KBOT(I) = min(KLCL(I)-1, ktopm(i)-1)
2435 KTOP(I) = min(KTOP(I)+1, ktopm(i))
2436 LSHC(I) = .FALSE.
2437 ENDDO
2438 DO K=1,KM-1
2439 KK = (K-1)*N2
2440 DO I=1,N2
2441 IF(K.GE.KBOT(I).AND.K.LT.KTOP(I)) THEN
2442 IK = KK + I
2443 IKU = IK + N2
2444 ELDQ = HVAP * (Q2(IK)-Q2(IKU))
2445 CPDT = CP * (T2(IK)-T2(IKU))
2446 RTDLS = (PRSL2(IK)-PRSL2(IKU)) / &
2447 & PRSI(index2(i),K+1)*RD*0.5*(T2(IK)+T2(IKU))
2448 DMSE = ELDQ + CPDT - RTDLS
2449 LSHC(I) = LSHC(I).OR.DMSE.GT.0.
2450 AU(IK) = G/RTDLS
2451 ENDIF
2452 ENDDO
2453 ENDDO
2454 K1=KM+1
2455 K2=0
2456 DO I=1,N2
2457 IF(.NOT.LSHC(I)) THEN
2458 KBOT(I) = KM+1
2459 KTOP(I) = 0
2460 ENDIF
2461 K1 = MIN(K1,KBOT(I))
2462 K2 = MAX(K2,KTOP(I))
2463 ENDDO
2464 KT = K2-K1+1
2465 IF(KT.LT.2) RETURN
2466 !-----------------------------------------------------------------------
2467 ! SET EDDY VISCOSITY COEFFICIENT CKU AT SIGMA INTERFACES.
2468 ! COMPUTE DIAGONALS AND RHS FOR TRIDIAGONAL MATRIX SOLVER.
2469 ! EXPAND FINAL FIELDS.
2470 KK = (K1-1) * N2
2471 DO I=1,N2
2472 IK = KK + I
2473 AD(IK) = 1.
2474 ENDDO
2475 !
2476 ! DTODSU=DT/DEL(K1)
2477 DO K=K1,K2-1
2478 ! DTODSL=DTODSU
2479 ! DTODSU= DT/DEL(K+1)
2480 ! DSIG=SL(K)-SL(K+1)
2481 KK = (K-1) * N2
2482 DO I=1,N2
2483 ii = index2(i)
2484 DTODSL = DT/DEL(II,K)
2485 DTODSU = DT/DEL(II,K+1)
2486 DSIG = PRSL(II,K) - PRSL(II,K+1)
2487 IK = KK + I
2488 IKU = IK + N2
2489 IF(K.EQ.KBOT(I)) THEN
2490 CK=1.5
2491 ELSEIF(K.EQ.KTOP(I)-1) THEN
2492 CK=1.
2493 ELSEIF(K.EQ.KTOP(I)-2) THEN
2494 CK=3.
2495 ELSEIF(K.GT.KBOT(I).AND.K.LT.KTOP(I)-2) THEN
2496 CK=5.
2497 ELSE
2498 CK=0.
2499 ENDIF
2500 DSDZ1 = CK*DSIG*AU(IK)*GOCP
2501 DSDZ2 = CK*DSIG*AU(IK)*AU(IK)
2502 AU(IK) = -DTODSL*DSDZ2
2503 AL(IK) = -DTODSU*DSDZ2
2504 AD(IK) = AD(IK)-AU(IK)
2505 AD(IKU) = 1.-AL(IK)
2506 T2(IK) = T2(IK)+DTODSL*DSDZ1
2507 T2(IKU) = T2(IKU)-DTODSU*DSDZ1
2508 ENDDO
2509 ENDDO
2510 IK1=(K1-1)*N2+1
2511 CALL TRIDI2T3(N2,KT,AL(IK1),AD(IK1),AU(IK1),Q2(IK1),T2(IK1), &
2512 & AU(IK1),Q2(IK1),T2(IK1))
2513 DO K=K1,K2
2514 KK = (K-1)*N2
2515 DO I=1,N2
2516 IK = KK + I
2517 Q(INDEX2(I),K) = Q2(IK)
2518 T(INDEX2(I),K) = T2(IK)
2519 ENDDO
2520 ENDDO
2521 !-----------------------------------------------------------------------
2522 RETURN
2523 END SUBROUTINE OLD_ARW_SHALCV
2525 !-----------------------------------------------------------------------
2526 SUBROUTINE TRIDI2T3(L,N,CL,CM,CU,R1,R2,AU,A1,A2)
2527 !yt INCLUDE DBTRIDI2;
2528 !!
2529 USE MODULE_GFS_MACHINE , ONLY : kind_phys
2530 implicit none
2531 integer k,n,l,i
2532 real(kind=kind_phys) fk
2533 !!
2534 real(kind=kind_phys) &
2535 & CL(L,2:N),CM(L,N),CU(L,N-1),R1(L,N),R2(L,N), &
2536 & AU(L,N-1),A1(L,N),A2(L,N)
2537 !-----------------------------------------------------------------------
2538 DO I=1,L
2539 FK=1./CM(I,1)
2540 AU(I,1)=FK*CU(I,1)
2541 A1(I,1)=FK*R1(I,1)
2542 A2(I,1)=FK*R2(I,1)
2543 ENDDO
2544 DO K=2,N-1
2545 DO I=1,L
2546 FK=1./(CM(I,K)-CL(I,K)*AU(I,K-1))
2547 AU(I,K)=FK*CU(I,K)
2548 A1(I,K)=FK*(R1(I,K)-CL(I,K)*A1(I,K-1))
2549 A2(I,K)=FK*(R2(I,K)-CL(I,K)*A2(I,K-1))
2550 ENDDO
2551 ENDDO
2552 DO I=1,L
2553 FK=1./(CM(I,N)-CL(I,N)*AU(I,N-1))
2554 A1(I,N)=FK*(R1(I,N)-CL(I,N)*A1(I,N-1))
2555 A2(I,N)=FK*(R2(I,N)-CL(I,N)*A2(I,N-1))
2556 ENDDO
2557 DO K=N-1,1,-1
2558 DO I=1,L
2559 A1(I,K)=A1(I,K)-AU(I,K)*A1(I,K+1)
2560 A2(I,K)=A2(I,K)-AU(I,K)*A2(I,K+1)
2561 ENDDO
2562 ENDDO
2563 !-----------------------------------------------------------------------
2564 RETURN
2565 END SUBROUTINE TRIDI2T3
2566 !-----------------------------------------------------------------------
2568 SUBROUTINE MSTADBT3(IM,KM,K1,K2,PRSL,PRSLK,TENV,QENV, &
2569 & KLCL,KBOT,KTOP,TCLD,QCLD)
2570 !yt INCLUDE DBMSTADB;
2571 !!
2572 USE MODULE_GFS_MACHINE, ONLY : kind_phys
2573 USE MODULE_GFS_FUNCPHYS, ONLY : FTDP, FTHE, FTLCL, STMA
2574 USE MODULE_GFS_PHYSCONS, EPS => con_eps, EPSM1 => con_epsm1, FV => con_FVirt
2576 implicit none
2577 !!
2578 ! include 'constant.h'
2579 !!
2580 integer k,k1,k2,km,i,im
2581 real(kind=kind_phys) pv,qma,slklcl,tdpd,thelcl,tlcl
2582 real(kind=kind_phys) tma,tvcld,tvenv
2583 !!
2584 real(kind=kind_phys) PRSL(IM,KM), PRSLK(IM,KM), TENV(IM,KM), &
2585 & QENV(IM,KM), TCLD(IM,KM), QCLD(IM,KM)
2586 INTEGER KLCL(IM), KBOT(IM), KTOP(IM)
2587 ! LOCAL ARRAYS
2588 real(kind=kind_phys) SLKMA(IM), THEMA(IM)
2589 !-----------------------------------------------------------------------
2590 ! DETERMINE WARMEST POTENTIAL WET-BULB TEMPERATURE BETWEEN K1 AND K2.
2591 ! COMPUTE ITS LIFTING CONDENSATION LEVEL.
2592 !
2593 DO I=1,IM
2594 SLKMA(I) = 0.
2595 THEMA(I) = 0.
2596 ENDDO
2597 DO K=K1,K2
2598 DO I=1,IM
2599 PV = 1000.0 * PRSL(I,K)*QENV(I,K)/(EPS-EPSM1*QENV(I,K))
2600 TDPD = TENV(I,K)-FTDP(PV)
2601 IF(TDPD.GT.0.) THEN
2602 TLCL = FTLCL(TENV(I,K),TDPD)
2603 SLKLCL = PRSLK(I,K)*TLCL/TENV(I,K)
2604 ELSE
2605 TLCL = TENV(I,K)
2606 SLKLCL = PRSLK(I,K)
2607 ENDIF
2608 THELCL=FTHE(TLCL,SLKLCL)
2609 IF(THELCL.GT.THEMA(I)) THEN
2610 SLKMA(I) = SLKLCL
2611 THEMA(I) = THELCL
2612 ENDIF
2613 ENDDO
2614 ENDDO
2615 !-----------------------------------------------------------------------
2616 ! SET CLOUD TEMPERATURES AND HUMIDITIES WHEREVER THE PARCEL LIFTED UP
2617 ! THE MOIST ADIABAT IS BUOYANT WITH RESPECT TO THE ENVIRONMENT.
2618 DO I=1,IM
2619 KLCL(I)=KM+1
2620 KBOT(I)=KM+1
2621 KTOP(I)=0
2622 ENDDO
2623 DO K=1,KM
2624 DO I=1,IM
2625 TCLD(I,K)=0.
2626 QCLD(I,K)=0.
2627 ENDDO
2628 ENDDO
2629 DO K=K1,KM
2630 DO I=1,IM
2631 IF(PRSLK(I,K).LE.SLKMA(I)) THEN
2632 KLCL(I)=MIN(KLCL(I),K)
2633 CALL STMA(THEMA(I),PRSLK(I,K),TMA,QMA)
2634 ! TMA=FTMA(THEMA(I),PRSLK(I,K),QMA)
2635 TVCLD=TMA*(1.+FV*QMA)
2636 TVENV=TENV(I,K)*(1.+FV*QENV(I,K))
2637 IF(TVCLD.GT.TVENV) THEN
2638 KBOT(I)=MIN(KBOT(I),K)
2639 KTOP(I)=MAX(KTOP(I),K)
2640 TCLD(I,K)=TMA-TENV(I,K)
2641 QCLD(I,K)=QMA-QENV(I,K)
2642 ENDIF
2643 ENDIF
2644 ENDDO
2645 ENDDO
2646 !-----------------------------------------------------------------------
2647 RETURN
2648 END SUBROUTINE MSTADBT3
2650 subroutine sascnvn(im,ix,km,jcap,delt,del,prsl,ps,phil,ql, &
2651 & q1,t1,u1,v1,rcs,cldwrk,rn,kbot,ktop,kcnv,slimsk, &
2652 & dot,ncloud,pgcon,sas_mass_flux)
2653 ! & dot,ncloud,ud_mf,dd_mf,dt_mf)
2654 ! & dot,ncloud,ud_mf,dd_mf,dt_mf,me)
2655 !
2656 ! use machine , only : kind_phys
2657 ! use funcphys , only : fpvs
2658 ! use physcons, grav => con_g, cp => con_cp, hvap => con_hvap &
2659 USE MODULE_GFS_MACHINE, ONLY : kind_phys
2660 USE MODULE_GFS_FUNCPHYS, ONLY : fpvs
2661 USE MODULE_GFS_PHYSCONS, grav => con_g, cp => con_cp &
2662 &, hvap => con_hvap &
2663 &, rv => con_rv, fv => con_fvirt, t0c => con_t0c &
2664 &, cvap => con_cvap, cliq => con_cliq &
2665 &, eps => con_eps, epsm1 => con_epsm1
2666 implicit none
2667 !
2668 integer im, ix, km, jcap, ncloud, &
2669 & kbot(im), ktop(im), kcnv(im)
2670 ! &, me
2671 real(kind=kind_phys) delt,sas_mass_flux
2672 real(kind=kind_phys) ps(im), del(ix,km), prsl(ix,km), &
2673 & ql(ix,km,2),q1(ix,km), t1(ix,km), &
2674 & u1(ix,km), v1(ix,km), rcs(im), &
2675 & cldwrk(im), rn(im), slimsk(im), &
2676 & dot(ix,km), phil(ix,km)
2677 ! hchuang code change mass flux output
2678 ! &, ud_mf(im,km),dd_mf(im,km),dt_mf(im,km)
2679 !
2680 integer i, j, indx, jmn, k, kk, latd, lond, km1
2681 !
2682 real(kind=kind_phys) clam, cxlamu, xlamde, xlamdd
2683 !
2684 real(kind=kind_phys) adw, aup, aafac, &
2685 & beta, betal, betas, &
2686 & c0, cpoel, dellat, delta, &
2687 & desdt, deta, detad, dg, &
2688 & dh, dhh, dlnsig, dp, &
2689 & dq, dqsdp, dqsdt, dt, &
2690 & dt2, dtmax, dtmin, dv1h, &
2691 & dv1q, dv2h, dv2q, dv1u, &
2692 & dv1v, dv2u, dv2v, dv3q, &
2693 & dv3h, dv3u, dv3v, &
2694 & dz, dz1, e1, edtmax, &
2695 & edtmaxl, edtmaxs, el2orc, elocp, &
2696 & es, etah, cthk, dthk, &
2697 & evef, evfact, evfactl, fact1, &
2698 & fact2, factor, fjcap, fkm, &
2699 & g, gamma, pprime, &
2700 & qlk, qrch, qs, c1, &
2701 & rain, rfact, shear, tem1, &
2702 & tem2, terr, val, val1, &
2703 & val2, w1, w1l, w1s, &
2704 & w2, w2l, w2s, w3, &
2705 & w3l, w3s, w4, w4l, &
2706 & w4s, xdby, xpw, xpwd, &
2707 & xqrch, mbdt, tem, &
2708 & ptem, ptem1
2709 !
2710 real(kind=kind_phys), intent(in) :: pgcon
2712 integer kb(im), kbcon(im), kbcon1(im), &
2713 & ktcon(im), ktcon1(im), &
2714 & jmin(im), lmin(im), kbmax(im), &
2715 & kbm(im), kmax(im)
2716 !
2717 real(kind=kind_phys) aa1(im), acrt(im), acrtfct(im), &
2718 & delhbar(im), delq(im), delq2(im), &
2719 & delqbar(im), delqev(im), deltbar(im), &
2720 & deltv(im), dtconv(im), edt(im), &
2721 & edto(im), edtx(im), fld(im), &
2722 & hcdo(im,km), hmax(im), hmin(im), &
2723 & ucdo(im,km), vcdo(im,km),aa2(im), &
2724 & pbcdif(im), pdot(im), po(im,km), &
2725 & pwavo(im), pwevo(im), xlamud(im), &
2726 & qcdo(im,km), qcond(im), qevap(im), &
2727 & rntot(im), vshear(im), xaa0(im), &
2728 & xk(im), xlamd(im), &
2729 & xmb(im), xmbmax(im), xpwav(im), &
2730 & xpwev(im), delubar(im),delvbar(im)
2731 !cj
2732 real(kind=kind_phys) cincr, cincrmax, cincrmin
2733 real(kind=kind_phys) xmbmx1
2734 !cj
2735 !c physical parameters
2736 parameter(g=grav)
2737 parameter(cpoel=cp/hvap,elocp=hvap/cp, &
2738 & el2orc=hvap*hvap/(rv*cp))
2739 parameter(terr=0.,c0=.002,c1=.002,delta=fv)
2740 parameter(fact1=(cvap-cliq)/rv,fact2=hvap/rv-fact1*t0c)
2741 parameter(cthk=150.,cincrmax=180.,cincrmin=120.,dthk=25.)
2742 !c local variables and arrays
2743 real(kind=kind_phys) pfld(im,km),to(im,km), qo(im,km), &
2744 & uo(im,km), vo(im,km), qeso(im,km)
2745 !c cloud water
2746 real(kind=kind_phys)qlko_ktcon(im),dellal(im,km),tvo(im,km), &
2747 & dbyo(im,km), zo(im,km), xlamue(im,km), &
2748 & fent1(im,km),fent2(im,km), frh(im,km), &
2749 & heo(im,km), heso(im,km), &
2750 & qrcd(im,km), dellah(im,km), dellaq(im,km), &
2751 & dellau(im,km),dellav(im,km), hcko(im,km), &
2752 & ucko(im,km), vcko(im,km), qcko(im,km), &
2753 & eta(im,km), etad(im,km), zi(im,km), &
2754 & qrcdo(im,km),pwo(im,km), pwdo(im,km), &
2755 & tx1(im), sumx(im)
2756 ! &, rhbar(im)
2757 !
2758 logical totflg, cnvflg(im), flg(im)
2759 !
2760 real(kind=kind_phys) pcrit(15), acritt(15), acrit(15)
2761 ! save pcrit, acritt
2762 data pcrit/850.,800.,750.,700.,650.,600.,550.,500.,450.,400.,&
2763 & 350.,300.,250.,200.,150./
2764 data acritt/.0633,.0445,.0553,.0664,.075,.1082,.1521,.2216, &
2765 & .3151,.3677,.41,.5255,.7663,1.1686,1.6851/
2766 !c gdas derived acrit
2767 !c data acritt/.203,.515,.521,.566,.625,.665,.659,.688,
2768 !c & .743,.813,.886,.947,1.138,1.377,1.896/
2769 real(kind=kind_phys) tf, tcr, tcrf
2770 parameter (tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf))
2771 !
2772 !c-----------------------------------------------------------------------
2773 !
2775 km1 = km - 1
2776 !c
2777 !c initialize arrays
2778 !c
2779 do i=1,im
2780 kcnv(i)=0
2781 cnvflg(i) = .true.
2782 rn(i)=0.
2783 kbot(i)=km+1
2784 ktop(i)=0
2785 kbcon(i)=km
2786 ktcon(i)=1
2787 dtconv(i) = 3600.
2788 cldwrk(i) = 0.
2789 pdot(i) = 0.
2790 pbcdif(i)= 0.
2791 lmin(i) = 1
2792 jmin(i) = 1
2793 qlko_ktcon(i) = 0.
2794 edt(i) = 0.
2795 edto(i) = 0.
2796 edtx(i) = 0.
2797 acrt(i) = 0.
2798 acrtfct(i) = 1.
2799 aa1(i) = 0.
2800 aa2(i) = 0.
2801 xaa0(i) = 0.
2802 pwavo(i)= 0.
2803 pwevo(i)= 0.
2804 xpwav(i)= 0.
2805 xpwev(i)= 0.
2806 vshear(i) = 0.
2807 enddo
2808 ! hchuang code change
2809 ! do k = 1, km
2810 ! do i = 1, im
2811 ! ud_mf(i,k) = 0.
2812 ! dd_mf(i,k) = 0.
2813 ! dt_mf(i,k) = 0.
2814 ! enddo
2815 ! enddo
2816 !c
2817 do k = 1, 15
2818 acrit(k) = acritt(k) * (975. - pcrit(k))
2819 enddo
2820 dt2 = delt
2821 val = 1200.
2822 dtmin = max(dt2, val )
2823 val = 3600.
2824 dtmax = max(dt2, val )
2825 !c model tunable parameters are all here
2826 mbdt = 10.
2827 edtmaxl = .3
2828 edtmaxs = .3
2829 clam = .1
2830 aafac = .1
2831 ! betal = .15
2832 ! betas = .15
2833 betal = .05
2834 betas = .05
2835 !c evef = 0.07
2836 evfact = 0.3
2837 evfactl = 0.3
2838 #if ( EM_CORE == 1 )
2839 ! HAWAII TEST - ZCX
2840 BETAl = .05
2841 betas = .05
2842 evfact = 0.5
2843 evfactl = 0.5
2844 #endif
2845 !
2846 cxlamu = 1.0e-4
2847 xlamde = 1.0e-4
2848 xlamdd = 1.0e-4
2849 !
2850 fjcap = (float(jcap) / 126.) ** 2
2851 val = 1.
2852 fjcap = max(fjcap,val)
2853 fkm = (float(km) / 28.) ** 2
2854 fkm = max(fkm,val)
2855 w1l = -8.e-3
2856 w2l = -4.e-2
2857 w3l = -5.e-3
2858 w4l = -5.e-4
2859 w1s = -2.e-4
2860 w2s = -2.e-3
2861 w3s = -1.e-3
2862 w4s = -2.e-5
2863 !c
2864 !c define top layer for search of the downdraft originating layer
2865 !c and the maximum thetae for updraft
2866 !c
2867 do i=1,im
2868 kbmax(i) = km
2869 kbm(i) = km
2870 kmax(i) = km
2871 tx1(i) = 1.0 / ps(i)
2872 enddo
2873 !
2874 do k = 1, km
2875 do i=1,im
2876 IF (prSL(I,K)*tx1(I) .GT. 0.04) KMAX(I) = MIN(KM,K + 1)
2877 !2011bugfix if (prsl(i,k)*tx1(i) .gt. 0.04) kmax(i) = k + 1
2878 if (prsl(i,k)*tx1(i) .gt. 0.45) kbmax(i) = k + 1
2879 if (prsl(i,k)*tx1(i) .gt. 0.70) kbm(i) = k + 1
2880 enddo
2881 enddo
2882 do i=1,im
2883 kbmax(i) = min(kbmax(i),kmax(i))
2884 kbm(i) = min(kbm(i),kmax(i))
2885 enddo
2886 !c
2887 !c hydrostatic height assume zero terr and initially assume
2888 !c updraft entrainment rate as an inverse function of height
2889 !c
2890 do k = 1, km
2891 do i=1,im
2892 zo(i,k) = phil(i,k) / g
2893 enddo
2894 enddo
2895 do k = 1, km1
2896 do i=1,im
2897 zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1))
2898 xlamue(i,k) = clam / zi(i,k)
2899 enddo
2900 enddo
2901 !c
2902 !c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2903 !c convert surface pressure to mb from cb
2904 !c
2905 do k = 1, km
2906 do i = 1, im
2907 if (k .le. kmax(i)) then
2908 pfld(i,k) = prsl(i,k) * 10.0
2909 eta(i,k) = 1.
2910 fent1(i,k)= 1.
2911 fent2(i,k)= 1.
2912 frh(i,k) = 0.
2913 hcko(i,k) = 0.
2914 qcko(i,k) = 0.
2915 ucko(i,k) = 0.
2916 vcko(i,k) = 0.
2917 etad(i,k) = 1.
2918 hcdo(i,k) = 0.
2919 qcdo(i,k) = 0.
2920 ucdo(i,k) = 0.
2921 vcdo(i,k) = 0.
2922 qrcd(i,k) = 0.
2923 qrcdo(i,k)= 0.
2924 dbyo(i,k) = 0.
2925 pwo(i,k) = 0.
2926 pwdo(i,k) = 0.
2927 dellal(i,k) = 0.
2928 to(i,k) = t1(i,k)
2929 qo(i,k) = q1(i,k)
2930 uo(i,k) = u1(i,k) * rcs(i)
2931 vo(i,k) = v1(i,k) * rcs(i)
2932 endif
2933 enddo
2934 enddo
2935 !c
2936 !c column variables
2937 !c p is pressure of the layer (mb)
2938 !c t is temperature at t-dt (k)..tn
2939 !c q is mixing ratio at t-dt (kg/kg)..qn
2940 !c to is temperature at t+dt (k)... this is after advection and turbulan
2941 !c qo is mixing ratio at t+dt (kg/kg)..q1
2942 !c
2943 do k = 1, km
2944 do i=1,im
2945 if (k .le. kmax(i)) then
2946 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
2947 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
2948 val1 = 1.e-8
2949 qeso(i,k) = max(qeso(i,k), val1)
2950 val2 = 1.e-10
2951 qo(i,k) = max(qo(i,k), val2 )
2952 ! qo(i,k) = min(qo(i,k),qeso(i,k))
2953 ! tvo(i,k) = to(i,k) + delta * to(i,k) * qo(i,k)
2954 endif
2955 enddo
2956 enddo
2957 !c
2958 !c compute moist static energy
2959 !c
2960 do k = 1, km
2961 do i=1,im
2962 if (k .le. kmax(i)) then
2963 ! tem = g * zo(i,k) + cp * to(i,k)
2964 tem = phil(i,k) + cp * to(i,k)
2965 heo(i,k) = tem + hvap * qo(i,k)
2966 heso(i,k) = tem + hvap * qeso(i,k)
2967 !c heo(i,k) = min(heo(i,k),heso(i,k))
2968 endif
2969 enddo
2970 enddo
2971 !c
2972 !c determine level with largest moist static energy
2973 !c this is the level where updraft starts
2974 !c
2975 do i=1,im
2976 hmax(i) = heo(i,1)
2977 kb(i) = 1
2978 enddo
2979 do k = 2, km
2980 do i=1,im
2981 if (k .le. kbm(i)) then
2982 if(heo(i,k).gt.hmax(i)) then
2983 kb(i) = k
2984 hmax(i) = heo(i,k)
2985 endif
2986 endif
2987 enddo
2988 enddo
2989 !c
2990 do k = 1, km1
2991 do i=1,im
2992 if (k .le. kmax(i)-1) then
2993 dz = .5 * (zo(i,k+1) - zo(i,k))
2994 dp = .5 * (pfld(i,k+1) - pfld(i,k))
2995 es = 0.01 * fpvs(to(i,k+1)) ! fpvs is in pa
2996 pprime = pfld(i,k+1) + epsm1 * es
2997 qs = eps * es / pprime
2998 dqsdp = - qs / pprime
2999 desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
3000 dqsdt = qs * pfld(i,k+1) * desdt / (es * pprime)
3001 gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
3002 dt = (g * dz + hvap * dqsdp * dp) / (cp * (1. + gamma))
3003 dq = dqsdt * dt + dqsdp * dp
3004 to(i,k) = to(i,k+1) + dt
3005 qo(i,k) = qo(i,k+1) + dq
3006 po(i,k) = .5 * (pfld(i,k) + pfld(i,k+1))
3007 endif
3008 enddo
3009 enddo
3010 !
3011 do k = 1, km1
3012 do i=1,im
3013 if (k .le. kmax(i)-1) then
3014 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
3015 qeso(i,k) = eps * qeso(i,k) / (po(i,k) + epsm1*qeso(i,k))
3016 val1 = 1.e-8
3017 qeso(i,k) = max(qeso(i,k), val1)
3018 val2 = 1.e-10
3019 qo(i,k) = max(qo(i,k), val2 )
3020 ! qo(i,k) = min(qo(i,k),qeso(i,k))
3021 val1 = 1.0
3022 frh(i,k) = 1. - min(qo(i,k)/qeso(i,k), val1)
3023 heo(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) + &
3024 & cp * to(i,k) + hvap * qo(i,k)
3025 heso(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) + &
3026 & cp * to(i,k) + hvap * qeso(i,k)
3027 uo(i,k) = .5 * (uo(i,k) + uo(i,k+1))
3028 vo(i,k) = .5 * (vo(i,k) + vo(i,k+1))
3029 endif
3030 enddo
3031 enddo
3032 !c
3033 !c look for the level of free convection as cloud base
3034 !c
3035 do i=1,im
3036 flg(i) = .true.
3037 kbcon(i) = kmax(i)
3038 enddo
3039 do k = 1, km1
3040 do i=1,im
3041 if (flg(i).and.k.le.kbmax(i)) then
3042 if(k.gt.kb(i).and.heo(i,kb(i)).gt.heso(i,k)) then
3043 kbcon(i) = k
3044 flg(i) = .false.
3045 endif
3046 endif
3047 enddo
3048 enddo
3049 !c
3050 do i=1,im
3051 if(kbcon(i).eq.kmax(i)) cnvflg(i) = .false.
3052 enddo
3053 !!
3054 totflg = .true.
3055 do i=1,im
3056 totflg = totflg .and. (.not. cnvflg(i))
3057 enddo
3058 if(totflg) return
3059 !!
3060 !c
3061 !c determine critical convective inhibition
3062 !c as a function of vertical velocity at cloud base.
3063 !c
3064 do i=1,im
3065 if(cnvflg(i)) then
3066 pdot(i) = 10.* dot(i,kbcon(i))
3067 endif
3068 enddo
3069 do i=1,im
3070 if(cnvflg(i)) then
3071 if(slimsk(i).eq.1.) then
3072 w1 = w1l
3073 w2 = w2l
3074 w3 = w3l
3075 w4 = w4l
3076 else
3077 w1 = w1s
3078 w2 = w2s
3079 w3 = w3s
3080 w4 = w4s
3081 endif
3082 if(pdot(i).le.w4) then
3083 tem = (pdot(i) - w4) / (w3 - w4)
3084 elseif(pdot(i).ge.-w4) then
3085 tem = - (pdot(i) + w4) / (w4 - w3)
3086 else
3087 tem = 0.
3088 endif
3089 val1 = -1.
3090 tem = max(tem,val1)
3091 val2 = 1.
3092 tem = min(tem,val2)
3093 tem = 1. - tem
3094 tem1= .5*(cincrmax-cincrmin)
3095 cincr = cincrmax - tem * tem1
3096 pbcdif(i) = pfld(i,kb(i)) - pfld(i,kbcon(i))
3097 if(pbcdif(i).gt.cincr) then
3098 cnvflg(i) = .false.
3099 endif
3100 endif
3101 enddo
3102 !!
3103 totflg = .true.
3104 do i=1,im
3105 totflg = totflg .and. (.not. cnvflg(i))
3106 enddo
3107 if(totflg) return
3108 !!
3109 !c
3110 !c assume that updraft entrainment rate above cloud base is
3111 !c same as that at cloud base
3112 !c
3113 do k = 2, km1
3114 do i=1,im
3115 if(cnvflg(i).and. &
3116 & (k.gt.kbcon(i).and.k.lt.kmax(i))) then
3117 xlamue(i,k) = xlamue(i,kbcon(i))
3118 endif
3119 enddo
3120 enddo
3121 !c
3122 !c assume the detrainment rate for the updrafts to be same as
3123 !c the entrainment rate at cloud base
3124 !c
3125 do i = 1, im
3126 if(cnvflg(i)) then
3127 xlamud(i) = xlamue(i,kbcon(i))
3128 endif
3129 enddo
3130 !c
3131 !c functions rapidly decreasing with height, mimicking a cloud ensemble
3132 !c (Bechtold et al., 2008)
3133 !c
3134 do k = 2, km1
3135 do i=1,im
3136 if(cnvflg(i).and. &
3137 & (k.gt.kbcon(i).and.k.lt.kmax(i))) then
3138 tem = qeso(i,k)/qeso(i,kbcon(i))
3139 fent1(i,k) = tem**2
3140 fent2(i,k) = tem**3
3141 endif
3142 enddo
3143 enddo
3144 !c
3145 !c final entrainment rate as the sum of turbulent part and organized entrainment
3146 !c depending on the environmental relative humidity
3147 !c (Bechtold et al., 2008)
3148 !c
3149 do k = 2, km1
3150 do i=1,im
3151 if(cnvflg(i).and. &
3152 & (k.ge.kbcon(i).and.k.lt.kmax(i))) then
3153 tem = cxlamu * frh(i,k) * fent2(i,k)
3154 xlamue(i,k) = xlamue(i,k)*fent1(i,k) + tem
3155 endif
3156 enddo
3157 enddo
3158 !c
3159 !c determine updraft mass flux for the subcloud layers
3160 !c
3161 do k = km1, 1, -1
3162 do i = 1, im
3163 if (cnvflg(i)) then
3164 if(k.lt.kbcon(i).and.k.ge.kb(i)) then
3165 dz = zi(i,k+1) - zi(i,k)
3166 ptem = 0.5*(xlamue(i,k)+xlamue(i,k+1))-xlamud(i)
3167 eta(i,k) = eta(i,k+1) / (1. + ptem * dz)
3168 endif
3169 endif
3170 enddo
3171 enddo
3172 !c
3173 !c compute mass flux above cloud base
3174 !c
3175 do k = 2, km1
3176 do i = 1, im
3177 if(cnvflg(i))then
3178 if(k.gt.kbcon(i).and.k.lt.kmax(i)) then
3179 dz = zi(i,k) - zi(i,k-1)
3180 ptem = 0.5*(xlamue(i,k)+xlamue(i,k-1))-xlamud(i)
3181 eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
3182 endif
3183 endif
3184 enddo
3185 enddo
3186 !c
3187 !c compute updraft cloud properties
3188 !c
3189 do i = 1, im
3190 if(cnvflg(i)) then
3191 indx = kb(i)
3192 hcko(i,indx) = heo(i,indx)
3193 ucko(i,indx) = uo(i,indx)
3194 vcko(i,indx) = vo(i,indx)
3195 pwavo(i) = 0.
3196 endif
3197 enddo
3198 !c
3199 !c cloud property is modified by the entrainment process
3200 !c
3201 do k = 2, km1
3202 do i = 1, im
3203 if (cnvflg(i)) then
3204 if(k.gt.kb(i).and.k.lt.kmax(i)) then
3205 dz = zi(i,k) - zi(i,k-1)
3206 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
3207 tem1 = 0.5 * xlamud(i) * dz
3208 factor = 1. + tem - tem1
3209 ptem = 0.5 * tem + pgcon
3210 ptem1= 0.5 * tem - pgcon
3211 hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5* &
3212 & (heo(i,k)+heo(i,k-1)))/factor
3213 ucko(i,k) = ((1.-tem1)*ucko(i,k-1)+ptem*uo(i,k) &
3214 & +ptem1*uo(i,k-1))/factor
3215 vcko(i,k) = ((1.-tem1)*vcko(i,k-1)+ptem*vo(i,k) &
3216 & +ptem1*vo(i,k-1))/factor
3217 dbyo(i,k) = hcko(i,k) - heso(i,k)
3218 endif
3219 endif
3220 enddo
3221 enddo
3222 !c
3223 !c taking account into convection inhibition due to existence of
3224 !c dry layers below cloud base
3225 !c
3226 do i=1,im
3227 flg(i) = cnvflg(i)
3228 kbcon1(i) = kmax(i)
3229 enddo
3230 do k = 2, km1
3231 do i=1,im
3232 if (flg(i).and.k.lt.kmax(i)) then
3233 if(k.ge.kbcon(i).and.dbyo(i,k).gt.0.) then
3234 kbcon1(i) = k
3235 flg(i) = .false.
3236 endif
3237 endif
3238 enddo
3239 enddo
3240 do i=1,im
3241 if(cnvflg(i)) then
3242 if(kbcon1(i).eq.kmax(i)) cnvflg(i) = .false.
3243 endif
3244 enddo
3245 do i=1,im
3246 if(cnvflg(i)) then
3247 tem = pfld(i,kbcon(i)) - pfld(i,kbcon1(i))
3248 if(tem.gt.dthk) then
3249 cnvflg(i) = .false.
3250 endif
3251 endif
3252 enddo
3253 !!
3254 totflg = .true.
3255 do i = 1, im
3256 totflg = totflg .and. (.not. cnvflg(i))
3257 enddo
3258 if(totflg) return
3259 !!
3260 !c
3261 !c determine first guess cloud top as the level of zero buoyancy
3262 !c
3263 do i = 1, im
3264 flg(i) = cnvflg(i)
3265 ktcon(i) = 1
3266 enddo
3267 do k = 2, km1
3268 do i = 1, im
3269 if (flg(i).and.k .lt. kmax(i)) then
3270 if(k.gt.kbcon1(i).and.dbyo(i,k).lt.0.) then
3271 ktcon(i) = k
3272 flg(i) = .false.
3273 endif
3274 endif
3275 enddo
3276 enddo
3277 !c
3278 do i = 1, im
3279 if(cnvflg(i)) then
3280 tem = pfld(i,kbcon(i))-pfld(i,ktcon(i))
3281 if(tem.lt.cthk) cnvflg(i) = .false.
3282 endif
3283 enddo
3284 !!
3285 totflg = .true.
3286 do i = 1, im
3287 totflg = totflg .and. (.not. cnvflg(i))
3288 enddo
3289 if(totflg) return
3290 !!
3291 !c
3292 !c search for downdraft originating level above theta-e minimum
3293 !c
3294 do i = 1, im
3295 if(cnvflg(i)) then
3296 hmin(i) = heo(i,kbcon1(i))
3297 lmin(i) = kbmax(i)
3298 jmin(i) = kbmax(i)
3299 endif
3300 enddo
3301 do k = 2, km1
3302 do i = 1, im
3303 if (cnvflg(i) .and. k .le. kbmax(i)) then
3304 if(k.gt.kbcon1(i).and.heo(i,k).lt.hmin(i)) then
3305 lmin(i) = k + 1
3306 hmin(i) = heo(i,k)
3307 endif
3308 endif
3309 enddo
3310 enddo
3311 !c
3312 !c make sure that jmin(i) is within the cloud
3313 !c
3314 do i = 1, im
3315 if(cnvflg(i)) then
3316 jmin(i) = min(lmin(i),ktcon(i)-1)
3317 jmin(i) = max(jmin(i),kbcon1(i)+1)
3318 if(jmin(i).ge.ktcon(i)) cnvflg(i) = .false.
3319 endif
3320 enddo
3321 !c
3322 !c specify upper limit of mass flux at cloud base
3323 !c
3324 do i = 1, im
3325 if(cnvflg(i)) then
3326 ! xmbmax(i) = .1
3327 !
3328 k = kbcon(i)
3329 dp = 1000. * del(i,k)
3330 xmbmax(i) = dp / (g * dt2)
3331 xmbmax(i) = min(sas_mass_flux,xmbmax(i))
3332 !
3333 ! tem = dp / (g * dt2)
3334 ! xmbmax(i) = min(tem, xmbmax(i))
3335 endif
3336 enddo
3337 !c
3338 !c compute cloud moisture property and precipitation
3339 !c
3340 do i = 1, im
3341 if (cnvflg(i)) then
3342 aa1(i) = 0.
3343 qcko(i,kb(i)) = qo(i,kb(i))
3344 ! rhbar(i) = 0.
3345 endif
3346 enddo
3347 do k = 2, km1
3348 do i = 1, im
3349 if (cnvflg(i)) then
3350 if(k.gt.kb(i).and.k.lt.ktcon(i)) then
3351 dz = zi(i,k) - zi(i,k-1)
3352 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
3353 qrch = qeso(i,k) &
3354 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
3355 !cj
3356 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
3357 tem1 = 0.5 * xlamud(i) * dz
3358 factor = 1. + tem - tem1
3359 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* &
3360 & (qo(i,k)+qo(i,k-1)))/factor
3361 !cj
3362 dq = eta(i,k) * (qcko(i,k) - qrch)
3363 !c
3364 ! rhbar(i) = rhbar(i) + qo(i,k) / qeso(i,k)
3365 !c
3366 !c check if there is excess moisture to release latent heat
3367 !c
3368 if(k.ge.kbcon(i).and.dq.gt.0.) then
3369 etah = .5 * (eta(i,k) + eta(i,k-1))
3370 if(ncloud.gt.0..and.k.gt.jmin(i)) then
3371 dp = 1000. * del(i,k)
3372 qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
3373 dellal(i,k) = etah * c1 * dz * qlk * g / dp
3374 else
3375 qlk = dq / (eta(i,k) + etah * c0 * dz)
3376 endif
3377 aa1(i) = aa1(i) - dz * g * qlk
3378 qcko(i,k) = qlk + qrch
3379 pwo(i,k) = etah * c0 * dz * qlk
3380 pwavo(i) = pwavo(i) + pwo(i,k)
3381 endif
3382 endif
3383 endif
3384 enddo
3385 enddo
3386 !c
3387 ! do i = 1, im
3388 ! if(cnvflg(i)) then
3389 ! indx = ktcon(i) - kb(i) - 1
3390 ! rhbar(i) = rhbar(i) / float(indx)
3391 ! endif
3392 ! enddo
3393 !c
3394 !c calculate cloud work function
3395 !c
3396 do k = 2, km1
3397 do i = 1, im
3398 if (cnvflg(i)) then
3399 if(k.ge.kbcon(i).and.k.lt.ktcon(i)) then
3400 dz1 = zo(i,k+1) - zo(i,k)
3401 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
3402 rfact = 1. + delta * cp * gamma &
3403 & * to(i,k) / hvap
3404 aa1(i) = aa1(i) + &
3405 & dz1 * (g / (cp * to(i,k))) &
3406 & * dbyo(i,k) / (1. + gamma) &
3407 & * rfact
3408 val = 0.
3409 aa1(i)=aa1(i)+ &
3410 & dz1 * g * delta * &
3411 & max(val,(qeso(i,k) - qo(i,k)))
3412 endif
3413 endif
3414 enddo
3415 enddo
3416 do i = 1, im
3417 if(cnvflg(i).and.aa1(i).le.0.) cnvflg(i) = .false.
3418 enddo
3419 !!
3420 totflg = .true.
3421 do i=1,im
3422 totflg = totflg .and. (.not. cnvflg(i))
3423 enddo
3424 if(totflg) return
3425 !!
3426 !c
3427 !c estimate the onvective overshooting as the level
3428 !c where the [aafac * cloud work function] becomes zero,
3429 !c which is the final cloud top
3430 !c
3431 do i = 1, im
3432 if (cnvflg(i)) then
3433 aa2(i) = aafac * aa1(i)
3434 endif
3435 enddo
3436 !c
3437 do i = 1, im
3438 flg(i) = cnvflg(i)
3439 ktcon1(i) = kmax(i) - 1
3440 enddo
3441 do k = 2, km1
3442 do i = 1, im
3443 if (flg(i)) then
3444 if(k.ge.ktcon(i).and.k.lt.kmax(i)) then
3445 dz1 = zo(i,k+1) - zo(i,k)
3446 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
3447 rfact = 1. + delta * cp * gamma &
3448 & * to(i,k) / hvap
3449 aa2(i) = aa2(i) + &
3450 & dz1 * (g / (cp * to(i,k))) &
3451 & * dbyo(i,k) / (1. + gamma) &
3452 & * rfact
3453 if(aa2(i).lt.0.) then
3454 ktcon1(i) = k
3455 flg(i) = .false.
3456 endif
3457 endif
3458 endif
3459 enddo
3460 enddo
3461 !c
3462 !c compute cloud moisture property, detraining cloud water
3463 !c and precipitation in overshooting layers
3464 !c
3465 do k = 2, km1
3466 do i = 1, im
3467 if (cnvflg(i)) then
3468 if(k.ge.ktcon(i).and.k.lt.ktcon1(i)) then
3469 dz = zi(i,k) - zi(i,k-1)
3470 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
3471 qrch = qeso(i,k) &
3472 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
3473 !cj
3474 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
3475 tem1 = 0.5 * xlamud(i) * dz
3476 factor = 1. + tem - tem1
3477 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* &
3478 & (qo(i,k)+qo(i,k-1)))/factor
3479 !cj
3480 dq = eta(i,k) * (qcko(i,k) - qrch)
3481 !c
3482 !c check if there is excess moisture to release latent heat
3483 !c
3484 if(dq.gt.0.) then
3485 etah = .5 * (eta(i,k) + eta(i,k-1))
3486 if(ncloud.gt.0.) then
3487 dp = 1000. * del(i,k)
3488 qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
3489 dellal(i,k) = etah * c1 * dz * qlk * g / dp
3490 else
3491 qlk = dq / (eta(i,k) + etah * c0 * dz)
3492 endif
3493 qcko(i,k) = qlk + qrch
3494 pwo(i,k) = etah * c0 * dz * qlk
3495 pwavo(i) = pwavo(i) + pwo(i,k)
3496 endif
3497 endif
3498 endif
3499 enddo
3500 enddo
3501 !c
3502 !c exchange ktcon with ktcon1
3503 !c
3504 do i = 1, im
3505 if(cnvflg(i)) then
3506 kk = ktcon(i)
3507 ktcon(i) = ktcon1(i)
3508 ktcon1(i) = kk
3509 endif
3510 enddo
3511 !c
3512 !c this section is ready for cloud water
3513 !c
3514 if(ncloud.gt.0) then
3515 !c
3516 !c compute liquid and vapor separation at cloud top
3517 !c
3518 do i = 1, im
3519 if(cnvflg(i)) then
3520 k = ktcon(i) - 1
3521 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
3522 qrch = qeso(i,k) &
3523 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
3524 dq = qcko(i,k) - qrch
3525 !c
3526 !c check if there is excess moisture to release latent heat
3527 !c
3528 if(dq.gt.0.) then
3529 qlko_ktcon(i) = dq
3530 qcko(i,k) = qrch
3531 endif
3532 endif
3533 enddo
3534 endif
3535 !c
3536 !ccccc if(lat.eq.latd.and.lon.eq.lond.and.cnvflg(i)) then
3537 !ccccc print *, ' aa1(i) before dwndrft =', aa1(i)
3538 !ccccc endif
3539 !c
3540 !c------- downdraft calculations
3541 !c
3542 !c--- compute precipitation efficiency in terms of windshear
3543 !c
3544 do i = 1, im
3545 if(cnvflg(i)) then
3546 vshear(i) = 0.
3547 endif
3548 enddo
3549 do k = 2, km
3550 do i = 1, im
3551 if (cnvflg(i)) then
3552 if(k.gt.kb(i).and.k.le.ktcon(i)) then
3553 shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2 &
3554 & + (vo(i,k)-vo(i,k-1)) ** 2)
3555 vshear(i) = vshear(i) + shear
3556 endif
3557 endif
3558 enddo
3559 enddo
3560 do i = 1, im
3561 if(cnvflg(i)) then
3562 vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i))-zi(i,kb(i)))
3563 e1=1.591-.639*vshear(i) &
3564 & +.0953*(vshear(i)**2)-.00496*(vshear(i)**3)
3565 edt(i)=1.-e1
3566 val = .9
3567 edt(i) = min(edt(i),val)
3568 val = .0
3569 edt(i) = max(edt(i),val)
3570 edto(i)=edt(i)
3571 edtx(i)=edt(i)
3572 endif
3573 enddo
3574 !c
3575 !c determine detrainment rate between 1 and kbcon
3576 !c
3577 do i = 1, im
3578 if(cnvflg(i)) then
3579 sumx(i) = 0.
3580 endif
3581 enddo
3582 do k = 1, km1
3583 do i = 1, im
3584 if(cnvflg(i).and.k.ge.1.and.k.lt.kbcon(i)) then
3585 dz = zi(i,k+1) - zi(i,k)
3586 sumx(i) = sumx(i) + dz
3587 endif
3588 enddo
3589 enddo
3590 do i = 1, im
3591 beta = betas
3592 if(slimsk(i).eq.1.) beta = betal
3593 if(cnvflg(i)) then
3594 dz = (sumx(i)+zi(i,1))/float(kbcon(i))
3595 tem = 1./float(kbcon(i))
3596 xlamd(i) = (1.-beta**tem)/dz
3597 endif
3598 enddo
3599 !c
3600 !c determine downdraft mass flux
3601 !c
3602 do k = km1, 1, -1
3603 do i = 1, im
3604 if (cnvflg(i) .and. k .le. kmax(i)-1) then
3605 if(k.lt.jmin(i).and.k.ge.kbcon(i)) then
3606 dz = zi(i,k+1) - zi(i,k)
3607 ptem = xlamdd - xlamde
3608 etad(i,k) = etad(i,k+1) * (1. - ptem * dz)
3609 else if(k.lt.kbcon(i)) then
3610 dz = zi(i,k+1) - zi(i,k)
3611 ptem = xlamd(i) + xlamdd - xlamde
3612 etad(i,k) = etad(i,k+1) * (1. - ptem * dz)
3613 endif
3614 endif
3615 enddo
3616 enddo
3617 !c
3618 !c--- downdraft moisture properties
3619 !c
3620 do i = 1, im
3621 if(cnvflg(i)) then
3622 jmn = jmin(i)
3623 hcdo(i,jmn) = heo(i,jmn)
3624 qcdo(i,jmn) = qo(i,jmn)
3625 qrcdo(i,jmn)= qeso(i,jmn)
3626 ucdo(i,jmn) = uo(i,jmn)
3627 vcdo(i,jmn) = vo(i,jmn)
3628 pwevo(i) = 0.
3629 endif
3630 enddo
3631 !cj
3632 do k = km1, 1, -1
3633 do i = 1, im
3634 if (cnvflg(i) .and. k.lt.jmin(i)) then
3635 dz = zi(i,k+1) - zi(i,k)
3636 if(k.ge.kbcon(i)) then
3637 tem = xlamde * dz
3638 tem1 = 0.5 * xlamdd * dz
3639 else
3640 tem = xlamde * dz
3641 tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
3642 endif
3643 factor = 1. + tem - tem1
3644 ptem = 0.5 * tem - pgcon
3645 ptem1= 0.5 * tem + pgcon
3646 hcdo(i,k) = ((1.-tem1)*hcdo(i,k+1)+tem*0.5* &
3647 & (heo(i,k)+heo(i,k+1)))/factor
3648 ucdo(i,k) = ((1.-tem1)*ucdo(i,k+1)+ptem*uo(i,k+1) &
3649 & +ptem1*uo(i,k))/factor
3650 vcdo(i,k) = ((1.-tem1)*vcdo(i,k+1)+ptem*vo(i,k+1) &
3651 & +ptem1*vo(i,k))/factor
3652 dbyo(i,k) = hcdo(i,k) - heso(i,k)
3653 endif
3654 enddo
3655 enddo
3656 !c
3657 do k = km1, 1, -1
3658 do i = 1, im
3659 if (cnvflg(i).and.k.lt.jmin(i)) then
3660 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
3661 qrcdo(i,k) = qeso(i,k)+ &
3662 & (1./hvap)*(gamma/(1.+gamma))*dbyo(i,k)
3663 ! detad = etad(i,k+1) - etad(i,k)
3664 !cj
3665 dz = zi(i,k+1) - zi(i,k)
3666 if(k.ge.kbcon(i)) then
3667 tem = xlamde * dz
3668 tem1 = 0.5 * xlamdd * dz
3669 else
3670 tem = xlamde * dz
3671 tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
3672 endif
3673 factor = 1. + tem - tem1
3674 qcdo(i,k) = ((1.-tem1)*qcdo(i,k+1)+tem*0.5* &
3675 & (qo(i,k)+qo(i,k+1)))/factor
3676 !cj
3677 ! pwdo(i,k) = etad(i,k+1) * qcdo(i,k+1) -
3678 ! & etad(i,k) * qrcdo(i,k)
3679 ! pwdo(i,k) = pwdo(i,k) - detad *
3680 ! & .5 * (qrcdo(i,k) + qrcdo(i,k+1))
3681 !cj
3682 pwdo(i,k) = etad(i,k+1) * (qcdo(i,k) - qrcdo(i,k))
3683 qcdo(i,k) = qrcdo(i,k)
3684 pwevo(i) = pwevo(i) + pwdo(i,k)
3685 endif
3686 enddo
3687 enddo
3688 !c
3689 !c--- final downdraft strength dependent on precip
3690 !c--- efficiency (edt), normalized condensate (pwav), and
3691 !c--- evaporate (pwev)
3692 !c
3693 do i = 1, im
3694 edtmax = edtmaxl
3695 if(slimsk(i).eq.0.) edtmax = edtmaxs
3696 if(cnvflg(i)) then
3697 if(pwevo(i).lt.0.) then
3698 edto(i) = -edto(i) * pwavo(i) / pwevo(i)
3699 edto(i) = min(edto(i),edtmax)
3700 else
3701 edto(i) = 0.
3702 endif
3703 endif
3704 enddo
3705 !c
3706 !c--- downdraft cloudwork functions
3707 !c
3708 do k = km1, 1, -1
3709 do i = 1, im
3710 if (cnvflg(i) .and. k .lt. jmin(i)) then
3711 gamma = el2orc * qeso(i,k) / to(i,k)**2
3712 dhh=hcdo(i,k)
3713 dt=to(i,k)
3714 dg=gamma
3715 dh=heso(i,k)
3716 dz=-1.*(zo(i,k+1)-zo(i,k))
3717 aa1(i)=aa1(i)+edto(i)*dz*(g/(cp*dt))*((dhh-dh)/(1.+dg)) &
3718 & *(1.+delta*cp*dg*dt/hvap)
3719 val=0.
3720 aa1(i)=aa1(i)+edto(i)* &
3721 & dz*g*delta*max(val,(qeso(i,k)-qo(i,k)))
3722 endif
3723 enddo
3724 enddo
3725 do i = 1, im
3726 if(cnvflg(i).and.aa1(i).le.0.) then
3727 cnvflg(i) = .false.
3728 endif
3729 enddo
3730 !!
3731 totflg = .true.
3732 do i=1,im
3733 totflg = totflg .and. (.not. cnvflg(i))
3734 enddo
3735 if(totflg) return
3736 !!
3737 !c
3738 !c--- what would the change be, that a cloud with unit mass
3739 !c--- will do to the environment?
3740 !c
3741 do k = 1, km
3742 do i = 1, im
3743 if(cnvflg(i) .and. k .le. kmax(i)) then
3744 dellah(i,k) = 0.
3745 dellaq(i,k) = 0.
3746 dellau(i,k) = 0.
3747 dellav(i,k) = 0.
3748 endif
3749 enddo
3750 enddo
3751 do i = 1, im
3752 if(cnvflg(i)) then
3753 dp = 1000. * del(i,1)
3754 dellah(i,1) = edto(i) * etad(i,1) * (hcdo(i,1) &
3755 & - heo(i,1)) * g / dp
3756 dellaq(i,1) = edto(i) * etad(i,1) * (qcdo(i,1) &
3757 & - qo(i,1)) * g / dp
3758 dellau(i,1) = edto(i) * etad(i,1) * (ucdo(i,1) &
3759 & - uo(i,1)) * g / dp
3760 dellav(i,1) = edto(i) * etad(i,1) * (vcdo(i,1) &
3761 & - vo(i,1)) * g / dp
3762 endif
3763 enddo
3764 !c
3765 !c--- changed due to subsidence and entrainment
3766 !c
3767 do k = 2, km1
3768 do i = 1, im
3769 if (cnvflg(i).and.k.lt.ktcon(i)) then
3770 aup = 1.
3771 if(k.le.kb(i)) aup = 0.
3772 adw = 1.
3773 if(k.gt.jmin(i)) adw = 0.
3774 dp = 1000. * del(i,k)
3775 dz = zi(i,k) - zi(i,k-1)
3776 !c
3777 dv1h = heo(i,k)
3778 dv2h = .5 * (heo(i,k) + heo(i,k-1))
3779 dv3h = heo(i,k-1)
3780 dv1q = qo(i,k)
3781 dv2q = .5 * (qo(i,k) + qo(i,k-1))
3782 dv3q = qo(i,k-1)
3783 dv1u = uo(i,k)
3784 dv2u = .5 * (uo(i,k) + uo(i,k-1))
3785 dv3u = uo(i,k-1)
3786 dv1v = vo(i,k)
3787 dv2v = .5 * (vo(i,k) + vo(i,k-1))
3788 dv3v = vo(i,k-1)
3789 !c
3790 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1))
3791 tem1 = xlamud(i)
3792 !c
3793 if(k.le.kbcon(i)) then
3794 ptem = xlamde
3795 ptem1 = xlamd(i)+xlamdd
3796 else
3797 ptem = xlamde
3798 ptem1 = xlamdd
3799 endif
3800 !cj
3801 dellah(i,k) = dellah(i,k) + &
3802 & ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1h &
3803 & - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3h &
3804 & - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2h*dz &
3805 & + aup*tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz &
3806 & + adw*edto(i)*ptem1*etad(i,k)*.5*(hcdo(i,k)+hcdo(i,k-1)) &
3807 & *dz) *g/dp
3808 !cj
3809 dellaq(i,k) = dellaq(i,k) + &
3810 & ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1q &
3811 & - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3q &
3812 & - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2q*dz &
3813 & + aup*tem1*eta(i,k-1)*.5*(qcko(i,k)+qcko(i,k-1))*dz &
3814 & + adw*edto(i)*ptem1*etad(i,k)*.5*(qrcdo(i,k)+qrcdo(i,k-1)) &
3815 & *dz) *g/dp
3816 !23456789012345678901234567890123456789012345678901234567890123456789012
3817 !cj
3818 dellau(i,k) = dellau(i,k) + &
3819 & ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1u &
3820 & - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3u &
3821 & - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2u*dz &
3822 & + aup*tem1*eta(i,k-1)*.5*(ucko(i,k)+ucko(i,k-1))*dz &
3823 & + adw*edto(i)*ptem1*etad(i,k)*.5*(ucdo(i,k)+ucdo(i,k-1))*dz &
3824 & - pgcon*(aup*eta(i,k-1)-adw*edto(i)*etad(i,k))*(dv1u-dv3u) &
3825 & ) *g/dp
3826 !cj
3827 dellav(i,k) = dellav(i,k) + &
3828 & ((aup*eta(i,k)-adw*edto(i)*etad(i,k))*dv1v &
3829 & - (aup*eta(i,k-1)-adw*edto(i)*etad(i,k-1))*dv3v &
3830 & - (aup*tem*eta(i,k-1)+adw*edto(i)*ptem*etad(i,k))*dv2v*dz &
3831 & + aup*tem1*eta(i,k-1)*.5*(vcko(i,k)+vcko(i,k-1))*dz &
3832 & + adw*edto(i)*ptem1*etad(i,k)*.5*(vcdo(i,k)+vcdo(i,k-1))*dz &
3833 & - pgcon*(aup*eta(i,k-1)-adw*edto(i)*etad(i,k))*(dv1v-dv3v) &
3834 & ) *g/dp
3835 !cj
3836 endif
3837 enddo
3838 enddo
3839 !c
3840 !c------- cloud top
3841 !c
3842 do i = 1, im
3843 if(cnvflg(i)) then
3844 indx = ktcon(i)
3845 dp = 1000. * del(i,indx)
3846 dv1h = heo(i,indx-1)
3847 dellah(i,indx) = eta(i,indx-1) * &
3848 & (hcko(i,indx-1) - dv1h) * g / dp
3849 dv1q = qo(i,indx-1)
3850 dellaq(i,indx) = eta(i,indx-1) * &
3851 & (qcko(i,indx-1) - dv1q) * g / dp
3852 dv1u = uo(i,indx-1)
3853 dellau(i,indx) = eta(i,indx-1) * &
3854 & (ucko(i,indx-1) - dv1u) * g / dp
3855 dv1v = vo(i,indx-1)
3856 dellav(i,indx) = eta(i,indx-1) * &
3857 & (vcko(i,indx-1) - dv1v) * g / dp
3858 !c
3859 !c cloud water
3860 !c
3861 dellal(i,indx) = eta(i,indx-1) * &
3862 & qlko_ktcon(i) * g / dp
3863 endif
3864 enddo
3865 !c
3866 !c------- final changed variable per unit mass flux
3867 !c
3868 do k = 1, km
3869 do i = 1, im
3870 if (cnvflg(i).and.k .le. kmax(i)) then
3871 if(k.gt.ktcon(i)) then
3872 qo(i,k) = q1(i,k)
3873 to(i,k) = t1(i,k)
3874 endif
3875 if(k.le.ktcon(i)) then
3876 qo(i,k) = dellaq(i,k) * mbdt + q1(i,k)
3877 dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp
3878 to(i,k) = dellat * mbdt + t1(i,k)
3879 val = 1.e-10
3880 qo(i,k) = max(qo(i,k), val )
3881 endif
3882 endif
3883 enddo
3884 enddo
3885 !c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3886 !c
3887 !c--- the above changed environment is now used to calulate the
3888 !c--- effect the arbitrary cloud (with unit mass flux)
3889 !c--- would have on the stability,
3890 !c--- which then is used to calculate the real mass flux,
3891 !c--- necessary to keep this change in balance with the large-scale
3892 !c--- destabilization.
3893 !c
3894 !c--- environmental conditions again, first heights
3895 !c
3896 do k = 1, km
3897 do i = 1, im
3898 if(cnvflg(i) .and. k .le. kmax(i)) then
3899 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
3900 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k)+epsm1*qeso(i,k))
3901 val = 1.e-8
3902 qeso(i,k) = max(qeso(i,k), val )
3903 ! tvo(i,k) = to(i,k) + delta * to(i,k) * qo(i,k)
3904 endif
3905 enddo
3906 enddo
3907 !c
3908 !c--- moist static energy
3909 !c
3910 do k = 1, km1
3911 do i = 1, im
3912 if(cnvflg(i) .and. k .le. kmax(i)-1) then
3913 dz = .5 * (zo(i,k+1) - zo(i,k))
3914 dp = .5 * (pfld(i,k+1) - pfld(i,k))
3915 es = 0.01 * fpvs(to(i,k+1)) ! fpvs is in pa
3916 pprime = pfld(i,k+1) + epsm1 * es
3917 qs = eps * es / pprime
3918 dqsdp = - qs / pprime
3919 desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
3920 dqsdt = qs * pfld(i,k+1) * desdt / (es * pprime)
3921 gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
3922 dt = (g * dz + hvap * dqsdp * dp) / (cp * (1. + gamma))
3923 dq = dqsdt * dt + dqsdp * dp
3924 to(i,k) = to(i,k+1) + dt
3925 qo(i,k) = qo(i,k+1) + dq
3926 po(i,k) = .5 * (pfld(i,k) + pfld(i,k+1))
3927 endif
3928 enddo
3929 enddo
3930 do k = 1, km1
3931 do i = 1, im
3932 if(cnvflg(i) .and. k .le. kmax(i)-1) then
3933 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
3934 qeso(i,k) = eps * qeso(i,k) / (po(i,k) + epsm1 * qeso(i,k))
3935 val1 = 1.e-8
3936 qeso(i,k) = max(qeso(i,k), val1)
3937 val2 = 1.e-10
3938 qo(i,k) = max(qo(i,k), val2 )
3939 ! qo(i,k) = min(qo(i,k),qeso(i,k))
3940 heo(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) + &
3941 & cp * to(i,k) + hvap * qo(i,k)
3942 heso(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) + &
3943 & cp * to(i,k) + hvap * qeso(i,k)
3944 endif
3945 enddo
3946 enddo
3947 do i = 1, im
3948 if(cnvflg(i)) then
3949 k = kmax(i)
3950 heo(i,k) = g * zo(i,k) + cp * to(i,k) + hvap * qo(i,k)
3951 heso(i,k) = g * zo(i,k) + cp * to(i,k) + hvap * qeso(i,k)
3952 !c heo(i,k) = min(heo(i,k),heso(i,k))
3953 endif
3954 enddo
3955 !c
3956 !c**************************** static control
3957 !c
3958 !c------- moisture and cloud work functions
3959 !c
3960 do i = 1, im
3961 if(cnvflg(i)) then
3962 xaa0(i) = 0.
3963 xpwav(i) = 0.
3964 endif
3965 enddo
3966 !c
3967 do i = 1, im
3968 if(cnvflg(i)) then
3969 indx = kb(i)
3970 hcko(i,indx) = heo(i,indx)
3971 qcko(i,indx) = qo(i,indx)
3972 endif
3973 enddo
3974 do k = 2, km1
3975 do i = 1, im
3976 if (cnvflg(i)) then
3977 if(k.gt.kb(i).and.k.le.ktcon(i)) then
3978 dz = zi(i,k) - zi(i,k-1)
3979 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
3980 tem1 = 0.5 * xlamud(i) * dz
3981 factor = 1. + tem - tem1
3982 hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5* &
3983 & (heo(i,k)+heo(i,k-1)))/factor
3984 endif
3985 endif
3986 enddo
3987 enddo
3988 do k = 2, km1
3989 do i = 1, im
3990 if (cnvflg(i)) then
3991 if(k.gt.kb(i).and.k.lt.ktcon(i)) then
3992 dz = zi(i,k) - zi(i,k-1)
3993 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
3994 xdby = hcko(i,k) - heso(i,k)
3995 xqrch = qeso(i,k) &
3996 & + gamma * xdby / (hvap * (1. + gamma))
3997 !cj
3998 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
3999 tem1 = 0.5 * xlamud(i) * dz
4000 factor = 1. + tem - tem1
4001 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* &
4002 & (qo(i,k)+qo(i,k-1)))/factor
4003 !cj
4004 dq = eta(i,k) * (qcko(i,k) - xqrch)
4005 !c
4006 if(k.ge.kbcon(i).and.dq.gt.0.) then
4007 etah = .5 * (eta(i,k) + eta(i,k-1))
4008 if(ncloud.gt.0..and.k.gt.jmin(i)) then
4009 qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
4010 else
4011 qlk = dq / (eta(i,k) + etah * c0 * dz)
4012 endif
4013 if(k.lt.ktcon1(i)) then
4014 xaa0(i) = xaa0(i) - dz * g * qlk
4015 endif
4016 qcko(i,k) = qlk + xqrch
4017 xpw = etah * c0 * dz * qlk
4018 xpwav(i) = xpwav(i) + xpw
4019 endif
4020 endif
4021 if(k.ge.kbcon(i).and.k.lt.ktcon1(i)) then
4022 dz1 = zo(i,k+1) - zo(i,k)
4023 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
4024 rfact = 1. + delta * cp * gamma &
4025 & * to(i,k) / hvap
4026 xaa0(i) = xaa0(i) &
4027 & + dz1 * (g / (cp * to(i,k))) &
4028 & * xdby / (1. + gamma) &
4029 & * rfact
4030 val=0.
4031 xaa0(i)=xaa0(i)+ &
4032 & dz1 * g * delta * &
4033 & max(val,(qeso(i,k) - qo(i,k)))
4034 endif
4035 endif
4036 enddo
4037 enddo
4038 !c
4039 !c------- downdraft calculations
4040 !c
4041 !c--- downdraft moisture properties
4042 !c
4043 do i = 1, im
4044 if(cnvflg(i)) then
4045 jmn = jmin(i)
4046 hcdo(i,jmn) = heo(i,jmn)
4047 qcdo(i,jmn) = qo(i,jmn)
4048 qrcd(i,jmn) = qeso(i,jmn)
4049 xpwev(i) = 0.
4050 endif
4051 enddo
4052 !cj
4053 do k = km1, 1, -1
4054 do i = 1, im
4055 if (cnvflg(i) .and. k.lt.jmin(i)) then
4056 dz = zi(i,k+1) - zi(i,k)
4057 if(k.ge.kbcon(i)) then
4058 tem = xlamde * dz
4059 tem1 = 0.5 * xlamdd * dz
4060 else
4061 tem = xlamde * dz
4062 tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
4063 endif
4064 factor = 1. + tem - tem1
4065 hcdo(i,k) = ((1.-tem1)*hcdo(i,k+1)+tem*0.5* &
4066 & (heo(i,k)+heo(i,k+1)))/factor
4067 endif
4068 enddo
4069 enddo
4070 !cj
4071 do k = km1, 1, -1
4072 do i = 1, im
4073 if (cnvflg(i) .and. k .lt. jmin(i)) then
4074 dq = qeso(i,k)
4075 dt = to(i,k)
4076 gamma = el2orc * dq / dt**2
4077 dh = hcdo(i,k) - heso(i,k)
4078 qrcd(i,k)=dq+(1./hvap)*(gamma/(1.+gamma))*dh
4079 ! detad = etad(i,k+1) - etad(i,k)
4080 !cj
4081 dz = zi(i,k+1) - zi(i,k)
4082 if(k.ge.kbcon(i)) then
4083 tem = xlamde * dz
4084 tem1 = 0.5 * xlamdd * dz
4085 else
4086 tem = xlamde * dz
4087 tem1 = 0.5 * (xlamd(i)+xlamdd) * dz
4088 endif
4089 factor = 1. + tem - tem1
4090 qcdo(i,k) = ((1.-tem1)*qcdo(i,k+1)+tem*0.5* &
4091 & (qo(i,k)+qo(i,k+1)))/factor
4092 !cj
4093 ! xpwd = etad(i,k+1) * qcdo(i,k+1) -
4094 ! & etad(i,k) * qrcd(i,k)
4095 ! xpwd = xpwd - detad *
4096 ! & .5 * (qrcd(i,k) + qrcd(i,k+1))
4097 !cj
4098 xpwd = etad(i,k+1) * (qcdo(i,k) - qrcd(i,k))
4099 qcdo(i,k)= qrcd(i,k)
4100 xpwev(i) = xpwev(i) + xpwd
4101 endif
4102 enddo
4103 enddo
4104 !c
4105 do i = 1, im
4106 edtmax = edtmaxl
4107 if(slimsk(i).eq.0.) edtmax = edtmaxs
4108 if(cnvflg(i)) then
4109 if(xpwev(i).ge.0.) then
4110 edtx(i) = 0.
4111 else
4112 edtx(i) = -edtx(i) * xpwav(i) / xpwev(i)
4113 edtx(i) = min(edtx(i),edtmax)
4114 endif
4115 endif
4116 enddo
4117 !c
4118 !c
4119 !c--- downdraft cloudwork functions
4120 !c
4121 !c
4122 do k = km1, 1, -1
4123 do i = 1, im
4124 if (cnvflg(i) .and. k.lt.jmin(i)) then
4125 gamma = el2orc * qeso(i,k) / to(i,k)**2
4126 dhh=hcdo(i,k)
4127 dt= to(i,k)
4128 dg= gamma
4129 dh= heso(i,k)
4130 dz=-1.*(zo(i,k+1)-zo(i,k))
4131 xaa0(i)=xaa0(i)+edtx(i)*dz*(g/(cp*dt))*((dhh-dh)/(1.+dg)) &
4132 & *(1.+delta*cp*dg*dt/hvap)
4133 val=0.
4134 xaa0(i)=xaa0(i)+edtx(i)* &
4135 & dz*g*delta*max(val,(qeso(i,k)-qo(i,k)))
4136 endif
4137 enddo
4138 enddo
4139 !c
4140 !c calculate critical cloud work function
4141 !c
4142 do i = 1, im
4143 if(cnvflg(i)) then
4144 if(pfld(i,ktcon(i)).lt.pcrit(15))then
4145 acrt(i)=acrit(15)*(975.-pfld(i,ktcon(i))) &
4146 & /(975.-pcrit(15))
4147 else if(pfld(i,ktcon(i)).gt.pcrit(1))then
4148 acrt(i)=acrit(1)
4149 else
4150 k = int((850. - pfld(i,ktcon(i)))/50.) + 2
4151 k = min(k,15)
4152 k = max(k,2)
4153 acrt(i)=acrit(k)+(acrit(k-1)-acrit(k))* &
4154 & (pfld(i,ktcon(i))-pcrit(k))/(pcrit(k-1)-pcrit(k))
4155 endif
4156 endif
4157 enddo
4158 do i = 1, im
4159 if(cnvflg(i)) then
4160 if(slimsk(i).eq.1.) then
4161 w1 = w1l
4162 w2 = w2l
4163 w3 = w3l
4164 w4 = w4l
4165 else
4166 w1 = w1s
4167 w2 = w2s
4168 w3 = w3s
4169 w4 = w4s
4170 endif
4171 !c
4172 !c modify critical cloud workfunction by cloud base vertical velocity
4173 !c
4174 if(pdot(i).le.w4) then
4175 acrtfct(i) = (pdot(i) - w4) / (w3 - w4)
4176 elseif(pdot(i).ge.-w4) then
4177 acrtfct(i) = - (pdot(i) + w4) / (w4 - w3)
4178 else
4179 acrtfct(i) = 0.
4180 endif
4181 val1 = -1.
4182 acrtfct(i) = max(acrtfct(i),val1)
4183 val2 = 1.
4184 acrtfct(i) = min(acrtfct(i),val2)
4185 acrtfct(i) = 1. - acrtfct(i)
4186 !c
4187 !c modify acrtfct(i) by colume mean rh if rhbar(i) is greater than 80 percent
4188 !c
4189 !c if(rhbar(i).ge..8) then
4190 !c acrtfct(i) = acrtfct(i) * (.9 - min(rhbar(i),.9)) * 10.
4191 !c endif
4192 !c
4193 !c modify adjustment time scale by cloud base vertical velocity
4194 !c
4195 val1=0.
4196 dtconv(i) = dt2 + max((1800. - dt2),val1) * &
4197 & (pdot(i) - w2) / (w1 - w2)
4198 !c dtconv(i) = max(dtconv(i), dt2)
4199 !c dtconv(i) = 1800. * (pdot(i) - w2) / (w1 - w2)
4200 dtconv(i) = max(dtconv(i),dtmin)
4201 dtconv(i) = min(dtconv(i),dtmax)
4202 !c
4203 endif
4204 enddo
4205 !c
4206 !c--- large scale forcing
4207 !c
4208 xmbmx1=-1.e20
4209 do i= 1, im
4210 if(cnvflg(i)) then
4211 fld(i)=(aa1(i)-acrt(i)* acrtfct(i))/dtconv(i)
4212 if(fld(i).le.0.) cnvflg(i) = .false.
4213 endif
4214 if(cnvflg(i)) then
4215 !c xaa0(i) = max(xaa0(i),0.)
4216 xk(i) = (xaa0(i) - aa1(i)) / mbdt
4217 if(xk(i).ge.0.) cnvflg(i) = .false.
4218 endif
4219 !c
4220 !c--- kernel, cloud base mass flux
4221 !c
4222 if(cnvflg(i)) then
4223 xmb(i) = -fld(i) / xk(i)
4224 xmb(i) = min(xmb(i),xmbmax(i))
4225 xmbmx1=max(xmbmx1,xmb(i))
4226 endif
4227 enddo
4228 ! if(xmbmx1.gt.0.4)print*,'qingfu test xmbmx1=',xmbmx1
4229 !!
4230 totflg = .true.
4231 do i=1,im
4232 totflg = totflg .and. (.not. cnvflg(i))
4233 enddo
4234 if(totflg) return
4235 !!
4236 !c
4237 !c restore to,qo,uo,vo to t1,q1,u1,v1 in case convection stops
4238 !c
4239 do k = 1, km
4240 do i = 1, im
4241 if (cnvflg(i) .and. k .le. kmax(i)) then
4242 to(i,k) = t1(i,k)
4243 qo(i,k) = q1(i,k)
4244 uo(i,k) = u1(i,k)
4245 vo(i,k) = v1(i,k)
4246 qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
4247 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
4248 val = 1.e-8
4249 qeso(i,k) = max(qeso(i,k), val )
4250 endif
4251 enddo
4252 enddo
4253 !c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4254 !c
4255 !c--- feedback: simply the changes from the cloud with unit mass flux
4256 !c--- multiplied by the mass flux necessary to keep the
4257 !c--- equilibrium with the larger-scale.
4258 !c
4259 do i = 1, im
4260 delhbar(i) = 0.
4261 delqbar(i) = 0.
4262 deltbar(i) = 0.
4263 delubar(i) = 0.
4264 delvbar(i) = 0.
4265 qcond(i) = 0.
4266 enddo
4267 do k = 1, km
4268 do i = 1, im
4269 if (cnvflg(i) .and. k .le. kmax(i)) then
4270 if(k.le.ktcon(i)) then
4271 dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp
4272 t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2
4273 q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2
4274 tem = 1./rcs(i)
4275 u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 * tem
4276 v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 * tem
4277 dp = 1000. * del(i,k)
4278 delhbar(i) = delhbar(i) + dellah(i,k)*xmb(i)*dp/g
4279 delqbar(i) = delqbar(i) + dellaq(i,k)*xmb(i)*dp/g
4280 deltbar(i) = deltbar(i) + dellat*xmb(i)*dp/g
4281 delubar(i) = delubar(i) + dellau(i,k)*xmb(i)*dp/g
4282 delvbar(i) = delvbar(i) + dellav(i,k)*xmb(i)*dp/g
4283 endif
4284 endif
4285 enddo
4286 enddo
4287 do k = 1, km
4288 do i = 1, im
4289 if (cnvflg(i) .and. k .le. kmax(i)) then
4290 if(k.le.ktcon(i)) then
4291 qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
4292 qeso(i,k) = eps * qeso(i,k)/(pfld(i,k) + epsm1*qeso(i,k))
4293 val = 1.e-8
4294 qeso(i,k) = max(qeso(i,k), val )
4295 endif
4296 endif
4297 enddo
4298 enddo
4299 !c
4300 do i = 1, im
4301 rntot(i) = 0.
4302 delqev(i) = 0.
4303 delq2(i) = 0.
4304 flg(i) = cnvflg(i)
4305 enddo
4306 do k = km, 1, -1
4307 do i = 1, im
4308 if (cnvflg(i) .and. k .le. kmax(i)) then
4309 if(k.lt.ktcon(i)) then
4310 aup = 1.
4311 if(k.le.kb(i)) aup = 0.
4312 adw = 1.
4313 if(k.ge.jmin(i)) adw = 0.
4314 rain = aup * pwo(i,k) + adw * edto(i) * pwdo(i,k)
4315 rntot(i) = rntot(i) + rain * xmb(i) * .001 * dt2
4316 endif
4317 endif
4318 enddo
4319 enddo
4320 do k = km, 1, -1
4321 do i = 1, im
4322 if (k .le. kmax(i)) then
4323 deltv(i) = 0.
4324 delq(i) = 0.
4325 qevap(i) = 0.
4326 if(cnvflg(i).and.k.lt.ktcon(i)) then
4327 aup = 1.
4328 if(k.le.kb(i)) aup = 0.
4329 adw = 1.
4330 if(k.ge.jmin(i)) adw = 0.
4331 rain = aup * pwo(i,k) + adw * edto(i) * pwdo(i,k)
4332 rn(i) = rn(i) + rain * xmb(i) * .001 * dt2
4333 endif
4334 if(flg(i).and.k.lt.ktcon(i)) then
4335 evef = edt(i) * evfact
4336 if(slimsk(i).eq.1.) evef=edt(i) * evfactl
4337 ! if(slimsk(i).eq.1.) evef=.07
4338 !c if(slimsk(i).ne.1.) evef = 0.
4339 qcond(i) = evef * (q1(i,k) - qeso(i,k)) &
4340 & / (1. + el2orc * qeso(i,k) / t1(i,k)**2)
4341 dp = 1000. * del(i,k)
4342 if(rn(i).gt.0..and.qcond(i).lt.0.) then
4343 qevap(i) = -qcond(i) * (1.-exp(-.32*sqrt(dt2*rn(i))))
4344 qevap(i) = min(qevap(i), rn(i)*1000.*g/dp)
4345 delq2(i) = delqev(i) + .001 * qevap(i) * dp / g
4346 endif
4347 if(rn(i).gt.0..and.qcond(i).lt.0..and. &
4348 & delq2(i).gt.rntot(i)) then
4349 qevap(i) = 1000.* g * (rntot(i) - delqev(i)) / dp
4350 flg(i) = .false.
4351 endif
4352 if(rn(i).gt.0..and.qevap(i).gt.0.) then
4353 q1(i,k) = q1(i,k) + qevap(i)
4354 t1(i,k) = t1(i,k) - elocp * qevap(i)
4355 rn(i) = rn(i) - .001 * qevap(i) * dp / g
4356 deltv(i) = - elocp*qevap(i)/dt2
4357 delq(i) = + qevap(i)/dt2
4358 delqev(i) = delqev(i) + .001*dp*qevap(i)/g
4359 endif
4360 dellaq(i,k) = dellaq(i,k) + delq(i) / xmb(i)
4361 delqbar(i) = delqbar(i) + delq(i)*dp/g
4362 deltbar(i) = deltbar(i) + deltv(i)*dp/g
4363 endif
4364 endif
4365 enddo
4366 enddo
4367 !cj
4368 ! do i = 1, im
4369 ! if(me.eq.31.and.cnvflg(i)) then
4370 ! if(cnvflg(i)) then
4371 ! print *, ' deep delhbar, delqbar, deltbar = ',
4372 ! & delhbar(i),hvap*delqbar(i),cp*deltbar(i)
4373 ! print *, ' deep delubar, delvbar = ',delubar(i),delvbar(i)
4374 ! print *, ' precip =', hvap*rn(i)*1000./dt2
4375 ! print*,'pdif= ',pfld(i,kbcon(i))-pfld(i,ktcon(i))
4376 ! endif
4377 ! enddo
4378 !c
4379 !c precipitation rate converted to actual precip
4380 !c in unit of m instead of kg
4381 !c
4382 do i = 1, im
4383 if(cnvflg(i)) then
4384 !c
4385 !c in the event of upper level rain evaporation and lower level downdraft
4386 !c moistening, rn can become negative, in this case, we back out of the
4387 !c heating and the moistening
4388 !c
4389 if(rn(i).lt.0..and..not.flg(i)) rn(i) = 0.
4390 if(rn(i).le.0.) then
4391 rn(i) = 0.
4392 else
4393 ktop(i) = ktcon(i)
4394 kbot(i) = kbcon(i)
4395 kcnv(i) = 1
4396 cldwrk(i) = aa1(i)
4397 endif
4398 endif
4399 enddo
4400 !c
4401 !c cloud water
4402 !c
4403 if (ncloud.gt.0) then
4404 !
4405 val1=1.0
4406 val2=0.0
4407 do k = 1, km
4408 do i = 1, im
4409 if (cnvflg(i) .and. rn(i).gt.0.) then
4410 if (k.gt.kb(i).and.k.le.ktcon(i)) then
4411 tem = dellal(i,k) * xmb(i) * dt2
4412 tem1 = max(val2, min(val1, (tcr-t1(i,k))*tcrf))
4413 if (ql(i,k,2) .gt. -999.0) then
4414 ql(i,k,1) = ql(i,k,1) + tem * tem1 ! ice
4415 ql(i,k,2) = ql(i,k,2) + tem *(1.0-tem1) ! water
4416 else
4417 ql(i,k,1) = ql(i,k,1) + tem
4418 endif
4419 endif
4420 endif
4421 enddo
4422 enddo
4423 !
4424 endif
4425 !c
4426 do k = 1, km
4427 do i = 1, im
4428 if(cnvflg(i).and.rn(i).le.0.) then
4429 if (k .le. kmax(i)) then
4430 t1(i,k) = to(i,k)
4431 q1(i,k) = qo(i,k)
4432 u1(i,k) = uo(i,k)
4433 v1(i,k) = vo(i,k)
4434 endif
4435 endif
4436 enddo
4437 enddo
4438 !
4439 ! hchuang code change
4440 !
4441 ! do k = 1, km
4442 ! do i = 1, im
4443 ! if(cnvflg(i).and.rn(i).gt.0.) then
4444 ! if(k.ge.kb(i) .and. k.lt.ktop(i)) then
4445 ! ud_mf(i,k) = eta(i,k) * xmb(i) * dt2
4446 ! endif
4447 ! endif
4448 ! enddo
4449 ! enddo
4450 ! do i = 1, im
4451 ! if(cnvflg(i).and.rn(i).gt.0.) then
4452 ! k = ktop(i)-1
4453 ! dt_mf(i,k) = ud_mf(i,k)
4454 ! endif
4455 ! enddo
4456 ! do k = 1, km
4457 ! do i = 1, im
4458 ! if(cnvflg(i).and.rn(i).gt.0.) then
4459 ! if(k.ge.1 .and. k.le.jmin(i)) then
4460 ! dd_mf(i,k) = edto(i) * etad(i,k) * xmb(i) * dt2
4461 ! endif
4462 ! endif
4463 ! enddo
4464 ! enddo
4465 !!
4466 return
4467 end subroutine sascnvn
4469 subroutine shalcnv(im,ix,km,jcap,delt,del,prsl,ps,phil,ql, &
4470 & q1,t1,u1,v1,rcs,rn,kbot,ktop,kcnv,slimsk, &
4471 & dot,ncloud,hpbl,heat,evap,pgcon)
4472 !
4473 use MODULE_GFS_machine , only : kind_phys
4474 use MODULE_GFS_funcphys , only : fpvs
4475 use MODULE_GFS_physcons, grav => con_g, cp => con_cp, hvap => con_hvap &
4476 &, rv => con_rv, fv => con_fvirt, t0c => con_t0c &
4477 &, rd => con_rd, cvap => con_cvap, cliq => con_cliq &
4478 &, eps => con_eps, epsm1 => con_epsm1
4479 implicit none
4480 !
4481 integer im, ix, km, jcap, ncloud, &
4482 & kbot(im), ktop(im), kcnv(im)
4483 real(kind=kind_phys) delt
4484 real(kind=kind_phys) ps(im), del(ix,km), prsl(ix,km), &
4485 & ql(ix,km,2),q1(ix,km), t1(ix,km), &
4486 & u1(ix,km), v1(ix,km), rcs(im), &
4487 & rn(im), slimsk(im), &
4488 & dot(ix,km), phil(ix,km), hpbl(im), &
4489 & heat(im), evap(im)
4490 ! &, ud_mf(im,km),dt_mf(im,km)
4492 real ud_mf(im,km),dt_mf(im,km)
4493 !
4494 integer i,j,indx, jmn, k, kk, latd, lond, km1
4495 integer kpbl(im)
4496 !
4497 real(kind=kind_phys) c0, cpoel, dellat, delta, &
4498 & desdt, deta, detad, dg, &
4499 & dh, dhh, dlnsig, dp, &
4500 & dq, dqsdp, dqsdt, dt, &
4501 & dt2, dtmax, dtmin, dv1h, &
4502 & dv1q, dv2h, dv2q, dv1u, &
4503 & dv1v, dv2u, dv2v, dv3q, &
4504 & dv3h, dv3u, dv3v, clam, &
4505 & dz, dz1, e1, &
4506 & el2orc, elocp, aafac, &
4507 & es, etah, h1, dthk, &
4508 & evef, evfact, evfactl, fact1, &
4509 & fact2, factor, fjcap, &
4510 & g, gamma, pprime, betaw, &
4511 & qlk, qrch, qs, c1, &
4512 & rain, rfact, shear, tem1, &
4513 & tem2, terr, val, val1, &
4514 & val2, w1, w1l, w1s, &
4515 & w2, w2l, w2s, w3, &
4516 & w3l, w3s, w4, w4l, &
4517 & w4s, tem, ptem, ptem1, &
4518 & pgcon
4519 !
4520 integer kb(im), kbcon(im), kbcon1(im), &
4521 & ktcon(im), ktcon1(im), &
4522 & kbm(im), kmax(im)
4523 !
4524 real(kind=kind_phys) aa1(im), &
4525 & delhbar(im), delq(im), delq2(im), &
4526 & delqbar(im), delqev(im), deltbar(im), &
4527 & deltv(im), edt(im), &
4528 & wstar(im), sflx(im), &
4529 & pdot(im), po(im,km), &
4530 & qcond(im), qevap(im), hmax(im), &
4531 & rntot(im), vshear(im), &
4532 & xlamud(im), xmb(im), xmbmax(im), &
4533 & delubar(im), delvbar(im)
4534 !c
4535 real(kind=kind_phys) cincr, cincrmax, cincrmin
4536 !cc
4537 !c physical parameters
4538 parameter(g=grav)
4539 parameter(cpoel=cp/hvap,elocp=hvap/cp, &
4540 & el2orc=hvap*hvap/(rv*cp))
4541 parameter(terr=0.,c0=.002,c1=5.e-4,delta=fv)
4542 parameter(fact1=(cvap-cliq)/rv,fact2=hvap/rv-fact1*t0c)
4543 parameter(cincrmax=180.,cincrmin=120.,dthk=25.)
4544 parameter(h1=0.33333333)
4545 !c local variables and arrays
4546 real(kind=kind_phys) pfld(im,km), to(im,km), qo(im,km), &
4547 & uo(im,km), vo(im,km), qeso(im,km)
4548 !c cloud water
4549 real(kind=kind_phys) qlko_ktcon(im), dellal(im,km), &
4550 & dbyo(im,km), zo(im,km), xlamue(im,km), &
4551 & heo(im,km), heso(im,km), &
4552 & dellah(im,km), dellaq(im,km), &
4553 & dellau(im,km), dellav(im,km), hcko(im,km), &
4554 & ucko(im,km), vcko(im,km), qcko(im,km), &
4555 & eta(im,km), zi(im,km), pwo(im,km), &
4556 & tx1(im)
4557 !
4558 logical totflg, cnvflg(im), flg(im)
4559 !
4560 real(kind=kind_phys) tf, tcr, tcrf
4561 parameter (tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf))
4562 !
4563 !c-----------------------------------------------------------------------
4564 !
4565 km1 = km - 1
4566 !c
4567 !c compute surface buoyancy flux
4568 !c
4569 do i=1,im
4570 sflx(i) = heat(i)+fv*t1(i,1)*evap(i)
4571 enddo
4572 !c
4573 !c initialize arrays
4574 !c
4575 do i=1,im
4576 cnvflg(i) = .true.
4577 if(kcnv(i).eq.1) cnvflg(i) = .false.
4578 if(sflx(i).le.0.) cnvflg(i) = .false.
4579 if(cnvflg(i)) then
4580 kbot(i)=km+1
4581 ktop(i)=0
4582 endif
4583 rn(i)=0.
4584 kbcon(i)=km
4585 ktcon(i)=1
4586 kb(i)=km
4587 pdot(i) = 0.
4588 qlko_ktcon(i) = 0.
4589 edt(i) = 0.
4590 aa1(i) = 0.
4591 vshear(i) = 0.
4592 enddo
4593 ! hchuang code change
4594 do k = 1, km
4595 do i = 1, im
4596 ud_mf(i,k) = 0.
4597 dt_mf(i,k) = 0.
4598 enddo
4599 enddo
4600 !!
4601 totflg = .true.
4602 do i=1,im
4603 totflg = totflg .and. (.not. cnvflg(i))
4604 enddo
4605 if(totflg) return
4606 !!
4607 !c
4608 dt2 = delt
4609 val = 1200.
4610 dtmin = max(dt2, val )
4611 val = 3600.
4612 dtmax = max(dt2, val )
4613 !c model tunable parameters are all here
4614 clam = .3
4615 aafac = .1
4616 betaw = .03
4617 !c evef = 0.07
4618 evfact = 0.3
4619 evfactl = 0.3
4620 !
4621 fjcap = (float(jcap) / 126.) ** 2
4622 val = 1.
4623 fjcap = max(fjcap,val)
4624 w1l = -8.e-3
4625 w2l = -4.e-2
4626 w3l = -5.e-3
4627 w4l = -5.e-4
4628 w1s = -2.e-4
4629 w2s = -2.e-3
4630 w3s = -1.e-3
4631 w4s = -2.e-5
4632 !c
4633 !c define top layer for search of the downdraft originating layer
4634 !c and the maximum thetae for updraft
4635 !c
4636 do i=1,im
4637 kbm(i) = km
4638 kmax(i) = km
4639 tx1(i) = 1.0 / ps(i)
4640 enddo
4641 !
4642 do k = 1, km
4643 do i=1,im
4644 if (prsl(i,k)*tx1(i) .gt. 0.70) kbm(i) = k + 1
4645 if (prsl(i,k)*tx1(i) .gt. 0.60) kmax(i) = k + 1
4646 enddo
4647 enddo
4648 do i=1,im
4649 kbm(i) = min(kbm(i),kmax(i))
4650 enddo
4651 !c
4652 !!c hydrostatic height assume zero terr and compute
4653 !c updraft entrainment rate as an inverse function of height
4654 !c
4655 do k = 1, km
4656 do i=1,im
4657 zo(i,k) = phil(i,k) / g
4658 enddo
4659 enddo
4660 do k = 1, km1
4661 do i=1,im
4662 zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1))
4663 xlamue(i,k) = clam / zi(i,k)
4664 enddo
4665 enddo
4666 do i=1,im
4667 xlamue(i,km) = xlamue(i,km1)
4668 enddo
4669 !c
4670 !c pbl height
4671 !c
4672 do i=1,im
4673 flg(i) = cnvflg(i)
4674 kpbl(i)= 1
4675 enddo
4676 do k = 2, km1
4677 do i=1,im
4678 if (flg(i).and.zo(i,k).le.hpbl(i)) then
4679 kpbl(i) = k
4680 else
4681 flg(i) = .false.
4682 endif
4683 enddo
4684 enddo
4685 do i=1,im
4686 kpbl(i)= min(kpbl(i),kbm(i))
4687 enddo
4688 !c
4689 !c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4690 !c convert surface pressure to mb from cb
4691 !c
4692 do k = 1, km
4693 do i = 1, im
4694 if (cnvflg(i) .and. k .le. kmax(i)) then
4695 pfld(i,k) = prsl(i,k) * 10.0
4696 eta(i,k) = 1.
4697 hcko(i,k) = 0.
4698 qcko(i,k) = 0.
4699 ucko(i,k) = 0.
4700 vcko(i,k) = 0.
4701 dbyo(i,k) = 0.
4702 pwo(i,k) = 0.
4703 dellal(i,k) = 0.
4704 to(i,k) = t1(i,k)
4705 qo(i,k) = q1(i,k)
4706 uo(i,k) = u1(i,k) * rcs(i)
4707 vo(i,k) = v1(i,k) * rcs(i)
4708 endif
4709 enddo
4710 enddo
4711 !c
4712 !c column variables
4713 !c p is pressure of the layer (mb)
4714 !c t is temperature at t-dt (k)..tn
4715 !c q is mixing ratio at t-dt (kg/kg)..qn
4716 !c to is temperature at t+dt (k)... this is after advection and turbulan
4717 !c qo is mixing ratio at t+dt (kg/kg)..q1
4718 !c
4719 do k = 1, km
4720 do i=1,im
4721 if (cnvflg(i) .and. k .le. kmax(i)) then
4722 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
4723 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
4724 val1 = 1.e-8
4725 qeso(i,k) = max(qeso(i,k), val1)
4726 val2 = 1.e-10
4727 qo(i,k) = max(qo(i,k), val2 )
4728 ! qo(i,k) = min(qo(i,k),qeso(i,k))
4729 ! tvo(i,k) = to(i,k) + delta * to(i,k) * qo(i,k)
4730 endif
4731 enddo
4732 enddo
4733 !c
4734 !c compute moist static energy
4735 !c
4736 do k = 1, km
4737 do i=1,im
4738 if (cnvflg(i) .and. k .le. kmax(i)) then
4739 ! tem = g * zo(i,k) + cp * to(i,k)
4740 tem = phil(i,k) + cp * to(i,k)
4741 heo(i,k) = tem + hvap * qo(i,k)
4742 heso(i,k) = tem + hvap * qeso(i,k)
4743 !c heo(i,k) = min(heo(i,k),heso(i,k))
4744 endif
4745 enddo
4746 enddo
4747 !c
4748 !c determine level with largest moist static energy within pbl
4749 !c this is the level where updraft starts
4750 !c
4751 do i=1,im
4752 if (cnvflg(i)) then
4753 hmax(i) = heo(i,1)
4754 kb(i) = 1
4755 endif
4756 enddo
4757 do k = 2, km
4758 do i=1,im
4759 if (cnvflg(i).and.k.le.kpbl(i)) then
4760 if(heo(i,k).gt.hmax(i)) then
4761 kb(i) = k
4762 hmax(i) = heo(i,k)
4763 endif
4764 endif
4765 enddo
4766 enddo
4767 !c
4768 do k = 1, km1
4769 do i=1,im
4770 if (cnvflg(i) .and. k .le. kmax(i)-1) then
4771 dz = .5 * (zo(i,k+1) - zo(i,k))
4772 dp = .5 * (pfld(i,k+1) - pfld(i,k))
4773 es = 0.01 * fpvs(to(i,k+1)) ! fpvs is in pa
4774 pprime = pfld(i,k+1) + epsm1 * es
4775 qs = eps * es / pprime
4776 dqsdp = - qs / pprime
4777 desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
4778 dqsdt = qs * pfld(i,k+1) * desdt / (es * pprime)
4779 gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
4780 dt = (g * dz + hvap * dqsdp * dp) / (cp * (1. + gamma))
4781 dq = dqsdt * dt + dqsdp * dp
4782 to(i,k) = to(i,k+1) + dt
4783 qo(i,k) = qo(i,k+1) + dq
4784 po(i,k) = .5 * (pfld(i,k) + pfld(i,k+1))
4785 endif
4786 enddo
4787 enddo
4788 !
4789 do k = 1, km1
4790 do i=1,im
4791 if (cnvflg(i) .and. k .le. kmax(i)-1) then
4792 qeso(i,k) = 0.01 * fpvs(to(i,k)) ! fpvs is in pa
4793 qeso(i,k) = eps * qeso(i,k) / (po(i,k) + epsm1*qeso(i,k))
4794 val1 = 1.e-8
4795 qeso(i,k) = max(qeso(i,k), val1)
4796 val2 = 1.e-10
4797 qo(i,k) = max(qo(i,k), val2 )
4798 ! qo(i,k) = min(qo(i,k),qeso(i,k))
4799 heo(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) + &
4800 & cp * to(i,k) + hvap * qo(i,k)
4801 heso(i,k) = .5 * g * (zo(i,k) + zo(i,k+1)) + &
4802 & cp * to(i,k) + hvap * qeso(i,k)
4803 uo(i,k) = .5 * (uo(i,k) + uo(i,k+1))
4804 vo(i,k) = .5 * (vo(i,k) + vo(i,k+1))
4805 endif
4806 enddo
4807 enddo
4808 !c
4809 !c look for the level of free convection as cloud base
4810 !c
4811 do i=1,im
4812 flg(i) = cnvflg(i)
4813 if(flg(i)) kbcon(i) = kmax(i)
4814 enddo
4815 do k = 2, km1
4816 do i=1,im
4817 if (flg(i).and.k.lt.kbm(i)) then
4818 if(k.gt.kb(i).and.heo(i,kb(i)).gt.heso(i,k)) then
4819 kbcon(i) = k
4820 flg(i) = .false.
4821 endif
4822 endif
4823 enddo
4824 enddo
4825 !c
4826 do i=1,im
4827 if(cnvflg(i)) then
4828 if(kbcon(i).eq.kmax(i)) cnvflg(i) = .false.
4829 endif
4830 enddo
4831 !!
4832 totflg = .true.
4833 do i=1,im
4834 totflg = totflg .and. (.not. cnvflg(i))
4835 enddo
4836 if(totflg) return
4837 !!
4838 !c
4839 !c determine critical convective inhibition
4840 !c as a function of vertical velocity at cloud base.
4841 !c
4842 do i=1,im
4843 if(cnvflg(i)) then
4844 pdot(i) = 10.* dot(i,kbcon(i))
4845 endif
4846 enddo
4847 do i=1,im
4848 if(cnvflg(i)) then
4849 if(slimsk(i).eq.1.) then
4850 w1 = w1l
4851 w2 = w2l
4852 w3 = w3l
4853 w4 = w4l
4854 else
4855 w1 = w1s
4856 w2 = w2s
4857 w3 = w3s
4858 w4 = w4s
4859 endif
4860 if(pdot(i).le.w4) then
4861 ptem = (pdot(i) - w4) / (w3 - w4)
4862 elseif(pdot(i).ge.-w4) then
4863 ptem = - (pdot(i) + w4) / (w4 - w3)
4864 else
4865 ptem = 0.
4866 endif
4867 val1 = -1.
4868 ptem = max(ptem,val1)
4869 val2 = 1.
4870 ptem = min(ptem,val2)
4871 ptem = 1. - ptem
4872 ptem1= .5*(cincrmax-cincrmin)
4873 cincr = cincrmax - ptem * ptem1
4874 tem1 = pfld(i,kb(i)) - pfld(i,kbcon(i))
4875 if(tem1.gt.cincr) then
4876 cnvflg(i) = .false.
4877 endif
4878 endif
4879 enddo
4880 !!
4881 totflg = .true.
4882 do i=1,im
4883 totflg = totflg .and. (.not. cnvflg(i))
4884 enddo
4885 if(totflg) return
4886 !!
4887 !c
4888 !c assume the detrainment rate for the updrafts to be same as
4889 !c the entrainment rate at cloud base
4890 !c
4891 do i = 1, im
4892 if(cnvflg(i)) then
4893 xlamud(i) = xlamue(i,kbcon(i))
4894 endif
4895 enddo
4896 !c
4897 !c determine updraft mass flux for the subcloud layers
4898 !c
4899 do k = km1, 1, -1
4900 do i = 1, im
4901 if (cnvflg(i)) then
4902 if(k.lt.kbcon(i).and.k.ge.kb(i)) then
4903 dz = zi(i,k+1) - zi(i,k)
4904 ptem = 0.5*(xlamue(i,k)+xlamue(i,k+1))-xlamud(i)
4905 eta(i,k) = eta(i,k+1) / (1. + ptem * dz)
4906 endif
4907 endif
4908 enddo
4909 enddo
4910 !c
4911 !c compute mass flux above cloud base
4912 !c
4913 do k = 2, km1
4914 do i = 1, im
4915 if(cnvflg(i))then
4916 if(k.gt.kbcon(i).and.k.lt.kmax(i)) then
4917 dz = zi(i,k) - zi(i,k-1)
4918 ptem = 0.5*(xlamue(i,k)+xlamue(i,k-1))-xlamud(i)
4919 eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
4920 endif
4921 endif
4922 enddo
4923 enddo
4924 !c
4925 !c compute updraft cloud property
4926 !c
4927 do i = 1, im
4928 if(cnvflg(i)) then
4929 indx = kb(i)
4930 hcko(i,indx) = heo(i,indx)
4931 ucko(i,indx) = uo(i,indx)
4932 vcko(i,indx) = vo(i,indx)
4933 endif
4934 enddo
4935 !c
4936 do k = 2, km1
4937 do i = 1, im
4938 if (cnvflg(i)) then
4939 if(k.gt.kb(i).and.k.lt.kmax(i)) then
4940 dz = zi(i,k) - zi(i,k-1)
4941 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
4942 tem1 = 0.5 * xlamud(i) * dz
4943 factor = 1. + tem - tem1
4944 ptem = 0.5 * tem + pgcon
4945 ptem1= 0.5 * tem - pgcon
4946 hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5* &
4947 & (heo(i,k)+heo(i,k-1)))/factor
4948 ucko(i,k) = ((1.-tem1)*ucko(i,k-1)+ptem*uo(i,k) &
4949 & +ptem1*uo(i,k-1))/factor
4950 vcko(i,k) = ((1.-tem1)*vcko(i,k-1)+ptem*vo(i,k) &
4951 & +ptem1*vo(i,k-1))/factor
4952 dbyo(i,k) = hcko(i,k) - heso(i,k)
4953 endif
4954 endif
4955 enddo
4956 enddo
4957 !c
4958 !c taking account into convection inhibition due to existence of
4959 !c dry layers below cloud base
4960 !c
4961 do i=1,im
4962 flg(i) = cnvflg(i)
4963 kbcon1(i) = kmax(i)
4964 enddo
4965 do k = 2, km1
4966 do i=1,im
4967 if (flg(i).and.k.lt.kbm(i)) then
4968 if(k.ge.kbcon(i).and.dbyo(i,k).gt.0.) then
4969 kbcon1(i) = k
4970 flg(i) = .false.
4971 endif
4972 endif
4973 enddo
4974 enddo
4975 do i=1,im
4976 if(cnvflg(i)) then
4977 if(kbcon1(i).eq.kmax(i)) cnvflg(i) = .false.
4978 endif
4979 enddo
4980 do i=1,im
4981 if(cnvflg(i)) then
4982 tem = pfld(i,kbcon(i)) - pfld(i,kbcon1(i))
4983 if(tem.gt.dthk) then
4984 cnvflg(i) = .false.
4985 endif
4986 endif
4987 enddo
4988 !!
4989 totflg = .true.
4990 do i = 1, im
4991 totflg = totflg .and. (.not. cnvflg(i))
4992 enddo
4993 if(totflg) return
4994 !!
4995 !c
4996 !c determine first guess cloud top as the level of zero buoyancy
4997 !c limited to the level of sigma=0.7
4998 !c
4999 do i = 1, im
5000 flg(i) = cnvflg(i)
5001 if(flg(i)) ktcon(i) = kbm(i)
5002 enddo
5003 do k = 2, km1
5004 do i=1,im
5005 if (flg(i).and.k .lt. kbm(i)) then
5006 if(k.gt.kbcon1(i).and.dbyo(i,k).lt.0.) then
5007 ktcon(i) = k
5008 flg(i) = .false.
5009 endif
5010 endif
5011 enddo
5012 enddo
5013 !c
5014 !c turn off shallow convection if cloud top is less than pbl top
5015 !c
5016 ! do i=1,im
5017 ! if(cnvflg(i)) then
5018 ! kk = kpbl(i)+1
5019 ! if(ktcon(i).le.kk) cnvflg(i) = .false.
5020 ! endif
5021 ! enddo
5022 !!
5023 ! totflg = .true.
5024 ! do i = 1, im
5025 ! totflg = totflg .and. (.not. cnvflg(i))
5026 ! enddo
5027 ! if(totflg) return
5028 !!
5029 !c
5030 !c specify upper limit of mass flux at cloud base
5031 !c
5032 do i = 1, im
5033 if(cnvflg(i)) then
5034 ! xmbmax(i) = .1
5035 !
5036 k = kbcon(i)
5037 dp = 1000. * del(i,k)
5038 xmbmax(i) = dp / (g * dt2)
5039 !
5040 ! tem = dp / (g * dt2)
5041 ! xmbmax(i) = min(tem, xmbmax(i))
5042 endif
5043 enddo
5044 !c
5045 !c compute cloud moisture property and precipitation
5046 !c
5047 do i = 1, im
5048 if (cnvflg(i)) then
5049 aa1(i) = 0.
5050 qcko(i,kb(i)) = qo(i,kb(i))
5051 endif
5052 enddo
5053 do k = 2, km1
5054 do i = 1, im
5055 if (cnvflg(i)) then
5056 if(k.gt.kb(i).and.k.lt.ktcon(i)) then
5057 dz = zi(i,k) - zi(i,k-1)
5058 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
5059 qrch = qeso(i,k) &
5060 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
5061 !cj
5062 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
5063 tem1 = 0.5 * xlamud(i) * dz
5064 factor = 1. + tem - tem1
5065 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* &
5066 & (qo(i,k)+qo(i,k-1)))/factor
5067 !cj
5068 dq = eta(i,k) * (qcko(i,k) - qrch)
5069 !c
5070 ! rhbar(i) = rhbar(i) + qo(i,k) / qeso(i,k)
5071 !c
5072 !c below lfc check if there is excess moisture to release latent heat
5073 !c
5074 if(k.ge.kbcon(i).and.dq.gt.0.) then
5075 etah = .5 * (eta(i,k) + eta(i,k-1))
5076 if(ncloud.gt.0.) then
5077 dp = 1000. * del(i,k)
5078 qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
5079 dellal(i,k) = etah * c1 * dz * qlk * g / dp
5080 else
5081 qlk = dq / (eta(i,k) + etah * c0 * dz)
5082 endif
5083 aa1(i) = aa1(i) - dz * g * qlk
5084 qcko(i,k)= qlk + qrch
5085 pwo(i,k) = etah * c0 * dz * qlk
5086 endif
5087 endif
5088 endif
5089 enddo
5090 enddo
5091 !c
5092 !c calculate cloud work function
5093 !c
5094 do k = 2, km1
5095 do i = 1, im
5096 if (cnvflg(i)) then
5097 if(k.ge.kbcon(i).and.k.lt.ktcon(i)) then
5098 dz1 = zo(i,k+1) - zo(i,k)
5099 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
5100 rfact = 1. + delta * cp * gamma &
5101 & * to(i,k) / hvap
5102 aa1(i) = aa1(i) + &
5103 & dz1 * (g / (cp * to(i,k))) &
5104 & * dbyo(i,k) / (1. + gamma) &
5105 & * rfact
5106 val = 0.
5107 aa1(i)=aa1(i)+ &
5108 & dz1 * g * delta * &
5109 & max(val,(qeso(i,k) - qo(i,k)))
5110 endif
5111 endif
5112 enddo
5113 enddo
5114 do i = 1, im
5115 if(cnvflg(i).and.aa1(i).le.0.) cnvflg(i) = .false.
5116 enddo
5117 !!
5118 totflg = .true.
5119 do i=1,im
5120 totflg = totflg .and. (.not. cnvflg(i))
5121 enddo
5122 if(totflg) return
5123 !!
5124 !c
5125 !c estimate the onvective overshooting as the level
5126 !c where the [aafac * cloud work function] becomes zero,
5127 !c which is the final cloud top
5128 !c limited to the level of sigma=0.7
5129 !c
5130 do i = 1, im
5131 if (cnvflg(i)) then
5132 aa1(i) = aafac * aa1(i)
5133 endif
5134 enddo
5135 !c
5136 do i = 1, im
5137 flg(i) = cnvflg(i)
5138 ktcon1(i) = kbm(i)
5139 enddo
5140 do k = 2, km1
5141 do i = 1, im
5142 if (flg(i)) then
5143 if(k.ge.ktcon(i).and.k.lt.kbm(i)) then
5144 dz1 = zo(i,k+1) - zo(i,k)
5145 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
5146 rfact = 1. + delta * cp * gamma &
5147 & * to(i,k) / hvap
5148 aa1(i) = aa1(i) + &
5149 & dz1 * (g / (cp * to(i,k))) &
5150 & * dbyo(i,k) / (1. + gamma) &
5151 & * rfact
5152 if(aa1(i).lt.0.) then
5153 ktcon1(i) = k
5154 flg(i) = .false.
5155 endif
5156 endif
5157 endif
5158 enddo
5159 enddo
5160 !c
5161 !c compute cloud moisture property, detraining cloud water
5162 !c and precipitation in overshooting layers
5163 !c
5164 do k = 2, km1
5165 do i = 1, im
5166 if (cnvflg(i)) then
5167 if(k.ge.ktcon(i).and.k.lt.ktcon1(i)) then
5168 dz = zi(i,k) - zi(i,k-1)
5169 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
5170 qrch = qeso(i,k) &
5171 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
5172 !cj
5173 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
5174 tem1 = 0.5 * xlamud(i) * dz
5175 factor = 1. + tem - tem1
5176 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* &
5177 & (qo(i,k)+qo(i,k-1)))/factor
5178 !cj
5179 dq = eta(i,k) * (qcko(i,k) - qrch)
5180 !c
5181 !c check if there is excess moisture to release latent heat
5182 !c
5183 if(dq.gt.0.) then
5184 etah = .5 * (eta(i,k) + eta(i,k-1))
5185 if(ncloud.gt.0.) then
5186 dp = 1000. * del(i,k)
5187 qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
5188 dellal(i,k) = etah * c1 * dz * qlk * g / dp
5189 else
5190 qlk = dq / (eta(i,k) + etah * c0 * dz)
5191 endif
5192 qcko(i,k) = qlk + qrch
5193 pwo(i,k) = etah * c0 * dz * qlk
5194 endif
5195 endif
5196 endif
5197 enddo
5198 enddo
5199 !c
5200 !c exchange ktcon with ktcon1
5201 !c
5202 do i = 1, im
5203 if(cnvflg(i)) then
5204 kk = ktcon(i)
5205 ktcon(i) = ktcon1(i)
5206 ktcon1(i) = kk
5207 endif
5208 enddo
5209 !c
5210 !c this section is ready for cloud water
5211 !c
5212 if(ncloud.gt.0) then
5213 !c
5214 !c compute liquid and vapor separation at cloud top
5215 !c
5216 do i = 1, im
5217 if(cnvflg(i)) then
5218 k = ktcon(i) - 1
5219 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
5220 qrch = qeso(i,k) &
5221 & + gamma * dbyo(i,k) / (hvap * (1. + gamma))
5222 dq = qcko(i,k) - qrch
5223 !c
5224 !c check if there is excess moisture to release latent heat
5225 !c
5226 if(dq.gt.0.) then
5227 qlko_ktcon(i) = dq
5228 qcko(i,k) = qrch
5229 endif
5230 endif
5231 enddo
5232 endif
5233 !!c
5234 !c--- compute precipitation efficiency in terms of windshear
5235 !c
5236 do i = 1, im
5237 if(cnvflg(i)) then
5238 vshear(i) = 0.
5239 endif
5240 enddo
5241 do k = 2, km
5242 do i = 1, im
5243 if (cnvflg(i)) then
5244 if(k.gt.kb(i).and.k.le.ktcon(i)) then
5245 shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2 &
5246 & + (vo(i,k)-vo(i,k-1)) ** 2)
5247 vshear(i) = vshear(i) + shear
5248 endif
5249 endif
5250 enddo
5251 enddo
5252 do i = 1, im
5253 if(cnvflg(i)) then
5254 vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i))-zi(i,kb(i)))
5255 e1=1.591-.639*vshear(i) &
5256 & +.0953*(vshear(i)**2)-.00496*(vshear(i)**3)
5257 edt(i)=1.-e1
5258 val = .9
5259 edt(i) = min(edt(i),val)
5260 val = .0
5261 edt(i) = max(edt(i),val)
5262 endif
5263 enddo
5264 !c
5265 !c--- what would the change be, that a cloud with unit mass
5266 !c--- will do to the environment?
5267 !c
5268 do k = 1, km
5269 do i = 1, im
5270 if(cnvflg(i) .and. k .le. kmax(i)) then
5271 dellah(i,k) = 0.
5272 dellaq(i,k) = 0.
5273 dellau(i,k) = 0.
5274 dellav(i,k) = 0.
5275 endif
5276 enddo
5277 enddo
5278 !c
5279 !c--- changed due to subsidence and entrainment
5280 !c
5281 do k = 2, km1
5282 do i = 1, im
5283 if (cnvflg(i)) then
5284 if(k.gt.kb(i).and.k.lt.ktcon(i)) then
5285 dp = 1000. * del(i,k)
5286 dz = zi(i,k) - zi(i,k-1)
5287 !c
5288 dv1h = heo(i,k)
5289 dv2h = .5 * (heo(i,k) + heo(i,k-1))
5290 dv3h = heo(i,k-1)
5291 dv1q = qo(i,k)
5292 dv2q = .5 * (qo(i,k) + qo(i,k-1))
5293 dv3q = qo(i,k-1)
5294 dv1u = uo(i,k)
5295 dv2u = .5 * (uo(i,k) + uo(i,k-1))
5296 dv3u = uo(i,k-1)
5297 dv1v = vo(i,k)
5298 dv2v = .5 * (vo(i,k) + vo(i,k-1))
5299 dv3v = vo(i,k-1)
5300 !c
5301 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1))
5302 tem1 = xlamud(i)
5303 !cj
5304 dellah(i,k) = dellah(i,k) + &
5305 & ( eta(i,k)*dv1h - eta(i,k-1)*dv3h &
5306 & - tem*eta(i,k-1)*dv2h*dz &
5307 & + tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz &
5308 & ) *g/dp
5309 !cj
5310 dellaq(i,k) = dellaq(i,k) + &
5311 & ( eta(i,k)*dv1q - eta(i,k-1)*dv3q &
5312 & - tem*eta(i,k-1)*dv2q*dz &
5313 & + tem1*eta(i,k-1)*.5*(qcko(i,k)+qcko(i,k-1))*dz &
5314 & ) *g/dp
5315 !cj
5316 dellau(i,k) = dellau(i,k) + &
5317 & ( eta(i,k)*dv1u - eta(i,k-1)*dv3u &
5318 & - tem*eta(i,k-1)*dv2u*dz &
5319 & + tem1*eta(i,k-1)*.5*(ucko(i,k)+ucko(i,k-1))*dz &
5320 & - pgcon*eta(i,k-1)*(dv1u-dv3u) &
5321 & ) *g/dp
5322 !cj
5323 dellav(i,k) = dellav(i,k) + &
5324 & ( eta(i,k)*dv1v - eta(i,k-1)*dv3v &
5325 & - tem*eta(i,k-1)*dv2v*dz &
5326 & + tem1*eta(i,k-1)*.5*(vcko(i,k)+vcko(i,k-1))*dz &
5327 & - pgcon*eta(i,k-1)*(dv1v-dv3v) &
5328 & ) *g/dp
5329 !cj
5330 endif
5331 endif
5332 enddo
5333 enddo
5334 !c
5335 !c------- cloud top
5336 !c
5337 do i = 1, im
5338 if(cnvflg(i)) then
5339 indx = ktcon(i)
5340 dp = 1000. * del(i,indx)
5341 dv1h = heo(i,indx-1)
5342 dellah(i,indx) = eta(i,indx-1) * &
5343 & (hcko(i,indx-1) - dv1h) * g / dp
5344 dv1q = qo(i,indx-1)
5345 dellaq(i,indx) = eta(i,indx-1) * &
5346 & (qcko(i,indx-1) - dv1q) * g / dp
5347 dv1u = uo(i,indx-1)
5348 dellau(i,indx) = eta(i,indx-1) * &
5349 & (ucko(i,indx-1) - dv1u) * g / dp
5350 dv1v = vo(i,indx-1)
5351 dellav(i,indx) = eta(i,indx-1) * &
5352 & (vcko(i,indx-1) - dv1v) * g / dp
5353 !c
5354 !c cloud water
5355 !c
5356 dellal(i,indx) = eta(i,indx-1) * &
5357 & qlko_ktcon(i) * g / dp
5358 endif
5359 enddo
5360 !c
5361 !c mass flux at cloud base for shallow convection
5362 !c (Grant, 2001)
5363 !c
5364 do i= 1, im
5365 if(cnvflg(i)) then
5366 k = kbcon(i)
5367 ! ptem = g*sflx(i)*zi(i,k)/t1(i,1)
5368 ptem = g*sflx(i)*hpbl(i)/t1(i,1)
5369 wstar(i) = ptem**h1
5370 tem = po(i,k)*100. / (rd*t1(i,k))
5371 xmb(i) = betaw*tem*wstar(i)
5372 xmb(i) = min(xmb(i),xmbmax(i))
5373 endif
5374 enddo
5375 !c
5376 !c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
5377 !c
5378 do k = 1, km
5379 do i = 1, im
5380 if (cnvflg(i) .and. k .le. kmax(i)) then
5381 qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
5382 qeso(i,k) = eps * qeso(i,k) / (pfld(i,k) + epsm1*qeso(i,k))
5383 val = 1.e-8
5384 qeso(i,k) = max(qeso(i,k), val )
5385 endif
5386 enddo
5387 enddo
5388 !c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
5389 !c
5390 do i = 1, im
5391 delhbar(i) = 0.
5392 delqbar(i) = 0.
5393 deltbar(i) = 0.
5394 delubar(i) = 0.
5395 delvbar(i) = 0.
5396 qcond(i) = 0.
5397 enddo
5398 do k = 1, km
5399 do i = 1, im
5400 if (cnvflg(i)) then
5401 if(k.gt.kb(i).and.k.le.ktcon(i)) then
5402 dellat = (dellah(i,k) - hvap * dellaq(i,k)) / cp
5403 t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2
5404 q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2
5405 tem = 1./rcs(i)
5406 u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 * tem
5407 v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 * tem
5408 dp = 1000. * del(i,k)
5409 delhbar(i) = delhbar(i) + dellah(i,k)*xmb(i)*dp/g
5410 delqbar(i) = delqbar(i) + dellaq(i,k)*xmb(i)*dp/g
5411 deltbar(i) = deltbar(i) + dellat*xmb(i)*dp/g
5412 delubar(i) = delubar(i) + dellau(i,k)*xmb(i)*dp/g
5413 delvbar(i) = delvbar(i) + dellav(i,k)*xmb(i)*dp/g
5414 endif
5415 endif
5416 enddo
5417 enddo
5418 do k = 1, km
5419 do i = 1, im
5420 if (cnvflg(i)) then
5421 if(k.gt.kb(i).and.k.le.ktcon(i)) then
5422 qeso(i,k) = 0.01 * fpvs(t1(i,k)) ! fpvs is in pa
5423 qeso(i,k) = eps * qeso(i,k)/(pfld(i,k) + epsm1*qeso(i,k))
5424 val = 1.e-8
5425 qeso(i,k) = max(qeso(i,k), val )
5426 endif
5427 endif
5428 enddo
5429 enddo
5430 !c
5431 do i = 1, im
5432 rntot(i) = 0.
5433 delqev(i) = 0.
5434 delq2(i) = 0.
5435 flg(i) = cnvflg(i)
5436 enddo
5437 do k = km, 1, -1
5438 do i = 1, im
5439 if (cnvflg(i)) then
5440 if(k.lt.ktcon(i).and.k.gt.kb(i)) then
5441 rntot(i) = rntot(i) + pwo(i,k) * xmb(i) * .001 * dt2
5442 endif
5443 endif
5444 enddo
5445 enddo
5446 !c
5447 !c evaporating rain
5448 !c
5449 do k = km, 1, -1
5450 do i = 1, im
5451 if (k .le. kmax(i)) then
5452 deltv(i) = 0.
5453 delq(i) = 0.
5454 qevap(i) = 0.
5455 if(cnvflg(i)) then
5456 if(k.lt.ktcon(i).and.k.gt.kb(i)) then
5457 rn(i) = rn(i) + pwo(i,k) * xmb(i) * .001 * dt2
5458 endif
5459 endif
5460 if(flg(i).and.k.lt.ktcon(i)) then
5461 evef = edt(i) * evfact
5462 if(slimsk(i).eq.1.) evef=edt(i) * evfactl
5463 ! if(slimsk(i).eq.1.) evef=.07
5464 !c if(slimsk(i).ne.1.) evef = 0.
5465 qcond(i) = evef * (q1(i,k) - qeso(i,k)) &
5466 & / (1. + el2orc * qeso(i,k) / t1(i,k)**2)
5467 dp = 1000. * del(i,k)
5468 if(rn(i).gt.0..and.qcond(i).lt.0.) then
5469 qevap(i) = -qcond(i) * (1.-exp(-.32*sqrt(dt2*rn(i))))
5470 qevap(i) = min(qevap(i), rn(i)*1000.*g/dp)
5471 delq2(i) = delqev(i) + .001 * qevap(i) * dp / g
5472 endif
5473 if(rn(i).gt.0..and.qcond(i).lt.0..and. &
5474 & delq2(i).gt.rntot(i)) then
5475 qevap(i) = 1000.* g * (rntot(i) - delqev(i)) / dp
5476 flg(i) = .false.
5477 endif
5478 if(rn(i).gt.0..and.qevap(i).gt.0.) then
5479 tem = .001 * dp / g
5480 tem1 = qevap(i) * tem
5481 if(tem1.gt.rn(i)) then
5482 qevap(i) = rn(i) / tem
5483 rn(i) = 0.
5484 else
5485 rn(i) = rn(i) - tem1
5486 endif
5487 q1(i,k) = q1(i,k) + qevap(i)
5488 t1(i,k) = t1(i,k) - elocp * qevap(i)
5489 deltv(i) = - elocp*qevap(i)/dt2
5490 delq(i) = + qevap(i)/dt2
5491 delqev(i) = delqev(i) + .001*dp*qevap(i)/g
5492 endif
5493 dellaq(i,k) = dellaq(i,k) + delq(i) / xmb(i)
5494 delqbar(i) = delqbar(i) + delq(i)*dp/g
5495 deltbar(i) = deltbar(i) + deltv(i)*dp/g
5496 endif
5497 endif
5498 enddo
5499 enddo
5500 !cj
5501 ! do i = 1, im
5502 ! if(me.eq.31.and.cnvflg(i)) then
5503 ! if(cnvflg(i)) then
5504 ! print *, ' shallow delhbar, delqbar, deltbar = ',
5505 ! & delhbar(i),hvap*delqbar(i),cp*deltbar(i)
5506 ! print *, ' shallow delubar, delvbar = ',delubar(i),delvbar(i)
5507 ! print *, ' precip =', hvap*rn(i)*1000./dt2
5508 ! print*,'pdif= ',pfld(i,kbcon(i))-pfld(i,ktcon(i))
5509 ! endif
5510 ! enddo
5511 !cj
5512 do i = 1, im
5513 if(cnvflg(i)) then
5514 if(rn(i).lt.0..or..not.flg(i)) rn(i) = 0.
5515 ktop(i) = ktcon(i)
5516 kbot(i) = kbcon(i)
5517 kcnv(i) = 0
5518 endif
5519 enddo
5520 !c
5521 !c cloud water
5522 !c
5523 if (ncloud.gt.0) then
5524 !
5525 val1 = 1.0
5526 val2 = 0.
5527 do k = 1, km1
5528 do i = 1, im
5529 if (cnvflg(i)) then
5530 if (k.gt.kb(i).and.k.le.ktcon(i)) then
5531 tem = dellal(i,k) * xmb(i) * dt2
5532 ! tem1 = max(0.0, min(1.0, (tcr-t1(i,k))*tcrf))
5533 tem1 = max(val2, min(val1, (tcr-t1(i,k))*tcrf))
5534 if (ql(i,k,2) .gt. -999.0) then
5535 ql(i,k,1) = ql(i,k,1) + tem * tem1 ! ice
5536 ql(i,k,2) = ql(i,k,2) + tem *(1.0-tem1) ! water
5537 else
5538 ql(i,k,1) = ql(i,k,1) + tem
5539 endif
5540 endif
5541 endif
5542 enddo
5543 enddo
5544 !
5545 endif
5546 !
5547 ! hchuang code change
5548 !
5549 do k = 1, km
5550 do i = 1, im
5551 if(cnvflg(i)) then
5552 if(k.ge.kb(i) .and. k.lt.ktop(i)) then
5553 ud_mf(i,k) = eta(i,k) * xmb(i) * dt2
5554 endif
5555 endif
5556 enddo
5557 enddo
5558 do i = 1, im
5559 if(cnvflg(i)) then
5560 k = ktop(i)-1
5561 dt_mf(i,k) = ud_mf(i,k)
5562 endif
5563 enddo
5564 !!
5565 return
5566 end subroutine shalcnv
5567 END MODULE module_cu_sas