/[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.12 by heimbach, Sun Mar 24 02:12:50 2002 UTC revision 1.17 by jmc, Sat Jan 3 00:46:53 2004 UTC
# Line 34  CBOP Line 34  CBOP
34  C !ROUTINE: GAD_ADVECTION  C !ROUTINE: GAD_ADVECTION
35    
36  C !INTERFACE: ==========================================================  C !INTERFACE: ==========================================================
37        SUBROUTINE GAD_ADVECTION(bi,bj,advectionScheme,tracerIdentity,        SUBROUTINE GAD_ADVECTION(
38       U                         Tracer,Gtracer,       I           implicitAdvection, advectionScheme, tracerIdentity,
39       I                         myTime,myIter,myThid)       I           uVel, vVel, wVel, tracer,
40         O           gTracer,
41         I           bi,bj, myTime,myIter,myThid)
42    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
43    
44  C !DESCRIPTION:  C !DESCRIPTION:
45  C Calculates the tendancy of a tracer due to advection.  C Calculates the tendancy of a tracer due to advection.
# Line 62  C !USES: =============================== Line 65  C !USES: ===============================
65  #include "SIZE.h"  #include "SIZE.h"
66  #include "EEPARAMS.h"  #include "EEPARAMS.h"
67  #include "PARAMS.h"  #include "PARAMS.h"
 #include "DYNVARS.h"  
68  #include "GRID.h"  #include "GRID.h"
69  #include "GAD.h"  #include "GAD.h"
70  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 71  C !USES: =============================== Line 73  C !USES: ===============================
73  #endif  #endif
74    
75  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
76  C  bi,bj                :: tile indices  C  implicitAdvection    :: vertical advection treated implicitly (later on)
77  C  advectionScheme      :: advection scheme to use  C  advectionScheme      :: advection scheme to use
78  C  tracerIdentity       :: identifier for the tracer (required only for OBCS)  C  tracerIdentity       :: identifier for the tracer (required only for OBCS)
79  C  Tracer               :: tracer field  C  uVel                 :: velocity, zonal component
80    C  vVel                 :: velocity, meridional component
81    C  wVel                 :: velocity, vertical component
82    C  tracer               :: tracer field
83    C  bi,bj                :: tile indices
84  C  myTime               :: current time  C  myTime               :: current time
85  C  myIter               :: iteration number  C  myIter               :: iteration number
86  C  myThid               :: thread number  C  myThid               :: thread number
87        INTEGER bi,bj        LOGICAL implicitAdvection
88        INTEGER advectionScheme        INTEGER advectionScheme
89        INTEGER tracerIdentity        INTEGER tracerIdentity
90        _RL Tracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL uVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
91          _RL vVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
92          _RL wVel  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
93          _RL tracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
94          INTEGER bi,bj
95        _RL myTime        _RL myTime
96        INTEGER myIter        INTEGER myIter
97        INTEGER myThid        INTEGER myThid
# Line 128  C  ipass                :: number of the Line 138  C  ipass                :: number of the
138  CEOP  CEOP
139    
140  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
141              act0 = tracerIdentity - 1
142              max0 = maxpass
143            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
144            max1 = myBxHi(myThid) - myBxLo(myThid) + 1            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
145            act2 = bj - myByLo(myThid)            act2 = bj - myByLo(myThid)
# Line 135  CEOP Line 147  CEOP
147            act3 = myThid - 1            act3 = myThid - 1
148            max3 = nTx*nTy            max3 = nTx*nTy
149            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
150            ikey = (act1 + 1) + act2*max1            igadkey = (act0 + 1)
151       &                      + act3*max1*max2       &                      + act1*max0
152       &                      + act4*max1*max2*max3       &                      + act2*max0*max1
153         &                      + act3*max0*max1*max2
154         &                      + act4*max0*max1*max2*max3
155              if (tracerIdentity.GT.maxpass) then
156                 print *, 'ph-pass gad_advection ', maxpass, tracerIdentity
157                 STOP 'maxpass seems smaller than tracerIdentity'
158              endif
159  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
160    
161  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
# Line 166  C     uninitialised but inert locations. Line 184  C     uninitialised but inert locations.
184  C--   Start of k loop for horizontal fluxes  C--   Start of k loop for horizontal fluxes
185        DO k=1,Nr        DO k=1,Nr
186  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
187           kkey = (ikey-1)*Nr + k           kkey = (igadkey-1)*Nr + k
188  CADJ STORE tracer(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE tracer(:,:,k,bi,bj) =
189    CADJ &     comlev1_bibj_k_gad, key=kkey, byte=isbyte
190  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
191    
192  C--   Get temporary terms used by tendency routines  C--   Get temporary terms used by tendency routines
# Line 192  C--   Make local copy of tracer array Line 211  C--   Make local copy of tracer array
211    
212        IF (useCubedSphereExchange) THEN        IF (useCubedSphereExchange) THEN
213         nipass=3         nipass=3
214    #ifdef ALLOW_AUTODIFF_TAMC
215           if ( nipass.GT.maxcube )
216         &      STOP 'maxcube needs to be = 3'
217    #endif
218        ELSE        ELSE
219         nipass=1         nipass=1
220        ENDIF        ENDIF
# Line 200  cph       nipass=1 Line 223  cph       nipass=1
223  C--   Multiple passes for different directions on different tiles  C--   Multiple passes for different directions on different tiles
224        DO ipass=1,nipass        DO ipass=1,nipass
225  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
226           passkey = ipass + (k-1)   *maxpass           passkey = ipass + (k-1)      *maxcube
227       &                   + (ikey-1)*maxpass*Nr       &                   + (igadkey-1)*maxcube*Nr
228           IF (nipass .GT. maxpass) THEN           IF (nipass .GT. maxpass) THEN
229            STOP 'GAD_ADVECTION: nipass > maxpass. check tamc.h'            STOP 'GAD_ADVECTION: nipass > maxcube. check tamc.h'
230           ENDIF           ENDIF
231  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
232    
# Line 252  C-    Advective flux in X Line 275  C-    Advective flux in X
275    
276  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
277  #ifndef DISABLE_MULTIDIM_ADVECTION  #ifndef DISABLE_MULTIDIM_ADVECTION
278  CADJ STORE localTij(:,:)  = comlev1_bibj_pass, key=passkey, byte=isbyte  CADJ STORE localTij(:,:)  =
279    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
280  #endif  #endif
281  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
282    
# Line 318  C-    Advective flux in Y Line 342  C-    Advective flux in Y
342    
343  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
344  #ifndef DISABLE_MULTIDIM_ADVECTION  #ifndef DISABLE_MULTIDIM_ADVECTION
345  CADJ STORE localTij(:,:)  = comlev1_bibj_pass, key=passkey, byte=isbyte  CADJ STORE localTij(:,:)  =
346    CADJ &     comlev1_bibj_k_gad_pass, key=passkey, byte=isbyte
347  #endif  #endif
348  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
349    
# Line 375  C--   End of K loop for horizontal fluxe Line 400  C--   End of K loop for horizontal fluxe
400  C--   Start of k loop for vertical flux  C--   Start of k loop for vertical flux
401        DO k=Nr,1,-1        DO k=Nr,1,-1
402  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
403           kkey = (ikey-1)*Nr + k           kkey = (igadkey-1)*Nr + k
404  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
405    
406  C--   kup    Cycles through 1,2 to point to w-layer above  C--   kup    Cycles through 1,2 to point to w-layer above
# Line 386  c     kp1=min(Nr,k+1) Line 411  c     kp1=min(Nr,k+1)
411        kp1Msk=1.        kp1Msk=1.
412        if (k.EQ.Nr) kp1Msk=0.        if (k.EQ.Nr) kp1Msk=0.
413    
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE localTijk(:,:,k)    
 CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE rTrans(:,:)    
 CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
414  C-- Compute Vertical transport  C-- Compute Vertical transport
415  C     Note: wVel needs to be masked  C     Note: wVel needs to be masked
416    
# Line 404  C- Surface interface : Line 422  C- Surface interface :
422           rTransKp1(i,j) = rTrans(i,j)           rTransKp1(i,j) = rTrans(i,j)
423           rTrans(i,j) = 0.           rTrans(i,j) = 0.
424           fVerT(i,j,kUp) = 0.           fVerT(i,j,kUp) = 0.
425             af(i,j) = 0.
426          ENDDO          ENDDO
427         ENDDO         ENDDO
428    
# Line 425  C--   Residual transp = Bolus transp + E Line 444  C--   Residual transp = Bolus transp + E
444       &                    rTrans, bi, bj, k, myThid)       &                    rTrans, bi, bj, k, myThid)
445  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
446    
447    #ifdef ALLOW_AUTODIFF_TAMC
448    CADJ STORE localTijk(:,:,k)  
449    CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
450    CADJ STORE rTrans(:,:)  
451    CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
452    #endif /* ALLOW_AUTODIFF_TAMC */
453    
454           IF ( .NOT.implicitAdvection ) THEN
455  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
456         IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
457          CALL GAD_FLUXLIMIT_ADV_R(           CALL GAD_FLUXLIMIT_ADV_R(
458       &       bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)       &        bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
459         ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
460          CALL GAD_DST3_ADV_R(           CALL GAD_DST3_ADV_R(
461       &       bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)       &        bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
462         ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
463          CALL GAD_DST3FL_ADV_R(           CALL GAD_DST3FL_ADV_R(
464       &       bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)       &        bi,bj,k,deltaTtracer,rTrans,wVel,localTijk,af,myThid)
465         ELSE          ELSE
466          STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'           STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
467            ENDIF
468         ENDIF         ENDIF
469  C-    add the advective flux to fVerT  C-    add the advective flux to fVerT
470         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
# Line 448  C-    add the advective flux to fVerT Line 476  C-    add the advective flux to fVerT
476  C- end Surface/Interior if bloc  C- end Surface/Interior if bloc
477        ENDIF        ENDIF
478    
479    #ifdef ALLOW_AUTODIFF_TAMC
480    CADJ STORE rTrans(:,:)  
481    CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
482    CADJ STORE rTranskp1(:,:)  
483    CADJ &     = comlev1_bibj_k_gad, key=kkey, byte=isbyte
484    #endif /* ALLOW_AUTODIFF_TAMC */
485    
486  C--   Divergence of fluxes  C--   Divergence of fluxes
487        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
488         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
Legend:
Removed from v.1.12  
changed lines
  Added in v.1.17
  ViewVC Help
Powered by ViewVC 1.1.22