Back to home page

MITgcm
 
 

    

File indexing completed on 2023-08-04 05:10:47 UTC

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

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!