/[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.3 by heimbach, Wed Sep 13 17:07:11 2000 UTC revision 1.21 by jmc, Sun Oct 17 23:05:09 2004 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "KPP_OPTIONS.h"  #include "KPP_OPTIONS.h"
5    
# Line 10  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 SMOOTH_HORIZ_RS - Apply horizontal smoothing to RS array.  C--  o KPP_SMOOTH_HORIZ - Apply horizontal smoothing to KPP array.
15  C--  o SMOOTH_HORIZ_RL - Apply horizontal smoothing to RL 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.
19    
20  c*************************************************************************  c***********************************************************************
21    
22        SUBROUTINE KPPMIX (        SUBROUTINE KPPMIX (
23       I     lri, kmtj, shsq, dvsq, ustar, bo, bosol       I       mytime, mythid
24       I     , dbloc, Ritop, coriol       I     , kmtj, shsq, dvsq, ustar
25       I     , ikey       I     , bo, bosol, dbloc, Ritop, coriol
26         I     , ikppkey
27       O     , diffus       O     , diffus
28       U     , ghat       U     , ghat
29       O     , hbl )       O     , hbl )
30    
31  c-------------------------------------------------------------------------  c-----------------------------------------------------------------------
32  c  c
33  c     Main driver subroutine for kpp vertical mixing scheme and  c     Main driver subroutine for kpp vertical mixing scheme and
34  c     interface to greater ocean model  c     interface to greater ocean model
# Line 50  c--------------------------------------- Line 52  c---------------------------------------
52  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
53    
54  c input  c input
55  c     lri             - mixing process switches  c     myTime          - current time in simulation
56    c     myThid          - thread number for this instance of the routine
57  c     kmtj   (imt)    - number of vertical layers on this row  c     kmtj   (imt)    - number of vertical layers on this row
58  c     shsq   (imt,Nr) - (local velocity shear)^2                      ((m/s)^2)  c     shsq   (imt,Nr) - (local velocity shear)^2                      ((m/s)^2)
59  c     dvsq   (imt,Nr) - (velocity shear re sfc)^2                     ((m/s)^2)  c     dvsq   (imt,Nr) - (velocity shear re sfc)^2                     ((m/s)^2)
# Line 67  c     note: there is a conversion from 2 Line 70  c     note: there is a conversion from 2
70  c           e.g., hbl(sNx,sNy) -> hbl(imt),  c           e.g., hbl(sNx,sNy) -> hbl(imt),
71  c           where hbl(i,j) -> hbl((j-1)*sNx+i)  c           where hbl(i,j) -> hbl((j-1)*sNx+i)
72    
73        logical lri        _RL     mytime
74          integer mythid
75        integer kmtj  (imt   )        integer kmtj  (imt   )
76        _RS    shsq   (imt,Nr)        _KPP_RL shsq  (imt,Nr)
77        _RS    dvsq   (imt,Nr)        _KPP_RL dvsq  (imt,Nr)
78        _RS    ustar  (imt   )        _KPP_RL ustar (imt   )
79        _RS    bo     (imt   )        _KPP_RL bo    (imt   )
80        _RS    bosol  (imt   )        _KPP_RL bosol (imt   )
81        _RS    dbloc  (imt,Nr)        _KPP_RL dbloc (imt,Nr)
82        _RS    Ritop  (imt,Nr)        _KPP_RL Ritop (imt,Nr)
83        _RS    coriol (imt   )        _KPP_RL coriol(imt   )
84    
85        integer ikey        integer ikppkey
86    
87  c output  c output
88  c     diffus (imt,1)  - vertical viscosity coefficient                  (m^2/s)  c     diffus (imt,1)  - vertical viscosity coefficient                  (m^2/s)
# Line 87  c     diffus (imt,3)  - vertical tempera Line 91  c     diffus (imt,3)  - vertical tempera
91  c     ghat   (imt)    - nonlocal transport coefficient                  (s/m^2)  c     ghat   (imt)    - nonlocal transport coefficient                  (s/m^2)
92  c     hbl    (imt)    - mixing layer depth                                  (m)  c     hbl    (imt)    - mixing layer depth                                  (m)
93    
94        _RS diffus(imt,0:Nrp1,mdiff)        _KPP_RL diffus(imt,0:Nrp1,mdiff)
95        _RS ghat  (imt,Nr)        _KPP_RL ghat  (imt,Nr)
96        _RS hbl   (imt)        _KPP_RL hbl   (imt)
97    
98  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
99    
# Line 103  c     blmc   (imt,Nr,mdiff) - boundary l Line 107  c     blmc   (imt,Nr,mdiff) - boundary l
107  c     sigma  (imt         ) - normalized depth (d / hbl)  c     sigma  (imt         ) - normalized depth (d / hbl)
108  c     Rib    (imt,Nr      ) - bulk Richardson number  c     Rib    (imt,Nr      ) - bulk Richardson number
109    
110        integer kbl (imt         )        integer kbl   (imt         )
111        _RS bfsfc   (imt         )        _KPP_RL bfsfc (imt         )
112        _RS casea   (imt         )        _KPP_RL casea (imt         )
113        _RS stable  (imt         )        _KPP_RL stable(imt         )
114        _RS dkm1    (imt,   mdiff)        _KPP_RL dkm1  (imt,   mdiff)
115        _RS blmc    (imt,Nr,mdiff)        _KPP_RL blmc  (imt,Nr,mdiff)
116        _RS sigma   (imt         )        _KPP_RL sigma (imt         )
117        _RS Rib     (imt,Nr      )        _KPP_RL Rib   (imt,Nr      )
118    
119        integer i, k, md        integer i, k, md
120    
# Line 121  c instability. Line 125  c instability.
125  c (ghat is temporary storage for horizontally smoothed dbloc)  c (ghat is temporary storage for horizontally smoothed dbloc)
126  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
127    
128    cph(
129    cph these storings avoid recomp. of Ri_iwmix
130    CADJ STORE ghat  = comlev1_kpp, key = ikppkey
131    CADJ STORE dbloc = comlev1_kpp, key = ikppkey
132    cph)
133        call Ri_iwmix (        call Ri_iwmix (
134       I       kmtj, shsq, dbloc, ghat       I       kmtj, shsq, dbloc, ghat
135       I     , ikey       I     , ikppkey
136       O     , diffus )       O     , diffus )
137    
138    cph(
139    cph these storings avoid recomp. of Ri_iwmix
140    cph DESPITE TAFs 'not necessary' warning!
141    CADJ STORE dbloc  = comlev1_kpp, key = ikppkey
142    CADJ STORE shsq   = comlev1_kpp, key = ikppkey
143    CADJ STORE ghat   = comlev1_kpp, key = ikppkey
144    CADJ STORE diffus = comlev1_kpp, key = ikppkey
145    cph)
146    
147  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
148  c set seafloor values to zero and fill extra "Nrp1" coefficients  c set seafloor values to zero and fill extra "Nrp1" coefficients
149  c for blmix  c for blmix
# Line 146  c diagnose the new boundary layer depth Line 164  c diagnose the new boundary layer depth
164  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
165    
166        call bldepth (        call bldepth (
167       I       kmtj       I       mytime, mythid
168         I     , kmtj
169       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol
170       I     , ikey       I     , ikppkey
171       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma
172       &     )       &     )
173    
174  CADJ STORE hbl,bfsfc,stable,casea,kbl = comlev1_kpp, key = ikey  CADJ STORE hbl,bfsfc,stable,casea,kbl = comlev1_kpp, key = ikppkey
175    
176  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
177  c compute boundary layer diffusivities  c compute boundary layer diffusivities
# Line 160  c--------------------------------------- Line 179  c---------------------------------------
179    
180        call blmix (        call blmix (
181       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl
182       O     , dkm1, blmc, ghat, sigma       O     , dkm1, blmc, ghat, sigma, ikppkey
183       &     )       &     )
184    cph(
185  CADJ STORE dkm1,blmc,ghat = comlev1_kpp, key = ikey  CADJ STORE dkm1,blmc,ghat = comlev1_kpp, key = ikppkey
186    CADJ STORE hbl, kbl, diffus, casea = comlev1_kpp, key = ikppkey
187    cph)
188    
189  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
190  c enhance diffusivity at interface kbl - 1  c enhance diffusivity at interface kbl - 1
# Line 174  c--------------------------------------- Line 195  c---------------------------------------
195       U     , ghat       U     , ghat
196       O     , blmc )       O     , blmc )
197    
198    cph(
199    cph avoids recomp. of enhance
200    CADJ STORE blmc = comlev1_kpp, key = ikppkey
201    cph)
202    
203  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
204  c combine interior and boundary layer coefficients and nonlocal term  c combine interior and boundary layer coefficients and nonlocal term
205    c !!!NOTE!!! In shallow (2-level) regions and for shallow mixed layers
206    c (< 1 level), diffusivity blmc can become negative.  The max-s below
207    c are a hack until this problem is properly diagnosed and fixed.
208  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
   
