pkg/generic_advdiff/gad_fluxlimit_adv_r.F

C $Header: /u/gcmpack/models/MITgcmUV/pkg/generic_advdiff/gad_fluxlimit_adv_r.F,v 1.5 2001/09/21 13:11:43 adcroft Exp $
C $Name:  $

#include "GAD_OPTIONS.h"

CBOP
C !ROUTINE: GAD_FLUXLIMIT_ADV_R

C !INTERFACE: ==========================================================
      SUBROUTINE GAD_FLUXLIMIT_ADV_R( 
     I           bi_arg,bj_arg,k,dTarg,
     I           rTrans, wVel,
     I           tracer,
     O           wT,
     I           myThid )

C !DESCRIPTION:
C Calculates the area integrated vertical flux due to advection of a tracer
C using second-order interpolation with a flux limiter:
C \begin{equation*}
C F^x_{adv} = W \overline{ \theta }^k
C - \frac{1}{2} \left(
C     [ 1 - \psi(C_r) ] |W|
C    + W \frac{w \Delta t}{\Delta r_c} \psi(C_r)
C              \right) \delta_k \theta
C \end{equation*}
C where the $\psi(C_r)$ is the limiter function and $C_r$ is
C the slope ratio.

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

C !INPUT PARAMETERS: ===================================================
C  bi_arg,bj_arg        :: tile indices
C  k                    :: vertical level
C  rTrans               :: vertical volume transport
C  wVel                 :: vertical flow
C  tracer               :: tracer field
C  myThid               :: thread number
      INTEGER bi_arg,bj_arg,k
      _RL dTarg
      _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL wVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
      _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
      INTEGER myThid

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

C !LOCAL VARIABLES: ====================================================
C  i,j                  :: loop indices
C  kp1                  :: =min( k+1 , Nr )
C  km1                  :: =max( k-1 , 1 )
C  km2                  :: =max( k-2 , 1 )
C  bi,bj                :: tile indices or (1,1) depending on use
C  Cr                   :: slope ratio
C  Rjm,Rj,Rjp           :: differences at i-1,i,i+1
      INTEGER i,j,kp1,km1,km2,bi,bj
      _RL Cr,Rjm,Rj,Rjp
C Statement function provides Limiter(Cr)
#include "GAD_FLUX_LIMITER.h"
CEOP

      IF (.NOT. multiDimAdvection) THEN
C      If using the standard time-stepping/advection schemes (ie. AB-II)
C      then the data-structures are all global arrays
       bi=bi_arg
       bj=bj_arg
      ELSE
C      otherwise if using the multi-dimensional advection schemes
C      then the data-structures are all local arrays except
C      for maskC(...) and wVel(...)
       bi=1
       bj=1
      ENDIF

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

      IF ( k.GT.Nr) THEN
       DO j=1-Oly,sNy+Oly
        DO i=1-Olx,sNx+Olx
         wT(i,j) = 0.
        ENDDO
       ENDDO
      ELSE
       DO j=1-Oly,sNy+Oly
        DO i=1-Olx,sNx+Olx
         Rjp=(tracer(i,j,kp1,bi,bj)-tracer(i,j,k,bi,bj))
     &        *maskC(i,j,kp1,bi_arg,bj_arg)
         Rj=(tracer(i,j,k,bi,bj)-tracer(i,j,kM1,bi,bj))
         Rjm=(tracer(i,j,km1,bi,bj)-tracer(i,j,kM2,bi,bj))
     &        *maskC(i,j,km2,bi_arg,bj_arg)
         IF (Rj.NE.0.) THEN
          IF (rTrans(i,j).LT.0.) THEN
            Cr=Rjm/Rj
          ELSE
            Cr=Rjp/Rj
          ENDIF
         ELSE
          IF (rTrans(i,j).LT.0.) THEN
            Cr=Rjm*1.E20
          ELSE
            Cr=Rjp*1.E20
          ENDIF
         ENDIF
         Cr=Limiter(Cr)
         wT(i,j) = maskC(i,j,kM1,bi_arg,bj_arg)*(
     &     rTrans(i,j)*
     &        (Tracer(i,j,k,bi,bj)+Tracer(i,j,kM1,bi,bj))*0.5 _d 0
     &    +(ABS(rTrans(i,j))*(1-Cr)
     &      +rTrans(i,j)*wVel(i,j,k,bi_arg,bj_arg)*dTarg*recip_drC(k)
     &                  *Cr
     &     )*Rj*0.5 _d 0                )
        ENDDO
       ENDDO
      ENDIF

      RETURN
      END
