                              SUBROUTINE CUCNVC
C     ******************************************************************
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .     
C SUBPROGRAM:    CUCNVC      CONVECTIVE PRECIPITATION PARAMETERIZATION
C   PRGRMMR: JANJIC          ORG: W/NP2      DATE: 93-11-02       
C     
C ABSTRACT:
C     CUCNVC CALCULATES THE SUB-GRID SCALE CONVECTION INCLUDING 
C     DEEP AND SHALLOW CONVECTIVE CLOUDS FOLLOWING THE SCHEME DESCRIBED
C     BY JANJIC (1994) BUT WITH SIGNIFICANT MODIFICATIONS.  IN ADDITION,
C     THE LATENT HEAT RELEASE AND MOISTURE CHANGE DUE TO PRECIPITATING
C     AND NON-PRECIPITATING CLOUDS ARE COMPUTED.
C     
C     
C PROGRAM HISTORY LOG:
C   87-09-??  JANJIC     - ORIGINATOR
C   90-11-21  JANJIC     - TWO SETS OF DSP PROFILES (FAST AND SLOW)
C                          REPLACE THE ORIGINAL ONE SET
C   95-03-25  BLACK      - CONVERSION FROM 1-D TO 2-D IN HORIZONTAL
C   96-03-28  BLACK      - ADDED EXTERNAL EDGE
C   98-11-02  BLACK      - MODIFIED FOR DISTRIBUTED MEMORY
C     
C USAGE: CALL CUCNVC FROM MAIN PROGRAM EBU
C
C   INPUT ARGUMENT LIST:
C     NONE     
C  
C   OUTPUT ARGUMENT LIST: 
C     NONE
C     
C   OUTPUT FILES:
C     NONE
C     
C   SUBPROGRAMS CALLED:
C  
C     UNIQUE:
C        TTBLEX
C  
C     LIBRARY:
C        NONE
C  
C   COMMON BLOCKS: CTLBLK
C                  LOOPS
C                  MASKS
C                  PHYS
C                  VRBLS
C                  CNVCLD
C                  PVRBLS
C                  ACMCLH
C                  INDX
C   
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 90
C   MACHINE : IBM SP
C$$$  
C     ******************************************************************
C     *                                                                *
C     *  REFERENCES:                                                   *
C     *                                                                *
C     *  JANJIC, Z.I., 1994:  THE STEP-MOUNTAIN ETA COORDINATE MODEL:  *
C     *    FURTHER DEVELOPMENTS OF THE CONVECTION, VISCOUS SUBLAYER    *
C     *    AND TURBULENCE CLOSURE SCHEMES.  MONTHLY WEATHER REVIEW,    *
C     *    VOL. 122, 927-945.                                          *
C     *                                                                *
C     ******************************************************************
C *** WARNING: THIS SUBROUTINE WILL NOT WORK IF LM.LT.12;
C              MUST BE CALLED IN THE SAME STEP WITH PROFQ2 BECAUSE PROFQ
C              DEFINES APE;
C----------------------------------------------------------------------
      INCLUDE "cuparm"
C----------------------------------------------------------------------
      INCLUDE "parmeta"
      INCLUDE "parm.tbl"
      INCLUDE "mpp.h"
C----------------------------------------------------------------------
                             P A R A M E T E R
CVVVVVVVVVV INSTABILITY FOR TOO LARGE LSH VVVVVVVVVVVVVVVVVVVVVVVVVVVVV
     & (KSMUD=0,NROW= 0)
C    & (KSMUD=0,NROW= 5)
CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
                             P A R A M E T E R
     & (IMJM=IM*JM-JM/2,JAM=6+2*(JM-10)
     &, IMJM_LOC=IDIM2*JDIM2
CVVVVVVVVVV INSTABILITY FOR TOO LARGE LSH VVVVVVVVVVVVVVVVVVVVVVVVVVVVV
C    &, LP1=LM+1,LM1=LM-1,LNO=1,LSH=LM/3-1,LSHU=LM/2-1,LQM=LM/5,KBM=3
C    &, LP1=LM+1,LM1=LM-1,LNO=1,LSH=LM/3  ,LSHU=LM/2-1,LQM=LM/5,KBM=3
C    &, LP1=LM+1,LM1=LM-1,LNO=3,LSH=LM/3  ,LSHU=LM/2-1,LQM=LM/5,KBM=3
C    &, LP1=LM+1,LM1=LM-1,LNO=2,LSH=LM/3-1,LSHU=LM/2-1,LQM=LM/5,KBM=3
C    &, LP1=LM+1,LM1=LM-1,LNO=2,LSH=LM/3-2,LSHU=LM/2-1,LQM=LM/5,KBM=3
     &, LP1=LM+1,LM1=LM-1)
CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
C----------------------------------------------------------------------
                             L O G I C A L
     & RUN,FIRST,RESTRT,SIGMA
C----------------------------------------------------------------------
      INCLUDE "CTLBLK.comm"
C-----------------------------------------------------------------------
      INCLUDE "LOOPS.comm"
C-----------------------------------------------------------------------
      INCLUDE "MASKS.comm"
C-----------------------------------------------------------------------
      INCLUDE "PHYS.comm"
C-----------------------------------------------------------------------
      INCLUDE "VRBLS.comm"
C-----------------------------------------------------------------------
      INCLUDE "CNVCLD.comm"
C-----------------------------------------------------------------------
      INCLUDE "PVRBLS.comm"
C-----------------------------------------------------------------------
      INCLUDE "ACMCLH.comm"
C-----------------------------------------------------------------------
      INCLUDE "INDX.comm"
CYL---------------------------------------------------------------------
!      INCLUDE "PPTASM.comm"
CYL---------------------------------------------------------------------
                             D I M E N S I O N
     & TREFK (LM),QREFK (LM),PK    (LM),APEK  (LM),TK    (LM)
     &,THSK  (LM),PSK   (LM),APESK (LM),QK    (LM),THERK (LM)
     &,THVREF(LM),THEVRF(LM),THVMOD(LM),DIFT  (LM),DIFQ  (LM)
     &,QSATK (LM),FPK   (LM)
     &,NTOPD (LM),NBOTD (LM),NTOPS (LM),NBOTS (LM)
     &,NDPTHD(LM),NDPTHS(LM)
C
                             D I M E N S I O N
     & LTOP  (IDIM1:IDIM2,JDIM1:JDIM2),LBOT  (IDIM1:IDIM2,JDIM1:JDIM2)
     &,PTOP  (IDIM1:IDIM2,JDIM1:JDIM2),PBOT  (IDIM1:IDIM2,JDIM1:JDIM2)
     &,IPTB  (IDIM1:IDIM2,JDIM1:JDIM2),ITHTB (IDIM1:IDIM2,JDIM1:JDIM2)
     &,PDSL  (IDIM1:IDIM2,JDIM1:JDIM2),APEBT (IDIM1:IDIM2,JDIM1:JDIM2)
     &,TBT   (IDIM1:IDIM2,JDIM1:JDIM2),Q2BT  (IDIM1:IDIM2,JDIM1:JDIM2)
     &,QQ    (IDIM1:IDIM2,JDIM1:JDIM2),PP    (IDIM1:IDIM2,JDIM1:JDIM2)
     &,PSP   (IDIM1:IDIM2,JDIM1:JDIM2),THBT  (IDIM1:IDIM2,JDIM1:JDIM2)
     &,THESP (IDIM1:IDIM2,JDIM1:JDIM2),P     (IDIM1:IDIM2,JDIM1:JDIM2)
     &,BTH   (IDIM1:IDIM2,JDIM1:JDIM2),STH   (IDIM1:IDIM2,JDIM1:JDIM2)
     &,T00   (IDIM1:IDIM2,JDIM1:JDIM2),T10   (IDIM1:IDIM2,JDIM1:JDIM2)
     &,T01   (IDIM1:IDIM2,JDIM1:JDIM2),T11   (IDIM1:IDIM2,JDIM1:JDIM2)
     &,WF1   (IDIM1:IDIM2,JDIM1:JDIM2),WF2   (IDIM1:IDIM2,JDIM1:JDIM2)
     &,WF3   (IDIM1:IDIM2,JDIM1:JDIM2),WF4   (IDIM1:IDIM2,JDIM1:JDIM2)
     &,PRECOL(IDIM1:IDIM2,JDIM1:JDIM2)
