/[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.39 by jmc, Mon Apr 30 13:49:40 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    CBOP
7    C !ROUTINE: KPP_CALC
8    
9    C !INTERFACE: ==========================================================
10        subroutine KPP_CALC(        subroutine KPP_CALC(
11       I     bi, bj, myTime, myThid )       I     bi, bj, myTime, myThid )
12    
13    C !DESCRIPTION: \bv
14  C     /==========================================================\  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                |
# 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 - array indices on which to apply calculations
124  c     myTime - Current time in simulation  c     myTime - Current time in simulation
# Line 121  c     myTime - Current time in simulatio Line 129  c     myTime - Current time in simulatio
129    
130  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
131    
132    C !LOCAL VARIABLES: ====================================================
133  c Local constants  c Local constants
134  c     minusone, p0, p5, p25, p125, p0625  c     minusone, p0, p5, p25, p125, p0625
135  c     imin, imax, jmin, jmax  - array computation indices  c     imin, imax, jmin, jmax  - array computation indices
136    
137        _RL        minusone        _RL        minusone
138        parameter( minusone=-1.0)        parameter( minusone=-1.0)
139        _KPP_RL    p0    , p5    , p25     , p125      , p0625        _RL        p0    , p5    , p25     , p125      , p0625
140        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 )
141        integer    imin      , imax        , jmin      , jmax        integer   imin      ,imax          ,jmin      ,jmax
142  #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  
143    
144  c Local arrays and variables  c Local arrays and variables
145  c     work?  (nx,ny)       - horizontal working arrays  c     work?  (nx,ny)       - horizontal working arrays
# Line 150  c                            for ri_iwmi Line 155  c                            for ri_iwmi
155  c     Ritop  (nx,ny,Nr)    - numerator of bulk richardson number  c     Ritop  (nx,ny,Nr)    - numerator of bulk richardson number
156  c                            at grid levels for bldepth  c                            at grid levels for bldepth
157  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)
158  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)
159  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)
160  c     ghat   (nx,ny,Nr)    - nonlocal transport coefficient           (s/m^2)  c     ghat   (nx,ny,Nr)    - nonlocal transport coefficient           (s/m^2)
161  c     hbl    (nx,ny)       - mixing layer depth                           (m)  c     hbl    (nx,ny)       - mixing layer depth                           (m)
162  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 165  c     zRef   (nx,ny)       - Reference d
165  c     uRef   (nx,ny)       - Reference zonal velocity                   (m/s)  c     uRef   (nx,ny)       - Reference zonal velocity                   (m/s)
166  c     vRef   (nx,ny)       - Reference meridional velocity              (m/s)  c     vRef   (nx,ny)       - Reference meridional velocity              (m/s)
167    
168        _RL     worka ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )        integer work1 ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy            )
169        integer work1 ( ibot:itop    , jbot:jtop                    )        _RL worka ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
170        _KPP_RL work2 ( ibot:itop    , jbot:jtop                    )        _RL work2 ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
171        _KPP_RL ustar ( ibot:itop    , jbot:jtop                    )        _RL work3 ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
172        _KPP_RL bo    ( ibot:itop    , jbot:jtop                    )        _RL ustar ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
173        _KPP_RL bosol ( ibot:itop    , jbot:jtop                    )        _RL bo    ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
174        _KPP_RL shsq  ( ibot:itop    , jbot:jtop    , Nr            )        _RL bosol ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
175        _KPP_RL dVsq  ( ibot:itop    , jbot:jtop    , Nr            )        _RL shsq  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
176        _KPP_RL dbloc ( ibot:itop    , jbot:jtop    , Nr            )        _RL dVsq  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
177        _KPP_RL Ritop ( ibot:itop    , jbot:jtop    , Nr            )        _RL dbloc ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
178        _KPP_RL vddiff( ibot:itop    , jbot:jtop    , 0:Nrp1, mdiff )        _RL Ritop ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
179        _KPP_RL ghat  ( ibot:itop    , jbot:jtop    , Nr            )        _RL vddiff( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, 0:Nrp1, mdiff )
180        _KPP_RL hbl   ( ibot:itop    , jbot:jtop                    )        _RL ghat  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr            )
181          _RL hbl   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
182    cph(
183          _RL TTALPHA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
184          _RL SSBETA ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
185    cph)
186  #ifdef KPP_ESTIMATE_UREF  #ifdef KPP_ESTIMATE_UREF
187        _KPP_RL z0    ( ibot:itop    , jbot:jtop                    )        _RL z0    ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
188        _KPP_RL zRef  ( ibot:itop    , jbot:jtop                    )        _RL zRef  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
189        _KPP_RL uRef  ( ibot:itop    , jbot:jtop                    )        _RL uRef  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
190        _KPP_RL vRef  ( ibot:itop    , jbot:jtop                    )        _RL vRef  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy                )
191  #endif /* KPP_ESTIMATE_UREF */  #endif /* KPP_ESTIMATE_UREF */
192                
193        _KPP_RL tempvar1, tempvar2        _RL tempvar2
194        integer i, j, k, kp1, im1, ip1, jm1, jp1        integer i, j, k, kp1, km1, im1, ip1, jm1, jp1
195    
196  #ifdef KPP_ESTIMATE_UREF  #ifdef KPP_ESTIMATE_UREF
197        _KPP_RL dBdz1, dBdz2, ustarX, ustarY        _RL tempvar1, dBdz1, dBdz2, ustarX, ustarY
198  #endif  #endif
199    
200    #ifdef ALLOW_AUTODIFF_TAMC
201              act1 = bi - myBxLo(myThid)
202              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
203              act2 = bj - myByLo(myThid)
204              max2 = myByHi(myThid) - myByLo(myThid) + 1
205              act3 = myThid - 1
206              max3 = nTx*nTy
207              act4 = ikey_dynamics - 1
208              ikppkey = (act1 + 1) + act2*max1
209         &                      + act3*max1*max2
210         &                      + act4*max1*max2*max3
211    #endif /* ALLOW_AUTODIFF_TAMC */
212    CEOP
213    
214  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?
215        IF ( DIFFERENT_MULTIPLE(kpp_freq,myTime,myTime-deltaTClock) .OR.        IF ( DIFFERENT_MULTIPLE(kpp_freq,myTime,deltaTClock)
216       1     myTime .EQ. startTime ) THEN       1     .OR. myTime .EQ. startTime ) THEN
217                    
218  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
219  c     prepare input arrays for subroutine "kppmix" to compute  c     prepare input arrays for subroutine "kppmix" to compute
# Line 227  c--------------------------------------- Line 251  c---------------------------------------
251    
252        CALL TIMER_START('STATEKPP      [KPP_CALC]', myThid)        CALL TIMER_START('STATEKPP      [KPP_CALC]', myThid)
253        CALL STATEKPP(        CALL STATEKPP(
254       I       bi, bj, myThid       I       ikppkey, bi, bj, myThid
255       O     , work2, dbloc, Ritop       O     , work2, dbloc, Ritop
256       O     , vddiff(ibot,jbot,1,1), vddiff(ibot,jbot,1,2)       O     , TTALPHA, SSBETA
257       &     )       &     )
258        CALL TIMER_STOP ('STATEKPP      [KPP_CALC]', myThid)        CALL TIMER_STOP ('STATEKPP      [KPP_CALC]', myThid)
259    
260        DO k = 1, Nr        DO k = 1, Nr
261           DO j = jbot, jtop           DO j = 1-OLy, sNy+OLy
262              DO i = ibot, itop              DO i = 1-OLx, sNx+OLx
263                 ghat(i,j,k) = dbloc(i,j,k)                 ghat(i,j,k) = dbloc(i,j,k)
264              ENDDO              ENDDO
265           ENDDO           ENDDO
# Line 248  c     dbloc(k) is buoyancy gradientnote Line 272  c     dbloc(k) is buoyancy gradientnote
272  c     levels therefore k+1 mask must be used  c     levels therefore k+1 mask must be used
273    
274        DO k = 1, Nr-1        DO k = 1, Nr-1
275           CALL KPP_SMOOTH_HORIZ (           CALL SMOOTH_HORIZ (
276       I        k+1, bi, bj,       I        k+1, bi, bj,
277       U        ghat (ibot,jbot,k) )       U        ghat (1-OLx,1-OLy,k) )
278        ENDDO        ENDDO
279    
280  #endif /* KPP_SMOOTH_DBLOC */  #endif /* KPP_SMOOTH_DBLOC */
281    
282  #ifdef KPP_SMOOTH_DENS  #ifdef KPP_SMOOTH_DENS
283  c     horizontally smooth density related quantities with 121 filters  c     horizontally smooth density related quantities with 121 filters
284        CALL KPP_SMOOTH_HORIZ (        CALL SMOOTH_HORIZ (
285       I     1, bi, bj,       I     1, bi, bj,
286       U     work2 )       U     work2 )
287        DO k = 1, Nr        DO k = 1, Nr
288           CALL KPP_SMOOTH_HORIZ (           CALL SMOOTH_HORIZ (
289       I        k+1, bi, bj,       I        k+1, bi, bj,
290       U        dbloc (ibot,jbot,k) )       U        dbloc (1-OLx,1-OLy,k) )
291           CALL KPP_SMOOTH_HORIZ (           CALL SMOOTH_HORIZ (
292       I        k, bi, bj,       I        k, bi, bj,
293       U        Ritop (ibot,jbot,k)  )       U        Ritop (1-OLx,1-OLy,k)  )
294           CALL KPP_SMOOTH_HORIZ (           CALL SMOOTH_HORIZ (
295       I        k, bi, bj,       I        k, bi, bj,
296       U        vddiff(ibot,jbot,k,1) )       U        TTALPHA(1-OLx,1-OLy,k) )
297           CALL KPP_SMOOTH_HORIZ (           CALL SMOOTH_HORIZ (
298       I        k, bi, bj,       I        k, bi, bj,
299       U        vddiff(ibot,jbot,k,2) )       U        SSBETA(1-OLx,1-OLy,k) )
300        ENDDO        ENDDO
301  #endif /* KPP_SMOOTH_DENS */  #endif /* KPP_SMOOTH_DENS */
302    
303        DO k = 1, Nr        DO k = 1, Nr
304           DO j = jbot, jtop           km1 = max(1,k-1)
305              DO i = ibot, itop           DO j = 1-OLy, sNy+OLy
306                DO i = 1-OLx, sNx+OLx
307    
308  c     zero out dbloc over land points (so that the convective  c     zero out dbloc over land points (so that the convective
309  c     part of the interior mixing can be diagnosed)  c     part of the interior mixing can be diagnosed)
310                 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)
311                 ghat(i,j,k)  = ghat(i,j,k)  * pMask(i,j,k,bi,bj)       &              * maskC(i,j,km1,bi,bj)
312                 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)
313         &              * maskC(i,j,km1,bi,bj)
314                   Ritop(i,j,k) = Ritop(i,j,k) * maskC(i,j,k,bi,bj)
315         &              * maskC(i,j,km1,bi,bj)
316                 if(k.eq.nzmax(i,j,bi,bj)) then                 if(k.eq.nzmax(i,j,bi,bj)) then
317                    dbloc(i,j,k) = p0                    dbloc(i,j,k) = p0
318                    ghat(i,j,k)  = p0                    ghat(i,j,k)  = p0
# Line 300  c     so that the subroutine "bldepth" w Line 328  c     so that the subroutine "bldepth" w
328           END DO           END DO
329        END DO        END DO
330    
331    cph(
332    cph  this avoids a single or double recomp./call of statekpp
333    CADJ store work2              = comlev1_kpp, key = ikppkey
334    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
335    CADJ store dbloc, Ritop, ghat = comlev1_kpp, key = ikppkey
336    CADJ store vddiff             = comlev1_kpp, key = ikppkey
337    CADJ store TTALPHA, SSBETA    = comlev1_kpp, key = ikppkey
338    #endif
339    cph)
340    
341  c------------------------------------------------------------------------  c------------------------------------------------------------------------
342  c     friction velocity, turbulent and radiative surface buoyancy forcing  c     friction velocity, turbulent and radiative surface buoyancy forcing
343  c     -------------------------------------------------------------------  c     -------------------------------------------------------------------
344  c     taux / rho = SurfaceTendencyU * delZ(1)                     (N/m^2)  c     taux / rho = surfaceForcingU                               (N/m^2)
345  c     tauy / rho = SurfaceTendencyV * delZ(1)                     (N/m^2)  c     tauy / rho = surfaceForcingV                               (N/m^2)
346  c     ustar = sqrt( sqrt( taux^2 + tauy^2 ) / rho )                 (m/s)  c     ustar = sqrt( sqrt( taux^2 + tauy^2 ) / rho )                (m/s)
347  c     bo    = - g * ( alpha*SurfaceTendencyT +  c     bo    = - g * ( alpha*surfaceForcingT +
348  c                     beta *SurfaceTendencyS ) * delZ(1) / rho  (m^2/s^3)  c                     beta *surfaceForcingS ) / rho            (m^2/s^3)
349  c     bosol = - g * alpha * Qsw * delZ(1) / rho                 (m^2/s^3)  c     bosol = - g * alpha * Qsw * drF(1) / rho                 (m^2/s^3)
350  c------------------------------------------------------------------------  c------------------------------------------------------------------------
351    
352  c initialize arrays to zero  c initialize arrays to zero
353        DO j = jbot, jtop        DO j = 1-OLy, sNy+OLy
354           DO i = ibot, itop           DO i = 1-OLx, sNx+OLx
355              ustar(i,j) = p0              ustar(i,j) = p0
356              bo   (I,J) = p0              bo   (I,J) = p0
357              bosol(I,J) = p0              bosol(I,J) = p0
# Line 324  c initialize arrays to zero Line 362  c initialize arrays to zero
362         jp1 = j + 1         jp1 = j + 1
363         DO i = imin, imax         DO i = imin, imax
364          ip1 = i+1          ip1 = i+1
365          tempVar1 =          work3(i,j) =
366       &   (SurfaceTendencyU(i,j,bi,bj) + SurfaceTendencyU(ip1,j,bi,bj)) *       &   (surfaceForcingU(i,j,bi,bj) + surfaceForcingU(ip1,j,bi,bj)) *
367       &   (SurfaceTendencyU(i,j,bi,bj) + SurfaceTendencyU(ip1,j,bi,bj)) +       &   (surfaceForcingU(i,j,bi,bj) + surfaceForcingU(ip1,j,bi,bj)) +
368       &   (SurfaceTendencyV(i,j,bi,bj) + SurfaceTendencyV(i,jp1,bi,bj)) *       &   (surfaceForcingV(i,j,bi,bj) + surfaceForcingV(i,jp1,bi,bj)) *
369       &   (SurfaceTendencyV(i,j,bi,bj) + SurfaceTendencyV(i,jp1,bi,bj))       &   (surfaceForcingV(i,j,bi,bj) + surfaceForcingV(i,jp1,bi,bj))
370          if ( tempVar1 .lt. (phepsi*phepsi) ) then         END DO
371             ustar(i,j) = SQRT( phepsi * p5 * delZ(1) )        END DO
372    cph(
373    CADJ store work3 = comlev1_kpp, key = ikppkey
374    cph)
375          DO j = jmin, jmax
376           jp1 = j + 1
377           DO i = imin, imax
378            ip1 = i+1
379    
380            if ( work3(i,j) .lt. (phepsi*phepsi*drF(1)*drF(1)) ) then
381               ustar(i,j) = SQRT( phepsi * p5 * drF(1) )
382          else          else
383             tempVar2 =  SQRT( tempVar1 ) * p5 * delZ(1)             tempVar2 =  SQRT( work3(i,j) ) * p5
384             ustar(i,j) = SQRT( tempVar2 )             ustar(i,j) = SQRT( tempVar2 )
385          endif          endif
386    
387           END DO
388          END DO
389    
390    CML#ifdef ALLOW_SHELFICE
391    CMLC     For the pbl parameterisation to work underneath the ice shelves
392    CMLC     it needs to know the surface (ice-ocean) fluxes. However, masking
393    CMLC     and indexing problems make this part of the code not work
394    CMLC     underneath the ice shelves and the following lines are only here
395    CMLC     to remind me that this still needs to be sorted out.
396    CML      IF ( useShelfIce ) THEN
397    CML      DO j = jmin, jmax
398    CML       jp1 = j + 1
399    CML       DO i = imin, imax
400    CML        ip1 = i+1
401    CML        k   = kTopC(I,J,bi,bj)
402    CML        bo(I,J) = - gravity *
403    CML     &       ( TTALPHA(I,J,K) * ( maskC(I,J,1,bi,bj)
404    CML     &       * ( surfaceForcingT(i,j,bi,bj)
405    CML     &       + surfaceForcingTice(i,j,bi,bj) )
406    CML     &       + shelficeForcingT(i,j,bi,bj) )
407    CML     &       + SSBETA(I,J,K) * ( maskC(I,J,1,bi,bj)
408    CML     &       * surfaceForcingS(i,j,bi,bj)
409    CML     &       + shelficeForcingS(i,j,bi,bj) ) )
410    CML     &       / work2(I,J)
411    CML
412    CML        bosol(I,J) = gravity * TTALPHA(I,J,1) * Qsw(i,j,bi,bj)
413    CML     &       * maskC(I,J,1,bi,bj) * recip_Cp*recip_rhoConst
414    CML     &       / work2(I,J)
415    CML
416    CML       END DO
417    CML      END DO
418    CML      ELSE
419    CML#endif /* ALLOW_SHELFICE */
420          DO j = jmin, jmax
421           jp1 = j + 1
422           DO i = imin, imax
423            ip1 = i+1
424    
425          bo(I,J) = - gravity *          bo(I,J) = - gravity *
426       &       ( vddiff(I,J,1,1) * SurfaceTendencyT(i,j,bi,bj) +       &       ( TTALPHA(I,J,1) * (surfaceForcingT(i,j,bi,bj)+
427       &         vddiff(I,J,1,2) * SurfaceTendencyS(i,j,bi,bj)       &       surfaceForcingTice(i,j,bi,bj)) +
428       &       ) *       &       SSBETA(I,J,1) * surfaceForcingS(i,j,bi,bj) )
429       &       delZ(1) / work2(I,J)       &       / work2(I,J)
430          bosol(I,J) = - gravity * vddiff(I,J,1,1) * Qsw(i,j,bi,bj) *  
431       &       delZ(1) / work2(I,J)          bosol(I,J) = gravity * TTALPHA(I,J,1) * Qsw(i,j,bi,bj) *
432         &       recip_Cp*recip_rhoConst
433         &       / work2(I,J)
434    
435         END DO         END DO
436        END DO        END DO
437    CML#ifdef ALLOW_SHELFICE
438    CML      ENDIF
439    CML#endif /* ALLOW_SHELFICE */
440          
441    cph(
442    CADJ store ustar = comlev1_kpp, key = ikppkey
443    cph)
444    
445  c------------------------------------------------------------------------  c------------------------------------------------------------------------
446  c     velocity shear  c     velocity shear
# Line 356  c--------------------------------------- Line 453  c---------------------------------------
453    
454  c initialize arrays to zero  c initialize arrays to zero
455        DO k = 1, Nr        DO k = 1, Nr
456           DO j = jbot, jtop           DO j = 1-OLy, sNy+OLy
457              DO i = ibot, itop              DO i = 1-OLx, sNx+OLx
458                 shsq(i,j,k) = p0                 shsq(i,j,k) = p0
459                 dVsq(i,j,k) = p0                 dVsq(i,j,k) = p0
460              END DO              END DO
# Line 384  c     dB/dz exceeds 5.2e-5 s^-2. Line 481  c     dB/dz exceeds 5.2e-5 s^-2.
481              work1(i,j) = nzmax(i,j,bi,bj)              work1(i,j) = nzmax(i,j,bi,bj)
482              DO k = 1, Nr              DO k = 1, Nr
483                 IF ( k .LT. nzmax(i,j,bi,bj) .AND.                 IF ( k .LT. nzmax(i,j,bi,bj) .AND.
484         &              maskC(I,J,k,bi,bj) .GT. 0. .AND.
485       &              dbloc(i,j,k) / drC(k+1) .GT. dB_dz )       &              dbloc(i,j,k) / drC(k+1) .GT. dB_dz )
486       &              work1(i,j) = k       &              work1(i,j) = k
487              END DO              END DO
# Line 436  c     Estimate reference velocity uRef a Line 534  c     Estimate reference velocity uRef a
534                 vRef(i,j) = p5 *                 vRef(i,j) = p5 *
535       &                     ( vVel(i,j,1,bi,bj) + vVel(i,jp1,1,bi,bj) )       &                     ( vVel(i,j,1,bi,bj) + vVel(i,jp1,1,bi,bj) )
536                 IF ( zRef(i,j) .LT. drF(1) ) THEN                 IF ( zRef(i,j) .LT. drF(1) ) THEN
537                    ustarX = ( SurfaceTendencyU(i,  j,bi,bj) +                    ustarX = ( surfaceForcingU(i,  j,bi,bj) +
538       &                       SurfaceTendencyU(ip1,j,bi,bj) ) * p5       &                       surfaceForcingU(ip1,j,bi,bj) ) * p5
539                    ustarY = ( SurfaceTendencyV(i,j,  bi,bj) +       &                   *recip_drF(1)
540       &                       SurfaceTendencyU(i,jp1,bi,bj) ) * p5                    ustarY = ( surfaceForcingV(i,j,  bi,bj) +
541         &                       surfaceForcingV(i,jp1,bi,bj) ) * p5
542         &                   *recip_drF(1)
543                    tempVar1 = ustarX * ustarX + ustarY * ustarY                    tempVar1 = ustarX * ustarX + ustarY * ustarY
544                    if ( tempVar1 .lt. (epsln*epsln) ) then                    if ( tempVar1 .lt. (epsln*epsln) ) then
545                       tempVar2 = epsln                       tempVar2 = epsln
# Line 580  c     shsq computation Line 680  c     shsq computation
680           END DO           END DO
681        END DO        END DO
682    
683    cph(
684    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
685    CADJ store dvsq, shsq = comlev1_kpp, key = ikppkey
686    #endif
687    cph)
688    
689  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
690  c     solve for viscosity, diffusivity, ghat, and hbl on "t-grid"  c     solve for viscosity, diffusivity, ghat, and hbl on "t-grid"
691  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
692    
693        DO j = jbot, jtop  c     precompute background vertical diffusivities, which are needed for
694           DO i = ibot, itop  c     matching diffusivities at bottom of KPP PBL
695          CALL CALC_3D_DIFFUSIVITY(
696         I        bi,bj,1-Olx,sNx+OLx,1-Oly,sNy+OLy,
697         I        GAD_SALINITY, .FALSE., .FALSE.,
698         O        KPPdiffKzS(1-Olx,1-Oly,1,bi,bj),
699         I        myThid)
700          CALL CALC_3D_DIFFUSIVITY(
701         I        bi,bj,1-Olx,sNx+OLx,1-Oly,sNy+OLy,
702         I        GAD_TEMPERATURE, .FALSE., .FALSE.,
703         O        KPPdiffKzT(1-Olx,1-Oly,1,bi,bj),
704         I        myThid)
705    
706          DO j = 1-OLy, sNy+OLy
707             DO i = 1-OLx, sNx+OLx
708              work1(i,j) = nzmax(i,j,bi,bj)              work1(i,j) = nzmax(i,j,bi,bj)
709              work2(i,j) = Fcori(i,j,bi,bj)              work2(i,j) = Fcori(i,j,bi,bj)
710           END DO           END DO
# Line 595  c--------------------------------------- Line 714  c---------------------------------------
714       I       mytime, mythid       I       mytime, mythid
715       I     , work1, shsq, dVsq, ustar       I     , work1, shsq, dVsq, ustar
716       I     , bo, bosol, dbloc, Ritop, work2       I     , bo, bosol, dbloc, Ritop, work2
717       I     , ikey       I     , KPPdiffKzS(1-Olx,1-Oly,1,bi,bj)
718         I     , KPPdiffKzT(1-Olx,1-Oly,1,bi,bj)
719         I     , ikppkey
720       O     , vddiff       O     , vddiff
721       U     , ghat       U     , ghat
722       O     , hbl )       O     , hbl )
   
