/[MITgcm]/MITgcm/pkg/kpp/kpp_calc.F
ViewVC logotype

Diff of /MITgcm/pkg/kpp/kpp_calc.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.4 by heimbach, Mon Nov 13 16:37:02 2000 UTC revision 1.42 by jmc, Thu May 3 21:37:34 2007 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "KPP_OPTIONS.h"  #include "KPP_OPTIONS.h"
5    
6        subroutine KPP_CALC(  CBOP
7       I     bi, bj, myTime, myThid )  C !ROUTINE: KPP_CALC
8  C     /==========================================================\  
9    C !INTERFACE: ==========================================================
10          SUBROUTINE KPP_CALC(
11         I     bi, bj, myTime, myIter, myThid )
12    
13    C !DESCRIPTION: \bv
14    C     *==========================================================*
15  C     | SUBROUTINE KPP_CALC                                      |  C     | SUBROUTINE KPP_CALC                                      |
16  C     | o Compute all KPP fields defined in KPP.h                |  C     | o Compute all KPP fields defined in KPP.h                |
17  C     |==========================================================|  C     *==========================================================*
18  C     | This subroutine serves as an interface between MITGCMUV  |  C     | This subroutine serves as an interface between MITGCMUV  |
19  C     | code and NCOM 1-D routines in kpp_routines.F             |  C     | code and NCOM 1-D routines in kpp_routines.F             |
20  C     \==========================================================/  C     *==========================================================*
21        IMPLICIT NONE        IMPLICIT NONE
22    
23  c=======================================================================  c=======================================================================
# Line 45  c         within the boundary layer, abo Line 52  c         within the boundary layer, abo
52  c         determined by turbulent surface fluxes, and interior mixing at  c         determined by turbulent surface fluxes, and interior mixing at
53  c         the lower boundary, i.e. at hbl.  c         the lower boundary, i.e. at hbl.
54  c      c    
55  c     this subroutine provides the interface between the MIT GCM UV and the  c     this subroutine provides the interface between the MITGCM and
56  c     subroutine "kppmix", where boundary layer depth, vertical  c     the routine "kppmix", where boundary layer depth, vertical
57  c     viscosity, vertical diffusivity, and counter gradient term (ghat)  c     viscosity, vertical diffusivity, and counter gradient term (ghat)
58  c     are computed slabwise.  c     are computed slabwise.
59  c     note: subroutine "kppmix" uses m-k-s units.  c     note: subroutine "kppmix" uses m-k-s units.
# Line 83  c     KPPfrac     - Fraction of short-wa Line 90  c     KPPfrac     - Fraction of short-wa
90    
91  c--   KPP_CALC computes vertical viscosity and diffusivity for region  c--   KPP_CALC computes vertical viscosity and diffusivity for region
92  c     (-2:sNx+3,-2:sNy+3) as required by CALC_DIFFUSIVITY and requires  c     (-2:sNx+3,-2:sNy+3) as required by CALC_DIFFUSIVITY and requires
93  c     values of uVel, vVel, SurfaceTendencyU, SurfaceTendencyV in the  c     values of uVel, vVel, surfaceForcingU, surfaceForcingV in the
94  c     region (-2:sNx+4,-2:sNy+4).  c     region (-2:sNx+4,-2:sNy+4).
95  c     Hence overlap region needs to be set OLx=4, OLy=4.  c     Hence overlap region needs to be set OLx=4, OLy=4.
96  c     When option FRUGAL_KPP is used, computation in overlap regions  c \ev
 c     is replaced with exchange calls hence reducing overlap requirements  
 c     to OLx=1, OLy=1.  
97    
98    C !USES: ===============================================================
99  #include "SIZE.h"  #include "SIZE.h"
100  #include "EEPARAMS.h"  #include "EEPARAMS.h"
101  #include "PARAMS.h"  #include "PARAMS.h"
# Line 98  c     to OLx=1, OLy=1. Line 104  c     to OLx=1, OLy=1.
104  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
105  #include "FFIELDS.h"  #include "FFIELDS.h"
106  #include "GRID.h"  #include "GRID.h"
107    #include "GAD.h"
108    #ifdef ALLOW_SHELFICE
109    # include "SHELFICE.h"
110    #endif /* ALLOW_SHELFICE */
111  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
112  #include "tamc.h"  #include "tamc.h"
113  #include "tamc_keys.h"  #include "tamc_keys.h"
       INTEGER    isbyte  
       PARAMETER( isbyte = 4 )  
114  #else /* ALLOW_AUTODIFF_TAMC */  #else /* ALLOW_AUTODIFF_TAMC */
115        integer ikey        integer ikppkey
116  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
117    
118        EXTERNAL DIFFERENT_MULTIPLE        EXTERNAL DIFFERENT_MULTIPLE
119        LOGICAL  DIFFERENT_MULTIPLE        LOGICAL  DIFFERENT_MULTIPLE
120    
121    C !INPUT PARAMETERS: ===================================================
122  c Routine arguments  c Routine arguments
123  c     bi, bj - array indices on which to apply calculations  c     bi, bj :: Current tile indices
124  c     myTime - Current time in simulation  c     myTime :: Current time in simulation
125    c     myIter :: Current iteration number in simulation
126    c     myThid :: My Thread Id. number
127    
128        INTEGER bi, bj        INTEGER bi, bj
       INTEGER myThid  
129        _RL     myTime        _RL     myTime
130          INTEGER myIter
131          INTEGER myThid
132    
133  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
134    
135    C !LOCAL VARIABLES: ====================================================
136  c Local constants  c Local constants
137  c     minusone, p0, p5, p25, p125, p0625  c     minusone, p0, p5, p25, p125, p0625
138  c     imin, imax, jmin, jmax  - array computation indices  c     imin, imax, jmin, jmax  - array computation indices
139    
140        _RL        minusone        _RL        minusone
141        parameter( minusone=-1.0)        parameter( minusone=-1.0)
142        _KPP_RL    p0    , p5    , p25     , p125      , p0625        _RL        p0    , p5    , p25     , p125      , p0625
143        parameter( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 )        parameter( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 )
144        integer    imin      , imax        , jmin      , jmax        integer   imin      ,imax          ,jmin      ,jmax
145  #ifdef FRUGAL_KPP        parameter(imin=2-OLx,imax=sNx+OLx-1,jmin=2-OLy,jmax=sNy+OLy-1)
       parameter( imin=1    , imax=sNx    , jmin=1    , jmax=sNy     )  
 #else  
       parameter( imin=-2   , imax=sNx+3  , jmin=-2   , jmax=sNy+3   )  
 #endif  
