/[MITgcm]/MITgcm/model/src/calc_3d_diffusivity.F
ViewVC logotype

Annotation of /MITgcm/model/src/calc_3d_diffusivity.F

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


Revision 1.20 - (hide annotations) (download)
Thu Dec 7 19:41:22 2017 UTC (6 years, 6 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint66o, checkpoint66n, checkpoint66m, HEAD
Changes since 1.19: +2 -2 lines
fix typo in pCelMix code (hFacC instead of hFacW)

1 jmc 1.20 C $Header: /u/gcmpack/MITgcm/model/src/calc_3d_diffusivity.F,v 1.19 2017/11/02 18:00:52 jmc Exp $
2 jmc 1.1 C $Name: $
3    
4     #include "PACKAGES_CONFIG.h"
5     #include "CPP_OPTIONS.h"
6    
7     CBOP
8     C !ROUTINE: CALC_3D_DIFFUSIVITY
9 jmc 1.11 C !INTERFACE:
10     SUBROUTINE CALC_3D_DIFFUSIVITY(
11 jmc 1.19 I bi, bj, iMin,iMax,jMin,jMax,
12 jmc 1.1 I trIdentity, trUseGMRedi, trUseKPP,
13     O KappaRTr,
14 jmc 1.19 I myThid )
15 jmc 1.1
16     C !DESCRIPTION: \bv
17     C *==========================================================*
18     C | SUBROUTINE CALC_3D_DIFFUSIVITY
19     C | o Calculate net (3D) vertical diffusivity for 1 tracer
20     C *==========================================================*
21     C | Combines spatially varying diffusion coefficients from
22     C | KPP and/or GM and/or convective stability test.
23     C *==========================================================*
24     C \ev
25    
26     C !USES:
27     IMPLICIT NONE
28     C == GLobal variables ==
29     #include "SIZE.h"
30     #include "EEPARAMS.h"
31     #include "PARAMS.h"
32     #include "DYNVARS.h"
33     #include "GRID.h"
34 jmc 1.3 #ifdef ALLOW_GENERIC_ADVDIFF
35 jmc 1.1 #include "GAD.h"
36 jmc 1.3 #endif
37 jmc 1.1 #ifdef ALLOW_PTRACERS
38     #include "PTRACERS_SIZE.h"
39 jmc 1.13 #include "PTRACERS_PARAMS.h"
40 jmc 1.1 #endif
41 jahn 1.15 #ifdef ALLOW_LONGSTEP
42     #include "LONGSTEP.h"
43     #endif
44 jmc 1.1
45     C !INPUT/OUTPUT PARAMETERS:
46     C == Routine arguments ==
47     C bi, bj :: tile indices
48     C iMin,iMax :: Range of points for which calculation is performed.
49     C jMin,jMax :: Range of points for which calculation is performed.
50     C trIdentity :: tracer identifier
51     C trUseGMRedi:: this tracer use GM-Redi
52     C trUseKPP :: this tracer use KPP
53     C myThid :: Instance number for this innvocation of CALC_3D_DIFFUSIVITY
54     C KappaRTr :: Net diffusivity for this tracer (trIdentity)
55     INTEGER bi,bj,iMin,iMax,jMin,jMax
56     INTEGER trIdentity
57     LOGICAL trUseGMRedi, trUseKPP
58 jmc 1.18 _RL KappaRTr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
59 jmc 1.1 INTEGER myThid
60    
61 jmc 1.3 #ifdef ALLOW_GENERIC_ADVDIFF
62 jmc 1.1 C !LOCAL VARIABLES:
63     C == Local variables ==
64     C i, j, k :: Loop counters
65     C iTr :: passive tracer index
66     C msgBuf :: message buffer
67     INTEGER i,j,k
68 jmc 1.11 _RL KbryanLewis79
69     #ifdef ALLOW_BL79_LAT_VARY
70     _RL KbryanLewisEQ
71     #endif
72 jmc 1.1 CHARACTER*(MAX_LEN_MBUF) msgBuf
73 jmc 1.4 #ifdef ALLOW_PTRACERS
74     INTEGER iTr
75     #endif
76 jmc 1.19 #ifndef EXCLUDE_PCELL_MIX_CODE
77     INTEGER km, mixSurf, mixBott
78     _RL pC_kFac
79     _RL tmpFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80     #endif
81 jmc 1.1 CEOP
82    
83 dimitri 1.10 IF ( .NOT. trUseKPP ) THEN
84 jmc 1.1 DO k = 1,Nr
85 jmc 1.8 KbryanLewis79=diffKrBL79surf+(diffKrBL79deep-diffKrBL79surf)
86 dimitri 1.10 & *(atan(-(rF(k)-diffKrBL79Ho)/diffKrBL79scl)/PI+0.5 _d 0)
87 dimitri 1.9 #ifdef ALLOW_BL79_LAT_VARY
88     KbryanLewisEQ=diffKrBLEQsurf+(diffKrBLEQdeep-diffKrBLEQsurf)
89 dimitri 1.10 & *(atan(-(rF(k)-diffKrBLEQHo)/diffKrBLEQscl)/PI+0.5 _d 0)
90 dimitri 1.9 #endif
91 jmc 1.18 DO j = 1-OLy,sNy+OLy
92     DO i = 1-OLx,sNx+OLx
93 jahn 1.15 #ifdef ALLOW_LONGSTEP
94     IF ( trIdentity .GE. GAD_TR1) THEN
95     KappaRTr(i,j,k) =
96     & LS_IVDConvCount(i,j,k,bi,bj)*ivdc_kappa
97     & + KbryanLewis79
98     #ifdef ALLOW_BL79_LAT_VARY
99     & + (KbryanLewisEQ-KbryanLewis79)*BL79LatArray(i,j,bi,bj)
100     #endif
101     ELSE
102     #else
103     IF ( .TRUE. ) THEN
104     #endif /* ALLOW_LONGSTEP */
105     KappaRTr(i,j,k) =
106 jmc 1.1 & IVDConvCount(i,j,k,bi,bj)*ivdc_kappa
107 dimitri 1.10 & + KbryanLewis79
108     #ifdef ALLOW_BL79_LAT_VARY
109     & + (KbryanLewisEQ-KbryanLewis79)*BL79LatArray(i,j,bi,bj)
110     #endif
111 jahn 1.15 ENDIF
112 dimitri 1.10 ENDDO
113     ENDDO
114     ENDDO
115     IF ( trIdentity.EQ.GAD_TEMPERATURE ) THEN
116     DO k = 1,Nr
117 jmc 1.18 DO j = 1-OLy,sNy+OLy
118     DO i = 1-OLx,sNx+OLx
119 jmc 1.11 KappaRTr(i,j,k) = KappaRTr(i,j,k)
120 jmc 1.17 #ifdef ALLOW_3D_DIFFKR
121 dimitri 1.10 & + diffKr(i,j,k,bi,bj)
122 jmc 1.1 #else
123 dimitri 1.10 & + diffKrNrT(k)
124 dimitri 1.9 #endif
125 dimitri 1.10 ENDDO
126 jmc 1.1 ENDDO
127     ENDDO
128 dimitri 1.10 ELSEIF ( trIdentity.EQ.GAD_SALINITY) THEN
129     DO k = 1,Nr
130 jmc 1.18 DO j = 1-OLy, sNy+OLy
131     DO i = 1-OLx, sNx+OLx
132 dimitri 1.10 KappaRTr(i,j,k) = KappaRTr(i,j,k)
133 jmc 1.17 #ifdef ALLOW_3D_DIFFKR
134 dimitri 1.10 & + diffKr(i,j,k,bi,bj)
135 jmc 1.1 #else
136 dimitri 1.10 & + diffKrNrS(k)
137 dimitri 1.9 #endif
138 dimitri 1.10 ENDDO
139 jmc 1.1 ENDDO
140     ENDDO
141     #ifdef ALLOW_PTRACERS
142 jahn 1.15 ELSEIF ( trIdentity.GE.GAD_TR1) THEN
143 jmc 1.1
144 dimitri 1.10 iTr = trIdentity - GAD_TR1 + 1
145     DO k = 1,Nr
146 jmc 1.18 DO j = 1-OLy, sNy+OLy
147     DO i = 1-OLx, sNx+OLx
148 dimitri 1.10 KappaRTr(i,j,k) = KappaRTr(i,j,k)
149 jmc 1.17 #ifdef ALLOW_3D_DIFFKR
150 dimitri 1.10 & + diffKr(i,j,k,bi,bj)
151 jmc 1.1 #else
152 dimitri 1.10 & + PTRACERS_diffKrNr(k,iTr)
153 dimitri 1.9 #endif
154 dimitri 1.10 ENDDO
155 jmc 1.1 ENDDO
156     ENDDO
157     #endif /* ALLOW_PTRACERS */
158 dimitri 1.10 ELSE
159 jmc 1.1 WRITE(msgBuf,'(A,I4)')
160 dimitri 1.10 & ' CALC_3D_DIFFUSIVITY: Invalid tracer Id: ',trIdentity
161 jmc 1.1 CALL PRINT_ERROR(msgBuf, myThid)
162     STOP 'ABNORMAL END: S/R CALC_3D_DIFFUSIVITY'
163 dimitri 1.10 ENDIF
164 jmc 1.1 ENDIF
165    
166     C-- Add physical pacakge contributions:
167    
168     #ifdef ALLOW_KPP
169     IF (trUseKPP) THEN
170 dimitri 1.12 C-- Set vertical diffusivity contribution from KPP
171 jmc 1.1 IF (trIdentity.EQ.GAD_TEMPERATURE) THEN
172     CALL KPP_CALC_DIFF_T(
173 jmc 1.2 I bi,bj,iMin,iMax,jMin,jMax,0,Nr,
174 dimitri 1.12 O KappaRTr,
175 jmc 1.1 I myThid)
176 jahn 1.15 ELSEIF (trIdentity.EQ.GAD_SALINITY) THEN
177 jmc 1.1 CALL KPP_CALC_DIFF_S(
178 jmc 1.2 I bi,bj,iMin,iMax,jMin,jMax,0,Nr,
179 dimitri 1.12 O KappaRTr,
180 jmc 1.1 I myThid)
181 jmc 1.17 #ifdef ALLOW_PTRACERS
182 jahn 1.15 ELSEIF ( trIdentity.GE.GAD_TR1) THEN
183 jmc 1.17 iTr = trIdentity - GAD_TR1 + 1
184 jahn 1.15 CALL KPP_CALC_DIFF_Ptr(
185     I bi,bj,iMin,iMax,jMin,jMax,0,Nr,
186     O KappaRTr,
187 jmc 1.17 I iTr, myThid )
188     #endif /* ALLOW_PTRACERS */
189     ELSE
190     WRITE(msgBuf,'(A,I4)')
191     & ' CALC_3D_DIFFUSIVITY: Invalid tracer Id: ',trIdentity
192     CALL PRINT_ERROR( msgBuf, myThid )
193     STOP 'ABNORMAL END: S/R CALC_3D_DIFFUSIVITY'
194 jmc 1.1 ENDIF
195 dimitri 1.10 ENDIF
196     #endif /* ALLOW_KPP */
197    
198     #ifdef ALLOW_GMREDI
199 jmc 1.11 IF (trUseGMRedi) THEN
200 dimitri 1.10 CALL GMREDI_CALC_DIFF(
201     I bi,bj,iMin,iMax,jMin,jMax,0,Nr,
202     U KappaRTr,
203 dfer 1.14 I trIdentity,myThid)
204 jmc 1.1 ENDIF
205     #endif
206    
207     #ifdef ALLOW_PP81
208     IF (usePP81) THEN
209     CALL PP81_CALC_DIFF(
210 jmc 1.2 I bi,bj,iMin,iMax,jMin,jMax,0,Nr,
211 jmc 1.1 U KappaRTr,
212     I myThid)
213     ENDIF
214     #endif
215    
216 jmc 1.18 #ifdef ALLOW_KL10
217     IF (useKL10) THEN
218     CALL KL10_CALC_DIFF(
219     I bi,bj,iMin,iMax,jMin,jMax,0,Nr,
220     U KappaRTr,
221     I myThid)
222     ENDIF
223     #endif
224    
225 jmc 1.1 #ifdef ALLOW_MY82
226     IF (useMY82) THEN
227     CALL MY82_CALC_DIFF(
228 jmc 1.2 I bi,bj,iMin,iMax,jMin,jMax,0,Nr,
229 jmc 1.1 U KappaRTr,
230     I myThid)
231     ENDIF
232     #endif
233 jmc 1.11
234 jmc 1.1 #ifdef ALLOW_GGL90
235     IF (useGGL90) THEN
236     CALL GGL90_CALC_DIFF(
237 jmc 1.2 I bi,bj,iMin,iMax,jMin,jMax,0,Nr,
238 jmc 1.1 O KappaRTr,
239     I myThid)
240     ENDIF
241     #endif
242 jmc 1.11
243 jmc 1.19 #ifndef EXCLUDE_PCELL_MIX_CODE
244     IF ( interDiffKr_pCell ) THEN
245     C-- This is a hack: alter vertical diffusivity (instead of changing many S/R)
246     C in order to account for missing hFac in diffusion term
247     DO k = 2,Nr
248     km = k - 1
249     C- account for true distance (including hFac) in vertical gradient
250     DO j = 2-OLy, sNy+OLy
251     DO i = 2-OLx, sNx+OLx
252     IF ( k.GT.kSurfC(i,j,bi,bj) .AND.
253     & k.LE.kLowC(i,j,bi,bj) ) THEN
254     KappaRTr(i,j,k) = KappaRTr(i,j,k)
255     & *twoRL/(hFacC(i,j,km,bi,bj)+hFacC(i,j,k,bi,bj))
256     ENDIF
257     ENDDO
258     ENDDO
259     ENDDO
260     ENDIF
261    
262     IF ( pCellMix_select.GT.0 ) THEN
263     C-- This is a hack: alter vertical diffusivity (instead of changing many S/R)
264     C in order to to increase mixing for too thin surface/bottom partial cell
265     mixSurf = pCellMix_select/10
266     mixBott = MOD(pCellMix_select,10)
267     DO k = 2,Nr
268     km = k - 1
269     pC_kFac = 1.
270     IF ( pCellMix_delR.LT.drF(k) )
271     & pC_kFac = pCellMix_delR*recip_drF(k)
272    
273     C- Increase KappaRTr above bottom level:
274     IF ( mixBott.GE.1 ) THEN
275     DO j = 2-OLy, sNy+OLy
276     DO i = 2-OLx, sNx+OLx
277     tmpFac(i,j) = 0. _d 0
278     IF ( k.EQ.kLowC(i,j,bi,bj) .AND.
279     & k.GT.kSurfC(i,j,bi,bj) ) THEN
280     tmpFac(i,j) = pC_kFac*_recip_hFacC(i,j,k,bi,bj)
281     ENDIF
282     ENDDO
283     ENDDO
284     ENDIF
285     IF ( mixBott.EQ.2 ) THEN
286     DO j = 2-OLy, sNy+OLy
287     DO i = 2-OLx, sNx+OLx
288     tmpFac(i,j) = tmpFac(i,j)*tmpFac(i,j)
289     ENDDO
290     ENDDO
291     ELSEIF ( mixBott.EQ.3 ) THEN
292     DO j = 2-OLy, sNy+OLy
293     DO i = 2-OLx, sNx+OLx
294     tmpFac(i,j) = tmpFac(i,j)*tmpFac(i,j)*tmpFac(i,j)
295     ENDDO
296     ENDDO
297     ELSEIF ( mixBott.EQ.4 ) THEN
298     DO j = 2-OLy, sNy+OLy
299     DO i = 2-OLx, sNx+OLx
300     tmpFac(i,j) = tmpFac(i,j)*tmpFac(i,j)
301     & * tmpFac(i,j)*tmpFac(i,j)
302     ENDDO
303     ENDDO
304     ENDIF
305     IF ( mixBott.GE.1 ) THEN
306     C- increase mixing above bottom (by ~(1/hFac)^mixBott) if too thin p-cell
307     DO j = 2-OLy, sNy+OLy
308     DO i = 2-OLx, sNx+OLx
309     tmpFac(i,j) = MIN( tmpFac(i,j), pCellMix_maxFac )
310     KappaRTr(i,j,k) = MAX( KappaRTr(i,j,k),
311     & pCellMix_diffKr(k)*tmpFac(i,j) )
312     ENDDO
313     ENDDO
314     ENDIF
315    
316     pC_kFac = 1.
317     IF ( pCellMix_delR.LT.drF(km) )
318     & pC_kFac = pCellMix_delR*recip_drF(km)
319    
320     C- Increase KappaRTr below surface level:
321     IF ( mixSurf.GE.1 ) THEN
322     DO j = 2-OLy, sNy+OLy
323     DO i = 2-OLx, sNx+OLx
324     tmpFac(i,j) = 0. _d 0
325     IF ( km.EQ.kSurfC(i,j,bi,bj) .AND.
326     & km.LT.kLowC(i,j,bi,bj) ) THEN
327 jmc 1.20 tmpFac(i,j) = pC_kFac*_recip_hFacC(i,j,km,bi,bj)
328 jmc 1.19 ENDIF
329     ENDDO
330     ENDDO
331     ENDIF
332     IF ( mixSurf.EQ.2 ) THEN
333     DO j = 2-OLy, sNy+OLy
334     DO i = 2-OLx, sNx+OLx
335     tmpFac(i,j) = tmpFac(i,j)*tmpFac(i,j)
336     ENDDO
337     ENDDO
338     ELSEIF ( mixSurf.EQ.3 ) THEN
339     DO j = 2-OLy, sNy+OLy
340     DO i = 2-OLx, sNx+OLx
341     tmpFac(i,j) = tmpFac(i,j)*tmpFac(i,j)*tmpFac(i,j)
342     ENDDO
343     ENDDO
344     ELSEIF ( mixSurf.EQ.4 ) THEN
345     DO j = 2-OLy, sNy+OLy
346     DO i = 2-OLx, sNx+OLx
347     tmpFac(i,j) = tmpFac(i,j)*tmpFac(i,j)
348     & * tmpFac(i,j)*tmpFac(i,j)
349     ENDDO
350     ENDDO
351     ENDIF
352     IF ( mixSurf.GE.1 ) THEN
353     C- increase mixing below surface (by ~(1/hFac)^mixSurf) if too thin p-cell
354     DO j = 2-OLy, sNy+OLy
355     DO i = 2-OLx, sNx+OLx
356     tmpFac(i,j) = MIN( tmpFac(i,j), pCellMix_maxFac )
357     KappaRTr(i,j,k) = MAX( KappaRTr(i,j,k),
358     & pCellMix_diffKr(k)*tmpFac(i,j) )
359     ENDDO
360     ENDDO
361     ENDIF
362    
363     C-- end of k loop
364     ENDDO
365     ENDIF
366     #endif /* ndef EXCLUDE_PCELL_MIX_CODE */
367    
368 jmc 1.11 C- Apply mask to vertical diffusivity
369 jmc 1.16 C jmc: do not have the impression that masking is needed
370     C but could be removed later if it is the case.
371 jmc 1.18 c DO j = 1-OLy, sNy+OLy
372     c DO i = 1-OLx, sNx+OLx
373 jmc 1.1 c KappaRTr(i,j,k) = maskUp(i,j)*KappaRTr(i,j,k)
374     c ENDDO
375     c ENDDO
376    
377 jmc 1.3 #endif /* ALLOW_GENERIC_ADVDIFF */
378    
379 jmc 1.1 RETURN
380     END

  ViewVC Help
Powered by ViewVC 1.1.22