FE-Project
Loading...
Searching...
No Matches
scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform Module Reference

module Atmosphere / Dynamics HEVI More...

Functions/Subroutines

subroutine, public atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform_init (mesh)
subroutine, public atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform_final ()
subroutine, public atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform_cal_tend (dens_dt, momx_dt, momy_dt, momz_dt, rhot_dt, ddens_, momx_, momy_, momz_, drhot_, dpres_, dens_hyd, pres_hyd, pres_hyd_ref, therm_hyd, coriolis, rtot, cvtot, cptot, dphyddx, dphyddy, element3d_operation, dx, dy, dz, sx, sy, sz, lift, lmesh, elem, lmesh2d, elem2d)
subroutine, public atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform_cal_vi (dens_dt, momx_dt, momy_dt, momz_dt, rhot_dt, ddens_, momx_, momy_, momz_, drhot_, dens_hyd, pres_hyd, ddens0_, momx0_, momy0_, momz0_, drhot0_, rtot, cvtot, cptot, element3d_operation, dz, lift, impl_fac, dt, lmesh, elem, lmesh2d, elem2d)

Detailed Description

module Atmosphere / Dynamics HEVI

Description
HEVI DGM scheme for Atmospheric dynamical process. To improve the numerical instability due to the aliasing errors, the split form based on Gassner et al. (2016, JCP) is used for advection terms.
Author
Yuta Kawai, Team SCALE

Function/Subroutine Documentation

◆ atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform_init()

subroutine, public scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform::atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform_init ( class(meshbase3d), intent(in) mesh)

Definition at line 80 of file scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform.F90.

83 implicit none
84 class(MeshBase3D), intent(in) :: mesh
85
86 integer :: p1, p2, p3, p_
87 type(ElementBase3D), pointer :: elem
88 !--------------------------------------------
89
90 call atm_dyn_dgm_nonhydro3d_common_init( mesh )
92
93 !-
94 elem => mesh%refElem3D
95
96 allocate( dxt1d_(elem%Nnode_h1D,elem%Nnode_h1D) )
97 allocate( dyt1d_(elem%Nnode_h1D,elem%Nnode_h1D) )
98 allocate( dzt1d_(elem%Nnode_v,elem%Nnode_v) )
99
100 do p1=1, elem%Nnode_h1D
101 dxt1d_(:,p1) = elem%Dx1(p1,1:elem%Nnode_h1D)
102 end do
103
104 do p2=1, elem%Nnode_h1D
105 do p_=1, elem%Nnode_h1D
106 dyt1d_(p_,p2) = elem%Dx2(1+(p2-1)*elem%Nnode_h1D,1+(p_-1)*elem%Nnode_h1D)
107 end do
108 end do
109
110 do p3=1, elem%Nnode_v
111 dzt1d_(:,p3) = elem%Dx3(elem%Colmask(p3,1),elem%Colmask(:,1))
112 end do
113
114 return
module FElib / Fluid dyn solver / Atmosphere / Nonhydrostatic model / HEVI / Common

References scale_atm_dyn_dgm_nonhydro3d_common::atm_dyn_dgm_nonhydro3d_common_init(), and scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common::atm_dyn_dgm_nonhydro3d_rhot_hevi_common_init().

◆ atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform_final()

subroutine, public scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform::atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform_final

Definition at line 118 of file scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform.F90.

121 implicit none
122 !--------------------------------------------
123
125 call atm_dyn_dgm_nonhydro3d_common_final()
126
127 deallocate( dxt1d_, dyt1d_, dzt1d_ )
128
129 return

References scale_atm_dyn_dgm_nonhydro3d_common::atm_dyn_dgm_nonhydro3d_common_final(), and scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common::atm_dyn_dgm_nonhydro3d_rhot_hevi_common_final().

◆ atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform_cal_tend()

subroutine, public scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform::atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform_cal_tend ( real(rp), dimension(elem%np,lmesh%nea), intent(out) dens_dt,
real(rp), dimension(elem%np,lmesh%nea), intent(out) momx_dt,
real(rp), dimension(elem%np,lmesh%nea), intent(out) momy_dt,
real(rp), dimension(elem%np,lmesh%nea), intent(out) momz_dt,
real(rp), dimension(elem%np,lmesh%nea), intent(out) rhot_dt,
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) dpres_,
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_hyd_ref,
real(rp), dimension(elem%np,lmesh%nea), intent(in) therm_hyd,
real(rp), dimension(elem2d%np,lmesh2d%nea), intent(in) coriolis,
real(rp), dimension(elem%np,lmesh%nea), intent(in) rtot,
real(rp), dimension(elem%np,lmesh%nea), intent(in) cvtot,
real(rp), dimension(elem%np,lmesh%nea), intent(in) cptot,
real(rp), dimension(elem%np,lmesh%nea), intent(in) dphyddx,
real(rp), dimension(elem%np,lmesh%nea), intent(in) dphyddy,
class(elementoperationbase3d), intent(in) element3d_operation,
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 )

Definition at line 134 of file scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform.F90.

