Back to home page

MITgcm
 
 

    

File indexing completed on 2023-09-03 05:10:27 UTC

view on githubraw file Latest commit 74487008 on 2023-09-03 01:50:18 UTC
6af7e0efdf Mart*0001 #include "KPP_OPTIONS.h"
1f89baba18 Patr*0002 #ifdef ALLOW_SALT_PLUME
                0003 #include "SALT_PLUME_OPTIONS.h"
                0004 #endif
6af7e0efdf Mart*0005 
                0006 CBOP
                0007 C !ROUTINE: KPP_FORCING_SURF
                0008 
                0009 C !INTERFACE: ==========================================================
                0010       SUBROUTINE KPP_FORCING_SURF(
                0011      I     rhoSurf, surfForcU, surfForcV,
29e16c9d38 Jean*0012      I     surfForcT, surfForcS, surfForcTice,
                0013      I     Qsw,
a67734cba4 Mart*0014 #ifdef KPP_ESTIMATE_UREF
                0015      I     dbloc,
                0016 #endif /* KPP_ESTIMATE_UREF */
63ceaaa79c Dimi*0017 #ifdef ALLOW_SALT_PLUME
1f89baba18 Patr*0018      I     SPforcS,SPforcT,
63ceaaa79c Dimi*0019 #endif /* ALLOW_SALT_PLUME */
29e16c9d38 Jean*0020      I     ttalpha, ssbeta,
                0021      O     ustar, bo, bosol,
63ceaaa79c Dimi*0022 #ifdef ALLOW_SALT_PLUME
                0023      O     boplume,
                0024 #endif /* ALLOW_SALT_PLUME */
                0025      O     dVsq,
edb6656069 Mart*0026      I     ikey, iMin, iMax, jMin, jMax, bi, bj, myTime, myThid )
6af7e0efdf Mart*0027 
                0028 C !DESCRIPTION: \bv
29e16c9d38 Jean*0029 C     *==========================================================*
6af7e0efdf Mart*0030 C     | SUBROUTINE KPP_FORCING_SURF                              |
                0031 C     | o Compute all surface related KPP fields:                |
                0032 C     |   - friction velocity ustar                              |
                0033 C     |   - turbulent and radiative surface buoyancy forcing,    |
63ceaaa79c Dimi*0034 C     |     bo and bosol, and surface haline buoyancy forcing    |
                0035 C     |     boplume                                              |
6af7e0efdf Mart*0036 C     |   - velocity shear relative to surface squared (this is  |
                0037 C     |     not really a surface affected quantity unless it is  |
                0038 C     |     computed with respect to some resolution independent |
                0039 C     |     reference level, that is KPP_ESTIMATE_UREF defined ) |
29e16c9d38 Jean*0040 C     *==========================================================*
6af7e0efdf Mart*0041       IMPLICIT NONE
                0042 
29e16c9d38 Jean*0043 C \ev
6af7e0efdf Mart*0044 
                0045 C !USES: ===============================================================
                0046 #include "SIZE.h"
                0047 #include "EEPARAMS.h"
                0048 #include "PARAMS.h"
                0049 #include "GRID.h"
                0050 #include "DYNVARS.h"
                0051 #include "KPP_PARAMS.h"
                0052 
                0053 C !INPUT PARAMETERS: ===================================================
                0054 C Routine arguments
edb6656069 Mart*0055 C     ikey   - key for storing trajectory for adjoint (taf)
29e16c9d38 Jean*0056 C     imin, imax, jmin, jmax  - array computation indices
6af7e0efdf Mart*0057 C     bi, bj - array indices on which to apply calculations
                0058 C     myTime - Current time in simulation
                0059 C     myThid - Current thread id
29e16c9d38 Jean*0060 C     rhoSurf- density of surface layer                            (kg/m^3)
6af7e0efdf Mart*0061 C     surfForcU     units are  r_unit.m/s^2 (=m^2/s^2 if r=z)
                0062 C     surfForcV     units are  r_unit.m/s^2 (=m^2/s^-2 if r=z)
