Context Navigation

module_cu_sas.F

Last change on this file was 11, checked in by aslmd, 14 years ago
spiga@svn-planeto:ajoute le modele meso-echelle martien
Property svn:executable set to ``*
File size: 83.9 KB

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

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

Download in other formats:

Original Format