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

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

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

revision 1.20 by jmc, Sun Jul 18 01:18:55 2004 UTC revision 1.29 by jmc, Mon Apr 30 13:49:40 2007 UTC
# Line 11  C--  o BLDEPTH     - Determine oceanic p Line 11  C--  o BLDEPTH     - Determine oceanic p
11  C--  o WSCALE      - Compute turbulent velocity scales.  C--  o WSCALE      - Compute turbulent velocity scales.
12  C--  o RI_IWMIX    - Compute interior viscosity diffusivity coefficients.  C--  o RI_IWMIX    - Compute interior viscosity diffusivity coefficients.
13  C--  o Z121        - Apply 121 vertical smoothing.  C--  o Z121        - Apply 121 vertical smoothing.
14  C--  o KPP_SMOOTH_HORIZ - Apply horizontal smoothing to KPP array.  C--  o SMOOTH_HORIZ- Apply horizontal smoothing to global array.
 C--  o SMOOTH_HORIZ     - Apply horizontal smoothing to global array.  
15  C--  o BLMIX       - Boundary layer mixing coefficients.  C--  o BLMIX       - Boundary layer mixing coefficients.
16  C--  o ENHANCE     - Enhance diffusivity at boundary layer interface.  C--  o ENHANCE     - Enhance diffusivity at boundary layer interface.
17  C--  o STATEKPP    - Compute buoyancy-related input arrays.  C--  o STATEKPP    - Compute buoyancy-related input arrays.
# Line 23  c*************************************** Line 22  c***************************************
22       I       mytime, mythid       I       mytime, mythid
23       I     , kmtj, shsq, dvsq, ustar       I     , kmtj, shsq, dvsq, ustar
24       I     , bo, bosol, dbloc, Ritop, coriol       I     , bo, bosol, dbloc, Ritop, coriol
25         I     , diffusKzS, diffusKzT
26       I     , ikppkey       I     , ikppkey
27       O     , diffus       O     , diffus
28       U     , ghat       U     , ghat
# Line 47  c--------------------------------------- Line 47  c---------------------------------------
47  #include "SIZE.h"  #include "SIZE.h"
48  #include "EEPARAMS.h"  #include "EEPARAMS.h"
49  #include "PARAMS.h"  #include "PARAMS.h"
 #include "DYNVARS.h"  
 #include "FFIELDS.h"  