209        do k = 1, Nr        do k = 1, Nr
210           do i = 1, imt           do i = 1, imt
211              if (k .lt. kbl(i)) then              if (k .lt. kbl(i)) then
212                 do md = 1, mdiff                 diffus(i,k,1) = max ( blmc(i,k,1), viscAr  )
213                    diffus(i,k,md) = blmc(i,k,md)                 diffus(i,k,2) = max ( blmc(i,k,2), diffKrNrS(k) )
214                 end do                 diffus(i,k,3) = max ( blmc(i,k,3), diffKrNrT(k) )
215              else              else
216                 ghat(i,k) = 0.                 ghat(i,k) = 0.
217              endif              endif
# Line 198  c--------------------------------------- Line 226  c---------------------------------------
226  c*************************************************************************  c*************************************************************************
227    
228        subroutine bldepth (        subroutine bldepth (
229       I       kmtj       I       mytime, mythid
230         I     , kmtj
231       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol
232       I     , ikey       I     , ikppkey
233       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma
234       &     )       &     )
235    
# Line 235  c Line 264  c
264    
265  c input  c input
266  c------  c------
267    c myTime    : current time in simulation
268    c myThid    : thread number for this instance of the routine
269  c kmtj      : number of vertical layers  c kmtj      : number of vertical layers
270  c dvsq      : (velocity shear re sfc)^2             ((m/s)^2)  c dvsq      : (velocity shear re sfc)^2             ((m/s)^2)
271  c dbloc     : local delta buoyancy across interfaces  (m/s^2)  c dbloc     : local delta buoyancy across interfaces  (m/s^2)
# Line 245  c ustar     : surface friction velocity Line 276  c ustar     : surface friction velocity
276  c bo        : surface turbulent buoyancy forcing    (m^2/s^3)  c bo        : surface turbulent buoyancy forcing    (m^2/s^3)
277  c bosol     : radiative buoyancy forcing            (m^2/s^3)  c bosol     : radiative buoyancy forcing            (m^2/s^3)
278  c coriol    : Coriolis parameter                        (1/s)  c coriol    : Coriolis parameter                        (1/s)
279          _RL     mytime
280          integer mythid
281        integer kmtj(imt)        integer kmtj(imt)
282        _RS dvsq   (imt,Nr)        _KPP_RL dvsq  (imt,Nr)
283        _RS dbloc  (imt,Nr)        _KPP_RL dbloc (imt,Nr)
284        _RS Ritop  (imt,Nr)        _KPP_RL Ritop (imt,Nr)
285        _RS ustar  (imt)        _KPP_RL ustar (imt)
286        _RS bo     (imt)        _KPP_RL bo    (imt)
287        _RS bosol  (imt)        _KPP_RL bosol (imt)
288        _RS coriol (imt)        _KPP_RL coriol(imt)
289        integer ikey        integer ikppkey, kkppkey
290    
291  c  output  c  output
292  c--------  c--------
# Line 264  c casea     : =1 in case A, =0 in case B Line 297  c casea     : =1 in case A, =0 in case B
297  c kbl       : -1 of first grid level below hbl  c kbl       : -1 of first grid level below hbl
298  c Rib       : Bulk Richardson number  c Rib       : Bulk Richardson number
299  c sigma     : normalized depth (d/hbl)  c sigma     : normalized depth (d/hbl)
300        _RS hbl    (imt)        _KPP_RL hbl   (imt)
301        _RS bfsfc  (imt)        _KPP_RL bfsfc (imt)
302        _RS stable (imt)        _KPP_RL stable(imt)
303        _RS casea  (imt)        _KPP_RL casea (imt)
304        integer kbl (imt)        integer kbl   (imt)
305        _RS Rib    (imt,Nr)        _KPP_RL Rib   (imt,Nr)
306        _RS sigma  (imt)        _KPP_RL sigma (imt)
307    
308  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
309    
310  c  local  c  local
311  c-------  c-------
312  c wm, ws    : turbulent velocity scales         (m/s)  c wm, ws    : turbulent velocity scales         (m/s)
313  c zwork     : depth work array        _KPP_RL wm(imt), ws(imt)
314        _RS wm(imt), ws(imt)        _RL     worka(imt)
       _RS zwork(imt)  
315    
316        _RS bvsq, vtsq, hekman, hmonob, hlimit        _KPP_RL bvsq, vtsq, hekman, hmonob, hlimit, tempVar1, tempVar2
317        integer i, kl        integer i, kl
318    
319        _RS        p5    , eins    , m1        _KPP_RL     p5    , eins
320        parameter (p5=0.5, eins=1.0, m1=-1.0 )        parameter ( p5=0.5, eins=1.0 )
321          _RL         minusone
322  crg intermediate result in RL to avoid overflow in adjoint        parameter ( minusone=-1.0 )
        _RL        dpshear2  
 #ifdef KPP_TEST_DENOM  
        _RL        dhelp  
 #endif  
323    
324  c find bulk Richardson number at every grid level until > Ricr  c find bulk Richardson number at every grid level until > Ricr
325  c  c
# Line 312  c     initialize hbl and kbl to bottomed Line 340  c     initialize hbl and kbl to bottomed
340    
341        do kl = 2, Nr        do kl = 2, Nr
342    
343    #ifdef ALLOW_AUTODIFF_TAMC
344             kkppkey = (ikppkey-1)*Nr + kl
345    #endif
346    
347  c     compute bfsfc = sw fraction at hbf * zgrid  c     compute bfsfc = sw fraction at hbf * zgrid
348    
349           do i = 1, imt           do i = 1, imt
350              zwork(i) = zgrid(kl)              worka(i) = zgrid(kl)
351           end do           end do
352    CADJ store worka = comlev1_kpp_k, key = kkppkey
353           call SWFRAC(           call SWFRAC(
354       I        imt, hbf, zwork,       I        imt, hbf,
355       O        bfsfc)       I        mytime, mythid,
356         U        worka )
357    CADJ store worka = comlev1_kpp_k, key = kkppkey
358    
359    
360           do i = 1, imt           do i = 1, imt
361    
# Line 329  c     use caseA as temporary array Line 365  c     use caseA as temporary array
365    
366  c     compute bfsfc= Bo + radiative contribution down to hbf * hbl  c     compute bfsfc= Bo + radiative contribution down to hbf * hbl
367    
368              bfsfc(i) = bo(i) + bosol(i)*(1. - bfsfc(i))              bfsfc(i) = bo(i) + bosol(i)*(1. - worka(i))
369              stable(i) = p5 + sign(p5,bfsfc(i))              stable(i) = p5 + sign(p5,bfsfc(i))
370              sigma(i) = stable(i) + (1. - stable(i)) * epsilon              sigma(i) = stable(i) + (1. - stable(i)) * epsilon
371    
# Line 342  c--------------------------------------- Line 378  c---------------------------------------
378           call wscale (           call wscale (
379       I        sigma, casea, ustar, bfsfc,       I        sigma, casea, ustar, bfsfc,
380       O        wm, ws )       O        wm, ws )
381    CADJ store ws = comlev1_kpp_k, key = kkppkey
382    
383           do i = 1, imt           do i = 1, imt
384    
# Line 370  c Line 407  c
407  c     rg: assignment to double precision variable to avoid overflow  c     rg: assignment to double precision variable to avoid overflow
408  c     ph: test for zero nominator  c     ph: test for zero nominator
409  c  c
410  #ifdef KPP_TEST_DENOM  
411  ctl replace              tempVar1  = dvsq(i,kl) + vtsq          
412              dpshear2  = dvsq(i,kl) + vtsq                        tempVar2 = max(tempVar1, phepsi)
413              dpshear2 = max(dpshear2, phepsi)              Rib(i,kl) = Ritop(i,kl) / tempVar2
414              Rib(i,kl) = Ritop(i,kl) / dpshear2  
 ctl end-replace  
 #else  
             dpshear2  = dvsq(i,kl) + vtsq + epsln  
             Rib(i,kl) = Ritop(i,kl) / dpshear2  
 #endif  
415           end do           end do
416        end do        end do
417          
418    cph(
419    cph  without this store, there is a recomputation error for
420    cph  rib in adbldepth (probably partial recomputation problem)    
421    CADJ store Rib = comlev1_kpp
422    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy),Nr /)
423    cph)
424    
425        do kl = 2, Nr        do kl = 2, Nr
426           do i = 1, imt           do i = 1, imt
427              if (kbl(i).eq.kmtj(i) .and. Rib(i,kl).gt.Ricr) kbl(i) = kl              if (kbl(i).eq.kmtj(i) .and. Rib(i,kl).gt.Ricr) kbl(i) = kl
# Line 390  ctl end-replace Line 429  ctl end-replace
429        end do        end do
430    
431  CADJ store kbl = comlev1_kpp  CADJ store kbl = comlev1_kpp
432  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
433    
434        do i = 1, imt        do i = 1, imt
435           kl = kbl(i)           kl = kbl(i)
   
436  c     linearly interpolate to find hbl where Rib = Ricr  c     linearly interpolate to find hbl where Rib = Ricr
   
437           if (kl.gt.1 .and. kl.lt.kmtj(i)) then           if (kl.gt.1 .and. kl.lt.kmtj(i)) then
438  #ifdef KPP_TEST_DENOM              tempVar1 = (Rib(i,kl)-Rib(i,kl-1))
             dhelp = (Rib(i,kl)-Rib(i,kl-1))  
             hbl(i) = -zgrid(kl-1) + (zgrid(kl-1)-zgrid(kl)) *  
      1           (Ricr - Rib(i,kl-1)) / dhelp  
 #else  
439              hbl(i) = -zgrid(kl-1) + (zgrid(kl-1)-zgrid(kl)) *              hbl(i) = -zgrid(kl-1) + (zgrid(kl-1)-zgrid(kl)) *
440       1           (Ricr - Rib(i,kl-1)) / (Rib(i,kl)-Rib(i,kl-1))       1           (Ricr - Rib(i,kl-1)) / tempVar1
 #endif  
