pkg/seaice/seaice_get_dynforcing.F

0dfb864288 Mart*0001 #include "SEAICE_OPTIONS.h"
8bf64a0174 Jean*0002 #ifdef ALLOW_EXF
                0003 # include "EXF_OPTIONS.h"
                0004 #endif
0dfb864288 Mart*0005
23cba332b3 Jean*0006 CBOP
8bf64a0174 Jean*0007 C     !ROUTINE: SEAICE_GET_DYNFORCING
23cba332b3 Jean*0008 C     !INTERFACE:
35fda33b05 Jean*0009       SUBROUTINE SEAICE_GET_DYNFORCING(
ec0d7df165 Mart*0010      I     uIce, vIce, icFrac, SIMaskU, SIMaskV,
8a6cc4f42a Jean*0011      O     taux, tauy,
0dfb864288 Mart*0012      I     myTime, myIter, myThid )
23cba332b3 Jean*0013 C     !DESCRIPTION: \bv
                0014 C     *==========================================================*
8a6cc4f42a Jean*0015 C     | SUBROUTINE SEAICE_GET_DYNFORCING
                0016 C     |   compute surface stress from atmopheric forcing fields
23cba332b3 Jean*0017 C     *==========================================================*
8a6cc4f42a Jean*0018 C     | started by Martin Losch, April 2007
23cba332b3 Jean*0019 C     *==========================================================*
                0020 C     \ev
0dfb864288 Mart*0021
23cba332b3 Jean*0022 C     !USES:
                0023       IMPLICIT NONE
0dfb864288 Mart*0024 C     === Global variables ===
                0025 #include "SIZE.h"
                0026 #include "EEPARAMS.h"
                0027 #include "PARAMS.h"
                0028 #include "GRID.h"
                0029 #include "FFIELDS.h"
7303eab4f2 Patr*0030 #include "SEAICE_SIZE.h"
0dfb864288 Mart*0031 #include "SEAICE_PARAMS.h"
855e83a213 Jean*0032 #ifdef HACK_FOR_GMAO_CPL
                0033 # include "SEAICE_LAYERS.h"
                0034 #endif
ae1fb66b64 Dimi*0035 #ifdef ALLOW_EXF
8bf64a0174 Jean*0036 # include "DYNVARS.h"
ae1fb66b64 Dimi*0037 # include "EXF_FIELDS.h"
58a71b9205 Dimi*0038 # include "EXF_PARAM.h"
ae1fb66b64 Dimi*0039 #endif
0dfb864288 Mart*0040
9da8ce8edb Jean*0041 #ifdef ALLOW_DIAGNOSTICS
                0042 C     !FUNCTIONS:
                0043       LOGICAL  DIAGNOSTICS_IS_ON
                0044       EXTERNAL DIAGNOSTICS_IS_ON
                0045 #endif /* ALLOW_DIAGNOSTICS */
                0046
23cba332b3 Jean*0047 C     !INPUT/OUTPUT PARAMETERS:
                0048 C   uIce   (inp) :: zonal      ice velocity (input)
                0049 C   vIce   (inp) :: meridional ice velocity (input)
9da8ce8edb Jean*0050 C   icFrac (inp) :: seaice fraction (input)
ec0d7df165 Mart*0051 C   SImaskU(inp) :: mask at U-point
                0052 C   SImaskV(inp) :: mask at V-point
23cba332b3 Jean*0053 C   taux   (out) :: zonal      wind stress over ice at U point
                0054 C   tauy   (out) :: meridional wind stress over ice at V point
                0055 C   myTime (inp) :: current time in simulation
                0056 C   myIter (inp) :: iteration number in simulation
                0057 C   myThid (inp) :: my Thread Id. number