C
     &,IBUOY (IMJM_LOC),JBUOY (IMJM_LOC)
     &,IDEEP (IMJM_LOC),JDEEP (IMJM_LOC)
     &,ISHAL (IMJM_LOC),JSHAL (IMJM_LOC)
     &,ILRES (IMJM_LOC),JLRES (IMJM_LOC)
     &,IHRES (IMJM_LOC),JHRES (IMJM_LOC)
C
     &,APE   (IDIM1:IDIM2,JDIM1:JDIM2,LM)
     &,TREF  (IDIM1:IDIM2,JDIM1:JDIM2,LM)
     &,TMOD  (IDIM1:IDIM2,JDIM1:JDIM2,LM)
     &,QMOD  (IDIM1:IDIM2,JDIM1:JDIM2,LM)
C
                             D I M E N S I O N
     & DSPB  (IDIM1:IDIM2,JDIM1:JDIM2),DSP0  (IDIM1:IDIM2,JDIM1:JDIM2)
     &,DSPT  (IDIM1:IDIM2,JDIM1:JDIM2)
     &,TL    (IDIM1:IDIM2,JDIM1:JDIM2),QL    (IDIM1:IDIM2,JDIM1:JDIM2)
     &,TNE   (IDIM1:IDIM2,JDIM1:JDIM2),TSE   (IDIM1:IDIM2,JDIM1:JDIM2)
     &,QNE   (IDIM1:IDIM2,JDIM1:JDIM2),QSE   (IDIM1:IDIM2,JDIM1:JDIM2)
C-----------------------------------------------------------------------
      CALL ZERO2(DSP0)
      CALL ZERO2(DSPB)
      CALL ZERO2(DSPT)
      CALL ZERO2(PSP)
C----------------------------------------------------------------------
      AVCNVC=AVCNVC+1.
      ACUTIM=ACUTIM+1.
      DTCNVC=NCNVC*DT
      RDTCNVC=1./DTCNVC
      TAUK=DTCNVC/TREL
c possible future change for shallow convection
c      TAUKSC=DTCNVC/(5.*TREL)   
C-----------------------------------------------------------------------
C************************** DIAGNOSTICS ********************************
C-----------------------------------------------------------------------
      DO L=1,LM
        NTOPD (L)=0
        NBOTD (L)=0
        NTOPS (L)=0
        NBOTS (L)=0
        NDPTHS(L)=0
        NDPTHD(L)=0
      ENDDO
C***********************************************************************
C---------------------------PREPARATIONS--------------------------------
C-----------------------------------------------------------------------
!$omp parallel do 
      DO 120 J=MYJS,MYJE
      DO 120 I=MYIS,MYIE
      LBOT (I,J)=LMH(I,J)
      THESP(I,J)=0.
      PDSL (I,J)=RES(I,J)*PD(I,J)
      PRECOL(I,J)=0.
      TREF(I,J,1)=T(I,J,1)
  120 CONTINUE
C-----------------------------------------------------------------------
C--------------PADDING SPECIFIC HUMIDITY IF TOO SMALL-------------------
C                  RESTORE APE TO SCRATCH ARRAY
C-----------------------------------------------------------------------
!$omp parallel do private(apests)
      DO 130 L=1,LM
      DO J=MYJS,MYJE
      DO I=MYIS,MYIE
        APESTS=PDSL(I,J)*AETA(L)+PT
        APE(I,J,L)=(1.E5/APESTS)**CAPA
        IF(Q(I,J,L).LT.EPSQ)Q(I,J,L)=HTM(I,J,L)*EPSQ
      ENDDO
      ENDDO
  130 CONTINUE
C-----------------------------------------------------------------------
C--------------SEARCH FOR MAXIMUM BUOYANCY LEVEL------------------------
C-----------------------------------------------------------------------
                             DO 170 KB=1,LM
C-----------------------------------------------------------------------
C--------------TRIAL MAXIMUM BUOYANCY LEVEL VARIABLES-------------------
C-----------------------------------------------------------------------
!$omp parallel do 
!$omp&  private(apesp,bq,bqs00k,bqs10k,iq,it,ittb,ittbk,iqtb,
!$omp&          lmhk,p00k,p01k,p10k,p11k,pkl,pp1,psfck,qbt,qq1,
!$omp&          sq,sqs00k,sqs10k,tpsp,tq,tth,tthbt,tthes)
      DO 155 J=MYJS2,MYJE2
      DO 150 I=MYIS1,MYIE1
C
      PKL=AETA(KB)*PDSL(I,J)+PT
      LMHK=LMH(I,J)
      PSFCK=AETA(LMHK)*PDSL(I,J)+PT
c now searching over a scaled depth in finding the parcel
c   with the max theta-e instead of the old 130 mb
c     IF(PKL.LT.PSFCK-PBM)GO TO 150
c     IF(KB.GT.LMHK)GO TO 150
      IF(KB.LE.LMHK .AND. PKL.GE.0.80*PSFCK) THEN
      QBT=Q(I,J,KB)
      TTHBT=T(I,J,KB)*APE(I,J,KB)
      TTH=(TTHBT-THL)*RDTH
      QQ1=TTH-AINT(TTH)
      ITTB=INT(TTH)+1
C--------------KEEPING INDICES WITHIN THE TABLE-------------------------
      IF(ITTB.LT.1)THEN
        ITTB=1
        QQ1=0.
      ENDIF
      IF(ITTB.GE.JTB)THEN
        ITTB=JTB-1
        QQ1=0.
      ENDIF
      CONTINUE
C--------------BASE AND SCALING FACTOR FOR SPEC. HUMIDITY---------------
      ITTBK=ITTB
      BQS00K=QS0(ITTBK)
      SQS00K=SQS(ITTBK)
      BQS10K=QS0(ITTBK+1)
      SQS10K=SQS(ITTBK+1)
C--------------SCALING SPEC. HUMIDITY & TABLE INDEX---------------------
      BQ  = (BQS10K-BQS00K)*QQ1+BQS00K
      SQ  = (SQS10K-SQS00K)*QQ1+SQS00K
      TQ  = (QBT-BQ)/SQ*RDQ
      PP1  =TQ - AINT(TQ)
      IQTB=INT(TQ)+1
C--------------KEEPING INDICES WITHIN THE TABLE-------------------------
      IF(IQTB.LT.1)THEN
        IQTB=1
        PP1=0.
      ENDIF
      IF(IQTB.GE.ITB)THEN
        IQTB=ITB-1
        PP1=0.
      ENDIF
C--------------SATURATION PRESSURE AT FOUR SURROUNDING TABLE PTS.-------
      IQ=IQTB
      IT=ITTB
      P00K=PTBL(IQ  ,IT  )
      P10K=PTBL(IQ+1,IT  )
      P01K=PTBL(IQ  ,IT+1)
      P11K=PTBL(IQ+1,IT+1)
C--------------SATURATION POINT VARIABLES AT THE BOTTOM-----------------
      TPSP=P00K+(P10K-P00K)*PP1+(P01K-P00K)*QQ1
     1        +(P00K-P10K-P01K+P11K)*PP1*QQ1
      APESP=(1.E5/TPSP)**CAPA
      TTHES=TTHBT*EXP(ELOCP*QBT*APESP/TTHBT)
