12#include "scaleFElib.h"
22 use scale_const,
only: &
25 cpdry => const_cpdry, &
26 cvdry => const_cvdry, &
63 integer,
private,
parameter :: VARS_DDENS_ID = 1
64 integer,
private,
parameter :: VARS_MOMX_ID = 2
65 integer,
private,
parameter :: VARS_MOMZ_ID = 3
66 integer,
private,
parameter :: VARS_DRHOT_ID = 4
67 integer,
private,
parameter :: PROG_VARS_NUM = 4
69 integer,
private,
parameter :: VARS_GxU_ID = 1
70 integer,
private,
parameter :: VARS_GzU_ID = 2
71 integer,
private,
parameter :: VARS_GxW_ID = 3
72 integer,
private,
parameter :: VARS_GzW_ID = 4
73 integer,
private,
parameter :: VARS_GxPT_ID = 5
74 integer,
private,
parameter :: VARS_GzPT_ID = 6
75 integer,
private,
parameter :: AUX_DIFFVARS_NUM = 6
77 real(RP),
private,
allocatable :: FilterMat(:,:)
78 real(RP),
private,
allocatable :: IntrpMat_VPOrdM1(:,:)
80 private :: cal_del_flux_dyn
81 private :: cal_del_graddiffvar
91 real(rp),
allocatable :: invv_vpordm1(:,:)
96 elem => mesh%refElem2D
98 allocate( invv_vpordm1(elem%Np,elem%Np) )
99 invv_vpordm1(:,:) = elem%invV
100 do p1=1, elem%PolyOrder+1
101 p_ = p1 + elem%PolyOrder*(elem%PolyOrder + 1)
102 invv_vpordm1(p_,:) = 0.0_rp
105 allocate( intrpmat_vpordm1(elem%Np,elem%Np) )
106 intrpmat_vpordm1(:,:) = matmul(elem%V, invv_vpordm1)
117 real(rp),
intent(in) :: etac
118 real(rp),
intent(in) :: alpha
119 integer,
intent(in) :: ord
121 real(rp) :: filter1d(elem%nfp)
129 eta = dble(p1-1)/dble(elem%PolyOrder)
130 if ( eta > etac .and. p1 /= 1)
then
131 filter1d(p1) = exp( - alpha*( ((eta - etac)/(1.0_rp - etac))**ord ))
135 allocate( filtermat(elem%Np,elem%Np) )
136 filtermat(:,:) = 0.0_rp
139 l = p1 + (p2-1)*elem%Nfp
140 filtermat(l,l) = filter1d(p1) * filter1d(p2)
143 filtermat(:,:) = matmul(filtermat, elem%invV)
144 filtermat(:,:) = matmul(elem%V, filtermat)
153 if(
allocated(intrpmat_vpordm1) )
deallocate( intrpmat_vpordm1 )
154 if(
allocated(filtermat) )
deallocate( filtermat )
162 DDENS_, MOMX_, MOMZ_, DRHOT_, lmesh, elem )
167 real(rp),
intent(inout) :: ddens_(elem%np,lmesh%nea)
168 real(rp),
intent(inout) :: momx_(elem%np,lmesh%nea)
169 real(rp),
intent(inout) :: momz_(elem%np,lmesh%nea)
170 real(rp),
intent(inout) :: drhot_(elem%np,lmesh%nea)
176 ddens_(:,k) = matmul(filtermat,ddens_(:,k))
177 momx_(:,k) = matmul(filtermat,momx_(:,k))
178 momz_(:,k) = matmul(filtermat,momz_(:,k))
179 drhot_(:,k) = matmul(filtermat,drhot_(:,k))
186 DENS_dt, MOMX_dt, MOMZ_dt, RHOT_dt, & ! (out)
187 ddens_, momx_, momz_, drhot_, dens_hyd, pres_hyd, &
188 gxu_, gzu_, gxw_, gzw_, gxpt_, gzpt_, &
189 visccoef_h, visccoef_v, diffcoef_h, diffcoef_v, &
190 dx, dz, sx, sz, lift, lmesh, elem)
196 type(
sparsemat),
intent(in) :: dx, dz, sx, sz, lift
197 real(rp),
intent(out) :: dens_dt(elem%np,lmesh%nea)
198 real(rp),
intent(out) :: momx_dt(elem%np,lmesh%nea)
199 real(rp),
intent(out) :: momz_dt(elem%np,lmesh%nea)
200 real(rp),
intent(out) :: rhot_dt(elem%np,lmesh%nea)
201 real(rp),
intent(in) :: ddens_(elem%np,lmesh%nea)
202 real(rp),
intent(in) :: momx_(elem%np,lmesh%nea)
203 real(rp),
intent(in) :: momz_(elem%np,lmesh%nea)
204 real(rp),
intent(in) :: drhot_(elem%np,lmesh%nea)
205 real(rp),
intent(in) :: dens_hyd(elem%np,lmesh%nea)
206 real(rp),
intent(in) :: pres_hyd(elem%np,lmesh%nea)
207 real(rp),
intent(in) :: gxu_(elem%np,lmesh%nea)
208 real(rp),
intent(in) :: gzu_(elem%np,lmesh%nea)
209 real(rp),
intent(in) :: gxw_(elem%np,lmesh%nea)
210 real(rp),
intent(in) :: gzw_(elem%np,lmesh%nea)
211 real(rp),
intent(in) :: gxpt_(elem%np,lmesh%nea)
212 real(rp),
intent(in) :: gzpt_(elem%np,lmesh%nea)
213 real(rp),
intent(in) :: visccoef_h
214 real(rp),
intent(in) :: visccoef_v
215 real(rp),
intent(in) :: diffcoef_h
216 real(rp),
intent(in) :: diffcoef_v
218 real(rp) :: fx(elem%np), fz(elem%np), liftdelflx(elem%np)
219 real(rp) :: del_flux(elem%nfptot,lmesh%ne,prog_vars_num)
220 real(rp) :: dens_(elem%np), rhot_(elem%np), dpres_(elem%np)
221 real(rp) :: pres_(elem%np), u_(elem%np), w_(elem%np)
226 call prof_rapstart(
'cal_dyn_tend_bndflux', 2)
227 call cal_del_flux_dyn( del_flux, &
228 ddens_, momx_, momz_, drhot_, dens_hyd, pres_hyd, &
229 gxu_, gzu_, gxw_, gzw_, gxpt_, gzpt_, &
230 visccoef_h, visccoef_v, &
231 diffcoef_h, diffcoef_v, &
232 lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), &
233 lmesh%vmapM, lmesh%vmapP, &
235 call prof_rapend(
'cal_dyn_tend_bndflux', 2)
238 call prof_rapstart(
'cal_dyn_tend_interior', 2)
239 do k = lmesh%NeS, lmesh%NeE
242 rhot_(:) = pres00/rdry * (pres_hyd(:,k)/pres00)**(cvdry/cpdry) + drhot_(:,k)
243 pres_(:) = pres00 * (rdry*rhot_(:)/pres00)**(cpdry/cvdry)
244 dpres_(:) = pres_(:) - pres_hyd(:,k)
245 dens_(:) = ddens_(:,k) + dens_hyd(:,k)
247 u_(:) = momx_(:,k)/dens_(:)
248 w_(:) = momz_(:,k)/dens_(:)
253 call sparsemat_matmul(lift, lmesh%Fscale(:,k)*del_flux(:,k,vars_ddens_id), liftdelflx)
256 lmesh%Escale(:,k,1,1) * fx(:) &
257 + lmesh%Escale(:,k,2,2) * fz(:) &
261 call sparsemat_matmul(dx, u_(:)*momx_(:,k) + dpres_(:) - visccoef_h*dens_(:)*gxu_(:,k), fx)
262 call sparsemat_matmul(dz, w_(:)*momx_(:,k) - visccoef_v*dens_(:)*gzu_(:,k) ,fz)
263 call sparsemat_matmul(lift, lmesh%Fscale(:,k)*del_flux(:,k,vars_momx_id), liftdelflx)
266 lmesh%Escale(:,k,1,1) * fx(:) &
267 + lmesh%Escale(:,k,2,2) * fz(:) &
271 call sparsemat_matmul(dx, u_(:)*momz_(:,k) - visccoef_h*dens_(:)*gxw_(:,k), fx)
272 call sparsemat_matmul(dz, w_(:)*momz_(:,k) + dpres_(:) - visccoef_v*dens_(:)*gzw_(:,k), fz)
273 call sparsemat_matmul(lift, lmesh%Fscale(:,k)*del_flux(:,k,vars_momz_id), liftdelflx)
276 lmesh%Escale(:,k,1,1) * fx(:) &
277 + lmesh%Escale(:,k,2,2) * fz(:) &
279 - matmul(intrpmat_vpordm1, ddens_(:,k)) * grav
284 call sparsemat_matmul(dx, u_(:)*rhot_(:) - diffcoef_h*dens_(:)*gxpt_(:,k), fx)
285 call sparsemat_matmul(dz, w_(:)*rhot_(:) - diffcoef_v*dens_(:)*gzpt_(:,k), fz)
286 call sparsemat_matmul(lift, lmesh%Fscale(:,k)*del_flux(:,k,vars_drhot_id), liftdelflx)
289 lmesh%Escale(:,k,1,1) * fx(:) &
290 + lmesh%Escale(:,k,2,2) * fz(:) &
294 call prof_rapend(
'cal_dyn_tend_interior', 2)
301 subroutine cal_del_flux_dyn( del_flux, &
302 DDENS_, MOMX_, MOMZ_, DRHOT_, DENS_hyd, PRES_hyd, &
303 GxU_, GzU_, GxW_, GzW_, GxPT_, GzPT_, &
304 viscCoef_h, viscCoef_v, &
305 diffCoef_h, diffCoef_v, &
306 nx, nz, vmapM, vmapP, lmesh, elem )
312 real(rp),
intent(out) :: del_flux(elem%nfptot*lmesh%ne,prog_vars_num)
313 real(rp),
intent(in) :: ddens_(elem%np*lmesh%nea)
314 real(rp),
intent(in) :: momx_(elem%np*lmesh%nea)
315 real(rp),
intent(in) :: momz_(elem%np*lmesh%nea)
316 real(rp),
intent(in) :: drhot_(elem%np*lmesh%nea)
317 real(rp),
intent(in) :: dens_hyd(elem%np*lmesh%nea)
318 real(rp),
intent(in) :: pres_hyd(elem%np*lmesh%nea)
319 real(rp),
intent(in) :: gxu_(elem%np*lmesh%nea)
320 real(rp),
intent(in) :: gzu_(elem%np*lmesh%nea)
321 real(rp),
intent(in) :: gxw_(elem%np*lmesh%nea)
322 real(rp),
intent(in) :: gzw_(elem%np*lmesh%nea)
323 real(rp),
intent(in) :: gxpt_(elem%np*lmesh%nea)
324 real(rp),
intent(in) :: gzpt_(elem%np*lmesh%nea)
325 real(rp),
intent(in) :: visccoef_h
326 real(rp),
intent(in) :: visccoef_v
327 real(rp),
intent(in) :: diffcoef_h
328 real(rp),
intent(in) :: diffcoef_v
329 real(rp),
intent(in) :: nx(elem%nfptot*lmesh%ne)
330 real(rp),
intent(in) :: nz(elem%nfptot*lmesh%ne)
331 integer,
intent(in) :: vmapm(elem%nfptot*lmesh%ne)
332 integer,
intent(in) :: vmapp(elem%nfptot*lmesh%ne)
335 real(rp) :: velp, velm, alpha
336 real(rp) :: um, up, wm, wp, presm, presp, dpresm, dpresp, densm, densp, rhotm, rhotp, rhot_hyd
337 real(rp) :: ddifffluxu, ddifffluxw, ddifffluxpt
338 real(rp) :: gamm, rgamm
339 real(rp) :: mu, visccoef, diffcoef
345 do i=1, elem%NfpTot*lmesh%Ne
346 im = vmapm(i); ip = vmapp(i)
348 rhot_hyd = pres00/rdry * (pres_hyd(im)/pres00)**rgamm
350 rhotm = rhot_hyd + drhot_(im)
351 presm = pres00 * (rdry*rhotm/pres00)**gamm
352 dpresm = presm - pres_hyd(im)
354 rhotp = rhot_hyd + drhot_(ip)
355 presp = pres00 * (rdry*rhotp/pres00)**gamm
356 dpresp = presp - pres_hyd(im)
358 densm = ddens_(im) + dens_hyd(im)
359 densp = ddens_(ip) + dens_hyd(im)
361 velm = (momx_(im)*nx(i) + momz_(im)*nz(i))/densm
362 velp = (momx_(ip)*nx(i) + momz_(ip)*nz(i))/densp
364 alpha = max( sqrt(gamm*presm/densm) + abs(velm), sqrt(gamm*presp/densp) + abs(velp) )
365 mu = (2.0_rp * dble((elem%PolyOrder+1)*(elem%PolyOrder+2)) / 2.0_rp / 600.0_rp)
366 visccoef = visccoef_h*abs(nx(i)) + visccoef_v*abs(nz(i))
367 diffcoef = diffcoef_h*abs(nx(i)) + diffcoef_v*abs(nz(i))
369 if (diffcoef > 0.0_rp)
then
370 ddifffluxu = visccoef*( &
371 (densp*gxu_(ip) - densm*gxu_(im))*nx(i) &
372 + (densp*gzu_(ip) - densm*gzu_(im))*nz(i) &
373 + mu*(densp + densm)*(momx_(ip)/densp - momx_(im)/densm) )
375 ddifffluxw = visccoef*( &
376 (densp*gxw_(ip) - densm*gxw_(im))*nx(i) &
377 + (densp*gzw_(ip) - densm*gzw_(im))*nz(i) &
378 + mu*(densp + densm)*(momz_(ip)/densp - momz_(im)/densm) )
380 ddifffluxpt = diffcoef*( &
381 (densp*gxpt_(ip) - densm*gxpt_(im))*nx(i) &
382 + (densp*gzpt_(ip) - densm*gzpt_(im))*nz(i) &
383 + mu*(densp + densm)*(rhotp/densp - rhotm/densm) )
390 del_flux(i,vars_ddens_id) = 0.5_rp*( &
391 ( densp*velp - densm*velm ) &
392 - alpha*(ddens_(ip) - ddens_(im)) )
394 del_flux(i,vars_momx_id) = 0.5_rp*( &
395 ( momx_(ip)*velp - momx_(im)*velm ) &
396 + ( dpresp - dpresm )*nx(i) &
397 - alpha*(momx_(ip) - momx_(im)) &
400 del_flux(i,vars_momz_id) = 0.5_rp*( &
401 ( momz_(ip)*velp - momz_(im)*velm) &
402 + ( dpresp - dpresm )*nz(i) &
403 - alpha*(momz_(ip) - momz_(im)) &
406 del_flux(i,vars_drhot_id) = 0.5_rp*( &
407 ( rhotp*velp - rhotm*velm ) &
408 - alpha*(drhot_(ip) - drhot_(im)) &
414 end subroutine cal_del_flux_dyn
417 DDENS_, MOMX_, MOMZ_, DRHOT_, DENS_hyd, PRES_hyd, &
418 Dx, Dz, Lift, lmesh, elem )
424 type(
sparsemat),
intent(in) :: dx, dz, lift
425 real(rp),
intent(out) :: gxu_(elem%np,lmesh%nea)
426 real(rp),
intent(out) :: gzu_(elem%np,lmesh%nea)
427 real(rp),
intent(out) :: gxw_(elem%np,lmesh%nea)
428 real(rp),
intent(out) :: gzw_(elem%np,lmesh%nea)
429 real(rp),
intent(out) :: gxpt_(elem%np,lmesh%nea)
430 real(rp),
intent(out) :: gzpt_(elem%np,lmesh%nea)
431 real(rp),
intent(in) :: ddens_(elem%np,lmesh%nea)
432 real(rp),
intent(in) :: momx_(elem%np,lmesh%nea)
433 real(rp),
intent(in) :: momz_(elem%np,lmesh%nea)
434 real(rp),
intent(in) :: drhot_(elem%np,lmesh%nea)
435 real(rp),
intent(in) :: dens_hyd(elem%np,lmesh%nea)
436 real(rp),
intent(in) :: pres_hyd(elem%np,lmesh%nea)
439 real(rp) :: dens_(elem%np), u_(elem%np), w_(elem%np), dtheta_(elem%np), rhot_(elem%np), rhot_hyd(elem%np)
440 real(rp) :: fx(elem%np), fz(elem%np), liftdelflx(elem%np)
441 real(rp) :: del_flux(elem%nfptot,lmesh%ne,aux_diffvars_num)
444 call cal_del_graddiffvar( del_flux, &
445 ddens_, momx_, momz_, drhot_, dens_hyd, pres_hyd, &
446 lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), &
447 lmesh%vmapM, lmesh%vmapP, &
450 do k=lmesh%NeS, lmesh%NeE
451 dens_(:) = ddens_(:,k) + dens_hyd(:,k)
452 rhot_hyd(:) = pres00/rdry * (pres_hyd(:,k)/pres00)**(cvdry/cpdry)
453 rhot_(:) = rhot_hyd(:) + drhot_(:,k)
455 u_(:) = momx_(:,k)/dens_(:)
456 w_(:) = momz_(:,k)/dens_(:)
457 dtheta_(:) = rhot_(:)/dens_(:) - rhot_hyd(:)/dens_hyd(:,k)
460 call sparsemat_matmul(lift, lmesh%Fscale(:,k)*del_flux(:,k,vars_gxu_id), liftdelflx)
461 gxu_(:,k) = lmesh%Escale(:,k,1,1)*fx(:) + liftdelflx(:)
464 call sparsemat_matmul(lift, lmesh%Fscale(:,k)*del_flux(:,k,vars_gzu_id), liftdelflx)
465 gzu_(:,k) = lmesh%Escale(:,k,2,2)*fz(:) + liftdelflx(:)
468 call sparsemat_matmul(lift, lmesh%Fscale(:,k)*del_flux(:,k,vars_gxw_id), liftdelflx)
469 gxw_(:,k) = lmesh%Escale(:,k,1,1)*fx(:) + liftdelflx(:)
472 call sparsemat_matmul(lift, lmesh%Fscale(:,k)*del_flux(:,k,vars_gzw_id), liftdelflx)
473 gzw_(:,k) = lmesh%Escale(:,k,2,2)*fz(:) + liftdelflx(:)
476 call sparsemat_matmul(lift, lmesh%Fscale(:,k)*del_flux(:,k,vars_gxpt_id), liftdelflx)
477 gxpt_(:,k) = lmesh%Escale(:,k,1,1)*fx(:) + liftdelflx(:)
480 call sparsemat_matmul(lift, lmesh%Fscale(:,k)*del_flux(:,k,vars_gzpt_id), liftdelflx)
481 gzpt_(:,k) = lmesh%Escale(:,k,2,2)*fz(:) + liftdelflx(:)
487 subroutine cal_del_graddiffvar( del_flux, &
488 DDENS_, MOMX_, MOMZ_, DRHOT_, DENS_hyd, PRES_hyd, nx, nz, vmapM, vmapP, lmesh, elem )
494 real(rp),
intent(out) :: del_flux(elem%nfptot*lmesh%ne,aux_diffvars_num)
495 real(rp),
intent(in) :: ddens_(elem%np*lmesh%nea)
496 real(rp),
intent(in) :: momx_(elem%np*lmesh%nea)
497 real(rp),
intent(in) :: momz_(elem%np*lmesh%nea)
498 real(rp),
intent(in) :: drhot_(elem%np*lmesh%nea)
499 real(rp),
intent(in) :: dens_hyd(elem%np*lmesh%nea)
500 real(rp),
intent(in) :: pres_hyd(elem%np*lmesh%nea)
501 real(rp),
intent(in) :: nx(elem%nfptot*lmesh%ne)
502 real(rp),
intent(in) :: nz(elem%nfptot*lmesh%ne)
503 integer,
intent(in) :: vmapm(elem%nfptot*lmesh%ne)
504 integer,
intent(in) :: vmapp(elem%nfptot*lmesh%ne)
507 real(rp) :: velp, velm, alpha
508 real(rp) :: densm, densp, rhot_hyd, rhotm, rhotp
509 real(rp) :: delu, delw, delpt
510 real(rp) :: momz_p, momx_p
511 real(rp) :: gamm, rgamm
517 do i=1, elem%NfpTot*lmesh%Ne
518 im = vmapm(i); ip = vmapp(i)
520 rhot_hyd = pres00/rdry * (pres_hyd(im)/pres00)**rgamm
521 rhotm = rhot_hyd + drhot_(im)
522 rhotp = rhot_hyd + drhot_(ip)
523 densm = ddens_(im) + dens_hyd(im)
524 densp = ddens_(ip) + dens_hyd(im)
526 delu = 0.5_rp*(momx_(ip)/densp - momx_(im)/densm)
527 delw = 0.5_rp*(momz_(ip)/densp - momz_(im)/densm)
528 delpt = 0.5_rp*(rhotp/densp - rhotm/densm)
530 del_flux(i,vars_gxu_id) = delu * nx(i)
531 del_flux(i,vars_gzu_id) = delu * nz(i)
532 del_flux(i,vars_gxw_id) = delw * nx(i)
533 del_flux(i,vars_gzw_id) = delw * nz(i)
534 del_flux(i,vars_gxpt_id) = delpt * nx(i)
535 del_flux(i,vars_gzpt_id) = delpt * nz(i)
539 end subroutine cal_del_graddiffvar
module FElib / Fluid dyn solver / Atmosphere / 2D nonhydrostatic model / HEVE
subroutine, public atm_dyn_dgm_nonhydro2d_filter_prgvar(ddens_, momx_, momz_, drhot_, lmesh, elem)
subroutine, public atm_dyn_dgm_nonhydro2d_cal_grad_diffvars(gxu_, gzu_, gxw_, gzw_, gxpt_, gzpt_, ddens_, momx_, momz_, drhot_, dens_hyd, pres_hyd, dx, dz, lift, lmesh, elem)
subroutine, public atm_dyn_dgm_nonhydro2d_cal_tend(dens_dt, momx_dt, momz_dt, rhot_dt, ddens_, momx_, momz_, drhot_, dens_hyd, pres_hyd, gxu_, gzu_, gxw_, gzw_, gxpt_, gzpt_, visccoef_h, visccoef_v, diffcoef_h, diffcoef_v, dx, dz, sx, sz, lift, lmesh, elem)
subroutine, public atm_dyn_dgm_nonhydro2d_final()
subroutine, public atm_dyn_dgm_nonhydro2d_prepair_expfilter(elem, etac, alpha, ord)
subroutine, public atm_dyn_dgm_nonhydro2d_init(mesh)
module FElib / Element / Base
module FElib / Element / Quadrilateral
module FElib / Mesh / Local 2D
module FElib / Data / base
module FElib / Mesh / Base 2D
module FElib / Data / base
module common / sparsemat