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