146    
147  c Local arrays and variables  c Local arrays and variables
148  c     work?  (nx,ny)       - horizontal working arrays  c     work?  (nx,ny)       - horizontal working arrays
# Line 150  c                            for ri_iwmi Line 158  c                            for ri_iwmi
158  c     Ritop  (nx,ny,Nr)    - numerator of bulk richardson number  c     Ritop  (nx,ny,Nr)    - numerator of bulk richardson number
159  c                            at grid levels for bldepth  c                            at grid levels for bldepth
160  c     vddiff (nx,ny,Nrp2,1)- vertical viscosity on "t-grid"           (m^2/s)  c     vddiff (nx,ny,Nrp2,1)- vertical viscosity on "t-grid"           (m^2/s)
161  c     vddiff (nx,ny,Nrp2,2)- vert. diff. on next row for temperature  (m^2/s)  c     vddiff (nx,ny,Nrp2,2)- vert. diff. on next row for salt&tracers (m^2/s)
162  c     vddiff (nx,ny,Nrp2,3)- vert. diff. on next row for salt&tracers (m^2/s)  c     vddiff (nx,ny,Nrp2,3)- vert. diff. on next row for temperature  (m^2/s)
163  c     ghat   (nx,ny,Nr)    - nonlocal transport coefficient           (s/m^2)  c     ghat   (nx,ny,Nr)    - nonlocal transport coefficient           (s/m^2)
164  c     hbl    (nx,ny)       - mixing layer depth                           (m)  c     hbl    (nx,ny)       - mixing layer depth                           (m)
165  c     kmtj   (nx,ny)       - maximum number of wet levels in each column  c     kmtj   (nx,ny)       - maximum number of wet levels in each column
# Line 160  c     zRef   (nx,ny)       - Reference d Line 168  c     zRef   (nx,ny)       - Reference d
168  c     uRef   (nx,ny)       - Reference zonal velocity                   (m/s)  c     uRef   (nx,ny)       - Reference zonal velocity                   (m/s)
169  c     vRef   (nx,ny)       - Reference meridional velocity              (m/s)  c     vRef   (nx,ny)       - Reference meridional velocity              (m/s)
170    
171        _RL     worka ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )        integer work1 ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy            )
172        integer work1 ( ibot:itop    , jbot:jtop                    )        _RL worka ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
173        _KPP_RL work2 ( ibot:itop    , jbot:jtop                    )        _RL work2 ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
174        _KPP_RL ustar ( ibot:itop    , jbot:jtop                    )        _RL work3 ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
175        _KPP_RL bo    ( ibot:itop    , jbot:jtop                    )        _RL ustar ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
176        _KPP_RL bosol ( ibot:itop    , jbot:jtop                    )        _RL bo    ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
177        _KPP_RL shsq  ( ibot:itop    , jbot:jtop    , Nr            )        _RL bosol ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
178        _KPP_RL dVsq  ( ibot:itop    , jbot:jtop    , Nr            )        _RL shsq  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
179        _KPP_RL dbloc ( ibot:itop    , jbot:jtop    , Nr            )        _RL dVsq  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
180        _KPP_RL Ritop ( ibot:itop    , jbot:jtop    , Nr            )        _RL dbloc ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
181        _KPP_RL vddiff( ibot:itop    , jbot:jtop    , 0:Nrp1, mdiff )        _RL Ritop ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
182        _KPP_RL ghat  ( ibot:itop    , jbot:jtop    , Nr            )        _RL vddiff( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, 0:Nrp1, mdiff )
183        _KPP_RL hbl   ( ibot:itop    , jbot:jtop                    )        _RL ghat  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
184          _RL hbl   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
185    cph(
186          _RL TTALPHA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
187          _RL SSBETA ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
188    cph)
189  #ifdef KPP_ESTIMATE_UREF  #ifdef KPP_ESTIMATE_UREF
190        _KPP_RL z0    ( ibot:itop    , jbot:jtop                    )        _RL z0    ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
191        _KPP_RL zRef  ( ibot:itop    , jbot:jtop                    )        _RL zRef  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
192        _KPP_RL uRef  ( ibot:itop    , jbot:jtop                    )        _RL uRef  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
193        _KPP_RL vRef  ( ibot:itop    , jbot:jtop                    )        _RL vRef  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
194  #endif /* KPP_ESTIMATE_UREF */  #endif /* KPP_ESTIMATE_UREF */
195          
196        _KPP_RL tempvar1, tempvar2        _RL tempvar2
197        integer i, j, k, kp1, im1, ip1, jm1, jp1        integer i, j, k, kp1, km1, im1, ip1, jm1, jp1
198    
199  #ifdef KPP_ESTIMATE_UREF  #ifdef KPP_ESTIMATE_UREF
200        _KPP_RL dBdz1, dBdz2, ustarX, ustarY        _RL tempvar1, dBdz1, dBdz2, ustarX, ustarY
201  #endif  #endif
202    
203    #ifdef ALLOW_AUTODIFF_TAMC
204              act1 = bi - myBxLo(myThid)
205              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
206              act2 = bj - myByLo(myThid)
207              max2 = myByHi(myThid) - myByLo(myThid) + 1
208              act3 = myThid - 1
209              max3 = nTx*nTy
210              act4 = ikey_dynamics - 1
211              ikppkey = (act1 + 1) + act2*max1
212         &                      + act3*max1*max2
213         &                      + act4*max1*max2*max3
214    #endif /* ALLOW_AUTODIFF_TAMC */
215    CEOP
216    
217  c     Check to see if new vertical mixing coefficient should be computed now?  c     Check to see if new vertical mixing coefficient should be computed now?
218        IF ( DIFFERENT_MULTIPLE(kpp_freq,myTime,myTime-deltaTClock) .OR.        IF ( DIFFERENT_MULTIPLE(kpp_freq,myTime,deltaTClock)
219       1     myTime .EQ. startTime ) THEN       1     .OR. myTime .EQ. startTime ) THEN
220            
221  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
222  c     prepare input arrays for subroutine "kppmix" to compute  c     prepare input arrays for subroutine "kppmix" to compute
223  c     viscosity and diffusivity and ghat.  c     viscosity and diffusivity and ghat.
# Line 227  c--------------------------------------- Line 254  c---------------------------------------
254    
255        CALL TIMER_START('STATEKPP      [KPP_CALC]', myThid)        CALL TIMER_START('STATEKPP      [KPP_CALC]', myThid)
256        CALL STATEKPP(        CALL STATEKPP(
257       I       bi, bj, myThid       O     work2, dbloc, Ritop,
258       O     , work2, dbloc, Ritop       O     TTALPHA, SSBETA,
259       O     , vddiff(ibot,jbot,1,1), vddiff(ibot,jbot,1,2)       I     ikppkey, bi, bj, myThid )
      &     )  
