C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_lsr.F,v 1.127 2017/05/23 16:22:45 mlosch Exp $
C $Name:  $

#include "SEAICE_OPTIONS.h"
#ifdef ALLOW_OBCS
# include "OBCS_OPTIONS.h"
#else
# define OBCS_UVICE_OLD
#endif
#ifdef ALLOW_AUTODIFF
# include "AUTODIFF_OPTIONS.h"
#endif

C--  File seaice_lsr.F: seaice LSR dynamical solver S/R:
C--   Contents
C--   o SEAICE_LSR
C--   o SEAICE_RESIDUAL
C--   o SEAICE_LSR_CALC_COEFFS
C--   o SEAICE_LSR_RHSU
C--   o SEAICE_LSR_RHSV
C--   o SEAICE_LSR_TRIDIAGU
C--   o SEAICE_LSR_TRIDIAGV

CBOP
C     !ROUTINE: SEAICE_LSR
C     !INTERFACE:
      SUBROUTINE SEAICE_LSR( myTime, myIter, myThid )

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R SEAICE_LSR
C     | o Solve ice momentum equation with an LSR dynamics solver
C     |   (see Zhang and Hibler,   JGR, 102, 8691-8702, 1997
C     |    and Zhang and Rothrock, MWR, 131,  845- 861, 2003)
C     |   Written by Jinlun Zhang, PSC/UW, Feb-2001
C     |                     zhang@apl.washington.edu
C     *==========================================================*
C     | C-grid version by Martin Losch
C     |   Since 2009/03/18: finite-Volume discretization of
C     |   stress divergence that includes all metric terms
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "GRID.h"
#include "SEAICE_SIZE.h"
#include "SEAICE_PARAMS.h"
#include "SEAICE.h"

#ifdef ALLOW_AUTODIFF_TAMC
# include "AUTODIFF_PARAMS.h"
# include "tamc.h"
#endif

C     !INPUT/OUTPUT PARAMETERS:
C     === Routine arguments ===
C     myTime     :: Simulation time
C     myIter     :: Simulation timestep number
C     myThid     :: my Thread Id. number
      _RL     myTime
      INTEGER myIter
      INTEGER myThid
CEOP

#ifdef SEAICE_CGRID
#ifdef SEAICE_ALLOW_DYNAMICS
C     !FUNCTIONS:
      LOGICAL  DIFFERENT_MULTIPLE
      EXTERNAL 

C     !LOCAL VARIABLES:
C     === Local variables ===
C     i,j,bi,bj  :: Loop counters

      INTEGER ipass, ipass0
      INTEGER i, j, m, bi, bj
      INTEGER iMin, iMax, jMin, jMax
      INTEGER ICOUNT1, ICOUNT2, linearIterLoc, nonLinIterLoc
      INTEGER kSrf
      INTEGER SOLV_MAX_TMP
      LOGICAL doIterate4u, doIterate4v
      CHARACTER*(MAX_LEN_MBUF) msgBuf

      _RL WFAU, WFAV, WFAU2, WFAV2
      _RL S1, S2, S1A, S2A
      _RL eplus, eminus
      _RL recip_deltaT

C     diagonals of coefficient matrices
      _RL AU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL BU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL CU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL AV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL BV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL CV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
C     RHS of TriDiag Linear system (both directions => 5 coef A,B,C & Rt1,2)
      _RL rhsU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL rhsV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
