/[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.17 by edhill, Mon Sep 29 19:24:31 2003 UTC revision 1.42 by heimbach, Fri Feb 13 21:58:36 2009 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4    #include "PACKAGES_CONFIG.h"
5  #include "KPP_OPTIONS.h"  #include "KPP_OPTIONS.h"
6    
7  C-- File kpp_routines.F: subroutines needed to implement  C-- File kpp_routines.F: subroutines needed to implement
# Line 11  C--  o BLDEPTH     - Determine oceanic p Line 12  C--  o BLDEPTH     - Determine oceanic p
12  C--  o WSCALE      - Compute turbulent velocity scales.  C--  o WSCALE      - Compute turbulent velocity scales.
13  C--  o RI_IWMIX    - Compute interior viscosity diffusivity coefficients.  C--  o RI_IWMIX    - Compute interior viscosity diffusivity coefficients.
14  C--  o Z121        - Apply 121 vertical smoothing.  C--  o Z121        - Apply 121 vertical smoothing.
15  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.  
16  C--  o BLMIX       - Boundary layer mixing coefficients.  C--  o BLMIX       - Boundary layer mixing coefficients.
17  C--  o ENHANCE     - Enhance diffusivity at boundary layer interface.  C--  o ENHANCE     - Enhance diffusivity at boundary layer interface.
18  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 20  C--  o STATEKPP    - Compute buoyancy-re
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     , boplume,SPDepth
27    #endif /* ALLOW_SALT_PLUME */
28         I     , dbloc, Ritop, coriol
29         I     , diffusKzS, diffusKzT
30       I     , ikppkey       I     , ikppkey
31       O     , diffus       O     , diffus
32       U     , ghat       U     , ghat
33       O     , hbl )       O     , hbl
34         I     , 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   myTime :: Current time in simulation
62  c     myThid          - thread number for this instance of the routine  c   myIter :: Current iteration number in simulation
63  c     kmtj   (imt)    - number of vertical layers on this row  c   myThid :: My Thread Id. number
64  c     shsq   (imt,Nr) - (local velocity shear)^2                      ((m/s)^2)  c     kmtj   (imt)     - number of vertical layers on this row
65  c     dvsq   (imt,Nr) - (velocity shear re sfc)^2                     ((m/s)^2)  c     msk    (imt)     - surface mask (=1 if water, =0 otherwise)
66  c     ustar  (imt)    - surface friction velocity                         (m/s)  c     shsq   (imt,Nr)  - (local velocity shear)^2                     ((m/s)^2)
67  c     bo     (imt)    - surface turbulent buoy. forcing               (m^2/s^3)  c     dvsq   (imt,Nr)  - (velocity shear re sfc)^2                    ((m/s)^2)
68  c     bosol  (imt)    - radiative buoyancy forcing                    (m^2/s^3)  c     ustar  (imt)     - surface friction velocity                        (m/s)
69  c     dbloc  (imt,Nr) - local delta buoyancy across interfaces          (m/s^2)  c     bo     (imt)     - surface turbulent buoy. forcing              (m^2/s^3)
70  c     dblocSm(imt,Nr) - horizontally smoothed dbloc                     (m/s^2)  c     bosol  (imt)     - radiative buoyancy forcing                   (m^2/s^3)
71  c                         stored in ghat to save space  c     boplume(imt)     - haline buoyancy forcing                      (m^2/s^3)
72  c     Ritop  (imt,Nr) - numerator of bulk Richardson Number  c     dbloc  (imt,Nr)  - local delta buoyancy across interfaces         (m/s^2)
73  c                         (zref-z) * delta buoyancy w.r.t. surface    ((m/s)^2)  c     dblocSm(imt,Nr)  - horizontally smoothed dbloc                    (m/s^2)
74  c     coriol (imt)    - Coriolis parameter                                (1/s)  c                          stored in ghat to save space
75    c     Ritop  (imt,Nr)  - numerator of bulk Richardson Number
76    c                          (zref-z) * delta buoyancy w.r.t. surface   ((m/s)^2)
77    c     coriol (imt)     - Coriolis parameter                               (1/s)
78    c     diffusKzS(imt,Nr)- background vertical diffusivity for scalars    (m^2/s)
79    c     diffusKzT(imt,Nr)- background vertical diffusivity for theta      (m^2/s)
80  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,
81  c           e.g., hbl(sNx,sNy) -> hbl(imt),  c           e.g., hbl(sNx,sNy) -> hbl(imt),
82  c           where hbl(i,j) -> hbl((j-1)*sNx+i)  c           where hbl(i,j) -> hbl((j-1)*sNx+i)
83    
84        _RL     mytime        _RL     myTime
85        integer mythid        integer myIter
86        integer kmtj  (imt   )        integer myThid
87        _KPP_RL shsq  (imt,Nr)        integer kmtj (imt   )
88        _KPP_RL dvsq  (imt,Nr)        _RL shsq     (imt,Nr)
89        _KPP_RL ustar (imt   )        _RL dvsq     (imt,Nr)
90        _KPP_RL bo    (imt   )        _RL ustar    (imt   )
91        _KPP_RL bosol (imt   )        _RL bo       (imt   )
92        _KPP_RL dbloc (imt,Nr)        _RL bosol    (imt   )
93        _KPP_RL Ritop (imt,Nr)  #ifdef ALLOW_SALT_PLUME
94        _KPP_RL coriol(imt   )        _RL boplume  (imt   )
95          _RL SPDepth  (imt   )
96    #endif /* ALLOW_SALT_PLUME */
97          _RL dbloc    (imt,Nr)
98          _RL Ritop    (imt,Nr)
99          _RL coriol   (imt   )
100          _RS msk      (imt   )
101          _RL diffusKzS(imt,Nr)
102          _RL diffusKzT(imt,Nr)
103    
104        integer ikppkey        integer ikppkey
105    
# Line 91  c     diffus (imt,3)  - vertical tempera Line 110  c     diffus (imt,3)  - vertical tempera
110  c     ghat   (imt)    - nonlocal transport coefficient                  (s/m^2)  c     ghat   (imt)    - nonlocal transport coefficient                  (s/m^2)
111  c     hbl    (imt)    - mixing layer depth                                  (m)  c     hbl    (imt)    - mixing layer depth                                  (m)
112    
113        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL diffus(imt,0:Nrp1,mdiff)
114        _KPP_RL ghat  (imt,Nr)        _RL ghat  (imt,Nr)
115        _KPP_RL hbl   (imt)        _RL hbl   (imt)
116    
117  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
118    
# Line 107  c     blmc   (imt,Nr,mdiff) - boundary l Line 126  c     blmc   (imt,Nr,mdiff) - boundary l
126  c     sigma  (imt         ) - normalized depth (d / hbl)  c     sigma  (imt         ) - normalized depth (d / hbl)
127  c     Rib    (imt,Nr      ) - bulk Richardson number  c     Rib    (imt,Nr      ) - bulk Richardson number
128    
129        integer kbl   (imt         )        integer kbl(imt         )
130        _KPP_RL bfsfc (imt         )        _RL bfsfc  (imt         )
131        _KPP_RL casea (imt         )        _RL casea  (imt         )
132        _KPP_RL stable(imt         )        _RL stable (imt         )
133        _KPP_RL dkm1  (imt,   mdiff)        _RL dkm1   (imt,   mdiff)
134        _KPP_RL blmc  (imt,Nr,mdiff)        _RL blmc   (imt,Nr,mdiff)
135        _KPP_RL sigma (imt         )        _RL sigma  (imt         )
136        _KPP_RL Rib   (imt,Nr      )        _RL Rib    (imt,Nr      )
137    
138        integer i, k, md        integer i, k, md
139    
# Line 127  c--------------------------------------- Line 146  c---------------------------------------
146    
147  cph(  cph(
148  cph these storings avoid recomp. of Ri_iwmix  cph these storings avoid recomp. of Ri_iwmix
149  CADJ STORE ghat  = comlev1_kpp, key = ikppkey  CADJ STORE ghat  = comlev1_kpp, key=ikppkey, kind=isbyte
150  CADJ STORE dbloc = comlev1_kpp, key = ikppkey  CADJ STORE dbloc = comlev1_kpp, key=ikppkey, kind=isbyte
151  cph)  cph)
152        call Ri_iwmix (        call Ri_iwmix (
153       I       kmtj, shsq, dbloc, ghat       I       kmtj, shsq, dbloc, ghat
154         I     , diffusKzS, diffusKzT
155       I     , ikppkey       I     , ikppkey
156       O     , diffus )       O     , diffus, myThid )
157    
158  cph(  cph(
159  cph these storings avoid recomp. of Ri_iwmix  cph these storings avoid recomp. of Ri_iwmix
160  cph DESPITE TAFs 'not necessary' warning!  cph DESPITE TAFs 'not necessary' warning!
161  CADJ STORE dbloc  = comlev1_kpp, key = ikppkey  CADJ STORE dbloc  = comlev1_kpp, key=ikppkey, kind=isbyte
162  CADJ STORE shsq   = comlev1_kpp, key = ikppkey  CADJ STORE shsq   = comlev1_kpp, key=ikppkey, kind=isbyte
163  CADJ STORE ghat   = comlev1_kpp, key = ikppkey  CADJ STORE ghat   = comlev1_kpp, key=ikppkey, kind=isbyte
164  CADJ STORE diffus = comlev1_kpp, key = ikppkey  CADJ STORE diffus = comlev1_kpp, key=ikppkey, kind=isbyte
165  cph)  cph)
166    
167  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 150  c for blmix Line 170  c for blmix
170  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
171    
172        do md = 1, mdiff        do md = 1, mdiff
173           do k=1,Nrp1
174           do i = 1,imt           do i = 1,imt
175              do k=kmtj(i),Nrp1               if(k.ge.kmtj(i))  diffus(i,k,md) = 0.0
                diffus(i,k,md) = 0.0  