C--------------CHECK FOR MAXIMUM BUOYANCY-------------------------------
       IF(TTHES.GT.THESP(I,J))THEN
        PSP  (I,J)=TPSP
        THBT (I,J)=TTHBT
        THESP(I,J)=TTHES
       ENDIF
      ENDIF
  150 CONTINUE
  155 CONTINUE
C-----------------------------------------------------------------------
  170 CONTINUE
C-----------------------------------------------------------------------
C---------CHOOSE CLOUD BASE AS MODEL LEVEL JUST BELOW PSP--------------
C-----------------------------------------------------------------------
      DO 240 L=1,LM1
      AETAL=AETA(L)
C
!$omp parallel do 
      DO J=MYJS2,MYJE2
      DO I=MYIS,MYIE
        P(I,J)=PDSL(I,J)*AETAL+PT
        IF(P(I,J).LT.PSP(I,J).AND.P(I,J).GE.PQM)LBOT(I,J)=L+1
      ENDDO
      ENDDO
C
  240 CONTINUE
C***
C*** WARNING: LBOT MUST NOT BE GT LMH(I,J)-1 IN SHALLOW CONVECTION
C*** MAKE SURE CLOUD BASE IS AT LEAST PONE ABOVE THE SURFACE
C***
!$omp parallel do private(lmhij,psfck)
      DO 250 J=MYJS2,MYJE2
      DO 250 I=MYIS,MYIE
      LMHIJ=LMH(I,J)
      PBOT(I,J)=AETA(LBOT(I,J))*PDSL(I,J)+PT
      PSFCK=AETA(LMHIJ)*PDSL(I,J)+PT
      IF(PBOT(I,J).GE.PSFCK-PONE.OR.LBOT(I,J).GE.LMHIJ)THEN
C
        DO L=1,LMHIJ-1
        P(I,J)=AETA(L)*PDSL(I,J)+PT
        IF(P(I,J).LT.PSFCK-PONE)LBOT(I,J)=L
        ENDDO
C
        PBOT(I,J)=AETA(LBOT(I,J))*PDSL(I,J)+PT
      ENDIF 
  250 CONTINUE
C--------------CLOUD TOP COMPUTATION------------------------------------
!$omp parallel do 
      DO J=MYJS,MYJE
      DO I=MYIS,MYIE
        LTOP(I,J)=LBOT(I,J)
        PTOP(I,J)=PBOT(I,J)
      ENDDO
      ENDDO
C-----------------------------------------------------------------------
!$omp parallel do 
!$omp&  private(ihres,ilres,iptb,ithtb,jhres,jlres,
!$omp&          knumh,knuml,pp,presk,qq)
C
      DO 290 L=LM,1,-1
C
C--------------SCALING PRESSURE & TT TABLE INDEX------------------------
      KNUML=0
      KNUMH=0
C
      DO 270 J=MYJS2,MYJE2
      DO 270 I=MYIS1,MYIE1
      PRESK=PDSL(I,J)*AETA(L)+PT
      IF(PRESK.LT.PLQ)THEN
        KNUML=KNUML+1
        ILRES(KNUML)=I
        JLRES(KNUML)=J
      ELSE
        KNUMH=KNUMH+1
        IHRES(KNUMH)=I
        JHRES(KNUMH)=J
      ENDIF
  270  CONTINUE
C***
C***  COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE<PL
C**
      IF(KNUML.GT.0)THEN
        CALL TTBLEX(TREF(IDIM1,JDIM1,L),TTBL,ITB,JTB,KNUML
     1,             ILRES,JLRES,PDSL,AETA(L),HTM(IDIM1,JDIM1,L)
     2,             PT,PL,QQ(IDIM1,JDIM1),PP(IDIM1,JDIM1)
     3,             RDP,THE0,STHE,RDTHE
     4,             THESP(IDIM1,JDIM1),IPTB(IDIM1,JDIM1)
     5,             ITHTB(IDIM1,JDIM1))
      ENDIF
C***
C***  COMPUTE PARCEL TEMPERATURE ALONG MOIST ADIABAT FOR PRESSURE>PL
C**
      IF(KNUMH.GT.0)THEN
        CALL TTBLEX(TREF(IDIM1,JDIM1,L),TTBLQ,ITBQ,JTBQ,KNUMH
     1,             IHRES,JHRES,PDSL,AETA(L),HTM(IDIM1,JDIM1,L)
     2,             PT,PLQ,QQ(IDIM1,JDIM1),PP(IDIM1,JDIM1)
     3,             RDPQ,THE0Q,STHEQ,RDTHEQ
     4,             THESP(IDIM1,JDIM1),IPTB(IDIM1,JDIM1)
     5,             ITHTB(IDIM1,JDIM1))
      ENDIF
  290  CONTINUE
C--------------BUOYANCY CHECK-------------------------------------------
      DO 295 L=LM,1,-1
!$omp parallel do 
      DO J=MYJS2,MYJE2
      DO I=MYIS1,MYIE1
        IF(TREF(I,J,L).GT.T(I,J,L)-DTTOP)LTOP(I,J)=L
      ENDDO
      ENDDO
C-----------------------------------------------------------------------
  295 CONTINUE
C-----------------CLOUD TOP PRESSURE------------------------------------
!$omp parallel do 
      DO J=MYJS2,MYJE2
      DO I=MYIS1,MYIE1
        PTOP(I,J)=AETA(LTOP(I,J))*PDSL(I,J)+PT
      ENDDO
      ENDDO
C-----------------------------------------------------------------------
C--------------DEFINE AND SMOOTH DSPS AND CLDEFI------------------------
C ************ UNIFIED OR SEPARATE LAND/SEA CONV ***********************
C-----------------------------------------------------------------------
C
      IF(UNIS)THEN
!$omp parallel do private(efi)
        DO J=MYJS,MYJE
        DO I=MYIS,MYIE
          EFI=CLDEFI(I,J)
          DSPB(I,J)=(EFI-EFIMN)*SLOPBS+DSPBSS
          DSP0(I,J)=(EFI-EFIMN)*SLOP0S+DSP0SS
          DSPT(I,J)=(EFI-EFIMN)*SLOPTS+DSPTSS
        ENDDO
        ENDDO
      ELSEIF(.NOT.UNIL)THEN
!$omp parallel do private(efi)
        DO J=MYJS,MYJE
        DO I=MYIS,MYIE
          EFI=CLDEFI(I,J)
          DSPB(I,J)=((EFI-EFIMN)*SLOPBS+DSPBSS)*SM(I,J)
     1             +((EFI-EFIMN)*SLOPBL+DSPBSL)*(1.-SM(I,J))
          DSP0(I,J)=((EFI-EFIMN)*SLOP0S+DSP0SS)*SM(I,J)
     1             +((EFI-EFIMN)*SLOP0L+DSP0SL)*(1.-SM(I,J))
          DSPT(I,J)=((EFI-EFIMN)*SLOPTS+DSPTSS)*SM(I,J)
     1             +((EFI-EFIMN)*SLOPTL+DSPTSL)*(1.-SM(I,J))
        ENDDO
        ENDDO
      ELSE
!$omp parallel do private(efi)
        DO J=MYJS,MYJE
        DO I=MYIS,MYIE
          EFI=CLDEFI(I,J)
          DSPB(I,J)=((EFI-EFIMN)*SLOPBL+DSPBSL)
          DSP0(I,J)=((EFI-EFIMN)*SLOP0L+DSP0SL)
          DSPT(I,J)=((EFI-EFIMN)*SLOPTL+DSPTSL) 
        ENDDO
        ENDDO
      ENDIF
C
C--------------EXTENDING SEA STRUCTURES INLAND ALONG COASTLINE----------
      IF(NROW.GT.0.AND..NOT.UNIS.AND..NOT.UNIL)THEN
