pkg/seaice/seaice_preconditioner.F

255a701086 Mart*0001 #include "SEAICE_OPTIONS.h"
                0002 #ifdef ALLOW_OBCS
                0003 # include "OBCS_OPTIONS.h"
                0004 #endif
                0005
5acccad966 Jean*0006 C--   File seaice_preconditioner.F:
                0007 C--   Contents
                0008 C--   o SEAICE_PRECONDITIONER
e780b26e64 Mart*0009 C--   o SEAICE_PRECOND_RHSU
                0010 C--   o SEAICE_PRECOND_RHSV
5acccad966 Jean*0011
                0012 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0013
255a701086 Mart*0014 CBOP
                0015 C     !ROUTINE: SEAICE_PRECONDITIONER
                0016 C     !INTERFACE:
5acccad966 Jean*0017       SUBROUTINE SEAICE_PRECONDITIONER(
                0018      U     duIce, dvIce,
48db1b87b2 Mart*0019      I     zetaPre, etaPre, etaZpre, zetaZpre, dwatPre,
255a701086 Mart*0020      I     newtonIter, krylovIter, myTime, myIter, myThid )
                0021
                0022 C     !DESCRIPTION: \bv
                0023 C     *==========================================================*
                0024 C     | SUBROUTINE SEAICE_PRECONDITIONER
                0025 C     | o Preconditioner for Jacobian-free Newton-Krylov solver,
                0026 C     |   compute improved first guess solution du/vIce, with
                0027 C     |   suboptimal solver, here LSOR
                0028 C     *==========================================================*
                0029 C     | written by Martin Losch, Oct 2012
                0030 C     *==========================================================*
                0031 C     \ev
                0032
                0033 C     !USES:
                0034       IMPLICIT NONE
                0035
                0036 C     === Global variables ===
                0037 #include "SIZE.h"
                0038 #include "EEPARAMS.h"
                0039 #include "PARAMS.h"
                0040 #include "DYNVARS.h"
                0041 #include "GRID.h"
                0042 #include "SEAICE_SIZE.h"
                0043 #include "SEAICE_PARAMS.h"
3f0f10fc37 Mart*0044 #include "SEAICE_GRID.h"
255a701086 Mart*0045 #include "SEAICE.h"
                0046
5acccad966 Jean*0047 C     !INPUT PARAMETERS:
255a701086 Mart*0048 C     === Routine arguments ===
                0049 C     myTime :: Simulation time
                0050 C     myIter :: Simulation timestep number
                0051 C     myThid :: my Thread Id. number
                0052 C     newtonIter :: current iterate of Newton iteration
                0053 C     krylovIter :: current iterate of Krylov iteration