176              end do              end do
177           end do           end do
178        end do        end do
# Line 164  c diagnose the new boundary layer depth Line 184  c diagnose the new boundary layer depth
184  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
185    
186        call bldepth (        call bldepth (
187       I       mytime, mythid       I       kmtj
188       I     , kmtj       I     , dvsq, dbloc, Ritop, ustar, bo, bosol
189       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol  #ifdef ALLOW_SALT_PLUME
190         I     , boplume,SPDepth
191    #endif /* ALLOW_SALT_PLUME */
192         I     , coriol
193       I     , ikppkey       I     , ikppkey
194       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma
195       &     )       I     , myTime, myIter, myThid )
196    
197  CADJ STORE hbl,bfsfc,stable,casea,kbl = comlev1_kpp, key = ikppkey  CADJ STORE hbl,bfsfc,stable,casea,kbl = comlev1_kpp,
198    CADJ &     key=ikppkey, kind=isbyte
199    
200  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
201  c compute boundary layer diffusivities  c compute boundary layer diffusivities
# Line 180  c--------------------------------------- Line 204  c---------------------------------------
204        call blmix (        call blmix (
205       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl
206       O     , dkm1, blmc, ghat, sigma, ikppkey       O     , dkm1, blmc, ghat, sigma, ikppkey
207       &     )       I     , myThid )
208  cph(  cph(
209  CADJ STORE dkm1,blmc,ghat = comlev1_kpp, key = ikppkey  CADJ STORE dkm1,blmc,ghat = comlev1_kpp,
210  CADJ STORE hbl, kbl, diffus, casea = comlev1_kpp, key = ikppkey  CADJ &     key=ikppkey, kind=isbyte
211    CADJ STORE hbl, kbl, diffus, casea = comlev1_kpp,
212    CADJ &     key=ikppkey, kind=isbyte
213  cph)  cph)
214    
215  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 193  c--------------------------------------- Line 219  c---------------------------------------
219        call enhance (        call enhance (
220       I       dkm1, hbl, kbl, diffus, casea       I       dkm1, hbl, kbl, diffus, casea
221       U     , ghat       U     , ghat
222       O     , blmc )       O     , blmc
223         I     , myThid )
224    
225  cph(  cph(
226  cph avoids recomp. of enhance  cph avoids recomp. of enhance
227  CADJ STORE blmc = comlev1_kpp, key = ikppkey  CADJ STORE blmc = comlev1_kpp, key=ikppkey, kind=isbyte
228  cph)  cph)
229    
230  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 209  c--------------------------------------- Line 236  c---------------------------------------
236        do k = 1, Nr        do k = 1, Nr
237           do i = 1, imt           do i = 1, imt
238              if (k .lt. kbl(i)) then              if (k .lt. kbl(i)) then
239    #ifdef ALLOW_SHELFICE
240    C     when there is shelfice on top (msk(i)=0), reset the boundary layer
241    C     mixing coefficients blmc to pure Ri-number based mixing
242                   blmc(i,k,1) = max ( blmc(i,k,1)*msk(i),
243         &              diffus(i,k,1) )
244                   blmc(i,k,2) = max ( blmc(i,k,2)*msk(i),
245         &              diffus(i,k,2) )
246                   blmc(i,k,3) = max ( blmc(i,k,3)*msk(i),
247         &              diffus(i,k,3) )
248    #endif /* not ALLOW_SHELFICE */
249                 diffus(i,k,1) = max ( blmc(i,k,1), viscAr  )                 diffus(i,k,1) = max ( blmc(i,k,1), viscAr  )
250                 diffus(i,k,2) = max ( blmc(i,k,2), diffKrS )                 diffus(i,k,2) = max ( blmc(i,k,2), diffusKzS(i,Nr) )
251                 diffus(i,k,3) = max ( blmc(i,k,3), diffKrT )                 diffus(i,k,3) = max ( blmc(i,k,3), diffusKzT(i,Nr) )
252              else              else
253                 ghat(i,k) = 0.                 ghat(i,k) = 0. _d 0
254              endif              endif
255           end do           end do
256        end do        end do
# Line 226  c--------------------------------------- Line 263  c---------------------------------------
263  c*************************************************************************  c*************************************************************************
264    
265        subroutine bldepth (        subroutine bldepth (
266       I       mytime, mythid       I       kmtj
267       I     , kmtj       I     , dvsq, dbloc, Ritop, ustar, bo, bosol
268       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol  #ifdef ALLOW_SALT_PLUME
269         I     , boplume,SPDepth
270    #endif /* ALLOW_SALT_PLUME */
271         I     , coriol
272       I     , ikppkey       I     , ikppkey
273       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma
274       &     )       I     , myTime, myIter, myThid )
275    
276  c     the oceanic planetary boundary layer depth, hbl, is determined as  c     the oceanic planetary boundary layer depth, hbl, is determined as
277  c     the shallowest depth where the bulk Richardson number is  c     the shallowest depth where the bulk Richardson number is
# Line 260  c Line 300  c
300  #include "EEPARAMS.h"  #include "EEPARAMS.h"
301  #include "PARAMS.h"  #include "PARAMS.h"
302  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
303  #include "FFIELDS.h"  #ifdef ALLOW_AUTODIFF
304    # include "tamc.h"
305    #endif
306    
307  c input  c input
308  c------  c------
309  c myTime    : current time in simulation  c   myTime :: Current time in simulation
310  c myThid    : thread number for this instance of the routine  c   myIter :: Current iteration number in simulation
311    c   myThid :: My Thread Id. number
312  c kmtj      : number of vertical layers  c kmtj      : number of vertical layers
313  c dvsq      : (velocity shear re sfc)^2             ((m/s)^2)  c dvsq      : (velocity shear re sfc)^2             ((m/s)^2)
314  c dbloc     : local delta buoyancy across interfaces  (m/s^2)  c dbloc     : local delta buoyancy across interfaces  (m/s^2)
# Line 275  c             buoyancy with respect to s Line 318  c             buoyancy with respect to s
318  c ustar     : surface friction velocity                 (m/s)  c ustar     : surface friction velocity                 (m/s)
319  c bo        : surface turbulent buoyancy forcing    (m^2/s^3)  c bo        : surface turbulent buoyancy forcing    (m^2/s^3)
320  c bosol     : radiative buoyancy forcing            (m^2/s^3)  c bosol     : radiative buoyancy forcing            (m^2/s^3)
321    c boplume   : haline buoyancy forcing               (m^2/s^3)
322  c coriol    : Coriolis parameter                        (1/s)  c coriol    : Coriolis parameter                        (1/s)
323        _RL     mytime        _RL     myTime
324        integer mythid        integer myIter
325          integer myThid
326        integer kmtj(imt)        integer kmtj(imt)
327        _KPP_RL dvsq  (imt,Nr)        _RL dvsq    (imt,Nr)
328        _KPP_RL dbloc (imt,Nr)        _RL dbloc   (imt,Nr)
329        _KPP_RL Ritop (imt,Nr)        _RL Ritop   (imt,Nr)
330        _KPP_RL ustar (imt)        _RL ustar   (imt)
331        _KPP_RL bo    (imt)        _RL bo      (imt)
332        _KPP_RL bosol (imt)        _RL bosol   (imt)
333        _KPP_RL coriol(imt)        _RL coriol  (imt)
334        integer ikppkey        integer ikppkey, kkppkey
335    #ifdef ALLOW_SALT_PLUME
336          _RL boplume (imt)
337          _RL SPDepth (imt)
338    #endif /* ALLOW_SALT_PLUME */
339    
340  c  output  c  output
341  c--------  c--------
# Line 297  c casea     : =1 in case A, =0 in case B Line 346  c casea     : =1 in case A, =0 in case B
346  c kbl       : -1 of first grid level below hbl  c kbl       : -1 of first grid level below hbl
347  c Rib       : Bulk Richardson number  c Rib       : Bulk Richardson number
348  c sigma     : normalized depth (d/hbl)  c sigma     : normalized depth (d/hbl)
349        _KPP_RL hbl   (imt)        _RL hbl    (imt)
350        _KPP_RL bfsfc (imt)        _RL bfsfc  (imt)
351        _KPP_RL stable(imt)        _RL stable (imt)
352        _KPP_RL casea (imt)        _RL casea  (imt)
353        integer kbl   (imt)        integer kbl(imt)
354        _KPP_RL Rib   (imt,Nr)        _RL Rib    (imt,Nr)
355        _KPP_RL sigma (imt)        _RL sigma  (imt)
356    
357  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
358    
359  c  local  c  local
360  c-------  c-------
361  c wm, ws    : turbulent velocity scales         (m/s)  c wm, ws    : turbulent velocity scales         (m/s)
362        _KPP_RL wm(imt), ws(imt)        _RL wm(imt), ws(imt)
363        _RL     worka(imt)        _RL worka(imt)
364          _RL bvsq, vtsq, hekman, hmonob, hlimit, tempVar1, tempVar2
       _KPP_RL bvsq, vtsq, hekman, hmonob, hlimit, tempVar1, tempVar2  
