pkg/kpp/kpp_forcing_surf.F

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)
00c7090dc0 Mart*0231        ENDDO
                0232       ENDDO
                0233       IF ( selectPenetratingSW .GE. 1 ) THEN
                0234        DO j = jmin, jmax
                0235         DO i = imin, imax
                0236          bosol(i,j) = gravity * TTALPHA(i,j,1) * Qsw(i,j,bi,bj)
                0237      &       *recip_Cp*recip_rhoConst
2b4c90c108 Mart*0238      &       / rhoSurf(i,j)
                0239 CML  &       * work3(i,j)
00c7090dc0 Mart*0240         ENDDO
29e16c9d38 Jean*0241        ENDDO
00c7090dc0 Mart*0242       ENDIF
30c6f5b1cd An T*0243
                0244 #ifdef ALLOW_SALT_PLUME
1f89baba18 Patr*0245 Catn: need check sign of SPforcT
                0246 Cnote: on input, if notdef salt_plume_volume, SPforc[S,T](k>1)=!0
30c6f5b1cd An T*0247       IF ( useSALT_PLUME ) THEN
1f89baba18 Patr*0248 #ifdef SALT_PLUME_VOLUME
a67734cba4 Mart*0249        DO k = 1, Nr
                0250         kp1 = k+1
                0251         DO j = jmin, jmax
                0252          DO i = imin, imax
2b4c90c108 Mart*0253           temparray(i,j) = - gravity *
                0254      &            ( SSBETA(i,j,k) * SPforcS(i,j,k) +
                0255      &              TTALPHA(i,j,k)* SPforcT(i,j,k) / HeatCapacity_Cp )
                0256      &              * recip_rhoConst / rhoSurf(i,j)
                0257 CML  &              * recip_rhoConst * work3(i,j)
                0258           boplume(i,j,kp1) = boplume(i,j,k)+temparray(i,j)
29e16c9d38 Jean*0259          ENDDO
a67734cba4 Mart*0260         ENDDO
                0261        ENDDO
1f89baba18 Patr*0262 #else /* SALT_PLUME_VOLUME */
a67734cba4 Mart*0263        DO j = jmin, jmax
                0264         DO i = imin, imax
                0265 C     initialization has been done already and does not have to be
                0266 C     repeated here (and definitely not inside the i/j loops)
                0267 CML         DO k = 2, Nrp1
2b4c90c108 Mart*0268 CML          boplume(i,j,k) = 0. _d 0
a67734cba4 Mart*0269 CML         ENDDO
2b4c90c108 Mart*0270          boplume(i,j,1) = - gravity * SSBETA(i,j,1)
1f89baba18 Patr*0271      &              * SPforcS(i,j,1)
2b4c90c108 Mart*0272      &              * recip_rhoConst / rhoSurf(i,j)
                0273 CML  &              * recip_rhoConst * work3(i,j)
a67734cba4 Mart*0274         ENDDO
                0275        ENDDO
1f89baba18 Patr*0276 #endif /* SALT_PLUME_VOLUME */
30c6f5b1cd An T*0277       ENDIF
                0278 #endif /* ALLOW_SALT_PLUME */
                0279
                0280 #ifdef ALLOW_DIAGNOSTICS
                0281       IF ( useDiagnostics ) THEN
a67734cba4 Mart*0282        CALL DIAGNOSTICS_FILL(bo     ,'KPPbo   ',0,1,2,bi,bj,myThid)
                0283        CALL DIAGNOSTICS_FILL(bosol  ,'KPPbosol',0,1,2,bi,bj,myThid)
30c6f5b1cd An T*0284 #ifdef ALLOW_SALT_PLUME
a67734cba4 Mart*0285        CALL DIAGNOSTICS_FILL(boplume,'KPPboplm',0,Nr,2,bi,bj,myThid)
30c6f5b1cd An T*0286 #endif /* ALLOW_SALT_PLUME */
                0287       ENDIF
                0288 #endif /* ALLOW_DIAGNOSTICS */
                0289
