C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_driving_stress_fd.F,v 1.2 2016/10/24 14:14:19 dgoldberg Exp $
C $Name:  $

#include "STREAMICE_OPTIONS.h"

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

CBOP
      SUBROUTINE STREAMICE_DRIVING_STRESS_FD( 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 dsdx
     & (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL dsdy
     & (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL face_factor, grd_below_sl, i_rhow
      

      i_rhow = 1./streamice_density_ocean_avg

      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO j=1-OLy+1,sNy+OLy-1
         DO i=1-OLx+1,sNx+OLx-1
          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 j=1,sNy
         DO i=1,sNx

          IF (streamice_hmask(i,j,bi,bj) .eq. 1) THEN

!    ds/dx ------------------------------------------------

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

            dsdx(i,j,bi,bj) = recip_dxC(i,j,bi,bj) *
     &       (surf_el_streamice(i,j,bi,bj) - 
     &        surf_el_streamice(i-1,j,bi,bj))

           ELSEIF (streamice_hmask(i-1,j,bi,bj).eq.0.or.
     &             streamice_hmask(i-1,j,bi,bj).eq.2) THEN

            IF (streamice_hmask(i+1,j,bi,bj).eq.1) THEN
             dsdx(i,j,bi,bj) = recip_dxC(i+1,j,bi,bj) *
     &       (surf_el_streamice(i+1,j,bi,bj) - 
     &        surf_el_streamice(i,j,bi,bj))
            ELSE
             dsdx(i,j,bi,bj) = 0.0
            ENDIF

           ELSE ! streamice_hmask(i-1,j,bi,bj).eq. -1
            dsdx(i,j,bi,bj) = 0.0
           ENDIF

!    ds/dy ------------------------------------------------

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

            dsdy(i,j,bi,bj) = recip_dyC(i,j,bi,bj) *
     &       (surf_el_streamice(i,j,bi,bj) - 
     &        surf_el_streamice(i,j-1,bi,bj))

           ELSEIF (streamice_hmask(i,j-1,bi,bj).eq.0.or.
     &             streamice_hmask(i,j-1,bi,bj).eq.2) THEN

            IF (streamice_hmask(i,j+1,bi,bj).eq.1) THEN
             dsdy(i,j,bi,bj) = recip_dyC(i,j+1,bi,bj) *
     &       (surf_el_streamice(i,j+1,bi,bj) - 
     &        surf_el_streamice(i,j,bi,bj))
            ELSE
             dsdy(i,j,bi,bj) = 0.0
            ENDIF

           ELSE ! streamice_hmask(i-1,j,bi,bj).eq. -1
            dsdx(i,j,bi,bj) = 0.0
           ENDIF

!    end ------------------------------------------------

          ELSEIF(streamice_hmask(i,j,bi,bj).eq.0.or.
     &            streamice_hmask(i,j,bi,bj).eq.2) THEN

!    ds/dx ------------------------------------------------

           IF(streamice_hmask(i-1,j,bi,bj).eq.1.and.
     &         streamice_hmask(i-2,j,bi,bj).eq.1) THEN

             dsdx(i,j,bi,bj) = recip_dxC(i-1,j,bi,bj) *
     &       (surf_el_streamice(i-1,j,bi,bj) - 
     &        surf_el_streamice(i-2,j,bi,bj))

           ELSE
            dsdx(i,j,bi,bj) = 0.0

           ENDIF

!    ds/dy ------------------------------------------------

           IF(streamice_hmask(i,j-1,bi,bj).eq.1.and.
     &         streamice_hmask(i,j-2,bi,bj).eq.1) THEN

             dsdy(i,j,bi,bj) = recip_dyC(i,j-1,bi,bj) *
     &       (surf_el_streamice(i,j-1,bi,bj) - 
     &        surf_el_streamice(i,j-2,bi,bj))

           ELSE
            dsdy(i,j,bi,bj) = 0.0

           ENDIF

!    end ------------------------------------------------
         
          ENDIF
         ENDDO
        ENDDO
       ENDDO
      ENDDO

      _EXCH_XY_RL (dsdy, mythid)
      _EXCH_XY_RL (dsdx, mythid)

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

!           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
          
           IF (streamice_umask(i,j,bi,bj).eq.1.0) THEN

            taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) + 
     &       gravity * streamice_density * 
     &       H_streamice (i,j,bi,bj) * dsdx (i,j,bi,bj) * 
     &       (0.5 * dyG(i,j,bi,bj) + 0.25 * dyF(i,j,bi,bj)) * 
     &       0.5 * dxG(i,j,bi,bj)

            taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) + 
     &       gravity * streamice_density * 
     &       H_streamice (i,j,bi,bj) * dsdx (i+1,j,bi,bj) * 
     &       0.25 * dyF(i,j,bi,bj) * 
     &       0.5 * dxG(i,j,bi,bj)

           ENDIF

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

            taudx_si(i+1,j,bi,bj) = taudx_si(i+1,j,bi,bj) + 
     &       gravity * streamice_density * 
     &       H_streamice (i,j,bi,bj) * dsdx (i+1,j,bi,bj) * 
     &       (0.5 * dyG(i+1,j,bi,bj) + 0.25 * dyF(i,j,bi,bj)) * 
     &       0.5 * dxG(i,j,bi,bj)

            taudx_si(i+1,j,bi,bj) = taudx_si(i+1,j,bi,bj) + 
     &       gravity * streamice_density * 
     &       H_streamice (i,j,bi,bj) * dsdx (i,j,bi,bj) * 
     &       0.25 * dyF(i,j,bi,bj) * 
     &       0.5 * dxG(i,j,bi,bj)

           ENDIF

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

            taudx_si(i,j+1,bi,bj) = taudx_si(i,j+1,bi,bj) + 
     &       gravity * streamice_density * 
     &       H_streamice (i,j,bi,bj) * dsdx (i,j,bi,bj) * 
     &       (0.5 * dyG(i,j,bi,bj) + 0.25 * dyF(i,j,bi,bj)) * 
     &       0.5 * dxG(i,j+1,bi,bj)

            taudx_si(i,j+1,bi,bj) = taudx_si(i,j+1,bi,bj) + 
     &       gravity * streamice_density * 
     &       H_streamice (i,j,bi,bj) * dsdx (i+1,j,bi,bj) * 
     &       0.25 * dyF(i,j,bi,bj) * 
     &       0.5 * dxG(i,j+1,bi,bj)

           ENDIF

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

            taudx_si(i+1,j+1,bi,bj) = taudx_si(i+1,j+1,bi,bj) + 
     &       gravity * streamice_density * 
     &       H_streamice (i,j,bi,bj) * dsdx (i+1,j,bi,bj) * 
     &       (0.5 * dyG(i+1,j,bi,bj) + 0.25 * dyF(i,j,bi,bj)) * 
     &       0.5 * dxG(i,j+1,bi,bj)

            taudx_si(i+1,j+1,bi,bj) = taudx_si(i+1,j+1,bi,bj) + 
     &       gravity * streamice_density * 
     &       H_streamice (i,j,bi,bj) * dsdx (i,j,bi,bj) * 
     &       0.25 * dyF(i,j,bi,bj) * 
     &       0.5 * dxG(i,j+1,bi,bj)

           ENDIF

#ifdef USE_ALT_RLOW
           IF (R_low_si(i,j,bi,bj) .lt. 0. _d 0) then
#else           
           IF (R_low(i,j,bi,bj) .lt. 0. _d 0) then
#endif             
              grd_below_sl = 1. _d 0
           else
              grd_below_sl = 0. _d 0
           endif


! check face to right
           IF (streamice_hmask(i+1,j,bi,bj).eq.0.or.
     &         streamice_hmask(i+1,j,bi,bj).eq.2.or.
     &         streamice_ufacemask(i+1,j,bi,bj).eq.2) THEN

            IF (float_frac_streamice(i,j,bi,bj) .eq. 1.0) THEN
     
             face_factor = 
     &        0.25 * dyG(i+1,j,bi,bj) *
     &        gravity * 
     &        streamice_density * H_streamice(i,j,bi,bj)**2-
#ifdef USE_ALT_RLOW
     &         streamice_density_ocean_avg * grd_below_sl * 
     &         R_low_si(i,j,bi,bj)**2
#else           
     &         streamice_density_ocean_avg * grd_below_sl * 
     &         R_low(i,j,bi,bj)**2
#endif

            ELSE

             face_factor = 
     &        0.25 * dyG(i+1,j,bi,bj) *
     &        streamice_density * gravity *
     &        (1-streamice_density*i_rhow) *
     &         H_streamice(i,j,bi,bj)**2 

            ENDIF

            taudx_si(i+1,j,bi,bj) = taudx_si(i+1,j,bi,bj) 
     &       - face_factor

            taudx_si(i+1,j+1,bi,bj) = taudx_si(i+1,j+1,bi,bj) 
     &       - face_factor

           ENDIF

! check face to left
           IF (streamice_hmask(i-1,j,bi,bj).eq.0.or.
     &         streamice_hmask(i-1,j,bi,bj).eq.2.or.
     &         streamice_ufacemask(i,j,bi,bj).eq.2) THEN

            IF (float_frac_streamice(i,j,bi,bj) .eq. 1.0) THEN
     
             face_factor = 
     &        0.25 * dyG(i,j,bi,bj) *
     &        gravity * 
     &        streamice_density * H_streamice(i,j,bi,bj)**2-
#ifdef USE_ALT_RLOW
     &         streamice_density_ocean_avg * grd_below_sl * 
     &         R_low_si(i,j,bi,bj)**2
#else           
     &         streamice_density_ocean_avg * grd_below_sl * 
     &         R_low(i,j,bi,bj)**2
#endif

            ELSE

             face_factor = 
     &        0.25 * dyG(i,j,bi,bj) *
     &        streamice_density * gravity *
     &        (1-streamice_density*i_rhow) *
     &         H_streamice(i,j,bi,bj)**2 

            ENDIF

            taudx_si(i,j,bi,bj) = taudx_si(i,j,bi,bj) 
     &       + face_factor

            taudx_si(i,j+1,bi,bj) = taudx_si(i,j+1,bi,bj) 
     &       + face_factor

           ENDIF

! Y FORCES           

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

            taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) + 
     &       gravity * streamice_density * 
     &       H_streamice (i,j,bi,bj) * dsdy (i,j,bi,bj) * 
     &       (0.5 * dxG(i,j,bi,bj) + 0.25 * dxF(i,j,bi,bj)) * 
     &       0.5 * dyG(i,j,bi,bj)

            taudy_si(i,j,bi,bj) = taudy_si(i,j,bi,bj) + 
     &       gravity * streamice_density * 
     &       H_streamice (i,j,bi,bj) * dsdy (i,j+1,bi,bj) * 
     &       0.25 * dxF(i,j,bi,bj) * 
     &       0.5 * dyG(i,j,bi,bj)

           ENDIF

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

            taudy_si(i,j+1,bi,bj) = taudy_si(i,j+1,bi,bj) + 
     &       gravity * streamice_density * 
     &       H_streamice (i,j,bi,bj) * dsdy (i,j+1,bi,bj) * 
     &       (0.5 * dxG(i,j+1,bi,bj) + 0.25 * dxF(i,j,bi,bj)) * 
     &       0.5 * dyG(i,j,bi,bj)

            taudy_si(i,j+1,bi,bj) = taudy_si(i,j+1,bi,bj) + 
     &       gravity * streamice_density * 
     &       H_streamice (i,j,bi,bj) * dsdy (i,j,bi,bj) * 
     &       0.25 * dxF(i,j,bi,bj) * 
     &       0.5 * dyG(i,j,bi,bj)

           ENDIF

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

            taudy_si(i+1,j,bi,bj) = taudy_si(i+1,j,bi,bj) + 
     &       gravity * streamice_density * 
     &       H_streamice (i,j,bi,bj) * dsdy (i,j,bi,bj) * 
     &       (0.5 * dxG(i,j,bi,bj) + 0.25 * dxF(i,j,bi,bj)) * 
     &       0.5 * dyG(i+1,j,bi,bj)

            taudy_si(i+1,j,bi,bj) = taudy_si(i+1,j,bi,bj) + 
     &       gravity * streamice_density * 
     &       H_streamice (i,j,bi,bj) * dsdy (i,j+1,bi,bj) * 
     &       0.25 * dxF(i,j,bi,bj) * 
     &       0.5 * dyG(i+1,j,bi,bj)

           ENDIF

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

            taudy_si(i+1,j+1,bi,bj) = taudy_si(i+1,j+1,bi,bj) + 
     &       gravity * streamice_density * 
     &       H_streamice (i,j,bi,bj) * dsdy (i,j+1,bi,bj) * 
     &       (0.5 * dxG(i,j+1,bi,bj) + 0.25 * dxF(i,j,bi,bj)) * 
     &       0.5 * dyG(i,j+1,bi,bj)

            taudy_si(i+1,j+1,bi,bj) = taudy_si(i+1,j+1,bi,bj) + 
     &       gravity * streamice_density * 
     &       H_streamice (i,j,bi,bj) * dsdy (i,j,bi,bj) * 
     &       0.25 * dxF(i,j,bi,bj) * 
     &       0.5 * dyG(i,j+1,bi,bj)

           ENDIF

