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