Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_cu_sas.F @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

Property svn:executable set to *

File size: 85.9 KB

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

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

Download in other formats:

Original Format