/[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.2 by heimbach, Tue Sep 12 18:14:32 2000 UTC revision 1.9 by heimbach, Tue Sep 18 20:30:59 2001 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     , bo, bosol, dbloc, Ritop, coriol
26       I     , ikey       I     , ikey
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 ikey
86    
# 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    CADJ STORE ghat = comlev1_kpp, key = ikey
129    
130        call Ri_iwmix (        call Ri_iwmix (
131       I       kmtj, shsq, dbloc, ghat       I       kmtj, shsq, dbloc, ghat
132       I     , ikey       I     , ikey
# Line 146  c diagnose the new boundary layer depth Line 152  c diagnose the new boundary layer depth
152  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
153    
154        call bldepth (        call bldepth (
155       I       kmtj       I       mytime, mythid
156         I     , kmtj
157       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol
158       I     , ikey       I     , ikey
159       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma
# Line 160  c--------------------------------------- Line 167  c---------------------------------------
167    
168        call blmix (        call blmix (
169       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl
170       O     , dkm1, blmc, ghat, sigma       O     , dkm1, blmc, ghat, sigma, ikey
171       &     )       &     )
172    
173  CADJ STORE dkm1,blmc,ghat = comlev1_kpp, key = ikey  CADJ STORE dkm1,blmc,ghat = comlev1_kpp, key = ikey
# Line 198  c--------------------------------------- Line 205  c---------------------------------------
205  c*************************************************************************  c*************************************************************************
206    
207        subroutine bldepth (        subroutine bldepth (
208       I       kmtj       I       mytime, mythid
209         I     , kmtj
210       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol
211       I     , ikey       I     , ikey
212       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma
# Line 235  c Line 243  c
243    
244  c input  c input
245  c------  c------
246    c myTime    : current time in simulation
247    c myThid    : thread number for this instance of the routine
248  c kmtj      : number of vertical layers  c kmtj      : number of vertical layers
249  c dvsq      : (velocity shear re sfc)^2             ((m/s)^2)  c dvsq      : (velocity shear re sfc)^2             ((m/s)^2)
250  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 255  c ustar     : surface friction velocity
255  c bo        : surface turbulent buoyancy forcing    (m^2/s^3)  c bo        : surface turbulent buoyancy forcing    (m^2/s^3)
256  c bosol     : radiative buoyancy forcing            (m^2/s^3)  c bosol     : radiative buoyancy forcing            (m^2/s^3)
257  c coriol    : Coriolis parameter                        (1/s)  c coriol    : Coriolis parameter                        (1/s)
258          _RL     mytime
259          integer mythid
260        integer kmtj(imt)        integer kmtj(imt)
261        _RS dvsq   (imt,Nr)        _KPP_RL dvsq  (imt,Nr)
262        _RS dbloc  (imt,Nr)        _KPP_RL dbloc (imt,Nr)
263        _RS Ritop  (imt,Nr)        _KPP_RL Ritop (imt,Nr)
264        _RS ustar  (imt)        _KPP_RL ustar (imt)
265        _RS bo     (imt)        _KPP_RL bo    (imt)
266        _RS bosol  (imt)        _KPP_RL bosol (imt)
267        _RS coriol (imt)        _KPP_RL coriol(imt)
268        integer ikey        integer ikey
269    
270  c  output  c  output
# Line 264  c casea     : =1 in case A, =0 in case B Line 276  c casea     : =1 in case A, =0 in case B
276  c kbl       : -1 of first grid level below hbl  c kbl       : -1 of first grid level below hbl
277  c Rib       : Bulk Richardson number  c Rib       : Bulk Richardson number
278  c sigma     : normalized depth (d/hbl)  c sigma     : normalized depth (d/hbl)
279        _RS hbl    (imt)        _KPP_RL hbl   (imt)
280        _RS bfsfc  (imt)        _KPP_RL bfsfc (imt)
281        _RS stable (imt)        _KPP_RL stable(imt)
282        _RS casea  (imt)        _KPP_RL casea (imt)
283        integer kbl (imt)        integer kbl   (imt)
284        _RS Rib    (imt,Nr)        _KPP_RL Rib   (imt,Nr)
285        _RS sigma  (imt)        _KPP_RL sigma (imt)
286    
287  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
288    
289  c  local  c  local
290  c-------  c-------
291  c wm, ws    : turbulent velocity scales         (m/s)  c wm, ws    : turbulent velocity scales         (m/s)
292  c zwork     : depth work array        _KPP_RL wm(imt), ws(imt)
293        _RS wm(imt), ws(imt)        _RL     worka(imt)
       _RS zwork(imt)  
294    
295        _RS bvsq, vtsq, hekman, hmonob, hlimit        _KPP_RL bvsq, vtsq, hekman, hmonob, hlimit, tempVar1, tempVar2
296        integer i, kl        integer i, kl
297    
298        _RS        p5    , eins    , m1        _KPP_RL     p5    , eins
299        parameter (p5=0.5, eins=1.0, m1=-1.0 )        parameter ( p5=0.5, eins=1.0 )
300          _RL         minusone
301  crg intermediate result in RL to avoid overflow in adjoint        parameter ( minusone=-1.0 )
        _RL        dpshear2  
 #ifdef KPP_TEST_DENOM  
        _RL        dhelp  
 #endif  
302    
303  c find bulk Richardson number at every grid level until > Ricr  c find bulk Richardson number at every grid level until > Ricr
304  c  c
# Line 315  c     initialize hbl and kbl to bottomed Line 322  c     initialize hbl and kbl to bottomed
322  c     compute bfsfc = sw fraction at hbf * zgrid  c     compute bfsfc = sw fraction at hbf * zgrid
323    
324           do i = 1, imt           do i = 1, imt
325              zwork(i) = zgrid(kl)              worka(i) = zgrid(kl)
326           end do           end do
327           call SWFRAC(           call SWFRAC(
328       I        imt, hbf, zwork,       I        imt, hbf,
329       O        bfsfc)       I        mytime, mythid,
330         U        worka )
331    
332           do i = 1, imt           do i = 1, imt
333    
# Line 329  c     use caseA as temporary array Line 337  c     use caseA as temporary array
337    
338  c     compute bfsfc= Bo + radiative contribution down to hbf * hbl  c     compute bfsfc= Bo + radiative contribution down to hbf * hbl
339    
340              bfsfc(i) = bo(i) + bosol(i)*(1. - bfsfc(i))              bfsfc(i) = bo(i) + bosol(i)*(1. - worka(i))
341              stable(i) = p5 + sign(p5,bfsfc(i))              stable(i) = p5 + sign(p5,bfsfc(i))
342              sigma(i) = stable(i) + (1. - stable(i)) * epsilon              sigma(i) = stable(i) + (1. - stable(i)) * epsilon
343    
# Line 370  c Line 378  c
378  c     rg: assignment to double precision variable to avoid overflow  c     rg: assignment to double precision variable to avoid overflow
379  c     ph: test for zero nominator  c     ph: test for zero nominator
380  c  c
381  #ifdef KPP_TEST_DENOM  
382  ctl replace              tempVar1  = dvsq(i,kl) + vtsq          
383              dpshear2  = dvsq(i,kl) + vtsq                        tempVar2 = max(tempVar1, phepsi)
384              dpshear2 = max(dpshear2, phepsi)              Rib(i,kl) = Ritop(i,kl) / tempVar2
385              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  
386           end do           end do
387        end do        end do
388                
# Line 394  CADJ &   , key = ikey, shape = (/ (sNx+2 Line 397  CADJ &   , key = ikey, shape = (/ (sNx+2
397    
398        do i = 1, imt        do i = 1, imt
399           kl = kbl(i)           kl = kbl(i)
   
400  c     linearly interpolate to find hbl where Rib = Ricr  c     linearly interpolate to find hbl where Rib = Ricr
   
401           if (kl.gt.1 .and. kl.lt.kmtj(i)) then           if (kl.gt.1 .and. kl.lt.kmtj(i)) then
402  #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  
403              hbl(i) = -zgrid(kl-1) + (zgrid(kl-1)-zgrid(kl)) *              hbl(i) = -zgrid(kl-1) + (zgrid(kl-1)-zgrid(kl)) *
404       1           (Ricr - Rib(i,kl-1)) / (Rib(i,kl)-Rib(i,kl-1))       1           (Ricr - Rib(i,kl-1)) / tempVar1
 #endif  
405           endif           endif
406        end do        end do
407    
# Line 416  c--------------------------------------- Line 412  c---------------------------------------
412  c     find stability and buoyancy forcing for boundary layer  c     find stability and buoyancy forcing for boundary layer
413  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
414    
415          do i = 1, imt
416             worka(i) = hbl(i)
417          end do
418        call SWFRAC(        call SWFRAC(
419       I     imt, m1, hbl,       I     imt, minusone,
420       O     bfsfc)       I     mytime, mythid,
421         U     worka )
422    
423          do i = 1, imt
424             bfsfc(i)  = bo(i) + bosol(i) * (1. - worka(i))
425          end do
426  CADJ store bfsfc = comlev1_kpp  CADJ store bfsfc = comlev1_kpp
427  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
428    
429    c--   ensure bfsfc is never 0
430        do i = 1, imt        do i = 1, imt
          bfsfc(i)  = bo(i) + bosol(i) * (1. - bfsfc(i))  
431           stable(i) = p5 + sign( p5, bfsfc(i) )           stable(i) = p5 + sign( p5, bfsfc(i) )
432             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  
   
433        end do        end do
434    
435    CADJ store bfsfc = comlev1_kpp
436    CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
437    
438  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
439  c check hbl limits for hekman or hmonob  c check hbl limits for hekman or hmonob
440  c     ph: test for zero nominator  c     ph: test for zero nominator
# Line 445  c--------------------------------------- Line 442  c---------------------------------------
442    
443        do i = 1, imt        do i = 1, imt
444           if (bfsfc(i) .gt. 0.0) then           if (bfsfc(i) .gt. 0.0) then
 #ifdef KPP_TEST_DENOM  
 ctl replace  
445              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  
446              hmonob = cmonob * ustar(i)*ustar(i)*ustar(i)              hmonob = cmonob * ustar(i)*ustar(i)*ustar(i)
447       &           / vonk / bfsfc(i)       &           / vonk / bfsfc(i)
448              hlimit = stable(i) * min(hekman,hmonob)              hlimit = stable(i) * min(hekman,hmonob)
449       &             + (stable(i)-1.) * zgrid(Nr)       &             + (stable(i)-1.) * zgrid(Nr)
450              hbl(i) = min(hbl(i),hlimit)              hbl(i) = min(hbl(i),hlimit)
451              hbl(i) = max(hbl(i),minKPPhbl)           end if
452           endif        end do
453    CADJ store hbl = comlev1_kpp
454    CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
455    
456          do i = 1, imt
457             hbl(i) = max(hbl(i),minKPPhbl)
458           kbl(i) = kmtj(i)           kbl(i) = kmtj(i)
459        end do        end do
460    
# Line 481  c--------------------------------------- Line 477  c---------------------------------------
477  c      find stability and buoyancy forcing for final hbl values  c      find stability and buoyancy forcing for final hbl values
478  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
479    
480          do i = 1, imt
481             worka(i) = hbl(i)
482          end do
483        call SWFRAC(        call SWFRAC(
484       I     imt, m1, hbl,       I     imt, minusone,
485       O     bfsfc)       I     mytime, mythid,
486         U     worka )
487    
488          do i = 1, imt
489             bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i))
490          end do
491  CADJ store bfsfc = comlev1_kpp  CADJ store bfsfc = comlev1_kpp
492  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
493    
494    c--   ensures bfsfc is never 0
495        do i = 1, imt        do i = 1, imt
          bfsfc(i) = bo(i) + bosol(i) * (1. - bfsfc(i))  
496           stable(i) = p5 + sign( p5, bfsfc(i) )           stable(i) = p5 + sign( p5, bfsfc(i) )
497  #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  
498        end do        end do
499    
500  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 540  c sigma   : normalized depth (d/hbl) Line 536  c sigma   : normalized depth (d/hbl)
536  c hbl     : boundary layer depth (m)  c hbl     : boundary layer depth (m)
537  c ustar   : surface friction velocity         (m/s)  c ustar   : surface friction velocity         (m/s)
538  c bfsfc   : total surface buoyancy flux       (m^2/s^3)  c bfsfc   : total surface buoyancy flux       (m^2/s^3)
539        _RS sigma(imt)        _KPP_RL sigma(imt)
540        _RS hbl  (imt)        _KPP_RL hbl  (imt)
541        _RS ustar(imt)        _KPP_RL ustar(imt)
542        _RS bfsfc(imt)        _KPP_RL bfsfc(imt)
543    
544  c  output  c  output
545  c--------  c--------
546  c wm, ws  : turbulent velocity scales at sigma  c wm, ws  : turbulent velocity scales at sigma
547        _RS wm(imt), ws(imt)        _KPP_RL wm(imt), ws(imt)
548    
549  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
550    
551  c local  c local
552  c------  c------
553  c zehat   : = zeta *  ustar**3  c zehat   : = zeta *  ustar**3
554        _RS zehat        _KPP_RL zehat
555    
556        integer iz, izp1, ju, i, jup1        integer iz, izp1, ju, i, jup1
557        _RS udiff, zdiff, zfrac, ufrac, fzfrac, wam, wbm, was, wbs, u3        _KPP_RL udiff, zdiff, zfrac, ufrac, fzfrac, wam
558        _RL dum        _KPP_RL wbm, was, wbs, u3, tempVar
559    
560  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
561  c use lookup table for zehat < zmax only; otherwise use  c use lookup table for zehat < zmax only; otherwise use
# Line 576  c--------------------------------------- Line 572  c---------------------------------------
572              iz    = min( iz, nni )              iz    = min( iz, nni )
573              iz    = max( iz, 0 )              iz    = max( iz, 0 )
574              izp1  = iz + 1              izp1  = iz + 1
575    
576              udiff = ustar(i) - umin              udiff = ustar(i) - umin
577              ju    = int( udiff / deltau )              ju    = int( udiff / deltau )
578              ju    = min( ju, nnj )              ju    = min( ju, nnj )
579              ju    = max( ju, 0 )              ju    = max( ju, 0 )
580              jup1  = ju + 1              jup1  = ju + 1
581    
582              zfrac = zdiff / deltaz - float(iz)              zfrac = zdiff / deltaz - float(iz)
583              ufrac = udiff / deltau - float(ju)              ufrac = udiff / deltau - float(ju)
584    
585              fzfrac= 1. - zfrac              fzfrac= 1. - zfrac
586              wam   = fzfrac     * wmt(iz,jup1) + zfrac * wmt(izp1,jup1)              wam   = fzfrac     * wmt(iz,jup1) + zfrac * wmt(izp1,jup1)
587              wbm   = fzfrac     * wmt(iz,ju  ) + zfrac * wmt(izp1,ju  )              wbm   = fzfrac     * wmt(iz,ju  ) + zfrac * wmt(izp1,ju  )
588              wm(i) = (1.-ufrac) * wbm          + ufrac * wam              wm(i) = (1.-ufrac) * wbm          + ufrac * wam
589    
590              was   = fzfrac     * wst(iz,jup1) + zfrac * wst(izp1,jup1)              was   = fzfrac     * wst(iz,jup1) + zfrac * wst(izp1,jup1)
591              wbs   = fzfrac     * wst(iz,ju  ) + zfrac * wst(izp1,ju  )              wbs   = fzfrac     * wst(iz,ju  ) + zfrac * wst(izp1,ju  )
592              ws(i) = (1.-ufrac) * wbs          + ufrac * was              ws(i) = (1.-ufrac) * wbs          + ufrac * was
593    
594           else           else
595    
596              u3    = ustar(i) * ustar(i) * ustar(i)              u3 = ustar(i) * ustar(i) * ustar(i)
597              dum   = u3 + conc1 * zehat              tempVar = u3 + conc1 * zehat
598              wm(i) = vonk     * ustar(i) * u3 / dum              wm(i) = vonk * ustar(i) * u3 / tempVar
599              ws(i) = wm(i)              ws(i) = wm(i)
600    
601           endif           endif
602    
603        end do        end do
# Line 636  c     shsq   (imt,Nr)      (local veloci Line 632  c     shsq   (imt,Nr)      (local veloci
632  c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)  c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)
633  c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)  c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)
634        integer kmtj   (imt)        integer kmtj   (imt)
635        _RS     shsq   (imt,Nr)        _KPP_RL shsq   (imt,Nr)
636        _RS     dbloc  (imt,Nr)        _KPP_RL dbloc  (imt,Nr)
637        _RS     dblocSm(imt,Nr)        _KPP_RL dblocSm(imt,Nr)
638        integer ikey        integer ikey
639    
640  c output  c output
641  c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)  c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)
642  c     diffus(imt,0:Nrp1,2)  vertical scalar diffusivity             (m^2/s)  c     diffus(imt,0:Nrp1,2)  vertical scalar diffusivity             (m^2/s)
643  c     diffus(imt,0:Nrp1,3)  vertical temperature diffusivity        (m^2/s)  c     diffus(imt,0:Nrp1,3)  vertical temperature diffusivity        (m^2/s)
644        _RS diffus(imt,0:Nrp1,3)        _KPP_RL diffus(imt,0:Nrp1,3)
645    
646  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
647    
648  c local variables  c local variables
649  c     Rig                   local Richardson number  c     Rig                   local Richardson number
650  c     fRi, fcon             function of Rig  c     fRi, fcon             function of Rig
651        _RS Rig        _KPP_RL Rig
652        _RS fRi, fcon        _KPP_RL fRi, fcon
653        _RS ratio        _KPP_RL ratio
654        integer i, ki, mr        integer i, ki, mr
655        _RS c1, c0        _KPP_RL c1, c0
656    
657    #ifdef ALLOW_KPP_VERTICALLY_SMOOTH
658    CADJ INIT kpp_ri_tape_mr = common, 1
659    #endif
660    
661  c constants  c constants
662        c1 = 1.0        c1 = 1.0
# Line 667  c     compute interior gradient Ri at al Line 667  c     compute interior gradient Ri at al
667  c     use diffus(*,*,1) as temporary storage of Ri to be smoothed  c     use diffus(*,*,1) as temporary storage of Ri to be smoothed
668  c     use diffus(*,*,2) as temporary storage for Brunt-Vaisala squared  c     use diffus(*,*,2) as temporary storage for Brunt-Vaisala squared
669  c     set values at bottom and below to nearest value above bottom  c     set values at bottom and below to nearest value above bottom
670    #ifdef ALLOW_AUTODIFF_TAMC
671    C     break data flow dependence on diffus
672          diffus(1,1,1) = 0.0
673    #endif
674    
675        do ki = 1, Nr        do ki = 1, Nr
676           do i = 1, imt           do i = 1, imt
677              if     (kmtj(i) .EQ. 0      ) then              if     (kmtj(i) .LE. 1      ) then
678                 diffus(i,ki,1) = 0.                 diffus(i,ki,1) = 0.
679                 diffus(i,ki,2) = 0.                 diffus(i,ki,2) = 0.
680              elseif (ki      .GE. kmtj(i)) then              elseif (ki      .GE. kmtj(i)) then
# Line 678  c     set values at bottom and below to Line 682  c     set values at bottom and below to
682                 diffus(i,ki,2) = diffus(i,ki-1,2)                 diffus(i,ki,2) = diffus(i,ki-1,2)
683              else              else
684                 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  
685       &              / max( Shsq(i,ki), phepsi )       &              / max( Shsq(i,ki), phepsi )
 #else  
      &              / ( shsq(i,ki) + epsln )  
 #endif  
