FE-Project
Loading...
Searching...
No Matches
Functions/Subroutines
mod_mkinit_util Module Reference

module Utility for mkinit More...

Functions/Subroutines

subroutine, public mkinitutil_gen_gpmat (gpmat, elem_intrp, elem)
subroutine, public mkinitutil_gen_vm1mat (vm1mat, elem_intrp, elem)
subroutine, public mkinitutil_calc_cosinebell (q, qmax, rx, ry, rz, xc, yc, zc, x, y, z, lcmesh3d, elem, intrppolyorder_h, intrppolyorder_v, z_func_type, z_func_params, cosbell_exponent)
 Calculate the distribution function of a cosine bell in regional domain.
subroutine, public mkinitutil_calc_cosinebell_global (q, qmax, rh, lonc, latc, rplanet, x, y, z, lcmesh3d, elem, intrppolyorder_h, intrppolyorder_v, z_func_type, z_func_params, cosbell_exponent)
 Calculate the distribution function of a cosine bell in global domain.
subroutine, public mkinitutil_galerkinprojection (q, func, intrppolyorder_h, intrppolyorder_v, lcmesh3d, elem)
 Apply the Galerkin projection to the user-defined function.
subroutine, public mkinitutil_galerkinprojection_global (q, func, intrppolyorder_h, intrppolyorder_v, lcmesh3d, elem, rplanet)
 Apply the Galerkin projection to the user-defined function (for global model)

Detailed Description

module Utility for mkinit

Description
subroutines useful to prepare initial data
Author
Yuta Kawai, Team SCALE

Function/Subroutine Documentation

◆ mkinitutil_gen_gpmat()

subroutine, public mod_mkinit_util::mkinitutil_gen_gpmat ( real(rp), dimension(elem%np,elem_intrp%np), intent(out) gpmat,
class(elementbase3d), intent(in) elem_intrp,
class(elementbase3d), intent(in) elem )

Definition at line 51 of file mod_mkinit_util.F90.

53 implicit none
54
55 class(ElementBase3D), intent(in) :: elem_intrp
56 class(ElementBase3D), intent(in) :: elem
57 real(RP), intent(out) :: GPMat(elem%Np,elem_intrp%Np)
58
59 integer :: p1, p2, p3, p_
60 integer :: p_intrp
61
62 real(RP) :: InvV_intrp(elem%Np,elem_intrp%Np)
63 !---------------------------------------------
64
65 invv_intrp(:,:) = 0.0_rp
66 do p3=1, elem%PolyOrder_v+1
67 do p2=1, elem%PolyOrder_h+1
68 do p1=1, elem%PolyOrder_h+1
69 p_ = p1 + (p2-1)*(elem%PolyOrder_h + 1) + (p3-1)*(elem%PolyOrder_h + 1)**2
70 p_intrp = p1 + (p2-1)*(elem_intrp%PolyOrder_h + 1) + (p3-1)*(elem_intrp%PolyOrder_h + 1)**2
71 invv_intrp(p_,:) = elem_intrp%invV(p_intrp,:)
72 end do
73 end do
74 end do
75 gpmat(:,:) = matmul(elem%V, invv_intrp)
76
77 return

Referenced by mkinitutil_calc_cosinebell(), mkinitutil_calc_cosinebell_global(), mkinitutil_galerkinprojection(), and mkinitutil_galerkinprojection_global().

◆ mkinitutil_gen_vm1mat()

subroutine, public mod_mkinit_util::mkinitutil_gen_vm1mat ( real(rp), dimension(elem%np,elem_intrp%np), intent(out) vm1mat,
class(elementbase3d), intent(in) elem_intrp,
class(elementbase3d), intent(in) elem )

Definition at line 81 of file mod_mkinit_util.F90.