5acccad966 Jean*0054 C     *Pre are precomputed and held fixed during the Krylov iteration
48db1b87b2 Mart*0055       _RL   zetaPre(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0056       _RL  zetaZPre(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0057       _RL    etaPre(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0058       _RL   etaZPre(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0059       _RL   dwatPre(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
5acccad966 Jean*0060       INTEGER newtonIter
                0061       INTEGER krylovIter
255a701086 Mart*0062       _RL     myTime
                0063       INTEGER myIter
                0064       INTEGER myThid
5acccad966 Jean*0065
                0066 C     !OUTPUT PARAMETERS:
255a701086 Mart*0067 C     du/vIce :: solution vector
                0068       _RL duIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0069       _RL dvIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
5acccad966 Jean*0070 CEOP
255a701086 Mart*0071
45315406aa Mart*0072 #if ( defined SEAICE_CGRID && \
                0073       ( defined SEAICE_ALLOW_JFNK || defined SEAICE_ALLOW_KRYLOV ) )
255a701086 Mart*0074 C     !FUNCTIONS:
                0075       LOGICAL  DIFFERENT_MULTIPLE
                0076       EXTERNAL DIFFERENT_MULTIPLE
                0077
                0078 C     !LOCAL VARIABLES:
                0079 C     === Local variables ===
                0080 C     i,j,bi,bj  :: Loop counters
                0081
bf019d885c Mart*0082       INTEGER i, j, m, bi, bj
255a701086 Mart*0083       INTEGER k
e780b26e64 Mart*0084       INTEGER iMin, iMax, jMin, jMax
255a701086 Mart*0085       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0086
cea0be0fc5 Mart*0087       _RL WFAU, WFAV
255a701086 Mart*0088
                0089 C     diagonals of coefficient matrices
                0090       _RL AU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0091       _RL BU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0092       _RL CU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0093       _RL AV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0094       _RL BV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0095       _RL CV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0096 C     RHS
                0097       _RL rhsU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0098       _RL rhsV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
cea0be0fc5 Mart*0099       _RL rhsU0(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0100       _RL rhsV0(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
255a701086 Mart*0101 C     coefficients for lateral points, u(j+/-1)
                0102       _RL uRt1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0103       _RL uRt2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0104 C     coefficients for lateral points, v(i+/-1)
                0105       _RL vRt1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0106       _RL vRt2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0107 C     abbreviations
                0108       _RL etaPlusZeta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0109       _RL zetaMinusEta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
df1dac8b7b Mart*0110 C     symmetric drag coefficient
                0111       _RL dragSym(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
255a701086 Mart*0112 C     auxillary fields
                0113       _RL uTmp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0114       _RL vTmp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
5acccad966 Jean*0115       _RS SINWAT
255a701086 Mart*0116       _RL COSWAT
cea0be0fc5 Mart*0117       _RL coriFac
                0118       _RL fricFac
255a701086 Mart*0119       LOGICAL printResidual
                0120       _RL residUini, residVini, residUend, residVend
e780b26e64 Mart*0121 C
255a701086 Mart*0122 CEOP
                0123
                0124 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0125
093db464af Mart*0126       printResidual = debugLevel.GE.debLevC
255a701086 Mart*0127      &  .AND. DIFFERENT_MULTIPLE( SEAICE_monFreq, myTime, deltaTClock )
                0128
e780b26e64 Mart*0129 C     extra overlap for (restricted) additive Schwarz method
                0130       jMin = 1-SEAICE_OLy
                0131       jMax = sNy+SEAICE_OLy
                0132       iMin = 1-SEAICE_OLx
                0133       iMax = sNx+SEAICE_OLx
cea0be0fc5 Mart*0134 C     convergence is affected with coriFac = fricFac = 1
                0135       coriFac = 0. _d 0
                0136       fricFac = coriFac
255a701086 Mart*0137 C     surface level
                0138       k = 1
                0139 C--   introduce turning angles
                0140       SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)
                0141       COSWAT=COS(SEAICE_waterTurnAngle*deg2rad)
                0142
e45202e340 Mart*0143 C     copy relaxation parameters
                0144       WFAU=SEAICE_LSRrelaxU
                0145       WFAV=SEAICE_LSRrelaxV
3f31c7d5de Mart*0146 C
cea0be0fc5 Mart*0147 C     Initialise
3f31c7d5de Mart*0148 C
255a701086 Mart*0149       DO bj=myByLo(myThid),myByHi(myThid)
                0150        DO bi=myBxLo(myThid),myBxHi(myThid)
                0151         DO j=1-OLy,sNy+OLy
                0152          DO i=1-OLx,sNx+OLx
cea0be0fc5 Mart*0153           rhsU (I,J,bi,bj) = 0. _d 0
                0154           rhsV (I,J,bi,bj) = 0. _d 0
                0155           rhsU0(I,J,bi,bj) = duIce(I,J,bi,bj)
                0156           rhsV0(I,J,bi,bj) = dvIce(I,J,bi,bj)
255a701086 Mart*0157 C     first guess for the increment is 0.
cea0be0fc5 Mart*0158           duIce(I,J,bi,bj) = 0. _d 0
                0159           dvIce(I,J,bi,bj) = 0. _d 0
df1dac8b7b Mart*0160 C     this is only the symmetric part of the drag
                0161           dragSym(I,J,bi,bj) = dwatPre(I,J,bi,bj)*COSWAT
255a701086 Mart*0162          ENDDO
                0163         ENDDO
                0164        ENDDO
                0165       ENDDO
                0166 C
                0167 C     some abbreviations
                0168 C
                0169       DO bj=myByLo(myThid),myByHi(myThid)
                0170        DO bi=myBxLo(myThid),myBxHi(myThid)
e780b26e64 Mart*0171         DO J=jMin-1,jMax
                0172          DO I=iMin-1,iMax
255a701086 Mart*0173           etaPlusZeta (I,J,bi,bj)= etaPre(I,J,bi,bj)+zetaPre(I,J,bi,bj)
                0174           zetaMinusEta(I,J,bi,bj)=zetaPre(I,J,bi,bj)- etaPre(I,J,bi,bj)
                0175          ENDDO
                0176         ENDDO
                0177        ENDDO
                0178       ENDDO
438e56056c Mart*0179 C
                0180 C     calculate coefficients of tridiagonal matrices for both u- and
                0181 C     v-equations
                0182 C
76afb58b68 Mart*0183       CALL SEAICE_LSR_CALC_COEFFS(
df1dac8b7b Mart*0184      I     etaPlusZeta, zetaMinusEta, etaZpre, zetaZpre, dragSym,
438e56056c Mart*0185      O     AU, BU, CU, AV, BV, CV, uRt1, uRt2, vRt1, vRt2,
e780b26e64 Mart*0186      I     iMin, iMax, jMin, jMax, myTime, myIter, myThid )
255a701086 Mart*0187
438e56056c Mart*0188 #ifndef OBCS_UVICE_OLD
                0189 C--     prevent tri-diagonal solver from modifying OB values:
255a701086 Mart*0190       DO bj=myByLo(myThid),myByHi(myThid)
                0191        DO bi=myBxLo(myThid),myBxHi(myThid)
e780b26e64 Mart*0192         DO J=jMin,jMax
                0193          DO I=iMin,iMax
255a701086 Mart*0194           IF ( maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj) .EQ. 0. ) THEN
438e56056c Mart*0195            AU(I,J,bi,bj)   = ZERO
                0196            BU(I,J,bi,bj)   = ONE
                0197            CU(I,J,bi,bj)   = ZERO
                0198            uRt1(I,J,bi,bj) = ZERO
                0199            uRt2(I,J,bi,bj) = ZERO
255a701086 Mart*0200           ENDIF
                0201           IF ( maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj) .EQ. 0. ) THEN
438e56056c Mart*0202            AV(I,J,bi,bj)   = ZERO
                0203            BV(I,J,bi,bj)   = ONE
                0204            CV(I,J,bi,bj)   = ZERO
20f7ac573b Mart*0205            vRt1(I,J,bi,bj) = ZERO
                0206            vRt2(I,J,bi,bj) = ZERO
255a701086 Mart*0207           ENDIF
                0208          ENDDO
                0209         ENDDO
                0210        ENDDO
                0211       ENDDO
438e56056c Mart*0212 #endif /* OBCS_UVICE_OLD */
255a701086 Mart*0213
                0214 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0215
                0216 #ifdef ALLOW_DEBUG
                0217       IF ( debugLevel .GE. debLevD ) THEN
                0218         WRITE(msgBuf,'(A,I3,A,I3,A)')
5acccad966 Jean*0219      &        'Uice pre iter (SEAICE_PRECONDITIONER',
255a701086 Mart*0220      &      newtonIter, ',', krylovIter, ')'
                0221         CALL DEBUG_STATS_RL( 1, UICE, msgBuf, myThid )
                0222         WRITE(msgBuf,'(A,I3,A,I3,A)')
5acccad966 Jean*0223      &        'Vice pre iter (SEAICE_PRECONDITIONER',
255a701086 Mart*0224      &      newtonIter, ',', krylovIter, ')'
                0225         CALL DEBUG_STATS_RL( 1, VICE, msgBuf, myThid )
                0226       ENDIF
                0227 #endif /* ALLOW_DEBUG */
                0228
                0229 C--   Calculate initial residual of the linearised system
304ceb7c08 Mart*0230       IF ( printResidual ) THEN
20f7ac573b Mart*0231 C     set up right-hand side now (will be redone in each iteration)
                0232        DO bj=myByLo(myThid),myByHi(myThid)
                0233         DO bi=myBxLo(myThid),myBxHi(myThid)
e780b26e64 Mart*0234          DO j=jMin,jMax
                0235           DO i=iMin,iMax
20f7ac573b Mart*0236            rhsU(I,J,bi,bj) = rhsU0(I,J,bi,bj)
                0237            rhsV(I,J,bi,bj) = rhsV0(I,J,bi,bj)
                0238           ENDDO
                0239          ENDDO
e780b26e64 Mart*0240          CALL SEAICE_PRECOND_RHSU (
48db1b87b2 Mart*0241      I        zetaMinusEta, etaPlusZeta, etaZpre, zetaZpre,
20f7ac573b Mart*0242      I        dwatPre, coriFac, fricFac, SINWAT, COSWAT,
48db1b87b2 Mart*0243      I        duIce, dvIce,
20f7ac573b Mart*0244      O        rhsU,
e780b26e64 Mart*0245      I        iMin,iMax,jMin,jMax,bi,bj,myThid )
                0246          CALL SEAICE_PRECOND_RHSV (
48db1b87b2 Mart*0247      I        zetaMinusEta, etaPlusZeta, etaZpre, zetaZpre,
20f7ac573b Mart*0248      I        dwatPre, coriFac, fricFac, SINWAT, COSWAT,
48db1b87b2 Mart*0249      I        duIce, dvIce,
20f7ac573b Mart*0250      O        rhsV,
e780b26e64 Mart*0251      I        iMin,iMax,jMin,jMax,bi,bj,myThid )
20f7ac573b Mart*0252 #ifndef OBCS_UVICE_OLD
e780b26e64 Mart*0253          DO J=jMin,jMax
                0254           DO I=iMin,iMax
20f7ac573b Mart*0255            IF ( maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj) .EQ. 0. ) THEN
                0256             rhsU(I,J,bi,bj) = duIce(I,J,bi,bj)
                0257            ENDIF
1d7787dd60 Mart*0258            IF ( maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj) .EQ. 0. ) THEN
                0259             rhsV(I,J,bi,bj) = dvIce(I,J,bi,bj)
                0260            ENDIF
20f7ac573b Mart*0261           ENDDO
                0262          ENDDO
                0263 #endif /* OBCS_UVICE_OLD */
                0264         ENDDO
                0265        ENDDO
                0266        CALL SEAICE_RESIDUAL(
255a701086 Mart*0267      I                  rhsU, rhsV, uRt1, uRt2, vRt1, vRt2,
304ceb7c08 Mart*0268      I                  AU, BU, CU, AV, BV, CV, duIce, dvIce,
255a701086 Mart*0269      O                  residUini, residVini, uTmp, vTmp,
                0270      I                  printResidual, myIter, myThid )
                0271       ENDIF
                0272
                0273 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0274
                0275 C NOW DO ITERATION
                0276
                0277 C ITERATION START -----------------------------------------------------
                0278
79df32c3f1 Mart*0279       DO m = 1, SEAICEpreconLinIter
255a701086 Mart*0280
                0281        DO bj=myByLo(myThid),myByHi(myThid)
                0282         DO bi=myBxLo(myThid),myBxHi(myThid)
5acccad966 Jean*0283
8df0ccd026 Jean*0284 C     save du/vIce prior to iteration
                0285          DO j=1-OLy,sNy+OLy
                0286           DO i=1-OLx,sNx+OLx
                0287            uTmp(I,J,bi,bj)=duIce(I,J,bi,bj)
                0288            vTmp(I,J,bi,bj)=dvIce(I,J,bi,bj)
                0289           ENDDO
                0290          ENDDO
                0291
20f7ac573b Mart*0292 C     set up right-hand sides for u- and v-equations
e780b26e64 Mart*0293          DO j=jMin,jMax
                0294           DO i=iMin,iMax
cea0be0fc5 Mart*0295            rhsU(I,J,bi,bj) = rhsU0(I,J,bi,bj)
e780b26e64 Mart*0296 #ifndef SEAICE_PRECOND_EXTRA_EXCHANGE
20f7ac573b Mart*0297            rhsV(I,J,bi,bj) = rhsV0(I,J,bi,bj)
e780b26e64 Mart*0298 #endif /* SEAICE_PRECOND_EXTRA_EXCHANGE */
cea0be0fc5 Mart*0299           ENDDO
                0300          ENDDO
e780b26e64 Mart*0301          CALL SEAICE_PRECOND_RHSU (
48db1b87b2 Mart*0302      I        zetaMinusEta, etaPlusZeta, etaZpre, zetaZPre,
cea0be0fc5 Mart*0303      I        dwatPre, coriFac, fricFac, SINWAT, COSWAT,
48db1b87b2 Mart*0304      I        duIce, dvIce,
cea0be0fc5 Mart*0305      U        rhsU,
e780b26e64 Mart*0306      I        iMin,iMax,jMin,jMax,bi,bj,myThid )
                0307 #ifndef SEAICE_PRECOND_EXTRA_EXCHANGE
                0308          CALL SEAICE_PRECOND_RHSV (
48db1b87b2 Mart*0309      I        zetaMinusEta, etaPlusZeta, etaZpre, zetaZpre,
20f7ac573b Mart*0310      I        dwatPre, coriFac, fricFac, SINWAT, COSWAT,
48db1b87b2 Mart*0311      I        duIce, dvIce,
20f7ac573b Mart*0312      U        rhsV,
e780b26e64 Mart*0313      I        iMin,iMax,jMin,jMax,bi,bj,myThid )
                0314 #endif /* SEAICE_PRECOND_EXTRA_EXCHANGE */
cea0be0fc5 Mart*0315 #ifndef OBCS_UVICE_OLD
48db1b87b2 Mart*0316 C--     prevent tri-diagonal solver from modifying OB values:
e780b26e64 Mart*0317          DO J=jMin,jMax
                0318           DO I=iMin,iMax
cea0be0fc5 Mart*0319            IF ( maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj) .EQ. 0. ) THEN
                0320             rhsU(I,J,bi,bj) = duIce(I,J,bi,bj)
                0321            ENDIF
e780b26e64 Mart*0322 #ifndef SEAICE_PRECOND_EXTRA_EXCHANGE
20f7ac573b Mart*0323            IF ( maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj) .EQ. 0. ) THEN
                0324             rhsV(I,J,bi,bj) = dvIce(I,J,bi,bj)
                0325            ENDIF
