/[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.7 by cnh, Sun Feb 4 14:38:50 2001 UTC revision 1.29 by jmc, Mon Apr 30 13:49:40 2007 UTC
# Line 11  C--  o BLDEPTH     - Determine oceanic p Line 11  C--  o BLDEPTH     - Determine oceanic p
11  C--  o WSCALE      - Compute turbulent velocity scales.  C--  o WSCALE      - Compute turbulent velocity scales.
12  C--  o RI_IWMIX    - Compute interior viscosity diffusivity coefficients.  C--  o RI_IWMIX    - Compute interior viscosity diffusivity coefficients.
13  C--  o Z121        - Apply 121 vertical smoothing.  C--  o Z121        - Apply 121 vertical smoothing.
14  C--  o KPP_SMOOTH_HORIZ - Apply horizontal smoothing to KPP array.  C--  o SMOOTH_HORIZ- Apply horizontal smoothing to global array.
 C--  o SMOOTH_HORIZ     - Apply horizontal smoothing to global array.  
15  C--  o BLMIX       - Boundary layer mixing coefficients.  C--  o BLMIX       - Boundary layer mixing coefficients.
16  C--  o ENHANCE     - Enhance diffusivity at boundary layer interface.  C--  o ENHANCE     - Enhance diffusivity at boundary layer interface.
17  C--  o STATEKPP    - Compute buoyancy-related input arrays.  C--  o STATEKPP    - Compute buoyancy-related input arrays.
# Line 23  c*************************************** Line 22  c***************************************
22       I       mytime, mythid       I       mytime, mythid
23       I     , kmtj, shsq, dvsq, ustar       I     , kmtj, shsq, dvsq, ustar
24       I     , bo, bosol, dbloc, Ritop, coriol       I     , bo, bosol, dbloc, Ritop, coriol
25       I     , ikey       I     , diffusKzS, diffusKzT
26         I     , ikppkey
27       O     , diffus       O     , diffus
28       U     , ghat       U     , ghat
29       O     , hbl )       O     , hbl )
# Line 47  c--------------------------------------- Line 47  c---------------------------------------
47  #include "SIZE.h"  #include "SIZE.h"
48  #include "EEPARAMS.h"  #include "EEPARAMS.h"
49  #include "PARAMS.h"  #include "PARAMS.h"
 #include "DYNVARS.h"  
 #include "FFIELDS.h"  
50  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
51    
52  c input  c input
53  c     myTime          - current time in simulation  c     myTime           - current time in simulation
54  c     myThid          - thread number for this instance of the routine  c     myThid           - thread number for this instance of the routine
55  c     kmtj   (imt)    - number of vertical layers on this row  c     kmtj   (imt)     - number of vertical layers on this row
56  c     shsq   (imt,Nr) - (local velocity shear)^2                      ((m/s)^2)  c     shsq   (imt,Nr)  - (local velocity shear)^2                     ((m/s)^2)
57  c     dvsq   (imt,Nr) - (velocity shear re sfc)^2                     ((m/s)^2)  c     dvsq   (imt,Nr)  - (velocity shear re sfc)^2                    ((m/s)^2)
58  c     ustar  (imt)    - surface friction velocity                         (m/s)  c     ustar  (imt)     - surface friction velocity                        (m/s)
59  c     bo     (imt)    - surface turbulent buoy. forcing               (m^2/s^3)  c     bo     (imt)     - surface turbulent buoy. forcing              (m^2/s^3)
60  c     bosol  (imt)    - radiative buoyancy forcing                    (m^2/s^3)  c     bosol  (imt)     - radiative buoyancy forcing                   (m^2/s^3)
61  c     dbloc  (imt,Nr) - local delta buoyancy across interfaces          (m/s^2)  c     dbloc  (imt,Nr)  - local delta buoyancy across interfaces         (m/s^2)
62  c     dblocSm(imt,Nr) - horizontally smoothed dbloc                     (m/s^2)  c     dblocSm(imt,Nr)  - horizontally smoothed dbloc                    (m/s^2)
63  c                         stored in ghat to save space  c                          stored in ghat to save space
64  c     Ritop  (imt,Nr) - numerator of bulk Richardson Number  c     Ritop  (imt,Nr)  - numerator of bulk Richardson Number
65  c                         (zref-z) * delta buoyancy w.r.t. surface    ((m/s)^2)  c                          (zref-z) * delta buoyancy w.r.t. surface   ((m/s)^2)
66  c     coriol (imt)    - Coriolis parameter                                (1/s)  c     coriol (imt)     - Coriolis parameter                               (1/s)
67    c     diffusKzS(imt,Nr)- background vertical diffusivity for scalars    (m^2/s)
68    c     diffusKzT(imt,Nr)- background vertical diffusivity for theta      (m^2/s)
69  c     note: there is a conversion from 2-D to 1-D for input output variables,  c     note: there is a conversion from 2-D to 1-D for input output variables,
70  c           e.g., hbl(sNx,sNy) -> hbl(imt),  c           e.g., hbl(sNx,sNy) -> hbl(imt),
71  c           where hbl(i,j) -> hbl((j-1)*sNx+i)  c           where hbl(i,j) -> hbl((j-1)*sNx+i)
72    
73        _RL     mytime        _RL     mytime
74        integer mythid        integer mythid
75        integer kmtj  (imt   )        integer kmtj (imt   )
76        _KPP_RL shsq  (imt,Nr)        _RL shsq     (imt,Nr)
77        _KPP_RL dvsq  (imt,Nr)        _RL dvsq     (imt,Nr)
78        _KPP_RL ustar (imt   )        _RL ustar    (imt   )
79        _KPP_RL bo    (imt   )        _RL bo       (imt   )
80        _KPP_RL bosol (imt   )        _RL bosol    (imt   )
81        _KPP_RL dbloc (imt,Nr)        _RL dbloc    (imt,Nr)
82        _KPP_RL Ritop (imt,Nr)        _RL Ritop    (imt,Nr)
83        _KPP_RL coriol(imt   )        _RL coriol   (imt   )
84          _RL diffusKzS(imt,Nr)
85          _RL diffusKzT(imt,Nr)
86    
87        integer ikey        integer ikppkey
88    
89  c output  c output
90  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 93  c     diffus (imt,3)  - vertical tempera
93  c     ghat   (imt)    - nonlocal transport coefficient                  (s/m^2)  c     ghat   (imt)    - nonlocal transport coefficient                  (s/m^2)
94  c     hbl    (imt)    - mixing layer depth                                  (m)  c     hbl    (imt)    - mixing layer depth                                  (m)
95    
96        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL diffus(imt,0:Nrp1,mdiff)
97        _KPP_RL ghat  (imt,Nr)        _RL ghat  (imt,Nr)
98        _KPP_RL hbl   (imt)        _RL hbl   (imt)
99    
100  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
101    
# Line 107  c     blmc   (imt,Nr,mdiff) - boundary l Line 109  c     blmc   (imt,Nr,mdiff) - boundary l
109  c     sigma  (imt         ) - normalized depth (d / hbl)  c     sigma  (imt         ) - normalized depth (d / hbl)
110  c     Rib    (imt,Nr      ) - bulk Richardson number  c     Rib    (imt,Nr      ) - bulk Richardson number
111    
112        integer kbl   (imt         )        integer kbl(imt         )
113        _KPP_RL bfsfc (imt         )        _RL bfsfc  (imt         )
114        _KPP_RL casea (imt         )        _RL casea  (imt         )
115        _KPP_RL stable(imt         )        _RL stable (imt         )
116        _KPP_RL dkm1  (imt,   mdiff)        _RL dkm1   (imt,   mdiff)
117        _KPP_RL blmc  (imt,Nr,mdiff)        _RL blmc   (imt,Nr,mdiff)
118        _KPP_RL sigma (imt         )        _RL sigma  (imt         )
119        _KPP_RL Rib   (imt,Nr      )        _RL Rib    (imt,Nr      )
120    
121        integer i, k, md        integer i, k, md
122    
# Line 125  c instability. Line 127  c instability.
127  c (ghat is temporary storage for horizontally smoothed dbloc)  c (ghat is temporary storage for horizontally smoothed dbloc)
128  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
129    
130  CADJ STORE ghat = comlev1_kpp, key = ikey  cph(
131    cph these storings avoid recomp. of Ri_iwmix
132    CADJ STORE ghat  = comlev1_kpp, key = ikppkey
133    CADJ STORE dbloc = comlev1_kpp, key = ikppkey
134    cph)
135        call Ri_iwmix (        call Ri_iwmix (
136       I       kmtj, shsq, dbloc, ghat       I       kmtj, shsq, dbloc, ghat
137       I     , ikey       I     , diffusKzS, diffusKzT
138         I     , ikppkey
139       O     , diffus )       O     , diffus )
140    
141    cph(
142    cph these storings avoid recomp. of Ri_iwmix
143    cph DESPITE TAFs 'not necessary' warning!
144    CADJ STORE dbloc  = comlev1_kpp, key = ikppkey
145    CADJ STORE shsq   = comlev1_kpp, key = ikppkey
146    CADJ STORE ghat   = comlev1_kpp, key = ikppkey
147    CADJ STORE diffus = comlev1_kpp, key = ikppkey
148    cph)
149    
150  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
151  c set seafloor values to zero and fill extra "Nrp1" coefficients  c set seafloor values to zero and fill extra "Nrp1" coefficients
152  c for blmix  c for blmix
# Line 155  c--------------------------------------- Line 170  c---------------------------------------
170       I       mytime, mythid       I       mytime, mythid
171       I     , kmtj       I     , kmtj
172       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol
173       I     , ikey       I     , ikppkey
174       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma
175       &     )       &     )
176    
177  CADJ STORE hbl,bfsfc,stable,casea,kbl = comlev1_kpp, key = ikey  CADJ STORE hbl,bfsfc,stable,casea,kbl = comlev1_kpp, key = ikppkey
178    
179  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
180  c compute boundary layer diffusivities  c compute boundary layer diffusivities
# Line 167  c--------------------------------------- Line 182  c---------------------------------------
182    
183        call blmix (        call blmix (
184       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl
185       O     , dkm1, blmc, ghat, sigma, ikey       O     , dkm1, blmc, ghat, sigma, ikppkey
186       &     )       &     )
187    cph(
188  CADJ STORE dkm1,blmc,ghat = comlev1_kpp, key = ikey  CADJ STORE dkm1,blmc,ghat = comlev1_kpp, key = ikppkey
189    CADJ STORE hbl, kbl, diffus, casea = comlev1_kpp, key = ikppkey
190    cph)
191    
192  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
193  c enhance diffusivity at interface kbl - 1  c enhance diffusivity at interface kbl - 1
# Line 181  c--------------------------------------- Line 198  c---------------------------------------
198       U     , ghat       U     , ghat
199       O     , blmc )       O     , blmc )
200    
201    cph(
202    cph avoids recomp. of enhance
203    CADJ STORE blmc = comlev1_kpp, key = ikppkey
204    cph)
205    
206  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
207  c combine interior and boundary layer coefficients and nonlocal term  c combine interior and boundary layer coefficients and nonlocal term
208    c !!!NOTE!!! In shallow (2-level) regions and for shallow mixed layers
209    c (< 1 level), diffusivity blmc can become negative.  The max-s below
210    c are a hack until this problem is properly diagnosed and fixed.
211  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
   
