347
353
354 implicit none
355
356 class(LocalMesh3D), intent(in) :: lmesh
357 class(ElementBase3D), intent(in) :: elem
358 class(LocalMesh2D), intent(in) :: lmesh2D
359 class(ElementBase2D), intent(in) :: elem2D
360 real(RP), intent(out) :: DENS_dt(elem%Np,lmesh%NeA)
361 real(RP), intent(out) :: MOMX_dt(elem%Np,lmesh%NeA)
362 real(RP), intent(out) :: MOMY_dt(elem%Np,lmesh%NeA)
363 real(RP), intent(out) :: MOMZ_dt(elem%Np,lmesh%NeA)
364 real(RP), intent(out) :: ETOT_dt(elem%Np,lmesh%NeA)
365 real(RP), intent(in) :: DDENS_(elem%Np,lmesh%NeA)
366 real(RP), intent(in) :: MOMX_(elem%Np,lmesh%NeA)
367 real(RP), intent(in) :: MOMY_(elem%Np,lmesh%NeA)
368 real(RP), intent(in) :: MOMZ_(elem%Np,lmesh%NeA)
369 real(RP), intent(in) :: ETOT_(elem%Np,lmesh%NeA)
370 real(RP), intent(in) :: DENS_hyd(elem%Np,lmesh%NeA)
371 real(RP), intent(in) :: PRES_hyd(elem%Np,lmesh%NeA)
372 real(RP), intent(in) :: DDENS0_(elem%Np,lmesh%NeA)
373 real(RP), intent(in) :: MOMX0_(elem%Np,lmesh%NeA)
374 real(RP), intent(in) :: MOMY0_(elem%Np,lmesh%NeA)
375 real(RP), intent(in) :: MOMZ0_(elem%Np,lmesh%NeA)
376 real(RP), intent(in) :: ETOT0_(elem%Np,lmesh%NeA)
377 real(RP), intent(in) :: Rtot(elem%Np,lmesh%NeA)
378 real(RP), intent(in) :: CVtot(elem%Np,lmesh%NeA)
379 real(RP), intent(in) :: CPtot(elem%Np,lmesh%NeA)
380 class(ElementOperationBase3D), intent(in) :: element3D_operation
381 class(SparseMat), intent(in) :: Dz, Lift
382 real(RP), intent(in) :: impl_fac
383 real(RP), intent(in) :: dt
384
385 real(RP) :: PROG_VARS (elem%Np,lmesh%NeZ,PRGVAR_NUM,lmesh%NeX*lmesh%NeY)
386 real(RP) :: DPRES (elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
387 real(RP) :: PROG_VARS0(elem%Np,lmesh%NeZ,PRGVAR_NUM,lmesh%NeX*lmesh%NeY)
388 real(RP) :: DPRES0 (elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
389 real(RP) :: b1D(3,elem%Nnode_v,lmesh%NeZ,elem%Nnode_h1D**2,lmesh%NeX*lmesh%NeY)
390 real(RP) :: GeoPot (elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
391 real(RP) :: KinHovDENS(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
392 integer :: ipiv(elem%Nnode_v*3*lmesh%NeZ,elem%Nnode_h1D**2)
393 real(RP) :: b1D_uv(elem%Nnode_v,lmesh%NeZ,2,elem%Nnode_h1D**2,lmesh%NeX*lmesh%NeY)
394 integer :: ipiv_uv(elem%Nnode_v*1*lmesh%NeZ,elem%Nnode_h1D**2)
395 real(RP) :: alph(elem%NfpTot,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
396 real(RP) :: Rtot_z(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
397 real(RP) :: CPtot_ov_CVtot(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
398 real(RP) :: DENS_hyd_z(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
399 real(RP) :: PRES_hyd_z(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
400 real(RP) :: GnnM_z(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
401 real(RP) :: G13_z(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
402 real(RP) :: G23_z(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
403 real(RP) :: GsqrtV_z(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
404 real(RP) :: nz(elem%NfpTot,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
405 integer :: vmapM(elem%NfpTot,lmesh%NeZ)
406 integer :: vmapP(elem%NfpTot,lmesh%NeZ)
407 integer :: ColMask(elem%Nnode_v)
408 integer :: ke_xy, ke_z, ke, ke2D, v
409 integer :: itr_nlin
410 integer :: kl, ku, nz_1D
411 integer :: kl_uv, ku_uv, nz_1D_uv
412 integer :: ij
413 logical :: is_converged
414
415 real(RP), allocatable :: PmatBnd(:,:,:)
416 real(RP), allocatable :: PmatBnd_uv(:,:,:)
417 integer :: info, info_uv
418
419 real(RP) :: DENS_(elem%Np)
420
421
422 call prof_rapstart( 'hevi_cal_vi_prep', 3)
423
424 nz_1d = elem%Nnode_v * 3 * lmesh%NeZ
425 kl = ( elem%Nnode_v + 1 ) * 3 - 1
426 ku = kl
427 nz_1d_uv = elem%Nnode_v * 1 * lmesh%NeZ
428 kl_uv = elem%Nnode_v
429 ku_uv = kl_uv
430 allocate( pmatbnd(2*kl+ku+1,nz_1d,elem%Nnode_h1D**2) )
431 allocate( pmatbnd_uv(2*kl_uv+ku_uv+1,nz_1d_uv,elem%Nnode_h1D**2) )
432
433 call lmesh%GetVmapZ1D( vmapm, vmapp )
434
435
436
437
438
439 do ke_xy=1, lmesh%NeX*lmesh%NeY
440 do ke_z=1, lmesh%NeZ
441 ke = ke_xy + (ke_z-1)*lmesh%NeX*lmesh%NeY
442 ke2d = lmesh%EMap3Dto2D(ke)
443
444 prog_vars(:,ke_z,dens_vid,ke_xy) = ddens0_(:,ke)
445 prog_vars(:,ke_z,momx_vid,ke_xy) = momx0_(:,ke)
446 prog_vars(:,ke_z,momy_vid,ke_xy) = momy0_(:,ke)
447 prog_vars(:,ke_z,momz_vid,ke_xy) = momz0_(:,ke)
448 prog_vars(:,ke_z,etot_vid,ke_xy) = etot0_(:,ke)
449
450 dens_hyd_z(:,ke_z,ke_xy) = dens_hyd(:,ke)
451 pres_hyd_z(:,ke_z,ke_xy) = pres_hyd(:,ke)
452 geopot(:,ke_z,ke_xy) = grav * lmesh%zlev(:,ke)
453
454 dens_(:) = dens_hyd(:,ke) + ddens0_(:,ke)
455 kinhovdens(:,ke_z,ke_xy) = 0.5_rp * ( &
456 momx0_(:,ke) * ( lmesh%G_ij(elem%IndexH2Dto3D,ke2d,1,1) * momx0_(:,ke) + lmesh%G_ij(elem%IndexH2Dto3D,ke2d,2,1) * momy0_(:,ke) ) &
457 + momy0_(:,ke) * ( lmesh%G_ij(elem%IndexH2Dto3D,ke2d,2,1) * momx0_(:,ke) + lmesh%G_ij(elem%IndexH2Dto3D,ke2d,2,2) * momy0_(:,ke) ) &
458 ) / dens_(:)**2
459
460 rtot_z(:,ke_z,ke_xy) = rtot(:,ke)
461 cptot_ov_cvtot(:,ke_z,ke_xy) = cptot(:,ke) / cvtot(:,ke)
462
463 dpres(:,ke_z,ke_xy) = &
464 ( cptot_ov_cvtot(:,ke_z,ke_xy) - 1.0_rp ) &
465 * ( etot0_(:,ke) - ( dens_(:) * ( kinhovdens(:,ke_z,ke_xy) + geopot(:,ke_z,ke_xy) ) + 0.5_rp * momz0_(:,ke)**2 / dens_(:) ) ) &
466 - pres_hyd(:,ke)
467
468 nz(:,ke_z,ke_xy) = lmesh%normal_fn(:,ke,3)
469 g13_z(:,ke_z,ke_xy) = lmesh%GI3(:,ke,1)
470 g23_z(:,ke_z,ke_xy) = lmesh%GI3(:,ke,2)
471 gsqrtv_z(:,ke_z,ke_xy) = lmesh%Gsqrt(:,ke) / lmesh%GsqrtH(elem%IndexH2Dto3D,ke2d)
472
473 gnnm_z(:,ke_z,ke_xy) = ( &
474 1.0_rp / gsqrtv_z(:,ke_z,ke_xy)**2 &
475 + g13_z(:,ke_z,ke_xy) * ( lmesh%GIJ(elem%IndexH2Dto3D,ke2d,1,1) * g13_z(:,ke_z,ke_xy) &
476 + lmesh%GIJ(elem%IndexH2Dto3D,ke2d,1,2) * g23_z(:,ke_z,ke_xy) ) &
477 + g23_z(:,ke_z,ke_xy) * ( lmesh%GIJ(elem%IndexH2Dto3D,ke2d,1,2) * g13_z(:,ke_z,ke_xy) &
478 + lmesh%GIJ(elem%IndexH2Dto3D,ke2d,2,2) * g23_z(:,ke_z,ke_xy) ) )
479 end do
480 end do
481
482
483 prog_vars0(:,:,:,:) = prog_vars(:,:,:,:)
484 dpres0(:,:,:) = dpres(:,:,:)
485
486
487
488 call prof_rapend( 'hevi_cal_vi_prep', 3)
489
490
491
492 if ( abs(impl_fac) > 0.0_rp ) then
493 call prof_rapstart( 'hevi_cal_vi_itr', 3)
494
495
496
497 do itr_nlin = 1, 1
498 call prof_rapstart( 'hevi_cal_vi_ax', 3)
499
500 call vi_eval_ax_uv( &
501 momx_dt(:,:), momy_dt(:,:), alph(:,:,:), &
502 prog_vars, dpres, prog_vars0, dpres0, &
503 ddens_, momx_, momy_, momz_, etot_, &
504 dens_hyd_z, pres_hyd_z, &
505 rtot_z, cptot_ov_cvtot, &
506 dz, lift, intrpmat_vpordm1, &
507 gnnm_z, g13_z, g23_z, gsqrtv_z, &
508 impl_fac, dt, &
509 lmesh, elem, nz, vmapm, vmapp, &
510 b1d_uv(:,:,:,:,:) )
511
512 call prof_rapend( 'hevi_cal_vi_ax', 3)
513
514 do ke_xy=1, lmesh%NeX * lmesh%NeY
515 call prof_rapstart( 'hevi_cal_vi_matbnd', 3)
516 call vi_construct_matbnd_uv( pmatbnd_uv(:,:,:), &
517 kl_uv, ku_uv, nz_1d_uv, &
518 prog_vars(:,:,:,ke_xy), kinhovdens(:,:,ke_xy), &
519 dens_hyd_z(:,:,ke_xy), pres_hyd_z(:,:,ke_xy), &
520 g13_z(:,:,ke_xy), g23_z(:,:,ke_xy), gsqrtv_z(:,:,ke_xy), &
521 alph(:,:,ke_xy), &
522 rtot_z(:,:,ke_xy), cptot_ov_cvtot(:,:,ke_xy), &
523 geopot(:,:,ke_xy), &
524 dz, lift, intrpmat_vpordm1, &
525 impl_fac, dt, &
526 lmesh, elem, nz(:,:,ke_xy), vmapm, vmapp, ke_xy, 1 )
527 call prof_rapend( 'hevi_cal_vi_matbnd', 3)
528
529 call prof_rapstart( 'hevi_cal_vi_lin', 3)
530
531
532 do ij=1, elem%Nnode_h1D**2
533 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)
534
535 colmask(:) = elem%Colmask(:,ij)
536 do ke_z=1, lmesh%NeZ
537 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)
538 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)
539 end do
540 end do
541
542
543 call prof_rapend( 'hevi_cal_vi_lin', 3)
544
545 end do
546
547 call prof_rapstart( 'hevi_cal_vi_ax', 3)
548 call vi_eval_ax( &
549 dens_dt(:,:), momz_dt(:,:), etot_dt(:,:), &
550 alph(:,:,:), &
551 prog_vars, dpres, prog_vars0, dpres0, &
552 ddens_, momx_, momy_, momz_, etot_, &
553 dens_hyd_z, pres_hyd_z, &
554 rtot_z, cptot_ov_cvtot, &
555 dz, lift, intrpmat_vpordm1, &
556 gnnm_z, g13_z, g23_z, gsqrtv_z, &
557 impl_fac, dt, &
558 lmesh, elem, nz, vmapm, vmapp, &
559 b1d(:,:,:,:,:) )
560 call prof_rapend( 'hevi_cal_vi_ax', 3)
561
562 do ke_xy=1, lmesh%NeX * lmesh%NeY
563 call prof_rapstart( 'hevi_cal_vi_matbnd', 3)
564 call vi_construct_matbnd( pmatbnd(:,:,:), &
565 kl, ku, nz_1d, &
566 prog_vars(:,:,:,ke_xy), kinhovdens(:,:,ke_xy), &
567 dens_hyd_z(:,:,ke_xy), pres_hyd_z(:,:,ke_xy), &
568 g13_z(:,:,ke_xy), g23_z(:,:,ke_xy), gsqrtv_z(:,:,ke_xy), &
569 alph(:,:,ke_xy), &
570 rtot_z(:,:,ke_xy), cptot_ov_cvtot(:,:,ke_xy), &
571 geopot(:,:,ke_xy), &
572 dz, lift, intrpmat_vpordm1, &
573 impl_fac, dt, &
574 lmesh, elem, nz(:,:,ke_xy), vmapm, vmapp, ke_xy, 1 )
575 call prof_rapend( 'hevi_cal_vi_matbnd', 3)
576
577 call prof_rapstart( 'hevi_cal_vi_lin', 3)
578
579
580 do ij=1, elem%Nnode_h1D**2
581 call dgbsv( nz_1d, kl, ku, 1, pmatbnd(:,:,ij), 2*kl+ku+1, ipiv(:,ij), b1d(:,:,:,ij,ke_xy), nz_1d, info)
582
583 colmask(:) = elem%Colmask(:,ij)
584 do ke_z=1, lmesh%NeZ
585 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)
586 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)
587 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)
588 end do
589 end do
590
591
592 do ke_z=1, lmesh%NeZ
593 dens_(:) = dens_hyd_z(:,ke_z,ke_xy) + prog_vars(:,ke_z,dens_vid,ke_xy)
594 dpres(:,ke_z,ke_xy) = &
595 ( cptot_ov_cvtot(:,ke_z,ke_xy) - 1.0_rp ) &
596 * ( prog_vars(:,ke_z,etot_vid,ke_xy) &
597 - ( dens_(:) * ( kinhovdens(:,ke_z,ke_xy) + geopot(:,ke_z,ke_xy) ) + 0.5_rp * prog_vars(:,ke_z,momz_vid,ke_xy)**2 / dens_(:) ) ) &
598 - pres_hyd_z(:,ke_z,ke_xy)
599 end do
600
601 call prof_rapend( 'hevi_cal_vi_lin', 3)
602
603 end do
604 end do
605
606 call prof_rapend( 'hevi_cal_vi_itr', 3)
607 end if
608
609 call prof_rapstart( 'hevi_cal_vi_retrun_var', 3)
610 if ( abs(impl_fac) > 0.0_rp) then
611
612 do ke_xy=1, lmesh%NeX * lmesh%NeY
613 do ke_z=1, lmesh%NeZ
614 ke = ke_xy + (ke_z-1)*lmesh%NeX*lmesh%NeY
615 dens_dt(:,ke) = ( prog_vars(:,ke_z,dens_vid,ke_xy) - ddens_(:,ke) ) / impl_fac
616 momx_dt(:,ke) = ( prog_vars(:,ke_z,momx_vid,ke_xy) - momx_(:,ke) ) / impl_fac
617 momy_dt(:,ke) = ( prog_vars(:,ke_z,momy_vid,ke_xy) - momy_(:,ke) ) / impl_fac
618 momz_dt(:,ke) = ( prog_vars(:,ke_z,momz_vid,ke_xy) - momz_(:,ke) ) / impl_fac
619 etot_dt(:,ke) = ( prog_vars(:,ke_z,etot_vid,ke_xy) - etot_(:,ke) ) / impl_fac
620 end do
621 end do
622 else
623 call vi_eval_ax_uv( &
624 momx_dt(:,:), momy_dt(:,:), alph(:,:,:), &
625 prog_vars, dpres, prog_vars0, dpres0, &
626 ddens_, momx_, momy_, momz_, etot_, &
627 dens_hyd_z, pres_hyd_z, &
628 rtot_z, cptot_ov_cvtot, &
629 dz, lift, intrpmat_vpordm1, &
630 gnnm_z, g13_z, g23_z, gsqrtv_z, &
631 impl_fac, dt, &
632 lmesh, elem, nz, vmapm, vmapp )
633
634 call vi_eval_ax( &
635 dens_dt(:,:), momz_dt(:,:), etot_dt(:,:), &
636 alph(:,:,:), &
637 prog_vars, dpres, prog_vars0, dpres0, &
638 ddens_, momx_, momy_, momz_, etot_, &
639 dens_hyd_z, pres_hyd_z, &
640 rtot_z, cptot_ov_cvtot, &
641 dz, lift, intrpmat_vpordm1, &
642 gnnm_z, g13_z, g23_z, gsqrtv_z, &
643 impl_fac, dt, &
644 lmesh, elem, nz, vmapm, vmapp )
645 end if
646 call prof_rapend( 'hevi_cal_vi_retrun_var', 3)
647
648
649 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)