New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

plm.h90 in vendors/PPR/src – NEMO

Context Navigation

source: vendors/PPR/src/plm.h90 @ 14212

Last change on this file since 14212 was 14212, checked in by jchanut, 4 years ago
#2222, update PPR library to circumvent issues found with cpp directives: remove timing instructions. Changed main module extension in .F90 to comply with FCM rules (no effect).
File size: 12.6 KB

Line
1
2	!
3	! This program may be freely redistributed under the
4	! condition that the copyright notices (including this
5	! entire header) are not removed, and no compensation
6	! is received through use of the software. Private,
7	! research, and institutional use is free. You may
8	! distribute modified versions of this code UNDER THE
9	! CONDITION THAT THIS CODE AND ANY MODIFICATIONS MADE
10	! TO IT IN THE SAME FILE REMAIN UNDER COPYRIGHT OF THE
11	! ORIGINAL AUTHOR, BOTH SOURCE AND OBJECT CODE ARE
12	! MADE FREELY AVAILABLE WITHOUT CHARGE, AND CLEAR
13	! NOTICE IS GIVEN OF THE MODIFICATIONS. Distribution
14	! of this code as part of a commercial system is
15	! permissible ONLY BY DIRECT ARRANGEMENT WITH THE
16	! AUTHOR. (If you are not directly supplying this
17	! code to a customer, and you are instead telling them
18	! how they can obtain it for free, then you are not
19	! required to make any arrangement with me.)
20	!
21	! Disclaimer: Neither I nor: Columbia University, the
22	! National Aeronautics and Space Administration, nor
23	! the Massachusetts Institute of Technology warrant
24	! or certify this code in any way whatsoever. This
25	! code is provided "as-is" to be used at your own risk.
26	!
27	!
28
29	!
30	! PLM.h90: a 1d, slope-limited piecewise linear method.
31	!
32	! Darren Engwirda
33	! 08-Sep-2016
34	! de2363 [at] columbia [dot] edu
35	!
36	!
37
38	pure subroutine plm(npos,nvar,ndof,delx, &
39	& fdat,fhat,dmin,ilim)
40
41	!
42	! NPOS no. edges over grid.
43	! NVAR no. state variables.
44	! NDOF no. degrees-of-freedom per grid-cell .
45	! DELX grid-cell spacing array. LENGTH(DELX) == +1 if
46	! spacing is uniform .
47	! FDAT grid-cell moments array. FDAT is an array with
48	! SIZE = NDOF-by-NVAR-by-NPOS-1 .
49	! FHAT grid-cell re-con. array. FHAT is an array with
50	! SIZE = MDOF-by-NVAR-by-NPOS-1 .
51	! DMIN min. grid-cell spacing thresh .
52	! ILIM cell slope-limiting selection .
53	!
54
55	implicit none
56
57	!------------------------------------------- arguments !
58	integer, intent( in) :: npos,nvar
59	integer, intent( in) :: ndof,ilim
60	real*8 , intent( in) :: dmin
61	real*8 , intent( in) :: delx(:)
62	real*8 , intent(out) :: fhat(:,:,:)
63	real*8 , intent( in) :: fdat(:,:,:)
64
65	if (size(delx).gt.+1) then
66
67	!------------------------------- variable grid-spacing !
68
69	call plmv(npos,nvar,ndof,delx,&
70	& fdat,fhat,&
71	& dmin,ilim )
72
73	else
74
75	!------------------------------- constant grid-spacing !
76
77	call plmc(npos,nvar,ndof,delx,&
78	& fdat,fhat,&
79	& dmin,ilim )
80
81	end if
82
83	return
84
85	end subroutine
86
87	!------------------------- assemble PLM reconstruction !
88
89	pure subroutine plmv(npos,nvar,ndof,delx, &
90	& fdat,fhat,dmin,ilim)
91
92	!
93	! *this is the variable grid-spacing variant .
94	!
95	! NPOS no. edges over grid.
96	! NVAR no. state variables.
97	! NDOF no. degrees-of-freedom per grid-cell .
98	! DELX grid-cell spacing array. LENGTH(DELX) == +1 if
99	! spacing is uniform .
100	! FDAT grid-cell moments array. FDAT is an array with
101	! SIZE = NDOF-by-NVAR-by-NPOS-1 .
102	! FHAT grid-cell re-con. array. FHAT is an array with
103	! SIZE = MDOF-by-NVAR-by-NPOS-1 .
104	! DMIN min. grid-cell spacing thresh .
105	! ILIM cell slope-limiting selection .
106	!
107
108	implicit none
109
110	!------------------------------------------- arguments !
111	integer, intent( in) :: npos,nvar
112	integer, intent( in) :: ndof,ilim
113	real*8 , intent( in) :: dmin
114	real*8 , intent( in) :: delx(:)
115	real*8 , intent(out) :: fhat(:,:,:)
116	real*8 , intent( in) :: fdat(:,:,:)
117
118	!------------------------------------------- variables !
119	integer :: ipos,ivar,head,tail
120	real*8 :: dfds(-1:+1)
121
122	head = +1; tail = npos - 1
123
124	if (npos.eq.2) then
125	!----------------------- reduce order if small stencil !
126	do ivar = +1, nvar
127	fhat(1,ivar,1) = &
128	& fdat(1,ivar,1)
129	fhat(2,ivar,1) = 0.d+0
130	end do
131	end if
132
133	if (npos.le.2) return
134
135	!-------------------------------------- lower-endpoint !
136
137	do ivar = +1 , nvar-0
138
139	call plsv( &
140	& fdat(1,ivar,head+0) , &
141	& delx(head+0), &
142	& fdat(1,ivar,head+0) , &
143	& delx(head+0), &
144	& fdat(1,ivar,head+1) , &
145	& delx(head+1), dfds)
146
147	fhat(1,ivar,head) = &
148	& fdat(1,ivar,head)
149	fhat(2,ivar,head) = dfds(0)
150
151	end do
152
153	!-------------------------------------- upper-endpoint !
154
155	do ivar = +1 , nvar-0
156
157	call plsv( &
158	& fdat(1,ivar,tail-1) , &
159	& delx(tail-1), &
160	& fdat(1,ivar,tail+0) , &
161	& delx(tail+0), &
162	& fdat(1,ivar,tail+0) , &
163	& delx(tail+0), dfds)
164
165	fhat(1,ivar,tail) = &
166	& fdat(1,ivar,tail)
167	fhat(2,ivar,tail) = dfds(0)
168
169	end do
170
171	!-------------------------------------- interior cells !
172
173	do ipos = +2 , npos-2
174	do ivar = +1 , nvar-0
175
176	call plsv( &
177	& fdat(1,ivar,ipos-1) , &
178	& delx(ipos-1), &
179	& fdat(1,ivar,ipos+0) , &
180	& delx(ipos+0), &
181	& fdat(1,ivar,ipos+1) , &
182	& delx(ipos+1), dfds)
183
184	fhat(1,ivar,ipos) = &
185	& fdat(1,ivar,ipos)
186	fhat(2,ivar,ipos) = dfds(0)
187
188	end do
189	end do
190
191	return
192
193	end subroutine
194
195	!------------------------- assemble PLM reconstruction !
196
197	pure subroutine plmc(npos,nvar,ndof,delx, &
198	& fdat,fhat,dmin,ilim)
199
200	!
201	! *this is the constant grid-spacing variant .
202	!
203	! NPOS no. edges over grid.
204	! NVAR no. state variables.
205	! NDOF no. degrees-of-freedom per grid-cell .
206	! DELX grid-cell spacing array. LENGTH(DELX) == +1 if
207	! spacing is uniform .
208	! FDAT grid-cell moments array. FDAT is an array with
209	! SIZE = NDOF-by-NVAR-by-NPOS-1 .
210	! FHAT grid-cell re-con. array. FHAT is an array with
211	! SIZE = MDOF-by-NVAR-by-NPOS-1 .
212	! DMIN min. grid-cell spacing thresh .
213	! ILIM cell slope-limiting selection .
214	!
215
216	implicit none
217
218	!------------------------------------------- arguments !
219	integer, intent( in) :: npos,nvar
220	integer, intent( in) :: ndof,ilim
221	real*8 , intent( in) :: dmin
222	real*8 , intent( in) :: delx(1)
223	real*8 , intent(out) :: fhat(:,:,:)
224	real*8 , intent( in) :: fdat(:,:,:)
225
226	!------------------------------------------- variables !
227	integer :: ipos,ivar,head,tail
228	real*8 :: dfds(-1:+1)
229
230	head = +1; tail = npos - 1
231
232	if (npos.eq.2) then
233	!----------------------- reduce order if small stencil !
234	do ivar = +1, nvar
235	fhat(1,ivar,1) = &
236	& fdat(1,ivar,1)
237	fhat(2,ivar,1) = 0.d+0
238	end do
239	end if
240
241	if (npos.le.2) return
242
243	!-------------------------------------- lower-endpoint !
244
245	do ivar = +1 , nvar-0
246
247	call plsc( &
248	& fdat(1,ivar,head+0) , &
249	& fdat(1,ivar,head+0) , &
250	& fdat(1,ivar,head+1) , &
251	& dfds)
252
253	fhat(1,ivar,head) = &
254	& fdat(1,ivar,head)
255	fhat(2,ivar,head) = dfds(0)
256
257	end do
258
259	!-------------------------------------- upper-endpoint !
260
261	do ivar = +1 , nvar-0
262
263	call plsc( &
264	& fdat(1,ivar,tail-1) , &
265	& fdat(1,ivar,tail+0) , &
266	& fdat(1,ivar,tail+0) , &
267	& dfds)
268
269	fhat(1,ivar,tail) = &
270	& fdat(1,ivar,tail)
271	fhat(2,ivar,tail) = dfds(0)
272
273	end do
274
275	!-------------------------------------- interior cells !
276
277	do ipos = +2 , npos-2
278	do ivar = +1 , nvar-0
279
280	call plsc( &
281	& fdat(1,ivar,ipos-1) , &
282	& fdat(1,ivar,ipos+0) , &
283	& fdat(1,ivar,ipos+1) , &
284	& dfds)
285
286	fhat(1,ivar,ipos) = &
287	& fdat(1,ivar,ipos)
288	fhat(2,ivar,ipos) = dfds(0)
289
290	end do
291	end do
292
293	return
294
295	end subroutine
296
297	!------------------------------- assemble PLM "slopes" !
298
299	pure subroutine plsv(ffll,hhll,ff00,hh00,&
300	& ffrr,hhrr,dfds)
301
302	!
303	! *this is the variable grid-spacing variant .
304	!
305	! FFLL left -biased grid-cell mean.
306	! HHLL left -biased grid-cell spac.
307	! FF00 centred grid-cell mean.
308	! HH00 centred grid-cell spac.
309	! FFRR right-biased grid-cell mean.
310	! HHRR right-biased grid-cell spac.
311	! DFDS piecewise linear gradients in local co-ord.'s.
312	! DFDS(+0) is a centred, slope-limited estimate,
313	! DFDS(-1), DFDS(+1) are left- and right-biased
314	! estimates (unlimited).
315	!
316
317	implicit none
318
319	!------------------------------------------- arguments !
320	real*8 , intent( in) :: ffll,ff00,ffrr
321	real*8 , intent( in) :: hhll,hh00,hhrr
322	real*8 , intent(out) :: dfds(-1:+1)
323
324	!------------------------------------------- variables !
325	real*8 :: fell,ferr,scal
326
327	real*8 , parameter :: ZERO = 1.d-14
328
329	!---------------------------- 2nd-order approximations !
330
331	dfds(-1) = ff00-ffll
332	dfds(+1) = ffrr-ff00
333
334	if (dfds(-1) * &
335	& dfds(+1) .gt. 0.0d+0) then
336
337	!---------------------------- calc. ll//rr edge values !
338
339	fell = (hh00ffll+hhllff00) &
340	& / (hhll+hh00)
341	ferr = (hhrrff00+hh00ffrr) &
342	& / (hh00+hhrr)
343
344	!---------------------------- calc. centred derivative !
345
346	dfds(+0) = &
347	& 0.5d+0 * (ferr - fell)
348
349	!---------------------------- monotonic slope-limiting !
350
351	scal = min(abs(dfds(-1)), &
352	& abs(dfds(+1))) &
353	& / max(abs(dfds(+0)), &
354	ZERO)
355	scal = min(scal,+1.0d+0)
356
357	dfds(+0) = scal * dfds(+0)
358
359	else
360
361	!---------------------------- flatten if local extrema !
362
363	dfds(+0) = +0.0d+0
364
365	end if
366
367	!---------------------------- scale onto local co-ord. !
368
369	dfds(-1) = dfds(-1) &
370	& / (hhll + hh00) * hh00
371	dfds(+1) = dfds(+1) &
372	& / (hh00 + hhrr) * hh00
373
374	return
375
376	end subroutine
377
378	!------------------------------- assemble PLM "slopes" !
379
380	pure subroutine plsc(ffll,ff00,ffrr,dfds)
381
382	!
383	! *this is the constant grid-spacing variant .
384	!
385	! FFLL left -biased grid-cell mean.
386	! FF00 centred grid-cell mean.
387	! FFRR right-biased grid-cell mean.
388	! DFDS piecewise linear gradients in local co-ord.'s.
389	! DFDS(+0) is a centred, slope-limited estimate,
390	! DFDS(-1), DFDS(+1) are left- and right-biased
391	! estimates (unlimited).
392	!
393
394	implicit none
395
396	!------------------------------------------- arguments !
397	real*8 , intent( in) :: ffll,ff00,ffrr
398	real*8 , intent(out) :: dfds(-1:+1)
399
400	!------------------------------------------- variables !
401	real*8 :: fell,ferr,scal
402
403	real*8 , parameter :: ZERO = 1.d-14
404
405	!---------------------------- 2nd-order approximations !
406
407	dfds(-1) = ff00-ffll
408	dfds(+1) = ffrr-ff00
409
410	if (dfds(-1) * &
411	& dfds(+1) .gt. 0.0d+0) then
412
413	!---------------------------- calc. ll//rr edge values !
414
415	fell = (ffll+ff00) * .5d+0
416	ferr = (ff00+ffrr) * .5d+0
417
418	!---------------------------- calc. centred derivative !
419
420	dfds(+0) = &
421	& 0.5d+0 * (ferr - fell)
422
423	!---------------------------- monotonic slope-limiting !
424
425	scal = min(abs(dfds(-1)), &
426	& abs(dfds(+1))) &
427	& / max(abs(dfds(+0)), &
428	ZERO)
429	scal = min(scal,+1.0d+0)
430
431	dfds(+0) = scal * dfds(+0)
432
433	else
434
435	!---------------------------- flatten if local extrema !
436
437	dfds(+0) = +0.0d+0
438
439	end if
440
441	!---------------------------- scale onto local co-ord. !
442
443	dfds(-1) = + 0.5d+0 * dfds(-1)
444	dfds(+1) = + 0.5d+0 * dfds(+1)
445
446	return
447
448	end subroutine
449
450
451

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format