/[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.12 by mlosch, Wed Sep 25 19:36:50 2002 UTC revision 1.29 by jmc, Sun Jul 18 15:29:52 2004 UTC
# Line 3  C $Name$ Line 3  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 84  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     When option FRUGAL_KPP is used, computation in overlap regions
97  c     is replaced with exchange calls hence reducing overlap requirements  c     is replaced with exchange calls hence reducing overlap requirements
98  c     to OLx=1, OLy=1.  c     to OLx=1, OLy=1.
99    c \ev
100    
101    C !USES: ===============================================================
102  #include "SIZE.h"  #include "SIZE.h"
103  #include "EEPARAMS.h"  #include "EEPARAMS.h"
104  #include "PARAMS.h"  #include "PARAMS.h"
# Line 99  c     to OLx=1, OLy=1. Line 107  c     to OLx=1, OLy=1.
107  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
108  #include "FFIELDS.h"  #include "FFIELDS.h"
109  #include "GRID.h"  #include "GRID.h"
   
110  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
111  #include "tamc.h"  #include "tamc.h"
112  #include "tamc_keys.h"  #include "tamc_keys.h"
113  #else /* ALLOW_AUTODIFF_TAMC */  #else /* ALLOW_AUTODIFF_TAMC */
114        integer ikey        integer ikppkey
115  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
116    
117        EXTERNAL DIFFERENT_MULTIPLE        EXTERNAL DIFFERENT_MULTIPLE
118        LOGICAL  DIFFERENT_MULTIPLE        LOGICAL  DIFFERENT_MULTIPLE
119    
120    C !INPUT PARAMETERS: ===================================================
121  c Routine arguments  c Routine arguments
122  c     bi, bj - array indices on which to apply calculations  c     bi, bj - array indices on which to apply calculations
123  c     myTime - Current time in simulation  c     myTime - Current time in simulation
# Line 120  c     myTime - Current time in simulatio Line 128  c     myTime - Current time in simulatio
128    
129  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
130    
131    C !LOCAL VARIABLES: ====================================================
132  c Local constants  c Local constants
133  c     minusone, p0, p5, p25, p125, p0625  c     minusone, p0, p5, p25, p125, p0625
134  c     imin, imax, jmin, jmax  - array computation indices  c     imin, imax, jmin, jmax  - array computation indices
# Line 128  c     imin, imax, jmin, jmax  - array co Line 137  c     imin, imax, jmin, jmax  - array co
137        parameter( minusone=-1.0)        parameter( minusone=-1.0)
138        _KPP_RL    p0    , p5    , p25     , p125      , p0625        _KPP_RL    p0    , p5    , p25     , p125      , p0625
139        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 )
140        integer    imin      , imax        , jmin      , jmax        integer   imin      ,imax          ,jmin      ,jmax
141  #ifdef FRUGAL_KPP  #ifdef FRUGAL_KPP
142        parameter( imin=1    , imax=sNx    , jmin=1    , jmax=sNy     )        parameter(imin=1    ,imax=sNx      ,jmin=1    ,jmax=sNy      )
143  #else  #else
144        parameter( imin=-2   , imax=sNx+3  , jmin=-2   , jmax=sNy+3   )        parameter(imin=2-OLx,imax=sNx+OLx-1,jmin=2-OLy,jmax=sNy+OLy-1)
145  #endif  #endif
146    
147  c Local arrays and variables  c Local arrays and variables
# Line 149  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 180  c     vRef   (nx,ny)       - Reference m Line 189  c     vRef   (nx,ny)       - Reference m
189        _KPP_RL vRef  ( ibot:itop    , jbot:jtop                    )        _KPP_RL vRef  ( ibot:itop    , jbot:jtop                    )
190  #endif /* KPP_ESTIMATE_UREF */  #endif /* KPP_ESTIMATE_UREF */
191                
192        _KPP_RL tempvar1, tempvar2        _KPP_RL tempvar2
193        integer i, j, k, kp1, im1, ip1, jm1, jp1        integer i, j, k, kp1, im1, ip1, jm1, jp1
194    
195  #ifdef KPP_ESTIMATE_UREF  #ifdef KPP_ESTIMATE_UREF
196        _KPP_RL dBdz1, dBdz2, ustarX, ustarY        _KPP_RL tempvar1, dBdz1, dBdz2, ustarX, ustarY
197  #endif  #endif
198    
199    #ifdef ALLOW_AUTODIFF_TAMC
200              act1 = bi - myBxLo(myThid)
201              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
202              act2 = bj - myByLo(myThid)
203              max2 = myByHi(myThid) - myByLo(myThid) + 1
204              act3 = myThid - 1
205              max3 = nTx*nTy
206              act4 = ikey_dynamics - 1
207              ikppkey = (act1 + 1) + act2*max1
208         &                      + act3*max1*max2
209         &                      + act4*max1*max2*max3
210    #endif /* ALLOW_AUTODIFF_TAMC */
211    CEOP
212    
213  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?
214        IF ( DIFFERENT_MULTIPLE(kpp_freq,myTime,myTime-deltaTClock) .OR.        IF ( DIFFERENT_MULTIPLE(kpp_freq,myTime,myTime-deltaTClock) .OR.
215       1     myTime .EQ. startTime ) THEN       1     myTime .EQ. startTime ) THEN
# Line 227  c--------------------------------------- Line 250  c---------------------------------------
250    
251        CALL TIMER_START('STATEKPP      [KPP_CALC]', myThid)        CALL TIMER_START('STATEKPP      [KPP_CALC]', myThid)
252        CALL STATEKPP(        CALL STATEKPP(
253       I       bi, bj, myThid       I       ikppkey, bi, bj, myThid
254       O     , work2, dbloc, Ritop       O     , work2, dbloc, Ritop
255       O     , vddiff(ibot,jbot,1,1), vddiff(ibot,jbot,1,2)       O     , vddiff(ibot,jbot,1,1), vddiff(ibot,jbot,1,2)
256       &     )       &     )
# Line 282  c     horizontally smooth density relate Line 305  c     horizontally smooth density relate
305    
306  c     zero out dbloc over land points (so that the convective  c     zero out dbloc over land points (so that the convective
307  c     part of the interior mixing can be diagnosed)  c     part of the interior mixing can be diagnosed)
308                 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)
309                 ghat(i,j,k)  = ghat(i,j,k)  * pMask(i,j,k,bi,bj)                 ghat(i,j,k)  = ghat(i,j,k)  * maskC(i,j,k,bi,bj)
310                 Ritop(i,j,k) = Ritop(i,j,k) * pMask(i,j,k,bi,bj)                 Ritop(i,j,k) = Ritop(i,j,k) * maskC(i,j,k,bi,bj)
311                 if(k.eq.nzmax(i,j,bi,bj)) then                 if(k.eq.nzmax(i,j,bi,bj)) then
312                    dbloc(i,j,k) = p0                    dbloc(i,j,k) = p0
313                    ghat(i,j,k)  = p0                    ghat(i,j,k)  = p0
# Line 302  c     so that the subroutine "bldepth" w Line 325  c     so that the subroutine "bldepth" w
325    
326  cph(  cph(
327  cph  this avoids a single or double recomp./call of statekpp  cph  this avoids a single or double recomp./call of statekpp
328  CADJ store work2              = comlev1_kpp, key = ikey  CADJ store work2              = comlev1_kpp, key = ikppkey
329  #ifdef ALLOW_AUTODIFF_KPP_EXTENSIVE_STORE  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
330  CADJ store dbloc, Ritop, ghat = comlev1_kpp, key = ikey  CADJ store dbloc, Ritop, ghat = comlev1_kpp, key = ikppkey
331  CADJ store vddiff             = comlev1_kpp, key = ikey  CADJ store vddiff             = comlev1_kpp, key = ikppkey
332  #endif  #endif
333  cph)  cph)
334    
335  c------------------------------------------------------------------------  c------------------------------------------------------------------------
336  c     friction velocity, turbulent and radiative surface buoyancy forcing  c     friction velocity, turbulent and radiative surface buoyancy forcing
337  c     -------------------------------------------------------------------  c     -------------------------------------------------------------------
338  c     taux / rho = SurfaceTendencyU * drF(1)                     (N/m^2)  c     taux / rho = surfaceForcingU                               (N/m^2)
339  c     tauy / rho = SurfaceTendencyV * drF(1)                     (N/m^2)  c     tauy / rho = surfaceForcingV                               (N/m^2)
340  c     ustar = sqrt( sqrt( taux^2 + tauy^2 ) / rho )                (m/s)  c     ustar = sqrt( sqrt( taux^2 + tauy^2 ) / rho )                (m/s)
341  c     bo    = - g * ( alpha*SurfaceTendencyT +  c     bo    = - g * ( alpha*surfaceForcingT +
342  c                     beta *SurfaceTendencyS ) * drF(1) / rho  (m^2/s^3)  c                     beta *surfaceForcingS ) / rho            (m^2/s^3)
343  c     bosol = - g * alpha * Qsw * drF(1) / rho                 (m^2/s^3)  c     bosol = - g * alpha * Qsw * drF(1) / rho                 (m^2/s^3)
344  c------------------------------------------------------------------------  c------------------------------------------------------------------------
345    
# Line 334  c initialize arrays to zero Line 357  c initialize arrays to zero
357         DO i = imin, imax         DO i = imin, imax
358          ip1 = i+1          ip1 = i+1
359          work3(i,j) =          work3(i,j) =
360       &   (SurfaceTendencyU(i,j,bi,bj) + SurfaceTendencyU(ip1,j,bi,bj)) *       &   (surfaceForcingU(i,j,bi,bj) + surfaceForcingU(ip1,j,bi,bj)) *
361       &   (SurfaceTendencyU(i,j,bi,bj) + SurfaceTendencyU(ip1,j,bi,bj)) +       &   (surfaceForcingU(i,j,bi,bj) + surfaceForcingU(ip1,j,bi,bj)) +
362       &   (SurfaceTendencyV(i,j,bi,bj) + SurfaceTendencyV(i,jp1,bi,bj)) *       &   (surfaceForcingV(i,j,bi,bj) + surfaceForcingV(i,jp1,bi,bj)) *
363       &   (SurfaceTendencyV(i,j,bi,bj) + SurfaceTendencyV(i,jp1,bi,bj))       &   (surfaceForcingV(i,j,bi,bj) + surfaceForcingV(i,jp1,bi,bj))
364         END DO         END DO
365        END DO        END DO
366  cph(  cph(
367  CADJ store work3 = comlev1_kpp, key = ikey  CADJ store work3 = comlev1_kpp, key = ikppkey
368  cph)  cph)
369        DO j = jmin, jmax        DO j = jmin, jmax
370         jp1 = j + 1         jp1 = j + 1
371         DO i = imin, imax         DO i = imin, imax
372          ip1 = i+1          ip1 = i+1
373          if ( work3(i,j) .lt. (phepsi*phepsi) ) then  
374            if ( work3(i,j) .lt. (phepsi*phepsi*drF(1)*drF(1)) ) then
375             ustar(i,j) = SQRT( phepsi * p5 * drF(1) )             ustar(i,j) = SQRT( phepsi * p5 * drF(1) )
376          else          else
377             tempVar2 =  SQRT( work3(i,j) ) * p5 * drF(1)             tempVar2 =  SQRT( work3(i,j) ) * p5
378             ustar(i,j) = SQRT( tempVar2 )             ustar(i,j) = SQRT( tempVar2 )
379          endif          endif
380    
381          bo(I,J) = - gravity *          bo(I,J) = - gravity *
382       &       ( vddiff(I,J,1,1) * SurfaceTendencyT(i,j,bi,bj) +       &       ( vddiff(I,J,1,1) * (surfaceForcingT(i,j,bi,bj)+
383       &         vddiff(I,J,1,2) * SurfaceTendencyS(i,j,bi,bj)       &       surfaceForcingTice(i,j,bi,bj)) +
384       &       ) *       &       vddiff(I,J,1,2) * surfaceForcingS(i,j,bi,bj) )
385       &       drF(1) / work2(I,J)       &       / work2(I,J)
386    
387          bosol(I,J) = gravity * vddiff(I,J,1,1) * Qsw(i,j,bi,bj) *          bosol(I,J) = gravity * vddiff(I,J,1,1) * Qsw(i,j,bi,bj) *
388       &       recip_Cp*recip_rhoConst*recip_dRf(1) *       &       recip_Cp*recip_rhoConst
389       &       drF(1) / work2(I,J)       &       / work2(I,J)
390    
391         END DO         END DO
392        END DO        END DO
393    
394  cph(  cph(
395  CADJ store ustar = comlev1_kpp, key = ikey  CADJ store ustar = comlev1_kpp, key = ikppkey
396  cph)  cph)
397    
398  c------------------------------------------------------------------------  c------------------------------------------------------------------------
# Line 459  c     Estimate reference velocity uRef a Line 486  c     Estimate reference velocity uRef a
486                 vRef(i,j) = p5 *                 vRef(i,j) = p5 *
487       &                     ( vVel(i,j,1,bi,bj) + vVel(i,jp1,1,bi,bj) )       &                     ( vVel(i,j,1,bi,bj) + vVel(i,jp1,1,bi,bj) )
488                 IF ( zRef(i,j) .LT. drF(1) ) THEN                 IF ( zRef(i,j) .LT. drF(1) ) THEN
489                    ustarX = ( SurfaceTendencyU(i,  j,bi,bj) +                    ustarX = ( surfaceForcingU(i,  j,bi,bj) +
490       &                       SurfaceTendencyU(ip1,j,bi,bj) ) * p5       &                       surfaceForcingU(ip1,j,bi,bj) ) * p5
491                    ustarY = ( SurfaceTendencyV(i,j,  bi,bj) +       &                   *recip_drF(1)
492       &                       SurfaceTendencyU(i,jp1,bi,bj) ) * p5                    ustarY = ( surfaceForcingV(i,j,  bi,bj) +
493         &                       surfaceForcingV(i,jp1,bi,bj) ) * p5
494         &                   *recip_drF(1)
495                    tempVar1 = ustarX * ustarX + ustarY * ustarY                    tempVar1 = ustarX * ustarX + ustarY * ustarY
496                    if ( tempVar1 .lt. (epsln*epsln) ) then                    if ( tempVar1 .lt. (epsln*epsln) ) then
497                       tempVar2 = epsln                       tempVar2 = epsln
# Line 604  c     shsq computation Line 633  c     shsq computation
633        END DO        END DO
634    
635  cph(  cph(
636  #ifdef ALLOW_AUTODIFF_KPP_EXTENSIVE_STORE  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
637  CADJ store dvsq, shsq = comlev1_kpp, key = ikey  CADJ store dvsq, shsq = comlev1_kpp, key = ikppkey
638  #endif  #endif
639  cph)  cph)
640    
# Line 624  c--------------------------------------- Line 653  c---------------------------------------
653       I       mytime, mythid       I       mytime, mythid
654       I     , work1, shsq, dVsq, ustar       I     , work1, shsq, dVsq, ustar
655       I     , bo, bosol, dbloc, Ritop, work2       I     , bo, bosol, dbloc, Ritop, work2
656       I     , ikey       I     , ikppkey
657       O     , vddiff       O     , vddiff
658       U     , ghat       U     , ghat
659       O     , hbl )       O     , hbl )
# Line 638  c--------------------------------------- Line 667  c---------------------------------------
667        DO j = jmin, jmax        DO j = jmin, jmax
668         DO i = imin, imax         DO i = imin, imax
669          DO k = 1, Nr          DO k = 1, Nr
670           KPPviscAz(i,j,k,bi,bj) = vddiff(i,j,k-1,1) * pMask(i,j,k,bi,bj)           KPPviscAz(i,j,k,bi,bj) = vddiff(i,j,k-1,1) * maskC(i,j,k,bi,bj)
671           KPPdiffKzS(i,j,k,bi,bj)= vddiff(i,j,k-1,2) * pMask(i,j,k,bi,bj)           KPPdiffKzS(i,j,k,bi,bj)= vddiff(i,j,k-1,2) * maskC(i,j,k,bi,bj)
672           KPPdiffKzT(i,j,k,bi,bj)= vddiff(i,j,k-1,3) * pMask(i,j,k,bi,bj)           KPPdiffKzT(i,j,k,bi,bj)= vddiff(i,j,k-1,3) * maskC(i,j,k,bi,bj)
673           KPPghat(i,j,k,bi,bj)   = ghat(i,j,k)       * pMask(i,j,k,bi,bj)           KPPghat(i,j,k,bi,bj)   = ghat(i,j,k)       * maskC(i,j,k,bi,bj)
674          END DO          END DO
675          KPPhbl(i,j,bi,bj) = hbl(i,j) * pMask(i,j,1,bi,bj)          KPPhbl(i,j,bi,bj) = hbl(i,j) * maskC(i,j,1,bi,bj)
676         END DO         END DO
677        END DO        END DO
678  #ifdef FRUGAL_KPP  #ifdef FRUGAL_KPP
# Line 680  c     horizontal smoothing of vertical d Line 709  c     horizontal smoothing of vertical d
709    
710  cph(  cph(
711  cph  crucial: this avoids full recomp./call of kppmix  cph  crucial: this avoids full recomp./call of kppmix
712  CADJ store KPPhbl = comlev1_kpp, key = ikey  CADJ store KPPhbl = comlev1_kpp, key = ikppkey
713  cph)  cph)
714    
715  C     Compute fraction of solar short-wave flux penetrating to  C     Compute fraction of solar short-wave flux penetrating to
# Line 745  c Local constants Line 774  c Local constants
774              KPPfrac(i,j,bi,bj) = 0.0              KPPfrac(i,j,bi,bj) = 0.0
775              DO k = 1,Nr              DO k = 1,Nr
776                 KPPghat   (i,j,k,bi,bj) = 0.0                 KPPghat   (i,j,k,bi,bj) = 0.0
777                 KPPviscAz (i,j,k,bi,bj) = viscAz                 KPPviscAz (i,j,k,bi,bj) = viscAr
778                 KPPdiffKzT(i,j,k,bi,bj) = diffKzT                 KPPdiffKzT(i,j,k,bi,bj) = diffKrT
779                 KPPdiffKzS(i,j,k,bi,bj) = diffKzS                 KPPdiffKzS(i,j,k,bi,bj) = diffKrS
780              ENDDO              ENDDO
781           ENDDO           ENDDO
782        ENDDO        ENDDO
Legend:
Removed from v.1.12  
changed lines
  Added in v.1.29
  ViewVC Help
Powered by ViewVC 1.1.22