module_cu_bmj.F
References to this file elsewhere.
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)*G*DZ8W(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)*DTCNVC*ABS(DQDT(K))
254 DELT=DELT+DPCOL(K)*DTCNVC*ABS(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,.01*PLYR &
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,DTCNVC*DTDT(K),1000.*DTCNVC*DQDT(K),.01*DPCOL(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=DSPBFL*FSL,DSP0SL=DSP0FL*FSL &
366 & ,DSPTSL=DSPTFL*FSL &
367 & ,DSPBSS=DSPBFS*FSS,DSP0SS=DSP0FS*FSS &
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/(ROW*G*ELWV)
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)*PP1*QQ1
482 APESP=(1.E5/TPSP)**CAPA
483 TTHES=TTHBT*EXP(ELOCP*QBT*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=FFUP*FPRS*FPRS*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)=((-FFUP*DPLO+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*(QQBT*0.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*(QQBT*0.608+1.) &
623 ! & -TSP*(QSP*0.608+1.))*0.5 &
624 ! & /(TSP*(QSP*0.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/PUP*EXP(A2*(TUP-A3)/(TUP-A4))
632 ! DP=PLO-PUP
633 ! TRMUP=(TUP*(QUP*0.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/PUP*EXP(A2*(TUP-A3)/(TUP-A4))
645 ! DP=PLO-PUP
646 ! TRMUP=(TUP*(QUP*0.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=AVGEFI*SM+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=PSH*PSFC*RSFCP
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)*DTHEM*RDP0T)/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=PFRZ*PSFC*RSFCP
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+DSPBSS*PBOTFC)*SM &
807 & +((EFI-EFIMN)*SLOPBL+DSPBSL*PBOTFC)*SM1
808 DSP0K=((EFI-EFIMN)*SLOP0S+DSP0SS*PBOTFC)*SM &
809 & +((EFI-EFIMN)*SLOP0L+DSP0SL*PBOTFC)*SM1
810 DSPTK=((EFI-EFIMN)*SLOPTS+DSPTSS*PBOTFC)*SM &
811 & +((EFI-EFIMN)*SLOPTL+DSPTSL*PBOTFC)*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=DIFQL*EL(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=(PRECK*SM+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)*FEFI*RDTCNVC
950 DQDT(L)=DIFQ(L)*FEFI*RDTCNVC
951 ENDDO
952 !
953 ELSE
954 !
955 !-----------------------------------------------------------------------
956 !*** REDUCE THE CLOUD TOP
957 !-----------------------------------------------------------------------
958 !
959 ! LTOP=LTOP+3
960 ! PTOP=PRSMID(LTOP)
961 ! DEPMIN=PSH*PSFC*RSFCP
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=EFIMN*SM+STEFI*(1.-SM)
972 !***
973 !*** SEARCH FOR SHALLOW CLOUD TOP
974 !***
975 LTSH=LBOT
976 LBM1=LBOT-1
977 PBTK=PK(LBOT)
978 DEPMIN=PSH*PSFC*RSFCP
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=PSH*PSFC*RSFCP
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 =TKL*APEKL*(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))*PPK*QQK
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 *RTBAR*DPRS(L)+POTSUM
1235 QOTSUM=QK(L) *RTBAR*DPRS(L)+QOTSUM
1236 DST =(TREFK(L)-TK(L))*RTBAR*DPRS(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 =DST*ROTSUM*CP
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))*TAUKSC*RDTCNVC
1356 DQDT(L)=(QREFK(L)-QK(L))*TAUKSC*RDTCNVC
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)*PP*QQ)
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/P*EXP(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/P*EXP(A2*(TH-A3*APE)/DENOM)
1612 ELSE
1613 QS=0.
1614 ENDIF
1615 ! QS=PQ0/P*EXP(A2*(TH-A3*APE)/(TH-A4*APE))
1616 TOLD(KTH)=TH/APE
1617 THEOLD(KTH)=TH*EXP(ELIWV*QS/(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/P*EXP(A2*(TH-A3*APE)/DENOM)
1681 ELSE
1682 QS=0.
1683 ENDIF
1684 ! QS=PQ0/P*EXP(A2*(TH-A3*APE)/(TH-A4*APE))
1685 TOLDQ(KTH)=TH/APE
1686 THEOLDQ(KTH)=TH*EXP(ELIWV*QS/(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=.1666667*RDX*(Y2KP1-Y2K)
1824 BK=0.5*Y2K
1825 CK=RDX*(YOLD(K+1)-YOLD(K))-.1666667*DX*(Y2KP1+Y2K+Y2K)
1826 !
1827 550 X=XK-XOLD(K)
1828 XSQ=X*X
1829 !
1830 YNEW(K1)=AK*XSQ*X+BK*XSQ+CK*X+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 !-----------------------------------------------------------------------