/[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.52 by heimbach, Wed May 21 10:45:33 2014 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "KPP_OPTIONS.h"  #include "KPP_OPTIONS.h"
5    #ifdef ALLOW_SALT_PLUME
6    #include "SALT_PLUME_OPTIONS.h"
7    #endif
8    
9  C-- File kpp_routines.F: subroutines needed to implement  C-- File kpp_routines.F: subroutines needed to implement
10  C--                      KPP vertical mixing scheme  C--                      KPP vertical mixing scheme
# Line 11  C--  o BLDEPTH     - Determine oceanic p Line 14  C--  o BLDEPTH     - Determine oceanic p
14  C--  o WSCALE      - Compute turbulent velocity scales.  C--  o WSCALE      - Compute turbulent velocity scales.
15  C--  o RI_IWMIX    - Compute interior viscosity diffusivity coefficients.  C--  o RI_IWMIX    - Compute interior viscosity diffusivity coefficients.
16  C--  o Z121        - Apply 121 vertical smoothing.  C--  o Z121        - Apply 121 vertical smoothing.
17  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.  
18  C--  o BLMIX       - Boundary layer mixing coefficients.  C--  o BLMIX       - Boundary layer mixing coefficients.
19  C--  o ENHANCE     - Enhance diffusivity at boundary layer interface.  C--  o ENHANCE     - Enhance diffusivity at boundary layer interface.
20  C--  o STATEKPP    - Compute buoyancy-related input arrays.  C--  o STATEKPP    - Compute buoyancy-related input arrays.
21    C--  o KPP_DOUBLEDIFF - Compute double diffusive contribution to diffusivities
22    
23  c***********************************************************************  c***********************************************************************
24    
25        SUBROUTINE KPPMIX (        SUBROUTINE KPPMIX (
26       I       mytime, mythid       I       kmtj, shsq, dvsq, ustar, msk
27       I     , kmtj, shsq, dvsq, ustar       I     , bo, bosol
28       I     , bo, bosol, dbloc, Ritop, coriol  #ifdef ALLOW_SALT_PLUME
29       I     , ikey       I     , boplume,SPDepth
30    #ifdef SALT_PLUME_SPLIT_BASIN
31         I     , lon,lat
32    #endif /* SALT_PLUME_SPLIT_BASIN */
33    #endif /* ALLOW_SALT_PLUME */
34         I     , dbloc, Ritop, coriol
35         I     , diffusKzS, diffusKzT
36         I     , ikppkey
37       O     , diffus       O     , diffus
38       U     , ghat       U     , ghat
39       O     , hbl )       O     , hbl
40         I     , bi, bj, myTime, myIter, myThid )
41    
42  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
43  c  c
# Line 47  c--------------------------------------- Line 58  c---------------------------------------
58  #include "SIZE.h"  #include "SIZE.h"
59  #include "EEPARAMS.h"  #include "EEPARAMS.h"
60  #include "PARAMS.h"  #include "PARAMS.h"
 #include "DYNVARS.h"  
 #include "FFIELDS.h"  
61  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
62    #ifdef ALLOW_AUTODIFF
63    # include "tamc.h"
64    #endif
65    
66  c input  c input
67  c     myTime          - current time in simulation  c   bi, bj :: Array indices on which to apply calculations
68  c     myThid          - thread number for this instance of the routine  c   myTime :: Current time in simulation
69  c     kmtj   (imt)    - number of vertical layers on this row  c   myIter :: Current iteration number in simulation
70  c     shsq   (imt,Nr) - (local velocity shear)^2                      ((m/s)^2)  c   myThid :: My Thread Id. number
71  c     dvsq   (imt,Nr) - (velocity shear re sfc)^2                     ((m/s)^2)  c     kmtj   (imt)     - number of vertical layers on this row
72  c     ustar  (imt)    - surface friction velocity                         (m/s)  c     msk    (imt)     - surface mask (=1 if water, =0 otherwise)
73  c     bo     (imt)    - surface turbulent buoy. forcing               (m^2/s^3)  c     shsq   (imt,Nr)  - (local velocity shear)^2                     ((m/s)^2)
74  c     bosol  (imt)    - radiative buoyancy forcing                    (m^2/s^3)  c     dvsq   (imt,Nr)  - (velocity shear re sfc)^2                    ((m/s)^2)
75  c     dbloc  (imt,Nr) - local delta buoyancy across interfaces          (m/s^2)  c     ustar  (imt)     - surface friction velocity                        (m/s)
76  c     dblocSm(imt,Nr) - horizontally smoothed dbloc                     (m/s^2)  c     bo     (imt)     - surface turbulent buoy. forcing              (m^2/s^3)
77  c                         stored in ghat to save space  c     bosol  (imt)     - radiative buoyancy forcing                   (m^2/s^3)
78  c     Ritop  (imt,Nr) - numerator of bulk Richardson Number  c     boplume(imt,Nrp1)- haline buoyancy forcing                      (m^2/s^3)
79  c                         (zref-z) * delta buoyancy w.r.t. surface    ((m/s)^2)  c     dbloc  (imt,Nr)  - local delta buoyancy across interfaces         (m/s^2)
80  c     coriol (imt)    - Coriolis parameter                                (1/s)  c     dblocSm(imt,Nr)  - horizontally smoothed dbloc                    (m/s^2)
81    c                          stored in ghat to save space
82    c     Ritop  (imt,Nr)  - numerator of bulk Richardson Number
83    c                          (zref-z) * delta buoyancy w.r.t. surface   ((m/s)^2)
84    c     coriol (imt)     - Coriolis parameter                               (1/s)
85    c     diffusKzS(imt,Nr)- background vertical diffusivity for scalars    (m^2/s)
86    c     diffusKzT(imt,Nr)- background vertical diffusivity for theta      (m^2/s)
87  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,
88  c           e.g., hbl(sNx,sNy) -> hbl(imt),  c           e.g., hbl(sNx,sNy) -> hbl(imt),
89  c           where hbl(i,j) -> hbl((j-1)*sNx+i)  c           where hbl(i,j) -> hbl((j-1)*sNx+i)
90          INTEGER bi, bj
91          _RL     myTime
92          integer myIter
93          integer myThid
94          integer kmtj (imt   )
95          _RL shsq     (imt,Nr)
96          _RL dvsq     (imt,Nr)
97          _RL ustar    (imt   )
98          _RL bo       (imt   )
99          _RL bosol    (imt   )
100    #ifdef ALLOW_SALT_PLUME
101          _RL boplume  (imt,Nrp1)
102          _RL SPDepth  (imt   )
103    #ifdef SALT_PLUME_SPLIT_BASIN
104          _RL lon  (imt   )
105          _RL lat  (imt   )
106    #endif /* SALT_PLUME_SPLIT_BASIN */
107    #endif /* ALLOW_SALT_PLUME */
108          _RL dbloc    (imt,Nr)
109          _RL Ritop    (imt,Nr)
110          _RL coriol   (imt   )
111          _RS msk      (imt   )
112          _RL diffusKzS(imt,Nr)
113          _RL diffusKzT(imt,Nr)
114    
115        _RL     mytime        integer ikppkey
       integer mythid  
       integer kmtj  (imt   )  
       _KPP_RL shsq  (imt,Nr)  
       _KPP_RL dvsq  (imt,Nr)  
       _KPP_RL ustar (imt   )  
       _KPP_RL bo    (imt   )  
       _KPP_RL bosol (imt   )  
       _KPP_RL dbloc (imt,Nr)  
       _KPP_RL Ritop (imt,Nr)  
       _KPP_RL coriol(imt   )  
   
       integer ikey  
116    
117  c output  c output
118  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 121  c     diffus (imt,3)  - vertical tempera
121  c     ghat   (imt)    - nonlocal transport coefficient                  (s/m^2)  c     ghat   (imt)    - nonlocal transport coefficient                  (s/m^2)
122  c     hbl    (imt)    - mixing layer depth                                  (m)  c     hbl    (imt)    - mixing layer depth                                  (m)
123    
124        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL diffus(imt,0:Nrp1,mdiff)
125        _KPP_RL ghat  (imt,Nr)        _RL ghat  (imt,Nr)
126        _KPP_RL hbl   (imt)        _RL hbl   (imt)
127    
128  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
129    
# Line 107  c     blmc   (imt,Nr,mdiff) - boundary l Line 137  c     blmc   (imt,Nr,mdiff) - boundary l
137  c     sigma  (imt         ) - normalized depth (d / hbl)  c     sigma  (imt         ) - normalized depth (d / hbl)
138  c     Rib    (imt,Nr      ) - bulk Richardson number  c     Rib    (imt,Nr      ) - bulk Richardson number
139    
140        integer kbl   (imt         )        integer kbl(imt         )
141        _KPP_RL bfsfc (imt         )        _RL bfsfc  (imt         )
142        _KPP_RL casea (imt         )        _RL casea  (imt         )
143        _KPP_RL stable(imt         )        _RL stable (imt         )
144        _KPP_RL dkm1  (imt,   mdiff)        _RL dkm1   (imt,   mdiff)
145        _KPP_RL blmc  (imt,Nr,mdiff)        _RL blmc   (imt,Nr,mdiff)
146        _KPP_RL sigma (imt         )        _RL sigma  (imt         )
147        _KPP_RL Rib   (imt,Nr      )        _RL Rib    (imt,Nr      )
148    
149        integer i, k, md        integer i, k, md
150    
# Line 125  c instability. Line 155  c instability.
155  c (ghat is temporary storage for horizontally smoothed dbloc)  c (ghat is temporary storage for horizontally smoothed dbloc)
156  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
157    
158  CADJ STORE ghat = comlev1_kpp, key = ikey  cph(
159    cph these storings avoid recomp. of Ri_iwmix
160    CADJ STORE ghat  = comlev1_kpp, key=ikppkey, kind=isbyte
161    CADJ STORE dbloc = comlev1_kpp, key=ikppkey, kind=isbyte
162    cph)
163        call Ri_iwmix (        call Ri_iwmix (
164       I       kmtj, shsq, dbloc, ghat       I       kmtj, shsq, dbloc, ghat
165       I     , ikey       I     , diffusKzS, diffusKzT
166       O     , diffus )       I     , ikppkey
167         O     , diffus, myThid )
168    
169    cph(
170    cph these storings avoid recomp. of Ri_iwmix
171    cph DESPITE TAFs 'not necessary' warning!
172    CADJ STORE dbloc  = comlev1_kpp, key=ikppkey, kind=isbyte
173    CADJ STORE shsq   = comlev1_kpp, key=ikppkey, kind=isbyte
174    CADJ STORE ghat   = comlev1_kpp, key=ikppkey, kind=isbyte
175    CADJ STORE diffus = comlev1_kpp, key=ikppkey, kind=isbyte
176    cph)
177    
178  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
179  c set seafloor values to zero and fill extra "Nrp1" coefficients  c set seafloor values to zero and fill extra "Nrp1" coefficients
# Line 138  c for blmix Line 181  c for blmix
181  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
182    
183        do md = 1, mdiff        do md = 1, mdiff
184           do k=1,Nrp1
185           do i = 1,imt           do i = 1,imt
186              do k=kmtj(i),Nrp1               if(k.ge.kmtj(i))  diffus(i,k,md) = 0.0
                diffus(i,k,md) = 0.0  
187              end do              end do
188           end do           end do
189        end do        end do
# Line 152  c diagnose the new boundary layer depth Line 195  c diagnose the new boundary layer depth
195  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
196    
197        call bldepth (        call bldepth (
198       I       mytime, mythid       I       kmtj
199       I     , kmtj       I     , dvsq, dbloc, Ritop, ustar, bo, bosol
200       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol  #ifdef ALLOW_SALT_PLUME
201       I     , ikey       I     , boplume,SPDepth
202    #ifdef SALT_PLUME_SPLIT_BASIN
203         I     , lon,lat
204    #endif /* SALT_PLUME_SPLIT_BASIN */
205    #endif /* ALLOW_SALT_PLUME */
206         I     , coriol
207         I     , ikppkey
208       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma
209       &     )       I     , bi, bj, myTime, myIter, myThid )
210    
211  CADJ STORE hbl,bfsfc,stable,casea,kbl = comlev1_kpp, key = ikey  CADJ STORE hbl,bfsfc,stable,casea,kbl = comlev1_kpp,
212    CADJ &     key=ikppkey, kind=isbyte
213    
214  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
215  c compute boundary layer diffusivities  c compute boundary layer diffusivities
# Line 167  c--------------------------------------- Line 217  c---------------------------------------
217    
218        call blmix (        call blmix (
219       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl
220       O     , dkm1, blmc, ghat, sigma, ikey       O     , dkm1, blmc, ghat, sigma, ikppkey
221       &     )       I     , myThid )
222    cph(
223  CADJ STORE dkm1,blmc,ghat = comlev1_kpp, key = ikey  CADJ STORE dkm1,blmc,ghat = comlev1_kpp,
224    CADJ &     key=ikppkey, kind=isbyte
225    CADJ STORE hbl, kbl, diffus, casea = comlev1_kpp,
226    CADJ &     key=ikppkey, kind=isbyte
227    cph)
228    
229  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
230  c enhance diffusivity at interface kbl - 1  c enhance diffusivity at interface kbl - 1
# Line 179  c--------------------------------------- Line 233  c---------------------------------------
233        call enhance (        call enhance (
234       I       dkm1, hbl, kbl, diffus, casea       I       dkm1, hbl, kbl, diffus, casea
235       U     , ghat       U     , ghat
236       O     , blmc )       O     , blmc
237         I     , myThid )
238    
239    cph(
240    cph avoids recomp. of enhance
241    CADJ STORE blmc = comlev1_kpp, key=ikppkey, kind=isbyte
242    cph)
243    
244  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
245  c combine interior and boundary layer coefficients and nonlocal term  c combine interior and boundary layer coefficients and nonlocal term
246    c !!!NOTE!!! In shallow (2-level) regions and for shallow mixed layers
247    c (< 1 level), diffusivity blmc can become negative.  The max-s below
248    c are a hack until this problem is properly diagnosed and fixed.
249  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
   
250        do k = 1, Nr        do k = 1, Nr
251           do i = 1, imt           do i = 1, imt
252              if (k .lt. kbl(i)) then              if (k .lt. kbl(i)) then
253                 do md = 1, mdiff  #ifdef ALLOW_SHELFICE
254                    diffus(i,k,md) = blmc(i,k,md)  C     when there is shelfice on top (msk(i)=0), reset the boundary layer
255                 end do  C     mixing coefficients blmc to pure Ri-number based mixing
256                   blmc(i,k,1) = max ( blmc(i,k,1)*msk(i),
257         &              diffus(i,k,1) )
258                   blmc(i,k,2) = max ( blmc(i,k,2)*msk(i),
259         &              diffus(i,k,2) )
260                   blmc(i,k,3) = max ( blmc(i,k,3)*msk(i),
261         &              diffus(i,k,3) )
262    #endif /* not ALLOW_SHELFICE */
263                   diffus(i,k,1) = max ( blmc(i,k,1), viscArNr(1) )
264                   diffus(i,k,2) = max ( blmc(i,k,2), diffusKzS(i,Nr) )
265                   diffus(i,k,3) = max ( blmc(i,k,3), diffusKzT(i,Nr) )
266              else              else
267                 ghat(i,k) = 0.                 ghat(i,k) = 0. _d 0
268              endif              endif
269           end do           end do
270        end do        end do
271    
272  #endif /* ALLOW_KPP */  #endif /* ALLOW_KPP */
273    
274        return        return
275        end        end
276    
277  c*************************************************************************  c*************************************************************************
278    
279        subroutine bldepth (        subroutine bldepth (
280       I       mytime, mythid       I       kmtj
281       I     , kmtj       I     , dvsq, dbloc, Ritop, ustar, bo, bosol
282       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol  #ifdef ALLOW_SALT_PLUME
283       I     , ikey       I     , boplume,SPDepth
284    #ifdef SALT_PLUME_SPLIT_BASIN
285         I     , lon,lat
286    #endif /* SALT_PLUME_SPLIT_BASIN */
287    #endif /* ALLOW_SALT_PLUME */
288         I     , coriol
289         I     , ikppkey
290       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma
291       &     )       I     , bi, bj, myTime, myIter, myThid )
292    
293  c     the oceanic planetary boundary layer depth, hbl, is determined as  c     the oceanic planetary boundary layer depth, hbl, is determined as
294  c     the shallowest depth where the bulk Richardson number is  c     the shallowest depth where the bulk Richardson number is
# Line 232  c     The water column and the surface f Line 310  c     The water column and the surface f
310  c     stable/ustable forcing conditions, and where hbl is relative  c     stable/ustable forcing conditions, and where hbl is relative
311  c     to grid points (caseA), so that conditional branches can be  c     to grid points (caseA), so that conditional branches can be
312  c     avoided in later subroutines.  c     avoided in later subroutines.
313  c  c
314        IMPLICIT NONE        IMPLICIT NONE
315    
316  #include "SIZE.h"  #include "SIZE.h"
317  #include "EEPARAMS.h"  #include "EEPARAMS.h"
318  #include "PARAMS.h"  #include "PARAMS.h"
319  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
320  #include "FFIELDS.h"  #ifdef ALLOW_AUTODIFF
321    # include "tamc.h"
322    #endif
323    
324  c input  c input
325  c------  c------
326  c myTime    : current time in simulation  c   bi, bj :: Array indices on which to apply calculations
327  c myThid    : thread number for this instance of the routine  c   myTime :: Current time in simulation
328    c   myIter :: Current iteration number in simulation
329    c   myThid :: My Thread Id. number
330  c kmtj      : number of vertical layers  c kmtj      : number of vertical layers
331  c dvsq      : (velocity shear re sfc)^2             ((m/s)^2)  c dvsq      : (velocity shear re sfc)^2             ((m/s)^2)
332  c dbloc     : local delta buoyancy across interfaces  (m/s^2)  c dbloc     : local delta buoyancy across interfaces  (m/s^2)
# Line 254  c             buoyancy with respect to s Line 336  c             buoyancy with respect to s
336  c ustar     : surface friction velocity                 (m/s)  c ustar     : surface friction velocity                 (m/s)
337  c bo        : surface turbulent buoyancy forcing    (m^2/s^3)  c bo        : surface turbulent buoyancy forcing    (m^2/s^3)
338  c bosol     : radiative buoyancy forcing            (m^2/s^3)  c bosol     : radiative buoyancy forcing            (m^2/s^3)
339    c boplume   : haline buoyancy forcing               (m^2/s^3)
340  c coriol    : Coriolis parameter                        (1/s)  c coriol    : Coriolis parameter                        (1/s)
341        _RL     mytime        INTEGER bi, bj
342        integer mythid        _RL     myTime
343          integer myIter
344          integer myThid
345        integer kmtj(imt)        integer kmtj(imt)
346        _KPP_RL dvsq  (imt,Nr)        _RL dvsq    (imt,Nr)
347        _KPP_RL dbloc (imt,Nr)        _RL dbloc   (imt,Nr)
348        _KPP_RL Ritop (imt,Nr)        _RL Ritop   (imt,Nr)
349        _KPP_RL ustar (imt)        _RL ustar   (imt)
350        _KPP_RL bo    (imt)        _RL bo      (imt)
351        _KPP_RL bosol (imt)        _RL bosol   (imt)
352        _KPP_RL coriol(imt)        _RL coriol  (imt)
353        integer ikey        integer ikppkey
354    #ifdef ALLOW_SALT_PLUME
355          _RL boplume (imt,Nrp1)
356          _RL SPDepth (imt)
357    #ifdef SALT_PLUME_SPLIT_BASIN
358          _RL lon (imt)
359          _RL lat (imt)
360    #endif /* SALT_PLUME_SPLIT_BASIN */
361    #endif /* ALLOW_SALT_PLUME */
362    
363  c  output  c  output
364  c--------  c--------
# Line 276  c casea     : =1 in case A, =0 in case B Line 369  c casea     : =1 in case A, =0 in case B
369  c kbl       : -1 of first grid level below hbl  c kbl       : -1 of first grid level below hbl
370  c Rib       : Bulk Richardson number  c Rib       : Bulk Richardson number
371  c sigma     : normalized depth (d/hbl)  c sigma     : normalized depth (d/hbl)
372        _KPP_RL hbl   (imt)        _RL hbl    (imt)
373        _KPP_RL bfsfc (imt)        _RL bfsfc  (imt)
374        _KPP_RL stable(imt)        _RL stable (imt)
375        _KPP_RL casea (imt)        _RL casea  (imt)
376        integer kbl   (imt)        integer kbl(imt)
377        _KPP_RL Rib   (imt,Nr)        _RL Rib    (imt,Nr)
378        _KPP_RL sigma (imt)        _RL sigma  (imt)
379    
380  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
381    
382  c  local  c  local
383  c-------  c-------
384  c wm, ws    : turbulent velocity scales         (m/s)  c wm, ws    : turbulent velocity scales         (m/s)
385        _KPP_RL wm(imt), ws(imt)        _RL wm(imt), ws(imt)
386        _RL     worka(imt)        _RL worka(imt)
387          _RL bvsq, vtsq, hekman, hmonob, hlimit, tempVar1, tempVar2
388        _KPP_RL bvsq, vtsq, hekman, hmonob, hlimit, tempVar1, tempVar2        integer i, k, kl, km, km1
389        integer i, kl        _RL temp
390    
391        _KPP_RL     p5    , eins        _RL         p5    , eins
392        parameter ( p5=0.5, eins=1.0 )        parameter ( p5=0.5, eins=1.0 )
393        _RL         minusone        _RL         minusone
394        parameter ( minusone=-1.0 )        parameter ( minusone=-1.0 )
395    #ifdef ALLOW_AUTODIFF_TAMC
396          integer kkppkey
397    #endif
398    
399    #ifdef ALLOW_DIAGNOSTICS
400    c     KPPBFSFC - Bo+radiation absorbed to d=hbf*hbl + plume (m^2/s^3)
401          _RL KPPBFSFC(imt,Nr)
402          _RL KPPRi(imt,Nr)
403    #endif /* ALLOW_DIAGNOSTICS */
404    
405  c find bulk Richardson number at every grid level until > Ricr  c find bulk Richardson number at every grid level until > Ricr
406  c  c
# Line 312  c       kbl(i)=kmtj(i) and hbl(i)=-zgrid Line 414  c       kbl(i)=kmtj(i) and hbl(i)=-zgrid
414  c     initialize hbl and kbl to bottomed out values  c     initialize hbl and kbl to bottomed out values
415    
416        do i = 1, imt        do i = 1, imt
417           Rib(i,1) = 0.0           Rib(i,1) = 0. _d 0
418           kbl(i) = max(kmtj(i),1)           kbl(i) = max(kmtj(i),1)
419           hbl(i) = -zgrid(kbl(i))           hbl(i) = -zgrid(kbl(i))
420        end do        end do
421    
422    #ifdef ALLOW_DIAGNOSTICS
423          do kl = 1, Nr
424             do i = 1, imt
425                KPPBFSFC(i,kl) = 0. _d 0
426                KPPRi(i,kl) = 0. _d 0
427             enddo
428          enddo
429    #endif /* ALLOW_DIAGNOSTICS */
430    
431        do kl = 2, Nr        do kl = 2, Nr
432    
433    #ifdef ALLOW_AUTODIFF_TAMC
434             kkppkey = (ikppkey-1)*Nr + kl
435    #endif
436    
437  c     compute bfsfc = sw fraction at hbf * zgrid  c     compute bfsfc = sw fraction at hbf * zgrid
438    
439           do i = 1, imt           do i = 1, imt
440              worka(i) = zgrid(kl)              worka(i) = zgrid(kl)
441           end do           end do
442    CADJ store worka = comlev1_kpp_k, key = kkppkey, kind=isbyte
443           call SWFRAC(           call SWFRAC(
444       I        imt, hbf,       I       imt, hbf,
445       I        mytime, mythid,       U       worka,
446       U        worka )       I       myTime, myIter, myThid )
447    CADJ store worka = comlev1_kpp_k, key = kkppkey, kind=isbyte
448    
449           do i = 1, imt           do i = 1, imt
450    
# Line 338  c     use caseA as temporary array Line 455  c     use caseA as temporary array
455  c     compute bfsfc= Bo + radiative contribution down to hbf * hbl  c     compute bfsfc= Bo + radiative contribution down to hbf * hbl
456    
457              bfsfc(i) = bo(i) + bosol(i)*(1. - worka(i))              bfsfc(i) = bo(i) + bosol(i)*(1. - worka(i))
             stable(i) = p5 + sign(p5,bfsfc(i))  
             sigma(i) = stable(i) + (1. - stable(i)) * epsilon  
458    
459           end do           end do
460    #ifdef ALLOW_SALT_PLUME
461    c     compute bfsfc = plume fraction at hbf * zgrid
462             IF ( useSALT_PLUME ) THEN
463    #ifndef SALT_PLUME_VOLUME
464               do i = 1, imt
465                  worka(i) = zgrid(kl)
466               enddo
467    Ccatn: in original way: accumulate all fractions of boplume above zgrid(kl)
468               call SALT_PLUME_FRAC(
469         I         imt, hbf,SPDepth,
470    #ifdef SALT_PLUME_SPLIT_BASIN
471         I         lon,lat,
472    #endif /* SALT_PLUME_SPLIT_BASIN */
473         U         worka,
474         I         myTime, myIter, myThid)
475               do i = 1, imt
476                  bfsfc(i) = bfsfc(i) + boplume(i,1)*(worka(i))
477    C            km=max(1,kbl(i)-1)
478    C            temp = (plumefrac(i,km)+plumefrac(i,kbl(i)))/2.0
479    C            bfsfc(i) = bfsfc(i) + boplume(i,1)*temp
480               enddo
481    #else /* def SALT_PLUME_VOLUME */
482    catn: in vol way: need to integrate down to hbl, so first locate
483    c     k level associated with this hbl, then sum up all SPforc[T,S]
484               DO i = 1, imt
485                km =max(1,kbl(i)-1)
486                km1=max(1,kbl(i))
487                temp = (boplume(i,km)+boplume(i,km1))/2.0
488                bfsfc(i) = bfsfc(i) + temp
489               ENDDO
490    #endif /* ndef SALT_PLUME_VOLUME */
491             ENDIF
492    #endif /* ALLOW_SALT_PLUME */
493    
494    #ifdef ALLOW_DIAGNOSTICS
495             do i = 1, imt
496                KPPBFSFC(i,kl) = bfsfc(i)
497             enddo
498    #endif /* ALLOW_DIAGNOSTICS */
499    
500             do i = 1, imt
501                stable(i) = p5 + sign(p5,bfsfc(i))
502                sigma(i) = stable(i) + (1. - stable(i)) * epsilon
503             enddo
504    
505  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
506  c     compute velocity scales at sigma, for hbl= caseA = -zgrid(kl)  c     compute velocity scales at sigma, for hbl= caseA = -zgrid(kl)
507  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
508    
509           call wscale (           call wscale (
510       I        sigma, casea, ustar, bfsfc,       I        sigma, casea, ustar, bfsfc,
511       O        wm, ws )       O        wm, ws, myThid )
512    CADJ store ws = comlev1_kpp_k, key = kkppkey, kind=isbyte
513    
514           do i = 1, imt           do i = 1, imt
515    
# Line 361  c--------------------------------------- Line 521  c---------------------------------------
521       1           ( dbloc(i,kl-1) / (zgrid(kl-1)-zgrid(kl  ))+       1           ( dbloc(i,kl-1) / (zgrid(kl-1)-zgrid(kl  ))+
522       2             dbloc(i,kl  ) / (zgrid(kl  )-zgrid(kl+1)))       2             dbloc(i,kl  ) / (zgrid(kl  )-zgrid(kl+1)))
523    
524              if (bvsq .eq. 0.) then              if (bvsq .eq. 0. _d 0) then
525                vtsq = 0.0                vtsq = 0. _d 0
526              else              else
527                vtsq = -zgrid(kl) * ws(i) * sqrt(abs(bvsq)) * Vtc                vtsq = -zgrid(kl) * ws(i) * sqrt(abs(bvsq)) * Vtc
528              endif              endif
# Line 379  c     rg: assignment to double precision Line 539  c     rg: assignment to double precision
539  c     ph: test for zero nominator  c     ph: test for zero nominator
540  c  c
541    
542              tempVar1  = dvsq(i,kl) + vtsq                        tempVar1  = dvsq(i,kl) + vtsq
543              tempVar2 = max(tempVar1, phepsi)              tempVar2 = max(tempVar1, phepsi)
544              Rib(i,kl) = Ritop(i,kl) / tempVar2              Rib(i,kl) = Ritop(i,kl) / tempVar2
545    #ifdef ALLOW_DIAGNOSTICS
546                KPPRi(i,kl) = Rib(i,kl)
547    #endif
548    
549           end do           end do
550        end do        end do
551          
552    #ifdef ALLOW_DIAGNOSTICS
553          IF ( useDiagnostics ) THEN
554             CALL DIAGNOSTICS_FILL(KPPBFSFC,'KPPbfsfc',0,Nr,2,bi,bj,myThid)
555             CALL DIAGNOSTICS_FILL(KPPRi   ,'KPPRi   ',0,Nr,2,bi,bj,myThid)
556          ENDIF
557    #endif /* ALLOW_DIAGNOSTICS */
558    
559    cph(
560    cph  without this store, there is a recomputation error for
561    cph  rib in adbldepth (probably partial recomputation problem)
562    CADJ store Rib = comlev1_kpp
563    CADJ &   , key=ikppkey, kind=isbyte,
564    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy),Nr /)
565    cph)
566    
567        do kl = 2, Nr        do kl = 2, Nr
568           do i = 1, imt           do i = 1, imt
569              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 571  c
571        end do        end do
572    
573  CADJ store kbl = comlev1_kpp  CADJ store kbl = comlev1_kpp
574  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key=ikppkey, kind=isbyte,
575    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
576    
577        do i = 1, imt        do i = 1, imt
578           kl = kbl(i)           kl = kbl(i)
# Line 406  c     linearly interpolate to find hbl w Line 585  c     linearly interpolate to find hbl w
585        end do        end do
586    
587  CADJ store hbl = comlev1_kpp  CADJ store hbl = comlev1_kpp
588  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key=ikppkey, kind=isbyte,
589    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
590    
591  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
592  c     find stability and buoyancy forcing for boundary layer  c     find stability and buoyancy forcing for boundary layer
# Line 415  c--------------------------------------- Line 595  c---------------------------------------
595        do i = 1, imt        do i = 1, imt
596           worka(i) = hbl(i)           worka(i) = hbl(i)
597        end do        end do
598    CADJ store worka = comlev1_kpp
599    CADJ &   , key=ikppkey, kind=isbyte,
600    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
601        call SWFRAC(        call SWFRAC(
602       I     imt, minusone,       I       imt, minusone,
603       I     mytime, mythid,       U       worka,
604       U     worka )       I       myTime, myIter, myThid )
605    CADJ store worka = comlev1_kpp
606    CADJ &   , key=ikppkey, kind=isbyte,
607    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
608    
609        do i = 1, imt        do i = 1, imt
610           bfsfc(i)  = bo(i) + bosol(i) * (1. - worka(i))           bfsfc(i)  = bo(i) + bosol(i) * (1. - worka(i))
611        end do        end do
612    
613    #ifdef ALLOW_SALT_PLUME
614          IF ( useSALT_PLUME ) THEN
615    #ifndef SALT_PLUME_VOLUME
616            do i = 1, imt
617               worka(i) = hbl(i)
618            enddo
619            call SALT_PLUME_FRAC(
620         I         imt,minusone,SPDepth,
621    #ifdef SALT_PLUME_SPLIT_BASIN
622         I         lon,lat,
623    #endif /* SALT_PLUME_SPLIT_BASIN */
624         U         worka,
625         I         myTime, myIter, myThid )
626            do i = 1, imt
627               bfsfc(i) = bfsfc(i) + boplume(i,1) * (worka(i))
628    C            km=max(1,kbl(i)-1)
629    C            temp = (plumefrac(i,km)+plumefrac(i,kbl(i)))/2.0
630    C            bfsfc(i) = bfsfc(i) + boplume(i,1)*temp
631            enddo
632    #else /* def SALT_PLUME_VOLUME */
633            DO i = 1, imt
634                km =max(1,kbl(i)-1)
635                km1=max(1,kbl(i))
636                temp = (boplume(i,km)+boplume(i,km1))/2.0
637                bfsfc(i) = bfsfc(i) + temp
638            ENDDO
639    #endif /* ndef SALT_PLUME_VOLUME */
640          ENDIF
641    #endif /* ALLOW_SALT_PLUME */
642  CADJ store bfsfc = comlev1_kpp  CADJ store bfsfc = comlev1_kpp
643  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key=ikppkey, kind=isbyte,
644    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
645    
646  c--   ensure bfsfc is never 0  c--   ensure bfsfc is never 0
647        do i = 1, imt        do i = 1, imt
# Line 432  c--   ensure bfsfc is never 0 Line 649  c--   ensure bfsfc is never 0
649           bfsfc(i) = sign(eins,bfsfc(i))*max(phepsi,abs(bfsfc(i)))           bfsfc(i) = sign(eins,bfsfc(i))*max(phepsi,abs(bfsfc(i)))
650        end do        end do
651    
652  CADJ store bfsfc = comlev1_kpp  cph(
653  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  cph  added stable to store list to avoid extensive recomp.
654    CADJ store bfsfc, stable = comlev1_kpp
655    CADJ &   , key=ikppkey, kind=isbyte,
656    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
657    cph)
658    
659  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
660  c check hbl limits for hekman or hmonob  c check hbl limits for hekman or hmonob
661  c     ph: test for zero nominator  c     ph: test for zero nominator
662  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
663    
664          IF ( LimitHblStable ) THEN
665        do i = 1, imt        do i = 1, imt
666           if (bfsfc(i) .gt. 0.0) then           if (bfsfc(i) .gt. 0.0) then
667              hekman = cekman * ustar(i) / max(abs(Coriol(i)),phepsi)              hekman = cekman * ustar(i) / max(abs(Coriol(i)),phepsi)
# Line 450  c--------------------------------------- Line 672  c---------------------------------------
672              hbl(i) = min(hbl(i),hlimit)              hbl(i) = min(hbl(i),hlimit)
673           end if           end if
674        end do        end do
675          ENDIF
676    
677  CADJ store hbl = comlev1_kpp  CADJ store hbl = comlev1_kpp
678  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key=ikppkey, kind=isbyte,
679    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
680    
681        do i = 1, imt        do i = 1, imt
682           hbl(i) = max(hbl(i),minKPPhbl)           hbl(i) = max(hbl(i),minKPPhbl)
# Line 459  CADJ &   , key = ikey, shape = (/ (sNx+2 Line 684  CADJ &   , key = ikey, shape = (/ (sNx+2
684        end do        end do
685    
686  CADJ store hbl = comlev1_kpp  CADJ store hbl = comlev1_kpp
687  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key=ikppkey, kind=isbyte,
688    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
689    
690  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
691  c      find new kbl  c      find new kbl
# Line 480  c--------------------------------------- Line 706  c---------------------------------------
706        do i = 1, imt        do i = 1, imt
707           worka(i) = hbl(i)           worka(i) = hbl(i)
708        end do        end do
709    CADJ store worka = comlev1_kpp
710    CADJ &   , key=ikppkey, kind=isbyte,
711    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
712        call SWFRAC(        call SWFRAC(
713       I     imt, minusone,       I     imt, minusone,
714       I     mytime, mythid,       U     worka,
715       U     worka )       I     myTime, myIter, myThid )
716    CADJ store worka = comlev1_kpp
717    CADJ &   , key=ikppkey, kind=isbyte,
718    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
719    
720        do i = 1, imt        do i = 1, imt
721           bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i))           bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i))
722        end do        end do
723    
724    #ifdef ALLOW_SALT_PLUME
725          IF ( useSALT_PLUME ) THEN
726    #ifndef SALT_PLUME_VOLUME
727            do i = 1, imt
728               worka(i) = hbl(i)
729            enddo
730            call SALT_PLUME_FRAC(
731         I         imt,minusone,SPDepth,
732    #ifdef SALT_PLUME_SPLIT_BASIN
733         I         lon,lat,
734    #endif /* SALT_PLUME_SPLIT_BASIN */
735         U         worka,
736         I         myTime, myIter, myThid )
737            do i = 1, imt
738               bfsfc(i) = bfsfc(i) + boplume(i,1) * (worka(i))
739    C            km=max(1,kbl(i)-1)
740    C            temp = (plumefrac(i,km)+plumefrac(i,kbl(i)))/2.0
741    C            bfsfc(i) = bfsfc(i) + boplume(i,1)*temp
742            enddo
743    #else /* def SALT_PLUME_VOLUME */
744            DO i = 1, imt
745                km =max(1,kbl(i)-1)
746                km1=max(1,kbl(i)-0)
747                temp = (boplume(i,km)+boplume(i,km1))/2.0
748                bfsfc(i) = bfsfc(i) + temp
749            ENDDO
750    #endif /* ndef SALT_PLUME_VOLUME */
751          ENDIF
752    #endif /* ALLOW_SALT_PLUME */
753  CADJ store bfsfc = comlev1_kpp  CADJ store bfsfc = comlev1_kpp
754  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key=ikppkey, kind=isbyte,
755    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
756    
757  c--   ensures bfsfc is never 0  c--   ensures bfsfc is never 0
758        do i = 1, imt        do i = 1, imt
# Line 515  c*************************************** Line 778  c***************************************
778    
779        subroutine wscale (        subroutine wscale (
780       I     sigma, hbl, ustar, bfsfc,       I     sigma, hbl, ustar, bfsfc,
781       O     wm, ws )       O     wm, ws,
782         I     myThid )
783    
784  c     compute turbulent velocity scales.  c     compute turbulent velocity scales.
785  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 524  c Line 788  c
788  c     note: the lookup table is only used for unstable conditions  c     note: the lookup table is only used for unstable conditions
789  c     (zehat.le.0), in the stable domain wm (=ws) gets computed  c     (zehat.le.0), in the stable domain wm (=ws) gets computed
790  c     directly.  c     directly.
791  c  c
792        IMPLICIT NONE        IMPLICIT NONE
793    
794  #include "SIZE.h"  #include "SIZE.h"
# Line 536  c sigma   : normalized depth (d/hbl) Line 800  c sigma   : normalized depth (d/hbl)
800  c hbl     : boundary layer depth (m)  c hbl     : boundary layer depth (m)
801  c ustar   : surface friction velocity         (m/s)  c ustar   : surface friction velocity         (m/s)
802  c bfsfc   : total surface buoyancy flux       (m^2/s^3)  c bfsfc   : total surface buoyancy flux       (m^2/s^3)
803        _KPP_RL sigma(imt)  c myThid  : thread number for this instance of the routine
804        _KPP_RL hbl  (imt)        integer myThid
805        _KPP_RL ustar(imt)        _RL sigma(imt)
806        _KPP_RL bfsfc(imt)        _RL hbl  (imt)
807          _RL ustar(imt)
808          _RL bfsfc(imt)
809    
810  c  output  c  output
811  c--------  c--------
812  c wm, ws  : turbulent velocity scales at sigma  c wm, ws  : turbulent velocity scales at sigma
813        _KPP_RL wm(imt), ws(imt)        _RL wm(imt), ws(imt)
814    
815  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
816    
817  c local  c local
818  c------  c------
819  c zehat   : = zeta *  ustar**3  c zehat   : = zeta *  ustar**3
820        _KPP_RL zehat        _RL zehat
821    
822        integer iz, izp1, ju, i, jup1        integer iz, izp1, ju, i, jup1
823        _KPP_RL udiff, zdiff, zfrac, ufrac, fzfrac, wam        _RL udiff, zdiff, zfrac, ufrac, fzfrac, wam
824        _KPP_RL wbm, was, wbs, u3, tempVar        _RL wbm, was, wbs, u3, tempVar
825    
826  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
827  c use lookup table for zehat < zmax only; otherwise use  c use lookup table for zehat < zmax only; otherwise use
# Line 599  c--------------------------------------- Line 865  c---------------------------------------
865              ws(i) = wm(i)              ws(i) = wm(i)
866    
867           endif           endif
868    
869        end do        end do
870    
871  #endif /* ALLOW_KPP */  #endif /* ALLOW_KPP */
872    
873        return        return
874        end        end
875    
876  c*************************************************************************  c*************************************************************************
877    
878        subroutine Ri_iwmix (        subroutine Ri_iwmix (
879       I       kmtj, shsq, dbloc, dblocSm       I       kmtj, shsq, dbloc, dblocSm,
880       I     , ikey       I       diffusKzS, diffusKzT,
881       O     , diffus )       I       ikppkey,
882         O       diffus,
883         I       myThid )
884    
885  c     compute interior viscosity diffusivity coefficients due  c     compute interior viscosity diffusivity coefficients due
886  c     to shear instability (dependent on a local Richardson number),  c     to shear instability (dependent on a local Richardson number),
# Line 625  c     to static instability (local Richa Line 893  c     to static instability (local Richa
893  #include "EEPARAMS.h"  #include "EEPARAMS.h"
894  #include "PARAMS.h"  #include "PARAMS.h"
895  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
896    #ifdef ALLOW_AUTODIFF
897    # include "AUTODIFF_PARAMS.h"
898    # include "tamc.h"
899    #endif
900    
901  c  input  c  input
902  c     kmtj   (imt)         number of vertical layers on this row  c     kmtj   (imt)         number of vertical layers on this row
903  c     shsq   (imt,Nr)      (local velocity shear)^2               ((m/s)^2)  c     shsq   (imt,Nr)      (local velocity shear)^2               ((m/s)^2)
904  c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)  c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)
905  c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)  c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)
906        integer kmtj   (imt)  c     diffusKzS(imt,Nr)- background vertical diffusivity for scalars    (m^2/s)
907        _KPP_RL shsq   (imt,Nr)  c     diffusKzT(imt,Nr)- background vertical diffusivity for theta      (m^2/s)
908        _KPP_RL dbloc  (imt,Nr)  c     myThid :: My Thread Id. number
909        _KPP_RL dblocSm(imt,Nr)        integer kmtj (imt)
910        integer ikey        _RL shsq     (imt,Nr)
911          _RL dbloc    (imt,Nr)
912          _RL dblocSm  (imt,Nr)
913          _RL diffusKzS(imt,Nr)
914          _RL diffusKzT(imt,Nr)
915          integer ikppkey
916          integer myThid
917    
918  c output  c output
919  c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)  c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)
920  c     diffus(imt,0:Nrp1,2)  vertical scalar diffusivity             (m^2/s)  c     diffus(imt,0:Nrp1,2)  vertical scalar diffusivity             (m^2/s)
921  c     diffus(imt,0:Nrp1,3)  vertical temperature diffusivity        (m^2/s)  c     diffus(imt,0:Nrp1,3)  vertical temperature diffusivity        (m^2/s)
922        _KPP_RL diffus(imt,0:Nrp1,3)        _RL diffus(imt,0:Nrp1,3)
923    
924  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
925    
926  c local variables  c local variables
927  c     Rig                   local Richardson number  c     Rig                   local Richardson number
928  c     fRi, fcon             function of Rig  c     fRi, fcon             function of Rig
929        _KPP_RL Rig        _RL Rig
930        _KPP_RL fRi, fcon        _RL fRi, fcon
931        _KPP_RL ratio        _RL ratio
932        integer i, ki, mr        integer i, ki, kp1
933        _KPP_RL c1, c0        _RL c1, c0
934    
935  #ifdef ALLOW_KPP_VERTICALLY_SMOOTH  #ifdef ALLOW_KPP_VERTICALLY_SMOOTH
936          integer mr
937  CADJ INIT kpp_ri_tape_mr = common, 1  CADJ INIT kpp_ri_tape_mr = common, 1
938  #endif  #endif
939    
940  c constants  c constants
941        c1 = 1.0        c1 = 1. _d 0
942        c0 = 0.0        c0 = 0. _d 0
943    
944  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
945  c     compute interior gradient Ri at all interfaces ki=1,Nr, (not surface)  c     compute interior gradient Ri at all interfaces ki=1,Nr, (not surface)
# Line 670  c     set values at bottom and below to Line 949  c     set values at bottom and below to
949  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
950  C     break data flow dependence on diffus  C     break data flow dependence on diffus
951        diffus(1,1,1) = 0.0        diffus(1,1,1) = 0.0
952    
953          do ki = 1, Nr
954             do i = 1, imt
955                diffus(i,ki,1) = 0.
956                diffus(i,ki,2) = 0.
957                diffus(i,ki,3) = 0.
958             enddo
959          enddo
960  #endif  #endif
961    
962        do ki = 1, Nr        do ki = 1, Nr
963           do i = 1, imt           do i = 1, imt
964              if     (kmtj(i) .EQ. 0      ) then              if     (kmtj(i) .LE. 1      ) then
965                 diffus(i,ki,1) = 0.                 diffus(i,ki,1) = 0.
966                 diffus(i,ki,2) = 0.                 diffus(i,ki,2) = 0.
967              elseif (ki      .GE. kmtj(i)) then              elseif (ki      .GE. kmtj(i)) then
# Line 687  C     break data flow dependence on diff Line 974  C     break data flow dependence on diff
974              endif              endif
975           end do           end do
976        end do        end do
977    CADJ store diffus = comlev1_kpp, key=ikppkey, kind=isbyte
978    
979  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
980  c     vertically smooth Ri  c     vertically smooth Ri
# Line 697  CADJ store diffus(:,:,1) = kpp_ri_tape_m Line 985  CADJ store diffus(:,:,1) = kpp_ri_tape_m
985  CADJ &  , key=mr, shape=(/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2 /)  CADJ &  , key=mr, shape=(/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2 /)
986    
987           call z121 (           call z121 (
988       U     diffus(1,0,1))       U     diffus(1,0,1),
989         I     myThid )
990        end do        end do
991  #endif  #endif
992    
# Line 706  c                           after smooth Line 995  c                           after smooth
995    
996        do ki = 1, Nr        do ki = 1, Nr
997           do i = 1, imt           do i = 1, imt
998    
999  c  evaluate f of Brunt-Vaisala squared for convection, store in fcon  c  evaluate f of Brunt-Vaisala squared for convection, store in fcon
1000    
1001              Rig   = max ( diffus(i,ki,2) , BVSQcon )              Rig   = max ( diffus(i,ki,2) , BVSQcon )
1002              ratio = min ( (BVSQcon - Rig) / BVSQcon, c1 )              ratio = min ( (BVSQcon - Rig) / BVSQcon, c1 )
1003              fcon  = c1 - ratio * ratio              fcon  = c1 - ratio * ratio
1004              fcon  = fcon * fcon * fcon              fcon  = fcon * fcon * fcon
1005    
1006  c  evaluate f of smooth Ri for shear instability, store in fRi  c  evaluate f of smooth Ri for shear instability, store in fRi
1007    
1008              Rig  = max ( diffus(i,ki,1), c0 )              Rig  = max ( diffus(i,ki,1), c0 )
1009              ratio = min ( Rig / Riinfty , c1 )              ratio = min ( Rig / Riinfty , c1 )
1010              fRi   = c1 - ratio * ratio              fRi   = c1 - ratio * ratio
1011              fRi   = fRi * fRi * fRi              fRi   = fRi * fRi * fRi
1012    
1013  c ----------------------------------------------------------------------  c ----------------------------------------------------------------------
1014  c            evaluate diffusivities and viscosity  c            evaluate diffusivities and viscosity
1015  c    mixing due to internal waves, and shear and static instability  c    mixing due to internal waves, and shear and static instability
   
             diffus(i,ki,1) = viscAr + fcon * difmcon + fRi * difm0  
             diffus(i,ki,2) = diffKrS + fcon * difscon + fRi * difs0  
             diffus(i,ki,3) = diffKrT + fcon * difscon + fRi * difs0  