83 implicit none
84
85 class(ElementBase3D), intent(in) :: elem_intrp
86 class(ElementBase3D), intent(in) :: elem
87 real(RP), intent(out) :: Vm1Mat(elem%Np,elem_intrp%Np)
88
89 integer :: p1, p2, p3, p_
90 integer :: p_intrp
91
92 real(RP) :: InvV_intrpVm1(elem%Np,elem_intrp%Np)
93 !---------------------------------------------
94
95 invv_intrpvm1(:,:) = 0.0_rp
96 do p3=1, elem%PolyOrder_v
97 do p2=1, elem%PolyOrder_h+1
98 do p1=1, elem%PolyOrder_h+1
99 p_ = p1 + (p2-1)*(elem%PolyOrder_h + 1) + (p3-1)*(elem%PolyOrder_h + 1)**2
100 p_intrp = p1 + (p2-1)*(elem_intrp%PolyOrder_h + 1) + (p3-1)*(elem_intrp%PolyOrder_h + 1)**2
101 invv_intrpvm1(p_,:) = elem_intrp%invV(p_intrp,:)
102 end do
103 end do
104 end do
105 vm1mat(:,:) = matmul(elem%V, invv_intrpvm1)
106
107 return

◆ mkinitutil_calc_cosinebell()

subroutine, public mod_mkinit_util::mkinitutil_calc_cosinebell ( real(rp), dimension(elem%np,lcmesh3d%nea), intent(out) q,
real(rp), intent(in) qmax,
real(rp), intent(in) rx,
real(rp), intent(in) ry,
real(rp), intent(in) rz,
real(rp), intent(in) xc,
real(rp), intent(in) yc,
real(rp), intent(in) zc,
real(rp), dimension(elem%np,lcmesh3d%ne), intent(in) x,
real(rp), dimension(elem%np,lcmesh3d%ne), intent(in) y,
real(rp), dimension(elem%np,lcmesh3d%ne), intent(in) z,
class(localmesh3d), intent(in) lcmesh3d,
class(elementbase3d), intent(in) elem,
integer, intent(in) intrppolyorder_h,
integer, intent(in) intrppolyorder_v,
character(len=*), intent(in), optional z_func_type,
real(rp), dimension(:), intent(in), optional z_func_params,
integer, intent(in), optional cosbell_exponent )

Calculate the distribution function of a cosine bell in regional domain.

If the vertical dependence is considered, specify z_func_type and z_func_params. For z_func_type = 'sin', the values of z_func_params is 1: the vertical model, 2: the half of wavelength

Parameters
[in]cosbell_exponentparameter to ensure 2*cosbell_exponent-1 continuous derivatives

Definition at line 117 of file mod_mkinit_util.F90.

