C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_calc_residual.F,v 1.9 2016/05/17 15:26:46 mlosch Exp $
C $Name:  $

#include "SEAICE_OPTIONS.h"

CBOP
C     !ROUTINE: SEAICE_CALC_RESIDUAL
C     !INTERFACE:
      SUBROUTINE SEAICE_CALC_RESIDUAL( 
     I     uIceLoc, vIceLoc, 
     O     uIceRes, vIceRes, 
     I     newtonIter, krylovIter, myTime, myIter, myThid )

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE SEAICE_CALC_RESIDUAL
C     | o For Jacobian-free Newton-Krylov solver compute
C     |   the residual of the momentum equations
C     *==========================================================*
C     | written by Martin Losch, Oct 2012
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "SEAICE_SIZE.h"
#ifdef SEAICE_ALLOW_BOTTOMDRAG
#include "SEAICE_PARAMS.h"
#endif /* SEAICE_ALLOW_BOTTOMDRAG */ 
#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
C     newtonIter :: current iterate of Newton iteration
C     krylovIter :: current iterate of Krylov iteration
      _RL     myTime
      INTEGER myIter
      INTEGER myThid
      INTEGER newtonIter
      INTEGER krylovIter
C     u/vIceLoc :: local copies of the current ice velocity
      _RL uIceLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL vIceLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
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)

#ifdef SEAICE_ALLOW_JFNK
C     u/vIceLHS :: left hand side of momentum equations 
      _RL uIceLHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL vIceLHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
C     u/vIceRHS :: righ hand side of momentum equations
      _RL uIceRHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL vIceRHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)

C     i,j,bi,bj :: loop indices
      INTEGER i,j,bi,bj
CEOP

C     Initialise
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO J=1-Oly,sNy+Oly
         DO I=1-Olx,sNx+Olx
          uIceLHS(I,J,bi,bj) = 0. _d 0
          vIceLHS(I,J,bi,bj) = 0. _d 0
          uIceRHS(I,J,bi,bj) = 0. _d 0
          vIceRHS(I,J,bi,bj) = 0. _d 0
         ENDDO
        ENDDO
       ENDDO
      ENDDO
C     u/vIceLoc have changed so that new drag coefficients and
C     viscosities are required
      CALL SEAICE_OCEANDRAG_COEFFS(
     I     uIceLoc, vIceLoc,
     O     DWATN,
     I     krylovIter, myTime, myIter, myThid )
#ifdef SEAICE_ALLOW_BOTTOMDRAG
      IF (SEAICEbasalDragK2.GT.0. _d 0) CALL SEAICE_BOTTOMDRAG_COEFFS(
     I     uIceLoc, vIceLoc,
#ifdef SEAICE_ITD
     I     HEFFITD, AREAITD, AREA,
#else
     I     HEFF, AREA,
#endif      
     O     CbotC,
     I     krylovIter, myTime, myIter, myThid )
#endif /* SEAICE_ALLOW_BOTTOMDRAG */
      CALL SEAICE_CALC_STRAINRATES(
     I     uIceLoc, vIceLoc,
     O     e11, e22, e12,
     I     krylovIter, myTime, myIter, myThid )
      CALL SEAICE_CALC_VISCOSITIES(
     I     e11, e22, e12, zMin, zMax, hEffM, press0, tensileStrFac,
     O     eta, etaZ, zeta, zetaZ, press, deltaC,
     I     krylovIter, myTime, myIter, myThid )

C     The scheme is backward Euler in time, i.e. the rhs-vector contains 
C     only terms that are independent of u/vIce, except for the time 
C     derivative part mass*(u/vIce-u/vIceNm1)/deltaT

C     compute new right hand side (depends to DWATN=Cdrag)
C     sea-surface tilt and wind stress: FORCEX0, FORCEY0
C     + mass*(u/vIceNm1)/deltaT
C     + Cdrag*(uVel*cosWat - vVel*sinWat)/(vVel*cosWat + uVel*sinWat)
      CALL SEAICE_CALC_RHS(
     O      uIceRHS, vIceRHS,
     I      newtonIter, krylovIter, myTime, myIter, myThid )

C     Left-hand side contributions:
C     + mass*(u/vIce)/deltaT
C     + Cdrag*(uIce*cosWat - vIce*sinWat)/(vIce*cosWat + uIce*sinWat)
C     + CdragBot*uIce/vIce
C     - mass*f*vIce/+mass*f*uIce
C     - dsigma/dx / -dsigma/dy, eta and zeta are only computed once per
C     Newton iterate
       CALL SEAICE_CALC_LHS(
     I      uIceLoc, vIceLoc,
     O      uIceLHS, vIceLHS,
     I      newtonIter, myTime, myIter, myThid )

C     Right-hand side contributions only need to be computed once per
C     time step, therefore we will put them into a separate routine
C     and call them elsewhere to save floating point operations

C     Calculate the residual
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         DO J=1,sNy
          DO I=1,sNx
           uIceRes(I,J,bi,bj) = uIceLHS(I,J,bi,bj) - uIceRHS(I,J,bi,bj)
           vIceRes(I,J,bi,bj) = vIceLHS(I,J,bi,bj) - vIceRHS(I,J,bi,bj)
          ENDDO
         ENDDO
        ENDDO
       ENDDO

#endif /* SEAICE_ALLOW_JFNK */

      RETURN
      END