1016    
1017                kp1 = MIN(ki+1,Nr)
1018    #ifdef EXCLUDE_KPP_SHEAR_MIX
1019                diffus(i,ki,1) = viscArNr(1)
1020                diffus(i,ki,2) = diffusKzS(i,kp1)
1021                diffus(i,ki,3) = diffusKzT(i,kp1)
1022    #else /* EXCLUDE_KPP_SHEAR_MIX */
1023    # ifdef ALLOW_AUTODIFF
1024                if ( inAdMode ) then
1025                  diffus(i,ki,1) = viscArNr(1)
1026                  diffus(i,ki,2) = diffusKzS(i,kp1)
1027                  diffus(i,ki,3) = diffusKzT(i,kp1)
1028                else
1029    # else /* ALLOW_AUTODIFF */
1030                if ( .TRUE. ) then
1031    # endif /* ALLOW_AUTODIFF */
1032                  diffus(i,ki,1) = viscArNr(1) + fcon*difmcon + fRi*difm0
1033                  diffus(i,ki,2) = diffusKzS(i,kp1)+fcon*difscon+fRi*difs0
1034                  diffus(i,ki,3) = diffusKzT(i,kp1)+fcon*diftcon+fRi*dift0
1035                endif
1036    #endif /* EXCLUDE_KPP_SHEAR_MIX */
1037           end do           end do
1038        end do        end do
1039    
1040  c ------------------------------------------------------------------------  c ------------------------------------------------------------------------
1041  c         set surface values to 0.0  c         set surface values to 0.0
1042    
1043        do i = 1, imt        do i = 1, imt
1044           diffus(i,0,1) = c0           diffus(i,0,1) = c0
1045           diffus(i,0,2) = c0           diffus(i,0,2) = c0
# Line 742  c         set surface values to 0.0 Line 1047  c         set surface values to 0.0
1047        end do        end do
1048    
1049  #endif /* ALLOW_KPP */  #endif /* ALLOW_KPP */
1050    
1051        return        return
1052        end        end
1053    
1054  c*************************************************************************  c*************************************************************************
1055    
1056        subroutine z121 (        subroutine z121 (
1057       U     v )       U     v,
1058         I     myThid )
1059    
1060  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)
1061  c     top (0) value is used as a dummy  c     top (0) value is used as a dummy
# Line 764  c     penetrative convection. Line 1070  c     penetrative convection.
1070  #include "SIZE.h"  #include "SIZE.h"
1071  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
1072    
1073  c input/output  c input/output
1074  c-------------  c-------------
1075  c v     : 2-D array to be smoothed in Nrp1 direction  c v     : 2-D array to be smoothed in Nrp1 direction
1076        _KPP_RL v(imt,0:Nrp1)  c myThid: thread number for this instance of the routine
1077          integer myThid
1078          _RL v(imt,0:Nrp1)
1079    
1080  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1081    
1082  c local  c local
1083        _KPP_RL zwork, zflag        _RL zwork, zflag
1084        _KPP_RL KRi_range(1:Nrp1)        _RL KRi_range(1:Nrp1)
1085        integer i, k, km1, kp1        integer i, k, km1, kp1
1086    
1087        _KPP_RL     p0      , p25       , p5      , p2        _RL         p0      , p25       , p5      , p2
1088        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 )
1089    
1090        KRi_range(Nrp1) = p0        KRi_range(Nrp1) = p0
# Line 830  CADJ STORE v(i,k), zwork = z121tape Line 1138  CADJ STORE v(i,k), zwork = z121tape
1138    
1139  c*************************************************************************  c*************************************************************************
1140    
       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*************************************************************************  
   
