FE-Project
Loading...
Searching...
No Matches
scale_atm_phy_mp_dgm_common.F90
Go to the documentation of this file.
1!> module FElib / Atmosphere / Physics cloud microphysics / common
2!!
3!! @par Description
4!! cloud microphysics process
5!! common subroutines
6!!
7!! To preserve nonnegativity in precipitation process,
8!! a limiter proposed by Light and Durran (2016, MWR) is used
9!!
10!! @par Reference
11!! - Light and Durran 2016:
12!! Preserving Nonnegativity in Discontinuous Galerkin Approximations to Scalar Transport via Truncation and Mass Aware Rescaling (TMAR).
13!! Monthly Weather Review, 144(12), 4771–4786.
14!!
15!! @author Yuta Kawai, Team SCALE
16!!
17!-------------------------------------------------------------------------------
18#include "scaleFElib.h"
20 !-----------------------------------------------------------------------------
21 !
22 !++ Used modules
23 !
24 use scale_precision
25 use scale_io
26 use scale_prc
27 use scale_prof
28 use scale_const, only: &
29 undef => const_undef8, &
30 grav => const_grav, &
31 pres00 => const_pre00
32
34 use scale_element_base, only: &
36 use scale_localmesh_3d, only: &
39
40 !-----------------------------------------------------------------------------
41 implicit none
42 private
43 !-----------------------------------------------------------------------------
44 !
45 !++ Public type & procedure
46 !
47
52
53 !-----------------------------------------------------------------------------
54 !++ Public parameters & variables
55 !-----------------------------------------------------------------------------
56
57 !-----------------------------------------------------------------------------
58 !
59 !++ Private procedure
60 !
61
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
65
66 !-----------------------------------------------------------------------------
67 !
68 !++ Private parameters & variables
69 !
70 !-----------------------------------------------------------------------------
71
72
73contains
74
75!OCL SERIAL
77 intWeight, & ! (out)
78 lcmesh ) ! (in)
80 implicit none
81 class(localmesh3d), target :: lcmesh
82 real(rp), intent(out) :: intweight(lcmesh%refelem3d%nfaces,lcmesh%refelem3d%nfptot)
83
84 class(elementbase3d), pointer :: elem
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(:)
89
90 integer :: f
91 integer :: i, j, k, l
92 integer :: is, ie
93 !--------------------------------------------
94
95 elem => lcmesh%refElem3D
96 intweight(:,:) = 0.0_rp
97
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) )
102
103 intweight_lgl1dpts_h(:) = polynomial_gengausslobattoptintweight(elem%PolyOrder_h)
104 intweight_lgl1dpts_v(:) = polynomial_gengausslobattoptintweight(elem%PolyOrder_v)
105
106 do f=1, elem%Nfaces_h
107 do k=1, elem%Nnode_v
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)
111 end do
112 end do
113
114 is = (f-1)*elem%Nfp_h + 1
115 ie = is + elem%Nfp_h - 1
116 intweight(f,is:ie) = intweight_h(:)
117 end do
118
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)
124 end do
125 end do
126
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(:)
130 end do
131
132 return
134
135!OCL SERIAL
137 DENS, RHOQ, CPtot, CVtot, RHOE, & ! (inout)
138 flx_hydro, sflx_rain, sflx_snow, esflx, & ! (inout)
139 temp, vterm, dt, rnstep, & ! (in)
140 dz, lift, nz, vmapm, vmapp, intweight, & ! (in)
141 qha, qla, qia, lcmesh, elem ) ! (in)
142
143 use scale_atmos_hydrometeor, only: &
144 cv_water, &
145 cp_water, &
146 cv_ice, &
147 cp_ice
148
149 implicit none
150
151 class(localmesh3d), intent(in) :: lcmesh
152 class(elementbase3d), intent(in) :: elem
153 integer, intent(in) :: qha !< hydrometeor (water + ice)
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
167 type(sparsemat), intent(in) :: dz
168 type(sparsemat), intent(in) :: lift
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
174
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)
184 real(rp) :: cp(qha)
185 real(rp) :: cv(qha)
186
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)
191
192 real(rp) :: fz(elem%np), liftdelflx(elem%np)
193 real(rp) :: rhoq_save(elem%np)
194
195 integer :: ke2d
196 integer :: ke_z
197 integer :: ke
198 integer :: iq
199
200 real(rp) :: q
201 real(rp) :: delz
202 !-------------------------------------------------------
203
204 do iq = 1, qha
205 if ( iq > qla + qia ) then
206 cp(iq) = undef
207 cv(iq) = undef
208 else if ( iq > qla ) then ! ice water
209 cp(iq) = cp_ice
210 cv(iq) = cv_ice
211 else ! liquid water
212 cp(iq) = cp_water
213 cv(iq) = cv_water
214 end if
215 end do
216
217 !$omp parallel do collapse(2)
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)
223 end do
224 end do
225
226 do iq = 1, qha
227 call atm_phy_mp_dgm_precipitation_get_delflux_dq( &
228 del_flux(:,:,:,:), & ! (out)
229 dens0(:,:,:), rhoq(:,:,:,iq), temp(:,:,:), cv(iq), nz(:,:,:), vmapm(:,:), vmapp(:,:), & ! (in)
230 lcmesh, elem ) ! (in)
231
232 !$omp parallel do private( &
233 !$omp ke2D, ke_z, ke, delz, Fz, LiftDelFlx )
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 )
238
239 call sparsemat_matmul( dz, rhoq(:,ke_z,ke2d,iq), fz )
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(:)
242
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(:)
246
247 ndcoefeuler(:,ke_z,ke2d) = 0.5_rp * delz * abs(vterm(:,ke_z,ke2d,iq))
248 end do
249 end do
250
251 call atm_phy_mp_dgm_netoutwardflux( &
252 netoutwardflux(:,:), & ! (out)
253 rhoq(:,:,:,iq), vterm(:,:,:,iq), dzrhoq(:,:,:), ndcoefeuler(:,:,:), & ! (in)
254 lcmesh%J(:,:), lcmesh%Fscale(:,:), & ! (in)
255 nz(:,:,:), vmapm(:,:), vmapp(:,:), lcmesh%VMapM(:,:), intweight(:,:), & ! (in)
256 lcmesh, elem ) ! (in)
257
258 !$omp parallel do collapse(2) private(ke, Q)
259 do ke2d = 1, lcmesh%Ne2D
260 do ke_z = 1, lcmesh%NeZ
261 ke = ke2d + (ke_z-1)*lcmesh%Ne2D
262
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 ) ) )
265 end do ! end loop for ke_z
266 end do ! end loop for ke2D
267
268 call atm_phy_mp_dgm_precipitation_get_delflux( &
269 del_flux(:,:,:,:), & ! (out)
270 dens0(:,:,:), rhoq(:,:,:,iq), temp(:,:,:), vterm(:,:,:,iq), & ! (in)
271 dzrhoq(:,:,:), dzrhoe(:,:,:), ndcoefeuler(:,:,:), & ! (in)
272 fct_coef(:,:,:), & ! (in)
273 cv(iq), lcmesh%J(:,:), lcmesh%Fscale(:,:), nz(:,:,:), & ! (in)
274 vmapm(:,:), vmapp(:,:), lcmesh%vmapM(:,:), intweight(:,:), & ! (in)
275 lcmesh, elem ) ! (in)
276
277 !$omp parallel do collapse(2) private( &
278 !$omp ke2D, ke_z, ke, &
279 !$omp qflx, eflx, dDENS, RHOQ_tmp, RHOQ0, RHOQ1, &
280 !$omp RHOQ_save, &
281 !$omp Fz, LiftDelFlx )
282 do ke2d = 1, lcmesh%Ne2D
283 do ke_z = 1, lcmesh%NeZ
284 ke = ke2d + (ke_z-1)*lcmesh%Ne2D
285
286 !--- update falling tracer
287
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)
291
292 call sparsemat_matmul( dz, qflx(:), fz )
293 call sparsemat_matmul( lift, lcmesh%Fscale(:,ke) * del_flux(:,ke_z,ke2d,1), liftdelflx )
294
295 ddens(:) = - dt * ( &
296 lcmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) )
297 rhoq_tmp(:) = max( 0.0_rp, rhoq(:,ke_z,ke2d,iq) + ddens(:) )
298! RHOQ_tmp(:) = RHOQ(:,ke_z,ke2D,iq) + dDENS(:)
299
300 !
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(:) )
303
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(:)
307
308
309 ! QTRC(iq; iq>QLA+QLI) is not mass tracer, such as number density
310 if ( iq > qla + qia ) cycle
311
312 flx_hydro(:,ke_z,ke2d) = flx_hydro(:,ke_z,ke2d) &
313 + qflx(:) * rnstep
314 if ( ke_z == 1 ) then
315 if ( iq > qla ) then ! ice water
316 sflx_snow(:,ke2d) = sflx_snow(:,ke2d) &
317 + qflx(elem%Hslice(:,1)) * rnstep
318 else ! liquid water
319 sflx_rain(:,ke2d) = sflx_rain(:,ke2d) &
320 + qflx(elem%Hslice(:,1)) * rnstep
321 end if
322 end if
323
324 !--- update density
325
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(:)
329
330 !--- update internal energy
331
332 eflx(:) = vterm(:,ke_z,ke2d,iq) * rhoq_save(:) * temp(:,ke_z,ke2d) * cv(iq) &
333 - ndcoefeuler(:,ke_z,ke2d) * dzrhoe(:,ke_z,ke2d)
334
335 call sparsemat_matmul( dz, eflx(:), fz )
336 call sparsemat_matmul( lift, lcmesh%Fscale(:,ke) * del_flux(:,ke_z,ke2d,2), liftdelflx )
337
338 rhoe(:,ke_z,ke2d) = rhoe(:,ke_z,ke2d) - dt * ( &
339 + lcmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) &
340 + qflx(:) * grav )
341
342 if ( ke_z == 1 ) then
343 esflx(:,ke2d) = esflx(:,ke2d) &
344 + eflx(elem%Hslice(:,1)) * rnstep
345 end if
346
347 end do ! end loop for ke_z
348 end do ! end loop for ke2D
349 end do ! end loop for iq
350
351 !$omp parallel do collapse(2)
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)
356 end do
357 end do
358
359 return
361
362!OCL SERIAL
364 MOMU_t, MOMV_t, MOMZ_t, & ! (out)
365 dens, momu, momv, momz, mflx, & ! (in)
366 dz, lift, nz, vmapm, vmapp, & ! (in)
367 lcmesh, elem ) ! (in)
368 implicit none
369
370 class(localmesh3d), intent(in) :: lcmesh
371 class(elementbase3d), intent(in) :: elem
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)
380 type(sparsemat), intent(in) :: dz
381 type(sparsemat), intent(in) :: lift
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)
385
386 integer :: ke2d
387 integer :: ke_z
388 integer :: ke
389
390 real(rp) :: fz(elem%np), liftdelflx(elem%np)
391 real(rp) :: del_flux(elem%nfptot,lcmesh%nez,lcmesh%ne2d,3)
392
393 real(rp) :: rdens(elem%np)
394 !-------------------------------------------------------
395
396 call atm_phy_mp_dgm_precipitation_momentum_get_delflux( &
397 del_flux(:,:,:,:), & ! (out)
398 dens(:,:,:), momu(:,:,:), momv(:,:,:), momz(:,:,:), & ! (in)
399 mflx(:,:,:), & ! (in)
400 nz(:,:,:), vmapm(:,:), vmapp(:,:), & ! (in)
401 lcmesh, elem ) ! (in)
402
403 !$omp parallel do collapse(2) private( &
404 !$omp ke2D, ke_z, ke, RDENS, Fz, LiftDelFlx )
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)
409
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(:) )
413
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(:) )
417
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(:) )
421 end do
422 end do
423
424 return
426
427
428!OCL SERIAL
430 QTRC, DDENS, PRES, &
431 CVtot, CPtot, Rtot, &
432 DENS_hyd, PRES_hyd, &
433 dt, lmesh, elem, QA, QLA, QIA, &
434 DRHOT )
435
436 use scale_const, only: &
437 cvdry => const_cvdry, &
438 cpdry => const_cpdry, &
439 rdry => const_rdry
440
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
446 implicit none
447
448 class(localmesh3d), intent(in) :: lmesh
449 class(elementbase3d), intent(in) :: elem
450 integer, intent(in) :: qa
451 type(localmeshfieldbaselist), intent(inout) :: qtrc(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)
462
463 integer :: ke
464 integer :: iq
465
466 real(rp) :: int_w(elem%np)
467
468 real(rp) :: dens(elem%np)
469 real(rp) :: ddens0(elem%np)
470
471 real(rp) :: trcmass0(elem%np), trcmass1(elem%np,qa)
472 real(rp) :: mass0_elem, mass1_elem
473 real(rp) :: inten0_elem, inten_elem
474
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)
479
480#ifdef SINGLE
481 real(rp), parameter :: trc_eps = 1e-32_rp
482#else
483 real(rp), parameter :: trc_eps = 1e-128_rp
484#endif
485 !------------------------------------------------
486
487 !$omp parallel do private( &
488 !$omp ke, iq, DENS, DDENS0, InternalEn, InternalEn0, TEMP, QTRC_tmp, &
489 !$omp TRCMASS0, TRCMASS1, MASS0_elem, MASS1_elem, IntEn0_elem, IntEn_elem, &
490 !$omp Qdry, CVtot_old, CPtot_old, Rtot_old, &
491 !$omp int_w, RHOT_hyd )
492 do ke = lmesh%NeS, lmesh%NeE
493
494 do iq = 1, qa
495 qtrc_tmp(:,iq) = qtrc(iq)%ptr%val(:,ke)
496 end do
497 call atmos_thermodyn_specific_heat( &
498 elem%Np, 1, elem%Np, qa, & ! (in)
499 qtrc_tmp, tracer_mass(:), tracer_r(:), tracer_cv(:), tracer_cp(:), & ! (in)
500 qdry, rtot_old, cvtot_old, cptot_old ) ! (out)
501
502 dens(:) = dens_hyd(:,ke) + ddens(:,ke)
503 ddens0(:) = ddens(:,ke)
504
505 ! RHOT_hyd(:) = PRES00 / Rdry * ( PRES_hyd(:,ke) / PRES00 )**( CVdry / CPdry )
506 ! ( Internal energy ) = Cvtot * RHO * T = CVtot * RHO * ( PT * EXNER )
507 ! InternalEn0(:) = CVtot_old(:) * ( RHOT_hyd(:) + DRHOT(:,ke) ) &
508 ! * ( Rtot_old(:) * ( RHOT_hyd(:) + DRHOT(:,ke) ) / PRES00 )**( Rtot_old(:) / CVtot_old(:) )
509 !InternalEn0(:) = CVtot_old(:) * PRES(:,ke) / Rtot_old(:)
510 internalen0(:) = cvtot(:,ke) * pres(:,ke) / rtot(:,ke)
511 !TEMP(:) = InternalEn0(:) / ( DENS(:) * CVtot_old(:) )
512 temp(:) = internalen0(:) / ( dens(:) * cvtot(:,ke) )
513 internalen(:) = internalen0(:)
514
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(:) )
519
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)
523! TRCMASS1(:,iq) = MASS0_elem / MASS1_elem * TRCMASS1(:,iq)
524
525 ddens(:,ke) = ddens(:,ke) + ( trcmass1(:,iq) - trcmass0(:) )
526 internalen(:) = internalen(:) + ( trcmass1(:,iq) - trcmass0(:) ) * tracer_cv(iq) * temp(:)
527 end do
528
529 !--
530
531 dens(:) = dens_hyd(:,ke) + ddens(:,ke)
532 do iq = 1, qa
533 qtrc_tmp(:,iq) = trcmass1(:,iq) / dens(:)
534 qtrc(iq)%ptr%val(:,ke) = qtrc_tmp(:,iq)
535 end do
536
537 inten0_elem = sum( int_w(:) * internalen0(:) )
538 inten_elem = sum( int_w(:) * internalen(:) )
539 internalen(:) = inten0_elem / inten_elem * internalen(:)
540
541 call atmos_thermodyn_specific_heat( &
542 elem%Np, 1, elem%Np, qa, & ! (in)
543 qtrc_tmp, tracer_mass(:), tracer_r(:), tracer_cv(:), tracer_cp(:), & ! (in)
544 qdry, rtot(:,ke), cvtot(:,ke), cptot(:,ke) ) ! (out)
545
546 internalen(:) = internalen(:) - ( ddens(:,ke) - ddens0(:) ) * grav * lmesh%zlev(:,ke)
547 pres(:,ke) = internalen(:) * rtot(:,ke) / cvtot(:,ke)
548
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 )
552 end if
553 end do
554
555 return
557
558!- private --------------------------------
559
560!OCL SERIAL
561 subroutine atm_phy_mp_dgm_netoutwardflux( &
562 net_outward_flux, &
563 RHOQ_, vterm_, &
564 DzRHOQ_, NDcoefEuler_, &
565 J, Fscale, &
566 nz, vmapM, vmapP, vmapM3D, IntWeight, &
567 lmesh, elem )
568 implicit none
569
570 class(localmesh3d), intent(in) :: lmesh
571 class(elementbase3d), intent(in) :: elem
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)
584
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)
590
591 integer :: ke
592 integer :: ke_z, ke2d
593 integer :: ip(elem%nfptot), im(elem%nfptot)
594 integer :: im3d(elem%nfptot)
595 !------------------------------------------------------------------------
596
597 !$omp parallel do collapse(2) private( &
598 !$omp ke2D, ke_z, ke, iM3D, iM, iP, &
599 !$omp RHOQ_M, RHOQ_P, velM, velP, alpha, &
600 !$omp numflux, outward_flux_tmp )
601 do ke2d=1, lmesh%Ne2D
602 do ke_z=1, lmesh%NeZ
603 ke = ke2d + (ke_z-1)*lmesh%Ne2D
604
605 im3d(:) = vmapm3d(:,ke)
606 im(:) = vmapm(:,ke_z); ip(:) = vmapp(:,ke_z)
607
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(:)) )
613
614 where (nz(:,ke_z,ke2d) > 1.0e-10 .and. ip(:) == im(:) )
615 velp(:) = - velm(:)
616 end where
617
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(:) ) )
621
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(:) ) )
624 end do
625 end do
626
627 return
628 end subroutine atm_phy_mp_dgm_netoutwardflux
629
630!OCL SERIAL
631 subroutine atm_phy_mp_dgm_precipitation_get_delflux_dq( &
632 del_flux, &
633 DENS_, RHOQ_,TEMP_, CV, nz, vmapM, vmapP, &
634 lmesh, elem )
635
636 implicit none
637
638 class(localmesh3d), intent(in) :: lmesh
639 class(elementbase3d), intent(in) :: elem
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)
648
649 integer :: ke
650 integer :: ke_z, ke2d
651 real(rp) :: rhoq_p(elem%nfptot), rhoq_m(elem%nfptot)
652 integer :: im(elem%nfptot), ip(elem%nfptot)
653 !-----------------------------------------
654
655 !$omp parallel do collapse(2) private( &
656 !$omp ke2D, ke_z, ke, iM, iP, RHOQ_M, RHOQ_P )
657 do ke2d=1, lmesh%Ne2D
658 do ke_z=1, lmesh%NeZ
659 ke = ke2d + (ke_z-1)*lmesh%Ne2D
660 im(:) = vmapm(:,ke_z); ip(:) = vmapp(:,ke_z)
661
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)
666 end do
667 end do
668
669 return
670 end subroutine atm_phy_mp_dgm_precipitation_get_delflux_dq
671
672
673!OCL SERIAL
674 subroutine atm_phy_mp_dgm_precipitation_get_delflux( &
675 del_flux, &
676 DENS_, RHOQ_, TEMP_, vterm_, &
677 DzRHOQ_, DzRHOE_, NDcoefEuler_, &
678 fct_coef_, CV, &
679 J, Fscale, nz, vmapM, vmapP, vmapM3D, IntWeight, &
680 lmesh, elem )
681
682 implicit none
683
684 class(localmesh3d), intent(in) :: lmesh
685 class(elementbase3d), intent(in) :: elem
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)
703
704 integer :: ke
705 integer :: ke_z, ke2d
706 integer :: f, p, fp
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)
712
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)
719 real(rp) :: r
720 !-----------------------------------------
721
722 !$omp parallel do collapse(2) private( &
723 !$omp ke2D, ke_z, ke, iM3D, iM, iP, &
724 !$omp R_M, R_P, RHOQ_M, RHOQ_P, TEMP_M, TEMP_P, &
725 !$omp velM, velP, alpha, numflux, numflux_ei, &
726 !$omp NDcoef_M, NDcoef_P, &
727 !$omp outward_flux_tmp, &
728 !$omp f, p, fp, R )
729 do ke2d=1, lmesh%Ne2D
730 do ke_z=1, lmesh%NeZ
731 ke = ke2d + (ke_z-1)*lmesh%Ne2D
732
733 im3d(:) = vmapm3d(:,ke)
734 im(:) = vmapm(:,ke_z); ip(:) = vmapp(:,ke_z)
735
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)
742
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(:)) )
746
747 where (nz(:,ke_z,ke2d) > 1.0e-10 .and. ip(:) == im(:) )
748 velp(:) = - velm(:)
749 end where
750
751 ndcoef_m(:) = ndcoefeuler_(im(:),ke2d)
752 ndcoef_p(:) = ndcoefeuler_(ip(:),ke2d)
753
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(:) ) )
757
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(:) ) )
761
762
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
767 do p=1, elem%Nfp_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)
776 end do
777 end do
778
779 end do
780 end do
781
782 return
783 end subroutine atm_phy_mp_dgm_precipitation_get_delflux
784
785!OCL SERIAL
786 subroutine atm_phy_mp_dgm_precipitation_momentum_get_delflux( &
787 del_flux, & ! (out)
788 dens_, momu_, momv_, momz_, mflx_, & ! (in)
789 nz, vmapm, vmapp, lmesh, elem ) ! (in)
790
791 implicit none
792
793 class(localmesh3d), intent(in) :: lmesh
794 class(elementbase3d), intent(in) :: 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)
804
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)
810 !-----------------------------------------
811
812 !$omp parallel do collapse(2) private( &
813 !$omp ke2D, ke_z, iM, iP, alpha, &
814 !$omp densM, densP, VelM, VelP )
815 do ke2d=1, lmesh%Ne2D
816 do ke_z=1, lmesh%NeZ
817 im(:) = vmapm(:,ke_z); ip(:) = vmapp(:,ke_z)
818
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(:)) )
824
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) ) )
828
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) ) )
832
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) ) )
836 end do
837 end do
838
839 return
840 end subroutine atm_phy_mp_dgm_precipitation_momentum_get_delflux
841
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 / 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.