/[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.1 by adcroft, Wed Jun 21 19:45:47 2000 UTC revision 1.43 by heimbach, Mon Jun 4 18:47:35 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, fu, fv in the region (-2:sNx+4,-2:sNy+4).  c     values of uVel, vVel, surfaceForcingU, surfaceForcingV in the
94    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 97  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"
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
137    c     minusone, p0, p5, p25, p125, p0625
138    c     imin, imax, jmin, jmax  - array computation indices
139    
140          _RL        minusone
141          parameter( minusone=-1.0)
142          _RL        p0    , p5    , p25     , p125      , p0625
143          parameter( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 )
144          integer   imin      ,imax          ,jmin      ,jmax
145          parameter(imin=2-OLx,imax=sNx+OLx-1,jmin=2-OLy,jmax=sNy+OLy-1)
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
149  c     ustar  (nx,ny)       - surface friction velocity                  (m/s)  c     ustar  (nx,ny)       - surface friction velocity                  (m/s)
# Line 132  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 142  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        _RS worka (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        integer work1 ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy            )
172        _RS workb (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL worka ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
173  #ifdef FRUGAL_KPP        _RL work2 ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
174        integer work1(sNx,sNy)        _RL work3 ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
175        _RS work2    (sNx,sNy)        _RL ustar ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
176        _RS ustar    (sNx,sNy)        _RL bo    ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
177        _RS bo       (sNx,sNy)        _RL bosol ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
178        _RS bosol    (sNx,sNy)        _RL shsq  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
179        _RS shsq     (sNx,sNy,Nr)        _RL dVsq  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
180        _RS dVsq     (sNx,sNy,Nr)        _RL dbloc ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
181        _RS dbloc    (sNx,sNy,Nr)        _RL Ritop ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
182        _RS Ritop    (sNx,sNy,Nr)        _RL vddiff( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, 0:Nrp1, mdiff )
183        _RS vddiff   (sNx,sNy,0:Nrp1,mdiff)        _RL ghat  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
184        _RS ghat     (sNx,sNy,Nr)        _RL hbl   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
185        _RS hbl      (sNx,sNy)  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        _RS z0       (sNx,sNy)        _RL z0    ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
191        _RS zRef     (sNx,sNy)        _RL zRef  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
192        _RS uRef     (sNx,sNy)        _RL uRef  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
193        _RS vRef     (sNx,sNy)        _RL vRef  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
194  #endif /* KPP_ESTIMATE_UREF */  #endif /* KPP_ESTIMATE_UREF */
 #else /* FRUGAL_KPP */  
       integer work1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS work2    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS ustar    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS bo       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS bosol    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS shsq     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RS dVsq     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RS dbloc    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RS Ritop    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RS vddiff   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,0:Nrp1,mdiff)  
       _RS ghat     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RS hbl      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
 #ifdef KPP_ESTIMATE_UREF  
       _RS z0       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS zRef     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS uRef     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS vRef     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
 #endif /* KPP_ESTIMATE_UREF */  
 #endif /* FRUGAL_KPP */  
   
 c     imin,imax,jmin,jmax  - array indices  
       integer    imin   , imax      , jmin   , jmax  
       parameter( imin=-2, imax=sNx+3, jmin=-2, jmax=sNy+3 )  
195    
196  c     mixing process switches        _RL tempvar2
197        logical    lri        integer i, j, k, kp1, km1, im1, ip1, jm1, jp1
       parameter( lri = .true. )  
   
       _RS        p0    , p5    , p25     , p125      , p0625  
       parameter( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 )  
         
       _RS     tempVar  
       integer i, j, k, kp1, im1, ip1, jm1, jp1  
198    
199  #ifdef KPP_ESTIMATE_UREF  #ifdef KPP_ESTIMATE_UREF
200        _RS     dBdz1, dBdz2, ustarX, ustarY        _RL tempvar1, dBdz1, dBdz2, ustarX, ustarY
201  #endif  #endif
202    
203        IF (use_KPPmixing) THEN  #ifdef ALLOW_AUTODIFF_TAMC
204              act1 = bi - myBxLo(myThid)
205  CADJ STORE fu     (:,:  ,bi,bj)  = uvtape, key = ikey, byte = isbyte            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
206  CADJ STORE fv     (:,:  ,bi,bj)  = uvtape, key = ikey, byte = isbyte            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 247  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  #ifdef FRUGAL_KPP       I     ikppkey, bi, bj, myThid )
      O     , vddiff(1    ,1    ,1,1), vddiff(1    ,1    ,1,2)  
 #else  
      O     , vddiff(1-OLx,1-OLy,1,1), vddiff(1-OLx,1-OLy,1,2)  
 #endif  
      &     )  