#ifdef USE_ALT_RLOW
           IF (R_low_si(i,j,bi,bj) .lt. 0. _d 0) then
#else           
           IF (R_low(i,j,bi,bj) .lt. 0. _d 0) then
#endif             
              grd_below_sl = 1. _d 0
           else
              grd_below_sl = 0. _d 0
           endif


! check face to right
           IF (streamice_hmask(i,j+1,bi,bj).eq.0.or.
     &         streamice_hmask(i,j+1,bi,bj).eq.2.or.
     &         streamice_vfacemask(i,j+1,bi,bj).eq.2) THEN

            IF (float_frac_streamice(i,j,bi,bj) .eq. 1.0) THEN
     
             face_factor = 
     &        0.25 * dxG(i,j+1,bi,bj) *
     &        gravity * 
     &        streamice_density * H_streamice(i,j,bi,bj)**2-
#ifdef USE_ALT_RLOW
     &         streamice_density_ocean_avg * grd_below_sl * 
     &         R_low_si(i,j,bi,bj)**2
#else           
     &         streamice_density_ocean_avg * grd_below_sl * 
     &         R_low(i,j,bi,bj)**2
#endif

            ELSE

             face_factor = 
     &        0.25 * dxG(i,j+1,bi,bj) *
     &        streamice_density * gravity *
     &        (1-streamice_density*i_rhow) *
     &         H_streamice(i,j,bi,bj)**2 

            ENDIF

            taudy_si(i,j+1,bi,bj) = taudy_si(i,j+1,bi,bj) 
     &       - face_factor

            taudy_si(i+1,j+1,bi,bj) = taudy_si(i+1,j+1,bi,bj) 
     &       - face_factor

           ENDIF

