FE-Project
Loading...
Searching...
No Matches
scale_atm_dyn_dgm_nonhydro3d_rhot_heve_numflux.F90
Go to the documentation of this file.
1!-------------------------------------------------------------------------------
2!> module FElib / Fluid dyn solver / Atmosphere / Nonhydrostatic model / HEVE / Numflux
3!!
4!! @par Description
5!! HEVE DGM scheme for Atmospheric dynamical process.
6!!
7!! @author Yuta Kawai, Team SCALE
8!<
9!-------------------------------------------------------------------------------
10#include "scaleFElib.h"
12 !-----------------------------------------------------------------------------
13 !
14 !++ Used modules
15 !
16 use scale_precision
17 use scale_io
18 use scale_prc
19 use scale_prof
20 use scale_const, only: &
21 grav => const_grav, &
22 rdry => const_rdry, &
23 cpdry => const_cpdry, &
24 cvdry => const_cvdry, &
25 pres00 => const_pre00
26
28 use scale_element_base, only: &
36
38 dens_vid => prgvar_ddens_id, rhot_vid => prgvar_drhot_id, &
39 momx_vid => prgvar_momx_id, momy_vid => prgvar_momy_id, &
40 momz_vid => prgvar_momz_id, &
42
43 !-----------------------------------------------------------------------------
44 implicit none
45 private
46 !-----------------------------------------------------------------------------
47 !
48 !++ Public procedures
49 !
52! public :: atm_dyn_dgm_nonhydro3d_rhot_heve_numflux_get_generalvc_LMARS
56 ! public :: atm_dyn_dgm_nonhydro3d_rhot_heve_add_bnd_contrib2_generalvc
57
58 !-----------------------------------------------------------------------------
59 !
60 !++ Public parameters & variables
61 !
62
63 !-----------------------------------------------------------------------------
64 !
65 !++ Private procedures & variables
66 !
67 !-------------------
68
69contains
70
71!OCL SERIAL
73 del_flux, del_flux_hyd, & ! (out)
74 ddens_, momx_, momy_, momz_, drhot_, dpres_, dens_hyd, pres_hyd, & ! (in)
75 rtot, cvtot, cptot, & ! (in)
76 gsqrt, g13, g23, nx, ny, nz, & ! (in)
77 vmapm, vmapp, lmesh, elem, lmesh2d, elem2d ) ! (in)
78
79 implicit none
80
81 class(localmesh3d), intent(in) :: lmesh
82 class(elementbase3d), intent(in) :: elem
83 class(localmesh2d), intent(in) :: lmesh2d
84 class(elementbase2d), intent(in) :: elem2d
85 real(rp), intent(out) :: del_flux(elem%nfptot,lmesh%ne,prgvar_num)
86 real(rp), intent(out) :: del_flux_hyd(elem%nfptot,lmesh%ne,2)
87 real(rp), intent(in) :: ddens_(elem%np*lmesh%nea)
88 real(rp), intent(in) :: momx_(elem%np*lmesh%nea)
89 real(rp), intent(in) :: momy_(elem%np*lmesh%nea)
90 real(rp), intent(in) :: momz_(elem%np*lmesh%nea)
91 real(rp), intent(in) :: drhot_(elem%np*lmesh%nea)
92 real(rp), intent(in) :: dpres_(elem%np*lmesh%nea)
93 real(rp), intent(in) :: dens_hyd(elem%np*lmesh%nea)
94 real(rp), intent(in) :: pres_hyd(elem%np*lmesh%nea)
95 real(rp), intent(in) :: rtot (elem%np*lmesh%nea)
96 real(rp), intent(in) :: cvtot(elem%np*lmesh%nea)
97 real(rp), intent(in) :: cptot(elem%np*lmesh%nea)
98 real(rp), intent(in) :: gsqrt(elem%np*lmesh%nea)
99 real(rp), intent(in) :: g13(elem%np*lmesh%nea)
100 real(rp), intent(in) :: g23(elem%np*lmesh%nea)
101 real(rp), intent(in) :: nx(elem%nfptot,lmesh%ne)
102 real(rp), intent(in) :: ny(elem%nfptot,lmesh%ne)
103 real(rp), intent(in) :: nz(elem%nfptot,lmesh%ne)
104 integer, intent(in) :: vmapm(elem%nfptot,lmesh%ne)
105 integer, intent(in) :: vmapp(elem%nfptot,lmesh%ne)
106
107 integer :: ke, i, ip(elem%nfptot), im(elem%nfptot)
108 integer :: ke2d
109 real(rp) :: velp(elem%nfptot), velm(elem%nfptot), alpha(elem%nfptot)
110 real(rp) :: dpresp(elem%nfptot), dpresm(elem%nfptot)
111 real(rp) :: gsqrtdensm(elem%nfptot), gsqrtdensp(elem%nfptot)
112 real(rp) :: gsqrtrhotm(elem%nfptot), gsqrtrhotp(elem%nfptot)
113 real(rp) :: gsqrtddens_p(elem%nfptot), gsqrtddens_m(elem%nfptot)
114 real(rp) :: gsqrtmomx_p(elem%nfptot), gsqrtmomx_m(elem%nfptot)
115 real(rp) :: gsqrtmomy_p(elem%nfptot), gsqrtmomy_m(elem%nfptot)
116 real(rp) :: gsqrtmomz_p(elem%nfptot), gsqrtmomz_m(elem%nfptot)
117 real(rp) :: gsqrtdrhot_p(elem%nfptot), gsqrtdrhot_m(elem%nfptot)
118 real(rp) :: phyd_p(elem%nfptot), phyd_m(elem%nfptot)
119 real(rp) :: gsqrt_p(elem%nfptot), gsqrt_m(elem%nfptot)
120 real(rp) :: gsqrtv_p(elem%nfptot), gsqrtv_m(elem%nfptot)
121 real(rp) :: g13_m(elem%nfptot), g13_p(elem%nfptot)
122 real(rp) :: g23_m(elem%nfptot), g23_p(elem%nfptot)
123 real(rp) :: gnn_m(elem%nfptot), gnn_p(elem%nfptot)
124
125 real(rp) :: gamm, rgamm
126 real(rp) :: rp0
127 real(rp) :: rovp0, p0ovr
128 !------------------------------------------------------------------------
129
130 gamm = cpdry / cvdry
131 rgamm = cvdry / cpdry
132 rp0 = 1.0_rp / pres00
133 rovp0 = rdry * rp0
134 p0ovr = pres00 / rdry
135
136 !$omp parallel do private( &
137 !$omp ke, iM, iP, ke2D, &
138 !$omp alpha, VelM, VelP, &
139 !$omp dpresM, dpresP, GsqrtDensM, GsqrtDensP, GsqrtRhotM, GsqrtRhotP, &
140 !$omp GsqrtMOMX_M, GsqrtMOMX_P, GsqrtMOMY_M, GsqrtMOMY_P, GsqrtMOMZ_M, GsqrtMOMZ_P, &
141 !$omp GsqrtDDENS_M, GsqrtDDENS_P, GsqrtDRHOT_M, GsqrtDRHOT_P, &
142 !$omp Phyd_M, Phyd_P, &
143 !$omp Gsqrt_P, Gsqrt_M, GsqrtV_P, GsqrtV_M, G13_P, G13_M, G23_P, G23_M, &
144 !$omp Gnn_P, Gnn_M )
145 do ke=lmesh%NeS, lmesh%NeE
146 im(:) = vmapm(:,ke); ip(:) = vmapp(:,ke)
147 ke2d = lmesh%EMap3Dto2D(ke)
148
149 gsqrt_m(:) = gsqrt(im)
150 gsqrt_p(:) = gsqrt(ip)
151 gsqrtv_m(:) = gsqrt_m(:)
152 gsqrtv_p(:) = gsqrt_p(:)
153
154 g13_m(:) = g13(im)
155 g13_p(:) = g13(ip)
156 g23_m(:) = g23(im)
157 g23_p(:) = g23(ip)
158
159 gsqrtddens_m(:) = gsqrt_m(:) * ddens_(im)
160 gsqrtddens_p(:) = gsqrt_p(:) * ddens_(ip)
161 gsqrtmomx_m(:) = gsqrt_m(:) * momx_(im)
162 gsqrtmomx_p(:) = gsqrt_p(:) * momx_(ip)
163 gsqrtmomy_m(:) = gsqrt_m(:) * momy_(im)
164 gsqrtmomy_p(:) = gsqrt_p(:) * momy_(ip)
165 gsqrtmomz_m(:) = gsqrt_m(:) * momz_(im)
166 gsqrtmomz_p(:) = gsqrt_p(:) * momz_(ip)
167 gsqrtdrhot_m(:) = gsqrt_m(:) * drhot_(im)
168 gsqrtdrhot_p(:) = gsqrt_p(:) * drhot_(ip)
169 phyd_m(:) = pres_hyd(im)
170 phyd_p(:) = pres_hyd(ip)
171
172 gnn_m(:) = abs( nx(:,ke) ) + abs( ny(:,ke) ) &
173 + ( 1.0_rp / gsqrtv_m(:)**2 + g13_m(:)**2 + g23_m(:)**2 ) * abs( nz(:,ke) )
174 gnn_p(:) = abs( nx(:,ke) ) + abs( ny(:,ke) ) &
175 + ( 1.0_rp / gsqrtv_p(:)**2 + g13_p(:)**2 + g23_p(:)**2 ) * abs( nz(:,ke) )
176
177 gsqrtdensm(:) = gsqrtddens_m(:) + gsqrt_m(:) * dens_hyd(im)
178 gsqrtdensp(:) = gsqrtddens_p(:) + gsqrt_p(:) * dens_hyd(ip)
179
180 gsqrtrhotm(:) = gsqrt_m(:) * p0ovr * (phyd_m(:) * rp0)**rgamm + gsqrtdrhot_m(:)
181 gsqrtrhotp(:) = gsqrt_p(:) * p0ovr * (phyd_p(:) * rp0)**rgamm + gsqrtdrhot_p(:)
182
183 velm(:) = ( gsqrtmomx_m(:) * nx(:,ke) + gsqrtmomy_m(:) * ny(:,ke) &
184 + ( ( gsqrtmomz_m(:) / gsqrtv_m(:) &
185 + g13_m(:) * gsqrtmomx_m(:) + g23_m(:) * gsqrtmomy_m(:) ) * nz(:,ke) ) &
186 ) / gsqrtdensm(:)
187 velp(:) = ( gsqrtmomx_p(:) * nx(:,ke) + gsqrtmomy_p(:) * ny(:,ke) &
188 + ( ( gsqrtmomz_p(:) / gsqrtv_p(:) &
189 + g13_p(:) * gsqrtmomx_p(:) + g23_p(:) * gsqrtmomy_p(:) ) * nz(:,ke) ) &
190 ) / gsqrtdensp(:)
191
192
193 ! dpresM(:) = PRES00 * ( Rtot(iM) * rP0 * GsqrtRhotM(:) / Gsqrt_M(:) )**( CPtot(iM) / CVtot(iM) ) &
194 ! - Phyd_M(:)
195 ! dpresP(:) = PRES00 * ( Rtot(iP) * rP0 * GsqrtRhotP(:) / Gsqrt_P(:) )**( CPtot(iP) / CVtot(iP) ) &
196 ! - Phyd_P(:)
197 dpresm(:) = dpres_(im)
198 dpresp(:) = dpres_(ip)
199
200 alpha(:) = max( sqrt( gnn_m(:) * gamm * ( phyd_m(:) + dpresm(:) ) * gsqrt_m(:) / gsqrtdensm(:) ) + abs(velm(:)), &
201 sqrt( gnn_p(:) * gamm * ( phyd_p(:) + dpresp(:) ) * gsqrt_p(:) / gsqrtdensp(:) ) + abs(velp(:)) )
202
203 del_flux(:,ke,dens_vid) = 0.5_rp * ( &
204 ( gsqrtdensp(:) * velp(:) - gsqrtdensm(:) * velm(:) ) &
205 - alpha(:) * ( gsqrtddens_p(:) - gsqrtddens_m(:) ) )
206
207 del_flux(:,ke,momx_vid ) = 0.5_rp * ( &
208 ( gsqrtmomx_p(:) * velp(:) - gsqrtmomx_m(:) * velm(:) ) &
209 + ( gsqrt_p(:) * ( nx(:,ke) + g13_p(:) * nz(:,ke)) * dpresp(:) &
210 - gsqrt_m(:) * ( nx(:,ke) + g13_m(:) * nz(:,ke)) * dpresm(:) ) &
211 - alpha(:) * ( gsqrtmomx_p(:) - gsqrtmomx_m(:) ) )
212
213 del_flux(:,ke,momy_vid ) = 0.5_rp * ( &
214 ( gsqrtmomy_p(:) * velp(:) - gsqrtmomy_m(:) * velm(:) ) &
215 + ( gsqrt_p(:) * ( ny(:,ke) + g23_p(:) * nz(:,ke)) * dpresp(:) &
216 - gsqrt_m(:) * ( ny(:,ke) + g23_m(:) * nz(:,ke)) * dpresm(:) ) &
217 - alpha(:) * ( gsqrtmomy_p(:) - gsqrtmomy_m(:) ) )
218
219 del_flux(:,ke,momz_vid ) = 0.5_rp * ( &
220 ( gsqrtmomz_p(:) * velp(:) - gsqrtmomz_m(:) * velm(:) ) &
221 + ( gsqrt_p(:) * dpresp(:) / gsqrtv_p(:) &
222 - gsqrt_m(:) * dpresm(:) / gsqrtv_m(:) ) * nz(:,ke) &
223 - alpha(:) * ( gsqrtmomz_p(:) - gsqrtmomz_m(:) ) )
224
225 del_flux(:,ke,rhot_vid) = 0.5_rp * ( &
226 ( gsqrtrhotp(:) * velp(:) - gsqrtrhotm(:) * velm(:) ) &
227 - alpha(:) * ( gsqrtdrhot_p(:) - gsqrtdrhot_m(:) ) )
228
229 del_flux_hyd(:,ke,1) = 0.5_rp * ( &
230 gsqrtv_p(:) * ( nx(:,ke) + g13_p(:) * nz(:,ke) ) * phyd_p(:) &
231 - gsqrtv_m(:) * ( nx(:,ke) + g13_m(:) * nz(:,ke) ) * phyd_m(:) )
232
233 del_flux_hyd(:,ke,2) = 0.5_rp * ( &
234 gsqrtv_p(:) * ( ny(:,ke) + g23_p(:) * nz(:,ke) ) * phyd_p(:) &
235 - gsqrtv_m(:) * ( ny(:,ke) + g23_m(:) * nz(:,ke) ) * phyd_m(:) )
236 end do
237
238 return
240
241!OCL SERIAL
243 PRGVAR_dt, & ! (inout)
244 ddens_, momx_, momy_, momz_, drhot_, dpres, & ! (in)
245 dens_hyd, pres_hyd, therm_hyd, rtot, cvtot, cptot, & ! (in)
246 gsqrt, g13, g23, nx, ny, nz, & ! (in)
247 vmapm, vmapp, lmesh, elem, lmesh2d, elem2d, elem3d_optr ) ! (in)
249 implicit none
250
251 class(localmesh3d), intent(in) :: lmesh
252 class(elementbase3d), intent(in) :: elem
253 class(localmesh2d), intent(in) :: lmesh2d
254 class(elementbase2d), intent(in) :: elem2d
255 real(rp), intent(inout) :: prgvar_dt(elem%np,lmesh%nea,5)
256 real(rp), intent(in) :: ddens_(elem%np*lmesh%nea)
257 real(rp), intent(in) :: momx_(elem%np*lmesh%nea)
258 real(rp), intent(in) :: momy_(elem%np*lmesh%nea)
259 real(rp), intent(in) :: momz_(elem%np*lmesh%nea)
260 real(rp), intent(in) :: drhot_(elem%np*lmesh%nea)
261 real(rp), intent(inout) :: dpres(elem%np*lmesh%nea)
262 real(rp), intent(in) :: dens_hyd(elem%np*lmesh%nea)
263 real(rp), intent(in) :: pres_hyd(elem%np*lmesh%nea)
264 real(rp), intent(in) :: therm_hyd(elem%np*lmesh%nea)
265 real(rp), intent(in) :: rtot (elem%np*lmesh%nea)
266 real(rp), intent(in) :: cvtot(elem%np*lmesh%nea)
267 real(rp), intent(in) :: cptot(elem%np*lmesh%nea)
268 real(rp), intent(in) :: gsqrt(elem%np*lmesh%nea)
269 real(rp), intent(in) :: g13(elem%np*lmesh%nea)
270 real(rp), intent(in) :: g23(elem%np*lmesh%nea)
271 real(rp), intent(in) :: nx(elem%nfptot,lmesh%ne)
272 real(rp), intent(in) :: ny(elem%nfptot,lmesh%ne)
273 real(rp), intent(in) :: nz(elem%nfptot,lmesh%ne)
274 integer, intent(in) :: vmapm(elem%nfptot,lmesh%ne)
275 integer, intent(in) :: vmapp(elem%nfptot,lmesh%ne)
276 class(elementoperationbase3d), intent(in) :: elem3d_optr
277
278 integer :: ke, fp, i, ip(elem%nfptot), im(elem%nfptot)
279 integer :: ke2d
280 real(rp) :: vel(elem%nfptot,2), alpha(elem%nfptot)
281 real(rp) :: dpres_(elem%nfptot,2)
282 real(rp) :: gsqrtdens(elem%nfptot,2)
283 real(rp) :: gsqrtrhot(elem%nfptot,2)
284 real(rp) :: gsqrtddens(elem%nfptot,2)
285 real(rp) :: gsqrtmomx(elem%nfptot,2)
286 real(rp) :: gsqrtmomy(elem%nfptot,2)
287 real(rp) :: gsqrtmomz(elem%nfptot,2)
288 real(rp) :: gsqrtdrhot(elem%nfptot,2)
289 real(rp) :: cs2_(elem%nfptot,2)
290 real(rp) :: gsqrt_(elem%nfptot,2)
291 real(rp) :: gsqrtv_(elem%nfptot,2)
292 real(rp) :: rgsqrtv(elem%nfptot,2)
293 real(rp) :: g13_(elem%nfptot,2)
294 real(rp) :: g23_(elem%nfptot,2)
295 real(rp) :: gnn_m, gnn_p
296
297 real(rp) :: gamm, rgamm
298 real(rp) :: rp0
299 real(rp) :: rovp0, p0ovr
300
301 integer, parameter :: in = 1
302 integer, parameter :: ex = 2
303
304 real(rp) :: tmp1, tmp2, tmp3, tmp4, tmp01, tmp02, tmp03
305 real(rp) :: del_flux_tmp_mom(elem%nfptot,3)
306
307 integer :: iv
308 integer :: p
309 real(rp) :: del_flux(elem%nfptot,5)
310 real(rp) :: lift(elem%np,5)
311 real(rp) :: rgsqrt(elem%np)
312
313 integer :: iim, iip, fpp
314 integer, parameter :: nfp_s = 0
315 !------------------------------------------------------------------------
316
317 gamm = cpdry / cvdry
318 rgamm = cvdry / cpdry
319 rp0 = 1.0_rp / pres00
320 rovp0 = rdry * rp0
321 p0ovr = pres00 / rdry
322 call prof_rapstart('cal_dyn_tend_bndflux1', 3)
323
324 !$omp parallel private( &
325 !$omp ke, iM, iP, ke2D, fp, i, iim, iip, fpp, &
326 !$omp alpha, Vel, &
327 !$omp dpres_, GsqrtDens, GsqrtRhot, &
328 !$omp GsqrtDDENS, GsqrtMOMX, GsqrtMOMY, GsqrtMOMZ, GsqrtDRHOT, &
329 !$omp Cs2_, &
330 !$omp Gsqrt_, GsqrtV_, RGsqrtV, G13_, G23_, &
331 !$omp Gnn_P, Gnn_M, tmp1, tmp2, tmp3, tmp4, tmp01, tmp02, tmp03, del_flux_tmp_mom, &
332 !$omp del_flux, p, iv, Lift, RGsqrt )
333
334 !$omp do
335 do i=elem%Np*lmesh%NeE+1, elem%Np*lmesh%NeE+size(lmesh%VMapB)
336 dpres(i) = pres00 * ( rtot(i) * rp0 * ( therm_hyd(i) + drhot_(i) ) )**( cptot(i) / cvtot(i) ) &
337 - pres_hyd(i)
338 end do
339 !$omp end do
340!OCL PREFETCH
341 !$omp do
342 do ke=lmesh%NeS, lmesh%NeE
343 im(:) = vmapm(:,ke); ip(:) = vmapp(:,ke)
344 ke2d = lmesh%EMap3Dto2D(ke)
345
346 do fp=1, elem%NfpTot
347 iim = im(fp)
348 gsqrt_(fp,in) = gsqrt(iim)
349 gsqrtv_(fp,in) = gsqrt_(fp,in)
350 rgsqrtv(fp,in) = 1.0_rp / gsqrtv_(fp,in)
351
352 gsqrtddens(fp,in) = gsqrt_(fp,in) * ddens_(iim)
353 gsqrtdens(fp,in) = gsqrtddens(fp,in) + gsqrt_(fp,in) * dens_hyd(iim)
354
355 gsqrtdrhot(fp,in) = gsqrt_(fp,in) * drhot_(iim)
356 gsqrtrhot(fp,in) = gsqrtdrhot(fp,in) + gsqrt_(fp,in) * therm_hyd(iim)
357 end do
358 do fp=1, elem%NfpTot
359 iim = im(fp)
360
361 dpres_(fp,in) = dpres(iim)
362 cs2_(fp,in) = gamm * ( pres_hyd(iim) + dpres_(fp,in) ) * gsqrt_(fp,in) / gsqrtdens(fp,in)
363 end do
364 do fp=1, elem%NfpTot
365 iim = im(fp)
366 gsqrtmomx(fp,in) = gsqrt_(fp,in) * momx_(iim)
367 gsqrtmomy(fp,in) = gsqrt_(fp,in) * momy_(iim)
368 gsqrtmomz(fp,in) = gsqrt_(fp,in) * momz_(iim)
369
370 g13_(fp,in) = g13(iim)
371 g23_(fp,in) = g23(iim)
372 end do
373 do fp=1, elem%NfpTot
374 fpp = nfp_s + fp
375 vel(fp,in) = ( gsqrtmomx(fp,in) * nx(fpp,ke) + gsqrtmomy(fp,in) * ny(fpp,ke) &
376 + ( ( gsqrtmomz(fp,in) * rgsqrtv(fp,in) &
377 + g13_(fp,in) * gsqrtmomx(fp,in) + g23_(fp,in) * gsqrtmomy(fp,in) ) * nz(fpp,ke) ) &
378 ) / gsqrtdens(fp,in)
379 end do
380
381 !-
382 do fp=1, elem%NfpTot
383 iip = ip(fp)
384
385 gsqrt_(fp,ex) = gsqrt(iip)
386 gsqrtv_(fp,ex) = gsqrt_(fp,ex)
387 rgsqrtv(fp,ex) = 1.0_rp / gsqrtv_(fp,ex)
388
389
390 gsqrtddens(fp,ex) = gsqrt_(fp,ex) * ddens_(iip)
391 gsqrtdens(fp,ex) = gsqrtddens(fp,ex) + gsqrt_(fp,ex) * dens_hyd(iip)
392
393 gsqrtdrhot(fp,ex) = gsqrt_(fp,ex) * drhot_(iip)
394 gsqrtrhot(fp,ex) = gsqrtdrhot(fp,ex) + gsqrt_(fp,ex) * therm_hyd(iip)
395 end do
396 do fp=1, elem%NfpTot
397 iip = ip(fp)
398 dpres_(fp,ex) = dpres(iip)
399 cs2_(fp,ex) = gamm * ( pres_hyd(iip) + dpres_(fp,ex) ) * gsqrt_(fp,ex) / gsqrtdens(fp,ex)
400 end do
401 do fp=1, elem%NfpTot
402 iip = ip(fp)
403 gsqrtmomx(fp,ex) = gsqrt_(fp,ex) * momx_(iip)
404 gsqrtmomy(fp,ex) = gsqrt_(fp,ex) * momy_(iip)
405 gsqrtmomz(fp,ex) = gsqrt_(fp,ex) * momz_(iip)
406
407 g13_(fp,ex) = g13(iip)
408 g23_(fp,ex) = g23(iip)
409 end do
410 do fp=1, elem%NfpTot
411 fpp = nfp_s + fp
412 vel(fp,ex) = ( gsqrtmomx(fp,ex) * nx(fpp,ke) + gsqrtmomy(fp,ex) * ny(fpp,ke) &
413 + ( ( gsqrtmomz(fp,ex) * rgsqrtv(fp,ex) &
414 + g13_(fp,ex) * gsqrtmomx(fp,ex) + g23_(fp,ex) * gsqrtmomy(fp,ex) ) * nz(fpp,ke) ) &
415 ) / gsqrtdens(fp,ex)
416 end do
417
418 do fp=1, elem%NfpTot
419 tmp1 = abs( nx(fp,ke) ) + abs( ny(fp,ke) )
420 gnn_m = tmp1 &
421 + ( 1.0_rp * rgsqrtv(fp,in)**2 + g13_(fp,in)**2 + g23_(fp,in)**2 ) * abs( nz(fp,ke) )
422 gnn_p = tmp1 &
423 + ( 1.0_rp * rgsqrtv(fp,ex)**2 + g13_(fp,ex)**2 + g23_(fp,ex)**2 ) * abs( nz(fp,ke) )
424
425 alpha(fp) = max( sqrt( gnn_m * cs2_(fp,in) ) + abs(vel(fp,in)), &
426 sqrt( gnn_p * cs2_(fp,ex) ) + abs(vel(fp,ex)) )
427 end do
428 do fp=1, elem%NfpTot
429 tmp1 = lmesh%Fscale(fp,ke) * 0.5_rp
430
431 tmp2 = - alpha(fp) * ( gsqrtddens(fp,ex) - gsqrtddens(fp,in) )
432 del_flux(fp,dens_vid) = tmp1 * ( &
433 gsqrtdens(fp,ex) * vel(fp,ex) - gsqrtdens(fp,in) * vel(fp,in) &
434 + tmp2 )
435
436 tmp2 = - alpha(fp) * ( gsqrtdrhot(fp,ex) - gsqrtdrhot(fp,in) )
437 del_flux(fp,rhot_vid) = tmp1 * ( &
438 gsqrtrhot(fp,ex) * vel(fp,ex) - gsqrtrhot(fp,in) * vel(fp,in) &
439 + tmp2 )
440 end do
441
442 do fp=1, elem%NfpTot
443 tmp3 = gsqrt_(fp,ex) * dpres_(fp,ex)
444 tmp4 = gsqrt_(fp,in) * dpres_(fp,in)
445
446 del_flux_tmp_mom(fp,1) = &
447 ( tmp3 * rgsqrtv(fp,ex) &
448 - tmp4 * rgsqrtv(fp,in) ) * nz(fp,ke)
449
450 del_flux_tmp_mom(fp,2) = &
451 ( nx(fp,ke) + g13_(fp,ex) * nz(fp,ke) ) * tmp3 &
452 - ( nx(fp,ke) + g13_(fp,in) * nz(fp,ke) ) * tmp4
453
454 del_flux_tmp_mom(fp,3) = &
455 ( ny(fp,ke) + g23_(fp,ex) * nz(fp,ke) ) * tmp3 &
456 - ( ny(fp,ke) + g23_(fp,in) * nz(fp,ke) ) * tmp4
457 end do
458
459 do fp=1, elem%NfpTot
460 tmp1 = lmesh%Fscale(fp,ke) * 0.5_rp
461
462 tmp2 = - alpha(fp) * ( gsqrtmomz(fp,ex) - gsqrtmomz(fp,in) )
463 del_flux(fp,momz_vid) = tmp1 * ( &
464 gsqrtmomz(fp,ex) * vel(fp,ex) - gsqrtmomz(fp,in) * vel(fp,in) &
465 + del_flux_tmp_mom(fp,1) &
466 + tmp2 )
467 end do
468
469 !-
470 do fp=1, elem%NfpTot
471 tmp1 = lmesh%Fscale(fp,ke) * 0.5_rp
472
473 tmp2 = - alpha(fp) * ( gsqrtmomx(fp,ex) - gsqrtmomx(fp,in) )
474 del_flux(fp,momx_vid) = tmp1 * ( &
475 gsqrtmomx(fp,ex) * vel(fp,ex) - gsqrtmomx(fp,in) * vel(fp,in) &
476 + del_flux_tmp_mom(fp,2) &
477 + tmp2 )
478
479 tmp2 = - alpha(fp) * ( gsqrtmomy(fp,ex) - gsqrtmomy(fp,in) )
480 del_flux(fp,momy_vid) = tmp1 * ( &
481 gsqrtmomy(fp,ex) * vel(fp,ex) - gsqrtmomy(fp,in) * vel(fp,in) &
482 + del_flux_tmp_mom(fp,3) &
483 + tmp2 )
484 end do
485
486 call elem3d_optr%Lift_var5( del_flux, lift )
487 rgsqrt(:) = 1.0_rp / lmesh%Gsqrt(:,ke)
488 do iv=1, 5
489 prgvar_dt(:,ke,iv) = prgvar_dt(:,ke,iv) - rgsqrt(:) * lift(:,iv)
490 end do
491 end do
492 !$omp end do
493 !$omp end parallel
494 call prof_rapend('cal_dyn_tend_bndflux1', 3)
495 return
497
498
499! !OCL SERIAL
500! subroutine atm_dyn_dgm_nonhydro3d_rhot_heve_add_bnd_contrib2_generalvc( &
501! PRGVAR_dt, & ! (inout)
502! DDENS_, MOMX_, MOMY_, MOMZ_, DRHOT_, DPRES, & ! (in)
503! DENS_hyd, PRES_hyd, THERM_hyd, Rtot, CVtot, CPtot, & ! (in)
504! Gsqrt, G13, G23, nx, ny, nz, & ! (in)
505! vmapM, vmapP, lmesh, elem, lmesh2D, elem2D, elem3D_optr ) ! (in)
506! use scale_element_operation_base, only: ElementOperationBase3D
507! implicit none
508
509! class(LocalMesh3D), intent(in) :: lmesh
510! class(ElementBase3D), intent(in) :: elem
511! class(LocalMesh2D), intent(in) :: lmesh2D
512! class(ElementBase2D), intent(in) :: elem2D
513! real(RP), intent(inout) :: PRGVAR_dt(elem%Np,lmesh%NeA,5)
514! real(RP), intent(in) :: DDENS_(elem%Np*lmesh%NeA)
515! real(RP), intent(in) :: MOMX_(elem%Np*lmesh%NeA)
516! real(RP), intent(in) :: MOMY_(elem%Np*lmesh%NeA)
517! real(RP), intent(in) :: MOMZ_(elem%Np*lmesh%NeA)
518! real(RP), intent(in) :: DRHOT_(elem%Np*lmesh%NeA)
519! real(RP), intent(inout) :: DPRES(elem%Np*lmesh%NeA)
520! real(RP), intent(in) :: DENS_hyd(elem%Np*lmesh%NeA)
521! real(RP), intent(in) :: PRES_hyd(elem%Np*lmesh%NeA)
522! real(RP), intent(in) :: THERM_hyd(elem%Np*lmesh%NeA)
523! real(RP), intent(in) :: Rtot (elem%Np*lmesh%NeA)
524! real(RP), intent(in) :: CVtot(elem%Np*lmesh%NeA)
525! real(RP), intent(in) :: CPtot(elem%Np*lmesh%NeA)
526! real(RP), intent(in) :: Gsqrt(elem%Np*lmesh%NeA)
527! real(RP), intent(in) :: G13(elem%Np*lmesh%NeA)
528! real(RP), intent(in) :: G23(elem%Np*lmesh%NeA)
529! real(RP), intent(in) :: nx(elem%NfpTot,lmesh%Ne)
530! real(RP), intent(in) :: ny(elem%NfpTot,lmesh%Ne)
531! real(RP), intent(in) :: nz(elem%NfpTot,lmesh%Ne)
532! integer, intent(in) :: vmapM(elem%NfpTot,lmesh%Ne)
533! integer, intent(in) :: vmapP(elem%NfpTot,lmesh%Ne)
534! class(ElementOperationBase3D), intent(in) :: elem3D_optr
535
536! integer :: ke, fp, i, iP(elem%NfpTot), iM(elem%NfpTot)
537
538! integer :: iv
539! real(RP) :: del_flux(elem%NfpTot,5)
540! real(RP) :: del_flux_P(elem%NfpTot,5)
541! real(RP) :: del_flux_save_h(elem%Nfp_h,5,elem%Nfaces_h,lmesh%Ne)
542! real(RP) :: del_flux_save_v(elem%Nfp_v,5,elem%Nfaces_v,lmesh%Ne)
543! real(RP) :: Lift(elem%Np,5)
544! real(RP) :: RGsqrt(elem%Np)
545! real(RP) :: tmp1
546
547! integer :: color(lmesh%Ne)
548! integer :: j, k
549
550! integer :: ke2, f, f2, fpp
551! integer :: Nfp_s
552
553! real(RP) :: rP0
554! !------------------------------------------------------------------------
555
556! do k=1, lmesh%NeZ
557! do j=1, lmesh%NeY
558! do i=1, lmesh%NeX
559! ke = i + (j-1)*lmesh%NeX + (k-1)*lmesh%Ne2D
560! if (mod(i+j+k,2) == 0) then
561! color(ke) = 0
562! else
563! color(ke) = 1
564! end if
565! end do
566! end do
567! end do
568! call PROF_rapstart('cal_dyn_tend_bndflux1', 3)
569
570! rP0 = 1.0_RP / PRES00
571
572! !$omp parallel private( &
573! !$omp ke, iM, iP, fp, i, f,fpp,f2,ke2, Nfp_s, &
574! !$omp del_flux, del_flux_P, iv, Lift, RGsqrt, tmp1 )
575
576! !$omp do
577! do i=elem%Np*lmesh%NeE+1, elem%Np*lmesh%NeE+size(lmesh%VMapB)
578! DPRES(i) = PRES00 * ( Rtot(i) * rP0 * ( THERM_hyd(i) + DRHOT_(i) ) )**( CPtot(i) / CVtot(i) ) &
579! - PRES_hyd(i)
580! end do
581! !$omp end do
582
583! !$omp do
584! do ke=lmesh%NeS, lmesh%NeE
585! if (color(ke) == 1) cycle
586
587! iM(:) = vmapM(:,ke); iP(:) = vmapP(:,ke)
588! Nfp_s = 0
589! call calc_bnd_flux_core( del_flux, del_flux_P, & ! (out)
590! DDENS_, MOMX_, MOMY_, MOMZ_, DRHOT_, DPRES, & ! (in)
591! DENS_hyd, PRES_hyd, THERM_hyd, Rtot, CVtot, CPtot, & ! (in)
592! Gsqrt, G13, G23, & ! (in)
593! nx(:,ke), ny(:,ke), nz(:,ke), lmesh%Fscale(:,ke), & ! (in)
594! iM, iP, elem%NfpTot, Nfp_s, elem%NfpTot, lmesh, elem ) ! (in)
595
596! !--
597! do f=1, elem%Nfaces_h
598! ke2 = lmesh%EToE(ke,f); f2 = lmesh%EToF(ke,f)
599! Nfp_s = (f-1)*elem%Nfp_h
600! do fpp=1, elem%Nfp_h
601! fp = Nfp_s + fpp
602! del_flux_save_h(fpp,DENS_VID,f2,ke2) = del_flux_P(fp,DENS_VID)
603! del_flux_save_h(fpp,RHOT_VID,f2,ke2) = del_flux_P(fp,RHOT_VID)
604! del_flux_save_h(fpp,MOMZ_VID,f2,ke2) = del_flux_P(fp,MOMZ_VID)
605! del_flux_save_h(fpp,MOMX_VID,f2,ke2) = del_flux_P(fp,MOMX_VID)
606! del_flux_save_h(fpp,MOMY_VID,f2,ke2) = del_flux_P(fp,MOMY_VID)
607! end do
608
609! ! call calc_bnd_flux_core( del_flux, del_flux_save_h(:,:,f2,ke2), & ! (out)
610! ! DDENS_, MOMX_, MOMY_, MOMZ_, DRHOT_, DPRES, & ! (in)
611! ! DENS_hyd, PRES_hyd, THERM_hyd, Rtot, CVtot, CPtot, & ! (in)
612! ! Gsqrt, G13, G23, & ! (in)
613! ! nx(:,ke), ny(:,ke), nz(:,ke), lmesh%Fscale(:,ke), & ! (in)
614! ! iM(Nfp_s+1:Nfp_s+elem%Nfp_h), iP(Nfp_s+1:Nfp_s+elem%Nfp_h), & ! (in)
615! ! elem%Nfp_h, Nfp_s, elem%NfpTot, lmesh, elem ) ! (in)
616! end do
617! do f=1, elem%Nfaces_v
618! ke2 = lmesh%EToE(ke,f+elem%Nfaces_h); f2 = lmesh%EToF(ke,f+elem%Nfaces_h) - elem%Nfaces_h
619! Nfp_s = elem%Nfaces_h*elem%Nfp_h + (f-1)*elem%Nfp_v
620! do fpp=1, elem%Nfp_v
621! fp = Nfp_s + fpp
622! del_flux_save_v(fpp,DENS_VID,f2,ke2) = del_flux_P(fp,DENS_VID)
623! del_flux_save_v(fpp,RHOT_VID,f2,ke2) = del_flux_P(fp,RHOT_VID)
624! del_flux_save_v(fpp,MOMZ_VID,f2,ke2) = del_flux_P(fp,MOMZ_VID)
625! del_flux_save_v(fpp,MOMX_VID,f2,ke2) = del_flux_P(fp,MOMX_VID)
626! del_flux_save_v(fpp,MOMY_VID,f2,ke2) = del_flux_P(fp,MOMY_VID)
627! end do
628
629! ! call calc_bnd_flux_core( del_flux, del_flux_save_v(:,:,f2,ke2), & ! (out)
630! ! DDENS_, MOMX_, MOMY_, MOMZ_, DRHOT_, DPRES, & ! (in)
631! ! DENS_hyd, PRES_hyd, THERM_hyd, Rtot, CVtot, CPtot, & ! (in)
632! ! Gsqrt, G13, G23, & ! (in)
633! ! nx(:,ke), ny(:,ke), nz(:,ke), lmesh%Fscale(:,ke), & ! (in)
634! ! iM(Nfp_s+1:Nfp_s+elem%Nfp_v), iP(Nfp_s+1:Nfp_s+elem%Nfp_v), & ! (in)
635! ! elem%Nfp_v, Nfp_s, elem%NfpTot, lmesh, elem ) ! (in)
636! end do
637
638! !--
639! call elem3D_optr%Lift_var5( del_flux, Lift )
640! RGsqrt(:) = 1.0_RP / lmesh%Gsqrt(:,ke)
641! do iv=1, 5
642! PRGVAR_dt(:,ke,iv) = PRGVAR_dt(:,ke,iv) - RGsqrt(:) * Lift(:,iv)
643! end do
644! end do
645! !$omp end do
646
647! !$omp do
648! do ke=lmesh%NeS, lmesh%NeE
649! if (color(ke) == 0) cycle
650
651! iM(:) = vmapM(:,ke); iP(:) = vmapP(:,ke)
652! !--
653! do f=1, elem%Nfaces_h
654! ke2 = lmesh%EToE(ke,f); f2 = lmesh%EToF(ke,f)
655! if (ke2 == ke) then
656! Nfp_s = (f-1)*elem%Nfp_h
657! call calc_bnd_flux_core( del_flux, del_flux_save_h(:,:,f2,ke2), & ! (out)
658! DDENS_, MOMX_, MOMY_, MOMZ_, DRHOT_, DPRES, & ! (in)
659! DENS_hyd, PRES_hyd, THERM_hyd, Rtot, CVtot, CPtot, & ! (in)
660! Gsqrt, G13, G23, & ! (in)
661! nx(:,ke), ny(:,ke), nz(:,ke), lmesh%Fscale(:,ke), & ! (in)
662! iM(Nfp_s+1:Nfp_s+elem%Nfp_h), iP(Nfp_s+1:Nfp_s+elem%Nfp_h), & ! (in)
663! elem%Nfp_h, Nfp_s, elem%NfpTot, lmesh, elem ) ! (in)
664! else
665! do fpp=1, elem%Nfp_h
666! fp = fpp + (f-1)*elem%Nfp_h
667! tmp1 = 0.5_RP * lmesh%Fscale(fp,ke)
668! del_flux(fp,DENS_VID) = tmp1 * del_flux_save_h(fpp,DENS_VID,f,ke)
669! del_flux(fp,RHOT_VID) = tmp1 * del_flux_save_h(fpp,RHOT_VID,f,ke)
670! del_flux(fp,MOMZ_VID) = tmp1 * del_flux_save_h(fpp,MOMZ_VID,f,ke)
671! del_flux(fp,MOMX_VID) = tmp1 * del_flux_save_h(fpp,MOMX_VID,f,ke)
672! del_flux(fp,MOMY_VID) = tmp1 * del_flux_save_h(fpp,MOMY_VID,f,ke)
673! end do
674! end if
675! end do
676! do f=1, elem%Nfaces_v
677! ke2 = lmesh%EToE(ke,f+elem%Nfaces_h); f2 = lmesh%EToF(ke,f+elem%Nfaces_h) - elem%Nfaces_h
678! if (ke2 == ke) then
679! Nfp_s = elem%Nfaces_h*elem%Nfp_h + (f-1)*elem%Nfp_v
680! call calc_bnd_flux_core( del_flux, del_flux_save_v(:,:,f2,ke2), & ! (out)
681! DDENS_, MOMX_, MOMY_, MOMZ_, DRHOT_, DPRES, & ! (in)
682! DENS_hyd, PRES_hyd, THERM_hyd, Rtot, CVtot, CPtot, & ! (in)
683! Gsqrt, G13, G23, & ! (in)
684! nx(:,ke), ny(:,ke), nz(:,ke), lmesh%Fscale(:,ke), & ! (in)
685! iM(Nfp_s+1:Nfp_s+elem%Nfp_v), iP(Nfp_s+1:Nfp_s+elem%Nfp_v), & ! (in)
686! elem%Nfp_v, Nfp_s, elem%NfpTot, lmesh, elem ) ! (in)
687! else
688! do fpp=1, elem%Nfp_v
689! fp = fpp + (f-1)*elem%Nfp_v + elem%Nfaces_h*elem%Nfp_h
690! tmp1 = 0.5_RP * lmesh%Fscale(fp,ke)
691! del_flux(fp,DENS_VID) = tmp1 * del_flux_save_v(fpp,DENS_VID,f,ke)
692! del_flux(fp,RHOT_VID) = tmp1 * del_flux_save_v(fpp,RHOT_VID,f,ke)
693! del_flux(fp,MOMZ_VID) = tmp1 * del_flux_save_v(fpp,MOMZ_VID,f,ke)
694! del_flux(fp,MOMX_VID) = tmp1 * del_flux_save_v(fpp,MOMX_VID,f,ke)
695! del_flux(fp,MOMY_VID) = tmp1 * del_flux_save_v(fpp,MOMY_VID,f,ke)
696! end do
697! end if
698! end do
699
700! call elem3D_optr%Lift_var5( del_flux, Lift )
701! RGsqrt(:) = 1.0_RP / lmesh%Gsqrt(:,ke)
702! do iv=1, 5
703! PRGVAR_dt(:,ke,iv) = PRGVAR_dt(:,ke,iv) - RGsqrt(:) * Lift(:,iv)
704! end do
705! end do
706! !$omp end do
707! !$omp end parallel
708
709! call PROF_rapend('cal_dyn_tend_bndflux1', 3)
710! return
711! end subroutine atm_dyn_dgm_nonhydro3d_rhot_heve_add_bnd_contrib2_generalvc
712
713! !OCL SERIAL
714! subroutine calc_bnd_flux_core( del_flux, del_flux_save, &
715! DDENS_, MOMX_, MOMY_, MOMZ_, DRHOT_, DPRES, & ! (in)
716! DENS_hyd, PRES_hyd, THERM_hyd, Rtot, CVtot, CPtot, & ! (in)
717! Gsqrt, G13, G23, & ! (in)
718! nx, ny, nz, Fscale, iM, iP, Nfp, Nfp_s, NfpTot, lmesh, elem )
719! integer, intent(in) :: Nfp
720! integer, intent(in) :: NfpTot
721! class(LocalMesh3D), intent(in) :: lmesh
722! class(ElementBase3D), intent(in) :: elem
723! real(RP), intent(out) :: del_flux(NfpTot,5)
724! real(RP), intent(out) :: del_flux_save(Nfp,5)
725! real(RP), intent(in) :: DDENS_(elem%Np*lmesh%NeA)
726! real(RP), intent(in) :: MOMX_(elem%Np*lmesh%NeA)
727! real(RP), intent(in) :: MOMY_(elem%Np*lmesh%NeA)
728! real(RP), intent(in) :: MOMZ_(elem%Np*lmesh%NeA)
729! real(RP), intent(in) :: DRHOT_(elem%Np*lmesh%NeA)
730! real(RP), intent(inout) :: DPRES(elem%Np*lmesh%NeA)
731! real(RP), intent(in) :: DENS_hyd(elem%Np*lmesh%NeA)
732! real(RP), intent(in) :: PRES_hyd(elem%Np*lmesh%NeA)
733! real(RP), intent(in) :: THERM_hyd(elem%Np*lmesh%NeA)
734! real(RP), intent(in) :: Rtot (elem%Np*lmesh%NeA)
735! real(RP), intent(in) :: CVtot(elem%Np*lmesh%NeA)
736! real(RP), intent(in) :: CPtot(elem%Np*lmesh%NeA)
737! real(RP), intent(in) :: Gsqrt(elem%Np*lmesh%NeA)
738! real(RP), intent(in) :: G13(elem%Np*lmesh%NeA)
739! real(RP), intent(in) :: G23(elem%Np*lmesh%NeA)
740! real(RP), intent(in) :: nx(NfpTot)
741! real(RP), intent(in) :: ny(NfpTot)
742! real(RP), intent(in) :: nz(NfpTot)
743! real(RP), intent(in) :: Fscale(NfpTot)
744! integer, intent(in) :: iM(Nfp)
745! integer, intent(in) :: iP(Nfp)
746! integer, intent(in) :: Nfp_s
747! !-
748! real(RP) :: Vel(Nfp,2), alpha(Nfp)
749! real(RP) :: DPRES_(Nfp,2)
750! real(RP) :: GsqrtDens(Nfp,2)
751! real(RP) :: GsqrtRhot(Nfp,2)
752! real(RP) :: GsqrtDDENS(Nfp,2)
753! real(RP) :: GsqrtMOMX(Nfp,2)
754! real(RP) :: GsqrtMOMY(Nfp,2)
755! real(RP) :: GsqrtMOMZ(Nfp,2)
756! real(RP) :: GsqrtDRHOT(Nfp,2)
757! real(RP) :: Cs2_(Nfp,2)
758! real(RP) :: Gsqrt_(Nfp,2)
759! real(RP) :: GsqrtV_(Nfp,2)
760! real(RP) :: RGsqrtV(Nfp,2)
761! real(RP) :: G13_(Nfp,2)
762! real(RP) :: G23_(Nfp,2)
763! implicit none
764
765! real(RP) :: Gnn_M, Gnn_P
766
767! real(RP) :: gamm
768
769! integer, parameter :: IN = 1
770! integer, parameter :: EX = 2
771
772! real(RP) :: tmp1, tmp2, tmp3, tmp4, tmp01, tmp02, tmp03
773! real(RP) :: del_flux_tmp_mom(Nfp,3)
774
775! integer :: fp, fpp
776! integer :: iv
777
778! integer :: iim, iip
779! !------------------------
780
781! gamm = CPDry / CvDry
782
783! !-
784! !OCL PREFETCH
785! do fp=1, Nfp
786! iim = iM(fp)
787! Gsqrt_(fp,IN) = Gsqrt(iim)
788! GsqrtV_(fp,IN) = Gsqrt_(fp,IN)
789! RGsqrtV(fp,IN) = 1.0_RP / GsqrtV_(fp,IN)
790
791! GsqrtDDENS(fp,IN) = Gsqrt_(fp,IN) * DDENS_(iim)
792! GsqrtDens(fp,IN) = GsqrtDDENS(fp,IN) + Gsqrt_(fp,IN) * DENS_hyd(iim)
793
794! GsqrtDRHOT(fp,IN) = Gsqrt_(fp,IN) * DRHOT_(iim)
795! GsqrtRhot(fp,IN) = GsqrtDRHOT(fp,IN) + Gsqrt_(fp,IN) * THERM_hyd(iim)
796! end do
797! !OCL PREFETCH
798! do fp=1, Nfp
799! iim = iM(fp)
800
801! DPRES_(fp,IN) = DPRES(iim)
802! Cs2_(fp,IN) = gamm * ( PRES_hyd(iim) + DPRES_(fp,IN) ) * Gsqrt_(fp,IN) / GsqrtDens(fp,IN)
803! end do
804! !OCL PREFETCH
805! do fp=1, Nfp
806! iim = iM(fp)
807! GsqrtMOMX (fp,IN) = Gsqrt_(fp,IN) * MOMX_ (iim)
808! GsqrtMOMY (fp,IN) = Gsqrt_(fp,IN) * MOMY_ (iim)
809! GsqrtMOMZ (fp,IN) = Gsqrt_(fp,IN) * MOMZ_ (iim)
810
811! G13_(fp,IN) = G13(iim)
812! G23_(fp,IN) = G23(iim)
813! end do
814! do fp=1, Nfp
815! fpp = Nfp_s + fp
816! Vel(fp,IN) = ( GsqrtMOMX(fp,IN) * nx(fpp) + GsqrtMOMY(fp,IN) * ny(fpp) &
817! + ( ( GsqrtMOMZ(fp,IN) * RGsqrtV(fp,IN) &
818! + G13_(fp,IN) * GsqrtMOMX(fp,IN) + G23_(fp,IN) * GsqrtMOMY(fp,IN) ) * nz(fpp) ) &
819! ) / GsqrtDens(fp,IN)
820! end do
821
822! !-
823! !OCL PREFETCH
824! do fp=1, Nfp
825! iip = iP(fp)
826
827! Gsqrt_(fp,EX) = Gsqrt(iip)
828! GsqrtV_(fp,EX) = Gsqrt_(fp,EX)
829! RGsqrtV(fp,EX) = 1.0_RP / GsqrtV_(fp,EX)
830
831
832! GsqrtDDENS(fp,EX) = Gsqrt_(fp,EX) * DDENS_(iip)
833! GsqrtDens(fp,EX) = GsqrtDDENS(fp,EX) + Gsqrt_(fp,EX) * DENS_hyd(iip)
834
835! GsqrtDRHOT(fp,EX) = Gsqrt_(fp,EX) * DRHOT_(iip)
836! GsqrtRhot(fp,EX) = GsqrtDRHOT(fp,EX) + Gsqrt_(fp,EX) * THERM_hyd(iip)
837! end do
838! !OCL PREFETCH
839! do fp=1, Nfp
840! iip = iP(fp)
841! DPRES_(fp,EX) = DPRES(iip)
842! Cs2_(fp,EX) = gamm * ( PRES_hyd(iip) + DPRES_(fp,EX) ) * Gsqrt_(fp,EX) / GsqrtDens(fp,EX)
843! end do
844! !OCL PREFETCH
845! do fp=1, Nfp
846! iip = iP(fp)
847! GsqrtMOMX (fp,EX) = Gsqrt_(fp,EX) * MOMX_ (iip)
848! GsqrtMOMY (fp,EX) = Gsqrt_(fp,EX) * MOMY_ (iip)
849! GsqrtMOMZ (fp,EX) = Gsqrt_(fp,EX) * MOMZ_ (iip)
850
851! G13_(fp,EX) = G13(iip)
852! G23_(fp,EX) = G23(iip)
853! end do
854! do fp=1, Nfp
855! fpp = Nfp_s + fp
856! Vel(fp,EX) = ( GsqrtMOMX(fp,EX) * nx(fpp) + GsqrtMOMY(fp,EX) * ny(fpp) &
857! + ( ( GsqrtMOMZ(fp,EX) * RGsqrtV(fp,EX) &
858! + G13_(fp,EX) * GsqrtMOMX(fp,EX) + G23_(fp,EX) * GsqrtMOMY(fp,EX) ) * nz(fpp) ) &
859! ) / GsqrtDens(fp,EX)
860! end do
861
862! !-
863! do fp=1, Nfp
864! fpp = Nfp_s + fp
865! tmp1 = abs( nx(fpp) ) + abs( ny(fpp) )
866! Gnn_M = tmp1 &
867! + ( 1.0_RP * RGsqrtV(fp,IN)**2 + G13_(fp,IN)**2 + G23_(fp,IN)**2 ) * abs( nz(fpp) )
868! Gnn_P = tmp1 &
869! + ( 1.0_RP * RGsqrtV(fp,EX)**2 + G13_(fp,EX)**2 + G23_(fp,EX)**2 ) * abs( nz(fpp) )
870
871! alpha(fp) = max( sqrt( Gnn_M * Cs2_(fp,IN) ) + abs(Vel(fp,IN)), &
872! sqrt( Gnn_P * Cs2_(fp,EX) ) + abs(Vel(fp,EX)) )
873! end do
874! do fp=1, Nfp
875! fpp = Nfp_s + fp
876! tmp1 = Fscale(fpp) * 0.5_RP
877
878! !-
879! tmp2 = GsqrtDens(fp,EX) * Vel(fp,EX) - GsqrtDens(fp,IN) * Vel(fp,IN)
880! tmp3 = - alpha(fp) * ( GsqrtDDENS(fp,EX) - GsqrtDDENS(fp,IN) )
881
882! del_flux(fpp,DENS_VID) = tmp1 * ( tmp2 + tmp3 )
883! del_flux_save(fp,DENS_VID) = tmp2 - tmp3
884
885! !-
886! tmp2 = GsqrtRhot(fp,EX) * Vel(fp,EX) - GsqrtRhot(fp,IN) * Vel(fp,IN)
887! tmp3 = - alpha(fp) * ( GsqrtDRHOT(fp,EX) - GsqrtDRHOT(fp,IN) )
888
889! del_flux(fpp,RHOT_VID) = tmp1 * ( tmp2 + tmp3 )
890! del_flux_save(fp,RHOT_VID) = tmp2 - tmp3
891! end do
892
893! do fp=1, Nfp
894! fpp = Nfp_s + fp
895! tmp3 = Gsqrt_(fp,EX) * dpres_(fp,EX)
896! tmp4 = Gsqrt_(fp,IN) * dpres_(fp,IN)
897
898! del_flux_tmp_mom(fp,1) = &
899! ( tmp3 * RGsqrtV(fp,EX) &
900! - tmp4 * RGsqrtV(fp,IN) ) * nz(fpp)
901
902! del_flux_tmp_mom(fp,2) = &
903! ( nx(fpp) + G13_(fp,EX) * nz(fpp) ) * tmp3 &
904! - ( nx(fpp) + G13_(fp,IN) * nz(fpp) ) * tmp4
905
906! del_flux_tmp_mom(fp,3) = &
907! ( ny(fpp) + G23_(fp,EX) * nz(fpp) ) * tmp3 &
908! - ( ny(fpp) + G23_(fp,IN) * nz(fpp) ) * tmp4
909! end do
910
911! do fp=1, Nfp
912! fpp = Nfp_s + fp
913! tmp1 = Fscale(fpp) * 0.5_RP
914
915! tmp2 = GsqrtMOMZ(fp,EX) * Vel(fp,EX) - GsqrtMOMZ(fp,IN) * Vel(fp,IN)
916! tmp3 = - alpha(fp) * ( GsqrtMOMZ(fp,EX) - GsqrtMOMZ(fp,IN) )
917
918! del_flux(fpp,MOMZ_VID) = tmp1 * ( tmp2 + del_flux_tmp_mom(fp,1) + tmp3 )
919! del_flux_save(fp,MOMZ_VID) = tmp2 + del_flux_tmp_mom(fp,1) - tmp3
920! end do
921
922! !-
923! do fp=1, Nfp
924! fpp = Nfp_s + fp
925! tmp1 = Fscale(fpp) * 0.5_RP
926
927! !-
928! tmp2 = GsqrtMOMX(fp,EX) * Vel(fp,EX) - GsqrtMOMX(fp,IN) * Vel(fp,IN)
929! tmp3 = - alpha(fp) * ( GsqrtMOMX(fp,EX) - GsqrtMOMX(fp,IN) )
930
931! del_flux(fpp,MOMX_VID) = tmp1 * ( tmp2 + del_flux_tmp_mom(fp,2) + tmp3 )
932! del_flux_save(fp,MOMX_VID) = tmp2 + del_flux_tmp_mom(fp,2) - tmp3
933
934! !-
935! tmp2 = GsqrtMOMY(fp,EX) * Vel(fp,EX) - GsqrtMOMY(fp,IN) * Vel(fp,IN)
936! tmp3 = - alpha(fp) * ( GsqrtMOMY(fp,EX) - GsqrtMOMY(fp,IN) )
937
938! del_flux(fpp,MOMY_VID) = tmp1 * ( tmp2 + del_flux_tmp_mom(fp,3) + tmp3 )
939! del_flux_save(fp,MOMY_VID) = tmp2 + del_flux_tmp_mom(fp,3) - tmp3
940! end do
941
942! return
943! end subroutine calc_bnd_flux_core
944
945!OCL SERIAL
947 del_flux, & ! (out)
948 ddens_, momx_, momy_, momz_, drhot_, dpres, & ! (in)
949 dens_hyd, pres_hyd, therm_hyd, rtot, cvtot, cptot, & ! (in)
950 gsqrt, g13, g23, nx, ny, nz, & ! (in)
951 vmapm, vmapp, lmesh, elem, lmesh2d, elem2d ) ! (in)
952
953 implicit none
954
955 class(localmesh3d), intent(in) :: lmesh
956 class(elementbase3d), intent(in) :: elem
957 class(localmesh2d), intent(in) :: lmesh2d
958 class(elementbase2d), intent(in) :: elem2d
959 real(rp), intent(out) :: del_flux(elem%nfptot,prgvar_num,lmesh%ne)
960 real(rp), intent(in) :: ddens_(elem%np*lmesh%nea)
961 real(rp), intent(in) :: momx_(elem%np*lmesh%nea)
962 real(rp), intent(in) :: momy_(elem%np*lmesh%nea)
963 real(rp), intent(in) :: momz_(elem%np*lmesh%nea)
964 real(rp), intent(in) :: drhot_(elem%np*lmesh%nea)
965 real(rp), intent(in) :: dpres(elem%np*lmesh%nea)
966 real(rp), intent(in) :: dens_hyd(elem%np*lmesh%nea)
967 real(rp), intent(in) :: pres_hyd(elem%np*lmesh%nea)
968 real(rp), intent(in) :: therm_hyd(elem%np*lmesh%nea)
969 real(rp), intent(in) :: rtot (elem%np*lmesh%nea)
970 real(rp), intent(in) :: cvtot(elem%np*lmesh%nea)
971 real(rp), intent(in) :: cptot(elem%np*lmesh%nea)
972 real(rp), intent(in) :: gsqrt(elem%np*lmesh%nea)
973 real(rp), intent(in) :: g13(elem%np*lmesh%nea)
974 real(rp), intent(in) :: g23(elem%np*lmesh%nea)
975 real(rp), intent(in) :: nx(elem%nfptot,lmesh%ne)
976 real(rp), intent(in) :: ny(elem%nfptot,lmesh%ne)
977 real(rp), intent(in) :: nz(elem%nfptot,lmesh%ne)
978 integer, intent(in) :: vmapm(elem%nfptot,lmesh%ne)
979 integer, intent(in) :: vmapp(elem%nfptot,lmesh%ne)
980
981 integer :: ke, fp, i, ip(elem%nfptot), im(elem%nfptot)
982 integer :: ke2d
983 real(rp) :: vel(elem%nfptot,2), alpha(elem%nfptot)
984 real(rp) :: dpres_(elem%nfptot,2)
985 real(rp) :: gsqrtdens(elem%nfptot,2)
986 real(rp) :: gsqrtrhot(elem%nfptot,2)
987 real(rp) :: gsqrtddens(elem%nfptot,2)
988 real(rp) :: gsqrtmomx(elem%nfptot,2)
989 real(rp) :: gsqrtmomy(elem%nfptot,2)
990 real(rp) :: gsqrtmomz(elem%nfptot,2)
991 real(rp) :: gsqrtdrhot(elem%nfptot,2)
992 real(rp) :: phyd_(elem%nfptot,2)
993 real(rp) :: gsqrt_(elem%nfptot,2)
994 real(rp) :: gsqrtv_(elem%nfptot,2)
995 real(rp) :: rgsqrtv(elem%nfptot,2)
996 real(rp) :: g13_(elem%nfptot,2)
997 real(rp) :: g23_(elem%nfptot,2)
998 real(rp) :: gnn_m, gnn_p
999
1000 real(rp) :: gamm, rgamm
1001 real(rp) :: rp0
1002 real(rp) :: rovp0, p0ovr
1003
1004 integer, parameter :: in = 1
1005 integer, parameter :: ex = 2
1006
1007 real(rp) :: tmp1, tmp2, tmp3, tmp4, tmp01, tmp02, tmp03
1008 real(rp) :: del_flux_tmp_mom(elem%nfptot,3)
1009 !------------------------------------------------------------------------
1010
1011 gamm = cpdry / cvdry
1012 rgamm = cvdry / cpdry
1013 rp0 = 1.0_rp / pres00
1014 rovp0 = rdry * rp0
1015 p0ovr = pres00 / rdry
1016
1017 !$omp parallel do private( &
1018 !$omp ke, iM, iP, ke2D, fp, &
1019 !$omp alpha, Vel, &
1020 !$omp dpres_, GsqrtDens, GsqrtRhot, &
1021 !$omp GsqrtDDENS, GsqrtMOMX, GsqrtMOMY, GsqrtMOMZ, GsqrtDRHOT, &
1022 !$omp Phyd_, &
1023 !$omp Gsqrt_, GsqrtV_, RGsqrtV, G13_, G23_, &
1024 !$omp Gnn_P, Gnn_M, tmp1, tmp2, tmp3, tmp4, tmp01, tmp02, tmp03, del_flux_tmp_mom )
1025!OCL PREFETCH
1026 do ke=lmesh%NeS, lmesh%NeE
1027 im(:) = vmapm(:,ke); ip(:) = vmapp(:,ke)
1028 ke2d = lmesh%EMap3Dto2D(ke)
1029
1030 gsqrt_(:,in) = gsqrt(im)
1031 gsqrt_(:,ex) = gsqrt(ip)
1032 gsqrtv_(:,:) = gsqrt_(:,:)
1033 rgsqrtv(:,:) = 1.0_rp / gsqrtv_(:,:)
1034
1035 g13_(:,in) = g13(im)
1036 g13_(:,ex) = g13(ip)
1037 g23_(:,in) = g23(im)
1038 g23_(:,ex) = g23(ip)
1039
1040 gsqrtddens(:,in) = gsqrt_(:,in) * ddens_(im)
1041 gsqrtddens(:,ex) = gsqrt_(:,ex) * ddens_(ip)
1042 gsqrtmomx(:,in) = gsqrt_(:,in) * momx_(im)
1043 gsqrtmomx(:,ex) = gsqrt_(:,ex) * momx_(ip)
1044 gsqrtmomy(:,in) = gsqrt_(:,in) * momy_(im)
1045 gsqrtmomy(:,ex) = gsqrt_(:,ex) * momy_(ip)
1046 gsqrtmomz(:,in) = gsqrt_(:,in) * momz_(im)
1047 gsqrtmomz(:,ex) = gsqrt_(:,ex) * momz_(ip)
1048 gsqrtdrhot(:,in) = gsqrt_(:,in) * drhot_(im)
1049 gsqrtdrhot(:,ex) = gsqrt_(:,ex) * drhot_(ip)
1050
1051 phyd_(:,in) = pres_hyd(im)
1052 phyd_(:,ex) = pres_hyd(ip)
1053 dpres_(:,in) = dpres(im)
1054 dpres_(:,ex) = dpres(ip)
1055
1056 gsqrtdens(:,in) = gsqrtddens(:,in) + gsqrt_(:,in) * dens_hyd(im)
1057 gsqrtdens(:,ex) = gsqrtddens(:,ex) + gsqrt_(:,ex) * dens_hyd(ip)
1058
1059 gsqrtrhot(:,in) = gsqrt_(:,in) * therm_hyd(im) + gsqrtdrhot(:,in)
1060 gsqrtrhot(:,ex) = gsqrt_(:,ex) * therm_hyd(ip) + gsqrtdrhot(:,ex)
1061
1062 vel(:,in) = ( gsqrtmomx(:,in) * nx(:,ke) + gsqrtmomy(:,in) * ny(:,ke) &
1063 + ( ( gsqrtmomz(:,in) * rgsqrtv(:,in) &
1064 + g13_(:,in) * gsqrtmomx(:,in) + g23_(:,in) * gsqrtmomy(:,in) ) * nz(:,ke) ) &
1065 ) / gsqrtdens(:,in)
1066 vel(:,ex) = ( gsqrtmomx(:,ex) * nx(:,ke) + gsqrtmomy(:,ex) * ny(:,ke) &
1067 + ( ( gsqrtmomz(:,ex) * rgsqrtv(:,ex) &
1068 + g13_(:,ex) * gsqrtmomx(:,ex) + g23_(:,ex) * gsqrtmomy(:,ex) ) * nz(:,ke) ) &
1069 ) / gsqrtdens(:,ex)
1070
1071 do fp=1, elem%NfpTot
1072 tmp1 = abs( nx(fp,ke) ) + abs( ny(fp,ke) )
1073 gnn_m = tmp1 &
1074 + ( 1.0_rp * rgsqrtv(fp,in)**2 + g13_(fp,in)**2 + g23_(fp,in)**2 ) * abs( nz(fp,ke) )
1075 gnn_p = tmp1 &
1076 + ( 1.0_rp * rgsqrtv(fp,ex)**2 + g13_(fp,ex)**2 + g23_(fp,ex)**2 ) * abs( nz(fp,ke) )
1077
1078 alpha(fp) = max( sqrt( gnn_m * gamm * ( phyd_(fp,in) + dpres_(fp,in) ) * gsqrt_(fp,in) / gsqrtdens(fp,in) ) + abs(vel(fp,in)), &
1079 sqrt( gnn_p * gamm * ( phyd_(fp,ex) + dpres_(fp,ex) ) * gsqrt_(fp,ex) / gsqrtdens(fp,ex) ) + abs(vel(fp,ex)) )
1080 end do
1081 do fp=1, elem%NfpTot
1082 tmp1 = lmesh%Fscale(fp,ke) * 0.5_rp
1083
1084 tmp2 = - alpha(fp) * ( gsqrtddens(fp,ex) - gsqrtddens(fp,in) )
1085 del_flux(fp,dens_vid,ke) = tmp1 * ( &
1086 gsqrtdens(fp,ex) * vel(fp,ex) - gsqrtdens(fp,in) * vel(fp,in) &
1087 + tmp2 )
1088
1089 tmp2 = - alpha(fp) * ( gsqrtdrhot(fp,ex) - gsqrtdrhot(fp,in) )
1090 del_flux(fp,rhot_vid,ke) = tmp1 * ( &
1091 gsqrtrhot(fp,ex) * vel(fp,ex) - gsqrtrhot(fp,in) * vel(fp,in) &
1092 + tmp2 )
1093 end do
1094
1095 do fp=1, elem%NfpTot
1096 tmp3 = gsqrt_(fp,ex) * dpres_(fp,ex)
1097 tmp4 = gsqrt_(fp,in) * dpres_(fp,in)
1098
1099 del_flux_tmp_mom(fp,1) = &
1100 ( tmp3 * rgsqrtv(fp,ex) &
1101 - tmp4 * rgsqrtv(fp,in) ) * nz(fp,ke)
1102
1103 del_flux_tmp_mom(fp,2) = &
1104 ( nx(fp,ke) + g13_(fp,ex) * nz(fp,ke) ) * tmp3 &
1105 - ( nx(fp,ke) + g13_(fp,in) * nz(fp,ke) ) * tmp4
1106
1107 del_flux_tmp_mom(fp,3) = &
1108 ( ny(fp,ke) + g23_(fp,ex) * nz(fp,ke) ) * tmp3 &
1109 - ( ny(fp,ke) + g23_(fp,in) * nz(fp,ke) ) * tmp4
1110 end do
1111 do fp=1, elem%NfpTot
1112 tmp1 = lmesh%Fscale(fp,ke) * 0.5_rp
1113
1114 tmp2 = - alpha(fp) * ( gsqrtmomz(fp,ex) - gsqrtmomz(fp,in) )
1115 del_flux(fp,momz_vid,ke) = tmp1 * ( &
1116 gsqrtmomz(fp,ex) * vel(fp,ex) - gsqrtmomz(fp,in) * vel(fp,in) &
1117 + del_flux_tmp_mom(fp,1) &
1118 + tmp2 )
1119 end do
1120 do fp=1, elem%NfpTot
1121 tmp1 = lmesh%Fscale(fp,ke) * 0.5_rp
1122
1123 tmp2 = - alpha(fp) * ( gsqrtmomx(fp,ex) - gsqrtmomx(fp,in) )
1124 del_flux(fp,momx_vid,ke) = tmp1 * ( &
1125 gsqrtmomx(fp,ex) * vel(fp,ex) - gsqrtmomx(fp,in) * vel(fp,in) &
1126 + del_flux_tmp_mom(fp,2) &
1127 + tmp2 )
1128
1129 tmp2 = - alpha(fp) * ( gsqrtmomy(fp,ex) - gsqrtmomy(fp,in) )
1130 del_flux(fp,momy_vid,ke) = tmp1 * ( &
1131 gsqrtmomy(fp,ex) * vel(fp,ex) - gsqrtmomy(fp,in) * vel(fp,in) &
1132 + del_flux_tmp_mom(fp,3) &
1133 + tmp2 )
1134 end do
1135 end do
1136
1137 return
1139
1140!OCL SERIAL
1141 subroutine atm_dyn_dgm_nonhydro3d_rhot_heve_numflux_get_generalvc_lmars( &
1142 del_flux, & ! (out)
1143 ddens_, momx_, momy_, momz_, drhot_, dpres, dens_hyd, pres_hyd, & ! (in)
1144 rtot, cvtot, cptot, & ! (in)
1145 gsqrt, g13, g23, nx, ny, nz, & ! (in)
1146 vmapm, vmapp, lmesh, elem, lmesh2d, elem2d ) ! (in)
1147
1148 implicit none
1149
1150 class(localmesh3d), intent(in) :: lmesh
1151 class(elementbase3d), intent(in) :: elem
1152 class(localmesh2d), intent(in) :: lmesh2d
1153 class(elementbase2d), intent(in) :: elem2d
1154 real(rp), intent(out) :: del_flux(elem%nfptot,prgvar_num,lmesh%ne)
1155 real(rp), intent(in) :: ddens_(elem%np*lmesh%nea)
1156 real(rp), intent(in) :: momx_(elem%np*lmesh%nea)
1157 real(rp), intent(in) :: momy_(elem%np*lmesh%nea)
1158 real(rp), intent(in) :: momz_(elem%np*lmesh%nea)
1159 real(rp), intent(in) :: drhot_(elem%np*lmesh%nea)
1160 real(rp), intent(in) :: dpres(elem%np*lmesh%nea)
1161 real(rp), intent(in) :: dens_hyd(elem%np*lmesh%nea)
1162 real(rp), intent(in) :: pres_hyd(elem%np*lmesh%nea)
1163 real(rp), intent(in) :: rtot (elem%np*lmesh%nea)
1164 real(rp), intent(in) :: cvtot(elem%np*lmesh%nea)
1165 real(rp), intent(in) :: cptot(elem%np*lmesh%nea)
1166 real(rp), intent(in) :: gsqrt(elem%np*lmesh%nea)
1167 real(rp), intent(in) :: g13(elem%np*lmesh%nea)
1168 real(rp), intent(in) :: g23(elem%np*lmesh%nea)
1169 real(rp), intent(in) :: nx(elem%nfptot,lmesh%ne)
1170 real(rp), intent(in) :: ny(elem%nfptot,lmesh%ne)
1171 real(rp), intent(in) :: nz(elem%nfptot,lmesh%ne)
1172 integer, intent(in) :: vmapm(elem%nfptot,lmesh%ne)
1173 integer, intent(in) :: vmapp(elem%nfptot,lmesh%ne)
1174
1175 integer :: ke, fp, i, ip(elem%nfptot), im(elem%nfptot)
1176 integer :: ke2d
1177 real(rp) :: vel(elem%nfptot,2), alpha(elem%nfptot), alpha2(elem%nfptot)
1178 real(rp) :: dpres_(elem%nfptot,2)
1179 real(rp) :: gsqrtdens(elem%nfptot,2)
1180 real(rp) :: gsqrtrhot(elem%nfptot,2)
1181 real(rp) :: gsqrtddens(elem%nfptot,2)
1182 real(rp) :: gsqrtmomx(elem%nfptot,2)
1183 real(rp) :: gsqrtmomy(elem%nfptot,2)
1184 real(rp) :: gsqrtmomz(elem%nfptot,2)
1185 real(rp) :: gsqrtdrhot(elem%nfptot,2)
1186 real(rp) :: phyd_(elem%nfptot,2)
1187 real(rp) :: gsqrt_(elem%nfptot,2)
1188 real(rp) :: gsqrtv_(elem%nfptot,2)
1189 real(rp) :: rgsqrtv(elem%nfptot,2)
1190 real(rp) :: g13_(elem%nfptot,2)
1191 real(rp) :: g23_(elem%nfptot,2)
1192 real(rp) :: gnn_m, gnn_p
1193
1194 real(rp) :: gamm, rgamm
1195 real(rp) :: rp0
1196 real(rp) :: rovp0, p0ovr
1197
1198 integer, parameter :: in = 1
1199 integer, parameter :: ex = 2
1200
1201 real(rp) :: tmp1, tmp2, tmp3, tmp4, tmp01, tmp02, tmp03
1202 real(rp) :: del_flux_tmp_mom(elem%nfptot,3)
1203
1204 real(rp) :: vel_n(elem%nfptot)
1205 real(rp) :: cs(elem%nfptot,2)
1206 !------------------------------------------------------------------------
1207
1208 gamm = cpdry / cvdry
1209 rgamm = cvdry / cpdry
1210 rp0 = 1.0_rp / pres00
1211 rovp0 = rdry * rp0
1212 p0ovr = pres00 / rdry
1213
1214 !$omp parallel do private( &
1215 !$omp ke, iM, iP, ke2D, fp, &
1216 !$omp alpha, alpha2, Vel, Vel_n, Cs, &
1217 !$omp dpres_, GsqrtDens, GsqrtRhot, &
1218 !$omp GsqrtDDENS, GsqrtMOMX, GsqrtMOMY, GsqrtMOMZ, GsqrtDRHOT, &
1219 !$omp Phyd_, &
1220 !$omp Gsqrt_, GsqrtV_, RGsqrtV, G13_, G23_, &
1221 !$omp Gnn_P, Gnn_M, tmp1, tmp2, tmp3, tmp4, tmp01, tmp02, tmp03, del_flux_tmp_mom )
1222!OCL PREFETCH
1223 do ke=lmesh%NeS, lmesh%NeE
1224 im(:) = vmapm(:,ke); ip(:) = vmapp(:,ke)
1225 ke2d = lmesh%EMap3Dto2D(ke)
1226
1227 gsqrt_(:,in) = gsqrt(im)
1228 gsqrt_(:,ex) = gsqrt(ip)
1229 gsqrtv_(:,:) = gsqrt_(:,:)
1230 rgsqrtv(:,:) = 1.0_rp / gsqrtv_(:,:)
1231
1232 g13_(:,in) = g13(im)
1233 g13_(:,ex) = g13(ip)
1234 g23_(:,in) = g23(im)
1235 g23_(:,ex) = g23(ip)
1236
1237 gsqrtddens(:,in) = gsqrt_(:,in) * ddens_(im)
1238 gsqrtddens(:,ex) = gsqrt_(:,ex) * ddens_(ip)
1239 gsqrtmomx(:,in) = gsqrt_(:,in) * momx_(im)
1240 gsqrtmomx(:,ex) = gsqrt_(:,ex) * momx_(ip)
1241 gsqrtmomy(:,in) = gsqrt_(:,in) * momy_(im)
1242 gsqrtmomy(:,ex) = gsqrt_(:,ex) * momy_(ip)
1243 gsqrtmomz(:,in) = gsqrt_(:,in) * momz_(im)
1244 gsqrtmomz(:,ex) = gsqrt_(:,ex) * momz_(ip)
1245 gsqrtdrhot(:,in) = gsqrt_(:,in) * drhot_(im)
1246 gsqrtdrhot(:,ex) = gsqrt_(:,ex) * drhot_(ip)
1247
1248 phyd_(:,in) = pres_hyd(im)
1249 phyd_(:,ex) = pres_hyd(ip)
1250 dpres_(:,in) = dpres(im)
1251 dpres_(:,ex) = dpres(ip)
1252
1253 gsqrtdens(:,in) = gsqrtddens(:,in) + gsqrt_(:,in) * dens_hyd(im)
1254 gsqrtdens(:,ex) = gsqrtddens(:,ex) + gsqrt_(:,ex) * dens_hyd(ip)
1255
1256 gsqrtrhot(:,:) = gsqrt_(:,:) * p0ovr * (phyd_(:,:) * rp0)**rgamm + gsqrtdrhot(:,:)
1257
1258 vel(:,in) = ( gsqrtmomx(:,in) * nx(:,ke) + gsqrtmomy(:,in) * ny(:,ke) &
1259 + ( ( gsqrtmomz(:,in) * rgsqrtv(:,in) &
1260 + g13_(:,in) * gsqrtmomx(:,in) + g23_(:,in) * gsqrtmomy(:,in) ) * nz(:,ke) ) &
1261 ) / gsqrtdens(:,in)
1262 vel(:,ex) = ( gsqrtmomx(:,ex) * nx(:,ke) + gsqrtmomy(:,ex) * ny(:,ke) &
1263 + ( ( gsqrtmomz(:,ex) * rgsqrtv(:,ex) &
1264 + g13_(:,ex) * gsqrtmomx(:,ex) + g23_(:,ex) * gsqrtmomy(:,ex) ) * nz(:,ke) ) &
1265 ) / gsqrtdens(:,ex)
1266
1267 do fp=1, elem%NfpTot
1268 tmp1 = abs( nx(fp,ke) ) + abs( ny(fp,ke) )
1269 gnn_m = tmp1 &
1270 + ( 1.0_rp * rgsqrtv(fp,in)**2 + g13_(fp,in)**2 + g23_(fp,in)**2 ) * abs( nz(fp,ke) )
1271 gnn_p = tmp1 &
1272 + ( 1.0_rp * rgsqrtv(fp,ex)**2 + g13_(fp,ex)**2 + g23_(fp,ex)**2 ) * abs( nz(fp,ke) )
1273
1274 cs(fp,in) = sqrt( gnn_m * gamm * ( phyd_(fp,in) + dpres_(fp,in) ) * gsqrt_(fp,in) / gsqrtdens(fp,in) )
1275 cs(fp,ex) = sqrt( gnn_p * gamm * ( phyd_(fp,ex) + dpres_(fp,ex) ) * gsqrt_(fp,ex) / gsqrtdens(fp,ex) )
1276 alpha(fp) = max( cs(fp,in) + abs(vel(fp,in)), cs(fp,ex) + abs(vel(fp,ex)) )
1277 end do
1278 vel_n(:) = 0.5_rp * ( &
1279 ( vel(:,in) + vel(:,ex) ) &
1280 - 2.0_rp * ( dpres_(:,ex) - dpres_(:,in) ) / ( 0.5_rp * ( gsqrtdens(:,in) / gsqrt_(:,in) * cs(:,in) + gsqrtdens(:,ex) / gsqrt_(:,ex) * cs(:,ex) ) ) &
1281 )
1282
1283
1284 do fp=1, elem%NfpTot
1285 tmp1 = lmesh%Fscale(fp,ke) * 0.5_rp
1286 tmp3 = vel_n(fp) - 2.0_rp * vel(fp,in)
1287
1288 tmp2 = - abs( vel_n(fp) ) * ( gsqrtddens(fp,ex) - gsqrtddens(fp,in) )
1289 del_flux(fp,dens_vid,ke) = tmp1 * ( &
1290 gsqrtdens(fp,ex) * vel_n(fp) + gsqrtdens(fp,in) * tmp3 &
1291 + tmp2 )
1292
1293 tmp2 = - abs( vel_n(fp) ) * ( gsqrtdrhot(fp,ex) - gsqrtdrhot(fp,in) )
1294 del_flux(fp,rhot_vid,ke) = tmp1 * ( &
1295 gsqrtrhot(fp,ex) * vel_n(fp) + gsqrtrhot(fp,in) * tmp3 &
1296 + tmp2 )
1297 end do
1298
1299 do fp=1, elem%NfpTot
1300 tmp1 = 0.5_rp * ( cs(fp,in) + cs(fp,ex) )
1301 tmp3 = gsqrt_(fp,ex) * dpres_(fp,ex)
1302 tmp4 = gsqrt_(fp,in) * dpres_(fp,in)
1303
1304 del_flux_tmp_mom(fp,1) = &
1305 ( tmp3 * rgsqrtv(fp,ex) &
1306 - tmp4 * rgsqrtv(fp,in) ) * nz(fp,ke) &
1307 - tmp1 * ( gsqrtmomz(fp,ex) - gsqrtmomz(fp,in) ) !* abs(nz(fp,ke))
1308
1309 del_flux_tmp_mom(fp,2) = &
1310 ( nx(fp,ke) + g13_(fp,ex) * nz(fp,ke) ) * tmp3 &
1311 - ( nx(fp,ke) + g13_(fp,in) * nz(fp,ke) ) * tmp4 &
1312 - tmp1 * ( gsqrtmomx(fp,ex) - gsqrtmomx(fp,in) ) !* abs(nx(fp,ke))
1313
1314 del_flux_tmp_mom(fp,3) = &
1315 ( ny(fp,ke) + g23_(fp,ex) * nz(fp,ke) ) * tmp3 &
1316 - ( ny(fp,ke) + g23_(fp,in) * nz(fp,ke) ) * tmp4 &
1317 - tmp1 * ( gsqrtmomy(fp,ex) - gsqrtmomy(fp,in) ) !* abs(ny(fp,ke))
1318
1319 end do
1320 do fp=1, elem%NfpTot
1321 tmp1 = lmesh%Fscale(fp,ke) * 0.5_rp
1322 tmp3 = vel_n(fp) - 2.0_rp * vel(fp,in)
1323
1324 tmp2 = - abs( vel_n(fp) ) * ( gsqrtmomz(fp,ex) - gsqrtmomz(fp,in) )
1325 del_flux(fp,momz_vid,ke) = tmp1 * ( &
1326 gsqrtmomz(fp,ex) * vel_n(fp) + gsqrtmomz(fp,in) * tmp3 &
1327 + del_flux_tmp_mom(fp,1) &
1328 + tmp2 )
1329 end do
1330 do fp=1, elem%NfpTot
1331 tmp1 = lmesh%Fscale(fp,ke) * 0.5_rp
1332 tmp3 = vel_n(fp) - 2.0_rp * vel(fp,in)
1333
1334 tmp2 = - abs( vel_n(fp) ) * ( gsqrtmomx(fp,ex) - gsqrtmomx(fp,in) )
1335 del_flux(fp,momx_vid,ke) = tmp1 * ( &
1336 gsqrtmomx(fp,ex) * vel_n(fp) + gsqrtmomx(fp,in) * tmp3 &
1337 + del_flux_tmp_mom(fp,2) &
1338 + tmp2 )
1339
1340 tmp2 = - abs( vel_n(fp) ) * ( gsqrtmomy(fp,ex) - gsqrtmomy(fp,in) )
1341 del_flux(fp,momy_vid,ke) = tmp1 * ( &
1342 gsqrtmomy(fp,ex) * vel_n(fp) + gsqrtmomy(fp,in) * tmp3 &
1343 + del_flux_tmp_mom(fp,3) &
1344 + tmp2 )
1345 end do
1346 end do
1347
1348 return
1349 end subroutine atm_dyn_dgm_nonhydro3d_rhot_heve_numflux_get_generalvc_lmars
1350!OCL SERIAL
1352 del_flux, del_flux_hyd, & ! (out)
1353 ddens_, momx_, momy_, momz_, drhot_, dpres_, dens_hyd, pres_hyd, & ! (in)
1354 rtot, cvtot, cptot, & ! (in)
1355 gsqrt, g11, g12, g22, gsqrth, gam, g13, g23, nx, ny, nz, & ! (in)
1356 vmapm, vmapp, im2dto3d, lmesh, elem, lmesh2d, elem2d ) ! (in)
1357
1358 implicit none
1359
1360 class(localmesh3d), intent(in) :: lmesh
1361 class(elementbase3d), intent(in) :: elem
1362 class(localmesh2d), intent(in) :: lmesh2d
1363 class(elementbase2d), intent(in) :: elem2d
1364 real(rp), intent(out) :: del_flux(elem%nfptot,lmesh%ne,prgvar_num)
1365 real(rp), intent(out) :: del_flux_hyd(elem%nfptot,lmesh%ne,2)
1366 real(rp), intent(in) :: ddens_(elem%np*lmesh%nea)
1367 real(rp), intent(in) :: momx_(elem%np*lmesh%nea)
1368 real(rp), intent(in) :: momy_(elem%np*lmesh%nea)
1369 real(rp), intent(in) :: momz_(elem%np*lmesh%nea)
1370 real(rp), intent(in) :: drhot_(elem%np*lmesh%nea)
1371 real(rp), intent(in) :: dpres_(elem%np*lmesh%nea)
1372 real(rp), intent(in) :: dens_hyd(elem%np*lmesh%nea)
1373 real(rp), intent(in) :: pres_hyd(elem%np*lmesh%nea)
1374 real(rp), intent(in) :: rtot (elem%np*lmesh%nea)
1375 real(rp), intent(in) :: cvtot(elem%np*lmesh%nea)
1376 real(rp), intent(in) :: cptot(elem%np*lmesh%nea)
1377 real(rp), intent(in) :: gsqrt(elem%np*lmesh%nea)
1378 real(rp), intent(in) :: g11(elem2d%np,lmesh2d%ne)
1379 real(rp), intent(in) :: g12(elem2d%np,lmesh2d%ne)
1380 real(rp), intent(in) :: g22(elem2d%np,lmesh2d%ne)
1381 real(rp), intent(in) :: gsqrth(elem2d%np,lmesh2d%ne)
1382 real(rp), intent(in) :: gam(elem%np*lmesh%nea)
1383 real(rp), intent(in) :: g13(elem%np*lmesh%nea)
1384 real(rp), intent(in) :: g23(elem%np*lmesh%nea)
1385 real(rp), intent(in) :: nx(elem%nfptot,lmesh%ne)
1386 real(rp), intent(in) :: ny(elem%nfptot,lmesh%ne)
1387 real(rp), intent(in) :: nz(elem%nfptot,lmesh%ne)
1388 integer, intent(in) :: vmapm(elem%nfptot,lmesh%ne)
1389 integer, intent(in) :: vmapp(elem%nfptot,lmesh%ne)
1390 integer, intent(in) :: im2dto3d(elem%nfptot)
1391
1392 integer :: ke, ip(elem%nfptot), im(elem%nfptot)
1393 integer :: ke2d
1394 real(rp) :: velp(elem%nfptot), velm(elem%nfptot), alpha(elem%nfptot)
1395 real(rp) :: dpresp(elem%nfptot), dpresm(elem%nfptot)
1396 real(rp) :: gsqrtdensm(elem%nfptot), gsqrtdensp(elem%nfptot)
1397 real(rp) :: gsqrtrhotm(elem%nfptot), gsqrtrhotp(elem%nfptot)
1398 real(rp) :: gsqrtddens_p(elem%nfptot), gsqrtddens_m(elem%nfptot)
1399 real(rp) :: gsqrtmomx_p(elem%nfptot), gsqrtmomx_m(elem%nfptot)
1400 real(rp) :: gsqrtmomy_p(elem%nfptot), gsqrtmomy_m(elem%nfptot)
1401 real(rp) :: gsqrtmomz_p(elem%nfptot), gsqrtmomz_m(elem%nfptot)
1402 real(rp) :: gsqrtdrhot_p(elem%nfptot), gsqrtdrhot_m(elem%nfptot)
1403 real(rp) :: phyd_p(elem%nfptot), phyd_m(elem%nfptot)
1404 real(rp) :: gsqrt_p(elem%nfptot), gsqrt_m(elem%nfptot)
1405 real(rp) :: gsqrtv_p(elem%nfptot), gsqrtv_m(elem%nfptot)
1406 real(rp) :: g13_m(elem%nfptot), g13_p(elem%nfptot)
1407 real(rp) :: g23_m(elem%nfptot), g23_p(elem%nfptot)
1408 real(rp) :: g1n_m(elem%nfptot), g2n_m(elem%nfptot)
1409 real(rp) :: gnn_m(elem%nfptot), gnn_p(elem%nfptot)
1410 real(rp) :: gxz_m(elem%nfptot), gxz_p(elem%nfptot)
1411 real(rp) :: gyz_m(elem%nfptot), gyz_p(elem%nfptot)
1412 real(rp) :: rgam2_m(elem%nfptot), rgam2_p(elem%nfptot)
1413
1414 real(rp) :: gamm, rgamm
1415 real(rp) :: rp0
1416 real(rp) :: rovp0, p0ovr
1417
1418 !------------------------------------------------------------------------
1419
1420 gamm = cpdry / cvdry
1421 rgamm = cvdry / cpdry
1422 rp0 = 1.0_rp / pres00
1423 rovp0 = rdry * rp0
1424 p0ovr = pres00 / rdry
1425
1426 !$omp parallel do private( &
1427 !$omp ke, iM, iP, ke2d, &
1428 !$omp alpha, VelM, VelP, &
1429 !$omp dpresM, dpresP, GsqrtDensM, GsqrtDensP, GsqrtRhotM, GsqrtRhotP, &
1430 !$omp GsqrtMOMX_M, GsqrtMOMX_P, GsqrtMOMY_M, GsqrtMOMY_P, GsqrtMOMZ_M, GsqrtMOMZ_P, &
1431 !$omp GsqrtDDENS_M, GsqrtDDENS_P, GsqrtDRHOT_M, GsqrtDRHOT_P, &
1432 !$omp Phyd_M, Phyd_P, &
1433 !$omp Gsqrt_P, Gsqrt_M, GsqrtV_P, GsqrtV_M, G13_P, G13_M, G23_P, G23_M, &
1434 !$omp rgam2_M, rgam2_P, Gxz_P, Gxz_M, Gyz_P, Gyz_M, G1n_M, G2n_M, Gnn_P, Gnn_M )
1435 do ke=lmesh%NeS, lmesh%NeE
1436 im(:) = vmapm(:,ke); ip(:) = vmapp(:,ke)
1437 ke2d = lmesh%EMap3Dto2D(ke)
1438
1439 gsqrt_m(:) = gsqrt(im)
1440 gsqrt_p(:) = gsqrt(ip)
1441
1442 rgam2_m(:) = 1.0_rp / gam(im)**2
1443 rgam2_p(:) = 1.0_rp / gam(ip)**2
1444 gsqrtv_m(:) = gsqrt_m(:) * rgam2_m(:) / gsqrth(im2dto3d(:),ke2d)
1445 gsqrtv_p(:) = gsqrt_p(:) * rgam2_p(:) / gsqrth(im2dto3d(:),ke2d)
1446
1447 g13_m(:) = g13(im)
1448 g13_p(:) = g13(ip)
1449 g23_m(:) = g23(im)
1450 g23_p(:) = g23(ip)
1451
1452 gsqrtddens_m(:) = gsqrt_m(:) * ddens_(im)
1453 gsqrtddens_p(:) = gsqrt_p(:) * ddens_(ip)
1454 gsqrtmomx_m(:) = gsqrt_m(:) * momx_(im)
1455 gsqrtmomx_p(:) = gsqrt_p(:) * momx_(ip)
1456 gsqrtmomy_m(:) = gsqrt_m(:) * momy_(im)
1457 gsqrtmomy_p(:) = gsqrt_p(:) * momy_(ip)
1458 gsqrtmomz_m(:) = gsqrt_m(:) * momz_(im)
1459 gsqrtmomz_p(:) = gsqrt_p(:) * momz_(ip)
1460 gsqrtdrhot_m(:) = gsqrt_m(:) * drhot_(im)
1461 gsqrtdrhot_p(:) = gsqrt_p(:) * drhot_(ip)
1462 phyd_m(:) = pres_hyd(im)
1463 phyd_p(:) = pres_hyd(ip)
1464
1465 gxz_m(:) = rgam2_m(:) * ( g11(im2dto3d(:),ke2d) * g13_m(:) + g12(im2dto3d(:),ke2d) * g23_m(:) )
1466 gxz_p(:) = rgam2_p(:) * ( g11(im2dto3d(:),ke2d) * g13_p(:) + g12(im2dto3d(:),ke2d) * g23_p(:) )
1467
1468 gyz_m(:) = rgam2_m(:) * ( g12(im2dto3d(:),ke2d) * g13_m(:) + g22(im2dto3d(:),ke2d) * g23_m(:) )
1469 gyz_p(:) = rgam2_p(:) * ( g12(im2dto3d(:),ke2d) * g13_p(:) + g22(im2dto3d(:),ke2d) * g23_p(:) )
1470
1471 g1n_m(:) = rgam2_m(:) * ( g11(im2dto3d(:),ke2d) * nx(:,ke) + g12(im2dto3d(:),ke2d) * ny(:,ke) )
1472 g2n_m(:) = rgam2_p(:) * ( g12(im2dto3d(:),ke2d) * nx(:,ke) + g22(im2dto3d(:),ke2d) * ny(:,ke) )
1473
1474 gnn_m(:) = rgam2_m(:) * ( g11(im2dto3d(:),ke2d) * abs( nx(:,ke) ) + g22(im2dto3d(:),ke2d) * abs( ny(:,ke) ) ) &
1475 + ( 1.0_rp / gsqrtv_m(:)**2 + g13_m(:) * gxz_m(:) + g23_m(:) * gyz_m(:) ) * abs( nz(:,ke) )
1476 gnn_p(:) = rgam2_p(:) * ( g11(im2dto3d(:),ke2d) * abs( nx(:,ke) ) + g22(im2dto3d(:),ke2d) * abs( ny(:,ke) ) ) &
1477 + ( 1.0_rp / gsqrtv_p(:)**2 + g13_p(:) * gxz_p(:) + g23_p(:) * gyz_p(:) ) * abs( nz(:,ke) )
1478
1479 gsqrtdensm(:) = gsqrtddens_m(:) + gsqrt_m(:) * dens_hyd(im)
1480 gsqrtdensp(:) = gsqrtddens_p(:) + gsqrt_p(:) * dens_hyd(ip)
1481
1482 gsqrtrhotm(:) = gsqrt_m(:) * p0ovr * (phyd_m(:) * rp0)**rgamm + gsqrtdrhot_m(:)
1483 gsqrtrhotp(:) = gsqrt_p(:) * p0ovr * (phyd_p(:) * rp0)**rgamm + gsqrtdrhot_p(:)
1484
1485 velm(:) = ( gsqrtmomx_m(:) * nx(:,ke) + gsqrtmomy_m(:) * ny(:,ke) &
1486 + ( ( gsqrtmomz_m(:) / gsqrtv_m(:) &
1487 + g13_m(:) * gsqrtmomx_m(:) + g23_m(:) * gsqrtmomy_m(:) ) * nz(:,ke) ) &
1488 ) / gsqrtdensm(:)
1489 velp(:) = ( gsqrtmomx_p(:) * nx(:,ke) + gsqrtmomy_p(:) * ny(:,ke) &
1490 + ( ( gsqrtmomz_p(:) / gsqrtv_p(:) &
1491 + g13_p(:) * gsqrtmomx_p(:) + g23_p(:) * gsqrtmomy_p(:) ) * nz(:,ke) ) &
1492 ) / gsqrtdensp(:)
1493
1494 ! dpresM(:) = PRES00 * ( Rtot(iM) * rP0 * GsqrtRhotM(:) / Gsqrt_M(:) )**( CPtot(iM) / CVtot(iM) ) &
1495 ! - Phyd_M(:)
1496 ! dpresP(:) = PRES00 * ( Rtot(iP) * rP0 * GsqrtRhotP(:) / Gsqrt_P(:) )**( CPtot(iP) / CVtot(iP) ) &
1497 ! - Phyd_P(:)
1498 dpresm(:) = dpres_(im)
1499 dpresp(:) = dpres_(ip)
1500
1501 alpha(:) = max( sqrt( gnn_m(:) * gamm * ( phyd_m(:) + dpresm(:) ) * gsqrt_m(:) / gsqrtdensm(:) ) + abs(velm(:)), &
1502 sqrt( gnn_p(:) * gamm * ( phyd_p(:) + dpresp(:) ) * gsqrt_p(:) / gsqrtdensp(:) ) + abs(velp(:)) )
1503
1504 del_flux(:,ke,dens_vid) = 0.5_rp * ( &
1505 ( gsqrtdensp(:) * velp(:) - gsqrtdensm(:) * velm(:) ) &
1506 - alpha(:) * ( gsqrtddens_p(:) - gsqrtddens_m(:) ) )
1507
1508 del_flux(:,ke,momx_vid ) = 0.5_rp * ( &
1509 ( gsqrtmomx_p(:) * velp(:) - gsqrtmomx_m(:) * velm(:) ) &
1510 + ( gsqrt_p(:) * ( g1n_m(:) + gxz_p(:) * nz(:,ke)) * dpresp(:) &
1511 - gsqrt_m(:) * ( g1n_m(:) + gxz_m(:) * nz(:,ke)) * dpresm(:) ) &
1512 - alpha(:) * ( gsqrtmomx_p(:) - gsqrtmomx_m(:) ) )
1513
1514 del_flux(:,ke,momy_vid ) = 0.5_rp * ( &
1515 ( gsqrtmomy_p(:) * velp(:) - gsqrtmomy_m(:) * velm(:) ) &
1516 + ( gsqrt_p(:) * ( g2n_m(:) + gyz_p(:) * nz(:,ke) ) * dpresp(:) &
1517 - gsqrt_m(:) * ( g2n_m(:) + gyz_m(:) * nz(:,ke) ) * dpresm(:) ) &
1518 - alpha(:) * ( gsqrtmomy_p(:) - gsqrtmomy_m(:) ) )
1519
1520 del_flux(:,ke,momz_vid ) = 0.5_rp * ( &
1521 ( gsqrtmomz_p(:) * velp(:) - gsqrtmomz_m(:) * velm(:) ) &
1522 + ( gsqrt_p(:) * dpresp(:) / gsqrtv_p(:) &
1523 - gsqrt_m(:) * dpresm(:) / gsqrtv_m(:) ) * nz(:,ke) &
1524 - alpha(:) * ( gsqrtmomz_p(:) - gsqrtmomz_m(:) ) )
1525
1526 del_flux(:,ke,rhot_vid) = 0.5_rp * ( &
1527 ( gsqrtrhotp(:) * velp(:) - gsqrtrhotm(:) * velm(:) ) &
1528 - alpha(:) * ( gsqrtdrhot_p(:) - gsqrtdrhot_m(:) ) )
1529
1530 del_flux_hyd(:,ke,1) = 0.5_rp * ( &
1531 gsqrtv_p(:) * ( nx(:,ke) + g13_p(:) * nz(:,ke) ) * phyd_p(:) &
1532 - gsqrtv_m(:) * ( nx(:,ke) + g13_m(:) * nz(:,ke) ) * phyd_m(:) )
1533
1534 del_flux_hyd(:,ke,2) = 0.5_rp * ( &
1535 gsqrtv_p(:) * ( ny(:,ke) + g23_p(:) * nz(:,ke) ) * phyd_p(:) &
1536 - gsqrtv_m(:) * ( ny(:,ke) + g23_m(:) * nz(:,ke) ) * phyd_m(:) )
1537 end do
1538
1539 return
1541
1542!OCL SERIAL
1544 del_flux, & ! (out)
1545 ddens_, momx_, momy_, momz_, drhot_, dpres, & ! (in)
1546 dens_hyd, pres_hyd, therm_hyd, rtot, cvtot, cptot, & ! (in)
1547 gsqrt, g11, g12, g22, gsqrth, gam, g13, g23, nx, ny, nz, & ! (in)
1548 vmapm, vmapp, im2dto3d, lmesh, elem, lmesh2d, elem2d ) ! (in)
1549
1550 implicit none
1551
1552 class(localmesh3d), intent(in) :: lmesh
1553 class(elementbase3d), intent(in) :: elem
1554 class(localmesh2d), intent(in) :: lmesh2d
1555 class(elementbase2d), intent(in) :: elem2d
1556 real(rp), intent(out) :: del_flux(elem%nfptot,prgvar_num,lmesh%ne)
1557 real(rp), intent(in) :: ddens_(elem%np*lmesh%nea)
1558 real(rp), intent(in) :: momx_(elem%np*lmesh%nea)
1559 real(rp), intent(in) :: momy_(elem%np*lmesh%nea)
1560 real(rp), intent(in) :: momz_(elem%np*lmesh%nea)
1561 real(rp), intent(in) :: drhot_(elem%np*lmesh%nea)
1562 real(rp), intent(in) :: dpres(elem%np*lmesh%nea)
1563 real(rp), intent(in) :: dens_hyd(elem%np*lmesh%nea)
1564 real(rp), intent(in) :: pres_hyd(elem%np*lmesh%nea)
1565 real(rp), intent(in) :: therm_hyd(elem%np*lmesh%nea)
1566 real(rp), intent(in) :: rtot (elem%np*lmesh%nea)
1567 real(rp), intent(in) :: cvtot(elem%np*lmesh%nea)
1568 real(rp), intent(in) :: cptot(elem%np*lmesh%nea)
1569 real(rp), intent(in) :: gsqrt(elem%np*lmesh%nea)
1570 real(rp), intent(in) :: g11(elem2d%np,lmesh2d%ne)
1571 real(rp), intent(in) :: g12(elem2d%np,lmesh2d%ne)
1572 real(rp), intent(in) :: g22(elem2d%np,lmesh2d%ne)
1573 real(rp), intent(in) :: gsqrth(elem2d%np,lmesh2d%ne)
1574 real(rp), intent(in) :: gam(elem%np*lmesh%nea)
1575 real(rp), intent(in) :: g13(elem%np*lmesh%nea)
1576 real(rp), intent(in) :: g23(elem%np*lmesh%nea)
1577 real(rp), intent(in) :: nx(elem%nfptot,lmesh%ne)
1578 real(rp), intent(in) :: ny(elem%nfptot,lmesh%ne)
1579 real(rp), intent(in) :: nz(elem%nfptot,lmesh%ne)
1580 integer, intent(in) :: vmapm(elem%nfptot,lmesh%ne)
1581 integer, intent(in) :: vmapp(elem%nfptot,lmesh%ne)
1582 integer, intent(in) :: im2dto3d(elem%nfptot)
1583
1584 integer :: ke, fp, i, ip(elem%nfptot), im(elem%nfptot)
1585 integer :: ke2d
1586 real(rp) :: vel(elem%nfptot,2), alpha(elem%nfptot)
1587 real(rp) :: dpres_(elem%nfptot,2)
1588 real(rp) :: gsqrtdens(elem%nfptot,2)
1589 real(rp) :: gsqrtrhot(elem%nfptot,2)
1590 real(rp) :: gsqrtddens(elem%nfptot,2)
1591 real(rp) :: gsqrtmomx(elem%nfptot,2)
1592 real(rp) :: gsqrtmomy(elem%nfptot,2)
1593 real(rp) :: gsqrtmomz(elem%nfptot,2)
1594 real(rp) :: gsqrtdrhot(elem%nfptot,2)
1595 real(rp) :: phyd_(elem%nfptot,2)
1596 real(rp) :: gsqrt_(elem%nfptot,2)
1597 real(rp) :: gsqrtv_(elem%nfptot,2)
1598 real(rp) :: rgsqrtv(elem%nfptot,2)
1599 real(rp) :: g13_(elem%nfptot,2)
1600 real(rp) :: g23_(elem%nfptot,2)
1601 real(rp) :: gxz_(elem%nfptot,2)
1602 real(rp) :: gyz_(elem%nfptot,2)
1603 real(rp) :: g1n_(elem%nfptot,2)
1604 real(rp) :: g2n_(elem%nfptot,2)
1605
1606 real(rp) :: gnn_m, gnn_p
1607 real(rp) :: rgam2(elem%nfptot,2)
1608
1609 real(rp) :: gamm, rgamm
1610 real(rp) :: rp0
1611 real(rp) :: rovp0, p0ovr
1612
1613 integer, parameter :: in = 1
1614 integer, parameter :: ex = 2
1615
1616 real(rp) :: tmp1, tmp2, tmp3, tmp4, tmp01, tmp02, tmp03
1617 real(rp) :: del_flux_tmp_mom(elem%nfptot,3)
1618
1619 real(rp) :: g11_, g12_, g22_
1620 !------------------------------------------------------------------------
1621
1622 gamm = cpdry / cvdry
1623 rgamm = cvdry / cpdry
1624 rp0 = 1.0_rp / pres00
1625 rovp0 = rdry * rp0
1626 p0ovr = pres00 / rdry
1627
1628 !$omp parallel do private( &
1629 !$omp ke, iM, iP, ke2D, fp, &
1630 !$omp alpha, Vel, &
1631 !$omp dpres_, GsqrtDens, GsqrtRhot, &
1632 !$omp GsqrtDDENS, GsqrtMOMX, GsqrtMOMY, GsqrtMOMZ, GsqrtDRHOT, &
1633 !$omp Phyd_, &
1634 !$omp Gsqrt_, GsqrtV_, RGsqrtV, G11_, G12_, G22_, G13_, G23_, rgam2, &
1635 !$omp Gxz_, Gyz_, G1n_, G2n_, &
1636 !$omp Gnn_P, Gnn_M, tmp1, tmp2, tmp3, tmp4, tmp01, tmp02, tmp03, del_flux_tmp_mom )
1637!OCL PREFETCH
1638 do ke=lmesh%NeS, lmesh%NeE
1639 im(:) = vmapm(:,ke); ip(:) = vmapp(:,ke)
1640 ke2d = lmesh%EMap3Dto2D(ke)
1641
1642 gsqrt_(:,in) = gsqrt(im)
1643 gsqrt_(:,ex) = gsqrt(ip)
1644
1645 rgam2(:,in) = 1.0_rp / gam(im)**2
1646 rgam2(:,ex) = 1.0_rp / gam(ip)**2
1647 gsqrtv_(:,in) = gsqrt_(:,in) * rgam2(:,in) / gsqrth(im2dto3d(:),ke2d)
1648 gsqrtv_(:,ex) = gsqrt_(:,ex) * rgam2(:,ex) / gsqrth(im2dto3d(:),ke2d)
1649 rgsqrtv(:,:) = 1.0_rp / gsqrtv_(:,:)
1650
1651 g13_(:,in) = g13(im)
1652 g13_(:,ex) = g13(ip)
1653 g23_(:,in) = g23(im)
1654 g23_(:,ex) = g23(ip)
1655
1656 gsqrtddens(:,in) = gsqrt_(:,in) * ddens_(im)
1657 gsqrtddens(:,ex) = gsqrt_(:,ex) * ddens_(ip)
1658 gsqrtmomx(:,in) = gsqrt_(:,in) * momx_(im)
1659 gsqrtmomx(:,ex) = gsqrt_(:,ex) * momx_(ip)
1660 gsqrtmomy(:,in) = gsqrt_(:,in) * momy_(im)
1661 gsqrtmomy(:,ex) = gsqrt_(:,ex) * momy_(ip)
1662 gsqrtmomz(:,in) = gsqrt_(:,in) * momz_(im)
1663 gsqrtmomz(:,ex) = gsqrt_(:,ex) * momz_(ip)
1664 gsqrtdrhot(:,in) = gsqrt_(:,in) * drhot_(im)
1665 gsqrtdrhot(:,ex) = gsqrt_(:,ex) * drhot_(ip)
1666
1667 phyd_(:,in) = pres_hyd(im)
1668 phyd_(:,ex) = pres_hyd(ip)
1669 dpres_(:,in) = dpres(im)
1670 dpres_(:,ex) = dpres(ip)
1671
1672 gsqrtdens(:,in) = gsqrtddens(:,in) + gsqrt_(:,in) * dens_hyd(im)
1673 gsqrtdens(:,ex) = gsqrtddens(:,ex) + gsqrt_(:,ex) * dens_hyd(ip)
1674
1675 gsqrtrhot(:,in) = gsqrt_(:,in) * therm_hyd(im) + gsqrtdrhot(:,in)
1676 gsqrtrhot(:,ex) = gsqrt_(:,ex) * therm_hyd(ip) + gsqrtdrhot(:,ex)
1677
1678 vel(:,in) = ( gsqrtmomx(:,in) * nx(:,ke) + gsqrtmomy(:,in) * ny(:,ke) &
1679 + ( ( gsqrtmomz(:,in) * rgsqrtv(:,in) &
1680 + g13_(:,in) * gsqrtmomx(:,in) + g23_(:,in) * gsqrtmomy(:,in) ) * nz(:,ke) ) &
1681 ) / gsqrtdens(:,in)
1682 vel(:,ex) = ( gsqrtmomx(:,ex) * nx(:,ke) + gsqrtmomy(:,ex) * ny(:,ke) &
1683 + ( ( gsqrtmomz(:,ex) * rgsqrtv(:,ex) &
1684 + g13_(:,ex) * gsqrtmomx(:,ex) + g23_(:,ex) * gsqrtmomy(:,ex) ) * nz(:,ke) ) &
1685 ) / gsqrtdens(:,ex)
1686
1687 do fp=1, elem%NfpTot
1688 g11_ = g11(im2dto3d(fp),ke2d); g12_ = g12(im2dto3d(fp),ke2d); g22_ = g22(im2dto3d(fp),ke2d)
1689
1690 gxz_(fp,in) = rgam2(fp,in) * ( g11_ * g13_(fp,in) + g12_ * g23_(fp,in) )
1691 gxz_(fp,ex) = rgam2(fp,ex) * ( g11_ * g13_(fp,ex) + g12_ * g23_(fp,ex) )
1692
1693 gyz_(fp,in) = rgam2(fp,in) * ( g12_ * g13_(fp,in) + g22_ * g23_(fp,in) )
1694 gyz_(fp,ex) = rgam2(fp,ex) * ( g12_ * g13_(fp,ex) + g22_ * g23_(fp,ex) )
1695
1696 g1n_(fp,in) = rgam2(fp,in) * ( g11_ * nx(fp,ke) + g12_ * ny(fp,ke) )
1697 g1n_(fp,ex) = rgam2(fp,ex) * ( g11_ * nx(fp,ke) + g12_ * ny(fp,ke) )
1698
1699 g2n_(fp,in) = rgam2(fp,in) * ( g12_ * nx(fp,ke) + g22_ * ny(fp,ke) )
1700 g2n_(fp,ex) = rgam2(fp,ex) * ( g12_ * nx(fp,ke) + g22_ * ny(fp,ke) )
1701 end do
1702 do fp=1, elem%NfpTot
1703 g11_ = g11(im2dto3d(fp),ke2d); g22_ = g22(im2dto3d(fp),ke2d)
1704 tmp1 = abs( g11_ * nx(fp,ke) ) + abs( g22_ * ny(fp,ke) )
1705
1706 gnn_m = rgam2(fp,in) * tmp1 &
1707 + ( 1.0_rp * rgsqrtv(fp,in)**2 + g13_(fp,in) * gxz_(fp,in) + g23_(fp,in) * gyz_(fp,in) ) * abs( nz(fp,ke) )
1708 gnn_p = rgam2(fp,ex) * tmp1 &
1709 + ( 1.0_rp * rgsqrtv(fp,ex)**2 + g13_(fp,ex) * gxz_(fp,ex) + g23_(fp,ex) * gyz_(fp,ex) ) * abs( nz(fp,ke) )
1710
1711 alpha(fp) = max( sqrt( gnn_m * gamm * ( phyd_(fp,in) + dpres_(fp,in) ) * gsqrt_(fp,in) / gsqrtdens(fp,in) ) + abs(vel(fp,in)), &
1712 sqrt( gnn_p * gamm * ( phyd_(fp,ex) + dpres_(fp,ex) ) * gsqrt_(fp,ex) / gsqrtdens(fp,ex) ) + abs(vel(fp,ex)) )
1713 end do
1714
1715 do fp=1, elem%NfpTot
1716 tmp1 = lmesh%Fscale(fp,ke) * 0.5_rp
1717
1718 tmp2 = - alpha(fp) * ( gsqrtddens(fp,ex) - gsqrtddens(fp,in) )
1719 del_flux(fp,dens_vid,ke) = tmp1 * ( &
1720 gsqrtdens(fp,ex) * vel(fp,ex) - gsqrtdens(fp,in) * vel(fp,in) &
1721 + tmp2 )
1722
1723 tmp2 = - alpha(fp) * ( gsqrtdrhot(fp,ex) - gsqrtdrhot(fp,in) )
1724 del_flux(fp,rhot_vid,ke) = tmp1 * ( &
1725 gsqrtrhot(fp,ex) * vel(fp,ex) - gsqrtrhot(fp,in) * vel(fp,in) &
1726 + tmp2 )
1727 end do
1728
1729 do fp=1, elem%NfpTot
1730 tmp3 = gsqrt_(fp,ex) * dpres_(fp,ex)
1731 tmp4 = gsqrt_(fp,in) * dpres_(fp,in)
1732
1733 del_flux_tmp_mom(fp,1) = &
1734 ( tmp3 * rgsqrtv(fp,ex) &
1735 - tmp4 * rgsqrtv(fp,in) ) * nz(fp,ke)
1736
1737 del_flux_tmp_mom(fp,2) = &
1738 ( g1n_(fp,ex) + gxz_(fp,ex) * nz(fp,ke) ) * tmp3 &
1739 - ( g1n_(fp,in) + gxz_(fp,in) * nz(fp,ke) ) * tmp4
1740
1741 del_flux_tmp_mom(fp,3) = &
1742 ( g2n_(fp,ex) + gyz_(fp,ex) * nz(fp,ke) ) * tmp3 &
1743 - ( g2n_(fp,in) + gyz_(fp,in) * nz(fp,ke) ) * tmp4
1744 end do
1745 do fp=1, elem%NfpTot
1746 tmp1 = lmesh%Fscale(fp,ke) * 0.5_rp
1747
1748 tmp2 = - alpha(fp) * ( gsqrtmomz(fp,ex) - gsqrtmomz(fp,in) )
1749 del_flux(fp,momz_vid,ke) = tmp1 * ( &
1750 gsqrtmomz(fp,ex) * vel(fp,ex) - gsqrtmomz(fp,in) * vel(fp,in) &
1751 + del_flux_tmp_mom(fp,1) &
1752 + tmp2 )
1753 end do
1754 do fp=1, elem%NfpTot
1755 tmp1 = lmesh%Fscale(fp,ke) * 0.5_rp
1756
1757 tmp2 = - alpha(fp) * ( gsqrtmomx(fp,ex) - gsqrtmomx(fp,in) )
1758 del_flux(fp,momx_vid,ke) = tmp1 * ( &
1759 gsqrtmomx(fp,ex) * vel(fp,ex) - gsqrtmomx(fp,in) * vel(fp,in) &
1760 + del_flux_tmp_mom(fp,2) &
1761 + tmp2 )
1762
1763 tmp2 = - alpha(fp) * ( gsqrtmomy(fp,ex) - gsqrtmomy(fp,in) )
1764 del_flux(fp,momy_vid,ke) = tmp1 * ( &
1765 gsqrtmomy(fp,ex) * vel(fp,ex) - gsqrtmomy(fp,in) * vel(fp,in) &
1766 + del_flux_tmp_mom(fp,3) &
1767 + tmp2 )
1768 end do
1769 end do
1770
1771 return
1773
module FElib / Fluid dyn solver / Atmosphere / Nonhydrostatic model / Common
module FElib / Fluid dyn solver / Atmosphere / Nonhydrostatic model / HEVE / Numflux
subroutine, public atm_dyn_dgm_nonhydro3d_rhot_heve_numflux_get_generalvc(del_flux, ddens_, momx_, momy_, momz_, drhot_, dpres, dens_hyd, pres_hyd, therm_hyd, rtot, cvtot, cptot, gsqrt, g13, g23, nx, ny, nz, vmapm, vmapp, lmesh, elem, lmesh2d, elem2d)
subroutine, public atm_dyn_dgm_nonhydro3d_rhot_heve_numflux_get_generalvc_asis(del_flux, del_flux_hyd, ddens_, momx_, momy_, momz_, drhot_, dpres_, dens_hyd, pres_hyd, rtot, cvtot, cptot, gsqrt, g13, g23, nx, ny, nz, vmapm, vmapp, lmesh, elem, lmesh2d, elem2d)
subroutine, public atm_dyn_dgm_nonhydro3d_rhot_heve_add_bnd_contrib_generalvc(prgvar_dt, ddens_, momx_, momy_, momz_, drhot_, dpres, dens_hyd, pres_hyd, therm_hyd, rtot, cvtot, cptot, gsqrt, g13, g23, nx, ny, nz, vmapm, vmapp, lmesh, elem, lmesh2d, elem2d, elem3d_optr)
subroutine, public atm_dyn_dgm_nonhydro3d_rhot_heve_numflux_get_generalhvc_asis(del_flux, del_flux_hyd, ddens_, momx_, momy_, momz_, drhot_, dpres_, dens_hyd, pres_hyd, rtot, cvtot, cptot, gsqrt, g11, g12, g22, gsqrth, gam, g13, g23, nx, ny, nz, vmapm, vmapp, im2dto3d, lmesh, elem, lmesh2d, elem2d)
subroutine, public atm_dyn_dgm_nonhydro3d_rhot_heve_numflux_get_generalhvc(del_flux, ddens_, momx_, momy_, momz_, drhot_, dpres, dens_hyd, pres_hyd, therm_hyd, rtot, cvtot, cptot, gsqrt, g11, g12, g22, gsqrth, gam, g13, g23, nx, ny, nz, vmapm, vmapp, im2dto3d, lmesh, elem, lmesh2d, elem2d)
module FElib / Element / Base
module FElib / Element / hexahedron
module FElib / Element / Operation / Base
module FElib / Mesh / Local 2D
module FElib / Mesh / Local 3D
module FElib / Mesh / Base 3D
module FElib / Data / base
Module common / sparsemat.
Derived type representing a 2D reference element.
Derived type representing a 3D reference element.
Derived type representing a hexahedral element.
Derived type representing a local mesh for 2D domain.
Derived type to manage a local 3D computational domain.
Derived type representing a field with 3D local mesh.
Derived type to manage a computational mesh (base type for 3D domain)
Derived type representing a field with 3D mesh.