/[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.9 by heimbach, Tue Sep 18 20:30:59 2001 UTC revision 1.30 by mlosch, Tue May 1 04:09:25 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 20  C--  o STATEKPP    - Compute buoyancy-re Line 19  C--  o STATEKPP    - Compute buoyancy-re
19  c***********************************************************************  c***********************************************************************
20    
21        SUBROUTINE KPPMIX (        SUBROUTINE KPPMIX (
22       I       mytime, mythid       I       kmtj, shsq, dvsq, ustar, msk
      I     , kmtj, shsq, dvsq, ustar  
23       I     , bo, bosol, dbloc, Ritop, coriol       I     , bo, bosol, dbloc, Ritop, coriol
24       I     , ikey       I     , diffusKzS, diffusKzT
25         I     , ikppkey
26       O     , diffus       O     , diffus
27       U     , ghat       U     , ghat
28       O     , hbl )       O     , hbl
29         I     , mytime, mythid )
30    
31  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
32  c  c
# 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     msk    (imt)     - surface mask (=1 if water, =0 otherwise)
57  c     dvsq   (imt,Nr) - (velocity shear re sfc)^2                     ((m/s)^2)  c     shsq   (imt,Nr)  - (local velocity shear)^2                     ((m/s)^2)
58  c     ustar  (imt)    - surface friction velocity                         (m/s)  c     dvsq   (imt,Nr)  - (velocity shear re sfc)^2                    ((m/s)^2)
59  c     bo     (imt)    - surface turbulent buoy. forcing               (m^2/s^3)  c     ustar  (imt)     - surface friction velocity                        (m/s)
60  c     bosol  (imt)    - radiative buoyancy forcing                    (m^2/s^3)  c     bo     (imt)     - surface turbulent buoy. forcing              (m^2/s^3)
61  c     dbloc  (imt,Nr) - local delta buoyancy across interfaces          (m/s^2)  c     bosol  (imt)     - radiative buoyancy forcing                   (m^2/s^3)
62  c     dblocSm(imt,Nr) - horizontally smoothed dbloc                     (m/s^2)  c     dbloc  (imt,Nr)  - local delta buoyancy across interfaces         (m/s^2)
63  c                         stored in ghat to save space  c     dblocSm(imt,Nr)  - horizontally smoothed dbloc                    (m/s^2)
64  c     Ritop  (imt,Nr) - numerator of bulk Richardson Number  c                          stored in ghat to save space
65  c                         (zref-z) * delta buoyancy w.r.t. surface    ((m/s)^2)  c     Ritop  (imt,Nr)  - numerator of bulk Richardson Number
66  c     coriol (imt)    - Coriolis parameter                                (1/s)  c                          (zref-z) * delta buoyancy w.r.t. surface   ((m/s)^2)
67    c     coriol (imt)     - Coriolis parameter                               (1/s)
68    c     diffusKzS(imt,Nr)- background vertical diffusivity for scalars    (m^2/s)
69    c     diffusKzT(imt,Nr)- background vertical diffusivity for theta      (m^2/s)
70  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,
71  c           e.g., hbl(sNx,sNy) -> hbl(imt),  c           e.g., hbl(sNx,sNy) -> hbl(imt),
72  c           where hbl(i,j) -> hbl((j-1)*sNx+i)  c           where hbl(i,j) -> hbl((j-1)*sNx+i)
73    
74        _RL     mytime        _RL     mytime
75        integer mythid        integer mythid
76        integer kmtj  (imt   )        integer kmtj (imt   )
77        _KPP_RL shsq  (imt,Nr)        _RL shsq     (imt,Nr)
78        _KPP_RL dvsq  (imt,Nr)        _RL dvsq     (imt,Nr)
79        _KPP_RL ustar (imt   )        _RL ustar    (imt   )
80        _KPP_RL bo    (imt   )        _RL bo       (imt   )
81        _KPP_RL bosol (imt   )        _RL bosol    (imt   )
82        _KPP_RL dbloc (imt,Nr)        _RL dbloc    (imt,Nr)
83        _KPP_RL Ritop (imt,Nr)        _RL Ritop    (imt,Nr)
84        _KPP_RL coriol(imt   )        _RL coriol   (imt   )
85          _RS msk      (imt   )
86          _RL diffusKzS(imt,Nr)
87          _RL diffusKzT(imt,Nr)
88    
89        integer ikey        integer ikppkey
90    
91  c output  c output
92  c     diffus (imt,1)  - vertical viscosity coefficient                  (m^2/s)  c     diffus (imt,1)  - vertical viscosity coefficient                  (m^2/s)
# Line 91  c     diffus (imt,3)  - vertical tempera Line 95  c     diffus (imt,3)  - vertical tempera
95  c     ghat   (imt)    - nonlocal transport coefficient                  (s/m^2)  c     ghat   (imt)    - nonlocal transport coefficient                  (s/m^2)
96  c     hbl    (imt)    - mixing layer depth                                  (m)  c     hbl    (imt)    - mixing layer depth                                  (m)
97    
98        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL diffus(imt,0:Nrp1,mdiff)
99        _KPP_RL ghat  (imt,Nr)        _RL ghat  (imt,Nr)
100        _KPP_RL hbl   (imt)        _RL hbl   (imt)
101    
102  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
103    
# Line 107  c     blmc   (imt,Nr,mdiff) - boundary l Line 111  c     blmc   (imt,Nr,mdiff) - boundary l
111  c     sigma  (imt         ) - normalized depth (d / hbl)  c     sigma  (imt         ) - normalized depth (d / hbl)
112  c     Rib    (imt,Nr      ) - bulk Richardson number  c     Rib    (imt,Nr      ) - bulk Richardson number
113    
114        integer kbl   (imt         )        integer kbl(imt         )
115        _KPP_RL bfsfc (imt         )        _RL bfsfc  (imt         )
116        _KPP_RL casea (imt         )        _RL casea  (imt         )
117        _KPP_RL stable(imt         )        _RL stable (imt         )
118        _KPP_RL dkm1  (imt,   mdiff)        _RL dkm1   (imt,   mdiff)
119        _KPP_RL blmc  (imt,Nr,mdiff)        _RL blmc   (imt,Nr,mdiff)
120        _KPP_RL sigma (imt         )        _RL sigma  (imt         )
121        _KPP_RL Rib   (imt,Nr      )        _RL Rib    (imt,Nr      )
122    
123        integer i, k, md        integer i, k, md
124    
# Line 125  c instability. Line 129  c instability.
129  c (ghat is temporary storage for horizontally smoothed dbloc)  c (ghat is temporary storage for horizontally smoothed dbloc)
130  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
131    
132  CADJ STORE ghat = comlev1_kpp, key = ikey  cph(
133    cph these storings avoid recomp. of Ri_iwmix
134    CADJ STORE ghat  = comlev1_kpp, key = ikppkey
135    CADJ STORE dbloc = comlev1_kpp, key = ikppkey
136    cph)
137        call Ri_iwmix (        call Ri_iwmix (
138       I       kmtj, shsq, dbloc, ghat       I       kmtj, shsq, dbloc, ghat
139       I     , ikey       I     , diffusKzS, diffusKzT
140         I     , ikppkey
141       O     , diffus )       O     , diffus )
142    
143    cph(
144    cph these storings avoid recomp. of Ri_iwmix
145    cph DESPITE TAFs 'not necessary' warning!
146    CADJ STORE dbloc  = comlev1_kpp, key = ikppkey
147    CADJ STORE shsq   = comlev1_kpp, key = ikppkey
148    CADJ STORE ghat   = comlev1_kpp, key = ikppkey
149    CADJ STORE diffus = comlev1_kpp, key = ikppkey
150    cph)
151    
152  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
153  c set seafloor values to zero and fill extra "Nrp1" coefficients  c set seafloor values to zero and fill extra "Nrp1" coefficients
154  c for blmix  c for blmix
# Line 152  c diagnose the new boundary layer depth Line 169  c diagnose the new boundary layer depth
169  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
170    
171        call bldepth (        call bldepth (
172       I       mytime, mythid       I       kmtj
      I     , kmtj  
173       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol
174       I     , ikey       I     , ikppkey
175       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma
176       &     )       I     , mytime, mythid )
177    
178  CADJ STORE hbl,bfsfc,stable,casea,kbl = comlev1_kpp, key = ikey  CADJ STORE hbl,bfsfc,stable,casea,kbl = comlev1_kpp, key = ikppkey
179    
180  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
181  c compute boundary layer diffusivities  c compute boundary layer diffusivities
# Line 167  c--------------------------------------- Line 183  c---------------------------------------
183    
184        call blmix (        call blmix (
185       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl
186       O     , dkm1, blmc, ghat, sigma, ikey       O     , dkm1, blmc, ghat, sigma, ikppkey
187       &     )       I     , mythid )
188    cph(
189  CADJ STORE dkm1,blmc,ghat = comlev1_kpp, key = ikey  CADJ STORE dkm1,blmc,ghat = comlev1_kpp, key = ikppkey
190    CADJ STORE hbl, kbl, diffus, casea = comlev1_kpp, key = ikppkey
191    cph)
192    
193  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
194  c enhance diffusivity at interface kbl - 1  c enhance diffusivity at interface kbl - 1
# Line 179  c--------------------------------------- Line 197  c---------------------------------------
197        call enhance (        call enhance (
198       I       dkm1, hbl, kbl, diffus, casea       I       dkm1, hbl, kbl, diffus, casea
199       U     , ghat       U     , ghat
200       O     , blmc )       O     , blmc
201         I     , mythid )
202    
203    cph(
204    cph avoids recomp. of enhance
205    CADJ STORE blmc = comlev1_kpp, key = ikppkey
206    cph)
207    
208  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
209  c combine interior and boundary layer coefficients and nonlocal term  c combine interior and boundary layer coefficients and nonlocal term
210    c !!!NOTE!!! In shallow (2-level) regions and for shallow mixed layers
211    c (< 1 level), diffusivity blmc can become negative.  The max-s below
212    c are a hack until this problem is properly diagnosed and fixed.
213  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
   
214        do k = 1, Nr        do k = 1, Nr
215           do i = 1, imt           do i = 1, imt
216              if (k .lt. kbl(i)) then              if (k .lt. kbl(i)) then
217                 do md = 1, mdiff  #ifdef ALLOW_SHELFICE
218                    diffus(i,k,md) = blmc(i,k,md)  C     when there is shelfice on top (msk(i)=0), reset the boundary layer
219                 end do  C     mixing coefficients blmc to pure Ri-number based mixing
220                   blmc(i,k,1) = max ( blmc(i,k,1)*msk(i),
221         &              diffus(i,k,1) )
222                   blmc(i,k,2) = max ( blmc(i,k,2)*msk(i),
223         &              diffus(i,k,2) )
224                   blmc(i,k,3) = max ( blmc(i,k,3)*msk(i),
225         &              diffus(i,k,3) )
226    #endif /* not ALLOW_SHELFICE */
227                   diffus(i,k,1) = max ( blmc(i,k,1), viscAr  )
228                   diffus(i,k,2) = max ( blmc(i,k,2), diffusKzS(i,Nr) )
229                   diffus(i,k,3) = max ( blmc(i,k,3), diffusKzT(i,Nr) )
230              else              else
231                 ghat(i,k) = 0.                 ghat(i,k) = 0.
232              endif              endif
# Line 205  c--------------------------------------- Line 241  c---------------------------------------
241  c*************************************************************************  c*************************************************************************
242    
243        subroutine bldepth (        subroutine bldepth (
244       I       mytime, mythid       I       kmtj
      I     , kmtj  
245       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol
246       I     , ikey       I     , ikppkey
247       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma
248       &     )       I     , mytime, mythid )
249    
250  c     the oceanic planetary boundary layer depth, hbl, is determined as  c     the oceanic planetary boundary layer depth, hbl, is determined as
251  c     the shallowest depth where the bulk Richardson number is  c     the shallowest depth where the bulk Richardson number is
# Line 236  c Line 271  c
271        IMPLICIT NONE        IMPLICIT NONE
272    
273  #include "SIZE.h"  #include "SIZE.h"
 #include "EEPARAMS.h"  
 #include "PARAMS.h"  
274  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
 #include "FFIELDS.h"  
275    
276  c input  c input
277  c------  c------
# Line 258  c coriol    : Coriolis parameter Line 290  c coriol    : Coriolis parameter
290        _RL     mytime        _RL     mytime
291        integer mythid        integer mythid
292        integer kmtj(imt)        integer kmtj(imt)
293        _KPP_RL dvsq  (imt,Nr)        _RL dvsq    (imt,Nr)
294        _KPP_RL dbloc (imt,Nr)        _RL dbloc   (imt,Nr)
295        _KPP_RL Ritop (imt,Nr)        _RL Ritop   (imt,Nr)
296        _KPP_RL ustar (imt)        _RL ustar   (imt)
297        _KPP_RL bo    (imt)        _RL bo      (imt)
298        _KPP_RL bosol (imt)        _RL bosol   (imt)
299        _KPP_RL coriol(imt)        _RL coriol  (imt)
300        integer ikey        integer ikppkey, kkppkey
301    
302  c  output  c  output
303  c--------  c--------
# Line 276  c casea     : =1 in case A, =0 in case B Line 308  c casea     : =1 in case A, =0 in case B
308  c kbl       : -1 of first grid level below hbl  c kbl       : -1 of first grid level below hbl
309  c Rib       : Bulk Richardson number  c Rib       : Bulk Richardson number
310  c sigma     : normalized depth (d/hbl)  c sigma     : normalized depth (d/hbl)
311        _KPP_RL hbl   (imt)        _RL hbl    (imt)
312        _KPP_RL bfsfc (imt)        _RL bfsfc  (imt)
313        _KPP_RL stable(imt)        _RL stable (imt)
314        _KPP_RL casea (imt)        _RL casea  (imt)
315        integer kbl   (imt)        integer kbl(imt)
316        _KPP_RL Rib   (imt,Nr)        _RL Rib    (imt,Nr)
317        _KPP_RL sigma (imt)        _RL sigma  (imt)
318    
319  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
320    
321  c  local  c  local
322  c-------  c-------
323  c wm, ws    : turbulent velocity scales         (m/s)  c wm, ws    : turbulent velocity scales         (m/s)
324        _KPP_RL wm(imt), ws(imt)        _RL wm(imt), ws(imt)
325        _RL     worka(imt)        _RL worka(imt)
326    
327        _KPP_RL bvsq, vtsq, hekman, hmonob, hlimit, tempVar1, tempVar2        _RL bvsq, vtsq, hekman, hmonob, hlimit, tempVar1, tempVar2
328        integer i, kl        integer i, kl
329    
330        _KPP_RL     p5    , eins        _RL         p5    , eins
331        parameter ( p5=0.5, eins=1.0 )        parameter ( p5=0.5, eins=1.0 )
332        _RL         minusone        _RL         minusone
333        parameter ( minusone=-1.0 )        parameter ( minusone=-1.0 )
# Line 319  c     initialize hbl and kbl to bottomed Line 351  c     initialize hbl and kbl to bottomed
351    
352        do kl = 2, Nr        do kl = 2, Nr
353    
354    #ifdef ALLOW_AUTODIFF_TAMC
355             kkppkey = (ikppkey-1)*Nr + kl
356    #endif
357    
358  c     compute bfsfc = sw fraction at hbf * zgrid  c     compute bfsfc = sw fraction at hbf * zgrid
359    
360           do i = 1, imt           do i = 1, imt
361              worka(i) = zgrid(kl)              worka(i) = zgrid(kl)
362           end do           end do
363    CADJ store worka = comlev1_kpp_k, key = kkppkey
364           call SWFRAC(           call SWFRAC(
365       I        imt, hbf,       I        imt, hbf,
366       I        mytime, mythid,       I        mytime, mythid,
367       U        worka )       U        worka )
368    CADJ store worka = comlev1_kpp_k, key = kkppkey
369    
370    
371           do i = 1, imt           do i = 1, imt
372    
# Line 349  c--------------------------------------- Line 388  c---------------------------------------
388    
389           call wscale (           call wscale (
390       I        sigma, casea, ustar, bfsfc,       I        sigma, casea, ustar, bfsfc,
391       O        wm, ws )       O        wm, ws, myThid )
392    CADJ store ws = comlev1_kpp_k, key = kkppkey
393    
394           do i = 1, imt           do i = 1, imt
395    
# Line 385  c Line 425  c
425    
426           end do           end do
427        end do        end do
428          
429    cph(
430    cph  without this store, there is a recomputation error for
431    cph  rib in adbldepth (probably partial recomputation problem)    
432    CADJ store Rib = comlev1_kpp
433    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy),Nr /)
434    cph)
435    
436        do kl = 2, Nr        do kl = 2, Nr
437           do i = 1, imt           do i = 1, imt
438              if (kbl(i).eq.kmtj(i) .and. Rib(i,kl).gt.Ricr) kbl(i) = kl              if (kbl(i).eq.kmtj(i) .and. Rib(i,kl).gt.Ricr) kbl(i) = kl
# Line 393  c Line 440  c
440        end do        end do
441    
442  CADJ store kbl = comlev1_kpp  CADJ store kbl = comlev1_kpp
443  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
444    
445        do i = 1, imt        do i = 1, imt
446           kl = kbl(i)           kl = kbl(i)
# Line 406  c     linearly interpolate to find hbl w Line 453  c     linearly interpolate to find hbl w
453        end do        end do
454    
455  CADJ store hbl = comlev1_kpp  CADJ store hbl = comlev1_kpp
456  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
457    
458  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
459  c     find stability and buoyancy forcing for boundary layer  c     find stability and buoyancy forcing for boundary layer
# Line 415  c--------------------------------------- Line 462  c---------------------------------------
462        do i = 1, imt        do i = 1, imt
463           worka(i) = hbl(i)           worka(i) = hbl(i)
464        end do        end do
465    CADJ store worka = comlev1_kpp
466    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
467        call SWFRAC(        call SWFRAC(
468       I     imt, minusone,       I     imt, minusone,
469       I     mytime, mythid,       I     mytime, mythid,
470       U     worka )       U     worka )
471    CADJ store worka = comlev1_kpp
472    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
473    
474        do i = 1, imt        do i = 1, imt
475           bfsfc(i)  = bo(i) + bosol(i) * (1. - worka(i))           bfsfc(i)  = bo(i) + bosol(i) * (1. - worka(i))
476        end do        end do
477  CADJ store bfsfc = comlev1_kpp  CADJ store bfsfc = comlev1_kpp
478  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
479    
480  c--   ensure bfsfc is never 0  c--   ensure bfsfc is never 0
481        do i = 1, imt        do i = 1, imt
# Line 432  c--   ensure bfsfc is never 0 Line 483  c--   ensure bfsfc is never 0
483           bfsfc(i) = sign(eins,bfsfc(i))*max(phepsi,abs(bfsfc(i)))           bfsfc(i) = sign(eins,bfsfc(i))*max(phepsi,abs(bfsfc(i)))
484        end do        end do
485    
486  CADJ store bfsfc = comlev1_kpp  cph(
487  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  cph  added stable to store list to avoid extensive recomp.
488    CADJ store bfsfc, stable = comlev1_kpp
489    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
490    cph)
491    
492  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
493  c check hbl limits for hekman or hmonob  c check hbl limits for hekman or hmonob
# Line 451  c--------------------------------------- Line 505  c---------------------------------------
505           end if           end if
506        end do        end do
507  CADJ store hbl = comlev1_kpp  CADJ store hbl = comlev1_kpp
508  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
509    
510        do i = 1, imt        do i = 1, imt
511           hbl(i) = max(hbl(i),minKPPhbl)           hbl(i) = max(hbl(i),minKPPhbl)
# Line 459  CADJ &   , key = ikey, shape = (/ (sNx+2 Line 513  CADJ &   , key = ikey, shape = (/ (sNx+2
513        end do        end do
514    
515  CADJ store hbl = comlev1_kpp  CADJ store hbl = comlev1_kpp
516  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
517    
518  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
519  c      find new kbl  c      find new kbl
# Line 480  c--------------------------------------- Line 534  c---------------------------------------
534        do i = 1, imt        do i = 1, imt
535           worka(i) = hbl(i)           worka(i) = hbl(i)
536        end do        end do
537    CADJ store worka = comlev1_kpp
538    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
539        call SWFRAC(        call SWFRAC(
540       I     imt, minusone,       I     imt, minusone,
541       I     mytime, mythid,       I     mytime, mythid,
542       U     worka )       U     worka )
543    CADJ store worka = comlev1_kpp
544    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
545    
546        do i = 1, imt        do i = 1, imt
547           bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i))           bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i))
548        end do        end do
549  CADJ store bfsfc = comlev1_kpp  CADJ store bfsfc = comlev1_kpp
550  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
551    
552  c--   ensures bfsfc is never 0  c--   ensures bfsfc is never 0
553        do i = 1, imt        do i = 1, imt
# Line 515  c*************************************** Line 573  c***************************************
573    
574        subroutine wscale (        subroutine wscale (
575       I     sigma, hbl, ustar, bfsfc,       I     sigma, hbl, ustar, bfsfc,
576       O     wm, ws )       O     wm, ws,
577         I     myThid )
578    
579  c     compute turbulent velocity scales.  c     compute turbulent velocity scales.
580  c     use a 2D-lookup table for wm and ws as functions of ustar and  c     use a 2D-lookup table for wm and ws as functions of ustar and
# Line 536  c sigma   : normalized depth (d/hbl) Line 595  c sigma   : normalized depth (d/hbl)
595  c hbl     : boundary layer depth (m)  c hbl     : boundary layer depth (m)
596  c ustar   : surface friction velocity         (m/s)  c ustar   : surface friction velocity         (m/s)
597  c bfsfc   : total surface buoyancy flux       (m^2/s^3)  c bfsfc   : total surface buoyancy flux       (m^2/s^3)
598        _KPP_RL sigma(imt)  c myThid  : thread number for this instance of the routine
599        _KPP_RL hbl  (imt)        integer myThid
600        _KPP_RL ustar(imt)        _RL sigma(imt)
601        _KPP_RL bfsfc(imt)        _RL hbl  (imt)
602          _RL ustar(imt)
603          _RL bfsfc(imt)
604    
605  c  output  c  output
606  c--------  c--------
607  c wm, ws  : turbulent velocity scales at sigma  c wm, ws  : turbulent velocity scales at sigma
608        _KPP_RL wm(imt), ws(imt)        _RL wm(imt), ws(imt)
609    
610  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
611    
612  c local  c local
613  c------  c------
614  c zehat   : = zeta *  ustar**3  c zehat   : = zeta *  ustar**3
615        _KPP_RL zehat        _RL zehat
616    
617        integer iz, izp1, ju, i, jup1        integer iz, izp1, ju, i, jup1
618        _KPP_RL udiff, zdiff, zfrac, ufrac, fzfrac, wam        _RL udiff, zdiff, zfrac, ufrac, fzfrac, wam
619        _KPP_RL wbm, was, wbs, u3, tempVar        _RL wbm, was, wbs, u3, tempVar
620    
621  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
622  c use lookup table for zehat < zmax only; otherwise use  c use lookup table for zehat < zmax only; otherwise use
# Line 611  c*************************************** Line 672  c***************************************
672    
673        subroutine Ri_iwmix (        subroutine Ri_iwmix (
674       I       kmtj, shsq, dbloc, dblocSm       I       kmtj, shsq, dbloc, dblocSm
675       I     , ikey       I     , diffusKzS, diffusKzT
676         I     , ikppkey
677       O     , diffus )       O     , diffus )
678    
679  c     compute interior viscosity diffusivity coefficients due  c     compute interior viscosity diffusivity coefficients due
# Line 631  c     kmtj   (imt)         number of ver Line 693  c     kmtj   (imt)         number of ver
693  c     shsq   (imt,Nr)      (local velocity shear)^2               ((m/s)^2)  c     shsq   (imt,Nr)      (local velocity shear)^2               ((m/s)^2)
694  c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)  c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)
695  c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)  c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)
696        integer kmtj   (imt)  c     diffusKzS(imt,Nr)- background vertical diffusivity for scalars    (m^2/s)
697        _KPP_RL shsq   (imt,Nr)  c     diffusKzT(imt,Nr)- background vertical diffusivity for theta      (m^2/s)
698        _KPP_RL dbloc  (imt,Nr)        integer kmtj (imt)
699        _KPP_RL dblocSm(imt,Nr)        _RL shsq     (imt,Nr)
700        integer ikey        _RL dbloc    (imt,Nr)
701          _RL dblocSm  (imt,Nr)
702          _RL diffusKzS(imt,Nr)
703          _RL diffusKzT(imt,Nr)
704          integer ikppkey
705    
706  c output  c output
707  c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)  c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)
708  c     diffus(imt,0:Nrp1,2)  vertical scalar diffusivity             (m^2/s)  c     diffus(imt,0:Nrp1,2)  vertical scalar diffusivity             (m^2/s)
709  c     diffus(imt,0:Nrp1,3)  vertical temperature diffusivity        (m^2/s)  c     diffus(imt,0:Nrp1,3)  vertical temperature diffusivity        (m^2/s)
710        _KPP_RL diffus(imt,0:Nrp1,3)        _RL diffus(imt,0:Nrp1,3)
711    
712  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
713    
714  c local variables  c local variables
715  c     Rig                   local Richardson number  c     Rig                   local Richardson number
716  c     fRi, fcon             function of Rig  c     fRi, fcon             function of Rig
717        _KPP_RL Rig        _RL Rig
718        _KPP_RL fRi, fcon        _RL fRi, fcon
719        _KPP_RL ratio        _RL ratio
720        integer i, ki, mr        integer i, ki
721        _KPP_RL c1, c0        _RL c1, c0
722    
723  #ifdef ALLOW_KPP_VERTICALLY_SMOOTH  #ifdef ALLOW_KPP_VERTICALLY_SMOOTH
724          integer mr
725  CADJ INIT kpp_ri_tape_mr = common, 1  CADJ INIT kpp_ri_tape_mr = common, 1
726  #endif  #endif
727    
# Line 670  c     set values at bottom and below to Line 737  c     set values at bottom and below to
737  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
738  C     break data flow dependence on diffus  C     break data flow dependence on diffus
739        diffus(1,1,1) = 0.0        diffus(1,1,1) = 0.0
740    
741          do ki = 1, Nr
742             do i = 1, imt
743                diffus(i,ki,1) = 0.
744                diffus(i,ki,2) = 0.
745                diffus(i,ki,3) = 0.
746             enddo
747          enddo
748  #endif  #endif
749                
750    
751        do ki = 1, Nr        do ki = 1, Nr
752           do i = 1, imt           do i = 1, imt
# Line 687  C     break data flow dependence on diff Line 763  C     break data flow dependence on diff
763              endif              endif
764           end do           end do
765        end do        end do
766    CADJ store diffus = comlev1_kpp, key = ikppkey
767    
768  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
769  c     vertically smooth Ri  c     vertically smooth Ri
# Line 697  CADJ store diffus(:,:,1) = kpp_ri_tape_m Line 774  CADJ store diffus(:,:,1) = kpp_ri_tape_m
774  CADJ &  , key=mr, shape=(/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2 /)  CADJ &  , key=mr, shape=(/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2 /)
775    
776           call z121 (           call z121 (
777       U     diffus(1,0,1))       U     diffus(1,0,1)
778         I     myThid )
779        end do        end do
780  #endif  #endif
781    
# Line 724  c  evaluate f of smooth Ri for shear ins Line 802  c  evaluate f of smooth Ri for shear ins
802  c ----------------------------------------------------------------------  c ----------------------------------------------------------------------
803  c            evaluate diffusivities and viscosity  c            evaluate diffusivities and viscosity
804  c    mixing due to internal waves, and shear and static instability  c    mixing due to internal waves, and shear and static instability
805    
806              diffus(i,ki,1) = viscAr + fcon * difmcon + fRi * difm0  #ifndef EXCLUDE_KPP_SHEAR_MIX
807              diffus(i,ki,2) = diffKrS + fcon * difscon + fRi * difs0              if ( .NOT. inAdMode ) then
808              diffus(i,ki,3) = diffKrT + fcon * difscon + fRi * difs0                 diffus(i,ki,1) = viscAr + fcon * difmcon + fRi * difm0
809                   diffus(i,ki,2) = diffusKzS(i,ki)+fcon*difscon+fRi*difs0
810                   diffus(i,ki,3) = diffusKzT(i,ki)+fcon*diftcon+fRi*dift0
811                else
812                   diffus(i,ki,1) = viscAr
813                   diffus(i,ki,2) = diffusKzS(i,ki)
814                   diffus(i,ki,3) = diffusKzT(i,ki)
815                endif
816    #else
817                diffus(i,ki,1) = viscAr
818                diffus(i,ki,2) = diffusKzS(i,ki)
819                diffus(i,ki,3) = diffusKzT(i,ki)
820    #endif
821    
822           end do           end do
823        end do        end do
# Line 749  c         set surface values to 0.0 Line 839  c         set surface values to 0.0
839  c*************************************************************************  c*************************************************************************
840    
841        subroutine z121 (        subroutine z121 (
842       U     v )       U     v,
843         I     myThid )
844    
845  c     Apply 121 smoothing in k to 2-d array V(i,k=1,Nr)  c     Apply 121 smoothing in k to 2-d array V(i,k=1,Nr)
846  c     top (0) value is used as a dummy  c     top (0) value is used as a dummy
# Line 767  c     penetrative convection. Line 858  c     penetrative convection.
858  c input/output  c input/output
859  c-------------  c-------------
860  c v     : 2-D array to be smoothed in Nrp1 direction  c v     : 2-D array to be smoothed in Nrp1 direction
861        _KPP_RL v(imt,0:Nrp1)  c myThid: thread number for this instance of the routine
862          integer myThid
863          _RL v(imt,0:Nrp1)
864    
865  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
866    
867  c local  c local
868        _KPP_RL zwork, zflag        _RL zwork, zflag
869        _KPP_RL KRi_range(1:Nrp1)        _RL KRi_range(1:Nrp1)
870        integer i, k, km1, kp1        integer i, k, km1, kp1
871    
872        _KPP_RL     p0      , p25       , p5      , p2        _RL         p0      , p25       , p5      , p2
873        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 )
874    
875        KRi_range(Nrp1) = p0        KRi_range(Nrp1) = p0
# Line 830  CADJ STORE v(i,k), zwork = z121tape Line 923  CADJ STORE v(i,k), zwork = z121tape
923    
924  c*************************************************************************  c*************************************************************************
925    
       subroutine kpp_smooth_horiz (  
      I     k, bi, bj,  
      U     fld )  
   
 c     Apply horizontal smoothing to KPP array  
   
       IMPLICIT NONE  
 #include "SIZE.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   *   pMask(i  ,j  ,k,bi,bj)   +  
      &           p125  * ( pMask(im1,j  ,k,bi,bj)   +  
      &                     pMask(ip1,j  ,k,bi,bj)   +  
      &                     pMask(i  ,jm1,k,bi,bj)   +  
      &                     pMask(i  ,jp1,k,bi,bj) ) +  
      &           p0625 * ( pMask(im1,jm1,k,bi,bj)   +  
      &                     pMask(im1,jp1,k,bi,bj)   +  
      &                     pMask(ip1,jm1,k,bi,bj)   +  
      &                     pMask(ip1,jp1,k,bi,bj) )  
             IF ( tempVar .GE. p25 ) THEN  
                fld_tmp(i,j) = (  
      &              p25  * fld(i  ,j  )*pMask(i  ,j  ,k,bi,bj) +  
      &              p125 *(fld(im1,j  )*pMask(im1,j  ,k,bi,bj) +  
      &                     fld(ip1,j  )*pMask(ip1,j  ,k,bi,bj) +  
      &                     fld(i  ,jm1)*pMask(i  ,jm1,k,bi,bj) +  
      &                     fld(i  ,jp1)*pMask(i  ,jp1,k,bi,bj))+  
      &              p0625*(fld(im1,jm1)*pMask(im1,jm1,k,bi,bj) +  
      &                     fld(im1,jp1)*pMask(im1,jp1,k,bi,bj) +  
      &                     fld(ip1,jm1)*pMask(ip1,jm1,k,bi,bj) +  
      &                     fld(ip1,jp1)*pMask(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*************************************************************************  
   
926        subroutine smooth_horiz (        subroutine smooth_horiz (
927       I     k, bi, bj,       I     k, bi, bj,
928       U     fld )       U     fld,
929         I     myThid )
930    
931  c     Apply horizontal smoothing to global _RL 2-D array  c     Apply horizontal smoothing to global _RL 2-D array
932    
933        IMPLICIT NONE        IMPLICIT NONE
934  #include "SIZE.h"  #include "SIZE.h"
935    #include "GRID.h"
936  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
937    
938  c     input  c     input
939  c     bi, bj : array indices  c     bi, bj : array indices
940  c     k      : vertical index used for masking  c     k      : vertical index used for masking
941    c     myThid : thread number for this instance of the routine
942          INTEGER myThid
943        integer k, bi, bj        integer k, bi, bj
944    
945  c     input/output  c     input/output
# Line 949  c     local Line 966  c     local
966              im1 = i-1              im1 = i-1
967              ip1 = i+1              ip1 = i+1
968              tempVar =              tempVar =
969       &           p25   *   pMask(i  ,j  ,k,bi,bj)   +       &           p25   *   maskC(i  ,j  ,k,bi,bj)   +
970       &           p125  * ( pMask(im1,j  ,k,bi,bj)   +       &           p125  * ( maskC(im1,j  ,k,bi,bj)   +
971       &                     pMask(ip1,j  ,k,bi,bj)   +       &                     maskC(ip1,j  ,k,bi,bj)   +
972       &                     pMask(i  ,jm1,k,bi,bj)   +       &                     maskC(i  ,jm1,k,bi,bj)   +
973       &                     pMask(i  ,jp1,k,bi,bj) ) +       &                     maskC(i  ,jp1,k,bi,bj) ) +
974       &           p0625 * ( pMask(im1,jm1,k,bi,bj)   +       &           p0625 * ( maskC(im1,jm1,k,bi,bj)   +
975       &                     pMask(im1,jp1,k,bi,bj)   +       &                     maskC(im1,jp1,k,bi,bj)   +
976       &                     pMask(ip1,jm1,k,bi,bj)   +       &                     maskC(ip1,jm1,k,bi,bj)   +
977       &                     pMask(ip1,jp1,k,bi,bj) )       &                     maskC(ip1,jp1,k,bi,bj) )
978              IF ( tempVar .GE. p25 ) THEN              IF ( tempVar .GE. p25 ) THEN
979                 fld_tmp(i,j) = (                 fld_tmp(i,j) = (
980       &              p25  * fld(i  ,j  )*pMask(i  ,j  ,k,bi,bj) +       &              p25  * fld(i  ,j  )*maskC(i  ,j  ,k,bi,bj) +
981       &              p125 *(fld(im1,j  )*pMask(im1,j  ,k,bi,bj) +       &              p125 *(fld(im1,j  )*maskC(im1,j  ,k,bi,bj) +
982       &                     fld(ip1,j  )*pMask(ip1,j  ,k,bi,bj) +       &                     fld(ip1,j  )*maskC(ip1,j  ,k,bi,bj) +
983       &                     fld(i  ,jm1)*pMask(i  ,jm1,k,bi,bj) +       &                     fld(i  ,jm1)*maskC(i  ,jm1,k,bi,bj) +
984       &                     fld(i  ,jp1)*pMask(i  ,jp1,k,bi,bj))+       &                     fld(i  ,jp1)*maskC(i  ,jp1,k,bi,bj))+
985       &              p0625*(fld(im1,jm1)*pMask(im1,jm1,k,bi,bj) +       &              p0625*(fld(im1,jm1)*maskC(im1,jm1,k,bi,bj) +
986       &                     fld(im1,jp1)*pMask(im1,jp1,k,bi,bj) +       &                     fld(im1,jp1)*maskC(im1,jp1,k,bi,bj) +
987       &                     fld(ip1,jm1)*pMask(ip1,jm1,k,bi,bj) +       &                     fld(ip1,jm1)*maskC(ip1,jm1,k,bi,bj) +
988       &                     fld(ip1,jp1)*pMask(ip1,jp1,k,bi,bj)))       &                     fld(ip1,jp1)*maskC(ip1,jp1,k,bi,bj)))
989       &              / tempVar       &              / tempVar
990              ELSE              ELSE
991                 fld_tmp(i,j) = fld(i,j)                 fld_tmp(i,j) = fld(i,j)
# Line 992  c*************************************** Line 1009  c***************************************
1009    
1010        subroutine blmix (        subroutine blmix (
1011       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl
1012       O     , dkm1, blmc, ghat, sigma, ikey       O     , dkm1, blmc, ghat, sigma, ikppkey
1013       &     )       I     , mythid )
1014    
1015  c     mixing coefficients within boundary layer depend on surface  c     mixing coefficients within boundary layer depend on surface
1016  c     forcing and the magnitude and gradient of interior mixing below  c     forcing and the magnitude and gradient of interior mixing below
# Line 1015  c     hbl   (imt)                 bounda Line 1032  c     hbl   (imt)                 bounda
1032  c     stable(imt)                 = 1 in stable forcing  c     stable(imt)                 = 1 in stable forcing
1033  c     casea (imt)                 = 1 in case A  c     casea (imt)                 = 1 in case A
1034  c     diffus(imt,0:Nrp1,mdiff)    vertical diffusivities              (m^2/s)  c     diffus(imt,0:Nrp1,mdiff)    vertical diffusivities              (m^2/s)
1035  c     kbl(imt)                    -1 of first grid level below hbl  c     kbl   (imt)                 -1 of first grid level below hbl
1036        _KPP_RL ustar (imt)  c     myThid               thread number for this instance of the routine
1037        _KPP_RL bfsfc (imt)        integer myThid
1038        _KPP_RL hbl   (imt)        _RL ustar (imt)
1039        _KPP_RL stable(imt)        _RL bfsfc (imt)
1040        _KPP_RL casea (imt)        _RL hbl   (imt)
1041        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL stable(imt)
1042          _RL casea (imt)
1043          _RL diffus(imt,0:Nrp1,mdiff)
1044        integer kbl(imt)        integer kbl(imt)
1045    
1046  c output  c output
# Line 1029  c     dkm1 (imt,mdiff)            bounda Line 1048  c     dkm1 (imt,mdiff)            bounda
1048  c     blmc (imt,Nr,mdiff)         boundary layer mixing coefficients  (m^2/s)  c     blmc (imt,Nr,mdiff)         boundary layer mixing coefficients  (m^2/s)
1049  c     ghat (imt,Nr)               nonlocal scalar transport  c     ghat (imt,Nr)               nonlocal scalar transport
1050  c     sigma(imt)                  normalized depth (d / hbl)  c     sigma(imt)                  normalized depth (d / hbl)
1051        _KPP_RL dkm1 (imt,mdiff)        _RL dkm1 (imt,mdiff)
1052        _KPP_RL blmc (imt,Nr,mdiff)        _RL blmc (imt,Nr,mdiff)
1053        _KPP_RL ghat (imt,Nr)        _RL ghat (imt,Nr)
1054        _KPP_RL sigma(imt)        _RL sigma(imt)
1055        integer ikey        integer ikppkey, kkppkey
1056    
1057  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1058    
# Line 1041  c  local Line 1060  c  local
1060  c     gat1*(imt)                 shape function at sigma = 1  c     gat1*(imt)                 shape function at sigma = 1
1061  c     dat1*(imt)                 derivative of shape function at sigma = 1  c     dat1*(imt)                 derivative of shape function at sigma = 1
1062  c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)  c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)
1063        _KPP_RL gat1m(imt), gat1s(imt), gat1t(imt)        _RL gat1m(imt), gat1s(imt), gat1t(imt)
1064        _KPP_RL dat1m(imt), dat1s(imt), dat1t(imt)        _RL dat1m(imt), dat1s(imt), dat1t(imt)
1065        _KPP_RL ws(imt), wm(imt)        _RL ws(imt), wm(imt)
1066        integer i, kn, ki        integer i, kn, ki
1067        _KPP_RL R, dvdzup, dvdzdn, viscp        _RL R, dvdzup, dvdzdn, viscp
1068        _KPP_RL difsp, diftp, visch, difsh, difth        _RL difsp, diftp, visch, difsh, difth
1069        _KPP_RL f1, sig, a1, a2, a3, delhat        _RL f1, sig, a1, a2, a3, delhat
1070        _KPP_RL Gm, Gs, Gt        _RL Gm, Gs, Gt
1071        _KPP_RL tempVar        _RL tempVar
1072    
1073        _KPP_RL    p0    , eins        _RL    p0    , eins
1074        parameter (p0=0.0, eins=1.0)        parameter (p0=0.0, eins=1.0)
1075    
1076  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 1062  c--------------------------------------- Line 1081  c---------------------------------------
1081           sigma(i) = stable(i) * 1.0 + (1. - stable(i)) * epsilon           sigma(i) = stable(i) * 1.0 + (1. - stable(i)) * epsilon
1082        end do        end do
1083    
1084    CADJ STORE sigma = comlev1_kpp, key = ikppkey
1085        call wscale (        call wscale (
1086       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1087       O        wm, ws )       O        wm, ws, myThid )
1088    CADJ STORE wm = comlev1_kpp, key = ikppkey
1089    CADJ STORE ws = comlev1_kpp, key = ikppkey
1090    
1091        do i = 1, imt        do i = 1, imt
1092           wm(i) = sign(eins,wm(i))*max(phepsi,abs(wm(i)))           wm(i) = sign(eins,wm(i))*max(phepsi,abs(wm(i)))
1093           ws(i) = sign(eins,ws(i))*max(phepsi,abs(ws(i)))           ws(i) = sign(eins,ws(i))*max(phepsi,abs(ws(i)))
1094        end do        end do
1095  CADJ STORE wm = comlev1_kpp, key = ikey  CADJ STORE wm = comlev1_kpp, key = ikppkey
1096  CADJ STORE ws = comlev1_kpp, key = ikey  CADJ STORE ws = comlev1_kpp, key = ikppkey
1097    
1098        do i = 1, imt        do i = 1, imt
1099    
# Line 1107  c--------------------------------------- Line 1129  c---------------------------------------
1129       &        max(ustar(i)**4,phepsi)       &        max(ustar(i)**4,phepsi)
1130           gat1m(i) = visch / hbl(i) / wm(i)           gat1m(i) = visch / hbl(i) / wm(i)
1131           dat1m(i) = -viscp / wm(i) + f1 * visch           dat1m(i) = -viscp / wm(i) + f1 * visch
          dat1m(i) = min(dat1m(i),p0)  
1132    
1133           gat1s(i) = difsh  / hbl(i) / ws(i)           gat1s(i) = difsh  / hbl(i) / ws(i)
1134           dat1s(i) = -difsp / ws(i) + f1 * difsh           dat1s(i) = -difsp / ws(i) + f1 * difsh
          dat1s(i) = min(dat1s(i),p0)  
1135    
1136           gat1t(i) = difth /  hbl(i) / ws(i)           gat1t(i) = difth /  hbl(i) / ws(i)
1137           dat1t(i) = -diftp / ws(i) + f1 * difth           dat1t(i) = -diftp / ws(i) + f1 * difth
          dat1t(i) = min(dat1t(i),p0)  
1138    
1139        end do        end do
1140    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1141    CADJ STORE gat1m = comlev1_kpp, key = ikppkey
1142    CADJ STORE gat1s = comlev1_kpp, key = ikppkey
1143    CADJ STORE gat1t = comlev1_kpp, key = ikppkey
1144    CADJ STORE dat1m = comlev1_kpp, key = ikppkey
1145    CADJ STORE dat1s = comlev1_kpp, key = ikppkey
1146    CADJ STORE dat1t = comlev1_kpp, key = ikppkey
1147    #endif
1148          do i = 1, imt
1149             dat1m(i) = min(dat1m(i),p0)
1150             dat1s(i) = min(dat1s(i),p0)
1151             dat1t(i) = min(dat1t(i),p0)
1152          end do
1153    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1154    CADJ STORE dat1m = comlev1_kpp, key = ikppkey
1155    CADJ STORE dat1s = comlev1_kpp, key = ikppkey
1156    CADJ STORE dat1t = comlev1_kpp, key = ikppkey
1157    #endif
1158    
1159        do ki = 1, Nr        do ki = 1, Nr
1160    
1161    #ifdef ALLOW_AUTODIFF_TAMC
1162             kkppkey = (ikppkey-1)*Nr + ki
1163    #endif
1164    
1165  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1166  c     compute turbulent velocity scales on the interfaces  c     compute turbulent velocity scales on the interfaces
1167  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 1129  c--------------------------------------- Line 1170  c---------------------------------------
1170              sig      = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i)              sig      = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i)
1171              sigma(i) = stable(i)*sig + (1.-stable(i))*min(sig,epsilon)              sigma(i) = stable(i)*sig + (1.-stable(i))*min(sig,epsilon)
1172           end do           end do
1173    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1174    CADJ STORE wm = comlev1_kpp_k, key = kkppkey
1175    CADJ STORE ws = comlev1_kpp_k, key = kkppkey
1176    #endif
1177    CADJ STORE sigma = comlev1_kpp_k, key = kkppkey
1178           call wscale (           call wscale (
1179       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1180       O        wm, ws )       O        wm, ws, myThid )
1181    CADJ STORE wm = comlev1_kpp_k, key = kkppkey
1182    CADJ STORE ws = comlev1_kpp_k, key = kkppkey
1183    
1184  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1185  c     compute the dimensionless shape functions at the interfaces  c     compute the dimensionless shape functions at the interfaces
# Line 1175  c--------------------------------------- Line 1223  c---------------------------------------
1223       &            + (1. - stable(i)) * min(sig,epsilon)       &            + (1. - stable(i)) * min(sig,epsilon)
1224        end do        end do
1225    
1226    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1227    CADJ STORE wm = comlev1_kpp, key = ikppkey
1228    CADJ STORE ws = comlev1_kpp, key = ikppkey
1229    #endif
1230    CADJ STORE sigma = comlev1_kpp, key = ikppkey
1231        call wscale (        call wscale (
1232       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1233       O        wm, ws )       O        wm, ws, myThid )
1234    CADJ STORE wm = comlev1_kpp, key = ikppkey
1235    CADJ STORE ws = comlev1_kpp, key = ikppkey
1236    
1237        do i = 1, imt        do i = 1, imt
1238           sig = -zgrid(kbl(i)-1) / hbl(i)           sig = -zgrid(kbl(i)-1) / hbl(i)
# Line 1203  c*************************************** Line 1258  c***************************************
1258       I       dkm1, hbl, kbl, diffus, casea       I       dkm1, hbl, kbl, diffus, casea
1259       U     , ghat       U     , ghat
1260       O     , blmc       O     , blmc
1261       &     )       &     , myThid )
1262    
1263  c enhance the diffusivity at the kbl-.5 interface  c enhance the diffusivity at the kbl-.5 interface
1264    
# Line 1218  c     hbl(imt)                  boundary Line 1273  c     hbl(imt)                  boundary
1273  c     kbl(imt)                  grid above hbl  c     kbl(imt)                  grid above hbl
1274  c     diffus(imt,0:Nrp1,mdiff) vertical diffusivities           (m^2/s)  c     diffus(imt,0:Nrp1,mdiff) vertical diffusivities           (m^2/s)
1275  c     casea(imt)                = 1 in caseA, = 0 in case B  c     casea(imt)                = 1 in caseA, = 0 in case B
1276        _KPP_RL dkm1  (imt,mdiff)  c     myThid                    thread number for this instance of the routine
1277        _KPP_RL hbl   (imt)        integer   myThid
1278          _RL dkm1  (imt,mdiff)
1279          _RL hbl   (imt)
1280        integer kbl   (imt)        integer kbl   (imt)
1281        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL diffus(imt,0:Nrp1,mdiff)
1282        _KPP_RL casea (imt)        _RL casea (imt)
1283    
1284  c input/output  c input/output
1285  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)
1286        _KPP_RL ghat (imt,Nr)        _RL ghat (imt,Nr)
1287    
1288  c output  c output
1289  c     enhanced bound. layer mixing coeff.  c     enhanced bound. layer mixing coeff.
1290        _KPP_RL blmc  (imt,Nr,mdiff)        _RL blmc  (imt,Nr,mdiff)
1291    
1292  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1293    
1294  c local  c local
1295  c     fraction hbl lies beteen zgrid neighbors  c     fraction hbl lies beteen zgrid neighbors
1296        _KPP_RL delta        _RL delta
1297        integer ki, i, md        integer ki, i, md
1298        _KPP_RL dkmp5, dstar        _RL dkmp5, dstar
1299    
1300        do i = 1, imt        do i = 1, imt
1301           ki = kbl(i)-1           ki = kbl(i)-1
# Line 1264  c     fraction hbl lies beteen zgrid nei Line 1321  c     fraction hbl lies beteen zgrid nei
1321  c*************************************************************************  c*************************************************************************
1322    
1323        SUBROUTINE STATEKPP (        SUBROUTINE STATEKPP (
1324       I     bi, bj, myThid,       O     RHO1, DBLOC, DBSFC, TTALPHA, SSBETA,
1325       O     RHO1, DBLOC, DBSFC, TTALPHA, SSBETA)       I     ikppkey, bi, bj, myThid )
1326  c  c
1327  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1328  c     "statekpp" computes all necessary input arrays  c     "statekpp" computes all necessary input arrays
# Line 1290  c Line 1347  c
1347  c     written  by: jan morzel,   feb. 10, 1995 (converted from "sigma" version)  c     written  by: jan morzel,   feb. 10, 1995 (converted from "sigma" version)
1348  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
1349  c  c
1350    c     28 april 05: added computation of mixed layer depth KPPmld
1351    c     for a deltaRho equivalent to a deltaTheta of 0.8 deg C
1352    
1353  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1354    
1355        IMPLICIT NONE        IMPLICIT NONE
# Line 1299  c--------------------------------------- Line 1359  c---------------------------------------
1359  #include "PARAMS.h"  #include "PARAMS.h"
1360  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
1361  #include "DYNVARS.h"  #include "DYNVARS.h"
1362    #include "GRID.h"
1363    
1364  c-------------- Routine arguments -----------------------------------------  c-------------- Routine arguments -----------------------------------------
1365        INTEGER bi, bj, myThid        INTEGER bi, bj, myThid
1366        _KPP_RL RHO1   ( ibot:itop, jbot:jtop       )        _RL RHO1   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy       )
1367        _KPP_RL DBLOC  ( ibot:itop, jbot:jtop, Nr   )        _RL DBLOC  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
1368        _KPP_RL DBSFC  ( ibot:itop, jbot:jtop, Nr   )        _RL DBSFC  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
1369        _KPP_RL TTALPHA( ibot:itop, jbot:jtop, Nrp1 )        _RL TTALPHA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1370        _KPP_RL SSBETA ( ibot:itop, jbot:jtop, Nrp1 )        _RL SSBETA ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1371    
1372  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1373    
# Line 1317  c Line 1378  c
1378  c     rhok         - density of t(k  ) & s(k  ) at depth k  c     rhok         - density of t(k  ) & s(k  ) at depth k
1379  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
1380  c     rho1k        - density of t(1  ) & s(1  ) at depth k  c     rho1k        - density of t(1  ) & s(1  ) at depth k
1381  c     work1, work2 - work arrays for holding horizontal slabs  c     work1,2,3    - work arrays for holding horizontal slabs
1382    
1383        _RL RHOK  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL RHOK  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1384        _RL RHOKM1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL RHOKM1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 1325  c     work1, work2 - work arrays for hol Line 1386  c     work1, work2 - work arrays for hol
1386        _RL WORK1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1387        _RL WORK2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1388        _RL WORK3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1389    #ifdef ALLOW_DIAGNOSTICS
1390    c     KPPMLD - mixed layer depth based on density criterion
1391          _RL KPPMLD(1    :sNx    ,1    :sNy    )
1392    #endif /* ALLOW_DIAGNOSTICS */
1393    
1394        INTEGER I, J, K        INTEGER I, J, K
1395          INTEGER ikppkey, kkppkey
1396    
1397  c calculate density, alpha, beta in surface layer, and set dbsfc to zero  c calculate density, alpha, beta in surface layer, and set dbsfc to zero
1398    
1399          kkppkey = (ikppkey-1)*Nr + 1
1400    
1401    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1402    CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1403    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1404    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1405        call FIND_RHO(        call FIND_RHO(
1406       I     bi, bj, ibot, itop, jbot, jtop, 1, 1, eosType,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1407       I     theta, salt,       I     theta, salt,
1408       O     WORK1,       O     WORK1,
1409       I     myThid )       I     myThid )
1410    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1411    CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1412    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1413    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1414    
1415        call FIND_ALPHA(        call FIND_ALPHA(
1416       I     bi, bj, ibot, itop, jbot, jtop, 1, 1, eosType,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1417       O     WORK2 )       O     WORK2, myThid )
1418    
1419        call FIND_BETA(        call FIND_BETA(
1420       I     bi, bj, ibot, itop, jbot, jtop, 1, 1, eosType,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1421       O     WORK3 )       O     WORK3, myThid )
1422    
1423        DO J = jbot, jtop        DO J = 1-OLy, sNy+OLy
1424           DO I = ibot, itop           DO I = 1-OLx, sNx+OLx
1425              RHO1(I,J)      = WORK1(I,J) + rhonil              RHO1(I,J)      = WORK1(I,J) + rhoConst
1426              TTALPHA(I,J,1) = WORK2(I,J)              TTALPHA(I,J,1) = WORK2(I,J)
1427              SSBETA(I,J,1)  = WORK3(I,J)              SSBETA(I,J,1)  = WORK3(I,J)
1428              DBSFC(I,J,1)   = 0.              DBSFC(I,J,1)   = 0.
1429           END DO           END DO
1430        END DO        END DO
1431    
1432    #ifdef ALLOW_DIAGNOSTICS
1433    c     work3 - density of t(1)-.8 & s(1  ) at depth 1
1434          IF ( useDiagnostics ) THEN
1435             DO J = 1, sNy
1436                DO I = 1, sNx
1437                   KPPMLD(I,J) = ABS(R_low(I,J,bi,bj))
1438                   WORK3 (I,J) = WORK1(I,J) - 0.8 * WORK2(I,J)
1439                END DO
1440             END DO
1441          ENDIF
1442    #endif /* ALLOW_DIAGNOSTICS */
1443    
1444  c calculate alpha, beta, and gradients in interior layers  c calculate alpha, beta, and gradients in interior layers
1445    
1446  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)
1447        DO K = 2, Nr        DO K = 2, Nr
1448    
1449          kkppkey = (ikppkey-1)*Nr + k
1450    
1451    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1452    CADJ STORE theta(:,:,k,bi,bj) = comlev1_kpp_k, key=kkppkey
1453    CADJ STORE salt (:,:,k,bi,bj) = comlev1_kpp_k, key=kkppkey
1454    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1455           call FIND_RHO(           call FIND_RHO(
1456       I        bi, bj, ibot, itop, jbot, jtop, K, K, eosType,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1457       I        theta, salt,       I        theta, salt,
1458       O        RHOK,       O        RHOK,
1459       I        myThid )       I        myThid )
1460    
1461    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1462    CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_kpp_k, key=kkppkey
1463    CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_kpp_k, key=kkppkey
1464    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1465           call FIND_RHO(           call FIND_RHO(
1466       I        bi, bj, ibot, itop, jbot, jtop, K-1, K, eosType,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K-1, K,
1467       I        theta, salt,       I        theta, salt,
1468       O        RHOKM1,       O        RHOKM1,
1469       I        myThid )       I        myThid )
1470    
1471    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1472    CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1473    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1474    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1475           call FIND_RHO(           call FIND_RHO(
1476       I        bi, bj, ibot, itop, jbot, jtop, 1, K, eosType,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, K,
1477       I        theta, salt,       I        theta, salt,
1478       O        RHO1K,       O        RHO1K,
1479       I        myThid )       I        myThid )
1480    
1481    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1482    CADJ STORE rhok  (:,:)          = comlev1_kpp_k, key=kkppkey
1483    CADJ STORE rhokm1(:,:)          = comlev1_kpp_k, key=kkppkey
1484    CADJ STORE rho1k (:,:)          = comlev1_kpp_k, key=kkppkey
1485    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1486    
1487           call FIND_ALPHA(           call FIND_ALPHA(
1488       I        bi, bj, ibot, itop, jbot, jtop, K, K, eosType,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1489       O        WORK1 )       O        WORK1, myThid )
1490    
1491           call FIND_BETA(           call FIND_BETA(
1492       I        bi, bj, ibot, itop, jbot, jtop, K, K, eosType,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1493       O        WORK2 )       O        WORK2, myThid )
1494    
1495           DO J = jbot, jtop           DO J = 1-OLy, sNy+OLy
1496              DO I = ibot, itop              DO I = 1-OLx, sNx+OLx
1497                 TTALPHA(I,J,K) = WORK1 (I,J)                 TTALPHA(I,J,K) = WORK1 (I,J)
1498                 SSBETA(I,J,K)  = WORK2 (I,J)                 SSBETA(I,J,K)  = WORK2 (I,J)
1499                 DBLOC(I,J,K-1) = gravity * (RHOK(I,J) - RHOKM1(I,J)) /                 DBLOC(I,J,K-1) = gravity * (RHOK(I,J) - RHOKM1(I,J)) /
1500       &                                    (RHOK(I,J) + rhonil)       &                                    (RHOK(I,J) + rhoConst)
1501                 DBSFC(I,J,K)   = gravity * (RHOK(I,J) - RHO1K (I,J)) /                 DBSFC(I,J,K)   = gravity * (RHOK(I,J) - RHO1K (I,J)) /
1502       &                                    (RHOK(I,J) + rhonil)       &                                    (RHOK(I,J) + rhoConst)
1503              END DO              END DO
1504           END DO           END DO
1505    
1506    #ifdef ALLOW_DIAGNOSTICS
1507             IF ( useDiagnostics ) THEN
1508    c     work1 - density of t(k-1)  & s(k-1) at depth 1
1509    c     work2 - density of t(k  )  & s(k  ) at depth 1
1510    c     work3 - density of t(1)-.8 & s(1  ) at depth 1
1511                call FIND_RHO(
1512         I           bi, bj, 1, sNx, 1, sNy, K-1, 1, theta, salt,
1513         O           WORK1,
1514         I           myThid )
1515                call FIND_RHO(
1516         I           bi, bj, 1, sNx, 1, sNy, K  , 1, theta, salt,
1517         O           WORK2,
1518         I           myThid )
1519                DO J = 1, sNy
1520                   DO I = 1, sNx
1521                      IF ( k .LE. klowC(I,J,bi,bj) .AND.
1522         &                 WORK1(I,J) .LT. WORK3(I,J) .AND.
1523         &                 WORK2(I,J) .GE. WORK3(I,J) ) THEN
1524                         KPPMLD(I,J) = MIN ( KPPMLD(I,J),
1525         &                    ABS(((WORK3(I,J)-WORK1(I,J))*rC(k)+
1526         &                    (WORK2(I,J)-WORK3(I,J))*rC(k-1))/
1527         &                    (WORK2(I,J)-WORK1(I,J))))
1528                      ENDIF
1529                   END DO
1530                END DO
1531             ENDIF
1532    #endif /* ALLOW_DIAGNOSTICS */
1533    
1534        END DO        END DO
1535    
1536  c     compute arrays for K = Nrp1  c     compute arrays for K = Nrp1
1537        DO J = jbot, jtop        DO J = 1-OLy, sNy+OLy
1538           DO I = ibot, itop           DO I = 1-OLx, sNx+OLx
1539              TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)              TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)
1540              SSBETA(I,J,Nrp1)  = SSBETA(I,J,Nr)              SSBETA(I,J,Nrp1)  = SSBETA(I,J,Nr)
1541              DBLOC(I,J,Nr)     = 0.              DBLOC(I,J,Nr)     = 0.
1542           END DO           END DO
1543        END DO        END DO
1544    
1545    #ifdef ALLOW_DIAGNOSTICS
1546          IF ( useDiagnostics ) THEN
1547             CALL DIAGNOSTICS_FILL(KPPmld,'KPPmld  ',0,1,3,bi,bj,myThid)
1548          ENDIF
1549    #endif /* ALLOW_DIAGNOSTICS */
1550    
1551  #endif /* ALLOW_KPP */  #endif /* ALLOW_KPP */
1552    
1553        RETURN        RETURN
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.30
  ViewVC Help
Powered by ViewVC 1.1.22