FE-Project
Loading...
Searching...
No Matches
scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common.F90
Go to the documentation of this file.
1!-------------------------------------------------------------------------------
2!> module FElib / Fluid dyn solver / Atmosphere / Nonhydrostatic model / HEVI / Common
3!!
4!! @par Description
5!! HEVI DGM scheme for Atmospheric dynamical process
6!! As the thermodynamics equation, a density * potential temperature equation is used.
7!!
8!! @author Yuta Kawai, Team SCALE
9!<
10!-------------------------------------------------------------------------------
11#include "scaleFElib.h"
12module scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common
13 !-----------------------------------------------------------------------------
14 !
15 !++ Used modules
16 !
17 use scale_precision
18 use scale_io
19 use scale_prc
20 use scale_prof
21 use scale_const, only: &
22 grav => const_grav, &
23 rdry => const_rdry, &
24 cpdry => const_cpdry, &
25 cvdry => const_cvdry, &
26 pres00 => const_pre00
27
28 use scale_sparsemat
29 use scale_element_base, only: &
30 elementbase2d, elementbase3d
31 use scale_element_modalfilter, only: modalfilter
32 use scale_element_hexahedral, only: hexahedralelement
33 use scale_localmesh_2d, only: localmesh2d
34 use scale_localmesh_3d, only: localmesh3d
35 use scale_mesh_base3d, only: meshbase3d
36 use scale_localmeshfield_base, only: localmeshfield3d
37 use scale_meshfield_base, only: meshfield3d
38
39 use scale_atm_dyn_dgm_nonhydro3d_common, only: &
40 dens_vid => prgvar_ddens_id, rhot_vid => prgvar_drhot_id, &
41 momx_vid => prgvar_momx_id, momy_vid => prgvar_momy_id, &
42 momz_vid => prgvar_momz_id, &
43 prgvar_num
44
45 !-----------------------------------------------------------------------------
46 implicit none
47 private
48 !-----------------------------------------------------------------------------
49 !
50 !++ Public procedures
51 !
52 public :: atm_dyn_dgm_nonhydro3d_rhot_hevi_common_init
53 public :: atm_dyn_dgm_nonhydro3d_rhot_hevi_common_final
54
55 public :: atm_dyn_dgm_nonhydro3d_rhot_hevi_common_eval_ax
56 public :: atm_dyn_dgm_nonhydro3d_rhot_hevi_common_eval_ax_uv
57 public :: atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd
58 public :: atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd_uv
59
60 !-----------------------------------------------------------------------------
61 !
62 !++ Public parameters & variables
63 !
64 logical, public :: vi_use_lapack_flag = .true.
65
66 !-----------------------------------------------------------------------------
67 !
68 !++ Private procedures & variables
69 !
70 !-------------------
71
72 private :: vi_cal_del_flux_dyn
73 private :: vi_cal_del_flux_dyn_uv
74
75contains
76 !------------------------------------------------
77 subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_init
78 implicit none
79 namelist / param_atmos_dyn_nonhydro3d_rhot_hevi_common / &
80 vi_use_lapack_flag
81
82 integer :: ierr
83 !---------------------------------------------------------
84
85 !--- read namelist
86 rewind(io_fid_conf)
87 read(io_fid_conf,nml=param_atmos_dyn_nonhydro3d_rhot_hevi_common,iostat=ierr)
88 if( ierr < 0 ) then !--- missing
89 log_info("ATMOS_DYN_nonhydro3d_rhot_hevi_common_setup",*) 'Not found namelist. Default used.'
90 elseif( ierr > 0 ) then !--- fatal error
91 log_error("ATMOS_DYN_nonhydro3d_rhot_hevi_common_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_DYN_NONHYDRO3D_RHOT_HEVI_COMMON. Check!'
92 call prc_abort
93 endif
94 log_nml(param_atmos_dyn_nonhydro3d_rhot_hevi_common)
95
96 return
97 end subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_init
98
99 subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_final
100 implicit none
101 !---------------------------------------------------------
102 return
103 end subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_final
104
105!OCL SERIAL
106 subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_eval_ax( &
107 DENS_t, MOMZ_t, RHOT_t, & ! (out)
108 alph, & ! (in)
109 prog_vars, prog_vars0, & ! (in)
110 ddens00, momx00, momy00, momz00, drhot00, & ! (in)
111 dens_hyd, pres_hyd, & ! (in)
112 rtot, cptot_ov_cvtot, & ! (in)
113 dz, lift, intrpmat_vpordm1, & ! (in)
114 gnnm, g13, g23, gsqrtv, & ! (in)
115 impl_fac, dt, & ! (in)
116 lmesh, elem, & ! (in)
117 nz, vmapm, vmapp, & ! (in)
118 b1d_ij ) ! (out)
119
120 implicit none
121
122 class(localmesh3d), intent(in) :: lmesh
123 class(elementbase3d), intent(in) :: elem
124 real(rp), intent(out) :: dens_t(elem%np,lmesh%nea)
125 real(rp), intent(out) :: momz_t(elem%np,lmesh%nea)
126 real(rp), intent(out) :: rhot_t(elem%np,lmesh%nea)
127 real(rp), intent(in) :: alph(elem%nfptot,lmesh%nez,lmesh%nex*lmesh%ney)
128 real(rp), intent(in) :: prog_vars (elem%np,lmesh%nez,prgvar_num,lmesh%nex*lmesh%ney)
129 real(rp), intent(in) :: prog_vars0 (elem%np,lmesh%nez,prgvar_num,lmesh%nex*lmesh%ney)
130 real(rp), intent(in) :: ddens00(elem%np,lmesh%nea)
131 real(rp), intent(in) :: momx00 (elem%np,lmesh%nea)
132 real(rp), intent(in) :: momy00 (elem%np,lmesh%nea)
133 real(rp), intent(in) :: momz00 (elem%np,lmesh%nea)
134 real(rp), intent(in) :: drhot00(elem%np,lmesh%nea)
135 real(rp), intent(in) :: dens_hyd(elem%np,lmesh%nez,lmesh%nex*lmesh%ney)
136 real(rp), intent(in) :: pres_hyd(elem%np,lmesh%nez,lmesh%nex*lmesh%ney)
137 real(rp), intent(in) :: rtot(elem%np,lmesh%nez,lmesh%nex*lmesh%ney)
138 real(rp), intent(in) :: cptot_ov_cvtot(elem%np,lmesh%nez,lmesh%nex*lmesh%ney)
139 class(sparsemat), intent(in) :: dz, lift
140 real(rp), intent(in) :: intrpmat_vpordm1(elem%np,elem%np)
141 real(rp), intent(in) :: gnnm(elem%np,lmesh%nez,lmesh%nex*lmesh%ney)
142 real(rp), intent(in) :: g13 (elem%np,lmesh%nez,lmesh%nex*lmesh%ney)
143 real(rp), intent(in) :: g23 (elem%np,lmesh%nez,lmesh%nex*lmesh%ney)
144 real(rp), intent(in) :: gsqrtv(elem%np,lmesh%nez,lmesh%nex*lmesh%ney)
145 real(rp), intent(in) :: impl_fac
146 real(rp), intent(in) :: dt
147 real(rp), intent(in) :: nz(elem%nfptot,lmesh%nez,lmesh%nex*lmesh%ney)
148 integer, intent(in) :: vmapm(elem%nfptot,lmesh%nez)
149 integer, intent(in) :: vmapp(elem%nfptot,lmesh%nez)
150 real(rp), intent(out), optional :: b1d_ij(3,elem%nnode_v,lmesh%nez,elem%nnode_h1d**2,lmesh%nex*lmesh%ney)
151
152 real(rp) :: rgsqrtv(elem%np)
153 real(rp) :: fscale(elem%nfptot), escale33(elem%np)
154 real(rp) :: fz(elem%np), liftdelflx(elem%np)
155 real(rp) :: del_flux(elem%nfptot,lmesh%nez,lmesh%nex*lmesh%ney,prgvar_num)
156 real(rp) :: momz(elem%np), ddens(elem%np), dpres(elem%np), rhot(elem%np)
157 real(rp) :: momw(elem%np)
158 integer :: ke_xy, ke_z
159 integer :: ke, ke2d
160 integer :: v
161 integer :: ij
162 integer :: colmask(elem%nnode_v)
163 real(rp) :: rdt
164 real(rp) :: drho(elem%np)
165
166 integer :: kk, kkk, p, pp
167 real(rp) :: vmf_v
168
169 real(rp) :: gamm, rgamm, rp0
170 !--------------------------------------------------------
171
172 gamm = cpdry/cvdry
173 rgamm = cvdry/cpdry
174 rdt = 1.0_rp / dt
175 rp0 = 1.0_rp / pres00
176
177 call vi_cal_del_flux_dyn( del_flux, &
178 alph, & ! (in)
179 prog_vars, prog_vars0, & ! (in)
180 dens_hyd, pres_hyd, & ! (in)
181 rtot, cptot_ov_cvtot, & ! (in)
182 gnnm, g13, g23, gsqrtv, nz, vmapm, vmapp, & ! (in)
183 lmesh, elem ) ! (in)
184
185 !$omp parallel private( ke_xy, ke_z, ke, ke2d, ij, v, &
186 !$omp MOMZ, MOMW, DDENS, DPRES, RHOT, Fz, LiftDelFlx, &
187 !$omp RGsqrtV, ColMask, Fscale, Escale33, drho, &
188 !$omp kk, p, kkk, vmf_v, pp )
189
190 !$omp do collapse(2)
191 do ke_xy=1, lmesh%NeX*lmesh%NeY
192 do ke_z=1, lmesh%NeZ
193 ke = ke_xy + (ke_z-1)*lmesh%NeX*lmesh%NeY
194 ke2d = lmesh%EMap3Dto2D(ke)
195
196 ddens(:) = prog_vars(:,ke_z,dens_vid,ke_xy)
197 drho(:) = matmul(intrpmat_vpordm1, ddens(:))
198
199 momz(:) = prog_vars(:,ke_z,momz_vid,ke_xy)
200
201 rhot(:) = pres00/rdry * (pres_hyd(:,ke_z,ke_xy)/pres00)**rgamm &
202 + prog_vars(:,ke_z,rhot_vid,ke_xy)
203 dpres(:) = pres00 * ( rtot(:,ke_z,ke_xy) * rp0 * rhot(:) )**cptot_ov_cvtot(:,ke_z,ke_xy) &
204 - pres_hyd(:,ke_z,ke_xy)
205
206 rgsqrtv(:) = 1.0_rp / gsqrtv(:,ke_z,ke_xy)
207 fscale(:) = lmesh%Fscale(:,ke)
208 escale33(:) = lmesh%Escale(:,ke,3,3)
209
210 momw(:) = momz(:) &
211 + gsqrtv(:,ke_z,ke_xy) * g13(:,ke_z,ke_xy) * prog_vars(:,ke_z,momx_vid,ke_xy) &
212 + gsqrtv(:,ke_z,ke_xy) * g23(:,ke_z,ke_xy) * prog_vars(:,ke_z,momy_vid,ke_xy)
213
214 !- DENS
215 call sparsemat_matmul(dz, momw(:), fz)
216 call sparsemat_matmul(lift, fscale(:) * del_flux(:,ke_z,ke_xy,dens_vid), liftdelflx)
217 dens_t(:,ke) = - ( escale33(:) * fz(:) + liftdelflx(:) ) * rgsqrtv(:)
218
219 !-MOMZ
220! call sparsemat_matmul(Dz, MOMZ(:)**2 / ( DENS_hyd(:,ke_z,ke_xy) + DDENS(:) ) + DPRES(:), Fz) ! [<- MOMZ x MOMZ / DENS + DPRES ]
221 call sparsemat_matmul(dz, dpres(:), fz)
222 call sparsemat_matmul(lift, fscale(:) * del_flux(:,ke_z,ke_xy,momz_vid), liftdelflx)
223 momz_t(:,ke) = - ( escale33(:) * fz(:) + liftdelflx(:) ) * rgsqrtv(:) &
224 - grav * drho(:)
225
226 !-RHOT
227 call sparsemat_matmul(dz, rhot(:) * momw(:) / ( dens_hyd(:,ke_z,ke_xy) + ddens(:) ), fz)
228 call sparsemat_matmul(lift, fscale(:) * del_flux(:,ke_z,ke_xy,rhot_vid), liftdelflx)
229 rhot_t(:,ke) = - ( escale33(:) * fz(:) + liftdelflx(:) ) * rgsqrtv(:)
230
231 end do
232 end do
233 !$omp end do
234
235 if ( present( b1d_ij ) ) then
236 !$omp do collapse(2)
237 do ke_xy=1, lmesh%NeX*lmesh%NeY
238 do ke_z=1, lmesh%NeZ
239 ke = ke_xy + (ke_z-1)*lmesh%NeX*lmesh%NeY
240
241 do ij=1, elem%Nnode_h1D**2
242 colmask(:) = elem%Colmask(:,ij)
243 b1d_ij(1,:,ke_z,ij,ke_xy) = impl_fac * dens_t(colmask(:),ke) &
244 - prog_vars(colmask(:),ke_z,dens_vid,ke_xy) &
245 + ddens00(colmask(:),ke)
246 b1d_ij(2,:,ke_z,ij,ke_xy) = impl_fac * momz_t(colmask(:),ke) &
247 - prog_vars(colmask(:),ke_z,momz_vid,ke_xy) &
248 + momz00(colmask(:),ke)
249 b1d_ij(3,:,ke_z,ij,ke_xy) = impl_fac * rhot_t(colmask(:),ke) &
250 - prog_vars(colmask(:),ke_z,rhot_vid,ke_xy) &
251 + drhot00(colmask(:),ke)
252 end do
253 end do
254 end do
255 !$omp end do
256 end if
257 !$omp end parallel
258
259 return
260 end subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_eval_ax
261
262!OCL SERIAL
263 subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_eval_ax_uv( &
264 MOMX_t, MOMY_t, & ! (out)
265 alph, & ! (out)
266 prog_vars, prog_vars0, & ! (in)
267 ddens00, momx00, momy00, momz00, drhot00, & ! (in)
268 dens_hyd, pres_hyd, & ! (in)
269 rtot, cptot_ov_cvtot, & ! (in)
270 dz, lift, intrpmat_vpordm1, & ! (in)
271 gnnm, g13, g23, gsqrtv, & ! (in)
272 impl_fac, dt, & ! (in)
273 lmesh, elem, & ! (in)
274 nz, vmapm, vmapp, & ! (in)
275 b1d_ij_uv ) ! (out)
276
277 implicit none
278
279 class(localmesh3d), intent(in) :: lmesh
280 class(elementbase3d), intent(in) :: elem
281 real(rp), intent(out) :: momx_t(elem%np,lmesh%nea)
282 real(rp), intent(out) :: momy_t(elem%np,lmesh%nea)
283 real(rp), intent(out) :: alph(elem%nfptot,lmesh%nez,lmesh%nex*lmesh%ney)
284 real(rp), intent(in) :: prog_vars (elem%np,lmesh%nez,prgvar_num,lmesh%nex*lmesh%ney)
285 real(rp), intent(in) :: prog_vars0 (elem%np,lmesh%nez,prgvar_num,lmesh%nex*lmesh%ney)
286 real(rp), intent(in) :: ddens00(elem%np,lmesh%nea)
287 real(rp), intent(in) :: momx00 (elem%np,lmesh%nea)
288 real(rp), intent(in) :: momy00 (elem%np,lmesh%nea)
289 real(rp), intent(in) :: momz00 (elem%np,lmesh%nea)
290 real(rp), intent(in) :: drhot00(elem%np,lmesh%nea)
291 real(rp), intent(in) :: dens_hyd(elem%np,lmesh%nez,lmesh%nex*lmesh%ney)
292 real(rp), intent(in) :: pres_hyd(elem%np,lmesh%nez,lmesh%nex*lmesh%ney)
293 real(rp), intent(in) :: rtot(elem%np,lmesh%nez,lmesh%nex*lmesh%ney)
294 real(rp), intent(in) :: cptot_ov_cvtot(elem%np,lmesh%nez,lmesh%nex*lmesh%ney)
295 class(sparsemat), intent(in) :: dz, lift
296 real(rp), intent(in) :: intrpmat_vpordm1(elem%np,elem%np)
297 real(rp), intent(in) :: gnnm(elem%np,lmesh%nez,lmesh%nex*lmesh%ney)
298 real(rp), intent(in) :: g13 (elem%np,lmesh%nez,lmesh%nex*lmesh%ney)
299 real(rp), intent(in) :: g23 (elem%np,lmesh%nez,lmesh%nex*lmesh%ney)
300 real(rp), intent(in) :: gsqrtv(elem%np,lmesh%nez,lmesh%nex*lmesh%ney)
301 real(rp), intent(in) :: impl_fac
302 real(rp), intent(in) :: dt
303 real(rp), intent(in) :: nz(elem%nfptot,lmesh%nez,lmesh%nex*lmesh%ney)
304 integer, intent(in) :: vmapm(elem%nfptot,lmesh%nez)
305 integer, intent(in) :: vmapp(elem%nfptot,lmesh%nez)
306 real(rp), intent(out), optional :: b1d_ij_uv(elem%nnode_v,lmesh%nez,2,elem%nnode_h1d**2,lmesh%nex*lmesh%ney)
307
308 real(rp) :: rgsqrtv(elem%np)
309 real(rp) :: fscale(elem%nfptot), escale33(elem%np)
310 real(rp) :: fz(elem%np), liftdelflx(elem%np)
311 real(rp) :: del_flux(elem%nfptot,lmesh%nez,lmesh%nex*lmesh%ney,prgvar_num)
312 integer :: ke_xy, ke_z
313 integer :: ke
314 integer :: v
315 integer :: ij
316 integer :: colmask(elem%nnode_v)
317
318 integer :: kk, kkk, p, pp
319 !--------------------------------------------------------
320
321 call vi_cal_del_flux_dyn_uv( del_flux, alph, & ! (out)
322 prog_vars, prog_vars0, & ! (in)
323 dens_hyd, pres_hyd, & ! (in)
324 rtot, cptot_ov_cvtot, & ! (in)
325 gnnm, g13, g23, gsqrtv, nz, vmapm, vmapp, & ! (in)
326 lmesh, elem ) ! (in)
327
328 !$omp parallel private( ke_xy, ke_z, ke, ij, v, &
329 !$omp Fz, LiftDelFlx, Fscale, RGsqrtV, ColMask, &
330 !$omp kk, p, kkk, pp )
331
332 !$omp do collapse(2)
333 do ke_xy=1, lmesh%NeX*lmesh%NeY
334 do ke_z=1, lmesh%NeZ
335 ke = ke_xy + (ke_z-1)*lmesh%NeX*lmesh%NeY
336
337 rgsqrtv(:) = 1.0_rp / gsqrtv(:,ke_z,ke_xy)
338 fscale(:) = lmesh%Fscale(:,ke)
339
340 !- MOMX
341 call sparsemat_matmul(lift, fscale(:) * del_flux(:,ke_z,ke_xy,momx_vid), liftdelflx)
342 momx_t(:,ke) = - liftdelflx(:) * rgsqrtv(:)
343
344 !-MOMY
345 call sparsemat_matmul(lift, fscale(:) * del_flux(:,ke_z,ke_xy,momy_vid), liftdelflx)
346 momy_t(:,ke) = - liftdelflx(:) * rgsqrtv(:)
347 end do
348 end do
349 !$omp end do
350
351 if ( present( b1d_ij_uv ) ) then
352 !$omp do collapse(2)
353 do ke_xy=1, lmesh%NeX*lmesh%NeY
354 do ke_z=1, lmesh%NeZ
355 ke = ke_xy + (ke_z-1)*lmesh%NeX*lmesh%NeY
356
357 do ij=1, elem%Nnode_h1D**2
358 colmask(:) = elem%Colmask(:,ij)
359 b1d_ij_uv(:,ke_z,1,ij,ke_xy) = impl_fac * momx_t(colmask(:),ke) &
360 - prog_vars(colmask(:),ke_z,momx_vid,ke_xy) &
361 + momx00(colmask(:),ke)
362 b1d_ij_uv(:,ke_z,2,ij,ke_xy) = impl_fac * momy_t(colmask(:),ke) &
363 - prog_vars(colmask(:),ke_z,momy_vid,ke_xy) &
364 + momy00(colmask(:),ke)
365 end do
366 end do
367 end do
368 !$omp end do
369 end if
370 !$omp end parallel
371
372 return
373 end subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_eval_ax_uv
374
375!OCL SERIAL
376 subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd( &
377 PmatBnd, & ! (out)
378 kl, ku, nz_1d, & ! (in)
379 prog_vars0, dens_hyd, pres_hyd, & ! (in)
380 g13, g23, gsqrtv, alph, & ! (in)
381 rtot, cptot_ov_cvtot, & ! (in)
382 dz, lift, intrpmat_vpordm1, & ! (in)
383 impl_fac, dt, & ! (in)
384 lmesh, elem, & ! (in)
385 nz, vmapm, vmapp, ke_x, ke_y ) ! (in)
386
387 implicit none
388
389 class(localmesh3d), intent(in) :: lmesh
390 class(elementbase3d), intent(in) :: elem
391 integer, intent(in) :: kl, ku, nz_1d
392 real(rp), intent(out) :: pmatbnd(2*kl+ku+1,3,elem%nnode_v,lmesh%nez,elem%nnode_h1d**2)
393 real(rp), intent(in) :: prog_vars0(elem%np,lmesh%nez,prgvar_num)
394 real(rp), intent(in) :: dens_hyd(elem%np,lmesh%nez)
395 real(rp), intent(in) :: pres_hyd(elem%np,lmesh%nez)
396 real(rp), intent(in) :: g13(elem%np,lmesh%nez)
397 real(rp), intent(in) :: g23(elem%np,lmesh%nez)
398 real(rp), intent(in) :: gsqrtv(elem%np,lmesh%nez)
399 real(rp), intent(in) :: alph(elem%nfptot,lmesh%nez)
400 real(rp), intent(in) :: rtot(elem%np,lmesh%nez)
401 real(rp), intent(in) :: cptot_ov_cvtot(elem%np,lmesh%nez)
402 class(sparsemat), intent(in) :: dz, lift
403 real(rp), intent(in) :: intrpmat_vpordm1(elem%np,elem%np)
404 real(rp), intent(in) :: impl_fac
405 real(rp), intent(in) :: dt
406 real(rp), intent(in) :: nz(elem%nfptot,lmesh%nez)
407 integer, intent(in) :: vmapm(elem%nfptot,lmesh%nez)
408 integer, intent(in) :: vmapp(elem%nfptot,lmesh%nez)
409 integer, intent(in) :: ke_x, ke_y
410
411 real(rp) :: rhot_hyd(elem%nnode_v)
412 real(rp) :: pot0(elem%nnode_v,lmesh%nez,elem%nnode_h1d**2)
413 real(rp) :: w0(elem%nnode_v,lmesh%nez,elem%nnode_h1d**2)
414 real(rp) :: wt0(elem%nnode_v,lmesh%nez,elem%nnode_h1d**2)
415 real(rp) :: dens0(elem%nnode_v,lmesh%nez,elem%nnode_h1d**2)
416 real(rp) :: dpdrhot0(elem%nnode_v,lmesh%nez,elem%nnode_h1d**2)
417 integer :: ke_z, ke_z2
418 integer :: v, ke, p, f1, f2, fp, fp2, fmv
419 real(rp) :: gamm, rgamm
420 real(rp) :: fac_dz_p(elem%nnode_v)
421 real(rp) :: pmatd(elem%nnode_v,elem%nnode_v,3,3)
422 real(rp) :: pmatl(elem%nnode_v,elem%nnode_v,3,3)
423 real(rp) :: pmatu(elem%nnode_v,elem%nnode_v,3,3)
424
425 integer :: colmask(elem%nnode_v)
426 real(rp) :: id(elem%nnode_v,elem%nnode_v)
427 real(rp) :: dd(elem%nnode_v)
428 real(rp) :: tmp1(elem%nnode_v)
429 real(rp) :: fac
430 integer :: ij, v1, v2, pv1, pv2, g_kj, g_kjp1, g_kjm1, pb1
431 logical :: bc_flag
432 logical :: eval_flag(3,3)
433
434 integer, parameter :: dens_vid_lc = 1
435 integer, parameter :: momz_vid_lc = 2
436 integer, parameter :: rhot_vid_lc = 3
437
438 !--------------------------------------------------------
439
440 gamm = cpdry/cvdry
441 rgamm = cvdry/cpdry
442
443 eval_flag(:,:) = .false.
444 do v=1, 3
445 eval_flag(v,v) = .true.
446 end do
447 eval_flag(dens_vid_lc,momz_vid_lc) = .true.
448 eval_flag(momz_vid_lc,dens_vid_lc) = .true.
449 eval_flag(momz_vid_lc,rhot_vid_lc) = .true.
450 eval_flag(rhot_vid_lc,momz_vid_lc) = .true.
451 eval_flag(rhot_vid_lc,dens_vid_lc) = .true.
452
453 id(:,:) = 0.0_rp
454 do p=1, elem%Nnode_v
455 id(p,p) = 1.0_rp
456 end do
457
458 !$omp parallel private(RHOT_hyd, Colmask)
459 !$omp workshare
460 pmatd(:,:,:,:) = 0.0_rp
461 pmatl(:,:,:,:) = 0.0_rp
462 pmatu(:,:,:,:) = 0.0_rp
463 !$omp end workshare
464 !$omp do
465 do ij=1, elem%Nnode_h1D**2
466 pmatbnd(:,:,:,:,ij) = 0.0_rp
467 end do
468 !$omp do collapse(2)
469 do ij=1, elem%Nnode_h1D**2
470 do ke_z=1, lmesh%NeZ
471 colmask(:) = elem%Colmask(:,ij)
472
473 rhot_hyd(:) = pres00/rdry * (pres_hyd(colmask(:),ke_z)/pres00)**rgamm
474 dens0(:,ke_z,ij) = dens_hyd(colmask(:),ke_z) + prog_vars0(colmask(:),ke_z,dens_vid)
475 pot0(:,ke_z,ij) = ( rhot_hyd(:) + prog_vars0(colmask(:),ke_z,rhot_vid) ) / dens0(:,ke_z,ij)
476 w0(:,ke_z,ij) = prog_vars0(colmask(:),ke_z,momz_vid) / dens0(:,ke_z,ij)
477
478 wt0(:,ke_z,ij) = w0(:,ke_z,ij) + gsqrtv(colmask(:),ke_z) * ( &
479 g13(colmask(:),ke_z) * prog_vars0(colmask(:),ke_z,momx_vid) &
480 + g23(colmask(:),ke_z) * prog_vars0(colmask(:),ke_z,momy_vid) ) / dens0(:,ke_z,ij)
481
482 ! DPDRHOT0(:,ke_z,ij) = gamm * PRES_hyd(Colmask(:),ke_z) / RHOT_hyd(:) &
483 ! * ( 1.0_RP + PROG_VARS0(Colmask(:),ke_z,RHOT_VID) / RHOT_hyd(:) )**(gamm-1.0_RP)
484 dpdrhot0(:,ke_z,ij) = cptot_ov_cvtot(colmask(:),ke_z) &
485 * pres00 * ( rtot(colmask(:),ke_z) / pres00 * ( dens0(:,ke_z,ij) * pot0(:,ke_z,ij) ) )**cptot_ov_cvtot(colmask(:),ke_z) &
486 / ( dens0(:,ke_z,ij) * pot0(:,ke_z,ij) )
487 end do
488 end do
489 !$omp end parallel
490
491 !$omp parallel private(ke_z, ke, ColMask, p, fp, fp2, v, f1, f2, ke_z2, fac_dz_p, &
492 !$omp fac, tmp1, FmV, &
493 !$omp ij, v1, v2, pv1, pv2, pb1, g_kj, g_kjp1, g_kjm1, bc_flag, &
494 !$omp Dd ) &
495 !$omp firstprivate(PmatD, PmatL, PmatU)
496
497 !$omp do collapse(2)
498 do ij=1, elem%Nnode_h1D**2
499 do ke_z=1, lmesh%NeZ
500 ke = ke_x + (ke_y-1)*lmesh%NeX + (ke_z-1)*lmesh%NeX*lmesh%NeY
501 colmask(:) = elem%Colmask(:,ij)
502
503 !-----
504 do p=1, elem%Nnode_v
505 fac_dz_p(:) = impl_fac * lmesh%Escale(colmask(:),ke,3,3) / gsqrtv(colmask(:),ke_z) &
506 * elem%Dx3(colmask(:),colmask(p))
507
508 dd(:) = id(:,p)
509
510 ! DDENS
511 pmatd(:,p,dens_vid_lc,dens_vid_lc) = dd(:)
512 pmatd(:,p,dens_vid_lc,momz_vid_lc) = fac_dz_p(:)
513
514 ! MOMZ
515 pmatd(:,p,momz_vid_lc,momz_vid_lc) = dd(:) !&
516 !+ fac_dz_p(:) * 2.0_RP * W0(p,ke_z,ij) [ <- d_MOMZ ( MOMZ x MOMZ / DENS ) ]
517 pmatd(:,p,momz_vid_lc,dens_vid_lc) = impl_fac * grav * intrpmat_vpordm1(colmask(:),colmask(p)) !&
518 !- fac_dz_p(:) * W0(p,ke_z,ij)**2 [ <- d_DENS ( MOMZ x MOMZ / DENS ) ]
519 pmatd(:,p,momz_vid_lc,rhot_vid_lc) = fac_dz_p(:) * dpdrhot0(p,ke_z,ij)
520
521 !DRHOT
522 pmatd(:,p,rhot_vid_lc,dens_vid_lc) = - fac_dz_p(:) * pot0(p,ke_z,ij) * wt0(p,ke_z,ij)
523 pmatd(:,p,rhot_vid_lc,momz_vid_lc) = fac_dz_p(:) * pot0(p,ke_z,ij)
524 pmatd(:,p,rhot_vid_lc,rhot_vid_lc) = dd(:) + fac_dz_p(:) * wt0(p,ke_z,ij)
525 end do
526
527 do f1=1, 2
528 if (f1==1) then
529 ke_z2 = max(ke_z-1,1)
530 pv1 = 1; pv2 = elem%Nnode_v
531 f2 = 2
532 else
533 ke_z2 = min(ke_z+1,lmesh%NeZ)
534 pv1 = elem%Nnode_v; pv2 = 1
535 f2 = 1
536 end if
537 fac = 0.5_rp * impl_fac / gsqrtv(colmask(pv1),ke_z)
538 if ( (ke_z == 1 .and. f1==1) .or. (ke_z == lmesh%NeZ .and. f1==elem%Nfaces_v) ) then
539 bc_flag = .true.
540 pv2 = pv1; f2 = f1
541 else
542 bc_flag = .false.
543 end if
544
545 fmv = elem%Fmask_v(ij,f1)
546 fp = elem%Nfp_h * elem%Nfaces_h + (f1-1)*elem%Nfp_v + ij
547 fp2 = elem%Nfp_h * elem%Nfaces_h + (f2-1)*elem%Nfp_v + ij
548
549 !--
550 tmp1(:) = fac * elem%Lift(colmask(:),fp) * lmesh%Fscale(fp,ke) &
551 * max( alph(fp,ke_z), alph(fp2,ke_z2) )
552 if (bc_flag) then
553 pmatd(:,pv1,momz_vid_lc,momz_vid_lc) = pmatd(:,pv1,momz_vid_lc,momz_vid_lc) + 2.0_rp * tmp1(:)
554 else
555 do v=1, 3
556 pmatd(:,pv1,v,v) = pmatd(:,pv1,v,v) + tmp1(:)
557 if (f1 == 1) then
558 pmatl(:,pv2,v,v) = - tmp1(:)
559 else
560 pmatu(:,pv2,v,v) = - tmp1(:)
561 end if
562 end do
563 end if
564
565 !--
566 tmp1(:) = fac * elem%Lift(colmask(:),fp) * lmesh%Fscale(fp,ke) * nz(fp,ke_z)
567
568 if (bc_flag) then
569 pmatd(:,pv1,dens_vid_lc,momz_vid_lc) = pmatd(:,pv1,dens_vid_lc,momz_vid_lc) - 2.0_rp * tmp1(:)
570 pmatd(:,pv1,rhot_vid_lc,momz_vid_lc) = pmatd(:,pv1,rhot_vid_lc,momz_vid_lc) - 2.0_rp * tmp1(:) * pot0(pv1,ke_z,ij)
571 pmatd(:,pv1,rhot_vid_lc,dens_vid_lc) = pmatd(:,pv1,rhot_vid_lc,dens_vid_lc) + 2.0_rp * tmp1(:) * pot0(pv1,ke_z,ij) * wt0(pv1,ke_z,ij)
572 pmatd(:,pv1,rhot_vid_lc,rhot_vid_lc) = pmatd(:,pv1,rhot_vid_lc,rhot_vid_lc) - 2.0_rp * tmp1(:) * wt0(pv1,ke_z,ij)
573 else
574 pmatd(:,pv1,dens_vid_lc,momz_vid_lc) = pmatd(:,pv1,dens_vid_lc,momz_vid_lc) - tmp1(:)
575
576 ! PmatD(pv1,pv1,MOMZ_VID_LC,MOMZ_VID_LC) = PmatD(pv1,pv1,MOMZ_VID_LC,MOMZ_VID_LC) - tmp1 * 2.0_RP * W0(pv1,ke_z,ij) [ <- d_MOMZ ( MOMZ x MOMZ / DENS ) ]
577 ! PmatD(pv1,pv1,MOMZ_VID_LC,DENS_VID_LC) = PmatD(pv1,pv1,MOMZ_VID_LC,DENS_VID_LC) + tmp1 * W0(pv1,ke_z,ij)**2 [ <- d_DENS ( MOMZ x MOMZ / DENS ) ]
578 pmatd(:,pv1,momz_vid_lc,rhot_vid_lc) = pmatd(:,pv1,momz_vid_lc,rhot_vid_lc) - tmp1(:) * dpdrhot0(pv1,ke_z,ij)
579
580 pmatd(:,pv1,rhot_vid_lc,momz_vid_lc) = pmatd(:,pv1,rhot_vid_lc,momz_vid_lc) - tmp1(:) * pot0(pv1,ke_z,ij)
581 pmatd(:,pv1,rhot_vid_lc,dens_vid_lc) = pmatd(:,pv1,rhot_vid_lc,dens_vid_lc) + tmp1(:) * pot0(pv1,ke_z,ij) * wt0(pv1,ke_z,ij)
582 pmatd(:,pv1,rhot_vid_lc,rhot_vid_lc) = pmatd(:,pv1,rhot_vid_lc,rhot_vid_lc) - tmp1(:) * wt0(pv1,ke_z,ij)
583
584 if (f1 == 1) then
585 ! PmatL(pv1,pv2,DENS_VID_LC,MOMZ_VID_LC) = + tmp1
586 pmatl(:,pv2,dens_vid_lc,momz_vid_lc) = + tmp1(:)
587
588 ! PmatL(pv1,pv2,MOMZ_VID_LC,MOMZ_VID_LC) = PmatL(pv1,pv2,MOMZ_VID_LC,MOMZ_VID_LC) & [ <- d_MOMZ ( MOMZ x MOMZ / DENS ) ]
589 ! + tmp1 * 2.0_RP * W0(pv2,ke_z2,ij)
590 ! PmatL(pv1,pv2,MOMZ_VID_LC,DENS_VID_LC) = - tmp1 * W0(pv2,ke_z2,ij) * W0(pv2,ke_z2,ij) [ <- d_DENS ( MOMZ x MOMZ / DENS ) ]
591 pmatl(:,pv2,momz_vid_lc,rhot_vid_lc) = + tmp1(:) * dpdrhot0(pv2,ke_z2,ij)
592
593 pmatl(:,pv2,rhot_vid_lc,momz_vid_lc) = + tmp1(:) * pot0(pv2,ke_z2,ij)
594 pmatl(:,pv2,rhot_vid_lc,dens_vid_lc) = - tmp1(:) * pot0(pv2,ke_z2,ij) * wt0(pv2,ke_z2,ij)
595 pmatl(:,pv2,rhot_vid_lc,rhot_vid_lc) = pmatl(:,pv2,rhot_vid_lc,rhot_vid_lc) &
596 + tmp1(:) * wt0(pv2,ke_z2,ij)
597 else
598 pmatu(:,pv2,dens_vid_lc,momz_vid_lc) = + tmp1(:)
599
600 ! PmatU(pv1,pv2,MOMZ_VID_LC,MOMZ_VID_LC) = PmatU(pv1,pv2,MOMZ_VID_LC,MOMZ_VID_LC ) & [ <- d_MOMZ ( MOMZ x MOMZ / DENS ) ]
601 ! + tmp1 * 2.0_RP * W0(pv2,ke_z2,ij)
602 ! PmatU(pv1,pv2,MOMZ_VID_LC,DENS_VID_LC) = - tmp1 * W0(pv2,ke_z2,ij) * W0(pv2,ke_z2,ij) [ <- d_DENS ( MOMZ x MOMZ / DENS ) ]
603 pmatu(:,pv2,momz_vid_lc,rhot_vid_lc) = + tmp1(:) * dpdrhot0(pv2,ke_z2,ij)
604
605 pmatu(:,pv2,rhot_vid_lc,momz_vid_lc) = + tmp1(:) * pot0(pv2,ke_z2,ij)
606 pmatu(:,pv2,rhot_vid_lc,dens_vid_lc) = - tmp1(:) * pot0(pv2,ke_z2,ij) * wt0(pv2,ke_z2,ij)
607 pmatu(:,pv2,rhot_vid_lc,rhot_vid_lc) = pmatu(:,pv2,rhot_vid_lc,rhot_vid_lc) &
608 + tmp1(:) * wt0(pv2,ke_z2,ij)
609 end if
610 end if
611 end do
612
613 do v2=1, 3
614 do v1=1, 3
615 if ( eval_flag(v1,v2) ) then
616 do pv2=1, elem%Nnode_v
617 g_kj = v2 + (pv2-1)*3 + (ke_z-1)*elem%Nnode_v*3
618 g_kjm1 = v2 + (pv2-1)*3 + (ke_z-2)*elem%Nnode_v*3
619 g_kjp1 = v2 + (pv2-1)*3 + (ke_z )*elem%Nnode_v*3
620
621 do pv1=1, elem%Nnode_v
622 pb1 = v1 + (pv1-1)*3 + (ke_z-1)*elem%Nnode_v*3
623
624 if (ke_z > 1 .and. pv2 == elem%Nnode_v ) then
625 pmatbnd(kl+ku+1+pb1-g_kjm1, v2,pv2,ke_z-1, ij) = pmatl(pv1,pv2,v1,v2)
626 end if
627 pmatbnd(kl+ku+1+pb1-g_kj, v2,pv2,ke_z, ij) = pmatd(pv1,pv2,v1,v2)
628 if (ke_z < lmesh%NeZ .and. pv2 == 1 ) then
629 pmatbnd(kl+ku+1+pb1-g_kjp1, v2,pv2,ke_z+1, ij) = pmatu(pv1,pv2,v1,v2)
630 end if
631 end do
632 end do
633 end if
634 end do
635 end do
636
637 end do
638 end do
639 !$omp end do
640 !$omp end parallel
641
642 return
643 end subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd
644
645!OCL SERIAL
646 subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd_uv( &
647 PmatBnd_uv, & ! (out)
648 kl_uv, ku_uv, nz_1d_uv, & ! (in)
649 prog_vars0, dens_hyd, pres_hyd, & ! (in)
650 g13, g23, gsqrtv, alph, & ! (in)
651 rtot, cptot_ov_cvtot, & ! (in)
652 dz, lift, intrpmat_vpordm1, & ! (in)
653 impl_fac, dt, & ! (in)
654 lmesh, elem, & ! (in)
655 nz, vmapm, vmapp, ke_x, ke_y ) ! (in)
656
657 implicit none
658
659 class(localmesh3d), intent(in) :: lmesh
660 class(elementbase3d), intent(in) :: elem
661 integer, intent(in) :: kl_uv, ku_uv, nz_1d_uv
662 real(rp), intent(out) :: pmatbnd_uv(2*kl_uv+ku_uv+1,elem%nnode_v,1,lmesh%nez,elem%nnode_h1d**2)
663 real(rp), intent(in) :: prog_vars0(elem%np,lmesh%nez,prgvar_num)
664 real(rp), intent(in) :: dens_hyd(elem%np,lmesh%nez)
665 real(rp), intent(in) :: pres_hyd(elem%np,lmesh%nez)
666 real(rp), intent(in) :: g13(elem%np,lmesh%nez)
667 real(rp), intent(in) :: g23(elem%np,lmesh%nez)
668 real(rp), intent(in) :: gsqrtv(elem%np,lmesh%nez)
669 real(rp), intent(in) :: alph(elem%nfptot,lmesh%nez)
670 real(rp), intent(in) :: rtot(elem%np,lmesh%nez)
671 real(rp), intent(in) :: cptot_ov_cvtot(elem%np,lmesh%nez)
672 class(sparsemat), intent(in) :: dz, lift
673 real(rp), intent(in) :: intrpmat_vpordm1(elem%np,elem%np)
674 real(rp), intent(in) :: impl_fac
675 real(rp), intent(in) :: dt
676 real(rp), intent(in) :: nz(elem%nfptot,lmesh%nez)
677 integer, intent(in) :: vmapm(elem%nfptot,lmesh%nez)
678 integer, intent(in) :: vmapp(elem%nfptot,lmesh%nez)
679 integer, intent(in) :: ke_x, ke_y
680
681 real(rp) :: rhot_hyd(elem%nnode_v)
682 real(rp) :: pot0(elem%nnode_v,lmesh%nez,elem%nnode_h1d**2)
683 real(rp) :: w0(elem%nnode_v,lmesh%nez,elem%nnode_h1d**2)
684 real(rp) :: wt0(elem%nnode_v,lmesh%nez,elem%nnode_h1d**2)
685 real(rp) :: dens0(elem%nnode_v,lmesh%nez,elem%nnode_h1d**2)
686 real(rp) :: dpdrhot0(elem%nnode_v,lmesh%nez,elem%nnode_h1d**2)
687 integer :: ke_z, ke_z2
688 integer :: v, ke, p, f1, f2, fp, fp2, fmv
689
690 real(rp) :: pmatd_uv(elem%nnode_v,elem%nnode_v)
691 real(rp) :: pmatl_uv(elem%nnode_v,elem%nnode_v)
692 real(rp) :: pmatu_uv(elem%nnode_v,elem%nnode_v)
693
694 integer :: colmask(elem%nnode_v)
695 real(rp) :: id(elem%nnode_v,elem%nnode_v)
696 real(rp) :: dd(elem%nnode_v)
697 real(rp) :: tmp1(elem%nnode_v)
698 real(rp) :: fac
699 integer :: ij, v1, v2, pv1, pv2, g_kj, g_kjp1, g_kjm1, pb1
700 logical :: bc_flag
701 !--------------------------------------------------------
702
703 id(:,:) = 0.0_rp
704 do p=1, elem%Nnode_v
705 id(p,p) = 1.0_rp
706 end do
707
708 !$omp parallel private(RHOT_hyd, Colmask)
709 !$omp workshare
710 pmatd_uv(:,:) = 0.0_rp
711 pmatl_uv(:,:) = 0.0_rp
712 pmatu_uv(:,:) = 0.0_rp
713 !$omp end workshare
714 !$omp do
715 do ij=1, elem%Nnode_h1D**2
716 pmatbnd_uv(:,:,:,:,ij) = 0.0_rp
717 end do
718 !$omp end parallel
719
720 !$omp parallel private(ke_z, ke, ColMask, p, fp, fp2, v, f1, f2, ke_z2, &
721 !$omp fac, tmp1, FmV, &
722 !$omp ij, v1, v2, pv1, pv2, pb1, g_kj, g_kjp1, g_kjm1, bc_flag, &
723 !$omp Dd ) &
724 !$omp firstprivate(PmatD_uv, PmatL_uv, PmatU_uv)
725
726 !$omp do collapse(2)
727 do ij=1, elem%Nnode_h1D**2
728 do ke_z=1, lmesh%NeZ
729 ke = ke_x + (ke_y-1)*lmesh%NeX + (ke_z-1)*lmesh%NeX*lmesh%NeY
730 colmask(:) = elem%Colmask(:,ij)
731
732 !-----
733 do p=1, elem%Nnode_v
734 dd(:) = id(:,p)
735 ! MOMX, MOMY
736 pmatd_uv(:,p) = dd(:)
737 end do
738
739 do f1=1, 2
740 if (f1==1) then
741 ke_z2 = max(ke_z-1,1)
742 pv1 = 1; pv2 = elem%Nnode_v
743 f2 = 2
744 else
745 ke_z2 = min(ke_z+1,lmesh%NeZ)
746 pv1 = elem%Nnode_v; pv2 = 1
747 f2 = 1
748 end if
749 fac = 0.5_rp * impl_fac / gsqrtv(colmask(pv1),ke_z)
750 if ( (ke_z == 1 .and. f1==1) .or. (ke_z == lmesh%NeZ .and. f1==elem%Nfaces_v) ) then
751 bc_flag = .true.
752 pv2 = pv1; f2 = f1
753 else
754 bc_flag = .false.
755 end if
756
757 fmv = elem%Fmask_v(ij,f1)
758 fp = elem%Nfp_h * elem%Nfaces_h + (f1-1)*elem%Nfp_v + ij
759 fp2 = elem%Nfp_h * elem%Nfaces_h + (f2-1)*elem%Nfp_v + ij
760
761 !--
762 tmp1(:) = fac * elem%Lift(colmask(:),fp) * lmesh%Fscale(fp,ke) &
763 * max( alph(fp,ke_z), alph(fp2,ke_z2) )
764 if (bc_flag) then
765 else
766 pmatd_uv(:,pv1) = pmatd_uv(:,pv1) + tmp1(:)
767 if (f1 == 1) then
768 pmatl_uv(:,pv2) = - tmp1(:)
769 else
770 pmatu_uv(:,pv2) = - tmp1(:)
771 end if
772 end if
773 end do
774
775 ! uv
776 do pv2=1, elem%Nnode_v
777 g_kj = pv2 + (ke_z-1)*elem%Nnode_v
778 g_kjm1 = pv2 + (ke_z-2)*elem%Nnode_v
779 g_kjp1 = pv2 + (ke_z )*elem%Nnode_v
780
781 do pv1=1, elem%Nnode_v
782 pb1 = pv1 + (ke_z-1)*elem%Nnode_v
783 if (ke_z > 1 .and. pv2 == elem%Nnode_v ) then
784 pmatbnd_uv(kl_uv+ku_uv+1+pb1-g_kjm1, pv2,1,ke_z-1, ij) = pmatl_uv(pv1,pv2)
785 end if
786 pmatbnd_uv(kl_uv+ku_uv+1+pb1-g_kj, pv2,1,ke_z, ij) = pmatd_uv(pv1,pv2)
787 if (ke_z < lmesh%NeZ .and. pv2 == 1) then
788 pmatbnd_uv(kl_uv+ku_uv+1+pb1-g_kjp1, pv2,1,ke_z+1, ij) = pmatu_uv(pv1,pv2)
789 end if
790 end do
791 end do
792
793 end do
794 end do
795 !$omp end do
796 !$omp end parallel
797
798 return
799 end subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd_uv
800
801!-- private ----------------
802
803!OCL SERIAL
804 subroutine vi_cal_del_flux_dyn_uv( del_flux, alph, & ! (out)
805 pvars_, pvars0_, & ! (in)
806 dens_hyd, pres_hyd, rtot, cptot_ov_cvtot, & ! (in)
807 gnn_, g13_, g23_, gsqrtv_, nz, vmapm, vmapp, lmesh, elem ) ! (in)
808
809 implicit none
810
811 class(localmesh3d), intent(in) :: lmesh
812 class(elementbase3d), intent(in) :: elem
813 real(rp), intent(out) :: del_flux(elem%nfptot*lmesh%nez,lmesh%nex*lmesh%ney,prgvar_num)
814 real(rp), intent(out) :: alph(elem%nfptot*lmesh%nez,lmesh%nex*lmesh%ney)
815 real(rp), intent(in) :: pvars_ (elem%np*lmesh%nez,prgvar_num,lmesh%nex*lmesh%ney)
816 real(rp), intent(in) :: pvars0_(elem%np*lmesh%nez,prgvar_num,lmesh%nex*lmesh%ney)
817 real(rp), intent(in) :: dens_hyd(elem%np*lmesh%nez,lmesh%nex*lmesh%ney)
818 real(rp), intent(in) :: pres_hyd(elem%np*lmesh%nez,lmesh%nex*lmesh%ney)
819 real(rp), intent(in) :: rtot(elem%np*lmesh%nez,lmesh%nex*lmesh%ney)
820 real(rp), intent(in) :: cptot_ov_cvtot(elem%np*lmesh%nez,lmesh%nex*lmesh%ney)
821 real(rp), intent(in) :: gnn_(elem%np*lmesh%nez,lmesh%nex*lmesh%ney)
822 real(rp), intent(in) :: g13_(elem%np*lmesh%nez,lmesh%nex*lmesh%ney)
823 real(rp), intent(in) :: g23_(elem%np*lmesh%nez,lmesh%nex*lmesh%ney)
824 real(rp), intent(in) :: gsqrtv_(elem%np*lmesh%nez,lmesh%nex*lmesh%ney)
825 real(rp), intent(in) :: nz(elem%nfptot*lmesh%nez,lmesh%nex*lmesh%ney)
826 integer, intent(in) :: vmapm(elem%nfptot*lmesh%nez)
827 integer, intent(in) :: vmapp(elem%nfptot*lmesh%nez)
828
829 integer :: i, p, ke_z, ip, im
830 integer :: ij
831 real(rp) :: wt0m, wt0p
832 real(rp) :: rhot_hyd_m, rhot_hyd_p
833 real(rp) :: pres0m, pres0p, dens0m, dens0p
834
835 real(rp) :: gamm, rgamm, pres0ovrdry, rp0
836 !------------------------------------------------------------------------
837
838 gamm = cpdry / cvdry
839 rgamm = cvdry / cpdry
840 pres0ovrdry = pres00 / rdry
841 rp0 = 1.0_rp / pres00
842
843 !$omp parallel private( ij, ke_z, p, i, iM, iP, &
844 !$omp rhot_hyd_M, rhot_hyd_P, &
845 !$omp wt0M, wt0P, dens0M, dens0P, pres0M, pres0P )
846
847 !$omp do collapse(2)
848 do ij=1, lmesh%NeX*lmesh%NeY
849 do ke_z=1, lmesh%NeZ
850 do p=1, elem%NfpTot
851 i = p + (ke_z-1)*elem%NfpTot
852 im = vmapm(i); ip = vmapp(i)
853
854 !-
855 rhot_hyd_m = pres0ovrdry * (pres_hyd(im,ij)/pres00)**rgamm
856 rhot_hyd_p = pres0ovrdry * (pres_hyd(ip,ij)/pres00)**rgamm
857
858 !-
859 dens0m = dens_hyd(im,ij) + pvars0_(im,dens_vid,ij)
860 dens0p = dens_hyd(ip,ij) + pvars0_(ip,dens_vid,ij)
861
862 pres0m = pres00 * ( rtot(im,ij) * rp0 * ( rhot_hyd_m + pvars0_(im,rhot_vid,ij) ) )**cptot_ov_cvtot(im,ij)
863 pres0p = pres00 * ( rtot(ip,ij) * rp0 * ( rhot_hyd_p + pvars0_(ip,rhot_vid,ij) ) )**cptot_ov_cvtot(ip,ij)
864
865 wt0m = ( pvars0_(im,momz_vid,ij) / gsqrtv_(im,ij) + g13_(im,ij) * pvars0_(im,momx_vid,ij) &
866 + g23_(im,ij) * pvars0_(im,momy_vid,ij) ) / dens0m
867 wt0p = ( pvars0_(ip,momz_vid,ij) / gsqrtv_(ip,ij) + g13_(ip,ij) * pvars0_(ip,momx_vid,ij) &
868 + g23_(ip,ij) * pvars0_(ip,momy_vid,ij) ) / dens0p
869
870 alph(i,ij) = nz(i,ij)**2 * max( abs( wt0m ) + sqrt( gnn_(im,ij) * gamm * pres0m / dens0m ), &
871 abs( wt0p ) + sqrt( gnn_(ip,ij) * gamm * pres0p / dens0p ) )
872
873 !----
874 del_flux(i,ij,momx_vid) = 0.5_rp * ( &
875 - alph(i,ij) * ( pvars_(ip,momx_vid,ij) - pvars_(im,momx_vid,ij) ) )
876
877 del_flux(i,ij,momy_vid) = 0.5_rp * ( &
878 - alph(i,ij) * ( pvars_(ip,momy_vid,ij) - pvars_(im,momy_vid,ij) ) )
879 end do
880 end do
881 end do
882 !$omp end do
883 !$omp end parallel
884
885 return
886 end subroutine vi_cal_del_flux_dyn_uv
887
888!OCL SERIAL
889 subroutine vi_cal_del_flux_dyn( del_flux, alph, & ! (out)
890 pvars_, pvars0_, & ! (in)
891 dens_hyd, pres_hyd, rtot, cptot_ov_cvtot, & ! (in)
892 gnn_, g13_, g23_, gsqrtv_, nz, vmapm, vmapp, lmesh, elem ) ! (in)
893
894 implicit none
895
896 class(localmesh3d), intent(in) :: lmesh
897 class(elementbase3d), intent(in) :: elem
898 real(rp), intent(out) :: del_flux(elem%nfptot*lmesh%nez,lmesh%nex*lmesh%ney,prgvar_num)
899 real(rp), intent(in) :: alph(elem%nfptot*lmesh%nez,lmesh%nex*lmesh%ney)
900 real(rp), intent(in) :: pvars_ (elem%np*lmesh%nez,prgvar_num,lmesh%nex*lmesh%ney)
901 real(rp), intent(in) :: pvars0_(elem%np*lmesh%nez,prgvar_num,lmesh%nex*lmesh%ney)
902 real(rp), intent(in) :: dens_hyd(elem%np*lmesh%nez,lmesh%nex*lmesh%ney)
903 real(rp), intent(in) :: pres_hyd(elem%np*lmesh%nez,lmesh%nex*lmesh%ney)
904 real(rp), intent(in) :: rtot(elem%np*lmesh%nez,lmesh%nex*lmesh%ney)
905 real(rp), intent(in) :: cptot_ov_cvtot(elem%np*lmesh%nez,lmesh%nex*lmesh%ney)
906 real(rp), intent(in) :: gnn_(elem%np*lmesh%nez,lmesh%nex*lmesh%ney)
907 real(rp), intent(in) :: g13_(elem%np*lmesh%nez,lmesh%nex*lmesh%ney)
908 real(rp), intent(in) :: g23_(elem%np*lmesh%nez,lmesh%nex*lmesh%ney)
909 real(rp), intent(in) :: gsqrtv_(elem%np*lmesh%nez,lmesh%nex*lmesh%ney)
910 real(rp), intent(in) :: nz(elem%nfptot*lmesh%nez,lmesh%nex*lmesh%ney)
911 integer, intent(in) :: vmapm(elem%nfptot*lmesh%nez)
912 integer, intent(in) :: vmapp(elem%nfptot*lmesh%nez)
913
914 integer :: i, p, ke_z, ip, im
915 integer :: ij
916 integer :: ke
917 real(rp) :: momz_p
918 real(rp) :: rhot_hyd_m, rhot_hyd_p
919 real(rp) :: dpresm, dpresp, densm, densp, pottm, pottp
920 real(rp) :: momw_p, momw_m
921
922 real(rp) :: rgamm, pres0ovrdry, rp0
923 !------------------------------------------------------------------------
924
925 rgamm = cvdry / cpdry
926 pres0ovrdry = pres00 / rdry
927 rp0 = 1.0_rp / pres00
928
929 !$omp parallel private( ij, ke_z, ke, p, i, iM, iP, &
930 !$omp rhot_hyd_M, rhot_hyd_P, densM, densP, pottM, pottP, &
931 !$omp dpresM, dpresP, MOMZ_P, MOMW_M, MOMW_P )
932
933 !$omp do collapse(2)
934 do ij=1, lmesh%NeX*lmesh%NeY
935 do ke_z=1, lmesh%NeZ
936 do p=1, elem%NfpTot
937 i = p + (ke_z-1)*elem%NfpTot
938 im = vmapm(i); ip = vmapp(i)
939
940 !-
941 rhot_hyd_m = pres0ovrdry * (pres_hyd(im,ij)/pres00)**rgamm
942 rhot_hyd_p = pres0ovrdry * (pres_hyd(ip,ij)/pres00)**rgamm
943
944 densm = dens_hyd(im,ij) + pvars_(im,dens_vid,ij)
945 densp = dens_hyd(ip,ij) + pvars_(ip,dens_vid,ij)
946
947 pottm = ( rhot_hyd_m + pvars_(im,rhot_vid,ij) ) / densm
948 pottp = ( rhot_hyd_p + pvars_(ip,rhot_vid,ij) ) / densp
949
950 dpresm = pres00 * ( rtot(im,ij) * rp0 * densm * pottm )**cptot_ov_cvtot(im,ij) &
951 - pres_hyd(im,ij)
952 dpresp = pres00 * ( rtot(ip,ij) * rp0 * densp * pottp )**cptot_ov_cvtot(ip,ij) &
953 - pres_hyd(ip,ij)
954
955 !-
956 momw_m = pvars_(im,momz_vid,ij) &
957 + gsqrtv_(im,ij) * g13_(im,ij) * pvars_(im,momx_vid,ij) &
958 + gsqrtv_(im,ij) * g23_(im,ij) * pvars_(im,momy_vid,ij)
959 momw_p = pvars_(ip,momz_vid,ij) &
960 + gsqrtv_(ip,ij) * g13_(ip,ij) * pvars_(ip,momx_vid,ij) &
961 + gsqrtv_(ip,ij) * g23_(ip,ij) * pvars_(ip,momy_vid,ij)
962
963 !----
964 if ( im == ip .and. (ke_z == 1 .or. ke_z == lmesh%NeZ) ) then
965 momz_p = - pvars_(im,momz_vid,ij) &
966 - 2.0_rp * gsqrtv_(im,ij) * ( g13_(im,ij) * pvars_(im,momx_vid,ij) &
967 + g23_(im,ij) * pvars_(im,momy_vid,ij) )
968 momw_p = - momw_m
969 else
970 momz_p = pvars_(ip,momz_vid,ij)
971 end if
972
973 del_flux(i,ij,dens_vid) = 0.5_rp * ( &
974 + ( momw_p - momw_m ) * nz(i,ij) &
975 - alph(i,ij) * ( pvars_(ip,dens_vid,ij) - pvars_(im,dens_vid,ij) ) )
976
977
978 del_flux(i,ij,momz_vid) = 0.5_rp * ( &
979 + ( dpresp - dpresm ) * nz(i,ij) &
980 - alph(i,ij) * ( momz_p - pvars_(im,momz_vid,ij) ) )
981
982 del_flux(i,ij,rhot_vid) = 0.5_rp * ( &
983 + ( pottp * momw_p - pottm * momw_m ) * nz(i,ij) &
984 - alph(i,ij) * ( pvars_(ip,rhot_vid,ij) - pvars_(im,rhot_vid,ij) ) )
985 end do
986 end do
987 end do
988 !$omp end do
989 !$omp end parallel
990
991 return
992 end subroutine vi_cal_del_flux_dyn
993end module scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common
scale_sparsemat::sparsemat_matmul
Definition scale_sparsemat.F90:57
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_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_drhot_id
integer, parameter, public prgvar_drhot_id
Definition scale_atm_dyn_dgm_nonhydro3d_common.F90:70
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::intrpmat_vpordm1
real(rp), dimension(:,:), allocatable, public intrpmat_vpordm1
Definition scale_atm_dyn_dgm_nonhydro3d_common.F90:102
scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common
module FElib / Fluid dyn solver / Atmosphere / Nonhydrostatic model / HEVI / Common
Definition scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common.F90:12
scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common::atm_dyn_dgm_nonhydro3d_rhot_hevi_common_init
subroutine, public atm_dyn_dgm_nonhydro3d_rhot_hevi_common_init
Definition scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common.F90:78
scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common::atm_dyn_dgm_nonhydro3d_rhot_hevi_common_eval_ax
subroutine, public atm_dyn_dgm_nonhydro3d_rhot_hevi_common_eval_ax(dens_t, momz_t, rhot_t, alph, prog_vars, prog_vars0, ddens00, momx00, momy00, momz00, drhot00, dens_hyd, pres_hyd, rtot, cptot_ov_cvtot, dz, lift, intrpmat_vpordm1, gnnm, g13, g23, gsqrtv, impl_fac, dt, lmesh, elem, nz, vmapm, vmapp, b1d_ij)
Definition scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common.F90:119
scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common::atm_dyn_dgm_nonhydro3d_rhot_hevi_common_eval_ax_uv
subroutine, public atm_dyn_dgm_nonhydro3d_rhot_hevi_common_eval_ax_uv(momx_t, momy_t, alph, prog_vars, prog_vars0, ddens00, momx00, momy00, momz00, drhot00, dens_hyd, pres_hyd, rtot, cptot_ov_cvtot, dz, lift, intrpmat_vpordm1, gnnm, g13, g23, gsqrtv, impl_fac, dt, lmesh, elem, nz, vmapm, vmapp, b1d_ij_uv)
Definition scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common.F90:276
scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common::atm_dyn_dgm_nonhydro3d_rhot_hevi_common_final
subroutine, public atm_dyn_dgm_nonhydro3d_rhot_hevi_common_final
Definition scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common.F90:100
scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common::atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd
subroutine, public atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd(pmatbnd, kl, ku, nz_1d, prog_vars0, dens_hyd, pres_hyd, g13, g23, gsqrtv, alph, rtot, cptot_ov_cvtot, dz, lift, intrpmat_vpordm1, impl_fac, dt, lmesh, elem, nz, vmapm, vmapp, ke_x, ke_y)
Definition scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common.F90:386
scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common::vi_use_lapack_flag
logical, public vi_use_lapack_flag
Definition scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common.F90:64
scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common::atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd_uv
subroutine, public atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd_uv(pmatbnd_uv, kl_uv, ku_uv, nz_1d_uv, prog_vars0, dens_hyd, pres_hyd, g13, g23, gsqrtv, alph, rtot, cptot_ov_cvtot, dz, lift, intrpmat_vpordm1, impl_fac, dt, lmesh, elem, nz, vmapm, vmapp, ke_x, ke_y)
Definition scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common.F90:656
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_element_modalfilter
module FElib / Element/ ModalFilter
Definition scale_element_modalfilter.F90:11
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_meshfield_base
module FElib / Data / base
Definition scale_meshfield_base.F90:10
scale_sparsemat
Module common / sparsemat.
Definition scale_sparsemat.F90:11
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_element_modalfilter::modalfilter
Derived type representing a modal filter.
Definition scale_element_modalfilter.F90:42
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_meshfield_base::meshfield3d
Derived type representing a field with 3D mesh.
Definition scale_meshfield_base.F90:95
scale_sparsemat::sparsemat
Derived type to manage a sparse matrix.
Definition scale_sparsemat.F90:33
Generated by doxygen 1.14.0