Back to home page

MITgcm

 
 

    


File indexing completed on 2024-06-06 05:11:10 UTC

view on githubraw file Latest commit af61e5eb on 2024-06-06 03:30:35 UTC
6e2f4e58fa Mart*0001 #include "SEAICE_OPTIONS.h"
fec75090d3 Jean*0002 #ifdef ALLOW_OBCS
                0003 # include "OBCS_OPTIONS.h"
                0004 #else
                0005 # define OBCS_UVICE_OLD
                0006 #endif
772b2ed80e Gael*0007 #ifdef ALLOW_AUTODIFF
                0008 # include "AUTODIFF_OPTIONS.h"
                0009 #endif
6e2f4e58fa Mart*0010 
c9a6744012 Jean*0011 CBOP
                0012 C     !ROUTINE: SEAICE_EVP
                0013 C     !INTERFACE:
6e2f4e58fa Mart*0014       SUBROUTINE SEAICE_EVP( myTime, myIter, myThid )
c9a6744012 Jean*0015 
                0016 C     !DESCRIPTION: \bv
                0017 C     *==========================================================*
                0018 C     | SUBROUTINE SEAICE_EVP
                0019 C     | o Ice dynamics using an EVP solver following
                0020 C     |   E. C. Hunke and J. K. Dukowicz. An
                0021 C     |   Elastic-Viscous-Plastic Model for Sea Ice Dynamics,
                0022 C     |   J. Phys. Oceanogr., 27, 1849-1867 (1997).
                0023 C     *==========================================================*
                0024 C     | written by Martin Losch, March 2006
                0025 C     *==========================================================*
                0026 C     \ev
                0027 
                0028 C     !USES:
6e2f4e58fa Mart*0029       IMPLICIT NONE
                0030 
                0031 C     === Global variables ===
                0032 #include "SIZE.h"
                0033 #include "EEPARAMS.h"
                0034 #include "PARAMS.h"
bcf4255045 Dimi*0035 #include "DYNVARS.h"
6e2f4e58fa Mart*0036 #include "GRID.h"
ccaa3c61f4 Patr*0037 #include "SEAICE_SIZE.h"
6e2f4e58fa Mart*0038 #include "SEAICE_PARAMS.h"
ccaa3c61f4 Patr*0039 #include "SEAICE.h"
6e2f4e58fa Mart*0040 
                0041 #ifdef ALLOW_AUTODIFF_TAMC
                0042 # include "tamc.h"
                0043 #endif
                0044 
c9a6744012 Jean*0045 C     !INPUT/OUTPUT PARAMETERS:
6e2f4e58fa Mart*0046 C     === Routine arguments ===
c9a6744012 Jean*0047 C     myTime :: Simulation time
                0048 C     myIter :: Simulation timestep number
                0049 C     myThid :: my Thread Id. number
6e2f4e58fa Mart*0050       _RL     myTime
                0051       INTEGER myIter
                0052       INTEGER myThid
                0053 
45315406aa Mart*0054 #if ( defined SEAICE_CGRID &&  defined SEAICE_ALLOW_EVP )
6e2f4e58fa Mart*0055 
                0056 C     === Local variables ===
c9a6744012 Jean*0057 C     i,j,bi,bj      :: Loop counters
                0058 C     kSrf           :: vertical index of surface layer
                0059 C     nEVPstep       :: number of timesteps within the EVP solver
                0060 C     iEVPstep       :: Loop counter
                0061 C     SIN/COSWAT     :: sine/cosine of turning angle
a9588c2504 Mart*0062 C     (recip_)ecc2   :: (one over) eccentricity squared
                0063 C     recip_evpAlpha :: 1/SEAICE_evpAlpha
                0064 C     recip_deltaT   :: 1/SEAICE_deltaTdyn
                0065 C     evpStarFac     :: 1 if SEAICEuseEVPstar = .true., 0 otherwise
                0066 C     betaFac        :: SEAICE_evpBeta/SEAICE_deltaTdyn=1/SEAICE_deltaTevp
                0067 C     betaFacP1      :: betaFac + evpStarFac/SEAICE_deltaTdyn
c9a6744012 Jean*0068 C     e11,e12,e22    :: components of strain rate tensor
                0069 C     seaice_div     :: divergence strain rates at C-points times P
                0070 C                       /divided by Delta minus 1
                0071 C     seaice_tension :: tension    strain rates at C-points times P
                0072 C                       /divided by Delta
                0073 C     seaice_shear   :: shear      strain rates, defined at Z-points times P
                0074 C                       /divided by Delta
                0075 C     sig11, sig22   :: sum and difference of diagonal terms of stress tensor
e4120ef5ff Jean*0076 C     modification for adaptive alpha and beta
85f5225996 Mart*0077 C               (see Kimmritz, Danilov, Losch 2015 for gamma << alpha beta)
a6389affa6 Mart*0078 C     EVPcFac        :: SEAICE_deltaTdyn*SEAICEaEVPcStar*(SEAICEaEVPcoeff*PI)**2
                0079 C                        with
                0080 C     SEAICEaEVPcStar:: multiple of stabilty factor: alpha*beta = cstar*gamma
                0081 C     SEAICEaEVPcoeff:: largest stabilized frequency according to
85f5225996 Mart*0082 C                        gamma = zeta * (cfac/cellarea)*deltaT/m
                0083 C                                with   (cfac/cellarea) <= pi**2/cellarea
e4120ef5ff Jean*0084 C     evpAlphaC/Z    :: alpha field on C points and on Z points
85f5225996 Mart*0085 C                        := sqrt(cstar gamma)
e4120ef5ff Jean*0086 C     evpBetaU/V     :: beta field on u and on v points
85f5225996 Mart*0087 C                        := sqrt(cstar gamma)
e4120ef5ff Jean*0088 C     evpAlphaMin    :: lower limit of alpha and beta, regularisation
                0089 C                     to prevent singularities of system matrix,
85f5225996 Mart*0090 C                     e.g. when ice concentration is too low.
e4120ef5ff Jean*0091 C     betaFacP1U/V   :: = betaFacP1 in standard case,
                0092 C                          with varying beta in the adaptive case
85f5225996 Mart*0093 C                           on u and on v point
                0094 C     betaFacU/V     :: analog betaFacP1U/V
b4949dd6db Jean*0095 
6e2f4e58fa Mart*0096       INTEGER i, j, bi, bj
282fbc7b4c Mart*0097       INTEGER kSrf
6e2f4e58fa Mart*0098       INTEGER nEVPstep, iEVPstep
8377b8ee87 Mart*0099 #ifndef ALLOW_AUTODIFF
                0100       INTEGER nEVPstepMax
                0101 #endif
452bde753c Mart*0102 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0103 C     evpkey :: tape key (depends on evp iteration)
                0104 C     tkey   :: tape key (depends on tile indices and evpkey)
                0105       INTEGER evpkey, tkey
cf9fa44a59 Patr*0106 #endif
6e2f4e58fa Mart*0107 
35fda33b05 Jean*0108       _RL COSWAT
                0109       _RS SINWAT