#ifdef SEAICE_GLOBAL_3DIAG_SOLVER
C     RHS of TriDiag Linear system (1 direction only => 3 coef A,B,C)
      _RL rhsUx(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL rhsVy(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      INTEGER iSubIter
      INTEGER errCode
#endif
C     coefficients for lateral points, u(j+/-1)
      _RL uRt1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL uRt2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
C     coefficients for lateral points, v(i+/-1)
      _RL vRt1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL vRt2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
C     abbreviations
      _RL etaPlusZeta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL zetaMinusEta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
C     symmetric drag coefficient
      _RL dragSym     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
C     intermediate time level centered (C) between n and n-1
      _RL uIceC(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL vIceC(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
C     auxillary fields
      _RL uTmp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL vTmp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
C     fractional area at velocity points
      _RL areaW(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL areaS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
#ifdef SEAICE_ALLOW_MOM_ADVECTION
C     tendency due to advection of momentum
      _RL gUmom(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL gVmom(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
#endif /*  SEAICE_ALLOW_MOM_ADVECTION */
#ifdef ALLOW_AUTODIFF_TAMC
      INTEGER itmpkey, itmpkey2, itmpkey3
#endif
      _RL COSWAT
      _RS SINWAT
      _RL UERR
#ifdef SEAICE_ALLOW_CHECK_LSR_CONVERGENCE
      _RL resnorm, EKnorm, counter
#endif /* SEAICE_ALLOW_CHECK_LSR_CONVERGENCE */
      LOGICAL printResidual
      _RL residUini, residVini, residUend, residVend
#ifdef SEAICE_ALLOW_FREEDRIFT
      _RL uRes (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL vRes (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL errIni, errFD, errSum
      _RL residU_fd, residV_fd
      _RL residUmix, residVmix
#endif /* SEAICE_ALLOW_FREEDRIFT */

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

      printResidual = debugLevel.GE.debLevA
     &  .AND. DIFFERENT_MULTIPLE( SEAICE_monFreq, myTime, deltaTClock )

C     extra overlap for (restricted) additive Schwarz method
      jMin = 1   - SEAICE_OLy
      jMax = sNy + SEAICE_OLy
      iMin = 1   - SEAICE_OLx
      iMax = sNx + SEAICE_OLx
C     
      linearIterLoc = SEAICElinearIterMax
      nonLinIterLoc = SEAICEnonLinIterMax
      IF ( SEAICEusePicardAsPrecon ) THEN
       linearIterLoc = SEAICEpreconlinIter
       nonLinIterLoc = SEAICEpreconNL_Iter
      ENDIF
C     
      recip_deltaT = 1. _d 0 / SEAICE_deltaTdyn

#ifdef ALLOW_AUTODIFF_TAMC
cph break artificial dependencies
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          deltaC (i,j,bi,bj) = 0. _d 0
          press  (i,j,bi,bj) = 0. _d 0
          zeta   (i,j,bi,bj) = 0. _d 0
          zetaZ  (i,j,bi,bj) = 0. _d 0
          eta    (i,j,bi,bj) = 0. _d 0
          etaZ   (i,j,bi,bj) = 0. _d 0
          uIceC  (i,j,bi,bj) = 0. _d 0
          vIceC  (i,j,bi,bj) = 0. _d 0
          uIceNm1(i,j,bi,bj) = 0. _d 0
          vIceNm1(i,j,bi,bj) = 0. _d 0
         ENDDO
        ENDDO
       ENDDO
      ENDDO
#endif

C     initialise fractional areas at velocity points
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          areaW(I,J,bi,bj) = 1. _d 0
          areaS(I,J,bi,bj) = 1. _d 0
         ENDDO
        ENDDO
        IF ( SEAICEscaleSurfStress ) THEN
         DO J=jMin,jMax
          DO I=iMin,iMax
           areaW(I,J,bi,bj) =
     &          0.5 _d 0*(AREA(I,J,bi,bj)+AREA(I-1,J,bi,bj))
           areaS(I,J,bi,bj) =
     &          0.5 _d 0*(AREA(I,J,bi,bj)+AREA(I,J-1,bi,bj))
          ENDDO
         ENDDO
        ENDIF
       ENDDO
      ENDDO

      ipass0 = 2
C     modify time stepping a little if more than the default 
C     modified Eulerian time step is used.
      IF ( nonLinIterLoc.GT.2 ) ipass0 = 1
#ifdef ALLOW_AUTODIFF_TAMC
      DO ipass=1,MPSEUDOTIMESTEPS
       itmpkey = (ikey_dynamics-1)*MPSEUDOTIMESTEPS + ipass
#else
      DO ipass=1,nonLinIterLoc
#endif

#ifdef ALLOW_AUTODIFF_TAMC
# ifdef SEAICE_ALLOW_FREEDRIFT
CADJ STORE uice_fd, vice_fd = comlev1_dynsol, kind=isbyte,
CADJ &     key = itmpkey
CADJ STORE ures,vres = comlev1_dynsol, kind=isbyte,
CADJ &     key = itmpkey
# endif
CADJ STORE utmp,vtmp = comlev1_dynsol, kind=isbyte,
CADJ &     key = itmpkey
#endif /* ALLOW_AUTODIFF_TAMC */

      IF ( ipass .LE. nonLinIterLoc ) THEN
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE uice = comlev1_dynsol, kind=isbyte,
CADJ &     key = itmpkey
CADJ STORE vice = comlev1_dynsol, kind=isbyte,
CADJ &     key = itmpkey
CADJ STORE uicenm1 = comlev1_dynsol, kind=isbyte,
CADJ &     key = itmpkey
CADJ STORE vicenm1 = comlev1_dynsol, kind=isbyte,
CADJ &     key = itmpkey
#endif /* ALLOW_AUTODIFF_TAMC */

C     surface level
      kSrf = 1
C--   introduce turning angles
      SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)
      COSWAT=COS(SEAICE_waterTurnAngle*deg2rad)

C SET SOME VALUES
      WFAU=SEAICE_LSRrelaxU
      WFAV=SEAICE_LSRrelaxV
      WFAU2=ZERO
      WFAV2=ZERO
      S1 = 0.
      S2 = 0.
      S1A=0.80 _d 0
      S2A=0.80 _d 0

      IF ( ipass .EQ. 1 ) THEN
C NOW DO PREDICTOR TIME STEP
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           uIceNm1(I,J,bi,bj) = uIce(I,J,bi,bj)
           vIceNm1(I,J,bi,bj) = vIce(I,J,bi,bj)
           uIceC  (I,J,bi,bj) = uIce(I,J,bi,bj)
           vIceC  (I,J,bi,bj) = vIce(I,J,bi,bj)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ELSE
C NOW DO MODIFIED EULER STEP
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           uIce (I,J,bi,bj)=HALF*(uIce(I,J,bi,bj)+uIceNm1(i,j,bi,bj))
           vIce (I,J,bi,bj)=HALF*(vIce(I,J,bi,bj)+vIceNm1(i,j,bi,bj))
           uIceC(I,J,bi,bj)=uIce(I,J,bi,bj)
           vIceC(I,J,bi,bj)=vIce(I,J,bi,bj)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDIF
      IF ( ipass .GT. ipass0 ) THEN
C     for additional (pseudo-time)steps update u/vIceNm1
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           uIceNm1(i,j,bi,bj)=uIce(I,J,bi,bj)
           vIceNm1(i,j,bi,bj)=vIce(I,J,bi,bj)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDIF

#ifdef ALLOW_AUTODIFF_TAMC
C     That is an important one! Note, that
C     * lsr is called multiple times (twice by default), thus the 
C       ipass index is part of itmpkey
C     * this storing is still outside the lsr iteration loop
CADJ STORE uice = comlev1_dynsol, kind=isbyte, key = itmpkey
CADJ STORE vice = comlev1_dynsol, kind=isbyte, key = itmpkey
#endif /* ALLOW_AUTODIFF_TAMC */

      CALL SEAICE_CALC_STRAINRATES(
     I     uIceC, vIceC,
     O     e11, e22, e12,
     I     ipass, myTime, myIter, myThid )

      CALL SEAICE_CALC_VISCOSITIES(
     I     e11, e22, e12, zMin, zMax, hEffM, press0, tensileStrFac,
     O     eta, etaZ, zeta, zetaZ, press, deltaC,
     I     ipass, myTime, myIter, myThid )

      CALL SEAICE_OCEANDRAG_COEFFS(
     I     uIceC, vIceC,
     O     DWATN,
     I     ipass, myTime, myIter, myThid )

#ifdef SEAICE_ALLOW_BOTTOMDRAG
      IF (SEAICEbasalDragK2.GT.0. _d 0) CALL SEAICE_BOTTOMDRAG_COEFFS(
     I     uIceC, vIceC,
#ifdef SEAICE_ITD
     I     HEFFITD, AREAITD, AREA,
#else
     I     HEFF, AREA,
#endif      
     O     CbotC,
     I     ipass, myTime, myIter, myThid )
#endif /* SEAICE_ALLOW_BOTTOMDRAG */
C
C     some abbreviations
C
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO J=jMin-1,jMax
         DO I=iMin-1,iMax
          etaPlusZeta (I,J,bi,bj) = ETA (I,J,bi,bj)+ZETA(I,J,bi,bj)
          zetaMinusEta(I,J,bi,bj) = ZETA(I,J,bi,bj)-ETA (I,J,bi,bj)
         ENDDO
        ENDDO
        DO J=1-OLy,sNy+OLy
         DO I=1-OLx,sNx+OLx
          dragSym(I,J,bi,bj) = DWATN(I,J,bi,bj)*COSWAT
#ifdef SEAICE_ALLOW_BOTTOMDRAG
     &         +CbotC(I,J,bi,bj)
#endif /* SEAICE_ALLOW_BOTTOMDRAG */
         ENDDO
        ENDDO
       ENDDO
      ENDDO
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C
C     Set up of right hand sides of momentum equations starts here
C
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO J=jMin,jMax
         DO I=iMin,iMax
C     set up anti symmetric drag force and add in current force
C     ( remember to average to correct velocity points )
          FORCEX(I,J,bi,bj)=FORCEX0(I,J,bi,bj)+
     &         ( 0.5 _d 0 * ( DWATN(I,J,bi,bj)+DWATN(I-1,J,bi,bj) ) *
     &         COSWAT * uVel(I,J,kSrf,bi,bj)
     &         - SIGN(SINWAT, _fCori(I,J,bi,bj))* 0.5 _d 0 *
     &         ( DWATN(I  ,J,bi,bj) * 0.5 _d 0 *
     &          (vVel(I  ,J  ,kSrf,bi,bj)-vIceC(I  ,J  ,bi,bj)
     &          +vVel(I  ,J+1,kSrf,bi,bj)-vIceC(I  ,J+1,bi,bj))
     &         + DWATN(I-1,J,bi,bj) * 0.5 _d 0 *
     &          (vVel(I-1,J  ,kSrf,bi,bj)-vIceC(I-1,J  ,bi,bj)
     &          +vVel(I-1,J+1,kSrf,bi,bj)-vIceC(I-1,J+1,bi,bj))
     &         ) ) * areaW(I,J,bi,bj)
          FORCEY(I,J,bi,bj)=FORCEY0(I,J,bi,bj)+
     &         ( 0.5 _d 0 * ( DWATN(I,J,bi,bj)+DWATN(I,J-1,bi,bj) ) *
     &         COSWAT * vVel(I,J,kSrf,bi,bj)
     &         + SIGN(SINWAT, _fCori(I,J,bi,bj)) * 0.5 _d 0 *
     &         ( DWATN(I,J  ,bi,bj) * 0.5 _d 0 *
     &          (uVel(I  ,J  ,kSrf,bi,bj)-uIceC(I  ,J  ,bi,bj)
     &          +uVel(I+1,J  ,kSrf,bi,bj)-uIceC(I+1,J  ,bi,bj))
     &         + DWATN(I,J-1,bi,bj) * 0.5 _d 0 *
     &          (uVel(I  ,J-1,kSrf,bi,bj)-uIceC(I  ,J-1,bi,bj)
     &          +uVel(I+1,J-1,kSrf,bi,bj)-uIceC(I+1,J-1,bi,bj))
     &         ) ) * areaS(I,J,bi,bj)
         ENDDO
        ENDDO
        DO j=jMin,jMax
         DO i=iMin,iMax
C-    add Coriolis term
          FORCEX(I,J,bi,bj) = FORCEX(I,J,bi,bj) + HALF*
     &         ( seaiceMassC(I  ,J,bi,bj) * _fCori(I  ,J,bi,bj)
     &          *0.5 _d 0*( vIceC( i ,j,bi,bj)+vIceC( i ,j+1,bi,bj) )
     &         + seaiceMassC(I-1,J,bi,bj) * _fCori(I-1,J,bi,bj)
     &          *0.5 _d 0*( vIceC(i-1,j,bi,bj)+vIceC(i-1,j+1,bi,bj) ) )
          FORCEY(I,J,bi,bj) = FORCEY(I,J,bi,bj) - HALF*
     &         ( seaiceMassC(I,J  ,bi,bj) * _fCori(I,J  ,bi,bj)
     &         *0.5 _d 0*( uIceC(i  ,j  ,bi,bj)+uIceC(i+1,  j,bi,bj) )
     &         + seaiceMassC(I,J-1,bi,bj) * _fCori(I,J-1,bi,bj)
     &         *0.5 _d 0*( uIceC(i  ,j-1,bi,bj)+uIceC(i+1,j-1,bi,bj) ) )
         ENDDO
        ENDDO
C      this is the rhs contribution of the time derivative
        DO j=jMin,jMax
         DO i=iMin,iMax
          FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)
     &         +seaiceMassU(I,J,bi,bj)*recip_deltaT
     &         *uIceNm1(i,j,bi,bj)
          FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)
     &         +seaiceMassV(I,J,bi,bj)*recip_deltaT
     &         *vIceNm1(i,j,bi,bj)
          FORCEX(I,J,bi,bj)=FORCEX(I,J,bi,bj)* seaiceMaskU(I,J,bi,bj)
          FORCEY(I,J,bi,bj)=FORCEY(I,J,bi,bj)* seaiceMaskV(I,J,bi,bj)
         ENDDO
        ENDDO
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C
C     u-equation
C
        DO j=jMin,jMax
         DO i=iMin,iMax
          rhsU(I,J,bi,bj) = FORCEX(I,J,bi,bj) 
         ENDDO
        ENDDO
        CALL SEAICE_LSR_RHSU(
     I       zetaMinusEta, etaPlusZeta, etaZ, zetaZ, press,
     I       uIceC, vIceC, 
     U       rhsU,
     I       iMin, iMax, jMin, jMax, bi, bj, myThid )
C
C     v-equation
C
        DO j=jMin,jMax
         DO i=iMin,iMax
          rhsV(I,J,bi,bj) = FORCEY(I,J,bi,bj) 
         ENDDO
        ENDDO
        CALL SEAICE_LSR_RHSV(
     I       zetaMinusEta, etaPlusZeta, etaZ, zetaZ, press,
     I       uIceC, vIceC, 
     U       rhsV,
     I       iMin, iMax, jMin, jMax, bi, bj, myThid )
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef SEAICE_ALLOW_MOM_ADVECTION
        IF ( SEAICEmomAdvection ) THEN 
         DO J=1-Oly,sNy+Oly
          DO I=1-Olx,sNx+Olx
           gUmom(I,J) = 0. _d 0
           gVmom(I,J) = 0. _d 0
          ENDDO
         ENDDO
         CALL SEAICE_MOM_ADVECTION(
     I        bi,bj,iMin,iMax,jMin,jMax,
     I        uIceC, vIceC,
     O        gUmom, gVmom,
     I        myTime, myIter, myThid )
         DO J=jMin,jMax
          DO I=jMin,jMax
           rhsU(I,J,bi,bj) = rhsU(I,J,bi,bj) + gUmom(I,J)
           rhsV(I,J,bi,bj) = rhsV(I,J,bi,bj) + gVmom(I,J)
          ENDDO
         ENDDO
        ENDIF
#endif /* SEAICE_ALLOW_MOM_ADVECTION */
       ENDDO
      ENDDO
C
C     Set up of right hand sides of momentum equations ends here
C
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C
C     calculate coefficients of tridiagonal matrices for both u- and
C     v-equations
C
      CALL SEAICE_LSR_CALC_COEFFS( 
     I     etaPlusZeta, zetaMinusEta, etaZ, zetaZ, dragSym,
     O     AU, BU, CU, AV, BV, CV, uRt1, uRt2, vRt1, vRt2,
     I     iMin, iMax, jMin, jMax, myTime, myIter, myThid )

#ifndef OBCS_UVICE_OLD
C--     prevent tri-diagonal solver from modifying OB values:
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO J=jMin,jMax
         DO I=iMin,iMax
          IF ( maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj) .EQ. 0. ) THEN
           AU(I,J,bi,bj)   = ZERO
           BU(I,J,bi,bj)   = ONE
           CU(I,J,bi,bj)   = ZERO
           uRt1(I,J,bi,bj) = ZERO
           uRt2(I,J,bi,bj) = ZERO
           rhsU(I,J,bi,bj) = uIce(I,J,bi,bj)
          ENDIF
          IF ( maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj) .EQ. 0. ) THEN
           AV(I,J,bi,bj)   = ZERO
           BV(I,J,bi,bj)   = ONE
           CV(I,J,bi,bj)   = ZERO
           vRt1(I,J,bi,bj) = ZERO
           vRt2(I,J,bi,bj) = ZERO
           rhsV(I,J,bi,bj) = vIce(I,J,bi,bj)
          ENDIF
         ENDDO
        ENDDO
       ENDDO
      ENDDO
#endif /* OBCS_UVICE_OLD */

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

#ifdef ALLOW_DEBUG
      IF ( debugLevel .GE. debLevD ) THEN
        WRITE(msgBuf,'(A,I3,A)')
     &        'Uice pre  iter (SEAICE_LSR', MOD(ipass,1000), ')'
        CALL DEBUG_STATS_RL( 1, UICE, msgBuf, myThid )
        WRITE(msgBuf,'(A,I3,A)')
     &        'Vice pre  iter (SEAICE_LSR', MOD(ipass,1000), ')'
        CALL DEBUG_STATS_RL( 1, VICE, msgBuf, myThid )
      ENDIF
#endif /* ALLOW_DEBUG */

C--   Calculate initial residual of the linearised system
      IF ( printResidual .OR. LSR_mixIniGuess.GE.1 )
     & CALL SEAICE_RESIDUAL(
     I                  rhsU, rhsV, uRt1, uRt2, vRt1, vRt2,
     I                  AU, BU, CU, AV, BV, CV, uIce, vIce,
     O                  residUini, residVini, uTmp, vTmp,
     I                  printResidual, myIter, myThid )

#ifdef SEAICE_ALLOW_FREEDRIFT
      IF ( printResidual .OR. LSR_mixIniGuess.GE.1 ) THEN
       IF ( LSR_mixIniGuess.GE.1 ) THEN
        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          DO j=jMin,jMax
           DO i=iMin,iMax
            uIce_fd(i,j,bi,bj) = FORCEX(i,j,bi,bj)
     &         / ( 1. _d 0 - seaiceMaskU(i,j,bi,bj)
     &           + seaiceMassU(i,j,bi,bj)*recip_deltaT
     &           + HALF*( dragSym(i,j,bi,bj) + dragSym(i-1,j,bi,bj) )
     &                 * areaW(i,j,bi,bj)
     &           )
            vIce_fd(i,j,bi,bj) = FORCEY(i,j,bi,bj)
     &         / ( 1. _d 0 - seaiceMaskV(i,j,bi,bj)
     &           + seaiceMassV(i,j,bi,bj)*recip_deltaT
     &           + HALF*( dragSym(i,j,bi,bj) + dragSym(i,j-1,bi,bj) )
     &                 * areaS(i,j,bi,bj)
     &           )
           ENDDO
          ENDDO
         ENDDO
        ENDDO
        CALL EXCH_UV_XY_RL( uIce_fd, vIce_fd, .TRUE., myThid )
       ENDIF
       IF ( LSR_mixIniGuess.GE.0 )
     &  CALL SEAICE_RESIDUAL(
     I                  rhsU, rhsV, uRt1, uRt2, vRt1, vRt2,
     I                  AU, BU, CU, AV, BV, CV, uIce_fd, vIce_fd,
     O                  residU_fd, residV_fd, uRes, vRes,
     I                  printResidual, myIter, myThid )
      ENDIF
      IF ( LSR_mixIniGuess.GE.2 ) THEN
# ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE uice, vice = comlev1_dynsol, kind=isbyte,
CADJ &     key = itmpkey
CADJ STORE utmp, vtmp = comlev1_dynsol, kind=isbyte,
CADJ &     key = itmpkey
# endif /* ALLOW_AUTODIFF_TAMC */
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         DO j=jMin,jMax
          DO i=iMin,iMax
           errIni = uTmp(i,j,bi,bj)*uTmp(i,j,bi,bj)
           errFD  = uRes(i,j,bi,bj)*uRes(i,j,bi,bj)
           IF ( LSR_mixIniGuess.GE.4 ) THEN
            errIni = errIni*errIni
            errFD  = errFD *errFD
           ENDIF
           errSum = ( errIni + errFD )
     &          *maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj)
           IF ( errSum.GT.0. ) THEN
            uIce(i,j,bi,bj) = ( errFD *uIce(i,j,bi,bj)
     &                        + errIni*uIce_fd(i,j,bi,bj)
     &                        )/errSum
           ENDIF
           errIni = vTmp(i,j,bi,bj)*vTmp(i,j,bi,bj)
           errFD  = vRes(i,j,bi,bj)*vRes(i,j,bi,bj)
           IF ( LSR_mixIniGuess.GE.4 ) THEN
            errIni = errIni*errIni
            errFD  = errFD *errFD
           ENDIF
           errSum = ( errIni + errFD )
     &          *maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj)
           IF ( errSum.GT.0. ) THEN
            vIce(i,j,bi,bj) = ( errFD *vIce(i,j,bi,bj)
     &                        + errIni*vIce_fd(i,j,bi,bj)
     &                        )/errSum
           ENDIF
          ENDDO
         ENDDO
        ENDDO
       ENDDO
       CALL EXCH_UV_XY_RL( uIce, vIce, .TRUE., myThid )
#ifndef ALLOW_AUTODIFF_TAMC
C     Here residU/Vmix and u/vTmp are computed but never used for anything
C     but reporting to STDOUT. Still TAF thinks this requires recomputations
C     so we hide this part of the code from TAF.
       IF ( printResidual ) THEN
        CALL SEAICE_RESIDUAL(
     I                  rhsU, rhsV, uRt1, uRt2, vRt1, vRt2,
     I                  AU, BU, CU, AV, BV, CV, uIce, vIce,
     O                  residUmix, residVmix, uTmp, vTmp,
     I                  printResidual, myIter, myThid )
       ENDIF
#endif /* ALLOW_AUTODIFF_TAMC */
      ENDIF
#endif /* SEAICE_ALLOW_FREEDRIFT */

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C NOW DO ITERATION

cph--- iteration starts here
cph--- need to kick out goto
      doIterate4u = .TRUE.
      doIterate4v = .TRUE.

C ITERATION START -----------------------------------------------------
C     Ideally, we would like to use the loop-directive for this
C     iterative loop, but (unsurprisingly) it does not seem to work for 
C     this convoluted case of mixed iterations, so it is commented out
CML#ifdef ALLOW_AUTODIFF_TAMC
CMLCADJ LOOP = iteration uice, vice = comlev1_lsr
CML#endif /* ALLOW_AUTODIFF_TAMC */

#if (defined (ALLOW_AUTODIFF_TAMC)  defined (SEAICE_LSR_ADJOINT_ITER))
      IF ( inAdMode ) THEN
        SOLV_MAX_TMP = SOLV_MAX_FIXED
      ELSE
        SOLV_MAX_TMP = linearIterLoc
      ENDIF
#else
      SOLV_MAX_TMP = linearIterLoc
#endif

      ICOUNT1 = SOLV_MAX_TMP
      ICOUNT2 = SOLV_MAX_TMP

      DO m = 1, SOLV_MAX_TMP

#ifdef ALLOW_AUTODIFF_TAMC
# ifdef SEAICE_LSR_ADJOINT_ITER
       itmpkey2 = (itmpkey-1)*SOLV_MAX_FIXED + MIN(m,SOLV_MAX_FIXED)
 else
       itmpkey2 = itmpkey
# endif /* SEAICE_LSR_ADJOINT_ITER */
C     Storing these scalars and logicals is necessary to do an exact
C     backward iteration
CADJ STORE wfau, wfav = comlev1_lsr,  kind=isbyte, key = itmpkey2
CADJ STORE doIterate4u, doIterate4v = comlev1_lsr, key = itmpkey2
#endif /* ALLOW_AUTODIFF_TAMC */

       IF ( doIterate4u .OR. doIterate4v ) THEN

        IF ( useCubedSphereExchange ) THEN
          doIterate4u = .TRUE.
          doIterate4v = .TRUE.
        ENDIF

# ifdef ALLOW_AUTODIFF_TAMC
C     Note that itmpkey2 can get very large when SEAICE_LSR_ADJOINT_ITER
C     is defined for an accurate adjoint, otherwise it is the same for
C     each m, so that the tapes get overwritten for each iterate and
C     the adjoint is necessarily wrong.
CADJ STORE uice, vice = comlev1_lsr, kind=isbyte, key = itmpkey2
# endif /* ALLOW_AUTODIFF_TAMC */

#ifdef SEAICE_GLOBAL_3DIAG_SOLVER
        IF ( SEAICEuseMultiTileSolver ) THEN

         DO bj=myByLo(myThid),myByHi(myThid)
          DO bi=myBxLo(myThid),myBxHi(myThid)
           IF ( doIterate4u ) THEN
            DO j=jMin,jMax
             DO i=iMin,iMax
              uTmp(i,j,bi,bj)  = uIce(i,j,bi,bj)
              rhsUx(i,j,bi,bj) = ( rhsU(i,j,bi,bj)
     &             + uRt1(i,j,bi,bj)*uIce(i,j-1,bi,bj)
     &             + uRt2(i,j,bi,bj)*uIce(i,j+1,bi,bj)
     &                           )*seaiceMaskU(i,j,bi,bj)
             ENDDO
            ENDDO
           ENDIF
           IF ( doIterate4v ) THEN
            DO j=jMin,jMax
             DO i=iMin,iMax
              vTmp(i,j,bi,bj)  = vIce(i,j,bi,bj)
              rhsVy(i,j,bi,bj) = ( rhsV(i,j,bi,bj)
     &             + vRt1(i,j,bi,bj)*vIce(i-1,J,bi,bj)
     &             + vRt2(i,j,bi,bj)*vIce(i+1,J,bi,bj)
     &                           )*seaiceMaskU(i,j,bi,bj)
             ENDDO
            ENDDO
           ENDIF
C     end bi,bj-loops
          ENDDO
         ENDDO

         iSubIter = SOLV_MAX_TMP*(ipass-1) + m
         CALL SOLVE_UV_TRIDIAGO(
     I                 1, OLx, doIterate4u, doIterate4v,
     I                 AU, BU, CU, rhsUx,
     I                 AV, BV, CV, rhsVy,
     O                 uIce, vIce, errCode,
     I                 iSubIter, myIter, myThid )

         DO bj=myByLo(myThid),myByHi(myThid)
          DO bi=myBxLo(myThid),myBxHi(myThid)
           IF ( doIterate4u ) THEN
            DO j=jMin,jMax
             DO i=iMin,iMax
              uIce(i,j,bi,bj) = uTmp(i,j,bi,bj)
     &            + WFAU*( uIce(i,j,bi,bj)-uTmp(i,j,bi,bj) )
             ENDDO
            ENDDO
           ENDIF
           IF ( doIterate4v ) THEN
            DO j=jMin,jMax
             DO i=iMin,iMax
              vIce(i,j,bi,bj) = vTmp(i,j,bi,bj)
     &            + WFAV*( vIce(i,j,bi,bj)-vTmp(i,j,bi,bj) )
             ENDDO
            ENDDO
           ENDIF
C     end bi,bj-loops
          ENDDO
         ENDDO

        ELSE
#endif /* SEAICE_GLOBAL_3DIAG_SOLVER */

        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)

#ifdef ALLOW_AUTODIFF_TAMC
          act1 = bi - myBxLo(myThid)
          max1 = myBxHi(myThid) - myBxLo(myThid) + 1
          act2 = bj - myByLo(myThid)
          max2 = myByHi(myThid) - myByLo(myThid) + 1
          act3 = myThid - 1
          max3 = nTx*nTy
          act4 = itmpkey2 - 1
          itmpkey3 = (act1 + 1) + act2*max1
     &                          + act3*max1*max2
     &                          + act4*max1*max2*max3
CADJ STORE uice(:,:,bi,bj) = comlev1_bibj_lsr,kind=isbyte,key=itmpkey3
CADJ STORE vice(:,:,bi,bj) = comlev1_bibj_lsr,kind=isbyte,key=itmpkey3
#endif /* ALLOW_AUTODIFF_TAMC */
C-jmc: get less TAF warnings when always (no if doIterate) saving uIce,vIce:
C     save u/vIce prior to iteration
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
            uTmp(I,J,bi,bj)=uIce(I,J,bi,bj)
            vTmp(I,J,bi,bj)=vIce(I,J,bi,bj)
           ENDDO
          ENDDO

          IF ( doIterate4u ) THEN
C Solve for uIce :
           CALL SEAICE_LSR_TRIDIAGU(
     I          AU, BU, CU, uRt1, uRt2, rhsU, uTmp, seaiceMaskU, WFAU,
     U          uIce,
     I          iMin, iMax, jMin, jMax, bi, bj, myTime, myIter, myThid )
          ENDIF

          IF ( doIterate4v ) THEN
C Solve for vIce
           CALL SEAICE_LSR_TRIDIAGV(
     I          AV, BV, CV, vRt1, vRt2, rhsV, vTmp, seaiceMaskV, WFAV, 
     U          vIce,
     I          iMin, iMax, jMin, jMax, bi, bj, myTime, myIter, myThid )
          ENDIF

C     end bi,bj-loops
         ENDDO
        ENDDO

#ifdef SEAICE_GLOBAL_3DIAG_SOLVER
        ENDIF
#endif /* SEAICE_GLOBAL_3DIAG_SOLVER */

        IF ( doIterate4u.AND.MOD(m,SOLV_NCHECK).EQ.0) THEN
         S1=ZERO
         DO bj=myByLo(myThid),myByHi(myThid)
          DO bi=myBxLo(myThid),myBxHi(myThid)
           DO J=1,sNy
            DO I=1,sNx
             UERR=(uIce(I,J,bi,bj)-uTmp(I,J,bi,bj))
     &               * seaiceMaskU(I,J,bi,bj)
             S1=MAX(ABS(UERR),S1)
            ENDDO
           ENDDO
          ENDDO
         ENDDO
         _GLOBAL_MAX_RL( S1, myThid )
c        WRITE(standardMessageUnit,'(A,2I6,1P4E16.9)')
c    &   ' U iters,error,WF = ',ipass,M,S1,S1A,WFAU
C SAFEGUARD AGAINST BAD FORCING ETC
         IF(m.GT.1.AND.S1.GT.S1A) WFAU=WFAU2
         S1A=S1
         IF(S1.LT.LSR_ERROR) THEN
          ICOUNT1=m
          doIterate4u = .FALSE.
         ENDIF
        ENDIF

        IF ( doIterate4v.AND.MOD(m,SOLV_NCHECK).EQ.0) THEN
         S2=ZERO
         DO bj=myByLo(myThid),myByHi(myThid)
          DO bi=myBxLo(myThid),myBxHi(myThid)
           DO J=1,sNy
            DO I=1,sNx
             UERR=(vIce(I,J,bi,bj)-vTmp(I,J,bi,bj))
     &               * seaiceMaskV(I,J,bi,bj)
             S2=MAX(ABS(UERR),S2)
            ENDDO
           ENDDO
          ENDDO
         ENDDO
         _GLOBAL_MAX_RL( S2, myThid )
C SAFEGUARD AGAINST BAD FORCING ETC
         IF(m.GT.1.AND.S2.GT.S2A) WFAV=WFAV2
         S2A=S2
         IF(S2.LT.LSR_ERROR) THEN
          ICOUNT2=m
          doIterate4v = .FALSE.
         ENDIF
        ENDIF

        CALL EXCH_UV_XY_RL( uIce, vIce,.TRUE.,myThid)

C--    end doIterate4u or doIterate4v
       ENDIF
      ENDDO
C ITERATION END -----------------------------------------------------

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

#ifndef ALLOW_AUTODIFF_TAMC
C     Here residU/Vend and u/vTmp are computed but never used for anything
C     but reporting to STDOUT. Still TAF thinks this requires recomputations
C     so we hide this part of the code from TAF.
      IF ( printResidual ) THEN
C--   Calculate final residual of the linearised system
        CALL SEAICE_RESIDUAL(
     I                  rhsU, rhsV, uRt1, uRt2, vRt1, vRt2,
     I                  AU, BU, CU, AV, BV, CV, uIce, vIce,
     O                  residUend, residVend, uTmp, vTmp,
     I                  printResidual, myIter, myThid )
        _BEGIN_MASTER( myThid )
        WRITE(standardMessageUnit,'(A,1X,I4,1P3E16.8)')
     &    ' SEAICE_LSR: Residual Initial ipass,Uice,Vice=',
     &        ipass, residUini, residVini
#ifdef SEAICE_ALLOW_FREEDRIFT
        IF ( LSR_mixIniGuess.GE.0 )
     &  WRITE(standardMessageUnit,'(A,1P3E16.8)')
     &    ' SEAICE_LSR: Residual FrDrift U_fd,V_fd=',
     &        residU_fd, residV_fd
        IF ( LSR_mixIniGuess.GE.2 )
     &  WRITE(standardMessageUnit,'(A,1P3E16.8)')
     &    ' SEAICE_LSR: Residual Mix-ini Uice,Vice=',
     &        residUmix, residVmix
#endif /* SEAICE_ALLOW_FREEDRIFT */
        WRITE(standardMessageUnit,'(A,I4,A,I6,1P2E16.8)')
     &    ' SEAICE_LSR (ipass=',ipass,') iters,dU,Resid=',
     &      ICOUNT1, S1, residUend
        WRITE(standardMessageUnit,'(A,I4,A,I6,1P2E16.8)')
     &    ' SEAICE_LSR (ipass=',ipass,') iters,dV,Resid=',
     &      ICOUNT2, S2, residVend
        _END_MASTER( myThid )
      ENDIF
#ifdef ALLOW_DEBUG
      IF ( debugLevel .GE. debLevD ) THEN
        WRITE(msgBuf,'(A,I3,A)')
     &        'Uice post iter (SEAICE_LSR', MOD(ipass,1000), ')'
        CALL DEBUG_STATS_RL( 1, UICE, msgBuf, myThid )
        WRITE(msgBuf,'(A,I3,A)')
     &        'Vice post iter (SEAICE_LSR', MOD(ipass,1000), ')'
        CALL DEBUG_STATS_RL( 1, VICE, msgBuf, myThid )
      ENDIF
#endif /* ALLOW_DEBUG */
      IF ( doIterate4u .OR. doIterate4v ) THEN
        WRITE(msgBuf,'(2A,I10,A,I4,A)') '** WARNING ** SEAICE_LSR ',
     &    '(it=', myIter, ',', ipass,') did not converge :'
        CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &                      SQUEEZE_RIGHT, myThid )
        WRITE(msgBuf,'(2(A,I6,0PF6.3,1PE14.6))')
     &      ' nIt,wFU,dU=', ICOUNT1, WFAU, S1,
     &    ' ; nIt,wFV,dV=', ICOUNT2, WFAV, S2
        CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &                      SQUEEZE_RIGHT, myThid )
      ENDIF