1	adcroft	1.6	C $Header: /u/gcmpack/models/MITgcmUV/pkg/generic_advdiff/gad_fluxlimit_adv_r.F,v 1.5 2001/09/21 13:11:43 adcroft Exp $
2	jmc	1.2	C $Name: $
3	adcroft	1.1
4			#include "GAD_OPTIONS.h"
5
6	adcroft	1.5	CBOP
7			C !ROUTINE: GAD_FLUXLIMIT_ADV_R
8
9			C !INTERFACE: ==========================================================
10	adcroft	1.1	SUBROUTINE GAD_FLUXLIMIT_ADV_R(
11	adcroft	1.3	I bi_arg,bj_arg,k,dTarg,
12	adcroft	1.1	I rTrans, wVel,
13			I tracer,
14			O wT,
15			I myThid )
16	adcroft	1.5
17			C !DESCRIPTION:
18			C Calculates the area integrated vertical flux due to advection of a tracer
19			C using second-order interpolation with a flux limiter:
20			C \begin{equation*}
21			C F^x_{adv} = W \overline{ \theta }^k
22			C - \frac{1}{2} \left(
23			C [ 1 - \psi(C_r) ] \|W\|
24			C + W \frac{w \Delta t}{\Delta r_c} \psi(C_r)
25			C \right) \delta_k \theta
26			C \end{equation*}
27			C where the $\psi(C_r)$ is the limiter function and $C_r$ is
28			C the slope ratio.
29
30			C !USES: ===============================================================
31	adcroft	1.1	IMPLICIT NONE
32			#include "SIZE.h"
33			#include "GRID.h"
34	adcroft	1.3	#include "EEPARAMS.h"
35			#include "PARAMS.h"
36	adcroft	1.1
37	adcroft	1.5	C !INPUT PARAMETERS: ===================================================
38			C bi_arg,bj_arg :: tile indices
39			C k :: vertical level
40			C rTrans :: vertical volume transport
41			C wVel :: vertical flow
42			C tracer :: tracer field
43			C myThid :: thread number
44	adcroft	1.3	INTEGER bi_arg,bj_arg,k
45			_RL dTarg
46	adcroft	1.1	_RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
47			_RL wVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
48			_RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
49	adcroft	1.5	INTEGER myThid
50
51			C !OUTPUT PARAMETERS: ==================================================
52			C wT :: vertical advective flux
53	adcroft	1.1	_RL wT (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
54
55	adcroft	1.5	C !LOCAL VARIABLES: ====================================================
56			C i,j :: loop indices
57			C kp1 :: =min( k+1 , Nr )
58			C km1 :: =max( k-1 , 1 )
59			C km2 :: =max( k-2 , 1 )
60			C bi,bj :: tile indices or (1,1) depending on use
61			C Cr :: slope ratio
62			C Rjm,Rj,Rjp :: differences at i-1,i,i+1
63	adcroft	1.3	INTEGER i,j,kp1,km1,km2,bi,bj
64	adcroft	1.1	_RL Cr,Rjm,Rj,Rjp
65	adcroft	1.5	C Statement function provides Limiter(Cr)
66	adcroft	1.1	#include "GAD_FLUX_LIMITER.h"
67	adcroft	1.5	CEOP
68	adcroft	1.1
69	adcroft	1.3	IF (.NOT. multiDimAdvection) THEN
70			C If using the standard time-stepping/advection schemes (ie. AB-II)
71			C then the data-structures are all global arrays
72			bi=bi_arg
73			bj=bj_arg
74			ELSE
75			C otherwise if using the multi-dimensional advection schemes
76			C then the data-structures are all local arrays except
77			C for maskC(...) and wVel(...)
78			bi=1
79			bj=1
80			ENDIF
81
82	adcroft	1.1	km2=MAX(1,k-2)
83			km1=MAX(1,k-1)
84			kp1=MIN(Nr,k+1)
85
86			IF ( k.GT.Nr) THEN
87			DO j=1-Oly,sNy+Oly
88			DO i=1-Olx,sNx+Olx
89			wT(i,j) = 0.
90			ENDDO
91			ENDDO
92			ELSE
93			DO j=1-Oly,sNy+Oly
94			DO i=1-Olx,sNx+Olx
95			Rjp=(tracer(i,j,kp1,bi,bj)-tracer(i,j,k,bi,bj))
96	adcroft	1.6	& *maskC(i,j,kp1,bi_arg,bj_arg)
97	adcroft	1.1	Rj=(tracer(i,j,k,bi,bj)-tracer(i,j,kM1,bi,bj))
98			Rjm=(tracer(i,j,km1,bi,bj)-tracer(i,j,kM2,bi,bj))
99	adcroft	1.6	& *maskC(i,j,km2,bi_arg,bj_arg)
100	adcroft	1.1	IF (Rj.NE.0.) THEN
101	adcroft	1.4	IF (rTrans(i,j).LT.0.) THEN
102	adcroft	1.1	Cr=Rjm/Rj
103			ELSE
104			Cr=Rjp/Rj
105			ENDIF
106			ELSE
107	adcroft	1.4	IF (rTrans(i,j).LT.0.) THEN
108	adcroft	1.1	Cr=Rjm*1.E20
109			ELSE
110			Cr=Rjp*1.E20
111			ENDIF
112			ENDIF
113			Cr=Limiter(Cr)
114	adcroft	1.3	wT(i,j) = maskC(i,j,kM1,bi_arg,bj_arg)*(
115	jmc	1.2	& rTrans(i,j)*
116	adcroft	1.1	& (Tracer(i,j,k,bi,bj)+Tracer(i,j,kM1,bi,bj))*0.5 _d 0
117	jmc	1.2	& +(ABS(rTrans(i,j))*(1-Cr)
118	adcroft	1.3	& +rTrans(i,j)wVel(i,j,k,bi_arg,bj_arg)dTarg*recip_drC(k)
119	jmc	1.2	& *Cr
120			& )Rj0.5 _d 0 )
121	adcroft	1.1	ENDDO
122			ENDDO
123			ENDIF
124
125			RETURN
126			END