260        CALL TIMER_STOP ('STATEKPP      [KPP_CALC]', myThid)        CALL TIMER_STOP ('STATEKPP      [KPP_CALC]', myThid)
261    
262        DO k = 1, Nr        DO k = 1, Nr
263           DO j = jbot, jtop           DO j = 1-OLy, sNy+OLy
264              DO i = ibot, itop              DO i = 1-OLx, sNx+OLx
265                 ghat(i,j,k) = dbloc(i,j,k)                 ghat(i,j,k) = dbloc(i,j,k)
266              ENDDO              ENDDO
267           ENDDO           ENDDO
# Line 248  c     dbloc(k) is buoyancy gradientnote Line 274  c     dbloc(k) is buoyancy gradientnote
274  c     levels therefore k+1 mask must be used  c     levels therefore k+1 mask must be used
275    
276        DO k = 1, Nr-1        DO k = 1, Nr-1
277           CALL KPP_SMOOTH_HORIZ (           CALL SMOOTH_HORIZ (
278       I        k+1, bi, bj,       I        k+1, bi, bj,
279       U        ghat (ibot,jbot,k) )       U        ghat (1-OLx,1-OLy,k),
280         I        myThid )
281        ENDDO        ENDDO
282    
283  #endif /* KPP_SMOOTH_DBLOC */  #endif /* KPP_SMOOTH_DBLOC */
284    
285  #ifdef KPP_SMOOTH_DENS  #ifdef KPP_SMOOTH_DENS
286  c     horizontally smooth density related quantities with 121 filters  c     horizontally smooth density related quantities with 121 filters
287        CALL KPP_SMOOTH_HORIZ (        CALL SMOOTH_HORIZ (
288       I     1, bi, bj,       I     1, bi, bj,
289       U     work2 )       U     work2,
290         I     myThid )
291        DO k = 1, Nr        DO k = 1, Nr
292           CALL KPP_SMOOTH_HORIZ (           CALL SMOOTH_HORIZ (
293       I        k+1, bi, bj,       I        k+1, bi, bj,
294       U        dbloc (ibot,jbot,k) )       U        dbloc (1-OLx,1-OLy,k),
295           CALL KPP_SMOOTH_HORIZ (       I        myThid )
296             CALL SMOOTH_HORIZ (
297       I        k, bi, bj,       I        k, bi, bj,
298       U        Ritop (ibot,jbot,k)  )       U        Ritop (1-OLx,1-OLy,k),
299           CALL KPP_SMOOTH_HORIZ (       I        myThid )
300             CALL SMOOTH_HORIZ (
301       I        k, bi, bj,       I        k, bi, bj,
302       U        vddiff(ibot,jbot,k,1) )       U        TTALPHA(1-OLx,1-OLy,k),
303           CALL KPP_SMOOTH_HORIZ (       I        myThid )
304             CALL SMOOTH_HORIZ (
305       I        k, bi, bj,       I        k, bi, bj,
306       U        vddiff(ibot,jbot,k,2) )       U        SSBETA(1-OLx,1-OLy,k),
307         I        myThid )
308        ENDDO        ENDDO
309  #endif /* KPP_SMOOTH_DENS */  #endif /* KPP_SMOOTH_DENS */
310    
311        DO k = 1, Nr        DO k = 1, Nr
312           DO j = jbot, jtop           km1 = max(1,k-1)
313              DO i = ibot, itop           DO j = 1-OLy, sNy+OLy
314                DO i = 1-OLx, sNx+OLx
315    
316  c     zero out dbloc over land points (so that the convective  c     zero out dbloc over land points (so that the convective
317  c     part of the interior mixing can be diagnosed)  c     part of the interior mixing can be diagnosed)
318                 dbloc(i,j,k) = dbloc(i,j,k) * pMask(i,j,k,bi,bj)                 dbloc(i,j,k) = dbloc(i,j,k) * maskC(i,j,k,bi,bj)
319                 ghat(i,j,k)  = ghat(i,j,k)  * pMask(i,j,k,bi,bj)       &              * maskC(i,j,km1,bi,bj)
320                 Ritop(i,j,k) = Ritop(i,j,k) * pMask(i,j,k,bi,bj)                 ghat(i,j,k)  = ghat(i,j,k)  * maskC(i,j,k,bi,bj)
321         &              * maskC(i,j,km1,bi,bj)
322                   Ritop(i,j,k) = Ritop(i,j,k) * maskC(i,j,k,bi,bj)
323         &              * maskC(i,j,km1,bi,bj)
324                 if(k.eq.nzmax(i,j,bi,bj)) then                 if(k.eq.nzmax(i,j,bi,bj)) then
325                    dbloc(i,j,k) = p0                    dbloc(i,j,k) = p0
326                    ghat(i,j,k)  = p0                    ghat(i,j,k)  = p0
# Line 296  c     note: land and ocean bottom values Line 332  c     note: land and ocean bottom values
332  c     so that the subroutine "bldepth" works correctly  c     so that the subroutine "bldepth" works correctly
333                 Ritop(i,j,k) = (zgrid(1)-zgrid(k)) * Ritop(i,j,k)                 Ritop(i,j,k) = (zgrid(1)-zgrid(k)) * Ritop(i,j,k)
334    
335              END DO              ENDDO
336           END DO           ENDDO
337        END DO        ENDDO
338    
339    cph(
340    cph  this avoids a single or double recomp./call of statekpp
341    CADJ store work2              = comlev1_kpp, key = ikppkey
342    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
343    CADJ store dbloc, Ritop, ghat = comlev1_kpp, key = ikppkey
344    CADJ store vddiff             = comlev1_kpp, key = ikppkey
345    CADJ store TTALPHA, SSBETA    = comlev1_kpp, key = ikppkey
346    #endif
347    cph)
348    
349    CML#ifdef ALLOW_SHELFICE
350    CMLC     For the pbl parameterisation to work underneath the ice shelves
351    CMLC     it needs to know the surface (ice-ocean) fluxes. However, masking
352    CMLC     and indexing problems make this part of the code not work
353    CMLC     underneath the ice shelves and the following lines are only here
354    CMLC     to remind me that this still needs to be sorted out.
355    CML      shelfIceFac = 0. _d 0
356    CML      IF ( useShelfIce ) selfIceFac = 1. _d 0
357    CML      DO j = jmin, jmax
358    CML       DO i = imin, imax
359    CML        surfForcT = surfaceForcingT(i,j,bi,bj)
360    CML     &       + shelficeForcingT(i,j,bi,bj) * shelfIceFac
361    CML        surfForcS = surfaceForcingS(i,j,bi,bj)
362    CML     &       + shelficeForcingS(i,j,bi,bj) * shelfIceFac
363    CML       ENDDO
364    CML      ENDDO
365    CML#endif /* ALLOW_SHELFICE */
366    
367  c------------------------------------------------------------------------  c------------------------------------------------------------------------
368  c     friction velocity, turbulent and radiative surface buoyancy forcing  c     friction velocity, turbulent and radiative surface buoyancy forcing
369  c     -------------------------------------------------------------------  c     -------------------------------------------------------------------
370  c     taux / rho = SurfaceTendencyU * delZ(1)                     (N/m^2)  c     taux / rho = surfaceForcingU                               (N/m^2)
371  c     tauy / rho = SurfaceTendencyV * delZ(1)                     (N/m^2)  c     tauy / rho = surfaceForcingV                               (N/m^2)
372  c     ustar = sqrt( sqrt( taux^2 + tauy^2 ) / rho )                 (m/s)  c     ustar = sqrt( sqrt( taux^2 + tauy^2 ) / rho )                (m/s)
373  c     bo    = - g * ( alpha*SurfaceTendencyT +  c     bo    = - g * ( alpha*surfaceForcingT +
374  c                     beta *SurfaceTendencyS ) * delZ(1) / rho  (m^2/s^3)  c                     beta *surfaceForcingS ) / rho            (m^2/s^3)
375  c     bosol = - g * alpha * Qsw * delZ(1) / rho                 (m^2/s^3)  c     bosol = - g * alpha * Qsw * drF(1) / rho                 (m^2/s^3)
 c------------------------------------------------------------------------  
   
 c initialize arrays to zero  
       DO j = jbot, jtop  
          DO i = ibot, itop  
             ustar(i,j) = p0  
             bo   (I,J) = p0  
             bosol(I,J) = p0  
          END DO  
       END DO  
   
       DO j = jmin, jmax  
        jp1 = j + 1  
        DO i = imin, imax  
         ip1 = i+1  
         tempVar1 =  
      &   (SurfaceTendencyU(i,j,bi,bj) + SurfaceTendencyU(ip1,j,bi,bj)) *  
      &   (SurfaceTendencyU(i,j,bi,bj) + SurfaceTendencyU(ip1,j,bi,bj)) +  
      &   (SurfaceTendencyV(i,j,bi,bj) + SurfaceTendencyV(i,jp1,bi,bj)) *  
      &   (SurfaceTendencyV(i,j,bi,bj) + SurfaceTendencyV(i,jp1,bi,bj))  
         if ( tempVar1 .lt. (phepsi*phepsi) ) then  
            ustar(i,j) = SQRT( phepsi * p5 * delZ(1) )  
         else  
            tempVar2 =  SQRT( tempVar1 ) * p5 * delZ(1)  
            ustar(i,j) = SQRT( tempVar2 )  
         endif  
         bo(I,J) = - gravity *  
      &       ( vddiff(I,J,1,1) * SurfaceTendencyT(i,j,bi,bj) +  
      &         vddiff(I,J,1,2) * SurfaceTendencyS(i,j,bi,bj)  
      &       ) *  
      &       delZ(1) / work2(I,J)  
         bosol(I,J) = - gravity * vddiff(I,J,1,1) * Qsw(i,j,bi,bj) *  
      &       delZ(1) / work2(I,J)  
        END DO  
       END DO  
   