e780b26e64 Mart*0326 #endif /* SEAICE_PRECOND_EXTRA_EXCHANGE */
cea0be0fc5 Mart*0327           ENDDO
                0328          ENDDO
                0329 #endif /* OBCS_UVICE_OLD */
255a701086 Mart*0330
                0331 C Solve for uIce :
bf019d885c Mart*0332          CALL SEAICE_LSR_TRIDIAGU(
8df0ccd026 Jean*0333      I        AU, BU, CU, uRt1, uRt2, rhsU, uTmp, seaiceMaskU, WFAU,
bf019d885c Mart*0334      U        duIce,
                0335      I        imin, imax, jmin, jmax, bi, bj, myTime, myIter, myThid )
5acccad966 Jean*0336
e780b26e64 Mart*0337 #ifdef SEAICE_PRECOND_EXTRA_EXCHANGE
20f7ac573b Mart*0338         ENDDO
                0339        ENDDO
                0340 C     ideally one would like to get rid off this exchange
                0341        CALL EXCH_UV_XY_RL( duIce, dvIce, .TRUE., myThid )
                0342
                0343        DO bj=myByLo(myThid),myByHi(myThid)
                0344         DO bi=myBxLo(myThid),myBxHi(myThid)
cea0be0fc5 Mart*0345 C     set up right-hand-side for v-equation
e780b26e64 Mart*0346          DO j=jMin,jMax
                0347           DO i=iMin,iMax
