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

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

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

revision 1.2 by heimbach, Wed Jun 7 01:55:14 2006 UTC revision 1.5 by jmc, Mon Jun 19 14:40:43 2006 UTC
# Line 7  CBOP Line 7  CBOP
7  C !ROUTINE: GAD_DST2U1_ADV_Y  C !ROUTINE: GAD_DST2U1_ADV_Y
8    
9  C !INTERFACE: ==========================================================  C !INTERFACE: ==========================================================
10        SUBROUTINE GAD_DST2U1_ADV_Y(        SUBROUTINE GAD_DST2U1_ADV_Y(
11       I           bi,bj,k, advectionScheme, deltaTloc,       I           bi,bj,k, advectionScheme, deltaTloc,
12       I           vTrans, vVel,       I           vTrans, vFld,
13       I           tracer,       I           tracer,
14       O           vT,       O           vT,
15       I           myThid )       I           myThid )
16    
17  C !DESCRIPTION:  C !DESCRIPTION:
18  C  Calculates the area integrated meridional flux due to advection  C  Calculates the area integrated meridional 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 31  C  k                 :: vertical level Line 31  C  k                 :: vertical level
31  C  advectionScheme   :: advection scheme to use: either 2nd Order DST  C  advectionScheme   :: advection scheme to use: either 2nd Order DST
32  C                                                or 1rst Order Upwind  C                                                or 1rst Order Upwind
33  C  vTrans            :: meridional volume transport  C  vTrans            :: meridional volume transport
34  C  vVel              :: meridional flow  C  vFld              :: meridional flow
35  C  tracer            :: tracer field  C  tracer            :: tracer field
36  C  myThid            :: thread number  C  myThid            :: thread number
37        INTEGER bi,bj, k, advectionScheme        INTEGER bi,bj, k, advectionScheme
38        _RL deltaTloc        _RL deltaTloc
39        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
40        _RL vVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL vFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
41        _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
42        INTEGER myThid        INTEGER myThid
43    
# Line 48  C  vT                :: meridional advec Line 48  C  vT                :: meridional advec
48  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
49  C  i,j               :: loop indices  C  i,j               :: loop indices
50  C  yLimit            :: centered (vs upwind) fraction  C  yLimit            :: centered (vs upwind) fraction
51  C  vFld              :: velocity [m/s], meridional component  C  vLoc              :: velocity [m/s], meridional component
52  C  vCFL              :: Courant-Friedrich-Levy number  C  vCFL              :: Courant-Friedrich-Levy number
53        INTEGER i,j        INTEGER i,j
54        _RL vFld, vCFL, yLimit        _RL vLoc, vCFL, yLimit, vAbs
55  CEOP  CEOP
56    
57        yLimit = 0. _d 0        yLimit = 0. _d 0
# Line 63  CEOP Line 63  CEOP
63        DO j=1-Oly+1,sNy+Oly        DO j=1-Oly+1,sNy+Oly
64         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
65    
66  c       vFld = vVel(i,j,k,bi,bj)          vLoc = vFld(i,j)
67          vFld = vTrans(i,j)*recip_dxG(i,j,bi,bj)  c       vLoc = vTrans(i,j)*recip_dxG(i,j,bi,bj)
68       &        *recip_drF(k)*_recip_hFacS(i,j,k,bi,bj)  c    &        *recip_drF(k)*_recip_hFacS(i,j,k,bi,bj)
69          vCFL = ABS(vFld*deltaTloc*recip_dyC(i,j,bi,bj))          vCFL = ABS(vLoc*deltaTloc*recip_dyC(i,j,bi,bj))
70    
71          vT(i,j) =  c       vT(i,j) =
72       &     vTrans(i,j)*(Tracer(i,j-1)+Tracer(i,j))*0.5 _d 0  c    &     vTrans(i,j)*(tracer(i,j-1)+tracer(i,j))*0.5 _d 0
73       &   + ( 1. _d 0 - yLimit*(1. _d 0 - vCFL) )*ABS(vTrans(i,j))  c    &   + ( 1. _d 0 - yLimit*(1. _d 0 - vCFL) )*ABS(vTrans(i,j))
74       &                *(tracer(i,j-1)-tracer(i,j))*0.5 _d 0  c    &                *(tracer(i,j-1)-tracer(i,j))*0.5 _d 0
75    C-- above formulation produces large truncation error when:
76    C    1rst.O upWind and   v > 0 & |tracer(i,j-1)| << |tracer(i,j)|
77    C                   or   v < 0 & |tracer(i,j-1)| >> |tracer(i,j)|
78    C-- change to a more robust expression:
79            vAbs = ABS(vTrans(i,j))
80         &       *( 1. _d 0 - yLimit*(1. _d 0 - vCFL) )
81            vT(i,j) = ( vTrans(i,j)+vAbs )* 0.5 _d 0 * tracer(i,j-1)
82         &          + ( vTrans(i,j)-vAbs )* 0.5 _d 0 * tracer(i,j)
83    
84         ENDDO         ENDDO
85        ENDDO        ENDDO
86    
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.22