141
144
145 implicit none
146
147 class(LocalMesh3D), intent(in) :: lmesh
148 class(ElementBase3D), intent(in) :: elem
149 class(LocalMesh2D), intent(in) :: lmesh2D
150 class(ElementBase2D), intent(in) :: elem2D
151 class(ElementOperationBase3D), intent(in) :: element3D_operation
152 type(SparseMat), intent(in) :: Dx, Dy, Dz, Sx, Sy, Sz, Lift
153 real(RP), intent(out) :: DENS_dt(elem%Np,lmesh%NeA)
154 real(RP), intent(out) :: MOMX_dt(elem%Np,lmesh%NeA)
155 real(RP), intent(out) :: MOMY_dt(elem%Np,lmesh%NeA)
156 real(RP), intent(out) :: MOMZ_dt(elem%Np,lmesh%NeA)
157 real(RP), intent(out) :: RHOT_dt(elem%Np,lmesh%NeA)
158 real(RP), intent(in) :: DDENS_(elem%Np,lmesh%NeA)
159 real(RP), intent(in) :: MOMX_(elem%Np,lmesh%NeA)
160 real(RP), intent(in) :: MOMY_(elem%Np,lmesh%NeA)
161 real(RP), intent(in) :: MOMZ_(elem%Np,lmesh%NeA)
162 real(RP), intent(in) :: DRHOT_(elem%Np,lmesh%NeA)
163 real(RP), intent(in) :: DPRES_(elem%Np,lmesh%NeA)
164 real(RP), intent(in) :: DENS_hyd(elem%Np,lmesh%NeA)
165 real(RP), intent(in) :: PRES_hyd(elem%Np,lmesh%NeA)
166 real(RP), intent(in) :: PRES_hyd_ref(elem%Np,lmesh%NeA)
167 real(RP), intent(in) :: THERM_hyd(elem%Np,lmesh%NeA)
168 real(RP), intent(in) :: CORIOLIS(elem2D%Np,lmesh2D%NeA)
169 real(RP), intent(in) :: Rtot(elem%Np,lmesh%NeA)
170 real(RP), intent(in) :: CVtot(elem%Np,lmesh%NeA)
171 real(RP), intent(in) :: CPtot(elem%Np,lmesh%NeA)
172 real(RP), intent(in) :: DPhydDx(elem%Np,lmesh%NeA)
173 real(RP), intent(in) :: DPhydDy(elem%Np,lmesh%NeA)
174
175 real(RP) :: Fx(elem%Np), Fy(elem%Np), Fz(elem%Np), LiftDelFlx(elem%Np)
176 real(RP) :: Fx_sp(elem%Np), Fy_sp(elem%Np), Fz_sp(elem%Np)
177 real(RP) :: DPRES_hyd(elem%Np), GradPhyd_x(elem%Np), GradPhyd_y(elem%Np)
178 real(RP) :: del_flux(elem%NfpTot,lmesh%Ne,PRGVAR_NUM)
179 real(RP) :: del_flux_hyd(elem%NfpTot,lmesh%Ne,2)
180 real(RP) :: GsqrtDens_(elem%Np), rdens_(elem%Np), RHOT_hyd(elem%Np), RHOT_(elem%Np)
181 real(RP) :: u_(elem%Np), v_(elem%Np), w_(elem%Np), wt_(elem%Np), pot_(elem%Np)
182 real(RP) :: Cori(elem%Np)
183 real(RP) :: GsqrtV(elem%Np), RGsqrtV(elem%Np)
184
185 integer :: ke, ke2d
186
187 real(RP) :: gamm, rgamm
188 real(RP) :: rP0
189 real(RP) :: RovP0, P0ovR
190 !------------------------------------------------------------------------
191
192 call prof_rapstart( 'cal_dyn_tend_bndflux', 3)
193 call get_ebnd_flux( &
194 del_flux, del_flux_hyd, & ! (out)
195 ddens_, momx_, momy_, momz_, drhot_, dpres_, dens_hyd, pres_hyd, & ! (in)
196 rtot, cvtot, cptot, & ! (in)
197 lmesh%Gsqrt, lmesh%GI3(:,:,1), lmesh%GI3(:,:,2), & ! (in)
198 lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), & ! (in)
199 lmesh%vmapM, lmesh%vmapP, & ! (in)
200 lmesh, elem, lmesh2d, elem2d ) ! (in)
201 call prof_rapend( 'cal_dyn_tend_bndflux', 3)
202
203 !-----
204 call prof_rapstart( 'cal_dyn_tend_interior', 3)
205 gamm = cpdry / cvdry
206 rgamm = cvdry / cpdry
207 rp0 = 1.0_rp / pres00
208 rovp0 = rdry * rp0
209 p0ovr = pres00 / rdry
210
211 !$omp parallel do private( ke2d, Cori, &
212 !$omp RHOT_, GsqrtDens_, rdens_, u_, v_, w_, wt_, pot_, &
213 !$omp DPRES_hyd, GradPhyd_x, GradPhyd_y, &
214 !$omp GsqrtV, RGsqrtV, &
215 !$omp Fx, Fy, Fz, Fx_sp, Fy_sp, Fz_sp, LiftDelFlx )
216 do ke = lmesh%NeS, lmesh%NeE
217 !--
218 ke2d = lmesh%EMap3Dto2D(ke)
219 cori(:) = coriolis(elem%IndexH2Dto3D(:),ke2d)
220
221 gsqrtv(:) = lmesh%Gsqrt(:,ke) / lmesh%GsqrtH(elem%IndexH2Dto3D,ke2d)
222 rgsqrtv(:) = 1.0_rp / gsqrtv(:)
223
224 !--
225 rhot_(:) = p0ovr * (pres_hyd(:,ke) * rp0)**rgamm + drhot_(:,ke)
226 ! DPRES_(:) = PRES00 * ( Rtot(:,ke) * rP0 * RHOT_(:) )**( CPtot(:,ke) / CVtot(:,ke) ) &
227 ! - PRES_hyd(:,ke)
228
229 gsqrtdens_(:) = lmesh%Gsqrt(:,ke) * ( ddens_(:,ke) + dens_hyd(:,ke) )
230 rdens_(:) = 1.0_rp / gsqrtdens_(:)
231 u_(:) = momx_(:,ke) * rdens_(:)
232 v_(:) = momy_(:,ke) * rdens_(:)
233 w_(:) = momz_(:,ke) * rdens_(:)
234 wt_(:) = w_(:) * rgsqrtv(:) + lmesh%GI3(:,ke,1) * u_(:) + lmesh%GI3(:,ke,2) * v_(:)
235 pot_(:) = rhot_(:) * rdens_(:)
236
237 ke2d = lmesh%EMap3Dto2D(ke)
238 cori(:) = coriolis(elem%IndexH2Dto3D(:),ke2d)
239
240 !-- Gradient hydrostatic pressure
241
242 dpres_hyd(:) = pres_hyd(:,ke) - pres_hyd_ref(:,ke)
243
244 call sparsemat_matmul(dx, gsqrtv(:) * dpres_hyd(:), fx)
245 call sparsemat_matmul(dz, gsqrtv(:) * lmesh%GI3(:,ke,1) * dpres_hyd(:), fz)
246 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux_hyd(:,ke,1), liftdelflx)
247 gradphyd_x(:) = lmesh%Escale(:,ke,1,1) * fx(:) &
248 + lmesh%Escale(:,ke,3,3) * fz(:) &
249 + liftdelflx(:)
250
251 call sparsemat_matmul(dy, gsqrtv(:) * dpres_hyd(:), fy)
252 call sparsemat_matmul(dz, gsqrtv(:) * lmesh%GI3(:,ke,2) * dpres_hyd(:), fz)
253 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux_hyd(:,ke,2), liftdelflx)
254 gradphyd_y(:) = lmesh%Escale(:,ke,2,2) * fy(:) &
255 + lmesh%Escale(:,ke,3,3) * fz(:) &
256 + liftdelflx(:)
257
258 !-- DENS
259 call dx_ab( dxt1d_, gsqrtdens_(:), u_, elem%Nnode_h1D, elem%Nnode_v, fx_sp )
260 call dy_ab( dyt1d_, gsqrtdens_(:), v_, elem%Nnode_h1D, elem%Nnode_v, fy_sp )
261 call dz_ab( dzt1d_, gsqrtdens_(:), wt_(:) - w_(:) * rgsqrtv(:), elem%Nnode_h1D, elem%Nnode_v, fz_sp )
262 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,dens_vid), liftdelflx)
263
264 dens_dt(:,ke) = - ( &
265 lmesh%Escale(:,ke,1,1) * fx_sp(:) &
266 + lmesh%Escale(:,ke,2,2) * fy_sp(:) &
267 + lmesh%Escale(:,ke,3,3) * fz_sp(:) &
268 + liftdelflx(:) )
269
270 !-- MOMX
271 call dx_abc( dxt1d_, gsqrtdens_, u_, u_, elem%Nnode_h1D, elem%Nnode_v, fx_sp )
272 call dy_abc( dyt1d_, gsqrtdens_, u_, v_, elem%Nnode_h1D, elem%Nnode_v, fy_sp )
273 call dz_abc( dzt1d_, gsqrtdens_, u_, wt_, elem%Nnode_h1D, elem%Nnode_v, fz_sp )
274 call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * dpres_(:,ke) , fx)
275 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momx_vid), liftdelflx)
276
277 momx_dt(:,ke) = &
278 - ( lmesh%Escale(:,ke,1,1) * ( fx_sp(:) + fx(:) ) &
279 + lmesh%Escale(:,ke,2,2) * fy_sp(:) &
280 + lmesh%Escale(:,ke,3,3) * fz_sp(:) &
281 + liftdelflx(:) ) / lmesh%Gsqrt(:,ke) &
282 - gradphyd_x(:) * rgsqrtv(:) &
283 + cori(:) * momy_(:,ke)
284
285 !-- MOMY
286 call dx_abc( dxt1d_, gsqrtdens_, v_, u_, elem%Nnode_h1D, elem%Nnode_v, fx_sp )
287 call dy_abc( dyt1d_, gsqrtdens_, v_, v_, elem%Nnode_h1D, elem%Nnode_v, fy_sp )
288 call dz_abc( dzt1d_, gsqrtdens_, v_, wt_, elem%Nnode_h1D, elem%Nnode_v, fz_sp )
289 call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * dpres_(:,ke) , fy)
290 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momy_vid), liftdelflx)
291
292 momy_dt(:,ke) = &
293 - ( lmesh%Escale(:,ke,1,1) * fx_sp(:) &
294 + lmesh%Escale(:,ke,2,2) * ( fy_sp(:) + fy(:) ) &
295 + lmesh%Escale(:,ke,3,3) * fz_sp(:) &
296 + liftdelflx(:) ) / lmesh%Gsqrt(:,ke) &
297 - gradphyd_y(:) * rgsqrtv(:) &
298 - cori(:) * momx_(:,ke)
299
300 !-- MOMZ
301 call dx_abc( dxt1d_, gsqrtdens_, w_, u_, elem%Nnode_h1D, elem%Nnode_v, fx_sp )
302 call dy_abc( dyt1d_, gsqrtdens_, w_, v_, elem%Nnode_h1D, elem%Nnode_v, fy_sp )
303 call dz_abc( dzt1d_, gsqrtdens_, w_, wt_, elem%Nnode_h1D, elem%Nnode_v, fz_sp )
304 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momz_vid), liftdelflx)
305
306 momz_dt(:,ke) = - ( &
307 lmesh%Escale(:,ke,1,1) * fx_sp(:) &
308 + lmesh%Escale(:,ke,2,2) * fy_sp(:) &
309 + lmesh%Escale(:,ke,3,3) * fz_sp(:) &
310 + liftdelflx(:) )
311
312 !-- RHOT
313 call dx_abc( dxt1d_, gsqrtdens_, pot_, u_, elem%Nnode_h1D, elem%Nnode_v, fx_sp )
314 call dy_abc( dyt1d_, gsqrtdens_, pot_, v_, elem%Nnode_h1D, elem%Nnode_v, fy_sp )
315 call dz_abc( dzt1d_, gsqrtdens_, pot_, wt_(:) - w_(:) * rgsqrtv(:), elem%Nnode_h1D, elem%Nnode_v, fz_sp )
316 call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,rhot_vid), liftdelflx)
317
318 rhot_dt(:,ke) = - ( &
319 lmesh%Escale(:,ke,1,1) * fx_sp(:) &
320 + lmesh%Escale(:,ke,2,2) * fy_sp(:) &
321 + lmesh%Escale(:,ke,3,3) * fz_sp(:) &
322 + liftdelflx(:) )
323 end do
324 call prof_rapend( 'cal_dyn_tend_interior', 3)
325
326 return
module FElib / Fluid dyn solver / Atmosphere / Nonhydrostatic model / HEVI / Numflux
subroutine, public atm_dyn_dgm_nonhydro3d_rhot_hevi_numflux_get_generalvc_asis(del_flux, del_flux_hyd, ddens_, momx_, momy_, momz_, drhot_, dpres_, dens_hyd, pres_hyd, rtot, cvtot, cptot, gsqrt, g13, g23, nx, ny, nz, vmapm, vmapp, lmesh, elem, lmesh2d, elem2d)

References scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_numflux::atm_dyn_dgm_nonhydro3d_rhot_hevi_numflux_get_generalvc_asis().

◆ atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform_cal_vi()

subroutine, public scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform::atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform_cal_vi ( real(rp), dimension(elem%np,lmesh%nea), intent(out) dens_dt,
real(rp), dimension(elem%np,lmesh%nea), intent(out) momx_dt,
real(rp), dimension(elem%np,lmesh%nea), intent(out) momy_dt,
real(rp), dimension(elem%np,lmesh%nea), intent(out) momz_dt,
real(rp), dimension(elem%np,lmesh%nea), intent(out) rhot_dt,
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) ddens0_,
real(rp), dimension(elem%np,lmesh%nea), intent(in) momx0_,
real(rp), dimension(elem%np,lmesh%nea), intent(in) momy0_,
real(rp), dimension(elem%np,lmesh%nea), intent(in) momz0_,
real(rp), dimension(elem%np,lmesh%nea), intent(in) drhot0_,
real(rp), dimension(elem%np,lmesh%nea), intent(in) rtot,
real(rp), dimension(elem%np,lmesh%nea), intent(in) cvtot,
real(rp), dimension(elem%np,lmesh%nea), intent(in) cptot,
class(elementoperationbase3d), intent(in) element3d_operation,
class(sparsemat), intent(in) dz,
class(sparsemat), intent(in) lift,
real(rp), intent(in) impl_fac,
real(rp), intent(in) dt,
class(localmesh3d), intent(in) lmesh,
class(elementbase3d), intent(in) elem,
class(localmesh2d), intent(in) lmesh2d,
class(elementbase2d), intent(in) elem2d )

