FE-Project
Loading...
Searching...
No Matches
scale_atm_dyn_dgm_hydrostatic.F90
Go to the documentation of this file.
1!-------------------------------------------------------------------------------
2!> module FElib / Fluid dyn solver / Atmosphere / Common
3!!
4!! @par Description
5!! Construct hydrostatic state for atmospheric dynamical process.
6!!
7!! @author Yuta Kawai, Team SCALE
8!<
9!-------------------------------------------------------------------------------
10#include "scaleFElib.h"
11module scale_atm_dyn_dgm_hydrostatic
12 !-----------------------------------------------------------------------------
13 !
14 !++ Used modules
15 !
16 use scale_precision
17 use scale_io
18 use scale_prof
19 use scale_prc
20
21 use scale_const, only: &
22 pi => const_pi, &
23 grav => const_grav, &
24 rdry => const_rdry, &
25 cpdry => const_cpdry, &
26 cvdry => const_cvdry, &
27 pres00 => const_pre00
28
29 use scale_gmres, only: gmres
30 use scale_sparsemat, only: &
31 sparsemat, &
32 sparsemat_matmul
33
34 use scale_element_base, only: &
35 elementbase, elementbase3d
36 use scale_element_hexahedral, only: hexahedralelement
37 use scale_localmesh_base, only: localmeshbase
38 use scale_localmesh_3d, only: localmesh3d
39
40 !-----------------------------------------------------------------------------
41 implicit none
42 private
43 !-----------------------------------------------------------------------------
44 !
45 !++ Public procedures
46 !
47 public :: hydrostatic_calc_basicstate_constt
48 public :: hydrostatic_calc_basicstate_constpt
49 public :: hydrostatic_calc_basicstate_consttlaps
50 public :: hydrostatic_calc_basicstate_constptlaps
51 public :: hydrostatic_calc_basicstate_constbvfreq
52 public :: hydrostatic_build_rho_xyz
53
54 interface hydrostatic_build_rho_xyz
55 module procedure hydrostatic_build_rho_xyz_dry
56 module procedure hydrostatic_build_rho_xyz_moist
57 end interface
58
59 !-----------------------------------------------------------------------------
60 !
61 !++ Public parameters & variables
62 !
63
64 !-----------------------------------------------------------------------------
65 !
66 !++ Private procedures & variables
67 !
68 private :: gmres_hydro_core
69 private :: eval_ax, eval_ax_lin
70 private :: cal_del_flux, cal_del_flux_lin
71 private :: construct_pmatinv
72
73contains
74 !> Calculate density and pressure in hydrostatic balance with a constant temperature
75 !!
76 !! @param DENS_hyd hydrostatic density [kg/m3]
77 !! @param PRES_hyd hydrostatic pressure [Pa]
78 !! @param Temp0 Temperature of isothermal atmosphere [K]
79 !! @param PRES_sfc Surface pressure [Pa]
80 !! @param x x-coordinate
81 !! @param y y-coordinate
82 !! @param z z-coordinate [m]
83 !! @param lcmesh3D A object to manage a local 3D mesh
84 !! @param elem A object to manage a 3D finite element
85!OCL SERIAL
86 subroutine hydrostatic_calc_basicstate_constt( &
87 DENS_hyd, PRES_hyd, &
88 Temp0, PRES_sfc, x, y, z, lcmesh3D, elem )
89
90 implicit none
91
92 class(localmesh3d), intent(in) :: lcmesh3d
93 class(elementbase3d), intent(in) :: elem
94 real(rp), intent(out) :: dens_hyd(elem%np,lcmesh3d%nea)
95 real(rp), intent(out) :: pres_hyd(elem%np,lcmesh3d%nea)
96 real(rp), intent(in) :: x(elem%np,lcmesh3d%ne)
97 real(rp), intent(in) :: y(elem%np,lcmesh3d%ne)
98 real(rp), intent(in) :: z(elem%np,lcmesh3d%ne)
99 real(rp), intent(in) :: temp0
100 real(rp), intent(in) :: pres_sfc
101
102 integer :: ke, p
103 real(rp) :: h0
104 !-----------------------------------------------
105
106 h0 = rdry * temp0 / grav
107
108 !$omp parallel do
109 !$acc parallel loop gang present(DENS_hyd, PRES_hyd, z, lcmesh3D, elem)
110 do ke=lcmesh3d%NeS, lcmesh3d%NeE
111 !$acc loop vector
112 do p=1, elem%Np
113 pres_hyd(p,ke) = pres_sfc * exp( - z(p,ke) / h0 )
114 dens_hyd(p,ke) = pres_hyd(p,ke) / ( rdry * temp0 )
115 end do
116 end do
117
118 return
119 end subroutine hydrostatic_calc_basicstate_constt
120
121 !> Calculate density and pressure in hydrostatic balance with a constant potential temperature
122 !!
123 !! @param DENS_hyd hydrostatic density [kg/m3]
124 !! @param PRES_hyd hydrostatic pressure [Pa]
125 !! @param PotTemp0 Constant potential temperature of atmosphere [K]
126 !! @param PRES_sfc Surface pressure [Pa]
127 !! @param x x-coordinate
128 !! @param y y-coordinate
129 !! @param z z-coordinate [m]
130 !! @param lcmesh3D A object to manage a local 3D mesh
131 !! @param elem A object to manage a 3D finite element
132!OCL SERIAL
133 subroutine hydrostatic_calc_basicstate_constpt( &
134 DENS_hyd, PRES_hyd, &
135 PotTemp0, PRES_sfc, x, y, z, lcmesh3D, elem )
136
137 implicit none
138
139 class(localmesh3d), intent(in) :: lcmesh3d
140 class(elementbase3d), intent(in) :: elem
141 real(rp), intent(out) :: dens_hyd(elem%np,lcmesh3d%nea)
142 real(rp), intent(out) :: pres_hyd(elem%np,lcmesh3d%nea)
143 real(rp), intent(in) :: x(elem%np,lcmesh3d%ne)
144 real(rp), intent(in) :: y(elem%np,lcmesh3d%ne)
145 real(rp), intent(in) :: z(elem%np,lcmesh3d%ne)
146 real(rp), intent(in) :: pottemp0
147 real(rp), intent(in) :: pres_sfc
148
149 integer :: ke, p
150
151 real(rp) :: exner
152 real(rp) :: exner_sfc
153 real(rp) :: rovcp
154 real(rp) :: cpovr
155 !-----------------------------------------------
156
157 rovcp = rdry / cpdry
158 cpovr = cpdry / rdry
159 exner_sfc = (pres_sfc / pres00)**rovcp
160
161 !$omp parallel do private(exner)
162 !$acc parallel loop gang present(DENS_hyd, PRES_hyd, z, lcmesh3D, elem)
163 do ke=lcmesh3d%NeS, lcmesh3d%NeE
164 !$acc loop vector
165 do p=1, elem%Np
166 ! Cp * PT0 * d exner / dz = - g
167 exner = exner_sfc - grav / (cpdry * pottemp0) * z(p,ke)
168 pres_hyd(p,ke) = pres00 * exner**cpovr
169 dens_hyd(p,ke) = pres_hyd(p,ke) / ( rdry * exner * pottemp0 )
170 end do
171 end do
172
173 return
174 end subroutine hydrostatic_calc_basicstate_constpt
175
176 !> Calculate density and pressure in hydrostatic balance with a constant Brunt–Väisälä frequency
177 !!
178 !! @param DENS_hyd hydrostatic density [kg/m3]
179 !! @param PRES_hyd hydrostatic pressure [Pa]
180 !! @param BruntVaisalaFreq Constant Brunt–Väisälä frequency [s-1]
181 !! @param PotTemp0 Constant potential temperature of atmosphere [K]
182 !! @param PRES_sfc Surface pressure [Pa]
183 !! @param x x-coordinate
184 !! @param y y-coordinate
185 !! @param z z-coordinate [m]
186 !! @param lcmesh3D A object to manage a local 3D mesh
187 !! @param elem A object to manage a 3D finite element
188!OCL SERIAL
189 subroutine hydrostatic_calc_basicstate_constbvfreq( &
190 DENS_hyd, PRES_hyd, &
191 BruntVaisalaFreq, PotTemp0, PRES_sfc, x, y, z, lcmesh3D, elem )
192
193 implicit none
194
195 class(localmesh3d), intent(in) :: lcmesh3d
196 class(elementbase3d), intent(in) :: elem
197 real(rp), intent(out) :: dens_hyd(elem%np,lcmesh3d%nea)
198 real(rp), intent(out) :: pres_hyd(elem%np,lcmesh3d%nea)
199 real(rp), intent(in) :: x(elem%np,lcmesh3d%ne)
200 real(rp), intent(in) :: y(elem%np,lcmesh3d%ne)
201 real(rp), intent(in) :: z(elem%np,lcmesh3d%ne)
202 real(rp), intent(in) :: bruntvaisalafreq
203 real(rp), intent(in) :: pottemp0
204 real(rp), intent(in) :: pres_sfc
205
206 integer :: ke, p
207
208 real(rp) :: pt
209 real(rp) :: exner
210 real(rp) :: exner_sfc
211 real(rp) :: rovcp
212 real(rp) :: cpovr
213 !-----------------------------------------------
214
215 rovcp = rdry / cpdry
216 cpovr = cpdry / rdry
217 exner_sfc = (pres_sfc / pres00)**rovcp
218
219 !$omp parallel do private(PT, exner)
220 !$acc parallel loop collapse(2) present(DENS_hyd, PRES_hyd, z, lcmesh3D, elem)
221 do ke=lcmesh3d%NeS, lcmesh3d%NeE
222 do p=1, elem%Np
223 ! d exner / dz = - g / ( Cp * PT0 ) * exp (- N2/g * z)
224 ! exner = exner(zs) - g^2 / (Cp * N^2) [ 1/PT (z) - 1/PT(zs) ]
225 pt = pottemp0 * exp( bruntvaisalafreq**2 / grav * z(p,ke) )
226 exner = exner_sfc + grav**2 / ( cpdry * bruntvaisalafreq**2 ) * ( 1.0_rp / pt - 1.0_rp / pottemp0 )
227
228 pres_hyd(p,ke) = pres00 * exner**cpovr
229 dens_hyd(p,ke) = pres_hyd(p,ke) / ( rdry * exner * pt )
230 end do
231 end do
232
233 return
234 end subroutine hydrostatic_calc_basicstate_constbvfreq
235
236 !> Calculate density and pressure in hydrostatic balance with a constant lapse rate of temperature
237 !!
238 !! @param DENS_hyd hydrostatic density [kg/m3]
239 !! @param PRES_hyd hydrostatic pressure [Pa]
240 !! @param TLAPS Constant lapse rate of atmosphere [K/m]
241 !! @param Temp0 Temperature at z=0 [K]
242 !! @param PRES_sfc Surface pressure [Pa]
243 !! @param x x-coordinate
244 !! @param y y-coordinate
245 !! @param z z-coordinate [m]
246 !! @param lcmesh3D A object to manage a local 3D mesh
247 !! @param elem A object to manage a 3D finite element
248!OCL SERIAL
249 subroutine hydrostatic_calc_basicstate_consttlaps( &
250 DENS_hyd, PRES_hyd, &
251 TLAPS, Temp0, PRES_sfc, x, y, z, lcmesh3D, elem )
252
253 implicit none
254
255 class(localmesh3d), intent(in) :: lcmesh3d
256 class(elementbase3d), intent(in) :: elem
257 real(rp), intent(out) :: dens_hyd(elem%np,lcmesh3d%nea)
258 real(rp), intent(out) :: pres_hyd(elem%np,lcmesh3d%nea)
259 real(rp), intent(in) :: x(elem%np,lcmesh3d%ne)
260 real(rp), intent(in) :: y(elem%np,lcmesh3d%ne)
261 real(rp), intent(in) :: z(elem%np,lcmesh3d%ne)
262 real(rp), intent(in) :: tlaps
263 real(rp), intent(in) :: temp0
264 real(rp), intent(in) :: pres_sfc
265
266 integer :: ke, p
267
268 real(rp) :: fac
269 real(rp) :: temp
270 !-----------------------------------------------
271
272 fac = grav / ( rdry * tlaps )
273
274 !$omp parallel do private(TEMP)
275 !$acc parallel loop collapse(2) present(DENS_hyd, PRES_hyd, z, lcmesh3D, elem)
276 do ke=lcmesh3d%NeS, lcmesh3d%NeE
277 do p=1, elem%Np
278 temp = temp0 * ( 1.0_rp - tlaps / temp0 * z(p,ke) )
279 pres_hyd(p,ke) = pres_sfc * ( temp / temp0 )**fac
280 dens_hyd(p,ke) = pres_hyd(p,ke) / ( rdry * temp )
281 end do
282 end do
283
284 return
285 end subroutine hydrostatic_calc_basicstate_consttlaps
286
287 !> Calculate density and pressure in hydrostatic balance with a constant lapse rate of potential temperature
288 !!
289 !! @param DENS_hyd hydrostatic density [kg/m3]
290 !! @param PRES_hyd hydrostatic pressure [Pa]
291 !! @param PTLAPS Constant lapse rate of potential temperature in atmosphere [K/m]
292 !! @param PotTemp0 Potentital temperature at z=0 [K]
293 !! @param PRES_sfc Surface pressure [Pa]
294 !! @param x x-coordinate
295 !! @param y y-coordinate
296 !! @param z z-coordinate [m]
297 !! @param lcmesh3D A object to manage a local 3D mesh
298 !! @param elem A object to manage a 3D finite element
299!OCL SERIAL
300 subroutine hydrostatic_calc_basicstate_constptlaps( &
301 DENS_hyd, PRES_hyd, &
302 PTLAPS, PotTemp0, PRES_sfc, x, y, z, lcmesh3D, elem )
303
304 implicit none
305
306 class(localmesh3d), intent(in) :: lcmesh3d
307 class(elementbase3d), intent(in) :: elem
308 real(rp), intent(out) :: dens_hyd(elem%np,lcmesh3d%nea)
309 real(rp), intent(out) :: pres_hyd(elem%np,lcmesh3d%nea)
310 real(rp), intent(in) :: x(elem%np,lcmesh3d%ne)
311 real(rp), intent(in) :: y(elem%np,lcmesh3d%ne)
312 real(rp), intent(in) :: z(elem%np,lcmesh3d%ne)
313 real(rp), intent(in) :: ptlaps
314 real(rp), intent(in) :: pottemp0
315 real(rp), intent(in) :: pres_sfc
316
317 integer :: ke, p
318
319 real(rp) :: pt
320 real(rp) :: exner
321 real(rp) :: exner_sfc
322 real(rp) :: rovcp
323 real(rp) :: cpovr
324 !-----------------------------------------------
325
326 rovcp = rdry / cpdry
327 cpovr = cpdry / rdry
328 exner_sfc = (pres_sfc / pres00)**rovcp
329
330 if ( ptlaps == 0.0_rp ) then
331 call hydrostatic_calc_basicstate_constpt( &
332 dens_hyd, pres_hyd, &
333 pottemp0, pres_sfc, x, y, z, lcmesh3d, elem )
334 else
335 !$omp parallel do private(PT, exner)
336 !$acc parallel loop collapse(2) present(DENS_hyd,PRES_hyd,z)
337 do ke=lcmesh3d%NeS, lcmesh3d%NeE
338 do p=1, elem%Np
339 ! d exner / dz = - g / ( Cp * PT0 ) / (1 + PTLAPS/PT0 * z)
340 ! exner = exner(zs) - g / (Cp * PTLAPS ) * log[ 1 + PTLAPS/PT0 * z ]
341 pt = pottemp0 + ptlaps * z(p,ke)
342 exner = exner_sfc - grav / ( cpdry * ptlaps ) * log( 1.0_rp + ptlaps / pottemp0 * z(p,ke) )
343
344 pres_hyd(p,ke) = pres00 * exner**cpovr
345 dens_hyd(p,ke) = pres_hyd(p,ke) / ( rdry * exner * pt )
346 end do
347 end do
348 end if
349
350 return
351 end subroutine hydrostatic_calc_basicstate_constptlaps
352
353 !> Build density in hydrostatic balance state of dry atmosphere
354 !!
355 !! @param DDENS Density deviation satisfing a discrete hydrostatic balance state of dry atmosphere [kg/m3]
356 !! @param DENS_hyd hydrostatic density [kg/m3]
357 !! @param PRES_hyd hydrostatic pressure [Pa]
358 !! @param POT Potentital temperature [K]
359 !! @param x x-coordinate
360 !! @param y y-coordinate
361 !! @param z model vertical coordinate
362 !! @param lcmesh A object to manage a local 3D mesh
363 !! @param elem A object to manage a 3D finite element
364 !! @param bnd_SFC_PRES Surface pressure which is given as surface boundary condition [Pa]
365!OCL SERIAL
366 subroutine hydrostatic_build_rho_xyz_dry( &
367 DDENS, &
368 DENS_hyd, PRES_hyd, &
369 POT, &
370 x, y, z, lcmesh, elem, &
371 bnd_SFC_PRES )
372
373 implicit none
374
375 class(localmesh3d), intent(in) :: lcmesh
376 class(elementbase3d), intent(in) :: elem
377 real(RP), intent(out) :: DDENS(elem%Np,lcmesh%NeA)
378 real(RP), intent(in) :: DENS_hyd(elem%Np,lcmesh%NeA)
379 real(RP), intent(in) :: PRES_hyd(elem%Np,lcmesh%NeA)
380 real(RP), intent(in) :: x(elem%Np,lcmesh%Ne)
381 real(RP), intent(in) :: y(elem%Np,lcmesh%Ne)
382 real(RP), intent(in) :: z(elem%Np,lcmesh%Ne)
383 real(RP), intent(in) :: POT(elem%Np,lcmesh%NeZ,lcmesh%NeX,lcmesh%NeY)
384 real(RP), intent(in), optional :: bnd_SFC_PRES(lcmesh%lcmesh2D%refElem2D%Np,lcmesh%Ne2DA)
385
386 real(RP) :: Rtot (elem%Np,lcmesh%NeZ,lcmesh%NeX,lcmesh%NeY)
387 real(RP) :: CPtot_ov_CVtot(elem%Np,lcmesh%NeZ,lcmesh%NeX,lcmesh%NeY)
388
389 !-----------------------------------------------
390
391 !$omp parallel
392 !$omp workshare
393 rtot(:,:,:,:) = rdry
394 cptot_ov_cvtot(:,:,:,:) = cpdry / cvdry
395 !$omp end workshare
396 !$omp end parallel
397
398 call hydrostatic_build_rho_xyz_moist( &
399 ddens, &
400 dens_hyd, pres_hyd, &
401 pot, rtot, cptot_ov_cvtot, &
402 x, y, z, lcmesh, elem, &
403 bnd_sfc_pres )
404
405 return
406 end subroutine hydrostatic_build_rho_xyz_dry
407
408 !> Build density in hydrostatic balance state of moist atmosphere
409 !!
410 !! @param DDENS Density deviation satisfing a discrete hydrostatic balance state of dry atmosphere [kg/m3]
411 !! @param DENS_hyd hydrostatic density [kg/m3]
412 !! @param PRES_hyd hydrostatic pressure [Pa]
413 !! @param POT Potentital temperature [K]
414 !! @param Rtot Specific gas constant of moist atmosphere [J/kg/K]
415 !! @param CPtot_ov_CVtot Specific heat ratio of moist atmosphere
416 !! @param x x-coordinate
417 !! @param y y-coordinate
418 !! @param z model vertical coordinate
419 !! @param lcmesh A object to manage a local 3D mesh
420 !! @param elem A object to manage a 3D finite element
421 !! @param bnd_SFC_PRES Surface pressure which is given as surface boundary condition [Pa]
422!OCL SERIAL
423 subroutine hydrostatic_build_rho_xyz_moist( &
424 DDENS, &
425 DENS_hyd, PRES_hyd, &
426 POT, Rtot, CPtot_ov_CVtot, &
427 x, y, z, lcmesh, elem, &
428 bnd_SFC_PRES )
429
430 use scale_const, only: &
431 eps0 => const_eps
432 use scale_element_operation_general, only: &
433 elementoperationgeneral_generate_vpordm1
434 implicit none
435
436 class(localmesh3d), intent(in) :: lcmesh
437 class(elementbase3d), intent(in) :: elem
438 real(RP), intent(out) :: DDENS(elem%Np,lcmesh%NeA)
439 real(RP), intent(in) :: DENS_hyd(elem%Np,lcmesh%NeA)
440 real(RP), intent(in) :: PRES_hyd(elem%Np,lcmesh%NeA)
441 real(RP), intent(in) :: x(elem%Np,lcmesh%Ne)
442 real(RP), intent(in) :: y(elem%Np,lcmesh%Ne)
443 real(RP), intent(in) :: z(elem%Np,lcmesh%Ne)
444 real(RP), intent(in) :: POT(elem%Np,lcmesh%NeZ,lcmesh%NeX,lcmesh%NeY)
445 real(RP), intent(in) :: Rtot(elem%Np,lcmesh%NeZ,lcmesh%NeX,lcmesh%NeY)
446 real(RP), intent(in) :: CPtot_ov_CVtot(elem%Np,lcmesh%NeZ,lcmesh%NeX,lcmesh%NeY)
447 real(RP), intent(in), optional :: bnd_SFC_PRES(lcmesh%lcmesh2D%refElem2D%Np,lcmesh%Ne2DA)
448
449 integer :: ke, ke2D
450 integer :: ke_x, ke_y, ke_z
451 integer :: itr_lin
452 integer :: itr_nlin
453
454 real(RP), parameter :: EPS = 1.0e-12_rp
455
456 type(sparsemat) :: Dz, Lift
457
458 type(gmres) :: gmres_hydro
459 real(RP), allocatable :: wj(:)
460 real(RP), allocatable :: pinv_v(:)
461 integer :: N, m
462 integer :: vmapM_z1D(elem%NfpTot,lcmesh%NeZ)
463 integer :: vmapP_z1D(elem%NfpTot,lcmesh%NeZ)
464 real(RP) :: VARS (elem%Np,lcmesh%NeZ)
465 real(RP) :: VARS0 (elem%Np,lcmesh%NeZ)
466 real(RP) :: VAR_DEL(elem%Np,lcmesh%NeZ)
467 real(RP) :: b(elem%Np,lcmesh%NeZ)
468 real(RP) :: Ax(elem%Np,lcmesh%NeZ)
469 real(RP) :: nz(elem%NfpTot,lcmesh%NeZ)
470 real(RP) :: DENS_hyd_z(elem%Np,lcmesh%NeZ)
471 real(RP) :: PRES_hyd_z(elem%Np,lcmesh%NeZ)
472 real(RP) :: PmatDlu(elem%Np,elem%Np,lcmesh%NeZ)
473 integer :: PmatDlu_ipiv(elem%Np,lcmesh%NeZ)
474 real(RP) :: PmatL(elem%Np,elem%Np,lcmesh%NeZ)
475 real(RP) :: PmatU(elem%Np,elem%Np,lcmesh%NeZ)
476 real(RP) :: GsqrtV_z(elem%Np,lcmesh%NeZ)
477
478 logical :: is_converged
479
480 real(RP) :: IntrpMat_VPOrdM1(elem%Np,elem%Np)
481
482 real(RP) :: bnd_SFC_PRES_tmp(lcmesh%lcmesh2D%refElem2D%Np,lcmesh%Ne2DA)
483 !-----------------------------------------------
484
485
486 if ( present(bnd_sfc_pres) ) then
487 !$omp parallel do
488 do ke2d=lcmesh%lcmesh2D%NeS, lcmesh%lcmesh2D%NeE
489 bnd_sfc_pres_tmp(:,ke2d) = bnd_sfc_pres(:,ke2d)
490 end do
491 else
492 !$omp parallel do
493 do ke2d=lcmesh%lcmesh2D%NeS, lcmesh%lcmesh2D%NeE
494 bnd_sfc_pres_tmp(:,ke2d) = pres_hyd(elem%Hslice(:,1),ke2d)
495 end do
496 end if
497
498 !--
499
500 n = elem%Np * lcmesh%NeZ
501 m = min(n / elem%Nnode_h1D**2, 30)
502! m = min(N / elem%Nnode_h1D**2, 256)
503
504 call gmres_hydro%Init( n, m, eps, eps0 )
505 allocate( wj(n), pinv_v(n) )
506
507 call dz%Init( elem%Dx3, storage_format='ELL' )
508 call lift%Init( elem%Lift, storage_format='ELL' )
509
510 call lcmesh%GetVmapZ1D( vmapm_z1d, vmapp_z1d ) ! (out)
511
512 call elementoperationgeneral_generate_vpordm1( intrpmat_vpordm1, & ! (out)
513 elem )
514
515 !-----------
516 do ke_y=1, lcmesh%NeY
517 do ke_x=1, lcmesh%NeX
518 ke2d = ke_x + (ke_y-1)*lcmesh%NeX
519 do ke_z=1, lcmesh%NeZ
520 ke = ke2d + (ke_z-1)*lcmesh%NeX*lcmesh%NeY
521 vars(:,ke_z) = 0.0_rp
522
523 vars0(:,ke_z) = vars(:,ke_z)
524 nz(:,ke_z) = lcmesh%normal_fn(:,ke,3)
525
526 dens_hyd_z(:,ke_z) = dens_hyd(:,ke)
527 pres_hyd_z(:,ke_z) = pres_hyd(:,ke)
528 gsqrtv_z(:,ke_z) = lcmesh%Gsqrt(:,ke) / lcmesh%GsqrtH(elem%IndexH2Dto3D(:),ke2d)
529 end do
530
531 do itr_nlin=1, 10
532
533 do ke_z=1, lcmesh%NeZ
534 var_del(:,ke_z) = 0.0_rp
535 end do
536
537 call eval_ax( ax(:,:), &
538 vars, vars0, pot(:,:,ke_x,ke_y), rtot(:,:,ke_x,ke_y), &
539 cptot_ov_cvtot(:,:,ke_x,ke_y), dens_hyd_z, pres_hyd_z, &
540 bnd_sfc_pres_tmp(:,ke2d), &
541 dz, lift, intrpmat_vpordm1, lcmesh, elem, &
542 nz, vmapm_z1d, vmapp_z1d, ke_x, ke_y )
543
544 do ke_z=1, lcmesh%NeZ
545 b(:,ke_z) = - ax(:,ke_z)
546 end do
547 if (lcmesh%tileID==1) then
548 if (itr_nlin > 1) then
549 log_progress(*) ke_x, ke_y, "itr_lin=", itr_lin
550 log_progress(*) "-------------------------------------"
551 end if
552 log_progress(*) ke_x, ke_y, "itr_nlin:", itr_nlin, 0, ": VAR", vars(elem%Colmask(:,1),1)
553 log_progress(*) ke_x, ke_y, "itr_nlin:", itr_nlin, 0, ": b", b(elem%Colmask(:,1),1)
554 if( io_l ) call flush(io_fid_log)
555 end if
556
557 if ( maxval(abs(b(:,:))) < 1.0e-10_rp ) exit
558
559 call construct_pmatinv( pmatdlu, pmatdlu_ipiv, pmatl, pmatu, & ! (out)
560 vars0, pot(:,:,ke_x,ke_y), rtot(:,:,ke_x,ke_y), & ! (in)
561 cptot_ov_cvtot(:,:,ke_x,ke_y), dens_hyd_z, pres_hyd_z, & ! (in)
562 dz, lift, intrpmat_vpordm1, gsqrtv_z, lcmesh, elem, & ! (in)
563 nz, vmapm_z1d, vmapp_z1d, ke_x, ke_y )
564
565 do itr_lin=1, 2*int(n/m)
566 !
567 call gmres_hydro_core( gmres_hydro, var_del, wj, is_converged, & ! (out)
568 vars, b, n, m, & ! (in)
569 pmatdlu, pmatdlu_ipiv, pmatl, pmatu, pinv_v, & ! (in)
570 pot(:,:,ke_x,ke_y), rtot(:,:,ke_x,ke_y), & ! (in)
571 cptot_ov_cvtot(:,:,ke_x,ke_y), dens_hyd_z, pres_hyd_z, & ! (in)
572 dz, lift, intrpmat_vpordm1, lcmesh, elem, & ! (in)
573 nz, vmapm_z1d, vmapp_z1d, ke_x, ke_y )
574
575 ! LOG_PROGRESS(*) ke_x, ke_y, "itr_lin:", itr_lin, ": VAR_DEL", VAR_DEL(elem%Colmask(:,1),1)
576 ! if( IO_L ) call flush(IO_FID_LOG)
577 if (is_converged) exit
578 end do ! itr_lin
579 do ke_z=1, lcmesh%NeZ
580 vars(:,ke_z) = vars(:,ke_z) + var_del(:,ke_z)
581 vars0(:,ke_z) = vars(:,ke_z)
582 end do
583 end do ! itr_nlin
584
585 do ke_z=1, lcmesh%NeZ
586 ke = ke_x + (ke_y-1)*lcmesh%NeX + (ke_z-1)*lcmesh%NeX*lcmesh%NeY
587 ddens(:,ke) = vars(:,ke_z)
588 end do
589 end do
590 end do
591
592 !
593 call gmres_hydro%Final()
594 call dz%Final()
595 call lift%Final()
596
597 return
598 end subroutine hydrostatic_build_rho_xyz_moist
599
600!-- private ---------------------------------
601
602!OCL SERIAL
603 subroutine gmres_hydro_core( gmres_hydro, x, wj, is_converged, &
604 x0, b, N, m, &
605 PmatDlu, PmatDlu_ipiv, PmatL, PmatU, pinv_v, & ! (in)
606 pot, rtot, cptot_ov_cvtot, dens_hyd, pres_hyd, & ! (in)
607 dz, lift, intrpmat_vpordm1, lmesh, elem, & ! (in)
608 nz, vmapm, vmapp, ke_x, ke_y )
609
610 implicit none
611
612 class(localmesh3d), intent(in) :: lmesh
613 class(elementbase3d), intent(in) :: elem
614 integer, intent(in) :: N
615 integer, intent(in) :: m
616
617 class(gmres), intent(inout) :: gmres_hydro
618 real(RP), intent(inout) :: x(N)
619 real(RP), intent(inout) :: wj(N)
620 logical, intent(out) :: is_converged
621 real(RP), intent(in) :: x0(N)
622 real(RP), intent(in) :: b(N)
623 real(RP), intent(in) :: PmatDlu(elem%Np,elem%Np,lmesh%NeZ)
624 integer, intent(in) :: PmatDlu_ipiv(elem%Np,lmesh%NeZ)
625 real(RP), intent(in) :: PmatL(elem%Np,elem%Np,lmesh%NeZ)
626 real(RP), intent(in) :: PmatU(elem%Np,elem%Np,lmesh%NeZ)
627 real(RP), intent(inout) :: pinv_v(N)
628 !---
629 real(RP), intent(in) :: POT(elem%Np,lmesh%NeZ)
630 real(RP), intent(in) :: Rtot(elem%Np,lmesh%NeZ)
631 real(RP), intent(in) :: CPtot_ov_CVtot(elem%Np,lmesh%NeZ)
632 real(RP), intent(in) :: DENS_hyd(elem%Np,lmesh%NeZ)
633 real(RP), intent(in) :: PRES_hyd(elem%Np,lmesh%NeZ)
634 class(sparsemat), intent(in) :: Dz, Lift
635 real(RP), intent(in) :: IntrpMat_VPOrdM1(elem%Np,elem%Np)
636 real(RP), intent(in) :: nz(elem%NfpTot,lmesh%NeZ)
637 integer, intent(in) :: vmapM(elem%NfpTot,lmesh%NeZ)
638 integer, intent(in) :: vmapP(elem%NfpTot,lmesh%NeZ)
639 integer, intent(in) :: ke_x, ke_y
640
641 integer :: j
642 !-------------------------------------------------------
643
644 call eval_ax_lin( wj(:), & ! (out)
645 x, x0, pot, rtot, cptot_ov_cvtot, & ! (in)
646 dens_hyd, pres_hyd, & ! (in)
647 dz, lift, intrpmat_vpordm1, lmesh, elem, & ! (in)
648 nz, vmapm, vmapp, ke_x, ke_y ) ! (in)
649
650 call gmres_hydro%Iterate_pre( b, wj, is_converged )
651 if (is_converged) return
652
653 do j=1, min(m, n)
654 call matmul_pinv_v( pinv_v, pmatdlu, pmatdlu_ipiv, pmatl, pmatu, gmres_hydro%v(:,j) )
655
656 call eval_ax_lin( wj(:), & ! (out)
657 pinv_v, x0, pot, rtot, cptot_ov_cvtot, & ! (in)
658 dens_hyd, pres_hyd, & ! (in)
659 dz, lift, intrpmat_vpordm1, lmesh, elem, & ! (in)
660 nz, vmapm, vmapp, ke_x, ke_y ) ! (in)
661
662 call gmres_hydro%Iterate_step_j( j, wj, is_converged )
663
664 log_info("GMRES check**:",*) "j=", j, "g:", gmres_hydro%g(j+1), "r:", gmres_hydro%r(j,j), "hj(j+1):", gmres_hydro%hj(j+1)
665 if ( is_converged ) exit
666 end do
667
668 do j=1, n
669 wj(j) = 0.0_rp
670 end do
671 call gmres_hydro%Iterate_post( wj )
672 call matmul_pinv_v_plus_x0( x, pmatdlu, pmatdlu_ipiv, pmatl, pmatu, wj)
673
674 return
675 contains
676
677!OCL SERIAL
678 subroutine matmul_pinv_v( pinv_v_, pDlu_, PmatDlu_ipiv_, pL, pU, v)
679 implicit none
680 real(RP), intent(out) :: pinv_v_(elem%Np,lmesh%NeZ)
681 real(RP), intent(in) :: pDlu_(elem%Np,elem%Np,lmesh%NeZ)
682 integer, intent(in) :: PmatDlu_ipiv_(elem%Np,lmesh%NeZ)
683 real(RP), intent(in) :: pL(elem%Np,elem%Np,lmesh%NeZ)
684 real(RP), intent(in) :: pU(elem%Np,elem%Np,lmesh%NeZ)
685 real(RP), intent(in) :: v(elem%Np,lmesh%NeZ)
686
687 integer :: k, n
688 real(RP) :: tmp(elem%Np)
689 integer :: vs, ve
690 integer :: info
691 !------------------------------------
692
693 n = elem%Np
694
695 vs = 1
696 ve = vs + elem%Np - 1
697 pinv_v_(:,1) = v(vs:ve,1)
698 call dgetrs('N', n, 1, pdlu_(:,:,1), n, pmatdlu_ipiv_(:,1), pinv_v_(:,1), n, info)
699
700 do k=2, lmesh%NeZ
701 vs = 1; ve = elem%Np
702 pinv_v_(:,k) = v(vs:ve,k) &
703 - matmul( pl(:,:,k), pinv_v_(:,k-1) )
704
705 call dgetrs('N', n, 1, pdlu_(:,:,k), n, pmatdlu_ipiv_(:,k), pinv_v_(:,k), n, info)
706 end do
707
708 !
709 do k=lmesh%NeZ-1, 1, -1
710 vs = 1; ve = elem%Np
711 tmp(vs:ve) = matmul( pu(:,:,k), pinv_v_(:,k+1) )
712 call dgetrs('N', n, 1, pdlu_(:,:,k), n, pmatdlu_ipiv_(:,k), tmp(:), n, info)
713
714 vs = 1
715 ve = vs + elem%Np - 1
716 pinv_v_(:,k) = pinv_v_(:,k) - tmp(vs:ve)
717 end do
718
719 return
720 end subroutine matmul_pinv_v
721
722!OCL SERIAL
723 subroutine matmul_pinv_v_plus_x0( x_, pDlu_, PmatDlu_ipiv_, pL, pU, v)
724 implicit none
725 real(RP), intent(inout) :: x_(elem%Np,lmesh%NeZ)
726 real(RP), intent(in) :: pDlu_(elem%Np,elem%Np,lmesh%NeZ)
727 integer, intent(in) :: PmatDlu_ipiv_(elem%Np,lmesh%NeZ)
728 real(RP), intent(in) :: pL(elem%Np,elem%Np,lmesh%NeZ)
729 real(RP), intent(in) :: pU(elem%Np,elem%Np,lmesh%NeZ)
730 real(RP), intent(in) :: v(elem%Np,lmesh%NeZ)
731
732 integer :: k
733 real(RP) :: tmp(elem%Np,lmesh%NeZ)
734
735 !------------------------------------
736
737 call matmul_pinv_v( tmp, pdlu_, pmatdlu_ipiv_, pl, pu, v)
738 !$omp parallel do
739 do k=1, lmesh%NeZ
740 x_(:,k) = x_(:,k) + tmp(:,k)
741 end do
742
743 return
744 end subroutine matmul_pinv_v_plus_x0
745
746 end subroutine gmres_hydro_core
747
748!OCL SERIAL
749 subroutine eval_ax( Ax, &
750 DDENS, DENS0, POT, Rtot, CPtot_ov_CVtot, & ! (in)
751 dens_hyd, pres_hyd, bnd_sfc_pres, & ! (in)
752 dz, lift, intrpmat_vpordm1, lmesh, elem, & ! (in)
753 nz, vmapm, vmapp, ke_x, ke_y ) ! (in)
754
755 implicit none
756 class(localmesh3d), intent(in) :: lmesh
757 class(elementbase3d), intent(in) :: elem
758
759 real(RP), intent(out) :: Ax(elem%Np,lmesh%NeZ)
760 real(RP), intent(in) :: DDENS (elem%Np,lmesh%NeZ)
761 real(RP), intent(in) :: DENS0(elem%Np,lmesh%NeZ)
762 real(RP), intent(in) :: POT(elem%Np,lmesh%NeZ)
763 real(RP), intent(in) :: Rtot(elem%Np,lmesh%NeZ)
764 real(RP), intent(in) :: CPtot_ov_CVtot(elem%Np,lmesh%NeZ)
765 real(RP), intent(in) :: DENS_hyd(elem%Np,lmesh%NeZ)
766 real(RP), intent(in) :: PRES_hyd(elem%Np,lmesh%NeZ)
767 real(RP), intent(in) :: bnd_SFC_PRES(lmesh%lcmesh2D%refElem2D%Np)
768 class(sparsemat), intent(in) :: Dz, Lift
769 real(RP), intent(in) :: IntrpMat_VPOrdM1(elem%Np,elem%Np)
770
771 real(RP), intent(in) :: nz(elem%NfpTot,lmesh%NeZ)
772 integer, intent(in) :: vmapM(elem%NfpTot,lmesh%NeZ)
773 integer, intent(in) :: vmapP(elem%NfpTot,lmesh%NeZ)
774 integer, intent(in) :: ke_x, ke_y
775
776 real(RP) :: DPRES(elem%Np), DENS(elem%Np)
777 real(RP) :: Fz(elem%Np), LiftDelFlx(elem%Np)
778 real(RP) :: del_flux(elem%NfpTot,lmesh%NeZ)
779
780 integer :: ke_z
781 integer :: ke, ke2D
782
783 real(RP) :: gamm
784 real(RP) :: RdOvP00
785 real(RP) :: rP0
786
787 real(RP) :: GsqrtV(elem%Np)
788 !-------------------------------------------
789
790 gamm = cpdry / cvdry
791 rdovp00 = rdry / pres00
792 rp0 = 1.0_rp / pres00
793
794 call cal_del_flux( del_flux, & ! (out)
795 ddens, pot, rtot, cptot_ov_cvtot, & ! (in)
796 dens_hyd, pres_hyd, bnd_sfc_pres, & ! (in)
797 nz, vmapm, vmapp, lmesh, elem ) ! (in)
798
799 !$omp parallel do private(ke, ke2D, DPRES, DENS, Fz, LiftDelFlx, GsqrtV)
800 do ke_z=1, lmesh%NeZ
801 ke = ke_x + (ke_y-1)*lmesh%NeX + (ke_z-1)*lmesh%NeX*lmesh%NeY
802 ke2d = lmesh%EMap3Dto2D(ke)
803
804 gsqrtv(:) = lmesh%Gsqrt(:,ke) / lmesh%GsqrtH(elem%IndexH2Dto3D,ke2d)
805
806 dens(:) = dens_hyd(:,ke_z) + ddens(:,ke_z)
807 dpres(:) = pres00 * ( rtot(:,ke_z) * rp0 * dens(:) * pot(:,ke_z) )**cptot_ov_cvtot(:,ke_z) !&
808 !- PRES_hyd(:,ke_z)
809 call sparsemat_matmul(dz, dpres, fz)
810 call sparsemat_matmul(lift, lmesh%Fscale(:,ke)*del_flux(:,ke_z), liftdelflx)
811
812 ! Ax(:,ke_z) = ( lmesh%Escale(:,ke,3,3) * Fz(:) + LiftDelFlx(:) ) / GsqrtV(:) &
813 ! + Grav * matmul(IntrpMat_VPOrdM1, DDENS(:,ke_z))
814 ax(:,ke_z) = ( lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) ) / gsqrtv(:) &
815 + grav * matmul(intrpmat_vpordm1, dens(:)) !DDENS(:,ke_z))
816 end do
817
818 return
819 end subroutine eval_ax
820
821!OCL SERIAL
822 subroutine cal_del_flux( del_flux, &
823 DDENS_, POT_, Rtot_, CPtot_ov_CVtot_, &
824 DENS_hyd, PRES_hyd, bnd_SFC_PRES, &
825 nz, vmapM, vmapP, lmesh, elem )
826
827 implicit none
828
829 class(localmesh3d), intent(in) :: lmesh
830 class(elementbase3d), intent(in) :: elem
831 real(RP), intent(out) :: del_flux(elem%NfpTot*lmesh%NeZ)
832 real(RP), intent(in) :: DDENS_(elem%Np*lmesh%NeZ)
833 real(RP), intent(in) :: POT_(elem%Np*lmesh%NeZ)
834 real(RP), intent(in) :: Rtot_(elem%Np*lmesh%NeZ)
835 real(RP), intent(in) :: CPtot_ov_CVtot_(elem%Np*lmesh%NeZ)
836 real(RP), intent(in) :: DENS_hyd(elem%Np*lmesh%NeZ)
837 real(RP), intent(in) :: PRES_hyd(elem%Np*lmesh%NeZ)
838 real(RP), intent(in) :: bnd_SFC_PRES(lmesh%lcmesh2D%refElem2D%Np)
839 real(RP), intent(in) :: nz(elem%NfpTot*lmesh%NeZ)
840 integer, intent(in) :: vmapM(elem%NfpTot*lmesh%NeZ)
841 integer, intent(in) :: vmapP(elem%NfpTot*lmesh%NeZ)
842
843 integer :: i, iP, iM
844 integer :: p, p2D, f, ke_z
845
846 real(RP) :: dpresP, dpresM
847 real(RP) :: RtotOvP00M, RtotOvP00P
848 real(RP) :: rP0
849 real(RP) :: fac
850
851 !-------------------------------
852
853 rp0 = 1.0_rp / pres00
854
855 !$omp parallel private( &
856 !$omp ke_z, p, p2D, f, i, iM, iP, fac, &
857 !$omp dpresM, dpresP, RtotOvP00M, RtotOvP00P )
858 !$omp do
859 do i=1, elem%NfpTot*lmesh%NeZ
860 del_flux(i) = 0.0_rp
861 end do
862 !$omp end do
863 !$omp do collapse(3)
864 do ke_z=1, lmesh%NeZ
865 do f=1, elem%Nfaces_v
866 do p2d=1, elem%Nfp_v
867 p = p2d + (f-1)*elem%Nfp_v + elem%Nfaces_h * elem%Nfp_h
868 i = p + (ke_z-1)*elem%NfpTot
869 im = vmapm(i); ip = vmapp(i)
870
871 rtotovp00m = rtot_(im) * rp0
872 rtotovp00p = rtot_(ip) * rp0
873
874 fac = 0.5_rp * ( 1.0_rp - sign(1.0_rp,nz(i)) )
875 dpresm = pres00 * ( rtotovp00m * (dens_hyd(im) + ddens_(im)) * pot_(im) )**cptot_ov_cvtot_(im) !- PRES_hyd(iM)
876 dpresp = pres00 * ( rtotovp00p * (dens_hyd(ip) + ddens_(ip)) * pot_(ip) )**cptot_ov_cvtot_(ip) !- PRES_hyd(iP)
877
878 if ( ke_z==1 .and. im==ip ) dpresp = bnd_sfc_pres(p2d)
879
880 del_flux(i) = fac * ( dpresp - dpresm ) * nz(i)
881 end do
882 end do
883 end do
884 !$omp end do
885 !$omp end parallel
886
887 return
888 end subroutine cal_del_flux
889
890!OCL SERIAL
891 subroutine eval_ax_lin( Ax, &
892 DDENS, DDENS0, POT, Rtot, CPtot_ov_CVtot, & ! (in)
893 dens_hyd, pres_hyd, & ! (in)
894 dz, lift, intrpmat_vpordm1, lmesh, elem, & ! (in)
895 nz, vmapm, vmapp, ke_x, ke_y )
896
897 implicit none
898 class(localmesh3d), intent(in) :: lmesh
899 class(elementbase3d), intent(in) :: elem
900
901 real(RP), intent(out) :: Ax(elem%Np,lmesh%NeZ)
902 real(RP), intent(in) :: DDENS (elem%Np,lmesh%NeZ)
903 real(RP), intent(in) :: DDENS0(elem%Np,lmesh%NeZ)
904 real(RP), intent(in) :: POT(elem%Np,lmesh%NeZ)
905 real(RP), intent(in) :: Rtot(elem%Np,lmesh%NeZ)
906 real(RP), intent(in) :: CPtot_ov_CVtot(elem%Np,lmesh%NeZ)
907 real(RP), intent(in) :: DENS_hyd(elem%Np,lmesh%NeZ)
908 real(RP), intent(in) :: PRES_hyd(elem%Np,lmesh%NeZ)
909 class(sparsemat), intent(in) :: Dz, Lift
910 real(RP), intent(in) :: IntrpMat_VPOrdM1(elem%Np,elem%Np)
911 real(RP), intent(in) :: nz(elem%NfpTot,lmesh%NeZ)
912 integer, intent(in) :: vmapM(elem%NfpTot,lmesh%NeZ)
913 integer, intent(in) :: vmapP(elem%NfpTot,lmesh%NeZ)
914 integer, intent(in) :: ke_x, ke_y
915
916 real(RP) :: PRES(elem%Np), PRES0(elem%Np)
917 real(RP) :: DPRES(elem%Np)
918 real(RP) :: Fz(elem%Np), LiftDelFlx(elem%Np)
919 real(RP) :: del_flux(elem%NfpTot,lmesh%NeZ)
920
921 integer :: ke_z
922 integer :: ke, ke2D
923
924 real(RP) :: gamm
925 real(RP) :: rP0
926
927 real(RP) :: GsqrtV(elem%Np)
928 !-------------------------------------------
929
930 gamm = cpdry/cvdry
931 rp0 = 1.0_rp / pres00
932
933 call cal_del_flux_lin( del_flux, & ! (out)
934 ddens, ddens0, pot, rtot, cptot_ov_cvtot, & ! (in)
935 dens_hyd, pres_hyd, & ! (in)
936 nz, vmapm, vmapp, lmesh, elem ) ! (in)
937
938 !$omp parallel do private(ke, ke2D, PRES, PRES0, DPRES, Fz, LiftDelFlx, GsqrtV)
939 do ke_z=1, lmesh%NeZ
940 ke = ke_x + (ke_y-1)*lmesh%NeX + (ke_z-1)*lmesh%NeX*lmesh%NeY
941 ke2d = lmesh%EMap3Dto2D(ke)
942
943 gsqrtv(:) = lmesh%Gsqrt(:,ke) / lmesh%GsqrtH(elem%IndexH2Dto3D,ke2d)
944
945 pres0(:) = pres00 * ( rtot(:,ke_z) * rp0 * (dens_hyd(:,ke_z) + ddens0(:,ke_z)) * pot(:,ke_z) )**cptot_ov_cvtot(:,ke_z)
946 dpres(:) = cptot_ov_cvtot(:,ke_z) * pres0(:) / (dens_hyd(:,ke_z) + ddens0(:,ke_z)) * ddens(:,ke_z)
947
948 call sparsemat_matmul(dz, dpres(:), fz)
949 call sparsemat_matmul(lift, lmesh%Fscale(:,ke)*del_flux(:,ke_z), liftdelflx)
950
951 ax(:,ke_z) = ( lmesh%Escale(:,ke,3,3) * fz(:) + liftdelflx(:) ) / gsqrtv(:) &
952 + grav * matmul(intrpmat_vpordm1, ddens(:,ke_z))
953 end do
954
955 return
956 end subroutine eval_ax_lin
957
958!OCL SERIAL
959 subroutine cal_del_flux_lin( del_flux, &
960 DDENS_, DDENS0_, POT_, Rtot_, CPtot_ov_CVtot_, &
961 DENS_hyd_, PRES_hyd_, &
962 nz, vmapM, vmapP, lmesh, elem )
963
964 implicit none
965
966 class(localmesh3d), intent(in) :: lmesh
967 class(elementbase3d), intent(in) :: elem
968 real(RP), intent(out) :: del_flux(elem%NfpTot*lmesh%NeZ)
969 real(RP), intent(in) :: DDENS_(elem%Np*lmesh%NeZ)
970 real(RP), intent(in) :: DDENS0_(elem%Np*lmesh%NeZ)
971 real(RP), intent(in) :: Rtot_(elem%Np*lmesh%NeZ)
972 real(RP), intent(in) :: CPtot_ov_CVtot_(elem%Np*lmesh%NeZ)
973 real(RP), intent(in) :: DENS_hyd_(elem%Np*lmesh%NeZ)
974 real(RP), intent(in) :: PRES_hyd_(elem%Np*lmesh%NeZ)
975 real(RP), intent(in) :: POT_(elem%Np*lmesh%NeZ)
976 real(RP), intent(in) :: nz(elem%NfpTot*lmesh%NeZ)
977 integer, intent(in) :: vmapM(elem%NfpTot*lmesh%NeZ)
978 integer, intent(in) :: vmapP(elem%NfpTot*lmesh%NeZ)
979
980 integer :: i, iP, iM
981 integer :: p, p2D, f, ke_z
982
983 real(RP) :: dpresP, dpresM
984 real(RP) :: pres0P, pres0M
985 real(RP) :: gamm
986 real(RP) :: rP0
987 real(RP) :: RtotOvP00M
988 real(RP) :: RtotOvP00P
989
990 real(RP) :: fac
991 !-------------------------------
992
993 gamm = cpdry/cvdry
994 rp0 = 1.0_rp / pres00
995
996 !$omp parallel private( &
997 !$omp p, p2D, f, ke_z, i, iM, iP, fac, &
998 !$omp dpresM, dpresP, pres0M, pres0P, RtotOvP00M, RtotOvP00P )
999 !$omp do
1000 do i=1, elem%NfpTot*lmesh%NeZ
1001 del_flux(i) = 0.0_rp
1002 end do
1003 !$omp end do
1004 !$omp do collapse(3)
1005 do ke_z=1, lmesh%NeZ
1006 do f=1, elem%Nfaces_v
1007 do p2d=1, elem%Nfp_v
1008 p = p2d + (f-1)*elem%Nfp_v + elem%Nfaces_h * elem%Nfp_h
1009 i = p + (ke_z-1)*elem%NfpTot
1010
1011 im = vmapm(i); ip = vmapp(i)
1012
1013 rtotovp00m = rtot_(im) * rp0
1014 rtotovp00p = rtot_(ip) * rp0
1015
1016 pres0m = pres00 * ( rtotovp00m * (dens_hyd_(im) + ddens0_(im)) * pot_(im) )**cptot_ov_cvtot_(im)
1017 pres0p = pres00 * ( rtotovp00p * (dens_hyd_(ip) + ddens0_(ip)) * pot_(ip) )**cptot_ov_cvtot_(ip)
1018
1019 fac = 0.5_rp * ( 1.0_rp - sign(1.0_rp,nz(i)) )
1020 dpresm = cptot_ov_cvtot_(im) * pres0m / (dens_hyd_(im) + ddens0_(im)) * ddens_(im)
1021 dpresp = cptot_ov_cvtot_(ip) * pres0p / (dens_hyd_(ip) + ddens0_(ip)) * ddens_(ip)
1022
1023 if ( ke_z == 1 .and. im == ip ) dpresp = 0.0_rp
1024
1025 del_flux(i) = fac * ( dpresp - dpresm ) * nz(i)
1026 end do
1027 end do
1028 end do
1029 !$omp end do
1030 !$omp end parallel
1031
1032 return
1033 end subroutine cal_del_flux_lin
1034
1035!OCL SERIAL
1036 subroutine construct_pmatinv( PmatDlu, PmatDlu_ipiv, PmatL, PmatU, & ! (out)
1037 ddens0, pot, rtot, cptot_ov_cvtot, dens_hyd, pres_hyd, & ! (in)
1038 dz, lift, intrpmat_vpordm1, gsqrtv, lmesh, elem, & ! (in)
1039 nz, vmapm, vmapp, ke_x, ke_y )
1040
1041 use scale_linalgebra, only: linalgebra_lu
1042 implicit none
1043
1044 class(localmesh3d), intent(in) :: lmesh
1045 class(elementbase3d), intent(in) :: elem
1046 real(RP), intent(out) :: PmatDlu(elem%Np,elem%Np,lmesh%NeZ)
1047 integer, intent(out) :: PmatDlu_ipiv(elem%Np,lmesh%NeZ)
1048 real(RP), intent(out) :: PmatL(elem%Np,elem%Np,lmesh%NeZ)
1049 real(RP), intent(out) :: PmatU(elem%Np,elem%Np,lmesh%NeZ)
1050 real(RP), intent(in) :: DDENS0(elem%Np,lmesh%NeZ)
1051 real(RP), intent(in) :: POT(elem%Np,lmesh%NeZ)
1052 real(RP), intent(in) :: Rtot(elem%Np,lmesh%NeZ)
1053 real(RP), intent(in) :: CPtot_ov_CVtot(elem%Np,lmesh%NeZ)
1054 real(RP), intent(in) :: DENS_hyd(elem%Np,lmesh%NeZ)
1055 real(RP), intent(in) :: PRES_hyd(elem%Np,lmesh%NeZ)
1056 class(sparsemat), intent(in) :: Dz, Lift
1057 real(RP), intent(in) :: IntrpMat_VPOrdM1(elem%Np,elem%Np)
1058 real(RP), intent(in) :: GsqrtV(elem%Np,lmesh%NeZ)
1059 real(RP), intent(in) :: nz(elem%NfpTot,lmesh%NeZ)
1060 integer, intent(in) :: vmapM(elem%NfpTot,lmesh%NeZ)
1061 integer, intent(in) :: vmapP(elem%NfpTot,lmesh%NeZ)
1062 integer, intent(in) :: ke_x, ke_y
1063
1064 real(RP) :: DENS0(elem%Np,lmesh%NeZ)
1065 real(RP) :: PRES0(elem%Np,lmesh%NeZ)
1066 integer :: ke_z, ke_z2
1067 integer :: ke, p, fp, v
1068 real(RP) :: gamm, rgamm, rP0
1069 real(RP) :: dz_p(elem%Np)
1070 real(RP) :: PmatD(elem%Np,elem%Np)
1071
1072 integer :: f1, f2, fp_s, fp_e
1073 integer :: FmV(elem%Nfp_v)
1074 integer :: FmV2 (elem%Nfp_v)
1075 real(RP) :: lift_op(elem%Np,elem%NfpTot)
1076 real(RP) :: lift_(elem%Np,elem%Np)
1077 real(RP) :: lift_2(elem%Np,elem%Np)
1078 real(RP) :: tmp(elem%Nfp_v)
1079
1080 real(RP) :: fac
1081 !--------------------------------------------------------
1082
1083 gamm = cpdry/cvdry
1084 rgamm = cvdry/cpdry
1085 rp0 = 1.0_rp / pres00
1086
1087 lift_op(:,:) = elem%Lift
1088
1089 !$omp parallel do
1090 do ke_z=1, lmesh%NeZ
1091 dens0(:,ke_z) = dens_hyd(:,ke_z) + ddens0(:,ke_z)
1092 pres0(:,ke_z) = pres00 * ( rtot(:,ke_z) * rp0 * dens0(:,ke_z) * pot(:,ke_z) )**cptot_ov_cvtot(:,ke_z)
1093 end do
1094
1095! !$omp parallel do private(ke, p, fp, v, f1, f2, ke_z2, dz_p, &
1096! !$omp tmp, lift_, lift_2, fac, &
1097! !$omp PmatD, FmV, FmV2, fp_s, fp_e)
1098 do ke_z=1, lmesh%NeZ
1099 ke = ke_x + (ke_y-1)*lmesh%NeX + (ke_z-1)*lmesh%NeX*lmesh%NeY
1100
1101 !-----
1102 pmatd(:,:) = 0.0_rp
1103 pmatl(:,:,ke_z) = 0.0_rp
1104 pmatu(:,:,ke_z) = 0.0_rp
1105
1106 do p=1, elem%Np
1107 dz_p(:) = lmesh%Escale(p,ke,3,3) * elem%Dx3(p,:)
1108 pmatd(p,:) = dz_p(:) * cptot_ov_cvtot(:,ke_z) * pres0(:,ke_z) / ( dens0(:,ke_z) * gsqrtv(p,ke_z) ) &
1109 + grav * intrpmat_vpordm1(p,:)
1110 end do
1111
1112 do f1=1, 2
1113 if (f1==1) then
1114 f2 = 2; ke_z2 = max(ke_z-1, 1)
1115 else
1116 f2 = 1; ke_z2 = min(ke_z+1, lmesh%NeZ)
1117 end if
1118 if ( (ke_z == 1 .and. f1==1) .or. (ke_z == lmesh%NeZ .and. f1==elem%Nfaces_v) ) then
1119 f2 = f1
1120 end if
1121
1122 fmv(:) = elem%Fmask_v(:,f1)
1123 fmv2(:) = elem%Fmask_v(:,f2)
1124
1125 fp_s = elem%Nfp_h * elem%Nfaces_h + 1 + (f1-1)*elem%Nfp_v
1126 fp_e = fp_s + elem%Nfp_v - 1
1127
1128 !
1129 lift_(:,:) = 0.0_rp
1130 lift_2(:,:) = 0.0_rp
1131 do fp=fp_s, fp_e
1132 p = fp-fp_s+1
1133
1134 fac = 0.5_rp * ( 1.0_rp - sign(1.0_rp,nz(fp,ke_z)) )
1135 tmp(:) = lift_op(fmv,fp) * lmesh%Fscale(fp,ke) * nz(fp,ke_z) * fac
1136 lift_(fmv,fmv(p)) = tmp(:) * cptot_ov_cvtot(fmv(p),ke_z ) * pres0(fmv(p),ke_z ) / dens0(fmv(p),ke_z )
1137 lift_2(fmv,fmv2(p)) = tmp(:) * cptot_ov_cvtot(fmv2(p),ke_z2) * pres0(fmv2(p),ke_z2) / dens0(fmv2(p),ke_z2)
1138 end do
1139
1140 !----
1141
1142 if ( ke_z == 1 .and. f1==1 ) then
1143 pmatd(:,:) = pmatd(:,:) - lift_(:,:)
1144 else if ( (ke_z == lmesh%NeZ .and. f1==elem%Nfaces_v) ) then
1145 else
1146 pmatd(:,:) = pmatd(:,:) - lift_(:,:)
1147 if (f1 == 1) then
1148 pmatl(:,:,ke_z) = lift_2(:,:)
1149 else
1150 pmatu(:,:,ke_z) = lift_2(:,:)
1151 end if
1152 end if
1153
1154 end do
1155 call get_pmatd_lu( pmatdlu(:,:,ke_z), pmatd(:,:), pmatdlu_ipiv(:,ke_z), elem%Np )
1156 end do
1157
1158 return
1159
1160 contains
1161!OCL SERIAL
1162 subroutine get_pmatd_lu( pmatDlu_, pmatD_, pmatDlu_ipiv_, N)
1163 use scale_linalgebra, only: linalgebra_lu
1164 implicit none
1165 integer, intent(in) :: N
1166 real(RP), intent(out) :: pmatDlu_(N,N)
1167 real(RP), intent(in) :: pmatD_(N,N)
1168 integer, intent(out) :: pmatDlu_ipiv_(N)
1169 integer :: info
1170 !------------------------------------------
1171
1172 pmatdlu_(:,:) = pmatd_(:,:)
1173 call linalgebra_lu(pmatdlu_, pmatdlu_ipiv_)
1174 return
1175 end subroutine get_pmatd_lu
1176
1177!OCL SERIAL
1178 subroutine get_pmatd_inv( pmatDinv_, pmatD_, pmatDlu_ipiv_, N)
1179 use scale_linalgebra, only: linalgebra_inv
1180 implicit none
1181 integer, intent(in) :: N
1182 real(RP), intent(out) :: pmatDinv_(N,N)
1183 real(RP), intent(in) :: pmatD_(N,N)
1184 integer, intent(out) :: pmatDlu_ipiv_(N)
1185 !------------------------------------------
1186
1187 pmatdinv_(:,:) = linalgebra_inv(pmatd_)
1188 return
1189 end subroutine get_pmatd_inv
1190 end subroutine construct_pmatinv
1191
1192end module scale_atm_dyn_dgm_hydrostatic
scale_atm_dyn_dgm_hydrostatic::hydrostatic_build_rho_xyz
Definition scale_atm_dyn_dgm_hydrostatic.F90:54
scale_element_operation_general::elementoperationgeneral_generate_vpordm1
Definition scale_element_operation_general.F90:69
scale_sparsemat::sparsemat_matmul
Definition scale_sparsemat.F90:57
scale_atm_dyn_dgm_hydrostatic
module FElib / Fluid dyn solver / Atmosphere / Common
Definition scale_atm_dyn_dgm_hydrostatic.F90:11
scale_atm_dyn_dgm_hydrostatic::hydrostatic_calc_basicstate_constptlaps
subroutine, public hydrostatic_calc_basicstate_constptlaps(dens_hyd, pres_hyd, ptlaps, pottemp0, pres_sfc, x, y, z, lcmesh3d, elem)
Calculate density and pressure in hydrostatic balance with a constant lapse rate of potential tempera...
Definition scale_atm_dyn_dgm_hydrostatic.F90:303
scale_atm_dyn_dgm_hydrostatic::hydrostatic_calc_basicstate_consttlaps
subroutine, public hydrostatic_calc_basicstate_consttlaps(dens_hyd, pres_hyd, tlaps, temp0, pres_sfc, x, y, z, lcmesh3d, elem)
Calculate density and pressure in hydrostatic balance with a constant lapse rate of temperature.
Definition scale_atm_dyn_dgm_hydrostatic.F90:252
scale_atm_dyn_dgm_hydrostatic::hydrostatic_calc_basicstate_constbvfreq
subroutine, public hydrostatic_calc_basicstate_constbvfreq(dens_hyd, pres_hyd, bruntvaisalafreq, pottemp0, pres_sfc, x, y, z, lcmesh3d, elem)
Calculate density and pressure in hydrostatic balance with a constant Brunt–Väisälä frequency.
Definition scale_atm_dyn_dgm_hydrostatic.F90:192
scale_atm_dyn_dgm_hydrostatic::hydrostatic_calc_basicstate_constt
subroutine, public hydrostatic_calc_basicstate_constt(dens_hyd, pres_hyd, temp0, pres_sfc, x, y, z, lcmesh3d, elem)
Calculate density and pressure in hydrostatic balance with a constant temperature.
Definition scale_atm_dyn_dgm_hydrostatic.F90:89
scale_atm_dyn_dgm_hydrostatic::hydrostatic_calc_basicstate_constpt
subroutine, public hydrostatic_calc_basicstate_constpt(dens_hyd, pres_hyd, pottemp0, pres_sfc, x, y, z, lcmesh3d, elem)
Calculate density and pressure in hydrostatic balance with a constant potential temperature.
Definition scale_atm_dyn_dgm_hydrostatic.F90:136
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_general
module FElib / Element / Operation with arbitary elements
Definition scale_element_operation_general.F90:12
scale_gmres
Module common / GMRES.
Definition scale_gmres.F90:15
scale_linalgebra
Module common / Linear algebra.
Definition scale_linalgebra.F90:11
scale_linalgebra::linalgebra_inv
real(rp) function, dimension(size(a, 1), size(a, 2)), public linalgebra_inv(a)
Calculate a inversion of matrix A.
Definition scale_linalgebra.F90:77
scale_linalgebra::linalgebra_lu
subroutine, public linalgebra_lu(a_lu, ipiv)
Perform LU factorization.
Definition scale_linalgebra.F90:174
scale_localmesh_3d
module FElib / Mesh / Local 3D
Definition scale_localmesh_3d.F90:10
scale_localmesh_base
module FElib / Mesh / Local, Base
Definition scale_localmesh_base.F90:10
scale_sparsemat
Module common / sparsemat.
Definition scale_sparsemat.F90:11
get_pmatd_lu
subroutine get_pmatd_lu(pmatdlu_, pmatd_, pmatdlu_ipiv_, n)
Definition scale_atm_dyn_dgm_hydrostatic.F90:1163
scale_element_base::elementbase3d
Derived type representing a 3D reference element.
Definition scale_element_base.F90:111
scale_element_base::elementbase
Derived type representing an arbitrary finite element.
Definition scale_element_base.F90:30
scale_element_hexahedral::hexahedralelement
Derived type representing a hexahedral element.
Definition scale_element_hexahedral.F90:31
scale_gmres::gmres
Derived type to provide a iterative solver for system of linear equations using GMRES.
Definition scale_gmres.F90:33
scale_localmesh_3d::localmesh3d
Derived type to manage a local 3D computational domain.
Definition scale_localmesh_3d.F90:32
scale_localmesh_base::localmeshbase
Derived type to manage a local computational domain (base type)
Definition scale_localmesh_base.F90:31
scale_sparsemat::sparsemat
Derived type to manage a sparse matrix.
Definition scale_sparsemat.F90:33
Generated by doxygen 1.14.0