pkg/generic_advdiff/gad_ppm_flx_x.F

C $Header:  $
C $Name:  $

#     include "GAD_OPTIONS.h"

      SUBROUTINE GAD_PPM_FLX_X(bi,bj,kk,iy,
     &           calc_CFL,delT,uvel,
     &           ufac,fhat,flux,myThid)
C     |================================================================|
C     | PPM_FLX_X: 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       kk        :: r-index.
C       iy        :: y-index.
C       calc_CFL  :: TRUE to calc. CFL from vel.
C       delT      :: time-step.
C       uvel      :: vel.-comp in x-direction.
C       ufac      :: vel.-flux in x-direction.
C       fhat      :: row of poly. coeff.
C       flux      :: adv.-flux in x-direction.
C       myThid    :: thread number.
C     ================================================================
          integer bi,bj,kk,iy
          logical calc_CFL
          _RL delT
          _RL uvel(1-OLx:sNx+OLx,
     &             1-OLy:sNy+OLy)
          _RL ufac(1-OLx:sNx+OLx,
     &             1-OLy:sNy+OLy)
          _RL fhat(1:3,
     &             1-OLx:sNx+OLx)
          _RL flux(1-OLx:sNx+OLx,
     &             1-OLy:sNy+OLy)
          integer myThid

C     ================================================================
C       ix        :: x-indexing.
C       uCFL      :: CFL number.
C       intF      :: upwind tracer edge-value.
C       ss11,ss22 :: int. endpoints.
C       ivec      :: int. basis vec.
C     ================================================================
          integer ix
          _RL uCFL,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  ix = 1-OLx+3, sNx+OLx-2

              if (uvel(ix,iy) .eq. 0. _d 0) then

                  flux(ix,iy)    = 0. _d 0

              else

              if (uvel(ix,iy) .gt. 0. _d 0) then

C     ==================== integrate PPM profile over upwind cell IX-1
              if ( calc_CFL ) then
              uCFL = uvel(ix,iy) * delT
     &             * recip_dxF(ix-1,iy,bi,bj)
     &             * recip_deepFacC(kk)
              else
              uCFL = uvel(ix,iy)
              end if

              ss11 = +1. _d 0 - 2. _d 0 * uCFL
              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,ix-1)
     &             + ivec(2) * fhat(2,ix-1)
     &             + ivec(3) * fhat(3,ix-1)

              intF = intF / (ss22 - ss11)

              else

C     ==================== integrate PPM profile over upwind cell IX+0
              if ( calc_CFL ) then
              uCFL = uvel(ix,iy) * delT
     &             * recip_dxF(ix-0,iy,bi,bj)
     &             * recip_deepFacC(kk)
              else
              uCFL = uvel(ix,iy)
              end if

              ss11 = -1. _d 0 - 2. _d 0 * uCFL
              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,ix-0)
     &             + ivec(2) * fhat(2,ix-0)
     &             + ivec(3) * fhat(3,ix-0)

              intF = intF / (ss22 - ss11)

              end if

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

              end if

          end do

          return

c     end subroutine GAD_PPM_FLX_X
      end
