/[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.40 by dfer, Thu Oct 23 18:11:00 2008 UTC
# Line 11  C--  o BLDEPTH     - Determine oceanic p Line 11  C--  o BLDEPTH     - Determine oceanic p
11  C--  o WSCALE      - Compute turbulent velocity scales.  C--  o WSCALE      - Compute turbulent velocity scales.
12  C--  o RI_IWMIX    - Compute interior viscosity diffusivity coefficients.  C--  o RI_IWMIX    - Compute interior viscosity diffusivity coefficients.
13  C--  o Z121        - Apply 121 vertical smoothing.  C--  o Z121        - Apply 121 vertical smoothing.
14  C--  o KPP_SMOOTH_HORIZ - Apply horizontal smoothing to KPP array.  C--  o SMOOTH_HORIZ- Apply horizontal smoothing to global array.
 C--  o SMOOTH_HORIZ     - Apply horizontal smoothing to global array.  
15  C--  o BLMIX       - Boundary layer mixing coefficients.  C--  o BLMIX       - Boundary layer mixing coefficients.
16  C--  o ENHANCE     - Enhance diffusivity at boundary layer interface.  C--  o ENHANCE     - Enhance diffusivity at boundary layer interface.
17  C--  o STATEKPP    - Compute buoyancy-related input arrays.  C--  o STATEKPP    - Compute buoyancy-related input arrays.
# Line 20  C--  o STATEKPP    - Compute buoyancy-re Line 19  C--  o STATEKPP    - Compute buoyancy-re
19  c***********************************************************************  c***********************************************************************
20    
21        SUBROUTINE KPPMIX (        SUBROUTINE KPPMIX (
22       I       mytime, mythid       I       kmtj, shsq, dvsq, ustar, msk
23       I     , kmtj, shsq, dvsq, ustar       I     , bo, bosol
24       I     , bo, bosol, dbloc, Ritop, coriol  #ifdef ALLOW_SALT_PLUME
25       I     , ikey       I     , boplume,SPDepth
26    #endif /* ALLOW_SALT_PLUME */
27         I     , dbloc, Ritop, coriol
28         I     , diffusKzS, diffusKzT
29         I     , ikppkey
30       O     , diffus       O     , diffus
31       U     , ghat       U     , ghat
32       O     , hbl )       O     , hbl
33         I     , myTime, myIter, myThid )
34    
35  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
36  c  c
# Line 47  c--------------------------------------- Line 51  c---------------------------------------
51  #include "SIZE.h"  #include "SIZE.h"
52  #include "EEPARAMS.h"  #include "EEPARAMS.h"
53  #include "PARAMS.h"  #include "PARAMS.h"
 #include "DYNVARS.h"  
 #include "FFIELDS.h"  
54  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
55    
56  c input  c input
57  c     myTime          - current time in simulation  c   myTime :: Current time in simulation
58  c     myThid          - thread number for this instance of the routine  c   myIter :: Current iteration number in simulation
59  c     kmtj   (imt)    - number of vertical layers on this row  c   myThid :: My Thread Id. number
60  c     shsq   (imt,Nr) - (local velocity shear)^2                      ((m/s)^2)  c     kmtj   (imt)     - number of vertical layers on this row
61  c     dvsq   (imt,Nr) - (velocity shear re sfc)^2                     ((m/s)^2)  c     msk    (imt)     - surface mask (=1 if water, =0 otherwise)
62  c     ustar  (imt)    - surface friction velocity                         (m/s)  c     shsq   (imt,Nr)  - (local velocity shear)^2                     ((m/s)^2)
63  c     bo     (imt)    - surface turbulent buoy. forcing               (m^2/s^3)  c     dvsq   (imt,Nr)  - (velocity shear re sfc)^2                    ((m/s)^2)
64  c     bosol  (imt)    - radiative buoyancy forcing                    (m^2/s^3)  c     ustar  (imt)     - surface friction velocity                        (m/s)
65  c     dbloc  (imt,Nr) - local delta buoyancy across interfaces          (m/s^2)  c     bo     (imt)     - surface turbulent buoy. forcing              (m^2/s^3)
66  c     dblocSm(imt,Nr) - horizontally smoothed dbloc                     (m/s^2)  c     bosol  (imt)     - radiative buoyancy forcing                   (m^2/s^3)
67  c                         stored in ghat to save space  c     boplume(imt)     - haline buoyancy forcing                      (m^2/s^3)
68  c     Ritop  (imt,Nr) - numerator of bulk Richardson Number  c     dbloc  (imt,Nr)  - local delta buoyancy across interfaces         (m/s^2)
69  c                         (zref-z) * delta buoyancy w.r.t. surface    ((m/s)^2)  c     dblocSm(imt,Nr)  - horizontally smoothed dbloc                    (m/s^2)
70  c     coriol (imt)    - Coriolis parameter                                (1/s)  c                          stored in ghat to save space
71    c     Ritop  (imt,Nr)  - numerator of bulk Richardson Number
72    c                          (zref-z) * delta buoyancy w.r.t. surface   ((m/s)^2)
73    c     coriol (imt)     - Coriolis parameter                               (1/s)
74    c     diffusKzS(imt,Nr)- background vertical diffusivity for scalars    (m^2/s)
75    c     diffusKzT(imt,Nr)- background vertical diffusivity for theta      (m^2/s)
76  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,
77  c           e.g., hbl(sNx,sNy) -> hbl(imt),  c           e.g., hbl(sNx,sNy) -> hbl(imt),
78  c           where hbl(i,j) -> hbl((j-1)*sNx+i)  c           where hbl(i,j) -> hbl((j-1)*sNx+i)
79    
80        _RL     mytime        _RL     myTime
81        integer mythid        integer myIter
82        integer kmtj  (imt   )        integer myThid
83        _KPP_RL shsq  (imt,Nr)        integer kmtj (imt   )
84        _KPP_RL dvsq  (imt,Nr)        _RL shsq     (imt,Nr)
85        _KPP_RL ustar (imt   )        _RL dvsq     (imt,Nr)
86        _KPP_RL bo    (imt   )        _RL ustar    (imt   )
87        _KPP_RL bosol (imt   )        _RL bo       (imt   )
88        _KPP_RL dbloc (imt,Nr)        _RL bosol    (imt   )
89        _KPP_RL Ritop (imt,Nr)  #ifdef ALLOW_SALT_PLUME
90        _KPP_RL coriol(imt   )        _RL boplume  (imt   )
91          _RL SPDepth  (imt   )
92    #endif /* ALLOW_SALT_PLUME */
93          _RL dbloc    (imt,Nr)
94          _RL Ritop    (imt,Nr)
95          _RL coriol   (imt   )
96          _RS msk      (imt   )
97          _RL diffusKzS(imt,Nr)
98          _RL diffusKzT(imt,Nr)
99    
100        integer ikey        integer ikppkey
101    
102  c output  c output
103  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 106  c     diffus (imt,3)  - vertical tempera
106  c     ghat   (imt)    - nonlocal transport coefficient                  (s/m^2)  c     ghat   (imt)    - nonlocal transport coefficient                  (s/m^2)
107  c     hbl    (imt)    - mixing layer depth                                  (m)  c     hbl    (imt)    - mixing layer depth                                  (m)
108    
109        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL diffus(imt,0:Nrp1,mdiff)
110        _KPP_RL ghat  (imt,Nr)        _RL ghat  (imt,Nr)
111        _KPP_RL hbl   (imt)        _RL hbl   (imt)
112    
113  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
114    
# Line 107  c     blmc   (imt,Nr,mdiff) - boundary l Line 122  c     blmc   (imt,Nr,mdiff) - boundary l
122  c     sigma  (imt         ) - normalized depth (d / hbl)  c     sigma  (imt         ) - normalized depth (d / hbl)
123  c     Rib    (imt,Nr      ) - bulk Richardson number  c     Rib    (imt,Nr      ) - bulk Richardson number
124    
125        integer kbl   (imt         )        integer kbl(imt         )
126        _KPP_RL bfsfc (imt         )        _RL bfsfc  (imt         )
127        _KPP_RL casea (imt         )        _RL casea  (imt         )
128        _KPP_RL stable(imt         )        _RL stable (imt         )
129        _KPP_RL dkm1  (imt,   mdiff)        _RL dkm1   (imt,   mdiff)
130        _KPP_RL blmc  (imt,Nr,mdiff)        _RL blmc   (imt,Nr,mdiff)
131        _KPP_RL sigma (imt         )        _RL sigma  (imt         )
132        _KPP_RL Rib   (imt,Nr      )        _RL Rib    (imt,Nr      )
133    
134        integer i, k, md        integer i, k, md
135    
# Line 125  c instability. Line 140  c instability.
140  c (ghat is temporary storage for horizontally smoothed dbloc)  c (ghat is temporary storage for horizontally smoothed dbloc)
141  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
142    
143  CADJ STORE ghat = comlev1_kpp, key = ikey  cph(
144    cph these storings avoid recomp. of Ri_iwmix
145    CADJ STORE ghat  = comlev1_kpp, key = ikppkey
146    CADJ STORE dbloc = comlev1_kpp, key = ikppkey
147    cph)
148        call Ri_iwmix (        call Ri_iwmix (
149       I       kmtj, shsq, dbloc, ghat       I       kmtj, shsq, dbloc, ghat
150       I     , ikey       I     , diffusKzS, diffusKzT
151       O     , diffus )       I     , ikppkey
152         O     , diffus, myThid )
153    
154    cph(
155    cph these storings avoid recomp. of Ri_iwmix
156    cph DESPITE TAFs 'not necessary' warning!
157    CADJ STORE dbloc  = comlev1_kpp, key = ikppkey
158    CADJ STORE shsq   = comlev1_kpp, key = ikppkey
159    CADJ STORE ghat   = comlev1_kpp, key = ikppkey
160    CADJ STORE diffus = comlev1_kpp, key = ikppkey
161    cph)
162    
163  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
164  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 166  c for blmix
166  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
167    
168        do md = 1, mdiff        do md = 1, mdiff
169           do k=1,Nrp1
170           do i = 1,imt           do i = 1,imt
171              do k=kmtj(i),Nrp1               if(k.ge.kmtj(i))  diffus(i,k,md) = 0.0
                diffus(i,k,md) = 0.0  
172              end do              end do
173           end do           end do
174        end do        end do
# Line 152  c diagnose the new boundary layer depth Line 180  c diagnose the new boundary layer depth
180  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
181    
182        call bldepth (        call bldepth (
183       I       mytime, mythid       I       kmtj
184       I     , kmtj       I     , dvsq, dbloc, Ritop, ustar, bo, bosol
185       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol  #ifdef ALLOW_SALT_PLUME
186       I     , ikey       I     , boplume,SPDepth
187    #endif /* ALLOW_SALT_PLUME */
188         I     , coriol
189         I     , ikppkey
190       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma
191       &     )       I     , myTime, myIter, myThid )
192    
193  CADJ STORE hbl,bfsfc,stable,casea,kbl = comlev1_kpp, key = ikey  CADJ STORE hbl,bfsfc,stable,casea,kbl = comlev1_kpp, key = ikppkey
194    
195  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
196  c compute boundary layer diffusivities  c compute boundary layer diffusivities
# Line 167  c--------------------------------------- Line 198  c---------------------------------------
198    
199        call blmix (        call blmix (
200       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl
201       O     , dkm1, blmc, ghat, sigma, ikey       O     , dkm1, blmc, ghat, sigma, ikppkey
202       &     )       I     , myThid )
203    cph(
204  CADJ STORE dkm1,blmc,ghat = comlev1_kpp, key = ikey  CADJ STORE dkm1,blmc,ghat = comlev1_kpp, key = ikppkey
205    CADJ STORE hbl, kbl, diffus, casea = comlev1_kpp, key = ikppkey
206    cph)
207    
208  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
209  c enhance diffusivity at interface kbl - 1  c enhance diffusivity at interface kbl - 1
# Line 179  c--------------------------------------- Line 212  c---------------------------------------
212        call enhance (        call enhance (
213       I       dkm1, hbl, kbl, diffus, casea       I       dkm1, hbl, kbl, diffus, casea
214       U     , ghat       U     , ghat
215       O     , blmc )       O     , blmc
216         I     , myThid )
217    
218    cph(
219    cph avoids recomp. of enhance
220    CADJ STORE blmc = comlev1_kpp, key = ikppkey
221    cph)
222    
223  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
224  c combine interior and boundary layer coefficients and nonlocal term  c combine interior and boundary layer coefficients and nonlocal term
225    c !!!NOTE!!! In shallow (2-level) regions and for shallow mixed layers
226    c (< 1 level), diffusivity blmc can become negative.  The max-s below
227    c are a hack until this problem is properly diagnosed and fixed.
228  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
   
229        do k = 1, Nr        do k = 1, Nr
230           do i = 1, imt           do i = 1, imt
231              if (k .lt. kbl(i)) then              if (k .lt. kbl(i)) then
232                 do md = 1, mdiff  #ifdef ALLOW_SHELFICE
233                    diffus(i,k,md) = blmc(i,k,md)  C     when there is shelfice on top (msk(i)=0), reset the boundary layer
234                 end do  C     mixing coefficients blmc to pure Ri-number based mixing
235                   blmc(i,k,1) = max ( blmc(i,k,1)*msk(i),
236         &              diffus(i,k,1) )
237                   blmc(i,k,2) = max ( blmc(i,k,2)*msk(i),
238         &              diffus(i,k,2) )
239                   blmc(i,k,3) = max ( blmc(i,k,3)*msk(i),
240         &              diffus(i,k,3) )
241    #endif /* not ALLOW_SHELFICE */
242                   diffus(i,k,1) = max ( blmc(i,k,1), viscAr  )
243                   diffus(i,k,2) = max ( blmc(i,k,2), diffusKzS(i,Nr) )
244                   diffus(i,k,3) = max ( blmc(i,k,3), diffusKzT(i,Nr) )
245              else              else
246                 ghat(i,k) = 0.                 ghat(i,k) = 0. _d 0
247              endif              endif
248           end do           end do
249        end do        end do
# Line 205  c--------------------------------------- Line 256  c---------------------------------------
256  c*************************************************************************  c*************************************************************************
257    
258        subroutine bldepth (        subroutine bldepth (
259       I       mytime, mythid       I       kmtj
260       I     , kmtj       I     , dvsq, dbloc, Ritop, ustar, bo, bosol
261       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol  #ifdef ALLOW_SALT_PLUME
262       I     , ikey       I     , boplume,SPDepth
263    #endif /* ALLOW_SALT_PLUME */
264         I     , coriol
265         I     , ikppkey
266       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma
267       &     )       I     , myTime, myIter, myThid )
268    
269  c     the oceanic planetary boundary layer depth, hbl, is determined as  c     the oceanic planetary boundary layer depth, hbl, is determined as
270  c     the shallowest depth where the bulk Richardson number is  c     the shallowest depth where the bulk Richardson number is
# Line 239  c Line 293  c
293  #include "EEPARAMS.h"  #include "EEPARAMS.h"
294  #include "PARAMS.h"  #include "PARAMS.h"
295  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
 #include "FFIELDS.h"  
296    
297  c input  c input
298  c------  c------
299  c myTime    : current time in simulation  c   myTime :: Current time in simulation
300  c myThid    : thread number for this instance of the routine  c   myIter :: Current iteration number in simulation
301    c   myThid :: My Thread Id. number
302  c kmtj      : number of vertical layers  c kmtj      : number of vertical layers
303  c dvsq      : (velocity shear re sfc)^2             ((m/s)^2)  c dvsq      : (velocity shear re sfc)^2             ((m/s)^2)
304  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 308  c             buoyancy with respect to s
308  c ustar     : surface friction velocity                 (m/s)  c ustar     : surface friction velocity                 (m/s)
309  c bo        : surface turbulent buoyancy forcing    (m^2/s^3)  c bo        : surface turbulent buoyancy forcing    (m^2/s^3)
310  c bosol     : radiative buoyancy forcing            (m^2/s^3)  c bosol     : radiative buoyancy forcing            (m^2/s^3)
311    c boplume   : haline buoyancy forcing               (m^2/s^3)
312  c coriol    : Coriolis parameter                        (1/s)  c coriol    : Coriolis parameter                        (1/s)
313        _RL     mytime        _RL     myTime
314        integer mythid        integer myIter
315          integer myThid
316        integer kmtj(imt)        integer kmtj(imt)
317        _KPP_RL dvsq  (imt,Nr)        _RL dvsq    (imt,Nr)
318        _KPP_RL dbloc (imt,Nr)        _RL dbloc   (imt,Nr)
319        _KPP_RL Ritop (imt,Nr)        _RL Ritop   (imt,Nr)
320        _KPP_RL ustar (imt)        _RL ustar   (imt)
321        _KPP_RL bo    (imt)        _RL bo      (imt)
322        _KPP_RL bosol (imt)        _RL bosol   (imt)
323        _KPP_RL coriol(imt)        _RL coriol  (imt)
324        integer ikey        integer ikppkey, kkppkey
325    #ifdef ALLOW_SALT_PLUME
326          _RL boplume (imt)
327          _RL SPDepth (imt)
328    #endif /* ALLOW_SALT_PLUME */
329    
330  c  output  c  output
331  c--------  c--------
# Line 276  c casea     : =1 in case A, =0 in case B Line 336  c casea     : =1 in case A, =0 in case B
336  c kbl       : -1 of first grid level below hbl  c kbl       : -1 of first grid level below hbl
337  c Rib       : Bulk Richardson number  c Rib       : Bulk Richardson number
338  c sigma     : normalized depth (d/hbl)  c sigma     : normalized depth (d/hbl)
339        _KPP_RL hbl   (imt)        _RL hbl    (imt)
340        _KPP_RL bfsfc (imt)        _RL bfsfc  (imt)
341        _KPP_RL stable(imt)        _RL stable (imt)
342        _KPP_RL casea (imt)        _RL casea  (imt)
343        integer kbl   (imt)        integer kbl(imt)
344        _KPP_RL Rib   (imt,Nr)        _RL Rib    (imt,Nr)
345        _KPP_RL sigma (imt)        _RL sigma  (imt)
346    
347  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
348    
349  c  local  c  local
350  c-------  c-------
351  c wm, ws    : turbulent velocity scales         (m/s)  c wm, ws    : turbulent velocity scales         (m/s)
352        _KPP_RL wm(imt), ws(imt)        _RL wm(imt), ws(imt)
353        _RL     worka(imt)        _RL worka(imt)
354          _RL bvsq, vtsq, hekman, hmonob, hlimit, tempVar1, tempVar2
       _KPP_RL bvsq, vtsq, hekman, hmonob, hlimit, tempVar1, tempVar2  
355        integer i, kl        integer i, kl
356    
357        _KPP_RL     p5    , eins        _RL         p5    , eins
358        parameter ( p5=0.5, eins=1.0 )        parameter ( p5=0.5, eins=1.0 )
359        _RL         minusone        _RL         minusone
360        parameter ( minusone=-1.0 )        parameter ( minusone=-1.0 )
361    
362    #ifdef ALLOW_DIAGNOSTICS
363    c     KPPBFSFC - Bo+radiation absorbed to d=hbf*hbl + plume (m^2/s^3)
364          _RL KPPBFSFC(imt,Nr)
365          _RL KPPRi(imt,Nr)
366    #endif /* ALLOW_DIAGNOSTICS */
367    
368  c find bulk Richardson number at every grid level until > Ricr  c find bulk Richardson number at every grid level until > Ricr
369  c  c
370  c note: the reference depth is -epsilon/2.*zgrid(k), but the reference  c note: the reference depth is -epsilon/2.*zgrid(k), but the reference
# Line 312  c       kbl(i)=kmtj(i) and hbl(i)=-zgrid Line 377  c       kbl(i)=kmtj(i) and hbl(i)=-zgrid
377  c     initialize hbl and kbl to bottomed out values  c     initialize hbl and kbl to bottomed out values
378    
379        do i = 1, imt        do i = 1, imt
380           Rib(i,1) = 0.0           Rib(i,1) = 0. _d 0
381           kbl(i) = max(kmtj(i),1)           kbl(i) = max(kmtj(i),1)
382           hbl(i) = -zgrid(kbl(i))           hbl(i) = -zgrid(kbl(i))
383        end do        end do
384    
385    #ifdef ALLOW_DIAGNOSTICS
386          do kl = 1, Nr
387             do i = 1, imt
388                KPPBFSFC(i,kl) = 0. _d 0
389                KPPRi(i,kl) = 0. _d 0
390             enddo
391          enddo
392    #endif /* ALLOW_DIAGNOSTICS */
393    
394        do kl = 2, Nr        do kl = 2, Nr
395    
396    #ifdef ALLOW_AUTODIFF_TAMC
397             kkppkey = (ikppkey-1)*Nr + kl
398    #endif
399    
400  c     compute bfsfc = sw fraction at hbf * zgrid  c     compute bfsfc = sw fraction at hbf * zgrid
401    
402           do i = 1, imt           do i = 1, imt
403              worka(i) = zgrid(kl)              worka(i) = zgrid(kl)
404           end do           end do
405    CADJ store worka = comlev1_kpp_k, key = kkppkey
406           call SWFRAC(           call SWFRAC(
407       I        imt, hbf,       I       imt, hbf,
408       I        mytime, mythid,       U       worka,
409       U        worka )       I       myTime, myIter, myThid )
410    CADJ store worka = comlev1_kpp_k, key = kkppkey
411    
412           do i = 1, imt           do i = 1, imt
413    
# Line 338  c     use caseA as temporary array Line 418  c     use caseA as temporary array
418  c     compute bfsfc= Bo + radiative contribution down to hbf * hbl  c     compute bfsfc= Bo + radiative contribution down to hbf * hbl
419    
420              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  
421    
422           end do           end do
423    #ifdef ALLOW_SALT_PLUME
424    c     compute bfsfc = plume fraction at hbf * zgrid
425             IF ( useSALT_PLUME ) THEN
426               do i = 1, imt
427                  worka(i) = zgrid(kl)
428               enddo
429               call SALT_PLUME_FRAC(
430         I         imt, hbf,SPDepth,
431         U         worka,
432         I         myTime, myIter, myThid)
433               do i = 1, imt
434                  bfsfc(i) = bfsfc(i) + boplume(i)*(1. - worka(i))
435               enddo
436             ENDIF
437    #endif /* ALLOW_SALT_PLUME */
438    
439    #ifdef ALLOW_DIAGNOSTICS
440             do i = 1, imt
441                KPPBFSFC(i,kl) = bfsfc(i)
442             enddo
443    #endif /* ALLOW_DIAGNOSTICS */
444    
445             do i = 1, imt
446                stable(i) = p5 + sign(p5,bfsfc(i))
447                sigma(i) = stable(i) + (1. - stable(i)) * epsilon
448             enddo
449    
450  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
451  c     compute velocity scales at sigma, for hbl= caseA = -zgrid(kl)  c     compute velocity scales at sigma, for hbl= caseA = -zgrid(kl)
# Line 349  c--------------------------------------- Line 453  c---------------------------------------
453    
454           call wscale (           call wscale (
455       I        sigma, casea, ustar, bfsfc,       I        sigma, casea, ustar, bfsfc,
456       O        wm, ws )       O        wm, ws, myThid )
457    CADJ store ws = comlev1_kpp_k, key = kkppkey
458    
459           do i = 1, imt           do i = 1, imt
460    
# Line 361  c--------------------------------------- Line 466  c---------------------------------------
466       1           ( dbloc(i,kl-1) / (zgrid(kl-1)-zgrid(kl  ))+       1           ( dbloc(i,kl-1) / (zgrid(kl-1)-zgrid(kl  ))+
467       2             dbloc(i,kl  ) / (zgrid(kl  )-zgrid(kl+1)))       2             dbloc(i,kl  ) / (zgrid(kl  )-zgrid(kl+1)))
468    
469              if (bvsq .eq. 0.) then              if (bvsq .eq. 0. _d 0) then
470                vtsq = 0.0                vtsq = 0. _d 0
471              else              else
472                vtsq = -zgrid(kl) * ws(i) * sqrt(abs(bvsq)) * Vtc                vtsq = -zgrid(kl) * ws(i) * sqrt(abs(bvsq)) * Vtc
473              endif              endif
# Line 382  c Line 487  c
487              tempVar1  = dvsq(i,kl) + vtsq                        tempVar1  = dvsq(i,kl) + vtsq          
488              tempVar2 = max(tempVar1, phepsi)              tempVar2 = max(tempVar1, phepsi)
489              Rib(i,kl) = Ritop(i,kl) / tempVar2              Rib(i,kl) = Ritop(i,kl) / tempVar2
490    #ifdef ALLOW_DIAGNOSTICS
491                KPPRi(i,kl) = Rib(i,kl)
492    #endif
493    
494           end do           end do
495        end do        end do
496          
497    #ifdef ALLOW_DIAGNOSTICS
498             CALL DIAGNOSTICS_FILL(KPPBFSFC,'KPPbfsfc',0,Nr,0,1,1,myThid)
499             CALL DIAGNOSTICS_FILL(KPPRi   ,'KPPRi   ',0,Nr,0,1,1,myThid)
500    #endif /* ALLOW_DIAGNOSTICS */
501    
502    cph(
503    cph  without this store, there is a recomputation error for
504    cph  rib in adbldepth (probably partial recomputation problem)    
505    CADJ store Rib = comlev1_kpp
506    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy),Nr /)
507    cph)
508    
509        do kl = 2, Nr        do kl = 2, Nr
510           do i = 1, imt           do i = 1, imt
511              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 513  c
513        end do        end do
514    
515  CADJ store kbl = comlev1_kpp  CADJ store kbl = comlev1_kpp
516  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
517    
518        do i = 1, imt        do i = 1, imt
519           kl = kbl(i)           kl = kbl(i)
# Line 406  c     linearly interpolate to find hbl w Line 526  c     linearly interpolate to find hbl w
526        end do        end do
527    
528  CADJ store hbl = comlev1_kpp  CADJ store hbl = comlev1_kpp
529  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
530    
531  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
532  c     find stability and buoyancy forcing for boundary layer  c     find stability and buoyancy forcing for boundary layer
# Line 415  c--------------------------------------- Line 535  c---------------------------------------
535        do i = 1, imt        do i = 1, imt
536           worka(i) = hbl(i)           worka(i) = hbl(i)
537        end do        end do
538    CADJ store worka = comlev1_kpp
539    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
540        call SWFRAC(        call SWFRAC(
541       I     imt, minusone,       I       imt, minusone,
542       I     mytime, mythid,       U       worka,
543       U     worka )       I       myTime, myIter, myThid )
544    CADJ store worka = comlev1_kpp
545    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
546    
547        do i = 1, imt        do i = 1, imt
548           bfsfc(i)  = bo(i) + bosol(i) * (1. - worka(i))           bfsfc(i)  = bo(i) + bosol(i) * (1. - worka(i))
549        end do        end do
550    
551    #ifdef ALLOW_SALT_PLUME
552          IF ( useSALT_PLUME ) THEN
553            do i = 1, imt
554               worka(i) = hbl(i)
555            enddo
556            call SALT_PLUME_FRAC(
557         I         imt,minusone,SPDepth,
558         U         worka,
559         I         myTime, myIter, myThid )
560            do i = 1, imt
561               bfsfc(i) = bfsfc(i) + boplume(i) * (1. - worka(i))
562            enddo
563          ENDIF
564    #endif /* ALLOW_SALT_PLUME */
565  CADJ store bfsfc = comlev1_kpp  CADJ store bfsfc = comlev1_kpp
566  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
567    
568  c--   ensure bfsfc is never 0  c--   ensure bfsfc is never 0
569        do i = 1, imt        do i = 1, imt
# Line 432  c--   ensure bfsfc is never 0 Line 571  c--   ensure bfsfc is never 0
571           bfsfc(i) = sign(eins,bfsfc(i))*max(phepsi,abs(bfsfc(i)))           bfsfc(i) = sign(eins,bfsfc(i))*max(phepsi,abs(bfsfc(i)))
572        end do        end do
573    
574  CADJ store bfsfc = comlev1_kpp  cph(
575  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  cph  added stable to store list to avoid extensive recomp.
576    CADJ store bfsfc, stable = comlev1_kpp
577    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
578    cph)
579    
580  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
581  c check hbl limits for hekman or hmonob  c check hbl limits for hekman or hmonob
# Line 451  c--------------------------------------- Line 593  c---------------------------------------
593           end if           end if
594        end do        end do
595  CADJ store hbl = comlev1_kpp  CADJ store hbl = comlev1_kpp
596  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
597    
598        do i = 1, imt        do i = 1, imt
599           hbl(i) = max(hbl(i),minKPPhbl)           hbl(i) = max(hbl(i),minKPPhbl)
# Line 459  CADJ &   , key = ikey, shape = (/ (sNx+2 Line 601  CADJ &   , key = ikey, shape = (/ (sNx+2
601        end do        end do
602    
603  CADJ store hbl = comlev1_kpp  CADJ store hbl = comlev1_kpp
604  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
605    
606  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
607  c      find new kbl  c      find new kbl
# Line 480  c--------------------------------------- Line 622  c---------------------------------------
622        do i = 1, imt        do i = 1, imt
623           worka(i) = hbl(i)           worka(i) = hbl(i)
624        end do        end do
625    CADJ store worka = comlev1_kpp
626    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
627        call SWFRAC(        call SWFRAC(
628       I     imt, minusone,       I     imt, minusone,
629       I     mytime, mythid,       U     worka,
630       U     worka )       I     myTime, myIter, myThid )
631    CADJ store worka = comlev1_kpp
632    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
633    
634        do i = 1, imt        do i = 1, imt
635           bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i))           bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i))
636        end do        end do
637    
638    #ifdef ALLOW_SALT_PLUME
639          IF ( useSALT_PLUME ) THEN
640            do i = 1, imt
641               worka(i) = hbl(i)
642            enddo
643            call SALT_PLUME_FRAC(
644         I         imt,minusone,SPDepth,
645         U         worka,
646         I         myTime, myIter, myThid )
647            do i = 1, imt
648               bfsfc(i) = bfsfc(i) + boplume(i) * (1. - worka(i))
649            enddo
650          ENDIF
651    #endif /* ALLOW_SALT_PLUME */
652  CADJ store bfsfc = comlev1_kpp  CADJ store bfsfc = comlev1_kpp
653  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
654    
655  c--   ensures bfsfc is never 0  c--   ensures bfsfc is never 0
656        do i = 1, imt        do i = 1, imt
# Line 515  c*************************************** Line 676  c***************************************
676    
677        subroutine wscale (        subroutine wscale (
678       I     sigma, hbl, ustar, bfsfc,       I     sigma, hbl, ustar, bfsfc,
679       O     wm, ws )       O     wm, ws,
680         I     myThid )
681    
682  c     compute turbulent velocity scales.  c     compute turbulent velocity scales.
683  c     use a 2D-lookup table for wm and ws as functions of ustar and  c     use a 2D-lookup table for wm and ws as functions of ustar and
# Line 536  c sigma   : normalized depth (d/hbl) Line 698  c sigma   : normalized depth (d/hbl)
698  c hbl     : boundary layer depth (m)  c hbl     : boundary layer depth (m)
699  c ustar   : surface friction velocity         (m/s)  c ustar   : surface friction velocity         (m/s)
700  c bfsfc   : total surface buoyancy flux       (m^2/s^3)  c bfsfc   : total surface buoyancy flux       (m^2/s^3)
701        _KPP_RL sigma(imt)  c myThid  : thread number for this instance of the routine
702        _KPP_RL hbl  (imt)        integer myThid
703        _KPP_RL ustar(imt)        _RL sigma(imt)
704        _KPP_RL bfsfc(imt)        _RL hbl  (imt)
705          _RL ustar(imt)
706          _RL bfsfc(imt)
707    
708  c  output  c  output
709  c--------  c--------
710  c wm, ws  : turbulent velocity scales at sigma  c wm, ws  : turbulent velocity scales at sigma
711        _KPP_RL wm(imt), ws(imt)        _RL wm(imt), ws(imt)
712    
713  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
714    
715  c local  c local
716  c------  c------
717  c zehat   : = zeta *  ustar**3  c zehat   : = zeta *  ustar**3
718        _KPP_RL zehat        _RL zehat
719    
720        integer iz, izp1, ju, i, jup1        integer iz, izp1, ju, i, jup1
721        _KPP_RL udiff, zdiff, zfrac, ufrac, fzfrac, wam        _RL udiff, zdiff, zfrac, ufrac, fzfrac, wam
722        _KPP_RL wbm, was, wbs, u3, tempVar        _RL wbm, was, wbs, u3, tempVar
723    
724  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
725  c use lookup table for zehat < zmax only; otherwise use  c use lookup table for zehat < zmax only; otherwise use
# Line 610  c--------------------------------------- Line 774  c---------------------------------------
774  c*************************************************************************  c*************************************************************************
775    
776        subroutine Ri_iwmix (        subroutine Ri_iwmix (
777       I       kmtj, shsq, dbloc, dblocSm       I       kmtj, shsq, dbloc, dblocSm,
778       I     , ikey       I       diffusKzS, diffusKzT,
779       O     , diffus )       I       ikppkey,
780         O       diffus,
781         I       myThid )
782    
783  c     compute interior viscosity diffusivity coefficients due  c     compute interior viscosity diffusivity coefficients due
784  c     to shear instability (dependent on a local Richardson number),  c     to shear instability (dependent on a local Richardson number),
# Line 631  c     kmtj   (imt)         number of ver Line 797  c     kmtj   (imt)         number of ver
797  c     shsq   (imt,Nr)      (local velocity shear)^2               ((m/s)^2)  c     shsq   (imt,Nr)      (local velocity shear)^2               ((m/s)^2)
798  c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)  c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)
799  c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)  c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)
800        integer kmtj   (imt)  c     diffusKzS(imt,Nr)- background vertical diffusivity for scalars    (m^2/s)
801        _KPP_RL shsq   (imt,Nr)  c     diffusKzT(imt,Nr)- background vertical diffusivity for theta      (m^2/s)
802        _KPP_RL dbloc  (imt,Nr)  c     myThid :: My Thread Id. number
803        _KPP_RL dblocSm(imt,Nr)        integer kmtj (imt)
804        integer ikey        _RL shsq     (imt,Nr)
805          _RL dbloc    (imt,Nr)
806          _RL dblocSm  (imt,Nr)
807          _RL diffusKzS(imt,Nr)
808          _RL diffusKzT(imt,Nr)
809          integer ikppkey
810          integer myThid
811    
812  c output  c output
813  c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)  c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)
814  c     diffus(imt,0:Nrp1,2)  vertical scalar diffusivity             (m^2/s)  c     diffus(imt,0:Nrp1,2)  vertical scalar diffusivity             (m^2/s)
815  c     diffus(imt,0:Nrp1,3)  vertical temperature diffusivity        (m^2/s)  c     diffus(imt,0:Nrp1,3)  vertical temperature diffusivity        (m^2/s)
816        _KPP_RL diffus(imt,0:Nrp1,3)        _RL diffus(imt,0:Nrp1,3)
817    
818  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
819    
820  c local variables  c local variables
821  c     Rig                   local Richardson number  c     Rig                   local Richardson number
822  c     fRi, fcon             function of Rig  c     fRi, fcon             function of Rig
823        _KPP_RL Rig        _RL Rig
824        _KPP_RL fRi, fcon        _RL fRi, fcon
825        _KPP_RL ratio        _RL ratio
826        integer i, ki, mr        integer i, ki
827        _KPP_RL c1, c0        _RL c1, c0
828    
829  #ifdef ALLOW_KPP_VERTICALLY_SMOOTH  #ifdef ALLOW_KPP_VERTICALLY_SMOOTH
830          integer mr
831  CADJ INIT kpp_ri_tape_mr = common, 1  CADJ INIT kpp_ri_tape_mr = common, 1
832  #endif  #endif
833    
834  c constants  c constants
835        c1 = 1.0        c1 = 1. _d 0
836        c0 = 0.0        c0 = 0. _d 0
837    
838  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
839  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 843  c     set values at bottom and below to
843  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
844  C     break data flow dependence on diffus  C     break data flow dependence on diffus
845        diffus(1,1,1) = 0.0        diffus(1,1,1) = 0.0
846    
847          do ki = 1, Nr
848             do i = 1, imt
849                diffus(i,ki,1) = 0.
850                diffus(i,ki,2) = 0.
851                diffus(i,ki,3) = 0.
852             enddo
853          enddo
854  #endif  #endif
855                
856    
857        do ki = 1, Nr        do ki = 1, Nr
858           do i = 1, imt           do i = 1, imt
859              if     (kmtj(i) .EQ. 0      ) then              if     (kmtj(i) .LE. 1      ) then
860                 diffus(i,ki,1) = 0.                 diffus(i,ki,1) = 0.
861                 diffus(i,ki,2) = 0.                 diffus(i,ki,2) = 0.
862              elseif (ki      .GE. kmtj(i)) then              elseif (ki      .GE. kmtj(i)) then
# Line 687  C     break data flow dependence on diff Line 869  C     break data flow dependence on diff
869              endif              endif
870           end do           end do
871        end do        end do
872    CADJ store diffus = comlev1_kpp, key = ikppkey
873    
874  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
875  c     vertically smooth Ri  c     vertically smooth Ri
# Line 697  CADJ store diffus(:,:,1) = kpp_ri_tape_m Line 880  CADJ store diffus(:,:,1) = kpp_ri_tape_m
880  CADJ &  , key=mr, shape=(/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2 /)  CADJ &  , key=mr, shape=(/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2 /)
881    
882           call z121 (           call z121 (
883       U     diffus(1,0,1))       U     diffus(1,0,1),
884         I     myThid )
885        end do        end do
886  #endif  #endif
887    
# Line 724  c  evaluate f of smooth Ri for shear ins Line 908  c  evaluate f of smooth Ri for shear ins
908  c ----------------------------------------------------------------------  c ----------------------------------------------------------------------
909  c            evaluate diffusivities and viscosity  c            evaluate diffusivities and viscosity
910  c    mixing due to internal waves, and shear and static instability  c    mixing due to internal waves, and shear and static instability
911    
912              diffus(i,ki,1) = viscAr + fcon * difmcon + fRi * difm0  #ifndef EXCLUDE_KPP_SHEAR_MIX
913              diffus(i,ki,2) = diffKrS + fcon * difscon + fRi * difs0              if ( .NOT. inAdMode ) then
914              diffus(i,ki,3) = diffKrT + fcon * difscon + fRi * difs0                 diffus(i,ki,1) = viscAr + fcon * difmcon + fRi * difm0
915                   diffus(i,ki,2) = diffusKzS(i,ki)+fcon*difscon+fRi*difs0
916                   diffus(i,ki,3) = diffusKzT(i,ki)+fcon*diftcon+fRi*dift0
917                else
918                   diffus(i,ki,1) = viscAr
919                   diffus(i,ki,2) = diffusKzS(i,ki)
920                   diffus(i,ki,3) = diffusKzT(i,ki)
921                endif
922    #else
923                diffus(i,ki,1) = viscAr
924                diffus(i,ki,2) = diffusKzS(i,ki)
925                diffus(i,ki,3) = diffusKzT(i,ki)
926    #endif
927    
928           end do           end do
929        end do        end do
# Line 749  c         set surface values to 0.0 Line 945  c         set surface values to 0.0
945  c*************************************************************************  c*************************************************************************
946    
947        subroutine z121 (        subroutine z121 (
948       U     v )       U     v,
949         I     myThid )
950    
951  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)
952  c     top (0) value is used as a dummy  c     top (0) value is used as a dummy
# Line 767  c     penetrative convection. Line 964  c     penetrative convection.
964  c input/output  c input/output
965  c-------------  c-------------
966  c v     : 2-D array to be smoothed in Nrp1 direction  c v     : 2-D array to be smoothed in Nrp1 direction
967        _KPP_RL v(imt,0:Nrp1)  c myThid: thread number for this instance of the routine
968          integer myThid
969          _RL v(imt,0:Nrp1)
970    
971  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
972    
973  c local  c local
974        _KPP_RL zwork, zflag        _RL zwork, zflag
975        _KPP_RL KRi_range(1:Nrp1)        _RL KRi_range(1:Nrp1)
976        integer i, k, km1, kp1        integer i, k, km1, kp1
977    
978        _KPP_RL     p0      , p25       , p5      , p2        _RL         p0      , p25       , p5      , p2
979        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 )
980    
981        KRi_range(Nrp1) = p0        KRi_range(Nrp1) = p0
# Line 830  CADJ STORE v(i,k), zwork = z121tape Line 1029  CADJ STORE v(i,k), zwork = z121tape
1029    
1030  c*************************************************************************  c*************************************************************************
1031    
       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*************************************************************************  
   