cea0be0fc5 Mart*0348            rhsV(I,J,bi,bj) = rhsV0(I,J,bi,bj)
                0349           ENDDO
                0350          ENDDO
e780b26e64 Mart*0351          CALL SEAICE_PRECOND_RHSV (
48db1b87b2 Mart*0352      I        zetaMinusEta, etaPlusZeta, etaZpre, zetaZpre,
cea0be0fc5 Mart*0353      I        dwatPre, coriFac, fricFac, SINWAT, COSWAT,
48db1b87b2 Mart*0354      I        duIce, dvIce,
cea0be0fc5 Mart*0355      U        rhsV,
e780b26e64 Mart*0356      I        iMin,iMax,jMin,jMax,bi,bj,myThid )
cea0be0fc5 Mart*0357 #ifndef OBCS_UVICE_OLD
48db1b87b2 Mart*0358 C--     prevent tri-diagonal solver from modifying OB values:
e780b26e64 Mart*0359          DO J=jMin,jMax
                0360           DO I=iMin,iMax
cea0be0fc5 Mart*0361            IF ( maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj) .EQ. 0. ) THEN
                0362             rhsV(I,J,bi,bj) = dvIce(I,J,bi,bj)
                0363            ENDIF
                0364           ENDDO
                0365          ENDDO
                0366 #endif /* OBCS_UVICE_OLD */
