pkg/seaice/seaice_ocean_stress.F

53092bcb42 Mart*0001 #include "SEAICE_OPTIONS.h"
                0002
17bd6b9422 jm-c 0003 CBOP
                0004 C     !ROUTINE: SEAICE_OCEAN_STRESS
                0005 C     !INTERFACE:
b4949dd6db Jean*0006       SUBROUTINE SEAICE_OCEAN_STRESS(
17bd6b9422 jm-c 0007      I                  windTauX, windTauY,
                0008      I                  myTime, myIter, myThid )
                0009 C     !DESCRIPTION: \bv
03c669d1ab Jean*0010 C     *==========================================================*
53092bcb42 Mart*0011 C     | SUBROUTINE SEAICE_OCEAN_STRESS                           |
                0012 C     | o Calculate ocean surface stresses                       |
                0013 C     |   - C-grid version                                       |
03c669d1ab Jean*0014 C     *==========================================================*
17bd6b9422 jm-c 0015 C     \ev
                0016
                0017 C     !USES:
53092bcb42 Mart*0018       IMPLICIT NONE
                0019
                0020 C     === Global variables ===
                0021 #include "SIZE.h"
                0022 #include "EEPARAMS.h"
                0023 #include "PARAMS.h"
bcf4255045 Dimi*0024 #include "DYNVARS.h"
3ed8349f04 Mart*0025 #include "GRID.h"
53092bcb42 Mart*0026 #include "FFIELDS.h"
03c669d1ab Jean*0027 #include "SEAICE_SIZE.h"
53092bcb42 Mart*0028 #include "SEAICE_PARAMS.h"
03c669d1ab Jean*0029 #include "SEAICE.h"
53092bcb42 Mart*0030
17bd6b9422 jm-c 0031 C     !INPUT/OUTPUT PARAMETERS:
                0032 C     windTauX :: X-direction wind stress over seaice at U point
                0033 C     windTauY :: Y-direction wind stress over seaice at V point
                0034 C     myTime   :: Current time in simulation
                0035 C     myIter   :: Current iteration number
                0036 C     myThid   :: my Thread Id number
                0037       _RL     windTauX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0038       _RL     windTauY(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
53092bcb42 Mart*0039       _RL     myTime
                0040       INTEGER myIter
                0041       INTEGER myThid
                0042
17bd6b9422 jm-c 0043 C     !LOCAL VARIABLES:
                0044 C     i,j,bi,bj :: Loop counters
                0045 C     kSrf      :: vertical index of surface layer
53092bcb42 Mart*0046       INTEGER i, j, bi, bj
486cba1d02 Mart*0047       INTEGER kSrf
35fda33b05 Jean*0048       _RL  COSWAT
                0049       _RS  SINWAT
8b339058e0 Mart*0050       _RL  fuIceLoc, fvIceLoc
4d5b524a4b Mart*0051       _RL  areaW, areaS
17bd6b9422 jm-c 0052 CEOP
53092bcb42 Mart*0053
486cba1d02 Mart*0054 C     surrface level
0320e25227 Mart*0055       IF ( usingPCoords ) THEN
                0056        kSrf = Nr
                0057       ELSE
                0058        kSrf = 1
                0059       ENDIF
486cba1d02 Mart*0060 C     introduce turning angle (default is zero)
53092bcb42 Mart*0061       SINWAT=SIN(SEAICE_waterTurnAngle*deg2rad)
                0062       COSWAT=COS(SEAICE_waterTurnAngle*deg2rad)
                0063
3ed8349f04 Mart*0064       IF ( useHB87StressCoupling ) THEN
                0065 C
b4949dd6db Jean*0066 C     use an intergral over ice and ocean surface layer to define
3ed8349f04 Mart*0067 C     surface stresses on ocean following Hibler and Bryan (1987, JPO)
b4949dd6db Jean*0068 C
3ed8349f04 Mart*0069        DO bj=myByLo(myThid),myByHi(myThid)
                0070         DO bi=myBxLo(myThid),myBxHi(myThid)
8b339058e0 Mart*0071          DO J=1,sNy
                0072           DO I=1,sNx
b4949dd6db Jean*0073 C     average wind stress over ice and ocean and apply averaged wind
e81deec45f Mart*0074 C     stress and internal ice stresses to surface layer of ocean
772590b63c Mart*0075            areaW = 0.5 * (AREA(I,J,bi,bj) + AREA(I-1,J,bi,bj))
8b339058e0 Mart*0076      &          * SEAICEstressFactor
                0077            fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)
17bd6b9422 jm-c 0078      &          + areaW*windTauX(I,J,bi,bj)
45315406aa Mart*0079 #ifdef SEAICE_CGRID
8b339058e0 Mart*0080      &          + stressDivergenceX(I,J,bi,bj) * SEAICEstressFactor
45315406aa Mart*0081 #endif
f38acc47ff Mart*0082           ENDDO
                0083          ENDDO
225a351aca Mart*0084 C     This loop separation makes the adjoint code vectorize
f38acc47ff Mart*0085          DO J=1,sNy
                0086           DO I=1,sNx
                0087            areaS = 0.5 * (AREA(I,J,bi,bj) + AREA(I,J-1,bi,bj))
                0088      &          * SEAICEstressFactor
