Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_cu_bmj.F @ 14

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

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

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 67.4 KB

Line
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