a9588c2504 Mart*0110       _RL ecc2, recip_ecc2, recip_evpAlpha, recip_deltaT
0912318009 Mart*0111       _RL betaFacP1, betaFac, evpStarFac, evpRevFac, recip_evpRevFac
6e2f4e58fa Mart*0112 
b7030e3a29 Mart*0113       _RL seaice_div    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0114       _RL seaice_tension(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0115       _RL seaice_shear  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68ff576152 Mart*0116       _RL sig11         (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0117       _RL sig22         (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70e078b38a Mart*0118 C     fractional area at velocity points
4d65bb090e Mart*0119       _RL areaW         (1:sNx,1:sNy,nSx,nSy)
                0120       _RL areaS         (1:sNx,1:sNy,nSx,nSy)
037351a1a6 Mart*0121 C     auxilliary variables
da4a30d88f Mart*0122       _RL ep            (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0123       _RL em            (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
d2082aecba Mart*0124       _RL e12Csq        (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1dd53d2d6a Mart*0125       _RL pressC        (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0126       _RL zetaC         (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0127       _RL deltaZ        (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
0fd6b2de1a Mart*0128 C     _RL zetaZ         (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
8ce59fd29e Mart*0129 #ifdef SEAICE_ALLOW_MOM_ADVECTION
                0130 C     tendency due to advection of momentum
                0131       _RL gUmom         (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0132       _RL gVmom         (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0133 #endif /*  SEAICE_ALLOW_MOM_ADVECTION */
e4120ef5ff Jean*0134       _RL deltaCreg, deltaSq, deltaMinSq, tmp
452bde753c Mart*0135 #ifdef SEAICE_ALLOW_TEM
b8e39f4242 Mart*0136       _RL etaDenC, zetaMaxC, etaDenZ, zetaMaxZ
452bde753c Mart*0137 #endif /* SEAICE_ALLOW_TEM */
fc27138d73 Mart*0138 #ifdef SEAICE_ALLOW_CLIPZETA
                0139       _RL zMaxZ, zMinZ, fac
                0140 #endif /* SEAICE_ALLOW_CLIPZETA */
e4120ef5ff Jean*0141       _RL denom1   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
a6389affa6 Mart*0142       _RL denom2   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
967b0b9f87 Mart*0143       _RL sumNorm, denomU, denomV
fec75090d3 Jean*0144       _RL locMaskU, locMaskV
a6389affa6 Mart*0145       _RL EVPcFac
85f5225996 Mart*0146       _RL evpAlphaC(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
a6389affa6 Mart*0147       _RL evpAlphaZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0148       _RL evpBetaU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0149       _RL evpBetaV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
85f5225996 Mart*0150       _RL betaFacP1U,  betaFacP1V
                0151       _RL betaFacU,    betaFacV
a6389affa6 Mart*0152       LOGICAL useAdaptiveEVP
85f5225996 Mart*0153 
e4120ef5ff Jean*0154 #ifdef ALLOW_SEAICE_EVP_RESIDUAL
85f5225996 Mart*0155       _RL resTile(nSx,nSy)
                0156       _RL resLoc
4669aaa4e4 Mart*0157       _RL uIcePm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0158       _RL vIcePm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0159       _RL sig11pm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0160       _RL sig22pm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0161       _RL sig12pm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
85f5225996 Mart*0162 
                0163       LOGICAL printResidual
                0164 C     CHARACTER*(10) suff
                0165 C     !FUNCTIONS:
                0166       LOGICAL  DIFFERENT_MULTIPLE
                0167       EXTERNAL DIFFERENT_MULTIPLE
                0168 #endif /* ALLOW_SEAICE_EVP_RESIDUAL */
                0169 CEOP
6e2f4e58fa Mart*0170 
a6389affa6 Mart*0171 C     set tuning parameters for adaptive EVP
                0172       useAdaptiveEVP = .FALSE.
                0173       IF ( SEAICEaEvpCoeff .NE. UNSET_RL ) useAdaptiveEVP = .TRUE.
                0174       EVPcFac = 0. _d 0
                0175       IF ( useAdaptiveEVP )
                0176      &     EVPcFac = SEAICE_deltaTdyn*SEAICEaEVPcStar
                0177      &                * (SEAICEaEvpCoeff * PI)**2
85f5225996 Mart*0178 
e4120ef5ff Jean*0179 #ifdef ALLOW_SEAICE_EVP_RESIDUAL
85f5225996 Mart*0180       printResidual = debugLevel.GE.debLevA
                0181      &  .AND. DIFFERENT_MULTIPLE( SEAICE_monFreq, myTime, deltaTClock )
                0182 #endif /* ALLOW_SEAICE_EVP_RESIDUAL */
c9a6744012 Jean*0183 C     surface level
0320e25227 Mart*0184       IF ( usingPCoords ) THEN
                0185        kSrf = Nr
                0186       ELSE
                0187        kSrf = 1
                0188       ENDIF
6e2f4e58fa Mart*0189 C--   introduce turning angles
                0190       SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)
                0191       COSWAT=COS(SEAICE_waterTurnAngle*deg2rad)
                0192 
                0193 C     abbreviation eccentricity squared
                0194       ecc2=SEAICE_eccen**2
037351a1a6 Mart*0195       recip_ecc2 = 0. _d 0
                0196       IF ( ecc2 .NE. 0. _d 0 ) recip_ecc2=ONE/ecc2
3daf25222c Mart*0197       deltaMinSq = SEAICE_deltaMin**2
dc26f158aa Mart*0198 C     copy number of internal time steps from previously defined parameter
                0199       nEVPstep = SEAICEnEVPstarSteps
fb2ddc4f23 Mart*0200 C     SEAICE_evpAlpha = 2. * SEAICE_evpTauRelax/SEAICE_deltaTevp
a6389affa6 Mart*0201       DO bj=myByLo(myThid),myByHi(myThid)
                0202        DO bi=myBxLo(myThid),myBxHi(myThid)
8377b8ee87 Mart*0203         DO j=1-OLy,sNy+OLy
                0204          DO i=1-OLx,sNx+OLx
                0205           denom1(i,j,bi,bj) = 1. _d 0 / ( SEAICE_evpAlpha + 1. _d 0 )
                0206           denom2(i,j,bi,bj) = 1. _d 0 / ( SEAICE_evpAlpha + ecc2 )
a6389affa6 Mart*0207          ENDDO
                0208         ENDDO
                0209        ENDDO
                0210       ENDDO
2bf4ed5c39 Mart*0211       recip_deltaT = 1. _d 0 / SEAICE_deltaTdyn
                0212       recip_evpAlpha = 0. _d 0
e4120ef5ff Jean*0213       IF ( SEAICE_evpAlpha .GT. 0. _d 0 )
2bf4ed5c39 Mart*0214      &     recip_evpAlpha = 1. _d 0 / SEAICE_evpAlpha
                0215       evpStarFac = 0. _d 0
0912318009 Mart*0216       evpRevFac  = 0. _d 0
                0217       recip_evpRevFac = 1. _d 0
2bf4ed5c39 Mart*0218       IF ( SEAICEuseEVPstar ) evpStarFac = 1. _d 0
0912318009 Mart*0219       IF ( SEAICEuseEVPrev  ) THEN
e4120ef5ff Jean*0220 C     the Bouillon et al. (2013) discretization in time has  more
d2082aecba Mart*0221 C     explicit terms
0912318009 Mart*0222        evpRevFac       = 1. _d 0
                0223        evpStarFac      = 1. _d 0
                0224        recip_evpRevFac = recip_ecc2
a6389affa6 Mart*0225        DO bj=myByLo(myThid),myByHi(myThid)
                0226         DO bi=myBxLo(myThid),myBxHi(myThid)
8377b8ee87 Mart*0227          DO j=1-OLy,sNy+OLy
                0228           DO i=1-OLx,sNx+OLx
                0229            denom1(i,j,bi,bj)     = 1. _d 0 / SEAICE_evpAlpha
                0230            denom2(i,j,bi,bj)     = denom1(i,j,bi,bj)
a6389affa6 Mart*0231           ENDDO
                0232          ENDDO
                0233         ENDDO
                0234        ENDDO
0912318009 Mart*0235       ENDIF
85f5225996 Mart*0236 C
                0237       betaFac    = SEAICE_evpBeta*recip_deltaT
                0238       betaFacU   = betaFac
                0239       betaFacV   = betaFac
                0240 C
                0241       betaFacP1  = betaFac + evpStarFac*recip_deltaT
                0242       betaFacP1U = betaFacP1
                0243       betaFacP1V = betaFacP1
8377b8ee87 Mart*0244 #ifndef ALLOW_AUTODIFF
cf9fa44a59 Patr*0245       nEVPstepMax = nEVPstep
                0246 #endif
6e2f4e58fa Mart*0247 
                0248       DO bj=myByLo(myThid),myByHi(myThid)
                0249        DO bi=myBxLo(myThid),myBxHi(myThid)
9e706e0219 Jean*0250         DO j=1-OLy,sNy+OLy
                0251          DO i=1-OLx,sNx+OLx
a3f6dcaab3 Mart*0252 C     use u/vIce as work arrays: copy previous time step to u/vIceNm1
8377b8ee87 Mart*0253           uIceNm1(i,j,bi,bj)   = uIce(i,j,bi,bj)
                0254           vIceNm1(i,j,bi,bj)   = vIce(i,j,bi,bj)
c9a6744012 Jean*0255 C     initialise strain rates
8377b8ee87 Mart*0256           e11  (i,j,bi,bj)     = 0. _d 0
                0257           e22  (i,j,bi,bj)     = 0. _d 0
                0258           e12  (i,j,bi,bj)     = 0. _d 0
a6389affa6 Mart*0259 C     initialise adative-EVP-specific fields
8377b8ee87 Mart*0260           evpAlphaC(i,j,bi,bj) = SEAICE_evpAlpha
                0261           evpAlphaZ(i,j,bi,bj) = SEAICE_evpAlpha
                0262           evpBetaU (i,j,bi,bj) = SEAICE_evpBeta
                0263           evpBetaV (i,j,bi,bj) = SEAICE_evpBeta
6e2f4e58fa Mart*0264          ENDDO
                0265         ENDDO
4d65bb090e Mart*0266 C     initialise fractional areas at velocity points
                0267         IF ( SEAICEscaleSurfStress ) THEN
8377b8ee87 Mart*0268          DO j=1,sNy
                0269           DO i=1,sNx
                0270            areaW(i,j,bi,bj) =
                0271      &          0.5 _d 0*(AREA(i,j,bi,bj)+AREA(i-1,j,bi,bj))
                0272            areaS(i,j,bi,bj) =
                0273      &          0.5 _d 0*(AREA(i,j,bi,bj)+AREA(i,j-1,bi,bj))
4d65bb090e Mart*0274           ENDDO
                0275          ENDDO
                0276         ELSE
8377b8ee87 Mart*0277          DO j=1,sNy
                0278           DO i=1,sNx
                0279            areaW(i,j,bi,bj) = 1. _d 0
                0280            areaS(i,j,bi,bj) = 1. _d 0
4d65bb090e Mart*0281           ENDDO
                0282          ENDDO
                0283         ENDIF
cf9fa44a59 Patr*0284        ENDDO
                0285       ENDDO
fc27138d73 Mart*0286 #ifdef SEAICE_ALLOW_CLIPZETA
8a89485793 Mart*0287 C     damping constraint (Hunke, J.Comp.Phys.,2001)
                0288       IF ( SEAICE_evpDampC .GT. 0. _d 0 ) THEN
2bf4ed5c39 Mart*0289 CML       fac = HALF * SEAICE_evpDampC * SEAICE_evpTauRelax
                0290 CML     &      /SEAICE_deltaTevp**2
                0291        fac = 0.25 _d 0 * SEAICE_evpDampC * SEAICE_evpAlpha
                0292      &      /SEAICE_deltaTevp
8a89485793 Mart*0293        DO bj=myByLo(myThid),myByHi(myThid)
                0294         DO bi=myBxLo(myThid),myBxHi(myThid)
9e706e0219 Jean*0295          DO j=1-OLy,sNy+OLy
                0296           DO i=1-OLx,sNx+OLx
8e32c48b8f Mart*0297            SEAICE_zMax(i,j,bi,bj) = _rA(i,j,bi,bj) * fac
8a89485793 Mart*0298           ENDDO
                0299          ENDDO
                0300         ENDDO
                0301        ENDDO
                0302       ENDIF
fc27138d73 Mart*0303 #endif /* SEAICE_ALLOW_CLIPZETA */
b4949dd6db Jean*0304 C
6e2f4e58fa Mart*0305 C     start of the main time loop
cf9fa44a59 Patr*0306       DO iEVPstep = 1, nEVPstepMax
                0307        IF (iEVPstep.LE.nEVPstep) THEN
f6a7d7e353 Patr*0308 C
6e2f4e58fa Mart*0309 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0310         evpkey = iEVPstep + (ikey_dynamics-1)*nEVPstepMax
                0311 CADJ STORE uice           = comlev1_evp, key = evpkey, byte = isbyte
                0312 CADJ STORE vice           = comlev1_evp, key = evpkey, byte = isbyte
                0313 CADJ STORE seaice_sigma1  = comlev1_evp, key = evpkey, byte = isbyte
                0314 CADJ STORE seaice_sigma2  = comlev1_evp, key = evpkey, byte = isbyte
                0315 CADJ STORE seaice_sigma12 = comlev1_evp, key = evpkey, byte = isbyte
6e2f4e58fa Mart*0316 #endif /* ALLOW_AUTODIFF_TAMC */
                0317 C
                0318 C     first calculate strain rates and bulk moduli/viscosities
b4949dd6db Jean*0319 C
                0320         CALL SEAICE_CALC_STRAINRATES(
a3f6dcaab3 Mart*0321      I       uIce, vIce,
037351a1a6 Mart*0322      O       e11, e22, e12,
2e75568507 Mart*0323      I       iEVPstep, myTime, myIter, myThid )
037351a1a6 Mart*0324 
88601fb2cb Gael*0325 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0326 CADJ STORE e22 = comlev1_evp,key = evpkey, byte = isbyte
88601fb2cb Gael*0327 #endif /* ALLOW_AUTODIFF_TAMC */
                0328 
037351a1a6 Mart*0329         DO bj=myByLo(myThid),myByHi(myThid)
                0330          DO bi=myBxLo(myThid),myBxHi(myThid)
77fc19a7e9 Mart*0331 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0332           tkey = bi + (bj-1)*nSx + (evpkey-1)*nSx*nSy
77fc19a7e9 Mart*0333 #endif /* ALLOW_AUTODIFF_TAMC */
85f5225996 Mart*0334 #ifdef ALLOW_SEAICE_EVP_RESIDUAL
                0335 C     save previous (p-1) iteration
                0336           IF ( printResidual ) THEN
                0337            DO j=1,sNy
                0338             DO i=1,sNx
8377b8ee87 Mart*0339              sig11Pm1(i,j,bi,bj) = seaice_sigma1(i,j,bi,bj)
                0340              sig22Pm1(i,j,bi,bj) = seaice_sigma2(i,j,bi,bj)
                0341              sig12Pm1(i,j,bi,bj) = seaice_sigma12(i,j,bi,bj)
                0342              uIcePm1 (i,j,bi,bj) = uIce(i,j,bi,bj)
                0343              vIcePm1 (i,j,bi,bj) = vIce(i,j,bi,bj)
85f5225996 Mart*0344             ENDDO
                0345            ENDDO
                0346           ENDIF
                0347 #endif /* ALLOW_SEAICE_EVP_RESIDUAL */
9e706e0219 Jean*0348           DO j=1-OLy,sNy+OLy
                0349            DO i=1-OLx,sNx+OLx
8377b8ee87 Mart*0350             seaice_div    (i,j) = 0. _d 0
                0351             seaice_tension(i,j) = 0. _d 0
                0352             seaice_shear  (i,j) = 0. _d 0
                0353             pressC        (i,j) = 0. _d 0
                0354             e12Csq        (i,j) = 0. _d 0
                0355             zetaC         (i,j) = 0. _d 0
                0356             deltaZ        (i,j) = 0. _d 0
                0357             zetaZ   (i,j,bi,bj) = 0. _d 0
                0358             deltaC  (i,j,bi,bj) = 0. _d 0
b7030e3a29 Mart*0359            ENDDO
                0360           ENDDO
9e706e0219 Jean*0361           DO j=1-OLy,sNy+OLy
                0362            DO i=1-OLx,sNx+OLx
da4a30d88f Mart*0363             ep(i,j) = e11(i,j,bi,bj) + e22(i,j,bi,bj)
                0364             em(i,j) = e11(i,j,bi,bj) - e22(i,j,bi,bj)
                0365            ENDDO
                0366           ENDDO
d2082aecba Mart*0367 C     need to do this beforehand for easier vectorization after
                0368 C     TAFization
                0369 C     average strain rates to C points
                0370           IF ( SEAICEetaZmethod .EQ. 0 ) THEN
                0371            DO j=1-OLy+1,sNy+OLy-1
                0372             DO i=1-OLx+1,sNx+OLx-1
                0373              tmp = 0.25 *
8377b8ee87 Mart*0374      &            ( e12(i,j  ,bi,bj) + e12(i+1,j  ,bi,bj)
                0375      &            + e12(i,j+1,bi,bj) + e12(i+1,j+1,bi,bj) )
d2082aecba Mart*0376              e12Csq(i,j) = tmp*tmp
                0377             ENDDO
                0378            ENDDO
                0379           ELSEIF ( SEAICEetaZmethod .EQ. 3 ) THEN
                0380            DO j=1-OLy+1,sNy+OLy-1
                0381             DO i=1-OLx+1,sNx+OLx-1
                0382 C     area weighted average of the squares of e12 is more accurate
85f5225996 Mart*0383 C     (and energy conserving) according to Bouillon et al. 2013, eq (11)
8377b8ee87 Mart*0384              e12Csq(i,j) = 0.25 _d 0 * recip_rA(i,j,bi,bj) *
                0385      &            ( rAz(i  ,j  ,bi,bj)*e12(i  ,j  ,bi,bj)**2
                0386      &            + rAz(i+1,j  ,bi,bj)*e12(i+1,j  ,bi,bj)**2
                0387      &            + rAz(i  ,j+1,bi,bj)*e12(i  ,j+1,bi,bj)**2
                0388      &            + rAz(i+1,j+1,bi,bj)*e12(i+1,j+1,bi,bj)**2 )
d2082aecba Mart*0389             ENDDO
                0390            ENDDO
                0391           ENDIF
68ff576152 Mart*0392           DO j=0,sNy+1
                0393            DO i=0,sNx+1
8377b8ee87 Mart*0394             deltaSq = ep(i,j)**2 + recip_ecc2 * em(i,j)**2
                0395      &           + recip_ecc2 * 4. _d 0 * e12Csq(i,j)
                0396 #ifdef ALLOW_AUTODIFF
d2082aecba Mart*0397 C     avoid sqrt of 0
8377b8ee87 Mart*0398             deltaC(i,j,bi,bj) = 0. _d 0
d2082aecba Mart*0399             IF ( deltaSq .GT. 0. _d 0 )
8377b8ee87 Mart*0400      &           deltaC(i,j,bi,bj) = SQRT(deltaSq)
d2082aecba Mart*0401 #else
8377b8ee87 Mart*0402             deltaC(i,j,bi,bj) = SQRT(deltaSq)
                0403 #endif /* ALLOW_AUTODIFF */
d2082aecba Mart*0404 #ifdef SEAICE_DELTA_SMOOTHREG
                0405 C     smooth regularization (without max-function) of delta for
                0406 C     better differentiability
85f5225996 Mart*0407 CML            deltaCreg  = SQRT(deltaSq + deltaMinSq)
8377b8ee87 Mart*0408             deltaCreg  = deltaC(i,j,bi,bj) + SEAICE_deltaMin
d2082aecba Mart*0409 #else
8377b8ee87 Mart*0410             deltaCreg  = MAX(deltaC(i,j,bi,bj),SEAICE_deltaMin)
d2082aecba Mart*0411 #endif /* SEAICE_DELTA_SMOOTHREG */
8377b8ee87 Mart*0412             zetaC(i,j) = HALF*( press0(i,j,bi,bj)
                0413      &           * ( 1. _d 0 + tensileStrFac(i,j,bi,bj) )
45c6bb8414 Mart*0414      &           )/deltaCreg
d2082aecba Mart*0415            ENDDO
                0416           ENDDO
a6389affa6 Mart*0417           IF ( useAdaptiveEVP ) THEN
                0418            DO j=0,sNy+1
                0419             DO i=0,sNx+1
e4120ef5ff Jean*0420 CML   I do not like these hidden regularisations, why do we need to
85f5225996 Mart*0421 CML   divide by mass?
af61e5eb16 Mart*0422 #ifdef ALLOW_AUTODIFF
                0423              evpAlphaC(i,j,bi,bj) = SEAICE_evpAlpha
                0424              IF ( zetaC(i,j) .GT. 0. _d 0 ) THEN
                0425 #endif
8377b8ee87 Mart*0426              evpAlphaC(i,j,bi,bj) = SQRT(zetaC(i,j)
                0427      &            * EVPcFac / MAX(seaiceMassC(i,j,bi,bj), 1.D-04)
                0428      &            * recip_rA(i,j,bi,bj) ) * HEFFM(i,j,bi,bj)
af61e5eb16 Mart*0429 #ifdef ALLOW_AUTODIFF
                0430              ENDIF
                0431 #endif
8377b8ee87 Mart*0432              evpAlphaC(i,j,bi,bj) =
                0433      &            MAX(evpAlphaC(i,j,bi,bj),SEAICEaEVPalphaMin)
a6389affa6 Mart*0434             ENDDO
85f5225996 Mart*0435            ENDDO
a6389affa6 Mart*0436           ENDIF
e4120ef5ff Jean*0437 C     compute zetaZ and deltaZ by simple averaging (following
d2082aecba Mart*0438 C     Bouillon et al., 2013)
8377b8ee87 Mart*0439           DO j=1,sNy+1
                0440            DO i=1,sNx+1
                0441             sumNorm = HEFFM(i,j,  bi,bj)+HEFFM(i-1,j,  bi,bj)
                0442      &              + HEFFM(i,j-1,bi,bj)+HEFFM(i-1,j-1,bi,bj)
967b0b9f87 Mart*0443             IF ( sumNorm.GT.0. _d 0 ) sumNorm = 1. _d 0 / sumNorm
8377b8ee87 Mart*0444             zetaZ(i,j,bi,bj) = sumNorm *
                0445      &           ( zetaC(i,  j) + zetaC(i-1,j-1)
                0446      &           + zetaC(i-1,j) + zetaC(i,  j-1) )
                0447             deltaZ(i,j) = sumNorm *
                0448      &           ( deltaC(i,  j,bi,bj) + deltaC(i-1,j-1,bi,bj)
                0449      &           + deltaC(i-1,j,bi,bj) + deltaC(i,  j-1,bi,bj) )
d2082aecba Mart*0450            ENDDO
                0451           ENDDO
fc27138d73 Mart*0452 #ifdef SEAICE_ALLOW_CLIPZETA
4619635089 Mart*0453 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0454 CADJ STORE zetac            = comlev1_bibj_evp, key = tkey, byte=isbyte
                0455 CADJ STORE zetaz(:,:,bi,bj) = comlev1_bibj_evp, key = tkey, byte=isbyte
4619635089 Mart*0456 #endif /* ALLOW_AUTODIFF_TAMC */
d8ae85c68c Mart*0457 C     regularize zeta if necessary
                0458           DO j=0,sNy+1
                0459            DO i=0,sNx+1
8e32c48b8f Mart*0460             zetaC(i,j)  = MAX( SEAICE_zMin(i,j,bi,bj),
                0461      &           MIN( SEAICE_zMax(i,j,bi,bj), zetaC(i,j) ) )
ec0d7df165 Mart*0462 CML            zetaC(I,J)   = zetaC(I,J)*HEFFM(I,J,bi,bj)
d8ae85c68c Mart*0463 C
282fbc7b4c Mart*0464 C     zMin, zMax are defined at C-points, make sure that they are not
                0465 C     masked by boundaries/land points
1cd6ce9d27 Mart*0466             zMaxZ       = MAX(
8e32c48b8f Mart*0467      &        MAX(SEAICE_zMax(i,  j,bi,bj),SEAICE_zMax(i,  j-1,bi,bj)),
                0468      &        MAX(SEAICE_zMax(i-1,j,bi,bj),SEAICE_zMax(i-1,j-1,bi,bj)) )
1cd6ce9d27 Mart*0469             zMinZ       = MAX(
8e32c48b8f Mart*0470      &        MAX(SEAICE_zMin(i,  j,bi,bj),SEAICE_zMin(i,  j-1,bi,bj)),
                0471      &        MAX(SEAICE_zMin(i-1,j,bi,bj),SEAICE_zMin(i-1,j-1,bi,bj)) )
8377b8ee87 Mart*0472             zetaZ(i,j,bi,bj) = MAX(zMinZ,MIN(zMaxZ,zetaZ(i,j,bi,bj)))
1cd6ce9d27 Mart*0473            ENDDO
                0474           ENDDO
fc27138d73 Mart*0475 #endif /* SEAICE_ALLOW_CLIPZETA */
d8ae85c68c Mart*0476 C     recompute pressure
                0477           DO j=0,sNy+1
                0478            DO i=0,sNx+1
8377b8ee87 Mart*0479             pressC(i,j) =
                0480      &           ( press0(i,j,bi,bj) * ( 1. _d 0 - SEAICEpressReplFac )
                0481      &           + TWO*zetaC(i,j)*deltaC(i,j,bi,bj)*SEAICEpressReplFac
                0482      &             /(1. _d 0 + tensileStrFac(i,j,bi,bj))
                0483      &           ) * (1. _d 0 - tensileStrFac(i,j,bi,bj))
d8ae85c68c Mart*0484            ENDDO
                0485           ENDDO
0beea9c759 Mart*0486 #ifdef ALLOW_DIAGNOSTICS
1cd6ce9d27 Mart*0487           IF ( useDiagnostics ) THEN
a3f6dcaab3 Mart*0488 C     save eta, zeta, and pressure for diagnostics
1cd6ce9d27 Mart*0489            DO j=1,sNy
                0490             DO i=1,sNx
8377b8ee87 Mart*0491              press(i,j,bi,bj) = pressC(i,j)
                0492              zeta (i,j,bi,bj) = zetaC(i,j)
                0493              eta  (i,j,bi,bj) = zetaC(i,j)*recip_ecc2
1cd6ce9d27 Mart*0494             ENDDO
                0495            ENDDO
                0496           ENDIF
0beea9c759 Mart*0497 #endif /* ALLOW_DIAGNOSTICS */
b7030e3a29 Mart*0498 C     Calculate the RHS of the stress equations. Do this now in order to
                0499 C     avoid multiple divisions by Delta
                0500 C     P * ( D_d/Delta - 1 ) = 2*zeta*D_d - P = 2*zeta*D_d - 2*zeta*Delta
                0501 C     P * ( D_t/Delta     ) = 2*zeta*D_t
                0502 C     P * ( D_s/Delta     ) = 2*zeta*D_s
dadd13178c Mart*0503 #ifdef SEAICE_ALLOW_TEM
1cd6ce9d27 Mart*0504           IF ( SEAICEuseTEM ) THEN
d2082aecba Mart*0505            DO j=0,sNy
                0506             DO i=0,sNx
8377b8ee87 Mart*0507              etaDenC   = em(i,j)**2 + 4. _d 0 * e12Csq(i,j)
d2082aecba Mart*0508              etaDenC  = SQRT(MAX(deltaMinSq,etaDenC))
8377b8ee87 Mart*0509              zetaMaxC = ecc2*zetaC(i,j)
                0510      &            *(deltaC(i,j,bi,bj)-ep(i,j))/etaDenC
d2082aecba Mart*0511 #ifdef ALLOW_DIAGNOSTICS
                0512 C     save new eta for diagnostics
8377b8ee87 Mart*0513              eta(i,j,bi,bj) = MIN(zetaC(i,j),zetaMaxC)*recip_ecc2
d2082aecba Mart*0514 #endif /* ALLOW_DIAGNOSTICS */
8377b8ee87 Mart*0515              seaice_div    (i,j) =
                0516      &            ( 2. _d 0 *zetaC(i,j)*ep(i,j) - pressC(i,j)
                0517      &            ) * HEFFM(i,j,bi,bj)
                0518              seaice_tension(i,j) = 2. _d 0*MIN(zetaC(i,j),zetaMaxC)
                0519      &            * em(i,j) * HEFFM(i,j,bi,bj)
d2082aecba Mart*0520             ENDDO
                0521            ENDDO
                0522            DO j=1,sNy+1
                0523             DO i=1,sNx+1
967b0b9f87 Mart*0524              sumNorm = 0.25 _d 0
                0525 CML            sumNorm = 1.0 _d 0
ec0d7df165 Mart*0526 CML     &           / MAX(1.D0,HEFFM(I,  J,  bi,bj)
                0527 CML     &           +          HEFFM(I-1,J,  bi,bj)
                0528 CML     &           +          HEFFM(I,  J-1,bi,bj)
                0529 CML     &           +          HEFFM(I-1,J-1,bi,bj) )
da4a30d88f Mart*0530 C     Averaging the squares gives more accurate viscous-plastic behavior
                0531 C     than squaring the averages
c9a6744012 Jean*0532              etaDenZ  =
8377b8ee87 Mart*0533      &            sumNorm * recip_rAz(i,j,bi,bj) *
                0534      &                    ( _rA(i  ,j  ,bi,bj) * em(i,  j  )**2
                0535      &                    + _rA(i-1,j-1,bi,bj) * em(i-1,j-1)**2
                0536      &                    + _rA(i-1,j  ,bi,bj) * em(i-1,j  )**2
                0537      &                    + _rA(i  ,j-1,bi,bj) * em(i,  j-1)**2 )
                0538      &                    + 4. _d 0*e12(i,j,bi,bj)**2
3daf25222c Mart*0539              etaDenZ  = SQRT(MAX(deltaMinSq,etaDenZ))
8377b8ee87 Mart*0540              zetaMaxZ = ecc2*zetaZ(i,j,bi,bj) * ( deltaZ(i,j)
                0541      &            - sumNorm * ( ep(i,j  ) + ep(i-1,j  )
                0542      &                        + ep(i,j-1) + ep(i-1,j-1) )
da4a30d88f Mart*0543      &            )/etaDenZ
8377b8ee87 Mart*0544              seaice_shear  (i,j) =
                0545      &            2. _d 0*MIN(zetaZ(i,j,bi,bj),zetaMaxZ)
                0546      &            * 2. _d 0*e12(i,j,bi,bj)
1cd6ce9d27 Mart*0547             ENDDO
                0548            ENDDO
                0549           ELSE
                0550 #else
                0551           IF (.TRUE. ) THEN
dadd13178c Mart*0552 #endif /* SEAICE_ALLOW_TEM */
8377b8ee87 Mart*0553            DO j=0,sNy
                0554             DO i=0,sNx
                0555              seaice_div    (i,j) =
                0556      &            ( 2. _d 0 *zetaC(i,j)*ep(i,j) - pressC(i,j)
                0557      &            ) * HEFFM(i,j,bi,bj)
                0558              seaice_tension(i,j) = 2. _d 0*zetaC(i,j)
                0559      &            * em(i,j) * HEFFM(i,j,bi,bj)
d2082aecba Mart*0560             ENDDO
                0561            ENDDO
8377b8ee87 Mart*0562            DO j=1,sNy+1
                0563             DO i=1,sNx+1
                0564              seaice_shear  (i,j) =
                0565      &            2. _d 0*zetaZ(i,j,bi,bj)*e12(i,j,bi,bj)
1cd6ce9d27 Mart*0566             ENDDO
037351a1a6 Mart*0567            ENDDO
1cd6ce9d27 Mart*0568           ENDIF
af61e5eb16 Mart*0569 #ifdef ALLOW_AUTODIFF_TAMC
                0570 CADJ STORE evpAlphaC(:,:,bi,bj) =comlev1_bibj_evp,key=tkey,byte=isbyte
                0571 CADJ STORE seaice_div        = comlev1_bibj_evp, key=tkey, byte=isbyte
                0572 CADJ STORE seaice_tension    = comlev1_bibj_evp, key=tkey, byte=isbyte
                0573 CADJ STORE seaice_shear      = comlev1_bibj_evp, key=tkey, byte=isbyte
                0574 #endif /* ALLOW_AUTODIFF_TAMC */
1fb8ac1204 Mart*0575 C
b7030e3a29 Mart*0576 C     first step stress equations
1fb8ac1204 Mart*0577 C
a6389affa6 Mart*0578           IF ( useAdaptiveEVP ) THEN
                0579            DO j=0,sNy
                0580             DO i=0,sNx
8377b8ee87 Mart*0581              denom1(i,j,bi,bj) = 1. _d 0 /  evpAlphaC(i,j,bi,bj)
                0582              denom2(i,j,bi,bj) = denom1(i,j,bi,bj)
a6389affa6 Mart*0583             ENDDO
                0584            ENDDO
                0585           ENDIF
af61e5eb16 Mart*0586 #ifdef ALLOW_AUTODIFF_TAMC
                0587 CADJ STORE denom1(:,:,bi,bj) = comlev1_bibj_evp, key=tkey, byte=isbyte
                0588 CADJ STORE denom2(:,:,bi,bj) = comlev1_bibj_evp, key=tkey, byte=isbyte
                0589 #endif /* ALLOW_AUTODIFF_TAMC */
1fb8ac1204 Mart*0590           DO j=0,sNy
                0591            DO i=0,sNx
b7030e3a29 Mart*0592 C     sigma1 and sigma2 are computed on C points
8377b8ee87 Mart*0593             seaice_sigma1 (i,j,bi,bj) = ( seaice_sigma1 (i,j,bi,bj)
                0594      &           * ( evpAlphaC(i,j,bi,bj) - evpRevFac )
                0595      &           + seaice_div(i,j)
                0596      &           ) * denom1(i,j,bi,bj)
                0597      &           *HEFFM(i,j,bi,bj)
                0598             seaice_sigma2 (i,j,bi,bj) = ( seaice_sigma2 (i,j,bi,bj)
                0599      &           * ( evpAlphaC(i,j,bi,bj) - evpRevFac )
                0600      &           + seaice_tension(i,j)*recip_evpRevFac
                0601      &           ) * denom2(i,j,bi,bj)
                0602      &         *HEFFM(i,j,bi,bj)
dc7b968f4a Mart*0603 #ifdef SEAICE_EVP_ELIMINATE_UNDERFLOWS
                0604 C     Code to avoid very small numbers that can degrade performance.
9e706e0219 Jean*0605 C     Many compilers can handle this more efficiently with the help of
dc7b968f4a Mart*0606 C     a flag (copied from CICE after correspondence with Elizabeth Hunke)
8377b8ee87 Mart*0607             seaice_sigma1(i,j,bi,bj) = SIGN(MAX(
                0608      &           ABS( seaice_sigma1(i,j,bi,bj) ), SEAICE_EPS ),
                0609      &           seaice_sigma1(i,j,bi,bj) )
                0610             seaice_sigma2(i,j,bi,bj) = SIGN(MAX(
                0611      &           ABS( seaice_sigma2(i,j,bi,bj) ), SEAICE_EPS ),
                0612      &           seaice_sigma2(i,j,bi,bj) )
dc7b968f4a Mart*0613 #endif /* SEAICE_EVP_ELIMINATE_UNDERFLOWS */
b7030e3a29 Mart*0614 C     recover sigma11 and sigma22
8377b8ee87 Mart*0615             sig11(i,j) = 0.5 _d 0 *
                0616      &           ( seaice_sigma1(i,j,bi,bj)+seaice_sigma2(i,j,bi,bj) )
                0617             sig22(i,j) = 0.5 _d 0 *
                0618      &           ( seaice_sigma1(i,j,bi,bj)-seaice_sigma2(i,j,bi,bj) )
b7030e3a29 Mart*0619            ENDDO
1fb8ac1204 Mart*0620           ENDDO
                0621 
                0622 C     sigma12 is computed on Z points
a6389affa6 Mart*0623           IF ( useAdaptiveEVP ) THEN
                0624            DO j=1,sNy+1
                0625             DO i=1,sNx+1
8377b8ee87 Mart*0626              evpAlphaZ(i,j,bi,bj) = 0.25 _d 0 *
                0627      &            ( evpAlphaC(i,  j,bi,bj)+evpAlphaC(i-1,j-1,bi,bj)
                0628      &            + evpAlphaC(i-1,j,bi,bj)+evpAlphaC(i,  j-1,bi,bj) )
                0629              denom2(i,j,bi,bj) = 1. _d 0 /  evpAlphaZ(i,j,bi,bj)
a6389affa6 Mart*0630             ENDDO
                0631            ENDDO
                0632           ENDIF
af61e5eb16 Mart*0633 #ifdef ALLOW_AUTODIFF_TAMC
                0634 CADJ STORE denom2   (:,:,bi,bj) = comlev1_bibj_evp,key=tkey,byte=isbyte
                0635 CADJ STORE evpAlphaZ(:,:,bi,bj) = comlev1_bibj_evp,key=tkey,byte=isbyte
                0636 #endif /* ALLOW_AUTODIFF_TAMC */
1fb8ac1204 Mart*0637           DO j=1,sNy+1
                0638            DO i=1,sNx+1
8377b8ee87 Mart*0639             seaice_sigma12(i,j,bi,bj) = ( seaice_sigma12(i,j,bi,bj)
                0640      &           * ( evpAlphaZ(i,j,bi,bj) - evpRevFac )
                0641      &           + seaice_shear(i,j)*recip_evpRevFac
                0642      &           ) * denom2(i,j,bi,bj)
dc7b968f4a Mart*0643 #ifdef SEAICE_EVP_ELIMINATE_UNDERFLOWS
8377b8ee87 Mart*0644             seaice_sigma12(i,j,bi,bj) = SIGN(MAX(
                0645      &           ABS( seaice_sigma12(i,j,bi,bj) ), SEAICE_EPS ),
                0646      &           seaice_sigma12(i,j,bi,bj) )
dc7b968f4a Mart*0647 #endif /* SEAICE_EVP_ELIMINATE_UNDERFLOWS */
1fb8ac1204 Mart*0648            ENDDO
b4949dd6db Jean*0649           ENDDO
b7030e3a29 Mart*0650 C
                0651 C     compute divergence of stress tensor
                0652 C
8377b8ee87 Mart*0653           DO j=1,sNy
                0654            DO i=1,sNx
                0655             stressDivergenceX(i,j,bi,bj) =
                0656      &           ( sig11(i  ,j  ) * _dyF(i  ,j,bi,bj)
                0657      &           - sig11(i-1,j  ) * _dyF(i-1,j,bi,bj)
                0658      &           + seaice_sigma12(i,j+1,bi,bj) * _dxV(i,j+1,bi,bj)
                0659      &           - seaice_sigma12(i,j  ,bi,bj) * _dxV(i,j  ,bi,bj)
                0660      &           ) * recip_rAw(i,j,bi,bj)
                0661             stressDivergenceY(i,j,bi,bj) =
                0662      &           ( sig22(i,j  ) * _dxF(i,j  ,bi,bj)
                0663      &           - sig22(i,j-1) * _dxF(i,j-1,bi,bj)
                0664      &           + seaice_sigma12(i+1,j,bi,bj) * _dyU(i+1,j,bi,bj)
                0665      &           - seaice_sigma12(i  ,j,bi,bj) * _dyU(i  ,j,bi,bj)
                0666      &           ) * recip_rAs(i,j,bi,bj)
b7030e3a29 Mart*0667            ENDDO
b4949dd6db Jean*0668           ENDDO
4669aaa4e4 Mart*0669          ENDDO
                0670         ENDDO
                0671 
e4120ef5ff Jean*0672 #ifdef ALLOW_SEAICE_EVP_RESIDUAL
4669aaa4e4 Mart*0673         IF ( printResidual ) THEN
                0674          DO bj=myByLo(myThid),myByHi(myThid)
                0675           DO bi=myBxLo(myThid),myBxHi(myThid)
85f5225996 Mart*0676            resTile(bi,bj) = 0. _d 0
                0677            DO j=1,sNy
                0678             DO i=1,sNx
e4120ef5ff Jean*0679              sig11Pm1(i,j,bi,bj)  =
4669aaa4e4 Mart*0680      &            seaice_sigma1(i,j,bi,bj)-sig11pm1(i,j,bi,bj)
e4120ef5ff Jean*0681              sig22Pm1(i,j,bi,bj)  =
4669aaa4e4 Mart*0682      &            seaice_sigma2(i,j,bi,bj)-sig22pm1(i,j,bi,bj)
e4120ef5ff Jean*0683              sig12Pm1(i,j,bi,bj)  =
4669aaa4e4 Mart*0684      &            seaice_sigma12(i,j,bi,bj)-sig12Pm1(i,j,bi,bj)
e4120ef5ff Jean*0685              sig11Pm1(i,j,bi,bj)  =
8377b8ee87 Mart*0686      &            evpAlphaC(i,j,bi,bj) * sig11Pm1(i,j,bi,bj)
e4120ef5ff Jean*0687              sig22Pm1(i,j,bi,bj)  =
8377b8ee87 Mart*0688      &            evpAlphaC(i,j,bi,bj) * sig22Pm1(i,j,bi,bj)
e4120ef5ff Jean*0689              sig12Pm1(i,j,bi,bj)  =
8377b8ee87 Mart*0690      &            evpAlphaZ(i,j,bi,bj) * sig12Pm1(i,j,bi,bj)
85f5225996 Mart*0691             ENDDO
                0692            ENDDO
859b587af3 Mart*0693            IF ( .NOT. SEAICEscaleSurfStress ) THEN
e4120ef5ff Jean*0694 C     multiply with mask (concentration) to only count ice contributions
859b587af3 Mart*0695             DO j=1,sNy
                0696              DO i=1,sNx
                0697               resTile(bi,bj) = resTile(bi,bj) + AREA(i,j,bi,bj) *
                0698      &             ( sig11Pm1(i,j,bi,bj)*sig11Pm1(i,j,bi,bj)
                0699      &             + sig22Pm1(i,j,bi,bj)*sig22Pm1(i,j,bi,bj)
                0700      &             + sig12Pm1(i,j,bi,bj)*sig12Pm1(i,j,bi,bj) )
                0701              ENDDO
                0702             ENDDO
                0703            ELSE
                0704 C     in this case the scaling with AREA is already done
                0705             DO j=1,sNy
                0706              DO i=1,sNx
                0707               resTile(bi,bj) = resTile(bi,bj)
                0708      &             + sig11Pm1(i,j,bi,bj)*sig11Pm1(i,j,bi,bj)
                0709      &             + sig22Pm1(i,j,bi,bj)*sig22Pm1(i,j,bi,bj)
                0710      &             + sig12Pm1(i,j,bi,bj)*sig12Pm1(i,j,bi,bj)
                0711              ENDDO
                0712             ENDDO
                0713            ENDIF
4669aaa4e4 Mart*0714           ENDDO
88601fb2cb Gael*0715          ENDDO
85f5225996 Mart*0716          resloc = 0. _d 0
                0717          CALL GLOBAL_SUM_TILE_RL( resTile, resloc, myThid )
                0718          resloc = SQRT(resloc)
                0719          WRITE(standardMessageUnit,'(A,1X,I6,1PE16.8)')
a6389affa6 Mart*0720      &   ' SEAICE_EVP: iEVPstep, residual sigma = ', iEVPstep, resLoc
85f5225996 Mart*0721         ENDIF
                0722 #endif /* ALLOW_SEAICE_EVP_RESIDUAL */
88601fb2cb Gael*0723 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0724 CADJ STORE stressDivergenceX = comlev1_evp, key = evpkey, byte=isbyte
                0725 CADJ STORE stressDivergenceY = comlev1_evp, key = evpkey, byte=isbyte
a6389affa6 Mart*0726 # ifdef SEAICE_DYN_STABLE_ADJOINT
c9a6744012 Jean*0727 Cgf zero out adjoint fields to stabilize pkg/seaice dyna. adjoint
4f95e6bec9 Gael*0728       CALL ZERO_ADJ( 1, stressDivergenceX, myThid)
                0729       CALL ZERO_ADJ( 1, stressDivergenceY, myThid)
a6389affa6 Mart*0730 # endif /* SEAICE_DYN_STABLE_ADJOINT */
4f95e6bec9 Gael*0731 #endif /* ALLOW_AUTODIFF_TAMC */
                0732 
6e2f4e58fa Mart*0733 C
037351a1a6 Mart*0734 C     set up rhs for stepping the velocity field
6e2f4e58fa Mart*0735 C
3f31c7d5de Mart*0736         CALL SEAICE_OCEANDRAG_COEFFS(
ec0d7df165 Mart*0737      I       uIce, vIce, HEFFM,
3f31c7d5de Mart*0738      O       DWATN,
                0739      I       iEVPstep, myTime, myIter, myThid )
df1dac8b7b Mart*0740 #ifdef SEAICE_ALLOW_BOTTOMDRAG
d5254b4e3d Mart*0741         CALL SEAICE_BOTTOMDRAG_COEFFS(
ec0d7df165 Mart*0742      I       uIce, vIce, HEFFM,
df1dac8b7b Mart*0743 #ifdef SEAICE_ITD
                0744      I       HEFFITD, AREAITD, AREA,
                0745 #else
                0746      I       HEFF, AREA,
e4120ef5ff Jean*0747 #endif
df1dac8b7b Mart*0748      O       CbotC,
                0749      I       iEVPstep, myTime, myIter, myThid )
                0750 #endif /* SEAICE_ALLOW_BOTTOMDRAG */
af61e5eb16 Mart*0751 #ifdef ALLOW_AUTODIFF_TAMC
                0752 CADJ STORE DWATN = comlev1_evp, key = evpkey, byte=isbyte
                0753 # ifdef SEAICE_ALLOW_BOTTOMDRAG
                0754 CADJ STORE CbotC = comlev1_evp, key = evpkey, byte=isbyte
                0755 # endif
                0756 #endif
88601fb2cb Gael*0757         DO bj=myByLo(myThid),myByHi(myThid)
                0758          DO bi=myBxLo(myThid),myBxHi(myThid)
af61e5eb16 Mart*0759 #ifdef ALLOW_AUTODIFF_TAMC
                0760           tkey = bi + (bj-1)*nSx + (evpkey-1)*nSx*nSy
                0761 CADJ STORE uIce(:,:,bi,bj) = comlev1_bibj_evp, key=tkey, byte=isbyte
                0762 CADJ STORE vIce(:,:,bi,bj) = comlev1_bibj_evp, key=tkey, byte=isbyte
                0763 #endif /* ALLOW_AUTODIFF_TAMC */
8377b8ee87 Mart*0764           DO j=1,sNy
                0765            DO i=1,sNx
e8b26f9f39 Mart*0766 C     over open water, all terms that contain sea ice mass drop out and
                0767 C     the balance is determined by the atmosphere-ice and ice-ocean stress;
e4120ef5ff Jean*0768 C     the staggering of uIce and vIce can cause stripes in the open ocean
                0769 C     solution when the turning angle is non-zero (SINWAT.NE.0);
e8b26f9f39 Mart*0770 C     we mask this term here in order to avoid the stripes and because
                0771 C     over open ocean, u/vIce do not advect anything, so that the associated
b8e39f4242 Mart*0772 C     error is small and most likely only confined to the ice edge but may
e8b26f9f39 Mart*0773 C     propagate into the ice covered regions.
8377b8ee87 Mart*0774             locMaskU = seaiceMassU(i,j,bi,bj)
                0775             locMaskV = seaiceMassV(i,j,bi,bj)
e8b26f9f39 Mart*0776             IF ( locMaskU .NE. 0. _d 0 ) locMaskU = 1. _d 0
                0777             IF ( locMaskV .NE. 0. _d 0 ) locMaskV = 1. _d 0
                0778 C     to recover old results replace the above lines with the line below
                0779 C           locMaskU = 1. _d 0
                0780 C           locMaskV = 1. _d 0
037351a1a6 Mart*0781 C     set up anti symmetric drag force and add in ice ocean stress
6e2f4e58fa Mart*0782 C     ( remember to average to correct velocity points )
8377b8ee87 Mart*0783             FORCEX(i,j,bi,bj)=FORCEX0(i,j,bi,bj)+
                0784      &           ( 0.5 _d 0 * ( DWATN(i,j,bi,bj)+DWATN(i-1,j,bi,bj) ) *
                0785      &           COSWAT * uVel(i,j,kSrf,bi,bj)
                0786      &           - SIGN(SINWAT, _fCori(i,j,bi,bj))* 0.5 _d 0 *
                0787      &           ( DWATN(i  ,j,bi,bj) * 0.5 _d 0 *
                0788      &            (vVel(i  ,j  ,kSrf,bi,bj)-vIce(i  ,j  ,bi,bj)
                0789      &            +vVel(i  ,j+1,kSrf,bi,bj)-vIce(i  ,j+1,bi,bj))
                0790      &           + DWATN(i-1,j,bi,bj) * 0.5 _d 0 *
                0791      &            (vVel(i-1,j  ,kSrf,bi,bj)-vIce(i-1,j  ,bi,bj)
                0792      &            +vVel(i-1,j+1,kSrf,bi,bj)-vIce(i-1,j+1,bi,bj))
                0793      &           )*locMaskU ) * areaW(i,j,bi,bj)
                0794             FORCEY(i,j,bi,bj)=FORCEY0(i,j,bi,bj)+
                0795      &           ( 0.5 _d 0 * ( DWATN(i,j,bi,bj)+DWATN(i,j-1,bi,bj) ) *
                0796      &           COSWAT * vVel(i,j,kSrf,bi,bj)
                0797      &           + SIGN(SINWAT, _fCori(i,j,bi,bj)) * 0.5 _d 0 *
                0798      &           ( DWATN(i,j  ,bi,bj) * 0.5 _d 0 *
                0799      &            (uVel(i  ,j  ,kSrf,bi,bj)-uIce(i  ,j  ,bi,bj)
                0800      &            +uVel(i+1,j  ,kSrf,bi,bj)-uIce(i+1,j  ,bi,bj))
                0801      &           + DWATN(i,j-1,bi,bj) * 0.5 _d 0 *
                0802      &            (uVel(i  ,j-1,kSrf,bi,bj)-uIce(i  ,j-1,bi,bj)
                0803      &            +uVel(i+1,j-1,kSrf,bi,bj)-uIce(i+1,j-1,bi,bj))
                0804      &           )*locMaskV ) * areaS(i,j,bi,bj)
037351a1a6 Mart*0805 C     coriols terms
8377b8ee87 Mart*0806             FORCEX(i,j,bi,bj)=FORCEX(i,j,bi,bj) + 0.5 _d 0*(
                0807      &             seaiceMassC(i  ,j,bi,bj) * _fCori(i  ,j,bi,bj)
                0808      &           * 0.5 _d 0*( vIce(i  ,j,bi,bj)+vIce(i  ,j+1,bi,bj) )
                0809      &           + seaiceMassC(i-1,j,bi,bj) * _fCori(i-1,j,bi,bj)
                0810      &           * 0.5 _d 0*( vIce(i-1,j,bi,bj)+vIce(i-1,j+1,bi,bj) )
b7030e3a29 Mart*0811      &           )
8377b8ee87 Mart*0812             FORCEY(i,j,bi,bj)=FORCEY(i,j,bi,bj) - 0.5 _d 0*(
                0813      &             seaiceMassC(i,j  ,bi,bj) * _fCori(i,j  ,bi,bj)
                0814      &           * 0.5 _d 0*( uIce(i,j  ,bi,bj)+uIce(i+1,  j,bi,bj) )
                0815      &           + seaiceMassC(i,j-1,bi,bj) * _fCori(i,j-1,bi,bj)
                0816      &           * 0.5 _d 0*( uIce(i,j-1,bi,bj)+uIce(i+1,j-1,bi,bj) )
b7030e3a29 Mart*0817      &           )
452bde753c Mart*0818            ENDDO
                0819           ENDDO
af61e5eb16 Mart*0820 #ifdef ALLOW_AUTODIFF_TAMC
                0821 CADJ STORE FORCEX(:,:,bi,bj) = comlev1_bibj_evp, key=tkey, byte=isbyte
                0822 CADJ STORE FORCEY(:,:,bi,bj) = comlev1_bibj_evp, key=tkey, byte=isbyte
                0823 #endif /* ALLOW_AUTODIFF_TAMC */
8ce59fd29e Mart*0824 #ifdef SEAICE_ALLOW_MOM_ADVECTION
ec0d7df165 Mart*0825           IF ( SEAICEmomAdvection ) THEN
8377b8ee87 Mart*0826            DO j=1-OLy,sNy+OLy
                0827             DO i=1-OLx,sNx+OLx
                0828              gUmom(i,j) = 0. _d 0
                0829              gVmom(i,j) = 0. _d 0
8ce59fd29e Mart*0830             ENDDO
                0831            ENDDO
                0832            CALL SEAICE_MOM_ADVECTION(
                0833      I          bi,bj,1,sNx,1,sNy,
                0834      I          uIce, vIce,
                0835      O          gUmom, gVmom,
                0836      I          myTime, myIter, myThid )
8377b8ee87 Mart*0837            DO j=1,sNy
                0838             DO i=1,sNx
                0839              FORCEX(i,j,bi,bj) = FORCEX(i,j,bi,bj) + gUmom(i,j)
                0840              FORCEY(i,j,bi,bj) = FORCEY(i,j,bi,bj) + gVmom(i,j)
8ce59fd29e Mart*0841             ENDDO
                0842            ENDDO
                0843           ENDIF
                0844 #endif /* SEAICE_ALLOW_MOM_ADVECTION */
6e2f4e58fa Mart*0845 C
                0846 C     step momentum equations with ice-ocean stress term treated implicitly
                0847 C
a6389affa6 Mart*0848           IF ( useAdaptiveEVP ) THEN
8377b8ee87 Mart*0849            DO j=1,sNy
                0850             DO i=1,sNx
85f5225996 Mart*0851 C     compute and adjust parameters that are constant otherwise
8377b8ee87 Mart*0852              evpBetaU(i,j,bi,bj) = 0.5 _d 0*(evpAlphaC(i-1,j,bi,bj)
                0853      &                                      +evpAlphaC(i,  j,bi,bj))
                0854              evpBetaV(i,j,bi,bj) = 0.5 _d 0*(evpAlphaC(i,j-1,bi,bj)
                0855      &                                      +evpAlphaC(i,j,  bi,bj))
85f5225996 Mart*0856 
a6389affa6 Mart*0857             ENDDO
                0858            ENDDO
                0859           ENDIF
af61e5eb16 Mart*0860 #ifdef ALLOW_AUTODIFF_TAMC
                0861 CADJ STORE evpBetaU(:,:,bi,bj) = comlev1_bibj_evp,key=tkey,byte=isbyte
                0862 CADJ STORE evpBetaV(:,:,bi,bj) = comlev1_bibj_evp,key=tkey,byte=isbyte
                0863 CADJ STORE uIce(:,:,bi,bj) = comlev1_bibj_evp, key=tkey, byte=isbyte
                0864 CADJ STORE vIce(:,:,bi,bj) = comlev1_bibj_evp, key=tkey, byte=isbyte
                0865 #endif /* ALLOW_AUTODIFF_TAMC */
8377b8ee87 Mart*0866           DO j=1,sNy
                0867            DO i=1,sNx
                0868             betaFacU   = evpBetaU(i,j,bi,bj)*recip_deltaT
                0869             betaFacV   = evpBetaV(i,j,bi,bj)*recip_deltaT
85f5225996 Mart*0870             tmp        = evpStarFac*recip_deltaT
                0871             betaFacP1V = betaFacV + tmp
                0872             betaFacP1U = betaFacU + tmp
8377b8ee87 Mart*0873             denomU = seaiceMassU(i,j,bi,bj)*betaFacP1U
                0874      &           + 0.5 _d 0*( DWATN(i,j,bi,bj) + DWATN(i-1,j,bi,bj) )
                0875      &           * COSWAT * areaW(i,j,bi,bj)
                0876             denomV = seaiceMassV(i,j,bi,bj)*betaFacP1V
                0877      &           + 0.5 _d 0*( DWATN(i,j,bi,bj) + DWATN(i,j-1,bi,bj) )
                0878      &           * COSWAT * areaS(i,j,bi,bj)
df1dac8b7b Mart*0879 #ifdef SEAICE_ALLOW_BOTTOMDRAG
8377b8ee87 Mart*0880             denomU = denomU + areaW(i,j,bi,bj)
                0881      &           * 0.5 _d 0*( CbotC(i,j,bi,bj) + CbotC(i-1,j,bi,bj) )
                0882             denomV = denomV + areaS(i,j,bi,bj)
                0883      &           * 0.5 _d 0*( CbotC(i,j,bi,bj) + CbotC(i,j-1,bi,bj) )
df1dac8b7b Mart*0884 #endif /* SEAICE_ALLOW_BOTTOMDRAG */
70e078b38a Mart*0885             IF ( denomU .EQ. 0. _d 0 ) denomU = 1. _d 0
                0886             IF ( denomV .EQ. 0. _d 0 ) denomV = 1. _d 0
8377b8ee87 Mart*0887             uIce(i,j,bi,bj) = seaiceMaskU(i,j,bi,bj) *
                0888      &           ( seaiceMassU(i,j,bi,bj)*betaFacU
                0889      &           * uIce(i,j,bi,bj)
                0890      &           + seaiceMassU(i,j,bi,bj)*recip_deltaT*evpStarFac
                0891      &           * uIceNm1(i,j,bi,bj)
                0892      &           + FORCEX(i,j,bi,bj)
                0893      &           + stressDivergenceX(i,j,bi,bj) ) / denomU
                0894             vIce(i,j,bi,bj) = seaiceMaskV(i,j,bi,bj) *
                0895      &           ( seaiceMassV(i,j,bi,bj)*betaFacV
                0896      &           * vIce(i,j,bi,bj)
                0897      &           + seaiceMassV(i,j,bi,bj)*recip_deltaT*evpStarFac
                0898      &           * vIceNm1(i,j,bi,bj)
                0899      &           + FORCEY(i,j,bi,bj)
                0900      &           + stressDivergenceY(i,j,bi,bj) ) / denomV
c9a6744012 Jean*0901 C--   to change results  of lab_sea.hb87 test exp. (only preserve 2 digits for cg2d)
a3f6dcaab3 Mart*0902 c           uIce(i,j,bi,bj) = uIceNm1(i,j,bi,bj)
                0903 c    &                       +( uIce(i,j,bi,bj) - uIceNm1(i,j,bi,bj) )
                0904 c           vIce(i,j,bi,bj) = vIceNm1(i,j,bi,bj)
                0905 c    &                       +( vIce(i,j,bi,bj) - vIceNm1(i,j,bi,bj) )
037351a1a6 Mart*0906            ENDDO
                0907           ENDDO
fec75090d3 Jean*0908 #ifndef OBCS_UVICE_OLD
af61e5eb16 Mart*0909 #ifdef ALLOW_AUTODIFF_TAMC
                0910 CADJ STORE uIce(:,:,bi,bj) = comlev1_bibj_evp, key=tkey, byte=isbyte
                0911 CADJ STORE vIce(:,:,bi,bj) = comlev1_bibj_evp, key=tkey, byte=isbyte
                0912 #endif /* ALLOW_AUTODIFF_TAMC */
fec75090d3 Jean*0913           DO j=1,sNy
                0914            DO i=1,sNx
                0915             locMaskU = maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj)
                0916             locMaskV = maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj)
a3f6dcaab3 Mart*0917             uIce(i,j,bi,bj) = uIce(i,j,bi,bj)*locMaskU
                0918      &                       + uIceNm1(i,j,bi,bj)*(ONE-locMaskU)
                0919             vIce(i,j,bi,bj) = vIce(i,j,bi,bj)*locMaskV
                0920      &                       + vIceNm1(i,j,bi,bj)*(ONE-locMaskV)
fec75090d3 Jean*0921            ENDDO
                0922           ENDDO
                0923 #endif /* OBCS_UVICE_OLD */
6e2f4e58fa Mart*0924          ENDDO
                0925         ENDDO
b4949dd6db Jean*0926 
a3f6dcaab3 Mart*0927         CALL EXCH_UV_XY_RL(uIce,vIce,.TRUE.,myThid)
b4949dd6db Jean*0928 
e4120ef5ff Jean*0929 #ifdef ALLOW_SEAICE_EVP_RESIDUAL
85f5225996 Mart*0930         IF ( printResidual ) THEN
                0931          DO bj=myByLo(myThid),myByHi(myThid)
                0932           DO bi=myBxLo(myThid),myBxHi(myThid)
                0933            resTile(bi,bj) = 0. _d 0
8377b8ee87 Mart*0934            DO j=1,sNy
                0935             DO i=1,sNx
                0936              uIcePm1(i,j,bi,bj) = seaiceMaskU(i,j,bi,bj) *
                0937      &            ( uIce(i,j,bi,bj)-uIcePm1(i,j,bi,bj) )
                0938      &            * evpBetaU(i,j,bi,bj)
                0939              vIcePm1(i,j,bi,bj) = seaiceMaskV(i,j,bi,bj) *
                0940      &            ( vIce(i,j,bi,bj)-vIcePm1(i,j,bi,bj) )
                0941      &            * evpBetaV(i,j,bi,bj)
85f5225996 Mart*0942             ENDDO
                0943            ENDDO
859b587af3 Mart*0944            IF ( .NOT. SEAICEscaleSurfStress ) THEN
e4120ef5ff Jean*0945 C     multiply with mask (concentration) to only count ice contributions
859b587af3 Mart*0946             DO j=1,sNy
                0947              DO i=1,sNx
                0948               resTile(bi,bj) = resTile(bi,bj) + AREA(i,j,bi,bj) *
8377b8ee87 Mart*0949      &             ( uIcePm1(i,j,bi,bj)*uIcePm1(i,j,bi,bj)
                0950      &             + vIcePm1(i,j,bi,bj)*vIcePm1(i,j,bi,bj) )
859b587af3 Mart*0951              ENDDO
                0952             ENDDO
                0953            ELSE
                0954             DO j=1,sNy
                0955              DO i=1,sNx
                0956               resTile(bi,bj) = resTile(bi,bj)
8377b8ee87 Mart*0957      &             + uIcePm1(i,j,bi,bj)*uIcePm1(i,j,bi,bj)
                0958      &             + vIcePm1(i,j,bi,bj)*vIcePm1(i,j,bi,bj)
859b587af3 Mart*0959              ENDDO
                0960             ENDDO
                0961            ENDIF
85f5225996 Mart*0962           ENDDO
                0963          ENDDO
                0964          resloc = 0. _d 0
                0965          CALL GLOBAL_SUM_TILE_RL( resTile, resloc, myThid )
                0966          resloc = SQRT(resloc)
                0967          WRITE(standardMessageUnit,'(A,1X,I6,1PE16.8)')
a6389affa6 Mart*0968      &        ' SEAICE_EVP: iEVPstep, residual U = ', iEVPstep, resLoc
85f5225996 Mart*0969         ENDIF
                0970 CML        WRITE(suff,'(I10.10)') myIter*100000+iEvpStep
                0971 CML        CALL WRITE_FLD_XY_RL( 'DELTA.',suff,deltaC,
                0972 CML     &       myIter*100000+iEvpStep,myThid)
                0973 CML        CALL WRITE_FLD_XY_RL( 'RSIG1.',suff,sig11pm1,
                0974 CML     &       myIter*100000+iEvpStep,myThid)
                0975 CML        CALL WRITE_FLD_XY_RL( 'RSIG2.',suff,sig22pm1,
                0976 CML     &       myIter*100000+iEvpStep,myThid)
                0977 CML        CALL WRITE_FLD_XY_RL( 'RSIG12.',suff,sig12pm1,
                0978 CML     &       myIter*100000+iEvpStep,myThid)
                0979 CML        CALL WRITE_FLD_XY_RL( 'RUICE.',suff,uIcePm1,
                0980 CML     &       myIter*100000+iEvpStep,myThid)
                0981 CML        CALL WRITE_FLD_XY_RL( 'RVICE.',suff,vIcePm1,
                0982 CML     &       myIter*100000+iEvpStep,myThid)
                0983 
                0984 #endif /* ALLOW_SEAICE_EVP_RESIDUAL */
                0985 
cf9fa44a59 Patr*0986        ENDIF
6e2f4e58fa Mart*0987 C     end of the main time loop
b4949dd6db Jean*0988       ENDDO
6e2f4e58fa Mart*0989 
45315406aa Mart*0990 #endif /* SEAICE_CGRID and SEAICE_ALLOW_EVP */
6e2f4e58fa Mart*0991 
                0992       RETURN
                0993       END