/[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.1 - (hide annotations) (download)
Sun Apr 20 04:03:07 2014 UTC (11 years, 3 months ago) by atn
Branch: MAIN
salt plume volume fraction code mod in progress, installment 01

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

  ViewVC Help
Powered by ViewVC 1.1.22