ba0b047096 Mart*0063 C     surfForcS     units are  r_unit.g/kg/s (=g/kg.m/s if r=z)
6af7e0efdf Mart*0064 C            - EmPmR * S_surf plus salinity relaxation*drF(1)
                0065 C     surfForcT     units are  r_unit.Kelvin/s (=Kelvin.m/s if r=z)
                0066 C            - Qnet (+Qsw) plus temp. relaxation*drF(1)
                0067 C                -> calculate        -lambda*(T(model)-T(clim))
                0068 C            Qnet assumed to be net heat flux including ShortWave rad.
                0069 C     surfForcTice
                0070 C            - equivalent Temperature flux in the top level that corresponds
                0071 C              to the melting or freezing of sea-ice.
                0072 C              Note that the surface level temperature is modified
                0073 C              directly by the sea-ice model in order to maintain
                0074 C              water temperature under sea-ice at the freezing
                0075 C              point.  But we need to keep track of the
                0076 C              equivalent amount of heat that this surface-level
                0077 C              temperature change implies because it is used by
                0078 C              the KPP package (kpp_calc.F and kpp_transport_t.F).
                0079 C              Units are r_unit.K/s (=Kelvin.m/s if r=z) (>0 for ocean warming).
                0080 C
                0081 C     Qsw     - surface shortwave radiation (upwards positive)
63ceaaa79c Dimi*0082 C     saltPlumeFlux - salt rejected during freezing (downward = positive)
6af7e0efdf Mart*0083 C     ttalpha - thermal expansion coefficient without 1/rho factor
                0084 C               d(rho{k,k})/d(T(k))                           (kg/m^3/C)
                0085 C     ssbeta  - salt expansion coefficient without 1/rho factor
ba0b047096 Mart*0086 C               d(rho{k,k})/d(S(k))                    ((kg/m^3)/(g/kg))
29e16c9d38 Jean*0087 C !OUTPUT PARAMETERS:
6af7e0efdf Mart*0088 C     ustar  (nx,ny)       - surface friction velocity                  (m/s)
                0089 C     bo     (nx,ny)       - surface turbulent buoyancy forcing     (m^2/s^3)
                0090 C     bosol  (nx,ny)       - surface radiative buoyancy forcing     (m^2/s^3)
1f89baba18 Patr*0091 C     boplume(nx,ny,Nr+1)  - surface haline buoyancy forcing        (m^2/s^3)
6af7e0efdf Mart*0092 C     dVsq   (nx,ny,Nr)    - velocity shear re surface squared
                0093 C                            at grid levels for bldepth             (m^2/s^2)
                0094 
edb6656069 Mart*0095       INTEGER ikey
6af7e0efdf Mart*0096       INTEGER iMin, iMax, jMin, jMax
                0097       INTEGER bi, bj
                0098       INTEGER myThid
                0099       _RL     myTime
                0100 
                0101       _RL rhoSurf     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0102       _RL surfForcU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0103       _RL surfForcV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0104       _RL surfForcT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0105       _RL surfForcS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0106       _RL surfForcTice(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
ac5511770f Jean*0107       _RS Qsw         (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
6af7e0efdf Mart*0108       _RL TTALPHA     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nrp1)
                0109       _RL SSBETA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nrp1)
                0110 
                0111       _RL ustar ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
                0112       _RL bo    ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
                0113       _RL bosol ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
a67734cba4 Mart*0114 #ifdef KPP_ESTIMATE_UREF
                0115       _RL dbloc ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
                0116 #endif /* KPP_ESTIMATE_UREF */
63ceaaa79c Dimi*0117 #ifdef ALLOW_SALT_PLUME
1f89baba18 Patr*0118       _RL SPforcS     (1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr       )
                0119       _RL SPforcT     (1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr       )
                0120       _RL boplume     (1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1     )
63ceaaa79c Dimi*0121 #endif /* ALLOW_SALT_PLUME */
6af7e0efdf Mart*0122       _RL dVsq  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
                0123 
                0124 C !LOCAL VARIABLES: ====================================================
29e16c9d38 Jean*0125 C Local constants
                0126       _RL        p0    , p5    , p125
                0127       PARAMETER( p0=0.0, p5=0.5, p125=0.125 )
                0128       INTEGER i, j, k, im1, ip1, jm1, jp1
a67734cba4 Mart*0129       _RL tempvar2, epsLoc, epsLocSq
7fb657f65b Patr*0130       _RL recip_Cp
6af7e0efdf Mart*0131 
                0132       _RL work3 ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
                0133 
                0134 #ifdef KPP_ESTIMATE_UREF