!$omp parallel do 
        DO J=MYJS,MYJE
        DO I=MYIS,MYIE
          WF1(I,J)=0.
          WF2(I,J)=0.
          WF3(I,J)=0.
          WF4(I,J)=0.
        ENDDO
        ENDDO
C
        KROW=NROW
C
        DO 350 KOUNT=1,KROW
!$omp parallel do 
        DO 345 J=MYJS2,MYJE2
        DO 345 I=MYIS1,MYIE1
        WF1(I,J)=(DSPB(I+IHW(J),J-1)+DSPB(I+IHE(J),J-1)
     1           +DSPB(I+IHW(J),J+1)+DSPB(I+IHE(J),J+1)+4.*DSPB(I,J))
     2           *HBM2(I,J)*0.125
        WF2(I,J)=(DSP0(I+IHW(J),J-1)+DSP0(I+IHE(J),J-1)
     1           +DSP0(I+IHW(J),J+1)+DSP0(I+IHE(J),J+1)+4.*DSP0(I,J))
     2           *HBM2(I,J)*0.125
        WF3(I,J)=(DSPT(I+IHW(J),J-1)+DSPT(I+IHE(J),J-1)
     1           +DSPT(I+IHW(J),J+1)+DSPT(I+IHE(J),J+1)+4.*DSPT(I,J))
     2           *HBM2(I,J)*0.125
        WF4(I,J)=(CLDEFI(I+IHW(J),J-1)+CLDEFI(I+IHE(J),J-1)
     1     +CLDEFI(I+IHW(J),J+1)+CLDEFI(I+IHE(J),J+1)+4.*CLDEFI(I,J))
     2     *HBM2(I,J)*0.125
  345   CONTINUE 
C
!$omp parallel do private(rsmk,smk)
        DO J=MYJS,MYJE
        DO I=MYIS,MYIE
          SMK =SM(I,J)
          RSMK=1.-SMK
          DSPB  (I,J)=DSPB  (I,J)*SMK+WF1(I,J)*RSMK
          DSP0  (I,J)=DSP0  (I,J)*SMK+WF2(I,J)*RSMK
          DSPT  (I,J)=DSPT  (I,J)*SMK+WF3(I,J)*RSMK 
          CLDEFI(I,J)=CLDEFI(I,J)*SMK+WF4(I,J)*RSMK
        ENDDO
        ENDDO
C
  350   CONTINUE
C-----------------------------------------------------------------------
      ENDIF
C--------------INITIALIZE CHANGES OF T AND Q DUE TO CONVECTION----------
!$omp parallel do 
      DO 360 L=1,LM
      DO J=MYJS,MYJE
      DO I=MYIS,MYIE
        TMOD(I,J,L)=0.
        QMOD(I,J,L)=0.
      ENDDO
      ENDDO
  360 CONTINUE
C--------------CLEAN UP AND GATHER DEEP CONVECTION POINTS---------------
!$omp parallel do 
      DO 380 J=MYJS2,MYJE2
      DO 380 I=MYIS1,MYIE1
      IF(LTOP(I,J).GE.LBOT(I,J))THEN
        LBOT(I,J)=0
        LTOP(I,J)=LBOT(I,J)
        PTOP(I,J)=PBOT(I,J)
      ENDIF
      IF(HBM2(I,J).LT.0.90)THEN
        LBOT(I,J)=0
        LTOP(I,J)=LBOT(I,J)
        PTOP(I,J)=PBOT(I,J)
      ENDIF
      IF(PTOP(I,J).GT.PBOT(I,J)-PNO.OR.LTOP(I,J).GT.LBOT(I,J)-2)
     1 CLDEFI(I,J)=AVGEFI*SM(I,J)+STEFI*(1.-SM(I,J))
  380 CONTINUE
C
      KHDEEP=0
      PSHNEW=20000.
      DO J=MYJS2,MYJE2
       DO I=MYIS1,MYIE1
        PSFCIJ=PD(I,J)+PT
c depth of cloud required to make the point a deep convection point
c  is now a scaled value of the PSFC instead of 290 mb everywhere 
        DEPMIN=PSHNEW*PSFCIJ*1.E-5
        DEPTH=PBOT(I,J)-PTOP(I,J)
        IF(DEPTH .GE. DEPMIN) THEN
c        IF(PTOP(I,J).LT.PBOT(I,J)-PSH)THEN
          KHDEEP=KHDEEP+1
          IDEEP(KHDEEP)=I
          JDEEP(KHDEEP)=J
         ENDIF
       ENDDO
      ENDDO
C************* HORIZONTAL LOOP FOR DEEP CONVECTION *********************
!$omp parallel do
!$omp&  private(apek,apekl,apekxx,apekxy,apesk,avrgt,avtgtl,dentpy,
!$omp&          depmin,depth,depwl,dhdt,difq,difql,dift,diftl,drheat,
!$omp&          drhdp,dsp,dsp0k,dspbk,dsptk,dthem,efi,fefi,hcorr,
!$omp&          i,j,l0,l0m1,lb,lbm1,lbtk,lcor,lqm,lshu,ltp1,ltp2,
!$omp&          ltpk,ltsh,pbtk,pk,pk0,pkb,pkl,pkt,preck,psfcij,psk,
!$omp&          pthrs,ptpk,qk,qkl,qrefk,qsatk,rdp0t,rhh,rhl,rhmax,
!$omp&          sumde,sumdp,therk,therkx,therky,thsk,thskl,tk,tkl,
!$omp&          trefk,trefkx,tskl)
C***********************************************************************
      DO 600 N=1,KHDEEP
C***********************************************************************
      I=IDEEP(N)
      J=JDEEP(N)
      PSFCIJ=PD(I,J)+PT
      LTPK=LTOP(I,J)
      LBTK=LBOT(I,J)
CDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCD
CDCDCDCDCDCDC  DEEP CONVECTION   DCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCD
CDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCD
      LB   =LBTK
      EFI  =CLDEFI(I,J)
      DSPBK=DSPB  (I,J)
      DSP0K=DSP0  (I,J)
      DSPTK=DSPT  (I,J)
C--------------INITIALIZE VARIABLES IN THE CONVECTIVE COLUMN------------
C***
C***  ONE SHOULD NOTE THAT THE VALUES ASSIGNED TO THE ARRAY TREFK
C***  IN THE 410 LOOP ARE REALLY ONLY RELEVANT IN ANCHORING THE
C***  REFERENCE TEMPERATURE PROFILE AT LEVEL LB.  WHEN BUILDING THE
C***  REFERENCE PROFILE FROM CLOUD BASE, THEN ASSIGNING THE
C***  AMBIENT TEMPERATURE TO TREFK IS ACCEPTABLE.  HOWEVER, WHEN
C***  BUILDING THE REFERENCE PROFILE FROM SOME OTHER LEVEL (SUCH AS
C***  ONE LEVEL ABOVE THE GROUND), THEN TREFK SHOULD BE FILLED WITH
C***  THE TEMPERATURES IN TREF(I,J,L) WHICH ARE THE TEMPERATURES OF
C***  THE MOIST ADIABAT THROUGH CLOUD BASE.  BY THE TIME THE LINE 
C***  NUMBERED 450 HAS BEEN REACHED, TREFK ACTUALLY DOES HOLD THE
C***  REFERENCE TEMPERATURE PROFILE.
C***
      DO 410 L=1,LM
      DIFT  (L)=0.
      DIFQ  (L)=0.
      TKL      =T(I,J,L)
      TK    (L)=TKL
      TREFK (L)=TKL
      QKL      =Q(I,J,L)
      QK    (L)=QKL
      QREFK (L)=QKL
      PKL      =AETA(L)*PDSL(I,J)+PT
      PK    (L)=PKL
      PSK   (L)=PKL
      APEKL    =APE(I,J,L)
      APEK  (L)=APEKL
      THERK (L)=TREF(I,J,L)*APEKL