Definition at line 331 of file scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform.F90.

339
347
348 implicit none
349
350 class(LocalMesh3D), intent(in) :: lmesh
351 class(ElementBase3D), intent(in) :: elem
352 class(LocalMesh2D), intent(in) :: lmesh2D
353 class(ElementBase2D), intent(in) :: elem2D
354 real(RP), intent(out) :: DENS_dt(elem%Np,lmesh%NeA)
355 real(RP), intent(out) :: MOMX_dt(elem%Np,lmesh%NeA)
356 real(RP), intent(out) :: MOMY_dt(elem%Np,lmesh%NeA)
357 real(RP), intent(out) :: MOMZ_dt(elem%Np,lmesh%NeA)
358 real(RP), intent(out) :: RHOT_dt(elem%Np,lmesh%NeA)
359 real(RP), intent(in) :: DDENS_(elem%Np,lmesh%NeA)
360 real(RP), intent(in) :: MOMX_(elem%Np,lmesh%NeA)
361 real(RP), intent(in) :: MOMY_(elem%Np,lmesh%NeA)
362 real(RP), intent(in) :: MOMZ_(elem%Np,lmesh%NeA)
363 real(RP), intent(in) :: DRHOT_(elem%Np,lmesh%NeA)
364 real(RP), intent(in) :: DENS_hyd(elem%Np,lmesh%NeA)
365 real(RP), intent(in) :: PRES_hyd(elem%Np,lmesh%NeA)
366 real(RP), intent(in) :: DDENS0_(elem%Np,lmesh%NeA)
367 real(RP), intent(in) :: MOMX0_(elem%Np,lmesh%NeA)
368 real(RP), intent(in) :: MOMY0_(elem%Np,lmesh%NeA)
369 real(RP), intent(in) :: MOMZ0_(elem%Np,lmesh%NeA)
370 real(RP), intent(in) :: DRHOT0_(elem%Np,lmesh%NeA)
371 real(RP), intent(in) :: Rtot(elem%Np,lmesh%NeA)
372 real(RP), intent(in) :: CVtot(elem%Np,lmesh%NeA)
373 real(RP), intent(in) :: CPtot(elem%Np,lmesh%NeA)
374 class(ElementOperationBase3D), intent(in) :: element3D_operation
375 class(SparseMat), intent(in) :: Dz, Lift
376 real(RP), intent(in) :: impl_fac
377 real(RP), intent(in) :: dt
378
379 real(RP) :: PROG_VARS (elem%Np,lmesh%NeZ,PRGVAR_NUM,lmesh%NeX*lmesh%NeY)
380 real(RP) :: PROG_VARS0(elem%Np,lmesh%NeZ,PRGVAR_NUM,lmesh%NeX*lmesh%NeY)
381 real(RP) :: b1D(elem%Nnode_v,3,lmesh%NeZ,elem%Nnode_h1D**2,lmesh%NeX*lmesh%NeY)
382 integer :: ipiv(elem%Nnode_v*3*lmesh%NeZ,elem%Nnode_h1D**2)
383 real(RP) :: b1D_uv(elem%Nnode_v,lmesh%NeZ,2,elem%Nnode_h1D**2,lmesh%NeX*lmesh%NeY)
384 integer :: ipiv_uv(elem%Nnode_v*1*lmesh%NeZ,elem%Nnode_h1D**2)
385 real(RP) :: alph(elem%NfpTot,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
386 real(RP) :: Rtot_z(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
387 real(RP) :: CPtot_ov_CVtot(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
388 real(RP) :: DENS_hyd_z(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
389 real(RP) :: PRES_hyd_z(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
390 real(RP) :: GnnM_z(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
391 real(RP) :: G13_z(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
392 real(RP) :: G23_z(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
393 real(RP) :: GsqrtV_z(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
394 real(RP) :: nz(elem%NfpTot,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
395 integer :: vmapM(elem%NfpTot,lmesh%NeZ)
396 integer :: vmapP(elem%NfpTot,lmesh%NeZ)
397 integer :: ColMask(elem%Nnode_v)
398 integer :: ke_xy, ke_z, ke, ke2D, v
399 integer :: itr_nlin
400 integer :: kl, ku, nz_1D
401 integer :: kl_uv, ku_uv, nz_1D_uv
402 integer :: ij, info
403 logical :: is_converged
404
405 real(RP), allocatable :: PmatBnd(:,:,:)
406 real(RP), allocatable :: PmatBnd_uv(:,:,:)
407 !------------------------------------------------------------------------
408
409
410 call prof_rapstart( 'hevi_cal_vi_prep', 3)
411
412 nz_1d = elem%Nnode_v * 3 * lmesh%NeZ
413 kl = ( elem%Nnode_v + 1 ) * 3 - 1
414 ku = kl
415 nz_1d_uv = elem%Nnode_v * 1 * lmesh%NeZ
416 kl_uv = elem%Nnode_v
417 ku_uv = kl_uv
418 allocate( pmatbnd(2*kl+ku+1,nz_1d,elem%Nnode_h1D**2) )
419 allocate( pmatbnd_uv(2*kl_uv+ku_uv+1,nz_1d_uv,elem%Nnode_h1D**2) )
420
421 call lmesh%GetVmapZ1D( vmapm, vmapp ) ! (out)
422
423 !-
424
425 !$omp parallel private( ke_xy, ke_z, ke, ke2D )
426 !$omp do collapse(2)
427 do ke_xy=1, lmesh%NeX*lmesh%NeY
428 do ke_z=1, lmesh%NeZ
429 ke = ke_xy + (ke_z-1)*lmesh%NeX*lmesh%NeY
430 ke2d = lmesh%EMap3Dto2D(ke)
431
432 prog_vars(:,ke_z,dens_vid,ke_xy) = ddens0_(:,ke)
433 prog_vars(:,ke_z,momx_vid,ke_xy) = momx0_(:,ke)
434 prog_vars(:,ke_z,momy_vid,ke_xy) = momy0_(:,ke)
435 prog_vars(:,ke_z,momz_vid,ke_xy) = momz0_(:,ke)
436 prog_vars(:,ke_z,rhot_vid,ke_xy) = drhot0_(:,ke)
437
438 dens_hyd_z(:,ke_z,ke_xy) = dens_hyd(:,ke)
439 pres_hyd_z(:,ke_z,ke_xy) = pres_hyd(:,ke)
440
441 rtot_z(:,ke_z,ke_xy) = rtot(:,ke)
442 cptot_ov_cvtot(:,ke_z,ke_xy) = cptot(:,ke) / cvtot(:,ke)
443
444 nz(:,ke_z,ke_xy) = lmesh%normal_fn(:,ke,3)
445 g13_z(:,ke_z,ke_xy) = lmesh%GI3(:,ke,1)
446 g23_z(:,ke_z,ke_xy) = lmesh%GI3(:,ke,2)
447 gsqrtv_z(:,ke_z,ke_xy) = lmesh%Gsqrt(:,ke) / lmesh%GsqrtH(elem%IndexH2Dto3D,ke2d)
448
449 gnnm_z(:,ke_z,ke_xy) = ( 1.0_rp / gsqrtv_z(:,ke_z,ke_xy)**2 &
450 + g13_z(:,ke_z,ke_xy)**2 + g23_z(:,ke_z,ke_xy) )
451 end do
452 end do
453 !$omp end do
454 !$omp workshare
455 prog_vars0(:,:,:,:) = prog_vars(:,:,:,:)
456 !$omp end workshare
457 !$omp end parallel
458
459 call prof_rapend( 'hevi_cal_vi_prep', 3)
460
461 !--
462
463 if ( abs(impl_fac) > 0.0_rp ) then
464 call prof_rapstart( 'hevi_cal_vi_itr', 3)
465
466 ! G = (q^n+1 - q^n*) + impl_fac * A(q^n+1) = 0
467 ! dG/dq^n+1 del[q] = - G(q^n*)
468 do itr_nlin = 1, 1
469 call prof_rapstart( 'hevi_cal_vi_ax', 3)
470
471 call vi_eval_ax_uv( &
472 momx_dt(:,:), momy_dt(:,:), alph(:,:,:), & ! (out)
473 prog_vars, prog_vars0, & ! (in)
474 ddens_, momx_, momy_, momz_, drhot_, & ! (in)
475 dens_hyd_z, pres_hyd_z, & ! (in)
476 rtot_z, cptot_ov_cvtot, & ! (in)
477 dz, lift, intrpmat_vpordm1, & ! (in)
478 gnnm_z, g13_z, g23_z, gsqrtv_z, & ! (in)
479 impl_fac, dt, & ! (in)
480 lmesh, elem, nz, vmapm, vmapp, & ! (in)
481 b1d_uv(:,:,:,:,:) ) ! (out)
482
483 call prof_rapend( 'hevi_cal_vi_ax', 3)
484
485 do ke_xy=1, lmesh%NeX * lmesh%NeY
486 call prof_rapstart( 'hevi_cal_vi_matbnd', 3)
487
488 call vi_construct_matbnd_uv( pmatbnd_uv(:,:,:), & ! (out)
489 kl_uv, ku_uv, nz_1d_uv, & ! (in)
490 prog_vars(:,:,:,ke_xy), & ! (in)
491 dens_hyd_z(:,:,ke_xy), pres_hyd_z(:,:,ke_xy), & ! (in)
492 g13_z(:,:,ke_xy), g23_z(:,:,ke_xy), gsqrtv_z(:,:,ke_xy), & ! (in)
493 alph(:,:,ke_xy), & ! (in)
494 rtot_z(:,:,ke_xy), cptot_ov_cvtot(:,:,ke_xy), & ! (in)
495 dz, lift, intrpmat_vpordm1, & ! (in)
496 impl_fac, dt, & ! (in)
497 lmesh, elem, nz(:,:,ke_xy), vmapm, vmapp, ke_xy, 1 ) ! (in)
498
499 call prof_rapend( 'hevi_cal_vi_matbnd', 3)
500
501 call prof_rapstart( 'hevi_cal_vi_lin', 3)
502 !$omp parallel private(ij, v, ke_z, info, ColMask)
503 !$omp do
504 do ij=1, elem%Nnode_h1D**2
505! call dgbsv( nz_1D_uv, kl_uv, ku_uv, 2, PmatBnd_uv(:,:,ij), 2*kl_uv+ku_uv+1, ipiv_uv(:,ij), b1D_uv(:,:,:,ij,ke_xy), nz_1D_uv, info)
506 call linalgebra_solvelineq_bndmat( pmatbnd_uv(:,:,ij), b1d_uv(:,:,:,ij,ke_xy), ipiv_uv(:,ij), nz_1d_uv, kl_uv, ku_uv, 2, vi_use_lapack_flag )
507
508 colmask(:) = elem%Colmask(:,ij)
509 do ke_z=1, lmesh%NeZ
510 prog_vars(colmask(:),ke_z,momx_vid,ke_xy) = prog_vars(colmask(:),ke_z,momx_vid,ke_xy) + b1d_uv(:,ke_z,1,ij,ke_xy)
511 prog_vars(colmask(:),ke_z,momy_vid,ke_xy) = prog_vars(colmask(:),ke_z,momy_vid,ke_xy) + b1d_uv(:,ke_z,2,ij,ke_xy)
512 end do
513 end do ! for ij
514 !$omp end do
515 !$omp end parallel
516 call prof_rapend( 'hevi_cal_vi_lin', 3)
517
518 end do ! for ke_xy
519
520 call prof_rapstart( 'hevi_cal_vi_ax', 3)
521 call vi_eval_ax( &
522 dens_dt(:,:), momz_dt(:,:), rhot_dt(:,:), & ! (out, dummy)
523 alph(:,:,:), & ! (in)
524 prog_vars, prog_vars0, & ! (in)
525 ddens_, momx_, momy_, momz_, drhot_, & ! (in)
526 dens_hyd_z, pres_hyd_z, & ! (in)
527 rtot_z, cptot_ov_cvtot, & ! (in)
528 dz, lift, intrpmat_vpordm1, & ! (in)
529 gnnm_z, g13_z, g23_z, gsqrtv_z, & ! (in)
530 impl_fac, dt, & ! (in)
531 lmesh, elem, nz, vmapm, vmapp, & ! (in)
532 b1d(:,:,:,:,:) ) ! (out)
533 call prof_rapend( 'hevi_cal_vi_ax', 3)
534
535 do ke_xy=1, lmesh%NeX * lmesh%NeY
536 call prof_rapstart( 'hevi_cal_vi_matbnd', 3)
537 call vi_construct_matbnd( pmatbnd(:,:,:), & ! (out)
538 kl, ku, nz_1d, & ! (in)
539 prog_vars(:,:,:,ke_xy), & ! (in)
540 dens_hyd_z(:,:,ke_xy), pres_hyd_z(:,:,ke_xy), & ! (in)
541 g13_z(:,:,ke_xy), g23_z(:,:,ke_xy), gsqrtv_z(:,:,ke_xy), & ! (in)
542 alph(:,:,ke_xy), & ! (in)
543 rtot_z(:,:,ke_xy), cptot_ov_cvtot(:,:,ke_xy), & ! (in)
544 dz, lift, intrpmat_vpordm1, & ! (in)
545 impl_fac, dt, & ! (in)
546 lmesh, elem, nz(:,:,ke_xy), vmapm, vmapp, ke_xy, 1 ) ! (in)
547
548 call prof_rapend( 'hevi_cal_vi_matbnd', 3)
549
550 call prof_rapstart( 'hevi_cal_vi_lin', 3)
551 !$omp parallel private(ij, v, ke_z, info, ColMask)
552 !$omp do
553 do ij=1, elem%Nnode_h1D**2
554! call dgbsv( nz_1D, kl, ku, 1, PmatBnd(:,:,ij), 2*kl+ku+1, ipiv(:,ij), b1D(:,:,:,ij,ke_xy), nz_1D, info)
555 call linalgebra_solvelineq_bndmat( pmatbnd(:,:,ij), b1d(:,:,:,ij,ke_xy), ipiv(:,ij), nz_1d, kl, ku, 1, vi_use_lapack_flag )
556
557 colmask(:) = elem%Colmask(:,ij)
558 do ke_z=1, lmesh%NeZ
559 prog_vars(colmask(:),ke_z,dens_vid,ke_xy) = prog_vars(colmask(:),ke_z,dens_vid,ke_xy) + b1d(1,:,ke_z,ij,ke_xy)
560 prog_vars(colmask(:),ke_z,momz_vid,ke_xy) = prog_vars(colmask(:),ke_z,momz_vid,ke_xy) + b1d(2,:,ke_z,ij,ke_xy)
561 prog_vars(colmask(:),ke_z,rhot_vid,ke_xy) = prog_vars(colmask(:),ke_z,rhot_vid,ke_xy) + b1d(3,:,ke_z,ij,ke_xy)
562 end do
563 end do ! for ij
564 !$omp end do
565 !$omp end parallel
566 call prof_rapend( 'hevi_cal_vi_lin', 3)
567
568 end do ! for ke_xy
569 end do ! itr nlin
570
571 call prof_rapend( 'hevi_cal_vi_itr', 3)
572 end if
573
574
575 call prof_rapstart( 'hevi_cal_vi_retrun_var', 3)
576 if ( abs(impl_fac) > 0.0_rp) then
577 !$omp parallel do collapse(2) private(ke_xy, ke_z, ke)
578 do ke_xy=1, lmesh%NeX * lmesh%NeY
579 do ke_z=1, lmesh%NeZ
580 ke = ke_xy + (ke_z-1)*lmesh%NeX*lmesh%NeY
581 dens_dt(:,ke) = ( prog_vars(:,ke_z,dens_vid,ke_xy) - ddens_(:,ke) ) / impl_fac
582 momx_dt(:,ke) = ( prog_vars(:,ke_z,momx_vid,ke_xy) - momx_(:,ke) ) / impl_fac
583 momy_dt(:,ke) = ( prog_vars(:,ke_z,momy_vid,ke_xy) - momy_(:,ke) ) / impl_fac
584 momz_dt(:,ke) = ( prog_vars(:,ke_z,momz_vid,ke_xy) - momz_(:,ke) ) / impl_fac
585 rhot_dt(:,ke) = ( prog_vars(:,ke_z,rhot_vid,ke_xy) - drhot_(:,ke) ) / impl_fac
586 end do
587 end do
588 else
589 call vi_eval_ax_uv( &
590 momx_dt(:,:), momy_dt(:,:), & ! (out)
591 alph(:,:,:), & ! (out, dummy)
592 prog_vars, prog_vars0, & ! (in)
593 ddens_, momx_, momy_, momz_, drhot_, & ! (in)
594 dens_hyd_z, pres_hyd_z, & ! (in)
595 rtot_z, cptot_ov_cvtot, & ! (in)
596 dz, lift, intrpmat_vpordm1, & ! (in)
597 gnnm_z, g13_z, g23_z, gsqrtv_z, & ! (in)
598 impl_fac, dt, & ! (in)
599 lmesh, elem, nz, vmapm, vmapp ) ! (in)
600
601 call vi_eval_ax( &
602 dens_dt(:,:), momz_dt(:,:), rhot_dt(:,:), & ! (out)
603 alph(:,:,:), & ! (in, dummy)
604 prog_vars, prog_vars0, & ! (in)
605 ddens_, momx_, momy_, momz_, drhot_, & ! (in)
606 dens_hyd_z, pres_hyd_z, & ! (in)
607 rtot_z, cptot_ov_cvtot, & ! (in)
608 dz, lift, intrpmat_vpordm1, & ! (in)
609 gnnm_z, g13_z, g23_z, gsqrtv_z, & ! (in)
610 impl_fac, dt, & ! (in)
611 lmesh, elem, nz, vmapm, vmapp ) ! (in)
612 end if
613 call prof_rapend( 'hevi_cal_vi_retrun_var', 3)
614
615 return
subroutine, public atm_dyn_dgm_nonhydro3d_rhot_hevi_common_eval_ax(dens_t, momz_t, rhot_t, alph, prog_vars, prog_vars0, ddens00, momx00, momy00, momz00, drhot00, dens_hyd, pres_hyd, rtot, cptot_ov_cvtot, dz, lift, intrpmat_vpordm1, gnnm, g13, g23, gsqrtv, impl_fac, dt, lmesh, elem, nz, vmapm, vmapp, b1d_ij)
subroutine, public atm_dyn_dgm_nonhydro3d_rhot_hevi_common_eval_ax_uv(momx_t, momy_t, alph, prog_vars, prog_vars0, ddens00, momx00, momy00, momz00, drhot00, dens_hyd, pres_hyd, rtot, cptot_ov_cvtot, dz, lift, intrpmat_vpordm1, gnnm, g13, g23, gsqrtv, impl_fac, dt, lmesh, elem, nz, vmapm, vmapp, b1d_ij_uv)
subroutine, public atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd(pmatbnd, kl, ku, nz_1d, prog_vars0, dens_hyd, pres_hyd, g13, g23, gsqrtv, alph, rtot, cptot_ov_cvtot, dz, lift, intrpmat_vpordm1, impl_fac, dt, lmesh, elem, nz, vmapm, vmapp, ke_x, ke_y)
subroutine, public atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd_uv(pmatbnd_uv, kl_uv, ku_uv, nz_1d_uv, prog_vars0, dens_hyd, pres_hyd, g13, g23, gsqrtv, alph, rtot, cptot_ov_cvtot, dz, lift, intrpmat_vpordm1, impl_fac, dt, lmesh, elem, nz, vmapm, vmapp, ke_x, ke_y)
Module common / Linear algebra.
subroutine, public linalgebra_solvelineq_bndmat(a, b, ipiv, n, kl, ku, nrhs, use_lapack)
Calculate the solution of linear equations with a band matrix.

References scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common::atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd(), scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common::atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd_uv(), scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common::atm_dyn_dgm_nonhydro3d_rhot_hevi_common_eval_ax(), scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common::atm_dyn_dgm_nonhydro3d_rhot_hevi_common_eval_ax_uv(), scale_atm_dyn_dgm_nonhydro3d_common::intrpmat_vpordm1, scale_linalgebra::linalgebra_solvelineq_bndmat(), and scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common::vi_use_lapack_flag.