#endif /* ALLOW_AUTODIFF_TAMC */

C     APPLY MASKS
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO J=1-OLy,sNy+OLy
         DO I=1-OLx,sNx+OLx
          uIce(I,J,bi,bj)=uIce(I,J,bi,bj)* seaiceMaskU(I,J,bi,bj)
          vIce(I,J,bi,bj)=vIce(I,J,bi,bj)* seaiceMaskV(I,J,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
CML      CALL EXCH_UV_XY_RL( uIce, vIce,.TRUE.,myThid)

#ifdef SEAICE_ALLOW_CHECK_LSR_CONVERGENCE
      IF ( debugLevel .GE. debLevC ) THEN
       resnorm = 0. _d 0
       EKnorm  = 0. _d 0
       counter = 0. _d 0
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
C--   Compute residual norm and print
         DO j=1,sNy
          DO i=1,sNx
           resnorm = resnorm + 0.5 _d 0 *
     &          ( ( (uIceNm1(i,j,bi,bj)+uIceNm1(i+1,j,bi,bj))
     &            - (uIce(i,j,bi,bj)+uIce(i+1,j,bi,bj)) )**2
     &          + ( (vIceNm1(i,j,bi,bj)+vIceNm1(i,j+1,bi,bj))
     &            - (vIce(i,j,bi,bj)+vIce(i,j+1,bi,bj)) )**2 )
           IF ( area(i,j,bi,bj) .gt. 0.5 _d 0 ) THEN
            EKnorm = EKnorm + 0.5 _d 0 * heff(i,j,bi,bj) *
     &           ( ( (uIce(i,j,bi,bj)+uIce(i+1,j,bi,bj)) )**2
     &           + ( (vIce(i,j,bi,bj)+vIce(i,j+1,bi,bj)) )**2 )
            counter = counter + 1. _d 0
           ENDIF
          ENDDO
         ENDDO
        ENDDO
       ENDDO
       _GLOBAL_SUM_RL( resnorm, myThid )
       _GLOBAL_SUM_RL( EKnorm, myThid )
       _GLOBAL_SUM_RL( counter, myThid )
       IF ( counter .gt. 0. _d 0 ) EKnorm = EKnorm/counter
       _BEGIN_MASTER( myThid )
       WRITE(standardMessageUnit,'(A,I7,1X,2E22.14)')
     &      'S/R SEAICE_LSR: IPSEUDO, RESNORM, EKNORM = ',
     &      ipass, sqrt(resnorm), EKnorm
       _END_MASTER( myThid )
      ENDIF
#endif /* SEAICE_ALLOW_CHECK_LSR_CONVERGENCE */

      ENDIF
C     end outer ipass pseudo-time stepping loop
      ENDDO

      IF ( useHB87StressCoupling ) THEN
# ifdef ALLOW_AUTODIFF_TAMC
C     These directives do not remove any recomputation warnings but avoid
C     recomputing the entire LSR loop before evaluating this code block
CADJ STORE uice, vice, eta, etaz, zeta = comlev1_dynsol, 
CADJ &     kind=isbyte, key = ikey_dynamics
# endif /* ALLOW_AUTODIFF_TAMC */
C     compute the divergence of stress here to be used later
C
C     compute strain rate from latest velocities
       CALL SEAICE_CALC_STRAINRATES(
     I       uIce, vIce,
     O       e11, e22, e12,
     I       3, myTime, myIter, myThid )
C     compute internal stresses with updated ice velocities
C     and evaluate divergence of stress and apply to forcing
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         CALL SEAICE_CALC_STRESSDIV( 
     I        e11, e22, e12, press, zeta, eta, etaZ,
     O        stressDivergenceX, stressDivergenceY,
     I        bi, bj, myTime, myIter, myThid )
        ENDDO
       ENDDO
C     endif  useHB87StressCoupling
      ENDIF

#endif /* SEAICE_ALLOW_DYNAMICS */
#endif /* SEAICE_CGRID */

      RETURN
      END


C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP C !ROUTINE: SEAICE_RESIDUAL C !INTERFACE: SUBROUTINE SEAICE_RESIDUAL( I rhsU, rhsV, uRt1, uRt2, vRt1, vRt2, I AU, BU, CU, AV, BV, CV, uFld, vFld, O residU, residV, uRes, vRes, I calcMeanResid, myIter, myThid ) C !DESCRIPTION: \bv C *==========================================================* C | S/R SEAICE_RESIDUAL C | o Compute Linear solver residual C *==========================================================* C \ev C !USES: IMPLICIT NONE C === Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "GRID.h" #include "SEAICE_SIZE.h" #include "SEAICE.h" C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === C calcMeanResid :: also compute mean residual (=RMS) C myIter :: Simulation timestep number C myThid :: my Thread Id. number C-- RHS _RL rhsU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL rhsV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C- coefficients for lateral points, u(j+/-1) _RL uRt1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL uRt2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C- coefficients for lateral points, v(i+/-1) _RL vRt1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL vRt2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C- diagonals of coefficient matrices _RL AU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL BU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL CU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL AV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL BV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL CV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C-- seaice velocity _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C- residual (output) _RL residU, residV _RL uRes (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL vRes (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) LOGICAL calcMeanResid INTEGER myIter INTEGER myThid CEOP #ifdef SEAICE_CGRID #ifdef SEAICE_ALLOW_DYNAMICS C !LOCAL VARIABLES: C === Local variables === C i,j,bi,bj :: Loop counters INTEGER i, j, bi, bj C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C-- Calculate residual of the linearised system residU = 0. residV = 0. DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1,sNy DO i=1,sNx uRes(i,j,bi,bj) = rhsU(i,j,bi,bj) & + uRt1(i,j,bi,bj)*uFld(i,j-1,bi,bj) & + uRt2(i,j,bi,bj)*uFld(i,j+1,bi,bj) & - ( AU(i,j,bi,bj)*uFld(i-1,j,bi,bj) & + BU(i,j,bi,bj)*uFld( i ,j,bi,bj) & + CU(i,j,bi,bj)*uFld(i+1,j,bi,bj) & ) vRes(i,j,bi,bj) = rhsV(i,j,bi,bj) & + vRt1(i,j,bi,bj)*vFld(i-1,j,bi,bj) & + vRt2(i,j,bi,bj)*vFld(i+1,j,bi,bj) & - ( AV(i,j,bi,bj)*vFld(i,j-1,bi,bj) & + BV(i,j,bi,bj)*vFld(i, j ,bi,bj) & + CV(i,j,bi,bj)*vFld(i,j+1,bi,bj) & ) ENDDO ENDDO IF ( calcMeanResid ) THEN DO j=1,sNy DO i=1,sNx residU = residU + uRes(i,j,bi,bj)*uRes(i,j,bi,bj) & *rAw(i,j,bi,bj)*maskInW(i,j,bi,bj) #ifndef OBCS_UVICE_OLD & *maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj) #endif residV = residV + vRes(i,j,bi,bj)*vRes(i,j,bi,bj) & *rAs(i,j,bi,bj)*maskInS(i,j,bi,bj) #ifndef OBCS_UVICE_OLD & *maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj) #endif ENDDO ENDDO ENDIF ENDDO ENDDO IF ( calcMeanResid ) THEN _GLOBAL_SUM_RL( residU, myThid ) _GLOBAL_SUM_RL( residV, myThid ) C scale residuals by globalArea so that they do not get ridiculously large IF ( residU.GT.0. ) residU = SQRT(residU/globalArea) IF ( residV.GT.0. ) residV = SQRT(residV/globalArea) ENDIF C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| #endif /* SEAICE_ALLOW_DYNAMICS */ #endif /* SEAICE_CGRID */ RETURN END


CBOP C !ROUTINE: SEAICE_LSR_CALC_COEFFS C !INTERFACE: SUBROUTINE SEAICE_LSR_CALC_COEFFS( I etaPlusZeta, zetaMinusEta, etaZloc, zetaZloc, dragSym, O AU, BU, CU, AV, BV, CV, uRt1, uRt2, vRt1, vRt2, I iMin, iMax, jMin, jMax, myTime, myIter, myThid ) C !DESCRIPTION: \bv C *==========================================================* C | S/R SEAICE_LSR_CALC_COEFFS C | o Calculate coefficient matrix for LSR solver C *==========================================================* C \ev C !USES: IMPLICIT NONE C === Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "SEAICE_SIZE.h" #include "SEAICE_PARAMS.h" #include "SEAICE.h" C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === C myTime :: Simulation time C myIter :: Simulation timestep number C myThid :: my Thread Id. number _RL myTime INTEGER myIter INTEGER myThid INTEGER iMin, iMax, jMin, jMax C combinations of bulk and shear viscosity at C points C abbreviations _RL etaPlusZeta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL zetaMinusEta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C shear and bulk viscosity at Z points _RL etaZloc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL zetaZloc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C symmetric ice-ocean stress coefficient _RL dragSym (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C coefficients of ice velocities in coefficient matrix C for both U and V-equation C diagonals of coefficient matrices _RL AU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL BU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL CU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL AV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL BV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL CV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C coefficients for lateral points, u(j+/-1) _RL uRt1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL uRt2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C coefficients for lateral points, v(i+/-1) _RL vRt1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL vRt2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) CEOP #ifdef SEAICE_CGRID #ifdef SEAICE_ALLOW_DYNAMICS C !LOCAL VARIABLES: C === Local variables === C i,j,bi,bj :: Loop counters INTEGER i, j, bi, bj _RL hFacM, hFacP C backward difference extrapolation factor (includes 1/deltaT) _RL bdfAlphaOverDt C strongly implicit coupling factor _RL strImpCplFac C fractional area at velocity points _RL areaW(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL areaS(1-OLx:sNx+OLx,1-OLy:sNy+OLy) C coefficients of ice velocities in coefficient matrix C for both U and V-equation C XX: double derivative in X C YY: double derivative in Y C XM: metric term with derivative in X C YM: metric term with derivative in Y _RL UXX (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL UYY (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL UXM (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL UYM (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL VXX (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL VYY (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL VXM (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL VYM (1-OLx:sNx+OLx,1-OLy:sNy+OLy) C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C backward difference extrapolation factor bdfAlphaOverDt = 1. _d 0 IF ( SEAICEuseBDF2 ) THEN IF ( myIter.EQ.nIter0 .AND. SEAICEmomStartBDF.EQ.0 ) THEN bdfAlphaOverDt = 1. _d 0 ELSE bdfAlphaOverDt = 1.5 _d 0 ENDIF ENDIF C bdfAlphaOverDt = bdfAlphaOverDt / SEAICE_deltaTdyn C strImpCplFac = 0. _d 0 IF ( SEAICEuseStrImpCpl ) strImpCplFac = 1. _d 0 DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) IF ( SEAICEscaleSurfStress ) THEN DO J=jMin,jMax DO I=iMin,iMax areaW(I,J) = 0.5 _d 0*(AREA(I,J,bi,bj)+AREA(I-1,J,bi,bj)) areaS(I,J) = 0.5 _d 0*(AREA(I,J,bi,bj)+AREA(I,J-1,bi,bj)) ENDDO ENDDO ELSE DO J=jMin,jMax DO I=iMin,iMax areaW(I,J) = 1. _d 0 areaS(I,J) = 1. _d 0 ENDDO ENDDO ENDIF C C coefficients of uIce(I,J) and vIce(I,J) belonging to ... C DO J=jMin,jMax DO I=iMin-1,iMax C ... d/dx (eta+zeta) d/dx u UXX(I,J) = _dyF(I,J,bi,bj) * etaPlusZeta(I,J,bi,bj) & * _recip_dxF(I,J,bi,bj) C ... d/dx (zeta-eta) k1 u UXM(I,J) = _dyF(I,J,bi,bj) * zetaMinusEta(I,J,bi,bj) & * k1AtC(I,J,bi,bj) * 0.5 _d 0 ENDDO ENDDO DO J=jMin,jMax+1 DO I=iMin,iMax C ... d/dy eta d/dy u UYY(I,J) = _dxV(I,J,bi,bj) * & ( etaZloc(I,J,bi,bj) + strImpCplFac*zetaZloc(I,J,bi,bj) ) & * _recip_dyU(I,J,bi,bj) C ... d/dy eta k2 u UYM(I,J) = _dxV(I,J,bi,bj) * etaZloc(I,J,bi,bj) & * k2AtZ(I,J,bi,bj) * 0.5 _d 0 ENDDO ENDDO DO J=jMin,jMax DO I=iMin,iMax+1 C ... d/dx eta dv/dx VXX(I,J) = _dyU(I,J,bi,bj) * & ( etaZloc(I,J,bi,bj) + strImpCplFac*zetaZloc(I,J,bi,bj) ) & * _recip_dxV(I,J,bi,bj) C ... d/dx eta k1 v VXM(I,J) = _dyU(I,J,bi,bj) * etaZloc(I,J,bi,bj) & * k1AtZ(I,J,bi,bj) * 0.5 _d 0 ENDDO ENDDO DO J=jMin-1,jMax DO I=iMin,iMax C ... d/dy eta+zeta dv/dy VYY(I,J) = _dxF(I,J,bi,bj) * etaPlusZeta(I,J,bi,bj) & * _recip_dyF(I,J,bi,bj) C ... d/dy (zeta-eta) k2 v VYM(I,J) = _dxF(I,J,bi,bj) * zetaMinusEta(I,J,bi,bj) & * k2AtC(I,J,bi,bj) * 0.5 _d 0 ENDDO ENDDO C-- Coefficients for solving uIce : C assemble coefficient matrix, beware of sign convention: because this C is the left hand side we calculate -grad(sigma), but the coefficients C of U(I,J+/-1) are counted on the right hand side DO J=jMin,jMax DO I=iMin,iMax C coefficients for UICE(I-1,J) AU(I,J,bi,bj)= ( - UXX(I-1,J) + UXM(I-1,J) ) & * seaiceMaskU(I,J,bi,bj) C coefficients for UICE(I+1,J) CU(I,J,bi,bj)= ( - UXX(I ,J) - UXM(I ,J) ) & * seaiceMaskU(I,J,bi,bj) C coefficients for UICE(I,J) BU(I,J,bi,bj)=(ONE - seaiceMaskU(I,J,bi,bj)) + & ( UXX(I-1,J) + UXX(I,J) + UYY(I,J+1) + UYY(I,J) & + UXM(I-1,J) - UXM(I,J) + UYM(I,J+1) - UYM(I,J) & ) * seaiceMaskU(I,J,bi,bj) C coefficients of uIce(I,J-1) uRt1(I,J,bi,bj)= UYY(I,J ) + UYM(I,J ) C coefficients of uIce(I,J+1) uRt2(I,J,bi,bj)= UYY(I,J+1) - UYM(I,J+1) ENDDO ENDDO C apply boundary conditions according to slip factor C for no slip, set u on boundary to zero: u(j+/-1)=-u(j) C for the free slip case sigma_12 = 0 DO J=jMin,jMax DO I=iMin,iMax hFacM = seaiceMaskU(I,J-1,bi,bj) hFacP = seaiceMaskU(I,J+1,bi,bj) C copy contributions to coefficient of U(I,J) C beware of sign convection: uRt1/2 have the opposite sign convention C than BU, hence the minus sign BU(I,J,bi,bj)=BU(I,J,bi,bj) + seaiceMaskU(I,J,bi,bj) * & ( ( 1. _d 0 - hFacM ) * ( UYY(I ,J ) + UYM(I ,J ) ) & + ( 1. _d 0 - hFacP ) * ( UYY(I ,J+1) - UYM(I ,J+1) ) ) C reset coefficients of U(I,J-1) and U(I,J+1) uRt1(I,J,bi,bj) = uRt1(I,J,bi,bj) * hFacM uRt2(I,J,bi,bj) = uRt2(I,J,bi,bj) * hFacP ENDDO ENDDO C now we need to normalize everything by the grid cell area DO J=jMin,jMax DO I=iMin,iMax AU(I,J,bi,bj) = AU(I,J,bi,bj) * recip_rAw(I,J,bi,bj) CU(I,J,bi,bj) = CU(I,J,bi,bj) * recip_rAw(I,J,bi,bj) C here we need to add the contribution from the time derivative (in C bdfAlphaOverDt) and the symmetric drag term; must be done after C normalizing BU(I,J,bi,bj) = BU(I,J,bi,bj) * recip_rAw(I,J,bi,bj) & + seaiceMaskU(I,J,bi,bj) * & ( bdfAlphaOverDt*seaiceMassU(I,J,bi,bj) & + 0.5 _d 0 * ( dragSym(I, J,bi,bj) & + dragSym(I-1,J,bi,bj) )*areaW(I,J) & ) uRt1(I,J,bi,bj) = uRt1(I,J,bi,bj) * recip_rAw(I,J,bi,bj) uRt2(I,J,bi,bj) = uRt2(I,J,bi,bj) * recip_rAw(I,J,bi,bj) ENDDO ENDDO C-- Coefficients for solving uIce : C assemble coefficient matrix, beware of sign convention: because this C is the left hand side we calculate -grad(sigma), but the coefficients C of V(I+/-1,J) are counted on the right hand side DO J=jMin,jMax DO I=iMin,iMax C coefficients for VICE(I,J-1) AV(I,J,bi,bj)=( - VYY(I,J-1) + VYM(I,J-1) & ) * seaiceMaskV(I,J,bi,bj) C coefficients for VICE(I,J+1) CV(I,J,bi,bj)=( - VYY(I,J ) - VYM(I,J ) & ) * seaiceMaskV(I,J,bi,bj) C coefficients for VICE(I,J) BV(I,J,bi,bj)= (ONE - seaiceMaskV(I,J,bi,bj)) + & ( VXX(I,J) + VXX(I+1,J) + VYY(I,J) + VYY(I,J-1) & - VXM(I,J) + VXM(I+1,J) - VYM(I,J) + VYM(I,J-1) & ) * seaiceMaskV(I,J,bi,bj) C coefficients for V(I-1,J) vRt1(I,J,bi,bj) = VXX(I ,J) + VXM(I ,J) C coefficients for V(I+1,J) vRt2(I,J,bi,bj) = VXX(I+1,J) - VXM(I+1,J) ENDDO ENDDO C apply boundary conditions according to slip factor C for no slip, set u on boundary to zero: v(i+/-1)=-v(i) C for the free slip case sigma_12 = 0 DO J=jMin,jMax DO I=iMin,iMax hFacM = seaiceMaskV(i-1,j,bi,bj) hFacP = seaiceMaskV(i+1,j,bi,bj) C copy contributions to coefficient of V(I,J) C beware of sign convection: vRt1/2 have the opposite sign convention C than BV, hence the minus sign BV(I,J,bi,bj)=BV(I,J,bi,bj) + seaiceMaskV(I,J,bi,bj) * & ( ( 1. _d 0 - hFacM ) * ( VXX(I ,J) + VXM(I ,J) ) & + ( 1. _d 0 - hFacP ) * ( VXX(I+1,J) - VXM(I+1,J) ) ) C reset coefficients of V(I-1,J) and V(I+1,J) vRt1(I,J,bi,bj) = vRt1(I,J,bi,bj) * hFacM vRt2(I,J,bi,bj) = vRt2(I,J,bi,bj) * hFacP ENDDO ENDDO C now we need to normalize everything by the grid cell area DO J=jMin,jMax DO I=iMin,iMax AV(I,J,bi,bj) = AV(I,J,bi,bj) * recip_rAs(I,J,bi,bj) CV(I,J,bi,bj) = CV(I,J,bi,bj) * recip_rAs(I,J,bi,bj) C here we need add the contribution from the time derivative (in C bdfAlphaOverDt) and the symmetric drag term; must be done after C normalizing BV(I,J,bi,bj) = BV(I,J,bi,bj) * recip_rAs(I,J,bi,bj) & + seaiceMaskV(I,J,bi,bj) * & ( bdfAlphaOverDt*seaiceMassV(I,J,bi,bj) & + 0.5 _d 0 * ( dragSym(I,J, bi,bj) & + dragSym(I,J-1,bi,bj) )*areaS(I,J) & ) vRt1(I,J,bi,bj) = vRt1(I,J,bi,bj) * recip_rAs(I,J,bi,bj) vRt2(I,J,bi,bj) = vRt2(I,J,bi,bj) * recip_rAs(I,J,bi,bj) ENDDO ENDDO IF ( ( useCubedSphereExchange .AND. & (SEAICE_OLx.GT.0 .OR. SEAICE_OLy.GT.0) ) .OR. & SEAICEscaleSurfStress ) THEN C this is clearly a hack: make sure that diagonals BU/BV are non-zero. C In my experience we only need to catch the case of SEAICE_OLx/y=OLx/y-2, C but for safety lets call this routine whenever SEAICE_OLx/y > 0 C & (SEAICE_OLx.EQ.OLx-2 .OR. SEAICE_OLy.EQ.OLy-2) ) THEN C When scaling the surface ice-ocean stress by AREA, then there will C be many cases of zero diagonal entries. DO J=jMin,jMax DO I=iMin,iMax IF (BU(I,J,bi,bj).EQ.0 _d 0) BU(I,J,bi,bj) = 1. _d 0 IF (BV(I,J,bi,bj).EQ.0 _d 0) BV(I,J,bi,bj) = 1. _d 0 ENDDO ENDDO ENDIF C bi/bj-loops ENDDO ENDDO RETURN END


C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP C !ROUTINE: SEAICE_LSR_RHSU C !INTERFACE: SUBROUTINE SEAICE_LSR_RHSU( I zetaMinusEta, etaPlusZeta, etaZloc, zetaZloc, pressLoc, I uIceC, vIceC, U rhsU, I iMin, iMax, jMin, jMax, bi, bj, myThid ) C !DESCRIPTION: \bv C *==========================================================* C | S/R SEAICE_LSR_RHSU C | o Set up stress contribution to rhs of u-momentum eq. C *==========================================================* C \ev C !USES: IMPLICIT NONE #include "SIZE.h" #include "EEPARAMS.h" #include "GRID.h" #include "SEAICE_SIZE.h" #include "SEAICE_PARAMS.h" #include "SEAICE.h" C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === C myThid :: my Thread Id. number INTEGER myThid INTEGER iMin, iMax, jMin, jMax, bi, bj C _RL zetaMinusEta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL etaPlusZeta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL etaZloc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL zetaZloc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL pressloc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL uIceC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL vIceC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL rhsU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C !LOCAL VARIABLES: C === Local variables === C i,j :: Loop counters INTEGER I,J INTEGER kSrf _RL hFacM _RL sig11(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL sig12(1-OLx:sNx+OLx,1-OLy:sNy+OLy) CEOP kSrf = 1 DO J=1-OLy,sNy+OLy DO I=1-OLx,sNx+OLx sig11(I,J) = 0. _d 0 sig12(I,J) = 0. _d 0 ENDDO ENDDO C contribution of sigma11 to rhs DO J=jMin,jMax DO I=iMin-1,iMax sig11(I,J) = zetaMinusEta(I,J,bi,bj) & * ( vIceC(I,J+1,bi,bj) - vIceC(I,J,bi,bj) ) & * _recip_dyF(I,J,bi,bj) & + etaPlusZeta(I,J,bi,bj) * k2AtC(I,J,bi,bj) & * 0.5 _d 0 * ( vIceC(I,J+1,bi,bj) + vIceC(I,J,bi,bj) ) & - 0.5 _d 0 * pressLoc(I,J,bi,bj) ENDDO ENDDO C contribution of sigma12 to rhs of u-equation DO J=jMin,jMax+1 DO I=iMin,iMax hFacM = seaiceMaskV(I,J,bi,bj) - seaiceMaskV(I-1,J,bi,bj) sig12(I,J) = etaZloc(I,J,bi,bj) * ( & ( vIceC(I,J,bi,bj) - vIceC(I-1,J,bi,bj) ) & * _recip_dxV(I,J,bi,bj) & - k1AtZ(I,J,bi,bj) & * 0.5 _d 0 * ( vIceC(I,J,bi,bj) + vIceC(I-1,J,bi,bj) ) & ) C free slip conditions (sig12=0) are taken care of by masking sig12 & *maskC(I ,J ,kSrf,bi,bj)*maskC(I-1,J ,kSrf,bi,bj) & *maskC(I ,J-1,kSrf,bi,bj)*maskC(I-1,J-1,kSrf,bi,bj) C no slip boundary conditions (v(i-1)=-v(i)) C v(i)+v(i-1) = 0 is also taken care of by masking sig12, so that we C only need to deal with v(i)-v(i-1) & + etaZloc(I,J,bi,bj) * _recip_dxV(I,J,bi,bj) & * ( vIceC(I,J,bi,bj) + vIceC(I-1,J,bi,bj) ) & * hFacM * 2. _d 0 ENDDO ENDDO IF ( SEAICEuseStrImpCpl ) THEN C strictly speaking, this is not a contribution to sig12, but for C convenience we add the explicit term -zetaZ*du/dy here; C we do not have to add any metric terms, because they cancel. DO J=jMin,jMax+1 DO I=iMin,iMax hFacM = seaiceMaskV(I,J,bi,bj) - seaiceMaskV(I-1,J,bi,bj) sig12(I,J) = sig12(I,J) - zetaZloc(I,J,bi,bj) * ( & ( uIceC(I,J,bi,bj) - uIceC(I,J-1,bi,bj) ) & * _recip_dyU(I,J,bi,bj) & ) C free slip conditions (sig12=0) are taken care of by masking sig12 & *maskC(I ,J ,kSrf,bi,bj)*maskC(I-1,J ,kSrf,bi,bj) & *maskC(I ,J-1,kSrf,bi,bj)*maskC(I-1,J-1,kSrf,bi,bj) C no slip boundary conditions (u(j-1)=-u(j)) C u(j)+u(j-1) = 0 is also taken care of by masking sig12, so that we C only need to deal with u(j)-u(j-1) & - zetaZloc(I,J,bi,bj) * _recip_dyU(I,J,bi,bj) & * ( uIceC(I,J,bi,bj) + uIceC(I,J-1,bi,bj) ) & * hFacM * 2. _d 0 ENDDO ENDDO ENDIF DO J=jMin,jMax DO I=iMin,iMax C contribution to the rhs part of grad(sigma)_x rhsU(I,J,bi,bj) = rhsU(I,J,bi,bj) & + recip_rAw(I,J,bi,bj) * seaiceMaskU(I,J,bi,bj) * & ( _dyF(I ,J ,bi,bj)*sig11(I ,J ) & - _dyF(I-1,J ,bi,bj)*sig11(I-1,J ) & + _dxV(I ,J+1,bi,bj)*sig12(I ,J+1) & - _dxV(I ,J ,bi,bj)*sig12(I ,J ) ) ENDDO ENDDO RETURN END


C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP C !ROUTINE: SEAICE_LSR_RHSV C !INTERFACE: SUBROUTINE SEAICE_LSR_RHSV( I zetaMinusEta, etaPlusZeta, etaZloc, zetaZloc, pressLoc, I uIceC, vIceC, U rhsV, I iMin, iMax, jMin, jMax, bi, bj, myThid ) C !DESCRIPTION: \bv C *==========================================================* C | S/R SEAICE_LSR_RHSV C | o Set up stress contribution to rhs of v-momentum eq. C *==========================================================* C \ev C !USES: IMPLICIT NONE #include "SIZE.h" #include "EEPARAMS.h" #include "GRID.h" #include "SEAICE_SIZE.h" #include "SEAICE_PARAMS.h" #include "SEAICE.h" C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === C myThid :: my Thread Id. number INTEGER myThid INTEGER iMin, iMax, jMin, jMax, bi, bj C _RL zetaMinusEta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL etaPlusZeta (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL etaZloc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL zetaZloc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL pressloc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL uIceC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL vIceC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL rhsV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C !LOCAL VARIABLES: C === Local variables === C i,j :: Loop counters INTEGER I,J INTEGER kSrf _RL hFacM _RL sig22(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL sig12(1-OLx:sNx+OLx,1-OLy:sNy+OLy) CEOP kSrf = 1 DO J = 1-Oly,sNy+Oly DO I = 1-Olx,sNx+Olx sig22(I,J) = 0. _d 0 sig12(I,J) = 0. _d 0 ENDDO ENDDO C contribution of sigma22 to rhs DO J=jMin-1,jMax DO I=iMin,iMax sig22(I,J) = zetaMinusEta(I,J,bi,bj) & * ( uIceC(I+1,J,bi,bj) - uIceC(I,J,bi,bj) ) & * _recip_dxF(I,J,bi,bj) & + etaPlusZeta(I,J,bi,bj) * k1AtC(I,J,bi,bj) & * 0.5 _d 0 * ( uIceC(I+1,J,bi,bj) + uIceC(I,J,bi,bj) ) & - 0.5 _d 0 * pressLoc(I,J,bi,bj) ENDDO ENDDO C contribution of sigma12 to rhs of v-equation DO J=jMin,jMax DO I=iMin,iMax+1 hFacM = seaiceMaskU(i,j,bi,bj) - seaiceMaskU(i,j-1,bi,bj) sig12(I,J) = etaZloc(I,J,bi,bj) * ( & ( uIceC(I,J,bi,bj) - uIceC(I,J-1,bi,bj) ) & * _recip_dyU(I,J,bi,bj) & - k2AtZ(I,J,bi,bj) & * 0.5 _d 0 * ( uIceC(I,J,bi,bj) + uIceC(I,J-1,bi,bj) ) & ) C free slip conditions (sig12=0) are taken care of by masking sig12, & *maskC(I ,J ,kSrf,bi,bj)*maskC(I-1,J ,kSrf,bi,bj) & *maskC(I ,J-1,kSrf,bi,bj)*maskC(I-1,J-1,kSrf,bi,bj) C no slip boundary conditions (u(j-1)=-u(j)) C u(j)+u(j-1) = 0 is also taken care of by masking sig12, so that we C only need to deal with u(j)-u(j-1) & + etaZloc(I,J,bi,bj) * _recip_dyU(I,J,bi,bj) & * ( uIceC(I,J,bi,bj) + uIceC(I,J-1,bi,bj) ) & * hFacM * 2. _d 0 ENDDO ENDDO IF ( SEAICEuseStrImpCpl ) THEN C strictly speaking, this is not a contribution to sig12, but for C convenience we add the explicit term -zetaZ*dv/dx here; C we do not have to add any metric terms, because they cancel. DO J=jMin,jMax DO I=iMin,iMax+1 hFacM = seaiceMaskU(i,j,bi,bj) - seaiceMaskU(i,j-1,bi,bj) sig12(I,J) = sig12(I,J) - zetaZloc(I,J,bi,bj) * ( & ( vIceC(I,J,bi,bj) - vIceC(I-1,J,bi,bj) ) & * _recip_dxV(I,J,bi,bj) & ) C free slip conditions (sig12=0) are taken care of by masking sig12 & *maskC(I ,J ,kSrf,bi,bj)*maskC(I-1,J ,kSrf,bi,bj) & *maskC(I ,J-1,kSrf,bi,bj)*maskC(I-1,J-1,kSrf,bi,bj) C no slip boundary conditions (v(i-1)=-v(i)) C v(i)+v(i-1) = 0 is also taken care of by masking sig12, so that we C only need to deal with v(i)-v(i-1) & - zetaZloc(I,J,bi,bj) * _recip_dxV(I,J,bi,bj) & * ( vIceC(I,J,bi,bj) + vIceC(I-1,J,bi,bj) ) & * hFacM * 2. _d 0 ENDDO ENDDO ENDIF DO J=jMin,jMax DO I=iMin,iMax C contribution to the rhs part of grad(sigma)_y rhsV(I,J,bi,bj) = rhsV(I,J,bi,bj) & + recip_rAs(I,J,bi,bj) * seaiceMaskV(I,J,bi,bj) * & ( _dyU(I+1,J ,bi,bj) * sig12(I+1,J ) & - _dyU(I ,J ,bi,bj) * sig12(I ,J ) & + _dxF(I ,J ,bi,bj) * sig22(I ,J ) & - _dxF(I ,J-1,bi,bj) * sig22(I ,J-1) ) ENDDO ENDDO RETURN END


C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP C !ROUTINE: SEAICE_LSR_TRIDIAGU C !INTERFACE: SUBROUTINE SEAICE_LSR_TRIDIAGU( I AU, BU, CU, uRt1, uRt2, rhsU, uTmp, seaiceMaskU, WFAU, U uIce, I iMin, iMax, jMin, jMax, bi, bj, myTime, myIter, myThid ) C !DESCRIPTION: \bv C *==========================================================* C | S/R SEAICE_LSR_TRIDIAGU C | o Solve tridiagonal problem in uIce for LSR solver C *==========================================================* C \ev C !USES: IMPLICIT NONE #include "SIZE.h" C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === C myTime :: Simulation time C myIter :: Simulation timestep number C myThid :: my Thread Id. number _RL myTime INTEGER myIter INTEGER myThid INTEGER iMin, iMax, jMin, jMax, bi, bj C diagonals of coefficient matrices _RL AU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL BU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL CU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C coefficients for lateral points, v(j+/-1) _RL uRt1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL uRt2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C right hand side _RL rhsU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C velocity of previous iteration step _RL uTmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C on entry: first guess and on exit: solution _RL uIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C mask _RL seaiceMaskU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C relaxation factor _RL WFAU CEOP C !LOCAL VARIABLES: C === Local variables === C i,j :: Loop counters INTEGER I,J,IM INTEGER jMinLoc, jStep #ifdef SEAICE_LSR_ZEBRA INTEGER K PARAMETER ( jStep = 2 ) #else PARAMETER ( jStep = 1 ) #endif /* SEAICE_LSR_ZEBRA */ _RL AA3, bet _RL URT(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL CUU(1-OLx:sNx+OLx,1-OLy:sNy+OLy) C Need to initialize local arrays DO J = 1-OLy,sNy+OLy DO I = 1-OLx,sNx+OLx URT(I,J) = 0. _d 0 CUU(I,J) = 0. _d 0 ENDDO ENDDO #ifdef SEAICE_LSR_ZEBRA DO K=0,1 jMinLoc = jMin+K #else jMinLoc = jMin #endif /* SEAICE_LSR_ZEBRA */ DO J=jMinLoc,jMax,jStep DO I=iMin,iMax AA3 = 0. _d 0 IF (I.EQ.iMin) AA3 = AA3 - AU(I,J,bi,bj)*uIce(I-1,J,bi,bj) IF (I.EQ.iMax) AA3 = AA3 - CU(I,J,bi,bj)*uIce(I+1,J,bi,bj) URT(I,J)=rhsU(I,J,bi,bj) & + AA3 & + uRt1(I,J,bi,bj)*uIce(I,J-1,bi,bj) & + uRt2(I,J,bi,bj)*uIce(I,J+1,bi,bj) URT(I,J)=URT(I,J)* seaiceMaskU(I,J,bi,bj) ENDDO C beginning of tridiagnoal solver DO I=iMin,iMax CUU(I,J)=CU(I,J,bi,bj) ENDDO CUU(iMin,J)=CUU(iMin,J)/BU(iMin,J,bi,bj) URT(iMin,J)=URT(iMin,J)/BU(iMin,J,bi,bj) #ifdef SEAICE_VECTORIZE_LSR ENDDO C start a new loop with reversed order to support automatic vectorization CADJ loop = sequential DO I=iMin+1,iMax IM=I-1 DO J=jMinLoc,jMax,jStep #else /* do not SEAICE_VECTORIZE_LSR */ CADJ loop = sequential DO I=iMin+1,iMax IM=I-1 #endif /* SEAICE_VECTORIZE_LSR */ bet = BU(I,J,bi,bj)-AU(I,J,bi,bj)*CUU(IM,J) CUU(I,J) = CUU(I,J)/bet URT(I,J) = (URT(I,J)-AU(I,J,bi,bj)*URT(IM,J))/bet ENDDO #ifdef SEAICE_VECTORIZE_LSR ENDDO CADJ loop = sequential DO I=iMin,iMax-1 IM=sNx-I DO J=jMinLoc,jMax,jStep #else CADJ loop = sequential DO I=iMin,iMax-1 IM=sNx-I #endif /* SEAICE_VECTORIZE_LSR */ URT(IM,J)=URT(IM,J)-CUU(IM,J)*URT(IM+1,J) ENDDO #ifdef SEAICE_VECTORIZE_LSR ENDDO #endif /* SEAICE_VECTORIZE_LSR */ C end of tridiagnoal solver, solution is URT C relaxation with WFAU #ifdef SEAICE_VECTORIZE_LSR C go back to original order DO J=jMinLoc,jMax,jStep #endif /* SEAICE_VECTORIZE_LSR */ DO I=iMin,iMax uIce(I,J,bi,bj)=uTmp(I,J,bi,bj) & +WFAU*(URT(I,J)-uTmp(I,J,bi,bj)) ENDDO ENDDO #ifdef SEAICE_LSR_ZEBRA ENDDO #endif /* SEAICE_LSR_ZEBRA */ RETURN END


C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP C !ROUTINE: SEAICE_LSR_TRIDIAGV C !INTERFACE: SUBROUTINE SEAICE_LSR_TRIDIAGV( I AV, BV, CV, vRt1, vRt2, rhsV, vTmp, seaiceMaskV, WFAV, U vIce, I iMin, iMax, jMin, jMax, bi, bj, myTime, myIter, myThid ) C !DESCRIPTION: \bv C *==========================================================* C | S/R SEAICE_LSR_TRIDIAGV C | o Solve tridiagonal problem in vIce for LSR solver C *==========================================================* C \ev C !USES: IMPLICIT NONE #include "SIZE.h" C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === C myTime :: Simulation time C myIter :: Simulation timestep number C myThid :: my Thread Id. number _RL myTime INTEGER myIter INTEGER myThid INTEGER iMin, iMax, jMin, jMax, bi, bj C diagonals of coefficient matrices _RL AV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL BV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL CV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C coefficients for lateral points, v(j+/-1) _RL vRt1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL vRt2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C right hand side _RL rhsV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C velocity of previous iteration step _RL vTmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C on entry: first guess and on exit: solution _RL vIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C mask _RL seaiceMaskV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) C relaxation factor _RL WFAV CEOP C !LOCAL VARIABLES: C === Local variables === C i,j :: Loop counters INTEGER I,J,JM INTEGER iMinLoc, iStep #ifdef SEAICE_LSR_ZEBRA INTEGER K PARAMETER ( iStep = 2 ) #else PARAMETER ( iStep = 1 ) #endif /* SEAICE_LSR_ZEBRA */ _RL AA3, bet _RL VRT(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL CVV(1-OLx:sNx+OLx,1-OLy:sNy+OLy) C Need to initialize local arrays DO J = 1-OLy,sNy+OLy DO I = 1-OLx,sNx+OLx VRT(I,J) = 0. _d 0 CVV(I,J) = 0. _d 0 ENDDO ENDDO #ifdef SEAICE_LSR_ZEBRA DO K=0,1 iMinLoc = iMin+K #else iMinLoc = iMin #endif /* SEAICE_LSR_ZEBRA */ #ifdef SEAICE_VECTORIZE_LSR DO J=jMin,jMax DO I=iMinLoc,iMax,iStep #else DO I=iMinLoc,iMax,iStep DO J=jMin,jMax #endif /* SEAICE_VECTORIZE_LSR */ AA3 = 0. _d 0 IF (J.EQ.jMin) AA3 = AA3 - AV(I,J,bi,bj)*vIce(I,J-1,bi,bj) IF (J.EQ.jMax) AA3 = AA3 - CV(I,J,bi,bj)*vIce(I,J+1,bi,bj) VRT(I,J)=rhsV(I,J,bi,bj) & + AA3 & + vRt1(I,J,bi,bj)*vIce(I-1,J,bi,bj) & + vRt2(I,J,bi,bj)*vIce(I+1,J,bi,bj) VRT(I,J)=VRT(I,J)* seaiceMaskV(I,J,bi,bj) C beginning of tridiagnoal solver CVV(I,J)=CV(I,J,bi,bj) ENDDO #ifdef SEAICE_VECTORIZE_LSR ENDDO DO I=iMinLoc,iMax,iStep #endif /* SEAICE_VECTORIZE_LSR */ CVV(I,jMin)=CVV(I,jMin)/BV(I,jMin,bi,bj) VRT(I,jMin)=VRT(I,jMin)/BV(I,jMin,bi,bj) #ifdef SEAICE_VECTORIZE_LSR ENDDO C start a new loop with reversed order to support automatic vectorization CADJ loop = sequential DO J=jMin+1,jMax JM=J-1 DO I=iMinLoc,iMax,iStep #else /* do not SEAICE_VECTORIZE_LSR */ CADJ loop = sequential DO J=jMin+1,jMax JM=J-1 #endif /* SEAICE_VECTORIZE_LSR */ bet = BV(I,J,bi,bj)-AV(I,J,bi,bj)*CVV(I,JM) CVV(I,J) = CVV(I,J)/bet VRT(I,J) = (VRT(I,J)-AV(I,J,bi,bj)*VRT(I,JM))/bet ENDDO #ifdef SEAICE_VECTORIZE_LSR ENDDO C go back to original order CADJ loop = sequential DO J=jMin,jMax-1 JM=sNy-J DO I=iMinLoc,iMax,iStep #else CADJ loop = sequential DO J=jMin,jMax-1 JM=sNy-J #endif /* SEAICE_VECTORIZE_LSR */ VRT(I,JM)=VRT(I,JM)-CVV(I,JM)*VRT(I,JM+1) ENDDO #ifdef SEAICE_VECTORIZE_LSR ENDDO C end of tridiagnoal solver, solution is VRT C relaxation with WFAV DO J=jMin,jMax DO I=iMinLoc,iMax,iStep #else DO J=jMin,jMax #endif /* SEAICE_VECTORIZE_LSR */ vIce(I,J,bi,bj)=vTmp(I,J,bi,bj) & +WFAV*(VRT(I,J)-vTmp(I,J,bi,bj)) ENDDO ENDDO #ifdef SEAICE_LSR_ZEBRA ENDDO #endif /* SEAICE_LSR_ZEBRA */ #endif /* SEAICE_ALLOW_DYNAMICS */ #endif /* SEAICE_CGRID */ RETURN END