Back to home page

MITgcm

 
 

    


File indexing completed on 2026-05-05 05:09:01 UTC

view on githubraw file Latest commit 3f0f10fc on 2026-05-04 14:55:37 UTC
925df4f4b9 Mart*0001 #include "SEAICE_OPTIONS.h"
                0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
                0005 
                0006 C--  File seaice_calc_stressdiv.F
                0007 C--   Contents
                0008 C--   o SEAICE_CALC_STRESSDIV
                0009 C--   o SEAICE_CALC_STRESS
                0010 
                0011 CBOP
                0012 C !ROUTINE: SEAICE_CALC_STRESSDIV
                0013 C !INTERFACE: ==========================================================
cbf501ab81 Jean*0014       SUBROUTINE SEAICE_CALC_STRESSDIV(
925df4f4b9 Mart*0015      I     e11, e22, e12, press, zeta, eta, etaZ,
                0016      O     stressDivergenceX, stressDivergenceY,
                0017      I     bi, bj, myTime, myIter, myThid )
                0018 
                0019 C !DESCRIPTION: \bv
                0020 C     *===========================================================*
                0021 C     | SUBROUTINE SEAICE_CALC_STRESSDIV
                0022 C     | o compute ice internal stress divergence
                0023 C     |
                0024 C     | Martin Losch, May 2017, Martin.Losch@awi.de
                0025 C     *===========================================================*
                0026 C \ev
                0027 
                0028 C !USES: ===============================================================
                0029       IMPLICIT NONE
                0030 
                0031 #include "SIZE.h"
                0032 #include "EEPARAMS.h"
                0033 #include "GRID.h"
3f0f10fc37 Mart*0034 #include "SEAICE_GRID.h"
925df4f4b9 Mart*0035 
                0036 C     !INPUT/OUTPUT PARAMETERS:
                0037 C     === Routine arguments ===
                0038 C     e11/e22/e12 :: strain rate tensor components
                0039 C     press  :: maximal compressive strength
                0040 C     zeta   :: bulk viscosity
cbf501ab81 Jean*0041 C     eta    :: shear viscosity
925df4f4b9 Mart*0042 C     etaZ   :: shear viscosity at vorticity points
                0043 C     stressDivergenceX/Y :: x/y-component of stress divergence
cbf501ab81 Jean*0044 C     myTime :: Simulation time
                0045 C     myIter :: Simulation timestep number
                0046 C     myThid :: my Thread Id. number
