/[MITgcm]/MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/kpp_calc.F
ViewVC logotype

Annotation of /MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/kpp_calc.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.2 - (hide annotations) (download)
Tue Apr 29 06:49:39 2014 UTC (11 years, 2 months ago) by atn
Branch: MAIN
Changes since 1.1: +4 -1 lines 
in progress:
1. add SALT_PLUME_OPTIONS.h to several files;
2. replace SPalpha (vol) with SPbrineSconst (salinity);
3. move diagnostics outside bi,bj loop.
1 atn 1.2 C $Header: /u/gcmpack/MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/kpp_calc.F,v 1.1 2014/04/20 04:03:07 atn Exp $
2 atn 1.1 C $Name: $
3    
4     #include "KPP_OPTIONS.h"
5 atn 1.2 #ifdef ALLOW_SALT_PLUME
6     #include "SALT_PLUME_OPTIONS.h"
7     #endif
8 atn 1.1
9     CBOP
10     C !ROUTINE: KPP_CALC
11    
12     C !INTERFACE: ==========================================================
13     SUBROUTINE KPP_CALC(
14     I bi, bj, myTime, myIter, myThid )
15    
16     C !DESCRIPTION: \bv
17     C *==========================================================*
18     C | SUBROUTINE KPP_CALC |
19     C | o Compute all KPP fields defined in KPP.h |
20     C *==========================================================*
21     C | This subroutine serves as an interface between MITGCMUV |
22     C | code and NCOM 1-D routines in kpp_routines.F |
23     C *==========================================================*
24     IMPLICIT NONE
25    
26     c=======================================================================
27     c
28     c written by : jan morzel, august 11, 1994
29     c modified by : jan morzel, january 25, 1995 : "dVsq" and 1d code
30     c detlef stammer, august, 1997 : for MIT GCM Classic
31     c d. menemenlis, july, 1998 : for MIT GCM UV
32     c
33     c compute vertical mixing coefficients based on the k-profile
34     c and oceanic planetary boundary layer scheme by large & mcwilliams.
35     c
36     c summary:
37     c - compute interior mixing everywhere:
38     c interior mixing gets computed at all interfaces due to constant
39     c internal wave background activity ("fkpm" and "fkph"), which
40     c is enhanced in places of static instability (local richardson
41     c number < 0).
42     c Additionally, mixing can be enhanced by adding contribution due
43     c to shear instability which is a function of the local richardson
44     c number
45     c - double diffusivity:
46     c interior mixing can be enhanced by double diffusion due to salt
47     c fingering and diffusive convection (ifdef "kmixdd").
48     c - kpp scheme in the boundary layer:
49     c
50     c a.boundary layer depth:
51     c at every gridpoint the depth of the oceanic boundary layer
52     c ("hbl") gets computed by evaluating bulk richardson numbers.
53     c b.boundary layer mixing:
54     c within the boundary layer, above hbl, vertical mixing is
55     c determined by turbulent surface fluxes, and interior mixing at
56     c the lower boundary, i.e. at hbl.
57     c
58     c this subroutine provides the interface between the MITGCM and
59     c the routine "kppmix", where boundary layer depth, vertical
60     c viscosity, vertical diffusivity, and counter gradient term (ghat)
61     c are computed slabwise.
62     c note: subroutine "kppmix" uses m-k-s units.
63     c
64     c time level:
65     c input tracer and velocity profiles are evaluated at time level
66     c tau, surface fluxes come from tau or tau-1.
67     c
68     c grid option:
69     c in this "1-grid" implementation, diffusivity and viscosity
70     c profiles are computed on the "t-grid" (by using velocity shear
71     c profiles averaged from the "u,v-grid" onto the "t-grid"; note, that
72     c the averaging includes zero values on coastal and seafloor grid
73     c points). viscosity on the "u,v-grid" is computed by averaging the
74     c "t-grid" viscosity values onto the "u,v-grid".
75     c
76     c vertical grid:
77     c mixing coefficients get evaluated at the bottom of the lowest
78     c layer, i.e., at depth zw(Nr). these values are only useful when
79     c the model ocean domain does not include the entire ocean down to
80     c the seafloor ("upperocean" setup) and allows flux through the
81     c bottom of the domain. for full-depth runs, these mixing
82     c coefficients are being zeroed out before leaving this subroutine.
83     c
84     c-------------------------------------------------------------------------
85    
86     c global parameters updated by kpp_calc
87     c KPPviscAz - KPP eddy viscosity coefficient (m^2/s)
88     c KPPdiffKzT - KPP diffusion coefficient for temperature (m^2/s)
89     c KPPdiffKzS - KPP diffusion coefficient for salt and tracers (m^2/s)
90     c KPPghat - Nonlocal transport coefficient (s/m^2)
91     c KPPhbl - Boundary layer depth on "t-grid" (m)
92     c KPPfrac - Fraction of short-wave flux penetrating mixing layer
93     c KPPplumefrac- Fraction of saltplume (flux) penetrating mixing layer
94    
95     c-- KPP_CALC computes vertical viscosity and diffusivity for region
96     c (-2:sNx+3,-2:sNy+3) as required by CALC_DIFFUSIVITY and requires
97     c values of uVel, vVel, surfaceForcingU, surfaceForcingV in the
98     c region (-2:sNx+4,-2:sNy+4).
99     c Hence overlap region needs to be set OLx=4, OLy=4.
100     c \ev
101    
102     C !USES: ===============================================================
103     #include "SIZE.h"
104     #include "EEPARAMS.h"
105     #include "PARAMS.h"
106     #include "DYNVARS.h"
107     #include "KPP.h"
108     #include "KPP_PARAMS.h"
109     #include "FFIELDS.h"
110     #include "GRID.h"
111     #include "GAD.h"
112     #ifdef ALLOW_SALT_PLUME
113     # include "SALT_PLUME.h"
114     #endif /* ALLOW_SALT_PLUME */
115     #ifdef ALLOW_SHELFICE
116     # include "SHELFICE.h"
117     #endif /* ALLOW_SHELFICE */
118     #ifdef ALLOW_AUTODIFF_TAMC
119     # include "tamc.h"
120     # include "tamc_keys.h"
121     #else /* ALLOW_AUTODIFF_TAMC */
122     #endif /* ALLOW_AUTODIFF_TAMC */
123    
124     EXTERNAL DIFFERENT_MULTIPLE
125     LOGICAL DIFFERENT_MULTIPLE
126    
127     C !INPUT PARAMETERS: ===================================================
128     c Routine arguments
129     c bi, bj :: Current tile indices
130     c myTime :: Current time in simulation
131     c myIter :: Current iteration number in simulation
132     c myThid :: My Thread Id. number
133    
134     INTEGER bi, bj
135     _RL myTime
136     INTEGER myIter
137     INTEGER myThid
138    
139     #ifdef ALLOW_KPP
140    
141     C !LOCAL VARIABLES: ====================================================
142     c Local constants
143     c minusone, p0, p5, p25, p125, p0625
144     c imin, imax, jmin, jmax - array computation indices
145    
146     _RL minusone
147     parameter( minusone=-1.0)
148     _RL p0 , p5 , p25 , p125 , p0625
149     parameter( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 )
150     integer imin ,imax ,jmin ,jmax
151     parameter(imin=2-OLx,imax=sNx+OLx-1,jmin=2-OLy,jmax=sNy+OLy-1)
152    
153     c Local arrays and variables
154     c work? (nx,ny) - horizontal working arrays
155     c ustar (nx,ny) - surface friction velocity (m/s)
156     c bo (nx,ny) - surface turbulent buoyancy forcing (m^2/s^3)
157     c bosol (nx,ny) - surface radiative buoyancy forcing (m^2/s^3)
158     c boplume(nx,ny) - surface haline buoyancy forcing (m^2/s^3)
159     c shsq (nx,ny,Nr) - local velocity shear squared
160     c at interfaces for ri_iwmix (m^2/s^2)
161     c dVsq (nx,ny,Nr) - velocity shear re surface squared
162     c at grid levels for bldepth (m^2/s^2)
163     c dbloc (nx,ny,Nr) - local delta buoyancy at interfaces
164     c for ri_iwmix and bldepth (m/s^2)
165     c Ritop (nx,ny,Nr) - numerator of bulk richardson number
166     c at grid levels for bldepth
167     c vddiff (nx,ny,Nrp2,1)- vertical viscosity on "t-grid" (m^2/s)
168     c vddiff (nx,ny,Nrp2,2)- vert. diff. on next row for salt&tracers (m^2/s)
169     c vddiff (nx,ny,Nrp2,3)- vert. diff. on next row for temperature (m^2/s)
170     c ghat (nx,ny,Nr) - nonlocal transport coefficient (s/m^2)
171     c hbl (nx,ny) - mixing layer depth (m)
172     c kmtj (nx,ny) - maximum number of wet levels in each column
173     c z0 (nx,ny) - Roughness length (m)
174     c zRef (nx,ny) - Reference depth: Hmix * epsilon (m)
175     c uRef (nx,ny) - Reference zonal velocity (m/s)
176     c vRef (nx,ny) - Reference meridional velocity (m/s)
177    
178     integer work1 ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
179     _RL worka ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
180     _RL work2 ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
181     _RL ustar ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
182     _RL bo ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
183     _RL bosol ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
184     #ifdef ALLOW_SALT_PLUME
185     _RL boplume ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
186     #endif /* ALLOW_SALT_PLUME */
187     _RL shsq ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr )
188     _RL dVsq ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr )
189     _RL dbloc ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr )
190     _RL Ritop ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr )
191     _RL vddiff( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, 0:Nrp1, mdiff )
192     _RL ghat ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr )
193     _RL hbl ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
194     cph(
195     _RL TTALPHA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
196     _RL SSBETA ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
197     cph)
198     #ifdef KPP_ESTIMATE_UREF
199     _RL z0 ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
200     _RL zRef ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
201     _RL uRef ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
202     _RL vRef ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
203     #endif /* KPP_ESTIMATE_UREF */
204    
205     integer i, j, k, kp1, km1, im1, ip1, jm1, jp1
206     integer ikppkey
207    
208     #ifdef KPP_ESTIMATE_UREF
209     _RL tempvar1, dBdz1, dBdz2, ustarX, ustarY
210     #endif
211    
212     #ifdef ALLOW_AUTODIFF_TAMC
213     act1 = bi - myBxLo(myThid)
214     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
215     act2 = bj - myByLo(myThid)
216     max2 = myByHi(myThid) - myByLo(myThid) + 1
217     act3 = myThid - 1
218     max3 = nTx*nTy
219     act4 = ikey_dynamics - 1
220     ikppkey = (act1 + 1) + act2*max1
221     & + act3*max1*max2
222     & + act4*max1*max2*max3
223     #endif /* ALLOW_AUTODIFF_TAMC */
224     CEOP
225    
226     c Check to see if new vertical mixing coefficient should be computed now?
227     IF ( DIFFERENT_MULTIPLE(kpp_freq,myTime,deltaTClock)
228     1 .OR. myTime .EQ. startTime ) THEN
229    
230     c-----------------------------------------------------------------------
231     c prepare input arrays for subroutine "kppmix" to compute
232     c viscosity and diffusivity and ghat.
233     c All input arrays need to be in m-k-s units.
234     c
235     c note: for the computation of the bulk richardson number in the
236     c "bldepth" subroutine, gradients of velocity and buoyancy are
237     c required at every depth. in the case of very fine vertical grids
238     c (thickness of top layer < 2m), the surface reference depth must
239     c be set to zref=epsilon/2*zgrid(k), and the reference value
240     c of velocity and buoyancy must be computed as vertical average
241     c between the surface and 2*zref. in the case of coarse vertical
242     c grids zref is zgrid(1)/2., and the surface reference value is
243     c simply the surface value at zgrid(1).
244     c-----------------------------------------------------------------------
245    
246     c------------------------------------------------------------------------
247     c density related quantities
248     c --------------------------
249     c
250     c work2 - density of surface layer (kg/m^3)
251     c dbloc - local buoyancy gradient at Nr interfaces
252     c g/rho{k+1,k+1} * [ drho{k,k+1}-drho{k+1,k+1} ] (m/s^2)
253     c dbsfc (stored in Ritop to conserve stack memory)
254     c - buoyancy difference with respect to the surface
255     c g * [ drho{1,k}/rho{1,k} - drho{k,k}/rho{k,k} ] (m/s^2)
256     c ttalpha (stored in vddiff(:,:,:,1) to conserve stack memory)
257     c - thermal expansion coefficient without 1/rho factor
258     c d(rho{k,k})/d(T(k)) (kg/m^3/C)
259     c ssbeta (stored in vddiff(:,:,:,2) to conserve stack memory)
260     c - salt expansion coefficient without 1/rho factor
261     c d(rho{k,k})/d(S(k)) (kg/m^3/PSU)
262     c------------------------------------------------------------------------
263    
264     CALL STATEKPP(
265     O work2, dbloc, Ritop,
266     O TTALPHA, SSBETA,
267     I ikppkey, bi, bj, myThid )
268    
269     DO k = 1, Nr
270     DO j = 1-OLy, sNy+OLy
271     DO i = 1-OLx, sNx+OLx
272     ghat(i,j,k) = dbloc(i,j,k)
273     ENDDO
274     ENDDO
275     ENDDO
276    
277     #ifdef KPP_SMOOTH_DBLOC
278     c horizontally smooth dbloc with a 121 filter
279     c smooth dbloc stored in ghat to save space
280     c dbloc(k) is buoyancy gradientnote between k and k+1
281     c levels therefore k+1 mask must be used
282    
283     DO k = 1, Nr-1
284     CALL SMOOTH_HORIZ (
285     I k+1, bi, bj,
286     U ghat (1-OLx,1-OLy,k),
287     I myThid )
288     ENDDO
289    
290     #endif /* KPP_SMOOTH_DBLOC */
291    
292     #ifdef KPP_SMOOTH_DENS
293     c horizontally smooth density related quantities with 121 filters
294     CALL SMOOTH_HORIZ (
295     I 1, bi, bj,
296     U work2,
297     I myThid )
298     DO k = 1, Nr
299     CALL SMOOTH_HORIZ (
300     I k+1, bi, bj,
301     U dbloc (1-OLx,1-OLy,k),
302     I myThid )
303     CALL SMOOTH_HORIZ (
304     I k, bi, bj,
305     U Ritop (1-OLx,1-OLy,k),
306     I myThid )
307     CALL SMOOTH_HORIZ (
308     I k, bi, bj,
309     U TTALPHA(1-OLx,1-OLy,k),
310     I myThid )
311     CALL SMOOTH_HORIZ (
312     I k, bi, bj,
313     U SSBETA(1-OLx,1-OLy,k),
314     I myThid )
315     ENDDO
316     #endif /* KPP_SMOOTH_DENS */
317    
318     DO k = 1, Nr
319     km1 = max(1,k-1)
320     DO j = 1-OLy, sNy+OLy
321     DO i = 1-OLx, sNx+OLx
322    
323     c zero out dbloc over land points (so that the convective
324     c part of the interior mixing can be diagnosed)
325     dbloc(i,j,k) = dbloc(i,j,k) * maskC(i,j,k,bi,bj)
326     & * maskC(i,j,km1,bi,bj)
327     ghat(i,j,k) = ghat(i,j,k) * maskC(i,j,k,bi,bj)
328     & * maskC(i,j,km1,bi,bj)
329     Ritop(i,j,k) = Ritop(i,j,k) * maskC(i,j,k,bi,bj)
330     & * maskC(i,j,km1,bi,bj)
331     if(k.eq.nzmax(i,j,bi,bj)) then
332     dbloc(i,j,k) = p0
333     ghat(i,j,k) = p0
334     Ritop(i,j,k) = p0
335     endif
336    
337     c numerator of bulk richardson number on grid levels
338     c note: land and ocean bottom values need to be set to zero
339     c so that the subroutine "bldepth" works correctly
340     Ritop(i,j,k) = (zgrid(1)-zgrid(k)) * Ritop(i,j,k)
341    
342     ENDDO
343     ENDDO
344     ENDDO
345    
346     cph(
347     cph this avoids a single or double recomp./call of statekpp
348     CADJ store work2 = comlev1_kpp, key = ikppkey
349     #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
350     CADJ store dbloc, Ritop, ghat = comlev1_kpp, key = ikppkey
351     CADJ store vddiff = comlev1_kpp, key = ikppkey
352     CADJ store TTALPHA, SSBETA = comlev1_kpp, key = ikppkey
353     #endif
354     cph)
355    
356     CML#ifdef ALLOW_SHELFICE
357     CMLC For the pbl parameterisation to work underneath the ice shelves
358     CMLC it needs to know the surface (ice-ocean) fluxes. However, masking
359     CMLC and indexing problems make this part of the code not work
360     CMLC underneath the ice shelves and the following lines are only here
361     CMLC to remind me that this still needs to be sorted out.
362     CML shelfIceFac = 0. _d 0
363     CML IF ( useShelfIce ) selfIceFac = 1. _d 0
364     CML DO j = jmin, jmax
365     CML DO i = imin, imax
366     CML surfForcT = surfaceForcingT(i,j,bi,bj)
367     CML & + shelficeForcingT(i,j,bi,bj) * shelfIceFac
368     CML surfForcS = surfaceForcingS(i,j,bi,bj)
369     CML & + shelficeForcingS(i,j,bi,bj) * shelfIceFac
370     CML ENDDO
371     CML ENDDO
372     CML#endif /* ALLOW_SHELFICE */
373    
374     c------------------------------------------------------------------------
375     c friction velocity, turbulent and radiative surface buoyancy forcing
376     c -------------------------------------------------------------------
377     c taux / rho = surfaceForcingU (N/m^2)
378     c tauy / rho = surfaceForcingV (N/m^2)
379     c ustar = sqrt( sqrt( taux^2 + tauy^2 ) / rho ) (m/s)
380     c bo = - g * ( alpha*surfaceForcingT +
381     c beta *surfaceForcingS ) / rho (m^2/s^3)
382     c bosol = - g * alpha * Qsw * drF(1) / rho (m^2/s^3)
383     c boplume = g * (beta * saltPlumeFlux/rhoConst ) /rho (m^2/s^3)
384     c------------------------------------------------------------------------
385     c velocity shear
386     c --------------
387     c Get velocity shear squared, averaged from "u,v-grid"
388     c onto "t-grid" (in (m/s)**2):
389     c dVsq(k)=(Uref-U(k))**2+(Vref-V(k))**2 at grid levels
390     c shsq(k)=(U(k)-U(k+1))**2+(V(k)-V(k+1))**2 at interfaces
391     c
392     c note: Vref can depend on the surface fluxes that is why we compute
393     c dVsq in the subroutine that does the surface related stuff
394     c (admittedly this is a bit messy)
395     c------------------------------------------------------------------------
396    
397     CALL KPP_FORCING_SURF(
398     I work2, surfaceForcingU, surfaceForcingV,
399     I surfaceForcingT, surfaceForcingS, surfaceForcingTice,
400     I Qsw,
401     #ifdef ALLOW_SALT_PLUME
402     #ifndef SALT_PLUME_VOLUME
403     I saltPlumeFlux,
404     #endif /* SALT_PLUME_VOLUME */
405     #endif /* ALLOW_SALT_PLUME */
406     I ttalpha, ssbeta,
407     O ustar, bo, bosol,
408     #ifdef ALLOW_SALT_PLUME
409     #ifndef SALT_PLUME_VOLUME
410     O boplume,
411     #endif /* SALT_PLUME_VOLUME */
412     #endif /* ALLOW_SALT_PLUME */
413     O dVsq,
414     I ikppkey, iMin, iMax, jMin, jMax, bi, bj, myTime, myThid )
415    
416     CMLcph(
417     CMLCADJ store ustar = comlev1_kpp, key = ikppkey
418     CMLcph)
419    
420     c initialize arrays to zero
421     DO k = 1, Nr
422     DO j = 1-OLy, sNy+OLy
423     DO i = 1-OLx, sNx+OLx
424     shsq(i,j,k) = p0
425     ENDDO
426     ENDDO
427     ENDDO
428    
429     c shsq computation
430     DO k = 1, Nrm1
431     kp1 = k + 1
432     DO j = jmin, jmax
433     jm1 = j - 1
434     jp1 = j + 1
435     DO i = imin, imax
436     im1 = i - 1
437     ip1 = i + 1
438     shsq(i,j,k) = p5 * (
439     & (uVel(i, j, k,bi,bj)-uVel(i, j, kp1,bi,bj)) *
440     & (uVel(i, j, k,bi,bj)-uVel(i, j, kp1,bi,bj)) +
441     & (uVel(ip1,j, k,bi,bj)-uVel(ip1,j, kp1,bi,bj)) *
442     & (uVel(ip1,j, k,bi,bj)-uVel(ip1,j, kp1,bi,bj)) +
443     & (vVel(i, j, k,bi,bj)-vVel(i, j, kp1,bi,bj)) *
444     & (vVel(i, j, k,bi,bj)-vVel(i, j, kp1,bi,bj)) +
445     & (vVel(i, jp1,k,bi,bj)-vVel(i, jp1,kp1,bi,bj)) *
446     & (vVel(i, jp1,k,bi,bj)-vVel(i, jp1,kp1,bi,bj)) )
447     #ifdef KPP_SMOOTH_SHSQ
448     shsq(i,j,k) = p5 * shsq(i,j,k) + p125 * (
449     & (uVel(i, jm1,k,bi,bj)-uVel(i, jm1,kp1,bi,bj)) *
450     & (uVel(i, jm1,k,bi,bj)-uVel(i, jm1,kp1,bi,bj)) +
451     & (uVel(ip1,jm1,k,bi,bj)-uVel(ip1,jm1,kp1,bi,bj)) *
452     & (uVel(ip1,jm1,k,bi,bj)-uVel(ip1,jm1,kp1,bi,bj)) +
453     & (uVel(i, jp1,k,bi,bj)-uVel(i, jp1,kp1,bi,bj)) *
454     & (uVel(i, jp1,k,bi,bj)-uVel(i, jp1,kp1,bi,bj)) +
455     & (uVel(ip1,jp1,k,bi,bj)-uVel(ip1,jp1,kp1,bi,bj)) *
456     & (uVel(ip1,jp1,k,bi,bj)-uVel(ip1,jp1,kp1,bi,bj)) +
457     & (vVel(im1,j, k,bi,bj)-vVel(im1,j, kp1,bi,bj)) *
458     & (vVel(im1,j, k,bi,bj)-vVel(im1,j, kp1,bi,bj)) +
459     & (vVel(im1,jp1,k,bi,bj)-vVel(im1,jp1,kp1,bi,bj)) *
460     & (vVel(im1,jp1,k,bi,bj)-vVel(im1,jp1,kp1,bi,bj)) +
461     & (vVel(ip1,j, k,bi,bj)-vVel(ip1,j, kp1,bi,bj)) *
462     & (vVel(ip1,j, k,bi,bj)-vVel(ip1,j, kp1,bi,bj)) +
463     & (vVel(ip1,jp1,k,bi,bj)-vVel(ip1,jp1,kp1,bi,bj)) *
464     & (vVel(ip1,jp1,k,bi,bj)-vVel(ip1,jp1,kp1,bi,bj)) )
465     #endif
466     ENDDO
467     ENDDO
468     ENDDO
469    
470     cph(
471     #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
472     CADJ store dvsq, shsq = comlev1_kpp, key = ikppkey
473     #endif
474     cph)
475    
476     c-----------------------------------------------------------------------
477     c solve for viscosity, diffusivity, ghat, and hbl on "t-grid"
478     c-----------------------------------------------------------------------
479    
480     c precompute background vertical diffusivities, which are needed for
481     c matching diffusivities at bottom of KPP PBL
482     CALL CALC_3D_DIFFUSIVITY(
483     I bi,bj,1-Olx,sNx+OLx,1-Oly,sNy+OLy,
484     I GAD_SALINITY, .FALSE., .FALSE.,
485     O KPPdiffKzS(1-Olx,1-Oly,1,bi,bj),
486     I myThid)
487     CALL CALC_3D_DIFFUSIVITY(
488     I bi,bj,1-Olx,sNx+OLx,1-Oly,sNy+OLy,
489     I GAD_TEMPERATURE, .FALSE., .FALSE.,
490     O KPPdiffKzT(1-Olx,1-Oly,1,bi,bj),
491     I myThid)
492     #ifndef EXCLUDE_KPP_DOUBLEDIFF
493     IF ( KPPuseDoubleDiff ) THEN
494     C Add the contribution of double diffusive effects (salt fingering
495     C and diffusive convection) here. It would be more logical to add
496     C them right after Ri_iwmix within kppmix, but ttalpha, ssbeta, theta
497     C and salt are not passed to kppmix and are thus not available there.
498     CALL KPP_DOUBLEDIFF(
499     I TTALPHA, SSBETA,
500     U KPPdiffKzT(1-Olx,1-Oly,1,bi,bj),
501     U KPPdiffKzS(1-Olx,1-Oly,1,bi,bj),
502     I ikppkey,1-Olx,sNx+OLx,1-Oly,sNy+OLy,bi,bj,myThid)
503     ENDIF
504     #endif /* ndef EXCLUDE_KPP_DOUBLEDIFF */
505    
506     DO j = 1-OLy, sNy+OLy
507     DO i = 1-OLx, sNx+OLx
508     work1(i,j) = nzmax(i,j,bi,bj)
509     work2(i,j) = Fcori(i,j,bi,bj)
510     ENDDO
511     ENDDO
512     CALL KPPMIX (
513     I work1, shsq, dVsq, ustar
514     I , maskC(1-Olx,1-Oly,1,bi,bj)
515     I , bo, bosol
516     #ifdef ALLOW_SALT_PLUME
517     #ifndef SALT_PLUME_VOLUME
518     I , boplume, SaltPlumeDepth(1-Olx,1-Oly,bi,bj)
519     #endif /* SALT_PLUME_VOLUME */
520     #endif /* ALLOW_SALT_PLUME */
521     I , dbloc, Ritop, work2
522     I , KPPdiffKzS(1-Olx,1-Oly,1,bi,bj)
523     I , KPPdiffKzT(1-Olx,1-Oly,1,bi,bj)
524     I , ikppkey
525     O , vddiff
526     U , ghat
527     O , hbl
528     I , bi, bj, mytime, myIter, mythid )
529    
530     c-----------------------------------------------------------------------
531     c zero out land values and transfer to global variables
532     c-----------------------------------------------------------------------
533    
534     DO j = jmin, jmax
535     DO i = imin, imax
536     DO k = 1, Nr
537     km1 = max(1,k-1)
538     KPPviscAz(i,j,k,bi,bj) = vddiff(i,j,k-1,1) * maskC(i,j,k,bi,bj)
539     & * maskC(i,j,km1,bi,bj)
540     KPPdiffKzS(i,j,k,bi,bj)= vddiff(i,j,k-1,2) * maskC(i,j,k,bi,bj)
541     & * maskC(i,j,km1,bi,bj)
542     KPPdiffKzT(i,j,k,bi,bj)= vddiff(i,j,k-1,3) * maskC(i,j,k,bi,bj)
543     & * maskC(i,j,km1,bi,bj)
544     KPPghat(i,j,k,bi,bj) = ghat(i,j,k) * maskC(i,j,k,bi,bj)
545     & * maskC(i,j,km1,bi,bj)
546     ENDDO
547     k = 1
548     #ifdef ALLOW_SHELFICE
549     if ( useShelfIce ) k = kTopC(i,j,bi,bj)
550     #endif /* ALLOW_SHELFICE */
551     KPPhbl(i,j,bi,bj) = hbl(i,j) * maskC(i,j,k,bi,bj)
552    
553     ENDDO
554     ENDDO
555    
556     #ifdef KPP_SMOOTH_VISC
557     c horizontal smoothing of vertical viscosity
558     DO k = 1, Nr
559     CALL SMOOTH_HORIZ (
560     I k, bi, bj,
561     U KPPviscAz(1-OLx,1-OLy,k,bi,bj),
562     I myThid )
563     ENDDO
564     C jmc: No EXCH inside bi,bj loop !!!
565     c _EXCH_XYZ_RL(KPPviscAz , myThid )
566     #endif /* KPP_SMOOTH_VISC */
567    
568     #ifdef KPP_SMOOTH_DIFF
569     c horizontal smoothing of vertical diffusivity
570     DO k = 1, Nr
571     CALL SMOOTH_HORIZ (
572     I k, bi, bj,
573     U KPPdiffKzS(1-OLx,1-OLy,k,bi,bj),
574     I myThid )
575     CALL SMOOTH_HORIZ (
576     I k, bi, bj,
577     U KPPdiffKzT(1-OLx,1-OLy,k,bi,bj),
578     I myThid )
579     ENDDO
580     #endif /* KPP_SMOOTH_DIFF */
581    
582     cph(
583     cph crucial: this avoids full recomp./call of kppmix
584     CADJ store KPPhbl = comlev1_kpp, key = ikppkey
585     cph)
586    
587     C Compute fraction of solar short-wave flux penetrating to
588     C the bottom of the mixing layer.
589     DO j=1-OLy,sNy+OLy
590     DO i=1-OLx,sNx+OLx
591     worka(i,j) = KPPhbl(i,j,bi,bj)
592     ENDDO
593     ENDDO
594     CALL SWFRAC(
595     I (sNx+2*OLx)*(sNy+2*OLy), minusone,
596     U worka,
597     I myTime, myIter, myThid )
598     DO j=1-OLy,sNy+OLy
599     DO i=1-OLx,sNx+OLx
600     KPPfrac(i,j,bi,bj) = worka(i,j)
601     ENDDO
602     ENDDO
603    
604     #ifdef ALLOW_SALT_PLUME
605     #ifndef SALT_PLUME_VOLUME
606     C Compute fraction of saltplume (flux) penetrating to
607     C the bottom of the mixing layer.
608     IF ( useSALT_PLUME ) THEN
609     DO j=1-OLy,sNy+OLy
610     DO i=1-OLx,sNx+OLx
611     work2(i,j) = SaltPlumeDepth(i,j,bi,bj)
612     worka(i,j) = KPPhbl(i,j,bi,bj)
613     ENDDO
614     ENDDO
615     CALL SALT_PLUME_FRAC(
616     I (sNx+2*OLx)*(sNy+2*OLy), minusone, work2,
617     U worka,
618     I myTime, myIter, myThid )
619     DO j=1-OLy,sNy+OLy
620     DO i=1-OLx,sNx+OLx
621     KPPplumefrac(i,j,bi,bj) = 1. _d 0 - worka(i,j)
622     ENDDO
623     ENDDO
624     ENDIF
625     #endif /* ndef SALT_PLUME_VOLUME */
626     #endif /* ALLOW_SALT_PLUME */
627    
628     ENDIF
629    
630     #endif /* ALLOW_KPP */
631    
632     RETURN
633     END
634    
635     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
636    
637     SUBROUTINE KPP_CALC_DUMMY(
638     I bi, bj, myTime, myIter, myThid )
639     C *==========================================================*
640     C | SUBROUTINE KPP_CALC_DUMMY |
641     C | o Compute all KPP fields defined in KPP.h |
642     C | o Dummy routine for TAMC
643     C *==========================================================*
644     C | This subroutine serves as an interface between MITGCMUV |
645     C | code and NCOM 1-D routines in kpp_routines.F |
646     C *==========================================================*
647     IMPLICIT NONE
648    
649     #include "SIZE.h"
650     #include "EEPARAMS.h"
651     #include "PARAMS.h"
652     #include "KPP.h"
653     #include "KPP_PARAMS.h"
654     #include "GRID.h"
655     #include "GAD.h"
656    
657     c Routine arguments
658     c bi, bj :: Current tile indices
659     c myTime :: Current time in simulation
660     c myIter :: Current iteration number in simulation
661     c myThid :: My Thread Id. number
662    
663     INTEGER bi, bj
664     _RL myTime
665     INTEGER myIter
666     INTEGER myThid
667    
668     #ifdef ALLOW_KPP
669    
670     c Local constants
671     integer i, j, k
672    
673     DO j=1-OLy,sNy+OLy
674     DO i=1-OLx,sNx+OLx
675     KPPhbl (i,j,bi,bj) = 1.0
676     KPPfrac(i,j,bi,bj) = 0.0
677     #ifdef ALLOW_SALT_PLUME
678     KPPplumefrac(i,j,bi,bj) = 0.0
679     #endif /* ALLOW_SALT_PLUME */
680     DO k = 1,Nr
681     KPPghat (i,j,k,bi,bj) = 0.0
682     KPPviscAz (i,j,k,bi,bj) = viscArNr(1)
683     ENDDO
684     ENDDO
685     ENDDO
686    
687     CALL CALC_3D_DIFFUSIVITY(
688     I bi,bj,1-Olx,sNx+OLx,1-Oly,sNy+OLy,
689     I GAD_SALINITY, .FALSE., .FALSE.,
690     O KPPdiffKzS(1-Olx,1-Oly,1,bi,bj),
691     I myThid)
692     CALL CALC_3D_DIFFUSIVITY(
693     I bi,bj,1-Olx,sNx+OLx,1-Oly,sNy+OLy,
694     I GAD_TEMPERATURE, .FALSE., .FALSE.,
695     O KPPdiffKzT(1-Olx,1-Oly,1,bi,bj),
696     I myThid)
697    
698     #endif
699     RETURN
700     END
  ViewVC Help
Powered by ViewVC 1.1.22