pkg/seaice/seaice_calc_rhs.F

2fd913c523 Mart*0001 #include "SEAICE_OPTIONS.h"
                0002
                0003 CBOP
                0004 C     !ROUTINE: SEAICE_CALC_RHS
                0005 C     !INTERFACE:
                0006       SUBROUTINE SEAICE_CALC_RHS(
                0007      O     uIceRHS, vIceRHS,
                0008      I     newtonIter, krylovIter, myTime, myIter, myThid )
                0009
                0010 C     !DESCRIPTION: \bv
                0011 C     *==========================================================*
                0012 C     | SUBROUTINE SEAICE_CALC_RHS
                0013 C     | o Right-hand side of momentum equations, i.e. all terms
5acccad966 Jean*0014 C     |   that do not depend on the ice velocities of the current
2fd913c523 Mart*0015 C     |   iterate of the Newton-Krylov iteration
                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 "PARAMS.h"
                0028 #include "DYNVARS.h"
                0029 #include "GRID.h"
                0030 #include "SEAICE_SIZE.h"
                0031 #include "SEAICE_PARAMS.h"
                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/vIceRHS :: RHS of momentum equations
                0047       _RL uIceRHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0048       _RL vIceRHS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0049
45315406aa Mart*0050 #if ( defined SEAICE_CGRID \
                0051       && ( defined SEAICE_ALLOW_JFNK || defined SEAICE_ALLOW_KRYLOV ) )
2fd913c523 Mart*0052 C     i,j,bi,bj,kSrf :: loop indices
                0053       INTEGER i,j,bi,bj
                0054       INTEGER kSrf
5acccad966 Jean*0055       _RS     SINWAT
                0056       _RL     COSWAT
70e078b38a Mart*0057 C     fractional area at velocity points
                0058       _RL areaW(1:sNx,1:sNy)
                0059       _RL areaS(1:sNx,1:sNy)
2fd913c523 Mart*0060 CEOP
                0061
0320e25227 Mart*0062       IF ( usingPCoords ) THEN
                0063        kSrf = Nr
                0064       ELSE
                0065        kSrf = 1
                0066       ENDIF
2fd913c523 Mart*0067 C--   introduce turning angles
                0068       SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)
                0069       COSWAT=COS(SEAICE_waterTurnAngle*deg2rad)
                0070
70e078b38a Mart*0071 C     initialise fractional areas at velocity points
                0072       DO J=1,sNy
                0073        DO I=1,sNx
                0074         areaW(I,J) = 1. _d 0
                0075         areaS(I,J) = 1. _d 0
                0076        ENDDO
                0077       ENDDO
                0078
2fd913c523 Mart*0079       DO bj=myByLo(myThid),myByHi(myThid)
                0080        DO bi=myBxLo(myThid),myBxHi(myThid)
70e078b38a Mart*0081         IF ( SEAICEscaleSurfStress ) THEN
                0082          DO J=1,sNy
                0083           DO I=1,sNx
                0084            areaW(I,J) = 0.5 _d 0*(AREA(I,J,bi,bj)+AREA(I-1,J,bi,bj))
d088941345 Mart*0085            areaS(I,J) = 0.5 _d 0*(AREA(I,J,bi,bj)+AREA(I,J-1,bi,bj))
70e078b38a Mart*0086           ENDDO
                0087          ENDDO
                0088         ENDIF
2fd913c523 Mart*0089         DO J=1,sNy
                0090          DO I=1,sNx
                0091 C     ice-velocity independent contribution to drag terms
                0092 C     - Cdrag*(uVel*cosWat - vVel*sinWat)/(vVel*cosWat + uVel*sinWat)
                0093 C     ( remember to average to correct velocity points )
                0094           uIceRHS(I,J,bi,bj) = FORCEX(I,J,bi,bj) + (
                0095      &         0.5 _d 0 * ( DWATN(I,J,bi,bj)+DWATN(I-1,J,bi,bj) ) *
                0096      &         COSWAT * uVel(I,J,kSrf,bi,bj)
                0097      &         - SIGN(SINWAT, _fCori(I,J,bi,bj))* 0.5 _d 0 *
                0098      &         ( DWATN(I  ,J,bi,bj) * 0.5 _d 0 *
                0099      &          ( vVel(I  ,J,kSrf,bi,bj)+vVel(I  ,J+1,kSrf,bi,bj) )
                0100      &         + DWATN(I-1,J,bi,bj) * 0.5 _d 0 *
                0101      &          ( vVel(I-1,J,kSrf,bi,bj)+vVel(I-1,J+1,kSrf,bi,bj) )
70e078b38a Mart*0102      &         ) )*areaW(I,J)
2fd913c523 Mart*0103           vIceRHS(I,J,bi,bj) = FORCEY(I,J,bi,bj) + (
                0104      &         0.5 _d 0 * ( DWATN(I,J,bi,bj)+DWATN(I,J-1,bi,bj) ) *
                0105      &         COSWAT * vVel(I,J,kSrf,bi,bj)
                0106      &         + SIGN(SINWAT, _fCori(I,J,bi,bj)) * 0.5 _d 0 *
                0107      &         ( DWATN(I,J  ,bi,bj) * 0.5 _d 0 *
                0108      &          ( uVel(I,J  ,kSrf,bi,bj)+uVel(I+1,J  ,kSrf,bi,bj))
                0109      &         + DWATN(I,J-1,bi,bj) * 0.5 _d 0 *
                0110      &          ( uVel(I,J-1,kSrf,bi,bj)+uVel(I+1,J-1,kSrf,bi,bj))
70e078b38a Mart*0111      &         ) )*areaS(I,J)
ad40ebac66 Mart*0112 C     apply masks for interior (important when we have open boundaries)
                0113           uIceRHS(I,J,bi,bj) = uIceRHS(I,J,bi,bj)*maskinW(I,J,bi,bj)
                0114           vIceRHS(I,J,bi,bj) = vIceRHS(I,J,bi,bj)*maskinS(I,J,bi,bj)
2fd913c523 Mart*0115          ENDDO
                0116         ENDDO
                0117        ENDDO
                0118       ENDDO
                0119
45315406aa Mart*0120 #endif /* SEAICE_CGRID, SEAICE_ALLOW_JFNK and KRYLOV */
2fd913c523 Mart*0121
                0122       RETURN
                0123       END