FE-Project
Loading...
Searching...
No Matches
scale_atm_dyn_dgm_globalnonhydro3d_rhot_heve.F90
Go to the documentation of this file.
1!-------------------------------------------------------------------------------
2!> module FElib / Fluid dyn solver / Atmosphere / Global nonhydrostatic model / HEVE
3!!
4!! @par Description
5!! HEVE DGM scheme for Global Atmospheric Dynamical process.
6!! The governing equations is a fully compressible nonhydrostatic equations,
7!! which consist of mass, momentum, and thermodynamics (density * potential temperature conservation) equations.
8!!
9!! @author Yuta Kawai, Xuanzhengbo Ren, and Team SCALE
10!<
11!-------------------------------------------------------------------------------
12#include "scaleFElib.h"
14 !-----------------------------------------------------------------------------
15 !
16 !++ Used modules
17 !
18 use scale_precision
19 use scale_io
20 use scale_prc
21 use scale_prof
22 use scale_const, only: &
23 grav => const_grav, &
24 rdry => const_rdry, &
25 cpdry => const_cpdry, &
26 cvdry => const_cvdry, &
27 pres00 => const_pre00, &
28 rplanet => const_radius
29
31 use scale_element_base, only: &
41
45 dens_vid => prgvar_ddens_id, rhot_vid => prgvar_drhot_id, &
46 momx_vid => prgvar_momx_id, momy_vid => prgvar_momy_id, &
47 momz_vid => prgvar_momz_id, &
49
50 !-----------------------------------------------------------------------------
51 implicit none
52 private
53 !-----------------------------------------------------------------------------
54 !
55 !++ Public procedures
56 !
62
63 !-----------------------------------------------------------------------------
64 !
65 !++ Public parameters & variables
66 !
67
68 !-----------------------------------------------------------------------------
69 !
70 !++ Private procedures & variables
71 !
72 !-------------------
73
74contains
75!OCL SERIAL
77 implicit none
78 class(meshbase3d), intent(in) :: mesh
79 !--------------------------------------------
80
82
83 return
85
86!OCL SERIAL
88 implicit none
89 !--------------------------------------------
90
92 return
94
95 !-------------------------------
96
97!OCL SERIAL
99 DENS_dt, MOMX_dt, MOMY_dt, MOMZ_dt, RHOT_dt, & ! (out)
100 ddens_, momx_, momy_, momz_, drhot_, dpres_, & ! (in)
101 dens_hyd, pres_hyd, pres_hyd_ref, therm_hyd, & ! (in)
102 coriolis, rtot, cvtot, cptot, dphyddx, dphyddy, & ! (in)
103 element3d_operation, dx, dy, dz, sx, sy, sz, lift, & ! (in)
104 lmesh, elem, lmesh2d, elem2d ) ! (in)
105
108 use scale_const, only: &
109 ohm => const_ohm
110 implicit none
111
112 class(localmesh3d), intent(in) :: lmesh
113 class(elementbase3d), intent(in) :: elem
114 class(localmesh2d), intent(in) :: lmesh2d
115 class(elementbase2d), intent(in) :: elem2d
116 class(elementoperationbase3d), intent(in) :: element3d_operation
117 type(sparsemat), intent(in) :: dx, dy, dz, sx, sy, sz, lift
118 real(rp), intent(out) :: dens_dt(elem%np,lmesh%nea)
119 real(rp), intent(out) :: momx_dt(elem%np,lmesh%nea)
120 real(rp), intent(out) :: momy_dt(elem%np,lmesh%nea)
121 real(rp), intent(out) :: momz_dt(elem%np,lmesh%nea)
122 real(rp), intent(out) :: rhot_dt(elem%np,lmesh%nea)
123 real(rp), intent(in) :: ddens_(elem%np,lmesh%nea)
124 real(rp), intent(in) :: momx_(elem%np,lmesh%nea)
125 real(rp), intent(in) :: momy_(elem%np,lmesh%nea)
126 real(rp), intent(in) :: momz_(elem%np,lmesh%nea)
127 real(rp), intent(in) :: drhot_(elem%np,lmesh%nea)
128 real(rp), intent(in) :: dpres_(elem%np,lmesh%nea)
129 real(rp), intent(in) :: dens_hyd(elem%np,lmesh%nea)
130 real(rp), intent(in) :: pres_hyd(elem%np,lmesh%nea)
131 real(rp), intent(in) :: pres_hyd_ref(elem%np,lmesh%nea)
132 real(rp), intent(in) :: therm_hyd(elem%np,lmesh%nea)
133 real(rp), intent(in) :: coriolis(elem2d%np,lmesh2d%nea)
134 real(rp), intent(in) :: rtot (elem%np,lmesh%nea)
135 real(rp), intent(in) :: cvtot(elem%np,lmesh%nea)
136 real(rp), intent(in) :: cptot(elem%np,lmesh%nea)
137 real(rp), intent(in) :: dphyddx(elem%np,lmesh%nea)
138 real(rp), intent(in) :: dphyddy(elem%np,lmesh%nea)
139
140 real(rp) :: fx(elem%np), fy(elem%np), fz(elem%np), liftdelflx(elem%np)
141 real(rp) :: dpres_hyd(elem%np), gradphyd_x(elem%np), gradphyd_y(elem%np)
142 real(rp) :: del_flux(elem%nfptot,lmesh%ne,prgvar_num)
143 real(rp) :: del_flux_hyd(elem%nfptot,lmesh%ne,2)
144 real(rp) :: rhot_(elem%np)
145 real(rp) :: rdens_(elem%np), u_(elem%np), v_(elem%np), w_(elem%np), wt_(elem%np), drho(elem%np)
146
147 real(rp) :: g11(elem%np), g12(elem%np), g22(elem%np)
148 real(rp) :: gsqrtv(elem%np), rgsqrtv(elem%np)
149 real(rp) :: x2d(elem2d%np,lmesh2d%ne), y2d(elem2d%np,lmesh2d%ne)
150 real(rp) :: x(elem%np), y(elem%np), twoovdel2(elem%np)
151 real(rp) :: cori(elem%np,2)
152 logical :: is_panel1to4
153 real(rp) :: s
154
155 integer :: ke, ke2d
156 integer :: p, p12, p3
157
158 real(rp) :: gamm, rgamm
159 real(rp) :: rp0
160 real(rp) :: rovp0, p0ovr
161 !------------------------------------------------------------------------
162
163 call prof_rapstart('cal_dyn_tend_bndflux', 3)
164 call get_ebnd_flux( &
165 del_flux, del_flux_hyd, & ! (out)
166 ddens_, momx_, momy_, momz_, drhot_, dpres_, dens_hyd, pres_hyd, & ! (in)
167 rtot, cvtot, cptot, & ! (in)
168 lmesh%Gsqrt, lmesh%GIJ(:,:,1,1), lmesh%GIJ(:,:,1,2), lmesh%GIJ(:,:,2,2), & ! (in)
169 lmesh%GsqrtH, lmesh%gam, lmesh%GI3(:,:,1), lmesh%GI3(:,:,2), & ! (in)
170 lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), & ! (in)
171 lmesh%vmapM, lmesh%vmapP, elem%IndexH2Dto3D_bnd, & ! (in)
172 lmesh, elem, lmesh2d, elem2d ) ! (in)
173 call prof_rapend('cal_dyn_tend_bndflux', 3)
174
175 !-----
176 call prof_rapstart('cal_dyn_tend_interior', 3)
177 gamm = cpdry / cvdry
178 rgamm = cvdry / cpdry
179 rp0 = 1.0_rp / pres00
180 rovp0 = rdry * rp0
181 p0ovr = pres00 / rdry
182
183 s = 1.0_rp
184 is_panel1to4 = .true.
185 if ( lmesh%panelID == 5 ) then
186 is_panel1to4 = .false.
187 else if ( lmesh%panelID == 6 ) then
188 is_panel1to4 = .false.
189 s = - 1.0_rp
190 end if
191
192 !$omp parallel private( &
193 !$omp RHOT_, rdens_, u_, v_, w_, wt_, &
194 !$omp Fx, Fy, Fz, LiftDelFlx, &
195 !$omp drho, DPRES_hyd, GradPhyd_x, GradPhyd_y, &
196 !$omp G11, G12, G22, GsqrtV, RGsqrtV, &
197 !$omp X, Y, twoOVdel2, &
198 !$omp CORI, ke, ke2D )
199
200 !$omp do
201 do ke2d = lmesh2d%NeS, lmesh2d%NeE
202 x2d(:,ke2d) = tan(lmesh2d%pos_en(:,ke2d,1))
203 y2d(:,ke2d) = tan(lmesh2d%pos_en(:,ke2d,2))
204 end do
205
206 !$omp do
207 do ke = lmesh%NeS, lmesh%NeE
208 !--
209 ke2d = lmesh%EMap3Dto2D(ke)
210 g11(:) = lmesh%GIJ(elem%IndexH2Dto3D,ke2d,1,1)
211 g12(:) = lmesh%GIJ(elem%IndexH2Dto3D,ke2d,1,2)
212 g22(:) = lmesh%GIJ(elem%IndexH2Dto3D,ke2d,2,2)
213 gsqrtv(:) = lmesh%Gsqrt(:,ke) / lmesh%GsqrtH(elem%IndexH2Dto3D,ke2d)
214 rgsqrtv(:) = 1.0_rp / gsqrtv(:)
215
216 !--
217 rhot_(:) = p0ovr * ( pres_hyd(:,ke) * rp0 )**rgamm + drhot_(:,ke)
218 ! DPRES_(:) = PRES00 * ( Rtot(:,ke) * rP0 * RHOT_(:) )**( CPtot(:,ke) / CVtot(:,ke) ) &
219 ! - PRES_hyd(:,ke)
220
221 rdens_(:) = 1.0_rp / ( ddens_(:,ke) + dens_hyd(:,ke) )
222 u_(:) = momx_(:,ke) * rdens_(:)
223 v_(:) = momy_(:,ke) * rdens_(:)
224 w_(:) = momz_(:,ke) * rdens_(:)
225 wt_(:) = w_(:) * rgsqrtv(:) + lmesh%GI3(:,ke,1) * u_(:) + lmesh%GI3(:,ke,2) * v_(:)
226
227 x(:) = x2d(elem%IndexH2Dto3D,ke2d)
228 y(:) = y2d(elem%IndexH2Dto3D,ke2d)
229 twoovdel2(:) = 2.0_rp / ( 1.0_rp + x(:)**2 + y(:)**2 )
230
231 cori(:,1) = s * ohm * twoovdel2(:) * ( - x(:) * y(:) * momx_(:,ke) + ( 1.0_rp + y(:)**2 ) * momy_(:,ke) )
232 cori(:,2) = s * ohm * twoovdel2(:) * ( - ( 1.0_rp + x(:)**2 ) * momx_(:,ke) + x(:) * y(:) * momy_(:,ke) )
233 if ( is_panel1to4 ) then
234 cori(:,1) = s * y(:) * cori(:,1)
235 cori(:,2) = s * y(:) * cori(:,2)
236 end if
237
238 drho(:) = matmul(intrpmat_vpordm1, ddens_(:,ke))
239
240 !-- Gradient hydrostatic pressure
241
242 dpres_hyd(:) = pres_hyd(:,ke) - pres_hyd_ref(:,ke)
243
244 call sparsemat_matmul(dx, gsqrtv(:) * dpres_hyd(:), fx)
245 call sparsemat_matmul(dz, gsqrtv(:) * lmesh%GI3(:,ke,1) * dpres_hyd(:), fz)
246 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux_hyd(:,ke,1), liftdelflx)
247 gradphyd_x(:) = lmesh%Escale(:,ke,1,1) * fx(:) &
248 + lmesh%Escale(:,ke,3,3) * fz(:) &
249 + liftdelflx(:)
250
251 call sparsemat_matmul(dy, gsqrtv(:) * dpres_hyd(:), fy)
252 call sparsemat_matmul(dz, gsqrtv(:) * lmesh%GI3(:,ke,2) * dpres_hyd(:), fz)
253 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux_hyd(:,ke,2), liftdelflx)
254 gradphyd_y(:) = lmesh%Escale(:,ke,2,2) * fy(:) &
255 + lmesh%Escale(:,ke,3,3) * fz(:) &
256 + liftdelflx(:)
257
258 !-- DENS
259 call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * momx_(:,ke), fx)
260 call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * momy_(:,ke), fy)
261 call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( ddens_(:,ke) + dens_hyd(:,ke) ) * wt_(:), fz)
262 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,dens_vid), liftdelflx)
263
264 dens_dt(:,ke) = - ( &
265 lmesh%Escale(:,ke,1,1) * fx(:) &
266 + lmesh%Escale(:,ke,2,2) * fy(:) &
267 + lmesh%Escale(:,ke,3,3) * fz(:) &
268 + liftdelflx(:) ) / lmesh%Gsqrt(:,ke)
269
270 !-- MOMX
271 call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * ( u_(:) * momx_(:,ke) + g11(:) * dpres_(:,ke) ), fx)
272 call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * ( v_(:) * momx_(:,ke) + g12(:) * dpres_(:,ke) ), fy)
273 call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( wt_(:) * momx_(:,ke) &
274 + ( lmesh%GI3(:,ke,1) * g11(:) + lmesh%GI3(:,ke,2) * g12(:) ) * dpres_(:,ke) ), fz)
275 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momx_vid), liftdelflx)
276
277 momx_dt(:,ke) = &
278 - ( lmesh%Escale(:,ke,1,1) * fx(:) &
279 + lmesh%Escale(:,ke,2,2) * fy(:) &
280 + lmesh%Escale(:,ke,3,3) * fz(:) &
281 + liftdelflx(:) ) / lmesh%Gsqrt(:,ke) &
282 - twoovdel2(:) * y(:) * &
283 ( x(:) * y(:) * u_(:) - (1.0_rp + y(:)**2) * v_(:) ) * momx_(:,ke) &
284 - ( g11(:) * gradphyd_x(:) + g12(:) * gradphyd_y(:) ) * rgsqrtv(:) &
285 + cori(:,1)
286
287 !-- MOMY
288 call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * ( u_(:) * momy_(:,ke) + g12(:) * dpres_(:,ke) ), fx)
289 call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * ( v_(:) * momy_(:,ke) + g22(:) * dpres_(:,ke) ), fy)
290 call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( wt_(:) * momy_(:,ke) &
291 + ( lmesh%GI3(:,ke,1) * g12(:) + lmesh%GI3(:,ke,2) * g22(:) ) * dpres_(:,ke) ), fz)
292 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momy_vid), liftdelflx)
293
294 momy_dt(:,ke) = &
295 - ( lmesh%Escale(:,ke,1,1) * fx(:) &
296 + lmesh%Escale(:,ke,2,2) * fy(:) &
297 + lmesh%Escale(:,ke,3,3) * fz(:) &
298 + liftdelflx(:) ) / lmesh%Gsqrt(:,ke) &
299 - twoovdel2(:) * x(:) * &
300 ( - (1.0_rp + x(:)**2) * u_(:) + x(:) * y(:) * v_(:) ) * momy_(:,ke) &
301 - ( g12(:) * gradphyd_x(:) + g22(:) * gradphyd_y(:) ) * rgsqrtv(:) &
302 + cori(:,2)
303
304 !-- MOMZ
305 call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * u_(:) * momz_(:,ke), fx)
306 call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * v_(:) * momz_(:,ke), fy)
307 call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( wt_(:) * momz_(:,ke) + rgsqrtv(:) * dpres_(:,ke) ), fz)
308 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momz_vid), liftdelflx)
309
310 momz_dt(:,ke) = &
311 - ( lmesh%Escale(:,ke,1,1) * fx(:) &
312 + lmesh%Escale(:,ke,2,2) * fy(:) &
313 + lmesh%Escale(:,ke,3,3) * fz(:) &
314 + liftdelflx(:) ) / lmesh%Gsqrt(:,ke) &
315 - grav * drho(:)
316
317 !-- RHOT
318 call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * u_(:) * rhot_(:), fx)
319 call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * v_(:) * rhot_(:), fy)
320 call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * wt_(:) * rhot_(:), fz)
321 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,rhot_vid), liftdelflx)
322
323 rhot_dt(:,ke) = &
324 - ( lmesh%Escale(:,ke,1,1) * fx(:) &
325 + lmesh%Escale(:,ke,2,2) * fy(:) &
326 + lmesh%Escale(:,ke,3,3) * fz(:) &
327 + liftdelflx(:) ) / lmesh%Gsqrt(:,ke)
328
329 end do
330 !$omp end do
331 !$omp end parallel
332 call prof_rapend('cal_dyn_tend_interior', 3)
333
334 return
336
337 !OCL SERIAL
339 DENS_dt, MOMX_dt, MOMY_dt, MOMZ_dt, RHOT_dt, & ! (out)
340 ddens_, momx_, momy_, momz_, drhot_, dpres_, & ! (in)
341 dens_hyd, pres_hyd, pres_hyd_ref, therm_hyd, & ! (in)
342 coriolis, rtot, cvtot, cptot, dphyddx, dphyddy, & ! (in)
343 element3d_operation, dx, dy, dz, sx, sy, sz, lift, & ! (in)
344 lmesh, elem, lmesh2d, elem2d ) ! (in)
345
348 use scale_const, only: &
349 ohm => const_ohm
350 implicit none
351
352 class(localmesh3d), intent(in) :: lmesh
353 class(elementbase3d), intent(in) :: elem
354 class(localmesh2d), intent(in) :: lmesh2d
355 class(elementbase2d), intent(in) :: elem2d
356 class(elementoperationbase3d), intent(in) :: element3d_operation
357 type(sparsemat), intent(in) :: dx, dy, dz, sx, sy, sz, lift
358 real(rp), intent(out) :: dens_dt(elem%np,lmesh%nea)
359 real(rp), intent(out) :: momx_dt(elem%np,lmesh%nea)
360 real(rp), intent(out) :: momy_dt(elem%np,lmesh%nea)
361 real(rp), intent(out) :: momz_dt(elem%np,lmesh%nea)
362 real(rp), intent(out) :: rhot_dt(elem%np,lmesh%nea)
363 real(rp), intent(in) :: ddens_(elem%np,lmesh%nea)
364 real(rp), intent(in) :: momx_(elem%np,lmesh%nea)
365 real(rp), intent(in) :: momy_(elem%np,lmesh%nea)
366 real(rp), intent(in) :: momz_(elem%np,lmesh%nea)
367 real(rp), intent(in) :: drhot_(elem%np,lmesh%nea)
368 real(rp), intent(in) :: dpres_(elem%np,lmesh%nea)
369 real(rp), intent(in) :: dens_hyd(elem%np,lmesh%nea)
370 real(rp), intent(in) :: pres_hyd(elem%np,lmesh%nea)
371 real(rp), intent(in) :: pres_hyd_ref(elem%np,lmesh%nea)
372 real(rp), intent(in) :: therm_hyd(elem%np,lmesh%nea)
373 real(rp), intent(in) :: coriolis(elem2d%np,lmesh2d%nea)
374 real(rp), intent(in) :: rtot (elem%np,lmesh%nea)
375 real(rp), intent(in) :: cvtot(elem%np,lmesh%nea)
376 real(rp), intent(in) :: cptot(elem%np,lmesh%nea)
377 real(rp), intent(in) :: dphyddx(elem%np,lmesh%nea)
378 real(rp), intent(in) :: dphyddy(elem%np,lmesh%nea)
379
380 real(rp) :: flux(elem%np,3,5), dflux(elem%np,4,5)
381 real(rp) :: del_flux(elem%nfptot,prgvar_num,lmesh%ne)
382 real(rp) :: u_, v_, w_, pt_
383 real(rp) :: drho(elem%np)
384 real(rp) :: rdens_(elem%np), gsqrtv(elem%np), rgsqrtv(elem%np), rgsqrt(elem%np)
385 real(rp) :: gsqrt_, gsqrtdpres_, e11, e22, e33
386
387 real(rp) :: g11(elem%np), g12(elem%np), g22(elem%np)
388 real(rp) :: x2d(elem2d%np,lmesh2d%ne), y2d(elem2d%np,lmesh2d%ne)
389 real(rp) :: x(elem%np), y(elem%np), twoovdel2(elem%np)
390 real(rp) :: cori(elem%np,2)
391 logical :: is_panel1to4
392 real(rp) :: s
393
394 integer :: ke, ke2d
395 integer :: p
396
397 real(rp) :: gamm, rgamm
398 real(rp) :: rp0
399 real(rp) :: rovp0, p0ovr
400 !------------------------------------------------------------------------
401
402 call prof_rapstart('cal_dyn_tend_bndflux', 3)
403 call get_ebnd_flux( &
404 del_flux, & ! (out)
405 ddens_, momx_, momy_, momz_, drhot_, dpres_, & ! (in)
406 dens_hyd, pres_hyd, therm_hyd, rtot, cvtot, cptot, & ! (in)
407 lmesh%Gsqrt, lmesh%GIJ(:,:,1,1), lmesh%GIJ(:,:,1,2), lmesh%GIJ(:,:,2,2), & ! (in)
408 lmesh%GsqrtH, lmesh%gam, lmesh%GI3(:,:,1), lmesh%GI3(:,:,2), & ! (in)
409 lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), & ! (in)
410 lmesh%vmapM, lmesh%vmapP, elem%IndexH2Dto3D_bnd, & ! (in)
411 lmesh, elem, lmesh2d, elem2d ) ! (in)
412 call prof_rapend('cal_dyn_tend_bndflux', 3)
413
414 !-----
415 call prof_rapstart('cal_dyn_tend_interior', 3)
416 gamm = cpdry / cvdry
417 rgamm = cvdry / cpdry
418 rp0 = 1.0_rp / pres00
419 rovp0 = rdry * rp0
420 p0ovr = pres00 / rdry
421
422 s = 1.0_rp
423 is_panel1to4 = .true.
424 if ( lmesh%panelID == 5 ) then
425 is_panel1to4 = .false.
426 else if ( lmesh%panelID == 6 ) then
427 is_panel1to4 = .false.
428 s = - 1.0_rp
429 end if
430
431 !$omp parallel private( ke, ke2d, p, &
432 !$omp u_, v_, w_, pt_, &
433 !$omp drho, CORI, &
434 !$omp RDENS_, RGsqrt, GsqrtV, RGsqrtV, Gsqrt_, GsqrtDPRES_, &
435 !$omp G11, G12, G22, X, Y, twoOVdel2, &
436 !$omp E11, E22, E33, DFlux, Flux )
437
438 !$omp do
439 do ke2d = lmesh2d%NeS, lmesh2d%NeE
440 x2d(:,ke2d) = tan(lmesh2d%pos_en(:,ke2d,1))
441 y2d(:,ke2d) = tan(lmesh2d%pos_en(:,ke2d,2))
442 end do
443
444 !$omp do
445 do ke = lmesh%NeS, lmesh%NeE
446 !--
447 ke2d = lmesh%EMap3Dto2D(ke)
448
449 do p=1, elem%Np
450 g11(p) = lmesh%GIJ(elem%IndexH2Dto3D(p),ke2d,1,1)
451 g12(p) = lmesh%GIJ(elem%IndexH2Dto3D(p),ke2d,1,2)
452 g22(p) = lmesh%GIJ(elem%IndexH2Dto3D(p),ke2d,2,2)
453 end do
454 do p=1, elem%Np
455 gsqrtv(p) = lmesh%Gsqrt(p,ke) / lmesh%GsqrtH(elem%IndexH2Dto3D(p),ke2d)
456 rgsqrtv(p) = 1.0_rp / gsqrtv(p)
457 rgsqrt(p) = 1.0_rp / lmesh%Gsqrt(p,ke)
458 rdens_(p) = 1.0_rp / ( ddens_(p,ke) + dens_hyd(p,ke) )
459 end do
460
461 !--
462
463 do p=1, elem%Np
464 gsqrt_ = lmesh%Gsqrt(p,ke)
465 flux(p,1,dens_vid) = gsqrt_ * momx_(p,ke)
466 flux(p,2,dens_vid) = gsqrt_ * momy_(p,ke)
467 flux(p,3,dens_vid) = gsqrt_ * ( &
468 momz_(p,ke) * rgsqrtv(p) &
469 + lmesh%GI3(p,ke,1) * momx_(p,ke) &
470 + lmesh%GI3(p,ke,2) * momy_(p,ke) )
471 end do
472 do p=1, elem%Np
473 pt_ = ( therm_hyd(p,ke) + drhot_(p,ke) ) * rdens_(p)
474
475 flux(p,1,rhot_vid) = flux(p,1,dens_vid) * pt_
476 flux(p,2,rhot_vid) = flux(p,2,dens_vid) * pt_
477 flux(p,3,rhot_vid) = flux(p,3,dens_vid) * pt_
478
479 w_ = momz_(p,ke) * rdens_(p)
480 flux(p,1,momz_vid) = flux(p,1,dens_vid) * w_
481 flux(p,2,momz_vid) = flux(p,2,dens_vid) * w_
482 flux(p,3,momz_vid) = flux(p,3,dens_vid) * w_ + lmesh%Gsqrt(p,ke) * rgsqrtv(p) * dpres_(p,ke)
483 end do
484 do p=1, elem%Np
485 gsqrtdpres_ = lmesh%Gsqrt(p,ke) * dpres_(p,ke)
486
487 u_ = momx_(p,ke) * rdens_(p)
488 v_ = momy_(p,ke) * rdens_(p)
489
490 flux(p,1,momx_vid) = flux(p,1,dens_vid) * u_ + g11(p) * gsqrtdpres_
491 flux(p,2,momx_vid) = flux(p,2,dens_vid) * u_ + g12(p) * gsqrtdpres_
492 flux(p,3,momx_vid) = flux(p,3,dens_vid) * u_ + gsqrtdpres_ * ( g11(p) * lmesh%GI3(p,ke,1) + g12(p) * lmesh%GI3(p,ke,2) )
493
494 flux(p,1,momy_vid) = flux(p,1,dens_vid) * v_ + g12(p) * gsqrtdpres_
495 flux(p,2,momy_vid) = flux(p,2,dens_vid) * v_ + g22(p) * gsqrtdpres_
496 flux(p,3,momy_vid) = flux(p,3,dens_vid) * v_ + gsqrtdpres_ * ( g12(p) * lmesh%GI3(p,ke,1) + g22(p) * lmesh%GI3(p,ke,2) )
497 end do
498
499 call element3d_operation%Div_var5( &
500 flux, del_flux(:,:,ke), & ! (in)
501 dflux ) ! (out)
502
503
504 do p=1, elem%Np
505 e11 = lmesh%Escale(p,ke,1,1)
506 e22 = lmesh%Escale(p,ke,2,2)
507 e33 = lmesh%Escale(p,ke,3,3)
508
509 dens_dt(p,ke) = - ( &
510 e11 * dflux(p,1,dens_vid) &
511 + e22 * dflux(p,2,dens_vid) &
512 + e33 * dflux(p,3,dens_vid) &
513 + dflux(p,4,dens_vid) ) * rgsqrt(p)
514
515 rhot_dt(p,ke) = - ( &
516 e11 * dflux(p,1,rhot_vid) &
517 + e22 * dflux(p,2,rhot_vid) &
518 + e33 * dflux(p,3,rhot_vid) &
519 + dflux(p,4,rhot_vid) ) * rgsqrt(p)
520 end do
521
522 call element3d_operation%VFilterPM1( ddens_(:,ke), & ! (in)
523 drho ) ! (out)
524
525 do p=1, elem%Np
526 e11 = lmesh%Escale(p,ke,1,1)
527 e22 = lmesh%Escale(p,ke,2,2)
528 e33 = lmesh%Escale(p,ke,3,3)
529
530 momz_dt(p,ke) = - ( &
531 e11 * dflux(p,1,momz_vid) &
532 + e22 * dflux(p,2,momz_vid) &
533 + e33 * dflux(p,3,momz_vid) &
534 + dflux(p,4,momz_vid) ) * rgsqrt(p) &
535 - grav * drho(p)
536 end do
537
538 !--
539 do p=1, elem%Np
540 x(p) = x2d(elem%IndexH2Dto3D(p),ke2d)
541 y(p) = y2d(elem%IndexH2Dto3D(p),ke2d)
542 twoovdel2(p) = 2.0_rp / ( 1.0_rp + x(p)**2 + y(p)**2 )
543 end do
544
545 do p=1, elem%Np
546 cori(p,1) = s * ohm * twoovdel2(p) * ( - x(p) * y(p) * momx_(p,ke) + ( 1.0_rp + y(p)**2 ) * momy_(p,ke) )
547 cori(p,2) = s * ohm * twoovdel2(p) * ( - ( 1.0_rp + x(p)**2 ) * momx_(p,ke) + x(p) * y(p) * momy_(p,ke) )
548 end do
549 if ( is_panel1to4 ) then
550 do p=1, elem%Np
551 cori(p,1) = s * y(p) * cori(p,1)
552 cori(p,2) = s * y(p) * cori(p,2)
553 end do
554 end if
555
556 do p=1, elem%Np
557 u_ = momx_(p,ke) * rdens_(p)
558 v_ = momy_(p,ke) * rdens_(p)
559
560 momx_dt(p,ke) = - ( g11(p) * dphyddx(p,ke) + g12(p) * dphyddy(p,ke) ) &
561 - twoovdel2(p) * y(p) * &
562 ( x(p) * y(p) * u_ - (1.0_rp + y(p)**2) * v_ ) * momx_(p,ke) &
563 + cori(p,1)
564
565 momy_dt(p,ke) = - ( g12(p) * dphyddx(p,ke) + g22(p) * dphyddy(p,ke) ) &
566 - twoovdel2(p) * x(p) * &
567 ( - (1.0_rp + x(p)**2) * u_ + x(p) * y(p) * v_ ) * momy_(p,ke) &
568 + cori(p,2)
569 end do
570
571 do p=1, elem%Np
572 e11 = lmesh%Escale(p,ke,1,1)
573 e22 = lmesh%Escale(p,ke,2,2)
574 e33 = lmesh%Escale(p,ke,3,3)
575
576 momx_dt(p,ke) = momx_dt(p,ke) - ( &
577 e11 * dflux(p,1,momx_vid) &
578 + e22 * dflux(p,2,momx_vid) &
579 + e33 * dflux(p,3,momx_vid) &
580 + dflux(p,4,momx_vid) ) * rgsqrt(p)
581
582 momy_dt(p,ke) = momy_dt(p,ke) - ( &
583 e11 * dflux(p,1,momy_vid) &
584 + e22 * dflux(p,2,momy_vid) &
585 + e33 * dflux(p,3,momy_vid) &
586 + dflux(p,4,momy_vid) ) * rgsqrt(p)
587 end do
588 end do
589 !$omp end do
590 !$omp end parallel
591 call prof_rapend('cal_dyn_tend_interior', 3)
592
593 return
595
596!OCL SERIAL
598 DENS_dt, MOMX_dt, MOMY_dt, MOMZ_dt, RHOT_dt, & ! (out)
599 ddens_, momx_, momy_, momz_, drhot_, dpres_, & ! (in)
600 dens_hyd, pres_hyd, pres_hyd_ref, therm_hyd, & ! (in)
601 coriolis, rtot, cvtot, cptot, dphyddx, dphyddy, & ! (in)
602 element3d_operation, dx, dy, dz, sx, sy, sz, lift, & ! (in)
603 lmesh, elem, lmesh2d, elem2d ) ! (in)
604
607 use scale_const, only: &
608 ohm => const_ohm
609 implicit none
610
611 class(localmesh3d), intent(in) :: lmesh
612 class(elementbase3d), intent(in) :: elem
613 class(localmesh2d), intent(in) :: lmesh2d
614 class(elementbase2d), intent(in) :: elem2d
615 class(elementoperationbase3d), intent(in) :: element3d_operation
616 type(sparsemat), intent(in) :: dx, dy, dz, sx, sy, sz, lift
617 real(rp), intent(out) :: dens_dt(elem%np,lmesh%nea)
618 real(rp), intent(out) :: momx_dt(elem%np,lmesh%nea)
619 real(rp), intent(out) :: momy_dt(elem%np,lmesh%nea)
620 real(rp), intent(out) :: momz_dt(elem%np,lmesh%nea)
621 real(rp), intent(out) :: rhot_dt(elem%np,lmesh%nea)
622 real(rp), intent(in) :: ddens_(elem%np,lmesh%nea)
623 real(rp), intent(in) :: momx_(elem%np,lmesh%nea)
624 real(rp), intent(in) :: momy_(elem%np,lmesh%nea)
625 real(rp), intent(in) :: momz_(elem%np,lmesh%nea)
626 real(rp), intent(in) :: drhot_(elem%np,lmesh%nea)
627 real(rp), intent(in) :: dpres_(elem%np,lmesh%nea)
628 real(rp), intent(in) :: dens_hyd(elem%np,lmesh%nea)
629 real(rp), intent(in) :: pres_hyd(elem%np,lmesh%nea)
630 real(rp), intent(in) :: pres_hyd_ref(elem%np,lmesh%nea)
631 real(rp), intent(in) :: therm_hyd(elem%np,lmesh%nea)
632 real(rp), intent(in) :: coriolis(elem2d%np,lmesh2d%nea)
633 real(rp), intent(in) :: rtot (elem%np,lmesh%nea)
634 real(rp), intent(in) :: cvtot(elem%np,lmesh%nea)
635 real(rp), intent(in) :: cptot(elem%np,lmesh%nea)
636 real(rp), intent(in) :: dphyddx(elem%np,lmesh%nea)
637 real(rp), intent(in) :: dphyddy(elem%np,lmesh%nea)
638
639 real(rp) :: fx(elem%np), fy(elem%np), fz(elem%np), liftdelflx(elem%np)
640 real(rp) :: dpres_hyd(elem%np), gradphyd_x(elem%np), gradphyd_y(elem%np)
641 real(rp) :: del_flux(elem%nfptot,lmesh%ne,prgvar_num)
642 real(rp) :: del_flux_hyd(elem%nfptot,lmesh%ne,2)
643 real(rp) :: rhot_(elem%np)
644 real(rp) :: rdens_(elem%np), u_(elem%np), v_(elem%np), w_(elem%np), wt_(elem%np), drho(elem%np)
645
646 real(rp) :: g11(elem%np), g12(elem%np), g22(elem%np)
647 real(rp) :: gsqrtv(elem%np), rgsqrtv(elem%np), rgam2(elem%np)
648 real(rp) :: x2d(elem2d%np,lmesh2d%ne), y2d(elem2d%np,lmesh2d%ne)
649 real(rp) :: x(elem%np), y(elem%np), twoovdel2(elem%np)
650 real(rp) :: om1(elem%np), om2(elem%np), om3(elem%np), del(elem%np), r(elem%np)
651 logical :: is_panel1to4
652 real(rp) :: s
653
654 integer :: ke, ke2d
655 integer :: p
656
657 real(rp) :: rgamm
658 real(rp) :: rp0
659 real(rp) :: p0ovr
660 !------------------------------------------------------------------------
661
662 call prof_rapstart('cal_dyn_tend_bndflux', 3)
663 call get_ebnd_flux( &
664 del_flux, del_flux_hyd, & ! (out)
665 ddens_, momx_, momy_, momz_, drhot_, dpres_, dens_hyd, pres_hyd, & ! (in)
666 rtot, cvtot, cptot, & ! (in)
667 lmesh%Gsqrt, lmesh%GIJ(:,:,1,1), lmesh%GIJ(:,:,1,2), lmesh%GIJ(:,:,2,2), & ! (in)
668 lmesh%GsqrtH, lmesh%gam, lmesh%GI3(:,:,1), lmesh%GI3(:,:,2), & ! (in)
669 lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), & ! (in)
670 lmesh%vmapM, lmesh%vmapP, elem%IndexH2Dto3D_bnd, & ! (in)
671 lmesh, elem, lmesh2d, elem2d ) ! (in)
672 call prof_rapend('cal_dyn_tend_bndflux', 3)
673
674 !-----
675 call prof_rapstart('cal_dyn_tend_interior', 3)
676 rgamm = cvdry / cpdry
677 rp0 = 1.0_rp / pres00
678 p0ovr = pres00 / rdry
679
680 s = 2.0_rp * ohm
681 is_panel1to4 = .true.
682 if ( lmesh%panelID == 5 ) then
683 is_panel1to4 = .false.
684 else if ( lmesh%panelID == 6 ) then
685 is_panel1to4 = .false.
686 s = - s
687 end if
688
689 !$omp parallel private( &
690 !$omp RHOT_, rdens_, u_, v_, w_, wt_, &
691 !$omp Fx, Fy, Fz, LiftDelFlx, &
692 !$omp drho, DPRES_hyd, GradPhyd_x, GradPhyd_y, &
693 !$omp G11, G12, G22, Rgam2, GsqrtV, RGsqrtV, &
694 !$omp X, Y, twoOVdel2, &
695 !$omp OM1, OM2, OM3, DEL, R, ke, ke2D )
696
697 !$omp do
698 do ke2d = lmesh2d%NeS, lmesh2d%NeE
699 x2d(:,ke2d) = tan(lmesh2d%pos_en(:,ke2d,1))
700 y2d(:,ke2d) = tan(lmesh2d%pos_en(:,ke2d,2))
701 end do
702
703 !$omp do
704 do ke = lmesh%NeS, lmesh%NeE
705 !--
706 ke2d = lmesh%EMap3Dto2D(ke)
707 rgam2(:) = 1.0_rp / lmesh%gam(:,ke)**2
708 g11(:) = lmesh%GIJ(elem%IndexH2Dto3D,ke2d,1,1) * rgam2(:)
709 g12(:) = lmesh%GIJ(elem%IndexH2Dto3D,ke2d,1,2) * rgam2(:)
710 g22(:) = lmesh%GIJ(elem%IndexH2Dto3D,ke2d,2,2) * rgam2(:)
711 gsqrtv(:) = lmesh%Gsqrt(:,ke) * rgam2(:) / lmesh%GsqrtH(elem%IndexH2Dto3D,ke2d)
712 rgsqrtv(:) = 1.0_rp / gsqrtv(:)
713
714 !--
715 rhot_(:) = p0ovr * ( pres_hyd(:,ke) * rp0 )**rgamm + drhot_(:,ke)
716 ! DPRES_(:) = PRES00 * ( Rtot(:,ke) * rP0 * RHOT_(:) )**( CPtot(:,ke) / CVtot(:,ke) ) &
717 ! - PRES_hyd(:,ke)
718
719 rdens_(:) = 1.0_rp / ( ddens_(:,ke) + dens_hyd(:,ke) )
720 u_(:) = momx_(:,ke) * rdens_(:)
721 v_(:) = momy_(:,ke) * rdens_(:)
722 w_(:) = momz_(:,ke) * rdens_(:)
723 wt_(:) = w_(:) * rgsqrtv(:) + lmesh%GI3(:,ke,1) * u_(:) + lmesh%GI3(:,ke,2) * v_(:)
724
725 x(:) = x2d(elem%IndexH2Dto3D,ke2d)
726 y(:) = y2d(elem%IndexH2Dto3D,ke2d)
727 del(:) = sqrt( 1.0_rp + x(:)**2 + y(:)**2 )
728 twoovdel2(:) = 2.0_rp / ( 1.0_rp + x(:)**2 + y(:)**2 )
729
730 r(:) = rplanet * lmesh%gam(:,ke)
731
732 !- pnl=1~4: OM1: 0, OM2: s del / (r (1+Y^2)) , s OM3 : Y /del
733 !- pnl=5,6: OM1: - s X del/(r (1+X^2)), OM2: - s Y del/(r(1+Y^2)), OM3 : s/del
734 if ( is_panel1to4 ) then
735 om1(:) = 0.0_rp
736 om2(:) = s * del(:) / ( r(:) * ( 1.0_rp + y(:)**2 ) )
737 om3(:) = s * y(:) / del(:)
738 else
739 om1(:) = - s * x(:) * del(:) / ( r(:) * ( 1.0_rp + x(:)**2 ) )
740 om2(:) = - s * y(:) * del(:) / ( r(:) * ( 1.0_rp + y(:)**2 ) )
741 om3(:) = s / del(:)
742 end if
743
744 drho(:) = matmul(intrpmat_vpordm1, ddens_(:,ke))
745
746 !-- Gradient hydrostatic pressure
747
748 dpres_hyd(:) = pres_hyd(:,ke) - pres_hyd_ref(:,ke)
749
750 call sparsemat_matmul(dx, gsqrtv(:) * dpres_hyd(:), fx)
751 call sparsemat_matmul(dz, gsqrtv(:) * lmesh%GI3(:,ke,1) * dpres_hyd(:), fz)
752 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux_hyd(:,ke,1), liftdelflx)
753 gradphyd_x(:) = lmesh%Escale(:,ke,1,1) * fx(:) &
754 + lmesh%Escale(:,ke,3,3) * fz(:) &
755 + liftdelflx(:)
756
757 call sparsemat_matmul(dy, gsqrtv(:) * dpres_hyd(:), fy)
758 call sparsemat_matmul(dz, gsqrtv(:) * lmesh%GI3(:,ke,2) * dpres_hyd(:), fz)
759 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux_hyd(:,ke,2), liftdelflx)
760 gradphyd_y(:) = lmesh%Escale(:,ke,2,2) * fy(:) &
761 + lmesh%Escale(:,ke,3,3) * fz(:) &
762 + liftdelflx(:)
763
764 !-- DENS
765 call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * momx_(:,ke), fx)
766 call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * momy_(:,ke), fy)
767 call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( ddens_(:,ke) + dens_hyd(:,ke) ) * wt_(:), fz)
768 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,dens_vid), liftdelflx)
769
770 dens_dt(:,ke) = - ( &
771 lmesh%Escale(:,ke,1,1) * fx(:) &
772 + lmesh%Escale(:,ke,2,2) * fy(:) &
773 + lmesh%Escale(:,ke,3,3) * fz(:) &
774 + liftdelflx(:) ) / lmesh%Gsqrt(:,ke)
775
776 !-- MOMX
777 call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * ( u_(:) * momx_(:,ke) + g11(:) * dpres_(:,ke) ), fx)
778 call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * ( v_(:) * momx_(:,ke) + g12(:) * dpres_(:,ke) ), fy)
779 call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( wt_(:) * momx_(:,ke) &
780 + ( lmesh%GI3(:,ke,1) * g11(:) + lmesh%GI3(:,ke,2) * g12(:) ) * dpres_(:,ke) ), fz)
781 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momx_vid), liftdelflx)
782
783 momx_dt(:,ke) = &
784 - ( lmesh%Escale(:,ke,1,1) * fx(:) &
785 + lmesh%Escale(:,ke,2,2) * fy(:) &
786 + lmesh%Escale(:,ke,3,3) * fz(:) &
787 + liftdelflx(:) ) / lmesh%Gsqrt(:,ke) &
788 - twoovdel2(:) * y(:) * & !-> metric terms
789 ( x(:) * y(:) * u_(:) - ( 1.0_rp + y(:)**2 ) * v_(:) ) * momx_(:,ke) & !
790 - 2.0_rp * u_(:) * momz_(:,ke) / r(:) & !<-
791 - ( g11(:) * gradphyd_x(:) + g12(:) * gradphyd_y(:) ) * rgsqrtv(:) & !-> gradient hydrostatic pressure
792 - lmesh%Gsqrt(:,ke) * ( g11(:) * ( om2(:) * momz_(:,ke) - om3(:) * momy_(:,ke) ) & !-> Coriolis term
793 - g12(:) * ( om1(:) * momz_(:,ke) - om3(:) * momx_(:,ke) ) )
794
795 !-- MOMY
796 call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * ( u_(:) * momy_(:,ke) + g12(:) * dpres_(:,ke) ), fx)
797 call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * ( v_(:) * momy_(:,ke) + g22(:) * dpres_(:,ke) ), fy)
798 call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( wt_(:) * momy_(:,ke) &
799 + ( lmesh%GI3(:,ke,1) * g12(:) + lmesh%GI3(:,ke,2) * g22(:) ) * dpres_(:,ke) ), fz)
800 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momy_vid), liftdelflx)
801
802 momy_dt(:,ke) = &
803 - ( lmesh%Escale(:,ke,1,1) * fx(:) &
804 + lmesh%Escale(:,ke,2,2) * fy(:) &
805 + lmesh%Escale(:,ke,3,3) * fz(:) &
806 + liftdelflx(:) ) / lmesh%Gsqrt(:,ke) &
807 - twoovdel2(:) * x(:) * & !-> metric terms
808 ( - (1.0_rp + x(:)**2) * u_(:) + x(:) * y(:) * v_(:) ) * momy_(:,ke) & !
809 - 2.0_rp * v_(:) * momz_(:,ke) / r(:) & !<-
810 - ( g12(:) * gradphyd_x(:) + g22(:) * gradphyd_y(:) ) * rgsqrtv(:) & !-> gradient hydrostatic pressure
811 - lmesh%Gsqrt(:,ke) * ( g12(:) * ( om2(:) * momz_(:,ke) - om3(:) * momy_(:,ke) ) & !-> Coriolis term
812 - g22(:) * ( om1(:) * momz_(:,ke) - om3(:) * momx_(:,ke) ) )
813
814 !-- MOMZ
815 call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * u_(:) * momz_(:,ke), fx)
816 call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * v_(:) * momz_(:,ke), fy)
817 call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( wt_(:) * momz_(:,ke) + rgsqrtv(:) * dpres_(:,ke) ), fz)
818 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momz_vid), liftdelflx)
819
820 momz_dt(:,ke) = &
821 - ( lmesh%Escale(:,ke,1,1) * fx(:) &
822 + lmesh%Escale(:,ke,2,2) * fy(:) &
823 + lmesh%Escale(:,ke,3,3) * fz(:) &
824 + liftdelflx(:) ) / lmesh%Gsqrt(:,ke) &
825 - 0.25_rp * r(:) * twoovdel2(:)**2 * ( 1.0_rp * x(:)**2 ) * ( 1.0_rp * y(:)**2 ) & !-> metric terms
826 * ( - ( 1.0_rp + x(:)**2 ) * momx_(:,ke) * u_(:) & !
827 + 2.0_rp * x(:) * y(:) * momx_(:,ke) * v_(:) & !
828 - ( 1.0_rp + y(:)**2 ) * momy_(:,ke) * v_(:) ) & !<-
829 + 2.0_rp * dpres_(:,ke) / r(:) & !-> metric term with gradient of pressure deviaition
830 - lmesh%Gsqrt(:,ke) * ( om1(:) * momy_(:,ke) - om2(:) * momx_(:,ke) ) & !-> Coriolis term
831 - grav * rgam2(:) * drho(:) !-> buoyancy term
832
833 !-- RHOT
834 call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * u_(:) * rhot_(:), fx)
835 call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * v_(:) * rhot_(:), fy)
836 call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * wt_(:) * rhot_(:), fz)
837 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,rhot_vid), liftdelflx)
838
839 rhot_dt(:,ke) = &
840 - ( lmesh%Escale(:,ke,1,1) * fx(:) &
841 + lmesh%Escale(:,ke,2,2) * fy(:) &
842 + lmesh%Escale(:,ke,3,3) * fz(:) &
843 + liftdelflx(:) ) / lmesh%Gsqrt(:,ke)
844 end do
845 !$omp end do
846 !$omp end parallel
847 call prof_rapend('cal_dyn_tend_interior', 3)
848
849 return
851
module FElib / Fluid dyn solver / Atmosphere / Global nonhydrostatic model / HEVE
subroutine, public atm_dyn_dgm_globalnh3d_rhot_heve_cal_tend_shallow_atm_asis(dens_dt, momx_dt, momy_dt, momz_dt, rhot_dt, ddens_, momx_, momy_, momz_, drhot_, dpres_, dens_hyd, pres_hyd, pres_hyd_ref, therm_hyd, coriolis, rtot, cvtot, cptot, dphyddx, dphyddy, element3d_operation, dx, dy, dz, sx, sy, sz, lift, lmesh, elem, lmesh2d, elem2d)
subroutine, public atm_dyn_dgm_globalnonhydro3d_rhot_heve_cal_tend_deep_atm(dens_dt, momx_dt, momy_dt, momz_dt, rhot_dt, ddens_, momx_, momy_, momz_, drhot_, dpres_, dens_hyd, pres_hyd, pres_hyd_ref, therm_hyd, coriolis, rtot, cvtot, cptot, dphyddx, dphyddy, element3d_operation, dx, dy, dz, sx, sy, sz, lift, lmesh, elem, lmesh2d, elem2d)
subroutine, public atm_dyn_dgm_globalnonhydro3d_rhot_heve_cal_tend_shallow_atm(dens_dt, momx_dt, momy_dt, momz_dt, rhot_dt, ddens_, momx_, momy_, momz_, drhot_, dpres_, dens_hyd, pres_hyd, pres_hyd_ref, therm_hyd, coriolis, rtot, cvtot, cptot, dphyddx, dphyddy, element3d_operation, dx, dy, dz, sx, sy, sz, lift, lmesh, elem, lmesh2d, elem2d)
module FElib / Fluid dyn solver / Atmosphere / Nonhydrostatic model / Common
subroutine, public atm_dyn_dgm_nonhydro3d_common_init(mesh)
Initialize a common module for atmospheric nonhydrostatic dynamical core.
real(rp), dimension(:,:), allocatable, public intrpmat_vpordm1
subroutine, public atm_dyn_dgm_nonhydro3d_common_final()
Finalize a common module for atmospheric nonhydrostatic dynamical core.
module FElib / Fluid dyn solver / Atmosphere / Nonhydrostatic model / HEVE / Numflux
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 2D
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 2D domain)
Derived type to manage a computational mesh (base type for 3D domain)
Derived type representing a field with 3D mesh.
Derived type to manage a sparse matrix.