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