925df4f4b9 Mart*0047       _RL e11  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0048       _RL e22  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0049       _RL e12  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0050       _RL press(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0051       _RL zeta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0052       _RL eta  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0053       _RL etaZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0054       _RL stressDivergenceX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0055       _RL stressDivergenceY(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
cbf501ab81 Jean*0056       INTEGER bi, bj
                0057       _RL     myTime
                0058       INTEGER myIter
                0059       INTEGER myThid
925df4f4b9 Mart*0060 CEOP
                0061 
                0062 #ifdef SEAICE_CGRID
                0063 C !LOCAL VARIABLES: ====================================================
                0064 C     === Local variables ===
                0065 C     i,j       :: inner loop counters
cbf501ab81 Jean*0066       _RL sig11(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0067       _RL sig22(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0068       _RL sig12(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
925df4f4b9 Mart*0069       INTEGER i, j
                0070 
cbf501ab81 Jean*0071       CALL SEAICE_CALC_STRESS(
925df4f4b9 Mart*0072      I     e11, e22, e12, press, zeta, eta, etaZ,
                0073      O     sig11, sig22, sig12,
                0074      I     bi, bj, myTime, myIter, myThid )
                0075 
                0076 C     compute divergence of stress tensor
cbf501ab81 Jean*0077 
5fe78992ba Mart*0078       DO j=1,sNy+1
                0079        DO i=1,sNx+1
925df4f4b9 Mart*0080         stressDivergenceX(i,j,bi,bj) =
                0081      &       ( sig11(i  ,j  ) * _dyF(i  ,j,bi,bj)
5fe78992ba Mart*0082      &       - sig11(i-1,j  ) * _dyF(i-1,j,bi,bj)
925df4f4b9 Mart*0083      &       + sig12(i  ,j+1) * _dxV(i,j+1,bi,bj)
                0084      &       - sig12(i  ,j  ) * _dxV(i,j  ,bi,bj)
                0085      &       ) * recip_rAw(i,j,bi,bj)
3f0f10fc37 Mart*0086 C     extra metric terms
                0087      &       + k2AtU(i,j,bi,bj) * halfRL*( sig12(i,j) + sig12(i,j+1) )
                0088      &       - k1AtU(i,j,bi,bj) * halfRL*( sig22(i,j) + sig22(i-1,j) )
925df4f4b9 Mart*0089         stressDivergenceY(i,j,bi,bj) =
                0090      &       ( sig22(i  ,j  ) * _dxF(i,j  ,bi,bj)
                0091      &       - sig22(i  ,j-1) * _dxF(i,j-1,bi,bj)
5fe78992ba Mart*0092      &       + sig12(i+1,j  ) * _dyU(i+1,j,bi,bj)
925df4f4b9 Mart*0093      &       - sig12(i  ,j  ) * _dyU(i  ,j,bi,bj)
                0094      &       ) * recip_rAs(i,j,bi,bj)
3f0f10fc37 Mart*0095 C     extra metric terms
                0096      &       + k1AtV(i,j,bi,bj) * halfRL*( sig12(i,j) + sig12(i+1,j) )
                0097      &       - k2AtV(i,j,bi,bj) * halfRL*( sig11(i,j) + sig11(i,j-1) )
925df4f4b9 Mart*0098        ENDDO
                0099       ENDDO
                0100 
                0101       RETURN
                0102       END
                0103 
                0104 CBOP
                0105 C !ROUTINE: SEAICE_CALC_STRESS
                0106 C !INTERFACE: ==========================================================
cbf501ab81 Jean*0107       SUBROUTINE SEAICE_CALC_STRESS(
925df4f4b9 Mart*0108      I     e11, e22, e12, press, zeta, eta, etaZ,
                0109      O     sig11, sig22, sig12,
                0110      I     bi, bj, myTime, myIter, myThid )
                0111 
                0112 C !DESCRIPTION: \bv
                0113 C     *===========================================================*
                0114 C     | SUBROUTINE SEAICE_CALC_STRESS
                0115 C     | o compute ice internal stress according to rheology
                0116 C     |
                0117 C     | Martin Losch, May 2017, Martin.Losch@awi.de
                0118 C     *===========================================================*
                0119 C \ev
                0120 
                0121 C !USES: ===============================================================
                0122       IMPLICIT NONE
                0123 
                0124 #include "SIZE.h"
                0125 #include "EEPARAMS.h"
                0126 #include "GRID.h"
                0127 
                0128 C     !INPUT/OUTPUT PARAMETERS:
                0129 C     === Routine arguments ===
                0130 C     myTime :: Simulation time
                0131 C     myIter :: Simulation timestep number
                0132 C     myThid :: my Thread Id. number
                0133 C     e11/e22/e12 :: strain rate tensor components
                0134 C     press  :: maximal compressive strength
                0135 C     zeta   :: bulk viscosity
cbf501ab81 Jean*0136 C     eta    :: shear viscosity
925df4f4b9 Mart*0137 C     etaZ   :: shear viscosity at vorticity points
                0138 C     sig11/sig22/sig12   :: stress tensor components
                0139       _RL     myTime
                0140       INTEGER myIter
                0141       INTEGER myThid
                0142       INTEGER bi, bj
                0143       _RL e11  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0144       _RL e22  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0145       _RL e12  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0146       _RL press(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0147       _RL zeta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0148       _RL eta  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0149       _RL etaZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0150       _RL sig11(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0151       _RL sig22(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0152       _RL sig12(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0153 CEOP
                0154 
                0155 C !LOCAL VARIABLES: ====================================================
                0156 C     === Local variables ===
                0157 C     i,j           :: inner loop counters
                0158 C     eplus, eminus :: convenient abbreviations for e11+e22, e11-e22
                0159       INTEGER i, j
                0160       _RL     eplus, eminus
                0161 
                0162 C     compute components of stress tensor from current strain rate tensor
                0163       DO j=1-OLy,sNy+OLy
                0164        DO i=1-OLx,sNx+OLx
                0165         sig11(i,j) = 0. _d 0
                0166         sig22(i,j) = 0. _d 0
                0167         sig12(i,j) = 0. _d 0
                0168        ENDDO
                0169       ENDDO
                0170 
5fe78992ba Mart*0171       DO j=0,sNy+1
                0172        DO i=0,sNx+1
925df4f4b9 Mart*0173         eplus = e11(i,j,bi,bj) + e22(i,j,bi,bj)
                0174         eminus= e11(i,j,bi,bj) - e22(i,j,bi,bj)
                0175         sig11(i,j) = zeta(i,j,bi,bj)*eplus + eta(i,j,bi,bj)*eminus
                0176      &       - 0.5 _d 0 * press(i,j,bi,bj)
                0177         sig22(i,j) = zeta(i,j,bi,bj)*eplus - eta(i,j,bi,bj)*eminus
                0178      &       - 0.5 _d 0 * press(i,j,bi,bj)
                0179        ENDDO
                0180       ENDDO
cbf501ab81 Jean*0181 
5fe78992ba Mart*0182       DO j=1,sNy+2
                0183        DO i=1,sNx+2
925df4f4b9 Mart*0184         sig12(i,j) = 2. _d 0 * e12(i,j,bi,bj) * etaZ(i,j,bi,bj)
                0185        ENDDO
                0186       ENDDO
                0187 #endif /* SEAICE_CGRID */
                0188       RETURN
                0189       END