Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_cu_bmj.F @ 3094

Last change on this file since 3094 was 2759, checked in by aslmd, 2 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
File size: 67.4 KB

Rev	Line
[2759]	1	!-----------------------------------------------------------------------
	2	!
	3	!WRF:MODEL_LAYER:PHYSICS
	4	!
	5	!-----------------------------------------------------------------------
	6	!
	7	MODULE MODULE_CU_BMJ
	8	!
	9	!-----------------------------------------------------------------------
	10	USE MODULE_MODEL_CONSTANTS
	11	!-----------------------------------------------------------------------
	12	!
	13	REAL,PARAMETER :: &
	14	& DSPC=-3000. &
	15	& ,DTTOP=0.,EFIFC=5.0,EFIMN=0.20,EFMNT=0.70 &
	16	& ,ELIWV=2.683E6,ENPLO=20000.,ENPUP=15000. &
	17	& ,EPSDN=1.05,EPSDT=0. &
	18	& ,EPSNTP=.0001,EPSNTT=.0001,EPSPR=1.E-7 &
	19	& ,EPSUP=1.00 &
	20	& ,FR=1.00,FSL=0.85,FSS=0.85 &
	21	& ,FUP=0. &
	22	& ,PBM=13000.,PFRZ=15000.,PNO=1000. &
	23	& ,PONE=2500.,PQM=20000. &
	24	& ,PSH=20000.,PSHU=45000. &
	25	& ,RENDP=1./(ENPLO-ENPUP) &
	26	& ,RHLSC=0.00,RHHSC=1.10 &
	27	& ,ROW=1.E3 &
	28	& ,STABDF=0.90,STABDS=0.90 &
	29	& ,STABS=1.0,STRESH=1.10 &
	30	& ,DTSHAL=-1.0,TREL=2400.
	31	!
	32	REAL,PARAMETER :: DTtrigr=-0.0 &
	33	,DTPtrigr=DTtrigr*PONE !<-- Average parcel virtual temperature deficit over depth PONE.
	34	!<-- NOTE: CAPEtrigr is scaled by the cloud base temperature (see below)
	35	!
	36	REAL,PARAMETER :: DSPBFL=-3875.*FR &
	37	& ,DSP0FL=-5875.*FR &
	38	& ,DSPTFL=-1875.*FR &
	39	& ,DSPBFS=-3875. &
	40	& ,DSP0FS=-5875. &
	41	& ,DSPTFS=-1875.
	42	!
	43	REAL,PARAMETER :: PL=2500.,PLQ=70000.,PH=105000. &
	44	& ,THL=210.,THH=365.,THHQ=325.
	45	!
	46	INTEGER,PARAMETER :: ITB=76,JTB=134,ITBQ=152,JTBQ=440
	47	!
	48	INTEGER,PARAMETER :: ITREFI_MAX=3
	49	!
	50	!*** ARRAYS FOR LOOKUP TABLES
	51	!
	52	REAL,DIMENSION(ITB),PRIVATE,SAVE :: STHE,THE0
	53	REAL,DIMENSION(JTB),PRIVATE,SAVE :: QS0,SQS
	54	REAL,DIMENSION(ITBQ),PRIVATE,SAVE :: STHEQ,THE0Q
	55	REAL,DIMENSION(ITB,JTB),PRIVATE,SAVE :: PTBL
	56	REAL,DIMENSION(JTB,ITB),PRIVATE,SAVE :: TTBL
	57	REAL,DIMENSION(JTBQ,ITBQ),PRIVATE,SAVE :: TTBLQ
	58
	59	!*** SHARE COPIES FOR MODULE_BL_MYJPBL
	60	!
	61	REAL,DIMENSION(JTB) :: QS0_EXP,SQS_EXP
	62	REAL,DIMENSION(ITB,JTB) :: PTBL_EXP
	63	!
	64	REAL,PARAMETER :: RDP=(ITB-1.)/(PH-PL),RDPQ=(ITBQ-1.)/(PH-PLQ) &
	65	& ,RDQ=ITB-1,RDTH=(JTB-1.)/(THH-THL) &
	66	& ,RDTHE=JTB-1.,RDTHEQ=JTBQ-1. &
	67	& ,RSFCP=1./101300.
	68	!
	69	REAL,PARAMETER :: AVGEFI=(EFIMN+1.)*0.5
	70	!
	71	!-----------------------------------------------------------------------
	72	!
	73	CONTAINS
	74	!
	75	!-----------------------------------------------------------------------
	76	SUBROUTINE BMJDRV( &
	77	& IDS,IDE,JDS,JDE,KDS,KDE &
	78	& ,IMS,IME,JMS,JME,KMS,KME &
	79	& ,ITS,ITE,JTS,JTE,KTS,KTE &
	80	& ,DT,ITIMESTEP,STEPCU &
	81	& ,CUDT, CURR_SECS, ADAPT_STEP_FLAG &
	82	& ,RAINCV,PRATEC,CUTOP,CUBOT,KPBL &
	83	& ,TH,T,QV &
	84	& ,PINT,PMID,PI,RHO,DZ8W &
	85	& ,CP,R,ELWV,ELIV,G,TFRZ,D608 &
	86	& ,CLDEFI,LOWLYR,XLAND,CU_ACT_FLAG &
	87	! optional
	88	& ,RTHCUTEN, RQVCUTEN &
	89	& )
	90	!-----------------------------------------------------------------------
	91	IMPLICIT NONE
	92	!-----------------------------------------------------------------------
	93	INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
	94	& ,IMS,IME,JMS,JME,KMS,KME &
	95	& ,ITS,ITE,JTS,JTE,KTS,KTE
	96	!
	97	INTEGER,INTENT(IN) :: ITIMESTEP,STEPCU
	98	!
	99	INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: KPBL,LOWLYR
	100	!
	101	REAL,INTENT(IN) :: CP,DT,ELIV,ELWV,G,R,TFRZ,D608
	102	!
	103	REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: XLAND
	104	!
	105	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: DZ8W &
	106	& ,PI,PINT &
	107	& ,PMID,QV &
	108	& ,RHO,T,TH
	109	!
	110	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME) &
	111	& ,OPTIONAL &
	112	& ,INTENT(INOUT) :: RQVCUTEN,RTHCUTEN
	113	!
	114	REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: CLDEFI,RAINCV, &
	115	PRATEC
	116	!
	117	REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: CUBOT,CUTOP
	118	!
	119	LOGICAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: CU_ACT_FLAG
	120
	121	! Adaptive time-step variables
	122	REAL, INTENT(IN ) :: CUDT
	123	REAL, INTENT(IN ) :: CURR_SECS
	124	LOGICAL,INTENT(IN ) :: ADAPT_STEP_FLAG
	125	!
	126	!-----------------------------------------------------------------------
	127	!***
	128	!*** LOCAL VARIABLES
	129	!***
	130	!-----------------------------------------------------------------------
	131	INTEGER :: LBOT,LPBL,LTOP
	132	!
	133	REAL :: DTCNVC,LANDMASK,PCPCOL,PSFC,PTOP
	134	!
	135	REAL,DIMENSION(KTS:KTE) :: DPCOL,DQDT,DTDT,PCOL,QCOL,TCOL
	136	!
	137	INTEGER :: I,J,K,KFLIP,LMH
	138
	139	LOGICAL :: run_param
	140	!
	141	!*** Begin debugging convection
	142	REAL :: DELQ,DELT,PLYR
	143	INTEGER :: IMD,JMD
	144	LOGICAL :: PRINT_DIAG
	145	!*** End debugging convection
	146	!
	147	!-----------------------------------------------------------------------
	148	!***********************************************************************
	149	!-----------------------------------------------------------------------
	150	!
	151	!*** PREPARE TO CALL BMJ CONVECTION SCHEME
	152	!
	153	!-----------------------------------------------------------------------
	154	!
	155	!*** CHECK TO SEE IF THIS IS A CONVECTION TIMESTEP
	156	!
	157	if (adapt_step_flag) then
	158	if ( (ITIMESTEP .eq. 0) .or. (cudt .eq. 0) .or. &
	159	( CURR_SECS + dt >= ( int( CURR_SECS / ( cudt * 60 ) ) + 1 ) * cudt * 60 ) ) then
	160	run_param = .TRUE.
	161	else
	162	run_param = .FALSE.
	163	endif
	164
	165	else
	166	if (MOD(ITIMESTEP,STEPCU) .EQ. 0 .or. ITIMESTEP .eq. 0) then
	167	run_param = .TRUE.
	168	else
	169	run_param = .FALSE.
	170	endif
	171	endif
	172
	173	!-----------------------------------------------------------------------
	174	!
	175	!*** COMPUTE CONVECTION EVERY STEPCU*DT/60.0 MINUTES
	176	!
	177	!*** Begin debugging convection
	178	IMD=(IMS+IME)/2
	179	JMD=(JMS+JME)/2
	180	PRINT_DIAG=.FALSE.
	181	!*** End debugging convection
	182
	183	IF(run_param)THEN
	184	!
	185	DO J=JTS,JTE
	186	DO I=ITS,ITE
	187	CU_ACT_FLAG(I,J)=.TRUE.
	188	ENDDO
	189	ENDDO
	190
	191	!
	192	DTCNVC=DT*STEPCU
	193	!
	194	DO J=JTS,JTE
	195	DO I=ITS,ITE
	196	!
	197	DO K=KTS,KTE
	198	DQDT(K)=0.
	199	DTDT(K)=0.
	200	ENDDO
	201	!
	202	RAINCV(I,J)=0.
	203	PRATEC(I,J)=0.
	204	PCPCOL=0.
	205	PSFC=PINT(I,LOWLYR(I,J),J)
	206	PTOP=PINT(I,KTE+1,J) ! KTE+1=KME
	207	!
	208	!*** CONVERT TO BMJ LAND MASK (1.0 FOR SEA; 0.0 FOR LAND)
	209	!
	210	LANDMASK=XLAND(I,J)-1.
	211	!
	212	!*** FILL 1-D VERTICAL ARRAYS
	213	!*** AND FLIP DIRECTION SINCE BMJ SCHEME
	214	!*** COUNTS DOWNWARD FROM THE DOMAIN'S TOP
	215	!
	216	DO K=KTS,KTE
	217	KFLIP=KTE+1-K
	218	!
	219	!*** CONVERT FROM MIXING RATIO TO SPECIFIC HUMIDITY
	220	!
	221	QCOL(K)=MAX(EPSQ,QV(I,KFLIP,J)/(1.+QV(I,KFLIP,J)))
	222	TCOL(K)=T(I,KFLIP,J)
	223	PCOL(K)=PMID(I,KFLIP,J)
	224	! DPCOL(K)=PINT(I,KFLIP,J)-PINT(I,KFLIP+1,J)
	225	DPCOL(K)=RHO(I,KFLIP,J)GDZ8W(I,KFLIP,J)
	226	ENDDO
	227	!
	228	!*** LOWEST LAYER ABOVE GROUND MUST ALSO BE FLIPPED
	229	!
	230	LMH=KTE+1-LOWLYR(I,J)
	231	LPBL=KTE+1-KPBL(I,J)
	232	!-----------------------------------------------------------------------
	233	!***
	234	!*** CALL CONVECTION
	235	!***
	236	!-----------------------------------------------------------------------
	237	!*** Begin debugging convection
	238	! PRINT_DIAG=.FALSE.
	239	! IF(I==IMD.AND.J==JMD)PRINT_DIAG=.TRUE.
	240	!*** End debugging convection
	241	!-----------------------------------------------------------------------
	242	CALL BMJ(ITIMESTEP,I,J,DTCNVC,LMH,LANDMASK,CLDEFI(I,J) &
	243	& ,DPCOL,PCOL,QCOL,TCOL,PSFC,PTOP &
	244	& ,DQDT,DTDT,PCPCOL,LBOT,LTOP,LPBL &
	245	& ,CP,R,ELWV,ELIV,G,TFRZ,D608 &
	246	& ,PRINT_DIAG &
	247	& ,IDS,IDE,JDS,JDE,KDS,KDE &
	248	& ,IMS,IME,JMS,JME,KMS,KME &
	249	& ,ITS,ITE,JTS,JTE,KTS,KTE)
	250	!-----------------------------------------------------------------------
	251	!
	252	!*** COMPUTE HEATING AND MOISTENING TENDENCIES
	253	!
	254	IF ( PRESENT( RTHCUTEN ) .AND. PRESENT( RQVCUTEN )) THEN
	255	DO K=KTS,KTE
	256	KFLIP=KTE+1-K
	257	RTHCUTEN(I,K,J)=DTDT(KFLIP)/PI(I,K,J)
	258	!
	259	!*** CONVERT FROM SPECIFIC HUMIDTY BACK TO MIXING RATIO
	260	!
	261	RQVCUTEN(I,K,J)=DQDT(KFLIP)/(1.-QCOL(KFLIP))**2
	262	ENDDO
	263	ENDIF
	264	!
	265	!*** ALL UNITS IN BMJ SCHEME ARE MKS, THUS CONVERT PRECIP FROM METERS
	266	!*** TO MILLIMETERS PER STEP FOR OUTPUT.
	267	!
	268	RAINCV(I,J)=PCPCOL*1.E3/STEPCU
	269	PRATEC(I,J)=PCPCOL1.E3/(STEPCU DT)
	270	!
	271	!*** CONVECTIVE CLOUD TOP AND BOTTOM FROM THIS CALL
	272	!
	273	CUTOP(I,J)=REAL(KTE+1-LTOP)
	274	CUBOT(I,J)=REAL(KTE+1-LBOT)
	275	!
	276	!-----------------------------------------------------------------------
	277	!*** Begin debugging convection
	278	IF(PRINT_DIAG)THEN
	279	DELT=0.
	280	DELQ=0.
	281	PLYR=0.
	282	IF(LBOT>0.AND.LTOP<LBOT)THEN
	283	DO K=LTOP,LBOT
	284	PLYR=PLYR+DPCOL(K)
	285	DELQ=DELQ+DPCOL(K)DTCNVCABS(DQDT(K))
	286	DELT=DELT+DPCOL(K)DTCNVCABS(DTDT(K))
	287	ENDDO
	288	DELQ=DELQ/PLYR
	289	DELT=DELT/PLYR
	290	ENDIF
	291	!
	292	WRITE(6,"(2a,2i4,3e12.4,f7.2,4i3)") &
	293	'{cu3 i,j,PCPCOL,DTavg,DQavg,PLYR,' &
	294	,'LBOT,LTOP,CUBOT,CUTOP = ' &
	295	,i,j, PCPCOL,DELT,1000.DELQ,.01PLYR &
	296	,LBOT,LTOP,NINT(CUBOT(I,J)),NINT(CUTOP(I,J))
	297	!
	298	IF(PLYR> 0.)THEN
	299	DO K=LBOT,LTOP,-1
	300	KFLIP=KTE+1-K
	301	WRITE(6,"(a,i3,2e12.4,f7.2)") '{cu3a KFLIP,DT,DQ,DP = ' &
	302	,KFLIP,DTCNVCDTDT(K),1000.DTCNVCDQDT(K),.01DPCOL(K)
	303	ENDDO
	304	ENDIF
	305	ENDIF
	306	!*** End debugging convection
	307	!-----------------------------------------------------------------------
	308	!
	309	ENDDO
	310	ENDDO
	311	!
	312	ENDIF
	313	!
	314	END SUBROUTINE BMJDRV
	315	!-----------------------------------------------------------------------
	316	!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
	317	!-----------------------------------------------------------------------
	318	SUBROUTINE BMJ &
	319	!-----------------------------------------------------------------------
	320	& (ITIMESTEP,I,J,DTCNVC,LMH,SM,CLDEFI &
	321	& ,DPRS,PRSMID,Q,T,PSFC,PT &
	322	& ,DQDT,DTDT,PCPCOL,LBOT,LTOP,LPBL &
	323	& ,CP,R,ELWV,ELIV,G,TFRZ,D608 &
	324	& ,PRINT_DIAG &
	325	& ,IDS,IDE,JDS,JDE,KDS,KDE &
	326	& ,IMS,IME,JMS,JME,KMS,KME &
	327	& ,ITS,ITE,JTS,JTE,KTS,KTE)
	328	!-----------------------------------------------------------------------
	329	IMPLICIT NONE
	330	!-----------------------------------------------------------------------
	331	INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
	332	,IMS,IME,JMS,JME,KMS,KME &
	333	,ITS,ITE,JTS,JTE,KTS,KTE &
	334	,I,J,ITIMESTEP
	335	!
	336	INTEGER,INTENT(IN) :: LMH,LPBL
	337	!
	338	INTEGER,INTENT(OUT) :: LBOT,LTOP
	339	!
	340	REAL,INTENT(IN) :: CP,D608,DTCNVC,ELIV,ELWV,G,PSFC,PT,R,SM,TFRZ
	341	!
	342	REAL,DIMENSION(KTS:KTE),INTENT(IN) :: DPRS,PRSMID,Q,T
	343	!
	344	REAL,INTENT(INOUT) :: CLDEFI,PCPCOL
	345	!
	346	REAL,DIMENSION(KTS:KTE),INTENT(INOUT) :: DQDT,DTDT
	347	!
	348	!-----------------------------------------------------------------------
	349	!*** DEFINE LOCAL VARIABLES
	350	!-----------------------------------------------------------------------
	351	!
	352	REAL,DIMENSION(KTS:KTE) :: APEK,APESK,EL,FPK &
	353	,PK,PSK,QK,QREFK,QSATK &
	354	,THERK,THEVRF,THSK &
	355	,THVMOD,THVREF,TK,TREFK
	356	REAL,DIMENSION(KTS:KTE) :: APE,DIFQ,DIFT,THEE,THES,TREF
	357	!
	358	REAL,DIMENSION(KTS:KTE) :: CPE,CPEcnv,DTV,DTVcnv,THEScnv !<-- CPE for shallow convection buoyancy check (24 Aug 2006)
	359	!
	360	LOGICAL :: DEEP,SHALLOW
	361	!
	362	!*** Begin debugging convection
	363	LOGICAL :: PRINT_DIAG
	364	!*** End debugging convection
	365	!
	366	!-----------------------------------------------------------------------
	367	!***
	368	!*** LOCAL SCALARS
	369	!***
	370	!-----------------------------------------------------------------------
	371	REAL :: APEKL,APEKXX,APEKXY,APES,APESTS &
	372	& ,AVRGT,AVRGTL,BQ,BQK,BQS00K,BQS10K &
	373	& ,CAPA,CUP,DEN,DENTPY,DEPMIN,DEPTH &
	374	& ,DEPWL,DHDT,DIFQL,DIFTL,DP,DPKL,DPLO,DPMIX,DPOT &
	375	& ,DPUP,DQREF,DRHDP,DRHEAT,DSP &
	376	& ,DSP0,DSP0K,DSPB,DSPBK,DSPT,DSPTK &
	377	& ,DSQ,DST,DSTQ,DTHEM,DTDP,EFI &
	378	& ,FEFI,FFUP,FPRS,FPTK,HCORR &
	379	& ,OTSUM,P,P00K,P01K,P10K,P11K &
	380	& ,PART1,PART2,PART3,PBOT,PBOTFC,PBTK &
	381	& ,PK0,PKB,PKL,PKT,PKXXXX,PKXXXY &
	382	& ,PLMH,PELEVFC,PBTmx,plo,POTSUM,PP1,PPK,PRECK &
	383	& ,PRESK,PSP,PSUM,PTHRS,PTOP,PTPK,PUP &
	384	& ,QBT,QKL,QNEW,QOTSUM,QQ1,QQK,QRFKL &
	385	& ,QRFTP,QSP,QSUM,QUP,RDP0T &
	386	& ,RDPSUM,RDTCNVC,RHH,RHL,RHMAX,ROTSUM,RTBAR,RHAVG &
	387	& ,SM1,SMIX,SQ,SQK,SQS00K,SQS10K,STABDL,SUMDE,SUMDP &
	388	& ,SUMDT,TAUK,TAUKSC,TCORR,THBT,THERKX,THERKY &
	389	& ,THESP,THSKL,THTPK,THVMKL,TKL,TNEW &
	390	& ,TQ,TQK,TREFKX,TRFKL,trmlo,trmup,TSKL,tsp,TTH &
	391	& ,TTHK,TUP &
	392	& ,CAPEcnv,PSPcnv,THBTcnv,CAPEtrigr,CAPE &
	393	& ,TLEV2,QSAT1,QSAT2,RHSHmax
	394	!
	395	INTEGER :: IQ,IQTB,IT,ITER,ITREFI,ITTB,ITTBK,KB,KNUMH,KNUML &
	396	& ,L,L0,L0M1,LB,LBM1,LCOR,LPT1 &
	397	& ,LQM,LSHU,LTP1,LTP2,LTSH, LBOTcnv,LTOPcnv,LMID
	398	!-----------------------------------------------------------------------
	399	!
	400	REAL,PARAMETER :: DSPBSL=DSPBFLFSL,DSP0SL=DSP0FLFSL &
	401	& ,DSPTSL=DSPTFL*FSL &
	402	& ,DSPBSS=DSPBFSFSS,DSP0SS=DSP0FSFSS &
	403	& ,DSPTSS=DSPTFS*FSS
	404	!
	405	REAL,PARAMETER :: ELEVFC=0.6,STEFI=1.
	406	!
	407	REAL,PARAMETER :: SLOPBL=(DSPBFL-DSPBSL)/(1.-EFIMN) &
	408	& ,SLOP0L=(DSP0FL-DSP0SL)/(1.-EFIMN) &
	409	& ,SLOPTL=(DSPTFL-DSPTSL)/(1.-EFIMN) &
	410	& ,SLOPBS=(DSPBFS-DSPBSS)/(1.-EFIMN) &
	411	& ,SLOP0S=(DSP0FS-DSP0SS)/(1.-EFIMN) &
	412	& ,SLOPTS=(DSPTFS-DSPTSS)/(1.-EFIMN) &
	413	& ,SLOPST=(STABDF-STABDS)/(1.-EFIMN) &
	414	& ,SLOPE=(1.-EFMNT)/(1.-EFIMN)
	415	!
	416	REAL :: A23M4L,CPRLG,ELOCP,RCP,QWAT
	417	!-----------------------------------------------------------------------
	418	!***********************************************************************
	419	!-----------------------------------------------------------------------
	420	CAPA=R/CP
	421	CPRLG=CP/(ROWGELWV)
	422	ELOCP=ELIWV/CP
	423	RCP=1./CP
	424	A23M4L=A2(A3-A4)ELWV
	425	RDTCNVC=1./DTCNVC
	426	DEPMIN=PSHPSFCRSFCP
	427	!
	428	DEEP=.FALSE.
	429	SHALLOW=.FALSE.
	430	!
	431	DSP0=0.
	432	DSPB=0.
	433	DSPT=0.
	434	!-----------------------------------------------------------------------
	435	TAUK=DTCNVC/TREL
	436	TAUKSC=DTCNVC/(1.0*TREL)
	437	!-----------------------------------------------------------------------
	438	!-----------------------------PREPARATIONS------------------------------
	439	!-----------------------------------------------------------------------
	440	LBOT=LMH
	441	PCPCOL=0.
	442	TREF(KTS)=T(KTS)
	443	!
	444	DO L=KTS,LMH
	445	APESTS=PRSMID(L)
	446	APE(L)=(1.E5/APESTS)**CAPA
	447	CPEcnv(L)=0.
	448	DTVcnv(L)=0.
	449	THEScnv(L)=0.
	450	ENDDO
	451	!
	452	!-----------------------------------------------------------------------
	453	!----------------SEARCH FOR MAXIMUM BUOYANCY LEVEL----------------------
	454	!-----------------------------------------------------------------------
	455	!
	456	PLMH=PRSMID(LMH)
	457	PELEVFC=PLMH*ELEVFC
	458	PBTmx=PRSMID(LMH)-PONE
	459	CAPEcnv=0.
	460	PSPcnv=0.
	461	THBTcnv=0.
	462	LBOTcnv=LBOT
	463	LTOPcnv=LBOT
	464	!
	465	!-----------------------------------------------------------------------
	466	!----------------TRIAL MAXIMUM BUOYANCY LEVEL VARIABLES-----------------
	467	!-----------------------------------------------------------------------
	468	!
	469	max_buoy_loop: DO KB=LMH,1,-1
	470	!
	471	!-----------------------------------------------------------------------
	472	!
	473	PKL=PRSMID(KB)
	474	! IF (PKL<PELEVFC .OR. T(KB)<=TFRZ) EXIT
	475	IF (PKL<PELEVFC) EXIT
	476	LBOT=LMH
	477	LTOP=LMH
	478	!
	479	!-----------------------------------------------------------------------
	480	!*** SEARCH OVER A SCALED DEPTH IN FINDING THE PARCEL
	481	!*** WITH THE MAX THETA-E
	482	!-----------------------------------------------------------------------
	483	!
	484	QBT=Q(KB)
	485	THBT=T(KB)*APE(KB)
	486	TTH=(THBT-THL)*RDTH
	487	QQ1=TTH-AINT(TTH)
	488	ITTB=INT(TTH)+1
	489	!----------------KEEPING INDICES WITHIN THE TABLE-----------------------
	490	IF(ITTB<1)THEN
	491	ITTB=1
	492	QQ1=0.
	493	ELSE IF(ITTB>=JTB)THEN
	494	ITTB=JTB-1
	495	QQ1=0.
	496	ENDIF
	497	!--------------BASE AND SCALING FACTOR FOR SPEC. HUMIDITY---------------
	498	ITTBK=ITTB
	499	BQS00K=QS0(ITTBK)
	500	SQS00K=SQS(ITTBK)
	501	BQS10K=QS0(ITTBK+1)
	502	SQS10K=SQS(ITTBK+1)
	503	!--------------SCALING SPEC. HUMIDITY & TABLE INDEX---------------------
	504	BQ=(BQS10K-BQS00K)*QQ1+BQS00K
	505	SQ=(SQS10K-SQS00K)*QQ1+SQS00K
	506	TQ=(QBT-BQ)/SQ*RDQ
	507	PP1=TQ-AINT(TQ)
	508	IQTB=INT(TQ)+1
	509	!----------------KEEPING INDICES WITHIN THE TABLE-----------------------
	510	IF(IQTB<1)THEN
	511	IQTB=1
	512	PP1=0.
	513	ELSE IF(IQTB>=ITB)THEN
	514	IQTB=ITB-1
	515	PP1=0.
	516	ENDIF
	517	!--------------SATURATION PRESSURE AT FOUR SURROUNDING TABLE PTS.-------
	518	IQ=IQTB
	519	IT=ITTB
	520	P00K=PTBL(IQ ,IT )
	521	P10K=PTBL(IQ+1,IT )
	522	P01K=PTBL(IQ ,IT+1)
	523	P11K=PTBL(IQ+1,IT+1)
	524	!
	525	!--------------SATURATION POINT VARIABLES AT THE BOTTOM-----------------
	526	!
	527	PSP=P00K+(P10K-P00K)PP1+(P01K-P00K)QQ1 &
	528	& +(P00K-P10K-P01K+P11K)PP1QQ1
	529	APES=(1.E5/PSP)**CAPA
	530	THESP=THBTEXP(ELOCPQBT*APES/THBT)
	531	!
	532	!-----------------------------------------------------------------------
	533	!-----------CHOOSE CLOUD BASE AS MODEL LEVEL JUST BELOW PSP-------------
	534	!-----------------------------------------------------------------------
	535	!
	536	DO L=KTS,LMH-1
	537	P=PRSMID(L)
	538	IF(P<PSP.AND.P>=PQM)LBOT=L+1
	539	ENDDO
	540	!***
	541	!*** WARNING: LBOT MUST NOT BE > LMH-1 IN SHALLOW CONVECTION
	542	!*** MAKE SURE CLOUD BASE IS AT LEAST PONE ABOVE THE SURFACE
	543	!***
	544	PBOT=PRSMID(LBOT)
	545	IF(PBOT>=PBTmx.OR.LBOT>=LMH)THEN
	546	DO L=KTS,LMH-1
	547	P=PRSMID(L)
	548	IF(P<PBTmx)LBOT=L
	549	ENDDO
	550	PBOT=PRSMID(LBOT)
	551	ENDIF
	552	!
	553	!-----------------------------------------------------------------------
	554	!----------------CLOUD TOP COMPUTATION----------------------------------
	555	!-----------------------------------------------------------------------
	556	!
	557	LTOP=LBOT
	558	PTOP=PBOT
	559	DO L=LMH,KTS,-1
	560	THES(L)=THESP
	561	ENDDO
	562	!
	563	!-----------------------------------------------------------------------
	564	!### BEGIN: ########### WARNING ########### WARNING ###########
	565	!-----------------------------------------------------------------------
	566	!
	567	!### IMPORTANT: THIS "DO KB=LMH,1,-1" loop must be broken up into two
	568	! separate loops in order for entrainment as programmed below to work
	569	! properly.
	570	!
	571	!--------------- ENTRAINMENT DURING PARCEL ASCENT --------------------
	572	!
	573	! DO L=LMH,KB,-1
	574	! THES(L)=THESP
	575	! ENDDO
	576	!
	577	! DO L=KTS,LMH
	578	! THEE(L)=THES(L)
	579	! ENDDO
	580	!!
	581	! FEFI=(CLDEFI-EFIMN)*SLOPE+EFMNT
	582	! FFUP=FUP/(FEFI*FEFI)
	583	!!
	584	! IF(PBOT>ENPLO)THEN
	585	! FPRS=1.
	586	! ELSEIF(PBOT>ENPUP)THEN
	587	! FPRS=(PBOT-ENPUP)*RENDP
	588	! ELSE
	589	! FPRS=0.
	590	! ENDIF
	591	!!
	592	! FFUP=FFUPFPRSFPRS*0.5
	593	! DPUP=DPRS(KB)
	594	!!
	595	! DO L=KB-1,KTS,-1
	596	! DPLO=DPUP
	597	! DPUP=DPRS(L)
	598	!!
	599	! THES(L)=((-FFUPDPLO+1.)THES(L+1) &
	600	! & +(THEE(L)DPUP+THEE(L+1)DPLO)*FFUP) &
	601	! & /(FFUP*DPUP+1.)
	602	! ENDDO
	603	!
	604	!-----------------------------------------------------------------------
	605	!### END: ########### WARNING ########### WARNING ###########
	606	!-----------------------------------------------------------------------
	607	!!
	608	!-----------------------------------------------------------------------
	609	!!*** COMPUTE PARCEL TEMPERATURE ALONG THE ASCENT TRAJECTORY
	610	!!*** SCALING PRESSURE & TT TABLE INDEX.
	611	!-----------------------------------------------------------------------
	612	!!
	613	!!
	614	! DO L=LMH,KTS,-1
	615	!!
	616	! PRESK=PRSMID(L)
	617	!!
	618	! IF(PRESK<PLQ)THEN
	619	! CALL TTBLEX(ITB,JTB,PL,PRESK,RDP,RDTHE &
	620	! & ,STHE,THE0,THES(L),TTBL,TREF(L))
	621	! ELSE
	622	! CALL TTBLEX(ITBQ,JTBQ,PLQ,PRESK,RDPQ,RDTHEQ &
	623	! & ,STHEQ,THE0Q,THES(L),TTBLQ,TREF(L))
	624	! ENDIF
	625	!!
	626	! ENDDO
	627	!!
	628	!!-----------------------------------------------------------------------
	629	!!----------------BUOYANCY CHECK-----------------------------------------
	630	!!-----------------------------------------------------------------------
	631	!!
	632	! DO L=LMH,KTS,-1
	633	! IF(TREF(L)>T(L)-DTTOP)LTOP=L
	634	! ENDDO
	635	!!
	636	!!*** CLOUD TOP PRESSURE
	637	!!
	638	! PTOP=PRSMID(LTOP)
	639	!
	640	!------------------FIRST ENTROPY CHECK----------------------------------
	641	!
	642	DO L=KTS,LMH
	643	CPE(L)=0.
	644	DTV(L)=0.
	645	ENDDO
	646	!-----------------------------------------------------------------------
	647	! lbot_ltop: IF(LBOT>LTOP)THEN
	648	!-----------------------------------------------------------------------
	649	!-- Begin: Buoyancy check including deep convection (24 Aug 2006)
	650	!-----------------------------------------------------------------------
	651	DENTPY=0.
	652	L=KB
	653	PLO=PRSMID(L)
	654	TRMLO=0.
	655	CAPEtrigr=DTPtrigr/T(LBOT)
	656	!
	657	!--- Below cloud base
	658	!
	659	IF(KB>LBOT) THEN
	660	DO L=KB-1,LBOT+1,-1
	661	PUP=PRSMID(L)
	662	TUP=THBT/APE(L)
	663	DP=PLO-PUP
	664	TRMUP=(TUP(QBT0.608+1.) &
	665	& -T(L)(Q(L)0.608+1.))*0.5 &
	666	& /(T(L)(Q(L)0.608+1.))
	667	DTV(L)=TRMLO+TRMUP
	668	DENTPY=DTV(L)*DP+DENTPY
	669	CPE(L)=DENTPY
	670	IF (DENTPY < CAPEtrigr) GO TO 170
	671	PLO=PUP
	672	TRMLO=TRMUP
	673	ENDDO
	674	ELSE
	675	L=LBOT+1
	676	PLO=PRSMID(L)
	677	TUP=THBT/APE(L)
	678	TRMLO=(TUP(QBT0.608+1.) &
	679	& -T(L)(Q(L)0.608+1.))*0.5 &
	680	& /(T(L)(Q(L)0.608+1.))
	681	ENDIF ! IF(KB>LBOT) THEN
	682	!
	683	!--- At cloud base
	684	!
	685	L=LBOT
	686	PUP=PSP
	687	TUP=THBT/APES
	688	TSP=(T(L+1)-T(L))/(PLO-PBOT) &
	689	& *(PUP-PBOT)+T(L)
	690	QSP=(Q(L+1)-Q(L))/(PLO-PBOT) &
	691	& *(PUP-PBOT)+Q(L)
	692	DP=PLO-PUP
	693	TRMUP=(TUP(QBT0.608+1.) &
	694	& -TSP(QSP0.608+1.))*0.5 &
	695	& /(TSP(QSP0.608+1.))
	696	DTV(L)=TRMLO+TRMUP
	697	DENTPY=DTV(L)*DP+DENTPY
	698	CPE(L)=DENTPY
	699	DTV(L)=DTV(L)*DP
	700	PLO=PUP
	701	TRMLO=TRMUP
	702	PUP=PRSMID(L)
	703	!
	704	!--- Calculate updraft temperature along moist adiabat (TUP)
	705	!
	706	IF(PUP<PLQ)THEN
	707	CALL TTBLEX(ITB,JTB,PL,PUP,RDP,RDTHE &
	708	& ,STHE,THE0,THES(L),TTBL,TUP)
	709	ELSE
	710	CALL TTBLEX(ITBQ,JTBQ,PLQ,PUP,RDPQ,RDTHEQ &
	711	& ,STHEQ,THE0Q,THES(L),TTBLQ,TUP)
	712	ENDIF
	713	!
	714	QUP=PQ0/PUPEXP(A2(TUP-A3)/(TUP-A4))
	715	QWAT=QBT-QUP !-- Water loading effects, reversible adiabat
	716	DP=PLO-PUP
	717	TRMUP=(TUP(QUP0.608+1.-QWAT) &
	718	& -T(L)(Q(L)0.608+1.))*0.5 &
	719	& /(T(L)(Q(L)0.608+1.))
	720	DENTPY=(TRMLO+TRMUP)*DP+DENTPY
	721	CPE(L)=DENTPY
	722	DTV(L)=(DTV(L)+(TRMLO+TRMUP)*DP)/(PRSMID(LBOT+1)-PRSMID(LBOT))
	723	!
	724	IF (DENTPY < CAPEtrigr) GO TO 170
	725	!
	726	PLO=PUP
	727	TRMLO=TRMUP
	728	!
	729	!-----------------------------------------------------------------------
	730	!--- In cloud above cloud base
	731	!-----------------------------------------------------------------------
	732	!
	733	DO L=LBOT-1,KTS,-1
	734	PUP=PRSMID(L)
	735	!
	736	!--- Calculate updraft temperature along moist adiabat (TUP)
	737	!
	738	IF(PUP<PLQ)THEN
	739	CALL TTBLEX(ITB,JTB,PL,PUP,RDP,RDTHE &
	740	& ,STHE,THE0,THES(L),TTBL,TUP)
	741	ELSE
	742	CALL TTBLEX(ITBQ,JTBQ,PLQ,PUP,RDPQ,RDTHEQ &
	743	& ,STHEQ,THE0Q,THES(L),TTBLQ,TUP)
	744	ENDIF
	745	!
	746	QUP=PQ0/PUPEXP(A2(TUP-A3)/(TUP-A4))
	747	QWAT=QBT-QUP !-- Water loading effects, reversible adiabat
	748	DP=PLO-PUP
	749	TRMUP=(TUP(QUP0.608+1.-QWAT) &
	750	& -T(L)(Q(L)0.608+1.))*0.5 &
	751	& /(T(L)(Q(L)0.608+1.))
	752	DTV(L)=TRMLO+TRMUP
	753	DENTPY=DTV(L)*DP+DENTPY
	754	CPE(L)=DENTPY
	755	!
	756	IF (DENTPY < CAPEtrigr) GO TO 170
	757	!
	758	PLO=PUP
	759	TRMLO=TRMUP
	760	ENDDO
	761	!
	762	!-----------------------------------------------------------------------
	763	!
	764	170 LTP1=KB
	765	CAPE=0.
	766	!
	767	!-----------------------------------------------------------------------
	768	!--- Cloud top level (LTOP) is located where CAPE is a maximum
	769	!--- Exit cloud-top search if CAPE < CAPEtrigr
	770	!-----------------------------------------------------------------------
	771	!
	772	DO L=KB,KTS,-1
	773	IF (CPE(L) < CAPEtrigr) THEN
	774	EXIT
	775	ELSE IF (CPE(L) > CAPE) THEN
	776	LTP1=L
	777	CAPE=CPE(L)
	778	ENDIF
	779	ENDDO !-- End DO L=KB,KTS,-1
	780	!
	781	LTOP=MIN(LTP1,LBOT)
	782	!
	783	!-----------------------------------------------------------------------
	784	!--------------- CHECK FOR MAXIMUM INSTABILITY ------------------------
	785	!-----------------------------------------------------------------------
	786	IF(CAPE > CAPEcnv) THEN
	787	CAPEcnv=CAPE
	788	PSPcnv=PSP
	789	THBTcnv=THBT
	790	LBOTcnv=LBOT
	791	LTOPcnv=LTOP
	792	DO L=LMH,KTS,-1
	793	CPEcnv(L)=CPE(L)
	794	DTVcnv(L)=DTV(L)
	795	THEScnv(L)=THES(L)
	796	ENDDO
	797	ENDIF ! End IF(CAPE > CAPEcnv) THEN
	798	!
	799	! ENDIF lbot_ltop
	800	!
	801	!-----------------------------------------------------------------------
	802	!
	803	ENDDO max_buoy_loop
	804	!
	805	!-----------------------------------------------------------------------
	806	!------------------------ MAXIMUM INSTABILITY ------------------------
	807	!-----------------------------------------------------------------------
	808	!
	809	IF(CAPEcnv > 0.) THEN
	810	PSP=PSPcnv
	811	THBT=THBTcnv
	812	LBOT=LBOTcnv
	813	LTOP=LTOPcnv
	814	PBOT=PRSMID(LBOT)
	815	PTOP=PRSMID(LTOP)
	816	!
	817	DO L=LMH,KTS,-1
	818	CPE(L)=CPEcnv(L)
	819	DTV(L)=DTVcnv(L)
	820	THES(L)=THEScnv(L)
	821	ENDDO
	822	!
	823	ENDIF
	824	!
	825	!-----------------------------------------------------------------------
	826	!----- Quick exit if cloud is too thin or no CAPE is present ---------
	827	!-----------------------------------------------------------------------
	828	!
	829	IF(PTOP>PBOT-PNO.OR.LTOP>LBOT-2.OR.CAPEcnv<=0.)THEN
	830	LBOT=0
	831	LTOP=KTE
	832	PBOT=PRSMID(LMH)
	833	PTOP=PBOT
	834	CLDEFI=AVGEFISM+STEFI(1.-SM)
	835	GO TO 800
	836	ENDIF
	837	!
	838	!*** DEPTH OF CLOUD REQUIRED TO MAKE THE POINT A DEEP CONVECTION POINT
	839	!*** IS A SCALED VALUE OF PSFC.
	840	!
	841	DEPTH=PBOT-PTOP
	842	!
	843	IF(DEPTH>=DEPMIN) THEN
	844	DEEP=.TRUE.
	845	ELSE
	846	SHALLOW=.TRUE.
	847	GO TO 600
	848	ENDIF
	849	!
	850	!-----------------------------------------------------------------------
	851	!DCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCD
	852	!DCDCDCDCDCDCDCDCDCDCDC DEEP CONVECTION DCDCDCDCDCDCDCDCDCDCDCDCDCD
	853	!DCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCD
	854	!-----------------------------------------------------------------------
	855	!
	856	300 CONTINUE
	857	!
	858	LB =LBOT
	859	EFI=CLDEFI
	860	!-----------------------------------------------------------------------
	861	!--------------INITIALIZE VARIABLES IN THE CONVECTIVE COLUMN------------
	862	!-----------------------------------------------------------------------
	863	!***
	864	!*** ONE SHOULD NOTE THAT THE VALUES ASSIGNED TO THE ARRAY TREFK
	865	!*** IN THE FOLLOWING LOOP ARE REALLY ONLY RELEVANT IN ANCHORING THE
	866	!*** REFERENCE TEMPERATURE PROFILE AT LEVEL LB. WHEN BUILDING THE
	867	!*** REFERENCE PROFILE FROM CLOUD BASE, THEN ASSIGNING THE
	868	!*** AMBIENT TEMPERATURE TO TREFK IS ACCEPTABLE. HOWEVER, WHEN
	869	!*** BUILDING THE REFERENCE PROFILE FROM SOME OTHER LEVEL (SUCH AS
	870	!*** ONE LEVEL ABOVE THE GROUND), THEN TREFK SHOULD BE FILLED WITH
	871	!*** THE TEMPERATURES IN TREF(L) WHICH ARE THE TEMPERATURES OF
	872	!*** THE MOIST ADIABAT THROUGH CLOUD BASE. BY THE TIME THE LINE
	873	!*** NUMBERED 450 HAS BEEN REACHED, TREFK ACTUALLY DOES HOLD THE
	874	!*** REFERENCE TEMPERATURE PROFILE.
	875	!***
	876	DO L=KTS,LMH
	877	DIFT (L)=0.
	878	DIFQ (L)=0.
	879	TKL =T(L)
	880	TK (L)=TKL
	881	TREFK (L)=TKL
	882	QKL =Q(L)
	883	QK (L)=QKL
	884	QREFK (L)=QKL
	885	PKL =PRSMID(L)
	886	PK (L)=PKL
	887	PSK (L)=PKL
	888	APEKL =APE(L)
	889	APEK (L)=APEKL
	890	!
	891	!--- Calculate temperature along moist adiabat (TREF)
	892	!
	893	IF(PKL<PLQ)THEN
	894	CALL TTBLEX(ITB,JTB,PL,PKL,RDP,RDTHE &
	895	& ,STHE,THE0,THES(L),TTBL,TREF(L))
	896	ELSE
	897	CALL TTBLEX(ITBQ,JTBQ,PLQ,PKL,RDPQ,RDTHEQ &
	898	& ,STHEQ,THE0Q,THES(L),TTBLQ,TREF(L))
	899	ENDIF
	900	THERK (L)=TREF(L)*APEKL
	901	ENDDO
	902	!
	903	!------------DEEP CONVECTION REFERENCE TEMPERATURE PROFILE------------
	904	!
	905	LTP1=LTOP+1
	906	LBM1=LB-1
	907	PKB=PK(LB)
	908	PKT=PK(LTOP)
	909	STABDL=(EFI-EFIMN)*SLOPST+STABDS
	910	!
	911	!------------TEMPERATURE REFERENCE PROFILE BELOW FREEZING LEVEL-------
	912	!
	913	EL(LB) = ELWV
	914	L0=LB
	915	PK0=PK(LB)
	916	TREFKX=TREFK(LB)
	917	THERKX=THERK(LB)
	918	APEKXX=APEK(LB)
	919	THERKY=THERK(LBM1)
	920	APEKXY=APEK(LBM1)
	921	!
	922	DO L=LBM1,LTOP,-1
	923	IF(T(L+1)<TFRZ)GO TO 430
	924	TREFKX=((THERKY-THERKX)*STABDL &
	925	& +TREFKX*APEKXX)/APEKXY
	926	TREFK(L)=TREFKX
	927	EL(L)=ELWV
	928	APEKXX=APEKXY
	929	THERKX=THERKY
	930	APEKXY=APEK(L-1)
	931	THERKY=THERK(L-1)
	932	L0=L
	933	PK0=PK(L0)
	934	ENDDO
	935	!
	936	!--------------FREEZING LEVEL AT OR ABOVE THE CLOUD TOP-----------------
	937	!
	938	GO TO 450
	939	!
	940	!--------------TEMPERATURE REFERENCE PROFILE ABOVE FREEZING LEVEL-------
	941	!
	942	430 L0M1=L0-1
	943	RDP0T=1./(PK0-PKT)
	944	DTHEM=THERK(L0)-TREFK(L0)*APEK(L0)
	945	!
	946	DO L=LTOP,L0M1
	947	TREFK(L)=(THERK(L)-(PK(L)-PKT)DTHEMRDP0T)/APEK(L)
	948	EL(L)=ELWV !ELIV
	949	ENDDO
	950	!
	951	!-----------------------------------------------------------------------
	952	!--------------DEEP CONVECTION REFERENCE HUMIDITY PROFILE---------------
	953	!-----------------------------------------------------------------------
	954	!
	955	!*** DEPWL IS THE PRESSURE DIFFERENCE BETWEEN CLOUD BASE AND
	956	!*** THE FREEZING LEVEL
	957	!
	958	450 CONTINUE
	959	DEPWL=PKB-PK0
	960	DEPTH=PFRZPSFCRSFCP
	961	SM1=1.-SM
	962	PBOTFC=1.
	963	!
	964	!-------------FIRST ADJUSTMENT OF TEMPERATURE PROFILE-------------------
	965	!!
	966	! SUMDT=0.
	967	! SUMDP=0.
	968	!!
	969	! DO L=LTOP,LB
	970	! SUMDT=(TK(L)-TREFK(L))*DPRS(L)+SUMDT
	971	! SUMDP=SUMDP+DPRS(L)
	972	! ENDDO
	973	!!
	974	! TCORR=SUMDT/SUMDP
	975	!!
	976	! DO L=LTOP,LB
	977	! TREFK(L)=TREFK(L)+TCORR
	978	! ENDDO
	979	!!
	980	!-----------------------------------------------------------------------
	981	!--------------- ITERATION LOOP FOR CLOUD EFFICIENCY -------------------
	982	!-----------------------------------------------------------------------
	983	!
	984	cloud_efficiency : DO ITREFI=1,ITREFI_MAX
	985	!
	986	!-----------------------------------------------------------------------
	987	DSPBK=((EFI-EFIMN)SLOPBS+DSPBSSPBOTFC)*SM &
	988	& +((EFI-EFIMN)SLOPBL+DSPBSLPBOTFC)*SM1
	989	DSP0K=((EFI-EFIMN)SLOP0S+DSP0SSPBOTFC)*SM &
	990	& +((EFI-EFIMN)SLOP0L+DSP0SLPBOTFC)*SM1
	991	DSPTK=((EFI-EFIMN)SLOPTS+DSPTSSPBOTFC)*SM &
	992	& +((EFI-EFIMN)SLOPTL+DSPTSLPBOTFC)*SM1
	993	!
	994	!-----------------------------------------------------------------------
	995	!
	996	DO L=LTOP,LB
	997	!
	998	!***
	999	!*** SATURATION PRESSURE DIFFERENCE
	1000	!***
	1001	IF(DEPWL>=DEPTH)THEN
	1002	IF(L<L0)THEN
	1003	DSP=((PK0-PK(L))DSPTK+(PK(L)-PKT)DSP0K)/(PK0-PKT)
	1004	ELSE
	1005	DSP=((PKB-PK(L))DSP0K+(PK(L)-PK0)DSPBK)/(PKB-PK0)
	1006	ENDIF
	1007	ELSE
	1008	DSP=DSP0K
	1009	IF(L<L0)THEN
	1010	DSP=((PK0-PK(L))DSPTK+(PK(L)-PKT)DSP0K)/(PK0-PKT)
	1011	ENDIF
	1012	ENDIF
	1013	!***
	1014	!*** HUMIDITY PROFILE
	1015	!***
	1016	IF(PK(L)>PQM)THEN
	1017	PSK(L)=PK(L)+DSP
	1018	APESK(L)=(1.E5/PSK(L))**CAPA
	1019	THSK(L)=TREFK(L)*APEK(L)
	1020	QREFK(L)=PQ0/PSK(L)EXP(A2(THSK(L)-A3*APESK(L)) &
	1021	& /(THSK(L)-A4*APESK(L)))
	1022	ELSE
	1023	QREFK(L)=QK(L)
	1024	ENDIF
	1025	!
	1026	ENDDO
	1027	!-----------------------------------------------------------------------
	1028	!***
	1029	!*** ENTHALPY CONSERVATION INTEGRAL
	1030	!***
	1031	!-----------------------------------------------------------------------
	1032	enthalpy_conservation : DO ITER=1,2
	1033	!
	1034	SUMDE=0.
	1035	SUMDP=0.
	1036	!
	1037	DO L=LTOP,LB
	1038	SUMDE=((TK(L)-TREFK(L))CP+(QK(L)-QREFK(L))EL(L))*DPRS(L) &
	1039	& +SUMDE
	1040	SUMDP=SUMDP+DPRS(L)
	1041	ENDDO
	1042	!
	1043	HCORR=SUMDE/(SUMDP-DPRS(LTOP))
	1044	LCOR=LTOP+1
	1045	!***
	1046	!*** FIND LQM
	1047	!***
	1048	LQM=1
	1049	DO L=KTS,LB
	1050	IF(PK(L)<=PQM)LQM=L
	1051	ENDDO
	1052	!***
	1053	!*** ABOVE LQM CORRECT TEMPERATURE ONLY
	1054	!***
	1055	IF(LCOR<=LQM)THEN
	1056	DO L=LCOR,LQM
	1057	TREFK(L)=TREFK(L)+HCORR*RCP
	1058	ENDDO
	1059	LCOR=LQM+1
	1060	ENDIF
	1061	!***
	1062	!*** BELOW LQM CORRECT BOTH TEMPERATURE AND MOISTURE
	1063	!***
	1064	DO L=LCOR,LB
	1065	TSKL=TREFK(L)*APEK(L)/APESK(L)
	1066	DHDT=QREFK(L)A23M4L/(TSKL-A4)*2+CP
	1067	TREFK(L)=HCORR/DHDT+TREFK(L)
	1068	THSKL=TREFK(L)*APEK(L)
	1069	QREFK(L)=PQ0/PSK(L)EXP(A2(THSKL-A3*APESK(L)) &
	1070	& /(THSKL-A4*APESK(L)))
	1071	ENDDO
	1072	!
	1073	ENDDO enthalpy_conservation
	1074	!-----------------------------------------------------------------------
	1075	!
	1076	!*** HEATING, MOISTENING, PRECIPITATION
	1077	!
	1078	!-----------------------------------------------------------------------
	1079	AVRGT=0.
	1080	PRECK=0.
	1081	DSQ=0.
	1082	DST=0.
	1083	!
	1084	DO L=LTOP,LB
	1085	TKL=TK(L)
	1086	DIFTL=(TREFK(L)-TKL )*TAUK
	1087	DIFQL=(QREFK(L)-QK(L))*TAUK
	1088	AVRGTL=(TKL+TKL+DIFTL)
	1089	DPOT=DPRS(L)/AVRGTL
	1090	DST=DIFTL*DPOT+DST
	1091	DSQ=DIFQLEL(L)DPOT+DSQ
	1092	AVRGT=AVRGTL*DPRS(L)+AVRGT
	1093	PRECK=DIFTL*DPRS(L)+PRECK
	1094	DIFT(L)=DIFTL
	1095	DIFQ(L)=DIFQL
	1096	ENDDO
	1097	!
	1098	DST=(DST+DST)*CP
	1099	DSQ=DSQ+DSQ
	1100	DENTPY=DST+DSQ
	1101	AVRGT=AVRGT/(SUMDP+SUMDP)
	1102	!
	1103	! DRHEAT=PRECK*CP/AVRGT
	1104	DRHEAT=(PRECKSM+MAX(1.E-7,PRECK)(1.-SM))*CP/AVRGT !As in Eta!
	1105	DRHEAT=MAX(DRHEAT,1.E-20)
	1106	EFI=EFIFC*DENTPY/DRHEAT
	1107	!-----------------------------------------------------------------------
	1108	EFI=MIN(EFI,1.)
	1109	EFI=MAX(EFI,EFIMN)
	1110	!-----------------------------------------------------------------------
	1111	!
	1112	ENDDO cloud_efficiency
	1113	!
	1114	!-----------------------------------------------------------------------
	1115	!
	1116	!-----------------------------------------------------------------------
	1117	!---------------------- DEEP CONVECTION --------------------------------
	1118	!-----------------------------------------------------------------------
	1119	!
	1120	IF(DENTPY>=EPSNTP.AND.PRECK>EPSPR)THEN
	1121	!
	1122	CLDEFI=EFI
	1123	FEFI=EFMNT+SLOPE*(EFI-EFIMN)
	1124	FEFI=(DENTPY-EPSNTP)*FEFI/DENTPY
	1125	PRECK=PRECK*FEFI
	1126	!
	1127	!*** UPDATE PRECIPITATION AND TENDENCIES OF TEMPERATURE AND MOISTURE
	1128	!
	1129	CUP=PRECK*CPRLG
	1130	PCPCOL=CUP
	1131	!
	1132	DO L=LTOP,LB
	1133	DTDT(L)=DIFT(L)FEFIRDTCNVC
	1134	DQDT(L)=DIFQ(L)FEFIRDTCNVC
	1135	ENDDO
	1136	!
	1137	ELSE
	1138	!
	1139	!-----------------------------------------------------------------------
	1140	!*** REDUCE THE CLOUD TOP
	1141	!-----------------------------------------------------------------------
	1142	!
	1143	! LTOP=LTOP+3
	1144	! PTOP=PRSMID(LTOP)
	1145	! DEPMIN=PSHPSFCRSFCP
	1146	! DEPTH=PBOT-PTOP
	1147	!***
	1148	!*** ITERATE DEEP CONVECTION PROCEDURE IF NEEDED
	1149	!***
	1150	! IF(DEPTH>=DEPMIN)THEN
	1151	! GO TO 300
	1152	! ENDIF
	1153	!
	1154	! CLDEFI=AVGEFI
	1155	CLDEFI=EFIMNSM+STEFI(1.-SM)
	1156	!***
	1157	!*** SEARCH FOR SHALLOW CLOUD TOP
	1158	!***
	1159	! LTSH=LBOT
	1160	! LBM1=LBOT-1
	1161	! PBTK=PK(LBOT)
	1162	! DEPMIN=PSHPSFCRSFCP
	1163	! PTPK=PBTK-DEPMIN
	1164	PTPK=MAX(PSHU, PK(LBOT)-DEPMIN)
	1165	!***
	1166	!*** CLOUD TOP IS THE LEVEL JUST BELOW PBTK-PSH or JUST BELOW PSHU
	1167	!***
	1168	DO L=KTS,LMH
	1169	IF(PK(L)<=PTPK)LTOP=L+1
	1170	ENDDO
	1171	!
	1172	! PTPK=PK(LTOP)
	1173	!!***
	1174	!!*** HIGHEST LEVEL ALLOWED IS LEVEL JUST BELOW PSHU
	1175	!!***
	1176	! IF(PTPK<=PSHU)THEN
	1177	!!
	1178	! DO L=KTS,LMH
	1179	! IF(PK(L)<=PSHU)LSHU=L+1
	1180	! ENDDO
	1181	!!
	1182	! LTOP=LSHU
	1183	! PTPK=PK(LTOP)
	1184	! ENDIF
	1185	!
	1186	! if(ltop>=lbot)then
	1187	!!!!!! lbot=0
	1188	! ltop=lmh
	1189	!!!!!! pbot=pk(lbot)
	1190	! ptop=pk(ltop)
	1191	! pbot=ptop
	1192	! go to 600
	1193	! endif
	1194	!
	1195	! LTP1=LTOP+1
	1196	! LTP2=LTOP+2
	1197	!!
	1198	! DO L=LTOP,LBOT
	1199	! QSATK(L)=PQ0/PK(L)EXP(A2(TK(L)-A3)/(TK(L)-A4))
	1200	! ENDDO
	1201	!!
	1202	! RHH=QK(LTOP)/QSATK(LTOP)
	1203	! RHMAX=0.
	1204	! LTSH=LTOP
	1205	!!
	1206	! DO L=LTP1,LBM1
	1207	! RHL=QK(L)/QSATK(L)
	1208	! DRHDP=(RHH-RHL)/(PK(L-1)-PK(L))
	1209	!!
	1210	! IF(DRHDP>RHMAX)THEN
	1211	! LTSH=L-1
	1212	! RHMAX=DRHDP
	1213	! ENDIF
	1214	!!
	1215	! RHH=RHL
	1216	! ENDDO
	1217	!
	1218	!-----------------------------------------------------------------------
	1219	!-- Make shallow cloud top a function of virtual temperature excess (DTV)
	1220	!-----------------------------------------------------------------------
	1221	!
	1222	LTP1=LBOT
	1223	DO L=LBOT-1,LTOP,-1
	1224	IF (DTV(L) > 0.) THEN
	1225	LTP1=L
	1226	ELSE
	1227	EXIT
	1228	ENDIF
	1229	ENDDO
	1230	LTOP=MIN(LTP1,LBOT)
	1231	!***
	1232	!*** CLOUD MUST BE AT LEAST TWO LAYERS THICK
	1233	!***
	1234	! IF(LBOT-LTOP<2)LTOP=LBOT-2 (eliminate this criterion)
	1235	!
	1236	!-- End: Buoyancy check (24 Aug 2006)
	1237	!
	1238	PTOP=PK(LTOP)
	1239	SHALLOW=.TRUE.
	1240	DEEP=.FALSE.
	1241	!
	1242	ENDIF
	1243	!DCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCD
	1244	!DCDCDCDCDCDCDC END OF DEEP CONVECTION DCDCDCDCDCDCD
	1245	!DCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCD
	1246	!
	1247	!-----------------------------------------------------------------------
	1248	!-----------------------------------------------------------------------
	1249	600 CONTINUE
	1250	!-----------------------------------------------------------------------
	1251	!-----------------------------------------------------------------------
	1252	!
	1253	!----------------GATHER SHALLOW CONVECTION POINTS-----------------------
	1254	!
	1255	! IF(PTOP<=PBOT-PNO.AND.LTOP<=LBOT-2)THEN
	1256	! DEPMIN=PSHPSFCRSFCP
	1257	!!
	1258	!! if(lpbl<lbot)lbot=lpbl
	1259	!! if(lbot>lmh-1)lbot=lmh-1
	1260	!! pbot=prsmid(lbot)
	1261	!!
	1262	! IF(PTOP+1.>=PBOT-DEPMIN)SHALLOW=.TRUE.
	1263	! ELSE
	1264	! LBOT=0
	1265	! LTOP=KTE
	1266	! ENDIF
	1267	!
	1268	!***********************************************************************
	1269	!-----------------------------------------------------------------------
	1270	!*** Begin debugging convection
	1271	IF(PRINT_DIAG)THEN
	1272	WRITE(6,"(a,2i3,L2,3e12.4)") &
	1273	'{cu2a lbot,ltop,shallow,pbot,ptop,depmin = ' &
	1274	,lbot,ltop,shallow,pbot,ptop,depmin
	1275	ENDIF
	1276	!*** End debugging convection
	1277	!-----------------------------------------------------------------------
	1278	!
	1279	IF(.NOT.SHALLOW)GO TO 800
	1280	!
	1281	!-----------------------------------------------------------------------
	1282	!***********************************************************************
	1283	!-----------------------------------------------------------------------
	1284	!SCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCS
	1285	!SCSCSCSCSCSCSC SHALLOW CONVECTION CSCSCSCSCSCSCSCSCSCS
	1286	!SCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCS
	1287	!-----------------------------------------------------------------------
	1288	DO L=KTS,LMH
	1289	TKL =T(L)
	1290	TK (L) =TKL
	1291	TREFK(L) =TKL
	1292	QKL =Q(L)
	1293	QK (L) =QKL
	1294	QREFK(L) =QKL
	1295	PKL =PRSMID(L)
	1296	PK (L) =PKL
	1297	QSATK(L) =PQ0/PK(L)EXP(A2(TK(L)-A3)/(TK(L)-A4))
	1298	APEKL =APE(L)
	1299	APEK (L) =APEKL
	1300	THVMKL =TKLAPEKL(QKL*D608+1.)
	1301	THVREF(L)=THVMKL
	1302	!
	1303	! IF(TKL>=TFRZ)THEN
	1304	EL(L)=ELWV
	1305	! ELSE
	1306	! EL(L)=ELIV
	1307	! ENDIF
	1308	ENDDO
	1309	!
	1310	!-----------------------------------------------------------------------
	1311	!-- Begin: Raise cloud top if avg RH>RHSHmax and CAPE>0
	1312	! RHSHmax=RH at cloud base associated with a DSP of PONE
	1313	!-----------------------------------------------------------------------
	1314	!
	1315	TLEV2=T(LBOT)((PK(LBOT)-PONE)/PK(LBOT))*CAPA
	1316	QSAT1=PQ0/PK(LBOT)EXP(A2(T(LBOT)-A3)/(TK(LBOT)-A4))
	1317	QSAT2=PQ0/(PK(LBOT)-PONE)EXP(A2(TLEV2-A3)/(TLEV2-A4))
	1318	RHSHmax=QSAT2/QSAT1
	1319	SUMDP=0.
	1320	RHAVG=0.
	1321	!
	1322	DO L=LBOT,LTOP,-1
	1323	RHAVG=RHAVG+DPRS(L)*QK(L)/QSATK(L)
	1324	SUMDP=SUMDP+DPRS(L)
	1325	ENDDO
	1326	!
	1327	IF (RHAVG/SUMDP > RHSHmax) THEN
	1328	LTSH=LTOP
	1329	DO L=LTOP-1,KTS,-1
	1330	RHAVG=RHAVG+DPRS(L)*QK(L)/QSATK(L)
	1331	SUMDP=SUMDP+DPRS(L)
	1332	IF (CPE(L) > 0.) THEN
	1333	LTSH=L
	1334	ELSE
	1335	EXIT
	1336	ENDIF
	1337	IF (RHAVG/SUMDP <= RHSHmax) EXIT
	1338	IF (PK(L) <= PSHU) EXIT
	1339	ENDDO
	1340	LTOP=LTSH
	1341	ENDIF
	1342	!
	1343	!-- End: Raise cloud top if avg RH>RHSHmax and CAPE>0
	1344	!
	1345	!---------------------------SHALLOW CLOUD TOP---------------------------
	1346	LBM1=LBOT-1
	1347	PTPK=PTOP
	1348	LTP1=LTOP-1
	1349	DEPTH=PBOT-PTOP
	1350	!-----------------------------------------------------------------------
	1351	!*** Begin debugging convection
	1352	IF(PRINT_DIAG)THEN
	1353	WRITE(6,"(a,4e12.4)") '{cu2b PBOT,PTOP,DEPTH,DEPMIN= ' &
	1354	,PBOT,PTOP,DEPTH,DEPMIN
	1355	ENDIF
	1356	!*** End debugging convection
	1357	!-----------------------------------------------------------------------
	1358	!
	1359	!BSF IF(DEPTH<DEPMIN)THEN
	1360	!BSF GO TO 800
	1361	!BSF ENDIF
	1362	!-----------------------------------------------------------------------
	1363	IF(PTOP>PBOT-PNO.OR.LTOP>LBOT-2)THEN
	1364	LBOT=0
	1365	!!! LTOP=LBOT
	1366	LTOP=KTE
	1367	PTOP=PBOT
	1368	GO TO 800
	1369	ENDIF
	1370	!
	1371	!--------------SCALING POTENTIAL TEMPERATURE & TABLE INDEX AT TOP-------
	1372	!
	1373	THTPK=T(LTP1)*APE(LTP1)
	1374	!
	1375	TTHK=(THTPK-THL)*RDTH
	1376	QQK =TTHK-AINT(TTHK)
	1377	IT =INT(TTHK)+1
	1378	!
	1379	IF(IT<1)THEN
	1380	IT=1
	1381	QQK=0.
	1382	ENDIF
	1383	!
	1384	IF(IT>=JTB)THEN
	1385	IT=JTB-1
	1386	QQK=0.
	1387	ENDIF
	1388	!
	1389	!--------------BASE AND SCALING FACTOR FOR SPEC. HUMIDITY AT TOP--------
	1390	!
	1391	BQS00K=QS0(IT)
	1392	SQS00K=SQS(IT)
	1393	BQS10K=QS0(IT+1)
	1394	SQS10K=SQS(IT+1)
	1395	!
	1396	!--------------SCALING SPEC. HUMIDITY & TABLE INDEX AT TOP--------------
	1397	!
	1398	BQK=(BQS10K-BQS00K)*QQK+BQS00K
	1399	SQK=(SQS10K-SQS00K)*QQK+SQS00K
	1400	!
	1401	! TQK=(Q(LTOP)-BQK)/SQK*RDQ
	1402	TQK=(Q(LTP1)-BQK)/SQK*RDQ
	1403	!
	1404	PPK=TQK-AINT(TQK)
	1405	IQ =INT(TQK)+1
	1406	!
	1407	IF(IQ<1)THEN
	1408	IQ=1
	1409	PPK=0.
	1410	ENDIF
	1411	!
	1412	IF(IQ>=ITB)THEN
	1413	IQ=ITB-1
	1414	PPK=0.
	1415	ENDIF
	1416	!
	1417	!----------------CLOUD TOP SATURATION POINT PRESSURE--------------------
	1418	!
	1419	PART1=(PTBL(IQ+1,IT)-PTBL(IQ,IT))*PPK
	1420	PART2=(PTBL(IQ,IT+1)-PTBL(IQ,IT))*QQK
	1421	PART3=(PTBL(IQ ,IT )-PTBL(IQ+1,IT ) &
	1422	& -PTBL(IQ ,IT+1)+PTBL(IQ+1,IT+1))PPKQQK
	1423	PTPK=PTBL(IQ,IT)+PART1+PART2+PART3
	1424	!-----------------------------------------------------------------------
	1425	DPMIX=PTPK-PSP
	1426	IF(ABS(DPMIX).LT.3000.)DPMIX=-3000.
	1427	!
	1428	!----------------TEMPERATURE PROFILE SLOPE------------------------------
	1429	!
	1430	SMIX=(THTPK-THBT)/DPMIX*STABS
	1431	!
	1432	TREFKX=TREFK(LBOT+1)
	1433	PKXXXX=PK(LBOT+1)
	1434	PKXXXY=PK(LBOT)
	1435	APEKXX=APEK(LBOT+1)
	1436	APEKXY=APEK(LBOT)
	1437	!
	1438	LMID=.5*(LBOT+LTOP)
	1439	!
	1440	DO L=LBOT,LTOP,-1
	1441	TREFKX=((PKXXXY-PKXXXX)*SMIX &
	1442	& +TREFKX*APEKXX)/APEKXY
	1443	TREFK(L)=TREFKX
	1444	IF(L<=LMID) TREFK(L)=MAX(TREFK(L), TK(L)+DTSHAL)
	1445	APEKXX=APEKXY
	1446	PKXXXX=PKXXXY
	1447	APEKXY=APEK(L-1)
	1448	PKXXXY=PK(L-1)
	1449	ENDDO
	1450	!
	1451	!----------------TEMPERATURE REFERENCE PROFILE CORRECTION---------------
	1452	!
	1453	SUMDT=0.
	1454	SUMDP=0.
	1455	!
	1456	DO L=LTOP,LBOT
	1457	SUMDT=(TK(L)-TREFK(L))*DPRS(L)+SUMDT
	1458	SUMDP=SUMDP+DPRS(L)
	1459	ENDDO
	1460	!
	1461	RDPSUM=1./SUMDP
	1462	FPK(LBOT)=TREFK(LBOT)
	1463	!
	1464	TCORR=SUMDT*RDPSUM
	1465	!
	1466	DO L=LTOP,LBOT
	1467	TRFKL =TREFK(L)+TCORR
	1468	TREFK(L)=TRFKL
	1469	FPK (L)=TRFKL
	1470	ENDDO
	1471	!
	1472	!----------------HUMIDITY PROFILE EQUATIONS-----------------------------
	1473	!
	1474	PSUM =0.
	1475	QSUM =0.
	1476	POTSUM=0.
	1477	QOTSUM=0.
	1478	OTSUM =0.
	1479	DST =0.
	1480	FPTK =FPK(LTOP)
	1481	!
	1482	DO L=LTOP,LBOT
	1483	DPKL =FPK(L)-FPTK
	1484	PSUM =DPKL *DPRS(L)+PSUM
	1485	QSUM =QK(L)*DPRS(L)+QSUM
	1486	RTBAR =2./(TREFK(L)+TK(L))
	1487	OTSUM =DPRS(L)*RTBAR+OTSUM
	1488	POTSUM=DPKL RTBARDPRS(L)+POTSUM
	1489	QOTSUM=QK(L) RTBARDPRS(L)+QOTSUM
	1490	DST =(TREFK(L)-TK(L))RTBARDPRS(L)/EL(L)+DST
	1491	ENDDO
	1492	!
	1493	PSUM =PSUM*RDPSUM
	1494	QSUM =QSUM*RDPSUM
	1495	ROTSUM=1./OTSUM
	1496	POTSUM=POTSUM*ROTSUM
	1497	QOTSUM=QOTSUM*ROTSUM
	1498	DST =DSTROTSUMCP
	1499	!
	1500	!-----------------------------------------------------------------------
	1501	!*** Begin debugging convection
	1502	IF(PRINT_DIAG)THEN
	1503	WRITE(6,"(a,5e12.4)") '{cu2c DST,PSUM,QSUM,POTSUM,QOTSUM = ' &
	1504	,DST,PSUM,QSUM,POTSUM,QOTSUM
	1505	ENDIF
	1506	!*** End debugging convection
	1507	!-----------------------------------------------------------------------
	1508	!
	1509	!----------------ENSURE POSITIVE ENTROPY CHANGE-------------------------
	1510	!
	1511	IF(DST>0.)THEN
	1512	! dstq=dst*epsup
	1513	LBOT=0
	1514	!!!! LTOP=LBOT
	1515	LTOP=KTE
	1516	PTOP=PBOT
	1517	GO TO 800
	1518	ELSE
	1519	DSTQ=DST*EPSDN
	1520	ENDIF
	1521	!
	1522	!----------------CHECK FOR ISOTHERMAL ATMOSPHERE------------------------
	1523	!
	1524	DEN=POTSUM-PSUM
	1525	!
	1526	IF(-DEN/PSUM<5.E-5)THEN
	1527	LBOT=0
	1528	!!!! LTOP=LBOT
	1529	LTOP=KTE
	1530	PTOP=PBOT
	1531	GO TO 800
	1532	!
	1533	!----------------SLOPE OF THE REFERENCE HUMIDITY PROFILE----------------
	1534	!
	1535	ELSE
	1536	DQREF=(QOTSUM-DSTQ-QSUM)/DEN
	1537	ENDIF
	1538	!
	1539	!-------------- HUMIDITY DOES NOT INCREASE WITH HEIGHT------------------
	1540	!
	1541	IF(DQREF<0.)THEN
	1542	LBOT=0
	1543	!!!! LTOP=LBOT
	1544	LTOP=KTE
	1545	PTOP=PBOT
	1546	GO TO 800
	1547	ENDIF
	1548	!
	1549	!----------------HUMIDITY AT THE CLOUD TOP------------------------------
	1550	!
	1551	QRFTP=QSUM-DQREF*PSUM
	1552	!
	1553	!----------------HUMIDITY PROFILE---------------------------------------
	1554	!
	1555	DO L=LTOP,LBOT
	1556	QRFKL=(FPK(L)-FPTK)*DQREF+QRFTP
	1557	!
	1558	!*** TOO DRY CLOUDS NOT ALLOWED
	1559	!
	1560	TNEW=(TREFK(L)-TK(L))*TAUKSC+TK(L)
	1561	QSATK(L)=PQ0/PK(L)EXP(A2(TNEW-A3)/(TNEW-A4))
	1562	QNEW=(QRFKL-QK(L))*TAUKSC+QK(L)
	1563	!
	1564	IF(QNEW<QSATK(L)*RHLSC)THEN
	1565	LBOT=0
	1566	!!!! LTOP=LBOT
	1567	LTOP=KTE
	1568	PTOP=PBOT
	1569	GO TO 800
	1570	ENDIF
	1571	!
	1572	!-------------TOO MOIST CLOUDS NOT ALLOWED------------------------------
	1573	!
	1574	IF(QNEW>QSATK(L)*RHHSC)THEN
	1575	LBOT=0
	1576	!!!! LTOP=LBOT
	1577	LTOP=KTE
	1578	PTOP=PBOT
	1579	GO TO 800
	1580	ENDIF
	1581
	1582	!
	1583	THVREF(L)=TREFK(L)APEK(L)(QRFKL*D608+1.)
	1584	QREFK(L)=QRFKL
	1585	ENDDO
	1586	!
	1587	!------------------ ELIMINATE CLOUDS WITH BOTTOMS TOO DRY --------------
	1588	!!
	1589	! qnew=(qrefk(lbot)-qk(lbot))*tauksc+qk(lbot)
	1590	!!
	1591	! if(qnew<qk(lbot+1)*stresh)then !!?? stresh too large!!
	1592	! lbot=0
	1593	!!!!!! ltop=lbot
	1594	! ltop=kte
	1595	! ptop=pbot
	1596	! go to 800
	1597	! endif
	1598	!!
	1599	!-------------- ELIMINATE IMPOSSIBLE SLOPES (BETTS,DTHETA/DQ)------------
	1600	!
	1601	DO L=LTOP,LBOT
	1602	DTDP=(THVREF(L-1)-THVREF(L))/(PRSMID(L)-PRSMID(L-1))
	1603	!
	1604	IF(DTDP<EPSDT)THEN
	1605	LBOT=0
	1606	!!!!! LTOP=LBOT
	1607	LTOP=KTE
	1608	PTOP=PBOT
	1609	GO TO 800
	1610	ENDIF
	1611	!
	1612	ENDDO
	1613	!-----------------------------------------------------------------------
	1614	!--------------RELAXATION TOWARD REFERENCE PROFILES---------------------
	1615	!-----------------------------------------------------------------------
	1616	!
	1617	DO L=LTOP,LBOT
	1618	DTDT(L)=(TREFK(L)-TK(L))TAUKSCRDTCNVC
	1619	DQDT(L)=(QREFK(L)-QK(L))TAUKSCRDTCNVC
	1620	ENDDO
	1621	!
	1622	!-----------------------------------------------------------------------
	1623	!*** Begin debugging convection
	1624	IF(PRINT_DIAG)THEN
	1625	DO L=LBOT,LTOP,-1
	1626	WRITE(6,"(a,i3,4e12.4)") '{cu2 KFLIP,DT,DTDT,DQ,DQDT = ' &
	1627	,KTE+1-L,TREFK(L)-TK(L),DTDT(L),QREFK(L)-QK(L),DQDT(L)
	1628	ENDDO
	1629	ENDIF
	1630	!*** End debugging convection
	1631	!-----------------------------------------------------------------------
	1632	!
	1633	!-----------------------------------------------------------------------
	1634	!SCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCS
	1635	!SCSCSCSCSCSCSC END OF SHALLOW CONVECTION SCSCSCSCSCSCSCS
	1636	!SCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCS
	1637	!-----------------------------------------------------------------------
	1638	800 CONTINUE
	1639	!-----------------------------------------------------------------------
	1640	END SUBROUTINE BMJ
	1641	!-----------------------------------------------------------------------
	1642	!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
	1643	!-----------------------------------------------------------------------
	1644	SUBROUTINE TTBLEX &
	1645	& (ITBX,JTBX,PLX,PRSMID,RDPX,RDTHEX,STHE &
	1646	& ,THE0,THESP,TTBL,TREF)
	1647	!-----------------------------------------------------------------------
	1648	! ******************************************************************
	1649	! * *
	1650	! * EXTRACT TEMPERATURE OF THE MOIST ADIABAT FROM *
	1651	! * THE APPROPRIATE TTBL *
	1652	! * *
	1653	! ******************************************************************
	1654	!-----------------------------------------------------------------------
	1655	IMPLICIT NONE
	1656	!-----------------------------------------------------------------------
	1657	INTEGER,INTENT(IN) :: ITBX,JTBX
	1658	!
	1659	REAL,INTENT(IN) :: PLX,PRSMID,RDPX,RDTHEX,THESP
	1660	!
	1661	REAL,DIMENSION(ITBX),INTENT(IN) :: STHE,THE0
	1662	!
	1663	REAL,DIMENSION(JTBX,ITBX),INTENT(IN) :: TTBL
	1664	!
	1665	REAL,INTENT(OUT) :: TREF
	1666	!-----------------------------------------------------------------------
	1667	REAL :: BTHE00K,BTHE10K,BTHK,PK,PP,QQ,STHE00K,STHE10K,STHK &
	1668	& ,T00K,T01K,T10K,T11K,TPK,TTHK
	1669	!
	1670	INTEGER :: IPTB,ITHTB
	1671	!-----------------------------------------------------------------------
	1672	!----------------SCALING PRESSURE & TT TABLE INDEX----------------------
	1673	!-----------------------------------------------------------------------
	1674	PK=PRSMID
	1675	TPK=(PK-PLX)*RDPX
	1676	QQ=TPK-AINT(TPK)
	1677	IPTB=INT(TPK)+1
	1678	!----------------KEEPING INDICES WITHIN THE TABLE-----------------------
	1679	IF(IPTB<1)THEN
	1680	IPTB=1
	1681	QQ=0.
	1682	ENDIF
	1683	!
	1684	IF(IPTB>=ITBX)THEN
	1685	IPTB=ITBX-1
	1686	QQ=0.
	1687	ENDIF
	1688	!----------------BASE AND SCALING FACTOR FOR THETAE---------------------
	1689	BTHE00K=THE0(IPTB)
	1690	STHE00K=STHE(IPTB)
	1691	BTHE10K=THE0(IPTB+1)
	1692	STHE10K=STHE(IPTB+1)
	1693	!----------------SCALING THE & TT TABLE INDEX---------------------------
	1694	BTHK=(BTHE10K-BTHE00K)*QQ+BTHE00K
	1695	STHK=(STHE10K-STHE00K)*QQ+STHE00K
	1696	TTHK=(THESP-BTHK)/STHK*RDTHEX
	1697	PP=TTHK-AINT(TTHK)
	1698	ITHTB=INT(TTHK)+1
	1699	!----------------KEEPING INDICES WITHIN THE TABLE-----------------------
	1700	IF(ITHTB<1)THEN
	1701	ITHTB=1
	1702	PP=0.
	1703	ENDIF
	1704	!
	1705	IF(ITHTB>=JTBX)THEN
	1706	ITHTB=JTBX-1
	1707	PP=0.
	1708	ENDIF
	1709	!----------------TEMPERATURE AT FOUR SURROUNDING TT TABLE PTS.----------
	1710	T00K=TTBL(ITHTB,IPTB)
	1711	T10K=TTBL(ITHTB+1,IPTB)
	1712	T01K=TTBL(ITHTB,IPTB+1)
	1713	T11K=TTBL(ITHTB+1,IPTB+1)
	1714	!-----------------------------------------------------------------------
	1715	!----------------PARCEL TEMPERATURE-------------------------------------
	1716	!-----------------------------------------------------------------------
	1717	TREF=(T00K+(T10K-T00K)PP+(T01K-T00K)QQ &
	1718	& +(T00K-T10K-T01K+T11K)PPQQ)
	1719	!-----------------------------------------------------------------------
	1720	END SUBROUTINE TTBLEX
	1721	!-----------------------------------------------------------------------
	1722	!-----------------------------------------------------------------------
	1723	SUBROUTINE BMJINIT(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN &
	1724	& ,CLDEFI,LOWLYR,CP,RD,RESTART &
	1725	& ,ALLOWED_TO_READ &
	1726	& ,IDS,IDE,JDS,JDE,KDS,KDE &
	1727	& ,IMS,IME,JMS,JME,KMS,KME &
	1728	& ,ITS,ITE,JTS,JTE,KTS,KTE)
	1729	!-----------------------------------------------------------------------
	1730	IMPLICIT NONE
	1731	!-----------------------------------------------------------------------
	1732	LOGICAL,INTENT(IN) :: RESTART,ALLOWED_TO_READ
	1733	INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
	1734	& ,IMS,IME,JMS,JME,KMS,KME &
	1735	& ,ITS,ITE,JTS,JTE,KTS,KTE
	1736	!
	1737	REAL,INTENT(IN) :: CP,RD
	1738	!
	1739	REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: &
	1740	& RTHCUTEN &
	1741	& ,RQVCUTEN &
	1742	& ,RQCCUTEN &
	1743	& ,RQRCUTEN
	1744	!
	1745	REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: CLDEFI
	1746
	1747	INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: LOWLYR
	1748	!
	1749	REAL,PARAMETER :: EPS=1.E-9
	1750	!
	1751	REAL, DIMENSION(JTB) :: APP,APT,AQP,AQT,PNEW,POLD,QSNEW,QSOLD &
	1752	& ,THENEW,THEOLD,TNEW,TOLD,Y2P,Y2T
	1753	!
	1754	REAL,DIMENSION(JTBQ) :: APTQ,AQTQ,THENEWQ,THEOLDQ &
	1755	& ,TNEWQ,TOLDQ,Y2TQ
	1756	!
	1757	INTEGER :: I,J,K,ITF,JTF,KTF
	1758	INTEGER :: KTH,KTHM,KTHM1,KP,KPM,KPM1
	1759	!
	1760	REAL :: APE,DP,DQS,DTH,DTHE,P,QS,QS0K,SQSK,STHEK &
	1761	& ,TH,THE0K,DENOM,ELOCP
	1762	!-----------------------------------------------------------------------
	1763
	1764	ELOCP=ELIWV/CP
	1765	JTF=MIN0(JTE,JDE-1)
	1766	KTF=MIN0(KTE,KDE-1)
	1767	ITF=MIN0(ITE,IDE-1)
	1768	!
	1769	IF(.NOT.RESTART)THEN
	1770	DO J=JTS,JTF
	1771	DO K=KTS,KTF
	1772	DO I=ITS,ITF
	1773	RTHCUTEN(I,K,J)=0.
	1774	RQVCUTEN(I,K,J)=0.
	1775	RQCCUTEN(I,K,J)=0.
	1776	RQRCUTEN(I,K,J)=0.
	1777	ENDDO
	1778	ENDDO
	1779	ENDDO
	1780	!
	1781	DO J=JTS,JTF
	1782	DO I=ITS,ITF
	1783	CLDEFI(I,J)=AVGEFI
	1784	ENDDO
	1785	ENDDO
	1786	ENDIF
	1787	!
	1788	!*** FOR NOW, ASSUME SIGMA MODE FOR LOWEST MODEL LAYER
	1789	!
	1790	DO J=JTS,JTF
	1791	DO I=ITS,ITF
	1792	LOWLYR(I,J)=1
	1793	ENDDO
	1794	ENDDO
	1795	!-----------------------------------------------------------------------
	1796	!----------------COARSE LOOK-UP TABLE FOR SATURATION POINT--------------
	1797	!-----------------------------------------------------------------------
	1798	!
	1799	KTHM=JTB
	1800	KPM=ITB
	1801	KTHM1=KTHM-1
	1802	KPM1=KPM-1
	1803	!
	1804	DTH=(THH-THL)/REAL(KTHM-1)
	1805	DP =(PH -PL )/REAL(KPM -1)
	1806	!
	1807	TH=THL-DTH
	1808	!-----------------------------------------------------------------------
	1809	DO 100 KTH=1,KTHM
	1810	!
	1811	TH=TH+DTH
	1812	P=PL-DP
	1813	!
	1814	DO KP=1,KPM
	1815	P=P+DP
	1816	APE=(100000./P)**(RD/CP)
	1817	DENOM=TH-A4*APE
	1818	IF (DENOM>EPS) THEN
	1819	QSOLD(KP)=PQ0/PEXP(A2(TH-A3*APE)/DENOM)
	1820	ELSE
	1821	QSOLD(KP)=0.
	1822	ENDIF
	1823	POLD(KP)=P
	1824	ENDDO
	1825	!
	1826	QS0K=QSOLD(1)
	1827	SQSK=QSOLD(KPM)-QSOLD(1)
	1828	QSOLD(1 )=0.
	1829	QSOLD(KPM)=1.
	1830	!
	1831	DO KP=2,KPM1
	1832	QSOLD(KP)=(QSOLD(KP)-QS0K)/SQSK
	1833	IF((QSOLD(KP)-QSOLD(KP-1)).LT.EPS)QSOLD(KP)=QSOLD(KP-1)+EPS
	1834	ENDDO
	1835	!
	1836	QS0(KTH)=QS0K
	1837	QS0_EXP(KTH)=QS0K
	1838	SQS(KTH)=SQSK
	1839	SQS_EXP(KTH)=SQSK
	1840	!-----------------------------------------------------------------------
	1841	QSNEW(1 )=0.
	1842	QSNEW(KPM)=1.
	1843	DQS=1./REAL(KPM-1)
	1844	!
	1845	DO KP=2,KPM1
	1846	QSNEW(KP)=QSNEW(KP-1)+DQS
	1847	ENDDO
	1848	!
	1849	Y2P(1 )=0.
	1850	Y2P(KPM )=0.
	1851	!
	1852	CALL SPLINE(JTB,KPM,QSOLD,POLD,Y2P,KPM,QSNEW,PNEW,APP,AQP)
	1853	!
	1854	DO KP=1,KPM
	1855	PTBL(KP,KTH)=PNEW(KP)
	1856	PTBL_EXP(KP,KTH)=PNEW(KP)
	1857	ENDDO
	1858	!-----------------------------------------------------------------------
	1859	100 CONTINUE
	1860	!-----------------------------------------------------------------------
	1861	!------------COARSE LOOK-UP TABLE FOR T(P) FROM CONSTANT THE------------
	1862	!-----------------------------------------------------------------------
	1863	P=PL-DP
	1864	!
	1865	DO 200 KP=1,KPM
	1866	!
	1867	P=P+DP
	1868	TH=THL-DTH
	1869	!
	1870	DO KTH=1,KTHM
	1871	TH=TH+DTH
	1872	APE=(1.E5/P)**(RD/CP)
	1873	DENOM=TH-A4*APE
	1874	IF (DENOM>EPS) THEN
	1875	QS=PQ0/PEXP(A2(TH-A3*APE)/DENOM)
	1876	ELSE
	1877	QS=0.
	1878	ENDIF
	1879	! QS=PQ0/PEXP(A2(TH-A3APE)/(TH-A4APE))
	1880	TOLD(KTH)=TH/APE
	1881	THEOLD(KTH)=THEXP(ELOCPQS/TOLD(KTH))
	1882	ENDDO
	1883	!
	1884	THE0K=THEOLD(1)
	1885	STHEK=THEOLD(KTHM)-THEOLD(1)
	1886	THEOLD(1 )=0.
	1887	THEOLD(KTHM)=1.
	1888	!
	1889	DO KTH=2,KTHM1
	1890	THEOLD(KTH)=(THEOLD(KTH)-THE0K)/STHEK
	1891	IF((THEOLD(KTH)-THEOLD(KTH-1)).LT.EPS) &
	1892	& THEOLD(KTH)=THEOLD(KTH-1) + EPS
	1893	ENDDO
	1894	!
	1895	THE0(KP)=THE0K
	1896	STHE(KP)=STHEK
	1897	!-----------------------------------------------------------------------
	1898	THENEW(1 )=0.
	1899	THENEW(KTHM)=1.
	1900	DTHE=1./REAL(KTHM-1)
	1901	!
	1902	DO KTH=2,KTHM1
	1903	THENEW(KTH)=THENEW(KTH-1)+DTHE
	1904	ENDDO
	1905	!
	1906	Y2T(1 )=0.
	1907	Y2T(KTHM)=0.
	1908	!
	1909	CALL SPLINE(JTB,KTHM,THEOLD,TOLD,Y2T,KTHM,THENEW,TNEW,APT,AQT)
	1910	!
	1911	DO KTH=1,KTHM
	1912	TTBL(KTH,KP)=TNEW(KTH)
	1913	ENDDO
	1914	!-----------------------------------------------------------------------
	1915	200 CONTINUE
	1916	!-----------------------------------------------------------------------
	1917	!
	1918	!-----------------------------------------------------------------------
	1919	!---------------FINE LOOK-UP TABLE FOR SATURATION POINT-----------------
	1920	!-----------------------------------------------------------------------
	1921	KTHM=JTBQ
	1922	KPM=ITBQ
	1923	KTHM1=KTHM-1
	1924	KPM1=KPM-1
	1925	!
	1926	DTH=(THHQ-THL)/REAL(KTHM-1)
	1927	DP=(PH-PLQ)/REAL(KPM-1)
	1928	!
	1929	TH=THL-DTH
	1930	P=PLQ-DP
	1931	!-----------------------------------------------------------------------
	1932	!---------------FINE LOOK-UP TABLE FOR T(P) FROM CONSTANT THE-----------
	1933	!-----------------------------------------------------------------------
	1934	DO 300 KP=1,KPM
	1935	!
	1936	P=P+DP
	1937	TH=THL-DTH
	1938	!
	1939	DO KTH=1,KTHM
	1940	TH=TH+DTH
	1941	APE=(1.E5/P)**(RD/CP)
	1942	DENOM=TH-A4*APE
	1943	IF (DENOM>EPS) THEN
	1944	QS=PQ0/PEXP(A2(TH-A3*APE)/DENOM)
	1945	ELSE
	1946	QS=0.
	1947	ENDIF
	1948	! QS=PQ0/PEXP(A2(TH-A3APE)/(TH-A4APE))
	1949	TOLDQ(KTH)=TH/APE
	1950	THEOLDQ(KTH)=THEXP(ELOCPQS/TOLDQ(KTH))
	1951	ENDDO
	1952	!
	1953	THE0K=THEOLDQ(1)
	1954	STHEK=THEOLDQ(KTHM)-THEOLDQ(1)
	1955	THEOLDQ(1 )=0.
	1956	THEOLDQ(KTHM)=1.
	1957	!
	1958	DO KTH=2,KTHM1
	1959	THEOLDQ(KTH)=(THEOLDQ(KTH)-THE0K)/STHEK
	1960	IF((THEOLDQ(KTH)-THEOLDQ(KTH-1))<EPS) &
	1961	& THEOLDQ(KTH)=THEOLDQ(KTH-1)+EPS
	1962	ENDDO
	1963	!
	1964	THE0Q(KP)=THE0K
	1965	STHEQ(KP)=STHEK
	1966	!-----------------------------------------------------------------------
	1967	THENEWQ(1 )=0.
	1968	THENEWQ(KTHM)=1.
	1969	DTHE=1./REAL(KTHM-1)
	1970	!
	1971	DO KTH=2,KTHM1
	1972	THENEWQ(KTH)=THENEWQ(KTH-1)+DTHE
	1973	ENDDO
	1974	!
	1975	Y2TQ(1 )=0.
	1976	Y2TQ(KTHM)=0.
	1977	!
	1978	CALL SPLINE(JTBQ,KTHM,THEOLDQ,TOLDQ,Y2TQ,KTHM &
	1979	& ,THENEWQ,TNEWQ,APTQ,AQTQ)
	1980	!
	1981	DO KTH=1,KTHM
	1982	TTBLQ(KTH,KP)=TNEWQ(KTH)
	1983	ENDDO
	1984	!-----------------------------------------------------------------------
	1985	300 CONTINUE
	1986	!-----------------------------------------------------------------------
	1987	END SUBROUTINE BMJINIT
	1988	!-----------------------------------------------------------------------
	1989	!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
	1990	!-----------------------------------------------------------------------
	1991	SUBROUTINE SPLINE(JTBX,NOLD,XOLD,YOLD,Y2,NNEW,XNEW,YNEW,P,Q)
	1992	! ********************************************************************
	1993	! * *
	1994	! * THIS IS A ONE-DIMENSIONAL CUBIC SPLINE FITTING ROUTINE *
	1995	! * PROGRAMED FOR A SMALL SCALAR MACHINE. *
	1996	! * *
	1997	! * PROGRAMER Z. JANJIC *
	1998	! * *
	1999	! * NOLD - NUMBER OF GIVEN VALUES OF THE FUNCTION. MUST BE GE 3. *
	2000	! * XOLD - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE *
	2001	! * FUNCTION ARE GIVEN. MUST BE IN ASCENDING ORDER. *
	2002	! * YOLD - THE GIVEN VALUES OF THE FUNCTION AT THE POINTS XOLD. *
	2003	! * Y2 - THE SECOND DERIVATIVES AT THE POINTS XOLD. IF NATURAL *
	2004	! * SPLINE IS FITTED Y2(1)=0. AND Y2(NOLD)=0. MUST BE *
	2005	! * SPECIFIED. *
	2006	! * NNEW - NUMBER OF VALUES OF THE FUNCTION TO BE CALCULATED. *
	2007	! * XNEW - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE *
	2008	! * FUNCTION ARE CALCULATED. XNEW(K) MUST BE GE XOLD(1) *
	2009	! * AND LE XOLD(NOLD). *
	2010	! * YNEW - THE VALUES OF THE FUNCTION TO BE CALCULATED. *
	2011	! * P, Q - AUXILIARY VECTORS OF THE LENGTH NOLD-2. *
	2012	! * *
	2013	! ********************************************************************
	2014	!-----------------------------------------------------------------------
	2015	IMPLICIT NONE
	2016	!-----------------------------------------------------------------------
	2017	INTEGER,INTENT(IN) :: JTBX,NNEW,NOLD
	2018	REAL,DIMENSION(JTBX),INTENT(IN) :: XNEW,XOLD,YOLD
	2019	REAL,DIMENSION(JTBX),INTENT(INOUT) :: P,Q,Y2
	2020	REAL,DIMENSION(JTBX),INTENT(OUT) :: YNEW
	2021	!
	2022	INTEGER :: K,K1,K2,KOLD,NOLDM1
	2023	REAL :: AK,BK,CK,DEN,DX,DXC,DXL,DXR,DYDXL,DYDXR &
	2024	& ,RDX,RTDXC,X,XK,XSQ,Y2K,Y2KP1
	2025	!-----------------------------------------------------------------------
	2026	NOLDM1=NOLD-1
	2027	!
	2028	DXL=XOLD(2)-XOLD(1)
	2029	DXR=XOLD(3)-XOLD(2)
	2030	DYDXL=(YOLD(2)-YOLD(1))/DXL
	2031	DYDXR=(YOLD(3)-YOLD(2))/DXR
	2032	RTDXC=0.5/(DXL+DXR)
	2033	!
	2034	P(1)= RTDXC(6.(DYDXR-DYDXL)-DXL*Y2(1))
	2035	Q(1)=-RTDXC*DXR
	2036	!
	2037	IF(NOLD==3)GO TO 150
	2038	!-----------------------------------------------------------------------
	2039	K=3
	2040	!
	2041	100 DXL=DXR
	2042	DYDXL=DYDXR
	2043	DXR=XOLD(K+1)-XOLD(K)
	2044	DYDXR=(YOLD(K+1)-YOLD(K))/DXR
	2045	DXC=DXL+DXR
	2046	DEN=1./(DXL*Q(K-2)+DXC+DXC)
	2047	!
	2048	P(K-1)= DEN(6.(DYDXR-DYDXL)-DXL*P(K-2))
	2049	Q(K-1)=-DEN*DXR
	2050	!
	2051	K=K+1
	2052	IF(K<NOLD)GO TO 100
	2053	!-----------------------------------------------------------------------
	2054	150 K=NOLDM1
	2055	!
	2056	200 Y2(K)=P(K-1)+Q(K-1)*Y2(K+1)
	2057	!
	2058	K=K-1
	2059	IF(K>1)GO TO 200
	2060	!-----------------------------------------------------------------------
	2061	K1=1
	2062	!
	2063	300 XK=XNEW(K1)
	2064	!
	2065	DO 400 K2=2,NOLD
	2066	!
	2067	IF(XOLD(K2)>XK)THEN
	2068	KOLD=K2-1
	2069	GO TO 450
	2070	ENDIF
	2071	!
	2072	400 CONTINUE
	2073	!
	2074	YNEW(K1)=YOLD(NOLD)
	2075	GO TO 600
	2076	!
	2077	450 IF(K1==1)GO TO 500
	2078	IF(K==KOLD)GO TO 550
	2079	!
	2080	500 K=KOLD
	2081	!
	2082	Y2K=Y2(K)
	2083	Y2KP1=Y2(K+1)
	2084	DX=XOLD(K+1)-XOLD(K)
	2085	RDX=1./DX
	2086	!
	2087	AK=.1666667RDX(Y2KP1-Y2K)
	2088	BK=0.5*Y2K
	2089	CK=RDX(YOLD(K+1)-YOLD(K))-.1666667DX*(Y2KP1+Y2K+Y2K)
	2090	!
	2091	550 X=XK-XOLD(K)
	2092	XSQ=X*X
	2093	!
	2094	YNEW(K1)=AKXSQX+BKXSQ+CKX+YOLD(K)
	2095	!
	2096	600 K1=K1+1
	2097	IF(K1<=NNEW)GO TO 300
	2098	!-----------------------------------------------------------------------
	2099	END SUBROUTINE SPLINE
	2100	!-----------------------------------------------------------------------
	2101	!
	2102	END MODULE MODULE_CU_BMJ
	2103	!
	2104	!-----------------------------------------------------------------------

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

Download in other formats:

Original Format