pkg/seaice/seaice_calc_lhs.F

2fd913c523 Mart*0001 #include "SEAICE_OPTIONS.h"
                0002
                0003 CBOP
                0004 C     !ROUTINE: SEAICE_CALC_LHS
                0005 C     !INTERFACE:
                0006       SUBROUTINE SEAICE_CALC_LHS(
                0007      I     uIceLoc, vIceLoc,
                0008      O     uIceLHS, vIceLHS,
                0009      I     newtonIter, myTime, myIter, myThid )
                0010
                0011 C     !DESCRIPTION: \bv
                0012 C     *==========================================================*
                0013 C     | SUBROUTINE SEAICE_CALC_LHS
                0014 C     | o Left-hand side of momentum equations, i.e. all terms
5acccad966 Jean*0015 C     |   that depend on the ice velocities of the current
2fd913c523 Mart*0016 C     |   iterate of the Newton-Krylov iteration
                0017 C     |
5acccad966 Jean*0018 C     | o The scheme is backward Euler in time, i.e. the
2fd913c523 Mart*0019 C     |   rhs-vector contains only terms that are independent
5acccad966 Jean*0020 C     |   of u/vIce, except for the time derivative part
2fd913c523 Mart*0021 C     |   mass*(u/vIce-u/vIceNm1)/deltaT
                0022 C     | o Left-hand side contributions
                0023 C     |   + mass*(u/vIce)/deltaT
                0024 C     |   + Cdrag*(uIce*cosWat - vIce*sinWat)
                0025 C     |          /(vIce*cosWat + uIce*sinWat)
                0026 C     |   - mass*f*vIce/+mass*f*uIce
                0027 C     |   - dsigma/dx / -dsigma/dy, eta and zeta are
                0028 C     |                   computed only once per Newton iterate
                0029 C     *==========================================================*
                0030 C     | written by Martin Losch, Oct 2012
                0031 C     *==========================================================*
                0032 C     \ev
                0033
                0034 C     !USES:
                0035       IMPLICIT NONE
                0036
                0037 C     === Global variables ===
                0038 #include "SIZE.h"
                0039 #include "EEPARAMS.h"
                0040 #include "PARAMS.h"
                0041 #include "GRID.h"
                0042 #include "SEAICE_SIZE.h"
                0043 #include "SEAICE_PARAMS.h"
                0044 #include "SEAICE.h"
                0045
                0046 C     !INPUT/OUTPUT PARAMETERS:
                0047 C     === Routine arguments ===
                0048 C     myTime :: Simulation time
                0049 C     myIter :: Simulation timestep number
                0050 C     myThid :: my Thread Id. number
                0051 C     newtonIter :: current iterate of Newton iteration
                0052       _RL     myTime
                0053       INTEGER myIter
                0054       INTEGER myThid
                0055       INTEGER newtonIter
                0056 C     u/vIceLoc :: local copies of the current ice velocity
                0057       _RL uIceLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0058       _RL vIceLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0059 C     u/vIceLHS :: LHS of momentum equations
                0060       _RL uIceLHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0061       _RL vIceLHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0062
45315406aa Mart*0063 #if ( defined SEAICE_CGRID \
                0064       && ( defined SEAICE_ALLOW_JFNK || defined SEAICE_ALLOW_KRYLOV ) )
2fd913c523 Mart*0065 C     i,j,bi,bj,k :: loop indices
                0066       INTEGER i,j,bi,bj
                0067       INTEGER k
5acccad966 Jean*0068       _RS     SINWAT
cbf501ab81 Jean*0069       _RL     COSWAT, recip_deltaT
e501eee760 Mart*0070 C     backward difference extrapolation factor
                0071       _RL bdfAlpha