260        CALL TIMER_STOP ('STATEKPP      [KPP_CALC]', myThid)        CALL TIMER_STOP ('STATEKPP      [KPP_CALC]', myThid)
261    
 #ifdef KPP_SMOOTH_DBLOC  
 c     horizontally smooth dbloc with a 121 filter  
 c     (stored in ghat to save space)  
   
       DO k = 1, Nr  
          CALL SMOOTH_HORIZ_RS (  
      I        k, bi, bj,  
      I        dbloc(1-OLx,1-OLy,k),  
      O        ghat (1-OLx,1-OLy,k) )  
       ENDDO  
   
 #else /* KPP_SMOOTH_DBLOC */  
   
262        DO k = 1, Nr        DO k = 1, Nr
 #ifdef FRUGAL_KPP  
          DO j = 1, sNy  
             DO i = 1, sNx  
 #else  
263           DO j = 1-OLy, sNy+OLy           DO j = 1-OLy, sNy+OLy
264              DO i = imin, imax              DO i = 1-OLx, sNx+OLx
 #endif  
265                 ghat(i,j,k) = dbloc(i,j,k)                 ghat(i,j,k) = dbloc(i,j,k)
266              ENDDO              ENDDO
267           ENDDO           ENDDO
268        ENDDO        ENDDO
269    
270    #ifdef KPP_SMOOTH_DBLOC
271    c     horizontally smooth dbloc with a 121 filter
272    c     smooth dbloc stored in ghat to save space
273    c     dbloc(k) is buoyancy gradientnote between k and k+1
274    c     levels therefore k+1 mask must be used
275    
276          DO k = 1, Nr-1
277             CALL SMOOTH_HORIZ (
278         I        k+1, bi, bj,
279         U        ghat (1-OLx,1-OLy,k),
280         I        myThid )
281          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 SMOOTH_HORIZ_RS (        CALL SMOOTH_HORIZ (
288       I     k, bi, bj,       I     1, bi, bj,
289       I     work2,       U     work2,
290       O     work2 )       I     myThid )
291        DO k = 1, Nr        DO k = 1, Nr
292           CALL SMOOTH_HORIZ_RS (           CALL SMOOTH_HORIZ (
293         I        k+1, bi, bj,
294         U        dbloc (1-OLx,1-OLy,k),
295         I        myThid )
296             CALL SMOOTH_HORIZ (
297       I        k, bi, bj,       I        k, bi, bj,
298       I        dbloc (1-OLx,1-OLy,k)  ,       U        Ritop (1-OLx,1-OLy,k),
299       O        dbloc (1-OLx,1-OLy,k)   )       I        myThid )
300           CALL SMOOTH_HORIZ_RS (           CALL SMOOTH_HORIZ (
301       I        k, bi, bj,       I        k, bi, bj,
302       I        Ritop (1-OLx,1-OLy,k)  ,       U        TTALPHA(1-OLx,1-OLy,k),
303       O        Ritop (1-OLx,1-OLy,k)   )       I        myThid )
304           CALL SMOOTH_HORIZ_RS (           CALL SMOOTH_HORIZ (
305       I        k, bi, bj,       I        k, bi, bj,
306       I        vddiff(1-OLx,1-OLy,k,1),       U        SSBETA(1-OLx,1-OLy,k),
307       O        vddiff(1-OLx,1-OLy,k,1) )       I        myThid )
          CALL SMOOTH_HORIZ_RS (  
      I        k, bi, bj,  
      I        vddiff(1-OLx,1-OLy,k,2),  
      O        vddiff(1-OLx,1-OLy,k,2) )  