686                 diffus(i,ki,2) = dbloc(i,ki)   / (zgrid(ki)-zgrid(ki+1))                 diffus(i,ki,2) = dbloc(i,ki)   / (zgrid(ki)-zgrid(ki+1))
687              endif              endif
688           end do           end do
# Line 690  c     set values at bottom and below to Line 690  c     set values at bottom and below to
690    
691  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
692  c     vertically smooth Ri  c     vertically smooth Ri
693    #ifdef ALLOW_KPP_VERTICALLY_SMOOTH
694        do mr = 1, num_v_smooth_Ri        do mr = 1, num_v_smooth_Ri
695    
696  CADJ store diffus(:,:,1) = comlev1_kpp_sm  CADJ store diffus(:,:,1) = kpp_ri_tape_mr
697  CADJ &  , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2 /)  CADJ &  , key=mr, shape=(/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2 /)
698    
699           call z121 (           call z121 (
700       U     diffus(1,0,1))       U     diffus(1,0,1))
701        end do        end do
702    #endif
 CADJ store diffus = comlev1_kpp  
 CADJ &  , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2,3 /)  
703    
704  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
705  c                           after smoothing loop  c                           after smoothing loop
# Line 769  c     penetrative convection. Line 767  c     penetrative convection.
767  c input/output  c input/output
768  c-------------  c-------------
769  c v     : 2-D array to be smoothed in Nrp1 direction  c v     : 2-D array to be smoothed in Nrp1 direction
770        _RS v(imt,0:Nrp1)        _KPP_RL v(imt,0:Nrp1)
771    
772  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
773    
774  c local  c local
775        _RS zwork, zflag        _KPP_RL zwork, zflag
776        _RS KRi_range(1:Nrp1)        _KPP_RL KRi_range(1:Nrp1)
777        integer i, k, km1, kp1        integer i, k, km1, kp1
778    
779        _RS         p0      , p25       , p5      , p2        _KPP_RL     p0      , p25       , p5      , p2
780        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 )
781    
782        KRi_range(Nrp1) = p0        KRi_range(Nrp1) = p0
# Line 789  C--   dummy assignment to end declaratio Line 787  C--   dummy assignment to end declaratio
787    
788  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
789  CHPF$ INDEPENDENT  CHPF$ INDEPENDENT
790    CADJ INIT z121tape = common, Nr
791  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
792    
793        do i = 1, imt        do i = 1, imt
794    
795             k = 1
796    CADJ STORE v(i,k) = z121tape
797           v(i,Nrp1) = v(i,Nr)           v(i,Nrp1) = v(i,Nr)
798    
799           do k = 1, Nr           do k = 1, Nr
# Line 806  CHPF$ INDEPENDENT Line 808  CHPF$ INDEPENDENT
808           zflag  = p2 + KRi_range(1) * KRi_range(2)           zflag  = p2 + KRi_range(1) * KRi_range(2)
809           v(i,1) = v(i,1) / zflag           v(i,1) = v(i,1) / zflag
810    
 CADJ INIT z121tape = common, Nr  
