/[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.51 by atn, Mon Mar 31 20:47:32 2014 UTC
# Line 11  C--  o BLDEPTH     - Determine oceanic p Line 11  C--  o BLDEPTH     - Determine oceanic p
11  C--  o WSCALE      - Compute turbulent velocity scales.  C--  o WSCALE      - Compute turbulent velocity scales.
12  C--  o RI_IWMIX    - Compute interior viscosity diffusivity coefficients.  C--  o RI_IWMIX    - Compute interior viscosity diffusivity coefficients.
13  C--  o Z121        - Apply 121 vertical smoothing.  C--  o Z121        - Apply 121 vertical smoothing.
14  C--  o KPP_SMOOTH_HORIZ - Apply horizontal smoothing to KPP array.  C--  o SMOOTH_HORIZ- Apply horizontal smoothing to global array.
 C--  o SMOOTH_HORIZ     - Apply horizontal smoothing to global array.  
15  C--  o BLMIX       - Boundary layer mixing coefficients.  C--  o BLMIX       - Boundary layer mixing coefficients.
16  C--  o ENHANCE     - Enhance diffusivity at boundary layer interface.  C--  o ENHANCE     - Enhance diffusivity at boundary layer interface.
17  C--  o STATEKPP    - Compute buoyancy-related input arrays.  C--  o STATEKPP    - Compute buoyancy-related input arrays.
18    C--  o KPP_DOUBLEDIFF - Compute double diffusive contribution to diffusivities
19    
20  c***********************************************************************  c***********************************************************************
21    
22        SUBROUTINE KPPMIX (        SUBROUTINE KPPMIX (
23       I       mytime, mythid       I       kmtj, shsq, dvsq, ustar, msk
24       I     , kmtj, shsq, dvsq, ustar       I     , bo, bosol
25       I     , bo, bosol, dbloc, Ritop, coriol  #ifdef ALLOW_SALT_PLUME
26       I     , ikey       I     , boplume,SPDepth
27    #endif /* ALLOW_SALT_PLUME */
28         I     , dbloc, Ritop, coriol
29         I     , diffusKzS, diffusKzT
30         I     , ikppkey
31       O     , diffus       O     , diffus
32       U     , ghat       U     , ghat
33       O     , hbl )       O     , hbl
34         I     , bi, bj, myTime, myIter, myThid )
35    
36  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
37  c  c
# Line 47  c--------------------------------------- Line 52  c---------------------------------------
52  #include "SIZE.h"  #include "SIZE.h"
53  #include "EEPARAMS.h"  #include "EEPARAMS.h"
54  #include "PARAMS.h"  #include "PARAMS.h"
 #include "DYNVARS.h"  
 #include "FFIELDS.h"  
55  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
56    #ifdef ALLOW_AUTODIFF
57    # include "tamc.h"
58    #endif
59    
60  c input  c input
61  c     myTime          - current time in simulation  c   bi, bj :: Array indices on which to apply calculations
62  c     myThid          - thread number for this instance of the routine  c   myTime :: Current time in simulation
63  c     kmtj   (imt)    - number of vertical layers on this row  c   myIter :: Current iteration number in simulation
64  c     shsq   (imt,Nr) - (local velocity shear)^2                      ((m/s)^2)  c   myThid :: My Thread Id. number
65  c     dvsq   (imt,Nr) - (velocity shear re sfc)^2                     ((m/s)^2)  c     kmtj   (imt)     - number of vertical layers on this row
66  c     ustar  (imt)    - surface friction velocity                         (m/s)  c     msk    (imt)     - surface mask (=1 if water, =0 otherwise)
67  c     bo     (imt)    - surface turbulent buoy. forcing               (m^2/s^3)  c     shsq   (imt,Nr)  - (local velocity shear)^2                     ((m/s)^2)
68  c     bosol  (imt)    - radiative buoyancy forcing                    (m^2/s^3)  c     dvsq   (imt,Nr)  - (velocity shear re sfc)^2                    ((m/s)^2)
69  c     dbloc  (imt,Nr) - local delta buoyancy across interfaces          (m/s^2)  c     ustar  (imt)     - surface friction velocity                        (m/s)
70  c     dblocSm(imt,Nr) - horizontally smoothed dbloc                     (m/s^2)  c     bo     (imt)     - surface turbulent buoy. forcing              (m^2/s^3)
71  c                         stored in ghat to save space  c     bosol  (imt)     - radiative buoyancy forcing                   (m^2/s^3)
72  c     Ritop  (imt,Nr) - numerator of bulk Richardson Number  c     boplume(imt)     - haline buoyancy forcing                      (m^2/s^3)
73  c                         (zref-z) * delta buoyancy w.r.t. surface    ((m/s)^2)  c     dbloc  (imt,Nr)  - local delta buoyancy across interfaces         (m/s^2)
74  c     coriol (imt)    - Coriolis parameter                                (1/s)  c     dblocSm(imt,Nr)  - horizontally smoothed dbloc                    (m/s^2)
75    c                          stored in ghat to save space
76    c     Ritop  (imt,Nr)  - numerator of bulk Richardson Number
77    c                          (zref-z) * delta buoyancy w.r.t. surface   ((m/s)^2)
78    c     coriol (imt)     - Coriolis parameter                               (1/s)
79    c     diffusKzS(imt,Nr)- background vertical diffusivity for scalars    (m^2/s)
80    c     diffusKzT(imt,Nr)- background vertical diffusivity for theta      (m^2/s)
81  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,
82  c           e.g., hbl(sNx,sNy) -> hbl(imt),  c           e.g., hbl(sNx,sNy) -> hbl(imt),
83  c           where hbl(i,j) -> hbl((j-1)*sNx+i)  c           where hbl(i,j) -> hbl((j-1)*sNx+i)
84          INTEGER bi, bj
85          _RL     myTime
86          integer myIter
87          integer myThid
88          integer kmtj (imt   )
89          _RL shsq     (imt,Nr)
90          _RL dvsq     (imt,Nr)
91          _RL ustar    (imt   )
92          _RL bo       (imt   )
93          _RL bosol    (imt   )
94    #ifdef ALLOW_SALT_PLUME
95          _RL boplume  (imt   )
96          _RL SPDepth  (imt   )
97    #endif /* ALLOW_SALT_PLUME */
98          _RL dbloc    (imt,Nr)
99          _RL Ritop    (imt,Nr)
100          _RL coriol   (imt   )
101          _RS msk      (imt   )
102          _RL diffusKzS(imt,Nr)
103          _RL diffusKzT(imt,Nr)
104    
105        _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  
106    
107  c output  c output
108  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 111  c     diffus (imt,3)  - vertical tempera
111  c     ghat   (imt)    - nonlocal transport coefficient                  (s/m^2)  c     ghat   (imt)    - nonlocal transport coefficient                  (s/m^2)
112  c     hbl    (imt)    - mixing layer depth                                  (m)  c     hbl    (imt)    - mixing layer depth                                  (m)
113    
114        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL diffus(imt,0:Nrp1,mdiff)
115        _KPP_RL ghat  (imt,Nr)        _RL ghat  (imt,Nr)
116        _KPP_RL hbl   (imt)        _RL hbl   (imt)
117    
118  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
119    
# Line 107  c     blmc   (imt,Nr,mdiff) - boundary l Line 127  c     blmc   (imt,Nr,mdiff) - boundary l
127  c     sigma  (imt         ) - normalized depth (d / hbl)  c     sigma  (imt         ) - normalized depth (d / hbl)
128  c     Rib    (imt,Nr      ) - bulk Richardson number  c     Rib    (imt,Nr      ) - bulk Richardson number
129    
130        integer kbl   (imt         )        integer kbl(imt         )
131        _KPP_RL bfsfc (imt         )        _RL bfsfc  (imt         )
132        _KPP_RL casea (imt         )        _RL casea  (imt         )
133        _KPP_RL stable(imt         )        _RL stable (imt         )
134        _KPP_RL dkm1  (imt,   mdiff)        _RL dkm1   (imt,   mdiff)
135        _KPP_RL blmc  (imt,Nr,mdiff)        _RL blmc   (imt,Nr,mdiff)
136        _KPP_RL sigma (imt         )        _RL sigma  (imt         )
137        _KPP_RL Rib   (imt,Nr      )        _RL Rib    (imt,Nr      )
138    
139        integer i, k, md        integer i, k, md
140    
# Line 125  c instability. Line 145  c instability.
145  c (ghat is temporary storage for horizontally smoothed dbloc)  c (ghat is temporary storage for horizontally smoothed dbloc)
146  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
147    
148  CADJ STORE ghat = comlev1_kpp, key = ikey  cph(
149    cph these storings avoid recomp. of Ri_iwmix
150    CADJ STORE ghat  = comlev1_kpp, key=ikppkey, kind=isbyte
151    CADJ STORE dbloc = comlev1_kpp, key=ikppkey, kind=isbyte
152    cph)
153        call Ri_iwmix (        call Ri_iwmix (
154       I       kmtj, shsq, dbloc, ghat       I       kmtj, shsq, dbloc, ghat
155       I     , ikey       I     , diffusKzS, diffusKzT
156       O     , diffus )       I     , ikppkey
157         O     , diffus, myThid )
158    
159    cph(
160    cph these storings avoid recomp. of Ri_iwmix
161    cph DESPITE TAFs 'not necessary' warning!
162    CADJ STORE dbloc  = comlev1_kpp, key=ikppkey, kind=isbyte
163    CADJ STORE shsq   = comlev1_kpp, key=ikppkey, kind=isbyte
164    CADJ STORE ghat   = comlev1_kpp, key=ikppkey, kind=isbyte
165    CADJ STORE diffus = comlev1_kpp, key=ikppkey, kind=isbyte
166    cph)
167    
168  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
169  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 171  c for blmix
171  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
172    
173        do md = 1, mdiff        do md = 1, mdiff
174           do k=1,Nrp1
175           do i = 1,imt           do i = 1,imt
176              do k=kmtj(i),Nrp1               if(k.ge.kmtj(i))  diffus(i,k,md) = 0.0
                diffus(i,k,md) = 0.0  
177              end do              end do
178           end do           end do
179        end do        end do
# Line 152  c diagnose the new boundary layer depth Line 185  c diagnose the new boundary layer depth
185  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
186    
187        call bldepth (        call bldepth (
188       I       mytime, mythid       I       kmtj
189       I     , kmtj       I     , dvsq, dbloc, Ritop, ustar, bo, bosol
190       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol  #ifdef ALLOW_SALT_PLUME
191       I     , ikey       I     , boplume,SPDepth
192    #endif /* ALLOW_SALT_PLUME */
193         I     , coriol
194         I     , ikppkey
195       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma
196       &     )       I     , bi, bj, myTime, myIter, myThid )
197    
198  CADJ STORE hbl,bfsfc,stable,casea,kbl = comlev1_kpp, key = ikey  CADJ STORE hbl,bfsfc,stable,casea,kbl = comlev1_kpp,
199    CADJ &     key=ikppkey, kind=isbyte
200    
201  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
202  c compute boundary layer diffusivities  c compute boundary layer diffusivities
# Line 167  c--------------------------------------- Line 204  c---------------------------------------
204    
205        call blmix (        call blmix (
206       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl
207       O     , dkm1, blmc, ghat, sigma, ikey       O     , dkm1, blmc, ghat, sigma, ikppkey
208       &     )       I     , myThid )
209    cph(
210  CADJ STORE dkm1,blmc,ghat = comlev1_kpp, key = ikey  CADJ STORE dkm1,blmc,ghat = comlev1_kpp,
211    CADJ &     key=ikppkey, kind=isbyte
212    CADJ STORE hbl, kbl, diffus, casea = comlev1_kpp,
213    CADJ &     key=ikppkey, kind=isbyte
214    cph)
215    
216  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
217  c enhance diffusivity at interface kbl - 1  c enhance diffusivity at interface kbl - 1
# Line 179  c--------------------------------------- Line 220  c---------------------------------------
220        call enhance (        call enhance (
221       I       dkm1, hbl, kbl, diffus, casea       I       dkm1, hbl, kbl, diffus, casea
222       U     , ghat       U     , ghat
223       O     , blmc )       O     , blmc
224         I     , myThid )
225    
226    cph(
227    cph avoids recomp. of enhance
228    CADJ STORE blmc = comlev1_kpp, key=ikppkey, kind=isbyte
229    cph)
230    
231  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
232  c combine interior and boundary layer coefficients and nonlocal term  c combine interior and boundary layer coefficients and nonlocal term
233    c !!!NOTE!!! In shallow (2-level) regions and for shallow mixed layers
234    c (< 1 level), diffusivity blmc can become negative.  The max-s below
235    c are a hack until this problem is properly diagnosed and fixed.
236  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
   
237        do k = 1, Nr        do k = 1, Nr
238           do i = 1, imt           do i = 1, imt
239              if (k .lt. kbl(i)) then              if (k .lt. kbl(i)) then
240                 do md = 1, mdiff  #ifdef ALLOW_SHELFICE
241                    diffus(i,k,md) = blmc(i,k,md)  C     when there is shelfice on top (msk(i)=0), reset the boundary layer
242                 end do  C     mixing coefficients blmc to pure Ri-number based mixing
243                   blmc(i,k,1) = max ( blmc(i,k,1)*msk(i),
244         &              diffus(i,k,1) )
245                   blmc(i,k,2) = max ( blmc(i,k,2)*msk(i),
246         &              diffus(i,k,2) )
247                   blmc(i,k,3) = max ( blmc(i,k,3)*msk(i),
248         &              diffus(i,k,3) )
249    #endif /* not ALLOW_SHELFICE */
250                   diffus(i,k,1) = max ( blmc(i,k,1), viscArNr(1) )
251                   diffus(i,k,2) = max ( blmc(i,k,2), diffusKzS(i,Nr) )
252                   diffus(i,k,3) = max ( blmc(i,k,3), diffusKzT(i,Nr) )
253              else              else
254                 ghat(i,k) = 0.                 ghat(i,k) = 0. _d 0
255              endif              endif
256           end do           end do
257        end do        end do
258    
259  #endif /* ALLOW_KPP */  #endif /* ALLOW_KPP */
260    
261        return        return
262        end        end
263    
264  c*************************************************************************  c*************************************************************************
265    
266        subroutine bldepth (        subroutine bldepth (
267       I       mytime, mythid       I       kmtj
268       I     , kmtj       I     , dvsq, dbloc, Ritop, ustar, bo, bosol
269       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol  #ifdef ALLOW_SALT_PLUME
270       I     , ikey       I     , boplume,SPDepth
271    #endif /* ALLOW_SALT_PLUME */
272         I     , coriol
273         I     , ikppkey
274       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma
275       &     )       I     , bi, bj, myTime, myIter, myThid )
276    
277  c     the oceanic planetary boundary layer depth, hbl, is determined as  c     the oceanic planetary boundary layer depth, hbl, is determined as
278  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 294  c     The water column and the surface f
294  c     stable/ustable forcing conditions, and where hbl is relative  c     stable/ustable forcing conditions, and where hbl is relative
295  c     to grid points (caseA), so that conditional branches can be  c     to grid points (caseA), so that conditional branches can be
296  c     avoided in later subroutines.  c     avoided in later subroutines.
297  c  c
298        IMPLICIT NONE        IMPLICIT NONE
299    
300  #include "SIZE.h"  #include "SIZE.h"
301  #include "EEPARAMS.h"  #include "EEPARAMS.h"
302  #include "PARAMS.h"  #include "PARAMS.h"
303  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
304  #include "FFIELDS.h"  #ifdef ALLOW_AUTODIFF
305    # include "tamc.h"
306    #endif
307    
308  c input  c input
309  c------  c------
310  c myTime    : current time in simulation  c   bi, bj :: Array indices on which to apply calculations
311  c myThid    : thread number for this instance of the routine  c   myTime :: Current time in simulation
312    c   myIter :: Current iteration number in simulation
313    c   myThid :: My Thread Id. number
314  c kmtj      : number of vertical layers  c kmtj      : number of vertical layers
315  c dvsq      : (velocity shear re sfc)^2             ((m/s)^2)  c dvsq      : (velocity shear re sfc)^2             ((m/s)^2)
316  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 320  c             buoyancy with respect to s
320  c ustar     : surface friction velocity                 (m/s)  c ustar     : surface friction velocity                 (m/s)
321  c bo        : surface turbulent buoyancy forcing    (m^2/s^3)  c bo        : surface turbulent buoyancy forcing    (m^2/s^3)
322  c bosol     : radiative buoyancy forcing            (m^2/s^3)  c bosol     : radiative buoyancy forcing            (m^2/s^3)
323    c boplume   : haline buoyancy forcing               (m^2/s^3)
324  c coriol    : Coriolis parameter                        (1/s)  c coriol    : Coriolis parameter                        (1/s)
325        _RL     mytime        INTEGER bi, bj
326        integer mythid        _RL     myTime
327          integer myIter
328          integer myThid
329        integer kmtj(imt)        integer kmtj(imt)
330        _KPP_RL dvsq  (imt,Nr)        _RL dvsq    (imt,Nr)
331        _KPP_RL dbloc (imt,Nr)        _RL dbloc   (imt,Nr)
332        _KPP_RL Ritop (imt,Nr)        _RL Ritop   (imt,Nr)
333        _KPP_RL ustar (imt)        _RL ustar   (imt)
334        _KPP_RL bo    (imt)        _RL bo      (imt)
335        _KPP_RL bosol (imt)        _RL bosol   (imt)
336        _KPP_RL coriol(imt)        _RL coriol  (imt)
337        integer ikey        integer ikppkey
338    #ifdef ALLOW_SALT_PLUME
339          _RL boplume (imt)
340          _RL SPDepth (imt)
341    #endif /* ALLOW_SALT_PLUME */
342    
343  c  output  c  output
344  c--------  c--------
# Line 276  c casea     : =1 in case A, =0 in case B Line 349  c casea     : =1 in case A, =0 in case B
349  c kbl       : -1 of first grid level below hbl  c kbl       : -1 of first grid level below hbl
350  c Rib       : Bulk Richardson number  c Rib       : Bulk Richardson number
351  c sigma     : normalized depth (d/hbl)  c sigma     : normalized depth (d/hbl)
352        _KPP_RL hbl   (imt)        _RL hbl    (imt)
353        _KPP_RL bfsfc (imt)        _RL bfsfc  (imt)
354        _KPP_RL stable(imt)        _RL stable (imt)
355        _KPP_RL casea (imt)        _RL casea  (imt)
356        integer kbl   (imt)        integer kbl(imt)
357        _KPP_RL Rib   (imt,Nr)        _RL Rib    (imt,Nr)
358        _KPP_RL sigma (imt)        _RL sigma  (imt)
359    
360  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
361    
362  c  local  c  local
363  c-------  c-------
364  c wm, ws    : turbulent velocity scales         (m/s)  c wm, ws    : turbulent velocity scales         (m/s)
365        _KPP_RL wm(imt), ws(imt)        _RL wm(imt), ws(imt)
366        _RL     worka(imt)        _RL worka(imt)
367          _RL bvsq, vtsq, hekman, hmonob, hlimit, tempVar1, tempVar2
       _KPP_RL bvsq, vtsq, hekman, hmonob, hlimit, tempVar1, tempVar2  
368        integer i, kl        integer i, kl
369    
370        _KPP_RL     p5    , eins        _RL         p5    , eins
371        parameter ( p5=0.5, eins=1.0 )        parameter ( p5=0.5, eins=1.0 )
372        _RL         minusone        _RL         minusone
373        parameter ( minusone=-1.0 )        parameter ( minusone=-1.0 )
374    #ifdef ALLOW_AUTODIFF_TAMC
375          integer kkppkey
376    #endif
377    
378    #ifdef ALLOW_DIAGNOSTICS
379    c     KPPBFSFC - Bo+radiation absorbed to d=hbf*hbl + plume (m^2/s^3)
380          _RL KPPBFSFC(imt,Nr)
381          _RL KPPRi(imt,Nr)
382    #endif /* ALLOW_DIAGNOSTICS */
383    
384  c find bulk Richardson number at every grid level until > Ricr  c find bulk Richardson number at every grid level until > Ricr
385  c  c
# Line 312  c       kbl(i)=kmtj(i) and hbl(i)=-zgrid Line 393  c       kbl(i)=kmtj(i) and hbl(i)=-zgrid
393  c     initialize hbl and kbl to bottomed out values  c     initialize hbl and kbl to bottomed out values
394    
395        do i = 1, imt        do i = 1, imt
396           Rib(i,1) = 0.0           Rib(i,1) = 0. _d 0
397           kbl(i) = max(kmtj(i),1)           kbl(i) = max(kmtj(i),1)
398           hbl(i) = -zgrid(kbl(i))           hbl(i) = -zgrid(kbl(i))
399        end do        end do
400    
401    #ifdef ALLOW_DIAGNOSTICS
402          do kl = 1, Nr
403             do i = 1, imt
404                KPPBFSFC(i,kl) = 0. _d 0
405                KPPRi(i,kl) = 0. _d 0
406             enddo
407          enddo
408    #endif /* ALLOW_DIAGNOSTICS */
409    
410        do kl = 2, Nr        do kl = 2, Nr
411    
412    #ifdef ALLOW_AUTODIFF_TAMC
413             kkppkey = (ikppkey-1)*Nr + kl
414    #endif
415    
416  c     compute bfsfc = sw fraction at hbf * zgrid  c     compute bfsfc = sw fraction at hbf * zgrid
417    
418           do i = 1, imt           do i = 1, imt
419              worka(i) = zgrid(kl)              worka(i) = zgrid(kl)
420           end do           end do
421    CADJ store worka = comlev1_kpp_k, key = kkppkey, kind=isbyte
422           call SWFRAC(           call SWFRAC(
423       I        imt, hbf,       I       imt, hbf,
424       I        mytime, mythid,       U       worka,
425       U        worka )       I       myTime, myIter, myThid )
426    CADJ store worka = comlev1_kpp_k, key = kkppkey, kind=isbyte
427    
428           do i = 1, imt           do i = 1, imt
429    
# Line 338  c     use caseA as temporary array Line 434  c     use caseA as temporary array
434  c     compute bfsfc= Bo + radiative contribution down to hbf * hbl  c     compute bfsfc= Bo + radiative contribution down to hbf * hbl
435    
436              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  
437    
438           end do           end do
439    #ifdef ALLOW_SALT_PLUME
440    c     compute bfsfc = plume fraction at hbf * zgrid
441             IF ( useSALT_PLUME ) THEN
442               do i = 1, imt
443                  worka(i) = zgrid(kl)
444               enddo
445               call SALT_PLUME_FRAC(
446         I         imt, hbf,SPDepth,
447         U         worka,
448         I         myTime, myIter, myThid)
449               do i = 1, imt
450                  bfsfc(i) = bfsfc(i) + boplume(i)*(worka(i))
451               enddo
452             ENDIF
453    #endif /* ALLOW_SALT_PLUME */
454    
455    #ifdef ALLOW_DIAGNOSTICS
456             do i = 1, imt
457                KPPBFSFC(i,kl) = bfsfc(i)
458             enddo
459    #endif /* ALLOW_DIAGNOSTICS */
460    
461             do i = 1, imt
462                stable(i) = p5 + sign(p5,bfsfc(i))
463                sigma(i) = stable(i) + (1. - stable(i)) * epsilon
464             enddo
465    
466  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
467  c     compute velocity scales at sigma, for hbl= caseA = -zgrid(kl)  c     compute velocity scales at sigma, for hbl= caseA = -zgrid(kl)
468  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
469    
470           call wscale (           call wscale (
471       I        sigma, casea, ustar, bfsfc,       I        sigma, casea, ustar, bfsfc,
472       O        wm, ws )       O        wm, ws, myThid )
473    CADJ store ws = comlev1_kpp_k, key = kkppkey, kind=isbyte
474    
475           do i = 1, imt           do i = 1, imt
476    
# Line 361  c--------------------------------------- Line 482  c---------------------------------------
482       1           ( dbloc(i,kl-1) / (zgrid(kl-1)-zgrid(kl  ))+       1           ( dbloc(i,kl-1) / (zgrid(kl-1)-zgrid(kl  ))+
483       2             dbloc(i,kl  ) / (zgrid(kl  )-zgrid(kl+1)))       2             dbloc(i,kl  ) / (zgrid(kl  )-zgrid(kl+1)))
484    
485              if (bvsq .eq. 0.) then              if (bvsq .eq. 0. _d 0) then
486                vtsq = 0.0                vtsq = 0. _d 0
487              else              else
488                vtsq = -zgrid(kl) * ws(i) * sqrt(abs(bvsq)) * Vtc                vtsq = -zgrid(kl) * ws(i) * sqrt(abs(bvsq)) * Vtc
489              endif              endif
# Line 379  c     rg: assignment to double precision Line 500  c     rg: assignment to double precision
500  c     ph: test for zero nominator  c     ph: test for zero nominator
501  c  c
502    
503              tempVar1  = dvsq(i,kl) + vtsq                        tempVar1  = dvsq(i,kl) + vtsq
504              tempVar2 = max(tempVar1, phepsi)              tempVar2 = max(tempVar1, phepsi)
505              Rib(i,kl) = Ritop(i,kl) / tempVar2              Rib(i,kl) = Ritop(i,kl) / tempVar2
506    #ifdef ALLOW_DIAGNOSTICS
507                KPPRi(i,kl) = Rib(i,kl)
508    #endif
509    
510           end do           end do
511        end do        end do
512          
513    #ifdef ALLOW_DIAGNOSTICS
514          IF ( useDiagnostics ) THEN
515             CALL DIAGNOSTICS_FILL(KPPBFSFC,'KPPbfsfc',0,Nr,2,bi,bj,myThid)
516             CALL DIAGNOSTICS_FILL(KPPRi   ,'KPPRi   ',0,Nr,2,bi,bj,myThid)
517          ENDIF
518    #endif /* ALLOW_DIAGNOSTICS */
519    
520    cph(
521    cph  without this store, there is a recomputation error for
522    cph  rib in adbldepth (probably partial recomputation problem)
523    CADJ store Rib = comlev1_kpp
524    CADJ &   , key=ikppkey, kind=isbyte,
525    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy),Nr /)
526    cph)
527    
528        do kl = 2, Nr        do kl = 2, Nr
529           do i = 1, imt           do i = 1, imt
530              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 532  c
532        end do        end do
533    
534  CADJ store kbl = comlev1_kpp  CADJ store kbl = comlev1_kpp
535  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key=ikppkey, kind=isbyte,
536    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
537    
538        do i = 1, imt        do i = 1, imt
539           kl = kbl(i)           kl = kbl(i)
# Line 406  c     linearly interpolate to find hbl w Line 546  c     linearly interpolate to find hbl w
546        end do        end do
547    
548  CADJ store hbl = comlev1_kpp  CADJ store hbl = comlev1_kpp
549  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key=ikppkey, kind=isbyte,
550    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
551    
552  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
553  c     find stability and buoyancy forcing for boundary layer  c     find stability and buoyancy forcing for boundary layer
# Line 415  c--------------------------------------- Line 556  c---------------------------------------
556        do i = 1, imt        do i = 1, imt
557           worka(i) = hbl(i)           worka(i) = hbl(i)
558        end do        end do
559    CADJ store worka = comlev1_kpp
560    CADJ &   , key=ikppkey, kind=isbyte,
561    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
562        call SWFRAC(        call SWFRAC(
563       I     imt, minusone,       I       imt, minusone,
564       I     mytime, mythid,       U       worka,
565       U     worka )       I       myTime, myIter, myThid )
566    CADJ store worka = comlev1_kpp
567    CADJ &   , key=ikppkey, kind=isbyte,
568    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
569    
570        do i = 1, imt        do i = 1, imt
571           bfsfc(i)  = bo(i) + bosol(i) * (1. - worka(i))           bfsfc(i)  = bo(i) + bosol(i) * (1. - worka(i))
572        end do        end do
573    
574    #ifdef ALLOW_SALT_PLUME
575          IF ( useSALT_PLUME ) THEN
576            do i = 1, imt
577               worka(i) = hbl(i)
578            enddo
579            call SALT_PLUME_FRAC(
580         I         imt,minusone,SPDepth,
581         U         worka,
582         I         myTime, myIter, myThid )
583            do i = 1, imt
584               bfsfc(i) = bfsfc(i) + boplume(i) * (worka(i))
585            enddo
586          ENDIF
587    #endif /* ALLOW_SALT_PLUME */
588  CADJ store bfsfc = comlev1_kpp  CADJ store bfsfc = comlev1_kpp
589  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key=ikppkey, kind=isbyte,
590    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
591    
592  c--   ensure bfsfc is never 0  c--   ensure bfsfc is never 0
593        do i = 1, imt        do i = 1, imt
# Line 432  c--   ensure bfsfc is never 0 Line 595  c--   ensure bfsfc is never 0
595           bfsfc(i) = sign(eins,bfsfc(i))*max(phepsi,abs(bfsfc(i)))           bfsfc(i) = sign(eins,bfsfc(i))*max(phepsi,abs(bfsfc(i)))
596        end do        end do
597    
598  CADJ store bfsfc = comlev1_kpp  cph(
599  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  cph  added stable to store list to avoid extensive recomp.
600    CADJ store bfsfc, stable = comlev1_kpp
601    CADJ &   , key=ikppkey, kind=isbyte,
602    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
603    cph)
604    
605  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
606  c check hbl limits for hekman or hmonob  c check hbl limits for hekman or hmonob
607  c     ph: test for zero nominator  c     ph: test for zero nominator
608  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
609    
610          IF ( LimitHblStable ) THEN
611        do i = 1, imt        do i = 1, imt
612           if (bfsfc(i) .gt. 0.0) then           if (bfsfc(i) .gt. 0.0) then
613              hekman = cekman * ustar(i) / max(abs(Coriol(i)),phepsi)              hekman = cekman * ustar(i) / max(abs(Coriol(i)),phepsi)
# Line 450  c--------------------------------------- Line 618  c---------------------------------------
618              hbl(i) = min(hbl(i),hlimit)              hbl(i) = min(hbl(i),hlimit)
619           end if           end if
620        end do        end do
621          ENDIF
622    
623  CADJ store hbl = comlev1_kpp  CADJ store hbl = comlev1_kpp
624  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key=ikppkey, kind=isbyte,
625    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
626    
627        do i = 1, imt        do i = 1, imt
628           hbl(i) = max(hbl(i),minKPPhbl)           hbl(i) = max(hbl(i),minKPPhbl)
# Line 459  CADJ &   , key = ikey, shape = (/ (sNx+2 Line 630  CADJ &   , key = ikey, shape = (/ (sNx+2
630        end do        end do
631    
632  CADJ store hbl = comlev1_kpp  CADJ store hbl = comlev1_kpp
633  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key=ikppkey, kind=isbyte,
634    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
635    
636  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
637  c      find new kbl  c      find new kbl
# Line 480  c--------------------------------------- Line 652  c---------------------------------------
652        do i = 1, imt        do i = 1, imt
653           worka(i) = hbl(i)           worka(i) = hbl(i)
654        end do        end do
655    CADJ store worka = comlev1_kpp
656    CADJ &   , key=ikppkey, kind=isbyte,
657    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
658        call SWFRAC(        call SWFRAC(
659       I     imt, minusone,       I     imt, minusone,
660       I     mytime, mythid,       U     worka,
661       U     worka )       I     myTime, myIter, myThid )
662    CADJ store worka = comlev1_kpp
663    CADJ &   , key=ikppkey, kind=isbyte,
664    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
665    
666        do i = 1, imt        do i = 1, imt
667           bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i))           bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i))
668        end do        end do
669    
670    #ifdef ALLOW_SALT_PLUME
671          IF ( useSALT_PLUME ) THEN
672            do i = 1, imt
673               worka(i) = hbl(i)
674            enddo
675            call SALT_PLUME_FRAC(
676         I         imt,minusone,SPDepth,
677         U         worka,
678         I         myTime, myIter, myThid )
679            do i = 1, imt
680               bfsfc(i) = bfsfc(i) + boplume(i) * (worka(i))
681            enddo
682          ENDIF
683    #endif /* ALLOW_SALT_PLUME */
684  CADJ store bfsfc = comlev1_kpp  CADJ store bfsfc = comlev1_kpp
685  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key=ikppkey, kind=isbyte,
686    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
687    
688  c--   ensures bfsfc is never 0  c--   ensures bfsfc is never 0
689        do i = 1, imt        do i = 1, imt
# Line 515  c*************************************** Line 709  c***************************************
709    
710        subroutine wscale (        subroutine wscale (
711       I     sigma, hbl, ustar, bfsfc,       I     sigma, hbl, ustar, bfsfc,
712       O     wm, ws )       O     wm, ws,
713         I     myThid )
714    
715  c     compute turbulent velocity scales.  c     compute turbulent velocity scales.
716  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 719  c
719  c     note: the lookup table is only used for unstable conditions  c     note: the lookup table is only used for unstable conditions
720  c     (zehat.le.0), in the stable domain wm (=ws) gets computed  c     (zehat.le.0), in the stable domain wm (=ws) gets computed
721  c     directly.  c     directly.
722  c  c
723        IMPLICIT NONE        IMPLICIT NONE
724    
725  #include "SIZE.h"  #include "SIZE.h"
# Line 536  c sigma   : normalized depth (d/hbl) Line 731  c sigma   : normalized depth (d/hbl)
731  c hbl     : boundary layer depth (m)  c hbl     : boundary layer depth (m)
732  c ustar   : surface friction velocity         (m/s)  c ustar   : surface friction velocity         (m/s)
733  c bfsfc   : total surface buoyancy flux       (m^2/s^3)  c bfsfc   : total surface buoyancy flux       (m^2/s^3)
734        _KPP_RL sigma(imt)  c myThid  : thread number for this instance of the routine
735        _KPP_RL hbl  (imt)        integer myThid
736        _KPP_RL ustar(imt)        _RL sigma(imt)
737        _KPP_RL bfsfc(imt)        _RL hbl  (imt)
738          _RL ustar(imt)
739          _RL bfsfc(imt)
740    
741  c  output  c  output
742  c--------  c--------
743  c wm, ws  : turbulent velocity scales at sigma  c wm, ws  : turbulent velocity scales at sigma
744        _KPP_RL wm(imt), ws(imt)        _RL wm(imt), ws(imt)
745    
746  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
747    
748  c local  c local
749  c------  c------
750  c zehat   : = zeta *  ustar**3  c zehat   : = zeta *  ustar**3
751        _KPP_RL zehat        _RL zehat
752    
753        integer iz, izp1, ju, i, jup1        integer iz, izp1, ju, i, jup1
754        _KPP_RL udiff, zdiff, zfrac, ufrac, fzfrac, wam        _RL udiff, zdiff, zfrac, ufrac, fzfrac, wam
755        _KPP_RL wbm, was, wbs, u3, tempVar        _RL wbm, was, wbs, u3, tempVar
756    
757  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
758  c use lookup table for zehat < zmax only; otherwise use  c use lookup table for zehat < zmax only; otherwise use
# Line 599  c--------------------------------------- Line 796  c---------------------------------------
796              ws(i) = wm(i)              ws(i) = wm(i)
797    
798           endif           endif
799    
800        end do        end do
801    
802  #endif /* ALLOW_KPP */  #endif /* ALLOW_KPP */
803    
804        return        return
805        end        end
806    
807  c*************************************************************************  c*************************************************************************
808    
809        subroutine Ri_iwmix (        subroutine Ri_iwmix (
810       I       kmtj, shsq, dbloc, dblocSm       I       kmtj, shsq, dbloc, dblocSm,
811       I     , ikey       I       diffusKzS, diffusKzT,
812       O     , diffus )       I       ikppkey,
813         O       diffus,
814         I       myThid )
815    
816  c     compute interior viscosity diffusivity coefficients due  c     compute interior viscosity diffusivity coefficients due
817  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 824  c     to static instability (local Richa
824  #include "EEPARAMS.h"  #include "EEPARAMS.h"
825  #include "PARAMS.h"  #include "PARAMS.h"
826  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
827    #ifdef ALLOW_AUTODIFF
828    # include "AUTODIFF_PARAMS.h"
829    # include "tamc.h"
830    #endif
831    
832  c  input  c  input
833  c     kmtj   (imt)         number of vertical layers on this row  c     kmtj   (imt)         number of vertical layers on this row
834  c     shsq   (imt,Nr)      (local velocity shear)^2               ((m/s)^2)  c     shsq   (imt,Nr)      (local velocity shear)^2               ((m/s)^2)
835  c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)  c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)
836  c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)  c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)
837        integer kmtj   (imt)  c     diffusKzS(imt,Nr)- background vertical diffusivity for scalars    (m^2/s)
838        _KPP_RL shsq   (imt,Nr)  c     diffusKzT(imt,Nr)- background vertical diffusivity for theta      (m^2/s)
839        _KPP_RL dbloc  (imt,Nr)  c     myThid :: My Thread Id. number
840        _KPP_RL dblocSm(imt,Nr)        integer kmtj (imt)
841        integer ikey        _RL shsq     (imt,Nr)
842          _RL dbloc    (imt,Nr)
843          _RL dblocSm  (imt,Nr)
844          _RL diffusKzS(imt,Nr)
845          _RL diffusKzT(imt,Nr)
846          integer ikppkey
847          integer myThid
848    
849  c output  c output
850  c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)  c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)
851  c     diffus(imt,0:Nrp1,2)  vertical scalar diffusivity             (m^2/s)  c     diffus(imt,0:Nrp1,2)  vertical scalar diffusivity             (m^2/s)
852  c     diffus(imt,0:Nrp1,3)  vertical temperature diffusivity        (m^2/s)  c     diffus(imt,0:Nrp1,3)  vertical temperature diffusivity        (m^2/s)
853        _KPP_RL diffus(imt,0:Nrp1,3)        _RL diffus(imt,0:Nrp1,3)
854    
855  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
856    
857  c local variables  c local variables
858  c     Rig                   local Richardson number  c     Rig                   local Richardson number
859  c     fRi, fcon             function of Rig  c     fRi, fcon             function of Rig
860        _KPP_RL Rig        _RL Rig
861        _KPP_RL fRi, fcon        _RL fRi, fcon
862        _KPP_RL ratio        _RL ratio
863        integer i, ki, mr        integer i, ki, kp1
864        _KPP_RL c1, c0        _RL c1, c0
865    
866  #ifdef ALLOW_KPP_VERTICALLY_SMOOTH  #ifdef ALLOW_KPP_VERTICALLY_SMOOTH
867          integer mr
868  CADJ INIT kpp_ri_tape_mr = common, 1  CADJ INIT kpp_ri_tape_mr = common, 1
869  #endif  #endif
870    
871  c constants  c constants
872        c1 = 1.0        c1 = 1. _d 0
873        c0 = 0.0        c0 = 0. _d 0
874    
875  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
876  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 880  c     set values at bottom and below to
880  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
881  C     break data flow dependence on diffus  C     break data flow dependence on diffus
882        diffus(1,1,1) = 0.0        diffus(1,1,1) = 0.0
883    
884          do ki = 1, Nr
885             do i = 1, imt
886                diffus(i,ki,1) = 0.
887                diffus(i,ki,2) = 0.
888                diffus(i,ki,3) = 0.
889             enddo
890          enddo
891  #endif  #endif
892    
893        do ki = 1, Nr        do ki = 1, Nr
894           do i = 1, imt           do i = 1, imt
895              if     (kmtj(i) .EQ. 0      ) then              if     (kmtj(i) .LE. 1      ) then
896                 diffus(i,ki,1) = 0.                 diffus(i,ki,1) = 0.
897                 diffus(i,ki,2) = 0.                 diffus(i,ki,2) = 0.
898              elseif (ki      .GE. kmtj(i)) then              elseif (ki      .GE. kmtj(i)) then
# Line 687  C     break data flow dependence on diff Line 905  C     break data flow dependence on diff
905              endif              endif
906           end do           end do
907        end do        end do
908    CADJ store diffus = comlev1_kpp, key=ikppkey, kind=isbyte
909    
910  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
911  c     vertically smooth Ri  c     vertically smooth Ri
# Line 697  CADJ store diffus(:,:,1) = kpp_ri_tape_m Line 916  CADJ store diffus(:,:,1) = kpp_ri_tape_m
916  CADJ &  , key=mr, shape=(/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2 /)  CADJ &  , key=mr, shape=(/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2 /)
917    
918           call z121 (           call z121 (
919       U     diffus(1,0,1))       U     diffus(1,0,1),
920         I     myThid )
921        end do        end do
922  #endif  #endif
923    
# Line 706  c                           after smooth Line 926  c                           after smooth
926    
927        do ki = 1, Nr        do ki = 1, Nr
928           do i = 1, imt           do i = 1, imt
929    
930  c  evaluate f of Brunt-Vaisala squared for convection, store in fcon  c  evaluate f of Brunt-Vaisala squared for convection, store in fcon
931    
932              Rig   = max ( diffus(i,ki,2) , BVSQcon )              Rig   = max ( diffus(i,ki,2) , BVSQcon )
933              ratio = min ( (BVSQcon - Rig) / BVSQcon, c1 )              ratio = min ( (BVSQcon - Rig) / BVSQcon, c1 )
934              fcon  = c1 - ratio * ratio              fcon  = c1 - ratio * ratio
935              fcon  = fcon * fcon * fcon              fcon  = fcon * fcon * fcon
936    
937  c  evaluate f of smooth Ri for shear instability, store in fRi  c  evaluate f of smooth Ri for shear instability, store in fRi
938    
939              Rig  = max ( diffus(i,ki,1), c0 )              Rig  = max ( diffus(i,ki,1), c0 )
940              ratio = min ( Rig / Riinfty , c1 )              ratio = min ( Rig / Riinfty , c1 )
941              fRi   = c1 - ratio * ratio              fRi   = c1 - ratio * ratio
942              fRi   = fRi * fRi * fRi              fRi   = fRi * fRi * fRi
943    
944  c ----------------------------------------------------------------------  c ----------------------------------------------------------------------
945  c            evaluate diffusivities and viscosity  c            evaluate diffusivities and viscosity
946  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  
947    
948                kp1 = MIN(ki+1,Nr)
949    #ifdef EXCLUDE_KPP_SHEAR_MIX
950                diffus(i,ki,1) = viscArNr(1)
951                diffus(i,ki,2) = diffusKzS(i,kp1)
952                diffus(i,ki,3) = diffusKzT(i,kp1)
953    #else /* EXCLUDE_KPP_SHEAR_MIX */
954    # ifdef ALLOW_AUTODIFF
955                if ( inAdMode ) then
956                  diffus(i,ki,1) = viscArNr(1)
957                  diffus(i,ki,2) = diffusKzS(i,kp1)
958                  diffus(i,ki,3) = diffusKzT(i,kp1)
959                else
960    # else /* ALLOW_AUTODIFF */
961                if ( .TRUE. ) then
962    # endif /* ALLOW_AUTODIFF */
963                  diffus(i,ki,1) = viscArNr(1) + fcon*difmcon + fRi*difm0
964                  diffus(i,ki,2) = diffusKzS(i,kp1)+fcon*difscon+fRi*difs0
965                  diffus(i,ki,3) = diffusKzT(i,kp1)+fcon*diftcon+fRi*dift0
966                endif
967    #endif /* EXCLUDE_KPP_SHEAR_MIX */
968           end do           end do
969        end do        end do
970    
971  c ------------------------------------------------------------------------  c ------------------------------------------------------------------------
972  c         set surface values to 0.0  c         set surface values to 0.0
973    
974        do i = 1, imt        do i = 1, imt
975           diffus(i,0,1) = c0           diffus(i,0,1) = c0
976           diffus(i,0,2) = c0           diffus(i,0,2) = c0
# Line 742  c         set surface values to 0.0 Line 978  c         set surface values to 0.0
978        end do        end do
979    
980  #endif /* ALLOW_KPP */  #endif /* ALLOW_KPP */
981    
982        return        return
983        end        end
984    
985  c*************************************************************************  c*************************************************************************
986    
987        subroutine z121 (        subroutine z121 (
988       U     v )       U     v,
989         I     myThid )
990    
991  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)
992  c     top (0) value is used as a dummy  c     top (0) value is used as a dummy
# Line 764  c     penetrative convection. Line 1001  c     penetrative convection.
1001  #include "SIZE.h"  #include "SIZE.h"
1002  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
1003    
1004  c input/output  c input/output
1005  c-------------  c-------------
1006  c v     : 2-D array to be smoothed in Nrp1 direction  c v     : 2-D array to be smoothed in Nrp1 direction
1007        _KPP_RL v(imt,0:Nrp1)  c myThid: thread number for this instance of the routine
1008          integer myThid
1009          _RL v(imt,0:Nrp1)
1010    
1011  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1012    
1013  c local  c local
1014        _KPP_RL zwork, zflag        _RL zwork, zflag
1015        _KPP_RL KRi_range(1:Nrp1)        _RL KRi_range(1:Nrp1)
1016        integer i, k, km1, kp1        integer i, k, km1, kp1
1017    
1018        _KPP_RL     p0      , p25       , p5      , p2        _RL         p0      , p25       , p5      , p2
1019        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 )
1020    
1021        KRi_range(Nrp1) = p0        KRi_range(Nrp1) = p0
# Line 830  CADJ STORE v(i,k), zwork = z121tape Line 1069  CADJ STORE v(i,k), zwork = z121tape
1069    
1070  c*************************************************************************  c*************************************************************************
1071    
       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*************************************************************************  
   