376  c------------------------------------------------------------------------  c------------------------------------------------------------------------
377  c     velocity shear  c     velocity shear
378  c     --------------  c     --------------
# Line 352  c     Get velocity shear squared, averag Line 380  c     Get velocity shear squared, averag
380  c     onto "t-grid" (in (m/s)**2):  c     onto "t-grid" (in (m/s)**2):
381  c     dVsq(k)=(Uref-U(k))**2+(Vref-V(k))**2      at grid levels  c     dVsq(k)=(Uref-U(k))**2+(Vref-V(k))**2      at grid levels
382  c     shsq(k)=(U(k)-U(k+1))**2+(V(k)-V(k+1))**2  at interfaces  c     shsq(k)=(U(k)-U(k+1))**2+(V(k)-V(k+1))**2  at interfaces
383    c    
384    c     note: Vref can depend on the surface fluxes that is why we compute
385    c     dVsq in the subroutine that does the surface related stuff
386    c     (admittedly this is a bit messy)
387  c------------------------------------------------------------------------  c------------------------------------------------------------------------
388    
389  c initialize arrays to zero        CALL KPP_FORCING_SURF(
390        DO k = 1, Nr       I     work2, surfaceForcingU, surfaceForcingV,
391           DO j = jbot, jtop       I     surfaceForcingT, surfaceForcingS, surfaceForcingTice,
392              DO i = ibot, itop       I     Qsw, ttalpha, ssbeta,
393                 shsq(i,j,k) = p0       O     ustar, bo, bosol, dVsq,
394                 dVsq(i,j,k) = p0       I     ikppkey, iMin, iMax, jMin, jMax, bi, bj, myTime, myThid )
395              END DO  
396           END DO  CMLcph(
397        END DO  CMLCADJ store ustar = comlev1_kpp, key = ikppkey
398    CMLcph)
 c     dVsq computation  
   
 #ifdef KPP_ESTIMATE_UREF  
   
 c     Get rid of vertical resolution dependence of dVsq term by  
 c     estimating a surface velocity that is independent of first level  
 c     thickness in the model.  First determine mixed layer depth hMix.  
 c     Second zRef = espilon * hMix.  Third determine roughness length  
 c     scale z0.  Third estimate reference velocity.  
   
       DO j = jmin, jmax  
          jp1 = j + 1  
          DO i = imin, imax  
             ip1 = i + 1  
   
 c     Determine mixed layer depth hMix as the shallowest depth at which  
 c     dB/dz exceeds 5.2e-5 s^-2.  
             work1(i,j) = nzmax(i,j,bi,bj)  
             DO k = 1, Nr  
                IF ( k .LT. nzmax(i,j,bi,bj) .AND.  
      &              dbloc(i,j,k) / drC(k+1) .GT. dB_dz )  
      &              work1(i,j) = k  
             END DO  
   
 c     Linearly interpolate to find hMix.  
             k = work1(i,j)  
             IF ( k .EQ. 0 .OR. nzmax(i,j,bi,bj) .EQ. 1 ) THEN  
                zRef(i,j) = p0  
             ELSEIF ( k .EQ. 1) THEN  
                dBdz2 = dbloc(i,j,1) / drC(2)  
                zRef(i,j) = drF(1) * dB_dz / dBdz2  
             ELSEIF ( k .LT. nzmax(i,j,bi,bj) ) THEN  
                dBdz1 = dbloc(i,j,k-1) / drC(k  )  
                dBdz2 = dbloc(i,j,k  ) / drC(k+1)  
                zRef(i,j) = rF(k) + drF(k) * (dB_dz - dBdz1) /  
      &                     MAX ( phepsi, dBdz2 - dBdz1 )  
             ELSE  
                zRef(i,j) = rF(k+1)  
             ENDIF  
   
 c     Compute roughness length scale z0 subject to 0 < z0  
                tempVar1 = p5 * (  
      &              (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  2,bi,bj)) *  
      &              (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  2,bi,bj)) +  
      &              (uVel(ip1,j,  1,bi,bj)-uVel(ip1,j,  2,bi,bj)) *  
      &              (uVel(ip1,j,  1,bi,bj)-uVel(ip1,j,  2,bi,bj)) +  
      &              (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  2,bi,bj)) *  
      &              (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  2,bi,bj)) +  
      &              (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,2,bi,bj)) *  
      &              (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,2,bi,bj)) )  
                if ( tempVar1 .lt. (epsln*epsln) ) then  
                   tempVar2 = epsln  
                else  
                   tempVar2 = SQRT ( tempVar1 )  
                endif  
                z0(i,j) = rF(2) *  
      &                   ( rF(3) * LOG ( rF(3) / rF(2) ) /  
      &                     ( rF(3) - rF(2) ) -  
      &                     tempVar2 * vonK /  
      &                     MAX ( ustar(i,j), phepsi ) )  
                z0(i,j) = MAX ( z0(i,j), phepsi )  
   
 c     zRef is set to 0.1 * hMix subject to z0 <= zRef <= drF(1)  
                zRef(i,j) = MAX ( epsilon * zRef(i,j), z0(i,j) )  
                zRef(i,j) = MIN ( zRef(i,j), drF(1) )  
   
 c     Estimate reference velocity uRef and vRef.  
                uRef(i,j) = p5 *  
      &                     ( uVel(i,j,1,bi,bj) + uVel(ip1,j,1,bi,bj) )  
                vRef(i,j) = p5 *  
      &                     ( vVel(i,j,1,bi,bj) + vVel(i,jp1,1,bi,bj) )  
                IF ( zRef(i,j) .LT. drF(1) ) THEN  
                   ustarX = ( SurfaceTendencyU(i,  j,bi,bj) +  
      &                       SurfaceTendencyU(ip1,j,bi,bj) ) * p5  
                   ustarY = ( SurfaceTendencyV(i,j,  bi,bj) +  
      &                       SurfaceTendencyU(i,jp1,bi,bj) ) * p5  
                   tempVar1 = ustarX * ustarX + ustarY * ustarY  
                   if ( tempVar1 .lt. (epsln*epsln) ) then  
                      tempVar2 = epsln  
                   else  
                      tempVar2 = SQRT ( tempVar1 )  
                   endif  
                   tempVar2 = ustar(i,j) *  
      &                 ( LOG ( zRef(i,j) / rF(2) ) +  
      &                 z0(i,j) / zRef(i,j) - z0(i,j) / rF(2) ) /  
      &                 vonK / tempVar2  
                   uRef(i,j) = uRef(i,j) + ustarX * tempVar2  
                   vRef(i,j) = vRef(i,j) + ustarY * tempVar2  
                ENDIF  
   
          END DO  
       END DO  
   
       DO k = 1, Nr  
          DO j = jmin, jmax  
             jm1 = j - 1  
             jp1 = j + 1  
             DO i = imin, imax  
                im1 = i - 1  
                ip1 = i + 1  
                dVsq(i,j,k) = p5 * (  
      $              (uRef(i,j) - uVel(i,  j,  k,bi,bj)) *  
      $              (uRef(i,j) - uVel(i,  j,  k,bi,bj)) +  
      $              (uRef(i,j) - uVel(ip1,j,  k,bi,bj)) *  
      $              (uRef(i,j) - uVel(ip1,j,  k,bi,bj)) +  
      $              (vRef(i,j) - vVel(i,  j,  k,bi,bj)) *  
      $              (vRef(i,j) - vVel(i,  j,  k,bi,bj)) +  
      $              (vRef(i,j) - vVel(i,  jp1,k,bi,bj)) *  
      $              (vRef(i,j) - vVel(i,  jp1,k,bi,bj)) )  
 #ifdef KPP_SMOOTH_DVSQ  
                dVsq(i,j,k) = p5 * dVsq(i,j,k) + p125 * (  
      $              (uRef(i,j) - uVel(i,  jm1,k,bi,bj)) *  
      $              (uRef(i,j) - uVel(i,  jm1,k,bi,bj)) +  
      $              (uRef(i,j) - uVel(ip1,jm1,k,bi,bj)) *  
      $              (uRef(i,j) - uVel(ip1,jm1,k,bi,bj)) +  
      $              (uRef(i,j) - uVel(i,  jp1,k,bi,bj)) *  
      $              (uRef(i,j) - uVel(i,  jp1,k,bi,bj)) +  
      $              (uRef(i,j) - uVel(ip1,jp1,k,bi,bj)) *  
      $              (uRef(i,j) - uVel(ip1,jp1,k,bi,bj)) +  
      $              (vRef(i,j) - vVel(im1,j,  k,bi,bj)) *  
      $              (vRef(i,j) - vVel(im1,j,  k,bi,bj)) +  
      $              (vRef(i,j) - vVel(im1,jp1,k,bi,bj)) *  
      $              (vRef(i,j) - vVel(im1,jp1,k,bi,bj)) +  
      $              (vRef(i,j) - vVel(ip1,j,  k,bi,bj)) *  
      $              (vRef(i,j) - vVel(ip1,j,  k,bi,bj)) +  
      $              (vRef(i,j) - vVel(ip1,jp1,k,bi,bj)) *  
      $              (vRef(i,j) - vVel(ip1,jp1,k,bi,bj)) )  
 #endif /* KPP_SMOOTH_DVSQ */  
             END DO  
          END DO  
       END DO  
   
 #else /* KPP_ESTIMATE_UREF */  