811           do k = 2, Nr           do k = 2, Nr
812  CADJ STORE v(i,k), zwork = z121tape  CADJ STORE v(i,k), zwork = z121tape
813              km1 = k - 1              km1 = k - 1
# Line 829  CADJ STORE v(i,k), zwork = z121tape Line 830  CADJ STORE v(i,k), zwork = z121tape
830    
831  c*************************************************************************  c*************************************************************************
832    
833        subroutine smooth_horiz_rs (        subroutine kpp_smooth_horiz (
834       I     k, bi, bj,       I     k, bi, bj,
835       I     fld_in,       U     fld )
      O     fld_out )  
836    
837  c     Apply horizontal smoothing to RS 2-D array  c     Apply horizontal smoothing to KPP array
838    
839        IMPLICIT NONE        IMPLICIT NONE
840  #include "SIZE.h"  #include "SIZE.h"
841  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
842    
843  c     input  c     input
844  c     k, bi, bj : array indices  c     bi, bj : array indices
845  c     fld_in    : 2-D array to be smoothed  c     k      : vertical index used for masking
846        integer k, bi, bj        integer k, bi, bj
       _RS fld_in(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
847    
848  c     output  c     input/output
849  c     fld_out   : smoothed 2-D array  c     fld    : 2-D array to be smoothed
850        _RS fld_out(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _KPP_RL fld( ibot:itop, jbot:jtop )
851    
852  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
853    
854  c     local  c     local
855        integer i, j, im1, ip1, jm1, jp1        integer i, j, im1, ip1, jm1, jp1
856        _RS tempVar        _KPP_RL tempVar
857        _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 )  
858    
859        integer    jmin        , jmax        integer    imin       , imax       , jmin       , jmax
860        parameter( jmin=1-OLy+1, jmax=sNy+OLy-1 )        parameter( imin=ibot+1, imax=itop-1, jmin=jbot+1, jmax=jtop-1 )
861    
862        _RS        p0    , p5    , p25     , p125      , p0625        _KPP_RL    p0    , p5    , p25     , p125      , p0625
863        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 )
864    
865        DO j = jmin, jmax        DO j = jmin, jmax
# Line 884  c     local Line 880  c     local
880       &                     pMask(ip1,jp1,k,bi,bj) )       &                     pMask(ip1,jp1,k,bi,bj) )
881              IF ( tempVar .GE. p25 ) THEN              IF ( tempVar .GE. p25 ) THEN
882                 fld_tmp(i,j) = (                 fld_tmp(i,j) = (
883       &              p25  * fld_in(i  ,j  )*pMask(i  ,j  ,k,bi,bj) +       &              p25  * fld(i  ,j  )*pMask(i  ,j  ,k,bi,bj) +
884       &              p125 *(fld_in(im1,j  )*pMask(im1,j  ,k,bi,bj) +       &              p125 *(fld(im1,j  )*pMask(im1,j  ,k,bi,bj) +
885       &                     fld_in(ip1,j  )*pMask(ip1,j  ,k,bi,bj) +       &                     fld(ip1,j  )*pMask(ip1,j  ,k,bi,bj) +
886       &                     fld_in(i  ,jm1)*pMask(i  ,jm1,k,bi,bj) +       &                     fld(i  ,jm1)*pMask(i  ,jm1,k,bi,bj) +
887       &                     fld_in(i  ,jp1)*pMask(i  ,jp1,k,bi,bj))+       &                     fld(i  ,jp1)*pMask(i  ,jp1,k,bi,bj))+
888       &              p0625*(fld_in(im1,jm1)*pMask(im1,jm1,k,bi,bj) +       &              p0625*(fld(im1,jm1)*pMask(im1,jm1,k,bi,bj) +
889       &                     fld_in(im1,jp1)*pMask(im1,jp1,k,bi,bj) +       &                     fld(im1,jp1)*pMask(im1,jp1,k,bi,bj) +
890       &                     fld_in(ip1,jm1)*pMask(ip1,jm1,k,bi,bj) +       &                     fld(ip1,jm1)*pMask(ip1,jm1,k,bi,bj) +
891       &                     fld_in(ip1,jp1)*pMask(ip1,jp1,k,bi,bj)))       &                     fld(ip1,jp1)*pMask(ip1,jp1,k,bi,bj)))
892       &              / tempVar       &              / tempVar
893              ELSE              ELSE
894                 fld_tmp(i,j) = fld_in(i,j)                 fld_tmp(i,j) = fld(i,j)
895              ENDIF              ENDIF
896           ENDDO           ENDDO
897        ENDDO        ENDDO
# Line 903  c     local Line 899  c     local
899  c     transfer smoothed field to output array  c     transfer smoothed field to output array
900        DO j = jmin, jmax        DO j = jmin, jmax
901           DO i = imin, imax           DO i = imin, imax
902              fld_out(i,j) = fld_tmp(i,j)              fld(i,j) = fld_tmp(i,j)
903           ENDDO           ENDDO
904        ENDDO        ENDDO
905    
 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  
   
