18#include "scaleFElib.h"
28 use scale_const,
only: &
29 undef => const_undef8, &
62 private :: atm_phy_mp_dgm_netoutwardflux
63 private :: atm_phy_mp_dgm_precipitation_get_delflux
64 private :: atm_phy_mp_dgm_precipitation_momentum_get_delflux
82 real(rp),
intent(out) :: intweight(lcmesh%refelem3d%nfaces,lcmesh%refelem3d%nfptot)
85 real(rp),
allocatable :: intweight_lgl1dpts_h(:)
86 real(rp),
allocatable :: intweight_lgl1dpts_v(:)
87 real(rp),
allocatable :: intweight_h(:)
88 real(rp),
allocatable :: intweight_v(:)
95 elem => lcmesh%refElem3D
96 intweight(:,:) = 0.0_rp
98 allocate( intweight_lgl1dpts_h(elem%Nnode_h1D) )
99 allocate( intweight_lgl1dpts_v(elem%Nnode_v) )
100 allocate( intweight_h(elem%Nnode_h1D*elem%Nnode_v) )
101 allocate( intweight_v(elem%Nnode_h1D**2) )
106 do f=1, elem%Nfaces_h
108 do i=1, elem%Nnode_h1D
109 l = i + (k-1)*elem%Nnode_h1D
110 intweight_h(l) = intweight_lgl1dpts_h(i) * intweight_lgl1dpts_v(k)
114 is = (f-1)*elem%Nfp_h + 1
115 ie = is + elem%Nfp_h - 1
116 intweight(f,is:ie) = intweight_h(:)
119 do f=1, elem%Nfaces_v
120 do j=1, elem%Nnode_h1D
121 do i=1, elem%Nnode_h1D
122 l = i + (j-1)*elem%Nnode_h1D
123 intweight_v(l) = intweight_lgl1dpts_h(i) * intweight_lgl1dpts_h(j)
127 is = elem%Nfaces_h*elem%Nfp_h + (f-1)*elem%Nfp_v + 1
128 ie = is + elem%Nfp_v - 1
129 intweight(elem%Nfaces_h+f,is:ie) = intweight_v(:)
137 DENS, RHOQ, CPtot, CVtot, RHOE, & ! (inout)
138 flx_hydro, sflx_rain, sflx_snow, esflx, &
139 temp, vterm, dt, rnstep, &
140 dz, lift, nz, vmapm, vmapp, intweight, &
141 qha, qla, qia, lcmesh, elem )
143 use scale_atmos_hydrometeor,
only: &
153 integer,
intent(in) :: qha
154 real(rp),
intent(inout) :: dens (elem%np,lcmesh%nez,lcmesh%ne2d)
155 real(rp),
intent(inout) :: rhoq (elem%np,lcmesh%nez,lcmesh%ne2d,qha)
156 real(rp),
intent(inout) :: cptot(elem%np,lcmesh%nez,lcmesh%ne2d)
157 real(rp),
intent(inout) :: cvtot(elem%np,lcmesh%nez,lcmesh%ne2d)
158 real(rp),
intent(inout) :: rhoe (elem%np,lcmesh%nez,lcmesh%ne2d)
159 real(rp),
intent(inout) :: flx_hydro(elem%np,lcmesh%nez,lcmesh%ne2d)
160 real(rp),
intent(inout) :: sflx_rain(elem%nfp_v,lcmesh%ne2da)
161 real(rp),
intent(inout) :: sflx_snow(elem%nfp_v,lcmesh%ne2da)
162 real(rp),
intent(inout) :: esflx (elem%nfp_v,lcmesh%ne2da)
163 real(rp),
intent(in) :: temp (elem%np,lcmesh%nez,lcmesh%ne2d)
164 real(rp),
intent(in) :: vterm(elem%np,lcmesh%nez,lcmesh%ne2d,qha)
165 real(rp),
intent(in) :: dt
166 real(rp),
intent(in) :: rnstep
169 real(rp),
intent(in) :: nz(elem%nfptot,lcmesh%nez,lcmesh%ne2d)
170 integer,
intent(in) :: vmapm(elem%nfptot,lcmesh%nez)
171 integer,
intent(in) :: vmapp(elem%nfptot,lcmesh%nez)
172 real(rp),
intent(in) :: intweight(elem%nfaces,elem%nfptot)
173 integer,
intent(in) :: qla, qia
175 real(rp) :: qflx(elem%np)
176 real(rp) :: eflx(elem%np)
177 real(rp) :: dens0(elem%np,lcmesh%nez,lcmesh%ne2d)
178 real(rp) :: rhocp(elem%np,lcmesh%nez,lcmesh%ne2d)
179 real(rp) :: rhocv(elem%np,lcmesh%nez,lcmesh%ne2d)
180 real(rp) :: ndcoefeuler(elem%np,lcmesh%nez,lcmesh%ne2d)
181 real(rp) :: dzrhoq(elem%np,lcmesh%nez,lcmesh%ne2d)
182 real(rp) :: dzrhoe(elem%np,lcmesh%nez,lcmesh%ne2d)
183 real(rp) :: ddens(elem%np)
187 real(rp) :: fct_coef(elem%np,lcmesh%nez,lcmesh%ne2d)
188 real(rp) :: rhoq0, rhoq1, rhoq_tmp(elem%np)
189 real(rp) :: netoutwardflux(lcmesh%nez,lcmesh%ne2d)
190 real(rp) :: del_flux(elem%nfptot,lcmesh%nez,lcmesh%ne2d,2)
192 real(rp) :: fz(elem%np), liftdelflx(elem%np)
193 real(rp) :: rhoq_save(elem%np)
205 if ( iq > qla + qia )
then
208 else if ( iq > qla )
then
218 do ke2d = 1, lcmesh%Ne2D
219 do ke_z = 1, lcmesh%NeZ
220 dens0(:,ke_z,ke2d) = dens(:,ke_z,ke2d)
221 rhocp(:,ke_z,ke2d) = cptot(:,ke_z,ke2d) * dens(:,ke_z,ke2d)
222 rhocv(:,ke_z,ke2d) = cvtot(:,ke_z,ke2d) * dens(:,ke_z,ke2d)
227 call atm_phy_mp_dgm_precipitation_get_delflux_dq( &
229 dens0(:,:,:), rhoq(:,:,:,iq), temp(:,:,:), cv(iq), nz(:,:,:), vmapm(:,:), vmapp(:,:), &
234 do ke2d = 1, lcmesh%Ne2D
235 do ke_z = 1, lcmesh%NeZ
236 ke = ke2d + (ke_z-1)*lcmesh%Ne2D
237 delz = ( lcmesh%pos_ev(lcmesh%EToV(ke,5),3) - lcmesh%pos_ev(lcmesh%EToV(ke,1),3) ) / dble( elem%Nnode_v )
240 call sparsemat_matmul( lift, lcmesh%Fscale(:,ke) * del_flux(:,ke_z,ke2d,1), liftdelflx )
241 dzrhoq(:,ke_z,ke2d) = lcmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:)
243 call sparsemat_matmul( dz, rhoq(:,ke_z,ke2d,iq) * cv(iq) * temp(:,ke_z,ke2d), fz )
244 call sparsemat_matmul( lift, lcmesh%Fscale(:,ke) * del_flux(:,ke_z,ke2d,2), liftdelflx )
245 dzrhoe(:,ke_z,ke2d) = lcmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:)
247 ndcoefeuler(:,ke_z,ke2d) = 0.5_rp * delz * abs(vterm(:,ke_z,ke2d,iq))
251 call atm_phy_mp_dgm_netoutwardflux( &
252 netoutwardflux(:,:), &
253 rhoq(:,:,:,iq), vterm(:,:,:,iq), dzrhoq(:,:,:), ndcoefeuler(:,:,:), &
254 lcmesh%J(:,:), lcmesh%Fscale(:,:), &
255 nz(:,:,:), vmapm(:,:), vmapp(:,:), lcmesh%VMapM(:,:), intweight(:,:), &
259 do ke2d = 1, lcmesh%Ne2D
260 do ke_z = 1, lcmesh%NeZ
261 ke = ke2d + (ke_z-1)*lcmesh%Ne2D
263 q = sum( lcmesh%J(:,ke) * elem%IntWeight_lgl(:) * rhoq(:,ke_z,ke2d,iq) ) / dt
264 fct_coef(:,ke_z,ke2d) = max( 0.0_rp, min( 1.0_rp, q / ( netoutwardflux(ke_z,ke2d) + 1.0e-10_rp ) ) )
268 call atm_phy_mp_dgm_precipitation_get_delflux( &
270 dens0(:,:,:), rhoq(:,:,:,iq), temp(:,:,:), vterm(:,:,:,iq), &
271 dzrhoq(:,:,:), dzrhoe(:,:,:), ndcoefeuler(:,:,:), &
273 cv(iq), lcmesh%J(:,:), lcmesh%Fscale(:,:), nz(:,:,:), &
274 vmapm(:,:), vmapp(:,:), lcmesh%vmapM(:,:), intweight(:,:), &
282 do ke2d = 1, lcmesh%Ne2D
283 do ke_z = 1, lcmesh%NeZ
284 ke = ke2d + (ke_z-1)*lcmesh%Ne2D
288 rhoq_save(:) = rhoq(:,ke_z,ke2d,iq)
289 qflx(:) = vterm(:,ke_z,ke2d,iq) * rhoq(:,ke_z,ke2d,iq) &
290 - ndcoefeuler(:,ke_z,ke2d) * dzrhoq(:,ke_z,ke2d)
293 call sparsemat_matmul( lift, lcmesh%Fscale(:,ke) * del_flux(:,ke_z,ke2d,1), liftdelflx )
295 ddens(:) = - dt * ( &
296 lcmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) )
297 rhoq_tmp(:) = max( 0.0_rp, rhoq(:,ke_z,ke2d,iq) + ddens(:) )
301 rhoq0 = sum( lcmesh%Gsqrt(:,ke) * lcmesh%J(:,ke) * elem%IntWeight_lgl(:) * ( rhoq(:,ke_z,ke2d,iq) + ddens(:) ) )
302 rhoq1 = sum( lcmesh%Gsqrt(:,ke) * lcmesh%J(:,ke) * elem%IntWeight_lgl(:) * rhoq_tmp(:) )
304 ddens(:) = rhoq0 / ( rhoq1 + 1.0e-32_rp ) * rhoq_tmp(:) &
305 - rhoq(:,ke_z,ke2d,iq)
306 rhoq(:,ke_z,ke2d,iq) = rhoq(:,ke_z,ke2d,iq) + ddens(:)
310 if ( iq > qla + qia ) cycle
312 flx_hydro(:,ke_z,ke2d) = flx_hydro(:,ke_z,ke2d) &
314 if ( ke_z == 1 )
then
316 sflx_snow(:,ke2d) = sflx_snow(:,ke2d) &
317 + qflx(elem%Hslice(:,1)) * rnstep
319 sflx_rain(:,ke2d) = sflx_rain(:,ke2d) &
320 + qflx(elem%Hslice(:,1)) * rnstep
326 rhocp(:,ke_z,ke2d) = rhocp(:,ke_z,ke2d) + cp(iq) * ddens(:)
327 rhocv(:,ke_z,ke2d) = rhocv(:,ke_z,ke2d) + cv(iq) * ddens(:)
328 dens(:,ke_z,ke2d) = dens(:,ke_z,ke2d) + ddens(:)
332 eflx(:) = vterm(:,ke_z,ke2d,iq) * rhoq_save(:) * temp(:,ke_z,ke2d) * cv(iq) &
333 - ndcoefeuler(:,ke_z,ke2d) * dzrhoe(:,ke_z,ke2d)
336 call sparsemat_matmul( lift, lcmesh%Fscale(:,ke) * del_flux(:,ke_z,ke2d,2), liftdelflx )
338 rhoe(:,ke_z,ke2d) = rhoe(:,ke_z,ke2d) - dt * ( &
339 + lcmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) &
342 if ( ke_z == 1 )
then
343 esflx(:,ke2d) = esflx(:,ke2d) &
344 + eflx(elem%Hslice(:,1)) * rnstep
352 do ke2d = 1, lcmesh%Ne2D
353 do ke_z = 1, lcmesh%NeZ
354 cptot(:,ke_z,ke2d) = rhocp(:,ke_z,ke2d) / dens(:,ke_z,ke2d)
355 cvtot(:,ke_z,ke2d) = rhocv(:,ke_z,ke2d) / dens(:,ke_z,ke2d)
364 MOMU_t, MOMV_t, MOMZ_t, & ! (out)
365 dens, momu, momv, momz, mflx, &
366 dz, lift, nz, vmapm, vmapp, &
372 real(rp),
intent(out) :: momu_t(elem%np,lcmesh%nea)
373 real(rp),
intent(out) :: momv_t(elem%np,lcmesh%nea)
374 real(rp),
intent(out) :: momz_t(elem%np,lcmesh%nea)
375 real(rp),
intent(in) :: dens(elem%np,lcmesh%nez,lcmesh%ne2d)
376 real(rp),
intent(in) :: momu(elem%np,lcmesh%nez,lcmesh%ne2d)
377 real(rp),
intent(in) :: momv(elem%np,lcmesh%nez,lcmesh%ne2d)
378 real(rp),
intent(in) :: momz(elem%np,lcmesh%nez,lcmesh%ne2d)
379 real(rp),
intent(in) :: mflx(elem%np,lcmesh%nez,lcmesh%ne2d)
382 real(rp),
intent(in) :: nz(elem%nfptot,lcmesh%nez,lcmesh%ne2d)
383 integer,
intent(in) :: vmapm(elem%nfptot,lcmesh%nez)
384 integer,
intent(in) :: vmapp(elem%nfptot,lcmesh%nez)
390 real(rp) :: fz(elem%np), liftdelflx(elem%np)
391 real(rp) :: del_flux(elem%nfptot,lcmesh%nez,lcmesh%ne2d,3)
393 real(rp) :: rdens(elem%np)
396 call atm_phy_mp_dgm_precipitation_momentum_get_delflux( &
398 dens(:,:,:), momu(:,:,:), momv(:,:,:), momz(:,:,:), &
400 nz(:,:,:), vmapm(:,:), vmapp(:,:), &
405 do ke2d = 1, lcmesh%Ne2D
406 do ke_z = 1, lcmesh%NeZ
407 ke = ke2d + (ke_z-1)*lcmesh%Ne2D
408 rdens(:) = 1.0_rp / dens(:,ke_z,ke2d)
410 call sparsemat_matmul( dz, mflx(:,ke_z,ke2d) * momu(:,ke_z,ke2d) * rdens(:), fz )
411 call sparsemat_matmul( lift, lcmesh%Fscale(:,ke) * del_flux(:,ke_z,ke2d,1), liftdelflx )
412 momu_t(:,ke) = - ( lcmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) )
414 call sparsemat_matmul( dz, mflx(:,ke_z,ke2d) * momv(:,ke_z,ke2d) * rdens(:), fz )
415 call sparsemat_matmul( lift, lcmesh%Fscale(:,ke) * del_flux(:,ke_z,ke2d,2), liftdelflx )
416 momv_t(:,ke) = - ( lcmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) )
418 call sparsemat_matmul( dz, mflx(:,ke_z,ke2d) * momz(:,ke_z,ke2d) * rdens(:), fz )
419 call sparsemat_matmul( lift, lcmesh%Fscale(:,ke) * del_flux(:,ke_z,ke2d,3), liftdelflx )
420 momz_t(:,ke) = - ( lcmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) )
431 CVtot, CPtot, Rtot, &
432 DENS_hyd, PRES_hyd, &
433 dt, lmesh, elem, QA, QLA, QIA, &
436 use scale_const,
only: &
437 cvdry => const_cvdry, &
438 cpdry => const_cpdry, &
441 use scale_tracer,
only: &
442 tracer_mass, tracer_r, tracer_cv, tracer_cp
443 use scale_atmos_thermodyn,
only: &
444 atmos_thermodyn_specific_heat
450 integer,
intent(in) :: qa
452 real(rp),
intent(inout) :: ddens(elem%np,lmesh%nea)
453 real(rp),
intent(inout) :: pres(elem%np,lmesh%nea)
454 real(rp),
intent(inout) :: cvtot(elem%np,lmesh%nea)
455 real(rp),
intent(inout) :: cptot(elem%np,lmesh%nea)
456 real(rp),
intent(inout) :: rtot(elem%np,lmesh%nea)
457 real(rp),
intent(in) :: dens_hyd(elem%np,lmesh%nea)
458 real(rp),
intent(in) :: pres_hyd(elem%np,lmesh%nea)
459 real(rp),
intent(in) :: dt
460 integer,
intent(in) :: qla, qia
461 real(rp),
intent(inout),
optional :: drhot(elem%np,lmesh%nea)
466 real(rp) :: int_w(elem%np)
468 real(rp) :: dens(elem%np)
469 real(rp) :: ddens0(elem%np)
471 real(rp) :: trcmass0(elem%np), trcmass1(elem%np,qa)
472 real(rp) :: mass0_elem, mass1_elem
473 real(rp) :: inten0_elem, inten_elem
475 real(rp) :: qtrc_tmp(elem%np,qa), qdry(elem%np)
476 real(rp) :: cvtot_old(elem%np), cptot_old(elem%np), rtot_old(elem%np)
477 real(rp) :: internalen(elem%np), internalen0(elem%np), temp(elem%np)
478 real(rp) :: rhot_hyd(elem%np)
481 real(rp),
parameter :: trc_eps = 1e-32_rp
483 real(rp),
parameter :: trc_eps = 1e-128_rp
492 do ke = lmesh%NeS, lmesh%NeE
495 qtrc_tmp(:,iq) = qtrc(iq)%ptr%val(:,ke)
497 call atmos_thermodyn_specific_heat( &
498 elem%Np, 1, elem%Np, qa, &
499 qtrc_tmp, tracer_mass(:), tracer_r(:), tracer_cv(:), tracer_cp(:), &
500 qdry, rtot_old, cvtot_old, cptot_old )
502 dens(:) = dens_hyd(:,ke) + ddens(:,ke)
503 ddens0(:) = ddens(:,ke)
510 internalen0(:) = cvtot(:,ke) * pres(:,ke) / rtot(:,ke)
512 temp(:) = internalen0(:) / ( dens(:) * cvtot(:,ke) )
513 internalen(:) = internalen0(:)
515 int_w(:) = lmesh%Gsqrt(:,ke) * lmesh%J(:,ke) * elem%IntWeight_lgl(:)
516 do iq = 1, 1 + qla + qia
517 trcmass0(:) = dens(:) * qtrc_tmp(:,iq)
518 trcmass1(:,iq) = max( trc_eps, trcmass0(:) )
520 mass0_elem = sum( int_w(:) * trcmass0(:) )
521 mass1_elem = sum( int_w(:) * trcmass1(:,iq) )
522 trcmass1(:,iq) = max(mass0_elem, 0.0e0_rp) / mass1_elem * trcmass1(:,iq)
525 ddens(:,ke) = ddens(:,ke) + ( trcmass1(:,iq) - trcmass0(:) )
526 internalen(:) = internalen(:) + ( trcmass1(:,iq) - trcmass0(:) ) * tracer_cv(iq) * temp(:)
531 dens(:) = dens_hyd(:,ke) + ddens(:,ke)
533 qtrc_tmp(:,iq) = trcmass1(:,iq) / dens(:)
534 qtrc(iq)%ptr%val(:,ke) = qtrc_tmp(:,iq)
537 inten0_elem = sum( int_w(:) * internalen0(:) )
538 inten_elem = sum( int_w(:) * internalen(:) )
539 internalen(:) = inten0_elem / inten_elem * internalen(:)
541 call atmos_thermodyn_specific_heat( &
542 elem%Np, 1, elem%Np, qa, &
543 qtrc_tmp, tracer_mass(:), tracer_r(:), tracer_cv(:), tracer_cp(:), &
544 qdry, rtot(:,ke), cvtot(:,ke), cptot(:,ke) )
546 internalen(:) = internalen(:) - ( ddens(:,ke) - ddens0(:) ) * grav * lmesh%zlev(:,ke)
547 pres(:,ke) = internalen(:) * rtot(:,ke) / cvtot(:,ke)
549 if (
present(drhot) )
then
550 drhot(:,ke) = pres00 / rtot(:,ke) * ( pres(:,ke) / pres00 )**( cvtot(:,ke) / cptot(:,ke) ) &
551 - pres00 / rdry * ( pres_hyd(:,ke) / pres00 )**( cvdry / cpdry )
561 subroutine atm_phy_mp_dgm_netoutwardflux( &
564 DzRHOQ_, NDcoefEuler_, &
566 nz, vmapM, vmapP, vmapM3D, IntWeight, &
572 real(rp),
intent(out) :: net_outward_flux(lmesh%nez,lmesh%ne2d)
573 real(rp),
intent(in) :: rhoq_(elem%np*lmesh%nez,lmesh%ne2d)
574 real(rp),
intent(in) :: vterm_(elem%np*lmesh%nez,lmesh%ne2d)
575 real(rp),
intent(in) :: dzrhoq_(elem%np*lmesh%nez,lmesh%ne2d)
576 real(rp),
intent(in) :: ndcoefeuler_(elem%np*lmesh%nez,lmesh%ne2d)
577 real(rp),
intent(in) :: j(elem%np*lmesh%ne)
578 real(rp),
intent(in) :: fscale(elem%nfptot,lmesh%ne)
579 real(rp),
intent(in) :: nz(elem%nfptot,lmesh%nez,lmesh%ne2d)
580 integer,
intent(in) :: vmapm(elem%nfptot,lmesh%nez)
581 integer,
intent(in) :: vmapp(elem%nfptot,lmesh%nez)
582 integer,
intent(in) :: vmapm3d(elem%nfptot,lmesh%ne)
583 real(rp),
intent(in) :: intweight(elem%nfaces,elem%nfptot)
585 real(rp) :: numflux(elem%nfptot)
586 real(rp) :: outward_flux_tmp(elem%nfaces)
587 real(rp) :: alpha(elem%nfptot)
588 real(rp) :: velm(elem%nfptot), velp(elem%nfptot)
589 real(rp) :: rhoq_m(elem%nfptot), rhoq_p(elem%nfptot)
592 integer :: ke_z, ke2d
593 integer :: ip(elem%nfptot), im(elem%nfptot)
594 integer :: im3d(elem%nfptot)
601 do ke2d=1, lmesh%Ne2D
603 ke = ke2d + (ke_z-1)*lmesh%Ne2D
605 im3d(:) = vmapm3d(:,ke)
606 im(:) = vmapm(:,ke_z); ip(:) = vmapp(:,ke_z)
608 rhoq_m(:) = rhoq_(im(:),ke2d)
609 rhoq_p(:) = rhoq_(ip(:),ke2d)
610 velm(:) = vterm_(im(:),ke2d) * nz(:,ke_z,ke2d)
611 velp(:) = vterm_(ip(:),ke2d) * nz(:,ke_z,ke2d)
612 alpha(:) = nz(:,ke_z,ke2d)**2 * max( abs(velm(:)), abs(velp(:)) )
614 where (nz(:,ke_z,ke2d) > 1.0e-10 .and. ip(:) == im(:) )
618 numflux(:) = 0.5_rp * ( rhoq_p(:) * velp(:) + rhoq_m(:) * velm(:) &
619 - ( ndcoefeuler_(ip(:),ke2d) * dzrhoq_(ip(:),ke2d) + ndcoefeuler_(im(:),ke2d) * dzrhoq_(im(:),ke2d) ) * nz(:,ke_z,ke2d) &
620 - alpha(:) * ( rhoq_p(:) - rhoq_m(:) ) )
622 outward_flux_tmp(:) = matmul( intweight(:,:), j(im3d(:)) * fscale(:,ke) * numflux(:) )
623 net_outward_flux(ke_z,ke2d) = sum( max( 0.0_rp, outward_flux_tmp(:) ) )
628 end subroutine atm_phy_mp_dgm_netoutwardflux
631 subroutine atm_phy_mp_dgm_precipitation_get_delflux_dq( &
633 DENS_, RHOQ_,TEMP_, CV, nz, vmapM, vmapP, &
640 real(rp),
intent(out) :: del_flux(elem%nfptot,lmesh%nez,lmesh%ne2d,2)
641 real(rp),
intent(in) :: dens_(elem%np*lmesh%nez,lmesh%ne2d)
642 real(rp),
intent(in) :: rhoq_(elem%np*lmesh%nez,lmesh%ne2d)
643 real(rp),
intent(in) :: temp_(elem%np*lmesh%nez,lmesh%ne2d)
644 real(rp),
intent(in) :: cv
645 real(rp),
intent(in) :: nz(elem%nfptot,lmesh%nez,lmesh%ne2d)
646 integer,
intent(in) :: vmapm(elem%nfptot,lmesh%nez)
647 integer,
intent(in) :: vmapp(elem%nfptot,lmesh%nez)
650 integer :: ke_z, ke2d
651 real(rp) :: rhoq_p(elem%nfptot), rhoq_m(elem%nfptot)
652 integer :: im(elem%nfptot), ip(elem%nfptot)
657 do ke2d=1, lmesh%Ne2D
659 ke = ke2d + (ke_z-1)*lmesh%Ne2D
660 im(:) = vmapm(:,ke_z); ip(:) = vmapp(:,ke_z)
662 rhoq_m(:) = rhoq_(im(:),ke2d)
663 rhoq_p(:) = rhoq_(ip(:),ke2d)
664 del_flux(:,ke_z,ke2d,1) = 0.5_rp * ( rhoq_p(:) - rhoq_m(:) ) * nz(:,ke_z,ke2d)
665 del_flux(:,ke_z,ke2d,2) = 0.5_rp * cv * ( rhoq_p(:) * temp_(ip(:),ke2d) - rhoq_m(:) * temp_(im(:),ke2d) ) * nz(:,ke_z,ke2d)
670 end subroutine atm_phy_mp_dgm_precipitation_get_delflux_dq
674 subroutine atm_phy_mp_dgm_precipitation_get_delflux( &
676 DENS_, RHOQ_, TEMP_, vterm_, &
677 DzRHOQ_, DzRHOE_, NDcoefEuler_, &
679 J, Fscale, nz, vmapM, vmapP, vmapM3D, IntWeight, &
686 real(rp),
intent(out) :: del_flux(elem%nfptot,lmesh%nez,lmesh%ne2d,2)
687 real(rp),
intent(in) :: dens_(elem%np*lmesh%nez,lmesh%ne2d)
688 real(rp),
intent(in) :: rhoq_(elem%np*lmesh%nez,lmesh%ne2d)
689 real(rp),
intent(in) :: temp_(elem%np*lmesh%nez,lmesh%ne2d)
690 real(rp),
intent(in) :: vterm_(elem%np*lmesh%nez,lmesh%ne2d)
691 real(rp),
intent(in) :: dzrhoq_(elem%np*lmesh%nez,lmesh%ne2d)
692 real(rp),
intent(in) :: dzrhoe_(elem%np*lmesh%nez,lmesh%ne2d)
693 real(rp),
intent(in) :: ndcoefeuler_(elem%np*lmesh%nez,lmesh%ne2d)
694 real(rp),
intent(in) :: fct_coef_(elem%np*lmesh%nez,lmesh%ne2d)
695 real(rp),
intent(in) :: cv
696 real(rp),
intent(in) :: j(elem%np*lmesh%ne)
697 real(rp),
intent(in) :: fscale(elem%nfptot,lmesh%ne)
698 real(rp),
intent(in) :: nz(elem%nfptot,lmesh%nez,lmesh%ne2d)
699 integer,
intent(in) :: vmapm(elem%nfptot,lmesh%nez)
700 integer,
intent(in) :: vmapp(elem%nfptot,lmesh%nez)
701 integer,
intent(in) :: vmapm3d(elem%nfptot,lmesh%ne)
702 real(rp),
intent(in) :: intweight(elem%nfaces,elem%nfptot)
705 integer :: ke_z, ke2d
707 integer :: ip(elem%nfptot), im(elem%nfptot)
708 real(rp) :: alpha(elem%nfptot)
709 real(rp) :: velm(elem%nfptot), velp(elem%nfptot)
710 real(rp) :: rhoq_m(elem%nfptot), rhoq_p(elem%nfptot)
711 real(rp) :: temp_m(elem%nfptot), temp_p(elem%nfptot)
713 integer :: im3d(elem%nfptot)
714 real(rp) :: r_m(elem%nfptot), r_p(elem%nfptot)
715 real(rp) :: ndcoef_m(elem%nfptot), ndcoef_p(elem%nfptot)
716 real(rp) :: numflux (elem%nfptot)
717 real(rp) :: numflux_ei(elem%nfptot)
718 real(rp) :: outward_flux_tmp(elem%nfaces)
729 do ke2d=1, lmesh%Ne2D
731 ke = ke2d + (ke_z-1)*lmesh%Ne2D
733 im3d(:) = vmapm3d(:,ke)
734 im(:) = vmapm(:,ke_z); ip(:) = vmapp(:,ke_z)
736 r_m(:) = fct_coef_(im(:),ke2d)
737 r_p(:) = fct_coef_(ip(:),ke2d)
738 rhoq_m(:) = rhoq_(im(:),ke2d)
739 rhoq_p(:) = rhoq_(ip(:),ke2d)
740 temp_m(:) = temp_(im(:),ke2d)
741 temp_p(:) = temp_(ip(:),ke2d)
743 velm(:) = vterm_(im(:),ke2d) * nz(:,ke_z,ke2d)
744 velp(:) = vterm_(ip(:),ke2d) * nz(:,ke_z,ke2d)
745 alpha(:) = nz(:,ke_z,ke2d)**2 * max( abs(velm(:)), abs(velp(:)) )
747 where (nz(:,ke_z,ke2d) > 1.0e-10 .and. ip(:) == im(:) )
751 ndcoef_m(:) = ndcoefeuler_(im(:),ke2d)
752 ndcoef_p(:) = ndcoefeuler_(ip(:),ke2d)
754 numflux(:) = 0.5_rp * ( rhoq_p(:) * velp(:) + rhoq_m(:) * velm(:) &
755 - ( ndcoef_p(:) * dzrhoq_(ip(:),ke2d) + ndcoef_m(:) * dzrhoq_(im(:),ke2d) ) * nz(:,ke_z,ke2d) &
756 - alpha(:) * ( rhoq_p(:) - rhoq_m(:) ) )
758 numflux_ei(:) = 0.5_rp * ( cv * ( rhoq_p(:) * temp_p(:) * velp(:) + rhoq_m(:) * temp_m(:) * velm(:) ) &
759 - ( ndcoef_p(:) * dzrhoe_(ip(:),ke2d) + ndcoef_m(:) * dzrhoe_(im(:),ke2d) ) * nz(:,ke_z,ke2d) &
760 - alpha(:) * cv * ( rhoq_p(:) * temp_p(:) - rhoq_m(:) * temp_m(:) ) )
763 del_flux(:,ke_z,ke2d,1) = 0.0_rp
764 del_flux(:,ke_z,ke2d,2) = 0.0_rp
765 outward_flux_tmp(:) = matmul( intweight(:,:), j(im3d(:)) * fscale(:,ke) * numflux(:) )
766 do f=1, elem%Nfaces_v
768 fp = p + (f-1)*elem%Nfp_v + elem%Nfaces_h * elem%Nfp_h
769 r = 0.5_rp * ( r_p(fp) + r_m(fp) - ( r_p(fp) - r_m(fp) ) * sign( 1.0_rp, outward_flux_tmp(elem%Nfaces_h+f) ) )
770 del_flux(fp,ke_z,ke2d,1) = numflux(fp) * r &
771 - rhoq_m(fp) * velm(fp) &
772 + ndcoef_m(fp) * dzrhoq_(im(fp),ke2d) * nz(fp,ke_z,ke2d)
773 del_flux(fp,ke_z,ke2d,2) = numflux_ei(fp) * r &
774 - rhoq_m(fp) * velm(fp) * cv * temp_m(fp) &
775 + ndcoef_m(fp) * dzrhoe_(im(fp),ke2d) * nz(fp,ke_z,ke2d)
783 end subroutine atm_phy_mp_dgm_precipitation_get_delflux
786 subroutine atm_phy_mp_dgm_precipitation_momentum_get_delflux( &
788 dens_, momu_, momv_, momz_, mflx_, &
789 nz, vmapm, vmapp, lmesh, elem )
795 real(rp),
intent(out) :: del_flux(elem%nfptot,lmesh%nez,lmesh%ne2d,3)
796 real(rp),
intent(in) :: dens_(elem%np*lmesh%nez,lmesh%ne2d)
797 real(rp),
intent(in) :: momu_(elem%np*lmesh%nez,lmesh%ne2d)
798 real(rp),
intent(in) :: momv_(elem%np*lmesh%nez,lmesh%ne2d)
799 real(rp),
intent(in) :: momz_(elem%np*lmesh%nez,lmesh%ne2d)
800 real(rp),
intent(in) :: mflx_(elem%np*lmesh%nez,lmesh%ne2d)
801 real(rp),
intent(in) :: nz(elem%nfptot,lmesh%nez,lmesh%ne2d)
802 integer,
intent(in) :: vmapm(elem%nfptot,lmesh%nez)
803 integer,
intent(in) :: vmapp(elem%nfptot,lmesh%nez)
805 integer :: ke_z, ke2d
806 integer :: ip(elem%nfptot), im(elem%nfptot)
807 real(rp) :: alpha(elem%nfptot)
808 real(rp) :: densm(elem%nfptot), densp(elem%nfptot)
809 real(rp) :: velm(elem%nfptot), velp(elem%nfptot)
815 do ke2d=1, lmesh%Ne2D
817 im(:) = vmapm(:,ke_z); ip(:) = vmapp(:,ke_z)
819 densm(:) = dens_(im(:),ke2d)
820 densp(:) = dens_(ip(:),ke2d)
821 velm(:) = mflx_(im(:),ke2d) * nz(:,ke_z,ke2d) / densm(:)
822 velp(:) = mflx_(ip(:),ke2d) * nz(:,ke_z,ke2d) / densp(:)
823 alpha(:) = nz(:,ke_z,ke2d)**2 * max( abs(velm(:)), abs(velp(:)) )
825 del_flux(:,ke_z,ke2d,1) = 0.5_rp * ( &
826 + ( momu_(ip(:),ke2d) * velp - momu_(im(:),ke2d) * velm ) &
827 - alpha(:) * ( momu_(ip(:),ke2d) - momu_(im(:),ke2d) ) )
829 del_flux(:,ke_z,ke2d,2) = 0.5_rp * ( &
830 + ( momv_(ip(:),ke2d) * velp - momv_(im(:),ke2d) * velm ) &
831 - alpha(:) * ( momv_(ip(:),ke2d) - momv_(im(:),ke2d) ) )
833 del_flux(:,ke_z,ke2d,3) = 0.5_rp * ( &
834 + ( momz_(ip(:),ke2d) * velp - momz_(im(:),ke2d) * velm ) &
835 - alpha(:) * ( momz_(ip(:),ke2d) - momz_(im(:),ke2d) ) )
840 end subroutine atm_phy_mp_dgm_precipitation_momentum_get_delflux
module FElib / Atmosphere / Physics cloud microphysics / common
subroutine, public atm_phy_mp_dgm_common_gen_intweight(intweight, lcmesh)
subroutine, public atm_phy_mp_dgm_common_negative_fixer(qtrc, ddens, pres, cvtot, cptot, rtot, dens_hyd, pres_hyd, dt, lmesh, elem, qa, qla, qia, drhot)
subroutine, public atm_phy_mp_dgm_common_precipitation(dens, rhoq, cptot, cvtot, rhoe, flx_hydro, sflx_rain, sflx_snow, esflx, temp, vterm, dt, rnstep, dz, lift, nz, vmapm, vmapp, intweight, qha, qla, qia, lcmesh, elem)
subroutine, public atm_phy_mp_dgm_common_precipitation_momentum(momu_t, momv_t, momz_t, dens, momu, momv, momz, mflx, dz, lift, nz, vmapm, vmapp, lcmesh, elem)
module FElib / Element / Base
module FElib / Mesh / Local 3D
module FElib / Data / base
module FElib / Mesh / Base 3D
Module common / Polynomial.
real(rp) function, dimension(nord+1), public polynomial_gengausslobattoptintweight(nord)
A function to calculate the Gauss-Lobbato weights.
Module common / sparsemat.
Derived type representing a 1D reference element.
Derived type representing a 2D reference element.
Derived type representing a 3D reference element.
Derived type to manage a local 3D computational domain.
Derived type to manage a computational mesh (base type for 3D domain)
Derived type to manage a sparse matrix.