pkg/generic_advdiff/gad_dst3_adv_r.F

C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_dst3_adv_r.F,v 1.9 2006/10/22 01:08:04 jmc Exp $
C $Name:  $

#include "GAD_OPTIONS.h"

CBOP
C !ROUTINE: GAD_DST3_ADV_R

C !INTERFACE: ==========================================================
      SUBROUTINE GAD_DST3_ADV_R(
     I           bi,bj,k,dTarg,
     I           rTrans, wFld,
     I           tracer,
     O           wT,
     I           myThid )

C !DESCRIPTION:
C  Calculates the area integrated vertical flux due to advection of a tracer
C  using 3rd-order Direct Space and Time (DST-3) Advection Scheme

C !USES: ===============================================================
      IMPLICIT NONE

C     == GLobal variables ==
#include "SIZE.h"
#ifdef OLD_DST3_FORMULATION
#include "EEPARAMS.h"
#include "PARAMS.h"
#endif
#include "GRID.h"
#include "GAD.h"

C     == Routine arguments ==
C !INPUT PARAMETERS: ===================================================
C  bi,bj             :: tile indices
C  k                 :: vertical level
C  deltaTloc         :: local time-step (s)
C  rTrans            :: vertical volume transport
C  wFld              :: vertical flow
C  tracer            :: tracer field
C  myThid            :: thread number
      INTEGER bi,bj,k
      _RL dTarg
      _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL wFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      INTEGER myThid

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

C     == Local variables ==
C !LOCAL VARIABLES: ====================================================
C  i,j               :: loop indices
C  km1               :: =max( k-1 , 1 )
C  wLoc              :: velocity, vertical component
C  wCFL              :: Courant-Friedrich-Levy number
      INTEGER i,j,kp1,km1,km2
      _RL wLoc
      _RL Rjm,Rj,Rjp,cfl,d0,d1
#ifdef OLD_DST3_FORMULATION
      _RL psiP,psiM,thetaP,thetaM
      _RL smallNo

      IF (inAdMode) THEN
       smallNo = 1.0D-20
      ELSE
       smallNo = 1.0D-20
      ENDIF
#endif

      km2=MAX(1,k-2)
      km1=MAX(1,k-1)
      kp1=MIN(Nr,k+1)

      DO j=1-Oly,sNy+Oly
       DO i=1-Olx,sNx+Olx
        Rjp=(tracer(i,j,k)-tracer(i,j,kp1))
     &         *maskC(i,j,kp1,bi,bj)
        Rj =(tracer(i,j,km1)-tracer(i,j,k))
     &         *maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
        Rjm=(tracer(i,j,km2)-tracer(i,j,km1))
     &         *maskC(i,j,km1,bi,bj)

        wLoc = wFld(i,j)
c       wLoc = rTrans(i,j)*recip_rA(i,j,bi,bj)
        cfl=ABS(wLoc*dTarg*recip_drC(k))
        d0=(2.-cfl)*(1.-cfl)*oneSixth
        d1=(1.-cfl*cfl)*oneSixth
#ifdef OLD_DST3_FORMULATION
        IF ( ABS(Rj).LT.smallNo .OR.
     &       ABS(Rjm).LT.smallNo ) THEN
         thetaP=0.
         psiP=0.
        ELSE
         thetaP=(Rjm+smallNo)/(smallNo+Rj)
         psiP=d0+d1*thetaP
        ENDIF
        IF ( ABS(Rj).LT.smallNo .OR.
     &       ABS(Rjp).LT.smallNo ) THEN
         thetaM=0.
         psiM=0.
        ELSE
         thetaM=(Rjp+smallNo)/(smallNo+Rj)
         psiM=d0+d1*thetaM
        ENDIF
         wT(i,j)=
     &    0.5*(rTrans(i,j)+ABS(rTrans(i,j)))
     &       *( tracer(i,j, k ) + psiM*Rj )
     &   +0.5*(rTrans(i,j)-ABS(rTrans(i,j)))
     &       *( tracer(i,j,km1) - psiP*Rj )
#else /* OLD_DST3_FORMULATION */
        wT(i,j)=
     &    0.5*(rTrans(i,j)+ABS(rTrans(i,j)))
     &       *( tracer(i,j, k ) + (d0*Rj+d1*Rjp) )
     &   +0.5*(rTrans(i,j)-ABS(rTrans(i,j)))
     &       *( tracer(i,j,km1) - (d0*Rj+d1*Rjm) )
#endif /* OLD_DST3_FORMULATION */

       ENDDO
      ENDDO

      RETURN
      END
