/[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.10 by jmc, Tue Mar 5 15:17:45 2002 UTC revision 1.11 by jmc, Wed Mar 6 02:01:54 2002 UTC
# Line 99  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 110  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 152  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 173  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 373  C--   kup    Cycles through 1,2 to point Line 383  C--   kup    Cycles through 1,2 to point
383  C--   kDown  Cycles through 2,1 to point to w-layer below  C--   kDown  Cycles through 2,1 to point to w-layer below
384        kup  = 1+MOD(k+1,2)        kup  = 1+MOD(k+1,2)
385        kDown= 1+MOD(k,2)        kDown= 1+MOD(k,2)
386    c     kp1=min(Nr,k+1)
387  C--   Get temporary terms used by tendency routines        kp1Msk=1.
388        CALL CALC_COMMON_FACTORS (        if (k.EQ.Nr) kp1Msk=0.
      I         bi,bj,iMin,iMax,jMin,jMax,k,  
      O         xA,yA,uTrans,vTrans,rTrans,maskUp,  
      I         myThid)  
   
 C-    Advective flux in R  
       DO j=1-Oly,sNy+Oly  
        DO i=1-Olx,sNx+Olx  
         af(i,j) = 0.  
        ENDDO  
       ENDDO  
389    
390  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
391  CADJ STORE localTijk(:,:,k)    CADJ STORE localTijk(:,:,k)  
392  CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
393  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
394    
395    C-- Compute Vertical transport
396  C     Note: wVel needs to be masked  C     Note: wVel needs to be masked
397        IF (K.GE.2) THEN  
398          IF (k.EQ.1) THEN
399    C- Surface interface :
400    
401           DO j=1-Oly,sNy+Oly
402            DO i=1-Olx,sNx+Olx
403             rTransKp1(i,j) = rTrans(i,j)
404             rTrans(i,j) = 0.
405             fVerT(i,j,kUp) = 0.
406            ENDDO
407           ENDDO
408    
409          ELSE
410    C- Interior interface :
411           DO j=1-Oly,sNy+Oly
412            DO i=1-Olx,sNx+Olx
413             rTransKp1(i,j) = kp1Msk*rTrans(i,j)
414             rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
415         &               *maskC(i,j,k-1,bi,bj)
416             af(i,j) = 0.
417            ENDDO
418           ENDDO
419    
420    #ifdef ALLOW_GMREDI
421    C--   Residual transp = Bolus transp + Eulerian transp
422           IF (useGMRedi)
423         &   CALL GMREDI_CALC_WFLOW(
424         &                    rTrans, bi, bj, k, myThid)
425    #endif /* ALLOW_GMREDI */
426    
427  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
428         IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN         IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
429          CALL GAD_FLUXLIMIT_ADV_R(          CALL GAD_FLUXLIMIT_ADV_R(
# Line 407  C-    Compute vertical advective flux in Line 437  C-    Compute vertical advective flux in
437         ELSE         ELSE
438          STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'          STOP 'GAD_ADVECTION: adv. scheme incompatibale with mutli-dim'
439         ENDIF         ENDIF
440  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)*tracer(i,j,k,bi,bj)  
 c    &             rTrans(i,j)*localTijk(i,j,k)  
         ENDDO  
        ENDDO  
       ELSE  
 C-    Surface "correction" term at k=1 :  
441         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
442          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
443           af(i,j) = rTrans(i,j)*tracer(i,j,k,bi,bj)           fVerT(i,j,kUp) = af(i,j)
 c        af(i,j) = rTrans(i,j)*localTijk(i,j,k)  
444          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)  
445         ENDDO         ENDDO
446        ENDDO  
447    C- end Surface/Interior if bloc
448          ENDIF
449    
450  C--   Divergence of fluxes  C--   Divergence of fluxes
       kp1=min(Nr,k+1)  
       kp1Msk=1.  
       if (k.EQ.Nr) kp1Msk=0.  
451        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
452         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
453          localTij(i,j)=localTijk(i,j,k)-deltaTtracer*          localTij(i,j)=localTijk(i,j,k)-deltaTtracer*
454       &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
455       &    *recip_rA(i,j,bi,bj)       &    *recip_rA(i,j,bi,bj)
456       &    *( fVerT(i,j,kUp)-fVerT(i,j,kDown)       &    *( fVerT(i,j,kUp)-fVerT(i,j,kDown)
457       &      -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))  
458       &     )*rkFac       &     )*rkFac
459          gTracer(i,j,k,bi,bj)=          gTracer(i,j,k,bi,bj)=
460       &   (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.10  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.22