308        ENDDO        ENDDO
309  #endif /* KPP_SMOOTH_DENS */  #endif /* KPP_SMOOTH_DENS */
310    
311        DO k = 1, Nr        DO k = 1, Nr
312  #ifdef FRUGAL_KPP           km1 = max(1,k-1)
          DO j = 1, sNy  
             DO i = 1, sNx  
 #else  
313           DO j = 1-OLy, sNy+OLy           DO j = 1-OLy, sNy+OLy
314              DO i = 1-OLx, sNx+OLx              DO i = 1-OLx, sNx+OLx
 #endif  
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 336  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 = fu * delZ(1)                                   (N/m^2)  c     taux / rho = surfaceForcingU                               (N/m^2)
371  c     tauy / rho = fv * 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*Qnet + beta*EmPmR) * delZ(1) / rho   (m^2/s^3)  c     bo    = - g * ( alpha*surfaceForcingT +
374  c     bosol = - g * alpha * Qsw * delZ(1) / rho                 (m^2/s^3)  c                     beta *surfaceForcingS ) / rho            (m^2/s^3)
375  c------------------------------------------------------------------------  c     bosol = - g * alpha * Qsw * drF(1) / rho                 (m^2/s^3)
   
 #ifdef FRUGAL_KPP  
       DO j = 1, sNy  
          jp1 = j + 1  
          DO i = 1, sNx  
 #else  
       DO j = jmin, jmax  
          jp1 = j + 1  
          DO i = imin, imax  
 #endif  
             ip1 = i+1  
             ustar(i,j) = SQRT( SQRT(  
      &           (fu(i,j,bi,bj) + fu(ip1,j,bi,bj)) * p5 * delZ(1) *  
      &           (fu(i,j,bi,bj) + fu(ip1,j,bi,bj)) * p5 * delZ(1) +  
      &           (fv(i,j,bi,bj) + fv(i,jp1,bi,bj)) * p5 * delZ(1) *  
      &           (fv(i,j,bi,bj) + fv(i,jp1,bi,bj)) * p5 * delZ(1) ))  
             bo(I,J) = - gravity *  
      &           ( vddiff(I,J,1,1) * Qnet(i,j,bi,bj) +  
      &             vddiff(I,J,1,2) * EmPmR(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  
   
 #ifndef FRUGAL_KPP  
 c     set array edges to zero  
       DO j = jmin, jmax  
          DO i = 1-OLx, imin-1  
             ustar(i,j) = p0  
             bo   (I,J) = p0  
             bosol(I,J) = p0  
          END DO  
          DO i = imax+1, sNx+OLx  
             ustar(i,j) = p0  
             bo   (I,J) = p0  
             bosol(I,J) = p0  
          END DO  
       END DO  
       DO i = 1-OLx, sNx+OLx  
          DO j = 1-OLy, jmin-1  
             ustar(i,j) = p0  
             bo   (I,J) = p0  
             bosol(I,J) = p0  
          END DO  
          DO j = jmax+1, sNy+OLy  
             ustar(i,j) = p0  
             bo   (I,J) = p0  
             bosol(I,J) = p0  
          END DO  
       END DO  
 #endif  
   
376  c------------------------------------------------------------------------  c------------------------------------------------------------------------
377  c     velocity shear  c     velocity shear
378  c     --------------  c     --------------
# Line 410  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     dVsq computation        CALL KPP_FORCING_SURF(
390         I     work2, surfaceForcingU, surfaceForcingV,
391  #ifdef KPP_ESTIMATE_UREF       I     surfaceForcingT, surfaceForcingS, surfaceForcingTice,
392         I     Qsw, ttalpha, ssbeta,
393  c     Get rid of vertical resolution dependence of dVsq term by       O     ustar, bo, bosol, dVsq,
394  c     estimating a surface velocity that is independent of first level       I     ikppkey, iMin, iMax, jMin, jMax, bi, bj, myTime, myThid )
395  c     thickness in the model.  First determine mixed layer depth hMix.  
396  c     Second zRef = espilon * hMix.  Third determine roughness length  CMLcph(
397  c     scale z0.  Third estimate reference velocity.  CMLCADJ store ustar = comlev1_kpp, key = ikppkey
398    CMLcph)
 #ifdef FRUGAL_KPP  
       DO j = 1, sNy  
          jp1 = j + 1  
          DO i = 1, sNx  
 #else  
       DO j = jmin, jmax  
          jp1 = j + 1  
          DO i = imin, imax  
 #endif /* FRUGAL_KPP */  
             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  
                tempVar = SQRT ( 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)) ) )  
                z0(i,j) = rF(2) *  
      &                   ( rF(3) * LOG ( rF(3) / rF(2) ) /  
      &                     ( rF(3) - rF(2) ) -  
      &                     tempVar * 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 = ( fu(i,j,bi,bj) + fu(ip1,j,bi,bj) ) * p5  
                   ustarY = ( fv(i,j,bi,bj) + fu(i,jp1,bi,bj) ) * p5  
                   tempVar = MAX ( phepsi,  
      &                 SQRT ( ustarX * ustarX + ustarY * ustarY ) )  
                   tempVar = ustar(i,j) *  
      &                 ( LOG ( zRef(i,j) / rF(2) ) +  
      &                 z0(i,j) / zRef(i,j) - z0(i,j) / rF(2) ) /  
      &                 vonK / tempVar  
                   uRef(i,j) = uRef(i,j) + ustarX * tempVar  
                   vRef(i,j) = vRef(i,j) + ustarY * tempVar  
                ENDIF  
   
          END DO  
       END DO  
   
       IF (KPPmixingMaps) THEN  
 #ifdef FRUGAL_KPP  
          CALL PRINT_MAPRS(  
      I        zRef, 'zRef', PRINT_MAP_XY,  
      I        1, sNx, 1, sNy, 1, 1, 1, 1,  
      I        1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 )  
          CALL PRINT_MAPRS(  
      I        z0, 'z0', PRINT_MAP_XY,  
      I        1, sNx, 1, sNy, 1, 1, 1, 1,  
      I        1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 )  
          CALL PRINT_MAPRS(  
      I        uRef, 'uRef', PRINT_MAP_XY,  
      I        1, sNx, 1, sNy, 1, 1, 1, 1,  
      I        1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 )  
          CALL PRINT_MAPRS(  
      I        vRef, 'vRef', PRINT_MAP_XY,  
      I        1, sNx, 1, sNy, 1, 1, 1, 1,  
      I        1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 )  
 #else  
          CALL PRINT_MAPRS(  
      I        zRef, 'zRef', PRINT_MAP_XY,  
      I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1, 1, 1,  
      I        1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 )  
          CALL PRINT_MAPRS(  
      I        z0, 'z0', PRINT_MAP_XY,  
      I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1, 1, 1,  
      I        1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 )  
          CALL PRINT_MAPRS(  
      I        uRef, 'uRef', PRINT_MAP_XY,  
      I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1, 1, 1,  
      I        1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 )  
          CALL PRINT_MAPRS(  
      I        vRef, 'vRef', PRINT_MAP_XY,  
      I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1, 1, 1,  
      I        1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 )  
 #endif  
       ENDIF  
   
       DO k = 1, Nr  
 #ifdef FRUGAL_KPP  
          DO j = 1, sNy  
             jm1 = j - 1  
             jp1 = j + 1  
             DO i = 1, sNx  
 #else  
          DO j = jmin, jmax  
             jm1 = j - 1  
             jp1 = j + 1  
             DO i = imin, imax  
 #endif /* FRUGAL_KPP */  
                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  #ifdef FRUGAL_KPP         DO j = 1-OLy, sNy+OLy
403           DO j = 1, sNy          DO i = 1-OLx, sNx+OLx
404              jm1 = j - 1           shsq(i,j,k) = p0
405              jp1 = j + 1          ENDDO
406              DO i = 1, sNx         ENDDO
407  #else        ENDDO
          DO j = jmin, jmax  
             jm1 = j - 1  
             jp1 = j + 1  
             DO i = imin, imax  
 #endif /* FRUGAL_KPP */  
                im1 = i - 1  
                ip1 = i + 1  
                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  #ifdef FRUGAL_KPP         DO j = jmin, jmax
413           DO j = 1, sNy          jm1 = j - 1
414              jm1 = j - 1          jp1 = j + 1
415              jp1 = j + 1          DO i = imin, imax
416              DO i = 1, sNx           im1 = i - 1
417  #else           ip1 = i + 1
418           DO j = jmin, jmax           shsq(i,j,k) = p5 * (
419              jm1 = j - 1       &        (uVel(i,  j,  k,bi,bj)-uVel(i,  j,  kp1,bi,bj)) *
420              jp1 = j + 1       &        (uVel(i,  j,  k,bi,bj)-uVel(i,  j,  kp1,bi,bj)) +
421              DO i = imin, imax       &        (uVel(ip1,j,  k,bi,bj)-uVel(ip1,j,  kp1,bi,bj)) *
422  #endif /* FRUGAL_KPP */       &        (uVel(ip1,j,  k,bi,bj)-uVel(ip1,j,  kp1,bi,bj)) +
423                 im1 = i - 1       &        (vVel(i,  j,  k,bi,bj)-vVel(i,  j,  kp1,bi,bj)) *
424                 ip1 = i + 1       &        (vVel(i,  j,  k,bi,bj)-vVel(i,  j,  kp1,bi,bj)) +
425                 shsq(i,j,k) = p5 * (       &        (vVel(i,  jp1,k,bi,bj)-vVel(i,  jp1,kp1,bi,bj)) *
426       $              (uVel(i,  j,  k,bi,bj)-uVel(i,  j,  kp1,bi,bj)) *       &        (vVel(i,  jp1,k,bi,bj)-vVel(i,  jp1,kp1,bi,bj)) )
      $              (uVel(i,  j,  k,bi,bj)-uVel(i,  j,  kp1,bi,bj)) +  
      $              (uVel(ip1,j,  k,bi,bj)-uVel(ip1,j,  kp1,bi,bj)) *  
      $              (uVel(ip1,j,  k,bi,bj)-uVel(ip1,j,  kp1,bi,bj)) +  
      $              (vVel(i,  j,  k,bi,bj)-vVel(i,  j,  kp1,bi,bj)) *  
      $              (vVel(i,  j,  k,bi,bj)-vVel(i,  j,  kp1,bi,bj)) +  
      $              (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)) )
 #endif  
             END DO  
          END DO  
       END DO  
   
 c     shsq @ Nr computation  
 #ifdef FRUGAL_KPP  
       DO j = 1, sNy  
          DO i = 1, sNx  
 #else  
       DO j = jmin, jmax  
          DO i = imin, imax  
