/[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.23 by dimitri, Fri Apr 29 18:47:02 2005 UTC revision 1.31 by jmc, Thu May 3 21:36:26 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
      I     , kmtj, shsq, dvsq, ustar  
23       I     , bo, bosol, dbloc, Ritop, coriol       I     , bo, bosol, dbloc, Ritop, coriol
24         I     , diffusKzS, diffusKzT
25       I     , ikppkey       I     , ikppkey
26       O     , diffus       O     , diffus
27       U     , ghat       U     , ghat
28       O     , hbl )       O     , hbl
29         I     , myTime, myIter, myThid )
30    
31  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
32  c  c
# Line 47  c--------------------------------------- Line 47  c---------------------------------------
47  #include "SIZE.h"  #include "SIZE.h"
48  #include "EEPARAMS.h"  #include "EEPARAMS.h"
49  #include "PARAMS.h"  #include "PARAMS.h"
 #include "DYNVARS.h"  
 #include "FFIELDS.h"  
50  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
51    
52  c input  c input
53  c     myTime          - current time in simulation  c   myTime :: Current time in simulation
54  c     myThid          - thread number for this instance of the routine  c   myIter :: Current iteration number in simulation
55  c     kmtj   (imt)    - number of vertical layers on this row  c   myThid :: My Thread Id. number
56  c     shsq   (imt,Nr) - (local velocity shear)^2                      ((m/s)^2)  c     kmtj   (imt)     - number of vertical layers on this row
57  c     dvsq   (imt,Nr) - (velocity shear re sfc)^2                     ((m/s)^2)  c     msk    (imt)     - surface mask (=1 if water, =0 otherwise)
58  c     ustar  (imt)    - surface friction velocity                         (m/s)  c     shsq   (imt,Nr)  - (local velocity shear)^2                     ((m/s)^2)
59  c     bo     (imt)    - surface turbulent buoy. forcing               (m^2/s^3)  c     dvsq   (imt,Nr)  - (velocity shear re sfc)^2                    ((m/s)^2)
60  c     bosol  (imt)    - radiative buoyancy forcing                    (m^2/s^3)  c     ustar  (imt)     - surface friction velocity                        (m/s)
61  c     dbloc  (imt,Nr) - local delta buoyancy across interfaces          (m/s^2)  c     bo     (imt)     - surface turbulent buoy. forcing              (m^2/s^3)
62  c     dblocSm(imt,Nr) - horizontally smoothed dbloc                     (m/s^2)  c     bosol  (imt)     - radiative buoyancy forcing                   (m^2/s^3)
63  c                         stored in ghat to save space  c     dbloc  (imt,Nr)  - local delta buoyancy across interfaces         (m/s^2)
64  c     Ritop  (imt,Nr) - numerator of bulk Richardson Number  c     dblocSm(imt,Nr)  - horizontally smoothed dbloc                    (m/s^2)
65  c                         (zref-z) * delta buoyancy w.r.t. surface    ((m/s)^2)  c                          stored in ghat to save space
66  c     coriol (imt)    - Coriolis parameter                                (1/s)  c     Ritop  (imt,Nr)  - numerator of bulk Richardson Number
67    c                          (zref-z) * delta buoyancy w.r.t. surface   ((m/s)^2)
68    c     coriol (imt)     - Coriolis parameter                               (1/s)
69    c     diffusKzS(imt,Nr)- background vertical diffusivity for scalars    (m^2/s)
70    c     diffusKzT(imt,Nr)- background vertical diffusivity for theta      (m^2/s)
71  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,
72  c           e.g., hbl(sNx,sNy) -> hbl(imt),  c           e.g., hbl(sNx,sNy) -> hbl(imt),
73  c           where hbl(i,j) -> hbl((j-1)*sNx+i)  c           where hbl(i,j) -> hbl((j-1)*sNx+i)
74    
75        _RL     mytime        _RL     myTime
76        integer mythid        integer myIter
77        integer kmtj  (imt   )        integer myThid
78        _KPP_RL shsq  (imt,Nr)        integer kmtj (imt   )
79        _KPP_RL dvsq  (imt,Nr)        _RL shsq     (imt,Nr)
80        _KPP_RL ustar (imt   )        _RL dvsq     (imt,Nr)
81        _KPP_RL bo    (imt   )        _RL ustar    (imt   )
82        _KPP_RL bosol (imt   )        _RL bo       (imt   )
83        _KPP_RL dbloc (imt,Nr)        _RL bosol    (imt   )
84        _KPP_RL Ritop (imt,Nr)        _RL dbloc    (imt,Nr)
85        _KPP_RL coriol(imt   )        _RL Ritop    (imt,Nr)
86          _RL coriol   (imt   )
87          _RS msk      (imt   )
88          _RL diffusKzS(imt,Nr)
89          _RL diffusKzT(imt,Nr)
90    
91        integer ikppkey        integer ikppkey
92    
# Line 91  c     diffus (imt,3)  - vertical tempera Line 97  c     diffus (imt,3)  - vertical tempera
97  c     ghat   (imt)    - nonlocal transport coefficient                  (s/m^2)  c     ghat   (imt)    - nonlocal transport coefficient                  (s/m^2)
98  c     hbl    (imt)    - mixing layer depth                                  (m)  c     hbl    (imt)    - mixing layer depth                                  (m)
99    
100        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL diffus(imt,0:Nrp1,mdiff)
101        _KPP_RL ghat  (imt,Nr)        _RL ghat  (imt,Nr)
102        _KPP_RL hbl   (imt)        _RL hbl   (imt)
103    
104  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
105    
# Line 107  c     blmc   (imt,Nr,mdiff) - boundary l Line 113  c     blmc   (imt,Nr,mdiff) - boundary l
113  c     sigma  (imt         ) - normalized depth (d / hbl)  c     sigma  (imt         ) - normalized depth (d / hbl)
114  c     Rib    (imt,Nr      ) - bulk Richardson number  c     Rib    (imt,Nr      ) - bulk Richardson number
115    
116        integer kbl   (imt         )        integer kbl(imt         )
117        _KPP_RL bfsfc (imt         )        _RL bfsfc  (imt         )
118        _KPP_RL casea (imt         )        _RL casea  (imt         )
119        _KPP_RL stable(imt         )        _RL stable (imt         )
120        _KPP_RL dkm1  (imt,   mdiff)        _RL dkm1   (imt,   mdiff)
121        _KPP_RL blmc  (imt,Nr,mdiff)        _RL blmc   (imt,Nr,mdiff)
122        _KPP_RL sigma (imt         )        _RL sigma  (imt         )
123        _KPP_RL Rib   (imt,Nr      )        _RL Rib    (imt,Nr      )
124    
125        integer i, k, md        integer i, k, md
126    
# Line 132  CADJ STORE dbloc = comlev1_kpp, key = ik Line 138  CADJ STORE dbloc = comlev1_kpp, key = ik
138  cph)  cph)
139        call Ri_iwmix (        call Ri_iwmix (
140       I       kmtj, shsq, dbloc, ghat       I       kmtj, shsq, dbloc, ghat
141         I     , diffusKzS, diffusKzT
142       I     , ikppkey       I     , ikppkey
143       O     , diffus )       O     , diffus )
144    
# Line 164  c diagnose the new boundary layer depth Line 171  c diagnose the new boundary layer depth
171  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
172    
173        call bldepth (        call bldepth (
174       I       mytime, mythid       I       kmtj
      I     , kmtj  
175       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol
176       I     , ikppkey       I     , ikppkey
177       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma
178       &     )       I     , myTime, myIter, myThid )
179    
180  CADJ STORE hbl,bfsfc,stable,casea,kbl = comlev1_kpp, key = ikppkey  CADJ STORE hbl,bfsfc,stable,casea,kbl = comlev1_kpp, key = ikppkey
181    
# Line 180  c--------------------------------------- Line 186  c---------------------------------------
186        call blmix (        call blmix (
187       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl
188       O     , dkm1, blmc, ghat, sigma, ikppkey       O     , dkm1, blmc, ghat, sigma, ikppkey
189       &     )       I     , myThid )
190  cph(  cph(
191  CADJ STORE dkm1,blmc,ghat = comlev1_kpp, key = ikppkey  CADJ STORE dkm1,blmc,ghat = comlev1_kpp, key = ikppkey
192  CADJ STORE hbl, kbl, diffus, casea = comlev1_kpp, key = ikppkey  CADJ STORE hbl, kbl, diffus, casea = comlev1_kpp, key = ikppkey
# Line 193  c--------------------------------------- Line 199  c---------------------------------------
199        call enhance (        call enhance (
200       I       dkm1, hbl, kbl, diffus, casea       I       dkm1, hbl, kbl, diffus, casea
201       U     , ghat       U     , ghat
202       O     , blmc )       O     , blmc
203         I     , myThid )
204    
205  cph(  cph(
206  cph avoids recomp. of enhance  cph avoids recomp. of enhance
# Line 209  c--------------------------------------- Line 216  c---------------------------------------
216        do k = 1, Nr        do k = 1, Nr
217           do i = 1, imt           do i = 1, imt
218              if (k .lt. kbl(i)) then              if (k .lt. kbl(i)) then
219    #ifdef ALLOW_SHELFICE
220    C     when there is shelfice on top (msk(i)=0), reset the boundary layer
221    C     mixing coefficients blmc to pure Ri-number based mixing
222                   blmc(i,k,1) = max ( blmc(i,k,1)*msk(i),
223         &              diffus(i,k,1) )
224                   blmc(i,k,2) = max ( blmc(i,k,2)*msk(i),
225         &              diffus(i,k,2) )
226                   blmc(i,k,3) = max ( blmc(i,k,3)*msk(i),
227         &              diffus(i,k,3) )
228    #endif /* not ALLOW_SHELFICE */
229                 diffus(i,k,1) = max ( blmc(i,k,1), viscAr  )                 diffus(i,k,1) = max ( blmc(i,k,1), viscAr  )
230                 diffus(i,k,2) = max ( blmc(i,k,2), diffKrNrS(k) )                 diffus(i,k,2) = max ( blmc(i,k,2), diffusKzS(i,Nr) )
231                 diffus(i,k,3) = max ( blmc(i,k,3), diffKrNrT(k) )                 diffus(i,k,3) = max ( blmc(i,k,3), diffusKzT(i,Nr) )
232              else              else
233                 ghat(i,k) = 0.                 ghat(i,k) = 0.
234              endif              endif
# Line 226  c--------------------------------------- Line 243  c---------------------------------------
243  c*************************************************************************  c*************************************************************************
244    
245        subroutine bldepth (        subroutine bldepth (
246       I       mytime, mythid       I       kmtj
      I     , kmtj  
247       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol       I     , dvsq, dbloc, Ritop, ustar, bo, bosol, coriol
248       I     , ikppkey       I     , ikppkey
249       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma       O     , hbl, bfsfc, stable, casea, kbl, Rib, sigma
250       &     )       I     , myTime, myIter, myThid )
251    
252  c     the oceanic planetary boundary layer depth, hbl, is determined as  c     the oceanic planetary boundary layer depth, hbl, is determined as
253  c     the shallowest depth where the bulk Richardson number is  c     the shallowest depth where the bulk Richardson number is
# Line 257  c Line 273  c
273        IMPLICIT NONE        IMPLICIT NONE
274    
275  #include "SIZE.h"  #include "SIZE.h"
 #include "EEPARAMS.h"  
 #include "PARAMS.h"  
276  #include "KPP_PARAMS.h"  #include "KPP_PARAMS.h"
 #include "FFIELDS.h"  
277    
278  c input  c input
279  c------  c------
280  c myTime    : current time in simulation  c   myTime :: Current time in simulation
281  c myThid    : thread number for this instance of the routine  c   myIter :: Current iteration number in simulation
282    c   myThid :: My Thread Id. number
283  c kmtj      : number of vertical layers  c kmtj      : number of vertical layers
284  c dvsq      : (velocity shear re sfc)^2             ((m/s)^2)  c dvsq      : (velocity shear re sfc)^2             ((m/s)^2)
285  c dbloc     : local delta buoyancy across interfaces  (m/s^2)  c dbloc     : local delta buoyancy across interfaces  (m/s^2)
# Line 276  c ustar     : surface friction velocity Line 290  c ustar     : surface friction velocity
290  c bo        : surface turbulent buoyancy forcing    (m^2/s^3)  c bo        : surface turbulent buoyancy forcing    (m^2/s^3)
291  c bosol     : radiative buoyancy forcing            (m^2/s^3)  c bosol     : radiative buoyancy forcing            (m^2/s^3)
292  c coriol    : Coriolis parameter                        (1/s)  c coriol    : Coriolis parameter                        (1/s)
293        _RL     mytime        _RL     myTime
294        integer mythid        integer myIter
295          integer myThid
296        integer kmtj(imt)        integer kmtj(imt)
297        _KPP_RL dvsq  (imt,Nr)        _RL dvsq    (imt,Nr)
298        _KPP_RL dbloc (imt,Nr)        _RL dbloc   (imt,Nr)
299        _KPP_RL Ritop (imt,Nr)        _RL Ritop   (imt,Nr)
300        _KPP_RL ustar (imt)        _RL ustar   (imt)
301        _KPP_RL bo    (imt)        _RL bo      (imt)
302        _KPP_RL bosol (imt)        _RL bosol   (imt)
303        _KPP_RL coriol(imt)        _RL coriol  (imt)
304        integer ikppkey, kkppkey        integer ikppkey, kkppkey
305    
306  c  output  c  output
# Line 297  c casea     : =1 in case A, =0 in case B Line 312  c casea     : =1 in case A, =0 in case B
312  c kbl       : -1 of first grid level below hbl  c kbl       : -1 of first grid level below hbl
313  c Rib       : Bulk Richardson number  c Rib       : Bulk Richardson number
314  c sigma     : normalized depth (d/hbl)  c sigma     : normalized depth (d/hbl)
315        _KPP_RL hbl   (imt)        _RL hbl    (imt)
316        _KPP_RL bfsfc (imt)        _RL bfsfc  (imt)
317        _KPP_RL stable(imt)        _RL stable (imt)
318        _KPP_RL casea (imt)        _RL casea  (imt)
319        integer kbl   (imt)        integer kbl(imt)
320        _KPP_RL Rib   (imt,Nr)        _RL Rib    (imt,Nr)
321        _KPP_RL sigma (imt)        _RL sigma  (imt)
322    
323  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
324    
325  c  local  c  local
326  c-------  c-------
327  c wm, ws    : turbulent velocity scales         (m/s)  c wm, ws    : turbulent velocity scales         (m/s)
328        _KPP_RL wm(imt), ws(imt)        _RL wm(imt), ws(imt)
329        _RL     worka(imt)        _RL worka(imt)
330    
331        _KPP_RL bvsq, vtsq, hekman, hmonob, hlimit, tempVar1, tempVar2        _RL bvsq, vtsq, hekman, hmonob, hlimit, tempVar1, tempVar2
332        integer i, kl        integer i, kl
333    
334        _KPP_RL     p5    , eins        _RL         p5    , eins
335        parameter ( p5=0.5, eins=1.0 )        parameter ( p5=0.5, eins=1.0 )
336        _RL         minusone        _RL         minusone
337        parameter ( minusone=-1.0 )        parameter ( minusone=-1.0 )
# Line 351  c     compute bfsfc = sw fraction at hbf Line 366  c     compute bfsfc = sw fraction at hbf
366           end do           end do
367  CADJ store worka = comlev1_kpp_k, key = kkppkey  CADJ store worka = comlev1_kpp_k, key = kkppkey
368           call SWFRAC(           call SWFRAC(
369       I        imt, hbf,       I       imt, hbf,
370       I        mytime, mythid,       U       worka,
371       U        worka )       I       myTime, myIter, myThid )
372  CADJ store worka = comlev1_kpp_k, key = kkppkey  CADJ store worka = comlev1_kpp_k, key = kkppkey
373    
374    
# Line 377  c--------------------------------------- Line 392  c---------------------------------------
392    
393           call wscale (           call wscale (
394       I        sigma, casea, ustar, bfsfc,       I        sigma, casea, ustar, bfsfc,
395       O        wm, ws )       O        wm, ws, myThid )
396  CADJ store ws = comlev1_kpp_k, key = kkppkey  CADJ store ws = comlev1_kpp_k, key = kkppkey
397    
398           do i = 1, imt           do i = 1, imt
# Line 454  c--------------------------------------- Line 469  c---------------------------------------
469  CADJ store worka = comlev1_kpp  CADJ store worka = comlev1_kpp
470  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
471        call SWFRAC(        call SWFRAC(
472       I     imt, minusone,       I       imt, minusone,
473       I     mytime, mythid,       U       worka,
474       U     worka )       I       myTime, myIter, myThid )
475  CADJ store worka = comlev1_kpp  CADJ store worka = comlev1_kpp
476  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
477    
# Line 527  CADJ store worka = comlev1_kpp Line 542  CADJ store worka = comlev1_kpp
542  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
543        call SWFRAC(        call SWFRAC(
544       I     imt, minusone,       I     imt, minusone,
545       I     mytime, mythid,       U     worka,
546       U     worka )       I     myTime, myIter, myThid )
547  CADJ store worka = comlev1_kpp  CADJ store worka = comlev1_kpp
548  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)  CADJ &   , key = ikppkey, shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
549    
# Line 562  c*************************************** Line 577  c***************************************
577    
578        subroutine wscale (        subroutine wscale (
579       I     sigma, hbl, ustar, bfsfc,       I     sigma, hbl, ustar, bfsfc,
580       O     wm, ws )       O     wm, ws,
581         I     myThid )
582    
583  c     compute turbulent velocity scales.  c     compute turbulent velocity scales.
584  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 599  c sigma   : normalized depth (d/hbl)
599  c hbl     : boundary layer depth (m)  c hbl     : boundary layer depth (m)
600  c ustar   : surface friction velocity         (m/s)  c ustar   : surface friction velocity         (m/s)
601  c bfsfc   : total surface buoyancy flux       (m^2/s^3)  c bfsfc   : total surface buoyancy flux       (m^2/s^3)
602        _KPP_RL sigma(imt)  c myThid  : thread number for this instance of the routine
603        _KPP_RL hbl  (imt)        integer myThid
604        _KPP_RL ustar(imt)        _RL sigma(imt)
605        _KPP_RL bfsfc(imt)        _RL hbl  (imt)
606          _RL ustar(imt)
607          _RL bfsfc(imt)
608    
609  c  output  c  output
610  c--------  c--------
611  c wm, ws  : turbulent velocity scales at sigma  c wm, ws  : turbulent velocity scales at sigma
612        _KPP_RL wm(imt), ws(imt)        _RL wm(imt), ws(imt)
613    
614  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
615    
616  c local  c local
617  c------  c------
618  c zehat   : = zeta *  ustar**3  c zehat   : = zeta *  ustar**3
619        _KPP_RL zehat        _RL zehat
620    
621        integer iz, izp1, ju, i, jup1        integer iz, izp1, ju, i, jup1
622        _KPP_RL udiff, zdiff, zfrac, ufrac, fzfrac, wam        _RL udiff, zdiff, zfrac, ufrac, fzfrac, wam
623        _KPP_RL wbm, was, wbs, u3, tempVar        _RL wbm, was, wbs, u3, tempVar
624    
625  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
626  c use lookup table for zehat < zmax only; otherwise use  c use lookup table for zehat < zmax only; otherwise use
# Line 657  c--------------------------------------- Line 675  c---------------------------------------
675  c*************************************************************************  c*************************************************************************
676    
677        subroutine Ri_iwmix (        subroutine Ri_iwmix (
678       I       kmtj, shsq, dbloc, dblocSm       I       kmtj, shsq, dbloc, dblocSm,
679       I     , ikppkey       I       diffusKzS, diffusKzT,
680       O     , diffus )       I       ikppkey,
681         O       diffus,
682         I       myThid )
683    
684  c     compute interior viscosity diffusivity coefficients due  c     compute interior viscosity diffusivity coefficients due
685  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 698  c     kmtj   (imt)         number of ver
698  c     shsq   (imt,Nr)      (local velocity shear)^2               ((m/s)^2)  c     shsq   (imt,Nr)      (local velocity shear)^2               ((m/s)^2)
699  c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)  c     dbloc  (imt,Nr)      local delta buoyancy                     (m/s^2)
700  c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)  c     dblocSm(imt,Nr)      horizontally smoothed dbloc              (m/s^2)
701        integer kmtj   (imt)  c     diffusKzS(imt,Nr)- background vertical diffusivity for scalars    (m^2/s)
702        _KPP_RL shsq   (imt,Nr)  c     diffusKzT(imt,Nr)- background vertical diffusivity for theta      (m^2/s)
703        _KPP_RL dbloc  (imt,Nr)  c     myThid :: My Thread Id. number
704        _KPP_RL dblocSm(imt,Nr)        integer kmtj (imt)
705          _RL shsq     (imt,Nr)
706          _RL dbloc    (imt,Nr)
707          _RL dblocSm  (imt,Nr)
708          _RL diffusKzS(imt,Nr)
709          _RL diffusKzT(imt,Nr)
710        integer ikppkey        integer ikppkey
711          integer myThid
712    
713  c output  c output
714  c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)  c     diffus(imt,0:Nrp1,1)  vertical viscosivity coefficient        (m^2/s)
715  c     diffus(imt,0:Nrp1,2)  vertical scalar diffusivity             (m^2/s)  c     diffus(imt,0:Nrp1,2)  vertical scalar diffusivity             (m^2/s)
716  c     diffus(imt,0:Nrp1,3)  vertical temperature diffusivity        (m^2/s)  c     diffus(imt,0:Nrp1,3)  vertical temperature diffusivity        (m^2/s)
717        _KPP_RL diffus(imt,0:Nrp1,3)        _RL diffus(imt,0:Nrp1,3)
718    
719  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
720    
721  c local variables  c local variables
722  c     Rig                   local Richardson number  c     Rig                   local Richardson number
723  c     fRi, fcon             function of Rig  c     fRi, fcon             function of Rig
724        _KPP_RL Rig        _RL Rig
725        _KPP_RL fRi, fcon        _RL fRi, fcon
726        _KPP_RL ratio        _RL ratio
727        integer i, ki        integer i, ki
728        _KPP_RL c1, c0        _RL c1, c0
729    
730  #ifdef ALLOW_KPP_VERTICALLY_SMOOTH  #ifdef ALLOW_KPP_VERTICALLY_SMOOTH
731        integer mr        integer mr
# Line 755  CADJ store diffus(:,:,1) = kpp_ri_tape_m Line 781  CADJ store diffus(:,:,1) = kpp_ri_tape_m
781  CADJ &  , key=mr, shape=(/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2 /)  CADJ &  , key=mr, shape=(/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2 /)
782    
783           call z121 (           call z121 (
784       U     diffus(1,0,1))       U     diffus(1,0,1)
785         I     myThid )
786        end do        end do
787  #endif  #endif
788    
# Line 786  c    mixing due to internal waves, and s Line 813  c    mixing due to internal waves, and s
813  #ifndef EXCLUDE_KPP_SHEAR_MIX  #ifndef EXCLUDE_KPP_SHEAR_MIX
814              if ( .NOT. inAdMode ) then              if ( .NOT. inAdMode ) then
815                 diffus(i,ki,1) = viscAr + fcon * difmcon + fRi * difm0                 diffus(i,ki,1) = viscAr + fcon * difmcon + fRi * difm0
816                 diffus(i,ki,2) = diffKrNrS(ki)+ fcon*difscon + fRi*difs0                 diffus(i,ki,2) = diffusKzS(i,ki)+fcon*difscon+fRi*difs0
817                 diffus(i,ki,3) = diffKrNrT(ki)+ fcon*difscon + fRi*difs0                 diffus(i,ki,3) = diffusKzT(i,ki)+fcon*diftcon+fRi*dift0
818              else              else
819                 diffus(i,ki,1) = viscAr                 diffus(i,ki,1) = viscAr
820                 diffus(i,ki,2) = diffKrNrS(ki)                 diffus(i,ki,2) = diffusKzS(i,ki)
821                 diffus(i,ki,3) = diffKrNrT(ki)                 diffus(i,ki,3) = diffusKzT(i,ki)
822              endif              endif
823  #else  #else
824              diffus(i,ki,1) = viscAr              diffus(i,ki,1) = viscAr
825              diffus(i,ki,2) = diffKrNrS(ki)              diffus(i,ki,2) = diffusKzS(i,ki)
826              diffus(i,ki,3) = diffKrNrT(ki)              diffus(i,ki,3) = diffusKzT(i,ki)
827  #endif  #endif
828    
829           end do           end do
# Line 819  c         set surface values to 0.0 Line 846  c         set surface values to 0.0
846  c*************************************************************************  c*************************************************************************
847    
848        subroutine z121 (        subroutine z121 (
849       U     v )       U     v,
850         I     myThid )
851    
852  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)
853  c     top (0) value is used as a dummy  c     top (0) value is used as a dummy
# Line 837  c     penetrative convection. Line 865  c     penetrative convection.
865  c input/output  c input/output
866  c-------------  c-------------
867  c v     : 2-D array to be smoothed in Nrp1 direction  c v     : 2-D array to be smoothed in Nrp1 direction
868        _KPP_RL v(imt,0:Nrp1)  c myThid: thread number for this instance of the routine
869          integer myThid
870          _RL v(imt,0:Nrp1)
871    
872  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
873    
874  c local  c local
875        _KPP_RL zwork, zflag        _RL zwork, zflag
876        _KPP_RL KRi_range(1:Nrp1)        _RL KRi_range(1:Nrp1)
877        integer i, k, km1, kp1        integer i, k, km1, kp1
878    
879        _KPP_RL     p0      , p25       , p5      , p2        _RL         p0      , p25       , p5      , p2
880        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 )
881    
882        KRi_range(Nrp1) = p0        KRi_range(Nrp1) = p0
# Line 900  CADJ STORE v(i,k), zwork = z121tape Line 930  CADJ STORE v(i,k), zwork = z121tape
930    
931  c*************************************************************************  c*************************************************************************
932    
       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*************************************************************************  
   