8b339058e0 Mart*0089            fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)
17bd6b9422 jm-c 0090      &          + areaS*windTauY(I,J,bi,bj)
45315406aa Mart*0091 #ifdef SEAICE_CGRID
8b339058e0 Mart*0092      &          + stressDivergenceY(I,J,bi,bj) * SEAICEstressFactor
45315406aa Mart*0093 #endif
e81deec45f Mart*0094           ENDDO
8b339058e0 Mart*0095          ENDDO
3ed8349f04 Mart*0096         ENDDO
                0097        ENDDO
8a6cc4f42a Jean*0098
3ed8349f04 Mart*0099       ELSE
8a6cc4f42a Jean*0100 C     else: useHB87StressCoupling=F
3ed8349f04 Mart*0101
b4949dd6db Jean*0102 C--   Compute ice-affected wind stress (interpolate to U/V-points)
                0103 C     by averaging wind stress and ice-ocean stress according to
3ed8349f04 Mart*0104 C     ice cover
53092bcb42 Mart*0105       DO bj=myByLo(myThid),myByHi(myThid)
                0106        DO bi=myBxLo(myThid),myBxHi(myThid)
                0107         DO j=1,sNy
                0108          DO i=1,sNx
4b491e6964 Mart*0109           fuIceLoc=HALF*( DWATN(I,J,bi,bj)+DWATN(I-1,J,bi,bj) )*
b4949dd6db Jean*0110      &         COSWAT *
486cba1d02 Mart*0111      &         ( uIce(I,J,bi,bj)-uVel(I,J,kSrf,bi,bj) )
b4949dd6db Jean*0112      &         - SIGN(SINWAT, _fCori(I,J,bi,bj)) * 0.5 _d 0 *
ff10a08150 Mart*0113      &         ( DWATN(I  ,J,bi,bj) *
486cba1d02 Mart*0114      &         0.5 _d 0*(vIce(I  ,J  ,bi,bj)-vVel(I  ,J  ,kSrf,bi,bj)
                0115      &                  +vIce(I  ,J+1,bi,bj)-vVel(I  ,J+1,kSrf,bi,bj))
ff10a08150 Mart*0116      &         + DWATN(I-1,J,bi,bj) *
486cba1d02 Mart*0117      &         0.5 _d 0*(vIce(I-1,J  ,bi,bj)-vVel(I-1,J  ,kSrf,bi,bj)
                0118      &                  +vIce(I-1,J+1,bi,bj)-vVel(I-1,J+1,kSrf,bi,bj))
53092bcb42 Mart*0119      &         )
4b491e6964 Mart*0120           fvIceLoc=HALF*( DWATN(I,J,bi,bj)+DWATN(I,J-1,bi,bj) )*
ff10a08150 Mart*0121      &         COSWAT *
486cba1d02 Mart*0122      &         ( vIce(I,J,bi,bj)-vVel(I,J,kSrf,bi,bj) )
ff10a08150 Mart*0123      &         + SIGN(SINWAT,  _fCori(I,J,bi,bj)) * 0.5 _d 0 *
                0124      &         ( DWATN(I,J  ,bi,bj) *
486cba1d02 Mart*0125      &         0.5 _d 0*(uIce(I  ,J  ,bi,bj)-uVel(I  ,J  ,kSrf,bi,bj)
                0126      &                  +uIce(I+1,J  ,bi,bj)-uVel(I+1,J  ,kSrf,bi,bj))
ff10a08150 Mart*0127      &         + DWATN(I,J-1,bi,bj) *
486cba1d02 Mart*0128      &         0.5 _d 0*(uIce(I  ,J-1,bi,bj)-uVel(I  ,J-1,kSrf,bi,bj)
                0129      &                  +uIce(I+1,J-1,bi,bj)-uVel(I+1,J-1,kSrf,bi,bj))
53092bcb42 Mart*0130      &         )
772590b63c Mart*0131           areaW = 0.5 _d 0 * (AREA(I,J,bi,bj) + AREA(I-1,J,bi,bj))
936a01b1f8 Mart*0132      &         * SEAICEstressFactor
772590b63c Mart*0133           areaS = 0.5 _d 0 * (AREA(I,J,bi,bj) + AREA(I,J-1,bi,bj))
936a01b1f8 Mart*0134      &         * SEAICEstressFactor
0dfb864288 Mart*0135           fu(I,J,bi,bj)=(ONE-areaW)*fu(I,J,bi,bj)+areaW*fuIceLoc
                0136           fv(I,J,bi,bj)=(ONE-areaS)*fv(I,J,bi,bj)+areaS*fvIceLoc
53092bcb42 Mart*0137          ENDDO
                0138         ENDDO
                0139        ENDDO
                0140       ENDDO
3ed8349f04 Mart*0141       ENDIF
53092bcb42 Mart*0142       CALL EXCH_UV_XY_RS(fu, fv, .TRUE., myThid)
93ea7880e3 Mart*0143
53092bcb42 Mart*0144       RETURN
                0145       END