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