50  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
51    
52  c input  c input
53  c     myTime          - current time in simulation  c     myTime           - current time in simulation
54  c     myThid          - thread number for this instance of the routine  c     myThid           - thread number for this instance of the routine
55  c     kmtj   (imt)    - number of vertical layers on this row  c     kmtj   (imt)     - number of vertical layers on this row
56  c     shsq   (imt,Nr) - (local velocity shear)^2                      ((m/s)^2)  c     shsq   (imt,Nr)  - (local velocity shear)^2                     ((m/s)^2)
57  c     dvsq   (imt,Nr) - (velocity shear re sfc)^2                     ((m/s)^2)  c     dvsq   (imt,Nr)  - (velocity shear re sfc)^2                    ((m/s)^2)
58  c     ustar  (imt)    - surface friction velocity                         (m/s)  c     ustar  (imt)     - surface friction velocity                        (m/s)
59  c     bo     (imt)    - surface turbulent buoy. forcing               (m^2/s^3)  c     bo     (imt)     - surface turbulent buoy. forcing              (m^2/s^3)
60  c     bosol  (imt)    - radiative buoyancy forcing                    (m^2/s^3)  c     bosol  (imt)     - radiative buoyancy forcing                   (m^2/s^3)
61  c     dbloc  (imt,Nr) - local delta buoyancy across interfaces          (m/s^2)  c     dbloc  (imt,Nr)  - local delta buoyancy across interfaces         (m/s^2)
62  c     dblocSm(imt,Nr) - horizontally smoothed dbloc                     (m/s^2)  c     dblocSm(imt,Nr)  - horizontally smoothed dbloc                    (m/s^2)
63  c                         stored in ghat to save space  c                          stored in ghat to save space
64  c     Ritop  (imt,Nr) - numerator of bulk Richardson Number  c     Ritop  (imt,Nr)  - numerator of bulk Richardson Number
65  c                         (zref-z) * delta buoyancy w.r.t. surface    ((m/s)^2)  c                          (zref-z) * delta buoyancy w.r.t. surface   ((m/s)^2)
66  c     coriol (imt)    - Coriolis parameter                                (1/s)  c     coriol (imt)     - Coriolis parameter                               (1/s)
67    c     diffusKzS(imt,Nr)- background vertical diffusivity for scalars    (m^2/s)
68    c     diffusKzT(imt,Nr)- background vertical diffusivity for theta      (m^2/s)
69  c     note: there is a conversion from 2-D to 1-D for input output variables,  c     note: there is a conversion from 2-D to 1-D for input output variables,
70  c           e.g., hbl(sNx,sNy) -> hbl(imt),  c           e.g., hbl(sNx,sNy) -> hbl(imt),
71  c           where hbl(i,j) -> hbl((j-1)*sNx+i)  c           where hbl(i,j) -> hbl((j-1)*sNx+i)
72    
73        _RL     mytime        _RL     mytime
74        integer mythid        integer mythid
75        integer kmtj  (imt   )        integer kmtj (imt   )
76        _KPP_RL shsq  (imt,Nr)        _RL shsq     (imt,Nr)
77        _KPP_RL dvsq  (imt,Nr)        _RL dvsq     (imt,Nr)
78        _KPP_RL ustar (imt   )        _RL ustar    (imt   )
79        _KPP_RL bo    (imt   )        _RL bo       (imt   )
80        _KPP_RL bosol (imt   )        _RL bosol    (imt   )
81        _KPP_RL dbloc (imt,Nr)        _RL dbloc    (imt,Nr)
82        _KPP_RL Ritop (imt,Nr)        _RL Ritop    (imt,Nr)
83        _KPP_RL coriol(imt   )        _RL coriol   (imt   )
84          _RL diffusKzS(imt,Nr)
85          _RL diffusKzT(imt,Nr)
86    
87        integer ikppkey        integer ikppkey
88    
# Line 91  c     diffus (imt,3)  - vertical tempera Line 93  c     diffus (imt,3)  - vertical tempera
93  c     ghat   (imt)    - nonlocal transport coefficient                  (s/m^2)  c     ghat   (imt)    - nonlocal transport coefficient                  (s/m^2)
94  c     hbl    (imt)    - mixing layer depth                                  (m)  c     hbl    (imt)    - mixing layer depth                                  (m)
95    
96        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL diffus(imt,0:Nrp1,mdiff)
97        _KPP_RL ghat  (imt,Nr)        _RL ghat  (imt,Nr)
98        _KPP_RL hbl   (imt)        _RL hbl   (imt)
99    
100  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
101    
# Line 107  c     blmc   (imt,Nr,mdiff) - boundary l Line 109  c     blmc   (imt,Nr,mdiff) - boundary l
109  c     sigma  (imt         ) - normalized depth (d / hbl)  c     sigma  (imt         ) - normalized depth (d / hbl)
110  c     Rib    (imt,Nr      ) - bulk Richardson number  c     Rib    (imt,Nr      ) - bulk Richardson number
111    
112        integer kbl   (imt         )        integer kbl(imt         )
113        _KPP_RL bfsfc (imt         )        _RL bfsfc  (imt         )
114        _KPP_RL casea (imt         )        _RL casea  (imt         )
115        _KPP_RL stable(imt         )        _RL stable (imt         )
116        _KPP_RL dkm1  (imt,   mdiff)        _RL dkm1   (imt,   mdiff)
117        _KPP_RL blmc  (imt,Nr,mdiff)        _RL blmc   (imt,Nr,mdiff)
118        _KPP_RL sigma (imt         )        _RL sigma  (imt         )
119        _KPP_RL Rib   (imt,Nr      )        _RL Rib    (imt,Nr      )
120    
121        integer i, k, md        integer i, k, md
122    
# Line 132  CADJ STORE dbloc = comlev1_kpp, key = ik Line 134  CADJ STORE dbloc = comlev1_kpp, key = ik
134  cph)  cph)
135        call Ri_iwmix (        call Ri_iwmix (
136       I       kmtj, shsq, dbloc, ghat       I       kmtj, shsq, dbloc, ghat
137         I     , diffusKzS, diffusKzT
138       I     , ikppkey       I     , ikppkey
139       O     , diffus )       O     , diffus )
140    
# Line 210  c--------------------------------------- Line 213  c---------------------------------------
213           do i = 1, imt           do i = 1, imt
214              if (k .lt. kbl(i)) then              if (k .lt. kbl(i)) then
215                 diffus(i,k,1) = max ( blmc(i,k,1), viscAr  )                 diffus(i,k,1) = max ( blmc(i,k,1), viscAr  )
216                 diffus(i,k,2) = max ( blmc(i,k,2), diffKrS )                 diffus(i,k,2) = max ( blmc(i,k,2), diffusKzS(i,Nr) )
217                 diffus(i,k,3) = max ( blmc(i,k,3), diffKrT )                 diffus(i,k,3) = max ( blmc(i,k,3), diffusKzT(i,Nr) )
218              else              else
219                 ghat(i,k) = 0.                 ghat(i,k) = 0.
220              endif              endif
# Line 257  c Line 260  c
260        IMPLICIT NONE        IMPLICIT NONE
261    
262  #include "SIZE.h"  #include "SIZE.h"
 #include "EEPARAMS.h"  
 #include "PARAMS.h"  
263  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
 #include "FFIELDS.h"  
