--- MITgcm/pkg/seaice/seaice_jfnk.F 2013/03/02 04:35:05 1.20 +++ MITgcm/pkg/seaice/seaice_jfnk.F 2014/10/20 03:20:57 1.27 @@ -1,7 +1,10 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_jfnk.F,v 1.20 2013/03/02 04:35:05 jmc Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_jfnk.F,v 1.27 2014/10/20 03:20:57 gforget 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 @@ -53,9 +56,7 @@ INTEGER myIter INTEGER myThid -#if ( (defined SEAICE_CGRID) && \ - (defined SEAICE_ALLOW_JFNK) && \ - (defined SEAICE_ALLOW_DYNAMICS) ) +#ifdef SEAICE_ALLOW_JFNK C !FUNCTIONS: LOGICAL DIFFERENT_MULTIPLE EXTERNAL DIFFERENT_MULTIPLE @@ -75,10 +76,12 @@ _RL JFNKresidual _RL JFNKresidualKm1 C parameters to compute convergence criterion - _RL phi_e, alp_e, JFNKgamma_lin + _RL JFNKgamma_lin _RL FGMRESeps _RL JFNKtol - +C backward differences extrapolation factors + _RL bdfFac, bdfAlpha +C _RL recip_deltaT LOGICAL JFNKconverged, krylovConverged LOGICAL writeNow @@ -87,6 +90,9 @@ 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) @@ -115,6 +121,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 @@ -123,21 +140,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 @@ -169,10 +198,9 @@ JFNKgamma_lin = JFNKgamma_lin_max IF ( newtonIter.GT.1.AND.newtonIter.LE.SEAICE_JFNK_tolIter & .AND.JFNKresidual.LT.JFNKres_t ) THEN -C Eisenstat, 1996, equ.(2.6) - phi_e = 1. _d 0 - alp_e = 1. _d 0 - JFNKgamma_lin = phi_e*( JFNKresidual/JFNKresidualKm1 )**alp_e +C Eisenstat and Walker (1996), eq.(2.6) + JFNKgamma_lin = SEAICE_JFNKphi + & *( JFNKresidual/JFNKresidualKm1 )**SEAICE_JFNKalpha JFNKgamma_lin = min(JFNKgamma_lin_max, JFNKgamma_lin) JFNKgamma_lin = max(JFNKgamma_lin_min, JFNKgamma_lin) ENDIF @@ -529,7 +557,7 @@ C This is the new residual JFNKresidual = resLoc -#endif /* SEAICE_ALLOW_DYNAMICS and SEAICE_CGRID and SEAICE_ALLOW_JFNK */ +#endif /* SEAICE_ALLOW_JFNK */ RETURN END