pkg/seaice/seaice_evp.F

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
                0316 CADJ STORE evpAlphaC      = comlev1_evp, key = evpkey, byte = isbyte
                0317 CADJ STORE evpBetaU       = comlev1_evp, key = evpkey, byte = isbyte
                0318 CADJ STORE evpBetaV       = comlev1_evp, key = evpkey, byte = isbyte
6e2f4e58fa Mart*0319 #endif /* ALLOW_AUTODIFF_TAMC */
                0320 C
                0321 C     first calculate strain rates and bulk moduli/viscosities
b4949dd6db Jean*0322 C
                0323         CALL SEAICE_CALC_STRAINRATES(
a3f6dcaab3 Mart*0324      I       uIce, vIce,
037351a1a6 Mart*0325      O       e11, e22, e12,
2e75568507 Mart*0326      I       iEVPstep, myTime, myIter, myThid )
037351a1a6 Mart*0327
88601fb2cb Gael*0328 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0329 CADJ STORE e11 = comlev1_evp,key = evpkey, byte = isbyte
                0330 CADJ STORE e12 = comlev1_evp,key = evpkey, byte = isbyte
                0331 CADJ STORE e22 = comlev1_evp,key = evpkey, byte = isbyte
88601fb2cb Gael*0332 #endif /* ALLOW_AUTODIFF_TAMC */
                0333
037351a1a6 Mart*0334         DO bj=myByLo(myThid),myByHi(myThid)
                0335          DO bi=myBxLo(myThid),myBxHi(myThid)
77fc19a7e9 Mart*0336 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0337           tkey = bi + (bj-1)*nSx + (evpkey-1)*nSx*nSy
77fc19a7e9 Mart*0338 #endif /* ALLOW_AUTODIFF_TAMC */
85f5225996 Mart*0339 #ifdef ALLOW_SEAICE_EVP_RESIDUAL
                0340 C     save previous (p-1) iteration
                0341           IF ( printResidual ) THEN
                0342            DO j=1,sNy
                0343             DO i=1,sNx
8377b8ee87 Mart*0344              sig11Pm1(i,j,bi,bj) = seaice_sigma1(i,j,bi,bj)
                0345              sig22Pm1(i,j,bi,bj) = seaice_sigma2(i,j,bi,bj)
                0346              sig12Pm1(i,j,bi,bj) = seaice_sigma12(i,j,bi,bj)
                0347              uIcePm1 (i,j,bi,bj) = uIce(i,j,bi,bj)
                0348              vIcePm1 (i,j,bi,bj) = vIce(i,j,bi,bj)
85f5225996 Mart*0349             ENDDO
                0350            ENDDO
                0351           ENDIF
                0352 #endif /* ALLOW_SEAICE_EVP_RESIDUAL */
9e706e0219 Jean*0353           DO j=1-OLy,sNy+OLy
                0354            DO i=1-OLx,sNx+OLx
8377b8ee87 Mart*0355             seaice_div    (i,j) = 0. _d 0
                0356             seaice_tension(i,j) = 0. _d 0
                0357             seaice_shear  (i,j) = 0. _d 0
                0358             pressC        (i,j) = 0. _d 0
                0359             e12Csq        (i,j) = 0. _d 0
                0360             zetaC         (i,j) = 0. _d 0
                0361             deltaZ        (i,j) = 0. _d 0
                0362             zetaZ   (i,j,bi,bj) = 0. _d 0
                0363             deltaC  (i,j,bi,bj) = 0. _d 0
b7030e3a29 Mart*0364            ENDDO
                0365           ENDDO
9e706e0219 Jean*0366           DO j=1-OLy,sNy+OLy
                0367            DO i=1-OLx,sNx+OLx
da4a30d88f Mart*0368             ep(i,j) = e11(i,j,bi,bj) + e22(i,j,bi,bj)
                0369             em(i,j) = e11(i,j,bi,bj) - e22(i,j,bi,bj)
                0370            ENDDO
                0371           ENDDO
d2082aecba Mart*0372 C     need to do this beforehand for easier vectorization after
                0373 C     TAFization
                0374 C     average strain rates to C points
                0375           IF ( SEAICEetaZmethod .EQ. 0 ) THEN
                0376            DO j=1-OLy+1,sNy+OLy-1
                0377             DO i=1-OLx+1,sNx+OLx-1
                0378              tmp = 0.25 *
8377b8ee87 Mart*0379      &            ( e12(i,j  ,bi,bj) + e12(i+1,j  ,bi,bj)
                0380      &            + e12(i,j+1,bi,bj) + e12(i+1,j+1,bi,bj) )
d2082aecba Mart*0381              e12Csq(i,j) = tmp*tmp
                0382             ENDDO
                0383            ENDDO
                0384           ELSEIF ( SEAICEetaZmethod .EQ. 3 ) THEN
                0385            DO j=1-OLy+1,sNy+OLy-1
                0386             DO i=1-OLx+1,sNx+OLx-1
                0387 C     area weighted average of the squares of e12 is more accurate
85f5225996 Mart*0388 C     (and energy conserving) according to Bouillon et al. 2013, eq (11)
8377b8ee87 Mart*0389              e12Csq(i,j) = 0.25 _d 0 * recip_rA(i,j,bi,bj) *
                0390      &            ( rAz(i  ,j  ,bi,bj)*e12(i  ,j  ,bi,bj)**2
                0391      &            + rAz(i+1,j  ,bi,bj)*e12(i+1,j  ,bi,bj)**2
                0392      &            + rAz(i  ,j+1,bi,bj)*e12(i  ,j+1,bi,bj)**2
                0393      &            + rAz(i+1,j+1,bi,bj)*e12(i+1,j+1,bi,bj)**2 )
d2082aecba Mart*0394             ENDDO
                0395            ENDDO
                0396           ENDIF
68ff576152 Mart*0397           DO j=0,sNy+1
                0398            DO i=0,sNx+1
8377b8ee87 Mart*0399             deltaSq = ep(i,j)**2 + recip_ecc2 * em(i,j)**2
                0400      &           + recip_ecc2 * 4. _d 0 * e12Csq(i,j)
                0401 #ifdef ALLOW_AUTODIFF
