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.
                0012 C--  o KPP_CALC_DUMMY   - Dummy routine for TAMC
                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)
                0097 c     KPPfrac     - Fraction of short-wave flux penetrating 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
                0156 c     work?  (nx,ny)       - horizontal working arrays
1f89baba18 Patr*0157 c     temp?  (nx,ny,Nr)    - 3d working arrays
2fa42a6013 Alis*0158 c     ustar  (nx,ny)       - surface friction velocity                  (m/s)
                0159 c     bo     (nx,ny)       - surface turbulent buoyancy forcing     (m^2/s^3)
                0160 c     bosol  (nx,ny)       - surface radiative buoyancy forcing     (m^2/s^3)
1f89baba18 Patr*0161 c     boplume(nx,ny,Nrp1)  - surface haline buoyancy forcing        (m^2/s^3)
2fa42a6013 Alis*0162 c     shsq   (nx,ny,Nr)    - local velocity shear squared
                0163 c                            at interfaces for ri_iwmix             (m^2/s^2)
                0164 c     dVsq   (nx,ny,Nr)    - velocity shear re surface squared
                0165 c                            at grid levels for bldepth             (m^2/s^2)
                0166 c     dbloc  (nx,ny,Nr)    - local delta buoyancy at interfaces
                0167 c                            for ri_iwmix and bldepth                 (m/s^2)
                0168 c     Ritop  (nx,ny,Nr)    - numerator of bulk richardson number
                0169 c                            at grid levels for bldepth
                0170 c     vddiff (nx,ny,Nrp2,1)- vertical viscosity on "t-grid"           (m^2/s)
42e0e34daa Dimi*0171 c     vddiff (nx,ny,Nrp2,2)- vert. diff. on next row for salt&tracers (m^2/s)
                0172 c     vddiff (nx,ny,Nrp2,3)- vert. diff. on next row for temperature  (m^2/s)
2fa42a6013 Alis*0173 c     ghat   (nx,ny,Nr)    - nonlocal transport coefficient           (s/m^2)
                0174 c     hbl    (nx,ny)       - mixing layer depth                           (m)
                0175 c     kmtj   (nx,ny)       - maximum number of wet levels in each column
2b4c90c108 Mart*0176       INTEGER work1 ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy            )
fe66051ebd Dimi*0177       _RL worka ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
                0178       _RL work2 ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
                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                )
63ceaaa79c Dimi*0182 #ifdef ALLOW_SALT_PLUME
1f89baba18 Patr*0183       _RL temp1   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr          )
                0184       _RL temp2   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr          )
                0185       _RL boplume ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1        )
                0186 #ifdef SALT_PLUME_SPLIT_BASIN
                0187       _RL lon ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                  )
                0188       _RL lat ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                  )
                0189 #endif /* SALT_PLUME_SPLIT_BASIN */
63ceaaa79c Dimi*0190 #endif /* ALLOW_SALT_PLUME */
fe66051ebd Dimi*0191       _RL shsq  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
                0192       _RL dVsq  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
                0193       _RL dbloc ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
                0194       _RL Ritop ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
                0195       _RL vddiff( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, 0:Nrp1, mdiff )
                0196       _RL ghat  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
                0197       _RL hbl   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
277d95ff7e Patr*0198 cph(
fe66051ebd Dimi*0199       _RL TTALPHA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
                0200       _RL SSBETA ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
277d95ff7e Patr*0201 cph)
7c50f07931 Mart*0202       integer i, j, kSurf, k, kp1, km1, im1, ip1, jm1, jp1
edb6656069 Mart*0203 C     ikey :: tape key (depends on tiles)
                0204       integer ikey
853c5e0e2c Jean*0205 CEOP
2fa42a6013 Alis*0206
ab9b91ea7a Patr*0207 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0208       ikey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
853c5e0e2c Jean*0209 #else /* ALLOW_AUTODIFF_TAMC */
edb6656069 Mart*0210       ikey = 0
ab9b91ea7a Patr*0211 #endif /* ALLOW_AUTODIFF_TAMC */
                0212
2fa42a6013 Alis*0213 c     Check to see if new vertical mixing coefficient should be computed now?
94a46dfe0d Jean*0214       IF ( DIFFERENT_MULTIPLE(kpp_freq,myTime,deltaTClock)
adf557b193 Jean*0215      1     .OR. myTime .EQ. startTime ) THEN
c417c38d79 Jean*0216
2fa42a6013 Alis*0217 c-----------------------------------------------------------------------
                0218 c     prepare input arrays for subroutine "kppmix" to compute
                0219 c     viscosity and diffusivity and ghat.
                0220 c     All input arrays need to be in m-k-s units.
                0221 c
                0222 c     note: for the computation of the bulk richardson number in the
                0223 c     "bldepth" subroutine, gradients of velocity and buoyancy are
                0224 c     required at every depth. in the case of very fine vertical grids
                0225 c     (thickness of top layer < 2m), the surface reference depth must
                0226 c     be set to zref=epsilon/2*zgrid(k), and the reference value
                0227 c     of velocity and buoyancy must be computed as vertical average
                0228 c     between the surface and 2*zref.  in the case of coarse vertical
                0229 c     grids zref is zgrid(1)/2., and the surface reference value is
                0230 c     simply the surface value at zgrid(1).
                0231 c-----------------------------------------------------------------------
                0232
                0233 c------------------------------------------------------------------------
                0234 c     density related quantities
                0235 c     --------------------------
                0236 c
                0237 c      work2   - density of surface layer                        (kg/m^3)
                0238 c      dbloc   - local buoyancy gradient at Nr interfaces
                0239 c                g/rho{k+1,k+1} * [ drho{k,k+1}-drho{k+1,k+1} ]   (m/s^2)
                0240 c      dbsfc (stored in Ritop to conserve stack memory)
                0241 c              - buoyancy difference with respect to the surface
                0242 c                g * [ drho{1,k}/rho{1,k} - drho{k,k}/rho{k,k} ]  (m/s^2)
                0243 c      ttalpha (stored in vddiff(:,:,:,1) to conserve stack memory)
                0244 c              - thermal expansion coefficient without 1/rho factor
                0245 c                d(rho{k,k})/d(T(k))                           (kg/m^3/C)
                0246 c      ssbeta (stored in vddiff(:,:,:,2) to conserve stack memory)
                0247 c              - salt expansion coefficient without 1/rho factor
