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