906  #endif /* ALLOW_KPP */  #endif /* ALLOW_KPP */
907    
908        return        return
# Line 932  c     set output array edges to input fi Line 910  c     set output array edges to input fi
910    
911  c*************************************************************************  c*************************************************************************
912    
913        subroutine smooth_horiz_rl (        subroutine smooth_horiz (
914       I     k, bi, bj,       I     k, bi, bj,
915       I     fld_in,       U     fld )
      O     fld_out )  
916    
917  c     Apply horizontal smoothing to RL 2-D array  c     Apply horizontal smoothing to global _RL 2-D array
918    
919        IMPLICIT NONE        IMPLICIT NONE
920  #include "SIZE.h"  #include "SIZE.h"
921  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
922    
923  c     input  c     input
924  c     k, bi, bj : array indices  c     bi, bj : array indices
925  c     fld_in    : 2-D array to be smoothed  c     k      : vertical index used for masking
926        integer k, bi, bj        integer k, bi, bj
       _RL fld_in(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
927    
928  c     output  c     input/output
929  c     fld_out   : smoothed 2-D array  c     fld    : 2-D array to be smoothed
930        _RL fld_out(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fld( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
931    
932  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
933    
934  c     local  c     local
935        integer i, j, im1, ip1, jm1, jp1        integer i, j, im1, ip1, jm1, jp1
936        _RL tempVar        _RL tempVar
937        _RL fld_tmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fld_tmp( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
938    
939        integer    imin        , imax        integer   imin      , imax          , jmin      , jmax
940        parameter( imin=1-OLx+1, imax=sNx+OLx-1 )        parameter(imin=2-OLx, imax=sNx+OLx-1, jmin=2-OLy, jmax=sNy+OLy-1)
941    
942        integer    jmin        , jmax        _RL        p0    , p5    , p25     , p125      , p0625
       parameter( jmin=1-OLy+1, jmax=sNy+OLy-1 )  
   
       _RS        p0    , p5    , p25     , p125      , p0625  
943        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 )
944    
945        DO j = jmin, jmax        DO j = jmin, jmax
# Line 987  c     local Line 960  c     local
960       &                     pMask(ip1,jp1,k,bi,bj) )       &                     pMask(ip1,jp1,k,bi,bj) )
961              IF ( tempVar .GE. p25 ) THEN              IF ( tempVar .GE. p25 ) THEN
962                 fld_tmp(i,j) = (                 fld_tmp(i,j) = (
963       &              p25  * fld_in(i  ,j  )*pMask(i  ,j  ,k,bi,bj) +       &              p25  * fld(i  ,j  )*pMask(i  ,j  ,k,bi,bj) +
964       &              p125 *(fld_in(im1,j  )*pMask(im1,j  ,k,bi,bj) +       &              p125 *(fld(im1,j  )*pMask(im1,j  ,k,bi,bj) +
965       &                     fld_in(ip1,j  )*pMask(ip1,j  ,k,bi,bj) +       &                     fld(ip1,j  )*pMask(ip1,j  ,k,bi,bj) +
966       &                     fld_in(i  ,jm1)*pMask(i  ,jm1,k,bi,bj) +       &                     fld(i  ,jm1)*pMask(i  ,jm1,k,bi,bj) +
967       &                     fld_in(i  ,jp1)*pMask(i  ,jp1,k,bi,bj))+       &                     fld(i  ,jp1)*pMask(i  ,jp1,k,bi,bj))+
968       &              p0625*(fld_in(im1,jm1)*pMask(im1,jm1,k,bi,bj) +       &              p0625*(fld(im1,jm1)*pMask(im1,jm1,k,bi,bj) +
969       &                     fld_in(im1,jp1)*pMask(im1,jp1,k,bi,bj) +       &                     fld(im1,jp1)*pMask(im1,jp1,k,bi,bj) +
970       &                     fld_in(ip1,jm1)*pMask(ip1,jm1,k,bi,bj) +       &                     fld(ip1,jm1)*pMask(ip1,jm1,k,bi,bj) +
971       &                     fld_in(ip1,jp1)*pMask(ip1,jp1,k,bi,bj)))       &                     fld(ip1,jp1)*pMask(ip1,jp1,k,bi,bj)))
972       &              / tempVar       &              / tempVar
973              ELSE              ELSE
974                 fld_tmp(i,j) = fld_in(i,j)                 fld_tmp(i,j) = fld(i,j)
975              ENDIF              ENDIF
976           ENDDO           ENDDO
977        ENDDO        ENDDO
# Line 1006  c     local Line 979  c     local
979  c     transfer smoothed field to output array  c     transfer smoothed field to output array
980        DO j = jmin, jmax        DO j = jmin, jmax
981           DO i = imin, imax           DO i = imin, imax
982              fld_out(i,j) = fld_tmp(i,j)              fld(i,j) = fld_tmp(i,j)
983           ENDDO           ENDDO
984        ENDDO        ENDDO
985    
 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  
   
986  #endif /* ALLOW_KPP */  #endif /* ALLOW_KPP */
987    
988        return        return
# Line 1037  c*************************************** Line 992  c***************************************
992    
993        subroutine blmix (        subroutine blmix (
994       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl
995       O     , dkm1, blmc, ghat, sigma       O     , dkm1, blmc, ghat, sigma, ikey
996       &     )       &     )
997    
998  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 1016  c     stable(imt)                 = 1 in
1016  c     casea (imt)                 = 1 in case A  c     casea (imt)                 = 1 in case A
1017  c     diffus(imt,0:Nrp1,mdiff)    vertical diffusivities              (m^2/s)  c     diffus(imt,0:Nrp1,mdiff)    vertical diffusivities              (m^2/s)
1018  c     kbl(imt)                    -1 of first grid level below hbl  c     kbl(imt)                    -1 of first grid level below hbl
1019        _RS ustar (imt)        _KPP_RL ustar (imt)
1020        _RS bfsfc (imt)        _KPP_RL bfsfc (imt)
1021        _RS hbl   (imt)        _KPP_RL hbl   (imt)
1022        _RS stable(imt)        _KPP_RL stable(imt)
1023        _RS casea (imt)        _KPP_RL casea (imt)
1024        _RS diffus(imt,0:Nrp1,mdiff)        _KPP_RL diffus(imt,0:Nrp1,mdiff)
1025        integer kbl(imt)        integer kbl(imt)
1026    
1027  c output  c output
# Line 1074  c     dkm1 (imt,mdiff)            bounda Line 1029  c     dkm1 (imt,mdiff)            bounda
1029  c     blmc (imt,Nr,mdiff)         boundary layer mixing coefficients  (m^2/s)  c     blmc (imt,Nr,mdiff)         boundary layer mixing coefficients  (m^2/s)
1030  c     ghat (imt,Nr)               nonlocal scalar transport  c     ghat (imt,Nr)               nonlocal scalar transport
1031  c     sigma(imt)                  normalized depth (d / hbl)  c     sigma(imt)                  normalized depth (d / hbl)
1032        _RS dkm1 (imt,mdiff)        _KPP_RL dkm1 (imt,mdiff)
1033        _RS blmc (imt,Nr,mdiff)        _KPP_RL blmc (imt,Nr,mdiff)
1034        _RS ghat (imt,Nr)        _KPP_RL ghat (imt,Nr)
1035        _RS sigma(imt)        _KPP_RL sigma(imt)
1036          integer ikey
1037    
1038  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1039    
# Line 1085  c  local Line 1041  c  local
1041  c     gat1*(imt)                 shape function at sigma = 1  c     gat1*(imt)                 shape function at sigma = 1
1042  c     dat1*(imt)                 derivative of shape function at sigma = 1  c     dat1*(imt)                 derivative of shape function at sigma = 1
1043  c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)  c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)
1044        _RS gat1m(imt), gat1s(imt), gat1t(imt)        _KPP_RL gat1m(imt), gat1s(imt), gat1t(imt)
1045        _RS dat1m(imt), dat1s(imt), dat1t(imt)        _KPP_RL dat1m(imt), dat1s(imt), dat1t(imt)
1046        _RS ws(imt), wm(imt)        _KPP_RL ws(imt), wm(imt)
1047        integer i, kn, ki        integer i, kn, ki
1048        _RS     R, dvdzup, dvdzdn, viscp        _KPP_RL R, dvdzup, dvdzdn, viscp
1049        _RS     difsp, diftp, visch, difsh, difth        _KPP_RL difsp, diftp, visch, difsh, difth
1050        _RS     f1, sig, a1, a2, a3, delhat        _KPP_RL f1, sig, a1, a2, a3, delhat
1051        _RS     Gm, Gs, Gt        _KPP_RL Gm, Gs, Gt
1052        _RL     dum        _KPP_RL tempVar
1053    
1054        _RS        p0    , eins        _KPP_RL    p0    , eins
1055        parameter (p0=0.0, eins=1.0)        parameter (p0=0.0, eins=1.0)
1056    
1057  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 1111  c--------------------------------------- Line 1067  c---------------------------------------
1067       O        wm, ws )       O        wm, ws )
1068    
1069        do i = 1, imt        do i = 1, imt
1070             wm(i) = sign(eins,wm(i))*max(phepsi,abs(wm(i)))
1071             ws(i) = sign(eins,ws(i))*max(phepsi,abs(ws(i)))
1072          end do
1073    CADJ STORE wm = comlev1_kpp, key = ikey
1074    CADJ STORE ws = comlev1_kpp, key = ikey
1075    
1076          do i = 1, imt
1077    
1078           kn = int(caseA(i)+phepsi) *(kbl(i) -1) +           kn = int(caseA(i)+phepsi) *(kbl(i) -1) +
1079       $        (1 - int(caseA(i)+phepsi)) * kbl(i)       $        (1 - int(caseA(i)+phepsi)) * kbl(i)
# Line 1140  c--------------------------------------- Line 1103  c---------------------------------------
1103           difsh  = diffus(i,kn,2) + difsp * delhat           difsh  = diffus(i,kn,2) + difsp * delhat
1104           difth  = diffus(i,kn,3) + diftp * delhat           difth  = diffus(i,kn,3) + diftp * delhat
1105    
 #ifdef KPP_TEST_DENOM  
 ctl replace (Important!!! not phepsi**4 !!!)  
