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.

pbc.h90 in NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/ext/PPR/src – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/ext/PPR/src/pbc.h90 @ 13926

Last change on this file since 13926 was 13926, checked in by jchanut, 4 years ago
#2222, add Piecewise Polynomial Reconstruction library
File size: 22.0 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	! PBC.f90: setup polynomial B.C.'s at domain endpoints.
31	!
32	! Darren Engwirda
33	! 09-Sep-2016
34	! de2363 [at] columbia [dot] edu
35	!
36	!
37
38	subroutine pbc(npos,nvar,ndof,delx, &
39	& fdat,bcon,edge,dfdx, &
40	& iend,dmin)
41
42	!
43	! NPOS no. edges over grid.
44	! NVAR no. state variables.
45	! NDOF no. degrees-of-freedom per grid-cell.
46	! DELX grid-cell spacing array. LENGTH(DELX) == +1 if
47	! spacing is uniform .
48	! FDAT grid-cell moments array. FDAT is an array with
49	! SIZE = NDOF-by-NVAR-by-NPOS-1 .
50	! BCON boundary condition data for endpoint .
51	! EDGE edge-centred interp. for function-value. EDGE
52	! is an array with SIZE = NVAR-by-NPOS .
53	! DFDX edge-centred interp. for 1st-derivative. DFDX
54	! is an array with SIZE = NVAR-by-NPOS .
55	! IEND domain endpoint, IEND < +0 for lower end-point
56	! and IEND > +0 for upper endpoint .
57	! DMIN min. grid-cell spacing thresh .
58	!
59
60	implicit none
61
62	!------------------------------------------- arguments !
63	integer, intent( in) :: npos,nvar,ndof
64	real*8 , intent( in) :: delx(:)
65	real*8 , intent( in) :: fdat(:,:,:)
66	real*8 , intent(out) :: edge(:,:)
67	real*8 , intent(out) :: dfdx(:,:)
68	integer, intent( in) :: iend
69	real*8 , intent( in) :: dmin
70	type(rcon_ends), intent(in) :: bcon(:)
71
72	!------------------------------------------- variables !
73	integer :: ivar,nlse,nval,nslp
74
75	nlse = 0 ; nval = 0 ; nslp = 0
76
77	do ivar = +1, nvar
78
79	select case (bcon(ivar)%bcopt)
80	!------------------------------------------- find BC's !
81	case(bcon_loose)
82	nlse = nlse + 1
83
84	case(bcon_value)
85	nval = nval + 1
86
87	case(bcon_slope)
88	nslp = nslp + 1
89
90	end select
91
92	end do
93
94	!---------------------------- setup "lower" conditions !
95
96	if (iend.lt.+0) then
97
98	if (nlse.gt.+0) then
99	!---------------------------- setup "unset" conditions !
100	call lbc(npos,nvar,ndof, &
101	& delx,fdat,bcon, &
102	& bcon_loose , &
103	& edge,dfdx,dmin)
104
105	end if
106
107	if (nval.gt.+0) then
108	!---------------------------- setup "value" conditions !
109	call lbc(npos,nvar,ndof, &
110	& delx,fdat,bcon, &
111	& bcon_value , &
112	& edge,dfdx,dmin)
113
114	end if
115
116	if (nslp.gt.+0) then
117	!---------------------------- setup "slope" conditions !
118	call lbc(npos,nvar,ndof, &
119	& delx,fdat,bcon, &
120	& bcon_slope , &
121	& edge,dfdx,dmin)
122
123	end if
124
125	end if
126
127	!---------------------------- setup "upper" conditions !
128
129	if (iend.gt.+0) then
130
131	if (nlse.gt.+0) then
132	!---------------------------- setup "unset" conditions !
133	call ubc(npos,nvar,ndof, &
134	& delx,fdat,bcon, &
135	& bcon_loose , &
136	& edge,dfdx,dmin)
137
138	end if
139
140	if (nval.gt.+0) then
141	!---------------------------- setup "value" conditions !
142	call ubc(npos,nvar,ndof, &
143	& delx,fdat,bcon, &
144	& bcon_value , &
145	& edge,dfdx,dmin)
146
147	end if
148
149	if (nslp.gt.+0) then
150	!---------------------------- setup "slope" conditions !
151	call ubc(npos,nvar,ndof, &
152	& delx,fdat,bcon, &
153	& bcon_slope , &
154	& edge,dfdx,dmin)
155
156	end if
157
158	end if
159
160	return
161
162	end subroutine
163
164	! LBC: impose a single B.C.-type at the lower endpoint !
165
166	subroutine lbc(npos,nvar,ndof,delx, &
167	& fdat,bcon,bopt,edge, &
168	& dfdx,dmin)
169
170	!
171	! NPOS no. edges over grid.
172	! NVAR no. state variables.
173	! NDOF no. degrees-of-freedom per grid-cell.
174	! DELX grid-cell spacing array. LENGTH(DELX) == +1 if
175	! spacing is uniform .
176	! FDAT grid-cell moments array. FDAT is an array with
177	! SIZE = NDOF-by-NVAR-by-NPOS-1 .
178	! BCON boundary condition data for endpoint .
179	! EDGE edge-centred interp. for function-value. EDGE
180	! is an array with SIZE = NVAR-by-NPOS .
181	! DFDX edge-centred interp. for 1st-derivative. DFDX
182	! is an array with SIZE = NVAR-by-NPOS .
183	! DMIN min. grid-cell spacing thresh .
184	!
185
186	implicit none
187
188	!------------------------------------------- arguments !
189	integer, intent( in) :: npos,nvar,ndof
190	integer, intent( in) :: bopt
191	real*8 , intent( in) :: delx(:)
192	real*8 , intent( in) :: fdat(:,:,:)
193	real*8 , intent(out) :: edge(:,:)
194	real*8 , intent(out) :: dfdx(:,:)
195	real*8 , intent( in) :: dmin
196	type(rcon_ends), intent(in) :: bcon(:)
197
198	!------------------------------------------- variables !
199	integer :: ivar,idof,isel, &
200	& head,tail,nsel
201	logical :: okay
202	real*8 :: xhat
203	real*8 :: delh(-1:+1)
204	real*8 :: xmap(-1:+2)
205	real*8 :: bvec(+3,-1:+2)
206	real*8 :: gvec(+3,-1:+2)
207	real*8 :: cmat(+3,+3)
208	real*8 :: fhat(+3, nvar)
209	real*8 :: eval(-1:+2)
210	real*8 :: gval(-1:+2)
211
212	integer, parameter :: NSIZ = +3
213	real*8 , parameter :: ZERO = +1.d-14
214
215	head = +2; tail = npos - 2
216
217	if (size(delx).gt.+1) then
218
219	!------------------ mean grid spacing about ii-th cell !
220
221	xhat = max(delx(head),dmin) * 0.5d+0
222
223	!------------------ grid spacing for all stencil cells !
224
225	delh(-1) = delx(head-1)
226	delh(+0) = delx(head+0)
227	delh(+1) = delx(head+1)
228
229	else
230
231	!------------------ mean grid spacing about ii-th cell !
232
233	xhat = max(delx( +1),dmin) * 0.5d+0
234
235	!------------------ grid spacing for all stencil cells !
236
237	delh(-1) = delx( +1)
238	delh(+0) = delx( +1)
239	delh(+1) = delx( +1)
240
241	end if
242
243	!---------- local coordinate mapping for stencil edges !
244
245	xmap(-1) =-(delh(-1) + &
246	& delh(+0)*0.5d0)/xhat
247	xmap(+0) = -1.d0
248	xmap(+1) = +1.d0
249	xmap(+2) = (delh(+1) + &
250	& delh(+0)*0.5d0)/xhat
251
252	!------------ linear system: lhs reconstruction matrix !
253
254	select case(bopt )
255	case( bcon_loose )
256
257	call bfun1d(-1,+3,xmap(-1),bvec(:,-1))
258	call bfun1d(-1,+3,xmap(+0),bvec(:,+0))
259	call bfun1d(-1,+3,xmap(+1),bvec(:,+1))
260	call bfun1d(-1,+3,xmap(+2),bvec(:,+2))
261
262	do idof = +1 , +3
263
264	cmat(1,idof) = bvec(idof,+0) &
265	& - bvec(idof,-1)
266	cmat(2,idof) = bvec(idof,+1) &
267	& - bvec(idof,+0)
268	cmat(3,idof) = bvec(idof,+2) &
269	& - bvec(idof,+1)
270
271	end do
272
273	case( bcon_value )
274
275	call bfun1d(+0,+3,xmap(-1),gvec(:,-1))
276
277	call bfun1d(-1,+3,xmap(-1),bvec(:,-1))
278	call bfun1d(-1,+3,xmap(+0),bvec(:,+0))
279	call bfun1d(-1,+3,xmap(+1),bvec(:,+1))
280
281	do idof = +1 , +3
282
283	cmat(1,idof) = bvec(idof,+0) &
284	& - bvec(idof,-1)
285	cmat(2,idof) = bvec(idof,+1) &
286	& - bvec(idof,+0)
287
288	cmat(3,idof) = gvec(idof,-1)
289
290	end do
291
292	case( bcon_slope )
293
294	call bfun1d(+1,+3,xmap(-1),gvec(:,-1))
295
296	call bfun1d(-1,+3,xmap(-1),bvec(:,-1))
297	call bfun1d(-1,+3,xmap(+0),bvec(:,+0))
298	call bfun1d(-1,+3,xmap(+1),bvec(:,+1))
299
300	do idof = +1 , +3
301
302	cmat(1,idof) = bvec(idof,+0) &
303	& - bvec(idof,-1)
304	cmat(2,idof) = bvec(idof,+1) &
305	& - bvec(idof,+0)
306
307	cmat(3,idof) = gvec(idof,-1)
308
309	end do
310
311	end select
312
313	!------------ linear system: rhs reconstruction vector !
314
315	isel = 0 ; nsel = 0
316
317	select case( bopt )
318	case ( bcon_loose )
319
320	do ivar = +1, nvar
321
322	if (bcon(ivar)%bcopt.eq.bcon_loose) then
323
324	isel = isel + 1
325	nsel = nsel + 1
326
327	fhat(1,isel) = delh(-1) * &
328	& fdat(1,ivar,head-1) / xhat
329	fhat(2,isel) = delh(+0) * &
330	& fdat(1,ivar,head+0) / xhat
331	fhat(3,isel) = delh(+1) * &
332	& fdat(1,ivar,head+1) / xhat
333
334	end if
335
336	end do
337
338	case ( bcon_value )
339
340	do ivar = +1, nvar
341
342	if (bcon(ivar)%bcopt.eq.bcon_value) then
343
344	isel = isel + 1
345	nsel = nsel + 1
346
347	fhat(1,isel) = delh(-1) * &
348	& fdat(1,ivar,head-1) / xhat
349	fhat(2,isel) = delh(+0) * &
350	& fdat(1,ivar,head+0) / xhat
351
352	fhat(3,isel) = bcon(ivar)%value
353
354	end if
355
356	end do
357
358	case ( bcon_slope )
359
360	do ivar = +1, nvar
361
362	if (bcon(ivar)%bcopt.eq.bcon_slope) then
363
364	isel = isel + 1
365	nsel = nsel + 1
366
367	fhat(1,isel) = delh(-1) * &
368	& fdat(1,ivar,head-1) / xhat
369	fhat(2,isel) = delh(+0) * &
370	& fdat(1,ivar,head+0) / xhat
371
372	fhat(3,isel) = &
373	& bcon(ivar)%slope * xhat
374
375	end if
376
377	end do
378
379	end select
380
381	!------------------------- factor/solve linear systems !
382
383	call slv_3x3(cmat,NSIZ,fhat , &
384	& NSIZ,nvar, &
385	& ZERO*dmin,okay)
386
387	if (okay .eqv..false.) then
388
389	# ifdef __PPR_WARNMAT__
390
391	write(,) &
392	& "WARNING::LBC-matrix-is-singular!"
393
394	# endif
395
396	end if
397
398	if (okay .eqv. .true.) then
399
400	!------------- extrapolate values/slopes at lower edge !
401
402	isel = +0
403
404	call bfun1d(+0,+3,xmap(-1),bvec(:,-1))
405	call bfun1d(+0,+3,xmap(+0),bvec(:,+0))
406	call bfun1d(+0,+3,xmap(+1),bvec(:,+1))
407
408	call bfun1d(+1,+3,xmap(-1),gvec(:,-1))
409	call bfun1d(+1,+3,xmap(+0),gvec(:,+0))
410	call bfun1d(+1,+3,xmap(+1),gvec(:,+1))
411
412	do ivar = +1, nvar
413
414	if (bcon(ivar)%bcopt.eq.bopt) then
415
416	isel = isel + 1
417
418	eval(-1) = dot_product( &
419	& bvec(:,-1),fhat(:,isel))
420	eval(+0) = dot_product( &
421	& bvec(:,+0),fhat(:,isel))
422	eval(+1) = dot_product( &
423	& bvec(:,+1),fhat(:,isel))
424
425	gval(-1) = dot_product( &
426	& gvec(:,-1),fhat(:,isel))
427	gval(+0) = dot_product( &
428	& gvec(:,+0),fhat(:,isel))
429	gval(+1) = dot_product( &
430	& gvec(:,+1),fhat(:,isel))
431
432	edge(ivar,head-1) = eval(-1)
433	edge(ivar,head+0) = eval(+0)
434	edge(ivar,head+1) = eval(+1)
435
436	dfdx(ivar,head-1) = gval(-1) &
437	& / xhat
438	dfdx(ivar,head+0) = gval(+0) &
439	& / xhat
440	dfdx(ivar,head+1) = gval(+1) &
441	& / xhat
442
443	end if
444
445	end do
446
447	else
448
449	!------------- low-order if re-con. matrix is singular !
450
451	do ivar = +1, nvar
452
453	if (bcon(ivar)%bcopt.eq.bopt) then
454
455	eval(-1) = &
456	& fdat(1,ivar,head-1) * 1.d0
457	eval(+0) = &
458	& fdat(1,ivar,head-1) * .5d0 + &
459	& fdat(1,ivar,head+0) * .5d0
460	eval(+1) = &
461	& fdat(1,ivar,head+0) * .5d0 + &
462	& fdat(1,ivar,head+1) * .5d0
463
464	gval(-1) = &
465	& fdat(1,ivar,head+0) * .5d0 - &
466	& fdat(1,ivar,head-1) * .5d0
467	gval(+0) = &
468	& fdat(1,ivar,head+0) * .5d0 - &
469	& fdat(1,ivar,head-1) * .5d0
470	gval(+1) = &
471	& fdat(1,ivar,head+1) * .5d0 - &
472	& fdat(1,ivar,head+0) * .5d0
473
474	edge(ivar,head-1) = eval(-1)
475	edge(ivar,head+0) = eval(+0)
476	edge(ivar,head+1) = eval(+1)
477
478	dfdx(ivar,head-1) = gval(-1) &
479	& / xhat
480	dfdx(ivar,head+0) = gval(+0) &
481	& / xhat
482	dfdx(ivar,head+1) = gval(+1) &
483	& / xhat
484
485	end if
486
487	end do
488
489	end if
490
491	return
492
493	end subroutine
494
495	! UBC: impose a single B.C.-type at the upper endpoint !
496
497	subroutine ubc(npos,nvar,ndof,delx, &
498	& fdat,bcon,bopt,edge, &
499	& dfdx,dmin)
500
501	!
502	! NPOS no. edges over grid.
503	! NVAR no. state variables.
504	! NDOF no. degrees-of-freedom per grid-cell.
505	! DELX grid-cell spacing array. LENGTH(DELX) == +1 if
506	! spacing is uniform .
507	! FDAT grid-cell moments array. FDAT is an array with
508	! SIZE = NDOF-by-NVAR-by-NPOS-1 .
509	! BCON boundary condition data for endpoint .
510	! EDGE edge-centred interp. for function-value. EDGE
511	! is an array with SIZE = NVAR-by-NPOS .
512	! DFDX edge-centred interp. for 1st-derivative. DFDX
513	! is an array with SIZE = NVAR-by-NPOS .
514	! DMIN min. grid-cell spacing thresh .
515	!
516
517	implicit none
518
519	!------------------------------------------- arguments !
520	integer, intent( in) :: npos,nvar,ndof
521	integer, intent( in) :: bopt
522	real*8 , intent( in) :: delx(:)
523	real*8 , intent( in) :: fdat(:,:,:)
524	real*8 , intent(out) :: edge(:,:)
525	real*8 , intent(out) :: dfdx(:,:)
526	real*8 , intent( in) :: dmin
527	type(rcon_ends), intent(in) :: bcon(:)
528
529	!------------------------------------------- variables !
530	integer :: ivar,idof,isel, &
531	& head,tail,nsel
532	logical :: okay
533	real*8 :: xhat
534	real*8 :: delh(-1:+1)
535	real*8 :: xmap(-1:+2)
536	real*8 :: bvec(+3,-1:+2)
537	real*8 :: gvec(+3,-1:+2)
538	real*8 :: cmat(+3,+3)
539	real*8 :: fhat(+3, nvar)
540	real*8 :: eval(-1:+2)
541	real*8 :: gval(-1:+2)
542
543	integer, parameter :: NSIZ = +3
544	real*8 , parameter :: ZERO = +1.d-14
545
546	head = +2; tail = npos - 2
547
548	if (size(delx).gt.+1) then
549
550	!------------------ mean grid spacing about ii-th cell !
551
552	xhat = max(delx(tail),dmin) * 0.5d+0
553
554	!------------------ grid spacing for all stencil cells !
555
556	delh(-1) = delx(tail-1)
557	delh(+0) = delx(tail+0)
558	delh(+1) = delx(tail+1)
559
560	else
561
562	!------------------ mean grid spacing about ii-th cell !
563
564	xhat = max(delx( +1),dmin) * 0.5d+0
565
566	!------------------ grid spacing for all stencil cells !
567
568	delh(-1) = delx( +1)
569	delh(+0) = delx( +1)
570	delh(+1) = delx( +1)
571
572	end if
573
574	!---------- local coordinate mapping for stencil edges !
575
576	xmap(-1) =-(delh(-1) + &
577	& delh(+0)*0.5d0)/xhat
578	xmap(+0) = -1.d0
579	xmap(+1) = +1.d0
580	xmap(+2) = (delh(+1) + &
581	& delh(+0)*0.5d0)/xhat
582
583	!------------ linear system: lhs reconstruction matrix !
584
585	select case(bopt )
586	case( bcon_loose )
587
588	call bfun1d(-1,+3,xmap(-1),bvec(:,-1))
589	call bfun1d(-1,+3,xmap(+0),bvec(:,+0))
590	call bfun1d(-1,+3,xmap(+1),bvec(:,+1))
591	call bfun1d(-1,+3,xmap(+2),bvec(:,+2))
592
593	do idof = +1 , +3
594
595	cmat(1,idof) = bvec(idof,+0) &
596	& - bvec(idof,-1)
597	cmat(2,idof) = bvec(idof,+1) &
598	& - bvec(idof,+0)
599	cmat(3,idof) = bvec(idof,+2) &
600	& - bvec(idof,+1)
601
602	end do
603
604	case( bcon_value )
605
606	call bfun1d(+0,+3,xmap(+2),gvec(:,+2))
607
608	call bfun1d(-1,+3,xmap(+0),bvec(:,+0))
609	call bfun1d(-1,+3,xmap(+1),bvec(:,+1))
610	call bfun1d(-1,+3,xmap(+2),bvec(:,+2))
611
612	do idof = +1 , +3
613
614	cmat(1,idof) = bvec(idof,+1) &
615	& - bvec(idof,+0)
616	cmat(2,idof) = bvec(idof,+2) &
617	& - bvec(idof,+1)
618
619	cmat(3,idof) = gvec(idof,+2)
620
621	end do
622
623	case( bcon_slope )
624
625	call bfun1d(+1,+3,xmap(+2),gvec(:,+2))
626
627	call bfun1d(-1,+3,xmap(+0),bvec(:,+0))
628	call bfun1d(-1,+3,xmap(+1),bvec(:,+1))
629	call bfun1d(-1,+3,xmap(+2),bvec(:,+2))
630
631	do idof = +1 , +3
632
633	cmat(1,idof) = bvec(idof,+1) &
634	& - bvec(idof,+0)
635	cmat(2,idof) = bvec(idof,+2) &
636	& - bvec(idof,+1)
637
638	cmat(3,idof) = gvec(idof,+2)
639
640	end do
641
642	end select
643
644	!------------ linear system: rhs reconstruction vector !
645
646	isel = 0 ; nsel = 0
647
648	select case( bopt )
649	case ( bcon_loose )
650
651	do ivar = +1, nvar
652
653	if (bcon(ivar)%bcopt.eq.bcon_loose) then
654
655	isel = isel + 1
656	nsel = nsel + 1
657
658	fhat(1,isel) = delh(-1) * &
659	& fdat(1,ivar,tail-1) / xhat
660	fhat(2,isel) = delh(+0) * &
661	& fdat(1,ivar,tail+0) / xhat
662	fhat(3,isel) = delh(+1) * &
663	& fdat(1,ivar,tail+1) / xhat
664
665	end if
666
667	end do
668
669	case ( bcon_value )
670
671	do ivar = +1, nvar
672
673	if (bcon(ivar)%bcopt.eq.bcon_value) then
674
675	isel = isel + 1
676	nsel = nsel + 1
677
678	fhat(1,isel) = delh(+0) * &
679	& fdat(1,ivar,tail+0) / xhat
680	fhat(2,isel) = delh(+1) * &
681	& fdat(1,ivar,tail+1) / xhat
682
683	fhat(3,isel) = bcon(ivar)%value
684
685	end if
686
687	end do
688
689	case ( bcon_slope )
690
691	do ivar = +1, nvar
692
693	if (bcon(ivar)%bcopt.eq.bcon_slope) then
694
695	isel = isel + 1
696	nsel = nsel + 1
697
698	fhat(1,isel) = delh(+0) * &
699	& fdat(1,ivar,tail+0) / xhat
700	fhat(2,isel) = delh(+1) * &
701	& fdat(1,ivar,tail+1) / xhat
702
703	fhat(3,isel) = &
704	& bcon(ivar)%slope * xhat
705
706	end if
707
708	end do
709
710	end select
711
712	!------------------------- factor/solve linear systems !
713
714	call slv_3x3(cmat,NSIZ,fhat , &
715	& NSIZ,nvar, &
716	& ZERO*dmin,okay)
717
718	if (okay .eqv..false.) then
719
720	# ifdef __PPR_WARNMAT__
721
722	write(,) &
723	& "WARNING::UBC-matrix-is-singular!"
724
725	# endif
726
727	end if
728
729	if (okay .eqv. .true.) then
730
731	!------------- extrapolate values/slopes at lower edge !
732
733	isel = +0
734
735	call bfun1d(+0,+3,xmap(+0),bvec(:,+0))
736	call bfun1d(+0,+3,xmap(+1),bvec(:,+1))
737	call bfun1d(+0,+3,xmap(+2),bvec(:,+2))
738
739	call bfun1d(+1,+3,xmap(+0),gvec(:,+0))
740	call bfun1d(+1,+3,xmap(+1),gvec(:,+1))
741	call bfun1d(+1,+3,xmap(+2),gvec(:,+2))
742
743	do ivar = +1, nvar
744
745	if (bcon(ivar)%bcopt.eq.bopt) then
746
747	isel = isel + 1
748
749	eval(+0) = dot_product( &
750	& bvec(:,+0),fhat(:,isel))
751	eval(+1) = dot_product( &
752	& bvec(:,+1),fhat(:,isel))
753	eval(+2) = dot_product( &
754	& bvec(:,+2),fhat(:,isel))
755
756	gval(+0) = dot_product( &
757	& gvec(:,+0),fhat(:,isel))
758	gval(+1) = dot_product( &
759	& gvec(:,+1),fhat(:,isel))
760	gval(+2) = dot_product( &
761	& gvec(:,+2),fhat(:,isel))
762
763	edge(ivar,tail+0) = eval(+0)
764	edge(ivar,tail+1) = eval(+1)
765	edge(ivar,tail+2) = eval(+2)
766
767	dfdx(ivar,tail+0) = gval(+0) &
768	& / xhat
769	dfdx(ivar,tail+1) = gval(+1) &
770	& / xhat
771	dfdx(ivar,tail+2) = gval(+2) &
772	& / xhat
773
774	end if
775
776	end do
777
778	else
779
780	!------------- low-order if re-con. matrix is singular !
781
782	do ivar = +1, nvar
783
784	if (bcon(ivar)%bcopt.eq.bopt) then
785
786	eval(+0) = &
787	& fdat(1,ivar,tail-1) * .5d0 + &
788	& fdat(1,ivar,tail+0) * .5d0
789	eval(+1) = &
790	& fdat(1,ivar,tail+0) * .5d0 + &
791	& fdat(1,ivar,tail+1) * .5d0
792	eval(+2) = &
793	& fdat(1,ivar,tail+1) * 1.d0
794
795	gval(+0) = &
796	& fdat(1,ivar,tail+0) * .5d0 - &
797	& fdat(1,ivar,tail-1) * .5d0
798	gval(+1) = &
799	& fdat(1,ivar,tail+1) * .5d0 - &
800	& fdat(1,ivar,tail+0) * .5d0
801	gval(+2) = &
802	& fdat(1,ivar,tail+1) * .5d0 - &
803	& fdat(1,ivar,tail+0) * .5d0
804
805	edge(ivar,tail+0) = eval(+0)
806	edge(ivar,tail+1) = eval(+1)
807	edge(ivar,tail+2) = eval(+2)
808
809	dfdx(ivar,tail+0) = gval(+0) &
810	& / xhat
811	dfdx(ivar,tail+1) = gval(+1) &
812	& / xhat
813	dfdx(ivar,tail+2) = gval(+2) &
814	& / xhat
815
816	end if
817
818	end do
819
820	end if
821
822	return
823
824	end subroutine
825
826
827

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

Download in other formats:

Original Format