pkg/generic_advdiff/gad_ppm_flx_r.F

C $Header:  $
C $Name:  $

#     include "GAD_OPTIONS.h"

      SUBROUTINE GAD_PPM_FLX_R(bi,bj,ix,iy,
     &           delT,wvel,wfac,fhat,
     &           flux,myThid )
C     |================================================================|
C     | PPM_FLX_R: evaluate PPM flux on grid-cell edges.               |
C     |================================================================|

          implicit none

C     =============================================== global variables
#         include "SIZE.h"
#         include "GRID.h"
#         include "GAD.h"

C     ================================================================
C       bi,bj     :: tile indexing.
C       ix        :: x-index.
C       iy        :: y-index.
C       delT      :: time-step.
C       wvel      :: vel.-comp in r-direction.
C       wfac      :: vel.-flux in r-direction.
C       fhat      :: row of poly. coeff.
C       flux      :: adv.-flux in r-direction.
C       myThid    :: thread number.
C     ================================================================
          integer bi,bj,ix,iy
          _RL delT(1:Nr)
          _RL wvel(1-OLx:sNx+OLx,
     &             1-OLy:sNy+OLy,
     &             1:Nr)
          _RL wfac(1-OLx:sNx+OLx,
     &             1-OLy:sNy+OLy,
     &             1:Nr)
          _RL fhat(1:3 ,
     &             1:Nr)
          _RL flux(1-OLx:sNx+OLx,
     &             1-OLy:sNy+OLy,
     &             1:Nr)
          integer myThid

C     ================================================================
C       ir        :: r-indexing.
C       wCFL      :: CFL number.
C       intF      :: upwind tracer edge-value.
C       ss11,ss22 :: int. endpoints.
C       ivec      :: int. basis vec.
C     ================================================================
          integer ir
          _RL wCFL,intF
          _RL ss11,ss22
          _RL ivec(1:3)

C     ================================================================
C       (1): calc. "departure-points" for each grid-cell edge by int-
C            egrating edge position backward in time over one single
C            time-step. This is a "single-cell" implementation: requ-
C            ires CFL <= 1.0.
C       (2): calc. flux as the integral of the upwind grid-cell poly-
C            nomial over the deformation interval found in (1).
C     ================================================================

          do  ir = +2, Nr

              if (wvel(ix,iy,ir) .eq. 0. _d 0) then

                  flux(ix,iy,ir)    = 0. _d 0

              else

              if (wvel(ix,iy,ir) .lt. 0. _d 0) then

C     ==================== integrate PQM profile over upwind cell IR-1
              wCFL = wvel(ix,iy,ir)
     &             * delT(ir-1)*recip_drF(ir-1)

              ss11 = +1. _d 0 + 2. _d 0 * wCFL
              ss22 = +1. _d 0

C     ==================== integrate profile over region swept by face
              ivec(1) = ss22 - ss11
              ivec(2) =(ss22 ** 2
     &                - ss11 ** 2)*(1. _d 0 / 2. _d 0)
              ivec(3) =(ss22 ** 3
     &                - ss11 ** 3)*(1. _d 0 / 3. _d 0)

              intF = ivec(1) * fhat(1,ir-1)
     &             + ivec(2) * fhat(2,ir-1)
     &             + ivec(3) * fhat(3,ir-1)

              intF = intF / (ss22 - ss11)

              else

C     ==================== integrate PQM profile over upwind cell IR+0
              wCFL = wvel(ix,iy,ir)
     &             * delT(ir-0)*recip_drF(ir-0)

              ss11 = -1. _d 0 + 2. _d 0 * wCFL
              ss22 = -1. _d 0

C     ==================== integrate profile over region swept by face
              ivec(1) = ss22 - ss11
              ivec(2) =(ss22 ** 2
     &                - ss11 ** 2)*(1. _d 0 / 2. _d 0)
              ivec(3) =(ss22 ** 3
     &                - ss11 ** 3)*(1. _d 0 / 3. _d 0)

              intF = ivec(1) * fhat(1,ir-0)
     &             + ivec(2) * fhat(2,ir-0)
     &             + ivec(3) * fhat(3,ir-0)

              intF = intF / (ss22 - ss11)

              end if

C     ==================== calc. flux = upwind tracer * face-transport
              flux(ix,iy,ir) = + wfac(ix,iy,ir) * intF

              end if

          end do

          return

c     end subroutine GAD_PPM_FLX_R
      end