1	jmc	1.10	C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_dst3_adv_r.F,v 1.9 2006/10/22 01:08:04 jmc Exp $
2	adcroft	1.2	C $Name: $
3	adcroft	1.1
4			#include "GAD_OPTIONS.h"
5
6	jmc	1.6	CBOP
7			C !ROUTINE: GAD_DST3_ADV_R
8
9			C !INTERFACE: ==========================================================
10	jmc	1.7	SUBROUTINE GAD_DST3_ADV_R(
11	jmc	1.6	I bi,bj,k,dTarg,
12	jmc	1.7	I rTrans, wFld,
13	adcroft	1.1	I tracer,
14			O wT,
15			I myThid )
16	jmc	1.6
17			C !DESCRIPTION:
18			C Calculates the area integrated vertical flux due to advection of a tracer
19			C using 3rd-order Direct Space and Time (DST-3) Advection Scheme
20
21			C !USES: ===============================================================
22	adcroft	1.1	IMPLICIT NONE
23
24			C == GLobal variables ==
25			#include "SIZE.h"
26	jmc	1.10	#ifdef OLD_DST3_FORMULATION
27	adcroft	1.1	#include "EEPARAMS.h"
28			#include "PARAMS.h"
29	jmc	1.10	#endif
30			#include "GRID.h"
31	adcroft	1.1	#include "GAD.h"
32
33			C == Routine arguments ==
34	jmc	1.6	C !INPUT PARAMETERS: ===================================================
35			C bi,bj :: tile indices
36			C k :: vertical level
37			C deltaTloc :: local time-step (s)
38			C rTrans :: vertical volume transport
39	jmc	1.7	C wFld :: vertical flow
40	jmc	1.6	C tracer :: tracer field
41			C myThid :: thread number
42			INTEGER bi,bj,k
43	adcroft	1.1	_RL dTarg
44			_RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
45	jmc	1.7	_RL wFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
46	jmc	1.6	_RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
47			INTEGER myThid
48
49			C !OUTPUT PARAMETERS: ==================================================
50			C wT :: vertical advective flux
51	adcroft	1.1	_RL wT (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
52
53			C == Local variables ==
54	jmc	1.6	C !LOCAL VARIABLES: ====================================================
55			C i,j :: loop indices
56			C km1 :: =max( k-1 , 1 )
57	jmc	1.7	C wLoc :: velocity, vertical component
58	jmc	1.6	C wCFL :: Courant-Friedrich-Levy number
59			INTEGER i,j,kp1,km1,km2
60	jmc	1.9	_RL wLoc
61	adcroft	1.1	_RL Rjm,Rj,Rjp,cfl,d0,d1
62	jmc	1.9	#ifdef OLD_DST3_FORMULATION
63	adcroft	1.1	_RL psiP,psiM,thetaP,thetaM
64	heimbach	1.4	_RL smallNo
65	adcroft	1.1
66	heimbach	1.4	IF (inAdMode) THEN
67			smallNo = 1.0D-20
68			ELSE
69			smallNo = 1.0D-20
70			ENDIF
71	jmc	1.9	#endif
72	heimbach	1.4
73	adcroft	1.1	km2=MAX(1,k-2)
74			km1=MAX(1,k-1)
75			kp1=MIN(Nr,k+1)
76
77			DO j=1-Oly,sNy+Oly
78			DO i=1-Olx,sNx+Olx
79	jmc	1.6	Rjp=(tracer(i,j,k)-tracer(i,j,kp1))
80			& *maskC(i,j,kp1,bi,bj)
81			Rj =(tracer(i,j,km1)-tracer(i,j,k))
82			& maskC(i,j,k,bi,bj)maskC(i,j,km1,bi,bj)
83			Rjm=(tracer(i,j,km2)-tracer(i,j,km1))
84			& *maskC(i,j,km1,bi,bj)
85	adcroft	1.1
86	jmc	1.8	wLoc = wFld(i,j)
87			c wLoc = rTrans(i,j)*recip_rA(i,j,bi,bj)
88	jmc	1.9	cfl=ABS(wLocdTargrecip_drC(k))
89	adcroft	1.1	d0=(2.-cfl)(1.-cfl)oneSixth
90			d1=(1.-cflcfl)oneSixth
91	jmc	1.9	#ifdef OLD_DST3_FORMULATION
92	heimbach	1.4	IF ( ABS(Rj).LT.smallNo .OR.
93			& ABS(Rjm).LT.smallNo ) THEN
94			thetaP=0.
95			psiP=0.
96			ELSE
97			thetaP=(Rjm+smallNo)/(smallNo+Rj)
98			psiP=d0+d1*thetaP
99			ENDIF
100			IF ( ABS(Rj).LT.smallNo .OR.
101			& ABS(Rjp).LT.smallNo ) THEN
102			thetaM=0.
103			psiM=0.
104			ELSE
105			thetaM=(Rjp+smallNo)/(smallNo+Rj)
106			psiM=d0+d1*thetaM
107			ENDIF
108			wT(i,j)=
109	jmc	1.9	& 0.5*(rTrans(i,j)+ABS(rTrans(i,j)))
110	jmc	1.6	& ( tracer(i,j, k ) + psiMRj )
111	jmc	1.9	& +0.5*(rTrans(i,j)-ABS(rTrans(i,j)))
112	jmc	1.6	& ( tracer(i,j,km1) - psiPRj )
113	jmc	1.9	#else /* OLD_DST3_FORMULATION */
114			wT(i,j)=
115			& 0.5*(rTrans(i,j)+ABS(rTrans(i,j)))
116			& ( tracer(i,j, k ) + (d0Rj+d1*Rjp) )
117			& +0.5*(rTrans(i,j)-ABS(rTrans(i,j)))
118			& ( tracer(i,j,km1) - (d0Rj+d1*Rjm) )
119			#endif /* OLD_DST3_FORMULATION */
120	jmc	1.5
121	adcroft	1.1	ENDDO
122			ENDDO
123
124			RETURN
125			END