1072        subroutine smooth_horiz (        subroutine smooth_horiz (
1073       I     k, bi, bj,       I     k, bi, bj,
1074       U     fld )       U     fld,
1075         I     myThid )
1076    
1077  c     Apply horizontal smoothing to global _RL 2-D array  c     Apply horizontal smoothing to global _RL 2-D array
1078    
1079        IMPLICIT NONE        IMPLICIT NONE
1080  #include "SIZE.h"  #include "SIZE.h"
1081    #include "GRID.h"
1082  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
1083    
1084  c     input  c     input
1085  c     bi, bj : array indices  c     bi, bj : array indices
1086  c     k      : vertical index used for masking  c     k      : vertical index used for masking
1087    c     myThid : thread number for this instance of the routine
1088          INTEGER myThid
1089        integer k, bi, bj        integer k, bi, bj
1090    
1091  c     input/output  c     input/output
# Line 949  c     local Line 1112  c     local
1112              im1 = i-1              im1 = i-1
1113              ip1 = i+1              ip1 = i+1
1114              tempVar =              tempVar =
1115       &           p25   *   pMask(i  ,j  ,k,bi,bj)   +       &           p25   *   maskC(i  ,j  ,k,bi,bj)   +
1116       &           p125  * ( pMask(im1,j  ,k,bi,bj)   +       &           p125  * ( maskC(im1,j  ,k,bi,bj)   +
1117       &                     pMask(ip1,j  ,k,bi,bj)   +       &                     maskC(ip1,j  ,k,bi,bj)   +
1118       &                     pMask(i  ,jm1,k,bi,bj)   +       &                     maskC(i  ,jm1,k,bi,bj)   +
1119       &                     pMask(i  ,jp1,k,bi,bj) ) +       &                     maskC(i  ,jp1,k,bi,bj) ) +
1120       &           p0625 * ( pMask(im1,jm1,k,bi,bj)   +       &           p0625 * ( maskC(im1,jm1,k,bi,bj)   +
1121       &                     pMask(im1,jp1,k,bi,bj)   +       &                     maskC(im1,jp1,k,bi,bj)   +
1122       &                     pMask(ip1,jm1,k,bi,bj)   +       &                     maskC(ip1,jm1,k,bi,bj)   +
1123       &                     pMask(ip1,jp1,k,bi,bj) )       &                     maskC(ip1,jp1,k,bi,bj) )
1124              IF ( tempVar .GE. p25 ) THEN              IF ( tempVar .GE. p25 ) THEN
1125                 fld_tmp(i,j) = (                 fld_tmp(i,j) = (
1126       &              p25  * fld(i  ,j  )*pMask(i  ,j  ,k,bi,bj) +       &              p25  * fld(i  ,j  )*maskC(i  ,j  ,k,bi,bj) +
1127       &              p125 *(fld(im1,j  )*pMask(im1,j  ,k,bi,bj) +       &              p125 *(fld(im1,j  )*maskC(im1,j  ,k,bi,bj) +
1128       &                     fld(ip1,j  )*pMask(ip1,j  ,k,bi,bj) +       &                     fld(ip1,j  )*maskC(ip1,j  ,k,bi,bj) +
1129       &                     fld(i  ,jm1)*pMask(i  ,jm1,k,bi,bj) +       &                     fld(i  ,jm1)*maskC(i  ,jm1,k,bi,bj) +
1130       &                     fld(i  ,jp1)*pMask(i  ,jp1,k,bi,bj))+       &                     fld(i  ,jp1)*maskC(i  ,jp1,k,bi,bj))+
1131       &              p0625*(fld(im1,jm1)*pMask(im1,jm1,k,bi,bj) +       &              p0625*(fld(im1,jm1)*maskC(im1,jm1,k,bi,bj) +
1132       &                     fld(im1,jp1)*pMask(im1,jp1,k,bi,bj) +       &                     fld(im1,jp1)*maskC(im1,jp1,k,bi,bj) +
1133       &                     fld(ip1,jm1)*pMask(ip1,jm1,k,bi,bj) +       &                     fld(ip1,jm1)*maskC(ip1,jm1,k,bi,bj) +
1134       &                     fld(ip1,jp1)*pMask(ip1,jp1,k,bi,bj)))       &                     fld(ip1,jp1)*maskC(ip1,jp1,k,bi,bj)))
1135       &              / tempVar       &              / tempVar
1136              ELSE              ELSE
1137                 fld_tmp(i,j) = fld(i,j)                 fld_tmp(i,j) = fld(i,j)
# Line 992  c*************************************** Line 1155  c***************************************
1155    
1156        subroutine blmix (        subroutine blmix (
1157       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl
1158       O     , dkm1, blmc, ghat, sigma, ikey       O     , dkm1, blmc, ghat, sigma, ikppkey
1159       &     )       I     , myThid )
1160    
1161  c     mixing coefficients within boundary layer depend on surface  c     mixing coefficients within boundary layer depend on surface
1162  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 1170  c
1170    
1171  #include "SIZE.h"  #include "SIZE.h"
1172  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
1173    #ifdef ALLOW_AUTODIFF
1174    # include "tamc.h"
1175    #endif
1176    
1177  c input  c input
1178  c     ustar (imt)                 surface friction velocity             (m/s)  c     ustar (imt)                 surface friction velocity             (m/s)
# Line 1015  c     hbl   (imt)                 bounda Line 1181  c     hbl   (imt)                 bounda
1181  c     stable(imt)                 = 1 in stable forcing  c     stable(imt)                 = 1 in stable forcing
1182  c     casea (imt)                 = 1 in case A  c     casea (imt)                 = 1 in case A
1183  c     diffus(imt,0:Nrp1,mdiff)    vertical diffusivities              (m^2/s)  c     diffus(imt,0:Nrp1,mdiff)    vertical diffusivities              (m^2/s)
1184  c     kbl(imt)                    -1 of first grid level below hbl  c     kbl   (imt)                 -1 of first grid level below hbl
1185        _KPP_RL ustar (imt)  c     myThid               thread number for this instance of the routine
1186        _KPP_RL bfsfc (imt)        integer myThid
1187        _KPP_RL hbl   (imt)        _RL ustar (imt)
1188        _KPP_RL stable(imt)        _RL bfsfc (imt)
1189        _KPP_RL casea (imt)        _RL hbl   (imt)
1190        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL stable(imt)
1191          _RL casea (imt)
1192          _RL diffus(imt,0:Nrp1,mdiff)
1193        integer kbl(imt)        integer kbl(imt)
1194    
1195  c output  c output
# Line 1029  c     dkm1 (imt,mdiff)            bounda Line 1197  c     dkm1 (imt,mdiff)            bounda
1197  c     blmc (imt,Nr,mdiff)         boundary layer mixing coefficients  (m^2/s)  c     blmc (imt,Nr,mdiff)         boundary layer mixing coefficients  (m^2/s)
1198  c     ghat (imt,Nr)               nonlocal scalar transport  c     ghat (imt,Nr)               nonlocal scalar transport
1199  c     sigma(imt)                  normalized depth (d / hbl)  c     sigma(imt)                  normalized depth (d / hbl)
1200        _KPP_RL dkm1 (imt,mdiff)        _RL dkm1 (imt,mdiff)
1201        _KPP_RL blmc (imt,Nr,mdiff)        _RL blmc (imt,Nr,mdiff)
1202        _KPP_RL ghat (imt,Nr)        _RL ghat (imt,Nr)
1203        _KPP_RL sigma(imt)        _RL sigma(imt)
1204        integer ikey        integer ikppkey
1205    
1206  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1207    
# Line 1041  c  local Line 1209  c  local
1209  c     gat1*(imt)                 shape function at sigma = 1  c     gat1*(imt)                 shape function at sigma = 1
1210  c     dat1*(imt)                 derivative of shape function at sigma = 1  c     dat1*(imt)                 derivative of shape function at sigma = 1
1211  c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)  c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)
1212        _KPP_RL gat1m(imt), gat1s(imt), gat1t(imt)        _RL gat1m(imt), gat1s(imt), gat1t(imt)
1213        _KPP_RL dat1m(imt), dat1s(imt), dat1t(imt)        _RL dat1m(imt), dat1s(imt), dat1t(imt)
1214        _KPP_RL ws(imt), wm(imt)        _RL ws(imt), wm(imt)
1215        integer i, kn, ki        integer i, kn, ki
1216        _KPP_RL R, dvdzup, dvdzdn, viscp        _RL R, dvdzup, dvdzdn, viscp
1217        _KPP_RL difsp, diftp, visch, difsh, difth        _RL difsp, diftp, visch, difsh, difth
1218        _KPP_RL f1, sig, a1, a2, a3, delhat        _RL f1, sig, a1, a2, a3, delhat
1219        _KPP_RL Gm, Gs, Gt        _RL Gm, Gs, Gt
1220        _KPP_RL tempVar        _RL tempVar
1221    
1222        _KPP_RL    p0    , eins        _RL    p0    , eins
1223        parameter (p0=0.0, eins=1.0)        parameter (p0=0.0, eins=1.0)
1224    #ifdef ALLOW_AUTODIFF_TAMC
1225          integer kkppkey
1226    #endif
1227    
1228  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1229  c compute velocity scales at hbl  c compute velocity scales at hbl
# Line 1062  c--------------------------------------- Line 1233  c---------------------------------------
1233           sigma(i) = stable(i) * 1.0 + (1. - stable(i)) * epsilon           sigma(i) = stable(i) * 1.0 + (1. - stable(i)) * epsilon
1234        end do        end do
1235    
1236    CADJ STORE sigma = comlev1_kpp, key=ikppkey, kind=isbyte
1237        call wscale (        call wscale (
1238       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1239       O        wm, ws )       O        wm, ws, myThid )
1240    CADJ STORE wm = comlev1_kpp, key=ikppkey, kind=isbyte
1241    CADJ STORE ws = comlev1_kpp, key=ikppkey, kind=isbyte
1242    
1243        do i = 1, imt        do i = 1, imt
1244           wm(i) = sign(eins,wm(i))*max(phepsi,abs(wm(i)))           wm(i) = sign(eins,wm(i))*max(phepsi,abs(wm(i)))
1245           ws(i) = sign(eins,ws(i))*max(phepsi,abs(ws(i)))           ws(i) = sign(eins,ws(i))*max(phepsi,abs(ws(i)))
1246        end do        end do
1247  CADJ STORE wm = comlev1_kpp, key = ikey  CADJ STORE wm = comlev1_kpp, key=ikppkey, kind=isbyte
1248  CADJ STORE ws = comlev1_kpp, key = ikey  CADJ STORE ws = comlev1_kpp, key=ikppkey, kind=isbyte
1249    
1250        do i = 1, imt        do i = 1, imt
1251    
# Line 1103  c--------------------------------------- Line 1277  c---------------------------------------
1277           difsh  = diffus(i,kn,2) + difsp * delhat           difsh  = diffus(i,kn,2) + difsp * delhat
1278           difth  = diffus(i,kn,3) + diftp * delhat           difth  = diffus(i,kn,3) + diftp * delhat
1279    
1280           f1 = stable(i) * conc1 * bfsfc(i) /           f1 = stable(i) * conc1 * bfsfc(i) /
1281       &        max(ustar(i)**4,phepsi)       &        max(ustar(i)**4,phepsi)
1282           gat1m(i) = visch / hbl(i) / wm(i)           gat1m(i) = visch / hbl(i) / wm(i)
1283           dat1m(i) = -viscp / wm(i) + f1 * visch           dat1m(i) = -viscp / wm(i) + f1 * visch
1284           dat1m(i) = min(dat1m(i),p0)  
   