1	jmc	1.1	C $Header: $
2			C $Name: $
3
4			# include "GAD_OPTIONS.h"
5
6			SUBROUTINE GAD_PPM_FLX_R(bi,bj,ix,iy,
7			& delT,wvel,wfac,fhat,
8			& flux,myThid )
9			C \|================================================================\|
10			C \| PPM_FLX_R: evaluate PPM flux on grid-cell edges. \|
11			C \|================================================================\|
12
13			implicit none
14
15			C =============================================== global variables
16			# include "SIZE.h"
17			# include "GRID.h"
18			# include "GAD.h"
19
20			C ================================================================
21			C bi,bj :: tile indexing.
22			C ix :: x-index.
23			C iy :: y-index.
24			C delT :: time-step.
25			C wvel :: vel.-comp in r-direction.
26			C wfac :: vel.-flux in r-direction.
27			C fhat :: row of poly. coeff.
28			C flux :: adv.-flux in r-direction.
29			C myThid :: thread number.
30			C ================================================================
31			integer bi,bj,ix,iy
32			_RL delT(1:Nr)
33			_RL wvel(1-OLx:sNx+OLx,
34			& 1-OLy:sNy+OLy,
35			& 1:Nr)
36			_RL wfac(1-OLx:sNx+OLx,
37			& 1-OLy:sNy+OLy,
38			& 1:Nr)
39			_RL fhat(1:3 ,
40			& 1:Nr)
41			_RL flux(1-OLx:sNx+OLx,
42			& 1-OLy:sNy+OLy,
43			& 1:Nr)
44			integer myThid
45
46			C ================================================================
47			C ir :: r-indexing.
48			C wCFL :: CFL number.
49			C intF :: upwind tracer edge-value.
50			C ss11,ss22 :: int. endpoints.
51			C ivec :: int. basis vec.
52			C ================================================================
53			integer ir
54			_RL wCFL,intF
55			_RL ss11,ss22
56			_RL ivec(1:3)
57
58			C ================================================================
59			C (1): calc. "departure-points" for each grid-cell edge by int-
60			C egrating edge position backward in time over one single
61			C time-step. This is a "single-cell" implementation: requ-
62			C ires CFL <= 1.0.
63			C (2): calc. flux as the integral of the upwind grid-cell poly-
64			C nomial over the deformation interval found in (1).
65			C ================================================================
66
67			do ir = +2, Nr
68
69			if (wvel(ix,iy,ir) .eq. 0. _d 0) then
70
71			flux(ix,iy,ir) = 0. _d 0
72
73			else
74
75			if (wvel(ix,iy,ir) .lt. 0. _d 0) then
76
77			C ==================== integrate PQM profile over upwind cell IR-1
78			wCFL = wvel(ix,iy,ir)
79			& * delT(ir-1)*recip_drF(ir-1)
80
81			ss11 = +1. _d 0 + 2. _d 0 * wCFL
82			ss22 = +1. _d 0
83
84			C ==================== integrate profile over region swept by face
85			ivec(1) = ss22 - ss11
86			ivec(2) =(ss22 ** 2
87			& - ss11 ** 2)*(1. _d 0 / 2. _d 0)
88			ivec(3) =(ss22 ** 3
89			& - ss11 ** 3)*(1. _d 0 / 3. _d 0)
90
91			intF = ivec(1) * fhat(1,ir-1)
92			& + ivec(2) * fhat(2,ir-1)
93			& + ivec(3) * fhat(3,ir-1)
94
95			intF = intF / (ss22 - ss11)
96
97			else
98
99			C ==================== integrate PQM profile over upwind cell IR+0
100			wCFL = wvel(ix,iy,ir)
101			& * delT(ir-0)*recip_drF(ir-0)
102
103			ss11 = -1. _d 0 + 2. _d 0 * wCFL
104			ss22 = -1. _d 0
105
106			C ==================== integrate profile over region swept by face
107			ivec(1) = ss22 - ss11
108			ivec(2) =(ss22 ** 2
109			& - ss11 ** 2)*(1. _d 0 / 2. _d 0)
110			ivec(3) =(ss22 ** 3
111			& - ss11 ** 3)*(1. _d 0 / 3. _d 0)
112
113			intF = ivec(1) * fhat(1,ir-0)
114			& + ivec(2) * fhat(2,ir-0)
115			& + ivec(3) * fhat(3,ir-0)
116
117			intF = intF / (ss22 - ss11)
118
119			end if
120
121			C ==================== calc. flux = upwind tracer * face-transport
122			flux(ix,iy,ir) = + wfac(ix,iy,ir) * intF
123
124			end if
125
126			end do
127
128			return
129
130			c end subroutine GAD_PPM_FLX_R
131			end