pkg/seaice/seaice_krylov.F

c1615e0916 Mart*0001 #include "SEAICE_OPTIONS.h"
                0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
                0005
                0006 C--  File seaice_krylov.F: seaice krylov dynamical solver S/R:
                0007
                0008 CBOP
                0009 C     !ROUTINE: SEAICE_KRYLOV
                0010 C     !INTERFACE:
                0011       SUBROUTINE SEAICE_KRYLOV( myTime, myIter, myThid )
                0012
                0013 C     !DESCRIPTION: \bv
                0014 C     *==========================================================*
                0015 C     | SUBROUTINE SEAICE_KRYLOV
                0016 C     | o Picard solver for ice dynamics using a preconditioned
ec0d7df165 Mart*0017 C     |   KRYLOV (Generalized Minimum RESidual=GMRES) method for
                0018 C     |   solving the linearised system following J.-F. Lemieux
c1615e0916 Mart*0019 C     |   et al., JGR 113, doi:10.1029/2007JC004680, 2008.
                0020 C     *==========================================================*
                0021 C     | written by Martin Losch, Jan 2016
                0022 C     *==========================================================*
                0023 C     \ev
                0024
                0025 C     !USES:
                0026       IMPLICIT NONE
                0027
                0028 C     === Global variables ===
                0029 #include "SIZE.h"
                0030 #include "EEPARAMS.h"
                0031 #include "PARAMS.h"
                0032 #include "DYNVARS.h"
                0033 #include "GRID.h"
                0034 #include "SEAICE_SIZE.h"
                0035 #include "SEAICE_PARAMS.h"
3f0f10fc37 Mart*0036 #include "SEAICE_GRID.h"
c1615e0916 Mart*0037 #include "SEAICE.h"
                0038
                0039 C     !INPUT/OUTPUT PARAMETERS:
                0040 C     === Routine arguments ===
                0041 C     myTime :: Simulation time
                0042 C     myIter :: Simulation timestep number
                0043 C     myThid :: my Thread Id. number
                0044       _RL     myTime
                0045       INTEGER myIter
                0046       INTEGER myThid
                0047
45315406aa Mart*0048 #if ( defined SEAICE_CGRID && defined SEAICE_ALLOW_KRYLOV )
c1615e0916 Mart*0049 C     !FUNCTIONS:
                0050       LOGICAL  DIFFERENT_MULTIPLE
                0051       EXTERNAL DIFFERENT_MULTIPLE
                0052
                0053 C     !LOCAL VARIABLES:
                0054 C     === Local variables ===
                0055 C     i,j,bi,bj :: loop indices
                0056       INTEGER i,j,bi,bj
                0057 C     loop indices
                0058       INTEGER picardIter
                0059       INTEGER krylovIter, krylovFails
                0060       INTEGER krylovIterMax, picardIterMax
                0061       INTEGER totalKrylovItersLoc, totalPicardItersLoc
                0062 C     FGMRES parameters
                0063 C     im      :: size of Krylov space
                0064 C     ifgmres :: interation counter
                0065       INTEGER im
                0066       PARAMETER ( im = 50 )
                0067       INTEGER ifgmres
                0068 C     FGMRES flag that determines amount of output messages of fgmres
                0069       INTEGER iOutFGMRES
                0070 C     FGMRES flag that indicates what fgmres wants us to do next
                0071       INTEGER iCode
                0072       _RL     picardResidual
                0073       _RL     picardResidualKm1
                0074 C     parameters to compute convergence criterion
                0075       _RL     krylovLinTol
                0076       _RL     FGMRESeps
                0077       _RL     picardTol
                0078 C     backward differences extrapolation factors
                0079       _RL bdfFac, bdfAlpha
                0080 C
                0081       _RL     recip_deltaT
                0082       LOGICAL picardConverged, krylovConverged
                0083       LOGICAL writeNow
                0084       CHARACTER*(MAX_LEN_MBUF) msgBuf
