Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_cu_sas.F @ 3026

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

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

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

Download in other formats:

Original Format