Line data Source code
1 : !MODULE FVM_ANALYTIC_MOD--------------------------------------------CE-for FVM!
2 : ! AUTHOR: CHRISTOPH ERATH, 17.October 2011 !
3 : ! This module contains all analytical terms for fvm !
4 : !-----------------------------------------------------------------------------!
5 : module fvm_analytic_mod
6 : use shr_kind_mod, only: r8=>shr_kind_r8
7 : use control_mod, only : north, south, east, west, neast, nwest, seast, swest
8 : use cam_abortutils, only: endrun
9 :
10 : implicit none
11 : private
12 :
13 : public :: get_high_order_weights_over_areas, compute_reconstruct_matrix
14 : public :: compute_halo_vars, init_flux_orient
15 : public :: I_00, I_10, I_01, I_20, I_02, I_11, gauss_points
16 : public :: F_00, F_10, F_01, F_20, F_02, F_11
17 : public :: create_interpolation_points, compute_basic_coordinate_vars
18 :
19 : CONTAINS
20 :
21 21600 : subroutine compute_basic_coordinate_vars(elem,&
22 21600 : nc,irecons,dalpha,dbeta,vtx_cart,center_cart,area_sphere,spherecentroid)
23 : use coordinate_systems_mod, only: cart2spherical
24 : use element_mod, only: element_t
25 : use coordinate_systems_mod, only: spherical_polar_t
26 :
27 : type (element_t), intent(in ) :: elem
28 : integer, intent(in) :: nc,irecons
29 :
30 : real (kind=r8), intent(out) :: dalpha, dbeta
31 : real (kind=r8), intent(out) :: vtx_cart (4,2,nc,nc)
32 : real (kind=r8), intent(out) :: area_sphere(nc,nc)
33 : real (kind=r8), intent(out) :: spherecentroid(irecons-1,nc,nc)
34 : type (spherical_polar_t), intent(out) :: center_cart(nc,nc) ! Spherical coordinates of fvm grid
35 :
36 : integer :: i,j
37 : real (kind=r8) :: centerx,centery
38 43200 : real (kind=r8) :: acartx(nc+1), acarty(nc+1)
39 :
40 21600 : dalpha=abs(elem%corners(1)%x-elem%corners(2)%x)/nc
41 21600 : dbeta =abs(elem%corners(1)%y-elem%corners(4)%y)/nc
42 :
43 108000 : do i=1,nc+1
44 86400 : acartx(i) = tan(elem%corners(1)%x+(i-1)*dalpha)
45 108000 : acarty(i) = tan(elem%corners(1)%y+(i-1)*dbeta)
46 : end do
47 :
48 86400 : do j=1,nc
49 280800 : do i=1,nc
50 194400 : centerx = tan(elem%corners(1)%x+(i-0.5_r8)*dalpha)
51 194400 : centery = tan(elem%corners(1)%y+(j-0.5_r8)*dbeta)
52 259200 : center_cart(i,j) = cart2spherical(centerx,centery,elem%FaceNum)
53 : enddo
54 : enddo
55 :
56 2224800 : vtx_cart = -9D9
57 86400 : do j=1,nc
58 280800 : do i=1,nc
59 194400 : vtx_cart(1,1,i,j) = acartx(i )
60 194400 : vtx_cart(1,2,i,j) = acarty(j )
61 :
62 194400 : vtx_cart(2,1,i,j) = acartx(i+1)
63 194400 : vtx_cart(2,2,i,j) = acarty(j )
64 :
65 194400 : vtx_cart(3,1,i,j) = acartx(i+1)
66 194400 : vtx_cart(3,2,i,j) = acarty(j+1)
67 :
68 194400 : vtx_cart(4,1,i,j) = acartx(i )
69 259200 : vtx_cart(4,2,i,j) = acarty(j+1)
70 : end do
71 : end do
72 : ! compute area and centroid for the interior and halo zone of interior elements
73 21600 : call moment_onsphere(nc,irecons,area_sphere,vtx_cart,.true.,spherecentroid)
74 21600 : end subroutine compute_basic_coordinate_vars
75 :
76 21600 : subroutine compute_halo_vars(faceno,cubeboundary,nc,nhc,nhe,&
77 21600 : jx_min,jx_max,jy_min,jy_max,flux_orient, ifct, rot_matrix)
78 : use control_mod, only : north, south, east, west, neast, nwest, seast, swest
79 :
80 : integer, intent(in) :: faceno,nc,nhc,nhe,cubeboundary
81 :
82 : integer, intent(out) :: jx_min(3),jx_max(3),jy_min(3),jy_max(3)
83 : real (kind=r8), intent(out) :: flux_orient(2, 1-nhc:nc+nhc,1-nhc:nc+nhc)
84 : integer, intent(out) :: ifct (1-nhc:nc+nhc,1-nhc:nc+nhc)
85 : integer, intent(out) :: rot_matrix(2,2,1-nhc:nc+nhc,1-nhc:nc+nhc)
86 :
87 : integer :: i,j
88 : integer :: rot90_matrix(2,2)
89 : integer :: ishft
90 :
91 :
92 21600 : jx_min(2) = 0; jx_max(2) = -1; jy_min(2) = 0; jy_max(2) = -1
93 21600 : jx_min(3) = 0; jx_max(3) = -1; jy_min(3) = 0; jy_max(3) = -1
94 :
95 40416 : select case (cubeboundary)
96 : case (0)
97 18816 : jx_min(1)=1-nhe; jx_max(1)=nc+1+nhe; jy_min(1)=1-nhe; jy_max(1)=nc+1+nhe
98 : case (west)
99 672 : jx_min(1)=1 ; jx_max(1)=nc+1+nhe; jy_min(1)=1-nhe; jy_max(1)=nc+1+nhe
100 672 : jx_min(2)=1-nhe; jx_max(2)=1 ; jy_min(2)=1-nhe; jy_max(2)=nc+1+nhe
101 : case(east)
102 672 : jx_min(1)=1-nhe; jx_max(1)=nc+1 ; jy_min(1)=1-nhe; jy_max(1)=nc+1+nhe
103 672 : jx_min(2)=nc+1 ; jx_max(2)=nc+1+nhe; jy_min(2)=1-nhe; jy_max(2)=nc+1+nhe
104 : case(north)
105 672 : jx_min(1)=1-nhe; jx_max(1)=nc+1+nhe; jy_min(1)=1-nhe; jy_max(1)=nc+1
106 672 : jx_min(2)=1-nhe; jx_max(2)=nc+1+nhe; jy_min(2)=nc+1 ; jy_max(2)=nc+1+nhe
107 : case(south)
108 672 : jx_min(1)=1-nhe; jx_max(1)=nc+1+nhe; jy_min(1)=1 ; jy_max(1)=nc+1+nhe
109 672 : jx_min(2)=1-nhe; jx_max(2)=nc+1+nhe; jy_min(2)=1-nhe; jy_max(2)=1
110 : case(swest)
111 24 : jx_min(1)=1 ; jx_max(1)=nc+1+nhe; jy_min(1)=1 ; jy_max(1)=nc+1+nhe
112 24 : jx_min(2)=1 ; jx_max(2)=nc+1+nhe; jy_min(2)=1-nhe; jy_max(2)=1
113 24 : jx_min(3)=1-nhe; jx_max(3)=1 ; jy_min(3)=1 ; jy_max(3)=nc+1+nhe
114 : case(seast)
115 24 : jx_min(1)=1-nhe; jx_max(1)=nc+1 ; jy_min(1)=1 ; jy_max(1)=nc+1+nhe
116 24 : jx_min(2)=1-nhe; jx_max(2)=nc+1 ; jy_min(2)=1-nhe; jy_max(2)=1
117 24 : jx_min(3)=nc+1 ; jx_max(3)=nc+1+nhe; jy_min(3)=1 ; jy_max(3)=nc+1+nhe
118 : case(neast)
119 24 : jx_min(1)=1-nhe; jx_max(1)=nc+1 ; jy_min(1)=1-nhe; jy_max(1)=nc+1
120 24 : jx_min(2)=1-nhe; jx_max(2)=nc+1 ; jy_min(2)=nc+1 ; jy_max(2)=nc+1+nhe
121 24 : jx_min(3)=nc+1 ; jx_max(3)=nc+1+nhe; jy_min(3)=1-nhe; jy_max(3)=nc+1
122 : case(nwest)
123 24 : jx_min(1)=1 ; jx_max(1)=nc+1+nhe; jy_min(1)=1-nhe; jy_max(1)=nc+1
124 24 : jx_min(2)=1 ; jx_max(2)=nc+1+nhe; jy_min(2)=nc+1 ; jy_max(2)=nc+1+nhe
125 24 : jx_min(3)=1-nhe; jx_max(3)=1 ; jy_min(3)=1-nhe; jy_max(3)=nc+1
126 :
127 : case default
128 0 : print *, 'Fatal Error in fvm_line_integrals_mod.F90.'
129 21600 : call endrun('Selected case for cubeboundary does not exists!')
130 : end select
131 : !
132 : ! init location of flux-sides
133 : !
134 21600 : call init_flux_orient(flux_orient,ifct,nc,nhc,cubeboundary,faceno)
135 3909600 : rot_matrix(1,1,:,:) = 1; rot_matrix(1,2,:,:) = 0;
136 3909600 : rot_matrix(2,1,:,:) = 0; rot_matrix(2,2,:,:) = 1;
137 :
138 21600 : if (cubeboundary>0) then
139 : !
140 : ! clockwise 90 rotation of vectors
141 : !
142 2784 : rot90_matrix(1,1) = 0; rot90_matrix(2,1) = -1;
143 2784 : rot90_matrix(1,2) = 1; rot90_matrix(2,2) = 0;
144 27840 : do j=1-nhc,nc+nhc
145 253344 : do i=1-nhc,nc+nhc
146 1076544 : do ishft=1,4-nint(flux_orient(2,i,j))
147 6007392 : rot_matrix(:,:,i,j) = MATMUL(rot90_matrix,rot_matrix(:,:,i,j))
148 : end do
149 : enddo
150 : enddo
151 : end if
152 21600 : end subroutine compute_halo_vars
153 :
154 :
155 : ! ----------------------------------------------------------------------------------!
156 : !SUBROUTINE MOMENT_ONSPHERE-----------------------------------------------CE-for FVM!
157 : ! AUTHOR: CHRISTOPH ERATH, 20.July 2011 !
158 : ! DESCRIPTION: Compute area and centroids/moments via line integrals !
159 : ! !
160 : ! INPUT: x ... x cartesian coordinats of the arrival grid on the cube !
161 : ! y ... y cartesian coordinats of the arrival grid on the cube !
162 : ! ... cell boundaries in x and y directions !
163 : ! INPUT/OUTPUT: !
164 : ! area ... area of cells on the sphere !
165 : ! centroid ... x,y,x^2,y^2,xy !
166 : !-----------------------------------------------------------------------------------!
167 21600 : subroutine moment_onsphere(nc,irecons,area,vtx_cart,lanalytic,spherecentroid)
168 : use dimensions_mod, only: ngpc
169 :
170 : integer, intent(in) :: nc,irecons
171 : real (kind=r8), dimension(nc,nc) , intent(out) :: area
172 : real (kind=r8), dimension(irecons-1,nc,nc), intent(out) :: spherecentroid
173 : real (kind=r8), dimension(4,2,nc,nc) , intent(in) :: vtx_cart
174 : logical, optional, intent(in) :: lanalytic
175 : integer :: i,j
176 : !
177 : ! variables for call to get_high_order_weights_over_areas
178 : !
179 : integer, parameter :: num_area=1, num_seg_max=2
180 : REAL(KIND=r8), dimension(2,num_seg_max,num_area) :: xx, dxx
181 : integer , dimension(num_area ), parameter :: num_seg=2
182 43200 : REAL(KIND=r8), dimension(irecons,num_area):: weights
183 43200 : real (kind=r8), dimension(nc+1) :: x, y
184 :
185 :
186 : real (kind=r8), dimension(ngpc):: gsweights, gspts
187 : !
188 : ! initialize quadrature weights for get_high_order_weights_over_areas
189 : !
190 21600 : call gauss_points(ngpc,gsweights,gspts) !set gauss points/weights
191 86400 : gspts = 0.5_r8*(gspts+1.0_r8) !shift location so in [0:1] instead of [-1:1]
192 :
193 86400 : x(1:nc) = vtx_cart(1,1,1:nc,1 )
194 86400 : y(1:nc) = vtx_cart(1,2,1 ,1:nc)
195 21600 : x(nc+1) = vtx_cart(2,1, nc,1 )
196 21600 : y(nc+1) = vtx_cart(3,2,1 ,nc )
197 :
198 0 : select case (irecons)
199 : case(1)
200 0 : if (present(lanalytic)) then
201 0 : do j=1,nc
202 0 : do i=1,nc
203 0 : area(i,j) = (I_00(x(i+1),y(j+1)) - I_00(x(i),y(j+1)) + &
204 0 : I_00(x(i),y(j)) - I_00(x(i+1),y(j)))
205 : end do
206 : end do
207 : else
208 0 : call endrun("non-analytic moments not coded for irecons=1")
209 : end if
210 :
211 : case(3)
212 0 : if (present(lanalytic)) then
213 0 : do j=1,nc
214 0 : do i=1,nc
215 0 : area(i,j) = (I_00(x(i+1),y(j+1)) - I_00(x(i),y(j+1)) + &
216 0 : I_00(x(i),y(j)) - I_00(x(i+1),y(j)))
217 : ! Compute centroids via line integrals
218 0 : spherecentroid(1,i,j) = (I_10(x(i+1),y(j+1)) - I_10(x(i),y(j+1)) + &
219 0 : I_10(x(i),y(j)) - I_10(x(i+1),y(j))) / area(i,j)
220 0 : spherecentroid(2,i,j) = (I_01(x(i+1),y(j+1)) - I_01(x(i),y(j+1)) + &
221 0 : I_01(x(i),y(j)) - I_01(x(i+1),y(j))) / area(i,j)
222 : end do
223 : end do
224 : else
225 0 : call endrun("non-analytic moments not coded for irecons=3")
226 : end if
227 :
228 :
229 : case(6)
230 21600 : if (present(lanalytic)) then
231 86400 : do j=1,nc
232 280800 : do i=1,nc
233 : ! area(i,j) = surfareaxy(x(i),x(i+1),y(j),y(j+1))
234 777600 : area(i,j) = (I_00(x(i+1),y(j+1)) - I_00(x(i),y(j+1)) + &
235 972000 : I_00(x(i),y(j)) - I_00(x(i+1),y(j)))
236 : ! Compute centroids via line integrals
237 194400 : spherecentroid(1,i,j) = (I_10(x(i+1),y(j+1)) - I_10(x(i),y(j+1)) + &
238 194400 : I_10(x(i),y(j)) - I_10(x(i+1),y(j))) / area(i,j)
239 194400 : spherecentroid(2,i,j) = (I_01(x(i+1),y(j+1)) - I_01(x(i),y(j+1)) + &
240 194400 : I_01(x(i),y(j)) - I_01(x(i+1),y(j))) / area(i,j)
241 : ! TAN(alpha)^2 component
242 194400 : spherecentroid(3,i,j) = (I_20(x(i+1),y(j+1)) - I_20(x(i),y(j+1)) + &
243 194400 : I_20(x(i),y(j)) - I_20(x(i+1),y(j))) / area(i,j)
244 : ! TAN(beta)^2 component
245 194400 : spherecentroid(4,i,j) = (I_02(x(i+1),y(j+1)) - I_02(x(i),y(j+1)) + &
246 194400 : I_02(x(i),y(j)) - I_02(x(i+1),y(j))) / area(i,j)
247 : ! TAN(alpha) TAN(beta) component
248 194400 : spherecentroid(5,i,j) = (I_11(x(i+1),y(j+1)) - I_11(x(i),y(j+1)) + &
249 259200 : I_11(x(i),y(j)) - I_11(x(i+1),y(j))) / area(i,j)
250 : end do
251 : end do
252 : else
253 0 : do j=1,nc
254 0 : do i=1,nc
255 :
256 0 : xx (1,1,1) = x(i) ; xx (2,1,1) = y(j+1);
257 0 : dxx(1,1,1) = x(i+1)-x(i); dxx(2,1,1) = 0.0_r8 ;
258 :
259 0 : xx (1,2,1) = x(i+1) ; xx (2,2,1) = y(j) ;
260 0 : dxx(1,2,1) = x(i)-x(i+1); dxx(2,2,1) = 0.0_r8 ;
261 :
262 0 : call get_high_order_weights_over_areas(xx,dxx,num_seg,num_seg_max,num_area,weights,ngpc,gsweights,gspts,irecons)
263 :
264 0 : area(i,j) = weights(1,1)
265 :
266 0 : spherecentroid(1:5,i,j) = weights(2:6,1)/area(i,j)
267 : end do
268 : end do
269 : end if
270 : case default
271 21600 : call endrun('SUBROUTINE moment_on_sphere: irecons out of range')
272 : end select
273 21600 : end subroutine moment_onsphere
274 :
275 :
276 : ! ----------------------------------------------------------------------------------!
277 : !SUBROUTINES I_00, I_01, I_20, I_02, I11----------------------------------CE-for FVM!
278 : ! AUTHOR: CHRISTOPH ERATH, 17.October 2011 !
279 : ! DESCRIPTION: calculates the exact integrals !
280 : ! !
281 : ! CALLS: none !
282 : ! INPUT: x ... x coordinate of the evaluation point (Cartesian on the cube) !
283 : ! y ... y coordinate of the evaluation point (Cartesian on the cube) !
284 : ! OUTPUT: I_00, I_01, I_20, I_02, I11 !
285 : !-----------------------------------------------------------------------------------!
286 1555200 : function I_00(x,y)
287 : implicit none
288 : real (kind=r8) :: I_00
289 : real (kind=r8), intent(in) :: x,y
290 :
291 1555200 : I_00 = ATAN(x*y/SQRT(1.0_r8+x*x+y*y))
292 1555200 : end function I_00
293 :
294 1555200 : function I_10(x,y)
295 : implicit none
296 : real (kind=r8) :: I_10
297 : real (kind=r8), intent(in) :: x,y
298 : real (kind=r8) :: tmp
299 :
300 : ! tmp = ATAN(x)
301 : ! I_10 = -ASINH(y*COS(tmp))
302 1555200 : tmp = y*COS(ATAN(x))
303 1555200 : I_10 = -log(tmp+sqrt(tmp*tmp+1))
304 1555200 : end function I_10
305 :
306 :
307 1555200 : function I_01(x,y)
308 : implicit none
309 : real (kind=r8) :: I_01
310 : real (kind=r8), intent(in) :: x,y
311 : real (kind=r8) :: tmp
312 :
313 : ! I_01 = -ASINH(x/SQRT(1+y*y))
314 1555200 : tmp=x/SQRT(1+y*y)
315 1555200 : I_01 = -log(tmp+sqrt(tmp*tmp+1))
316 1555200 : end function I_01
317 :
318 1555200 : function I_20(x,y)
319 : implicit none
320 : real (kind=r8) :: I_20
321 : real (kind=r8), intent(in) :: x,y
322 : real (kind=r8) :: tmp,tmp1
323 :
324 1555200 : tmp = 1.0_r8+y*y
325 1555200 : tmp1=x/SQRT(tmp)
326 1555200 : I_20 = y*log(tmp1+sqrt(tmp1*tmp1+1))+ACOS(x*y/(SQRT((1.0_r8+x*x)*tmp)))
327 1555200 : end function I_20
328 :
329 1555200 : function I_02(x,y)
330 : implicit none
331 : real (kind=r8) :: I_02
332 : real (kind=r8), intent(in) :: x,y
333 : real (kind=r8) :: tmp,tmp1
334 :
335 : ! tmp=1.0_r8+x*x
336 : ! I_02 = x*ASINH(y/SQRT(tmp))+ACOS(x*y/SQRT(tmp*(1+y*y)))
337 1555200 : tmp=1.0_r8+x*x
338 1555200 : tmp1=y/SQRT(tmp)
339 :
340 1555200 : I_02 = x*log(tmp1+sqrt(tmp1*tmp1+1))+ACOS(x*y/SQRT(tmp*(1+y*y)))
341 :
342 1555200 : end function I_02
343 :
344 1555200 : function I_11(x,y)
345 : implicit none
346 : real (kind=r8) :: I_11
347 : real (kind=r8), intent(in) :: x,y
348 :
349 1555200 : I_11 = -SQRT(1+x*x+y*y)
350 1555200 : end function I_11
351 : !END SUBROUTINES I_00, I_01, I_20, I_02, I11------------------------------CE-for FVM!
352 :
353 :
354 0 : real (kind=r8) function F_00(x_in,y_in)
355 : implicit none
356 : real (kind=r8), intent(in) :: x_in,y_in
357 : real (kind=r8) :: x,y
358 : !
359 0 : x = x_in
360 0 : y = y_in
361 0 : F_00 =y/((1.0_r8+x*x)*SQRT(1.0_r8+x*x+y*y))
362 0 : end function F_00
363 :
364 0 : real (kind=r8) function F_10(x_in,y_in)
365 : implicit none
366 : real (kind=r8), intent(in) :: x_in,y_in
367 : real (kind=r8) :: x,y
368 :
369 0 : x = x_in
370 0 : y = y_in
371 :
372 0 : F_10 =x*y/((1.0_r8+x*x)*SQRT(1.0_r8+x*x+y*y))
373 0 : end function F_10
374 :
375 0 : real (kind=r8) function F_01(x_in,y_in)
376 : implicit none
377 : real (kind=r8), intent(in) :: x_in,y_in
378 : real (kind=r8) :: x,y
379 :
380 0 : x = x_in
381 0 : y = y_in
382 :
383 0 : F_01 =-1.0_r8/(SQRT(1.0_r8+x*x+y*y))
384 0 : end function F_01
385 :
386 0 : real (kind=r8) function F_20(x_in,y_in)
387 : implicit none
388 : real (kind=r8), intent(in) :: x_in,y_in
389 : real (kind=r8) :: x,y
390 :
391 0 : x = x_in
392 0 : y = y_in
393 :
394 0 : F_20 =x*x*y/((1.0_r8+x*x)*SQRT(1.0_r8+x*x+y*y))
395 0 : end function F_20
396 :
397 0 : real (kind=r8) function F_02(x_in,y_in)
398 : implicit none
399 : real (kind=r8), intent(in) :: x_in,y_in
400 : real (kind=r8) :: x,y,alpha,tmp
401 :
402 0 : x = x_in
403 0 : y = y_in
404 :
405 : alpha = ATAN(x)
406 : ! F_02 =-y/SQRT(1.0_r8+x*x+y*y)+ASINH(y*COS(alpha))
407 0 : tmp=y*COS(alpha)
408 0 : F_02 =-y/SQRT(1.0_r8+x*x+y*y)+log(tmp+sqrt(tmp*tmp+1))
409 :
410 : !
411 : ! cos(alpha) = 1/sqrt(1+x*x)
412 : !
413 0 : end function F_02
414 :
415 0 : real (kind=r8) function F_11(x_in,y_in)
416 : implicit none
417 : real (kind=r8), intent(in) :: x_in,y_in
418 : real (kind=r8) :: x,y
419 :
420 0 : x = x_in
421 0 : y = y_in
422 :
423 0 : F_11 =-x/(SQRT(1.0_r8+x*x+y*y))
424 0 : end function F_11
425 :
426 :
427 :
428 : !
429 : ! matrix version of reconstruct_cubic_onface
430 : !
431 21600 : subroutine compute_reconstruct_matrix(nc,nhe,nhc,irecons,dalpha,dbeta,spherecentroid,vtx_cart,&
432 21600 : centroid_stretch,vertex_recons_weights,recons_metrics,recons_metrics_integral)
433 : implicit none
434 : integer , intent(in) :: nc,nhe,irecons,nhc
435 : real (kind=r8), intent(in) :: dalpha,dbeta
436 : real (kind=r8), dimension(irecons-1,1-nhc:nc+nhc,1-nhc:nc+nhc), intent(in) :: spherecentroid
437 : real (kind=r8), dimension(4,2,1-nhc:nc+nhc,1-nhc:nc+nhc) , intent(in) :: vtx_cart
438 :
439 : real (kind=r8), dimension(7,1-nhe:nc+nhe,1-nhe:nc+nhe) , intent(out):: centroid_stretch
440 : real (kind=r8), dimension(4,1:irecons-1,1-nhe:nc+nhe,1-nhe:nc+nhe), intent(out):: vertex_recons_weights
441 : real (kind=r8), dimension(3,1-nhe:nc+nhe,1-nhe:nc+nhe) , intent(out):: recons_metrics
442 : real (kind=r8), dimension(3,1-nhe:nc+nhe,1-nhe:nc+nhe) , intent(out):: recons_metrics_integral
443 :
444 : !
445 : integer :: i, j, count, m, n
446 : real (kind=r8) :: coef,tmp,cartx,carty
447 : !
448 : ! pre-compute variables for reconstruction
449 : !
450 : select case (irecons)
451 : case(3)
452 0 : do j= 1-nhe,nc+nhe
453 0 : do i=1-nhe,nc+nhe
454 0 : count = 1
455 0 : do n = j, j+1
456 0 : do m = i, i+1
457 0 : cartx = vtx_cart(count,1,i,j); carty = vtx_cart(count,2,i,j);
458 :
459 0 : vertex_recons_weights(count,1,i,j) = cartx - spherecentroid(1,i,j)
460 0 : vertex_recons_weights(count,2,i,j) = carty - spherecentroid(2,i,j)
461 :
462 0 : count=count+1
463 : end do
464 : enddo
465 : end do
466 : end do
467 0 : call endrun("recons_metrics and recons_metrics_integral not initialize")
468 : !
469 : ! for reconstruction
470 : !
471 0 : do j= 1-nhe,nc+nhe
472 0 : do i=1-nhe,nc+nhe
473 : !
474 : !***************
475 : !* dfdx *
476 : !***************
477 : !
478 0 : coef = 1.0_r8/(12.0_r8 * dalpha) !finite difference coefficient
479 0 : coef = coef /( 1.0_r8 + spherecentroid(1,i,j)**2) !stretching coefficient
480 :
481 0 : centroid_stretch(1,i,j) = coef
482 : !
483 : !***************
484 : !* dfdy *
485 : !***************
486 : !
487 0 : coef = 1.0_r8/(12.0_r8 * dbeta) !finite difference coefficient
488 0 : coef = coef /( 1.0_r8 + spherecentroid(2,i,j)**2) !stretching coefficient
489 :
490 0 : centroid_stretch(2,i,j) = coef
491 : end do
492 : end do
493 : case(6)
494 108000 : do j= 1-nhe,nc+nhe
495 475200 : do i=1-nhe,nc+nhe
496 1922400 : do count=1,4
497 1468800 : cartx = vtx_cart(count,1,i,j); carty = vtx_cart(count,2,i,j);
498 :
499 1468800 : vertex_recons_weights(count,1,i,j) = cartx - spherecentroid(1,i,j)
500 1468800 : vertex_recons_weights(count,2,i,j) = carty - spherecentroid(2,i,j)
501 :
502 1468800 : vertex_recons_weights(count,3,i,j) = (spherecentroid(1,i,j)**2 - &
503 : spherecentroid(3,i,j)) + &
504 1468800 : (cartx - spherecentroid(1,i,j))**2
505 1468800 : vertex_recons_weights(count,4,i,j) = (spherecentroid(2,i,j)**2 - &
506 : spherecentroid(4,i,j)) + &
507 1468800 : (carty - spherecentroid(2,i,j))**2
508 :
509 1468800 : vertex_recons_weights(count,5,i,j) = (cartx - spherecentroid(1,i,j))* &
510 : (carty - spherecentroid(2,i,j))+ &
511 : (spherecentroid(1,i,j) * &
512 : spherecentroid(2,i,j) - &
513 1836000 : spherecentroid(5,i,j))
514 : end do
515 : end do
516 : end do
517 :
518 108000 : do j= 1-nhe,nc+nhe
519 475200 : do i=1-nhe,nc+nhe
520 367200 : recons_metrics(1,i,j) = spherecentroid(1,i,j)**2 -spherecentroid(3,i,j)
521 367200 : recons_metrics(2,i,j) = spherecentroid(2,i,j)**2 -spherecentroid(4,i,j)
522 : recons_metrics(3,i,j) = spherecentroid(1,i,j)*spherecentroid(2,i,j)-&
523 367200 : spherecentroid(5,i,j)
524 :
525 : recons_metrics_integral(1,i,j) = &
526 367200 : 2.0_r8*spherecentroid(1,i,j)**2 -spherecentroid(3,i,j)
527 : recons_metrics_integral(2,i,j) = &
528 367200 : 2.0_r8*spherecentroid(2,i,j)**2 -spherecentroid(4,i,j)
529 : recons_metrics_integral(3,i,j) = &
530 : 2.0_r8*spherecentroid(1,i,j)*spherecentroid(2,i,j)-&
531 453600 : spherecentroid(5,i,j)
532 : end do
533 : end do
534 :
535 :
536 :
537 : !
538 : ! pre-compute variables for reconstruction
539 : !
540 108000 : do j= 1-nhe,nc+nhe
541 475200 : do i=1-nhe,nc+nhe
542 : !
543 : !***************
544 : !* dfdx *
545 : !***************
546 : !
547 367200 : coef = 1.0_r8/(12.0_r8 * dalpha) !finite difference coefficient
548 367200 : coef = coef /( 1.0_r8 + spherecentroid(1,i,j)**2) !stretching coefficient
549 :
550 367200 : centroid_stretch(1,i,j) = coef
551 : !
552 : !***************
553 : !* dfdy *
554 : !***************
555 : !
556 367200 : coef = 1.0_r8/(12.0_r8 * dbeta) !finite difference coefficient
557 367200 : coef = coef /( 1.0_r8 + spherecentroid(2,i,j)**2) !stretching coefficient
558 :
559 367200 : centroid_stretch(2,i,j) = coef
560 :
561 : !*****************
562 : !* d2fdx2 *
563 : !*****************
564 : !
565 367200 : coef = 1.0_r8 / (12.0_r8 * dalpha**2) !finite difference coefficient
566 : !
567 : ! stretching coefficient part 2
568 : ! recons(3,i,j) = (a * recons(1,i,j)+ recons(3,i,j))*b
569 : !
570 367200 : tmp = 0.5_r8/((1.0_r8 + spherecentroid(1,i,j)**2)**2)
571 :
572 367200 : centroid_stretch(3,i,j) = coef*tmp
573 367200 : centroid_stretch(6,i,j) = -spherecentroid(1,i,j)/(1.0_r8 + spherecentroid(1,i,j)**2)
574 :
575 : !
576 : !*****************
577 : !* d2fdy2 *
578 : !*****************
579 : !
580 : !
581 367200 : coef = 1.0_r8 / (12.0_r8 * dbeta**2) !finite difference coefficient
582 : !
583 : ! stretching coefficient part 2
584 : !
585 : ! recons(4,i,j) = (a * recons(1,i,j)+ recons(4,i,j))*b
586 : !
587 367200 : tmp =0.5_r8/((1.0_r8 + spherecentroid(2,i,j)**2)**2)
588 :
589 367200 : centroid_stretch(4,i,j) = coef*tmp
590 367200 : centroid_stretch(7,i,j) = -spherecentroid(2,i,j)/(1.0_r8 + spherecentroid(2,i,j)**2)
591 : !
592 : !*****************
593 : !* d2fdxdy *
594 : !*****************
595 : !
596 : !
597 367200 : coef = 1.0_r8 / (4.0_r8 * dalpha * dbeta) !finite difference coefficient
598 : coef = coef / ((1.0_r8 + spherecentroid(1,i,j)**2) * &
599 367200 : (1.0_r8 + spherecentroid(2,i,j)**2)) !stretching coefficient
600 :
601 453600 : centroid_stretch(5,i,j) = coef
602 : enddo
603 : enddo
604 : case default
605 21600 : call endrun('SUBROUTINE compute_reconstruct_matrix: irecons out of range')
606 : end select
607 21600 : end subroutine compute_reconstruct_matrix
608 :
609 :
610 11594793592 : subroutine get_high_order_weights_over_areas(x,dx,num_seg,num_seg_max,num_area,weights,ngpc,gsweights, gspts,irecons)
611 : implicit none
612 : integer , intent(in) :: num_area, num_seg_max, irecons
613 : REAL(KIND=r8), dimension(2,num_seg_max,num_area ), intent(inout) :: x, dx
614 : integer , intent(in) :: ngpc
615 : integer , dimension(num_area ), intent(in) :: num_seg
616 : REAL(KIND=r8), dimension(irecons,num_area), intent(out) :: weights
617 :
618 23189587184 : real (kind=r8), dimension(ngpc,num_seg_max ) :: xq,yq !quadrature points along line segments
619 23189587184 : real (kind=r8), dimension(ngpc,num_seg_max,irecons) :: F !potentials
620 23189587184 : real (kind=r8), dimension( irecons) :: weights_area
621 23189587184 : real (kind=r8), dimension(ngpc,num_seg_max) :: xq2, yrh, rho, tmp !intermediate variables for optimization
622 11594793592 : REAL(KIND=r8) , dimension(ngpc,num_seg_max) :: xq2ir, xq2i, rhoi !intermediate variables for optimization
623 :
624 : integer :: iseg,iarea,i,j,k
625 :
626 : real (kind=r8), dimension(ngpc) :: gsweights, gspts
627 :
628 >41741*10^7 : weights(1:irecons,1:num_area) = 0.0_r8 !may not be necessary dbgxxx
629 69568761552 : do iarea=1,num_area
630 >14061*10^7 : do iseg=1,num_seg(iarea)
631 >33058*10^7 : xq(:,iseg) = x(1,iseg,iarea)+dx(1,iseg,iarea)*gspts(:)
632 >38855*10^7 : yq(:,iseg) = x(2,iseg,iarea)+dx(2,iseg,iarea)*gspts(:)
633 : end do
634 : !
635 : ! potentials (equation's 23-28 in CSLAM paper; Lauritzen et al., 2010):
636 : !
637 : ! (Rory Kelly optimization)
638 : !
639 >14061*10^7 : do j=1,num_seg(iarea)
640 : !DIR$ SIMD
641 >33058*10^7 : do i=1,ngpc
642 >24793*10^7 : xq2(i,j) = xq(i,j)*xq(i,j)
643 >24793*10^7 : xq2i(i,j) = 1.0_r8/(1.0_r8+xq2(i,j))
644 >24793*10^7 : xq2ir(i,j) = SQRT(xq2i(i,j))
645 >24793*10^7 : rho(i,j) = SQRT(1.0_r8+xq2(i,j)+yq(i,j)*yq(i,j))
646 >24793*10^7 : rhoi(i,j) = 1.0_r8/rho(i,j)
647 >24793*10^7 : yrh(i,j) = yq(i,j)*rhoi(i,j)
648 >24793*10^7 : tmp(i,j) = yq(i,j)*xq2ir(i,j)
649 >24793*10^7 : F(i,j,1) = yrh(i,j)*xq2i(i,j) !F_00 !F_00
650 >24793*10^7 : F(i,j,2) = xq(i,j)*yrh(i,j)*xq2i(i,j) !F_10 !F_10
651 >24793*10^7 : F(i,j,3) = -1.0_r8*rhoi(i,j) !F_01 !F_01
652 >24793*10^7 : F(i,j,4) = xq2(i,j)*yrh(i,j)*xq2i(i,j) !F_20 !F_20
653 >33058*10^7 : F(i,j,6) = -xq(i,j)*rhoi(i,j) !F_11 !F_11
654 : enddo
655 : !
656 : ! take F(i,j,5) out of loop above since it prevents vectorization
657 : !
658 >38855*10^7 : do i=1,ngpc
659 >33058*10^7 : F(i,j,5) = -yq(i,j)*rhoi(i,j)+log(tmp(i,j)+rho(i,j)*xq2ir(i,j)) !F_02 !F_02
660 : end do
661 : enddo
662 >40581*10^7 : weights_area = 0.0_r8
663 >40581*10^7 : do k=1,irecons
664 >90169*10^7 : do iseg=1,num_seg(iarea)
665 >23313*10^8 : weights_area(k) = weights_area(k) + sum(gsweights(:)*F(:,iseg,k))*0.5_r8*dx(1,iseg,iarea)
666 : end do
667 : end do
668 >41741*10^7 : weights(1:irecons,iarea) = weights_area(1:irecons)
669 : end do
670 11594793592 : end subroutine get_high_order_weights_over_areas
671 :
672 :
673 : !********************************************************************************
674 : !
675 : ! Gauss-Legendre quadrature
676 : !
677 : ! Tabulated values
678 : !
679 : !********************************************************************************
680 760416 : subroutine gauss_points(n,weights,points)
681 : implicit none
682 : integer, intent(in ) :: n
683 : real (kind=r8), dimension(:), intent(out) :: weights, points !dimension(n)
684 :
685 760416 : select case (n)
686 : ! CASE(1)
687 : ! abscissae(1) = 0.0_r8
688 : ! weights(1) = 2.0_r8
689 : case(2)
690 0 : points(1) = -sqrt(1.0_r8/3.0_r8)
691 0 : points(2) = sqrt(1.0_r8/3.0_r8)
692 0 : weights(1) = 1.0_r8
693 0 : weights(2) = 1.0_r8
694 : case(3)
695 760416 : points(1) = -0.774596669241483377035853079956_r8
696 760416 : points(2) = 0.0_r8
697 760416 : points(3) = 0.774596669241483377035853079956_r8
698 760416 : weights(1) = 0.555555555555555555555555555556_r8
699 760416 : weights(2) = 0.888888888888888888888888888889_r8
700 760416 : weights(3) = 0.555555555555555555555555555556_r8
701 : case(4)
702 0 : points(1) = -0.861136311594052575223946488893_r8
703 0 : points(2) = -0.339981043584856264802665659103_r8
704 0 : points(3) = 0.339981043584856264802665659103_r8
705 0 : points(4) = 0.861136311594052575223946488893_r8
706 0 : weights(1) = 0.347854845137453857373063949222_r8
707 0 : weights(2) = 0.652145154862546142626936050778_r8
708 0 : weights(3) = 0.652145154862546142626936050778_r8
709 0 : weights(4) = 0.347854845137453857373063949222_r8
710 : case(5)
711 0 : points(1) = -(1.0_r8/3.0_r8)*sqrt(5.0_r8+2.0_r8*sqrt(10.0_r8/7.0_r8))
712 0 : points(2) = -(1.0_r8/3.0_r8)*sqrt(5.0_r8-2.0_r8*sqrt(10.0_r8/7.0_r8))
713 0 : points(3) = 0.0_r8
714 0 : points(4) = (1.0_r8/3.0_r8)*sqrt(5.0_r8-2.0_r8*sqrt(10.0_r8/7.0_r8))
715 0 : points(5) = (1.0_r8/3.0_r8)*sqrt(5.0_r8+2.0_r8*sqrt(10.0_r8/7.0_r8))
716 0 : weights(1) = (322.0_r8-13.0_r8*sqrt(70.0_r8))/900.0_r8
717 0 : weights(2) = (322.0_r8+13.0_r8*sqrt(70.0_r8))/900.0_r8
718 0 : weights(3) = 128.0_r8/225.0_r8
719 0 : weights(4) = (322.0_r8+13.0_r8*sqrt(70.0_r8))/900.0_r8
720 0 : weights(5) = (322.0_r8-13.0_r8*sqrt(70.0_r8))/900.0_r8
721 : case default
722 760416 : call endrun('SUBROUTINE gauss_points: n out of range in (0<n<5)')
723 : end select
724 :
725 760416 : end subroutine gauss_points
726 :
727 :
728 :
729 21600 : subroutine init_flux_orient(flux_orient,ifct,nc,nhc,cubeboundary,faceno)
730 : implicit none
731 : integer , intent(in) :: cubeboundary,faceno,nc,nhc
732 : real (kind=r8), intent(out) :: flux_orient(2 ,1-nhc:nc+nhc,1-nhc:nc+nhc)
733 : integer , intent(out) :: ifct (1-nhc:nc+nhc,1-nhc:nc+nhc)
734 :
735 : integer :: ib
736 :
737 5464800 : flux_orient = 99.9D9 !for debugging
738 : !
739 : ! halo of flux_orient will be filled through boundary exchange
740 : !
741 280800 : flux_orient (1,1:nc,1:nc) = dble(faceno)
742 1965600 : flux_orient (2,:,:) = 0.0_r8
743 1965600 : ifct(:,:) = 1
744 21600 : if (cubeboundary>0) then
745 :
746 : !
747 : ! cshift (permute) value needed to be applied to vertex number so that they match orientation
748 : ! of the interior of the panel
749 : !
750 : !
751 2784 : ib = cubeboundary
752 2784 : if (faceno==2) then
753 4184 : if (ib==north.or.ib==nwest.or.ib==neast) flux_orient(2,1-nhc:nc+nhc,nc+1 :nc+nhc) = 1
754 4184 : if (ib==south.or.ib==swest.or.ib==seast) flux_orient(2,1-nhc:nc+nhc,1-nhc:0 ) = 3
755 : end if
756 2784 : if (faceno==3) then
757 4184 : if (ib==north.or.ib==nwest.or.ib==neast) flux_orient (2,1-nhc:nc+nhc,nc+1 :nc+nhc) = 2
758 4184 : if (ib==south.or.ib==swest.or.ib==seast) flux_orient (2,1-nhc:nc+nhc,1-nhc:0 ) = 2
759 : end if
760 2784 : if (faceno==4) then
761 4184 : if (ib==north.or.ib==nwest.or.ib==neast) flux_orient (2,1-nhc:nc+nhc,nc+1 :nc+nhc) = 3
762 4184 : if (ib==south.or.ib==swest.or.ib==seast) flux_orient (2,1-nhc:nc+nhc,1-nhc:0 ) = 1
763 : end if
764 2784 : if (faceno==5) then
765 4184 : if (ib==south.or.ib==swest.or.ib==seast) flux_orient (2,1-nhc:nc+nhc,1-nhc:0 ) = 2
766 4904 : if (ib== west.or.ib==swest.or.ib==nwest) flux_orient (2,1-nhc:0 ,1-nhc:nc+nhc) = 3
767 4904 : if (ib== east.or.ib==seast.or.ib==neast) flux_orient (2, nc+1:nc+nhc,1-nhc:nc+nhc) = 1
768 : end if
769 :
770 2784 : if (faceno==6) then
771 4184 : if (ib==north.or.ib==nwest.or.ib==neast ) flux_orient (2,1-nhc:nc+nhc,nc+1 :nc+nhc) = 2
772 4904 : if (ib==west .or.ib==swest.or.ib==nwest ) flux_orient (2,1-nhc:0 ,1-nhc:nc+nhc) = 1
773 4904 : if (ib==east .or.ib==seast.or.ib==neast ) flux_orient (2,nc+1:nc+nhc,1-nhc:nc+nhc) = 3
774 : end if
775 : !
776 : ! non-existent cells in physical space
777 : !
778 2784 : if (cubeboundary==nwest) then
779 336 : flux_orient(2,1-nhc:0 ,nc+1 :nc+nhc) = 0
780 312 : ifct ( 1-nhc:0 ,nc+1 :nc+nhc) = 0
781 2760 : else if (cubeboundary==swest) then
782 336 : flux_orient (2,1-nhc:0 ,1-nhc:0 ) = 0
783 312 : ifct ( 1-nhc:0 ,1-nhc:0 ) = 0
784 2736 : else if (cubeboundary==neast) then
785 336 : flux_orient (2,nc+1 :nc+nhc,nc+1 :nc+nhc) = 0
786 312 : ifct ( nc+1 :nc+nhc,nc+1 :nc+nhc) = 0
787 2712 : else if (cubeboundary==seast) then
788 336 : flux_orient (2,nc+1 :nc+nhc,1-nhc:0 ) = 0
789 312 : ifct ( nc+1 :nc+nhc,1-nhc:0 ) = 0
790 : end if
791 : end if
792 :
793 21600 : end subroutine init_flux_orient
794 :
795 : !
796 : !
797 : !
798 :
799 : ! ----------------------------------------------------------------------------------!
800 : !SUBROUTINE CREATE_INTERPOLATIION_POINTS----------------------------------CE-for FVM!
801 : ! AUTHOR: CHRISTOPH ERATH, 17.October 2011 !
802 : ! DESCRIPTION: for elements, which share a cube edge, we have to do some !
803 : ! interpolation on different cubic faces, also in the halo region: !
804 : ! because we also need the reconstruction coefficients in the halo zone, !
805 : ! which is basically calculated twice, on the original cell of an element !
806 : ! on face A and on a cell in the halo region of an element of face B !
807 : ! The crux is, that the interpolation has to be the same to ensure !
808 : ! conservation of the scheme !
809 : ! SYMMETRY of the CUBE is used for calucaltion the interpolation_point !
810 : ! !
811 : ! CALLS: interpolation_point !
812 : ! INPUT/OUTPUT: !
813 : ! elem ... element structure from HOMME !
814 : ! fvm ... structure !
815 : !-----------------------------------------------------------------------------------!
816 21600 : subroutine create_interpolation_points(elem,&
817 : nc,nhc,nhr,ns,nh,cubeboundary,&
818 21600 : dalpha,dbeta,ibase,halo_interp_weight)
819 : use element_mod , only: element_t
820 : use coordinate_systems_mod, only: cartesian2D_t
821 : use control_mod , only: north, south, east, west, neast, nwest, seast, swest
822 : use cube_mod , only: cube_xstart, cube_xend, cube_ystart, cube_yend
823 :
824 : implicit none
825 : type (element_t), intent(in) :: elem
826 :
827 : integer , intent(in) :: nc,nhc,nhr,ns,nh,cubeboundary
828 : integer , intent(out) :: ibase(1-nh:nc+nh,1:nhr,2)
829 : real (kind=r8), intent(out) :: halo_interp_weight(1:ns,1-nh:nc+nh,1:nhr,2)
830 : !
831 : ! pre-compute weight/index matrices
832 : !
833 : integer :: imin,imax,jmin,jmax,iinterp
834 : real (kind=r8), intent(in) :: dalpha,dbeta
835 :
836 43200 : real (kind=r8), dimension(1-nhc:nc+nhc) :: gnomxstart, gnomxend, gnomystart, gnomyend
837 : integer :: i, halo, ida, ide, iref1
838 : type (cartesian2D_t) :: tmpgnom
839 43200 : real (kind=r8) :: interp(1-nh:nc+nh,1:nhr,2)
840 : integer ::ibaseref
841 21600 : integer :: ibase_tmp(1-nh:nc+nh,1:nhr,2)
842 :
843 756000 : ibase = 99999 !dbg
844 2570400 : halo_interp_weight(:,:,:,:) = 9.99E9_r8 !dbg
845 :
846 : ! element is not on a corner, but shares a cube edge (call of subroutine)
847 21600 : if(cubeboundary <= 4) then
848 21504 : gnomxstart(1-nhc)=elem%corners(1)%x-(nhc-0.5_r8)*dalpha
849 21504 : gnomystart(1-nhc)=elem%corners(1)%y-(nhc-0.5_r8)*dbeta
850 193536 : do i=2-nhc,nc+nhc
851 172032 : gnomxstart(i)=gnomxstart(i-1)+dalpha
852 193536 : gnomystart(i)=gnomystart(i-1)+dbeta
853 : end do
854 21504 : ida=1-nhc !lower bound
855 21504 : ide=nc+nhc !upper bound
856 : select case (cubeboundary)
857 : !INTERIOR element
858 : case(0)
859 : ! nothing to do!
860 : !CASE WEST
861 : case(west)
862 2016 : do halo=1,nhr
863 : ! iref1=ida
864 1344 : tmpgnom%x=cube_xstart-(halo-0.5_r8)*dalpha
865 10080 : do i=halo-nh,nc+nh-(halo-1) !see fvm_reconstruction to understand these boundaries
866 8064 : iref1=ida
867 8064 : tmpgnom%y=gnomystart(i)
868 8064 : call interpolation_point(nc,ns,tmpgnom,gnomystart,1,4,1,interp(i,halo,1),&
869 17472 : ida,ide,iref1,ibase_tmp(i,halo,1))
870 : end do
871 : end do
872 :
873 : !CASE EAST
874 : case(east)
875 : ! east zone
876 2016 : do halo=1,nhr
877 1344 : iref1=ida
878 1344 : tmpgnom%x=cube_xend+(halo-0.5_r8)*dalpha
879 10080 : do i=halo-nh,nc+nh-(halo-1)
880 8064 : tmpgnom%y=gnomystart(i)
881 8064 : call interpolation_point(nc,ns,tmpgnom,gnomystart,1,2,1,interp(i,halo,1),&
882 17472 : ida,ide,iref1,ibase_tmp(i,halo,1))
883 : end do
884 : end do
885 :
886 : !CASE NORTH
887 : case(north)
888 : ! north zone
889 2016 : do halo=1,nhr
890 1344 : tmpgnom%y=cube_yend+(halo-0.5_r8)*dbeta
891 1344 : iref1=ida
892 10080 : do i=halo-nh,nc+nh-(halo-1)
893 8064 : tmpgnom%x=gnomxstart(i)
894 : !
895 : ! dbg - change to interp(i,halo,1) instead of interp(i,halo,2)
896 : ! so that I can get rid of iinterp = 1 in fvm_reconstruction_mod
897 : !
898 8064 : call interpolation_point(nc,ns,tmpgnom,gnomxstart,1,6,0,interp(i,halo,2),&
899 17472 : ida,ide,iref1,ibase_tmp(i,halo,2))
900 : end do
901 : end do
902 : !CASE SOUTH
903 : case(south)
904 : !south zone
905 2016 : do halo=1,nhr
906 1344 : iref1=ida
907 1344 : tmpgnom%y=cube_ystart-(halo-0.5_r8)*dbeta
908 10080 : do i=halo-nh,nc+nh-(halo-1)
909 8064 : tmpgnom%x=gnomxstart(i)
910 8064 : call interpolation_point(nc,ns,tmpgnom,gnomxstart,1,5,0,interp(i,halo,2),&
911 17472 : ida,ide,iref1,ibase_tmp(i,halo,2))
912 : end do
913 : end do
914 :
915 : !
916 : !THIS CASE SHOULD NOT HAPPEN!
917 : case default
918 0 : print *,'Fatal Error in first select statement:'
919 21504 : call endrun('fvm_reconstruction_mod.F90 subroutine fillhalo_cubic!' )
920 : end select
921 : !CORNER TREATMENT
922 : else
923 96 : gnomxstart(1-nhc)=cube_xstart-(nhc-0.5_r8)*dalpha
924 96 : gnomxend(nc+nhc)=cube_xend+(nhc-0.5_r8)*dalpha
925 96 : gnomystart(1-nhc)=cube_ystart-(nhc-0.5_r8)*dbeta
926 96 : gnomyend(nc+nhc)=cube_yend+(nhc-0.5_r8)*dbeta
927 864 : do i=2-nhc,nc+nhc
928 768 : gnomxstart(i)=gnomxstart(i-1)+dalpha
929 768 : gnomxend(nc+1-i)=gnomxend(nc+2-i)-dalpha
930 768 : gnomystart(i)=gnomystart(i-1)+dbeta
931 864 : gnomyend(nc+1-i)=gnomyend(nc+2-i)-dbeta
932 : end do
933 :
934 : select case (cubeboundary)
935 : !CASE SOUTH WEST
936 : case(swest)
937 : ! west zone
938 72 : do halo=1,nhr
939 48 : tmpgnom%x=cube_xstart-(halo-0.5_r8)*dalpha
940 48 : ida=1
941 48 : ide=nc+nc
942 48 : iref1=ida
943 336 : do i=0,nc+nh-(halo-1)
944 264 : tmpgnom%y=gnomystart(i)
945 264 : call interpolation_point(nc,ns,tmpgnom,gnomystart,1,4,1,interp(i,halo,1),&
946 576 : ida,ide,iref1,ibase_tmp(i,halo,1))
947 : end do
948 : end do
949 : !CASE SOUTH EAST
950 : case(seast)
951 : ! east zone
952 72 : do halo=1,nhr
953 48 : tmpgnom%x=cube_xend+(halo-0.5_r8)*dalpha
954 48 : ida=1
955 48 : ide=nc+nc
956 48 : iref1=ida
957 336 : do i=0,nc+nh-(halo-1)
958 264 : tmpgnom%y=gnomystart(i)
959 264 : call interpolation_point(nc,ns,tmpgnom,gnomystart,1,2,1, interp(i,halo,1),&
960 576 : ida,ide,iref1,ibase_tmp(i,halo,1))
961 : end do
962 : end do
963 : !CASE NORTH EAST
964 : case(neast)
965 : ! east zone
966 72 : do halo=1,nhr
967 48 : tmpgnom%x=cube_xend+(halo-0.5_r8)*dalpha
968 48 : ida=1-nc
969 48 : ide=nc
970 48 : iref1=ida
971 336 : do i=halo-nh,nc+1
972 264 : tmpgnom%y=gnomyend(i)
973 264 : call interpolation_point(nc,ns,tmpgnom,gnomyend,1,2,1, interp(i,halo,1),&
974 576 : ida,ide,iref1,ibase_tmp(i,halo,1))
975 : end do
976 : end do
977 : !CASE NORTH WEST
978 : case(nwest)
979 : ! west zone
980 72 : do halo=1,2
981 48 : tmpgnom%x=cube_xstart-(halo-0.5_r8)*dalpha
982 48 : ida=1-nc
983 48 : ide=nc
984 48 : iref1=ida
985 336 : do i=halo-nh,nc+1
986 264 : tmpgnom%y=gnomyend(i)
987 264 : call interpolation_point(nc,ns,tmpgnom,gnomyend,1,4,1, interp(i,halo,1),&
988 576 : ida,ide,iref1,ibase_tmp(i,halo,1))
989 : end do
990 : end do
991 : !THIS CASE SHOULD NOT HAPPEN!
992 : case default
993 0 : print *,'Fatal Error in second select statement:'
994 96 : call endrun('fvm_reconstruction_mod.F90 subroutine create_interpolationpoint!')
995 : end select
996 : endif
997 :
998 : !**************************
999 : !
1000 : ! compute haloe weights and indices
1001 : !
1002 21600 : if (cubeboundary>0) then
1003 2784 : if (cubeboundary<5) then
1004 : !
1005 : ! element is located at a panel side but is not a corner element
1006 : ! (west,east,south,north) = (1,2,3,4)
1007 : !
1008 2688 : if (cubeboundary==west .or.cubeboundary==east ) then
1009 : iinterp = 1
1010 : end if
1011 2688 : if (cubeboundary==north.or.cubeboundary==south) iinterp = 2
1012 8064 : do halo=1,nhr
1013 40320 : do i=halo-nh,nc+nh-(halo-1)
1014 32256 : ibaseref=ibase_tmp(i,halo,iinterp)
1015 32256 : ibase(i,halo,1) = ibaseref
1016 : call get_equispace_weights(dbeta, interp(i,halo,iinterp),&
1017 37632 : halo_interp_weight(:,i,halo,1),ns)
1018 : end do
1019 : end do
1020 : else
1021 : !
1022 : ! element is located at a cube corner
1023 : ! (swest,seast,nwest,neast)=(5,6,7,8)
1024 : !
1025 288 : do halo=1,nhr
1026 192 : if (cubeboundary==swest .or.cubeboundary==seast) then
1027 96 : imin = 0 ; imax = nc+nh-(halo-1);
1028 96 : jmin = halo-nh; jmax = nc+1;
1029 : else
1030 96 : jmin = 0 ; jmax = nc+nh-(halo-1);
1031 96 : imin = halo-nh; imax = nc+1;
1032 : end if
1033 1248 : do i=imin,imax
1034 1056 : ibaseref=ibase_tmp(i,halo,1)
1035 1056 : ibase(i,halo,1) = ibaseref
1036 1248 : call get_equispace_weights(dbeta, interp(i,halo,1),halo_interp_weight(:,i,halo,1),ns)
1037 : end do
1038 : !
1039 : ! reverse weights/indices for fotherpanel (see details on reconstruct_matrix)
1040 : !
1041 4416 : halo_interp_weight(1:ns,jmin:jmax,halo,2) = halo_interp_weight(ns:1:-1,imax:imin:-1,halo,1)
1042 1344 : ibase (jmin:jmax,halo ,2) = nc+1-(ns-1)-ibase(imax:imin:-1,halo ,1)
1043 : end do
1044 : end if
1045 :
1046 : end if
1047 :
1048 :
1049 21600 : end subroutine create_interpolation_points
1050 :
1051 :
1052 :
1053 :
1054 : !END SUBROUTINE CREATE_INTERPOLATION_POINTS-------------------------------CE-for FVM!
1055 :
1056 :
1057 :
1058 : ! ----------------------------------------------------------------------------------!
1059 : !SUBROUTINE INTERPOLATION_POINT-------------------------------------------CE-for FVM!
1060 : ! AUTHOR: CHRISTOPH ERATH, 14.November 2011 !
1061 : ! DESCRIPTION: calculates the interpolation point on from face 1 in face 2 in !
1062 : ! alpha/beta coordinates, only 1D !
1063 : ! !
1064 : ! CALLS: cubedsphere2cart, cart2cubedsphere !
1065 : ! INPUT: gnom... 1D coordinates !
1066 : ! gnom1d... 1d coordinates !
1067 : ! face1... orginal face !
1068 : ! face2... target face (where the interpolation has to be done) !
1069 : ! xy ... 0 for alpha coordinate, any other for beta !
1070 : ! except.which type, interior, left edge (-1), right edge (1) !
1071 : ! point... interpolation point !
1072 : ! ida ... begin of interpval !
1073 : ! ide ... end of interpval !
1074 :
1075 :
1076 : ! INPUT/OUTPUT/RETURN: !
1077 : ! iref ... where we start the search, is also an OUTPUT, so we know for the !
1078 : ! next point where to start !
1079 : !-----------------------------------------------------------------------------------!
1080 : !
1081 : ! DESCRIPTION: searchs where the interpolation point has to be (iref), two values !
1082 : ! of interpval on the left and on the right, except if we are out of range !
1083 : ! which is indicated through ia and ie, respectively !
1084 : ! It is a 1D interpolation, use alpha/beta coordinates!!! !
1085 : ! !
1086 : ! CALLS: cubic_equispace_interp !
1087 : ! INPUT: iref ... where we start the search, is also an OUTPUT, so we know for the !
1088 : ! next point where to start !
1089 : ! ibaseref ... startindex of the four tracer value for the reconstruction !
1090 : ! point ... provides the difference of the interpolation point to use it !
1091 : ! directly in CUBIC_EQUISPACE_INTERP !
1092 : !-----------------------------------------------------------------------------------!
1093 33312 : function get_gno_point(gnom,face1,face2,xy) result(point)
1094 : use coordinate_systems_mod, only : cubedsphere2cart, cart2cubedsphere, &
1095 : cartesian2D_t,cartesian3D_t
1096 : implicit none
1097 : type (cartesian2D_t), intent(in) :: gnom
1098 : integer, intent(in) :: face1, face2, xy
1099 : real (kind=r8) :: point
1100 :
1101 : type(cartesian3D_t) :: tmpcart3d
1102 : type (cartesian2D_t) :: tmpgnom
1103 :
1104 33312 : tmpcart3d=cubedsphere2cart(gnom,face1)
1105 33312 : tmpgnom=cart2cubedsphere(tmpcart3d,face2)
1106 33312 : if(xy==0) then
1107 16128 : point=tmpgnom%x
1108 : else
1109 17184 : point=tmpgnom%y
1110 : end if
1111 33312 : end function get_gno_point
1112 :
1113 33312 : subroutine interpolation_point(nc,ns,gnom,gnom1d,face1,face2,xy,point,ida,ide,iref,ibaseref)
1114 : use coordinate_systems_mod, only : cartesian2D_t
1115 : implicit none
1116 : integer , intent(in) :: nc,ns
1117 : type (cartesian2D_t), intent(in) :: gnom
1118 : real (kind=r8), dimension(1-nc:), intent(in) :: gnom1d !dimension(1-nhc:nc+nhc)
1119 : integer, intent(in) :: face1, face2, xy
1120 : integer,intent(in) :: ida, ide
1121 : integer,intent(inout) :: iref,ibaseref
1122 : real (kind=r8), intent(inout) :: point
1123 :
1124 : ! type(cartesian3D_t) :: tmpcart3d
1125 : ! type (cartesian2D_t) :: tmpgnom
1126 :
1127 33312 : point = get_gno_point(gnom,face1,face2,xy)
1128 :
1129 : ! tmpcart3d=cubedsphere2cart(gnom,face1)
1130 : ! tmpgnom=cart2cubedsphere(tmpcart3d,face2)
1131 : ! if(xy==0) then
1132 : ! point=tmpgnom%x
1133 : ! else
1134 : ! point=tmpgnom%y
1135 : ! end if
1136 : !
1137 : ! in which cell is interpolation point located? gno(iref) is location of point to the right that is closest
1138 : !
1139 : ! |----------|---------|------x---|----------|------|------
1140 : ! gno(iref-1) gno(iref)
1141 : !
1142 33312 : iref=ida
1143 181584 : do while (point>gnom1d(iref))
1144 148320 : iref = iref + 1
1145 148320 : if (iref>ide+1) then
1146 0 : call endrun("error in search - ABORT; probably invalid ns-nc combination")
1147 : end if
1148 181584 : if (iref>ide) then
1149 : ! write(*,*) "extrapolation in interpolation_point",iref,ide
1150 48 : iref=ide
1151 48 : exit
1152 : endif
1153 : end do
1154 : !
1155 : ! this routine works for ns=1 and ns even
1156 : !
1157 33312 : if (MOD(ns,2)==1) then
1158 33312 : iref = max(iref,ida+1)!make sure gnom1d does not go out of bounds for extrapolation
1159 33312 : if (gnom1d(iref)-point>point-gnom1d(iref-1)) iref=iref-1
1160 33312 : iref=iref-((ns-1)/2)
1161 33312 : ibaseref = min(max(iref,ida),ide-(ns-1))!extrapolation
1162 33312 : point=point-gnom1d(ibaseref)
1163 0 : else if (MOD(ns, 2)==0) then
1164 : !
1165 : ! this code is only coded for ns even
1166 : !
1167 : ! ibaseref is the left most index used for 1D interpolation
1168 : ! (hence iref = iref-ns/2 except near corners)
1169 : !
1170 0 : iref = iref-ns/2
1171 0 : ibaseref = min(max(iref,ida),ide-(ns-1))
1172 0 : point=point-gnom1d(ibaseref)
1173 : end if
1174 33312 : end subroutine interpolation_point
1175 : !END SUBROUTINE INTERPOLATION_POINT---------------------------------------CE-for FVM!
1176 : ! ---------------------------------------------------------------------!
1177 : ! !
1178 : ! Precompute weights for Lagrange interpolation !
1179 : ! for equi-distant source grid values !
1180 : ! !
1181 : !----------------------------------------------------------------------!
1182 :
1183 33312 : subroutine get_equispace_weights(dx, x, w,ns)
1184 : !
1185 : ! Coordinate system for Lagrange interpolation:
1186 : !
1187 : ! |------|------|------|------|
1188 : ! 0 dx 2*dx 3*dx ns*dx
1189 : !
1190 : implicit none
1191 : real (kind=r8),intent(in) :: dx ! spacing of points, alpha/beta
1192 : real (kind=r8),intent(in) :: x ! X coordinate where interpolation is to be applied
1193 : real (kind=r8),dimension(:),intent(out) :: w ! dimension(ns)
1194 : integer ,intent(in) :: ns
1195 : !
1196 : integer :: j,k
1197 : !
1198 : ! use Lagrange interpolation formulae, e.g.,:
1199 : !
1200 : ! http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html
1201 : !
1202 133248 : w = 1.0_r8
1203 33312 : if (ns.ne.1) then
1204 133248 : do j=1,ns
1205 433056 : do k=1,ns
1206 399744 : if (k.ne.j) then
1207 199872 : w(j)=w(j)*(x-dble(k-1)*dx)/(dble(j-1)*dx-dble(k-1)*dx)
1208 : end if
1209 : end do
1210 : end do
1211 : end if
1212 33312 : end subroutine get_equispace_weights
1213 :
1214 : end module fvm_analytic_mod
|