e780b26e64 Mart*0367 #endif /* SEAICE_PRECOND_EXTRA_EXCHANGE */
cea0be0fc5 Mart*0368
255a701086 Mart*0369 C Solve for dvIce
bf019d885c Mart*0370          CALL SEAICE_LSR_TRIDIAGV(
778d8fff73 Jean*0371      I        AV, BV, CV, vRt1, vRt2, rhsV, vTmp, seaiceMaskV, WFAV,
bf019d885c Mart*0372      U        dvIce,
                0373      I        imin, imax, jmin, jmax, bi, bj, myTime, myIter, myThid )
255a701086 Mart*0374
                0375 C     end bi,bj-loops
                0376         ENDDO
                0377        ENDDO
5acccad966 Jean*0378
255a701086 Mart*0379        CALL EXCH_UV_XY_RL( duIce, dvIce, .TRUE., myThid )
                0380
                0381       ENDDO
                0382 C ITERATION END -----------------------------------------------------
                0383
                0384 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0385
                0386       IF ( printResidual ) THEN
                0387 C--   Calculate final residual of the linearised system
                0388         CALL SEAICE_RESIDUAL(
                0389      I                  rhsU, rhsV, uRt1, uRt2, vRt1, vRt2,
                0390      I                  AU, BU, CU, AV, BV, CV, duIce, dvIce,
                0391      O                  residUend, residVend, uTmp, vTmp,
                0392      I                  printResidual, myIter, myThid )
                0393         _BEGIN_MASTER( myThid )
304ceb7c08 Mart*0394         WRITE(standardMessageUnit,'(A,A,1X,1P2E16.8)')
                0395      &       ' SEAICE_PRECONDITIONER: Residual Initial Uice,Vice     =',
                0396      &       '     ', residUini, residVini
255a701086 Mart*0397         WRITE(standardMessageUnit,'(A,I4,A,I4,A,I6,1P2E16.8)')
                0398      &       ' SEAICE_PRECONDITIONER (iter=',newtonIter,',',
5acccad966 Jean*0399      &       krylovIter, ') iters, U/VResid=',
79df32c3f1 Mart*0400      &       SEAICEpreconLinIter, residUend, residVend
255a701086 Mart*0401         _END_MASTER( myThid )
                0402       ENDIF
                0403 #ifdef ALLOW_DEBUG
                0404       IF ( debugLevel .GE. debLevD ) THEN
                0405         WRITE(msgBuf,'(A,I3,A,I3,A)')
5acccad966 Jean*0406      &        'Uice post iter (SEAICE_PRECONDITIONER',
255a701086 Mart*0407      &      newtonIter, ',', krylovIter, ')'
                0408         CALL DEBUG_STATS_RL( 1, UICE, msgBuf, myThid )
                0409         WRITE(msgBuf,'(A,I3,A,I3,A)')