365        integer i, kl        integer i, kl
366    
367        _KPP_RL     p5    , eins        _RL         p5    , eins
368        parameter ( p5=0.5, eins=1.0 )        parameter ( p5=0.5, eins=1.0 )
369        _RL         minusone        _RL         minusone
370        parameter ( minusone=-1.0 )        parameter ( minusone=-1.0 )
371    
372    #ifdef ALLOW_DIAGNOSTICS
373    c     KPPBFSFC - Bo+radiation absorbed to d=hbf*hbl + plume (m^2/s^3)
374          _RL KPPBFSFC(imt,Nr)
375          _RL KPPRi(imt,Nr)
376    #endif /* ALLOW_DIAGNOSTICS */
377    
378  c find bulk Richardson number at every grid level until > Ricr  c find bulk Richardson number at every grid level until > Ricr
379  c  c
380  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 333  c       kbl(i)=kmtj(i) and hbl(i)=-zgrid Line 387  c       kbl(i)=kmtj(i) and hbl(i)=-zgrid
387  c     initialize hbl and kbl to bottomed out values  c     initialize hbl and kbl to bottomed out values
388    
389        do i = 1, imt        do i = 1, imt
390           Rib(i,1) = 0.0           Rib(i,1) = 0. _d 0
391           kbl(i) = max(kmtj(i),1)           kbl(i) = max(kmtj(i),1)
392           hbl(i) = -zgrid(kbl(i))           hbl(i) = -zgrid(kbl(i))
393        end do        end do
394    
395    #ifdef ALLOW_DIAGNOSTICS
396          do kl = 1, Nr
397             do i = 1, imt
398                KPPBFSFC(i,kl) = 0. _d 0
399                KPPRi(i,kl) = 0. _d 0
400             enddo
401          enddo
402    #endif /* ALLOW_DIAGNOSTICS */
403    
404        do kl = 2, Nr        do kl = 2, Nr
405    
406    #ifdef ALLOW_AUTODIFF_TAMC
407             kkppkey = (ikppkey-1)*Nr + kl
408    #endif
409    
410  c     compute bfsfc = sw fraction at hbf * zgrid  c     compute bfsfc = sw fraction at hbf * zgrid
411    
412           do i = 1, imt           do i = 1, imt
413              worka(i) = zgrid(kl)              worka(i) = zgrid(kl)
414           end do           end do
415    CADJ store worka = comlev1_kpp_k, key = kkppkey, kind=isbyte
416           call SWFRAC(           call SWFRAC(
417       I        imt, hbf,       I       imt, hbf,
418       I        mytime, mythid,       U       worka,
419       U        worka )       I       myTime, myIter, myThid )
420    CADJ store worka = comlev1_kpp_k, key = kkppkey, kind=isbyte
421    
422           do i = 1, imt           do i = 1, imt
423    
# Line 359  c     use caseA as temporary array Line 428  c     use caseA as temporary array
428  c     compute bfsfc= Bo + radiative contribution down to hbf * hbl  c     compute bfsfc= Bo + radiative contribution down to hbf * hbl
429    
430              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  
431    
432           end do           end do
433    #ifdef ALLOW_SALT_PLUME
434    c     compute bfsfc = plume fraction at hbf * zgrid
435             IF ( useSALT_PLUME ) THEN
436               do i = 1, imt
437                  worka(i) = zgrid(kl)
438               enddo
439               call SALT_PLUME_FRAC(
440         I         imt, hbf,SPDepth,
441         U         worka,
442         I         myTime, myIter, myThid)
443               do i = 1, imt
444                  bfsfc(i) = bfsfc(i) + boplume(i)*(1. - worka(i))
445               enddo
446             ENDIF
447    #endif /* ALLOW_SALT_PLUME */
448    
449    #ifdef ALLOW_DIAGNOSTICS
450             do i = 1, imt
451                KPPBFSFC(i,kl) = bfsfc(i)
452             enddo
453    #endif /* ALLOW_DIAGNOSTICS */
454    
455             do i = 1, imt
456                stable(i) = p5 + sign(p5,bfsfc(i))
457                sigma(i) = stable(i) + (1. - stable(i)) * epsilon
458             enddo
459    
460  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
461  c     compute velocity scales at sigma, for hbl= caseA = -zgrid(kl)  c     compute velocity scales at sigma, for hbl= caseA = -zgrid(kl)
# Line 370  c--------------------------------------- Line 463  c---------------------------------------
463    
464           call wscale (           call wscale (
465       I        sigma, casea, ustar, bfsfc,       I        sigma, casea, ustar, bfsfc,
466       O        wm, ws )       O        wm, ws, myThid )
467    CADJ store ws = comlev1_kpp_k, key = kkppkey, kind=isbyte
468    
469           do i = 1, imt           do i = 1, imt
470    
# Line 382  c--------------------------------------- Line 476  c---------------------------------------
476       1           ( dbloc(i,kl-1) / (zgrid(kl-1)-zgrid(kl  ))+       1           ( dbloc(i,kl-1) / (zgrid(kl-1)-zgrid(kl  ))+
477       2             dbloc(i,kl  ) / (zgrid(kl  )-zgrid(kl+1)))       2             dbloc(i,kl  ) / (zgrid(kl  )-zgrid(kl+1)))
478    
479              if (bvsq .eq. 0.) then              if (bvsq .eq. 0. _d 0) then
480                vtsq = 0.0                vtsq = 0. _d 0
481              else              else
482                vtsq = -zgrid(kl) * ws(i) * sqrt(abs(bvsq)) * Vtc                vtsq = -zgrid(kl) * ws(i) * sqrt(abs(bvsq)) * Vtc
483              endif              endif
# Line 403  c Line 497  c
497              tempVar1  = dvsq(i,kl) + vtsq                        tempVar1  = dvsq(i,kl) + vtsq          
498              tempVar2 = max(tempVar1, phepsi)              tempVar2 = max(tempVar1, phepsi)
499              Rib(i,kl) = Ritop(i,kl) / tempVar2              Rib(i,kl) = Ritop(i,kl) / tempVar2
500    #ifdef ALLOW_DIAGNOSTICS
501                KPPRi(i,kl) = Rib(i,kl)
502    #endif
503    
504           end do           end do
505        end do        end do
506    
507    #ifdef ALLOW_DIAGNOSTICS
508             CALL DIAGNOSTICS_FILL(KPPBFSFC,'KPPbfsfc',0,Nr,0,1,1,myThid)
509             CALL DIAGNOSTICS_FILL(KPPRi   ,'KPPRi   ',0,Nr,0,1,1,myThid)
510    #endif /* ALLOW_DIAGNOSTICS */
511    
512  cph(  cph(
513  cph  without this store, there is a recomputation error for  cph  without this store, there is a recomputation error for
514  cph  rib in adbldepth (probably partial recomputation problem)      cph  rib in adbldepth (probably partial recomputation problem)    
515  CADJ store Rib = comlev1_kpp  CADJ store Rib = comlev1_kpp
516  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy),Nr /)  CADJ &   , key=ikppkey, kind=isbyte,
517    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy),Nr /)
518  cph)  cph)
519    
520        do kl = 2, Nr        do kl = 2, Nr
# Line 421  cph) Line 524  cph)
524        end do        end do
525    
526  CADJ store kbl = comlev1_kpp  CADJ store kbl = comlev1_kpp
527  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key=ikppkey, kind=isbyte,
528    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
529    
530        do i = 1, imt        do i = 1, imt
531           kl = kbl(i)           kl = kbl(i)
# Line 434  c     linearly interpolate to find hbl w Line 538  c     linearly interpolate to find hbl w
538        end do        end do
539    
540  CADJ store hbl = comlev1_kpp  CADJ store hbl = comlev1_kpp
541  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key=ikppkey, kind=isbyte,
542    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
543    
544  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
545  c     find stability and buoyancy forcing for boundary layer  c     find stability and buoyancy forcing for boundary layer
# Line 443  c--------------------------------------- Line 548  c---------------------------------------
548        do i = 1, imt        do i = 1, imt
549           worka(i) = hbl(i)           worka(i) = hbl(i)
550        end do        end do
551    CADJ store worka = comlev1_kpp
552    CADJ &   , key=ikppkey, kind=isbyte,
553    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
554        call SWFRAC(        call SWFRAC(
555       I     imt, minusone,       I       imt, minusone,
556       I     mytime, mythid,       U       worka,
557       U     worka )       I       myTime, myIter, myThid )
558    CADJ store worka = comlev1_kpp
559    CADJ &   , key=ikppkey, kind=isbyte,
560    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
561    
562        do i = 1, imt        do i = 1, imt
563           bfsfc(i)  = bo(i) + bosol(i) * (1. - worka(i))           bfsfc(i)  = bo(i) + bosol(i) * (1. - worka(i))
564        end do        end do
565    
566    #ifdef ALLOW_SALT_PLUME
567          IF ( useSALT_PLUME ) THEN
568            do i = 1, imt
569               worka(i) = hbl(i)
570            enddo
571            call SALT_PLUME_FRAC(
572         I         imt,minusone,SPDepth,
573         U         worka,
574         I         myTime, myIter, myThid )
575            do i = 1, imt
576               bfsfc(i) = bfsfc(i) + boplume(i) * (1. - worka(i))
577            enddo
578          ENDIF
579    #endif /* ALLOW_SALT_PLUME */
580  CADJ store bfsfc = comlev1_kpp  CADJ store bfsfc = comlev1_kpp
581  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key=ikppkey, kind=isbyte,
582    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
583    
584  c--   ensure bfsfc is never 0  c--   ensure bfsfc is never 0
585        do i = 1, imt        do i = 1, imt
# Line 463  c--   ensure bfsfc is never 0 Line 590  c--   ensure bfsfc is never 0
590  cph(  cph(
591  cph  added stable to store list to avoid extensive recomp.  cph  added stable to store list to avoid extensive recomp.
592  CADJ store bfsfc, stable = comlev1_kpp  CADJ store bfsfc, stable = comlev1_kpp
593  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key=ikppkey, kind=isbyte,
594    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
595  cph)  cph)
596    
597  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 482  c--------------------------------------- Line 610  c---------------------------------------
610           end if           end if
611        end do        end do
612  CADJ store hbl = comlev1_kpp  CADJ store hbl = comlev1_kpp
613  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key=ikppkey, kind=isbyte,
614    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
615    
616        do i = 1, imt        do i = 1, imt
617           hbl(i) = max(hbl(i),minKPPhbl)           hbl(i) = max(hbl(i),minKPPhbl)
# Line 490  CADJ &   , key = ikppkey, shape = (/ (sN Line 619  CADJ &   , key = ikppkey, shape = (/ (sN
619        end do        end do
620    
621  CADJ store hbl = comlev1_kpp  CADJ store hbl = comlev1_kpp
622  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key=ikppkey, kind=isbyte,
623    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
624    
625  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
626  c      find new kbl  c      find new kbl
# Line 511  c--------------------------------------- Line 641  c---------------------------------------
641        do i = 1, imt        do i = 1, imt
642           worka(i) = hbl(i)           worka(i) = hbl(i)
643        end do        end do
644    CADJ store worka = comlev1_kpp
645    CADJ &   , key=ikppkey, kind=isbyte,
646    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
647        call SWFRAC(        call SWFRAC(
648       I     imt, minusone,       I     imt, minusone,
649       I     mytime, mythid,       U     worka,
650       U     worka )       I     myTime, myIter, myThid )
651    CADJ store worka = comlev1_kpp
652    CADJ &   , key=ikppkey, kind=isbyte,
653    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
654    
655        do i = 1, imt        do i = 1, imt
656           bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i))           bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i))
657        end do        end do
658    
659    #ifdef ALLOW_SALT_PLUME
660          IF ( useSALT_PLUME ) THEN
661            do i = 1, imt
662               worka(i) = hbl(i)
663            enddo
664            call SALT_PLUME_FRAC(
665         I         imt,minusone,SPDepth,
666         U         worka,
667         I         myTime, myIter, myThid )
668            do i = 1, imt
669               bfsfc(i) = bfsfc(i) + boplume(i) * (1. - worka(i))
670            enddo
671          ENDIF
672    #endif /* ALLOW_SALT_PLUME */
673  CADJ store bfsfc = comlev1_kpp  CADJ store bfsfc = comlev1_kpp
674  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key=ikppkey, kind=isbyte,
675    CADJ &     shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
676    
677  c--   ensures bfsfc is never 0  c--   ensures bfsfc is never 0
678        do i = 1, imt        do i = 1, imt
# Line 546  c*************************************** Line 698  c***************************************
698    
699        subroutine wscale (        subroutine wscale (
700       I     sigma, hbl, ustar, bfsfc,       I     sigma, hbl, ustar, bfsfc,
701       O     wm, ws )       O     wm, ws,
702         I     myThid )
703    
704  c     compute turbulent velocity scales.  c     compute turbulent velocity scales.
705  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 567  c sigma   : normalized depth (d/hbl) Line 720  c sigma   : normalized depth (d/hbl)
720  c hbl     : boundary layer depth (m)  c hbl     : boundary layer depth (m)
721  c ustar   : surface friction velocity         (m/s)  c ustar   : surface friction velocity         (m/s)
722  c bfsfc   : total surface buoyancy flux       (m^2/s^3)  c bfsfc   : total surface buoyancy flux       (m^2/s^3)
723        _KPP_RL sigma(imt)  c myThid  : thread number for this instance of the routine
724        _KPP_RL hbl  (imt)        integer myThid
725        _KPP_RL ustar(imt)        _RL sigma(imt)
726        _KPP_RL bfsfc(imt)        _RL hbl  (imt)
727          _RL ustar(imt)
728          _RL bfsfc(imt)
729    
730  c  output  c  output
731  c--------  c--------
732  c wm, ws  : turbulent velocity scales at sigma  c wm, ws  : turbulent velocity scales at sigma
733        _KPP_RL wm(imt), ws(imt)        _RL wm(imt), ws(imt)
734    
735  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
736    
737  c local  c local
738  c------  c------
739  c zehat   : = zeta *  ustar**3  c zehat   : = zeta *  ustar**3
740        _KPP_RL zehat        _RL zehat
741    
742        integer iz, izp1, ju, i, jup1        integer iz, izp1, ju, i, jup1
743        _KPP_RL udiff, zdiff, zfrac, ufrac, fzfrac, wam        _RL udiff, zdiff, zfrac, ufrac, fzfrac, wam
744        _KPP_RL wbm, was, wbs, u3, tempVar        _RL wbm, was, wbs, u3, tempVar
745    
746  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
747  c use lookup table for zehat < zmax only; otherwise use  c use lookup table for zehat < zmax only; otherwise use
# Line 641  c--------------------------------------- Line 796  c---------------------------------------
796  c*************************************************************************  c*************************************************************************
797    
798        subroutine Ri_iwmix (        subroutine Ri_iwmix (
799       I       kmtj, shsq, dbloc, dblocSm       I       kmtj, shsq, dbloc, dblocSm,
800       I     , ikppkey       I       diffusKzS, diffusKzT,
801       O     , diffus )       I       ikppkey,
802         O       diffus,
803         I       myThid )
804    
805  c     compute interior viscosity diffusivity coefficients due  c     compute interior viscosity diffusivity coefficients due
806  c     to shear instability (dependent on a local Richardson number),  c     to shear instability (dependent on a local Richardson number),
# Line 656  c     to static instability (local Richa Line 813  c     to static instability (local Richa
813  #include "EEPARAMS.h"  #include "EEPARAMS.h"
814  #include "PARAMS.h"  #include "PARAMS.h"
815  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
816    #ifdef ALLOW_AUTODIFF
817    # include "tamc.h"
818    #endif
819    
820  c  input  c  input
821  c     kmtj   (imt)         number of vertical layers on this row  c     kmtj   (imt)         number of vertical layers on this row
822  c     shsq   (imt,Nr)      (local velocity shear)^2               ((m/s)^2)  c     shsq   (imt,Nr)      (local velocity shear)^2               ((m/s)^2)
823  c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)  c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)
824  c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)  c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)
825        integer kmtj   (imt)  c     diffusKzS(imt,Nr)- background vertical diffusivity for scalars    (m^2/s)
826        _KPP_RL shsq   (imt,Nr)  c     diffusKzT(imt,Nr)- background vertical diffusivity for theta      (m^2/s)
827        _KPP_RL dbloc  (imt,Nr)  c     myThid :: My Thread Id. number
828        _KPP_RL dblocSm(imt,Nr)        integer kmtj (imt)
829          _RL shsq     (imt,Nr)
830          _RL dbloc    (imt,Nr)
831          _RL dblocSm  (imt,Nr)
832          _RL diffusKzS(imt,Nr)
833          _RL diffusKzT(imt,Nr)
834        integer ikppkey        integer ikppkey
835          integer myThid
836    
837  c output  c output
838  c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)  c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)
839  c     diffus(imt,0:Nrp1,2)  vertical scalar diffusivity             (m^2/s)  c     diffus(imt,0:Nrp1,2)  vertical scalar diffusivity             (m^2/s)
840  c     diffus(imt,0:Nrp1,3)  vertical temperature diffusivity        (m^2/s)  c     diffus(imt,0:Nrp1,3)  vertical temperature diffusivity        (m^2/s)
841        _KPP_RL diffus(imt,0:Nrp1,3)        _RL diffus(imt,0:Nrp1,3)
842    
843  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
844    
845  c local variables  c local variables
846  c     Rig                   local Richardson number  c     Rig                   local Richardson number
847  c     fRi, fcon             function of Rig  c     fRi, fcon             function of Rig
848        _KPP_RL Rig        _RL Rig
849        _KPP_RL fRi, fcon        _RL fRi, fcon
850        _KPP_RL ratio        _RL ratio
851        integer i, ki        integer i, ki
852        _KPP_RL c1, c0        _RL c1, c0
853    
854  #ifdef ALLOW_KPP_VERTICALLY_SMOOTH  #ifdef ALLOW_KPP_VERTICALLY_SMOOTH
855        integer mr        integer mr
# Line 691  CADJ INIT kpp_ri_tape_mr = common, 1 Line 857  CADJ INIT kpp_ri_tape_mr = common, 1
857  #endif  #endif
858    
859  c constants  c constants
860        c1 = 1.0        c1 = 1. _d 0
861        c0 = 0.0        c0 = 0. _d 0
862    
863  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
864  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 702  c     set values at bottom and below to Line 868  c     set values at bottom and below to
868  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
869  C     break data flow dependence on diffus  C     break data flow dependence on diffus
870        diffus(1,1,1) = 0.0        diffus(1,1,1) = 0.0
871    
872          do ki = 1, Nr
873             do i = 1, imt
874                diffus(i,ki,1) = 0.
875                diffus(i,ki,2) = 0.
876                diffus(i,ki,3) = 0.
877             enddo
878          enddo
879  #endif  #endif
880                
881    
882        do ki = 1, Nr        do ki = 1, Nr
883           do i = 1, imt           do i = 1, imt
# Line 719  C     break data flow dependence on diff Line 894  C     break data flow dependence on diff
894              endif              endif
895           end do           end do
896        end do        end do
897    CADJ store diffus = comlev1_kpp, key=ikppkey, kind=isbyte
898    
899  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
900  c     vertically smooth Ri  c     vertically smooth Ri
# Line 729  CADJ store diffus(:,:,1) = kpp_ri_tape_m Line 905  CADJ store diffus(:,:,1) = kpp_ri_tape_m
905  CADJ &  , key=mr, shape=(/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2 /)  CADJ &  , key=mr, shape=(/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2 /)
906    
907           call z121 (           call z121 (
908       U     diffus(1,0,1))       U     diffus(1,0,1),
909         I     myThid )
910        end do        end do
911  #endif  #endif
912    
# Line 756  c  evaluate f of smooth Ri for shear ins Line 933  c  evaluate f of smooth Ri for shear ins
933  c ----------------------------------------------------------------------  c ----------------------------------------------------------------------
934  c            evaluate diffusivities and viscosity  c            evaluate diffusivities and viscosity
935  c    mixing due to internal waves, and shear and static instability  c    mixing due to internal waves, and shear and static instability
936    
937              diffus(i,ki,1) = viscAr + fcon * difmcon + fRi * difm0  #ifndef EXCLUDE_KPP_SHEAR_MIX
938              diffus(i,ki,2) = diffKrS + fcon * difscon + fRi * difs0              if ( .NOT. inAdMode ) then
939              diffus(i,ki,3) = diffKrT + fcon * difscon + fRi * difs0                 diffus(i,ki,1) = viscAr + fcon * difmcon + fRi * difm0
940                   diffus(i,ki,2) = diffusKzS(i,ki)+fcon*difscon+fRi*difs0
941                   diffus(i,ki,3) = diffusKzT(i,ki)+fcon*diftcon+fRi*dift0
942                else
943                   diffus(i,ki,1) = viscAr
944                   diffus(i,ki,2) = diffusKzS(i,ki)
945                   diffus(i,ki,3) = diffusKzT(i,ki)
946                endif
947    #else
948                diffus(i,ki,1) = viscAr
949                diffus(i,ki,2) = diffusKzS(i,ki)
950                diffus(i,ki,3) = diffusKzT(i,ki)
951    #endif
952    
953           end do           end do
954        end do        end do
# Line 781  c         set surface values to 0.0 Line 970  c         set surface values to 0.0
970  c*************************************************************************  c*************************************************************************
971    
972        subroutine z121 (        subroutine z121 (
973       U     v )       U     v,
974         I     myThid )
975    
976  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)
977  c     top (0) value is used as a dummy  c     top (0) value is used as a dummy
# Line 799  c     penetrative convection. Line 989  c     penetrative convection.
989  c input/output  c input/output
990  c-------------  c-------------
991  c v     : 2-D array to be smoothed in Nrp1 direction  c v     : 2-D array to be smoothed in Nrp1 direction
992        _KPP_RL v(imt,0:Nrp1)  c myThid: thread number for this instance of the routine
993          integer myThid
994          _RL v(imt,0:Nrp1)
995    
996  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
997    
998  c local  c local
999        _KPP_RL zwork, zflag        _RL zwork, zflag
1000        _KPP_RL KRi_range(1:Nrp1)        _RL KRi_range(1:Nrp1)
1001        integer i, k, km1, kp1        integer i, k, km1, kp1
1002    
1003        _KPP_RL     p0      , p25       , p5      , p2        _RL         p0      , p25       , p5      , p2
1004        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 )
1005    
1006        KRi_range(Nrp1) = p0        KRi_range(Nrp1) = p0
# Line 862  CADJ STORE v(i,k), zwork = z121tape Line 1054  CADJ STORE v(i,k), zwork = z121tape
1054    
1055  c*************************************************************************  c*************************************************************************
1056    
       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*************************************************************************  
   
