/[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.5 by adcroft, Fri Sep 21 13:11:43 2001 UTC revision 1.12 by heimbach, Sun Mar 24 02:12:50 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 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 96  C  kdown                :: index into 2 Line 99  C  kdown                :: index into 2
99  C  kp1                  :: =k+1 for k<Nr, =Nr for k=Nr  C  kp1                  :: =k+1 for k<Nr, =Nr for k=Nr
100  C  xA,yA                :: areas of X and Y face of tracer cells  C  xA,yA                :: areas of X and Y face of tracer cells
101  C  uTrans,vTrans,rTrans :: 2-D arrays of volume transports at U,V and W points  C  uTrans,vTrans,rTrans :: 2-D arrays of volume transports at U,V and W points
102    C  rTransKp1            :: vertical volume transport at interface k+1
103  C  af                   :: 2-D array for horizontal advective flux  C  af                   :: 2-D array for horizontal advective flux
104  C  fVerT                :: 2 1/2D arrays for vertical advective flux  C  fVerT                :: 2 1/2D arrays for vertical advective flux
105  C  localTij             :: 2-D array used as temporary local copy of tracer fld  C  localTij             :: 2-D array used as temporary local copy of tracer fld
# Line 107  C  nipass               :: number of pas Line 111  C  nipass               :: number of pas
111  C  ipass                :: number of the current pass being made  C  ipass                :: number of the current pass being made
112        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113        INTEGER iMin,iMax,jMin,jMax        INTEGER iMin,iMax,jMin,jMax
114        INTEGER i,j,k,kup,kDown,kp1        INTEGER i,j,k,kup,kDown
115        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
117        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
119        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
120          _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121        _RL af      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL af      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
123        _RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 122  C  ipass                :: number of the Line 127  C  ipass                :: number of the
127        INTEGER nipass,ipass        INTEGER nipass,ipass
128  CEOP  CEOP
129    
130    #ifdef ALLOW_AUTODIFF_TAMC
131              act1 = bi - myBxLo(myThid)
132              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
133              act2 = bj - myByLo(myThid)
134              max2 = myByHi(myThid) - myByLo(myThid) + 1
135              act3 = myThid - 1
136              max3 = nTx*nTy
137              act4 = ikey_dynamics - 1
138              ikey = (act1 + 1) + act2*max1
139         &                      + act3*max1*max2
140         &                      + act4*max1*max2*max3
141    #endif /* ALLOW_AUTODIFF_TAMC */
142    
143  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
144  C     These inital values do not alter the numerical results. They  C     These inital values do not alter the numerical results. They
145  C     just ensure that all memory references are to valid floating  C     just ensure that all memory references are to valid floating
# Line 136  C     uninitialised but inert locations. Line 154  C     uninitialised but inert locations.
154          rTrans(i,j)  = 0. _d 0          rTrans(i,j)  = 0. _d 0
155          fVerT(i,j,1) = 0. _d 0          fVerT(i,j,1) = 0. _d 0
156          fVerT(i,j,2) = 0. _d 0          fVerT(i,j,2) = 0. _d 0
157            rTransKp1(i,j)= 0. _d 0
158         ENDDO         ENDDO
159        ENDDO        ENDDO
160    
# Line 146  C     uninitialised but inert locations. Line 165  C     uninitialised but inert locations.
165    
166  C--   Start of k loop for horizontal fluxes  C--   Start of k loop for horizontal fluxes
167        DO k=1,Nr        DO k=1,Nr
168    #ifdef ALLOW_AUTODIFF_TAMC
169             kkey = (ikey-1)*Nr + k
170    CADJ STORE tracer(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
171    #endif /* ALLOW_AUTODIFF_TAMC */
172    
173  C--   Get temporary terms used by tendency routines  C--   Get temporary terms used by tendency routines
174        CALL CALC_COMMON_FACTORS (        CALL CALC_COMMON_FACTORS (
# Line 153  C--   Get temporary terms used by tenden Line 176  C--   Get temporary terms used by tenden
176       O         xA,yA,uTrans,vTrans,rTrans,maskUp,       O         xA,yA,uTrans,vTrans,rTrans,maskUp,
177       I         myThid)       I         myThid)
178    
179    #ifdef ALLOW_GMREDI
180    C--   Residual transp = Bolus transp + Eulerian transp
181           IF (useGMRedi)
182         &   CALL GMREDI_CALC_UVFLOW(
183         &            uTrans, vTrans, bi, bj, k, myThid)
184    #endif /* ALLOW_GMREDI */
185    
186  C--   Make local copy of tracer array  C--   Make local copy of tracer array
187        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
188         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
# Line 165  C--   Make local copy of tracer array Line 195  C--   Make local copy of tracer array
195        ELSE        ELSE
196         nipass=1         nipass=1
197        ENDIF        ENDIF
198         nipass=1  cph       nipass=1
199    
200  C--   Multiple passes for different directions on different tiles  C--   Multiple passes for different directions on different tiles
201        DO ipass=1,nipass        DO ipass=1,nipass
202    #ifdef ALLOW_AUTODIFF_TAMC
203             passkey = ipass + (k-1)   *maxpass
204         &                   + (ikey-1)*maxpass*Nr
205             IF (nipass .GT. maxpass) THEN
206              STOP 'GAD_ADVECTION: nipass > maxpass. check tamc.h'
207             ENDIF
208    #endif /* ALLOW_AUTODIFF_TAMC */
209    
210        IF (nipass.EQ.3) THEN        IF (nipass.EQ.3) THEN
211         calc_fluxes_X=.FALSE.         calc_fluxes_X=.FALSE.
# Line 212  C-    Advective flux in X Line 249  C-    Advective flux in X
249          af(i,j) = 0.          af(i,j) = 0.
250         ENDDO         ENDDO
251        ENDDO        ENDDO
252    
253    #ifdef ALLOW_AUTODIFF_TAMC
254    #ifndef DISABLE_MULTIDIM_ADVECTION
255    CADJ STORE localTij(:,:)  = comlev1_bibj_pass, key=passkey, byte=isbyte
256    #endif
257    #endif /* ALLOW_AUTODIFF_TAMC */
258    
259        IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN        IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
260         CALL GAD_FLUXLIMIT_ADV_X(         CALL GAD_FLUXLIMIT_ADV_X(
261       &      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 266  C-    Advective flux in X
266         CALL GAD_DST3FL_ADV_X(         CALL GAD_DST3FL_ADV_X(
267       &       bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)       &       bi,bj,k,deltaTtracer,uTrans,uVel,localTij,af,myThid)
268        ELSE        ELSE
269         STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'         STOP 'GAD_ADVECTION: adv. scheme incompatibale with multi-dim'
270        ENDIF        ENDIF
271    
272        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
273         DO i=1-Olx,sNx+Olx-1         DO i=1-Olx,sNx+Olx-1
274          localTij(i,j)=localTij(i,j)-deltaTtracer*          localTij(i,j)=localTij(i,j)-deltaTtracer*
# Line 270  C-    Advective flux in Y Line 315  C-    Advective flux in Y
315          af(i,j) = 0.          af(i,j) = 0.
316         ENDDO         ENDDO
317        ENDDO        ENDDO
318    
319    #ifdef ALLOW_AUTODIFF_TAMC
320    #ifndef DISABLE_MULTIDIM_ADVECTION
321    CADJ STORE localTij(:,:)  = comlev1_bibj_pass, key=passkey, byte=isbyte
322    #endif
323    #endif /* ALLOW_AUTODIFF_TAMC */
324    
325        IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN        IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
326         CALL GAD_FLUXLIMIT_ADV_Y(         CALL GAD_FLUXLIMIT_ADV_Y(
327       &       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 334  C-    Advective flux in Y
334        ELSE        ELSE
335         STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'         STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
336        ENDIF        ENDIF
337    
338        DO j=1-Oly,sNy+Oly-1        DO j=1-Oly,sNy+Oly-1
339         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
340          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 374  C--   End of K loop for horizontal fluxe
374    
375  C--   Start of k loop for vertical flux  C--   Start of k loop for vertical flux
376        DO k=Nr,1,-1        DO k=Nr,1,-1
377    #ifdef ALLOW_AUTODIFF_TAMC
378             kkey = (ikey-1)*Nr + k
379    #endif /* ALLOW_AUTODIFF_TAMC */
380    
381  C--   kup    Cycles through 1,2 to point to w-layer above  C--   kup    Cycles through 1,2 to point to w-layer above
382  C--   kDown  Cycles through 2,1 to point to w-layer below  C--   kDown  Cycles through 2,1 to point to w-layer below
383        kup  = 1+MOD(k+1,2)        kup  = 1+MOD(k+1,2)
384        kDown= 1+MOD(k,2)        kDown= 1+MOD(k,2)
385    c     kp1=min(Nr,k+1)
386          kp1Msk=1.
387          if (k.EQ.Nr) kp1Msk=0.
388    
389  C--   Get temporary terms used by tendency routines  #ifdef ALLOW_AUTODIFF_TAMC
390        CALL CALC_COMMON_FACTORS (  CADJ STORE localTijk(:,:,k)  
391       I         bi,bj,iMin,iMax,jMin,jMax,k,  CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
392       O         xA,yA,uTrans,vTrans,rTrans,maskUp,  CADJ STORE rTrans(:,:)  
393       I         myThid)  CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
394    #endif /* ALLOW_AUTODIFF_TAMC */
395    
396  C-    Advective flux in R  C-- Compute Vertical transport
397        DO j=1-Oly,sNy+Oly  C     Note: wVel needs to be masked
398         DO i=1-Olx,sNx+Olx  
399          af(i,j) = 0.        IF (k.EQ.1) THEN
400    C- Surface interface :
401    
402           DO j=1-Oly,sNy+Oly
403            DO i=1-Olx,sNx+Olx
404             rTransKp1(i,j) = rTrans(i,j)
405             rTrans(i,j) = 0.
406             fVerT(i,j,kUp) = 0.
407            ENDDO
408         ENDDO         ENDDO
       ENDDO  
409    
410  C     Note: wVel needs to be masked        ELSE
411        IF (K.GE.2) THEN  C- Interior interface :
412           DO j=1-Oly,sNy+Oly
413            DO i=1-Olx,sNx+Olx
414             rTransKp1(i,j) = kp1Msk*rTrans(i,j)
415             rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
416         &               *maskC(i,j,k-1,bi,bj)
417             af(i,j) = 0.
418            ENDDO
419           ENDDO
420    
421    #ifdef ALLOW_GMREDI
422    C--   Residual transp = Bolus transp + Eulerian transp
423           IF (useGMRedi)
424         &   CALL GMREDI_CALC_WFLOW(
425         &                    rTrans, bi, bj, k, myThid)
426    #endif /* ALLOW_GMREDI */
427    
428  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
429         IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN         IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
430          CALL GAD_FLUXLIMIT_ADV_R(          CALL GAD_FLUXLIMIT_ADV_R(
# Line 355  C-    Compute vertical advective flux in Line 438  C-    Compute vertical advective flux in
438         ELSE         ELSE
439          STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'          STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
440         ENDIF         ENDIF
441  C-    Surface "correction" term at k>1 :  C-    add the advective flux to fVerT
        DO j=1-Oly,sNy+Oly  
         DO i=1-Olx,sNx+Olx  
          af(i,j) = af(i,j)  
      &           + (maskC(i,j,k,bi,bj)-maskC(i,j,k-1,bi,bj))*  
      &             rTrans(i,j)*localTijk(i,j,k)  
         ENDDO  
        ENDDO  
       ELSE  
 C-    Surface "correction" term at k=1 :  
442         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
443          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
444           af(i,j) = rTrans(i,j)*localTijk(i,j,k)           fVerT(i,j,kUp) = af(i,j)
445          ENDDO          ENDDO
        ENDDO  
       ENDIF  
 C-    add the advective flux to fVerT  
       DO j=1-Oly,sNy+Oly  
        DO i=1-Olx,sNx+Olx  
         fVerT(i,j,kUp) = af(i,j)  
446         ENDDO         ENDDO
447        ENDDO  
448    C- end Surface/Interior if bloc
449          ENDIF
450    
451  C--   Divergence of fluxes  C--   Divergence of fluxes
       kp1=min(Nr,k+1)  
       kp1Msk=1.  
       if (k.EQ.Nr) kp1Msk=0.  
452        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
453         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
454          localTij(i,j)=localTijk(i,j,k)-deltaTtracer*          localTij(i,j)=localTijk(i,j,k)-deltaTtracer*
455       &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
456       &    *recip_rA(i,j,bi,bj)       &    *recip_rA(i,j,bi,bj)
457       &    *( fVerT(i,j,kUp)-fVerT(i,j,kDown)       &    *( fVerT(i,j,kUp)-fVerT(i,j,kDown)
458       &      -tracer(i,j,k,bi,bj)*rA(i,j,bi,bj)*       &      -tracer(i,j,k,bi,bj)*(rTrans(i,j)-rTransKp1(i,j))
      &        (wVel(i,j,k,bi,bj)-kp1Msk*wVel(i,j,kp1,bi,bj))  
459       &     )*rkFac       &     )*rkFac
460          gTracer(i,j,k,bi,bj)=          gTracer(i,j,k,bi,bj)=
461       &   (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer       &   (localTij(i,j)-tracer(i,j,k,bi,bj))/deltaTtracer
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.22