--- MITgcm/pkg/seaice/seaice_jfnk.F 2014/02/04 18:30:31 1.24 +++ MITgcm/pkg/seaice/seaice_jfnk.F 2014/12/01 12:31:36 1.28 @@ -1,7 +1,10 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_jfnk.F,v 1.24 2014/02/04 18:30:31 mlosch Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_jfnk.F,v 1.28 2014/12/01 12:31:36 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 @@ -96,6 +99,7 @@ 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) @@ -150,10 +154,11 @@ vIceNm1(I,J,bi,bj) = vIce(I,J,bi,bj) ENDDO ENDDO - IF ( .NOT.SEAICEuseIMEX ) THEN +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*(abA*u/vIceNm1+abB*(u/vIceNm1-u/vIceNm2))/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) @@ -162,7 +167,7 @@ & + seaiceMassV(I,J,bi,bj)*dvIcNm1(I,J,bi,bj)*recip_deltaT ENDDO ENDDO - ENDIF +CML ENDIF ENDDO ENDDO C Start nonlinear Newton iteration: outer loop iteration @@ -183,6 +188,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) @@ -239,7 +245,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