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

Contents of /MITgcm/pkg/generic_advdiff/gad_ppm_adv_x.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:03 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_X(meth,bi,bj,kk,
7 I calc_CFL,delT,uvel,ufac,fbar,
8 O flux,myThid )
9 C |================================================================|
10 C | PPM_ADV_X: evaluate grid-cell advective flux in X. |
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 kk :: r-index.
25 C calc_CFL :: TRUE to calc. CFL from vel.
26 C delT :: time-step.
27 C uvel :: vel.-comp in x-direction.
28 C ufac :: vel.-flux in x-direction.
29 C fbar :: grid-cell values.
30 C flux :: adv.-flux in x-direction.
31 C myThid :: thread number.
32 C ================================================================
33 integer meth
34 integer bi,bj,kk
35 logical calc_CFL
36 _RL delT
37 _RL uvel(1-OLx:sNx+OLx,
38 & 1-OLy:sNy+OLy)
39 _RL ufac(1-OLx:sNx+OLx,
40 & 1-OLy:sNy+OLy)
41 _RL fbar(1-OLx:sNx+OLx,
42 & 1-OLy:sNy+OLy)
43 _RL flux(1-OLx:sNx+OLx,
44 & 1-OLy:sNy+OLy)
45 integer myThid
46
47 C ================================================================
48 C ix,iy,ir :: grid indexing.
49 C floc :: row of grid-cell values.
50 C mloc :: row of grid-cell mask values.
51 C fhat :: row of poly. coeff.
52 C - FHAT(:,I) = PQM coeff.
53 C edge :: row of edge-wise values/slopes.
54 C - EDGE(1,:) = VALUE.
55 C - EDGE(2,:) = DF/DY.
56 C ohat :: row of oscl. coeff.
57 C - OHAT(1,:) = D^1F/DS^1.
58 C - OHAT(2,:) = D^2F/DS^2.
59 C ================================================================
60 integer ix,iy
61 _RL mloc(1-OLx:sNx+OLx)
62 _RL floc(1-OLx:sNx+OLx)
63 _RL fhat(1:3,
64 & 1-OLx:sNx+OLx)
65 _RL edge(1-OLx:sNx+OLx)
66 _RL ohat(1:2,
67 & 1-OLx:sNx+OLx)
68 _RL vsum
69
70 do iy = 1-OLy+0, sNy+OLy-0
71 C ==================== zero stencil "ghost" cells along boundaries
72 flux( +1-OLx+0,iy) = 0. _d 0
73 flux( +1-OLx+1,iy) = 0. _d 0
74 flux( +1-OLx+2,iy) = 0. _d 0
75 flux( +1-OLx+3,iy) = 0. _d 0
76 flux(sNx+OLx-0,iy) = 0. _d 0
77 flux(sNx+OLx-1,iy) = 0. _d 0
78 flux(sNx+OLx-2,iy) = 0. _d 0
79 end do
80
81 C ================================================================
82 C (1): copy a single row of data onto contiguous storage, treat
83 C as a set of one-dimensional problems.
84 C (2): calc. "oscillation-indicators" for each grid-cell if ad-
85 C vection scheme is WENO-class.
86 C (3): calc. edge-centred values/slopes by high-order interpol-
87 C ation.
88 C (4): calc. cell-centred polynomial profiles with appropriate
89 C slope-limiting.
90 C (5): calc. fluxes using a local, semi-lagrangian integration.
91 C ================================================================
92
93 do iy = 1-OLy+0, sNy+OLy-0
94
95 vsum = 0.0 _d 0
96 do ix = 1-OLx+0, sNx+OLx-0
97 C ================================== quick break on zero transport
98 vsum = vsum
99 & + abs(ufac(ix,iy))
100 end do
101
102 if (vsum .gt. 0. _d 0) then
103
104 do ix = 1-OLx+0, sNx+OLx-0
105 C ================================== make local unit-stride copies
106 floc(ix) = fbar (ix,iy)
107 mloc(ix) =
108 & maskC(ix,iy,kk,bi,bj)
109 end do
110
111 C ==================== reconstruct derivatives for WENO indicators
112 if (meth.eq.ENUM_PPM_WENO_LIMIT) then
113 CALL GAD_OSC_HAT_X(bi,bj,kk,iy,
114 & mloc,floc,
115 & ohat,myThid)
116 end if
117
118 C ==================== reconstruct 3rd--order accurate edge values
119 CALL GAD_PPM_P3E_X(bi,bj,kk,iy,
120 & mloc,floc,
121 & edge,myThid)
122
123 C ==================== reconstruct coeff. for grid-cell poynomials
124 CALL GAD_PPM_HAT_X(bi,bj,kk,iy,
125 & meth,
126 & mloc,floc,
127 & edge,ohat,
128 & fhat,myThid)
129
130 C ==================== evaluate integral fluxes on grid-cell edges
131 CALL GAD_PPM_FLX_X(bi,bj,kk,iy,
132 & calc_CFL,
133 & delT,uvel,
134 & ufac,fhat,
135 & flux,myThid)
136
137 else
138
139 do ix = 1-OLx+3, sNx+OLx-2
140 C ================================== "null" flux on zero transport
141 flux(ix,iy) = 0.0 _d 0
142 end do
143
144 end if
145
146 end do
147
148 return
149
150 c end subroutine GAD_PPM_ADV_X
151 end

  ViewVC Help
Powered by ViewVC 1.1.22