1032        subroutine smooth_horiz (        subroutine smooth_horiz (
1033       I     k, bi, bj,       I     k, bi, bj,
1034       U     fld )       U     fld,
1035         I     myThid )
1036    
1037  c     Apply horizontal smoothing to global _RL 2-D array  c     Apply horizontal smoothing to global _RL 2-D array
1038    
1039        IMPLICIT NONE        IMPLICIT NONE
1040  #include "SIZE.h"  #include "SIZE.h"
1041    #include "GRID.h"
1042  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
1043    
1044  c     input  c     input
1045  c     bi, bj : array indices  c     bi, bj : array indices
1046  c     k      : vertical index used for masking  c     k      : vertical index used for masking
1047    c     myThid : thread number for this instance of the routine
1048          INTEGER myThid
1049        integer k, bi, bj        integer k, bi, bj
1050    
1051  c     input/output  c     input/output
# Line 949  c     local Line 1072  c     local
1072              im1 = i-1              im1 = i-1
1073              ip1 = i+1              ip1 = i+1
1074              tempVar =              tempVar =
1075       &           p25   *   pMask(i  ,j  ,k,bi,bj)   +       &           p25   *   maskC(i  ,j  ,k,bi,bj)   +
1076       &           p125  * ( pMask(im1,j  ,k,bi,bj)   +       &           p125  * ( maskC(im1,j  ,k,bi,bj)   +
1077       &                     pMask(ip1,j  ,k,bi,bj)   +       &                     maskC(ip1,j  ,k,bi,bj)   +
1078       &                     pMask(i  ,jm1,k,bi,bj)   +       &                     maskC(i  ,jm1,k,bi,bj)   +
1079       &                     pMask(i  ,jp1,k,bi,bj) ) +       &                     maskC(i  ,jp1,k,bi,bj) ) +
1080       &           p0625 * ( pMask(im1,jm1,k,bi,bj)   +       &           p0625 * ( maskC(im1,jm1,k,bi,bj)   +
1081       &                     pMask(im1,jp1,k,bi,bj)   +       &                     maskC(im1,jp1,k,bi,bj)   +
1082       &                     pMask(ip1,jm1,k,bi,bj)   +       &                     maskC(ip1,jm1,k,bi,bj)   +
1083       &                     pMask(ip1,jp1,k,bi,bj) )       &                     maskC(ip1,jp1,k,bi,bj) )
1084              IF ( tempVar .GE. p25 ) THEN              IF ( tempVar .GE. p25 ) THEN
1085                 fld_tmp(i,j) = (                 fld_tmp(i,j) = (
1086       &              p25  * fld(i  ,j  )*pMask(i  ,j  ,k,bi,bj) +       &              p25  * fld(i  ,j  )*maskC(i  ,j  ,k,bi,bj) +
1087       &              p125 *(fld(im1,j  )*pMask(im1,j  ,k,bi,bj) +       &              p125 *(fld(im1,j  )*maskC(im1,j  ,k,bi,bj) +
1088       &                     fld(ip1,j  )*pMask(ip1,j  ,k,bi,bj) +       &                     fld(ip1,j  )*maskC(ip1,j  ,k,bi,bj) +
1089       &                     fld(i  ,jm1)*pMask(i  ,jm1,k,bi,bj) +       &                     fld(i  ,jm1)*maskC(i  ,jm1,k,bi,bj) +
1090       &                     fld(i  ,jp1)*pMask(i  ,jp1,k,bi,bj))+       &                     fld(i  ,jp1)*maskC(i  ,jp1,k,bi,bj))+
1091       &              p0625*(fld(im1,jm1)*pMask(im1,jm1,k,bi,bj) +       &              p0625*(fld(im1,jm1)*maskC(im1,jm1,k,bi,bj) +
1092       &                     fld(im1,jp1)*pMask(im1,jp1,k,bi,bj) +       &                     fld(im1,jp1)*maskC(im1,jp1,k,bi,bj) +
1093       &                     fld(ip1,jm1)*pMask(ip1,jm1,k,bi,bj) +       &                     fld(ip1,jm1)*maskC(ip1,jm1,k,bi,bj) +
1094       &                     fld(ip1,jp1)*pMask(ip1,jp1,k,bi,bj)))       &                     fld(ip1,jp1)*maskC(ip1,jp1,k,bi,bj)))
1095       &              / tempVar       &              / tempVar
1096              ELSE              ELSE
1097                 fld_tmp(i,j) = fld(i,j)                 fld_tmp(i,j) = fld(i,j)
# Line 992  c*************************************** Line 1115  c***************************************
1115    
1116        subroutine blmix (        subroutine blmix (
1117       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl
1118       O     , dkm1, blmc, ghat, sigma, ikey       O     , dkm1, blmc, ghat, sigma, ikppkey
1119       &     )       I     , myThid )
1120    
1121  c     mixing coefficients within boundary layer depend on surface  c     mixing coefficients within boundary layer depend on surface
1122  c     forcing and the magnitude and gradient of interior mixing below  c     forcing and the magnitude and gradient of interior mixing below
# Line 1015  c     hbl   (imt)                 bounda Line 1138  c     hbl   (imt)                 bounda
1138  c     stable(imt)                 = 1 in stable forcing  c     stable(imt)                 = 1 in stable forcing
1139  c     casea (imt)                 = 1 in case A  c     casea (imt)                 = 1 in case A
1140  c     diffus(imt,0:Nrp1,mdiff)    vertical diffusivities              (m^2/s)  c     diffus(imt,0:Nrp1,mdiff)    vertical diffusivities              (m^2/s)
1141  c     kbl(imt)                    -1 of first grid level below hbl  c     kbl   (imt)                 -1 of first grid level below hbl
1142        _KPP_RL ustar (imt)  c     myThid               thread number for this instance of the routine
1143        _KPP_RL bfsfc (imt)        integer myThid
1144        _KPP_RL hbl   (imt)        _RL ustar (imt)
1145        _KPP_RL stable(imt)        _RL bfsfc (imt)
1146        _KPP_RL casea (imt)        _RL hbl   (imt)
1147        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL stable(imt)
1148          _RL casea (imt)
1149          _RL diffus(imt,0:Nrp1,mdiff)
1150        integer kbl(imt)        integer kbl(imt)
1151    
1152  c output  c output
# Line 1029  c     dkm1 (imt,mdiff)            bounda Line 1154  c     dkm1 (imt,mdiff)            bounda
1154  c     blmc (imt,Nr,mdiff)         boundary layer mixing coefficients  (m^2/s)  c     blmc (imt,Nr,mdiff)         boundary layer mixing coefficients  (m^2/s)
1155  c     ghat (imt,Nr)               nonlocal scalar transport  c     ghat (imt,Nr)               nonlocal scalar transport
1156  c     sigma(imt)                  normalized depth (d / hbl)  c     sigma(imt)                  normalized depth (d / hbl)
1157        _KPP_RL dkm1 (imt,mdiff)        _RL dkm1 (imt,mdiff)
1158        _KPP_RL blmc (imt,Nr,mdiff)        _RL blmc (imt,Nr,mdiff)
1159        _KPP_RL ghat (imt,Nr)        _RL ghat (imt,Nr)
1160        _KPP_RL sigma(imt)        _RL sigma(imt)
1161        integer ikey        integer ikppkey, kkppkey
1162    
1163  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1164    
# Line 1041  c  local Line 1166  c  local
1166  c     gat1*(imt)                 shape function at sigma = 1  c     gat1*(imt)                 shape function at sigma = 1
1167  c     dat1*(imt)                 derivative of shape function at sigma = 1  c     dat1*(imt)                 derivative of shape function at sigma = 1
1168  c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)  c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)
1169        _KPP_RL gat1m(imt), gat1s(imt), gat1t(imt)        _RL gat1m(imt), gat1s(imt), gat1t(imt)
1170        _KPP_RL dat1m(imt), dat1s(imt), dat1t(imt)        _RL dat1m(imt), dat1s(imt), dat1t(imt)
1171        _KPP_RL ws(imt), wm(imt)        _RL ws(imt), wm(imt)
1172        integer i, kn, ki        integer i, kn, ki
1173        _KPP_RL R, dvdzup, dvdzdn, viscp        _RL R, dvdzup, dvdzdn, viscp
1174        _KPP_RL difsp, diftp, visch, difsh, difth        _RL difsp, diftp, visch, difsh, difth
1175        _KPP_RL f1, sig, a1, a2, a3, delhat        _RL f1, sig, a1, a2, a3, delhat
1176        _KPP_RL Gm, Gs, Gt        _RL Gm, Gs, Gt
1177        _KPP_RL tempVar        _RL tempVar
1178    
1179        _KPP_RL    p0    , eins        _RL    p0    , eins
1180        parameter (p0=0.0, eins=1.0)        parameter (p0=0.0, eins=1.0)
1181    
1182  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 1062  c--------------------------------------- Line 1187  c---------------------------------------
1187           sigma(i) = stable(i) * 1.0 + (1. - stable(i)) * epsilon           sigma(i) = stable(i) * 1.0 + (1. - stable(i)) * epsilon
1188        end do        end do
1189    
1190    CADJ STORE sigma = comlev1_kpp, key = ikppkey
1191        call wscale (        call wscale (
1192       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1193       O        wm, ws )       O        wm, ws, myThid )
1194    CADJ STORE wm = comlev1_kpp, key = ikppkey
1195    CADJ STORE ws = comlev1_kpp, key = ikppkey
1196    
1197        do i = 1, imt        do i = 1, imt
1198           wm(i) = sign(eins,wm(i))*max(phepsi,abs(wm(i)))           wm(i) = sign(eins,wm(i))*max(phepsi,abs(wm(i)))
1199           ws(i) = sign(eins,ws(i))*max(phepsi,abs(ws(i)))           ws(i) = sign(eins,ws(i))*max(phepsi,abs(ws(i)))
1200        end do        end do
1201  CADJ STORE wm = comlev1_kpp, key = ikey  CADJ STORE wm = comlev1_kpp, key = ikppkey
1202  CADJ STORE ws = comlev1_kpp, key = ikey  CADJ STORE ws = comlev1_kpp, key = ikppkey
1203    
1204        do i = 1, imt        do i = 1, imt
1205    
# Line 1107  c--------------------------------------- Line 1235  c---------------------------------------
1235       &        max(ustar(i)**4,phepsi)       &        max(ustar(i)**4,phepsi)
1236           gat1m(i) = visch / hbl(i) / wm(i)           gat1m(i) = visch / hbl(i) / wm(i)
1237           dat1m(i) = -viscp / wm(i) + f1 * visch           dat1m(i) = -viscp / wm(i) + f1 * visch
          dat1m(i) = min(dat1m(i),p0)  
