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

Contents of /MITgcm/pkg/generic_advdiff/gad_pqm_hat_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_PQM_HAT_X(bi,bj,kk,iy,
7 & method,mask,fbar,edge,
8 & ohat,fhat,myThid)
9 C |================================================================|
10 C | PQM_HAT_X: reconstruct grid-cell PQM polynomials. |
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 ====================================================== arguments
21 C bi,bj :: tile indices.
22 C kk :: r-index.
23 C iy :: y-index.
24 C method :: advection scheme.
25 C mask :: row of cell-wise mask values.
26 C fbar :: row of cell-wise values.
27 C edge :: row of edge-wise values/slopes.
28 C - EDGE(1,:) = VALUE.
29 C - EDGE(2,:) = DF/DX.
30 C ohat :: row of oscl. coeff.
31 C - OHAT(1,:) = D^1F/DS^1.
32 C - OHAT(2,:) = D^2F/DS^2.
33 C fhat :: row of poly. coeff.
34 C - FHAT(:,I) = PQM coeff.
35 C myThid :: thread number.
36 C ================================================================
37 integer bi,bj,kk,iy
38 integer method
39 _RL mask(1-OLx:sNx+OLx)
40 _RL fbar(1-OLx:sNx+OLx)
41 _RL edge(1:2,
42 & 1-OLx:sNx+OLx)
43 _RL ohat(1:2,
44 & 1-OLx:sNx+OLx)
45 _RL fhat(1:5,
46 & 1-OLx:sNx+OLx)
47 integer myThid
48
49 C ====================================================== variables
50 C ii,ix :: local x-indexing.
51 C ff00 :: centre-biased cell value.
52 C ffll,ffrr :: left-, and right-biased cell values.
53 C xhat :: local coord. scaling.
54 C fell,ferr :: left-, and right-biased edge values.
55 C dell,derr :: left-, and right-biased edge slopes.
56 C dfds :: linear slope estimates.
57 C uhat :: "NULL" limited poly. coeff.
58 C lhat :: "MONO" limited poly. coeff.
59 C scal :: "WENO" weights.
60 C fmag,fdel :: local perturbation indicators.
61 C ================================================================
62 integer ii,ix
63 _RL ff00
64 _RL ffll,ffrr
65 _RL xhat
66 _RL fell,ferr
67 _RL dell,derr
68 _RL dfds(-1:+1)
69 _RL uhat(+1:+5)
70 _RL lhat(+1:+5)
71 _RL scal(+1:+2)
72 _RL fmag,fdel
73 integer mono
74
75 do ix = 1-OLx+3, sNx+OLx-3
76
77 if (mask(ix) .gt. 0. _d 0) then
78
79 C =============================== scale to local grid-cell co-ords
80 xhat = dxF(ix,iy,bi,bj) * 0.5 _d 0
81
82 C =============================== assemble cell mean + edge values
83 ff00 = fbar(ix+0)
84 ffll = ff00
85 & + mask(ix-1)*(fbar(ix-1)-ff00)
86 ffrr = ff00
87 & + mask(ix+1)*(fbar(ix+1)-ff00)
88
89 fell = edge(+1,ix-0)
90 ferr = edge(+1,ix+1)
91
92 dell = edge(+2,ix-0)
93 derr = edge(+2,ix+1)
94
95 dell = dell * xhat
96 derr = derr * xhat
97
98 c select case(method)
99 c case(ENUM_PQM_NULL_LIMIT)
100 if ( method.eq.ENUM_PQM_NULL_LIMIT ) then
101 C =============================== "NULL" limited grid-cell profile
102 CALL GAD_PQM_FUN_NULL ( ff00,
103 & fell,ferr,dell,derr,lhat,mono)
104
105 c case(ENUM_PQM_MONO_LIMIT)
106 elseif ( method.eq.ENUM_PQM_MONO_LIMIT ) then
107 C =============================== "MONO" limited grid-cell profile
108 CALL GAD_PLM_FUN_U(ffll,ff00,ffrr,dfds)
109
110 CALL GAD_PQM_FUN_MONO ( ff00,ffll,ffrr,
111 & fell,ferr,dell,derr,dfds,lhat,
112 & mono)
113
114 c case(ENUM_PQM_WENO_LIMIT)
115 elseif ( method.eq.ENUM_PQM_WENO_LIMIT ) then
116 C =============================== "WENO" limited grid-cell profile
117 CALL GAD_PLM_FUN_U(ffll,ff00,ffrr,dfds)
118
119 CALL GAD_PQM_FUN_NULL ( ff00,
120 & fell,ferr,dell,derr,uhat,mono)
121
122 CALL GAD_PQM_FUN_MONO ( ff00,ffll,ffrr,
123 & fell,ferr,dell,derr,dfds,lhat,
124 & mono)
125
126 if ( mono .gt. 0) then
127
128 C =============================== only apply WENO if it is worth it
129 fdel = abs(ffrr-ff00)+abs(ff00-ffll)
130 fmag = abs(ffll)+abs(ff00)+abs(ffrr)
131
132 if (fdel .gt. 1. _d -6 * fmag) then
133
134 C =============================== calc. WENO oscillation weighting
135 CALL GAD_OSC_MUL_X(ix,+2,mask,
136 & ohat,scal)
137
138 do ii = +1, +5
139 C =============================== blend limited/un-limited profile
140 lhat(ii) = scal(1) * uhat(ii)
141 & + scal(2) * lhat(ii)
142 end do
143
144 end if
145
146 end if
147
148 c end select
149 endif
150
151 do ii = +1, +5
152 C =============================== copy polynomial onto output data
153 fhat(ii,ix) = lhat(ii)
154 end do
155
156 else
157
158 do ii = +1, +5
159 fhat(ii,ix) = 0.0 _d 0
160 end do
161
162 end if
163
164 end do
165
166 return
167
168 c end subroutine GAD_PQM_HAT_X
169 end

  ViewVC Help
Powered by ViewVC 1.1.22