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

Annotation 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 - (hide annotations) (download)
Fri May 2 05:35:22 2014 UTC (11 years, 3 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 atn 1.3 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 atn 1.1 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 atn 1.2 #ifdef ALLOW_SALT_PLUME
13     # include "SALT_PLUME_OPTIONS.h"
14     #endif
15 atn 1.1
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 atn 1.3 I salt,1,
388 atn 1.1 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