933        subroutine smooth_horiz (        subroutine smooth_horiz (
934       I     k, bi, bj,       I     k, bi, bj,
935       U     fld )       U     fld,
936         I     myThid )
937    
938  c     Apply horizontal smoothing to global _RL 2-D array  c     Apply horizontal smoothing to global _RL 2-D array
939    
# Line 995  c     Apply horizontal smoothing to glob Line 945  c     Apply horizontal smoothing to glob
945  c     input  c     input
946  c     bi, bj : array indices  c     bi, bj : array indices
947  c     k      : vertical index used for masking  c     k      : vertical index used for masking
948    c     myThid : thread number for this instance of the routine
949          INTEGER myThid
950        integer k, bi, bj        integer k, bi, bj
951    
952  c     input/output  c     input/output
# Line 1065  c*************************************** Line 1017  c***************************************
1017        subroutine blmix (        subroutine blmix (
1018       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl       I       ustar, bfsfc, hbl, stable, casea, diffus, kbl
1019       O     , dkm1, blmc, ghat, sigma, ikppkey       O     , dkm1, blmc, ghat, sigma, ikppkey
1020       &     )       I     , myThid )
1021    
1022  c     mixing coefficients within boundary layer depend on surface  c     mixing coefficients within boundary layer depend on surface
1023  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 1039  c     hbl   (imt)                 bounda
1039  c     stable(imt)                 = 1 in stable forcing  c     stable(imt)                 = 1 in stable forcing
1040  c     casea (imt)                 = 1 in case A  c     casea (imt)                 = 1 in case A
1041  c     diffus(imt,0:Nrp1,mdiff)    vertical diffusivities              (m^2/s)  c     diffus(imt,0:Nrp1,mdiff)    vertical diffusivities              (m^2/s)
1042  c     kbl(imt)                    -1 of first grid level below hbl  c     kbl   (imt)                 -1 of first grid level below hbl
1043        _KPP_RL ustar (imt)  c     myThid               thread number for this instance of the routine
1044        _KPP_RL bfsfc (imt)        integer myThid
1045        _KPP_RL hbl   (imt)        _RL ustar (imt)
1046        _KPP_RL stable(imt)        _RL bfsfc (imt)
1047        _KPP_RL casea (imt)        _RL hbl   (imt)
1048        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL stable(imt)
1049          _RL casea (imt)
1050          _RL diffus(imt,0:Nrp1,mdiff)
1051        integer kbl(imt)        integer kbl(imt)
1052    
1053  c output  c output
# Line 1101  c     dkm1 (imt,mdiff)            bounda Line 1055  c     dkm1 (imt,mdiff)            bounda
1055  c     blmc (imt,Nr,mdiff)         boundary layer mixing coefficients  (m^2/s)  c     blmc (imt,Nr,mdiff)         boundary layer mixing coefficients  (m^2/s)
1056  c     ghat (imt,Nr)               nonlocal scalar transport  c     ghat (imt,Nr)               nonlocal scalar transport
1057  c     sigma(imt)                  normalized depth (d / hbl)  c     sigma(imt)                  normalized depth (d / hbl)
1058        _KPP_RL dkm1 (imt,mdiff)        _RL dkm1 (imt,mdiff)
1059        _KPP_RL blmc (imt,Nr,mdiff)        _RL blmc (imt,Nr,mdiff)
1060        _KPP_RL ghat (imt,Nr)        _RL ghat (imt,Nr)
1061        _KPP_RL sigma(imt)        _RL sigma(imt)
1062        integer ikppkey, kkppkey        integer ikppkey, kkppkey
1063    
1064  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
# Line 1113  c  local Line 1067  c  local
1067  c     gat1*(imt)                 shape function at sigma = 1  c     gat1*(imt)                 shape function at sigma = 1
1068  c     dat1*(imt)                 derivative of shape function at sigma = 1  c     dat1*(imt)                 derivative of shape function at sigma = 1
1069  c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)  c     ws(imt), wm(imt)           turbulent velocity scales             (m/s)
1070        _KPP_RL gat1m(imt), gat1s(imt), gat1t(imt)        _RL gat1m(imt), gat1s(imt), gat1t(imt)
1071        _KPP_RL dat1m(imt), dat1s(imt), dat1t(imt)        _RL dat1m(imt), dat1s(imt), dat1t(imt)
1072        _KPP_RL ws(imt), wm(imt)        _RL ws(imt), wm(imt)
1073        integer i, kn, ki        integer i, kn, ki
1074        _KPP_RL R, dvdzup, dvdzdn, viscp        _RL R, dvdzup, dvdzdn, viscp
1075        _KPP_RL difsp, diftp, visch, difsh, difth        _RL difsp, diftp, visch, difsh, difth
1076        _KPP_RL f1, sig, a1, a2, a3, delhat        _RL f1, sig, a1, a2, a3, delhat
1077        _KPP_RL Gm, Gs, Gt        _RL Gm, Gs, Gt
1078        _KPP_RL tempVar        _RL tempVar
1079    
1080        _KPP_RL    p0    , eins        _RL    p0    , eins
1081        parameter (p0=0.0, eins=1.0)        parameter (p0=0.0, eins=1.0)
1082    
1083  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
# Line 1137  c--------------------------------------- Line 1091  c---------------------------------------
1091  CADJ STORE sigma = comlev1_kpp, key = ikppkey  CADJ STORE sigma = comlev1_kpp, key = ikppkey
1092        call wscale (        call wscale (
1093       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1094       O        wm, ws )       O        wm, ws, myThid )
1095  CADJ STORE wm = comlev1_kpp, key = ikppkey  CADJ STORE wm = comlev1_kpp, key = ikppkey
1096  CADJ STORE ws = comlev1_kpp, key = ikppkey  CADJ STORE ws = comlev1_kpp, key = ikppkey
1097    
# Line 1230  CADJ STORE ws = comlev1_kpp_k, key = kkp Line 1184  CADJ STORE ws = comlev1_kpp_k, key = kkp
1184  CADJ STORE sigma = comlev1_kpp_k, key = kkppkey  CADJ STORE sigma = comlev1_kpp_k, key = kkppkey
1185           call wscale (           call wscale (
1186       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1187       O        wm, ws )       O        wm, ws, myThid )
1188  CADJ STORE wm = comlev1_kpp_k, key = kkppkey  CADJ STORE wm = comlev1_kpp_k, key = kkppkey
1189  CADJ STORE ws = comlev1_kpp_k, key = kkppkey  CADJ STORE ws = comlev1_kpp_k, key = kkppkey
1190    
# Line 1283  CADJ STORE ws = comlev1_kpp, key = ikppk Line 1237  CADJ STORE ws = comlev1_kpp, key = ikppk
1237  CADJ STORE sigma = comlev1_kpp, key = ikppkey  CADJ STORE sigma = comlev1_kpp, key = ikppkey
1238        call wscale (        call wscale (
1239       I        sigma, hbl, ustar, bfsfc,       I        sigma, hbl, ustar, bfsfc,
1240       O        wm, ws )       O        wm, ws, myThid )
1241  CADJ STORE wm = comlev1_kpp, key = ikppkey  CADJ STORE wm = comlev1_kpp, key = ikppkey
1242  CADJ STORE ws = comlev1_kpp, key = ikppkey  CADJ STORE ws = comlev1_kpp, key = ikppkey
1243    
# Line 1311  c*************************************** Line 1265  c***************************************
1265       I       dkm1, hbl, kbl, diffus, casea       I       dkm1, hbl, kbl, diffus, casea
1266       U     , ghat       U     , ghat
1267       O     , blmc       O     , blmc
1268       &     )       &     , myThid )
1269    
1270  c enhance the diffusivity at the kbl-.5 interface  c enhance the diffusivity at the kbl-.5 interface
1271    
# Line 1326  c     hbl(imt)                  boundary Line 1280  c     hbl(imt)                  boundary
1280  c     kbl(imt)                  grid above hbl  c     kbl(imt)                  grid above hbl
1281  c     diffus(imt,0:Nrp1,mdiff) vertical diffusivities           (m^2/s)  c     diffus(imt,0:Nrp1,mdiff) vertical diffusivities           (m^2/s)
1282  c     casea(imt)                = 1 in caseA, = 0 in case B  c     casea(imt)                = 1 in caseA, = 0 in case B
1283        _KPP_RL dkm1  (imt,mdiff)  c     myThid                    thread number for this instance of the routine
1284        _KPP_RL hbl   (imt)        integer   myThid
1285          _RL dkm1  (imt,mdiff)
1286          _RL hbl   (imt)
1287        integer kbl   (imt)        integer kbl   (imt)
1288        _KPP_RL diffus(imt,0:Nrp1,mdiff)        _RL diffus(imt,0:Nrp1,mdiff)
1289        _KPP_RL casea (imt)        _RL casea (imt)
1290    
1291  c input/output  c input/output
1292  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)
1293        _KPP_RL ghat (imt,Nr)        _RL ghat (imt,Nr)
1294    
1295  c output  c output
1296  c     enhanced bound. layer mixing coeff.  c     enhanced bound. layer mixing coeff.
1297        _KPP_RL blmc  (imt,Nr,mdiff)        _RL blmc  (imt,Nr,mdiff)
1298    
1299  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1300    
1301  c local  c local
1302  c     fraction hbl lies beteen zgrid neighbors  c     fraction hbl lies beteen zgrid neighbors
1303        _KPP_RL delta        _RL delta
1304        integer ki, i, md        integer ki, i, md
1305        _KPP_RL dkmp5, dstar        _RL dkmp5, dstar
1306    
1307        do i = 1, imt        do i = 1, imt
1308           ki = kbl(i)-1           ki = kbl(i)-1
# Line 1372  c     fraction hbl lies beteen zgrid nei Line 1328  c     fraction hbl lies beteen zgrid nei
1328  c*************************************************************************  c*************************************************************************
1329    
1330        SUBROUTINE STATEKPP (        SUBROUTINE STATEKPP (
1331       I     ikppkey, bi, bj, myThid,       O     RHO1, DBLOC, DBSFC, TTALPHA, SSBETA,
1332       O     RHO1, DBLOC, DBSFC, TTALPHA, SSBETA)       I     ikppkey, bi, bj, myThid )
1333  c  c
1334  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
1335  c     "statekpp" computes all necessary input arrays  c     "statekpp" computes all necessary input arrays
# Line 1414  c--------------------------------------- Line 1370  c---------------------------------------
1370    
1371  c-------------- Routine arguments -----------------------------------------  c-------------- Routine arguments -----------------------------------------
1372        INTEGER bi, bj, myThid        INTEGER bi, bj, myThid
1373        _KPP_RL RHO1   ( ibot:itop, jbot:jtop       )        _RL RHO1   ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy       )
1374        _KPP_RL DBLOC  ( ibot:itop, jbot:jtop, Nr   )        _RL DBLOC  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
1375        _KPP_RL DBSFC  ( ibot:itop, jbot:jtop, Nr   )        _RL DBSFC  ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr   )
1376        _KPP_RL TTALPHA( ibot:itop, jbot:jtop, Nrp1 )        _RL TTALPHA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1377        _KPP_RL SSBETA ( ibot:itop, jbot:jtop, Nrp1 )        _RL SSBETA ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1378    
1379  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
1380    
# Line 1454  CADJ STORE theta(:,:,1,bi,bj) = comlev1_ Line 1410  CADJ STORE theta(:,:,1,bi,bj) = comlev1_
1410  CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey  CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1411  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1412        call FIND_RHO(        call FIND_RHO(
1413       I     bi, bj, ibot, itop, jbot, jtop, 1, 1,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1414       I     theta, salt,       I     theta, salt,
1415       O     WORK1,       O     WORK1,
1416       I     myThid )       I     myThid )
# Line 1464  CADJ STORE salt (:,:,1,bi,bj) = comlev1_ Line 1420  CADJ STORE salt (:,:,1,bi,bj) = comlev1_
1420  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1421    
1422        call FIND_ALPHA(        call FIND_ALPHA(
1423       I     bi, bj, ibot, itop, jbot, jtop, 1, 1,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1424       O     WORK2 )       O     WORK2, myThid )
1425    
1426        call FIND_BETA(        call FIND_BETA(
1427       I     bi, bj, ibot, itop, jbot, jtop, 1, 1,       I     bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1428       O     WORK3 )       O     WORK3, myThid )
1429    
1430        DO J = jbot, jtop        DO J = 1-OLy, sNy+OLy
1431           DO I = ibot, itop           DO I = 1-OLx, sNx+OLx
1432              RHO1(I,J)      = WORK1(I,J) + rhoConst              RHO1(I,J)      = WORK1(I,J) + rhoConst
1433              TTALPHA(I,J,1) = WORK2(I,J)              TTALPHA(I,J,1) = WORK2(I,J)
1434              SSBETA(I,J,1)  = WORK3(I,J)              SSBETA(I,J,1)  = WORK3(I,J)
# Line 1504  CADJ STORE theta(:,:,k,bi,bj) = comlev1_ Line 1460  CADJ STORE theta(:,:,k,bi,bj) = comlev1_
1460  CADJ STORE salt (:,:,k,bi,bj) = comlev1_kpp_k, key=kkppkey  CADJ STORE salt (:,:,k,bi,bj) = comlev1_kpp_k, key=kkppkey
1461  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1462           call FIND_RHO(           call FIND_RHO(
1463       I        bi, bj, ibot, itop, jbot, jtop, K, K,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1464       I        theta, salt,       I        theta, salt,
1465       O        RHOK,       O        RHOK,
1466       I        myThid )       I        myThid )
# Line 1514  CADJ STORE theta(:,:,k-1,bi,bj) = comlev Line 1470  CADJ STORE theta(:,:,k-1,bi,bj) = comlev
1470  CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_kpp_k, key=kkppkey  CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_kpp_k, key=kkppkey
1471  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1472           call FIND_RHO(           call FIND_RHO(
1473       I        bi, bj, ibot, itop, jbot, jtop, K-1, K,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K-1, K,
1474       I        theta, salt,       I        theta, salt,
1475       O        RHOKM1,       O        RHOKM1,
1476       I        myThid )       I        myThid )
# Line 1524  CADJ STORE theta(:,:,1,bi,bj) = comlev1_ Line 1480  CADJ STORE theta(:,:,1,bi,bj) = comlev1_
1480  CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey  CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k, key=kkppkey
1481  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1482           call FIND_RHO(           call FIND_RHO(
1483       I        bi, bj, ibot, itop, jbot, jtop, 1, K,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, K,
1484       I        theta, salt,       I        theta, salt,
1485       O        RHO1K,       O        RHO1K,
1486       I        myThid )       I        myThid )
# Line 1536  CADJ STORE rho1k (:,:)          = comlev Line 1492  CADJ STORE rho1k (:,:)          = comlev
1492  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */  #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1493    
1494           call FIND_ALPHA(           call FIND_ALPHA(
1495       I        bi, bj, ibot, itop, jbot, jtop, K, K,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1496       O        WORK1 )       O        WORK1, myThid )
1497    
1498           call FIND_BETA(           call FIND_BETA(
1499       I        bi, bj, ibot, itop, jbot, jtop, K, K,       I        bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1500       O        WORK2 )       O        WORK2, myThid )
1501    
1502           DO J = jbot, jtop           DO J = 1-OLy, sNy+OLy
1503              DO I = ibot, itop              DO I = 1-OLx, sNx+OLx
1504                 TTALPHA(I,J,K) = WORK1 (I,J)                 TTALPHA(I,J,K) = WORK1 (I,J)
1505                 SSBETA(I,J,K)  = WORK2 (I,J)                 SSBETA(I,J,K)  = WORK2 (I,J)
1506                 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 1585  c     work3 - density of t(1)-.8 & s(1 Line 1541  c     work3 - density of t(1)-.8 & s(1
1541        END DO        END DO
1542    
1543  c     compute arrays for K = Nrp1  c     compute arrays for K = Nrp1
1544        DO J = jbot, jtop        DO J = 1-OLy, sNy+OLy
1545           DO I = ibot, itop           DO I = 1-OLx, sNx+OLx
1546              TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)              TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)
1547              SSBETA(I,J,Nrp1)  = SSBETA(I,J,Nr)              SSBETA(I,J,Nrp1)  = SSBETA(I,J,Nr)
1548              DBLOC(I,J,Nr)     = 0.              DBLOC(I,J,Nr)     = 0.
Legend:
Removed from v.1.23  
changed lines
  Added in v.1.31
  ViewVC Help
Powered by ViewVC 1.1.22