pkg/generic_advdiff/gad_fluxlimit_adv_r.F

C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_fluxlimit_adv_r.F,v 1.11 2006/06/19 14:40:43 jmc Exp $
C $Name:  $

#include "GAD_OPTIONS.h"

CBOP
C !ROUTINE: GAD_FLUXLIMIT_ADV_R

C !INTERFACE: ==========================================================
      SUBROUTINE GAD_FLUXLIMIT_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 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,bj              :: tile indices
C  k                  :: vertical level
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  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
C  wLoc               :: velocity, vertical component
      INTEGER i,j,kp1,km1,km2
      _RL Cr,Rjm,Rj,Rjp
      _RL wLoc, wCFL
C Statement function provides Limiter(Cr)
#include "GAD_FLUX_LIMITER.h"
CEOP

      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

         wLoc = wFld(i,j)
         wCFL = ABS( wLoc*dTarg*recip_drC(k) )
         Rjp=(tracer(i,j,kp1)-tracer(i,j,k))
     &        *maskC(i,j,kp1,bi,bj)
         Rj= (tracer(i,j,k)  -tracer(i,j,kM1))
         Rjm=(tracer(i,j,km1)-tracer(i,j,kM2))
     &        *maskC(i,j,km2,bi,bj)

         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,bj)*(
     &     rTrans(i,j)*
     &        (tracer(i,j,k)+tracer(i,j,kM1))*0.5 _d 0
     &    +ABS(rTrans(i,j))*((1.-Cr)+wCFL*Cr)
     &                     *Rj*0.5 _d 0 )
        ENDDO
       ENDDO
      ENDIF

      RETURN
      END
1	jmc	1.12	C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_fluxlimit_adv_r.F,v 1.11 2006/06/19 14:40:43 jmc 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	jmc	1.10	SUBROUTINE GAD_FLUXLIMIT_ADV_R(
11	jmc	1.9	I bi,bj,k,dTarg,
12	jmc	1.10	I rTrans, wFld,
13	adcroft	1.1	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	jmc	1.9	C bi,bj :: tile indices
39	edhill	1.8	C k :: vertical level
40			C rTrans :: vertical volume transport
41	jmc	1.10	C wFld :: vertical flow
42	edhill	1.8	C tracer :: tracer field
43			C myThid :: thread number
44	jmc	1.9	INTEGER bi,bj,k
45	adcroft	1.3	_RL dTarg
46	adcroft	1.1	_RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
47	jmc	1.10	_RL wFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
48	jmc	1.9	_RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
49	adcroft	1.5	INTEGER myThid
50
51			C !OUTPUT PARAMETERS: ==================================================
52	edhill	1.8	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	edhill	1.8	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	jmc	1.10	C wLoc :: velocity, vertical component
64	jmc	1.9	INTEGER i,j,kp1,km1,km2
65	adcroft	1.1	_RL Cr,Rjm,Rj,Rjp
66	jmc	1.12	_RL wLoc, wCFL
67	adcroft	1.5	C Statement function provides Limiter(Cr)
68	adcroft	1.1	#include "GAD_FLUX_LIMITER.h"
69	adcroft	1.5	CEOP
70	adcroft	1.1
71			km2=MAX(1,k-2)
72			km1=MAX(1,k-1)
73			kp1=MIN(Nr,k+1)
74
75			IF ( k.GT.Nr) THEN
76			DO j=1-Oly,sNy+Oly
77			DO i=1-Olx,sNx+Olx
78			wT(i,j) = 0.
79			ENDDO
80			ENDDO
81			ELSE
82			DO j=1-Oly,sNy+Oly
83			DO i=1-Olx,sNx+Olx
84	jmc	1.7
85	jmc	1.11	wLoc = wFld(i,j)
86	jmc	1.12	wCFL = ABS( wLocdTargrecip_drC(k) )
87	jmc	1.9	Rjp=(tracer(i,j,kp1)-tracer(i,j,k))
88			& *maskC(i,j,kp1,bi,bj)
89			Rj= (tracer(i,j,k) -tracer(i,j,kM1))
90			Rjm=(tracer(i,j,km1)-tracer(i,j,kM2))
91			& *maskC(i,j,km2,bi,bj)
92	jmc	1.7
93	adcroft	1.1	IF (Rj.NE.0.) THEN
94	adcroft	1.4	IF (rTrans(i,j).LT.0.) THEN
95	adcroft	1.1	Cr=Rjm/Rj
96			ELSE
97			Cr=Rjp/Rj
98			ENDIF
99			ELSE
100	adcroft	1.4	IF (rTrans(i,j).LT.0.) THEN
101	adcroft	1.1	Cr=Rjm*1.E20
102			ELSE
103			Cr=Rjp*1.E20
104			ENDIF
105			ENDIF
106			Cr=Limiter(Cr)
107	jmc	1.9	wT(i,j) = maskC(i,j,kM1,bi,bj)*(
108	jmc	1.2	& rTrans(i,j)*
109	jmc	1.9	& (tracer(i,j,k)+tracer(i,j,kM1))*0.5 _d 0
110	jmc	1.12	& +ABS(rTrans(i,j))((1.-Cr)+wCFLCr)
111			& Rj0.5 _d 0 )
112	adcroft	1.1	ENDDO
113			ENDDO
114			ENDIF
115
116			RETURN
117			END