a67734cba4 Mart*0135 c     z0     (nx,ny)       - Roughness length                             (m)
                0136 c     zRef   (nx,ny)       - Reference depth: Hmix * epsilon              (m)
                0137 c     uRef   (nx,ny)       - Reference zonal velocity                   (m/s)
                0138 c     vRef   (nx,ny)       - Reference meridional velocity              (m/s)
                0139       INTEGER kTmp
                0140       _RL tempvar1, dBdz1, dBdz2, ustarX, ustarY, zFac
6af7e0efdf Mart*0141       _RL z0    ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
                0142       _RL zRef  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
                0143       _RL uRef  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
                0144       _RL vRef  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
                0145 #endif
29e16c9d38 Jean*0146 #ifdef ALLOW_SALT_PLUME
                0147 #ifdef SALT_PLUME_VOLUME
                0148       INTEGER kp1
                0149       _RL temparray   (1-OLx:sNx+OLx, 1-OLy:sNy+OLy           )
                0150 #endif /* SALT_PLUME_VOLUME */
                0151 #endif /* ALLOW_SALT_PLUME */
6af7e0efdf Mart*0152 CEOP
                0153 
29e16c9d38 Jean*0154 C------------------------------------------------------------------------
                0155 C     friction velocity, turbulent and radiative surface buoyancy forcing
                0156 C     -------------------------------------------------------------------
                0157 C     taux / rho = surfForcU                               (N/m^2)
                0158 C     tauy / rho = surfForcV                               (N/m^2)
                0159 C     ustar = sqrt( sqrt( taux^2 + tauy^2 ) / rho )        (m/s)
                0160 C     bo    = - g * ( alpha*surfForcT +
                0161 C                     beta *surfForcS ) / rho              (m^2/s^3)
                0162 C     bosol = - g * alpha * Qsw * drF(1) / rho             (m^2/s^3)
                0163 C     boplume =-g * ( beta *saltPlumeFlux/rhoConst )/rho   (m^2/s^3)
                0164 C             =-g * ( beta *SPforcS      /rhoConst )/rho   (m^2/s^3)
                0165 C              -g * (alpha *SPforcT/Cp   /rhoConst )/rho   (m^2/s^3)
                0166 C------------------------------------------------------------------------
6af7e0efdf Mart*0167 
7fb657f65b Patr*0168       recip_Cp = 1. _d 0 / HeatCapacity_Cp
29e16c9d38 Jean*0169 C initialize arrays to zero
6af7e0efdf Mart*0170       DO j = 1-OLy, sNy+OLy
a67734cba4 Mart*0171        DO i = 1-OLx, sNx+OLx
                0172         ustar(i,j) = p0
2b4c90c108 Mart*0173         bo   (i,j) = p0
                0174         bosol(i,j) = p0
a67734cba4 Mart*0175        ENDDO
                0176       ENDDO
63ceaaa79c Dimi*0177 #ifdef ALLOW_SALT_PLUME
a67734cba4 Mart*0178       DO k = 1, Nrp1
                0179        DO j = 1-OLy, sNy+OLy
                0180         DO i = 1-OLx, sNx+OLx
2b4c90c108 Mart*0181          boplume(i,j,k) = p0
a67734cba4 Mart*0182         ENDDO
                0183        ENDDO
29e16c9d38 Jean*0184       ENDDO
a67734cba4 Mart*0185 #endif /* ALLOW_SALT_PLUME */
6af7e0efdf Mart*0186 
                0187       DO j = jmin, jmax
                0188        jp1 = j + 1
                0189        DO i = imin, imax
                0190         ip1 = i+1
                0191         work3(i,j) =
                0192      &   (surfForcU(i,j,bi,bj) + surfForcU(ip1,j,bi,bj)) *
                0193      &   (surfForcU(i,j,bi,bj) + surfForcU(ip1,j,bi,bj)) +
                0194      &   (surfForcV(i,j,bi,bj) + surfForcV(i,jp1,bi,bj)) *
                0195      &   (surfForcV(i,j,bi,bj) + surfForcV(i,jp1,bi,bj))
29e16c9d38 Jean*0196        ENDDO
                0197       ENDDO
a67734cba4 Mart*0198       epsLocSq = phepsi*phepsi*drF(1)*drF(1)
                0199       epsLoc   = SQRT ( p5*phepsi*drF(1) )
6af7e0efdf Mart*0200       DO j = jmin, jmax
                0201        DO i = imin, imax
                0202 
