/[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.4 by heimbach, Mon Nov 13 16:37:02 2000 UTC
# Line 10  C--  o BLDEPTH     - Determine oceanic p Line 10  C--  o BLDEPTH     - Determine oceanic p
10  C--  o WSCALE      - Compute turbulent velocity scales.  C--  o WSCALE      - Compute turbulent velocity scales.
11  C--  o RI_IWMIX    - Compute interior viscosity diffusivity coefficients.  C--  o RI_IWMIX    - Compute interior viscosity diffusivity coefficients.
12  C--  o Z121        - Apply 121 vertical smoothing.  C--  o Z121        - Apply 121 vertical smoothing.
13  C--  o SMOOTH_HORIZ_RS - Apply horizontal smoothing to RS array.  C--  o KPP_SMOOTH_HORIZ - Apply horizontal smoothing to KPP array.
14  C--  o SMOOTH_HORIZ_RL - Apply horizontal smoothing to RL array.  C--  o SMOOTH_HORIZ     - Apply horizontal smoothing to global array.
15  C--  o BLMIX       - Boundary layer mixing coefficients.  C--  o BLMIX       - Boundary layer mixing coefficients.
16  C--  o ENHANCE     - Enhance diffusivity at boundary layer interface.  C--  o ENHANCE     - Enhance diffusivity at boundary layer interface.
17  C--  o STATEKPP    - Compute buoyancy-related input arrays.  C--  o STATEKPP    - Compute buoyancy-related input arrays.
18    
19  c*************************************************************************  c***********************************************************************
20    
21        SUBROUTINE KPPMIX (        SUBROUTINE KPPMIX (
22       I     lri, kmtj, shsq, dvsq, ustar, bo, bosol       I       mytime, mythid
23       I     , dbloc, Ritop, coriol       I     , kmtj, shsq, dvsq, ustar
24         I     , bo, bosol, dbloc, Ritop, coriol
25       I     , ikey       I     , ikey
26       O     , diffus       O     , diffus
27       U     , ghat       U     , ghat
28       O     , hbl )       O     , hbl )
29    
30  c-------------------------------------------------------------------------  c-----------------------------------------------------------------------
31  c  c
32  c     Main driver subroutine for kpp vertical mixing scheme and  c     Main driver subroutine for kpp vertical mixing scheme and
33  c     interface to greater ocean model  c     interface to greater ocean model
# Line 50  c--------------------------------------- Line 51  c---------------------------------------
51  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
52    
53  c input  c input
54  c     lri             - mixing process switches  c     myTime          - current time in simulation
55    c     myThid          - thread number for this instance of the routine
56  c     kmtj   (imt)    - number of vertical layers on this row  c     kmtj   (imt)    - number of vertical layers on this row
57  c     shsq   (imt,Nr) - (local velocity shear)^2                      ((m/s)^2)  c     shsq   (imt,Nr) - (local velocity shear)^2                      ((m/s)^2)
58  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 69  c     note: there is a conversion from 2
69  c           e.g., hbl(sNx,sNy) -> hbl(imt),  c           e.g., hbl(sNx,sNy) -> hbl(imt),
70  c           where hbl(i,j) -> hbl((j-1)*sNx+i)  c           where hbl(i,j) -> hbl((j-1)*sNx+i)
71    
72        logical lri        _RL     mytime
73          integer mythid
74        integer kmtj  (imt   )        integer kmtj  (imt   )
75        _RS    shsq   (imt,Nr)        _KPP_RL shsq  (imt,Nr)
76        _RS    dvsq   (imt,Nr)        _KPP_RL dvsq  (imt,Nr)
77        _RS    ustar  (imt   )        _KPP_RL ustar (imt   )
78        _RS    bo     (imt   )        _KPP_RL bo    (imt   )
79        _RS    bosol  (imt   )        _KPP_RL bosol (imt   )
80        _RS    dbloc  (imt,Nr)        _KPP_RL dbloc (imt,Nr)
81        _RS    Ritop  (imt,Nr)        _KPP_RL Ritop (imt,Nr)
82        _RS    coriol (imt   )        _KPP_RL coriol(imt   )
83    
84        integer ikey        integer ikey
85    
# Line 87  c     diffus (imt,3)  - vertical tempera Line 90  c     diffus (imt,3)  - vertical tempera
90  c     ghat   (imt)    - nonlocal transport coefficient                  (s/m^2)  c     ghat   (imt)    - nonlocal transport coefficient                  (s/m^2)
91  c     hbl    (imt)    - mixing layer depth                                  (m)  c     hbl    (imt)    - mixing layer depth                                  (m)
92    
93        _RS diffus(imt,0:Nrp1,mdiff)        _KPP_RL diffus(imt,0:Nrp1,mdiff)
94        _RS ghat  (imt,Nr)        _KPP_RL ghat  (imt,Nr)
95        _RS hbl   (imt)        _KPP_RL hbl   (imt)
96    
97  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
98    
# Line 103  c     blmc   (imt,Nr,mdiff) - boundary l Line 106  c     blmc   (imt,Nr,mdiff) - boundary l
106  c     sigma  (imt         ) - normalized depth (d / hbl)  c     sigma  (imt         ) - normalized depth (d / hbl)
107  c     Rib    (imt,Nr      ) - bulk Richardson number  c     Rib    (imt,Nr      ) - bulk Richardson number
108    
109        integer kbl (imt         )        integer kbl   (imt         )
110        _RS bfsfc   (imt         )        _KPP_RL bfsfc (imt         )
111        _RS casea   (imt         )        _KPP_RL casea (imt         )
112        _RS stable  (imt         )        _KPP_RL stable(imt         )
113        _RS dkm1    (imt,   mdiff)        _KPP_RL dkm1  (imt,   mdiff)
114        _RS blmc    (imt,Nr,mdiff)        _KPP_RL blmc  (imt,Nr,mdiff)
115        _RS sigma   (imt         )        _KPP_RL sigma (imt         )
116        _RS Rib     (imt,Nr      )        _KPP_RL Rib   (imt,Nr      )
117    
118        integer i, k, md        integer i, k, md
119    
# Line 121  c instability. Line 124  c instability.
124  c (ghat is temporary storage for horizontally smoothed dbloc)  c (ghat is temporary storage for horizontally smoothed dbloc)
125  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
126    
127    CADJ STORE ghat = comlev1_kpp, key = ikey
128    
129        call Ri_iwmix (        call Ri_iwmix (
130       I       kmtj, shsq, dbloc, ghat       I       kmtj, shsq, dbloc, ghat
131       I     , ikey       I     , ikey
# Line 146  c diagnose the new boundary layer depth Line 151  c diagnose the new boundary layer depth
151  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
152    
153        call bldepth (        call bldepth (
154       I       kmtj       I       mytime, mythid
155         I     , kmtj
156       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol
157       I     , ikey       I     , ikey
158       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma
# Line 160  c--------------------------------------- Line 166  c---------------------------------------
166    
167        call blmix (        call blmix (
168       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl
169       O     , dkm1, blmc, ghat, sigma       O     , dkm1, blmc, ghat, sigma, ikey
170       &     )       &     )
171    
172  CADJ STORE dkm1,blmc,ghat = comlev1_kpp, key = ikey  CADJ STORE dkm1,blmc,ghat = comlev1_kpp, key = ikey
# Line 198  c--------------------------------------- Line 204  c---------------------------------------
204  c*************************************************************************  c*************************************************************************
205    
206        subroutine bldepth (        subroutine bldepth (
207       I       kmtj       I       mytime, mythid
208         I     , kmtj
209       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol
210       I     , ikey       I     , ikey
211       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma
# Line 235  c Line 242  c
242    
243  c input  c input
244  c------  c------
245    c myTime    : current time in simulation
246    c myThid    : thread number for this instance of the routine
247  c kmtj      : number of vertical layers  c kmtj      : number of vertical layers
248  c dvsq      : (velocity shear re sfc)^2             ((m/s)^2)  c dvsq      : (velocity shear re sfc)^2             ((m/s)^2)
249  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 254  c ustar     : surface friction velocity
254  c bo        : surface turbulent buoyancy forcing    (m^2/s^3)  c bo        : surface turbulent buoyancy forcing    (m^2/s^3)
255  c bosol     : radiative buoyancy forcing            (m^2/s^3)  c bosol     : radiative buoyancy forcing            (m^2/s^3)
256  c coriol    : Coriolis parameter                        (1/s)  c coriol    : Coriolis parameter                        (1/s)
257          _RL     mytime
258          integer mythid
259        integer kmtj(imt)        integer kmtj(imt)
260        _RS dvsq   (imt,Nr)        _KPP_RL dvsq  (imt,Nr)
261        _RS dbloc  (imt,Nr)        _KPP_RL dbloc (imt,Nr)
262        _RS Ritop  (imt,Nr)        _KPP_RL Ritop (imt,Nr)
263        _RS ustar  (imt)        _KPP_RL ustar (imt)
264        _RS bo     (imt)        _KPP_RL bo    (imt)
265        _RS bosol  (imt)        _KPP_RL bosol (imt)
266        _RS coriol (imt)        _KPP_RL coriol(imt)
267        integer ikey        integer ikey
268    
269  c  output  c  output
# Line 264  c casea     : =1 in case A, =0 in case B Line 275  c casea     : =1 in case A, =0 in case B
275  c kbl       : -1 of first grid level below hbl  c kbl       : -1 of first grid level below hbl
276  c Rib       : Bulk Richardson number  c Rib       : Bulk Richardson number
277  c sigma     : normalized depth (d/hbl)  c sigma     : normalized depth (d/hbl)
278        _RS hbl    (imt)        _KPP_RL hbl   (imt)
279        _RS bfsfc  (imt)        _KPP_RL bfsfc (imt)
280        _RS stable (imt)        _KPP_RL stable(imt)
281        _RS casea  (imt)        _KPP_RL casea (imt)
282        integer kbl (imt)        integer kbl   (imt)
283        _RS Rib    (imt,Nr)        _KPP_RL Rib   (imt,Nr)
284        _RS sigma  (imt)        _KPP_RL sigma (imt)
285    
286  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
287    
288  c  local  c  local
289  c-------  c-------
290  c wm, ws    : turbulent velocity scales         (m/s)  c wm, ws    : turbulent velocity scales         (m/s)
291  c zwork     : depth work array        _KPP_RL wm(imt), ws(imt)
292        _RS wm(imt), ws(imt)        _RL     worka(imt)
       _RS zwork(imt)  
293    
294        _RS bvsq, vtsq, hekman, hmonob, hlimit        _KPP_RL bvsq, vtsq, hekman, hmonob, hlimit, tempVar1, tempVar2
295        integer i, kl        integer i, kl
296    
297        _RS        p5    , eins    , m1        _KPP_RL     p5    , eins
298        parameter (p5=0.5, eins=1.0, m1=-1.0 )        parameter ( p5=0.5, eins=1.0 )
299          _RL         minusone
300  crg intermediate result in RL to avoid overflow in adjoint        parameter ( minusone=-1.0 )
        _RL        dpshear2  
 #ifdef KPP_TEST_DENOM  
        _RL        dhelp  
 #endif  
301    
302  c find bulk Richardson number at every grid level until > Ricr  c find bulk Richardson number at every grid level until > Ricr
303  c  c
# Line 315  c     initialize hbl and kbl to bottomed Line 321  c     initialize hbl and kbl to bottomed
321  c     compute bfsfc = sw fraction at hbf * zgrid  c     compute bfsfc = sw fraction at hbf * zgrid
322    
323           do i = 1, imt           do i = 1, imt
324              zwork(i) = zgrid(kl)              worka(i) = zgrid(kl)
325           end do           end do
326           call SWFRAC(           call SWFRAC(
327       I        imt, hbf, zwork,       I        imt, hbf,
328       O        bfsfc)       I        mytime, mythid,
329         U        worka )
330    
331           do i = 1, imt           do i = 1, imt
332    
# Line 329  c     use caseA as temporary array Line 336  c     use caseA as temporary array
336    
337  c     compute bfsfc= Bo + radiative contribution down to hbf * hbl  c     compute bfsfc= Bo + radiative contribution down to hbf * hbl
338    
339              bfsfc(i) = bo(i) + bosol(i)*(1. - bfsfc(i))              bfsfc(i) = bo(i) + bosol(i)*(1. - worka(i))
340              stable(i) = p5 + sign(p5,bfsfc(i))              stable(i) = p5 + sign(p5,bfsfc(i))
341              sigma(i) = stable(i) + (1. - stable(i)) * epsilon              sigma(i) = stable(i) + (1. - stable(i)) * epsilon
342    
# Line 370  c Line 377  c
377  c     rg: assignment to double precision variable to avoid overflow  c     rg: assignment to double precision variable to avoid overflow
378  c     ph: test for zero nominator  c     ph: test for zero nominator
379  c  c
380  #ifdef KPP_TEST_DENOM  
381  ctl replace              tempVar1  = dvsq(i,kl) + vtsq          
382              dpshear2  = dvsq(i,kl) + vtsq                        tempVar2 = max(tempVar1, phepsi)
383              dpshear2 = max(dpshear2, phepsi)              Rib(i,kl) = Ritop(i,kl) / tempVar2
384              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  
385           end do           end do
386        end do        end do
387                
# Line 394  CADJ &   , key = ikey, shape = (/ (sNx+2 Line 396  CADJ &   , key = ikey, shape = (/ (sNx+2
396    
397        do i = 1, imt        do i = 1, imt
398           kl = kbl(i)           kl = kbl(i)
   
399  c     linearly interpolate to find hbl where Rib = Ricr  c     linearly interpolate to find hbl where Rib = Ricr
   
400           if (kl.gt.1 .and. kl.lt.kmtj(i)) then           if (kl.gt.1 .and. kl.lt.kmtj(i)) then
401  #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  
402              hbl(i) = -zgrid(kl-1) + (zgrid(kl-1)-zgrid(kl)) *              hbl(i) = -zgrid(kl-1) + (zgrid(kl-1)-zgrid(kl)) *
403       1           (Ricr - Rib(i,kl-1)) / (Rib(i,kl)-Rib(i,kl-1))       1           (Ricr - Rib(i,kl-1)) / tempVar1
 #endif  
404           endif           endif
405        end do        end do
406    
# Line 416  c--------------------------------------- Line 411  c---------------------------------------
411  c     find stability and buoyancy forcing for boundary layer  c     find stability and buoyancy forcing for boundary layer
412  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
413    
414          do i = 1, imt
415             worka(i) = hbl(i)
416          end do
417        call SWFRAC(        call SWFRAC(
418       I     imt, m1, hbl,       I     imt, minusone,
419       O     bfsfc)       I     mytime, mythid,
420         U     worka )
421    
422          do i = 1, imt
423             bfsfc(i)  = bo(i) + bosol(i) * (1. - worka(i))
424          end do
425  CADJ store bfsfc = comlev1_kpp  CADJ store bfsfc = comlev1_kpp
426  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
427    
428    c--   ensure bfsfc is never 0
429        do i = 1, imt        do i = 1, imt
          bfsfc(i)  = bo(i) + bosol(i) * (1. - bfsfc(i))  
430           stable(i) = p5 + sign( p5, bfsfc(i) )           stable(i) = p5 + sign( p5, bfsfc(i) )
431             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  
   
432        end do        end do
433    
434    CADJ store bfsfc = comlev1_kpp
435    CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
436    
437  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
438  c check hbl limits for hekman or hmonob  c check hbl limits for hekman or hmonob
439  c     ph: test for zero nominator  c     ph: test for zero nominator
# Line 445  c--------------------------------------- Line 441  c---------------------------------------
441    
442        do i = 1, imt        do i = 1, imt
443           if (bfsfc(i) .gt. 0.0) then           if (bfsfc(i) .gt. 0.0) then
 #ifdef KPP_TEST_DENOM  
 ctl replace  
444              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  
445              hmonob = cmonob * ustar(i)*ustar(i)*ustar(i)              hmonob = cmonob * ustar(i)*ustar(i)*ustar(i)
446       &           / vonk / bfsfc(i)       &           / vonk / bfsfc(i)
447              hlimit = stable(i) * min(hekman,hmonob)              hlimit = stable(i) * min(hekman,hmonob)
448       &             + (stable(i)-1.) * zgrid(Nr)       &             + (stable(i)-1.) * zgrid(Nr)
449              hbl(i) = min(hbl(i),hlimit)              hbl(i) = min(hbl(i),hlimit)
450              hbl(i) = max(hbl(i),minKPPhbl)           end if
451           endif        end do
452    CADJ store hbl = comlev1_kpp
453    CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
454    
455          do i = 1, imt
456             hbl(i) = max(hbl(i),minKPPhbl)
457           kbl(i) = kmtj(i)           kbl(i) = kmtj(i)
458        end do        end do
459    
# Line 481  c--------------------------------------- Line 476  c---------------------------------------
476  c      find stability and buoyancy forcing for final hbl values  c      find stability and buoyancy forcing for final hbl values
477  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
478    
479          do i = 1, imt
480             worka(i) = hbl(i)
481          end do
482        call SWFRAC(        call SWFRAC(
483       I     imt, m1, hbl,       I     imt, minusone,
484       O     bfsfc)       I     mytime, mythid,
485         U     worka )
486    
487          do i = 1, imt
488             bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i))
489          end do
490  CADJ store bfsfc = comlev1_kpp  CADJ store bfsfc = comlev1_kpp
491  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
492    
493    c--   ensures bfsfc is never 0
494        do i = 1, imt        do i = 1, imt
          bfsfc(i) = bo(i) + bosol(i) * (1. - bfsfc(i))  
495           stable(i) = p5 + sign( p5, bfsfc(i) )           stable(i) = p5 + sign( p5, bfsfc(i) )
496  #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  
497        end do        end do
498    
499  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 540  c sigma   : normalized depth (d/hbl) Line 535  c sigma   : normalized depth (d/hbl)
535  c hbl     : boundary layer depth (m)  c hbl     : boundary layer depth (m)
536  c ustar   : surface friction velocity         (m/s)  c ustar   : surface friction velocity         (m/s)
537  c bfsfc   : total surface buoyancy flux       (m^2/s^3)  c bfsfc   : total surface buoyancy flux       (m^2/s^3)
538        _RS sigma(imt)        _KPP_RL sigma(imt)
539        _RS hbl  (imt)        _KPP_RL hbl  (imt)
540        _RS ustar(imt)        _KPP_RL ustar(imt)
541        _RS bfsfc(imt)        _KPP_RL bfsfc(imt)
542    
543  c  output  c  output
544  c--------  c--------
545  c wm, ws  : turbulent velocity scales at sigma  c wm, ws  : turbulent velocity scales at sigma
546        _RS wm(imt), ws(imt)        _KPP_RL wm(imt), ws(imt)
547    
548  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
549    
550  c local  c local
551  c------  c------
552  c zehat   : = zeta *  ustar**3  c zehat   : = zeta *  ustar**3
553        _RS zehat        _KPP_RL zehat
554    
555        integer iz, izp1, ju, i, jup1        integer iz, izp1, ju, i, jup1
556        _RS udiff, zdiff, zfrac, ufrac, fzfrac, wam, wbm, was, wbs, u3        _KPP_RL udiff, zdiff, zfrac, ufrac, fzfrac, wam
557        _RL dum        _KPP_RL wbm, was, wbs, u3, tempVar
558    
559  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
560  c use lookup table for zehat < zmax only; otherwise use  c use lookup table for zehat < zmax only; otherwise use
# Line 576  c--------------------------------------- Line 571  c---------------------------------------
571              iz    = min( iz, nni )              iz    = min( iz, nni )
572              iz    = max( iz, 0 )              iz    = max( iz, 0 )
573              izp1  = iz + 1              izp1  = iz + 1
574    
575              udiff = ustar(i) - umin              udiff = ustar(i) - umin
576              ju    = int( udiff / deltau )              ju    = int( udiff / deltau )
577              ju    = min( ju, nnj )              ju    = min( ju, nnj )
578              ju    = max( ju, 0 )              ju    = max( ju, 0 )
579              jup1  = ju + 1              jup1  = ju + 1
580    
581              zfrac = zdiff / deltaz - float(iz)              zfrac = zdiff / deltaz - float(iz)
582              ufrac = udiff / deltau - float(ju)              ufrac = udiff / deltau - float(ju)
583    
584              fzfrac= 1. - zfrac              fzfrac= 1. - zfrac
585              wam   = fzfrac     * wmt(iz,jup1) + zfrac * wmt(izp1,jup1)              wam   = fzfrac     * wmt(iz,jup1) + zfrac * wmt(izp1,jup1)
586              wbm   = fzfrac     * wmt(iz,ju  ) + zfrac * wmt(izp1,ju  )              wbm   = fzfrac     * wmt(iz,ju  ) + zfrac * wmt(izp1,ju  )
587              wm(i) = (1.-ufrac) * wbm          + ufrac * wam              wm(i) = (1.-ufrac) * wbm          + ufrac * wam
588    
589              was   = fzfrac     * wst(iz,jup1) + zfrac * wst(izp1,jup1)              was   = fzfrac     * wst(iz,jup1) + zfrac * wst(izp1,jup1)
590              wbs   = fzfrac     * wst(iz,ju  ) + zfrac * wst(izp1,ju  )              wbs   = fzfrac     * wst(iz,ju  ) + zfrac * wst(izp1,ju  )
591              ws(i) = (1.-ufrac) * wbs          + ufrac * was              ws(i) = (1.-ufrac) * wbs          + ufrac * was
592    
593           else           else
594    
595              u3    = ustar(i) * ustar(i) * ustar(i)              u3 = ustar(i) * ustar(i) * ustar(i)
596              dum   = u3 + conc1 * zehat              tempVar = u3 + conc1 * zehat
597              wm(i) = vonk     * ustar(i) * u3 / dum              wm(i) = vonk * ustar(i) * u3 / tempVar
598              ws(i) = wm(i)              ws(i) = wm(i)
599    
600           endif           endif
601    
602        end do        end do
# Line 636  c     shsq   (imt,Nr)      (local veloci Line 631  c     shsq   (imt,Nr)      (local veloci
631  c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)  c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)
632  c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)  c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)
633        integer kmtj   (imt)        integer kmtj   (imt)
634        _RS     shsq   (imt,Nr)        _KPP_RL shsq   (imt,Nr)
635        _RS     dbloc  (imt,Nr)        _KPP_RL dbloc  (imt,Nr)
636        _RS     dblocSm(imt,Nr)        _KPP_RL dblocSm(imt,Nr)
637        integer ikey        integer ikey
638    
639  c output  c output
640  c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)  c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)
641  c     diffus(imt,0:Nrp1,2)  vertical scalar diffusivity             (m^2/s)  c     diffus(imt,0:Nrp1,2)  vertical scalar diffusivity             (m^2/s)
642  c     diffus(imt,0:Nrp1,3)  vertical temperature diffusivity        (m^2/s)  c     diffus(imt,0:Nrp1,3)  vertical temperature diffusivity        (m^2/s)
643        _RS diffus(imt,0:Nrp1,3)        _KPP_RL diffus(imt,0:Nrp1,3)
644    
645  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
646    
647  c local variables  c local variables
648  c     Rig                   local Richardson number  c     Rig                   local Richardson number
649  c     fRi, fcon             function of Rig  c     fRi, fcon             function of Rig
650        _RS Rig        _KPP_RL Rig
651        _RS fRi, fcon        _KPP_RL fRi, fcon
652        _RS ratio        _KPP_RL ratio
653        integer i, ki, mr        integer i, ki, mr
654        _RS c1, c0        _KPP_RL c1, c0
655    
656  c constants  c constants
657        c1 = 1.0        c1 = 1.0
# Line 678  c     set values at bottom and below to Line 673  c     set values at bottom and below to
673                 diffus(i,ki,2) = diffus(i,ki-1,2)                 diffus(i,ki,2) = diffus(i,ki-1,2)
674              else              else
675                 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  
676       &              / max( Shsq(i,ki), phepsi )       &              / max( Shsq(i,ki), phepsi )
 #else  
      &              / ( shsq(i,ki) + epsln )  
 #endif  
