/[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.13 by dimitri, Sat Dec 28 10:11:11 2002 UTC revision 1.31 by jmc, Thu May 3 21:36:26 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, myIter, 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   myIter :: Current iteration number in simulation
55  c     kmtj   (imt)    - number of vertical layers on this row  c   myThid :: My Thread Id. number
56  c     shsq   (imt,Nr) - (local velocity shear)^2                      ((m/s)^2)  c     kmtj   (imt)     - number of vertical layers on this row
57  c     dvsq   (imt,Nr) - (velocity shear re sfc)^2                     ((m/s)^2)  c     msk    (imt)     - surface mask (=1 if water, =0 otherwise)
58  c     ustar  (imt)    - surface friction velocity                         (m/s)  c     shsq   (imt,Nr)  - (local velocity shear)^2                     ((m/s)^2)
59  c     bo     (imt)    - surface turbulent buoy. forcing               (m^2/s^3)  c     dvsq   (imt,Nr)  - (velocity shear re sfc)^2                    ((m/s)^2)
60  c     bosol  (imt)    - radiative buoyancy forcing                    (m^2/s^3)  c     ustar  (imt)     - surface friction velocity                        (m/s)
61  c     dbloc  (imt,Nr) - local delta buoyancy across interfaces          (m/s^2)  c     bo     (imt)     - surface turbulent buoy. forcing              (m^2/s^3)
62  c     dblocSm(imt,Nr) - horizontally smoothed dbloc                     (m/s^2)  c     bosol  (imt)     - radiative buoyancy forcing                   (m^2/s^3)
63  c                         stored in ghat to save space  c     dbloc  (imt,Nr)  - local delta buoyancy across interfaces         (m/s^2)
64  c     Ritop  (imt,Nr) - numerator of bulk Richardson Number  c     dblocSm(imt,Nr)  - horizontally smoothed dbloc                    (m/s^2)
65  c                         (zref-z) * delta buoyancy w.r.t. surface    ((m/s)^2)  c                          stored in ghat to save space
66  c     coriol (imt)    - Coriolis parameter                                (1/s)  c     Ritop  (imt,Nr)  - numerator of bulk Richardson Number
67    c                          (zref-z) * delta buoyancy w.r.t. surface   ((m/s)^2)
68    c     coriol (imt)     - Coriolis parameter                               (1/s)
69    c     diffusKzS(imt,Nr)- background vertical diffusivity for scalars    (m^2/s)
70    c     diffusKzT(imt,Nr)- background vertical diffusivity for theta      (m^2/s)
71  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,
72  c           e.g., hbl(sNx,sNy) -> hbl(imt),  c           e.g., hbl(sNx,sNy) -> hbl(imt),
73  c           where hbl(i,j) -> hbl((j-1)*sNx+i)  c           where hbl(i,j) -> hbl((j-1)*sNx+i)
74    
75        _RL     mytime        _RL     myTime
76        integer mythid        integer myIter
77        integer kmtj  (imt   )        integer myThid
78        _KPP_RL shsq  (imt,Nr)        integer kmtj (imt   )
79        _KPP_RL dvsq  (imt,Nr)        _RL shsq     (imt,Nr)
80        _KPP_RL ustar (imt   )        _RL dvsq     (imt,Nr)
81        _KPP_RL bo    (imt   )        _RL ustar    (imt   )
82        _KPP_RL bosol (imt   )        _RL bo       (imt   )
83        _KPP_RL dbloc (imt,Nr)        _RL bosol    (imt   )
84        _KPP_RL Ritop (imt,Nr)        _RL dbloc    (imt,Nr)
85        _KPP_RL coriol(imt   )        _RL Ritop    (imt,Nr)
86          _RL coriol   (imt   )
87          _RS msk      (imt   )
88          _RL diffusKzS(imt,Nr)
89          _RL diffusKzT(imt,Nr)
90    
91        integer ikey        integer ikppkey
92    
93  c output  c output
94  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 97  c     diffus (imt,3)  - vertical tempera
97  c     ghat   (imt)    - nonlocal transport coefficient                  (s/m^2)  c     ghat   (imt)    - nonlocal transport coefficient                  (s/m^2)
98  c     hbl    (imt)    - mixing layer depth                                  (m)  c     hbl    (imt)    - mixing layer depth                                  (m)
99    
100        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL diffus(imt,0:Nrp1,mdiff)
101        _KPP_RL ghat  (imt,Nr)        _RL ghat  (imt,Nr)
102        _KPP_RL hbl   (imt)        _RL hbl   (imt)
103    
104  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
105    
# Line 107  c     blmc   (imt,Nr,mdiff) - boundary l Line 113  c     blmc   (imt,Nr,mdiff) - boundary l
113  c     sigma  (imt         ) - normalized depth (d / hbl)  c     sigma  (imt         ) - normalized depth (d / hbl)
114  c     Rib    (imt,Nr      ) - bulk Richardson number  c     Rib    (imt,Nr      ) - bulk Richardson number
115    
116        integer kbl   (imt         )        integer kbl(imt         )
117        _KPP_RL bfsfc (imt         )        _RL bfsfc  (imt         )
118        _KPP_RL casea (imt         )        _RL casea  (imt         )
119        _KPP_RL stable(imt         )        _RL stable (imt         )
120        _KPP_RL dkm1  (imt,   mdiff)        _RL dkm1   (imt,   mdiff)
121        _KPP_RL blmc  (imt,Nr,mdiff)        _RL blmc   (imt,Nr,mdiff)
122        _KPP_RL sigma (imt         )        _RL sigma  (imt         )
123        _KPP_RL Rib   (imt,Nr      )        _RL Rib    (imt,Nr      )
124    
125        integer i, k, md        integer i, k, md
126    
# Line 125  c instability. Line 131  c instability.
131  c (ghat is temporary storage for horizontally smoothed dbloc)  c (ghat is temporary storage for horizontally smoothed dbloc)
132  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
133    
134  CADJ STORE ghat = comlev1_kpp, key = ikey  cph(
135    cph these storings avoid recomp. of Ri_iwmix
136    CADJ STORE ghat  = comlev1_kpp, key = ikppkey
137    CADJ STORE dbloc = comlev1_kpp, key = ikppkey
138    cph)
139        call Ri_iwmix (        call Ri_iwmix (
140       I       kmtj, shsq, dbloc, ghat       I       kmtj, shsq, dbloc, ghat
141       I     , ikey       I     , diffusKzS, diffusKzT
142         I     , ikppkey
143       O     , diffus )       O     , diffus )
144    
145    cph(
146    cph these storings avoid recomp. of Ri_iwmix
147    cph DESPITE TAFs 'not necessary' warning!
148    CADJ STORE dbloc  = comlev1_kpp, key = ikppkey
149    CADJ STORE shsq   = comlev1_kpp, key = ikppkey
150    CADJ STORE ghat   = comlev1_kpp, key = ikppkey
151    CADJ STORE diffus = comlev1_kpp, key = ikppkey
152    cph)
153    
154  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
155  c set seafloor values to zero and fill extra "Nrp1" coefficients  c set seafloor values to zero and fill extra "Nrp1" coefficients
156  c for blmix  c for blmix
# Line 152  c diagnose the new boundary layer depth Line 171  c diagnose the new boundary layer depth
171  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
172    
173        call bldepth (        call bldepth (
174       I       mytime, mythid       I       kmtj
      I     , kmtj  
175       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol
176       I     , ikey       I     , ikppkey
177       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma
178       &     )       I     , myTime, myIter, myThid )
179    
180  CADJ STORE hbl,bfsfc,stable,casea,kbl = comlev1_kpp, key = ikey  CADJ STORE hbl,bfsfc,stable,casea,kbl = comlev1_kpp, key = ikppkey
181    
182  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
183  c compute boundary layer diffusivities  c compute boundary layer diffusivities
# Line 167  c--------------------------------------- Line 185  c---------------------------------------
185    
186        call blmix (        call blmix (
187       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl
188       O     , dkm1, blmc, ghat, sigma, ikey       O     , dkm1, blmc, ghat, sigma, ikppkey
189       &     )       I     , myThid )
190    cph(
191  CADJ STORE dkm1,blmc,ghat = comlev1_kpp, key = ikey  CADJ STORE dkm1,blmc,ghat = comlev1_kpp, key = ikppkey
192    CADJ STORE hbl, kbl, diffus, casea = comlev1_kpp, key = ikppkey
193    cph)
194    
195  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
196  c enhance diffusivity at interface kbl - 1  c enhance diffusivity at interface kbl - 1
# Line 179  c--------------------------------------- Line 199  c---------------------------------------
199        call enhance (        call enhance (
200       I       dkm1, hbl, kbl, diffus, casea       I       dkm1, hbl, kbl, diffus, casea
201       U     , ghat       U     , ghat
202       O     , blmc )       O     , blmc
203         I     , myThid )
204    
205    cph(
206    cph avoids recomp. of enhance
207    CADJ STORE blmc = comlev1_kpp, key = ikppkey
208    cph)
209    
210  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
211  c combine interior and boundary layer coefficients and nonlocal term  c combine interior and boundary layer coefficients and nonlocal term
212    c !!!NOTE!!! In shallow (2-level) regions and for shallow mixed layers
213    c (< 1 level), diffusivity blmc can become negative.  The max-s below
214    c are a hack until this problem is properly diagnosed and fixed.
215  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
   
216        do k = 1, Nr        do k = 1, Nr
217           do i = 1, imt           do i = 1, imt
218              if (k .lt. kbl(i)) then              if (k .lt. kbl(i)) then
219                 do md = 1, mdiff  #ifdef ALLOW_SHELFICE
220                    diffus(i,k,md) = blmc(i,k,md)  C     when there is shelfice on top (msk(i)=0), reset the boundary layer
221                 end do  C     mixing coefficients blmc to pure Ri-number based mixing
222                   blmc(i,k,1) = max ( blmc(i,k,1)*msk(i),
223         &              diffus(i,k,1) )
224                   blmc(i,k,2) = max ( blmc(i,k,2)*msk(i),
225         &              diffus(i,k,2) )
226                   blmc(i,k,3) = max ( blmc(i,k,3)*msk(i),
227         &              diffus(i,k,3) )
228    #endif /* not ALLOW_SHELFICE */
229                   diffus(i,k,1) = max ( blmc(i,k,1), viscAr  )
230                   diffus(i,k,2) = max ( blmc(i,k,2), diffusKzS(i,Nr) )
231                   diffus(i,k,3) = max ( blmc(i,k,3), diffusKzT(i,Nr) )
232              else              else
233                 ghat(i,k) = 0.                 ghat(i,k) = 0.
234              endif              endif
# Line 205  c--------------------------------------- Line 243  c---------------------------------------
243  c*************************************************************************  c*************************************************************************
244    
245        subroutine bldepth (        subroutine bldepth (
246       I       mytime, mythid       I       kmtj
      I     , kmtj  
247       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol
248       I     , ikey       I     , ikppkey
249       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma
250       &     )       I     , myTime, myIter, myThid )
251    
252  c     the oceanic planetary boundary layer depth, hbl, is determined as  c     the oceanic planetary boundary layer depth, hbl, is determined as
253  c     the shallowest depth where the bulk Richardson number is  c     the shallowest depth where the bulk Richardson number is
# Line 236  c Line 273  c
273        IMPLICIT NONE        IMPLICIT NONE
274    
275  #include "SIZE.h"  #include "SIZE.h"
 #include "EEPARAMS.h"  
 #include "PARAMS.h"  
276  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
 #include "FFIELDS.h"  
277    
278  c input  c input
279  c------  c------
280  c myTime    : current time in simulation  c   myTime :: Current time in simulation
281  c myThid    : thread number for this instance of the routine  c   myIter :: Current iteration number in simulation
282    c   myThid :: My Thread Id. number
283  c kmtj      : number of vertical layers  c kmtj      : number of vertical layers
284  c dvsq      : (velocity shear re sfc)^2             ((m/s)^2)  c dvsq      : (velocity shear re sfc)^2             ((m/s)^2)
285  c dbloc     : local delta buoyancy across interfaces  (m/s^2)  c dbloc     : local delta buoyancy across interfaces  (m/s^2)
# Line 255  c ustar     : surface friction velocity Line 290  c ustar     : surface friction velocity
290  c bo        : surface turbulent buoyancy forcing    (m^2/s^3)  c bo        : surface turbulent buoyancy forcing    (m^2/s^3)
291  c bosol     : radiative buoyancy forcing            (m^2/s^3)  c bosol     : radiative buoyancy forcing            (m^2/s^3)
292  c coriol    : Coriolis parameter                        (1/s)  c coriol    : Coriolis parameter                        (1/s)
293        _RL     mytime        _RL     myTime
294        integer mythid        integer myIter
295          integer myThid
296        integer kmtj(imt)        integer kmtj(imt)
297        _KPP_RL dvsq  (imt,Nr)        _RL dvsq    (imt,Nr)
298        _KPP_RL dbloc (imt,Nr)        _RL dbloc   (imt,Nr)
299        _KPP_RL Ritop (imt,Nr)        _RL Ritop   (imt,Nr)
300        _KPP_RL ustar (imt)        _RL ustar   (imt)
301        _KPP_RL bo    (imt)        _RL bo      (imt)
302        _KPP_RL bosol (imt)        _RL bosol   (imt)
303        _KPP_RL coriol(imt)        _RL coriol  (imt)
304        integer ikey        integer ikppkey, kkppkey
305    
306  c  output  c  output
307  c--------  c--------
# Line 276  c casea     : =1 in case A, =0 in case B Line 312  c casea     : =1 in case A, =0 in case B
312  c kbl       : -1 of first grid level below hbl  c kbl       : -1 of first grid level below hbl
313  c Rib       : Bulk Richardson number  c Rib       : Bulk Richardson number
314  c sigma     : normalized depth (d/hbl)  c sigma     : normalized depth (d/hbl)
315        _KPP_RL hbl   (imt)        _RL hbl    (imt)
316        _KPP_RL bfsfc (imt)        _RL bfsfc  (imt)
317        _KPP_RL stable(imt)        _RL stable (imt)
318        _KPP_RL casea (imt)        _RL casea  (imt)
319        integer kbl   (imt)        integer kbl(imt)
320        _KPP_RL Rib   (imt,Nr)        _RL Rib    (imt,Nr)
321        _KPP_RL sigma (imt)        _RL sigma  (imt)
322    
323  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
324    
325  c  local  c  local
326  c-------  c-------
327  c wm, ws    : turbulent velocity scales         (m/s)  c wm, ws    : turbulent velocity scales         (m/s)
328        _KPP_RL wm(imt), ws(imt)        _RL wm(imt), ws(imt)
329        _RL     worka(imt)        _RL worka(imt)
330    
331        _KPP_RL bvsq, vtsq, hekman, hmonob, hlimit, tempVar1, tempVar2        _RL bvsq, vtsq, hekman, hmonob, hlimit, tempVar1, tempVar2
332        integer i, kl        integer i, kl
333    
334        _KPP_RL     p5    , eins        _RL         p5    , eins
335        parameter ( p5=0.5, eins=1.0 )        parameter ( p5=0.5, eins=1.0 )
336        _RL         minusone        _RL         minusone
337        parameter ( minusone=-1.0 )        parameter ( minusone=-1.0 )
# Line 319  c     initialize hbl and kbl to bottomed Line 355  c     initialize hbl and kbl to bottomed
355    
356        do kl = 2, Nr        do kl = 2, Nr
357    
358    #ifdef ALLOW_AUTODIFF_TAMC
359             kkppkey = (ikppkey-1)*Nr + kl
360    #endif
361    
362  c     compute bfsfc = sw fraction at hbf * zgrid  c     compute bfsfc = sw fraction at hbf * zgrid
363    
364           do i = 1, imt           do i = 1, imt
365              worka(i) = zgrid(kl)              worka(i) = zgrid(kl)
366           end do           end do
367    CADJ store worka = comlev1_kpp_k, key = kkppkey
368           call SWFRAC(           call SWFRAC(
369       I        imt, hbf,       I       imt, hbf,
370       I        mytime, mythid,       U       worka,
371       U        worka )       I       myTime, myIter, myThid )
372    CADJ store worka = comlev1_kpp_k, key = kkppkey
373    
374    
375           do i = 1, imt           do i = 1, imt
376    
# Line 349  c--------------------------------------- Line 392  c---------------------------------------
392    
393           call wscale (           call wscale (
394       I        sigma, casea, ustar, bfsfc,       I        sigma, casea, ustar, bfsfc,
395       O        wm, ws )       O        wm, ws, myThid )
396    CADJ store ws = comlev1_kpp_k, key = kkppkey
397    
398           do i = 1, imt           do i = 1, imt
399    
# Line 387  c Line 431  c
431        end do        end do
432    
433  cph(  cph(
434  cph  without this store, there's a recomputation error for  cph  without this store, there is a recomputation error for
435  cph  rib in adbldepth (probably partial recomputation problem)      cph  rib in adbldepth (probably partial recomputation problem)    
436  CADJ store Rib = comlev1_kpp  CADJ store Rib = comlev1_kpp
437  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy),Nr /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy),Nr /)
438  cph)  cph)
439    
440        do kl = 2, Nr        do kl = 2, Nr
# Line 400  cph) Line 444  cph)
444        end do        end do
445    
446  CADJ store kbl = comlev1_kpp  CADJ store kbl = comlev1_kpp
447  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
448    
449        do i = 1, imt        do i = 1, imt
450           kl = kbl(i)           kl = kbl(i)
# Line 413  c     linearly interpolate to find hbl w Line 457  c     linearly interpolate to find hbl w
457        end do        end do
458    
459  CADJ store hbl = comlev1_kpp  CADJ store hbl = comlev1_kpp
460  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
461    
462  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
463  c     find stability and buoyancy forcing for boundary layer  c     find stability and buoyancy forcing for boundary layer
# Line 422  c--------------------------------------- Line 466  c---------------------------------------
466        do i = 1, imt        do i = 1, imt
467           worka(i) = hbl(i)           worka(i) = hbl(i)
468        end do        end do
469    CADJ store worka = comlev1_kpp
470    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
471        call SWFRAC(        call SWFRAC(
472       I     imt, minusone,       I       imt, minusone,
473       I     mytime, mythid,       U       worka,
474       U     worka )       I       myTime, myIter, myThid )
475    CADJ store worka = comlev1_kpp
476    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
477    
478        do i = 1, imt        do i = 1, imt
479           bfsfc(i)  = bo(i) + bosol(i) * (1. - worka(i))           bfsfc(i)  = bo(i) + bosol(i) * (1. - worka(i))
480        end do        end do
481  CADJ store bfsfc = comlev1_kpp  CADJ store bfsfc = comlev1_kpp
482  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
483    
484  c--   ensure bfsfc is never 0  c--   ensure bfsfc is never 0
485        do i = 1, imt        do i = 1, imt
# Line 442  c--   ensure bfsfc is never 0 Line 490  c--   ensure bfsfc is never 0
490  cph(  cph(
491  cph  added stable to store list to avoid extensive recomp.  cph  added stable to store list to avoid extensive recomp.
492  CADJ store bfsfc, stable = comlev1_kpp  CADJ store bfsfc, stable = comlev1_kpp
493  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
494  cph)  cph)
495    
496  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 461  c--------------------------------------- Line 509  c---------------------------------------
509           end if           end if
510        end do        end do
511  CADJ store hbl = comlev1_kpp  CADJ store hbl = comlev1_kpp
512  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
513    
514        do i = 1, imt        do i = 1, imt
515           hbl(i) = max(hbl(i),minKPPhbl)           hbl(i) = max(hbl(i),minKPPhbl)
# Line 469  CADJ &   , key = ikey, shape = (/ (sNx+2 Line 517  CADJ &   , key = ikey, shape = (/ (sNx+2
517        end do        end do
518    
519  CADJ store hbl = comlev1_kpp  CADJ store hbl = comlev1_kpp
520  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
521    
522  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
523  c      find new kbl  c      find new kbl
# Line 490  c--------------------------------------- Line 538  c---------------------------------------
538        do i = 1, imt        do i = 1, imt
539           worka(i) = hbl(i)           worka(i) = hbl(i)
540        end do        end do
541    CADJ store worka = comlev1_kpp
542    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
543        call SWFRAC(        call SWFRAC(
544       I     imt, minusone,       I     imt, minusone,
545       I     mytime, mythid,       U     worka,
546       U     worka )       I     myTime, myIter, myThid )
547    CADJ store worka = comlev1_kpp
548    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
549    
550        do i = 1, imt        do i = 1, imt
551           bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i))           bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i))
552        end do        end do
553  CADJ store bfsfc = comlev1_kpp  CADJ store bfsfc = comlev1_kpp
554  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
555    
556  c--   ensures bfsfc is never 0  c--   ensures bfsfc is never 0
557        do i = 1, imt        do i = 1, imt
# Line 525  c*************************************** Line 577  c***************************************
577    
578        subroutine wscale (        subroutine wscale (
579       I     sigma, hbl, ustar, bfsfc,       I     sigma, hbl, ustar, bfsfc,
580       O     wm, ws )       O     wm, ws,
581         I     myThid )
582    
583  c     compute turbulent velocity scales.  c     compute turbulent velocity scales.
584  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 546  c sigma   : normalized depth (d/hbl) Line 599  c sigma   : normalized depth (d/hbl)
599  c hbl     : boundary layer depth (m)  c hbl     : boundary layer depth (m)
600  c ustar   : surface friction velocity         (m/s)  c ustar   : surface friction velocity         (m/s)
601  c bfsfc   : total surface buoyancy flux       (m^2/s^3)  c bfsfc   : total surface buoyancy flux       (m^2/s^3)
602        _KPP_RL sigma(imt)  c myThid  : thread number for this instance of the routine
603        _KPP_RL hbl  (imt)        integer myThid
604        _KPP_RL ustar(imt)        _RL sigma(imt)
605        _KPP_RL bfsfc(imt)        _RL hbl  (imt)
606          _RL ustar(imt)
607          _RL bfsfc(imt)
608    
609  c  output  c  output
610  c--------  c--------
611  c wm, ws  : turbulent velocity scales at sigma  c wm, ws  : turbulent velocity scales at sigma
612        _KPP_RL wm(imt), ws(imt)        _RL wm(imt), ws(imt)
613    
614  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
615    
616  c local  c local
617  c------  c------
618  c zehat   : = zeta *  ustar**3  c zehat   : = zeta *  ustar**3
619        _KPP_RL zehat        _RL zehat
620    
621        integer iz, izp1, ju, i, jup1        integer iz, izp1, ju, i, jup1
622        _KPP_RL udiff, zdiff, zfrac, ufrac, fzfrac, wam        _RL udiff, zdiff, zfrac, ufrac, fzfrac, wam
623        _KPP_RL wbm, was, wbs, u3, tempVar        _RL wbm, was, wbs, u3, tempVar
624    
625  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
626  c use lookup table for zehat < zmax only; otherwise use  c use lookup table for zehat < zmax only; otherwise use
# Line 620  c--------------------------------------- Line 675  c---------------------------------------
675  c*************************************************************************  c*************************************************************************
676    
677        subroutine Ri_iwmix (        subroutine Ri_iwmix (
678       I       kmtj, shsq, dbloc, dblocSm       I       kmtj, shsq, dbloc, dblocSm,
679       I     , ikey       I       diffusKzS, diffusKzT,
680       O     , diffus )       I       ikppkey,
681         O       diffus,
682         I       myThid )
683    
684  c     compute interior viscosity diffusivity coefficients due  c     compute interior viscosity diffusivity coefficients due
685  c     to shear instability (dependent on a local Richardson number),  c     to shear instability (dependent on a local Richardson number),
# Line 641  c     kmtj   (imt)         number of ver Line 698  c     kmtj   (imt)         number of ver
698  c     shsq   (imt,Nr)      (local velocity shear)^2               ((m/s)^2)  c     shsq   (imt,Nr)      (local velocity shear)^2               ((m/s)^2)
699  c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)  c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)
700  c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)  c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)
701        integer kmtj   (imt)  c     diffusKzS(imt,Nr)- background vertical diffusivity for scalars    (m^2/s)
702        _KPP_RL shsq   (imt,Nr)  c     diffusKzT(imt,Nr)- background vertical diffusivity for theta      (m^2/s)
703        _KPP_RL dbloc  (imt,Nr)  c     myThid :: My Thread Id. number
704        _KPP_RL dblocSm(imt,Nr)        integer kmtj (imt)
705        integer ikey        _RL shsq     (imt,Nr)
706          _RL dbloc    (imt,Nr)
707          _RL dblocSm  (imt,Nr)
708          _RL diffusKzS(imt,Nr)
709          _RL diffusKzT(imt,Nr)
710          integer ikppkey
711          integer myThid
712    
713  c output  c output
714  c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)  c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)
715  c     diffus(imt,0:Nrp1,2)  vertical scalar diffusivity             (m^2/s)  c     diffus(imt,0:Nrp1,2)  vertical scalar diffusivity             (m^2/s)
716  c     diffus(imt,0:Nrp1,3)  vertical temperature diffusivity        (m^2/s)  c     diffus(imt,0:Nrp1,3)  vertical temperature diffusivity        (m^2/s)
717        _KPP_RL diffus(imt,0:Nrp1,3)        _RL diffus(imt,0:Nrp1,3)
718    
719  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
720    
721  c local variables  c local variables
722  c     Rig                   local Richardson number  c     Rig                   local Richardson number
723  c     fRi, fcon             function of Rig  c     fRi, fcon             function of Rig
724        _KPP_RL Rig        _RL Rig
725        _KPP_RL fRi, fcon        _RL fRi, fcon
726        _KPP_RL ratio        _RL ratio
727        integer i, ki        integer i, ki
728        _KPP_RL c1, c0        _RL c1, c0
729    
730  #ifdef ALLOW_KPP_VERTICALLY_SMOOTH  #ifdef ALLOW_KPP_VERTICALLY_SMOOTH
731        integer mr        integer mr
# Line 681  c     set values at bottom and below to Line 744  c     set values at bottom and below to
744  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
745  C     break data flow dependence on diffus  C     break data flow dependence on diffus
746        diffus(1,1,1) = 0.0        diffus(1,1,1) = 0.0
747    
748          do ki = 1, Nr
749             do i = 1, imt
750                diffus(i,ki,1) = 0.
751                diffus(i,ki,2) = 0.
752                diffus(i,ki,3) = 0.
753             enddo
754          enddo
755  #endif  #endif
756                
757    
758        do ki = 1, Nr        do ki = 1, Nr
759           do i = 1, imt           do i = 1, imt
# Line 698  C     break data flow dependence on diff Line 770  C     break data flow dependence on diff
770              endif              endif
771           end do           end do
772        end do        end do
773    CADJ store diffus = comlev1_kpp, key = ikppkey
774    
775  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
776  c     vertically smooth Ri  c     vertically smooth Ri
# Line 708  CADJ store diffus(:,:,1) = kpp_ri_tape_m Line 781  CADJ store diffus(:,:,1) = kpp_ri_tape_m
781  CADJ &  , key=mr, shape=(/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2 /)  CADJ &  , key=mr, shape=(/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2 /)
782    
783           call z121 (           call z121 (
784       U     diffus(1,0,1))       U     diffus(1,0,1)
785         I     myThid )
786        end do        end do
787  #endif  #endif
788    
# Line 735  c  evaluate f of smooth Ri for shear ins Line 809  c  evaluate f of smooth Ri for shear ins
809  c ----------------------------------------------------------------------  c ----------------------------------------------------------------------
810  c            evaluate diffusivities and viscosity  c            evaluate diffusivities and viscosity
811  c    mixing due to internal waves, and shear and static instability  c    mixing due to internal waves, and shear and static instability
812    
813              diffus(i,ki,1) = viscAr + fcon * difmcon + fRi * difm0  #ifndef EXCLUDE_KPP_SHEAR_MIX
814              diffus(i,ki,2) = diffKrS + fcon * difscon + fRi * difs0              if ( .NOT. inAdMode ) then
815              diffus(i,ki,3) = diffKrT + fcon * difscon + fRi * difs0                 diffus(i,ki,1) = viscAr + fcon * difmcon + fRi * difm0
816                   diffus(i,ki,2) = diffusKzS(i,ki)+fcon*difscon+fRi*difs0
817                   diffus(i,ki,3) = diffusKzT(i,ki)+fcon*diftcon+fRi*dift0
818                else
819                   diffus(i,ki,1) = viscAr
820                   diffus(i,ki,2) = diffusKzS(i,ki)
821                   diffus(i,ki,3) = diffusKzT(i,ki)
822                endif
823    #else
824                diffus(i,ki,1) = viscAr
825                diffus(i,ki,2) = diffusKzS(i,ki)
826                diffus(i,ki,3) = diffusKzT(i,ki)
827    #endif
828    
829           end do           end do
830        end do        end do
# Line 760  c         set surface values to 0.0 Line 846  c         set surface values to 0.0
846  c*************************************************************************  c*************************************************************************
847    
848        subroutine z121 (        subroutine z121 (
849       U     v )       U     v,
850         I     myThid )
851    
852  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)
853  c     top (0) value is used as a dummy  c     top (0) value is used as a dummy
# Line 778  c     penetrative convection. Line 865  c     penetrative convection.
865  c input/output  c input/output
866  c-------------  c-------------
867  c v     : 2-D array to be smoothed in Nrp1 direction  c v     : 2-D array to be smoothed in Nrp1 direction
868        _KPP_RL v(imt,0:Nrp1)  c myThid: thread number for this instance of the routine
869          integer myThid
870          _RL v(imt,0:Nrp1)
871    
872  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
873    
874  c local  c local
875        _KPP_RL zwork, zflag        _RL zwork, zflag
876        _KPP_RL KRi_range(1:Nrp1)        _RL KRi_range(1:Nrp1)
877        integer i, k, km1, kp1        integer i, k, km1, kp1
878    
879        _KPP_RL     p0      , p25       , p5      , p2        _RL         p0      , p25       , p5      , p2
880        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 )
881    
882        KRi_range(Nrp1) = p0        KRi_range(Nrp1) = p0
# Line 841  CADJ STORE v(i,k), zwork = z121tape Line 930  CADJ STORE v(i,k), zwork = z121tape
930    
931  c*************************************************************************  c*************************************************************************
932    
       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*************************************************************************  
   