df1dac8b7b Mart*0072 C     symmetric drag coefficient
                0073       _RL dragSym(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70e078b38a Mart*0074 C     fractional area at velocity points
                0075       _RL areaW(1:sNx,1:sNy)
                0076       _RL areaS(1:sNx,1:sNy)
8ce59fd29e Mart*0077 #ifdef SEAICE_ALLOW_MOM_ADVECTION
                0078 C     tendency due to advection of momentum
                0079       _RL gUmom(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0080       _RL gVmom(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0081 #endif /*  SEAICE_ALLOW_MOM_ADVECTION */
2fd913c523 Mart*0082 CEOP
                0083
                0084       k=1
                0085       recip_deltaT = 1. _d 0 / SEAICE_deltaTdyn
                0086 C--   introduce turning angles
                0087       SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)
                0088       COSWAT=COS(SEAICE_waterTurnAngle*deg2rad)
e501eee760 Mart*0089 C     backward difference extrapolation factor
                0090       bdfAlpha = 1. _d 0
                0091       IF ( SEAICEuseBDF2 ) THEN
                0092        IF ( myIter.EQ.nIter0 .AND. SEAICEmomStartBDF.EQ.0 ) THEN
                0093         bdfAlpha = 1. _d 0
6cbc659de0 Mart*0094        ELSE
e501eee760 Mart*0095         bdfAlpha = 1.5 _d 0
6cbc659de0 Mart*0096        ENDIF
                0097       ENDIF
2fd913c523 Mart*0098
70e078b38a Mart*0099 C     initialise fractional areas at velocity points
                0100       DO J=1,sNy
                0101        DO I=1,sNx
                0102         areaW(I,J) = 1. _d 0
                0103         areaS(I,J) = 1. _d 0
                0104        ENDDO
                0105       ENDDO
                0106
2fd913c523 Mart*0107       DO bj=myByLo(myThid),myByHi(myThid)
                0108        DO bi=myBxLo(myThid),myBxHi(myThid)
df1dac8b7b Mart*0109 C     symmetric drag coefficient may include bottomdrag for grounded ice
                0110         DO j=1-OLy,sNy+OLy
                0111          DO i=1-OLx,sNx+OLx
                0112           dragSym(I,J) = DWATN(I,J,bi,bj)*COSWAT
                0113 #ifdef SEAICE_ALLOW_BOTTOMDRAG
                0114      &         +CbotC(I,J,bi,bj)
                0115 #endif /* SEAICE_ALLOW_BOTTOMDRAG */
                0116          ENDDO
                0117         ENDDO
cbf501ab81 Jean*0118 C
2fd913c523 Mart*0119 C     compute components of stress tensor from current velocity field
925df4f4b9 Mart*0120 C     and compute divergence of stress tensor
2fd913c523 Mart*0121 C
cbf501ab81 Jean*0122         CALL SEAICE_CALC_STRESSDIV(
925df4f4b9 Mart*0123      I       e11, e22, e12, press, zeta, eta, etaZ,
                0124      O       stressDivergenceX, stressDivergenceY,
                0125      I       bi, bj, myTime, myIter, myThid )
2fd913c523 Mart*0126 C     compute lhs side of momentum equations
70e078b38a Mart*0127         IF ( SEAICEscaleSurfStress ) THEN
                0128          DO J=1,sNy
                0129           DO I=1,sNx
                0130            areaW(I,J) = 0.5 _d 0*(AREA(I,J,bi,bj)+AREA(I-1,J,bi,bj))
d088941345 Mart*0131            areaS(I,J) = 0.5 _d 0*(AREA(I,J,bi,bj)+AREA(I,J-1,bi,bj))
70e078b38a Mart*0132           ENDDO
                0133          ENDDO
                0134         ENDIF
2fd913c523 Mart*0135         DO J=1,sNy
                0136          DO I=1,sNx
5acccad966 Jean*0137 C     mass*(uIce)/deltaT - dsigma/dx
cbf501ab81 Jean*0138           uIceLHS(I,J,bi,bj) =
e501eee760 Mart*0139      &         bdfAlpha*seaiceMassU(I,J,bi,bj)*recip_deltaT
2fd913c523 Mart*0140      &         *uIceLoc(I,J,bi,bj) - stressDivergenceX(I,J,bi,bj)
                0141 C     mass*(vIce)/deltaT - dsigma/dy
cbf501ab81 Jean*0142           vIceLHS(I,J,bi,bj) =
e501eee760 Mart*0143      &         bdfAlpha*seaiceMassV(I,J,bi,bj)*recip_deltaT
2fd913c523 Mart*0144      &         *vIceLoc(I,J,bi,bj) - stressDivergenceY(I,J,bi,bj)
                0145 C     coriols terms: - mass*f*vIce
                0146           uIceLHS(I,J,bi,bj) = uIceLHS(I,J,bi,bj) - 0.5 _d 0*(
                0147      &         seaiceMassC(I  ,J,bi,bj) * _fCori(I  ,J,bi,bj)
                0148      &       * 0.5 _d 0*( vIceLoc(I  ,J,bi,bj)+vIceLoc(I  ,J+1,bi,bj) )
                0149      &       + seaiceMassC(I-1,J,bi,bj) * _fCori(I-1,J,bi,bj)
                0150      &       * 0.5 _d 0*( vIceLoc(I-1,J,bi,bj)+vIceLoc(I-1,J+1,bi,bj) )
                0151      &           )
                0152 C                    + mass*f*uIce
                0153           vIceLHS(I,J,bi,bj) = vIceLHS(I,J,bi,bj) + 0.5 _d 0*(
                0154      &         seaiceMassC(I,J  ,bi,bj) * _fCori(I,J  ,bi,bj)
                0155      &       * 0.5 _d 0*( uIceLoc(I,J  ,bi,bj)+uIceLoc(I+1,  J,bi,bj) )
                0156      &       + seaiceMassC(I,J-1,bi,bj) * _fCori(I,J-1,bi,bj)
                0157      &       * 0.5 _d 0*( uIceLoc(I,J-1,bi,bj)+uIceLoc(I+1,J-1,bi,bj) )
                0158      &           )
df1dac8b7b Mart*0159 C     ocean-ice and bottom drag terms: + (Cdrag*cosWat+Cb)*uIce - vIce*sinWat)
bce2e149b6 Mart*0160           uIceLHS(I,J,bi,bj) = uIceLHS(I,J,bi,bj) + (
df1dac8b7b Mart*0161      &         0.5 _d 0 * ( dragSym(I,J)+dragSym(I-1,J) )
                0162      &         * uIceLoc(I,J,bi,bj)
2fd913c523 Mart*0163      &         - SIGN(SINWAT, _fCori(I,J,bi,bj))* 0.5 _d 0 *
                0164      &         (  DWATN(I  ,J,bi,bj) * 0.5 _d 0 *
                0165      &         (vIceLoc(I  ,J,bi,bj)+vIceLoc(I  ,J+1,bi,bj))
                0166      &         +  DWATN(I-1,J,bi,bj) * 0.5 _d 0 *
                0167      &         (vIceLoc(I-1,J,bi,bj)+vIceLoc(I-1,J+1,bi,bj))
bce2e149b6 Mart*0168      &         ) ) * areaW(I,J)
df1dac8b7b Mart*0169 C                                      + (Cdrag*cosWat+Cb)*uIce + uIce*sinWat)
bce2e149b6 Mart*0170           vIceLHS(I,J,bi,bj) = vIceLHS(I,J,bi,bj) + (
df1dac8b7b Mart*0171      &         0.5 _d 0 * ( dragSym(I,J)+dragSym(I,J-1) )
                0172      &         * vIceLoc(I,J,bi,bj)
2fd913c523 Mart*0173      &         + SIGN(SINWAT, _fCori(I,J,bi,bj)) * 0.5 _d 0 *
                0174      &         (  DWATN(I,J  ,bi,bj) * 0.5 _d 0 *
                0175      &         (uIceLoc(I,J  ,bi,bj)+uIceLoc(I+1,J  ,bi,bj))
                0176      &         +  DWATN(I,J-1,bi,bj) * 0.5 _d 0 *
                0177      &         (uIceLoc(I,J-1,bi,bj)+uIceLoc(I+1,J-1,bi,bj))
bce2e149b6 Mart*0178      &         ) ) * areaS(I,J)
ad40ebac66 Mart*0179 C     apply masks for interior (important when we have open boundaries)
                0180           uIceLHS(I,J,bi,bj) = uIceLHS(I,J,bi,bj)*maskinW(I,J,bi,bj)
                0181           vIceLHS(I,J,bi,bj) = vIceLHS(I,J,bi,bj)*maskinS(I,J,bi,bj)
