FE-Project
Loading...
Searching...
No Matches
scale_atm_phy_bl_dgm_common.F90
Go to the documentation of this file.
1!> module FElib / Atmosphere / Physics / boundary layer turbulence
2!!
3!! @par Description
4!! Boundary layer turbulence process
5!!
6!! @author Yuta Kawai, Team SCALE
7!<
8!-------------------------------------------------------------------------------
9#include "scaleFElib.h"
10module scale_atm_phy_bl_dgm_common
11 !-----------------------------------------------------------------------------
12 !
13 !++ Used modules
14 !
15 use scale_precision
16 use scale_io
17 use scale_prc
18 use scale_prof
19 use scale_const, only: &
20 eps => const_eps
21 use scale_tracer, only: qa
22
23 use scale_sparsemat
24 use scale_element_base, only: &
25 elementbase1d, elementbase2d, elementbase3d
26 use scale_element_hexahedral, only: hexahedralelement
27 use scale_localmesh_2d, only: localmesh2d
28 use scale_localmesh_3d, only: localmesh3d
29 use scale_mesh_base3d, only: meshbase3d
30 use scale_localmeshfield_base, only: localmeshfield3d, localmeshfieldbaselist
31 use scale_meshfield_base, only: meshfield3d
32
33 use scale_element_operation_base, only: elementoperationbase3d
34
35 !-----------------------------------------------------------------------------
36 implicit none
37 private
38 !-----------------------------------------------------------------------------
39 !
40 !++ Public procedures
41 !
42 public :: atm_phy_bl_dgm_common_calc_tendency
43
44 !-----------------------------------------------------------------------------
45 !
46 !++ Private procedure
47 !
48
49 !-----------------------------------------------------------------------------
50 !
51 !++ Private parameters & variables
52 !
53contains
54 !> Calculate tendency with PBL turbulence models
55!OCL SERIAL
56 subroutine atm_phy_bl_dgm_common_calc_tendency( &
57 RHOU_tp, RHOV_tp, DRHOT_tp, RHOQ_tp_list, & ! (out)
58 ddens_, momx_, momy_, drhot_, qtrc_list, & ! (in)
59 pt_, dens_hyd, pres_hyd, nu, kh, & ! (in)
60 element3d_operation, c_ip, dtsec, & ! (in)
61 lmesh, elem, elem1d, is_bound, & ! (in)
62 use_delta_form ) ! (in)
63 use scale_atm_dyn_dgm_hevi_common_linalgebra, only: &
64 atm_dyn_dgm_hevi_common_linalgebra_get_param
65 implicit none
66 class(localmesh3d), intent(in), target :: lmesh
67 class(elementbase3d), intent(in) :: elem
68 class(elementbase1d), intent(in) :: elem1d
69 real(rp), intent(out) :: rhou_tp(elem%np,lmesh%nea)
70 real(rp), intent(out) :: rhov_tp(elem%np,lmesh%nea)
71 real(rp), intent(out) :: drhot_tp(elem%np,lmesh%nea)
72 type(localmeshfieldbaselist), intent(inout) :: rhoq_tp_list(qa)
73 real(rp), intent(in) :: ddens_(elem%np,lmesh%nea)
74 real(rp), intent(in) :: momx_(elem%np,lmesh%nea)
75 real(rp), intent(in) :: momy_(elem%np,lmesh%nea)
76 real(rp), intent(in) :: drhot_(elem%np,lmesh%nea)
77 type(localmeshfieldbaselist), intent(in) :: qtrc_list(qa)
78 real(rp), intent(in) :: pt_(elem%np,lmesh%nea)
79 real(rp), intent(in) :: dens_hyd(elem%np,lmesh%nea)
80 real(rp), intent(in) :: pres_hyd(elem%np,lmesh%nea)
81 real(rp), intent(in) :: nu(elem%np,lmesh%nea)
82 real(rp), intent(in) :: kh(elem%np,lmesh%nea)
83 class(elementoperationbase3d), intent(in) :: element3d_operation
84 real(rp), intent(in) :: c_ip
85 real(rp), intent(in) :: dtsec
86 logical, intent(in) :: is_bound(elem%nfptot,lmesh%ne)
87 logical, intent(in) :: use_delta_form
88
89 class(localmesh2d), pointer :: lmesh2d
90 class(elementbase2d), pointer :: elem2d
91
92 integer :: iq
93 real(rp) :: qtrc00_(elem%np,qa,lmesh%ne)
94
95 real(rp) :: prog_vars (elem%np,lmesh%nex*lmesh%ney,lmesh%nez,3+qa)
96 real(rp) :: alph_m(elem%nfptot,lmesh%ne)
97 real(rp) :: alph_h(elem%nfptot,lmesh%ne)
98 real(rp) :: gsqrtv(elem%np,lmesh%ne)
99
100 integer :: vmapm(elem%nfptot,lmesh%ne)
101 integer :: vmapp(elem%nfptot,lmesh%ne)
102 integer :: ke_xy, ke_z, ke, ke2d
103 integer :: p
104
105 integer :: im, jm
106
107 real(rp) :: dens(elem%np,lmesh%ne)
108 real(rp), allocatable :: b1d_ij(:,:,:,:,:)
109 real(rp), allocatable :: bndmatl(:,:,:,:,:)
110 real(rp), allocatable :: bndmatd(:,:,:,:,:)
111 real(rp), allocatable :: g(:,:,:,:,:,:)
112
113 real(rp) :: impl_fac, r_impl_fac
114 !------------------------------------------------------------------------
115
116 lmesh2d => lmesh%lcmesh2D
117 elem2d => lmesh2d%refElem2D
118 impl_fac = 1.0_rp * dtsec
119
120 call lmesh%GetVmapZ3D( vmapm, vmapp ) ! (out)
121 call atm_dyn_dgm_hevi_common_linalgebra_get_param( im, jm, & ! (out)
122 elem%Nnode_v )
123
124 allocate( b1d_ij(im*elem%Nnode_v,3+qa,jm,lmesh%Ne2D,lmesh%NeZ) )
125 allocate( bndmatd(im*elem%Nnode_v,elem%Nnode_v,jm,lmesh%Ne2D,2) )
126 allocate( bndmatl(im*elem%Nnode_v,elem%Nnode_v,jm,lmesh%Ne2D,2) )
127 allocate( g(im*elem%Nnode_v,elem%Nnode_v,jm,lmesh%Ne2D,lmesh%NeZ,2) )
128
129 !$omp parallel private( ke, ke2D )
130 !$omp do collapse(2)
131 do iq = 1, qa
132 do ke=lmesh%NeS, lmesh%NeE
133 qtrc00_(:,iq,ke) = qtrc_list(iq)%ptr%val(:,ke)
134 end do
135 end do
136 !$omp end do
137
138 !$omp do collapse(2)
139 do ke_z =1, lmesh%NeZ
140 do ke_xy=1, lmesh%NeX * lmesh%NeY
141 ke = ke_xy + (ke_z-1)*lmesh%Ne2D
142 ke2d = lmesh%EMap3Dto2D(ke)
143
144 dens(:,ke) = dens_hyd(:,ke) + ddens_(:,ke)
145
146 prog_vars(:,ke_xy,ke_z,1) = momx_(:,ke)
147 prog_vars(:,ke_xy,ke_z,2) = momy_(:,ke)
148 prog_vars(:,ke_xy,ke_z,3) = dens(:,ke) * pt_(:,ke)
149 do iq = 1, qa
150 prog_vars(:,ke_xy,ke_z,3+iq) = dens(:,ke) * qtrc00_(:,iq,ke)
151 end do
152
153 do p=1, elem%Np
154 gsqrtv(p,ke) = lmesh%Gsqrt(p,ke) / lmesh%GsqrtH(elem%IndexH2Dto3D(p),ke2d)
155 end do
156 end do
157 end do
158 !$omp end do
159 !$omp end parallel
160
161 call eval_ax( rhou_tp, rhov_tp, drhot_tp, rhoq_tp_list, alph_m, alph_h, & !(out)
162 prog_vars, momx_, momy_, pt_, qtrc00_, nu, kh, dens, gsqrtv, & !(in)
163 impl_fac, dtsec, lmesh, elem, vmapm, vmapp, is_bound, & !(in)
164 element3d_operation, c_ip, im, jm, b1d_ij, use_delta_form ) !(in)
165
166 call vi_solve( prog_vars, & ! (inout)
167 bndmatl, bndmatd, g, b1d_ij, & ! (inout)
168 dens, nu, kh, gsqrtv, c_ip, dtsec, impl_fac, & ! (in)
169 im, jm, lmesh, elem, elem1d, use_delta_form ) ! (in)
170
171 !---
172 r_impl_fac = 1.0_rp / impl_fac
173
174 !$omp parallel do collapse(2) private( ke, iq )
175 do ke_z =1, lmesh%NeZ
176 do ke_xy=1, lmesh%NeX * lmesh%NeY
177 ke = ke_xy + (ke_z-1)*lmesh%Ne2D
178 rhou_tp(:,ke) = ( prog_vars(:,ke_xy,ke_z,1) - momx_(:,ke) ) * r_impl_fac
179 rhov_tp(:,ke) = ( prog_vars(:,ke_xy,ke_z,2) - momy_(:,ke) ) * r_impl_fac
180 drhot_tp(:,ke) = ( prog_vars(:,ke_xy,ke_z,3) - dens(:,ke) * pt_(:,ke) ) * r_impl_fac
181
182 do iq=1, qa
183 rhoq_tp_list(iq)%ptr%val(:,ke) = ( prog_vars(:,ke_xy,ke_z,3+iq) - dens(:,ke) * qtrc00_(:,iq,ke) ) * r_impl_fac
184 end do
185 end do
186 end do
187
188 return
189 end subroutine atm_phy_bl_dgm_common_calc_tendency
190
191!- Private subroutines -----------------------------
192
193 !> Solve the block tridiagonal system which is generated by the vertical diffusion equation for PBL turbulence models
194!OCL SERIAL
195 subroutine vi_solve( PROG_VARS, & ! (inout)
196 bndmatl, bndmatd, g, b1d_ij, & ! (inout)
197 dens, nu, kh, gsqrtv, c_ip, dtsec, impl_fac, & ! (in)
198 im, jm, lmesh, elem, elem1d, use_delta_form ) ! (in)
199 implicit none
200 integer, intent(in) :: im, jm
201 class(localmesh3d), intent(in) :: lmesh
202 class(elementbase3d), intent(in) :: elem
203 class(elementbase1d), intent(in) :: elem1d
204 real(rp), intent(inout) :: prog_vars(elem%nnode_h1d**2,elem%nnode_v,lmesh%ne2d,lmesh%nez,3+qa)
205 real(rp), intent(inout) :: bndmatl(im*elem%nnode_v,elem%nnode_v,jm,lmesh%ne2d,2)
206 real(rp), intent(inout) :: bndmatd(im*elem%nnode_v,elem%nnode_v,jm,lmesh%ne2d,2)
207 real(rp), intent(inout) :: g(im*elem%nnode_v,elem%nnode_v,jm,lmesh%ne2d,lmesh%nez,2)
208 real(rp), intent(inout) :: b1d_ij(im*elem%nnode_v,3+qa,jm,lmesh%ne2d,lmesh%nez)
209 real(rp), intent(in) :: dens(elem%np,lmesh%ne)
210 real(rp), intent(in) :: nu(elem%np,lmesh%nea)
211 real(rp), intent(in) :: kh(elem%np,lmesh%nea)
212 real(rp), intent(in) :: gsqrtv(elem%np,lmesh%ne)
213 real(rp), intent(in) :: c_ip
214 real(rp), intent(in) :: dtsec
215 real(rp), intent(in) :: impl_fac
216 logical, intent(in) :: use_delta_form
217
218 integer :: ke_z, ke_xy
219 integer :: i, j, ij
220 integer :: pv, pv1, pv2, pp, p2
221 integer :: iv
222
223 real(rp) :: tmp(im*elem%nnode_v,2)
224 real(rp) :: tmp_b(im*elem%nnode_v,2)
225 real(rp) :: tmp_b2(im*elem%nnode_v)
226 !-------------------------------------------------------
227
228 do ke_z=1, lmesh%NeZ
229 call construct_matbnd_sip( &
230 bndmatl(:,:,:,:,1), bndmatd(:,:,:,:,1), g(:,:,:,:,ke_z,1), & ! (out)
231 dens, nu, gsqrtv, elem1d%Dx1, elem1d%M, elem1d%invM, & ! (in)
232 impl_fac, dtsec, c_ip, lmesh, elem, im, jm, ke_z ) ! (in)
233
234 call construct_matbnd_sip( &
235 bndmatl(:,:,:,:,2), bndmatd(:,:,:,:,2), g(:,:,:,:,ke_z,2), & ! (out)
236 dens, kh, gsqrtv, elem1d%Dx1, elem1d%M, elem1d%invM, & ! (in)
237 impl_fac, dtsec, c_ip, lmesh, elem, im, jm, ke_z ) ! (in)
238
239 if ( ke_z > 1 ) then
240 !$omp parallel do collapse(2) private( pp, p2, tmp, tmp_b,tmp_b2 )
241 do ke_xy=1, lmesh%NeX * lmesh%NeY
242 do j=1, jm
243 !* D_k <- D_k - L_k * G_{k-1} ------------------
244 do pv2=1, elem%Nnode_v
245 tmp(:,:) = 0.0_rp
246 do pv=1, elem%Nnode_v
247 do pv1=1, elem%Nnode_v
248 do i=1, im
249 pp = i + (pv1-1)*im; p2 = i + (pv -1)*im
250 tmp(pp,1) = tmp(pp,1) + bndmatl(pp,pv,j,ke_xy,1) * g(p2,pv2,j,ke_xy,ke_z-1,1)
251 tmp(pp,2) = tmp(pp,2) + bndmatl(pp,pv,j,ke_xy,2) * g(p2,pv2,j,ke_xy,ke_z-1,2)
252 end do
253 end do
254 end do
255 bndmatd(:,pv2,j,ke_xy,1) = bndmatd(:,pv2,j,ke_xy,1) - tmp(:,1)
256 bndmatd(:,pv2,j,ke_xy,2) = bndmatd(:,pv2,j,ke_xy,2) - tmp(:,2)
257 end do ! loop for pv2
258
259 !* b_k <- b_k - L_k * b_{k-1} ------------------
260 tmp_b(:,:) = 0.0_rp
261 do pv=1, elem%Nnode_v
262 do pv1=1, elem%Nnode_v
263 do i=1, im
264 pp = i + (pv1-1)*im; p2 = i + (pv -1)*im
265 tmp_b(pp,1) = tmp_b(pp,1) + bndmatl(pp,pv,j,ke_xy,1) * b1d_ij(p2,1,j,ke_xy,ke_z-1)
266 tmp_b(pp,2) = tmp_b(pp,2) + bndmatl(pp,pv,j,ke_xy,1) * b1d_ij(p2,2,j,ke_xy,ke_z-1)
267 end do
268 end do
269 end do
270 b1d_ij(:,1,j,ke_xy,ke_z) = b1d_ij(:,1,j,ke_xy,ke_z) - tmp_b(:,1)
271 b1d_ij(:,2,j,ke_xy,ke_z) = b1d_ij(:,2,j,ke_xy,ke_z) - tmp_b(:,2)
272
273 do iv=3, 3+qa
274 tmp_b2(:) = 0.0_rp
275 do pv=1, elem%Nnode_v
276 do pv1=1, elem%Nnode_v
277 do i=1, im
278 pp = i + (pv1-1)*im; p2 = i + (pv -1)*im
279 tmp_b2(pp) = tmp_b2(pp) + bndmatl(pp,pv,j,ke_xy,2) * b1d_ij(p2,iv,j,ke_xy,ke_z-1)
280 end do
281 end do
282 end do
283 b1d_ij(:,iv,j,ke_xy,ke_z) = b1d_ij(:,iv,j,ke_xy,ke_z) - tmp_b2(:)
284 end do
285
286 end do
287 end do
288 end if
289
290 call linkernel_solve_sip( b1d_ij(:,:,:,:,ke_z), g(:,:,:,:,ke_z,1), g(:,:,:,:,ke_z,2),&
291 bndmatd(:,:,:,:,1), bndmatd(:,:,:,:,2), elem%Nnode_v, im, jm, lmesh%Ne2D, ke_z == lmesh%NeZ )
292
293 end do
294 do ke_z=lmesh%NeZ-1, 1, -1
295 !$omp parallel do collapse(2) private( pp, p2, tmp_b,tmp_b2 )
296 do ke_xy=1, lmesh%NeX * lmesh%NeY
297 do j=1, jm
298 tmp_b(:,:) = 0.0_rp
299 do pv2=1, elem%Nnode_v
300 do pv =1, elem%Nnode_v
301 do i=1, im
302 pp = i + (pv -1)*im; p2 = i + (pv2-1)*im
303 tmp_b(pp,1) = tmp_b(pp,1) + g(pp,pv2,j,ke_xy,ke_z,1) * b1d_ij(p2,1,j,ke_xy,ke_z+1)
304 tmp_b(pp,2) = tmp_b(pp,2) + g(pp,pv2,j,ke_xy,ke_z,1) * b1d_ij(p2,2,j,ke_xy,ke_z+1)
305 end do
306 end do
307 end do
308 b1d_ij(:,1,j,ke_xy,ke_z) = b1d_ij(:,1,j,ke_xy,ke_z) - tmp_b(:,1)
309 b1d_ij(:,2,j,ke_xy,ke_z) = b1d_ij(:,2,j,ke_xy,ke_z) - tmp_b(:,2)
310
311 do iv=3, 3+qa
312 tmp_b2(:) = 0.0_rp
313 do pv2=1, elem%Nnode_v
314 do pv =1, elem%Nnode_v
315 do i=1, im
316 pp = i + (pv -1)*im; p2 = i + (pv2-1)*im
317 tmp_b2(pp) = tmp_b2(pp) + g(pp,pv2,j,ke_xy,ke_z,2) * b1d_ij(p2,iv,j,ke_xy,ke_z+1)
318 end do
319 end do
320 end do
321 b1d_ij(:,iv,j,ke_xy,ke_z) = b1d_ij(:,iv,j,ke_xy,ke_z) - tmp_b2(:)
322 end do
323
324 end do
325 end do
326 end do
327
328 !$omp parallel do collapse(2) private( p2, pp )
329 do ke_z=1, lmesh%NeZ
330 do ke_xy=1, lmesh%NeX * lmesh%NeY
331 do j=1, jm
332 if ( use_delta_form ) then
333 do iv=1, 3+qa
334 do pv=1, elem%Nnode_v
335 do i=1, im
336 p2 = i + (pv-1)*im; pp = i + (j-1)*im
337 prog_vars(pp,pv,ke_xy,ke_z,iv) = prog_vars(pp,pv,ke_xy,ke_z,iv) + b1d_ij(p2,iv,j,ke_xy,ke_z)
338 end do
339 end do
340 end do
341 else
342 do iv=1, 3+qa
343 do pv=1, elem%Nnode_v
344 do i=1, im
345 p2 = i + (pv-1)*im; pp = i + (j-1)*im
346 prog_vars(pp,pv,ke_xy,ke_z,iv) = b1d_ij(p2,iv,j,ke_xy,ke_z)
347 end do
348 end do
349 end do
350 end if
351 end do
352 end do
353 end do
354
355 return
356 end subroutine vi_solve
357
358 !> Solve the linear system associated with the diagonal block of the block tridiagonal system
359!OCL SERIAL
360 subroutine linkernel_solve_sip( b, G1, G2, & ! (inout)
361 bndmatd1, bndmatd2, nv, im, jm, ne2d, top_flag ) ! (in)
362 implicit none
363
364 integer, intent(in) :: nv
365 integer, intent(in) :: im
366 integer, intent(in) :: jm
367 integer, intent(in) :: ne2d
368 real(rp), intent(inout) :: b(im*nv,3+qa,jm,ne2d)
369 real(rp), intent(inout) :: g1(im*nv,nv,jm,ne2d)
370 real(rp), intent(inout) :: g2(im*nv,nv,jm,ne2d)
371 real(rp), intent(in) :: bndmatd1(im*nv,nv,jm,ne2d)
372 real(rp), intent(in) :: bndmatd2(im*nv,nv,jm,ne2d)
373 logical, intent(in) :: top_flag
374
375 integer :: ke_xy
376 integer :: i, j
377 integer :: pv1, pv2
378 integer :: iv
379 integer :: pp
380
381 integer :: info
382 integer :: nrhs1, nrhs2
383 integer :: ipiv1(nv), ipiv2(nv)
384
385 real(rp) :: amat1(nv,nv)
386 real(rp) :: rhs1(nv,nv+2)
387 real(rp) :: amat2(nv,nv)
388 real(rp) :: rhs2(nv,nv+1+qa)
389 !------------------------------------------------------------
390
391 !$omp parallel do collapse(2) &
392 !$omp private(ke_xy,j,i,pv1,pv2,iv,pp,info, &
393 !$omp nrhs1,nrhs2,ipiv1,ipiv2,Amat1,RHS1,Amat2,RHS2 )
394 do ke_xy=1, ne2d
395 do j=1, jm
396 do i=1, im
397
398 do pv2=1, nv
399 do pv1=1, nv
400 pp = i + (pv1-1)*im
401 amat1(pv1,pv2) = bndmatd1(pp,pv2,j,ke_xy)
402 amat2(pv1,pv2) = bndmatd2(pp,pv2,j,ke_xy)
403 end do
404 end do
405
406 rhs1(:,:) = 0.0_rp
407 rhs2(:,:) = 0.0_rp
408
409 if ( top_flag ) then
410 nrhs1 = 2
411 do iv=1, 2
412 do pv1=1, nv
413 pp = i + (pv1-1)*im
414 rhs1(pv1,iv) = b(pp,iv,j,ke_xy)
415 end do
416 end do
417
418 nrhs2 = 1 + qa
419 do iv=1, 1+qa
420 do pv1=1, nv
421 pp = i + (pv1-1)*im
422 rhs2(pv1,iv) = b(pp,2+iv,j,ke_xy)
423 end do
424 end do
425 else
426 nrhs1 = nv + 2
427 nrhs2 = nv + 1 + qa
428
429 do pv2=1, nv
430 do pv1=1, nv
431 pp = i + (pv1-1)*im
432 rhs1(pv1,pv2) = g1(pp,pv2,j,ke_xy)
433 rhs2(pv1,pv2) = g2(pp,pv2,j,ke_xy)
434 end do
435 end do
436 do iv=1, 2
437 do pv1=1, nv
438 pp = i + (pv1-1)*im
439 rhs1(pv1,nv+iv) = b(pp,iv,j,ke_xy)
440 end do
441 end do
442 do iv=1, 1+qa
443 do pv1=1, nv
444 pp = i + (pv1-1)*im
445 rhs2(pv1,nv+iv) = b(pp,2+iv,j,ke_xy)
446 end do
447 end do
448
449 end if
450
451 ! LU factorization:
452 call dgetrf( nv, nv, amat1, nv, ipiv1, info )
453 if ( info /= 0 ) then
454 log_error('linkernel_solve_sip',*) 'NU, DGETRF failed: info=', info, ', i=', i, ', j=', j, ', ke_xy=', ke_xy
455 call prc_abort
456 end if
457 call dgetrf( nv, nv, amat2, nv, ipiv2, info )
458 if ( info /= 0 ) then
459 log_error('linkernel_solve_sip',*) 'KH, DGETRF failed: info=', info, ', i=', i, ', j=', j, ', ke_xy=', ke_xy
460 call prc_abort
461 end if
462
463 ! Solve all right-hand sides using the same LU factors.
464 call dgetrs( 'N', nv, nrhs1, amat1, nv, ipiv1, rhs1, nv, info )
465 if ( info /= 0 ) then
466 log_error('linkernel_solve_sip',*) 'NU, DGETRS failed: info=', info, ', i=', i, ', j=', j, ', ke_xy=', ke_xy
467 call prc_abort
468 end if
469 call dgetrs( 'N', nv, nrhs2, amat2, nv, ipiv2, rhs2, nv, info )
470 if ( info /= 0 ) then
471 log_error('linkernel_solve_sip',*) 'KH, DGETRS failed: info=', info, ', i=', i, ', j=', j, ', ke_xy=', ke_xy
472 call prc_abort
473 end if
474
475 if ( top_flag ) then
476 ! At the top level RHS(:,1:2) contains D^{-1} b.
477 do iv=1, 2
478 do pv1=1, nv
479 pp = i + (pv1-1)*im
480 b(pp,iv,j,ke_xy) = rhs1(pv1,iv)
481 end do
482 end do
483 do iv=1, 1+qa
484 do pv1=1, nv
485 pp = i + (pv1-1)*im
486 b(pp,2+iv,j,ke_xy) = rhs2(pv1,iv)
487 end do
488 end do
489
490 ! Not used in backward substitution.
491 do pv2=1, nv
492 do pv1=1, nv
493 pp = i + (pv1-1)*im
494 g1(pp,pv2,j,ke_xy) = 0.0_rp
495 g2(pp,pv2,j,ke_xy) = 0.0_rp
496 end do
497 end do
498
499 else
500 ! RHS(:,1:nv) now contains D^{-1} U.
501 do pv2=1, nv
502 do pv1=1, nv
503 pp = i + (pv1-1)*im
504 g1(pp,pv2,j,ke_xy) = rhs1(pv1,pv2)
505 g2(pp,pv2,j,ke_xy) = rhs2(pv1,pv2)
506 end do
507 end do
508 ! RHS1(:,nv+1:nv+2) contains D^{-1} b.
509 do iv=1, 2
510 do pv1=1, nv
511 pp = i + (pv1-1)*im
512 b(pp,iv,j,ke_xy) = rhs1(pv1,nv+iv)
513 end do
514 end do
515 ! RHS2(:,nv+1:nv+1+QA) contains D^{-1} b.
516 do iv=1, 1+qa
517 do pv1=1, nv
518 pp = i + (pv1-1)*im
519 b(pp,2+iv,j,ke_xy) = rhs2(pv1,nv+iv)
520 end do
521 end do
522 end if
523
524 end do ! i
525 end do ! j
526 end do ! ke_xy
527 !$omp end parallel do
528
529 return
530 end subroutine linkernel_solve_sip
531
532 !> Construct the block tridiagonal matrix for the vertical diffusion equation
533!OCL SERIAL
534 subroutine construct_matbnd_sip( BndMatL, BndMatD, BndMatU, & ! (out)
535 rho, kdiff, gsqrtv, dx1d, m1d,invm1d, & ! (in)
536 impl_fac, dt, penalty_fac, lmesh, elem, im, jm, ke_z ) ! (in)
537 implicit none
538 class(localmesh3d), intent(in) :: lmesh
539 class(elementbase3d), intent(in) :: elem
540 integer, intent(in) :: im, jm
541 real(rp), intent(out) :: bndmatl(im,elem%nnode_v,elem%nnode_v,jm,lmesh%ne2d)
542 real(rp), intent(out) :: bndmatd(im,elem%nnode_v,elem%nnode_v,jm,lmesh%ne2d)
543 real(rp), intent(out) :: bndmatu(im,elem%nnode_v,elem%nnode_v,jm,lmesh%ne2d)
544 real(rp), intent(in) :: rho(elem%np,lmesh%ne)
545 real(rp), intent(in) :: kdiff(elem%np,lmesh%nea)
546 real(rp), intent(in) :: gsqrtv(elem%np,lmesh%ne)
547 real(rp), intent(in) :: dx1d(elem%nnode_v,elem%nnode_v)
548 real(rp), intent(in) :: m1d(elem%nnode_v,elem%nnode_v)
549 real(rp), intent(in) :: invm1d(elem%nnode_v,elem%nnode_v)
550 real(rp), intent(in) :: impl_fac
551 real(rp), intent(in) :: dt
552 real(rp), intent(in) :: penalty_fac
553 integer, intent(in) :: ke_z
554
555 integer :: ke2d, ke, p
556 integer :: i, j
557 integer :: pv, pv1, pv2
558 integer :: f
559 integer :: ke_nb, ke_z_nb
560 integer :: pvm, pvp
561
562 real(rp) :: lambda
563
564 real(rp) :: mu_loc (elem%nnode_v)
565 real(rp) :: rinv_loc(elem%nnode_v)
566 real(rp) :: rgsqrtv_loc(elem%nnode_v)
567
568 real(rp) :: mu_nb (elem%nnode_v)
569 real(rp) :: rinv_nb(elem%nnode_v)
570
571 real(rp) :: avol (elem%nnode_v,elem%nnode_v)
572 real(rp) :: affmm(elem%nnode_v,elem%nnode_v)
573 real(rp) :: affmp(elem%nnode_v,elem%nnode_v)
574
575 real(rp) :: minvaloc(elem%nnode_v,elem%nnode_v)
576 real(rp) :: minvanb (elem%nnode_v,elem%nnode_v)
577
578 real(rp) :: dz(elem%nnode_v,elem%nnode_v)
579 real(rp) :: dz_loc(elem%nnode_v), dz_nb(elem%nnode_v)
580
581 real(rp) :: mu_face_m, mu_face_p
582 real(rp) :: sigma
583 real(rp) :: fscale_m
584
585 logical :: boundary_flag
586
587 integer :: nnode_v
588 !---------------------------------------------------------------------------
589
590 call prof_rapstart('phy_bl_cal_vi_matbnd_sip', 3)
591
592 lambda = impl_fac
593 nnode_v = elem%Nnode_v
594
595 !$omp parallel do collapse(2) private( ke, p, ke_nb, ke_z_nb, pvM, pvP, &
596 !$omp mu_loc, rinv_loc, rGsqrtV_loc, mu_nb, rinv_nb, &
597 !$omp Avol, AffMM, AffMP, MinvALoc, MinvANb, Dz, Dz_loc, Dz_nb, &
598 !$omp mu_face_M, mu_face_P, sigma, Fscale_M, boundary_flag )
599 do ke2d = 1, lmesh%Ne2D
600 do j = 1, jm
601 ke = ke2d + (ke_z-1)*lmesh%Ne2D
602
603 bndmatl(:,:,:,j,ke2d) = 0.0_rp
604 bndmatd(:,:,:,j,ke2d) = 0.0_rp
605 bndmatu(:,:,:,j,ke2d) = 0.0_rp
606
607 do i=1, im
608 do pv=1, nnode_v
609 p = i + (j-1)*im + (pv-1)*im*jm
610
611 mu_loc(pv) = rho(p,ke) * kdiff(p,ke)
612 rinv_loc(pv) = 1.0_rp / rho(p,ke)
613 rgsqrtv_loc(pv) = 1.0_rp / gsqrtv(p,ke)
614
615 avol(:,pv) = 0.0_rp
616 bndmatd(i,pv,pv,j,ke2d) = 1.0_rp
617 end do
618 do pv2=1, nnode_v
619 do pv1=1, nnode_v
620 p = i + (j-1)*im + (pv1-1)*im*jm
621 dz(pv1,pv2) = lmesh%Escale(p,ke,3,3) * dx1d(pv1,pv2)
622 end do
623 end do
624
625 !- Volume contribution ---
626
627 ! Avol = Dz^T M diag(mu) Dz diag(1/rho)
628 do pv2=1, nnode_v
629 do pv1=1, nnode_v
630 do pv=1, nnode_v
631 avol(pv1,pv2) = avol(pv1,pv2) &
632 + dz(pv,pv1) * m1d(pv,pv) * mu_loc(pv) * rgsqrtv_loc(pv) * dz(pv,pv2) * rinv_loc(pv2)
633 end do
634 end do
635 end do
636
637 ! Minv Avol
638 minvaloc(:,:) = 0.0_rp
639 do pv2=1, nnode_v
640 do pv=1, nnode_v
641 do pv1=1, nnode_v
642 minvaloc(pv1,pv2) = minvaloc(pv1,pv2) + invm1d(pv1,pv) * avol(pv,pv2)
643 end do
644 end do
645 minvaloc(pv1,pv2) = rgsqrtv_loc(pv1) * minvaloc(pv1,pv2)
646 end do
647
648 do pv2=1, nnode_v
649 do pv1=1, nnode_v
650 bndmatd(i,pv1,pv2,j,ke2d) = bndmatd(i,pv1,pv2,j,ke2d) + lambda * minvaloc(pv1,pv2)
651 end do
652 end do
653
654 !- Bottom and top faces
655
656 do f=1, 2
657 if (f == 1) then
658 ke_z_nb = max(ke_z-1,1)
659 pvm = 1; pvp = nnode_v
660 boundary_flag = (ke_z == 1)
661 else
662 ke_z_nb = min(ke_z+1,lmesh%NeZ)
663 pvm = nnode_v; pvp = 1
664 boundary_flag = (ke_z == lmesh%NeZ)
665 end if
666 ke_nb = ke2d + (ke_z_nb-1)*lmesh%Ne2D
667
668 ! Homogeneous Neumann boundary condition:
669 ! mu dphi/dn = 0
670 ! No SIP boundary contribution is added here.
671
672 if (boundary_flag) cycle
673
674 do pv=1, nnode_v
675 p = i + (j-1)*im + (pv-1)*im*jm
676 mu_nb(pv) = rho(p,ke_nb) * kdiff(p,ke_nb)
677 rinv_nb(pv) = 1.0_rp / rho(p,ke_nb)
678 end do
679
680 mu_face_m = mu_loc(pvm)
681 mu_face_p = mu_nb(pvp)
682 sigma = penalty_fac * real(nnode_v, kind=rp)**2 * max(mu_face_m, mu_face_p)
683
684
685 !--
686
687 dz_loc(:) = dz(pvm,:)
688 p = i + (j-1)*im + (pvp-1)*im*jm
689 do pv=1, nnode_v
690 dz_nb(pv) = lmesh%Escale(p,ke_nb,3,3) * dx1d(pvp,pv)
691 end do
692
693 call construct_sip_face_blocks_lgl( affmm, affmp, & ! (out)
694 dz_loc, dz_nb, mu_face_m, mu_face_p, rinv_loc, rinv_nb, & ! (in)
695 sigma, pvm, pvp, f, elem%Nnode_v ) ! (in)
696
697 fscale_m = lmesh%Fscale(elem%Nfp_h*elem%Nfaces_h+1,ke)
698 affmm(:,:) = fscale_m * affmm(:,:)
699 affmp(:,:) = fscale_m * affmp(:,:)
700
701 ! M^-1 AffMM, M^-1 AffMP
702 minvaloc(:,:) = 0.0_rp; minvanb(:,:) = 0.0_rp
703 do pv2=1, nnode_v
704 do pv=1, nnode_v
705 do pv1=1, nnode_v
706 minvaloc(pv1,pv2) = minvaloc(pv1,pv2) + invm1d(pv1,pv) * affmm(pv,pv2)
707 minvanb(pv1,pv2) = minvanb(pv1,pv2) + invm1d(pv1,pv) * affmp(pv,pv2)
708 end do
709 end do
710 end do
711
712 do pv2=1, nnode_v
713 do pv1=1, nnode_v
714 bndmatd(i,pv1,pv2,j,ke2d) = bndmatd(i,pv1,pv2,j,ke2d) + lambda * rgsqrtv_loc(pv1) * minvaloc(pv1,pv2)
715 if ( f == 1 ) then
716 bndmatl(i,pv1,pv2,j,ke2d) = bndmatl(i,pv1,pv2,j,ke2d) + lambda * rgsqrtv_loc(pv1) * minvanb(pv1,pv2)
717 else
718 bndmatu(i,pv1,pv2,j,ke2d) = bndmatu(i,pv1,pv2,j,ke2d) + lambda * rgsqrtv_loc(pv1) * minvanb(pv1,pv2)
719 end if
720 end do
721 end do
722
723 end do ! end face loop
724 end do
725 end do
726 end do
727
728 call prof_rapend('phy_bl_cal_vi_matbnd_sip', 3)
729 return
730 end subroutine construct_matbnd_sip
731
732!OCL SERIAL
733 subroutine construct_sip_face_blocks_lgl( &
734 AffMM, AffMP, & ! (out)
735 dm, dp, mum, mup, rinvm, rinvp, & ! (in)
736 sigma, pvm, pvp, face_id, nv ) ! (in)
737 implicit none
738 integer, intent(in) :: nv
739 real(rp), intent(out) :: affmm(nv,nv)
740 real(rp), intent(out) :: affmp(nv,nv)
741 real(rp), intent(in) :: dm(nv), dp(nv)
742 real(rp), intent(in) :: mum, mup
743 real(rp), intent(in) :: rinvm(nv), rinvp(nv)
744 real(rp), intent(in) :: sigma
745 integer, intent(in) :: pvm, pvp
746 integer, intent(in) :: face_id
747
748 integer :: i, j
749 real(rp) :: ncom
750 !----------------------------------------------------
751
752 if (face_id == 1) then
753 ncom = -1.0_rp
754 else
755 ncom = +1.0_rp
756 end if
757
758 affmm(:,:) = 0.0_rp; affmp(:,:) = 0.0_rp
759 do i=1, nv
760 do j=1, nv
761 if (i == pvm) then
762 affmm(i,j) = affmm(i,j) &
763 - 0.5_rp * mum * ncom * dm(j) * rinvm(j)
764 end if
765 if (j == pvm) then
766 affmm(i,j) = affmm(i,j) &
767 - 0.5_rp * mum * ncom * dm(i) * rinvm(j)
768 end if
769 if (i == pvm .and. j == pvm) then
770 affmm(i,j) = affmm(i,j) + sigma * rinvm(j)
771 end if
772 !-
773 if ( i == pvm ) then
774 affmp(i,j) = affmp(i,j) &
775 - 0.5_rp * mup * ncom * dp(j) * rinvp(j)
776 end if
777 if ( j == pvp ) then
778 affmp(i,j) = affmp(i,j) &
779 + 0.5_rp * mum * ncom * dm(i) * rinvp(j)
780 end if
781 if ( i == pvm .and. j == pvp ) then
782 affmp(i,j) = affmp(i,j) - sigma * rinvp(j)
783 end if
784 end do
785 end do
786 return
787 end subroutine construct_sip_face_blocks_lgl
788
789!OCL SERIAL
790 subroutine eval_ax( MOMX_t, MOMY_t, DRHOT_t, RHOQ_t_list, alph_M, alph_H, &
791 PROG_VARS, MOMX00, MOMY00, PT00, QTRC00, NU, KH, DENS, GsqrtV, impl_fac, dt, &
792 lmesh, elem, vmapM, vmapP, is_bound, element3D_operation, C_IP, &
793 im, jm, b, use_delta_form )
794 implicit none
795 class(localmesh3d), intent(in) :: lmesh
796 class(elementbase3d), intent(in) :: elem
797 real(rp), intent(out) :: momx_t(elem%np,lmesh%ne)
798 real(rp), intent(out) :: momy_t(elem%np,lmesh%ne)
799 real(rp), intent(out) :: drhot_t(elem%np,lmesh%ne)
800 type(localmeshfieldbaselist), intent(inout) :: rhoq_t_list(qa)
801 real(rp), intent(out) :: alph_m(elem%nfptot,lmesh%ne)
802 real(rp), intent(out) :: alph_h(elem%nfptot,lmesh%ne)
803 real(rp), intent(in) :: prog_vars(elem%np,lmesh%nex*lmesh%ney*lmesh%nez,3+qa)
804 real(rp), intent(in) :: momx00(elem%np,lmesh%nea)
805 real(rp), intent(in) :: momy00(elem%np,lmesh%nea)
806 real(rp), intent(in) :: pt00(elem%np,lmesh%nea)
807 real(rp), intent(in) :: qtrc00(elem%np,qa,lmesh%ne)
808 real(rp), intent(in) :: nu(elem%np,lmesh%nea)
809 real(rp), intent(in) :: kh(elem%np,lmesh%nea)
810 real(rp), intent(in) :: dens(elem%np,lmesh%ne)
811 real(rp), intent(in) :: gsqrtv(elem%np,lmesh%ne)
812 real(rp), intent(in) :: impl_fac
813 real(rp), intent(in) :: dt
814 integer, intent(in) :: vmapm(elem%nfptot,lmesh%ne)
815 integer, intent(in) :: vmapp(elem%nfptot,lmesh%ne)
816 logical, intent(in) :: is_bound(elem%nfptot,lmesh%ne)
817 class(elementoperationbase3d), intent(in) :: element3d_operation
818 real(rp), intent(in) :: c_ip
819 integer, intent(in) :: im, jm
820 real(rp), intent(out), optional :: b(im,elem%nnode_v,3+qa,jm,lmesh%ne)
821 logical, intent(in), optional :: use_delta_form
822
823 real(rp) :: flux(elem%np,3), dflux(elem%np,2,3)
824 real(rp) :: flux_q(elem%np), dflux_q(elem%np,2)
825 real(rp) :: rdens(elem%np)
826 real(rp) :: rgsqrtv
827 real(rp) :: e33
828
829 real(rp) :: diff_flux_z_broken(elem%np,lmesh%nea,3+qa)
830 real(rp) :: diff_flux_z(elem%np,lmesh%nea,3+qa)
831
832 real(rp) :: del_flux(elem%nfptot,3,lmesh%ne)
833 real(rp) :: del_flux_q(elem%nfptot,qa,lmesh%ne)
834
835 integer :: ke_xy, ke_z
836 integer :: ke, ke2d
837 integer :: p, fp
838 integer :: iv, iq
839
840 integer :: i, j
841 integer :: pv
842 logical :: flag_cal_b
843 logical :: flag_use_delta_form
844 !---------------------------------------------------------
845
846 if ( present(b) .and. present(use_delta_form) ) then
847 flag_cal_b = .true.
848 flag_use_delta_form = use_delta_form
849 else
850 flag_cal_b = .false.
851 flag_use_delta_form = .true.
852 end if
853
854 if ( flag_use_delta_form ) then
855 call cal_grad_del_flux( del_flux, del_flux_q, & ! (out)
856 prog_vars, dens, & ! (in)
857 lmesh%normal_fn(:,:,3), lmesh%Fscale, vmapm, vmapp, & ! (in)
858 lmesh, elem, lmesh%lcmesh2D, lmesh%lcmesh2D%refElem2D ) ! (in)
859
860 !$omp parallel do collapse(2) private( ke,ke2D,iv, DFlux,DFlux_q, RDENS, RGsqrtV, E33 )
861 do ke_z=1, lmesh%NeZ
862 do ke_xy=1, lmesh%NeX*lmesh%NeY
863 ke = ke_xy + (ke_z-1)*lmesh%NeX*lmesh%NeY
864 ke2d = lmesh%EMap3Dto2D(ke)
865
866 rdens(:) = 1.0_rp / dens(:,ke)
867 do iv=1, 3
868 call element3d_operation%Dz( prog_vars(:,ke,iv) * rdens(:), dflux(:,1,iv) )
869 call element3d_operation%Lift( del_flux(:,iv,ke), dflux(:,2,iv) )
870 end do
871
872 do p=1, elem%Np
873 rgsqrtv = 1.0_rp / gsqrtv(p,ke)
874 e33 = lmesh%Escale(p,ke,3,3)
875 diff_flux_z_broken(p,ke,1) = dens(p,ke) * nu(p,ke) * e33 * dflux(p,1,1) * rgsqrtv
876 diff_flux_z_broken(p,ke,2) = dens(p,ke) * nu(p,ke) * e33 * dflux(p,1,2) * rgsqrtv
877 diff_flux_z_broken(p,ke,3) = dens(p,ke) * kh(p,ke) * e33 * dflux(p,1,3) * rgsqrtv
878
879 diff_flux_z(p,ke,1) = dens(p,ke) * nu(p,ke) * ( e33 * dflux(p,1,1) + dflux(p,2,1) ) * rgsqrtv
880 diff_flux_z(p,ke,2) = dens(p,ke) * nu(p,ke) * ( e33 * dflux(p,1,2) + dflux(p,2,2) ) * rgsqrtv
881 diff_flux_z(p,ke,3) = dens(p,ke) * kh(p,ke) * ( e33 * dflux(p,1,3) + dflux(p,2,3) ) * rgsqrtv
882 end do
883
884 do iq=1, qa
885 iv = 3 + iq
886 call element3d_operation%Dz( prog_vars(:,ke,iv) * rdens(:), dflux_q(:,1) )
887 call element3d_operation%Lift( del_flux_q(:,iq,ke), dflux_q(:,2) )
888 do p=1, elem%Np
889 rgsqrtv = 1.0_rp / gsqrtv(p,ke)
890 e33 = lmesh%Escale(p,ke,3,3)
891 diff_flux_z_broken(p,ke,iv) = dens(p,ke) * kh(p,ke) * e33 * dflux_q(p,1) * rgsqrtv
892 diff_flux_z(p,ke,iv) = dens(p,ke) * kh(p,ke) * ( e33 * dflux_q(p,1) + dflux_q(p,2) ) * rgsqrtv
893 end do
894 end do
895
896 end do
897 end do
898
899 !---------------------------------------------------------
900
901 call cal_del_flux( del_flux, del_flux_q, alph_m, alph_h, & ! (out)
902 diff_flux_z, diff_flux_z_broken, prog_vars, dens, nu, kh, c_ip, & ! (in)
903 lmesh%normal_fn(:,:,3), lmesh%Fscale, vmapm, vmapp, is_bound, & ! (in)
904 lmesh, elem, lmesh%lcmesh2D, lmesh%lcmesh2D%refElem2D ) ! (in)
905
906 !$omp parallel do collapse(2) private( ke,iv, ke2D, DFlux,DFlux_q, RGsqrtV, E33 )
907 do ke_z=1, lmesh%NeZ
908 do ke_xy=1, lmesh%NeX*lmesh%NeY
909 ke = ke_xy + (ke_z-1)*lmesh%NeX*lmesh%NeY
910 ke2d = lmesh%EMap3Dto2D(ke)
911 do iv=1, 3
912 call element3d_operation%Dz( diff_flux_z(:,ke,iv), dflux(:,1,iv) )
913 call element3d_operation%Lift( del_flux(:,iv,ke), dflux(:,2,iv) )
914 end do
915
916 do p=1, elem%Np
917 rgsqrtv = 1.0_rp / gsqrtv(p,ke)
918 e33 = lmesh%Escale(p,ke,3,3)
919 momx_t(p,ke) = ( e33 * dflux(p,1,1) + dflux(p,2,1) ) * rgsqrtv
920 momy_t(p,ke) = ( e33 * dflux(p,1,2) + dflux(p,2,2) ) * rgsqrtv
921 drhot_t(p,ke) = ( e33 * dflux(p,1,3) + dflux(p,2,3) ) * rgsqrtv
922 end do
923
924 do iq=1, qa
925 iv = 3 + iq
926 call element3d_operation%Dz( diff_flux_z(:,ke,iv), dflux_q(:,1) )
927 call element3d_operation%Lift( del_flux_q(:,iq,ke), dflux_q(:,2) )
928 do p=1, elem%Np
929 rgsqrtv = 1.0_rp / gsqrtv(p,ke)
930 e33 = lmesh%Escale(p,ke,3,3)
931 rhoq_t_list(iq)%ptr%val(p,ke) = ( e33 * dflux_q(p,1) + dflux_q(p,2) ) * rgsqrtv
932 end do
933 end do
934
935 end do
936 end do
937 end if
938
939 if ( flag_cal_b ) then
940 if ( use_delta_form ) then
941 !$omp parallel do private(p,iv)
942 do ke=1, lmesh%Ne
943 do pv=1, elem%Nnode_v
944
945 do j=1, jm
946 do i=1, im
947 p = i + (j-1)*im + (pv-1)*im*jm
948 b(i,pv,1,j,ke) = impl_fac * momx_t(p,ke) &
949 - prog_vars(p,ke,1) &
950 + momx00(p,ke)
951 b(i,pv,2,j,ke) = impl_fac * momy_t(p,ke) &
952 - prog_vars(p,ke,2) &
953 + momy00(p,ke)
954 b(i,pv,3,j,ke) = impl_fac * drhot_t(p,ke) &
955 - prog_vars(p,ke,3) &
956 + dens(p,ke) * pt00(p,ke)
957 end do
958 end do
959 !-
960 do iq=1, qa
961 iv = 3 + iq
962 do j=1, jm
963 do i=1, im
964 p = i + (j-1)*im + (pv-1)*im*jm
965 b(i,pv,iv,j,ke) = impl_fac * rhoq_t_list(iq)%ptr%val(p,ke) &
966 - prog_vars(p,ke,iv) &
967 + dens(p,ke) * qtrc00(p,iq,ke)
968 end do
969 end do
970 end do
971
972 end do
973 end do
974 else
975 !$omp parallel do private(p,iv)
976 do ke=1, lmesh%Ne
977 do pv=1, elem%Nnode_v
978
979 do j=1, jm
980 do i=1, im
981 p = i + (j-1)*im + (pv-1)*im*jm
982 b(i,pv,1,j,ke) = momx00(p,ke)
983 b(i,pv,2,j,ke) = momy00(p,ke)
984 b(i,pv,3,j,ke) = dens(p,ke) * pt00(p,ke)
985 end do
986 end do
987 !-
988 do iq=1, qa
989 iv = 3 + iq
990 do j=1, jm
991 do i=1, im
992 p = i + (j-1)*im + (pv-1)*im*jm
993 b(i,pv,iv,j,ke) = dens(p,ke) * qtrc00(p,iq,ke)
994 end do
995 end do
996 end do
997
998 end do
999 end do
1000 end if
1001 end if
1002 return
1003 end subroutine eval_ax
1004
1005!OCL SERIAL
1006 subroutine cal_del_flux( del_flux, del_flux_q, alph_M, alph_H, &
1007 DIFF_flux_z, DIFF_flux_z_broken, PROG_VARS, DENS, NU, KH, C_IP, &
1008 nz, Fscale, vmapM, vmapP, is_bound, lmesh, elem, lmesh2D, elem2D )
1009 implicit none
1010 class(localmesh3d), intent(in) :: lmesh
1011 class(elementbase3d), intent(in) :: elem
1012 class(localmesh2d), intent(in) :: lmesh2d
1013 class(elementbase2d), intent(in) :: elem2d
1014 real(rp), intent(out) :: del_flux(elem%nfptot,3,lmesh%ne)
1015 real(rp), intent(out) :: del_flux_q(elem%nfptot,qa,lmesh%ne)
1016 real(rp), intent(out) :: alph_m(elem%nfptot,lmesh%ne)
1017 real(rp), intent(out) :: alph_h(elem%nfptot,lmesh%ne)
1018 real(rp), intent(in) :: diff_flux_z(elem%np*lmesh%nea,3+qa)
1019 real(rp), intent(in) :: diff_flux_z_broken(elem%np*lmesh%nea,3+qa)
1020 real(rp), intent(in) :: prog_vars(elem%np*lmesh%nex*lmesh%ney*lmesh%nez,3+qa)
1021 real(rp), intent(in) :: dens(elem%np*lmesh%nea)
1022 real(rp), intent(in) :: nu(elem%np*lmesh%nea)
1023 real(rp), intent(in) :: kh(elem%np*lmesh%nea)
1024 real(rp), intent(in) :: c_ip
1025 real(rp), intent(in) :: nz(elem%nfptot,lmesh%ne)
1026 real(rp), intent(in) :: fscale(elem%nfptot,lmesh%ne)
1027 integer, intent(in) :: vmapm(elem%nfptot,lmesh%ne)
1028 integer, intent(in) :: vmapp(elem%nfptot,lmesh%ne)
1029 logical, intent(in) :: is_bound(elem%nfptot,lmesh%ne)
1030
1031 integer :: ke, ke_z, ke_xy
1032 integer :: iq
1033
1034 integer :: ip(elem%nfptot), im(elem%nfptot)
1035 real(rp) :: diff_flux_z_p(elem%nfptot,3)
1036
1037 real(rp) :: coef(elem%nfptot)
1038 real(rp) :: rdens_m(elem%nfptot), rdens_p(elem%nfptot)
1039 real(rp) :: numflux(elem%nfptot,3)
1040 real(rp) :: numflux_q(elem%nfptot)
1041 !------------------------------------------------
1042
1043 !$omp parallel do collapse(2) &
1044 !$omp private( ke,iq, iM,iP, coef, RDENS_M, RDENS_P, numflux,numflux_q )
1045 do ke_z=1, lmesh%NeZ
1046 do ke_xy=1, lmesh%NeX*lmesh%NeY
1047 ke = ke_xy + (ke_z-1)*lmesh%Ne2D
1048 im(:) = vmapm(:,ke); ip(:) = vmapp(:,ke)
1049
1050 alph_m(:,ke) = c_ip * (elem%Nnode_v)**2 * max( dens(im)*nu(im), dens(ip)*nu(ip) )
1051 alph_h(:,ke) = c_ip * (elem%Nnode_v)**2 * max( dens(im)*kh(im), dens(ip)*kh(ip) )
1052
1053 coef(:) = nz(:,ke) * fscale(:,ke)
1054 rdens_m(:) = 1.0_rp / dens(im)
1055 rdens_p(:) = 1.0_rp / dens(ip)
1056
1057 where ( is_bound(:,ke) )
1058 numflux(:,1) = 0.0_rp
1059 numflux(:,2) = 0.0_rp
1060 numflux(:,3) = 0.0_rp
1061 elsewhere
1062 numflux(:,1) = 0.5_rp * ( diff_flux_z_broken(ip,1) + diff_flux_z(im,1) )
1063 numflux(:,2) = 0.5_rp * ( diff_flux_z_broken(ip,2) + diff_flux_z(im,2) )
1064 numflux(:,3) = 0.5_rp * ( diff_flux_z_broken(ip,3) + diff_flux_z(im,3) )
1065 end where
1066
1067 !-
1068 del_flux(:,1,ke) = coef(:) * ( numflux(:,1) - diff_flux_z(im,1) ) &
1069 + alph_m(:,ke) * fscale(:,ke) * ( prog_vars(ip,1) * rdens_p(:)- prog_vars(im,1) * rdens_m(:) )
1070
1071 del_flux(:,2,ke) = coef(:) * ( numflux(:,2) - diff_flux_z(im,2) ) &
1072 + alph_m(:,ke) * fscale(:,ke) * ( prog_vars(ip,2) * rdens_p(:) - prog_vars(im,2) * rdens_m(:) )
1073
1074 del_flux(:,3,ke) = coef(:) * ( numflux(:,3) - diff_flux_z(im,3) ) &
1075 + alph_h(:,ke) * fscale(:,ke) * ( prog_vars(ip,3) * rdens_p(:) - prog_vars(im,3) * rdens_m(:) )
1076
1077 do iq=1, qa
1078 where ( is_bound(:,ke) )
1079 numflux_q(:) = 0.0_rp
1080 elsewhere
1081 numflux_q(:) = 0.5_rp * ( diff_flux_z_broken(ip,3+iq) + diff_flux_z(im,3+iq) )
1082 end where
1083
1084 del_flux_q(:,iq,ke) = coef(:) * ( numflux_q(:) - diff_flux_z(im,3+iq) ) &
1085 + alph_h(:,ke) * fscale(:,ke) * ( prog_vars(ip,3+iq) * rdens_p(:) - prog_vars(im,3+iq) * rdens_m(:) )
1086 end do
1087 end do
1088 end do
1089
1090 return
1091 end subroutine cal_del_flux
1092
1093!OCL SERIAL
1094 subroutine cal_grad_del_flux( del_flux, del_flux_q, &
1095 PROG_VARS, DENS, nz, Fscale,vmapM, vmapP, &
1096 lmesh, elem, lmesh2D, elem2D )
1097 implicit none
1098 class(localmesh3d), intent(in) :: lmesh
1099 class(elementbase3d), intent(in) :: elem
1100 class(localmesh2d), intent(in) :: lmesh2d
1101 class(elementbase2d), intent(in) :: elem2d
1102 real(rp), intent(out) :: del_flux(elem%nfptot,3,lmesh%ne)
1103 real(rp), intent(out) :: del_flux_q(elem%nfptot,qa,lmesh%ne)
1104 real(rp), intent(in) :: prog_vars(elem%np*lmesh%nex*lmesh%ney*lmesh%nez,3+qa)
1105 real(rp), intent(in) :: dens(elem%np*lmesh%ne)
1106 real(rp), intent(in) :: nz(elem%nfptot,lmesh%ne)
1107 real(rp), intent(in) :: fscale(elem%nfptot,lmesh%ne)
1108 integer, intent(in) :: vmapm(elem%nfptot,lmesh%ne)
1109 integer, intent(in) :: vmapp(elem%nfptot,lmesh%ne)
1110
1111 integer :: ke, ke_z, ke_xy
1112 integer :: iq
1113 integer :: ip(elem%nfptot), im(elem%nfptot)
1114
1115 real(rp) :: coef(elem%nfptot)
1116 real(rp) :: rdens_m(elem%nfptot), rdens_p(elem%nfptot)
1117 !------------------------------------------------
1118
1119 !$omp parallel do collapse(2) &
1120 !$omp private( ke,iq, iM,iP, coef, RDENS_M, RDENS_P )
1121 do ke_z=1, lmesh%NeZ
1122 do ke_xy=1, lmesh%NeX*lmesh%NeY
1123 ke = ke_xy + (ke_z-1)*lmesh%Ne2D
1124 im(:) = vmapm(:,ke); ip(:) = vmapp(:,ke)
1125
1126 coef(:) = 0.5_rp * nz(:,ke) * fscale(:,ke)
1127 rdens_m(:) = 1.0_rp / dens(im)
1128 rdens_p(:) = 1.0_rp / dens(ip)
1129
1130 del_flux(:,1,ke) = coef(:) * ( prog_vars(ip,1) * rdens_p(:) - prog_vars(im,1) * rdens_m(:) )
1131 del_flux(:,2,ke) = coef(:) * ( prog_vars(ip,2) * rdens_p(:) - prog_vars(im,2) * rdens_m(:) )
1132 del_flux(:,3,ke) = coef(:) * ( prog_vars(ip,3) * rdens_p(:) - prog_vars(im,3) * rdens_m(:) )
1133
1134 do iq=1, qa
1135 del_flux_q(:,iq,ke) = coef(:) * ( prog_vars(ip,3+iq) * rdens_p(:) - prog_vars(im,3+iq) * rdens_m(:) )
1136 end do
1137 end do
1138 end do
1139
1140 return
1141 end subroutine cal_grad_del_flux
1142end module scale_atm_phy_bl_dgm_common
scale_atm_dyn_dgm_hevi_common_linalgebra
module FElib / Fluid dyn solver / Atmosphere / HEVI / Common
Definition scale_atm_dyn_dgm_hevi_common_linalgebra.F90:17
scale_atm_dyn_dgm_hevi_common_linalgebra::atm_dyn_dgm_hevi_common_linalgebra_get_param
subroutine, public atm_dyn_dgm_hevi_common_linalgebra_get_param(im, jm, nnode_h1d)
Definition scale_atm_dyn_dgm_hevi_common_linalgebra.F90:50
scale_atm_phy_bl_dgm_common
module FElib / Atmosphere / Physics / boundary layer turbulence
Definition scale_atm_phy_bl_dgm_common.F90:10
scale_atm_phy_bl_dgm_common::atm_phy_bl_dgm_common_calc_tendency
subroutine, public atm_phy_bl_dgm_common_calc_tendency(rhou_tp, rhov_tp, drhot_tp, rhoq_tp_list, ddens_, momx_, momy_, drhot_, qtrc_list, pt_, dens_hyd, pres_hyd, nu, kh, element3d_operation, c_ip, dtsec, lmesh, elem, elem1d, is_bound, use_delta_form)
Calculate tendency with PBL turbulence models.
Definition scale_atm_phy_bl_dgm_common.F90:63
scale_element_base
module FElib / Element / Base
Definition scale_element_base.F90:11
scale_element_hexahedral
module FElib / Element / hexahedron
Definition scale_element_hexahedral.F90:10
scale_element_operation_base
module FElib / Element / Operation / Base
Definition scale_element_operation_base.F90:12
scale_localmesh_2d
module FElib / Mesh / Local 2D
Definition scale_localmesh_2d.F90:10
scale_localmesh_3d
module FElib / Mesh / Local 3D
Definition scale_localmesh_3d.F90:10
scale_localmeshfield_base
module FElib / Data / base
Definition scale_localmeshfield_base.F90:10
scale_mesh_base3d
module FElib / Mesh / Base 3D
Definition scale_mesh_base3d.F90:10
scale_meshfield_base
module FElib / Data / base
Definition scale_meshfield_base.F90:10
scale_sparsemat
Module common / sparsemat.
Definition scale_sparsemat.F90:11
scale_element_base::elementbase1d
Derived type representing a 1D reference element.
Definition scale_element_base.F90:53
scale_element_base::elementbase2d
Derived type representing a 2D reference element.
Definition scale_element_base.F90:71
scale_element_base::elementbase3d
Derived type representing a 3D reference element.
Definition scale_element_base.F90:111
scale_element_hexahedral::hexahedralelement
Derived type representing a hexahedral element.
Definition scale_element_hexahedral.F90:31
scale_element_operation_base::elementoperationbase3d
Base type for elementwise operations.
Definition scale_element_operation_base.F90:33
scale_localmesh_2d::localmesh2d
Derived type representing a local mesh for 2D domain.
Definition scale_localmesh_2d.F90:30
scale_localmesh_3d::localmesh3d
Derived type to manage a local 3D computational domain.
Definition scale_localmesh_3d.F90:32
scale_localmeshfield_base::localmeshfield3d
Derived type representing a field with 3D local mesh.
Definition scale_localmeshfield_base.F90:69
scale_localmeshfield_base::localmeshfieldbaselist
Definition scale_localmeshfield_base.F90:48
scale_mesh_base3d::meshbase3d
Derived type to manage a computational mesh (base type for 3D domain)
Definition scale_mesh_base3d.F90:37
scale_meshfield_base::meshfield3d
Derived type representing a field with 3D mesh.
Definition scale_meshfield_base.F90:95
Generated by doxygen 1.14.0