1141        subroutine smooth_horiz (        subroutine smooth_horiz (
1142       I     k, bi, bj,       I     k, bi, bj,
1143       U     fld )       U     fld,
1144         I     myThid )
1145    
1146  c     Apply horizontal smoothing to global _RL 2-D array  c     Apply horizontal smoothing to global _RL 2-D array
1147    
1148        IMPLICIT NONE        IMPLICIT NONE
1149  #include "SIZE.h"  #include "SIZE.h"
1150    #include "GRID.h"
1151  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
1152    
1153  c     input  c     input
1154  c     bi, bj : array indices  c     bi, bj : array indices
1155  c     k      : vertical index used for masking  c     k      : vertical index used for masking
1156    c     myThid : thread number for this instance of the routine
1157          INTEGER myThid
1158        integer k, bi, bj        integer k, bi, bj
1159    
1160  c     input/output  c     input/output
# Line 949  c     local Line 1181  c     local
1181              im1 = i-1              im1 = i-1
1182              ip1 = i+1              ip1 = i+1
1183              tempVar =              tempVar =
1184       &           p25   *   pMask(i  ,j  ,k,bi,bj)   +       &           p25   *   maskC(i  ,j  ,k,bi,bj)   +
1185       &           p125  * ( pMask(im1,j  ,k,bi,bj)   +       &           p125  * ( maskC(im1,j  ,k,bi,bj)   +
1186       &                     pMask(ip1,j  ,k,bi,bj)   +       &                     maskC(ip1,j  ,k,bi,bj)   +
1187       &                     pMask(i  ,jm1,k,bi,bj)   +       &                     maskC(i  ,jm1,k,bi,bj)   +
1188       &                     pMask(i  ,jp1,k,bi,bj) ) +       &                     maskC(i  ,jp1,k,bi,bj) ) +
1189       &           p0625 * ( pMask(im1,jm1,k,bi,bj)   +       &           p0625 * ( maskC(im1,jm1,k,bi,bj)   +
1190       &                     pMask(im1,jp1,k,bi,bj)   +       &                     maskC(im1,jp1,k,bi,bj)   +
1191       &                     pMask(ip1,jm1,k,bi,bj)   +       &                     maskC(ip1,jm1,k,bi,bj)   +
1192       &                     pMask(ip1,jp1,k,bi,bj) )       &                     maskC(ip1,jp1,k,bi,bj) )
1193              IF ( tempVar .GE. p25 ) THEN              IF ( tempVar .GE. p25 ) THEN
1194                 fld_tmp(i,j) = (                 fld_tmp(i,j) = (
1195       &              p25  * fld(i  ,j  )*pMask(i  ,j  ,k,bi,bj) +       &              p25  * fld(i  ,j  )*maskC(i  ,j  ,k,bi,bj) +
1196       &              p125 *(fld(im1,j  )*pMask(im1,j  ,k,bi,bj) +       &              p125 *(fld(im1,j  )*maskC(im1,j  ,k,bi,bj) +
1197       &                     fld(ip1,j  )*pMask(ip1,j  ,k,bi,bj) +       &                     fld(ip1,j  )*maskC(ip1,j  ,k,bi,bj) +
1198       &                     fld(i  ,jm1)*pMask(i  ,jm1,k,bi,bj) +       &                     fld(i  ,jm1)*maskC(i  ,jm1,k,bi,bj) +
1199       &                     fld(i  ,jp1)*pMask(i  ,jp1,k,bi,bj))+       &                     fld(i  ,jp1)*maskC(i  ,jp1,k,bi,bj))+
1200       &              p0625*(fld(im1,jm1)*pMask(im1,jm1,k,bi,bj) +       &              p0625*(fld(im1,jm1)*maskC(im1,jm1,k,bi,bj) +
1201       &                     fld(im1,jp1)*pMask(im1,jp1,k,bi,bj) +       &                     fld(im1,jp1)*maskC(im1,jp1,k,bi,bj) +
1202       &                     fld(ip1,jm1)*pMask(ip1,jm1,k,bi,bj) +       &                     fld(ip1,jm1)*maskC(ip1,jm1,k,bi,bj) +
1203       &                     fld(ip1,jp1)*pMask(ip1,jp1,k,bi,bj)))       &                     fld(ip1,jp1)*maskC(ip1,jp1,k,bi,bj)))
1204       &              / tempVar       &              / tempVar
1205              ELSE              ELSE
1206                 fld_tmp(i,j) = fld(i,j)                 fld_tmp(i,j) = fld(i,j)
# Line 992  c*************************************** Line 1224  c***************************************
1224    
1225        subroutine blmix (        subroutine blmix (
1226       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl
1227       O     , dkm1, blmc, ghat, sigma, ikey       O     , dkm1, blmc, ghat, sigma, ikppkey
1228       &     )       I     , myThid )
1229    
1230  c     mixing coefficients within boundary layer depend on surface  c     mixing coefficients within boundary layer depend on surface
1231  c     forcing and the magnitude and gradient of interior mixing below  c     forcing and the magnitude and gradient of interior mixing below
# Line 1007  c Line 1239  c
1239    
1240  #include "SIZE.h"  #include "SIZE.h"
1241  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
1242    #ifdef ALLOW_AUTODIFF
1243    # include "tamc.h"
1244    #endif
1245    
1246  c input  c input
1247  c     ustar (imt)                 surface friction velocity             (m/s)  c     ustar (imt)                 surface friction velocity             (m/s)
# Line 1015  c     hbl   (imt)                 bounda Line 1250  c     hbl   (imt)                 bounda
1250  c     stable(imt)                 = 1 in stable forcing  c     stable(imt)                 = 1 in stable forcing
1251  c     casea (imt)                 = 1 in case A  c     casea (imt)                 = 1 in case A
1252  c     diffus(imt,0:Nrp1,mdiff)    vertical diffusivities              (m^2/s)  c     diffus(imt,0:Nrp1,mdiff)    vertical diffusivities              (m^2/s)
1253  c     kbl(imt)                    -1 of first grid level below hbl  c     kbl   (imt)                 -1 of first grid level below hbl
1254        _KPP_RL ustar (imt)  c     myThid               thread number for this instance of the routine
1255        _KPP_RL bfsfc (imt)        integer myThid
1256        _KPP_RL hbl   (imt)        _RL ustar (imt)
1257        _KPP_RL stable(imt)        _RL bfsfc (imt)
1258        _KPP_RL casea (imt)        _RL hbl   (imt)
1259        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL stable(imt)
1260          _RL casea (imt)
1261          _RL diffus(imt,0:Nrp1,mdiff)
1262        integer kbl(imt)        integer kbl(imt)
1263    
1264  c output  c output
# Line 1029  c     dkm1 (imt,mdiff)            bounda Line 1266  c     dkm1 (imt,mdiff)            bounda
1266  c     blmc (imt,Nr,mdiff)         boundary layer mixing coefficients  (m^2/s)  c     blmc (imt,Nr,mdiff)         boundary layer mixing coefficients  (m^2/s)
1267  c     ghat (imt,Nr)               nonlocal scalar transport  c     ghat (imt,Nr)               nonlocal scalar transport
1268  c     sigma(imt)                  normalized depth (d / hbl)  c     sigma(imt)                  normalized depth (d / hbl)
1269        _KPP_RL dkm1 (imt,mdiff)        _RL dkm1 (imt,mdiff)
1270        _KPP_RL blmc (imt,Nr,mdiff)        _RL blmc (imt,Nr,mdiff)
1271        _KPP_RL ghat (imt,Nr)        _RL ghat (imt,Nr)
1272        _KPP_RL sigma(imt)        _RL sigma(imt)
1273        integer ikey        integer ikppkey
1274    
1275  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1276    
# Line 1041  c  local Line 1278  c  local
1278  c     gat1*(imt)                 shape function at sigma = 1  c     gat1*(imt)                 shape function at sigma = 1
1279  c     dat1*(imt)                 derivative of shape function at sigma = 1  c     dat1*(imt)                 derivative of shape function at sigma = 1
1280  c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)  c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)
1281        _KPP_RL gat1m(imt), gat1s(imt), gat1t(imt)        _RL gat1m(imt), gat1s(imt), gat1t(imt)
1282        _KPP_RL dat1m(imt), dat1s(imt), dat1t(imt)        _RL dat1m(imt), dat1s(imt), dat1t(imt)
1283        _KPP_RL ws(imt), wm(imt)        _RL ws(imt), wm(imt)
1284        integer i, kn, ki        integer i, kn, ki
1285        _KPP_RL R, dvdzup, dvdzdn, viscp        _RL R, dvdzup, dvdzdn, viscp
1286        _KPP_RL difsp, diftp, visch, difsh, difth        _RL difsp, diftp, visch, difsh, difth
1287        _KPP_RL f1, sig, a1, a2, a3, delhat        _RL f1, sig, a1, a2, a3, delhat
1288        _KPP_RL Gm, Gs, Gt        _RL Gm, Gs, Gt
1289        _KPP_RL tempVar        _RL tempVar
1290    
1291        _KPP_RL    p0    , eins        _RL    p0    , eins
1292        parameter (p0=0.0, eins=1.0)        parameter (p0=0.0, eins=1.0)
1293    #ifdef ALLOW_AUTODIFF_TAMC
1294          integer kkppkey
1295    #endif
1296    
1297  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1298  c compute velocity scales at hbl  c compute velocity scales at hbl
# Line 1062  c--------------------------------------- Line 1302  c---------------------------------------
1302           sigma(i) = stable(i) * 1.0 + (1. - stable(i)) * epsilon           sigma(i) = stable(i) * 1.0 + (1. - stable(i)) * epsilon
1303        end do        end do
1304    
1305    CADJ STORE sigma = comlev1_kpp, key=ikppkey, kind=isbyte
1306        call wscale (        call wscale (
1307       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1308       O        wm, ws )       O        wm, ws, myThid )
1309    CADJ STORE wm = comlev1_kpp, key=ikppkey, kind=isbyte
1310    CADJ STORE ws = comlev1_kpp, key=ikppkey, kind=isbyte
1311    
1312        do i = 1, imt        do i = 1, imt
1313           wm(i) = sign(eins,wm(i))*max(phepsi,abs(wm(i)))           wm(i) = sign(eins,wm(i))*max(phepsi,abs(wm(i)))
1314           ws(i) = sign(eins,ws(i))*max(phepsi,abs(ws(i)))           ws(i) = sign(eins,ws(i))*max(phepsi,abs(ws(i)))
1315        end do        end do
1316  CADJ STORE wm = comlev1_kpp, key = ikey  CADJ STORE wm = comlev1_kpp, key=ikppkey, kind=isbyte
1317  CADJ STORE ws = comlev1_kpp, key = ikey  CADJ STORE ws = comlev1_kpp, key=ikppkey, kind=isbyte
1318    
1319        do i = 1, imt        do i = 1, imt
1320    
# Line 1103  c--------------------------------------- Line 1346  c---------------------------------------
1346           difsh  = diffus(i,kn,2) + difsp * delhat           difsh  = diffus(i,kn,2) + difsp * delhat
1347           difth  = diffus(i,kn,3) + diftp * delhat           difth  = diffus(i,kn,3) + diftp * delhat
1348    
1349           f1 = stable(i) * conc1 * bfsfc(i) /           f1 = stable(i) * conc1 * bfsfc(i) /
1350       &        max(ustar(i)**4,phepsi)       &        max(ustar(i)**4,phepsi)
1351           gat1m(i) = visch / hbl(i) / wm(i)           gat1m(i) = visch / hbl(i) / wm(i)
1352           dat1m(i) = -viscp / wm(i) + f1 * visch           dat1m(i) = -viscp / wm(i) + f1 * visch
1353           dat1m(i) = min(dat1m(i),p0)  
   
