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