212        do k = 1, Nr        do k = 1, Nr
213           do i = 1, imt           do i = 1, imt
214              if (k .lt. kbl(i)) then              if (k .lt. kbl(i)) then
215                 do md = 1, mdiff                 diffus(i,k,1) = max ( blmc(i,k,1), viscAr  )
216                    diffus(i,k,md) = blmc(i,k,md)                 diffus(i,k,2) = max ( blmc(i,k,2), diffusKzS(i,Nr) )
217                 end do                 diffus(i,k,3) = max ( blmc(i,k,3), diffusKzT(i,Nr) )
218              else              else
219                 ghat(i,k) = 0.                 ghat(i,k) = 0.
220              endif              endif
# Line 208  c*************************************** Line 232  c***************************************
232       I       mytime, mythid       I       mytime, mythid
233       I     , kmtj       I     , kmtj
234       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol
235       I     , ikey       I     , ikppkey
236       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma
237       &     )       &     )
238    
# Line 236  c Line 260  c
260        IMPLICIT NONE        IMPLICIT NONE
261    
262  #include "SIZE.h"  #include "SIZE.h"
 #include "EEPARAMS.h"  
 #include "PARAMS.h"  
263  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
 #include "FFIELDS.h"  
264    
265  c input  c input
266  c------  c------
# Line 258  c coriol    : Coriolis parameter Line 279  c coriol    : Coriolis parameter
279        _RL     mytime        _RL     mytime
280        integer mythid        integer mythid
281        integer kmtj(imt)        integer kmtj(imt)
282        _KPP_RL dvsq  (imt,Nr)        _RL dvsq    (imt,Nr)
283        _KPP_RL dbloc (imt,Nr)        _RL dbloc   (imt,Nr)
284        _KPP_RL Ritop (imt,Nr)        _RL Ritop   (imt,Nr)
285        _KPP_RL ustar (imt)        _RL ustar   (imt)
286        _KPP_RL bo    (imt)        _RL bo      (imt)
287        _KPP_RL bosol (imt)        _RL bosol   (imt)
288        _KPP_RL coriol(imt)        _RL coriol  (imt)
289        integer ikey        integer ikppkey, kkppkey
290    
291  c  output  c  output
292  c--------  c--------
# Line 276  c casea     : =1 in case A, =0 in case B Line 297  c casea     : =1 in case A, =0 in case B
297  c kbl       : -1 of first grid level below hbl  c kbl       : -1 of first grid level below hbl
298  c Rib       : Bulk Richardson number  c Rib       : Bulk Richardson number
299  c sigma     : normalized depth (d/hbl)  c sigma     : normalized depth (d/hbl)
300        _KPP_RL hbl   (imt)        _RL hbl    (imt)
301        _KPP_RL bfsfc (imt)        _RL bfsfc  (imt)
302        _KPP_RL stable(imt)        _RL stable (imt)
303        _KPP_RL casea (imt)        _RL casea  (imt)
304        integer kbl   (imt)        integer kbl(imt)
305        _KPP_RL Rib   (imt,Nr)        _RL Rib    (imt,Nr)
306        _KPP_RL sigma (imt)        _RL sigma  (imt)
307    
308  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
309    
310  c  local  c  local
311  c-------  c-------
312  c wm, ws    : turbulent velocity scales         (m/s)  c wm, ws    : turbulent velocity scales         (m/s)
313        _KPP_RL wm(imt), ws(imt)        _RL wm(imt), ws(imt)
314        _RL     worka(imt)        _RL worka(imt)
315    
316        _KPP_RL bvsq, vtsq, hekman, hmonob, hlimit, tempVar1, tempVar2        _RL bvsq, vtsq, hekman, hmonob, hlimit, tempVar1, tempVar2
317        integer i, kl        integer i, kl
318    
319        _KPP_RL     p5    , eins        _RL         p5    , eins
320        parameter ( p5=0.5, eins=1.0 )        parameter ( p5=0.5, eins=1.0 )
321        _RL         minusone        _RL         minusone
322        parameter ( minusone=-1.0 )        parameter ( minusone=-1.0 )
# Line 319  c     initialize hbl and kbl to bottomed Line 340  c     initialize hbl and kbl to bottomed
340    
341        do kl = 2, Nr        do kl = 2, Nr
342    
343    #ifdef ALLOW_AUTODIFF_TAMC
344             kkppkey = (ikppkey-1)*Nr + kl
345    #endif
346    
347  c     compute bfsfc = sw fraction at hbf * zgrid  c     compute bfsfc = sw fraction at hbf * zgrid
348    
349           do i = 1, imt           do i = 1, imt
350              worka(i) = zgrid(kl)              worka(i) = zgrid(kl)
351           end do           end do
352    CADJ store worka = comlev1_kpp_k, key = kkppkey
353           call SWFRAC(           call SWFRAC(
354       I        imt, hbf,       I        imt, hbf,
355       I        mytime, mythid,       I        mytime, mythid,
356       U        worka )       U        worka )
357    CADJ store worka = comlev1_kpp_k, key = kkppkey
358    
359    
360           do i = 1, imt           do i = 1, imt
361    
# Line 350  c--------------------------------------- Line 378  c---------------------------------------
378           call wscale (           call wscale (
379       I        sigma, casea, ustar, bfsfc,       I        sigma, casea, ustar, bfsfc,
380       O        wm, ws )       O        wm, ws )
381    CADJ store ws = comlev1_kpp_k, key = kkppkey
382    
383           do i = 1, imt           do i = 1, imt
384    
# Line 385  c Line 414  c
414    
415           end do           end do
416        end do        end do
417          
418    cph(
419    cph  without this store, there is a recomputation error for
420    cph  rib in adbldepth (probably partial recomputation problem)    
421    CADJ store Rib = comlev1_kpp
422    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy),Nr /)
423    cph)
424    
425        do kl = 2, Nr        do kl = 2, Nr
426           do i = 1, imt           do i = 1, imt
427              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 429  c
429        end do        end do
430    
431  CADJ store kbl = comlev1_kpp  CADJ store kbl = comlev1_kpp
432  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
433    
434        do i = 1, imt        do i = 1, imt
435           kl = kbl(i)           kl = kbl(i)
# Line 406  c     linearly interpolate to find hbl w Line 442  c     linearly interpolate to find hbl w
442        end do        end do
443    
444  CADJ store hbl = comlev1_kpp  CADJ store hbl = comlev1_kpp
445  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
446    
447  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
448  c     find stability and buoyancy forcing for boundary layer  c     find stability and buoyancy forcing for boundary layer
# Line 415  c--------------------------------------- Line 451  c---------------------------------------
451        do i = 1, imt        do i = 1, imt
452           worka(i) = hbl(i)           worka(i) = hbl(i)
453        end do        end do
454    CADJ store worka = comlev1_kpp
455    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
456        call SWFRAC(        call SWFRAC(
457       I     imt, minusone,       I     imt, minusone,
458       I     mytime, mythid,       I     mytime, mythid,
459       U     worka )       U     worka )
460    CADJ store worka = comlev1_kpp
461    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
462    
463        do i = 1, imt        do i = 1, imt
464           bfsfc(i)  = bo(i) + bosol(i) * (1. - worka(i))           bfsfc(i)  = bo(i) + bosol(i) * (1. - worka(i))
465        end do        end do
466  CADJ store bfsfc = comlev1_kpp  CADJ store bfsfc = comlev1_kpp
467  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
468    
469  c--   ensure bfsfc is never 0  c--   ensure bfsfc is never 0
470        do i = 1, imt        do i = 1, imt
# Line 432  c--   ensure bfsfc is never 0 Line 472  c--   ensure bfsfc is never 0
472           bfsfc(i) = sign(eins,bfsfc(i))*max(phepsi,abs(bfsfc(i)))           bfsfc(i) = sign(eins,bfsfc(i))*max(phepsi,abs(bfsfc(i)))
473        end do        end do
474    
475  CADJ store bfsfc = comlev1_kpp  cph(
476  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  cph  added stable to store list to avoid extensive recomp.
477    CADJ store bfsfc, stable = comlev1_kpp
478    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
479    cph)
480    
481  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
482  c check hbl limits for hekman or hmonob  c check hbl limits for hekman or hmonob
# Line 451  c--------------------------------------- Line 494  c---------------------------------------
494           end if           end if
495        end do        end do
496  CADJ store hbl = comlev1_kpp  CADJ store hbl = comlev1_kpp
497  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
498    
499        do i = 1, imt        do i = 1, imt
500           hbl(i) = max(hbl(i),minKPPhbl)           hbl(i) = max(hbl(i),minKPPhbl)
# Line 459  CADJ &   , key = ikey, shape = (/ (sNx+2 Line 502  CADJ &   , key = ikey, shape = (/ (sNx+2
502        end do        end do
503    
504  CADJ store hbl = comlev1_kpp  CADJ store hbl = comlev1_kpp
505  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
506    
507  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
508  c      find new kbl  c      find new kbl
# Line 480  c--------------------------------------- Line 523  c---------------------------------------
523        do i = 1, imt        do i = 1, imt
524           worka(i) = hbl(i)           worka(i) = hbl(i)
525        end do        end do
526    CADJ store worka = comlev1_kpp
527    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
528        call SWFRAC(        call SWFRAC(
529       I     imt, minusone,       I     imt, minusone,
530       I     mytime, mythid,       I     mytime, mythid,
531       U     worka )       U     worka )
532    CADJ store worka = comlev1_kpp
533    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
534    
535        do i = 1, imt        do i = 1, imt
536           bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i))           bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i))
537        end do        end do
538  CADJ store bfsfc = comlev1_kpp  CADJ store bfsfc = comlev1_kpp
539  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
540    
541  c--   ensures bfsfc is never 0  c--   ensures bfsfc is never 0
542        do i = 1, imt        do i = 1, imt
# Line 536  c sigma   : normalized depth (d/hbl) Line 583  c sigma   : normalized depth (d/hbl)
583  c hbl     : boundary layer depth (m)  c hbl     : boundary layer depth (m)
584  c ustar   : surface friction velocity         (m/s)  c ustar   : surface friction velocity         (m/s)
585  c bfsfc   : total surface buoyancy flux       (m^2/s^3)  c bfsfc   : total surface buoyancy flux       (m^2/s^3)
586        _KPP_RL sigma(imt)        _RL sigma(imt)
587        _KPP_RL hbl  (imt)        _RL hbl  (imt)
588        _KPP_RL ustar(imt)        _RL ustar(imt)
589        _KPP_RL bfsfc(imt)        _RL bfsfc(imt)
590    
591  c  output  c  output
592  c--------  c--------
593  c wm, ws  : turbulent velocity scales at sigma  c wm, ws  : turbulent velocity scales at sigma
594        _KPP_RL wm(imt), ws(imt)        _RL wm(imt), ws(imt)
595    
596  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
597    
598  c local  c local
599  c------  c------
600  c zehat   : = zeta *  ustar**3  c zehat   : = zeta *  ustar**3
601        _KPP_RL zehat        _RL zehat
602    
603        integer iz, izp1, ju, i, jup1        integer iz, izp1, ju, i, jup1
604        _KPP_RL udiff, zdiff, zfrac, ufrac, fzfrac, wam        _RL udiff, zdiff, zfrac, ufrac, fzfrac, wam
605        _KPP_RL wbm, was, wbs, u3, tempVar        _RL wbm, was, wbs, u3, tempVar
606    
607  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
608  c use lookup table for zehat < zmax only; otherwise use  c use lookup table for zehat < zmax only; otherwise use
# Line 611  c*************************************** Line 658  c***************************************
658    
659        subroutine Ri_iwmix (        subroutine Ri_iwmix (
660       I       kmtj, shsq, dbloc, dblocSm       I       kmtj, shsq, dbloc, dblocSm
661       I     , ikey       I     , diffusKzS, diffusKzT
662         I     , ikppkey
663       O     , diffus )       O     , diffus )
664    
665  c     compute interior viscosity diffusivity coefficients due  c     compute interior viscosity diffusivity coefficients due
# Line 631  c     kmtj   (imt)         number of ver Line 679  c     kmtj   (imt)         number of ver
679  c     shsq   (imt,Nr)      (local velocity shear)^2               ((m/s)^2)  c     shsq   (imt,Nr)      (local velocity shear)^2               ((m/s)^2)
680  c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)  c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)
681  c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)  c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)
682        integer kmtj   (imt)  c     diffusKzS(imt,Nr)- background vertical diffusivity for scalars    (m^2/s)
683        _KPP_RL shsq   (imt,Nr)  c     diffusKzT(imt,Nr)- background vertical diffusivity for theta      (m^2/s)
684        _KPP_RL dbloc  (imt,Nr)        integer kmtj (imt)
685        _KPP_RL dblocSm(imt,Nr)        _RL shsq     (imt,Nr)
686        integer ikey        _RL dbloc    (imt,Nr)
687          _RL dblocSm  (imt,Nr)
688          _RL diffusKzS(imt,Nr)
689          _RL diffusKzT(imt,Nr)
690          integer ikppkey
691    
692  c output  c output
693  c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)  c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)
694  c     diffus(imt,0:Nrp1,2)  vertical scalar diffusivity             (m^2/s)  c     diffus(imt,0:Nrp1,2)  vertical scalar diffusivity             (m^2/s)
695  c     diffus(imt,0:Nrp1,3)  vertical temperature diffusivity        (m^2/s)  c     diffus(imt,0:Nrp1,3)  vertical temperature diffusivity        (m^2/s)
696        _KPP_RL diffus(imt,0:Nrp1,3)        _RL diffus(imt,0:Nrp1,3)
697    
698  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
699    
700  c local variables  c local variables
701  c     Rig                   local Richardson number  c     Rig                   local Richardson number
702  c     fRi, fcon             function of Rig  c     fRi, fcon             function of Rig
703        _KPP_RL Rig        _RL Rig
704        _KPP_RL fRi, fcon        _RL fRi, fcon
705        _KPP_RL ratio        _RL ratio
706        integer i, ki, mr        integer i, ki
707        _KPP_RL c1, c0        _RL c1, c0
708    
709  #ifdef ALLOW_KPP_VERTICALLY_SMOOTH  #ifdef ALLOW_KPP_VERTICALLY_SMOOTH
710          integer mr
711  CADJ INIT kpp_ri_tape_mr = common, 1  CADJ INIT kpp_ri_tape_mr = common, 1
712  #endif  #endif
713    
# Line 670  c     set values at bottom and below to Line 723  c     set values at bottom and below to
723  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
724  C     break data flow dependence on diffus  C     break data flow dependence on diffus
725        diffus(1,1,1) = 0.0        diffus(1,1,1) = 0.0
726    
727          do ki = 1, Nr
728             do i = 1, imt
729                diffus(i,ki,1) = 0.
730                diffus(i,ki,2) = 0.
731                diffus(i,ki,3) = 0.
732             enddo
733          enddo
734  #endif  #endif
735                
736    
737        do ki = 1, Nr        do ki = 1, Nr
738           do i = 1, imt           do i = 1, imt
739              if     (kmtj(i) .EQ. 0      ) then              if     (kmtj(i) .LE. 1      ) then
740                 diffus(i,ki,1) = 0.                 diffus(i,ki,1) = 0.
741                 diffus(i,ki,2) = 0.                 diffus(i,ki,2) = 0.
742              elseif (ki      .GE. kmtj(i)) then              elseif (ki      .GE. kmtj(i)) then
# Line 687  C     break data flow dependence on diff Line 749  C     break data flow dependence on diff
749              endif              endif
750           end do           end do
751        end do        end do
752    CADJ store diffus = comlev1_kpp, key = ikppkey
753    
754  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
755  c     vertically smooth Ri  c     vertically smooth Ri
# Line 724  c  evaluate f of smooth Ri for shear ins Line 787  c  evaluate f of smooth Ri for shear ins
787  c ----------------------------------------------------------------------  c ----------------------------------------------------------------------
788  c            evaluate diffusivities and viscosity  c            evaluate diffusivities and viscosity
789  c    mixing due to internal waves, and shear and static instability  c    mixing due to internal waves, and shear and static instability
790    
791              diffus(i,ki,1) = viscAr + fcon * difmcon + fRi * difm0  #ifndef EXCLUDE_KPP_SHEAR_MIX
792              diffus(i,ki,2) = diffKrS + fcon * difscon + fRi * difs0              if ( .NOT. inAdMode ) then
793              diffus(i,ki,3) = diffKrT + fcon * difscon + fRi * difs0                 diffus(i,ki,1) = viscAr + fcon * difmcon + fRi * difm0
794                   diffus(i,ki,2) = diffusKzS(i,ki)+fcon*difscon+fRi*difs0
795                   diffus(i,ki,3) = diffusKzT(i,ki)+fcon*diftcon+fRi*dift0
796                else
797                   diffus(i,ki,1) = viscAr
798                   diffus(i,ki,2) = diffusKzS(i,ki)
799                   diffus(i,ki,3) = diffusKzT(i,ki)
800                endif
801    #else
802                diffus(i,ki,1) = viscAr
803                diffus(i,ki,2) = diffusKzS(i,ki)
804                diffus(i,ki,3) = diffusKzT(i,ki)
805    #endif
806    
807           end do           end do
808        end do        end do
# Line 767  c     penetrative convection. Line 842  c     penetrative convection.
842  c input/output  c input/output
843  c-------------  c-------------
844  c v     : 2-D array to be smoothed in Nrp1 direction  c v     : 2-D array to be smoothed in Nrp1 direction
845        _KPP_RL v(imt,0:Nrp1)        _RL v(imt,0:Nrp1)
846    
847  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
848    
849  c local  c local
850        _KPP_RL zwork, zflag        _RL zwork, zflag
851        _KPP_RL KRi_range(1:Nrp1)        _RL KRi_range(1:Nrp1)
852        integer i, k, km1, kp1        integer i, k, km1, kp1
853    
854        _KPP_RL     p0      , p25       , p5      , p2        _RL         p0      , p25       , p5      , p2
855        parameter ( p0 = 0.0, p25 = 0.25, p5 = 0.5, p2 = 2.0 )        parameter ( p0 = 0.0, p25 = 0.25, p5 = 0.5, p2 = 2.0 )
856    
857        KRi_range(Nrp1) = p0        KRi_range(Nrp1) = p0
# Line 830  CADJ STORE v(i,k), zwork = z121tape Line 905  CADJ STORE v(i,k), zwork = z121tape
905    
906  c*************************************************************************  c*************************************************************************
907    
       subroutine kpp_smooth_horiz (  
      I     k, bi, bj,  
      U     fld )  
   
 c     Apply horizontal smoothing to KPP array  
   
       IMPLICIT NONE  
 #include "SIZE.h"  
 #include "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*************************************************************************  
   