399    
400    c initialize arrays to zero
401        DO k = 1, Nr        DO k = 1, Nr
402           DO j = jmin, jmax         DO j = 1-OLy, sNy+OLy
403              jm1 = j - 1          DO i = 1-OLx, sNx+OLx
404              jp1 = j + 1           shsq(i,j,k) = p0
405              DO i = imin, imax          ENDDO
406                 im1 = i - 1         ENDDO
407                 ip1 = i + 1        ENDDO
                dVsq(i,j,k) = p5 * (  
      $              (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  k,bi,bj)) *  
      $              (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  k,bi,bj)) +  
      $              (uVel(ip1,j,  1,bi,bj)-uVel(ip1,j,  k,bi,bj)) *  
      $              (uVel(ip1,j,  1,bi,bj)-uVel(ip1,j,  k,bi,bj)) +  
      $              (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  k,bi,bj)) *  
      $              (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  k,bi,bj)) +  
      $              (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,k,bi,bj)) *  
      $              (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,k,bi,bj)) )  
 #ifdef KPP_SMOOTH_DVSQ  
                dVsq(i,j,k) = p5 * dVsq(i,j,k) + p125 * (  
      $              (uVel(i,  jm1,1,bi,bj)-uVel(i,  jm1,k,bi,bj)) *  
      $              (uVel(i,  jm1,1,bi,bj)-uVel(i,  jm1,k,bi,bj)) +  
      $              (uVel(ip1,jm1,1,bi,bj)-uVel(ip1,jm1,k,bi,bj)) *  
      $              (uVel(ip1,jm1,1,bi,bj)-uVel(ip1,jm1,k,bi,bj)) +  
      $              (uVel(i,  jp1,1,bi,bj)-uVel(i,  jp1,k,bi,bj)) *  
      $              (uVel(i,  jp1,1,bi,bj)-uVel(i,  jp1,k,bi,bj)) +  
      $              (uVel(ip1,jp1,1,bi,bj)-uVel(ip1,jp1,k,bi,bj)) *  
      $              (uVel(ip1,jp1,1,bi,bj)-uVel(ip1,jp1,k,bi,bj)) +  
      $              (vVel(im1,j,  1,bi,bj)-vVel(im1,j,  k,bi,bj)) *  
      $              (vVel(im1,j,  1,bi,bj)-vVel(im1,j,  k,bi,bj)) +  
      $              (vVel(im1,jp1,1,bi,bj)-vVel(im1,jp1,k,bi,bj)) *  
      $              (vVel(im1,jp1,1,bi,bj)-vVel(im1,jp1,k,bi,bj)) +  
      $              (vVel(ip1,j,  1,bi,bj)-vVel(ip1,j,  k,bi,bj)) *  
      $              (vVel(ip1,j,  1,bi,bj)-vVel(ip1,j,  k,bi,bj)) +  
      $              (vVel(ip1,jp1,1,bi,bj)-vVel(ip1,jp1,k,bi,bj)) *  
      $              (vVel(ip1,jp1,1,bi,bj)-vVel(ip1,jp1,k,bi,bj)) )  
 #endif /* KPP_SMOOTH_DVSQ */  
             END DO  
          END DO  
       END DO  
   
 #endif /* KPP_ESTIMATE_UREF */  
