FE-Project
Loading...
Searching...
No Matches
scale_atm_dyn_dgm_nonhydro2d.F90
Go to the documentation of this file.
1!-------------------------------------------------------------------------------
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
36
37
38 !-----------------------------------------------------------------------------
39 implicit none
40 private
41 !-----------------------------------------------------------------------------
42 !
43 !++ Public procedures
44 !
51
52 !-----------------------------------------------------------------------------
53 !
54 !++ Public parameters & variables
55 !
56
57 !-----------------------------------------------------------------------------
58 !
59 !++ Private procedures & variables
60 !
61 !-------------------
62
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
68
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
76
77 real(RP), private, allocatable :: FilterMat(:,:)
78 real(RP), private, allocatable :: IntrpMat_VPOrdM1(:,:)
79
80 private :: cal_del_flux_dyn
81 private :: cal_del_graddiffvar
82
83contains
84 subroutine atm_dyn_dgm_nonhydro2d_init( mesh )
85
86 implicit none
87 class(meshbase2d), intent(in), target :: mesh
88
89
90 integer :: p1, p_
91 real(rp), allocatable :: invv_vpordm1(:,:)
92
93 type(elementbase2d), pointer :: elem
94 !--------------------------------------------
95
96 elem => mesh%refElem2D
97
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
103 end do
104
105 allocate( intrpmat_vpordm1(elem%Np,elem%Np) )
106 intrpmat_vpordm1(:,:) = matmul(elem%V, invv_vpordm1)
107
108 return
109 end subroutine atm_dyn_dgm_nonhydro2d_init
110
112 elem, &
113 etac, alpha, ord )
114
115 implicit none
116 class(elementbase2d), intent(in) :: elem
117 real(rp), intent(in) :: etac
118 real(rp), intent(in) :: alpha
119 integer, intent(in) :: ord
120
121 real(rp) :: filter1d(elem%nfp)
122 real(rp) :: eta
123 integer :: p1, p2
124 integer :: l
125 !----------------------------------------------------
126
127 filter1d(:) = 1.0_rp
128 do p1=1, 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 ))
132 end if
133 end do
134
135 allocate( filtermat(elem%Np,elem%Np) )
136 filtermat(:,:) = 0.0_rp
137 do p2=1, elem%Nfp
138 do p1=1, elem%Nfp
139 l = p1 + (p2-1)*elem%Nfp
140 filtermat(l,l) = filter1d(p1) * filter1d(p2)
141 end do
142 end do
143 filtermat(:,:) = matmul(filtermat, elem%invV)
144 filtermat(:,:) = matmul(elem%V, filtermat)
145
146 return
148
150 implicit none
151 !--------------------------------------------
152
153 if( allocated(intrpmat_vpordm1) ) deallocate( intrpmat_vpordm1 )
154 if( allocated(filtermat) ) deallocate( filtermat )
155
156 return
157 end subroutine atm_dyn_dgm_nonhydro2d_final
158
159 !-------------------------------
160
162 DDENS_, MOMX_, MOMZ_, DRHOT_, lmesh, elem )
163 implicit none
164
165 class(localmesh2d), intent(in) :: lmesh
166 class(elementbase2d), intent(in) :: 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)
171
172 integer :: k
173 !------------------------------------
174
175 do k=1, lmesh%Ne
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))
180 end do
181
182 return
184
186 DENS_dt, MOMX_dt, MOMZ_dt, RHOT_dt, & ! (out)
187 ddens_, momx_, momz_, drhot_, dens_hyd, pres_hyd, & ! (in)
188 gxu_, gzu_, gxw_, gzw_, gxpt_, gzpt_, & ! (in)
189 visccoef_h, visccoef_v, diffcoef_h, diffcoef_v, & ! (in)
190 dx, dz, sx, sz, lift, lmesh, elem)
191
192 implicit none
193
194 class(localmesh2d), intent(in) :: lmesh
195 class(elementbase2d), intent(in) :: 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
217
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)
222
223 integer :: k
224 !------------------------------------------------------------------------
225
226 call prof_rapstart( 'cal_dyn_tend_bndflux', 2)
227 call cal_del_flux_dyn( del_flux, & ! (out)
228 ddens_, momx_, momz_, drhot_, dens_hyd, pres_hyd, & ! (in)
229 gxu_, gzu_, gxw_, gzw_, gxpt_, gzpt_, & ! (in)
230 visccoef_h, visccoef_v, & ! (in)
231 diffcoef_h, diffcoef_v, & ! (in)
232 lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), & ! (in)
233 lmesh%vmapM, lmesh%vmapP, & ! (in)
234 lmesh, elem ) ! (in)
235 call prof_rapend( 'cal_dyn_tend_bndflux', 2)
236
237 !-----
238 call prof_rapstart( 'cal_dyn_tend_interior', 2)
239 do k = lmesh%NeS, lmesh%NeE
240 !--
241
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)
246
247 u_(:) = momx_(:,k)/dens_(:)
248 w_(:) = momz_(:,k)/dens_(:)
249
250 !-- DENS
251 call sparsemat_matmul(dx, momx_(:,k), fx)
252 call sparsemat_matmul(dz, momz_(:,k), fz)
253 call sparsemat_matmul(lift, lmesh%Fscale(:,k)*del_flux(:,k,vars_ddens_id), liftdelflx)
254
255 dens_dt(:,k) = - ( &
256 lmesh%Escale(:,k,1,1) * fx(:) &
257 + lmesh%Escale(:,k,2,2) * fz(:) &
258 + liftdelflx(:) )
259
260 !-- MOMX
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)
264
265 momx_dt(:,k) = - ( &
266 lmesh%Escale(:,k,1,1) * fx(:) &
267 + lmesh%Escale(:,k,2,2) * fz(:) &
268 + liftdelflx(:) )
269
270 !-- MOMZ
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)
274
275 momz_dt(:,k) = - ( &
276 lmesh%Escale(:,k,1,1) * fx(:) &
277 + lmesh%Escale(:,k,2,2) * fz(:) &
278 + liftdelflx(:) ) &
279 - matmul(intrpmat_vpordm1, ddens_(:,k)) * grav
280 !- DDENS_(:,k)*Grav
281
282
283 !-- RHOT
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)
287
288 rhot_dt(:,k) = - ( &
289 lmesh%Escale(:,k,1,1) * fx(:) &
290 + lmesh%Escale(:,k,2,2) * fz(:) &
291 + liftdelflx )
292
293 end do
294 call prof_rapend( 'cal_dyn_tend_interior', 2)
295
296 return
298
299 !------
300
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 )
307
308 implicit none
309
310 class(localmesh2d), intent(in) :: lmesh
311 class(elementbase2d), intent(in) :: 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)
333
334 integer :: i, ip, im
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
340 !------------------------------------------------------------------------
341
342 gamm = cpdry/cvdry
343 rgamm = cvdry/cpdry
344
345 do i=1, elem%NfpTot*lmesh%Ne
346 im = vmapm(i); ip = vmapp(i)
347
348 rhot_hyd = pres00/rdry * (pres_hyd(im)/pres00)**rgamm
349
350 rhotm = rhot_hyd + drhot_(im)
351 presm = pres00 * (rdry*rhotm/pres00)**gamm
352 dpresm = presm - pres_hyd(im)
353
354 rhotp = rhot_hyd + drhot_(ip)
355 presp = pres00 * (rdry*rhotp/pres00)**gamm
356 dpresp = presp - pres_hyd(im)
357
358 densm = ddens_(im) + dens_hyd(im)
359 densp = ddens_(ip) + dens_hyd(im)
360
361 velm = (momx_(im)*nx(i) + momz_(im)*nz(i))/densm
362 velp = (momx_(ip)*nx(i) + momz_(ip)*nz(i))/densp
363
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))
368
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) )
374
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) )
379
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) )
384 else
385 ddifffluxu = 0.0_rp
386 ddifffluxw = 0.0_rp
387 ddifffluxpt = 0.0_rp
388 end if
389
390 del_flux(i,vars_ddens_id) = 0.5_rp*( &
391 ( densp*velp - densm*velm ) &
392 - alpha*(ddens_(ip) - ddens_(im)) )
393
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)) &
398 - ddifffluxu )
399
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)) &
404 - ddifffluxw )
405
406 del_flux(i,vars_drhot_id) = 0.5_rp*( &
407 ( rhotp*velp - rhotm*velm ) &
408 - alpha*(drhot_(ip) - drhot_(im)) &
409 - ddifffluxpt )
410
411 end do
412
413 return
414 end subroutine cal_del_flux_dyn
415
416 subroutine atm_dyn_dgm_nonhydro2d_cal_grad_diffvars( GxU_, GzU_, GxW_, GzW_, GxPT_, GzPT_, &
417 DDENS_, MOMX_, MOMZ_, DRHOT_, DENS_hyd, PRES_hyd, &
418 Dx, Dz, Lift, lmesh, elem )
419
420 implicit none
421
422 class(localmesh2d), intent(in) :: lmesh
423 class(elementbase2d), intent(in) :: 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)
437
438 integer :: k
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)
442 !------------------------------------------------------------------------------
443
444 call cal_del_graddiffvar( del_flux, & ! (out)
445 ddens_, momx_, momz_, drhot_, dens_hyd, pres_hyd, & ! (in)
446 lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), & ! (in)
447 lmesh%vmapM, lmesh%vmapP, & ! (in)
448 lmesh, elem ) ! (in)
449
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)
454
455 u_(:) = momx_(:,k)/dens_(:)
456 w_(:) = momz_(:,k)/dens_(:)
457 dtheta_(:) = rhot_(:)/dens_(:) - rhot_hyd(:)/dens_hyd(:,k)
458
459 call sparsemat_matmul(dx, u_, fx)
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(:)
462
463 call sparsemat_matmul(dz, u_, fz)
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(:)
466
467 call sparsemat_matmul(dx, w_, fx)
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(:)
470
471 call sparsemat_matmul(dz, w_, fz)
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(:)
474
475 call sparsemat_matmul(dx, dtheta_, fx)
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(:)
478
479 call sparsemat_matmul(dz, dtheta_, fz)
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(:)
482 end do
483
484 return
486
487 subroutine cal_del_graddiffvar( del_flux, &
488 DDENS_, MOMX_, MOMZ_, DRHOT_, DENS_hyd, PRES_hyd, nx, nz, vmapM, vmapP, lmesh, elem )
489
490 implicit none
491
492 class(localmesh2d), intent(in) :: lmesh
493 class(elementbase2d), intent(in) :: 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)
505
506 integer :: i, ip, im
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
512 !------------------------------------------------------------------------
513
514 gamm = cpdry/cvdry
515 rgamm = cvdry/cpdry
516
517 do i=1, elem%NfpTot*lmesh%Ne
518 im = vmapm(i); ip = vmapp(i)
519
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)
525
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)
529
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)
536 end do
537
538 return
539 end subroutine cal_del_graddiffvar
540
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_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 / Mesh / Base 2D
module FElib / Data / base
module common / sparsemat