pkg/kpp/kpp_calc.F

2fa42a6013 Alis*0001 #include "KPP_OPTIONS.h"
853c5e0e2c Jean*0002 #ifdef ALLOW_AUTODIFF
                0003 # include "AUTODIFF_OPTIONS.h"
                0004 #endif
1f89baba18 Patr*0005 #ifdef ALLOW_SALT_PLUME
853c5e0e2c Jean*0006 # include "SALT_PLUME_OPTIONS.h"
1f89baba18 Patr*0007 #endif
2fa42a6013 Alis*0008
2b4c90c108 Mart*0009 C-- File kpp_calc.F:
                0010 C--  Contents
                0011 C--  o KPP_CALC         - Main KPP interface routine.
b7b61e618a Mart*0012 C--  o KPP_CALC_DUMMY   - Dummy routine for TAF
2b4c90c108 Mart*0013
                0014 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
ef080e1d37 Dimi*0015 CBOP
                0016 C !ROUTINE: KPP_CALC
                0017
                0018 C !INTERFACE: ==========================================================
c417c38d79 Jean*0019       SUBROUTINE KPP_CALC(
                0020      I     bi, bj, myTime, myIter, myThid )
ef080e1d37 Dimi*0021
                0022 C !DESCRIPTION: \bv
c417c38d79 Jean*0023 C     *==========================================================*
2fa42a6013 Alis*0024 C     | SUBROUTINE KPP_CALC                                      |
                0025 C     | o Compute all KPP fields defined in KPP.h                |
c417c38d79 Jean*0026 C     *==========================================================*
2fa42a6013 Alis*0027 C     | This subroutine serves as an interface between MITGCMUV  |
                0028 C     | code and NCOM 1-D routines in kpp_routines.F             |
c417c38d79 Jean*0029 C     *==========================================================*
2fa42a6013 Alis*0030
                0031 c=======================================================================
                0032 c
                0033 c     written  by  : jan morzel, august  11, 1994
                0034 c     modified by  : jan morzel, january 25, 1995 : "dVsq" and 1d code
                0035 c                    detlef stammer, august, 1997 : for MIT GCM Classic
                0036 c                    d. menemenlis,    july, 1998 : for MIT GCM UV
                0037 c
                0038 c     compute vertical mixing coefficients based on the k-profile
                0039 c     and oceanic planetary boundary layer scheme by large & mcwilliams.
                0040 c
                0041 c     summary:
                0042 c     - compute interior mixing everywhere:
                0043 c       interior mixing gets computed at all interfaces due to constant
                0044 c       internal wave background activity ("fkpm" and "fkph"), which
                0045 c       is enhanced in places of static instability (local richardson
                0046 c       number < 0).
                0047 c       Additionally, mixing can be enhanced by adding contribution due
                0048 c       to shear instability which is a function of the local richardson
a9d2e4c565 Jean*0049 c       number
2fa42a6013 Alis*0050 c     - double diffusivity:
                0051 c       interior mixing can be enhanced by double diffusion due to salt
                0052 c       fingering and diffusive convection (ifdef "kmixdd").
                0053 c     - kpp scheme in the boundary layer:
a9d2e4c565 Jean*0054 c
2fa42a6013 Alis*0055 c       a.boundary layer depth:
a9d2e4c565 Jean*0056 c         at every gridpoint the depth of the oceanic boundary layer
2fa42a6013 Alis*0057 c         ("hbl") gets computed by evaluating bulk richardson numbers.
                0058 c       b.boundary layer mixing:
a9d2e4c565 Jean*0059 c         within the boundary layer, above hbl, vertical mixing is
2fa42a6013 Alis*0060 c         determined by turbulent surface fluxes, and interior mixing at
                0061 c         the lower boundary, i.e. at hbl.
a9d2e4c565 Jean*0062 c
c417c38d79 Jean*0063 c     this subroutine provides the interface between the MITGCM and
a9d2e4c565 Jean*0064 c     the routine "kppmix", where boundary layer depth, vertical
2fa42a6013 Alis*0065 c     viscosity, vertical diffusivity, and counter gradient term (ghat)
                0066 c     are computed slabwise.
                0067 c     note: subroutine "kppmix" uses m-k-s units.
                0068 c
                0069 c     time level:
a9d2e4c565 Jean*0070 c     input tracer and velocity profiles are evaluated at time level
2fa42a6013 Alis*0071 c     tau, surface fluxes come from tau or tau-1.
                0072 c
                0073 c     grid option:
                0074 c     in this "1-grid" implementation, diffusivity and viscosity
                0075 c     profiles are computed on the "t-grid" (by using velocity shear
                0076 c     profiles averaged from the "u,v-grid" onto the "t-grid"; note, that
a9d2e4c565 Jean*0077 c     the averaging includes zero values on coastal and seafloor grid
                0078 c     points).  viscosity on the "u,v-grid" is computed by averaging the
2fa42a6013 Alis*0079 c     "t-grid" viscosity values onto the "u,v-grid".
                0080 c
                0081 c     vertical grid:
a9d2e4c565 Jean*0082 c     mixing coefficients get evaluated at the bottom of the lowest
                0083 c     layer, i.e., at depth zw(Nr).  these values are only useful when
2fa42a6013 Alis*0084 c     the model ocean domain does not include the entire ocean down to
                0085 c     the seafloor ("upperocean" setup) and allows flux through the
a9d2e4c565 Jean*0086 c     bottom of the domain.  for full-depth runs, these mixing
2fa42a6013 Alis*0087 c     coefficients are being zeroed out before leaving this subroutine.
                0088 c
                0089 c-------------------------------------------------------------------------
                0090
                0091 c global parameters updated by kpp_calc
                0092 c     KPPviscAz   - KPP eddy viscosity coefficient                 (m^2/s)
                0093 c     KPPdiffKzT  - KPP diffusion coefficient for temperature      (m^2/s)
                0094 c     KPPdiffKzS  - KPP diffusion coefficient for salt and tracers (m^2/s)
                0095 c     KPPghat     - Nonlocal transport coefficient                 (s/m^2)
                0096 c     KPPhbl      - Boundary layer depth on "t-grid"                   (m)
00c7090dc0 Mart*0097 c     KPPfrac     - Fraction of short-wave flux heating the mixing layer
63ceaaa79c Dimi*0098 c     KPPplumefrac- Fraction of saltplume (flux) penetrating mixing layer
2fa42a6013 Alis*0099
                0100 c--   KPP_CALC computes vertical viscosity and diffusivity for region
                0101 c     (-2:sNx+3,-2:sNy+3) as required by CALC_DIFFUSIVITY and requires
