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 118 of file mod_mkinit_util.F90.

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

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

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 348 of file mod_mkinit_util.F90.

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

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

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

Generated by doxygen 1.10.0