408    
409  c     shsq computation  c     shsq computation
410        DO k = 1, Nrm1        DO k = 1, Nrm1
411           kp1 = k + 1         kp1 = k + 1
412           DO j = jmin, jmax         DO j = jmin, jmax
413              jm1 = j - 1          jm1 = j - 1
414              jp1 = j + 1          jp1 = j + 1
415              DO i = imin, imax          DO i = imin, imax
416                 im1 = i - 1           im1 = i - 1
417                 ip1 = i + 1           ip1 = i + 1
418                 shsq(i,j,k) = p5 * (           shsq(i,j,k) = p5 * (
419       $              (uVel(i,  j,  k,bi,bj)-uVel(i,  j,  kp1,bi,bj)) *       &        (uVel(i,  j,  k,bi,bj)-uVel(i,  j,  kp1,bi,bj)) *
420       $              (uVel(i,  j,  k,bi,bj)-uVel(i,  j,  kp1,bi,bj)) +       &        (uVel(i,  j,  k,bi,bj)-uVel(i,  j,  kp1,bi,bj)) +
421       $              (uVel(ip1,j,  k,bi,bj)-uVel(ip1,j,  kp1,bi,bj)) *       &        (uVel(ip1,j,  k,bi,bj)-uVel(ip1,j,  kp1,bi,bj)) *
422       $              (uVel(ip1,j,  k,bi,bj)-uVel(ip1,j,  kp1,bi,bj)) +       &        (uVel(ip1,j,  k,bi,bj)-uVel(ip1,j,  kp1,bi,bj)) +
423       $              (vVel(i,  j,  k,bi,bj)-vVel(i,  j,  kp1,bi,bj)) *       &        (vVel(i,  j,  k,bi,bj)-vVel(i,  j,  kp1,bi,bj)) *
424       $              (vVel(i,  j,  k,bi,bj)-vVel(i,  j,  kp1,bi,bj)) +       &        (vVel(i,  j,  k,bi,bj)-vVel(i,  j,  kp1,bi,bj)) +
425       $              (vVel(i,  jp1,k,bi,bj)-vVel(i,  jp1,kp1,bi,bj)) *       &        (vVel(i,  jp1,k,bi,bj)-vVel(i,  jp1,kp1,bi,bj)) *
426       $              (vVel(i,  jp1,k,bi,bj)-vVel(i,  jp1,kp1,bi,bj)) )       &        (vVel(i,  jp1,k,bi,bj)-vVel(i,  jp1,kp1,bi,bj)) )
427  #ifdef KPP_SMOOTH_SHSQ  #ifdef KPP_SMOOTH_SHSQ
428                 shsq(i,j,k) = p5 * shsq(i,j,k) + p125 * (           shsq(i,j,k) = p5 * shsq(i,j,k) + p125 * (
429       $              (uVel(i,  jm1,k,bi,bj)-uVel(i,  jm1,kp1,bi,bj)) *       &        (uVel(i,  jm1,k,bi,bj)-uVel(i,  jm1,kp1,bi,bj)) *
430       $              (uVel(i,  jm1,k,bi,bj)-uVel(i,  jm1,kp1,bi,bj)) +       &        (uVel(i,  jm1,k,bi,bj)-uVel(i,  jm1,kp1,bi,bj)) +
431       $              (uVel(ip1,jm1,k,bi,bj)-uVel(ip1,jm1,kp1,bi,bj)) *       &        (uVel(ip1,jm1,k,bi,bj)-uVel(ip1,jm1,kp1,bi,bj)) *
432       $              (uVel(ip1,jm1,k,bi,bj)-uVel(ip1,jm1,kp1,bi,bj)) +       &        (uVel(ip1,jm1,k,bi,bj)-uVel(ip1,jm1,kp1,bi,bj)) +
433       $              (uVel(i,  jp1,k,bi,bj)-uVel(i,  jp1,kp1,bi,bj)) *       &        (uVel(i,  jp1,k,bi,bj)-uVel(i,  jp1,kp1,bi,bj)) *
434       $              (uVel(i,  jp1,k,bi,bj)-uVel(i,  jp1,kp1,bi,bj)) +       &        (uVel(i,  jp1,k,bi,bj)-uVel(i,  jp1,kp1,bi,bj)) +
435       $              (uVel(ip1,jp1,k,bi,bj)-uVel(ip1,jp1,kp1,bi,bj)) *       &        (uVel(ip1,jp1,k,bi,bj)-uVel(ip1,jp1,kp1,bi,bj)) *
436       $              (uVel(ip1,jp1,k,bi,bj)-uVel(ip1,jp1,kp1,bi,bj)) +       &        (uVel(ip1,jp1,k,bi,bj)-uVel(ip1,jp1,kp1,bi,bj)) +
437       $              (vVel(im1,j,  k,bi,bj)-vVel(im1,j,  kp1,bi,bj)) *       &        (vVel(im1,j,  k,bi,bj)-vVel(im1,j,  kp1,bi,bj)) *
438       $              (vVel(im1,j,  k,bi,bj)-vVel(im1,j,  kp1,bi,bj)) +       &        (vVel(im1,j,  k,bi,bj)-vVel(im1,j,  kp1,bi,bj)) +
439       $              (vVel(im1,jp1,k,bi,bj)-vVel(im1,jp1,kp1,bi,bj)) *       &        (vVel(im1,jp1,k,bi,bj)-vVel(im1,jp1,kp1,bi,bj)) *
440       $              (vVel(im1,jp1,k,bi,bj)-vVel(im1,jp1,kp1,bi,bj)) +       &        (vVel(im1,jp1,k,bi,bj)-vVel(im1,jp1,kp1,bi,bj)) +
441       $              (vVel(ip1,j,  k,bi,bj)-vVel(ip1,j,  kp1,bi,bj)) *       &        (vVel(ip1,j,  k,bi,bj)-vVel(ip1,j,  kp1,bi,bj)) *
442       $              (vVel(ip1,j,  k,bi,bj)-vVel(ip1,j,  kp1,bi,bj)) +       &        (vVel(ip1,j,  k,bi,bj)-vVel(ip1,j,  kp1,bi,bj)) +
443       $              (vVel(ip1,jp1,k,bi,bj)-vVel(ip1,jp1,kp1,bi,bj)) *       &        (vVel(ip1,jp1,k,bi,bj)-vVel(ip1,jp1,kp1,bi,bj)) *
444       $              (vVel(ip1,jp1,k,bi,bj)-vVel(ip1,jp1,kp1,bi,bj)) )       &        (vVel(ip1,jp1,k,bi,bj)-vVel(ip1,jp1,kp1,bi,bj)) )
445  #endif  #endif
446              END DO          ENDDO
447           END DO         ENDDO
448        END DO        ENDDO
449    
450    cph(
451    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
452    CADJ store dvsq, shsq = comlev1_kpp, key = ikppkey
453    #endif
454    cph)
455    
456  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
457  c     solve for viscosity, diffusivity, ghat, and hbl on "t-grid"  c     solve for viscosity, diffusivity, ghat, and hbl on "t-grid"
458  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
459    
460        DO j = jbot, jtop  c     precompute background vertical diffusivities, which are needed for
461           DO i = ibot, itop  c     matching diffusivities at bottom of KPP PBL
462          CALL CALC_3D_DIFFUSIVITY(
463         I        bi,bj,1-Olx,sNx+OLx,1-Oly,sNy+OLy,
464         I        GAD_SALINITY, .FALSE., .FALSE.,
465         O        KPPdiffKzS(1-Olx,1-Oly,1,bi,bj),
466         I        myThid)
467          CALL CALC_3D_DIFFUSIVITY(
468         I        bi,bj,1-Olx,sNx+OLx,1-Oly,sNy+OLy,
469         I        GAD_TEMPERATURE, .FALSE., .FALSE.,
470         O        KPPdiffKzT(1-Olx,1-Oly,1,bi,bj),
471         I        myThid)
472    
473          DO j = 1-OLy, sNy+OLy
474             DO i = 1-OLx, sNx+OLx
475              work1(i,j) = nzmax(i,j,bi,bj)              work1(i,j) = nzmax(i,j,bi,bj)
476              work2(i,j) = Fcori(i,j,bi,bj)              work2(i,j) = Fcori(i,j,bi,bj)
477           END DO           ENDDO
478        END DO        ENDDO
479        CALL TIMER_START('KPPMIX [KPP_CALC]', myThid)        CALL TIMER_START('KPPMIX [KPP_CALC]', myThid)
480        CALL KPPMIX (        CALL KPPMIX (
481       I       mytime, mythid       I       work1, shsq, dVsq, ustar
482       I     , work1, shsq, dVsq, ustar       I     , maskC(1-Olx,1-Oly,1,bi,bj)
483       I     , bo, bosol, dbloc, Ritop, work2       I     , bo, bosol, dbloc, Ritop, work2
484       I     , ikey       I     , KPPdiffKzS(1-Olx,1-Oly,1,bi,bj)
485         I     , KPPdiffKzT(1-Olx,1-Oly,1,bi,bj)
486         I     , ikppkey
487       O     , vddiff       O     , vddiff
488       U     , ghat       U     , ghat
489       O     , hbl )       O     , hbl
490         I     , mytime, mythid )
491        CALL TIMER_STOP ('KPPMIX [KPP_CALC]', myThid)        CALL TIMER_STOP ('KPPMIX [KPP_CALC]', myThid)
492    
 #ifdef ALLOW_AUTODIFF_TAMC  
 cph( storing not necessary  
 cphCADJ STORE vddiff, ghat  = comlev1_kpp, key = ikey  
 cph)  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