933        subroutine smooth_horiz (        subroutine smooth_horiz (
934       I     k, bi, bj,       I     k, bi, bj,
935       U     fld )       U     fld,
936         I     myThid )
937    
938  c     Apply horizontal smoothing to global _RL 2-D array  c     Apply horizontal smoothing to global _RL 2-D array
939    
940        IMPLICIT NONE        IMPLICIT NONE
941  #include "SIZE.h"  #include "SIZE.h"
942    #include "GRID.h"
943  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
944    
945  c     input  c     input
946  c     bi, bj : array indices  c     bi, bj : array indices
947  c     k      : vertical index used for masking  c     k      : vertical index used for masking
948    c     myThid : thread number for this instance of the routine
949          INTEGER myThid
950        integer k, bi, bj        integer k, bi, bj
951    
952  c     input/output  c     input/output
# Line 960  c     local Line 973  c     local
973              im1 = i-1              im1 = i-1
974              ip1 = i+1              ip1 = i+1
975              tempVar =              tempVar =
976       &           p25   *   pMask(i  ,j  ,k,bi,bj)   +       &           p25   *   maskC(i  ,j  ,k,bi,bj)   +
977       &           p125  * ( pMask(im1,j  ,k,bi,bj)   +       &           p125  * ( maskC(im1,j  ,k,bi,bj)   +
978       &                     pMask(ip1,j  ,k,bi,bj)   +       &                     maskC(ip1,j  ,k,bi,bj)   +
979       &                     pMask(i  ,jm1,k,bi,bj)   +       &                     maskC(i  ,jm1,k,bi,bj)   +
980       &                     pMask(i  ,jp1,k,bi,bj) ) +       &                     maskC(i  ,jp1,k,bi,bj) ) +
981       &           p0625 * ( pMask(im1,jm1,k,bi,bj)   +       &           p0625 * ( maskC(im1,jm1,k,bi,bj)   +
982       &                     pMask(im1,jp1,k,bi,bj)   +       &                     maskC(im1,jp1,k,bi,bj)   +
983       &                     pMask(ip1,jm1,k,bi,bj)   +       &                     maskC(ip1,jm1,k,bi,bj)   +
984       &                     pMask(ip1,jp1,k,bi,bj) )       &                     maskC(ip1,jp1,k,bi,bj) )
985              IF ( tempVar .GE. p25 ) THEN              IF ( tempVar .GE. p25 ) THEN
986                 fld_tmp(i,j) = (                 fld_tmp(i,j) = (
987       &              p25  * fld(i  ,j  )*pMask(i  ,j  ,k,bi,bj) +       &              p25  * fld(i  ,j  )*maskC(i  ,j  ,k,bi,bj) +
988       &              p125 *(fld(im1,j  )*pMask(im1,j  ,k,bi,bj) +       &              p125 *(fld(im1,j  )*maskC(im1,j  ,k,bi,bj) +
989       &                     fld(ip1,j  )*pMask(ip1,j  ,k,bi,bj) +       &                     fld(ip1,j  )*maskC(ip1,j  ,k,bi,bj) +
990       &                     fld(i  ,jm1)*pMask(i  ,jm1,k,bi,bj) +       &                     fld(i  ,jm1)*maskC(i  ,jm1,k,bi,bj) +
991       &                     fld(i  ,jp1)*pMask(i  ,jp1,k,bi,bj))+       &                     fld(i  ,jp1)*maskC(i  ,jp1,k,bi,bj))+
992       &              p0625*(fld(im1,jm1)*pMask(im1,jm1,k,bi,bj) +       &              p0625*(fld(im1,jm1)*maskC(im1,jm1,k,bi,bj) +
993       &                     fld(im1,jp1)*pMask(im1,jp1,k,bi,bj) +       &                     fld(im1,jp1)*maskC(im1,jp1,k,bi,bj) +
994       &                     fld(ip1,jm1)*pMask(ip1,jm1,k,bi,bj) +       &                     fld(ip1,jm1)*maskC(ip1,jm1,k,bi,bj) +
995       &                     fld(ip1,jp1)*pMask(ip1,jp1,k,bi,bj)))       &                     fld(ip1,jp1)*maskC(ip1,jp1,k,bi,bj)))
996       &              / tempVar       &              / tempVar
997              ELSE              ELSE
998                 fld_tmp(i,j) = fld(i,j)                 fld_tmp(i,j) = fld(i,j)
# Line 1003  c*************************************** Line 1016  c***************************************
1016    
1017        subroutine blmix (        subroutine blmix (
1018       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl
1019       O     , dkm1, blmc, ghat, sigma, ikey       O     , dkm1, blmc, ghat, sigma, ikppkey
1020       &     )       I     , myThid )
1021    
1022  c     mixing coefficients within boundary layer depend on surface  c     mixing coefficients within boundary layer depend on surface
1023  c     forcing and the magnitude and gradient of interior mixing below  c     forcing and the magnitude and gradient of interior mixing below
# Line 1026  c     hbl   (imt)                 bounda Line 1039  c     hbl   (imt)                 bounda
1039  c     stable(imt)                 = 1 in stable forcing  c     stable(imt)                 = 1 in stable forcing
1040  c     casea (imt)                 = 1 in case A  c     casea (imt)                 = 1 in case A
1041  c     diffus(imt,0:Nrp1,mdiff)    vertical diffusivities              (m^2/s)  c     diffus(imt,0:Nrp1,mdiff)    vertical diffusivities              (m^2/s)
1042  c     kbl(imt)                    -1 of first grid level below hbl  c     kbl   (imt)                 -1 of first grid level below hbl
1043        _KPP_RL ustar (imt)  c     myThid               thread number for this instance of the routine
1044        _KPP_RL bfsfc (imt)        integer myThid
1045        _KPP_RL hbl   (imt)        _RL ustar (imt)
1046        _KPP_RL stable(imt)        _RL bfsfc (imt)
1047        _KPP_RL casea (imt)        _RL hbl   (imt)
1048        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL stable(imt)
1049          _RL casea (imt)
1050          _RL diffus(imt,0:Nrp1,mdiff)
1051        integer kbl(imt)        integer kbl(imt)
1052    
1053  c output  c output
# Line 1040  c     dkm1 (imt,mdiff)            bounda Line 1055  c     dkm1 (imt,mdiff)            bounda
1055  c     blmc (imt,Nr,mdiff)         boundary layer mixing coefficients  (m^2/s)  c     blmc (imt,Nr,mdiff)         boundary layer mixing coefficients  (m^2/s)
1056  c     ghat (imt,Nr)               nonlocal scalar transport  c     ghat (imt,Nr)               nonlocal scalar transport
1057  c     sigma(imt)                  normalized depth (d / hbl)  c     sigma(imt)                  normalized depth (d / hbl)
1058        _KPP_RL dkm1 (imt,mdiff)        _RL dkm1 (imt,mdiff)
1059        _KPP_RL blmc (imt,Nr,mdiff)        _RL blmc (imt,Nr,mdiff)
1060        _KPP_RL ghat (imt,Nr)        _RL ghat (imt,Nr)
1061        _KPP_RL sigma(imt)        _RL sigma(imt)
1062        integer ikey        integer ikppkey, kkppkey
1063    
1064  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1065    
# Line 1052  c  local Line 1067  c  local
1067  c     gat1*(imt)                 shape function at sigma = 1  c     gat1*(imt)                 shape function at sigma = 1
1068  c     dat1*(imt)                 derivative of shape function at sigma = 1  c     dat1*(imt)                 derivative of shape function at sigma = 1
1069  c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)  c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)
1070        _KPP_RL gat1m(imt), gat1s(imt), gat1t(imt)        _RL gat1m(imt), gat1s(imt), gat1t(imt)
1071        _KPP_RL dat1m(imt), dat1s(imt), dat1t(imt)        _RL dat1m(imt), dat1s(imt), dat1t(imt)
1072        _KPP_RL ws(imt), wm(imt)        _RL ws(imt), wm(imt)
1073        integer i, kn, ki        integer i, kn, ki
1074        _KPP_RL R, dvdzup, dvdzdn, viscp        _RL R, dvdzup, dvdzdn, viscp
1075        _KPP_RL difsp, diftp, visch, difsh, difth        _RL difsp, diftp, visch, difsh, difth
1076        _KPP_RL f1, sig, a1, a2, a3, delhat        _RL f1, sig, a1, a2, a3, delhat
1077        _KPP_RL Gm, Gs, Gt        _RL Gm, Gs, Gt
1078        _KPP_RL tempVar        _RL tempVar
1079    
1080        _KPP_RL    p0    , eins        _RL    p0    , eins
1081        parameter (p0=0.0, eins=1.0)        parameter (p0=0.0, eins=1.0)
1082    
1083  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 1073  c--------------------------------------- Line 1088  c---------------------------------------
1088           sigma(i) = stable(i) * 1.0 + (1. - stable(i)) * epsilon           sigma(i) = stable(i) * 1.0 + (1. - stable(i)) * epsilon
1089        end do        end do
1090    
1091    CADJ STORE sigma = comlev1_kpp, key = ikppkey
1092        call wscale (        call wscale (
1093       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1094       O        wm, ws )       O        wm, ws, myThid )
1095    CADJ STORE wm = comlev1_kpp, key = ikppkey
1096    CADJ STORE ws = comlev1_kpp, key = ikppkey
1097    
1098        do i = 1, imt        do i = 1, imt
1099           wm(i) = sign(eins,wm(i))*max(phepsi,abs(wm(i)))           wm(i) = sign(eins,wm(i))*max(phepsi,abs(wm(i)))
1100           ws(i) = sign(eins,ws(i))*max(phepsi,abs(ws(i)))           ws(i) = sign(eins,ws(i))*max(phepsi,abs(ws(i)))
1101        end do        end do
1102  CADJ STORE wm = comlev1_kpp, key = ikey  CADJ STORE wm = comlev1_kpp, key = ikppkey
1103  CADJ STORE ws = comlev1_kpp, key = ikey  CADJ STORE ws = comlev1_kpp, key = ikppkey
1104    
1105        do i = 1, imt        do i = 1, imt
1106    
# Line 1118  c--------------------------------------- Line 1136  c---------------------------------------
1136       &        max(ustar(i)**4,phepsi)       &        max(ustar(i)**4,phepsi)
1137           gat1m(i) = visch / hbl(i) / wm(i)           gat1m(i) = visch / hbl(i) / wm(i)
1138           dat1m(i) = -viscp / wm(i) + f1 * visch           dat1m(i) = -viscp / wm(i) + f1 * visch
          dat1m(i) = min(dat1m(i),p0)  