d2082aecba Mart*0402 C     avoid sqrt of 0
8377b8ee87 Mart*0403             deltaC(i,j,bi,bj) = 0. _d 0
d2082aecba Mart*0404             IF ( deltaSq .GT. 0. _d 0 )
8377b8ee87 Mart*0405      &           deltaC(i,j,bi,bj) = SQRT(deltaSq)
d2082aecba Mart*0406 #else
8377b8ee87 Mart*0407             deltaC(i,j,bi,bj) = SQRT(deltaSq)
                0408 #endif /* ALLOW_AUTODIFF */
d2082aecba Mart*0409 #ifdef SEAICE_DELTA_SMOOTHREG
                0410 C     smooth regularization (without max-function) of delta for
                0411 C     better differentiability
85f5225996 Mart*0412 CML            deltaCreg  = SQRT(deltaSq + deltaMinSq)
8377b8ee87 Mart*0413             deltaCreg  = deltaC(i,j,bi,bj) + SEAICE_deltaMin
d2082aecba Mart*0414 #else
8377b8ee87 Mart*0415             deltaCreg  = MAX(deltaC(i,j,bi,bj),SEAICE_deltaMin)
d2082aecba Mart*0416 #endif /* SEAICE_DELTA_SMOOTHREG */
8377b8ee87 Mart*0417             zetaC(i,j) = HALF*( press0(i,j,bi,bj)
                0418      &           * ( 1. _d 0 + tensileStrFac(i,j,bi,bj) )
45c6bb8414 Mart*0419      &           )/deltaCreg
d2082aecba Mart*0420            ENDDO
                0421           ENDDO
a6389affa6 Mart*0422           IF ( useAdaptiveEVP ) THEN
                0423            DO j=0,sNy+1
                0424             DO i=0,sNx+1
e4120ef5ff Jean*0425 CML   I do not like these hidden regularisations, why do we need to
85f5225996 Mart*0426 CML   divide by mass?
8377b8ee87 Mart*0427              evpAlphaC(i,j,bi,bj) = SQRT(zetaC(i,j)
                0428      &            * EVPcFac / MAX(seaiceMassC(i,j,bi,bj), 1.D-04)
                0429      &            * recip_rA(i,j,bi,bj) ) * HEFFM(i,j,bi,bj)
                0430              evpAlphaC(i,j,bi,bj) =
                0431      &            MAX(evpAlphaC(i,j,bi,bj),SEAICEaEVPalphaMin)
a6389affa6 Mart*0432             ENDDO
85f5225996 Mart*0433            ENDDO
a6389affa6 Mart*0434           ENDIF
e4120ef5ff Jean*0435 C     compute zetaZ and deltaZ by simple averaging (following
d2082aecba Mart*0436 C     Bouillon et al., 2013)
8377b8ee87 Mart*0437           DO j=1,sNy+1
                0438            DO i=1,sNx+1
                0439             sumNorm = HEFFM(i,j,  bi,bj)+HEFFM(i-1,j,  bi,bj)
                0440      &              + HEFFM(i,j-1,bi,bj)+HEFFM(i-1,j-1,bi,bj)
967b0b9f87 Mart*0441             IF ( sumNorm.GT.0. _d 0 ) sumNorm = 1. _d 0 / sumNorm
8377b8ee87 Mart*0442             zetaZ(i,j,bi,bj) = sumNorm *
                0443      &           ( zetaC(i,  j) + zetaC(i-1,j-1)
                0444      &           + zetaC(i-1,j) + zetaC(i,  j-1) )
                0445             deltaZ(i,j) = sumNorm *
                0446      &           ( deltaC(i,  j,bi,bj) + deltaC(i-1,j-1,bi,bj)
                0447      &           + deltaC(i-1,j,bi,bj) + deltaC(i,  j-1,bi,bj) )
d2082aecba Mart*0448            ENDDO
                0449           ENDDO
fc27138d73 Mart*0450 #ifdef SEAICE_ALLOW_CLIPZETA
4619635089 Mart*0451 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0452 CADJ STORE zetac            = comlev1_bibj_evp, key = tkey, byte=isbyte
                0453 CADJ STORE zetaz(:,:,bi,bj) = comlev1_bibj_evp, key = tkey, byte=isbyte
4619635089 Mart*0454 #endif /* ALLOW_AUTODIFF_TAMC */
d8ae85c68c Mart*0455 C     regularize zeta if necessary
                0456           DO j=0,sNy+1
                0457            DO i=0,sNx+1
8e32c48b8f Mart*0458             zetaC(i,j)  = MAX( SEAICE_zMin(i,j,bi,bj),
                0459      &           MIN( SEAICE_zMax(i,j,bi,bj), zetaC(i,j) ) )
ec0d7df165 Mart*0460 CML            zetaC(I,J)   = zetaC(I,J)*HEFFM(I,J,bi,bj)
d8ae85c68c Mart*0461 C
282fbc7b4c Mart*0462 C     zMin, zMax are defined at C-points, make sure that they are not
                0463 C     masked by boundaries/land points
1cd6ce9d27 Mart*0464             zMaxZ       = MAX(
8e32c48b8f Mart*0465      &        MAX(SEAICE_zMax(i,  j,bi,bj),SEAICE_zMax(i,  j-1,bi,bj)),
                0466      &        MAX(SEAICE_zMax(i-1,j,bi,bj),SEAICE_zMax(i-1,j-1,bi,bj)) )
1cd6ce9d27 Mart*0467             zMinZ       = MAX(
8e32c48b8f Mart*0468      &        MAX(SEAICE_zMin(i,  j,bi,bj),SEAICE_zMin(i,  j-1,bi,bj)),
                0469      &        MAX(SEAICE_zMin(i-1,j,bi,bj),SEAICE_zMin(i-1,j-1,bi,bj)) )
8377b8ee87 Mart*0470             zetaZ(i,j,bi,bj) = MAX(zMinZ,MIN(zMaxZ,zetaZ(i,j,bi,bj)))
1cd6ce9d27 Mart*0471            ENDDO
                0472           ENDDO
fc27138d73 Mart*0473 #endif /* SEAICE_ALLOW_CLIPZETA */
d8ae85c68c Mart*0474 C     recompute pressure
                0475           DO j=0,sNy+1
                0476            DO i=0,sNx+1