C 410 THVMOD(L)=(QKL*0.608+1.)*TKL*APEKL
  410     CONTINUE
C--------------DEEP CONVECTION REFERENCE TEMPERATURE PROFILE------------
      LTP1=LTPK+1
      LBM1=LB-1
      PKB=PK(LB)
      PKT=PK(LTPK)
C--------------TEMPERATURE REFERENCE PROFILE BELOW FREEZING LEVEL-------
      L0=LB
      PK0=PK(LB)
      TREFKX=TREFK(LB)
      THERKX=THERK(LB)
      APEKXX=APEK(LB)
      THERKY=THERK(LBM1)
      APEKXY=APEK(LBM1)
C
      DO 420 L=LBM1,LTPK,-1
      IF(T(I,J,L+1).LT.TFRZ)GO TO 430
      STABDL=STABD
      TREFKX=((THERKY-THERKX)*STABDL
     1        +TREFKX*APEKXX)/APEKXY
      TREFK(L)=TREFKX
      APEKXX=APEKXY
      THERKX=THERKY
      APEKXY=APEK(L-1)
      THERKY=THERK(L-1)
      L0=L
      PK0=PK(L0)
  420 CONTINUE
C--------------FREEZING LEVEL AT OR ABOVE THE CLOUD TOP-----------------
      L0M1=L0-1
      GO TO 450
C--------------TEMPERATURE REFERENCE PROFILE ABOVE FREEZING LEVEL-------
  430 L0M1=L0-1
      RDP0T=1./(PK0-PKT)
      DTHEM=THERK(L0)-TREFK(L0)*APEK(L0)
C
      DO L=LTPK,L0M1
        TREFK(L)=(THERK(L)-(PK(L)-PKT)*DTHEM*RDP0T)/APEK(L)
      ENDDO
C-----------------------------------------------------------------------
C--------------DEEP CONVECTION REFERENCE HUMIDITY PROFILE---------------
C-----------------------------------------------------------------------
c  reference profile had been too dry in the case where the cloud
c   base was close to the freezing level - this is now changed
  450 DEPTH=PFRZ*PSFCIJ*1.E-5
      DEPWL=PKB-PK0
      DO 460 L=LTPK,LB
C-----------------------------------------------------------------------
C--------------SATURATION PRESSURE DIFFERENCE---------------------------
C-----------------------------------------------------------------------
c     IF(PKB-PK0.GE.PFRZ)THEN
      IF(DEPWL .GE. DEPTH) THEN
        IF(L.LT.L0)THEN
          DSP=((PK0-PK(L))*DSPTK+(PK(L)-PKT)*DSP0K)/(PK0-PKT)
        ELSE
          DSP=((PKB-PK(L))*DSP0K+(PK(L)-PK0)*DSPBK)/(PKB-PK0)
        ENDIF
      ELSE
c       DSP=DSPC
        DSP=DSP0K
        IF(L.LT.L0)
     1 DSP=( (PK0-PK(L))*DSPTK+(PK(L)-PKT)*DSP0K)/(PK0-PKT)
      ENDIF

C--------------HUMIDITY PROFILE-----------------------------------------
      IF(PK(L).GT.PQM)THEN
        PSK(L)=PK(L)+DSP
        APESK(L)=(1.E5/PSK(L))**CAPA
        THSK(L)=TREFK(L)*APEK(L)
        QREFK(L)=PQ0/PSK(L)*EXP(A2*(THSK(L)-A3*APESK(L))
     1                            /(THSK(L)-A4*APESK(L)))
      ELSE
        QREFK(L)=Q(I,J,L)
      ENDIF
  460 CONTINUE
C--------------ENTHALPY CONSERVATION INTEGRAL--------------------------
                             DO 520 ITER=1,2
C-----------------------------------------------------------------------
      SUMDE=0.
      SUMDP=0.
C
      DO L=LTPK,LB
        SUMDE=((TK(L)-TREFK(L))*CP+(QK(L)-QREFK(L))*ELWV)*DETA(L)
     1        +SUMDE
        SUMDP=SUMDP+DETA(L)
      ENDDO
C
      HCORR=SUMDE/(SUMDP-DETA(LTPK))
      LCOR=LTPK+1
C-----------------------FIND LQM----------------------------------------
      DO L=1,LB
        IF(PK(L).LE.PQM)LQM=L
      ENDDO
C--------------ABOVE LQM CORRECT TEMPERATURE ONLY-----------------------
      IF(LCOR.LE.LQM)THEN
        DO L=LCOR,LQM
          TREFK(L)=TREFK(L)+HCORR*RCP
        ENDDO
        LCOR=LQM+1
      ENDIF
C--------------BELOW LQM CORRECT BOTH TEMPERATURE AND MOISTURE----------
C
      DO 510 L=LCOR,LB
      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))
     1                          /(THSKL-A4*APESK(L)))
  510 CONTINUE
C-----------------------------------------------------------------------
  520 CONTINUE
C--------------HEATING, MOISTENING, PRECIPITATION-----------------------
      DENTPY=0.
      AVRGT =0.
      PRECK =0.
C
      DO 530 L=LTPK,LB
      TKL    =TK(L)
      DIFTL  =(TREFK(L)-TKL  )*TAUK
      DIFQL  =(QREFK(L)-QK(L))*TAUK
      AVRGTL =(TKL+TKL+DIFTL)
      DENTPY =(DIFTL*CP+DIFQL*ELWV)*DETA(L)/AVRGTL+DENTPY
      AVRGT  =AVRGTL*DETA(L)+AVRGT
      PRECK  =DETA(L)*DIFTL+PRECK
      DIFT(L)=DIFTL
      DIFQ(L)=DIFQL
  530 CONTINUE
C
      DENTPY=DENTPY+DENTPY
      AVRGT =AVRGT/(SUMDP+SUMDP)
C-------------SWAP IF ENTROPY AND/OR PRECIP .LT. 0 ...------------------
      IF(DENTPY.LT.EPSNTP.OR.PRECK.LT.0.)THEN
        IF(OCT90)THEN
          CLDEFI(I,J)=EFIMN
        ELSE
          CLDEFI(I,J)=EFIMN*SM(I,J)+STEFI*(1.-SM(I,J))
        ENDIF
C
C--------------SEARCH FOR SHALLOW CLOUD TOP-----------------------------
        LBTK=LBOT(I,J)
        LTSH=LBTK
        LBM1=LBTK-1
        PBTK=PK(LBTK)
c use new threshold for cloud depth
c        PTPK=PBTK-PSH
        PSFCIJ=PD(I,J)+PT
        DEPMIN=PSHNEW*PSFCIJ*1.E-5
        PTPK=PBTK-DEPMIN
C-------------CLOUD TOP IS THE LEVEL JUST BELOW PBTK-PSH----------------
        DO L=1,LM
          IF(PK(L).LE.PTPK)LTPK=L+1
        ENDDO
        PTPK=PK(LTPK)
C--------------HIGHEST LEVEL ALLOWED IS LEVEL JUST BELOW PSHU-----------
        IF(PTPK.LE.PSHU)THEN
          DO L=1,LM
            IF(PK(L).LE.PSHU)LSHU=L+1
          ENDDO
          LTPK=LSHU
          PTPK=PK(LTPK)
        ENDIF
C
        LTP1=LTPK+1
        LTP2=LTPK+2
C-----------------------------------------------------------------------
        DO L=LTPK,LBTK
          QSATK(L)=PQ0/PK(L)*EXP(A2*(TK(L)-A3)/(TK(L)-A4))
        ENDDO
C-----------------------------------------------------------------------
        RHH=QK(LTPK)/QSATK(LTPK)
        RHMAX=0.
C
        DO 570 L=LTP1,LBM1
        RHL=QK(L)/QSATK(L)
        DRHDP=(RHH-RHL)/(PK(L-1)-PK(L))
        IF(DRHDP.GT.RHMAX)THEN
          LTSH=L-1
          RHMAX=DRHDP
        ENDIF
        RHH=RHL
  570   CONTINUE
