pkg/generic_advdiff/gad_dst2u1_adv_x.F

C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_dst2u1_adv_x.F,v 1.1 2005/10/22 19:56:33 jmc Exp $
C $Name:  $

#include "GAD_OPTIONS.h"

CBOP
C !ROUTINE: GAD_DST2U1_ADV_X

C !INTERFACE: ==========================================================
      SUBROUTINE GAD_DST2U1_ADV_X( 
     I           bi,bj,k, advectionScheme,
     I           deltaTloc, uTrans, uVel,
     I           tracer,
     O           uT,
     I           myThid )

C !DESCRIPTION:
C  Calculates the area integrated zonal flux due to advection
C  of a tracer using second-order Direct Space and Time (DST-2) 
C  interpolation (=Lax-Wendroff) or simple 1rst order upwind scheme.

C !USES: ===============================================================
      IMPLICIT NONE
#include "SIZE.h"
#include "GRID.h"
#include "GAD.h"

C !INPUT PARAMETERS: ===================================================
C  bi,bj             :: tile indices
C  k                 :: vertical level
C  advectionScheme   :: advection scheme to use: either 2nd Order DST 
C                                                or 1rst Order Upwind
C  uTrans            :: zonal volume transport
C  uVel              :: zonal flow
C  tracer            :: tracer field
C  myThid            :: thread number
      INTEGER bi,bj,k
      INTEGER advectionScheme
      _RL deltaTloc
      _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL uVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
      _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      INTEGER myThid

C !OUTPUT PARAMETERS: ==================================================
C  uT                :: zonal advective flux
      _RL uT    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)

C !LOCAL VARIABLES: ====================================================
C  i,j               :: loop indices
C  rLimit            :: centered (vs upwind) fraction
C  uFld              :: velocity [m/s], zonal component
C  uCFL              :: Courant-Friedrich-Levy number
      INTEGER i,j
      _RL uFld, uCFL, xLimit
CEOP

      xLimit = 0. _d 0
      IF ( advectionScheme.EQ.ENUM_DST2 ) xLimit = 1. _d 0

      DO j=1-Oly,sNy+Oly
       uT(1-Olx,j)=0.
       DO i=1-Olx+1,sNx+Olx

c       uFld = uVel(i,j,k,bi,bj)
        uFld = uTrans(i,j)*recip_dyG(i,j,bi,bj)
     &        *recip_drF(k)*_recip_hFacW(i,j,k,bi,bj)
        uCFL = ABS(uFld*deltaTloc*recip_dxC(i,j,bi,bj))

        uT(i,j) = 
     &     uTrans(i,j)*(Tracer(i-1,j)+Tracer(i,j))*0.5 _d 0
     &   + ( 1. _d 0 - xLimit*(1. _d 0 - uCFL) )*ABS(uTrans(i,j))
     &                *(tracer(i-1,j)-tracer(i,j))*0.5 _d 0
       ENDDO
      ENDDO

      RETURN
      END
1	C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_dst2u1_adv_x.F,v 1.1 2005/10/22 19:56:33 jmc Exp $
2	C $Name: $
3
4	#include "GAD_OPTIONS.h"
5
6	CBOP
7	C !ROUTINE: GAD_DST2U1_ADV_X
8
9	C !INTERFACE: ==========================================================
10	SUBROUTINE GAD_DST2U1_ADV_X(
11	I bi,bj,k, advectionScheme,
12	I deltaTloc, uTrans, uVel,
13	I tracer,
14	O uT,
15	I myThid )
16
17	C !DESCRIPTION:
18	C Calculates the area integrated zonal flux due to advection
19	C of a tracer using second-order Direct Space and Time (DST-2)
20	C interpolation (=Lax-Wendroff) or simple 1rst order upwind scheme.
21
22	C !USES: ===============================================================
23	IMPLICIT NONE
24	#include "SIZE.h"
25	#include "GRID.h"
26	#include "GAD.h"
27
28	C !INPUT PARAMETERS: ===================================================
29	C bi,bj :: tile indices
30	C k :: vertical level
31	C advectionScheme :: advection scheme to use: either 2nd Order DST
32	C or 1rst Order Upwind
33	C uTrans :: zonal volume transport
34	C uVel :: zonal flow
35	C tracer :: tracer field
36	C myThid :: thread number
37	INTEGER bi,bj,k
38	INTEGER advectionScheme
39	_RL deltaTloc
40	_RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
41	_RL uVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
42	_RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
43	INTEGER myThid
44
45	C !OUTPUT PARAMETERS: ==================================================
46	C uT :: zonal advective flux
47	_RL uT (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
48
49	C !LOCAL VARIABLES: ====================================================
50	C i,j :: loop indices
51	C rLimit :: centered (vs upwind) fraction
52	C uFld :: velocity [m/s], zonal component
53	C uCFL :: Courant-Friedrich-Levy number
54	INTEGER i,j
55	_RL uFld, uCFL, xLimit
56	CEOP
57
58	xLimit = 0. _d 0
59	IF ( advectionScheme.EQ.ENUM_DST2 ) xLimit = 1. _d 0
60
61	DO j=1-Oly,sNy+Oly
62	uT(1-Olx,j)=0.
63	DO i=1-Olx+1,sNx+Olx
64
65	c uFld = uVel(i,j,k,bi,bj)
66	uFld = uTrans(i,j)*recip_dyG(i,j,bi,bj)
67	& recip_drF(k)_recip_hFacW(i,j,k,bi,bj)
68	uCFL = ABS(uFlddeltaTlocrecip_dxC(i,j,bi,bj))
69
70	uT(i,j) =
71	& uTrans(i,j)(Tracer(i-1,j)+Tracer(i,j))0.5 _d 0
72	& + ( 1. _d 0 - xLimit(1. _d 0 - uCFL) )ABS(uTrans(i,j))
73	& (tracer(i-1,j)-tracer(i,j))0.5 _d 0
74	ENDDO
75	ENDDO
76
77	RETURN
78	END