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