1354           gat1s(i) = difsh  / hbl(i) / ws(i)           gat1s(i) = difsh  / hbl(i) / ws(i)
1355           dat1s(i) = -difsp / ws(i) + f1 * difsh           dat1s(i) = -difsp / ws(i) + f1 * difsh
1356           dat1s(i) = min(dat1s(i),p0)  
   
1357           gat1t(i) = difth /  hbl(i) / ws(i)           gat1t(i) = difth /  hbl(i) / ws(i)
1358           dat1t(i) = -diftp / ws(i) + f1 * difth           dat1t(i) = -diftp / ws(i) + f1 * difth
1359    
1360          end do
1361    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1362    CADJ STORE gat1m = comlev1_kpp, key=ikppkey, kind=isbyte
1363    CADJ STORE gat1s = comlev1_kpp, key=ikppkey, kind=isbyte
1364    CADJ STORE gat1t = comlev1_kpp, key=ikppkey, kind=isbyte
1365    CADJ STORE dat1m = comlev1_kpp, key=ikppkey, kind=isbyte
1366    CADJ STORE dat1s = comlev1_kpp, key=ikppkey, kind=isbyte
1367    CADJ STORE dat1t = comlev1_kpp, key=ikppkey, kind=isbyte
1368    #endif
1369          do i = 1, imt
1370             dat1m(i) = min(dat1m(i),p0)
1371             dat1s(i) = min(dat1s(i),p0)
1372           dat1t(i) = min(dat1t(i),p0)           dat1t(i) = min(dat1t(i),p0)
   
