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

Contents of /MITgcm/pkg/generic_advdiff/gad_pqm_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: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_PQM_ADV_R(meth,bi,bj,
7 I delT,wvel,wfac,fbar,
8 O flux,myThid )
9 C |================================================================|
10 C | PQM_ADV_R: evaluate grid-cell advective flux in R. |
11 C | Lagrangian-type Piecewise Quartic Method (PQM). |
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 wvel :: vel.-comp in r-direction.
26 C wfac :: vel.-flux in r-direction.
27 C fbar :: grid-cell values.
28 C flux :: adv.-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 wvel(1-OLx:sNx+OLx,
35 & 1-OLy:sNy+OLy,
36 & 1:Nr)
37 _RL wfac(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:5,1-0:Nr+0)
65 _RL edge(1:2,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 = +1, Nr
102 C ================================== quick break on zero transport
103 vsum = vsum
104 & + abs(wfac(ix,iy,ir))
105 C ================================== make local unit-stride copies
106 floc(ir) = fbar (ix,iy,ir)
107 mloc(ir) =
108 & maskC(ix,iy,ir,bi,bj)
109 end do
110
111 if (vsum .gt. 0. _d 0) then
112
113 C ================================== make mask boundary conditions
114 floc( -2) = floc(+1)
115 floc( -1) = floc(+1)
116 floc( +0) = floc(+1)
117 floc(Nr+1) = floc(Nr)
118 floc(Nr+2) = floc(Nr)
119 floc(Nr+3) = floc(Nr)
120
121 C ==================== reconstruct derivatives for WENO indicators
122 if (meth.eq.ENUM_PQM_WENO_LIMIT) then
123 CALL GAD_OSC_HAT_R(bi,bj,ix,iy,
124 & mloc,floc,
125 & ohat,myThid)
126 end if
127
128 C ==================== reconstruct 5th--order accurate edge values
129 CALL GAD_PQM_P5E_R(bi,bj,ix,iy,
130 & mloc,floc,
131 & edge,myThid)
132
133 C ==================== reconstruct coeff. for grid-cell poynomials
134 CALL GAD_PQM_HAT_R(bi,bj,ix,iy,
135 & meth,
136 & mloc,floc,
137 & edge,ohat,
138 & fhat,myThid)
139
140 C ==================== evaluate integral fluxes on grid-cell edges
141 CALL GAD_PQM_FLX_R(bi,bj,ix,iy,
142 & delT,wvel,
143 & wfac,fhat,
144 & flux,myThid)
145
146 else
147
148 do ir = 2, Nr
149 C ================================== "null" flux on zero transport
150 flux(ix,iy,ir) = 0. _d 0
151 end do
152
153 end if
154
155 end do
156 end do
157
158 return
159
160 c end subroutine GAD_PQM_ADV_R
161 end

  ViewVC Help
Powered by ViewVC 1.1.22