8377b8ee87 Mart*0477             pressC(i,j) =
                0478      &           ( press0(i,j,bi,bj) * ( 1. _d 0 - SEAICEpressReplFac )
                0479      &           + TWO*zetaC(i,j)*deltaC(i,j,bi,bj)*SEAICEpressReplFac
                0480      &             /(1. _d 0 + tensileStrFac(i,j,bi,bj))
                0481      &           ) * (1. _d 0 - tensileStrFac(i,j,bi,bj))
d8ae85c68c Mart*0482            ENDDO
                0483           ENDDO
0beea9c759 Mart*0484 #ifdef ALLOW_DIAGNOSTICS
1cd6ce9d27 Mart*0485           IF ( useDiagnostics ) THEN
a3f6dcaab3 Mart*0486 C     save eta, zeta, and pressure for diagnostics
1cd6ce9d27 Mart*0487            DO j=1,sNy
                0488             DO i=1,sNx
8377b8ee87 Mart*0489              press(i,j,bi,bj) = pressC(i,j)
                0490              zeta (i,j,bi,bj) = zetaC(i,j)
                0491              eta  (i,j,bi,bj) = zetaC(i,j)*recip_ecc2
1cd6ce9d27 Mart*0492             ENDDO
                0493            ENDDO
                0494           ENDIF
0beea9c759 Mart*0495 #endif /* ALLOW_DIAGNOSTICS */
b7030e3a29 Mart*0496 C     Calculate the RHS of the stress equations. Do this now in order to
                0497 C     avoid multiple divisions by Delta
                0498 C     P * ( D_d/Delta - 1 ) = 2*zeta*D_d - P = 2*zeta*D_d - 2*zeta*Delta
                0499 C     P * ( D_t/Delta     ) = 2*zeta*D_t
                0500 C     P * ( D_s/Delta     ) = 2*zeta*D_s
dadd13178c Mart*0501 #ifdef SEAICE_ALLOW_TEM
1cd6ce9d27 Mart*0502           IF ( SEAICEuseTEM ) THEN
d2082aecba Mart*0503            DO j=0,sNy
                0504             DO i=0,sNx
8377b8ee87 Mart*0505              etaDenC   = em(i,j)**2 + 4. _d 0 * e12Csq(i,j)
d2082aecba Mart*0506              etaDenC  = SQRT(MAX(deltaMinSq,etaDenC))
8377b8ee87 Mart*0507              zetaMaxC = ecc2*zetaC(i,j)
                0508      &            *(deltaC(i,j,bi,bj)-ep(i,j))/etaDenC
d2082aecba Mart*0509 #ifdef ALLOW_DIAGNOSTICS
                0510 C     save new eta for diagnostics
8377b8ee87 Mart*0511              eta(i,j,bi,bj) = MIN(zetaC(i,j),zetaMaxC)*recip_ecc2
d2082aecba Mart*0512 #endif /* ALLOW_DIAGNOSTICS */
8377b8ee87 Mart*0513              seaice_div    (i,j) =
                0514      &            ( 2. _d 0 *zetaC(i,j)*ep(i,j) - pressC(i,j)
                0515      &            ) * HEFFM(i,j,bi,bj)
                0516              seaice_tension(i,j) = 2. _d 0*MIN(zetaC(i,j),zetaMaxC)
                0517      &            * em(i,j) * HEFFM(i,j,bi,bj)
d2082aecba Mart*0518             ENDDO
                0519            ENDDO
                0520            DO j=1,sNy+1
                0521             DO i=1,sNx+1
967b0b9f87 Mart*0522              sumNorm = 0.25 _d 0
                0523 CML            sumNorm = 1.0 _d 0
ec0d7df165 Mart*0524 CML     &           / MAX(1.D0,HEFFM(I,  J,  bi,bj)
                0525 CML     &           +          HEFFM(I-1,J,  bi,bj)
                0526 CML     &           +          HEFFM(I,  J-1,bi,bj)
                0527 CML     &           +          HEFFM(I-1,J-1,bi,bj) )
da4a30d88f Mart*0528 C     Averaging the squares gives more accurate viscous-plastic behavior
                0529 C     than squaring the averages
c9a6744012 Jean*0530              etaDenZ  =
8377b8ee87 Mart*0531      &            sumNorm * recip_rAz(i,j,bi,bj) *
                0532      &                    ( _rA(i  ,j  ,bi,bj) * em(i,  j  )**2
                0533      &                    + _rA(i-1,j-1,bi,bj) * em(i-1,j-1)**2
                0534      &                    + _rA(i-1,j  ,bi,bj) * em(i-1,j  )**2
                0535      &                    + _rA(i  ,j-1,bi,bj) * em(i,  j-1)**2 )
                0536      &                    + 4. _d 0*e12(i,j,bi,bj)**2
3daf25222c Mart*0537              etaDenZ  = SQRT(MAX(deltaMinSq,etaDenZ))
8377b8ee87 Mart*0538              zetaMaxZ = ecc2*zetaZ(i,j,bi,bj) * ( deltaZ(i,j)
                0539      &            - sumNorm * ( ep(i,j  ) + ep(i-1,j  )
                0540      &                        + ep(i,j-1) + ep(i-1,j-1) )
da4a30d88f Mart*0541      &            )/etaDenZ
8377b8ee87 Mart*0542              seaice_shear  (i,j) =
                0543      &            2. _d 0*MIN(zetaZ(i,j,bi,bj),zetaMaxZ)
                0544      &            * 2. _d 0*e12(i,j,bi,bj)
1cd6ce9d27 Mart*0545             ENDDO
                0546            ENDDO
                0547           ELSE
                0548 #else
                0549           IF (.TRUE. ) THEN
dadd13178c Mart*0550 #endif /* SEAICE_ALLOW_TEM */
8377b8ee87 Mart*0551            DO j=0,sNy
                0552             DO i=0,sNx
                0553              seaice_div    (i,j) =
                0554      &            ( 2. _d 0 *zetaC(i,j)*ep(i,j) - pressC(i,j)
                0555      &            ) * HEFFM(i,j,bi,bj)
                0556              seaice_tension(i,j) = 2. _d 0*zetaC(i,j)
                0557      &            * em(i,j) * HEFFM(i,j,bi,bj)
d2082aecba Mart*0558             ENDDO
                0559            ENDDO