493  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
494  c     zero out land values and transfer to global variables  c     zero out land values and transfer to global variables
495  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 615  c--------------------------------------- Line 497  c---------------------------------------
497        DO j = jmin, jmax        DO j = jmin, jmax
498         DO i = imin, imax         DO i = imin, imax
499          DO k = 1, Nr          DO k = 1, Nr
500           KPPviscAz(i,j,k,bi,bj) = vddiff(i,j,k-1,1) * pMask(i,j,k,bi,bj)           km1 = max(1,k-1)
501           KPPdiffKzS(i,j,k,bi,bj)= vddiff(i,j,k-1,2) * pMask(i,j,k,bi,bj)           KPPviscAz(i,j,k,bi,bj) = vddiff(i,j,k-1,1) * maskC(i,j,k,bi,bj)
502           KPPdiffKzT(i,j,k,bi,bj)= vddiff(i,j,k-1,3) * pMask(i,j,k,bi,bj)       &        * maskC(i,j,km1,bi,bj)
503           KPPghat(i,j,k,bi,bj)   = ghat(i,j,k)       * pMask(i,j,k,bi,bj)           KPPdiffKzS(i,j,k,bi,bj)= vddiff(i,j,k-1,2) * maskC(i,j,k,bi,bj)
504          END DO       &        * maskC(i,j,km1,bi,bj)
505          KPPhbl(i,j,bi,bj) = hbl(i,j) * pMask(i,j,1,bi,bj)           KPPdiffKzT(i,j,k,bi,bj)= vddiff(i,j,k-1,3) * maskC(i,j,k,bi,bj)
506         END DO       &        * maskC(i,j,km1,bi,bj)
507        END DO           KPPghat(i,j,k,bi,bj)   = ghat(i,j,k)       * maskC(i,j,k,bi,bj)
508  #ifdef FRUGAL_KPP       &        * maskC(i,j,km1,bi,bj)
509        _EXCH_XYZ_R8(KPPviscAz  , myThid )          ENDDO
510        _EXCH_XYZ_R8(KPPdiffKzS , myThid )          k = 1
511        _EXCH_XYZ_R8(KPPdiffKzT , myThid )  #ifdef ALLOW_SHELFICE
512        _EXCH_XYZ_R8(KPPghat    , myThid )          if ( useShelfIce ) k = kTopC(i,j,bi,bj)
513        _EXCH_XY_R8 (KPPhbl     , myThid )  #endif /* ALLOW_SHELFICE */
514  #endif          KPPhbl(i,j,bi,bj) = hbl(i,j) * maskC(i,j,k,bi,bj)
515    
516           ENDDO
517          ENDDO
518    
519  #ifdef KPP_SMOOTH_VISC  #ifdef KPP_SMOOTH_VISC
520  c     horizontal smoothing of vertical viscosity  c     horizontal smoothing of vertical viscosity
521        DO k = 1, Nr        DO k = 1, Nr
522           CALL SMOOTH_HORIZ (           CALL SMOOTH_HORIZ (
523       I        k, bi, bj,       I        k, bi, bj,
524       U        KPPviscAz(1-OLx,1-OLy,k,bi,bj) )       U        KPPviscAz(1-OLx,1-OLy,k,bi,bj),
525        END DO       I        myThid )
526        _EXCH_XYZ_R8(KPPviscAz  , myThid )        ENDDO
527    C jmc: No EXCH inside bi,bj loop !!!
528    c     _EXCH_XYZ_R8(KPPviscAz  , myThid )
529  #endif /* KPP_SMOOTH_VISC */  #endif /* KPP_SMOOTH_VISC */
530    
531  #ifdef KPP_SMOOTH_DIFF  #ifdef KPP_SMOOTH_DIFF
# Line 646  c     horizontal smoothing of vertical d Line 533  c     horizontal smoothing of vertical d
533        DO k = 1, Nr        DO k = 1, Nr
534           CALL SMOOTH_HORIZ (           CALL SMOOTH_HORIZ (
535       I        k, bi, bj,       I        k, bi, bj,
536       U        KPPdiffKzS(1-OLx,1-OLy,k,bi,bj) )       U        KPPdiffKzS(1-OLx,1-OLy,k,bi,bj),
537         I        myThid )
538           CALL SMOOTH_HORIZ (           CALL SMOOTH_HORIZ (
539       I        k, bi, bj,       I        k, bi, bj,
540       U        KPPdiffKzT(1-OLx,1-OLy,k,bi,bj) )       U        KPPdiffKzT(1-OLx,1-OLy,k,bi,bj),
541        END DO       I        myThid )
542        _EXCH_XYZ_R8(KPPdiffKzS , myThid )        ENDDO
       _EXCH_XYZ_R8(KPPdiffKzT , myThid )  
