Back to home page

MITgcm
 
 

    

File indexing completed on 2024-03-02 06:10:28 UTC

view on githubraw file Latest commit 5cf43646 on 2024-03-01 18:50:49 UTC
3752238fd8 Patr*0001 #include "EXF_OPTIONS.h"
3a255f48df Gael*0002 #ifdef ALLOW_CTRL
                0003 # include "CTRL_OPTIONS.h"
                0004 #endif
                0005 #ifdef ALLOW_AUTODIFF
                0006 # include "AUTODIFF_OPTIONS.h"
                0007 #endif
3752238fd8 Patr*0008 
c5e3091043 Jean*0009       SUBROUTINE EXF_WIND( myTime, myIter, myThid )
3752238fd8 Patr*0010 
c5e3091043 Jean*0011 C     ==================================================================
                0012 C     SUBROUTINE exf_wind
                0013 C     ==================================================================
                0014 C
                0015 C     o Prepare wind speed and stress fields
                0016 C
                0017 C     ==================================================================
                0018 C     SUBROUTINE exf_wind
                0019 C     ==================================================================
3752238fd8 Patr*0020 
c5e3091043 Jean*0021       IMPLICIT NONE
3752238fd8 Patr*0022 
c5e3091043 Jean*0023 C     == global variables ==
3752238fd8 Patr*0024 
                0025 #include "SIZE.h"
50113055f4 Jean*0026 #include "EEPARAMS.h"
3752238fd8 Patr*0027 #include "PARAMS.h"
                0028 
082e18c36c Jean*0029 #include "EXF_PARAM.h"
                0030 #include "EXF_FIELDS.h"
                0031 #include "EXF_CONSTANTS.h"
56d13a40ed Mart*0032 #include "DYNVARS.h"
3752238fd8 Patr*0033 
                0034 #ifdef ALLOW_AUTODIFF_TAMC
                0035 #include "tamc.h"
                0036 #endif
423768d890 Jean*0037 #ifdef ALLOW_GENTIM2D_CONTROL
6a3fbcbfe3 Gael*0038 # include "CTRL_SIZE.h"
5cf4364659 Mart*0039 # include "CTRL.h"
6a3fbcbfe3 Gael*0040 # include "CTRL_GENARR.h"
                0041 #endif
3752238fd8 Patr*0042 
c5e3091043 Jean*0043 C     == routine arguments ==
3752238fd8 Patr*0044 
c5e3091043 Jean*0045       _RL     myTime
                0046       INTEGER myIter
                0047       INTEGER myThid
3752238fd8 Patr*0048 
c5e3091043 Jean*0049 C     == external functions ==
3752238fd8 Patr*0050 
c5e3091043 Jean*0051 C     == local variables ==
3752238fd8 Patr*0052 
c5e3091043 Jean*0053       INTEGER bi,bj
56d13a40ed Mart*0054       INTEGER i,j,ks
                0055       _RL     urelw(1:sNx,1:sNy)
                0056       _RL     vrelw(1:sNx,1:sNy)
                0057       _RL     wsLoc(1:sNx,1:sNy)
c5e3091043 Jean*0058       _RL     wsSq
423768d890 Jean*0059 #ifdef ALLOW_BULKFORMULAE
c5e3091043 Jean*0060       _RL     usSq, recip_sqrtRhoA, ustar
                0061       _RL     tmp1, tmp2, tmp3, tmp4
423768d890 Jean*0062 #endif /* ALLOW_BULKFORMULAE */
c5e3091043 Jean*0063       _RL     oneThirdRL
                0064       PARAMETER ( oneThirdRL = 1.d0 / 3.d0 )
423768d890 Jean*0065 #if !(defined ALLOW_BULKFORMULAE) || !(defined ALLOW_ATM_TEMP)
                0066       _RL wsm, tmpbulk
                0067 #endif
                0068 #ifdef ALLOW_GENTIM2D_CONTROL
6a3fbcbfe3 Gael*0069       INTEGER iarr
                0070 #endif
7c50f07931 Mart*0071 #ifdef ALLOW_AUTODIFF_TAMC
                0072       INTEGER ikey
                0073 #endif
c5e3091043 Jean*0074 C     == end of interface ==
3752238fd8 Patr*0075 
c5e3091043 Jean*0076 C--   Use atmospheric state to compute surface fluxes.
3752238fd8 Patr*0077 
56d13a40ed Mart*0078       ks = 1
                0079       IF ( usingPCoords ) ks = Nr
                0080 
c5e3091043 Jean*0081 C     Loop over tiles.
                0082       DO bj = myByLo(myThid),myByHi(myThid)
                0083        DO bi = myBxLo(myThid),myBxHi(myThid)