1139    
1140           gat1s(i) = difsh  / hbl(i) / ws(i)           gat1s(i) = difsh  / hbl(i) / ws(i)
1141           dat1s(i) = -difsp / ws(i) + f1 * difsh           dat1s(i) = -difsp / ws(i) + f1 * difsh
          dat1s(i) = min(dat1s(i),p0)  
1142    
1143           gat1t(i) = difth /  hbl(i) / ws(i)           gat1t(i) = difth /  hbl(i) / ws(i)
1144           dat1t(i) = -diftp / ws(i) + f1 * difth           dat1t(i) = -diftp / ws(i) + f1 * difth
          dat1t(i) = min(dat1t(i),p0)  
1145    
1146        end do        end do
1147    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1148    CADJ STORE gat1m = comlev1_kpp, key = ikppkey
1149    CADJ STORE gat1s = comlev1_kpp, key = ikppkey
1150    CADJ STORE gat1t = comlev1_kpp, key = ikppkey
1151    CADJ STORE dat1m = comlev1_kpp, key = ikppkey
1152    CADJ STORE dat1s = comlev1_kpp, key = ikppkey
1153    CADJ STORE dat1t = comlev1_kpp, key = ikppkey
1154    #endif
1155          do i = 1, imt
1156             dat1m(i) = min(dat1m(i),p0)
1157             dat1s(i) = min(dat1s(i),p0)
1158             dat1t(i) = min(dat1t(i),p0)
1159          end do
1160    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1161    CADJ STORE dat1m = comlev1_kpp, key = ikppkey
1162    CADJ STORE dat1s = comlev1_kpp, key = ikppkey
1163    CADJ STORE dat1t = comlev1_kpp, key = ikppkey
1164    #endif
1165    
1166        do ki = 1, Nr        do ki = 1, Nr
1167    
1168    #ifdef ALLOW_AUTODIFF_TAMC
1169             kkppkey = (ikppkey-1)*Nr + ki
1170    #endif
1171    
1172  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1173  c     compute turbulent velocity scales on the interfaces  c     compute turbulent velocity scales on the interfaces
1174  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 1140  c--------------------------------------- Line 1177  c---------------------------------------
1177              sig      = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i)              sig      = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i)
1178              sigma(i) = stable(i)*sig + (1.-stable(i))*min(sig,epsilon)              sigma(i) = stable(i)*sig + (1.-stable(i))*min(sig,epsilon)
1179           end do           end do
1180    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1181    CADJ STORE wm = comlev1_kpp_k, key = kkppkey
1182    CADJ STORE ws = comlev1_kpp_k, key = kkppkey
1183    #endif
1184    CADJ STORE sigma = comlev1_kpp_k, key = kkppkey
1185           call wscale (           call wscale (
1186       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1187       O        wm, ws )       O        wm, ws, myThid )
1188    CADJ STORE wm = comlev1_kpp_k, key = kkppkey
1189    CADJ STORE ws = comlev1_kpp_k, key = kkppkey
1190    
1191  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1192  c     compute the dimensionless shape functions at the interfaces  c     compute the dimensionless shape functions at the interfaces
# Line 1186  c--------------------------------------- Line 1230  c---------------------------------------
1230       &            + (1. - stable(i)) * min(sig,epsilon)       &            + (1. - stable(i)) * min(sig,epsilon)
1231        end do        end do
1232    
1233    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1234    CADJ STORE wm = comlev1_kpp, key = ikppkey
1235    CADJ STORE ws = comlev1_kpp, key = ikppkey
1236    #endif
1237    CADJ STORE sigma = comlev1_kpp, key = ikppkey
1238        call wscale (        call wscale (
1239       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1240       O        wm, ws )       O        wm, ws, myThid )
1241    CADJ STORE wm = comlev1_kpp, key = ikppkey
1242    CADJ STORE ws = comlev1_kpp, key = ikppkey
1243    
1244        do i = 1, imt        do i = 1, imt
1245           sig = -zgrid(kbl(i)-1) / hbl(i)           sig = -zgrid(kbl(i)-1) / hbl(i)
# Line 1214  c*************************************** Line 1265  c***************************************
1265       I       dkm1, hbl, kbl, diffus, casea       I       dkm1, hbl, kbl, diffus, casea
1266       U     , ghat       U     , ghat
1267       O     , blmc       O     , blmc
1268       &     )       &     , myThid )
1269    
1270  c enhance the diffusivity at the kbl-.5 interface  c enhance the diffusivity at the kbl-.5 interface
1271    
# Line 1229  c     hbl(imt)                  boundary Line 1280  c     hbl(imt)                  boundary
1280  c     kbl(imt)                  grid above hbl  c     kbl(imt)                  grid above hbl
1281  c     diffus(imt,0:Nrp1,mdiff) vertical diffusivities           (m^2/s)  c     diffus(imt,0:Nrp1,mdiff) vertical diffusivities           (m^2/s)
1282  c     casea(imt)                = 1 in caseA, = 0 in case B  c     casea(imt)                = 1 in caseA, = 0 in case B
1283        _KPP_RL dkm1  (imt,mdiff)  c     myThid                    thread number for this instance of the routine
1284        _KPP_RL hbl   (imt)        integer   myThid
1285          _RL dkm1  (imt,mdiff)
1286          _RL hbl   (imt)
1287        integer kbl   (imt)        integer kbl   (imt)
1288        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL diffus(imt,0:Nrp1,mdiff)
1289        _KPP_RL casea (imt)        _RL casea (imt)
1290    
1291  c input/output  c input/output
1292  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)
1293        _KPP_RL ghat (imt,Nr)        _RL ghat (imt,Nr)
1294    
1295  c output  c output
1296  c     enhanced bound. layer mixing coeff.  c     enhanced bound. layer mixing coeff.
1297        _KPP_RL blmc  (imt,Nr,mdiff)        _RL blmc  (imt,Nr,mdiff)
1298    
1299  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1300    
1301  c local  c local
1302  c     fraction hbl lies beteen zgrid neighbors  c     fraction hbl lies beteen zgrid neighbors
1303        _KPP_RL delta        _RL delta
1304        integer ki, i, md        integer ki, i, md
1305        _KPP_RL dkmp5, dstar        _RL dkmp5, dstar
1306    
1307        do i = 1, imt        do i = 1, imt
1308           ki = kbl(i)-1           ki = kbl(i)-1
# Line 1275  c     fraction hbl lies beteen zgrid nei Line 1328  c     fraction hbl lies beteen zgrid nei
1328  c*************************************************************************  c*************************************************************************
1329    
1330        SUBROUTINE STATEKPP (        SUBROUTINE STATEKPP (
1331       I     bi, bj, myThid,       O     RHO1, DBLOC, DBSFC, TTALPHA, SSBETA,
1332       O     RHO1, DBLOC, DBSFC, TTALPHA, SSBETA)       I     ikppkey, bi, bj, myThid )
1333  c  c
1334  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1335  c     "statekpp" computes all necessary input arrays  c     "statekpp" computes all necessary input arrays
# Line 1301  c Line 1354  c
1354  c     written  by: jan morzel,   feb. 10, 1995 (converted from "sigma" version)  c     written  by: jan morzel,   feb. 10, 1995 (converted from "sigma" version)
1355  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
1356  c  c
1357    c     28 april 05: added computation of mixed layer depth KPPmld
1358    c     for a deltaRho equivalent to a deltaTheta of 0.8 deg C
1359    
1360  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1361    
1362        IMPLICIT NONE        IMPLICIT NONE
# Line 1310  c--------------------------------------- Line 1366  c---------------------------------------
1366  #include "PARAMS.h"  #include "PARAMS.h"
1367  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
1368  #include "DYNVARS.h"  #include "DYNVARS.h"
1369    #include "GRID.h"
1370    
1371  c-------------- Routine arguments -----------------------------------------  c-------------- Routine arguments -----------------------------------------
1372        INTEGER bi, bj, myThid        INTEGER bi, bj, myThid
1373        _KPP_RL RHO1   ( ibot:itop, jbot:jtop       )        _RL RHO1   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy       )
1374        _KPP_RL DBLOC  ( ibot:itop, jbot:jtop, Nr   )        _RL DBLOC  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
1375        _KPP_RL DBSFC  ( ibot:itop, jbot:jtop, Nr   )        _RL DBSFC  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
1376        _KPP_RL TTALPHA( ibot:itop, jbot:jtop, Nrp1 )        _RL TTALPHA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1377        _KPP_RL SSBETA ( ibot:itop, jbot:jtop, Nrp1 )        _RL SSBETA ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1378    
1379  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1380    
# Line 1328  c Line 1385  c
1385  c     rhok         - density of t(k  ) & s(k  ) at depth k  c     rhok         - density of t(k  ) & s(k  ) at depth k
1386  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
1387  c     rho1k        - density of t(1  ) & s(1  ) at depth k  c     rho1k        - density of t(1  ) & s(1  ) at depth k
1388  c     work1, work2 - work arrays for holding horizontal slabs  c     work1,2,3    - work arrays for holding horizontal slabs
1389    
1390        _RL RHOK  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL RHOK  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1391        _RL RHOKM1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL RHOKM1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 1336  c     work1, work2 - work arrays for hol Line 1393  c     work1, work2 - work arrays for hol
1393        _RL WORK1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1394        _RL WORK2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1395        _RL WORK3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1396    #ifdef ALLOW_DIAGNOSTICS
1397    c     KPPMLD - mixed layer depth based on density criterion
1398          _RL KPPMLD(1    :sNx    ,1    :sNy    )
1399    #endif /* ALLOW_DIAGNOSTICS */
1400    
1401        INTEGER I, J, K        INTEGER I, J, K
1402          INTEGER ikppkey, kkppkey
1403    
1404  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
1405    
1406          kkppkey = (ikppkey-1)*Nr + 1
1407    
1408    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1409    CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1410    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1411    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1412        call FIND_RHO(        call FIND_RHO(
1413       I     bi, bj, ibot, itop, jbot, jtop, 1, 1,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1414       I     theta, salt,       I     theta, salt,
1415       O     WORK1,       O     WORK1,
1416       I     myThid )       I     myThid )
1417    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1418    CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1419    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1420    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1421    
1422        call FIND_ALPHA(        call FIND_ALPHA(
1423       I     bi, bj, ibot, itop, jbot, jtop, 1, 1,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1424       O     WORK2 )       O     WORK2, myThid )
1425    
1426        call FIND_BETA(        call FIND_BETA(
1427       I     bi, bj, ibot, itop, jbot, jtop, 1, 1,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1428       O     WORK3 )       O     WORK3, myThid )
1429    
1430        DO J = jbot, jtop        DO J = 1-OLy, sNy+OLy
1431           DO I = ibot, itop           DO I = 1-OLx, sNx+OLx
1432              RHO1(I,J)      = WORK1(I,J) + rhoConst              RHO1(I,J)      = WORK1(I,J) + rhoConst
1433              TTALPHA(I,J,1) = WORK2(I,J)              TTALPHA(I,J,1) = WORK2(I,J)
1434              SSBETA(I,J,1)  = WORK3(I,J)              SSBETA(I,J,1)  = WORK3(I,J)
# Line 1363  c calculate density, alpha, beta in surf Line 1436  c calculate density, alpha, beta in surf
1436           END DO           END DO
1437        END DO        END DO
1438    
1439    #ifdef ALLOW_DIAGNOSTICS
1440    c     work3 - density of t(1)-.8 & s(1  ) at depth 1
1441          IF ( useDiagnostics ) THEN
1442             DO J = 1, sNy
1443                DO I = 1, sNx
1444                   KPPMLD(I,J) = ABS(R_low(I,J,bi,bj))
1445                   WORK3 (I,J) = WORK1(I,J) - 0.8 * WORK2(I,J)
1446                END DO
1447             END DO
1448          ENDIF
1449    #endif /* ALLOW_DIAGNOSTICS */
1450    
1451  c calculate alpha, beta, and gradients in interior layers  c calculate alpha, beta, and gradients in interior layers
1452    
1453  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)
1454        DO K = 2, Nr        DO K = 2, Nr
1455    
1456          kkppkey = (ikppkey-1)*Nr + k
1457    
1458    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1459    CADJ STORE theta(:,:,k,bi,bj) = comlev1_kpp_k, key=kkppkey
1460    CADJ STORE salt (:,:,k,bi,bj) = comlev1_kpp_k, key=kkppkey
1461    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1462           call FIND_RHO(           call FIND_RHO(
1463       I        bi, bj, ibot, itop, jbot, jtop, K, K,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1464       I        theta, salt,       I        theta, salt,
1465       O        RHOK,       O        RHOK,
1466       I        myThid )       I        myThid )
1467    
1468    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1469    CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_kpp_k, key=kkppkey
1470    CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_kpp_k, key=kkppkey
1471    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1472           call FIND_RHO(           call FIND_RHO(
1473       I        bi, bj, ibot, itop, jbot, jtop, K-1, K,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K-1, K,
1474       I        theta, salt,       I        theta, salt,
1475       O        RHOKM1,       O        RHOKM1,
1476       I        myThid )       I        myThid )
1477    
1478    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1479    CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1480    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1481    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1482           call FIND_RHO(           call FIND_RHO(
1483       I        bi, bj, ibot, itop, jbot, jtop, 1, K,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, K,
1484       I        theta, salt,       I        theta, salt,
1485       O        RHO1K,       O        RHO1K,
1486       I        myThid )       I        myThid )
1487    
1488    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1489    CADJ STORE rhok  (:,:)          = comlev1_kpp_k, key=kkppkey
1490    CADJ STORE rhokm1(:,:)          = comlev1_kpp_k, key=kkppkey
1491    CADJ STORE rho1k (:,:)          = comlev1_kpp_k, key=kkppkey
1492    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1493    
1494           call FIND_ALPHA(           call FIND_ALPHA(
1495       I        bi, bj, ibot, itop, jbot, jtop, K, K,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1496       O        WORK1 )       O        WORK1, myThid )
1497    
1498           call FIND_BETA(           call FIND_BETA(
1499       I        bi, bj, ibot, itop, jbot, jtop, K, K,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1500       O        WORK2 )       O        WORK2, myThid )
1501    
1502           DO J = jbot, jtop           DO J = 1-OLy, sNy+OLy
1503              DO I = ibot, itop              DO I = 1-OLx, sNx+OLx
1504                 TTALPHA(I,J,K) = WORK1 (I,J)                 TTALPHA(I,J,K) = WORK1 (I,J)
1505                 SSBETA(I,J,K)  = WORK2 (I,J)                 SSBETA(I,J,K)  = WORK2 (I,J)
1506                 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 1405  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO Line 1510  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO
1510              END DO              END DO
1511           END DO           END DO
1512    
1513    #ifdef ALLOW_DIAGNOSTICS
1514             IF ( useDiagnostics ) THEN
1515    c     work1 - density of t(k-1)  & s(k-1) at depth 1
1516    c     work2 - density of t(k  )  & s(k  ) at depth 1
1517    c     work3 - density of t(1)-.8 & s(1  ) at depth 1
1518                call FIND_RHO(
1519         I           bi, bj, 1, sNx, 1, sNy, K-1, 1, theta, salt,
1520         O           WORK1,
1521         I           myThid )
1522                call FIND_RHO(
1523         I           bi, bj, 1, sNx, 1, sNy, K  , 1, theta, salt,
1524         O           WORK2,
1525         I           myThid )
1526                DO J = 1, sNy
1527                   DO I = 1, sNx
1528                      IF ( k .LE. klowC(I,J,bi,bj) .AND.
1529         &                 WORK1(I,J) .LT. WORK3(I,J) .AND.
1530         &                 WORK2(I,J) .GE. WORK3(I,J) ) THEN
1531                         KPPMLD(I,J) = MIN ( KPPMLD(I,J),
1532         &                    ABS(((WORK3(I,J)-WORK1(I,J))*rC(k)+
1533         &                    (WORK2(I,J)-WORK3(I,J))*rC(k-1))/
1534         &                    (WORK2(I,J)-WORK1(I,J))))
1535                      ENDIF
1536                   END DO
1537                END DO
1538             ENDIF
1539    #endif /* ALLOW_DIAGNOSTICS */
1540    
1541        END DO        END DO
1542    
1543  c     compute arrays for K = Nrp1  c     compute arrays for K = Nrp1
1544        DO J = jbot, jtop        DO J = 1-OLy, sNy+OLy
1545           DO I = ibot, itop           DO I = 1-OLx, sNx+OLx
1546              TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)              TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)
1547              SSBETA(I,J,Nrp1)  = SSBETA(I,J,Nr)              SSBETA(I,J,Nrp1)  = SSBETA(I,J,Nr)
1548              DBLOC(I,J,Nr)     = 0.              DBLOC(I,J,Nr)     = 0.
1549           END DO           END DO
1550        END DO        END DO
1551    
1552    #ifdef ALLOW_DIAGNOSTICS
1553          IF ( useDiagnostics ) THEN
1554             CALL DIAGNOSTICS_FILL(KPPmld,'KPPmld  ',0,1,3,bi,bj,myThid)
1555          ENDIF
1556    #endif /* ALLOW_DIAGNOSTICS */
1557    
1558  #endif /* ALLOW_KPP */  #endif /* ALLOW_KPP */
1559    
1560        RETURN        RETURN
Legend:
Removed from v.1.13  
changed lines
  Added in v.1.31
  ViewVC Help
Powered by ViewVC 1.1.22