    SUBROUTINE DIVHOASTQL
!     ******************************************************************
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!                .      .    .
! SUBPROGRAM:    DIVHOA      DIVERGENCE/HORIZONTAL OMEGA-ALPHA
!   PRGRMMR: JANJIC          ORG: W/NP22     DATE: 93-10-28

! ABSTRACT:
!     DIVHOA COMPUTES THE DIVERGENCE INCLUDING THE
!     MODIFICATION PREVENTING GRAVITY WAVE GRID SEPARATION, AND
!     CALCULATES THE HORIZONTAL PART OF THE OMEGA-ALPHA TERM
!     (THE PART PROPORTIONAL TO THE ADVECTION OF MASS ALONG
!     ETA/SIGMA SURFACES).

!     MODIFIED TO INCLUDE DIVERGENCE RESULTING FROM SLOPING STEPS
!     ("S": DIVHOAST)
!     EXPANDED to do also Slantwise T Advection ("T": DIVHOAST), along
!     with appropriate omega-alpha changes to the advected T.
!     Warning: this makes Subroutine SLADVT redundant
!        (needs to be commented out or removed)
!***********************************************************************
! PROGRAM HISTORY LOG:
!   87-06-??  JANJIC     - ORIGINATOR
!   95-03-25  BLACK      - CONVERSION FROM 1-D TO 2-D IN HORIZONTAL
!   96-03-29  BLACK      - ADDED EXTERNAL EDGE
!   97-03-17  MESINGER   - SPLIT FROM PFDHT
!   98-10-30  BLACK      - MODIFIED FOR DISTRIBUTED MEMORY
!   00-10-20  BLACK      - INCORPORATED PRESSURE GRADIENT METHOD
!                          FROM MESO MODEL

!   04  MESINGER and JOVIC - ADDED DIVERGENCE CALCULATION DUE TO SLOPES
!   May 2006  MESINGER, at CPTEC - EXPANDED TO INCLUDE SLANTWISE T ADV.
!   June 2006  MESINGER and JOVIC - REVISED OMEGA-ALPHA DUE TO SLOPES

! USAGE: CALL DIVHOA FROM MAIN PROGRAM EBU
!   INPUT ARGUMENT LIST:
!       NONE

!   OUTPUT ARGUMENT LIST:
!     NONE

!   OUTPUT FILES:
!     NONE

!   SUBPROGRAMS CALLED:

!     UNIQUE: NONE

!     LIBRARY: NONE

!   COMMON BLOCKS: CTLBLK
!                  MASKS
!                  LOOPS
!                  DYNAM
!                  VRBLS
!                  CONTIN
!                  INDX

! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!   MACHINE : IBM SP
!$$$
!***********************************************************************
!-----------------------------------------------------------------------
    INCLUDE "parmeta.f90"
    INCLUDE "mpp.h"
!-----------------------------------------------------------------------
    PARAMETER &
    (LP1=LM+1,JAM=6+2*(JM-10),RFCP=1.E0/(4.E0*1004.6E0))
!-----------------------------------------------------------------------
    LOGICAL :: &
    RUN,FIRST,RESTRT,SIGMA
!----------------------------------------------------------------------
    INCLUDE "COMM_CTLBLK.f90"
!-----------------------------------------------------------------------
    include "COMM_LOOPS.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_MASKS.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_INDX.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_DYNAM.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_VRBLS.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_CONTIN.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_NHYDRO.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_CLDWTR.f90"
!-----------------------------------------------------------------------
    INCLUDE "COMM_SLOPES.f90"