445  #endif  #endif
446              shsq(i,j,Nr) = p0          ENDDO
447           END DO         ENDDO
448        END DO        ENDDO
449    
450  #ifndef FRUGAL_KPP  cph(
451  c     set array edges to zero  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
452        DO k = 1, Nr  CADJ store dvsq, shsq = comlev1_kpp, key = ikppkey
          DO j = jmin, jmax  
             DO i = 1-OLx, imin-1  
                shsq(i,j,k) = p0  
                dVsq(i,j,k) = p0  
             END DO  
             DO i = imax+1, sNx+OLx  
                shsq(i,j,k) = p0  
                dVsq(i,j,k) = p0  
             END DO  
          END DO  
          DO i = 1-OLx, sNx+OLx  
             DO j = 1-OLy, jmin-1  
                shsq(i,j,k) = p0  
                dVsq(i,j,k) = p0  
             END DO  
             DO j = jmax+1, sNy+OLy  
                shsq(i,j,k) = p0  
                dVsq(i,j,k) = p0  
             END DO  
          END DO  
       END DO  
453  #endif  #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  #ifdef FRUGAL_KPP  c     precompute background vertical diffusivities, which are needed for
461        DO j = 1, sNy  c     matching diffusivities at bottom of KPP PBL
462           DO i = 1, sNx        CALL CALC_3D_DIFFUSIVITY(
463  #else       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        DO j = 1-OLy, sNy+OLy
474           DO i = 1-OLx, sNx+OLx           DO i = 1-OLx, sNx+OLx
 #endif  
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       lri, work1, shsq, dVsq, ustar       I       work1, shsq, dVsq, ustar
482         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, myIter, mythid )
   