5acccad966 Jean*0410      &        'Vice post iter (SEAICE_PRECONDITIONER',
255a701086 Mart*0411      &      newtonIter, ',', krylovIter, ')'
                0412         CALL DEBUG_STATS_RL( 1, VICE, msgBuf, myThid )
                0413       ENDIF
                0414 #endif /* ALLOW_DEBUG */
                0415
                0416 C     APPLY MASKS
                0417       DO bj=myByLo(myThid),myByHi(myThid)
                0418        DO bi=myBxLo(myThid),myBxHi(myThid)
                0419         DO J=1-OLy,sNy+OLy
                0420          DO I=1-OLx,sNx+OLx
                0421           duIce(I,J,bi,bj)=duIce(I,J,bi,bj)* seaiceMaskU(I,J,bi,bj)
                0422           dvIce(I,J,bi,bj)=dvIce(I,J,bi,bj)* seaiceMaskV(I,J,bi,bj)
                0423          ENDDO
                0424         ENDDO
                0425        ENDDO
                0426       ENDDO
                0427
5acccad966 Jean*0428       RETURN
cea0be0fc5 Mart*0429       END
                0430
1d7787dd60 Mart*0431 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0432
5acccad966 Jean*0433 CBOP
e780b26e64 Mart*0434 C     !ROUTINE: SEAICE_PRECOND_RHSU
5acccad966 Jean*0435 C     !INTERFACE:
e780b26e64 Mart*0436       SUBROUTINE SEAICE_PRECOND_RHSU (
48db1b87b2 Mart*0437      I     zetaMinusEta, etaPlusZeta, etaZpre, zetaZpre,
cea0be0fc5 Mart*0438      I     dwatPre, coriFac, fricFac, SINWAT, COSWAT,
48db1b87b2 Mart*0439      I     uIceLoc, vIceLoc,
cea0be0fc5 Mart*0440      U     rhsU,
e780b26e64 Mart*0441      I     iMin,iMax,jMin,jMax,bi,bj,myThid )
cea0be0fc5 Mart*0442
1d7787dd60 Mart*0443 C     !DESCRIPTION: \bv
                0444 C     *==========================================================*
e780b26e64 Mart*0445 C     | SUBROUTINE SEAICE_PRECOND_RHSU
5acccad966 Jean*0446 C     | o Calculate the right-hand-side of the u-momentum equation
1d7787dd60 Mart*0447 C     *==========================================================*
                0448 C     \ev
                0449
5acccad966 Jean*0450 C     !USES:
cea0be0fc5 Mart*0451       IMPLICIT NONE
5acccad966 Jean*0452
cea0be0fc5 Mart*0453 #include "SIZE.h"
                0454 #include "EEPARAMS.h"
                0455 #include "PARAMS.h"
                0456 #include "DYNVARS.h"
                0457 #include "GRID.h"
                0458 #include "SEAICE_SIZE.h"
48db1b87b2 Mart*0459 #include "SEAICE_PARAMS.h"
cea0be0fc5 Mart*0460 #include "SEAICE.h"
                0461
5acccad966 Jean*0462 C     !INPUT/OUTPUT PARAMETERS:
cea0be0fc5 Mart*0463       _RL zetaMinusEta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
5acccad966 Jean*0464       _RL etaPlusZeta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
48db1b87b2 Mart*0465       _RL  etaZpre    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0466       _RL zetaZpre    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0467       _RL uIceLoc     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
5acccad966 Jean*0468       _RL vIceLoc     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0469       _RL dwatPre     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0470       _RL coriFac, fricFac
                0471       _RS SINWAT
                0472       _RL COSWAT
                0473       _RL rhsU        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
