pkg/seaice/seaice_calc_residual.F

2fd913c523 Mart*0001 #include "SEAICE_OPTIONS.h"
                0002
                0003 CBOP
                0004 C     !ROUTINE: SEAICE_CALC_RESIDUAL
                0005 C     !INTERFACE:
ec0d7df165 Mart*0006       SUBROUTINE SEAICE_CALC_RESIDUAL(
                0007      I     uIceLoc, vIceLoc,
                0008      O     uIceRes, vIceRes,
2fd913c523 Mart*0009      I     newtonIter, krylovIter, myTime, myIter, myThid )
                0010
                0011 C     !DESCRIPTION: \bv
                0012 C     *==========================================================*
                0013 C     | SUBROUTINE SEAICE_CALC_RESIDUAL
                0014 C     | o For Jacobian-free Newton-Krylov solver compute
                0015 C     |   the residual of the momentum equations
                0016 C     *==========================================================*
                0017 C     | written by Martin Losch, Oct 2012
                0018 C     *==========================================================*
                0019 C     \ev
                0020
                0021 C     !USES:
                0022       IMPLICIT NONE
                0023
                0024 C     === Global variables ===
                0025 #include "SIZE.h"
                0026 #include "EEPARAMS.h"
                0027 #include "SEAICE_SIZE.h"
5bb179ddc2 Mart*0028 #if ( defined SEAICE_ALLOW_BOTTOMDRAG || defined SEAICE_ALLOW_SIDEDRAG )
                0029 # include "SEAICE_PARAMS.h"
                0030 #endif
3f0f10fc37 Mart*0031 #include "SEAICE_GRID.h"
2fd913c523 Mart*0032 #include "SEAICE.h"
                0033
                0034 C     !INPUT/OUTPUT PARAMETERS:
                0035 C     === Routine arguments ===
                0036 C     myTime :: Simulation time
                0037 C     myIter :: Simulation timestep number
                0038 C     myThid :: my Thread Id. number
                0039 C     newtonIter :: current iterate of Newton iteration
                0040 C     krylovIter :: current iterate of Krylov iteration
                0041       _RL     myTime
                0042       INTEGER myIter
                0043       INTEGER myThid
                0044       INTEGER newtonIter
                0045       INTEGER krylovIter
                0046 C     u/vIceLoc :: local copies of the current ice velocity
                0047       _RL uIceLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0048       _RL vIceLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0049 C     u/vIceRes :: residual of sea-ice momentum equations
                0050       _RL uIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0051       _RL vIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0052
45315406aa Mart*0053 #if ( defined SEAICE_CGRID && defined SEAICE_ALLOW_JFNK )
ec0d7df165 Mart*0054 C     u/vIceLHS :: left hand side of momentum equations
2fd913c523 Mart*0055       _RL uIceLHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0056       _RL vIceLHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0057 C     u/vIceRHS :: righ hand side of momentum equations
                0058       _RL uIceRHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0059       _RL vIceRHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0060
                0061 C     i,j,bi,bj :: loop indices
                0062       INTEGER i,j,bi,bj
                0063 CEOP
                0064
                0065 C     Initialise
                0066       DO bj=myByLo(myThid),myByHi(myThid)
                0067        DO bi=myBxLo(myThid),myBxHi(myThid)
ec0d7df165 Mart*0068         DO J=1-OLy,sNy+OLy
                0069          DO I=1-OLx,sNx+OLx
2fd913c523 Mart*0070           uIceLHS(I,J,bi,bj) = 0. _d 0
                0071           vIceLHS(I,J,bi,bj) = 0. _d 0
                0072           uIceRHS(I,J,bi,bj) = 0. _d 0
                0073           vIceRHS(I,J,bi,bj) = 0. _d 0
                0074          ENDDO
                0075         ENDDO
                0076        ENDDO
                0077       ENDDO
                0078 C     u/vIceLoc have changed so that new drag coefficients and
                0079 C     viscosities are required
                0080       CALL SEAICE_OCEANDRAG_COEFFS(
ec0d7df165 Mart*0081      I     uIceLoc, vIceLoc, HEFFM,
2fd913c523 Mart*0082      O     DWATN,
                0083      I     krylovIter, myTime, myIter, myThid )
