!########################################################################
!########################################################################
!#########                                                      #########
!#########     WRF:MODEL_LAYER:PHYSICS:MODULE module_cu_bmj     #########
!#########                                                      #########
!#########                      Adapted by                      #########
!#########     Center for Analysis and Prediction of Storms     #########
!#########                University of Oklahoma                #########
!#########                                                      #########
!########################################################################
!########################################################################


MODULE module_cu_bmj 2

!------------------------------------------------------------------------
!
! PURPOSE:
!
! Contains the routines and parameters for the Betts-Miller-Janjic 
! deep and shallow convective adjustment scheme.
!
! REFERENCES:
!
! Janjic, Z. I., 1994:  The step-mountain Eta coordinate model:  Further
!   developments of the convection, viscous sublayer, and turbulence
!   closure schemes.  _Mon. Wea. Rev._, 122, 927-945.
! Betts, A. K., and M. J. Miller, 1993:  The Betts-Miller Scheme.  _The
!   Representation of Cumulus Convection in Numerical Models_.  _Meteor.
!   Mono._, 24, 107-121.
! Janjic, Z. I., 1990:  The step-mountain coordinate:  Physical package.
!   _Mon. Wea. Rev._, 118, 1429-1443.
! Betts, A. K., 1986:  A new convective adjustment scheme.  Part I:  
!   Observational and theroetical basis.  _Quart. J. Roy. Meteor. Soc._,
!   112, 677-691.
! Betts, A. K., and M. J. Miller, 1986:  A new convective adjustment
!   scheme.  Part II:  Single column tests using GATE wave, BOMEX, ATEX
!   and arctic air-mass data sets.  _Quart. J. Roy. Meteor. Soc._, 112,
!   693-709.
! Betts, A. K., 1985:  Mixing line analysis of clouds and cloudy 
!   boundary layers.  _J. Atmos. Sci._, 42, 2751-2763.
! Betts, A. K., 1984:  Boundary layer thermodynamics of a High Plains
!   severe storm.  _Mon. Wea. Rev._, 112, 2199-2211.
! Betts, A. K., 1983a:  Thermodynamics of mixed stratocumulus layers:  
!   Saturation point budgets.  _J. Atmos. Sci._, 40, 2655-2670.
! Betts, A. K., 1983b:  Atmospheric convective structure and a convection
!   scheme based on saturation point adjustment.  Workshop on convection
!   in large-scale models, 20 Nov. to 1 Dec., ECMWF.
! Betts, A. K., 1982a:  Saturation point analysis of moist convective
!   overturning.  _J. Atmos. Sci._, 39, 1484-1505.
! Betts, A. K., 1982b:  Cloud thermodynamic models in saturation point
!   coordinates.  _J. Atmos. Sci._, 39, 2182-2191.
!
!------------------------------------------------------------------------
!
! AUTHOR: ???
!
! MODIFICATIONS:
!
! Eric Kemp, 21 September 2001
! Reformatted code.
!
!------------------------------------------------------------------------
!
! Declare module parameters and variables.
!
!------------------------------------------------------------------------

  LOGICAL :: UNIS=.TRUE.,UNIL=.FALSE.

  REAL,PARAMETER :: A2=17.2693882,A3=273.16,A4=35.86                     &
                       ,DSPC=-3000.                                      &
                       ,DTTOP=0.,EFIFC=5.0,EFIMN=.20,EFMNT=.70           &
                       ,ELIVW=2.72E6                                     &
                       ,EPSDN=1.05,EPSDT=0.,EPSNTP=1.E2                  &
                       ,EPSP=1.E-7,EPSQ=2.E-12,EPSUP=1.0                 &
                       ,FCC=.50,FSL=.850,FSS=.85                         &
                       ,PBM=13000.,PFRZ=15000.,PNO=1000.                 &
                       ,PONE=2500.,PQM=20000.,PQ0=379.90516              &
                       ,PSH=20000.,PSHU=45000.,RHF=0.10,ROW=1.E3         &
                       ,STABD=.90,STABFC=1.0                             &
                       ,STABS=1.0,STRESH=1.10                            &
                       ,T1=274.16,TREL=2400.

  REAL,PARAMETER :: DSPBFL=-4843.75,DSP0FL=-7050.0,DSPTFL=-2250.0        &
                       ,DSPBFS=-3875.,DSP0FS=-5875.,DSPTFS=-1875.

  REAL,PARAMETER :: PL=2500.,PLQ=70000.,PH=105000.                       &
                       ,THL=210.,THH=365.,THHQ=325.

  INTEGER,PARAMETER :: ITB=76,JTB=134,ITBQ=152,JTBQ=440

!------------------------------------------------------------------------
!
! Arrays for lookup tables.
!
!------------------------------------------------------------------------

  REAL,DIMENSION(ITB),PRIVATE,SAVE :: STHE,THE0
  REAL,DIMENSION(JTB),PRIVATE,SAVE :: QS0,SQS
  REAL,DIMENSION(ITBQ),PRIVATE,SAVE :: STHEQ,THE0Q
  REAL,DIMENSION(ITB,JTB),PRIVATE,SAVE :: PTBL
  REAL,DIMENSION(JTB,ITB),PRIVATE,SAVE :: TTBL
  REAL,DIMENSION(JTBQ,ITBQ),PRIVATE,SAVE :: TTBLQ

  REAL,PARAMETER :: RDP=(ITB-1.)/(PH-PL),RDPQ=(ITBQ-1.)/(PH-PLQ)        &
                       ,RDQ=ITB-1,RDTH=(JTB-1.)/(THH-THL)               &
                       ,RDTHE=JTB-1.,RDTHEQ=JTBQ-1.

CONTAINS ! Subroutine declarations

!########################################################################
!########################################################################
!#########                                                      #########
!#########                  SUBROUTINE bmjdrv                   #########
!#########                                                      #########
!#########                      Adapted by                      #########
!#########     Center for Analysis and Prediction of Storms     #########
!#########                University of Oklahoma                #########
!#########                                                      #########
!########################################################################
!########################################################################


SUBROUTINE BMJDRV(DT,ITIMESTEP,STEPCU                                   & 1,1
                       ,RR,RTHCUTEN,RQVCUTEN                            &
                       ,RAINCV                                          &
                       ,TH,T,QV,PINT,PMID,PI,RHO,DZ8W                   &
                       ,CP,R,RVOVRD,ELWV,G,TFRZ,D608                    &
                       ,CLDEFI,LOWLYR,XLAND                             &
                       ,IDS,IDE,JDS,JDE,KDS,KDE                         &
                       ,IMS,IME,JMS,JME,KMS,KME                         &
                       ,ITS,ITE,JTS,JTE,KTS,KTE)

!------------------------------------------------------------------------
!
! PURPOSE:
!
! Driver subroutine for Betts-Miller-Janjic deep and shallow convective
! adjustment scheme.
!
!------------------------------------------------------------------------
!
! AUTHOR:  ???
!
! MODIFICATIONS:
!
! Eric Kemp, 21 September 2001
! Reformatted code.
!
!------------------------------------------------------------------------
!
! Force explicit declarations.
!
!------------------------------------------------------------------------

  IMPLICIT NONE

!------------------------------------------------------------------------
!
! Declare arguments
!
!------------------------------------------------------------------------

!  INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE                        &
!                           ,IMS,IME,JMS,JME,KMS,KME                     & 
!                           ,ITS,ITE,JTS,JTE,KTS,KTE
  INTEGER,INTENT(IN) :: IDS, & ! start index for i in domain
                        IDE, & ! end index for i in domain
                        JDS, & ! start index for j in domain
                        JDE, & ! end index for j in domain
                        KDS, & ! start index for k in domain
                        KDE, & ! end index for k in domain
                        IMS, & ! start index for i in memory
                        IME, & ! end index for i in memory
                        JMS, & ! start index for j in memory
                        JME, & ! end index for j in memory
                        KMS, & ! start index for k in memory
                        KME, & ! end index for k in memory
                        ITS, & ! start index for i in tile
                        ITE, & ! end index for i in tile
                        JTS, & ! start index for j in tile
                        JTE, & ! end index for j in tile
                        KTS, & ! start index for k in tile
                        KTE    ! end index for k in tile

  INTEGER,INTENT(IN) :: ITIMESTEP, &  ! Number of time step (integer)
                        STEPCU        ! Number of fundamental timesteps 
                                      !   between convection calls.

  INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: &
                     LOWLYR ! Index of lowest model layer above the ground