a67734cba4 Mart*0203         IF ( work3(i,j) .lt. epsLocSq ) THEN
                0204            ustar(i,j) = epsLoc
                0205         ELSE
6af7e0efdf Mart*0206            tempVar2 =  SQRT( work3(i,j) ) * p5
                0207            ustar(i,j) = SQRT( tempVar2 )
a67734cba4 Mart*0208         ENDIF
6af7e0efdf Mart*0209 
29e16c9d38 Jean*0210        ENDDO
                0211       ENDDO
6af7e0efdf Mart*0212 
a67734cba4 Mart*0213 C     I would like to replace 4 instances of 1/rhoSurf by multiplication
                0214 C     but this changes results in lab_sea.hb87/salt_plume, natl_box, and
                0215 C     seaice_obcs, so maybe later
                0216 CML      DO j = jmin, jmax
                0217 CML       DO i = imin, imax
                0218 CMLC     rhoSurf = dRho + rhoConst > 0 always, so we are safe here
2b4c90c108 Mart*0219 CML        work3(i,j) = 1. _d 0 / rhoSurf(i,j)
a67734cba4 Mart*0220 CML       ENDDO
                0221 CML      ENDDO
                0222 
6af7e0efdf Mart*0223       DO j = jmin, jmax
                0224        DO i = imin, imax
2b4c90c108 Mart*0225         bo(i,j) = - gravity *
                0226      &       ( TTALPHA(i,j,1) * (surfForcT(i,j,bi,bj)+
6af7e0efdf Mart*0227      &       surfForcTice(i,j,bi,bj)) +
2b4c90c108 Mart*0228      &       SSBETA(i,j,1) * surfForcS(i,j,bi,bj) )
                0229      &       / rhoSurf(i,j)
                0230 CML  &       * work3(i,j)
                0231         bosol(i,j) = gravity * TTALPHA(i,j,1) * Qsw(i,j,bi,bj) *
7fb657f65b Patr*0232      &       recip_Cp*recip_rhoConst
2b4c90c108 Mart*0233      &       / rhoSurf(i,j)
                0234 CML  &       * work3(i,j)
29e16c9d38 Jean*0235        ENDDO
                0236       ENDDO
30c6f5b1cd An T*0237 
                0238 #ifdef ALLOW_SALT_PLUME
1f89baba18 Patr*0239 Catn: need check sign of SPforcT
                0240 Cnote: on input, if notdef salt_plume_volume, SPforc[S,T](k>1)=!0
30c6f5b1cd An T*0241       IF ( useSALT_PLUME ) THEN
1f89baba18 Patr*0242 #ifdef SALT_PLUME_VOLUME
a67734cba4 Mart*0243        DO k = 1, Nr
                0244         kp1 = k+1
                0245         DO j = jmin, jmax
                0246          DO i = imin, imax
2b4c90c108 Mart*0247           temparray(i,j) = - gravity *
                0248      &            ( SSBETA(i,j,k) * SPforcS(i,j,k) +
                0249      &              TTALPHA(i,j,k)* SPforcT(i,j,k) / HeatCapacity_Cp )
                0250      &              * recip_rhoConst / rhoSurf(i,j)
                0251 CML  &              * recip_rhoConst * work3(i,j)
                0252           boplume(i,j,kp1) = boplume(i,j,k)+temparray(i,j)
29e16c9d38 Jean*0253          ENDDO
a67734cba4 Mart*0254         ENDDO
                0255        ENDDO
1f89baba18 Patr*0256 #else /* SALT_PLUME_VOLUME */
a67734cba4 Mart*0257        DO j = jmin, jmax
                0258         DO i = imin, imax
                0259 C     initialization has been done already and does not have to be
                0260 C     repeated here (and definitely not inside the i/j loops)
                0261 CML         DO k = 2, Nrp1
2b4c90c108 Mart*0262 CML          boplume(i,j,k) = 0. _d 0
a67734cba4 Mart*0263 CML         ENDDO
2b4c90c108 Mart*0264          boplume(i,j,1) = - gravity * SSBETA(i,j,1)
1f89baba18 Patr*0265      &              * SPforcS(i,j,1)
2b4c90c108 Mart*0266      &              * recip_rhoConst / rhoSurf(i,j)
                0267 CML  &              * recip_rhoConst * work3(i,j)
a67734cba4 Mart*0268         ENDDO
                0269        ENDDO
