Back to home page

darwin3
 
 

    

File indexing completed on 2024-12-17 18:38:22 UTC

view on githubraw file Latest commit 5bb179dd on 2024-10-17 18:00:27 UTC
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)
5bb179ddc2 Mart*0179 #ifdef SEAICE_ALLOW_SIDEDRAG
                0180 C     lateral drag terms: + sideDragU*uIce
                0181           uIceLHS(I,J,bi,bj) = uIceLHS(I,J,bi,bj)
                0182      &         + sideDragU(I,J,bi,bj)*uIceLoc(I,J,bi,bj)
                0183 C                         + sideDragV*vIce
                0184           vIceLHS(I,J,bi,bj) = vIceLHS(I,J,bi,bj)
                0185      &         + sideDragV(I,J,bi,bj)*vIceLoc(I,J,bi,bj)
                0186 #endif /* SEAICE_ALLOW_SIDEDRAG */
ad40ebac66 Mart*0187 C     apply masks for interior (important when we have open boundaries)
                0188           uIceLHS(I,J,bi,bj) = uIceLHS(I,J,bi,bj)*maskinW(I,J,bi,bj)
                0189           vIceLHS(I,J,bi,bj) = vIceLHS(I,J,bi,bj)*maskinS(I,J,bi,bj)
2fd913c523 Mart*0190          ENDDO
                0191         ENDDO
8ce59fd29e Mart*0192 #ifdef SEAICE_ALLOW_MOM_ADVECTION
cbf501ab81 Jean*0193         IF ( SEAICEmomAdvection ) THEN
                0194          DO J=1-OLy,sNy+OLy
                0195           DO I=1-OLx,sNx+OLx
8ce59fd29e Mart*0196            gUmom(I,J) = 0. _d 0
                0197            gVmom(I,J) = 0. _d 0
                0198           ENDDO
                0199          ENDDO
                0200          CALL SEAICE_MOM_ADVECTION(
                0201      I        bi,bj,1,sNx,1,sNy,
                0202      I        uIceLoc, vIceLoc,
                0203      O        gUmom, gVmom,
                0204      I        myTime, myIter, myThid )
b3193ad894 Mart*0205 C     Beware of sign! gU/Vmom is computed for the rhs of the equation;
                0206 C     therefore, we need to substract gU/Vmom from the left hand side
8ce59fd29e Mart*0207          DO J=1,sNy
                0208           DO I=1,sNx
b3193ad894 Mart*0209            uIceLHS(I,J,bi,bj) = uIceLHS(I,J,bi,bj) - gUmom(I,J)
                0210            vIceLHS(I,J,bi,bj) = vIceLHS(I,J,bi,bj) - gVmom(I,J)
8ce59fd29e Mart*0211           ENDDO
                0212          ENDDO
                0213         ENDIF
                0214 #endif /* SEAICE_ALLOW_MOM_ADVECTION */
2fd913c523 Mart*0215        ENDDO
                0216       ENDDO
                0217 
45315406aa Mart*0218 #endif /* SEAICE_CGRID, SEAICE_ALLOW_JFNK and KRYLOV */
2fd913c523 Mart*0219 
                0220       RETURN
                0221       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/darwinproject/darwin3 by the 2.3.7-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!