Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_cu_sas.F @ 3094

Last change on this file since 3094 was 2759, checked in by aslmd, 2 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
Property svn:executable set to ``*
File size: 84.5 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format