123
124 implicit none
125
126 class(LocalMesh3D), intent(in) :: lcmesh3D
127 class(ElementBase3D), intent(in) :: elem
128 real(RP), intent(out) :: q(elem%Np,lcmesh3D%NeA)
129 real(RP), intent(in) :: qmax
130 real(RP), intent(in) :: rx, ry, rz
131 real(RP), intent(in) :: xc, yc, zc
132 real(RP), intent(in) :: x(elem%Np,lcmesh3D%Ne)
133 real(RP), intent(in) :: y(elem%Np,lcmesh3D%Ne)
134 real(RP), intent(in) :: z(elem%Np,lcmesh3D%Ne)
135 integer, intent(in) :: IntrpPolyOrder_h
136 integer, intent(in) :: IntrpPolyOrder_v
137 character(len=*), optional, intent(in) :: z_func_type
138 real(RP), optional, intent(in) :: z_func_params(:)
139 integer, intent(in), optional :: cosbell_exponent !< parameter to ensure 2*cosbell_exponent-1 continuous derivatives
140
141 integer :: ke
142
143 type(HexahedralElement) :: elem_intrp
144 real(RP), allocatable :: x_intrp(:,:), y_intrp(:,:), z_intrp(:,:)
145 real(RP), allocatable :: r_intrp(:)
146 real(RP), allocatable :: z_func(:,:)
147 real(RP) :: vx(elem%Nv), vy(elem%Nv), vz(elem%Nv)
148
149 real(RP), allocatable :: IntrpMat(:,:)
150 real(RP), allocatable :: q_intrp(:)
151
152 integer :: exponent
153 !-----------------------------------------------
154
155 if ( present(cosbell_exponent) ) then
156 exponent = cosbell_exponent
157 else
158 exponent = 1
159 end if
160
161 call elem_intrp%Init( intrppolyorder_h, intrppolyorder_v, .false. )
162
163 allocate( intrpmat(elem%Np,elem_intrp%Np) )
164 call mkinitutil_gen_gpmat( intrpmat, elem_intrp, elem )
165
166 allocate( x_intrp(elem_intrp%Np,lcmesh3d%Ne), y_intrp(elem_intrp%Np,lcmesh3d%Ne), z_intrp(elem_intrp%Np,lcmesh3d%Ne) )
167 allocate( z_func(elem_intrp%Np,lcmesh3d%Ne))
168 allocate( r_intrp(elem_intrp%Np) )
169 allocate( q_intrp(elem_intrp%Np) )
170
171 !$omp parallel private( &
172 !$omp q_intrp, vx, vy, vz, &
173 !$omp x_intrp, y_intrp, z_intrp, r_intrp )
174
175 !$omp do
176 do ke=lcmesh3d%NeS, lcmesh3d%NeE
177
178 vx(:) = lcmesh3d%pos_ev(lcmesh3d%EToV(ke,:),1)
179 vy(:) = lcmesh3d%pos_ev(lcmesh3d%EToV(ke,:),2)
180 vz(:) = lcmesh3d%pos_ev(lcmesh3d%EToV(ke,:),3)
181 x_intrp(:,ke) = vx(1) + 0.5_rp * ( elem_intrp%x1(:) + 1.0_rp ) * ( vx(2) - vx(1) )
182 y_intrp(:,ke) = vy(1) + 0.5_rp * ( elem_intrp%x2(:) + 1.0_rp ) * ( vy(4) - vy(1) )
183 z_intrp(:,ke) = vz(1) + 0.5_rp * ( elem_intrp%x3(:) + 1.0_rp ) * ( vz(5) - vz(1) )
184
185 z_func(:,ke) = 1.0_rp
186 end do
187 !$omp end do
188
189 ! Calculate the vertical function
190 if ( present(z_func_type) ) then
191 select case(z_func_type)
192 case ('sin')
193 !$omp do
194 do ke=lcmesh3d%NeS, lcmesh3d%NeE
195 z_func(:,ke) = sin( z_func_params(1) * pi * z_intrp(:,ke) / z_func_params(2) )
196 end do
197 !$omp end do
198 end select
199 end if
200
201 !$omp do
202 do ke=lcmesh3d%NeS, lcmesh3d%NeE
203 r_intrp(:) = sqrt( &
204 ( (x_intrp(:,ke) - xc) / rx )**2 &
205 + ( (y_intrp(:,ke) - yc) / ry )**2 &
206 + ( (z_intrp(:,ke) - zc) / rz )**2 )
207
208 where( r_intrp(:) <= 1.0_rp )
209 q_intrp(:) = qmax * ( 0.5_rp * (1.0_rp + cos( pi * r_intrp(:) ) ) )**exponent
210 elsewhere
211 q_intrp(:) = 0.0_rp
212 end where
213 q(:,ke) = matmul(intrpmat, q_intrp(:) * z_func(:,ke))
214 end do
215 !$omp end do
216 !$omp end parallel
217
218 call elem_intrp%Final()
219
220 return

References mkinitutil_gen_gpmat().

◆ mkinitutil_calc_cosinebell_global()

subroutine, public mod_mkinit_util::mkinitutil_calc_cosinebell_global ( real(rp), dimension(elem%np,lcmesh3d%nea), intent(out) q,
real(rp), intent(in) qmax,
real(rp), intent(in) rh,
real(rp), intent(in) lonc,
real(rp), intent(in) latc,
real(rp), intent(in) rplanet,
real(rp), dimension(elem%np,lcmesh3d%ne), intent(in) x,
real(rp), dimension(elem%np,lcmesh3d%ne), intent(in) y,
real(rp), dimension(elem%np,lcmesh3d%ne), intent(in) z,
class(localmesh3d), intent(in) lcmesh3d,
class(elementbase3d), intent(in) elem,
integer, intent(in) intrppolyorder_h,
integer, intent(in) intrppolyorder_v,
character(len=*), intent(in), optional z_func_type,
real(rp), dimension(:), intent(in), optional z_func_params,
integer, intent(in), optional cosbell_exponent )

Calculate the distribution function of a cosine bell in global domain.

