                             SUBROUTINE DIVHOAST
C     ******************************************************************
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .
C SUBPROGRAM:    DIVHOA      DIVERGENCE/HORIZONTAL OMEGA-ALPHA
C   PRGRMMR: JANJIC          ORG: W/NP22     DATE: 93-10-28
C
C ABSTRACT:
C     DIVHOA COMPUTES THE DIVERGENCE INCLUDING THE
C     MODIFICATION PREVENTING GRAVITY WAVE GRID SEPARATION, AND
C     CALCULATES THE HORIZONTAL PART OF THE OMEGA-ALPHA TERM
C     (THE PART PROPORTIONAL TO THE ADVECTION OF MASS ALONG
C     ETA/SIGMA SURFACES).
C  
C     MODIFIED TO INCLUDE DIVERGENCE RESULTING FROM SLOPING STEPS 
C     ("S": DIVHOAST)
C     EXPANDED to do also Slantwise T Advection ("T": DIVHOAST), along
C     with appropriate omega-alpha changes to the advected T.  
C     Warning: this makes Subroutine SLADVT redundant
C        (needs to be commented out or removed)
C***********************************************************************
C PROGRAM HISTORY LOG:
C   87-06-??  JANJIC     - ORIGINATOR
C   95-03-25  BLACK      - CONVERSION FROM 1-D TO 2-D IN HORIZONTAL
C   96-03-29  BLACK      - ADDED EXTERNAL EDGE
C   97-03-17  MESINGER   - SPLIT FROM PFDHT
C   98-10-30  BLACK      - MODIFIED FOR DISTRIBUTED MEMORY
C   00-10-20  BLACK      - INCORPORATED PRESSURE GRADIENT METHOD
C                          FROM MESO MODEL
C
C   04  MESINGER and JOVIC - ADDED DIVERGENCE CALCULATION DUE TO SLOPES
C   May 2006  MESINGER, at CPTEC - EXPANDED TO INCLUDE SLANTWISE T ADV.
C   June 2006  MESINGER and JOVIC - REVISED OMEGA-ALPHA DUE TO SLOPES 
C
C USAGE: CALL DIVHOA FROM MAIN PROGRAM EBU
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: NONE
C
C     LIBRARY: NONE
C
C   COMMON BLOCKS: CTLBLK
C                  MASKS
C                  LOOPS
C                  DYNAM
C                  VRBLS
C                  CONTIN
C                  INDX
C   
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 90
C   MACHINE : IBM SP
C$$$  
C***********************************************************************
C-----------------------------------------------------------------------
      INCLUDE "parmeta"
      INCLUDE "mpp.h"
C-----------------------------------------------------------------------
                             P A R A M E T E R
     & (LP1=LM+1,JAM=6+2*(JM-10),RFCP=1.E0/(4.E0*1004.6E0))
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 "INDX.comm"
c-----------------------------------------------------------------------
      INCLUDE "DYNAM.comm"
C-----------------------------------------------------------------------
      INCLUDE "VRBLS.comm"
C-----------------------------------------------------------------------
      INCLUDE "CONTIN.comm"
C-----------------------------------------------------------------------
      INCLUDE "NHYDRO.comm"
C-----------------------------------------------------------------------
      INCLUDE "SLOPES.comm"
C-----------------------------------------------------------------------
                             R E A L
     & PINTLG(IDIM1:IDIM2,JDIM1:JDIM2,LM+1)
C
                             R E A L
     & 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)
C
                             R E A L
     & 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)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
      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)
C-----------------------------------------------------------------------
C--------------PREPARATORY CALCULATIONS---------------------------------
C-----------------------------------------------------------------------
      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   CONTINUE
      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   CONTINUE
      ENDIF
C
      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
C
!$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.
C   DIVergence to Save is needed to collect contributions to DIV due to
C   slopes at h points of the layer being processed (L).  They can be
C   added only once the code moves one layer up.  The reason is that the 
C   collection is done for points also sideways of the h point that the
C   code is processing, so that if added immediatelly these would be 
C   lost when the code visits these sideways located points, and puts
C   the divergence "correction" into DIV (loop 270)

C   Accumulated T changes at L and at L+1
          ATC  (I,J)=0.0
          ATCP1(I,J)=0.0
        ENDDO
        ENDDO

C-----------------------------------------------------------------------
!$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
C-----------------------------------------------------------------------
C-------------------- MAIN VERTICAL INTEGRATION LOOP -------------------
C-----------------------------------------------------------------------
                             DO 400 L=LM,1,-1