1285           gat1s(i) = difsh  / hbl(i) / ws(i)           gat1s(i) = difsh  / hbl(i) / ws(i)
1286           dat1s(i) = -difsp / ws(i) + f1 * difsh           dat1s(i) = -difsp / ws(i) + f1 * difsh
1287           dat1s(i) = min(dat1s(i),p0)  
   
1288           gat1t(i) = difth /  hbl(i) / ws(i)           gat1t(i) = difth /  hbl(i) / ws(i)
1289           dat1t(i) = -diftp / ws(i) + f1 * difth           dat1t(i) = -diftp / ws(i) + f1 * difth
1290    
1291          end do
1292    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1293    CADJ STORE gat1m = comlev1_kpp, key=ikppkey, kind=isbyte
1294    CADJ STORE gat1s = comlev1_kpp, key=ikppkey, kind=isbyte
1295    CADJ STORE gat1t = comlev1_kpp, key=ikppkey, kind=isbyte
1296    CADJ STORE dat1m = comlev1_kpp, key=ikppkey, kind=isbyte
1297    CADJ STORE dat1s = comlev1_kpp, key=ikppkey, kind=isbyte
1298    CADJ STORE dat1t = comlev1_kpp, key=ikppkey, kind=isbyte
1299    #endif
1300          do i = 1, imt
1301             dat1m(i) = min(dat1m(i),p0)
1302             dat1s(i) = min(dat1s(i),p0)
1303           dat1t(i) = min(dat1t(i),p0)           dat1t(i) = min(dat1t(i),p0)
   