ec0d7df165 Mart*0085
c1615e0916 Mart*0086       _RL uIceLin(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0087       _RL vIceLin(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0088 C     extra time level required for backward difference time stepping
                0089       _RL duIcNm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0090       _RL dvIcNm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0091 C     u/vWork   :: work arrays
                0092       _RL uWork  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0093       _RL vWork  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0094 C     u/vIceLHS :: left hand side of momentum equation (A*x)
                0095       _RL uIceLHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0096       _RL vIceLHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0097 C     u/vIceRHS :: right hand side of momentum equation (b)
                0098       _RL uIceRHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0099       _RL vIceRHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0100 C     helper array
                0101       _RL resTmp (nVec,1,nSx,nSy)
                0102 C     work arrays
                0103       _RL rhs(nVec,nSx,nSy), sol(nVec,nSx,nSy)
                0104       _RL vv(nVec,im+1,nSx,nSy), w(nVec,im,nSx,nSy)
                0105       _RL wk1(nVec,nSx,nSy), wk2(nVec,nSx,nSy)
                0106 CEOP
                0107
                0108 C     Initialise
                0109       picardIter          = 0
                0110       krylovFails         = 0
                0111       totalKrylovItersLoc = 0
                0112       picardConverged     = .FALSE.
                0113       picardTol           = 0. _d 0
                0114       picardResidual      = 0. _d 0
                0115       picardResidualKm1   = 0. _d 0
                0116       FGMRESeps           = 0. _d 0
                0117       recip_deltaT        = 1. _d 0 / SEAICE_deltaTdyn
                0118
                0119       krylovIterMax       = SEAICElinearIterMax
                0120       picardIterMax       = SEAICEnonLinIterMax
                0121       IF ( SEAICEusePicardAsPrecon ) THEN
                0122        krylovIterMax       = SEAICEpreconlinIter
                0123        picardIterMax       = SEAICEpreconNL_Iter
                0124       ENDIF
                0125
                0126       iOutFGMRES=0
                0127 C     with iOutFgmres=1, seaice_fgmres prints the residual at each iteration
ec0d7df165 Mart*0128       IF ( debugLevel.GE.debLevC .AND.
c1615e0916 Mart*0129      &     .NOT.SEAICEusePicardAsPrecon .AND.
                0130      &     DIFFERENT_MULTIPLE( SEAICE_monFreq, myTime, deltaTClock ) )
                0131      &     iOutFGMRES=1
                0132
                0133 C     backward difference extrapolation factors
                0134       bdfFac = 0. _d 0
                0135       IF ( SEAICEuseBDF2 ) THEN
                0136        IF ( myIter.EQ.nIter0 .AND. SEAICEmomStartBDF.EQ.0 ) THEN
                0137         bdfFac = 0. _d 0
                0138        ELSE
                0139         bdfFac = 0.5 _d 0
                0140        ENDIF
                0141       ENDIF
ec0d7df165 Mart*0142       bdfAlpha = 1. _d 0 + bdfFac
c1615e0916 Mart*0143
                0144       DO bj=myByLo(myThid),myByHi(myThid)
                0145        DO bi=myBxLo(myThid),myBxHi(myThid)
                0146         DO J=1-OLy,sNy+OLy
                0147          DO I=1-OLx,sNx+OLx
                0148           uIceLHS(I,J,bi,bj) = 0. _d 0
                0149           vIceLHS(I,J,bi,bj) = 0. _d 0
                0150           uIceRHS(I,J,bi,bj) = 0. _d 0
                0151           vIceRHS(I,J,bi,bj) = 0. _d 0
                0152          ENDDO
                0153         ENDDO
                0154 C     cycle ice velocities
                0155         DO J=1-OLy,sNy+OLy
                0156          DO I=1-OLx,sNx+OLx
ec0d7df165 Mart*0157           duIcNm1(I,J,bi,bj) = uIce(I,J,bi,bj) * bdfAlpha
c1615e0916 Mart*0158      &         + ( uIce(I,J,bi,bj) - uIceNm1(I,J,bi,bj) ) * bdfFac
ec0d7df165 Mart*0159           dvIcNm1(I,J,bi,bj) = vIce(I,J,bi,bj) * bdfAlpha
c1615e0916 Mart*0160      &         + ( vIce(I,J,bi,bj) - vIceNm1(I,J,bi,bj) ) * bdfFac
                0161           uIceNm1(I,J,bi,bj) = uIce(I,J,bi,bj)
                0162           vIceNm1(I,J,bi,bj) = vIce(I,J,bi,bj)
                0163           uIceLin(I,J,bi,bj) = uIce(I,J,bi,bj)
                0164           vIceLin(I,J,bi,bj) = vIce(I,J,bi,bj)
                0165          ENDDO
                0166         ENDDO
                0167 C     Compute things that do no change during the OL iteration:
                0168 C     sea-surface tilt and wind stress:
                0169 C     FORCEX/Y0 - mass*(1.5*u/vIceNm1+0.5*(u/vIceNm1-u/vIceNm2))/deltaT
                0170         DO J=1-OLy,sNy+OLy
                0171          DO I=1-OLx,sNx+OLx
                0172           FORCEX(I,J,bi,bj) = FORCEX0(I,J,bi,bj)
                0173      &         + seaiceMassU(I,J,bi,bj)*duIcNm1(I,J,bi,bj)*recip_deltaT
                0174           FORCEY(I,J,bi,bj) = FORCEY0(I,J,bi,bj)
                0175      &         + seaiceMassV(I,J,bi,bj)*dvIcNm1(I,J,bi,bj)*recip_deltaT
                0176          ENDDO
                0177         ENDDO
                0178 CML        ENDIF
                0179        ENDDO
                0180       ENDDO
                0181 C     Start nonlinear Picard iteration: outer loop iteration
                0182       DO WHILE ( picardIter.LT.picardIterMax .AND.
                0183      &     .NOT.picardConverged )
                0184        picardIter = picardIter + 1
                0185 C     smooth ice velocities in time for computation of
ec0d7df165 Mart*0186 C     the non-linear drag coefficents
c1615e0916 Mart*0187        DO bj=myByLo(myThid),myByHi(myThid)
                0188         DO bi=myBxLo(myThid),myBxHi(myThid)
                0189          DO j=1-OLy,sNy+OLy
                0190           DO i=1-OLx,sNx+OLx
                0191            uIceLin(I,J,bi,bj) = 0.5 _d 0 *
                0192      &          (uIce(I,J,bi,bj) + uIceLin(I,J,bi,bj))
ec0d7df165 Mart*0193            vIceLin(I,J,bi,bj) = 0.5 _d 0 *
c1615e0916 Mart*0194      &          (vIce(I,J,bi,bj) + vIceLin(I,J,bi,bj))
                0195           ENDDO
                0196          ENDDO
                0197         ENDDO
                0198        ENDDO
ec0d7df165 Mart*0199 C     u/vIce have changed in Picard iteration so that new drag
c1615e0916 Mart*0200 C     coefficients and viscosities are required (that will not change in
                0201 C     the Krylov iteration)
                0202        CALL SEAICE_OCEANDRAG_COEFFS(
ec0d7df165 Mart*0203      I      uIceLin, vIceLin, HEFFM,
c1615e0916 Mart*0204      O      DWATN,
                0205      I      0, myTime, myIter, myThid )
df1dac8b7b Mart*0206 #ifdef SEAICE_ALLOW_BOTTOMDRAG
d5254b4e3d Mart*0207        CALL SEAICE_BOTTOMDRAG_COEFFS(
ec0d7df165 Mart*0208      I      uIceLin, vIceLin, HEFFM,
df1dac8b7b Mart*0209 #ifdef SEAICE_ITD
                0210      I      HEFFITD, AREAITD, AREA,
                0211 #else
                0212      I      HEFF, AREA,
ec0d7df165 Mart*0213 #endif
df1dac8b7b Mart*0214      O      CbotC,
                0215      I      0, myTime, myIter, myThid )
                0216 #endif /* SEAICE_ALLOW_BOTTOMDRAG */
c1615e0916 Mart*0217        CALL SEAICE_CALC_STRAINRATES(
                0218      I      uIceLin, vIceLin,
                0219      O      e11, e22, e12,
                0220      I      0, myTime, myIter, myThid )
                0221        CALL SEAICE_CALC_VISCOSITIES(
8e32c48b8f Mart*0222      I      e11, e22, e12, SEAICE_zMin, SEAICE_zMax, HEFFM, press0,
                0223      I      tensileStrFac,
c1615e0916 Mart*0224      O      eta, etaZ, zeta, zetaZ, press, deltaC,
                0225      I      0, myTime, myIter, myThid )
                0226 C     compute rhs that does not change during Krylov iteration
                0227        CALL SEAICE_CALC_RHS(
                0228      O      uIceRHS, vIceRHS,
                0229      I      picardIter, 0, myTime, myIter, myThid )
                0230 C     compute rhs for initial residual
                0231        CALL SEAICE_CALC_LHS(
                0232      I         uIce, vIce,
                0233      O         uIceLHS, vIceLHS,
                0234      I         picardIter, myTime, myIter, myThid )
                0235 C     Calculate the residual
                0236        DO bj=myByLo(myThid),myByHi(myThid)
                0237         DO bi=myBxLo(myThid),myBxHi(myThid)
                0238          DO J=1,sNy
                0239           DO I=1,sNx
                0240            uIceLHS(I,J,bi,bj) = uIceLHS(I,J,bi,bj) - uIceRHS(I,J,bi,bj)
                0241            vIceLHS(I,J,bi,bj) = vIceLHS(I,J,bi,bj) - vIceRHS(I,J,bi,bj)
                0242 C     save u/vIceLin as k-2nd step for linearization (does not work properly)
                0243 CML           uIceLin(I,J,bi,bj) = uIce(I,J,bi,bj)
                0244 CML           vIceLin(I,J,bi,bj) = vIce(I,J,bi,bj)
                0245           ENDDO
                0246          ENDDO
                0247         ENDDO
                0248        ENDDO
                0249 C     Important: Compute the norm of the residual using the same scalar
                0250 C     product that SEAICE_FGMRES does
                0251        CALL SEAICE_MAP2VEC(nVec,uIceLHS,vIceLHS,resTmp,.TRUE.,myThid)
                0252        CALL SEAICE_SCALPROD(nVec,1,1,1,resTmp,resTmp,
                0253      &      picardResidual,myThid)
                0254        picardResidual = SQRT(picardResidual)
                0255
                0256 C     compute convergence criterion for linear preconditioned FGMRES
                0257        krylovLinTol = JFNKgamma_lin_max
                0258 C     the best method is still not clear to me
                0259 C     this is described in Lemieux et al 2008
                0260 CML       IF ( picardIter .EQ. 1 ) krylovLinTol = 1./10.
                0261 CML       IF ( picardIter .EQ. 2 ) krylovLinTol = 1./20.
                0262 CML       IF ( picardIter .EQ. 3 ) krylovLinTol = 1./20.
                0263 CML       IF ( picardIter .EQ. 4 ) krylovLinTol = 1./30.
                0264 CML       IF ( picardIter .EQ. 5 ) krylovLinTol = 1./30.
                0265 CML       IF ( picardIter .EQ. 6 ) krylovLinTol = 1./30.
                0266 CML       IF ( picardIter .EQ. 7 ) krylovLinTol = 1./30.
                0267 CML       IF ( picardIter .EQ. 8 ) krylovLinTol = 1./40.
                0268 CML       IF ( picardIter .EQ. 9 ) krylovLinTol = 1./50.
                0269 CML       IF ( picardIter .GT. 9 ) krylovLinTol = 1./80.
                0270 C     this is used with the JFNK solver, but the Picard-Krylov solver
                0271 C     converges too slowly for this scheme
                0272 CML       IF ( picardIter.GT.1.AND.picardResidual.LT.JFNKres_t ) THEN
                0273 CMLC     Eisenstat and Walker (1996), eq.(2.6)
                0274 CML        krylovLinTol = SEAICE_JFNKphi
                0275 CML     &       *( picardResidual/picardResidualKm1 )**SEAICE_JFNKalpha
                0276 CML        krylovLinTol = min(JFNKgamma_lin_max, krylovLinTol)
                0277 CML        krylovLinTol = max(JFNKgamma_lin_min, krylovLinTol)
                0278 CML       ENDIF
                0279 CML       krylovLinTol = 1. _d -1
                0280 C     save the residual for the next iteration
                0281        picardResidualKm1 = picardResidual
                0282
                0283 C     The Krylov iteration uses FGMRES, the preconditioner is LSOR
                0284 C     for now. The code is adapted from SEAICE_LSR, but heavily stripped
                0285 C     down.
                0286 C     krylovIter is mapped into "its" in seaice_fgmres and is incremented
                0287 C     in that routine
                0288        krylovIter    = 0
                0289        iCode         = 0
                0290
                0291        picardConverged = picardResidual.LT.picardTol
0d02be7d13 Mart*0292      &      .OR.picardResidual.EQ.0. _d 0
c1615e0916 Mart*0293
                0294 C     do Krylov loop only if convergence is not reached
                0295
                0296        IF ( .NOT.picardConverged ) THEN
                0297
                0298 C     start Krylov iteration (FGMRES)
                0299
                0300         krylovConverged = .FALSE.
                0301         FGMRESeps = krylovLinTol * picardResidual
                0302         CALL SEAICE_MAP2VEC(nVec,uIce,vIce,sol,.TRUE.,myThid)
                0303         CALL SEAICE_MAP2VEC(nVec,uIceRHS,vIceRHS,rhs,.TRUE.,myThid)
                0304         DO bj=myByLo(myThid),myByHi(myThid)
                0305          DO bi=myBxLo(myThid),myBxHi(myThid)
                0306           DO j=1-OLy,sNy+OLy
                0307            DO i=1-OLx,sNx+OLx
                0308             uWork(i,j,bi,bj) = 0. _d 0
                0309             vWork(i,j,bi,bj) = 0. _d 0
                0310            ENDDO
                0311           ENDDO
                0312          ENDDO
                0313         ENDDO
                0314         DO WHILE ( .NOT.krylovConverged )
                0315 C     solution vector sol = u/vIce
                0316 C     residual vector (rhs) Fu = u/vIceRHS
                0317 C     output work vectors wk1, -> input work vector wk2
                0318
                0319 C     map results to wk2
                0320          CALL SEAICE_MAP2VEC(nVec,uWork,vWork,wk2,.TRUE.,myThid)
                0321
                0322          CALL SEAICE_FGMRES (nVec,im,rhs,sol,ifgmres,krylovIter,
                0323      U        vv,w,wk1,wk2,
                0324      I        FGMRESeps,krylovIterMax,iOutFGMRES,
                0325      U        iCode,
                0326      I        myThid)
                0327 C
                0328          IF ( iCode .EQ. 0 ) THEN
                0329 C     map sol(ution) vector to u/vIce
                0330           CALL SEAICE_MAP2VEC(nVec,uIce,vIce,sol,.FALSE.,myThid)
                0331           CALL EXCH_UV_XY_RL( uIce, vIce,.TRUE.,myThid)
                0332          ELSE
                0333 C     map work vector to du/vIce to either compute a preconditioner
                0334 C     solution (wk1=rhs) or a matrix times wk1
                0335           CALL SEAICE_MAP2VEC(nVec,uWork,vWork,wk1,.FALSE.,myThid)
                0336           CALL EXCH_UV_XY_RL( uWork, vWork,.TRUE.,myThid)
                0337          ENDIF
                0338
                0339 C     FGMRES returns iCode either asking for an new preconditioned vector
                0340 C     or product of matrix times vector. For iCode = 0, terminate
                0341 C     iteration
                0342          IF (iCode.EQ.1) THEN
                0343 C     Call preconditioner
                0344           IF ( SEAICEpreconLinIter .GT. 0 )
                0345      &         CALL SEAICE_PRECONDITIONER(
                0346      U         uWork, vWork,
                0347      I         zeta, eta, etaZ, zetaZ, dwatn,
                0348      I         picardIter, krylovIter, myTime, myIter, myThid )
                0349          ELSEIF (iCode.GE.2) THEN
                0350 C     Compute lhs of equations (A*x)
                0351           CALL SEAICE_CALC_STRAINRATES(
                0352      I         uWork, vWork,
                0353      O         e11, e22, e12,
                0354      I         krylovIter, myTime, myIter, myThid )
                0355           CALL SEAICE_CALC_LHS(
                0356      I         uWork, vWork,
                0357      O         uIceLHS, vIceLHS,
                0358      I         picardIter, myTime, myIter, myThid )
                0359           DO bj=myByLo(myThid),myByHi(myThid)
                0360            DO bi=myBxLo(myThid),myBxHi(myThid)
                0361             DO j=1-OLy,sNy+OLy
                0362              DO i=1-OLx,sNx+OLx
                0363               uWork(i,j,bi,bj) = uIceLHS(i,j,bi,bj)
                0364               vWork(i,j,bi,bj) = vIceLHS(i,j,bi,bj)
                0365              ENDDO
                0366             ENDDO
                0367            ENDDO
                0368           ENDDO
                0369          ENDIF
                0370          krylovConverged = iCode.EQ.0
                0371 C     End of Krylov iterate
                0372         ENDDO
                0373         totalKrylovItersLoc = totalKrylovItersLoc + krylovIter
                0374 C     some output diagnostics
                0375         IF ( debugLevel.GE.debLevA
                0376      &       .AND. .NOT.SEAICEusePicardAsPrecon ) THEN
                0377          _BEGIN_MASTER( myThid )
                0378          totalPicardItersLoc =
                0379      &        picardIterMax*(myIter-nIter0)+picardIter
871ac80009 Jean*0380          WRITE(msgBuf,'(2A,2(1X,I6),2E12.5)')
c1615e0916 Mart*0381      &        ' S/R SEAICE_KRYLOV: Picard iterate / total, ',
                0382      &        'KRYLOVgamma_lin, initial norm = ',
                0383      &        picardIter, totalPicardItersLoc,
                0384      &        krylovLinTol,picardResidual
                0385          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0386      &        SQUEEZE_RIGHT, myThid )
                0387          WRITE(msgBuf,'(3(A,I6))')
                0388      &        ' S/R SEAICE_KRYLOV: Picard iterate / total = ',
                0389      &        picardIter, ' / ', totalPicardItersLoc,
                0390      &        ', Nb. of FGMRES iterations = ', krylovIter
                0391          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0392      &        SQUEEZE_RIGHT, myThid )
                0393          _END_MASTER( myThid )
                0394         ENDIF
                0395         IF ( krylovIter.EQ.krylovIterMax ) THEN
                0396          krylovFails = krylovFails + 1
                0397         ENDIF
                0398 C     Set the stopping criterion for the Picard iteration and the
                0399 C     criterion for the transition from accurate to approximate FGMRES
                0400         IF ( picardIter .EQ. 1 ) THEN
                0401          picardTol=SEAICEnonLinTol*picardResidual
                0402          IF ( JFNKres_tFac .NE. UNSET_RL )
                0403      &        JFNKres_t = picardResidual * JFNKres_tFac
                0404         ENDIF
                0405        ENDIF
                0406 C     end of Picard iterate
                0407       ENDDO
                0408
                0409 C--   Output diagnostics
                0410
                0411       IF ( SEAICE_monFreq .GT. 0. _d 0 ) THEN