1057        subroutine smooth_horiz (        subroutine smooth_horiz (
1058       I     k, bi, bj,       I     k, bi, bj,
1059       U     fld )       U     fld,
1060         I     myThid )
1061    
1062  c     Apply horizontal smoothing to global _RL 2-D array  c     Apply horizontal smoothing to global _RL 2-D array
1063    
1064        IMPLICIT NONE        IMPLICIT NONE
1065  #include "SIZE.h"  #include "SIZE.h"
1066    #include "GRID.h"
1067  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
1068    
1069  c     input  c     input
1070  c     bi, bj : array indices  c     bi, bj : array indices
1071  c     k      : vertical index used for masking  c     k      : vertical index used for masking
1072    c     myThid : thread number for this instance of the routine
1073          INTEGER myThid
1074        integer k, bi, bj        integer k, bi, bj
1075    
1076  c     input/output  c     input/output
# Line 981  c     local Line 1097  c     local
1097              im1 = i-1              im1 = i-1
1098              ip1 = i+1              ip1 = i+1
1099              tempVar =              tempVar =
1100       &           p25   *   pMask(i  ,j  ,k,bi,bj)   +       &           p25   *   maskC(i  ,j  ,k,bi,bj)   +
1101       &           p125  * ( pMask(im1,j  ,k,bi,bj)   +       &           p125  * ( maskC(im1,j  ,k,bi,bj)   +
1102       &                     pMask(ip1,j  ,k,bi,bj)   +       &                     maskC(ip1,j  ,k,bi,bj)   +
1103       &                     pMask(i  ,jm1,k,bi,bj)   +       &                     maskC(i  ,jm1,k,bi,bj)   +
1104       &                     pMask(i  ,jp1,k,bi,bj) ) +       &                     maskC(i  ,jp1,k,bi,bj) ) +
1105       &           p0625 * ( pMask(im1,jm1,k,bi,bj)   +       &           p0625 * ( maskC(im1,jm1,k,bi,bj)   +
1106       &                     pMask(im1,jp1,k,bi,bj)   +       &                     maskC(im1,jp1,k,bi,bj)   +
1107       &                     pMask(ip1,jm1,k,bi,bj)   +       &                     maskC(ip1,jm1,k,bi,bj)   +
1108       &                     pMask(ip1,jp1,k,bi,bj) )       &                     maskC(ip1,jp1,k,bi,bj) )
1109              IF ( tempVar .GE. p25 ) THEN              IF ( tempVar .GE. p25 ) THEN
1110                 fld_tmp(i,j) = (                 fld_tmp(i,j) = (
1111       &              p25  * fld(i  ,j  )*pMask(i  ,j  ,k,bi,bj) +       &              p25  * fld(i  ,j  )*maskC(i  ,j  ,k,bi,bj) +
1112       &              p125 *(fld(im1,j  )*pMask(im1,j  ,k,bi,bj) +       &              p125 *(fld(im1,j  )*maskC(im1,j  ,k,bi,bj) +
1113       &                     fld(ip1,j  )*pMask(ip1,j  ,k,bi,bj) +       &                     fld(ip1,j  )*maskC(ip1,j  ,k,bi,bj) +
1114       &                     fld(i  ,jm1)*pMask(i  ,jm1,k,bi,bj) +       &                     fld(i  ,jm1)*maskC(i  ,jm1,k,bi,bj) +
1115       &                     fld(i  ,jp1)*pMask(i  ,jp1,k,bi,bj))+       &                     fld(i  ,jp1)*maskC(i  ,jp1,k,bi,bj))+
1116       &              p0625*(fld(im1,jm1)*pMask(im1,jm1,k,bi,bj) +       &              p0625*(fld(im1,jm1)*maskC(im1,jm1,k,bi,bj) +
1117       &                     fld(im1,jp1)*pMask(im1,jp1,k,bi,bj) +       &                     fld(im1,jp1)*maskC(im1,jp1,k,bi,bj) +
1118       &                     fld(ip1,jm1)*pMask(ip1,jm1,k,bi,bj) +       &                     fld(ip1,jm1)*maskC(ip1,jm1,k,bi,bj) +
1119       &                     fld(ip1,jp1)*pMask(ip1,jp1,k,bi,bj)))       &                     fld(ip1,jp1)*maskC(ip1,jp1,k,bi,bj)))
1120       &              / tempVar       &              / tempVar
1121              ELSE              ELSE
1122                 fld_tmp(i,j) = fld(i,j)                 fld_tmp(i,j) = fld(i,j)
# Line 1025  c*************************************** Line 1141  c***************************************
1141        subroutine blmix (        subroutine blmix (
1142       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl
1143       O     , dkm1, blmc, ghat, sigma, ikppkey       O     , dkm1, blmc, ghat, sigma, ikppkey
1144       &     )       I     , myThid )
1145    
1146  c     mixing coefficients within boundary layer depend on surface  c     mixing coefficients within boundary layer depend on surface
1147  c     forcing and the magnitude and gradient of interior mixing below  c     forcing and the magnitude and gradient of interior mixing below
# Line 1039  c Line 1155  c
1155    
1156  #include "SIZE.h"  #include "SIZE.h"
1157  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
1158    #ifdef ALLOW_AUTODIFF
1159    # include "tamc.h"
1160    #endif
1161    
1162  c input  c input
1163  c     ustar (imt)                 surface friction velocity             (m/s)  c     ustar (imt)                 surface friction velocity             (m/s)
# Line 1047  c     hbl   (imt)                 bounda Line 1166  c     hbl   (imt)                 bounda
1166  c     stable(imt)                 = 1 in stable forcing  c     stable(imt)                 = 1 in stable forcing
1167  c     casea (imt)                 = 1 in case A  c     casea (imt)                 = 1 in case A
1168  c     diffus(imt,0:Nrp1,mdiff)    vertical diffusivities              (m^2/s)  c     diffus(imt,0:Nrp1,mdiff)    vertical diffusivities              (m^2/s)
1169  c     kbl(imt)                    -1 of first grid level below hbl  c     kbl   (imt)                 -1 of first grid level below hbl
1170        _KPP_RL ustar (imt)  c     myThid               thread number for this instance of the routine
1171        _KPP_RL bfsfc (imt)        integer myThid
1172        _KPP_RL hbl   (imt)        _RL ustar (imt)
1173        _KPP_RL stable(imt)        _RL bfsfc (imt)
1174        _KPP_RL casea (imt)        _RL hbl   (imt)
1175        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL stable(imt)
1176          _RL casea (imt)
1177          _RL diffus(imt,0:Nrp1,mdiff)
1178        integer kbl(imt)        integer kbl(imt)
1179    
1180  c output  c output
# Line 1061  c     dkm1 (imt,mdiff)            bounda Line 1182  c     dkm1 (imt,mdiff)            bounda
1182  c     blmc (imt,Nr,mdiff)         boundary layer mixing coefficients  (m^2/s)  c     blmc (imt,Nr,mdiff)         boundary layer mixing coefficients  (m^2/s)
1183  c     ghat (imt,Nr)               nonlocal scalar transport  c     ghat (imt,Nr)               nonlocal scalar transport
1184  c     sigma(imt)                  normalized depth (d / hbl)  c     sigma(imt)                  normalized depth (d / hbl)
1185        _KPP_RL dkm1 (imt,mdiff)        _RL dkm1 (imt,mdiff)
1186        _KPP_RL blmc (imt,Nr,mdiff)        _RL blmc (imt,Nr,mdiff)
1187        _KPP_RL ghat (imt,Nr)        _RL ghat (imt,Nr)
1188        _KPP_RL sigma(imt)        _RL sigma(imt)
1189        integer ikppkey        integer ikppkey, kkppkey
1190    
1191  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1192    
# Line 1073  c  local Line 1194  c  local
1194  c     gat1*(imt)                 shape function at sigma = 1  c     gat1*(imt)                 shape function at sigma = 1
1195  c     dat1*(imt)                 derivative of shape function at sigma = 1  c     dat1*(imt)                 derivative of shape function at sigma = 1
1196  c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)  c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)
1197        _KPP_RL gat1m(imt), gat1s(imt), gat1t(imt)        _RL gat1m(imt), gat1s(imt), gat1t(imt)
1198        _KPP_RL dat1m(imt), dat1s(imt), dat1t(imt)        _RL dat1m(imt), dat1s(imt), dat1t(imt)
1199        _KPP_RL ws(imt), wm(imt)        _RL ws(imt), wm(imt)
1200        integer i, kn, ki        integer i, kn, ki
1201        _KPP_RL R, dvdzup, dvdzdn, viscp        _RL R, dvdzup, dvdzdn, viscp
1202        _KPP_RL difsp, diftp, visch, difsh, difth        _RL difsp, diftp, visch, difsh, difth
1203        _KPP_RL f1, sig, a1, a2, a3, delhat        _RL f1, sig, a1, a2, a3, delhat
1204        _KPP_RL Gm, Gs, Gt        _RL Gm, Gs, Gt
1205        _KPP_RL tempVar        _RL tempVar
1206    
1207        _KPP_RL    p0    , eins        _RL    p0    , eins
1208        parameter (p0=0.0, eins=1.0)        parameter (p0=0.0, eins=1.0)
1209    
1210  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 1094  c--------------------------------------- Line 1215  c---------------------------------------
1215           sigma(i) = stable(i) * 1.0 + (1. - stable(i)) * epsilon           sigma(i) = stable(i) * 1.0 + (1. - stable(i)) * epsilon
1216        end do        end do
1217    
1218    CADJ STORE sigma = comlev1_kpp, key=ikppkey, kind=isbyte
1219        call wscale (        call wscale (
1220       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1221       O        wm, ws )       O        wm, ws, myThid )
1222    CADJ STORE wm = comlev1_kpp, key=ikppkey, kind=isbyte
1223    CADJ STORE ws = comlev1_kpp, key=ikppkey, kind=isbyte
1224    
1225        do i = 1, imt        do i = 1, imt
1226           wm(i) = sign(eins,wm(i))*max(phepsi,abs(wm(i)))           wm(i) = sign(eins,wm(i))*max(phepsi,abs(wm(i)))
1227           ws(i) = sign(eins,ws(i))*max(phepsi,abs(ws(i)))           ws(i) = sign(eins,ws(i))*max(phepsi,abs(ws(i)))
1228        end do        end do
1229  CADJ STORE wm = comlev1_kpp, key = ikppkey  CADJ STORE wm = comlev1_kpp, key=ikppkey, kind=isbyte
1230  CADJ STORE ws = comlev1_kpp, key = ikppkey  CADJ STORE ws = comlev1_kpp, key=ikppkey, kind=isbyte
1231    
1232        do i = 1, imt        do i = 1, imt
1233    
# Line 1139  c--------------------------------------- Line 1263  c---------------------------------------
1263       &        max(ustar(i)**4,phepsi)       &        max(ustar(i)**4,phepsi)
1264           gat1m(i) = visch / hbl(i) / wm(i)           gat1m(i) = visch / hbl(i) / wm(i)
1265           dat1m(i) = -viscp / wm(i) + f1 * visch           dat1m(i) = -viscp / wm(i) + f1 * visch
          dat1m(i) = min(dat1m(i),p0)  