441           endif           endif
442        end do        end do
443    
444  CADJ store hbl = comlev1_kpp  CADJ store hbl = comlev1_kpp
445  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
446    
447  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
448  c     find stability and buoyancy forcing for boundary layer  c     find stability and buoyancy forcing for boundary layer
449  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
450    
451          do i = 1, imt
452             worka(i) = hbl(i)
453          end do
454    CADJ store worka = comlev1_kpp
455    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
456        call SWFRAC(        call SWFRAC(
457       I     imt, m1, hbl,       I     imt, minusone,
458       O     bfsfc)       I     mytime, mythid,
459         U     worka )
460    CADJ store worka = comlev1_kpp
461    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
462    
463          do i = 1, imt
464             bfsfc(i)  = bo(i) + bosol(i) * (1. - worka(i))
465          end do
466  CADJ store bfsfc = comlev1_kpp  CADJ store bfsfc = comlev1_kpp
467  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
468    
469    c--   ensure bfsfc is never 0
470        do i = 1, imt        do i = 1, imt
          bfsfc(i)  = bo(i) + bosol(i) * (1. - bfsfc(i))  
471           stable(i) = p5 + sign( p5, bfsfc(i) )           stable(i) = p5 + sign( p5, bfsfc(i) )
472             bfsfc(i) = sign(eins,bfsfc(i))*max(phepsi,abs(bfsfc(i)))
 #ifdef KPP_TEST_DENOM  
 ctl add  
         bfsfc(i) = sign(eins,bfsfc(i))*max(phepsi,abs(bfsfc(i)))  
 ctl end-add  
 #else  
 c--      ensures bfsfc is never 0  
          bfsfc(i)  = bfsfc(i) + stable(i) * epsln  
 #endif  
   
473        end do        end do
474    
475    cph(
476    cph  added stable to store list to avoid extensive recomp.
477    CADJ store bfsfc, stable = comlev1_kpp
478    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
479    cph)
480    
481  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
482  c check hbl limits for hekman or hmonob  c check hbl limits for hekman or hmonob
483  c     ph: test for zero nominator  c     ph: test for zero nominator
# Line 445  c--------------------------------------- Line 485  c---------------------------------------
485    
486        do i = 1, imt        do i = 1, imt
487           if (bfsfc(i) .gt. 0.0) then           if (bfsfc(i) .gt. 0.0) then
 #ifdef KPP_TEST_DENOM  
 ctl replace  
488              hekman = cekman * ustar(i) / max(abs(Coriol(i)),phepsi)              hekman = cekman * ustar(i) / max(abs(Coriol(i)),phepsi)
 ctl end-replace  
 #else  
             hekman = cekman * ustar(i) / (abs(coriol(i))+epsln)  
 #endif  
489              hmonob = cmonob * ustar(i)*ustar(i)*ustar(i)              hmonob = cmonob * ustar(i)*ustar(i)*ustar(i)
490       &           / vonk / bfsfc(i)       &           / vonk / bfsfc(i)
491              hlimit = stable(i) * min(hekman,hmonob)              hlimit = stable(i) * min(hekman,hmonob)
492       &             + (stable(i)-1.) * zgrid(Nr)       &             + (stable(i)-1.) * zgrid(Nr)
493              hbl(i) = min(hbl(i),hlimit)              hbl(i) = min(hbl(i),hlimit)
494              hbl(i) = max(hbl(i),minKPPhbl)           end if
495           endif        end do
496    CADJ store hbl = comlev1_kpp
497    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
498    
499          do i = 1, imt
500             hbl(i) = max(hbl(i),minKPPhbl)
501           kbl(i) = kmtj(i)           kbl(i) = kmtj(i)
502        end do        end do
503    
504  CADJ store hbl = comlev1_kpp  CADJ store hbl = comlev1_kpp
505  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
506    
507  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
508  c      find new kbl  c      find new kbl
# Line 481  c--------------------------------------- Line 520  c---------------------------------------
520  c      find stability and buoyancy forcing for final hbl values  c      find stability and buoyancy forcing for final hbl values
521  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
522    
523          do i = 1, imt
524             worka(i) = hbl(i)
525          end do
526    CADJ store worka = comlev1_kpp
527    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
528        call SWFRAC(        call SWFRAC(
529       I     imt, m1, hbl,       I     imt, minusone,
530       O     bfsfc)       I     mytime, mythid,
531         U     worka )
532    CADJ store worka = comlev1_kpp
533    CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
534    
535          do i = 1, imt
536             bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i))
537          end do
538  CADJ store bfsfc = comlev1_kpp  CADJ store bfsfc = comlev1_kpp
539  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
540    
541    c--   ensures bfsfc is never 0
542        do i = 1, imt        do i = 1, imt
          bfsfc(i) = bo(i) + bosol(i) * (1. - bfsfc(i))  
543           stable(i) = p5 + sign( p5, bfsfc(i) )           stable(i) = p5 + sign( p5, bfsfc(i) )
544  #ifdef KPP_TEST_DENOM           bfsfc(i) = sign(eins,bfsfc(i))*max(phepsi,abs(bfsfc(i)))
 ctl add  
         bfsfc(i) = sign(eins,bfsfc(i))*max(phepsi,abs(bfsfc(i)))  
 ctl end-add  
 #else  
 c--      ensures bfsfc is never 0  
          bfsfc(i)  = bfsfc(i) + stable(i) * epsln  
 #endif  