C-----------------------------------------------------------------------
C***
C***  INTEGRATE THE GEOPOTENTIAL
C***
Cmp
	FIM=0.
Cmp
C
!$omp parallel do private (alp1p,dfi,fiupk,rdpds)
C	write(6,*) 'FIM defined over I: ', MYIS_P5,MYIE_P5
C	write(6,*) 'FIM defined over J: ', MYJS_P5,MYJE_P5
      DO 125 J=MYJS_P5,MYJE_P5
      DO 125 I=MYIS_P5,MYIE_P5
C
      ALP1P=PINTLG(I,J,L)
C
      DFI=(Q(I,J,L)*0.608+1.)*T(I,J,L)*R*(ALP1(I,J)-ALP1P)/DWDT(I,J,L)
C
      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)) .le. 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
C
      FILO(I,J)=(FIUPK-DFL(L))*HTM(I,J,L)+DFL(L)
      ALP1(I,J)=ALP1P
  125 CONTINUE
C
C-----------------------------------------------------------------------
!$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)
C
      ALP1P =PINTLG(I,J,L)
      ALP1PL=PINTLG(I,J,L+1)
      ALP2P =ALP1P*ALP1P
      ALP2PL=ALP1PL*ALP1PL
C
      DFI=(Q(I,J,L)*0.608+1.)*T(I,J,L)*R*(ALP1PL-ALP1P)/DWDT(I,J,L)
      DFDZ(I,J)=DFI*DWDT(I,J,L)/(ALP2PL-ALP2P)
      APEL(I,J)=(ALP2PL+ALP2P)*0.5
  205 CONTINUE
C
!$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 CONTINUE
C
         IF (L.LT.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 CONTINUE
         END IF
C
!$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 CONTINUE
C
!$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 CONTINUE
C
C--------------DIAGONAL CONTRIBUTIONS TO PRESSURE GRADIENT FORCE--------
C
!$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))
     1        *(DWDT(I+IHE(J),J+1,L)+DWDT(I,J,L))
	if ( ABS(PNE(I,J)) .le. 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.
     1        *(APEL(I+IHE(J),J+1)-APEL(I,J))
      PCNE(I,J)=CNE(I,J)*ADPDNE(I,J)
  240 CONTINUE
C
!$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))
     1        *(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.
     1        *(APEL(I+IHE(J),J-1)-APEL(I,J))
      PCSE(I,J)=CSE(I,J)*ADPDSE(I,J)
  250 CONTINUE
C
C--------------CONTINUITY EQUATION MODIFICATION-------------------------
C
!$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)
     1         +CNE(I+IVW(J),J)+PSE(I+IVW(J),J)+CSE(I+IVW(J),J)
     2         -PNE(I,J-1)-CNE(I,J-1)-PSE(I,J+1)-CSE(I,J+1))
  260 CONTINUE
C-----------------------------------------------------------------------
      DO 270 J=MYJS2,MYJE2
      DO 270 I=MYIS1,MYIE1
      DIV(I,J,L)=DETA(L)*WPDAR(I,J)
     1         *(PCXC(I+IHE(J),J)-PCXC(I,J+1)
     2          +PCXC(I+IHW(J),J)-PCXC(I,J-1))
  270 CONTINUE
C
C--------------LAT & LONG PRESSURE FORCE COMPONENTS---------------------
C
!$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 CONTINUE
C
C--------------LAT & LON FLUXES & OMEGA-ALPHA COMPONENTS----------------
C
!$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 CONTINUE
C
C--------------DIAGONAL FLUXES AND DIAGONALLY AVERAGED WIND-------------
C
!$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 CONTINUE
C
!$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 CONTINUE
         IF (L.LT.LM-1) THEN  
      DO 335 J=MYJS2_P1,MYJE2_P1
      DO 335 I=MYIS1_P1,MYIE1_P1
         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)
         ATCP1(I,J)=ATC(I,J)
         ATC  (I,J)=0.0
  335 CONTINUE
         END IF 
C
C--------------HORIZONTAL PART OF OMEGA-ALPHA & DIVERGENCE--------------
C
!$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)
     1              +TNS(I,J-1)+TNE(I,J)+TNE(I+IHW(J),J-1)+TSE(I,J)
     2              +TSE(I+IHW(J),J+1))*RDPD(I,J)*FCP(I,J)*HM(I,J)
         IF (L.LT.3) THEN
      T(I,J,L)=OMGALF(I,J,L)+T(I,J,L)
         ELSE
      ATC(I,J)=ATC(I,J)+OMGALF(I,J,L)
         END IF