491        CALL TIMER_STOP ('KPPMIX [KPP_CALC]', myThid)        CALL TIMER_STOP ('KPPMIX [KPP_CALC]', myThid)
492    
       IF (KPPmixingMaps) THEN  
 #ifdef FRUGAL_KPP  
        CALL PRINT_MAPRS(  
      I      hbl, 'hbl', PRINT_MAP_XY,  
      I      1, sNx, 1, sNy, 1, 1, 1, 1,  
      I      1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 )  
 #else  
        CALL PRINT_MAPRS(  
      I      hbl, 'hbl', PRINT_MAP_XY,  
      I      1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1, 1, 1,  
      I      1, sNx, 1, sNy, 1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1 )  
 #endif  
        ENDIF  
   
 CADJ STORE vddiff, ghat  = uvtape, key = ikey  
   
493  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
494  c     zero out land values,  c     zero out land values and transfer to global variables
 c     make sure coefficients are within reasonable bounds,  
 c     and transfer to global variables  
495  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
496    
 #ifdef FRUGAL_KPP  
       DO j = 1, sNy  
          DO i = 1, sNx  
 #else  
497        DO j = jmin, jmax        DO j = jmin, jmax
498           DO i = imin, imax         DO i = imin, imax
499  #endif          DO k = 1, Nr
500              DO k = 1, Nr           km1 = max(1,k-1)
501  c KPPviscAz           KPPviscAz(i,j,k,bi,bj) = vddiff(i,j,k-1,1) * maskC(i,j,k,bi,bj)
502                 tempVar = min( maxKPPviscAz(k), vddiff(i,j,k-1,1) )       &        * maskC(i,j,km1,bi,bj)
503                 tempVar = max( minKPPviscAz, tempVar )           KPPdiffKzS(i,j,k,bi,bj)= vddiff(i,j,k-1,2) * maskC(i,j,k,bi,bj)
504                 KPPviscAz(i,j,k,bi,bj) = tempVar*pMask(i,j,k,bi,bj)       &        * maskC(i,j,km1,bi,bj)
505  c KPPdiffKzS           KPPdiffKzT(i,j,k,bi,bj)= vddiff(i,j,k-1,3) * maskC(i,j,k,bi,bj)
506                 tempVar = min( maxKPPdiffKzS, vddiff(i,j,k-1,2) )       &        * maskC(i,j,km1,bi,bj)
507                 tempVar = max( minKPPdiffKzS, tempVar )           KPPghat(i,j,k,bi,bj)   = ghat(i,j,k)       * maskC(i,j,k,bi,bj)
508                 KPPdiffKzS(i,j,k,bi,bj) = tempVar*pMask(i,j,k,bi,bj)       &        * maskC(i,j,km1,bi,bj)
509  c KPPdiffKzT          ENDDO
510                 tempVar = min( maxKPPdiffKzT, vddiff(i,j,k-1,3) )          k = 1
511                 tempVar = max( minKPPdiffKzT, tempVar )  #ifdef ALLOW_SHELFICE
512                 KPPdiffKzT(i,j,k,bi,bj) = tempVar*pMask(i,j,k,bi,bj)          if ( useShelfIce ) k = kTopC(i,j,bi,bj)
513  c KPPghat  #endif /* ALLOW_SHELFICE */
514                 tempVar = min( maxKPPghat, ghat(i,j,k) )          KPPhbl(i,j,bi,bj) = hbl(i,j) * maskC(i,j,k,bi,bj)
515                 tempVar = max( minKPPghat, tempVar )  
516                 KPPghat(i,j,k,bi,bj) = tempVar*pMask(i,j,k,bi,bj)         ENDDO
517              END DO        ENDDO
 c KPPhbl: set to -zgrid(1) over land  
             KPPhbl(i,j,bi,bj) = (hbl(i,j) + zgrid(1))  
      &           * pMask(i,j,1,bi,bj) -  
      &           zgrid(1)  
          END DO  
       END DO  
 #ifdef FRUGAL_KPP  
       _EXCH_XYZ_R8(KPPviscAz  , myThid )  
       _EXCH_XYZ_R8(KPPdiffKzS , myThid )  
       _EXCH_XYZ_R8(KPPdiffKzT , myThid )  
       _EXCH_XYZ_R8(KPPghat    , myThid )  
       _EXCH_XY_R8 (KPPhbl     , myThid )  
 #endif  
