phys/module_cu_bmj.F

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