If the vertical dependence is considered, specify z_func_type and z_func_params. For z_func_type = 'sin', the values of z_func_params is 1: the vertical model, 2: the half of wavelength

Parameters
[in]cosbell_exponentparameter to ensure 2*cosbell_exponent-1 continuous derivatives

Definition at line 230 of file mod_mkinit_util.F90.

236
237 use scale_cubedsphere_coord_cnv, only: &
238 cubedspherecoordcnv_cs2lonlatpos
239 implicit none
240
241 class(LocalMesh3D), intent(in) :: lcmesh3D
242 class(ElementBase3D), intent(in) :: elem
243 real(RP), intent(out) :: q(elem%Np,lcmesh3D%NeA)
244 real(RP), intent(in) :: qmax
245 real(RP), intent(in) :: rh
246 real(RP), intent(in) :: lonc, latc
247 real(RP), intent(in) :: rplanet
248 real(RP), intent(in) :: x(elem%Np,lcmesh3D%Ne)
249 real(RP), intent(in) :: y(elem%Np,lcmesh3D%Ne)
250 real(RP), intent(in) :: z(elem%Np,lcmesh3D%Ne)
251 integer, intent(in) :: IntrpPolyOrder_h
252 integer, intent(in) :: IntrpPolyOrder_v
253 character(len=*), optional, intent(in) :: z_func_type
254 real(RP), optional, intent(in) :: z_func_params(:)
255 integer, intent(in), optional :: cosbell_exponent !< parameter to ensure 2*cosbell_exponent-1 continuous derivatives
256
257 integer :: ke
258
259 type(HexahedralElement) :: elem_intrp
260 real(RP), allocatable :: x_intrp(:,:), y_intrp(:,:), z_intrp(:,:)
261 real(RP), allocatable :: gam_intrp(:,:)
262 real(RP), allocatable :: lon_intrp(:,:), lat_intrp(:,:)
263 real(RP), allocatable :: z_func(:,:)
264 real(RP), allocatable :: r_intrp(:)
265 real(RP) :: vx(elem%Nv), vy(elem%Nv), vz(elem%Nv)
266
267 real(RP), allocatable :: IntrpMat(:,:)
268 real(RP), allocatable :: q_intrp(:)
269
270 integer :: exponent
271
272 !-----------------------------------------------
273
274 if ( present(cosbell_exponent) ) then
275 exponent = cosbell_exponent
276 else
277 exponent = 1
278 end if
279
280 call elem_intrp%Init( intrppolyorder_h, intrppolyorder_v, .false. )
281
282 allocate( intrpmat(elem%Np,elem_intrp%Np) )
283 call mkinitutil_gen_gpmat( intrpmat, elem_intrp, elem )
284
285 allocate( x_intrp(elem_intrp%Np,lcmesh3d%Ne), y_intrp(elem_intrp%Np,lcmesh3d%Ne), z_intrp(elem_intrp%Np,lcmesh3d%Ne) )
286 allocate( gam_intrp(elem_intrp%Np,lcmesh3d%Ne) )
287 allocate( lon_intrp(elem_intrp%Np,lcmesh3d%Ne), lat_intrp(elem_intrp%Np,lcmesh3d%Ne) )
288 allocate( z_func(elem_intrp%Np,lcmesh3d%Ne) )
289 allocate( r_intrp(elem_intrp%Np) )
290 allocate( q_intrp(elem_intrp%Np) )
291
292
293 !$omp parallel do private(vx, vy, vz)
294 do ke=lcmesh3d%NeS, lcmesh3d%NeE
295 vx(:) = lcmesh3d%pos_ev(lcmesh3d%EToV(ke,:),1)
296 vy(:) = lcmesh3d%pos_ev(lcmesh3d%EToV(ke,:),2)
297 vz(:) = lcmesh3d%pos_ev(lcmesh3d%EToV(ke,:),3)
298 x_intrp(:,ke) = vx(1) + 0.5_rp * ( elem_intrp%x1(:) + 1.0_rp ) * ( vx(2) - vx(1) )
299 y_intrp(:,ke) = vy(1) + 0.5_rp * ( elem_intrp%x2(:) + 1.0_rp ) * ( vy(4) - vy(1) )
300 z_intrp(:,ke) = vz(1) + 0.5_rp * ( elem_intrp%x3(:) + 1.0_rp ) * ( vz(5) - vz(1) )
301 gam_intrp(:,ke) = 1.0_rp
302
303 z_func(:,ke) = 1.0_rp
304 end do
305
306 call cubedspherecoordcnv_cs2lonlatpos( lcmesh3d%panelID, x_intrp, y_intrp, gam_intrp, & ! (in)
307 elem_intrp%Np * lcmesh3d%Ne, & ! (in)
308 lon_intrp(:,:), lat_intrp(:,:) ) ! (out)
309
310 ! Calculate the vertical function
311 if ( present(z_func_type) ) then
312 select case(z_func_type)
313 case ('sin')
314 !$omp parallel do
315 do ke=lcmesh3d%NeS, lcmesh3d%NeE
316 z_func(:,ke) = sin( z_func_params(1) * pi * z_intrp(:,ke) / z_func_params(2) )
317 end do
318 end select
319 end if
320
321 !$omp parallel do private( r_intrp, q_intrp )
322 do ke=lcmesh3d%NeS, lcmesh3d%NeE
323
324 ! Calculate the horizontal function
325 r_intrp(:) = rplanet / rh * acos( sin(latc) * sin(lat_intrp(:,ke)) + cos(latc) * cos(lat_intrp(:,ke)) * cos(lon_intrp(:,ke) - lonc) )
326 where( r_intrp(:) <= 1.0_rp )
327 q_intrp(:) = qmax * ( 0.5_rp * (1.0_rp + cos( pi * r_intrp(:) ) ) )**exponent
328 elsewhere
329 q_intrp(:) = 0.0_rp
330 end where
331
332 ! Perform Galerkin projection
333 q(:,ke) = matmul(intrpmat, q_intrp(:) * z_func(:,ke))
334 end do
335
336 call elem_intrp%Final()
337
338 return
scale_cubedsphere_coord_cnv
Module common / Coordinate conversion with cubed-sphere projection.
Definition scale_cubedsphere_coord_cnv.F90:21
scale_cubedsphere_coord_cnv::cubedspherecoordcnv_cs2lonlatpos
subroutine, public cubedspherecoordcnv_cs2lonlatpos(panelid, alpha, beta, gam, np, lon, lat)
Calculate longitude and latitude coordinates from local coordinates using the central angles in an eq...
Definition scale_cubedsphere_coord_cnv.F90:88