c512e371cc drin*0412 C     Only Master Thread updates counters in common block:
                0413       _BEGIN_MASTER(myThid)
                0414 C     Cumulate some diagnostic counters for the Krylov solver
c1615e0916 Mart*0415        totalJFNKtimeSteps = totalJFNKtimeSteps + 1
                0416        totalNewtonIters   = totalNewtonIters + picardIter
                0417        totalKrylovIters   = totalKrylovIters + totalKrylovItersLoc
                0418 C     Record failure
                0419        totalKrylovFails   = totalKrylovFails + krylovFails
                0420        IF ( picardIter .EQ. picardIterMax ) THEN
                0421         totalNewtonfails = totalNewtonfails + 1
                0422        ENDIF
c512e371cc drin*0423        _END_MASTER( myThid )
c1615e0916 Mart*0424       ENDIF
                0425 C     Decide whether it is time to dump and reset the counter
                0426       writeNow = DIFFERENT_MULTIPLE(SEAICE_monFreq,
                0427      &     myTime+deltaTClock, deltaTClock)
                0428 #ifdef ALLOW_CAL
                0429       IF ( useCAL ) THEN
                0430        CALL CAL_TIME2DUMP(
                0431      I      zeroRL, SEAICE_monFreq,  deltaTClock,
                0432      U      writeNow,
ec0d7df165 Mart*0433      I      myTime+deltaTClock, myIter+1, myThid )
c1615e0916 Mart*0434       ENDIF
                0435 #endif
                0436       IF ( writeNow ) THEN
                0437        _BEGIN_MASTER( myThid )
                0438        WRITE(msgBuf,'(A)')
                0439      &' // ======================================================='
                0440        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0441      &      SQUEEZE_RIGHT, myThid )
                0442        WRITE(msgBuf,'(A)') ' // Begin KRYLOV statistics'
                0443        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0444      &      SQUEEZE_RIGHT, myThid )
                0445        WRITE(msgBuf,'(A)')
                0446      &' // ======================================================='
                0447        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0448      &      SQUEEZE_RIGHT, myThid )
                0449        WRITE(msgBuf,'(A,I10)')
                0450      &      ' %KRYLOV_MON: time step              = ', myIter+1
                0451        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0452      &      SQUEEZE_RIGHT, myThid )
                0453        WRITE(msgBuf,'(A,I10)')