C
        LTOP(I,J)=LTSH
C---------------CLOUD MUST BE AT LEAST TWO LAYERS THICK-----------------
        IF(LBOT(I,J)-LTOP(I,J).LT.2)LTOP(I,J)=LBOT(I,J)-2
C
        PTOP(I,J)=PK(LTOP(I,J))
        GO TO 600
      ENDIF
C-----------------------------------------------------------------------
C--------------... DEEP CONVECTION OTHERWISE----------------------------
C-----------------------------------------------------------------------
      DRHEAT=(PRECK*SM(I,J)+AMAX1(EPSP,PRECK)*(1.-SM(I,J)))
     1       *CP/AVRGT
      EFI=EFIFC*DENTPY/DRHEAT
C
C************** UNIFIED OR SEPARATE LAND/SEA CONV. **************
C
      IF(OCT90)THEN
        IF(UNIS)THEN
          EFI=CLDEFI(I,J)*FCB+EFI*FCC
        ELSEIF(.NOT.UNIL)THEN
          EFI=(CLDEFI(I,J)*FCB+EFI*FCC)*SM(I,J)+1.-SM(I,J)
        ELSE
          EFI=1.
        ENDIF
      ELSE
        EFI=CLDEFI(I,J)*FCB+EFI*FCC
      ENDIF
C
      IF(EFI.GT.1.   )    EFI=1.
      IF(EFI.LT.EFIMN)    EFI=EFIMN
      IF(PRECK.EQ.0.)     EFI=1.
      CLDEFI(I,J)=EFI
C
      FEFI=EFMNT+SLOPE*(EFI-EFIMN)
C     FEFI=AMAX1(EFI,EFMNT)
C
      PRECK=PRECK*FEFI
C
C--------------UPDATE PRECIPITATION, TEMPERATURE & MOISTURE-------------
C
      PRECOL(I,J)=PDSL(I,J)*PRECK*CPRLG
      PREC  (I,J)=PDSL(I,J)*PRECK*CPRLG+PREC  (I,J)
      CUPREC(I,J)=PDSL(I,J)*PRECK*CPRLG+CUPREC(I,J)
      ACPREC(I,J)=PDSL(I,J)*PRECK*CPRLG+ACPREC(I,J)
      CUPPT(I,J) =PDSL(I,J)*PRECK*CPRLG+CUPPT(I,J)
C
      DO L=LTPK,LB
        TMOD(I,J,L)=DIFT(L)*FEFI
        QMOD(I,J,L)=DIFQ(L)*FEFI
CYL
C        TLAT(I,J,L)=DIFT(L)*FEFI + TLAT(I,J,L)
CYL
      ENDDO
C
CDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCD
CDCDCDCDCDCDC          END OF DEEP CONVECTION            DCDCDCDCDCDCDCD
CDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCDCD
C
C-----------------------------------------------------------------------
  600 CONTINUE
C-----------------------------------------------------------------------
      NDEEP=0
C
      DO 620 J=MYJS2,MYJE2
      DO 620 I=MYIS,MYIE
      LTPK=LTOP(I,J)
      LBTK=LBOT(I,J)
      LB  =LMH(I,J)-1
      PSFCIJ=PD(I,J)+PT
      DEPMIN=PSHNEW*PSFCIJ*1.E-5
      IF(PTOP(I,J).LT.PBOT(I,J)-DEPMIN)THEN
        NDEEP=NDEEP+1
        NDEPTH=LB-LTPK
        NTOPD (LTPK  )=NTOPD (LTPK  )+1
        NBOTD (LB    )=NBOTD (LB    )+1
        NDPTHD(NDEPTH)=NDPTHD(NDEPTH)+1
      ENDIF
  620 CONTINUE
      NNEG=KHDEEP-NDEEP
C
C--------------GATHER SHALLOW CONVECTION POINTS-------------------------
C
      KHSHAL=0
      NDSTN =0
      NDSTP =0
C
      DO 630 J=MYJS2,MYJE2
      DO 630 I=MYIS,MYIE
      IF(PTOP(I,J).GT.PBOT(I,J)-PNO.OR.
     1   LTOP(I,J).GT.LBOT(I,J)-2)GO TO 630
      PSFCIJ=PD(I,J)+PT
      DEPMIN=PSHNEW*PSFCIJ*1.E-5
      IF(PTOP(I,J)+1..GE.PBOT(I,J)-DEPMIN)THEN
        KHSHAL=KHSHAL+1
        ISHAL(KHSHAL)=I
        JSHAL(KHSHAL)=J
      ENDIF
  630 CONTINUE
C
C************* HORIZONTAL LOOP FOR SHALLOW CONVECTION ******************
!$omp parallel do 
!$omp&  private(apek,apekl,apekxx,apekxy,bqk,bqs00k,bqs10k,den,dentpy,
!$omp&           dpkl,dpmix,dqref,dst,dstq,dtdeta,fpk,fptk,i,iq,it,j,
!$omp&           lbm1,lbtk,ltp1,ltpk,otsum,part1,part2,part3,pk,pkl,
!$omp&           pkxxxx,pkxxxy,potsum,ppk,psum,ptpk,pz0,qk,qkl,qnew,
!$omp&           qotsum,qqk,qrefk,qrfkl,qrftp,qsatk,qsum,rdpsum,rtbar,
!$omp&           smix,sqk,sqs00k,sqs10k,sumdp,sumdt,tcorr,thvmkl,
!$omp&           thvref,tk,tkl,tqk,trefk,trefkx,trfkl,tthk)
C
      DO 800 N=1,KHSHAL
C***********************************************************************
      I=ISHAL(N)
      J=JSHAL(N)
C-----------------------------------------------------------------------
CSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCS
CSCSCSCSCSCSC         SHALLOW CONVECTION          CSCSCSCSCSCSCSCSCSCSCS
CSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCS
C-----------------------------------------------------------------------
      PZ0=PD(I,J)+PT
      LLMH=LMH(I,J)
C
      DO 650 L=1,LLMH
      TKL      =T  (I,J,L)
      TK   (L) =TKL
      TREFK(L) =TKL
      QKL      =Q  (I,J,L)
      QK   (L) =QKL
      QREFK(L) =QKL
      PKL      =AETA(L)*PDSL(I,J)+PT
      PK   (L) =PKL
      QSATK(L) =PQ0/PK(L)*EXP(A2*(TK(L)-A3)/(TK(L)-A4))
      APEKL    =APE(I,J,L)
CVVVVVVVVVVVV CHOOSE THE PRESSURE FUNCTION VVVVVVVVVVVVVVVVVVVVVVVVVVVVV
C         FPK  (L) =ALOG(PKL)
C         FPK  (L) =PKL
C         FPK  (L) =-1./PKL
CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
      APEK (L) =APEKL
      THVMKL   =TKL*APEKL*(QKL*0.608+1.)
C     THVMOD(L)=THVMKL
      THVREF(L)=THVMKL
  650 CONTINUE
C--------------SHALLOW CLOUD TOP-----------------------------
      LBTK=LBOT(I,J)
      LBM1=LBTK-1
      PTPK=PTOP(I,J)
      LTP1=LTOP(I,J)
      LTPK=LTOP(I,J)-1
C-----------------------------------------------------------------------
      IF(PTOP(I,J).GT.PBOT(I,J)-PNO.OR.LTOP(I,J).GT.LBOT(I,J)-2)THEN
        LBOT(I,J)=0
        LTOP(I,J)=LBOT(I,J)
        PTOP(I,J)=PBOT(I,J)
        GO TO 800
      ENDIF
C--------------SCALING POTENTIAL TEMPERATURE & TABLE INDEX AT TOP-------
C
C     THTPK=T(I,J,LTP1)*APE(I,J,LTP1)
      THTPK=T(I,J,LTPK)*APE(I,J,LTPK)