8377b8ee87 Mart*0560            DO j=1,sNy+1
                0561             DO i=1,sNx+1
                0562              seaice_shear  (i,j) =
                0563      &            2. _d 0*zetaZ(i,j,bi,bj)*e12(i,j,bi,bj)
1cd6ce9d27 Mart*0564             ENDDO
037351a1a6 Mart*0565            ENDDO
1cd6ce9d27 Mart*0566           ENDIF
1fb8ac1204 Mart*0567 C
b7030e3a29 Mart*0568 C     first step stress equations
1fb8ac1204 Mart*0569 C
a6389affa6 Mart*0570 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0571 CADJ STORE denom1(:,:,bi,bj) = comlev1_bibj_evp, key=tkey, byte=isbyte
                0572 CADJ STORE denom2(:,:,bi,bj) = comlev1_bibj_evp, key=tkey, byte=isbyte
a6389affa6 Mart*0573 #endif /* ALLOW_AUTODIFF_TAMC */
                0574           IF ( useAdaptiveEVP ) THEN
                0575            DO j=0,sNy
                0576             DO i=0,sNx
8377b8ee87 Mart*0577              denom1(i,j,bi,bj) = 1. _d 0 /  evpAlphaC(i,j,bi,bj)
                0578              denom2(i,j,bi,bj) = denom1(i,j,bi,bj)
a6389affa6 Mart*0579             ENDDO
                0580            ENDDO
                0581           ENDIF
1fb8ac1204 Mart*0582           DO j=0,sNy
                0583            DO i=0,sNx
b7030e3a29 Mart*0584 C     sigma1 and sigma2 are computed on C points
8377b8ee87 Mart*0585             seaice_sigma1 (i,j,bi,bj) = ( seaice_sigma1 (i,j,bi,bj)
                0586      &           * ( evpAlphaC(i,j,bi,bj) - evpRevFac )
                0587      &           + seaice_div(i,j)
                0588      &           ) * denom1(i,j,bi,bj)
                0589      &           *HEFFM(i,j,bi,bj)
                0590             seaice_sigma2 (i,j,bi,bj) = ( seaice_sigma2 (i,j,bi,bj)
                0591      &           * ( evpAlphaC(i,j,bi,bj) - evpRevFac )
                0592      &           + seaice_tension(i,j)*recip_evpRevFac
                0593      &           ) * denom2(i,j,bi,bj)
                0594      &         *HEFFM(i,j,bi,bj)
dc7b968f4a Mart*0595 #ifdef SEAICE_EVP_ELIMINATE_UNDERFLOWS
                0596 C     Code to avoid very small numbers that can degrade performance.
9e706e0219 Jean*0597 C     Many compilers can handle this more efficiently with the help of
dc7b968f4a Mart*0598 C     a flag (copied from CICE after correspondence with Elizabeth Hunke)
8377b8ee87 Mart*0599             seaice_sigma1(i,j,bi,bj) = SIGN(MAX(
                0600      &           ABS( seaice_sigma1(i,j,bi,bj) ), SEAICE_EPS ),
                0601      &           seaice_sigma1(i,j,bi,bj) )
                0602             seaice_sigma2(i,j,bi,bj) = SIGN(MAX(
                0603      &           ABS( seaice_sigma2(i,j,bi,bj) ), SEAICE_EPS ),
                0604      &           seaice_sigma2(i,j,bi,bj) )
dc7b968f4a Mart*0605 #endif /* SEAICE_EVP_ELIMINATE_UNDERFLOWS */
b7030e3a29 Mart*0606 C     recover sigma11 and sigma22
8377b8ee87 Mart*0607             sig11(i,j) = 0.5 _d 0 *
                0608      &           ( seaice_sigma1(i,j,bi,bj)+seaice_sigma2(i,j,bi,bj) )
                0609             sig22(i,j) = 0.5 _d 0 *
                0610      &           ( seaice_sigma1(i,j,bi,bj)-seaice_sigma2(i,j,bi,bj) )
b7030e3a29 Mart*0611            ENDDO
1fb8ac1204 Mart*0612           ENDDO
                0613
77fc19a7e9 Mart*0614 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0615 CADJ STORE denom2   (:,:,bi,bj) = comlev1_bibj_evp,key=tkey,byte=isbyte
                0616 CADJ STORE evpAlphaZ(:,:,bi,bj) = comlev1_bibj_evp,key=tkey,byte=isbyte
77fc19a7e9 Mart*0617 #endif /* ALLOW_AUTODIFF_TAMC */
1fb8ac1204 Mart*0618 C     sigma12 is computed on Z points
a6389affa6 Mart*0619           IF ( useAdaptiveEVP ) THEN
                0620            DO j=1,sNy+1
                0621             DO i=1,sNx+1
8377b8ee87 Mart*0622              evpAlphaZ(i,j,bi,bj) = 0.25 _d 0 *
                0623      &            ( evpAlphaC(i,  j,bi,bj)+evpAlphaC(i-1,j-1,bi,bj)
                0624      &            + evpAlphaC(i-1,j,bi,bj)+evpAlphaC(i,  j-1,bi,bj) )
                0625              denom2(i,j,bi,bj) = 1. _d 0 /  evpAlphaZ(i,j,bi,bj)
a6389affa6 Mart*0626             ENDDO
                0627            ENDDO
                0628           ENDIF
1fb8ac1204 Mart*0629           DO j=1,sNy+1
                0630            DO i=1,sNx+1
8377b8ee87 Mart*0631             seaice_sigma12(i,j,bi,bj) = ( seaice_sigma12(i,j,bi,bj)
                0632      &           * ( evpAlphaZ(i,j,bi,bj) - evpRevFac )
                0633      &           + seaice_shear(i,j)*recip_evpRevFac
                0634      &           ) * denom2(i,j,bi,bj)
dc7b968f4a Mart*0635 #ifdef SEAICE_EVP_ELIMINATE_UNDERFLOWS
8377b8ee87 Mart*0636             seaice_sigma12(i,j,bi,bj) = SIGN(MAX(
                0637      &           ABS( seaice_sigma12(i,j,bi,bj) ), SEAICE_EPS ),
                0638      &           seaice_sigma12(i,j,bi,bj) )
dc7b968f4a Mart*0639 #endif /* SEAICE_EVP_ELIMINATE_UNDERFLOWS */
1fb8ac1204 Mart*0640            ENDDO
b4949dd6db Jean*0641           ENDDO
b7030e3a29 Mart*0642 C
                0643 C     compute divergence of stress tensor
                0644 C
