/[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.22 by heimbach, Thu Feb 10 05:41:44 2005 UTC revision 1.35 by jmc, Fri Oct 19 19:11:17 2007 UTC
# Line 11  C--  o BLDEPTH     - Determine oceanic p Line 11  C--  o BLDEPTH     - Determine oceanic p
11  C--  o WSCALE      - Compute turbulent velocity scales.  C--  o WSCALE      - Compute turbulent velocity scales.
12  C--  o RI_IWMIX    - Compute interior viscosity diffusivity coefficients.  C--  o RI_IWMIX    - Compute interior viscosity diffusivity coefficients.
13  C--  o Z121        - Apply 121 vertical smoothing.  C--  o Z121        - Apply 121 vertical smoothing.
14  C--  o KPP_SMOOTH_HORIZ - Apply horizontal smoothing to KPP array.  C--  o SMOOTH_HORIZ- Apply horizontal smoothing to global array.
 C--  o SMOOTH_HORIZ     - Apply horizontal smoothing to global array.  
15  C--  o BLMIX       - Boundary layer mixing coefficients.  C--  o BLMIX       - Boundary layer mixing coefficients.
16  C--  o ENHANCE     - Enhance diffusivity at boundary layer interface.  C--  o ENHANCE     - Enhance diffusivity at boundary layer interface.
17  C--  o STATEKPP    - Compute buoyancy-related input arrays.  C--  o STATEKPP    - Compute buoyancy-related input arrays.
# Line 20  C--  o STATEKPP    - Compute buoyancy-re Line 19  C--  o STATEKPP    - Compute buoyancy-re
19  c***********************************************************************  c***********************************************************************
20    
21        SUBROUTINE KPPMIX (        SUBROUTINE KPPMIX (
22       I       mytime, mythid       I       kmtj, shsq, dvsq, ustar, msk
23       I     , kmtj, shsq, dvsq, ustar       I     , bo, bosol
24       I     , bo, bosol, dbloc, Ritop, coriol  #ifdef ALLOW_SALT_PLUME
25         I     , boplume,SPDepth
26    #endif /* ALLOW_SALT_PLUME */
27         I     , dbloc, Ritop, coriol
28         I     , diffusKzS, diffusKzT
29       I     , ikppkey       I     , ikppkey
30       O     , diffus       O     , diffus
31       U     , ghat       U     , ghat
32       O     , hbl )       O     , hbl
33         I     , myTime, myIter, myThid )
34    
35  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
36  c  c
# Line 47  c--------------------------------------- Line 51  c---------------------------------------
51  #include "SIZE.h"  #include "SIZE.h"
52  #include "EEPARAMS.h"  #include "EEPARAMS.h"
53  #include "PARAMS.h"  #include "PARAMS.h"
 #include "DYNVARS.h"  
 #include "FFIELDS.h"  
54  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
55    
56  c input  c input
57  c     myTime          - current time in simulation  c   myTime :: Current time in simulation
58  c     myThid          - thread number for this instance of the routine  c   myIter :: Current iteration number in simulation
59  c     kmtj   (imt)    - number of vertical layers on this row  c   myThid :: My Thread Id. number
60  c     shsq   (imt,Nr) - (local velocity shear)^2                      ((m/s)^2)  c     kmtj   (imt)     - number of vertical layers on this row
61  c     dvsq   (imt,Nr) - (velocity shear re sfc)^2                     ((m/s)^2)  c     msk    (imt)     - surface mask (=1 if water, =0 otherwise)
62  c     ustar  (imt)    - surface friction velocity                         (m/s)  c     shsq   (imt,Nr)  - (local velocity shear)^2                     ((m/s)^2)
63  c     bo     (imt)    - surface turbulent buoy. forcing               (m^2/s^3)  c     dvsq   (imt,Nr)  - (velocity shear re sfc)^2                    ((m/s)^2)
64  c     bosol  (imt)    - radiative buoyancy forcing                    (m^2/s^3)  c     ustar  (imt)     - surface friction velocity                        (m/s)
65  c     dbloc  (imt,Nr) - local delta buoyancy across interfaces          (m/s^2)  c     bo     (imt)     - surface turbulent buoy. forcing              (m^2/s^3)
66  c     dblocSm(imt,Nr) - horizontally smoothed dbloc                     (m/s^2)  c     bosol  (imt)     - radiative buoyancy forcing                   (m^2/s^3)
67  c                         stored in ghat to save space  c     boplume(imt)     - haline buoyancy forcing                      (m^2/s^3)
68  c     Ritop  (imt,Nr) - numerator of bulk Richardson Number  c     dbloc  (imt,Nr)  - local delta buoyancy across interfaces         (m/s^2)
69  c                         (zref-z) * delta buoyancy w.r.t. surface    ((m/s)^2)  c     dblocSm(imt,Nr)  - horizontally smoothed dbloc                    (m/s^2)
70  c     coriol (imt)    - Coriolis parameter                                (1/s)  c                          stored in ghat to save space
71    c     Ritop  (imt,Nr)  - numerator of bulk Richardson Number
72    c                          (zref-z) * delta buoyancy w.r.t. surface   ((m/s)^2)
73    c     coriol (imt)     - Coriolis parameter                               (1/s)
74    c     diffusKzS(imt,Nr)- background vertical diffusivity for scalars    (m^2/s)
75    c     diffusKzT(imt,Nr)- background vertical diffusivity for theta      (m^2/s)
76  c     note: there is a conversion from 2-D to 1-D for input output variables,  c     note: there is a conversion from 2-D to 1-D for input output variables,
77  c           e.g., hbl(sNx,sNy) -> hbl(imt),  c           e.g., hbl(sNx,sNy) -> hbl(imt),
78  c           where hbl(i,j) -> hbl((j-1)*sNx+i)  c           where hbl(i,j) -> hbl((j-1)*sNx+i)
79    
80        _RL     mytime        _RL     myTime
81        integer mythid        integer myIter
82        integer kmtj  (imt   )        integer myThid
83        _KPP_RL shsq  (imt,Nr)        integer kmtj (imt   )
84        _KPP_RL dvsq  (imt,Nr)        _RL shsq     (imt,Nr)
85        _KPP_RL ustar (imt   )        _RL dvsq     (imt,Nr)
86        _KPP_RL bo    (imt   )        _RL ustar    (imt   )
87        _KPP_RL bosol (imt   )        _RL bo       (imt   )
88        _KPP_RL dbloc (imt,Nr)        _RL bosol    (imt   )
89        _KPP_RL Ritop (imt,Nr)  #ifdef ALLOW_SALT_PLUME
90        _KPP_RL coriol(imt   )        _RL boplume  (imt   )
91          _RL SPDepth  (imt   )
92    #endif /* ALLOW_SALT_PLUME */
93          _RL dbloc    (imt,Nr)
94          _RL Ritop    (imt,Nr)
95          _RL coriol   (imt   )
96          _RS msk      (imt   )
97          _RL diffusKzS(imt,Nr)
98          _RL diffusKzT(imt,Nr)
99    
100        integer ikppkey        integer ikppkey
101    
# Line 91  c     diffus (imt,3)  - vertical tempera Line 106  c     diffus (imt,3)  - vertical tempera
106  c     ghat   (imt)    - nonlocal transport coefficient                  (s/m^2)  c     ghat   (imt)    - nonlocal transport coefficient                  (s/m^2)
107  c     hbl    (imt)    - mixing layer depth                                  (m)  c     hbl    (imt)    - mixing layer depth                                  (m)
108    
109        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL diffus(imt,0:Nrp1,mdiff)
110        _KPP_RL ghat  (imt,Nr)        _RL ghat  (imt,Nr)
111        _KPP_RL hbl   (imt)        _RL hbl   (imt)
112    
113  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
114    
# Line 107  c     blmc   (imt,Nr,mdiff) - boundary l Line 122  c     blmc   (imt,Nr,mdiff) - boundary l
122  c     sigma  (imt         ) - normalized depth (d / hbl)  c     sigma  (imt         ) - normalized depth (d / hbl)
123  c     Rib    (imt,Nr      ) - bulk Richardson number  c     Rib    (imt,Nr      ) - bulk Richardson number
124    
125        integer kbl   (imt         )        integer kbl(imt         )
126        _KPP_RL bfsfc (imt         )        _RL bfsfc  (imt         )
127        _KPP_RL casea (imt         )        _RL casea  (imt         )
128        _KPP_RL stable(imt         )        _RL stable (imt         )
129        _KPP_RL dkm1  (imt,   mdiff)        _RL dkm1   (imt,   mdiff)
130        _KPP_RL blmc  (imt,Nr,mdiff)        _RL blmc   (imt,Nr,mdiff)
131        _KPP_RL sigma (imt         )        _RL sigma  (imt         )
132        _KPP_RL Rib   (imt,Nr      )        _RL Rib    (imt,Nr      )
133    
134        integer i, k, md        integer i, k, md
135    
# Line 132  CADJ STORE dbloc = comlev1_kpp, key = ik Line 147  CADJ STORE dbloc = comlev1_kpp, key = ik
147  cph)  cph)
148        call Ri_iwmix (        call Ri_iwmix (
149       I       kmtj, shsq, dbloc, ghat       I       kmtj, shsq, dbloc, ghat
150         I     , diffusKzS, diffusKzT
151       I     , ikppkey       I     , ikppkey
152       O     , diffus )       O     , diffus, myThid )
153    
154  cph(  cph(
155  cph these storings avoid recomp. of Ri_iwmix  cph these storings avoid recomp. of Ri_iwmix
# Line 150  c for blmix Line 166  c for blmix
166  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
167    
168        do md = 1, mdiff        do md = 1, mdiff
169           do k=1,Nrp1
170           do i = 1,imt           do i = 1,imt
171              do k=kmtj(i),Nrp1               if(k.ge.kmtj(i))  diffus(i,k,md) = 0.0
                diffus(i,k,md) = 0.0  
172              end do              end do
173           end do           end do
174        end do        end do
# Line 164  c diagnose the new boundary layer depth Line 180  c diagnose the new boundary layer depth
180  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
181    
182        call bldepth (        call bldepth (
183       I       mytime, mythid       I       kmtj
184       I     , kmtj       I     , dvsq, dbloc, Ritop, ustar, bo, bosol
185       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol  #ifdef ALLOW_SALT_PLUME
186         I     , boplume,SPDepth
187    #endif /* ALLOW_SALT_PLUME */
188         I     , coriol
189       I     , ikppkey       I     , ikppkey
190       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma
191       &     )       I     , myTime, myIter, myThid )
192    
193  CADJ STORE hbl,bfsfc,stable,casea,kbl = comlev1_kpp, key = ikppkey  CADJ STORE hbl,bfsfc,stable,casea,kbl = comlev1_kpp, key = ikppkey
194    
# Line 180  c--------------------------------------- Line 199  c---------------------------------------
199        call blmix (        call blmix (
200       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl
201       O     , dkm1, blmc, ghat, sigma, ikppkey       O     , dkm1, blmc, ghat, sigma, ikppkey
202       &     )       I     , myThid )
203  cph(  cph(
204  CADJ STORE dkm1,blmc,ghat = comlev1_kpp, key = ikppkey  CADJ STORE dkm1,blmc,ghat = comlev1_kpp, key = ikppkey
205  CADJ STORE hbl, kbl, diffus, casea = comlev1_kpp, key = ikppkey  CADJ STORE hbl, kbl, diffus, casea = comlev1_kpp, key = ikppkey
# Line 193  c--------------------------------------- Line 212  c---------------------------------------
212        call enhance (        call enhance (
213       I       dkm1, hbl, kbl, diffus, casea       I       dkm1, hbl, kbl, diffus, casea
214       U     , ghat       U     , ghat
215       O     , blmc )       O     , blmc
216         I     , myThid )
217    
218  cph(  cph(
219  cph avoids recomp. of enhance  cph avoids recomp. of enhance
# Line 209  c--------------------------------------- Line 229  c---------------------------------------
229        do k = 1, Nr        do k = 1, Nr
230           do i = 1, imt           do i = 1, imt
231              if (k .lt. kbl(i)) then              if (k .lt. kbl(i)) then
232    #ifdef ALLOW_SHELFICE
233    C     when there is shelfice on top (msk(i)=0), reset the boundary layer
234    C     mixing coefficients blmc to pure Ri-number based mixing
235                   blmc(i,k,1) = max ( blmc(i,k,1)*msk(i),
236         &              diffus(i,k,1) )
237                   blmc(i,k,2) = max ( blmc(i,k,2)*msk(i),
238         &              diffus(i,k,2) )
239                   blmc(i,k,3) = max ( blmc(i,k,3)*msk(i),
240         &              diffus(i,k,3) )
241    #endif /* not ALLOW_SHELFICE */
242                 diffus(i,k,1) = max ( blmc(i,k,1), viscAr  )                 diffus(i,k,1) = max ( blmc(i,k,1), viscAr  )
243                 diffus(i,k,2) = max ( blmc(i,k,2), diffKrNrS(k) )                 diffus(i,k,2) = max ( blmc(i,k,2), diffusKzS(i,Nr) )
244                 diffus(i,k,3) = max ( blmc(i,k,3), diffKrNrT(k) )                 diffus(i,k,3) = max ( blmc(i,k,3), diffusKzT(i,Nr) )
245              else              else
246                 ghat(i,k) = 0.                 ghat(i,k) = 0.
247              endif              endif
# Line 226  c--------------------------------------- Line 256  c---------------------------------------
256  c*************************************************************************  c*************************************************************************
257    
258        subroutine bldepth (        subroutine bldepth (
259       I       mytime, mythid       I       kmtj
260       I     , kmtj       I     , dvsq, dbloc, Ritop, ustar, bo, bosol
261       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol  #ifdef ALLOW_SALT_PLUME
262         I     , boplume,SPDepth
263    #endif /* ALLOW_SALT_PLUME */
264         I     , coriol
265       I     , ikppkey       I     , ikppkey
266       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma
267       &     )       I     , myTime, myIter, myThid )
268    
269  c     the oceanic planetary boundary layer depth, hbl, is determined as  c     the oceanic planetary boundary layer depth, hbl, is determined as
270  c     the shallowest depth where the bulk Richardson number is  c     the shallowest depth where the bulk Richardson number is
# Line 257  c Line 290  c
290        IMPLICIT NONE        IMPLICIT NONE
291    
292  #include "SIZE.h"  #include "SIZE.h"
 #include "EEPARAMS.h"  
 #include "PARAMS.h"  
293  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
 #include "FFIELDS.h"  
294    
295  c input  c input
296  c------  c------
297  c myTime    : current time in simulation  c   myTime :: Current time in simulation
298  c myThid    : thread number for this instance of the routine  c   myIter :: Current iteration number in simulation
299    c   myThid :: My Thread Id. number
300  c kmtj      : number of vertical layers  c kmtj      : number of vertical layers
301  c dvsq      : (velocity shear re sfc)^2             ((m/s)^2)  c dvsq      : (velocity shear re sfc)^2             ((m/s)^2)
302  c dbloc     : local delta buoyancy across interfaces  (m/s^2)  c dbloc     : local delta buoyancy across interfaces  (m/s^2)
# Line 275  c             buoyancy with respect to s Line 306  c             buoyancy with respect to s
306  c ustar     : surface friction velocity                 (m/s)  c ustar     : surface friction velocity                 (m/s)
307  c bo        : surface turbulent buoyancy forcing    (m^2/s^3)  c bo        : surface turbulent buoyancy forcing    (m^2/s^3)
308  c bosol     : radiative buoyancy forcing            (m^2/s^3)  c bosol     : radiative buoyancy forcing            (m^2/s^3)
309    c boplume   : haline buoyancy forcing               (m^2/s^3)
310  c coriol    : Coriolis parameter                        (1/s)  c coriol    : Coriolis parameter                        (1/s)
311        _RL     mytime        _RL     myTime
312        integer mythid        integer myIter
313          integer myThid
314        integer kmtj(imt)        integer kmtj(imt)
315        _KPP_RL dvsq  (imt,Nr)        _RL dvsq    (imt,Nr)
316        _KPP_RL dbloc (imt,Nr)        _RL dbloc   (imt,Nr)
317        _KPP_RL Ritop (imt,Nr)        _RL Ritop   (imt,Nr)
318        _KPP_RL ustar (imt)        _RL ustar   (imt)
319        _KPP_RL bo    (imt)        _RL bo      (imt)
320        _KPP_RL bosol (imt)        _RL bosol   (imt)
321        _KPP_RL coriol(imt)        _RL coriol  (imt)
322        integer ikppkey, kkppkey        integer ikppkey, kkppkey
323    #ifdef ALLOW_SALT_PLUME
324          _RL boplume (imt)
325          _RL SPDepth (imt)
326    #endif /* ALLOW_SALT_PLUME */
327    
328  c  output  c  output
329  c--------  c--------
# Line 297  c casea     : =1 in case A, =0 in case B Line 334  c casea     : =1 in case A, =0 in case B
334  c kbl       : -1 of first grid level below hbl  c kbl       : -1 of first grid level below hbl
335  c Rib       : Bulk Richardson number  c Rib       : Bulk Richardson number
336  c sigma     : normalized depth (d/hbl)  c sigma     : normalized depth (d/hbl)
337        _KPP_RL hbl   (imt)        _RL hbl    (imt)
338        _KPP_RL bfsfc (imt)        _RL bfsfc  (imt)
339        _KPP_RL stable(imt)        _RL stable (imt)
340        _KPP_RL casea (imt)        _RL casea  (imt)
341        integer kbl   (imt)        integer kbl(imt)
342        _KPP_RL Rib   (imt,Nr)        _RL Rib    (imt,Nr)
343        _KPP_RL sigma (imt)        _RL sigma  (imt)
344    
345  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
346    
347  c  local  c  local
348  c-------  c-------
349  c wm, ws    : turbulent velocity scales         (m/s)  c wm, ws    : turbulent velocity scales         (m/s)
350        _KPP_RL wm(imt), ws(imt)        _RL wm(imt), ws(imt)
351        _RL     worka(imt)        _RL worka(imt)
352          _RL bvsq, vtsq, hekman, hmonob, hlimit, tempVar1, tempVar2
       _KPP_RL bvsq, vtsq, hekman, hmonob, hlimit, tempVar1, tempVar2  
353        integer i, kl        integer i, kl
354    
355        _KPP_RL     p5    , eins        _RL         p5    , eins
356        parameter ( p5=0.5, eins=1.0 )        parameter ( p5=0.5, eins=1.0 )
357        _RL         minusone        _RL         minusone
358        parameter ( minusone=-1.0 )        parameter ( minusone=-1.0 )
# Line 351  c     compute bfsfc = sw fraction at hbf Line 387  c     compute bfsfc = sw fraction at hbf
387           end do           end do
388  CADJ store worka = comlev1_kpp_k, key = kkppkey  CADJ store worka = comlev1_kpp_k, key = kkppkey
389           call SWFRAC(           call SWFRAC(
390       I        imt, hbf,       I       imt, hbf,
391       I        mytime, mythid,       U       worka,
392       U        worka )       I       myTime, myIter, myThid )
393  CADJ store worka = comlev1_kpp_k, key = kkppkey  CADJ store worka = comlev1_kpp_k, key = kkppkey
394    
   
395           do i = 1, imt           do i = 1, imt
396    
397  c     use caseA as temporary array  c     use caseA as temporary array
# Line 366  c     use caseA as temporary array Line 401  c     use caseA as temporary array
401  c     compute bfsfc= Bo + radiative contribution down to hbf * hbl  c     compute bfsfc= Bo + radiative contribution down to hbf * hbl
402    
403              bfsfc(i) = bo(i) + bosol(i)*(1. - worka(i))              bfsfc(i) = bo(i) + bosol(i)*(1. - worka(i))
             stable(i) = p5 + sign(p5,bfsfc(i))  
             sigma(i) = stable(i) + (1. - stable(i)) * epsilon  
404    
405           end do           end do
406    #ifdef ALLOW_SALT_PLUME
407    c     compute bfsfc = plume fraction at hbf * zgrid
408             do i = 1, imt
409                worka(i) = zgrid(kl)
410             enddo
411             call PLUMEFRAC(
412         I       imt, hbf,SPDepth,
413         U       worka,
414         I       myTime, myIter, myThid)
415             do i = 1, imt
416                bfsfc(i) = bfsfc(i) + boplume(i)*(1. - worka(i))
417             enddo
418    #endif /* ALLOW_SALT_PLUME */
419    
420             do i = 1, imt
421                stable(i) = p5 + sign(p5,bfsfc(i))
422                sigma(i) = stable(i) + (1. - stable(i)) * epsilon
423             enddo
424    
425  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
426  c     compute velocity scales at sigma, for hbl= caseA = -zgrid(kl)  c     compute velocity scales at sigma, for hbl= caseA = -zgrid(kl)
# Line 377  c--------------------------------------- Line 428  c---------------------------------------
428    
429           call wscale (           call wscale (
430       I        sigma, casea, ustar, bfsfc,       I        sigma, casea, ustar, bfsfc,
431       O        wm, ws )       O        wm, ws, myThid )
432  CADJ store ws = comlev1_kpp_k, key = kkppkey  CADJ store ws = comlev1_kpp_k, key = kkppkey
433    
434           do i = 1, imt           do i = 1, imt
# Line 454  c--------------------------------------- Line 505  c---------------------------------------
505  CADJ store worka = comlev1_kpp  CADJ store worka = comlev1_kpp
506  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
507        call SWFRAC(        call SWFRAC(
508       I     imt, minusone,       I       imt, minusone,
509       I     mytime, mythid,       U       worka,
510       U     worka )       I       myTime, myIter, myThid )
511  CADJ store worka = comlev1_kpp  CADJ store worka = comlev1_kpp
512  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
513    
514        do i = 1, imt        do i = 1, imt
515           bfsfc(i)  = bo(i) + bosol(i) * (1. - worka(i))           bfsfc(i)  = bo(i) + bosol(i) * (1. - worka(i))
516        end do        end do
517    
518    #ifdef ALLOW_SALT_PLUME
519          do i = 1, imt
520             worka(i) = hbl(i)
521          enddo
522          call PLUMEFRAC(
523         I       imt,minusone,SPDepth,
524         U       worka,
525         I       myTime, myIter, myThid )
526          do i = 1, imt
527             bfsfc(i) = bfsfc(i) + boplume(i) * (1. - worka(i))
528          enddo
529    #endif /* ALLOW_SALT_PLUME */
530  CADJ store bfsfc = comlev1_kpp  CADJ store bfsfc = comlev1_kpp
531  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
532    
# Line 527  CADJ store worka = comlev1_kpp Line 591  CADJ store worka = comlev1_kpp
591  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
592        call SWFRAC(        call SWFRAC(
593       I     imt, minusone,       I     imt, minusone,
594       I     mytime, mythid,       U     worka,
595       U     worka )       I     myTime, myIter, myThid )
596  CADJ store worka = comlev1_kpp  CADJ store worka = comlev1_kpp
597  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
598    
599        do i = 1, imt        do i = 1, imt
600           bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i))           bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i))
601        end do        end do
602    
603    #ifdef ALLOW_SALT_PLUME
604          do i = 1, imt
605             worka(i) = hbl(i)
606          enddo
607          call PLUMEFRAC(
608         I       imt,minusone,SPDepth,
609         U       worka,
610         I       myTime, myIter, myThid )
611          do i = 1, imt
612             bfsfc(i) = bfsfc(i) + boplume(i) * (1. - worka(i))
613          enddo
614    #endif /* ALLOW_SALT_PLUME */
615  CADJ store bfsfc = comlev1_kpp  CADJ store bfsfc = comlev1_kpp
616  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
617    
# Line 562  c*************************************** Line 639  c***************************************
639    
640        subroutine wscale (        subroutine wscale (
641       I     sigma, hbl, ustar, bfsfc,       I     sigma, hbl, ustar, bfsfc,
642       O     wm, ws )       O     wm, ws,
643         I     myThid )
644    
645  c     compute turbulent velocity scales.  c     compute turbulent velocity scales.
646  c     use a 2D-lookup table for wm and ws as functions of ustar and  c     use a 2D-lookup table for wm and ws as functions of ustar and
# Line 583  c sigma   : normalized depth (d/hbl) Line 661  c sigma   : normalized depth (d/hbl)
661  c hbl     : boundary layer depth (m)  c hbl     : boundary layer depth (m)
662  c ustar   : surface friction velocity         (m/s)  c ustar   : surface friction velocity         (m/s)
663  c bfsfc   : total surface buoyancy flux       (m^2/s^3)  c bfsfc   : total surface buoyancy flux       (m^2/s^3)
664        _KPP_RL sigma(imt)  c myThid  : thread number for this instance of the routine
665        _KPP_RL hbl  (imt)        integer myThid
666        _KPP_RL ustar(imt)        _RL sigma(imt)
667        _KPP_RL bfsfc(imt)        _RL hbl  (imt)
668          _RL ustar(imt)
669          _RL bfsfc(imt)
670    
671  c  output  c  output
672  c--------  c--------
673  c wm, ws  : turbulent velocity scales at sigma  c wm, ws  : turbulent velocity scales at sigma
674        _KPP_RL wm(imt), ws(imt)        _RL wm(imt), ws(imt)
675    
676  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
677    
678  c local  c local
679  c------  c------
680  c zehat   : = zeta *  ustar**3  c zehat   : = zeta *  ustar**3
681        _KPP_RL zehat        _RL zehat
682    
683        integer iz, izp1, ju, i, jup1        integer iz, izp1, ju, i, jup1
684        _KPP_RL udiff, zdiff, zfrac, ufrac, fzfrac, wam        _RL udiff, zdiff, zfrac, ufrac, fzfrac, wam
685        _KPP_RL wbm, was, wbs, u3, tempVar        _RL wbm, was, wbs, u3, tempVar
686    
687  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
688  c use lookup table for zehat < zmax only; otherwise use  c use lookup table for zehat < zmax only; otherwise use
# Line 657  c--------------------------------------- Line 737  c---------------------------------------
737  c*************************************************************************  c*************************************************************************
738    
739        subroutine Ri_iwmix (        subroutine Ri_iwmix (
740       I       kmtj, shsq, dbloc, dblocSm       I       kmtj, shsq, dbloc, dblocSm,
741       I     , ikppkey       I       diffusKzS, diffusKzT,
742       O     , diffus )       I       ikppkey,
743         O       diffus,
744         I       myThid )
745    
746  c     compute interior viscosity diffusivity coefficients due  c     compute interior viscosity diffusivity coefficients due
747  c     to shear instability (dependent on a local Richardson number),  c     to shear instability (dependent on a local Richardson number),
# Line 678  c     kmtj   (imt)         number of ver Line 760  c     kmtj   (imt)         number of ver
760  c     shsq   (imt,Nr)      (local velocity shear)^2               ((m/s)^2)  c     shsq   (imt,Nr)      (local velocity shear)^2               ((m/s)^2)
761  c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)  c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)
762  c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)  c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)
763        integer kmtj   (imt)  c     diffusKzS(imt,Nr)- background vertical diffusivity for scalars    (m^2/s)
764        _KPP_RL shsq   (imt,Nr)  c     diffusKzT(imt,Nr)- background vertical diffusivity for theta      (m^2/s)
765        _KPP_RL dbloc  (imt,Nr)  c     myThid :: My Thread Id. number
766        _KPP_RL dblocSm(imt,Nr)        integer kmtj (imt)
767          _RL shsq     (imt,Nr)
768          _RL dbloc    (imt,Nr)
769          _RL dblocSm  (imt,Nr)
770          _RL diffusKzS(imt,Nr)
771          _RL diffusKzT(imt,Nr)
772        integer ikppkey        integer ikppkey
773          integer myThid
774    
775  c output  c output
776  c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)  c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)
777  c     diffus(imt,0:Nrp1,2)  vertical scalar diffusivity             (m^2/s)  c     diffus(imt,0:Nrp1,2)  vertical scalar diffusivity             (m^2/s)
778  c     diffus(imt,0:Nrp1,3)  vertical temperature diffusivity        (m^2/s)  c     diffus(imt,0:Nrp1,3)  vertical temperature diffusivity        (m^2/s)
779        _KPP_RL diffus(imt,0:Nrp1,3)        _RL diffus(imt,0:Nrp1,3)
780    
781  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
782    
783  c local variables  c local variables
784  c     Rig                   local Richardson number  c     Rig                   local Richardson number
785  c     fRi, fcon             function of Rig  c     fRi, fcon             function of Rig
786        _KPP_RL Rig        _RL Rig
787        _KPP_RL fRi, fcon        _RL fRi, fcon
788        _KPP_RL ratio        _RL ratio
789        integer i, ki        integer i, ki
790        _KPP_RL c1, c0        _RL c1, c0
791    
792  #ifdef ALLOW_KPP_VERTICALLY_SMOOTH  #ifdef ALLOW_KPP_VERTICALLY_SMOOTH
793        integer mr        integer mr
# Line 755  CADJ store diffus(:,:,1) = kpp_ri_tape_m Line 843  CADJ store diffus(:,:,1) = kpp_ri_tape_m
843  CADJ &  , key=mr, shape=(/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2 /)  CADJ &  , key=mr, shape=(/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2 /)
844    
845           call z121 (           call z121 (
846       U     diffus(1,0,1))       U     diffus(1,0,1)
847         I     myThid )
848        end do        end do
849  #endif  #endif
850    
# Line 786  c    mixing due to internal waves, and s Line 875  c    mixing due to internal waves, and s
875  #ifndef EXCLUDE_KPP_SHEAR_MIX  #ifndef EXCLUDE_KPP_SHEAR_MIX
876              if ( .NOT. inAdMode ) then              if ( .NOT. inAdMode ) then
877                 diffus(i,ki,1) = viscAr + fcon * difmcon + fRi * difm0                 diffus(i,ki,1) = viscAr + fcon * difmcon + fRi * difm0
878                 diffus(i,ki,2) = diffKrNrS(ki)+ fcon*difscon + fRi*difs0                 diffus(i,ki,2) = diffusKzS(i,ki)+fcon*difscon+fRi*difs0
879                 diffus(i,ki,3) = diffKrNrT(ki)+ fcon*difscon + fRi*difs0                 diffus(i,ki,3) = diffusKzT(i,ki)+fcon*diftcon+fRi*dift0
880              else              else
881                 diffus(i,ki,1) = viscAr                 diffus(i,ki,1) = viscAr
882                 diffus(i,ki,2) = diffKrNrS(ki)                 diffus(i,ki,2) = diffusKzS(i,ki)
883                 diffus(i,ki,3) = diffKrNrT(ki)                 diffus(i,ki,3) = diffusKzT(i,ki)
884              endif              endif
885  #else  #else
886              diffus(i,ki,1) = viscAr              diffus(i,ki,1) = viscAr
887              diffus(i,ki,2) = diffKrNrS(ki)              diffus(i,ki,2) = diffusKzS(i,ki)
888              diffus(i,ki,3) = diffKrNrT(ki)              diffus(i,ki,3) = diffusKzT(i,ki)
889  #endif  #endif
890    
891           end do           end do
# Line 819  c         set surface values to 0.0 Line 908  c         set surface values to 0.0
908  c*************************************************************************  c*************************************************************************
909    
910        subroutine z121 (        subroutine z121 (
911       U     v )       U     v,
912         I     myThid )
913    
914  c     Apply 121 smoothing in k to 2-d array V(i,k=1,Nr)  c     Apply 121 smoothing in k to 2-d array V(i,k=1,Nr)
915  c     top (0) value is used as a dummy  c     top (0) value is used as a dummy
# Line 837  c     penetrative convection. Line 927  c     penetrative convection.
927  c input/output  c input/output
928  c-------------  c-------------
929  c v     : 2-D array to be smoothed in Nrp1 direction  c v     : 2-D array to be smoothed in Nrp1 direction
930        _KPP_RL v(imt,0:Nrp1)  c myThid: thread number for this instance of the routine
931          integer myThid
932          _RL v(imt,0:Nrp1)
933    
934  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
935    
936  c local  c local
937        _KPP_RL zwork, zflag        _RL zwork, zflag
938        _KPP_RL KRi_range(1:Nrp1)        _RL KRi_range(1:Nrp1)
939        integer i, k, km1, kp1        integer i, k, km1, kp1
940    
941        _KPP_RL     p0      , p25       , p5      , p2        _RL         p0      , p25       , p5      , p2
942        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 )
943    
944        KRi_range(Nrp1) = p0        KRi_range(Nrp1) = p0
# Line 900  CADJ STORE v(i,k), zwork = z121tape Line 992  CADJ STORE v(i,k), zwork = z121tape
992    
993  c*************************************************************************  c*************************************************************************
994    
       subroutine kpp_smooth_horiz (  
      I     k, bi, bj,  
      U     fld )  
   
 c     Apply horizontal smoothing to KPP array  
   
       IMPLICIT NONE  
 #include "SIZE.h"  
 #include "GRID.h"  
 #include "KPP_PARAMS.h"  
   
 c     input  
 c     bi, bj : array indices  
 c     k      : vertical index used for masking  
       integer k, bi, bj  
   
 c     input/output  
 c     fld    : 2-D array to be smoothed  
       _KPP_RL fld( ibot:itop, jbot:jtop )  
   
 #ifdef ALLOW_KPP  
   
 c     local  
       integer i, j, im1, ip1, jm1, jp1  
       _KPP_RL tempVar  
       _KPP_RL fld_tmp( ibot:itop, jbot:jtop )  
   
       integer    imin       , imax       , jmin       , jmax  
       parameter( imin=ibot+1, imax=itop-1, jmin=jbot+1, jmax=jtop-1 )  
   
       _KPP_RL    p0    , p5    , p25     , p125      , p0625  
       parameter( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 )  
   
       DO j = jmin, jmax  
          jm1 = j-1  
          jp1 = j+1  
          DO i = imin, imax  
             im1 = i-1  
             ip1 = i+1  
             tempVar =  
      &           p25   *   maskC(i  ,j  ,k,bi,bj)   +  
      &           p125  * ( maskC(im1,j  ,k,bi,bj)   +  
      &                     maskC(ip1,j  ,k,bi,bj)   +  
      &                     maskC(i  ,jm1,k,bi,bj)   +  
      &                     maskC(i  ,jp1,k,bi,bj) ) +  
      &           p0625 * ( maskC(im1,jm1,k,bi,bj)   +  
      &                     maskC(im1,jp1,k,bi,bj)   +  
      &                     maskC(ip1,jm1,k,bi,bj)   +  
      &                     maskC(ip1,jp1,k,bi,bj) )  
             IF ( tempVar .GE. p25 ) THEN  
                fld_tmp(i,j) = (  
      &              p25  * fld(i  ,j  )*maskC(i  ,j  ,k,bi,bj) +  
      &              p125 *(fld(im1,j  )*maskC(im1,j  ,k,bi,bj) +  
      &                     fld(ip1,j  )*maskC(ip1,j  ,k,bi,bj) +  
      &                     fld(i  ,jm1)*maskC(i  ,jm1,k,bi,bj) +  
      &                     fld(i  ,jp1)*maskC(i  ,jp1,k,bi,bj))+  
      &              p0625*(fld(im1,jm1)*maskC(im1,jm1,k,bi,bj) +  
      &                     fld(im1,jp1)*maskC(im1,jp1,k,bi,bj) +  
      &                     fld(ip1,jm1)*maskC(ip1,jm1,k,bi,bj) +  
      &                     fld(ip1,jp1)*maskC(ip1,jp1,k,bi,bj)))  
      &              / tempVar  
             ELSE  
                fld_tmp(i,j) = fld(i,j)  
             ENDIF  
          ENDDO  
       ENDDO  
   
 c     transfer smoothed field to output array  
       DO j = jmin, jmax  
          DO i = imin, imax  
             fld(i,j) = fld_tmp(i,j)  
          ENDDO  
       ENDDO  
   
 #endif /* ALLOW_KPP */  
   
       return  
       end  
   
 c*************************************************************************  
   