29e16c9d38 Jean*0290 C------------------------------------------------------------------------
                0291 C     velocity shear
                0292 C     --------------
                0293 C     Get velocity shear squared, averaged from "u,v-grid"
                0294 C     onto "t-grid" (in (m/s)**2):
                0295 C     dVsq(k)=(Uref-U(k))**2+(Vref-V(k))**2      at grid levels
                0296 C------------------------------------------------------------------------
6af7e0efdf Mart*0297
29e16c9d38 Jean*0298 C initialize arrays to zero
6af7e0efdf Mart*0299       DO k = 1, Nr
                0300        DO j = 1-OLy, sNy+OLy
                0301         DO i = 1-OLx, sNx+OLx
                0302          dVsq(i,j,k) = p0
29e16c9d38 Jean*0303         ENDDO
                0304        ENDDO
                0305       ENDDO
6af7e0efdf Mart*0306
29e16c9d38 Jean*0307 C     dVsq computation
6af7e0efdf Mart*0308
                0309 #ifdef KPP_ESTIMATE_UREF
                0310
29e16c9d38 Jean*0311 C     Get rid of vertical resolution dependence of dVsq term by
                0312 C     estimating a surface velocity that is independent of first level
                0313 C     thickness in the model.  First determine mixed layer depth hMix.
                0314 C     Second zRef = espilon * hMix.  Third determine roughness length
                0315 C     scale z0.  Third estimate reference velocity.
a67734cba4 Mart*0316 CML   Comment: so far this code only handles the case of coarse resolution.
                0317 CML   High resolution, i.e. the case zref=hMix*epsilon > drF(1), is not
                0318 CML   handled by this code.
6af7e0efdf Mart*0319
a67734cba4 Mart*0320       zFac = ABS(rF(3)) * LOG ( rF(3) / rF(2) ) * recip_drF(2)
6af7e0efdf Mart*0321       DO j = jmin, jmax
                0322        jp1 = j + 1
                0323        DO i = imin, imax
                0324         ip1 = i + 1
                0325
29e16c9d38 Jean*0326 C     Determine mixed layer depth hMix as the shallowest depth at which
                0327 C     dB/dz exceeds 5.2e-5 s^-2.
a67734cba4 Mart*0328         kTmp = nzmax(i,j,bi,bj)
                0329         DO k = Nr, 1, -1
6af7e0efdf Mart*0330          IF ( k .LT. nzmax(i,j,bi,bj) .AND.
2b4c90c108 Mart*0331      &        maskC(i,j,k,bi,bj) .GT. 0. .AND.
a67734cba4 Mart*0332      &        dbloc(i,j,k) * recip_drC(k+1) .GT. dB_dz )
                0333      &        kTmp = k
6af7e0efdf Mart*0334         ENDDO
                0335
29e16c9d38 Jean*0336 C     Linearly interpolate to find hMix.
a67734cba4 Mart*0337         k = kTmp
6af7e0efdf Mart*0338         IF ( k .EQ. 0 .OR. nzmax(i,j,bi,bj) .EQ. 1 ) THEN
                0339          zRef(i,j) = p0
                0340         ELSEIF ( k .EQ. 1) THEN
a67734cba4 Mart*0341          dBdz2 = dbloc(i,j,1) *recip_drC(2)
6af7e0efdf Mart*0342          zRef(i,j) = drF(1) * dB_dz / dBdz2
                0343         ELSEIF ( k .LT. nzmax(i,j,bi,bj) ) THEN
a67734cba4 Mart*0344          dBdz1 = dbloc(i,j,k-1) * recip_drC(k  )
                0345          dBdz2 = dbloc(i,j,k  ) * recip_drC(k+1)
                0346          zRef(i,j) = ABS(rF(k)) + drF(k) * (dB_dz - dBdz1) /
6af7e0efdf Mart*0347      &        MAX ( phepsi, dBdz2 - dBdz1 )
                0348         ELSE
a67734cba4 Mart*0349          zRef(i,j) = ABS(rF(k+1))
6af7e0efdf Mart*0350         ENDIF
29e16c9d38 Jean*0351
                0352 C     Compute roughness length scale z0 subject to 0 < z0