8377b8ee87 Mart*0645           DO j=1,sNy
                0646            DO i=1,sNx
                0647             stressDivergenceX(i,j,bi,bj) =
                0648      &           ( sig11(i  ,j  ) * _dyF(i  ,j,bi,bj)
                0649      &           - sig11(i-1,j  ) * _dyF(i-1,j,bi,bj)
                0650      &           + seaice_sigma12(i,j+1,bi,bj) * _dxV(i,j+1,bi,bj)
                0651      &           - seaice_sigma12(i,j  ,bi,bj) * _dxV(i,j  ,bi,bj)
                0652      &           ) * recip_rAw(i,j,bi,bj)
                0653             stressDivergenceY(i,j,bi,bj) =
                0654      &           ( sig22(i,j  ) * _dxF(i,j  ,bi,bj)
                0655      &           - sig22(i,j-1) * _dxF(i,j-1,bi,bj)
                0656      &           + seaice_sigma12(i+1,j,bi,bj) * _dyU(i+1,j,bi,bj)
                0657      &           - seaice_sigma12(i  ,j,bi,bj) * _dyU(i  ,j,bi,bj)
                0658      &           ) * recip_rAs(i,j,bi,bj)
b7030e3a29 Mart*0659            ENDDO
b4949dd6db Jean*0660           ENDDO
4669aaa4e4 Mart*0661          ENDDO
                0662         ENDDO
                0663
e4120ef5ff Jean*0664 #ifdef ALLOW_SEAICE_EVP_RESIDUAL
4669aaa4e4 Mart*0665         IF ( printResidual ) THEN
                0666          DO bj=myByLo(myThid),myByHi(myThid)
                0667           DO bi=myBxLo(myThid),myBxHi(myThid)
85f5225996 Mart*0668            resTile(bi,bj) = 0. _d 0
                0669            DO j=1,sNy
                0670             DO i=1,sNx
e4120ef5ff Jean*0671              sig11Pm1(i,j,bi,bj)  =
4669aaa4e4 Mart*0672      &            seaice_sigma1(i,j,bi,bj)-sig11pm1(i,j,bi,bj)
e4120ef5ff Jean*0673              sig22Pm1(i,j,bi,bj)  =
4669aaa4e4 Mart*0674      &            seaice_sigma2(i,j,bi,bj)-sig22pm1(i,j,bi,bj)
e4120ef5ff Jean*0675              sig12Pm1(i,j,bi,bj)  =
4669aaa4e4 Mart*0676      &            seaice_sigma12(i,j,bi,bj)-sig12Pm1(i,j,bi,bj)
e4120ef5ff Jean*0677              sig11Pm1(i,j,bi,bj)  =
8377b8ee87 Mart*0678      &            evpAlphaC(i,j,bi,bj) * sig11Pm1(i,j,bi,bj)
e4120ef5ff Jean*0679              sig22Pm1(i,j,bi,bj)  =
8377b8ee87 Mart*0680      &            evpAlphaC(i,j,bi,bj) * sig22Pm1(i,j,bi,bj)
e4120ef5ff Jean*0681              sig12Pm1(i,j,bi,bj)  =
8377b8ee87 Mart*0682      &            evpAlphaZ(i,j,bi,bj) * sig12Pm1(i,j,bi,bj)
85f5225996 Mart*0683             ENDDO
                0684            ENDDO
859b587af3 Mart*0685            IF ( .NOT. SEAICEscaleSurfStress ) THEN
e4120ef5ff Jean*0686 C     multiply with mask (concentration) to only count ice contributions
859b587af3 Mart*0687             DO j=1,sNy
                0688              DO i=1,sNx
                0689               resTile(bi,bj) = resTile(bi,bj) + AREA(i,j,bi,bj) *
                0690      &             ( sig11Pm1(i,j,bi,bj)*sig11Pm1(i,j,bi,bj)
                0691      &             + sig22Pm1(i,j,bi,bj)*sig22Pm1(i,j,bi,bj)
                0692      &             + sig12Pm1(i,j,bi,bj)*sig12Pm1(i,j,bi,bj) )
                0693              ENDDO
                0694             ENDDO
                0695            ELSE
                0696 C     in this case the scaling with AREA is already done
                0697             DO j=1,sNy
                0698              DO i=1,sNx
                0699               resTile(bi,bj) = resTile(bi,bj)
                0700      &             + sig11Pm1(i,j,bi,bj)*sig11Pm1(i,j,bi,bj)
                0701      &             + sig22Pm1(i,j,bi,bj)*sig22Pm1(i,j,bi,bj)
                0702      &             + sig12Pm1(i,j,bi,bj)*sig12Pm1(i,j,bi,bj)
                0703              ENDDO
                0704             ENDDO
                0705            ENDIF
4669aaa4e4 Mart*0706           ENDDO
88601fb2cb Gael*0707          ENDDO
85f5225996 Mart*0708          resloc = 0. _d 0
                0709          CALL GLOBAL_SUM_TILE_RL( resTile, resloc, myThid )
                0710          resloc = SQRT(resloc)
                0711          WRITE(standardMessageUnit,'(A,1X,I6,1PE16.8)')
a6389affa6 Mart*0712      &   ' SEAICE_EVP: iEVPstep, residual sigma = ', iEVPstep, resLoc
85f5225996 Mart*0713         ENDIF
                0714 #endif /* ALLOW_SEAICE_EVP_RESIDUAL */
88601fb2cb Gael*0715 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0716 CADJ STORE stressDivergenceX = comlev1_evp, key = evpkey, byte=isbyte
                0717 CADJ STORE stressDivergenceY = comlev1_evp, key = evpkey, byte=isbyte
a6389affa6 Mart*0718 # ifdef SEAICE_DYN_STABLE_ADJOINT
c9a6744012 Jean*0719 Cgf zero out adjoint fields to stabilize pkg/seaice dyna. adjoint
4f95e6bec9 Gael*0720       CALL ZERO_ADJ( 1, stressDivergenceX, myThid)
                0721       CALL ZERO_ADJ( 1, stressDivergenceY, myThid)
a6389affa6 Mart*0722 # endif /* SEAICE_DYN_STABLE_ADJOINT */
4f95e6bec9 Gael*0723 #endif /* ALLOW_AUTODIFF_TAMC */
                0724
