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

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

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

revision 1.4 by jmc, Sat Dec 4 00:22:25 2004 UTC revision 1.5 by jmc, Sat Dec 4 18:50:34 2004 UTC
# Line 42  C     c3d          :: upper diagonal of Line 42  C     c3d          :: upper diagonal of
42  C     myThid       :: thread number  C     myThid       :: thread number
43        INTEGER bi,bj,k        INTEGER bi,bj,k
44        INTEGER iMin,iMax,jMin,jMax        INTEGER iMin,iMax,jMin,jMax
45        _RL deltaTarg        _RL deltaTarg(Nr)
46        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
47        _RL tFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL tFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
48        _RL a3d   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL a3d   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
# Line 65  C     rUpwind      :: upwind   contribut Line 65  C     rUpwind      :: upwind   contribut
65        _RL Cr,Rjm,Rj,Rjp, w_CFL        _RL Cr,Rjm,Rj,Rjp, w_CFL
66        _RL upwindFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL upwindFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
67        _RL rCenter, rUpwind        _RL rCenter, rUpwind
68          _RL deltaTcfl
69    
70  C Statement function provides Limiter(Cr)  C Statement function provides Limiter(Cr)
71  #include "GAD_FLUX_LIMITER.h"  #include "GAD_FLUX_LIMITER.h"
# Line 74  CEOP Line 75  CEOP
75        km1=MAX(1,k-1)        km1=MAX(1,k-1)
76        kp1=MIN(Nr,k+1)        kp1=MIN(Nr,k+1)
77    
78        IF ( k.GT.Nr .OR. k.LT.2 ) RETURN  C--   process interior interface only:
79          IF ( k.GT.1 .AND. k.LE.Nr ) THEN
80    
81  C--   Compute the upwind fraction:  C--   Compute the upwind fraction:
82           deltaTcfl = deltaTarg(k)
83         DO j=jMin,jMax         DO j=jMin,jMax
84          DO i=iMin,iMax          DO i=iMin,iMax
85           w_CFL = rTrans(i,j)*recip_rA(i,j,bi,bj)*deltaTarg*recip_drC(k)           w_CFL = deltaTcfl*rTrans(i,j)*recip_rA(i,j,bi,bj)*recip_drC(k)
86           Rjp=(tFld(i,j,kp1)-tFld(i,j,k)  )*maskC(i,j,kp1,bi,bj)           Rjp=(tFld(i,j,kp1)-tFld(i,j,k)  )*maskC(i,j,kp1,bi,bj)
87           Rj =(tFld(i,j,k)  -tFld(i,j,km1))           Rj =(tFld(i,j,k)  -tFld(i,j,km1))
88           Rjm=(tFld(i,j,km1)-tFld(i,j,km2))*maskC(i,j,km2,bi,bj)           Rjm=(tFld(i,j,km1)-tFld(i,j,km2))*maskC(i,j,km2,bi,bj)
# Line 102  C--   Compute the upwind fraction: Line 105  C--   Compute the upwind fraction:
105  C--    Add centered & upwind contributions  C--    Add centered & upwind contributions
106         DO j=jMin,jMax         DO j=jMin,jMax
107           DO i=iMin,iMax           DO i=iMin,iMax
108             rCenter = 0.5 _d 0 *deltaTarg*rTrans(i,j)             rCenter = 0.5 _d 0 *rTrans(i,j)*recip_rA(i,j,bi,bj)*rkFac
      &                        *recip_rA(i,j,bi,bj)*rkFac  
109             rUpwind = abs(rCenter)*upwindFac(i,j)             rUpwind = abs(rCenter)*upwindFac(i,j)
110             a3d(i,j,k)   = a3d(i,j,k)             a3d(i,j,k)   = a3d(i,j,k)
111       &                  + (rCenter-rUpwind)       &                  + (rCenter-rUpwind)*deltaTarg(k)
112       &                   *recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &                   *recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
113             b3d(i,j,k)   = b3d(i,j,k)             b3d(i,j,k)   = b3d(i,j,k)
114       &                  + (rCenter+rUpwind)       &                  + (rCenter+rUpwind)*deltaTarg(k)
115       &                    *recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &                    *recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
116             b3d(i,j,k-1) = b3d(i,j,k-1)             b3d(i,j,k-1) = b3d(i,j,k-1)
117       &                  - (rCenter-rUpwind)       &                  - (rCenter-rUpwind)*deltaTarg(k-1)
118       &                    *recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)       &                    *recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)
119             c3d(i,j,k-1) = c3d(i,j,k-1)             c3d(i,j,k-1) = c3d(i,j,k-1)
120       &                  - (rCenter+rUpwind)       &                  - (rCenter+rUpwind)*deltaTarg(k-1)
121       &                    *recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)       &                    *recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)
122           ENDDO           ENDDO
123         ENDDO         ENDDO
124    
125    C--   process interior interface only: end
126          ENDIF
127    
128        RETURN        RETURN
129        END        END
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.22