1f89baba18 Patr*0270 #endif /* SALT_PLUME_VOLUME */
30c6f5b1cd An T*0271       ENDIF
                0272 #endif /* ALLOW_SALT_PLUME */
                0273 
                0274 #ifdef ALLOW_DIAGNOSTICS
                0275       IF ( useDiagnostics ) THEN
a67734cba4 Mart*0276        CALL DIAGNOSTICS_FILL(bo     ,'KPPbo   ',0,1,2,bi,bj,myThid)
                0277        CALL DIAGNOSTICS_FILL(bosol  ,'KPPbosol',0,1,2,bi,bj,myThid)
30c6f5b1cd An T*0278 #ifdef ALLOW_SALT_PLUME
a67734cba4 Mart*0279        CALL DIAGNOSTICS_FILL(boplume,'KPPboplm',0,Nr,2,bi,bj,myThid)
30c6f5b1cd An T*0280 #endif /* ALLOW_SALT_PLUME */
                0281       ENDIF
                0282 #endif /* ALLOW_DIAGNOSTICS */
                0283 
29e16c9d38 Jean*0284 C------------------------------------------------------------------------
                0285 C     velocity shear
                0286 C     --------------
                0287 C     Get velocity shear squared, averaged from "u,v-grid"
                0288 C     onto "t-grid" (in (m/s)**2):
                0289 C     dVsq(k)=(Uref-U(k))**2+(Vref-V(k))**2      at grid levels
                0290 C------------------------------------------------------------------------
6af7e0efdf Mart*0291 
29e16c9d38 Jean*0292 C initialize arrays to zero
6af7e0efdf Mart*0293       DO k = 1, Nr
                0294        DO j = 1-OLy, sNy+OLy
                0295         DO i = 1-OLx, sNx+OLx
                0296          dVsq(i,j,k) = p0
29e16c9d38 Jean*0297         ENDDO
                0298        ENDDO
                0299       ENDDO
6af7e0efdf Mart*0300 
29e16c9d38 Jean*0301 C     dVsq computation
6af7e0efdf Mart*0302 
                0303 #ifdef KPP_ESTIMATE_UREF
                0304 
29e16c9d38 Jean*0305 C     Get rid of vertical resolution dependence of dVsq term by
                0306 C     estimating a surface velocity that is independent of first level
                0307 C     thickness in the model.  First determine mixed layer depth hMix.
                0308 C     Second zRef = espilon * hMix.  Third determine roughness length
                0309 C     scale z0.  Third estimate reference velocity.
a67734cba4 Mart*0310 CML   Comment: so far this code only handles the case of coarse resolution.
                0311 CML   High resolution, i.e. the case zref=hMix*epsilon > drF(1), is not
                0312 CML   handled by this code.
6af7e0efdf Mart*0313 
a67734cba4 Mart*0314       zFac = ABS(rF(3)) * LOG ( rF(3) / rF(2) ) * recip_drF(2)
6af7e0efdf Mart*0315       DO j = jmin, jmax
                0316        jp1 = j + 1
                0317        DO i = imin, imax
                0318         ip1 = i + 1
                0319 
29e16c9d38 Jean*0320 C     Determine mixed layer depth hMix as the shallowest depth at which
                0321 C     dB/dz exceeds 5.2e-5 s^-2.
a67734cba4 Mart*0322         kTmp = nzmax(i,j,bi,bj)
                0323         DO k = Nr, 1, -1
6af7e0efdf Mart*0324          IF ( k .LT. nzmax(i,j,bi,bj) .AND.
2b4c90c108 Mart*0325      &        maskC(i,j,k,bi,bj) .GT. 0. .AND.
a67734cba4 Mart*0326      &        dbloc(i,j,k) * recip_drC(k+1) .GT. dB_dz )
                0327      &        kTmp = k
6af7e0efdf Mart*0328         ENDDO
                0329 
29e16c9d38 Jean*0330 C     Linearly interpolate to find hMix.
a67734cba4 Mart*0331         k = kTmp
6af7e0efdf Mart*0332         IF ( k .EQ. 0 .OR. nzmax(i,j,bi,bj) .EQ. 1 ) THEN
                0333          zRef(i,j) = p0
                0334         ELSEIF ( k .EQ. 1) THEN
a67734cba4 Mart*0335          dBdz2 = dbloc(i,j,1) *recip_drC(2)
6af7e0efdf Mart*0336          zRef(i,j) = drF(1) * dB_dz / dBdz2
                0337         ELSEIF ( k .LT. nzmax(i,j,bi,bj) ) THEN