6e2f4e58fa Mart*0725 C
037351a1a6 Mart*0726 C     set up rhs for stepping the velocity field
6e2f4e58fa Mart*0727 C
3f31c7d5de Mart*0728         CALL SEAICE_OCEANDRAG_COEFFS(
ec0d7df165 Mart*0729      I       uIce, vIce, HEFFM,
3f31c7d5de Mart*0730      O       DWATN,
                0731      I       iEVPstep, myTime, myIter, myThid )
df1dac8b7b Mart*0732 #ifdef SEAICE_ALLOW_BOTTOMDRAG
d5254b4e3d Mart*0733         CALL SEAICE_BOTTOMDRAG_COEFFS(
ec0d7df165 Mart*0734      I       uIce, vIce, HEFFM,
df1dac8b7b Mart*0735 #ifdef SEAICE_ITD
                0736      I       HEFFITD, AREAITD, AREA,
                0737 #else
                0738      I       HEFF, AREA,
e4120ef5ff Jean*0739 #endif
df1dac8b7b Mart*0740      O       CbotC,
                0741      I       iEVPstep, myTime, myIter, myThid )
                0742 #endif /* SEAICE_ALLOW_BOTTOMDRAG */
3f31c7d5de Mart*0743
88601fb2cb Gael*0744         DO bj=myByLo(myThid),myByHi(myThid)
                0745          DO bi=myBxLo(myThid),myBxHi(myThid)
8377b8ee87 Mart*0746           DO j=1,sNy
                0747            DO i=1,sNx
e8b26f9f39 Mart*0748 C     over open water, all terms that contain sea ice mass drop out and
                0749 C     the balance is determined by the atmosphere-ice and ice-ocean stress;
e4120ef5ff Jean*0750 C     the staggering of uIce and vIce can cause stripes in the open ocean
                0751 C     solution when the turning angle is non-zero (SINWAT.NE.0);
e8b26f9f39 Mart*0752 C     we mask this term here in order to avoid the stripes and because
                0753 C     over open ocean, u/vIce do not advect anything, so that the associated
b8e39f4242 Mart*0754 C     error is small and most likely only confined to the ice edge but may
e8b26f9f39 Mart*0755 C     propagate into the ice covered regions.
8377b8ee87 Mart*0756             locMaskU = seaiceMassU(i,j,bi,bj)
                0757             locMaskV = seaiceMassV(i,j,bi,bj)
e8b26f9f39 Mart*0758             IF ( locMaskU .NE. 0. _d 0 ) locMaskU = 1. _d 0
                0759             IF ( locMaskV .NE. 0. _d 0 ) locMaskV = 1. _d 0
                0760 C     to recover old results replace the above lines with the line below
                0761 C           locMaskU = 1. _d 0
                0762 C           locMaskV = 1. _d 0
037351a1a6 Mart*0763 C     set up anti symmetric drag force and add in ice ocean stress
6e2f4e58fa Mart*0764 C     ( remember to average to correct velocity points )
8377b8ee87 Mart*0765             FORCEX(i,j,bi,bj)=FORCEX0(i,j,bi,bj)+
                0766      &           ( 0.5 _d 0 * ( DWATN(i,j,bi,bj)+DWATN(i-1,j,bi,bj) ) *
                0767      &           COSWAT * uVel(i,j,kSrf,bi,bj)
                0768      &           - SIGN(SINWAT, _fCori(i,j,bi,bj))* 0.5 _d 0 *
                0769      &           ( DWATN(i  ,j,bi,bj) * 0.5 _d 0 *
                0770      &            (vVel(i  ,j  ,kSrf,bi,bj)-vIce(i  ,j  ,bi,bj)
                0771      &            +vVel(i  ,j+1,kSrf,bi,bj)-vIce(i  ,j+1,bi,bj))
                0772      &           + DWATN(i-1,j,bi,bj) * 0.5 _d 0 *
                0773      &            (vVel(i-1,j  ,kSrf,bi,bj)-vIce(i-1,j  ,bi,bj)
                0774      &            +vVel(i-1,j+1,kSrf,bi,bj)-vIce(i-1,j+1,bi,bj))
                0775      &           )*locMaskU ) * areaW(i,j,bi,bj)
                0776             FORCEY(i,j,bi,bj)=FORCEY0(i,j,bi,bj)+
                0777      &           ( 0.5 _d 0 * ( DWATN(i,j,bi,bj)+DWATN(i,j-1,bi,bj) ) *
                0778      &           COSWAT * vVel(i,j,kSrf,bi,bj)
                0779      &           + SIGN(SINWAT, _fCori(i,j,bi,bj)) * 0.5 _d 0 *
                0780      &           ( DWATN(i,j  ,bi,bj) * 0.5 _d 0 *
                0781      &            (uVel(i  ,j  ,kSrf,bi,bj)-uIce(i  ,j  ,bi,bj)
                0782      &            +uVel(i+1,j  ,kSrf,bi,bj)-uIce(i+1,j  ,bi,bj))
                0783      &           + DWATN(i,j-1,bi,bj) * 0.5 _d 0 *
                0784      &            (uVel(i  ,j-1,kSrf,bi,bj)-uIce(i  ,j-1,bi,bj)
                0785      &            +uVel(i+1,j-1,kSrf,bi,bj)-uIce(i+1,j-1,bi,bj))
                0786      &           )*locMaskV ) * areaS(i,j,bi,bj)
037351a1a6 Mart*0787 C     coriols terms
8377b8ee87 Mart*0788             FORCEX(i,j,bi,bj)=FORCEX(i,j,bi,bj) + 0.5 _d 0*(
                0789      &             seaiceMassC(i  ,j,bi,bj) * _fCori(i  ,j,bi,bj)
                0790      &           * 0.5 _d 0*( vIce(i  ,j,bi,bj)+vIce(i  ,j+1,bi,bj) )
                0791      &           + seaiceMassC(i-1,j,bi,bj) * _fCori(i-1,j,bi,bj)
                0792      &           * 0.5 _d 0*( vIce(i-1,j,bi,bj)+vIce(i-1,j+1,bi,bj) )
b7030e3a29 Mart*0793      &           )
8377b8ee87 Mart*0794             FORCEY(i,j,bi,bj)=FORCEY(i,j,bi,bj) - 0.5 _d 0*(
                0795      &             seaiceMassC(i,j  ,bi,bj) * _fCori(i,j  ,bi,bj)
                0796      &           * 0.5 _d 0*( uIce(i,j  ,bi,bj)+uIce(i+1,  j,bi,bj) )
                0797      &           + seaiceMassC(i,j-1,bi,bj) * _fCori(i,j-1,bi,bj)
                0798      &           * 0.5 _d 0*( uIce(i,j-1,bi,bj)+uIce(i+1,j-1,bi,bj) )
b7030e3a29 Mart*0799      &           )
452bde753c Mart*0800            ENDDO
                0801           ENDDO
