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