2fd913c523 Mart*0182          ENDDO
                0183         ENDDO
8ce59fd29e Mart*0184 #ifdef SEAICE_ALLOW_MOM_ADVECTION
cbf501ab81 Jean*0185         IF ( SEAICEmomAdvection ) THEN
                0186          DO J=1-OLy,sNy+OLy
                0187           DO I=1-OLx,sNx+OLx
8ce59fd29e Mart*0188            gUmom(I,J) = 0. _d 0
                0189            gVmom(I,J) = 0. _d 0
                0190           ENDDO
                0191          ENDDO
                0192          CALL SEAICE_MOM_ADVECTION(
                0193      I        bi,bj,1,sNx,1,sNy,
                0194      I        uIceLoc, vIceLoc,
                0195      O        gUmom, gVmom,
                0196      I        myTime, myIter, myThid )
b3193ad894 Mart*0197 C     Beware of sign! gU/Vmom is computed for the rhs of the equation;
                0198 C     therefore, we need to substract gU/Vmom from the left hand side
8ce59fd29e Mart*0199          DO J=1,sNy
                0200           DO I=1,sNx
b3193ad894 Mart*0201            uIceLHS(I,J,bi,bj) = uIceLHS(I,J,bi,bj) - gUmom(I,J)
                0202            vIceLHS(I,J,bi,bj) = vIceLHS(I,J,bi,bj) - gVmom(I,J)
8ce59fd29e Mart*0203           ENDDO
                0204          ENDDO
                0205         ENDIF
                0206 #endif /* SEAICE_ALLOW_MOM_ADVECTION */
2fd913c523 Mart*0207        ENDDO
                0208       ENDDO
                0209
45315406aa Mart*0210 #endif /* SEAICE_CGRID, SEAICE_ALLOW_JFNK and KRYLOV */
2fd913c523 Mart*0211
                0212       RETURN
                0213       END