--- MITgcm/pkg/seaice/seaice_jfnk.F 2013/04/23 08:40:06 1.22 +++ MITgcm/pkg/seaice/seaice_jfnk.F 2016/01/27 14:03:34 1.29 @@ -1,7 +1,10 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_jfnk.F,v 1.22 2013/04/23 08:40:06 mlosch Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_jfnk.F,v 1.29 2016/01/27 14:03:34 mlosch Exp $ C $Name: $ #include "SEAICE_OPTIONS.h" +#ifdef ALLOW_AUTODIFF +# include "AUTODIFF_OPTIONS.h" +#endif C-- File seaice_jfnk.F: seaice jfnk dynamical solver S/R: C-- Contents @@ -60,12 +63,18 @@ C !LOCAL VARIABLES: C === Local variables === -C i,j,bi,bj :: loop indices - INTEGER i,j,bi,bj +C i,j,k,bi,bj :: loop indices + INTEGER i,j,k,bi,bj C loop indices INTEGER newtonIter INTEGER krylovIter, krylovFails INTEGER totalKrylovItersLoc, totalNewtonItersLoc +C FGMRES parameters +C im :: size of Krylov space +C ifgmres :: interation counter + INTEGER im + PARAMETER ( im = 50 ) + INTEGER ifgmres C FGMRES flag that determines amount of output messages of fgmres INTEGER iOutFGMRES C FGMRES flag that indicates what fgmres wants us to do next @@ -76,7 +85,9 @@ _RL JFNKgamma_lin _RL FGMRESeps _RL JFNKtol - +C backward differences extrapolation factors + _RL bdfFac, bdfAlpha +C _RL recip_deltaT LOGICAL JFNKconverged, krylovConverged LOGICAL writeNow @@ -85,15 +96,23 @@ C u/vIceRes :: residual of sea-ice momentum equations _RL uIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL vIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) +C extra time level required for backward difference time stepping + _RL duIcNm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL dvIcNm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C du/vIce :: ice velocity increment to be added to u/vIce _RL duIce (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL dvIce (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C precomputed (= constant per Newton iteration) versions of C zeta, eta, and DWATN, press _RL zetaPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) + _RL zetaZPre(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL etaPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL etaZPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL dwatPre (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) +C work arrays + _RL rhs(nVec,nSx,nSy), sol(nVec,nSx,nSy) + _RL vv(nVec,im+1,nSx,nSy), w(nVec,im,nSx,nSy) + _RL wk1(nVec,nSx,nSy), wk2(nVec,nSx,nSy) CEOP C Initialise @@ -113,6 +132,17 @@ & DIFFERENT_MULTIPLE( SEAICE_monFreq, myTime, deltaTClock ) ) & iOutFGMRES=1 +C backward difference extrapolation factors + bdfFac = 0. _d 0 + IF ( SEAICEuseBDF2 ) THEN + IF ( myIter.EQ.nIter0 .AND. SEAICEmomStartBDF.EQ.0 ) THEN + bdfFac = 0. _d 0 + ELSE + bdfFac = 0.5 _d 0 + ENDIF + ENDIF + bdfAlpha = 1. _d 0 + bdfFac + DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO J=1-OLy,sNy+OLy @@ -121,21 +151,33 @@ vIceRes(I,J,bi,bj) = 0. _d 0 duIce (I,J,bi,bj) = 0. _d 0 dvIce (I,J,bi,bj) = 0. _d 0 + ENDDO + ENDDO +C cycle ice velocities + DO J=1-OLy,sNy+OLy + DO I=1-OLx,sNx+OLx + duIcNm1(I,J,bi,bj) = uIce(I,J,bi,bj) * bdfAlpha + & + ( uIce(I,J,bi,bj) - uIceNm1(I,J,bi,bj) ) * bdfFac + dvIcNm1(I,J,bi,bj) = vIce(I,J,bi,bj) * bdfAlpha + & + ( vIce(I,J,bi,bj) - vIceNm1(I,J,bi,bj) ) * bdfFac uIceNm1(I,J,bi,bj) = uIce(I,J,bi,bj) vIceNm1(I,J,bi,bj) = vIce(I,J,bi,bj) ENDDO ENDDO +C As long as IMEX is not properly implemented leave this commented out +CML IF ( .NOT.SEAICEuseIMEX ) THEN C Compute things that do no change during the Newton iteration: C sea-surface tilt and wind stress: -C FORCEX/Y0 - mass*(u/vIceNm1)/deltaT +C FORCEX/Y0 - mass*(1.5*u/vIceNm1+0.5*(u/vIceNm1-u/vIceNm2))/deltaT DO J=1-OLy,sNy+OLy DO I=1-OLx,sNx+OLx FORCEX(I,J,bi,bj) = FORCEX0(I,J,bi,bj) - & + seaiceMassU(I,J,bi,bj)*uIceNm1(I,J,bi,bj)*recip_deltaT + & + seaiceMassU(I,J,bi,bj)*duIcNm1(I,J,bi,bj)*recip_deltaT FORCEY(I,J,bi,bj) = FORCEY0(I,J,bi,bj) - & + seaiceMassV(I,J,bi,bj)*vIceNm1(I,J,bi,bj)*recip_deltaT + & + seaiceMassV(I,J,bi,bj)*dvIcNm1(I,J,bi,bj)*recip_deltaT ENDDO ENDDO +CML ENDIF ENDDO ENDDO C Start nonlinear Newton iteration: outer loop iteration @@ -156,6 +198,7 @@ DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx zetaPre(I,J,bi,bj) = zeta(I,J,bi,bj) + zetaZPre(I,J,bi,bj)= zetaZ(I,J,bi,bj) etaPre(I,J,bi,bj) = eta(I,J,bi,bj) etaZPre(I,J,bi,bj) = etaZ(I,J,bi,bj) dwatPre(I,J,bi,bj) = DWATN(I,J,bi,bj) @@ -194,16 +237,36 @@ krylovConverged = .FALSE. FGMRESeps = JFNKgamma_lin * JFNKresidual +C map first guess sol; it is zero because the solution is a correction + CALL SEAICE_MAP2VEC(nVec,duIce,dvIce,sol,.TRUE.,myThid) +C map rhs and change its sign because we are solving J*u = -F + CALL SEAICE_MAP2VEC(nVec,-uIceRes,-vIceRes,rhs,.TRUE.,myThid) DO WHILE ( .NOT.krylovConverged ) C solution vector sol = du/vIce C residual vector (rhs) Fu = u/vIceRes C output work vectors wk1, -> input work vector wk2 - CALL SEAICE_FGMRES_DRIVER( - I uIceRes, vIceRes, - U duIce, dvIce, iCode, - I FGMRESeps, iOutFGMRES, - I newtonIter, krylovIter, myTime, myIter, myThid ) +C map preconditioner results or Jacobian times vector, +C stored in du/vIce to wk2, for iCode=0, wk2 is set to zero, +C because du/vIce = 0 + CALL SEAICE_MAP2VEC(nVec,duIce,dvIce,wk2,.TRUE.,myThid) +C + CALL SEAICE_FGMRES (nVec,im,rhs,sol,ifgmres,krylovIter, + U vv,w,wk1,wk2, + I FGMRESeps,SEAICEkrylovIterMax,iOutFGMRES, + U iCode, + I myThid) +C + IF ( iCode .EQ. 0 ) THEN +C map sol(ution) vector to du/vIce + CALL SEAICE_MAP2VEC(nVec,duIce,dvIce,sol,.FALSE.,myThid) + ELSE +C map work vector to du/vIce to either compute a preconditioner +C solution (wk1=rhs) or a Jacobian times wk1 + CALL SEAICE_MAP2VEC(nVec,duIce,dvIce,wk1,.FALSE.,myThid) + ENDIF +C Fill overlaps in updated fields + CALL EXCH_UV_XY_RL( duIce, dvIce,.TRUE.,myThid) C FGMRES returns iCode either asking for an new preconditioned vector C or product of matrix (Jacobian) times vector. For iCode = 0, terminate C iteration @@ -212,7 +275,7 @@ IF ( SOLV_MAX_ITERS .GT. 0 ) & CALL SEAICE_PRECONDITIONER( U duIce, dvIce, - I zetaPre, etaPre, etaZpre, dwatPre, + I zetaPre, etaPre, etaZpre, zetaZpre, dwatPre, I newtonIter, krylovIter, myTime, myIter, myThid ) ELSEIF (iCode.GE.2) THEN C Compute Jacobian times vector