References scale_cubedsphere_coord_cnv::cubedspherecoordcnv_cs2lonlatpos(), and mkinitutil_gen_gpmat().

◆ mkinitutil_galerkinprojection()

subroutine, public mod_mkinit_util::mkinitutil_galerkinprojection ( real(rp), dimension(elem%np,lcmesh3d%nea), intent(out) q,
external subroutine(real(rp), dimension(elem_intrp%np), intent(out) q_intrp, real(rp), dimension(elem_intrp%np), intent(in) x, real(rp), dimension(elem_intrp%np), intent(in) y, real(rp), dimension(elem_intrp%np), intent(in) z, class(elementbase3d), intent(in) elem_intrp) func,
integer, intent(in) intrppolyorder_h,
integer, intent(in) intrppolyorder_v,
class(localmesh3d), intent(in) lcmesh3d,
class(elementbase3d), intent(in) elem )

Apply the Galerkin projection to the user-defined function.

Definition at line 345 of file mod_mkinit_util.F90.

348
349 implicit none
350 class(LocalMesh3D), intent(in) :: lcmesh3D
351 class(ElementBase3D), intent(in) :: elem
352 real(RP), intent(out) :: q(elem%Np,lcmesh3D%NeA)
353 integer, intent(in) :: IntrpPolyOrder_h
354 integer, intent(in) :: IntrpPolyOrder_v
355
356 interface
357 subroutine func( q_intrp, &
358 x, y, z, elem_intrp )
359 import elementbase3d
360 import rp
361 class(ElementBase3D), intent(in) :: elem_intrp
362 real(RP), intent(out) :: q_intrp(elem_intrp%Np)
363 real(RP), intent(in) :: x(elem_intrp%Np)
364 real(RP), intent(in) :: y(elem_intrp%Np)
365 real(RP), intent(in) :: z(elem_intrp%Np)
366 end subroutine func
367 end interface
368
369 type(HexahedralElement) :: elem_intrp
370 real(RP), allocatable :: x_intrp(:,:), y_intrp(:,:), z_intrp(:,:)
371 real(RP) :: vx(elem%Nv), vy(elem%Nv), vz(elem%Nv)
372
373 real(RP), allocatable :: IntrpMat(:,:)
374 real(RP), allocatable :: q_intrp(:)
375
376 integer :: ke
377 !-----------------------------------------------
378
379 call elem_intrp%Init( intrppolyorder_h, intrppolyorder_v, .false. )
380
381 allocate( intrpmat(elem%Np,elem_intrp%Np) )
382 call mkinitutil_gen_gpmat( intrpmat, elem_intrp, elem )
383
384 allocate( x_intrp(elem_intrp%Np,lcmesh3d%Ne), y_intrp(elem_intrp%Np,lcmesh3d%Ne), z_intrp(elem_intrp%Np,lcmesh3d%Ne) )
385 allocate( q_intrp(elem_intrp%Np) )
386
387 !$omp parallel do private(vx, vy, vz)
388 do ke=lcmesh3d%NeS, lcmesh3d%NeE
389 vx(:) = lcmesh3d%pos_ev(lcmesh3d%EToV(ke,:),1)
390 vy(:) = lcmesh3d%pos_ev(lcmesh3d%EToV(ke,:),2)
391 vz(:) = lcmesh3d%pos_ev(lcmesh3d%EToV(ke,:),3)
392 x_intrp(:,ke) = vx(1) + 0.5_rp * ( elem_intrp%x1(:) + 1.0_rp ) * ( vx(2) - vx(1) )
393 y_intrp(:,ke) = vy(1) + 0.5_rp * ( elem_intrp%x2(:) + 1.0_rp ) * ( vy(4) - vy(1) )
394 z_intrp(:,ke) = vz(1) + 0.5_rp * ( elem_intrp%x3(:) + 1.0_rp ) * ( vz(5) - vz(1) )
395 end do
396
397 !$omp parallel do private( q_intrp )
398 do ke=lcmesh3d%NeS, lcmesh3d%NeE
399
400 call func( q_intrp, & ! (out)
401 x_intrp(:,ke), y_intrp(:,ke), z_intrp(:,ke), & ! (in)
402 elem_intrp ) ! (in)
403
404 ! Perform Galerkin projection
405 q(:,ke) = matmul( intrpmat, q_intrp )
406 end do
407
408 call elem_intrp%Final()
409
410 return