995        subroutine smooth_horiz (        subroutine smooth_horiz (
996       I     k, bi, bj,       I     k, bi, bj,
997       U     fld )       U     fld,
998         I     myThid )
999    
1000  c     Apply horizontal smoothing to global _RL 2-D array  c     Apply horizontal smoothing to global _RL 2-D array
1001    
# Line 995  c     Apply horizontal smoothing to glob Line 1007  c     Apply horizontal smoothing to glob
1007  c     input  c     input
1008  c     bi, bj : array indices  c     bi, bj : array indices
1009  c     k      : vertical index used for masking  c     k      : vertical index used for masking
1010    c     myThid : thread number for this instance of the routine
1011          INTEGER myThid
1012        integer k, bi, bj        integer k, bi, bj
1013    
1014  c     input/output  c     input/output
# Line 1065  c*************************************** Line 1079  c***************************************
1079        subroutine blmix (        subroutine blmix (
1080       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl
1081       O     , dkm1, blmc, ghat, sigma, ikppkey       O     , dkm1, blmc, ghat, sigma, ikppkey
1082       &     )       I     , myThid )
1083    
1084  c     mixing coefficients within boundary layer depend on surface  c     mixing coefficients within boundary layer depend on surface
1085  c     forcing and the magnitude and gradient of interior mixing below  c     forcing and the magnitude and gradient of interior mixing below
# Line 1087  c     hbl   (imt)                 bounda Line 1101  c     hbl   (imt)                 bounda
1101  c     stable(imt)                 = 1 in stable forcing  c     stable(imt)                 = 1 in stable forcing
1102  c     casea (imt)                 = 1 in case A  c     casea (imt)                 = 1 in case A
1103  c     diffus(imt,0:Nrp1,mdiff)    vertical diffusivities              (m^2/s)  c     diffus(imt,0:Nrp1,mdiff)    vertical diffusivities              (m^2/s)
1104  c     kbl(imt)                    -1 of first grid level below hbl  c     kbl   (imt)                 -1 of first grid level below hbl
1105        _KPP_RL ustar (imt)  c     myThid               thread number for this instance of the routine
1106        _KPP_RL bfsfc (imt)        integer myThid
1107        _KPP_RL hbl   (imt)        _RL ustar (imt)
1108        _KPP_RL stable(imt)        _RL bfsfc (imt)
1109        _KPP_RL casea (imt)        _RL hbl   (imt)
1110        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL stable(imt)
1111          _RL casea (imt)
1112          _RL diffus(imt,0:Nrp1,mdiff)
1113        integer kbl(imt)        integer kbl(imt)
1114    
1115  c output  c output
# Line 1101  c     dkm1 (imt,mdiff)            bounda Line 1117  c     dkm1 (imt,mdiff)            bounda
1117  c     blmc (imt,Nr,mdiff)         boundary layer mixing coefficients  (m^2/s)  c     blmc (imt,Nr,mdiff)         boundary layer mixing coefficients  (m^2/s)
1118  c     ghat (imt,Nr)               nonlocal scalar transport  c     ghat (imt,Nr)               nonlocal scalar transport
1119  c     sigma(imt)                  normalized depth (d / hbl)  c     sigma(imt)                  normalized depth (d / hbl)
1120        _KPP_RL dkm1 (imt,mdiff)        _RL dkm1 (imt,mdiff)
1121        _KPP_RL blmc (imt,Nr,mdiff)        _RL blmc (imt,Nr,mdiff)
1122        _KPP_RL ghat (imt,Nr)        _RL ghat (imt,Nr)
1123        _KPP_RL sigma(imt)        _RL sigma(imt)
1124        integer ikppkey, kkppkey        integer ikppkey, kkppkey
1125    
1126  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
# Line 1113  c  local Line 1129  c  local
1129  c     gat1*(imt)                 shape function at sigma = 1  c     gat1*(imt)                 shape function at sigma = 1
1130  c     dat1*(imt)                 derivative of shape function at sigma = 1  c     dat1*(imt)                 derivative of shape function at sigma = 1
1131  c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)  c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)
1132        _KPP_RL gat1m(imt), gat1s(imt), gat1t(imt)        _RL gat1m(imt), gat1s(imt), gat1t(imt)
1133        _KPP_RL dat1m(imt), dat1s(imt), dat1t(imt)        _RL dat1m(imt), dat1s(imt), dat1t(imt)
1134        _KPP_RL ws(imt), wm(imt)        _RL ws(imt), wm(imt)
1135        integer i, kn, ki        integer i, kn, ki
1136        _KPP_RL R, dvdzup, dvdzdn, viscp        _RL R, dvdzup, dvdzdn, viscp
1137        _KPP_RL difsp, diftp, visch, difsh, difth        _RL difsp, diftp, visch, difsh, difth
1138        _KPP_RL f1, sig, a1, a2, a3, delhat        _RL f1, sig, a1, a2, a3, delhat
1139        _KPP_RL Gm, Gs, Gt        _RL Gm, Gs, Gt
1140        _KPP_RL tempVar        _RL tempVar
1141    
1142        _KPP_RL    p0    , eins        _RL    p0    , eins
1143        parameter (p0=0.0, eins=1.0)        parameter (p0=0.0, eins=1.0)
1144    
1145  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 1137  c--------------------------------------- Line 1153  c---------------------------------------
1153  CADJ STORE sigma = comlev1_kpp, key = ikppkey  CADJ STORE sigma = comlev1_kpp, key = ikppkey
1154        call wscale (        call wscale (
1155       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1156       O        wm, ws )       O        wm, ws, myThid )
1157  CADJ STORE wm = comlev1_kpp, key = ikppkey  CADJ STORE wm = comlev1_kpp, key = ikppkey
1158  CADJ STORE ws = comlev1_kpp, key = ikppkey  CADJ STORE ws = comlev1_kpp, key = ikppkey
1159    
# Line 1230  CADJ STORE ws = comlev1_kpp_k, key = kkp Line 1246  CADJ STORE ws = comlev1_kpp_k, key = kkp
1246  CADJ STORE sigma = comlev1_kpp_k, key = kkppkey  CADJ STORE sigma = comlev1_kpp_k, key = kkppkey
1247           call wscale (           call wscale (
1248       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1249       O        wm, ws )       O        wm, ws, myThid )
1250  CADJ STORE wm = comlev1_kpp_k, key = kkppkey  CADJ STORE wm = comlev1_kpp_k, key = kkppkey
1251  CADJ STORE ws = comlev1_kpp_k, key = kkppkey  CADJ STORE ws = comlev1_kpp_k, key = kkppkey
1252    
# Line 1283  CADJ STORE ws = comlev1_kpp, key = ikppk Line 1299  CADJ STORE ws = comlev1_kpp, key = ikppk
1299  CADJ STORE sigma = comlev1_kpp, key = ikppkey  CADJ STORE sigma = comlev1_kpp, key = ikppkey
1300        call wscale (        call wscale (
1301       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1302       O        wm, ws )       O        wm, ws, myThid )
1303  CADJ STORE wm = comlev1_kpp, key = ikppkey  CADJ STORE wm = comlev1_kpp, key = ikppkey
1304  CADJ STORE ws = comlev1_kpp, key = ikppkey  CADJ STORE ws = comlev1_kpp, key = ikppkey
1305    
# Line 1311  c*************************************** Line 1327  c***************************************
1327       I       dkm1, hbl, kbl, diffus, casea       I       dkm1, hbl, kbl, diffus, casea
1328       U     , ghat       U     , ghat
1329       O     , blmc       O     , blmc
1330       &     )       &     , myThid )
1331    
1332  c enhance the diffusivity at the kbl-.5 interface  c enhance the diffusivity at the kbl-.5 interface
1333    
# Line 1326  c     hbl(imt)                  boundary Line 1342  c     hbl(imt)                  boundary
1342  c     kbl(imt)                  grid above hbl  c     kbl(imt)                  grid above hbl
1343  c     diffus(imt,0:Nrp1,mdiff) vertical diffusivities           (m^2/s)  c     diffus(imt,0:Nrp1,mdiff) vertical diffusivities           (m^2/s)
1344  c     casea(imt)                = 1 in caseA, = 0 in case B  c     casea(imt)                = 1 in caseA, = 0 in case B
1345        _KPP_RL dkm1  (imt,mdiff)  c     myThid                    thread number for this instance of the routine
1346        _KPP_RL hbl   (imt)        integer   myThid
1347          _RL dkm1  (imt,mdiff)
1348          _RL hbl   (imt)
1349        integer kbl   (imt)        integer kbl   (imt)
1350        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL diffus(imt,0:Nrp1,mdiff)
1351        _KPP_RL casea (imt)        _RL casea (imt)
1352    
1353  c input/output  c input/output
1354  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)
1355        _KPP_RL ghat (imt,Nr)        _RL ghat (imt,Nr)
1356    
1357  c output  c output
1358  c     enhanced bound. layer mixing coeff.  c     enhanced bound. layer mixing coeff.
1359        _KPP_RL blmc  (imt,Nr,mdiff)        _RL blmc  (imt,Nr,mdiff)
1360    
1361  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1362    
1363  c local  c local
1364  c     fraction hbl lies beteen zgrid neighbors  c     fraction hbl lies beteen zgrid neighbors
1365        _KPP_RL delta        _RL delta
1366        integer ki, i, md        integer ki, i, md
1367        _KPP_RL dkmp5, dstar        _RL dkmp5, dstar
1368    
1369        do i = 1, imt        do i = 1, imt
1370           ki = kbl(i)-1           ki = kbl(i)-1
# Line 1372  c     fraction hbl lies beteen zgrid nei Line 1390  c     fraction hbl lies beteen zgrid nei
1390  c*************************************************************************  c*************************************************************************
1391    
1392        SUBROUTINE STATEKPP (        SUBROUTINE STATEKPP (
1393       I     ikppkey, bi, bj, myThid,       O     RHO1, DBLOC, DBSFC, TTALPHA, SSBETA,
1394       O     RHO1, DBLOC, DBSFC, TTALPHA, SSBETA)       I     ikppkey, bi, bj, myThid )
1395  c  c
1396  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1397  c     "statekpp" computes all necessary input arrays  c     "statekpp" computes all necessary input arrays
# Line 1398  c Line 1416  c
1416  c     written  by: jan morzel,   feb. 10, 1995 (converted from "sigma" version)  c     written  by: jan morzel,   feb. 10, 1995 (converted from "sigma" version)
1417  c     modified by: d. menemenlis,    june 1998 : for use with MIT GCM UV  c     modified by: d. menemenlis,    june 1998 : for use with MIT GCM UV
1418  c  c
1419    c     28 april 05: added computation of mixed layer depth KPPmld
1420    c     for a deltaRho equivalent to a deltaTheta of 0.8 deg C
1421    
1422  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1423    
1424        IMPLICIT NONE        IMPLICIT NONE
# Line 1407  c--------------------------------------- Line 1428  c---------------------------------------
1428  #include "PARAMS.h"  #include "PARAMS.h"
1429  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
1430  #include "DYNVARS.h"  #include "DYNVARS.h"
1431    #include "GRID.h"
1432    
1433  c-------------- Routine arguments -----------------------------------------  c-------------- Routine arguments -----------------------------------------
1434        INTEGER bi, bj, myThid        INTEGER bi, bj, myThid
1435        _KPP_RL RHO1   ( ibot:itop, jbot:jtop       )        _RL RHO1   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy       )
1436        _KPP_RL DBLOC  ( ibot:itop, jbot:jtop, Nr   )        _RL DBLOC  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
1437        _KPP_RL DBSFC  ( ibot:itop, jbot:jtop, Nr   )        _RL DBSFC  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
1438        _KPP_RL TTALPHA( ibot:itop, jbot:jtop, Nrp1 )        _RL TTALPHA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1439        _KPP_RL SSBETA ( ibot:itop, jbot:jtop, Nrp1 )        _RL SSBETA ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1440    
1441  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1442    
# Line 1425  c Line 1447  c
1447  c     rhok         - density of t(k  ) & s(k  ) at depth k  c     rhok         - density of t(k  ) & s(k  ) at depth k
1448  c     rhokm1       - density of t(k-1) & s(k-1) at depth k  c     rhokm1       - density of t(k-1) & s(k-1) at depth k
1449  c     rho1k        - density of t(1  ) & s(1  ) at depth k  c     rho1k        - density of t(1  ) & s(1  ) at depth k
1450  c     work1, work2 - work arrays for holding horizontal slabs  c     work1,2,3    - work arrays for holding horizontal slabs
1451    
1452        _RL RHOK  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL RHOK  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1453        _RL RHOKM1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL RHOKM1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 1433  c     work1, work2 - work arrays for hol Line 1455  c     work1, work2 - work arrays for hol
1455        _RL WORK1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1456        _RL WORK2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1457        _RL WORK3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL WORK3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1458    #ifdef ALLOW_DIAGNOSTICS
1459    c     KPPMLD - mixed layer depth based on density criterion
1460          _RL KPPMLD(1    :sNx    ,1    :sNy    )
1461    #endif /* ALLOW_DIAGNOSTICS */
1462    
1463        INTEGER I, J, K        INTEGER I, J, K
1464        INTEGER ikppkey, kkppkey        INTEGER ikppkey, kkppkey
1465    
# Line 1445  CADJ STORE theta(:,:,1,bi,bj) = comlev1_ Line 1472  CADJ STORE theta(:,:,1,bi,bj) = comlev1_
1472  CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey  CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1473  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1474        call FIND_RHO(        call FIND_RHO(
1475       I     bi, bj, ibot, itop, jbot, jtop, 1, 1,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1476       I     theta, salt,       I     theta, salt,
1477       O     WORK1,       O     WORK1,
1478       I     myThid )       I     myThid )
# Line 1455  CADJ STORE salt (:,:,1,bi,bj) = comlev1_ Line 1482  CADJ STORE salt (:,:,1,bi,bj) = comlev1_
1482  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1483    
1484        call FIND_ALPHA(        call FIND_ALPHA(
1485       I     bi, bj, ibot, itop, jbot, jtop, 1, 1,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1486       O     WORK2 )       O     WORK2, myThid )
1487    
1488        call FIND_BETA(        call FIND_BETA(
1489       I     bi, bj, ibot, itop, jbot, jtop, 1, 1,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1490       O     WORK3 )       O     WORK3, myThid )
1491    
1492        DO J = jbot, jtop        DO J = 1-OLy, sNy+OLy
1493           DO I = ibot, itop           DO I = 1-OLx, sNx+OLx
1494              RHO1(I,J)      = WORK1(I,J) + rhoConst              RHO1(I,J)      = WORK1(I,J) + rhoConst
1495              TTALPHA(I,J,1) = WORK2(I,J)              TTALPHA(I,J,1) = WORK2(I,J)
1496              SSBETA(I,J,1)  = WORK3(I,J)              SSBETA(I,J,1)  = WORK3(I,J)
# Line 1471  CADJ STORE salt (:,:,1,bi,bj) = comlev1_ Line 1498  CADJ STORE salt (:,:,1,bi,bj) = comlev1_
1498           END DO           END DO
1499        END DO        END DO
1500    
1501    #ifdef ALLOW_DIAGNOSTICS
1502    c     work3 - density of t(1)-.8 & s(1  ) at depth 1
1503          IF ( useDiagnostics ) THEN
1504             DO J = 1, sNy
1505                DO I = 1, sNx
1506                   KPPMLD(I,J) = ABS(R_low(I,J,bi,bj))
1507                   WORK3 (I,J) = WORK1(I,J) - 0.8 _d 0 * WORK2(I,J)
1508                END DO
1509             END DO
1510          ENDIF
1511    #endif /* ALLOW_DIAGNOSTICS */
1512    
1513  c calculate alpha, beta, and gradients in interior layers  c calculate alpha, beta, and gradients in interior layers
1514    
1515  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)  CHPF$  INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)
# Line 1483  CADJ STORE theta(:,:,k,bi,bj) = comlev1_ Line 1522  CADJ STORE theta(:,:,k,bi,bj) = comlev1_
1522  CADJ STORE salt (:,:,k,bi,bj) = comlev1_kpp_k, key=kkppkey  CADJ STORE salt (:,:,k,bi,bj) = comlev1_kpp_k, key=kkppkey
1523  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1524           call FIND_RHO(           call FIND_RHO(
1525       I        bi, bj, ibot, itop, jbot, jtop, K, K,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1526       I        theta, salt,       I        theta, salt,
1527       O        RHOK,       O        RHOK,
1528       I        myThid )       I        myThid )
# Line 1493  CADJ STORE theta(:,:,k-1,bi,bj) = comlev Line 1532  CADJ STORE theta(:,:,k-1,bi,bj) = comlev
1532  CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_kpp_k, key=kkppkey  CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_kpp_k, key=kkppkey
1533  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1534           call FIND_RHO(           call FIND_RHO(
1535       I        bi, bj, ibot, itop, jbot, jtop, K-1, K,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K-1, K,
1536       I        theta, salt,       I        theta, salt,
1537       O        RHOKM1,       O        RHOKM1,
1538       I        myThid )       I        myThid )
# Line 1503  CADJ STORE theta(:,:,1,bi,bj) = comlev1_ Line 1542  CADJ STORE theta(:,:,1,bi,bj) = comlev1_
1542  CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey  CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1543  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1544           call FIND_RHO(           call FIND_RHO(
1545       I        bi, bj, ibot, itop, jbot, jtop, 1, K,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, K,
1546       I        theta, salt,       I        theta, salt,
1547       O        RHO1K,       O        RHO1K,
1548       I        myThid )       I        myThid )
# Line 1515  CADJ STORE rho1k (:,:)          = comlev Line 1554  CADJ STORE rho1k (:,:)          = comlev
1554  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1555    
1556           call FIND_ALPHA(           call FIND_ALPHA(
1557       I        bi, bj, ibot, itop, jbot, jtop, K, K,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1558       O        WORK1 )       O        WORK1, myThid )
1559    
1560           call FIND_BETA(           call FIND_BETA(
1561       I        bi, bj, ibot, itop, jbot, jtop, K, K,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1562       O        WORK2 )       O        WORK2, myThid )
1563    
1564           DO J = jbot, jtop           DO J = 1-OLy, sNy+OLy
1565              DO I = ibot, itop              DO I = 1-OLx, sNx+OLx
1566                 TTALPHA(I,J,K) = WORK1 (I,J)                 TTALPHA(I,J,K) = WORK1 (I,J)
1567                 SSBETA(I,J,K)  = WORK2 (I,J)                 SSBETA(I,J,K)  = WORK2 (I,J)
1568                 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 1533  CADJ STORE rho1k (:,:)          = comlev Line 1572  CADJ STORE rho1k (:,:)          = comlev
1572              END DO              END DO
1573           END DO           END DO
1574    
1575    #ifdef ALLOW_DIAGNOSTICS
1576             IF ( useDiagnostics ) THEN
1577    c     work1 - density of t(k-1)  & s(k-1) at depth 1
1578    c     work2 - density of t(k  )  & s(k  ) at depth 1
1579    c     work3 - density of t(1)-.8 & s(1  ) at depth 1
1580                call FIND_RHO(
1581         I           bi, bj, 1, sNx, 1, sNy, K-1, 1, theta, salt,
1582         O           WORK1,
1583         I           myThid )
1584                call FIND_RHO(
1585         I           bi, bj, 1, sNx, 1, sNy, K  , 1, theta, salt,
1586         O           WORK2,
1587         I           myThid )
1588                DO J = 1, sNy
1589                   DO I = 1, sNx
1590                      IF ( k .LE. klowC(I,J,bi,bj) .AND.
1591         &                 WORK1(I,J) .LT. WORK3(I,J) .AND.
1592         &                 WORK2(I,J) .GE. WORK3(I,J) ) THEN
1593                         KPPMLD(I,J) = MIN ( KPPMLD(I,J),
1594         &                    ABS(((WORK3(I,J)-WORK1(I,J))*rC(k)+
1595         &                    (WORK2(I,J)-WORK3(I,J))*rC(k-1))/
1596         &                    (WORK2(I,J)-WORK1(I,J))))
1597                      ENDIF
1598                   END DO
1599                END DO
1600             ENDIF
1601    #endif /* ALLOW_DIAGNOSTICS */
1602    
1603        END DO        END DO
1604    
1605  c     compute arrays for K = Nrp1  c     compute arrays for K = Nrp1
1606        DO J = jbot, jtop        DO J = 1-OLy, sNy+OLy
1607           DO I = ibot, itop           DO I = 1-OLx, sNx+OLx
1608              TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)              TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)
1609              SSBETA(I,J,Nrp1)  = SSBETA(I,J,Nr)              SSBETA(I,J,Nrp1)  = SSBETA(I,J,Nr)
1610              DBLOC(I,J,Nr)     = 0.              DBLOC(I,J,Nr)     = 0.
1611           END DO           END DO
1612        END DO        END DO
1613    
1614    #ifdef ALLOW_DIAGNOSTICS
1615          IF ( useDiagnostics ) THEN
1616             CALL DIAGNOSTICS_FILL(KPPmld,'KPPmld  ',0,1,3,bi,bj,myThid)
1617          ENDIF
1618    #endif /* ALLOW_DIAGNOSTICS */
1619    
1620  #endif /* ALLOW_KPP */  #endif /* ALLOW_KPP */
1621    
1622        RETURN        RETURN
Legend:
Removed from v.1.22  
changed lines
  Added in v.1.35
  ViewVC Help
Powered by ViewVC 1.1.22