908        subroutine smooth_horiz (        subroutine smooth_horiz (
909       I     k, bi, bj,       I     k, bi, bj,
910       U     fld )       U     fld )
# Line 918  c     Apply horizontal smoothing to glob Line 913  c     Apply horizontal smoothing to glob
913    
914        IMPLICIT NONE        IMPLICIT NONE
915  #include "SIZE.h"  #include "SIZE.h"
916    #include "GRID.h"
917  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
918    
919  c     input  c     input
# Line 949  c     local Line 945  c     local
945              im1 = i-1              im1 = i-1
946              ip1 = i+1              ip1 = i+1
947              tempVar =              tempVar =
948       &           p25   *   pMask(i  ,j  ,k,bi,bj)   +       &           p25   *   maskC(i  ,j  ,k,bi,bj)   +
949       &           p125  * ( pMask(im1,j  ,k,bi,bj)   +       &           p125  * ( maskC(im1,j  ,k,bi,bj)   +
950       &                     pMask(ip1,j  ,k,bi,bj)   +       &                     maskC(ip1,j  ,k,bi,bj)   +
951       &                     pMask(i  ,jm1,k,bi,bj)   +       &                     maskC(i  ,jm1,k,bi,bj)   +
952       &                     pMask(i  ,jp1,k,bi,bj) ) +       &                     maskC(i  ,jp1,k,bi,bj) ) +
953       &           p0625 * ( pMask(im1,jm1,k,bi,bj)   +       &           p0625 * ( maskC(im1,jm1,k,bi,bj)   +
954       &                     pMask(im1,jp1,k,bi,bj)   +       &                     maskC(im1,jp1,k,bi,bj)   +
955       &                     pMask(ip1,jm1,k,bi,bj)   +       &                     maskC(ip1,jm1,k,bi,bj)   +
956       &                     pMask(ip1,jp1,k,bi,bj) )       &                     maskC(ip1,jp1,k,bi,bj) )
957              IF ( tempVar .GE. p25 ) THEN              IF ( tempVar .GE. p25 ) THEN
958                 fld_tmp(i,j) = (                 fld_tmp(i,j) = (
959       &              p25  * fld(i  ,j  )*pMask(i  ,j  ,k,bi,bj) +       &              p25  * fld(i  ,j  )*maskC(i  ,j  ,k,bi,bj) +
960       &              p125 *(fld(im1,j  )*pMask(im1,j  ,k,bi,bj) +       &              p125 *(fld(im1,j  )*maskC(im1,j  ,k,bi,bj) +
961       &                     fld(ip1,j  )*pMask(ip1,j  ,k,bi,bj) +       &                     fld(ip1,j  )*maskC(ip1,j  ,k,bi,bj) +
962       &                     fld(i  ,jm1)*pMask(i  ,jm1,k,bi,bj) +       &                     fld(i  ,jm1)*maskC(i  ,jm1,k,bi,bj) +
963       &                     fld(i  ,jp1)*pMask(i  ,jp1,k,bi,bj))+       &                     fld(i  ,jp1)*maskC(i  ,jp1,k,bi,bj))+
964       &              p0625*(fld(im1,jm1)*pMask(im1,jm1,k,bi,bj) +       &              p0625*(fld(im1,jm1)*maskC(im1,jm1,k,bi,bj) +
965       &                     fld(im1,jp1)*pMask(im1,jp1,k,bi,bj) +       &                     fld(im1,jp1)*maskC(im1,jp1,k,bi,bj) +
966       &                     fld(ip1,jm1)*pMask(ip1,jm1,k,bi,bj) +       &                     fld(ip1,jm1)*maskC(ip1,jm1,k,bi,bj) +
967       &                     fld(ip1,jp1)*pMask(ip1,jp1,k,bi,bj)))       &                     fld(ip1,jp1)*maskC(ip1,jp1,k,bi,bj)))
968       &              / tempVar       &              / tempVar
969              ELSE              ELSE
970                 fld_tmp(i,j) = fld(i,j)                 fld_tmp(i,j) = fld(i,j)
# Line 992  c*************************************** Line 988  c***************************************
988    
989        subroutine blmix (        subroutine blmix (
990       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl
991       O     , dkm1, blmc, ghat, sigma, ikey       O     , dkm1, blmc, ghat, sigma, ikppkey
992       &     )       &     )
993    
994  c     mixing coefficients within boundary layer depend on surface  c     mixing coefficients within boundary layer depend on surface
# Line 1015  c     hbl   (imt)                 bounda Line 1011  c     hbl   (imt)                 bounda
1011  c     stable(imt)                 = 1 in stable forcing  c     stable(imt)                 = 1 in stable forcing
1012  c     casea (imt)                 = 1 in case A  c     casea (imt)                 = 1 in case A
1013  c     diffus(imt,0:Nrp1,mdiff)    vertical diffusivities              (m^2/s)  c     diffus(imt,0:Nrp1,mdiff)    vertical diffusivities              (m^2/s)
1014  c     kbl(imt)                    -1 of first grid level below hbl  c     kbl   (imt)                 -1 of first grid level below hbl
1015        _KPP_RL ustar (imt)        _RL ustar (imt)
1016        _KPP_RL bfsfc (imt)        _RL bfsfc (imt)
1017        _KPP_RL hbl   (imt)        _RL hbl   (imt)
1018        _KPP_RL stable(imt)        _RL stable(imt)
1019        _KPP_RL casea (imt)        _RL casea (imt)
1020        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL diffus(imt,0:Nrp1,mdiff)
1021        integer kbl(imt)        integer kbl(imt)
1022    
1023  c output  c output
# Line 1029  c     dkm1 (imt,mdiff)            bounda Line 1025  c     dkm1 (imt,mdiff)            bounda
1025  c     blmc (imt,Nr,mdiff)         boundary layer mixing coefficients  (m^2/s)  c     blmc (imt,Nr,mdiff)         boundary layer mixing coefficients  (m^2/s)
1026  c     ghat (imt,Nr)               nonlocal scalar transport  c     ghat (imt,Nr)               nonlocal scalar transport
1027  c     sigma(imt)                  normalized depth (d / hbl)  c     sigma(imt)                  normalized depth (d / hbl)
1028        _KPP_RL dkm1 (imt,mdiff)        _RL dkm1 (imt,mdiff)
1029        _KPP_RL blmc (imt,Nr,mdiff)        _RL blmc (imt,Nr,mdiff)
1030        _KPP_RL ghat (imt,Nr)        _RL ghat (imt,Nr)
1031        _KPP_RL sigma(imt)        _RL sigma(imt)
1032        integer ikey        integer ikppkey, kkppkey
1033    
1034  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1035    
# Line 1041  c  local Line 1037  c  local
1037  c     gat1*(imt)                 shape function at sigma = 1  c     gat1*(imt)                 shape function at sigma = 1
1038  c     dat1*(imt)                 derivative of shape function at sigma = 1  c     dat1*(imt)                 derivative of shape function at sigma = 1
1039  c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)  c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)
1040        _KPP_RL gat1m(imt), gat1s(imt), gat1t(imt)        _RL gat1m(imt), gat1s(imt), gat1t(imt)
1041        _KPP_RL dat1m(imt), dat1s(imt), dat1t(imt)        _RL dat1m(imt), dat1s(imt), dat1t(imt)
1042        _KPP_RL ws(imt), wm(imt)        _RL ws(imt), wm(imt)
1043        integer i, kn, ki        integer i, kn, ki
1044        _KPP_RL R, dvdzup, dvdzdn, viscp        _RL R, dvdzup, dvdzdn, viscp
1045        _KPP_RL difsp, diftp, visch, difsh, difth        _RL difsp, diftp, visch, difsh, difth
1046        _KPP_RL f1, sig, a1, a2, a3, delhat        _RL f1, sig, a1, a2, a3, delhat
1047        _KPP_RL Gm, Gs, Gt        _RL Gm, Gs, Gt
1048        _KPP_RL tempVar        _RL tempVar
1049    
1050        _KPP_RL    p0    , eins        _RL    p0    , eins
1051        parameter (p0=0.0, eins=1.0)        parameter (p0=0.0, eins=1.0)
1052    
1053  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 1062  c--------------------------------------- Line 1058  c---------------------------------------
1058           sigma(i) = stable(i) * 1.0 + (1. - stable(i)) * epsilon           sigma(i) = stable(i) * 1.0 + (1. - stable(i)) * epsilon
1059        end do        end do
1060    
1061    CADJ STORE sigma = comlev1_kpp, key = ikppkey
1062        call wscale (        call wscale (
1063       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1064       O        wm, ws )       O        wm, ws )
1065    CADJ STORE wm = comlev1_kpp, key = ikppkey
1066    CADJ STORE ws = comlev1_kpp, key = ikppkey
1067    
1068        do i = 1, imt        do i = 1, imt
1069           wm(i) = sign(eins,wm(i))*max(phepsi,abs(wm(i)))           wm(i) = sign(eins,wm(i))*max(phepsi,abs(wm(i)))
1070           ws(i) = sign(eins,ws(i))*max(phepsi,abs(ws(i)))           ws(i) = sign(eins,ws(i))*max(phepsi,abs(ws(i)))
1071        end do        end do
1072  CADJ STORE wm = comlev1_kpp, key = ikey  CADJ STORE wm = comlev1_kpp, key = ikppkey
1073  CADJ STORE ws = comlev1_kpp, key = ikey  CADJ STORE ws = comlev1_kpp, key = ikppkey
1074    
1075        do i = 1, imt        do i = 1, imt
1076    
# Line 1107  c--------------------------------------- Line 1106  c---------------------------------------
1106       &        max(ustar(i)**4,phepsi)       &        max(ustar(i)**4,phepsi)
1107           gat1m(i) = visch / hbl(i) / wm(i)           gat1m(i) = visch / hbl(i) / wm(i)
1108           dat1m(i) = -viscp / wm(i) + f1 * visch           dat1m(i) = -viscp / wm(i) + f1 * visch
          dat1m(i) = min(dat1m(i),p0)  