1266    
1267           gat1s(i) = difsh  / hbl(i) / ws(i)           gat1s(i) = difsh  / hbl(i) / ws(i)
1268           dat1s(i) = -difsp / ws(i) + f1 * difsh           dat1s(i) = -difsp / ws(i) + f1 * difsh
          dat1s(i) = min(dat1s(i),p0)  
1269    
1270           gat1t(i) = difth /  hbl(i) / ws(i)           gat1t(i) = difth /  hbl(i) / ws(i)
1271           dat1t(i) = -diftp / ws(i) + f1 * difth           dat1t(i) = -diftp / ws(i) + f1 * difth
          dat1t(i) = min(dat1t(i),p0)  
1272    
1273        end do        end do
1274    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1275    CADJ STORE gat1m = comlev1_kpp, key=ikppkey, kind=isbyte
1276    CADJ STORE gat1s = comlev1_kpp, key=ikppkey, kind=isbyte
1277    CADJ STORE gat1t = comlev1_kpp, key=ikppkey, kind=isbyte
1278    CADJ STORE dat1m = comlev1_kpp, key=ikppkey, kind=isbyte
1279    CADJ STORE dat1s = comlev1_kpp, key=ikppkey, kind=isbyte
1280    CADJ STORE dat1t = comlev1_kpp, key=ikppkey, kind=isbyte
1281    #endif
1282          do i = 1, imt
1283             dat1m(i) = min(dat1m(i),p0)
1284             dat1s(i) = min(dat1s(i),p0)
1285             dat1t(i) = min(dat1t(i),p0)
1286          end do
1287    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1288    CADJ STORE dat1m = comlev1_kpp, key=ikppkey, kind=isbyte
1289    CADJ STORE dat1s = comlev1_kpp, key=ikppkey, kind=isbyte
1290    CADJ STORE dat1t = comlev1_kpp, key=ikppkey, kind=isbyte
1291    #endif
1292    
1293        do ki = 1, Nr        do ki = 1, Nr
1294    
1295    #ifdef ALLOW_AUTODIFF_TAMC
1296             kkppkey = (ikppkey-1)*Nr + ki
1297    #endif
1298    
1299  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1300  c     compute turbulent velocity scales on the interfaces  c     compute turbulent velocity scales on the interfaces
1301  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 1161  c--------------------------------------- Line 1304  c---------------------------------------
1304              sig      = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i)              sig      = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i)
1305              sigma(i) = stable(i)*sig + (1.-stable(i))*min(sig,epsilon)              sigma(i) = stable(i)*sig + (1.-stable(i))*min(sig,epsilon)
1306           end do           end do
1307    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1308    CADJ STORE wm = comlev1_kpp_k, key = kkppkey
1309    CADJ STORE ws = comlev1_kpp_k, key = kkppkey
1310    #endif
1311    CADJ STORE sigma = comlev1_kpp_k, key = kkppkey
1312           call wscale (           call wscale (
1313       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1314       O        wm, ws )       O        wm, ws, myThid )
1315    CADJ STORE wm = comlev1_kpp_k, key = kkppkey
1316    CADJ STORE ws = comlev1_kpp_k, key = kkppkey
1317    
1318  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1319  c     compute the dimensionless shape functions at the interfaces  c     compute the dimensionless shape functions at the interfaces
# Line 1207  c--------------------------------------- Line 1357  c---------------------------------------
1357       &            + (1. - stable(i)) * min(sig,epsilon)       &            + (1. - stable(i)) * min(sig,epsilon)
1358        end do        end do
1359    
1360    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1361    CADJ STORE wm = comlev1_kpp, key=ikppkey, kind=isbyte
1362    CADJ STORE ws = comlev1_kpp, key=ikppkey, kind=isbyte
1363    #endif
1364    CADJ STORE sigma = comlev1_kpp, key=ikppkey, kind=isbyte
1365        call wscale (        call wscale (
1366       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1367       O        wm, ws )       O        wm, ws, myThid )
1368    CADJ STORE wm = comlev1_kpp, key=ikppkey, kind=isbyte
1369    CADJ STORE ws = comlev1_kpp, key=ikppkey, kind=isbyte
1370    
1371        do i = 1, imt        do i = 1, imt
1372           sig = -zgrid(kbl(i)-1) / hbl(i)           sig = -zgrid(kbl(i)-1) / hbl(i)
# Line 1235  c*************************************** Line 1392  c***************************************
1392       I       dkm1, hbl, kbl, diffus, casea       I       dkm1, hbl, kbl, diffus, casea
1393       U     , ghat       U     , ghat
1394       O     , blmc       O     , blmc
1395       &     )       &     , myThid )
1396    
1397  c enhance the diffusivity at the kbl-.5 interface  c enhance the diffusivity at the kbl-.5 interface
1398    
# Line 1250  c     hbl(imt)                  boundary Line 1407  c     hbl(imt)                  boundary
1407  c     kbl(imt)                  grid above hbl  c     kbl(imt)                  grid above hbl
1408  c     diffus(imt,0:Nrp1,mdiff) vertical diffusivities           (m^2/s)  c     diffus(imt,0:Nrp1,mdiff) vertical diffusivities           (m^2/s)
1409  c     casea(imt)                = 1 in caseA, = 0 in case B  c     casea(imt)                = 1 in caseA, = 0 in case B
1410        _KPP_RL dkm1  (imt,mdiff)  c     myThid                    thread number for this instance of the routine
1411        _KPP_RL hbl   (imt)        integer   myThid
1412          _RL dkm1  (imt,mdiff)
1413          _RL hbl   (imt)
1414        integer kbl   (imt)        integer kbl   (imt)
1415        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL diffus(imt,0:Nrp1,mdiff)
1416        _KPP_RL casea (imt)        _RL casea (imt)
1417    
1418  c input/output  c input/output
1419  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)
1420        _KPP_RL ghat (imt,Nr)        _RL ghat (imt,Nr)
1421    
1422  c output  c output
1423  c     enhanced bound. layer mixing coeff.  c     enhanced bound. layer mixing coeff.
1424        _KPP_RL blmc  (imt,Nr,mdiff)        _RL blmc  (imt,Nr,mdiff)
1425    
1426  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1427    
1428  c local  c local
1429  c     fraction hbl lies beteen zgrid neighbors  c     fraction hbl lies beteen zgrid neighbors
1430        _KPP_RL delta        _RL delta
1431        integer ki, i, md        integer ki, i, md
1432        _KPP_RL dkmp5, dstar        _RL dkmp5, dstar
1433    
1434        do i = 1, imt        do i = 1, imt
1435           ki = kbl(i)-1           ki = kbl(i)-1
# Line 1296  c     fraction hbl lies beteen zgrid nei Line 1455  c     fraction hbl lies beteen zgrid nei
1455  c*************************************************************************  c*************************************************************************
1456    
1457        SUBROUTINE STATEKPP (        SUBROUTINE STATEKPP (
1458       I     bi, bj, myThid,       O     RHO1, DBLOC, DBSFC, TTALPHA, SSBETA,
1459       O     RHO1, DBLOC, DBSFC, TTALPHA, SSBETA)       I     ikppkey, bi, bj, myThid )
1460  c  c
1461  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1462  c     "statekpp" computes all necessary input arrays  c     "statekpp" computes all necessary input arrays
# Line 1322  c Line 1481  c
1481  c     written  by: jan morzel,   feb. 10, 1995 (converted from "sigma" version)  c     written  by: jan morzel,   feb. 10, 1995 (converted from "sigma" version)
1482  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
1483  c  c
1484    
1485  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1486    
1487        IMPLICIT NONE        IMPLICIT NONE
# Line 1331  c--------------------------------------- Line 1491  c---------------------------------------
1491  #include "PARAMS.h"  #include "PARAMS.h"
1492  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
1493  #include "DYNVARS.h"  #include "DYNVARS.h"
1494    #include "GRID.h"
1495    #ifdef ALLOW_AUTODIFF
1496    # include "tamc.h"
1497    #endif
1498    
1499  c-------------- Routine arguments -----------------------------------------  c-------------- Routine arguments -----------------------------------------
1500        INTEGER bi, bj, myThid        INTEGER bi, bj, myThid
1501        _KPP_RL RHO1   ( ibot:itop, jbot:jtop       )        _RL RHO1   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy       )
1502        _KPP_RL DBLOC  ( ibot:itop, jbot:jtop, Nr   )        _RL DBLOC  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
1503        _KPP_RL DBSFC  ( ibot:itop, jbot:jtop, Nr   )        _RL DBSFC  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
1504        _KPP_RL TTALPHA( ibot:itop, jbot:jtop, Nrp1 )        _RL TTALPHA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1505        _KPP_RL SSBETA ( ibot:itop, jbot:jtop, Nrp1 )        _RL SSBETA ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1506    
1507  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1508    
# Line 1349  c Line 1513  c
1513  c     rhok         - density of t(k  ) & s(k  ) at depth k  c     rhok         - density of t(k  ) & s(k  ) at depth k
1514  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
1515  c     rho1k        - density of t(1  ) & s(1  ) at depth k  c     rho1k        - density of t(1  ) & s(1  ) at depth k
1516  c     work1, work2 - work arrays for holding horizontal slabs  c     work1,2,3    - work arrays for holding horizontal slabs
1517    
1518        _RL RHOK  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL RHOK  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1519        _RL RHOKM1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL RHOKM1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 1357  c     work1, work2 - work arrays for hol Line 1521  c     work1, work2 - work arrays for hol
1521        _RL WORK1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1522        _RL WORK2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1523        _RL WORK3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1524    
1525        INTEGER I, J, K        INTEGER I, J, K
1526          INTEGER ikppkey, kkppkey
1527    
1528  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
1529    
1530        call FIND_RHO(        kkppkey = (ikppkey-1)*Nr + 1
1531       I     bi, bj, ibot, itop, jbot, jtop, 1, 1,  
1532       I     theta, salt,  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1533    CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k,
1534    CADJ &     key=kkppkey, kind=isbyte
1535    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k,
1536    CADJ &     key=kkppkey, kind=isbyte
1537    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1538          CALL FIND_RHO_2D(
1539         I     1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1,
1540         I     theta(1-OLx,1-OLy,1,bi,bj), salt(1-OLx,1-OLy,1,bi,bj),
1541       O     WORK1,       O     WORK1,
1542       I     myThid )       I     1, bi, bj, myThid )
1543    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1544    CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k,
1545    CADJ &     key=kkppkey, kind=isbyte
1546    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k,
1547    CADJ &     key=kkppkey, kind=isbyte
1548    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1549    
1550        call FIND_ALPHA(        call FIND_ALPHA(
1551       I     bi, bj, ibot, itop, jbot, jtop, 1, 1,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1552       O     WORK2 )       O     WORK2, myThid )
1553    
1554        call FIND_BETA(        call FIND_BETA(
1555       I     bi, bj, ibot, itop, jbot, jtop, 1, 1,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1556       O     WORK3 )       O     WORK3, myThid )
1557    
1558        DO J = jbot, jtop        DO J = 1-OLy, sNy+OLy
1559           DO I = ibot, itop           DO I = 1-OLx, sNx+OLx
1560              RHO1(I,J)      = WORK1(I,J) + rhoConst              RHO1(I,J)      = WORK1(I,J) + rhoConst
1561              TTALPHA(I,J,1) = WORK2(I,J)              TTALPHA(I,J,1) = WORK2(I,J)
1562              SSBETA(I,J,1)  = WORK3(I,J)              SSBETA(I,J,1)  = WORK3(I,J)
# Line 1389  c calculate alpha, beta, and gradients i Line 1569  c calculate alpha, beta, and gradients i
1569  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)
1570        DO K = 2, Nr        DO K = 2, Nr
1571    
1572           call FIND_RHO(        kkppkey = (ikppkey-1)*Nr + k
1573       I        bi, bj, ibot, itop, jbot, jtop, K, K,  
1574       I        theta, salt,  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1575    CADJ STORE theta(:,:,k,bi,bj) = comlev1_kpp_k,
1576    CADJ &     key=kkppkey, kind=isbyte
1577    CADJ STORE salt (:,:,k,bi,bj) = comlev1_kpp_k,
1578    CADJ &     key=kkppkey, kind=isbyte
1579    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1580             CALL FIND_RHO_2D(
1581         I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, k,
1582         I        theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
1583       O        RHOK,       O        RHOK,
1584       I        myThid )       I        k, bi, bj, myThid )
1585    
1586           call FIND_RHO(  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1587       I        bi, bj, ibot, itop, jbot, jtop, K-1, K,  CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_kpp_k,
1588       I        theta, salt,  CADJ &     key=kkppkey, kind=isbyte
1589    CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_kpp_k,
1590    CADJ &     key=kkppkey, kind=isbyte
1591    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1592             CALL FIND_RHO_2D(
1593         I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, k,
1594         I        theta(1-OLx,1-OLy,k-1,bi,bj),salt(1-OLx,1-OLy,k-1,bi,bj),
1595       O        RHOKM1,       O        RHOKM1,
1596       I        myThid )       I        k-1, bi, bj, myThid )
1597    
1598           call FIND_RHO(  #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1599       I        bi, bj, ibot, itop, jbot, jtop, 1, K,  CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k,
1600       I        theta, salt,  CADJ &     key=kkppkey, kind=isbyte
1601    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k,
1602    CADJ &     key=kkppkey, kind=isbyte
1603    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1604             CALL FIND_RHO_2D(
1605         I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, k,
1606         I        theta(1-OLx,1-OLy,1,bi,bj), salt(1-OLx,1-OLy,1,bi,bj),
1607       O        RHO1K,       O        RHO1K,
1608       I        myThid )       I        1, bi, bj, myThid )
1609    
1610    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1611    CADJ STORE rhok  (:,:)          = comlev1_kpp_k,
1612    CADJ &     key=kkppkey, kind=isbyte
1613    CADJ STORE rhokm1(:,:)          = comlev1_kpp_k,
1614    CADJ &     key=kkppkey, kind=isbyte
1615    CADJ STORE rho1k (:,:)          = comlev1_kpp_k,
1616    CADJ &     key=kkppkey, kind=isbyte
1617    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1618    
1619           call FIND_ALPHA(           call FIND_ALPHA(
1620       I        bi, bj, ibot, itop, jbot, jtop, K, K,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1621       O        WORK1 )       O        WORK1, myThid )
1622    
1623           call FIND_BETA(           call FIND_BETA(
1624       I        bi, bj, ibot, itop, jbot, jtop, K, K,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1625       O        WORK2 )       O        WORK2, myThid )
1626    
1627           DO J = jbot, jtop           DO J = 1-OLy, sNy+OLy
1628              DO I = ibot, itop              DO I = 1-OLx, sNx+OLx
1629                 TTALPHA(I,J,K) = WORK1 (I,J)                 TTALPHA(I,J,K) = WORK1 (I,J)
1630                 SSBETA(I,J,K)  = WORK2 (I,J)                 SSBETA(I,J,K)  = WORK2 (I,J)
1631                 DBLOC(I,J,K-1) = gravity * (RHOK(I,J) - RHOKM1(I,J)) /                 DBLOC(I,J,K-1) = gravity * (RHOK(I,J) - RHOKM1(I,J)) /
# Line 1429  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO Line 1638  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO
1638        END DO        END DO
1639    
1640  c     compute arrays for K = Nrp1  c     compute arrays for K = Nrp1
1641        DO J = jbot, jtop        DO J = 1-OLy, sNy+OLy
1642           DO I = ibot, itop           DO I = 1-OLx, sNx+OLx
1643              TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)              TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)
1644              SSBETA(I,J,Nrp1)  = SSBETA(I,J,Nr)              SSBETA(I,J,Nrp1)  = SSBETA(I,J,Nr)
1645              DBLOC(I,J,Nr)     = 0.              DBLOC(I,J,Nr)     = 0.
1646           END DO           END DO
1647        END DO        END DO
1648    
1649    #ifdef ALLOW_DIAGNOSTICS
1650          IF ( useDiagnostics ) THEN
1651             CALL DIAGNOSTICS_FILL(DBSFC ,'KPPdbsfc',0,Nr,0,1,1,myThid)
1652             CALL DIAGNOSTICS_FILL(DBLOC ,'KPPdbloc',0,Nr,0,1,1,myThid)
1653          ENDIF
1654    #endif /* ALLOW_DIAGNOSTICS */
1655    
1656  #endif /* ALLOW_KPP */  #endif /* ALLOW_KPP */
1657    
1658        RETURN        RETURN
Legend:
Removed from v.1.17  
changed lines
  Added in v.1.42
  ViewVC Help
Powered by ViewVC 1.1.22