545        end do        end do
546    
547  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 540  c sigma   : normalized depth (d/hbl) Line 583  c sigma   : normalized depth (d/hbl)
583  c hbl     : boundary layer depth (m)  c hbl     : boundary layer depth (m)
584  c ustar   : surface friction velocity         (m/s)  c ustar   : surface friction velocity         (m/s)
585  c bfsfc   : total surface buoyancy flux       (m^2/s^3)  c bfsfc   : total surface buoyancy flux       (m^2/s^3)
586        _RS sigma(imt)        _KPP_RL sigma(imt)
587        _RS hbl  (imt)        _KPP_RL hbl  (imt)
588        _RS ustar(imt)        _KPP_RL ustar(imt)
589        _RS bfsfc(imt)        _KPP_RL bfsfc(imt)
590    
591  c  output  c  output
592  c--------  c--------
593  c wm, ws  : turbulent velocity scales at sigma  c wm, ws  : turbulent velocity scales at sigma
594        _RS wm(imt), ws(imt)        _KPP_RL wm(imt), ws(imt)
595    
596  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
597    
598  c local  c local
599  c------  c------
600  c zehat   : = zeta *  ustar**3  c zehat   : = zeta *  ustar**3
601        _RS zehat        _KPP_RL zehat
602    
603        integer iz, izp1, ju, i, jup1        integer iz, izp1, ju, i, jup1
604        _RS udiff, zdiff, zfrac, ufrac, fzfrac, wam, wbm, was, wbs, u3        _KPP_RL udiff, zdiff, zfrac, ufrac, fzfrac, wam
605        _RL dum        _KPP_RL wbm, was, wbs, u3, tempVar
606    
607  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
608  c use lookup table for zehat < zmax only; otherwise use  c use lookup table for zehat < zmax only; otherwise use
# Line 576  c--------------------------------------- Line 619  c---------------------------------------
619              iz    = min( iz, nni )              iz    = min( iz, nni )
620              iz    = max( iz, 0 )              iz    = max( iz, 0 )
621              izp1  = iz + 1              izp1  = iz + 1
622    
623              udiff = ustar(i) - umin              udiff = ustar(i) - umin
624              ju    = int( udiff / deltau )              ju    = int( udiff / deltau )
625              ju    = min( ju, nnj )              ju    = min( ju, nnj )
626              ju    = max( ju, 0 )              ju    = max( ju, 0 )
627              jup1  = ju + 1              jup1  = ju + 1
628    
629              zfrac = zdiff / deltaz - float(iz)              zfrac = zdiff / deltaz - float(iz)
630              ufrac = udiff / deltau - float(ju)              ufrac = udiff / deltau - float(ju)
631    
632              fzfrac= 1. - zfrac              fzfrac= 1. - zfrac
633              wam   = fzfrac     * wmt(iz,jup1) + zfrac * wmt(izp1,jup1)              wam   = fzfrac     * wmt(iz,jup1) + zfrac * wmt(izp1,jup1)
634              wbm   = fzfrac     * wmt(iz,ju  ) + zfrac * wmt(izp1,ju  )              wbm   = fzfrac     * wmt(iz,ju  ) + zfrac * wmt(izp1,ju  )
635              wm(i) = (1.-ufrac) * wbm          + ufrac * wam              wm(i) = (1.-ufrac) * wbm          + ufrac * wam
636    
637              was   = fzfrac     * wst(iz,jup1) + zfrac * wst(izp1,jup1)              was   = fzfrac     * wst(iz,jup1) + zfrac * wst(izp1,jup1)
638              wbs   = fzfrac     * wst(iz,ju  ) + zfrac * wst(izp1,ju  )              wbs   = fzfrac     * wst(iz,ju  ) + zfrac * wst(izp1,ju  )
639              ws(i) = (1.-ufrac) * wbs          + ufrac * was              ws(i) = (1.-ufrac) * wbs          + ufrac * was
640    
641           else           else
642    
643              u3    = ustar(i) * ustar(i) * ustar(i)              u3 = ustar(i) * ustar(i) * ustar(i)
644              dum   = u3 + conc1 * zehat              tempVar = u3 + conc1 * zehat
645              wm(i) = vonk     * ustar(i) * u3 / dum              wm(i) = vonk * ustar(i) * u3 / tempVar
646              ws(i) = wm(i)              ws(i) = wm(i)
647    
648           endif           endif
649    
650        end do        end do
# Line 615  c*************************************** Line 658  c***************************************
658    
659        subroutine Ri_iwmix (        subroutine Ri_iwmix (
660       I       kmtj, shsq, dbloc, dblocSm       I       kmtj, shsq, dbloc, dblocSm
661       I     , ikey       I     , ikppkey
662       O     , diffus )       O     , diffus )
663    
664  c     compute interior viscosity diffusivity coefficients due  c     compute interior viscosity diffusivity coefficients due
# Line 636  c     shsq   (imt,Nr)      (local veloci Line 679  c     shsq   (imt,Nr)      (local veloci
679  c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)  c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)
680  c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)  c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)
681        integer kmtj   (imt)        integer kmtj   (imt)
682        _RS     shsq   (imt,Nr)        _KPP_RL shsq   (imt,Nr)
683        _RS     dbloc  (imt,Nr)        _KPP_RL dbloc  (imt,Nr)
684        _RS     dblocSm(imt,Nr)        _KPP_RL dblocSm(imt,Nr)
685        integer ikey        integer ikppkey
686    
687  c output  c output
688  c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)  c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)
689  c     diffus(imt,0:Nrp1,2)  vertical scalar diffusivity             (m^2/s)  c     diffus(imt,0:Nrp1,2)  vertical scalar diffusivity             (m^2/s)
690  c     diffus(imt,0:Nrp1,3)  vertical temperature diffusivity        (m^2/s)  c     diffus(imt,0:Nrp1,3)  vertical temperature diffusivity        (m^2/s)
691        _RS diffus(imt,0:Nrp1,3)        _KPP_RL diffus(imt,0:Nrp1,3)
692    
693  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
694    
695  c local variables  c local variables
696  c     Rig                   local Richardson number  c     Rig                   local Richardson number
697  c     fRi, fcon             function of Rig  c     fRi, fcon             function of Rig
698        _RS Rig        _KPP_RL Rig
699        _RS fRi, fcon        _KPP_RL fRi, fcon
700        _RS ratio        _KPP_RL ratio
701        integer i, ki, mr        integer i, ki
702        _RS c1, c0        _KPP_RL c1, c0
703    
704    #ifdef ALLOW_KPP_VERTICALLY_SMOOTH
705          integer mr
706    CADJ INIT kpp_ri_tape_mr = common, 1
707    #endif
708    
709  c constants  c constants
710        c1 = 1.0        c1 = 1.0
# Line 667  c     compute interior gradient Ri at al Line 715  c     compute interior gradient Ri at al
715  c     use diffus(*,*,1) as temporary storage of Ri to be smoothed  c     use diffus(*,*,1) as temporary storage of Ri to be smoothed
716  c     use diffus(*,*,2) as temporary storage for Brunt-Vaisala squared  c     use diffus(*,*,2) as temporary storage for Brunt-Vaisala squared
717  c     set values at bottom and below to nearest value above bottom  c     set values at bottom and below to nearest value above bottom
718    #ifdef ALLOW_AUTODIFF_TAMC
719    C     break data flow dependence on diffus
720          diffus(1,1,1) = 0.0
721    
722        do ki = 1, Nr        do ki = 1, Nr
723           do i = 1, imt           do i = 1, imt
724              if     (kmtj(i) .EQ. 0      ) then              diffus(i,ki,1) = 0.
725                diffus(i,ki,2) = 0.
726                diffus(i,ki,3) = 0.
727             enddo
728          enddo
729    #endif
730                
731    
732          do ki = 1, Nr
733             do i = 1, imt
734                if     (kmtj(i) .LE. 1      ) then
735                 diffus(i,ki,1) = 0.                 diffus(i,ki,1) = 0.
736                 diffus(i,ki,2) = 0.                 diffus(i,ki,2) = 0.
737              elseif (ki      .GE. kmtj(i)) then              elseif (ki      .GE. kmtj(i)) then
# Line 678  c     set values at bottom and below to Line 739  c     set values at bottom and below to
739                 diffus(i,ki,2) = diffus(i,ki-1,2)                 diffus(i,ki,2) = diffus(i,ki-1,2)
740              else              else
741                 diffus(i,ki,1) = dblocSm(i,ki) * (zgrid(ki)-zgrid(ki+1))                 diffus(i,ki,1) = dblocSm(i,ki) * (zgrid(ki)-zgrid(ki+1))
 #ifdef KPP_TEST_DENOM  
742       &              / max( Shsq(i,ki), phepsi )       &              / max( Shsq(i,ki), phepsi )
 #else  
      &              / ( shsq(i,ki) + epsln )  
 #endif  
743                 diffus(i,ki,2) = dbloc(i,ki)   / (zgrid(ki)-zgrid(ki+1))                 diffus(i,ki,2) = dbloc(i,ki)   / (zgrid(ki)-zgrid(ki+1))
744              endif              endif
745           end do           end do
746        end do        end do
747    CADJ store diffus = comlev1_kpp, key = ikppkey
748    
749  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
750  c     vertically smooth Ri  c     vertically smooth Ri
751    #ifdef ALLOW_KPP_VERTICALLY_SMOOTH
752        do mr = 1, num_v_smooth_Ri        do mr = 1, num_v_smooth_Ri
753    
754  CADJ store diffus(:,:,1) = comlev1_kpp_sm  CADJ store diffus(:,:,1) = kpp_ri_tape_mr
755  CADJ &  , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2 /)  CADJ &  , key=mr, shape=(/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2 /)
756    
757           call z121 (           call z121 (
758       U     diffus(1,0,1))       U     diffus(1,0,1))
759        end do        end do
760    #endif
 CADJ store diffus = comlev1_kpp  
 CADJ &  , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2,3 /)  
761    
762  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
763  c                           after smoothing loop  c                           after smoothing loop
# Line 726  c  evaluate f of smooth Ri for shear ins Line 782  c  evaluate f of smooth Ri for shear ins
782  c ----------------------------------------------------------------------  c ----------------------------------------------------------------------
783  c            evaluate diffusivities and viscosity  c            evaluate diffusivities and viscosity
784  c    mixing due to internal waves, and shear and static instability  c    mixing due to internal waves, and shear and static instability
785    
786    #ifndef EXCLUDE_KPP_SHEAR_MIX
787              diffus(i,ki,1) = viscAr + fcon * difmcon + fRi * difm0              diffus(i,ki,1) = viscAr + fcon * difmcon + fRi * difm0
788              diffus(i,ki,2) = diffKrS + fcon * difscon + fRi * difs0              diffus(i,ki,2) = diffKrNrS(ki)+ fcon*difscon + fRi*difs0
789              diffus(i,ki,3) = diffKrT + fcon * difscon + fRi * difs0              diffus(i,ki,3) = diffKrNrT(ki)+ fcon*difscon + fRi*difs0
790    #else
791                diffus(i,ki,1) = viscAr
792                diffus(i,ki,2) = diffKrNrS(ki)
793                diffus(i,ki,3) = diffKrNrT(ki)
794    #endif
795    
796           end do           end do
797        end do        end do
# Line 769  c     penetrative convection. Line 831  c     penetrative convection.
831  c input/output  c input/output
832  c-------------  c-------------
833  c v     : 2-D array to be smoothed in Nrp1 direction  c v     : 2-D array to be smoothed in Nrp1 direction
834        _RS v(imt,0:Nrp1)        _KPP_RL v(imt,0:Nrp1)
835    
836  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
837    
838  c local  c local
839        _RS zwork, zflag        _KPP_RL zwork, zflag
840        _RS KRi_range(1:Nrp1)        _KPP_RL KRi_range(1:Nrp1)
841        integer i, k, km1, kp1        integer i, k, km1, kp1
842    
843        _RS         p0      , p25       , p5      , p2        _KPP_RL     p0      , p25       , p5      , p2
844        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 )
845    
846        KRi_range(Nrp1) = p0        KRi_range(Nrp1) = p0
# Line 789  C--   dummy assignment to end declaratio Line 851  C--   dummy assignment to end declaratio
851    
852  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
853  CHPF$ INDEPENDENT  CHPF$ INDEPENDENT
854    CADJ INIT z121tape = common, Nr
855  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
856    
857        do i = 1, imt        do i = 1, imt
858    
859             k = 1
860    CADJ STORE v(i,k) = z121tape
861           v(i,Nrp1) = v(i,Nr)           v(i,Nrp1) = v(i,Nr)
862    
863           do k = 1, Nr           do k = 1, Nr
# Line 806  CHPF$ INDEPENDENT Line 872  CHPF$ INDEPENDENT
872           zflag  = p2 + KRi_range(1) * KRi_range(2)           zflag  = p2 + KRi_range(1) * KRi_range(2)
873           v(i,1) = v(i,1) / zflag           v(i,1) = v(i,1) / zflag
874    
 CADJ INIT z121tape = common, Nr  
