/[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.1 - (hide annotations) (download)
Tue Apr 22 04:32:25 2014 UTC (11 years, 3 months ago) by atn
Branch: MAIN
in progress:
1. sort out bi,bj loop
2. skip remove saltplumeflux in external_forcing_surf, (thus) skip kpp
3. move saltplume outside k-loop in [salt,temp]_integrate.F
4. add diagnostics

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

  ViewVC Help
Powered by ViewVC 1.1.22