c572a3ecb8 Jean*0102 c     values of uVel, vVel, surfaceForcingU, surfaceForcingV in the
1d478690dc Patr*0103 c     region (-2:sNx+4,-2:sNy+4).
2fa42a6013 Alis*0104 c     Hence overlap region needs to be set OLx=4, OLy=4.
ef080e1d37 Dimi*0105 c \ev
2fa42a6013 Alis*0106
ef080e1d37 Dimi*0107 C !USES: ===============================================================
2b4c90c108 Mart*0108       IMPLICIT NONE
2fa42a6013 Alis*0109 #include "SIZE.h"
                0110 #include "EEPARAMS.h"
                0111 #include "PARAMS.h"
                0112 #include "DYNVARS.h"
                0113 #include "KPP.h"
                0114 #include "KPP_PARAMS.h"
                0115 #include "FFIELDS.h"
                0116 #include "GRID.h"
059d9fc14f Dimi*0117 #include "GAD.h"
8749f0b4bc Dimi*0118 #ifdef ALLOW_SALT_PLUME
                0119 # include "SALT_PLUME.h"
                0120 #endif /* ALLOW_SALT_PLUME */
36aa272a96 Mart*0121 #ifdef ALLOW_SHELFICE
                0122 # include "SHELFICE.h"
                0123 #endif /* ALLOW_SHELFICE */
2fa42a6013 Alis*0124 #ifdef ALLOW_AUTODIFF_TAMC
8749f0b4bc Dimi*0125 # include "tamc.h"
2fa42a6013 Alis*0126 #endif /* ALLOW_AUTODIFF_TAMC */
                0127
ef080e1d37 Dimi*0128 C !INPUT PARAMETERS: ===================================================
2fa42a6013 Alis*0129 c Routine arguments
c417c38d79 Jean*0130 c     bi, bj :: Current tile indices
                0131 c     myTime :: Current time in simulation
                0132 c     myIter :: Current iteration number in simulation
                0133 c     myThid :: My Thread Id. number
2fa42a6013 Alis*0134       INTEGER bi, bj
                0135       _RL     myTime
c417c38d79 Jean*0136       INTEGER myIter
                0137       INTEGER myThid
2fa42a6013 Alis*0138
                0139 #ifdef ALLOW_KPP
2b4c90c108 Mart*0140 C !FUNCTIONS:
                0141       LOGICAL  DIFFERENT_MULTIPLE
                0142       EXTERNAL DIFFERENT_MULTIPLE
2fa42a6013 Alis*0143
ef080e1d37 Dimi*0144 C !LOCAL VARIABLES: ====================================================
956c0a5824 Patr*0145 c Local constants
                0146 c     minusone, p0, p5, p25, p125, p0625
7c50f07931 Mart*0147 c     iMin, iMax, jMin, jMax  - array computation indices
956c0a5824 Patr*0148       _RL        minusone
2b4c90c108 Mart*0149       PARAMETER( minusone=-1.0)
fe66051ebd Dimi*0150       _RL        p0    , p5    , p25     , p125      , p0625
2b4c90c108 Mart*0151       PARAMETER( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 )
7c50f07931 Mart*0152       INTEGER   iMin      ,iMax          ,jMin      ,jMax
                0153       PARAMETER(iMin=2-OLx,iMax=sNx+OLx-1,jMin=2-OLy,jMax=sNy+OLy-1)