264    
265  c input  c input
266  c------  c------
# Line 279  c coriol    : Coriolis parameter Line 279  c coriol    : Coriolis parameter
279        _RL     mytime        _RL     mytime
280        integer mythid        integer mythid
281        integer kmtj(imt)        integer kmtj(imt)
282        _KPP_RL dvsq  (imt,Nr)        _RL dvsq    (imt,Nr)
283        _KPP_RL dbloc (imt,Nr)        _RL dbloc   (imt,Nr)
284        _KPP_RL Ritop (imt,Nr)        _RL Ritop   (imt,Nr)
285        _KPP_RL ustar (imt)        _RL ustar   (imt)
286        _KPP_RL bo    (imt)        _RL bo      (imt)
287        _KPP_RL bosol (imt)        _RL bosol   (imt)
288        _KPP_RL coriol(imt)        _RL coriol  (imt)
289        integer ikppkey, kkppkey        integer ikppkey, kkppkey
290    
291  c  output  c  output
# Line 297  c casea     : =1 in case A, =0 in case B Line 297  c casea     : =1 in case A, =0 in case B
297  c kbl       : -1 of first grid level below hbl  c kbl       : -1 of first grid level below hbl
298  c Rib       : Bulk Richardson number  c Rib       : Bulk Richardson number
299  c sigma     : normalized depth (d/hbl)  c sigma     : normalized depth (d/hbl)
300        _KPP_RL hbl   (imt)        _RL hbl    (imt)
301        _KPP_RL bfsfc (imt)        _RL bfsfc  (imt)
302        _KPP_RL stable(imt)        _RL stable (imt)
303        _KPP_RL casea (imt)        _RL casea  (imt)
304        integer kbl   (imt)        integer kbl(imt)
305        _KPP_RL Rib   (imt,Nr)        _RL Rib    (imt,Nr)
306        _KPP_RL sigma (imt)        _RL sigma  (imt)
307    
308  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
309    
310  c  local  c  local
311  c-------  c-------
312  c wm, ws    : turbulent velocity scales         (m/s)  c wm, ws    : turbulent velocity scales         (m/s)
313        _KPP_RL wm(imt), ws(imt)        _RL wm(imt), ws(imt)
314        _RL     worka(imt)        _RL worka(imt)
315    
316        _KPP_RL bvsq, vtsq, hekman, hmonob, hlimit, tempVar1, tempVar2        _RL bvsq, vtsq, hekman, hmonob, hlimit, tempVar1, tempVar2
317        integer i, kl        integer i, kl
318    
319        _KPP_RL     p5    , eins        _RL         p5    , eins
320        parameter ( p5=0.5, eins=1.0 )        parameter ( p5=0.5, eins=1.0 )
321        _RL         minusone        _RL         minusone
322        parameter ( minusone=-1.0 )        parameter ( minusone=-1.0 )
# Line 583  c sigma   : normalized depth (d/hbl) Line 583  c sigma   : normalized depth (d/hbl)
583  c hbl     : boundary layer depth (m)  c hbl     : boundary layer depth (m)
584  c ustar   : surface friction velocity         (m/s)  c ustar   : surface friction velocity         (m/s)
585  c bfsfc   : total surface buoyancy flux       (m^2/s^3)  c bfsfc   : total surface buoyancy flux       (m^2/s^3)
586        _KPP_RL sigma(imt)        _RL sigma(imt)
587        _KPP_RL hbl  (imt)        _RL hbl  (imt)
588        _KPP_RL ustar(imt)        _RL ustar(imt)
589        _KPP_RL bfsfc(imt)        _RL bfsfc(imt)
590    
591  c  output  c  output
592  c--------  c--------
593  c wm, ws  : turbulent velocity scales at sigma  c wm, ws  : turbulent velocity scales at sigma
594        _KPP_RL wm(imt), ws(imt)        _RL wm(imt), ws(imt)
595    
596  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
597    
598  c local  c local
599  c------  c------
600  c zehat   : = zeta *  ustar**3  c zehat   : = zeta *  ustar**3
601        _KPP_RL zehat        _RL zehat
602    
603        integer iz, izp1, ju, i, jup1        integer iz, izp1, ju, i, jup1
604        _KPP_RL udiff, zdiff, zfrac, ufrac, fzfrac, wam        _RL udiff, zdiff, zfrac, ufrac, fzfrac, wam
605        _KPP_RL wbm, was, wbs, u3, tempVar        _RL wbm, was, wbs, u3, tempVar
606    
607  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
608  c use lookup table for zehat < zmax only; otherwise use  c use lookup table for zehat < zmax only; otherwise use
# Line 658  c*************************************** Line 658  c***************************************
658    
659        subroutine Ri_iwmix (        subroutine Ri_iwmix (
660       I       kmtj, shsq, dbloc, dblocSm       I       kmtj, shsq, dbloc, dblocSm
661         I     , diffusKzS, diffusKzT
662       I     , ikppkey       I     , ikppkey
663       O     , diffus )       O     , diffus )
664    
# Line 678  c     kmtj   (imt)         number of ver Line 679  c     kmtj   (imt)         number of ver
679  c     shsq   (imt,Nr)      (local velocity shear)^2               ((m/s)^2)  c     shsq   (imt,Nr)      (local velocity shear)^2               ((m/s)^2)
680  c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)  c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)
681  c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)  c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)
682        integer kmtj   (imt)  c     diffusKzS(imt,Nr)- background vertical diffusivity for scalars    (m^2/s)
683        _KPP_RL shsq   (imt,Nr)  c     diffusKzT(imt,Nr)- background vertical diffusivity for theta      (m^2/s)
684        _KPP_RL dbloc  (imt,Nr)        integer kmtj (imt)
685        _KPP_RL dblocSm(imt,Nr)        _RL shsq     (imt,Nr)
686          _RL dbloc    (imt,Nr)
687          _RL dblocSm  (imt,Nr)
688          _RL diffusKzS(imt,Nr)
689          _RL diffusKzT(imt,Nr)
690        integer ikppkey        integer ikppkey
691    
692  c output  c output
693  c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)  c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)
694  c     diffus(imt,0:Nrp1,2)  vertical scalar diffusivity             (m^2/s)  c     diffus(imt,0:Nrp1,2)  vertical scalar diffusivity             (m^2/s)
695  c     diffus(imt,0:Nrp1,3)  vertical temperature diffusivity        (m^2/s)  c     diffus(imt,0:Nrp1,3)  vertical temperature diffusivity        (m^2/s)
696        _KPP_RL diffus(imt,0:Nrp1,3)        _RL diffus(imt,0:Nrp1,3)
697    
698  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
699    
700  c local variables  c local variables
701  c     Rig                   local Richardson number  c     Rig                   local Richardson number
702  c     fRi, fcon             function of Rig  c     fRi, fcon             function of Rig
703        _KPP_RL Rig        _RL Rig
704        _KPP_RL fRi, fcon        _RL fRi, fcon
705        _KPP_RL ratio        _RL ratio
706        integer i, ki        integer i, ki
707        _KPP_RL c1, c0        _RL c1, c0
708    
709  #ifdef ALLOW_KPP_VERTICALLY_SMOOTH  #ifdef ALLOW_KPP_VERTICALLY_SMOOTH
710        integer mr        integer mr
# Line 784  c            evaluate diffusivities and Line 789  c            evaluate diffusivities and
789  c    mixing due to internal waves, and shear and static instability  c    mixing due to internal waves, and shear and static instability
790    
791  #ifndef EXCLUDE_KPP_SHEAR_MIX  #ifndef EXCLUDE_KPP_SHEAR_MIX
792              diffus(i,ki,1) = viscAr + fcon * difmcon + fRi * difm0              if ( .NOT. inAdMode ) then
793              diffus(i,ki,2) = diffKrS + fcon * difscon + fRi * difs0                 diffus(i,ki,1) = viscAr + fcon * difmcon + fRi * difm0
794              diffus(i,ki,3) = diffKrT + fcon * difscon + fRi * difs0                 diffus(i,ki,2) = diffusKzS(i,ki)+fcon*difscon+fRi*difs0
795                   diffus(i,ki,3) = diffusKzT(i,ki)+fcon*diftcon+fRi*dift0
796                else
797                   diffus(i,ki,1) = viscAr
798                   diffus(i,ki,2) = diffusKzS(i,ki)
799                   diffus(i,ki,3) = diffusKzT(i,ki)
800                endif
801  #else  #else
802              diffus(i,ki,1) = viscAr              diffus(i,ki,1) = viscAr
803              diffus(i,ki,2) = diffKrS              diffus(i,ki,2) = diffusKzS(i,ki)
804              diffus(i,ki,3) = diffKrT              diffus(i,ki,3) = diffusKzT(i,ki)
805  #endif  #endif
806    
807           end do           end do
# Line 831  c     penetrative convection. Line 842  c     penetrative convection.
842  c input/output  c input/output
843  c-------------  c-------------
844  c v     : 2-D array to be smoothed in Nrp1 direction  c v     : 2-D array to be smoothed in Nrp1 direction
845        _KPP_RL v(imt,0:Nrp1)        _RL v(imt,0:Nrp1)
846    
847  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
848    
849  c local  c local
850        _KPP_RL zwork, zflag        _RL zwork, zflag
851        _KPP_RL KRi_range(1:Nrp1)        _RL KRi_range(1:Nrp1)
852        integer i, k, km1, kp1        integer i, k, km1, kp1
853    
854        _KPP_RL     p0      , p25       , p5      , p2        _RL         p0      , p25       , p5      , p2
855        parameter ( p0 = 0.0, p25 = 0.25, p5 = 0.5, p2 = 2.0 )        parameter ( p0 = 0.0, p25 = 0.25, p5 = 0.5, p2 = 2.0 )
856    
857        KRi_range(Nrp1) = p0        KRi_range(Nrp1) = p0
# Line 894  CADJ STORE v(i,k), zwork = z121tape Line 905  CADJ STORE v(i,k), zwork = z121tape
905    
906  c*************************************************************************  c*************************************************************************
907    
       subroutine kpp_smooth_horiz (  
      I     k, bi, bj,  
      U     fld )  
   
 c     Apply horizontal smoothing to KPP array  
   
       IMPLICIT NONE  
 #include "SIZE.h"  
 #include "GRID.h"  
 #include "KPP_PARAMS.h"  
   
 c     input  
 c     bi, bj : array indices  
 c     k      : vertical index used for masking  
       integer k, bi, bj  
   
 c     input/output  
 c     fld    : 2-D array to be smoothed  
       _KPP_RL fld( ibot:itop, jbot:jtop )  
   
 #ifdef ALLOW_KPP  
   
 c     local  
       integer i, j, im1, ip1, jm1, jp1  
       _KPP_RL tempVar  
       _KPP_RL fld_tmp( ibot:itop, jbot:jtop )  
   
       integer    imin       , imax       , jmin       , jmax  
       parameter( imin=ibot+1, imax=itop-1, jmin=jbot+1, jmax=jtop-1 )  
   
       _KPP_RL    p0    , p5    , p25     , p125      , p0625  
       parameter( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 )  
   
       DO j = jmin, jmax  
          jm1 = j-1  
          jp1 = j+1  
          DO i = imin, imax  
             im1 = i-1  
             ip1 = i+1  
             tempVar =  
      &           p25   *   maskC(i  ,j  ,k,bi,bj)   +  
      &           p125  * ( maskC(im1,j  ,k,bi,bj)   +  
      &                     maskC(ip1,j  ,k,bi,bj)   +  
      &                     maskC(i  ,jm1,k,bi,bj)   +  
      &                     maskC(i  ,jp1,k,bi,bj) ) +  
      &           p0625 * ( maskC(im1,jm1,k,bi,bj)   +  
      &                     maskC(im1,jp1,k,bi,bj)   +  
      &                     maskC(ip1,jm1,k,bi,bj)   +  
      &                     maskC(ip1,jp1,k,bi,bj) )  
             IF ( tempVar .GE. p25 ) THEN  
                fld_tmp(i,j) = (  
      &              p25  * fld(i  ,j  )*maskC(i  ,j  ,k,bi,bj) +  
      &              p125 *(fld(im1,j  )*maskC(im1,j  ,k,bi,bj) +  
      &                     fld(ip1,j  )*maskC(ip1,j  ,k,bi,bj) +  
      &                     fld(i  ,jm1)*maskC(i  ,jm1,k,bi,bj) +  
      &                     fld(i  ,jp1)*maskC(i  ,jp1,k,bi,bj))+  
      &              p0625*(fld(im1,jm1)*maskC(im1,jm1,k,bi,bj) +  
      &                     fld(im1,jp1)*maskC(im1,jp1,k,bi,bj) +  
      &                     fld(ip1,jm1)*maskC(ip1,jm1,k,bi,bj) +  
      &                     fld(ip1,jp1)*maskC(ip1,jp1,k,bi,bj)))  
      &              / tempVar  
             ELSE  
                fld_tmp(i,j) = fld(i,j)  
             ENDIF  
          ENDDO  
       ENDDO  
   
 c     transfer smoothed field to output array  
       DO j = jmin, jmax  
          DO i = imin, imax  
             fld(i,j) = fld_tmp(i,j)  
          ENDDO  
       ENDDO  
   
 #endif /* ALLOW_KPP */  
   
       return  
       end  
   
 c*************************************************************************  
   
