297
303
304 implicit none
305
306 class(LocalMesh3D), intent(in) :: lmesh
307 class(ElementBase3D), intent(in) :: elem
308 class(LocalMesh2D), intent(in) :: lmesh2D
309 class(ElementBase2D), intent(in) :: elem2D
310 real(RP), intent(out) :: DENS_dt(elem%Np,lmesh%NeA)
311 real(RP), intent(out) :: MOMX_dt(elem%Np,lmesh%NeA)
312 real(RP), intent(out) :: MOMY_dt(elem%Np,lmesh%NeA)
313 real(RP), intent(out) :: MOMZ_dt(elem%Np,lmesh%NeA)
314 real(RP), intent(out) :: ETOT_dt(elem%Np,lmesh%NeA)
315 real(RP), intent(in) :: DDENS_(elem%Np,lmesh%NeA)
316 real(RP), intent(in) :: MOMX_(elem%Np,lmesh%NeA)
317 real(RP), intent(in) :: MOMY_(elem%Np,lmesh%NeA)
318 real(RP), intent(in) :: MOMZ_(elem%Np,lmesh%NeA)
319 real(RP), intent(in) :: ETOT_(elem%Np,lmesh%NeA)
320 real(RP), intent(in) :: DENS_hyd(elem%Np,lmesh%NeA)
321 real(RP), intent(in) :: PRES_hyd(elem%Np,lmesh%NeA)
322 real(RP), intent(in) :: DDENS0_(elem%Np,lmesh%NeA)
323 real(RP), intent(in) :: MOMX0_(elem%Np,lmesh%NeA)
324 real(RP), intent(in) :: MOMY0_(elem%Np,lmesh%NeA)
325 real(RP), intent(in) :: MOMZ0_(elem%Np,lmesh%NeA)
326 real(RP), intent(in) :: ETOT0_(elem%Np,lmesh%NeA)
327 real(RP), intent(in) :: Rtot(elem%Np,lmesh%NeA)
328 real(RP), intent(in) :: CVtot(elem%Np,lmesh%NeA)
329 real(RP), intent(in) :: CPtot(elem%Np,lmesh%NeA)
330 class(ElementOperationBase3D), intent(in) :: element3D_operation
331 class(SparseMat), intent(in) :: Dz, Lift
332 real(RP), intent(in) :: impl_fac
333 real(RP), intent(in) :: dt
334
335 real(RP) :: PROG_VARS (elem%Np,lmesh%NeZ,PRGVAR_NUM,lmesh%NeX*lmesh%NeY)
336 real(RP) :: DPRES (elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
337 real(RP) :: PROG_VARS0(elem%Np,lmesh%NeZ,PRGVAR_NUM,lmesh%NeX*lmesh%NeY)
338 real(RP) :: DPRES0 (elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
339 real(RP) :: b1D(3,elem%Nnode_v,lmesh%NeZ,elem%Nnode_h1D**2,lmesh%NeX*lmesh%NeY)
340 real(RP) :: GeoPot (elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
341 real(RP) :: KinHovDENS(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
342 integer :: ipiv(elem%Nnode_v*3*lmesh%NeZ,elem%Nnode_h1D**2)
343 real(RP) :: b1D_uv(elem%Nnode_v,lmesh%NeZ,2,elem%Nnode_h1D**2,lmesh%NeX*lmesh%NeY)
344 integer :: ipiv_uv(elem%Nnode_v*1*lmesh%NeZ,elem%Nnode_h1D**2)
345 real(RP) :: alph(elem%NfpTot,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
346 real(RP) :: Rtot_z(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
347 real(RP) :: CPtot_ov_CVtot(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
348 real(RP) :: DENS_hyd_z(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
349 real(RP) :: PRES_hyd_z(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
350 real(RP) :: GnnM_z(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
351 real(RP) :: G13_z(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
352 real(RP) :: G23_z(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
353 real(RP) :: GsqrtV_z(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
354 real(RP) :: nz(elem%NfpTot,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
355 integer :: vmapM(elem%NfpTot,lmesh%NeZ)
356 integer :: vmapP(elem%NfpTot,lmesh%NeZ)
357 integer :: ColMask(elem%Nnode_v)
358 integer :: ke_xy, ke_z, ke, ke2D, v
359 integer :: itr_nlin
360 integer :: kl, ku, nz_1D
361 integer :: kl_uv, ku_uv, nz_1D_uv
362 integer :: ij, info
363 logical :: is_converged
364
365 real(RP), allocatable :: PmatBnd(:,:,:)
366 real(RP), allocatable :: PmatBnd_uv(:,:,:)
367
368 real(RP) :: DENS_(elem%Np)
369
370
371 call prof_rapstart( 'hevi_cal_vi_prep', 3)
372
373 nz_1d = elem%Nnode_v * 3 * lmesh%NeZ
374 kl = ( elem%Nnode_v + 1 ) * 3 - 1
375 ku = kl
376 nz_1d_uv = elem%Nnode_v * 1 * lmesh%NeZ
377 kl_uv = elem%Nnode_v
378 ku_uv = kl_uv
379 allocate( pmatbnd(2*kl+ku+1,nz_1d,elem%Nnode_h1D**2) )
380 allocate( pmatbnd_uv(2*kl_uv+ku_uv+1,nz_1d_uv,elem%Nnode_h1D**2) )
381
382 call lmesh%GetVmapZ1D( vmapm, vmapp )
383
384
385
386
387
388 do ke_xy=1, lmesh%NeX*lmesh%NeY
389 do ke_z=1, lmesh%NeZ
390 ke = ke_xy + (ke_z-1)*lmesh%NeX*lmesh%NeY
391 ke2d = lmesh%EMap3Dto2D(ke)
392
393 prog_vars(:,ke_z,dens_vid,ke_xy) = ddens0_(:,ke)
394 prog_vars(:,ke_z,momx_vid,ke_xy) = momx0_(:,ke)
395 prog_vars(:,ke_z,momy_vid,ke_xy) = momy0_(:,ke)
396 prog_vars(:,ke_z,momz_vid,ke_xy) = momz0_(:,ke)
397 prog_vars(:,ke_z,etot_vid,ke_xy) = etot0_(:,ke)
398
399 dens_hyd_z(:,ke_z,ke_xy) = dens_hyd(:,ke)
400 pres_hyd_z(:,ke_z,ke_xy) = pres_hyd(:,ke)
401 geopot(:,ke_z,ke_xy) = grav * lmesh%zlev(:,ke)
402
403 dens_(:) = dens_hyd(:,ke) + ddens0_(:,ke)
404 kinhovdens(:,ke_z,ke_xy) = 0.5_rp * ( &
405 momx0_(:,ke) * ( lmesh%G_ij(elem%IndexH2Dto3D,ke2d,1,1) * momx0_(:,ke) + lmesh%G_ij(elem%IndexH2Dto3D,ke2d,2,1) * momy0_(:,ke) ) &
406 + momy0_(:,ke) * ( lmesh%G_ij(elem%IndexH2Dto3D,ke2d,2,1) * momx0_(:,ke) + lmesh%G_ij(elem%IndexH2Dto3D,ke2d,2,2) * momy0_(:,ke) ) &
407 ) / dens_(:)**2
408
409
410 rtot_z(:,ke_z,ke_xy) = rtot(:,ke)
411 cptot_ov_cvtot(:,ke_z,ke_xy) = cptot(:,ke) / cvtot(:,ke)
412
413 dpres(:,ke_z,ke_xy) = &
414 ( cptot_ov_cvtot(:,ke_z,ke_xy) - 1.0_rp ) &
415 * ( etot0_(:,ke) - ( dens_(:) * ( kinhovdens(:,ke_z,ke_xy) + geopot(:,ke_z,ke_xy) ) + 0.5_rp * momz0_(:,ke)**2 / dens_(:) ) ) &
416 - pres_hyd(:,ke)
417
418 nz(:,ke_z,ke_xy) = lmesh%normal_fn(:,ke,3)
419 g13_z(:,ke_z,ke_xy) = lmesh%GI3(:,ke,1)
420 g23_z(:,ke_z,ke_xy) = lmesh%GI3(:,ke,2)
421 gsqrtv_z(:,ke_z,ke_xy) = lmesh%Gsqrt(:,ke) / lmesh%GsqrtH(elem%IndexH2Dto3D,ke2d)
422
423 gnnm_z(:,ke_z,ke_xy) = ( 1.0_rp / gsqrtv_z(:,ke_z,ke_xy)**2 &
424 + g13_z(:,ke_z,ke_xy)**2 + g23_z(:,ke_z,ke_xy) )
425 end do
426 end do
427
428
429 prog_vars0(:,:,:,:) = prog_vars(:,:,:,:)
430 dpres0(:,:,:) = dpres(:,:,:)
431
432
433
434 call prof_rapend( 'hevi_cal_vi_prep', 3)
435
436
437
438 if ( abs(impl_fac) > 0.0_rp ) then
439 call prof_rapstart( 'hevi_cal_vi_itr', 3)
440
441
442
443 do itr_nlin = 1, 1
444 call prof_rapstart( 'hevi_cal_vi_ax', 3)
445
446 call vi_eval_ax_uv( &
447 momx_dt(:,:), momy_dt(:,:), alph(:,:,:), &
448 prog_vars, dpres, prog_vars0, dpres0, &
449 ddens_, momx_, momy_, momz_, etot_, &
450 dens_hyd_z, pres_hyd_z, &
451 rtot_z, cptot_ov_cvtot, &
452 dz, lift, intrpmat_vpordm1, &
453 gnnm_z, g13_z, g23_z, gsqrtv_z, &
454 impl_fac, dt, &
455 lmesh, elem, nz, vmapm, vmapp, &
456 b1d_uv(:,:,:,:,:) )
457
458 call prof_rapend( 'hevi_cal_vi_ax', 3)
459
460 do ke_xy=1, lmesh%NeX * lmesh%NeY
461 call prof_rapstart( 'hevi_cal_vi_matbnd', 3)
462 call vi_construct_matbnd_uv( pmatbnd_uv(:,:,:), &
463 kl_uv, ku_uv, nz_1d_uv, &
464 prog_vars(:,:,:,ke_xy), kinhovdens(:,:,ke_xy), &
465 dens_hyd_z(:,:,ke_xy), pres_hyd_z(:,:,ke_xy), &
466 g13_z(:,:,ke_xy), g23_z(:,:,ke_xy), gsqrtv_z(:,:,ke_xy), &
467 alph(:,:,ke_xy), &
468 rtot_z(:,:,ke_xy), cptot_ov_cvtot(:,:,ke_xy), &
469 geopot(:,:,ke_xy), &
470 dz, lift, intrpmat_vpordm1, &
471 impl_fac, dt, &
472 lmesh, elem, nz(:,:,ke_xy), vmapm, vmapp, ke_xy, 1 )
473 call prof_rapend( 'hevi_cal_vi_matbnd', 3)
474
475 call prof_rapstart( 'hevi_cal_vi_lin', 3)
476
477
478 do ij=1, elem%Nnode_h1D**2
479 call dgbsv( nz_1d_uv, kl_uv, ku_uv, 2, pmatbnd_uv(:,:,ij), 2*kl_uv+ku_uv+1, ipiv_uv(:,ij), b1d_uv(:,:,:,ij,ke_xy), nz_1d_uv, info)
480
481 colmask(:) = elem%Colmask(:,ij)
482 do ke_z=1, lmesh%NeZ
483 prog_vars(colmask(:),ke_z,momx_vid,ke_xy) = prog_vars(colmask(:),ke_z,momx_vid,ke_xy) + b1d_uv(:,ke_z,1,ij,ke_xy)
484 prog_vars(colmask(:),ke_z,momy_vid,ke_xy) = prog_vars(colmask(:),ke_z,momy_vid,ke_xy) + b1d_uv(:,ke_z,2,ij,ke_xy)
485 end do
486 end do
487
488
489 call prof_rapend( 'hevi_cal_vi_lin', 3)
490
491 end do
492
493 call prof_rapstart( 'hevi_cal_vi_ax', 3)
494 call vi_eval_ax( &
495 dens_dt(:,:), momz_dt(:,:), etot_dt(:,:), &
496 alph(:,:,:), &
497 prog_vars, dpres, prog_vars0, dpres0, &
498 ddens_, momx_, momy_, momz_, etot_, &
499 dens_hyd_z, pres_hyd_z, &
500 rtot_z, cptot_ov_cvtot, &
501 dz, lift, intrpmat_vpordm1, &
502 gnnm_z, g13_z, g23_z, gsqrtv_z, &
503 impl_fac, dt, &
504 lmesh, elem, nz, vmapm, vmapp, &
505 b1d(:,:,:,:,:) )
506 call prof_rapend( 'hevi_cal_vi_ax', 3)
507
508 do ke_xy=1, lmesh%NeX * lmesh%NeY
509 call prof_rapstart( 'hevi_cal_vi_matbnd', 3)
510 call vi_construct_matbnd( pmatbnd(:,:,:), &
511 kl, ku, nz_1d, &
512 prog_vars(:,:,:,ke_xy), kinhovdens(:,:,ke_xy), &
513 dens_hyd_z(:,:,ke_xy), pres_hyd_z(:,:,ke_xy), &
514 g13_z(:,:,ke_xy), g23_z(:,:,ke_xy), gsqrtv_z(:,:,ke_xy), &
515 alph(:,:,ke_xy), &
516 rtot_z(:,:,ke_xy), cptot_ov_cvtot(:,:,ke_xy), &
517 geopot(:,:,ke_xy), &
518 dz, lift, intrpmat_vpordm1, &
519 impl_fac, dt, &
520 lmesh, elem, nz(:,:,ke_xy), vmapm, vmapp, ke_xy, 1 )
521 call prof_rapend( 'hevi_cal_vi_matbnd', 3)
522
523 call prof_rapstart( 'hevi_cal_vi_lin', 3)
524
525
526 do ij=1, elem%Nnode_h1D**2
527 call dgbsv( nz_1d, kl, ku, 1, pmatbnd(:,:,ij), 2*kl+ku+1, ipiv(:,ij), b1d(:,:,:,ij,ke_xy), nz_1d, info)
528
529 colmask(:) = elem%Colmask(:,ij)
530 do ke_z=1, lmesh%NeZ
531 prog_vars(colmask(:),ke_z,dens_vid,ke_xy) = prog_vars(colmask(:),ke_z,dens_vid,ke_xy) + b1d(1,:,ke_z,ij,ke_xy)
532 prog_vars(colmask(:),ke_z,momz_vid,ke_xy) = prog_vars(colmask(:),ke_z,momz_vid,ke_xy) + b1d(2,:,ke_z,ij,ke_xy)
533 prog_vars(colmask(:),ke_z,etot_vid,ke_xy) = prog_vars(colmask(:),ke_z,etot_vid,ke_xy) + b1d(3,:,ke_z,ij,ke_xy)
534 end do
535 end do
536
537
538 do ke_z=1, lmesh%NeZ
539 dens_(:) = dens_hyd_z(:,ke_z,ke_xy) + prog_vars(:,ke_z,dens_vid,ke_xy)
540 dpres(:,ke_z,ke_xy) = &
541 ( cptot_ov_cvtot(:,ke_z,ke_xy) - 1.0_rp ) &
542 * ( prog_vars(:,ke_z,etot_vid,ke_xy) &
543 - ( dens_(:) * ( kinhovdens(:,ke_z,ke_xy) + geopot(:,ke_z,ke_xy) ) + 0.5_rp * prog_vars(:,ke_z,momz_vid,ke_xy)**2 / dens_(:) ) ) &
544 - pres_hyd_z(:,ke_z,ke_xy)
545 end do
546
547 call prof_rapend( 'hevi_cal_vi_lin', 3)
548
549 end do
550 end do
551
552 call prof_rapend( 'hevi_cal_vi_itr', 3)
553 end if
554
555 call prof_rapstart( 'hevi_cal_vi_retrun_var', 3)
556 if ( abs(impl_fac) > 0.0_rp) then
557
558 do ke_xy=1, lmesh%NeX * lmesh%NeY
559 do ke_z=1, lmesh%NeZ
560 ke = ke_xy + (ke_z-1)*lmesh%NeX*lmesh%NeY
561 dens_dt(:,ke) = ( prog_vars(:,ke_z,dens_vid,ke_xy) - ddens_(:,ke) ) / impl_fac
562 momx_dt(:,ke) = ( prog_vars(:,ke_z,momx_vid,ke_xy) - momx_(:,ke) ) / impl_fac
563 momy_dt(:,ke) = ( prog_vars(:,ke_z,momy_vid,ke_xy) - momy_(:,ke) ) / impl_fac
564 momz_dt(:,ke) = ( prog_vars(:,ke_z,momz_vid,ke_xy) - momz_(:,ke) ) / impl_fac
565 etot_dt(:,ke) = ( prog_vars(:,ke_z,etot_vid,ke_xy) - etot_(:,ke) ) / impl_fac
566 end do
567 end do
568 else
569 call vi_eval_ax_uv( &
570 momx_dt(:,:), momy_dt(:,:), alph(:,:,:), &
571 prog_vars, dpres, prog_vars0, dpres0, &
572 ddens_, momx_, momy_, momz_, etot_, &
573 dens_hyd_z, pres_hyd_z, &
574 rtot_z, cptot_ov_cvtot, &
575 dz, lift, intrpmat_vpordm1, &
576 gnnm_z, g13_z, g23_z, gsqrtv_z, &
577 impl_fac, dt, &
578 lmesh, elem, nz, vmapm, vmapp )
579
580 call vi_eval_ax( &
581 dens_dt(:,:), momz_dt(:,:), etot_dt(:,:), &
582 alph(:,:,:), &
583 prog_vars, dpres, prog_vars0, dpres0, &
584 ddens_, momx_, momy_, momz_, etot_, &
585 dens_hyd_z, pres_hyd_z, &
586 rtot_z, cptot_ov_cvtot, &
587 dz, lift, intrpmat_vpordm1, &
588 gnnm_z, g13_z, g23_z, gsqrtv_z, &
589 impl_fac, dt, &
590 lmesh, elem, nz, vmapm, vmapp )
591 end if
592 call prof_rapend( 'hevi_cal_vi_retrun_var', 3)
593
594 return
module FElib / Fluid dyn solver / Atmosphere / Nonhydrostatic model / HEVI / Common
subroutine, public atm_dyn_dgm_nonhydro3d_etot_hevi_common_construct_matbnd(pmatbnd, kl, ku, nz_1d, prog_vars0, kinhovdens00, dens_hyd, pres_hyd, g13, g23, gsqrtv, alph, rtot, cptot_ov_cvtot, geopot, dz, lift, intrpmat_vpordm1, impl_fac, dt, lmesh, elem, nz, vmapm, vmapp, ke_x, ke_y)
subroutine, public atm_dyn_dgm_nonhydro3d_etot_hevi_common_eval_ax(dens_t, momz_t, etot_t, alph, prog_vars, dpres, prog_vars0, dpres0, ddens00, momx00, momy00, momz00, entot00, dens_hyd, pres_hyd, rtot, cptot_ov_cvtot, dz, lift, intrpmat_vpordm1, gnnm, g13, g23, gsqrtv, impl_fac, dt, lmesh, elem, nz, vmapm, vmapp, b1d_ij)
subroutine, public atm_dyn_dgm_nonhydro3d_etot_hevi_common_eval_ax_uv(momx_t, momy_t, alph, prog_vars, dpres, prog_vars0, dpres0, ddens00, momx00, momy00, momz00, entot00, dens_hyd, pres_hyd, rtot, cptot_ov_cvtot, dz, lift, intrpmat_vpordm1, gnnm, g13, g23, gsqrtv, impl_fac, dt, lmesh, elem, nz, vmapm, vmapp, b1d_ij_uv)
subroutine, public atm_dyn_dgm_nonhydro3d_etot_hevi_common_construct_matbnd_uv(pmatbnd_uv, kl_uv, ku_uv, nz_1d_uv, prog_vars0, kinhovdens00, dens_hyd, pres_hyd, g13, g23, gsqrtv, alph, rtot, cptot_ov_cvtot, geopot, dz, lift, intrpmat_vpordm1, impl_fac, dt, lmesh, elem, nz, vmapm, vmapp, ke_x, ke_y)