boundary_custom.f90
Go to the documentation of this file.
1 !#############################################################################
2 !# #
3 !# fosite - 3D hydrodynamical simulation program #
4 !# module: boundary_custom.f90 #
5 !# #
6 !# Copyright (C) 2006-2018 #
7 !# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
8 !# Jannes Klee <jklee@astrophysik.uni-kiel.de> #
9 !# #
10 !# This program is free software; you can redistribute it and/or modify #
11 !# it under the terms of the GNU General Public License as published by #
12 !# the Free Software Foundation; either version 2 of the License, or (at #
13 !# your option) any later version. #
14 !# #
15 !# This program is distributed in the hope that it will be useful, but #
16 !# WITHOUT ANY WARRANTY; without even the implied warranty of #
17 !# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
18 !# NON INFRINGEMENT. See the GNU General Public License for more #
19 !# details. #
20 !# #
21 !# You should have received a copy of the GNU General Public License #
22 !# along with this program; if not, write to the Free Software #
23 !# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
24 !# #
25 !#############################################################################
26 
27 !----------------------------------------------------------------------------!
35 !----------------------------------------------------------------------------!
36 MODULE boundary_custom_mod
37  USE boundary_base_mod
38  USE boundary_fixed_mod
39  USE marray_compound_mod
40  USE mesh_base_mod
41  USE physics_base_mod
42  USE common_dict
43  IMPLICIT NONE
44  !--------------------------------------------------------------------------!
45  PRIVATE
46  TYPE, EXTENDS(boundary_fixed) :: boundary_custom
47  INTEGER, DIMENSION(:), ALLOCATABLE :: cbtype
48  REAL, DIMENSION(:,:,:), ALLOCATABLE :: rscale, & !< radial scaling constants
49  invrscale
50  REAL, DIMENSION(:,:,:), POINTER :: radius
51  CONTAINS
52  PROCEDURE :: initboundary_custom
53  PROCEDURE :: finalize
54  PROCEDURE :: setboundarydata
55  PROCEDURE :: setcustomboundaries
56  END TYPE
57  CHARACTER(LEN=32), PARAMETER :: boundcond_name = "custom"
58  !--------------------------------------------------------------------------!
59  ! create bit masks for custom boundary conditions
60  enum, bind(c)
61  enumerator :: custom_undefined = 0, &
62  custom_nograd = 1, & ! no gradients (default)
63  custom_period = 2, & ! periodic
64  custom_reflect = 3, & ! reflecting
65  custom_reflneg = 4, & ! reflect and change sign
66  custom_extrapol = 5, & ! linear extrapolation
67  custom_fixed = 6, & ! set fixed boundary data
68  custom_logexpol = 7, & ! extrapolation of log values
69  custom_outflow = 8, & ! nograd/reflect depending on flow direction
70  custom_kepler = 9, & ! extrapolate according to Keplers law
71  custom_angkepler = 10 ! same but for angular momentum instead of
72  END enum
73  !--------------------------------------------------------------------------!
74  PUBLIC :: &
75  boundary_custom, &
76  ! constants
77  custom_nograd, custom_period, custom_reflect, custom_reflneg, &
78  custom_extrapol, custom_fixed, custom_logexpol, &
79  custom_outflow, custom_kepler, custom_angkepler
80  !--------------------------------------------------------------------------!
81 
82 CONTAINS
83 
85  SUBROUTINE initboundary_custom(this,Mesh,Physics,dir,config)
86  IMPLICIT NONE
87  !------------------------------------------------------------------------!
88  CLASS(boundary_custom), INTENT(INOUT) :: this
89  CLASS(mesh_base), INTENT(IN) :: Mesh
90  CLASS(physics_base), INTENT(IN) :: Physics
91  TYPE(Dict_TYP), POINTER :: config
92  INTEGER, INTENT(IN) :: dir
93  !------------------------------------------------------------------------!
94  INTEGER :: err = 0
95  INTEGER :: i,j,k
96  !------------------------------------------------------------------------!
97  CALL this%InitBoundary(mesh,physics,custom,boundcond_name,dir,config)
98 
99  ! allocate memory for boundary data and mask
100  ALLOCATE(this%cbtype(physics%VNUM),stat=err)
101  IF (err.NE.0) &
102  CALL this%Error("InitBoundary_custom", "Unable to allocate memory.")
103 
104  ! allocate memory for boundary data and mask
105  ! this array contains the boundary condition for each primitive variable;
106  ! the user must call SetCustomBoundaries() after initialization to specifiy
107  ! the actual boundary condition for each variable and supply user defined
108  ! data arrays if necessary
109  this%cbtype(:) = custom_nograd
110  END SUBROUTINE initboundary_custom
111 
112 
114  PURE SUBROUTINE setboundarydata(this,Mesh,Physics,time,pvar)
115  IMPLICIT NONE
116  !------------------------------------------------------------------------!
117  CLASS(boundary_custom), INTENT(INOUT) :: this
118  CLASS(mesh_base), INTENT(IN) :: mesh
119  CLASS(physics_base), INTENT(IN) :: physics
120  REAL, INTENT(IN) :: time
121  CLASS(marray_compound), INTENT(INOUT) :: pvar
122  !------------------------------------------------------------------------!
123  INTEGER :: i,j,k,m
124  !------------------------------------------------------------------------!
125  SELECT CASE(this%GetDirection())
126  CASE(west)
127  DO m=1,physics%VNUM
128  SELECT CASE(this%cbtype(m))
129  CASE(custom_nograd)
130 !NEC$ SHORTLOOP
131  DO i=1,mesh%GINUM
132  pvar%data4d(mesh%IMIN-i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
133  = pvar%data4d(mesh%IMIN,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m)
134  END DO
135  CASE(custom_period)
136 !NEC$ SHORTLOOP
137  DO i=1,mesh%GINUM
138  pvar%data4d(mesh%IMIN-i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
139  = pvar%data4d(mesh%IMAX-i+1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m)
140  END DO
141  CASE(custom_reflect)
142 !NEC$ SHORTLOOP
143  DO i=1,mesh%GINUM
144  pvar%data4d(mesh%IMIN-i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
145  = pvar%data4d(mesh%IMIN+i-1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m)
146  END DO
147  CASE(custom_reflneg)
148 !NEC$ SHORTLOOP
149  DO i=1,mesh%GINUM
150  pvar%data4d(mesh%IMIN-i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
151  = -pvar%data4d(mesh%IMIN+i-1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m)
152  END DO
153  CASE(custom_extrapol)
154 !NEC$ SHORTLOOP
155  DO i=1,mesh%GINUM
156  pvar%data4d(mesh%IMIN-i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
157  = (i+1)*pvar%data4d(mesh%IMIN,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
158  - i*pvar%data4d(mesh%IMIN+1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m)
159  END DO
160  CASE(custom_fixed)
161 !NEC$ SHORTLOOP
162  DO i=1,mesh%GINUM
163  pvar%data4d(mesh%IMIN-i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
164  = this%data(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m)
165  END DO
166  CASE(custom_logexpol)
167 !NEC$ SHORTLOOP
168  DO i=1,mesh%GINUM
169  pvar%data4d(mesh%IMIN-i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
170  = pvar%data4d(mesh%IMIN,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
171  * abs(pvar%data4d(mesh%IMIN,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
172  / pvar%data4d(mesh%IMIN+1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m))**i
173  END DO
174  CASE(custom_outflow)
175 ! REMARK: this would work for any GINUM, but with lower performance
176 ! !NEC$ SHORTLOOP
177 ! DO i=1,Mesh%GINUM
178 ! WHERE (pvar%data4d(Mesh%IMIN,Mesh%JMIN:Mesh%JMAX,Mesh%KMIN:Mesh%KMAX,m).GE.0.0 )
179 ! pvar%data4d(Mesh%IMIN-i,Mesh%JMIN:Mesh%JMAX,Mesh%KMIN:Mesh%KMAX,m) &
180 ! = -pvar%data4d(Mesh%IMIN+i-1,Mesh%JMIN:Mesh%JMAX,Mesh%KMIN:Mesh%KMAX,m)
181 ! ELSEWHERE
182 ! pvar%data4d(Mesh%IMIN-i,Mesh%JMIN:Mesh%JMAX,Mesh%KMIN:Mesh%KMAX,m) &
183 ! = pvar%data4d(Mesh%IMIN,Mesh%JMIN:Mesh%JMAX,Mesh%KMIN:Mesh%KMAX,m)
184 ! END WHERE
185 ! END DO
186  WHERE (pvar%data4d(mesh%IMIN,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m).GE.0.0 )
187  pvar%data4d(mesh%IGMIN,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
188  = -pvar%data4d(mesh%IMIN+1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m)
189  pvar%data4d(mesh%IGMIN+1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
190  = -pvar%data4d(mesh%IMIN,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m)
191  ELSEWHERE
192  pvar%data4d(mesh%IGMIN,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
193  = pvar%data4d(mesh%IMIN,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m)
194  pvar%data4d(mesh%IGMIN+1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
195  = pvar%data4d(mesh%IMIN,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m)
196  END WHERE
197  CASE(custom_kepler)
198 !NEC$ SHORTLOOP
199  DO i=1,mesh%GINUM
200  pvar%data4d(mesh%IMIN-i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
201  = (pvar%data4d(mesh%IMIN,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
202  + mesh%radius%bcenter(mesh%IMIN,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX)*mesh%OMEGA) &
203  * this%invRscale(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX) &
204  - mesh%radius%bcenter(mesh%IMIN-i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX)*mesh%OMEGA
205  END DO
206  CASE(custom_angkepler)
207 !NEC$ SHORTLOOP
208  DO i=1,mesh%GINUM
209  pvar%data4d(mesh%IMIN-i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
210  = (pvar%data4d(mesh%IMIN,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
211  + mesh%radius%bcenter(mesh%IMIN,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX)**2*mesh%OMEGA) &
212  * this%Rscale(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX) &
213  - mesh%radius%bcenter(mesh%IMIN-i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX)**2*mesh%OMEGA
214  END DO
215  CASE DEFAULT
216  this%err = ior(custom,z'0100')
217  END SELECT
218  END DO
219  CASE(east)
220  DO m=1,physics%VNUM
221  SELECT CASE(this%cbtype(m))
222  CASE(custom_nograd)
223 !NEC$ SHORTLOOP
224  DO i=1,mesh%GINUM
225  pvar%data4d(mesh%IMAX+i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
226  = pvar%data4d(mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m)
227  END DO
228  CASE(custom_period)
229 !NEC$ SHORTLOOP
230  DO i=1,mesh%GINUM
231  pvar%data4d(mesh%IMAX+i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
232  = pvar%data4d(mesh%IMIN+i-1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m)
233  END DO
234  CASE(custom_reflect)
235 !NEC$ SHORTLOOP
236  DO i=1,mesh%GINUM
237  pvar%data4d(mesh%IMAX+i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
238  = pvar%data4d(mesh%IMAX-i+1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m)
239  END DO
240  CASE(custom_reflneg)
241 !NEC$ SHORTLOOP
242  DO i=1,mesh%GINUM
243  pvar%data4d(mesh%IMAX+i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
244  = -pvar%data4d(mesh%IMAX-i+1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m)
245  END DO
246  CASE(custom_extrapol)
247 !NEC$ SHORTLOOP
248  DO i=1,mesh%GINUM
249  pvar%data4d(mesh%IMAX+i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
250  = (i+1)*pvar%data4d(mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
251  - i*pvar%data4d(mesh%IMAX-1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m)
252  END DO
253  CASE(custom_fixed)
254 !NEC$ SHORTLOOP
255  DO i=1,mesh%GINUM
256  pvar%data4d(mesh%IMAX+i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
257  = this%data(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m)
258  END DO
259  CASE(custom_logexpol)
260 !NEC$ SHORTLOOP
261  DO i=1,mesh%GINUM
262  pvar%data4d(mesh%IMAX+i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
263  = pvar%data4d(mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
264  * abs(pvar%data4d(mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
265  / pvar%data4d(mesh%IMAX-1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m))**i
266  END DO
267  CASE(custom_outflow)
268 ! REMARK: this would work for any GINUM, but with lower performance
269 ! !NEC$ SHORTLOOP
270 ! DO i=1,Mesh%GINUM
271 ! WHERE (pvar%data4d(Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,Mesh%KMIN:Mesh%KMAX,m).LE.0.0 )
272 ! pvar%data4d(Mesh%IMAX+i,Mesh%JMIN:Mesh%JMAX,Mesh%KMIN:Mesh%KMAX,m) &
273 ! = -pvar%data4d(Mesh%IMAX-i+1,Mesh%JMIN:Mesh%JMAX,Mesh%KMIN:Mesh%KMAX,m)
274 ! ELSEWHERE
275 ! pvar%data4d(Mesh%IMAX+i,Mesh%JMIN:Mesh%JMAX,Mesh%KMIN:Mesh%KMAX,m) &
276 ! = pvar%data4d(Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,Mesh%KMIN:Mesh%KMAX,m)
277 ! END WHERE
278 ! END DO
279  WHERE (pvar%data4d(mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m).LE.0.0 )
280  pvar%data4d(mesh%IGMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
281  = -pvar%data4d(mesh%IMAX-1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m)
282  pvar%data4d(mesh%IGMAX-1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
283  = -pvar%data4d(mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m)
284  ELSEWHERE
285  pvar%data4d(mesh%IGMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
286  = pvar%data4d(mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m)
287  pvar%data4d(mesh%IGMAX-1,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
288  = pvar%data4d(mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m)
289  END WHERE
290  CASE(custom_kepler)
291 !NEC$ SHORTLOOP
292  DO i=1,mesh%GINUM
293  pvar%data4d(mesh%IMAX+i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
294  = (pvar%data4d(mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
295  + mesh%radius%bcenter(mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX)*mesh%OMEGA)&
296  * this%invRscale(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX) &
297  - mesh%radius%bcenter(mesh%IMAX+i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX)*mesh%OMEGA
298  END DO
299  CASE(custom_angkepler)
300 !NEC$ SHORTLOOP
301  DO i=1,mesh%GINUM
302  pvar%data4d(mesh%IMAX+i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
303  = (pvar%data4d(mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
304  + mesh%radius%bcenter(mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX)**2*mesh%OMEGA)&
305  * this%Rscale(i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX) &
306  - mesh%radius%bcenter(mesh%IMAX+i,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX)**2*mesh%OMEGA
307  END DO
308  CASE DEFAULT
309  this%err = ior(custom,z'0200')
310  END SELECT
311  END DO
312  CASE(south)
313  DO m=1,physics%VNUM
314  SELECT CASE(this%cbtype(m))
315  CASE(custom_nograd)
316 !NEC$ SHORTLOOP
317  DO j=1,mesh%GJNUM
318  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX,m) &
319  = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN,mesh%KMIN:mesh%KMAX,m)
320  END DO
321  CASE(custom_period)
322 !NEC$ SHORTLOOP
323  DO j=1,mesh%GJNUM
324  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX,m) &
325  = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX-j+1,mesh%KMIN:mesh%KMAX,m)
326  END DO
327  CASE(custom_reflect)
328 !NEC$ SHORTLOOP
329  DO j=1,mesh%GJNUM
330  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX,m) &
331  = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN+j-1,mesh%KMIN:mesh%KMAX,m)
332  END DO
333  CASE(custom_reflneg)
334 !NEC$ SHORTLOOP
335  DO j=1,mesh%GJNUM
336  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX,m) &
337  = -pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN+j-1,mesh%KMIN:mesh%KMAX,m)
338  END DO
339  CASE(custom_extrapol)
340 !NEC$ SHORTLOOP
341  DO j=1,mesh%GJNUM
342  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX,m) &
343  = (j+1)*pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN,mesh%KMIN:mesh%KMAX,m) &
344  - j*pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN+1,mesh%KMIN:mesh%KMAX,m)
345  END DO
346  CASE(custom_fixed)
347 !NEC$ SHORTLOOP
348  DO j=1,mesh%GJNUM
349  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX,m) &
350  = this%data(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX,m)
351  END DO
352  CASE(custom_logexpol)
353 !NEC$ SHORTLOOP
354  DO j=1,mesh%GJNUM
355  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX,m) &
356  = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN,mesh%KMIN:mesh%KMAX,m) &
357  * abs(pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN,mesh%KMIN:mesh%KMAX,m) &
358  / pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN+1,mesh%KMIN:mesh%KMAX,m))**j
359  END DO
360  CASE(custom_outflow)
361 ! REMARK: this would work for any GJNUM, but with lower performance
362 ! !NEC$ SHORTLOOP
363 ! DO j=1,Mesh%GJNUM
364 ! WHERE (pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN,Mesh%KMIN:Mesh%KMAX,m).GE.0.0 )
365 ! pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN-j,Mesh%KMIN:Mesh%KMAX,m) &
366 ! = -pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN+j-1,Mesh%KMIN:Mesh%KMAX,m)
367 ! ELSEWHERE
368 ! pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN-j,Mesh%KMIN:Mesh%KMAX,m) &
369 ! = pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN,Mesh%KMIN:Mesh%KMAX,m)
370 ! END WHERE
371 ! END DO
372  WHERE (pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN,mesh%KMIN:mesh%KMAX,m).GE.0.0 )
373  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JGMIN,mesh%KMIN:mesh%KMAX,m) &
374  = -pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN+1,mesh%KMIN:mesh%KMAX,m)
375  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JGMIN+1,mesh%KMIN:mesh%KMAX,m) &
376  = -pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN,mesh%KMIN:mesh%KMAX,m)
377  ELSEWHERE
378  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JGMIN,mesh%KMIN:mesh%KMAX,m) &
379  = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN,mesh%KMIN:mesh%KMAX,m)
380  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JGMIN+1,mesh%KMIN:mesh%KMAX,m) &
381  = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN,mesh%KMIN:mesh%KMAX,m)
382  END WHERE
383  CASE(custom_kepler)
384 !NEC$ SHORTLOOP
385  DO j=1,mesh%GJNUM
386  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX,m) &
387  = (pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN,mesh%KMIN:mesh%KMAX,m) &
388  + mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN,mesh%KMIN:mesh%KMAX)*mesh%OMEGA)&
389  * this%invRscale(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX) &
390  - mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX)*mesh%OMEGA
391  END DO
392  CASE(custom_angkepler)
393 !NEC$ SHORTLOOP
394  DO j=1,mesh%GJNUM
395  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX,m) &
396  = (pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN,mesh%KMIN:mesh%KMAX,m) &
397  + mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN,mesh%KMIN:mesh%KMAX)**2*mesh%OMEGA)&
398  * this%Rscale(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX) &
399  - mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN-j,mesh%KMIN:mesh%KMAX)**2*mesh%OMEGA
400  END DO
401  CASE DEFAULT
402  this%err = ior(custom,z'0300')
403  END SELECT
404  END DO
405  CASE(north)
406  DO m=1,physics%VNUM
407  SELECT CASE(this%cbtype(m))
408  CASE(custom_nograd)
409 !NEC$ SHORTLOOP
410  DO j=1,mesh%GJNUM
411  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX,m) &
412  = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX,mesh%KMIN:mesh%KMAX,m)
413  END DO
414  CASE(custom_period)
415 !NEC$ SHORTLOOP
416  DO j=1,mesh%GJNUM
417  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX,m) &
418  = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN+j-1,mesh%KMIN:mesh%KMAX,m)
419  END DO
420  CASE(custom_reflect)
421 !NEC$ SHORTLOOP
422  DO j=1,mesh%GJNUM
423  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX,m) &
424  = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX-j+1,mesh%KMIN:mesh%KMAX,m)
425  END DO
426  CASE(custom_reflneg)
427 !NEC$ SHORTLOOP
428  DO j=1,mesh%GJNUM
429  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX,m) &
430  = -pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX-j+1,mesh%KMIN:mesh%KMAX,m)
431  END DO
432  CASE(custom_extrapol)
433 !NEC$ SHORTLOOP
434  DO j=1,mesh%GJNUM
435  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX,m) &
436  = (j+1)*pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
437  - j*pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX-1,mesh%KMIN:mesh%KMAX,m)
438  END DO
439  CASE(custom_fixed)
440 !NEC$ SHORTLOOP
441  DO j=1,mesh%GJNUM
442  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX,m) &
443  = this%data(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX,m)
444  END DO
445  CASE(custom_logexpol)
446 !NEC$ SHORTLOOP
447  DO j=1,mesh%GJNUM
448  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX,m) &
449  = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
450  * abs(pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
451  / pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX-1,mesh%KMIN:mesh%KMAX,m))**j
452  END DO
453  CASE(custom_outflow)
454 ! REMARK: this would work for any GJNUM, but with lower performance
455 ! !NEC$ SHORTLOOP
456 ! DO j=1,Mesh%GJNUM
457 ! WHERE (pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMAX,Mesh%KMIN:Mesh%KMAX,m).LE.0.0 )
458 ! pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMAX+j,Mesh%KMIN:Mesh%KMAX,m) &
459 ! = -pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMAX-j+1,Mesh%KMIN:Mesh%KMAX,m)
460 ! ELSEWHERE
461 ! pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMAX+j,Mesh%KMIN:Mesh%KMAX,m) &
462 ! = pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMAX,Mesh%KMIN:Mesh%KMAX,m)
463 ! END WHERE
464 ! END DO
465  WHERE (pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX,mesh%KMIN:mesh%KMAX,m).LE.0.0 )
466  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JGMAX,mesh%KMIN:mesh%KMAX,m) &
467  = -pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX-1,mesh%KMIN:mesh%KMAX,m)
468  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JGMAX-1,mesh%KMIN:mesh%KMAX,m) &
469  = -pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX,mesh%KMIN:mesh%KMAX,m)
470  ELSEWHERE
471  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JGMAX,mesh%KMIN:mesh%KMAX,m) &
472  = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX,mesh%KMIN:mesh%KMAX,m)
473  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JGMAX-1,mesh%KMIN:mesh%KMAX,m) &
474  = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX,mesh%KMIN:mesh%KMAX,m)
475  END WHERE
476  CASE(custom_kepler)
477 !NEC$ SHORTLOOP
478  DO j=1,mesh%GJNUM
479  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX,m) &
480  = (pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
481  + mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMAX,mesh%KMIN:mesh%KMAX)*mesh%OMEGA)&
482  * this%invRscale(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX) &
483  - mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX)*mesh%OMEGA
484  END DO
485  CASE(custom_angkepler)
486 !NEC$ SHORTLOOP
487  DO j=1,mesh%GJNUM
488  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX,m) &
489  = (pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMAX,mesh%KMIN:mesh%KMAX,m) &
490  + mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMAX,mesh%KMIN:mesh%KMAX)**2*mesh%OMEGA)&
491  * this%Rscale(mesh%IMIN:mesh%IMAX,j,mesh%KMIN:mesh%KMAX) &
492  - mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMAX+j,mesh%KMIN:mesh%KMAX)**2*mesh%OMEGA
493  END DO
494  CASE DEFAULT
495  this%err = ior(custom,z'0400')
496  END SELECT
497  END DO
498  CASE(bottom)
499  DO m=1,physics%VNUM
500  SELECT CASE(this%cbtype(m))
501  CASE(custom_nograd)
502 !NEC$ SHORTLOOP
503  DO k=1,mesh%GKNUM
504  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN-k,m) &
505  = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN,m)
506  END DO
507  CASE(custom_period)
508 !NEC$ SHORTLOOP
509  DO k=1,mesh%GKNUM
510  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN-k,m) &
511  = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX-k+1,m)
512  END DO
513  CASE(custom_reflect)
514 !NEC$ SHORTLOOP
515  DO k=1,mesh%GKNUM
516  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN-k,m) &
517  = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN+k-1,m)
518  END DO
519  CASE(custom_reflneg)
520 !NEC$ SHORTLOOP
521  DO k=1,mesh%GKNUM
522  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN-k,m) &
523  = -pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN+k-1,m)
524  END DO
525  CASE(custom_extrapol)
526 !NEC$ SHORTLOOP
527  DO k=1,mesh%GKNUM
528  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN-k,m) &
529  = (k+1)*pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN,m) &
530  - k*pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN+1,m)
531  END DO
532  CASE(custom_fixed)
533 !NEC$ SHORTLOOP
534  DO k=1,mesh%GKNUM
535  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN-k,m) &
536  = this%data(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k,m)
537  END DO
538  CASE(custom_logexpol)
539 !NEC$ SHORTLOOP
540  DO k=1,mesh%GKNUM
541  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN-k,m) &
542  = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN,m) &
543  * abs(pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN,m) &
544  / pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN+1,m))**k
545  END DO
546  CASE(custom_outflow)
547 ! REMARK: this would work for any GKNUM, but with lower performance
548 ! !NEC$ SHORTLOOP
549 ! DO k=1,Mesh%GKNUM
550 ! WHERE (pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,Mesh%KMIN,m).GE.0.0 )
551 ! pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,Mesh%KMIN-k,m) &
552 ! = -pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,Mesh%KMIN+k-1,m)
553 ! ELSEWHERE
554 ! pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,Mesh%KMIN-k,m) &
555 ! = pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,Mesh%KMIN,m)
556 ! END WHERE
557 ! END DO
558 !NEC$ SHORTLOOP
559  DO k=1,mesh%GKNUM
560  WHERE (pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN,m).GE.0.0 )
561  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KGMIN,m) &
562  = -pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN+1,m)
563  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KGMIN+1,m) &
564  = -pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN,m)
565  ELSEWHERE
566  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KGMIN,m) &
567  = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN,m)
568  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KGMIN+1,m) &
569  = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN,m)
570  END WHERE
571  END DO
572  CASE(custom_kepler)
573 !NEC$ SHORTLOOP
574  DO k=1,mesh%GKNUM
575  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN-k,m) &
576  = (pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN,m) &
577  + mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN)*mesh%OMEGA)&
578  * this%invRscale(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k) &
579  - mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN-k)*mesh%OMEGA
580  END DO
581  CASE(custom_angkepler)
582 !NEC$ SHORTLOOP
583  DO k=1,mesh%GKNUM
584  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN-k,m) &
585  = (pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN,m) &
586  + mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN)**2*mesh%OMEGA)&
587  * this%invRscale(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k) &
588  - mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN-k)**2*mesh%OMEGA
589  END DO
590  CASE DEFAULT
591  this%err = ior(custom,z'0500')
592  END SELECT
593  END DO
594  CASE(top)
595  DO m=1,physics%VNUM
596  SELECT CASE(this%cbtype(m))
597  CASE(custom_nograd)
598 !NEC$ SHORTLOOP
599  DO k=1,mesh%GKNUM
600  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX+k,m) &
601  = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX,m)
602  END DO
603  CASE(custom_period)
604 !NEC$ SHORTLOOP
605  DO k=1,mesh%GKNUM
606  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX+k,m) &
607  = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN+k-1,m)
608  END DO
609  CASE(custom_reflect)
610 !NEC$ SHORTLOOP
611  DO k=1,mesh%GKNUM
612  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX+k,m) &
613  = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX-k+1,m)
614  END DO
615  CASE(custom_reflneg)
616 !NEC$ SHORTLOOP
617  DO k=1,mesh%GKNUM
618  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX+k,m) &
619  = -pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX-k+1,m)
620  END DO
621  CASE(custom_extrapol)
622 !NEC$ SHORTLOOP
623  DO k=1,mesh%GKNUM
624  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX+k,m) &
625  = (k+1)*pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX,m) &
626  - k*pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX-1,m)
627  END DO
628  CASE(custom_fixed)
629 !NEC$ SHORTLOOP
630  DO k=1,mesh%GKNUM
631  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX+k,m) &
632  = this%data(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k,m)
633  END DO
634  CASE(custom_logexpol)
635 !NEC$ SHORTLOOP
636  DO k=1,mesh%GKNUM
637  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX+k,m) &
638  = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX,m) &
639  * abs(pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX,m) &
640  / pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX-1,m))**k
641  END DO
642  CASE(custom_outflow)
643 ! REMARK: this would work for any GKNUM, but with lower performance
644 ! !NEC$ SHORTLOOP
645 ! DO k=1,Mesh%GKNUM
646 ! WHERE (pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,Mesh%KMAX,m).LE.0.0)
647 ! pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,Mesh%KMAX+k,m) &
648 ! = -pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,Mesh%KMAX-k+1,m)
649 ! ELSEWHERE
650 ! pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,Mesh%KMAX+k,m) &
651 ! = pvar%data4d(Mesh%IMIN:Mesh%IMAX,Mesh%JMIN:Mesh%JMAX,Mesh%KMAX,m)
652 ! END WHERE
653 !NEC$ SHORTLOOP
654  DO k=1,mesh%GKNUM
655  WHERE (pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX,m).LE.0.0)
656  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KGMAX,m) &
657  = -pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX-1,m)
658  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KGMAX-1,m) &
659  = -pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX,m)
660  ELSEWHERE
661  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KGMAX,m) &
662  = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX,m)
663  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KGMAX-1,m) &
664  = pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX,m)
665  END WHERE
666  END DO
667  CASE(custom_kepler)
668 !NEC$ SHORTLOOP
669  DO k=1,mesh%GKNUM
670  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX+k,m) &
671  = (pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX,m) &
672  + mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX)*mesh%OMEGA)&
673  * this%invRscale(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k) &
674  - mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX+k)*mesh%OMEGA
675  END DO
676  CASE(custom_angkepler)
677 !NEC$ SHORTLOOP
678  DO k=1,mesh%GKNUM
679  pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX+k,m) &
680  = (pvar%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX,m) &
681  + mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX)**2*mesh%OMEGA)&
682  * this%invRscale(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,k) &
683  - mesh%radius%bcenter(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMAX+k)**2*mesh%OMEGA
684  END DO
685  CASE DEFAULT
686  this%err = ior(custom,z'0600')
687  END SELECT
688  END DO
689  END SELECT
690  END SUBROUTINE setboundarydata
691 
692  SUBROUTINE setcustomboundaries(this,Mesh,Physics,cbtype,kepler_radius)
693  IMPLICIT NONE
694  !------------------------------------------------------------------------!
695  CLASS(boundary_custom), INTENT(INOUT) :: this
696  CLASS(mesh_base), INTENT(IN) :: Mesh
697  CLASS(physics_base), INTENT(IN) :: Physics
698  INTEGER, DIMENSION(1:Physics%VNUM) :: cbtype
699  REAL, DIMENSION(:,:,:), POINTER, OPTIONAL :: kepler_radius
700  !------------------------------------------------------------------------!
701  INTEGER :: i,j,k,err = 0
702  !------------------------------------------------------------------------!
703  this%cbtype(:) = cbtype(:)
704  IF (any(this%cbtype(:).EQ.custom_fixed)) THEN
705  SELECT CASE(this%GetDirection())
706  CASE(west,east)
707  ALLOCATE(this%data(mesh%GINUM,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,physics%VNUM), &
708  stat=err)
709  CASE(south,north)
710  ALLOCATE(this%data(mesh%IMIN:mesh%IMAX,mesh%GJNUM,mesh%KMIN:mesh%KMAX,physics%VNUM), &
711  stat=err)
712  CASE(bottom,top)
713  ALLOCATE(this%data(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%GKNUM,physics%VNUM), &
714  stat=err)
715  END SELECT
716  IF (err.NE.0) CALL this%Error("SetCustomBoundaries","Unable to allocate memory.")
717  this%data(:,:,:,:) = 0.0
718  END IF
719 
720  IF (any(this%cbtype(:).EQ.custom_outflow)) THEN
721  SELECT CASE(this%GetDirection())
722  CASE(west,east)
723  IF (mesh%GINUM.NE.2) CALL this%Error("SetCustomBoundaries","GINUM must be 2 for outflow boundaries")
724  CASE(south,north)
725  IF (mesh%GJNUM.NE.2) CALL this%Error("SetCustomBoundaries","GJNUM must be 2 for outflow boundaries")
726  CASE(bottom,top)
727  IF (mesh%GKNUM.NE.2) CALL this%Error("SetCustomBoundaries","GKNUM must be 2 for outflow boundaries")
728  END SELECT
729  END IF
730 
731  IF (any(this%cbtype(:).EQ.custom_kepler)) THEN
732  SELECT CASE(this%GetDirection())
733  CASE(west,east)
734  ALLOCATE(this%invRscale(mesh%GINUM,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
735  stat=err)
736  CASE(south,north)
737  ALLOCATE(this%invRscale(mesh%IMIN:mesh%IMAX,mesh%GJNUM,mesh%KMIN:mesh%KMAX), &
738  stat=err)
739  CASE(bottom,top)
740  ALLOCATE(this%invRscale(mesh%IMIN:mesh%IMAX,mesh%KGMIN:mesh%JMAX,mesh%GKNUM), &
741  stat=err)
742  END SELECT
743  IF (err.NE.0) CALL this%Error("SetCustomBoundaries","Unable to allocate memory.")
744 
745  IF (PRESENT(kepler_radius)) THEN
746  this%radius => kepler_radius
747  IF (mesh%OMEGA.GT.tiny(mesh%OMEGA)) &
748  CALL this%Warning("SetCustomBoundaries", &
749  "user supplied radius and rotating frame may yield unexpected/wrong results")
750  ELSE
751  this%radius => mesh%radius%bcenter
752  END IF
753 
754  SELECT CASE(this%GetDirection())
755  CASE(west)
756  FORALL (i=1:mesh%GINUM,j=mesh%JMIN:mesh%JMAX,k=mesh%KMIN:mesh%KMAX) &
757  this%invRscale(i,j,k) = this%radius(mesh%IMIN,j,k) / this%radius(mesh%IMIN-i,j,k)
758  CASE(east)
759  FORALL (i=1:mesh%GINUM,j=mesh%JMIN:mesh%JMAX,k=mesh%KMIN:mesh%KMAX) &
760  this%invRscale(i,j,k) = this%radius(mesh%IMAX,j,k) / this%radius(mesh%IMAX+i,j,k)
761  CASE(south)
762  FORALL (i=mesh%IMIN:mesh%IMAX,j=1:mesh%GjNUM,k=mesh%KMIN:mesh%KMAX) &
763  this%invRscale(i,j,k) = this%radius(i,mesh%JMIN,k) / this%radius(i,mesh%JMIN-j,k)
764  CASE(north)
765  FORALL (i=mesh%IMIN:mesh%IMAX,j=1:mesh%GjNUM,k=mesh%KMIN:mesh%KMAX) &
766  this%invRscale(i,j,k) = this%radius(i,mesh%JMAX,k) / this%radius(i,mesh%JMAX+j,k)
767  CASE(bottom)
768  FORALL (i=mesh%IMIN:mesh%IMAX,j=mesh%JMIN:mesh%JMAX,k=1:mesh%GKNUM) &
769  this%invRscale(i,j,k) = this%radius(i,j,mesh%KMIN) / this%radius(i,j,mesh%KMIN-k)
770  CASE(top)
771  FORALL (i=mesh%IMIN:mesh%IMAX,j=mesh%JMIN:mesh%JMAX,k=1:mesh%GKNUM) &
772  this%invRscale(i,j,k) = this%radius(i,j,mesh%KMAX) / this%radius(i,j,mesh%KMAX+k)
773  END SELECT
774  this%invRscale(:,:,:) = sqrt(this%invRscale(:,:,:))
775  END IF
776 
777  IF (any(this%cbtype(:).EQ.custom_angkepler)) THEN
778  SELECT CASE(this%GetDirection())
779  CASE(west,east)
780  ALLOCATE(this%Rscale(mesh%GINUM,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX), &
781  stat=err)
782  CASE(south,north)
783  ALLOCATE(this%Rscale(mesh%IMIN:mesh%IMAX,mesh%GJNUM,mesh%KMIN:mesh%KMAX), &
784  stat=err)
785  CASE(bottom,top)
786  ALLOCATE(this%Rscale(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%GKNUM), &
787  stat=err)
788  END SELECT
789 
790  IF (PRESENT(kepler_radius)) THEN
791  this%radius => kepler_radius
792  IF (mesh%OMEGA.GT.tiny(mesh%OMEGA)) &
793  CALL this%Warning("SetCustomBoundaries", &
794  "user supplied radius and rotating frame may yield unexpected/wrong results")
795  ELSE
796  this%radius => mesh%radius%bcenter
797  END IF
798 
799  IF (err.NE.0) CALL this%Error("SetCustomBoundaries","Unable to allocate memory.")
800  SELECT CASE(this%GetDirection())
801  CASE(west)
802  FORALL (i=1:mesh%GINUM,j=mesh%JMIN:mesh%JMAX,k=mesh%KMIN:mesh%KMAX) &
803  this%Rscale(i,j,k) = this%radius(mesh%IMIN-i,j,k) / this%radius(mesh%IMIN,j,k)
804  CASE(east)
805  FORALL (i=1:mesh%GINUM,j=mesh%JMIN:mesh%JMAX,k=mesh%KMIN:mesh%KMAX) &
806  this%Rscale(i,j,k) = this%radius(mesh%IMAX+i,j,k) / this%radius(mesh%IMAX,j,k)
807  CASE(south)
808  FORALL (i=mesh%IMIN:mesh%IMAX,j=1:mesh%GjNUM,k=mesh%KMIN:mesh%KMAX) &
809  this%Rscale(i,j,k) = this%radius(i,mesh%JMIN-j,k) / this%radius(i,mesh%JMIN,k)
810  CASE(north)
811  FORALL (i=mesh%IMIN:mesh%IMAX,j=1:mesh%GjNUM,k=mesh%KMIN:mesh%KMAX) &
812  this%Rscale(i,j,k) = this%radius(i,mesh%JMAX+j,k) / this%radius(i,mesh%JMAX,k)
813  CASE(bottom)
814  FORALL (i=mesh%IMIN:mesh%IMAX,j=mesh%JMIN:mesh%JMAX,k=1:mesh%GKNUM) &
815  this%Rscale(i,j,k) = this%radius(i,j,mesh%KMIN-k) / this%radius(i,j,mesh%KMIN)
816  CASE(top)
817  FORALL (i=mesh%IMIN:mesh%IMAX,j=mesh%JMIN:mesh%JMAX,k=1:mesh%GKNUM) &
818  this%Rscale(i,j,k) = this%radius(i,j,mesh%KMAX+k) / this%radius(i,j,mesh%KMAX)
819  END SELECT
820  this%Rscale(:,:,:) = sqrt(this%Rscale(:,:,:))
821  END IF
822  END SUBROUTINE
823 
824  SUBROUTINE finalize(this)
825  IMPLICIT NONE
826  !------------------------------------------------------------------------!
827  CLASS(boundary_custom), INTENT(INOUT) :: this
828  !------------------------------------------------------------------------!
829  IF (ALLOCATED(this%data)) DEALLOCATE(this%data)
830  IF (ALLOCATED(this%cbtype)) DEALLOCATE(this%cbtype)
831  IF (ALLOCATED(this%Rscale)) DEALLOCATE(this%Rscale)
832  IF (ALLOCATED(this%invRscale)) DEALLOCATE(this%invRscale)
833 
834  CALL this%Finalize_base()
835  END SUBROUTINE finalize
836 END MODULE boundary_custom_mod
boundary_fixed_mod::boundary_fixed
Definition: boundary_fixed.f90:46
mesh_base
boundary_custom_mod::setboundarydata
pure subroutine setboundarydata(this, Mesh, Physics, time, pvar)
Applies the custom boundary conditions.
Definition: boundary_custom.f90:115
boundary_custom_mod::boundary_custom
Definition: boundary_custom.f90:46
marray_compound_mod
derived class for compound of mesh arrays
Definition: marray_compound.f90:35
physics_base_mod
Basic physics module.
Definition: physics_base.f90:55
boundary_custom_mod::initboundary_custom
subroutine initboundary_custom(this, Mesh, Physics, dir, config)
Constructor for custom boundary conditions.
Definition: boundary_custom.f90:86
common_dict
Dictionary for generic data types.
Definition: common_dict.f90:61
boundary_custom_mod::setcustomboundaries
subroutine setcustomboundaries(this, Mesh, Physics, cbtype, kepler_radius)
Definition: boundary_custom.f90:693
boundary_custom_mod::finalize
subroutine finalize(this)
Definition: boundary_custom.f90:825
boundary_custom_mod::boundcond_name
character(len=32), parameter boundcond_name
Definition: boundary_custom.f90:57
boundary_custom_mod
Boundary module for custom conditions.
Definition: boundary_custom.f90:36
boundary_fixed_mod
Boundary module for fixed in/outflow conditions.
Definition: boundary_fixed.f90:37
boundary_base_mod
Definition: boundary_base.f90:40
Generated by   doxygen 1.8.14