/[MITgcm]/MITgcm/pkg/generic_advdiff/gad_ppm_adv_r.F
ViewVC logotype

Contents of /MITgcm/pkg/generic_advdiff/gad_ppm_adv_r.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (show annotations) (download)
Sun Mar 13 01:44:02 2016 UTC (8 years, 1 month ago) by jmc
Branch: MAIN
CVS Tags: checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65v, checkpoint65w, checkpoint65u, HEAD
- from Darren: add PPM and PQM advection schemes (number 40-42 and 50-52)
  with 2 types of limiter (see: Engwirda & Kelley, submit. to JCP);
  Note (from Darren): unlimited PPM/PQM scheme (40 & 50) are just for
  testing and not for actual use.

1 C $Header: $
2 C $Name: $
3
4 # include "GAD_OPTIONS.h"
5
6 SUBROUTINE GAD_PPM_ADV_R(meth,bi,bj,
7 I delT,velR,facR,fbar,
8 O flux,myThid )
9 C |================================================================|
10 C | PPM_ADV_R: evaluate grid-cell advective flux in R. |
11 C | Lagrangian-type Piecewise Parabolic Method (PPM). |
12 C |================================================================|
13
14 implicit none
15
16 C =============================================== global variables
17 # include "SIZE.h"
18 # include "GRID.h"
19 # include "GAD.h"
20
21 C ================================================================
22 C meth :: advection method.
23 C bi,bj :: tile indexing
24 C delT :: level-wise time-steps.
25 C velR :: vel. field in r-direction.
26 C facR :: grid-areas in r-direction.
27 C fbar :: grid-cell values.
28 C flux :: trac.-flux in r-direction.
29 C myThid :: thread number.
30 C ================================================================
31 integer meth
32 integer bi,bj
33 _RL delT(1:Nr)
34 _RL velR(1-OLx:sNx+OLx,
35 & 1-OLy:sNy+OLy,
36 & 1:Nr)
37 _RL facR(1-OLx:sNx+OLx,
38 & 1-OLy:sNy+OLy,
39 & 1:Nr)
40 _RL fbar(1-OLx:sNx+OLx,
41 & 1-OLy:sNy+OLy,
42 & 1:Nr)
43 _RL flux(1-OLx:sNx+OLx,
44 & 1-OLy:sNy+OLy,
45 & 1:Nr)
46 integer myThid
47
48 C ================================================================
49 C ix,iy,ir :: grid indexing.
50 C floc :: col. of grid-cell values.
51 C mloc :: col. of grid-cell mask values.
52 C fhat :: col. of poly. coeff.
53 C - FHAT(:,I) = PQM coeff.
54 C edge :: col. of edge-wise values/slopes.
55 C - EDGE(1,:) = VALUE.
56 C - EDGE(2,:) = DF/DR.
57 C ohat :: col. of oscl. coeff.
58 C - OHAT(1,:) = D^1F/DS^1.
59 C - OHAT(2,:) = D^2F/DS^2.
60 C ================================================================
61 integer ix,iy,ir
62 _RL floc( 1-3:Nr+3)
63 _RL mloc( 1-3:Nr+3)
64 _RL fhat(1:3,1-0:Nr+0)
65 _RL edge( 1-0:Nr+1)
66 _RL ohat(1:2,1-3:Nr+3)
67 _RL vsum
68
69 C ======================================= mask boundary conditions
70 mloc( -2) = 0.0 _d 0
71 mloc( -1) = 0.0 _d 0
72 mloc( +0) = 0.0 _d 0
73 mloc(Nr+1) = 0.0 _d 0
74 mloc(Nr+2) = 0.0 _d 0
75 mloc(Nr+3) = 0.0 _d 0
76
77 C ================================================================
78 C (1): copy a single row of data onto contiguous storage, treat
79 C as a set of one-dimensional problems.
80 C (2): calc. "oscillation-indicators" for each grid-cell if ad-
81 C vection scheme is WENO-class.
82 C (3): calc. edge-centred values/slopes by high-order interpol-
83 C ation.
84 C (4): calc. cell-centred polynomial profiles with appropriate
85 C slope-limiting.
86 C (5): calc. fluxes using a local, semi-lagrangian integration.
87 C ================================================================
88
89 do iy = 1-OLy+0, sNy+OLy-0
90 do ix = 1-OLx+0, sNx+OLx-0
91 C ======================================= no flux through surf. bc
92 flux(ix,iy,+1) = 0.0 _d +0
93 end do
94 end do
95
96 C ==================== calculate transport for interior grid-cells
97 do iy = 1-OLy+0, sNy+OLy-0
98 do ix = 1-OLx+0, sNx+OLx-0
99
100 vsum = 0.0 _d 0
101 do ir = 2, Nr
102 C ================================== quick break on zero transport
103 vsum = vsum
104 & + abs(velR(ix,iy,ir))
105 end do
106
107 if (vsum .gt. 0. _d 0) then
108
109 do ir = 1, Nr
110 C ================================== make local unit-stride copies
111 floc(ir) = fbar(ix,iy,ir)
112 mloc(ir) =
113 & maskC(ix,iy,ir,bi,bj)
114 end do
115
116 C ================================== make mask boundary conditions
117 floc( -2) = floc(+1)
118 floc( -1) = floc(+1)
119 floc( +0) = floc(+1)
120 floc(Nr+1) = floc(Nr)
121 floc(Nr+2) = floc(Nr)
122 floc(Nr+3) = floc(Nr)
123
124 C ==================== reconstruct derivatives for WENO indicators
125 if (meth.eq.ENUM_PPM_WENO_LIMIT) then
126 CALL GAD_OSC_HAT_R(bi,bj,ix,iy,
127 & mloc,floc,
128 & ohat,myThid)
129 end if
130
131 C ==================== reconstruct 3rd--order accurate edge values
132 CALL GAD_PPM_P3E_R(bi,bj,ix,iy,
133 & mloc,floc,
134 & edge,myThid)
135
136 C ==================== reconstruct coeff. for grid-cell poynomials
137 CALL GAD_PPM_HAT_R(bi,bj,ix,iy,
138 & meth,
139 & mloc,floc,
140 & edge,ohat,
141 & fhat,myThid)
142
143 C ==================== evaluate integral fluxes on grid-cell edges
144 CALL GAD_PPM_FLX_R(bi,bj,ix,iy,
145 & delT,velR,
146 & facR,fhat,
147 & flux,myThid)
148
149 else
150
151 do ir = 2, Nr
152 C ================================== "null" flux on zero transport
153 flux(ix,iy,ir) = 0.0 _d +0
154 end do
155
156 end if
157
158 end do
159 end do
160
161 return
162
163 c end subroutine GAD_PPM_ADV_R
164 end

  ViewVC Help
Powered by ViewVC 1.1.22