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

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

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


Revision 1.4 - (hide annotations) (download)
Thu May 1 21:30:48 2014 UTC (11 years, 3 months ago) by atn
Branch: MAIN
Changes since 1.3: +11 -11 lines
bug fixing

1 atn 1.4 C $Header: /u/gcmpack/MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/kpp_routines.F,v 1.3 2014/05/01 08:09:30 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     C-- File kpp_routines.F: subroutines needed to implement
10     C-- KPP vertical mixing scheme
11     C-- Contents
12     C-- o KPPMIX - Main driver and interface routine.
13     C-- o BLDEPTH - Determine oceanic planetary boundary layer depth.
14     C-- o WSCALE - Compute turbulent velocity scales.
15     C-- o RI_IWMIX - Compute interior viscosity diffusivity coefficients.
16     C-- o Z121 - Apply 121 vertical smoothing.
17     C-- o SMOOTH_HORIZ- Apply horizontal smoothing to global array.
18     C-- o BLMIX - Boundary layer mixing coefficients.
19     C-- o ENHANCE - Enhance diffusivity at boundary layer interface.
20     C-- o STATEKPP - Compute buoyancy-related input arrays.
21     C-- o KPP_DOUBLEDIFF - Compute double diffusive contribution to diffusivities
22    
23     c***********************************************************************
24    
25     SUBROUTINE KPPMIX (
26     I kmtj, shsq, dvsq, ustar, msk
27     I , bo, bosol
28     #ifdef ALLOW_SALT_PLUME
29     I , boplume,SPDepth
30     #endif /* ALLOW_SALT_PLUME */
31     I , dbloc, Ritop, coriol
32     I , diffusKzS, diffusKzT
33     I , ikppkey
34     O , diffus
35     U , ghat
36     O , hbl
37     I , bi, bj, myTime, myIter, myThid )
38    
39     c-----------------------------------------------------------------------
40     c
41     c Main driver subroutine for kpp vertical mixing scheme and
42     c interface to greater ocean model
43     c
44     c written by: bill large, june 6, 1994
45     c modified by: jan morzel, june 30, 1994
46     c bill large, august 11, 1994
47     c bill large, january 25, 1995 : "dVsq" and 1d code
48     c detlef stammer, august 1997 : for use with MIT GCM Classic
49     c d. menemenlis, june 1998 : for use with MIT GCM UV
50     c
51     c-----------------------------------------------------------------------
52    
53     IMPLICIT NONE
54    
55     #include "SIZE.h"
56     #include "EEPARAMS.h"
57     #include "PARAMS.h"
58     #include "KPP_PARAMS.h"
59     #ifdef ALLOW_AUTODIFF
60     # include "tamc.h"
61     #endif
62    
63     c input
64     c bi, bj :: Array indices on which to apply calculations
65     c myTime :: Current time in simulation
66     c myIter :: Current iteration number in simulation
67     c myThid :: My Thread Id. number
68     c kmtj (imt) - number of vertical layers on this row
69     c msk (imt) - surface mask (=1 if water, =0 otherwise)
70     c shsq (imt,Nr) - (local velocity shear)^2 ((m/s)^2)
71     c dvsq (imt,Nr) - (velocity shear re sfc)^2 ((m/s)^2)
72     c ustar (imt) - surface friction velocity (m/s)
73     c bo (imt) - surface turbulent buoy. forcing (m^2/s^3)
74     c bosol (imt) - radiative buoyancy forcing (m^2/s^3)
75 atn 1.4 c boplume(imt,Nr) - haline buoyancy forcing (m^2/s^3)
76 atn 1.1 c dbloc (imt,Nr) - local delta buoyancy across interfaces (m/s^2)
77     c dblocSm(imt,Nr) - horizontally smoothed dbloc (m/s^2)
78     c stored in ghat to save space
79     c Ritop (imt,Nr) - numerator of bulk Richardson Number
80     c (zref-z) * delta buoyancy w.r.t. surface ((m/s)^2)
81     c coriol (imt) - Coriolis parameter (1/s)
82     c diffusKzS(imt,Nr)- background vertical diffusivity for scalars (m^2/s)
83     c diffusKzT(imt,Nr)- background vertical diffusivity for theta (m^2/s)
84     c note: there is a conversion from 2-D to 1-D for input output variables,
85     c e.g., hbl(sNx,sNy) -> hbl(imt),
86     c where hbl(i,j) -> hbl((j-1)*sNx+i)
87     INTEGER bi, bj
88     _RL myTime
89     integer myIter
90     integer myThid
91     integer kmtj (imt )
92     _RL shsq (imt,Nr)
93     _RL dvsq (imt,Nr)
94     _RL ustar (imt )
95     _RL bo (imt )
96     _RL bosol (imt )
97     #ifdef ALLOW_SALT_PLUME
98 atn 1.4 _RL boplume (imt,Nr)
99 atn 1.1 _RL SPDepth (imt )
100     #endif /* ALLOW_SALT_PLUME */
101     _RL dbloc (imt,Nr)
102     _RL Ritop (imt,Nr)
103     _RL coriol (imt )
104     _RS msk (imt )
105     _RL diffusKzS(imt,Nr)
106     _RL diffusKzT(imt,Nr)
107    
108     integer ikppkey
109    
110     c output
111     c diffus (imt,1) - vertical viscosity coefficient (m^2/s)
112     c diffus (imt,2) - vertical scalar diffusivity (m^2/s)
113     c diffus (imt,3) - vertical temperature diffusivity (m^2/s)
114     c ghat (imt) - nonlocal transport coefficient (s/m^2)
115     c hbl (imt) - mixing layer depth (m)
116    
117     _RL diffus(imt,0:Nrp1,mdiff)
118     _RL ghat (imt,Nr)
119     _RL hbl (imt)
120    
121     #ifdef ALLOW_KPP
122    
123     c local
124     c kbl (imt ) - index of first grid level below hbl
125     c bfsfc (imt ) - surface buoyancy forcing (m^2/s^3)
126     c casea (imt ) - 1 in case A; 0 in case B
127     c stable (imt ) - 1 in stable forcing; 0 if unstable
128     c dkm1 (imt, mdiff) - boundary layer diffusivity at kbl-1 level
129     c blmc (imt,Nr,mdiff) - boundary layer mixing coefficients
130     c sigma (imt ) - normalized depth (d / hbl)
131     c Rib (imt,Nr ) - bulk Richardson number
132    
133     integer kbl(imt )
134     _RL bfsfc (imt )
135     _RL casea (imt )
136     _RL stable (imt )
137     _RL dkm1 (imt, mdiff)
138     _RL blmc (imt,Nr,mdiff)
139     _RL sigma (imt )
140     _RL Rib (imt,Nr )
141    
142     integer i, k, md
143    
144     c-----------------------------------------------------------------------
145     c compute interior mixing coefficients everywhere, due to constant
146     c internal wave activity, static instability, and local shear
147     c instability.
148     c (ghat is temporary storage for horizontally smoothed dbloc)
149     c-----------------------------------------------------------------------
150    
151     cph(
152     cph these storings avoid recomp. of Ri_iwmix
153     CADJ STORE ghat = comlev1_kpp, key=ikppkey, kind=isbyte
154     CADJ STORE dbloc = comlev1_kpp, key=ikppkey, kind=isbyte
155     cph)
156     call Ri_iwmix (
157     I kmtj, shsq, dbloc, ghat
158     I , diffusKzS, diffusKzT
159     I , ikppkey
160     O , diffus, myThid )
161    
162     cph(
163     cph these storings avoid recomp. of Ri_iwmix
164     cph DESPITE TAFs 'not necessary' warning!
165     CADJ STORE dbloc = comlev1_kpp, key=ikppkey, kind=isbyte
166     CADJ STORE shsq = comlev1_kpp, key=ikppkey, kind=isbyte
167     CADJ STORE ghat = comlev1_kpp, key=ikppkey, kind=isbyte
168     CADJ STORE diffus = comlev1_kpp, key=ikppkey, kind=isbyte
169     cph)
170    
171     c-----------------------------------------------------------------------
172     c set seafloor values to zero and fill extra "Nrp1" coefficients
173     c for blmix
174     c-----------------------------------------------------------------------
175    
176     do md = 1, mdiff
177     do k=1,Nrp1
178     do i = 1,imt
179     if(k.ge.kmtj(i)) diffus(i,k,md) = 0.0
180     end do
181     end do
182     end do
183    
184     c-----------------------------------------------------------------------
185     c compute boundary layer mixing coefficients:
186     c
187     c diagnose the new boundary layer depth
188     c-----------------------------------------------------------------------
189    
190     call bldepth (
191     I kmtj
192     I , dvsq, dbloc, Ritop, ustar, bo, bosol
193     #ifdef ALLOW_SALT_PLUME
194     I , boplume,SPDepth
195     #endif /* ALLOW_SALT_PLUME */
196     I , coriol
197     I , ikppkey
198     O , hbl, bfsfc, stable, casea, kbl, Rib, sigma
199     I , bi, bj, myTime, myIter, myThid )
200    
201     CADJ STORE hbl,bfsfc,stable,casea,kbl = comlev1_kpp,
202     CADJ & key=ikppkey, kind=isbyte
203    
204     c-----------------------------------------------------------------------
205     c compute boundary layer diffusivities
206     c-----------------------------------------------------------------------
207    
208     call blmix (
209     I ustar, bfsfc, hbl, stable, casea, diffus, kbl
210     O , dkm1, blmc, ghat, sigma, ikppkey
211     I , myThid )
212     cph(
213     CADJ STORE dkm1,blmc,ghat = comlev1_kpp,
214     CADJ & key=ikppkey, kind=isbyte
215     CADJ STORE hbl, kbl, diffus, casea = comlev1_kpp,
216     CADJ & key=ikppkey, kind=isbyte
217     cph)
218    
219     c-----------------------------------------------------------------------
220     c enhance diffusivity at interface kbl - 1
221     c-----------------------------------------------------------------------
222    
223     call enhance (
224     I dkm1, hbl, kbl, diffus, casea
225     U , ghat
226     O , blmc
227     I , myThid )
228    
229     cph(
230     cph avoids recomp. of enhance
231     CADJ STORE blmc = comlev1_kpp, key=ikppkey, kind=isbyte
232     cph)
233    
234     c-----------------------------------------------------------------------
235     c combine interior and boundary layer coefficients and nonlocal term
236     c !!!NOTE!!! In shallow (2-level) regions and for shallow mixed layers
237     c (< 1 level), diffusivity blmc can become negative. The max-s below
238     c are a hack until this problem is properly diagnosed and fixed.
239     c-----------------------------------------------------------------------
240     do k = 1, Nr
241     do i = 1, imt
242     if (k .lt. kbl(i)) then
243     #ifdef ALLOW_SHELFICE
244     C when there is shelfice on top (msk(i)=0), reset the boundary layer
245     C mixing coefficients blmc to pure Ri-number based mixing
246     blmc(i,k,1) = max ( blmc(i,k,1)*msk(i),
247     & diffus(i,k,1) )
248     blmc(i,k,2) = max ( blmc(i,k,2)*msk(i),
249     & diffus(i,k,2) )
250     blmc(i,k,3) = max ( blmc(i,k,3)*msk(i),
251     & diffus(i,k,3) )
252     #endif /* not ALLOW_SHELFICE */
253     diffus(i,k,1) = max ( blmc(i,k,1), viscArNr(1) )
254     diffus(i,k,2) = max ( blmc(i,k,2), diffusKzS(i,Nr) )
255     diffus(i,k,3) = max ( blmc(i,k,3), diffusKzT(i,Nr) )
256     else
257     ghat(i,k) = 0. _d 0
258     endif
259     end do
260     end do
261    
262     #endif /* ALLOW_KPP */
263    
264     return
265     end
266    
267     c*************************************************************************
268    
269     subroutine bldepth (
270     I kmtj
271     I , dvsq, dbloc, Ritop, ustar, bo, bosol
272     #ifdef ALLOW_SALT_PLUME
273     I , boplume,SPDepth
274     #endif /* ALLOW_SALT_PLUME */
275     I , coriol
276     I , ikppkey
277     O , hbl, bfsfc, stable, casea, kbl, Rib, sigma
278     I , bi, bj, myTime, myIter, myThid )
279    
280     c the oceanic planetary boundary layer depth, hbl, is determined as
281     c the shallowest depth where the bulk Richardson number is
282     c equal to the critical value, Ricr.
283     c
284     c bulk Richardson numbers are evaluated by computing velocity and
285     c buoyancy differences between values at zgrid(kl) < 0 and surface
286     c reference values.
287     c in this configuration, the reference values are equal to the
288     c values in the surface layer.
289     c when using a very fine vertical grid, these values should be
290     c computed as the vertical average of velocity and buoyancy from
291     c the surface down to epsilon*zgrid(kl).
292     c
293     c when the bulk Richardson number at k exceeds Ricr, hbl is
294     c linearly interpolated between grid levels zgrid(k) and zgrid(k-1).
295     c
296     c The water column and the surface forcing are diagnosed for
297     c stable/ustable forcing conditions, and where hbl is relative
298     c to grid points (caseA), so that conditional branches can be
299     c avoided in later subroutines.
300     c
301     IMPLICIT NONE
302    
303     #include "SIZE.h"
304     #include "EEPARAMS.h"
305     #include "PARAMS.h"
306     #include "KPP_PARAMS.h"
307     #ifdef ALLOW_AUTODIFF
308     # include "tamc.h"
309     #endif
310    
311     c input
312     c------
313     c bi, bj :: Array indices on which to apply calculations
314     c myTime :: Current time in simulation
315     c myIter :: Current iteration number in simulation
316     c myThid :: My Thread Id. number
317     c kmtj : number of vertical layers
318     c dvsq : (velocity shear re sfc)^2 ((m/s)^2)
319     c dbloc : local delta buoyancy across interfaces (m/s^2)
320     c Ritop : numerator of bulk Richardson Number
321     c =(z-zref)*dbsfc, where dbsfc=delta
322     c buoyancy with respect to surface ((m/s)^2)
323     c ustar : surface friction velocity (m/s)
324     c bo : surface turbulent buoyancy forcing (m^2/s^3)
325     c bosol : radiative buoyancy forcing (m^2/s^3)
326     c boplume : haline buoyancy forcing (m^2/s^3)
327     c coriol : Coriolis parameter (1/s)
328     INTEGER bi, bj
329     _RL myTime
330     integer myIter
331     integer myThid
332     integer kmtj(imt)
333     _RL dvsq (imt,Nr)
334     _RL dbloc (imt,Nr)
335     _RL Ritop (imt,Nr)
336     _RL ustar (imt)
337     _RL bo (imt)
338     _RL bosol (imt)
339     _RL coriol (imt)
340     integer ikppkey
341     #ifdef ALLOW_SALT_PLUME
342 atn 1.4 _RL boplume (imt,Nr)
343 atn 1.1 _RL SPDepth (imt)
344     #endif /* ALLOW_SALT_PLUME */
345    
346     c output
347     c--------
348     c hbl : boundary layer depth (m)
349     c bfsfc : Bo+radiation absorbed to d=hbf*hbl (m^2/s^3)
350     c stable : =1 in stable forcing; =0 unstable
351     c casea : =1 in case A, =0 in case B
352     c kbl : -1 of first grid level below hbl
353     c Rib : Bulk Richardson number
354     c sigma : normalized depth (d/hbl)
355     _RL hbl (imt)
356     _RL bfsfc (imt)
357     _RL stable (imt)
358     _RL casea (imt)
359     integer kbl(imt)
360     _RL Rib (imt,Nr)
361     _RL sigma (imt)
362    
363     #ifdef ALLOW_KPP
364    
365     c local
366     c-------
367     c wm, ws : turbulent velocity scales (m/s)
368     _RL wm(imt), ws(imt)
369     _RL worka(imt)
370     _RL bvsq, vtsq, hekman, hmonob, hlimit, tempVar1, tempVar2
371 atn 1.4 integer i, k, kl
372 atn 1.1
373     _RL p5 , eins
374     parameter ( p5=0.5, eins=1.0 )
375     _RL minusone
376     parameter ( minusone=-1.0 )
377     #ifdef ALLOW_AUTODIFF_TAMC
378     integer kkppkey
379     #endif
380    
381     #ifdef ALLOW_DIAGNOSTICS
382     c KPPBFSFC - Bo+radiation absorbed to d=hbf*hbl + plume (m^2/s^3)
383     _RL KPPBFSFC(imt,Nr)
384     _RL KPPRi(imt,Nr)
385     #endif /* ALLOW_DIAGNOSTICS */
386    
387     c find bulk Richardson number at every grid level until > Ricr
388     c
389     c note: the reference depth is -epsilon/2.*zgrid(k), but the reference
390     c u,v,t,s values are simply the surface layer values,
391     c and not the averaged values from 0 to 2*ref.depth,
392     c which is necessary for very fine grids(top layer < 2m thickness)
393     c note: max values when Ricr never satisfied are
394     c kbl(i)=kmtj(i) and hbl(i)=-zgrid(kmtj(i))
395    
396     c initialize hbl and kbl to bottomed out values
397    
398     do i = 1, imt
399     Rib(i,1) = 0. _d 0
400     kbl(i) = max(kmtj(i),1)
401     hbl(i) = -zgrid(kbl(i))
402     end do
403    
404     #ifdef ALLOW_DIAGNOSTICS
405     do kl = 1, Nr
406     do i = 1, imt
407     KPPBFSFC(i,kl) = 0. _d 0
408     KPPRi(i,kl) = 0. _d 0
409     enddo
410     enddo
411     #endif /* ALLOW_DIAGNOSTICS */
412    
413     do kl = 2, Nr
414    
415     #ifdef ALLOW_AUTODIFF_TAMC
416     kkppkey = (ikppkey-1)*Nr + kl
417     #endif
418    
419     c compute bfsfc = sw fraction at hbf * zgrid
420    
421     do i = 1, imt
422     worka(i) = zgrid(kl)
423     end do
424     CADJ store worka = comlev1_kpp_k, key = kkppkey, kind=isbyte
425     call SWFRAC(
426     I imt, hbf,
427     U worka,
428     I myTime, myIter, myThid )
429     CADJ store worka = comlev1_kpp_k, key = kkppkey, kind=isbyte
430    
431     do i = 1, imt
432    
433     c use caseA as temporary array
434    
435     casea(i) = -zgrid(kl)
436    
437     c compute bfsfc= Bo + radiative contribution down to hbf * hbl
438    
439     bfsfc(i) = bo(i) + bosol(i)*(1. - worka(i))
440    
441     end do
442     #ifdef ALLOW_SALT_PLUME
443     c compute bfsfc = plume fraction at hbf * zgrid
444     IF ( useSALT_PLUME ) THEN
445     do i = 1, imt
446     worka(i) = zgrid(kl)
447     enddo
448 atn 1.3 #ifndef SALT_PLUME_VOLUME
449 atn 1.4 catn: in original way: accumulate all fractions of boplume above zgrid(kl)
450 atn 1.1 call SALT_PLUME_FRAC(
451     I imt, hbf,SPDepth,
452     U worka,
453     I myTime, myIter, myThid)
454     do i = 1, imt
455 atn 1.4 bfsfc(i) = bfsfc(i) + boplume(i,1)*(worka(i))
456 atn 1.1 enddo
457 atn 1.3 #else /* def SALT_PLUME_VOLUME */
458     catn: in vol way: need to integrate down to hbl, so first locate
459     c k level associated with this hbl, then sum up all SPforc[T,S]
460     DO i = 1, imt
461 atn 1.4 DO k = 1, kl
462     IF (abs(worka(i)).GE.abs( zgrid(k)-hwide(k)/2.0 )) THEN
463 atn 1.3 bfsfc(i) = bfsfc(i) + boplume(i,k)
464     ENDIF
465     ENDDO
466     ENDDO
467     #endif /* ndef SALT_PLUME_VOLUME */
468 atn 1.1 ENDIF
469     #endif /* ALLOW_SALT_PLUME */
470    
471     #ifdef ALLOW_DIAGNOSTICS
472     do i = 1, imt
473     KPPBFSFC(i,kl) = bfsfc(i)
474     enddo
475     #endif /* ALLOW_DIAGNOSTICS */
476    
477     do i = 1, imt
478     stable(i) = p5 + sign(p5,bfsfc(i))
479     sigma(i) = stable(i) + (1. - stable(i)) * epsilon
480     enddo
481    
482     c-----------------------------------------------------------------------
483     c compute velocity scales at sigma, for hbl= caseA = -zgrid(kl)
484     c-----------------------------------------------------------------------
485    
486     call wscale (
487     I sigma, casea, ustar, bfsfc,
488     O wm, ws, myThid )
489     CADJ store ws = comlev1_kpp_k, key = kkppkey, kind=isbyte
490    
491     do i = 1, imt
492    
493     c-----------------------------------------------------------------------
494     c compute the turbulent shear contribution to Rib
495     c-----------------------------------------------------------------------
496    
497     bvsq = p5 *
498     1 ( dbloc(i,kl-1) / (zgrid(kl-1)-zgrid(kl ))+
499     2 dbloc(i,kl ) / (zgrid(kl )-zgrid(kl+1)))
500    
501     if (bvsq .eq. 0. _d 0) then
502     vtsq = 0. _d 0
503     else
504     vtsq = -zgrid(kl) * ws(i) * sqrt(abs(bvsq)) * Vtc
505     endif
506    
507     c compute bulk Richardson number at new level
508     c note: Ritop needs to be zero on land and ocean bottom
509     c points so that the following if statement gets triggered
510     c correctly; otherwise, hbl might get set to (big) negative
511     c values, that might exceed the limit for the "exp" function
512     c in "SWFRAC"
513    
514     c
515     c rg: assignment to double precision variable to avoid overflow
516     c ph: test for zero nominator
517     c
518    
519     tempVar1 = dvsq(i,kl) + vtsq
520     tempVar2 = max(tempVar1, phepsi)
521     Rib(i,kl) = Ritop(i,kl) / tempVar2
522     #ifdef ALLOW_DIAGNOSTICS
523     KPPRi(i,kl) = Rib(i,kl)
524     #endif
525    
526     end do
527     end do
528    
529     #ifdef ALLOW_DIAGNOSTICS
530     IF ( useDiagnostics ) THEN
531     CALL DIAGNOSTICS_FILL(KPPBFSFC,'KPPbfsfc',0,Nr,2,bi,bj,myThid)
532     CALL DIAGNOSTICS_FILL(KPPRi ,'KPPRi ',0,Nr,2,bi,bj,myThid)
533     ENDIF
534     #endif /* ALLOW_DIAGNOSTICS */
535    
536     cph(
537     cph without this store, there is a recomputation error for
538     cph rib in adbldepth (probably partial recomputation problem)
539     CADJ store Rib = comlev1_kpp
540     CADJ & , key=ikppkey, kind=isbyte,
541     CADJ & shape = (/ (sNx+2*OLx)*(sNy+2*OLy),Nr /)
542     cph)
543    
544     do kl = 2, Nr
545     do i = 1, imt
546     if (kbl(i).eq.kmtj(i) .and. Rib(i,kl).gt.Ricr) kbl(i) = kl
547     end do
548     end do
549    
550     CADJ store kbl = comlev1_kpp
551     CADJ & , key=ikppkey, kind=isbyte,
552     CADJ & shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
553    
554     do i = 1, imt
555     kl = kbl(i)
556     c linearly interpolate to find hbl where Rib = Ricr
557     if (kl.gt.1 .and. kl.lt.kmtj(i)) then
558     tempVar1 = (Rib(i,kl)-Rib(i,kl-1))
559     hbl(i) = -zgrid(kl-1) + (zgrid(kl-1)-zgrid(kl)) *
560     1 (Ricr - Rib(i,kl-1)) / tempVar1
561     endif
562     end do
563    
564     CADJ store hbl = comlev1_kpp
565     CADJ & , key=ikppkey, kind=isbyte,
566     CADJ & shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
567    
568     c-----------------------------------------------------------------------
569     c find stability and buoyancy forcing for boundary layer
570     c-----------------------------------------------------------------------
571    
572     do i = 1, imt
573     worka(i) = hbl(i)
574     end do
575     CADJ store worka = comlev1_kpp
576     CADJ & , key=ikppkey, kind=isbyte,
577     CADJ & shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
578     call SWFRAC(
579     I imt, minusone,
580     U worka,
581     I myTime, myIter, myThid )
582     CADJ store worka = comlev1_kpp
583     CADJ & , key=ikppkey, kind=isbyte,
584     CADJ & shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
585    
586     do i = 1, imt
587     bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i))
588     end do
589    
590     #ifdef ALLOW_SALT_PLUME
591 atn 1.3 IF ( useSALT_PLUME ) THEN
592 atn 1.1 #ifndef SALT_PLUME_VOLUME
593     do i = 1, imt
594     worka(i) = hbl(i)
595     enddo
596     call SALT_PLUME_FRAC(
597     I imt,minusone,SPDepth,
598     U worka,
599     I myTime, myIter, myThid )
600     do i = 1, imt
601     bfsfc(i) = bfsfc(i) + boplume(i) * (worka(i))
602     enddo
603 atn 1.3 #else /* def SALT_PLUME_VOLUME */
604     DO i = 1, imt
605     DO k = 1, Nr
606 atn 1.4 IF (hbl(i).GE.abs( zgrid(k)-hwide(k)/2.0 )) THEN
607 atn 1.3 bfsfc(i) = bfsfc(i) + boplume(i,k)
608     ENDIF
609     ENDDO
610     ENDDO
611     #endif /* ndef SALT_PLUME_VOLUME */
612 atn 1.1 ENDIF
613     #endif /* ALLOW_SALT_PLUME */
614     CADJ store bfsfc = comlev1_kpp
615     CADJ & , key=ikppkey, kind=isbyte,
616     CADJ & shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
617    
618     c-- ensure bfsfc is never 0
619     do i = 1, imt
620     stable(i) = p5 + sign( p5, bfsfc(i) )
621     bfsfc(i) = sign(eins,bfsfc(i))*max(phepsi,abs(bfsfc(i)))
622     end do
623    
624     cph(
625     cph added stable to store list to avoid extensive recomp.
626     CADJ store bfsfc, stable = comlev1_kpp
627     CADJ & , key=ikppkey, kind=isbyte,
628     CADJ & shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
629     cph)
630    
631     c-----------------------------------------------------------------------
632     c check hbl limits for hekman or hmonob
633     c ph: test for zero nominator
634     c-----------------------------------------------------------------------
635    
636     IF ( LimitHblStable ) THEN
637     do i = 1, imt
638     if (bfsfc(i) .gt. 0.0) then
639     hekman = cekman * ustar(i) / max(abs(Coriol(i)),phepsi)
640     hmonob = cmonob * ustar(i)*ustar(i)*ustar(i)
641     & / vonk / bfsfc(i)
642     hlimit = stable(i) * min(hekman,hmonob)
643     & + (stable(i)-1.) * zgrid(Nr)
644     hbl(i) = min(hbl(i),hlimit)
645     end if
646     end do
647     ENDIF
648    
649     CADJ store hbl = comlev1_kpp
650     CADJ & , key=ikppkey, kind=isbyte,
651     CADJ & shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
652    
653     do i = 1, imt
654     hbl(i) = max(hbl(i),minKPPhbl)
655     kbl(i) = kmtj(i)
656     end do
657    
658     CADJ store hbl = comlev1_kpp
659     CADJ & , key=ikppkey, kind=isbyte,
660     CADJ & shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
661    
662     c-----------------------------------------------------------------------
663     c find new kbl
664     c-----------------------------------------------------------------------
665    
666     do kl = 2, Nr
667     do i = 1, imt
668     if ( kbl(i).eq.kmtj(i) .and. (-zgrid(kl)).gt.hbl(i) ) then
669     kbl(i) = kl
670     endif
671     end do
672     end do
673    
674     c-----------------------------------------------------------------------
675     c find stability and buoyancy forcing for final hbl values
676     c-----------------------------------------------------------------------
677    
678     do i = 1, imt
679     worka(i) = hbl(i)
680     end do
681     CADJ store worka = comlev1_kpp
682     CADJ & , key=ikppkey, kind=isbyte,
683     CADJ & shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
684     call SWFRAC(
685     I imt, minusone,
686     U worka,
687     I myTime, myIter, myThid )
688     CADJ store worka = comlev1_kpp
689     CADJ & , key=ikppkey, kind=isbyte,
690     CADJ & shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
691    
692     do i = 1, imt
693     bfsfc(i) = bo(i) + bosol(i) * (1. - worka(i))
694     end do
695    
696     #ifdef ALLOW_SALT_PLUME
697 atn 1.3 IF ( useSALT_PLUME ) THEN
698 atn 1.1 #ifndef SALT_PLUME_VOLUME
699     do i = 1, imt
700     worka(i) = hbl(i)
701     enddo
702     call SALT_PLUME_FRAC(
703     I imt,minusone,SPDepth,
704     U worka,
705     I myTime, myIter, myThid )
706     do i = 1, imt
707     bfsfc(i) = bfsfc(i) + boplume(i) * (worka(i))
708     enddo
709 atn 1.3 #else /* def SALT_PLUME_VOLUME */
710     DO i = 1, imt
711     DO k = 1, Nr
712 atn 1.4 IF (hbl(i).GE.abs( zgrid(k)-hwide(k)/2.0 )) THEN
713 atn 1.3 bfsfc(i) = bfsfc(i) + boplume(i,k)
714     ENDIF
715     ENDDO
716     ENDDO
717     #endif /* ndef SALT_PLUME_VOLUME */
718 atn 1.1 ENDIF
719     #endif /* ALLOW_SALT_PLUME */
720     CADJ store bfsfc = comlev1_kpp
721     CADJ & , key=ikppkey, kind=isbyte,
722     CADJ & shape = (/ (sNx+2*OLx)*(sNy+2*OLy) /)
723    
724     c-- ensures bfsfc is never 0
725     do i = 1, imt
726     stable(i) = p5 + sign( p5, bfsfc(i) )
727     bfsfc(i) = sign(eins,bfsfc(i))*max(phepsi,abs(bfsfc(i)))
728     end do
729    
730     c-----------------------------------------------------------------------
731     c determine caseA and caseB
732     c-----------------------------------------------------------------------
733    
734     do i = 1, imt
735     casea(i) = p5 +
736     1 sign(p5, -zgrid(kbl(i)) - p5*hwide(kbl(i)) - hbl(i))
737     end do
738    
739     #endif /* ALLOW_KPP */
740    
741     return
742     end
743    
744     c*************************************************************************
745    
746     subroutine wscale (
747     I sigma, hbl, ustar, bfsfc,
748     O wm, ws,
749     I myThid )
750    
751     c compute turbulent velocity scales.
752     c use a 2D-lookup table for wm and ws as functions of ustar and
753     c zetahat (=vonk*sigma*hbl*bfsfc).
754     c
755     c note: the lookup table is only used for unstable conditions
756     c (zehat.le.0), in the stable domain wm (=ws) gets computed
757     c directly.
758     c
759     IMPLICIT NONE
760    
761     #include "SIZE.h"
762     #include "KPP_PARAMS.h"
763    
764     c input
765     c------
766     c sigma : normalized depth (d/hbl)
767     c hbl : boundary layer depth (m)
768     c ustar : surface friction velocity (m/s)
769     c bfsfc : total surface buoyancy flux (m^2/s^3)
770     c myThid : thread number for this instance of the routine
771     integer myThid
772     _RL sigma(imt)
773     _RL hbl (imt)
774     _RL ustar(imt)
775     _RL bfsfc(imt)
776    
777     c output
778     c--------
779     c wm, ws : turbulent velocity scales at sigma
780     _RL wm(imt), ws(imt)
781    
782     #ifdef ALLOW_KPP
783    
784     c local
785     c------
786     c zehat : = zeta * ustar**3
787     _RL zehat
788    
789     integer iz, izp1, ju, i, jup1
790     _RL udiff, zdiff, zfrac, ufrac, fzfrac, wam
791     _RL wbm, was, wbs, u3, tempVar
792    
793     c-----------------------------------------------------------------------
794     c use lookup table for zehat < zmax only; otherwise use
795     c stable formulae
796     c-----------------------------------------------------------------------
797    
798     do i = 1, imt
799     zehat = vonk*sigma(i)*hbl(i)*bfsfc(i)
800    
801     if (zehat .le. zmax) then
802    
803     zdiff = zehat - zmin
804     iz = int( zdiff / deltaz )
805     iz = min( iz, nni )
806     iz = max( iz, 0 )
807     izp1 = iz + 1
808    
809     udiff = ustar(i) - umin
810     ju = int( udiff / deltau )
811     ju = min( ju, nnj )
812     ju = max( ju, 0 )
813     jup1 = ju + 1
814    
815     zfrac = zdiff / deltaz - float(iz)
816     ufrac = udiff / deltau - float(ju)
817    
818     fzfrac= 1. - zfrac
819     wam = fzfrac * wmt(iz,jup1) + zfrac * wmt(izp1,jup1)
820     wbm = fzfrac * wmt(iz,ju ) + zfrac * wmt(izp1,ju )
821     wm(i) = (1.-ufrac) * wbm + ufrac * wam
822    
823     was = fzfrac * wst(iz,jup1) + zfrac * wst(izp1,jup1)
824     wbs = fzfrac * wst(iz,ju ) + zfrac * wst(izp1,ju )
825     ws(i) = (1.-ufrac) * wbs + ufrac * was
826    
827     else
828    
829     u3 = ustar(i) * ustar(i) * ustar(i)
830     tempVar = u3 + conc1 * zehat
831     wm(i) = vonk * ustar(i) * u3 / tempVar
832     ws(i) = wm(i)
833    
834     endif
835    
836     end do
837    
838     #endif /* ALLOW_KPP */
839    
840     return
841     end
842    
843     c*************************************************************************
844    
845     subroutine Ri_iwmix (
846     I kmtj, shsq, dbloc, dblocSm,
847     I diffusKzS, diffusKzT,
848     I ikppkey,
849     O diffus,
850     I myThid )
851    
852     c compute interior viscosity diffusivity coefficients due
853     c to shear instability (dependent on a local Richardson number),
854     c to background internal wave activity, and
855     c to static instability (local Richardson number < 0).
856    
857     IMPLICIT NONE
858    
859     #include "SIZE.h"
860     #include "EEPARAMS.h"
861     #include "PARAMS.h"
862     #include "KPP_PARAMS.h"
863     #ifdef ALLOW_AUTODIFF
864     # include "AUTODIFF_PARAMS.h"
865     # include "tamc.h"
866     #endif
867    
868     c input
869     c kmtj (imt) number of vertical layers on this row
870     c shsq (imt,Nr) (local velocity shear)^2 ((m/s)^2)
871     c dbloc (imt,Nr) local delta buoyancy (m/s^2)
872     c dblocSm(imt,Nr) horizontally smoothed dbloc (m/s^2)
873     c diffusKzS(imt,Nr)- background vertical diffusivity for scalars (m^2/s)
874     c diffusKzT(imt,Nr)- background vertical diffusivity for theta (m^2/s)
875     c myThid :: My Thread Id. number
876     integer kmtj (imt)
877     _RL shsq (imt,Nr)
878     _RL dbloc (imt,Nr)
879     _RL dblocSm (imt,Nr)
880     _RL diffusKzS(imt,Nr)
881     _RL diffusKzT(imt,Nr)
882     integer ikppkey
883     integer myThid
884    
885     c output
886     c diffus(imt,0:Nrp1,1) vertical viscosivity coefficient (m^2/s)
887     c diffus(imt,0:Nrp1,2) vertical scalar diffusivity (m^2/s)
888     c diffus(imt,0:Nrp1,3) vertical temperature diffusivity (m^2/s)
889     _RL diffus(imt,0:Nrp1,3)
890    
891     #ifdef ALLOW_KPP
892    
893     c local variables
894     c Rig local Richardson number
895     c fRi, fcon function of Rig
896     _RL Rig
897     _RL fRi, fcon
898     _RL ratio
899     integer i, ki, kp1
900     _RL c1, c0
901    
902     #ifdef ALLOW_KPP_VERTICALLY_SMOOTH
903     integer mr
904     CADJ INIT kpp_ri_tape_mr = common, 1
905     #endif
906    
907     c constants
908     c1 = 1. _d 0
909     c0 = 0. _d 0
910    
911     c-----------------------------------------------------------------------
912     c compute interior gradient Ri at all interfaces ki=1,Nr, (not surface)
913     c use diffus(*,*,1) as temporary storage of Ri to be smoothed
914     c use diffus(*,*,2) as temporary storage for Brunt-Vaisala squared
915     c set values at bottom and below to nearest value above bottom
916     #ifdef ALLOW_AUTODIFF_TAMC
917     C break data flow dependence on diffus
918     diffus(1,1,1) = 0.0
919    
920     do ki = 1, Nr
921     do i = 1, imt
922     diffus(i,ki,1) = 0.
923     diffus(i,ki,2) = 0.
924     diffus(i,ki,3) = 0.
925     enddo
926     enddo
927     #endif
928    
929     do ki = 1, Nr
930     do i = 1, imt
931     if (kmtj(i) .LE. 1 ) then
932     diffus(i,ki,1) = 0.
933     diffus(i,ki,2) = 0.
934     elseif (ki .GE. kmtj(i)) then
935     diffus(i,ki,1) = diffus(i,ki-1,1)
936     diffus(i,ki,2) = diffus(i,ki-1,2)
937     else
938     diffus(i,ki,1) = dblocSm(i,ki) * (zgrid(ki)-zgrid(ki+1))
939     & / max( Shsq(i,ki), phepsi )
940     diffus(i,ki,2) = dbloc(i,ki) / (zgrid(ki)-zgrid(ki+1))
941     endif
942     end do
943     end do
944     CADJ store diffus = comlev1_kpp, key=ikppkey, kind=isbyte
945    
946     c-----------------------------------------------------------------------
947     c vertically smooth Ri
948     #ifdef ALLOW_KPP_VERTICALLY_SMOOTH
949     do mr = 1, num_v_smooth_Ri
950    
951     CADJ store diffus(:,:,1) = kpp_ri_tape_mr
952     CADJ & , key=mr, shape=(/ (sNx+2*OLx)*(sNy+2*OLy),Nr+2 /)
953    
954     call z121 (
955     U diffus(1,0,1),
956     I myThid )
957     end do
958     #endif
959    
960     c-----------------------------------------------------------------------
961     c after smoothing loop
962    
963     do ki = 1, Nr
964     do i = 1, imt
965    
966     c evaluate f of Brunt-Vaisala squared for convection, store in fcon
967    
968     Rig = max ( diffus(i,ki,2) , BVSQcon )
969     ratio = min ( (BVSQcon - Rig) / BVSQcon, c1 )
970     fcon = c1 - ratio * ratio
971     fcon = fcon * fcon * fcon
972    
973     c evaluate f of smooth Ri for shear instability, store in fRi
974    
975     Rig = max ( diffus(i,ki,1), c0 )
976     ratio = min ( Rig / Riinfty , c1 )
977     fRi = c1 - ratio * ratio
978     fRi = fRi * fRi * fRi
979    
980     c ----------------------------------------------------------------------
981     c evaluate diffusivities and viscosity
982     c mixing due to internal waves, and shear and static instability
983    
984     kp1 = MIN(ki+1,Nr)
985     #ifdef EXCLUDE_KPP_SHEAR_MIX
986     diffus(i,ki,1) = viscArNr(1)
987     diffus(i,ki,2) = diffusKzS(i,kp1)
988     diffus(i,ki,3) = diffusKzT(i,kp1)
989     #else /* EXCLUDE_KPP_SHEAR_MIX */
990     # ifdef ALLOW_AUTODIFF
991     if ( inAdMode ) then
992     diffus(i,ki,1) = viscArNr(1)
993     diffus(i,ki,2) = diffusKzS(i,kp1)
994     diffus(i,ki,3) = diffusKzT(i,kp1)
995     else
996     # else /* ALLOW_AUTODIFF */
997     if ( .TRUE. ) then
998     # endif /* ALLOW_AUTODIFF */
999     diffus(i,ki,1) = viscArNr(1) + fcon*difmcon + fRi*difm0
1000     diffus(i,ki,2) = diffusKzS(i,kp1)+fcon*difscon+fRi*difs0
1001     diffus(i,ki,3) = diffusKzT(i,kp1)+fcon*diftcon+fRi*dift0
1002     endif
1003     #endif /* EXCLUDE_KPP_SHEAR_MIX */
1004     end do
1005     end do
1006    
1007     c ------------------------------------------------------------------------
1008     c set surface values to 0.0
1009    
1010     do i = 1, imt
1011     diffus(i,0,1) = c0
1012     diffus(i,0,2) = c0
1013     diffus(i,0,3) = c0
1014     end do
1015    
1016     #endif /* ALLOW_KPP */
1017    
1018     return
1019     end
1020    
1021     c*************************************************************************
1022    
1023     subroutine z121 (
1024     U v,
1025     I myThid )
1026    
1027     c Apply 121 smoothing in k to 2-d array V(i,k=1,Nr)
1028     c top (0) value is used as a dummy
1029     c bottom (Nrp1) value is set to input value from above.
1030    
1031     c Note that it is important to exclude from the smoothing any points
1032     c that are outside the range of the K(Ri) scheme, ie. >0.8, or <0.0.
1033     c Otherwise, there is interference with other physics, especially
1034     c penetrative convection.
1035    
1036     IMPLICIT NONE
1037     #include "SIZE.h"
1038     #include "KPP_PARAMS.h"
1039    
1040     c input/output
1041     c-------------
1042     c v : 2-D array to be smoothed in Nrp1 direction
1043     c myThid: thread number for this instance of the routine
1044     integer myThid
1045     _RL v(imt,0:Nrp1)
1046    
1047     #ifdef ALLOW_KPP
1048    
1049     c local
1050     _RL zwork, zflag
1051     _RL KRi_range(1:Nrp1)
1052     integer i, k, km1, kp1
1053    
1054     _RL p0 , p25 , p5 , p2
1055     parameter ( p0 = 0.0, p25 = 0.25, p5 = 0.5, p2 = 2.0 )
1056    
1057     KRi_range(Nrp1) = p0
1058    
1059     #ifdef ALLOW_AUTODIFF_TAMC
1060     C-- dummy assignment to end declaration part for TAMC
1061     i = 0
1062    
1063     C-- HPF directive to help TAMC
1064     CHPF$ INDEPENDENT
1065     CADJ INIT z121tape = common, Nr
1066     #endif /* ALLOW_AUTODIFF_TAMC */
1067    
1068     do i = 1, imt
1069    
1070     k = 1
1071     CADJ STORE v(i,k) = z121tape
1072     v(i,Nrp1) = v(i,Nr)
1073    
1074     do k = 1, Nr
1075     KRi_range(k) = p5 + SIGN(p5,v(i,k))
1076     KRi_range(k) = KRi_range(k) *
1077     & ( p5 + SIGN(p5,(Riinfty-v(i,k))) )
1078     end do
1079    
1080     zwork = KRi_range(1) * v(i,1)
1081     v(i,1) = p2 * v(i,1) +
1082     & KRi_range(1) * KRi_range(2) * v(i,2)
1083     zflag = p2 + KRi_range(1) * KRi_range(2)
1084     v(i,1) = v(i,1) / zflag
1085    
1086     do k = 2, Nr
1087     CADJ STORE v(i,k), zwork = z121tape
1088     km1 = k - 1
1089     kp1 = k + 1
1090     zflag = v(i,k)
1091     v(i,k) = p2 * v(i,k) +
1092     & KRi_range(k) * KRi_range(kp1) * v(i,kp1) +
1093     & KRi_range(k) * zwork
1094     zwork = KRi_range(k) * zflag
1095     zflag = p2 + KRi_range(k)*(KRi_range(kp1)+KRi_range(km1))
1096     v(i,k) = v(i,k) / zflag
1097     end do
1098    
1099     end do
1100    
1101     #endif /* ALLOW_KPP */
1102    
1103     return
1104     end
1105    
1106     c*************************************************************************
1107    
1108     subroutine smooth_horiz (
1109     I k, bi, bj,
1110     U fld,
1111     I myThid )
1112    
1113     c Apply horizontal smoothing to global _RL 2-D array
1114    
1115     IMPLICIT NONE
1116     #include "SIZE.h"
1117     #include "GRID.h"
1118     #include "KPP_PARAMS.h"
1119    
1120     c input
1121     c bi, bj : array indices
1122     c k : vertical index used for masking
1123     c myThid : thread number for this instance of the routine
1124     INTEGER myThid
1125     integer k, bi, bj
1126    
1127     c input/output
1128     c fld : 2-D array to be smoothed
1129     _RL fld( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
1130    
1131     #ifdef ALLOW_KPP
1132    
1133     c local
1134     integer i, j, im1, ip1, jm1, jp1
1135     _RL tempVar
1136     _RL fld_tmp( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
1137    
1138     integer imin , imax , jmin , jmax
1139     parameter(imin=2-OLx, imax=sNx+OLx-1, jmin=2-OLy, jmax=sNy+OLy-1)
1140    
1141     _RL p0 , p5 , p25 , p125 , p0625
1142     parameter( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 )
1143    
1144     DO j = jmin, jmax
1145     jm1 = j-1
1146     jp1 = j+1
1147     DO i = imin, imax
1148     im1 = i-1
1149     ip1 = i+1
1150     tempVar =
1151     & p25 * maskC(i ,j ,k,bi,bj) +
1152     & p125 * ( maskC(im1,j ,k,bi,bj) +
1153     & maskC(ip1,j ,k,bi,bj) +
1154     & maskC(i ,jm1,k,bi,bj) +
1155     & maskC(i ,jp1,k,bi,bj) ) +
1156     & p0625 * ( maskC(im1,jm1,k,bi,bj) +
1157     & maskC(im1,jp1,k,bi,bj) +
1158     & maskC(ip1,jm1,k,bi,bj) +
1159     & maskC(ip1,jp1,k,bi,bj) )
1160     IF ( tempVar .GE. p25 ) THEN
1161     fld_tmp(i,j) = (
1162     & p25 * fld(i ,j )*maskC(i ,j ,k,bi,bj) +
1163     & p125 *(fld(im1,j )*maskC(im1,j ,k,bi,bj) +
1164     & fld(ip1,j )*maskC(ip1,j ,k,bi,bj) +
1165     & fld(i ,jm1)*maskC(i ,jm1,k,bi,bj) +
1166     & fld(i ,jp1)*maskC(i ,jp1,k,bi,bj))+
1167     & p0625*(fld(im1,jm1)*maskC(im1,jm1,k,bi,bj) +
1168     & fld(im1,jp1)*maskC(im1,jp1,k,bi,bj) +
1169     & fld(ip1,jm1)*maskC(ip1,jm1,k,bi,bj) +
1170     & fld(ip1,jp1)*maskC(ip1,jp1,k,bi,bj)))
1171     & / tempVar
1172     ELSE
1173     fld_tmp(i,j) = fld(i,j)
1174     ENDIF
1175     ENDDO
1176     ENDDO
1177    
1178     c transfer smoothed field to output array
1179     DO j = jmin, jmax
1180     DO i = imin, imax
1181     fld(i,j) = fld_tmp(i,j)
1182     ENDDO
1183     ENDDO
1184    
1185     #endif /* ALLOW_KPP */
1186    
1187     return
1188     end
1189    
1190     c*************************************************************************
1191    
1192     subroutine blmix (
1193     I ustar, bfsfc, hbl, stable, casea, diffus, kbl
1194     O , dkm1, blmc, ghat, sigma, ikppkey
1195     I , myThid )
1196    
1197     c mixing coefficients within boundary layer depend on surface
1198     c forcing and the magnitude and gradient of interior mixing below
1199     c the boundary layer ("matching").
1200     c
1201     c caution: if mixing bottoms out at hbl = -zgrid(Nr) then
1202     c fictitious layer at Nrp1 is needed with small but finite width
1203     c hwide(Nrp1) (eg. epsln = 1.e-20).
1204     c
1205     IMPLICIT NONE
1206    
1207     #include "SIZE.h"
1208     #include "KPP_PARAMS.h"
1209     #ifdef ALLOW_AUTODIFF
1210     # include "tamc.h"
1211     #endif
1212    
1213     c input
1214     c ustar (imt) surface friction velocity (m/s)
1215     c bfsfc (imt) surface buoyancy forcing (m^2/s^3)
1216     c hbl (imt) boundary layer depth (m)
1217     c stable(imt) = 1 in stable forcing
1218     c casea (imt) = 1 in case A
1219     c diffus(imt,0:Nrp1,mdiff) vertical diffusivities (m^2/s)
1220     c kbl (imt) -1 of first grid level below hbl
1221     c myThid thread number for this instance of the routine
1222     integer myThid
1223     _RL ustar (imt)
1224     _RL bfsfc (imt)
1225     _RL hbl (imt)
1226     _RL stable(imt)
1227     _RL casea (imt)
1228     _RL diffus(imt,0:Nrp1,mdiff)
1229     integer kbl(imt)
1230    
1231     c output
1232     c dkm1 (imt,mdiff) boundary layer difs at kbl-1 level
1233     c blmc (imt,Nr,mdiff) boundary layer mixing coefficients (m^2/s)
1234     c ghat (imt,Nr) nonlocal scalar transport
1235     c sigma(imt) normalized depth (d / hbl)
1236     _RL dkm1 (imt,mdiff)
1237     _RL blmc (imt,Nr,mdiff)
1238     _RL ghat (imt,Nr)
1239     _RL sigma(imt)
1240     integer ikppkey
1241    
1242     #ifdef ALLOW_KPP
1243    
1244     c local
1245     c gat1*(imt) shape function at sigma = 1
1246     c dat1*(imt) derivative of shape function at sigma = 1
1247     c ws(imt), wm(imt) turbulent velocity scales (m/s)
1248     _RL gat1m(imt), gat1s(imt), gat1t(imt)
1249     _RL dat1m(imt), dat1s(imt), dat1t(imt)
1250     _RL ws(imt), wm(imt)
1251     integer i, kn, ki
1252     _RL R, dvdzup, dvdzdn, viscp
1253     _RL difsp, diftp, visch, difsh, difth
1254     _RL f1, sig, a1, a2, a3, delhat
1255     _RL Gm, Gs, Gt
1256     _RL tempVar
1257    
1258     _RL p0 , eins
1259     parameter (p0=0.0, eins=1.0)
1260     #ifdef ALLOW_AUTODIFF_TAMC
1261     integer kkppkey
1262     #endif
1263    
1264     c-----------------------------------------------------------------------
1265     c compute velocity scales at hbl
1266     c-----------------------------------------------------------------------
1267    
1268     do i = 1, imt
1269     sigma(i) = stable(i) * 1.0 + (1. - stable(i)) * epsilon
1270     end do
1271    
1272     CADJ STORE sigma = comlev1_kpp, key=ikppkey, kind=isbyte
1273     call wscale (
1274     I sigma, hbl, ustar, bfsfc,
1275     O wm, ws, myThid )
1276     CADJ STORE wm = comlev1_kpp, key=ikppkey, kind=isbyte
1277     CADJ STORE ws = comlev1_kpp, key=ikppkey, kind=isbyte
1278    
1279     do i = 1, imt
1280     wm(i) = sign(eins,wm(i))*max(phepsi,abs(wm(i)))
1281     ws(i) = sign(eins,ws(i))*max(phepsi,abs(ws(i)))
1282     end do
1283     CADJ STORE wm = comlev1_kpp, key=ikppkey, kind=isbyte
1284     CADJ STORE ws = comlev1_kpp, key=ikppkey, kind=isbyte
1285    
1286     do i = 1, imt
1287    
1288     kn = int(caseA(i)+phepsi) *(kbl(i) -1) +
1289     $ (1 - int(caseA(i)+phepsi)) * kbl(i)
1290    
1291     c-----------------------------------------------------------------------
1292     c find the interior viscosities and derivatives at hbl(i)
1293     c-----------------------------------------------------------------------
1294    
1295     delhat = 0.5*hwide(kn) - zgrid(kn) - hbl(i)
1296     R = 1.0 - delhat / hwide(kn)
1297     dvdzup = (diffus(i,kn-1,1) - diffus(i,kn ,1)) / hwide(kn)
1298     dvdzdn = (diffus(i,kn ,1) - diffus(i,kn+1,1)) / hwide(kn+1)
1299     viscp = 0.5 * ( (1.-R) * (dvdzup + abs(dvdzup)) +
1300     1 R * (dvdzdn + abs(dvdzdn)) )
1301    
1302     dvdzup = (diffus(i,kn-1,2) - diffus(i,kn ,2)) / hwide(kn)
1303     dvdzdn = (diffus(i,kn ,2) - diffus(i,kn+1,2)) / hwide(kn+1)
1304     difsp = 0.5 * ( (1.-R) * (dvdzup + abs(dvdzup)) +
1305     1 R * (dvdzdn + abs(dvdzdn)) )
1306    
1307     dvdzup = (diffus(i,kn-1,3) - diffus(i,kn ,3)) / hwide(kn)
1308     dvdzdn = (diffus(i,kn ,3) - diffus(i,kn+1,3)) / hwide(kn+1)
1309     diftp = 0.5 * ( (1.-R) * (dvdzup + abs(dvdzup)) +
1310     1 R * (dvdzdn + abs(dvdzdn)) )
1311    
1312     visch = diffus(i,kn,1) + viscp * delhat
1313     difsh = diffus(i,kn,2) + difsp * delhat
1314     difth = diffus(i,kn,3) + diftp * delhat
1315    
1316     f1 = stable(i) * conc1 * bfsfc(i) /
1317     & max(ustar(i)**4,phepsi)
1318     gat1m(i) = visch / hbl(i) / wm(i)
1319     dat1m(i) = -viscp / wm(i) + f1 * visch
1320    
1321     gat1s(i) = difsh / hbl(i) / ws(i)
1322     dat1s(i) = -difsp / ws(i) + f1 * difsh
1323    
1324     gat1t(i) = difth / hbl(i) / ws(i)
1325     dat1t(i) = -diftp / ws(i) + f1 * difth
1326    
1327     end do
1328     #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1329     CADJ STORE gat1m = comlev1_kpp, key=ikppkey, kind=isbyte
1330     CADJ STORE gat1s = comlev1_kpp, key=ikppkey, kind=isbyte
1331     CADJ STORE gat1t = comlev1_kpp, key=ikppkey, kind=isbyte
1332     CADJ STORE dat1m = comlev1_kpp, key=ikppkey, kind=isbyte
1333     CADJ STORE dat1s = comlev1_kpp, key=ikppkey, kind=isbyte
1334     CADJ STORE dat1t = comlev1_kpp, key=ikppkey, kind=isbyte
1335     #endif
1336     do i = 1, imt
1337     dat1m(i) = min(dat1m(i),p0)
1338     dat1s(i) = min(dat1s(i),p0)
1339     dat1t(i) = min(dat1t(i),p0)
1340     end do
1341     #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1342     CADJ STORE dat1m = comlev1_kpp, key=ikppkey, kind=isbyte
1343     CADJ STORE dat1s = comlev1_kpp, key=ikppkey, kind=isbyte
1344     CADJ STORE dat1t = comlev1_kpp, key=ikppkey, kind=isbyte
1345     #endif
1346    
1347     do ki = 1, Nr
1348    
1349     #ifdef ALLOW_AUTODIFF_TAMC
1350     kkppkey = (ikppkey-1)*Nr + ki
1351     #endif
1352    
1353     c-----------------------------------------------------------------------
1354     c compute turbulent velocity scales on the interfaces
1355     c-----------------------------------------------------------------------
1356    
1357     do i = 1, imt
1358     sig = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i)
1359     sigma(i) = stable(i)*sig + (1.-stable(i))*min(sig,epsilon)
1360     end do
1361     #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1362     CADJ STORE wm = comlev1_kpp_k, key = kkppkey
1363     CADJ STORE ws = comlev1_kpp_k, key = kkppkey
1364     #endif
1365     CADJ STORE sigma = comlev1_kpp_k, key = kkppkey
1366     call wscale (
1367     I sigma, hbl, ustar, bfsfc,
1368     O wm, ws, myThid )
1369     CADJ STORE wm = comlev1_kpp_k, key = kkppkey
1370     CADJ STORE ws = comlev1_kpp_k, key = kkppkey
1371    
1372     c-----------------------------------------------------------------------
1373     c compute the dimensionless shape functions at the interfaces
1374     c-----------------------------------------------------------------------
1375    
1376     do i = 1, imt
1377     sig = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl(i)
1378     a1 = sig - 2.
1379     a2 = 3. - 2. * sig
1380     a3 = sig - 1.
1381    
1382     Gm = a1 + a2 * gat1m(i) + a3 * dat1m(i)
1383     Gs = a1 + a2 * gat1s(i) + a3 * dat1s(i)
1384     Gt = a1 + a2 * gat1t(i) + a3 * dat1t(i)
1385    
1386     c-----------------------------------------------------------------------
1387     c compute boundary layer diffusivities at the interfaces
1388     c-----------------------------------------------------------------------
1389    
1390     blmc(i,ki,1) = hbl(i) * wm(i) * sig * (1. + sig * Gm)
1391     blmc(i,ki,2) = hbl(i) * ws(i) * sig * (1. + sig * Gs)
1392     blmc(i,ki,3) = hbl(i) * ws(i) * sig * (1. + sig * Gt)
1393    
1394     c-----------------------------------------------------------------------
1395     c nonlocal transport term = ghat * <ws>o
1396     c-----------------------------------------------------------------------
1397    
1398     tempVar = ws(i) * hbl(i)
1399     ghat(i,ki) = (1.-stable(i)) * cg / max(phepsi,tempVar)
1400    
1401     end do
1402     end do
1403    
1404     c-----------------------------------------------------------------------
1405     c find diffusivities at kbl-1 grid level
1406     c-----------------------------------------------------------------------
1407    
1408     do i = 1, imt
1409     sig = -zgrid(kbl(i)-1) / hbl(i)
1410     sigma(i) = stable(i) * sig
1411     & + (1. - stable(i)) * min(sig,epsilon)
1412     end do
1413    
1414     #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1415     CADJ STORE wm = comlev1_kpp, key=ikppkey, kind=isbyte
1416     CADJ STORE ws = comlev1_kpp, key=ikppkey, kind=isbyte
1417     #endif
1418     CADJ STORE sigma = comlev1_kpp, key=ikppkey, kind=isbyte
1419     call wscale (
1420     I sigma, hbl, ustar, bfsfc,
1421     O wm, ws, myThid )
1422     CADJ STORE wm = comlev1_kpp, key=ikppkey, kind=isbyte
1423     CADJ STORE ws = comlev1_kpp, key=ikppkey, kind=isbyte
1424    
1425     do i = 1, imt
1426     sig = -zgrid(kbl(i)-1) / hbl(i)
1427     a1 = sig - 2.
1428     a2 = 3. - 2. * sig
1429     a3 = sig - 1.
1430     Gm = a1 + a2 * gat1m(i) + a3 * dat1m(i)
1431     Gs = a1 + a2 * gat1s(i) + a3 * dat1s(i)
1432     Gt = a1 + a2 * gat1t(i) + a3 * dat1t(i)
1433     dkm1(i,1) = hbl(i) * wm(i) * sig * (1. + sig * Gm)
1434     dkm1(i,2) = hbl(i) * ws(i) * sig * (1. + sig * Gs)
1435     dkm1(i,3) = hbl(i) * ws(i) * sig * (1. + sig * Gt)
1436     end do
1437    
1438     #endif /* ALLOW_KPP */
1439    
1440     return
1441     end
1442    
1443     c*************************************************************************
1444    
1445     subroutine enhance (
1446     I dkm1, hbl, kbl, diffus, casea
1447     U , ghat
1448     O , blmc
1449     & , myThid )
1450    
1451     c enhance the diffusivity at the kbl-.5 interface
1452    
1453     IMPLICIT NONE
1454    
1455     #include "SIZE.h"
1456     #include "KPP_PARAMS.h"
1457    
1458     c input
1459     c dkm1(imt,mdiff) bl diffusivity at kbl-1 grid level
1460     c hbl(imt) boundary layer depth (m)
1461     c kbl(imt) grid above hbl
1462     c diffus(imt,0:Nrp1,mdiff) vertical diffusivities (m^2/s)
1463     c casea(imt) = 1 in caseA, = 0 in case B
1464     c myThid thread number for this instance of the routine
1465     integer myThid
1466     _RL dkm1 (imt,mdiff)
1467     _RL hbl (imt)
1468     integer kbl (imt)
1469     _RL diffus(imt,0:Nrp1,mdiff)
1470     _RL casea (imt)
1471    
1472     c input/output
1473     c nonlocal transport, modified ghat at kbl(i)-1 interface (s/m**2)
1474     _RL ghat (imt,Nr)
1475    
1476     c output
1477     c enhanced bound. layer mixing coeff.
1478     _RL blmc (imt,Nr,mdiff)
1479    
1480     #ifdef ALLOW_KPP
1481    
1482     c local
1483     c fraction hbl lies beteen zgrid neighbors
1484     _RL delta
1485     integer ki, i, md
1486     _RL dkmp5, dstar
1487    
1488     do i = 1, imt
1489     ki = kbl(i)-1
1490     if ((ki .ge. 1) .and. (ki .lt. Nr)) then
1491     delta = (hbl(i) + zgrid(ki)) / (zgrid(ki) - zgrid(ki+1))
1492     do md = 1, mdiff
1493     dkmp5 = casea(i) * diffus(i,ki,md) +
1494     1 (1.- casea(i)) * blmc (i,ki,md)
1495     dstar = (1.- delta)**2 * dkm1(i,md)
1496     & + delta**2 * dkmp5
1497     blmc(i,ki,md) = (1.- delta)*diffus(i,ki,md)
1498     & + delta*dstar
1499     end do
1500     ghat(i,ki) = (1.- casea(i)) * ghat(i,ki)
1501     endif
1502     end do
1503    
1504     #endif /* ALLOW_KPP */
1505    
1506     return
1507     end
1508    
1509     c*************************************************************************
1510    
1511     SUBROUTINE STATEKPP (
1512     O RHO1, DBLOC, DBSFC, TTALPHA, SSBETA,
1513     I ikppkey, bi, bj, myThid )
1514     c
1515     c-----------------------------------------------------------------------
1516     c "statekpp" computes all necessary input arrays
1517     c for the kpp mixing scheme
1518     c
1519     c input:
1520     c bi, bj = array indices on which to apply calculations
1521     c
1522     c output:
1523     c rho1 = potential density of surface layer (kg/m^3)
1524     c dbloc = local buoyancy gradient at Nr interfaces
1525     c g/rho{k+1,k+1} * [ drho{k,k+1}-drho{k+1,k+1} ] (m/s^2)
1526     c dbsfc = buoyancy difference with respect to the surface
1527     c g * [ drho{1,k}/rho{1,k} - drho{k,k}/rho{k,k} ] (m/s^2)
1528     c ttalpha= thermal expansion coefficient without 1/rho factor
1529     c d(rho) / d(potential temperature) (kg/m^3/C)
1530     c ssbeta = salt expansion coefficient without 1/rho factor
1531     c d(rho) / d(salinity) (kg/m^3/PSU)
1532     c
1533     c see also subroutines find_rho.F find_alpha.F find_beta.F
1534     c
1535     c written by: jan morzel, feb. 10, 1995 (converted from "sigma" version)
1536     c modified by: d. menemenlis, june 1998 : for use with MIT GCM UV
1537     c
1538    
1539     c-----------------------------------------------------------------------
1540    
1541     IMPLICIT NONE
1542    
1543     #include "SIZE.h"
1544     #include "EEPARAMS.h"
1545     #include "PARAMS.h"
1546     #include "KPP_PARAMS.h"
1547     #include "DYNVARS.h"
1548     #include "GRID.h"
1549     #ifdef ALLOW_AUTODIFF
1550     # include "tamc.h"
1551     #endif
1552    
1553     c-------------- Routine arguments -----------------------------------------
1554     INTEGER bi, bj, myThid
1555     _RL RHO1 ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
1556     _RL DBLOC ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr )
1557     _RL DBSFC ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr )
1558     _RL TTALPHA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1559     _RL SSBETA ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1560    
1561     #ifdef ALLOW_KPP
1562    
1563     c--------------------------------------------------------------------------
1564     c
1565     c local arrays:
1566     c
1567     c rhok - density of t(k ) & s(k ) at depth k
1568     c rhokm1 - density of t(k-1) & s(k-1) at depth k
1569     c rho1k - density of t(1 ) & s(1 ) at depth k
1570     c work1,2,3 - work arrays for holding horizontal slabs
1571    
1572     _RL RHOK (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1573     _RL RHOKM1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1574     _RL RHO1K (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1575     _RL WORK1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1576     _RL WORK2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1577     _RL WORK3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1578    
1579     INTEGER I, J, K
1580     INTEGER ikppkey, kkppkey
1581    
1582     c calculate density, alpha, beta in surface layer, and set dbsfc to zero
1583    
1584     kkppkey = (ikppkey-1)*Nr + 1
1585    
1586     #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1587     CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k,
1588     CADJ & key=kkppkey, kind=isbyte
1589     CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k,
1590     CADJ & key=kkppkey, kind=isbyte
1591     #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1592     CALL FIND_RHO_2D(
1593     I 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1,
1594     I theta(1-OLx,1-OLy,1,bi,bj), salt(1-OLx,1-OLy,1,bi,bj),
1595     O WORK1,
1596     I 1, bi, bj, myThid )
1597     #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1598     CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k,
1599     CADJ & key=kkppkey, kind=isbyte
1600     CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k,
1601     CADJ & key=kkppkey, kind=isbyte
1602     #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1603    
1604     call FIND_ALPHA(
1605     I bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1606     O WORK2, myThid )
1607    
1608     call FIND_BETA(
1609     I bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
1610     O WORK3, myThid )
1611    
1612     DO J = 1-OLy, sNy+OLy
1613     DO I = 1-OLx, sNx+OLx
1614     RHO1(I,J) = WORK1(I,J) + rhoConst
1615     TTALPHA(I,J,1) = WORK2(I,J)
1616     SSBETA(I,J,1) = WORK3(I,J)
1617     DBSFC(I,J,1) = 0.
1618     END DO
1619     END DO
1620    
1621     c calculate alpha, beta, and gradients in interior layers
1622    
1623     CHPF$ INDEPENDENT, NEW (RHOK,RHOKM1,RHO1K,WORK1,WORK2)
1624     DO K = 2, Nr
1625    
1626     kkppkey = (ikppkey-1)*Nr + k
1627    
1628     #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1629     CADJ STORE theta(:,:,k,bi,bj) = comlev1_kpp_k,
1630     CADJ & key=kkppkey, kind=isbyte
1631     CADJ STORE salt (:,:,k,bi,bj) = comlev1_kpp_k,
1632     CADJ & key=kkppkey, kind=isbyte
1633     #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1634     CALL FIND_RHO_2D(
1635     I 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, k,
1636     I theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
1637     O RHOK,
1638     I k, bi, bj, myThid )
1639    
1640     #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1641     CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_kpp_k,
1642     CADJ & key=kkppkey, kind=isbyte
1643     CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_kpp_k,
1644     CADJ & key=kkppkey, kind=isbyte
1645     #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1646     CALL FIND_RHO_2D(
1647     I 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, k,
1648     I theta(1-OLx,1-OLy,k-1,bi,bj),salt(1-OLx,1-OLy,k-1,bi,bj),
1649     O RHOKM1,
1650     I k-1, bi, bj, myThid )
1651    
1652     #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1653     CADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k,
1654     CADJ & key=kkppkey, kind=isbyte
1655     CADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k,
1656     CADJ & key=kkppkey, kind=isbyte
1657     #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1658     CALL FIND_RHO_2D(
1659     I 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, k,
1660     I theta(1-OLx,1-OLy,1,bi,bj), salt(1-OLx,1-OLy,1,bi,bj),
1661     O RHO1K,
1662     I 1, bi, bj, myThid )
1663    
1664     #ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1665     CADJ STORE rhok (:,:) = comlev1_kpp_k,
1666     CADJ & key=kkppkey, kind=isbyte
1667     CADJ STORE rhokm1(:,:) = comlev1_kpp_k,
1668     CADJ & key=kkppkey, kind=isbyte
1669     CADJ STORE rho1k (:,:) = comlev1_kpp_k,
1670     CADJ & key=kkppkey, kind=isbyte
1671     #endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1672    
1673     call FIND_ALPHA(
1674     I bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1675     O WORK1, myThid )
1676    
1677     call FIND_BETA(
1678     I bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, K,
1679     O WORK2, myThid )
1680    
1681     DO J = 1-OLy, sNy+OLy
1682     DO I = 1-OLx, sNx+OLx
1683     TTALPHA(I,J,K) = WORK1 (I,J)
1684     SSBETA(I,J,K) = WORK2 (I,J)
1685     DBLOC(I,J,K-1) = gravity * (RHOK(I,J) - RHOKM1(I,J)) /
1686     & (RHOK(I,J) + rhoConst)
1687     DBSFC(I,J,K) = gravity * (RHOK(I,J) - RHO1K (I,J)) /
1688     & (RHOK(I,J) + rhoConst)
1689     END DO
1690     END DO
1691    
1692     END DO
1693    
1694     c compute arrays for K = Nrp1
1695     DO J = 1-OLy, sNy+OLy
1696     DO I = 1-OLx, sNx+OLx
1697     TTALPHA(I,J,Nrp1) = TTALPHA(I,J,Nr)
1698     SSBETA(I,J,Nrp1) = SSBETA(I,J,Nr)
1699     DBLOC(I,J,Nr) = 0.
1700     END DO
1701     END DO
1702    
1703     #ifdef ALLOW_DIAGNOSTICS
1704     IF ( useDiagnostics ) THEN
1705     CALL DIAGNOSTICS_FILL(DBSFC ,'KPPdbsfc',0,Nr,2,bi,bj,myThid)
1706     CALL DIAGNOSTICS_FILL(DBLOC ,'KPPdbloc',0,Nr,2,bi,bj,myThid)
1707     ENDIF
1708     #endif /* ALLOW_DIAGNOSTICS */
1709    
1710     #endif /* ALLOW_KPP */
1711    
1712     RETURN
1713     END
1714    
1715     c*************************************************************************
1716    
1717     SUBROUTINE KPP_DOUBLEDIFF (
1718     I TTALPHA, SSBETA,
1719     U kappaRT,
1720     U kappaRS,
1721     I ikppkey, imin, imax, jmin, jmax, bi, bj, myThid )
1722     c
1723     c-----------------------------------------------------------------------
1724     c "KPP_DOUBLEDIFF" adds the double diffusive contributions
1725     C as Rrho-dependent parameterizations to kappaRT and kappaRS
1726     c
1727     c input:
1728     c bi, bj = array indices on which to apply calculations
1729     c imin, imax, jmin, jmax = array boundaries
1730     c ikppkey = key for TAMC/TAF automatic differentiation
1731     c myThid = thread id
1732     c
1733     c ttalpha= thermal expansion coefficient without 1/rho factor
1734     c d(rho) / d(potential temperature) (kg/m^3/C)
1735     c ssbeta = salt expansion coefficient without 1/rho factor
1736     c d(rho) / d(salinity) (kg/m^3/PSU)
1737     c output: updated
1738     c kappaRT/S :: background diffusivities for temperature and salinity
1739     c
1740     c written by: martin losch, sept. 15, 2009
1741     c
1742    
1743     c-----------------------------------------------------------------------
1744    
1745     IMPLICIT NONE
1746    
1747     #include "SIZE.h"
1748     #include "EEPARAMS.h"
1749     #include "PARAMS.h"
1750     #include "KPP_PARAMS.h"
1751     #include "DYNVARS.h"
1752     #include "GRID.h"
1753     #ifdef ALLOW_AUTODIFF
1754     # include "tamc.h"
1755     #endif
1756    
1757     c-------------- Routine arguments -----------------------------------------
1758     INTEGER ikppkey, imin, imax, jmin, jmax, bi, bj, myThid
1759    
1760     _RL TTALPHA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1761     _RL SSBETA ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nrp1 )
1762     _RL KappaRT( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr )
1763     _RL KappaRS( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, Nr )
1764    
1765     #ifdef ALLOW_KPP
1766    
1767     C--------------------------------------------------------------------------
1768     C
1769     C local variables
1770     C I,J,K :: loop indices
1771     C kkppkey :: key for TAMC/TAF automatic differentiation
1772     C
1773     INTEGER I, J, K
1774     INTEGER kkppkey
1775     C alphaDT :: d\rho/d\theta * d\theta
1776     C betaDS :: d\rho/dsalt * dsalt
1777     C Rrho :: "density ratio" R_{\rho} = \alpha dT/dz / \beta dS/dz
1778     C nuddt/s :: double diffusive diffusivities
1779     C numol :: molecular diffusivity
1780     C rFac :: abbreviation for 1/(R_{\rho0}-1)
1781    
1782     _RL alphaDT ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
1783     _RL betaDS ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
1784     _RL nuddt ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
1785     _RL nudds ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
1786     _RL Rrho
1787     _RL numol, rFac, nutmp
1788     INTEGER Km1
1789    
1790     C set some constants here
1791     numol = 1.5 _d -06
1792     rFac = 1. _d 0 / (Rrho0 - 1. _d 0 )
1793     C
1794     kkppkey = (ikppkey-1)*Nr + 1
1795    
1796     CML#ifdef KPP_AUTODIFF_EXCESSIVE_STORE
1797     CMLCADJ STORE theta(:,:,1,bi,bj) = comlev1_kpp_k,
1798     CMLCADJ & key=kkppkey, kind=isbyte
1799     CMLCADJ STORE salt (:,:,1,bi,bj) = comlev1_kpp_k,
1800     CMLCADJ & key=kkppkey, kind=isbyte
1801     CML#endif /* KPP_AUTODIFF_EXCESSIVE_STORE */
1802    
1803     DO K = 1, Nr
1804     Km1 = MAX(K-1,1)
1805     DO J = 1-OLy, sNy+OLy
1806     DO I = 1-OLx, sNx+OLx
1807     alphaDT(I,J) = ( theta(I,J,Km1,bi,bj)-theta(I,J,K,bi,bj) )
1808     & * 0.5 _d 0 * ABS( TTALPHA(I,J,Km1) + TTALPHA(I,J,K) )
1809     betaDS(I,J) = ( salt(I,J,Km1,bi,bj)-salt(I,J,K,bi,bj) )
1810     & * 0.5 _d 0 * ( SSBETA(I,J,Km1) + SSBETA(I,J,K) )
1811     nuddt(I,J) = 0. _d 0
1812     nudds(I,J) = 0. _d 0
1813     ENDDO
1814     ENDDO
1815     IF ( K .GT. 1 ) THEN
1816     DO J = jMin, jMax
1817     DO I = iMin, iMax
1818     Rrho = 0. _d 0
1819     C Now we have many different cases
1820     C a. alphaDT > 0 and betaDS > 0 => salt fingering
1821     C (salinity destabilizes)
1822     IF ( alphaDT(I,J) .GT. betaDS(I,J)
1823     & .AND. betaDS(I,J) .GT. 0. _d 0 ) THEN
1824     Rrho = MIN( alphaDT(I,J)/betaDS(I,J), Rrho0 )
1825     C Large et al. 1994, eq. 31a
1826     C nudds(I,J) = dsfmax * ( 1. _d 0 - (Rrho - 1. _d 0) * rFac )**3
1827     nutmp = ( 1. _d 0 - (Rrho - 1. _d 0) * rFac )
1828     nudds(I,J) = dsfmax * nutmp * nutmp * nutmp
1829     C Large et al. 1994, eq. 31c
1830     nuddt(I,J) = 0.7 _d 0 * nudds(I,J)
1831     ELSEIF ( alphaDT(I,J) .LT. 0. _d 0
1832     & .AND. betaDS(I,J) .LT. 0. _d 0
1833     & .AND.alphaDT(I,J) .GT. betaDS(I,J) ) THEN
1834     C b. alphaDT < 0 and betaDS < 0 => semi-convection, diffusive convection
1835     C (temperature destabilizes)
1836     C for Rrho >= 1 the water column is statically unstable and we never
1837     C reach this point
1838     Rrho = alphaDT(I,J)/betaDS(I,J)
1839     C Large et al. 1994, eq. 32
1840     nuddt(I,J) = numol * 0.909 _d 0
1841     & * exp ( 4.6 _d 0 * exp (
1842     & - 5.4 _d 0 * ( 1. _d 0/Rrho - 1. _d 0 ) ) )
1843     CMLC or
1844     CMLC Large et al. 1994, eq. 33
1845     CML nuddt(I,J) = numol * 8.7 _d 0 * Rrho**1.1
1846     C Large et al. 1994, eqs. 34
1847     nudds(I,J) = nuddt(I,J) * MAX( 0.15 _d 0 * Rrho,
1848     & 1.85 _d 0 * Rrho - 0.85 _d 0 )
1849     ELSE
1850     C Do nothing, because in this case the water colume is unstable
1851     C => double diffusive processes are negligible and mixing due
1852     C to shear instability will dominate
1853     ENDIF
1854     ENDDO
1855     ENDDO
1856     C ENDIF ( K .GT. 1 )
1857     ENDIF
1858     C
1859     DO J = 1-OLy, sNy+OLy
1860     DO I = 1-OLx, sNx+OLx
1861     kappaRT(I,J,K) = kappaRT(I,J,K) + nuddt(I,J)
1862     kappaRS(I,J,K) = kappaRS(I,J,K) + nudds(I,J)
1863     ENDDO
1864     ENDDO
1865     #ifdef ALLOW_DIAGNOSTICS
1866     IF ( useDiagnostics ) THEN
1867     CALL DIAGNOSTICS_FILL(nuddt,'KPPnuddt',k,1,2,bi,bj,myThid)
1868     CALL DIAGNOSTICS_FILL(nudds,'KPPnudds',k,1,2,bi,bj,myThid)
1869     ENDIF
1870     #endif /* ALLOW_DIAGNOSTICS */
1871     C end of K-loop
1872     ENDDO
1873     #endif /* ALLOW_KPP */
1874    
1875     RETURN
1876     END

  ViewVC Help
Powered by ViewVC 1.1.22