pkg/seaice/seaice_lsr.F

53092bcb42 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
53092bcb42 Mart*0010
efcb5efb58 Jean*0011 C--  File seaice_lsr.F: seaice LSR dynamical solver S/R:
                0012 C--   Contents
                0013 C--   o SEAICE_LSR
                0014 C--   o SEAICE_RESIDUAL
5fb12930f9 Mart*0015 C--   o SEAICE_LSR_CALC_COEFFS
bfe6f7447e Mart*0016 C--   o SEAICE_LSR_RHSU
                0017 C--   o SEAICE_LSR_RHSV
5fb12930f9 Mart*0018 C--   o SEAICE_LSR_TRIDIAGU
                0019 C--   o SEAICE_LSR_TRIDIAGV
efcb5efb58 Jean*0020
86fa82a64d Jean*0021 CBOP
e7ae128053 Jean*0022 C     !ROUTINE: SEAICE_LSR
86fa82a64d Jean*0023 C     !INTERFACE:
bb3037b865 Gael*0024       SUBROUTINE SEAICE_LSR( myTime, myIter, myThid )
86fa82a64d Jean*0025
                0026 C     !DESCRIPTION: \bv
                0027 C     *==========================================================*
                0028 C     | S/R SEAICE_LSR
                0029 C     | o Solve ice momentum equation with an LSR dynamics solver
                0030 C     |   (see Zhang and Hibler,   JGR, 102, 8691-8702, 1997
                0031 C     |    and Zhang and Rothrock, MWR, 131,  845- 861, 2003)
                0032 C     |   Written by Jinlun Zhang, PSC/UW, Feb-2001
                0033 C     |                     zhang@apl.washington.edu
efcb5efb58 Jean*0034 C     *==========================================================*
86fa82a64d Jean*0035 C     | C-grid version by Martin Losch
                0036 C     |   Since 2009/03/18: finite-Volume discretization of
                0037 C     |   stress divergence that includes all metric terms
efcb5efb58 Jean*0038 C     *==========================================================*
86fa82a64d Jean*0039 C     \ev
                0040
                0041 C     !USES:
589eca759f Mart*0042       IMPLICIT NONE
                0043
                0044 C     === Global variables ===
                0045 #include "SIZE.h"
                0046 #include "EEPARAMS.h"
                0047 #include "PARAMS.h"
                0048 #include "DYNVARS.h"
                0049 #include "GRID.h"
9b4636ff01 Patr*0050 #include "SEAICE_SIZE.h"
589eca759f Mart*0051 #include "SEAICE_PARAMS.h"
3f0f10fc37 Mart*0052 #include "SEAICE_GRID.h"
03c669d1ab Jean*0053 #include "SEAICE.h"
589eca759f Mart*0054
                0055 #ifdef ALLOW_AUTODIFF_TAMC
                0056 # include "tamc.h"
                0057 #endif
                0058
86fa82a64d Jean*0059 C     !INPUT/OUTPUT PARAMETERS:
589eca759f Mart*0060 C     === Routine arguments ===
86fa82a64d Jean*0061 C     myTime     :: Simulation time
                0062 C     myIter     :: Simulation timestep number
                0063 C     myThid     :: my Thread Id. number
589eca759f Mart*0064       _RL     myTime
                0065       INTEGER myIter
                0066       INTEGER myThid
86fa82a64d Jean*0067 CEOP
589eca759f Mart*0068
                0069 #ifdef SEAICE_CGRID
7b7a25aa58 Jean*0070 C     !FUNCTIONS:
                0071       LOGICAL  DIFFERENT_MULTIPLE
                0072       EXTERNAL DIFFERENT_MULTIPLE
589eca759f Mart*0073
86fa82a64d Jean*0074 C     !LOCAL VARIABLES:
589eca759f Mart*0075 C     === Local variables ===
86fa82a64d Jean*0076 C     i,j,bi,bj  :: Loop counters
589eca759f Mart*0077
4a08d54d3a Mart*0078       INTEGER ipass
5fb12930f9 Mart*0079       INTEGER i, j, m, bi, bj
                0080       INTEGER iMin, iMax, jMin, jMax
79df32c3f1 Mart*0081       INTEGER ICOUNT1, ICOUNT2, linearIterLoc, nonLinIterLoc
486cba1d02 Mart*0082       INTEGER kSrf
7b7a25aa58 Jean*0083       LOGICAL doIterate4u, doIterate4v
143d9ce879 Dami*0084       LOGICAL doNonLinLoop
86fa82a64d Jean*0085       CHARACTER*(MAX_LEN_MBUF) msgBuf
589eca759f Mart*0086
143d9ce879 Dami*0087 #ifdef ALLOW_DIAGNOSTICS
                0088       LOGICAL  DIAGNOSTICS_IS_ON
                0089       EXTERNAL DIAGNOSTICS_IS_ON
                0090       _RL RTmp (1:sNx,1:sNy)
                0091 #endif /* ALLOW_DIAGNOSTICS */
                0092
e45202e340 Mart*0093       _RL WFAU, WFAV, WFAU2, WFAV2
5fb12930f9 Mart*0094       _RL S1, S2, S1A, S2A
e232688b80 Mart*0095       _RL recip_deltaT
589eca759f Mart*0096
143d9ce879 Dami*0097 #ifdef SEAICE_ALLOW_LSR_FLEX
                0098       _RL residIni, residIniNonLin
                0099       _RL LSRflexFac, LSR_ERRORfu, LSR_ERRORfv
                0100       _RL residkm1, residkm2
                0101 #endif /* SEAICE_ALLOW_LSR_FLEX */
                0102