ba0b047096 Mart*0248 c                d(rho{k,k})/d(S(k))                    ((kg/m^3)/(g/kg))
2fa42a6013 Alis*0249 c------------------------------------------------------------------------
                0250
                0251       CALL STATEKPP(
25c3477f99 Mart*0252      O     work2, dbloc, Ritop,
                0253      O     TTALPHA, SSBETA,
edb6656069 Mart*0254      I     ikey, bi, bj, myThid )
2fa42a6013 Alis*0255
                0256       DO k = 1, Nr
a67734cba4 Mart*0257        DO j = 1-OLy, sNy+OLy
                0258         DO i = 1-OLx, sNx+OLx
                0259          ghat(i,j,k) = dbloc(i,j,k)
                0260         ENDDO
                0261        ENDDO
2fa42a6013 Alis*0262       ENDDO
                0263
956c0a5824 Patr*0264 #ifdef KPP_SMOOTH_DBLOC
                0265 c     horizontally smooth dbloc with a 121 filter
                0266 c     smooth dbloc stored in ghat to save space
                0267 c     dbloc(k) is buoyancy gradientnote between k and k+1
                0268 c     levels therefore k+1 mask must be used
                0269
                0270       DO k = 1, Nr-1
a67734cba4 Mart*0271        CALL SMOOTH_HORIZ (
                0272      I      k+1, bi, bj,
                0273      U      ghat (1-OLx,1-OLy,k),
                0274      I      myThid )
956c0a5824 Patr*0275       ENDDO
                0276
2fa42a6013 Alis*0277 #endif /* KPP_SMOOTH_DBLOC */
                0278
                0279 #ifdef KPP_SMOOTH_DENS
                0280 c     horizontally smooth density related quantities with 121 filters
fe66051ebd Dimi*0281       CALL SMOOTH_HORIZ (
956c0a5824 Patr*0282      I     1, bi, bj,
c417c38d79 Jean*0283      U     work2,
25c3477f99 Mart*0284      I     myThid )
2fa42a6013 Alis*0285       DO k = 1, Nr
fe66051ebd Dimi*0286          CALL SMOOTH_HORIZ (
956c0a5824 Patr*0287      I        k+1, bi, bj,
c417c38d79 Jean*0288      U        dbloc (1-OLx,1-OLy,k),
25c3477f99 Mart*0289      I        myThid )
fe66051ebd Dimi*0290          CALL SMOOTH_HORIZ (
2fa42a6013 Alis*0291      I        k, bi, bj,
c417c38d79 Jean*0292      U        Ritop (1-OLx,1-OLy,k),
25c3477f99 Mart*0293      I        myThid )
fe66051ebd Dimi*0294          CALL SMOOTH_HORIZ (
2fa42a6013 Alis*0295      I        k, bi, bj,
c417c38d79 Jean*0296      U        TTALPHA(1-OLx,1-OLy,k),
25c3477f99 Mart*0297      I        myThid )
fe66051ebd Dimi*0298          CALL SMOOTH_HORIZ (
2fa42a6013 Alis*0299      I        k, bi, bj,
c417c38d79 Jean*0300      U        SSBETA(1-OLx,1-OLy,k),
25c3477f99 Mart*0301      I        myThid )
2fa42a6013 Alis*0302       ENDDO
                0303 #endif /* KPP_SMOOTH_DENS */
                0304
7c50f07931 Mart*0305 c     surface index
                0306       kSurf = 1
2fa42a6013 Alis*0307       DO k = 1, Nr
a67734cba4 Mart*0308        km1 = max(1,k-1)
2b4c90c108 Mart*0309        kp1 = min(Nr,k+1)
a67734cba4 Mart*0310        DO j = 1-OLy, sNy+OLy
                0311         DO i = 1-OLx, sNx+OLx
2fa42a6013 Alis*0312
                0313 c     zero out dbloc over land points (so that the convective
2b4c90c108 Mart*0314 c     part of the interior mixing can be diagnosed).
                0315 c     ML: Note that dbloc and ghat (=smoothed dbloc) are defined at
                0316 c     the bottom of the cell rather than the (MITgcm-standard) top
                0317 c     of the cell; therefore the proper masking is maskC(k)*maskC(k+1)
