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

module FElib / Atmosphere / Physics turbulence More...

Functions/Subroutines

subroutine, public atm_phy_tb_dgm_smg_init (mesh)
 
subroutine, public atm_phy_tb_dgm_smg_final ()
 
subroutine, public atm_phy_tb_dgm_smg_cal_grad (t11, t12, t13, t21, t22, t23, t31, t32, t33, df1, df2, df3, tke, nu, kh, ddens_, momx_, momy_, momz_, drhot_, dens_hyd, pres_hyd, pres, pt, dx, dy, dz, sx, sy, sz, lift, lmesh, elem, lmesh2d, elem2d, is_bound)
 Calculate parameterized stress tensor and eddy heat flux with turbulent model.
 

Detailed Description

module FElib / Atmosphere / Physics turbulence

Description
Sub-grid scale turbulnce process Smagorinsky-type
Author
Yuta Kawai, Team SCALE
Reference
  • Brown et al., 1994: Large-eddy simulaition of stable atmospheric boundary layers with a revised stochastic subgrid model. Roy. Meteor. Soc., 120, 1485-1512
  • Scotti et al., 1993: Generalized Smagorinsky model for anisotropic grids. Phys. Fluids A, 5, 2306-2308
  • Nishizawa et al., 2015: Influence of grid aspect ratio on planetary boundary layer turbulence in large-eddy simulations Geosci. Model Dev., 8, 3393–3419
NAMELIST
  • PARAM_ATMOS_PHY_TB_DGM_SMG
    nametypedefault valuecomment
    CS real(RP) 0.13_RP Smagorinsky constant (Scotti et al. 1993)
    NU_MAX real(RP) 10000.0_RP
    FILTER_FAC real(RP) 2.0_RP
    CONSISTENT_TKE logical .true.

History Output
No history output

Function/Subroutine Documentation

◆ atm_phy_tb_dgm_smg_init()

subroutine, public scale_atm_phy_tb_dgm_smg::atm_phy_tb_dgm_smg_init ( class(meshbase3d), intent(in) mesh)

Definition at line 100 of file scale_atm_phy_tb_dgm_smg.F90.