518    
519  #ifdef KPP_SMOOTH_VISC  #ifdef KPP_SMOOTH_VISC
520  c     horizontal smoothing of vertical viscosity  c     horizontal smoothing of vertical viscosity
 c     as coded requires FRUGAL_KPP and OLx=4, OLy=4  
 c     alternatively could recode with OLx=5, OLy=5  
   
521        DO k = 1, Nr        DO k = 1, Nr
522           CALL SMOOTH_HORIZ_RL (           CALL SMOOTH_HORIZ (
523       I        k, bi, bj,       I        k, bi, bj,
524       I        KPPviscAz(1-OLx,1-OLy,k,bi,bj),       U        KPPviscAz(1-OLx,1-OLy,k,bi,bj),
525       O        KPPviscAz(1-OLx,1-OLy,k,bi,bj) )       I        myThid )
526        END DO        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
532  c     horizontal smoothing of vertical diffusivity  c     horizontal smoothing of vertical diffusivity
 c     as coded requires FRUGAL_KPP and OLx=4, OLy=4  
 c     alternatively could recode with OLx=5, OLy=5  
   
533        DO k = 1, Nr        DO k = 1, Nr
534           CALL SMOOTH_HORIZ_RL (           CALL SMOOTH_HORIZ (
535       I        k, bi, bj,       I        k, bi, bj,
536       I        KPPdiffKzS(1-OLx,1-OLy,k,bi,bj),       U        KPPdiffKzS(1-OLx,1-OLy,k,bi,bj),
537       O        KPPdiffKzS(1-OLx,1-OLy,k,bi,bj) )       I        myThid )
538           CALL SMOOTH_HORIZ_RL (           CALL SMOOTH_HORIZ (
539       I        k, bi, bj,       I        k, bi, bj,
540       I        KPPdiffKzT(1-OLx,1-OLy,k,bi,bj),       U        KPPdiffKzT(1-OLx,1-OLy,k,bi,bj),
541       O        KPPdiffKzT(1-OLx,1-OLy,k,bi,bj) )       I        myThid )
542        END DO        ENDDO
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.
# Line 847  C     the bottom of the mixing layer. Line 555  C     the bottom of the mixing layer.
555           ENDDO           ENDDO
556        ENDDO        ENDDO
557        CALL SWFRAC(        CALL SWFRAC(
558       I     (sNx+2*OLx)*(sNy+2*OLy), -1., worka,       I     (sNx+2*OLx)*(sNy+2*OLy), minusone,
559       O     workb )       U     worka,
560         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) = workb(i,j)              KPPfrac(i,j,bi,bj) = worka(i,j)
564           ENDDO           ENDDO
565        ENDDO        ENDDO
           
 CADJ STORE KPPhbl    (:,:  ,bi,bj)  = uvtape, key = ikey, byte = isbyte  
 CADJ STORE KPPghat   (:,:,:,bi,bj)  = uvtape, key = ikey, byte = isbyte  
 CADJ STORE KPPviscAz (:,:,:,bi,bj)  = uvtape, key = ikey, byte = isbyte  
 CADJ STORE KPPdiffKzT(:,:,:,bi,bj)  = uvtape, key = ikey, byte = isbyte  
 CADJ STORE KPPdiffKzS(:,:,:,bi,bj)  = uvtape, key = ikey, byte = isbyte  
566    
567        ENDIF        ENDIF
       ENDIF  
568    
569  #endif ALLOW_KPP  #endif /* ALLOW_KPP */
570    
571        RETURN        RETURN
572        END        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
634          END
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.43
  ViewVC Help
Powered by ViewVC 1.1.22