8ce59fd29e Mart*0802 #ifdef SEAICE_ALLOW_MOM_ADVECTION
ec0d7df165 Mart*0803           IF ( SEAICEmomAdvection ) THEN
8377b8ee87 Mart*0804            DO j=1-OLy,sNy+OLy
                0805             DO i=1-OLx,sNx+OLx
                0806              gUmom(i,j) = 0. _d 0
                0807              gVmom(i,j) = 0. _d 0
8ce59fd29e Mart*0808             ENDDO
                0809            ENDDO
                0810            CALL SEAICE_MOM_ADVECTION(
                0811      I          bi,bj,1,sNx,1,sNy,
                0812      I          uIce, vIce,
                0813      O          gUmom, gVmom,
                0814      I          myTime, myIter, myThid )
8377b8ee87 Mart*0815            DO j=1,sNy
                0816             DO i=1,sNx
                0817              FORCEX(i,j,bi,bj) = FORCEX(i,j,bi,bj) + gUmom(i,j)
                0818              FORCEY(i,j,bi,bj) = FORCEY(i,j,bi,bj) + gVmom(i,j)
8ce59fd29e Mart*0819             ENDDO
                0820            ENDDO
                0821           ENDIF
                0822 #endif /* SEAICE_ALLOW_MOM_ADVECTION */
6e2f4e58fa Mart*0823 C
                0824 C     step momentum equations with ice-ocean stress term treated implicitly
                0825 C
a6389affa6 Mart*0826           IF ( useAdaptiveEVP ) THEN
8377b8ee87 Mart*0827            DO j=1,sNy
                0828             DO i=1,sNx
85f5225996 Mart*0829 C     compute and adjust parameters that are constant otherwise
8377b8ee87 Mart*0830              evpBetaU(i,j,bi,bj) = 0.5 _d 0*(evpAlphaC(i-1,j,bi,bj)
                0831      &                                      +evpAlphaC(i,  j,bi,bj))
                0832              evpBetaV(i,j,bi,bj) = 0.5 _d 0*(evpAlphaC(i,j-1,bi,bj)
                0833      &                                      +evpAlphaC(i,j,  bi,bj))
85f5225996 Mart*0834
a6389affa6 Mart*0835             ENDDO
                0836            ENDDO
                0837           ENDIF
8377b8ee87 Mart*0838           DO j=1,sNy
                0839            DO i=1,sNx
                0840             betaFacU   = evpBetaU(i,j,bi,bj)*recip_deltaT
                0841             betaFacV   = evpBetaV(i,j,bi,bj)*recip_deltaT
85f5225996 Mart*0842             tmp        = evpStarFac*recip_deltaT
                0843             betaFacP1V = betaFacV + tmp
                0844             betaFacP1U = betaFacU + tmp
8377b8ee87 Mart*0845             denomU = seaiceMassU(i,j,bi,bj)*betaFacP1U
                0846      &           + 0.5 _d 0*( DWATN(i,j,bi,bj) + DWATN(i-1,j,bi,bj) )
                0847      &           * COSWAT * areaW(i,j,bi,bj)
                0848             denomV = seaiceMassV(i,j,bi,bj)*betaFacP1V
                0849      &           + 0.5 _d 0*( DWATN(i,j,bi,bj) + DWATN(i,j-1,bi,bj) )
                0850      &           * COSWAT * areaS(i,j,bi,bj)
df1dac8b7b Mart*0851 #ifdef SEAICE_ALLOW_BOTTOMDRAG
8377b8ee87 Mart*0852             denomU = denomU + areaW(i,j,bi,bj)
                0853      &           * 0.5 _d 0*( CbotC(i,j,bi,bj) + CbotC(i-1,j,bi,bj) )
                0854             denomV = denomV + areaS(i,j,bi,bj)
                0855      &           * 0.5 _d 0*( CbotC(i,j,bi,bj) + CbotC(i,j-1,bi,bj) )
df1dac8b7b Mart*0856 #endif /* SEAICE_ALLOW_BOTTOMDRAG */
70e078b38a Mart*0857             IF ( denomU .EQ. 0. _d 0 ) denomU = 1. _d 0
                0858             IF ( denomV .EQ. 0. _d 0 ) denomV = 1. _d 0
8377b8ee87 Mart*0859             uIce(i,j,bi,bj) = seaiceMaskU(i,j,bi,bj) *
                0860      &           ( seaiceMassU(i,j,bi,bj)*betaFacU
                0861      &           * uIce(i,j,bi,bj)
                0862      &           + seaiceMassU(i,j,bi,bj)*recip_deltaT*evpStarFac
                0863      &           * uIceNm1(i,j,bi,bj)
                0864      &           + FORCEX(i,j,bi,bj)
                0865      &           + stressDivergenceX(i,j,bi,bj) ) / denomU
                0866             vIce(i,j,bi,bj) = seaiceMaskV(i,j,bi,bj) *
                0867      &           ( seaiceMassV(i,j,bi,bj)*betaFacV
                0868      &           * vIce(i,j,bi,bj)
                0869      &           + seaiceMassV(i,j,bi,bj)*recip_deltaT*evpStarFac
                0870      &           * vIceNm1(i,j,bi,bj)
                0871      &           + FORCEY(i,j,bi,bj)
                0872      &           + stressDivergenceY(i,j,bi,bj) ) / denomV
