FE-Project
Loading...
Searching...
No Matches
Functions/Subroutines
scale_atm_phy_mp_dgm_common Module Reference

module FElib / Atmosphere / Physics cloud microphysics / common More...

Functions/Subroutines

subroutine, public atm_phy_mp_dgm_common_gen_intweight (intweight, lcmesh)
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)
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)

Detailed Description

module FElib / Atmosphere / Physics cloud microphysics / common

Description
cloud microphysics process common subroutines

To preserve nonnegativity in precipitation process, a limiter proposed by Light and Durran (2016, MWR) is used

Reference
  • Light and Durran 2016: Preserving Nonnegativity in Discontinuous Galerkin Approximations to Scalar Transport via Truncation and Mass Aware Rescaling (TMAR). Monthly Weather Review, 144(12), 4771–4786.
Author
Yuta Kawai, Team SCALE

Function/Subroutine Documentation

◆ atm_phy_mp_dgm_common_gen_intweight()

subroutine, public scale_atm_phy_mp_dgm_common::atm_phy_mp_dgm_common_gen_intweight ( real(rp), dimension(lcmesh%refelem3d%nfaces,lcmesh%refelem3d%nfptot), intent(out) intweight,
class(localmesh3d), target lcmesh )

Definition at line 76 of file scale_atm_phy_mp_dgm_common.F90.

79 use scale_polynominal, only: polynominal_gengausslobattoptintweight
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(:) = polynominal_gengausslobattoptintweight(elem%PolyOrder_h)
104 intweight_lgl1dpts_v(:) = polynominal_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
scale_polynominal
Module common / Polynominal.
Definition scale_polynominal.F90:11
scale_polynominal::polynominal_gengausslobattoptintweight
real(rp) function, dimension(nord+1), public polynominal_gengausslobattoptintweight(nord)
A function to calculate the Gauss-Lobbato weights.
Definition scale_polynominal.F90:272

References scale_polynominal::polynominal_gengausslobattoptintweight().

◆ atm_phy_mp_dgm_common_precipitation()

subroutine, public scale_atm_phy_mp_dgm_common::atm_phy_mp_dgm_common_precipitation ( real(rp), dimension (elem%np,lcmesh%nez,lcmesh%ne2d), intent(inout) dens,
real(rp), dimension (elem%np,lcmesh%nez,lcmesh%ne2d,qha), intent(inout) rhoq,
real(rp), dimension(elem%np,lcmesh%nez,lcmesh%ne2d), intent(inout) cptot,
real(rp), dimension(elem%np,lcmesh%nez,lcmesh%ne2d), intent(inout) cvtot,
real(rp), dimension (elem%np,lcmesh%nez,lcmesh%ne2d), intent(inout) rhoe,
real(rp), dimension(elem%np,lcmesh%nez,lcmesh%ne2d), intent(inout) flx_hydro,
real(rp), dimension(elem%nfp_v,lcmesh%ne2da), intent(inout) sflx_rain,
real(rp), dimension(elem%nfp_v,lcmesh%ne2da), intent(inout) sflx_snow,
real(rp), dimension (elem%nfp_v,lcmesh%ne2da), intent(inout) esflx,
real(rp), dimension (elem%np,lcmesh%nez,lcmesh%ne2d), intent(in) temp,
real(rp), dimension(elem%np,lcmesh%nez,lcmesh%ne2d,qha), intent(in) vterm,
real(rp), intent(in) dt,
real(rp), intent(in) rnstep,
type(sparsemat), intent(in) dz,
type(sparsemat), intent(in) lift,
real(rp), dimension(elem%nfptot,lcmesh%nez,lcmesh%ne2d), intent(in) nz,
integer, dimension(elem%nfptot,lcmesh%nez), intent(in) vmapm,
integer, dimension(elem%nfptot,lcmesh%nez), intent(in) vmapp,
real(rp), dimension(elem%nfaces,elem%nfptot), intent(in) intweight,
integer, intent(in) qha,
integer, intent(in) qla,
integer, intent(in) qia,
class(localmesh3d), intent(in) lcmesh,
class(elementbase3d), intent(in) elem )
Parameters
[in]qhahydrometeor (water + ice)

Definition at line 136 of file scale_atm_phy_mp_dgm_common.F90.

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

◆ atm_phy_mp_dgm_common_precipitation_momentum()

subroutine, public scale_atm_phy_mp_dgm_common::atm_phy_mp_dgm_common_precipitation_momentum ( real(rp), dimension(elem%np,lcmesh%nea), intent(out) momu_t,
real(rp), dimension(elem%np,lcmesh%nea), intent(out) momv_t,
real(rp), dimension(elem%np,lcmesh%nea), intent(out) momz_t,
real(rp), dimension(elem%np,lcmesh%nez,lcmesh%ne2d), intent(in) dens,
real(rp), dimension(elem%np,lcmesh%nez,lcmesh%ne2d), intent(in) momu,
real(rp), dimension(elem%np,lcmesh%nez,lcmesh%ne2d), intent(in) momv,
real(rp), dimension(elem%np,lcmesh%nez,lcmesh%ne2d), intent(in) momz,
real(rp), dimension(elem%np,lcmesh%nez,lcmesh%ne2d), intent(in) mflx,
type(sparsemat), intent(in) dz,
type(sparsemat), intent(in) lift,
real(rp), dimension(elem%nfptot,lcmesh%nez,lcmesh%ne2d), intent(in) nz,
integer, dimension(elem%nfptot,lcmesh%nez), intent(in) vmapm,
integer, dimension(elem%nfptot,lcmesh%nez), intent(in) vmapp,
class(localmesh3d), intent(in) lcmesh,
class(elementbase3d), intent(in) elem )

Definition at line 363 of file scale_atm_phy_mp_dgm_common.F90.

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

◆ atm_phy_mp_dgm_common_negative_fixer()

subroutine, public scale_atm_phy_mp_dgm_common::atm_phy_mp_dgm_common_negative_fixer ( type(localmeshfieldbaselist), dimension(qa), intent(inout) qtrc,
real(rp), dimension(elem%np,lmesh%nea), intent(inout) ddens,
real(rp), dimension(elem%np,lmesh%nea), intent(inout) pres,
real(rp), dimension(elem%np,lmesh%nea), intent(inout) cvtot,
real(rp), dimension(elem%np,lmesh%nea), intent(inout) cptot,
real(rp), dimension(elem%np,lmesh%nea), intent(inout) rtot,
real(rp), dimension(elem%np,lmesh%nea), intent(in) dens_hyd,
real(rp), dimension(elem%np,lmesh%nea), intent(in) pres_hyd,
real(rp), intent(in) dt,
class(localmesh3d), intent(in) lmesh,
class(elementbase3d), intent(in) elem,
integer, intent(in) qa,
integer, intent(in) qla,
integer, intent(in) qia,
real(rp), dimension(elem%np,lmesh%nea), intent(inout), optional drhot )

Definition at line 429 of file scale_atm_phy_mp_dgm_common.F90.

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
445 use scale_localmeshfield_base, only: localmeshfieldbaselist
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
scale_localmeshfield_base
module FElib / Data / base
Definition scale_localmeshfield_base.F90:10
scale_localmeshfield_base::localmeshfieldbaselist
Definition scale_localmeshfield_base.F90:48
Generated by doxygen 1.14.0