1373        end do        end do
1374    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1375    CADJ STORE dat1m = comlev1_kpp, key=ikppkey, kind=isbyte
1376    CADJ STORE dat1s = comlev1_kpp, key=ikppkey, kind=isbyte
1377    CADJ STORE dat1t = comlev1_kpp, key=ikppkey, kind=isbyte
1378    #endif
1379    
1380        do ki = 1, Nr        do ki = 1, Nr
1381    
1382    #ifdef ALLOW_AUTODIFF_TAMC
1383             kkppkey = (ikppkey-1)*Nr + ki
1384    #endif
1385    
1386  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1387  c     compute turbulent velocity scales on the interfaces  c     compute turbulent velocity scales on the interfaces
1388  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 1129  c--------------------------------------- Line 1391  c---------------------------------------
1391              sig      = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i)              sig      = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i)
1392              sigma(i) = stable(i)*sig + (1.-stable(i))*min(sig,epsilon)              sigma(i) = stable(i)*sig + (1.-stable(i))*min(sig,epsilon)
1393           end do           end do
1394    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1395    CADJ STORE wm = comlev1_kpp_k, key = kkppkey
1396    CADJ STORE ws = comlev1_kpp_k, key = kkppkey
1397    #endif
1398    CADJ STORE sigma = comlev1_kpp_k, key = kkppkey
1399           call wscale (           call wscale (
1400       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1401       O        wm, ws )       O        wm, ws, myThid )
1402    CADJ STORE wm = comlev1_kpp_k, key = kkppkey
1403    CADJ STORE ws = comlev1_kpp_k, key = kkppkey
1404    
1405  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1406  c     compute the dimensionless shape functions at the interfaces  c     compute the dimensionless shape functions at the interfaces
# Line 1161  c--------------------------------------- Line 1430  c---------------------------------------
1430    
1431              tempVar = ws(i) * hbl(i)              tempVar = ws(i) * hbl(i)
1432              ghat(i,ki) = (1.-stable(i)) * cg / max(phepsi,tempVar)              ghat(i,ki) = (1.-stable(i)) * cg / max(phepsi,tempVar)
1433    
1434           end do           end do
1435        end do        end do
1436    
# Line 1171  c--------------------------------------- Line 1440  c---------------------------------------
1440    
1441        do i = 1, imt        do i = 1, imt
1442           sig      = -zgrid(kbl(i)-1) / hbl(i)           sig      = -zgrid(kbl(i)-1) / hbl(i)
1443           sigma(i) = stable(i) * sig           sigma(i) = stable(i) * sig
1444       &            + (1. - stable(i)) * min(sig,epsilon)       &            + (1. - stable(i)) * min(sig,epsilon)
1445        end do        end do
1446    
1447    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1448    CADJ STORE wm = comlev1_kpp, key=ikppkey, kind=isbyte
1449    CADJ STORE ws = comlev1_kpp, key=ikppkey, kind=isbyte
1450    #endif
1451    CADJ STORE sigma = comlev1_kpp, key=ikppkey, kind=isbyte
1452        call wscale (        call wscale (
1453       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1454       O        wm, ws )       O        wm, ws, myThid )
1455    CADJ STORE wm = comlev1_kpp, key=ikppkey, kind=isbyte
1456    CADJ STORE ws = comlev1_kpp, key=ikppkey, kind=isbyte
1457    
1458        do i = 1, imt        do i = 1, imt
1459           sig = -zgrid(kbl(i)-1) / hbl(i)           sig = -zgrid(kbl(i)-1) / hbl(i)
# Line 1194  c--------------------------------------- Line 1470  c---------------------------------------
1470    
1471  #endif /* ALLOW_KPP */  #endif /* ALLOW_KPP */
1472    
1473        return        return
1474        end        end
1475    
1476  c*************************************************************************  c*************************************************************************
1477    
1478        subroutine enhance (        subroutine enhance (
1479       I       dkm1, hbl, kbl, diffus, casea       I       dkm1, hbl, kbl, diffus, casea
1480       U     , ghat       U     , ghat
1481       O     , blmc       O     , blmc
1482       &     )       &     , myThid )
1483    
1484  c enhance the diffusivity at the kbl-.5 interface  c enhance the diffusivity at the kbl-.5 interface
1485    
# Line 1218  c     hbl(imt)                  boundary Line 1494  c     hbl(imt)                  boundary
1494  c     kbl(imt)                  grid above hbl  c     kbl(imt)                  grid above hbl
1495  c     diffus(imt,0:Nrp1,mdiff) vertical diffusivities           (m^2/s)  c     diffus(imt,0:Nrp1,mdiff) vertical diffusivities           (m^2/s)
1496  c     casea(imt)                = 1 in caseA, = 0 in case B  c     casea(imt)                = 1 in caseA, = 0 in case B
1497        _KPP_RL dkm1  (imt,mdiff)  c     myThid                    thread number for this instance of the routine
1498        _KPP_RL hbl   (imt)        integer   myThid
1499          _RL dkm1  (imt,mdiff)
1500          _RL hbl   (imt)
1501        integer kbl   (imt)        integer kbl   (imt)
1502        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL diffus(imt,0:Nrp1,mdiff)
1503        _KPP_RL casea (imt)        _RL casea (imt)
1504    
1505  c input/output  c input/output
1506  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)
1507        _KPP_RL ghat (imt,Nr)        _RL ghat (imt,Nr)
1508    
1509  c output  c output
1510  c     enhanced bound. layer mixing coeff.  c     enhanced bound. layer mixing coeff.
1511        _KPP_RL blmc  (imt,Nr,mdiff)        _RL blmc  (imt,Nr,mdiff)
1512    
1513  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1514    
1515  c local  c local
1516  c     fraction hbl lies beteen zgrid neighbors  c     fraction hbl lies beteen zgrid neighbors
1517        _KPP_RL delta        _RL delta
1518        integer ki, i, md        integer ki, i, md
1519        _KPP_RL dkmp5, dstar        _RL dkmp5, dstar
1520    
1521        do i = 1, imt        do i = 1, imt
1522           ki = kbl(i)-1           ki = kbl(i)-1
# Line 1258  c     fraction hbl lies beteen zgrid nei Line 1536  c     fraction hbl lies beteen zgrid nei
1536    
1537  #endif /* ALLOW_KPP */  #endif /* ALLOW_KPP */
1538    
1539        return        return
1540        end        end
1541    
1542  c*************************************************************************  c*************************************************************************
1543    
1544        SUBROUTINE STATEKPP (        SUBROUTINE STATEKPP (
1545       I     bi, bj, myThid,       O     RHO1, DBLOC, DBSFC, TTALPHA, SSBETA,
1546       O     RHO1, DBLOC, DBSFC, TTALPHA, SSBETA)       I     ikppkey, bi, bj, myThid )
1547  c  c
1548  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1549  c     "statekpp" computes all necessary input arrays  c     "statekpp" computes all necessary input arrays
# Line 1290  c Line 1568  c
1568  c     written  by: jan morzel,   feb. 10, 1995 (converted from "sigma" version)  c     written  by: jan morzel,   feb. 10, 1995 (converted from "sigma" version)
1569  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
1570  c  c
1571    
1572  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1573    
1574        IMPLICIT NONE        IMPLICIT NONE
# Line 1299  c--------------------------------------- Line 1578  c---------------------------------------
1578  #include "PARAMS.h"  #include "PARAMS.h"
1579  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
1580  #include "DYNVARS.h"  #include "DYNVARS.h"
1581    #include "GRID.h"
1582    #ifdef ALLOW_AUTODIFF
1583    # include "tamc.h"
1584    #endif
1585    
1586  c-------------- Routine arguments -----------------------------------------  c-------------- Routine arguments -----------------------------------------
1587        INTEGER bi, bj, myThid        INTEGER bi, bj, myThid
1588        _KPP_RL RHO1   ( ibot:itop, jbot:jtop       )        _RL RHO1   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy       )
1589        _KPP_RL DBLOC  ( ibot:itop, jbot:jtop, Nr   )        _RL DBLOC  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
1590        _KPP_RL DBSFC  ( ibot:itop, jbot:jtop, Nr   )        _RL DBSFC  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
1591        _KPP_RL TTALPHA( ibot:itop, jbot:jtop, Nrp1 )        _RL TTALPHA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1592        _KPP_RL SSBETA ( ibot:itop, jbot:jtop, Nrp1 )        _RL SSBETA ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1593    
1594  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1595    
# Line 1317  c Line 1600  c
1600  c     rhok         - density of t(k  ) & s(k  ) at depth k  c     rhok         - density of t(k  ) & s(k  ) at depth k
1601  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
1602  c     rho1k        - density of t(1  ) & s(1  ) at depth k  c     rho1k        - density of t(1  ) & s(1  ) at depth k
1603  c     work1, work2 - work arrays for holding horizontal slabs  c     work1,2,3    - work arrays for holding horizontal slabs
1604    
1605        _RL RHOK  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL RHOK  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1606        _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 1608  c     work1, work2 - work arrays for hol
1608        _RL WORK1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1609        _RL WORK2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1610        _RL WORK3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1611    
1612        INTEGER I, J, K        INTEGER I, J, K
1613          INTEGER ikppkey, kkppkey
1614    
1615  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
1616    
1617        call FIND_RHO(        kkppkey = (ikppkey-1)*Nr + 1
1618       I     bi, bj, ibot, itop, jbot, jtop, 1, 1, eosType,  
1619       I     theta, salt,  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1620    CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k,
1621    CADJ &     key=kkppkey, kind=isbyte
1622    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k,
1623    CADJ &     key=kkppkey, kind=isbyte
1624    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1625          CALL FIND_RHO_2D(
1626         I     1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1,
1627         I     theta(1-OLx,1-OLy,1,bi,bj), salt(1-OLx,1-OLy,1,bi,bj),
1628       O     WORK1,       O     WORK1,
1629       I     myThid )       I     1, bi, bj, myThid )
1630    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1631    CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k,
1632    CADJ &     key=kkppkey, kind=isbyte
1633    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k,
1634    CADJ &     key=kkppkey, kind=isbyte
1635    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1636    
1637        call FIND_ALPHA(        call FIND_ALPHA(
1638       I     bi, bj, ibot, itop, jbot, jtop, 1, 1, eosType,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1639       O     WORK2 )       O     WORK2, myThid )
1640    
1641        call FIND_BETA(        call FIND_BETA(
1642       I     bi, bj, ibot, itop, jbot, jtop, 1, 1, eosType,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1643       O     WORK3 )       O     WORK3, myThid )
1644    
1645        DO J = jbot, jtop        DO J = 1-OLy, sNy+OLy
1646           DO I = ibot, itop           DO I = 1-OLx, sNx+OLx
1647              RHO1(I,J)      = WORK1(I,J) + rhonil              RHO1(I,J)      = WORK1(I,J) + rhoConst
1648              TTALPHA(I,J,1) = WORK2(I,J)              TTALPHA(I,J,1) = WORK2(I,J)
1649              SSBETA(I,J,1)  = WORK3(I,J)              SSBETA(I,J,1)  = WORK3(I,J)
1650              DBSFC(I,J,1)   = 0.              DBSFC(I,J,1)   = 0.
# Line 1357  c calculate alpha, beta, and gradients i Line 1656  c calculate alpha, beta, and gradients i
1656  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)
1657        DO K = 2, Nr        DO K = 2, Nr
1658    
1659           call FIND_RHO(        kkppkey = (ikppkey-1)*Nr + k
1660       I        bi, bj, ibot, itop, jbot, jtop, K, K, eosType,  
1661       I        theta, salt,  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1662    CADJ STORE theta(:,:,k,bi,bj) = comlev1_kpp_k,
1663    CADJ &     key=kkppkey, kind=isbyte
1664    CADJ STORE salt (:,:,k,bi,bj) = comlev1_kpp_k,
1665    CADJ &     key=kkppkey, kind=isbyte
1666    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1667             CALL FIND_RHO_2D(
1668         I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, k,
1669         I        theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
1670       O        RHOK,       O        RHOK,
1671       I        myThid )       I        k, bi, bj, myThid )
1672    
1673           call FIND_RHO(  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1674       I        bi, bj, ibot, itop, jbot, jtop, K-1, K, eosType,  CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_kpp_k,
1675       I        theta, salt,  CADJ &     key=kkppkey, kind=isbyte
1676    CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_kpp_k,
1677    CADJ &     key=kkppkey, kind=isbyte
1678    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1679             CALL FIND_RHO_2D(
1680         I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, k,
1681         I        theta(1-OLx,1-OLy,k-1,bi,bj),salt(1-OLx,1-OLy,k-1,bi,bj),
1682       O        RHOKM1,       O        RHOKM1,
1683       I        myThid )       I        k-1, bi, bj, myThid )
1684    
1685           call FIND_RHO(  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1686       I        bi, bj, ibot, itop, jbot, jtop, 1, K, eosType,  CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k,
1687       I        theta, salt,  CADJ &     key=kkppkey, kind=isbyte
1688    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k,
1689    CADJ &     key=kkppkey, kind=isbyte
1690    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1691             CALL FIND_RHO_2D(
1692         I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, k,
1693         I        theta(1-OLx,1-OLy,1,bi,bj), salt(1-OLx,1-OLy,1,bi,bj),
1694       O        RHO1K,       O        RHO1K,
1695       I        myThid )       I        1, bi, bj, myThid )
1696    
1697    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1698    CADJ STORE rhok  (:,:)          = comlev1_kpp_k,
1699    CADJ &     key=kkppkey, kind=isbyte
1700    CADJ STORE rhokm1(:,:)          = comlev1_kpp_k,
1701    CADJ &     key=kkppkey, kind=isbyte
1702    CADJ STORE rho1k (:,:)          = comlev1_kpp_k,
1703    CADJ &     key=kkppkey, kind=isbyte
1704    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1705    
1706           call FIND_ALPHA(           call FIND_ALPHA(
1707       I        bi, bj, ibot, itop, jbot, jtop, K, K, eosType,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1708       O        WORK1 )       O        WORK1, myThid )
1709    
1710           call FIND_BETA(           call FIND_BETA(
1711       I        bi, bj, ibot, itop, jbot, jtop, K, K, eosType,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1712       O        WORK2 )       O        WORK2, myThid )
1713    
1714           DO J = jbot, jtop           DO J = 1-OLy, sNy+OLy
1715              DO I = ibot, itop              DO I = 1-OLx, sNx+OLx
1716                 TTALPHA(I,J,K) = WORK1 (I,J)                 TTALPHA(I,J,K) = WORK1 (I,J)
1717                 SSBETA(I,J,K)  = WORK2 (I,J)                 SSBETA(I,J,K)  = WORK2 (I,J)
1718                 DBLOC(I,J,K-1) = gravity * (RHOK(I,J) - RHOKM1(I,J)) /                 DBLOC(I,J,K-1) = gravity * (RHOK(I,J) - RHOKM1(I,J)) /
1719       &                                    (RHOK(I,J) + rhonil)       &                                    (RHOK(I,J) + rhoConst)
1720                 DBSFC(I,J,K)   = gravity * (RHOK(I,J) - RHO1K (I,J)) /                 DBSFC(I,J,K)   = gravity * (RHOK(I,J) - RHO1K (I,J)) /
1721       &                                    (RHOK(I,J) + rhonil)       &                                    (RHOK(I,J) + rhoConst)
1722              END DO              END DO
1723           END DO           END DO
1724    
1725        END DO        END DO
1726    
1727  c     compute arrays for K = Nrp1  c     compute arrays for K = Nrp1
1728        DO J = jbot, jtop        DO J = 1-OLy, sNy+OLy
1729           DO I = ibot, itop           DO I = 1-OLx, sNx+OLx
1730              TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)              TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)
1731              SSBETA(I,J,Nrp1)  = SSBETA(I,J,Nr)              SSBETA(I,J,Nrp1)  = SSBETA(I,J,Nr)
1732              DBLOC(I,J,Nr)     = 0.              DBLOC(I,J,Nr)     = 0.
1733           END DO           END DO
1734        END DO        END DO
1735    
1736    #ifdef ALLOW_DIAGNOSTICS
1737          IF ( useDiagnostics ) THEN
1738             CALL DIAGNOSTICS_FILL(DBSFC ,'KPPdbsfc',0,Nr,2,bi,bj,myThid)
1739             CALL DIAGNOSTICS_FILL(DBLOC ,'KPPdbloc',0,Nr,2,bi,bj,myThid)
1740          ENDIF
1741    #endif /* ALLOW_DIAGNOSTICS */
1742    
1743    #endif /* ALLOW_KPP */
1744    
1745          RETURN
1746          END
1747    
1748    c*************************************************************************
1749    
1750          SUBROUTINE KPP_DOUBLEDIFF (
1751         I     TTALPHA, SSBETA,
1752         U     kappaRT,
1753         U     kappaRS,
1754         I     ikppkey, imin, imax, jmin, jmax, bi, bj, myThid )
1755    c
1756    c-----------------------------------------------------------------------
1757    c     "KPP_DOUBLEDIFF" adds the double diffusive contributions
1758    C     as Rrho-dependent parameterizations to kappaRT and kappaRS
1759    c
1760    c     input:
1761    c     bi, bj  = array indices on which to apply calculations
1762    c     imin, imax, jmin, jmax = array boundaries
1763    c     ikppkey = key for TAMC/TAF automatic differentiation
1764    c     myThid  = thread id
1765    c
1766    c      ttalpha= thermal expansion coefficient without 1/rho factor
1767    c               d(rho) / d(potential temperature)                    (kg/m^3/C)
1768    c      ssbeta = salt expansion coefficient without 1/rho factor
1769    c               d(rho) / d(salinity)                               (kg/m^3/PSU)
1770    c     output: updated
1771    c     kappaRT/S :: background diffusivities for temperature and salinity
1772    c
1773    c     written  by: martin losch,   sept. 15, 2009
1774    c
1775    
1776    c-----------------------------------------------------------------------
1777    
1778          IMPLICIT NONE
1779    
1780    #include "SIZE.h"
1781    #include "EEPARAMS.h"
1782    #include "PARAMS.h"
1783    #include "KPP_PARAMS.h"
1784    #include "DYNVARS.h"
1785    #include "GRID.h"
1786    #ifdef ALLOW_AUTODIFF
1787    # include "tamc.h"
1788    #endif
1789    
1790    c-------------- Routine arguments -----------------------------------------
1791          INTEGER ikppkey, imin, imax, jmin, jmax, bi, bj, myThid
1792    
1793          _RL TTALPHA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1794          _RL SSBETA ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1795          _RL KappaRT( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
1796          _RL KappaRS( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
1797    
1798    #ifdef ALLOW_KPP
1799    
1800    C--------------------------------------------------------------------------
1801    C
1802    C     local variables
1803    C     I,J,K :: loop indices
1804    C     kkppkey :: key for TAMC/TAF automatic differentiation
1805    C
1806          INTEGER I, J, K
1807          INTEGER kkppkey
1808    C     alphaDT   :: d\rho/d\theta * d\theta
1809    C     betaDS    :: d\rho/dsalt * dsalt
1810    C     Rrho      :: "density ratio" R_{\rho} = \alpha dT/dz / \beta dS/dz
1811    C     nuddt/s   :: double diffusive diffusivities
1812    C     numol     :: molecular diffusivity
1813    C     rFac      :: abbreviation for 1/(R_{\rho0}-1)
1814    
1815          _RL alphaDT ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
1816          _RL betaDS  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
1817          _RL nuddt   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
1818          _RL nudds   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
1819          _RL Rrho
1820          _RL numol, rFac, nutmp
1821          INTEGER Km1
1822    
1823    C     set some constants here
1824          numol = 1.5 _d -06
1825          rFac  = 1. _d 0 / (Rrho0 - 1. _d 0 )
1826    C
1827          kkppkey = (ikppkey-1)*Nr + 1
1828    
1829    CML#ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1830    CMLCADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k,
1831    CMLCADJ &     key=kkppkey, kind=isbyte
1832    CMLCADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k,
1833    CMLCADJ &     key=kkppkey, kind=isbyte
1834    CML#endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1835    
1836          DO K = 1, Nr
1837           Km1 = MAX(K-1,1)
1838           DO J = 1-OLy, sNy+OLy
1839            DO I = 1-OLx, sNx+OLx
1840             alphaDT(I,J) = ( theta(I,J,Km1,bi,bj)-theta(I,J,K,bi,bj) )
1841         &        * 0.5 _d 0 * ABS( TTALPHA(I,J,Km1) + TTALPHA(I,J,K) )
1842             betaDS(I,J)  = ( salt(I,J,Km1,bi,bj)-salt(I,J,K,bi,bj) )
1843         &        * 0.5 _d 0 * ( SSBETA(I,J,Km1) + SSBETA(I,J,K) )
1844             nuddt(I,J) = 0. _d 0
1845             nudds(I,J) = 0. _d 0
1846            ENDDO
1847           ENDDO
1848           IF ( K .GT. 1 ) THEN
1849            DO J = jMin, jMax
1850             DO I = iMin, iMax
1851              Rrho  = 0. _d 0
1852    C     Now we have many different cases
1853    C     a. alphaDT > 0 and betaDS > 0 => salt fingering
1854    C        (salinity destabilizes)
1855              IF (      alphaDT(I,J) .GT. betaDS(I,J)
1856         &         .AND. betaDS(I,J) .GT. 0. _d 0 ) THEN
1857               Rrho = MIN( alphaDT(I,J)/betaDS(I,J), Rrho0 )
1858    C     Large et al. 1994, eq. 31a
1859    C          nudds(I,J) = dsfmax * ( 1. _d 0 - (Rrho - 1. _d 0) * rFac )**3
1860               nutmp      =          ( 1. _d 0 - (Rrho - 1. _d 0) * rFac )
1861               nudds(I,J) = dsfmax * nutmp * nutmp * nutmp
1862    C     Large et al. 1994, eq. 31c
1863               nuddt(I,J) = 0.7 _d 0 * nudds(I,J)
1864              ELSEIF (   alphaDT(I,J) .LT. 0. _d 0
1865         &          .AND. betaDS(I,J) .LT. 0. _d 0
1866         &          .AND.alphaDT(I,J) .GT. betaDS(I,J) ) THEN
1867    C     b. alphaDT < 0 and betaDS < 0 => semi-convection, diffusive convection
1868    C        (temperature destabilizes)
1869    C     for Rrho >= 1 the water column is statically unstable and we never
1870    C     reach this point
1871               Rrho = alphaDT(I,J)/betaDS(I,J)
1872    C     Large et al. 1994, eq. 32
1873               nuddt(I,J) = numol * 0.909 _d 0
1874         &          * exp ( 4.6 _d 0 * exp (
1875         &          - 5.4 _d 0 * ( 1. _d 0/Rrho - 1. _d 0 ) ) )
1876    CMLC     or
1877    CMLC     Large et al. 1994, eq. 33
1878    CML         nuddt(I,J) = numol * 8.7 _d 0 * Rrho**1.1
1879    C     Large et al. 1994, eqs. 34
1880               nudds(I,J) = nuddt(I,J) * MAX( 0.15 _d 0 * Rrho,
1881         &          1.85 _d 0 * Rrho - 0.85 _d 0 )
1882              ELSE
1883    C     Do nothing, because in this case the water colume is unstable
1884    C     => double diffusive processes are negligible and mixing due
1885    C     to shear instability will dominate
1886              ENDIF
1887             ENDDO
1888            ENDDO
1889    C     ENDIF ( K .GT. 1 )
1890           ENDIF
1891    C
1892           DO J = 1-OLy, sNy+OLy
1893            DO I = 1-OLx, sNx+OLx
1894             kappaRT(I,J,K) = kappaRT(I,J,K) + nuddt(I,J)
1895             kappaRS(I,J,K) = kappaRS(I,J,K) + nudds(I,J)
1896            ENDDO
1897           ENDDO
1898    #ifdef ALLOW_DIAGNOSTICS
1899           IF ( useDiagnostics ) THEN
1900            CALL DIAGNOSTICS_FILL(nuddt,'KPPnuddt',k,1,2,bi,bj,myThid)
1901            CALL DIAGNOSTICS_FILL(nudds,'KPPnudds',k,1,2,bi,bj,myThid)
1902           ENDIF
1903    #endif /* ALLOW_DIAGNOSTICS */
1904    C     end of K-loop
1905          ENDDO
1906  #endif /* ALLOW_KPP */  #endif /* ALLOW_KPP */
1907    
1908        RETURN        RETURN
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.52
  ViewVC Help
Powered by ViewVC 1.1.22