df1dac8b7b Mart*0084 #ifdef SEAICE_ALLOW_BOTTOMDRAG
d5254b4e3d Mart*0085       CALL SEAICE_BOTTOMDRAG_COEFFS(
ec0d7df165 Mart*0086      I     uIceLoc, vIceLoc, HEFFM,
df1dac8b7b Mart*0087 #ifdef SEAICE_ITD
                0088      I     HEFFITD, AREAITD, AREA,
                0089 #else
                0090      I     HEFF, AREA,
ec0d7df165 Mart*0091 #endif
df1dac8b7b Mart*0092      O     CbotC,
                0093      I     krylovIter, myTime, myIter, myThid )
                0094 #endif /* SEAICE_ALLOW_BOTTOMDRAG */
5bb179ddc2 Mart*0095 #ifdef SEAICE_ALLOW_SIDEDRAG
                0096         IF ( SEAICEsideDrag .NE. 0. _d 0 ) CALL SEAICE_SIDEDRAG_STRESS(
                0097      I       uIceLoc, vIceLoc, coastRoughU, coastRoughV, AREA,
                0098      O       sideDragU, sideDragV,
                0099      I       krylovIter, myTime, myIter, myThid )
                0100 #endif /* SEAICE_ALLOW_SIDEDRAG */
2fd913c523 Mart*0101       CALL SEAICE_CALC_STRAINRATES(
                0102      I     uIceLoc, vIceLoc,
                0103      O     e11, e22, e12,
                0104      I     krylovIter, myTime, myIter, myThid )
                0105       CALL SEAICE_CALC_VISCOSITIES(
8e32c48b8f Mart*0106      I     e11, e22, e12, SEAICE_zMin, SEAICE_zMax, HEFFM, press0,
                0107      I     tensileStrFac,
0fd6b2de1a Mart*0108      O     eta, etaZ, zeta, zetaZ, press, deltaC,
2fd913c523 Mart*0109      I     krylovIter, myTime, myIter, myThid )
                0110
ec0d7df165 Mart*0111 C     The scheme is backward Euler in time, i.e. the rhs-vector contains
                0112 C     only terms that are independent of u/vIce, except for the time
2fd913c523 Mart*0113 C     derivative part mass*(u/vIce-u/vIceNm1)/deltaT
                0114
                0115 C     compute new right hand side (depends to DWATN=Cdrag)
                0116 C     sea-surface tilt and wind stress: FORCEX0, FORCEY0
                0117 C     + mass*(u/vIceNm1)/deltaT
                0118 C     + Cdrag*(uVel*cosWat - vVel*sinWat)/(vVel*cosWat + uVel*sinWat)
                0119       CALL SEAICE_CALC_RHS(
                0120      O      uIceRHS, vIceRHS,
                0121      I      newtonIter, krylovIter, myTime, myIter, myThid )
                0122
                0123 C     Left-hand side contributions:
                0124 C     + mass*(u/vIce)/deltaT
                0125 C     + Cdrag*(uIce*cosWat - vIce*sinWat)/(vIce*cosWat + uIce*sinWat)
df1dac8b7b Mart*0126 C     + CdragBot*uIce/vIce
2fd913c523 Mart*0127 C     - mass*f*vIce/+mass*f*uIce
                0128 C     - dsigma/dx / -dsigma/dy, eta and zeta are only computed once per
                0129 C     Newton iterate
                0130        CALL SEAICE_CALC_LHS(
                0131      I      uIceLoc, vIceLoc,
                0132      O      uIceLHS, vIceLHS,
                0133      I      newtonIter, myTime, myIter, myThid )
                0134
                0135 C     Right-hand side contributions only need to be computed once per
                0136 C     time step, therefore we will put them into a separate routine
                0137 C     and call them elsewhere to save floating point operations
                0138
                0139 C     Calculate the residual
                0140        DO bj=myByLo(myThid),myByHi(myThid)
                0141         DO bi=myBxLo(myThid),myBxHi(myThid)
                0142          DO J=1,sNy
                0143           DO I=1,sNx
                0144            uIceRes(I,J,bi,bj) = uIceLHS(I,J,bi,bj) - uIceRHS(I,J,bi,bj)
                0145            vIceRes(I,J,bi,bj) = vIceLHS(I,J,bi,bj) - vIceRHS(I,J,bi,bj)
                0146           ENDDO
                0147          ENDDO
                0148         ENDDO
                0149        ENDDO
                0150
45315406aa Mart*0151 #endif /* SEAICE_CGRID and SEAICE_ALLOW_JFNK */
2fd913c523 Mart*0152
                0153       RETURN
                0154       END