6af7e0efdf Mart*0353         tempVar1 = p5 * (
                0354      &       (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  2,bi,bj)) *
                0355      &       (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  2,bi,bj)) +
                0356      &       (uVel(ip1,j,  1,bi,bj)-uVel(ip1,j,  2,bi,bj)) *
                0357      &       (uVel(ip1,j,  1,bi,bj)-uVel(ip1,j,  2,bi,bj)) +
                0358      &       (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  2,bi,bj)) *
29e16c9d38 Jean*0359      &       (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  2,bi,bj)) +
6af7e0efdf Mart*0360      &       (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,2,bi,bj)) *
                0361      &       (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,2,bi,bj)) )
                0362         IF ( tempVar1 .lt. (epsln*epsln) ) THEN
                0363          tempVar2 = epsln
                0364         ELSE
                0365          tempVar2 = SQRT ( tempVar1 )
                0366         ENDIF
a67734cba4 Mart*0367 C     ustar > sqrt(0.5*drF(1)*phepsi) by definition, so I do not know
                0368 C     why this extra max function call is required. Too small ustar
                0369 C     will lead to negative z0, which is caught in the next line
                0370         z0(i,j) = drF(1) * ( zFac - tempVar2 * vonK / ustar(i,j) )
                0371 CML     &       ( zFac - tempVar2 * vonK / MAX ( ustar(i,j), phepsi ) )
6af7e0efdf Mart*0372         z0(i,j) = MAX ( z0(i,j), phepsi )
                0373
29e16c9d38 Jean*0374 C     zRef is set to 0.1 * hMix subject to z0 <= zRef <= drF(1)
6af7e0efdf Mart*0375         zRef(i,j) = MAX ( epsilon * zRef(i,j), z0(i,j) )
2b4c90c108 Mart*0376 C     with this statement the high resolution case zRef > drF(1) is ruled out
                0377 C     so we comment it out
                0378 CML        zRef(i,j) = MIN ( zRef(i,j), drF(1) )
29e16c9d38 Jean*0379
                0380 C     Estimate reference velocity uRef and vRef.
6af7e0efdf Mart*0381         uRef(i,j) = p5 * ( uVel(i,j,1,bi,bj) + uVel(ip1,j,1,bi,bj) )
                0382         vRef(i,j) = p5 * ( vVel(i,j,1,bi,bj) + vVel(i,jp1,1,bi,bj) )
                0383         IF ( zRef(i,j) .LT. drF(1) ) THEN
29e16c9d38 Jean*0384          ustarX = ( surfForcU(i,  j,bi,bj) +
6af7e0efdf Mart*0385      &        surfForcU(ip1,j,bi,bj) ) * p5 *recip_drF(1)
                0386          ustarY = ( surfForcV(i,j,  bi,bj) +
                0387      &        surfForcV(i,jp1,bi,bj) ) * p5 *recip_drF(1)
                0388          tempVar1 = ustarX * ustarX + ustarY * ustarY
                0389          if ( tempVar1 .lt. (epsln*epsln) ) then
                0390           tempVar2 = epsln
                0391          else
                0392           tempVar2 = SQRT ( tempVar1 )
                0393          endif
                0394          tempVar2 = ustar(i,j) *
a67734cba4 Mart*0395      &        ( LOG ( zRef(i,j) * recip_drF(1) ) +
                0396      &        z0(i,j) / zRef(i,j) - z0(i,j) * recip_drF(1) ) /
6af7e0efdf Mart*0397      &        vonK / tempVar2
                0398          uRef(i,j) = uRef(i,j) + ustarX * tempVar2
                0399          vRef(i,j) = vRef(i,j) + ustarY * tempVar2