a67734cba4 Mart*0318          dbloc(i,j,k) = dbloc(i,j,k) * maskC(i,j,k,bi,bj)
2b4c90c108 Mart*0319      &        * maskC(i,j,kp1,bi,bj)
a67734cba4 Mart*0320          ghat(i,j,k)  = ghat(i,j,k)  * maskC(i,j,k,bi,bj)
2b4c90c108 Mart*0321      &        * maskC(i,j,kp1,bi,bj)
7c50f07931 Mart*0322 #ifdef ALLOW_SHELFICE
                0323          IF ( useShelfIce ) kSurf = MAX(1,kTopC(i,j,bi,bj))
                0324 #endif /* ALLOW_SHELFICE */
a67734cba4 Mart*0325          Ritop(i,j,k) = Ritop(i,j,k) * maskC(i,j,k,bi,bj)
2b4c90c108 Mart*0326      &        * maskC(i,j,kSurf,bi,bj)
                0327          IF(k.EQ.nzmax(i,j,bi,bj)) THEN
a67734cba4 Mart*0328           dbloc(i,j,k) = p0
                0329           ghat(i,j,k)  = p0
                0330           Ritop(i,j,k) = p0
2b4c90c108 Mart*0331          ENDIF
2fa42a6013 Alis*0332
                0333 c     numerator of bulk richardson number on grid levels
                0334 c     note: land and ocean bottom values need to be set to zero
                0335 c     so that the subroutine "bldepth" works correctly
a67734cba4 Mart*0336          Ritop(i,j,k) = (zgrid(1)-zgrid(k)) * Ritop(i,j,k)
2fa42a6013 Alis*0337
a67734cba4 Mart*0338         ENDDO
                0339        ENDDO
c417c38d79 Jean*0340       ENDDO
2fa42a6013 Alis*0341
32e2940ee8 Patr*0342 cph(
                0343 cph  this avoids a single or double recomp./call of statekpp
853c5e0e2c Jean*0344 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0345 CADJ STORE work2              = comlev1_kpp, key = ikey
3067745c9b Patr*0346 #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
edb6656069 Mart*0347 CADJ STORE dbloc, Ritop, ghat = comlev1_kpp, key = ikey
                0348 CADJ STORE vddiff             = comlev1_kpp, key = ikey
                0349 CADJ STORE TTALPHA, SSBETA    = comlev1_kpp, key = ikey
32e2940ee8 Patr*0350 #endif
853c5e0e2c Jean*0351 #endif /* ALLOW_AUTODIFF_TAMC */
32e2940ee8 Patr*0352 cph)
                0353
36aa272a96 Mart*0354 CML#ifdef ALLOW_SHELFICE
a9d2e4c565 Jean*0355 CMLC     For the pbl parameterisation to work underneath the ice shelves
36aa272a96 Mart*0356 CMLC     it needs to know the surface (ice-ocean) fluxes. However, masking
a9d2e4c565 Jean*0357 CMLC     and indexing problems make this part of the code not work
36aa272a96 Mart*0358 CMLC     underneath the ice shelves and the following lines are only here
                0359 CMLC     to remind me that this still needs to be sorted out.
6af7e0efdf Mart*0360 CML      shelfIceFac = 0. _d 0
                0361 CML      IF ( useShelfIce ) selfIceFac = 1. _d 0
7c50f07931 Mart*0362 CML      DO j = jMin, jMax
                0363 CML       DO i = iMin, iMax
a9d2e4c565 Jean*0364 CML        surfForcT = surfaceForcingT(i,j,bi,bj)
6af7e0efdf Mart*0365 CML     &       + shelficeForcingT(i,j,bi,bj) * shelfIceFac
a9d2e4c565 Jean*0366 CML        surfForcS = surfaceForcingS(i,j,bi,bj)
6af7e0efdf Mart*0367 CML     &       + shelficeForcingS(i,j,bi,bj) * shelfIceFac
c417c38d79 Jean*0368 CML       ENDDO
                0369 CML      ENDDO
36aa272a96 Mart*0370 CML#endif /* ALLOW_SHELFICE */
32e2940ee8 Patr*0371
2fa42a6013 Alis*0372 c------------------------------------------------------------------------
6af7e0efdf Mart*0373 c     friction velocity, turbulent and radiative surface buoyancy forcing
                0374 c     -------------------------------------------------------------------
                0375 c     taux / rho = surfaceForcingU                               (N/m^2)
                0376 c     tauy / rho = surfaceForcingV                               (N/m^2)
                0377 c     ustar = sqrt( sqrt( taux^2 + tauy^2 ) / rho )                (m/s)
                0378 c     bo    = - g * ( alpha*surfaceForcingT +
                0379 c                     beta *surfaceForcingS ) / rho            (m^2/s^3)
                0380 c     bosol = - g * alpha * Qsw * drF(1) / rho                 (m^2/s^3)
