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