1304        end do        end do
1305    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1306    CADJ STORE dat1m = comlev1_kpp, key=ikppkey, kind=isbyte
1307    CADJ STORE dat1s = comlev1_kpp, key=ikppkey, kind=isbyte
1308    CADJ STORE dat1t = comlev1_kpp, key=ikppkey, kind=isbyte
1309    #endif
1310    
1311        do ki = 1, Nr        do ki = 1, Nr
1312    
1313    #ifdef ALLOW_AUTODIFF_TAMC
1314             kkppkey = (ikppkey-1)*Nr + ki
1315    #endif
1316    
1317  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1318  c     compute turbulent velocity scales on the interfaces  c     compute turbulent velocity scales on the interfaces
1319  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 1129  c--------------------------------------- Line 1322  c---------------------------------------
1322              sig      = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i)              sig      = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i)
1323              sigma(i) = stable(i)*sig + (1.-stable(i))*min(sig,epsilon)              sigma(i) = stable(i)*sig + (1.-stable(i))*min(sig,epsilon)
1324           end do           end do
1325    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1326    CADJ STORE wm = comlev1_kpp_k, key = kkppkey
1327    CADJ STORE ws = comlev1_kpp_k, key = kkppkey
1328    #endif
1329    CADJ STORE sigma = comlev1_kpp_k, key = kkppkey
1330           call wscale (           call wscale (
1331       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1332       O        wm, ws )       O        wm, ws, myThid )
1333    CADJ STORE wm = comlev1_kpp_k, key = kkppkey
1334    CADJ STORE ws = comlev1_kpp_k, key = kkppkey
1335    
1336  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1337  c     compute the dimensionless shape functions at the interfaces  c     compute the dimensionless shape functions at the interfaces
# Line 1161  c--------------------------------------- Line 1361  c---------------------------------------
1361    
1362              tempVar = ws(i) * hbl(i)              tempVar = ws(i) * hbl(i)
1363              ghat(i,ki) = (1.-stable(i)) * cg / max(phepsi,tempVar)              ghat(i,ki) = (1.-stable(i)) * cg / max(phepsi,tempVar)
1364    
1365           end do           end do
1366        end do        end do
1367    
# Line 1171  c--------------------------------------- Line 1371  c---------------------------------------
1371    
1372        do i = 1, imt        do i = 1, imt
1373           sig      = -zgrid(kbl(i)-1) / hbl(i)           sig      = -zgrid(kbl(i)-1) / hbl(i)
1374           sigma(i) = stable(i) * sig           sigma(i) = stable(i) * sig
1375       &            + (1. - stable(i)) * min(sig,epsilon)       &            + (1. - stable(i)) * min(sig,epsilon)
1376        end do        end do
1377    
1378    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1379    CADJ STORE wm = comlev1_kpp, key=ikppkey, kind=isbyte
1380    CADJ STORE ws = comlev1_kpp, key=ikppkey, kind=isbyte
1381    #endif
1382    CADJ STORE sigma = comlev1_kpp, key=ikppkey, kind=isbyte
1383        call wscale (        call wscale (
1384       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1385       O        wm, ws )       O        wm, ws, myThid )
1386    CADJ STORE wm = comlev1_kpp, key=ikppkey, kind=isbyte
1387    CADJ STORE ws = comlev1_kpp, key=ikppkey, kind=isbyte
1388    
1389        do i = 1, imt        do i = 1, imt
1390           sig = -zgrid(kbl(i)-1) / hbl(i)           sig = -zgrid(kbl(i)-1) / hbl(i)
# Line 1194  c--------------------------------------- Line 1401  c---------------------------------------
1401    
1402  #endif /* ALLOW_KPP */  #endif /* ALLOW_KPP */
1403    
1404        return        return
1405        end        end
1406    
1407  c*************************************************************************  c*************************************************************************
1408    
1409        subroutine enhance (        subroutine enhance (
1410       I       dkm1, hbl, kbl, diffus, casea       I       dkm1, hbl, kbl, diffus, casea
1411       U     , ghat       U     , ghat
1412       O     , blmc       O     , blmc
1413       &     )       &     , myThid )
1414    
1415  c enhance the diffusivity at the kbl-.5 interface  c enhance the diffusivity at the kbl-.5 interface
1416    
# Line 1218  c     hbl(imt)                  boundary Line 1425  c     hbl(imt)                  boundary
1425  c     kbl(imt)                  grid above hbl  c     kbl(imt)                  grid above hbl
1426  c     diffus(imt,0:Nrp1,mdiff) vertical diffusivities           (m^2/s)  c     diffus(imt,0:Nrp1,mdiff) vertical diffusivities           (m^2/s)
1427  c     casea(imt)                = 1 in caseA, = 0 in case B  c     casea(imt)                = 1 in caseA, = 0 in case B
1428        _KPP_RL dkm1  (imt,mdiff)  c     myThid                    thread number for this instance of the routine
1429        _KPP_RL hbl   (imt)        integer   myThid
1430          _RL dkm1  (imt,mdiff)
1431          _RL hbl   (imt)
1432        integer kbl   (imt)        integer kbl   (imt)
1433        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL diffus(imt,0:Nrp1,mdiff)
1434        _KPP_RL casea (imt)        _RL casea (imt)
1435    
1436  c input/output  c input/output
1437  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)
1438        _KPP_RL ghat (imt,Nr)        _RL ghat (imt,Nr)
1439    
1440  c output  c output
1441  c     enhanced bound. layer mixing coeff.  c     enhanced bound. layer mixing coeff.
1442        _KPP_RL blmc  (imt,Nr,mdiff)        _RL blmc  (imt,Nr,mdiff)
1443    
1444  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1445    
1446  c local  c local
1447  c     fraction hbl lies beteen zgrid neighbors  c     fraction hbl lies beteen zgrid neighbors
1448        _KPP_RL delta        _RL delta
1449        integer ki, i, md        integer ki, i, md
1450        _KPP_RL dkmp5, dstar        _RL dkmp5, dstar
1451    
1452        do i = 1, imt        do i = 1, imt
1453           ki = kbl(i)-1           ki = kbl(i)-1
# Line 1258  c     fraction hbl lies beteen zgrid nei Line 1467  c     fraction hbl lies beteen zgrid nei
1467    
1468  #endif /* ALLOW_KPP */  #endif /* ALLOW_KPP */
1469    
1470        return        return
1471        end        end
1472    
1473  c*************************************************************************  c*************************************************************************
1474    
1475        SUBROUTINE STATEKPP (        SUBROUTINE STATEKPP (
1476       I     bi, bj, myThid,       O     RHO1, DBLOC, DBSFC, TTALPHA, SSBETA,
1477       O     RHO1, DBLOC, DBSFC, TTALPHA, SSBETA)       I     ikppkey, bi, bj, myThid )
1478  c  c
1479  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1480  c     "statekpp" computes all necessary input arrays  c     "statekpp" computes all necessary input arrays
# Line 1290  c Line 1499  c
1499  c     written  by: jan morzel,   feb. 10, 1995 (converted from "sigma" version)  c     written  by: jan morzel,   feb. 10, 1995 (converted from "sigma" version)
1500  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
1501  c  c
1502    
1503  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1504    
1505        IMPLICIT NONE        IMPLICIT NONE
# Line 1299  c--------------------------------------- Line 1509  c---------------------------------------
1509  #include "PARAMS.h"  #include "PARAMS.h"
1510  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
1511  #include "DYNVARS.h"  #include "DYNVARS.h"
1512    #include "GRID.h"
1513    #ifdef ALLOW_AUTODIFF
1514    # include "tamc.h"
1515    #endif
1516    
1517  c-------------- Routine arguments -----------------------------------------  c-------------- Routine arguments -----------------------------------------
1518        INTEGER bi, bj, myThid        INTEGER bi, bj, myThid
1519        _KPP_RL RHO1   ( ibot:itop, jbot:jtop       )        _RL RHO1   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy       )
1520        _KPP_RL DBLOC  ( ibot:itop, jbot:jtop, Nr   )        _RL DBLOC  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
1521        _KPP_RL DBSFC  ( ibot:itop, jbot:jtop, Nr   )        _RL DBSFC  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
1522        _KPP_RL TTALPHA( ibot:itop, jbot:jtop, Nrp1 )        _RL TTALPHA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1523        _KPP_RL SSBETA ( ibot:itop, jbot:jtop, Nrp1 )        _RL SSBETA ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1524    
1525  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1526    
# Line 1317  c Line 1531  c
1531  c     rhok         - density of t(k  ) & s(k  ) at depth k  c     rhok         - density of t(k  ) & s(k  ) at depth k
1532  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
1533  c     rho1k        - density of t(1  ) & s(1  ) at depth k  c     rho1k        - density of t(1  ) & s(1  ) at depth k
1534  c     work1, work2 - work arrays for holding horizontal slabs  c     work1,2,3    - work arrays for holding horizontal slabs
1535    
1536        _RL RHOK  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL RHOK  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1537        _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 1539  c     work1, work2 - work arrays for hol
1539        _RL WORK1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1540        _RL WORK2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1541        _RL WORK3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1542    
1543        INTEGER I, J, K        INTEGER I, J, K
1544          INTEGER ikppkey, kkppkey
1545    
1546  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
1547    
1548        call FIND_RHO(        kkppkey = (ikppkey-1)*Nr + 1
1549       I     bi, bj, ibot, itop, jbot, jtop, 1, 1, eosType,  
1550       I     theta, salt,  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1551    CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k,
1552    CADJ &     key=kkppkey, kind=isbyte
1553    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k,
1554    CADJ &     key=kkppkey, kind=isbyte
1555    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1556          CALL FIND_RHO_2D(
1557         I     1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1,
1558         I     theta(1-OLx,1-OLy,1,bi,bj), salt(1-OLx,1-OLy,1,bi,bj),
1559       O     WORK1,       O     WORK1,
1560       I     myThid )       I     1, bi, bj, myThid )
1561    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1562    CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k,
1563    CADJ &     key=kkppkey, kind=isbyte
1564    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k,
1565    CADJ &     key=kkppkey, kind=isbyte
1566    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1567    
1568        call FIND_ALPHA(        call FIND_ALPHA(
1569       I     bi, bj, ibot, itop, jbot, jtop, 1, 1, eosType,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1570       O     WORK2 )       O     WORK2, myThid )
1571    
1572        call FIND_BETA(        call FIND_BETA(
1573       I     bi, bj, ibot, itop, jbot, jtop, 1, 1, eosType,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1574       O     WORK3 )       O     WORK3, myThid )
1575    
1576        DO J = jbot, jtop        DO J = 1-OLy, sNy+OLy
1577           DO I = ibot, itop           DO I = 1-OLx, sNx+OLx
1578              RHO1(I,J)      = WORK1(I,J) + rhonil              RHO1(I,J)      = WORK1(I,J) + rhoConst
1579              TTALPHA(I,J,1) = WORK2(I,J)              TTALPHA(I,J,1) = WORK2(I,J)
1580              SSBETA(I,J,1)  = WORK3(I,J)              SSBETA(I,J,1)  = WORK3(I,J)
1581              DBSFC(I,J,1)   = 0.              DBSFC(I,J,1)   = 0.
# Line 1357  c calculate alpha, beta, and gradients i Line 1587  c calculate alpha, beta, and gradients i
1587  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)
1588        DO K = 2, Nr        DO K = 2, Nr
1589    
1590           call FIND_RHO(        kkppkey = (ikppkey-1)*Nr + k
1591       I        bi, bj, ibot, itop, jbot, jtop, K, K, eosType,  
1592       I        theta, salt,  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1593    CADJ STORE theta(:,:,k,bi,bj) = comlev1_kpp_k,
1594    CADJ &     key=kkppkey, kind=isbyte
1595    CADJ STORE salt (:,:,k,bi,bj) = comlev1_kpp_k,
1596    CADJ &     key=kkppkey, kind=isbyte
1597    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1598             CALL FIND_RHO_2D(
1599         I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, k,
1600         I        theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
1601       O        RHOK,       O        RHOK,
1602       I        myThid )       I        k, bi, bj, myThid )
1603    
1604           call FIND_RHO(  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1605       I        bi, bj, ibot, itop, jbot, jtop, K-1, K, eosType,  CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_kpp_k,
1606       I        theta, salt,  CADJ &     key=kkppkey, kind=isbyte
1607    CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_kpp_k,
1608    CADJ &     key=kkppkey, kind=isbyte
1609    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1610             CALL FIND_RHO_2D(
1611         I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, k,
1612         I        theta(1-OLx,1-OLy,k-1,bi,bj),salt(1-OLx,1-OLy,k-1,bi,bj),
1613       O        RHOKM1,       O        RHOKM1,
1614       I        myThid )       I        k-1, bi, bj, myThid )
1615    
1616           call FIND_RHO(  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1617       I        bi, bj, ibot, itop, jbot, jtop, 1, K, eosType,  CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k,
1618       I        theta, salt,  CADJ &     key=kkppkey, kind=isbyte
1619    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k,
1620    CADJ &     key=kkppkey, kind=isbyte
1621    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1622             CALL FIND_RHO_2D(
1623         I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, k,
1624         I        theta(1-OLx,1-OLy,1,bi,bj), salt(1-OLx,1-OLy,1,bi,bj),
1625       O        RHO1K,       O        RHO1K,
1626       I        myThid )       I        1, bi, bj, myThid )
1627    
1628    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1629    CADJ STORE rhok  (:,:)          = comlev1_kpp_k,
1630    CADJ &     key=kkppkey, kind=isbyte
1631    CADJ STORE rhokm1(:,:)          = comlev1_kpp_k,
1632    CADJ &     key=kkppkey, kind=isbyte
1633    CADJ STORE rho1k (:,:)          = 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, K, K, eosType,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1639       O        WORK1 )       O        WORK1, myThid )
1640    
1641           call FIND_BETA(           call FIND_BETA(
1642       I        bi, bj, ibot, itop, jbot, jtop, K, K, eosType,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1643       O        WORK2 )       O        WORK2, 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                 TTALPHA(I,J,K) = WORK1 (I,J)                 TTALPHA(I,J,K) = WORK1 (I,J)
1648                 SSBETA(I,J,K)  = WORK2 (I,J)                 SSBETA(I,J,K)  = WORK2 (I,J)
1649                 DBLOC(I,J,K-1) = gravity * (RHOK(I,J) - RHOKM1(I,J)) /                 DBLOC(I,J,K-1) = gravity * (RHOK(I,J) - RHOKM1(I,J)) /
1650       &                                    (RHOK(I,J) + rhonil)       &                                    (RHOK(I,J) + rhoConst)
1651                 DBSFC(I,J,K)   = gravity * (RHOK(I,J) - RHO1K (I,J)) /                 DBSFC(I,J,K)   = gravity * (RHOK(I,J) - RHO1K (I,J)) /
1652       &                                    (RHOK(I,J) + rhonil)       &                                    (RHOK(I,J) + rhoConst)
1653              END DO              END DO
1654           END DO           END DO
1655    
1656        END DO        END DO
1657    
1658  c     compute arrays for K = Nrp1  c     compute arrays for K = Nrp1
1659        DO J = jbot, jtop        DO J = 1-OLy, sNy+OLy
1660           DO I = ibot, itop           DO I = 1-OLx, sNx+OLx
1661              TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)              TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)
1662              SSBETA(I,J,Nrp1)  = SSBETA(I,J,Nr)              SSBETA(I,J,Nrp1)  = SSBETA(I,J,Nr)
1663              DBLOC(I,J,Nr)     = 0.              DBLOC(I,J,Nr)     = 0.
1664           END DO           END DO
1665        END DO        END DO
1666    
1667    #ifdef ALLOW_DIAGNOSTICS
1668          IF ( useDiagnostics ) THEN
1669             CALL DIAGNOSTICS_FILL(DBSFC ,'KPPdbsfc',0,Nr,2,bi,bj,myThid)
1670             CALL DIAGNOSTICS_FILL(DBLOC ,'KPPdbloc',0,Nr,2,bi,bj,myThid)
1671          ENDIF
1672    #endif /* ALLOW_DIAGNOSTICS */
1673    
1674    #endif /* ALLOW_KPP */
1675    
1676          RETURN
1677          END
1678    
1679    c*************************************************************************
1680    
1681          SUBROUTINE KPP_DOUBLEDIFF (
1682         I     TTALPHA, SSBETA,
1683         U     kappaRT,
1684         U     kappaRS,
1685         I     ikppkey, imin, imax, jmin, jmax, bi, bj, myThid )
1686    c
1687    c-----------------------------------------------------------------------
1688    c     "KPP_DOUBLEDIFF" adds the double diffusive contributions
1689    C     as Rrho-dependent parameterizations to kappaRT and kappaRS
1690    c
1691    c     input:
1692    c     bi, bj  = array indices on which to apply calculations
1693    c     imin, imax, jmin, jmax = array boundaries
1694    c     ikppkey = key for TAMC/TAF automatic differentiation
1695    c     myThid  = thread id
1696    c
1697    c      ttalpha= thermal expansion coefficient without 1/rho factor
1698    c               d(rho) / d(potential temperature)                    (kg/m^3/C)
1699    c      ssbeta = salt expansion coefficient without 1/rho factor
1700    c               d(rho) / d(salinity)                               (kg/m^3/PSU)
1701    c     output: updated
1702    c     kappaRT/S :: background diffusivities for temperature and salinity
1703    c
1704    c     written  by: martin losch,   sept. 15, 2009
1705    c
1706    
1707    c-----------------------------------------------------------------------
1708    
1709          IMPLICIT NONE
1710    
1711    #include "SIZE.h"
1712    #include "EEPARAMS.h"
1713    #include "PARAMS.h"
1714    #include "KPP_PARAMS.h"
1715    #include "DYNVARS.h"
1716    #include "GRID.h"
1717    #ifdef ALLOW_AUTODIFF
1718    # include "tamc.h"
1719    #endif
1720    
1721    c-------------- Routine arguments -----------------------------------------
1722          INTEGER ikppkey, imin, imax, jmin, jmax, bi, bj, myThid
1723    
1724          _RL TTALPHA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1725          _RL SSBETA ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1726          _RL KappaRT( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
1727          _RL KappaRS( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
1728    
1729    #ifdef ALLOW_KPP
1730    
1731    C--------------------------------------------------------------------------
1732    C
1733    C     local variables
1734    C     I,J,K :: loop indices
1735    C     kkppkey :: key for TAMC/TAF automatic differentiation
1736    C
1737          INTEGER I, J, K
1738          INTEGER kkppkey
1739    C     alphaDT   :: d\rho/d\theta * d\theta
1740    C     betaDS    :: d\rho/dsalt * dsalt
1741    C     Rrho      :: "density ratio" R_{\rho} = \alpha dT/dz / \beta dS/dz
1742    C     nuddt/s   :: double diffusive diffusivities
1743    C     numol     :: molecular diffusivity
1744    C     rFac      :: abbreviation for 1/(R_{\rho0}-1)
1745    
1746          _RL alphaDT ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
1747          _RL betaDS  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
1748          _RL nuddt   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
1749          _RL nudds   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
1750          _RL Rrho
1751          _RL numol, rFac, nutmp
1752          INTEGER Km1
1753    
1754    C     set some constants here
1755          numol = 1.5 _d -06
1756          rFac  = 1. _d 0 / (Rrho0 - 1. _d 0 )
1757    C
1758          kkppkey = (ikppkey-1)*Nr + 1
1759    
1760    CML#ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1761    CMLCADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k,
1762    CMLCADJ &     key=kkppkey, kind=isbyte
1763    CMLCADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k,
1764    CMLCADJ &     key=kkppkey, kind=isbyte
1765    CML#endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1766    
1767          DO K = 1, Nr
1768           Km1 = MAX(K-1,1)
1769           DO J = 1-OLy, sNy+OLy
1770            DO I = 1-OLx, sNx+OLx
1771             alphaDT(I,J) = ( theta(I,J,Km1,bi,bj)-theta(I,J,K,bi,bj) )
1772         &        * 0.5 _d 0 * ABS( TTALPHA(I,J,Km1) + TTALPHA(I,J,K) )
1773             betaDS(I,J)  = ( salt(I,J,Km1,bi,bj)-salt(I,J,K,bi,bj) )
1774         &        * 0.5 _d 0 * ( SSBETA(I,J,Km1) + SSBETA(I,J,K) )
1775             nuddt(I,J) = 0. _d 0
1776             nudds(I,J) = 0. _d 0
1777            ENDDO
1778           ENDDO
1779           IF ( K .GT. 1 ) THEN
1780            DO J = jMin, jMax
1781             DO I = iMin, iMax
1782              Rrho  = 0. _d 0
1783    C     Now we have many different cases
1784    C     a. alphaDT > 0 and betaDS > 0 => salt fingering
1785    C        (salinity destabilizes)
1786              IF (      alphaDT(I,J) .GT. betaDS(I,J)
1787         &         .AND. betaDS(I,J) .GT. 0. _d 0 ) THEN
1788               Rrho = MIN( alphaDT(I,J)/betaDS(I,J), Rrho0 )
1789    C     Large et al. 1994, eq. 31a
1790    C          nudds(I,J) = dsfmax * ( 1. _d 0 - (Rrho - 1. _d 0) * rFac )**3
1791               nutmp      =          ( 1. _d 0 - (Rrho - 1. _d 0) * rFac )
1792               nudds(I,J) = dsfmax * nutmp * nutmp * nutmp
1793    C     Large et al. 1994, eq. 31c
1794               nuddt(I,J) = 0.7 _d 0 * nudds(I,J)
1795              ELSEIF (   alphaDT(I,J) .LT. 0. _d 0
1796         &          .AND. betaDS(I,J) .LT. 0. _d 0
1797         &          .AND.alphaDT(I,J) .GT. betaDS(I,J) ) THEN
1798    C     b. alphaDT < 0 and betaDS < 0 => semi-convection, diffusive convection
1799    C        (temperature destabilizes)
1800    C     for Rrho >= 1 the water column is statically unstable and we never
1801    C     reach this point
1802               Rrho = alphaDT(I,J)/betaDS(I,J)
1803    C     Large et al. 1994, eq. 32
1804               nuddt(I,J) = numol * 0.909 _d 0
1805         &          * exp ( 4.6 _d 0 * exp (
1806         &          - 5.4 _d 0 * ( 1. _d 0/Rrho - 1. _d 0 ) ) )
1807    CMLC     or
1808    CMLC     Large et al. 1994, eq. 33
1809    CML         nuddt(I,J) = numol * 8.7 _d 0 * Rrho**1.1
1810    C     Large et al. 1994, eqs. 34
1811               nudds(I,J) = nuddt(I,J) * MAX( 0.15 _d 0 * Rrho,
1812         &          1.85 _d 0 * Rrho - 0.85 _d 0 )
1813              ELSE
1814    C     Do nothing, because in this case the water colume is unstable
1815    C     => double diffusive processes are negligible and mixing due
1816    C     to shear instability will dominate
1817              ENDIF
1818             ENDDO
1819            ENDDO
1820    C     ENDIF ( K .GT. 1 )
1821           ENDIF
1822    C
1823           DO J = 1-OLy, sNy+OLy
1824            DO I = 1-OLx, sNx+OLx
1825             kappaRT(I,J,K) = kappaRT(I,J,K) + nuddt(I,J)
1826             kappaRS(I,J,K) = kappaRS(I,J,K) + nudds(I,J)
1827            ENDDO
1828           ENDDO
1829    #ifdef ALLOW_DIAGNOSTICS
1830           IF ( useDiagnostics ) THEN
1831            CALL DIAGNOSTICS_FILL(nuddt,'KPPnuddt',k,1,2,bi,bj,myThid)
1832            CALL DIAGNOSTICS_FILL(nudds,'KPPnudds',k,1,2,bi,bj,myThid)
1833           ENDIF
1834    #endif /* ALLOW_DIAGNOSTICS */
1835    C     end of K-loop
1836          ENDDO
1837  #endif /* ALLOW_KPP */  #endif /* ALLOW_KPP */
1838    
1839        RETURN        RETURN
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.51
  ViewVC Help
Powered by ViewVC 1.1.22