589eca759f Mart*0103 C     diagonals of coefficient matrices
03c669d1ab Jean*0104       _RL AU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0105       _RL BU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0106       _RL CU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0107       _RL AV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0108       _RL BV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0109       _RL CV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
2d6335c1f7 Jean*0110 C     RHS of TriDiag Linear system (both directions => 5 coef A,B,C & Rt1,2)
03c669d1ab Jean*0111       _RL rhsU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0112       _RL rhsV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
2d6335c1f7 Jean*0113 #ifdef SEAICE_GLOBAL_3DIAG_SOLVER
                0114 C     RHS of TriDiag Linear system (1 direction only => 3 coef A,B,C)
                0115       _RL rhsUx(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0116       _RL rhsVy(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0117       INTEGER iSubIter
b2e217c119 Jean*0118       INTEGER errCode
2d6335c1f7 Jean*0119 #endif
7b7a25aa58 Jean*0120 C     coefficients for lateral points, u(j+/-1)
03c669d1ab Jean*0121       _RL uRt1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0122       _RL uRt2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
7b7a25aa58 Jean*0123 C     coefficients for lateral points, v(i+/-1)
03c669d1ab Jean*0124       _RL vRt1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0125       _RL vRt2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
589eca759f Mart*0126 C     abbreviations
03c669d1ab Jean*0127       _RL etaPlusZeta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0128       _RL zetaMinusEta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
df1dac8b7b Mart*0129 C     symmetric drag coefficient
                0130       _RL dragSym     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
bfe6f7447e Mart*0131 C     intermediate time level centered (C) between n and n-1
                0132       _RL uIceC(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0133       _RL vIceC(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
589eca759f Mart*0134 C     auxillary fields
03c669d1ab Jean*0135       _RL uTmp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0136       _RL vTmp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
4a08d54d3a Mart*0137 C     auxillary fields to store intermediate FORCEX/Y0
                0138       _RL fxTmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0139       _RL fyTmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
70e078b38a Mart*0140 C     fractional area at velocity points
733ca03edd Mart*0141       _RL areaW(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0142       _RL areaS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
8ce59fd29e Mart*0143 #ifdef SEAICE_ALLOW_MOM_ADVECTION
                0144 C     tendency due to advection of momentum
                0145       _RL gUmom(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0146       _RL gVmom(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0147 #endif /*  SEAICE_ALLOW_MOM_ADVECTION */
d24acf3cfc Jean*0148 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0149 C     nlkey :: tape key (depends on nonlinear iteration)
                0150 C     lnkey :: tape key (depends on linear and nonlinear iteration)
                0151 C     tkey  :: tape key (depends on tile indices and lnkey)
                0152       INTEGER nlkey, lnkey, tkey
d24acf3cfc Jean*0153 #endif
35fda33b05 Jean*0154       _RL COSWAT
                0155       _RS SINWAT
589eca759f Mart*0156       _RL UERR
63a26f3b3e Mart*0157 #ifdef SEAICE_ALLOW_CHECK_LSR_CONVERGENCE
                0158       _RL resnorm, EKnorm, counter
                0159 #endif /* SEAICE_ALLOW_CHECK_LSR_CONVERGENCE */
efcb5efb58 Jean*0160       LOGICAL printResidual
                0161       _RL residUini, residVini, residUend, residVend
                0162 #ifdef SEAICE_ALLOW_FREEDRIFT
03c669d1ab Jean*0163       _RL uRes (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0164       _RL vRes (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
efcb5efb58 Jean*0165       _RL errIni, errFD, errSum
                0166       _RL residU_fd, residV_fd
                0167       _RL residUmix, residVmix
                0168 #endif /* SEAICE_ALLOW_FREEDRIFT */
0c560a7486 Mart*0169       _RL areaMinLoc
efcb5efb58 Jean*0170
                0171 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0172
                0173       printResidual = debugLevel.GE.debLevA
                0174      &  .AND. DIFFERENT_MULTIPLE( SEAICE_monFreq, myTime, deltaTClock )
589eca759f Mart*0175
a68248c150 Mart*0176 C     extra overlap for (restricted) additive Schwarz method
                0177       jMin = 1   - SEAICE_OLy
                0178       jMax = sNy + SEAICE_OLy
                0179       iMin = 1   - SEAICE_OLx
                0180       iMax = sNx + SEAICE_OLx
cbf501ab81 Jean*0181 C
79df32c3f1 Mart*0182       linearIterLoc = SEAICElinearIterMax
                0183       nonLinIterLoc = SEAICEnonLinIterMax
                0184       IF ( SEAICEusePicardAsPrecon ) THEN
                0185        linearIterLoc = SEAICEpreconlinIter
                0186        nonLinIterLoc = SEAICEpreconNL_Iter
                0187       ENDIF
143d9ce879 Dami*0188
e232688b80 Mart*0189       recip_deltaT = 1. _d 0 / SEAICE_deltaTdyn
5fb12930f9 Mart*0190
143d9ce879 Dami*0191 C     surface level
                0192       IF ( usingPCoords ) THEN
                0193        kSrf = Nr
                0194       ELSE
                0195        kSrf = 1
                0196       ENDIF
                0197 C--   introduce turning angles
                0198       SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)
                0199       COSWAT=COS(SEAICE_waterTurnAngle*deg2rad)
                0200
8377b8ee87 Mart*0201 #ifdef ALLOW_AUTODIFF
d976cd8bb3 Patr*0202 cph break artificial dependencies
                0203       DO bj=myByLo(myThid),myByHi(myThid)
                0204        DO bi=myBxLo(myThid),myBxHi(myThid)
                0205         DO j=1-OLy,sNy+OLy
                0206          DO i=1-OLx,sNx+OLx
a16f976091 Mart*0207           deltaC (i,j,bi,bj) = 0. _d 0
                0208           press  (i,j,bi,bj) = 0. _d 0
7c0da9bbfa Mart*0209           zeta   (i,j,bi,bj) = 0. _d 0
0fd6b2de1a Mart*0210           zetaZ  (i,j,bi,bj) = 0. _d 0
7c0da9bbfa Mart*0211           eta    (i,j,bi,bj) = 0. _d 0
                0212           etaZ   (i,j,bi,bj) = 0. _d 0
                0213           uIceC  (i,j,bi,bj) = 0. _d 0
                0214           vIceC  (i,j,bi,bj) = 0. _d 0
                0215           uIceNm1(i,j,bi,bj) = 0. _d 0
                0216           vIceNm1(i,j,bi,bj) = 0. _d 0
d976cd8bb3 Patr*0217          ENDDO
                0218         ENDDO
                0219        ENDDO
                0220       ENDDO
                0221 #endif
                0222
70e078b38a Mart*0223 C     initialise fractional areas at velocity points
733ca03edd Mart*0224       DO bj=myByLo(myThid),myByHi(myThid)
                0225        DO bi=myBxLo(myThid),myBxHi(myThid)
                0226         DO j=1-OLy,sNy+OLy
                0227          DO i=1-OLx,sNx+OLx
4a08d54d3a Mart*0228           areaW(i,j,bi,bj) = 1. _d 0
                0229           areaS(i,j,bi,bj) = 1. _d 0
                0230           ENDDO
733ca03edd Mart*0231          ENDDO
0c560a7486 Mart*0232         areaMinLoc = 0. _d 0
733ca03edd Mart*0233         IF ( SEAICEscaleSurfStress ) THEN
0c560a7486 Mart*0234          areaMinLoc = 1. _d -10
                0235 CML         areaMinLoc = SEAICE_EPS
4a08d54d3a Mart*0236          DO j=jMin,jMax
                0237           DO i=iMin,iMax
                0238            areaW(i,j,bi,bj) =
                0239      &          0.5 _d 0*(AREA(i,j,bi,bj)+AREA(i-1,j,bi,bj))
                0240            areaS(i,j,bi,bj) =
                0241      &          0.5 _d 0*(AREA(i,j,bi,bj)+AREA(i,j-1,bi,bj))
733ca03edd Mart*0242           ENDDO
                0243          ENDDO
                0244         ENDIF
70e078b38a Mart*0245        ENDDO
                0246       ENDDO
                0247
4a08d54d3a Mart*0248 #ifdef ALLOW_AUTODIFF_TAMC
                0249 CADJ STORE uice    = comlev1, key=ikey_dynamics, kind=isbyte
                0250 CADJ STORE vice    = comlev1, key=ikey_dynamics, kind=isbyte
                0251 #endif /*  ALLOW_AUTODIFF_TAMC */
                0252 C     This is the rhs contribution of the time derivative.
                0253 C     We add it here, because it does not change in the iteration.
                0254 C     We store this update on a local field to avoid touching FORCEX/Y0.
                0255 C     As a side effect, this avoids TAF recomputation warnings.
                0256       DO bj=myByLo(myThid),myByHi(myThid)
                0257        DO bi=myBxLo(myThid),myBxHi(myThid)
                0258         DO j=1-OLy,sNy+OLy
                0259          DO i=1-OLx,sNx+OLx
                0260           uIceNm1(i,j,bi,bj) = uIce(i,j,bi,bj)
                0261           vIceNm1(i,j,bi,bj) = vIce(i,j,bi,bj)
                0262           fxTmp  (i,j,bi,bj) = FORCEX0(i,j,bi,bj)
                0263      &         +seaiceMassU(i,j,bi,bj)*recip_deltaT
                0264      &         *uIceNm1(i,j,bi,bj)
                0265           fyTmp  (i,j,bi,bj) = FORCEY0(i,j,bi,bj)
                0266      &         +seaiceMassV(i,j,bi,bj)*recip_deltaT
                0267      &         *vIceNm1(i,j,bi,bj)
                0268          ENDDO
                0269         ENDDO
                0270        ENDDO
                0271       ENDDO
0d75a51072 Mart*0272
143d9ce879 Dami*0273 #ifdef SEAICE_ALLOW_LSR_FLEX
                0274       IF ( SEAICEuseLSRflex ) THEN
                0275        residkm1       = 1. _d 0
                0276        residkm2       = 1. _d 0
                0277        residIniNonLin = 1. _d 0
                0278       ENDIF
                0279 #endif /* SEAICE_ALLOW_LSR_FLEX */
                0280
                0281       doNonLinLoop = .TRUE.
                0282
                0283 C     Start of non-linear iteration ====================================
                0284
                0285 C     Note: in S/R SEAICE_CHECK, we make sure that nonLinIterLoc
                0286 C     (=SEAICEnonLinIterMax) is not larger than MPSEUDOTIMESTEPS (which
                0287 C     defines the length of the TAF tape, so it is safe to use the
                0288 C     runtime parameter value here.
79df32c3f1 Mart*0289       DO ipass=1,nonLinIterLoc
9b4636ff01 Patr*0290
                0291 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0292        nlkey = (ikey_dynamics-1)*MPSEUDOTIMESTEPS + ipass
9b4636ff01 Patr*0293 # ifdef SEAICE_ALLOW_FREEDRIFT
edb6656069 Mart*0294 CADJ STORE uice_fd, vice_fd = comlev1_dynsol, key = nlkey, kind=isbyte
                0295 CADJ STORE ures,vres        = comlev1_dynsol, key = nlkey, kind=isbyte
9b4636ff01 Patr*0296 # endif
edb6656069 Mart*0297 CADJ STORE utmp,vtmp        = comlev1_dynsol, key = nlkey, kind=isbyte
9b4636ff01 Patr*0298 #endif /* ALLOW_AUTODIFF_TAMC */
                0299
143d9ce879 Dami*0300       IF ( doNonLinLoop ) THEN
bb3037b865 Gael*0301 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0302 CADJ STORE uice    = comlev1_dynsol, key = nlkey, kind=isbyte
                0303 CADJ STORE vice    = comlev1_dynsol, key = nlkey, kind=isbyte
                0304 CADJ STORE uicenm1 = comlev1_dynsol, key = nlkey, kind=isbyte
                0305 CADJ STORE vicenm1 = comlev1_dynsol, key = nlkey, kind=isbyte
bb3037b865 Gael*0306 #endif /* ALLOW_AUTODIFF_TAMC */
                0307
589eca759f Mart*0308 C SET SOME VALUES
e45202e340 Mart*0309       WFAU=SEAICE_LSRrelaxU
                0310       WFAV=SEAICE_LSRrelaxV
589eca759f Mart*0311       WFAU2=ZERO
                0312       WFAV2=ZERO
7b7a25aa58 Jean*0313       S1 = 0.
                0314       S2 = 0.
e45202e340 Mart*0315       S1A=0.80 _d 0
                0316       S2A=0.80 _d 0
589eca759f Mart*0317
4a08d54d3a Mart*0318       DO bj=myByLo(myThid),myByHi(myThid)
                0319        DO bi=myBxLo(myThid),myBxHi(myThid)
                0320         IF ( ipass .EQ. 1 ) THEN
                0321 C     This is the predictor time step
a10601a61f Mart*0322          DO j=1-OLy,sNy+OLy
                0323           DO i=1-OLx,sNx+OLx
4a08d54d3a Mart*0324            uIceC(i,j,bi,bj) = uIce(i,j,bi,bj)
                0325            vIceC(i,j,bi,bj) = vIce(i,j,bi,bj)
a10601a61f Mart*0326           ENDDO
                0327          ENDDO
4a08d54d3a Mart*0328         ELSEIF ( ipass .EQ. 2 .AND. SEAICEnonLinIterMax .LE. 2 ) THEN
                0329 C     This is the modified Euler step
a10601a61f Mart*0330          DO j=1-OLy,sNy+OLy
                0331           DO i=1-OLx,sNx+OLx
4a08d54d3a Mart*0332            uIce (i,j,bi,bj)=HALF*( uIce(i,j,bi,bj)+uIceNm1(i,j,bi,bj) )
                0333            vIce (i,j,bi,bj)=HALF*( vIce(i,j,bi,bj)+vIceNm1(i,j,bi,bj) )
                0334            uIceC(i,j,bi,bj)=uIce(i,j,bi,bj)
                0335            vIceC(i,j,bi,bj)=vIce(i,j,bi,bj)
a10601a61f Mart*0336           ENDDO
                0337          ENDDO
4a08d54d3a Mart*0338         ELSE
                0339 C     This is the case for SEAICEnonLinIterMax > 2, and here we use
                0340 C     a different iterative scheme. u/vIceC = u/vIce is unstable, and
                0341 C     different stabilisation methods are possible.
de36f9bd7d Mart*0342          DO j=1-OLy,sNy+OLy
                0343           DO i=1-OLx,sNx+OLx
4a08d54d3a Mart*0344 C     This is stable but slow to converge.
                0345 C          uIceC(i,j,bi,bj) = HALF*(uIce(i,j,bi,bj)+uIceNm1(i,j,bi,bj))
                0346 C          vIceC(i,j,bi,bj) = HALF*(vIce(i,j,bi,bj)+vIceNm1(i,j,bi,bj))
                0347 C     This converges slightly faster than the previous lines.
                0348            uIceC(i,j,bi,bj) = HALF*( uIce(i,j,bi,bj)+uIceC(i,j,bi,bj) )
                0349            vIceC(i,j,bi,bj) = HALF*( vIce(i,j,bi,bj)+vIceC(i,j,bi,bj) )
de36f9bd7d Mart*0350           ENDDO
                0351          ENDDO
4a08d54d3a Mart*0352         ENDIF
                0353 C     bi/bj-loops
de36f9bd7d Mart*0354        ENDDO
4a08d54d3a Mart*0355       ENDDO
a10601a61f Mart*0356
589eca759f Mart*0357 #ifdef ALLOW_AUTODIFF_TAMC
3e8bf7e743 Mart*0358 C     That is an important one! Note, that
143d9ce879 Dami*0359 C     * ipass index is part of nlkey because we are inside
                0360 C       the nonlinear iteration
3e8bf7e743 Mart*0361 C     * this storing is still outside the lsr iteration loop
edb6656069 Mart*0362 CADJ STORE uice  = comlev1_dynsol, kind=isbyte, key = nlkey
                0363 CADJ STORE vice  = comlev1_dynsol, kind=isbyte, key = nlkey
                0364 CADJ STORE uicec = comlev1_dynsol, kind=isbyte, key = nlkey
                0365 CADJ STORE vicec = comlev1_dynsol, kind=isbyte, key = nlkey
589eca759f Mart*0366 #endif /* ALLOW_AUTODIFF_TAMC */
                0367
                0368       CALL SEAICE_CALC_STRAINRATES(
                0369      I     uIceC, vIceC,
                0370      O     e11, e22, e12,
3e8bf7e743 Mart*0371      I     ipass, myTime, myIter, myThid )
589eca759f Mart*0372
                0373       CALL SEAICE_CALC_VISCOSITIES(
8e32c48b8f Mart*0374      I     e11, e22, e12, SEAICE_zMin, SEAICE_zMax, HEFFM, press0,
                0375      I     tensileStrFac,
0fd6b2de1a Mart*0376      O     eta, etaZ, zeta, zetaZ, press, deltaC,
3e8bf7e743 Mart*0377      I     ipass, myTime, myIter, myThid )
589eca759f Mart*0378
3f31c7d5de Mart*0379       CALL SEAICE_OCEANDRAG_COEFFS(
ec0d7df165 Mart*0380      I     uIceC, vIceC, HEFFM,
3f31c7d5de Mart*0381      O     DWATN,
3e8bf7e743 Mart*0382      I     ipass, myTime, myIter, myThid )
df1dac8b7b Mart*0383
                0384 #ifdef SEAICE_ALLOW_BOTTOMDRAG
d5254b4e3d Mart*0385       CALL SEAICE_BOTTOMDRAG_COEFFS(
ec0d7df165 Mart*0386      I     uIceC, vIceC, HEFFM,
df1dac8b7b Mart*0387 #ifdef SEAICE_ITD
                0388      I     HEFFITD, AREAITD, AREA,
                0389 #else
                0390      I     HEFF, AREA,
cbf501ab81 Jean*0391 #endif
df1dac8b7b Mart*0392      O     CbotC,
                0393      I     ipass, myTime, myIter, myThid )
                0394 #endif /* SEAICE_ALLOW_BOTTOMDRAG */
5bb179ddc2 Mart*0395
                0396 #ifdef SEAICE_ALLOW_SIDEDRAG
                0397       IF ( SEAICEsideDrag .NE. 0. _d 0 ) CALL SEAICE_SIDEDRAG_STRESS(
                0398      I     uIceC, vIceC, coastRoughU, coastRoughV, AREA,
                0399      O     sideDragU, sideDragV,
                0400      I     ipass, myTime, myIter, myThid )
                0401 #ifdef ALLOW_AUTODIFF_TAMC
                0402 CADJ STORE sideDragU = comlev1_dynsol, kind=isbyte, key = nlkey
                0403 CADJ STORE sideDragV = comlev1_dynsol, kind=isbyte, key = nlkey
                0404 #endif /* ALLOW_AUTODIFF_TAMC */
                0405 #endif /* SEAICE_ALLOW_SIDEDRAG */
d681af2948 Mart*0406 C
                0407 C     some abbreviations
                0408 C
                0409       DO bj=myByLo(myThid),myByHi(myThid)
                0410        DO bi=myBxLo(myThid),myBxHi(myThid)
4a08d54d3a Mart*0411         DO j=jMin-1,jMax
                0412          DO i=iMin-1,iMax
                0413           etaPlusZeta (i,j,bi,bj) = ETA (i,j,bi,bj)+ZETA(i,j,bi,bj)
                0414           zetaMinusEta(i,j,bi,bj) = ZETA(i,j,bi,bj)-ETA (i,j,bi,bj)
d681af2948 Mart*0415          ENDDO
                0416         ENDDO
4a08d54d3a Mart*0417         DO j=1-OLy,sNy+OLy
                0418          DO i=1-OLx,sNx+OLx
                0419           dragSym(i,j,bi,bj) = DWATN(i,j,bi,bj)*COSWAT
df1dac8b7b Mart*0420 #ifdef SEAICE_ALLOW_BOTTOMDRAG
4a08d54d3a Mart*0421      &         +CbotC(i,j,bi,bj)
df1dac8b7b Mart*0422 #endif /* SEAICE_ALLOW_BOTTOMDRAG */
                0423          ENDDO
                0424         ENDDO
d681af2948 Mart*0425        ENDDO
                0426       ENDDO
                0427 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0428 C
                0429 C     Set up of right hand sides of momentum equations starts here
                0430 C
589eca759f Mart*0431       DO bj=myByLo(myThid),myByHi(myThid)
                0432        DO bi=myBxLo(myThid),myBxHi(myThid)
4a08d54d3a Mart*0433         DO j=jMin,jMax
                0434          DO i=iMin,iMax
efcb5efb58 Jean*0435 C     set up anti symmetric drag force and add in current force
589eca759f Mart*0436 C     ( remember to average to correct velocity points )
4a08d54d3a Mart*0437           FORCEX(i,j,bi,bj) = fxTmp(i,j,bi,bj)+
                0438      &         ( 0.5 _d 0 * ( DWATN(i,j,bi,bj)+DWATN(i-1,j,bi,bj) ) *
                0439      &         COSWAT * uVel(i,j,kSrf,bi,bj)
                0440      &         - SIGN(SINWAT, _fCori(i,j,bi,bj))* 0.5 _d 0 *
                0441      &         ( DWATN(i  ,j,bi,bj) * 0.5 _d 0 *
                0442      &          (vVel(i  ,j  ,kSrf,bi,bj)-vIceC(i  ,j  ,bi,bj)
                0443      &          +vVel(i  ,j+1,kSrf,bi,bj)-vIceC(i  ,j+1,bi,bj))
                0444      &         + DWATN(i-1,j,bi,bj) * 0.5 _d 0 *
                0445      &          (vVel(i-1,j  ,kSrf,bi,bj)-vIceC(i-1,j  ,bi,bj)
                0446      &          +vVel(i-1,j+1,kSrf,bi,bj)-vIceC(i-1,j+1,bi,bj))
                0447      &         ) ) * areaW(i,j,bi,bj)
                0448           FORCEY(i,j,bi,bj) = fyTmp(i,j,bi,bj)+
                0449      &         ( 0.5 _d 0 * ( DWATN(i,j,bi,bj)+DWATN(i,j-1,bi,bj) ) *
                0450      &         COSWAT * vVel(i,j,kSrf,bi,bj)
                0451      &         + SIGN(SINWAT, _fCori(i,j,bi,bj)) * 0.5 _d 0 *
                0452      &         ( DWATN(i,j  ,bi,bj) * 0.5 _d 0 *
                0453      &          (uVel(i  ,j  ,kSrf,bi,bj)-uIceC(i  ,j  ,bi,bj)
                0454      &          +uVel(i+1,j  ,kSrf,bi,bj)-uIceC(i+1,j  ,bi,bj))
                0455      &         + DWATN(i,j-1,bi,bj) * 0.5 _d 0 *
                0456      &          (uVel(i  ,j-1,kSrf,bi,bj)-uIceC(i  ,j-1,bi,bj)
                0457      &          +uVel(i+1,j-1,kSrf,bi,bj)-uIceC(i+1,j-1,bi,bj))
                0458      &         ) ) * areaS(i,j,bi,bj)
589eca759f Mart*0459          ENDDO
                0460         ENDDO
5fb12930f9 Mart*0461         DO j=jMin,jMax
                0462          DO i=iMin,iMax
efcb5efb58 Jean*0463 C-    add Coriolis term
4a08d54d3a Mart*0464           FORCEX(i,j,bi,bj) = FORCEX(i,j,bi,bj) + HALF*
                0465      &         ( seaiceMassC(i  ,j,bi,bj) * _fCori(i  ,j,bi,bj)
efcb5efb58 Jean*0466      &          *0.5 _d 0*( vIceC( i ,j,bi,bj)+vIceC( i ,j+1,bi,bj) )
4a08d54d3a Mart*0467      &         + seaiceMassC(i-1,j,bi,bj) * _fCori(i-1,j,bi,bj)
efcb5efb58 Jean*0468      &          *0.5 _d 0*( vIceC(i-1,j,bi,bj)+vIceC(i-1,j+1,bi,bj) ) )
4a08d54d3a Mart*0469           FORCEY(i,j,bi,bj) = FORCEY(i,j,bi,bj) - HALF*
                0470      &         ( seaiceMassC(i,j  ,bi,bj) * _fCori(i,j  ,bi,bj)
efcb5efb58 Jean*0471      &         *0.5 _d 0*( uIceC(i  ,j  ,bi,bj)+uIceC(i+1,  j,bi,bj) )
4a08d54d3a Mart*0472      &         + seaiceMassC(i,j-1,bi,bj) * _fCori(i,j-1,bi,bj)
efcb5efb58 Jean*0473      &         *0.5 _d 0*( uIceC(i  ,j-1,bi,bj)+uIceC(i+1,j-1,bi,bj) ) )
                0474          ENDDO
                0475         ENDDO
4a08d54d3a Mart*0476 C      mask the forcing so far
5fb12930f9 Mart*0477         DO j=jMin,jMax
                0478          DO i=iMin,iMax
4a08d54d3a Mart*0479           FORCEX(i,j,bi,bj) = FORCEX(i,j,bi,bj)*seaiceMaskU(i,j,bi,bj)
                0480           FORCEY(i,j,bi,bj) = FORCEY(i,j,bi,bj)*seaiceMaskV(i,j,bi,bj)
589eca759f Mart*0481          ENDDO
                0482         ENDDO
7b7a25aa58 Jean*0483 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
b7c2bb63b7 Mart*0484 C
                0485 C     u-equation
                0486 C
bfe6f7447e Mart*0487         DO j=jMin,jMax
                0488          DO i=iMin,iMax
4a08d54d3a Mart*0489           rhsU(i,j,bi,bj) = FORCEX(i,j,bi,bj)
3ff8c9e192 Jean*0490          ENDDO
                0491         ENDDO
bfe6f7447e Mart*0492         CALL SEAICE_LSR_RHSU(
0fd6b2de1a Mart*0493      I       zetaMinusEta, etaPlusZeta, etaZ, zetaZ, press,
cbf501ab81 Jean*0494      I       uIceC, vIceC,
bfe6f7447e Mart*0495      U       rhsU,
                0496      I       iMin, iMax, jMin, jMax, bi, bj, myThid )
b7c2bb63b7 Mart*0497 C
                0498 C     v-equation
                0499 C
bfe6f7447e Mart*0500         DO j=jMin,jMax
                0501          DO i=iMin,iMax
4a08d54d3a Mart*0502           rhsV(i,j,bi,bj) = FORCEY(i,j,bi,bj)
3ff8c9e192 Jean*0503          ENDDO
                0504         ENDDO
bfe6f7447e Mart*0505         CALL SEAICE_LSR_RHSV(
0fd6b2de1a Mart*0506      I       zetaMinusEta, etaPlusZeta, etaZ, zetaZ, press,
cbf501ab81 Jean*0507      I       uIceC, vIceC,
bfe6f7447e Mart*0508      U       rhsV,
                0509      I       iMin, iMax, jMin, jMax, bi, bj, myThid )
8ce59fd29e Mart*0510 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0511 #ifdef SEAICE_ALLOW_MOM_ADVECTION
cbf501ab81 Jean*0512         IF ( SEAICEmomAdvection ) THEN
4a08d54d3a Mart*0513          DO j=1-OLy,sNy+OLy
                0514           DO i=1-OLx,sNx+OLx
                0515            gUmom(i,j) = 0. _d 0
                0516            gVmom(i,j) = 0. _d 0
8ce59fd29e Mart*0517           ENDDO
                0518          ENDDO
                0519          CALL SEAICE_MOM_ADVECTION(
                0520      I        bi,bj,iMin,iMax,jMin,jMax,
                0521      I        uIceC, vIceC,
                0522      O        gUmom, gVmom,
                0523      I        myTime, myIter, myThid )
4a08d54d3a Mart*0524          DO j=jMin,jMax
                0525           DO i=jMin,jMax
                0526            rhsU(i,j,bi,bj) = rhsU(i,j,bi,bj) + gUmom(i,j)
                0527            rhsV(i,j,bi,bj) = rhsV(i,j,bi,bj) + gVmom(i,j)
8ce59fd29e Mart*0528           ENDDO
                0529          ENDDO
                0530         ENDIF
                0531 #endif /* SEAICE_ALLOW_MOM_ADVECTION */
b7c2bb63b7 Mart*0532        ENDDO
                0533       ENDDO
d681af2948 Mart*0534 C
                0535 C     Set up of right hand sides of momentum equations ends here
                0536 C
                0537 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0538 C
                0539 C     calculate coefficients of tridiagonal matrices for both u- and
                0540 C     v-equations
                0541 C
cbf501ab81 Jean*0542       CALL SEAICE_LSR_CALC_COEFFS(
df1dac8b7b Mart*0543      I     etaPlusZeta, zetaMinusEta, etaZ, zetaZ, dragSym,
d681af2948 Mart*0544      O     AU, BU, CU, AV, BV, CV, uRt1, uRt2, vRt1, vRt2,
                0545      I     iMin, iMax, jMin, jMax, myTime, myIter, myThid )
589eca759f Mart*0546
5bb179ddc2 Mart*0547 #ifdef SEAICE_ALLOW_SIDEDRAG
                0548       IF ( SEAICEsideDrag .NE. 0. _d 0 ) THEN
                0549        DO bj=myByLo(myThid),myByHi(myThid)
                0550         DO bi=myBxLo(myThid),myBxHi(myThid)
                0551          DO j=jMin,jMax
                0552           DO i=iMin,iMax
                0553            BU(i,j,bi,bj) = BU(i,j,bi,bj)
                0554      &          + seaiceMaskU(i,j,bi,bj) * sideDragU(i,j,bi,bj)
                0555            BV(i,j,bi,bj) = BV(i,j,bi,bj)
                0556      &          + seaiceMaskV(i,j,bi,bj) * sideDragV(i,j,bi,bj)
                0557           ENDDO
                0558          ENDDO
                0559         ENDDO
                0560        ENDDO
                0561       ENDIF
                0562 #endif /* SEAICE_ALLOW_SIDEDRAG */
                0563
fec75090d3 Jean*0564 #ifndef OBCS_UVICE_OLD
b7c2bb63b7 Mart*0565 C--     prevent tri-diagonal solver from modifying OB values:
                0566       DO bj=myByLo(myThid),myByHi(myThid)
                0567        DO bi=myBxLo(myThid),myBxHi(myThid)
4a08d54d3a Mart*0568         DO j=jMin,jMax
                0569          DO i=iMin,iMax
b7c2bb63b7 Mart*0570           IF ( maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj) .EQ. 0. ) THEN
4a08d54d3a Mart*0571            AU(i,j,bi,bj)   = ZERO
                0572            BU(i,j,bi,bj)   = ONE
                0573            CU(i,j,bi,bj)   = ZERO
                0574            uRt1(i,j,bi,bj) = ZERO
                0575            uRt2(i,j,bi,bj) = ZERO
                0576            rhsU(i,j,bi,bj) = uIce(i,j,bi,bj)
b7c2bb63b7 Mart*0577           ENDIF
fec75090d3 Jean*0578           IF ( maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj) .EQ. 0. ) THEN
4a08d54d3a Mart*0579            AV(i,j,bi,bj)   = ZERO
                0580            BV(i,j,bi,bj)   = ONE
                0581            CV(i,j,bi,bj)   = ZERO
                0582            vRt1(i,j,bi,bj) = ZERO
                0583            vRt2(i,j,bi,bj) = ZERO
                0584            rhsV(i,j,bi,bj) = vIce(i,j,bi,bj)
fec75090d3 Jean*0585           ENDIF
                0586          ENDDO
                0587         ENDDO
589eca759f Mart*0588        ENDDO
                0589       ENDDO
b7c2bb63b7 Mart*0590 #endif /* OBCS_UVICE_OLD */
589eca759f Mart*0591
7b7a25aa58 Jean*0592 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0593
589eca759f Mart*0594 #ifdef ALLOW_DEBUG
86fa82a64d Jean*0595       IF ( debugLevel .GE. debLevD ) THEN
                0596         WRITE(msgBuf,'(A,I3,A)')
3e8bf7e743 Mart*0597      &        'Uice pre  iter (SEAICE_LSR', MOD(ipass,1000), ')'
7b7a25aa58 Jean*0598         CALL DEBUG_STATS_RL( 1, UICE, msgBuf, myThid )
                0599         WRITE(msgBuf,'(A,I3,A)')
3e8bf7e743 Mart*0600      &        'Vice pre  iter (SEAICE_LSR', MOD(ipass,1000), ')'
86fa82a64d Jean*0601         CALL DEBUG_STATS_RL( 1, VICE, msgBuf, myThid )
589eca759f Mart*0602       ENDIF
                0603 #endif /* ALLOW_DEBUG */
                0604
efcb5efb58 Jean*0605 C--   Calculate initial residual of the linearised system
143d9ce879 Dami*0606       IF ( printResidual .OR. LSR_mixIniGuess.GE.1
                0607      &     .OR. SEAICEuseLSRflex )
                0608      &    CALL SEAICE_RESIDUAL(
d681af2948 Mart*0609      I                  rhsU, rhsV, uRt1, uRt2, vRt1, vRt2,
                0610      I                  AU, BU, CU, AV, BV, CV, uIce, vIce,
                0611      O                  residUini, residVini, uTmp, vTmp,
                0612      I                  printResidual, myIter, myThid )
                0613
efcb5efb58 Jean*0614 #ifdef SEAICE_ALLOW_FREEDRIFT
d681af2948 Mart*0615       IF ( printResidual .OR. LSR_mixIniGuess.GE.1 ) THEN
                0616        IF ( LSR_mixIniGuess.GE.1 ) THEN
                0617         DO bj=myByLo(myThid),myByHi(myThid)
                0618          DO bi=myBxLo(myThid),myBxHi(myThid)
                0619           DO j=jMin,jMax
                0620            DO i=iMin,iMax
                0621             uIce_fd(i,j,bi,bj) = FORCEX(i,j,bi,bj)
                0622      &         / ( 1. _d 0 - seaiceMaskU(i,j,bi,bj)
e232688b80 Mart*0623      &           + seaiceMassU(i,j,bi,bj)*recip_deltaT
df1dac8b7b Mart*0624      &           + HALF*( dragSym(i,j,bi,bj) + dragSym(i-1,j,bi,bj) )
0c560a7486 Mart*0625      &                 * MAX(areaW(i,j,bi,bj),areaMinLoc)
d681af2948 Mart*0626      &           )
                0627             vIce_fd(i,j,bi,bj) = FORCEY(i,j,bi,bj)
                0628      &         / ( 1. _d 0 - seaiceMaskV(i,j,bi,bj)
e232688b80 Mart*0629      &           + seaiceMassV(i,j,bi,bj)*recip_deltaT
df1dac8b7b Mart*0630      &           + HALF*( dragSym(i,j,bi,bj) + dragSym(i,j-1,bi,bj) )
0c560a7486 Mart*0631      &                 * MAX(areaS(i,j,bi,bj),areaMinLoc)
d681af2948 Mart*0632      &           )
                0633            ENDDO
                0634           ENDDO
                0635          ENDDO
                0636         ENDDO
                0637         CALL EXCH_UV_XY_RL( uIce_fd, vIce_fd, .TRUE., myThid )
                0638        ENDIF
                0639        IF ( LSR_mixIniGuess.GE.0 )
efcb5efb58 Jean*0640      &  CALL SEAICE_RESIDUAL(
                0641      I                  rhsU, rhsV, uRt1, uRt2, vRt1, vRt2,
                0642      I                  AU, BU, CU, AV, BV, CV, uIce_fd, vIce_fd,
                0643      O                  residU_fd, residV_fd, uRes, vRes,
                0644      I                  printResidual, myIter, myThid )
                0645       ENDIF
                0646       IF ( LSR_mixIniGuess.GE.2 ) THEN
9b4636ff01 Patr*0647 # ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0648 CADJ STORE uice, vice = comlev1_dynsol, key = nlkey, kind=isbyte
                0649 CADJ STORE utmp, vtmp = comlev1_dynsol, key = nlkey, kind=isbyte
9b4636ff01 Patr*0650 # endif /* ALLOW_AUTODIFF_TAMC */
efcb5efb58 Jean*0651        DO bj=myByLo(myThid),myByHi(myThid)
                0652         DO bi=myBxLo(myThid),myBxHi(myThid)
d681af2948 Mart*0653          DO j=jMin,jMax
                0654           DO i=iMin,iMax
                0655            errIni = uTmp(i,j,bi,bj)*uTmp(i,j,bi,bj)
                0656            errFD  = uRes(i,j,bi,bj)*uRes(i,j,bi,bj)
                0657            IF ( LSR_mixIniGuess.GE.4 ) THEN
                0658             errIni = errIni*errIni
                0659             errFD  = errFD *errFD
                0660            ENDIF
                0661            errSum = ( errIni + errFD )
                0662      &          *maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj)
                0663            IF ( errSum.GT.0. ) THEN
                0664             uIce(i,j,bi,bj) = ( errFD *uIce(i,j,bi,bj)
                0665      &                        + errIni*uIce_fd(i,j,bi,bj)
                0666      &                        )/errSum
                0667            ENDIF
                0668            errIni = vTmp(i,j,bi,bj)*vTmp(i,j,bi,bj)
                0669            errFD  = vRes(i,j,bi,bj)*vRes(i,j,bi,bj)
                0670            IF ( LSR_mixIniGuess.GE.4 ) THEN
                0671             errIni = errIni*errIni
                0672             errFD  = errFD *errFD
                0673            ENDIF
                0674            errSum = ( errIni + errFD )
                0675      &          *maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj)
                0676            IF ( errSum.GT.0. ) THEN
                0677             vIce(i,j,bi,bj) = ( errFD *vIce(i,j,bi,bj)
                0678      &                        + errIni*vIce_fd(i,j,bi,bj)
                0679      &                        )/errSum
                0680            ENDIF
efcb5efb58 Jean*0681           ENDDO
d681af2948 Mart*0682          ENDDO
efcb5efb58 Jean*0683         ENDDO
                0684        ENDDO
                0685        CALL EXCH_UV_XY_RL( uIce, vIce, .TRUE., myThid )
dfc97bd020 Mart*0686 #ifndef ALLOW_AUTODIFF_TAMC
0e2f719d67 Mart*0687 C     Here residU/Vmix and u/vTmp are computed but never used for anything
                0688 C     but reporting to STDOUT. Still TAF thinks this requires recomputations
                0689 C     so we hide this part of the code from TAF.
efcb5efb58 Jean*0690        IF ( printResidual ) THEN
                0691         CALL SEAICE_RESIDUAL(
                0692      I                  rhsU, rhsV, uRt1, uRt2, vRt1, vRt2,
                0693      I                  AU, BU, CU, AV, BV, CV, uIce, vIce,
                0694      O                  residUmix, residVmix, uTmp, vTmp,
                0695      I                  printResidual, myIter, myThid )
                0696        ENDIF
dfc97bd020 Mart*0697 #endif /* ALLOW_AUTODIFF_TAMC */
efcb5efb58 Jean*0698       ENDIF
                0699 #endif /* SEAICE_ALLOW_FREEDRIFT */
                0700
                0701 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0702
589eca759f Mart*0703 C NOW DO ITERATION
                0704
7b7a25aa58 Jean*0705       doIterate4u = .TRUE.
                0706       doIterate4v = .TRUE.
589eca759f Mart*0707
                0708 C ITERATION START -----------------------------------------------------
d681af2948 Mart*0709 C     Ideally, we would like to use the loop-directive for this
cbf501ab81 Jean*0710 C     iterative loop, but (unsurprisingly) it does not seem to work for
964a494b4b Mart*0711 C     this convoluted case of mixed iterations, so it is commented out
                0712 CML#ifdef ALLOW_AUTODIFF_TAMC
0e2f719d67 Mart*0713 CMLCADJ LOOP = iteration uice, vice = comlev1_lsr
964a494b4b Mart*0714 CML#endif /* ALLOW_AUTODIFF_TAMC */
9b4636ff01 Patr*0715
0d75a51072 Mart*0716       ICOUNT1 = linearIterLoc
                0717       ICOUNT2 = linearIterLoc
9b4636ff01 Patr*0718
143d9ce879 Dami*0719 #ifdef SEAICE_ALLOW_LSR_FLEX
                0720 C     first IF doNonLinLoop
                0721       ENDIF
                0722
                0723       IF ( SEAICEuseLSRflex ) THEN
                0724 C     Do not do anything if the initial residuals are zero
                0725        IF ( residUini .EQ. 0. _d 0) THEN
                0726         doIterate4u = .FALSE.
                0727         ICOUNT1 = 0
                0728         S1      = residUini
                0729        ENDIF
                0730        IF ( residVini .EQ. 0. _d 0) THEN
                0731         doIterate4v = .FALSE.
                0732         ICOUNT2 = 0
                0733         S2      = residVini
                0734        ENDIF
                0735
                0736        residIni = SQRT(residUini*residUini+residVini*residVini)
                0737
                0738 C     If you want to use SEAICEnonLinTol to be an absolute criterion,
                0739 C     comment out the following line.
                0740        IF ( ipass .EQ. 1 ) residIniNonLin = residIni
                0741
                0742 C     Skip the the rest of the nonlinear loop if there is nothing to do.
                0743        doNonLinLoop = doNonLinLoop .AND.( doIterate4u .OR. doIterate4v )
                0744 C     Check convergence of nonlinear loop, but do at least two passes
                0745        doNonLinLoop = doNonLinLoop .AND. .NOT.
                0746      &      ( ipass .GT. 2 .AND.
                0747      &      ( residIni .LT. SEAICEnonLinTol*residIniNonlin ) )
                0748
                0749 C     Catch the case of zero residual
                0750        IF  ( residIni .EQ. 0. _d 0 ) residIni = 1. _d -20
                0751 C     From Lemieux et al. 2010
                0752 C       LSRflexFac = MIN(LSRflexFac/residkm2*0.9,LSRflexFac)
                0753 C     Other possibilities
                0754        LSRflexFac = 1. _d 0 / (1. _d 0 + ABS( LOG10(residIni) ) )
                0755 C       LSRflexFac = ABS(LOG10(residkm2))/ABS(LOG10(residkm1))
                0756 C       LSRflexFac = (residkm1/residkm2)/(1+ABS(LOG10(residIni)))
                0757        LSRflexFac = MIN(LSRflexFac, 0.99 _d 0 )
                0758 C      LSRflexFac = MAX(LSRflexFac, 0.3 _d 0 )
                0759
                0760 C       IF ( ipass .LE. 3 ) THEN
                0761 C        LSRflexFac = 0.99 _d 0
                0762 C       ENDIF
                0763
                0764        LSR_ERRORfu = residUini*LSRflexFac
                0765        LSR_ERRORfv = residVini*LSRflexFac
                0766
                0767        IF ( printResidual ) THEN
                0768         _BEGIN_MASTER( myThid )
                0769         WRITE(standardMessageUnit,'(A,1X,I5,1X,F10.6,3(1X,E12.6))')
                0770      &       ' SEAICE_LSR: ipass, FLEX_FACTOR,'//
                0771      &       ' LSR_ERRORfu, LSR_ERRORfv, residIni =',
                0772      &       ipass, LSRflexFac, LSR_ERRORfu, LSR_ERRORfv, residIni
                0773         _END_MASTER( myThid )
                0774        ENDIF
                0775
                0776       ENDIF
                0777
                0778 C     test doNonLinLoop a second time, as it may have changed
                0779       IF ( doNonLinLoop ) THEN
                0780 #endif /* SEAICE_ALLOW_LSR_FLEX */
                0781
0d75a51072 Mart*0782       DO m = 1, linearIterLoc
9b4636ff01 Patr*0783
0e2f719d67 Mart*0784 #ifdef ALLOW_AUTODIFF_TAMC
                0785 # ifdef SEAICE_LSR_ADJOINT_ITER
edb6656069 Mart*0786        lnkey = (nlkey-1)*SOLV_MAX_FIXED + MIN(m,SOLV_MAX_FIXED)
0e2f719d67 Mart*0787 # else
edb6656069 Mart*0788        lnkey = nlkey
0e2f719d67 Mart*0789 # endif /* SEAICE_LSR_ADJOINT_ITER */
3e8bf7e743 Mart*0790 C     Storing these scalars and logicals is necessary to do an exact
                0791 C     backward iteration
edb6656069 Mart*0792 CADJ STORE wfau, wfav = comlev1_lsr,  kind=isbyte, key = lnkey
                0793 CADJ STORE doIterate4u, doIterate4v = comlev1_lsr, key = lnkey
0e2f719d67 Mart*0794 #endif /* ALLOW_AUTODIFF_TAMC */
9b4636ff01 Patr*0795
7b7a25aa58 Jean*0796        IF ( doIterate4u .OR. doIterate4v ) THEN
589eca759f Mart*0797
7b7a25aa58 Jean*0798         IF ( useCubedSphereExchange ) THEN
                0799           doIterate4u = .TRUE.
                0800           doIterate4v = .TRUE.
                0801         ENDIF
589eca759f Mart*0802
9b4636ff01 Patr*0803 # ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0804 C     Note that lnkey can get very large when SEAICE_LSR_ADJOINT_ITER
3e8bf7e743 Mart*0805 C     is defined for an accurate adjoint, otherwise it is the same for
                0806 C     each m, so that the tapes get overwritten for each iterate and
                0807 C     the adjoint is necessarily wrong.
edb6656069 Mart*0808 CADJ STORE uice, vice = comlev1_lsr, kind=isbyte, key = lnkey
9b4636ff01 Patr*0809 # endif /* ALLOW_AUTODIFF_TAMC */
                0810
2d6335c1f7 Jean*0811 #ifdef SEAICE_GLOBAL_3DIAG_SOLVER
                0812         IF ( SEAICEuseMultiTileSolver ) THEN
                0813
                0814          DO bj=myByLo(myThid),myByHi(myThid)
                0815           DO bi=myBxLo(myThid),myBxHi(myThid)
                0816            IF ( doIterate4u ) THEN
5fb12930f9 Mart*0817             DO j=jMin,jMax
                0818              DO i=iMin,iMax
2d6335c1f7 Jean*0819               uTmp(i,j,bi,bj)  = uIce(i,j,bi,bj)
                0820               rhsUx(i,j,bi,bj) = ( rhsU(i,j,bi,bj)
                0821      &             + uRt1(i,j,bi,bj)*uIce(i,j-1,bi,bj)
                0822      &             + uRt2(i,j,bi,bj)*uIce(i,j+1,bi,bj)
                0823      &                           )*seaiceMaskU(i,j,bi,bj)
                0824              ENDDO
                0825             ENDDO
                0826            ENDIF
                0827            IF ( doIterate4v ) THEN
5fb12930f9 Mart*0828             DO j=jMin,jMax
                0829              DO i=iMin,iMax
2d6335c1f7 Jean*0830               vTmp(i,j,bi,bj)  = vIce(i,j,bi,bj)
                0831               rhsVy(i,j,bi,bj) = ( rhsV(i,j,bi,bj)
4a08d54d3a Mart*0832      &             + vRt1(i,j,bi,bj)*vIce(i-1,j,bi,bj)
                0833      &             + vRt2(i,j,bi,bj)*vIce(i+1,j,bi,bj)
2d6335c1f7 Jean*0834      &                           )*seaiceMaskU(i,j,bi,bj)
                0835              ENDDO
                0836             ENDDO
                0837            ENDIF
                0838 C     end bi,bj-loops
                0839           ENDDO
                0840          ENDDO
                0841
0d75a51072 Mart*0842          iSubIter = linearIterLoc*(ipass-1) + m
2d6335c1f7 Jean*0843          CALL SOLVE_UV_TRIDIAGO(
                0844      I                 1, OLx, doIterate4u, doIterate4v,
                0845      I                 AU, BU, CU, rhsUx,
                0846      I                 AV, BV, CV, rhsVy,
                0847      O                 uIce, vIce, errCode,
                0848      I                 iSubIter, myIter, myThid )
                0849
                0850          DO bj=myByLo(myThid),myByHi(myThid)
                0851           DO bi=myBxLo(myThid),myBxHi(myThid)
                0852            IF ( doIterate4u ) THEN
5fb12930f9 Mart*0853             DO j=jMin,jMax
                0854              DO i=iMin,iMax
2d6335c1f7 Jean*0855               uIce(i,j,bi,bj) = uTmp(i,j,bi,bj)
                0856      &            + WFAU*( uIce(i,j,bi,bj)-uTmp(i,j,bi,bj) )
                0857              ENDDO
                0858             ENDDO
                0859            ENDIF
                0860            IF ( doIterate4v ) THEN
5fb12930f9 Mart*0861             DO j=jMin,jMax
                0862              DO i=iMin,iMax
2d6335c1f7 Jean*0863               vIce(i,j,bi,bj) = vTmp(i,j,bi,bj)
                0864      &            + WFAV*( vIce(i,j,bi,bj)-vTmp(i,j,bi,bj) )
                0865              ENDDO
                0866             ENDDO
                0867            ENDIF
                0868 C     end bi,bj-loops
                0869           ENDDO
                0870          ENDDO
                0871
                0872         ELSE
                0873 #endif /* SEAICE_GLOBAL_3DIAG_SOLVER */
                0874
7b7a25aa58 Jean*0875         DO bj=myByLo(myThid),myByHi(myThid)
                0876          DO bi=myBxLo(myThid),myBxHi(myThid)
                0877
f4df908d9c Mart*0878 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0879           tkey = bi + (bj-1)*nSx + (lnkey-1)*nSx*nSy
                0880 CADJ STORE uice(:,:,bi,bj) = comlev1_bibj_lsr,kind=isbyte,key=tkey
                0881 CADJ STORE vice(:,:,bi,bj) = comlev1_bibj_lsr,kind=isbyte,key=tkey
f4df908d9c Mart*0882 #endif /* ALLOW_AUTODIFF_TAMC */
4a08d54d3a Mart*0883 C-jmc: get fewer TAF warnings when always (no if doIterate) saving uIce,vIce:
3e8bf7e743 Mart*0884 C     save u/vIce prior to iteration
8df0ccd026 Jean*0885           DO j=1-OLy,sNy+OLy
                0886            DO i=1-OLx,sNx+OLx
4a08d54d3a Mart*0887             uTmp(i,j,bi,bj)=uIce(i,j,bi,bj)
                0888             vTmp(i,j,bi,bj)=vIce(i,j,bi,bj)
8df0ccd026 Jean*0889            ENDDO
                0890           ENDDO
                0891
7b7a25aa58 Jean*0892           IF ( doIterate4u ) THEN
                0893 C Solve for uIce :
5fb12930f9 Mart*0894            CALL SEAICE_LSR_TRIDIAGU(
8df0ccd026 Jean*0895      I          AU, BU, CU, uRt1, uRt2, rhsU, uTmp, seaiceMaskU, WFAU,
5fb12930f9 Mart*0896      U          uIce,
8d9be9cde5 Mart*0897      I          iMin, iMax, jMin, jMax, bi, bj, myTime, myIter, myThid )
3ff8c9e192 Jean*0898           ENDIF
589eca759f Mart*0899
7b7a25aa58 Jean*0900           IF ( doIterate4v ) THEN
                0901 C Solve for vIce
5fb12930f9 Mart*0902            CALL SEAICE_LSR_TRIDIAGV(
cbf501ab81 Jean*0903      I          AV, BV, CV, vRt1, vRt2, rhsV, vTmp, seaiceMaskV, WFAV,
5fb12930f9 Mart*0904      U          vIce,
8d9be9cde5 Mart*0905      I          iMin, iMax, jMin, jMax, bi, bj, myTime, myIter, myThid )
7b7a25aa58 Jean*0906           ENDIF
589eca759f Mart*0907
7b7a25aa58 Jean*0908 C     end bi,bj-loops
3ff8c9e192 Jean*0909          ENDDO
589eca759f Mart*0910         ENDDO
                0911
2d6335c1f7 Jean*0912 #ifdef SEAICE_GLOBAL_3DIAG_SOLVER
                0913         ENDIF
                0914 #endif /* SEAICE_GLOBAL_3DIAG_SOLVER */
                0915
143d9ce879 Dami*0916 #ifdef SEAICE_ALLOW_LSR_FLEX
                0917         IF ( SEAICEuseLSRflex ) THEN
                0918          IF ( MOD(m,SOLV_NCHECK) .EQ. 0
                0919      &        .AND. (doIterate4u.OR.doIterate4v) ) THEN
                0920           CALL SEAICE_RESIDUAL(
                0921      I         rhsU, rhsV, uRt1, uRt2, vRt1, vRt2,
                0922      I         AU, BU, CU, AV, BV, CV, uIce, vIce,
                0923      O         S1, S2, uTmp, vTmp,
                0924      I         printResidual, myIter, myThid )
                0925           IF( S1.LT.LSR_ERRORfu ) THEN
                0926            ICOUNT1=m
                0927            doIterate4u = .FALSE.
                0928           ENDIF
                0929           IF( S2.LT.LSR_ERRORfv ) THEN
                0930            ICOUNT2=m
                0931            doIterate4v = .FALSE.
                0932           ENDIF
                0933          ENDIF
                0934         ELSE
                0935 #endif
7b7a25aa58 Jean*0936         IF ( doIterate4u.AND.MOD(m,SOLV_NCHECK).EQ.0) THEN
                0937          S1=ZERO
                0938          DO bj=myByLo(myThid),myByHi(myThid)
                0939           DO bi=myBxLo(myThid),myBxHi(myThid)
4a08d54d3a Mart*0940            DO j=1,sNy
                0941             DO i=1,sNx
                0942              UERR=(uIce(i,j,bi,bj)-uTmp(i,j,bi,bj))
                0943      &               * seaiceMaskU(i,j,bi,bj)
7b7a25aa58 Jean*0944              S1=MAX(ABS(UERR),S1)
                0945             ENDDO
                0946            ENDDO
3ff8c9e192 Jean*0947           ENDDO
                0948          ENDDO
7b7a25aa58 Jean*0949          _GLOBAL_MAX_RL( S1, myThid )
                0950 c        WRITE(standardMessageUnit,'(A,2I6,1P4E16.9)')
3e8bf7e743 Mart*0951 c    &   ' U iters,error,WF = ',ipass,M,S1,S1A,WFAU
589eca759f Mart*0952 C SAFEGUARD AGAINST BAD FORCING ETC
7b7a25aa58 Jean*0953          IF(m.GT.1.AND.S1.GT.S1A) WFAU=WFAU2
                0954          S1A=S1
                0955          IF(S1.LT.LSR_ERROR) THEN
                0956           ICOUNT1=m
                0957           doIterate4u = .FALSE.
                0958          ENDIF
                0959         ENDIF
589eca759f Mart*0960
7b7a25aa58 Jean*0961         IF ( doIterate4v.AND.MOD(m,SOLV_NCHECK).EQ.0) THEN
                0962          S2=ZERO
                0963          DO bj=myByLo(myThid),myByHi(myThid)
                0964           DO bi=myBxLo(myThid),myBxHi(myThid)
4a08d54d3a Mart*0965            DO j=1,sNy
                0966             DO i=1,sNx
                0967              UERR=(vIce(i,j,bi,bj)-vTmp(i,j,bi,bj))
                0968      &               * seaiceMaskV(i,j,bi,bj)
7b7a25aa58 Jean*0969              S2=MAX(ABS(UERR),S2)
                0970             ENDDO
                0971            ENDDO
                0972           ENDDO
                0973          ENDDO
                0974          _GLOBAL_MAX_RL( S2, myThid )
                0975 C SAFEGUARD AGAINST BAD FORCING ETC
                0976          IF(m.GT.1.AND.S2.GT.S2A) WFAV=WFAV2
                0977          S2A=S2
                0978          IF(S2.LT.LSR_ERROR) THEN
                0979           ICOUNT2=m
                0980           doIterate4v = .FALSE.
                0981          ENDIF
                0982         ENDIF
143d9ce879 Dami*0983 #ifdef SEAICE_ALLOW_LSR_FLEX
                0984         ENDIF
                0985 #endif
7b7a25aa58 Jean*0986
                0987         CALL EXCH_UV_XY_RL( uIce, vIce,.TRUE.,myThid)
                0988
                0989 C--    end doIterate4u or doIterate4v
                0990        ENDIF
589eca759f Mart*0991       ENDDO
                0992 C ITERATION END -----------------------------------------------------
                0993
143d9ce879 Dami*0994 #ifdef SEAICE_ALLOW_LSR_FLEX
                0995       IF ( SEAICEuseLSRflex ) THEN
                0996        residkm2 = residkm1
                0997        residkm1 = residIni
                0998       ENDIF
                0999 #endif /* SEAICE_ALLOW_LSR_FLEX */
                1000
efcb5efb58 Jean*1001 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                1002
dfc97bd020 Mart*1003 #ifndef ALLOW_AUTODIFF_TAMC
0e2f719d67 Mart*1004 C     Here residU/Vend and u/vTmp are computed but never used for anything
                1005 C     but reporting to STDOUT. Still TAF thinks this requires recomputations
                1006 C     so we hide this part of the code from TAF.
efcb5efb58 Jean*1007       IF ( printResidual ) THEN
                1008 C--   Calculate final residual of the linearised system
                1009         CALL SEAICE_RESIDUAL(
                1010      I                  rhsU, rhsV, uRt1, uRt2, vRt1, vRt2,
                1011      I                  AU, BU, CU, AV, BV, CV, uIce, vIce,
                1012      O                  residUend, residVend, uTmp, vTmp,
                1013      I                  printResidual, myIter, myThid )
143d9ce879 Dami*1014
86fa82a64d Jean*1015         _BEGIN_MASTER( myThid )
c512e371cc drin*1016         WRITE(standardMessageUnit,'(A,1X,I9,1P3E16.8)')
e116529eb1 Mart*1017      &    ' SEAICE_LSR: Residual Initial ipass,Uice,Vice=',
3e8bf7e743 Mart*1018      &        ipass, residUini, residVini
efcb5efb58 Jean*1019 #ifdef SEAICE_ALLOW_FREEDRIFT
                1020         IF ( LSR_mixIniGuess.GE.0 )
                1021      &  WRITE(standardMessageUnit,'(A,1P3E16.8)')
                1022      &    ' SEAICE_LSR: Residual FrDrift U_fd,V_fd=',
                1023      &        residU_fd, residV_fd
                1024         IF ( LSR_mixIniGuess.GE.2 )
                1025      &  WRITE(standardMessageUnit,'(A,1P3E16.8)')
                1026      &    ' SEAICE_LSR: Residual Mix-ini Uice,Vice=',
                1027      &        residUmix, residVmix
                1028 #endif /* SEAICE_ALLOW_FREEDRIFT */
c512e371cc drin*1029         WRITE(standardMessageUnit,'(A,I9,A,I9,1P2E16.8)')
3e8bf7e743 Mart*1030      &    ' SEAICE_LSR (ipass=',ipass,') iters,dU,Resid=',
efcb5efb58 Jean*1031      &      ICOUNT1, S1, residUend
c512e371cc drin*1032         WRITE(standardMessageUnit,'(A,I9,A,I9,1P2E16.8)')
3e8bf7e743 Mart*1033      &    ' SEAICE_LSR (ipass=',ipass,') iters,dV,Resid=',
efcb5efb58 Jean*1034      &      ICOUNT2, S2, residVend
86fa82a64d Jean*1035         _END_MASTER( myThid )
7b7a25aa58 Jean*1036       ENDIF
efcb5efb58 Jean*1037 #ifdef ALLOW_DEBUG
7b7a25aa58 Jean*1038       IF ( debugLevel .GE. debLevD ) THEN
                1039         WRITE(msgBuf,'(A,I3,A)')
3e8bf7e743 Mart*1040      &        'Uice post iter (SEAICE_LSR', MOD(ipass,1000), ')'
7b7a25aa58 Jean*1041         CALL DEBUG_STATS_RL( 1, UICE, msgBuf, myThid )
                1042         WRITE(msgBuf,'(A,I3,A)')
3e8bf7e743 Mart*1043      &        'Vice post iter (SEAICE_LSR', MOD(ipass,1000), ')'
7b7a25aa58 Jean*1044         CALL DEBUG_STATS_RL( 1, VICE, msgBuf, myThid )
589eca759f Mart*1045       ENDIF
                1046 #endif /* ALLOW_DEBUG */
bc6ed4e27a Jean*1047       IF ( doIterate4u .OR. doIterate4v ) THEN
                1048         WRITE(msgBuf,'(2A,I10,A,I4,A)') '** WARNING ** SEAICE_LSR ',
3e8bf7e743 Mart*1049      &    '(it=', myIter, ',', ipass,') did not converge :'
bc6ed4e27a Jean*1050         CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
                1051      &                      SQUEEZE_RIGHT, myThid )
                1052         WRITE(msgBuf,'(2(A,I6,0PF6.3,1PE14.6))')
                1053      &      ' nIt,wFU,dU=', ICOUNT1, WFAU, S1,
                1054      &    ' ; nIt,wFV,dV=', ICOUNT2, WFAV, S2
                1055         CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
                1056      &                      SQUEEZE_RIGHT, myThid )
                1057       ENDIF
dfc97bd020 Mart*1058 #endif /* ALLOW_AUTODIFF_TAMC */
589eca759f Mart*1059
                1060 C     APPLY MASKS
                1061       DO bj=myByLo(myThid),myByHi(myThid)
                1062        DO bi=myBxLo(myThid),myBxHi(myThid)
4a08d54d3a Mart*1063         DO j=1-OLy,sNy+OLy
                1064          DO i=1-OLx,sNx+OLx
                1065           uIce(i,j,bi,bj)=uIce(i,j,bi,bj)* seaiceMaskU(i,j,bi,bj)
                1066           vIce(i,j,bi,bj)=vIce(i,j,bi,bj)* seaiceMaskV(i,j,bi,bj)
3ff8c9e192 Jean*1067          ENDDO
                1068         ENDDO
589eca759f Mart*1069        ENDDO
                1070       ENDDO
772590b63c Mart*1071 CML      CALL EXCH_UV_XY_RL( uIce, vIce,.TRUE.,myThid)
589eca759f Mart*1072
992ae64189 Mart*1073 #ifdef SEAICE_ALLOW_CHECK_LSR_CONVERGENCE
86fa82a64d Jean*1074       IF ( debugLevel .GE. debLevC ) THEN
992ae64189 Mart*1075        resnorm = 0. _d 0
                1076        EKnorm  = 0. _d 0
                1077        counter = 0. _d 0
                1078        DO bj=myByLo(myThid),myByHi(myThid)
                1079         DO bi=myBxLo(myThid),myBxHi(myThid)
                1080 C--   Compute residual norm and print
                1081          DO j=1,sNy
                1082           DO i=1,sNx
                1083            resnorm = resnorm + 0.5 _d 0 *
                1084      &          ( ( (uIceNm1(i,j,bi,bj)+uIceNm1(i+1,j,bi,bj))
                1085      &            - (uIce(i,j,bi,bj)+uIce(i+1,j,bi,bj)) )**2
                1086      &          + ( (vIceNm1(i,j,bi,bj)+vIceNm1(i,j+1,bi,bj))
                1087      &            - (vIce(i,j,bi,bj)+vIce(i,j+1,bi,bj)) )**2 )
8377b8ee87 Mart*1088            IF ( area(i,j,bi,bj) .GT. 0.5 _d 0 ) THEN
992ae64189 Mart*1089             EKnorm = EKnorm + 0.5 _d 0 * heff(i,j,bi,bj) *
                1090      &           ( ( (uIce(i,j,bi,bj)+uIce(i+1,j,bi,bj)) )**2
                1091      &           + ( (vIce(i,j,bi,bj)+vIce(i,j+1,bi,bj)) )**2 )
                1092             counter = counter + 1. _d 0
                1093            ENDIF
                1094           ENDDO
                1095          ENDDO
                1096         ENDDO
                1097        ENDDO
                1098        _GLOBAL_SUM_RL( resnorm, myThid )
                1099        _GLOBAL_SUM_RL( EKnorm, myThid )
                1100        _GLOBAL_SUM_RL( counter, myThid )
8377b8ee87 Mart*1101        IF ( counter .GT. 0. _d 0 ) EKnorm = EKnorm/counter
86fa82a64d Jean*1102        _BEGIN_MASTER( myThid )
                1103        WRITE(standardMessageUnit,'(A,I7,1X,2E22.14)')
                1104      &      'S/R SEAICE_LSR: IPSEUDO, RESNORM, EKNORM = ',
8377b8ee87 Mart*1105      &      ipass, SQRT(resnorm), EKnorm
86fa82a64d Jean*1106        _END_MASTER( myThid )
992ae64189 Mart*1107       ENDIF
                1108 #endif /* SEAICE_ALLOW_CHECK_LSR_CONVERGENCE */
                1109
143d9ce879 Dami*1110 C     second IF doNonLinLoop
992ae64189 Mart*1111       ENDIF
3e8bf7e743 Mart*1112 C     end outer ipass pseudo-time stepping loop
992ae64189 Mart*1113       ENDDO
                1114
143d9ce879 Dami*1115 #ifdef ALLOW_DIAGNOSTICS
                1116       IF ( DIAGNOSTICS_IS_ON('SIlsrRe ',myThid) ) THEN
                1117 C     S/R SEAICE_CALC_RESIDUAL does not compute residual velocities on
                1118 C     halos, so we need an exchange here to get the correct diagnostics
                1119 C     at cell centers
                1120        CALL EXCH_UV_XY_RL( uTmp, vTmp,.TRUE.,myThid)
                1121        DO bj = myByLo(myThid), myByHi(myThid)
                1122         DO bi = myBxLo(myThid), myBxHi(myThid)
                1123          DO j=1,sNy
                1124           DO i=1,sNx
                1125            RTmp(i,j) = SQRT( HALF * (
                1126      &            uTmp(i,  j,bi,bj) * uTmp(i,  j,bi,bj)
                1127      &          + uTmp(i+1,j,bi,bj) * uTmp(i+1,j,bi,bj)
                1128      &          + vTmp(i,j,  bi,bj) * vTmp(i,j,  bi,bj)
                1129      &          + vTmp(i,j+1,bi,bj) * vTmp(i,j+1,bi,bj)
                1130      &          ) )
                1131           ENDDO
                1132          ENDDO
                1133          CALL DIAGNOSTICS_FILL(RTmp,'SIlsrRe ',0,1,3,bi,bj,myThid)
                1134         ENDDO
                1135        ENDDO
                1136       ENDIF
                1137 #endif   /* ALLOW_DIAGNOSTICS */
                1138
992ae64189 Mart*1139       IF ( useHB87StressCoupling ) THEN
dfc97bd020 Mart*1140 # ifdef ALLOW_AUTODIFF_TAMC
                1141 C     These directives do not remove any recomputation warnings but avoid
                1142 C     recomputing the entire LSR loop before evaluating this code block
cbf501ab81 Jean*1143 CADJ STORE uice, vice, eta, etaz, zeta = comlev1_dynsol,
dfc97bd020 Mart*1144 CADJ &     kind=isbyte, key = ikey_dynamics
                1145 # endif /* ALLOW_AUTODIFF_TAMC */
8b339058e0 Mart*1146 C     compute the divergence of stress here to be used later
35fda33b05 Jean*1147 C
8b339058e0 Mart*1148 C     compute strain rate from latest velocities
                1149        CALL SEAICE_CALC_STRAINRATES(
                1150      I       uIce, vIce,
                1151      O       e11, e22, e12,
2e75568507 Mart*1152      I       3, myTime, myIter, myThid )
8b339058e0 Mart*1153 C     compute internal stresses with updated ice velocities
925df4f4b9 Mart*1154 C     and evaluate divergence of stress and apply to forcing
8b339058e0 Mart*1155        DO bj=myByLo(myThid),myByHi(myThid)
                1156         DO bi=myBxLo(myThid),myBxHi(myThid)
cbf501ab81 Jean*1157          CALL SEAICE_CALC_STRESSDIV(
925df4f4b9 Mart*1158      I        e11, e22, e12, press, zeta, eta, etaZ,
                1159      O        stressDivergenceX, stressDivergenceY,
                1160      I        bi, bj, myTime, myIter, myThid )
8b339058e0 Mart*1161         ENDDO
964a494b4b Mart*1162        ENDDO
35fda33b05 Jean*1163 C     endif  useHB87StressCoupling
8b339058e0 Mart*1164       ENDIF
                1165
53092bcb42 Mart*1166       RETURN
                1167       END
efcb5efb58 Jean*1168
                1169 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                1170 CBOP
                1171 C     !ROUTINE: SEAICE_RESIDUAL
                1172 C     !INTERFACE:
                1173       SUBROUTINE SEAICE_RESIDUAL(
                1174      I                  rhsU, rhsV, uRt1, uRt2, vRt1, vRt2,
                1175      I                  AU, BU, CU, AV, BV, CV, uFld, vFld,
                1176      O                  residU, residV, uRes, vRes,
                1177      I                  calcMeanResid, myIter, myThid )
                1178
                1179 C     !DESCRIPTION: \bv
                1180 C     *==========================================================*
                1181 C     | S/R SEAICE_RESIDUAL
                1182 C     | o Compute Linear solver residual
                1183 C     *==========================================================*
                1184 C     \ev
                1185
                1186 C     !USES:
                1187       IMPLICIT NONE
                1188
                1189 C     === Global variables ===
                1190 #include "SIZE.h"
                1191 #include "EEPARAMS.h"
                1192 #include "GRID.h"
03c669d1ab Jean*1193 #include "SEAICE_SIZE.h"
efcb5efb58 Jean*1194 #include "SEAICE.h"
                1195
                1196 C     !INPUT/OUTPUT PARAMETERS:
                1197 C     === Routine arguments ===
                1198 C     calcMeanResid :: also compute mean residual (=RMS)
                1199 C     myIter        :: Simulation timestep number
                1200 C     myThid        :: my Thread Id. number
                1201 C--   RHS
03c669d1ab Jean*1202       _RL rhsU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1203       _RL rhsV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
efcb5efb58 Jean*1204 C-    coefficients for lateral points, u(j+/-1)
03c669d1ab Jean*1205       _RL uRt1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1206       _RL uRt2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
efcb5efb58 Jean*1207 C-    coefficients for lateral points, v(i+/-1)
03c669d1ab Jean*1208       _RL vRt1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1209       _RL vRt2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
efcb5efb58 Jean*1210 C-    diagonals of coefficient matrices
03c669d1ab Jean*1211       _RL AU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1212       _RL BU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1213       _RL CU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1214       _RL AV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1215       _RL BV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1216       _RL CV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
efcb5efb58 Jean*1217 C--   seaice velocity
03c669d1ab Jean*1218       _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1219       _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
efcb5efb58 Jean*1220 C-    residual (output)
                1221       _RL residU, residV
03c669d1ab Jean*1222       _RL uRes (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1223       _RL vRes (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
efcb5efb58 Jean*1224       LOGICAL calcMeanResid
                1225       INTEGER myIter
                1226       INTEGER myThid
                1227 CEOP
                1228
                1229 C     !LOCAL VARIABLES:
                1230 C     === Local variables ===
                1231 C     i,j,bi,bj  :: Loop counters
                1232
                1233       INTEGER i, j, bi, bj
                1234
                1235 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                1236 C--   Calculate residual of the linearised system
                1237       residU = 0.
                1238       residV = 0.
                1239
                1240       DO bj=myByLo(myThid),myByHi(myThid)
                1241        DO bi=myBxLo(myThid),myBxHi(myThid)
                1242          DO j=1,sNy
                1243           DO i=1,sNx
                1244            uRes(i,j,bi,bj) = rhsU(i,j,bi,bj)
                1245      &            + uRt1(i,j,bi,bj)*uFld(i,j-1,bi,bj)
                1246      &            + uRt2(i,j,bi,bj)*uFld(i,j+1,bi,bj)
                1247      &            - ( AU(i,j,bi,bj)*uFld(i-1,j,bi,bj)
                1248      &              + BU(i,j,bi,bj)*uFld( i ,j,bi,bj)
                1249      &              + CU(i,j,bi,bj)*uFld(i+1,j,bi,bj)
                1250      &              )
                1251            vRes(i,j,bi,bj) = rhsV(i,j,bi,bj)
                1252      &            + vRt1(i,j,bi,bj)*vFld(i-1,j,bi,bj)
                1253      &            + vRt2(i,j,bi,bj)*vFld(i+1,j,bi,bj)
                1254      &            - ( AV(i,j,bi,bj)*vFld(i,j-1,bi,bj)
                1255      &              + BV(i,j,bi,bj)*vFld(i, j ,bi,bj)
                1256      &              + CV(i,j,bi,bj)*vFld(i,j+1,bi,bj)
                1257      &              )
                1258           ENDDO
                1259          ENDDO
56d0d0c8e2 Mart*1260          IF ( calcMeanResid ) THEN
                1261           DO j=1,sNy
                1262            DO i=1,sNx
efcb5efb58 Jean*1263             residU = residU + uRes(i,j,bi,bj)*uRes(i,j,bi,bj)
                1264      &                       *rAw(i,j,bi,bj)*maskInW(i,j,bi,bj)
5adf1d4e76 Mart*1265 #ifndef OBCS_UVICE_OLD
                1266      &                       *maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj)
                1267 #endif
efcb5efb58 Jean*1268             residV = residV + vRes(i,j,bi,bj)*vRes(i,j,bi,bj)
                1269      &                       *rAs(i,j,bi,bj)*maskInS(i,j,bi,bj)
5adf1d4e76 Mart*1270 #ifndef OBCS_UVICE_OLD
                1271      &                       *maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj)
                1272 #endif
56d0d0c8e2 Mart*1273            ENDDO
efcb5efb58 Jean*1274           ENDDO
56d0d0c8e2 Mart*1275          ENDIF
efcb5efb58 Jean*1276        ENDDO
                1277       ENDDO
                1278       IF ( calcMeanResid ) THEN
56d0d0c8e2 Mart*1279        _GLOBAL_SUM_RL( residU, myThid )
                1280        _GLOBAL_SUM_RL( residV, myThid )
                1281 C     scale residuals by globalArea so that they do not get ridiculously large
                1282        IF ( residU.GT.0. ) residU = SQRT(residU/globalArea)
                1283        IF ( residV.GT.0. ) residV = SQRT(residV/globalArea)
efcb5efb58 Jean*1284       ENDIF
                1285
                1286 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                1287
                1288       RETURN
                1289       END
b7c2bb63b7 Mart*1290
                1291 CBOP
5fb12930f9 Mart*1292 C     !ROUTINE: SEAICE_LSR_CALC_COEFFS
b7c2bb63b7 Mart*1293 C     !INTERFACE:
cbf501ab81 Jean*1294       SUBROUTINE SEAICE_LSR_CALC_COEFFS(
df1dac8b7b Mart*1295      I     etaPlusZeta, zetaMinusEta, etaZloc, zetaZloc, dragSym,
b7c2bb63b7 Mart*1296      O     AU, BU, CU, AV, BV, CV, uRt1, uRt2, vRt1, vRt2,
d585a5aee9 Mart*1297      I     iMin, iMax, jMin, jMax, myTime, myIter, myThid )
b7c2bb63b7 Mart*1298
                1299 C     !DESCRIPTION: \bv
                1300 C     *==========================================================*
5fb12930f9 Mart*1301 C     | S/R SEAICE_LSR_CALC_COEFFS
b7c2bb63b7 Mart*1302 C     | o Calculate coefficient matrix for LSR solver
                1303 C     *==========================================================*
                1304 C     \ev
                1305
                1306 C     !USES:
                1307       IMPLICIT NONE
                1308
                1309 C     === Global variables ===
                1310 #include "SIZE.h"
                1311 #include "EEPARAMS.h"
                1312 #include "PARAMS.h"
                1313 #include "GRID.h"
                1314 #include "SEAICE_SIZE.h"
                1315 #include "SEAICE_PARAMS.h"
3f0f10fc37 Mart*1316 #include "SEAICE_GRID.h"
b7c2bb63b7 Mart*1317 #include "SEAICE.h"
                1318
                1319 C     !INPUT/OUTPUT PARAMETERS:
                1320 C     === Routine arguments ===
                1321 C     myTime     :: Simulation time
                1322 C     myIter     :: Simulation timestep number
                1323 C     myThid     :: my Thread Id. number
                1324       _RL     myTime
                1325       INTEGER myIter
                1326       INTEGER myThid
d585a5aee9 Mart*1327       INTEGER iMin, iMax, jMin, jMax
b7c2bb63b7 Mart*1328 C     combinations of bulk and shear viscosity at C points
                1329 C     abbreviations
                1330       _RL etaPlusZeta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1331       _RL zetaMinusEta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0fd6b2de1a Mart*1332 C     shear and bulk viscosity at Z points
                1333       _RL  etaZloc    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1334       _RL zetaZloc    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
df1dac8b7b Mart*1335 C     symmetric ice-ocean stress coefficient
                1336       _RL dragSym     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
b7c2bb63b7 Mart*1337 C     coefficients of ice velocities in coefficient matrix
                1338 C     for both U and V-equation
                1339 C     diagonals of coefficient matrices
                1340       _RL AU          (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1341       _RL BU          (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1342       _RL CU          (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1343       _RL AV          (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1344       _RL BV          (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1345       _RL CV          (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1346 C     coefficients for lateral points, u(j+/-1)
                1347       _RL uRt1        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1348       _RL uRt2        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1349 C     coefficients for lateral points, v(i+/-1)
                1350       _RL vRt1        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1351       _RL vRt2        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1352 CEOP
                1353
                1354 C     !LOCAL VARIABLES:
                1355 C     === Local variables ===
                1356 C     i,j,bi,bj  :: Loop counters
                1357
                1358       INTEGER i, j, bi, bj
                1359       _RL hFacM, hFacP
e232688b80 Mart*1360 C     backward difference extrapolation factor (includes 1/deltaT)
                1361       _RL bdfAlphaOverDt
0fd6b2de1a Mart*1362 C     strongly implicit coupling factor
                1363       _RL strImpCplFac
b7c2bb63b7 Mart*1364
70e078b38a Mart*1365 C     fractional area at velocity points
                1366       _RL areaW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                1367       _RL areaS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
b7c2bb63b7 Mart*1368 C     coefficients of ice velocities in coefficient matrix
                1369 C     for both U and V-equation
                1370 C     XX: double derivative in X
                1371 C     YY: double derivative in Y
                1372 C     XM: metric term with derivative in X
                1373 C     YM: metric term with derivative in Y
bfe6f7447e Mart*1374       _RL UXX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                1375       _RL UYY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                1376       _RL UXM  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                1377       _RL UYM  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                1378       _RL VXX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                1379       _RL VYY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                1380       _RL VXM  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                1381       _RL VYM  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
b7c2bb63b7 Mart*1382 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                1383
e501eee760 Mart*1384 C     backward difference extrapolation factor
e232688b80 Mart*1385       bdfAlphaOverDt = 1. _d 0
e501eee760 Mart*1386       IF ( SEAICEuseBDF2 ) THEN
                1387        IF ( myIter.EQ.nIter0 .AND. SEAICEmomStartBDF.EQ.0 ) THEN
e232688b80 Mart*1388         bdfAlphaOverDt = 1. _d 0
6cbc659de0 Mart*1389        ELSE
e232688b80 Mart*1390         bdfAlphaOverDt = 1.5 _d 0
6cbc659de0 Mart*1391        ENDIF
                1392       ENDIF
e232688b80 Mart*1393 C
                1394       bdfAlphaOverDt = bdfAlphaOverDt / SEAICE_deltaTdyn
                1395 C
0fd6b2de1a Mart*1396       strImpCplFac = 0. _d 0
                1397       IF ( SEAICEuseStrImpCpl ) strImpCplFac = 1. _d 0
                1398
b7c2bb63b7 Mart*1399       DO bj=myByLo(myThid),myByHi(myThid)
                1400        DO bi=myBxLo(myThid),myBxHi(myThid)
70e078b38a Mart*1401         IF ( SEAICEscaleSurfStress ) THEN
4a08d54d3a Mart*1402          DO j=jMin,jMax
                1403           DO i=iMin,iMax
                1404            areaW(i,j) = 0.5 _d 0*(AREA(i,j,bi,bj)+AREA(i-1,j,bi,bj))
                1405            areaS(i,j) = 0.5 _d 0*(AREA(i,j,bi,bj)+AREA(i,j-1,bi,bj))
70e078b38a Mart*1406           ENDDO
                1407          ENDDO
                1408         ELSE
4a08d54d3a Mart*1409          DO j=jMin,jMax
                1410           DO i=iMin,iMax
                1411            areaW(i,j) = 1. _d 0
                1412            areaS(i,j) = 1. _d 0
70e078b38a Mart*1413           ENDDO
                1414          ENDDO
                1415         ENDIF
b7c2bb63b7 Mart*1416 C
4a08d54d3a Mart*1417 C     coefficients of uIce(i,j) and vIce(i,j) belonging to ...
b7c2bb63b7 Mart*1418 C
4a08d54d3a Mart*1419         DO j=jMin,jMax
                1420          DO i=iMin-1,iMax
b7c2bb63b7 Mart*1421 C     ... d/dx (eta+zeta) d/dx u
4a08d54d3a Mart*1422           UXX(i,j) = _dyF(i,j,bi,bj) * etaPlusZeta(i,j,bi,bj)
                1423      &         * _recip_dxF(i,j,bi,bj)
b7c2bb63b7 Mart*1424 C     ... d/dx (zeta-eta) k1 u
4a08d54d3a Mart*1425           UXM(i,j) = _dyF(i,j,bi,bj) * zetaMinusEta(i,j,bi,bj)
                1426      &         * k1AtC(i,j,bi,bj) * 0.5 _d 0
b7c2bb63b7 Mart*1427          ENDDO
                1428         ENDDO
4a08d54d3a Mart*1429         DO j=jMin,jMax+1
                1430          DO i=iMin,iMax
b7c2bb63b7 Mart*1431 C     ... d/dy eta d/dy u
4a08d54d3a Mart*1432           UYY(i,j) = _dxV(i,j,bi,bj) *
                1433      &         ( etaZloc(i,j,bi,bj) + strImpCplFac*zetaZloc(i,j,bi,bj) )
                1434      &         * _recip_dyU(i,j,bi,bj)
b7c2bb63b7 Mart*1435 C     ... d/dy eta k2 u
4a08d54d3a Mart*1436           UYM(i,j) = _dxV(i,j,bi,bj) * etaZloc(i,j,bi,bj)
                1437      &         * k2AtZ(i,j,bi,bj) * 0.5 _d 0
b7c2bb63b7 Mart*1438          ENDDO
                1439         ENDDO
4a08d54d3a Mart*1440         DO j=jMin,jMax
                1441          DO i=iMin,iMax+1
b7c2bb63b7 Mart*1442 C     ... d/dx eta dv/dx
4a08d54d3a Mart*1443           VXX(i,j) = _dyU(i,j,bi,bj) *
                1444      &         ( etaZloc(i,j,bi,bj) + strImpCplFac*zetaZloc(i,j,bi,bj) )
                1445      &         * _recip_dxV(i,j,bi,bj)
b7c2bb63b7 Mart*1446 C     ... d/dx eta k1 v
4a08d54d3a Mart*1447           VXM(i,j) = _dyU(i,j,bi,bj) * etaZloc(i,j,bi,bj)
                1448      &         * k1AtZ(i,j,bi,bj) * 0.5 _d 0
b7c2bb63b7 Mart*1449          ENDDO
                1450         ENDDO
4a08d54d3a Mart*1451         DO j=jMin-1,jMax
                1452          DO i=iMin,iMax
b7c2bb63b7 Mart*1453 C     ... d/dy eta+zeta dv/dy
4a08d54d3a Mart*1454           VYY(i,j) = _dxF(i,j,bi,bj) * etaPlusZeta(i,j,bi,bj)
                1455      &         * _recip_dyF(i,j,bi,bj)
b7c2bb63b7 Mart*1456 C     ... d/dy (zeta-eta) k2 v
4a08d54d3a Mart*1457           VYM(i,j) = _dxF(i,j,bi,bj) * zetaMinusEta(i,j,bi,bj)
                1458      &         * k2AtC(i,j,bi,bj) * 0.5 _d 0
b7c2bb63b7 Mart*1459          ENDDO
                1460         ENDDO
                1461
                1462 C--   Coefficients for solving uIce :
                1463
                1464 C     assemble coefficient matrix, beware of sign convention: because this
                1465 C     is the left hand side we calculate -grad(sigma), but the coefficients
4a08d54d3a Mart*1466 C     of U(i,j+/-1) are counted on the right hand side
                1467         DO j=jMin,jMax
                1468          DO i=iMin,iMax
                1469 C     coefficients for UICE(i-1,j)
                1470           AU(i,j,bi,bj)= ( - UXX(i-1,j) + UXM(i-1,j) )
                1471      &         * seaiceMaskU(i,j,bi,bj)
                1472 C     coefficients for UICE(i+1,j)
                1473           CU(i,j,bi,bj)= ( - UXX(i  ,j) - UXM(i  ,j) )
                1474      &         * seaiceMaskU(i,j,bi,bj)
                1475 C     coefficients for UICE(i,j)
                1476           BU(i,j,bi,bj)=(ONE - seaiceMaskU(i,j,bi,bj)) +
                1477      &         ( UXX(i-1,j) + UXX(i,j) + UYY(i,j+1) + UYY(i,j)
                1478      &         + UXM(i-1,j) - UXM(i,j) + UYM(i,j+1) - UYM(i,j)
                1479      &         ) * seaiceMaskU(i,j,bi,bj)
                1480 C     coefficients of uIce(i,j-1)
                1481           uRt1(i,j,bi,bj)= UYY(i,j  ) + UYM(i,j  )
                1482 C     coefficients of uIce(i,j+1)
                1483           uRt2(i,j,bi,bj)= UYY(i,j+1) - UYM(i,j+1)
b7c2bb63b7 Mart*1484          ENDDO
                1485         ENDDO
                1486
                1487 C     apply boundary conditions according to slip factor
                1488 C     for no slip, set u on boundary to zero: u(j+/-1)=-u(j)
                1489 C     for the free slip case sigma_12 = 0
4a08d54d3a Mart*1490         DO j=jMin,jMax
                1491          DO i=iMin,iMax
                1492           hFacM = seaiceMaskU(i,j-1,bi,bj)
                1493           hFacP = seaiceMaskU(i,j+1,bi,bj)
                1494 C     copy contributions to coefficient of U(i,j)
b7c2bb63b7 Mart*1495 C     beware of sign convection: uRt1/2 have the opposite sign convention
                1496 C     than BU, hence the minus sign
4a08d54d3a Mart*1497           BU(i,j,bi,bj)=BU(i,j,bi,bj) + seaiceMaskU(i,j,bi,bj) *
                1498      &         ( ( 1. _d 0 - hFacM ) * ( UYY(i  ,j  ) + UYM(i  ,j  ) )
                1499      &         + ( 1. _d 0 - hFacP ) * ( UYY(i  ,j+1) - UYM(i  ,j+1) ) )
                1500 C     reset coefficients of U(i,j-1) and U(i,j+1)
                1501           uRt1(i,j,bi,bj) = uRt1(i,j,bi,bj) * hFacM
                1502           uRt2(i,j,bi,bj) = uRt2(i,j,bi,bj) * hFacP
b7c2bb63b7 Mart*1503          ENDDO
                1504         ENDDO
                1505
                1506 C     now we need to normalize everything by the grid cell area
4a08d54d3a Mart*1507         DO j=jMin,jMax
                1508          DO i=iMin,iMax
                1509           AU(i,j,bi,bj)   = AU(i,j,bi,bj)   * recip_rAw(i,j,bi,bj)
                1510           CU(i,j,bi,bj)   = CU(i,j,bi,bj)   * recip_rAw(i,j,bi,bj)
cbf501ab81 Jean*1511 C     here we need to add the contribution from the time derivative (in
                1512 C     bdfAlphaOverDt) and the symmetric drag term; must be done after
e232688b80 Mart*1513 C     normalizing
4a08d54d3a Mart*1514           BU(i,j,bi,bj)    = BU(i,j,bi,bj)  * recip_rAw(i,j,bi,bj)
                1515      &         + seaiceMaskU(i,j,bi,bj) *
                1516      &         ( bdfAlphaOverDt*seaiceMassU(i,j,bi,bj)
8377b8ee87 Mart*1517      &         + 0.5 _d 0 * ( dragSym(i,  j,bi,bj)
4a08d54d3a Mart*1518      &                      + dragSym(i-1,j,bi,bj) )*areaW(i,j)
df1dac8b7b Mart*1519      &         )
4a08d54d3a Mart*1520           uRt1(i,j,bi,bj) = uRt1(i,j,bi,bj) * recip_rAw(i,j,bi,bj)
                1521           uRt2(i,j,bi,bj) = uRt2(i,j,bi,bj) * recip_rAw(i,j,bi,bj)
b7c2bb63b7 Mart*1522          ENDDO
                1523         ENDDO
                1524
                1525 C--   Coefficients for solving uIce :
                1526
                1527 C     assemble coefficient matrix, beware of sign convention: because this
                1528 C     is the left hand side we calculate -grad(sigma), but the coefficients
4a08d54d3a Mart*1529 C     of V(i+/-1,j) are counted on the right hand side
                1530         DO j=jMin,jMax
                1531          DO i=iMin,iMax
                1532 C     coefficients for VICE(i,j-1)
                1533           AV(i,j,bi,bj)=( - VYY(i,j-1) + VYM(i,j-1)
                1534      &         ) * seaiceMaskV(i,j,bi,bj)
                1535 C     coefficients for VICE(i,j+1)
                1536           CV(i,j,bi,bj)=( - VYY(i,j  ) - VYM(i,j  )
                1537      &         ) * seaiceMaskV(i,j,bi,bj)
                1538 C     coefficients for VICE(i,j)
                1539           BV(i,j,bi,bj)= (ONE - seaiceMaskV(i,j,bi,bj)) +
                1540      &         ( VXX(i,j) + VXX(i+1,j) + VYY(i,j) + VYY(i,j-1)
                1541      &         - VXM(i,j) + VXM(i+1,j) - VYM(i,j) + VYM(i,j-1)
                1542      &         ) * seaiceMaskV(i,j,bi,bj)
                1543 C     coefficients for V(i-1,j)
                1544           vRt1(i,j,bi,bj) = VXX(i  ,j) + VXM(i  ,j)
                1545 C     coefficients for V(i+1,j)
                1546           vRt2(i,j,bi,bj) = VXX(i+1,j) - VXM(i+1,j)
b7c2bb63b7 Mart*1547          ENDDO
                1548         ENDDO
                1549
                1550 C     apply boundary conditions according to slip factor
                1551 C     for no slip, set u on boundary to zero: v(i+/-1)=-v(i)
                1552 C     for the free slip case sigma_12 = 0
4a08d54d3a Mart*1553         DO j=jMin,jMax
                1554          DO i=iMin,iMax
b7c2bb63b7 Mart*1555           hFacM = seaiceMaskV(i-1,j,bi,bj)
                1556           hFacP = seaiceMaskV(i+1,j,bi,bj)
4a08d54d3a Mart*1557 C     copy contributions to coefficient of V(i,j)
b7c2bb63b7 Mart*1558 C     beware of sign convection: vRt1/2 have the opposite sign convention
                1559 C     than BV, hence the minus sign
4a08d54d3a Mart*1560           BV(i,j,bi,bj)=BV(i,j,bi,bj) + seaiceMaskV(i,j,bi,bj) *
                1561      &         ( ( 1. _d 0 - hFacM ) * ( VXX(i  ,j) + VXM(i  ,j) )
                1562      &         + ( 1. _d 0 - hFacP ) * ( VXX(i+1,j) - VXM(i+1,j) ) )
                1563 C     reset coefficients of V(i-1,j) and V(i+1,j)
                1564           vRt1(i,j,bi,bj) = vRt1(i,j,bi,bj) * hFacM
                1565           vRt2(i,j,bi,bj) = vRt2(i,j,bi,bj) * hFacP
b7c2bb63b7 Mart*1566          ENDDO
                1567         ENDDO
                1568
                1569 C     now we need to normalize everything by the grid cell area
4a08d54d3a Mart*1570         DO j=jMin,jMax
                1571          DO i=iMin,iMax
                1572           AV(i,j,bi,bj)   = AV(i,j,bi,bj)   * recip_rAs(i,j,bi,bj)
                1573           CV(i,j,bi,bj)   = CV(i,j,bi,bj)   * recip_rAs(i,j,bi,bj)
e232688b80 Mart*1574 C     here we need add the contribution from the time derivative (in
cbf501ab81 Jean*1575 C     bdfAlphaOverDt) and the symmetric drag term; must be done after
e232688b80 Mart*1576 C     normalizing
4a08d54d3a Mart*1577           BV(i,j,bi,bj)   = BV(i,j,bi,bj)   * recip_rAs(i,j,bi,bj)
                1578      &         + seaiceMaskV(i,j,bi,bj) *
                1579      &         ( bdfAlphaOverDt*seaiceMassV(i,j,bi,bj)
                1580      &         + 0.5 _d 0 * ( dragSym(i,j,  bi,bj)
                1581      &                      + dragSym(i,j-1,bi,bj) )*areaS(i,j)
df1dac8b7b Mart*1582      &         )
4a08d54d3a Mart*1583           vRt1(i,j,bi,bj) = vRt1(i,j,bi,bj) * recip_rAs(i,j,bi,bj)
                1584           vRt2(i,j,bi,bj) = vRt2(i,j,bi,bj) * recip_rAs(i,j,bi,bj)
b7c2bb63b7 Mart*1585          ENDDO
                1586         ENDDO
                1587
70e078b38a Mart*1588         IF ( ( useCubedSphereExchange .AND.
                1589      &       (SEAICE_OLx.GT.0 .OR. SEAICE_OLy.GT.0) ) .OR.
                1590      &       SEAICEscaleSurfStress ) THEN
05180ad1ef Mart*1591 C     this is clearly a hack: make sure that diagonals BU/BV are non-zero.
                1592 C     In my experience we only need to catch the case of SEAICE_OLx/y=OLx/y-2,
                1593 C     but for safety lets call this routine whenever SEAICE_OLx/y > 0
                1594 C     &       (SEAICE_OLx.EQ.OLx-2 .OR. SEAICE_OLy.EQ.OLy-2) ) THEN
70e078b38a Mart*1595 C     When scaling the surface ice-ocean stress by AREA, then there will
                1596 C     be many cases of zero diagonal entries.
4a08d54d3a Mart*1597          DO j=jMin,jMax
                1598           DO i=iMin,iMax
                1599            IF (BU(i,j,bi,bj).EQ.0 _d 0) BU(i,j,bi,bj) = 1. _d 0
                1600            IF (BV(i,j,bi,bj).EQ.0 _d 0) BV(i,j,bi,bj) = 1. _d 0
05180ad1ef Mart*1601           ENDDO
                1602          ENDDO
                1603         ENDIF
b7c2bb63b7 Mart*1604 C     bi/bj-loops
                1605        ENDDO
                1606       ENDDO
                1607
5fb12930f9 Mart*1608       RETURN
                1609       END
                1610
                1611 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
bfe6f7447e Mart*1612
                1613 CBOP
                1614 C     !ROUTINE: SEAICE_LSR_RHSU
                1615 C     !INTERFACE:
                1616       SUBROUTINE SEAICE_LSR_RHSU(
cbf501ab81 Jean*1617      I     zetaMinusEta, etaPlusZeta, etaZloc, zetaZloc, pressLoc,
                1618      I     uIceC, vIceC,
bfe6f7447e Mart*1619      U     rhsU,
                1620      I     iMin, iMax, jMin, jMax, bi, bj, myThid )
                1621
                1622 C     !DESCRIPTION: \bv
                1623 C     *==========================================================*
                1624 C     | S/R SEAICE_LSR_RHSU
                1625 C     | o Set up stress contribution to rhs of u-momentum eq.
                1626 C     *==========================================================*
                1627 C     \ev
                1628
                1629 C     !USES:
                1630       IMPLICIT NONE
                1631 #include "SIZE.h"
0fd6b2de1a Mart*1632 #include "EEPARAMS.h"
bfe6f7447e Mart*1633 #include "GRID.h"
                1634 #include "SEAICE_SIZE.h"
0fd6b2de1a Mart*1635 #include "SEAICE_PARAMS.h"
3f0f10fc37 Mart*1636 #include "SEAICE_GRID.h"
bfe6f7447e Mart*1637 #include "SEAICE.h"
                1638
                1639 C     !INPUT/OUTPUT PARAMETERS:
                1640 C     === Routine arguments ===
                1641 C     myThid     :: my Thread Id. number
                1642       INTEGER myThid
                1643       INTEGER iMin, iMax, jMin, jMax, bi, bj
                1644 C
                1645       _RL zetaMinusEta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1646       _RL etaPlusZeta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1647       _RL etaZloc     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0fd6b2de1a Mart*1648       _RL zetaZloc    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
bfe6f7447e Mart*1649       _RL pressloc    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0fd6b2de1a Mart*1650       _RL uIceC       (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
bfe6f7447e Mart*1651       _RL vIceC       (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1652       _RL rhsU        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1653
                1654 C     !LOCAL VARIABLES:
                1655 C     === Local variables ===
                1656 C     i,j  :: Loop counters
3f0f10fc37 Mart*1657 C     eplus, eminus :: convenient abbreviations for e11+e22, e11-e22
4a08d54d3a Mart*1658       INTEGER i,j
3f0f10fc37 Mart*1659       _RL eplus, eminus
bfe6f7447e Mart*1660       _RL hFacM
                1661       _RL sig11(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                1662       _RL sig12(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                1663 CEOP
                1664
4a08d54d3a Mart*1665       DO j=1-OLy,sNy+OLy
                1666        DO i=1-OLx,sNx+OLx
                1667         sig11(i,j) = 0. _d 0
                1668         sig12(i,j) = 0. _d 0
bfe6f7447e Mart*1669        ENDDO
                1670       ENDDO
                1671 C     contribution of sigma11 to rhs
4a08d54d3a Mart*1672       DO j=jMin,jMax
                1673        DO i=iMin-1,iMax
                1674         sig11(i,j) = zetaMinusEta(i,j,bi,bj)
                1675      &       * ( vIceC(i,j+1,bi,bj) - vIceC(i,j,bi,bj) )
                1676      &       * _recip_dyF(i,j,bi,bj)
                1677      &       + etaPlusZeta(i,j,bi,bj) * k2AtC(i,j,bi,bj)
                1678      &       * 0.5 _d 0 * ( vIceC(i,j+1,bi,bj) + vIceC(i,j,bi,bj) )
                1679      &       - 0.5 _d 0 * pressLoc(i,j,bi,bj)
bfe6f7447e Mart*1680        ENDDO
                1681       ENDDO
                1682 C     contribution of sigma12 to rhs of u-equation
4a08d54d3a Mart*1683       DO j=jMin,jMax+1
                1684        DO i=iMin,iMax
                1685         hFacM = seaiceMaskV(i,j,bi,bj) - seaiceMaskV(i-1,j,bi,bj)
                1686         sig12(i,j) = etaZloc(i,j,bi,bj) * (
                1687      &       ( vIceC(i,j,bi,bj) - vIceC(i-1,j,bi,bj) )
                1688      &       * _recip_dxV(i,j,bi,bj)
                1689      &       - k1AtZ(i,j,bi,bj)
                1690      &       * 0.5 _d 0 * ( vIceC(i,j,bi,bj) + vIceC(i-1,j,bi,bj) )
bfe6f7447e Mart*1691      &       )
                1692 C     free slip conditions (sig12=0) are taken care of by masking sig12
4a08d54d3a Mart*1693      &       *HEFFM(i  ,j  ,bi,bj)*HEFFM(i-1,j  ,bi,bj)
                1694      &       *HEFFM(i  ,j-1,bi,bj)*HEFFM(i-1,j-1,bi,bj)
bfe6f7447e Mart*1695 C     no slip boundary conditions (v(i-1)=-v(i))
                1696 C     v(i)+v(i-1) = 0 is also taken care of by masking sig12, so that we
                1697 C     only need to deal with v(i)-v(i-1)
4a08d54d3a Mart*1698      &       + etaZloc(i,j,bi,bj) * _recip_dxV(i,j,bi,bj)
                1699      &       * ( vIceC(i,j,bi,bj) + vIceC(i-1,j,bi,bj) )
bfe6f7447e Mart*1700      &       * hFacM * 2. _d 0
                1701        ENDDO
                1702       ENDDO
                1703
0fd6b2de1a Mart*1704       IF ( SEAICEuseStrImpCpl ) THEN
cbf501ab81 Jean*1705 C     strictly speaking, this is not a contribution to sig12, but for
0fd6b2de1a Mart*1706 C     convenience we add the explicit term -zetaZ*du/dy here;
                1707 C     we do not have to add any metric terms, because they cancel.
4a08d54d3a Mart*1708        DO j=jMin,jMax+1
                1709         DO i=iMin,iMax
                1710          hFacM = seaiceMaskV(i,j,bi,bj) - seaiceMaskV(i-1,j,bi,bj)
                1711          sig12(i,j) = sig12(i,j) - zetaZloc(i,j,bi,bj) * (
                1712      &        ( uIceC(i,j,bi,bj) - uIceC(i,j-1,bi,bj) )
                1713      &        * _recip_dyU(i,j,bi,bj)
0fd6b2de1a Mart*1714      &        )
                1715 C     free slip conditions (sig12=0) are taken care of by masking sig12
4a08d54d3a Mart*1716      &        *HEFFM(i  ,j  ,bi,bj)*HEFFM(i-1,j  ,bi,bj)
                1717      &        *HEFFM(i  ,j-1,bi,bj)*HEFFM(i-1,j-1,bi,bj)
0fd6b2de1a Mart*1718 C     no slip boundary conditions (u(j-1)=-u(j))
                1719 C     u(j)+u(j-1) = 0 is also taken care of by masking sig12, so that we
                1720 C     only need to deal with u(j)-u(j-1)
4a08d54d3a Mart*1721      &        - zetaZloc(i,j,bi,bj) * _recip_dyU(i,j,bi,bj)
                1722      &        * ( uIceC(i,j,bi,bj) + uIceC(i,j-1,bi,bj) )
0fd6b2de1a Mart*1723      &        * hFacM * 2. _d 0
                1724         ENDDO
                1725        ENDDO
                1726       ENDIF
cbf501ab81 Jean*1727
4a08d54d3a Mart*1728       DO j=jMin,jMax
                1729        DO i=iMin,iMax
bfe6f7447e Mart*1730 C     contribution to the rhs part of grad(sigma)_x
4a08d54d3a Mart*1731         rhsU(i,j,bi,bj) = rhsU(i,j,bi,bj)
                1732      &       + recip_rAw(i,j,bi,bj) * seaiceMaskU(i,j,bi,bj) *
                1733      &       ( _dyF(i  ,j  ,bi,bj)*sig11(i  ,j  )
                1734      &       - _dyF(i-1,j  ,bi,bj)*sig11(i-1,j  )
                1735      &       + _dxV(i  ,j+1,bi,bj)*sig12(i  ,j+1)
                1736      &       - _dxV(i  ,j  ,bi,bj)*sig12(i  ,j  ) )
bfe6f7447e Mart*1737        ENDDO
                1738       ENDDO
                1739
3f0f10fc37 Mart*1740 C     add extra metric terms
                1741       IF ( SEAICEselectMetricTerms .GE. 2 ) THEN
                1742 C     to save memory we re-use sig11 as sig22, hopefully without confusion
                1743        DO j=1-OLy,sNy+OLy
                1744         DO i=1-OLx,sNx+OLx
                1745          eplus = e11(i,j,bi,bj) + e22(i,j,bi,bj)
                1746          eminus= e11(i,j,bi,bj) - e22(i,j,bi,bj)
                1747          sig11(i,j) = zeta(i,j,bi,bj)*eplus - eta(i,j,bi,bj)*eminus
                1748      &        - 0.5 _d 0 * press(i,j,bi,bj)
                1749          sig12(i,j) = 2. _d 0 * e12(i,j,bi,bj) * etaZ(i,j,bi,bj)
                1750         ENDDO
                1751        ENDDO
                1752
                1753        DO j=jMin,jMax
                1754         DO i=iMin,iMax
                1755          rhsU(i,j,bi,bj) = rhsU(i,j,bi,bj)
                1756      &        + seaiceMaskU(i,j,bi,bj) * halfRL *
                1757      &        ( k2AtU(i,j,bi,bj) * ( sig12(i,j) + sig12(i,j+1) )
                1758      &        - k1AtU(i,j,bi,bj) * ( sig11(i,j) + sig11(i-1,j) ) )
                1759         ENDDO
                1760        ENDDO
                1761       ENDIF
                1762
bfe6f7447e Mart*1763       RETURN
                1764       END
                1765
                1766 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                1767
                1768 CBOP
                1769 C     !ROUTINE: SEAICE_LSR_RHSV
                1770 C     !INTERFACE:
                1771       SUBROUTINE SEAICE_LSR_RHSV(
0fd6b2de1a Mart*1772      I     zetaMinusEta, etaPlusZeta, etaZloc, zetaZloc, pressLoc,
                1773      I     uIceC, vIceC,
bfe6f7447e Mart*1774      U     rhsV,
                1775      I     iMin, iMax, jMin, jMax, bi, bj, myThid )
                1776
                1777 C     !DESCRIPTION: \bv
                1778 C     *==========================================================*
                1779 C     | S/R SEAICE_LSR_RHSV
                1780 C     | o Set up stress contribution to rhs of v-momentum eq.
                1781 C     *==========================================================*
                1782 C     \ev
                1783
                1784 C     !USES:
                1785       IMPLICIT NONE
                1786 #include "SIZE.h"
0fd6b2de1a Mart*1787 #include "EEPARAMS.h"
bfe6f7447e Mart*1788 #include "GRID.h"
                1789 #include "SEAICE_SIZE.h"
0fd6b2de1a Mart*1790 #include "SEAICE_PARAMS.h"
3f0f10fc37 Mart*1791 #include "SEAICE_GRID.h"
bfe6f7447e Mart*1792 #include "SEAICE.h"
                1793
                1794 C     !INPUT/OUTPUT PARAMETERS:
                1795 C     === Routine arguments ===
                1796 C     myThid     :: my Thread Id. number
                1797       INTEGER myThid
                1798       INTEGER iMin, iMax, jMin, jMax, bi, bj
                1799 C
                1800       _RL zetaMinusEta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1801       _RL etaPlusZeta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1802       _RL etaZloc     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0fd6b2de1a Mart*1803       _RL zetaZloc    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
bfe6f7447e Mart*1804       _RL pressloc    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1805       _RL uIceC       (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0fd6b2de1a Mart*1806       _RL vIceC       (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
bfe6f7447e Mart*1807       _RL rhsV        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1808
                1809 C     !LOCAL VARIABLES:
                1810 C     === Local variables ===
                1811 C     i,j  :: Loop counters
3f0f10fc37 Mart*1812 C     eplus, eminus :: convenient abbreviations for e11+e22, e11-e22
4a08d54d3a Mart*1813       INTEGER i,j
3f0f10fc37 Mart*1814       _RL eplus, eminus
bfe6f7447e Mart*1815       _RL hFacM
                1816       _RL sig22(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                1817       _RL sig12(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                1818 CEOP
                1819
4a08d54d3a Mart*1820       DO j=1-OLy,sNy+OLy
                1821        DO i=1-OLx,sNx+OLx
                1822         sig22(i,j) = 0. _d 0
                1823         sig12(i,j) = 0. _d 0
bfe6f7447e Mart*1824        ENDDO
                1825       ENDDO
                1826
                1827 C     contribution of sigma22 to rhs
4a08d54d3a Mart*1828       DO j=jMin-1,jMax
                1829        DO i=iMin,iMax
                1830         sig22(i,j) = zetaMinusEta(i,j,bi,bj)
                1831      &       * ( uIceC(i+1,j,bi,bj) - uIceC(i,j,bi,bj) )
                1832      &       * _recip_dxF(i,j,bi,bj)
                1833      &       + etaPlusZeta(i,j,bi,bj) * k1AtC(i,j,bi,bj)
                1834      &       * 0.5 _d 0 * ( uIceC(i+1,j,bi,bj) + uIceC(i,j,bi,bj) )
                1835      &       - 0.5 _d 0 * pressLoc(i,j,bi,bj)
bfe6f7447e Mart*1836        ENDDO
                1837       ENDDO
                1838 C     contribution of sigma12 to rhs of v-equation
4a08d54d3a Mart*1839       DO j=jMin,jMax
                1840        DO i=iMin,iMax+1
bfe6f7447e Mart*1841         hFacM = seaiceMaskU(i,j,bi,bj) - seaiceMaskU(i,j-1,bi,bj)
4a08d54d3a Mart*1842         sig12(i,j) = etaZloc(i,j,bi,bj) * (
                1843      &       ( uIceC(i,j,bi,bj) - uIceC(i,j-1,bi,bj) )
                1844      &       * _recip_dyU(i,j,bi,bj)
                1845      &       - k2AtZ(i,j,bi,bj)
                1846      &       * 0.5 _d 0 * ( uIceC(i,j,bi,bj) + uIceC(i,j-1,bi,bj) )
bfe6f7447e Mart*1847      &       )
                1848 C     free slip conditions (sig12=0) are taken care of by masking sig12,
4a08d54d3a Mart*1849      &       *HEFFM(i  ,j  ,bi,bj)*HEFFM(i-1,j  ,bi,bj)
                1850      &       *HEFFM(i  ,j-1,bi,bj)*HEFFM(i-1,j-1,bi,bj)
bfe6f7447e Mart*1851 C     no slip boundary conditions (u(j-1)=-u(j))
                1852 C     u(j)+u(j-1) = 0 is also taken care of by masking sig12, so that we
                1853 C     only need to deal with u(j)-u(j-1)
4a08d54d3a Mart*1854      &       + etaZloc(i,j,bi,bj) * _recip_dyU(i,j,bi,bj)
                1855      &       * ( uIceC(i,j,bi,bj) + uIceC(i,j-1,bi,bj) )
bfe6f7447e Mart*1856      &         * hFacM * 2. _d 0
                1857        ENDDO
                1858       ENDDO
                1859
0fd6b2de1a Mart*1860       IF ( SEAICEuseStrImpCpl ) THEN
cbf501ab81 Jean*1861 C     strictly speaking, this is not a contribution to sig12, but for
0fd6b2de1a Mart*1862 C     convenience we add the explicit term -zetaZ*dv/dx here;
                1863 C     we do not have to add any metric terms, because they cancel.
4a08d54d3a Mart*1864        DO j=jMin,jMax
                1865         DO i=iMin,iMax+1
0fd6b2de1a Mart*1866          hFacM = seaiceMaskU(i,j,bi,bj) - seaiceMaskU(i,j-1,bi,bj)
4a08d54d3a Mart*1867          sig12(i,j) = sig12(i,j) - zetaZloc(i,j,bi,bj) * (
                1868      &        ( vIceC(i,j,bi,bj) - vIceC(i-1,j,bi,bj) )
                1869      &        * _recip_dxV(i,j,bi,bj)
0fd6b2de1a Mart*1870      &        )
                1871 C     free slip conditions (sig12=0) are taken care of by masking sig12
4a08d54d3a Mart*1872      &        *HEFFM(i  ,j  ,bi,bj)*HEFFM(i-1,j  ,bi,bj)
                1873      &        *HEFFM(i  ,j-1,bi,bj)*HEFFM(i-1,j-1,bi,bj)
0fd6b2de1a Mart*1874 C     no slip boundary conditions (v(i-1)=-v(i))
                1875 C     v(i)+v(i-1) = 0 is also taken care of by masking sig12, so that we
                1876 C     only need to deal with v(i)-v(i-1)
4a08d54d3a Mart*1877      &        - zetaZloc(i,j,bi,bj) * _recip_dxV(i,j,bi,bj)
                1878      &        * ( vIceC(i,j,bi,bj) + vIceC(i-1,j,bi,bj) )
0fd6b2de1a Mart*1879      &        * hFacM * 2. _d 0
                1880         ENDDO
                1881        ENDDO
                1882       ENDIF
                1883
4a08d54d3a Mart*1884       DO j=jMin,jMax
                1885        DO i=iMin,iMax
bfe6f7447e Mart*1886 C     contribution to the rhs part of grad(sigma)_y
4a08d54d3a Mart*1887         rhsV(i,j,bi,bj) = rhsV(i,j,bi,bj)
                1888      &       + recip_rAs(i,j,bi,bj) * seaiceMaskV(i,j,bi,bj) *
                1889      &       ( _dyU(i+1,j  ,bi,bj) * sig12(i+1,j  )
                1890      &       - _dyU(i  ,j  ,bi,bj) * sig12(i  ,j  )
                1891      &       + _dxF(i  ,j  ,bi,bj) * sig22(i  ,j  )
                1892      &       - _dxF(i  ,j-1,bi,bj) * sig22(i  ,j-1) )
bfe6f7447e Mart*1893        ENDDO
                1894       ENDDO
                1895
3f0f10fc37 Mart*1896       IF ( SEAICEselectMetricTerms .GE. 2 ) THEN
                1897 C     to save memory we re-use sig22 as sig11, hopefully without confusion
                1898        DO j=1-OLy,sNy+OLy
                1899         DO i=1-OLx,sNx+OLx
                1900          eplus = e11(i,j,bi,bj) + e22(i,j,bi,bj)
                1901          eminus= e11(i,j,bi,bj) - e22(i,j,bi,bj)
                1902          sig22(i,j) = zeta(i,j,bi,bj)*eplus + eta(i,j,bi,bj)*eminus
                1903      &        - 0.5 _d 0 * press(i,j,bi,bj)
                1904          sig12(i,j) = 2. _d 0 * e12(i,j,bi,bj) * etaZ(i,j,bi,bj)
                1905         ENDDO
                1906        ENDDO
                1907
                1908        DO j=jMin,jMax
                1909         DO i=iMin,iMax
                1910          rhsV(i,j,bi,bj) = rhsV(i,j,bi,bj)
                1911      &        + seaiceMaskU(i,j,bi,bj) * halfRL *
                1912      &        ( k1AtV(i,j,bi,bj) * ( sig12(i,j) + sig12(i+1,j) )
                1913      &        - k2AtV(i,j,bi,bj) * ( sig22(i,j) + sig22(i,j-1) ) )
                1914         ENDDO
                1915        ENDDO
                1916       ENDIF
                1917
bfe6f7447e Mart*1918       RETURN
                1919       END
                1920
                1921 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
5fb12930f9 Mart*1922
                1923 CBOP
                1924 C     !ROUTINE: SEAICE_LSR_TRIDIAGU
                1925 C     !INTERFACE:
                1926       SUBROUTINE SEAICE_LSR_TRIDIAGU(
8df0ccd026 Jean*1927      I     AU, BU, CU, uRt1, uRt2, rhsU, uTmp, seaiceMaskU, WFAU,
5fb12930f9 Mart*1928      U     uIce,
8d9be9cde5 Mart*1929      I     iMin, iMax, jMin, jMax, bi, bj, myTime, myIter, myThid )
5fb12930f9 Mart*1930
                1931 C     !DESCRIPTION: \bv
                1932 C     *==========================================================*
8d9be9cde5 Mart*1933 C     | S/R SEAICE_LSR_TRIDIAGU
5fb12930f9 Mart*1934 C     | o Solve tridiagonal problem in uIce for LSR solver
                1935 C     *==========================================================*
                1936 C     \ev
                1937
                1938 C     !USES:
                1939       IMPLICIT NONE
                1940 #include "SIZE.h"
                1941
                1942 C     !INPUT/OUTPUT PARAMETERS:
                1943 C     === Routine arguments ===
                1944 C     myTime     :: Simulation time
                1945 C     myIter     :: Simulation timestep number
                1946 C     myThid     :: my Thread Id. number
                1947       _RL     myTime
                1948       INTEGER myIter
                1949       INTEGER myThid
                1950       INTEGER iMin, iMax, jMin, jMax, bi, bj
                1951 C     diagonals of coefficient matrices
                1952       _RL AU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1953       _RL BU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1954       _RL CU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1955 C     coefficients for lateral points, v(j+/-1)
                1956       _RL uRt1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1957       _RL uRt2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1958 C     right hand side
                1959       _RL rhsU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
8df0ccd026 Jean*1960 C     velocity of previous iteration step
                1961       _RL uTmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
5fb12930f9 Mart*1962 C     on entry: first guess and on exit: solution
                1963       _RL uIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1964 C     mask
                1965       _RL seaiceMaskU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                1966 C     relaxation factor
                1967       _RL WFAU
                1968 CEOP
                1969
                1970 C     !LOCAL VARIABLES:
                1971 C     === Local variables ===
                1972 C     i,j  :: Loop counters
4a08d54d3a Mart*1973       INTEGER i,j,iM
032b010f0c Mart*1974       INTEGER jMinLoc, jStep
12dace3fd6 Mart*1975 #ifdef SEAICE_LSR_ZEBRA
4a08d54d3a Mart*1976       INTEGER k
032b010f0c Mart*1977       PARAMETER ( jStep = 2 )
                1978 #else
                1979       PARAMETER ( jStep = 1 )
12dace3fd6 Mart*1980 #endif /* SEAICE_LSR_ZEBRA */
7c0da9bbfa Mart*1981       _RL AA3, bet
5fb12930f9 Mart*1982       _RL URT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                1983       _RL CUU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
cbf501ab81 Jean*1984
dadeed8422 Mart*1985 C     Need to initialize local arrays
4a08d54d3a Mart*1986       DO j=1-OLy,sNy+OLy
                1987        DO i=1-OLx,sNx+OLx
                1988         URT(i,j) = 0. _d 0
                1989         CUU(i,j) = 0. _d 0
dadeed8422 Mart*1990        ENDDO
                1991       ENDDO
12dace3fd6 Mart*1992 #ifdef SEAICE_LSR_ZEBRA
4a08d54d3a Mart*1993       DO k=0,1
                1994       jMinLoc = jMin+k
8d9be9cde5 Mart*1995 #else
032b010f0c Mart*1996       jMinLoc = jMin
12dace3fd6 Mart*1997 #endif /* SEAICE_LSR_ZEBRA */
4a08d54d3a Mart*1998       DO j=jMinLoc,jMax,jStep
                1999        DO i=iMin,iMax
5fb12930f9 Mart*2000         AA3 = 0. _d 0
4a08d54d3a Mart*2001         IF (i.EQ.iMin) AA3 = AA3 - AU(i,j,bi,bj)*uIce(i-1,j,bi,bj)
                2002         IF (i.EQ.iMax) AA3 = AA3 - CU(i,j,bi,bj)*uIce(i+1,j,bi,bj)
cbf501ab81 Jean*2003
4a08d54d3a Mart*2004         URT(i,j)=rhsU(i,j,bi,bj)
5fb12930f9 Mart*2005      &       + AA3
4a08d54d3a Mart*2006      &       + uRt1(i,j,bi,bj)*uIce(i,j-1,bi,bj)
                2007      &       + uRt2(i,j,bi,bj)*uIce(i,j+1,bi,bj)
                2008         URT(i,j)=URT(i,j)* seaiceMaskU(i,j,bi,bj)
5fb12930f9 Mart*2009        ENDDO
032b010f0c Mart*2010 C     beginning of tridiagnoal solver
4a08d54d3a Mart*2011        DO i=iMin,iMax
                2012         CUU(i,j)=CU(i,j,bi,bj)
5fb12930f9 Mart*2013        ENDDO
4a08d54d3a Mart*2014        CUU(iMin,j)=CUU(iMin,j)/BU(iMin,j,bi,bj)
                2015        URT(iMin,j)=URT(iMin,j)/BU(iMin,j,bi,bj)
5fb12930f9 Mart*2016 #ifdef SEAICE_VECTORIZE_LSR
                2017       ENDDO
                2018 C     start a new loop with reversed order to support automatic vectorization
032b010f0c Mart*2019 CADJ loop = sequential
4a08d54d3a Mart*2020       DO i=iMin+1,iMax
                2021        iM=i-1
                2022        DO j=jMinLoc,jMax,jStep
5fb12930f9 Mart*2023 #else /* do not SEAICE_VECTORIZE_LSR */
032b010f0c Mart*2024 CADJ loop = sequential
4a08d54d3a Mart*2025        DO i=iMin+1,iMax
                2026         iM=i-1
5fb12930f9 Mart*2027 #endif /* SEAICE_VECTORIZE_LSR */
4a08d54d3a Mart*2028         bet      = BU(i,j,bi,bj)-AU(i,j,bi,bj)*CUU(iM,j)
                2029         CUU(i,j) = CUU(i,j)/bet
                2030         URT(i,j) = (URT(i,j)-AU(i,j,bi,bj)*URT(iM,j))/bet
5fb12930f9 Mart*2031        ENDDO
                2032 #ifdef SEAICE_VECTORIZE_LSR
                2033       ENDDO
032b010f0c Mart*2034 CADJ loop = sequential
4a08d54d3a Mart*2035       DO i=iMin,iMax-1
8377b8ee87 Mart*2036        iM=sNx-i
4a08d54d3a Mart*2037        DO j=jMinLoc,jMax,jStep
8d9be9cde5 Mart*2038 #else
032b010f0c Mart*2039 CADJ loop = sequential
4a08d54d3a Mart*2040        DO i=iMin,iMax-1
8377b8ee87 Mart*2041         iM=sNx-i
032b010f0c Mart*2042 #endif /* SEAICE_VECTORIZE_LSR */
4a08d54d3a Mart*2043         URT(iM,j)=URT(iM,j)-CUU(iM,j)*URT(iM+1,j)
5fb12930f9 Mart*2044        ENDDO
032b010f0c Mart*2045 #ifdef SEAICE_VECTORIZE_LSR
                2046       ENDDO
                2047 #endif /* SEAICE_VECTORIZE_LSR */
                2048 C     end of tridiagnoal solver, solution is URT
                2049 C     relaxation with WFAU
                2050 #ifdef SEAICE_VECTORIZE_LSR
                2051 C     go back to original order
4a08d54d3a Mart*2052       DO j=jMinLoc,jMax,jStep
032b010f0c Mart*2053 #endif /* SEAICE_VECTORIZE_LSR */
4a08d54d3a Mart*2054        DO i=iMin,iMax
                2055         uIce(i,j,bi,bj)=uTmp(i,j,bi,bj)
                2056      &       +WFAU*(URT(i,j)-uTmp(i,j,bi,bj))
5fb12930f9 Mart*2057        ENDDO
                2058       ENDDO
12dace3fd6 Mart*2059 #ifdef SEAICE_LSR_ZEBRA
8d9be9cde5 Mart*2060       ENDDO
12dace3fd6 Mart*2061 #endif /* SEAICE_LSR_ZEBRA */
5fb12930f9 Mart*2062
                2063       RETURN
                2064       END
                2065
                2066 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                2067
                2068 CBOP
                2069 C     !ROUTINE: SEAICE_LSR_TRIDIAGV
                2070 C     !INTERFACE:
                2071       SUBROUTINE SEAICE_LSR_TRIDIAGV(
cbf501ab81 Jean*2072      I     AV, BV, CV, vRt1, vRt2, rhsV, vTmp, seaiceMaskV, WFAV,
5fb12930f9 Mart*2073      U     vIce,
8d9be9cde5 Mart*2074      I     iMin, iMax, jMin, jMax, bi, bj, myTime, myIter, myThid )
5fb12930f9 Mart*2075
                2076 C     !DESCRIPTION: \bv
                2077 C     *==========================================================*
                2078 C     | S/R SEAICE_LSR_TRIDIAGV
                2079 C     | o Solve tridiagonal problem in vIce for LSR solver
                2080 C     *==========================================================*
                2081 C     \ev
                2082
                2083 C     !USES:
                2084       IMPLICIT NONE
                2085 #include "SIZE.h"
                2086
                2087 C     !INPUT/OUTPUT PARAMETERS:
                2088 C     === Routine arguments ===
                2089 C     myTime     :: Simulation time
                2090 C     myIter     :: Simulation timestep number
                2091 C     myThid     :: my Thread Id. number
                2092       _RL     myTime
                2093       INTEGER myIter
                2094       INTEGER myThid
                2095       INTEGER iMin, iMax, jMin, jMax, bi, bj
                2096 C     diagonals of coefficient matrices
                2097       _RL AV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                2098       _RL BV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                2099       _RL CV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                2100 C     coefficients for lateral points, v(j+/-1)
                2101       _RL vRt1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                2102       _RL vRt2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                2103 C     right hand side
                2104       _RL rhsV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
8df0ccd026 Jean*2105 C     velocity of previous iteration step
                2106       _RL vTmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
5fb12930f9 Mart*2107 C     on entry: first guess and on exit: solution
                2108       _RL vIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                2109 C     mask
                2110       _RL seaiceMaskV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                2111 C     relaxation factor
                2112       _RL WFAV
                2113 CEOP
                2114
                2115 C     !LOCAL VARIABLES:
                2116 C     === Local variables ===
                2117 C     i,j  :: Loop counters
4a08d54d3a Mart*2118       INTEGER i,j,jM
032b010f0c Mart*2119       INTEGER iMinLoc, iStep
12dace3fd6 Mart*2120 #ifdef SEAICE_LSR_ZEBRA
4a08d54d3a Mart*2121       INTEGER k
032b010f0c Mart*2122       PARAMETER ( iStep = 2 )
                2123 #else
                2124       PARAMETER ( iStep = 1 )
12dace3fd6 Mart*2125 #endif /* SEAICE_LSR_ZEBRA */
7c0da9bbfa Mart*2126       _RL AA3, bet
5fb12930f9 Mart*2127       _RL VRT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                2128       _RL CVV(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                2129
dadeed8422 Mart*2130 C     Need to initialize local arrays
4a08d54d3a Mart*2131       DO j=1-OLy,sNy+OLy
                2132        DO i=1-OLx,sNx+OLx
                2133         VRT(i,j) = 0. _d 0
                2134         CVV(i,j) = 0. _d 0
dadeed8422 Mart*2135        ENDDO
                2136       ENDDO
12dace3fd6 Mart*2137 #ifdef SEAICE_LSR_ZEBRA
4a08d54d3a Mart*2138       DO k=0,1
                2139       iMinLoc = iMin+k
032b010f0c Mart*2140 #else
                2141       iMinLoc = iMin
12dace3fd6 Mart*2142 #endif /*  SEAICE_LSR_ZEBRA */
8d9be9cde5 Mart*2143 #ifdef SEAICE_VECTORIZE_LSR
4a08d54d3a Mart*2144       DO j=jMin,jMax
                2145        DO i=iMinLoc,iMax,iStep
8d9be9cde5 Mart*2146 #else
4a08d54d3a Mart*2147       DO i=iMinLoc,iMax,iStep
                2148        DO j=jMin,jMax
8d9be9cde5 Mart*2149 #endif /* SEAICE_VECTORIZE_LSR */
5fb12930f9 Mart*2150         AA3 = 0. _d 0
4a08d54d3a Mart*2151         IF (j.EQ.jMin) AA3 = AA3 - AV(i,j,bi,bj)*vIce(i,j-1,bi,bj)
                2152         IF (j.EQ.jMax) AA3 = AA3 - CV(i,j,bi,bj)*vIce(i,j+1,bi,bj)
5fb12930f9 Mart*2153
4a08d54d3a Mart*2154         VRT(i,j)=rhsV(i,j,bi,bj)
5fb12930f9 Mart*2155      &       + AA3
4a08d54d3a Mart*2156      &       + vRt1(i,j,bi,bj)*vIce(i-1,j,bi,bj)
                2157      &       + vRt2(i,j,bi,bj)*vIce(i+1,j,bi,bj)
                2158         VRT(i,j)=VRT(i,j)* seaiceMaskV(i,j,bi,bj)
5fb12930f9 Mart*2159
032b010f0c Mart*2160 C     beginning of tridiagnoal solver
4a08d54d3a Mart*2161         CVV(i,j)=CV(i,j,bi,bj)
5fb12930f9 Mart*2162        ENDDO
8d9be9cde5 Mart*2163 #ifdef SEAICE_VECTORIZE_LSR
                2164       ENDDO
4a08d54d3a Mart*2165       DO i=iMinLoc,iMax,iStep
8d9be9cde5 Mart*2166 #endif /* SEAICE_VECTORIZE_LSR */
4a08d54d3a Mart*2167        CVV(i,jMin)=CVV(i,jMin)/BV(i,jMin,bi,bj)
                2168        VRT(i,jMin)=VRT(i,jMin)/BV(i,jMin,bi,bj)
8d9be9cde5 Mart*2169 #ifdef SEAICE_VECTORIZE_LSR
                2170       ENDDO
                2171 C     start a new loop with reversed order to support automatic vectorization
032b010f0c Mart*2172 CADJ loop = sequential
4a08d54d3a Mart*2173       DO j=jMin+1,jMax
                2174        jM=j-1
                2175        DO i=iMinLoc,iMax,iStep
8d9be9cde5 Mart*2176 #else /* do not SEAICE_VECTORIZE_LSR */
032b010f0c Mart*2177 CADJ loop = sequential
4a08d54d3a Mart*2178        DO j=jMin+1,jMax
                2179         jM=j-1
8d9be9cde5 Mart*2180 #endif /* SEAICE_VECTORIZE_LSR */
4a08d54d3a Mart*2181         bet      = BV(i,j,bi,bj)-AV(i,j,bi,bj)*CVV(i,jM)
                2182         CVV(i,j) = CVV(i,j)/bet
                2183         VRT(i,j) = (VRT(i,j)-AV(i,j,bi,bj)*VRT(i,jM))/bet
5fb12930f9 Mart*2184        ENDDO
8d9be9cde5 Mart*2185 #ifdef SEAICE_VECTORIZE_LSR
                2186       ENDDO
                2187 C     go back to original order
032b010f0c Mart*2188 CADJ loop = sequential
4a08d54d3a Mart*2189       DO j=jMin,jMax-1
                2190        jM=sNy-j
                2191        DO i=iMinLoc,iMax,iStep
032b010f0c Mart*2192 #else
                2193 CADJ loop = sequential
4a08d54d3a Mart*2194        DO j=jMin,jMax-1
                2195         jM=sNy-j
032b010f0c Mart*2196 #endif /* SEAICE_VECTORIZE_LSR */
4a08d54d3a Mart*2197         VRT(i,jM)=VRT(i,jM)-CVV(i,jM)*VRT(i,jM+1)
5fb12930f9 Mart*2198        ENDDO
032b010f0c Mart*2199 #ifdef SEAICE_VECTORIZE_LSR
                2200       ENDDO
                2201 C     end of tridiagnoal solver, solution is VRT
                2202 C     relaxation with WFAV
4a08d54d3a Mart*2203       DO j=jMin,jMax
                2204        DO i=iMinLoc,iMax,iStep
032b010f0c Mart*2205 #else
4a08d54d3a Mart*2206        DO j=jMin,jMax
032b010f0c Mart*2207 #endif /* SEAICE_VECTORIZE_LSR */
4a08d54d3a Mart*2208         vIce(i,j,bi,bj)=vTmp(i,j,bi,bj)
                2209      &       +WFAV*(VRT(i,j)-vTmp(i,j,bi,bj))
5fb12930f9 Mart*2210        ENDDO
                2211       ENDDO
12dace3fd6 Mart*2212 #ifdef SEAICE_LSR_ZEBRA
8d9be9cde5 Mart*2213       ENDDO
12dace3fd6 Mart*2214 #endif /* SEAICE_LSR_ZEBRA */
5fb12930f9 Mart*2215
b7c2bb63b7 Mart*2216 #endif /* SEAICE_CGRID */
                2217
                2218       RETURN
                2219       END