3752238fd8 Patr*0084 
ddea17ebdb Gael*0085 #ifdef ALLOW_AUTODIFF_TAMC
20dee61641 Mart*0086         ikey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
ddea17ebdb Gael*0087 #endif /* ALLOW_AUTODIFF_TAMC */
                0088 
c5e3091043 Jean*0089 C--   Initialise
                0090         DO j = 1,sNy
                0091          DO i = 1,sNx
                0092           wsLoc(i,j) = 0. _d 0
56d13a40ed Mart*0093           urelw(i,j) = uwind(i,j,bi,bj)
                0094           vrelw(i,j) = vwind(i,j,bi,bj)
3752238fd8 Patr*0095           cw(i,j,bi,bj) = 0. _d 0
                0096           sw(i,j,bi,bj) = 0. _d 0
                0097           sh(i,j,bi,bj) = 0. _d 0
c5e3091043 Jean*0098           wStress(i,j,bi,bj) = 0. _d 0
                0099          ENDDO
                0100         ENDDO
                0101 
56d13a40ed Mart*0102         IF (useRelativeWind) THEN
                0103 C     Subtract uVel and vVel from uwind and vwind.
                0104          DO j = 1,sNy
                0105           DO i = 1,sNx
                0106            urelw(i,j) = uwind(i,j,bi,bj) - 0.5 _d 0
                0107      &          * (uVel(i,j,ks,bi,bj)+uVel(i+1,j,ks,bi,bj))
                0108            vrelw(i,j) = vwind(i,j,bi,bj) - 0.5 _d 0
                0109      &          * (vVel(i,j,ks,bi,bj)+vVel(i,j+1,ks,bi,bj))
                0110           ENDDO
                0111          ENDDO
                0112 #ifdef ALLOW_DIAGNOSTICS
                0113          IF ( useDiagnostics ) THEN
                0114           CALL DIAGNOSTICS_FILL( urelw,'EXFurelw',0,1,3,bi,bj,myThid )
                0115           CALL DIAGNOSTICS_FILL( vrelw,'EXFvrelw',0,1,3,bi,bj,myThid )
                0116          ENDIF
                0117 #endif
                0118         ENDIF
                0119 
358649780a Gael*0120 #ifdef ALLOW_AUTODIFF_TAMC
56d13a40ed Mart*0121 CADJ STORE urelw (:,:)       = comlev1_bibj, key=ikey, byte=isbyte
                0122 CADJ STORE vrelw (:,:)       = comlev1_bibj, key=ikey, byte=isbyte