723        CALL TIMER_STOP ('KPPMIX [KPP_CALC]', myThid)        CALL TIMER_STOP ('KPPMIX [KPP_CALC]', myThid)
724    
 #ifdef ALLOW_AUTODIFF_TAMC  
 cph( storing not necessary  
 cphCADJ STORE vddiff, ghat  = comlev1_kpp, key = ikey  
 cph)  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
725  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
726  c     zero out land values and transfer to global variables  c     zero out land values and transfer to global variables
727  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 615  c--------------------------------------- Line 729  c---------------------------------------
729        DO j = jmin, jmax        DO j = jmin, jmax
730         DO i = imin, imax         DO i = imin, imax
731          DO k = 1, Nr          DO k = 1, Nr
732           KPPviscAz(i,j,k,bi,bj) = vddiff(i,j,k-1,1) * pMask(i,j,k,bi,bj)           km1 = max(1,k-1)
733           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)
734           KPPdiffKzT(i,j,k,bi,bj)= vddiff(i,j,k-1,3) * pMask(i,j,k,bi,bj)       &        * maskC(i,j,km1,bi,bj)
735           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)
736         &        * maskC(i,j,km1,bi,bj)
737             KPPdiffKzT(i,j,k,bi,bj)= vddiff(i,j,k-1,3) * maskC(i,j,k,bi,bj)
738         &        * maskC(i,j,km1,bi,bj)
739             KPPghat(i,j,k,bi,bj)   = ghat(i,j,k)       * maskC(i,j,k,bi,bj)
740         &        * maskC(i,j,km1,bi,bj)
741          END DO          END DO
742          KPPhbl(i,j,bi,bj) = hbl(i,j) * pMask(i,j,1,bi,bj)          k = 1
743    #ifdef ALLOW_SHELFICE
744            if ( useShelfIce ) k = kTopC(i,j,bi,bj)
745    #endif /* ALLOW_SHELFICE */
746            KPPhbl(i,j,bi,bj) = hbl(i,j) * maskC(i,j,k,bi,bj)
747    
748         END DO         END DO
749        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  
750    
751  #ifdef KPP_SMOOTH_VISC  #ifdef KPP_SMOOTH_VISC
752  c     horizontal smoothing of vertical viscosity  c     horizontal smoothing of vertical viscosity
# Line 655  c     horizontal smoothing of vertical d Line 772  c     horizontal smoothing of vertical d
772        _EXCH_XYZ_R8(KPPdiffKzT , myThid )        _EXCH_XYZ_R8(KPPdiffKzT , myThid )
773  #endif /* KPP_SMOOTH_DIFF */  #endif /* KPP_SMOOTH_DIFF */
774    
775    cph(
776    cph  crucial: this avoids full recomp./call of kppmix
777    CADJ store KPPhbl = comlev1_kpp, key = ikppkey
778    cph)
779    
780  C     Compute fraction of solar short-wave flux penetrating to  C     Compute fraction of solar short-wave flux penetrating to
781  C     the bottom of the mixing layer.  C     the bottom of the mixing layer.
782        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
# Line 674  C     the bottom of the mixing layer. Line 796  C     the bottom of the mixing layer.
796    
797        ENDIF        ENDIF
798    
799  #endif ALLOW_KPP  #endif /* ALLOW_KPP */
800    
801        RETURN        RETURN
802        END        END
803    
804          subroutine KPP_CALC_DUMMY(
805         I     bi, bj, myTime, myThid )
806    C     /==========================================================\
807    C     | SUBROUTINE KPP_CALC_DUMMY                                |
808    C     | o Compute all KPP fields defined in KPP.h                |
809    C     | o Dummy routine for TAMC
810    C     |==========================================================|
811    C     | This subroutine serves as an interface between MITGCMUV  |
812    C     | code and NCOM 1-D routines in kpp_routines.F             |
813    C     \==========================================================/
814          IMPLICIT NONE
815    
816    #include "SIZE.h"
817    #include "EEPARAMS.h"
818    #include "PARAMS.h"
819    #include "KPP.h"
820    #include "KPP_PARAMS.h"
821    #include "GRID.h"
822    #include "GAD.h"
823    
824    c Routine arguments
825    c     bi, bj - array indices on which to apply calculations
826    c     myTime - Current time in simulation
827    
828          INTEGER bi, bj
829          INTEGER myThid
830          _RL     myTime
831    
832    #ifdef ALLOW_KPP
833    
834    c Local constants
835          integer i, j, k
836    
837          DO j=1-OLy,sNy+OLy
838             DO i=1-OLx,sNx+OLx
839                KPPhbl (i,j,bi,bj) = 1.0
840                KPPfrac(i,j,bi,bj) = 0.0
841                DO k = 1,Nr
842                   KPPghat   (i,j,k,bi,bj) = 0.0
843                   KPPviscAz (i,j,k,bi,bj) = viscAr
844                ENDDO
845             ENDDO
846          ENDDO
847    
848          CALL CALC_3D_DIFFUSIVITY(
849         I     bi,bj,1-Olx,sNx+OLx,1-Oly,sNy+OLy,
850         I     GAD_SALINITY, .FALSE., .FALSE.,
851         O     KPPdiffKzS(1-Olx,1-Oly,1,bi,bj),
852         I     myThid)
853          CALL CALC_3D_DIFFUSIVITY(
854         I     bi,bj,1-Olx,sNx+OLx,1-Oly,sNy+OLy,
855         I     GAD_TEMPERATURE, .FALSE., .FALSE.,
856         O     KPPdiffKzT(1-Olx,1-Oly,1,bi,bj),
857         I     myThid)
858    
859    #endif
860          RETURN
861          END
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.39
  ViewVC Help
Powered by ViewVC 1.1.22