908        subroutine smooth_horiz (        subroutine smooth_horiz (
909       I     k, bi, bj,       I     k, bi, bj,
910       U     fld )       U     fld )
# Line 1081  c     hbl   (imt)                 bounda Line 1011  c     hbl   (imt)                 bounda
1011  c     stable(imt)                 = 1 in stable forcing  c     stable(imt)                 = 1 in stable forcing
1012  c     casea (imt)                 = 1 in case A  c     casea (imt)                 = 1 in case A
1013  c     diffus(imt,0:Nrp1,mdiff)    vertical diffusivities              (m^2/s)  c     diffus(imt,0:Nrp1,mdiff)    vertical diffusivities              (m^2/s)
1014  c     kbl(imt)                    -1 of first grid level below hbl  c     kbl   (imt)                 -1 of first grid level below hbl
1015        _KPP_RL ustar (imt)        _RL ustar (imt)
1016        _KPP_RL bfsfc (imt)        _RL bfsfc (imt)
1017        _KPP_RL hbl   (imt)        _RL hbl   (imt)
1018        _KPP_RL stable(imt)        _RL stable(imt)
1019        _KPP_RL casea (imt)        _RL casea (imt)
1020        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL diffus(imt,0:Nrp1,mdiff)
1021        integer kbl(imt)        integer kbl(imt)
1022    
1023  c output  c output
# Line 1095  c     dkm1 (imt,mdiff)            bounda Line 1025  c     dkm1 (imt,mdiff)            bounda
1025  c     blmc (imt,Nr,mdiff)         boundary layer mixing coefficients  (m^2/s)  c     blmc (imt,Nr,mdiff)         boundary layer mixing coefficients  (m^2/s)
1026  c     ghat (imt,Nr)               nonlocal scalar transport  c     ghat (imt,Nr)               nonlocal scalar transport
1027  c     sigma(imt)                  normalized depth (d / hbl)  c     sigma(imt)                  normalized depth (d / hbl)
1028        _KPP_RL dkm1 (imt,mdiff)        _RL dkm1 (imt,mdiff)
1029        _KPP_RL blmc (imt,Nr,mdiff)        _RL blmc (imt,Nr,mdiff)
1030        _KPP_RL ghat (imt,Nr)        _RL ghat (imt,Nr)
1031        _KPP_RL sigma(imt)        _RL sigma(imt)
1032        integer ikppkey, kkppkey        integer ikppkey, kkppkey
1033    
1034  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
# Line 1107  c  local Line 1037  c  local
1037  c     gat1*(imt)                 shape function at sigma = 1  c     gat1*(imt)                 shape function at sigma = 1
1038  c     dat1*(imt)                 derivative of shape function at sigma = 1  c     dat1*(imt)                 derivative of shape function at sigma = 1
1039  c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)  c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)
1040        _KPP_RL gat1m(imt), gat1s(imt), gat1t(imt)        _RL gat1m(imt), gat1s(imt), gat1t(imt)
1041        _KPP_RL dat1m(imt), dat1s(imt), dat1t(imt)        _RL dat1m(imt), dat1s(imt), dat1t(imt)
1042        _KPP_RL ws(imt), wm(imt)        _RL ws(imt), wm(imt)
1043        integer i, kn, ki        integer i, kn, ki
1044        _KPP_RL R, dvdzup, dvdzdn, viscp        _RL R, dvdzup, dvdzdn, viscp
1045        _KPP_RL difsp, diftp, visch, difsh, difth        _RL difsp, diftp, visch, difsh, difth
1046        _KPP_RL f1, sig, a1, a2, a3, delhat        _RL f1, sig, a1, a2, a3, delhat
1047        _KPP_RL Gm, Gs, Gt        _RL Gm, Gs, Gt
1048        _KPP_RL tempVar        _RL tempVar
1049    
1050        _KPP_RL    p0    , eins        _RL    p0    , eins
1051        parameter (p0=0.0, eins=1.0)        parameter (p0=0.0, eins=1.0)
1052    
1053  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 1320  c     hbl(imt)                  boundary Line 1250  c     hbl(imt)                  boundary
1250  c     kbl(imt)                  grid above hbl  c     kbl(imt)                  grid above hbl
1251  c     diffus(imt,0:Nrp1,mdiff) vertical diffusivities           (m^2/s)  c     diffus(imt,0:Nrp1,mdiff) vertical diffusivities           (m^2/s)
1252  c     casea(imt)                = 1 in caseA, = 0 in case B  c     casea(imt)                = 1 in caseA, = 0 in case B
1253        _KPP_RL dkm1  (imt,mdiff)        _RL dkm1  (imt,mdiff)
1254        _KPP_RL hbl   (imt)        _RL hbl   (imt)
1255        integer kbl   (imt)        integer kbl   (imt)
1256        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL diffus(imt,0:Nrp1,mdiff)
1257        _KPP_RL casea (imt)        _RL casea (imt)
1258    
1259  c input/output  c input/output
1260  c     nonlocal transport, modified ghat at kbl(i)-1 interface    (s/m**2)  c     nonlocal transport, modified ghat at kbl(i)-1 interface    (s/m**2)
1261        _KPP_RL ghat (imt,Nr)        _RL ghat (imt,Nr)
1262    
1263  c output  c output
1264  c     enhanced bound. layer mixing coeff.  c     enhanced bound. layer mixing coeff.
1265        _KPP_RL blmc  (imt,Nr,mdiff)        _RL blmc  (imt,Nr,mdiff)
1266    
1267  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1268    
1269  c local  c local
1270  c     fraction hbl lies beteen zgrid neighbors  c     fraction hbl lies beteen zgrid neighbors
1271        _KPP_RL delta        _RL delta
1272        integer ki, i, md        integer ki, i, md
1273        _KPP_RL dkmp5, dstar        _RL dkmp5, dstar
1274    
1275        do i = 1, imt        do i = 1, imt
1276           ki = kbl(i)-1           ki = kbl(i)-1
# Line 1392  c Line 1322  c
1322  c     written  by: jan morzel,   feb. 10, 1995 (converted from "sigma" version)  c     written  by: jan morzel,   feb. 10, 1995 (converted from "sigma" version)
1323  c     modified by: d. menemenlis,    june 1998 : for use with MIT GCM UV  c     modified by: d. menemenlis,    june 1998 : for use with MIT GCM UV
1324  c  c
1325    c     28 april 05: added computation of mixed layer depth KPPmld
1326    c     for a deltaRho equivalent to a deltaTheta of 0.8 deg C
1327    
1328  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1329    
1330        IMPLICIT NONE        IMPLICIT NONE
# Line 1401  c--------------------------------------- Line 1334  c---------------------------------------
1334  #include "PARAMS.h"  #include "PARAMS.h"
1335  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
1336  #include "DYNVARS.h"  #include "DYNVARS.h"
1337    #include "GRID.h"
1338    
1339  c-------------- Routine arguments -----------------------------------------  c-------------- Routine arguments -----------------------------------------
1340        INTEGER bi, bj, myThid        INTEGER bi, bj, myThid
1341        _KPP_RL RHO1   ( ibot:itop, jbot:jtop       )        _RL RHO1   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy       )
1342        _KPP_RL DBLOC  ( ibot:itop, jbot:jtop, Nr   )        _RL DBLOC  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
1343        _KPP_RL DBSFC  ( ibot:itop, jbot:jtop, Nr   )        _RL DBSFC  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
1344        _KPP_RL TTALPHA( ibot:itop, jbot:jtop, Nrp1 )        _RL TTALPHA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1345        _KPP_RL SSBETA ( ibot:itop, jbot:jtop, Nrp1 )        _RL SSBETA ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1346    
1347  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1348    
# Line 1419  c Line 1353  c
1353  c     rhok         - density of t(k  ) & s(k  ) at depth k  c     rhok         - density of t(k  ) & s(k  ) at depth k
1354  c     rhokm1       - density of t(k-1) & s(k-1) at depth k  c     rhokm1       - density of t(k-1) & s(k-1) at depth k
1355  c     rho1k        - density of t(1  ) & s(1  ) at depth k  c     rho1k        - density of t(1  ) & s(1  ) at depth k
1356  c     work1, work2 - work arrays for holding horizontal slabs  c     work1,2,3    - work arrays for holding horizontal slabs
1357    
1358        _RL RHOK  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL RHOK  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1359        _RL RHOKM1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL RHOKM1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 1427  c     work1, work2 - work arrays for hol Line 1361  c     work1, work2 - work arrays for hol
1361        _RL WORK1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1362        _RL WORK2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1363        _RL WORK3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1364    #ifdef ALLOW_DIAGNOSTICS
1365    c     KPPMLD - mixed layer depth based on density criterion
1366          _RL KPPMLD(1    :sNx    ,1    :sNy    )
1367    #endif /* ALLOW_DIAGNOSTICS */
1368    
1369        INTEGER I, J, K        INTEGER I, J, K
1370        INTEGER ikppkey, kkppkey        INTEGER ikppkey, kkppkey
1371    
# Line 1439  CADJ STORE theta(:,:,1,bi,bj) = comlev1_ Line 1378  CADJ STORE theta(:,:,1,bi,bj) = comlev1_
1378  CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey  CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1379  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1380        call FIND_RHO(        call FIND_RHO(
1381       I     bi, bj, ibot, itop, jbot, jtop, 1, 1,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1382       I     theta, salt,       I     theta, salt,
1383       O     WORK1,       O     WORK1,
1384       I     myThid )       I     myThid )
# Line 1449  CADJ STORE salt (:,:,1,bi,bj) = comlev1_ Line 1388  CADJ STORE salt (:,:,1,bi,bj) = comlev1_
1388  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1389    
1390        call FIND_ALPHA(        call FIND_ALPHA(
1391       I     bi, bj, ibot, itop, jbot, jtop, 1, 1,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1392       O     WORK2 )       O     WORK2 )
1393    
1394        call FIND_BETA(        call FIND_BETA(
1395       I     bi, bj, ibot, itop, jbot, jtop, 1, 1,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1396       O     WORK3 )       O     WORK3 )
1397    
1398        DO J = jbot, jtop        DO J = 1-OLy, sNy+OLy
1399           DO I = ibot, itop           DO I = 1-OLx, sNx+OLx
1400              RHO1(I,J)      = WORK1(I,J) + rhoConst              RHO1(I,J)      = WORK1(I,J) + rhoConst
1401              TTALPHA(I,J,1) = WORK2(I,J)              TTALPHA(I,J,1) = WORK2(I,J)
1402              SSBETA(I,J,1)  = WORK3(I,J)              SSBETA(I,J,1)  = WORK3(I,J)
# Line 1465  CADJ STORE salt (:,:,1,bi,bj) = comlev1_ Line 1404  CADJ STORE salt (:,:,1,bi,bj) = comlev1_
1404           END DO           END DO
1405        END DO        END DO
1406    
1407    #ifdef ALLOW_DIAGNOSTICS
1408    c     work3 - density of t(1)-.8 & s(1  ) at depth 1
1409          IF ( useDiagnostics ) THEN
1410             DO J = 1, sNy
1411                DO I = 1, sNx
1412                   KPPMLD(I,J) = ABS(R_low(I,J,bi,bj))
1413                   WORK3 (I,J) = WORK1(I,J) - 0.8 * WORK2(I,J)
1414                END DO
1415             END DO
1416          ENDIF
1417    #endif /* ALLOW_DIAGNOSTICS */
1418    
1419  c calculate alpha, beta, and gradients in interior layers  c calculate alpha, beta, and gradients in interior layers
1420    
1421  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)
# Line 1477  CADJ STORE theta(:,:,k,bi,bj) = comlev1_ Line 1428  CADJ STORE theta(:,:,k,bi,bj) = comlev1_
1428  CADJ STORE salt (:,:,k,bi,bj) = comlev1_kpp_k, key=kkppkey  CADJ STORE salt (:,:,k,bi,bj) = comlev1_kpp_k, key=kkppkey
1429  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1430           call FIND_RHO(           call FIND_RHO(
1431       I        bi, bj, ibot, itop, jbot, jtop, K, K,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1432       I        theta, salt,       I        theta, salt,
1433       O        RHOK,       O        RHOK,
1434       I        myThid )       I        myThid )
# Line 1487  CADJ STORE theta(:,:,k-1,bi,bj) = comlev Line 1438  CADJ STORE theta(:,:,k-1,bi,bj) = comlev
1438  CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_kpp_k, key=kkppkey  CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_kpp_k, key=kkppkey
1439  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1440           call FIND_RHO(           call FIND_RHO(
1441       I        bi, bj, ibot, itop, jbot, jtop, K-1, K,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K-1, K,
1442       I        theta, salt,       I        theta, salt,
1443       O        RHOKM1,       O        RHOKM1,
1444       I        myThid )       I        myThid )
# Line 1497  CADJ STORE theta(:,:,1,bi,bj) = comlev1_ Line 1448  CADJ STORE theta(:,:,1,bi,bj) = comlev1_
1448  CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey  CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1449  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1450           call FIND_RHO(           call FIND_RHO(
1451       I        bi, bj, ibot, itop, jbot, jtop, 1, K,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, K,
1452       I        theta, salt,       I        theta, salt,
1453       O        RHO1K,       O        RHO1K,
1454       I        myThid )       I        myThid )
# Line 1509  CADJ STORE rho1k (:,:)          = comlev Line 1460  CADJ STORE rho1k (:,:)          = comlev
1460  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1461    
1462           call FIND_ALPHA(           call FIND_ALPHA(
1463       I        bi, bj, ibot, itop, jbot, jtop, K, K,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1464       O        WORK1 )       O        WORK1 )
1465    
1466           call FIND_BETA(           call FIND_BETA(
1467       I        bi, bj, ibot, itop, jbot, jtop, K, K,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1468       O        WORK2 )       O        WORK2 )
1469    
1470           DO J = jbot, jtop           DO J = 1-OLy, sNy+OLy
1471              DO I = ibot, itop              DO I = 1-OLx, sNx+OLx
1472                 TTALPHA(I,J,K) = WORK1 (I,J)                 TTALPHA(I,J,K) = WORK1 (I,J)
1473                 SSBETA(I,J,K)  = WORK2 (I,J)                 SSBETA(I,J,K)  = WORK2 (I,J)
1474                 DBLOC(I,J,K-1) = gravity * (RHOK(I,J) - RHOKM1(I,J)) /                 DBLOC(I,J,K-1) = gravity * (RHOK(I,J) - RHOKM1(I,J)) /
# Line 1527  CADJ STORE rho1k (:,:)          = comlev Line 1478  CADJ STORE rho1k (:,:)          = comlev
1478              END DO              END DO
1479           END DO           END DO
1480    
1481    #ifdef ALLOW_DIAGNOSTICS
1482             IF ( useDiagnostics ) THEN
1483    c     work1 - density of t(k-1)  & s(k-1) at depth 1
1484    c     work2 - density of t(k  )  & s(k  ) at depth 1
1485    c     work3 - density of t(1)-.8 & s(1  ) at depth 1
1486                call FIND_RHO(
1487         I           bi, bj, 1, sNx, 1, sNy, K-1, 1, theta, salt,
1488         O           WORK1,
1489         I           myThid )
1490                call FIND_RHO(
1491         I           bi, bj, 1, sNx, 1, sNy, K  , 1, theta, salt,
1492         O           WORK2,
1493         I           myThid )
1494                DO J = 1, sNy
1495                   DO I = 1, sNx
1496                      IF ( k .LE. klowC(I,J,bi,bj) .AND.
1497         &                 WORK1(I,J) .LT. WORK3(I,J) .AND.
1498         &                 WORK2(I,J) .GE. WORK3(I,J) ) THEN
1499                         KPPMLD(I,J) = MIN ( KPPMLD(I,J),
1500         &                    ABS(((WORK3(I,J)-WORK1(I,J))*rC(k)+
1501         &                    (WORK2(I,J)-WORK3(I,J))*rC(k-1))/
1502         &                    (WORK2(I,J)-WORK1(I,J))))
1503                      ENDIF
1504                   END DO
1505                END DO
1506             ENDIF
1507    #endif /* ALLOW_DIAGNOSTICS */
1508    
1509        END DO        END DO
1510    
1511  c     compute arrays for K = Nrp1  c     compute arrays for K = Nrp1
1512        DO J = jbot, jtop        DO J = 1-OLy, sNy+OLy
1513           DO I = ibot, itop           DO I = 1-OLx, sNx+OLx
1514              TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)              TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)
1515              SSBETA(I,J,Nrp1)  = SSBETA(I,J,Nr)              SSBETA(I,J,Nrp1)  = SSBETA(I,J,Nr)
1516              DBLOC(I,J,Nr)     = 0.              DBLOC(I,J,Nr)     = 0.
1517           END DO           END DO
1518        END DO        END DO
1519    
1520    #ifdef ALLOW_DIAGNOSTICS
1521          IF ( useDiagnostics ) THEN
1522             CALL DIAGNOSTICS_FILL(KPPmld,'KPPmld  ',0,1,3,bi,bj,myThid)
1523          ENDIF
1524    #endif /* ALLOW_DIAGNOSTICS */
1525    
1526  #endif /* ALLOW_KPP */  #endif /* ALLOW_KPP */
1527    
1528        RETURN        RETURN
Legend:
Removed from v.1.20  
changed lines
  Added in v.1.29
  ViewVC Help
Powered by ViewVC 1.1.22