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

Contents of /MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/salt_integrate.F

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


Revision 1.3 - (show annotations) (download)
Fri May 2 05:35:22 2014 UTC (11 years, 2 months ago) by atn
Branch: MAIN
CVS Tags: HEAD
Changes since 1.2: +2 -2 lines
bug fix. move #include back inside autodiff bracket.

1 C $Header: /u/gcmpack/MITgcm_contrib/atnguyen/code_21Dec2012_saltplume/salt_integrate.F,v 1.2 2014/04/29 06:49:40 atn Exp $
2 C $Name: $
3
4 #include "PACKAGES_CONFIG.h"
5 #include "CPP_OPTIONS.h"
6 #ifdef ALLOW_AUTODIFF
7 # include "AUTODIFF_OPTIONS.h"
8 #endif
9 #ifdef ALLOW_GENERIC_ADVDIFF
10 # include "GAD_OPTIONS.h"
11 #endif
12 #ifdef ALLOW_SALT_PLUME
13 # include "SALT_PLUME_OPTIONS.h"
14 #endif
15
16 CBOP
17 C !ROUTINE: SALT_INTEGRATE
18 C !INTERFACE:
19 SUBROUTINE SALT_INTEGRATE(
20 I bi, bj, recip_hFac,
21 I uFld, vFld, wFld,
22 U KappaRk,
23 I myTime, myIter, myThid )
24 C !DESCRIPTION: \bv
25 C *==========================================================*
26 C | SUBROUTINE SALT_INTEGRATE
27 C | o Calculate tendency for salt
28 C | and integrates forward in time.
29 C *==========================================================*
30 C | A procedure called EXTERNAL_FORCING_S is called from
31 C | here. These procedures can be used to add per problem
32 C | E-P flux source terms.
33 C | Note: Although it is slightly counter-intuitive the
34 C | EXTERNAL_FORCING routine is not the place to put
35 C | file I/O. Instead files that are required to
36 C | calculate the external source terms are generally
37 C | read during the model main loop. This makes the
38 C | logistics of multi-processing simpler and also
39 C | makes the adjoint generation simpler. It also
40 C | allows for I/O to overlap computation where that
41 C | is supported by hardware.
42 C | Aside from the problem specific term the code here
43 C | forms the tendency terms due to advection and mixing
44 C | The baseline implementation here uses a centered
45 C | difference form for the advection term and a tensorial
46 C | divergence of a flux form for the diffusive term. The
47 C | diffusive term is formulated so that isopycnal mixing
48 C | and GM-style subgrid-scale terms can be incorporated by
49 C | simply setting the diffusion tensor terms appropriately.
50 C *==========================================================*
51 C \ev
52
53 C !USES:
54 IMPLICIT NONE
55 C == GLobal variables ==
56 #include "SIZE.h"
57 #include "EEPARAMS.h"
58 #include "PARAMS.h"
59 #include "GRID.h"
60 #include "DYNVARS.h"
61 #include "RESTART.h"
62 #ifdef ALLOW_GENERIC_ADVDIFF
63 # include "GAD.h"
64 # include "GAD_SOM_VARS.h"
65 #endif
66 #ifdef ALLOW_AUTODIFF
67 # include "tamc.h"
68 # include "tamc_keys.h"
69 #endif
70
71 C !INPUT/OUTPUT PARAMETERS:
72 C == Routine arguments ==
73 C bi, bj, :: tile indices
74 C recip_hFac :: reciprocal of cell open-depth factor (@ next iter)
75 C uFld,vFld :: Local copy of horizontal velocity field
76 C wFld :: Local copy of vertical velocity field
77 C KappaRk :: Vertical diffusion for Salinity
78 C myTime :: current time
79 C myIter :: current iteration number
80 C myThid :: my Thread Id. number
81 INTEGER bi, bj
82 _RS recip_hFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
83 _RL uFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
84 _RL vFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
85 _RL wFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
86 _RL KappaRk (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
87 _RL myTime
88 INTEGER myIter
89 INTEGER myThid
90 CEOP
91
92 #ifdef ALLOW_GENERIC_ADVDIFF
93 C !LOCAL VARIABLES:
94 C iMin, iMax :: 1rst index loop range
95 C jMin, jMax :: 2nd index loop range
96 C k :: vertical index
97 C kM1 :: =k-1 for k>1, =1 for k=1
98 C kUp :: index into 2 1/2D array, toggles between 1|2
99 C kDown :: index into 2 1/2D array, toggles between 2|1
100 C xA :: Tracer cell face area normal to X
101 C yA :: Tracer cell face area normal to X
102 C maskUp :: Land/water mask for Wvel points (interface k)
103 C uTrans :: Zonal volume transport through cell face
104 C vTrans :: Meridional volume transport through cell face
105 C rTrans :: Vertical volume transport at interface k
106 C rTransKp :: Vertical volume transport at inteface k+1
107 C fZon :: Flux of salt (S) in the zonal direction
108 C fMer :: Flux of salt (S) in the meridional direction
109 C fVer :: Flux of salt (S) in the vertical direction
110 C at the upper(U) and lower(D) faces of a cell.
111 INTEGER iMin, iMax, jMin, jMax
112 INTEGER i, j, k
113 INTEGER kUp, kDown, kM1
114 _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115 _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116 _RS maskUp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
117 _RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118 _RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
119 _RL rTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
120 _RL rTransKp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121 _RL fZon (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122 _RL fMer (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
123 _RL fVer (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
124 _RL gs_AB (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
125 LOGICAL calcAdvection
126 INTEGER iterNb
127 #ifdef ALLOW_ADAMSBASHFORTH_3
128 INTEGER m1, m2
129 #endif
130 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
131
132 C- Loop ranges for daughter routines
133 iMin = 1-OLx+2
134 iMax = sNx+OLx-1
135 jMin = 1-OLy+2
136 jMax = sNy+OLy-1
137
138 iterNb = myIter
139 IF (staggerTimeStep) iterNb = myIter - 1
140
141 #ifdef ALLOW_AUTODIFF_TAMC
142 act1 = bi - myBxLo(myThid)
143 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
144 act2 = bj - myByLo(myThid)
145 max2 = myByHi(myThid) - myByLo(myThid) + 1
146 act3 = myThid - 1
147 max3 = nTx*nTy
148 act4 = ikey_dynamics - 1
149 itdkey = (act1 + 1) + act2*max1
150 & + act3*max1*max2
151 & + act4*max1*max2*max3
152 #endif /* ALLOW_AUTODIFF_TAMC */
153
154 C- Tracer tendency needs to be set to zero (moved here from gad_calc_rhs):
155 DO k=1,Nr
156 DO j=1-OLy,sNy+OLy
157 DO i=1-OLx,sNx+OLx
158 gS(i,j,k,bi,bj) = 0. _d 0
159 ENDDO
160 ENDDO
161 ENDDO
162 DO j=1-OLy,sNy+OLy
163 DO i=1-OLx,sNx+OLx
164 fVer(i,j,1) = 0. _d 0
165 fVer(i,j,2) = 0. _d 0
166 ENDDO
167 ENDDO
168 #ifdef ALLOW_AUTODIFF
169 DO k=1,Nr
170 DO j=1-OLy,sNy+OLy
171 DO i=1-OLx,sNx+OLx
172 kappaRk(i,j,k) = 0. _d 0
173 ENDDO
174 ENDDO
175 ENDDO
176 CADJ STORE salt(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
177 CADJ STORE wFld(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
178 #endif /* ALLOW_AUTODIFF */
179
180 #ifdef INCLUDE_CALC_DIFFUSIVITY_CALL
181 CALL CALC_3D_DIFFUSIVITY(
182 I bi, bj, iMin, iMax, jMin, jMax,
183 I GAD_SALINITY, useGMredi, useKPP,
184 O kappaRk,
185 I myThid )
186 #endif /* INCLUDE_CALC_DIFFUSIVITY_CALL */
187
188 #ifndef DISABLE_MULTIDIM_ADVECTION
189 C-- Some advection schemes are better calculated using a multi-dimensional
190 C method in the absence of any other terms and, if used, is done here.
191 C
192 C The CPP flag DISABLE_MULTIDIM_ADVECTION is currently unset in GAD_OPTIONS.h
193 C The default is to use multi-dimensinal advection for non-linear advection
194 C schemes. However, for the sake of efficiency of the adjoint it is necessary
195 C to be able to exclude this scheme to avoid excessive storage and
196 C recomputation. It *is* differentiable, if you need it.
197 C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
198 C disable this section of code.
199 #ifdef GAD_ALLOW_TS_SOM_ADV
200 # ifdef ALLOW_AUTODIFF_TAMC
201 CADJ STORE som_S = comlev1_bibj, key=itdkey, byte=isbyte
202 # endif
203 IF ( saltSOM_Advection ) THEN
204 # ifdef ALLOW_DEBUG
205 IF (debugMode) CALL DEBUG_CALL('GAD_SOM_ADVECT',myThid)
206 # endif
207 CALL GAD_SOM_ADVECT(
208 I saltImplVertAdv, saltAdvScheme, saltVertAdvScheme,
209 I GAD_SALINITY, dTtracerLev,
210 I uFld, vFld, wFld, salt,
211 U som_S,
212 O gS,
213 I bi, bj, myTime, myIter, myThid )
214 ELSEIF (saltMultiDimAdvec) THEN
215 #else /* GAD_ALLOW_TS_SOM_ADV */
216 IF (saltMultiDimAdvec) THEN
217 #endif /* GAD_ALLOW_TS_SOM_ADV */
218 # ifdef ALLOW_DEBUG
219 IF (debugMode) CALL DEBUG_CALL('GAD_ADVECTION',myThid)
220 # endif
221 CALL GAD_ADVECTION(
222 I saltImplVertAdv, saltAdvScheme, saltVertAdvScheme,
223 I GAD_SALINITY, dTtracerLev,
224 I uFld, vFld, wFld, salt,
225 O gS,
226 I bi, bj, myTime, myIter, myThid )
227 ENDIF
228 #endif /* DISABLE_MULTIDIM_ADVECTION */
229
230 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
231
232 C- Start vertical index (k) loop (Nr:1)
233 calcAdvection = saltAdvection .AND. .NOT.saltMultiDimAdvec
234 DO k=Nr,1,-1
235 #ifdef ALLOW_AUTODIFF_TAMC
236 kkey = (itdkey-1)*Nr + k
237 #endif
238 kM1 = MAX(1,k-1)
239 kUp = 1+MOD(k+1,2)
240 kDown= 1+MOD(k,2)
241
242 #ifdef ALLOW_AUTODIFF_TAMC
243 CADJ STORE fVer(:,:,:) = comlev1_bibj_k, key=kkey,
244 CADJ & byte=isbyte, kind = isbyte
245 CADJ STORE gS(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey,
246 CADJ & byte=isbyte, kind = isbyte
247 # ifdef ALLOW_ADAMSBASHFORTH_3
248 CADJ STORE gsNm(:,:,k,bi,bj,1) = comlev1_bibj_k, key=kkey,
249 CADJ & byte=isbyte, kind = isbyte
250 CADJ STORE gsNm(:,:,k,bi,bj,2) = comlev1_bibj_k, key=kkey,
251 CADJ & kind = isbyte
252 # else
253 CADJ STORE gsNm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey,
254 CADJ & byte=isbyte, kind = isbyte
255 # endif
256 #endif /* ALLOW_AUTODIFF_TAMC */
257 CALL CALC_ADV_FLOW(
258 I uFld, vFld, wFld,
259 U rTrans,
260 O uTrans, vTrans, rTransKp,
261 O maskUp, xA, yA,
262 I k, bi, bj, myThid )
263
264 #ifdef ALLOW_ADAMSBASHFORTH_3
265 m1 = 1 + MOD(iterNb+1,2)
266 m2 = 1 + MOD( iterNb ,2)
267 CALL GAD_CALC_RHS(
268 I bi, bj, iMin,iMax,jMin,jMax, k, kM1, kUp, kDown,
269 I xA, yA, maskUp, uFld(1-OLx,1-OLy,k),
270 I vFld(1-OLx,1-OLy,k), wFld(1-OLx,1-OLy,k),
271 I uTrans, vTrans, rTrans, rTransKp,
272 I diffKhS, diffK4S, KappaRk(1-OLx,1-OLy,k), diffKr4S,
273 I gsNm(1-OLx,1-OLy,1,1,1,m2), salt, dTtracerLev,
274 I GAD_SALINITY, saltAdvScheme, saltVertAdvScheme,
275 I calcAdvection, saltImplVertAdv, AdamsBashforth_S,
276 I saltVertDiff4, useGMRedi, useKPP,
277 O fZon, fMer,
278 U fVer, gS,
279 I myTime, myIter, myThid )
280 #else /* ALLOW_ADAMSBASHFORTH_3 */
281 CALL GAD_CALC_RHS(
282 I bi, bj, iMin,iMax,jMin,jMax, k, kM1, kUp, kDown,
283 I xA, yA, maskUp, uFld(1-OLx,1-OLy,k),
284 I vFld(1-OLx,1-OLy,k), wFld(1-OLx,1-OLy,k),
285 I uTrans, vTrans, rTrans, rTransKp,
286 I diffKhS, diffK4S, KappaRk(1-OLx,1-OLy,k), diffKr4S,
287 I gsNm1, salt, dTtracerLev,
288 I GAD_SALINITY, saltAdvScheme, saltVertAdvScheme,
289 I calcAdvection, saltImplVertAdv, AdamsBashforth_S,
290 I saltVertDiff4, useGMRedi, useKPP,
291 O fZon, fMer,
292 U fVer, gS,
293 I myTime, myIter, myThid )
294 #endif /* ALLOW_ADAMSBASHFORTH_3 */
295
296 C-- External salinity forcing term(s) inside Adams-Bashforth:
297 IF ( saltForcing .AND. tracForcingOutAB.NE.1 )
298 & CALL EXTERNAL_FORCING_S(
299 I iMin, iMax, jMin, jMax, bi, bj, k,
300 I myTime, myThid )
301
302 IF ( AdamsBashforthGs ) THEN
303 #ifdef ALLOW_ADAMSBASHFORTH_3
304 CALL ADAMS_BASHFORTH3(
305 I bi, bj, k, Nr,
306 U gS, gsNm, gs_AB,
307 I saltStartAB, iterNb, myThid )
308 #else
309 CALL ADAMS_BASHFORTH2(
310 I bi, bj, k, Nr,
311 U gS, gsNm1, gs_AB,
312 I saltStartAB, iterNb, myThid )
313 #endif
314 #ifdef ALLOW_DIAGNOSTICS
315 IF ( useDiagnostics ) THEN
316 CALL DIAGNOSTICS_FILL(gs_AB,'AB_gS ',k,1,2,bi,bj,myThid)
317 ENDIF
318 #endif /* ALLOW_DIAGNOSTICS */
319 ENDIF
320
321 C-- External salinity forcing term(s) outside Adams-Bashforth:
322 IF ( saltForcing .AND. tracForcingOutAB.EQ.1 )
323 & CALL EXTERNAL_FORCING_S(
324 I iMin, iMax, jMin, jMax, bi, bj, k,
325 I myTime, myThid )
326
327 #ifdef NONLIN_FRSURF
328 IF (nonlinFreeSurf.GT.0) THEN
329 CALL FREESURF_RESCALE_G(
330 I bi, bj, k,
331 U gS,
332 I myThid )
333 IF ( AdamsBashforthGs ) THEN
334 #ifdef ALLOW_ADAMSBASHFORTH_3
335 # ifdef ALLOW_AUTODIFF_TAMC
336 CADJ STORE gsNm(:,:,k,bi,bj,1) = comlev1_bibj_k, key=kkey,
337 CADJ & byte=isbyte, kind = isbyte
338 CADJ STORE gsNm(:,:,k,bi,bj,2) = comlev1_bibj_k, key=kkey,
339 CADJ & kind = isbyte
340 # endif
341 CALL FREESURF_RESCALE_G(
342 I bi, bj, k,
343 U gsNm(1-OLx,1-OLy,1,1,1,1),
344 I myThid )
345 CALL FREESURF_RESCALE_G(
346 I bi, bj, k,
347 U gsNm(1-OLx,1-OLy,1,1,1,2),
348 I myThid )
349 #else
350 CALL FREESURF_RESCALE_G(
351 I bi, bj, k,
352 U gsNm1,
353 I myThid )
354 #endif
355 ENDIF
356 ENDIF
357 #endif /* NONLIN_FRSURF */
358
359 #ifdef ALLOW_ADAMSBASHFORTH_3
360 IF ( AdamsBashforth_S ) THEN
361 CALL TIMESTEP_TRACER(
362 I bi, bj, k, dTtracerLev(k),
363 I gsNm(1-OLx,1-OLy,1,1,1,m2),
364 U gS,
365 I myIter, myThid )
366 ELSE
367 #endif
368 CALL TIMESTEP_TRACER(
369 I bi, bj, k, dTtracerLev(k),
370 I salt,
371 U gS,
372 I myIter, myThid )
373 #ifdef ALLOW_ADAMSBASHFORTH_3
374 ENDIF
375 #endif
376
377 C- end of vertical index (k) loop (Nr:1)
378 ENDDO
379
380 #ifdef ALLOW_SALT_PLUME
381 #ifdef SALT_PLUME_VOLUME
382 IF ( useSALT_PLUME ) THEN
383 IF ( usingZCoords ) THEN
384 CALL SALT_PLUME_APPLY(
385 I 2, bi, bj,
386 I recip_hFacC,
387 I salt,1,
388 U gS,
389 I myTime, myIter, myThid )
390 ENDIF
391 ENDIF
392 #endif /* SALT_PLUME_VOLUME */
393 #endif /* ALLOW_SALT_PLUME */
394
395 #ifdef ALLOW_DOWN_SLOPE
396 IF ( useDOWN_SLOPE ) THEN
397 IF ( usingPCoords ) THEN
398 CALL DWNSLP_APPLY(
399 I GAD_SALINITY, bi, bj, kSurfC,
400 I recip_drF, recip_hFacC, recip_rA,
401 I dTtracerLev,
402 I salt,
403 U gS,
404 I myTime, myIter, myThid )
405 ELSE
406 CALL DWNSLP_APPLY(
407 I GAD_SALINITY, bi, bj, kLowC,
408 I recip_drF, recip_hFacC, recip_rA,
409 I dTtracerLev,
410 I salt,
411 U gS,
412 I myTime, myIter, myThid )
413 ENDIF
414 ENDIF
415 #endif /* ALLOW_DOWN_SLOPE */
416
417 iMin = 0
418 iMax = sNx+1
419 jMin = 0
420 jMax = sNy+1
421
422 C-- Implicit vertical advection & diffusion
423
424 #ifdef INCLUDE_IMPLVERTADV_CODE
425 IF ( saltImplVertAdv ) THEN
426 #ifdef ALLOW_AUTODIFF_TAMC
427 CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
428 CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
429 CADJ STORE wFld(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
430 CADJ STORE salt(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
431 CADJ STORE recip_hFac(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
432 #endif /* ALLOW_AUTODIFF_TAMC */
433 CALL GAD_IMPLICIT_R(
434 I saltImplVertAdv, saltVertAdvScheme, GAD_SALINITY,
435 I dTtracerLev,
436 I kappaRk, recip_hFac, wFld, salt,
437 U gS,
438 I bi, bj, myTime, myIter, myThid )
439 ELSEIF ( implicitDiffusion ) THEN
440 #else /* INCLUDE_IMPLVERTADV_CODE */
441 IF ( implicitDiffusion ) THEN
442 #endif /* INCLUDE_IMPLVERTADV_CODE */
443 #ifdef ALLOW_AUTODIFF_TAMC
444 CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
445 CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
446 #endif /* ALLOW_AUTODIFF_TAMC */
447 CALL IMPLDIFF(
448 I bi, bj, iMin, iMax, jMin, jMax,
449 I GAD_SALINITY, kappaRk, recip_hFac,
450 U gS,
451 I myThid )
452 ENDIF
453
454 #endif /* ALLOW_GENERIC_ADVDIFF */
455
456 RETURN
457 END

  ViewVC Help
Powered by ViewVC 1.1.22