a67734cba4 Mart*0338          dBdz1 = dbloc(i,j,k-1) * recip_drC(k  )
                0339          dBdz2 = dbloc(i,j,k  ) * recip_drC(k+1)
                0340          zRef(i,j) = ABS(rF(k)) + drF(k) * (dB_dz - dBdz1) /
6af7e0efdf Mart*0341      &        MAX ( phepsi, dBdz2 - dBdz1 )
                0342         ELSE
a67734cba4 Mart*0343          zRef(i,j) = ABS(rF(k+1))
6af7e0efdf Mart*0344         ENDIF
29e16c9d38 Jean*0345 
                0346 C     Compute roughness length scale z0 subject to 0 < z0
6af7e0efdf Mart*0347         tempVar1 = p5 * (
                0348      &       (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  2,bi,bj)) *
                0349      &       (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  2,bi,bj)) +
                0350      &       (uVel(ip1,j,  1,bi,bj)-uVel(ip1,j,  2,bi,bj)) *
                0351      &       (uVel(ip1,j,  1,bi,bj)-uVel(ip1,j,  2,bi,bj)) +
                0352      &       (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  2,bi,bj)) *
29e16c9d38 Jean*0353      &       (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  2,bi,bj)) +
6af7e0efdf Mart*0354      &       (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,2,bi,bj)) *
                0355      &       (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,2,bi,bj)) )
                0356         IF ( tempVar1 .lt. (epsln*epsln) ) THEN
                0357          tempVar2 = epsln
                0358         ELSE
                0359          tempVar2 = SQRT ( tempVar1 )
                0360         ENDIF
a67734cba4 Mart*0361 C     ustar > sqrt(0.5*drF(1)*phepsi) by definition, so I do not know
                0362 C     why this extra max function call is required. Too small ustar
                0363 C     will lead to negative z0, which is caught in the next line
                0364         z0(i,j) = drF(1) * ( zFac - tempVar2 * vonK / ustar(i,j) )
                0365 CML     &       ( zFac - tempVar2 * vonK / MAX ( ustar(i,j), phepsi ) )
6af7e0efdf Mart*0366         z0(i,j) = MAX ( z0(i,j), phepsi )
                0367 
29e16c9d38 Jean*0368 C     zRef is set to 0.1 * hMix subject to z0 <= zRef <= drF(1)
6af7e0efdf Mart*0369         zRef(i,j) = MAX ( epsilon * zRef(i,j), z0(i,j) )
2b4c90c108 Mart*0370 C     with this statement the high resolution case zRef > drF(1) is ruled out
                0371 C     so we comment it out
                0372 CML        zRef(i,j) = MIN ( zRef(i,j), drF(1) )
29e16c9d38 Jean*0373 
                0374 C     Estimate reference velocity uRef and vRef.
6af7e0efdf Mart*0375         uRef(i,j) = p5 * ( uVel(i,j,1,bi,bj) + uVel(ip1,j,1,bi,bj) )
                0376         vRef(i,j) = p5 * ( vVel(i,j,1,bi,bj) + vVel(i,jp1,1,bi,bj) )
                0377         IF ( zRef(i,j) .LT. drF(1) ) THEN
29e16c9d38 Jean*0378          ustarX = ( surfForcU(i,  j,bi,bj) +
6af7e0efdf Mart*0379      &        surfForcU(ip1,j,bi,bj) ) * p5 *recip_drF(1)
                0380          ustarY = ( surfForcV(i,j,  bi,bj) +
                0381      &        surfForcV(i,jp1,bi,bj) ) * p5 *recip_drF(1)
                0382          tempVar1 = ustarX * ustarX + ustarY * ustarY
                0383          if ( tempVar1 .lt. (epsln*epsln) ) then
                0384           tempVar2 = epsln
                0385          else
                0386           tempVar2 = SQRT ( tempVar1 )
                0387          endif
                0388          tempVar2 = ustar(i,j) *
a67734cba4 Mart*0389      &        ( LOG ( zRef(i,j) * recip_drF(1) ) +
                0390      &        z0(i,j) / zRef(i,j) - z0(i,j) * recip_drF(1) ) /
6af7e0efdf Mart*0391      &        vonK / tempVar2
                0392          uRef(i,j) = uRef(i,j) + ustarX * tempVar2
                0393          vRef(i,j) = vRef(i,j) + ustarY * tempVar2