c9a6744012 Jean*0873 C--   to change results  of lab_sea.hb87 test exp. (only preserve 2 digits for cg2d)
a3f6dcaab3 Mart*0874 c           uIce(i,j,bi,bj) = uIceNm1(i,j,bi,bj)
                0875 c    &                       +( uIce(i,j,bi,bj) - uIceNm1(i,j,bi,bj) )
                0876 c           vIce(i,j,bi,bj) = vIceNm1(i,j,bi,bj)
                0877 c    &                       +( vIce(i,j,bi,bj) - vIceNm1(i,j,bi,bj) )
037351a1a6 Mart*0878            ENDDO
                0879           ENDDO
fec75090d3 Jean*0880 #ifndef OBCS_UVICE_OLD
                0881           DO j=1,sNy
                0882            DO i=1,sNx
                0883             locMaskU = maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj)
                0884             locMaskV = maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj)
a3f6dcaab3 Mart*0885             uIce(i,j,bi,bj) = uIce(i,j,bi,bj)*locMaskU
                0886      &                       + uIceNm1(i,j,bi,bj)*(ONE-locMaskU)
                0887             vIce(i,j,bi,bj) = vIce(i,j,bi,bj)*locMaskV
                0888      &                       + vIceNm1(i,j,bi,bj)*(ONE-locMaskV)
fec75090d3 Jean*0889            ENDDO
                0890           ENDDO
                0891 #endif /* OBCS_UVICE_OLD */
6e2f4e58fa Mart*0892          ENDDO
                0893         ENDDO
b4949dd6db Jean*0894
88601fb2cb Gael*0895 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0896 CADJ STORE uIce = comlev1_evp, key = evpkey, byte = isbyte
                0897 CADJ STORE vIce = comlev1_evp, key = evpkey, byte = isbyte
88601fb2cb Gael*0898 #endif /* ALLOW_AUTODIFF_TAMC */
                0899
a3f6dcaab3 Mart*0900         CALL EXCH_UV_XY_RL(uIce,vIce,.TRUE.,myThid)
b4949dd6db Jean*0901
e4120ef5ff Jean*0902 #ifdef ALLOW_SEAICE_EVP_RESIDUAL
85f5225996 Mart*0903         IF ( printResidual ) THEN
                0904          DO bj=myByLo(myThid),myByHi(myThid)
                0905           DO bi=myBxLo(myThid),myBxHi(myThid)
                0906            resTile(bi,bj) = 0. _d 0
8377b8ee87 Mart*0907            DO j=1,sNy
                0908             DO i=1,sNx
                0909              uIcePm1(i,j,bi,bj) = seaiceMaskU(i,j,bi,bj) *
                0910      &            ( uIce(i,j,bi,bj)-uIcePm1(i,j,bi,bj) )
                0911      &            * evpBetaU(i,j,bi,bj)
                0912              vIcePm1(i,j,bi,bj) = seaiceMaskV(i,j,bi,bj) *
                0913      &            ( vIce(i,j,bi,bj)-vIcePm1(i,j,bi,bj) )
                0914      &            * evpBetaV(i,j,bi,bj)
85f5225996 Mart*0915             ENDDO
                0916            ENDDO
859b587af3 Mart*0917            IF ( .NOT. SEAICEscaleSurfStress ) THEN
e4120ef5ff Jean*0918 C     multiply with mask (concentration) to only count ice contributions
859b587af3 Mart*0919             DO j=1,sNy
                0920              DO i=1,sNx
                0921               resTile(bi,bj) = resTile(bi,bj) + AREA(i,j,bi,bj) *
8377b8ee87 Mart*0922      &             ( uIcePm1(i,j,bi,bj)*uIcePm1(i,j,bi,bj)
                0923      &             + vIcePm1(i,j,bi,bj)*vIcePm1(i,j,bi,bj) )
859b587af3 Mart*0924              ENDDO
                0925             ENDDO
                0926            ELSE
                0927             DO j=1,sNy
                0928              DO i=1,sNx
                0929               resTile(bi,bj) = resTile(bi,bj)
8377b8ee87 Mart*0930      &             + uIcePm1(i,j,bi,bj)*uIcePm1(i,j,bi,bj)
                0931      &             + vIcePm1(i,j,bi,bj)*vIcePm1(i,j,bi,bj)
859b587af3 Mart*0932              ENDDO
                0933             ENDDO
                0934            ENDIF
85f5225996 Mart*0935           ENDDO
                0936          ENDDO
                0937          resloc = 0. _d 0
                0938          CALL GLOBAL_SUM_TILE_RL( resTile, resloc, myThid )
                0939          resloc = SQRT(resloc)
                0940          WRITE(standardMessageUnit,'(A,1X,I6,1PE16.8)')
a6389affa6 Mart*0941      &        ' SEAICE_EVP: iEVPstep, residual U = ', iEVPstep, resLoc
85f5225996 Mart*0942         ENDIF
                0943 CML        WRITE(suff,'(I10.10)') myIter*100000+iEvpStep
                0944 CML        CALL WRITE_FLD_XY_RL( 'DELTA.',suff,deltaC,
                0945 CML     &       myIter*100000+iEvpStep,myThid)
                0946 CML        CALL WRITE_FLD_XY_RL( 'RSIG1.',suff,sig11pm1,
                0947 CML     &       myIter*100000+iEvpStep,myThid)
                0948 CML        CALL WRITE_FLD_XY_RL( 'RSIG2.',suff,sig22pm1,
                0949 CML     &       myIter*100000+iEvpStep,myThid)
                0950 CML        CALL WRITE_FLD_XY_RL( 'RSIG12.',suff,sig12pm1,
                0951 CML     &       myIter*100000+iEvpStep,myThid)
                0952 CML        CALL WRITE_FLD_XY_RL( 'RUICE.',suff,uIcePm1,
                0953 CML     &       myIter*100000+iEvpStep,myThid)
                0954 CML        CALL WRITE_FLD_XY_RL( 'RVICE.',suff,vIcePm1,
                0955 CML     &       myIter*100000+iEvpStep,myThid)
                0956
                0957 #endif /* ALLOW_SEAICE_EVP_RESIDUAL */
                0958
cf9fa44a59 Patr*0959        ENDIF
6e2f4e58fa Mart*0960 C     end of the main time loop
b4949dd6db Jean*0961       ENDDO
6e2f4e58fa Mart*0962
45315406aa Mart*0963 #endif /* SEAICE_CGRID and SEAICE_ALLOW_EVP */
6e2f4e58fa Mart*0964
                0965       RETURN
                0966       END