ec0d7df165 Mart*0454      &      ' %KRYLOV_MON: Nb. of time steps      = ',
c1615e0916 Mart*0455      &      totalJFNKtimeSteps
                0456        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0457      &      SQUEEZE_RIGHT, myThid )
                0458        WRITE(msgBuf,'(A,I10)')
                0459      &      ' %KRYLOV_MON: Nb. of Picard steps    = ', totalNewtonIters
                0460        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0461      &      SQUEEZE_RIGHT, myThid )
                0462        WRITE(msgBuf,'(A,I10)')
                0463      &      ' %KRYLOV_MON: Nb. of Krylov steps    = ', totalKrylovIters
                0464        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0465      &      SQUEEZE_RIGHT, myThid )
                0466        WRITE(msgBuf,'(A,I10)')
                0467      &      ' %KRYLOV_MON: Nb. of Picard failures = ', totalNewtonfails
                0468        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0469      &      SQUEEZE_RIGHT, myThid )
                0470        WRITE(msgBuf,'(A,I10)')
                0471      &      ' %KRYLOV_MON: Nb. of Krylov failures = ', totalKrylovFails
                0472        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0473      &      SQUEEZE_RIGHT, myThid )
                0474        WRITE(msgBuf,'(A)')
                0475      &' // ======================================================='
                0476        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0477      &      SQUEEZE_RIGHT, myThid )
                0478        WRITE(msgBuf,'(A)') ' // End KRYLOV statistics'
                0479        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0480      &      SQUEEZE_RIGHT, myThid )
                0481        WRITE(msgBuf,'(A)')
                0482      &' // ======================================================='
                0483        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0484      &      SQUEEZE_RIGHT, myThid )
