Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_cu_bmj.F @ 3094

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

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

Download in other formats:

Original Format