2b4c90c108 Mart*0394         ELSEIF ( zRef(i,j) .GE. drF(1) ) THEN
                0395 C     if zref > drF(1), average over zref
                0396          uRef(i,j) = uRef(i,j)*drF(1)
                0397          vRef(i,j) = vRef(i,j)*drF(1)
                0398          k = 2
                0399          DO WHILE ( ABS(rF(k+1)) .LE. zRef(i,j) )
                0400           uRef(i,j) = uRef(i,j) + drF(k) * p5
                0401      &         * ( uVel(i,j,k,bi,bj) + uVel(ip1,j,k,bi,bj) )
                0402           vRef(i,j) = vRef(i,j) + drF(k) * p5
                0403      &         * ( vVel(i,j,k,bi,bj) + vVel(i,jp1,k,bi,bj) )
                0404           k = k+1
                0405          ENDDO
                0406 C     now k should be the index of the layer into which zref reaches
                0407          uRef(i,j) = uRef(i,j) + MAX(0.,zref(i,j)-ABS(rF(k))) * p5
                0408      &        * ( uVel(i,j,k,bi,bj) + uVel(ip1,j,k,bi,bj) )
                0409          vRef(i,j) = vRef(i,j) + MAX(0.,zref(i,j)-ABS(rF(k))) * p5
                0410      &        * ( vVel(i,j,k,bi,bj) + vVel(i,jp1,k,bi,bj) )
                0411          uRef(i,j) = uRef(i,j)/zRef(i,j)
                0412          vRef(i,j) = vRef(i,j)/zRef(i,j)
6af7e0efdf Mart*0413         ENDIF
29e16c9d38 Jean*0414 
6af7e0efdf Mart*0415        ENDDO
                0416       ENDDO
                0417 
                0418       DO k = 1, Nr
                0419        DO j = jmin, jmax
                0420         jm1 = j - 1
                0421         jp1 = j + 1
                0422         DO i = imin, imax
                0423          im1 = i - 1
                0424          ip1 = i + 1
                0425          dVsq(i,j,k) = p5 * (
29e16c9d38 Jean*0426      &        (uRef(i,j) - uVel(i,  j,  k,bi,bj)) *
                0427      &        (uRef(i,j) - uVel(i,  j,  k,bi,bj)) +
                0428      &        (uRef(i,j) - uVel(ip1,j,  k,bi,bj)) *
                0429      &        (uRef(i,j) - uVel(ip1,j,  k,bi,bj)) +
                0430      &        (vRef(i,j) - vVel(i,  j,  k,bi,bj)) *
                0431      &        (vRef(i,j) - vVel(i,  j,  k,bi,bj)) +
                0432      &        (vRef(i,j) - vVel(i,  jp1,k,bi,bj)) *
                0433      &        (vRef(i,j) - vVel(i,  jp1,k,bi,bj)) )
6af7e0efdf Mart*0434 #ifdef KPP_SMOOTH_DVSQ
                0435          dVsq(i,j,k) = p5 * dVsq(i,j,k) + p125 * (
29e16c9d38 Jean*0436      &        (uRef(i,j) - uVel(i,  jm1,k,bi,bj)) *
                0437      &        (uRef(i,j) - uVel(i,  jm1,k,bi,bj)) +
                0438      &        (uRef(i,j) - uVel(ip1,jm1,k,bi,bj)) *
                0439      &        (uRef(i,j) - uVel(ip1,jm1,k,bi,bj)) +
                0440      &        (uRef(i,j) - uVel(i,  jp1,k,bi,bj)) *
                0441      &        (uRef(i,j) - uVel(i,  jp1,k,bi,bj)) +
                0442      &        (uRef(i,j) - uVel(ip1,jp1,k,bi,bj)) *
                0443      &        (uRef(i,j) - uVel(ip1,jp1,k,bi,bj)) +
                0444      &        (vRef(i,j) - vVel(im1,j,  k,bi,bj)) *
                0445      &        (vRef(i,j) - vVel(im1,j,  k,bi,bj)) +
                0446      &        (vRef(i,j) - vVel(im1,jp1,k,bi,bj)) *
                0447      &        (vRef(i,j) - vVel(im1,jp1,k,bi,bj)) +
                0448      &        (vRef(i,j) - vVel(ip1,j,  k,bi,bj)) *
                0449      &        (vRef(i,j) - vVel(ip1,j,  k,bi,bj)) +
                0450      &        (vRef(i,j) - vVel(ip1,jp1,k,bi,bj)) *
                0451      &        (vRef(i,j) - vVel(ip1,jp1,k,bi,bj)) )