2b4c90c108 Mart*0400         ELSEIF ( zRef(i,j) .GE. drF(1) ) THEN
                0401 C     if zref > drF(1), average over zref
                0402          uRef(i,j) = uRef(i,j)*drF(1)
                0403          vRef(i,j) = vRef(i,j)*drF(1)
                0404          k = 2
                0405          DO WHILE ( ABS(rF(k+1)) .LE. zRef(i,j) )
                0406           uRef(i,j) = uRef(i,j) + drF(k) * p5
                0407      &         * ( uVel(i,j,k,bi,bj) + uVel(ip1,j,k,bi,bj) )
                0408           vRef(i,j) = vRef(i,j) + drF(k) * p5
                0409      &         * ( vVel(i,j,k,bi,bj) + vVel(i,jp1,k,bi,bj) )
                0410           k = k+1
                0411          ENDDO
                0412 C     now k should be the index of the layer into which zref reaches
                0413          uRef(i,j) = uRef(i,j) + MAX(0.,zref(i,j)-ABS(rF(k))) * p5
                0414      &        * ( uVel(i,j,k,bi,bj) + uVel(ip1,j,k,bi,bj) )
                0415          vRef(i,j) = vRef(i,j) + MAX(0.,zref(i,j)-ABS(rF(k))) * p5
                0416      &        * ( vVel(i,j,k,bi,bj) + vVel(i,jp1,k,bi,bj) )
                0417          uRef(i,j) = uRef(i,j)/zRef(i,j)
                0418          vRef(i,j) = vRef(i,j)/zRef(i,j)
6af7e0efdf Mart*0419         ENDIF
29e16c9d38 Jean*0420
6af7e0efdf Mart*0421        ENDDO
                0422       ENDDO
                0423
                0424       DO k = 1, Nr
                0425        DO j = jmin, jmax
                0426         jm1 = j - 1
                0427         jp1 = j + 1
                0428         DO i = imin, imax
                0429          im1 = i - 1
                0430          ip1 = i + 1
                0431          dVsq(i,j,k) = p5 * (
29e16c9d38 Jean*0432      &        (uRef(i,j) - uVel(i,  j,  k,bi,bj)) *
                0433      &        (uRef(i,j) - uVel(i,  j,  k,bi,bj)) +
                0434      &        (uRef(i,j) - uVel(ip1,j,  k,bi,bj)) *
                0435      &        (uRef(i,j) - uVel(ip1,j,  k,bi,bj)) +
                0436      &        (vRef(i,j) - vVel(i,  j,  k,bi,bj)) *
                0437      &        (vRef(i,j) - vVel(i,  j,  k,bi,bj)) +
                0438      &        (vRef(i,j) - vVel(i,  jp1,k,bi,bj)) *
                0439      &        (vRef(i,j) - vVel(i,  jp1,k,bi,bj)) )
6af7e0efdf Mart*0440 #ifdef KPP_SMOOTH_DVSQ
                0441          dVsq(i,j,k) = p5 * dVsq(i,j,k) + p125 * (
29e16c9d38 Jean*0442      &        (uRef(i,j) - uVel(i,  jm1,k,bi,bj)) *
                0443      &        (uRef(i,j) - uVel(i,  jm1,k,bi,bj)) +
                0444      &        (uRef(i,j) - uVel(ip1,jm1,k,bi,bj)) *
                0445      &        (uRef(i,j) - uVel(ip1,jm1,k,bi,bj)) +
                0446      &        (uRef(i,j) - uVel(i,  jp1,k,bi,bj)) *
                0447      &        (uRef(i,j) - uVel(i,  jp1,k,bi,bj)) +
                0448      &        (uRef(i,j) - uVel(ip1,jp1,k,bi,bj)) *
                0449      &        (uRef(i,j) - uVel(ip1,jp1,k,bi,bj)) +
                0450      &        (vRef(i,j) - vVel(im1,j,  k,bi,bj)) *
                0451      &        (vRef(i,j) - vVel(im1,j,  k,bi,bj)) +
                0452      &        (vRef(i,j) - vVel(im1,jp1,k,bi,bj)) *
                0453      &        (vRef(i,j) - vVel(im1,jp1,k,bi,bj)) +
                0454      &        (vRef(i,j) - vVel(ip1,j,  k,bi,bj)) *
                0455      &        (vRef(i,j) - vVel(ip1,j,  k,bi,bj)) +
                0456      &        (vRef(i,j) - vVel(ip1,jp1,k,bi,bj)) *
                0457      &        (vRef(i,j) - vVel(ip1,jp1,k,bi,bj)) )
