FE-Project
Loading...
Searching...
No Matches
scale_atm_dyn_dgm_nonhydro3d_hevi_gmres.F90
Go to the documentation of this file.
1!-------------------------------------------------------------------------------
2!> module FElib / Fluid dyn solver / Atmosphere / Regional nonhydrostatic model / HEVI
3!!
4!! @par Description
5!! HEVI DGM scheme for Atmospheric dynamical process in which GMRES is used for implicit solver.
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, Team SCALE
10!<
11!-------------------------------------------------------------------------------
12#include "scaleFElib.h"
13module scale_atm_dyn_dgm_nonhydro3d_hevi_gmres
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
29 use scale_gmres, only: gmres
30
31 use scale_sparsemat
32 use scale_element_base, only: &
33 elementbase2d, elementbase3d
34 use scale_element_hexahedral, only: hexahedralelement
35 use scale_localmesh_2d, only: localmesh2d
36 use scale_localmesh_3d, only: localmesh3d
37 use scale_mesh_base3d, only: meshbase3d
38 use scale_localmeshfield_base, only: localmeshfield3d
39 use scale_meshfield_base, only: meshfield3d
40
41
42 !-----------------------------------------------------------------------------
43 implicit none
44 private
45 !-----------------------------------------------------------------------------
46 !
47 !++ Public procedures
48 !
49 public :: atm_dyn_dgm_nonhydro3d_hevi_init
50 public :: atm_dyn_dgm_nonhydro3d_hevi_final
51 public :: atm_dyn_dgm_nonhydro3d_hevi_cal_tend
52 public :: atm_dyn_dgm_nonhydro3d_hevi_cal_grad_diffvars
53 public :: atm_dyn_dgm_nonhydro3d_hevi_cal_vi
54
55 !-----------------------------------------------------------------------------
56 !
57 !++ Public parameters & variables
58 !
59
60 !-----------------------------------------------------------------------------
61 !
62 !++ Private procedures & variables
63 !
64 !-------------------
65
66 integer, private, parameter :: DDENS_VID = 1
67 integer, private, parameter :: MOMX_VID = 2
68 integer, private, parameter :: MOMY_VID = 3
69 integer, private, parameter :: MOMZ_VID = 4
70 integer, private, parameter :: DRHOT_VID = 5
71 integer, private, parameter :: PROG_VARS_NUM = 5
72
73 integer, private, parameter :: VARS_GxU_ID = 1
74 integer, private, parameter :: VARS_GyU_ID = 2
75 integer, private, parameter :: VARS_GzU_ID = 3
76 integer, private, parameter :: VARS_GxV_ID = 4
77 integer, private, parameter :: VARS_GyV_ID = 5
78 integer, private, parameter :: VARS_GzV_ID = 6
79 integer, private, parameter :: VARS_GxW_ID = 7
80 integer, private, parameter :: VARS_GyW_ID = 8
81 integer, private, parameter :: VARS_GzW_ID = 9
82 integer, private, parameter :: VARS_GxPT_ID = 10
83 integer, private, parameter :: VARS_GyPT_ID = 11
84 integer, private, parameter :: VARS_GzPT_ID = 12
85 integer, private, parameter :: AUX_DIFFVARS_NUM = 12
86
87 real(RP), private, allocatable :: IntrpMat_VPOrdM1(:,:)
88
89 private :: cal_del_flux_dyn
90 private :: cal_del_graddiffvar
91
92contains
93 subroutine atm_dyn_dgm_nonhydro3d_hevi_init( mesh )
94
95 implicit none
96 class(meshbase3d), intent(in) :: mesh
97
98 integer :: p1, p2, p3, p_
99 type(elementbase3d), pointer :: elem
100 real(rp) :: invv_pordm1(mesh%refelem3d%np,mesh%refelem3d%np)
101 !--------------------------------------------
102
103 elem => mesh%refElem3D
104 allocate( intrpmat_vpordm1(elem%Np,elem%Np) )
105
106 invv_pordm1(:,:) = elem%invV
107 do p2=1, elem%Nnode_h1D
108 do p1=1, elem%Nnode_h1D
109 p_ = p1 + (p2-1)*elem%Nnode_h1D + (elem%Nnode_v-1)*elem%Nnode_h1D**2
110 invv_pordm1(p_,:) = 0.0_rp
111 end do
112 end do
113 intrpmat_vpordm1(:,:) = matmul(elem%V, invv_pordm1)
114
115 return
116 end subroutine atm_dyn_dgm_nonhydro3d_hevi_init
117
118 subroutine atm_dyn_dgm_nonhydro3d_hevi_final()
119 implicit none
120 !--------------------------------------------
121
122 deallocate( intrpmat_vpordm1 )
123
124 return
125 end subroutine atm_dyn_dgm_nonhydro3d_hevi_final
126
127 !-------------------------------
128
129 subroutine atm_dyn_dgm_nonhydro3d_hevi_cal_tend( &
130 DENS_dt, MOMX_dt, MOMY_dt, MOMZ_dt, RHOT_dt, & ! (out)
131 ddens_, momx_, momy_, momz_, drhot_, dens_hyd, pres_hyd, coriolis, & ! (in)
132 gxu_, gyu_, gzu_, gxv_, gyv_, gzv_, gxw_, gyw_, gzw_, gxpt_, gypt_, gzpt_, & ! (in)
133 visccoef_h, visccoef_v, diffcoef_h, diffcoef_v, & ! (in)
134 dx, dy, dz, sx, sy, sz, lift, lmesh, elem, lmesh2d, elem2d )
135
136 implicit none
137
138 class(localmesh3d), intent(in) :: lmesh
139 class(elementbase3d), intent(in) :: elem
140 class(localmesh2d), intent(in) :: lmesh2d
141 class(elementbase2d), intent(in) :: elem2d
142 type(sparsemat), intent(in) :: dx, dy, dz, sx, sy, sz, lift
143 real(rp), intent(out) :: dens_dt(elem%np,lmesh%nea)
144 real(rp), intent(out) :: momx_dt(elem%np,lmesh%nea)
145 real(rp), intent(out) :: momy_dt(elem%np,lmesh%nea)
146 real(rp), intent(out) :: momz_dt(elem%np,lmesh%nea)
147 real(rp), intent(out) :: rhot_dt(elem%np,lmesh%nea)
148 real(rp), intent(in) :: ddens_(elem%np,lmesh%nea)
149 real(rp), intent(in) :: momx_(elem%np,lmesh%nea)
150 real(rp), intent(in) :: momy_(elem%np,lmesh%nea)
151 real(rp), intent(in) :: momz_(elem%np,lmesh%nea)
152 real(rp), intent(in) :: drhot_(elem%np,lmesh%nea)
153 real(rp), intent(in) :: dens_hyd(elem%np,lmesh%nea)
154 real(rp), intent(in) :: pres_hyd(elem%np,lmesh%nea)
155 real(rp), intent(in) :: coriolis(elem2d%np,lmesh2d%nea)
156 real(rp), intent(in) :: gxu_(elem%np,lmesh%nea)
157 real(rp), intent(in) :: gyu_(elem%np,lmesh%nea)
158 real(rp), intent(in) :: gzu_(elem%np,lmesh%nea)
159 real(rp), intent(in) :: gxv_(elem%np,lmesh%nea)
160 real(rp), intent(in) :: gyv_(elem%np,lmesh%nea)
161 real(rp), intent(in) :: gzv_(elem%np,lmesh%nea)
162 real(rp), intent(in) :: gxw_(elem%np,lmesh%nea)
163 real(rp), intent(in) :: gyw_(elem%np,lmesh%nea)
164 real(rp), intent(in) :: gzw_(elem%np,lmesh%nea)
165 real(rp), intent(in) :: gxpt_(elem%np,lmesh%nea)
166 real(rp), intent(in) :: gypt_(elem%np,lmesh%nea)
167 real(rp), intent(in) :: gzpt_(elem%np,lmesh%nea)
168 real(rp), intent(in) :: visccoef_h
169 real(rp), intent(in) :: visccoef_v
170 real(rp), intent(in) :: diffcoef_h
171 real(rp), intent(in) :: diffcoef_v
172
173 real(rp) :: fx(elem%np), fy(elem%np), fz(elem%np), liftdelflx(elem%np)
174 real(rp) :: del_flux(elem%nfptot,lmesh%ne,prog_vars_num)
175 real(rp) :: dens_(elem%np), rhot_hyd(elem%np), rhot_(elem%np), dpres_(elem%np)
176 real(rp) :: pres_(elem%np), u_(elem%np), v_(elem%np), w_(elem%np), pot_(elem%np)
177 real(rp) :: cori(elem%np)
178
179 real(rp) :: tmp(elem%np)
180 integer :: ke, ke2d
181 real(rp) :: gamm, rgamm
182 !------------------------------------------------------------------------
183
184 call prof_rapstart( 'cal_dyn_tend_bndflux', 3)
185 call cal_del_flux_dyn( del_flux, & ! (out)
186 ddens_, momx_, momy_, momz_, drhot_, dens_hyd, pres_hyd, & ! (in)
187 gxu_, gyu_, gzu_, gxv_, gyv_, gzv_, gxw_, gyw_, gzw_, & ! (in)
188 gxpt_, gypt_, gzpt_, & ! (in)
189 visccoef_h, visccoef_v, & ! (in)
190 diffcoef_h, diffcoef_v, & ! (in)
191 lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), & ! (in)
192 lmesh%vmapM, lmesh%vmapP, & ! (in)
193 lmesh, elem ) ! (in)
194 call prof_rapend( 'cal_dyn_tend_bndflux', 3)
195
196 !-----
197 call prof_rapstart( 'cal_dyn_tend_interior', 3)
198 gamm = cpdry / cvdry
199 rgamm = cvdry / cpdry
200
201 !$omp parallel do private(RHOT_hyd,RHOT_,pres_,dpres_,dens_,u_,v_,w_,pot_,ke2d,Cori,Fx,Fy,Fz,LiftDelFlx,tmp)
202 do ke = lmesh%NeS, lmesh%NeE
203 !--
204 rhot_hyd(:) = pres00/rdry * (pres_hyd(:,ke)/pres00)**rgamm
205 rhot_(:) = rhot_hyd(:) + drhot_(:,ke)
206 pres_(:) = pres_hyd(:,ke) * (1.0_rp + drhot_(:,ke)/rhot_hyd(:))**gamm
207 dpres_(:) = pres_(:) - pres_hyd(:,ke)
208 dens_(:) = ddens_(:,ke) + dens_hyd(:,ke)
209
210 u_(:) = momx_(:,ke)/dens_(:)
211 v_(:) = momy_(:,ke)/dens_(:)
212 w_(:) = momz_(:,ke)/dens_(:)
213 pot_(:) = rhot_(:)/dens_(:)
214
215 ke2d = lmesh%EMap3Dto2D(ke)
216 cori(:) = coriolis(elem%IndexH2Dto3D(:),ke2d)
217
218 !-- DENS
219 call sparsemat_matmul(dx, momx_(:,ke), fx)
220 call sparsemat_matmul(dy, momy_(:,ke), fy)
221 call sparsemat_matmul(lift, lmesh%Fscale(:,ke)*del_flux(:,ke,ddens_vid), liftdelflx)
222
223 dens_dt(:,ke) = - ( &
224 lmesh%Escale(:,ke,1,1) * fx(:) &
225 + lmesh%Escale(:,ke,2,2) * fy(:) &
226 + liftdelflx(:) )
227
228 !-- MOMX
229 call sparsemat_matmul(dx, u_(:)*momx_(:,ke) + pres_(:) - visccoef_h*dens_(:)*gxu_(:,ke), fx)
230 call sparsemat_matmul(dy, v_(:)*momx_(:,ke) - visccoef_h*dens_(:)*gyu_(:,ke), fy)
231 call sparsemat_matmul(dz, w_(:)*momx_(:,ke) - visccoef_v*dens_(:)*gzu_(:,ke) ,fz)
232 call sparsemat_matmul(lift, lmesh%Fscale(:,ke)*del_flux(:,ke,momx_vid), liftdelflx)
233
234 momx_dt(:,ke) = - ( &
235 lmesh%Escale(:,ke,1,1) * fx(:) &
236 + lmesh%Escale(:,ke,2,2) * fy(:) &
237 + lmesh%Escale(:,ke,3,3) * fz(:) &
238 + liftdelflx(:) &
239 - cori(:)*momy_(:,ke) &
240 )
241
242 !-- MOMY
243 call sparsemat_matmul(dx, u_(:)*momy_(:,ke) - visccoef_h*dens_(:)*gxv_(:,ke), fx)
244 call sparsemat_matmul(dy, v_(:)*momy_(:,ke) + pres_(:) - visccoef_h*dens_(:)*gyv_(:,ke), fy)
245 call sparsemat_matmul(dz, w_(:)*momy_(:,ke) - visccoef_v*dens_(:)*gzv_(:,ke) ,fz)
246 call sparsemat_matmul(lift, lmesh%Fscale(:,ke)*del_flux(:,ke,momy_vid), liftdelflx)
247
248 momy_dt(:,ke) = - ( &
249 lmesh%Escale(:,ke,1,1) * fx(:) &
250 + lmesh%Escale(:,ke,2,2) * fy(:) &
251 + lmesh%Escale(:,ke,3,3) * fz(:) &
252 + liftdelflx(:) &
253 + cori(:)*momx_(:,ke) &
254 )
255
256 !-- MOMZ
257 call sparsemat_matmul(dx, u_(:)*momz_(:,ke) - visccoef_h*dens_(:)*gxw_(:,ke), fx)
258 call sparsemat_matmul(dy, v_(:)*momz_(:,ke) - visccoef_h*dens_(:)*gyw_(:,ke), fy)
259 call sparsemat_matmul(dz, w_(:)*momz_(:,ke) - visccoef_v*dens_(:)*gzw_(:,ke), fz)
260 call sparsemat_matmul(lift, lmesh%Fscale(:,ke)*del_flux(:,ke,momz_vid), liftdelflx)
261
262 momz_dt(:,ke) = - ( &
263 lmesh%Escale(:,ke,1,1) * fx(:) &
264 + lmesh%Escale(:,ke,2,2) * fy(:) &
265 + lmesh%Escale(:,ke,3,3) * fz(:) &
266 + liftdelflx(:) )
267
268 !-- RHOT
269 call sparsemat_matmul(dx, pot_(:)*momx_(:,ke) - diffcoef_h*dens_(:)*gxpt_(:,ke), fx)
270 call sparsemat_matmul(dy, pot_(:)*momy_(:,ke) - diffcoef_h*dens_(:)*gypt_(:,ke), fy)
271 call sparsemat_matmul(dz, - diffcoef_v*dens_(:)*gzpt_(:,ke), fz)
272 call sparsemat_matmul(lift, lmesh%Fscale(:,ke)*del_flux(:,ke,drhot_vid), liftdelflx)
273
274 rhot_dt(:,ke) = - ( &
275 lmesh%Escale(:,ke,1,1) * fx(:) &
276 + lmesh%Escale(:,ke,2,2) * fy(:) &
277 + lmesh%Escale(:,ke,3,3) * fz(:) &
278 + liftdelflx(:) )
279 end do
280 call prof_rapend( 'cal_dyn_tend_interior', 3)
281
282 return
283 end subroutine atm_dyn_dgm_nonhydro3d_hevi_cal_tend
284
285 !------
286
287 subroutine cal_del_flux_dyn( del_flux, &
288 DDENS_, MOMX_, MOMY_, MOMZ_, DRHOT_, DENS_hyd, PRES_hyd, &
289 GxU_, GyU_, GzU_, GxV_, GyV_, GzV_, GxW_, GyW_, GzW_, &
290 GxPT_, GyPT_, GzPT_, &
291 viscCoef_h, viscCoef_v, &
292 diffCoef_h, diffCoef_v, &
293 nx, ny, nz, vmapM, vmapP, lmesh, elem )
294
295 implicit none
296
297 class(localmesh3d), intent(in) :: lmesh
298 class(elementbase3d), intent(in) :: elem
299 real(rp), intent(out) :: del_flux(elem%nfptot*lmesh%ne,prog_vars_num)
300 real(rp), intent(in) :: ddens_(elem%np*lmesh%nea)
301 real(rp), intent(in) :: momx_(elem%np*lmesh%nea)
302 real(rp), intent(in) :: momy_(elem%np*lmesh%nea)
303 real(rp), intent(in) :: momz_(elem%np*lmesh%nea)
304 real(rp), intent(in) :: drhot_(elem%np*lmesh%nea)
305 real(rp), intent(in) :: dens_hyd(elem%np*lmesh%nea)
306 real(rp), intent(in) :: pres_hyd(elem%np*lmesh%nea)
307 real(rp), intent(in) :: gxu_(elem%np*lmesh%nea)
308 real(rp), intent(in) :: gyu_(elem%np*lmesh%nea)
309 real(rp), intent(in) :: gzu_(elem%np*lmesh%nea)
310 real(rp), intent(in) :: gxv_(elem%np*lmesh%nea)
311 real(rp), intent(in) :: gyv_(elem%np*lmesh%nea)
312 real(rp), intent(in) :: gzv_(elem%np*lmesh%nea)
313 real(rp), intent(in) :: gxw_(elem%np*lmesh%nea)
314 real(rp), intent(in) :: gyw_(elem%np*lmesh%nea)
315 real(rp), intent(in) :: gzw_(elem%np*lmesh%nea)
316 real(rp), intent(in) :: gxpt_(elem%np*lmesh%nea)
317 real(rp), intent(in) :: gypt_(elem%np*lmesh%nea)
318 real(rp), intent(in) :: gzpt_(elem%np*lmesh%nea)
319 real(rp), intent(in) :: visccoef_h
320 real(rp), intent(in) :: visccoef_v
321 real(rp), intent(in) :: diffcoef_h
322 real(rp), intent(in) :: diffcoef_v
323 real(rp), intent(in) :: nx(elem%nfptot*lmesh%ne)
324 real(rp), intent(in) :: ny(elem%nfptot*lmesh%ne)
325 real(rp), intent(in) :: nz(elem%nfptot*lmesh%ne)
326 integer, intent(in) :: vmapm(elem%nfptot*lmesh%ne)
327 integer, intent(in) :: vmapp(elem%nfptot*lmesh%ne)
328
329 integer :: i, ip, im
330 real(rp) :: velp, velm, momz_p, alpha, swv
331 real(rp) :: presm, presp, dpresm, dpresp, densm, densp, rhotm, rhotp, rhot_hyd_m, rhot_hyd_p
332 real(rp) :: ddifffluxu, ddifffluxv, ddifffluxw, ddifffluxpt
333 real(rp) :: gamm, rgamm
334 real(rp) :: mu
335 logical :: visc_flag, diff_flag
336
337 !------------------------------------------------------------------------
338
339 gamm = cpdry/cvdry
340 rgamm = cvdry/cpdry
341
342 if (visccoef_h > 0.0_rp .or. visccoef_v > 0.0_rp) then
343 visc_flag = .true.
344 else
345 visc_flag = .false.
346 ddifffluxu = 0.0_rp
347 ddifffluxv = 0.0_rp
348 ddifffluxw = 0.0_rp
349 end if
350
351 if (diffcoef_h > 0.0_rp .or. diffcoef_v > 0.0_rp) then
352 diff_flag = .true.
353 else
354 diff_flag = .false.
355 ddifffluxpt = 0.0_rp
356 end if
357
358 !$omp parallel do &
359 !$omp private( iM, iP, VelP, VelM, MOMZ_P, alpha, swV, &
360 !$omp presM, presP, dpresM, dpresP, densM, densP, rhotM, rhotP, rhot_hyd_M, rhot_hyd_P ) &
361 !$omp firstprivate( dDiffFluxU, dDiffFluxV, dDiffFluxW, dDiffFluxPT, mu )
362 do i=1, elem%NfpTot*lmesh%Ne
363 im = vmapm(i); ip = vmapp(i)
364
365 rhot_hyd_m = pres00/rdry * (pres_hyd(im)/pres00)**rgamm
366 rhot_hyd_p = pres00/rdry * (pres_hyd(ip)/pres00)**rgamm
367
368 rhotm = rhot_hyd_m + drhot_(im)
369 presm = pres_hyd(im) * (1.0_rp + drhot_(im)/rhot_hyd_m)**gamm
370 dpresm = presm - pres_hyd(im)*abs(nz(i))
371
372 rhotp = rhot_hyd_p + drhot_(ip)
373 presp = pres_hyd(ip) * (1.0_rp + drhot_(ip)/rhot_hyd_p)**gamm
374 dpresp = presp - pres_hyd(ip)*abs(nz(i))
375
376 densm = ddens_(im) + dens_hyd(im)
377 densp = ddens_(ip) + dens_hyd(ip)
378
379 swv = 1.0_rp - nz(i)**2
380 velm = (momx_(im)*nx(i) + momy_(im)*ny(i) + momz_(im)*nz(i))/densm
381 velp = (momx_(ip)*nx(i) + momy_(ip)*ny(i) + momz_(ip)*nz(i))/densp
382 momz_p = momz_(ip)
383
384 alpha = swv*max( sqrt(gamm*presm/densm) + abs(velm), &
385 sqrt(gamm*presp/densp) + abs(velp) )
386
387 mu = (2.0_rp * dble((elem%PolyOrder_h+1)*(elem%PolyOrder_h+2)) / 2.0_rp / 600.0_rp)
388
389
390 if ( visc_flag ) then
391 ddifffluxu = ( &
392 visccoef_h*(densp*gxu_(ip) - densm*gxu_(im))*nx(i) &
393 + visccoef_h*(densp*gyu_(ip) - densm*gyu_(im))*ny(i) &
394 + visccoef_v*(densp*gzu_(ip) - densm*gzu_(im))*nz(i) &
395 + mu*(densp + densm)*(momx_(ip)/densp - momx_(im)/densm) )
396
397 ddifffluxv = ( &
398 visccoef_h*(densp*gxv_(ip) - densm*gxv_(im))*nx(i) &
399 + visccoef_h*(densp*gyv_(ip) - densm*gyv_(im))*ny(i) &
400 + visccoef_v*(densp*gzv_(ip) - densm*gzv_(im))*nz(i) &
401 + mu*(densp + densm)*(momy_(ip)/densp - momy_(im)/densm) )
402
403 ddifffluxw = ( &
404 visccoef_h*(densp*gxw_(ip) - densm*gxw_(im))*nx(i) &
405 + visccoef_h*(densp*gyw_(ip) - densm*gyw_(im))*ny(i) &
406 + visccoef_v*(densp*gzw_(ip) - densm*gzw_(im))*nz(i) &
407 + mu*(densp + densm)*(momz_(ip)/densp - momz_(im)/densm) )
408 end if
409 if ( diff_flag ) then
410 ddifffluxpt = ( &
411 diffcoef_h*(densp*gxpt_(ip) - densm*gxpt_(im))*nx(i) &
412 + diffcoef_h*(densp*gypt_(ip) - densm*gypt_(im))*ny(i) &
413 + diffcoef_v*(densp*gzpt_(ip) - densm*gzpt_(im))*nz(i) &
414 + mu*(densp + densm)*(rhotp/densp - rhotm/densm) )
415 end if
416
417 del_flux(i,ddens_vid) = 0.5_rp*( &
418 (momx_(ip) - momx_(im))*nx(i) &
419 + (momy_(ip) - momy_(im))*ny(i) &
420 - alpha * (ddens_(ip) - ddens_(im)) )
421
422 del_flux(i,momx_vid) = 0.5_rp*( &
423 ( momx_(ip)*velp - momx_(im)*velm ) &
424 + ( dpresp - dpresm )*nx(i) &
425 - alpha * (momx_(ip) - momx_(im)) &
426 - ddifffluxu )
427
428 del_flux(i,momy_vid) = 0.5_rp*( &
429 ( momy_(ip)*velp - momy_(im)*velm ) &
430 + ( dpresp - dpresm )*ny(i) &
431 - alpha * (momy_(ip) - momy_(im)) &
432 - ddifffluxv )
433
434 del_flux(i,momz_vid) = 0.5_rp*( &
435 ( momz_p*velp - momz_(im)*velm) &
436 - alpha * (momz_(ip) - momz_(im)) &
437 - ddifffluxw )
438
439 del_flux(i,drhot_vid) = 0.5_rp*( &
440 swv*( rhotp*velp - rhotm*velm ) &
441 - alpha *(drhot_(ip) - drhot_(im)) &
442 - ddifffluxpt )
443 end do
444
445 return
446 end subroutine cal_del_flux_dyn
447
448 subroutine cal_del_flux_dyn_ausmup( del_flux, &
449 DDENS_, MOMX_, MOMY_, MOMZ_, DRHOT_, DENS_hyd, PRES_hyd, &
450 GxU_, GyU_, GzU_, GxV_, GyV_, GzV_, GxW_, GyW_, GzW_, &
451 GxPT_, GyPT_, GzPT_, &
452 viscCoef_h, viscCoef_v, &
453 diffCoef_h, diffCoef_v, &
454 nx, ny, nz, vmapM, vmapP, lmesh, elem )
455
456 implicit none
457
458 class(localmesh3d), intent(in) :: lmesh
459 class(elementbase3d), intent(in) :: elem
460 real(rp), intent(out) :: del_flux(elem%nfptot*lmesh%ne,prog_vars_num)
461 real(rp), intent(in) :: ddens_(elem%np*lmesh%nea)
462 real(rp), intent(in) :: momx_(elem%np*lmesh%nea)
463 real(rp), intent(in) :: momy_(elem%np*lmesh%nea)
464 real(rp), intent(in) :: momz_(elem%np*lmesh%nea)
465 real(rp), intent(in) :: drhot_(elem%np*lmesh%nea)
466 real(rp), intent(in) :: dens_hyd(elem%np*lmesh%nea)
467 real(rp), intent(in) :: pres_hyd(elem%np*lmesh%nea)
468 real(rp), intent(in) :: gxu_(elem%np*lmesh%nea)
469 real(rp), intent(in) :: gyu_(elem%np*lmesh%nea)
470 real(rp), intent(in) :: gzu_(elem%np*lmesh%nea)
471 real(rp), intent(in) :: gxv_(elem%np*lmesh%nea)
472 real(rp), intent(in) :: gyv_(elem%np*lmesh%nea)
473 real(rp), intent(in) :: gzv_(elem%np*lmesh%nea)
474 real(rp), intent(in) :: gxw_(elem%np*lmesh%nea)
475 real(rp), intent(in) :: gyw_(elem%np*lmesh%nea)
476 real(rp), intent(in) :: gzw_(elem%np*lmesh%nea)
477 real(rp), intent(in) :: gxpt_(elem%np*lmesh%nea)
478 real(rp), intent(in) :: gypt_(elem%np*lmesh%nea)
479 real(rp), intent(in) :: gzpt_(elem%np*lmesh%nea)
480 real(rp), intent(in) :: visccoef_h
481 real(rp), intent(in) :: visccoef_v
482 real(rp), intent(in) :: diffcoef_h
483 real(rp), intent(in) :: diffcoef_v
484 real(rp), intent(in) :: nx(elem%nfptot*lmesh%ne)
485 real(rp), intent(in) :: ny(elem%nfptot*lmesh%ne)
486 real(rp), intent(in) :: nz(elem%nfptot*lmesh%ne)
487 integer, intent(in) :: vmapm(elem%nfptot*lmesh%ne)
488 integer, intent(in) :: vmapp(elem%nfptot*lmesh%ne)
489
490 integer :: i, ip, im
491 real(rp) :: velp, velm, momz_p, alpha, swv
492 real(rp) :: presm, presp, dpresm, dpresp, densm, densp, rhotm, rhotp, rhot_hyd_m, rhot_hyd_p
493 real(rp) :: mam, map, ma, mflx, csp, csm, betap, betam, chi, dpres_num
494 real(rp) :: ddifffluxu, ddifffluxv, ddifffluxw, ddifffluxpt
495 real(rp) :: gamm, rgamm
496 real(rp) :: mu
497 logical :: visc_flag, diff_flag
498
499 !------------------------------------------------------------------------
500
501 gamm = cpdry/cvdry
502 rgamm = cvdry/cpdry
503
504 if (visccoef_h > 0.0_rp .or. visccoef_v > 0.0_rp) then
505 visc_flag = .true.
506 else
507 visc_flag = .false.
508 ddifffluxu = 0.0_rp
509 ddifffluxv = 0.0_rp
510 ddifffluxw = 0.0_rp
511 end if
512
513 if (diffcoef_h > 0.0_rp .or. diffcoef_v > 0.0_rp) then
514 diff_flag = .true.
515 else
516 diff_flag = .false.
517 ddifffluxpt = 0.0_rp
518 end if
519
520 !$omp parallel do &
521 !$omp private( iM, iP, VelP, VelM, MOMZ_P, alpha, swV, &
522 !$omp presM, presP, dpresM, dpresP, densM, densP, rhotM, rhotP, rhot_hyd_M, rhot_hyd_P, &
523 !$omp MaM, MaP, Ma, mflx, CsP, CsM, betaP, betaM, chi, dpres_num ) &
524 !$omp firstprivate( dDiffFluxU, dDiffFluxV, dDiffFluxW, dDiffFluxPT, mu )
525 do i=1, elem%NfpTot*lmesh%Ne
526 im = vmapm(i); ip = vmapp(i)
527
528 rhot_hyd_m = pres00/rdry * (pres_hyd(im)/pres00)**rgamm
529 rhot_hyd_p = pres00/rdry * (pres_hyd(ip)/pres00)**rgamm
530
531 rhotm = rhot_hyd_m + drhot_(im)
532 presm = pres_hyd(im) * (1.0_rp + drhot_(im)/rhot_hyd_m)**gamm
533 dpresm = presm - pres_hyd(im)*abs(nz(i))
534
535 rhotp = rhot_hyd_p + drhot_(ip)
536 presp = pres_hyd(ip) * (1.0_rp + drhot_(ip)/rhot_hyd_p)**gamm
537 dpresp = presp - pres_hyd(ip)*abs(nz(i))
538
539 densm = ddens_(im) + dens_hyd(im)
540 densp = ddens_(ip) + dens_hyd(ip)
541
542 swv = 1.0_rp - nz(i)**2
543 velm = (momx_(im)*nx(i) + momy_(im)*ny(i) + momz_(im)*nz(i))/densm
544 velp = (momx_(ip)*nx(i) + momy_(ip)*ny(i) + momz_(ip)*nz(i))/densp
545 momz_p = momz_(ip)
546
547 csm = sqrt(gamm*presm/densm)
548 csp = sqrt(gamm*presp/densp)
549 mam = 2.0_rp * velm / (csm + csp)
550 map = 2.0_rp * velm / (csm + csp)
551 betam = 0.25_rp * ( 2.0_rp + mam ) * ( mam - 1.0_rp )**2
552 betap = 0.25_rp * ( 2.0_rp - map ) * ( map - 1.0_rp )**2
553 chi = min( 1.0_rp, &
554 2.0_rp/(csm + csp) * sqrt( 0.5_rp*( &
555 (momx_(im)**2 + momy_(im)**2 + momz_(im)**2)/densm**2 &
556 + (momx_(ip)**2 + momy_(ip)**2 + momz_(ip)**2)/densp**2 ) ) )
557 alpha = swv * ( densm * abs(velm) + densp * abs(velp) ) / ( densm + densp )
558
559 mflx = 0.5_rp * ( &
560 densm * velm + densp * velp &
561 - alpha*( ddens_(ip) - ddens_(im) ) &
562 - swv * 2.0_rp * (dpresp - dpresm)/(csm + csp) )
563
564 dpres_num = 0.5_rp * ( &
565 (1.0_rp + (1.0_rp - chi)*(betam + betap - 1.0_rp)) * (dpresm + dpresp) &
566 - (betap - betam)*(dpresp - dpresm) )
567
568 mu = (2.0_rp * dble((elem%PolyOrder_h+1)*(elem%PolyOrder_h+2)) / 2.0_rp / 600.0_rp)
569
570
571 if ( visc_flag ) then
572 ddifffluxu = ( &
573 visccoef_h*(densp*gxu_(ip) - densm*gxu_(im))*nx(i) &
574 + visccoef_h*(densp*gyu_(ip) - densm*gyu_(im))*ny(i) &
575 + visccoef_v*(densp*gzu_(ip) - densm*gzu_(im))*nz(i) &
576 + mu*(densp + densm)*(momx_(ip)/densp - momx_(im)/densm) )
577
578 ddifffluxv = ( &
579 visccoef_h*(densp*gxv_(ip) - densm*gxv_(im))*nx(i) &
580 + visccoef_h*(densp*gyv_(ip) - densm*gyv_(im))*ny(i) &
581 + visccoef_v*(densp*gzv_(ip) - densm*gzv_(im))*nz(i) &
582 + mu*(densp + densm)*(momy_(ip)/densp - momy_(im)/densm) )
583
584 ddifffluxw = ( &
585 visccoef_h*(densp*gxw_(ip) - densm*gxw_(im))*nx(i) &
586 + visccoef_h*(densp*gyw_(ip) - densm*gyw_(im))*ny(i) &
587 + visccoef_v*(densp*gzw_(ip) - densm*gzw_(im))*nz(i) &
588 + mu*(densp + densm)*(momz_(ip)/densp - momz_(im)/densm) )
589 end if
590 if ( diff_flag ) then
591 ddifffluxpt = ( &
592 diffcoef_h*(densp*gxpt_(ip) - densm*gxpt_(im))*nx(i) &
593 + diffcoef_h*(densp*gypt_(ip) - densm*gypt_(im))*ny(i) &
594 + diffcoef_v*(densp*gzpt_(ip) - densm*gzpt_(im))*nz(i) &
595 + mu*(densp + densm)*(rhotp/densp - rhotm/densm) )
596 end if
597
598 del_flux(i,ddens_vid) = swv * ( &
599 mflx &
600 - densm * velm )
601
602 del_flux(i,momx_vid) = 0.5_rp*( &
603 (mflx + abs(mflx)) * momx_(im) / densm &
604 + (mflx - abs(mflx)) * momx_(ip) / densp &
605 - 2.0_rp * momx_(im) * velm &
606 + 2.0_rp * (dpres_num - dpresm ) * nx(i) &
607 - ddifffluxu )
608
609 del_flux(i,momy_vid) = 0.5_rp*( &
610 (mflx + abs(mflx)) * momy_(im) / densm &
611 + (mflx - abs(mflx)) * momy_(ip) / densp &
612 - 2.0_rp * momy_(im) * velm &
613 + 2.0_rp * (dpres_num - dpresm ) * ny(i) &
614 - ddifffluxv )
615
616 del_flux(i,momz_vid) = 0.5_rp*( &
617 (mflx + abs(mflx)) * momz_(im) / densm &
618 + (mflx - abs(mflx)) * momz_(ip) / densp &
619 - 2.0_rp * momz_(im) * velm &
620 - ddifffluxw )
621
622 del_flux(i,drhot_vid) = 0.5_rp*( &
623 (mflx + abs(mflx)) * rhotm / densm &
624 + (mflx - abs(mflx)) * rhotp / densp &
625 - 2.0_rp * rhotm * velm &
626 - ddifffluxpt )
627 end do
628
629 return
630 end subroutine cal_del_flux_dyn_ausmup
631
632 subroutine atm_dyn_dgm_nonhydro3d_hevi_cal_grad_diffvars( &
633 GxU_, GyU_, GzU_, GxV_, GyV_, GzV_, GxW_, GyW_, GzW_, GxPT_, GyPT_, GzPT_, &
634 DDENS_, MOMX_, MOMY_, MOMZ_, DRHOT_, DENS_hyd, PRES_hyd, &
635 Dx, Dy, Dz, Lift, lmesh, elem )
636
637 implicit none
638
639 class(localmesh3d), intent(in) :: lmesh
640 class(elementbase3d), intent(in) :: elem
641 type(sparsemat), intent(in) :: dx, dy, dz, lift
642 real(rp), intent(out) :: gxu_(elem%np,lmesh%nea)
643 real(rp), intent(out) :: gyu_(elem%np,lmesh%nea)
644 real(rp), intent(out) :: gzu_(elem%np,lmesh%nea)
645 real(rp), intent(out) :: gxv_(elem%np,lmesh%nea)
646 real(rp), intent(out) :: gyv_(elem%np,lmesh%nea)
647 real(rp), intent(out) :: gzv_(elem%np,lmesh%nea)
648 real(rp), intent(out) :: gxw_(elem%np,lmesh%nea)
649 real(rp), intent(out) :: gyw_(elem%np,lmesh%nea)
650 real(rp), intent(out) :: gzw_(elem%np,lmesh%nea)
651 real(rp), intent(out) :: gxpt_(elem%np,lmesh%nea)
652 real(rp), intent(out) :: gypt_(elem%np,lmesh%nea)
653 real(rp), intent(out) :: gzpt_(elem%np,lmesh%nea)
654 real(rp), intent(in) :: ddens_(elem%np,lmesh%nea)
655 real(rp), intent(in) :: momx_(elem%np,lmesh%nea)
656 real(rp), intent(in) :: momy_(elem%np,lmesh%nea)
657 real(rp), intent(in) :: momz_(elem%np,lmesh%nea)
658 real(rp), intent(in) :: drhot_(elem%np,lmesh%nea)
659 real(rp), intent(in) :: dens_hyd(elem%np,lmesh%nea)
660 real(rp), intent(in) :: pres_hyd(elem%np,lmesh%nea)
661
662 real(rp) :: dens_(elem%np), u_(elem%np), v_(elem%np), w_(elem%np)
663 real(rp) :: dtheta_(elem%np), rhot_(elem%np), rhot_hyd(elem%np)
664 real(rp) :: fx(elem%np), fy(elem%np), fz(elem%np), liftdelflx(elem%np)
665 real(rp) :: del_flux(elem%nfptot,lmesh%ne,aux_diffvars_num)
666
667 integer :: ke
668 !------------------------------------------------------------------------------
669
670 call cal_del_graddiffvar( del_flux, & ! (out)
671 ddens_, momx_, momy_, momz_, drhot_, dens_hyd, pres_hyd, & ! (in)
672 lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), & ! (in)
673 lmesh%vmapM, lmesh%vmapP, & ! (in)
674 lmesh, elem ) ! (in)
675
676 do ke=lmesh%NeS, lmesh%NeE
677 dens_(:) = ddens_(:,ke) + dens_hyd(:,ke)
678 rhot_hyd(:) = pres00/rdry * (pres_hyd(:,ke)/pres00)**(cvdry/cpdry)
679 rhot_(:) = rhot_hyd(:) + drhot_(:,ke)
680
681 u_(:) = momx_(:,ke)/dens_(:)
682 v_(:) = momy_(:,ke)/dens_(:)
683 w_(:) = momz_(:,ke)/dens_(:)
684 dtheta_(:) = rhot_(:)/dens_(:) - rhot_hyd(:)/dens_hyd(:,ke)
685
686 !- U
687
688 call sparsemat_matmul(dx, u_, fx)
689 call sparsemat_matmul(lift, lmesh%Fscale(:,ke)*del_flux(:,ke,vars_gxu_id), liftdelflx)
690 gxu_(:,ke) = lmesh%Escale(:,ke,1,1)*fx(:) + liftdelflx(:)
691
692 call sparsemat_matmul(dy, u_, fy)
693 call sparsemat_matmul(lift, lmesh%Fscale(:,ke)*del_flux(:,ke,vars_gyu_id), liftdelflx)
694 gyu_(:,ke) = lmesh%Escale(:,ke,2,2)*fx(:) + liftdelflx(:)
695
696 call sparsemat_matmul(dz, u_, fz)
697 call sparsemat_matmul(lift, lmesh%Fscale(:,ke)*del_flux(:,ke,vars_gzu_id), liftdelflx)
698 gzu_(:,ke) = lmesh%Escale(:,ke,3,3)*fz(:) + liftdelflx(:)
699
700 !- V
701
702 call sparsemat_matmul(dx, v_, fx)
703 call sparsemat_matmul(lift, lmesh%Fscale(:,ke)*del_flux(:,ke,vars_gxv_id), liftdelflx)
704 gxv_(:,ke) = lmesh%Escale(:,ke,1,1)*fx(:) + liftdelflx(:)
705
706 call sparsemat_matmul(dy, v_, fy)
707 call sparsemat_matmul(lift, lmesh%Fscale(:,ke)*del_flux(:,ke,vars_gyv_id), liftdelflx)
708 gyv_(:,ke) = lmesh%Escale(:,ke,2,2)*fx(:) + liftdelflx(:)
709
710 call sparsemat_matmul(dz, v_, fz)
711 call sparsemat_matmul(lift, lmesh%Fscale(:,ke)*del_flux(:,ke,vars_gzv_id), liftdelflx)
712 gzv_(:,ke) = lmesh%Escale(:,ke,3,3)*fz(:) + liftdelflx(:)
713
714 !- W
715
716 call sparsemat_matmul(dx, w_, fx)
717 call sparsemat_matmul(lift, lmesh%Fscale(:,ke)*del_flux(:,ke,vars_gxw_id), liftdelflx)
718 gxw_(:,ke) = lmesh%Escale(:,ke,1,1)*fx(:) + liftdelflx(:)
719
720 call sparsemat_matmul(dy, w_, fx)
721 call sparsemat_matmul(lift, lmesh%Fscale(:,ke)*del_flux(:,ke,vars_gyw_id), liftdelflx)
722 gyw_(:,ke) = lmesh%Escale(:,ke,2,2)*fy(:) + liftdelflx(:)
723
724 call sparsemat_matmul(dz, w_, fz)
725 call sparsemat_matmul(lift, lmesh%Fscale(:,ke)*del_flux(:,ke,vars_gzw_id), liftdelflx)
726 gzw_(:,ke) = lmesh%Escale(:,ke,3,3)*fz(:) + liftdelflx(:)
727
728 !- PT
729
730 call sparsemat_matmul(dx, dtheta_, fx)
731 call sparsemat_matmul(lift, lmesh%Fscale(:,ke)*del_flux(:,ke,vars_gxpt_id), liftdelflx)
732 gxpt_(:,ke) = lmesh%Escale(:,ke,1,1)*fx(:) + liftdelflx(:)
733
734 call sparsemat_matmul(dy, dtheta_, fy)
735 call sparsemat_matmul(lift, lmesh%Fscale(:,ke)*del_flux(:,ke,vars_gypt_id), liftdelflx)
736 gypt_(:,ke) = lmesh%Escale(:,ke,2,2)*fz(:) + liftdelflx(:)
737
738 call sparsemat_matmul(dz, dtheta_, fz)
739 call sparsemat_matmul(lift, lmesh%Fscale(:,ke)*del_flux(:,ke,vars_gzpt_id), liftdelflx)
740 gzpt_(:,ke) = lmesh%Escale(:,ke,3,3)*fz(:) + liftdelflx(:)
741
742 end do
743
744 return
745 end subroutine atm_dyn_dgm_nonhydro3d_hevi_cal_grad_diffvars
746
747 subroutine cal_del_graddiffvar( del_flux, &
748 DDENS_, MOMX_, MOMY_, MOMZ_, DRHOT_, DENS_hyd, PRES_hyd, &
749 nx, ny, nz, vmapM, vmapP, lmesh, elem )
750
751 implicit none
752
753 class(localmesh3d), intent(in) :: lmesh
754 class(elementbase3d), intent(in) :: elem
755 real(rp), intent(out) :: del_flux(elem%nfptot*lmesh%ne,aux_diffvars_num)
756 real(rp), intent(in) :: ddens_(elem%np*lmesh%nea)
757 real(rp), intent(in) :: momx_(elem%np*lmesh%nea)
758 real(rp), intent(in) :: momy_(elem%np*lmesh%nea)
759 real(rp), intent(in) :: momz_(elem%np*lmesh%nea)
760 real(rp), intent(in) :: drhot_(elem%np*lmesh%nea)
761 real(rp), intent(in) :: dens_hyd(elem%np*lmesh%nea)
762 real(rp), intent(in) :: pres_hyd(elem%np*lmesh%nea)
763 real(rp), intent(in) :: nx(elem%nfptot*lmesh%ne)
764 real(rp), intent(in) :: ny(elem%nfptot*lmesh%ne)
765 real(rp), intent(in) :: nz(elem%nfptot*lmesh%ne)
766 integer, intent(in) :: vmapm(elem%nfptot*lmesh%ne)
767 integer, intent(in) :: vmapp(elem%nfptot*lmesh%ne)
768
769 integer :: i, ip, im
770 real(rp) :: velp, velm, alpha
771 real(rp) :: densm, densp, rhot_hyd, rhotm, rhotp
772 real(rp) :: delu, delv, delw, delpt
773 real(rp) :: momx_p, momy_p, momz_p
774 real(rp) :: gamm, rgamm
775 !------------------------------------------------------------------------
776
777 gamm = cpdry/cvdry
778 rgamm = cvdry/cpdry
779
780 do i=1, elem%NfpTot*lmesh%Ne
781 im = vmapm(i); ip = vmapp(i)
782
783 rhot_hyd = pres00/rdry * (pres_hyd(im)/pres00)**rgamm
784 rhotm = rhot_hyd + drhot_(im)
785 rhotp = rhot_hyd + drhot_(ip)
786 densm = ddens_(im) + dens_hyd(im)
787 densp = ddens_(ip) + dens_hyd(im)
788
789 delu = 0.5_rp*(momx_(ip)/densp - momx_(im)/densm)
790 delv = 0.5_rp*(momy_(ip)/densp - momy_(im)/densm)
791 delw = 0.5_rp*(momz_(ip)/densp - momz_(im)/densm)
792 delpt = 0.5_rp*(rhotp/densp - rhotm/densm)
793
794 !- U
795 del_flux(i,vars_gxu_id) = delu * nx(i)
796 del_flux(i,vars_gyu_id) = delu * ny(i)
797 del_flux(i,vars_gzu_id) = delu * nz(i)
798
799 !- V
800 del_flux(i,vars_gxv_id) = delv * nx(i)
801 del_flux(i,vars_gyv_id) = delv * ny(i)
802 del_flux(i,vars_gzv_id) = delv * nz(i)
803
804 !- W
805 del_flux(i,vars_gxw_id) = delw * nx(i)
806 del_flux(i,vars_gyw_id) = delw * ny(i)
807 del_flux(i,vars_gzw_id) = delw * nz(i)
808
809 !- PT
810 del_flux(i,vars_gxpt_id) = delpt * nx(i)
811 del_flux(i,vars_gypt_id) = delpt * ny(i)
812 del_flux(i,vars_gzpt_id) = delpt * nz(i)
813 end do
814
815 return
816 end subroutine cal_del_graddiffvar
817
818 subroutine atm_dyn_dgm_nonhydro3d_hevi_cal_vi( &
819 DENS_dt, MOMX_dt, MOMY_dt, MOMZ_dt, RHOT_dt, & ! (out)
820 ddens_, momx_, momy_, momz_, drhot_, dens_hyd, pres_hyd, & ! (in)
821 dz, lift, impl_fac, lmesh, elem, lmesh2d, elem2d )
822
823 implicit none
824
825 class(localmesh3d), intent(in) :: lmesh
826 class(elementbase3d), intent(in) :: elem
827 class(localmesh2d), intent(in) :: lmesh2d
828 class(elementbase2d), intent(in) :: elem2d
829 real(rp), intent(out) :: dens_dt(elem%np,lmesh%nea)
830 real(rp), intent(out) :: momx_dt(elem%np,lmesh%nea)
831 real(rp), intent(out) :: momy_dt(elem%np,lmesh%nea)
832 real(rp), intent(out) :: momz_dt(elem%np,lmesh%nea)
833 real(rp), intent(out) :: rhot_dt(elem%np,lmesh%nea)
834 real(rp), intent(in) :: ddens_(elem%np,lmesh%nea)
835 real(rp), intent(in) :: momx_(elem%np,lmesh%nea)
836 real(rp), intent(in) :: momy_(elem%np,lmesh%nea)
837 real(rp), intent(in) :: momz_(elem%np,lmesh%nea)
838 real(rp), intent(in) :: drhot_(elem%np,lmesh%nea)
839 real(rp), intent(in) :: dens_hyd(elem%np,lmesh%nea)
840 real(rp), intent(in) :: pres_hyd(elem%np,lmesh%nea)
841 class(sparsemat), intent(in) :: dz, lift
842 real(rp), intent(in) :: impl_fac
843
844 real(rp) :: prog_vars(elem%np,prog_vars_num,lmesh%nez)
845 real(rp) :: prog_vars0(elem%np,prog_vars_num,lmesh%nez)
846 real(rp) :: prog_vars00(elem%np,prog_vars_num,lmesh%nez)
847 real(rp) :: prog_del(elem%np,prog_vars_num,lmesh%nez)
848 real(rp) :: b(elem%np,prog_vars_num,lmesh%nez)
849 real(rp) :: ax(elem%np,prog_vars_num,lmesh%nez)
850 real(rp) :: tend(elem%np,prog_vars_num,lmesh%nez)
851 real(rp) :: dens_hyd_z(elem%np,lmesh%nez)
852 real(rp) :: pres_hyd_z(elem%np,lmesh%nez)
853 real(rp) :: nz(elem%nfptot,lmesh%nez)
854 integer :: vmapm(elem%nfptot,lmesh%nez)
855 integer :: vmapp(elem%nfptot,lmesh%nez)
856 integer :: ke_x, ke_y, ke_z, ke, p, v
857 integer :: itr_lin, itr_nlin
858 integer :: m, n
859 integer :: f, vs, ve
860 logical :: is_converged
861
862
863 type(gmres) :: gmres_hevi
864 real(rp), allocatable :: wj(:), pinv_v(:)
865 real(rp) :: pmatdlu(elem%np*prog_vars_num,elem%np*prog_vars_num,lmesh%nez)
866 integer :: pmatdlu_ipiv(elem%np*prog_vars_num,lmesh%nez)
867 real(rp) :: pmatl(elem%np,elem%np,prog_vars_num,prog_vars_num,lmesh%nez)
868 real(rp) :: pmatu(elem%np,elem%np,prog_vars_num,prog_vars_num,lmesh%nez)
869 real(rp), parameter :: eps0 = 1.0e-16_rp
870 real(rp), parameter :: eps = 1.0e-16_rp
871
872 !------------------------------------------------------------------------
873
874 call prof_rapstart( 'hevi_cal_vi', 3)
875 call prof_rapstart( 'hevi_cal_vi_prep', 3)
876
877 n = elem%Np * prog_vars_num * lmesh%NeZ
878 m = n / prog_vars_num / elem%Nnode_h1D**2
879
880 call gmres_hevi%Init( n, m, eps, eps )
881 allocate( wj(n), pinv_v(n) )
882
883 do ke_z=1, lmesh%NeZ
884 do f=1, elem%Nfaces_h
885 vs = 1 + (f-1)*elem%Nfp_h
886 ve = vs + elem%Nfp_h - 1
887 vmapm(vs:ve,ke_z) = elem%Fmask_h(:,f) + (ke_z-1)*elem%Np
888 end do
889 do f=1, elem%Nfaces_v
890 vs = elem%Nfp_h*elem%Nfaces_h + 1 + (f-1)*elem%Nfp_v
891 ve = vs + elem%Nfp_v - 1
892 vmapm(vs:ve,ke_z) = elem%Fmask_v(:,f) + (ke_z-1)*elem%Np
893 end do
894 vmapp(:,ke_z) = vmapm(:,ke_z)
895 end do
896
897 do ke_z=1, lmesh%NeZ
898 vs = elem%Nfp_h*elem%Nfaces_h + 1
899 ve = vs + elem%Nfp_v - 1
900 if (ke_z > 1) &
901 vmapp(vs:ve,ke_z) = elem%Fmask_v(:,2) + (ke_z-2)*elem%Np
902 ! if (ke_z > 1) then
903 ! vmapP(vs:ve,ke_z) = elem%Fmask_v(:,2) + (ke_z-2)*elem%Np
904 ! else
905 ! vmapP(vs:ve,ke_z) = elem%Fmask_v(:,2) + (lmesh%NeZ-1)*elem%Np
906 ! end if
907
908 vs = elem%Nfp_h*elem%Nfaces_h + elem%Nfp_v + 1
909 ve = vs + elem%Nfp_v - 1
910 if (ke_z < lmesh%NeZ) &
911 vmapp(vs:ve,ke_z) = elem%Fmask_v(:,1) + ke_z*elem%Np
912 ! if (ke_z < lmesh%NeZ) then
913 ! vmapP(vs:ve,ke_z) = elem%Fmask_v(:,1) + ke_z*elem%Np
914 ! else
915 ! vmapP(vs:ve,ke_z) = elem%Fmask_v(:,1)
916 ! end if
917 end do
918 call prof_rapend( 'hevi_cal_vi_prep', 3)
919
920 do ke_y=1, lmesh%NeY
921 do ke_x=1, lmesh%NeX
922
923 call prof_rapstart( 'hevi_cal_vi_get_var', 3)
924
925 do ke_z=1, lmesh%NeZ
926 ke = ke_x + (ke_y-1)*lmesh%NeX + (ke_z-1)*lmesh%NeX*lmesh%NeY
927
928 prog_vars(:,ddens_vid,ke_z) = ddens_(:,ke)
929 prog_vars(:,momx_vid,ke_z) = momx_(:,ke)
930 prog_vars(:,momy_vid,ke_z) = momy_(:,ke)
931 prog_vars(:,momz_vid,ke_z) = momz_(:,ke)
932 prog_vars(:,drhot_vid,ke_z) = drhot_(:,ke)
933 dens_hyd_z(:,ke_z) = dens_hyd(:,ke)
934 pres_hyd_z(:,ke_z) = pres_hyd(:,ke)
935
936 prog_vars0(:,:,ke_z) = prog_vars(:,:,ke_z)
937 prog_vars00(:,:,ke_z) = prog_vars(:,:,ke_z)
938 nz(:,ke_z) = lmesh%normal_fn(:,ke,3)
939 end do
940 call prof_rapend( 'hevi_cal_vi_get_var', 3)
941
942 if ( impl_fac > 0.0_rp ) then
943 call prof_rapstart( 'hevi_cal_vi_itr', 3)
944
945 do itr_nlin = 1, 1
946
947 do ke_z=1, lmesh%NeZ
948 prog_del(:,:,ke_z) = 0.0_rp
949 end do
950
951 call vi_eval_ax( ax(:,:,:), & ! (out)
952 prog_vars, prog_vars0, dens_hyd_z, pres_hyd_z, & ! (in)
953 dz, lift, impl_fac, lmesh, elem, & ! (in)
954 nz, vmapm, vmapp, ke_x, ke_y, .false. )
955
956 do ke_z=1, lmesh%NeZ
957 b(:,:,ke_z) = - ax(:,:,ke_z) + prog_vars00(:,:,ke_z)
958 end do
959
960 call prof_rapstart( 'hevi_cal_vi_pmatinv', 3)
961
962 call vi_construct_pmatinv( pmatdlu, pmatdlu_ipiv, pmatl, pmatu, & ! (out)
963 prog_vars0, dens_hyd_z, pres_hyd_z, & ! (in)
964 dz, lift, impl_fac, lmesh, elem, & ! (in)
965 nz, vmapm, vmapp, ke_x, ke_y )
966
967 call prof_rapend( 'hevi_cal_vi_pmatinv', 3)
968
969 call prof_rapstart( 'hevi_cal_vi_itr_lin', 3)
970 do itr_lin=1, 2*int(n/m)
971
972 call vi_gmres_core( gmres_hevi, prog_del(:,:,:), wj, & ! (inout)
973 is_converged, & ! (out)
974 prog_vars(:,:,:), b(:,:,:), n, m, pmatdlu, pmatdlu_ipiv, pmatl, pmatu, pinv_v, & ! (in)
975 dens_hyd_z, pres_hyd_z, & ! (in)
976 dz, lift, impl_fac, lmesh, elem, & ! (in)
977 nz, vmapm, vmapp, ke_x, ke_y )
978
979 if (is_converged) exit
980 end do ! itr lin
981! write(*,*) lmesh%tileID, ke_x, ke_y, itr_lin, gmres_hevi%m_out
982 call prof_rapend( 'hevi_cal_vi_itr_lin', 3)
983
984 do ke_z=1, lmesh%NeZ
985 prog_vars(:,:,ke_z) = prog_vars(:,:,ke_z) + prog_del(:,:,ke_z)
986 prog_vars0(:,:,ke_z) = prog_vars(:,:,ke_z)
987 end do
988 end do ! itr nlin
989
990 call prof_rapend( 'hevi_cal_vi_itr', 3)
991 end if
992
993 call prof_rapstart( 'hevi_cal_vi_retrun_var', 3)
994 call vi_eval_ax( tend(:,:,:), & ! (out)
995 prog_vars, prog_vars, dens_hyd_z, pres_hyd_z, & ! (in)
996 dz, lift, impl_fac, lmesh, elem, & ! (in)
997 nz, vmapm, vmapp, ke_x, ke_y, .true. )
998
999 !tend(:,:,:) = 0.0_RP
1000 do ke_z=1, lmesh%NeZ
1001 ke = ke_x + (ke_y-1)*lmesh%NeX + (ke_z-1)*lmesh%NeX*lmesh%NeY
1002 dens_dt(:,ke) = - tend(:,ddens_vid,ke_z)
1003 momx_dt(:,ke) = - tend(:,momx_vid,ke_z)
1004 momy_dt(:,ke) = - tend(:,momy_vid,ke_z)
1005 momz_dt(:,ke) = - tend(:,momz_vid,ke_z)
1006 rhot_dt(:,ke) = - tend(:,drhot_vid,ke_z)
1007 end do
1008 !tend(:,:,:) = 0.0_RP
1009 call prof_rapend( 'hevi_cal_vi_retrun_var', 3)
1010 end do
1011 end do
1012
1013 call gmres_hevi%Final()
1014
1015 call prof_rapend( 'hevi_cal_vi', 3)
1016
1017 return
1018 end subroutine atm_dyn_dgm_nonhydro3d_hevi_cal_vi
1019
1020 !------------------------------------------------
1021
1022 subroutine vi_gmres_core( gmres_hevi, x, wj, is_converged, & ! (inout)
1023 x0, b, n, m, pmatdlu, pmatdlu_ipiv, pmatl, pmatu, pinv_v, & ! (in)
1024 dens_hyd, pres_hyd, & ! (in)
1025 dz, lift, impl_fac, lmesh, elem, & ! (in)
1026 nz, vmapm, vmapp, ke_x, ke_y )
1027
1028 implicit none
1029
1030 class(localmesh3d), intent(in) :: lmesh
1031 class(elementbase3d), intent(in) :: elem
1032 integer, intent(in) :: n
1033 integer, intent(in) :: m
1034
1035 class(gmres), intent(inout) :: gmres_hevi
1036 real(rp), intent(inout) :: x(n)
1037 real(rp), intent(inout) :: wj(n)
1038 logical, intent(out) :: is_converged
1039 real(rp), intent(in) :: x0(n)
1040 real(rp), intent(in) :: b(n)
1041 real(rp), intent(in) :: pmatdlu(elem%np*prog_vars_num,elem%np*prog_vars_num,lmesh%nez)
1042 integer, intent(in) :: pmatdlu_ipiv(elem%np*prog_vars_num,lmesh%nez)
1043 real(rp), intent(in) :: pmatl(elem%np,elem%np,prog_vars_num,prog_vars_num,lmesh%nez)
1044 real(rp), intent(in) :: pmatu(elem%np,elem%np,prog_vars_num,prog_vars_num,lmesh%nez)
1045 real(rp), intent(inout) :: pinv_v(n)
1046 !---
1047 real(rp), intent(in) :: dens_hyd(elem%np,lmesh%nez)
1048 real(rp), intent(in) :: pres_hyd(elem%np,lmesh%nez)
1049 class(sparsemat), intent(in) :: dz, lift
1050 real(rp), intent(in) :: impl_fac
1051 real(rp), intent(in) :: nz(elem%nfptot,lmesh%nez)
1052 integer, intent(in) :: vmapm(elem%nfptot,lmesh%nez)
1053 integer, intent(in) :: vmapp(elem%nfptot,lmesh%nez)
1054 integer, intent(in) :: ke_x, ke_y
1055
1056 integer :: j
1057
1058 !--------------------------------------
1059
1060! call PROF_rapstart('vi_lin_core_eval_Ax', 3)
1061 call vi_eval_ax_lin( wj(:), & ! (out)
1062 x, x0, dens_hyd, pres_hyd, & ! (in)
1063 dz, lift, impl_fac, lmesh, elem, & ! (in)
1064 nz, vmapm, vmapp, ke_x, ke_y, .false. ) ! (in)
1065! call PROF_rapend('vi_lin_core_eval_Ax', 3)
1066
1067 ! call PROF_rapstart('vi_lin_core_pre', 3)
1068 call gmres_hevi%Iterate_pre( b, wj, is_converged)
1069 if (is_converged) return
1070 ! call PROF_rapend('vi_lin_core_pre', 3)
1071
1072 do j=1, min(m, n)
1073 call prof_rapstart('vi_lin_core_pinv_v', 3)
1074 call matmul_pinv_v( pinv_v, pmatdlu, pmatdlu_ipiv, pmatl, pmatu, gmres_hevi%v(:,j) )
1075 call prof_rapend('vi_lin_core_pinv_v', 3)
1076
1077 call prof_rapstart('vi_lin_core_eval_Ax', 3)
1078 call vi_eval_ax_lin( wj(:), & ! (out)
1079 pinv_v, x0, dens_hyd, pres_hyd, & ! (in)
1080 dz, lift, impl_fac, lmesh, elem, & ! (in)
1081 nz, vmapm, vmapp, ke_x, ke_y, .false. ) ! (in)
1082 call prof_rapend('vi_lin_core_eval_Ax', 3)
1083
1084 call prof_rapstart('vi_lin_core_itr', 3)
1085 call gmres_hevi%Iterate_step_j( j, wj, is_converged)
1086 call prof_rapend('vi_lin_core_itr', 3)
1087 if (is_converged) exit
1088 end do
1089
1090 call prof_rapstart('vi_lin_core_post', 3)
1091 do j=1, n
1092 wj(j) = 0.0_rp
1093 end do
1094 call gmres_hevi%Iterate_post( wj )
1095 call prof_rapstart('vi_lin_core_pinv_v_px0', 3)
1096 call matmul_pinv_v_plus_x0( x, pmatdlu, pmatdlu_ipiv, pmatl, pmatu, wj)
1097 call prof_rapend('vi_lin_core_pinv_v_px0', 3)
1098 call prof_rapend('vi_lin_core_post', 3)
1099
1100 return
1101 contains
1102 subroutine matmul_pinv_v( pinv_v_, pDlu_, PmatDlu_ipiv_, pL, pU, v)
1103 implicit none
1104 real(rp), intent(out) :: pinv_v_(elem%np,prog_vars_num,lmesh%nez)
1105 real(rp), intent(in) :: pdlu_(elem%np*prog_vars_num,elem%np*prog_vars_num,lmesh%nez)
1106 integer, intent(in) :: pmatdlu_ipiv_(elem%np*prog_vars_num,lmesh%nez)
1107 real(rp), intent(in) :: pl(elem%np,elem%np,prog_vars_num,prog_vars_num,lmesh%nez)
1108 real(rp), intent(in) :: pu(elem%np,elem%np,prog_vars_num,prog_vars_num,lmesh%nez)
1109 real(rp), intent(in) :: v(elem%np*prog_vars_num,lmesh%nez)
1110
1111 integer :: k, n
1112 real(rp) :: tmp(elem%np*prog_vars_num)
1113 integer :: vid, vs, ve
1114 integer :: info
1115 !------------------------------------
1116
1117 n = elem%Np * prog_vars_num
1118
1119 do vid = 1, prog_vars_num
1120 vs = 1 + (vid-1) * elem%Np
1121 ve = vs + elem%Np - 1
1122 pinv_v_(:,vid,1) = v(vs:ve,1)
1123 end do
1124 call dgetrs('N', n, 1, pdlu_(:,:,1), n, pmatdlu_ipiv_(:,1), pinv_v_(:,:,1), n, info)
1125
1126! tmp(:) = matmul( pDlu_(:,:,1), v(:,1) )
1127! do vid = 1, PROG_VARS_NUM
1128! ! pinv_v_(:,vid,1) = matmul( pDlu_(:,:,vid,1), v(:,1) )
1129! vs = 1 + (vid-1) * elem%Np
1130! ve = vs + elem%Np - 1
1131! pinv_v_(:,vid,1) = tmp(vs:ve)
1132! end do
1133
1134 do k=2, lmesh%NeZ
1135 vs = 1; ve = elem%Np
1136 pinv_v_(:,ddens_vid,k) = v(vs:ve,k) &
1137 - matmul( pl(:,:,ddens_vid,ddens_vid,k), pinv_v_(:,ddens_vid,k-1) ) &
1138 - matmul( pl(:,:,momz_vid ,ddens_vid,k), pinv_v_(:,momz_vid ,k-1) )
1139
1140 vs = ve+1; ve = vs + elem%Np - 1
1141 pinv_v_(:,momx_vid,k) = v(vs:ve,k) &
1142 - matmul( pl(:,:,momx_vid,momx_vid,k), pinv_v_(:,momx_vid,k-1) )
1143
1144 vs = ve+1; ve = vs + elem%Np - 1
1145 pinv_v_(:,momy_vid,k) = v(vs:ve,k) &
1146 - matmul( pl(:,:,momy_vid,momy_vid,k), pinv_v_(:,momy_vid,k-1) )
1147
1148 vs = ve+1; ve = vs + elem%Np - 1
1149 pinv_v_(:,momz_vid,k) = v(vs:ve,k) &
1150 - matmul( pl(:,:,ddens_vid,momz_vid,k), pinv_v_(:,ddens_vid,k-1) ) &
1151 - matmul( pl(:,:,momz_vid ,momz_vid,k), pinv_v_(:,momz_vid ,k-1) ) &
1152 - matmul( pl(:,:,drhot_vid,momz_vid,k), pinv_v_(:,drhot_vid,k-1) )
1153
1154 vs = ve+1; ve = vs + elem%Np - 1
1155 pinv_v_(:,drhot_vid,k) = v(vs:ve,k) &
1156 - matmul( pl(:,:,ddens_vid ,drhot_vid,k), pinv_v_(:,ddens_vid,k-1) ) &
1157 - matmul( pl(:,:,momz_vid ,drhot_vid,k), pinv_v_(:,momz_vid ,k-1) ) &
1158 - matmul( pl(:,:,drhot_vid ,drhot_vid,k), pinv_v_(:,drhot_vid,k-1) )
1159
1160 call dgetrs('N', n, 1, pdlu_(:,:,k), n, pmatdlu_ipiv_(:,k), pinv_v_(:,:,k), n, info)
1161
1162 ! !$omp parallel do
1163 ! do vid = 1, PROG_VARS_NUM
1164 ! pinv_v_(:,vid,k) = matmul( pDlu_(:,:,vid,k), tmp(:) )
1165 ! end do
1166
1167 ! do vid = 1, PROG_VARS_NUM
1168 ! vs = 1 + (vid-1) * elem%Np; ve = vs + elem%Np - 1
1169 ! tmp(vs:ve) = pinv_v_(:,vid,k)
1170 ! end do
1171 ! tmp(:) = matmul( pDlu_(:,:,k), tmp(:) )
1172 ! do vid = 1, PROG_VARS_NUM
1173 ! vs = 1 + (vid-1) * elem%Np; ve = vs + elem%Np - 1
1174 ! pinv_v_(:,vid,k) = tmp(vs:ve)
1175 ! end do
1176 end do
1177
1178 !
1179 do k=lmesh%NeZ-1, 1, -1
1180 vs = 1; ve = elem%Np
1181 tmp(vs:ve) = &
1182 matmul( pu(:,:,ddens_vid,ddens_vid,k), pinv_v_(:,ddens_vid,k+1) ) &
1183 + matmul( pu(:,:,momz_vid ,ddens_vid,k), pinv_v_(:,momz_vid ,k+1) )
1184
1185 vs = ve+1; ve = vs + elem%Np - 1
1186 tmp(vs:ve) = &
1187 matmul( pu(:,:,momx_vid,momx_vid,k), pinv_v_(:,momx_vid,k+1) )
1188
1189 vs = ve+1; ve = vs + elem%Np - 1
1190 tmp(vs:ve) = &
1191 matmul( pu(:,:,momy_vid,momy_vid,k), pinv_v_(:,momy_vid,k+1) )
1192
1193 vs = ve+1; ve = vs + elem%Np - 1
1194 tmp(vs:ve) = &
1195 matmul( pu(:,:,ddens_vid,momz_vid,k), pinv_v_(:,ddens_vid,k+1) ) &
1196 + matmul( pu(:,:,momz_vid ,momz_vid,k), pinv_v_(:,momz_vid ,k+1) ) &
1197 + matmul( pu(:,:,drhot_vid,momz_vid,k), pinv_v_(:,drhot_vid,k+1) )
1198
1199 vs = ve+1; ve = vs + elem%Np - 1
1200 tmp(vs:ve) = &
1201 matmul( pu(:,:,ddens_vid ,drhot_vid,k), pinv_v_(:,ddens_vid,k+1) ) &
1202 + matmul( pu(:,:,momz_vid ,drhot_vid,k), pinv_v_(:,momz_vid ,k+1) ) &
1203 + matmul( pu(:,:,drhot_vid ,drhot_vid,k), pinv_v_(:,drhot_vid,k+1) )
1204
1205 call dgetrs('N', n, 1, pdlu_(:,:,k), n, pmatdlu_ipiv_(:,k), tmp(:), n, info)
1206
1207 !! $omp parallel do
1208 do vid = 1, prog_vars_num
1209 vs = 1 + (vid-1) * elem%Np
1210 ve = vs + elem%Np - 1
1211 pinv_v_(:,vid,k) = pinv_v_(:,vid,k) - tmp(vs:ve)
1212 end do
1213
1214 ! tmp(:) = matmul( pDlu_(:,:,k), tmp(:) )
1215 ! do vid = 1, PROG_VARS_NUM
1216 ! vs = 1 + (vid-1) * elem%Np; ve = vs + elem%Np - 1
1217 ! pinv_v_(:,vid,k) = pinv_v_(:,vid,k) - tmp(vs:ve)
1218 ! end do
1219
1220 end do
1221
1222 return
1223 end subroutine matmul_pinv_v
1224
1225 subroutine matmul_pinv_v_plus_x0( x_, pDlu_, PmatDlu_ipiv_, pL, pU, v)
1226 implicit none
1227 real(rp), intent(inout) :: x_(elem%np*prog_vars_num,lmesh%nez)
1228 real(rp), intent(in) :: pdlu_(elem%np*prog_vars_num,elem%np*prog_vars_num,lmesh%nez)
1229 integer, intent(in) :: pmatdlu_ipiv_(elem%np*prog_vars_num,lmesh%nez)
1230 real(rp), intent(in) :: pl(elem%np*prog_vars_num,elem%np*prog_vars_num,lmesh%nez)
1231 real(rp), intent(in) :: pu(elem%np*prog_vars_num,elem%np*prog_vars_num,lmesh%nez)
1232 real(rp), intent(in) :: v(elem%np*prog_vars_num,lmesh%nez)
1233
1234 integer :: k
1235 real(rp) :: tmp(elem%np*prog_vars_num,lmesh%nez)
1236
1237 !------------------------------------
1238
1239 call matmul_pinv_v( tmp, pdlu_, pmatdlu_ipiv_, pl, pu, v)
1240 !$omp parallel do
1241 do k=1, lmesh%NeZ
1242 x_(:,k) = x_(:,k) + tmp(:,k)
1243 end do
1244
1245 return
1246 end subroutine matmul_pinv_v_plus_x0
1247 end subroutine vi_gmres_core
1248
1249 !---
1250 subroutine vi_eval_ax( Ax, &
1251 PROG_VARS, PROG_VARS0, DENS_hyd, PRES_hyd, & ! (in)
1252 dz, lift, impl_fac, lmesh, elem, & ! (in)
1253 nz, vmapm, vmapp, ke_x, ke_y, cal_tend_flag )
1254
1255 implicit none
1256
1257 class(localmesh3d), intent(in) :: lmesh
1258 class(elementbase3d), intent(in) :: elem
1259 real(rp), intent(out) :: ax(elem%np,prog_vars_num,lmesh%nez)
1260 real(rp), intent(in) :: prog_vars(elem%np,prog_vars_num,lmesh%nez)
1261 real(rp), intent(in) :: prog_vars0(elem%np,prog_vars_num,lmesh%nez)
1262 real(rp), intent(in) :: dens_hyd(elem%np,lmesh%nez)
1263 real(rp), intent(in) :: pres_hyd(elem%np,lmesh%nez)
1264 class(sparsemat), intent(in) :: dz, lift
1265 real(rp), intent(in) :: impl_fac
1266 real(rp), intent(in) :: nz(elem%nfptot,lmesh%nez)
1267 integer, intent(in) :: vmapm(elem%nfptot,lmesh%nez)
1268 integer, intent(in) :: vmapp(elem%nfptot,lmesh%nez)
1269 integer, intent(in) :: ke_x, ke_y
1270 logical, intent(in) :: cal_tend_flag
1271
1272 real(rp) :: fz(elem%np), liftdelflx(elem%np)
1273 real(rp) :: del_flux(elem%nfptot,lmesh%nez,prog_vars_num)
1274 real(rp) :: rhot_hyd(elem%np), pot(elem%np)
1275 real(rp) :: dpres(elem%np)
1276 integer :: ke_z
1277 integer :: ke
1278 real(rp) :: gamm, rgamm
1279
1280 !--------------------------------------------------------
1281
1282 gamm = cpdry/cvdry
1283 rgamm = cvdry/cpdry
1284
1285 call vi_cal_del_flux_dyn( del_flux, & ! (out)
1286 prog_vars(:,ddens_vid,:), prog_vars(:,momx_vid,:), & ! (in)
1287 prog_vars(:,momy_vid ,:), prog_vars(:,momz_vid,:), & ! (in)
1288 prog_vars(:,drhot_vid,:), & ! (in)
1289 prog_vars0(:,ddens_vid,:), prog_vars0(:,momx_vid,:), & ! (in)
1290 prog_vars0(:,momy_vid ,:), prog_vars0(:,momz_vid,:), & ! (in)
1291 prog_vars0(:,drhot_vid,:), & ! (in)
1292 dens_hyd, pres_hyd, nz, vmapm, vmapp, & ! (in)
1293 lmesh, elem ) ! (in)
1294
1295 !$omp parallel do private(ke, RHOT_hyd, DPRES, POT, Fz, LiftDelFlx)
1296 do ke_z=1, lmesh%NeZ
1297 ke = ke_x + (ke_y-1)*lmesh%NeX + (ke_z-1)*lmesh%NeX*lmesh%NeY
1298
1299 rhot_hyd(:) = pres00/rdry * (pres_hyd(:,ke_z)/pres00)**rgamm
1300
1301 dpres(:) = pres_hyd(:,ke_z) * ((1.0_rp + prog_vars(:,drhot_vid,ke_z)/rhot_hyd(:))**gamm - 1.0_rp)
1302 pot(:) = (rhot_hyd(:) + prog_vars(:,drhot_vid,ke_z))/(dens_hyd(:,ke_z) + prog_vars(:,ddens_vid,ke_z))
1303
1304 !- DENS
1305 call sparsemat_matmul(dz, prog_vars(:,momz_vid,ke_z), fz)
1306 call sparsemat_matmul(lift, lmesh%Fscale(:,ke)*del_flux(:,ke_z,ddens_vid), liftdelflx)
1307 ax(:,ddens_vid,ke_z) = lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:)
1308
1309 !- MOMX
1310 call sparsemat_matmul(lift, lmesh%Fscale(:,ke)*del_flux(:,ke_z,momx_vid), liftdelflx)
1311 ax(:,momx_vid,ke_z) = liftdelflx(:)
1312
1313 !-MOMY
1314 call sparsemat_matmul(lift, lmesh%Fscale(:,ke)*del_flux(:,ke_z,momy_vid), liftdelflx)
1315 ax(:,momy_vid,ke_z) = liftdelflx(:)
1316
1317 !-MOMZ
1318 call sparsemat_matmul(dz, dpres(:), fz)
1319 call sparsemat_matmul(lift, lmesh%Fscale(:,ke)*del_flux(:,ke_z,momz_vid), liftdelflx)
1320 ax(:,momz_vid,ke_z) = lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) &
1321 + grav * matmul(intrpmat_vpordm1, prog_vars(:,ddens_vid,ke_z))
1322
1323 !-RHOT
1324 call sparsemat_matmul(dz, pot(:)*prog_vars(:,momz_vid,ke_z), fz)
1325 call sparsemat_matmul(lift, lmesh%Fscale(:,ke)*del_flux(:,ke_z,drhot_vid), liftdelflx)
1326 ax(:,drhot_vid,ke_z) = lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:)
1327
1328 !--
1329 if ( .not. cal_tend_flag ) then
1330 ax(:,:,ke_z) = prog_vars(:,:,ke_z) + impl_fac * ax(:,:,ke_z)
1331 end if
1332
1333 end do
1334
1335 return
1336 end subroutine vi_eval_ax
1337
1338 subroutine vi_cal_del_flux_dyn( del_flux, &
1339 DDENS_, MOMX_, MOMY_, MOMZ_, DRHOT_, &
1340 DDENS0_, MOMX0_, MOMY0_, MOMZ0_, DRHOT0_, &
1341 DENS_hyd, PRES_hyd, nz, vmapM, vmapP, lmesh, elem )
1342
1343 implicit none
1344
1345 class(localmesh3d), intent(in) :: lmesh
1346 class(elementbase3d), intent(in) :: elem
1347 real(rp), intent(out) :: del_flux(elem%nfptot*lmesh%nez,prog_vars_num)
1348 real(rp), intent(in) :: ddens_(elem%np*lmesh%nez)
1349 real(rp), intent(in) :: momx_(elem%np*lmesh%nez)
1350 real(rp), intent(in) :: momy_(elem%np*lmesh%nez)
1351 real(rp), intent(in) :: momz_(elem%np*lmesh%nez)
1352 real(rp), intent(in) :: drhot_(elem%np*lmesh%nez)
1353 real(rp), intent(in) :: ddens0_(elem%np*lmesh%nez)
1354 real(rp), intent(in) :: momx0_(elem%np*lmesh%nez)
1355 real(rp), intent(in) :: momy0_(elem%np*lmesh%nez)
1356 real(rp), intent(in) :: momz0_(elem%np*lmesh%nez)
1357 real(rp), intent(in) :: drhot0_(elem%np*lmesh%nez)
1358 real(rp), intent(in) :: dens_hyd(elem%np*lmesh%nez)
1359 real(rp), intent(in) :: pres_hyd(elem%np*lmesh%nez)
1360 real(rp), intent(in) :: nz(elem%nfptot*lmesh%nez)
1361 integer, intent(in) :: vmapm(elem%nfptot*lmesh%nez)
1362 integer, intent(in) :: vmapp(elem%nfptot*lmesh%nez)
1363
1364 integer :: i, p, ke_z, ip, im
1365 real(rp) :: alpha0, swv
1366 real(rp) :: momz_p
1367 real(rp) :: rhot_hyd_m, rhot_hyd_p
1368 real(rp) :: dpresm, dpresp, densm, densp, pottm, pottp
1369 real(rp) :: pres0m, pres0p, dens0m, dens0p
1370 real(rp) :: gamm, rgamm
1371 !------------------------------------------------------------------------
1372
1373 gamm = cpdry/cvdry
1374 rgamm = cvdry/cpdry
1375
1376 !$omp parallel do private( p, i, iM, iP, &
1377 !$omp rhot_hyd_M, rhot_hyd_P, densM, densP, pottM, pottP, &
1378 !$omp dpresM, dpresP, MOMZ_P, &
1379 !$omp dens0M, dens0P, pres0M, pres0P, &
1380 !$omp swV, alpha0 )
1381 do ke_z=1, lmesh%NeZ
1382 do p=1, elem%NfpTot
1383 i = p + (ke_z-1)*elem%NfpTot
1384 im = vmapm(i); ip = vmapp(i)
1385
1386 !-
1387 rhot_hyd_m = pres00/rdry * (pres_hyd(im)/pres00)**rgamm
1388 rhot_hyd_p = pres00/rdry * (pres_hyd(ip)/pres00)**rgamm
1389
1390 densm = dens_hyd(im) + ddens_(im)
1391 densp = dens_hyd(ip) + ddens_(ip)
1392
1393 pottm = (rhot_hyd_m + drhot_(im)) / densm
1394 pottp = (rhot_hyd_p + drhot_(ip)) / densp
1395
1396 dpresm = pres_hyd(im) * ((1.0_rp + drhot_(im)/rhot_hyd_m)**gamm - 1.0_rp)
1397 dpresp = pres_hyd(ip) * ((1.0_rp + drhot_(ip)/rhot_hyd_p)**gamm - 1.0_rp)
1398
1399 !-
1400 dens0m = dens_hyd(im) + ddens0_(im)
1401 dens0p = dens_hyd(ip) + ddens0_(ip)
1402
1403 pres0m = pres_hyd(im) * (1.0_rp + drhot0_(im)/rhot_hyd_m)**gamm
1404 pres0p = pres_hyd(ip) * (1.0_rp + drhot0_(ip)/rhot_hyd_p)**gamm
1405
1406 swv = nz(i)**2
1407 alpha0 = swv * max( abs(momz0_(im)/dens0m) + sqrt(gamm * pres0m/dens0m), &
1408 abs(momz0_(ip)/dens0p) + sqrt(gamm * pres0p/dens0p) )
1409
1410 if (im==ip .and. (ke_z == 1 .or. ke_z == lmesh%NeZ)) then
1411 momz_p = - momz_(im)
1412 else
1413 momz_p = momz_(ip)
1414 end if
1415
1416 del_flux(i,ddens_vid) = 0.5_rp * ( &
1417 + ( momz_p - momz_(im) ) * nz(i) &
1418 - alpha0 * ( ddens_(ip) - ddens_(im) ) )
1419
1420 del_flux(i,momx_vid) = 0.5_rp * ( &
1421 - alpha0 * ( momx_(ip) - momx_(im) ) )
1422
1423 del_flux(i,momy_vid) = 0.5_rp * ( &
1424 - alpha0 * ( momy_(ip) - momy_(im) ) )
1425
1426 del_flux(i,momz_vid) = 0.5_rp * ( &
1427 + ( dpresp - dpresm ) * nz(i) &
1428 - alpha0 * ( momz_p - momz_(im) ) )
1429
1430 del_flux(i,drhot_vid) = 0.5_rp * ( &
1431 + ( pottp * momz_p - pottm * momz_(im) ) * nz(i) &
1432 - alpha0 * ( drhot_(ip) - drhot_(im) ) )
1433 end do
1434 end do
1435
1436 return
1437 end subroutine vi_cal_del_flux_dyn
1438!-
1439 subroutine vi_eval_ax_lin( Ax, &
1440 PROG_VARS, PROG_VARS0, DENS_hyd, PRES_hyd, & ! (in)
1441 dz, lift, impl_fac, lmesh, elem, & ! (in)
1442 nz, vmapm, vmapp, ke_x, ke_y, cal_tend_flag )
1443
1444 implicit none
1445
1446 class(localmesh3d), intent(in) :: lmesh
1447 class(elementbase3d), intent(in) :: elem
1448 real(rp), intent(out) :: ax(elem%np,prog_vars_num,lmesh%nez)
1449 real(rp), intent(in) :: prog_vars(elem%np,prog_vars_num,lmesh%nez)
1450 real(rp), intent(in) :: prog_vars0(elem%np,prog_vars_num,lmesh%nez)
1451 real(rp), intent(in) :: dens_hyd(elem%np,lmesh%nez)
1452 real(rp), intent(in) :: pres_hyd(elem%np,lmesh%nez)
1453 class(sparsemat), intent(in) :: dz, lift
1454 real(rp), intent(in) :: impl_fac
1455 real(rp), intent(in) :: nz(elem%nfptot,lmesh%nez)
1456 integer, intent(in) :: vmapm(elem%nfptot,lmesh%nez)
1457 integer, intent(in) :: vmapp(elem%nfptot,lmesh%nez)
1458 integer, intent(in) :: ke_x, ke_y
1459 logical, intent(in) :: cal_tend_flag
1460
1461 real(rp) :: fz(elem%np), liftdelflx(elem%np)
1462 real(rp) :: del_flux(elem%nfptot,lmesh%nez,prog_vars_num)
1463 real(rp) :: rhot_hyd(elem%np)
1464 real(rp) :: pot0(elem%np), dens0(elem%np)
1465 real(rp) :: dpres(elem%np)
1466 integer :: ke_z
1467 integer :: ke
1468 real(rp) :: gamm, rgamm
1469
1470 !--------------------------------------------------------
1471
1472 gamm = cpdry/cvdry
1473 rgamm = cvdry/cpdry
1474
1475 call vi_cal_del_flux_dyn_lin( del_flux, & ! (out)
1476 prog_vars(:,ddens_vid,:), prog_vars(:,momx_vid,:), & ! (in)
1477 prog_vars(:,momy_vid ,:), prog_vars(:,momz_vid,:), & ! (in)
1478 prog_vars(:,drhot_vid,:), & ! (in)
1479 prog_vars0(:,ddens_vid,:), prog_vars0(:,momx_vid,:), & ! (in)
1480 prog_vars0(:,momy_vid ,:), prog_vars0(:,momz_vid,:), & ! (in)
1481 prog_vars0(:,drhot_vid,:), & ! (in)
1482 dens_hyd, pres_hyd, nz, vmapm, vmapp, & ! (in)
1483 lmesh, elem ) ! (in)
1484
1485
1486 !$omp parallel do private(ke, RHOT_hyd, DPRES, DENS0, POT0, Fz, LiftDelFlx)
1487 do ke_z=1, lmesh%NeZ
1488 ke = ke_x + (ke_y-1)*lmesh%NeX + (ke_z-1)*lmesh%NeX*lmesh%NeY
1489
1490 rhot_hyd(:) = pres00/rdry * (pres_hyd(:,ke_z)/pres00)**rgamm
1491
1492 dpres(:) = gamm * pres_hyd(:,ke_z) / rhot_hyd(:) &
1493 * ( 1.0_rp + prog_vars0(:,drhot_vid,ke_z) / rhot_hyd(:) )**(gamm-1) &
1494 * prog_vars(:,drhot_vid,ke_z)
1495
1496 dens0(:) = dens_hyd(:,ke_z) + prog_vars0(:,ddens_vid,ke_z)
1497 pot0(:) = ( rhot_hyd(:) + prog_vars0(:,drhot_vid,ke_z) ) / dens0(:)
1498
1499 !- DENS
1500 call sparsemat_matmul(dz, prog_vars(:,momz_vid,ke_z), fz)
1501 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke_z,ddens_vid), liftdelflx)
1502 ax(:,ddens_vid,ke_z) = lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:)
1503
1504 !- MOMX
1505 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke_z,momx_vid), liftdelflx)
1506 ax(:,momx_vid,ke_z) = liftdelflx(:)
1507
1508 !-MOMY
1509 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke_z,momy_vid), liftdelflx)
1510 ax(:,momy_vid,ke_z) = liftdelflx(:)
1511
1512 !-MOMZ
1513 call sparsemat_matmul(dz, dpres(:), fz)
1514 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke_z,momz_vid), liftdelflx)
1515 ax(:,momz_vid,ke_z) = lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) &
1516 + grav * matmul(intrpmat_vpordm1, prog_vars(:,ddens_vid,ke_z))
1517
1518 !-RHOT
1519 call sparsemat_matmul(dz, pot0(:) * prog_vars(:,momz_vid,ke_z) &
1520 + prog_vars0(:,momz_vid,ke_z) / dens0(:) * prog_vars(:,drhot_vid,ke_z) &
1521 - pot0(:) * prog_vars0(:,momz_vid,ke_z) / dens0(:) * prog_vars(:,ddens_vid,ke_z), &
1522 fz )
1523
1524 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke_z,drhot_vid), liftdelflx)
1525 ax(:,drhot_vid,ke_z) = lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:)
1526
1527 !--
1528 if ( .not. cal_tend_flag ) then
1529 ax(:,:,ke_z) = prog_vars(:,:,ke_z) + impl_fac * ax(:,:,ke_z)
1530 end if
1531
1532 end do
1533
1534
1535 return
1536 end subroutine vi_eval_ax_lin
1537
1538 subroutine vi_cal_del_flux_dyn_lin( del_flux, &
1539 DDENS_, MOMX_, MOMY_, MOMZ_, DRHOT_, &
1540 DDENS0_, MOMX0_, MOMY0_, MOMZ0_, DRHOT0_, &
1541 DENS_hyd, PRES_hyd, nz, vmapM, vmapP, lmesh, elem )
1542
1543 implicit none
1544
1545 class(localmesh3d), intent(in) :: lmesh
1546 class(elementbase3d), intent(in) :: elem
1547 real(rp), intent(out) :: del_flux(elem%nfptot*lmesh%nez,prog_vars_num)
1548 real(rp), intent(in) :: ddens_(elem%np*lmesh%nez)
1549 real(rp), intent(in) :: momx_(elem%np*lmesh%nez)
1550 real(rp), intent(in) :: momy_(elem%np*lmesh%nez)
1551 real(rp), intent(in) :: momz_(elem%np*lmesh%nez)
1552 real(rp), intent(in) :: drhot_(elem%np*lmesh%nez)
1553 real(rp), intent(in) :: ddens0_(elem%np*lmesh%nez)
1554 real(rp), intent(in) :: momx0_(elem%np*lmesh%nez)
1555 real(rp), intent(in) :: momy0_(elem%np*lmesh%nez)
1556 real(rp), intent(in) :: momz0_(elem%np*lmesh%nez)
1557 real(rp), intent(in) :: drhot0_(elem%np*lmesh%nez)
1558 real(rp), intent(in) :: dens_hyd(elem%np*lmesh%nez)
1559 real(rp), intent(in) :: pres_hyd(elem%np*lmesh%nez)
1560 real(rp), intent(in) :: nz(elem%nfptot*lmesh%nez)
1561 integer, intent(in) :: vmapm(elem%nfptot*lmesh%nez)
1562 integer, intent(in) :: vmapp(elem%nfptot*lmesh%nez)
1563
1564 integer :: i, p, ke_z, ip, im
1565 real(rp) :: alpha0, swv
1566 real(rp) :: rhot_hyd_m, rhot_hyd_p
1567 real(rp) :: dpresm, dpresp, densm, desnp, momz_p
1568 real(rp) :: pres0m, pres0p, dens0m, dens0p, pott0m, pott0p, momz0_p
1569
1570 real(rp) :: gamm, rgamm
1571 !------------------------------------------------------------------------
1572
1573 gamm = cpdry/cvdry
1574 rgamm = cvdry/cpdry
1575
1576 !$omp parallel do private( p, i, iM, iP, &
1577 !$omp rhot_hyd_M, rhot_hyd_P, dpresM, dpresP, MOMZ_P, &
1578 !$omp dens0M, dens0P, pott0M, pott0P, pres0M, pres0P, MOMZ0_P, &
1579 !$omp swV, alpha0 )
1580 do ke_z=1, lmesh%NeZ
1581 do p=1, elem%NfpTot
1582 i = p + (ke_z-1)*elem%NfpTot
1583 im = vmapm(i); ip = vmapp(i)
1584
1585 rhot_hyd_m = pres00/rdry * (pres_hyd(im)/pres00)**rgamm
1586 rhot_hyd_p = pres00/rdry * (pres_hyd(ip)/pres00)**rgamm
1587
1588 dpresm = gamm * pres_hyd(im) / rhot_hyd_m * ( 1.0_rp + drhot0_(im) / rhot_hyd_m )**(gamm-1) * drhot_(im)
1589 dpresp = gamm * pres_hyd(ip) / rhot_hyd_p * ( 1.0_rp + drhot0_(ip) / rhot_hyd_p )**(gamm-1) * drhot_(ip)
1590
1591 !-
1592 dens0m = dens_hyd(im) + ddens0_(im)
1593 dens0p = dens_hyd(ip) + ddens0_(ip)
1594
1595 pott0m = ( rhot_hyd_m + drhot0_(im) ) / dens0m
1596 pott0p = ( rhot_hyd_p + drhot0_(ip) ) / dens0p
1597
1598 pres0m = pres_hyd(im) * ( 1.0_rp + drhot0_(im) / rhot_hyd_m )**gamm
1599 pres0p = pres_hyd(ip) * ( 1.0_rp + drhot0_(ip) / rhot_hyd_p )**gamm
1600
1601 swv = nz(i)**2
1602 alpha0 = swv * max( abs(momz0_(im)/dens0m) + sqrt(gamm*pres0m/dens0m), &
1603 abs(momz0_(ip)/dens0p) + sqrt(gamm*pres0p/dens0p) )
1604
1605 if (im==ip .and. (ke_z == 1 .or. ke_z == lmesh%NeZ)) then
1606 momz_p = - momz_(im)
1607 momz0_p = - momz0_(im)
1608 else
1609 momz_p = momz_(ip)
1610 momz0_p = momz0_(ip)
1611 end if
1612
1613 del_flux(i,ddens_vid) = 0.5_rp * ( &
1614 + ( momz_p - momz_(im) ) * nz(i) &
1615 - alpha0 * ( ddens_(ip) - ddens_(im) ) )
1616
1617 del_flux(i,momx_vid) = 0.5_rp * ( &
1618 - alpha0 * ( momx_(ip) - momx_(im) ) )
1619
1620 del_flux(i,momy_vid) = 0.5_rp * ( &
1621 - alpha0 * ( momy_(ip) - momy_(im) ) )
1622
1623 del_flux(i,momz_vid) = 0.5_rp * ( &
1624 + ( dpresp - dpresm ) * nz(i) &
1625 - alpha0 * ( momz_p - momz_(im) ) )
1626
1627 del_flux(i,drhot_vid) = 0.5_rp * ( &
1628 ( pott0p * momz_p &
1629 - pott0m * momz_(im) &
1630 + momz0_p / dens0p * drhot_(ip) &
1631 - momz0_(im) / dens0m * drhot_(im) &
1632 - pott0p * momz0_p / dens0p * ddens_(ip) &
1633 + pott0m * momz0_(im) / dens0m * ddens_(im) &
1634 ) * nz(i) &
1635 - alpha0 * ( drhot_(ip) - drhot_(im) ) )
1636
1637 end do
1638 end do
1639
1640 return
1641 end subroutine vi_cal_del_flux_dyn_lin
1642
1643!-
1644 subroutine vi_construct_pmatinv( PmatDlu, PmatDlu_ipiv, PmatL, PmatU, &
1645 PROG_VARS0, DENS_hyd, PRES_hyd, & ! (in)
1646 dz, lift, impl_fac, lmesh, elem, & ! (in)
1647 nz, vmapm, vmapp, ke_x, ke_y )
1648
1649 use scale_linalgebra, only: linalgebra_lu
1650 implicit none
1651
1652 class(localmesh3d), intent(in) :: lmesh
1653 class(elementbase3d), intent(in) :: elem
1654 real(rp), intent(out) :: pmatdlu(elem%np*prog_vars_num,elem%np*prog_vars_num,lmesh%nez)
1655 integer, intent(out) :: pmatdlu_ipiv(elem%np*prog_vars_num,lmesh%nez)
1656 real(rp), intent(out) :: pmatl(elem%np,elem%np,prog_vars_num,prog_vars_num,lmesh%nez)
1657 real(rp), intent(out) :: pmatu(elem%np,elem%np,prog_vars_num,prog_vars_num,lmesh%nez)
1658 real(rp), intent(in) :: prog_vars0(elem%np,prog_vars_num,lmesh%nez)
1659 real(rp), intent(in) :: dens_hyd(elem%np,lmesh%nez)
1660 real(rp), intent(in) :: pres_hyd(elem%np,lmesh%nez)
1661 class(sparsemat), intent(in) :: dz, lift
1662 real(rp), intent(in) :: impl_fac
1663 real(rp), intent(in) :: nz(elem%nfptot,lmesh%nez)
1664 integer, intent(in) :: vmapm(elem%nfptot,lmesh%nez)
1665 integer, intent(in) :: vmapp(elem%nfptot,lmesh%nez)
1666 integer, intent(in) :: ke_x, ke_y
1667
1668 real(rp) :: rhot_hyd(elem%np)
1669 real(rp) :: pot0(elem%np,lmesh%nez)
1670 real(rp) :: dens0(elem%np,lmesh%nez)
1671 real(rp) :: pres0(elem%np,lmesh%nez)
1672 real(rp) :: w0(elem%np,lmesh%nez)
1673 real(rp) :: dpdrhot0(elem%np,lmesh%nez)
1674 integer :: ke_z, ke_z2
1675 integer :: ke, p, fp, v
1676 real(rp) :: gamm, rgamm
1677 real(rp) :: dz_p(elem%np)
1678 real(rp) :: pmatd(elem%np,prog_vars_num,elem%np,prog_vars_num)
1679
1680 integer :: f1, f2, fp_s, fp_e
1681 integer :: fmv(elem%nfp_v)
1682 integer :: fmv2(elem%nfp_v)
1683 real(rp) :: lift_op(elem%np,elem%nfptot)
1684 real(rp) :: id(elem%np,elem%np)
1685 real(rp) :: lift_(elem%np,elem%np)
1686 real(rp) :: lift_2(elem%np,elem%np)
1687 real(rp) :: lift_s(elem%np,elem%np)
1688 real(rp) :: lift_s2(elem%np,elem%np)
1689 real(rp) :: lift_p(elem%np,elem%np)
1690 real(rp) :: lift_p2(elem%np,elem%np)
1691 real(rp) :: lift_pt(elem%np,elem%np)
1692 real(rp) :: lift_pt2(elem%np,elem%np)
1693 real(rp) :: lift_pt_rho(elem%np,elem%np)
1694 real(rp) :: lift_pt_rho2(elem%np,elem%np)
1695 real(rp) :: lift_pt_rhot(elem%np,elem%np)
1696 real(rp) :: lift_pt_rhot2(elem%np,elem%np)
1697 real(rp) :: tmp(elem%nfp_v)
1698 real(rp) :: fac, facw
1699
1700 !--------------------------------------------------------
1701
1702 gamm = cpdry/cvdry
1703 rgamm = cvdry/cpdry
1704
1705 lift_op(:,:) = elem%Lift
1706
1707 id(:,:) = 0.0_rp
1708 do p=1, elem%Np
1709 id(p,p) = 1.0_rp
1710 end do
1711
1712 !$omp parallel do private(RHOT_hyd)
1713 do ke_z=1, lmesh%NeZ
1714 rhot_hyd(:) = pres00/rdry * (pres_hyd(:,ke_z)/pres00)**rgamm
1715
1716 dpdrhot0(:,ke_z) = gamm * pres_hyd(:,ke_z) / rhot_hyd(:) &
1717 * ( 1.0_rp + prog_vars0(:,drhot_vid,ke_z) / rhot_hyd(:) )**(gamm-1)
1718
1719 dens0(:,ke_z) = dens_hyd(:,ke_z) + prog_vars0(:,ddens_vid,ke_z)
1720 pot0(:,ke_z) = ( rhot_hyd(:) + prog_vars0(:,drhot_vid,ke_z) ) / dens0(:,ke_z)
1721
1722 w0(:,ke_z) = prog_vars0(:,momz_vid,ke_z) / dens0(:,ke_z)
1723 pres0(:,ke_z) = pres_hyd(:,ke_z) * ( 1.0_rp + prog_vars0(:,drhot_vid,ke_z) / rhot_hyd(:) )**gamm
1724 end do
1725
1726 !$omp parallel do private(ke, p, fp, v, f1, f2, ke_z2, dz_p, &
1727 !$omp fac, facw, tmp, lift_, lift_2, lift_s, lift_s2, lift_p, lift_p2, lift_pt, lift_pt2, lift_pt_rho, lift_pt_rho2, lift_pt_rhot, lift_pt_rhot2, &
1728 !$omp PmatD, FmV, FmV2, fp_s, fp_e)
1729 do ke_z=1, lmesh%NeZ
1730 ke = ke_x + (ke_y-1)*lmesh%NeX + (ke_z-1)*lmesh%NeX*lmesh%NeY
1731
1732 !-----
1733 pmatd(:,:,:,:) = 0.0_rp
1734 pmatl(:,:,:,:,ke_z) = 0.0_rp
1735 pmatu(:,:,:,:,ke_z) = 0.0_rp
1736
1737 do v=1, prog_vars_num
1738 pmatd(:,v,:,v) = id(:,:)
1739 end do
1740
1741 do p=1, elem%Np
1742 dz_p(:) = lmesh%Escale(p,ke,3,3) * elem%Dx3(p,:)
1743
1744 ! DDENS
1745 pmatd(p,ddens_vid,:,momz_vid) = impl_fac * dz_p(:)
1746
1747 ! MOMZ
1748 pmatd(p,momz_vid,:,ddens_vid) = impl_fac * grav * intrpmat_vpordm1(p,:)
1749 pmatd(p,momz_vid,:,drhot_vid) = impl_fac * dz_p(:) * dpdrhot0(:,ke_z)
1750
1751 !DRHOT
1752 pmatd(p,drhot_vid,:,ddens_vid) = - impl_fac * pot0(:,ke_z) * w0(:,ke_z) * dz_p(:)
1753 pmatd(p,drhot_vid,:,momz_vid) = impl_fac * pot0(:,ke_z) * dz_p(:)
1754 pmatd(p,drhot_vid,:,drhot_vid) = pmatd(p,drhot_vid,:,drhot_vid) &
1755 + impl_fac * w0(:,ke_z) * dz_p(:)
1756 end do
1757
1758 do f1=1, 2
1759 if (f1==1) then
1760 f2 = 2; ke_z2 = max(ke_z-1, 1)
1761 else
1762 f2 = 1; ke_z2 = min(ke_z+1, lmesh%NeZ)
1763 end if
1764 fac = 0.5_rp * impl_fac
1765 facw = fac
1766 if ( (ke_z == 1 .and. f1==1) .or. (ke_z == lmesh%NeZ .and. f1==elem%Nfaces_v) ) then
1767 facw = impl_fac
1768 f2 = f1
1769 end if
1770
1771 fmv(:) = elem%Fmask_v(:,f1)
1772 fmv2(:) = elem%Fmask_v(:,f2)
1773
1774 fp_s = elem%Nfp_h * elem%Nfaces_h + 1 + (f1-1)*elem%Nfp_v
1775 fp_e = fp_s + elem%Nfp_v - 1
1776
1777 !
1778 lift_(:,:) = 0.0_rp
1779 lift_2(:,:) = 0.0_rp
1780 !
1781 lift_s(:,:) = 0.0_rp
1782 lift_s2(:,:) = 0.0_rp
1783 !
1784 lift_p(:,:) = 0.0_rp
1785 lift_p2(:,:) = 0.0_rp
1786 !
1787 lift_pt(:,:) = 0.0_rp
1788 lift_pt2(:,:) = 0.0_rp
1789 !
1790 lift_pt_rho(:,:) = 0.0_rp
1791 lift_pt_rho2(:,:) = 0.0_rp
1792 !
1793 lift_pt_rhot(:,:) = 0.0_rp
1794 lift_pt_rhot2(:,:) = 0.0_rp
1795
1796 !--
1797 do fp=fp_s, fp_e
1798 p = fp-fp_s+1
1799 tmp(:) = lift_op(fmv,fp) * lmesh%Fscale(fp,ke) * nz(fp,ke_z)
1800 lift_(fmv,fmv(p)) = tmp(:)
1801 lift_2(fmv,fmv2(p)) = tmp(:)
1802 end do
1803
1804 !---
1805 do fp=fp_s, fp_e
1806 p = fp-fp_s+1
1807 tmp(:) = lift_op(fmv,fp) * lmesh%Fscale(fp,ke) * &
1808 max( abs( w0(fmv(p),ke_z ) ) + sqrt( gamm * pres0(fmv(p),ke_z ) / dens0(fmv(p),ke_z ) ), &
1809 abs( w0(fmv2(p),ke_z2) ) + sqrt( gamm * pres0(fmv2(p),ke_z2) / dens0(fmv2(p),ke_z2) ) )
1810
1811 lift_s(fmv,fmv(p)) = tmp(:)
1812 lift_s2(fmv,fmv2(p)) = tmp(:)
1813 end do
1814
1815 !--
1816 do fp=fp_s, fp_e
1817 p = fp-fp_s+1
1818 tmp(:) = lift_op(fmv,fp) * lmesh%Fscale(fp,ke) * nz(fp,ke_z)
1819 lift_p(fmv,fmv(p)) = tmp(:) * dpdrhot0(fmv(p),ke_z )
1820 lift_p2(fmv,fmv2(p)) = tmp(:) * dpdrhot0(fmv2(p),ke_z2)
1821 end do
1822
1823 !--
1824 do fp=fp_s, fp_e
1825 p = fp-fp_s+1
1826 tmp(:) = lift_op(fmv,fp) * lmesh%Fscale(fp,ke) * nz(fp,ke_z)
1827 lift_pt(fmv,fmv(p)) = tmp(:) * pot0(fmv(p),ke_z )
1828 lift_pt2(fmv,fmv2(p)) = tmp(:) * pot0(fmv2(p),ke_z2)
1829 end do
1830
1831 !--
1832
1833 do fp=fp_s, fp_e
1834 p = fp-fp_s+1
1835 tmp(:) = lift_op(fmv,fp) * lmesh%Fscale(fp,ke) * nz(fp,ke_z)
1836 lift_pt_rho(fmv,fmv(p)) = - tmp(:) * pot0(fmv(p),ke_z ) * w0(fmv(p),ke_z )
1837 lift_pt_rho2(fmv,fmv2(p)) = - tmp(:) * pot0(fmv2(p),ke_z2) * w0(fmv2(p),ke_z2)
1838 end do
1839
1840 !--
1841
1842 do fp=fp_s, fp_e
1843 p = fp-fp_s+1
1844 tmp(:) = lift_op(fmv,fp) * lmesh%Fscale(fp,ke) * nz(fp,ke_z)
1845 lift_pt_rhot(fmv,fmv(p)) = tmp(:) * w0(fmv(p),ke_z )
1846 lift_pt_rhot2(fmv,fmv2(p)) = tmp(:) * w0(fmv2(p),ke_z2)
1847 end do
1848
1849 !----
1850
1851 if ( (ke_z == 1 .and. f1==1) .or. (ke_z == lmesh%NeZ .and. f1==elem%Nfaces_v) ) then
1852 pmatd(:,ddens_vid,:,momz_vid) = pmatd(:,ddens_vid,:,momz_vid) &
1853 - facw * lift_(:,:)
1854
1855 pmatd(:,drhot_vid,:,momz_vid) = pmatd(:,drhot_vid,:,momz_vid) &
1856 - facw * lift_pt(:,:)
1857
1858 pmatd(:,momz_vid,:,momz_vid) = pmatd(:,momz_vid,:,momz_vid) &
1859 + facw * lift_s(:,:)
1860 else
1861
1862 do v=1, prog_vars_num
1863 pmatd(:,v,:,v) = pmatd(:,v,:,v) + fac * lift_s(:,:)
1864 if (f1 == 1) then
1865 pmatl(:,:,v,v,ke_z) = - fac * lift_s2(:,:)
1866 else
1867 pmatu(:,:,v,v,ke_z) = - fac * lift_s2(:,:)
1868 end if
1869 end do
1870
1871 !-
1872 pmatd(:,ddens_vid,:,momz_vid) = pmatd(:,ddens_vid,:,momz_vid) &
1873 - fac * lift_(:,:)
1874 if (f1 == 1) then
1875 pmatl(:,:,momz_vid,ddens_vid,ke_z) = fac * lift_2(:,:)
1876 else
1877 pmatu(:,:,momz_vid,ddens_vid,ke_z) = fac * lift_2(:,:)
1878 end if
1879
1880 !-
1881 pmatd(:,momz_vid,:,drhot_vid) = pmatd(:,momz_vid,:,drhot_vid) &
1882 - fac * lift_p(:,:)
1883 if (f1 == 1) then
1884 pmatl(:,:,drhot_vid,momz_vid,ke_z) = fac * lift_p2(:,:)
1885 else
1886 pmatu(:,:,drhot_vid,momz_vid,ke_z) = fac * lift_p2(:,:)
1887 end if
1888
1889 pmatd(:,drhot_vid,:,momz_vid) = pmatd(:,drhot_vid,:,momz_vid) &
1890 - fac * lift_pt(:,:)
1891 pmatd(:,drhot_vid,:,ddens_vid) = pmatd(:,drhot_vid,:,ddens_vid) &
1892 - fac * lift_pt_rho(:,:)
1893 pmatd(:,drhot_vid,:,drhot_vid) = pmatd(:,drhot_vid,:,drhot_vid) &
1894 - fac * lift_pt_rhot(:,:)
1895 if (f1 == 1) then
1896 pmatl(:,:,momz_vid ,drhot_vid,ke_z) = fac * lift_pt2(:,:)
1897 pmatl(:,:,ddens_vid,drhot_vid,ke_z) = fac * lift_pt_rho2(:,:)
1898 pmatl(:,:,drhot_vid,drhot_vid,ke_z) = pmatl(:,:,drhot_vid,drhot_vid,ke_z) &
1899 + fac * lift_pt_rhot2(:,:)
1900 else
1901 pmatu(:,:,momz_vid ,drhot_vid,ke_z) = fac * lift_pt2(:,:)
1902 pmatu(:,:,ddens_vid,drhot_vid,ke_z) = fac * lift_pt_rho2(:,:)
1903 pmatu(:,:,drhot_vid,drhot_vid,ke_z) = pmatu(:,:,drhot_vid,drhot_vid,ke_z) &
1904 + fac * lift_pt_rhot2(:,:)
1905 end if
1906
1907 end if
1908 end do
1909
1910 call prof_rapstart( 'hevi_cal_vi_cal_pinv', 3)
1911 call get_pmatd_lu( pmatdlu(:,:,ke_z), pmatd, pmatdlu_ipiv(:,ke_z), elem%Np * prog_vars_num )
1912 !call get_PmatD_inv( PmatDlu(:,:,ke_z), PmatD, PmatDlu_ipiv(:,ke_z), elem%Np * PROG_VARS_NUM )
1913 call prof_rapend( 'hevi_cal_vi_cal_pinv', 3)
1914 end do
1915
1916 return
1917
1918 contains
1919 subroutine get_pmatd_lu( pmatDlu_, pmatD_, pmatDlu_ipiv_, N)
1920 use scale_linalgebra, only: linalgebra_lu
1921 implicit none
1922 integer, intent(in) :: N
1923 real(RP), intent(out) :: pmatDlu_(N,N)
1924 real(RP), intent(in) :: pmatD_(N,N)
1925 integer, intent(out) :: pmatDlu_ipiv_(N)
1926 integer :: info
1927 !------------------------------------------
1928
1929 pmatdlu_(:,:) = pmatd_(:,:)
1930 call linalgebra_lu(pmatdlu_, pmatdlu_ipiv_)
1931 return
1932 end subroutine get_pmatd_lu
1933
1934 subroutine get_pmatd_inv( pmatDinv_, pmatD_, pmatDlu_ipiv_, N)
1935 use scale_linalgebra, only: linalgebra_inv
1936 implicit none
1937 integer, intent(in) :: N
1938 real(RP), intent(out) :: pmatDinv_(N,N)
1939 real(RP), intent(in) :: pmatD_(N,N)
1940 integer, intent(out) :: pmatDlu_ipiv_(N)
1941 !------------------------------------------
1942
1943 pmatdinv_(:,:) = linalgebra_inv(pmatd_)
1944 return
1945 end subroutine get_pmatd_inv
1946 end subroutine vi_construct_pmatinv
1947
1948
1949end module scale_atm_dyn_dgm_nonhydro3d_hevi_gmres
scale_sparsemat::sparsemat_matmul
Definition scale_sparsemat.F90:57
scale_atm_dyn_dgm_nonhydro3d_hevi_gmres
module FElib / Fluid dyn solver / Atmosphere / Regional nonhydrostatic model / HEVI
Definition scale_atm_dyn_dgm_nonhydro3d_hevi_gmres.F90:13
scale_atm_dyn_dgm_nonhydro3d_hevi_gmres::atm_dyn_dgm_nonhydro3d_hevi_final
subroutine, public atm_dyn_dgm_nonhydro3d_hevi_final()
Definition scale_atm_dyn_dgm_nonhydro3d_hevi_gmres.F90:119
scale_atm_dyn_dgm_nonhydro3d_hevi_gmres::atm_dyn_dgm_nonhydro3d_hevi_init
subroutine, public atm_dyn_dgm_nonhydro3d_hevi_init(mesh)
Definition scale_atm_dyn_dgm_nonhydro3d_hevi_gmres.F90:94
scale_atm_dyn_dgm_nonhydro3d_hevi_gmres::atm_dyn_dgm_nonhydro3d_hevi_cal_grad_diffvars
subroutine, public atm_dyn_dgm_nonhydro3d_hevi_cal_grad_diffvars(gxu_, gyu_, gzu_, gxv_, gyv_, gzv_, gxw_, gyw_, gzw_, gxpt_, gypt_, gzpt_, ddens_, momx_, momy_, momz_, drhot_, dens_hyd, pres_hyd, dx, dy, dz, lift, lmesh, elem)
Definition scale_atm_dyn_dgm_nonhydro3d_hevi_gmres.F90:636
scale_atm_dyn_dgm_nonhydro3d_hevi_gmres::atm_dyn_dgm_nonhydro3d_hevi_cal_tend
subroutine, public atm_dyn_dgm_nonhydro3d_hevi_cal_tend(dens_dt, momx_dt, momy_dt, momz_dt, rhot_dt, ddens_, momx_, momy_, momz_, drhot_, dens_hyd, pres_hyd, coriolis, gxu_, gyu_, gzu_, gxv_, gyv_, gzv_, gxw_, gyw_, gzw_, gxpt_, gypt_, gzpt_, visccoef_h, visccoef_v, diffcoef_h, diffcoef_v, dx, dy, dz, sx, sy, sz, lift, lmesh, elem, lmesh2d, elem2d)
Definition scale_atm_dyn_dgm_nonhydro3d_hevi_gmres.F90:135
scale_atm_dyn_dgm_nonhydro3d_hevi_gmres::atm_dyn_dgm_nonhydro3d_hevi_cal_vi
subroutine, public atm_dyn_dgm_nonhydro3d_hevi_cal_vi(dens_dt, momx_dt, momy_dt, momz_dt, rhot_dt, ddens_, momx_, momy_, momz_, drhot_, dens_hyd, pres_hyd, dz, lift, impl_fac, lmesh, elem, lmesh2d, elem2d)
Definition scale_atm_dyn_dgm_nonhydro3d_hevi_gmres.F90:822
scale_element_base
module FElib / Element / Base
Definition scale_element_base.F90:11
scale_element_hexahedral
module FElib / Element / hexahedron
Definition scale_element_hexahedral.F90:10
scale_gmres
Module common / GMRES.
Definition scale_gmres.F90:15
scale_linalgebra
Module common / Linear algebra.
Definition scale_linalgebra.F90:11
scale_linalgebra::linalgebra_inv
real(rp) function, dimension(size(a, 1), size(a, 2)), public linalgebra_inv(a)
Calculate a inversion of matrix A.
Definition scale_linalgebra.F90:77
scale_linalgebra::linalgebra_lu
subroutine, public linalgebra_lu(a_lu, ipiv)
Perform LU factorization.
Definition scale_linalgebra.F90:174
scale_localmesh_2d
module FElib / Mesh / Local 2D
Definition scale_localmesh_2d.F90:10
scale_localmesh_3d
module FElib / Mesh / Local 3D
Definition scale_localmesh_3d.F90:10
scale_localmeshfield_base
module FElib / Data / base
Definition scale_localmeshfield_base.F90:10
scale_mesh_base3d
module FElib / Mesh / Base 3D
Definition scale_mesh_base3d.F90:10
scale_meshfield_base
module FElib / Data / base
Definition scale_meshfield_base.F90:10
scale_sparsemat
Module common / sparsemat.
Definition scale_sparsemat.F90:11
get_pmatd_lu
subroutine get_pmatd_lu(pmatdlu_, pmatd_, pmatdlu_ipiv_, n)
Definition scale_atm_dyn_dgm_nonhydro3d_hevi_gmres.F90:1920
scale_element_base::elementbase2d
Derived type representing a 2D reference element.
Definition scale_element_base.F90:71
scale_element_base::elementbase3d
Derived type representing a 3D reference element.
Definition scale_element_base.F90:111
scale_element_hexahedral::hexahedralelement
Derived type representing a hexahedral element.
Definition scale_element_hexahedral.F90:31
scale_gmres::gmres
Derived type to provide a iterative solver for system of linear equations using GMRES.
Definition scale_gmres.F90:33
scale_localmesh_2d::localmesh2d
Derived type representing a local mesh for 2D domain.
Definition scale_localmesh_2d.F90:30
scale_localmesh_3d::localmesh3d
Derived type to manage a local 3D computational domain.
Definition scale_localmesh_3d.F90:32
scale_localmeshfield_base::localmeshfield3d
Derived type representing a field with 3D local mesh.
Definition scale_localmeshfield_base.F90:69
scale_mesh_base3d::meshbase3d
Derived type to manage a computational mesh (base type for 3D domain)
Definition scale_mesh_base3d.F90:37
scale_meshfield_base::meshfield3d
Derived type representing a field with 3D mesh.
Definition scale_meshfield_base.F90:95
scale_sparsemat::sparsemat
Derived type to manage a sparse matrix.
Definition scale_sparsemat.F90:33
Generated by doxygen 1.14.0