875           do k = 2, Nr           do k = 2, Nr
876  CADJ STORE v(i,k), zwork = z121tape  CADJ STORE v(i,k), zwork = z121tape
877              km1 = k - 1              km1 = k - 1
# Line 829  CADJ STORE v(i,k), zwork = z121tape Line 894  CADJ STORE v(i,k), zwork = z121tape
894    
895  c*************************************************************************  c*************************************************************************
896    
897        subroutine smooth_horiz_rs (        subroutine kpp_smooth_horiz (
898       I     k, bi, bj,       I     k, bi, bj,
899       I     fld_in,       U     fld )
      O     fld_out )  
900    
901  c     Apply horizontal smoothing to RS 2-D array  c     Apply horizontal smoothing to KPP array
902    
903        IMPLICIT NONE        IMPLICIT NONE
904  #include "SIZE.h"  #include "SIZE.h"
905    #include "GRID.h"
906  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
907    
908  c     input  c     input
909  c     k, bi, bj : array indices  c     bi, bj : array indices
910  c     fld_in    : 2-D array to be smoothed  c     k      : vertical index used for masking
911        integer k, bi, bj        integer k, bi, bj
       _RS fld_in(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
912    
913  c     output  c     input/output
914  c     fld_out   : smoothed 2-D array  c     fld    : 2-D array to be smoothed
915        _RS fld_out(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _KPP_RL fld( ibot:itop, jbot:jtop )
916    
917  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
918    
919  c     local  c     local
920        integer i, j, im1, ip1, jm1, jp1        integer i, j, im1, ip1, jm1, jp1
921        _RS tempVar        _KPP_RL tempVar
922        _RS fld_tmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _KPP_RL fld_tmp( ibot:itop, jbot:jtop )
   
       integer    imin        , imax  
       parameter( imin=1-OLx+1, imax=sNx+OLx-1 )  
923    
924        integer    jmin        , jmax        integer    imin       , imax       , jmin       , jmax
925        parameter( jmin=1-OLy+1, jmax=sNy+OLy-1 )        parameter( imin=ibot+1, imax=itop-1, jmin=jbot+1, jmax=jtop-1 )
926    
927        _RS        p0    , p5    , p25     , p125      , p0625        _KPP_RL    p0    , p5    , p25     , p125      , p0625
928        parameter( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 )        parameter( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 )
929    
930        DO j = jmin, jmax        DO j = jmin, jmax
# Line 873  c     local Line 934  c     local
934              im1 = i-1              im1 = i-1
935              ip1 = i+1              ip1 = i+1
936              tempVar =              tempVar =
937       &           p25   *   pMask(i  ,j  ,k,bi,bj)   +       &           p25   *   maskC(i  ,j  ,k,bi,bj)   +
938       &           p125  * ( pMask(im1,j  ,k,bi,bj)   +       &           p125  * ( maskC(im1,j  ,k,bi,bj)   +
939       &                     pMask(ip1,j  ,k,bi,bj)   +       &                     maskC(ip1,j  ,k,bi,bj)   +
940       &                     pMask(i  ,jm1,k,bi,bj)   +       &                     maskC(i  ,jm1,k,bi,bj)   +
941       &                     pMask(i  ,jp1,k,bi,bj) ) +       &                     maskC(i  ,jp1,k,bi,bj) ) +
942       &           p0625 * ( pMask(im1,jm1,k,bi,bj)   +       &           p0625 * ( maskC(im1,jm1,k,bi,bj)   +
943       &                     pMask(im1,jp1,k,bi,bj)   +       &                     maskC(im1,jp1,k,bi,bj)   +
944       &                     pMask(ip1,jm1,k,bi,bj)   +       &                     maskC(ip1,jm1,k,bi,bj)   +
945       &                     pMask(ip1,jp1,k,bi,bj) )       &                     maskC(ip1,jp1,k,bi,bj) )
946              IF ( tempVar .GE. p25 ) THEN              IF ( tempVar .GE. p25 ) THEN
947                 fld_tmp(i,j) = (                 fld_tmp(i,j) = (
948       &              p25  * fld_in(i  ,j  )*pMask(i  ,j  ,k,bi,bj) +       &              p25  * fld(i  ,j  )*maskC(i  ,j  ,k,bi,bj) +
949       &              p125 *(fld_in(im1,j  )*pMask(im1,j  ,k,bi,bj) +       &              p125 *(fld(im1,j  )*maskC(im1,j  ,k,bi,bj) +
950       &                     fld_in(ip1,j  )*pMask(ip1,j  ,k,bi,bj) +       &                     fld(ip1,j  )*maskC(ip1,j  ,k,bi,bj) +
951       &                     fld_in(i  ,jm1)*pMask(i  ,jm1,k,bi,bj) +       &                     fld(i  ,jm1)*maskC(i  ,jm1,k,bi,bj) +
952       &                     fld_in(i  ,jp1)*pMask(i  ,jp1,k,bi,bj))+       &                     fld(i  ,jp1)*maskC(i  ,jp1,k,bi,bj))+
953       &              p0625*(fld_in(im1,jm1)*pMask(im1,jm1,k,bi,bj) +       &              p0625*(fld(im1,jm1)*maskC(im1,jm1,k,bi,bj) +
954       &                     fld_in(im1,jp1)*pMask(im1,jp1,k,bi,bj) +       &                     fld(im1,jp1)*maskC(im1,jp1,k,bi,bj) +
955       &                     fld_in(ip1,jm1)*pMask(ip1,jm1,k,bi,bj) +       &                     fld(ip1,jm1)*maskC(ip1,jm1,k,bi,bj) +
956       &                     fld_in(ip1,jp1)*pMask(ip1,jp1,k,bi,bj)))       &                     fld(ip1,jp1)*maskC(ip1,jp1,k,bi,bj)))
957       &              / tempVar       &              / tempVar
958              ELSE              ELSE
959                 fld_tmp(i,j) = fld_in(i,j)                 fld_tmp(i,j) = fld(i,j)
960              ENDIF              ENDIF
961           ENDDO           ENDDO
962        ENDDO        ENDDO
# Line 903  c     local Line 964  c     local
964  c     transfer smoothed field to output array  c     transfer smoothed field to output array
965        DO j = jmin, jmax        DO j = jmin, jmax
966           DO i = imin, imax           DO i = imin, imax
967              fld_out(i,j) = fld_tmp(i,j)              fld(i,j) = fld_tmp(i,j)
968           ENDDO           ENDDO
969        ENDDO        ENDDO
970    
 c     set output array edges to input field values  
       DO j = jmin, jmax  
          DO i = 1-OLx, imin-1  
             fld_out(i,j) = fld_in(i,j)  
          END DO  
          DO i = imax+1, sNx+OLx  
             fld_out(i,j) = fld_in(i,j)  
          END DO  
       END DO  
       DO i = 1-OLx, sNx+OLx  
          DO j = 1-OLy, jmin-1  
             fld_out(i,j) = fld_in(i,j)  
          END DO  
          DO j = jmax+1, sNy+OLy  
             fld_out(i,j) = fld_in(i,j)  
          END DO  
       END DO  
   
971  #endif /* ALLOW_KPP */  #endif /* ALLOW_KPP */
972    
973        return        return
# Line 932  c     set output array edges to input fi Line 975  c     set output array edges to input fi
975    
976  c*************************************************************************  c*************************************************************************
977    
978        subroutine smooth_horiz_rl (        subroutine smooth_horiz (
979       I     k, bi, bj,       I     k, bi, bj,
980       I     fld_in,       U     fld )
      O     fld_out )  
981    
982  c     Apply horizontal smoothing to RL 2-D array  c     Apply horizontal smoothing to global _RL 2-D array
983    
984        IMPLICIT NONE        IMPLICIT NONE
985  #include "SIZE.h"  #include "SIZE.h"
986    #include "GRID.h"
987  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
988    
989  c     input  c     input
990  c     k, bi, bj : array indices  c     bi, bj : array indices
991  c     fld_in    : 2-D array to be smoothed  c     k      : vertical index used for masking
992        integer k, bi, bj        integer k, bi, bj
       _RL fld_in(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
993    
994  c     output  c     input/output
995  c     fld_out   : smoothed 2-D array  c     fld    : 2-D array to be smoothed
996        _RL fld_out(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fld( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
997    
998  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
999    
1000  c     local  c     local
1001        integer i, j, im1, ip1, jm1, jp1        integer i, j, im1, ip1, jm1, jp1
1002        _RL tempVar        _RL tempVar
1003        _RL fld_tmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fld_tmp( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
   
       integer    imin        , imax  
       parameter( imin=1-OLx+1, imax=sNx+OLx-1 )  
1004    
1005        integer    jmin        , jmax        integer   imin      , imax          , jmin      , jmax
1006        parameter( jmin=1-OLy+1, jmax=sNy+OLy-1 )        parameter(imin=2-OLx, imax=sNx+OLx-1, jmin=2-OLy, jmax=sNy+OLy-1)
1007    
1008        _RS        p0    , p5    , p25     , p125      , p0625        _RL        p0    , p5    , p25     , p125      , p0625
1009        parameter( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 )        parameter( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 )
1010    
1011        DO j = jmin, jmax        DO j = jmin, jmax
# Line 976  c     local Line 1015  c     local
1015              im1 = i-1              im1 = i-1
1016              ip1 = i+1              ip1 = i+1
1017              tempVar =              tempVar =
1018       &           p25   *   pMask(i  ,j  ,k,bi,bj)   +       &           p25   *   maskC(i  ,j  ,k,bi,bj)   +
1019       &           p125  * ( pMask(im1,j  ,k,bi,bj)   +       &           p125  * ( maskC(im1,j  ,k,bi,bj)   +
1020       &                     pMask(ip1,j  ,k,bi,bj)   +       &                     maskC(ip1,j  ,k,bi,bj)   +
1021       &                     pMask(i  ,jm1,k,bi,bj)   +       &                     maskC(i  ,jm1,k,bi,bj)   +
1022       &                     pMask(i  ,jp1,k,bi,bj) ) +       &                     maskC(i  ,jp1,k,bi,bj) ) +
1023       &           p0625 * ( pMask(im1,jm1,k,bi,bj)   +       &           p0625 * ( maskC(im1,jm1,k,bi,bj)   +
1024       &                     pMask(im1,jp1,k,bi,bj)   +       &                     maskC(im1,jp1,k,bi,bj)   +
1025       &                     pMask(ip1,jm1,k,bi,bj)   +       &                     maskC(ip1,jm1,k,bi,bj)   +
1026       &                     pMask(ip1,jp1,k,bi,bj) )       &                     maskC(ip1,jp1,k,bi,bj) )
1027              IF ( tempVar .GE. p25 ) THEN              IF ( tempVar .GE. p25 ) THEN
1028                 fld_tmp(i,j) = (                 fld_tmp(i,j) = (
1029       &              p25  * fld_in(i  ,j  )*pMask(i  ,j  ,k,bi,bj) +       &              p25  * fld(i  ,j  )*maskC(i  ,j  ,k,bi,bj) +
1030       &              p125 *(fld_in(im1,j  )*pMask(im1,j  ,k,bi,bj) +       &              p125 *(fld(im1,j  )*maskC(im1,j  ,k,bi,bj) +
1031       &                     fld_in(ip1,j  )*pMask(ip1,j  ,k,bi,bj) +       &                     fld(ip1,j  )*maskC(ip1,j  ,k,bi,bj) +
1032       &                     fld_in(i  ,jm1)*pMask(i  ,jm1,k,bi,bj) +       &                     fld(i  ,jm1)*maskC(i  ,jm1,k,bi,bj) +
1033       &                     fld_in(i  ,jp1)*pMask(i  ,jp1,k,bi,bj))+       &                     fld(i  ,jp1)*maskC(i  ,jp1,k,bi,bj))+
1034       &              p0625*(fld_in(im1,jm1)*pMask(im1,jm1,k,bi,bj) +       &              p0625*(fld(im1,jm1)*maskC(im1,jm1,k,bi,bj) +
1035       &                     fld_in(im1,jp1)*pMask(im1,jp1,k,bi,bj) +       &                     fld(im1,jp1)*maskC(im1,jp1,k,bi,bj) +
1036       &                     fld_in(ip1,jm1)*pMask(ip1,jm1,k,bi,bj) +       &                     fld(ip1,jm1)*maskC(ip1,jm1,k,bi,bj) +
1037       &                     fld_in(ip1,jp1)*pMask(ip1,jp1,k,bi,bj)))       &                     fld(ip1,jp1)*maskC(ip1,jp1,k,bi,bj)))
1038       &              / tempVar       &              / tempVar
1039              ELSE              ELSE
1040                 fld_tmp(i,j) = fld_in(i,j)                 fld_tmp(i,j) = fld(i,j)
1041              ENDIF              ENDIF
1042           ENDDO           ENDDO
1043        ENDDO        ENDDO
# Line 1006  c     local Line 1045  c     local
1045  c     transfer smoothed field to output array  c     transfer smoothed field to output array
1046        DO j = jmin, jmax        DO j = jmin, jmax
1047           DO i = imin, imax           DO i = imin, imax
1048              fld_out(i,j) = fld_tmp(i,j)              fld(i,j) = fld_tmp(i,j)
1049           ENDDO           ENDDO
1050        ENDDO        ENDDO
1051    
 c     set output array edges to input field values  
       DO j = jmin, jmax  
          DO i = 1-OLx, imin-1  
             fld_out(i,j) = fld_in(i,j)  
          END DO  
          DO i = imax+1, sNx+OLx  
             fld_out(i,j) = fld_in(i,j)  
          END DO  
       END DO  
       DO i = 1-OLx, sNx+OLx  
          DO j = 1-OLy, jmin-1  
             fld_out(i,j) = fld_in(i,j)  
          END DO  
          DO j = jmax+1, sNy+OLy  
             fld_out(i,j) = fld_in(i,j)  
          END DO  
       END DO  
   
1052  #endif /* ALLOW_KPP */  #endif /* ALLOW_KPP */
1053    
1054        return        return
# Line 1037  c*************************************** Line 1058  c***************************************
1058    
1059        subroutine blmix (        subroutine blmix (
1060       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl
1061       O     , dkm1, blmc, ghat, sigma       O     , dkm1, blmc, ghat, sigma, ikppkey
1062       &     )       &     )
1063    
1064  c     mixing coefficients within boundary layer depend on surface  c     mixing coefficients within boundary layer depend on surface
# Line 1061  c     stable(imt)                 = 1 in Line 1082  c     stable(imt)                 = 1 in
1082  c     casea (imt)                 = 1 in case A  c     casea (imt)                 = 1 in case A
1083  c     diffus(imt,0:Nrp1,mdiff)    vertical diffusivities              (m^2/s)  c     diffus(imt,0:Nrp1,mdiff)    vertical diffusivities              (m^2/s)
1084  c     kbl(imt)                    -1 of first grid level below hbl  c     kbl(imt)                    -1 of first grid level below hbl
1085        _RS ustar (imt)        _KPP_RL ustar (imt)
1086        _RS bfsfc (imt)        _KPP_RL bfsfc (imt)
1087        _RS hbl   (imt)        _KPP_RL hbl   (imt)
1088        _RS stable(imt)        _KPP_RL stable(imt)
1089        _RS casea (imt)        _KPP_RL casea (imt)
1090        _RS diffus(imt,0:Nrp1,mdiff)        _KPP_RL diffus(imt,0:Nrp1,mdiff)
1091        integer kbl(imt)        integer kbl(imt)
1092    
1093  c output  c output
# Line 1074  c     dkm1 (imt,mdiff)            bounda Line 1095  c     dkm1 (imt,mdiff)            bounda
1095  c     blmc (imt,Nr,mdiff)         boundary layer mixing coefficients  (m^2/s)  c     blmc (imt,Nr,mdiff)         boundary layer mixing coefficients  (m^2/s)
1096  c     ghat (imt,Nr)               nonlocal scalar transport  c     ghat (imt,Nr)               nonlocal scalar transport
1097  c     sigma(imt)                  normalized depth (d / hbl)  c     sigma(imt)                  normalized depth (d / hbl)
1098        _RS dkm1 (imt,mdiff)        _KPP_RL dkm1 (imt,mdiff)
1099        _RS blmc (imt,Nr,mdiff)        _KPP_RL blmc (imt,Nr,mdiff)
1100        _RS ghat (imt,Nr)        _KPP_RL ghat (imt,Nr)
1101        _RS sigma(imt)        _KPP_RL sigma(imt)
1102          integer ikppkey, kkppkey
1103    
1104  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1105    
# Line 1085  c  local Line 1107  c  local
1107  c     gat1*(imt)                 shape function at sigma = 1  c     gat1*(imt)                 shape function at sigma = 1
1108  c     dat1*(imt)                 derivative of shape function at sigma = 1  c     dat1*(imt)                 derivative of shape function at sigma = 1
1109  c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)  c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)
1110        _RS gat1m(imt), gat1s(imt), gat1t(imt)        _KPP_RL gat1m(imt), gat1s(imt), gat1t(imt)
1111        _RS dat1m(imt), dat1s(imt), dat1t(imt)        _KPP_RL dat1m(imt), dat1s(imt), dat1t(imt)
1112        _RS ws(imt), wm(imt)        _KPP_RL ws(imt), wm(imt)
1113        integer i, kn, ki        integer i, kn, ki
1114        _RS     R, dvdzup, dvdzdn, viscp        _KPP_RL R, dvdzup, dvdzdn, viscp
1115        _RS     difsp, diftp, visch, difsh, difth        _KPP_RL difsp, diftp, visch, difsh, difth
1116        _RS     f1, sig, a1, a2, a3, delhat        _KPP_RL f1, sig, a1, a2, a3, delhat
1117        _RS     Gm, Gs, Gt        _KPP_RL Gm, Gs, Gt
1118        _RL     dum        _KPP_RL tempVar
1119    
1120        _RS        p0    , eins        _KPP_RL    p0    , eins
1121        parameter (p0=0.0, eins=1.0)        parameter (p0=0.0, eins=1.0)
1122    
1123  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 1106  c--------------------------------------- Line 1128  c---------------------------------------
1128           sigma(i) = stable(i) * 1.0 + (1. - stable(i)) * epsilon           sigma(i) = stable(i) * 1.0 + (1. - stable(i)) * epsilon
1129        end do        end do
1130    
1131    CADJ STORE sigma = comlev1_kpp, key = ikppkey
1132        call wscale (        call wscale (
1133       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1134       O        wm, ws )       O        wm, ws )
1135    CADJ STORE wm = comlev1_kpp, key = ikppkey
1136    CADJ STORE ws = comlev1_kpp, key = ikppkey
1137    
1138          do i = 1, imt
1139             wm(i) = sign(eins,wm(i))*max(phepsi,abs(wm(i)))
1140             ws(i) = sign(eins,ws(i))*max(phepsi,abs(ws(i)))
1141          end do
1142    CADJ STORE wm = comlev1_kpp, key = ikppkey
1143    CADJ STORE ws = comlev1_kpp, key = ikppkey
1144    
1145        do i = 1, imt        do i = 1, imt
1146    
# Line 1140  c--------------------------------------- Line 1172  c---------------------------------------
1172           difsh  = diffus(i,kn,2) + difsp * delhat           difsh  = diffus(i,kn,2) + difsp * delhat
1173           difth  = diffus(i,kn,3) + diftp * delhat           difth  = diffus(i,kn,3) + diftp * delhat
1174    
 #ifdef KPP_TEST_DENOM  
 ctl replace (Important!!! not phepsi**4 !!!)  
1175           f1 = stable(i) * conc1 * bfsfc(i) /           f1 = stable(i) * conc1 * bfsfc(i) /
1176       &        max(ustar(i)**4,phepsi)       &        max(ustar(i)**4,phepsi)
          wm(i) = sign(eins,wm(i))*max(phepsi,abs(wm(i)))  
1177           gat1m(i) = visch / hbl(i) / wm(i)           gat1m(i) = visch / hbl(i) / wm(i)
1178           dat1m(i) = -viscp / wm(i) + f1 * visch           dat1m(i) = -viscp / wm(i) + f1 * visch
 #else  
          f1 = stable(i) * conc1 * bfsfc(i) / (ustar(i)**4+epsln)  
          gat1m(i) = visch / hbl(i) / (wm(i)+epsln)  
          dat1m(i) = -viscp / (wm(i)+epsln) + f1 * visch  
 #endif  
          dat1m(i) = min(dat1m(i),p0)  
1179    
 #ifdef KPP_TEST_DENOM  
          ws(i) = sign(eins,ws(i))*max(phepsi,abs(ws(i)))  
1180           gat1s(i) = difsh  / hbl(i) / ws(i)           gat1s(i) = difsh  / hbl(i) / ws(i)
1181           dat1s(i) = -difsp / ws(i) + f1 * difsh           dat1s(i) = -difsp / ws(i) + f1 * difsh
 #else  
          gat1s(i) = difsh / hbl(i) / (ws(i)+epsln)  
          dat1s(i) = -difsp / (ws(i)+epsln) + f1 * difsh  
 #endif  
          dat1s(i) = min(dat1s(i),p0)  
1182    
 #ifdef KPP_TEST_DENOM  
1183           gat1t(i) = difth /  hbl(i) / ws(i)           gat1t(i) = difth /  hbl(i) / ws(i)
1184           dat1t(i) = -diftp / ws(i) + f1 * difth           dat1t(i) = -diftp / ws(i) + f1 * difth
 #else  
          gat1t(i) = difth / hbl(i) / (ws(i)+epsln)  
          dat1t(i) = -diftp / (ws(i)+epsln) + f1 * difth  
 #endif  
          dat1t(i) = min(dat1t(i),p0)  
1185    
1186        end do        end do
1187    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1188    CADJ STORE gat1m = comlev1_kpp, key = ikppkey
1189    CADJ STORE gat1s = comlev1_kpp, key = ikppkey
1190    CADJ STORE gat1t = comlev1_kpp, key = ikppkey
1191    CADJ STORE dat1m = comlev1_kpp, key = ikppkey
1192    CADJ STORE dat1s = comlev1_kpp, key = ikppkey
1193    CADJ STORE dat1t = comlev1_kpp, key = ikppkey
1194    #endif
1195          do i = 1, imt
1196             dat1m(i) = min(dat1m(i),p0)
1197             dat1s(i) = min(dat1s(i),p0)
1198             dat1t(i) = min(dat1t(i),p0)
1199          end do
1200    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1201    CADJ STORE dat1m = comlev1_kpp, key = ikppkey
1202    CADJ STORE dat1s = comlev1_kpp, key = ikppkey
1203    CADJ STORE dat1t = comlev1_kpp, key = ikppkey
1204    #endif
1205    
1206        do ki = 1, Nr        do ki = 1, Nr
1207    
1208    #ifdef ALLOW_AUTODIFF_TAMC
1209             kkppkey = (ikppkey-1)*Nr + ki
1210    #endif
1211    
1212  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1213  c     compute turbulent velocity scales on the interfaces  c     compute turbulent velocity scales on the interfaces
1214  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 1185  c--------------------------------------- Line 1217  c---------------------------------------
1217              sig      = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i)              sig      = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i)
1218              sigma(i) = stable(i)*sig + (1.-stable(i))*min(sig,epsilon)              sigma(i) = stable(i)*sig + (1.-stable(i))*min(sig,epsilon)
1219           end do           end do
1220    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1221    CADJ STORE wm = comlev1_kpp_k, key = kkppkey
1222    CADJ STORE ws = comlev1_kpp_k, key = kkppkey
1223    #endif
1224    CADJ STORE sigma = comlev1_kpp_k, key = kkppkey
1225           call wscale (           call wscale (
1226       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1227       O        wm, ws )       O        wm, ws )
1228    CADJ STORE wm = comlev1_kpp_k, key = kkppkey
1229    CADJ STORE ws = comlev1_kpp_k, key = kkppkey
1230    
1231  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1232  c     compute the dimensionless shape functions at the interfaces  c     compute the dimensionless shape functions at the interfaces
# Line 1214  c--------------------------------------- Line 1253  c---------------------------------------
1253  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1254  c     nonlocal transport term = ghat * <ws>o  c     nonlocal transport term = ghat * <ws>o
1255  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1256  #ifdef KPP_TEST_DENOM  
1257              dum = ws(i) * hbl(i)              tempVar = ws(i) * hbl(i)
1258  ctl replace              ghat(i,ki) = (1.-stable(i)) * cg / max(phepsi,tempVar)
             ghat(i,ki) = (1.-stable(i)) * cg / max(phepsi,dum)  
 ctl end-replace  
 #else  
             dum = ws(i) * hbl(i) + epsln  
             ghat(i,ki) = (1. - stable(i)) * cg / dum  
 #endif  
1259    
1260           end do           end do
1261        end do        end do
# Line 1237  c--------------------------------------- Line 1270  c---------------------------------------
1270       &            + (1. - stable(i)) * min(sig,epsilon)       &            + (1. - stable(i)) * min(sig,epsilon)
1271        end do        end do
1272    
1273    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1274    CADJ STORE wm = comlev1_kpp, key = ikppkey
1275    CADJ STORE ws = comlev1_kpp, key = ikppkey
1276    #endif
1277    CADJ STORE sigma = comlev1_kpp, key = ikppkey
1278        call wscale (        call wscale (
1279       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1280       O        wm, ws )       O        wm, ws )
1281    CADJ STORE wm = comlev1_kpp, key = ikppkey
1282    CADJ STORE ws = comlev1_kpp, key = ikppkey
1283    
1284        do i = 1, imt        do i = 1, imt
1285           sig = -zgrid(kbl(i)-1) / hbl(i)           sig = -zgrid(kbl(i)-1) / hbl(i)
# Line 1280  c     hbl(imt)                  boundary Line 1320  c     hbl(imt)                  boundary
1320  c     kbl(imt)                  grid above hbl  c     kbl(imt)                  grid above hbl
1321  c     diffus(imt,0:Nrp1,mdiff) vertical diffusivities           (m^2/s)  c     diffus(imt,0:Nrp1,mdiff) vertical diffusivities           (m^2/s)
1322  c     casea(imt)                = 1 in caseA, = 0 in case B  c     casea(imt)                = 1 in caseA, = 0 in case B
1323        _RS     dkm1  (imt,mdiff)        _KPP_RL dkm1  (imt,mdiff)
1324        _RS     hbl   (imt)        _KPP_RL hbl   (imt)
1325        integer kbl   (imt)        integer kbl   (imt)
1326        _RS     diffus(imt,0:Nrp1,mdiff)        _KPP_RL diffus(imt,0:Nrp1,mdiff)
1327        _RS     casea (imt)        _KPP_RL casea (imt)
1328    
1329  c input/output  c input/output
1330  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)
1331        _RS     ghat (imt,Nr)        _KPP_RL ghat (imt,Nr)
1332    
1333  c output  c output
1334  c     enhanced bound. layer mixing coeff.  c     enhanced bound. layer mixing coeff.
1335        _RS     blmc  (imt,Nr,mdiff)        _KPP_RL blmc  (imt,Nr,mdiff)
1336    
1337  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1338    
1339  c local  c local
1340  c     fraction hbl lies beteen zgrid neighbors  c     fraction hbl lies beteen zgrid neighbors
1341        _RS delta        _KPP_RL delta
1342        integer ki, i, md        integer ki, i, md
1343        _RS dkmp5, dstar        _KPP_RL dkmp5, dstar
1344    
1345        do i = 1, imt        do i = 1, imt
1346           ki = kbl(i)-1           ki = kbl(i)-1
# Line 1326  c     fraction hbl lies beteen zgrid nei Line 1366  c     fraction hbl lies beteen zgrid nei
1366  c*************************************************************************  c*************************************************************************
1367    
1368        SUBROUTINE STATEKPP (        SUBROUTINE STATEKPP (
1369       I     bi, bj, myThid,       I     ikppkey, bi, bj, myThid,
1370       O     RHO1, DBLOC, DBSFC, TTALPHA, SSBETA)       O     RHO1, DBLOC, DBSFC, TTALPHA, SSBETA)
1371  c  c
1372  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 1360  c--------------------------------------- Line 1400  c---------------------------------------
1400  #include "EEPARAMS.h"  #include "EEPARAMS.h"
1401  #include "PARAMS.h"  #include "PARAMS.h"
1402  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
1403    #include "DYNVARS.h"
1404    
1405  c-------------- Routine arguments -----------------------------------------  c-------------- Routine arguments -----------------------------------------
1406        INTEGER bi, bj, myThid        INTEGER bi, bj, myThid
1407  #ifdef FRUGAL_KPP        _KPP_RL RHO1   ( ibot:itop, jbot:jtop       )
1408        _RS RHO1   (sNx,sNy)        _KPP_RL DBLOC  ( ibot:itop, jbot:jtop, Nr   )
1409        _RS DBLOC  (sNx,sNy,Nr)        _KPP_RL DBSFC  ( ibot:itop, jbot:jtop, Nr   )
1410        _RS DBSFC  (sNx,sNy,Nr)        _KPP_RL TTALPHA( ibot:itop, jbot:jtop, Nrp1 )
1411        _RS TTALPHA(sNx,sNy,Nrp1)        _KPP_RL SSBETA ( ibot:itop, jbot:jtop, Nrp1 )
       _RS SSBETA (sNx,sNy,Nrp1)  
 #else  
       _RS RHO1   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS DBLOC  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RS DBSFC  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RS TTALPHA(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nrp1)  
       _RS SSBETA (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nrp1)  
 #endif  
1412    
1413  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1414    
# Line 1395  c     work1, work2 - work arrays for hol Line 1428  c     work1, work2 - work arrays for hol
1428        _RL WORK2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1429        _RL WORK3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1430        INTEGER I, J, K        INTEGER I, J, K
1431          INTEGER ikppkey, kkppkey
1432    
1433  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
1434    
1435          kkppkey = (ikppkey-1)*Nr + 1
1436    
1437    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1438    CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1439    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1440    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1441        call FIND_RHO(        call FIND_RHO(
1442  #ifdef FRUGAL_KPP       I     bi, bj, ibot, itop, jbot, jtop, 1, 1,
1443       I     bi, bj, 1, sNx, 1, sNy,  1, 1, eosType,       I     theta, salt,
 #else  
      I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy,  1, 1, eosType,  
 #endif  
1444       O     WORK1,       O     WORK1,
1445       I     myThid )       I     myThid )
1446    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1447    CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1448    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1449    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1450    
1451        call FIND_ALPHA(        call FIND_ALPHA(
1452  #ifdef FRUGAL_KPP       I     bi, bj, ibot, itop, jbot, jtop, 1, 1,
      I     bi, bj, 1, sNx, 1, sNy,  1, 1, eosType,  
 #else  
      I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy,  1, 1, eosType,  
 #endif  
1453       O     WORK2 )       O     WORK2 )
1454    
1455        call FIND_BETA(        call FIND_BETA(
1456  #ifdef FRUGAL_KPP       I     bi, bj, ibot, itop, jbot, jtop, 1, 1,
      I     bi, bj, 1, sNx, 1, sNy,  1, 1, eosType,  
 #else  
      I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy,  1, 1, eosType,  
 #endif  
1457       O     WORK3 )       O     WORK3 )
1458    
1459  #ifdef FRUGAL_KPP        DO J = jbot, jtop
1460        DO J = 1, sNy           DO I = ibot, itop
1461           DO I = 1, sNx              RHO1(I,J)      = WORK1(I,J) + rhoConst
 #else  
       DO J = 1-OLy, sNy+OLy  
          DO I = 1-OLx, sNx+OLx  
 #endif  
             RHO1(I,J)      = WORK1(I,J) + rhonil  
1462              TTALPHA(I,J,1) = WORK2(I,J)              TTALPHA(I,J,1) = WORK2(I,J)
1463              SSBETA(I,J,1)  = WORK3(I,J)              SSBETA(I,J,1)  = WORK3(I,J)
1464              DBSFC(I,J,1)   = 0.              DBSFC(I,J,1)   = 0.
# Line 1442  c calculate alpha, beta, and gradients i Line 1470  c calculate alpha, beta, and gradients i
1470  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)
1471        DO K = 2, Nr        DO K = 2, Nr
1472    
1473          kkppkey = (ikppkey-1)*Nr + k
1474    
1475    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1476    CADJ STORE theta(:,:,k,bi,bj) = comlev1_kpp_k, key=kkppkey
1477    CADJ STORE salt (:,:,k,bi,bj) = comlev1_kpp_k, key=kkppkey
1478    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1479           call FIND_RHO(           call FIND_RHO(
1480  #ifdef FRUGAL_KPP       I        bi, bj, ibot, itop, jbot, jtop, K, K,
1481       I        bi, bj, 1, sNx, 1, sNy,  K,   K, eosType,       I        theta, salt,
 #else  
      I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy,  K,   K, eosType,  
 #endif  
1482       O        RHOK,       O        RHOK,
1483       I        myThid )       I        myThid )
1484    
1485    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1486    CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_kpp_k, key=kkppkey
1487    CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_kpp_k, key=kkppkey
1488    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1489           call FIND_RHO(           call FIND_RHO(
1490  #ifdef FRUGAL_KPP       I        bi, bj, ibot, itop, jbot, jtop, K-1, K,
1491       I        bi, bj, 1, sNx, 1, sNy,  K-1, K, eosType,       I        theta, salt,
 #else  
      I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy,  K-1, K, eosType,  
 #endif  
1492       O        RHOKM1,       O        RHOKM1,
1493       I        myThid )       I        myThid )
1494    
1495    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1496    CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1497    CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1498    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1499           call FIND_RHO(           call FIND_RHO(
1500  #ifdef FRUGAL_KPP       I        bi, bj, ibot, itop, jbot, jtop, 1, K,
1501       I        bi, bj, 1, sNx, 1, sNy,  1,   K, eosType,       I        theta, salt,
 #else  
      I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy,  1,   K, eosType,  
 #endif  
1502       O        RHO1K,       O        RHO1K,
1503       I        myThid )       I        myThid )
1504    
1505    #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1506    CADJ STORE rhok  (:,:)          = comlev1_kpp_k, key=kkppkey
1507    CADJ STORE rhokm1(:,:)          = comlev1_kpp_k, key=kkppkey
1508    CADJ STORE rho1k (:,:)          = comlev1_kpp_k, key=kkppkey
1509    #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1510    
1511           call FIND_ALPHA(           call FIND_ALPHA(
1512  #ifdef FRUGAL_KPP       I        bi, bj, ibot, itop, jbot, jtop, K, K,
      I        bi, bj, 1, sNx, 1, sNy,  K,   K, eosType,  
 #else  
      I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy,  K,   K, eosType,  
 #endif  
1513       O        WORK1 )       O        WORK1 )
1514    
1515           call FIND_BETA(           call FIND_BETA(
1516  #ifdef FRUGAL_KPP       I        bi, bj, ibot, itop, jbot, jtop, K, K,
      I        bi, bj, 1, sNx, 1, sNy,  K,   K, eosType,  
 #else  
      I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy,  K,   K, eosType,  
 #endif  
1517       O        WORK2 )       O        WORK2 )
1518    
1519  #ifdef FRUGAL_KPP           DO J = jbot, jtop
1520           DO J = 1, sNy              DO I = ibot, itop
             DO I = 1, sNx  
 #else  
          DO J = 1-OLy, sNy+OLy  
             DO I = 1-OLx, sNx+OLx  
 #endif  