6af7e0efdf Mart*0458 #endif /* KPP_SMOOTH_DVSQ */
                0459         ENDDO
                0460        ENDDO
                0461       ENDDO
                0462
                0463 #else /* KPP_ESTIMATE_UREF */
                0464
                0465       DO k = 1, Nr
                0466        DO j = jmin, jmax
                0467         jm1 = j - 1
                0468         jp1 = j + 1
                0469         DO i = imin, imax
                0470          im1 = i - 1
                0471          ip1 = i + 1
                0472          dVsq(i,j,k) = p5 * (
29e16c9d38 Jean*0473      &        (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  k,bi,bj)) *
                0474      &        (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  k,bi,bj)) +
                0475      &        (uVel(ip1,j,  1,bi,bj)-uVel(ip1,j,  k,bi,bj)) *
                0476      &        (uVel(ip1,j,  1,bi,bj)-uVel(ip1,j,  k,bi,bj)) +
                0477      &        (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  k,bi,bj)) *
                0478      &        (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  k,bi,bj)) +
                0479      &        (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,k,bi,bj)) *
                0480      &        (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,k,bi,bj)) )
6af7e0efdf Mart*0481 #ifdef KPP_SMOOTH_DVSQ
                0482          dVsq(i,j,k) = p5 * dVsq(i,j,k) + p125 * (
29e16c9d38 Jean*0483      &        (uVel(i,  jm1,1,bi,bj)-uVel(i,  jm1,k,bi,bj)) *
                0484      &        (uVel(i,  jm1,1,bi,bj)-uVel(i,  jm1,k,bi,bj)) +
                0485      &        (uVel(ip1,jm1,1,bi,bj)-uVel(ip1,jm1,k,bi,bj)) *
                0486      &        (uVel(ip1,jm1,1,bi,bj)-uVel(ip1,jm1,k,bi,bj)) +
                0487      &        (uVel(i,  jp1,1,bi,bj)-uVel(i,  jp1,k,bi,bj)) *
                0488      &        (uVel(i,  jp1,1,bi,bj)-uVel(i,  jp1,k,bi,bj)) +
                0489      &        (uVel(ip1,jp1,1,bi,bj)-uVel(ip1,jp1,k,bi,bj)) *
                0490      &        (uVel(ip1,jp1,1,bi,bj)-uVel(ip1,jp1,k,bi,bj)) +
                0491      &        (vVel(im1,j,  1,bi,bj)-vVel(im1,j,  k,bi,bj)) *
                0492      &        (vVel(im1,j,  1,bi,bj)-vVel(im1,j,  k,bi,bj)) +
                0493      &        (vVel(im1,jp1,1,bi,bj)-vVel(im1,jp1,k,bi,bj)) *
                0494      &        (vVel(im1,jp1,1,bi,bj)-vVel(im1,jp1,k,bi,bj)) +
                0495      &        (vVel(ip1,j,  1,bi,bj)-vVel(ip1,j,  k,bi,bj)) *
                0496      &        (vVel(ip1,j,  1,bi,bj)-vVel(ip1,j,  k,bi,bj)) +
                0497      &        (vVel(ip1,jp1,1,bi,bj)-vVel(ip1,jp1,k,bi,bj)) *
                0498      &        (vVel(ip1,jp1,1,bi,bj)-vVel(ip1,jp1,k,bi,bj)) )
6af7e0efdf Mart*0499 #endif /* KPP_SMOOTH_DVSQ */
                0500         ENDDO
                0501        ENDDO
                0502       ENDDO
                0503
                0504 #endif /* KPP_ESTIMATE_UREF */
                0505
2b4c90c108 Mart*0506 #ifdef ALLOW_DIAGNOSTICS
                0507       IF ( useDiagnostics ) THEN
                0508        CALL DIAGNOSTICS_FILL(dVsq,  'KPPdVsq ',0,Nr,2,bi,bj,myThid)
                0509       ENDIF
                0510 #endif /* ALLOW_DIAGNOSTICS */
                0511
6af7e0efdf Mart*0512       RETURN
                0513       END