C
      TTHK =(THTPK-THL)*RDTH
      QQK  =TTHK-AINT(TTHK)
      IT   =INT(TTHK)+1
      IF(IT.LT.1)THEN
        IT=1
        QQK=0.
      ENDIF
      IF(IT.GE.JTB)THEN
        IT=JTB-1
        QQK=0.
      ENDIF
C--------------BASE AND SCALING FACTOR FOR SPEC. HUMIDITY AT TOP--------
      BQS00K=QS0(IT)
      SQS00K=SQS(IT)
      BQS10K=QS0(IT+1)
      SQS10K=SQS(IT+1)
C--------------SCALING SPEC. HUMIDITY & TABLE INDEX AT TOP--------------
      BQK=(BQS10K-BQS00K)*QQK+BQS00K
      SQK=(SQS10K-SQS00K)*QQK+SQS00K
C
C     TQK=(Q(I,J,LTP1)-BQK)/SQK*RDQ
      TQK=(Q(I,J,LTPK)-BQK)/SQK*RDQ
C
      PPK=TQK-AINT(TQK)
      IQ =INT(TQK)+1
      IF(IQ.LT.1)THEN
        IQ=1
        PPK=0.
      ENDIF
      IF(IQ.GE.ITB)THEN
        IQ=ITB-1
        PPK=0.
      ENDIF
C--------------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  )
     1      -PTBL(IQ  ,IT+1)+PTBL(IQ+1,IT+1))*PPK*QQK
      PTPK=PTBL(IQ,IT)+PART1+PART2+PART3
C-----------------------------------------------------------------------
      DPMIX=PTPK-PSP(I,J)
      IF(ABS(DPMIX).LT.3000.)DPMIX=-3000.
C--------------TEMPERATURE PROFILE SLOPE--------------------------------
      SMIX=(THTPK-THBT(I,J))/DPMIX*STABS
C
      TREFKX=TREFK(LBTK+1)
      PKXXXX=PK(LBTK+1)
      PKXXXY=PK(LBTK)
      APEKXX=APEK(LBTK+1)
      APEKXY=APEK(LBTK)
C
      DO 670 L=LBTK,LTP1,-1
      TREFKX=((PKXXXY-PKXXXX)*SMIX
     1        +TREFKX*APEKXX)/APEKXY
      TREFK(L)=TREFKX
      APEKXX=APEKXY
      PKXXXX=PKXXXY
      APEKXY=APEK(L-1)
      PKXXXY=PK(L-1)
  670 CONTINUE
C--------------TEMPERATURE REFERENCE PROFILE CORRECTION-----------------
      SUMDT=0.
      SUMDP=0.
C
      DO L=LTP1,LBTK
        SUMDT=(TK(L)-TREFK(L))*DETA(L)+SUMDT
        SUMDP=SUMDP+DETA(L)
      ENDDO
C
      RDPSUM=1./SUMDP
      FPK(LBTK)=TREFK(LBTK)
C
      TCORR=SUMDT*RDPSUM
C
      DO L=LTP1,LBTK
C     TCORR=SUMDT/(SUMDP-DETA(LBTK))
C     DO L=LTP1,LBM1
C
        TRFKL   =TREFK(L)+TCORR
        TREFK(L)=TRFKL
        FPK  (L)=TRFKL
      ENDDO
C--------------HUMIDITY PROFILE EQUATIONS-------------------------------
      PSUM  =0.
      QSUM  =0.
      POTSUM=0.
      QOTSUM=0.
      OTSUM =0.
      DST   =0.
      FPTK  =FPK(LTP1)
C
      DO 700 L=LTP1,LBTK
      DPKL  =FPK(L)-FPTK
      PSUM  =DPKL *DETA(L)+PSUM
      QSUM  =QK(L)*DETA(L)+QSUM
      RTBAR =2./(TREFK(L)+TK(L))
      OTSUM =DETA(L)*RTBAR+OTSUM
      POTSUM=DPKL   *RTBAR*DETA(L)+POTSUM
      QOTSUM=QK(L)  *RTBAR*DETA(L)+QOTSUM
      DST   =(TREFK(L)-TK(L))*RTBAR*DETA(L)+DST
  700 CONTINUE
C
      PSUM  =PSUM*RDPSUM
      QSUM  =QSUM*RDPSUM
      ROTSUM=1./OTSUM
      POTSUM=POTSUM*ROTSUM
      QOTSUM=QOTSUM*ROTSUM
      DST   =DST   *ROTSUM*CP/ELWV
C--------------ENSURE POSITIVE ENTROPY CHANGE---------------------------
      IF(DST.GT.0.)THEN
C
C       DSTQ=DST*EPSUP
        LBOT(I,J)=0
        LTOP(I,J)=LBOT(I,J)
        PTOP(I,J)=PBOT(I,J)
cccc    NDSTP=NDSTP+1
        GO TO 800
C
      ELSE
        DSTQ=DST*EPSDN
      ENDIF
C--------------CHECK FOR ISOTHERMAL ATMOSPHERE--------------------------
      DEN=POTSUM-PSUM
C
      IF(-DEN/PSUM.LT.5.E-5)THEN
        LBOT(I,J)=0
        LTOP(I,J)=LBOT(I,J)
        PTOP(I,J)=PBOT(I,J)
        GO TO 800
C
C--------------SLOPE OF THE REFERENCE HUMIDITY PROFILE------------------
      ELSE
        DQREF=(QOTSUM-DSTQ-QSUM)/DEN
      ENDIF
C------------ HUMIDITY DOES NOT INCREASE WITH HEIGHT--------------------
      IF(DQREF.LT.0.)THEN
        LBOT(I,J)=0
        LTOP(I,J)=LBOT(I,J)
        PTOP(I,J)=PBOT(I,J)
        GO TO 800
      ENDIF
C--------------HUMIDITY AT THE CLOUD TOP--------------------------------
      QRFTP=QSUM-DQREF*PSUM
C--------------HUMIDITY PROFILE-----------------------------------------
C
      DO 720 L=LTP1,LBTK
      QRFKL=(FPK(L)-FPTK)*DQREF+QRFTP
C
C--------------SUPERSATURATION OR NEGATIVE Q NOT ALLOWED----------------
C
      QNEW =(QRFKL-QK(L))*TAUK+QK(L)
      IF(QNEW.GT.QSATK(L)*STRESH.OR.QNEW.LT.0.)THEN
        LBOT(I,J)=0
        LTOP(I,J)=LBOT(I,J)
        PTOP(I,J)=PBOT(I,J)
        GO TO 800
      ENDIF
C
C-----------------------------------------------------------------------
      THVREF(L)=TREFK(L)*APEK(L)*(QRFKL*0.608+1.)
      QREFK(L)=QRFKL
  720 CONTINUE
C--------------ELIMINATE IMPOSSIBLE SLOPES (BETTS, DTHETA/DQ)-----------
      DO 730 L=LTP1,LBTK
      DTDETA=(THVREF(L-1)-THVREF(L))/(AETA(L)-AETA(L-1))
      IF(DTDETA.LT.EPSTH)THEN
        LBOT(I,J)=0
        LTOP(I,J)=LBOT(I,J)
        PTOP(I,J)=PBOT(I,J)
        GO TO 800
      ENDIF
  730 CONTINUE
C*************************** DIAGNOSTICS ******************************
cccc  IF(DST.GT.0.)THEN
cccc    NDSTP=NDSTP+1
cccc  ELSE
cccc    NDSTN=NDSTN+1
cccc  ENDIF
C**********************************************************************
      DENTPY=0.
C
      DO L=LTP1,LBTK
        DENTPY=((TREFK(L)-TK(L))*CP+(QREFK(L)-QK(L))*ELWV)
     1         /(TK(L)+TREFK(L))*DETA(L)+DENTPY
      ENDDO
