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
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 */
5bb179ddc2 Mart*0094 #ifdef SEAICE_ALLOW_SIDEDRAG
                0095         IF ( SEAICEsideDrag .NE. 0. _d 0 ) CALL SEAICE_SIDEDRAG_STRESS(
                0096      I       uIceLoc, vIceLoc, coastRoughU, coastRoughV, AREA,
                0097      O       sideDragU, sideDragV,
                0098      I       krylovIter, myTime, myIter, myThid )
                0099 #endif /* SEAICE_ALLOW_SIDEDRAG */
2fd913c523 Mart*0100       CALL SEAICE_CALC_STRAINRATES(
                0101      I     uIceLoc, vIceLoc,
                0102      O     e11, e22, e12,
                0103      I     krylovIter, myTime, myIter, myThid )
                0104       CALL SEAICE_CALC_VISCOSITIES(
8e32c48b8f Mart*0105      I     e11, e22, e12, SEAICE_zMin, SEAICE_zMax, HEFFM, press0,
                0106      I     tensileStrFac,
0fd6b2de1a Mart*0107      O     eta, etaZ, zeta, zetaZ, press, deltaC,
2fd913c523 Mart*0108      I     krylovIter, myTime, myIter, myThid )
                0109
ec0d7df165 Mart*0110 C     The scheme is backward Euler in time, i.e. the rhs-vector contains
                0111 C     only terms that are independent of u/vIce, except for the time
2fd913c523 Mart*0112 C     derivative part mass*(u/vIce-u/vIceNm1)/deltaT
                0113
                0114 C     compute new right hand side (depends to DWATN=Cdrag)
                0115 C     sea-surface tilt and wind stress: FORCEX0, FORCEY0
                0116 C     + mass*(u/vIceNm1)/deltaT
                0117 C     + Cdrag*(uVel*cosWat - vVel*sinWat)/(vVel*cosWat + uVel*sinWat)
                0118       CALL SEAICE_CALC_RHS(
                0119      O      uIceRHS, vIceRHS,
                0120      I      newtonIter, krylovIter, myTime, myIter, myThid )
                0121
                0122 C     Left-hand side contributions:
                0123 C     + mass*(u/vIce)/deltaT
                0124 C     + Cdrag*(uIce*cosWat - vIce*sinWat)/(vIce*cosWat + uIce*sinWat)
df1dac8b7b Mart*0125 C     + CdragBot*uIce/vIce
2fd913c523 Mart*0126 C     - mass*f*vIce/+mass*f*uIce
                0127 C     - dsigma/dx / -dsigma/dy, eta and zeta are only computed once per
                0128 C     Newton iterate
                0129        CALL SEAICE_CALC_LHS(
                0130      I      uIceLoc, vIceLoc,
                0131      O      uIceLHS, vIceLHS,
                0132      I      newtonIter, myTime, myIter, myThid )
                0133
                0134 C     Right-hand side contributions only need to be computed once per
                0135 C     time step, therefore we will put them into a separate routine
                0136 C     and call them elsewhere to save floating point operations
                0137
                0138 C     Calculate the residual
                0139        DO bj=myByLo(myThid),myByHi(myThid)
                0140         DO bi=myBxLo(myThid),myBxHi(myThid)
                0141          DO J=1,sNy
                0142           DO I=1,sNx
                0143            uIceRes(I,J,bi,bj) = uIceLHS(I,J,bi,bj) - uIceRHS(I,J,bi,bj)
                0144            vIceRes(I,J,bi,bj) = vIceLHS(I,J,bi,bj) - vIceRHS(I,J,bi,bj)
                0145           ENDDO
                0146          ENDDO
                0147         ENDDO
                0148        ENDDO
                0149
45315406aa Mart*0150 #endif /* SEAICE_CGRID and SEAICE_ALLOW_JFNK */
2fd913c523 Mart*0151
                0152       RETURN
                0153       END