! check face to left
           IF (streamice_hmask(i,j-1,bi,bj).eq.0.or.
     &         streamice_hmask(i,j-1,bi,bj).eq.2.or.
     &         streamice_vfacemask(i,j,bi,bj).eq.2) THEN

            IF (float_frac_streamice(i,j,bi,bj) .eq. 1.0) THEN
     
             face_factor = 
     &        0.25 * dxG(i,j,bi,bj) *
     &        gravity * 
     &        streamice_density * H_streamice(i,j,bi,bj)**2-
#ifdef USE_ALT_RLOW
     &         streamice_density_ocean_avg * grd_below_sl * 
     &         R_low_si(i,j,bi,bj)**2
#else           
     &         streamice_density_ocean_avg * grd_below_sl * 
     &         R_low(i,j,bi,bj)**2
#endif

            ELSE

             face_factor = 
     &        0.25 * dxG(i,j,bi,bj) *
     &        streamice_density * gravity *
     &        (1-streamice_density*i_rhow) *
     &         H_streamice(i,j,bi,bj)**2 

            ENDIF

            taudy_si(i,j,bi,bj) = taudy_SI(i,j,bi,bj) 
     &       + face_factor

            taudy_si(i+1,j,bi,bj) = taudy_si(i+1,j,bi,bj) 
     &       + face_factor

           ENDIF


          ENDIF

         ENDDO
        ENDDO
       ENDDO
      ENDDO

      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO j=1-OLy+1,sNy+OLy-1
         DO i=1-OLx+1,sNx+OLx-1
          taudx_SI(i,j,bi,bj) = -1.*taudx_SI(i,j,bi,bj)
          taudy_SI(i,j,bi,bj) = -1.*taudy_SI(i,j,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
!       call write_fld_xy_rl ('driving_taux','',taudx_si,0,mythid)

#endif
      RETURN
      END