!-----------------------------------------------------------------------
    REAL :: &
    PINTLG(IDIM1:IDIM2,JDIM1:JDIM2,LM+1)

    REAL :: &
    FIM   (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,FILO  (IDIM1:IDIM2,JDIM1:JDIM2),RDPD  (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,ADPDX (IDIM1:IDIM2,JDIM1:JDIM2),RDPDX (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,ADPDY (IDIM1:IDIM2,JDIM1:JDIM2),RDPDY (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,ADPDNE(IDIM1:IDIM2,JDIM1:JDIM2),ADPDSE(IDIM1:IDIM2,JDIM1:JDIM2) &
    ,PEW   (IDIM1:IDIM2,JDIM1:JDIM2),PNS   (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,PCEW  (IDIM1:IDIM2,JDIM1:JDIM2),PCNS  (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,DPFEW (IDIM1:IDIM2,JDIM1:JDIM2),DPFNS (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,FNS   (IDIM1:IDIM2,JDIM1:JDIM2),TNS   (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,HM    (IDIM1:IDIM2,JDIM1:JDIM2),VM    (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,EDIV  (IDIM1:IDIM2,JDIM1:JDIM2),DIVL  (IDIM1:IDIM2,JDIM1:JDIM2)

    REAL :: &
    DPDE  (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,APEL  (IDIM1:IDIM2,JDIM1:JDIM2),PCXC  (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,ALP1  (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,DFDZ  (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,UDY   (IDIM1:IDIM2,JDIM1:JDIM2),VDX   (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,TEW   (IDIM1:IDIM2,JDIM1:JDIM2),FEW   (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,TNE   (IDIM1:IDIM2,JDIM1:JDIM2),TSE   (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,FNE   (IDIM1:IDIM2,JDIM1:JDIM2),FSE   (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,PNE   (IDIM1:IDIM2,JDIM1:JDIM2),PSE   (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,CNE   (IDIM1:IDIM2,JDIM1:JDIM2),CSE   (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,PPNE  (IDIM1:IDIM2,JDIM1:JDIM2),PPSE  (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,PCNE  (IDIM1:IDIM2,JDIM1:JDIM2),PCSE  (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,DIVS  (IDIM1:IDIM2,JDIM1:JDIM2),DPDEP1(IDIM1:IDIM2,JDIM1:JDIM2) &
    ,ATC   (IDIM1:IDIM2,JDIM1:JDIM2),ATCP1 (IDIM1:IDIM2,JDIM1:JDIM2) &
! Chou: including specific humidity and liquid water contributions
    ,AQC   (IDIM1:IDIM2,JDIM1:JDIM2),AQCP1 (IDIM1:IDIM2,JDIM1:JDIM2) &
    ,ACWMC   (IDIM1:IDIM2,JDIM1:JDIM2),ACWMCP1 (IDIM1:IDIM2,JDIM1:JDIM2)
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
    CALL ZERO2(ALP1)
    CALL ZERO2(DPDE)
    CALL ZERO2(APEL)
    CALL ZERO2(PCXC)
    CALL ZERO2(DFDZ)
    CALL ZERO2(UDY)
    CALL ZERO2(VDX)
    CALL ZERO2(TEW)
    CALL ZERO2(FEW)
    CALL ZERO2(TNE)
    CALL ZERO2(TSE)
    CALL ZERO2(FNE)
    CALL ZERO2(FSE)
    CALL ZERO2(PNE)
    CALL ZERO2(PSE)
    CALL ZERO2(CNE)
    CALL ZERO2(CSE)
    CALL ZERO2(PPNE)
    CALL ZERO2(PPSE)
    CALL ZERO2(PCNE)
    CALL ZERO2(PCSE)
!-----------------------------------------------------------------------
!--------------PREPARATORY CALCULATIONS---------------------------------
!-----------------------------------------------------------------------
    IF(SIGMA)THEN
    ! omp parallel do
        DO 50 J=MYJS_P5,MYJE_P5
            DO 50 I=MYIS_P5,MYIE_P5
                FILO(I,J)=FIS(I,J)
                PDSL(I,J)=PD(I,J)
        50 END DO
    ELSE
    ! omp parallel do
        DO 100 J=MYJS_P5,MYJE_P5
            DO 100 I=MYIS_P5,MYIE_P5
                FILO(I,J)=0.0
                PDSL(I,J)=RES(I,J)*PD(I,J)
        100 END DO
    ENDIF

    IF(HYDRO)THEN
    ! omp parallel do
        DO L=1,LM+1
            DO J=MYJS_P5,MYJE_P5
                DO I=MYIS_P5,MYIE_P5
                    PINTLG(I,J,L)=ALOG(ETA(L)*PDSL(I,J)+PT)
                ENDDO
            ENDDO
        ENDDO
    ELSE
    ! omp parallel do
        DO L=1,LM+1
            DO J=MYJS_P5,MYJE_P5
                DO I=MYIS_P5,MYIE_P5
                    PINTLG(I,J,L)=ALOG(PINT(I,J,L))
                ENDDO
            ENDDO
        ENDDO
    ENDIF

! omp parallel do
    DO L=1,LM
        DO J=MYJS_P4,MYJE_P4
            DO I=MYIS_P4,MYIE_P4
                OMGALF(I,J,L)=0.
                DIV(I,J,L)=0.
            ENDDO
        ENDDO
    ENDDO

    DO J=MYJS_P4,MYJE_P4
        DO I=MYIS_P4,MYIE_P4
            DIVS(I,J)=0.
        !   DIVergence to Save is needed to collect contributions to DIV due to
        !   slopes at h points of the layer being processed (L).  They can be
        !   added only once the code moves one layer up.  The reason is that the
        !   collection is done for points also sideways of the h point that the
        !   code is processing, so that if added immediatelly these would be
        !   lost when the code visits these sideways located points, and puts
        !   the divergence "correction" into DIV (loop 270)

        !   Accumulated T changes at L and at L+1
            ATC  (I,J)=0.0
            ATCP1(I,J)=0.0
! Chou: including specific humidity and liquid water fluxes
            AQC  (I,J)=0.0
            AQCP1(I,J)=0.0
            ACWMC  (I,J)=0.0
            ACWMCP1(I,J)=0.0
        ENDDO
    ENDDO

!-----------------------------------------------------------------------
! omp parallel do private (alp1x)
    DO J=MYJS_P5,MYJE_P5
        DO I=MYIS_P5,MYIE_P5
            ALP1X=PINTLG(I,J,LM+1)
            ALP1(I,J)=ALP1X
        ENDDO
    ENDDO
!-----------------------------------------------------------------------
!-------------------- MAIN VERTICAL INTEGRATION LOOP -------------------
!-----------------------------------------------------------------------
    DO 400 L=LM,1,-1
    !-----------------------------------------------------------------------
    !***
    !***  INTEGRATE THE GEOPOTENTIAL
    !***
    ! p
        FIM=0.
    ! p
    
    ! omp parallel do private (alp1p,dfi,fiupk,rdpds)
    !	write(6,*) 'FIM defined over I: ', MYIS_P5,MYIE_P5
    !	write(6,*) 'FIM defined over J: ', MYJS_P5,MYJE_P5
        DO 125 J=MYJS_P5,MYJE_P5
            DO 125 I=MYIS_P5,MYIE_P5
            
                ALP1P=PINTLG(I,J,L)
            
                DFI=((Q(I,J,L)*0.608+1.)*T(I,J,L)*R*(ALP1(I,J)-ALP1P)/(1+CWM(I,J,L)))/DWDT(I,J,L)
            
                RDPDS=1./(DETA(L)*PDSL(I,J))
                RTOP(I,J,L)=RDPDS*DFI
                FIUPK=FILO(I,J)+DFI
                FIM(I,J)=FILO(I,J)+FIUPK
                if (abs(FIM(I,J)) <= 5.e+10) then
                else
                    write(6,*) 'bad FIM ', I,J,FIM(I,J),FILO(I,J),DFI
                    write(6,*) 'Q,T,ALP1,ALP1P,DWDT: ', &
                    Q(I,J,L),T(I,J,L),ALP1(I,J),ALP1P,DWDT(I,J,L)
                    STOP
                endif
            
                FILO(I,J)=(FIUPK-DFL(L))*HTM(I,J,L)+DFL(L)
                ALP1(I,J)=ALP1P
        125 END DO
    
    !-----------------------------------------------------------------------
    ! omp parallel do private (alp1p,alp1pl,alp2p,alp2pl,dfi)
        DO 205 J=MYJS_P5,MYJE_P5
            DO 205 I=MYIS_P5,MYIE_P5
                HM(I,J)=HTM(I,J,L)*HBM2(I,J)
                VM(I,J)=VTM(I,J,L)*VBM2(I,J)
            
                ALP1P =PINTLG(I,J,L)
                ALP1PL=PINTLG(I,J,L+1)
                ALP2P =ALP1P*ALP1P
                ALP2PL=ALP1PL*ALP1PL
            
                DFI=((Q(I,J,L)*0.608+1.-CWM(I,J,L))*T(I,J,L)*R*(ALP1PL-ALP1P)/(1+CWM(I,J,L)))/DWDT(I,J,L)
                DFDZ(I,J)=DFI*DWDT(I,J,L)/(ALP2PL-ALP2P)
                APEL(I,J)=(ALP2PL+ALP2P)*0.5
        205 END DO
    
    ! omp parallel do
        DO 210 J=MYJS_P4,MYJE_P4
            DO 210 I=MYIS_P4,MYIE_P4
                DPDE(I,J)=DETA(L)*PDSL(I,J)
                DIVL(I,J)=0.
                EDIV(I,J)=0.
        210 END DO
    
        IF (L < LM) THEN
        ! omp parallel do
            DO 211 J=MYJS_P4,MYJE_P4
                DO 211 I=MYIS_P4,MYIE_P4
                    DPDEP1(I,J)=DETA(L+1)*PDSL(I,J)
            211 END DO
        END IF
    
    ! omp parallel do
        DO 215 J=MYJS_P1,MYJE_P1
            DO 215 I=MYIS_P1,MYIE_P1
                RDPD(I,J)=1./DPDE(I,J)
        215 END DO
    
    ! omp parallel do
        DO 220 J=MYJS1_P3,MYJE1_P3
            DO 220 I=MYIS_P3,MYIE_P3
                ADPDX(I,J)=DPDE(I+IVW(J),J)+DPDE(I+IVE(J),J)
                ADPDY(I,J)=DPDE(I,J-1)+DPDE(I,J+1)
                RDPDX(I,J)=1./ADPDX(I,J)
                RDPDY(I,J)=1./ADPDY(I,J)
        220 END DO
    
    !--------------DIAGONAL CONTRIBUTIONS TO PRESSURE GRADIENT FORCE--------
    
    ! omp parallel do
        DO 240 J=MYJS_P4,MYJE_P4
            DO 240 I=MYIS_P4,MYIE_P4
                ADPDNE(I,J)=DPDE(I+IHE(J),J+1)+DPDE(I,J)
                PNE(I,J)=(FIM (I+IHE(J),J+1)-FIM (I,J)) &
                *(DWDT(I+IHE(J),J+1,L)+DWDT(I,J,L))
                if ( ABS(PNE(I,J)) <= 5.e10) then
                else
                    write(6,*) 'crazy PNE ',I,J,PNE(I,J)
                    write(6,*) 'pieces', I+IHE(J),J+1,FIM (I+IHE(J),J+1)
                endif
                PPNE(I,J)=PNE(I,J)*ADPDNE(I,J)
                CNE(I,J)=(DFDZ(I+IHE(J),J+1)+DFDZ(I,J))*2. &
                *(APEL(I+IHE(J),J+1)-APEL(I,J))
                PCNE(I,J)=CNE(I,J)*ADPDNE(I,J)
        240 END DO
    
    ! omp parallel do
        DO 250 J=MYJS1_P4,MYJE_P4
            DO 250 I=MYIS_P4,MYIE1_P4
                ADPDSE(I,J)=DPDE(I+IHE(J),J-1)+DPDE(I,J)
                PSE(I,J)=(FIM(I+IHE(J),J-1)-FIM(I,J)) &
                *(DWDT(I+IHE(J),J-1,L)+DWDT(I,J,L))
                PPSE(I,J)=PSE(I,J)*ADPDSE(I,J)
                CSE(I,J)=(DFDZ(I+IHE(J),J-1)+DFDZ(I,J))*2. &
                *(APEL(I+IHE(J),J-1)-APEL(I,J))
                PCSE(I,J)=CSE(I,J)*ADPDSE(I,J)
        250 END DO
    
    !--------------CONTINUITY EQUATION MODIFICATION-------------------------
    
    ! omp parallel do
        DO 260 J=MYJS1_P1,MYJE1_P1
            DO 260 I=MYIS_P1,MYIE_P1
                PCXC(I,J)=VBM3(I,J)*VTM(I,J,L)*(PNE(I+IVW(J),J) &
                +CNE(I+IVW(J),J)+PSE(I+IVW(J),J)+CSE(I+IVW(J),J) &
                -PNE(I,J-1)-CNE(I,J-1)-PSE(I,J+1)-CSE(I,J+1))
        260 END DO
    !-----------------------------------------------------------------------
        DO 270 J=MYJS2,MYJE2
            DO 270 I=MYIS1,MYIE1
                DIV(I,J,L)=DETA(L)*WPDAR(I,J) &
                *(PCXC(I+IHE(J),J)-PCXC(I,J+1) &
                +PCXC(I+IHW(J),J)-PCXC(I,J-1))
        270 END DO
    
    !--------------LAT & LONG PRESSURE FORCE COMPONENTS---------------------
    
    ! omp parallel do private (dcnek,dcsek,dpnek,dpsek)
        DO 280 J=MYJS1_P3,MYJE1_P3
            DO 280 I=MYIS_P3,MYIE_P3
                DPNEK=PNE(I+IVW(J),J)+PNE(I,J-1)
                DPSEK=PSE(I+IVW(J),J)+PSE(I,J+1)
                PEW(I,J)=DPNEK+DPSEK
                PNS(I,J)=DPNEK-DPSEK
                DCNEK=CNE(I+IVW(J),J)+CNE(I,J-1)
                DCSEK=CSE(I+IVW(J),J)+CSE(I,J+1)
                PCEW(I,J)=(DCNEK+DCSEK)*ADPDX(I,J)
                PCNS(I,J)=(DCNEK-DCSEK)*ADPDY(I,J)
        280 END DO
    
    !--------------LAT & LON FLUXES & OMEGA-ALPHA COMPONENTS----------------
    
    ! omp parallel do
        DO 310 J=MYJS1_P3,MYJE1_P3
            DO 310 I=MYIS_P3,MYIE_P3
                UDY(I,J)=DY*U(I,J,L)
                FEW(I,J)=UDY(I,J)*ADPDX(I,J)
                TEW(I,J)=UDY(I,J)*PCEW(I,J)
                VDX(I,J)=DX(I,J)*V(I,J,L)
                FNS(I,J)=VDX(I,J)*ADPDY(I,J)
                TNS(I,J)=VDX(I,J)*PCNS(I,J)
        310 END DO
    
    !--------------DIAGONAL FLUXES AND DIAGONALLY AVERAGED WIND-------------
    
    ! omp parallel do private (pvnek)
        DO 320 J=MYJS1_P2,MYJE2_P2
            DO 320 I=MYIS_P2,MYIE1_P2
                PVNEK=(UDY(I+IHE(J),J)+VDX(I+IHE(J),J))+(UDY(I,J+1)+VDX(I,J+1))
                FNE(I,J)=PVNEK*ADPDNE(I,J)
                TNE(I,J)=PVNEK*PCNE(I,J)*2.
        320 END DO
    
    ! omp parallel do private (pvsek)
        DO 330 J=MYJS2_P2,MYJE1_P2
            DO 330 I=MYIS_P2,MYIE1_P2
                PVSEK=(UDY(I+IHE(J),J)-VDX(I+IHE(J),J))+(UDY(I,J-1)-VDX(I,J-1))
                FSE(I,J)=PVSEK*ADPDSE(I,J)
                TSE(I,J)=PVSEK*PCSE(I,J)*2.
        330 END DO
        IF (L < LM-1) THEN
            DO 335 J=MYJS2_P2,MYJE2_P2
                DO 335 I=MYIS1_P2,MYIE1_P2
                    DIV(I,J,L+1)=DIV(I,J,L+1)+DIVS(I,J)
                    DIVS(I,J)=0.0
                    T(I,J,L+2)=T(I,J,L+2)+ATCP1(I,J)
                    Q(I,J,L+2)=Q(I,J,L+2)+AQCP1(I,J)         !chou
                    CWM(I,J,L+2)=CWM(I,J,L+2)+ACWMCP1(I,J)   !chou

                    ATCP1(I,J)=ATC(I,J)
                    ATC  (I,J)=0.0
! Chou: including specific humidity and liquid water contributions
                    AQCP1(I,J)=AQC(I,J)        !chou
                    AQC  (I,J)=0.0             !chou
                    ACWMCP1(I,J)=ACWMC(I,J)    !chou
                    ACWMC  (I,J)=0.0           !chou

            335 END DO
        END IF
    
    !--------------HORIZONTAL PART OF OMEGA-ALPHA & DIVERGENCE--------------
    
    ! omp parallel do
        DO 340 J=MYJS2_P1,MYJE2_P1
            DO 340 I=MYIS1_P1,MYIE1_P1
                OMGALF(I,J,L)=(TEW(I+IHE(J),J)+TEW(I+IHW(J),J)+TNS(I,J+1) &
                +TNS(I,J-1)+TNE(I,J)+TNE(I+IHW(J),J-1)+TSE(I,J) &
                +TSE(I+IHW(J),J+1))*RDPD(I,J)*FCP(I,J)*HM(I,J)
                IF (L < 3) THEN
                    T(I,J,L)=OMGALF(I,J,L)+T(I,J,L)
                ELSEIF (L < LM) THEN
                    ATC(I,J)=ATC(I,J)+OMGALF(I,J,L)
                END IF
            
                IF (L < LM .AND. (HTM(I,J,L+1)*HBM2(I,J)) > 0) THEN
                !   The T point below is in the atmosphere.  Should it receive
                !   omega-alpha contributions from fluxes of layer L?

                !   Check each of the four surrounding V points if they are
                !   the points just above the slope

                    TOASL=0.0

                    IF (VTMS(I+IHE(J),J,L+1) == 1) THEN
                        IF (ISLD(I+IHE(J),J) == 1) &
                        TOASL=TOASL+0.50*TEW(I+IHE(J),J)
                        IF (ISLD(I+IHE(J),J) == 8) &
                        TOASL=TOASL+0.50*TEW(I+IHE(J),J)+0.25*TSE(I,J)
                        IF (ISLD(I+IHE(J),J) == 2) &
                        TOASL=TOASL+0.50*TEW(I+IHE(J),J)+0.25*TNE(I,J)
                    ENDIF

                    IF (VTMS(I,J+1,L+1) == 1) THEN
                        IF (ISLD(I,J+1) == 3) &
                        TOASL=TOASL+0.50*TNS(I,J+1)
                        IF (ISLD(I,J+1) == 2) &
                        TOASL=TOASL+0.50*TNS(I,J+1)+0.25*TNE(I,J)
                        IF (ISLD(I,J+1) == 4) &
                        TOASL=TOASL+0.50*TNS(I,J+1)+0.25*TSE(I+IHW(J),J+1)
                    ENDIF

                    IF (VTMS(I+IHW(J),J,L+1) == 1) THEN
                        IF (ISLD(I+IHW(J),J) == 5) &
                        TOASL=TOASL+0.50*TEW(I+IHW(J),J)
                        IF (ISLD(I+IHW(J),J) == 4) &
                        TOASL=TOASL+0.50*TEW(I+IHW(J),J)+0.25*TSE(I+IHW(J),J+1)
                        IF (ISLD(I+IHW(J),J) == 6) &
                        TOASL=TOASL+0.50*TEW(I+IHW(J),J)+0.25*TNE(I+IHW(J),J-1)
                    ENDIF

                    IF (VTMS(I,J-1,L+1) == 1) THEN
                        IF (ISLD(I,J-1) == 7) &
                        TOASL=TOASL+0.50*TNS(I,J-1)
                        IF (ISLD(I,J-1) == 6) &
                        TOASL=TOASL+0.50*TNS(I,J-1)+0.25*TNE(I+IHW(J),J-1)
                        IF (ISLD(I,J-1) == 8) &
                        TOASL=TOASL+0.50*TNS(I,J-1)+0.25*TSE(I,J)
                    ENDIF

                    OAADP1 = TOASL*DETA(L+1)/DETA(L)*RDPD(I,J)*FCP(I,J)
                    OMGALF(I,J,L+1)= OMGALF(I,J,L+1) +OAADP1
                    IF (L < 3) THEN
                        T(I,J,L+1)=T(I,J,L+1) +OAADP1
                    ELSEIF (L == LM-1) THEN
                        ATCP1(I,J)=ATCP1(I,J) + OMGALF(I,J,L+1)
                    ELSE
                        ATCP1(I,J)=ATCP1(I,J) +OAADP1
                    END IF

                ENDIF
            
                EDIV(I,J)=((FEW(I+IHE(J),J)+FNS(I,J+1) &
                +FNE(I,J)+FSE(I,J)) &
                -(FEW(I+IHW(J),J)+FNS(I,J-1) &
                +FNE(I+IHW(J),J-1)+FSE(I+IHW(J),J+1)))*FDIV(I,J)
                DIVL(I,J)=EDIV(I,J)*HBM2(I,J)
                 
            !  Slantwise mass and T advection, with T increment to account for the
            !  omega-alpha impact, T change due to the movement to a different p

                IF (L < LM .AND. L == LMV(I+IHE(J),J) &
                 .AND. ISLD(I+IHE(J),J) > 0) THEN

                    TFSTR=0.5*DETA(L+1)/DETA(L)  * FDIV(I,J)*HM(I,J)
                    FSTR=0.5*TFSTR

                    STSEP1=FSTR*FSE(I+IHE(J),J+1)
                    STNE  =FSTR*FNE(I,J)
                    STSE  =FSTR*FSE(I,J)
                    STNEM1=FSTR*FNE(I+IHE(J),J-1)
                    STEW  =TFSTR*FEW(I+IHE(J),J)
                    STNS  =TFSTR*FNS(I+IHE(J),J)

                    IF      (ISLD(I+IHE(J),J) == 1) THEN
                        DIV(I+IHE(J),J+1,L+1) = DIV(I+IHE(J),J+1,L+1) + STSEP1
                        DIVS(I+1,J)           = DIVS(I+1,J)           - STSEP1
                        wdep=ABS( STSEP1)*DT
                        DPT=(RTOP(I+IHE(J),J+1,L+1)+RTOP(I+1,J,L))* &
                        (DPDEP1(I+IHE(J),J+1)    +DPDE(I+1,J))  *RFCP
                        WDUTD=wdep*SIGN(1.0, STSEP1) &
                        *(T(I+IHE(J),J+1,L+1)-DPT-T(I+1,J,L))

! Chou: including specific humidity and liquid water contributions
                        WDUQD=wdep*SIGN(1.0, STSEP1) *(Q(I+IHE(J),J+1,L+1)-Q(I+1,J,L))
                        WDUCWMD=wdep*SIGN(1.0, STSEP1) *(CWM(I+IHE(J),J+1,L+1)-CWM(I+1,J,L))

                        IF ( STSEP1 > 0.0) THEN
                            warr=DPDE(I+1,J)
                            ATC(I+1,J)=ATC(I+1,J)+WDUTD/(warr+wdep)
                            AQC(I+1,J)=AQC(I+1,J)+WDUQD/(warr+wdep)
                            ACWMC(I+1,J)=ACWMC(I+1,J)+WDUCWMD/(warr+wdep)
                        ELSE
                            warr=DPDEP1(I+IHE(J),J+1)
                            ATCP1(I+IHE(J),J+1)=ATCP1(I+IHE(J),J+1)+WDUTD/(warr+wdep)
                            AQCP1(I+IHE(J),J+1)=AQCP1(I+IHE(J),J+1)+WDUQD/(warr+wdep)
                            ACWMCP1(I+IHE(J),J+1)=ACWMCP1(I+IHE(J),J+1)+WDUCWMD/(warr+wdep)
                        END IF

                        DIV(I+IHE(J),J-1,L+1) = DIV(I+IHE(J),J-1,L+1) + STNEM1
                        DIVS(I+1,J)           = DIVS(I+1,J)           - STNEM1
                        wdep=ABS( STNEM1)*DT
                        DPT=(RTOP(I+IHE(J),J-1,L+1)+RTOP(I+1,J,L))* &
                        (DPDEP1(I+IHE(J),J-1)    +DPDE(I+1,J))  *RFCP
                        WDUTD=wdep*SIGN(1.0, STNEM1) &
                        *(T(I+IHE(J),J-1,L+1)-DPT-T(I+1,J,L))

! Chou: including specific humidity and liquid water contributions
                        WDUQD=wdep*SIGN(1.0,  STNEM1)*(Q(I+IHE(J),J-1,L+1)-Q(I+1,J,L))
                        WDUCWMD=wdep*SIGN(1.0,STNEM1)*(CWM(I+IHE(J),J-1,L+1)-CWM(I+1,J,L))

                        IF ( STNEM1 > 0.0) THEN
                            warr=DPDE(I+1,J)
                            ATC(I+1,J)=ATC(I+1,J)+WDUTD/(warr+wdep)
                            AQC(I+1,J)=AQC(I+1,J)+WDUQD/(warr+wdep)
                            ACWMC(I+1,J)=ACWMC(I+1,J)+WDUCWMD/(warr+wdep)
                        ELSE
                            warr=DPDEP1(I+IHE(J),J-1)
                            ATCP1(I+IHE(J),J-1)=ATCP1(I+IHE(J),J-1)+WDUTD/(warr+wdep)
                            AQCP1(I+IHE(J),J-1)=AQCP1(I+IHE(J),J-1)+WDUQD/(warr+wdep)
                            ACWMCP1(I+IHE(J),J-1)=ACWMCP1(I+IHE(J),J-1)+WDUCWMD/(warr+wdep)
                        END IF

                        DIV(I,J,L+1)          = DIV(I,J,L+1)          + STEW  !blue
                        DIVS(I+1,J)           = DIVS(I+1,J)           - STEW  !blue
                        wdep=ABS(   STEW)*DT
                        DPT=(RTOP(I,J,L+1)+RTOP(I+1,J,L))* &
                        (DPDEP1(I,J)    +DPDE(I+1,J))  *RFCP
                        WDUTD=wdep*SIGN(1.0,   STEW)*(T(I,J,L+1)-DPT-T(I+1,J,L))

! Chou: including specific humidity and liquid water contributions
                        WDUQD  =wdep*SIGN(1.0,   STEW)*(Q(I,J,L+1)-Q(I+1,J,L))
                        WDUCWMD=wdep*SIGN(1.0,   STEW)*(CWM(I,J,L+1)-CWM(I+1,J,L))
                        IF (   STEW > 0.0) THEN
                            warr=DPDE(I+1,J)
                            ATC(I+1,J)=ATC(I+1,J)+WDUTD/(warr+wdep)
                            AQC(I+1,J)=AQC(I+1,J)+WDUQD/(warr+wdep)
                            ACWMC(I+1,J)=ACWMC(I+1,J)+WDUCWMD/(warr+wdep)
                        ELSE
                            warr=DPDEP1(I,J)
                            ATCP1(I,J)=ATCP1(I,J)+WDUTD/(warr+wdep)
                            AQCP1(I,J)=AQCP1(I,J)+WDUQD/(warr+wdep)
                            ACWMCP1(I,J)=ACWMCP1(I,J)+WDUCWMD/(warr+wdep)
                        END IF

                    ELSE IF (ISLD(I+IHE(J),J) == 2) THEN
                        DIV(I,J,L+1)          = DIV(I,J,L+1)          + STNE
                        DIVS(I+IHE(J),J+1)    = DIVS(I+IHE(J),J+1)    - STNE
                        wdep=ABS(   STNE)*DT
                        DPT=(RTOP(I,J,L+1)+RTOP(I+IHE(J),J+1,L))* &
                        (DPDEP1(I,J)    +DPDE(I+IHE(J),J+1))  *RFCP
                        WDUTD=wdep*SIGN(1.0,   STNE)*(T(I,J,L+1)-DPT-T(I+IHE(J),J+1,L))

! Chou: including specific humidity and liquid water contributions
                        WDUQD=wdep*SIGN(1.0,   STNE)*(Q(I,J,L+1)-Q(I+IHE(J),J+1,L))
                        WDUCWMD=wdep*SIGN(1.0, STNE)*(CWM(I,J,L+1)-CWM(I+IHE(J),J+1,L))

                        IF (   STNE > 0.0) THEN
                            warr=DPDE(I+IHE(J),J+1)
                            ATC(I+IHE(J),J+1)=ATC(I+IHE(J),J+1)+WDUTD/(warr+wdep)
                            AQC(I+IHE(J),J+1)=AQC(I+IHE(J),J+1)+WDUQD/(warr+wdep)
                            ACWMC(I+IHE(J),J+1)=ACWMC(I+IHE(J),J+1)+WDUCWMD/(warr+wdep)
                        ELSE
                            warr=DPDEP1(I,J)
                            ATCP1(I,J)=ATCP1(I,J)+WDUTD/(warr+wdep)
                            AQCP1(I,J)=AQCP1(I,J)+WDUQD/(warr+wdep)
                            ACWMCP1(I,J)=ACWMCP1(I,J)+WDUCWMD/(warr+wdep)
                        END IF
                        		   
                        DIV(I+IHE(J),J-1,L+1) = DIV(I+IHE(J),J-1,L+1) + STNEM1
                        DIVS(I+1,J)           = DIVS(I+1,J)           - STNEM1
                        wdep=ABS( STNEM1)*DT
                        DPT=(RTOP(I+IHE(J),J-1,L+1)+RTOP(I+1,J,L))* &
                        (DPDEP1(I+IHE(J),J-1)    +DPDE(I+1,J))  *RFCP
                        WDUTD=wdep*SIGN(1.0, STNEM1) &
                        *(T(I+IHE(J),J-1,L+1)-DPT-T(I+1,J,L))

! Chou: including specific humidity and liquid water contributions
                        WDUQD=wdep*SIGN(1.0,  STNEM1)*(Q(I+IHE(J),J-1,L+1)-Q(I+1,J,L))
                        WDUCWMD=wdep*SIGN(1.0,STNEM1)*(CWM(I+IHE(J),J-1,L+1)-CWM(I+1,J,L))

                        IF ( STNEM1 > 0.0) THEN
                            warr=DPDE(I+1,J)
                            ATC(I+1,J)=ATC(I+1,J)+WDUTD/(warr+wdep)
                            AQC(I+1,J)=AQC(I+1,J)+WDUQD/(warr+wdep)
                            ACWMC(I+1,J)=ACWMC(I+1,J)+WDUCWMD/(warr+wdep)
                        ELSE
                            warr=DPDEP1(I+IHE(J),J-1)
                            ATCP1(I+IHE(J),J-1)=ATCP1(I+IHE(J),J-1)+WDUTD/(warr+wdep)
                            AQCP1(I+IHE(J),J-1)=AQCP1(I+IHE(J),J-1)+WDUQD/(warr+wdep)
                            ACWMCP1(I+IHE(J),J-1)=ACWMCP1(I+IHE(J),J-1)+WDUCWMD/(warr+wdep)
                        END IF

                        DIV(I,J,L+1)          = DIV(I,J,L+1)          + STEW  !blue
                        DIVS(I+1,J)           = DIVS(I+1,J)           - STEW  !blue
                        wdep=ABS(   STEW)*DT
                        DPT=(RTOP(I,J,L+1)+RTOP(I+1,J,L))* &
                        (DPDEP1(I,J)    +DPDE(I+1,J))  *RFCP
                        WDUTD=wdep*SIGN(1.0,   STEW)*(T(I,J,L+1)-DPT-T(I+1,J,L))

! Chou: including specific humidity and liquid water contributions
                        WDUQD=wdep*SIGN(1.0,   STEW)*(Q(I,J,L+1)-Q(I+1,J,L))
                        WDUCWMD=wdep*SIGN(1.0, STEW)*(CWM(I,J,L+1)-CWM(I+1,J,L))

                        IF (   STEW > 0.0) THEN
                            warr=DPDE(I+1,J)
                            ATC(I+1,J)=ATC(I+1,J)+WDUTD/(warr+wdep)
                            AQC(I+1,J)=AQC(I+1,J)+WDUQD/(warr+wdep)
                            ACWMC(I+1,J)=ACWMC(I+1,J)+WDUCWMD/(warr+wdep)
                        ELSE
                            warr=DPDEP1(I,J)
                            ATCP1(I,J)=ATCP1(I,J)+WDUTD/(warr+wdep)
                            AQCP1(I,J)=AQCP1(I,J)+WDUQD/(warr+wdep)
                            ACWMCP1(I,J)=ACWMCP1(I,J)+WDUCWMD/(warr+wdep)
                        END IF

                        DIV(I+IHE(J),J-1,L+1) = DIV(I+IHE(J),J-1,L+1) + STNS  !yellow
                        DIVS(I+IHE(J),J+1)    = DIVS(I+IHE(J),J+1)    - STNS  !yellow
                        wdep=ABS(  STNS)*DT
                        DPT=(RTOP(I+IHE(J),J-1,L+1)+RTOP(I+IHE(J),J+1,L))* &
                        (DPDEP1(I+IHE(J),J-1)    +DPDE(I+IHE(J),J+1))  *RFCP
                        WDUTD=wdep*SIGN(1.0,  STNS) &
                        *(T(I+IHE(J),J-1,L+1)-DPT-T(I+IHE(J),J+1,L))

! Chou: including specific humidity and liquid water contributions
                        WDUQD=wdep*SIGN(1.0,  STNS)*(Q(I+IHE(J),J-1,L+1)-Q(I+IHE(J),J+1,L))
                        WDUCWMD=wdep*SIGN(1.0,STNS)*(CWM(I+IHE(J),J-1,L+1)-CWM(I+IHE(J),J+1,L))

                        IF (  STNS > 0.0) THEN
                            warr=DPDE(I+IHE(J),J+1)
                            ATC(I+IHE(J),J+1)=ATC(I+IHE(J),J+1)+WDUTD/(warr+wdep)
                            AQC(I+IHE(J),J+1)=AQC(I+IHE(J),J+1)+WDUQD/(warr+wdep)
                            ACWMC(I+IHE(J),J+1)=ACWMC(I+IHE(J),J+1)+WDUCWMD/(warr+wdep)
                        ELSE
                            warr=DPDEP1(I+IHE(J),J-1)
                            ATCP1(I+IHE(J),J-1)=ATCP1(I+IHE(J),J-1)+WDUTD/(warr+wdep)
                            AQCP1(I+IHE(J),J-1)=AQCP1(I+IHE(J),J-1)+WDUQD/(warr+wdep)
                            ACWMCP1(I+IHE(J),J-1)=ACWMCP1(I+IHE(J),J-1)+WDUCWMD/(warr+wdep)
                        END IF

                    ELSE IF (ISLD(I+IHE(J),J) == 3) THEN
                        DIV(I+1,J,L+1)        = DIV(I+1,J,L+1)        - STSEP1
                        DIVS(I+IHE(J),J+1)    = DIVS(I+IHE(J),J+1)    + STSEP1
                        wdep=ABS(-STSEP1)*DT
                        DPT=(RTOP(I+1,J,L+1)+RTOP(I+IHE(J),J+1,L))* &
                        (DPDEP1(I+1,J)    +DPDE(I+IHE(J),J+1))  *RFCP
                        WDUTD=wdep*SIGN(1.0,-STSEP1) &
                        *(T(I+1,J,L+1)-DPT-T(I+IHE(J),J+1,L))

! Chou: including specific humidity and liquid water contributions
                        WDUQD=wdep*SIGN(1.0,  -STSEP1)*(Q(I+1,J,L+1)-Q(I+IHE(J),J+1,L))
                        WDUCWMD=wdep*SIGN(1.0,-STSEP1)*(CWM(I+1,J,L+1)-CWM(I+IHE(J),J+1,L))

                        IF (-STSEP1 > 0.0) THEN
                            warr=DPDE(I+IHE(J),J+1)
                            ATC(I+IHE(J),J+1)=ATC(I+IHE(J),J+1)+WDUTD/(warr+wdep)
                            AQC(I+IHE(J),J+1)=AQC(I+IHE(J),J+1)+WDUQD/(warr+wdep)
                            ACWMC(I+IHE(J),J+1)=ACWMC(I+IHE(J),J+1)+WDUCWMD/(warr+wdep)
                        ELSE
                            warr=DPDEP1(I+1,J)
                            ATCP1(I+1,J)=ATCP1(I+1,J)+WDUTD/(warr+wdep)
                            AQCP1(I+1,J)=AQCP1(I+1,J)+WDUQD/(warr+wdep)
                            ACWMCP1(I+1,J)=ACWMCP1(I+1,J)+WDUCWMD/(warr+wdep)
                        END IF
                        		   
                        DIV(I,J,L+1)          = DIV(I,J,L+1)          + STNE
                        DIVS(I+IHE(J),J+1)    = DIVS(I+IHE(J),J+1)    - STNE
                        wdep=ABS(   STNE)*DT
                        DPT=(RTOP(I,J,L+1)+RTOP(I+IHE(J),J+1,L))* &
                        (DPDEP1(I,J)    +DPDE(I+IHE(J),J+1))  *RFCP
                        WDUTD=wdep*SIGN(1.0,   STNE)*(T(I,J,L+1)-DPT-T(I+IHE(J),J+1,L))

! Chou: including specific humidity and liquid water contributions
                        WDUQD=wdep*SIGN(1.0,   STNE)*(Q(I,J,L+1)-Q(I+IHE(J),J+1,L))
                        WDUCWMD=wdep*SIGN(1.0, STNE)*(CWM(I,J,L+1)-CWM(I+IHE(J),J+1,L))

                        IF (   STNE > 0.0) THEN
                            warr=DPDE(I+IHE(J),J+1)
                            ATC(I+IHE(J),J+1)=ATC(I+IHE(J),J+1)+WDUTD/(warr+wdep)
                            AQC(I+IHE(J),J+1)=AQC(I+IHE(J),J+1)+WDUQD/(warr+wdep)
                            ACWMC(I+IHE(J),J+1)=ACWMC(I+IHE(J),J+1)+WDUCWMD/(warr+wdep)
                        ELSE
                            warr=DPDEP1(I,J)
                            ATCP1(I,J)=ATCP1(I,J)+WDUTD/(warr+wdep)
                            AQCP1(I,J)=AQCP1(I,J)+WDUQD/(warr+wdep)
                            ACWMCP1(I,J)=ACWMCP1(I,J)+WDUCWMD/(warr+wdep)
                        END IF

                        DIV(I+IHE(J),J-1,L+1) = DIV(I+IHE(J),J-1,L+1) + STNS  !yellow
                        DIVS(I+IHE(J),J+1)    = DIVS(I+IHE(J),J+1)    - STNS  !yellow
                        wdep=ABS(  STNS)*DT
                        DPT=(RTOP(I+IHE(J),J-1,L+1)+RTOP(I+IHE(J),J+1,L))* &
                        (DPDEP1(I+IHE(J),J-1)    +DPDE(I+IHE(J),J+1))  *RFCP
                        WDUTD=wdep*SIGN(1.0,  STNS) &
                        *(T(I+IHE(J),J-1,L+1)-DPT-T(I+IHE(J),J+1,L))

! Chou: including specific humidity and liquid water contributions
                        WDUQD=wdep*SIGN(1.0,  STNS)*(Q(I+IHE(J),J-1,L+1)-Q(I+IHE(J),J+1,L))
                        WDUCWMD=wdep*SIGN(1.0,STNS)*(CWM(I+IHE(J),J-1,L+1)-CWM(I+IHE(J),J+1,L))

                        IF (  STNS > 0.0) THEN
                            warr=DPDE(I+IHE(J),J+1)
                            ATC(I+IHE(J),J+1)=ATC(I+IHE(J),J+1)+WDUTD/(warr+wdep)
                            AQC(I+IHE(J),J+1)=AQC(I+IHE(J),J+1)+WDUQD/(warr+wdep)
                            ACWMC(I+IHE(J),J+1)=ACWMC(I+IHE(J),J+1)+WDUCWMD/(warr+wdep)
                        ELSE
                            warr=DPDEP1(I+IHE(J),J-1)
                            ATCP1(I+IHE(J),J-1)=ATCP1(I+IHE(J),J-1)+WDUTD/(warr+wdep)
                            AQCP1(I+IHE(J),J-1)=AQCP1(I+IHE(J),J-1)+WDUQD/(warr+wdep)
                            ACWMCP1(I+IHE(J),J-1)=ACWMCP1(I+IHE(J),J-1)+WDUCWMD/(warr+wdep)
                        END IF

                    ELSE IF (ISLD(I+IHE(J),J) == 4) THEN
                        DIV(I+1,J,L+1)        = DIV(I+1,J,L+1)        - STSEP1
                        DIVS(I+IHE(J),J+1)    = DIVS(I+IHE(J),J+1)    + STSEP1
                        wdep=ABS(-STSEP1)*DT
                        DPT=(RTOP(I+1,J,L+1)+RTOP(I+IHE(J),J+1,L))* &
                        (DPDEP1(I+1,J)    +DPDE(I+IHE(J),J+1))  *RFCP
                        WDUTD=wdep*SIGN(1.0,-STSEP1) &
                        *(T(I+1,J,L+1)-DPT-T(I+IHE(J),J+1,L))

! Chou: including specific humidity and liquid water contributions
                        WDUQD=wdep*SIGN(1.0,  -STSEP1)*(Q(I+1,J,L+1)-Q(I+IHE(J),J+1,L))
                        WDUCWMD=wdep*SIGN(1.0,-STSEP1)*(CWM(I+1,J,L+1)-CWM(I+IHE(J),J+1,L))

                        IF (-STSEP1 > 0.0) THEN
                            warr=DPDE(I+IHE(J),J+1)
                            ATC(I+IHE(J),J+1)=ATC(I+IHE(J),J+1)+WDUTD/(warr+wdep)
                            AQC(I+IHE(J),J+1)=AQC(I+IHE(J),J+1)+WDUQD/(warr+wdep)
                            ACWMC(I+IHE(J),J+1)=ACWMC(I+IHE(J),J+1)+WDUCWMD/(warr+wdep)
                        ELSE
                            warr=DPDEP1(I+1,J)
                            ATCP1(I+1,J)=ATCP1(I+1,J)+WDUTD/(warr+wdep)
                            AQCP1(I+1,J)=AQCP1(I+1,J)+WDUQD/(warr+wdep)
                            ACWMCP1(I+1,J)=ACWMCP1(I+1,J)+WDUCWMD/(warr+wdep)
                        END IF

                        DIV(I+IHE(J),J-1,L+1) = DIV(I+IHE(J),J-1,L+1) - STSE
                        DIVS(I,J)             = DIVS(I,J)             + STSE
                        wdep=ABS(   -STSE)*DT
                        DPT=(RTOP(I+IHE(J),J-1,L+1)+RTOP(I,J,L))* &
                        (DPDEP1(I+IHE(J),J-1)    +DPDE(I,J))  *RFCP
                        WDUTD=wdep*SIGN(1.0,   -STSE)*(T(I+IHE(J),J-1,L+1)-DPT-T(I,J,L))

! Chou: including specific humidity and liquid water contributions
                        WDUQD=wdep*SIGN(1.0,   -STSE)*(Q(I+IHE(J),J-1,L+1)-Q(I,J,L))
                        WDUCWMD=wdep*SIGN(1.0, -STSE)*(CWM(I+IHE(J),J-1,L+1)-CWM(I,J,L))

                        IF (   -STSE > 0.0) THEN
                            warr=DPDE(I,J)
                            ATC(I,J)=ATC(I,J)+WDUTD/(warr+wdep)
                            AQC(I,J)=AQC(I,J)+WDUQD/(warr+wdep)
                            ACWMC(I,J)=ACWMC(I,J)+WDUCWMD/(warr+wdep)
                        ELSE
                            warr=DPDEP1(I+IHE(J),J-1)
                            ATCP1(I+IHE(J),J-1)=ATCP1(I+IHE(J),J-1)+WDUTD/(warr+wdep)
                            AQCP1(I+IHE(J),J-1)=AQCP1(I+IHE(J),J-1)+WDUQD/(warr+wdep)
                            ACWMCP1(I+IHE(J),J-1)=ACWMCP1(I+IHE(J),J-1)+WDUCWMD/(warr+wdep)
                        END IF

                        DIV(I+IHE(J),J-1,L+1) = DIV(I+IHE(J),J-1,L+1) + STNS  !yellow
                        DIVS(I+IHE(J),J+1)    = DIVS(I+IHE(J),J+1)    - STNS  !yellow
                        wdep=ABS(  STNS)*DT
                        DPT=(RTOP(I+IHE(J),J-1,L+1)+RTOP(I+IHE(J),J+1,L))* &
                        (DPDEP1(I+IHE(J),J-1)    +DPDE(I+IHE(J),J+1))  *RFCP
                        WDUTD=wdep*SIGN(1.0,  STNS) &
                        *(T(I+IHE(J),J-1,L+1)-DPT-T(I+IHE(J),J+1,L))

! Chou: including specific humidity and liquid water contributions
                        WDUQD=wdep*SIGN(1.0,  STNS)*(Q(I+IHE(J),J-1,L+1)-Q(I+IHE(J),J+1,L))
                        WDUCWMD=wdep*SIGN(1.0,STNS)*(CWM(I+IHE(J),J-1,L+1)-CWM(I+IHE(J),J+1,L))

                        IF (  STNS > 0.0) THEN
                            warr=DPDE(I+IHE(J),J+1)
                            ATC(I+IHE(J),J+1)=ATC(I+IHE(J),J+1)+WDUTD/(warr+wdep)
                            AQC(I+IHE(J),J+1)=AQC(I+IHE(J),J+1)+WDUQD/(warr+wdep)
                            ACWMC(I+IHE(J),J+1)=ACWMC(I+IHE(J),J+1)+WDUCWMD/(warr+wdep)
                        ELSE
                            warr=DPDEP1(I+IHE(J),J-1)
                            ATCP1(I+IHE(J),J-1)=ATCP1(I+IHE(J),J-1)+WDUTD/(warr+wdep)
                            AQCP1(I+IHE(J),J-1)=AQCP1(I+IHE(J),J-1)+WDUQD/(warr+wdep)
                            ACWMCP1(I+IHE(J),J-1)=ACWMCP1(I+IHE(J),J-1)+WDUCWMD/(warr+wdep)
                        END IF

                        DIV(I+1,J,L+1)        = DIV(I+1,J,L+1)        - STEW  !green
                        DIVS(I,J)             = DIVS(I,J)             + STEW  !green
                        wdep=ABS(  -STEW)*DT
                        DPT=(RTOP(I+1,J,L+1)+RTOP(I,J,L))* &
                        (DPDEP1(I+1,J)    +DPDE(I,J))  *RFCP
                        WDUTD=wdep*SIGN(1.0,  -STEW)*(T(I+1,J,L+1)-DPT-T(I,J,L))

! Chou: including specific humidity and liquid water contributions
                        WDUQD=wdep*SIGN(1.0,  -STEW)*(Q(I+1,J,L+1)-Q(I,J,L))
                        WDUCWMD=wdep*SIGN(1.0,-STEW)*(CWM(I+1,J,L+1)-CWM(I,J,L))

                        IF (  -STEW > 0.0) THEN
                            warr=DPDE(I,J)
                            ATC(I,J)=ATC(I,J)+WDUTD/(warr+wdep)
                            AQC(I,J)=AQC(I,J)+WDUQD/(warr+wdep)
                            ACWMC(I,J)=ACWMC(I,J)+WDUCWMD/(warr+wdep)
                        ELSE
                            warr=DPDEP1(I+1,J)
                            ATCP1(I+1,J)=ATCP1(I+1,J)+WDUTD/(warr+wdep)
                            AQCP1(I+1,J)=AQCP1(I+1,J)+WDUQD/(warr+wdep)
                            ACWMCP1(I+1,J)=ACWMCP1(I+1,J)+WDUCWMD/(warr+wdep)
                        END IF

                    ELSE IF (ISLD(I+IHE(J),J) == 5) THEN
                        DIV(I+IHE(J),J+1,L+1) = DIV(I+IHE(J),J+1,L+1) - STNE
                        DIVS(I,J)             = DIVS(I,J)             + STNE
                        wdep=ABS(  -STNE)*DT
                        DPT=(RTOP(I+IHE(J),J+1,L+1)+RTOP(I,J,L))* &
                        (DPDEP1(I+IHE(J),J+1)    +DPDE(I,J))  *RFCP
                        WDUTD=wdep*SIGN(1.0,  -STNE)*(T(I+IHE(J),J+1,L+1)-DPT-T(I,J,L))

! Chou: including specific humidity and liquid water contributions
                        WDUQD=wdep*SIGN(1.0,  -STNE)*(Q(I+IHE(J),J+1,L+1)-Q(I,J,L))
                        WDUCWMD=wdep*SIGN(1.0,-STNE)*(CWM(I+IHE(J),J+1,L+1)-CWM(I,J,L))

                        IF (  -STNE > 0.0) THEN
                            warr=DPDE(I,J)
                            ATC(I,J)=ATC(I,J)+WDUTD/(warr+wdep)
                            AQC(I,J)=AQC(I,J)+WDUQD/(warr+wdep)
                            ACWMC(I,J)=ACWMC(I,J)+WDUCWMD/(warr+wdep)
                        ELSE
                            warr=DPDEP1(I+IHE(J),J+1)
                            ATCP1(I+IHE(J),J+1)=ATCP1(I+IHE(J),J+1)+WDUTD/(warr+wdep)
                            AQCP1(I+IHE(J),J+1)=AQCP1(I+IHE(J),J+1)+WDUQD/(warr+wdep)
                            ACWMCP1(I+IHE(J),J+1)=ACWMCP1(I+IHE(J),J+1)+WDUCWMD/(warr+wdep)
                        END IF

                        DIV(I+IHE(J),J-1,L+1) = DIV(I+IHE(J),J-1,L+1) - STSE
                        DIVS(I,J)             = DIVS(I,J)             + STSE
                        wdep=ABS(   -STSE)*DT
                        DPT=(RTOP(I+IHE(J),J-1,L+1)+RTOP(I,J,L))* &
                        (DPDEP1(I+IHE(J),J-1)    +DPDE(I,J))  *RFCP
                        WDUTD=wdep*SIGN(1.0,   -STSE)*(T(I+IHE(J),J-1,L+1)-DPT-T(I,J,L))

! Chou: including specific humidity and liquid water contributions
                        WDUQD=wdep*SIGN(1.0,   -STSE)*(Q(I+IHE(J),J-1,L+1)-Q(I,J,L))
                        WDUCWMD=wdep*SIGN(1.0, -STSE)*(CWM(I+IHE(J),J-1,L+1)-CWM(I,J,L))

                        IF (   -STSE > 0.0) THEN
                            warr=DPDE(I,J)
                            ATC(I,J)=ATC(I,J)+WDUTD/(warr+wdep)
                            AQC(I,J)=AQC(I,J)+WDUQD/(warr+wdep)
                            ACWMC(I,J)=ACWMC(I,J)+WDUCWMD/(warr+wdep)
                        ELSE
                            warr=DPDEP1(I+IHE(J),J-1)
                            ATCP1(I+IHE(J),J-1)=ATCP1(I+IHE(J),J-1)+WDUTD/(warr+wdep)
                            AQCP1(I+IHE(J),J-1)=AQCP1(I+IHE(J),J-1)+WDUQD/(warr+wdep)
                            ACWMCP1(I+IHE(J),J-1)=ACWMCP1(I+IHE(J),J-1)+WDUCWMD/(warr+wdep)
                        END IF

                        DIV(I+1,J,L+1)        = DIV(I+1,J,L+1)        - STEW  !green
                        DIVS(I,J)             = DIVS(I,J)             + STEW  !green
                        wdep=ABS(  -STEW)*DT
                        DPT=(RTOP(I+1,J,L+1)+RTOP(I,J,L))* &
                        (DPDEP1(I+1,J)    +DPDE(I,J))  *RFCP
                        WDUTD=wdep*SIGN(1.0,  -STEW)*(T(I+1,J,L+1)-DPT-T(I,J,L))

! Chou: including specific humidity and liquid water contributions
                        WDUQD=wdep*SIGN(1.0,  -STEW)*(Q(I+1,J,L+1)-Q(I,J,L))
                        WDUCWMD=wdep*SIGN(1.0,-STEW)*(CWM(I+1,J,L+1)-CWM(I,J,L))

                        IF (  -STEW > 0.0) THEN
                            warr=DPDE(I,J)
                            ATC(I,J)=ATC(I,J)+WDUTD/(warr+wdep)
                            AQC(I,J)=AQC(I,J)+WDUQD/(warr+wdep)
                            ACWMC(I,J)=ACWMC(I,J)+WDUCWMD/(warr+wdep)
                        ELSE
                            warr=DPDEP1(I+1,J)
                            ATCP1(I+1,J)=ATCP1(I+1,J)+WDUTD/(warr+wdep)
                            AQCP1(I+1,J)=AQCP1(I+1,J)+WDUQD/(warr+wdep)
                            ACWMCP1(I+1,J)=ACWMCP1(I+1,J)+WDUCWMD/(warr+wdep)
                        END IF

                    ELSE IF (ISLD(I+IHE(J),J) == 6) THEN
                        DIV(I+IHE(J),J+1,L+1) = DIV(I+IHE(J),J+1,L+1) - STNE
                        DIVS(I,J)             = DIVS(I,J)             + STNE
                        wdep=ABS(  -STNE)*DT
                        DPT=(RTOP(I+IHE(J),J+1,L+1)+RTOP(I,J,L))* &
                        (DPDEP1(I+IHE(J),J+1)    +DPDE(I,J))  *RFCP
                        WDUTD=wdep*SIGN(1.0,  -STNE)*(T(I+IHE(J),J+1,L+1)-DPT-T(I,J,L))

! Chou: including specific humidity and liquid water contributions
                        WDUQD=wdep*SIGN(1.0,  -STNE)*(Q(I+IHE(J),J+1,L+1)-Q(I,J,L))
                        WDUCWMD=wdep*SIGN(1.0,-STNE)*(CWM(I+IHE(J),J+1,L+1)-CWM(I,J,L))

                        IF (  -STNE > 0.0) THEN
                            warr=DPDE(I,J)
                            ATC(I,J)=ATC(I,J)+WDUTD/(warr+wdep)
                            AQC(I,J)=AQC(I,J)+WDUQD/(warr+wdep)
                            ACWMC(I,J)=ACWMC(I,J)+WDUCWMD/(warr+wdep)
                        ELSE
                            warr=DPDEP1(I+IHE(J),J+1)
                            ATCP1(I+IHE(J),J+1)=ATCP1(I+IHE(J),J+1)+WDUTD/(warr+wdep)
                            AQCP1(I+IHE(J),J+1)=AQCP1(I+IHE(J),J+1)+WDUQD/(warr+wdep)
                            ACWMCP1(I+IHE(J),J+1)=ACWMCP1(I+IHE(J),J+1)+WDUCWMD/(warr+wdep)
                        END IF

                        DIV(I+1,J,L+1)        = DIV(I+1,J,L+1)        - STNEM1
                        DIVS(I+IHE(J),J-1)    = DIVS(I+IHE(J),J-1)    + STNEM1
                        wdep=ABS(-STNEM1)*DT
                        DPT=(RTOP(I+1,J,L+1)+RTOP(I+IHE(J),J-1,L))* &
                        (DPDEP1(I+1,J)    +DPDE(I+IHE(J),J-1))  *RFCP
                        WDUTD=wdep*SIGN(1.0,-STNEM1)*(T(I+1,J,L+1)-DPT-T(I+IHE(J),J-1,L))

! Chou: including specific humidity and liquid water contributions
                        WDUQD=wdep*SIGN(1.0,-STNEM1)*(Q(I+1,J,L+1)-Q(I+IHE(J),J-1,L))
                        WDUCWMD=wdep*SIGN(1.0,-STNEM1)*(CWM(I+1,J,L+1)-CWM(I+IHE(J),J-1,L))

                        IF (-STNEM1 > 0.0) THEN
                            warr=DPDE(I+IHE(J),J-1)
                            ATC(I+IHE(J),J-1)=ATC(I+IHE(J),J-1)+WDUTD/(warr+wdep)
                            AQC(I+IHE(J),J-1)=AQC(I+IHE(J),J-1)+WDUQD/(warr+wdep)
                            ACWMC(I+IHE(J),J-1)=ACWMC(I+IHE(J),J-1)+WDUCWMD/(warr+wdep)
                        ELSE
                            warr=DPDEP1(I+1,J)
                            ATCP1(I+1,J)=ATCP1(I+1,J)+WDUTD/(warr+wdep)
                            AQCP1(I+1,J)=AQCP1(I+1,J)+WDUQD/(warr+wdep)
                            ACWMCP1(I+1,J)=ACWMCP1(I+1,J)+WDUCWMD/(warr+wdep)
                        END IF

                        DIV(I+1,J,L+1)        = DIV(I+1,J,L+1)        - STEW  !green
                        DIVS(I,J)             = DIVS(I,J)             + STEW  !green
                        wdep=ABS(  -STEW)*DT
                        DPT=(RTOP(I+1,J,L+1)+RTOP(I,J,L))* &
                        (DPDEP1(I+1,J)    +DPDE(I,J))  *RFCP
                        WDUTD=wdep*SIGN(1.0,  -STEW)*(T(I+1,J,L+1)-DPT-T(I,J,L))

! Chou: including specific humidity and liquid water contributions
                        WDUQD=wdep*SIGN(1.0,  -STEW)*(Q(I+1,J,L+1)-Q(I,J,L))
                        WDUCWMD=wdep*SIGN(1.0,-STEW)*(CWM(I+1,J,L+1)-CWM(I,J,L))

                        IF (  -STEW > 0.0) THEN
                            warr=DPDE(I,J)
                            ATC(I,J)=ATC(I,J)+WDUTD/(warr+wdep)
                            AQC(I,J)=AQC(I,J)+WDUQD/(warr+wdep)
                            ACWMC(I,J)=ACWMC(I,J)+WDUCWMD/(warr+wdep)
                        ELSE
                            warr=DPDEP1(I+1,J)
                            ATCP1(I+1,J)=ATCP1(I+1,J)+WDUTD/(warr+wdep)
                            AQCP1(I+1,J)=AQCP1(I+1,J)+WDUQD/(warr+wdep)
                            ACWMCP1(I+1,J)=ACWMCP1(I+1,J)+WDUCWMD/(warr+wdep)
                        END IF

                        DIV(I+IHE(J),J+1,L+1) = DIV(I+IHE(J),J+1,L+1) - STNS  !pink
                        DIVS(I+IHE(J),J-1)    = DIVS(I+IHE(J),J-1)    + STNS  !pink
                        wdep=ABS(  -STNS)*DT
                        DPT=(RTOP(I+IHE(J),J+1,L+1)+RTOP(I+IHE(J),J-1,L))* &
                        (DPDEP1(I+IHE(J),J+1)    +DPDE(I+IHE(J),J-1))  *RFCP
                        WDUTD=wdep*SIGN(1.0,  -STNS) &
                        *(T(I+IHE(J),J+1,L+1)-DPT-T(I+IHE(J),J-1,L))

! Chou: including specific humidity and liquid water contributions
                        WDUQD=wdep*SIGN(1.0,  -STNS)*(Q(I+IHE(J),J+1,L+1)-Q(I+IHE(J),J-1,L))
                        WDUCWMD=wdep*SIGN(1.0,-STNS)*(CWM(I+IHE(J),J+1,L+1)-CWM(I+IHE(J),J-1,L))

                        IF (  -STNS > 0.0) THEN
                            warr=DPDE(I+IHE(J),J-1)
                            ATC(I+IHE(J),J-1)=ATC(I+IHE(J),J-1)+WDUTD/(warr+wdep)
                            AQC(I+IHE(J),J-1)=AQC(I+IHE(J),J-1)+WDUQD/(warr+wdep)
                            ACWMC(I+IHE(J),J-1)=ACWMC(I+IHE(J),J-1)+WDUCWMD/(warr+wdep)
                        ELSE
                            warr=DPDEP1(I+IHE(J),J+1)
                            ATCP1(I+IHE(J),J+1)=ATCP1(I+IHE(J),J+1)+WDUTD/(warr+wdep)
                            AQCP1(I+IHE(J),J+1)=AQCP1(I+IHE(J),J+1)+WDUQD/(warr+wdep)
                            ACWMCP1(I+IHE(J),J+1)=ACWMCP1(I+IHE(J),J+1)+WDUCWMD/(warr+wdep)
                        END IF

                    ELSE IF (ISLD(I+IHE(J),J) == 7) THEN
                        DIV(I,J,L+1)          = DIV(I,J,L+1)          + STSE
                        DIVS(I+IHE(J),J-1)    = DIVS(I+IHE(J),J-1)    - STSE
                        wdep=ABS(   STSE)*DT
                        DPT=(RTOP(I,J,L+1)+RTOP(I+IHE(J),J-1,L))* &
                        (DPDEP1(I,J)    +DPDE(I+IHE(J),J-1))  *RFCP
                        WDUTD=wdep*SIGN(1.0,   STSE)*(T(I,J,L+1)-DPT-T(I+IHE(J),J-1,L))

! Chou: including specific humidity and liquid water contributions
                        WDUQD=wdep*SIGN(1.0,   STSE)*(Q(I,J,L+1)-Q(I+IHE(J),J-1,L))
                        WDUCWMD=wdep*SIGN(1.0, STSE)*(CWM(I,J,L+1)-CWM(I+IHE(J),J-1,L))

                        IF (   STSE > 0.0) THEN
                            warr=DPDE(I+IHE(J),J-1)
                            ATC(I+IHE(J),J-1)=ATC(I+IHE(J),J-1)+WDUTD/(warr+wdep)
                            AQC(I+IHE(J),J-1)=AQC(I+IHE(J),J-1)+WDUQD/(warr+wdep)
                            ACWMC(I+IHE(J),J-1)=ACWMC(I+IHE(J),J-1)+WDUCWMD/(warr+wdep)
                        ELSE
                            warr=DPDEP1(I,J)
                            ATCP1(I,J)=ATCP1(I,J)+WDUTD/(warr+wdep)
                            AQCP1(I,J)=AQCP1(I,J)+WDUQD/(warr+wdep)
                            ACWMCP1(I,J)=ACWMCP1(I,J)+WDUCWMD/(warr+wdep)
                        END IF

                        DIV(I+1,J,L+1)        = DIV(I+1,J,L+1)        - STNEM1
                        DIVS(I+IHE(J),J-1)    = DIVS(I+IHE(J),J-1)    + STNEM1
                        wdep=ABS(-STNEM1)*DT
                        DPT=(RTOP(I+1,J,L+1)+RTOP(I+IHE(J),J-1,L))* &
                        (DPDEP1(I+1,J)    +DPDE(I+IHE(J),J-1))  *RFCP
                        WDUTD=wdep*SIGN(1.0,-STNEM1)*(T(I+1,J,L+1)-DPT-T(I+IHE(J),J-1,L))

! Chou: including specific humidity and liquid water contributions
                        WDUQD=wdep*SIGN(1.0,  -STNEM1)*(Q(I+1,J,L+1)-Q(I+IHE(J),J-1,L))
                        WDUCWMD=wdep*SIGN(1.0,-STNEM1)*(CWM(I+1,J,L+1)-CWM(I+IHE(J),J-1,L))

                        IF (-STNEM1 > 0.0) THEN
                            warr=DPDE(I+IHE(J),J-1)
                            ATC(I+IHE(J),J-1)=ATC(I+IHE(J),J-1)+WDUTD/(warr+wdep)
                            AQC(I+IHE(J),J-1)=AQC(I+IHE(J),J-1)+WDUQD/(warr+wdep)
                            ACWMC(I+IHE(J),J-1)=ACWMC(I+IHE(J),J-1)+WDUCWMD/(warr+wdep)
                        ELSE
                            warr=DPDEP1(I+1,J)
                            ATCP1(I+1,J)=ATCP1(I+1,J)+WDUTD/(warr+wdep)
                            AQCP1(I+1,J)=AQCP1(I+1,J)+WDUQD/(warr+wdep)
                            ACWMCP1(I+1,J)=ACWMCP1(I+1,J)+WDUCWMD/(warr+wdep)
                        END IF

                        DIV(I+IHE(J),J+1,L+1) = DIV(I+IHE(J),J+1,L+1) - STNS  !pink
                        DIVS(I+IHE(J),J-1)    = DIVS(I+IHE(J),J-1)    + STNS  !pink
                        wdep=ABS(  -STNS)*DT
                        DPT=(RTOP(I+IHE(J),J+1,L+1)+RTOP(I+IHE(J),J-1,L))* &
                        (DPDEP1(I+IHE(J),J+1)    +DPDE(I+IHE(J),J-1))  *RFCP
                        WDUTD=wdep*SIGN(1.0,  -STNS) &
                        *(T(I+IHE(J),J+1,L+1)-DPT-T(I+IHE(J),J-1,L))

! Chou: including specific humidity and liquid water contributions
                        WDUQD=wdep*SIGN(1.0,  -STNS) *(Q(I+IHE(J),J+1,L+1)-Q(I+IHE(J),J-1,L))
                        WDUCWMD=wdep*SIGN(1.0,-STNS) *(CWM(I+IHE(J),J+1,L+1)-CWM(I+IHE(J),J-1,L))

                        IF (  -STNS > 0.0) THEN
                            warr=DPDE(I+IHE(J),J-1)
                            ATC(I+IHE(J),J-1)=ATC(I+IHE(J),J-1)+WDUTD/(warr+wdep)
                            AQC(I+IHE(J),J-1)=AQC(I+IHE(J),J-1)+WDUQD/(warr+wdep)
                            ACWMC(I+IHE(J),J-1)=ACWMC(I+IHE(J),J-1)+WDUCWMD/(warr+wdep)
                        ELSE
                            warr=DPDEP1(I+IHE(J),J+1)
                            ATCP1(I+IHE(J),J+1)=ATCP1(I+IHE(J),J+1)+WDUTD/(warr+wdep)
                            AQCP1(I+IHE(J),J+1)=AQCP1(I+IHE(J),J+1)+WDUQD/(warr+wdep)
                            ACWMCP1(I+IHE(J),J+1)=ACWMCP1(I+IHE(J),J+1)+WDUCWMD/(warr+wdep)
                        END IF

                    ELSE IF (ISLD(I+IHE(J),J) == 8) THEN
                        DIV(I+IHE(J),J+1,L+1) = DIV(I+IHE(J),J+1,L+1) + STSEP1
                        DIVS(I+1,J)           = DIVS(I+1,J)           - STSEP1
                        wdep=ABS( STSEP1)*DT
                        DPT=(RTOP(I+IHE(J),J+1,L+1)+RTOP(I+1,J,L))* &
                        (DPDEP1(I+IHE(J),J+1)    +DPDE(I+1,J))  *RFCP
                        WDUTD=wdep*SIGN(1.0, STSEP1) &
                        *(T(I+IHE(J),J+1,L+1)-DPT-T(I+1,J,L))

! Chou: including specific humidity and liquid water contributions
                        WDUQD=wdep*SIGN(1.0,  STSEP1)*(Q(I+IHE(J),J+1,L+1)-Q(I+1,J,L))
                        WDUCWMD=wdep*SIGN(1.0,STSEP1)*(CWM(I+IHE(J),J+1,L+1)-CWM(I+1,J,L))

                        IF ( STSEP1 > 0.0) THEN
                            warr=DPDE(I+1,J)
                            ATC(I+1,J)=ATC(I+1,J)+WDUTD/(warr+wdep)
                            AQC(I+1,J)=AQC(I+1,J)+WDUQD/(warr+wdep)
                            ACWMC(I+1,J)=ACWMC(I+1,J)+WDUCWMD/(warr+wdep)
                        ELSE
                            warr=DPDEP1(I+IHE(J),J+1)
                            ATCP1(I+IHE(J),J+1)=ATCP1(I+IHE(J),J+1)+WDUTD/(warr+wdep)
                            AQCP1(I+IHE(J),J+1)=AQCP1(I+IHE(J),J+1)+WDUQD/(warr+wdep)
                            ACWMCP1(I+IHE(J),J+1)=ACWMCP1(I+IHE(J),J+1)+WDUCWMD/(warr+wdep)
                        END IF
                        		   
                        DIV(I,J,L+1)          = DIV(I,J,L+1)          + STSE
                        DIVS(I+IHE(J),J-1)    = DIVS(I+IHE(J),J-1)    - STSE
                        wdep=ABS(   STSE)*DT
                        DPT=(RTOP(I,J,L+1)+RTOP(I+IHE(J),J-1,L))* &
                        (DPDEP1(I,J)    +DPDE(I+IHE(J),J-1))  *RFCP
                        WDUTD=wdep*SIGN(1.0,   STSE)*(T(I,J,L+1)-DPT-T(I+IHE(J),J-1,L))

! Chou: including specific humidity and liquid water contributions
                        WDUQD=wdep*SIGN(1.0,   STSE)*(Q(I,J,L+1)-Q(I+IHE(J),J-1,L))
                        WDUCWMD=wdep*SIGN(1.0, STSE)*(CWM(I,J,L+1)-CWM(I+IHE(J),J-1,L))

                        IF (   STSE > 0.0) THEN
                            warr=DPDE(I+IHE(J),J-1)
                            ATC(I+IHE(J),J-1)=ATC(I+IHE(J),J-1)+WDUTD/(warr+wdep)
                            AQC(I+IHE(J),J-1)=AQC(I+IHE(J),J-1)+WDUQD/(warr+wdep)
                            ACWMC(I+IHE(J),J-1)=ACWMC(I+IHE(J),J-1)+WDUCWMD/(warr+wdep)
                        ELSE
                            warr=DPDEP1(I,J)
                            ATCP1(I,J)=ATCP1(I,J)+WDUTD/(warr+wdep)
                            AQCP1(I,J)=AQCP1(I,J)+WDUQD/(warr+wdep)
                            ACWMCP1(I,J)=ACWMCP1(I,J)+WDUCWMD/(warr+wdep)
                        END IF

                        DIV(I,J,L+1)          = DIV(I,J,L+1)          + STEW  !blue
                        DIVS(I+1,J)           = DIVS(I+1,J)           - STEW  !blue
                        wdep=ABS(   STEW)*DT
                        DPT=(RTOP(I,J,L+1)+RTOP(I+1,J,L))* &
                        (DPDEP1(I,J)    +DPDE(I+1,J))  *RFCP
                        WDUTD=wdep*SIGN(1.0,   STEW)*(T(I,J,L+1)-DPT-T(I+1,J,L))

! Chou: including specific humidity and liquid water contributions
                        WDUQD=wdep*SIGN(1.0,   STEW)*(Q(I,J,L+1)-Q(I+1,J,L))
                        WDUCWMD=wdep*SIGN(1.0, STEW)*(CWM(I,J,L+1)-CWM(I+1,J,L))

                        IF (   STEW > 0.0) THEN
                            warr=DPDE(I+1,J)
                            ATC(I+1,J)=ATC(I+1,J)+WDUTD/(warr+wdep)
                            AQC(I+1,J)=AQC(I+1,J)+WDUQD/(warr+wdep)
                            ACWMC(I+1,J)=ACWMC(I+1,J)+WDUCWMD/(warr+wdep)
                        ELSE
                            warr=DPDEP1(I,J)
                            ATCP1(I,J)=ATCP1(I,J)+WDUTD/(warr+wdep)
                            AQCP1(I,J)=AQCP1(I,J)+WDUQD/(warr+wdep)
                            ACWMCP1(I,J)=ACWMCP1(I,J)+WDUCWMD/(warr+wdep)
                        END IF

                        DIV(I+IHE(J),J+1,L+1) = DIV(I+IHE(J),J+1,L+1) - STNS  !pink
                        DIVS(I+IHE(J),J-1)    = DIVS(I+IHE(J),J-1)    + STNS  !pink
                        wdep=ABS(  -STNS)*DT
                        DPT=(RTOP(I+IHE(J),J+1,L+1)+RTOP(I+IHE(J),J-1,L))* &
                        (DPDEP1(I+IHE(J),J+1)    +DPDE(I+IHE(J),J-1))  *RFCP
                        WDUTD=wdep*SIGN(1.0,  -STNS) &
                        *(T(I+IHE(J),J+1,L+1)-DPT-T(I+IHE(J),J-1,L))

! Chou: including specific humidity and liquid water contributions
                        WDUQD=wdep*SIGN(1.0,  -STNS)*(Q(I+IHE(J),J+1,L+1)-Q(I+IHE(J),J-1,L))
                        WDUCWMD=wdep*SIGN(1.0,-STNS)*(CWM(I+IHE(J),J+1,L+1)-CWM(I+IHE(J),J-1,L))

                        IF (  -STNS > 0.0) THEN
                            warr=DPDE(I+IHE(J),J-1)
                            ATC(I+IHE(J),J-1)=ATC(I+IHE(J),J-1)+WDUTD/(warr+wdep)
                            AQC(I+IHE(J),J-1)=AQC(I+IHE(J),J-1)+WDUQD/(warr+wdep)
                            ACWMC(I+IHE(J),J-1)=ACWMC(I+IHE(J),J-1)+WDUCWMD/(warr+wdep)
                        ELSE
                            warr=DPDEP1(I+IHE(J),J+1)
                            ATCP1(I+IHE(J),J+1)=ATCP1(I+IHE(J),J+1)+WDUTD/(warr+wdep)
                            AQCP1(I+IHE(J),J+1)=AQCP1(I+IHE(J),J+1)+WDUQD/(warr+wdep)
                            ACWMCP1(I+IHE(J),J+1)=ACWMCP1(I+IHE(J),J+1)+WDUCWMD/(warr+wdep)
                        END IF
                    END IF

                END IF !! L < LM
            
        340 END DO
    
    !------------------------------------------------------------------------
    
    ! omp parallel do
        DO 390 J=MYJS,MYJE
            DO 390 I=MYIS,MYIE
                DIV(I,J,L)=(DIV(I,J,L)+DIVL(I,J))*HM(I,J)
        390 END DO
    !-----------------------------------------------------------------------
    400 END DO
!-----------------------------------------------------------------------
    RETURN
    END SUBROUTINE DIVHOASTQL