101 implicit none
102 class(MeshBase3D), intent(in) :: mesh
103
104 logical :: consistent_tke = .true.
105
106 namelist / param_atmos_phy_tb_dgm_smg / &
107 cs, &
108 nu_max, &
109 filter_fac, &
110 consistent_tke
111
112 integer :: ierr
113 !--------------------------------------------------------------------
114
115 log_newline
116 log_info("ATMOS_PHY_TB_dgm_smg_setup",*) 'Setup'
117 log_info("ATMOS_PHY_TB_dgm_smg_setup",*) 'Smagorinsky-type Eddy Viscocity Model'
118
119 !--- read namelist
120 rewind(io_fid_conf)
121 read(io_fid_conf,nml=param_atmos_phy_tb_dgm_smg,iostat=ierr)
122 if( ierr < 0 ) then !--- missing
123 log_info("ATMOS_PHY_TB_dgm_smg_setup",*) 'Not found namelist. Default used.'
124 elseif( ierr > 0 ) then !--- fatal error
125 log_error("ATMOS_PHY_TB_dgm_smg_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_TB_DGM_SMG. Check!'
126 call prc_abort
127 endif
128 log_nml(param_atmos_phy_tb_dgm_smg)
129
130 rprn = 1.0_rp / prn
131 rric = 1.0_rp / ric
132 onemprnovric = ( 1.0_rp - prn ) * rric
133
134 if ( consistent_tke ) then
135 tke_fac = 1.0_rp
136 else
137 tke_fac = 0.0_rp
138 end if
139
140 return

◆ atm_phy_tb_dgm_smg_final()

subroutine, public scale_atm_phy_tb_dgm_smg::atm_phy_tb_dgm_smg_final

Definition at line 144 of file scale_atm_phy_tb_dgm_smg.F90.

145 implicit none
146 !--------------------------------------------------------------------
147
148 return

◆ atm_phy_tb_dgm_smg_cal_grad()

subroutine, public scale_atm_phy_tb_dgm_smg::atm_phy_tb_dgm_smg_cal_grad ( real(rp), dimension(elem%np,lmesh%nea), intent(out) t11,
real(rp), dimension(elem%np,lmesh%nea), intent(out) t12,
real(rp), dimension(elem%np,lmesh%nea), intent(out) t13,
real(rp), dimension(elem%np,lmesh%nea), intent(out) t21,
real(rp), dimension(elem%np,lmesh%nea), intent(out) t22,
real(rp), dimension(elem%np,lmesh%nea), intent(out) t23,
real(rp), dimension(elem%np,lmesh%nea), intent(out) t31,
real(rp), dimension(elem%np,lmesh%nea), intent(out) t32,
real(rp), dimension(elem%np,lmesh%nea), intent(out) t33,
real(rp), dimension(elem%np,lmesh%nea), intent(out) df1,
real(rp), dimension(elem%np,lmesh%nea), intent(out) df2,
real(rp), dimension(elem%np,lmesh%nea), intent(out) df3,
real(rp), dimension(elem%np,lmesh%nea), intent(out) tke,
real(rp), dimension(elem%np,lmesh%nea), intent(out) nu,
real(rp), dimension(elem%np,lmesh%nea), intent(out) kh,
real(rp), dimension(elem%np,lmesh%nea), intent(in) ddens_,
real(rp), dimension (elem%np,lmesh%nea), intent(in) momx_,
real(rp), dimension (elem%np,lmesh%nea), intent(in) momy_,
real(rp), dimension (elem%np,lmesh%nea), intent(in) momz_,
real(rp), dimension(elem%np,lmesh%nea), intent(in) drhot_,
real(rp), dimension(elem%np,lmesh%nea), intent(in) dens_hyd,
real(rp), dimension(elem%np,lmesh%nea), intent(in) pres_hyd,
real(rp), dimension(elem%np,lmesh%nea), intent(in) pres,
real(rp), dimension(elem%np,lmesh%nea), intent(in) pt,
type(sparsemat), intent(in) dx,
type(sparsemat), intent(in) dy,
type(sparsemat), intent(in) dz,
type(sparsemat), intent(in) sx,
type(sparsemat), intent(in) sy,
type(sparsemat), intent(in) sz,
type(sparsemat), intent(in) lift,
class(localmesh3d), intent(in) lmesh,
class(elementbase3d), intent(in) elem,
class(localmesh2d), intent(in) lmesh2d,
class(elementbase2d), intent(in) elem2d,
logical, dimension(elem%nfptot,lmesh%ne), intent(in) is_bound )

Calculate parameterized stress tensor and eddy heat flux with turbulent model.

Parameters
[out]t11(1,1) component of stress tensor
[out]t12(1,2) component of stress tensor
[out]t13(1,3) component of stress tensor
[out]t21(2,1) component of stress tensor
[out]t22(2,2) component of stress tensor
[out]t23(2,3) component of stress tensor
[out]t31(3,1) component of stress tensor
[out]t32(3,2) component of stress tensor
[out]t33(3,3) component of stress tensor
[out]df1Diffusive heat flux in x1 direction / density
[out]df2Diffusive heat flux in x2 direction / density
[out]df3Diffusive heat flux in x3 direction / density
[out]tkeParameterized turbulent kinetic energy
[out]nuEddy viscosity
[out]khEddy diffusivity
[in]ddens_Density perturbation
[in]momx_Momentum in x1 direction
[in]momy_Momentum in x2 direction
[in]momz_Momentum in x3 direction
[in]drhot_Density x potential temperature perturbation
[in]dens_hydReference pressure in hydrostatic balance
[in]pres_hydReference density in hydrostatic balance
[in]presPressure
[in]ptPotential temperature
[in]dzDifferential matrix managed by sparse matrix type
[in]szStiffness matrix managed by sparse matrix type
[in]liftLifting matrix managed by sparse matrix type
[in]is_boundFlag whether nodes are located at domain boundaries

Definition at line 154 of file scale_atm_phy_tb_dgm_smg.F90.

162
163 use scale_atm_phy_tb_dgm_common, only: &
164 atm_phy_tb_dgm_common_calc_lambda
165 implicit none
166
167 class(LocalMesh3D), intent(in) :: lmesh
168 class(ElementBase3D), intent(in) :: elem
169 class(LocalMesh2D), intent(in) :: lmesh2D
170 class(ElementBase2D), intent(in) :: elem2D
171 real(RP), intent(out) :: T11(elem%Np,lmesh%NeA)
172 real(RP), intent(out) :: T12(elem%Np,lmesh%NeA)
173 real(RP), intent(out) :: T13(elem%Np,lmesh%NeA)
174 real(RP), intent(out) :: T21(elem%Np,lmesh%NeA)
175 real(RP), intent(out) :: T22(elem%Np,lmesh%NeA)
176 real(RP), intent(out) :: T23(elem%Np,lmesh%NeA)
177 real(RP), intent(out) :: T31(elem%Np,lmesh%NeA)
178 real(RP), intent(out) :: T32(elem%Np,lmesh%NeA)
179 real(RP), intent(out) :: T33(elem%Np,lmesh%NeA)
180 real(RP), intent(out) :: DF1(elem%Np,lmesh%NeA)
181 real(RP), intent(out) :: DF2(elem%Np,lmesh%NeA)
182 real(RP), intent(out) :: DF3(elem%Np,lmesh%NeA)
183 real(RP), intent(out) :: TKE(elem%Np,lmesh%NeA)
184 real(RP), intent(out) :: Nu(elem%Np,lmesh%NeA)
185 real(RP), intent(out) :: Kh(elem%Np,lmesh%NeA)
186 real(RP), intent(in) :: DDENS_(elem%Np,lmesh%NeA)
187 real(RP), intent(in) :: MOMX_ (elem%Np,lmesh%NeA)
188 real(RP), intent(in) :: MOMY_ (elem%Np,lmesh%NeA)
189 real(RP), intent(in) :: MOMZ_ (elem%Np,lmesh%NeA)
190 real(RP), intent(in) :: DRHOT_(elem%Np,lmesh%NeA)
191 real(RP), intent(in) :: DENS_hyd(elem%Np,lmesh%NeA)
192 real(RP), intent(in) :: PRES_hyd(elem%Np,lmesh%NeA)
193 real(RP), intent(in) :: PRES(elem%Np,lmesh%NeA)
194 real(RP), intent(in) :: PT(elem%Np,lmesh%NeA)
195 type(SparseMat), intent(in) :: Dx, Dy, Dz
196 type(SparseMat), intent(in) :: Sx, Sy, Sz
197 type(SparseMat), intent(in) :: Lift
198 logical, intent(in) :: is_bound(elem%NfpTot,lmesh%Ne)
199
200 real(RP) :: Fx(elem%Np), Fy(elem%Np), Fz(elem%Np), LiftDelFlx(elem%Np)
201 real(RP) :: DENS(elem%Np), RDENS(elem%Np), RHOT(elem%Np), Q(elem%Np)
202 real(RP) :: DdensDxi(elem%Np,3)
203 real(RP) :: DVelDxi(elem%Np,3,3)
204 real(RP) :: del_flux_rho (elem%NfpTot,lmesh%Ne,3)
205 real(RP) :: del_flux_mom (elem%NfpTot,lmesh%Ne,3,3)
206 real(RP) :: del_flux_rhot(elem%NfpTot,lmesh%Ne,3)
207
208 real(RP) :: S11(elem%Np), S12(elem%Np), S22(elem%Np), S23(elem%Np), S31(elem%Np), S33(elem%Np)
209 real(RP) :: SkkOvThree
210 real(RP) :: TKEMulTwoOvThree
211 real(RP) :: coef
212
213 real(RP) :: Ri ! local gradient Richardson number
214 real(RP) :: S2 ! (2SijSij)^1/2
215 real(RP) :: fm ! factor in eddy viscosity which represents the stability dependence of
216 ! the Brown et al (1994)'s subgrid model
217 real(RP) :: Pr ! Parandtl number (=Nu/Kh= fm/fh)
218
219 real(RP) :: lambda (elem%Np,lmesh%Ne) ! basic mixing length
220 real(RP) :: lambda_r(elem%Np) ! characteristic subgrid length scale
221 real(RP) :: E(elem%Np) ! subgrid kinetic energy
222 real(RP) :: C1(elem%Np) ! factor in the relation with energy disspation rate, lambda_r, and E
223
224 integer :: ke
225 integer :: p
226 !--------------------------------------------------------------------
227
228 call cal_del_flux_grad( del_flux_rho, del_flux_mom, del_flux_rhot, & ! (out)
229 ddens_, momx_, momy_, momz_, drhot_, dens_hyd, pres_hyd, pt, & ! (in)
230 lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), & ! (in)
231 lmesh%vmapM, lmesh%vmapP, & ! (in)
232 lmesh, elem, is_bound ) ! (in)
233
234 call atm_phy_tb_dgm_common_calc_lambda( lambda, & ! (out)
235 cs, filter_fac, lmesh, elem, lmesh2d, elem2d ) ! (in)
236
237 !$omp parallel do private( &
238 !$omp Fx, Fy, Fz, LiftDelFlx, &
239 !$omp DENS, RHOT, RDENS, Q, DdensDxi, DVelDxi, &
240 !$omp S11, S12, S22, S23, S31, S33, &
241 !$omp coef, SkkOvThree, TKEMulTwoOvThree, &
242 !$omp p, Ri, S2, fm, Pr, lambda_r, E, C1 )
243 do ke=lmesh%NeS, lmesh%NeE
244 !---
245 dens(:) = dens_hyd(:,ke) + ddens_(:,ke)
246 rdens(:) = 1.0_rp / dens(:)
247 rhot(:) = dens(:) * pt(:,ke)
248
249 ! gradient of density
250 call sparsemat_matmul( dx, dens, fx )
251 call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_rho(:,ke,1), liftdelflx )
252 ddensdxi(:,1) = lmesh%Escale(:,ke,1,1) * fx(:) + liftdelflx(:)
253
254 call sparsemat_matmul( dy, dens, fy )
255 call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_rho(:,ke,2), liftdelflx )
256 ddensdxi(:,2) = lmesh%Escale(:,ke,2,2) * fy(:) + liftdelflx(:)
257
258 call sparsemat_matmul( dz, dens, fz )
259 call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_rho(:,ke,3), liftdelflx )
260 ddensdxi(:,3) = lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:)
261
262 ! gradient of u
263 q(:) = momx_(:,ke) * rdens(:)
264
265 call sparsemat_matmul( dx, momx_(:,ke), fx )
266 call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,1,1), liftdelflx )
267 dveldxi(:,1,1) = ( lmesh%Escale(:,ke,1,1) * fx(:) + liftdelflx(:) - q(:) * ddensdxi(:,1) ) * rdens(:)
268
269 call sparsemat_matmul( dy, momx_(:,ke), fy )
270 call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,2,1), liftdelflx )
271 dveldxi(:,2,1) = ( lmesh%Escale(:,ke,2,2) * fy(:) + liftdelflx(:) - q(:) * ddensdxi(:,2) ) * rdens(:)
272
273 call sparsemat_matmul( dz, momx_(:,ke), fz )
274 call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,3,1), liftdelflx )
275 dveldxi(:,3,1) = ( lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) - q(:) * ddensdxi(:,3) ) * rdens(:)
276
277 ! gradient of v
278 q(:) = momy_(:,ke) * rdens(:)
279
280 call sparsemat_matmul( dx, momy_(:,ke), fx )
281 call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,1,2), liftdelflx )
282 dveldxi(:,1,2) = ( lmesh%Escale(:,ke,1,1) * fx(:) + liftdelflx(:) - q(:) * ddensdxi(:,1) ) * rdens(:)
283
284 call sparsemat_matmul( dy, momy_(:,ke), fy )
285 call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,2,2), liftdelflx )
286 dveldxi(:,2,2) = ( lmesh%Escale(:,ke,2,2) * fy(:) + liftdelflx(:) - q(:) * ddensdxi(:,2) ) * rdens(:)
287
288 call sparsemat_matmul( dz, momy_(:,ke), fz )
289 call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,3,2), liftdelflx )
290 dveldxi(:,3,2) = ( lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) - q(:) * ddensdxi(:,3) ) * rdens(:)
291
292 ! gradient of w
293 q(:) = momz_(:,ke) * rdens(:)
294
295 call sparsemat_matmul( dx, momz_(:,ke), fx )
296 call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,1,3), liftdelflx )
297 dveldxi(:,1,3) = ( lmesh%Escale(:,ke,1,1) * fx(:) + liftdelflx(:) - q(:) * ddensdxi(:,1) ) * rdens(:)
298
299 call sparsemat_matmul( dy, momz_(:,ke), fy )
300 call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,2,3), liftdelflx )
301 dveldxi(:,2,3) = ( lmesh%Escale(:,ke,2,2) * fy(:) + liftdelflx(:) - q(:) * ddensdxi(:,2) ) * rdens(:)
302
303 call sparsemat_matmul( dz, momz_(:,ke), fz )
304 call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,3,3), liftdelflx )
305 dveldxi(:,3,3) = ( lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) - q(:) * ddensdxi(:,3) ) * rdens(:)
306
307 ! gradient of pt
308 q(:) = rhot(:) * rdens(:)
309
310 call sparsemat_matmul( dx, rhot, fx )
311 call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_rhot(:,ke,1), liftdelflx )
312 df1(:,ke) = ( lmesh%Escale(:,ke,1,1) * fx(:) + liftdelflx(:) - q(:) * ddensdxi(:,1) ) * rdens(:)
313
314 call sparsemat_matmul( dy, rhot, fy )
315 call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_rhot(:,ke,2), liftdelflx )
316 df2(:,ke) = ( lmesh%Escale(:,ke,2,2) * fy(:) + liftdelflx(:) - q(:) * ddensdxi(:,2) ) * rdens(:)
317
318 call sparsemat_matmul( dz, rhot, fz )
319 call sparsemat_matmul( lift, lmesh%Fscale(:,ke) * del_flux_rhot(:,ke,3), liftdelflx )
320 df3(:,ke) = ( lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) - q(:) * ddensdxi(:,3) ) * rdens(:)
321
322
323 ! Calculate the component of strain velocity tensor
324 s11(:) = dveldxi(:,1,1)
325 s12(:) = 0.5_rp * ( dveldxi(:,1,2) + dveldxi(:,2,1) )
326 s22(:) = dveldxi(:,2,2)
327 s23(:) = 0.5_rp * ( dveldxi(:,2,3) + dveldxi(:,3,2) )
328 s31(:) = 0.5_rp * ( dveldxi(:,1,3) + dveldxi(:,3,1) )
329 s33(:) = dveldxi(:,3,3)
330
331 ! Caclulate eddy viscosity & eddy diffusivity
332
333 do p=1, elem%Np
334 s2 = 2.0_rp * ( s11(p)**2 + s22(p)**2 + s33(p)**2 ) &
335 + 4.0_rp * ( s31(p)**2 + s12(p)**2 + s23(p)**2 )
336
337 ri = grav / pt(p,ke) * df3(p,ke) / max( s2, eps )
338
339 ! The Stability functions fm and fh are given by the appendix A of Brown et al. (1994).
340 if (ri < 0.0_rp ) then ! unstable
341 fm = sqrt( 1.0_rp - fmc * ri )
342 nu(p,ke) = lambda(p,ke)**2 * sqrt( s2 ) * fm
343 pr = fm / sqrt( 1.0_rp - fhb * ri ) * prn
344 else if ( ri < ric ) then ! stable
345 fm = ( 1.0_rp - ri * rric )**4
346 nu(p,ke) = lambda(p,ke)**2 * sqrt( s2 ) * fm
347 pr = prn / ( 1.0_rp - onemprnovric * ri )
348 else ! strongly stable
349 fm = 0.0_rp
350 nu(p,ke) = 0.0_rp
351 kh(p,ke) = 0.0_rp
352 pr = 1.0_rp
353 end if
354
355 if ( ri < ric ) then
356 kh(p,ke) = max( min( nu(p,ke) / pr, nu_max ), eps )
357 nu(p,ke) = max( min( nu(p,ke), nu_max ), eps )
358 pr = nu(p,ke) / kh(p,ke)
359 lambda_r(p) = lambda(p,ke) * sqrt( fm / sqrt( 1.0_rp - ri/pr ) )
360 else
361 lambda_r(p) = 0.0_rp
362 end if
363 end do
364
365! Nu(:,ke) = 0.0_RP
366
367 ! if ( backscatter ) then
368 ! else
369 e(:) = nu(:,ke)**3 / ( lambda_r(:)**4 + eps )
370 c1(:) = c1o
371 ! end if
372
373 ! TKE
374 tke(:,ke) = ( e(:) * lambda_r(:) / c1(:) )**twooverthree
375
376 !---
377
378 do p=1, elem%Np
379 tkemultwoovthree = twooverthree * tke(p,ke) * tke_fac
380 skkovthree = ( s11(p) + s22(p) + s33(p) ) * oneoverthree
381 coef = 2.0_rp * nu(p,ke)
382
383 t11(p,ke) = dens(p) * ( coef * ( s11(p) - skkovthree ) - tkemultwoovthree )
384 t12(p,ke) = dens(p) * coef * s12(p)
385 t13(p,ke) = dens(p) * coef * s31(p)
386
387 t21(p,ke) = dens(p) * coef * s12(p)
388 t22(p,ke) = dens(p) * ( coef * ( s22(p) - skkovthree ) - tkemultwoovthree )
389 t23(p,ke) = dens(p) * coef * s23(p)
390
391 t31(p,ke) = dens(p) * coef * s31(p)
392 t32(p,ke) = dens(p) * coef * s23(p)
393 t33(p,ke) = dens(p) * ( coef * ( s33(p) - skkovthree ) - tkemultwoovthree )
394 end do
395
396 df1(:,ke) = kh(:,ke) * df1(:,ke)
397 df2(:,ke) = kh(:,ke) * df2(:,ke)
398 df3(:,ke) = kh(:,ke) * df3(:,ke)
399 end do
400
401 return
scale_atm_phy_tb_dgm_common
module FElib / Fluid dyn solver / Atmosphere / Physics turbulence / Common
Definition scale_atm_phy_tb_dgm_common.F90:11
scale_atm_phy_tb_dgm_common::atm_phy_tb_dgm_common_calc_lambda
subroutine, public atm_phy_tb_dgm_common_calc_lambda(lambda, cs, filter_fac, lmesh, elem, lmesh2d, elem2d)
Definition scale_atm_phy_tb_dgm_common.F90:226

References scale_atm_phy_tb_dgm_common::atm_phy_tb_dgm_common_calc_lambda().

Generated by doxygen 1.10.0