9da8ce8edb Jean*0058       _RL uIce   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0059       _RL vIce   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0060       _RL icFrac (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
ec0d7df165 Mart*0061       _RL SImaskU(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0062       _RL SImaskV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
9da8ce8edb Jean*0063       _RL taux   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0064       _RL tauy   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0dfb864288 Mart*0065       _RL     myTime
                0066       INTEGER myIter
                0067       INTEGER myThid
23cba332b3 Jean*0068 CEOP
0dfb864288 Mart*0069
23cba332b3 Jean*0070 C     !LOCAL VARIABLES:
                0071 C     i,j,bi,bj :: Loop counters
                0072 C     ks        :: vertical index of surface layer
0dfb864288 Mart*0073       INTEGER bi, bj, i, j
23cba332b3 Jean*0074       INTEGER ks
35fda33b05 Jean*0075       _RL  COSWIN
                0076       _RS  SINWIN
8a6cc4f42a Jean*0077 C     CDAIR   :: local wind stress coefficient (used twice)
                0078 C     oceTauX :: wind-stress over open-ocean (on Arakawa A-grid), X direction
                0079 C     oceTauY :: wind-stress over open-ocean (on Arakawa A-grid), Y direction
9da8ce8edb Jean*0080       _RL CDAIR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
56d13a40ed Mart*0081 #ifdef ALLOW_EXF
                0082       _RL AAA
                0083       _RL uTmp   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0084       _RL vTmp   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0085 # ifndef SEAICE_EXTERNAL_FLUXES
                0086       _RL U1, V1
9da8ce8edb Jean*0087       _RL oceTauX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0088       _RL oceTauY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
56d13a40ed Mart*0089 # endif
                0090 #endif /* ALLOW_EXF */
9da8ce8edb Jean*0091 #ifdef ALLOW_DIAGNOSTICS
                0092       _RL locVar (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0093       _RL locFrac
                0094 #endif /* ALLOW_DIAGNOSTICS */
0dfb864288 Mart*0095
23cba332b3 Jean*0096 C--   surface level
                0097       ks = 1
56d13a40ed Mart*0098       IF ( usingPcoords ) ks = Nr
0dfb864288 Mart*0099 C--   introduce turning angle (default is zero)
                0100       SINWIN=SIN(SEAICE_airTurnAngle*deg2rad)
                0101       COSWIN=COS(SEAICE_airTurnAngle*deg2rad)
                0102
                0103 C--   NOW SET UP FORCING FIELDS
                0104
                0105       DO bj=myByLo(myThid),myByHi(myThid)
                0106        DO bi=myBxLo(myThid),myBxHi(myThid)
58a71b9205 Dimi*0107
8bf64a0174 Jean*0108 #ifdef ALLOW_EXF
                0109 C--   Wind stress is computed on center of C-grid cell
                0110 C     and interpolated to U and V points later
                0111
58a71b9205 Dimi*0112 #ifndef SEAICE_EXTERNAL_FLUXES
                0113 C--   First compute wind-stress over open ocean: this will results in
                0114 C     over-writing fu and fv that were computed or read-in by pkg/exf.
8bf64a0174 Jean*0115         DO j=1-OLy,sNy+OLy
                0116          DO i=1-OLx,sNx+OLx
203c016029 Jean*0117           U1=UWIND(i,j,bi,bj)
                0118           V1=VWIND(i,j,bi,bj)
0dfb864288 Mart*0119           AAA=U1**2+V1**2
                0120           IF ( AAA .LE. SEAICE_EPS_SQ ) THEN
                0121              AAA=SEAICE_EPS
                0122           ELSE
                0123              AAA=SQRT(AAA)
                0124           ENDIF
203c016029 Jean*0125           CDAIR(i,j)=SEAICE_rhoAir*OCEAN_drag
0dfb864288 Mart*0126      &         *(2.70 _d 0+0.142 _d 0*AAA+0.0764 _d 0*AAA*AAA)
203c016029 Jean*0127           oceTauX(i,j)=CDAIR(i,j)*
                0128      &         (COSWIN*U1-SIGN(SINWIN, _fCori(i,j,bi,bj))*V1)
                0129           oceTauY(i,j)=CDAIR(i,j)*
                0130      &         (SIGN(SINWIN, _fCori(i,j,bi,bj))*U1+COSWIN*V1)
0dfb864288 Mart*0131          ENDDO
                0132         ENDDO
8a6cc4f42a Jean*0133 C--   Interpolate wind stress over open ocean (N/m^2)
35fda33b05 Jean*0134 C     from A-grid to U and V points of C-grid
8bf64a0174 Jean*0135         DO j=1-OLy+1,sNy+OLy
                0136          DO i=1-OLx+1,sNx+OLx
203c016029 Jean*0137           fu(i,j,bi,bj) = 0.5 _d 0*( oceTauX(i,j) + oceTauX(i-1,j) )
ec0d7df165 Mart*0138      &                            * SIMaskU(i,j,bi,bj)
203c016029 Jean*0139           fv(i,j,bi,bj) = 0.5 _d 0*( oceTauY(i,j) + oceTauY(i,j-1) )
ec0d7df165 Mart*0140      &                            * SIMaskV(i,j,bi,bj)
8a6cc4f42a Jean*0141          ENDDO
                0142         ENDDO
                0143 #endif /* ndef SEAICE_EXTERNAL_FLUXES */
                0144
8bf64a0174 Jean*0145 # ifdef SEAICE_EXTERNAL_FLUXES
5650936ed9 Jean*0146         IF ( useEXF .AND. useAtmWind ) THEN
8bf64a0174 Jean*0147 # else  /* SEAICE_EXTERNAL_FLUXES */
5650936ed9 Jean*0148         IF ( useEXF ) THEN
8bf64a0174 Jean*0149 # endif /* SEAICE_EXTERNAL_FLUXES */
56d13a40ed Mart*0150 C     make local copies of wind field
                0151          DO j=1-OLy,sNy+OLy
                0152           DO i=1-OLx,sNx+OLx
                0153            uTmp(i,j)=uwind(i,j,bi,bj)
                0154            vTmp(i,j)=vwind(i,j,bi,bj)
                0155           ENDDO
                0156          ENDDO
                0157          IF ( useRelativeWind ) THEN
8377b8ee87 Mart*0158 C     subtract ice velocities from each component wind-speed
5650936ed9 Jean*0159           DO j=1-OLy,sNy+OLy-1
                0160            DO i=1-OLx,sNx+OLx-1
56d13a40ed Mart*0161             uTmp(i,j)=uwind(i,j,bi,bj)
203c016029 Jean*0162      &          - 0.5 _d 0 * (uIce(i,j,bi,bj)+uIce(i+1,j,bi,bj))
56d13a40ed Mart*0163             vTmp(i,j)=vwind(i,j,bi,bj)
203c016029 Jean*0164      &          - 0.5 _d 0 * (vIce(i,j,bi,bj)+vIce(i,j+1,bi,bj))
5650936ed9 Jean*0165            ENDDO
58a71b9205 Dimi*0166           ENDDO
5650936ed9 Jean*0167          ENDIF
56d13a40ed Mart*0168 C--   Now compute ice surface stress
                0169          DO j=1-OLy,sNy+OLy
                0170           DO i=1-OLx,sNx+OLx
                0171            AAA=uTmp(i,j)**2+vTmp(i,j)**2
                0172            IF ( AAA .LE. SEAICE_EPS_SQ ) THEN
                0173             AAA=SEAICE_EPS
                0174            ELSE
                0175             AAA=SQRT(AAA)
                0176            ENDIF
                0177            IF ( yC(i,j,bi,bj) .LT. ZERO ) THEN
                0178             CDAIR(i,j) = SEAICE_rhoAir*SEAICE_drag_south*AAA
                0179            ELSE
                0180             CDAIR(i,j) = SEAICE_rhoAir*SEAICE_drag*AAA
                0181            ENDIF
58a71b9205 Dimi*0182           ENDDO
56d13a40ed Mart*0183          ENDDO
                0184          DO j=1-OLy+1,sNy+OLy
                0185           DO i=1-OLx+1,sNx+OLx
58a71b9205 Dimi*0186 C     interpolate to U points
56d13a40ed Mart*0187            taux(i,j,bi,bj)=0.5 _d 0 *
23cba332b3 Jean*0188      &         (  CDAIR(i  ,j)*(
56d13a40ed Mart*0189      &          COSWIN                            *uTmp(i  ,j)
                0190      &          -SIGN(SINWIN, _fCori(i  ,j,bi,bj))*vTmp(i  ,j) )
203c016029 Jean*0191      &          + CDAIR(i-1,j)*(
56d13a40ed Mart*0192      &          COSWIN                            *uTmp(i-1,j)
                0193      &          -SIGN(SINWIN, _fCori(i-1,j,bi,bj))*vTmp(i-1,j) )
ec0d7df165 Mart*0194      &         )*SIMaskU(i,j,bi,bj)
58a71b9205 Dimi*0195 C     interpolate to V points
56d13a40ed Mart*0196            tauy(i,j,bi,bj)=0.5 _d 0 *
23cba332b3 Jean*0197      &         (  CDAIR(i,j  )*(
56d13a40ed Mart*0198      &          SIGN(SINWIN, _fCori(i,j  ,bi,bj))*uTmp(i,j  )
                0199      &          +COSWIN                          *vTmp(i,j  ) )
203c016029 Jean*0200      &          + CDAIR(i,j-1)*(
56d13a40ed Mart*0201      &          SIGN(SINWIN, _fCori(i,j-1,bi,bj))*uTmp(i,j-1)
                0202      &          +COSWIN                          *vTmp(i,j-1) )
ec0d7df165 Mart*0203      &         )*SIMaskV(i,j,bi,bj)
58a71b9205 Dimi*0204           ENDDO
56d13a40ed Mart*0205          ENDDO
58a71b9205 Dimi*0206
8bf64a0174 Jean*0207         ELSE
5650936ed9 Jean*0208 #else  /* ALLOW_EXF */
                0209         IF (.TRUE.) THEN
                0210 #endif /* ALLOW_EXF */
8bf64a0174 Jean*0211
855e83a213 Jean*0212 C--   Wind stress is available on U and V points, copy it to seaice variables.
                0213 #ifdef HACK_FOR_GMAO_CPL
                0214          DO j=2-OLy,sNy+OLy
                0215           DO i=2-OLx,sNx+OLx
ec0d7df165 Mart*0216            taux(i,j,bi,bj) = SIwindTauX(i,j,bi,bj)*SIMaskU(i,j,bi,bj)
                0217            tauy(i,j,bi,bj) = SIwindTauY(i,j,bi,bj)*SIMaskV(i,j,bi,bj)
855e83a213 Jean*0218           ENDDO
                0219          ENDDO
                0220 #else /* HACK_FOR_GMAO_CPL */
5650936ed9 Jean*0221          DO j=1-OLy,sNy+OLy
                0222           DO i=1-OLx,sNx+OLx
0dfb864288 Mart*0223 C now ice surface stress
5650936ed9 Jean*0224            IF ( yC(i,j,bi,bj) .LT. ZERO ) THEN
                0225             CDAIR(i,j) = SEAICE_drag_south/OCEAN_drag
                0226            ELSE
                0227             CDAIR(i,j) = SEAICE_drag      /OCEAN_drag
                0228            ENDIF
                0229            taux(i,j,bi,bj) = CDAIR(i,j)*fu(i,j,bi,bj)
ec0d7df165 Mart*0230      &                     *SIMaskU(i,j,bi,bj)
5650936ed9 Jean*0231            tauy(i,j,bi,bj) = CDAIR(i,j)*fv(i,j,bi,bj)
ec0d7df165 Mart*0232      &                     *SIMaskV(i,j,bi,bj)
5650936ed9 Jean*0233           ENDDO
0dfb864288 Mart*0234          ENDDO
855e83a213 Jean*0235 #endif /* HACK_FOR_GMAO_CPL */
8bf64a0174 Jean*0236
                0237         ENDIF
                0238
                0239 C--   end bi,bj loops
0dfb864288 Mart*0240        ENDDO
                0241       ENDDO
8bf64a0174 Jean*0242
9da8ce8edb Jean*0243 #ifdef ALLOW_DIAGNOSTICS
                0244       IF ( useDiagnostics ) THEN
                0245        CALL DIAGNOSTICS_FILL( taux, 'SItaux  ', 0,1, 0,1,1, myThid )
                0246        CALL DIAGNOSTICS_FILL( tauy, 'SItauy  ', 0,1, 0,1,1, myThid )
                0247        IF ( DIAGNOSTICS_IS_ON( 'SIatmTx ', myThid ) ) THEN
                0248         DO bj=myByLo(myThid),myByHi(myThid)
                0249          DO bi=myBxLo(myThid),myBxHi(myThid)
                0250            DO j=2-OLy,sNy+OLy
                0251             DO i=2-OLx,sNx+OLx
                0252               locFrac = ( icFrac( i, j,bi,bj)
                0253      &                  + icFrac(i-1,j,bi,bj) )*halfRL
                0254               locVar(i,j) = taux(i,j,bi,bj)*locFrac
                0255      &                      + fu(i,j,bi,bj)*(oneRL-locFrac)
                0256             ENDDO
                0257            ENDDO
                0258            CALL DIAGNOSTICS_FILL(locVar,'SIatmTx ',0,1,2,bi,bj,myThid)
                0259          ENDDO
                0260         ENDDO
                0261        ENDIF
                0262        IF ( DIAGNOSTICS_IS_ON( 'SIatmTy ', myThid ) ) THEN
                0263         DO bj=myByLo(myThid),myByHi(myThid)
                0264          DO bi=myBxLo(myThid),myBxHi(myThid)
                0265            DO j=2-OLy,sNy+OLy
                0266             DO i=2-OLx,sNx+OLx
                0267               locFrac = ( icFrac(i, j, bi,bj)
                0268      &                  + icFrac(i,j-1,bi,bj) )*halfRL
                0269               locVar(i,j) = tauy(i,j,bi,bj)*locFrac
                0270      &                      + fv(i,j,bi,bj)*(oneRL-locFrac)
                0271             ENDDO
                0272            ENDDO
                0273            CALL DIAGNOSTICS_FILL(locVar,'SIatmTy ',0,1,2,bi,bj,myThid)
                0274          ENDDO
                0275         ENDDO
                0276        ENDIF
                0277       ENDIF
                0278 #endif /* ALLOW_DIAGNOSTICS */
                0279
0dfb864288 Mart*0280       RETURN
                0281       END