/[MITgcm]/MITgcm/pkg/generic_advdiff/gad_advection.F
ViewVC logotype

Diff of /MITgcm/pkg/generic_advdiff/gad_advection.F

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

revision 1.4 by adcroft, Wed Sep 19 20:45:09 2001 UTC revision 1.9 by adcroft, Mon Mar 4 17:26:41 2002 UTC
# Line 4  C $Name$ Line 4  C $Name$
4  CBOI  CBOI
5  C !TITLE: pkg/generic\_advdiff  C !TITLE: pkg/generic\_advdiff
6  C !AUTHORS: adcroft@mit.edu  C !AUTHORS: adcroft@mit.edu
7  C !INTRODUCTION:  C !INTRODUCTION: Generic Advection Diffusion Package
 C \section{Generica Advection Diffusion Package}  
8  C  C
9  C Package "generic\_advdiff" provides a common set of routines for calculating  C Package "generic\_advdiff" provides a common set of routines for calculating
10  C advective/diffusive fluxes for tracers (cell centered quantities on a C-grid).  C advective/diffusive fluxes for tracers (cell centered quantities on a C-grid).
# Line 48  C Line 47  C
47  C The algorithm is as follows:  C The algorithm is as follows:
48  C \begin{itemize}  C \begin{itemize}
49  C \item{$\theta^{(n+1/3)} = \theta^{(n)}  C \item{$\theta^{(n+1/3)} = \theta^{(n)}
50  C      - \Delta t \partial_x (u\theta) + \theta \partial_x u$}  C      - \Delta t \partial_x (u\theta^{(n)}) + \theta^{(n)} \partial_x u$}
51  C \item{$\theta^{(n+2/3)} = \theta^{(n+1/3)}  C \item{$\theta^{(n+2/3)} = \theta^{(n+1/3)}
52  C      - \Delta t \partial_y (v\theta) + \theta \partial_y v$}  C      - \Delta t \partial_y (v\theta^{(n+1/3)}) + \theta^{(n)} \partial_y v$}
53  C \item{$\theta^{(n+3/3)} = \theta^{(n+2/3)}  C \item{$\theta^{(n+3/3)} = \theta^{(n+2/3)}
54  C      - \Delta t \partial_r (w\theta) + \theta \partial_r w$}  C      - \Delta t \partial_r (w\theta^{(n+2/3)}) + \theta^{(n)} \partial_r w$}
55  C \item{$G_\theta = ( \theta^{(n+3/3)} - \theta^{(n)} )/\Delta t$}  C \item{$G_\theta = ( \theta^{(n+3/3)} - \theta^{(n)} )/\Delta t$}
56  C \end{itemize}  C \end{itemize}
57  C  C
# Line 66  C !USES: =============================== Line 65  C !USES: ===============================
65  #include "DYNVARS.h"  #include "DYNVARS.h"
66  #include "GRID.h"  #include "GRID.h"
67  #include "GAD.h"  #include "GAD.h"
68    #ifdef ALLOW_AUTODIFF_TAMC
69    # include "tamc.h"
70    # include "tamc_keys.h"
71    #endif
72    
73  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
74  C  bi,bj                :: tile indices  C  bi,bj                :: tile indices
# Line 78  C  myThid               :: thread number Line 81  C  myThid               :: thread number
81        INTEGER bi,bj        INTEGER bi,bj
82        INTEGER advectionScheme        INTEGER advectionScheme
83        INTEGER tracerIdentity        INTEGER tracerIdentity
84        _RL Gtracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL Tracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
85        _RL myTime        _RL myTime
86        INTEGER myIter        INTEGER myIter
87        INTEGER myThid        INTEGER myThid
88    
89  C !OUTPUT PARAMETERS: ==================================================  C !OUTPUT PARAMETERS: ==================================================
90  C  gTracer              :: tendancy array  C  gTracer              :: tendancy array
91        _RL Tracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL gTracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
92    
93  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
94  C  maskUp               :: 2-D array for mask at W points  C  maskUp               :: 2-D array for mask at W points
# Line 122  C  ipass                :: number of the Line 125  C  ipass                :: number of the
125        INTEGER nipass,ipass        INTEGER nipass,ipass
126  CEOP  CEOP
127    
128    #ifdef ALLOW_AUTODIFF_TAMC
129              act1 = bi - myBxLo(myThid)
130              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
131              act2 = bj - myByLo(myThid)
132              max2 = myByHi(myThid) - myByLo(myThid) + 1
133              act3 = myThid - 1
134              max3 = nTx*nTy
135              act4 = ikey_dynamics - 1
136              ikey = (act1 + 1) + act2*max1
137         &                      + act3*max1*max2
138         &                      + act4*max1*max2*max3
139    #endif /* ALLOW_AUTODIFF_TAMC */
140    
141  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
142  C     These inital values do not alter the numerical results. They  C     These inital values do not alter the numerical results. They
143  C     just ensure that all memory references are to valid floating  C     just ensure that all memory references are to valid floating
# Line 146  C     uninitialised but inert locations. Line 162  C     uninitialised but inert locations.
162    
163  C--   Start of k loop for horizontal fluxes  C--   Start of k loop for horizontal fluxes
164        DO k=1,Nr        DO k=1,Nr
165    #ifdef ALLOW_AUTODIFF_TAMC
166             kkey = (ikey-1)*Nr + k
167    CADJ STORE tracer(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
168    #endif /* ALLOW_AUTODIFF_TAMC */
169    
170  C--   Get temporary terms used by tendency routines  C--   Get temporary terms used by tendency routines
171        CALL CALC_COMMON_FACTORS (        CALL CALC_COMMON_FACTORS (
# Line 165  C--   Make local copy of tracer array Line 185  C--   Make local copy of tracer array
185        ELSE        ELSE
186         nipass=1         nipass=1
187        ENDIF        ENDIF
188         nipass=1  cph       nipass=1
189    
190  C--   Multiple passes for different directions on different tiles  C--   Multiple passes for different directions on different tiles
191        DO ipass=1,nipass        DO ipass=1,nipass
192    #ifdef ALLOW_AUTODIFF_TAMC
193             passkey = ipass + (k-1)   *maxpass
194         &                   + (ikey-1)*maxpass*Nr
195             IF (nipass .GT. maxpass) THEN
196              STOP 'GAD_ADVECTION: nipass > maxpass. check tamc.h'
197             ENDIF
198    #endif /* ALLOW_AUTODIFF_TAMC */
199    
200        IF (nipass.EQ.3) THEN        IF (nipass.EQ.3) THEN
201         calc_fluxes_X=.FALSE.         calc_fluxes_X=.FALSE.
# Line 212  C-    Advective flux in X Line 239  C-    Advective flux in X
239          af(i,j) = 0.          af(i,j) = 0.
240         ENDDO         ENDDO
241        ENDDO        ENDDO
242    
243    #ifdef ALLOW_AUTODIFF_TAMC
244    #ifndef DISABLE_MULTIDIM_ADVECTION
245    CADJ STORE localTij(:,:)  = comlev1_bibj_pass, key=passkey, byte=isbyte
246    #endif
247    #endif /* ALLOW_AUTODIFF_TAMC */
248    
249        IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN        IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
250         CALL GAD_FLUXLIMIT_ADV_X(         CALL GAD_FLUXLIMIT_ADV_X(
251       &      bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)       &      bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)
# Line 222  C-    Advective flux in X Line 256  C-    Advective flux in X
256         CALL GAD_DST3FL_ADV_X(         CALL GAD_DST3FL_ADV_X(
257       &       bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)       &       bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)
258        ELSE        ELSE
259         STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'         write(0,*) advectionScheme
260           STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'
261        ENDIF        ENDIF
262    
263        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
264         DO i=1-Olx,sNx+Olx-1         DO i=1-Olx,sNx+Olx-1
265          localTij(i,j)=localTij(i,j)-deltaTtracer*          localTij(i,j)=localTij(i,j)-deltaTtracer*
# Line 270  C-    Advective flux in Y Line 306  C-    Advective flux in Y
306          af(i,j) = 0.          af(i,j) = 0.
307         ENDDO         ENDDO
308        ENDDO        ENDDO
309    
310    #ifdef ALLOW_AUTODIFF_TAMC
311    #ifndef DISABLE_MULTIDIM_ADVECTION
312    CADJ STORE localTij(:,:)  = comlev1_bibj_pass, key=passkey, byte=isbyte
313    #endif
314    #endif /* ALLOW_AUTODIFF_TAMC */
315    
316        IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN        IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
317         CALL GAD_FLUXLIMIT_ADV_Y(         CALL GAD_FLUXLIMIT_ADV_Y(
318       &       bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)       &       bi,bj,k,deltaTtracer,vTrans,vVel,localTij,af,myThid)
# Line 282  C-    Advective flux in Y Line 325  C-    Advective flux in Y
325        ELSE        ELSE
326         STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'         STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
327        ENDIF        ENDIF
328    
329        DO j=1-Oly,sNy+Oly-1        DO j=1-Oly,sNy+Oly-1
330         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
331          localTij(i,j)=localTij(i,j)-deltaTtracer*          localTij(i,j)=localTij(i,j)-deltaTtracer*
# Line 321  C--   End of K loop for horizontal fluxe Line 365  C--   End of K loop for horizontal fluxe
365    
366  C--   Start of k loop for vertical flux  C--   Start of k loop for vertical flux
367        DO k=Nr,1,-1        DO k=Nr,1,-1
368    #ifdef ALLOW_AUTODIFF_TAMC
369             kkey = (ikey-1)*Nr + k
370    #endif /* ALLOW_AUTODIFF_TAMC */
371    
372  C--   kup    Cycles through 1,2 to point to w-layer above  C--   kup    Cycles through 1,2 to point to w-layer above
373  C--   kDown  Cycles through 2,1 to point to w-layer below  C--   kDown  Cycles through 2,1 to point to w-layer below
# Line 340  C-    Advective flux in R Line 387  C-    Advective flux in R
387         ENDDO         ENDDO
388        ENDDO        ENDDO
389    
390    #ifdef ALLOW_AUTODIFF_TAMC
391    CADJ STORE localTijk(:,:,k)  
392    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
393    #endif /* ALLOW_AUTODIFF_TAMC */
394    
395  C     Note: wVel needs to be masked  C     Note: wVel needs to be masked
396        IF (K.GE.2) THEN        IF (K.GE.2) THEN
397  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.22