FE-Project
Loading...
Searching...
No Matches
scale_atm_dyn_dgm_nonhydro3d_numdiff.F90
Go to the documentation of this file.
1!-------------------------------------------------------------------------------
2!> module Atmosphere / Dynamics common
3!!
4!! @par Description
5!! Explicit numerical diffusion for Atmospheric dynamical process.
6!! For the discretization, the local DGM (e.g., Cockburn and Shu, 1998) is used.
7!!
8!! @author Yuta Kawai, Team SCALE
9!<
10!-------------------------------------------------------------------------------
11#include "scaleFElib.h"
12module scale_atm_dyn_dgm_nonhydro3d_numdiff
13 !-----------------------------------------------------------------------------
14 !
15 !++ Used modules
16 !
17 use scale_precision
18 use scale_io
19 use scale_prc
20 use scale_prof
21
22 use scale_sparsemat
23 use scale_element_base, only: &
24 elementbase2d, elementbase3d
25 use scale_element_hexahedral, only: hexahedralelement
26 use scale_localmesh_2d, only: localmesh2d
27 use scale_localmesh_3d, only: localmesh3d
28 use scale_mesh_base, only: meshbase
29 use scale_mesh_base3d, only: meshbase3d
30 use scale_localmeshfield_base, only: localmeshfield3d
31 use scale_meshfield_base, only: meshfield3d
32
33 use scale_model_var_manager, only: &
34 modelvarmanager, variableinfo
35 use scale_model_mesh_manager, only: modelmesh3d
36
37 use scale_atm_dyn_dgm_nonhydro3d_common, only: &
38 prgvar_num, &
39 dens_vid => prgvar_ddens_id, therm_vid => prgvar_therm_id,&
40 momx_vid => prgvar_momx_id, momy_vid => prgvar_momy_id, &
41 momz_vid => prgvar_momz_id, denshyd_vid => auxvar_denshydro_id
42
43 use scale_atm_dyn_dgm_bnd, only: atmdynbnd
44
45 !-----------------------------------------------------------------------------
46 implicit none
47 private
48 !-----------------------------------------------------------------------------
49 !
50 !++ Public procedures
51 !
52
53 type, public :: atmdyn_nonhydro3d_numdiff
54 integer :: nd_laplacian_num
55 real(rp) :: nd_coef_h
56 real(rp) :: nd_coef_v
57 real(rp) :: dtsec
58
59 type(meshfield3d), allocatable :: numdiff_flux_vars3d(:)
60 type(modelvarmanager) :: numdiff_flux_manager
61 integer :: numdiff_flux_commid
62
63 type(meshfield3d), allocatable :: numdiff_tend_vars3d(:)
64 type(modelvarmanager) :: numdiff_tend_manager
65 integer :: numdiff_tend_commid
66
67 contains
68 procedure :: init => atm_dyn_dgm_nonhydro3d_numdiff_init
69 procedure :: final => atm_dyn_dgm_nonhydro3d_numdiff_final
70 procedure :: apply => atm_dyn_dgm_nonhydro3d_numdiff_apply
71 end type atmdyn_nonhydro3d_numdiff
72
73 !-----------------------------------------------------------------------------
74 !
75 !++ Public parameters & variables
76 !
77 integer, public, parameter :: numdiff_flux_num = 3
78 integer, public, parameter :: numdiffflx_x_id = 1
79 integer, public, parameter :: numdiffflx_y_id = 2
80 integer, public, parameter :: numdiffflx_z_id = 3
81
82 type(variableinfo), public :: atmos_dyn_numdiff_flux_vinfo(numdiff_flux_num)
83 DATA atmos_dyn_numdiff_flux_vinfo / &
84 variableinfo( numdiffflx_x_id, 'DIFFFLX_X', 'flux in x-direction', &
85 '?.m/s', 3, 'XYZ', '' ), &
86 variableinfo( numdiffflx_y_id, 'DIFFFLX_Y', 'flux in y-direction', &
87 '?.m/s', 3, 'XYZ', '' ), &
88 variableinfo( numdiffflx_z_id, 'DIFFFLX_Z', 'flux in z-direction', &
89 '?.m/s', 3, 'XYZ', '' ) /
90
91 !-
92 integer, public, parameter :: numdiff_tend_num = 2
93 integer, public, parameter :: numdiff_laplah_id = 1
94 integer, public, parameter :: numdiff_laplav_id = 2
95
96 type(variableinfo), public :: atmos_dyn_numdiff_tend_vinfo(numdiff_tend_num)
97 DATA atmos_dyn_numdiff_tend_vinfo / &
98 variableinfo( numdiff_laplah_id, 'NUMDIFF_LAPLAH', 'tendency due to nundiff', &
99 '?/s', 3, 'XYZ', '' ), &
100 variableinfo( numdiff_laplav_id, 'NUMDIFF_LAPLAV', 'tendency due to nundiff', &
101 '?/s', 3, 'XYZ', '' ) /
102
103 !-----------------------------------------------------------------------------
104 !
105 !++ Private procedures & variables
106 !
107 !-------------------
108
109 private :: numdiff_tend
110 private :: numdiff_cal_laplacian
111 private :: numdiff_cal_flx
112
113 private :: cal_del_flux_lap
114 private :: cal_del_flux_lap_with_coef
115 private :: cal_del_graddiffvar
116
117contains
118!OCL SERIAL
119 subroutine atm_dyn_dgm_nonhydro3d_numdiff_init( this, model_mesh3D, dtsec )
120 use scale_prc, only: prc_abort
121 implicit none
122 class(atmdyn_nonhydro3d_numdiff), intent(inout) :: this
123 class(modelmesh3d), intent(inout), target :: model_mesh3d
124 real(rp), intent(in) :: dtsec
125
126 class(meshbase3d), pointer :: mesh3d
127 !--------------------------------------------
128
129 integer :: nd_laplacian_num = 1
130 real(rp) :: nd_coef_h = 0.0_rp
131 real(rp) :: nd_coef_v = 0.0_rp
132
133 namelist /param_atmos_dyn_numdiff/ &
134 nd_laplacian_num, &
135 nd_coef_h, nd_coef_v
136
137 integer :: ierr
138
139 integer :: v
140 !---------------------------------------------------------------
141
142 rewind(io_fid_conf)
143 read(io_fid_conf,nml=param_atmos_dyn_numdiff,iostat=ierr)
144 if( ierr < 0 ) then !--- missing
145 log_info("ATMOS_DYN_setup_numdiff",*) 'Not found namelist. Default used.'
146 elseif( ierr > 0 ) then !--- fatal error
147 log_error("ATMOS_DYN_setup_numdiff",*) 'Not appropriate names in namelist PARAM_ATMOS_DYN_NUMDIFF. Check!'
148 call prc_abort
149 endif
150 log_nml(param_atmos_dyn_numdiff)
151
152 this%ND_LAPLACIAN_NUM = nd_laplacian_num
153 this%ND_COEF_H = nd_coef_h
154 this%ND_COEF_v = nd_coef_v
155
156 this%dtsec = dtsec
157
158 !------------
159 mesh3d => model_mesh3d%ptr_mesh
160
161 call this%NUMDIFF_FLUX_manager%Init()
162 allocate( this%NUMDIFF_FLUX_VARS3D(numdiff_flux_num) )
163
164 do v = 1, numdiff_flux_num
165 call this%NUMDIFF_FLUX_manager%Regist( &
166 atmos_dyn_numdiff_flux_vinfo(v), mesh3d, & ! (in)
167 this%NUMDIFF_FLUX_VARS3D(v), & ! (inout)
168 .false., fill_zero=.true. ) ! (in)
169 end do
170
171 call model_mesh3d%Create_communicator( &
172 numdiff_flux_num, 0, 0, & ! (in)
173 this%NUMDIFF_FLUX_manager, & ! (inout)
174 this%NUMDIFF_FLUX_VARS3D(:), & ! (in)
175 this%NUMDIFF_FLUX_commid ) ! (out)
176
177 !-
178 call this%NUMDIFF_TEND_manager%Init()
179 allocate( this%NUMDIFF_TEND_VARS3D(numdiff_tend_num) )
180
181 do v = 1, numdiff_tend_num
182 call this%NUMDIFF_TEND_manager%Regist( &
183 atmos_dyn_numdiff_tend_vinfo(v), mesh3d, & ! (in)
184 this%NUMDIFF_TEND_VARS3D(v), & ! (inout)
185 .false., fill_zero=.true. ) ! (in)
186 end do
187
188 call model_mesh3d%Create_communicator( &
189 numdiff_tend_num, 0, 0, & ! (in)
190 this%NUMDIFF_TEND_manager, & ! (inout)
191 this%NUMDIFF_TEND_VARS3D(:), & ! (in)
192 this%NUMDIFF_TEND_commid ) ! (out)
193
194 return
195 end subroutine atm_dyn_dgm_nonhydro3d_numdiff_init
196
197!OCL SERIAL
198 subroutine atm_dyn_dgm_nonhydro3d_numdiff_final( this )
199 implicit none
200
201 class(atmdyn_nonhydro3d_numdiff), intent(inout) :: this
202 !--------------------------------------------
203
204 call this%NUMDIFF_FLUX_manager%Final()
205 deallocate( this%NUMDIFF_FLUX_VARS3D )
206
207 call this%NUMDIFF_TEND_manager%Final()
208 deallocate( this%NUMDIFF_TEND_VARS3D )
209
210 return
211 end subroutine atm_dyn_dgm_nonhydro3d_numdiff_final
212
213!OCL SERIAL
214 subroutine atm_dyn_dgm_nonhydro3d_numdiff_apply( this, &
215 PROG_VARS, AUX_VARS, boundary_cond, &
216 Dx, Dy, Dz, Lift, mesh )
217
218 implicit none
219
220 class(atmdyn_nonhydro3d_numdiff), intent(inout) :: this
221 class(modelvarmanager), intent(inout) :: prog_vars
222 class(modelvarmanager), intent(inout) :: aux_vars
223 class(atmdynbnd), intent(in) :: boundary_cond
224 type(sparsemat), intent(in) :: dx, dy, dz, lift
225 class(meshbase3d), intent(in), target :: mesh
226
227 class(meshfield3d), pointer :: ddens, dens_hyd
228 !--------------------------------------------
229
230 call prog_vars%Get3D(dens_vid, ddens)
231 call aux_vars%Get3D(denshyd_vid, dens_hyd)
232
233 call prog_vars%MeshFieldComm_Exchange()
234
235 call apply_numfilter( this, therm_vid, prog_vars, boundary_cond, ddens, dens_hyd, dx, dy, dz, lift, mesh )
236 call apply_numfilter( this, momz_vid, prog_vars, boundary_cond, ddens, dens_hyd, dx, dy, dz, lift, mesh )
237 call apply_numfilter( this, momx_vid, prog_vars, boundary_cond, ddens, dens_hyd, dx, dy, dz, lift, mesh )
238 call apply_numfilter( this, momy_vid, prog_vars, boundary_cond, ddens, dens_hyd, dx, dy, dz, lift, mesh )
239 call apply_numfilter( this, dens_vid, prog_vars, boundary_cond, ddens, dens_hyd, dx, dy, dz, lift, mesh )
240
241 return
242 end subroutine atm_dyn_dgm_nonhydro3d_numdiff_apply
243
244 !- private ------------------------------
245
246!OCL SERIAL
247 subroutine apply_numfilter( this, varid, prgvars_list, boundary_cond, DDENS, DENS_hyd, &
248 Dx, Dy, Dz, Lift, mesh )
249
250 use scale_atm_dyn_dgm_nonhydro3d_common, only: &
251 prgvar_ddens_id
252
253 implicit none
254
255 class(atmdyn_nonhydro3d_numdiff), intent(inout) :: this
256 integer, intent(in) :: varid
257 class(modelvarmanager), intent(inout) :: prgvars_list
258 class(atmdynbnd), intent(in) :: boundary_cond
259 class(meshfield3d), intent(in) :: ddens
260 class(meshfield3d), intent(in) :: dens_hyd
261 type(sparsemat), intent(in) :: dx, dy, dz, lift
262 class(meshbase3d), intent(in), target :: mesh
263
264 class(meshfield3d), pointer :: var
265 class(meshfield3d), pointer :: nd_flx_x, nd_flx_y, nd_flx_z
266 class(meshfield3d), pointer :: nd_lapla_h, nd_lapla_v
267
268 class(localmesh3d), pointer :: lcmesh
269 integer :: n
270 integer :: ke
271
272 integer :: nd_itr
273 real(rp) :: nd_sign
274 logical :: dens_weight_flag
275 logical, allocatable :: is_bound(:,:)
276
277 real(rp), allocatable :: tmp_tend(:,:)
278 !-----------------------------------------
279
280 nd_sign = (-1)**(mod(this%ND_LAPLACIAN_NUM+1,2))
281 dens_weight_flag = (varid /= prgvar_ddens_id)
282
283 call prgvars_list%Get3D(varid, var)
284 call this%NUMDIFF_FLUX_manager%Get3D(numdiffflx_x_id, nd_flx_x)
285 call this%NUMDIFF_FLUX_manager%Get3D(numdiffflx_y_id, nd_flx_y)
286 call this%NUMDIFF_FLUX_manager%Get3D(numdiffflx_z_id, nd_flx_z)
287
288 call this%NUMDIFF_TEND_manager%Get3D(numdiff_laplah_id, nd_lapla_h)
289 call this%NUMDIFF_TEND_manager%Get3D(numdiff_laplav_id, nd_lapla_v)
290
291 do n=1, mesh%LOCAL_MESH_NUM
292 lcmesh => mesh%lcmesh_list(n)
293
294 allocate( is_bound(lcmesh%refElem%NfpTot,lcmesh%Ne) )
295 call boundary_cond%ApplyBC_numdiff_even_lc( var%local(n)%val, is_bound, varid, n, &
296 lcmesh%normal_fn(:,:,1), lcmesh%normal_fn(:,:,2), lcmesh%normal_fn(:,:,3), &
297 lcmesh%vmapM, lcmesh%vmapP, lcmesh%vmapB, lcmesh, lcmesh%refElem3D )
298
299 call numdiff_cal_flx( &
300 nd_flx_x%local(n)%val, nd_flx_y%local(n)%val, nd_flx_z%local(n)%val, &
301 var%local(n)%val, var%local(n)%val, ddens%local(n)%val, dens_hyd%local(n)%val, &
302 dx, dy, dz, lift, lcmesh, lcmesh%refElem3D, is_bound, dens_weight_flag )
303
304 deallocate( is_bound )
305 end do
306
307 !* Exchange halo data
308 call this%NUMDIFF_FLUX_manager%MeshFieldComm_Exchange()
309
310 do nd_itr=1, this%ND_LAPLACIAN_NUM-1
311 do n = 1, mesh%LOCAL_MESH_NUM
312 lcmesh => mesh%lcmesh_list(n)
313
314 allocate( is_bound(lcmesh%refElem%NfpTot,lcmesh%Ne) )
315 call boundary_cond%ApplyBC_numdiff_odd_lc( &
316 nd_flx_x%local(n)%val, nd_flx_y%local(n)%val, nd_flx_z%local(n)%val, is_bound, varid, n, &
317 lcmesh%normal_fn(:,:,1), lcmesh%normal_fn(:,:,2), lcmesh%normal_fn(:,:,3), &
318 lcmesh%vmapM, lcmesh%vmapP, lcmesh%vmapB, lcmesh, lcmesh%refElem3D )
319
320 call numdiff_cal_laplacian( &
321 nd_lapla_h%local(n)%val, nd_lapla_v%local(n)%val, &
322 nd_flx_x%local(n)%val, nd_flx_y%local(n)%val, nd_flx_z%local(n)%val, &
323 dx, dy, dz, lift, lcmesh, lcmesh%refElem3D, is_bound )
324
325 deallocate( is_bound )
326 end do
327 !* Exchange halo data
328 call this%NUMDIFF_TEND_manager%MeshFieldComm_Exchange()
329
330 do n = 1, mesh%LOCAL_MESH_NUM
331 lcmesh => mesh%lcmesh_list(n)
332
333 allocate( is_bound(lcmesh%refElem%NfpTot,lcmesh%Ne) )
334 call boundary_cond%ApplyBC_numdiff_even_lc( &
335 nd_lapla_h%local(n)%val, is_bound, varid, n, &
336 lcmesh%normal_fn(:,:,1), lcmesh%normal_fn(:,:,2), lcmesh%normal_fn(:,:,3), &
337 lcmesh%vmapM, lcmesh%vmapP, lcmesh%vmapB, lcmesh, lcmesh%refElem3D )
338
339 call numdiff_cal_flx( &
340 nd_flx_x%local(n)%val, nd_flx_y%local(n)%val, nd_flx_z%local(n)%val, &
341 nd_lapla_h%local(n)%val, nd_lapla_v%local(n)%val, &
342 ddens%local(n)%val, dens_hyd%local(n)%val, &
343 dx, dy, dz, lift, lcmesh, lcmesh%refElem3D, is_bound, .false. )
344
345 deallocate( is_bound )
346 end do
347 !* Exchange halo data
348 call this%NUMDIFF_FLUX_manager%MeshFieldComm_Exchange()
349 end do
350
351 do n = 1, mesh%LOCAL_MESH_NUM
352
353 allocate( is_bound(lcmesh%refElem%NfpTot,lcmesh%Ne), tmp_tend(lcmesh%refElem3D%Np,lcmesh%NeA) )
354
355 call boundary_cond%ApplyBC_numdiff_odd_lc( &
356 nd_flx_x%local(n)%val, nd_flx_y%local(n)%val, nd_flx_z%local(n)%val, &
357 is_bound, varid, n, &
358 lcmesh%normal_fn(:,:,1), lcmesh%normal_fn(:,:,2), lcmesh%normal_fn(:,:,3), &
359 lcmesh%vmapM, lcmesh%vmapP, lcmesh%vmapB, lcmesh, lcmesh%refElem3D )
360
361 call numdiff_tend( tmp_tend(:,:), &
362 nd_flx_x%local(n)%val, nd_flx_y%local(n)%val, nd_flx_z%local(n)%val, &
363 ddens%local(n)%val, dens_hyd%local(n)%val, &
364 nd_sign * this%ND_COEF_H, nd_sign * this%ND_COEF_V, &
365 dx, dy, dz, lift, lcmesh, lcmesh%refElem3D, is_bound, dens_weight_flag )
366
367 !$omp parallel do
368 do ke=lcmesh%NeS, lcmesh%NeE
369 var%local(n)%val(:,ke) = var%local(n)%val(:,ke) + this%dtsec * tmp_tend(:,ke)
370 end do
371
372 deallocate( is_bound, tmp_tend )
373 end do
374
375 return
376 end subroutine apply_numfilter
377
378!OCL SERIAL
379 subroutine numdiff_tend( &
380 tend_, & ! (out)
381 gxv_, gyv_, gzv_, & ! (in)
382 ddens_, dens_hyd, diffcoef_h, diffcoef_v, & ! (in)
383 dx, dy, dz, lift, lmesh, elem, is_bound, mul_dens_flag )
384
385 implicit none
386
387 class(localmesh3d), intent(in) :: lmesh
388 class(elementbase3d), intent(in) :: elem
389 type(sparsemat), intent(in) :: dx, dy, dz, lift
390 real(rp), intent(inout) :: tend_(elem%np,lmesh%nea)
391 real(rp), intent(in) :: gxv_(elem%np,lmesh%nea)
392 real(rp), intent(in) :: gyv_(elem%np,lmesh%nea)
393 real(rp), intent(in) :: gzv_(elem%np,lmesh%nea)
394 real(rp), intent(in) :: ddens_(elem%np,lmesh%nea)
395 real(rp), intent(in) :: dens_hyd(elem%np,lmesh%nea)
396 real(rp), intent(in) :: diffcoef_h
397 real(rp), intent(in) :: diffcoef_v
398 logical, intent(in) :: is_bound(elem%nfptot,lmesh%ne)
399 logical, intent(in) :: mul_dens_flag
400
401 real(rp) :: del_flux(elem%nfptot,lmesh%ne)
402 real(rp) :: coef_h(elem%np), coef_v(elem%np)
403 real(rp) :: fx(elem%np), fy(elem%np), fz(elem%np), liftdelflx(elem%np)
404
405 integer :: ke
406
407 !--------------------------------------------
408
409 call cal_del_flux_lap_with_coef( del_flux, & ! (out)
410 gxv_, gyv_, gzv_, & ! (in)
411 ddens_, dens_hyd, diffcoef_h, diffcoef_v, & ! (in)
412 lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), & ! (in)
413 lmesh%vmapM, lmesh%vmapP, & ! (in)
414 lmesh, elem, is_bound, mul_dens_flag ) ! (in)
415
416 !$omp parallel do private( &
417 !$omp coef_h, coef_v, Fx, Fy, Fz, LiftDelFlx )
418 do ke = lmesh%NeS, lmesh%NeE
419
420 if (mul_dens_flag) then
421 coef_h(:) = diffcoef_h * (dens_hyd(:,ke) + ddens_(:,ke))
422 coef_v(:) = diffcoef_v * (dens_hyd(:,ke) + ddens_(:,ke))
423 else
424 coef_h(:) = diffcoef_h
425 coef_v(:) = diffcoef_v
426 end if
427 call sparsemat_matmul(dx, coef_h(:) * gxv_(:,ke), fx)
428 call sparsemat_matmul(dy, coef_h(:) * gyv_(:,ke), fy)
429 call sparsemat_matmul(dz, coef_v(:) * gzv_(:,ke) ,fz)
430 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke), liftdelflx)
431
432 tend_(:,ke) = ( &
433 lmesh%Escale(:,ke,1,1) * fx(:) &
434 + lmesh%Escale(:,ke,2,2) * fy(:) &
435 + lmesh%Escale(:,ke,3,3) * fz(:) &
436 + liftdelflx(:) )
437 end do
438
439 return
440 end subroutine numdiff_tend
441
442!OCL SERIAL
443 subroutine cal_del_flux_lap_with_coef( del_flux, & ! (out)
444 gxv_, gyv_, gzv_, & ! (in)
445 ddens_, dens_hyd, coef_h, coef_v, & ! (in)
446 nx, ny, nz, vmapm, vmapp, lmesh, elem, & ! (in)
447 is_bound, mul_dens_flag ) ! (in)
448
449 implicit none
450
451 class(localmesh3d), intent(in) :: lmesh
452 class(elementbase3d), intent(in) :: elem
453 real(rp), intent(out) :: del_flux(elem%nfptot*lmesh%ne)
454 real(rp), intent(in) :: gxv_(elem%np*lmesh%nea)
455 real(rp), intent(in) :: gyv_(elem%np*lmesh%nea)
456 real(rp), intent(in) :: gzv_(elem%np*lmesh%nea)
457 real(rp), intent(in) :: ddens_(elem%np*lmesh%nea)
458 real(rp), intent(in) :: dens_hyd(elem%np*lmesh%nea)
459 real(rp), intent(in) :: coef_h
460 real(rp), intent(in) :: coef_v
461 real(rp), intent(in) :: nx(elem%nfptot*lmesh%ne)
462 real(rp), intent(in) :: ny(elem%nfptot*lmesh%ne)
463 real(rp), intent(in) :: nz(elem%nfptot*lmesh%ne)
464 integer, intent(in) :: vmapm(elem%nfptot*lmesh%ne)
465 integer, intent(in) :: vmapp(elem%nfptot*lmesh%ne)
466 logical, intent(in) :: is_bound(elem%nfptot*lmesh%ne)
467 logical, intent(in) :: mul_dens_flag
468
469 integer :: i, ip, im
470 real(rp) :: weightp, weightm
471 !------------------------------------------------------------------------
472
473 !$omp parallel do &
474 !$omp private( iM, iP, weightP, weightM )
475 do i=1, elem%NfpTot*lmesh%Ne
476 im = vmapm(i); ip = vmapp(i)
477
478 if (mul_dens_flag) then
479 weightm = 0.5_rp * (dens_hyd(im) + ddens_(im))
480 weightp = 0.5_rp * (dens_hyd(ip) + ddens_(ip))
481 else
482 weightm = 0.5_rp
483 weightp = 0.5_rp
484 end if
485
486 if (is_bound(i)) then
487 del_flux(i) = &
488 coef_h * ( weightp * gxv_(ip) - weightm * gxv_(im) ) * nx(i) &
489 + coef_h * ( weightp * gyv_(ip) - weightm * gyv_(im) ) * ny(i) &
490 + coef_v * ( weightp * gzv_(ip) - weightm * gzv_(im) ) * nz(i)
491 else
492 del_flux(i) = &
493 ( 1.0_rp + sign(1.0_rp,nx(i)) ) * coef_h * ( weightp * gxv_(ip) - weightm * gxv_(im) ) * nx(i) &
494 + ( 1.0_rp + sign(1.0_rp,ny(i)) ) * coef_h * ( weightp * gyv_(ip) - weightm * gyv_(im) ) * ny(i) &
495 + ( 1.0_rp + sign(1.0_rp,nz(i)) ) * coef_v * ( weightp * gzv_(ip) - weightm * gzv_(im) ) * nz(i)
496 end if
497 end do
498
499 return
500 end subroutine cal_del_flux_lap_with_coef
501
502 !--
503
504!OCL SERIAL
505 subroutine numdiff_cal_laplacian( &
506 lapla_h, lapla_v, & ! (out)
507 gxv_, gyv_, gzv_, & ! (in)
508 dx, dy, dz, lift, lmesh, elem, is_bound ) ! (in)
509
510 implicit none
511
512 class(localmesh3d), intent(in) :: lmesh
513 class(elementbase3d), intent(in) :: elem
514 type(sparsemat), intent(in) :: dx, dy, dz, lift
515 real(rp), intent(out) :: lapla_h(elem%np,lmesh%nea)
516 real(rp), intent(out) :: lapla_v(elem%np,lmesh%nea)
517 real(rp), intent(in) :: gxv_(elem%np,lmesh%nea)
518 real(rp), intent(in) :: gyv_(elem%np,lmesh%nea)
519 real(rp), intent(in) :: gzv_(elem%np,lmesh%nea)
520 logical, intent(in) :: is_bound(elem%nfptot,lmesh%ne)
521
522 real(rp) :: del_flux_h(elem%nfptot,lmesh%ne)
523 real(rp) :: del_flux_v(elem%nfptot,lmesh%ne)
524 real(rp) :: fx(elem%np), fy(elem%np), fz(elem%np), liftdelflx(elem%np)
525
526 integer :: ke
527 !--------------------------------------------
528
529 call cal_del_flux_lap( del_flux_h, del_flux_v, & ! (out)
530 gxv_, gyv_, gzv_, & ! (in)
531 lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), & ! (in)
532 lmesh%vmapM, lmesh%vmapP, & ! (in)
533 lmesh, elem, is_bound ) ! (in)
534
535 !$omp parallel do private( &
536 !$omp Fx, Fy, Fz, LiftDelFlx )
537 do ke = lmesh%NeS, lmesh%NeE
538
539 call sparsemat_matmul(dx, gxv_(:,ke), fx)
540 call sparsemat_matmul(dy, gyv_(:,ke), fy)
541 call sparsemat_matmul(lift, lmesh%Fscale(:,ke)*del_flux_h(:,ke), liftdelflx)
542
543 lapla_h(:,ke) = ( &
544 lmesh%Escale(:,ke,1,1) * fx(:) &
545 + lmesh%Escale(:,ke,2,2) * fy(:) &
546 + liftdelflx(:) &
547 )
548
549 call sparsemat_matmul(dz, gzv_(:,ke) ,fz)
550 call sparsemat_matmul(lift, lmesh%Fscale(:,ke)*del_flux_v(:,ke), liftdelflx)
551 lapla_v(:,ke) = ( &
552 lmesh%Escale(:,ke,3,3) * fz(:) &
553 + liftdelflx(:) &
554 )
555 end do
556
557 return
558 end subroutine numdiff_cal_laplacian
559
560!OCL SERIAL
561 subroutine cal_del_flux_lap( del_flux_h, del_flux_v, & ! (out)
562 gxv_, gyv_, gzv_, & ! (in)
563 nx, ny, nz, vmapm, vmapp, lmesh, elem, is_bound ) ! (in)
564
565 implicit none
566
567 class(localmesh3d), intent(in) :: lmesh
568 class(elementbase3d), intent(in) :: elem
569 real(rp), intent(out) :: del_flux_h(elem%nfptot*lmesh%ne)
570 real(rp), intent(out) :: del_flux_v(elem%nfptot*lmesh%ne)
571 real(rp), intent(in) :: gxv_(elem%np*lmesh%nea)
572 real(rp), intent(in) :: gyv_(elem%np*lmesh%nea)
573 real(rp), intent(in) :: gzv_(elem%np*lmesh%nea)
574 real(rp), intent(in) :: nx(elem%nfptot*lmesh%ne)
575 real(rp), intent(in) :: ny(elem%nfptot*lmesh%ne)
576 real(rp), intent(in) :: nz(elem%nfptot*lmesh%ne)
577 integer, intent(in) :: vmapm(elem%nfptot*lmesh%ne)
578 integer, intent(in) :: vmapp(elem%nfptot*lmesh%ne)
579 logical, intent(in) :: is_bound(elem%nfptot*lmesh%ne)
580
581 integer :: i, ip, im
582 !------------------------------------------------------------------------
583
584 !$omp parallel do &
585 !$omp private( iM, iP )
586 do i=1, elem%NfpTot*lmesh%Ne
587 im = vmapm(i); ip = vmapp(i)
588
589 if (is_bound(i)) then
590 del_flux_h(i) = 0.5_rp * ( ( gxv_(ip) - gxv_(im) ) * nx(i) &
591 + ( gyv_(ip) - gyv_(im) ) * ny(i) )
592 del_flux_v(i) = 0.5_rp * ( gzv_(ip) - gzv_(im) ) * nz(i)
593 else
594 del_flux_h(i) = 0.5_rp * ( ( 1.0_rp + sign(1.0_rp,nx(i)) ) * ( gxv_(ip) - gxv_(im) ) * nx(i) &
595 + ( 1.0_rp + sign(1.0_rp,ny(i)) ) * ( gyv_(ip) - gyv_(im) ) * ny(i) )
596 del_flux_v(i) = 0.5_rp * ( 1.0_rp + sign(1.0_rp,nz(i)) ) * ( gzv_(ip) - gzv_(im) ) * nz(i)
597 end if
598 end do
599
600 return
601 end subroutine cal_del_flux_lap
602
603 !-------------------------------------------------------
604
605!OCL SERIAL
606 subroutine numdiff_cal_flx( &
607 GxV_, GyV_, GzV_, & ! (out)
608 varh_, varv_, ddens_, dens_hyd_, & ! (in)
609 dx, dy, dz, lift, lmesh, elem, is_bound, divide_dens_flag ) ! (in)
610
611 implicit none
612
613 class(localmesh3d), intent(in) :: lmesh
614 class(elementbase3d), intent(in) :: elem
615 type(sparsemat), intent(in) :: dx, dy, dz, lift
616 real(rp), intent(out) :: gxv_(elem%np,lmesh%nea)
617 real(rp), intent(out) :: gyv_(elem%np,lmesh%nea)
618 real(rp), intent(out) :: gzv_(elem%np,lmesh%nea)
619 real(rp), intent(in) :: varh_(elem%np,lmesh%nea)
620 real(rp), intent(in) :: varv_(elem%np,lmesh%nea)
621 real(rp), intent(in) :: ddens_(elem%np,lmesh%nea)
622 real(rp), intent(in) :: dens_hyd_(elem%np,lmesh%nea)
623 logical, intent(in) :: is_bound(elem%nfptot,lmesh%ne)
624 logical, intent(in) :: divide_dens_flag
625
626 real(rp) :: fx(elem%np), fy(elem%np), fz(elem%np), liftdelflx(elem%np)
627 real(rp) :: vh(elem%np), vv(elem%np)
628 real(rp) :: del_flux(elem%nfptot,lmesh%ne,3)
629
630 integer :: ke
631
632 !------------------------------------------------------------------------------
633
634 call cal_del_graddiffvar( del_flux, & ! (out)
635 varh_, varv_, ddens_, dens_hyd_, & ! (in)
636 lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), & ! (in)
637 lmesh%vmapM, lmesh%vmapP, & ! (in)
638 lmesh, elem, is_bound, divide_dens_flag ) ! (in)
639
640 !$omp parallel do private(Fx, Fy, Fz, LiftDelFlx, vh, vv)
641 do ke=lmesh%NeS, lmesh%NeE
642
643 if (divide_dens_flag) then
644 vh(:) = varh_(:,ke) / (ddens_(:,ke) + dens_hyd_(:,ke))
645 vv(:) = varv_(:,ke) / (ddens_(:,ke) + dens_hyd_(:,ke))
646 else
647 vh(:) = varh_(:,ke)
648 vv(:) = varv_(:,ke)
649 end if
650
651 call sparsemat_matmul(dx, vh, fx)
652 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,1), liftdelflx)
653 gxv_(:,ke) = lmesh%Escale(:,ke,1,1) * fx(:) + liftdelflx(:)
654
655 call sparsemat_matmul(dy, vh, fy)
656 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,2), liftdelflx)
657 gyv_(:,ke) = lmesh%Escale(:,ke,2,2) * fy(:) + liftdelflx(:)
658
659 call sparsemat_matmul(dz, vv, fz)
660 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,3), liftdelflx)
661 gzv_(:,ke) = lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:)
662 end do
663
664 return
665 end subroutine numdiff_cal_flx
666
667!OCL SERIAL
668 subroutine cal_del_graddiffvar( del_flux, & ! (out)
669 varh_, varv_, ddens_, dens_hyd_, & ! (in)
670 nx, ny, nz, vmapm, vmapp, lmesh, elem, & ! (in)
671 is_bound, divide_dens_flag ) ! (in)
672
673 implicit none
674
675 class(localmesh3d), intent(in) :: lmesh
676 class(elementbase3d), intent(in) :: elem
677 real(rp), intent(out) :: del_flux(elem%nfptot*lmesh%ne,3)
678 real(rp), intent(in) :: varh_(elem%np*lmesh%nea)
679 real(rp), intent(in) :: varv_(elem%np*lmesh%nea)
680 real(rp), intent(in) :: ddens_(elem%np*lmesh%nea)
681 real(rp), intent(in) :: dens_hyd_(elem%np*lmesh%nea)
682 real(rp), intent(in) :: nx(elem%nfptot*lmesh%ne)
683 real(rp), intent(in) :: ny(elem%nfptot*lmesh%ne)
684 real(rp), intent(in) :: nz(elem%nfptot*lmesh%ne)
685 integer, intent(in) :: vmapm(elem%nfptot*lmesh%ne)
686 integer, intent(in) :: vmapp(elem%nfptot*lmesh%ne)
687 logical, intent(in) :: is_bound(elem%nfptot*lmesh%ne)
688 logical, intent(in) :: divide_dens_flag
689
690 integer :: i, ip, im
691 real(rp) :: delvarh, delvarv
692 real(rp) :: weight_p, weight_m
693 !------------------------------------------------------------------------
694
695 !$omp parallel do private(iM, iP, delVarh, delVarv, weight_P, weight_M)
696 do i=1, elem%NfpTot*lmesh%Ne
697 im = vmapm(i); ip = vmapp(i)
698
699 if (divide_dens_flag) then
700 weight_p = 1.0_rp / (ddens_(ip) + dens_hyd_(ip))
701 weight_m = 1.0_rp / (ddens_(im) + dens_hyd_(im))
702 else
703 weight_m = 1.0_rp
704 weight_p = 1.0_rp
705 end if
706
707 delvarh = 0.5_rp * (varh_(ip) * weight_p - varh_(im) * weight_m)
708 delvarv = 0.5_rp * (varv_(ip) * weight_p - varv_(im) * weight_m)
709 if (is_bound(i)) then
710 del_flux(i,1) = delvarh * nx(i)
711 del_flux(i,2) = delvarh * ny(i)
712 del_flux(i,3) = delvarv * nz(i)
713 else
714 del_flux(i,1) = ( 1.0_rp - sign(1.0_rp,nx(i)) ) * delvarh * nx(i)
715 del_flux(i,2) = ( 1.0_rp - sign(1.0_rp,ny(i)) ) * delvarh * ny(i)
716 del_flux(i,3) = ( 1.0_rp - sign(1.0_rp,nz(i)) ) * delvarv * nz(i)
717 end if
718 end do
719
720 return
721 end subroutine cal_del_graddiffvar
722
723 !------------------------------------------------
724end module scale_atm_dyn_dgm_nonhydro3d_numdiff
scale_sparsemat::sparsemat_matmul
Definition scale_sparsemat.F90:57
scale_atm_dyn_dgm_bnd
module FElib / Fluid dyn solver / Atmosphere / Boundary
Definition scale_atm_dyn_dgm_bnd.F90:11
scale_atm_dyn_dgm_nonhydro3d_common
module FElib / Fluid dyn solver / Atmosphere / Nonhydrostatic model / Common
Definition scale_atm_dyn_dgm_nonhydro3d_common.F90:11
scale_atm_dyn_dgm_nonhydro3d_common::prgvar_momy_id
integer, parameter, public prgvar_momy_id
Definition scale_atm_dyn_dgm_nonhydro3d_common.F90:74
scale_atm_dyn_dgm_nonhydro3d_common::prgvar_therm_id
integer, parameter, public prgvar_therm_id
Definition scale_atm_dyn_dgm_nonhydro3d_common.F90:69
scale_atm_dyn_dgm_nonhydro3d_common::prgvar_ddens_id
integer, parameter, public prgvar_ddens_id
Definition scale_atm_dyn_dgm_nonhydro3d_common.F90:68
scale_atm_dyn_dgm_nonhydro3d_common::prgvar_momz_id
integer, parameter, public prgvar_momz_id
Definition scale_atm_dyn_dgm_nonhydro3d_common.F90:72
scale_atm_dyn_dgm_nonhydro3d_common::prgvar_momx_id
integer, parameter, public prgvar_momx_id
Definition scale_atm_dyn_dgm_nonhydro3d_common.F90:73
scale_atm_dyn_dgm_nonhydro3d_common::prgvar_num
integer, parameter, public prgvar_num
Definition scale_atm_dyn_dgm_nonhydro3d_common.F90:77
scale_atm_dyn_dgm_nonhydro3d_common::auxvar_denshydro_id
integer, parameter, public auxvar_denshydro_id
Definition scale_atm_dyn_dgm_nonhydro3d_common.F90:89
scale_atm_dyn_dgm_nonhydro3d_numdiff
module Atmosphere / Dynamics common
Definition scale_atm_dyn_dgm_nonhydro3d_numdiff.F90:12
scale_atm_dyn_dgm_nonhydro3d_numdiff::atmos_dyn_numdiff_tend_vinfo
type(variableinfo), dimension(numdiff_tend_num), public atmos_dyn_numdiff_tend_vinfo
Definition scale_atm_dyn_dgm_nonhydro3d_numdiff.F90:96
scale_atm_dyn_dgm_nonhydro3d_numdiff::numdiff_laplah_id
integer, parameter, public numdiff_laplah_id
Definition scale_atm_dyn_dgm_nonhydro3d_numdiff.F90:93
scale_atm_dyn_dgm_nonhydro3d_numdiff::numdiff_laplav_id
integer, parameter, public numdiff_laplav_id
Definition scale_atm_dyn_dgm_nonhydro3d_numdiff.F90:94
scale_atm_dyn_dgm_nonhydro3d_numdiff::numdiffflx_x_id
integer, parameter, public numdiffflx_x_id
Definition scale_atm_dyn_dgm_nonhydro3d_numdiff.F90:78
scale_atm_dyn_dgm_nonhydro3d_numdiff::numdiffflx_y_id
integer, parameter, public numdiffflx_y_id
Definition scale_atm_dyn_dgm_nonhydro3d_numdiff.F90:79
scale_atm_dyn_dgm_nonhydro3d_numdiff::numdiff_flux_num
integer, parameter, public numdiff_flux_num
Definition scale_atm_dyn_dgm_nonhydro3d_numdiff.F90:77
scale_atm_dyn_dgm_nonhydro3d_numdiff::atmos_dyn_numdiff_flux_vinfo
type(variableinfo), dimension(numdiff_flux_num), public atmos_dyn_numdiff_flux_vinfo
Definition scale_atm_dyn_dgm_nonhydro3d_numdiff.F90:82
scale_atm_dyn_dgm_nonhydro3d_numdiff::numdiffflx_z_id
integer, parameter, public numdiffflx_z_id
Definition scale_atm_dyn_dgm_nonhydro3d_numdiff.F90:80
scale_atm_dyn_dgm_nonhydro3d_numdiff::numdiff_tend_num
integer, parameter, public numdiff_tend_num
Definition scale_atm_dyn_dgm_nonhydro3d_numdiff.F90:92
scale_element_base
module FElib / Element / Base
Definition scale_element_base.F90:11
scale_element_hexahedral
module FElib / Element / hexahedron
Definition scale_element_hexahedral.F90:10
scale_localmesh_2d
module FElib / Mesh / Local 2D
Definition scale_localmesh_2d.F90:10
scale_localmesh_3d
module FElib / Mesh / Local 3D
Definition scale_localmesh_3d.F90:10
scale_localmeshfield_base
module FElib / Data / base
Definition scale_localmeshfield_base.F90:10
scale_mesh_base3d
module FElib / Mesh / Base 3D
Definition scale_mesh_base3d.F90:10
scale_mesh_base
module FElib / Mesh / Base
Definition scale_mesh_base.F90:10
scale_meshfield_base
module FElib / Data / base
Definition scale_meshfield_base.F90:10
scale_model_mesh_manager
FElib / model framework / mesh manager.
Definition scale_model_mesh_manager.F90:12
scale_model_var_manager
FElib / model framework / variable manager.
Definition scale_model_var_manager.F90:12
scale_sparsemat
Module common / sparsemat.
Definition scale_sparsemat.F90:11
scale_atm_dyn_dgm_bnd::atmdynbnd
A derived type useful for apply boundary conditions.
Definition scale_atm_dyn_dgm_bnd.F90:57
scale_atm_dyn_dgm_nonhydro3d_numdiff::atmdyn_nonhydro3d_numdiff
Definition scale_atm_dyn_dgm_nonhydro3d_numdiff.F90:53
scale_element_base::elementbase2d
Derived type representing a 2D reference element.
Definition scale_element_base.F90:71
scale_element_base::elementbase3d
Derived type representing a 3D reference element.
Definition scale_element_base.F90:111
scale_element_hexahedral::hexahedralelement
Derived type representing a hexahedral element.
Definition scale_element_hexahedral.F90:31
scale_localmesh_2d::localmesh2d
Derived type representing a local mesh for 2D domain.
Definition scale_localmesh_2d.F90:30
scale_localmesh_3d::localmesh3d
Derived type to manage a local 3D computational domain.
Definition scale_localmesh_3d.F90:32
scale_localmeshfield_base::localmeshfield3d
Derived type representing a field with 3D local mesh.
Definition scale_localmeshfield_base.F90:69
scale_mesh_base3d::meshbase3d
Derived type to manage a computational mesh (base type for 3D domain)
Definition scale_mesh_base3d.F90:37
scale_mesh_base::meshbase
Base type to manage a computational mesh.
Definition scale_mesh_base.F90:40
scale_meshfield_base::meshfield3d
Derived type representing a field with 3D mesh.
Definition scale_meshfield_base.F90:95
scale_model_mesh_manager::modelmesh3d
Definition scale_model_mesh_manager.F90:57
scale_model_var_manager::modelvarmanager
Definition scale_model_var_manager.F90:41
scale_sparsemat::sparsemat
Derived type to manage a sparse matrix.
Definition scale_sparsemat.F90:33
Generated by doxygen 1.14.0