C
      IF (L.LT.LM .AND. (HTM(I,J,L+1)*HBM2(I,J)).GT.0) THEN
C   The T point below is in the atmosphere.  Should it receive
C   omega-alpha contributions from fluxes of layer L?

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

         TOASL=0.0

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

         IF (VTMS(I,J+1,L+1).EQ.1) THEN
            IF (ISLD(I,J+1).EQ.3)
     &         TOASL=TOASL+0.50*TNS(I,J+1)
            IF (ISLD(I,J+1).EQ.2)
     &         TOASL=TOASL+0.50*TNS(I,J+1)+0.25*TNE(I,J)
            IF (ISLD(I,J+1).EQ.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).EQ.1) THEN
            IF (ISLD(I+IHW(J),J).EQ.5)
     &         TOASL=TOASL+0.50*TEW(I+IHW(J),J)
            IF (ISLD(I+IHW(J),J).EQ.4)
     &         TOASL=TOASL+0.50*TEW(I+IHW(J),J)+0.25*TSE(I+IHW(J),J+1)
            IF (ISLD(I+IHW(J),J).EQ.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).EQ.1) THEN
            IF (ISLD(I,J-1).EQ.7)
     &         TOASL=TOASL+0.50*TNS(I,J-1)
            IF (ISLD(I,J-1).EQ.6)
     &         TOASL=TOASL+0.50*TNS(I,J-1)+0.25*TNE(I+IHW(J),J-1)
            IF (ISLD(I,J-1).EQ.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.LT.3) THEN
      T(I,J,L+1)=T(I,J,L+1) +OAADP1
         ELSE
      ATCP1(I,J)=ATCP1(I,J) +OAADP1
         END IF

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

      IF (L.LT.LM .AND. L.EQ.LMV(I+IHE(J),J)
     1   .AND. ISLD(I+IHE(J),J).GT.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).EQ.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))
           IF ( STSEP1.GT.0.0) THEN
        warr=DPDE(I+1,J)
        ATC(I+1,J)=ATC(I+1,J)+WDUTD/(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)
           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))
           IF ( STNEM1.GT.0.0) THEN
        warr=DPDE(I+1,J)
        ATC(I+1,J)=ATC(I+1,J)+WDUTD/(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)
           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))
           IF (   STEW.GT.0.0) THEN
        warr=DPDE(I+1,J)
        ATC(I+1,J)=ATC(I+1,J)+WDUTD/(warr+wdep)
           ELSE
        warr=DPDEP1(I,J)
        ATCP1(I,J)=ATCP1(I,J)+WDUTD/(warr+wdep)
           END IF

        ELSE IF (ISLD(I+IHE(J),J).EQ.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))
           IF (   STNE.GT.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)
           ELSE
        warr=DPDEP1(I,J)
        ATCP1(I,J)=ATCP1(I,J)+WDUTD/(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))
           IF ( STNEM1.GT.0.0) THEN
        warr=DPDE(I+1,J)
        ATC(I+1,J)=ATC(I+1,J)+WDUTD/(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)
           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))
           IF (   STEW.GT.0.0) THEN
        warr=DPDE(I+1,J)
        ATC(I+1,J)=ATC(I+1,J)+WDUTD/(warr+wdep)
           ELSE
        warr=DPDEP1(I,J)
        ATCP1(I,J)=ATCP1(I,J)+WDUTD/(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))
           IF (  STNS.GT.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)
           ELSE
        warr=DPDEP1(I+IHE(J),J-1)
        ATCP1(I+IHE(J),J-1)=ATCP1(I+IHE(J),J-1)+WDUTD/(warr+wdep)
           END IF

        ELSE IF (ISLD(I+IHE(J),J).EQ.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))
           IF (-STSEP1.GT.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)
           ELSE
        warr=DPDEP1(I+1,J)
        ATCP1(I+1,J)=ATCP1(I+1,J)+WDUTD/(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))
           IF (   STNE.GT.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)
           ELSE
        warr=DPDEP1(I,J)
        ATCP1(I,J)=ATCP1(I,J)+WDUTD/(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))
           IF (  STNS.GT.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)
           ELSE
        warr=DPDEP1(I+IHE(J),J-1)
        ATCP1(I+IHE(J),J-1)=ATCP1(I+IHE(J),J-1)+WDUTD/(warr+wdep)
           END IF

        ELSE IF (ISLD(I+IHE(J),J).EQ.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))
           IF (-STSEP1.GT.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)
           ELSE
        warr=DPDEP1(I+1,J)
        ATCP1(I+1,J)=ATCP1(I+1,J)+WDUTD/(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))
           IF (   -STSE.GT.0.0) THEN
        warr=DPDE(I,J)
        ATC(I,J)=ATC(I,J)+WDUTD/(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)
           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))
           IF (  STNS.GT.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)
           ELSE
        warr=DPDEP1(I+IHE(J),J-1)
        ATCP1(I+IHE(J),J-1)=ATCP1(I+IHE(J),J-1)+WDUTD/(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))
           IF (  -STEW.GT.0.0) THEN
        warr=DPDE(I,J)
        ATC(I,J)=ATC(I,J)+WDUTD/(warr+wdep)
           ELSE
        warr=DPDEP1(I+1,J)
        ATCP1(I+1,J)=ATCP1(I+1,J)+WDUTD/(warr+wdep)
           END IF

        ELSE IF (ISLD(I+IHE(J),J).EQ.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))
           IF (  -STNE.GT.0.0) THEN
        warr=DPDE(I,J)
        ATC(I,J)=ATC(I,J)+WDUTD/(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)
           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))
           IF (   -STSE.GT.0.0) THEN
        warr=DPDE(I,J)
        ATC(I,J)=ATC(I,J)+WDUTD/(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)
           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))
           IF (  -STEW.GT.0.0) THEN
        warr=DPDE(I,J)
        ATC(I,J)=ATC(I,J)+WDUTD/(warr+wdep)
           ELSE
        warr=DPDEP1(I+1,J)
        ATCP1(I+1,J)=ATCP1(I+1,J)+WDUTD/(warr+wdep)
           END IF

        ELSE IF (ISLD(I+IHE(J),J).EQ.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))
           IF (  -STNE.GT.0.0) THEN
        warr=DPDE(I,J)
        ATC(I,J)=ATC(I,J)+WDUTD/(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)
           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))
           IF (-STNEM1.GT.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)
           ELSE
        warr=DPDEP1(I+1,J)
        ATCP1(I+1,J)=ATCP1(I+1,J)+WDUTD/(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))
           IF (  -STEW.GT.0.0) THEN
        warr=DPDE(I,J)
        ATC(I,J)=ATC(I,J)+WDUTD/(warr+wdep)
           ELSE
        warr=DPDEP1(I+1,J)
        ATCP1(I+1,J)=ATCP1(I+1,J)+WDUTD/(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))
           IF (  -STNS.GT.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)
           ELSE
        warr=DPDEP1(I+IHE(J),J+1)
        ATCP1(I+IHE(J),J+1)=ATCP1(I+IHE(J),J+1)+WDUTD/(warr+wdep)
           END IF

        ELSE IF (ISLD(I+IHE(J),J).EQ.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))
           IF (   STSE.GT.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)
           ELSE
        warr=DPDEP1(I,J)
        ATCP1(I,J)=ATCP1(I,J)+WDUTD/(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))
           IF (-STNEM1.GT.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)
           ELSE
        warr=DPDEP1(I+1,J)
        ATCP1(I+1,J)=ATCP1(I+1,J)+WDUTD/(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))
           IF (  -STNS.GT.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)
           ELSE
        warr=DPDEP1(I+IHE(J),J+1)
        ATCP1(I+IHE(J),J+1)=ATCP1(I+IHE(J),J+1)+WDUTD/(warr+wdep)
           END IF

        ELSE IF (ISLD(I+IHE(J),J).EQ.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))
           IF ( STSEP1.GT.0.0) THEN
        warr=DPDE(I+1,J)
        ATC(I+1,J)=ATC(I+1,J)+WDUTD/(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)
           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))
           IF (   STSE.GT.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)
           ELSE
        warr=DPDEP1(I,J)
        ATCP1(I,J)=ATCP1(I,J)+WDUTD/(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))
           IF (   STEW.GT.0.0) THEN
        warr=DPDE(I+1,J)
        ATC(I+1,J)=ATC(I+1,J)+WDUTD/(warr+wdep)
           ELSE
        warr=DPDEP1(I,J)
        ATCP1(I,J)=ATCP1(I,J)+WDUTD/(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))
           IF (  -STNS.GT.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)
           ELSE
        warr=DPDEP1(I+IHE(J),J+1)
        ATCP1(I+IHE(J),J+1)=ATCP1(I+IHE(J),J+1)+WDUTD/(warr+wdep)
           END IF
        END IF

      END IF !! L.LT.LM
C
  340 CONTINUE
C
C------------------------------------------------------------------------
C
!$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 CONTINUE
C-----------------------------------------------------------------------
  400                        CONTINUE
C-----------------------------------------------------------------------
                             RETURN
                             END                                        