956c0a5824 Patr*0154
2fa42a6013 Alis*0155 c Local arrays and variables
00c7090dc0 Mart*0156 c     sdens  (nx,ny        - density of surface layer                (kg/m^3)
2fa42a6013 Alis*0157 c     ustar  (nx,ny)       - surface friction velocity                  (m/s)
                0158 c     bo     (nx,ny)       - surface turbulent buoyancy forcing     (m^2/s^3)
                0159 c     bosol  (nx,ny)       - surface radiative buoyancy forcing     (m^2/s^3)
                0160 c     shsq   (nx,ny,Nr)    - local velocity shear squared
                0161 c                            at interfaces for ri_iwmix             (m^2/s^2)
                0162 c     dVsq   (nx,ny,Nr)    - velocity shear re surface squared
                0163 c                            at grid levels for bldepth             (m^2/s^2)
                0164 c     dbloc  (nx,ny,Nr)    - local delta buoyancy at interfaces
                0165 c                            for ri_iwmix and bldepth                 (m/s^2)
                0166 c     Ritop  (nx,ny,Nr)    - numerator of bulk richardson number
                0167 c                            at grid levels for bldepth
                0168 c     vddiff (nx,ny,Nrp2,1)- vertical viscosity on "t-grid"           (m^2/s)
42e0e34daa Dimi*0169 c     vddiff (nx,ny,Nrp2,2)- vert. diff. on next row for salt&tracers (m^2/s)
                0170 c     vddiff (nx,ny,Nrp2,3)- vert. diff. on next row for temperature  (m^2/s)
2fa42a6013 Alis*0171 c     ghat   (nx,ny,Nr)    - nonlocal transport coefficient           (s/m^2)
                0172 c     hbl    (nx,ny)       - mixing layer depth                           (m)
00c7090dc0 Mart*0173 c     kbl    (nx,ny)       - k-cell index of mixing layer depth
2fa42a6013 Alis*0174 c     kmtj   (nx,ny)       - maximum number of wet levels in each column
00c7090dc0 Mart*0175 c     work?  (nx,ny)       - horizontal working arrays
                0176 c     temp?  (nx,ny,Nr)    - 3d working arrays
                0177 c     boplume(nx,ny,Nrp1)  - surface haline buoyancy forcing        (m^2/s^3)
                0178       _RL sdens ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
fe66051ebd Dimi*0179       _RL ustar ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
                0180       _RL bo    ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
                0181       _RL bosol ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
                0182       _RL shsq  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
                0183       _RL dVsq  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
                0184       _RL dbloc ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
                0185       _RL Ritop ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
                0186       _RL vddiff( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, 0:Nrp1, mdiff )
                0187       _RL ghat  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
                0188       _RL hbl   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
277d95ff7e Patr*0189 cph(
fe66051ebd Dimi*0190       _RL TTALPHA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
                0191       _RL SSBETA ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
277d95ff7e Patr*0192 cph)
00c7090dc0 Mart*0193 #if ( defined ALLOW_SALT_PLUME || defined SHORTWAVE_HEATING )
                0194       _RL worka ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
                0195 #endif
                0196 #ifdef ALLOW_SALT_PLUME
                0197       _RL temp1   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr          )
                0198       _RL temp2   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr          )
                0199       _RL boplume ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1        )
                0200 #ifdef SALT_PLUME_SPLIT_BASIN
                0201       _RL lon ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                  )
                0202       _RL lat ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                  )
                0203 #endif /* SALT_PLUME_SPLIT_BASIN */
                0204 #endif /* ALLOW_SALT_PLUME */
                0205 #ifdef SHORTWAVE_HEATING
                0206       INTEGER kbl(1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
                0207 C     interpolation factor 1 >= (KPPhbl+rF(kbl))/drF(kbl) >= 0
                0208       _RL rFac
                0209 #endif
                0210       INTEGER i, j, kSurf, k, kp1, km1, im1, ip1, jm1, jp1
edb6656069 Mart*0211 C     ikey :: tape key (depends on tiles)
00c7090dc0 Mart*0212       INTEGER ikey
853c5e0e2c Jean*0213 CEOP
2fa42a6013 Alis*0214
ab9b91ea7a Patr*0215 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0216       ikey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
853c5e0e2c Jean*0217 #else /* ALLOW_AUTODIFF_TAMC */
edb6656069 Mart*0218       ikey = 0
ab9b91ea7a Patr*0219 #endif /* ALLOW_AUTODIFF_TAMC */
                0220
2fa42a6013 Alis*0221 c     Check to see if new vertical mixing coefficient should be computed now?
94a46dfe0d Jean*0222       IF ( DIFFERENT_MULTIPLE(kpp_freq,myTime,deltaTClock)
adf557b193 Jean*0223      1     .OR. myTime .EQ. startTime ) THEN
c417c38d79 Jean*0224
2fa42a6013 Alis*0225 c-----------------------------------------------------------------------
                0226 c     prepare input arrays for subroutine "kppmix" to compute
                0227 c     viscosity and diffusivity and ghat.
                0228 c     All input arrays need to be in m-k-s units.
                0229 c
                0230 c     note: for the computation of the bulk richardson number in the
                0231 c     "bldepth" subroutine, gradients of velocity and buoyancy are
                0232 c     required at every depth. in the case of very fine vertical grids
                0233 c     (thickness of top layer < 2m), the surface reference depth must
                0234 c     be set to zref=epsilon/2*zgrid(k), and the reference value
                0235 c     of velocity and buoyancy must be computed as vertical average
                0236 c     between the surface and 2*zref.  in the case of coarse vertical
                0237 c     grids zref is zgrid(1)/2., and the surface reference value is
                0238 c     simply the surface value at zgrid(1).
                0239 c-----------------------------------------------------------------------
                0240
                0241 c------------------------------------------------------------------------
                0242 c     density related quantities
                0243 c     --------------------------
                0244 c
00c7090dc0 Mart*0245 c      sdens   - density of surface layer                        (kg/m^3)
2fa42a6013 Alis*0246 c      dbloc   - local buoyancy gradient at Nr interfaces
                0247 c                g/rho{k+1,k+1} * [ drho{k,k+1}-drho{k+1,k+1} ]   (m/s^2)
                0248 c      dbsfc (stored in Ritop to conserve stack memory)
                0249 c              - buoyancy difference with respect to the surface
                0250 c                g * [ drho{1,k}/rho{1,k} - drho{k,k}/rho{k,k} ]  (m/s^2)
00c7090dc0 Mart*0251 c      ttalpha - thermal expansion coefficient without 1/rho factor
2fa42a6013 Alis*0252 c                d(rho{k,k})/d(T(k))                           (kg/m^3/C)
00c7090dc0 Mart*0253 c      ssbeta  - salt expansion coefficient without 1/rho factor
ba0b047096 Mart*0254 c                d(rho{k,k})/d(S(k))                    ((kg/m^3)/(g/kg))
2fa42a6013 Alis*0255 c------------------------------------------------------------------------
                0256
                0257       CALL STATEKPP(
00c7090dc0 Mart*0258      O     sdens, dbloc, Ritop,
25c3477f99 Mart*0259      O     TTALPHA, SSBETA,
edb6656069 Mart*0260      I     ikey, bi, bj, myThid )
2fa42a6013 Alis*0261
                0262       DO k = 1, Nr
a67734cba4 Mart*0263        DO j = 1-OLy, sNy+OLy
                0264         DO i = 1-OLx, sNx+OLx
                0265          ghat(i,j,k) = dbloc(i,j,k)
                0266         ENDDO
                0267        ENDDO
2fa42a6013 Alis*0268       ENDDO
                0269
956c0a5824 Patr*0270 #ifdef KPP_SMOOTH_DBLOC
                0271 c     horizontally smooth dbloc with a 121 filter
                0272 c     smooth dbloc stored in ghat to save space
                0273 c     dbloc(k) is buoyancy gradientnote between k and k+1
                0274 c     levels therefore k+1 mask must be used
                0275
                0276       DO k = 1, Nr-1
a67734cba4 Mart*0277        CALL SMOOTH_HORIZ (
                0278      I      k+1, bi, bj,
                0279      U      ghat (1-OLx,1-OLy,k),
                0280      I      myThid )
956c0a5824 Patr*0281       ENDDO
                0282
2fa42a6013 Alis*0283 #endif /* KPP_SMOOTH_DBLOC */
                0284
                0285 #ifdef KPP_SMOOTH_DENS
                0286 c     horizontally smooth density related quantities with 121 filters
fe66051ebd Dimi*0287       CALL SMOOTH_HORIZ (
956c0a5824 Patr*0288      I     1, bi, bj,
00c7090dc0 Mart*0289      U     sdens,
25c3477f99 Mart*0290      I     myThid )
2fa42a6013 Alis*0291       DO k = 1, Nr
fe66051ebd Dimi*0292          CALL SMOOTH_HORIZ (
956c0a5824 Patr*0293      I        k+1, bi, bj,
c417c38d79 Jean*0294      U        dbloc (1-OLx,1-OLy,k),
25c3477f99 Mart*0295      I        myThid )
fe66051ebd Dimi*0296          CALL SMOOTH_HORIZ (
2fa42a6013 Alis*0297      I        k, bi, bj,
c417c38d79 Jean*0298      U        Ritop (1-OLx,1-OLy,k),
25c3477f99 Mart*0299      I        myThid )
fe66051ebd Dimi*0300          CALL SMOOTH_HORIZ (
2fa42a6013 Alis*0301      I        k, bi, bj,
c417c38d79 Jean*0302      U        TTALPHA(1-OLx,1-OLy,k),
25c3477f99 Mart*0303      I        myThid )
fe66051ebd Dimi*0304          CALL SMOOTH_HORIZ (
2fa42a6013 Alis*0305      I        k, bi, bj,
c417c38d79 Jean*0306      U        SSBETA(1-OLx,1-OLy,k),
25c3477f99 Mart*0307      I        myThid )
2fa42a6013 Alis*0308       ENDDO
                0309 #endif /* KPP_SMOOTH_DENS */
                0310
7c50f07931 Mart*0311 c     surface index
                0312       kSurf = 1
2fa42a6013 Alis*0313       DO k = 1, Nr
a67734cba4 Mart*0314        km1 = max(1,k-1)
2b4c90c108 Mart*0315        kp1 = min(Nr,k+1)
a67734cba4 Mart*0316        DO j = 1-OLy, sNy+OLy
                0317         DO i = 1-OLx, sNx+OLx
2fa42a6013 Alis*0318
                0319 c     zero out dbloc over land points (so that the convective
2b4c90c108 Mart*0320 c     part of the interior mixing can be diagnosed).
                0321 c     ML: Note that dbloc and ghat (=smoothed dbloc) are defined at
                0322 c     the bottom of the cell rather than the (MITgcm-standard) top
                0323 c     of the cell; therefore the proper masking is maskC(k)*maskC(k+1)
a67734cba4 Mart*0324          dbloc(i,j,k) = dbloc(i,j,k) * maskC(i,j,k,bi,bj)
2b4c90c108 Mart*0325      &        * maskC(i,j,kp1,bi,bj)
a67734cba4 Mart*0326          ghat(i,j,k)  = ghat(i,j,k)  * maskC(i,j,k,bi,bj)
2b4c90c108 Mart*0327      &        * maskC(i,j,kp1,bi,bj)
7c50f07931 Mart*0328 #ifdef ALLOW_SHELFICE
                0329          IF ( useShelfIce ) kSurf = MAX(1,kTopC(i,j,bi,bj))
                0330 #endif /* ALLOW_SHELFICE */
a67734cba4 Mart*0331          Ritop(i,j,k) = Ritop(i,j,k) * maskC(i,j,k,bi,bj)
2b4c90c108 Mart*0332      &        * maskC(i,j,kSurf,bi,bj)
                0333          IF(k.EQ.nzmax(i,j,bi,bj)) THEN
a67734cba4 Mart*0334           dbloc(i,j,k) = p0
                0335           ghat(i,j,k)  = p0
                0336           Ritop(i,j,k) = p0
2b4c90c108 Mart*0337          ENDIF
2fa42a6013 Alis*0338
                0339 c     numerator of bulk richardson number on grid levels
                0340 c     note: land and ocean bottom values need to be set to zero
                0341 c     so that the subroutine "bldepth" works correctly
a67734cba4 Mart*0342          Ritop(i,j,k) = (zgrid(1)-zgrid(k)) * Ritop(i,j,k)
2fa42a6013 Alis*0343
a67734cba4 Mart*0344         ENDDO
                0345        ENDDO
c417c38d79 Jean*0346       ENDDO
2fa42a6013 Alis*0347
32e2940ee8 Patr*0348 cph(
                0349 cph  this avoids a single or double recomp./call of statekpp
853c5e0e2c Jean*0350 #ifdef ALLOW_AUTODIFF_TAMC
00c7090dc0 Mart*0351 CADJ STORE sdens, Ritop       = comlev1_kpp, key = ikey
edb6656069 Mart*0352 CADJ STORE TTALPHA, SSBETA    = comlev1_kpp, key = ikey
853c5e0e2c Jean*0353 #endif /* ALLOW_AUTODIFF_TAMC */
32e2940ee8 Patr*0354 cph)
                0355
36aa272a96 Mart*0356 CML#ifdef ALLOW_SHELFICE
a9d2e4c565 Jean*0357 CMLC     For the pbl parameterisation to work underneath the ice shelves
36aa272a96 Mart*0358 CMLC     it needs to know the surface (ice-ocean) fluxes. However, masking
a9d2e4c565 Jean*0359 CMLC     and indexing problems make this part of the code not work
36aa272a96 Mart*0360 CMLC     underneath the ice shelves and the following lines are only here
                0361 CMLC     to remind me that this still needs to be sorted out.
6af7e0efdf Mart*0362 CML      shelfIceFac = 0. _d 0
                0363 CML      IF ( useShelfIce ) selfIceFac = 1. _d 0
7c50f07931 Mart*0364 CML      DO j = jMin, jMax
                0365 CML       DO i = iMin, iMax
a9d2e4c565 Jean*0366 CML        surfForcT = surfaceForcingT(i,j,bi,bj)
6af7e0efdf Mart*0367 CML     &       + shelficeForcingT(i,j,bi,bj) * shelfIceFac
a9d2e4c565 Jean*0368 CML        surfForcS = surfaceForcingS(i,j,bi,bj)
6af7e0efdf Mart*0369 CML     &       + shelficeForcingS(i,j,bi,bj) * shelfIceFac
c417c38d79 Jean*0370 CML       ENDDO
                0371 CML      ENDDO
36aa272a96 Mart*0372 CML#endif /* ALLOW_SHELFICE */
32e2940ee8 Patr*0373
2fa42a6013 Alis*0374 c------------------------------------------------------------------------
6af7e0efdf Mart*0375 c     friction velocity, turbulent and radiative surface buoyancy forcing
                0376 c     -------------------------------------------------------------------
                0377 c     taux / rho = surfaceForcingU                               (N/m^2)
                0378 c     tauy / rho = surfaceForcingV                               (N/m^2)
                0379 c     ustar = sqrt( sqrt( taux^2 + tauy^2 ) / rho )                (m/s)
                0380 c     bo    = - g * ( alpha*surfaceForcingT +
                0381 c                     beta *surfaceForcingS ) / rho            (m^2/s^3)
                0382 c     bosol = - g * alpha * Qsw * drF(1) / rho                 (m^2/s^3)
63ceaaa79c Dimi*0383 c     boplume = g * (beta * saltPlumeFlux/rhoConst ) /rho      (m^2/s^3)
1f89baba18 Patr*0384 c             = g * (beta * SPforcingS   /rhoConst ) /rho
                0385 c              +g * (alpha* SPforcingT   / ??
6af7e0efdf Mart*0386 c------------------------------------------------------------------------
2fa42a6013 Alis*0387 c     velocity shear
                0388 c     --------------
                0389 c     Get velocity shear squared, averaged from "u,v-grid"
                0390 c     onto "t-grid" (in (m/s)**2):
                0391 c     dVsq(k)=(Uref-U(k))**2+(Vref-V(k))**2      at grid levels
                0392 c     shsq(k)=(U(k)-U(k+1))**2+(V(k)-V(k+1))**2  at interfaces
a9d2e4c565 Jean*0393 c
6af7e0efdf Mart*0394 c     note: Vref can depend on the surface fluxes that is why we compute
                0395 c     dVsq in the subroutine that does the surface related stuff
                0396 c     (admittedly this is a bit messy)
2fa42a6013 Alis*0397 c------------------------------------------------------------------------
c417c38d79 Jean*0398
1f89baba18 Patr*0399 #ifdef ALLOW_SALT_PLUME
7c50f07931 Mart*0400       DO j=jMin,jMax
                0401        DO i=iMin,iMax
1f89baba18 Patr*0402 #ifndef SALT_PLUME_VOLUME
                0403         temp1(i,j,1) = saltPlumeFlux(i,j,bi,bj)
                0404         temp2(i,j,1) = 0. _d 0
                0405         DO k=2,Nr
                0406          temp1(i,j,k) = 0. _d 0
                0407          temp2(i,j,k) = 0. _d 0
                0408         ENDDO
                0409 #else /* def SALT_PLUME_VOLUME */
                0410         DO k=1,Nr
                0411          temp1(i,j,k) = SPforcingS(i,j,k,bi,bj)
                0412          temp2(i,j,k) = SPforcingT(i,j,k,bi,bj)
                0413         ENDDO
                0414 #endif /* SALT_PLUME_VOLUME */
                0415        ENDDO
                0416       ENDDO
                0417 #endif /* ALLOW_SALT_PLUME */
                0418
6af7e0efdf Mart*0419       CALL KPP_FORCING_SURF(
00c7090dc0 Mart*0420      I     sdens, surfaceForcingU, surfaceForcingV,
16fe9e8e73 Jean*0421      I     surfaceForcingT, surfaceForcingS, adjustColdSST_diag,
a9d2e4c565 Jean*0422      I     Qsw,
a67734cba4 Mart*0423 #ifdef KPP_ESTIMATE_UREF
                0424      I     dbloc,
                0425 #endif /* KPP_ESTIMATE_UREF */
63ceaaa79c Dimi*0426 #ifdef ALLOW_SALT_PLUME
1f89baba18 Patr*0427      I     temp1, temp2,
63ceaaa79c Dimi*0428 #endif /* ALLOW_SALT_PLUME */
                0429      I     ttalpha, ssbeta,
a9d2e4c565 Jean*0430      O     ustar, bo, bosol,
63ceaaa79c Dimi*0431 #ifdef ALLOW_SALT_PLUME
                0432      O     boplume,
                0433 #endif /* ALLOW_SALT_PLUME */
                0434      O     dVsq,
edb6656069 Mart*0435      I     ikey, iMin, iMax, jMin, jMax, bi, bj, myTime, myThid )
2fa42a6013 Alis*0436
00c7090dc0 Mart*0437 #ifdef ALLOW_AUTODIFF_TAMC
                0438 CADJ STORE ustar, bo, bosol, dVsq = comlev1_kpp, key = ikey
                0439 # ifdef ALLOW_SALT_PLUME
                0440 CADJ STORE boplume = comlev1_kpp, key = ikey
                0441 # endif
                0442 #endif
2fa42a6013 Alis*0443
6af7e0efdf Mart*0444 c initialize arrays to zero
2fa42a6013 Alis*0445       DO k = 1, Nr
6af7e0efdf Mart*0446        DO j = 1-OLy, sNy+OLy
                0447         DO i = 1-OLx, sNx+OLx
                0448          shsq(i,j,k) = p0
c417c38d79 Jean*0449         ENDDO
                0450        ENDDO
                0451       ENDDO
2fa42a6013 Alis*0452
                0453 c     shsq computation
                0454       DO k = 1, Nrm1
6af7e0efdf Mart*0455        kp1 = k + 1
7c50f07931 Mart*0456        DO j = jMin, jMax
6af7e0efdf Mart*0457         jm1 = j - 1
                0458         jp1 = j + 1
7c50f07931 Mart*0459         DO i = iMin, iMax
6af7e0efdf Mart*0460          im1 = i - 1
                0461          ip1 = i + 1
                0462          shsq(i,j,k) = p5 * (
c417c38d79 Jean*0463      &        (uVel(i,  j,  k,bi,bj)-uVel(i,  j,  kp1,bi,bj)) *
                0464      &        (uVel(i,  j,  k,bi,bj)-uVel(i,  j,  kp1,bi,bj)) +
                0465      &        (uVel(ip1,j,  k,bi,bj)-uVel(ip1,j,  kp1,bi,bj)) *
                0466      &        (uVel(ip1,j,  k,bi,bj)-uVel(ip1,j,  kp1,bi,bj)) +
                0467      &        (vVel(i,  j,  k,bi,bj)-vVel(i,  j,  kp1,bi,bj)) *
                0468      &        (vVel(i,  j,  k,bi,bj)-vVel(i,  j,  kp1,bi,bj)) +
                0469      &        (vVel(i,  jp1,k,bi,bj)-vVel(i,  jp1,kp1,bi,bj)) *
                0470      &        (vVel(i,  jp1,k,bi,bj)-vVel(i,  jp1,kp1,bi,bj)) )
2fa42a6013 Alis*0471 #ifdef KPP_SMOOTH_SHSQ
6af7e0efdf Mart*0472          shsq(i,j,k) = p5 * shsq(i,j,k) + p125 * (
c417c38d79 Jean*0473      &        (uVel(i,  jm1,k,bi,bj)-uVel(i,  jm1,kp1,bi,bj)) *
                0474      &        (uVel(i,  jm1,k,bi,bj)-uVel(i,  jm1,kp1,bi,bj)) +
                0475      &        (uVel(ip1,jm1,k,bi,bj)-uVel(ip1,jm1,kp1,bi,bj)) *
                0476      &        (uVel(ip1,jm1,k,bi,bj)-uVel(ip1,jm1,kp1,bi,bj)) +
                0477      &        (uVel(i,  jp1,k,bi,bj)-uVel(i,  jp1,kp1,bi,bj)) *
                0478      &        (uVel(i,  jp1,k,bi,bj)-uVel(i,  jp1,kp1,bi,bj)) +
                0479      &        (uVel(ip1,jp1,k,bi,bj)-uVel(ip1,jp1,kp1,bi,bj)) *
                0480      &        (uVel(ip1,jp1,k,bi,bj)-uVel(ip1,jp1,kp1,bi,bj)) +
                0481      &        (vVel(im1,j,  k,bi,bj)-vVel(im1,j,  kp1,bi,bj)) *
                0482      &        (vVel(im1,j,  k,bi,bj)-vVel(im1,j,  kp1,bi,bj)) +
                0483      &        (vVel(im1,jp1,k,bi,bj)-vVel(im1,jp1,kp1,bi,bj)) *
                0484      &        (vVel(im1,jp1,k,bi,bj)-vVel(im1,jp1,kp1,bi,bj)) +
                0485      &        (vVel(ip1,j,  k,bi,bj)-vVel(ip1,j,  kp1,bi,bj)) *
                0486      &        (vVel(ip1,j,  k,bi,bj)-vVel(ip1,j,  kp1,bi,bj)) +
                0487      &        (vVel(ip1,jp1,k,bi,bj)-vVel(ip1,jp1,kp1,bi,bj)) *
                0488      &        (vVel(ip1,jp1,k,bi,bj)-vVel(ip1,jp1,kp1,bi,bj)) )
2fa42a6013 Alis*0489 #endif
c417c38d79 Jean*0490         ENDDO
                0491        ENDDO
                0492       ENDDO
2fa42a6013 Alis*0493
2b4c90c108 Mart*0494 #ifdef ALLOW_DIAGNOSTICS
                0495       IF ( useDiagnostics ) THEN
                0496        CALL DIAGNOSTICS_FILL(shsq,  'KPPshsq ',0,Nr,2,bi,bj,myThid)
                0497       ENDIF
                0498 #endif /* ALLOW_DIAGNOSTICS */
32e2940ee8 Patr*0499 cph(
853c5e0e2c Jean*0500 #ifdef ALLOW_AUTODIFF_TAMC
3067745c9b Patr*0501 #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
00c7090dc0 Mart*0502 CADJ STORE shsq = comlev1_kpp, key = ikey
32e2940ee8 Patr*0503 #endif
853c5e0e2c Jean*0504 #endif /* ALLOW_AUTODIFF_TAMC */
32e2940ee8 Patr*0505 cph)
                0506
2fa42a6013 Alis*0507 c-----------------------------------------------------------------------
                0508 c     solve for viscosity, diffusivity, ghat, and hbl on "t-grid"
                0509 c-----------------------------------------------------------------------
                0510
af61e5eb16 Mart*0511 #if ( !defined EXCLUDE_PCELL_MIX_CODE && defined  ALLOW_AUTODIFF_TAMC )
                0512 CADJ STORE KPPdiffKzS(:,:,:,bi,bj) = comlev1_kpp, key = ikey
                0513 CADJ STORE KPPdiffKzT(:,:,:,bi,bj) = comlev1_kpp, key = ikey
                0514 #endif
059d9fc14f Dimi*0515 c     precompute background vertical diffusivities, which are needed for
                0516 c     matching diffusivities at bottom of KPP PBL
                0517       CALL CALC_3D_DIFFUSIVITY(
853c5e0e2c Jean*0518      I        bi,bj,1-OLx,sNx+OLx,1-OLy,sNy+OLy,
059d9fc14f Dimi*0519      I        GAD_SALINITY, .FALSE., .FALSE.,
853c5e0e2c Jean*0520      O        KPPdiffKzS(1-OLx,1-OLy,1,bi,bj),
059d9fc14f Dimi*0521      I        myThid)
                0522       CALL CALC_3D_DIFFUSIVITY(
853c5e0e2c Jean*0523      I        bi,bj,1-OLx,sNx+OLx,1-OLy,sNy+OLy,
059d9fc14f Dimi*0524      I        GAD_TEMPERATURE, .FALSE., .FALSE.,
853c5e0e2c Jean*0525      O        KPPdiffKzT(1-OLx,1-OLy,1,bi,bj),
059d9fc14f Dimi*0526      I        myThid)
e750a5e49e Mart*0527 #ifndef EXCLUDE_KPP_DOUBLEDIFF
                0528       IF ( KPPuseDoubleDiff ) THEN
a9d2e4c565 Jean*0529 C     Add the contribution of double diffusive effects (salt fingering
e750a5e49e Mart*0530 C     and diffusive convection) here. It would be more logical to add
a9d2e4c565 Jean*0531 C     them right after Ri_iwmix within kppmix, but ttalpha, ssbeta, theta
e750a5e49e Mart*0532 C     and salt are not passed to kppmix and are thus not available there.
                0533        CALL KPP_DOUBLEDIFF(
                0534      I      TTALPHA, SSBETA,
853c5e0e2c Jean*0535      U      KPPdiffKzT(1-OLx,1-OLy,1,bi,bj),
                0536      U      KPPdiffKzS(1-OLx,1-OLy,1,bi,bj),
edb6656069 Mart*0537      I      ikey,1-OLx,sNx+OLx,1-OLy,sNy+OLy,bi,bj,myThid)
e750a5e49e Mart*0538       ENDIF
                0539 #endif /* ndef EXCLUDE_KPP_DOUBLEDIFF */
00c7090dc0 Mart*0540 #ifdef ALLOW_AUTODIFF_TAMC
                0541 CADJ STORE KPPdiffKzS(:,:,:,bi,bj) = comlev1_kpp, key = ikey
                0542 CADJ STORE KPPdiffKzT(:,:,:,bi,bj) = comlev1_kpp, key = ikey
                0543 #endif
059d9fc14f Dimi*0544
2fa42a6013 Alis*0545       CALL KPPMIX (
00c7090dc0 Mart*0546      I       nzmax(1-OLx,1-OLy,bi,bj), shsq, dVsq, ustar,
edb6656069 Mart*0547      I       maskC(1-OLx,1-OLy,1,bi,bj),
                0548      I       bo, bosol,
63ceaaa79c Dimi*0549 #ifdef ALLOW_SALT_PLUME
edb6656069 Mart*0550      I       boplume, SaltPlumeDepth(1-OLx,1-OLy,bi,bj),
1f89baba18 Patr*0551 #ifdef SALT_PLUME_SPLIT_BASIN
edb6656069 Mart*0552      I       XC(1-OLx,1-OLy,bi,bj), YC(1-OLx,1-OLy,bi,bj),
1f89baba18 Patr*0553 #endif /* SALT_PLUME_SPLIT_BASIN */
63ceaaa79c Dimi*0554 #endif /* ALLOW_SALT_PLUME */
00c7090dc0 Mart*0555      I       dbloc, Ritop, fCori(1-OLx,1-OLy,bi,bj),
                0556 #ifdef SHORTWAVE_HEATING
                0557      I       SWFRac3D  (1-OLx,1-OLy,1,bi,bj),
                0558 #endif
edb6656069 Mart*0559      I       KPPdiffKzS(1-OLx,1-OLy,1,bi,bj),
                0560      I       KPPdiffKzT(1-OLx,1-OLy,1,bi,bj),
                0561      I       ikey,
                0562      O       vddiff,
                0563      U       ghat,
                0564      O       hbl,
00c7090dc0 Mart*0565 #ifdef SHORTWAVE_HEATING
                0566      O       kbl,
                0567 #endif
edb6656069 Mart*0568      I       bi, bj, myTime, myIter, myThid )
2fa42a6013 Alis*0569
                0570 c-----------------------------------------------------------------------
956c0a5824 Patr*0571 c     zero out land values and transfer to global variables
2fa42a6013 Alis*0572 c-----------------------------------------------------------------------
                0573
a67734cba4 Mart*0574       DO k = 1, Nr
                0575        km1 = max(1,k-1)
7c50f07931 Mart*0576        DO j = jMin, jMax
                0577         DO i = iMin, iMax
c572a3ecb8 Jean*0578          KPPviscAz(i,j,k,bi,bj) = vddiff(i,j,k-1,1) * maskC(i,j,k,bi,bj)
36aa272a96 Mart*0579      &        * maskC(i,j,km1,bi,bj)
c572a3ecb8 Jean*0580          KPPdiffKzS(i,j,k,bi,bj)= vddiff(i,j,k-1,2) * maskC(i,j,k,bi,bj)
36aa272a96 Mart*0581      &        * maskC(i,j,km1,bi,bj)
c572a3ecb8 Jean*0582          KPPdiffKzT(i,j,k,bi,bj)= vddiff(i,j,k-1,3) * maskC(i,j,k,bi,bj)
36aa272a96 Mart*0583      &        * maskC(i,j,km1,bi,bj)
c572a3ecb8 Jean*0584          KPPghat(i,j,k,bi,bj)   = ghat(i,j,k)       * maskC(i,j,k,bi,bj)
36aa272a96 Mart*0585      &        * maskC(i,j,km1,bi,bj)
c417c38d79 Jean*0586         ENDDO
a67734cba4 Mart*0587        ENDDO
                0588       ENDDO
7c50f07931 Mart*0589       kSurf = 1
                0590       DO j = jMin, jMax
                0591        DO i = iMin, iMax
36aa272a96 Mart*0592 #ifdef ALLOW_SHELFICE
7c50f07931 Mart*0593         IF ( useShelfIce ) kSurf = MAX(1,kTopC(i,j,bi,bj))
36aa272a96 Mart*0594 #endif /* ALLOW_SHELFICE */
2b4c90c108 Mart*0595         KPPhbl(i,j,bi,bj) = hbl(i,j) * maskC(i,j,kSurf,bi,bj)
c417c38d79 Jean*0596        ENDDO
                0597       ENDDO
2fa42a6013 Alis*0598
                0599 #ifdef KPP_SMOOTH_VISC
                0600 c     horizontal smoothing of vertical viscosity
                0601       DO k = 1, Nr
956c0a5824 Patr*0602          CALL SMOOTH_HORIZ (
2fa42a6013 Alis*0603      I        k, bi, bj,
c417c38d79 Jean*0604      U        KPPviscAz(1-OLx,1-OLy,k,bi,bj),
25c3477f99 Mart*0605      I        myThid )
c417c38d79 Jean*0606       ENDDO
                0607 C jmc: No EXCH inside bi,bj loop !!!
6637358eea Jean*0608 c     _EXCH_XYZ_RL(KPPviscAz  , myThid )
2fa42a6013 Alis*0609 #endif /* KPP_SMOOTH_VISC */
                0610
                0611 #ifdef KPP_SMOOTH_DIFF
                0612 c     horizontal smoothing of vertical diffusivity
                0613       DO k = 1, Nr
956c0a5824 Patr*0614          CALL SMOOTH_HORIZ (
2fa42a6013 Alis*0615      I        k, bi, bj,
c417c38d79 Jean*0616      U        KPPdiffKzS(1-OLx,1-OLy,k,bi,bj),
25c3477f99 Mart*0617      I        myThid )
956c0a5824 Patr*0618          CALL SMOOTH_HORIZ (
2fa42a6013 Alis*0619      I        k, bi, bj,
c417c38d79 Jean*0620      U        KPPdiffKzT(1-OLx,1-OLy,k,bi,bj),
25c3477f99 Mart*0621      I        myThid )
c417c38d79 Jean*0622       ENDDO
2fa42a6013 Alis*0623 #endif /* KPP_SMOOTH_DIFF */
                0624
32e2940ee8 Patr*0625 cph(
                0626 cph  crucial: this avoids full recomp./call of kppmix
853c5e0e2c Jean*0627 #ifdef ALLOW_AUTODIFF_TAMC
00c7090dc0 Mart*0628 CADJ STORE KPPhbl(:,:,bi,bj) = comlev1_kpp, key = ikey
                0629 # ifdef SHORTWAVE_HEATING
                0630 CADJ STORE kbl = comlev1_kpp, key = ikey
                0631 # endif
853c5e0e2c Jean*0632 #endif /* ALLOW_AUTODIFF_TAMC */
32e2940ee8 Patr*0633 cph)
                0634
00c7090dc0 Mart*0635 #ifdef SHORTWAVE_HEATING
                0636       IF ( selectPenetratingSW .GE. 1 ) THEN
                0637 #ifdef ALLOW_AUTODIFF_TAMC
                0638 C     Avoid uselessly recomputing this block because of re-using worka
                0639 CADJ INCOMPLETE worka
                0640 #endif
                0641        IF ( KPPuseSWfrac3D ) THEN
                0642 C     Determine the fraction of solar shortwave flux penetrating to the
                0643 C     bottom of the mixed layer by interpolation.
                0644         DO j=1-OLy,sNy+OLy
                0645          DO i=1-OLx,sNx+OLx
                0646           k = kbl(i,j)
                0647           rFac = MAX( (KPPhbl(i,j,bi,bj)+rF(k))*recip_drF(k), zeroRL )
                0648           worka(i,j) = SWFrac3D(i,j,k,bi,bj)
                0649      &         + rFac*(SWFrac3D(i,j,k+1,bi,bj) - SWFrac3D(i,j,k,bi,bj))
                0650          ENDDO
                0651         ENDDO
                0652        ELSE
2fa42a6013 Alis*0653 C     Compute fraction of solar short-wave flux penetrating to
                0654 C     the bottom of the mixing layer.
00c7090dc0 Mart*0655         DO j=1-OLy,sNy+OLy
2fa42a6013 Alis*0656          DO i=1-OLx,sNx+OLx
00c7090dc0 Mart*0657           worka(i,j) = KPPhbl(i,j,bi,bj)
2fa42a6013 Alis*0658          ENDDO
00c7090dc0 Mart*0659         ENDDO
                0660         CALL SWFRAC(
956c0a5824 Patr*0661      I     (sNx+2*OLx)*(sNy+2*OLy), minusone,
c417c38d79 Jean*0662      U     worka,
                0663      I     myTime, myIter, myThid )
00c7090dc0 Mart*0664        ENDIF
                0665        DO j=1-OLy,sNy+OLy
                0666         DO i=1-OLx,sNx+OLx
                0667          KPPfrac(i,j,bi,bj) = 1. _d 0 - worka(i,j)
                0668         ENDDO
                0669        ENDDO
                0670       ENDIF
                0671 #endif /* SHORTWAVE_HEATING */
2fa42a6013 Alis*0672
63ceaaa79c Dimi*0673 #ifdef ALLOW_SALT_PLUME
                0674 C     Compute fraction of saltplume (flux) penetrating to
                0675 C     the bottom of the mixing layer.
00c7090dc0 Mart*0676 #ifdef ALLOW_AUTODIFF_TAMC
                0677 C     Avoid uselessly recomputing the previous block because of worka
                0678 CADJ INCOMPLETE worka
                0679 #endif
30c6f5b1cd An T*0680       IF ( useSALT_PLUME ) THEN
1f89baba18 Patr*0681 #ifndef SALT_PLUME_VOLUME
a67734cba4 Mart*0682        DO j=1-OLy,sNy+OLy
                0683         DO i=1-OLx,sNx+OLx
                0684          worka(i,j) = KPPhbl(i,j,bi,bj)
1f89baba18 Patr*0685 #ifdef SALT_PLUME_SPLIT_BASIN
a67734cba4 Mart*0686          lon(i,j) = XC(i,j,bi,bj)
                0687          lat(i,j) = YC(i,j,bi,bj)
1f89baba18 Patr*0688 #endif /* SALT_PLUME_SPLIT_BASIN */
30c6f5b1cd An T*0689         ENDDO
a67734cba4 Mart*0690        ENDDO
                0691        CALL SALT_PLUME_FRAC(
00c7090dc0 Mart*0692      I      (sNx+2*OLx)*(sNy+2*OLy), minusone,
                0693      I      SaltPlumeDepth(1-OLx,1-OLy,bi,bj),
1f89baba18 Patr*0694 #ifdef SALT_PLUME_SPLIT_BASIN
a67734cba4 Mart*0695      I      lon,lat,
1f89baba18 Patr*0696 #endif /* SALT_PLUME_SPLIT_BASIN */
a67734cba4 Mart*0697      U      worka,
                0698      I      myTime, myIter, myThid )
                0699        DO j=1-OLy,sNy+OLy
                0700         DO i=1-OLx,sNx+OLx
                0701          KPPplumefrac(i,j,bi,bj) = 1. _d 0 - worka(i,j)
30c6f5b1cd An T*0702         ENDDO
a67734cba4 Mart*0703        ENDDO
1f89baba18 Patr*0704 #else /* SALT_PLUME_VOLUME */
a67734cba4 Mart*0705        DO j=1-OLy,sNy+OLy
                0706         DO i=1-OLx,sNx+OLx
                0707          KPPplumefrac(i,j,bi,bj) = 0. _d 0
1f89baba18 Patr*0708         ENDDO
a67734cba4 Mart*0709        ENDDO
1f89baba18 Patr*0710 #endif /* SALT_PLUME_VOLUME */
30c6f5b1cd An T*0711       ENDIF
63ceaaa79c Dimi*0712 #endif /* ALLOW_SALT_PLUME */
                0713
2fa42a6013 Alis*0714       ENDIF
                0715
3587afcb9b Alis*0716 #endif /* ALLOW_KPP */
2fa42a6013 Alis*0717
                0718       RETURN
                0719       END
b7230be7a5 Patr*0720
a9d2e4c565 Jean*0721 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
2b4c90c108 Mart*0722 CBOP
                0723 C !ROUTINE: KPP_CALC_DUMMY
a9d2e4c565 Jean*0724
2b4c90c108 Mart*0725 C !INTERFACE: ==========================================================
c417c38d79 Jean*0726       SUBROUTINE KPP_CALC_DUMMY(
                0727      I     bi, bj, myTime, myIter, myThid )
2b4c90c108 Mart*0728
                0729 C !DESCRIPTION: \bv
c417c38d79 Jean*0730 C     *==========================================================*
b7230be7a5 Patr*0731 C     | SUBROUTINE KPP_CALC_DUMMY                                |
                0732 C     | o Compute all KPP fields defined in KPP.h                |
b7b61e618a Mart*0733 C     | o Dummy routine for TAF                                  |
c417c38d79 Jean*0734 C     *==========================================================*
b7230be7a5 Patr*0735
2b4c90c108 Mart*0736 C !USES: ===============================================================
                0737       IMPLICIT NONE
b7230be7a5 Patr*0738 #include "SIZE.h"
                0739 #include "EEPARAMS.h"
                0740 #include "PARAMS.h"
                0741 #include "KPP.h"
                0742 #include "KPP_PARAMS.h"
                0743 #include "GRID.h"
059d9fc14f Dimi*0744 #include "GAD.h"
af61e5eb16 Mart*0745 #ifdef ALLOW_AUTODIFF_TAMC
                0746 # include "tamc.h"
                0747 #endif /* ALLOW_AUTODIFF_TAMC */
b7230be7a5 Patr*0748
2b4c90c108 Mart*0749 C !INPUT PARAMETERS: ===================================================
c417c38d79 Jean*0750 c     bi, bj :: Current tile indices
                0751 c     myTime :: Current time in simulation
                0752 c     myIter :: Current iteration number in simulation
                0753 c     myThid :: My Thread Id. number
b7230be7a5 Patr*0754       INTEGER bi, bj
                0755       _RL     myTime
c417c38d79 Jean*0756       INTEGER myIter
                0757       INTEGER myThid
b7230be7a5 Patr*0758
                0759 #ifdef ALLOW_KPP
2b4c90c108 Mart*0760 C !LOCAL VARIABLES: ====================================================
                0761       INTEGER i, j, k
af61e5eb16 Mart*0762 #if ( !defined EXCLUDE_PCELL_MIX_CODE && defined  ALLOW_AUTODIFF_TAMC )
                0763       INTEGER ikey
                0764 #endif
2b4c90c108 Mart*0765 CEOP
b7230be7a5 Patr*0766
                0767       DO j=1-OLy,sNy+OLy
                0768          DO i=1-OLx,sNx+OLx
                0769             KPPhbl (i,j,bi,bj) = 1.0
                0770             KPPfrac(i,j,bi,bj) = 0.0
63ceaaa79c Dimi*0771 #ifdef ALLOW_SALT_PLUME
                0772             KPPplumefrac(i,j,bi,bj) = 0.0
                0773 #endif /* ALLOW_SALT_PLUME */
b7230be7a5 Patr*0774             DO k = 1,Nr
                0775                KPPghat   (i,j,k,bi,bj) = 0.0
a9d2e4c565 Jean*0776                KPPviscAz (i,j,k,bi,bj) = viscArNr(1)
b7230be7a5 Patr*0777             ENDDO
                0778          ENDDO
                0779       ENDDO
059d9fc14f Dimi*0780
af61e5eb16 Mart*0781 #if ( !defined EXCLUDE_PCELL_MIX_CODE && defined  ALLOW_AUTODIFF_TAMC )
                0782       ikey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
                0783 CADJ STORE KPPdiffKzS(:,:,:,bi,bj) = comlev1_kpp, key = ikey
                0784 CADJ STORE KPPdiffKzT(:,:,:,bi,bj) = comlev1_kpp, key = ikey
                0785 #endif
059d9fc14f Dimi*0786       CALL CALC_3D_DIFFUSIVITY(
853c5e0e2c Jean*0787      I     bi,bj,1-OLx,sNx+OLx,1-OLy,sNy+OLy,
059d9fc14f Dimi*0788      I     GAD_SALINITY, .FALSE., .FALSE.,
853c5e0e2c Jean*0789      O     KPPdiffKzS(1-OLx,1-OLy,1,bi,bj),
059d9fc14f Dimi*0790      I     myThid)
                0791       CALL CALC_3D_DIFFUSIVITY(
853c5e0e2c Jean*0792      I     bi,bj,1-OLx,sNx+OLx,1-OLy,sNy+OLy,
059d9fc14f Dimi*0793      I     GAD_TEMPERATURE, .FALSE., .FALSE.,
853c5e0e2c Jean*0794      O     KPPdiffKzT(1-OLx,1-OLy,1,bi,bj),
059d9fc14f Dimi*0795      I     myThid)
                0796
b7230be7a5 Patr*0797 #endif
                0798       RETURN
                0799       END