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