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.
oscl1d.h90 in vendors/PPR/src – NEMO

source: vendors/PPR/src/oscl1d.h90 @ 14387

Last change on this file since 14387 was 14387, checked in by smasson, 3 years ago

PPR: suppress real*8 <-> real*16 compilation errors or implicit conversions, see #2617
File size: 8.3 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    ! OSCL1D.h90: "oscillation-indicators" for WENO interp.
31    !
32    ! Darren Engwirda
33    ! 08-Sep-2016
34    ! de2363 [at] columbia [dot] edu
35    !
36    !
37   
38    pure subroutine oscli (npos,nvar,ndof,delx,&
39        &                  fdat,oscl,dmin)
40
41    !
42    ! NPOS  no. edges over grid.
43    ! NVAR  no. state variables.
44    ! NDOF  no. degrees-of-freedom per grid-cell .
45    ! DELX  (constant) grid-cell spacing. LENGTH(DELX)==+1 .
46    ! FDAT  grid-cell moments array. FDAT is an array with
47    !       SIZE = NDOF-by-NVAR-by-NPOS-1 .
48    ! OSCL  grid-cell oscil. dof.'s. OSCL is an array with
49    !       SIZE = +2  -by-NVAR-by-NPOS-1 .
50    ! DMIN  min. grid-cell spacing thresh .
51    !
52
53        implicit none
54
55    !------------------------------------------- arguments !
56        integer, intent( in) :: npos,nvar,ndof
57        real*8 , intent( in) :: dmin
58        real*8 , intent( in) :: delx(:)
59        real*8 , intent( in) :: fdat(:,:,:)
60        real*8 , intent(out) :: oscl(:,:,:)
61
62    !------------------------------------------- variables !
63        integer :: ivar,ipos
64       
65        if (npos.lt.3) then
66    !------------------------------- at least 3 grid-cells !
67        do  ipos = +1, npos-1
68        do  ivar = +1, nvar-0
69            oscl(1,ivar,ipos) = +0.e0
70            oscl(2,ivar,ipos) = +0.e0
71        end do
72        end do
73        end if
74       
75        if (npos.lt.3) return
76        if (nvar.lt.1) return
77        if (ndof.lt.1) return
78
79        if (size(delx).gt.+1) then
80
81    !------------------------------- variable grid-spacing !
82
83            call osclv(npos,nvar,ndof,delx, &
84        &              fdat,oscl,dmin)
85       
86        else
87
88    !------------------------------- constant grid-spacing !
89       
90            call osclc(npos,nvar,ndof,delx, &
91        &              fdat,oscl,dmin)
92               
93        end if
94
95        return
96       
97    end  subroutine
98   
99    pure subroutine osclv (npos,nvar,ndof,delx,&
100        &                  fdat,oscl,dmin)
101
102    !
103    ! *this is the variable grid-spacing variant .
104    !
105    ! NPOS  no. edges over grid.
106    ! NVAR  no. state variables.
107    ! NDOF  no. degrees-of-freedom per grid-cell .
108    ! DELX  (variable) grid-cell spacing. LENGTH(DELX)!=+1 .
109    ! FDAT  grid-cell moments array. FDAT is an array with
110    !       SIZE = NDOF-by-NVAR-by-NPOS-1 .
111    ! OSCL  grid-cell oscil. dof.'s. OSCL is an array with
112    !       SIZE = +2  -by-NVAR-by-NPOS-1 .
113    ! DMIN  min. grid-cell spacing thresh .
114    !
115
116        implicit none
117
118    !------------------------------------------- arguments !
119        integer, intent( in) :: npos,nvar,ndof
120        real*8 , intent( in) :: dmin
121        real*8 , intent( in) :: delx(:)
122        real*8 , intent( in) :: fdat(:,:,:)
123        real*8 , intent(out) :: oscl(:,:,:)
124
125    !------------------------------------------- variables !
126        integer :: head,tail
127        integer :: ipos,ivar
128        real*8  :: hhll,hhcc,hhrr
129        real*8  :: hhmm,hhrc,hhlc
130        real*8  :: cmat(2,3)
131
132        head = +1 ; tail = npos-1
133
134    !--------------------------------------- centred point !
135
136        do  ipos = head+1, tail-1
137
138        hhll = max(delx(ipos-1),dmin)
139        hhcc = max(delx(ipos+0),dmin)
140        hhrr = max(delx(ipos+1),dmin)
141
142        hhrc = hhrr + hhcc
143        hhlc = hhll + hhcc
144        hhmm = hhll + hhcc + hhrr
145
146        cmat(1,1) = -(hhcc+2.e0*hhrr)/(hhlc*hhmm)
147        cmat(1,2) = -(hhll-hhrr)* &
148        &          (3.e0*hhcc+2.e0*(hhll+hhrr))/&
149        &            (hhlc*hhrc*hhmm)
150        cmat(1,3) = +(hhcc+2.e0*hhll)/(hhrc*hhmm)
151
152        cmat(2,1) = +3.e0/(hhlc*hhmm)
153        cmat(2,2) = -3.e0*(2.e0*hhcc+hhll+hhrr)/&
154        &            (hhlc*hhrc*hhmm)
155        cmat(2,3) = +3.e0/(hhrc*hhmm)
156
157        do  ivar = 1, nvar
158
159            oscl(1,ivar,ipos) = +1.e0 * ( &
160        & + cmat(1,1)*fdat(1,ivar,ipos-1) &
161        & + cmat(1,2)*fdat(1,ivar,ipos+0) &
162        & + cmat(1,3)*fdat(1,ivar,ipos+1) )
163
164            oscl(2,ivar,ipos) = +2.e0 * ( &
165        & + cmat(2,1)*fdat(1,ivar,ipos-1) &
166        & + cmat(2,2)*fdat(1,ivar,ipos+0) &
167        & + cmat(2,3)*fdat(1,ivar,ipos+1) )
168
169        end do
170           
171        end do
172
173    !-------------------------------------- lower endpoint !
174
175        hhll = max(delx(head+0),dmin)
176        hhcc = max(delx(head+1),dmin)
177        hhrr = max(delx(head+2),dmin)
178
179        cmat(1,1) = -2.e0 / (hhll+hhcc)
180        cmat(1,2) = +2.e0 / (hhll+hhcc)
181
182        do  ivar = 1, nvar
183
184            oscl(1,ivar,head) = &
185        & + cmat(1,1)*fdat(1,ivar,head+0) &
186        & + cmat(1,2)*fdat(1,ivar,head+1)
187             
188            oscl(2,ivar,head) = +0.e0
189
190        end do
191
192    !-------------------------------------- upper endpoint !
193
194        hhll = max(delx(tail-2),dmin)
195        hhcc = max(delx(tail-1),dmin)
196        hhrr = max(delx(tail-0),dmin)
197
198        cmat(1,2) = -2.e0 / (hhrr+hhcc)
199        cmat(1,3) = +2.e0 / (hhrr+hhcc)
200
201        do  ivar = 1, nvar
202
203            oscl(1,ivar,tail) = &
204        & + cmat(1,2)*fdat(1,ivar,tail-1) &
205        & + cmat(1,3)*fdat(1,ivar,tail+0)
206
207            oscl(2,ivar,tail) = +0.e0
208
209        end do
210       
211        return
212
213    end  subroutine
214   
215    pure subroutine osclc (npos,nvar,ndof,delx,&
216        &                  fdat,oscl,dmin)
217
218    !
219    ! *this is the constant grid-spacing variant .
220    !
221    ! NPOS  no. edges over grid.
222    ! NVAR  no. state variables.
223    ! NDOF  no. degrees-of-freedom per grid-cell .
224    ! DELX  (constant) grid-cell spacing. LENGTH(DELX)==+1 .
225    ! FDAT  grid-cell moments array. FDAT is an array with
226    !       SIZE = NDOF-by-NVAR-by-NPOS-1 .
227    ! OSCL  grid-cell oscil. dof.'s. OSCL is an array with
228    !       SIZE = +2  -by-NVAR-by-NPOS-1 .
229    ! DMIN  min. grid-cell spacing thresh .
230    !
231
232        implicit none
233
234    !------------------------------------------- arguments !
235        integer, intent( in) :: npos,nvar,ndof
236        real*8 , intent( in) :: dmin
237        real*8 , intent( in) :: delx(1)
238        real*8 , intent( in) :: fdat(:,:,:)
239        real*8 , intent(out) :: oscl(:,:,:)
240
241    !------------------------------------------- variables !
242        integer :: head,tail,ipos,ivar
243       
244        head = +1; tail = npos - 1
245
246    !-------------------------------------- centred points !
247
248        do  ipos = 2, npos-2
249        do  ivar = 1, nvar-0
250
251            oscl(1,ivar,ipos) = &
252        & + .25d+0 * fdat(1,ivar,ipos+1) &
253        & - .25d+0 * fdat(1,ivar,ipos-1)
254       
255            oscl(2,ivar,ipos) = &
256        & + .25d+0 * fdat(1,ivar,ipos+1) &
257        & - .50d+0 * fdat(1,ivar,ipos+0) &
258        & + .25d+0 * fdat(1,ivar,ipos-1)
259
260        end do
261        end do
262
263    !-------------------------------------- lower endpoint !
264
265        do  ivar = 1, nvar
266
267            oscl(1,ivar,head) = &
268        & + .50d+0 * fdat(1,ivar,head+1) &
269        & - .50d+0 * fdat(1,ivar,head+0)
270       
271            oscl(2,ivar,head) = +0.e0
272
273        end do
274
275    !-------------------------------------- upper endpoint !
276
277        do  ivar = 1, nvar
278
279            oscl(1,ivar,tail) = &
280        & + .50d+0 * fdat(1,ivar,tail+0) &
281        & - .50d+0 * fdat(1,ivar,tail-1)
282           
283            oscl(2,ivar,tail) = +0.e0
284
285        end do
286       
287        return
288
289    end  subroutine
290   
291   
292
Note: See TracBrowser for help on using the repository browser.