1109    
1110           gat1s(i) = difsh  / hbl(i) / ws(i)           gat1s(i) = difsh  / hbl(i) / ws(i)
1111           dat1s(i) = -difsp / ws(i) + f1 * difsh           dat1s(i) = -difsp / ws(i) + f1 * difsh
          dat1s(i) = min(dat1s(i),p0)  
1112    
1113           gat1t(i) = difth /  hbl(i) / ws(i)           gat1t(i) = difth /  hbl(i) / ws(i)
1114           dat1t(i) = -diftp / ws(i) + f1 * difth           dat1t(i) = -diftp / ws(i) + f1 * difth
          dat1t(i) = min(dat1t(i),p0)  
1115    
1116        end do        end do
1117    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1118    CADJ STORE gat1m = comlev1_kpp, key = ikppkey
1119    CADJ STORE gat1s = comlev1_kpp, key = ikppkey
1120    CADJ STORE gat1t = comlev1_kpp, key = ikppkey
1121    CADJ STORE dat1m = comlev1_kpp, key = ikppkey
1122    CADJ STORE dat1s = comlev1_kpp, key = ikppkey
1123    CADJ STORE dat1t = comlev1_kpp, key = ikppkey
1124    #endif
1125          do i = 1, imt
1126             dat1m(i) = min(dat1m(i),p0)
1127             dat1s(i) = min(dat1s(i),p0)
1128             dat1t(i) = min(dat1t(i),p0)
1129          end do
1130    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1131    CADJ STORE dat1m = comlev1_kpp, key = ikppkey
1132    CADJ STORE dat1s = comlev1_kpp, key = ikppkey
1133    CADJ STORE dat1t = comlev1_kpp, key = ikppkey
1134    #endif
1135    
1136        do ki = 1, Nr        do ki = 1, Nr
1137    
1138    #ifdef ALLOW_AUTODIFF_TAMC
1139             kkppkey = (ikppkey-1)*Nr + ki
1140    #endif
1141    
1142  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1143  c     compute turbulent velocity scales on the interfaces  c     compute turbulent velocity scales on the interfaces
1144  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 1129  c--------------------------------------- Line 1147  c---------------------------------------
1147              sig      = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i)              sig      = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i)
1148              sigma(i) = stable(i)*sig + (1.-stable(i))*min(sig,epsilon)              sigma(i) = stable(i)*sig + (1.-stable(i))*min(sig,epsilon)
1149           end do           end do
1150    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1151    CADJ STORE wm = comlev1_kpp_k, key = kkppkey
1152    CADJ STORE ws = comlev1_kpp_k, key = kkppkey
1153    #endif
1154    CADJ STORE sigma = comlev1_kpp_k, key = kkppkey
1155           call wscale (           call wscale (
1156       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1157       O        wm, ws )       O        wm, ws )
1158    CADJ STORE wm = comlev1_kpp_k, key = kkppkey
1159    CADJ STORE ws = comlev1_kpp_k, key = kkppkey
1160    
1161  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1162  c     compute the dimensionless shape functions at the interfaces  c     compute the dimensionless shape functions at the interfaces
# Line 1175  c--------------------------------------- Line 1200  c---------------------------------------
1200       &            + (1. - stable(i)) * min(sig,epsilon)       &            + (1. - stable(i)) * min(sig,epsilon)
1201        end do        end do
1202    
1203    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1204    CADJ STORE wm = comlev1_kpp, key = ikppkey
1205    CADJ STORE ws = comlev1_kpp, key = ikppkey
1206    #endif
1207    CADJ STORE sigma = comlev1_kpp, key = ikppkey
1208        call wscale (        call wscale (
1209       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1210       O        wm, ws )       O        wm, ws )
1211    CADJ STORE wm = comlev1_kpp, key = ikppkey
1212    CADJ STORE ws = comlev1_kpp, key = ikppkey
1213    
1214        do i = 1, imt        do i = 1, imt
1215           sig = -zgrid(kbl(i)-1) / hbl(i)           sig = -zgrid(kbl(i)-1) / hbl(i)
# Line 1218  c     hbl(imt)                  boundary Line 1250  c     hbl(imt)                  boundary
1250  c     kbl(imt)                  grid above hbl  c     kbl(imt)                  grid above hbl
1251  c     diffus(imt,0:Nrp1,mdiff) vertical diffusivities           (m^2/s)  c     diffus(imt,0:Nrp1,mdiff) vertical diffusivities           (m^2/s)
1252  c     casea(imt)                = 1 in caseA, = 0 in case B  c     casea(imt)                = 1 in caseA, = 0 in case B
1253        _KPP_RL dkm1  (imt,mdiff)        _RL dkm1  (imt,mdiff)
1254        _KPP_RL hbl   (imt)        _RL hbl   (imt)
1255        integer kbl   (imt)        integer kbl   (imt)
1256        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL diffus(imt,0:Nrp1,mdiff)
1257        _KPP_RL casea (imt)        _RL casea (imt)
1258    
1259  c input/output  c input/output
1260  c     nonlocal transport, modified ghat at kbl(i)-1 interface    (s/m**2)  c     nonlocal transport, modified ghat at kbl(i)-1 interface    (s/m**2)
1261        _KPP_RL ghat (imt,Nr)        _RL ghat (imt,Nr)
1262    
1263  c output  c output
1264  c     enhanced bound. layer mixing coeff.  c     enhanced bound. layer mixing coeff.
1265        _KPP_RL blmc  (imt,Nr,mdiff)        _RL blmc  (imt,Nr,mdiff)
1266    
1267  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1268    
1269  c local  c local
1270  c     fraction hbl lies beteen zgrid neighbors  c     fraction hbl lies beteen zgrid neighbors
1271        _KPP_RL delta        _RL delta
1272        integer ki, i, md        integer ki, i, md
1273        _KPP_RL dkmp5, dstar        _RL dkmp5, dstar
1274    
1275        do i = 1, imt        do i = 1, imt
1276           ki = kbl(i)-1           ki = kbl(i)-1
# Line 1264  c     fraction hbl lies beteen zgrid nei Line 1296  c     fraction hbl lies beteen zgrid nei
1296  c*************************************************************************  c*************************************************************************
1297    
1298        SUBROUTINE STATEKPP (        SUBROUTINE STATEKPP (
1299       I     bi, bj, myThid,       I     ikppkey, bi, bj, myThid,
1300       O     RHO1, DBLOC, DBSFC, TTALPHA, SSBETA)       O     RHO1, DBLOC, DBSFC, TTALPHA, SSBETA)
1301  c  c
1302  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 1290  c Line 1322  c
1322  c     written  by: jan morzel,   feb. 10, 1995 (converted from "sigma" version)  c     written  by: jan morzel,   feb. 10, 1995 (converted from "sigma" version)
1323  c     modified by: d. menemenlis,    june 1998 : for use with MIT GCM UV  c     modified by: d. menemenlis,    june 1998 : for use with MIT GCM UV
1324  c  c
1325    c     28 april 05: added computation of mixed layer depth KPPmld
1326    c     for a deltaRho equivalent to a deltaTheta of 0.8 deg C
1327    
1328  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1329    
1330        IMPLICIT NONE        IMPLICIT NONE
# Line 1299  c--------------------------------------- Line 1334  c---------------------------------------
1334  #include "PARAMS.h"  #include "PARAMS.h"
1335  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
1336  #include "DYNVARS.h"  #include "DYNVARS.h"
1337    #include "GRID.h"
1338    
1339  c-------------- Routine arguments -----------------------------------------  c-------------- Routine arguments -----------------------------------------
1340        INTEGER bi, bj, myThid        INTEGER bi, bj, myThid
1341        _KPP_RL RHO1   ( ibot:itop, jbot:jtop       )        _RL RHO1   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy       )
1342        _KPP_RL DBLOC  ( ibot:itop, jbot:jtop, Nr   )        _RL DBLOC  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
1343        _KPP_RL DBSFC  ( ibot:itop, jbot:jtop, Nr   )        _RL DBSFC  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
1344        _KPP_RL TTALPHA( ibot:itop, jbot:jtop, Nrp1 )        _RL TTALPHA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1345        _KPP_RL SSBETA ( ibot:itop, jbot:jtop, Nrp1 )        _RL SSBETA ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1346    
1347  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1348    
# Line 1317  c Line 1353  c
1353  c     rhok         - density of t(k  ) & s(k  ) at depth k  c     rhok         - density of t(k  ) & s(k  ) at depth k
1354  c     rhokm1       - density of t(k-1) & s(k-1) at depth k  c     rhokm1       - density of t(k-1) & s(k-1) at depth k
1355  c     rho1k        - density of t(1  ) & s(1  ) at depth k  c     rho1k        - density of t(1  ) & s(1  ) at depth k
1356  c     work1, work2 - work arrays for holding horizontal slabs  c     work1,2,3    - work arrays for holding horizontal slabs
1357    
1358        _RL RHOK  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL RHOK  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1359        _RL RHOKM1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL RHOKM1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 1325  c     work1, work2 - work arrays for hol Line 1361  c     work1, work2 - work arrays for hol
1361        _RL WORK1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1362        _RL WORK2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1363        _RL WORK3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1364    #ifdef ALLOW_DIAGNOSTICS
1365    c     KPPMLD - mixed layer depth based on density criterion
1366          _RL KPPMLD(1    :sNx    ,1    :sNy    )
1367    #endif /* ALLOW_DIAGNOSTICS */
1368    
1369        INTEGER I, J, K        INTEGER I, J, K
1370          INTEGER ikppkey, kkppkey
1371    
1372  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
1373    
1374          kkppkey = (ikppkey-1)*Nr + 1
1375    
1376    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1377    CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1378    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1379    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1380        call FIND_RHO(        call FIND_RHO(
1381       I     bi, bj, ibot, itop, jbot, jtop, 1, 1, eosType,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1382       I     theta, salt,       I     theta, salt,
1383       O     WORK1,       O     WORK1,
1384       I     myThid )       I     myThid )
1385    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1386    CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1387    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1388    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1389    
1390        call FIND_ALPHA(        call FIND_ALPHA(
1391       I     bi, bj, ibot, itop, jbot, jtop, 1, 1, eosType,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1392       O     WORK2 )       O     WORK2 )
1393    
1394        call FIND_BETA(        call FIND_BETA(
1395       I     bi, bj, ibot, itop, jbot, jtop, 1, 1, eosType,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1396       O     WORK3 )       O     WORK3 )
1397    
1398        DO J = jbot, jtop        DO J = 1-OLy, sNy+OLy
1399           DO I = ibot, itop           DO I = 1-OLx, sNx+OLx
1400              RHO1(I,J)      = WORK1(I,J) + rhonil              RHO1(I,J)      = WORK1(I,J) + rhoConst
1401              TTALPHA(I,J,1) = WORK2(I,J)              TTALPHA(I,J,1) = WORK2(I,J)
1402              SSBETA(I,J,1)  = WORK3(I,J)              SSBETA(I,J,1)  = WORK3(I,J)
1403              DBSFC(I,J,1)   = 0.              DBSFC(I,J,1)   = 0.
1404           END DO           END DO
1405        END DO        END DO
1406    
1407    #ifdef ALLOW_DIAGNOSTICS
1408    c     work3 - density of t(1)-.8 & s(1  ) at depth 1
1409          IF ( useDiagnostics ) THEN
1410             DO J = 1, sNy
1411                DO I = 1, sNx
1412                   KPPMLD(I,J) = ABS(R_low(I,J,bi,bj))
1413                   WORK3 (I,J) = WORK1(I,J) - 0.8 * WORK2(I,J)
1414                END DO
1415             END DO
1416          ENDIF
1417    #endif /* ALLOW_DIAGNOSTICS */
1418    
1419  c calculate alpha, beta, and gradients in interior layers  c calculate alpha, beta, and gradients in interior layers
1420    
1421  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)
1422        DO K = 2, Nr        DO K = 2, Nr
1423    
1424          kkppkey = (ikppkey-1)*Nr + k
1425    
1426    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1427    CADJ STORE theta(:,:,k,bi,bj) = comlev1_kpp_k, key=kkppkey
1428    CADJ STORE salt (:,:,k,bi,bj) = comlev1_kpp_k, key=kkppkey
1429    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1430           call FIND_RHO(           call FIND_RHO(
1431       I        bi, bj, ibot, itop, jbot, jtop, K, K, eosType,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1432       I        theta, salt,       I        theta, salt,
1433       O        RHOK,       O        RHOK,
1434       I        myThid )       I        myThid )
1435    
1436    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1437    CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_kpp_k, key=kkppkey
1438    CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_kpp_k, key=kkppkey
1439    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1440           call FIND_RHO(           call FIND_RHO(
1441       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,
1442       I        theta, salt,       I        theta, salt,
1443       O        RHOKM1,       O        RHOKM1,
1444       I        myThid )       I        myThid )
1445    
1446    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1447    CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1448    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1449    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1450           call FIND_RHO(           call FIND_RHO(
1451       I        bi, bj, ibot, itop, jbot, jtop, 1, K, eosType,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, K,
1452       I        theta, salt,       I        theta, salt,
1453       O        RHO1K,       O        RHO1K,
1454       I        myThid )       I        myThid )
1455    
1456    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1457    CADJ STORE rhok  (:,:)          = comlev1_kpp_k, key=kkppkey
1458    CADJ STORE rhokm1(:,:)          = comlev1_kpp_k, key=kkppkey
1459    CADJ STORE rho1k (:,:)          = comlev1_kpp_k, key=kkppkey
1460    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1461    
1462           call FIND_ALPHA(           call FIND_ALPHA(
1463       I        bi, bj, ibot, itop, jbot, jtop, K, K, eosType,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1464       O        WORK1 )       O        WORK1 )
1465    
1466           call FIND_BETA(           call FIND_BETA(
1467       I        bi, bj, ibot, itop, jbot, jtop, K, K, eosType,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1468       O        WORK2 )       O        WORK2 )
1469    
1470           DO J = jbot, jtop           DO J = 1-OLy, sNy+OLy
1471              DO I = ibot, itop              DO I = 1-OLx, sNx+OLx
1472                 TTALPHA(I,J,K) = WORK1 (I,J)                 TTALPHA(I,J,K) = WORK1 (I,J)
1473                 SSBETA(I,J,K)  = WORK2 (I,J)                 SSBETA(I,J,K)  = WORK2 (I,J)
1474                 DBLOC(I,J,K-1) = gravity * (RHOK(I,J) - RHOKM1(I,J)) /                 DBLOC(I,J,K-1) = gravity * (RHOK(I,J) - RHOKM1(I,J)) /
1475       &                                    (RHOK(I,J) + rhonil)       &                                    (RHOK(I,J) + rhoConst)
1476                 DBSFC(I,J,K)   = gravity * (RHOK(I,J) - RHO1K (I,J)) /                 DBSFC(I,J,K)   = gravity * (RHOK(I,J) - RHO1K (I,J)) /
1477       &                                    (RHOK(I,J) + rhonil)       &                                    (RHOK(I,J) + rhoConst)
1478              END DO              END DO
1479           END DO           END DO
1480    
1481    #ifdef ALLOW_DIAGNOSTICS
1482             IF ( useDiagnostics ) THEN
1483    c     work1 - density of t(k-1)  & s(k-1) at depth 1
1484    c     work2 - density of t(k  )  & s(k  ) at depth 1
1485    c     work3 - density of t(1)-.8 & s(1  ) at depth 1
1486                call FIND_RHO(
1487         I           bi, bj, 1, sNx, 1, sNy, K-1, 1, theta, salt,
1488         O           WORK1,
1489         I           myThid )
1490                call FIND_RHO(
1491         I           bi, bj, 1, sNx, 1, sNy, K  , 1, theta, salt,
1492         O           WORK2,
1493         I           myThid )
1494                DO J = 1, sNy
1495                   DO I = 1, sNx
1496                      IF ( k .LE. klowC(I,J,bi,bj) .AND.
1497         &                 WORK1(I,J) .LT. WORK3(I,J) .AND.
1498         &                 WORK2(I,J) .GE. WORK3(I,J) ) THEN
1499                         KPPMLD(I,J) = MIN ( KPPMLD(I,J),
1500         &                    ABS(((WORK3(I,J)-WORK1(I,J))*rC(k)+
1501         &                    (WORK2(I,J)-WORK3(I,J))*rC(k-1))/
1502         &                    (WORK2(I,J)-WORK1(I,J))))
1503                      ENDIF
1504                   END DO
1505                END DO
1506             ENDIF
1507    #endif /* ALLOW_DIAGNOSTICS */
1508    
1509        END DO        END DO
1510    
1511  c     compute arrays for K = Nrp1  c     compute arrays for K = Nrp1
1512        DO J = jbot, jtop        DO J = 1-OLy, sNy+OLy
1513           DO I = ibot, itop           DO I = 1-OLx, sNx+OLx
1514              TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)              TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)
1515              SSBETA(I,J,Nrp1)  = SSBETA(I,J,Nr)              SSBETA(I,J,Nrp1)  = SSBETA(I,J,Nr)
1516              DBLOC(I,J,Nr)     = 0.              DBLOC(I,J,Nr)     = 0.
1517           END DO           END DO
1518        END DO        END DO
1519    
1520    #ifdef ALLOW_DIAGNOSTICS
1521          IF ( useDiagnostics ) THEN
1522             CALL DIAGNOSTICS_FILL(KPPmld,'KPPmld  ',0,1,3,bi,bj,myThid)
1523          ENDIF
1524    #endif /* ALLOW_DIAGNOSTICS */
1525    
1526  #endif /* ALLOW_KPP */  #endif /* ALLOW_KPP */
1527    
1528        RETURN        RETURN
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.29
  ViewVC Help
Powered by ViewVC 1.1.22