References mkinitutil_gen_gpmat().

◆ mkinitutil_galerkinprojection_global()

subroutine, public mod_mkinit_util::mkinitutil_galerkinprojection_global ( real(rp), dimension(elem%np,lcmesh3d%nea), intent(out) q,
external subroutine(real(rp), dimension(elem_intrp%np), intent(out) q_intrp, real(rp), dimension(elem_intrp%np), intent(in) lon, real(rp), dimension(elem_intrp%np), intent(in) lat, real(rp), dimension(elem_intrp%np), intent(in) z, class(elementbase3d), intent(in) elem_intrp, real(rp), intent(in) rplanet_) func,
integer, intent(in) intrppolyorder_h,
integer, intent(in) intrppolyorder_v,
class(localmesh3d), intent(in) lcmesh3d,
class(elementbase3d), intent(in) elem,
real(rp), intent(in) rplanet )

Apply the Galerkin projection to the user-defined function (for global model)

Definition at line 416 of file mod_mkinit_util.F90.

419
420 use scale_cubedsphere_coord_cnv, only: &
421 cubedspherecoordcnv_cs2lonlatpos
422
423 implicit none
424 class(LocalMesh3D), intent(in) :: lcmesh3D
425 class(ElementBase3D), intent(in) :: elem
426 real(RP), intent(out) :: q(elem%Np,lcmesh3D%NeA)
427 integer, intent(in) :: IntrpPolyOrder_h
428 integer, intent(in) :: IntrpPolyOrder_v
429 real(RP), intent(in) :: rplanet
430
431 interface
432 subroutine func( q_intrp, &
433 lon, lat, z, elem_intrp, rplanet_ )
434 import elementbase3d
435 import rp
436 class(ElementBase3D), intent(in) :: elem_intrp
437 real(RP), intent(out) :: q_intrp(elem_intrp%Np)
438 real(RP), intent(in) :: lon(elem_intrp%Np)
439 real(RP), intent(in) :: lat(elem_intrp%Np)
440 real(RP), intent(in) :: z(elem_intrp%Np)
441 real(RP), intent(in) :: rplanet_
442 end subroutine func
443 end interface
444
445 type(HexahedralElement) :: elem_intrp
446 real(RP), allocatable :: x_intrp(:,:), y_intrp(:,:), z_intrp(:,:)
447 real(RP), allocatable :: gam_intrp(:,:)
448 real(RP), allocatable :: lon_intrp(:,:), lat_intrp(:,:)
449 real(RP) :: vx(elem%Nv), vy(elem%Nv), vz(elem%Nv)
450
451 real(RP), allocatable :: IntrpMat(:,:)
452 real(RP), allocatable :: q_intrp(:)
453
454 integer :: ke
455 !-----------------------------------------------
456
457 call elem_intrp%Init( intrppolyorder_h, intrppolyorder_v, .false. )
458
459 allocate( intrpmat(elem%Np,elem_intrp%Np) )
460 call mkinitutil_gen_gpmat( intrpmat, elem_intrp, elem )
461
462 allocate( x_intrp(elem_intrp%Np,lcmesh3d%Ne), y_intrp(elem_intrp%Np,lcmesh3d%Ne), z_intrp(elem_intrp%Np,lcmesh3d%Ne) )
463 allocate( gam_intrp(elem_intrp%Np,lcmesh3d%Ne) )
464 allocate( lon_intrp(elem_intrp%Np,lcmesh3d%Ne), lat_intrp(elem_intrp%Np,lcmesh3d%Ne) )
465 allocate( q_intrp(elem_intrp%Np) )
466
467 !$omp parallel do private(vx, vy, vz)
468 do ke=lcmesh3d%NeS, lcmesh3d%NeE
469 vx(:) = lcmesh3d%pos_ev(lcmesh3d%EToV(ke,:),1)
470 vy(:) = lcmesh3d%pos_ev(lcmesh3d%EToV(ke,:),2)
471 vz(:) = lcmesh3d%pos_ev(lcmesh3d%EToV(ke,:),3)
472 x_intrp(:,ke) = vx(1) + 0.5_rp * ( elem_intrp%x1(:) + 1.0_rp ) * ( vx(2) - vx(1) )
473 y_intrp(:,ke) = vy(1) + 0.5_rp * ( elem_intrp%x2(:) + 1.0_rp ) * ( vy(4) - vy(1) )
474 z_intrp(:,ke) = vz(1) + 0.5_rp * ( elem_intrp%x3(:) + 1.0_rp ) * ( vz(5) - vz(1) )
475
476 gam_intrp(:,ke) = 1.0_rp
477 end do
478
479 call cubedspherecoordcnv_cs2lonlatpos( lcmesh3d%panelID, x_intrp, y_intrp, gam_intrp, & ! (in)
480 elem_intrp%Np * lcmesh3d%Ne, & ! (in)
481 lon_intrp(:,:), lat_intrp(:,:) ) ! (out)
482
483 !$omp parallel do private( q_intrp )
484 do ke=lcmesh3d%NeS, lcmesh3d%NeE
485
486 call func( q_intrp, & ! (out)
487 lon_intrp(:,ke), lat_intrp(:,ke), z_intrp(:,ke), & ! (in)
488 elem_intrp, rplanet ) ! (in)
489
490 ! Perform Galerkin projection
491 q(:,ke) = matmul( intrpmat, q_intrp )
492 end do
493
494 call elem_intrp%Final()
495
496 return

References scale_cubedsphere_coord_cnv::cubedspherecoordcnv_cs2lonlatpos(), and mkinitutil_gen_gpmat().

Generated by doxygen 1.14.0