677                 diffus(i,ki,2) = dbloc(i,ki)   / (zgrid(ki)-zgrid(ki+1))                 diffus(i,ki,2) = dbloc(i,ki)   / (zgrid(ki)-zgrid(ki+1))
678              endif              endif
679           end do           end do
# Line 769  c     penetrative convection. Line 760  c     penetrative convection.
760  c input/output  c input/output
761  c-------------  c-------------
762  c v     : 2-D array to be smoothed in Nrp1 direction  c v     : 2-D array to be smoothed in Nrp1 direction
763        _RS v(imt,0:Nrp1)        _KPP_RL v(imt,0:Nrp1)
764    
765  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
766    
767  c local  c local
768        _RS zwork, zflag        _KPP_RL zwork, zflag
769        _RS KRi_range(1:Nrp1)        _KPP_RL KRi_range(1:Nrp1)
770        integer i, k, km1, kp1        integer i, k, km1, kp1
771    
772        _RS         p0      , p25       , p5      , p2        _KPP_RL     p0      , p25       , p5      , p2
773        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 )
774    
775        KRi_range(Nrp1) = p0        KRi_range(Nrp1) = p0
# Line 789  C--   dummy assignment to end declaratio Line 780  C--   dummy assignment to end declaratio
780    
781  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
782  CHPF$ INDEPENDENT  CHPF$ INDEPENDENT
783    CADJ INIT z121tape = common, Nr
784  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
785    
786        do i = 1, imt        do i = 1, imt
787    
788             k = 1
789    CADJ STORE v(i,k) = z121tape
790           v(i,Nrp1) = v(i,Nr)           v(i,Nrp1) = v(i,Nr)
791    
792           do k = 1, Nr           do k = 1, Nr
# Line 806  CHPF$ INDEPENDENT Line 801  CHPF$ INDEPENDENT
801           zflag  = p2 + KRi_range(1) * KRi_range(2)           zflag  = p2 + KRi_range(1) * KRi_range(2)
802           v(i,1) = v(i,1) / zflag           v(i,1) = v(i,1) / zflag
803    
 CADJ INIT z121tape = common, Nr  