1106           f1 = stable(i) * conc1 * bfsfc(i) /           f1 = stable(i) * conc1 * bfsfc(i) /
1107       &        max(ustar(i)**4,phepsi)       &        max(ustar(i)**4,phepsi)
          wm(i) = sign(eins,wm(i))*max(phepsi,abs(wm(i)))  
1108           gat1m(i) = visch / hbl(i) / wm(i)           gat1m(i) = visch / hbl(i) / wm(i)
1109           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  
1110           dat1m(i) = min(dat1m(i),p0)           dat1m(i) = min(dat1m(i),p0)
1111    
 #ifdef KPP_TEST_DENOM  
          ws(i) = sign(eins,ws(i))*max(phepsi,abs(ws(i)))  
1112           gat1s(i) = difsh  / hbl(i) / ws(i)           gat1s(i) = difsh  / hbl(i) / ws(i)
1113           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  
1114           dat1s(i) = min(dat1s(i),p0)           dat1s(i) = min(dat1s(i),p0)
1115    
 #ifdef KPP_TEST_DENOM  
1116           gat1t(i) = difth /  hbl(i) / ws(i)           gat1t(i) = difth /  hbl(i) / ws(i)
1117           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  
1118           dat1t(i) = min(dat1t(i),p0)           dat1t(i) = min(dat1t(i),p0)
1119    
1120        end do        end do
# Line 1214  c--------------------------------------- Line 1158  c---------------------------------------
1158  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1159  c     nonlocal transport term = ghat * <ws>o  c     nonlocal transport term = ghat * <ws>o
1160  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1161  #ifdef KPP_TEST_DENOM  
1162              dum = ws(i) * hbl(i)              tempVar = ws(i) * hbl(i)
1163  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  
1164    
1165           end do           end do
1166        end do        end do
# Line 1280  c     hbl(imt)                  boundary Line 1218  c     hbl(imt)                  boundary
1218  c     kbl(imt)                  grid above hbl  c     kbl(imt)                  grid above hbl
1219  c     diffus(imt,0:Nrp1,mdiff) vertical diffusivities           (m^2/s)  c     diffus(imt,0:Nrp1,mdiff) vertical diffusivities           (m^2/s)
1220  c     casea(imt)                = 1 in caseA, = 0 in case B  c     casea(imt)                = 1 in caseA, = 0 in case B
1221        _RS     dkm1  (imt,mdiff)        _KPP_RL dkm1  (imt,mdiff)
1222        _RS     hbl   (imt)        _KPP_RL hbl   (imt)
1223        integer kbl   (imt)        integer kbl   (imt)
1224        _RS     diffus(imt,0:Nrp1,mdiff)        _KPP_RL diffus(imt,0:Nrp1,mdiff)
1225        _RS     casea (imt)        _KPP_RL casea (imt)
1226    
1227  c input/output  c input/output
1228  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)
1229        _RS     ghat (imt,Nr)        _KPP_RL ghat (imt,Nr)
1230    
1231  c output  c output
1232  c     enhanced bound. layer mixing coeff.  c     enhanced bound. layer mixing coeff.
1233        _RS     blmc  (imt,Nr,mdiff)        _KPP_RL blmc  (imt,Nr,mdiff)
1234    
1235  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1236    
1237  c local  c local
1238  c     fraction hbl lies beteen zgrid neighbors  c     fraction hbl lies beteen zgrid neighbors
1239        _RS delta        _KPP_RL delta
1240        integer ki, i, md        integer ki, i, md
1241        _RS dkmp5, dstar        _KPP_RL dkmp5, dstar
1242    
1243        do i = 1, imt        do i = 1, imt
1244           ki = kbl(i)-1           ki = kbl(i)-1
# Line 1360  c--------------------------------------- Line 1298  c---------------------------------------
1298  #include "EEPARAMS.h"  #include "EEPARAMS.h"
1299  #include "PARAMS.h"  #include "PARAMS.h"
1300  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
1301    #include "DYNVARS.h"
1302    
1303  c-------------- Routine arguments -----------------------------------------  c-------------- Routine arguments -----------------------------------------
1304        INTEGER bi, bj, myThid        INTEGER bi, bj, myThid
1305  #ifdef FRUGAL_KPP        _KPP_RL RHO1   ( ibot:itop, jbot:jtop       )
1306        _RS RHO1   (sNx,sNy)        _KPP_RL DBLOC  ( ibot:itop, jbot:jtop, Nr   )
1307        _RS DBLOC  (sNx,sNy,Nr)        _KPP_RL DBSFC  ( ibot:itop, jbot:jtop, Nr   )
1308        _RS DBSFC  (sNx,sNy,Nr)        _KPP_RL TTALPHA( ibot:itop, jbot:jtop, Nrp1 )
1309        _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  
1310    
1311  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1312    
# Line 1399  c     work1, work2 - work arrays for hol Line 1330  c     work1, work2 - work arrays for hol
1330  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
1331    
1332        call FIND_RHO(        call FIND_RHO(
1333  #ifdef FRUGAL_KPP       I     bi, bj, ibot, itop, jbot, jtop, 1, 1, eosType,
1334       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  
1335       O     WORK1,       O     WORK1,
1336       I     myThid )       I     myThid )
1337    
1338        call FIND_ALPHA(        call FIND_ALPHA(
1339  #ifdef FRUGAL_KPP       I     bi, bj, ibot, itop, jbot, jtop, 1, 1, eosType,
      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  
1340       O     WORK2 )       O     WORK2 )
1341    
1342        call FIND_BETA(        call FIND_BETA(
1343  #ifdef FRUGAL_KPP       I     bi, bj, ibot, itop, jbot, jtop, 1, 1, eosType,
      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  
1344       O     WORK3 )       O     WORK3 )
1345    
1346  #ifdef FRUGAL_KPP        DO J = jbot, jtop
1347        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  
1348              RHO1(I,J)      = WORK1(I,J) + rhonil              RHO1(I,J)      = WORK1(I,J) + rhonil
1349              TTALPHA(I,J,1) = WORK2(I,J)              TTALPHA(I,J,1) = WORK2(I,J)
1350              SSBETA(I,J,1)  = WORK3(I,J)              SSBETA(I,J,1)  = WORK3(I,J)
# Line 1443  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO Line 1358  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO
1358        DO K = 2, Nr        DO K = 2, Nr
1359    
1360           call FIND_RHO(           call FIND_RHO(
1361  #ifdef FRUGAL_KPP       I        bi, bj, ibot, itop, jbot, jtop, K, K, eosType,
1362       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  
1363       O        RHOK,       O        RHOK,
1364       I        myThid )       I        myThid )
1365    
1366           call FIND_RHO(           call FIND_RHO(
1367  #ifdef FRUGAL_KPP       I        bi, bj, ibot, itop, jbot, jtop, K-1, K, eosType,
1368       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  
1369       O        RHOKM1,       O        RHOKM1,
1370       I        myThid )       I        myThid )
1371    
1372           call FIND_RHO(           call FIND_RHO(
1373  #ifdef FRUGAL_KPP       I        bi, bj, ibot, itop, jbot, jtop, 1, K, eosType,
1374       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  
1375       O        RHO1K,       O        RHO1K,
1376       I        myThid )       I        myThid )
1377    
1378           call FIND_ALPHA(           call FIND_ALPHA(
1379  #ifdef FRUGAL_KPP       I        bi, bj, ibot, itop, jbot, jtop, K, K, eosType,
      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  
1380       O        WORK1 )       O        WORK1 )
1381    
1382           call FIND_BETA(           call FIND_BETA(
1383  #ifdef FRUGAL_KPP       I        bi, bj, ibot, itop, jbot, jtop, K, K, eosType,
      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  
1384       O        WORK2 )       O        WORK2 )
1385    
1386  #ifdef FRUGAL_KPP           DO J = jbot, jtop
1387           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  
1388                 TTALPHA(I,J,K) = WORK1 (I,J)                 TTALPHA(I,J,K) = WORK1 (I,J)
1389                 SSBETA(I,J,K)  = WORK2 (I,J)                 SSBETA(I,J,K)  = WORK2 (I,J)
1390                 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 1504  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO Line 1397  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO
1397        END DO        END DO
1398    
1399  c     compute arrays for K = Nrp1  c     compute arrays for K = Nrp1
1400  #ifdef FRUGAL_KPP        DO J = jbot, jtop
1401        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  
1402              TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)              TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)
1403              SSBETA(I,J,Nrp1)  = SSBETA(I,J,Nr)              SSBETA(I,J,Nrp1)  = SSBETA(I,J,Nr)
1404              DBLOC(I,J,Nr)     = 0.              DBLOC(I,J,Nr)     = 0.
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.22