543  #endif /* KPP_SMOOTH_DIFF */  #endif /* KPP_SMOOTH_DIFF */
544    
545    cph(
546    cph  crucial: this avoids full recomp./call of kppmix
547    CADJ store KPPhbl = comlev1_kpp, key = ikppkey
548    cph)
549    
550  C     Compute fraction of solar short-wave flux penetrating to  C     Compute fraction of solar short-wave flux penetrating to
551  C     the bottom of the mixing layer.  C     the bottom of the mixing layer.
552        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
# Line 664  C     the bottom of the mixing layer. Line 556  C     the bottom of the mixing layer.
556        ENDDO        ENDDO
557        CALL SWFRAC(        CALL SWFRAC(
558       I     (sNx+2*OLx)*(sNy+2*OLy), minusone,       I     (sNx+2*OLx)*(sNy+2*OLy), minusone,
559       I     mytime, mythid,       U     worka,
560       U     worka )       I     myTime, myIter, myThid )
561        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
562           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
563              KPPfrac(i,j,bi,bj) = worka(i,j)              KPPfrac(i,j,bi,bj) = worka(i,j)
# Line 674  C     the bottom of the mixing layer. Line 566  C     the bottom of the mixing layer.
566    
567        ENDIF        ENDIF
568    
569  #endif ALLOW_KPP  #endif /* ALLOW_KPP */
570    
571          RETURN
572          END
573    
574          SUBROUTINE KPP_CALC_DUMMY(
575         I     bi, bj, myTime, myIter, myThid )
576    C     *==========================================================*
577    C     | SUBROUTINE KPP_CALC_DUMMY                                |
578    C     | o Compute all KPP fields defined in KPP.h                |
579    C     | o Dummy routine for TAMC
580    C     *==========================================================*
581    C     | This subroutine serves as an interface between MITGCMUV  |
582    C     | code and NCOM 1-D routines in kpp_routines.F             |
583    C     *==========================================================*
584          IMPLICIT NONE
585    
586    #include "SIZE.h"
587    #include "EEPARAMS.h"
588    #include "PARAMS.h"
589    #include "KPP.h"
590    #include "KPP_PARAMS.h"
591    #include "GRID.h"
592    #include "GAD.h"
593    
594    c Routine arguments
595    c     bi, bj :: Current tile indices
596    c     myTime :: Current time in simulation
597    c     myIter :: Current iteration number in simulation
598    c     myThid :: My Thread Id. number
599    
600          INTEGER bi, bj
601          _RL     myTime
602          INTEGER myIter
603          INTEGER myThid
604    
605    #ifdef ALLOW_KPP
606    
607    c Local constants
608          integer i, j, k
609    
610          DO j=1-OLy,sNy+OLy
611             DO i=1-OLx,sNx+OLx
612                KPPhbl (i,j,bi,bj) = 1.0
613                KPPfrac(i,j,bi,bj) = 0.0
614                DO k = 1,Nr
615                   KPPghat   (i,j,k,bi,bj) = 0.0
616                   KPPviscAz (i,j,k,bi,bj) = viscAr
617                ENDDO
618             ENDDO
619          ENDDO
620    
621          CALL CALC_3D_DIFFUSIVITY(
622         I     bi,bj,1-Olx,sNx+OLx,1-Oly,sNy+OLy,
623         I     GAD_SALINITY, .FALSE., .FALSE.,
624         O     KPPdiffKzS(1-Olx,1-Oly,1,bi,bj),
625         I     myThid)
626          CALL CALC_3D_DIFFUSIVITY(
627         I     bi,bj,1-Olx,sNx+OLx,1-Oly,sNy+OLy,
628         I     GAD_TEMPERATURE, .FALSE., .FALSE.,
629         O     KPPdiffKzT(1-Olx,1-Oly,1,bi,bj),
630         I     myThid)
631    
632    #endif
633        RETURN        RETURN
634        END        END
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.42
  ViewVC Help
Powered by ViewVC 1.1.22