358649780a Gael*0123 CADJ STORE wspeed(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
                0124 #endif
                0125 
                0126         IF ( useAtmWind ) THEN
c5e3091043 Jean*0127 
                0128 C--   Wind speed and direction.
423768d890 Jean*0129          DO j = 1,sNy
                0130           DO i = 1,sNx
56d13a40ed Mart*0131            wsSq = urelw(i,j)*urelw(i,j)
                0132      &          + vrelw(i,j)*vrelw(i,j)
c5e3091043 Jean*0133            IF ( wsSq .NE. 0. _d 0 ) THEN
                0134              wsLoc(i,j) = SQRT(wsSq)
56d13a40ed Mart*0135              cw(i,j,bi,bj) = urelw(i,j)/wsLoc(i,j)
                0136              sw(i,j,bi,bj) = vrelw(i,j)/wsLoc(i,j)
c5e3091043 Jean*0137            ELSE
                0138              wsLoc(i,j) = 0. _d 0
3752238fd8 Patr*0139              cw(i,j,bi,bj) = 0. _d 0
                0140              sw(i,j,bi,bj) = 0. _d 0
c5e3091043 Jean*0141            ENDIF
423768d890 Jean*0142           ENDDO
c5e3091043 Jean*0143          ENDDO
423768d890 Jean*0144          IF ( wspeedfile .EQ. ' ' ) THEN
c5e3091043 Jean*0145 C-    wind-speed is not loaded from file: save local array into common block
                0146           DO j = 1,sNy
                0147            DO i = 1,sNx
                0148              wspeed(i,j,bi,bj) = wsLoc(i,j)
                0149            ENDDO
                0150           ENDDO
423768d890 Jean*0151          ENDIF
c5e3091043 Jean*0152 
423768d890 Jean*0153 #ifdef ALLOW_BULKFORMULAE
358649780a Gael*0154         ELSE
423768d890 Jean*0155 C--     case useAtmWind=F
c5e3091043 Jean*0156 
                0157 C--   Wind stress and direction.
423768d890 Jean*0158          DO j = 1,sNy
                0159           DO i = 1,sNx
c5e3091043 Jean*0160            IF ( stressIsOnCgrid ) THEN
                0161              usSq = ( ustress(i,  j,bi,bj)*ustress(i  ,j,bi,bj)
                0162      &               +ustress(i+1,j,bi,bj)*ustress(i+1,j,bi,bj)
                0163      &               +vstress(i,j,  bi,bj)*vstress(i,j  ,bi,bj)
                0164      &               +vstress(i,j+1,bi,bj)*vstress(i,j+1,bi,bj)
                0165      &              )*0.5 _d 0
                0166            ELSE
                0167              usSq = ustress(i,j,bi,bj)*ustress(i,j,bi,bj)
                0168      &             +vstress(i,j,bi,bj)*vstress(i,j,bi,bj)
                0169            ENDIF
                0170            IF ( usSq .NE. 0. _d 0 ) THEN
                0171              wStress(i,j,bi,bj) = SQRT(usSq)
                0172 c            ustar = SQRT(usSq/atmrho)
                0173              cw(i,j,bi,bj) = ustress(i,j,bi,bj)/wStress(i,j,bi,bj)
                0174              sw(i,j,bi,bj) = vstress(i,j,bi,bj)/wStress(i,j,bi,bj)
                0175            ELSE
                0176              wStress(i,j,bi,bj) = 0. _d 0
                0177              cw(i,j,bi,bj)      = 0. _d 0
                0178              sw(i,j,bi,bj)      = 0. _d 0
                0179            ENDIF
423768d890 Jean*0180           ENDDO
c5e3091043 Jean*0181          ENDDO
                0182 
423768d890 Jean*0183          IF ( wspeedfile .EQ. ' ' ) THEN
0fcb643df7 Mart*0184 C--   wspeed is not loaded ; derive wind-speed by inversion of
c5e3091043 Jean*0185 C     wind-stress=fct(wind-speed) relation:
                0186 C             The variables us, sh and rdn have to be computed from
                0187 C             given wind stresses inverting relationship for neutral
                0188 C             drag coeff. cdn.
                0189 C             The inversion is based on linear and quadratic form of
                0190 C             cdn(umps); ustar can be directly computed from stress;
423768d890 Jean*0191           recip_sqrtRhoA = 1. _d 0 / SQRT(atmrho)
                0192           DO j = 1,sNy
                0193            DO i = 1,sNx
7c50f07931 Mart*0194 C     check for zero wStress to please AD tools
502f242cbd Mart*0195             IF ( wStress(i,j,bi,bj) .LE. 0. _d 0 ) THEN
                0196              ustar      = 0. _d 0
                0197              wsloc(i,j) = 0. _d 0
c5e3091043 Jean*0198             ELSE
502f242cbd Mart*0199              ustar = SQRT(wStress(i,j,bi,bj))*recip_sqrtRhoA
                0200              IF ( ustar .LT. ustofu11 ) THEN
                0201               tmp1 = -cquadrag_2/cquadrag_1*exf_half
                0202               tmp2 = SQRT(tmp1*tmp1 + ustar*ustar/cquadrag_1)
                0203               wsLoc(i,j) = SQRT(tmp1 + tmp2)
                0204              ELSE
                0205               tmp1 = clindrag_2/clindrag_1*oneThirdRL
                0206               tmp2 = ustar*ustar/clindrag_1*exf_half
                0207      &             - tmp1*tmp1*tmp1
                0208               tmp3 = SQRT( ustar*ustar/clindrag_1*
                0209      &             (ustar*ustar/clindrag_1*0.25 _d 0 - tmp1*tmp1*tmp1 )
                0210      &                    )
                0211               tmp4 = (tmp2 + tmp3)**oneThirdRL
                0212               wsLoc(i,j) = tmp4 + tmp1*tmp1 / tmp4 - tmp1
                0213 c             wsLoc(i,j) = (tmp2 + tmp3)**oneThirdRL +
                0214 c    &             tmp1*tmp1 * (tmp2 + tmp3)**(-oneThirdRL) - tmp1
                0215              ENDIF
c5e3091043 Jean*0216             ENDIF
423768d890 Jean*0217            ENDDO
c5e3091043 Jean*0218           ENDDO
                0219 C-    save local array wind-speed to common block
423768d890 Jean*0220           DO j = 1,sNy
                0221            DO i = 1,sNx
c5e3091043 Jean*0222             wspeed(i,j,bi,bj) = wsLoc(i,j)
423768d890 Jean*0223            ENDDO
c5e3091043 Jean*0224           ENDDO
423768d890 Jean*0225 C-    end if wspeedfile = empty
                0226          ENDIF
c5e3091043 Jean*0227 
358649780a Gael*0228 #ifdef ALLOW_AUTODIFF_TAMC
                0229 CADJ STORE wspeed(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
                0230 CADJ STORE uwind (:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
                0231 CADJ STORE vwind (:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
                0232 CADJ STORE cw    (:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
                0233 CADJ STORE sw    (:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
                0234 #endif
                0235 
c5e3091043 Jean*0236 C--   infer wind field from wind-speed & wind-stress direction
423768d890 Jean*0237          DO j = 1,sNy
                0238           DO i = 1,sNx
c5e3091043 Jean*0239            uwind(i,j,bi,bj) = wspeed(i,j,bi,bj)*cw(i,j,bi,bj)
                0240            vwind(i,j,bi,bj) = wspeed(i,j,bi,bj)*sw(i,j,bi,bj)
423768d890 Jean*0241           ENDDO
c5e3091043 Jean*0242          ENDDO
423768d890 Jean*0243 
                0244 #endif /* ALLOW_BULKFORMULAE */
                0245 C--   end if/else useAtmWind
358649780a Gael*0246         ENDIF
3752238fd8 Patr*0247 
423768d890 Jean*0248 #ifdef ALLOW_GENTIM2D_CONTROL
6a3fbcbfe3 Gael*0249         DO j = 1,sNy
                0250          DO i = 1,sNx
                0251            do iarr = 1, maxCtrlTim2D
d15205fe5b Gael*0252            if (xx_gentim2d_file(iarr)(1:9).EQ.'xx_wspeed')
6a3fbcbfe3 Gael*0253      &       wspeed(i,j,bi,bj)=wspeed(i,j,bi,bj)+
                0254      &                         xx_gentim2d(i,j,bi,bj,iarr)
                0255            enddo
                0256          ENDDO
                0257         ENDDO
                0258 #endif
                0259 
ddea17ebdb Gael*0260 #ifdef ALLOW_AUTODIFF_TAMC
                0261 CADJ STORE wspeed(:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
                0262 #endif
                0263 
c5e3091043 Jean*0264 C--   set wind-speed lower limit
                0265         DO j = 1,sNy
                0266          DO i = 1,sNx
                0267            sh(i,j,bi,bj) = MAX(wspeed(i,j,bi,bj),uMin)
                0268          ENDDO
                0269         ENDDO
                0270 
423768d890 Jean*0271 #if !(defined ALLOW_BULKFORMULAE) || !(defined ALLOW_ATM_TEMP)
                0272 C     Note: In case ALLOW_BULKFORMULAE and ALLOW_ATM_TEMP are defined,
                0273 C     wind-stress is computed (if needed) within S/R EXF_BULKFORMULAE
                0274         IF ( useAtmWind ) THEN
bab1d47ef6 Jean*0275 c#ifdef ALLOW_ATM_WIND
423768d890 Jean*0276 C--   Computes wind-stress:
                0277          DO j = 1,sNy
                0278           DO i = 1,sNx
                0279            wsm     = sh(i,j,bi,bj)
cc60455fbb Mart*0280 # ifdef  ALLOW_DRAG_LARGEYEAGER09
                0281 C     Large and Yeager (2009), Climate Dynamics, equation 11a/b
                0282            tmpbulk =  cdrag_1/wsm + cdrag_2 + cdrag_3*wsm
                0283      &          + cdrag_8 * wsm**6
                0284            tmpbulk = exf_scal_BulkCdn * (
                0285      &            ( halfRL - SIGN(halfRL, wsm-umax) )*tmpbulk
                0286      &          + ( halfRL + SIGN(halfRL, wsm-umax) )*cdragMax
                0287      &          )
                0288 # else
423768d890 Jean*0289            tmpbulk = exf_scal_BulkCdn
                0290      &             * ( cdrag_1/wsm + cdrag_2 + cdrag_3*wsm )
cc60455fbb Mart*0291 # endif
423768d890 Jean*0292            ustress(i,j,bi,bj) = atmrho*tmpbulk*wspeed(i,j,bi,bj)
56d13a40ed Mart*0293      &                        * urelw(i,j)
423768d890 Jean*0294            vstress(i,j,bi,bj) = atmrho*tmpbulk*wspeed(i,j,bi,bj)
56d13a40ed Mart*0295      &                        * vrelw(i,j)
423768d890 Jean*0296           ENDDO
                0297          ENDDO
bab1d47ef6 Jean*0298 c#else /* ALLOW_ATM_WIND */
                0299 c        STOP 'ABNORMAL END: S/R EXF_WIND: missing code for useAtmWind'
                0300 c#endif /* ALLOW_ATM_WIND */
423768d890 Jean*0301 C--   end if useAtmWind
                0302         ENDIF
                0303 #endif /* ndef ALLOW_BULKFORMULAE or ndef ALLOW_ATM_TEMP */
                0304 
c5e3091043 Jean*0305 C--   end bi,bj loops
                0306        ENDDO
                0307       ENDDO
                0308 
                0309       RETURN
                0310       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

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

The LXR team
Valid CSS 2.1! Valid HTML 4.01!