804           do k = 2, Nr           do k = 2, Nr
805  CADJ STORE v(i,k), zwork = z121tape  CADJ STORE v(i,k), zwork = z121tape
806              km1 = k - 1              km1 = k - 1
# Line 829  CADJ STORE v(i,k), zwork = z121tape Line 823  CADJ STORE v(i,k), zwork = z121tape
823    
824  c*************************************************************************  c*************************************************************************
825    
826        subroutine smooth_horiz_rs (        subroutine kpp_smooth_horiz (
827       I     k, bi, bj,       I     k, bi, bj,
828       I     fld_in,       U     fld )
      O     fld_out )  
829    
830  c     Apply horizontal smoothing to RS 2-D array  c     Apply horizontal smoothing to KPP array
831    
832        IMPLICIT NONE        IMPLICIT NONE
833  #include "SIZE.h"  #include "SIZE.h"
834  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
835    
836  c     input  c     input
837  c     k, bi, bj : array indices  c     bi, bj : array indices
838  c     fld_in    : 2-D array to be smoothed  c     k      : vertical index used for masking
839        integer k, bi, bj        integer k, bi, bj
       _RS fld_in(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
840    
841  c     output  c     input/output
842  c     fld_out   : smoothed 2-D array  c     fld    : 2-D array to be smoothed
843        _RS fld_out(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _KPP_RL fld( ibot:itop, jbot:jtop )
844    
845  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
846    
847  c     local  c     local
848        integer i, j, im1, ip1, jm1, jp1        integer i, j, im1, ip1, jm1, jp1
849        _RS tempVar        _KPP_RL tempVar
850        _RS fld_tmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _KPP_RL fld_tmp( ibot:itop, jbot:jtop )
851    
852        integer    imin        , imax        integer    imin       , imax       , jmin       , jmax
853        parameter( imin=1-OLx+1, imax=sNx+OLx-1 )        parameter( imin=ibot+1, imax=itop-1, jmin=jbot+1, jmax=jtop-1 )
854    
855        integer    jmin        , jmax        _KPP_RL    p0    , p5    , p25     , p125      , p0625
       parameter( jmin=1-OLy+1, jmax=sNy+OLy-1 )  
   
       _RS        p0    , p5    , p25     , p125      , p0625  
856        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 )
857    
858        DO j = jmin, jmax        DO j = jmin, jmax
# Line 884  c     local Line 873  c     local
873       &                     pMask(ip1,jp1,k,bi,bj) )       &                     pMask(ip1,jp1,k,bi,bj) )
874              IF ( tempVar .GE. p25 ) THEN              IF ( tempVar .GE. p25 ) THEN
875                 fld_tmp(i,j) = (                 fld_tmp(i,j) = (
876       &              p25  * fld_in(i  ,j  )*pMask(i  ,j  ,k,bi,bj) +       &              p25  * fld(i  ,j  )*pMask(i  ,j  ,k,bi,bj) +
877       &              p125 *(fld_in(im1,j  )*pMask(im1,j  ,k,bi,bj) +       &              p125 *(fld(im1,j  )*pMask(im1,j  ,k,bi,bj) +
878       &                     fld_in(ip1,j  )*pMask(ip1,j  ,k,bi,bj) +       &                     fld(ip1,j  )*pMask(ip1,j  ,k,bi,bj) +
879       &                     fld_in(i  ,jm1)*pMask(i  ,jm1,k,bi,bj) +       &                     fld(i  ,jm1)*pMask(i  ,jm1,k,bi,bj) +
880       &                     fld_in(i  ,jp1)*pMask(i  ,jp1,k,bi,bj))+       &                     fld(i  ,jp1)*pMask(i  ,jp1,k,bi,bj))+
881       &              p0625*(fld_in(im1,jm1)*pMask(im1,jm1,k,bi,bj) +       &              p0625*(fld(im1,jm1)*pMask(im1,jm1,k,bi,bj) +
882       &                     fld_in(im1,jp1)*pMask(im1,jp1,k,bi,bj) +       &                     fld(im1,jp1)*pMask(im1,jp1,k,bi,bj) +
883       &                     fld_in(ip1,jm1)*pMask(ip1,jm1,k,bi,bj) +       &                     fld(ip1,jm1)*pMask(ip1,jm1,k,bi,bj) +
884       &                     fld_in(ip1,jp1)*pMask(ip1,jp1,k,bi,bj)))       &                     fld(ip1,jp1)*pMask(ip1,jp1,k,bi,bj)))
885       &              / tempVar       &              / tempVar
886              ELSE              ELSE
887                 fld_tmp(i,j) = fld_in(i,j)                 fld_tmp(i,j) = fld(i,j)
888              ENDIF              ENDIF
889           ENDDO           ENDDO
890        ENDDO        ENDDO
# Line 903  c     local Line 892  c     local
892  c     transfer smoothed field to output array  c     transfer smoothed field to output array
893        DO j = jmin, jmax        DO j = jmin, jmax
894           DO i = imin, imax           DO i = imin, imax
895              fld_out(i,j) = fld_tmp(i,j)              fld(i,j) = fld_tmp(i,j)
896           ENDDO           ENDDO
897        ENDDO        ENDDO
898    
 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  
   
