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

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

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

revision 1.1 by jmc, Sat Oct 22 19:56:32 2005 UTC revision 1.4 by jmc, Mon Jun 19 14:40:43 2006 UTC
# Line 7  CBOP Line 7  CBOP
7  C !ROUTINE: GAD_DST2U1_ADV_R  C !ROUTINE: GAD_DST2U1_ADV_R
8    
9  C !INTERFACE: ==========================================================  C !INTERFACE: ==========================================================
10        SUBROUTINE GAD_DST2U1_ADV_R(        SUBROUTINE GAD_DST2U1_ADV_R(
11       I           bi,bj,k, advectionScheme,       I           bi,bj,k, advectionScheme,
12       I           deltaTloc, rTrans, wVel,       I           deltaTloc, rTrans, wFld,
13       I           tracer,       I           tracer,
14       O           wT,       O           wT,
15       I           myThid )       I           myThid )
16    
17  C !DESCRIPTION:  C !DESCRIPTION:
18  C  Calculates the area integrated vertical flux due to advection  C  Calculates the area integrated vertical flux due to advection
19  C  of a tracer using second-order Direct Space and Time (DST-2)  C  of a tracer using second-order Direct Space and Time (DST-2)
20  C  interpolation (=Lax-Wendroff) or simple 1rst order upwind scheme.  C  interpolation (=Lax-Wendroff) or simple 1rst order upwind scheme.
21    
22  C !USES: ===============================================================  C !USES: ===============================================================
# Line 32  C  advectionScheme   :: advection scheme Line 32  C  advectionScheme   :: advection scheme
32  C                                                or 1rst Order Upwind  C                                                or 1rst Order Upwind
33  C  deltaTloc         :: local time-step (s)  C  deltaTloc         :: local time-step (s)
34  C  rTrans            :: vertical volume transport  C  rTrans            :: vertical volume transport
35  C  wVel              :: vertical flow  C  wFld              :: vertical flow
36  C  tracer            :: tracer field  C  tracer            :: tracer field
37  C  myThid            :: thread number  C  myThid            :: thread number
38        INTEGER bi,bj,k        INTEGER bi,bj,k
39        INTEGER advectionScheme        INTEGER advectionScheme
40        _RL deltaTloc        _RL deltaTloc
41        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
42        _RL wVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL wFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
43        _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
44        INTEGER myThid        INTEGER myThid
45    
# Line 51  C !LOCAL VARIABLES: ==================== Line 51  C !LOCAL VARIABLES: ====================
51  C  i,j               :: loop indices  C  i,j               :: loop indices
52  C  km1               :: =max( k-1 , 1 )  C  km1               :: =max( k-1 , 1 )
53  C  rLimit            :: centered (vs upwind) fraction  C  rLimit            :: centered (vs upwind) fraction
54  C  wFld              :: velocity, vertical component  C  wLoc              :: velocity, vertical component
55  C  wCFL              :: Courant-Friedrich-Levy number  C  wCFL              :: Courant-Friedrich-Levy number
56        INTEGER i,j,km1        INTEGER i,j,km1
57        _RL wFld, wCFL, rLimit        _RL wLoc, wCFL, rLimit, wAbs
58  CEOP  CEOP
59    
60        rLimit = 0. _d 0        rLimit = 0. _d 0
# Line 72  CEOP Line 72  CEOP
72         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
73          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
74    
75  c        wFld = wVel(i,j,k,bi,bj)           wLoc = wFld(i,j)
76           wFld = rTrans(i,j)*recip_rA(i,j,bi,bj)  c        wLoc = rTrans(i,j)*recip_rA(i,j,bi,bj)
77           wCFL = ABS(wFld*deltaTloc*recip_drC(k))           wCFL = ABS(wLoc*deltaTloc*recip_drC(k))
78    
79           wT(i,j) = maskC(i,j,kM1,bi,bj)*(  c        wT(i,j) = maskC(i,j,km1,bi,bj)*(
80       &     rTrans(i,j)*(Tracer(i,j,km1)+Tracer(i,j,k))*0.5 _d 0  c    &     rTrans(i,j)*(tracer(i,j,km1)+tracer(i,j,k))*0.5 _d 0
81       &   + ( 1. _d 0 - rLimit*(1. _d 0 - wCFL) )*ABS(rTrans(i,j))  c    &   + ( 1. _d 0 - rLimit*(1. _d 0 - wCFL) )*ABS(rTrans(i,j))
82       &                *(tracer(i,j,km1)-tracer(i,j,k))*0.5 _d 0*rkSign  c    &                *(tracer(i,j,km1)-tracer(i,j,k))*0.5 _d 0*rkSign
83    c    &                                  )
84    C-- above formulation produces large truncation error when:
85    C    1rst.O upWind and   w*rkSign > 0 & |tracer(k-1)| << |tracer(k)|
86    C                   or   w*rkSign < 0 & |tracer(k-1)| >> |tracer(k)|
87    C-- change to a more robust expression:
88             wAbs = ABS(rTrans(i,j))*rkSign
89         &       *( 1. _d 0 - rLimit*(1. _d 0 - wCFL) )
90             wT(i,j) = maskC(i,j,km1,bi,bj)*(
91         &             ( rTrans(i,j)+wAbs )* 0.5 _d 0 * tracer(i,j,km1)
92         &           + ( rTrans(i,j)-wAbs )* 0.5 _d 0 * tracer(i,j,k)
93       &                                  )       &                                  )
94          ENDDO          ENDDO
95         ENDDO         ENDDO
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22