e780b26e64 Mart*0474       INTEGER iMin, iMax, jMin, jMax, bi, bj, myThid
5acccad966 Jean*0475 CEOP
cea0be0fc5 Mart*0476
5acccad966 Jean*0477 C     !LOCAL VARIABLES:
cea0be0fc5 Mart*0478       INTEGER I,J,K
48db1b87b2 Mart*0479       _RL zeros(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
70e078b38a Mart*0480       _RL areaW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
48db1b87b2 Mart*0481
cea0be0fc5 Mart*0482 C     surface level
                0483       k = 1
48db1b87b2 Mart*0484 C     set dummy pressure to zero
                0485       DO J=1-OLy,sNy+OLy
                0486        DO I=1-OLx,sNx+OLx
                0487         zeros(I,J,bi,bj) = 0. _d 0
cea0be0fc5 Mart*0488        ENDDO
                0489       ENDDO
48db1b87b2 Mart*0490       CALL SEAICE_LSR_RHSU(
                0491      I     zetaMinusEta, etaPlusZeta, etaZpre, zetaZpre, zeros,
778d8fff73 Jean*0492      I     uIceLoc, vIceLoc,
48db1b87b2 Mart*0493      U     rhsU,
                0494      I     iMin, iMax, jMin, jMax, bi, bj, myThid )
cea0be0fc5 Mart*0495
                0496 C     neglected for preconditioning step
965029703f Mart*0497       IF ( fricFac+coriFac .NE. 0. _d 0 ) THEN
70e078b38a Mart*0498        IF ( SEAICEscaleSurfStress ) THEN
                0499         DO J=jMin,jMax
                0500          DO I=iMin,iMax
                0501           areaW(I,J) = 0.5 _d 0*(AREA(I,J,bi,bj)+AREA(I-1,J,bi,bj))
                0502          ENDDO
                0503         ENDDO
                0504        ELSE
                0505         DO J=jMin,jMax
                0506          DO I=iMin,iMax
                0507           areaW(I,J) = 1. _d 0
                0508          ENDDO
                0509         ENDDO
                0510        ENDIF
e780b26e64 Mart*0511        DO J=jMin,jMax
                0512         DO I=iMin,iMax
cea0be0fc5 Mart*0513          rhsU(I,J,bi,bj) = rhsU(I,J,bi,bj)
                0514      &        - SIGN(SINWAT, _fCori(I,J,bi,bj))* 0.5 _d 0 *
                0515      &        ( dwatPre(I  ,J,bi,bj) * 0.5 _d 0 *
                0516      &        (vVel(I  ,J  ,k,bi,bj)-vIceLoc(I  ,J  ,bi,bj)
                0517      &        +vVel(I  ,J+1,k,bi,bj)-vIceLoc(I  ,J+1,bi,bj))
                0518      &        + dwatPre(I-1,J,bi,bj) * 0.5 _d 0 *
                0519      &        (vVel(I-1,J  ,k,bi,bj)-vIceLoc(I-1,J  ,bi,bj)
                0520      &        +vVel(I-1,J+1,k,bi,bj)-vIceLoc(I-1,J+1,bi,bj))
70e078b38a Mart*0521      &        ) * fricFac * areaW(I,J)
cea0be0fc5 Mart*0522 C-    add Coriolis term
                0523          rhsU(I,J,bi,bj) = rhsU(I,J,bi,bj) + 0.5 _d 0 *
                0524      &        ( seaiceMassC(I  ,J,bi,bj) * _fCori(I  ,J,bi,bj)
                0525      &        *0.5 _d 0*(vIceLoc( i ,j,bi,bj)+vIceLoc( i ,j+1,bi,bj))
                0526      &        + seaiceMassC(I-1,J,bi,bj) * _fCori(I-1,J,bi,bj)
5acccad966 Jean*0527      &        *0.5 _d 0*(vIceLoc(i-1,j,bi,bj)+vIceLoc(i-1,j+1,bi,bj))
cea0be0fc5 Mart*0528      &        ) * coriFac
                0529         ENDDO
                0530        ENDDO
                0531       ENDIF
5acccad966 Jean*0532
                0533       RETURN
cea0be0fc5 Mart*0534       END
                0535
1d7787dd60 Mart*0536 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0537
5acccad966 Jean*0538 CBOP
e780b26e64 Mart*0539 C     !ROUTINE: SEAICE_PRECOND_RHSV
5acccad966 Jean*0540 C     !INTERFACE:
e780b26e64 Mart*0541       SUBROUTINE SEAICE_PRECOND_RHSV (
48db1b87b2 Mart*0542      I     zetaMinusEta, etaPlusZeta, etaZpre, zetaZpre,
cea0be0fc5 Mart*0543      I     dwatPre, coriFac, fricFac, SINWAT, COSWAT,
48db1b87b2 Mart*0544      I     uIceLoc, vIceLoc,
cea0be0fc5 Mart*0545      U     rhsV,
e780b26e64 Mart*0546      I     iMin,iMax,jMin,jMax,bi,bj,myThid )
cea0be0fc5 Mart*0547
1d7787dd60 Mart*0548 C     !DESCRIPTION: \bv
                0549 C     *==========================================================*
e780b26e64 Mart*0550 C     | SUBROUTINE SEAICE_PRECOND_RHSV
1d7787dd60 Mart*0551 C     | o Calculate the right-hand-side of the v-momentum equation
                0552 C     *==========================================================*
                0553 C     \ev
                0554
5acccad966 Jean*0555 C     !USES:
cea0be0fc5 Mart*0556       IMPLICIT NONE
5acccad966 Jean*0557
cea0be0fc5 Mart*0558 #include "SIZE.h"
                0559 #include "EEPARAMS.h"
                0560 #include "PARAMS.h"
                0561 #include "DYNVARS.h"
                0562 #include "GRID.h"
                0563 #include "SEAICE_SIZE.h"
48db1b87b2 Mart*0564 #include "SEAICE_PARAMS.h"
cea0be0fc5 Mart*0565 #include "SEAICE.h"
                0566
5acccad966 Jean*0567 C     !INPUT/OUTPUT PARAMETERS:
                0568       _RL zetaMinusEta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0569       _RL etaPlusZeta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
48db1b87b2 Mart*0570       _RL  etaZpre    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0571       _RL zetaZpre    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
5acccad966 Jean*0572       _RL uIceLoc     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
48db1b87b2 Mart*0573       _RL vIceLoc     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
5acccad966 Jean*0574       _RL dwatPre     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
cea0be0fc5 Mart*0575       _RL coriFac, fricFac
5acccad966 Jean*0576       _RS SINWAT
cea0be0fc5 Mart*0577       _RL COSWAT
                0578       _RL rhsV        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
e780b26e64 Mart*0579       INTEGER iMin, iMax, jMin, jMax, bi, bj, myThid
5acccad966 Jean*0580 CEOP
cea0be0fc5 Mart*0581
5acccad966 Jean*0582 C     !LOCAL VARIABLES:
cea0be0fc5 Mart*0583       INTEGER I,J,K
48db1b87b2 Mart*0584       _RL zeros(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
70e078b38a Mart*0585       _RL areaS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
48db1b87b2 Mart*0586
cea0be0fc5 Mart*0587 C     surface level
                0588       k = 1
48db1b87b2 Mart*0589 C     set dummy pressure to zero
                0590       DO J=1-OLy,sNy+OLy
                0591        DO I=1-OLx,sNx+OLx
                0592         zeros(I,J,bi,bj) = 0. _d 0
cea0be0fc5 Mart*0593        ENDDO
                0594       ENDDO
48db1b87b2 Mart*0595       CALL SEAICE_LSR_RHSV(
                0596      I     zetaMinusEta, etaPlusZeta, etaZpre, zetaZpre, zeros,
778d8fff73 Jean*0597      I     uIceLoc, vIceLoc,
48db1b87b2 Mart*0598      U     rhsV,
                0599      I     iMin, iMax, jMin, jMax, bi, bj, myThid )
cea0be0fc5 Mart*0600
                0601 C     neglected for preconditioning step
                0602       IF ( fricFac+coriFac .NE. 0. _d 0 ) THEN
70e078b38a Mart*0603        IF ( SEAICEscaleSurfStress ) THEN
                0604         DO J=jMin,jMax
                0605          DO I=iMin,iMax
0e072cb3ca Mart*0606           areaS(I,J) = 0.5 _d 0*(AREA(I,J,bi,bj)+AREA(I,J-1,bi,bj))
70e078b38a Mart*0607          ENDDO
                0608         ENDDO
                0609        ELSE
                0610         DO J=jMin,jMax
                0611          DO I=iMin,iMax
                0612           areaS(I,J) = 1. _d 0
                0613          ENDDO
                0614         ENDDO
                0615        ENDIF
e780b26e64 Mart*0616        DO J=jMin,jMax
                0617         DO I=iMin,iMax
5acccad966 Jean*0618          rhsV(I,J,bi,bj) = rhsV(I,J,bi,bj)
cea0be0fc5 Mart*0619      &        + SIGN(SINWAT, _fCori(I,J,bi,bj)) * 0.5 _d 0 *
                0620      &        ( dwatPre(I,J  ,bi,bj) * 0.5 _d 0 *
                0621      &        (uVel(I  ,J  ,k,bi,bj)-uIceLoc(I  ,J  ,bi,bj)
                0622      &        +uVel(I+1,J  ,k,bi,bj)-uIceLoc(I+1,J  ,bi,bj))
                0623      &        + dwatPre(I,J-1,bi,bj) * 0.5 _d 0 *
                0624      &        (uVel(I  ,J-1,k,bi,bj)-uIceLoc(I  ,J-1,bi,bj)
                0625      &        +uVel(I+1,J-1,k,bi,bj)-uIceLoc(I+1,J-1,bi,bj))
70e078b38a Mart*0626      &        ) * fricFac * areaS(I,J)
cea0be0fc5 Mart*0627 C-    add Coriolis term
                0628          rhsV(I,J,bi,bj) = rhsV(I,J,bi,bj) - 0.5 _d 0 *
                0629      &        ( seaiceMassC(I,J  ,bi,bj) * _fCori(I,J  ,bi,bj)
                0630      &        *0.5 _d 0*(uIceLoc(i  ,j  ,bi,bj)+uIceLoc(i+1,  j,bi,bj))
                0631      &        + seaiceMassC(I,J-1,bi,bj) * _fCori(I,J-1,bi,bj)
                0632      &        *0.5 _d 0*(uIceLoc(i  ,j-1,bi,bj)+uIceLoc(i+1,j-1,bi,bj))
                0633      &        ) * coriFac
                0634         ENDDO
                0635        ENDDO
                0636       ENDIF
                0637
45315406aa Mart*0638 #endif /* SEAICE_CGRID, SEAICE_ALLOW_JFNK and KRYLOV */
255a701086 Mart*0639
                0640       RETURN
                0641       END