#include "STREAMICE_OPTIONS.h"

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

CBOP
      SUBROUTINE STREAMICE_DRIVING_STRESS_PPM( myThid ) 
!      O taudx, 
!      O taudy )

C     /============================================================\
C     | SUBROUTINE                                                 |   
C     | o                                                          |
C     |============================================================|
C     |                                                            |
C     \============================================================/
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "STREAMICE.h"
#include "STREAMICE_CG.h"

C     !INPUT/OUTPUT ARGUMENTS
      INTEGER myThid
!       _RL taudx (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
!       _RL taudx (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)

#ifdef ALLOW_STREAMICE


C     LOCAL VARIABLES
      INTEGER i, j, bi, bj, k, l,
     &        Gi, Gj
      LOGICAL at_west_bdry, at_east_bdry, 
     &        at_north_bdry, at_south_bdry
      _RL sx, sy, diffx, diffy, neu_val
      
      IF (myXGlobalLo.eq.1) at_west_bdry = .true.
      IF (myYGlobalLo.eq.1) at_south_bdry = .true.
      IF (myXGlobalLo-1+sNx*nSx.eq.Nx) 
     & at_east_bdry = .false.
      IF (myYGlobalLo-1+sNy*nSy.eq.Ny) 
     & at_north_bdry = .false.

      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          taudx_SI(i,j,bi,bj) = 0. _d 0
          taudy_SI(i,j,bi,bj) = 0. _d 0
         ENDDO
        ENDDO
       ENDDO
      ENDDO

      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        
        DO i=0,sNx+1
         DO j=0,sNy+1

          diffx = 0. _d 0
          diffy = 0. _d 0
          sx = 0. _d 0
          sy = 0. _d 0

          Gi = (myXGlobalLo-1)+(bi-1)*sNx+i
          Gj = (myYGlobalLo-1)+(bj-1)*sNy+j

          IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN

           ! we are in an "active" cell

           IF (Gi.eq.1.AND..NOT.STREAMICE_EW_periodic) THEN

            ! western boundary - only one sided possible

            IF (STREAMICE_hmask(i+1,j,bi,bj).eq.1.0) THEN

             ! cell to east is active

             sx = (surf_el_streamice(i+1,j,bi,bj)-
     &             surf_el_streamice(i,j,bi,bj))/dxC(i+1,j,bi,bj)
            ELSE

             ! cell to east is empty

             sx = 0. _d 0
            ENDIF

            DO k=0,1
             DO l=0,1
              IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
               taudx_SI(i+k,j+l,bi,bj) = taudx_SI(i+k,j+l,bi,bj) -
     &          0.25 * streamice_density * gravity * 
     &          (streamice_bg_surf_slope_x+sx) * 
     &          H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
              ENDIF
             ENDDO
            ENDDO

           ELSEIF (Gi.eq.Nx.AND..NOT.STREAMICE_EW_periodic) THEN

            ! eastern boundary - only one sided possible

            IF (STREAMICE_hmask(i-1,j,bi,bj).eq.1.0) THEN

             ! cell to west is active

             sx = (surf_el_streamice(i,j,bi,bj)-
     &             surf_el_streamice(i-1,j,bi,bj))/dxC(i,j,bi,bj)
            ELSE

             ! cell to west is inactive

             sx = 0. _d 0
            ENDIF

            DO k=0,1
             DO l=0,1
              IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
               taudx_SI(i+k,j+l,bi,bj) = taudx_SI(i+k,j+l,bi,bj) -
     &          0.25 * streamice_density * gravity * 
     &          (streamice_bg_surf_slope_x+sx) * 
     &          H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
              ENDIF
             ENDDO
            ENDDO

           ELSE

            ! interior (west-east) cell

            IF (STREAMICE_hmask(i+1,j,bi,bj).eq.1.0 .and.
     &          STREAMICE_hmask(i-1,j,bi,bj).eq.1.0) THEN

             k = 0
             DO l=0,1
              IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
               taudx_SI(i+k,j+l,bi,bj) = taudx_SI(i+k,j+l,bi,bj) -
     &          streamice_density * gravity * (1./6.) *
     &          (-2.*surf_el_streamice(i-1,j,bi,bj) + 
     &           surf_el_streamice(i,j,bi,bj) + 
     &           surf_el_streamice(i+1,j,bi,bj) +
     &           3.*streamice_bg_surf_slope_x * dxF(i,j,bi,bj)) *
     &          H_streamice(i,j,bi,bj) * .5 * dyF(i,j,bi,bj)
              ENDIF
             ENDDO

             k = 1
             DO l=0,1
              IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
               taudx_SI(i+k,j+l,bi,bj) = taudx_SI(i+k,j+l,bi,bj) -
     &          streamice_density * gravity * (1./6.) *
     &          (-surf_el_streamice(i-1,j,bi,bj) - 
     &           surf_el_streamice(i,j,bi,bj) + 
     &           2*surf_el_streamice(i+1,j,bi,bj) +
     &           3.*streamice_bg_surf_slope_x * dxF(i,j,bi,bj)) *
     &          H_streamice(i,j,bi,bj) * .5 * dyF(i,j,bi,bj)
              ENDIF
             ENDDO


            ELSE 

             IF (STREAMICE_hmask(i+1,j,bi,bj).eq.1.0) THEN

              sx = (surf_el_streamice(i+1,j,bi,bj)-
     &             surf_el_streamice(i,j,bi,bj))/dxC(i+1,j,bi,bj)

             ELSEIF (STREAMICE_hmask(i-1,j,bi,bj).eq.1.0) THEN

              sx = (surf_el_streamice(i,j,bi,bj)-
     &             surf_el_streamice(i-1,j,bi,bj))/dxC(i,j,bi,bj)

             ELSE

              sx = 0. _d 0

             ENDIF

             DO k=0,1
              DO l=0,1
               IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
                taudx_SI(i+k,j+l,bi,bj) = taudx_SI(i+k,j+l,bi,bj) -
     &           0.25 * streamice_density * gravity * 
     &           (streamice_bg_surf_slope_x+sx) * 
     &           H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
               ENDIF
              ENDDO
             ENDDO

            ENDIF

           ENDIF

!!!!!!!! DONE WITH X-GRADIENT

           IF (Gj.eq.1.AND..NOT.STREAMICE_NS_periodic) THEN

            ! western boundary - only one sided possible

            IF (STREAMICE_hmask(i,j+1,bi,bj).eq.1.0) THEN

             ! cell to east is active

             sy = (surf_el_streamice(i,j+1,bi,bj)-
     &             surf_el_streamice(i,j,bi,bj))/dyC(i,j+1,bi,bj)
            ELSE

             ! cell to east is empty

             sy = 0. _d 0
            ENDIF

            DO k=0,1
             DO l=0,1
              IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
               taudy_SI(i+k,j+l,bi,bj) = taudy_SI(i+k,j+l,bi,bj) -
     &          0.25 * streamice_density * gravity * 
     &          (streamice_bg_surf_slope_y+sy) * 
     &          H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
              ENDIF
             ENDDO
            ENDDO

           ELSEIF (Gj.eq.Ny.AND..NOT.STREAMICE_NS_periodic) THEN

            ! eastern boundary - only one sided possible

            IF (STREAMICE_hmask(i,j-1,bi,bj).eq.1.0) THEN

             ! cell to west is active

             sy = (surf_el_streamice(i,j,bi,bj)-
     &             surf_el_streamice(i,j-1,bi,bj))/dyC(i,j,bi,bj)

            ELSE

             ! cell to west is inactive

             sy = 0. _d 0
            ENDIF

            DO k=0,1
             DO l=0,1
              IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
               taudy_SI(i+k,j+l,bi,bj) = taudy_SI(i+k,j+l,bi,bj) -
     &          0.25 * streamice_density * gravity * 
     &          (streamice_bg_surf_slope_y+sy) * 
     &          H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
              ENDIF
             ENDDO
            ENDDO

           ELSE

            ! interior (west-east) cell

            IF (STREAMICE_hmask(i,j+1,bi,bj).eq.1.0 .and.
     &          STREAMICE_hmask(i,j-1,bi,bj).eq.1.0) THEN

             l = 0
             DO k=0,1
              IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
               taudy_SI(i+k,j+l,bi,bj) = taudy_SI(i+k,j+l,bi,bj) -
     &          streamice_density * gravity * (1./6.) *
     &          (-2.*surf_el_streamice(i,j-1,bi,bj) + 
     &           surf_el_streamice(i,j,bi,bj) + 
     &           surf_el_streamice(i,j+1,bi,bj) +
     &           3.*streamice_bg_surf_slope_y * dyF(i,j,bi,bj)) *
     &          H_streamice(i,j,bi,bj) * .5 * dxF(i,j,bi,bj)
              ENDIF
             ENDDO

             l = 1
             DO k=0,1
              IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
               taudy_SI(i+k,j+l,bi,bj) = taudy_SI(i+k,j+l,bi,bj) -
     &          streamice_density * gravity * (1./6.) *
     &          (-surf_el_streamice(i,j-1,bi,bj) - 
     &           surf_el_streamice(i,j,bi,bj) + 
     &           2*surf_el_streamice(i,j+1,bi,bj) +
     &           3.*streamice_bg_surf_slope_y * dyF(i,j,bi,bj)) *
     &          H_streamice(i,j,bi,bj) * .5 * dxF(i,j,bi,bj)
              ENDIF
             ENDDO


            ELSE 

             IF (STREAMICE_hmask(i,j+1,bi,bj).eq.1.0) THEN

              sy = (surf_el_streamice(i,j+1,bi,bj)-
     &             surf_el_streamice(i,j,bi,bj))/dxC(i,j+1,bi,bj)

             ELSEIF (STREAMICE_hmask(i,j-1,bi,bj).eq.1.0) THEN

              sy = (surf_el_streamice(i,j,bi,bj)-
     &             surf_el_streamice(i,j-1,bi,bj))/dxC(i,j,bi,bj)

             ELSE

              sy = 0. _d 0

             ENDIF

             DO k=0,1
              DO l=0,1
               IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
                taudy_SI(i+k,j+l,bi,bj) = taudy_SI(i+k,j+l,bi,bj) -
     &           0.25 * streamice_density * gravity * 
     &           (streamice_bg_surf_slope_y+sy) * 
     &           H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
               ENDIF
              ENDDO
             ENDDO

            ENDIF

           ENDIF

!            DO k=0,1
!             DO l=0,1
!              IF (STREAMICE_umask(i+k,j+l,bi,bj).eq.1.0) THEN
!               taudx_SI(i+k,j+l,bi,bj) = taudx_SI(i+k,j+l,bi,bj) -
!      &         0.25 * streamice_density * gravity * 
!      &         (streamice_bg_surf_slope_x+sx) * 
!      &         H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
!               taudy_SI(i+k,j+l,bi,bj) = taudy_SI(i+k,j+l,bi,bj) -
!      &         0.25 * streamice_density * gravity * 
!      &         (streamice_bg_surf_slope_y+sy) * 
!      &         H_streamice(i,j,bi,bj) * rA(i,j,bi,bj)
!               
!              ENDIF
!             ENDDO
!            ENDDO

           IF (float_frac_streamice(i,j,bi,bj) .eq. 1.0) then
#ifdef USE_ALT_RLOW
            neu_val = .5 * gravity * 
     &       (streamice_density * H_streamice (i,j,bi,bj) ** 2 - 
     &        streamice_density_ocean_avg * R_low_si(i,j,bi,bj) ** 2)
#else
            neu_val = .5 * gravity * 
     &       (streamice_density * H_streamice (i,j,bi,bj) ** 2 - 
     &        streamice_density_ocean_avg * R_low(i,j,bi,bj) ** 2)
#endif
           ELSE
            neu_val = .5 * gravity * 
     &       (1-streamice_density/streamice_density_ocean_avg) * 
     &        streamice_density * H_streamice(i,j,bi,bj) ** 2 
           ENDIF

           IF ((STREAMICE_ufacemask(i,j,bi,bj) .eq. 2) 
     &      .OR. (STREAMICE_hmask(i-1,j,bi,bj) .eq. 0)
     &      .OR. (STREAMICE_hmask(i-1,j,bi,bj) .eq. 2) ) THEN ! left face of the cell is at a stress boundary
          ! the depth-integrated longitudinal stress is equal to the difference of depth-integrated pressure on either side of the face
          ! on the ice side, it is rho g h^2 / 2
          ! on the ocean side, it is rhow g (delta OD)^2 / 2
          ! OD can be zero under the ice; but it is ASSUMED on the ice-free side of the face, topography elevation is not above the base of the 
          !     ice in the current cell
             
             taudx_SI(i,j,bi,bj) = taudx_SI(i,j,bi,bj) - 
     &        .5 * dyG(i,j,bi,bj)*(neu_val+streamice_addl_backstress)
             taudx_SI(i,j+1,bi,bj) = taudx_SI(i,j+1,bi,bj) - 
     &        .5 * dyG(i,j,bi,bj)*(neu_val+streamice_addl_backstress)
           ENDIF

           IF ((STREAMICE_ufacemask(i+1,j,bi,bj) .eq. 2) 
     &      .OR. (STREAMICE_hmask(i+1,j,bi,bj) .eq. 0)
     &      .OR. (STREAMICE_hmask(i+1,j,bi,bj) .eq. 2) ) THEN 
             
             taudx_SI(i+1,j,bi,bj) = taudx_SI(i+1,j,bi,bj) + 
     &        .5 * dyG(i+1,j,bi,bj)*(neu_val+streamice_addl_backstress)  ! note negative sign is due to direction of normal vector
             taudx_SI(i+1,j+1,bi,bj) = taudx_SI(i+1,j+1,bi,bj) +
     &        .5 * dyG(i+1,j,bi,bj)*(neu_val+streamice_addl_backstress)
           ENDIF

           IF ((STREAMICE_vfacemask(i,j,bi,bj) .eq. 2) 
     &      .OR. (STREAMICE_hmask(i,j-1,bi,bj) .eq. 0)
     &      .OR. (STREAMICE_hmask(i,j-1,bi,bj) .eq. 2) ) THEN 
             
             taudy_SI(i,j,bi,bj) = taudy_SI(i,j,bi,bj) - 
     &        .5 * dxG(i,j,bi,bj)*(neu_val+streamice_addl_backstress)
             taudy_SI(i+1,j,bi,bj) = taudy_SI(i+1,j,bi,bj) - 
     &        .5 * dxG(i,j,bi,bj)*(neu_val+streamice_addl_backstress)
           ENDIF

           IF ((STREAMICE_vfacemask(i,j+1,bi,bj) .eq. 2) 
     &      .OR. (STREAMICE_hmask(i,j+1,bi,bj) .eq. 0)
     &      .OR. (STREAMICE_hmask(i,j+1,bi,bj) .eq. 2) ) THEN 
             
             taudy_SI(i,j+1,bi,bj) = taudy_SI(i,j+1,bi,bj) +
     &        .5 * dxG(i,j+1,bi,bj)*(neu_val+streamice_addl_backstress)
             taudy_SI(i+1,j+1,bi,bj) = taudy_SI(i+1,j+1,bi,bj) +
     &        .5 * dxG(i,j+1,bi,bj)*(neu_val+streamice_addl_backstress)
           ENDIF

          ENDIF
         ENDDO
        ENDDO
       ENDDO
      ENDDO



#endif
      RETURN
      END
       
      