899  #endif /* ALLOW_KPP */  #endif /* ALLOW_KPP */
900    
901        return        return
# Line 932  c     set output array edges to input fi Line 903  c     set output array edges to input fi
903    
904  c*************************************************************************  c*************************************************************************
905    
906        subroutine smooth_horiz_rl (        subroutine smooth_horiz (
907       I     k, bi, bj,       I     k, bi, bj,
908       I     fld_in,       U     fld )
      O     fld_out )  
909    
910  c     Apply horizontal smoothing to RL 2-D array  c     Apply horizontal smoothing to global _RL 2-D array
911    
912        IMPLICIT NONE        IMPLICIT NONE
913  #include "SIZE.h"  #include "SIZE.h"
914  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
915    
916  c     input  c     input
917  c     k, bi, bj : array indices  c     bi, bj : array indices
918  c     fld_in    : 2-D array to be smoothed  c     k      : vertical index used for masking
919        integer k, bi, bj        integer k, bi, bj
       _RL fld_in(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
920    
921  c     output  c     input/output
922  c     fld_out   : smoothed 2-D array  c     fld    : 2-D array to be smoothed
923        _RL fld_out(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fld( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
924    
925  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
926    
927  c     local  c     local
928        integer i, j, im1, ip1, jm1, jp1        integer i, j, im1, ip1, jm1, jp1
929        _RL tempVar        _RL tempVar
930        _RL fld_tmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fld_tmp( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
931    
932        integer    imin        , imax        integer   imin      , imax          , jmin      , jmax
933        parameter( imin=1-OLx+1, imax=sNx+OLx-1 )        parameter(imin=2-OLx, imax=sNx+OLx-1, jmin=2-OLy, jmax=sNy+OLy-1)
934    
935        integer    jmin        , jmax        _RL        p0    , p5    , p25     , p125      , p0625
       parameter( jmin=1-OLy+1, jmax=sNy+OLy-1 )  
   
       _RS        p0    , p5    , p25     , p125      , p0625  
936        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 )
937    
938        DO j = jmin, jmax        DO j = jmin, jmax
# Line 987  c     local Line 953  c     local
953       &                     pMask(ip1,jp1,k,bi,bj) )       &                     pMask(ip1,jp1,k,bi,bj) )
954              IF ( tempVar .GE. p25 ) THEN              IF ( tempVar .GE. p25 ) THEN
955                 fld_tmp(i,j) = (                 fld_tmp(i,j) = (
956       &              p25  * fld_in(i  ,j  )*pMask(i  ,j  ,k,bi,bj) +       &              p25  * fld(i  ,j  )*pMask(i  ,j  ,k,bi,bj) +
957       &              p125 *(fld_in(im1,j  )*pMask(im1,j  ,k,bi,bj) +       &              p125 *(fld(im1,j  )*pMask(im1,j  ,k,bi,bj) +
958       &                     fld_in(ip1,j  )*pMask(ip1,j  ,k,bi,bj) +       &                     fld(ip1,j  )*pMask(ip1,j  ,k,bi,bj) +
959       &                     fld_in(i  ,jm1)*pMask(i  ,jm1,k,bi,bj) +       &                     fld(i  ,jm1)*pMask(i  ,jm1,k,bi,bj) +
960       &                     fld_in(i  ,jp1)*pMask(i  ,jp1,k,bi,bj))+       &                     fld(i  ,jp1)*pMask(i  ,jp1,k,bi,bj))+
961       &              p0625*(fld_in(im1,jm1)*pMask(im1,jm1,k,bi,bj) +       &              p0625*(fld(im1,jm1)*pMask(im1,jm1,k,bi,bj) +
962       &                     fld_in(im1,jp1)*pMask(im1,jp1,k,bi,bj) +       &                     fld(im1,jp1)*pMask(im1,jp1,k,bi,bj) +
963       &                     fld_in(ip1,jm1)*pMask(ip1,jm1,k,bi,bj) +       &                     fld(ip1,jm1)*pMask(ip1,jm1,k,bi,bj) +
964       &                     fld_in(ip1,jp1)*pMask(ip1,jp1,k,bi,bj)))       &                     fld(ip1,jp1)*pMask(ip1,jp1,k,bi,bj)))
965       &              / tempVar       &              / tempVar
966              ELSE              ELSE
967                 fld_tmp(i,j) = fld_in(i,j)                 fld_tmp(i,j) = fld(i,j)
968              ENDIF              ENDIF
969           ENDDO           ENDDO
970        ENDDO        ENDDO
# Line 1006  c     local Line 972  c     local
972  c     transfer smoothed field to output array  c     transfer smoothed field to output array
973        DO j = jmin, jmax        DO j = jmin, jmax
974           DO i = imin, imax           DO i = imin, imax
975              fld_out(i,j) = fld_tmp(i,j)              fld(i,j) = fld_tmp(i,j)
976           ENDDO           ENDDO
977        ENDDO        ENDDO
978    
 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  
   
979  #endif /* ALLOW_KPP */  #endif /* ALLOW_KPP */
980    
981        return        return
# Line 1037  c*************************************** Line 985  c***************************************
985    
986        subroutine blmix (        subroutine blmix (
987       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl
988       O     , dkm1, blmc, ghat, sigma       O     , dkm1, blmc, ghat, sigma, ikey
989       &     )       &     )
990    
991  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 1009  c     stable(imt)                 = 1 in
1009  c     casea (imt)                 = 1 in case A  c     casea (imt)                 = 1 in case A
1010  c     diffus(imt,0:Nrp1,mdiff)    vertical diffusivities              (m^2/s)  c     diffus(imt,0:Nrp1,mdiff)    vertical diffusivities              (m^2/s)
1011  c     kbl(imt)                    -1 of first grid level below hbl  c     kbl(imt)                    -1 of first grid level below hbl
1012        _RS ustar (imt)        _KPP_RL ustar (imt)
1013        _RS bfsfc (imt)        _KPP_RL bfsfc (imt)
1014        _RS hbl   (imt)        _KPP_RL hbl   (imt)
1015        _RS stable(imt)        _KPP_RL stable(imt)
1016        _RS casea (imt)        _KPP_RL casea (imt)
1017        _RS diffus(imt,0:Nrp1,mdiff)        _KPP_RL diffus(imt,0:Nrp1,mdiff)
1018        integer kbl(imt)        integer kbl(imt)
1019    
1020  c output  c output
# Line 1074  c     dkm1 (imt,mdiff)            bounda Line 1022  c     dkm1 (imt,mdiff)            bounda
1022  c     blmc (imt,Nr,mdiff)         boundary layer mixing coefficients  (m^2/s)  c     blmc (imt,Nr,mdiff)         boundary layer mixing coefficients  (m^2/s)
1023  c     ghat (imt,Nr)               nonlocal scalar transport  c     ghat (imt,Nr)               nonlocal scalar transport
1024  c     sigma(imt)                  normalized depth (d / hbl)  c     sigma(imt)                  normalized depth (d / hbl)
1025        _RS dkm1 (imt,mdiff)        _KPP_RL dkm1 (imt,mdiff)
1026        _RS blmc (imt,Nr,mdiff)        _KPP_RL blmc (imt,Nr,mdiff)
1027        _RS ghat (imt,Nr)        _KPP_RL ghat (imt,Nr)
1028        _RS sigma(imt)        _KPP_RL sigma(imt)
1029          integer ikey
1030    
1031  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1032    
# Line 1085  c  local Line 1034  c  local
1034  c     gat1*(imt)                 shape function at sigma = 1  c     gat1*(imt)                 shape function at sigma = 1
1035  c     dat1*(imt)                 derivative of shape function at sigma = 1  c     dat1*(imt)                 derivative of shape function at sigma = 1
1036  c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)  c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)
1037        _RS gat1m(imt), gat1s(imt), gat1t(imt)        _KPP_RL gat1m(imt), gat1s(imt), gat1t(imt)
1038        _RS dat1m(imt), dat1s(imt), dat1t(imt)        _KPP_RL dat1m(imt), dat1s(imt), dat1t(imt)
1039        _RS ws(imt), wm(imt)        _KPP_RL ws(imt), wm(imt)
1040        integer i, kn, ki        integer i, kn, ki
1041        _RS     R, dvdzup, dvdzdn, viscp        _KPP_RL R, dvdzup, dvdzdn, viscp
1042        _RS     difsp, diftp, visch, difsh, difth        _KPP_RL difsp, diftp, visch, difsh, difth
1043        _RS     f1, sig, a1, a2, a3, delhat        _KPP_RL f1, sig, a1, a2, a3, delhat
1044        _RS     Gm, Gs, Gt        _KPP_RL Gm, Gs, Gt
1045        _RL     dum        _KPP_RL tempVar
1046    
1047        _RS        p0    , eins        _KPP_RL    p0    , eins
1048        parameter (p0=0.0, eins=1.0)        parameter (p0=0.0, eins=1.0)
1049    
1050  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 1111  c--------------------------------------- Line 1060  c---------------------------------------
1060       O        wm, ws )       O        wm, ws )
1061    
1062        do i = 1, imt        do i = 1, imt
1063             wm(i) = sign(eins,wm(i))*max(phepsi,abs(wm(i)))
1064             ws(i) = sign(eins,ws(i))*max(phepsi,abs(ws(i)))
1065          end do
1066    CADJ STORE wm = comlev1_kpp, key = ikey
1067    CADJ STORE ws = comlev1_kpp, key = ikey
1068    
1069          do i = 1, imt
1070    
1071           kn = int(caseA(i)+phepsi) *(kbl(i) -1) +           kn = int(caseA(i)+phepsi) *(kbl(i) -1) +
1072       $        (1 - int(caseA(i)+phepsi)) * kbl(i)       $        (1 - int(caseA(i)+phepsi)) * kbl(i)
# Line 1140  c--------------------------------------- Line 1096  c---------------------------------------
1096           difsh  = diffus(i,kn,2) + difsp * delhat           difsh  = diffus(i,kn,2) + difsp * delhat
1097           difth  = diffus(i,kn,3) + diftp * delhat           difth  = diffus(i,kn,3) + diftp * delhat
1098    
 #ifdef KPP_TEST_DENOM  
 ctl replace (Important!!! not phepsi**4 !!!)  
