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