Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_cu_sas.F @ 13

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

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

WRF: version v3.3
LMDZ: version v1818

More details in:

Property svn:executable set to *

File size: 85.9 KB

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

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

Download in other formats:

Original Format