1099           f1 = stable(i) * conc1 * bfsfc(i) /           f1 = stable(i) * conc1 * bfsfc(i) /
1100       &        max(ustar(i)**4,phepsi)       &        max(ustar(i)**4,phepsi)
          wm(i) = sign(eins,wm(i))*max(phepsi,abs(wm(i)))  
1101           gat1m(i) = visch / hbl(i) / wm(i)           gat1m(i) = visch / hbl(i) / wm(i)
1102           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  
1103           dat1m(i) = min(dat1m(i),p0)           dat1m(i) = min(dat1m(i),p0)
1104    
 #ifdef KPP_TEST_DENOM  
          ws(i) = sign(eins,ws(i))*max(phepsi,abs(ws(i)))  
1105           gat1s(i) = difsh  / hbl(i) / ws(i)           gat1s(i) = difsh  / hbl(i) / ws(i)
1106           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  
1107           dat1s(i) = min(dat1s(i),p0)           dat1s(i) = min(dat1s(i),p0)
1108    
 #ifdef KPP_TEST_DENOM  
1109           gat1t(i) = difth /  hbl(i) / ws(i)           gat1t(i) = difth /  hbl(i) / ws(i)
1110           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  
1111           dat1t(i) = min(dat1t(i),p0)           dat1t(i) = min(dat1t(i),p0)
1112    
1113        end do        end do
# Line 1214  c--------------------------------------- Line 1151  c---------------------------------------
1151  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1152  c     nonlocal transport term = ghat * <ws>o  c     nonlocal transport term = ghat * <ws>o
1153  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1154  #ifdef KPP_TEST_DENOM  
1155              dum = ws(i) * hbl(i)              tempVar = ws(i) * hbl(i)
1156  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  
1157    
1158           end do           end do
1159        end do        end do
# Line 1280  c     hbl(imt)                  boundary Line 1211  c     hbl(imt)                  boundary
1211  c     kbl(imt)                  grid above hbl  c     kbl(imt)                  grid above hbl
1212  c     diffus(imt,0:Nrp1,mdiff) vertical diffusivities           (m^2/s)  c     diffus(imt,0:Nrp1,mdiff) vertical diffusivities           (m^2/s)
1213  c     casea(imt)                = 1 in caseA, = 0 in case B  c     casea(imt)                = 1 in caseA, = 0 in case B
1214        _RS     dkm1  (imt,mdiff)        _KPP_RL dkm1  (imt,mdiff)
1215        _RS     hbl   (imt)        _KPP_RL hbl   (imt)
1216        integer kbl   (imt)        integer kbl   (imt)
1217        _RS     diffus(imt,0:Nrp1,mdiff)        _KPP_RL diffus(imt,0:Nrp1,mdiff)
1218        _RS     casea (imt)        _KPP_RL casea (imt)
1219    
1220  c input/output  c input/output
1221  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)
1222        _RS     ghat (imt,Nr)        _KPP_RL ghat (imt,Nr)
1223    
1224  c output  c output
1225  c     enhanced bound. layer mixing coeff.  c     enhanced bound. layer mixing coeff.
1226        _RS     blmc  (imt,Nr,mdiff)        _KPP_RL blmc  (imt,Nr,mdiff)
1227    
1228  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1229    
1230  c local  c local
1231  c     fraction hbl lies beteen zgrid neighbors  c     fraction hbl lies beteen zgrid neighbors
1232        _RS delta        _KPP_RL delta
1233        integer ki, i, md        integer ki, i, md
1234        _RS dkmp5, dstar        _KPP_RL dkmp5, dstar
1235    
1236        do i = 1, imt        do i = 1, imt
1237           ki = kbl(i)-1           ki = kbl(i)-1
# Line 1363  c--------------------------------------- Line 1294  c---------------------------------------
1294    
1295  c-------------- Routine arguments -----------------------------------------  c-------------- Routine arguments -----------------------------------------
1296        INTEGER bi, bj, myThid        INTEGER bi, bj, myThid
1297  #ifdef FRUGAL_KPP        _KPP_RL RHO1   ( ibot:itop, jbot:jtop       )
1298        _RS RHO1   (sNx,sNy)        _KPP_RL DBLOC  ( ibot:itop, jbot:jtop, Nr   )
1299        _RS DBLOC  (sNx,sNy,Nr)        _KPP_RL DBSFC  ( ibot:itop, jbot:jtop, Nr   )
1300        _RS DBSFC  (sNx,sNy,Nr)        _KPP_RL TTALPHA( ibot:itop, jbot:jtop, Nrp1 )
1301        _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  
1302    
1303  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1304    
# Line 1399  c     work1, work2 - work arrays for hol Line 1322  c     work1, work2 - work arrays for hol
1322  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
1323    
1324        call FIND_RHO(        call FIND_RHO(
1325  #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  
1326       O     WORK1,       O     WORK1,
1327       I     myThid )       I     myThid )
1328    
1329        call FIND_ALPHA(        call FIND_ALPHA(
1330  #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  
1331       O     WORK2 )       O     WORK2 )
1332    
1333        call FIND_BETA(        call FIND_BETA(
1334  #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  
1335       O     WORK3 )       O     WORK3 )
1336    
1337  #ifdef FRUGAL_KPP        DO J = jbot, jtop
1338        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  
1339              RHO1(I,J)      = WORK1(I,J) + rhonil              RHO1(I,J)      = WORK1(I,J) + rhonil
1340              TTALPHA(I,J,1) = WORK2(I,J)              TTALPHA(I,J,1) = WORK2(I,J)
1341              SSBETA(I,J,1)  = WORK3(I,J)              SSBETA(I,J,1)  = WORK3(I,J)
# Line 1443  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO Line 1349  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO
1349        DO K = 2, Nr        DO K = 2, Nr
1350    
1351           call FIND_RHO(           call FIND_RHO(
1352  #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  
1353       O        RHOK,       O        RHOK,
1354       I        myThid )       I        myThid )
1355    
1356           call FIND_RHO(           call FIND_RHO(
1357  #ifdef FRUGAL_KPP       I        bi, bj, ibot, itop, jbot, jtop, K-1, K, eosType,
      I        bi, bj, 1, sNx, 1, sNy,  K-1, K, eosType,  
 #else  
      I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy,  K-1, K, eosType,  
 #endif  
1358       O        RHOKM1,       O        RHOKM1,
1359       I        myThid )       I        myThid )
1360    
1361           call FIND_RHO(           call FIND_RHO(
1362  #ifdef FRUGAL_KPP       I        bi, bj, ibot, itop, jbot, jtop, 1, K, eosType,
      I        bi, bj, 1, sNx, 1, sNy,  1,   K, eosType,  
 #else  
      I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy,  1,   K, eosType,  
 #endif  
1363       O        RHO1K,       O        RHO1K,
1364       I        myThid )       I        myThid )
1365    
1366           call FIND_ALPHA(           call FIND_ALPHA(
1367  #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  
1368       O        WORK1 )       O        WORK1 )
1369    
1370           call FIND_BETA(           call FIND_BETA(
1371  #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  
1372       O        WORK2 )       O        WORK2 )
1373    
1374  #ifdef FRUGAL_KPP           DO J = jbot, jtop
1375           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  
1376                 TTALPHA(I,J,K) = WORK1 (I,J)                 TTALPHA(I,J,K) = WORK1 (I,J)
1377                 SSBETA(I,J,K)  = WORK2 (I,J)                 SSBETA(I,J,K)  = WORK2 (I,J)
1378                 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 1385  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO
1385        END DO        END DO
1386    
1387  c     compute arrays for K = Nrp1  c     compute arrays for K = Nrp1
1388  #ifdef FRUGAL_KPP        DO J = jbot, jtop
1389        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  
1390              TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)              TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)
1391              SSBETA(I,J,Nrp1)  = SSBETA(I,J,Nr)              SSBETA(I,J,Nrp1)  = SSBETA(I,J,Nr)
1392              DBLOC(I,J,Nr)     = 0.              DBLOC(I,J,Nr)     = 0.
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22