63ceaaa79c Dimi*0381 c     boplume = g * (beta * saltPlumeFlux/rhoConst ) /rho      (m^2/s^3)
1f89baba18 Patr*0382 c             = g * (beta * SPforcingS   /rhoConst ) /rho
                0383 c              +g * (alpha* SPforcingT   / ??
6af7e0efdf Mart*0384 c------------------------------------------------------------------------
2fa42a6013 Alis*0385 c     velocity shear
                0386 c     --------------
                0387 c     Get velocity shear squared, averaged from "u,v-grid"
                0388 c     onto "t-grid" (in (m/s)**2):
                0389 c     dVsq(k)=(Uref-U(k))**2+(Vref-V(k))**2      at grid levels
                0390 c     shsq(k)=(U(k)-U(k+1))**2+(V(k)-V(k+1))**2  at interfaces
a9d2e4c565 Jean*0391 c
6af7e0efdf Mart*0392 c     note: Vref can depend on the surface fluxes that is why we compute
                0393 c     dVsq in the subroutine that does the surface related stuff
                0394 c     (admittedly this is a bit messy)
2fa42a6013 Alis*0395 c------------------------------------------------------------------------
c417c38d79 Jean*0396
1f89baba18 Patr*0397 #ifdef ALLOW_SALT_PLUME
7c50f07931 Mart*0398       DO j=jMin,jMax
                0399        DO i=iMin,iMax
1f89baba18 Patr*0400 #ifndef SALT_PLUME_VOLUME
                0401         temp1(i,j,1) = saltPlumeFlux(i,j,bi,bj)
                0402         temp2(i,j,1) = 0. _d 0
                0403         DO k=2,Nr
                0404          temp1(i,j,k) = 0. _d 0
                0405          temp2(i,j,k) = 0. _d 0
                0406         ENDDO
                0407 #else /* def SALT_PLUME_VOLUME */
                0408         DO k=1,Nr
                0409          temp1(i,j,k) = SPforcingS(i,j,k,bi,bj)
                0410          temp2(i,j,k) = SPforcingT(i,j,k,bi,bj)
                0411         ENDDO
                0412 #endif /* SALT_PLUME_VOLUME */
                0413        ENDDO
                0414       ENDDO
                0415 #endif /* ALLOW_SALT_PLUME */
                0416
6af7e0efdf Mart*0417       CALL KPP_FORCING_SURF(
                0418      I     work2, surfaceForcingU, surfaceForcingV,
16fe9e8e73 Jean*0419      I     surfaceForcingT, surfaceForcingS, adjustColdSST_diag,
a9d2e4c565 Jean*0420      I     Qsw,
a67734cba4 Mart*0421 #ifdef KPP_ESTIMATE_UREF
                0422      I     dbloc,
                0423 #endif /* KPP_ESTIMATE_UREF */
63ceaaa79c Dimi*0424 #ifdef ALLOW_SALT_PLUME
1f89baba18 Patr*0425      I     temp1, temp2,
63ceaaa79c Dimi*0426 #endif /* ALLOW_SALT_PLUME */
                0427      I     ttalpha, ssbeta,
a9d2e4c565 Jean*0428      O     ustar, bo, bosol,
63ceaaa79c Dimi*0429 #ifdef ALLOW_SALT_PLUME
                0430      O     boplume,
                0431 #endif /* ALLOW_SALT_PLUME */
                0432      O     dVsq,
edb6656069 Mart*0433      I     ikey, iMin, iMax, jMin, jMax, bi, bj, myTime, myThid )
2fa42a6013 Alis*0434
6af7e0efdf Mart*0435 CMLcph(
853c5e0e2c Jean*0436 CML#ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0437 CMLCADJ STORE ustar = comlev1_kpp, key = ikey
853c5e0e2c Jean*0438 CML#endif
6af7e0efdf Mart*0439 CMLcph)
2fa42a6013 Alis*0440
6af7e0efdf Mart*0441 c initialize arrays to zero
2fa42a6013 Alis*0442       DO k = 1, Nr
6af7e0efdf Mart*0443        DO j = 1-OLy, sNy+OLy
                0444         DO i = 1-OLx, sNx+OLx
                0445          shsq(i,j,k) = p0
c417c38d79 Jean*0446         ENDDO
                0447        ENDDO
                0448       ENDDO
2fa42a6013 Alis*0449
                0450 c     shsq computation
                0451       DO k = 1, Nrm1
6af7e0efdf Mart*0452        kp1 = k + 1
7c50f07931 Mart*0453        DO j = jMin, jMax
6af7e0efdf Mart*0454         jm1 = j - 1
                0455         jp1 = j + 1
7c50f07931 Mart*0456         DO i = iMin, iMax
6af7e0efdf Mart*0457          im1 = i - 1
                0458          ip1 = i + 1
                0459          shsq(i,j,k) = p5 * (
c417c38d79 Jean*0460      &        (uVel(i,  j,  k,bi,bj)-uVel(i,  j,  kp1,bi,bj)) *
                0461      &        (uVel(i,  j,  k,bi,bj)-uVel(i,  j,  kp1,bi,bj)) +
                0462      &        (uVel(ip1,j,  k,bi,bj)-uVel(ip1,j,  kp1,bi,bj)) *
                0463      &        (uVel(ip1,j,  k,bi,bj)-uVel(ip1,j,  kp1,bi,bj)) +
                0464      &        (vVel(i,  j,  k,bi,bj)-vVel(i,  j,  kp1,bi,bj)) *
                0465      &        (vVel(i,  j,  k,bi,bj)-vVel(i,  j,  kp1,bi,bj)) +
                0466      &        (vVel(i,  jp1,k,bi,bj)-vVel(i,  jp1,kp1,bi,bj)) *
                0467      &        (vVel(i,  jp1,k,bi,bj)-vVel(i,  jp1,kp1,bi,bj)) )
2fa42a6013 Alis*0468 #ifdef KPP_SMOOTH_SHSQ
6af7e0efdf Mart*0469          shsq(i,j,k) = p5 * shsq(i,j,k) + p125 * (
c417c38d79 Jean*0470      &        (uVel(i,  jm1,k,bi,bj)-uVel(i,  jm1,kp1,bi,bj)) *
                0471      &        (uVel(i,  jm1,k,bi,bj)-uVel(i,  jm1,kp1,bi,bj)) +
                0472      &        (uVel(ip1,jm1,k,bi,bj)-uVel(ip1,jm1,kp1,bi,bj)) *
                0473      &        (uVel(ip1,jm1,k,bi,bj)-uVel(ip1,jm1,kp1,bi,bj)) +
                0474      &        (uVel(i,  jp1,k,bi,bj)-uVel(i,  jp1,kp1,bi,bj)) *
                0475      &        (uVel(i,  jp1,k,bi,bj)-uVel(i,  jp1,kp1,bi,bj)) +
                0476      &        (uVel(ip1,jp1,k,bi,bj)-uVel(ip1,jp1,kp1,bi,bj)) *
                0477      &        (uVel(ip1,jp1,k,bi,bj)-uVel(ip1,jp1,kp1,bi,bj)) +
                0478      &        (vVel(im1,j,  k,bi,bj)-vVel(im1,j,  kp1,bi,bj)) *
                0479      &        (vVel(im1,j,  k,bi,bj)-vVel(im1,j,  kp1,bi,bj)) +
                0480      &        (vVel(im1,jp1,k,bi,bj)-vVel(im1,jp1,kp1,bi,bj)) *
                0481      &        (vVel(im1,jp1,k,bi,bj)-vVel(im1,jp1,kp1,bi,bj)) +
                0482      &        (vVel(ip1,j,  k,bi,bj)-vVel(ip1,j,  kp1,bi,bj)) *
                0483      &        (vVel(ip1,j,  k,bi,bj)-vVel(ip1,j,  kp1,bi,bj)) +
                0484      &        (vVel(ip1,jp1,k,bi,bj)-vVel(ip1,jp1,kp1,bi,bj)) *
                0485      &        (vVel(ip1,jp1,k,bi,bj)-vVel(ip1,jp1,kp1,bi,bj)) )
2fa42a6013 Alis*0486 #endif
c417c38d79 Jean*0487         ENDDO
                0488        ENDDO
                0489       ENDDO
2fa42a6013 Alis*0490
2b4c90c108 Mart*0491 #ifdef ALLOW_DIAGNOSTICS
                0492       IF ( useDiagnostics ) THEN
                0493        CALL DIAGNOSTICS_FILL(shsq,  'KPPshsq ',0,Nr,2,bi,bj,myThid)
                0494       ENDIF
                0495 #endif /* ALLOW_DIAGNOSTICS */
32e2940ee8 Patr*0496 cph(
853c5e0e2c Jean*0497 #ifdef ALLOW_AUTODIFF_TAMC
3067745c9b Patr*0498 #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
edb6656069 Mart*0499 CADJ STORE dvsq, shsq = comlev1_kpp, key = ikey
32e2940ee8 Patr*0500 #endif
853c5e0e2c Jean*0501 #endif /* ALLOW_AUTODIFF_TAMC */
32e2940ee8 Patr*0502 cph)
                0503
2fa42a6013 Alis*0504 c-----------------------------------------------------------------------
                0505 c     solve for viscosity, diffusivity, ghat, and hbl on "t-grid"
                0506 c-----------------------------------------------------------------------
                0507
af61e5eb16 Mart*0508 #if ( !defined EXCLUDE_PCELL_MIX_CODE && defined  ALLOW_AUTODIFF_TAMC )
                0509 CADJ STORE KPPdiffKzS(:,:,:,bi,bj) = comlev1_kpp, key = ikey
                0510 CADJ STORE KPPdiffKzT(:,:,:,bi,bj) = comlev1_kpp, key = ikey
                0511 #endif
059d9fc14f Dimi*0512 c     precompute background vertical diffusivities, which are needed for
                0513 c     matching diffusivities at bottom of KPP PBL
                0514       CALL CALC_3D_DIFFUSIVITY(
853c5e0e2c Jean*0515      I        bi,bj,1-OLx,sNx+OLx,1-OLy,sNy+OLy,
059d9fc14f Dimi*0516      I        GAD_SALINITY, .FALSE., .FALSE.,
853c5e0e2c Jean*0517      O        KPPdiffKzS(1-OLx,1-OLy,1,bi,bj),
059d9fc14f Dimi*0518      I        myThid)
                0519       CALL CALC_3D_DIFFUSIVITY(
853c5e0e2c Jean*0520      I        bi,bj,1-OLx,sNx+OLx,1-OLy,sNy+OLy,
059d9fc14f Dimi*0521      I        GAD_TEMPERATURE, .FALSE., .FALSE.,
853c5e0e2c Jean*0522      O        KPPdiffKzT(1-OLx,1-OLy,1,bi,bj),
059d9fc14f Dimi*0523      I        myThid)
e750a5e49e Mart*0524 #ifndef EXCLUDE_KPP_DOUBLEDIFF
                0525       IF ( KPPuseDoubleDiff ) THEN
a9d2e4c565 Jean*0526 C     Add the contribution of double diffusive effects (salt fingering
e750a5e49e Mart*0527 C     and diffusive convection) here. It would be more logical to add
a9d2e4c565 Jean*0528 C     them right after Ri_iwmix within kppmix, but ttalpha, ssbeta, theta
e750a5e49e Mart*0529 C     and salt are not passed to kppmix and are thus not available there.
                0530        CALL KPP_DOUBLEDIFF(
                0531      I      TTALPHA, SSBETA,
853c5e0e2c Jean*0532      U      KPPdiffKzT(1-OLx,1-OLy,1,bi,bj),
                0533      U      KPPdiffKzS(1-OLx,1-OLy,1,bi,bj),
edb6656069 Mart*0534      I      ikey,1-OLx,sNx+OLx,1-OLy,sNy+OLy,bi,bj,myThid)
e750a5e49e Mart*0535       ENDIF
                0536 #endif /* ndef EXCLUDE_KPP_DOUBLEDIFF */
059d9fc14f Dimi*0537
4fa51dc730 Dimi*0538       DO j = 1-OLy, sNy+OLy
                0539          DO i = 1-OLx, sNx+OLx
2fa42a6013 Alis*0540             work1(i,j) = nzmax(i,j,bi,bj)
                0541             work2(i,j) = Fcori(i,j,bi,bj)
c417c38d79 Jean*0542          ENDDO
                0543       ENDDO
2fa42a6013 Alis*0544       CALL KPPMIX (
edb6656069 Mart*0545      I       work1, shsq, dVsq, ustar,
                0546      I       maskC(1-OLx,1-OLy,1,bi,bj),
                0547      I       bo, bosol,
63ceaaa79c Dimi*0548 #ifdef ALLOW_SALT_PLUME
edb6656069 Mart*0549      I       boplume, SaltPlumeDepth(1-OLx,1-OLy,bi,bj),
1f89baba18 Patr*0550 #ifdef SALT_PLUME_SPLIT_BASIN
edb6656069 Mart*0551      I       XC(1-OLx,1-OLy,bi,bj), YC(1-OLx,1-OLy,bi,bj),
1f89baba18 Patr*0552 #endif /* SALT_PLUME_SPLIT_BASIN */
63ceaaa79c Dimi*0553 #endif /* ALLOW_SALT_PLUME */
edb6656069 Mart*0554      I       dbloc, Ritop, work2,
                0555      I       KPPdiffKzS(1-OLx,1-OLy,1,bi,bj),
                0556      I       KPPdiffKzT(1-OLx,1-OLy,1,bi,bj),
                0557      I       ikey,
                0558      O       vddiff,
                0559      U       ghat,
                0560      O       hbl,
                0561      I       bi, bj, myTime, myIter, myThid )
2fa42a6013 Alis*0562
                0563 c-----------------------------------------------------------------------
956c0a5824 Patr*0564 c     zero out land values and transfer to global variables
2fa42a6013 Alis*0565 c-----------------------------------------------------------------------
                0566
a67734cba4 Mart*0567       DO k = 1, Nr
                0568        km1 = max(1,k-1)
7c50f07931 Mart*0569        DO j = jMin, jMax
                0570         DO i = iMin, iMax
c572a3ecb8 Jean*0571          KPPviscAz(i,j,k,bi,bj) = vddiff(i,j,k-1,1) * maskC(i,j,k,bi,bj)
36aa272a96 Mart*0572      &        * maskC(i,j,km1,bi,bj)
c572a3ecb8 Jean*0573          KPPdiffKzS(i,j,k,bi,bj)= vddiff(i,j,k-1,2) * maskC(i,j,k,bi,bj)
36aa272a96 Mart*0574      &        * maskC(i,j,km1,bi,bj)
c572a3ecb8 Jean*0575          KPPdiffKzT(i,j,k,bi,bj)= vddiff(i,j,k-1,3) * maskC(i,j,k,bi,bj)
36aa272a96 Mart*0576      &        * maskC(i,j,km1,bi,bj)
c572a3ecb8 Jean*0577          KPPghat(i,j,k,bi,bj)   = ghat(i,j,k)       * maskC(i,j,k,bi,bj)
36aa272a96 Mart*0578      &        * maskC(i,j,km1,bi,bj)
c417c38d79 Jean*0579         ENDDO
a67734cba4 Mart*0580        ENDDO
                0581       ENDDO
7c50f07931 Mart*0582       kSurf = 1
                0583       DO j = jMin, jMax
                0584        DO i = iMin, iMax
36aa272a96 Mart*0585 #ifdef ALLOW_SHELFICE
7c50f07931 Mart*0586         IF ( useShelfIce ) kSurf = MAX(1,kTopC(i,j,bi,bj))
36aa272a96 Mart*0587 #endif /* ALLOW_SHELFICE */
2b4c90c108 Mart*0588         KPPhbl(i,j,bi,bj) = hbl(i,j) * maskC(i,j,kSurf,bi,bj)
c417c38d79 Jean*0589        ENDDO
                0590       ENDDO
2fa42a6013 Alis*0591
                0592 #ifdef KPP_SMOOTH_VISC
                0593 c     horizontal smoothing of vertical viscosity
                0594       DO k = 1, Nr
956c0a5824 Patr*0595          CALL SMOOTH_HORIZ (
2fa42a6013 Alis*0596      I        k, bi, bj,
c417c38d79 Jean*0597      U        KPPviscAz(1-OLx,1-OLy,k,bi,bj),
25c3477f99 Mart*0598      I        myThid )
c417c38d79 Jean*0599       ENDDO
                0600 C jmc: No EXCH inside bi,bj loop !!!
6637358eea Jean*0601 c     _EXCH_XYZ_RL(KPPviscAz  , myThid )
2fa42a6013 Alis*0602 #endif /* KPP_SMOOTH_VISC */
                0603
                0604 #ifdef KPP_SMOOTH_DIFF
                0605 c     horizontal smoothing of vertical diffusivity
                0606       DO k = 1, Nr
956c0a5824 Patr*0607          CALL SMOOTH_HORIZ (
2fa42a6013 Alis*0608      I        k, bi, bj,
c417c38d79 Jean*0609      U        KPPdiffKzS(1-OLx,1-OLy,k,bi,bj),
25c3477f99 Mart*0610      I        myThid )
956c0a5824 Patr*0611          CALL SMOOTH_HORIZ (
2fa42a6013 Alis*0612      I        k, bi, bj,
c417c38d79 Jean*0613      U        KPPdiffKzT(1-OLx,1-OLy,k,bi,bj),
25c3477f99 Mart*0614      I        myThid )
c417c38d79 Jean*0615       ENDDO
2fa42a6013 Alis*0616 #endif /* KPP_SMOOTH_DIFF */
                0617
32e2940ee8 Patr*0618 cph(
                0619 cph  crucial: this avoids full recomp./call of kppmix
853c5e0e2c Jean*0620 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0621 CADJ STORE KPPhbl = comlev1_kpp, key = ikey
853c5e0e2c Jean*0622 #endif /* ALLOW_AUTODIFF_TAMC */
32e2940ee8 Patr*0623 cph)
                0624
2fa42a6013 Alis*0625 C     Compute fraction of solar short-wave flux penetrating to
                0626 C     the bottom of the mixing layer.
                0627       DO j=1-OLy,sNy+OLy
                0628          DO i=1-OLx,sNx+OLx
                0629             worka(i,j) = KPPhbl(i,j,bi,bj)
                0630          ENDDO
                0631       ENDDO
                0632       CALL SWFRAC(
956c0a5824 Patr*0633      I     (sNx+2*OLx)*(sNy+2*OLy), minusone,
c417c38d79 Jean*0634      U     worka,
                0635      I     myTime, myIter, myThid )
2fa42a6013 Alis*0636       DO j=1-OLy,sNy+OLy
                0637          DO i=1-OLx,sNx+OLx
956c0a5824 Patr*0638             KPPfrac(i,j,bi,bj) = worka(i,j)
2fa42a6013 Alis*0639          ENDDO
                0640       ENDDO
                0641
63ceaaa79c Dimi*0642 #ifdef ALLOW_SALT_PLUME
                0643 C     Compute fraction of saltplume (flux) penetrating to
                0644 C     the bottom of the mixing layer.
30c6f5b1cd An T*0645       IF ( useSALT_PLUME ) THEN
1f89baba18 Patr*0646 #ifndef SALT_PLUME_VOLUME
a67734cba4 Mart*0647        DO j=1-OLy,sNy+OLy
                0648         DO i=1-OLx,sNx+OLx
                0649          work2(i,j) = SaltPlumeDepth(i,j,bi,bj)
                0650          worka(i,j) = KPPhbl(i,j,bi,bj)
1f89baba18 Patr*0651 #ifdef SALT_PLUME_SPLIT_BASIN
a67734cba4 Mart*0652          lon(i,j) = XC(i,j,bi,bj)
                0653          lat(i,j) = YC(i,j,bi,bj)
1f89baba18 Patr*0654 #endif /* SALT_PLUME_SPLIT_BASIN */
30c6f5b1cd An T*0655         ENDDO
a67734cba4 Mart*0656        ENDDO
                0657        CALL SALT_PLUME_FRAC(
                0658      I      (sNx+2*OLx)*(sNy+2*OLy), minusone, work2,
1f89baba18 Patr*0659 #ifdef SALT_PLUME_SPLIT_BASIN
a67734cba4 Mart*0660      I      lon,lat,
1f89baba18 Patr*0661 #endif /* SALT_PLUME_SPLIT_BASIN */
a67734cba4 Mart*0662      U      worka,
                0663      I      myTime, myIter, myThid )
                0664        DO j=1-OLy,sNy+OLy
                0665         DO i=1-OLx,sNx+OLx
                0666          KPPplumefrac(i,j,bi,bj) = 1. _d 0 - worka(i,j)
30c6f5b1cd An T*0667         ENDDO
a67734cba4 Mart*0668        ENDDO
1f89baba18 Patr*0669 #else /* SALT_PLUME_VOLUME */
a67734cba4 Mart*0670        DO j=1-OLy,sNy+OLy
                0671         DO i=1-OLx,sNx+OLx
                0672          KPPplumefrac(i,j,bi,bj) = 0. _d 0
1f89baba18 Patr*0673         ENDDO
a67734cba4 Mart*0674        ENDDO
1f89baba18 Patr*0675 #endif /* SALT_PLUME_VOLUME */
30c6f5b1cd An T*0676       ENDIF
63ceaaa79c Dimi*0677 #endif /* ALLOW_SALT_PLUME */
                0678
2fa42a6013 Alis*0679       ENDIF
                0680
3587afcb9b Alis*0681 #endif /* ALLOW_KPP */
2fa42a6013 Alis*0682
                0683       RETURN
                0684       END
b7230be7a5 Patr*0685
a9d2e4c565 Jean*0686 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
2b4c90c108 Mart*0687 CBOP
                0688 C !ROUTINE: KPP_CALC_DUMMY
a9d2e4c565 Jean*0689
2b4c90c108 Mart*0690 C !INTERFACE: ==========================================================
c417c38d79 Jean*0691       SUBROUTINE KPP_CALC_DUMMY(
                0692      I     bi, bj, myTime, myIter, myThid )
2b4c90c108 Mart*0693
                0694 C !DESCRIPTION: \bv
c417c38d79 Jean*0695 C     *==========================================================*
b7230be7a5 Patr*0696 C     | SUBROUTINE KPP_CALC_DUMMY                                |
                0697 C     | o Compute all KPP fields defined in KPP.h                |
                0698 C     | o Dummy routine for TAMC
c417c38d79 Jean*0699 C     *==========================================================*
b7230be7a5 Patr*0700
2b4c90c108 Mart*0701 C !USES: ===============================================================
                0702       IMPLICIT NONE
b7230be7a5 Patr*0703 #include "SIZE.h"
                0704 #include "EEPARAMS.h"
                0705 #include "PARAMS.h"
                0706 #include "KPP.h"
                0707 #include "KPP_PARAMS.h"
                0708 #include "GRID.h"
059d9fc14f Dimi*0709 #include "GAD.h"
af61e5eb16 Mart*0710 #ifdef ALLOW_AUTODIFF_TAMC
                0711 # include "tamc.h"
                0712 #endif /* ALLOW_AUTODIFF_TAMC */
b7230be7a5 Patr*0713
2b4c90c108 Mart*0714 C !INPUT PARAMETERS: ===================================================
c417c38d79 Jean*0715 c     bi, bj :: Current tile indices
                0716 c     myTime :: Current time in simulation
                0717 c     myIter :: Current iteration number in simulation
                0718 c     myThid :: My Thread Id. number
b7230be7a5 Patr*0719       INTEGER bi, bj
                0720       _RL     myTime
c417c38d79 Jean*0721       INTEGER myIter
                0722       INTEGER myThid
b7230be7a5 Patr*0723
                0724 #ifdef ALLOW_KPP
2b4c90c108 Mart*0725 C !LOCAL VARIABLES: ====================================================
                0726       INTEGER i, j, k
af61e5eb16 Mart*0727 #if ( !defined EXCLUDE_PCELL_MIX_CODE && defined  ALLOW_AUTODIFF_TAMC )
                0728       INTEGER ikey
                0729 #endif
2b4c90c108 Mart*0730 CEOP
b7230be7a5 Patr*0731
                0732       DO j=1-OLy,sNy+OLy
                0733          DO i=1-OLx,sNx+OLx
                0734             KPPhbl (i,j,bi,bj) = 1.0
                0735             KPPfrac(i,j,bi,bj) = 0.0
63ceaaa79c Dimi*0736 #ifdef ALLOW_SALT_PLUME
                0737             KPPplumefrac(i,j,bi,bj) = 0.0
                0738 #endif /* ALLOW_SALT_PLUME */
b7230be7a5 Patr*0739             DO k = 1,Nr
                0740                KPPghat   (i,j,k,bi,bj) = 0.0
a9d2e4c565 Jean*0741                KPPviscAz (i,j,k,bi,bj) = viscArNr(1)
b7230be7a5 Patr*0742             ENDDO
                0743          ENDDO
                0744       ENDDO
059d9fc14f Dimi*0745
af61e5eb16 Mart*0746 #if ( !defined EXCLUDE_PCELL_MIX_CODE && defined  ALLOW_AUTODIFF_TAMC )
                0747       ikey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
                0748 CADJ STORE KPPdiffKzS(:,:,:,bi,bj) = comlev1_kpp, key = ikey
                0749 CADJ STORE KPPdiffKzT(:,:,:,bi,bj) = comlev1_kpp, key = ikey
                0750 #endif
059d9fc14f Dimi*0751       CALL CALC_3D_DIFFUSIVITY(
853c5e0e2c Jean*0752      I     bi,bj,1-OLx,sNx+OLx,1-OLy,sNy+OLy,
059d9fc14f Dimi*0753      I     GAD_SALINITY, .FALSE., .FALSE.,
853c5e0e2c Jean*0754      O     KPPdiffKzS(1-OLx,1-OLy,1,bi,bj),
059d9fc14f Dimi*0755      I     myThid)
                0756       CALL CALC_3D_DIFFUSIVITY(
853c5e0e2c Jean*0757      I     bi,bj,1-OLx,sNx+OLx,1-OLy,sNy+OLy,
059d9fc14f Dimi*0758      I     GAD_TEMPERATURE, .FALSE., .FALSE.,
853c5e0e2c Jean*0759      O     KPPdiffKzT(1-OLx,1-OLy,1,bi,bj),
059d9fc14f Dimi*0760      I     myThid)
                0761
b7230be7a5 Patr*0762 #endif
                0763       RETURN
                0764       END