C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_calc_stressdiv.F,v 1.1 2017/05/23 16:22:45 mlosch Exp $ C $Name: $ #include "SEAICE_OPTIONS.h" #ifdef ALLOW_AUTODIFF # include "AUTODIFF_OPTIONS.h" #endif C-- File seaice_calc_stressdiv.F C-- Contents C-- o SEAICE_CALC_STRESSDIV C-- o SEAICE_CALC_STRESS CBOP C !ROUTINE: SEAICE_CALC_STRESSDIV C !INTERFACE: ========================================================== SUBROUTINE SEAICE_CALC_STRESSDIV( I e11, e22, e12, press, zeta, eta, etaZ, O stressDivergenceX, stressDivergenceY, I bi, bj, myTime, myIter, myThid ) C !DESCRIPTION: \bv C *===========================================================* C | SUBROUTINE SEAICE_CALC_STRESSDIV C | o compute ice internal stress divergence C | C | Martin Losch, May 2017, Martin.Losch@awi.de C *===========================================================* C \ev C !USES: =============================================================== IMPLICIT NONE #include "SIZE.h" #include "EEPARAMS.h" #include "GRID.h" C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === C myTime :: Simulation time C myIter :: Simulation timestep number C myThid :: my Thread Id. number C e11/e22/e12 :: strain rate tensor components C press :: maximal compressive strength C zeta :: bulk viscosity C eta :: shear viscosity C etaZ :: shear viscosity at vorticity points C stressDivergenceX/Y :: x/y-component of stress divergence _RL myTime INTEGER myIter INTEGER myThid INTEGER bi, bj _RL e11 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL e22 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL e12 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL press(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL zeta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL eta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL etaZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL sig11(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL sig22(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL sig12(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL stressDivergenceX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL stressDivergenceY(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) CEOP #ifdef SEAICE_CGRID C !LOCAL VARIABLES: ==================================================== C === Local variables === C i,j :: inner loop counters C INTEGER i, j CALL SEAICE_CALC_STRESS( I e11, e22, e12, press, zeta, eta, etaZ, O sig11, sig22, sig12, I bi, bj, myTime, myIter, myThid ) C C compute divergence of stress tensor C DO J=1,sNy DO I=1,sNx stressDivergenceX(i,j,bi,bj) = & ( sig11(i ,j ) * _dyF(i ,j,bi,bj) & - sig11(i-1,J ) * _dyF(i-1,J,bi,bj) & + sig12(i ,j+1) * _dxV(i,j+1,bi,bj) & - sig12(i ,j ) * _dxV(i,j ,bi,bj) & ) * recip_rAw(i,j,bi,bj) stressDivergenceY(i,j,bi,bj) = & ( sig22(i ,j ) * _dxF(i,j ,bi,bj) & - sig22(i ,j-1) * _dxF(i,j-1,bi,bj) & + sig12(i+1,J ) * _dyU(i+1,J,bi,bj) & - sig12(i ,j ) * _dyU(i ,j,bi,bj) & ) * recip_rAs(i,j,bi,bj) ENDDO ENDDO RETURN END
CBOP C !ROUTINE: SEAICE_CALC_STRESS C !INTERFACE: ========================================================== SUBROUTINE SEAICE_CALC_STRESS( I e11, e22, e12, press, zeta, eta, etaZ, O sig11, sig22, sig12, I bi, bj, myTime, myIter, myThid ) C !DESCRIPTION: \bv C *===========================================================* C | SUBROUTINE SEAICE_CALC_STRESS C | o compute ice internal stress according to rheology C | C | Martin Losch, May 2017, Martin.Losch@awi.de C *===========================================================* C \ev C !USES: =============================================================== IMPLICIT NONE #include "SIZE.h" #include "EEPARAMS.h" #include "GRID.h" C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === C myTime :: Simulation time C myIter :: Simulation timestep number C myThid :: my Thread Id. number C e11/e22/e12 :: strain rate tensor components C press :: maximal compressive strength C zeta :: bulk viscosity C eta :: shear viscosity C etaZ :: shear viscosity at vorticity points C sig11/sig22/sig12 :: stress tensor components _RL myTime INTEGER myIter INTEGER myThid INTEGER bi, bj _RL e11 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL e22 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL e12 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL press(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL zeta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL eta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL etaZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL sig11(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL sig22(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL sig12(1-OLx:sNx+OLx,1-OLy:sNy+OLy) CEOP C !LOCAL VARIABLES: ==================================================== C === Local variables === C i,j :: inner loop counters C eplus, eminus :: convenient abbreviations for e11+e22, e11-e22 INTEGER i, j _RL eplus, eminus C compute components of stress tensor from current strain rate tensor DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx sig11(i,j) = 0. _d 0 sig22(i,j) = 0. _d 0 sig12(i,j) = 0. _d 0 ENDDO ENDDO DO j=0,sNy DO i=0,sNx eplus = e11(i,j,bi,bj) + e22(i,j,bi,bj) eminus= e11(i,j,bi,bj) - e22(i,j,bi,bj) sig11(i,j) = zeta(i,j,bi,bj)*eplus + eta(i,j,bi,bj)*eminus & - 0.5 _d 0 * press(i,j,bi,bj) sig22(i,j) = zeta(i,j,bi,bj)*eplus - eta(i,j,bi,bj)*eminus & - 0.5 _d 0 * press(i,j,bi,bj) ENDDO ENDDO DO j=1,sNy+1 DO i=1,sNx+1 sig12(i,j) = 2. _d 0 * e12(i,j,bi,bj) * etaZ(i,j,bi,bj) ENDDO ENDDO #endif /* SEAICE_CGRID */ RETURN END