6af7e0efdf Mart*0452 #endif /* KPP_SMOOTH_DVSQ */
                0453         ENDDO
                0454        ENDDO
                0455       ENDDO
                0456 
                0457 #else /* KPP_ESTIMATE_UREF */
                0458 
                0459       DO k = 1, Nr
                0460        DO j = jmin, jmax
                0461         jm1 = j - 1
                0462         jp1 = j + 1
                0463         DO i = imin, imax
                0464          im1 = i - 1
                0465          ip1 = i + 1
                0466          dVsq(i,j,k) = p5 * (
29e16c9d38 Jean*0467      &        (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  k,bi,bj)) *
                0468      &        (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  k,bi,bj)) +
                0469      &        (uVel(ip1,j,  1,bi,bj)-uVel(ip1,j,  k,bi,bj)) *
                0470      &        (uVel(ip1,j,  1,bi,bj)-uVel(ip1,j,  k,bi,bj)) +
                0471      &        (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  k,bi,bj)) *
                0472      &        (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  k,bi,bj)) +
                0473      &        (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,k,bi,bj)) *
                0474      &        (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,k,bi,bj)) )
6af7e0efdf Mart*0475 #ifdef KPP_SMOOTH_DVSQ
                0476          dVsq(i,j,k) = p5 * dVsq(i,j,k) + p125 * (
29e16c9d38 Jean*0477      &        (uVel(i,  jm1,1,bi,bj)-uVel(i,  jm1,k,bi,bj)) *
                0478      &        (uVel(i,  jm1,1,bi,bj)-uVel(i,  jm1,k,bi,bj)) +
                0479      &        (uVel(ip1,jm1,1,bi,bj)-uVel(ip1,jm1,k,bi,bj)) *
                0480      &        (uVel(ip1,jm1,1,bi,bj)-uVel(ip1,jm1,k,bi,bj)) +
                0481      &        (uVel(i,  jp1,1,bi,bj)-uVel(i,  jp1,k,bi,bj)) *
                0482      &        (uVel(i,  jp1,1,bi,bj)-uVel(i,  jp1,k,bi,bj)) +
                0483      &        (uVel(ip1,jp1,1,bi,bj)-uVel(ip1,jp1,k,bi,bj)) *
                0484      &        (uVel(ip1,jp1,1,bi,bj)-uVel(ip1,jp1,k,bi,bj)) +
                0485      &        (vVel(im1,j,  1,bi,bj)-vVel(im1,j,  k,bi,bj)) *
                0486      &        (vVel(im1,j,  1,bi,bj)-vVel(im1,j,  k,bi,bj)) +
                0487      &        (vVel(im1,jp1,1,bi,bj)-vVel(im1,jp1,k,bi,bj)) *
                0488      &        (vVel(im1,jp1,1,bi,bj)-vVel(im1,jp1,k,bi,bj)) +
                0489      &        (vVel(ip1,j,  1,bi,bj)-vVel(ip1,j,  k,bi,bj)) *
                0490      &        (vVel(ip1,j,  1,bi,bj)-vVel(ip1,j,  k,bi,bj)) +
                0491      &        (vVel(ip1,jp1,1,bi,bj)-vVel(ip1,jp1,k,bi,bj)) *
                0492      &        (vVel(ip1,jp1,1,bi,bj)-vVel(ip1,jp1,k,bi,bj)) )
6af7e0efdf Mart*0493 #endif /* KPP_SMOOTH_DVSQ */
                0494         ENDDO
                0495        ENDDO
                0496       ENDDO
                0497 
                0498 #endif /* KPP_ESTIMATE_UREF */
                0499 
2b4c90c108 Mart*0500 #ifdef ALLOW_DIAGNOSTICS
                0501       IF ( useDiagnostics ) THEN
                0502        CALL DIAGNOSTICS_FILL(dVsq,  'KPPdVsq ',0,Nr,2,bi,bj,myThid)
                0503       ENDIF
                0504 #endif /* ALLOW_DIAGNOSTICS */
                0505 
6af7e0efdf Mart*0506       RETURN
                0507       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!