!  REAL,INTENT(IN) :: CP,DT,ELWV,G,R,RVOVRD,TFRZ,D608
  REAL,INTENT(IN) :: CP, &      ! specific heat at constant pressure 
                                !   (1004 J/k/kg)
                     DT, &      ! Time step (sec)
                     ELWV, &    ! Latent heat of vaporization at 0 deg C
                                !   (2.5e6 J kg**-1)
                     G,  &      ! Acceleration due to gravity 
                                !   (9.81 m s**-2)
                     R,  &      ! Dry gas constant (287 J K**-1 kg**-1)
                     RVOVRD, &  ! Ratio of gas constant for water vapor
                                !   over dry gas constant
                                !   (461.6/287 J K**-1 kg**-1
                     TFRZ, &    ! Freezing point (273.15 K)
                     D608       ! rvovrd - 1.

  REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: &
                     XLAND      ! Land-sea mask (1.0 for land; 2.0 for 
                                !   water)

  REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: &
                     TH, &      ! Potential temperature (K)
                     T, &       ! Temperature (K)
                     QV, &      ! Water vapor mixing ratio (kg/kg)
                     RR, &      ! Dry air density (kg/m^3)
                     PINT, &    ! Pressure at full levels (Pa)
                     PMID, &    ! Pressure (Pa)
                     PI, &      ! Exner function (dimensionless)
                     RHO, &     ! Density (kg/m^3)
                     DZ8W       ! dz between full levels (m)

  REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(INOUT) ::               &
                     RQVCUTEN,& ! Rho_dTheta_m tendency due to
                                !   cumulus scheme precipitation 
                                !   (kg/m^3 . K)
                     RTHCUTEN   ! Rho_dQv tendency due to
                                !   cumulus scheme precipitation 
                                !   (kg/m^3 . kg/kg)

  REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: &
                     RAINCV, &  ! Cumulus scheme precipitation (mm)
                     CLDEFI     ! Precipitation efficiency (for BMJ 
                                !   scheme) (dimensionless)

!------------------------------------------------------------------------
!
! Local variables.
!
!------------------------------------------------------------------------

  REAL :: CUBOT,CUTOP,DTCNVC,LANDMASK,PCPCOL,PSFC,PTOP

  REAL,DIMENSION(KTS:KTE) :: DPCOL,DQDT,DTDT,PCOL,QCOL,TCOL

  INTEGER :: I,J,K,KFLIP,ICLDCK,LBOT,LMH,LTOP

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
! 
! Beginning of executable code...
! 
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

!------------------------------------------------------------------------
!
! Check to see if this is a convection timestep.
!
!------------------------------------------------------------------------

  ICLDCK=MOD(ITIMESTEP,STEPCU)                                              

!------------------------------------------------------------------------
!
! Compute convection every stepcu*dt/60.0 minutes.
!
!------------------------------------------------------------------------

  IF (ICLDCK.EQ.0) THEN

    DTCNVC=DT*STEPCU

    DO J=JTS,JTE  
      DO I=ITS,ITE
        DO K=KTS,KTE
          DQDT(K)=0. 
          DTDT(K)=0.
        END DO ! DO k = kts,kte

        RAINCV(I,J)=0.
        PCPCOL=0.
        PSFC=PINT(I,LOWLYR(I,J),J)  ! Surface pressure (Pa)
        PTOP=PINT(I,KTE+1,J)        ! Top level pressure (Pa)  KTE+1=KME

!------------------------------------------------------------------------
!
!       Convert to BMJ land mask.  (1.0 for sea, 0.0 for land)
!
!------------------------------------------------------------------------

        LANDMASK=XLAND(I,J)-1.

!------------------------------------------------------------------------
!
!       Fill 1-D vertical arrays and flip direction since BMJ scheme
!       counts downward from the domain's top.
!
!------------------------------------------------------------------------

        DO K=KTS,KTE
          KFLIP=KTE+1-K
          TCOL(K)=T(I,KFLIP,J)  ! Local temperature (K)

! change 10/18/00 
! from mixing ratio to specific humidity (prevent -ve 4/23/01)

          QCOL(K)=MAX(0.,QV(I,KFLIP,J)/(1.+QV(I,KFLIP,J))) ! Spec. Humidity
          PCOL(K)=PMID(I,KFLIP,J)                      ! Pressure
          DPCOL(K)=PINT(I,KFLIP,J)-PINT(I,KFLIP+1,J)
          DPCOL(K)=RHO(I,KFLIP,J)*G*DZ8W(I,KFLIP,J)   ! Pressure thickness
                                                      ! of layer.
        END DO ! DO k=kts,kte

!------------------------------------------------------------------------
!
!       Lowest layer above ground must also be flipped.
!
!------------------------------------------------------------------------

        LMH=KTE+1-LOWLYR(I,J)

!------------------------------------------------------------------------
!
!       Run BMJ scheme.
!
!------------------------------------------------------------------------

        CALL BMJ(I,J,DTCNVC,LMH,LANDMASK,CLDEFI(I,J)                     &
                  ,DPCOL,PCOL,QCOL,TCOL,PSFC,PTOP                        &
                  ,DQDT,DTDT,PCPCOL,LBOT,LTOP                            &
                  ,CP,R,ELWV,G,TFRZ,D608                                 &
                  ,IDS,IDE,JDS,JDE,KDS,KDE                               &
                  ,IMS,IME,JMS,JME,KMS,KME                               &
                  ,ITS,ITE,JTS,JTE,KTS,KTE)

!------------------------------------------------------------------------
!
!       Compute heating and moistening tendencies.
!
!------------------------------------------------------------------------

        DO K=KTS,KTE
          KFLIP=KTE+1-K

          RTHCUTEN(I,K,J)=RR(I,K,J)*                                     &
                          ((1.+RVOVRD*QV(I,K,J))*DTDT(KFLIP)             &
                                                /PI(I,K,J)               &
                             +RVOVRD*TH(I,K,J)*DQDT(KFLIP))

! change 10/18/00 
! from specific humidity back to mixing ration

          RQVCUTEN(I,K,J)=RR(I,K,J)*DQDT(KFLIP)/(1.-QCOL(KFLIP)**2)
        END DO ! DO k = kts,kte

!------------------------------------------------------------------------
!
!       All units in BMJ scheme are MKS, so convert precip from meters
!       to millimeters per step.
!
!------------------------------------------------------------------------

        RAINCV(I,J)=PCPCOL*1.E3/STEPCU

!------------------------------------------------------------------------
!
!       Variables cubot and cutop are the bottom and top layers, 
!       respectively, of the convective cloud.  If need be, they can be 
!       used in the radiation computations. (Note:  They are not
!       used at this time! EMK)
!
!------------------------------------------------------------------------

        IF (CUBOT.NE.0) CUBOT=REAL(KTE+1-LBOT)
        IF (CUTOP.NE.0) CUTOP=REAL(KTE+1-LTOP)

      END DO ! DO i = its,ite
    END DO ! DO j = jts,jte

  END IF ! IF (ICLDCK.EQ.0) THEN

END SUBROUTINE BMJDRV

!########################################################################
!########################################################################
!#########                                                      #########
!#########                    SUBROUTINE bmj                    #########
!#########                                                      #########
!#########                      Adapted by                      #########
!#########     Center for Analysis and Prediction of Storms     #########
!#########                University of Oklahoma                #########
!#########                                                      #########
!########################################################################
!########################################################################


SUBROUTINE BMJ(I,J,DTCNVC,LMH,SM,CLDEFI                                 & 1,2
              ,DPRS,PRSMID,Q,T,PSFC,PT                                  &
              ,DQDT,DTDT,PCPCOL,LBOT,LTOP                               &
              ,CP,R,ELWV,G,TFRZ,D608                                    &
              ,IDS,IDE,JDS,JDE,KDS,KDE                                  &
              ,IMS,IME,JMS,JME,KMS,KME                                  &
              ,ITS,ITE,JTS,JTE,KTS,KTE)

!------------------------------------------------------------------------
!
! PURPOSE:
!
! Runs the Betts-Miller-Janjic deep and shallow convective adjustment
! scheme.
!
!------------------------------------------------------------------------
!
! AUTHOR: ???
!
! MODIFICATION:
!
! Eric Kemp, 21 September 2001
! Reformatted code.
!
!------------------------------------------------------------------------
!
! Force explicit declarations
!
!------------------------------------------------------------------------

  IMPLICIT NONE

!------------------------------------------------------------------------
!
! Declare arguments
!
!------------------------------------------------------------------------

  INTEGER,INTENT(IN) :: IDS, & ! start index for i in domain
                        IDE, & ! end index for i in domain
                        JDS, & ! start index for j in domain
                        JDE, & ! end index for j in domain
                        KDS, & ! start index for k in domain
                        KDE, & ! end index for k in domain
                        IMS, & ! start index for i in memory
                        IME, & ! end index for i in memory
                        JMS, & ! start index for j in memory
                        JME, & ! end index for j in memory
                        KMS, & ! start index for k in memory
                        KME, & ! end index for k in memory
                        ITS, & ! start index for i in tile
                        ITE, & ! end index for i in tile
                        JTS, & ! start index for j in tile
                        JTE, & ! end index for j in tile
                        KTS, & ! start index for k in tile
                        KTE, & ! end index for k in tile
                        I,   & ! I index of column
                        J      ! J index of column

  INTEGER,INTENT(IN) :: LMH    ! Lowest k index above ground

  INTEGER,INTENT(OUT) :: LBOT,&! K level of bottom of convection
                         LTOP  ! K level of top of convection

  REAL,INTENT(IN) :: CP, &     ! Specific heatat constant pressure
                               !   (1004 J/K/kg)
                     D608, &   ! (Rv/Rd) - 1.
                     DTCNVC, & ! Convection time step (stepcu*dt)
                     ELWV, &   ! Latent heat of vaporization at 0 deg C
                               !   (2.5e6 J/kg)
                     G, &      ! Gravitational acceleration (9.81 m/s/s)
                     PSFC, &   ! Surface (lowest layer) pressure (Pa)
                     PT, &     ! Model top pressure (Pa)
                     R, &      ! Dry gas constant (287 J/K/kg)
                     SM, &     ! Land mask (1.0 for sea, 0.0 for land)
                     TFRZ      ! Freezing point (273.15 K)

  REAL,DIMENSION(KTS:KTE),INTENT(IN) :: &
                     DPRS, &   ! Delta pressure across layer (Pa)
                     PRSMID, & ! Pressure at layer (Pa)
                     Q, &      ! Specific humidity (kg/kg)
                     T         ! Temperature (K)

  REAL,INTENT(INOUT) :: &
                     CLDEFI, & ! Precip. efficiency (dimensionless)
                     PCPCOL    ! Convective precip. (m)

  REAL,DIMENSION(KTS:KTE),INTENT(INOUT) :: &
                     DQDT, &   ! Specific humidity time tendency (kg/kg/s)
                     DTDT      ! Temperature time tendency (K/s)

!------------------------------------------------------------------------
!
! Declare local variable arrays.
!
!------------------------------------------------------------------------

  REAL,DIMENSION(KTS:KTE) :: APEK,APESK,FPK                              &
                            ,PK,PSK,QK,QREFK,QSATK                       &
                            ,THERK,THEVRF,THSK                           &
                            ,THVMOD,THVREF,TK,TREFK

  REAL,DIMENSION(KTS:KTE) :: APE,DIFQ,DIFT,TREF

!------------------------------------------------------------------------
!
! Declare local scalar variables.
!
!------------------------------------------------------------------------

 LOGICAL :: DEEP,SHALLOW

 REAL :: APEKL,APEKXX,APEKXY,APESP,APESTS                                &
             ,AVRGT,AVRGTL,BQ,BQK,BQS00K,BQS10K                          &
             ,CAPA,CTHRS,DEN,DENTPY,DEPMIN,DEPTH                         &
             ,DEPWL,DHDT,DIFQL,DIFTL,DPKL,DPMIX                          &
             ,DQREF,DRHDP,DRHEAT,DSP                                     &
             ,DSP0,DSP0K,DSPB,DSPBK,DSPT,DSPTK                           &
             ,DST,DSTQ,DTHEM,DTDP,EFI                                    &
             ,FEFI,FPTK,HCORR,OTSUM,P,P00K,P01K,P10K,P11K                &
             ,PART1,PART2,PART3,PBOT,PBTK                                &
             ,PK0,PKB,PKL,PKT,PKXXXX,PKXXXY                              &
             ,PLMH,POTSUM,PP1,PPK,PRECK                                  &
             ,PRESK,PSP,PSUM,PTHRS,PTOP,PTPK                             &
             ,QBT,QKL,QNEW,QOTSUM,QQ1,QQK,QRFKL,QRFTP,QSUM,RDP0T         &
             ,RDPSUM,RDTCNVC,RHH,RHL,RHMAX,ROTSUM,RTBAR                  &
             ,SMIX,SQ,SQK,SQS00K,SQS10K,STABDL,SUMDE,SUMDP               &
             ,SUMDT,TAUK,TCORR,THBT,THERKX,THERKY                        &
             ,THESP,THSKL,THTPK,THVMKL,TKL                               &
             ,TPSP,TQ,TQK,TREFKX,TRFKL,TSKL,TTH,TTHBT,TTHES,TTHK

  INTEGER :: IQ,IQTB,IT,ITER,ITTB,ITTBK,KB,KNUMH,KNUML                   &
                ,L,L0,L0M1,LB,LBM1,LCOR,LPT1                             &
                ,LQM,LSHU,LTP1,LTP2,LTSH

!------------------------------------------------------------------------
!
! Declare local parameters.
!
!------------------------------------------------------------------------

  REAL,PARAMETER :: DSPBSL=DSPBFL*FSL,DSP0SL=DSP0FL*FSL                  &
                       ,DSPTSL=DSPTFL*FSL                                &
                       ,DSPBSS=DSPBFS*FSS,DSP0SS=DSP0FS*FSS              &
                       ,DSPTSS=DSPTFS*FSS

  REAL,PARAMETER :: AVGEFI=(EFIMN+1.)*.5,STEFI=1.

  REAL,PARAMETER :: SLOPBL=(DSPBFL-DSPBSL)/(1.-EFIMN)                    &
                       ,SLOP0L=(DSP0FL-DSP0SL)/(1.-EFIMN)                &
                       ,SLOPTL=(DSPTFL-DSPTSL)/(1.-EFIMN)                &
                       ,SLOPBS=(DSPBFS-DSPBSS)/(1.-EFIMN)                &
                       ,SLOP0S=(DSP0FS-DSP0SS)/(1.-EFIMN)                &
                       ,SLOPTS=(DSPTFS-DSPTSS)/(1.-EFIMN)                &
                       ,SLOPE=(1.-EFMNT)/(1.-EFIMN)

  REAL :: A23M4L,CPRLG,ELOCP,FCB,RCP

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
! 
! Beginning of executable code...
! 
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

  CAPA=R/CP
  CPRLG=CP/(ROW*G*ELWV)
  ELOCP=ELIVW/CP
  RCP=1./CP
  A23M4L=A2*(A3-A4)*ELWV
  FCB=1.-FCC
  RDTCNVC=1./DTCNVC

  DEEP=.FALSE.
  SHALLOW=.FALSE.

  DSP0=0.
  DSPB=0.
  DSPT=0.
  PSP=0.

  TAUK=DTCNVC/TREL
  CTHRS=(0.006350/86400.)*DTCNVC/CPRLG

!------------------------------------------------------------------------
!
! Preparations. (???)
!
!------------------------------------------------------------------------

  LBOT=LMH
  THESP=0.
  PCPCOL=0.
  TREF(1)=T(1)

  DO L=1,LMH
    APESTS=PRSMID(L)
    APE(L)=(1.E5/APESTS)**CAPA
  END DO ! DO l = 1,lmh

!------------------------------------------------------------------------
!
! Search for maximum buoyancy level.
!
!------------------------------------------------------------------------

  DO 170 KB=1,LMH

!------------------------------------------------------------------------
!
!   Trial maximum buoyancy level variables.
!
!------------------------------------------------------------------------

    PKL=PRSMID(KB)
    PLMH=PRSMID(LMH)

!------------------------------------------------------------------------
!
!   Search over a scaled depth (between bottom pressure level and 80%
!   of the bottom pressure level) to find the parcel with the maximum
!   equivalent potential temperature.
!
!------------------------------------------------------------------------

    IF (PKL.GE.0.80*PLMH) THEN
      QBT=Q(KB)
      TTHBT=T(KB)*APE(KB)   ! T to Theta
      TTH=(TTHBT-THL)*RDTH   
      QQ1=TTH-AINT(TTH)     ! Remainder used with (bi)linear 
                            ! interpolation.
      ITTB=INT(TTH)+1       ! Look-up table index.

!------------------------------------------------------------------------
!
!     Keep indices within the table. (A look-up table is used for this
!     code --EMK)
!
!------------------------------------------------------------------------

      IF (ITTB.LT.1) THEN
        ITTB=1
        QQ1=0.
      END IF ! IF (ITTB.LT.1) THEN

      IF (ITTB.GE.JTB) THEN
        ITTB=JTB-1
        QQ1=0.
      END IF ! IF (ITTB.GE.JTB) THEN

!------------------------------------------------------------------------
!
!     Base and scaling factor for saturation specific humidity.
!
!------------------------------------------------------------------------

      ITTBK=ITTB
      BQS00K=QS0(ITTBK)
      SQS00K=SQS(ITTBK)
      BQS10K=QS0(ITTBK+1)
      SQS10K=SQS(ITTBK+1)

!------------------------------------------------------------------------
!
!     Scaling specific humidity and table index using linear 
!     interpolation.
!
!------------------------------------------------------------------------

      BQ=(BQS10K-BQS00K)*QQ1+BQS00K
      SQ=(SQS10K-SQS00K)*QQ1+SQS00K
      TQ=(QBT-BQ)/SQ*RDQ
      PP1=TQ-AINT(TQ)      ! Remainder used with bilinear interpolation.
      IQTB=INT(TQ)+1       ! Index in look-up table

!------------------------------------------------------------------------
!
!     Keep indices within the table. (Refers to look-up table -- EMK)
!
!------------------------------------------------------------------------

      IF (IQTB.LT.1) THEN
        IQTB=1
        PP1=0.
      END IF ! IF (IQTB.LT.1) THEN

      IF (IQTB.GE.ITB) THEN
        IQTB=ITB-1
        PP1=0.
      END IF ! IF (IQTB.GE.ITB) THEN

!------------------------------------------------------------------------
!
!     Saturation pressure at four surrounding table points. (Refers to
!     look-up table --EMK)
!
!------------------------------------------------------------------------

      IQ=IQTB
      IT=ITTB
      P00K=PTBL(IQ  ,IT  )
      P10K=PTBL(IQ+1,IT  )
      P01K=PTBL(IQ  ,IT+1)
      P11K=PTBL(IQ+1,IT+1)

!------------------------------------------------------------------------
!
!     Saturation point variables at the bottom.
!
!------------------------------------------------------------------------

      TPSP=P00K+(P10K-P00K)*PP1+(P01K-P00K)*QQ1                          &
               +(P00K-P10K-P01K+P11K)*PP1*QQ1 ! Saturation pressure
      APESP=(1.E5/TPSP)**CAPA                 
      TTHES=TTHBT*EXP(ELOCP*QBT*APESP/TTHBT)  ! Theta-e.

!------------------------------------------------------------------------
!
!     Check for maximum buoyancy.
!
!------------------------------------------------------------------------

      IF (TTHES.GT.THESP) THEN
        PSP  =TPSP
        THBT =TTHBT
        THESP=TTHES
      END IF ! IF (TTHES.GT.THESP) THEN

    END IF ! IF (ITTB.LT.1) THEN

  170 CONTINUE ! DO 170 KB=1,LMH

!------------------------------------------------------------------------
!
! Choose cloud base as model level just below psp.
!
!------------------------------------------------------------------------

  DO L=1,LMH-1
    P=PRSMID(L)
    IF (P.LT.PSP.AND.P.GE.PQM) LBOT=L+1 ! K level of bottom of convection
  END DO ! DO L=1,LMH-1

!------------------------------------------------------------------------
!
! WARNING:  LBOT must not be greater than LMH-1 in shallow convection
! scheme.  Make sure the cloud base is at least PONE above the surface. 
! PONE is currently set to 2500 Pa.
!
!------------------------------------------------------------------------

  PBOT=PRSMID(LBOT) ! Pressure at bottom of convection.
  PLMH=PRSMID(LMH)  ! Bottom level pressure.

  IF (PBOT.GE.PLMH-PONE.OR.LBOT.GE.LMH) THEN
    DO L=1,LMH-1
      P=PRSMID(L)
      IF(P.LT.PLMH-PONE) LBOT=L
    END DO ! DO L=1,LMH-1

    PBOT=PRSMID(LBOT)
  END IF ! IF (PBOT.GE.PLMH-PONE.OR.LBOT.GE.LMH) THEN

!------------------------------------------------------------------------
!
! Cloud top computation.
!
!------------------------------------------------------------------------

  LTOP=LBOT ! Initialize k-index of top of convection
  PTOP=PBOT ! Initialize pressure of top of convection

!------------------------------------------------------------------------
!
! Compute parcel temperature along moist adiabat, scaling pressure
! and TT table index.  (???)
!
!------------------------------------------------------------------------

  DO L=LMH,1,-1

    PRESK=PRSMID(L)

    IF (PRESK.LT.PLQ) THEN ! PLQ = 70000 Pa
      CALL TTBLEX(ITB,JTB,PL,PRESK,RDP,RDTHE                             &
                 ,STHE,THE0,THESP,TTBL,TREF(L))
    ELSE
      CALL TTBLEX(ITBQ,JTBQ,PLQ,PRESK,RDPQ,RDTHEQ                        &
                 ,STHEQ,THE0Q,THESP,TTBLQ,TREF(L))
    END IF ! IF (PRESK.LT.PLQ) THEN

  END DO ! DO L=LMH,1,-1

!------------------------------------------------------------------------
!
! Buoyancy check.
!
!------------------------------------------------------------------------

  DO L=LMH,1,-1
    IF (TREF(L).GT.T(L)-DTTOP) LTOP=L ! K-index of top of cloud
  END DO ! DO L=LMH,1,-1

!------------------------------------------------------------------------
!
! Cloud top pressure.
!
!------------------------------------------------------------------------

  PTOP=PRSMID(LTOP)

!------------------------------------------------------------------------
!
! Define and smooth dsps and cldefi.  Unified or separate land/sea
! conv 
!
! DSPB is the saturation pressure difference at cloud base.
! DSP0 is the saturation pressure difference at freezing level.
! DSPT is the saturation pressure difference at cloud top.
!
!------------------------------------------------------------------------

  EFI=CLDEFI

  IF (UNIS) THEN ! Currently, UNIS is set to true.

!------------------------------------------------------------------------
!
!   All profiles are sea profiles (unified sea) (???)
! 
!------------------------------------------------------------------------

    DSPB=(EFI-EFIMN)*SLOPBS+DSPBSS
    DSP0=(EFI-EFIMN)*SLOP0S+DSP0SS
    DSPT=(EFI-EFIMN)*SLOPTS+DSPTSS

  ELSE IF (.NOT.UNIL) THEN ! Currently, UNIL is set to false

!------------------------------------------------------------------------
! 
!   All profiles are not land profiles (not unified land) (???)
!
!------------------------------------------------------------------------

    DSPB=((EFI-EFIMN)*SLOPBS+DSPBSS)*SM                                  &
            +((EFI-EFIMN)*SLOPBL+DSPBSL)*(1.-SM)
    DSP0=((EFI-EFIMN)*SLOP0S+DSP0SS)*SM                                  &
            +((EFI-EFIMN)*SLOP0L+DSP0SL)*(1.-SM)
    DSPT=((EFI-EFIMN)*SLOPTS+DSPTSS)*SM                                  &
                 +((EFI-EFIMN)*SLOPTL+DSPTSL)*(1.-SM)
  ELSE

!------------------------------------------------------------------------
! 
!   All profiles are land profiles (unified land) (???)
!
!------------------------------------------------------------------------

    DSPB=((EFI-EFIMN)*SLOPBL+DSPBSL)
    DSP0=((EFI-EFIMN)*SLOP0L+DSP0SL)
    DSPT=((EFI-EFIMN)*SLOPTL+DSPTSL) 
  END IF ! IF (UNIS) ... ELSE IF (.NOT. UNIL) ...

!------------------------------------------------------------------------
! 
! Clean up and gather deep convection points.
!
!------------------------------------------------------------------------

  IF (LTOP.GE.LBOT) THEN
    LBOT=0
    LTOP=LBOT
    PTOP=PBOT
  END IF

  IF (PTOP.GT.PBOT-PNO.OR.LTOP.GT.LBOT-2)                                &
       CLDEFI=AVGEFI*SM+STEFI*(1.-SM) ! Recalculate precip. efficiency

!------------------------------------------------------------------------
! 
! Depth of cloud required to make the point a deep convection point is
! a scaled value of psfc.
!
!------------------------------------------------------------------------

  DEPMIN=PSH*PSFC*1.E-5
  DEPTH=PBOT-PTOP

  IF (DEPTH.GE.DEPMIN) THEN
    DEEP=.TRUE.
  END IF ! IF (DEPTH.GE.DEPMIN) THEN

!------------------------------------------------------------------------
!
! Deep convection.
!
!------------------------------------------------------------------------

  IF (.NOT.DEEP) GO TO 600 ! Skip to shallow convection

  LB   =LBOT
  EFI  =CLDEFI
  DSPBK=DSPB 
  DSP0K=DSP0
  DSPTK=DSPT

!------------------------------------------------------------------------
! 
! Initialize variables in the convective column.  One should note that
! the values assigned to the array trefk in the following loop are 
! really only relevant in anchoring the reference temperature profile
! at level lb.  When building the reference profile from cloud base,
! then assigning the ambient temperature to trefk is acceptable.  
! However, when building the reference profile from some other level
! (such as one level above the ground), then trefk should be filled with
! the temperatures in tref(l) which are the temperatures of the moist
! adiabat through cloud base.  By the time the line numbered 450 has 
! been reached, trefk actually does hold the reference temperature
! profile.
!
!------------------------------------------------------------------------

  DO L=1,LMH
    DIFT  (L)=0.
    DIFQ  (L)=0.
    TKL      =T(L)
    TK    (L)=TKL           ! Temperature
    TREFK (L)=TKL           ! Temperature
    QKL      =Q(L)
    QK    (L)=QKL           ! Specific humidity
    QREFK (L)=QKL           ! Specific humidity
    PKL      =PRSMID(L)
    PK    (L)=PKL           ! Pressure
    PSK   (L)=PKL           ! Pressure
    APEKL    =APE(L)
    APEK  (L)=APEKL
    THERK (L)=TREF(L)*APEKL  ! Temperature to Theta
  END DO ! DO L=1,LMH

!------------------------------------------------------------------------
!
! Deep convection reference temperature profile.
!
!------------------------------------------------------------------------
  
  LTP1=LTOP+1    ! K-index of one level below cloud top
  LBM1=LB-1      ! K-index of one level above cloud bottom
  PKB=PK(LB)     ! Pressure at cloud bottom
  PKT=PK(LTOP)   ! Pressure at cloud top

!------------------------------------------------------------------------
!
! Temperature reference profile below freezing level.
!
!------------------------------------------------------------------------

  L0=LB
  PK0=PK(LB)
  TREFKX=TREFK(LB)     ! Temperature at cloud bottom
  THERKX=THERK(LB)     ! Theta at cloud bottom
  APEKXX=APEK(LB)
  THERKY=THERK(LBM1)   ! Theta at one level above cloud bottom
  APEKXY=APEK(LBM1)

  DO L=LBM1,LTOP,-1
    IF (T(L+1).LT.TFRZ) GO TO 430 ! Skip to code for temperature profile
                                  ! above freezing level.
    STABDL=STABD
    TREFKX=((THERKY-THERKX)*STABDL                                       &
                +TREFKX*APEKXX)/APEKXY  ! Reference profile temperature
    TREFK(L)=TREFKX                     ! Reference profile temperature

    APEKXX=APEKXY
    THERKX=THERKY                       ! Change lower level theta to
                                        ! upper level theta.
    APEKXY=APEK(L-1)
    THERKY=THERK(L-1)                   ! Change upper level theta to
                                        ! theta from next highest level.
    L0=L
    PK0=PK(L0)                          ! Pressure
  END DO ! DO L=LBM1,LTOP,-1

!------------------------------------------------------------------------
!
! Freezing level at or above the cloud top.
!
!------------------------------------------------------------------------

  L0M1=L0-1 
  GO TO 450 ! Skip to deep convection reference humidity profile.

!------------------------------------------------------------------------
!
! Temperature reference profile above freezing level.
!
!------------------------------------------------------------------------

  430 L0M1=L0-1
  RDP0T=1./(PK0-PKT) ! Reciprocal of difference in pressure between 
                     ! current level and cloud top.
  DTHEM=THERK(L0)-TREFK(L0)*APEK(L0) ! Difference between theta and
                                     ! reference theta.

  DO L=LTOP,L0M1
    TREFK(L)=(THERK(L)-(PK(L)-PKT)*DTHEM*RDP0T)/APEK(L) ! Reference 
                                                        ! profile
                                                        ! temperature.
  END DO ! DO L=LTOP,L0M1

!------------------------------------------------------------------------
!
! Deep convection reference humidity profile.
!
! DEPWL is the pressure difference between cloud base and the freezing
! level.
!
!------------------------------------------------------------------------

  450 DEPWL=PKB-PK0
  DEPTH=PFRZ*PSFC*1.E-5 ! 15000 Pa * lowest layer pressure * 1.E-5

!------------------------------------------------------------------------
!
! Saturation pressure difference.
!
!------------------------------------------------------------------------

  DO 460 L=LTOP,LB ! Cloud top to cloud base.

    IF (DEPWL.GE.DEPTH) THEN
      IF (L.LT.L0) THEN ! L0 is k-index of freezing level
        DSP=((PK0-PK(L))*DSPTK+(PK(L)-PKT)*DSP0K)/(PK0-PKT)
      ELSE
        DSP=((PKB-PK(L))*DSP0K+(PK(L)-PK0)*DSPBK)/(PKB-PK0)
      END IF ! IF (L.LT.L0) THEN
    ELSE
      DSP=DSP0K
      IF (L.LT.L0)                                                       &
        DSP=((PK0-PK(L))*DSPTK+(PK(L)-PKT)*DSP0K)/(PK0-PKT)
    END IF ! IF (DEPWL.GE.DEPTH) THEN

!------------------------------------------------------------------------
!
!   Humidity profile.
!
!------------------------------------------------------------------------

    IF (PK(L).GT.PQM) THEN ! PQM = 20000 Pa
      PSK(L)=PK(L)+DSP
      APESK(L)=(1.E5/PSK(L))**CAPA
      THSK(L)=TREFK(L)*APEK(L)     ! Reference profile theta
      QREFK(L)=PQ0/PSK(L)*EXP(A2*(THSK(L)-A3*APESK(L))                   &
                                  /(THSK(L)-A4*APESK(L))) ! Reference
                                                          ! profile
                                                          ! sp. humidity.
    ELSE
      QREFK(L)=Q(L) ! Reference profile specific humidity.
    END IF ! IF (PK(L).GT.PQM) THEN

  460 CONTINUE ! DO 460 L=LTOP,LB

!------------------------------------------------------------------------
!
! Enthalpy conservation integral.
!
!------------------------------------------------------------------------

  LQM=0

  DO 520 ITER=1,2
    SUMDE=0. ! Integral of enthalpy multiplied by delta pressure
             ! across each layer.
    SUMDP=0. ! Integral of pressure from cloud base to cloud top. 
    DO L=LTOP,LB ! Cloud top to cloud base
      SUMDE=((TK(L)-TREFK(L))*CP+(QK(L)-QREFK(L))*ELWV)*DPRS(L)          &
              +SUMDE
      SUMDP=SUMDP+DPRS(L)
    END DO ! DO L=LTOP,LB

    HCORR=SUMDE/(SUMDP-DPRS(LTOP)) ! Integrated enthalpy.
    LCOR=LTOP+1 ! One level below cloud top height.

!------------------------------------------------------------------------
! 
!   Find LQM.
!       
!------------------------------------------------------------------------

    DO L=1,LB
      IF (PK(L).LE.PQM) LQM=L ! PQM = 20000 Pa
    END DO ! DO L=1,LB

!------------------------------------------------------------------------
! 
!   Below lqm (above 200 mb level) correct temperature only.
!       
!------------------------------------------------------------------------

    IF (LCOR.LE.LQM) THEN
      DO L=LCOR,LQM
        TREFK(L)=TREFK(L)+HCORR*RCP
      END DO ! DO L=LCOR,LQM
      LCOR=LQM+1
    END IF ! IF (LCOR.LE.LQM) THEN

!------------------------------------------------------------------------
! 
!   Above lqm (below 200 mb level) correct both temperature and moisture.
!       
!------------------------------------------------------------------------

    DO L=LCOR,LB ! LQM + 1 to cloud bottom               
      TSKL=TREFK(L)*APEK(L)/APESK(L)
      DHDT=QREFK(L)*A23M4L/(TSKL-A4)**2+CP
      TREFK(L)=HCORR/DHDT+TREFK(L)
      THSKL=TREFK(L)*APEK(L)
      QREFK(L)=PQ0/PSK(L)*EXP(A2*(THSKL-A3*APESK(L))                     &
                                  /(THSKL-A4*APESK(L)))
    END DO ! DO L=LCOR,LB

  520 CONTINUE ! DO 520 ITER=1,2

!------------------------------------------------------------------------
! 
! Heating, moistening, precipitation.  (???)
!       
!------------------------------------------------------------------------

  DENTPY=0. ! Integral of difference in entropy between actual and
            ! reference thermodynamic profile.
  AVRGT =0.
  PRECK =0. ! Condensate/precipitation

  DO L=LTOP,LB ! Cloud top to cloud base
    TKL    =TK(L)
    DIFTL  =(TREFK(L)-TKL  )*TAUK ! TAUK=DTCNVC/TREL
    DIFQL  =(QREFK(L)-QK(L))*TAUK ! TAUK=DTCNVC/TREL
    AVRGTL =(TKL+TKL+DIFTL)
    DENTPY =(DIFTL*CP+DIFQL*ELWV)*DPRS(L)/AVRGTL+DENTPY
    AVRGT  =AVRGTL*DPRS(L)+AVRGT
    PRECK  =DPRS(L)*DIFTL+PRECK
    DIFT(L)=DIFTL
    DIFQ(L)=DIFQL
  END DO ! DO L=LTOP,LB

  DENTPY=DENTPY+DENTPY
  AVRGT =AVRGT/(SUMDP+SUMDP)

!------------------------------------------------------------------------
!
! Swap if entropy and/or precip is less than zero.
!
!------------------------------------------------------------------------

  IF (DENTPY.LT.EPSNTP.OR.PRECK.LT.0.) THEN ! EPSNTP = 1.E2
    CLDEFI=EFIMN*SM+STEFI*(1.-SM) ! New precip. efficiency

!------------------------------------------------------------------------
!
!   Search for shallow cloud top.
!
!------------------------------------------------------------------------

    LTSH=LBOT                 ! K-index of cloud bottom
    LBM1=LBOT-1               ! K-index of level above cloud bottom
    PBTK=PK(LBOT)             ! Pressure at cloud bottom
    DEPMIN=PSH*PSFC*1.E-5
    PTPK=PBTK-DEPMIN          ! Est. pressure at shallow cloud top

!------------------------------------------------------------------------
!
!   Cloud top is the level just below pbtk-psh.
!
!------------------------------------------------------------------------

    DO L=1,LMH
      IF(PK(L).LE.PTPK)LTOP=L+1 ! K-index of shallow cloud top
    END DO ! DO L=1,LMH

    PTPK=PK(LTOP) ! Pressure of shallow cloud top

!------------------------------------------------------------------------
!
!   Highest level allowed is level just below pshu.
!
!------------------------------------------------------------------------

    IF (PTPK.LE.PSHU) THEN ! PSHU = 45000 Pa

      DO L=1,LMH
        IF (PK(L).LE.PSHU) LSHU=L+1
      END DO ! DO L=1,LMH

      LTOP=LSHU            ! K-index of shallow cloud top
      PTPK=PK(LTOP)        ! Pressure of shallow cloud top
    END IF ! IF (PTPK.LE.PSHU) THEN

    LTP1=LTOP+1    ! K-index of first level below shallow cloud top
    LTP2=LTOP+2    ! K-index of second level below shallow cloud top

    DO L=LTOP,LBOT
      QSATK(L)=PQ0/PK(L)*EXP(A2*(TK(L)-A3)/(TK(L)-A4)) ! Saturation
                                                       ! sp. humidity.
    END DO ! DO L=LTOP,LBOT

    RHH=QK(LTOP)/QSATK(LTOP) ! Relative humidity at shallow cloud top
    RHMAX=0.

    DO L=LTP1,LBM1 ! Cloud top to cloud bottom
      RHL=QK(L)/QSATK(L)     ! Relative humidity
      DRHDP=(RHH-RHL)/(PK(L-1)-PK(L)) ! Lapse rate of relative humidity

      IF (DRHDP.GT.RHMAX) THEN
        LTSH=L-1     ! K-index of maximum relative humidity
        RHMAX=DRHDP
      END IF ! IF (DRHDP.GT.RHMAX) THEN

      RHH=RHL
    END DO ! DO L=LTP1,LBM1

    LTOP=LTSH ! K-index of shallow cloud top

!------------------------------------------------------------------------
!
!   Cloud must be at least two layers thick.
!
!------------------------------------------------------------------------

    IF (LBOT-LTOP.LT.2) LTOP=LBOT-2

    PTOP=PK(LTOP) ! Pressure of shallow cloud top
    GO TO 600 ! Go to shallow convection
  END IF ! IF (DENTPY.LT.EPSNTP.OR.PRECK.LT.0.) THEN

!------------------------------------------------------------------------
!
! Deep convection otherwise.
!
! Keep the land value of efi equal to 1 until precip surpasses
! a threshold value, currently set to 0.25 in per 24 hrs.
!
!------------------------------------------------------------------------

  PTHRS=CTHRS
  DRHEAT=(PRECK*SM+AMAX1(EPSP,PRECK-PTHRS)*(1.-SM))                      &
             *CP/AVRGT ! Latent heat release
  EFI=EFIFC*DENTPY/DRHEAT ! Precip. efficiency

!------------------------------------------------------------------------
!
! Unified or separate land/sea conv. (???)
!
!------------------------------------------------------------------------

  EFI=CLDEFI*FCB+EFI*FCC ! Precip. efficiency, between 0.20 and 1.

  IF(EFI.GT.1.   )EFI=1.
  IF(EFI.LT.EFIMN)EFI=EFIMN
  IF(PRECK.EQ.0.) EFI=1.
  CLDEFI=EFI

  FEFI=EFMNT+SLOPE*(EFI-EFIMN)
!  FEFI=AMAX1(EFI,EFMNT)
!
  PRECK=PRECK*FEFI

!------------------------------------------------------------------------
!
! Precipitation, temperature and moisture changes.
!
!------------------------------------------------------------------------

  PCPCOL=PRECK*CPRLG ! Precipitation accumuluation.

  DO L=LTOP,LB
    DTDT(L)=DIFT(L)*FEFI*RDTCNVC  ! Temperature tendency
    DQDT(L)=DIFQ(L)*FEFI*RDTCNVC  ! Moisture tendency
  END DO ! DO L=LTOP,LB

  600 CONTINUE ! End of deep convection

!------------------------------------------------------------------------
!
! Gather shallow convection points.
!
!------------------------------------------------------------------------

  IF (PTOP.LE.PBOT-PNO.AND.LTOP.LE.LBOT-2) THEN
    DEPMIN=PSH*PSFC*1.E-5 ! Minimum depth requirement

    IF (PTOP+1..GE.PBOT-DEPMIN) THEN
      SHALLOW=.TRUE.
    END IF ! IF (PTOP+1..GE.PBOT-DEPMIN) THEN

  END IF ! IF (PTOP.LE.PBOT-PNO.AND.LTOP.LE.LBOT-2) THEN

  IF (.NOT.SHALLOW) GO TO 800 ! Skip shallow convection

!------------------------------------------------------------------------
!
! Shallow convection.
!
!------------------------------------------------------------------------

  DO L=1,LMH
    TKL      =T(L)  
    TK   (L) =TKL           ! Temperature
    TREFK(L) =TKL           ! Temperature
    QKL      =Q(L)          
    QK   (L) =QKL           ! Specific humidity
    QREFK(L) =QKL           ! Specific humidity
    PKL      =PRSMID(L)
    PK   (L) =PKL           ! Pressure
    QSATK(L) =PQ0/PK(L)*EXP(A2*(TK(L)-A3)/(TK(L)-A4))
                            ! Saturation specific humidity
    APEKL    =APE(L)
    APEK (L) =APEKL
    THVMKL   =TKL*APEKL*(QKL*D608+1.) ! Theta-v
    THVREF(L)=THVMKL                  ! Theta-v
  END DO ! DO L=1,LMH

!------------------------------------------------------------------------
!
! Shallow cloud top.
!
!------------------------------------------------------------------------

  LBM1=LBOT-1     ! One k-level above cloud bottom
  PTPK=PTOP       ! Pressure at cloud top
  LTP1=LTOP-1     ! One k-level above cloud top

! NOTE:  This appears to be redundant...EMK
  IF (PTOP.GT.PBOT-PNO.OR.LTOP.GT.LBOT-2) THEN
    LBOT=0
    LTOP=LBOT
    PTOP=PBOT
    GO TO 800 ! Skip shallow convection
  END IF ! IF (PTOP.GT.PBOT-PNO.OR.LTOP.GT.LBOT-2) THEN

!------------------------------------------------------------------------
!
! Scaling potential temperature and table index at top.
!
!------------------------------------------------------------------------

  THTPK=T(LTP1)*APE(LTP1) ! Theta at one level above cloud top

  TTHK=(THTPK-THL)*RDTH   
  QQK =TTHK-AINT(TTHK)    ! Remainder for linear interpolation.
  IT  =INT(TTHK)+1        ! Index for look-up table

  IF (IT.LT.1) THEN
    IT=1
    QQK=0.
  END IF ! IF (IT.LT.1) THEN

  IF (IT.GE.JTB) THEN
    IT=JTB-1
    QQK=0.
  END IF ! IF (IT.GE.JTB) THEN

!------------------------------------------------------------------------
!
! Base and scaling factor for saturation specific humidity at top.
!
!------------------------------------------------------------------------

  BQS00K=QS0(IT)
  SQS00K=SQS(IT)
  BQS10K=QS0(IT+1)
  SQS10K=SQS(IT+1)

!------------------------------------------------------------------------
!
! Scaling saturation specific humidity and table index at top.
!
!------------------------------------------------------------------------

  BQK=(BQS10K-BQS00K)*QQK+BQS00K
  SQK=(SQS10K-SQS00K)*QQK+SQS00K

!  TQK=(Q(LTOP)-BQK)/SQK*RDQ
  TQK=(Q(LTP1)-BQK)/SQK*RDQ

  PPK=TQK-AINT(TQK)             ! Remainder for linear interpolation.
  IQ =INT(TQK)+1                ! Index for look-up table.

  IF (IQ.LT.1) THEN
    IQ=1
    PPK=0.
  END IF ! IF (IQ.LT.1) THEN

  IF (IQ.GE.ITB) THEN
    IQ=ITB-1
    PPK=0.
  END IF ! IF (IQ.GE.ITB) THEN

!------------------------------------------------------------------------
!
! Cloud top saturation point pressure.
!
!------------------------------------------------------------------------

  PART1=(PTBL(IQ+1,IT)-PTBL(IQ,IT))*PPK
  PART2=(PTBL(IQ,IT+1)-PTBL(IQ,IT))*QQK
  PART3=(PTBL(IQ  ,IT  )-PTBL(IQ+1,IT  )                                 &
         -PTBL(IQ  ,IT+1)+PTBL(IQ+1,IT+1))*PPK*QQK
  PTPK=PTBL(IQ,IT)+PART1+PART2+PART3 ! Saturation point pressure


  DPMIX=PTPK-PSP ! Difference in saturation pressure along mixing line 
                 ! between cloud base and cloud top.
  IF (ABS(DPMIX).LT.3000.) DPMIX=-3000.

!------------------------------------------------------------------------
!
! Temperature profile slope.
!
!------------------------------------------------------------------------

  SMIX=(THTPK-THBT)/DPMIX*STABS

  TREFKX=TREFK(LBOT+1)  ! Reference profile temperature one level below
                        ! cloud bottom.
  PKXXXX=PK(LBOT+1)     ! Pressure one level below cloud bottom.
  PKXXXY=PK(LBOT)       ! Pressure at cloud bottom.
  APEKXX=APEK(LBOT+1)
  APEKXY=APEK(LBOT)

  DO L=LBOT,LTOP,-1
        TREFKX=((PKXXXY-PKXXXX)*SMIX                                     &
                +TREFKX*APEKXX)/APEKXY
    TREFK(L)=TREFKX ! Reference profile temperature.
    APEKXX=APEKXY
    PKXXXX=PKXXXY   ! Switch lower-level pressure with upper-level
                    ! pressure.
    APEKXY=APEK(L-1)
    PKXXXY=PK(L-1)  ! Switch upper-level pressure with pressure at
                    ! next highest level.
  END DO ! DO L=LBOT,LTOP,-1

!------------------------------------------------------------------------
!
! Temperature reference profile correction.
!
!------------------------------------------------------------------------

  SUMDT=0. ! Integral of temperature difference times pressure thickness
           ! of layers.
  SUMDP=0. ! Integral of thickness of pressure of layers.

  DO L=LTOP,LBOT
    SUMDT=(TK(L)-TREFK(L))*DPRS(L)+SUMDT
    SUMDP=SUMDP+DPRS(L)
  END DO ! DO L=LTOP,LBOT

  RDPSUM=1./SUMDP
  FPK(LBOT)=TREFK(LBOT)

  TCORR=SUMDT*RDPSUM ! Integral of temperature difference.

  DO L=LTOP,LBOT
    TRFKL   =TREFK(L)+TCORR
    TREFK(L)=TRFKL   ! Reference profile temperature.
    FPK  (L)=TRFKL
  END DO ! DO L=LTOP,LBOT

!------------------------------------------------------------------------
!
! Humidity profile equations.
!
!------------------------------------------------------------------------

  PSUM  =0. ! Integral of pressure differences 1 (pressure of layer and
            ! pressure of cloud top) times pressure difference 2 
            ! (thickness of layers).
  QSUM  =0. ! Integral of specific humidity times pressure thickness
            ! of layers.
  POTSUM=0.
  QOTSUM=0.
  OTSUM =0.
  DST   =0. ! Change in entropy.
  FPTK  =FPK(LTOP) ! Reference profile temperature at cloud top.

  DO L=LTOP,LBOT ! Cloud top to cloud base
    DPKL  =FPK(L)-FPTK
    PSUM  =DPKL *DPRS(L)+PSUM
    QSUM  =QK(L)*DPRS(L)+QSUM
    RTBAR =2./(TREFK(L)+TK(L))
    OTSUM =DPRS(L)*RTBAR+OTSUM
    POTSUM=DPKL    *RTBAR*DPRS(L)+POTSUM
    QOTSUM=QK(L)   *RTBAR*DPRS(L)+QOTSUM
    DST   =(TREFK(L)-TK(L))*RTBAR*DPRS(L)+DST
  END DO ! DO L=LTOP,LBOT

  PSUM  =PSUM*RDPSUM
  QSUM  =QSUM*RDPSUM
  ROTSUM=1./OTSUM
  POTSUM=POTSUM*ROTSUM
  QOTSUM=QOTSUM*ROTSUM
  DST   =DST*ROTSUM*CP/ELWV ! Change in entropy

!------------------------------------------------------------------------
!
! Ensure positive entropy change.
!
!------------------------------------------------------------------------

  IF (DST.GT.0.) THEN

!    DSTQ=DST*EPSUP
    LBOT=0
    LTOP=LBOT
    PTOP=PBOT
    GO TO 800 ! Skip shallow convection

  ELSE
    DSTQ=DST*EPSDN ! EPSDN = 1.05
  END IF ! IF (DST.GT.0.) THEN

!------------------------------------------------------------------------
!
! Check for isothermal atmosphere.
!
!------------------------------------------------------------------------

  DEN=POTSUM-PSUM

  IF (-DEN/PSUM.LT.5.E-5) THEN
    LBOT=0
    LTOP=LBOT
    PTOP=PBOT
    GO TO 800 ! Skip shallow convection

  ELSE

!------------------------------------------------------------------------
!
! Slope of the reference humidity profile.
!
!------------------------------------------------------------------------

    DQREF=(QOTSUM-DSTQ-QSUM)/DEN
  END IF ! IF (-DEN/PSUM.LT.5.E-5) THEN

!------------------------------------------------------------------------
!
! Humidity does not increase with height.
!
!------------------------------------------------------------------------

  IF (DQREF.LT.0.) THEN
    LBOT=0
    LTOP=LBOT
    PTOP=PBOT
    GO TO 800 ! Skip shallow convection
  END IF ! IF (DQREF.LT.0.) THEN

!------------------------------------------------------------------------
!
! Humidity at the cloud top.
!
!------------------------------------------------------------------------

  QRFTP=QSUM-DQREF*PSUM

!------------------------------------------------------------------------
!
! Humidity profile.
!
!------------------------------------------------------------------------

  DO L=LTOP,LBOT ! Cloud top to cloud base
    QRFKL=(FPK(L)-FPTK)*DQREF+QRFTP

!------------------------------------------------------------------------
!
!   Supersaturation or negative q not allowed.
!
!------------------------------------------------------------------------

    QNEW=(QRFKL-QK(L))*TAUK+QK(L) ! New specific humidity.

    IF (QNEW.GT.QSATK(L)*STRESH .OR. QNEW.LT.0.) THEN
      LBOT=0
      LTOP=LBOT
      PTOP=PBOT
      GO TO 800 ! Skip shallow convection
    END IF ! IF (QNEW.GT.QSATK(L)*STRESH .OR. QNEW.LT.0.) THEN

    THVREF(L)=TREFK(L)*APEK(L)*(QRFKL*D608+1.) ! New theta-v.
    QREFK(L)=QRFKL  ! New reference profile specific humidity.
  END DO ! DO L=LTOP,LBOT

!------------------------------------------------------------------------
!
! Eliminate impossible slopes.  (Betts, dtheta/dq)
!
!------------------------------------------------------------------------

  DO L=LTOP,LBOT ! Cloud top to cloud base.
    DTDP=(THVREF(L-1)-THVREF(L))/(PRSMID(L)-PRSMID(L-1)) ! Dthetav/DP

    IF (DTDP.LT.EPSDT) THEN ! EPSDT = 0.
      LBOT=0
      LTOP=LBOT
      PTOP=PBOT
      GO TO 800 ! Skip shallow convection
    END IF ! IF (DTDP.LT.EPSDT) THEN

  END DO ! DO L=LTOP,LBOT

  DENTPY=0.

  DO L=LTOP,LBOT ! Cloud top to cloud base.
    DENTPY=((TREFK(L)-TK(L))*CP+(QREFK(L)-QK(L))*ELWV)                   &
               /(TK(L)+TREFK(L))*DPRS(L)+DENTPY  ! Change in entropy
  END DO ! DO L=LTOP,LBOT 

!------------------------------------------------------------------------
!
! Relaxation toward reference profiles.
!
!------------------------------------------------------------------------

  DO L=LTOP,LBOT ! Cloud top to cloud base
    DTDT(L)=(TREFK(L)-TK(L))*TAUK*RDTCNVC  ! Temperature tendency
    DQDT(L)=(QREFK(L)-QK(L))*TAUK*RDTCNVC  ! Specific humidity tendency
  END DO ! DO L=LTOP,LBOT

!------------------------------------------------------------------------
!
! End of shallow convection.
!
!------------------------------------------------------------------------

  800 CONTINUE
END SUBROUTINE BMJ

!########################################################################
!########################################################################
!#########                                                      #########
!#########                  SUBROUTINE ttblex                   #########
!#########                                                      #########
!#########                      Adapted by                      #########
!#########     Center for Analysis and Prediction of Storms     #########
!#########                University of Oklahoma                #########
!#########                                                      #########
!########################################################################
!########################################################################


SUBROUTINE TTBLEX(ITBX,JTBX,PLX,PRSMID,RDPX,RDTHEX,STHE,THE0,THESP,TTBL & 4
                  ,TREF)

!------------------------------------------------------------------------
!
! PURPOSE:
!
! Extract temperature of the moist adiabat from the appropriate ttbl.
!
!------------------------------------------------------------------------
!
! AUTHOR: ???
!
! MODIFICATIONS:
!
! Eric Kemp, 24 September 2001
! Reformatted code.
!
!------------------------------------------------------------------------
!
! Force explicit declarations.
!
!------------------------------------------------------------------------

  IMPLICIT NONE

!------------------------------------------------------------------------
!
! Declare arguments.
!
!------------------------------------------------------------------------

  INTEGER,INTENT(IN) :: ITBX,JTBX  ! Dimensions of look-up tables

  REAL,INTENT(IN) :: PLX, &    ! Pressure reference used with look-up
                               ! table calculation.
                     PRSMID, & ! Grid point pressure (Pa).
                     RDPX, &   ! Remainder used with look-up table
                               ! calculation. 
                     RDTHEX, & ! Remainder used with look-up table
                               ! calculation.
                     THESP     ! Saturation point potential temperature 
                               ! (K).

  REAL,DIMENSION(ITBX),INTENT(IN) :: STHE, & ! Look-up table.
                                     THE0    ! Look-up table.

  REAL,DIMENSION(JTBX,ITBX),INTENT(IN) :: TTBL ! Look-up table.

  REAL,INTENT(OUT) :: TREF ! Parcel temperature (K)

!------------------------------------------------------------------------
!
! Internal variables.
!
!------------------------------------------------------------------------

  REAL :: BTHE00K,BTHE10K,BTHK,PK,PP,QQ,STHE00K,STHE10K,STHK             &
          ,T00K,T01K,T10K,T11K,TPK,TTHK

  INTEGER :: IPTB,ITHTB

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
! 
! Beginning of executable code...
! 
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

!------------------------------------------------------------------------
!
! Scaling pressure and TT table index.
!
! Determine pressure index for look-up table and remainder for (bi)linear
! interpolation.
!
!------------------------------------------------------------------------

  PK=PRSMID
  TPK=(PK-PLX)*RDPX
  QQ=TPK-AINT(TPK) ! Remainder for linear interpolation.
  IPTB=INT(TPK)+1  ! Index for look-up table.

!------------------------------------------------------------------------
!
! Keeping indices within the table.
!
!------------------------------------------------------------------------

  IF (IPTB.LT.1) THEN
    IPTB=1  ! Index for look-up table.
    QQ=0.   ! Remainder for linear interpolation.
  END IF ! IF (IPTB.LT.1) THEN

  IF (IPTB.GE.ITBX) THEN
    IPTB=ITBX-1  ! Index for look-up table.
    QQ=0.        ! Remainder for linear interpolation.
  END IF ! IF (IPTB.GE.ITBX) THEN

!------------------------------------------------------------------------
!
! Base and scaling factor for thetae.
!
!------------------------------------------------------------------------

  BTHE00K=THE0(IPTB)
  STHE00K=STHE(IPTB)
  BTHE10K=THE0(IPTB+1)
  STHE10K=STHE(IPTB+1)

!------------------------------------------------------------------------
!
! Scaling the and TT table index. (???)
!
! Use linear interpolation to determine the scale and base factors for
! thetae, then use them to determine theta-e index for look-up table
! and remainder for bilinear interpolation.
!
!------------------------------------------------------------------------

  BTHK=(BTHE10K-BTHE00K)*QQ+BTHE00K
  STHK=(STHE10K-STHE00K)*QQ+STHE00K
  TTHK=(THESP-BTHK)/STHK*RDTHEX
  PP=TTHK-AINT(TTHK)   ! Remainder for linear interpolation.
  ITHTB=INT(TTHK)+1    ! Index for look-up table.

!------------------------------------------------------------------------
!
! Keeping indices within the table.
!
!------------------------------------------------------------------------

  IF (ITHTB.LT.1) THEN
    ITHTB=1            ! Index for look-up table.
    PP=0.              ! Remainder for linear interpolation.
  END IF ! IF (ITHTB.LT.1) THEN

  IF (ITHTB.GE.JTBX) THEN
    ITHTB=JTBX-1       ! Index for look-up table.
    PP=0.              ! Remainder for linear interpolation.
  END IF ! IF (ITHTB.GE.JTBX) THEN

!------------------------------------------------------------------------
!
! Temperature at four surrounding TT table points.
!
!------------------------------------------------------------------------

  T00K=TTBL(ITHTB,IPTB)
  T10K=TTBL(ITHTB+1,IPTB)
  T01K=TTBL(ITHTB,IPTB+1)
  T11K=TTBL(ITHTB+1,IPTB+1)

!------------------------------------------------------------------------
!
! Parcel temperature, using bilinear interpolation.
!
!------------------------------------------------------------------------

  TREF=(T00K+(T10K-T00K)*PP+(T01K-T00K)*QQ                               &
       +(T00K-T10K-T01K+T11K)*PP*QQ)

END SUBROUTINE TTBLEX

!########################################################################
!########################################################################
!#########                                                      #########
!#########                  SUBROUTINE bmjinit                  #########
!#########                                                      #########
!#########                      Adapted by                      #########
!#########     Center for Analysis and Prediction of Storms     #########
!#########                University of Oklahoma                #########
!#########                                                      #########
!########################################################################
!########################################################################


SUBROUTINE BMJINIT(RTHCUTEN,RQVCUTEN,RQCCUTEN,RQRCUTEN                  & 1,3
                        ,CLDEFI,LOWLYR,CP,RD,restart                    &
                        ,IDS,IDE,JDS,JDE,KDS,KDE                        &
                        ,IMS,IME,JMS,JME,KMS,KME                        &
                        ,ITS,ITE,JTS,JTE,KTS,KTE)

!------------------------------------------------------------------------
!
! PURPOSE:
!
! Initialize variables and look-up tables used by Betts-Miller-Janjic 
! convective scheme.
!
!------------------------------------------------------------------------
!
! AUTHOR: ???
!
! MODIFICATIONS:
!
! Eric Kemp, 24 September 2001
! Reformatted code.
!
!------------------------------------------------------------------------
!
! Force explicit declarations.
!
!------------------------------------------------------------------------

  IMPLICIT NONE

!------------------------------------------------------------------------
!
! Declare arguments.
!
!------------------------------------------------------------------------

  LOGICAL , INTENT(IN) :: restart ! Restart flag.

  INTEGER,INTENT(IN) :: IDS, & ! start index for i in domain
                        IDE, & ! end index for i in domain
                        JDS, & ! start index for j in domain
                        JDE, & ! end index for j in domain
                        KDS, & ! start index for k in domain
                        KDE, & ! end index for k in domain
                        IMS, & ! start index for i in memory
                        IME, & ! end index for i in memory
                        JMS, & ! start index for j in memory
                        JME, & ! end index for j in memory
                        KMS, & ! start index for k in memory
                        KME, & ! end index for k in memory
                        ITS, & ! start index for i in tile
                        ITE, & ! end index for i in tile
                        JTS, & ! start index for j in tile
                        JTE, & ! end index for j in tile
                        KTS, & ! start index for k in tile
                        KTE    ! end index for k in tile

  REAL,INTENT(IN) :: CP, & ! specific heat at constant pressure
                           !   (1004 J/k/kg)
                     RD    ! Dry air gas constant.

  REAL,DIMENSION(IMS:,KMS:,JMS:),INTENT(OUT) ::                          &
            RTHCUTEN, & ! Rho_dTheta_m tendency due to  
                        !   cumulus scheme precipitation
                        !   (kg/m^3 . K)
            RQVCUTEN, & ! Rho_dQv tendency due to
                        !   cumulus scheme precipitation
                        !   (kg/m^3 . kg/kg)
            RQCCUTEN, & ! Rho_dQc tendency due to
                        !   cumulus scheme precipitation
                        !   (kg/m^3 . kg/kg)
            RQRCUTEN    ! Rho_dQr tendency due to
                        !   cumulus scheme precipitation
                        !   (kg/m^3 . kg/kg)

  REAL,DIMENSION(IMS:,JMS:),INTENT(OUT) :: &
            CLDEFI      ! Precipitation efficiency (for BMJ scheme) 
                        !   (dimensionless)

  INTEGER,DIMENSION(IMS:,JMS:),INTENT(OUT) :: &
            LOWLYR      ! Index of lowest model layer above the ground

!------------------------------------------------------------------------
!
! Declare internal parameters.
!
!------------------------------------------------------------------------

  REAL,PARAMETER :: ELIWV=2.683E6,EPS=1.E-9

!------------------------------------------------------------------------
!
! Declare internal variables.
!
!------------------------------------------------------------------------

  REAL, DIMENSION(JTB) :: APP,APT,AQP,AQT,PNEW,POLD,QSNEW,QSOLD          &
                             ,THENEW,THEOLD,TNEW,TOLD,Y2P,Y2T

  REAL,DIMENSION(JTBQ) :: APTQ,AQTQ,THENEWQ,THEOLDQ                      &
                             ,TNEWQ,TOLDQ,Y2TQ

  INTEGER :: I,J,K,ITF,JTF,KTF
  INTEGER :: KTH,KTHM,KTHM1,KP,KPM,KPM1

  REAL :: APE,DP,DQS,DTH,DTHE,P,QS,QS0K,SQSK,STHEK                       &
             ,TH,THE0K

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
! 
! Beginning of executable code...
! 
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
 
  JTF=MIN0(JTE,JDE-1)
  KTF=MIN0(KTE,KDE-1)
  ITF=MIN0(ITE,IDE-1)

  IF (.not.restart) THEN
    DO J=JTS,JTF
      DO K=KTS,KTF
        DO I=ITS,ITF
          RTHCUTEN(I,K,J)=0.
          RQVCUTEN(I,K,J)=0.
          RQCCUTEN(I,K,J)=0.
          RQRCUTEN(I,K,J)=0.
        END DO ! DO I=ITS,ITF
      END DO ! DO K=KTS,KTF
    END DO ! DO J=JTS,JTF

    DO J=JTS,JTF
      DO I=ITS,ITF
        CLDEFI(I,J)=1.
      END DO ! DO I=ITS,ITF
    END DO ! DO J=JTS,JTF
  END IF ! IF (.not.restart) THEN

!------------------------------------------------------------------------
!
! For now, assume sigma mode for lowest model layer.
!
!------------------------------------------------------------------------

  DO J=JTS,JTF
    DO I=ITS,ITF
      LOWLYR(I,J)=1
    END DO ! DO I=ITS,ITF
  END DO ! DO J=JTS,JTF

!------------------------------------------------------------------------
!
! Coarse look-up table for saturation point.
!
!------------------------------------------------------------------------

  KTHM=JTB
  KPM=ITB
  KTHM1=KTHM-1
  KPM1=KPM-1

  DTH=(THH-THL)/REAL(KTHM-1)
  DP =(PH -PL )/REAL(KPM -1)

  TH=THL-DTH

  DO 100 KTH=1,KTHM

    TH=TH+DTH
    P=PL-DP

    DO KP=1,KPM
      P=P+DP
      APE=(100000./P)**(RD/CP)
      QSOLD(KP)=PQ0/P*EXP(A2*(TH-A3*APE)/(TH-A4*APE)) ! Saturation 
                                                      ! specific humidity.
      POLD(KP)=P
    END DO ! DO KP=1,KPM

    QS0K=QSOLD(1)
    SQSK=QSOLD(KPM)-QSOLD(1) ! Scaling factor.
    QSOLD(1  )=0.
    QSOLD(KPM)=1.

    DO KP=2,KPM1
      QSOLD(KP)=(QSOLD(KP)-QS0K)/SQSK ! Base look-up table
      IF ((QSOLD(KP)-QSOLD(KP-1)).LT.EPS) QSOLD(KP)=QSOLD(KP-1)+EPS
    END DO ! DO KP=2,KPM1

    QS0(KTH)=QS0K  ! Base factor look-up table.
    SQS(KTH)=SQSK  ! Scaling factor look-up table.

    QSNEW(1  )=0.
    QSNEW(KPM)=1.
    DQS=1./REAL(KPM-1)

    DO KP=2,KPM1
      QSNEW(KP)=QSNEW(KP-1)+DQS
    END DO ! DO KP=2,KPM1

    Y2P(1   )=0.
    Y2P(KPM )=0.

    CALL SPLINE(JTB,KPM,QSOLD,POLD,Y2P,KPM,QSNEW,PNEW,APP,AQP)

    DO KP=1,KPM
      PTBL(KP,KTH)=PNEW(KP)
    END DO ! DO KP=1,KPM

  100 CONTINUE ! DO 100 KTH=1,KTHM

!------------------------------------------------------------------------
!
! Coarse look-up table for T(P) from constant THE.  (???)
!
!------------------------------------------------------------------------

  P=PL-DP

  DO 200 KP=1,KPM

    P=P+DP
    TH=THL-DTH

    DO KTH=1,KTHM
      TH=TH+DTH
      APE=(1.E5/P)**(RD/CP)
      QS=PQ0/P*EXP(A2*(TH-A3*APE)/(TH-A4*APE))
      TOLD(KTH)=TH/APE
      THEOLD(KTH)=TH*EXP(ELIWV*QS/(CP*TOLD(KTH)))
    END DO ! DO KTH=1,KTHM

    THE0K=THEOLD(1)
    STHEK=THEOLD(KTHM)-THEOLD(1)
    THEOLD(1   )=0.
    THEOLD(KTHM)=1.

    DO KTH=2,KTHM1
      THEOLD(KTH)=(THEOLD(KTH)-THE0K)/STHEK
      IF((THEOLD(KTH)-THEOLD(KTH-1)).LT.EPS)                             &
            THEOLD(KTH)=THEOLD(KTH-1)  +  EPS
    END DO ! DO KTH=2,KTHM1

    THE0(KP)=THE0K
    STHE(KP)=STHEK

    THENEW(1  )=0.
    THENEW(KTHM)=1.
    DTHE=1./REAL(KTHM-1)

    DO KTH=2,KTHM1
      THENEW(KTH)=THENEW(KTH-1)+DTHE
    END DO ! DO KTH=2,KTHM1

    Y2T(1   )=0.
    Y2T(KTHM)=0.

    CALL SPLINE(JTB,KTHM,THEOLD,TOLD,Y2T,KTHM,THENEW,TNEW,APT,AQT)

    DO KTH=1,KTHM
      TTBL(KTH,KP)=TNEW(KTH)
    END DO ! DO KTH=1,KTHM

  200 CONTINUE ! DO 200 KP=1,KPM

!------------------------------------------------------------------------
!
! Fine look-up table for saturation point.
!
!------------------------------------------------------------------------

  KTHM=JTBQ
  KPM=ITBQ
  KTHM1=KTHM-1
  KPM1=KPM-1

  DTH=(THHQ-THL)/REAL(KTHM-1)
  DP=(PH-PLQ)/REAL(KPM-1)

  TH=THL-DTH
  P=PLQ-DP

!------------------------------------------------------------------------
!
! Fine look-up table for T(P) from constant THE. (???)
!
!------------------------------------------------------------------------

  DO 300 KP=1,KPM

    P=P+DP
    TH=THL-DTH

    DO KTH=1,KTHM
      TH=TH+DTH
      APE=(1.E5/P)**(RD/CP)
      QS=PQ0/P*EXP(A2*(TH-A3*APE)/(TH-A4*APE))
      TOLDQ(KTH)=TH/APE
      THEOLDQ(KTH)=TH*EXP(ELIWV*QS/(CP*TOLDQ(KTH)))
    END DO ! DO KTH=1,KTHM

    THE0K=THEOLDQ(1)
    STHEK=THEOLDQ(KTHM)-THEOLDQ(1)
    THEOLDQ(1   )=0.
    THEOLDQ(KTHM)=1.

    DO KTH=2,KTHM1
      THEOLDQ(KTH)=(THEOLDQ(KTH)-THE0K)/STHEK
      IF ((THEOLDQ(KTH)-THEOLDQ(KTH-1)).LT.EPS)                          &
            THEOLDQ(KTH)=THEOLDQ(KTH-1)  +  EPS
    END DO ! DO KTH=2,KTHM1

    THE0Q(KP)=THE0K
    STHEQ(KP)=STHEK

    THENEWQ(1  )=0.
    THENEWQ(KTHM)=1.
    DTHE=1./REAL(KTHM-1)

    DO KTH=2,KTHM1
      THENEWQ(KTH)=THENEWQ(KTH-1)+DTHE
    END DO ! DO KTH=2,KTHM1

    Y2TQ(1   )=0.
    Y2TQ(KTHM)=0.

    CALL SPLINE(JTBQ,KTHM,THEOLDQ,TOLDQ,Y2TQ,KTHM                        &
                 ,THENEWQ,TNEWQ,APTQ,AQTQ)

    DO KTH=1,KTHM
      TTBLQ(KTH,KP)=TNEWQ(KTH)
    END DO ! DO KTH=1,KTHM

  300 CONTINUE

END SUBROUTINE BMJINIT

!########################################################################
!########################################################################
!#########                                                      #########
!#########                  SUBROUTINE spline                   #########
!#########                                                      #########
!#########                      Adapted by                      #########
!#########     Center for Analysis and Prediction of Storms     #########
!#########                University of Oklahoma                #########
!#########                                                      #########
!########################################################################
!########################################################################


SUBROUTINE SPLINE(JTBX,NOLD,XOLD,YOLD,Y2,NNEW,XNEW,YNEW,P,Q) 6

!------------------------------------------------------------------------
!
! PURPOSE:
!
! 1-D cubic spline fitting subroutine.
!
!------------------------------------------------------------------------
!
! AUTHOR: Z. Janjic
!
! MODIFICATIONS:
!
! Eric Kemp, 24 September 2001
! Reformatted code.  However, the overall "spagetti" structure of
! the code has been left undisturbed.
!
!------------------------------------------------------------------------
!
! Original documentation:
!
!   ******************************************************************
!   *                                                                *
!   *  THIS IS A ONE-DIMENSIONAL CUBIC SPLINE FITTING ROUTINE        *
!   *  PROGRAMED FOR A SMALL SCALAR MACHINE.                         *
!   *                                                                *
!   *  PROGRAMER Z. JANJIC                                           *
!   *                                                                *
!   *  NOLD - NUMBER OF GIVEN VALUES OF THE FUNCTION.  MUST BE GE 3. *
!   *  XOLD - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE     *
!   *         FUNCTION ARE GIVEN.  MUST BE IN ASCENDING ORDER.       *
!   *  YOLD - THE GIVEN VALUES OF THE FUNCTION AT THE POINTS XOLD.   *
!   *  Y2   - THE SECOND DERIVATIVES AT THE POINTS XOLD.  IF NATURAL *
!   *         SPLINE IS FITTED Y2(1)=0. AND Y2(NOLD)=0. MUST BE      *
!   *         SPECIFIED.                                             *
!   *  NNEW - NUMBER OF VALUES OF THE FUNCTION TO BE CALCULATED.     *
!   *  XNEW - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE     *
!   *         FUNCTION ARE CALCULATED.  XNEW(K) MUST BE GE XOLD(1)   *
!   *         AND LE XOLD(NOLD).                                     *
!   *  YNEW - THE VALUES OF THE FUNCTION TO BE CALCULATED.           *
!   *  P, Q - AUXILIARY VECTORS OF THE LENGTH NOLD-2.                *
!   *                                                                *
!   ******************************************************************
!
!------------------------------------------------------------------------
!
! Force explicit declarations.
!
!------------------------------------------------------------------------

  IMPLICIT NONE

!------------------------------------------------------------------------
!
! Declare arguments.
!
!------------------------------------------------------------------------

  INTEGER,INTENT(IN) :: JTBX,NNEW,NOLD
  REAL,DIMENSION(JTBX),INTENT(IN) :: XNEW,XOLD,YOLD
  REAL,DIMENSION(JTBX),INTENT(INOUT) :: P,Q,Y2
  REAL,DIMENSION(JTBX),INTENT(OUT) :: YNEW

!------------------------------------------------------------------------
!
! Declare internal variables.
!
!------------------------------------------------------------------------

  INTEGER :: K,K1,K2,KOLD,NOLDM1
  REAL :: AK,BK,CK,DEN,DX,DXC,DXL,DXR,DYDXL,DYDXR                        &
             ,RDX,RTDXC,X,XK,XSQ,Y2K,Y2KP1

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
! 
! Beginning of executable code...
! 
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

  NOLDM1=NOLD-1

  DXL=XOLD(2)-XOLD(1)
  DXR=XOLD(3)-XOLD(2)
  DYDXL=(YOLD(2)-YOLD(1))/DXL
  DYDXR=(YOLD(3)-YOLD(2))/DXR
  RTDXC=0.5/(DXL+DXR)

  P(1)= RTDXC*(6.*(DYDXR-DYDXL)-DXL*Y2(1))
  Q(1)=-RTDXC*DXR

  IF (NOLD.EQ.3) GO TO 150 ! Skip down

  K=3

  100 DXL=DXR
  DYDXL=DYDXR
  DXR=XOLD(K+1)-XOLD(K)
  DYDXR=(YOLD(K+1)-YOLD(K))/DXR
  DXC=DXL+DXR
  DEN=1./(DXL*Q(K-2)+DXC+DXC)

  P(K-1)= DEN*(6.*(DYDXR-DYDXL)-DXL*P(K-2))
  Q(K-1)=-DEN*DXR

  K=K+1
  IF (K.LT.NOLD) GO TO 100 ! Jump back up

  150 K=NOLDM1

  200 Y2(K)=P(K-1)+Q(K-1)*Y2(K+1)

  K=K-1
  IF (K.GT.1) GO TO 200 ! Jump back up

  K1=1

  300 XK=XNEW(K1)

  DO 400 K2=2,NOLD

    IF (XOLD(K2).GT.XK) THEN
      KOLD=K2-1
      GO TO 450 ! Exit loop, jump down
    END IF ! IF (XOLD(K2).GT.XK) THEN

  400 CONTINUE ! DO 400 K2=2,NOLD

  YNEW(K1)=YOLD(NOLD)
  GO TO 600 ! Jump down

  450 IF (K1.EQ.1) GO TO 500 ! Jump down
  IF (K.EQ.KOLD) GO TO 550 ! Jump down

  500 K=KOLD

  Y2K=Y2(K)
  Y2KP1=Y2(K+1)
  DX=XOLD(K+1)-XOLD(K)
  RDX=1./DX

  AK=.1666667*RDX*(Y2KP1-Y2K)
  BK=0.5*Y2K
  CK=RDX*(YOLD(K+1)-YOLD(K))-.1666667*DX*(Y2KP1+Y2K+Y2K)

  550 X=XK-XOLD(K)
  XSQ=X*X

  YNEW(K1)=AK*XSQ*X+BK*XSQ+CK*X+YOLD(K)

  600 K1=K1+1
  IF (K1.LE.NNEW) GO TO 300 ! Jump up

END SUBROUTINE SPLINE

!------------------------------------------------------------------------
!
! End subroutine declarations.
!
!------------------------------------------------------------------------

END MODULE MODULE_CU_BMJ