1521                 TTALPHA(I,J,K) = WORK1 (I,J)                 TTALPHA(I,J,K) = WORK1 (I,J)
1522                 SSBETA(I,J,K)  = WORK2 (I,J)                 SSBETA(I,J,K)  = WORK2 (I,J)
1523                 DBLOC(I,J,K-1) = gravity * (RHOK(I,J) - RHOKM1(I,J)) /                 DBLOC(I,J,K-1) = gravity * (RHOK(I,J) - RHOKM1(I,J)) /
1524       &                                    (RHOK(I,J) + rhonil)       &                                    (RHOK(I,J) + rhoConst)
1525                 DBSFC(I,J,K)   = gravity * (RHOK(I,J) - RHO1K (I,J)) /                 DBSFC(I,J,K)   = gravity * (RHOK(I,J) - RHO1K (I,J)) /
1526       &                                    (RHOK(I,J) + rhonil)       &                                    (RHOK(I,J) + rhoConst)
1527              END DO              END DO
1528           END DO           END DO
1529    
1530        END DO        END DO
1531    
1532  c     compute arrays for K = Nrp1  c     compute arrays for K = Nrp1
1533  #ifdef FRUGAL_KPP        DO J = jbot, jtop
1534        DO J = 1, sNy           DO I = ibot, itop
          DO I = 1, sNx  
 #else  
       DO J = 1-OLy, sNy+OLy  
          DO I = 1-OLx, sNx+OLx  
 #endif  
1535              TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)              TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)
1536              SSBETA(I,J,Nrp1)  = SSBETA(I,J,Nr)              SSBETA(I,J,Nrp1)  = SSBETA(I,J,Nr)
1537              DBLOC(I,J,Nr)     = 0.              DBLOC(I,J,Nr)     = 0.
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.21
  ViewVC Help
Powered by ViewVC 1.1.22