c512e371cc drin*0485 C     Reset and start again
c1615e0916 Mart*0486        totalJFNKtimeSteps = 0
                0487        totalNewtonIters   = 0
                0488        totalKrylovIters   = 0
                0489        totalKrylovFails   = 0
                0490        totalNewtonfails   = 0
c512e371cc drin*0491        _END_MASTER( myThid )
c1615e0916 Mart*0492       ENDIF
                0493
                0494 C     Print more debugging information
ec0d7df165 Mart*0495       IF ( debugLevel.GE.debLevA
c1615e0916 Mart*0496      &     .AND. .NOT.SEAICEusePicardAsPrecon ) THEN
                0497        IF ( picardIter .EQ. picardIterMax ) THEN
                0498         _BEGIN_MASTER( myThid )
                0499         WRITE(msgBuf,'(A,I10)')
                0500      &       ' S/R SEAICE_KRYLOV: Solver did not converge in timestep ',
                0501      &       myIter+1
                0502         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0503      &       SQUEEZE_RIGHT, myThid )
                0504         _END_MASTER( myThid )
                0505        ENDIF
                0506        IF ( krylovFails .GT. 0 ) THEN
                0507         _BEGIN_MASTER( myThid )
                0508         WRITE(msgBuf,'(A,I4,A,I10)')
                0509      &       ' S/R SEAICE_KRYLOV: FGMRES did not converge ',
                0510      &       krylovFails, ' times in timestep ', myIter+1
                0511         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0512      &       SQUEEZE_RIGHT, myThid )
                0513         _END_MASTER( myThid )
                0514        ENDIF
                0515        _BEGIN_MASTER( myThid )
                0516        WRITE(msgBuf,'(A,I6,A,I10)')
                0517      &      ' S/R SEAICE_KRYLOV: Total number FGMRES iterations = ',
                0518      &      totalKrylovItersLoc, ' in timestep ', myIter+1
                0519        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0520      &      SQUEEZE_RIGHT, myThid )
                0521        _END_MASTER( myThid )
                0522       ENDIF
                0523
45315406aa Mart*0524 #endif /* SEAICE_CGRID and SEAICE_ALLOW_KRYLOV */
c1615e0916 Mart*0525
                0526       RETURN
                0527       END