1	C $Header: $
2	C $Name: $
3
4	# include "GAD_OPTIONS.h"
5
6	SUBROUTINE GAD_PPM_FLX_X(bi,bj,kk,iy,
7	& calc_CFL,delT,uvel,
8	& ufac,fhat,flux,myThid)
9	C \|================================================================\|
10	C \| PPM_FLX_X: 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 kk :: r-index.
23	C iy :: y-index.
24	C calc_CFL :: TRUE to calc. CFL from vel.
25	C delT :: time-step.
26	C uvel :: vel.-comp in x-direction.
27	C ufac :: vel.-flux in x-direction.
28	C fhat :: row of poly. coeff.
29	C flux :: adv.-flux in x-direction.
30	C myThid :: thread number.
31	C ================================================================
32	integer bi,bj,kk,iy
33	logical calc_CFL
34	_RL delT
35	_RL uvel(1-OLx:sNx+OLx,
36	& 1-OLy:sNy+OLy)
37	_RL ufac(1-OLx:sNx+OLx,
38	& 1-OLy:sNy+OLy)
39	_RL fhat(1:3,
40	& 1-OLx:sNx+OLx)
41	_RL flux(1-OLx:sNx+OLx,
42	& 1-OLy:sNy+OLy)
43	integer myThid
44
45	C ================================================================
46	C ix :: x-indexing.
47	C uCFL :: CFL number.
48	C intF :: upwind tracer edge-value.
49	C ss11,ss22 :: int. endpoints.
50	C ivec :: int. basis vec.
51	C ================================================================
52	integer ix
53	_RL uCFL,intF
54	_RL ss11,ss22
55	_RL ivec(1:3)
56
57	C ================================================================
58	C (1): calc. "departure-points" for each grid-cell edge by int-
59	C egrating edge position backward in time over one single
60	C time-step. This is a "single-cell" implementation: requ-
61	C ires CFL <= 1.0.
62	C (2): calc. flux as the integral of the upwind grid-cell poly-
63	C nomial over the deformation interval found in (1).
64	C ================================================================
65
66	do ix = 1-OLx+3, sNx+OLx-2
67
68	if (uvel(ix,iy) .eq. 0. _d 0) then
69
70	flux(ix,iy) = 0. _d 0
71
72	else
73
74	if (uvel(ix,iy) .gt. 0. _d 0) then
75
76	C ==================== integrate PPM profile over upwind cell IX-1
77	if ( calc_CFL ) then
78	uCFL = uvel(ix,iy) * delT
79	& * recip_dxF(ix-1,iy,bi,bj)
80	& * recip_deepFacC(kk)
81	else
82	uCFL = uvel(ix,iy)
83	end if
84
85	ss11 = +1. _d 0 - 2. _d 0 * uCFL
86	ss22 = +1. _d 0
87
88	C ==================== integrate profile over region swept by face
89	ivec(1) = ss22 - ss11
90	ivec(2) =(ss22 ** 2
91	& - ss11 ** 2)*(1. _d 0 / 2. _d 0)
92	ivec(3) =(ss22 ** 3
93	& - ss11 ** 3)*(1. _d 0 / 3. _d 0)
94
95	intF = ivec(1) * fhat(1,ix-1)
96	& + ivec(2) * fhat(2,ix-1)
97	& + ivec(3) * fhat(3,ix-1)
98
99	intF = intF / (ss22 - ss11)
100
101	else
102
103	C ==================== integrate PPM profile over upwind cell IX+0
104	if ( calc_CFL ) then
105	uCFL = uvel(ix,iy) * delT
106	& * recip_dxF(ix-0,iy,bi,bj)
107	& * recip_deepFacC(kk)
108	else
109	uCFL = uvel(ix,iy)
110	end if
111
112	ss11 = -1. _d 0 - 2. _d 0 * uCFL
113	ss22 = -1. _d 0
114
115	C ==================== integrate profile over region swept by face
116	ivec(1) = ss22 - ss11
117	ivec(2) =(ss22 ** 2
118	& - ss11 ** 2)*(1. _d 0 / 2. _d 0)
119	ivec(3) =(ss22 ** 3
120	& - ss11 ** 3)*(1. _d 0 / 3. _d 0)
121
122	intF = ivec(1) * fhat(1,ix-0)
123	& + ivec(2) * fhat(2,ix-0)
124	& + ivec(3) * fhat(3,ix-0)
125
126	intF = intF / (ss22 - ss11)
127
128	end if
129
130	C ==================== calc. flux = upwind tracer * face-transport
131	flux(ix,iy) = + ufac(ix,iy) * intF
132
133	end if
134
135	end do
136
137	return
138
139	c end subroutine GAD_PPM_FLX_X
140	end