1238    
1239           gat1s(i) = difsh  / hbl(i) / ws(i)           gat1s(i) = difsh  / hbl(i) / ws(i)
1240           dat1s(i) = -difsp / ws(i) + f1 * difsh           dat1s(i) = -difsp / ws(i) + f1 * difsh
          dat1s(i) = min(dat1s(i),p0)  
1241    
1242           gat1t(i) = difth /  hbl(i) / ws(i)           gat1t(i) = difth /  hbl(i) / ws(i)
1243           dat1t(i) = -diftp / ws(i) + f1 * difth           dat1t(i) = -diftp / ws(i) + f1 * difth
          dat1t(i) = min(dat1t(i),p0)  
1244    
1245        end do        end do
1246    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1247    CADJ STORE gat1m = comlev1_kpp, key = ikppkey
1248    CADJ STORE gat1s = comlev1_kpp, key = ikppkey
1249    CADJ STORE gat1t = comlev1_kpp, key = ikppkey
1250    CADJ STORE dat1m = comlev1_kpp, key = ikppkey
1251    CADJ STORE dat1s = comlev1_kpp, key = ikppkey
1252    CADJ STORE dat1t = comlev1_kpp, key = ikppkey
1253    #endif
1254          do i = 1, imt
1255             dat1m(i) = min(dat1m(i),p0)
1256             dat1s(i) = min(dat1s(i),p0)
1257             dat1t(i) = min(dat1t(i),p0)
1258          end do
1259    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1260    CADJ STORE dat1m = comlev1_kpp, key = ikppkey
1261    CADJ STORE dat1s = comlev1_kpp, key = ikppkey
1262    CADJ STORE dat1t = comlev1_kpp, key = ikppkey
1263    #endif
1264    
1265        do ki = 1, Nr        do ki = 1, Nr
1266    
1267    #ifdef ALLOW_AUTODIFF_TAMC
1268             kkppkey = (ikppkey-1)*Nr + ki
1269    #endif
1270    
1271  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1272  c     compute turbulent velocity scales on the interfaces  c     compute turbulent velocity scales on the interfaces
1273  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 1129  c--------------------------------------- Line 1276  c---------------------------------------
1276              sig      = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i)              sig      = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i)
1277              sigma(i) = stable(i)*sig + (1.-stable(i))*min(sig,epsilon)              sigma(i) = stable(i)*sig + (1.-stable(i))*min(sig,epsilon)
1278           end do           end do
1279    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1280    CADJ STORE wm = comlev1_kpp_k, key = kkppkey
1281    CADJ STORE ws = comlev1_kpp_k, key = kkppkey
1282    #endif
1283    CADJ STORE sigma = comlev1_kpp_k, key = kkppkey
1284           call wscale (           call wscale (
1285       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1286       O        wm, ws )       O        wm, ws, myThid )
1287    CADJ STORE wm = comlev1_kpp_k, key = kkppkey
1288    CADJ STORE ws = comlev1_kpp_k, key = kkppkey
1289    
1290  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1291  c     compute the dimensionless shape functions at the interfaces  c     compute the dimensionless shape functions at the interfaces
# Line 1175  c--------------------------------------- Line 1329  c---------------------------------------
1329       &            + (1. - stable(i)) * min(sig,epsilon)       &            + (1. - stable(i)) * min(sig,epsilon)
1330        end do        end do
1331    
1332    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1333    CADJ STORE wm = comlev1_kpp, key = ikppkey
1334    CADJ STORE ws = comlev1_kpp, key = ikppkey
1335    #endif
1336    CADJ STORE sigma = comlev1_kpp, key = ikppkey
1337        call wscale (        call wscale (
1338       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1339       O        wm, ws )       O        wm, ws, myThid )
1340    CADJ STORE wm = comlev1_kpp, key = ikppkey
1341    CADJ STORE ws = comlev1_kpp, key = ikppkey
1342    
1343        do i = 1, imt        do i = 1, imt
1344           sig = -zgrid(kbl(i)-1) / hbl(i)           sig = -zgrid(kbl(i)-1) / hbl(i)
# Line 1203  c*************************************** Line 1364  c***************************************
1364       I       dkm1, hbl, kbl, diffus, casea       I       dkm1, hbl, kbl, diffus, casea
1365       U     , ghat       U     , ghat
1366       O     , blmc       O     , blmc
1367       &     )       &     , myThid )
1368    
1369  c enhance the diffusivity at the kbl-.5 interface  c enhance the diffusivity at the kbl-.5 interface
1370    
# Line 1218  c     hbl(imt)                  boundary Line 1379  c     hbl(imt)                  boundary
1379  c     kbl(imt)                  grid above hbl  c     kbl(imt)                  grid above hbl
1380  c     diffus(imt,0:Nrp1,mdiff) vertical diffusivities           (m^2/s)  c     diffus(imt,0:Nrp1,mdiff) vertical diffusivities           (m^2/s)
1381  c     casea(imt)                = 1 in caseA, = 0 in case B  c     casea(imt)                = 1 in caseA, = 0 in case B
1382        _KPP_RL dkm1  (imt,mdiff)  c     myThid                    thread number for this instance of the routine
1383        _KPP_RL hbl   (imt)        integer   myThid
1384          _RL dkm1  (imt,mdiff)
1385          _RL hbl   (imt)
1386        integer kbl   (imt)        integer kbl   (imt)
1387        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL diffus(imt,0:Nrp1,mdiff)
1388        _KPP_RL casea (imt)        _RL casea (imt)
1389    
1390  c input/output  c input/output
1391  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)
1392        _KPP_RL ghat (imt,Nr)        _RL ghat (imt,Nr)
1393    
1394  c output  c output
1395  c     enhanced bound. layer mixing coeff.  c     enhanced bound. layer mixing coeff.
1396        _KPP_RL blmc  (imt,Nr,mdiff)        _RL blmc  (imt,Nr,mdiff)
1397    
1398  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1399    
1400  c local  c local
1401  c     fraction hbl lies beteen zgrid neighbors  c     fraction hbl lies beteen zgrid neighbors
1402        _KPP_RL delta        _RL delta
1403        integer ki, i, md        integer ki, i, md
1404        _KPP_RL dkmp5, dstar        _RL dkmp5, dstar
1405    
1406        do i = 1, imt        do i = 1, imt
1407           ki = kbl(i)-1           ki = kbl(i)-1
# Line 1264  c     fraction hbl lies beteen zgrid nei Line 1427  c     fraction hbl lies beteen zgrid nei
1427  c*************************************************************************  c*************************************************************************
1428    
1429        SUBROUTINE STATEKPP (        SUBROUTINE STATEKPP (
1430       I     bi, bj, myThid,       O     RHO1, DBLOC, DBSFC, TTALPHA, SSBETA,
1431       O     RHO1, DBLOC, DBSFC, TTALPHA, SSBETA)       I     ikppkey, bi, bj, myThid )
1432  c  c
1433  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1434  c     "statekpp" computes all necessary input arrays  c     "statekpp" computes all necessary input arrays
# Line 1290  c Line 1453  c
1453  c     written  by: jan morzel,   feb. 10, 1995 (converted from "sigma" version)  c     written  by: jan morzel,   feb. 10, 1995 (converted from "sigma" version)
1454  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
1455  c  c
1456    c     28 april 05: added computation of mixed layer depth KPPmld
1457    c     for a deltaRho equivalent to a deltaTheta of 0.8 deg C
1458    
1459  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1460    
1461        IMPLICIT NONE        IMPLICIT NONE
# Line 1299  c--------------------------------------- Line 1465  c---------------------------------------
1465  #include "PARAMS.h"  #include "PARAMS.h"
1466  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
1467  #include "DYNVARS.h"  #include "DYNVARS.h"
1468    #include "GRID.h"
1469    
1470  c-------------- Routine arguments -----------------------------------------  c-------------- Routine arguments -----------------------------------------
1471        INTEGER bi, bj, myThid        INTEGER bi, bj, myThid
1472        _KPP_RL RHO1   ( ibot:itop, jbot:jtop       )        _RL RHO1   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy       )
1473        _KPP_RL DBLOC  ( ibot:itop, jbot:jtop, Nr   )        _RL DBLOC  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
1474        _KPP_RL DBSFC  ( ibot:itop, jbot:jtop, Nr   )        _RL DBSFC  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
1475        _KPP_RL TTALPHA( ibot:itop, jbot:jtop, Nrp1 )        _RL TTALPHA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1476        _KPP_RL SSBETA ( ibot:itop, jbot:jtop, Nrp1 )        _RL SSBETA ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1477    
1478  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1479    
# Line 1317  c Line 1484  c
1484  c     rhok         - density of t(k  ) & s(k  ) at depth k  c     rhok         - density of t(k  ) & s(k  ) at depth k
1485  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
1486  c     rho1k        - density of t(1  ) & s(1  ) at depth k  c     rho1k        - density of t(1  ) & s(1  ) at depth k
1487  c     work1, work2 - work arrays for holding horizontal slabs  c     work1,2,3    - work arrays for holding horizontal slabs
1488    
1489        _RL RHOK  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL RHOK  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1490        _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 1492  c     work1, work2 - work arrays for hol
1492        _RL WORK1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1493        _RL WORK2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1494        _RL WORK3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1495    #ifdef ALLOW_DIAGNOSTICS
1496    c     KPPMLD - mixed layer depth based on density criterion
1497          _RL KPPMLD(1    :sNx    ,1    :sNy    )
1498    #endif /* ALLOW_DIAGNOSTICS */
1499    
1500        INTEGER I, J, K        INTEGER I, J, K
1501          INTEGER ikppkey, kkppkey
1502    
1503  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
1504    
1505        call FIND_RHO(        kkppkey = (ikppkey-1)*Nr + 1
1506       I     bi, bj, ibot, itop, jbot, jtop, 1, 1, eosType,  
1507       I     theta, salt,  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1508    CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1509    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1510    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1511          CALL FIND_RHO_2D(
1512         I     1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1,
1513         I     theta(1-OLx,1-OLy,1,bi,bj), salt(1-OLx,1-OLy,1,bi,bj),
1514       O     WORK1,       O     WORK1,
1515       I     myThid )       I     1, bi, bj, myThid )
1516    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1517    CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1518    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1519    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1520    
1521        call FIND_ALPHA(        call FIND_ALPHA(
1522       I     bi, bj, ibot, itop, jbot, jtop, 1, 1, eosType,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1523       O     WORK2 )       O     WORK2, myThid )
1524    
1525        call FIND_BETA(        call FIND_BETA(
1526       I     bi, bj, ibot, itop, jbot, jtop, 1, 1, eosType,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1527       O     WORK3 )       O     WORK3, myThid )
1528    
1529        DO J = jbot, jtop        DO J = 1-OLy, sNy+OLy
1530           DO I = ibot, itop           DO I = 1-OLx, sNx+OLx
1531              RHO1(I,J)      = WORK1(I,J) + rhonil              RHO1(I,J)      = WORK1(I,J) + rhoConst
1532              TTALPHA(I,J,1) = WORK2(I,J)              TTALPHA(I,J,1) = WORK2(I,J)
1533              SSBETA(I,J,1)  = WORK3(I,J)              SSBETA(I,J,1)  = WORK3(I,J)
1534              DBSFC(I,J,1)   = 0.              DBSFC(I,J,1)   = 0.
1535           END DO           END DO
1536        END DO        END DO
1537    
1538    #ifdef ALLOW_DIAGNOSTICS
1539    c     work3 - density of t(1)-.8 & s(1  ) at depth 1
1540          IF ( useDiagnostics ) THEN
1541             DO J = 1, sNy
1542                DO I = 1, sNx
1543                   KPPMLD(I,J) = ABS(R_low(I,J,bi,bj))
1544                   WORK3 (I,J) = WORK1(I,J) - 0.8 _d 0 * WORK2(I,J)
1545                END DO
1546             END DO
1547          ENDIF
1548    #endif /* ALLOW_DIAGNOSTICS */
1549    
1550  c calculate alpha, beta, and gradients in interior layers  c calculate alpha, beta, and gradients in interior layers
1551    
1552  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)
1553        DO K = 2, Nr        DO K = 2, Nr
1554    
1555           call FIND_RHO(        kkppkey = (ikppkey-1)*Nr + k
1556       I        bi, bj, ibot, itop, jbot, jtop, K, K, eosType,  
1557       I        theta, salt,  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1558    CADJ STORE theta(:,:,k,bi,bj) = comlev1_kpp_k, key=kkppkey
1559    CADJ STORE salt (:,:,k,bi,bj) = comlev1_kpp_k, key=kkppkey
1560    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1561             CALL FIND_RHO_2D(
1562         I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, k,
1563         I        theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
1564       O        RHOK,       O        RHOK,
1565       I        myThid )       I        k, bi, bj, myThid )
1566    
1567           call FIND_RHO(  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1568       I        bi, bj, ibot, itop, jbot, jtop, K-1, K, eosType,  CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_kpp_k, key=kkppkey
1569       I        theta, salt,  CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_kpp_k, key=kkppkey
1570    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1571             CALL FIND_RHO_2D(
1572         I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, k,
1573         I        theta(1-OLx,1-OLy,k-1,bi,bj),salt(1-OLx,1-OLy,k-1,bi,bj),
1574       O        RHOKM1,       O        RHOKM1,
1575       I        myThid )       I        k-1, bi, bj, myThid )
1576    
1577           call FIND_RHO(  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1578       I        bi, bj, ibot, itop, jbot, jtop, 1, K, eosType,  CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1579       I        theta, salt,  CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1580    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1581             CALL FIND_RHO_2D(
1582         I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, k,
1583         I        theta(1-OLx,1-OLy,1,bi,bj), salt(1-OLx,1-OLy,1,bi,bj),
1584       O        RHO1K,       O        RHO1K,
1585       I        myThid )       I        1, bi, bj, myThid )
1586    
1587    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1588    CADJ STORE rhok  (:,:)          = comlev1_kpp_k, key=kkppkey
1589    CADJ STORE rhokm1(:,:)          = comlev1_kpp_k, key=kkppkey
1590    CADJ STORE rho1k (:,:)          = comlev1_kpp_k, key=kkppkey
1591    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1592    
1593           call FIND_ALPHA(           call FIND_ALPHA(
1594       I        bi, bj, ibot, itop, jbot, jtop, K, K, eosType,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1595       O        WORK1 )       O        WORK1, myThid )
1596    
1597           call FIND_BETA(           call FIND_BETA(
1598       I        bi, bj, ibot, itop, jbot, jtop, K, K, eosType,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1599       O        WORK2 )       O        WORK2, myThid )
1600    
1601           DO J = jbot, jtop           DO J = 1-OLy, sNy+OLy
1602              DO I = ibot, itop              DO I = 1-OLx, sNx+OLx
1603                 TTALPHA(I,J,K) = WORK1 (I,J)                 TTALPHA(I,J,K) = WORK1 (I,J)
1604                 SSBETA(I,J,K)  = WORK2 (I,J)                 SSBETA(I,J,K)  = WORK2 (I,J)
1605                 DBLOC(I,J,K-1) = gravity * (RHOK(I,J) - RHOKM1(I,J)) /                 DBLOC(I,J,K-1) = gravity * (RHOK(I,J) - RHOKM1(I,J)) /
1606       &                                    (RHOK(I,J) + rhonil)       &                                    (RHOK(I,J) + rhoConst)
1607                 DBSFC(I,J,K)   = gravity * (RHOK(I,J) - RHO1K (I,J)) /                 DBSFC(I,J,K)   = gravity * (RHOK(I,J) - RHO1K (I,J)) /
1608       &                                    (RHOK(I,J) + rhonil)       &                                    (RHOK(I,J) + rhoConst)
1609              END DO              END DO
1610           END DO           END DO
1611    
1612    #ifdef ALLOW_DIAGNOSTICS
1613             IF ( useDiagnostics ) THEN
1614    c     work1 - density of t(k-1)  & s(k-1) at depth 1
1615    c     work2 - density of t(k  )  & s(k  ) at depth 1
1616    c     work3 - density of t(1)-.8 & s(1  ) at depth 1
1617                CALL FIND_RHO_2D(
1618         I           1, sNx, 1, sNy, 1,
1619         I           theta(1-OLx,1-OLy,k-1,bi,bj),
1620         I           salt (1-OLx,1-OLy,k-1,bi,bj),
1621         O           WORK1,
1622         I           k-1, bi, bj, myThid )
1623                CALL FIND_RHO_2D(
1624         I           1, sNx, 1, sNy, 1,
1625         I           theta(1-OLx,1-OLy,k,bi,bj),
1626         I           salt (1-OLx,1-OLy,k,bi,bj),
1627         O           WORK2,
1628         I           k, bi, bj, myThid )
1629                DO J = 1, sNy
1630                   DO I = 1, sNx
1631                      IF ( k .LE. klowC(I,J,bi,bj) .AND.
1632         &                 WORK1(I,J) .LT. WORK3(I,J) .AND.
1633         &                 WORK2(I,J) .GE. WORK3(I,J) ) THEN
1634                         KPPMLD(I,J) = MIN ( KPPMLD(I,J),
1635         &                    ABS(((WORK3(I,J)-WORK1(I,J))*rC(k)+
1636         &                    (WORK2(I,J)-WORK3(I,J))*rC(k-1))/
1637         &                    (WORK2(I,J)-WORK1(I,J))))
1638                      ENDIF
1639                   END DO
1640                END DO
1641             ENDIF
1642    #endif /* ALLOW_DIAGNOSTICS */
1643    
1644        END DO        END DO
1645    
1646  c     compute arrays for K = Nrp1  c     compute arrays for K = Nrp1
1647        DO J = jbot, jtop        DO J = 1-OLy, sNy+OLy
1648           DO I = ibot, itop           DO I = 1-OLx, sNx+OLx
1649              TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)              TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)
1650              SSBETA(I,J,Nrp1)  = SSBETA(I,J,Nr)              SSBETA(I,J,Nrp1)  = SSBETA(I,J,Nr)
1651              DBLOC(I,J,Nr)     = 0.              DBLOC(I,J,Nr)     = 0.
1652           END DO           END DO
1653        END DO        END DO
1654    
1655    #ifdef ALLOW_DIAGNOSTICS
1656          IF ( useDiagnostics ) THEN
1657             CALL DIAGNOSTICS_FILL(KPPmld,'KPPmld  ',0,1,3,bi,bj,myThid)
1658             CALL DIAGNOSTICS_FILL(DBSFC ,'KPPdbsfc',0,Nr,0,1,1,myThid)
1659             CALL DIAGNOSTICS_FILL(DBLOC ,'KPPdbloc',0,Nr,0,1,1,myThid)
1660          ENDIF
1661    #endif /* ALLOW_DIAGNOSTICS */
1662    
1663  #endif /* ALLOW_KPP */  #endif /* ALLOW_KPP */
1664    
1665        RETURN        RETURN
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.40
  ViewVC Help
Powered by ViewVC 1.1.22