C
C-----------------------------------------------------------------------
C--------------RELAXATION TOWARDS REFERENCE PROFILES--------------------
C-----------------------------------------------------------------------
C
C
      DO 750 L=LTP1,LBTK
      TMOD(I,J,L)=(TREFK(L)-TK(L))*TAUK
      QMOD(I,J,L)=(QREFK(L)-QK(L))*TAUK
  750 CONTINUE
C-----------------------------------------------------------------------
CSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCS
CSCSCSCSCSCSC         END OF SHALLOW CONVECTION        SCSCSCSCSCSCSCSCS
CSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCSCS
C-----------------------------------------------------------------------
  800                        CONTINUE
C-----------------------------------------------------------------------
C
C************************** DIAGNOSTICS ********************************
      NSHAL=0
C
      DO 820 J=MYJS2,MYJE2
      DO 820 I=MYIS,MYIE
      LTPK=LTOP(I,J)
      LBTK=LBOT(I,J)
      PTPK=PTOP(I,J)
      PBTK=PBOT(I,J)
      IF(PTPK.GT.PBTK-PNO.OR.LTPK.GT.LBTK-2)GO TO 820
      PSFCIJ=PD(I,J)+PT
      DEPMIN=PSHNEW*PSFCIJ*1.E-5
C      IF(PTPK.GE.PBTK-PSH)THEN
      IF(PTPK.GE.PBTK-DEPMIN)THEN
        NSHAL=NSHAL+1
        NTOPS(LTPK)=NTOPS(LTPK)+1
        NBOTS(LBTK)=NBOTS(LBTK)+1
        NDEPTH=LBTK-LTPK
        NDPTHS(NDEPTH)=NDPTHS(NDEPTH)+1
      ENDIF
  820 CONTINUE
      NEGDS=KHSHAL-NSHAL
C***********************************************************************
C
C--------------SMOOTHING TEMPERATURE & HUMIDITY CORRECTIONS-------------
      IF(KSMUD.EQ.0)THEN
!$omp parallel do 
C
        DO 900 L=1,LM
C***
C***  UPDATE THE FUNDAMENTAL PROGNOSTIC ARRAYS
C***
        DO 830 J=MYJS,MYJE
        DO 830 I=MYIS,MYIE
        T(I,J,L)=T(I,J,L)+TMOD(I,J,L)
        Q(I,J,L)=Q(I,J,L)+QMOD(I,J,L)
C***
C***  ACCUMULATE LATENT HEATING DUE TO CONVECTION.
C***  SCALE BY THE RECIPROCAL OF THE PERIOD AT WHICH THIS ROUTINE
C***  IS CALLED.  THIS PERIOD IS THE CONVECTION TIMESTEP.
C***
        TCUCN(I,J,L)=TCUCN(I,J,L)+TMOD(I,J,L)*RDTCNVC
  830  CONTINUE
  900  CONTINUE
      ELSE
!$omp parallel do private(ql,qne,qse,tl,tne,tse)
C
        DO 910 L=1,LM
C
        CALL ZERO2(QL)
        CALL ZERO2(QNE)
        CALL ZERO2(QSE)
        CALL ZERO2(TL)
        CALL ZERO2(TNE)
        CALL ZERO2(TSE)
C-----------------------------------------------------------------------
        DO J=MYJS,MYJE
        DO I=MYIS,MYIE
          TL(I,J)=TMOD(I,J,L)
          QL(I,J)=QMOD(I,J,L)
        ENDDO
        ENDDO
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
        NSMUD=KSMUD
C
        DO 870 KS=1,NSMUD
C
        DO J=MYJS,MYJE1
        DO I=MYIS,MYIE
        TNE(I,J)=(TL(I+IHE(J),J+1)-TL(I,J))
     1           *HTM(I,J,L)*HTM(I+IHE(J),J+1,L)
        QNE(I,J)=(QL(I+IHE(J),J+1)-QL(I,J))
     1           *HTM(I,J,L)*HTM(I+IHE(J),J+1,L)
        ENDDO
        ENDDO
C
        DO J=MYJS1,MYJE
        DO I=MYIS,MYIE
        TSE(I,J)=(TL(I+IHE(J),J-1)-TL(I,J))
     1           *HTM(I+IHE(J),J-1,L)*HTM(I,J,L)
        QSE(I,J)=(QL(I+IHE(J),J-1)-QL(I,J))
     1           *HTM(I+IHE(J),J-1,L)*HTM(I,J,L)
        ENDDO
        ENDDO
C
        DO J=MYJS2,MYJE2
        DO I=MYIS,MYIE
        TL(I,J)=(TNE(I,J)-TNE(I+IHW(J),J-1)+TSE(I,J)-TSE(I+IHW(J),J+1))
     1          *HBM2(I,J)*0.125+TL(I,J)
        QL(I,J)=(QNE(I,J)-QNE(I+IHW(J),J-1)+QSE(I,J)-QSE(I+IHW(J),J+1))
     1          *HBM2(I,J)*0.125+QL(I,J)
        ENDDO
        ENDDO
C
  870   CONTINUE
C-----------------------------------------------------------------------
C***
C***  UPDATE THE FUNDAMENTAL PROGNOSTIC ARRAYS
C***
        DO J=MYJS,MYJE
        DO I=MYIS,MYIE
          T(I,J,L)=T(I,J,L)+TL(I,J)
          Q(I,J,L)=Q(I,J,L)+QL(I,J)
        ENDDO
        ENDDO
C***
C***  ACCUMULATE LATENT HEATING DUE TO CONVECTION.
C***  SCALE BY THE RECIPROCAL OF THE PERIOD AT WHICH THIS ROUTINE
C***  IS CALLED.  THIS PERIOD IS THE CONVECTION TIMESTEP.
C***
        DO J=MYJS,MYJE
        DO I=MYIS,MYIE
          TCUCN(I,J,L)=TCUCN(I,J,L)+TL(I,J)*RDTCNVC
        ENDDO
        ENDDO
C
  910   CONTINUE
C
      ENDIF
C--------------SAVE CLOUD TOP AND BOTTOM FOR RADIATION------------------
!$omp parallel do 
      DO J=MYJS,MYJE
      DO I=MYIS,MYIE
        HTOP(I,J)=MIN(FLOAT(LTOP(I,J)),HTOP(I,J))
        HBOT(I,J)=MAX(FLOAT(LBOT(I,J)),HBOT(I,J))
        CNVTOP(I,J)=MIN(FLOAT(LTOP(I,J)),CNVTOP(I,J))
        CNVBOT(I,J)=MAX(FLOAT(LBOT(I,J)),CNVBOT(I,J))
      ENDDO
      ENDDO
C-----------------------------------------------------------------------
C************************* DIAGNOSTICS *********************************
C
C     WRITE(LIST,950)NTSD,NSHAL,NDEEP,NNEG,NEGDS,NDSTN,NDSTP
C         DO 940 L=1,LM
C     WRITE(LIST,952)L
C     WRITE(LIST,954)NBOTS(L),NTOPS(L),NDPTHS(L)
C    1,              NBOTD(L),NTOPD(L),NDPTHD(L)
C 940 CONTINUE
C 950 FORMAT(' NTSD=',I3,I8,' SHALLOW, ',I4,' DEEP, ',
C    1 I4,' NEG., ',I4,' NEG. SHALL.,',I4,' DST.LT.0, ',I4,' DST.GT.0')
C 952 FORMAT(' LAYER (FROM TOP),',I2)
C 954 FORMAT('     NBOTS=',I4,'     NTOPS=',I4,'     NDPTHS=',I4,
C    1       '     NBOTD=',I4,'     NTOPD=',I4,'     NDPTHD=',I4)
C***********************************************************************
                             RETURN
                             END
