diff --git a/FElib/src/Makefile b/FElib/src/Makefile index 58101d55..eab10590 100644 --- a/FElib/src/Makefile +++ b/FElib/src/Makefile @@ -23,6 +23,7 @@ VPATH = \ fluid_dyn_solver: \ surface: \ turbulence: \ + bl_turbulence: \ microphysics: \ radiation: \ model_framework: @@ -171,6 +172,10 @@ OBJS_NAME_TURBULENCE = \ scale_atm_phy_tb_dgm_globalsmg.o \ scale_atm_phy_tb_dgm_common.o +OBJS_NAME_BL_TURBULENCE = \ + scale_atm_phy_bl_dgm_mynn_lv2.o \ + scale_atm_phy_bl_dgm_common.o + OBJS_NAME_MICROPHYS = \ scale_atm_phy_mp_dgm_common.o \ scale_atm_phy_mp_lscond.o @@ -197,6 +202,7 @@ OBJS_NAME = \ $(OBJS_NAME_FLUID_DYN_SOLVER) \ $(OBJS_NAME_SURFACE) \ $(OBJS_NAME_TURBULENCE) \ + $(OBJS_NAME_BL_TURBULENCE) \ $(OBJS_NAME_MICROPHYS) \ $(OBJS_NAME_RADIATION) \ $(OBJS_NAME_MODEL_FRAMEWORK) diff --git a/FElib/src/bl_turbulence/scale_atm_phy_bl_dgm_common.F90 b/FElib/src/bl_turbulence/scale_atm_phy_bl_dgm_common.F90 new file mode 100644 index 00000000..63ddfc04 --- /dev/null +++ b/FElib/src/bl_turbulence/scale_atm_phy_bl_dgm_common.F90 @@ -0,0 +1,998 @@ +!> module FElib / Atmosphere / Physics / boundary layer turbulence +!! +!! @par Description +!! Boundary layer turbulence process +!! +!! @author Yuta Kawai, Team SCALE +!< +!------------------------------------------------------------------------------- +#include "scaleFElib.h" +module scale_atm_phy_bl_dgm_common + !----------------------------------------------------------------------------- + ! + !++ Used modules + ! + use scale_precision + use scale_io + use scale_prc + use scale_prof + use scale_const, only: & + EPS => CONST_EPS + use scale_sparsemat + use scale_element_base, only: & + ElementBase1D, ElementBase2D, ElementBase3D + use scale_element_hexahedral, only: HexahedralElement + use scale_localmesh_2d, only: LocalMesh2D + use scale_localmesh_3d, only: LocalMesh3D + use scale_mesh_base3d, only: MeshBase3D + use scale_localmeshfield_base, only: LocalMeshField3D + use scale_meshfield_base, only: MeshField3D + + use scale_element_operation_base, only: ElementOperationBase3D + + !----------------------------------------------------------------------------- + implicit none + private + !----------------------------------------------------------------------------- + ! + !++ Public procedures + ! + public :: atm_phy_bl_dgm_common_calc_tendency + + !----------------------------------------------------------------------------- + ! + !++ Private procedure + ! + + !----------------------------------------------------------------------------- + ! + !++ Private parameters & variables + ! +contains +!OCL SERIAL + subroutine atm_phy_bl_dgm_common_calc_tendency( & + RHOU_tp, RHOV_tp, DRHOT_tp, & ! (out) + DDENS_, MOMX_, MOMY_, DRHOT_, & ! (in) + PT_, DENS_hyd, PRES_hyd, NU, KH, & ! (in) + element3D_operation, C_IP, dtsec, & ! (in) + lmesh, elem, elem1D, is_bound, & ! (in) + use_delta_form ) ! (in) + use scale_atm_dyn_dgm_hevi_common_linalgebra, only: & + atm_dyn_dgm_hevi_common_linalgebra_get_param + implicit none + class(LocalMesh3D), intent(in), target :: lmesh + class(ElementBase3D), intent(in) :: elem + class(ElementBase1D), intent(in) :: elem1D + real(RP), intent(out) :: RHOU_tp(elem%Np,lmesh%NeA) + real(RP), intent(out) :: RHOV_tp(elem%Np,lmesh%NeA) + real(RP), intent(out) :: DRHOT_tp(elem%Np,lmesh%NeA) + real(RP), intent(in) :: DDENS_(elem%Np,lmesh%NeA) + real(RP), intent(in) :: MOMX_(elem%Np,lmesh%NeA) + real(RP), intent(in) :: MOMY_(elem%Np,lmesh%NeA) + real(RP), intent(in) :: DRHOT_(elem%Np,lmesh%NeA) + real(RP), intent(in) :: PT_(elem%Np,lmesh%NeA) + real(RP), intent(in) :: DENS_hyd(elem%Np,lmesh%NeA) + real(RP), intent(in) :: PRES_hyd(elem%Np,lmesh%NeA) + real(RP), intent(in) :: NU(elem%Np,lmesh%NeA) + real(RP), intent(in) :: KH(elem%Np,lmesh%NeA) + class(ElementOperationBase3D), intent(in) :: element3D_operation + real(RP), intent(in) :: C_IP + real(RP), intent(in) :: dtsec + logical, intent(in) :: is_bound(elem%NfpTot,lmesh%Ne) + logical, intent(in) :: use_delta_form + + class(LocalMesh2D), pointer :: lmesh2D + class(ElementBase2D), pointer :: elem2D + + real(RP) :: PROG_VARS (elem%Np,lmesh%NeX*lmesh%NeY,lmesh%NeZ,3) + real(RP) :: alph_M(elem%NfpTot,lmesh%Ne) + real(RP) :: alph_H(elem%NfpTot,lmesh%Ne) + real(RP) :: GsqrtV(elem%Np,lmesh%Ne) + + integer :: vmapM(elem%NfpTot,lmesh%Ne) + integer :: vmapP(elem%NfpTot,lmesh%Ne) + integer :: ke_xy, ke_z, ke, ke2d + integer :: p + + integer :: im, jm + + real(RP) :: DENS(elem%Np,lmesh%Ne) + real(RP), allocatable :: b1D_ij(:,:,:,:,:) + real(RP), allocatable :: BndMatL(:,:,:,:,:) + real(RP), allocatable :: BndMatD(:,:,:,:,:) + real(RP), allocatable :: G(:,:,:,:,:,:) + + real(RP) :: impl_fac + !------------------------------------------------------------------------ + + lmesh2D => lmesh%lcmesh2D + elem2D => lmesh2D%refElem2D + impl_fac = 1.0_RP * dtsec + + call lmesh%GetVmapZ3D( vmapM, vmapP ) ! (out) + call atm_dyn_dgm_hevi_common_linalgebra_get_param( im, jm, & ! (out) + elem%Nnode_v ) + + allocate( b1D_ij(im*elem%Nnode_v,3,jm,lmesh%Ne2D,lmesh%NeZ) ) + allocate( BndMatD(im*elem%Nnode_v,elem%Nnode_v,jm,lmesh%Ne2D,2) ) + allocate( BndMatL(im*elem%Nnode_v,elem%Nnode_v,jm,lmesh%Ne2D,2) ) + allocate( G(im*elem%Nnode_v,elem%Nnode_v,jm,lmesh%Ne2D,lmesh%NeZ,2) ) + + !$omp parallel do collapse(2) private( ke, ke2D ) + do ke_z =1, lmesh%NeZ + do ke_xy=1, lmesh%NeX * lmesh%NeY + ke = ke_xy + (ke_z-1)*lmesh%Ne2D + ke2D = lmesh%EMap3Dto2D(ke) + + DENS(:,ke) = DENS_hyd(:,ke) + DDENS_(:,ke) + + PROG_VARS(:,ke_xy,ke_z,1) = MOMX_ (:,ke) + PROG_VARS(:,ke_xy,ke_z,2) = MOMY_ (:,ke) + PROG_VARS(:,ke_xy,ke_z,3) = DENS(:,ke) * PT_(:,ke) + + do p=1, elem%Np + GsqrtV(p,ke) = lmesh%Gsqrt(p,ke) / lmesh%GsqrtH(elem%IndexH2Dto3D(p),ke2D) + end do + end do + end do + + call eval_Ax( RHOU_tp, RHOV_tp, DRHOT_tp, alph_M, alph_H, & !(out) + PROG_VARS, MOMX_, MOMY_, PT_, NU, KH, DENS, GsqrtV, & !(in) + impl_fac, dtsec, lmesh, elem, vmapM, vmapP, is_bound, & !(in) + element3D_operation, C_IP, im, jm, b1D_ij, use_delta_form ) !(in) + + call vi_solve( PROG_VARS, & ! (inout) + BndMatL, BndMatD, G, b1D_ij, & ! (inout) + DENS, NU, KH, GsqrtV, C_IP, dtsec, impl_fac, & ! (in) + im, jm, lmesh, elem, elem1D, use_delta_form ) ! (in) + + !--- + !$omp parallel do collapse(2) private( ke ) + do ke_z =1, lmesh%NeZ + do ke_xy=1, lmesh%NeX * lmesh%NeY + ke = ke_xy + (ke_z-1)*lmesh%Ne2D + RHOU_tp (:,ke) = ( PROG_VARS(:,ke_xy,ke_z,1) - MOMX_ (:,ke) ) / impl_fac + RHOV_tp (:,ke) = ( PROG_VARS(:,ke_xy,ke_z,2) - MOMY_ (:,ke) ) / impl_fac + DRHOT_tp(:,ke) = ( PROG_VARS(:,ke_xy,ke_z,3) - DENS(:,ke) * PT_(:,ke) ) / impl_fac + end do + end do + + return + end subroutine atm_phy_bl_dgm_common_calc_tendency + +!- Private subroutines ----------------------------- + +!OCL SERIAL + subroutine vi_solve( PROG_VARS, & ! (inout) + BndMatL, BndMatD, G, b1D_ij, & ! (inout) + DENS, NU, KH, GsqrtV, C_IP, dtsec, impl_fac, & ! (in) + im, jm, lmesh, elem, elem1D, use_delta_form ) ! (in) + implicit none + integer, intent(in) :: im, jm + class(LocalMesh3D), intent(in) :: lmesh + class(ElementBase3D), intent(in) :: elem + class(ElementBase1D), intent(in) :: elem1D + real(RP), intent(inout) :: PROG_VARS(elem%Nnode_h1D**2,elem%Nnode_v,lmesh%Ne2D,lmesh%NeZ,3) + real(RP), intent(inout) :: BndMatL(im*elem%Nnode_v,elem%Nnode_v,jm,lmesh%Ne2D,2) + real(RP), intent(inout) :: BndMatD(im*elem%Nnode_v,elem%Nnode_v,jm,lmesh%Ne2D,2) + real(RP), intent(inout) :: G(im*elem%Nnode_v,elem%Nnode_v,jm,lmesh%Ne2D,lmesh%NeZ,2) + real(RP), intent(inout) :: b1D_ij(im*elem%Nnode_v,3,jm,lmesh%Ne2D,lmesh%NeZ) + real(RP), intent(in) :: DENS(elem%Np,lmesh%Ne) + real(RP), intent(in) :: NU(elem%Np,lmesh%NeA) + real(RP), intent(in) :: KH(elem%Np,lmesh%NeA) + real(RP), intent(in) :: GsqrtV(elem%Np,lmesh%Ne) + real(RP), intent(in) :: C_IP + real(RP), intent(in) :: dtsec + real(RP), intent(in) :: impl_fac + logical, intent(in) :: use_delta_form + + integer :: ke_z, ke_xy + integer :: i, j, ij + integer :: pv, pv1, pv2, pp, p2 + + real(RP) :: tmp(im*elem%Nnode_v,2) + real(RP) :: tmp_b(im*elem%Nnode_v,3) + !------------------------------------------------------- + + do ke_z=1, lmesh%NeZ + call construct_matbnd_sip( & + BndMatL(:,:,:,:,1), BndMatD(:,:,:,:,1), G(:,:,:,:,ke_z,1), & ! (out) + DENS, NU, GsqrtV, elem1D%Dx1, elem1D%M, elem1D%invM, & ! (in) + impl_fac, dtsec, C_IP, lmesh, elem, im, jm, ke_z ) ! (in) + + call construct_matbnd_sip( & + BndMatL(:,:,:,:,2), BndMatD(:,:,:,:,2), G(:,:,:,:,ke_z,2), & ! (out) + DENS, KH, GsqrtV, elem1D%Dx1, elem1D%M, elem1D%invM, & ! (in) + impl_fac, dtsec, C_IP, lmesh, elem, im, jm, ke_z ) ! (in) + + if ( ke_z > 1 ) then + !$omp parallel do collapse(2) private( pp, p2, tmp, tmp_b ) + do ke_xy=1, lmesh%NeX * lmesh%NeY + do j=1, jm + !* D_k <- D_k - L_k * G_{k-1} ------------------ + do pv2=1, elem%Nnode_v + tmp(:,:) = 0.0_RP + do pv=1, elem%Nnode_v + do pv1=1, elem%Nnode_v + do i=1, im + pp = i + (pv1-1)*im; p2 = i + (pv -1)*im + tmp(pp,1) = tmp(pp,1) + BndMatL(pp,pv,j,ke_xy,1) * G(p2,pv2,j,ke_xy,ke_z-1,1) + tmp(pp,2) = tmp(pp,2) + BndMatL(pp,pv,j,ke_xy,2) * G(p2,pv2,j,ke_xy,ke_z-1,2) + end do + end do + end do + BndMatD(:,pv2,j,ke_xy,1) = BndMatD(:,pv2,j,ke_xy,1) - tmp(:,1) + BndMatD(:,pv2,j,ke_xy,2) = BndMatD(:,pv2,j,ke_xy,2) - tmp(:,2) + end do ! loop for pv2 + + !* b_k <- b_k - L_k * b_{k-1} ------------------ + tmp_b(:,:) = 0.0_RP + do pv=1, elem%Nnode_v + do pv1=1, elem%Nnode_v + do i=1, im + pp = i + (pv1-1)*im; p2 = i + (pv -1)*im + tmp_b(pp,1) = tmp_b(pp,1) + BndMatL(pp,pv,j,ke_xy,1) * b1D_ij(p2,1,j,ke_xy,ke_z-1) + tmp_b(pp,2) = tmp_b(pp,2) + BndMatL(pp,pv,j,ke_xy,1) * b1D_ij(p2,2,j,ke_xy,ke_z-1) + tmp_b(pp,3) = tmp_b(pp,3) + BndMatL(pp,pv,j,ke_xy,2) * b1D_ij(p2,3,j,ke_xy,ke_z-1) + end do + end do + end do + b1D_ij(:,1,j,ke_xy,ke_z) = b1D_ij(:,1,j,ke_xy,ke_z) - tmp_b(:,1) + b1D_ij(:,2,j,ke_xy,ke_z) = b1D_ij(:,2,j,ke_xy,ke_z) - tmp_b(:,2) + b1D_ij(:,3,j,ke_xy,ke_z) = b1D_ij(:,3,j,ke_xy,ke_z) - tmp_b(:,3) + end do + end do + end if + + call linkernel_solve_sip( b1D_ij(:,:,:,:,ke_z), G(:,:,:,:,ke_z,1), G(:,:,:,:,ke_z,2),& + BndMatD(:,:,:,:,1), BndMatD(:,:,:,:,2), elem%Nnode_v, im, jm, lmesh%Ne2D, ke_z == lmesh%NeZ ) + + end do + do ke_z=lmesh%NeZ-1, 1, -1 + !$omp parallel do collapse(2) private( pp, p2, tmp_b ) + do ke_xy=1, lmesh%NeX * lmesh%NeY + do j=1, jm + tmp_b(:,:) = 0.0_RP + do pv2=1, elem%Nnode_v + do pv =1, elem%Nnode_v + do i=1, im + pp = i + (pv -1)*im; p2 = i + (pv2-1)*im + tmp_b(pp,1) = tmp_b(pp,1) + G(pp,pv2,j,ke_xy,ke_z,1) * b1D_ij(p2,1,j,ke_xy,ke_z+1) + tmp_b(pp,2) = tmp_b(pp,2) + G(pp,pv2,j,ke_xy,ke_z,1) * b1D_ij(p2,2,j,ke_xy,ke_z+1) + tmp_b(pp,3) = tmp_b(pp,3) + G(pp,pv2,j,ke_xy,ke_z,2) * b1D_ij(p2,3,j,ke_xy,ke_z+1) + end do + end do + end do + b1D_ij(:,1,j,ke_xy,ke_z) = b1D_ij(:,1,j,ke_xy,ke_z) - tmp_b(:,1) + b1D_ij(:,2,j,ke_xy,ke_z) = b1D_ij(:,2,j,ke_xy,ke_z) - tmp_b(:,2) + b1D_ij(:,3,j,ke_xy,ke_z) = b1D_ij(:,3,j,ke_xy,ke_z) - tmp_b(:,3) + end do + end do + end do + + !$omp parallel do collapse(2) private( p2, pp ) + do ke_z=1, lmesh%NeZ + do ke_xy=1, lmesh%NeX * lmesh%NeY + do j=1, jm + if ( use_delta_form ) then + do pv=1, elem%Nnode_v + do i=1, im + p2 = i + (pv-1)*im; pp = i + (j-1)*im + PROG_VARS(pp,pv,ke_xy,ke_z,1) = PROG_VARS(pp,pv,ke_xy,ke_z,1) + b1D_ij(p2,1,j,ke_xy,ke_z) + PROG_VARS(pp,pv,ke_xy,ke_z,2) = PROG_VARS(pp,pv,ke_xy,ke_z,2) + b1D_ij(p2,2,j,ke_xy,ke_z) + PROG_VARS(pp,pv,ke_xy,ke_z,3) = PROG_VARS(pp,pv,ke_xy,ke_z,3) + b1D_ij(p2,3,j,ke_xy,ke_z) + end do + end do + else + do pv=1, elem%Nnode_v + do i=1, im + p2 = i + (pv-1)*im; pp = i + (j-1)*im + PROG_VARS(pp,pv,ke_xy,ke_z,1) = b1D_ij(p2,1,j,ke_xy,ke_z) + PROG_VARS(pp,pv,ke_xy,ke_z,2) = b1D_ij(p2,2,j,ke_xy,ke_z) + PROG_VARS(pp,pv,ke_xy,ke_z,3) = b1D_ij(p2,3,j,ke_xy,ke_z) + end do + end do + end if + end do + end do + end do + + return + end subroutine vi_solve + +!OCL SERIAL + subroutine linkernel_solve_sip( b, G1, G2, & ! (inout) + BndMatD1, BndMatD2, nv, im, jm, Ne2D, top_flag ) ! (in) + implicit none + + integer, intent(in) :: nv + integer, intent(in) :: im + integer, intent(in) :: jm + integer, intent(in) :: Ne2D + real(RP), intent(inout) :: b(im*nv,3,jm,Ne2D) + real(RP), intent(inout) :: G1(im*nv,nv,jm,Ne2D) + real(RP), intent(inout) :: G2(im*nv,nv,jm,Ne2D) + real(RP), intent(in) :: BndMatD1(im*nv,nv,jm,Ne2D) + real(RP), intent(in) :: BndMatD2(im*nv,nv,jm,Ne2D) + logical, intent(in) :: top_flag + + integer :: ke_xy + integer :: i, j + integer :: pv1, pv2 + integer :: iv + integer :: pp + + integer :: info + integer :: nrhs1, nrhs2 + integer :: ipiv1(nv), ipiv2(nv) + + real(RP) :: Amat1(nv,nv) + real(RP) :: RHS1(nv,nv+2) + real(RP) :: Amat2(nv,nv) + real(RP) :: RHS2(nv,nv+2) + !------------------------------------------------------------ + + !$omp parallel do collapse(2) & + !$omp private(ke_xy,j,i,pv1,pv2,iv,pp,info, & + !$omp nrhs1,nrhs2,ipiv1,ipiv2,Amat1,RHS1,Amat2,RHS2 ) + do ke_xy=1, Ne2D + do j=1, jm + do i=1, im + + do pv2=1, nv + do pv1=1, nv + pp = i + (pv1-1)*im + Amat1(pv1,pv2) = BndMatD1(pp,pv2,j,ke_xy) + Amat2(pv1,pv2) = BndMatD2(pp,pv2,j,ke_xy) + end do + end do + + RHS1(:,:) = 0.0_RP + RHS2(:,:) = 0.0_RP + + if ( top_flag ) then + nrhs1 = 2 + do iv=1, 2 + do pv1=1, nv + pp = i + (pv1-1)*im + RHS1(pv1,iv) = b(pp,iv,j,ke_xy) + end do + end do + + nrhs2 = 1 + do pv1=1, nv + pp = i + (pv1-1)*im + RHS2(pv1,1) = b(pp,3,j,ke_xy) + end do + else + nrhs1 = nv + 2 + nrhs2 = nv + 1 + + do pv2=1, nv + do pv1=1, nv + pp = i + (pv1-1)*im + RHS1(pv1,pv2) = G1(pp,pv2,j,ke_xy) + RHS2(pv1,pv2) = G2(pp,pv2,j,ke_xy) + end do + end do + do iv=1, 2 + do pv1=1, nv + pp = i + (pv1-1)*im + RHS1(pv1,nv+iv) = b(pp,iv,j,ke_xy) + end do + end do + do pv1=1, nv + pp = i + (pv1-1)*im + RHS2(pv1,nv+1) = b(pp,3,j,ke_xy) + end do + + end if + + ! LU factorization: + call DGETRF( nv, nv, Amat1, nv, ipiv1, info ) + if ( info /= 0 ) then + LOG_ERROR('linkernel_solve_sip',*) 'NU, DGETRF failed: info=', info, ', i=', i, ', j=', j, ', ke_xy=', ke_xy + call PRC_abort + end if + call DGETRF( nv, nv, Amat2, nv, ipiv2, info ) + if ( info /= 0 ) then + LOG_ERROR('linkernel_solve_sip',*) 'KH, DGETRF failed: info=', info, ', i=', i, ', j=', j, ', ke_xy=', ke_xy + call PRC_abort + end if + + ! Solve all right-hand sides using the same LU factors. + call DGETRS( 'N', nv, nrhs1, Amat1, nv, ipiv1, RHS1, nv, info ) + if ( info /= 0 ) then + LOG_ERROR('linkernel_solve_sip',*) 'NU, DGETRS failed: info=', info, ', i=', i, ', j=', j, ', ke_xy=', ke_xy + call PRC_abort + end if + call DGETRS( 'N', nv, nrhs2, Amat2, nv, ipiv2, RHS2, nv, info ) + if ( info /= 0 ) then + LOG_ERROR('linkernel_solve_sip',*) 'KH, DGETRS failed: info=', info, ', i=', i, ', j=', j, ', ke_xy=', ke_xy + call PRC_abort + end if + + if ( top_flag ) then + ! At the top level RHS(:,1:2) contains D^{-1} b. + do iv=1, 2 + do pv1=1, nv + pp = i + (pv1-1)*im + b(pp,iv,j,ke_xy) = RHS1(pv1,iv) + end do + end do + do pv1=1, nv + pp = i + (pv1-1)*im + b(pp,3,j,ke_xy) = RHS2(pv1,1) + end do + + ! Not used in backward substitution. + do pv2=1, nv + do pv1=1, nv + pp = i + (pv1-1)*im + G1(pp,pv2,j,ke_xy) = 0.0_RP + G2(pp,pv2,j,ke_xy) = 0.0_RP + end do + end do + + else + ! RHS(:,1:nv) now contains D^{-1} U. + do pv2=1, nv + do pv1=1, nv + pp = i + (pv1-1)*im + G1(pp,pv2,j,ke_xy) = RHS1(pv1,pv2) + G2(pp,pv2,j,ke_xy) = RHS2(pv1,pv2) + end do + end do + ! RHS(:,nv+1:nv+2) contains D^{-1} b. + do iv=1, 2 + do pv1=1, nv + pp = i + (pv1-1)*im + b(pp,iv,j,ke_xy) = RHS1(pv1,nv+iv) + end do + end do + do pv1=1, nv + pp = i + (pv1-1)*im + b(pp,3,j,ke_xy) = RHS2(pv1,nv+1) + end do + end if + + end do ! i + end do ! j + end do ! ke_xy + !$omp end parallel do + + return + end subroutine linkernel_solve_sip + +!OCL SERIAL + subroutine construct_matbnd_sip( BndMatL, BndMatD, BndMatU, & ! (out) + RHO, KDIFF, GsqrtV, Dx1D, M1D,invM1D, & ! (in) + impl_fac, dt, penalty_fac, lmesh, elem, im, jm, ke_z ) ! (in) + implicit none + class(LocalMesh3D), intent(in) :: lmesh + class(ElementBase3D), intent(in) :: elem + integer, intent(in) :: im, jm + real(RP), intent(out) :: BndMatL(im,elem%Nnode_v,elem%Nnode_v,jm,lmesh%Ne2D) + real(RP), intent(out) :: BndMatD(im,elem%Nnode_v,elem%Nnode_v,jm,lmesh%Ne2D) + real(RP), intent(out) :: BndMatU(im,elem%Nnode_v,elem%Nnode_v,jm,lmesh%Ne2D) + real(RP), intent(in) :: RHO(elem%Np,lmesh%Ne) + real(RP), intent(in) :: KDIFF(elem%Np,lmesh%NeA) + real(RP), intent(in) :: GsqrtV(elem%Np,lmesh%Ne) + real(RP), intent(in) :: Dx1D(elem%Nnode_v,elem%Nnode_v) + real(RP), intent(in) :: M1D(elem%Nnode_v,elem%Nnode_v) + real(RP), intent(in) :: invM1D(elem%Nnode_v,elem%Nnode_v) + real(RP), intent(in) :: impl_fac + real(RP), intent(in) :: dt + real(RP), intent(in) :: penalty_fac + integer, intent(in) :: ke_z + + integer :: ke2D, ke, p + integer :: i, j + integer :: pv, pv1, pv2 + integer :: f + integer :: ke_nb, ke_z_nb + integer :: pvM, pvP + + real(RP) :: lambda + + real(RP) :: mu_loc (elem%Nnode_v) + real(RP) :: rinv_loc(elem%Nnode_v) + real(RP) :: rGsqrtV_loc(elem%Nnode_v) + + real(RP) :: mu_nb (elem%Nnode_v) + real(RP) :: rinv_nb(elem%Nnode_v) + + real(RP) :: Avol (elem%Nnode_v,elem%Nnode_v) + real(RP) :: AffMM(elem%Nnode_v,elem%Nnode_v) + real(RP) :: AffMP(elem%Nnode_v,elem%Nnode_v) + + real(RP) :: MinvALoc(elem%Nnode_v,elem%Nnode_v) + real(RP) :: MinvANb (elem%Nnode_v,elem%Nnode_v) + + real(RP) :: Dz(elem%Nnode_v,elem%Nnode_v) + real(RP) :: Dz_loc(elem%Nnode_v), Dz_nb(elem%Nnode_v) + + real(RP) :: mu_face_M, mu_face_P + real(RP) :: sigma + real(RP) :: Fscale_M + + logical :: boundary_flag + + integer :: Nnode_v + !--------------------------------------------------------------------------- + + call PROF_rapstart('phy_bl_cal_vi_matbnd_sip', 3) + + lambda = impl_fac + Nnode_v = elem%Nnode_v + + !$omp parallel do collapse(2) private( ke, p, ke_nb, ke_z_nb, pvM, pvP, & + !$omp mu_loc, rinv_loc, rGsqrtV_loc, mu_nb, rinv_nb, & + !$omp Avol, AffMM, AffMP, MinvALoc, MinvANb, Dz, Dz_loc, Dz_nb, & + !$omp mu_face_M, mu_face_P, sigma, Fscale_M, boundary_flag ) + do ke2D = 1, lmesh%Ne2D + do j = 1, jm + ke = ke2D + (ke_z-1)*lmesh%Ne2D + + BndMatL(:,:,:,j,ke2D) = 0.0_RP + BndMatD(:,:,:,j,ke2D) = 0.0_RP + BndMatU(:,:,:,j,ke2D) = 0.0_RP + + do i=1, im + do pv=1, Nnode_v + p = i + (j-1)*im + (pv-1)*im*jm + + mu_loc(pv) = RHO(p,ke) * KDIFF(p,ke) + rinv_loc(pv) = 1.0_RP / RHO(p,ke) + rGsqrtV_loc(pv) = 1.0_RP / GsqrtV(p,ke) + + Avol(:,pv) = 0.0_RP + BndMatD(i,pv,pv,j,ke2D) = 1.0_RP + end do + do pv2=1, Nnode_v + do pv1=1, Nnode_v + p = i + (j-1)*im + (pv1-1)*im*jm + Dz(pv1,pv2) = lmesh%Escale(p,ke,3,3) * Dx1D(pv1,pv2) + end do + end do + + !- Volume contribution --- + + ! Avol = Dz^T M diag(mu) Dz diag(1/rho) + do pv2=1, Nnode_v + do pv1=1, Nnode_v + do pv=1, Nnode_v + Avol(pv1,pv2) = Avol(pv1,pv2) & + + Dz(pv,pv1) * M1D(pv,pv) * mu_loc(pv) * rGsqrtV_loc(pv) * Dz(pv,pv2) * rinv_loc(pv2) + end do + end do + end do + + ! Minv Avol + MinvALoc(:,:) = 0.0_RP + do pv2=1, Nnode_v + do pv=1, Nnode_v + do pv1=1, Nnode_v + MinvALoc(pv1,pv2) = MinvALoc(pv1,pv2) + invM1D(pv1,pv) * Avol(pv,pv2) + end do + end do + MinvALoc(pv1,pv2) = rGsqrtV_loc(pv1) * MinvALoc(pv1,pv2) + end do + + do pv2=1, Nnode_v + do pv1=1, Nnode_v + BndMatD(i,pv1,pv2,j,ke2D) = BndMatD(i,pv1,pv2,j,ke2D) + lambda * MinvALoc(pv1,pv2) + end do + end do + + !- Bottom and top faces + + do f=1, 2 + if (f == 1) then + ke_z_nb = max(ke_z-1,1) + pvM = 1; pvP = Nnode_v + boundary_flag = (ke_z == 1) + else + ke_z_nb = min(ke_z+1,lmesh%NeZ) + pvM = Nnode_v; pvP = 1 + boundary_flag = (ke_z == lmesh%NeZ) + end if + ke_nb = ke2D + (ke_z_nb-1)*lmesh%Ne2D + + ! Homogeneous Neumann boundary condition: + ! mu dphi/dn = 0 + ! No SIP boundary contribution is added here. + + if (boundary_flag) cycle + + do pv=1, Nnode_v + p = i + (j-1)*im + (pv-1)*im*jm + mu_nb(pv) = RHO(p,ke_nb) * KDIFF(p,ke_nb) + rinv_nb(pv) = 1.0_RP / RHO(p,ke_nb) + end do + + mu_face_M = mu_loc(pvM) + mu_face_P = mu_nb (pvP) + sigma = penalty_fac * real(Nnode_v, kind=RP)**2 * max(mu_face_M, mu_face_P) + + + !-- + + Dz_loc(:) = Dz(pvM,:) + p = i + (j-1)*im + (pvP-1)*im*jm + do pv=1, Nnode_v + Dz_nb(pv) = lmesh%Escale(p,ke_nb,3,3) * Dx1D(pvP,pv) + end do + + call construct_sip_face_blocks_lgl( AffMM, AffMP, & ! (out) + Dz_loc, Dz_nb, mu_face_M, mu_face_P, rinv_loc, rinv_nb, & ! (in) + sigma, pvM, pvP, f, elem%Nnode_v ) ! (in) + + Fscale_M = lmesh%Fscale(elem%Nfp_h*elem%Nfaces_h+1,ke) + AffMM(:,:) = Fscale_M * AffMM(:,:) + AffMP(:,:) = Fscale_M * AffMP(:,:) + + ! M^-1 AffMM, M^-1 AffMP + MinvALoc(:,:) = 0.0_RP; MinvAnb(:,:) = 0.0_RP + do pv2=1, Nnode_v + do pv=1, Nnode_v + do pv1=1, Nnode_v + MinvALoc(pv1,pv2) = MinvALoc(pv1,pv2) + invM1D(pv1,pv) * AffMM(pv,pv2) + MinvANb (pv1,pv2) = MinvANb (pv1,pv2) + invM1D(pv1,pv) * AffMP(pv,pv2) + end do + end do + end do + + do pv2=1, Nnode_v + do pv1=1, Nnode_v + BndMatD(i,pv1,pv2,j,ke2D) = BndMatD(i,pv1,pv2,j,ke2D) + lambda * rGsqrtV_loc(pv1) * MinvALoc(pv1,pv2) + if ( f == 1 ) then + BndMatL(i,pv1,pv2,j,ke2D) = BndMatL(i,pv1,pv2,j,ke2D) + lambda * rGsqrtV_loc(pv1) * MinvANb(pv1,pv2) + else + BndMatU(i,pv1,pv2,j,ke2D) = BndMatU(i,pv1,pv2,j,ke2D) + lambda * rGsqrtV_loc(pv1) * MinvANb(pv1,pv2) + end if + end do + end do + + end do ! end face loop + end do + end do + end do + + call PROF_rapend('phy_bl_cal_vi_matbnd_sip', 3) + return + end subroutine construct_matbnd_sip + +!OCL SERIAL + subroutine construct_sip_face_blocks_lgl( & + AffMM, AffMP, & ! (out) + dM, dP, muM, muP, rinvM, rinvP, & ! (in) + sigma, pvM, pvP, face_id, nv ) ! (in) + implicit none + integer, intent(in) :: nv + real(RP), intent(out) :: AffMM(nv,nv) + real(RP), intent(out) :: AffMP(nv,nv) + real(RP), intent(in) :: dM(nv), dP(nv) + real(RP), intent(in) :: muM, muP + real(RP), intent(in) :: rinvM(nv), rinvP(nv) + real(RP), intent(in) :: sigma + integer, intent(in) :: pvM, pvP + integer, intent(in) :: face_id + + integer :: i, j + real(RP) :: ncom + !---------------------------------------------------- + + if (face_id == 1) then + ncom = -1.0_RP + else + ncom = +1.0_RP + end if + + AffMM(:,:) = 0.0_RP; AffMP(:,:) = 0.0_RP + do i=1, nv + do j=1, nv + if (i == pvM) then + AffMM(i,j) = AffMM(i,j) & + - 0.5_RP * muM * ncom * dM(j) * rinvM(j) + end if + if (j == pvM) then + AffMM(i,j) = AffMM(i,j) & + - 0.5_RP * muM * ncom * dM(i) * rinvM(j) + end if + if (i == pvM .and. j == pvM) then + AffMM(i,j) = AffMM(i,j) + sigma * rinvM(j) + end if + !- + if ( i == pvM ) then + AffMP(i,j) = AffMP(i,j) & + - 0.5_RP * muP * ncom * dP(j) * rinvP(j) + end if + if ( j == pvP ) then + AffMP(i,j) = AffMP(i,j) & + + 0.5_RP * muM * ncom * dM(i) * rinvP(j) + end if + if ( i == pvM .and. j == pvP ) then + AffMP(i,j) = AffMP(i,j) - sigma * rinvP(j) + end if + end do + end do + return + end subroutine construct_sip_face_blocks_lgl + +!OCL SERIAL + subroutine eval_Ax( MOMX_t, MOMY_t, DRHOT_t, alph_M, alph_H, & + PROG_VARS, MOMX00, MOMY00, PT00, NU, KH, DENS, GsqrtV, impl_fac, dt, & + lmesh, elem, vmapM, vmapP, is_bound, element3D_operation, C_IP, & + im, jm, b, use_delta_form ) + implicit none + class(LocalMesh3D), intent(in) :: lmesh + class(ElementBase3D), intent(in) :: elem + real(RP), intent(out) :: MOMX_t(elem%Np,lmesh%Ne) + real(RP), intent(out) :: MOMY_t(elem%Np,lmesh%Ne) + real(RP), intent(out) :: DRHOT_t(elem%Np,lmesh%Ne) + real(RP), intent(out) :: alph_M(elem%NfpTot,lmesh%Ne) + real(RP), intent(out) :: alph_H(elem%NfpTot,lmesh%Ne) + real(RP), intent(in) :: PROG_VARS(elem%Np,lmesh%NeX*lmesh%NeY*lmesh%NeZ,3) + real(RP), intent(in) :: MOMX00(elem%Np,lmesh%NeA) + real(RP), intent(in) :: MOMY00(elem%Np,lmesh%NeA) + real(RP), intent(in) :: PT00(elem%Np,lmesh%NeA) + real(RP), intent(in) :: NU(elem%Np,lmesh%NeA) + real(RP), intent(in) :: KH(elem%Np,lmesh%NeA) + real(RP), intent(in) :: DENS(elem%Np,lmesh%Ne) + real(RP), intent(in) :: GsqrtV(elem%Np,lmesh%Ne) + real(RP), intent(in) :: impl_fac + real(RP), intent(in) :: dt + integer, intent(in) :: vmapM(elem%NfpTot,lmesh%Ne) + integer, intent(in) :: vmapP(elem%NfpTot,lmesh%Ne) + logical, intent(in) :: is_bound(elem%NfpTot,lmesh%Ne) + class(ElementOperationBase3D), intent(in) :: element3D_operation + real(RP), intent(in) :: C_IP + integer, intent(in) :: im, jm + real(RP), intent(out), optional :: b(im,elem%Nnode_v,3,jm,lmesh%Ne) + logical, intent(in), optional :: use_delta_form + + real(RP) :: Flux(elem%Np,3), DFlux(elem%Np,2,3) + real(RP) :: RDENS(elem%Np) + real(RP) :: RGsqrtV + real(RP) :: E33 + + real(RP) :: DIFF_flux_z_broken(elem%Np,lmesh%NeA,3) + real(RP) :: DIFF_flux_z(elem%Np,lmesh%NeA,3) + + real(RP) :: del_flux(elem%NfpTot,3,lmesh%Ne) + + integer :: ke_xy, ke_z + integer :: ke, ke2D + integer :: p, fp + integer :: iv + + integer :: i, j + integer :: pv + logical :: flag_cal_b + logical :: flag_use_delta_form + !--------------------------------------------------------- + + if ( present(b) .and. present(use_delta_form) ) then + flag_cal_b = .true. + flag_use_delta_form = use_delta_form + else + flag_cal_b = .false. + flag_use_delta_form = .true. + end if + + if ( flag_use_delta_form ) then + call cal_grad_del_flux( del_flux, & ! (out) + PROG_VARS, DENS, & ! (in) + lmesh%normal_fn(:,:,3), lmesh%Fscale, vmapM, vmapP, & ! (in) + lmesh, elem, lmesh%lcmesh2D, lmesh%lcmesh2D%refElem2D ) ! (in) + + !$omp parallel do collapse(2) private( ke, ke2D, DFlux, RDENS, RGsqrtV, E33 ) + do ke_z=1, lmesh%NeZ + do ke_xy=1, lmesh%NeX*lmesh%NeY + ke = ke_xy + (ke_z-1)*lmesh%NeX*lmesh%NeY + ke2D = lmesh%EMap3Dto2D(ke) + + RDENS(:) = 1.0_RP / DENS(:,ke) + do iv=1, 3 + call element3D_operation%Dz( PROG_VARS(:,ke,iv) * RDENS(:), DFlux(:,1,iv) ) + call element3D_operation%Lift( del_flux(:,iv,ke), DFlux(:,2,iv) ) + end do + + do p=1, elem%Np + RGsqrtV = 1.0_RP / GsqrtV(p,ke) + E33 = lmesh%Escale(p,ke,3,3) + DIFF_flux_z_broken(p,ke,1) = DENS(p,ke) * NU(p,ke) * E33 * DFlux(p,1,1) * RGsqrtV + DIFF_flux_z_broken(p,ke,2) = DENS(p,ke) * NU(p,ke) * E33 * DFlux(p,1,2) * RGsqrtV + DIFF_flux_z_broken(p,ke,3) = DENS(p,ke) * KH(p,ke) * E33 * DFlux(p,1,3) * RGsqrtV + + DIFF_flux_z(p,ke,1) = DENS(p,ke) * NU(p,ke) * ( E33 * DFlux(p,1,1) + DFlux(p,2,1) ) * RGsqrtV + DIFF_flux_z(p,ke,2) = DENS(p,ke) * NU(p,ke) * ( E33 * DFlux(p,1,2) + DFlux(p,2,2) ) * RGsqrtV + DIFF_flux_z(p,ke,3) = DENS(p,ke) * KH(p,ke) * ( E33 * DFlux(p,1,3) + DFlux(p,2,3) ) * RGsqrtV + end do + end do + end do + + !--------------------------------------------------------- + + call cal_del_flux( del_flux, alph_M, alph_H, & ! (out) + DIFF_flux_z, DIFF_flux_z_broken, PROG_VARS, DENS, NU, KH, C_IP, & ! (in) + lmesh%normal_fn(:,:,3), lmesh%Fscale, vmapM, vmapP, is_bound, & ! (in) + lmesh, elem, lmesh%lcmesh2D, lmesh%lcmesh2D%refElem2D ) ! (in) + + !$omp parallel do collapse(2) private( ke, ke2D, DFlux, RGsqrtV, E33 ) + do ke_z=1, lmesh%NeZ + do ke_xy=1, lmesh%NeX*lmesh%NeY + ke = ke_xy + (ke_z-1)*lmesh%NeX*lmesh%NeY + ke2D = lmesh%EMap3Dto2D(ke) + do iv=1, 3 + call element3D_operation%Dz( DIFF_flux_z(:,ke,iv), DFlux(:,1,iv) ) + call element3D_operation%Lift( del_flux(:,iv,ke), DFlux(:,2,iv) ) + end do + + do p=1, elem%Np + RGsqrtV = 1.0_RP / GsqrtV(p,ke) + E33 = lmesh%Escale(p,ke,3,3) + MOMX_t (p,ke) = ( E33 * DFlux(p,1,1) + DFlux(p,2,1) ) * RGsqrtV + MOMY_t (p,ke) = ( E33 * DFlux(p,1,2) + DFlux(p,2,2) ) * RGsqrtV + DRHOT_t(p,ke) = ( E33 * DFlux(p,1,3) + DFlux(p,2,3) ) * RGsqrtV + end do + end do + end do + end if + + if ( flag_cal_b ) then + if ( use_delta_form ) then + !$omp parallel do private(p) + do ke=1, lmesh%Ne + do pv=1, elem%Nnode_v + do j=1, jm + do i=1, im + p = i + (j-1)*im + (pv-1)*im*jm + b(i,pv,1,j,ke) = impl_fac * MOMX_t(p,ke) & + - PROG_VARS (p,ke,1) & + + MOMX00(p,ke) + b(i,pv,2,j,ke) = impl_fac * MOMY_t(p,ke) & + - PROG_VARS (p,ke,2) & + + MOMY00(p,ke) + b(i,pv,3,j,ke) = impl_fac * DRHOT_t(p,ke) & + - PROG_VARS (p,ke,3) & + + DENS(p,ke) * PT00(p,ke) + end do + end do + end do + end do + else + !$omp parallel do private(p) + do ke=1, lmesh%Ne + do pv=1, elem%Nnode_v + do j=1, jm + do i=1, im + p = i + (j-1)*im + (pv-1)*im*jm + b(i,pv,1,j,ke) = MOMX00(p,ke) + b(i,pv,2,j,ke) = MOMY00(p,ke) + b(i,pv,3,j,ke) = DENS(p,ke) * PT00(p,ke) + end do + end do + end do + end do + end if + end if + return + end subroutine eval_Ax + +!OCL SERIAL + subroutine cal_del_flux( del_flux, alph_M, alph_H, & + DIFF_flux_z, DIFF_flux_z_broken, PROG_VARS, DENS, NU, KH, C_IP, & + nz, Fscale, vmapM, vmapP, is_bound, lmesh, elem, lmesh2D, elem2D ) + implicit none + class(LocalMesh3D), intent(in) :: lmesh + class(ElementBase3D), intent(in) :: elem + class(LocalMesh2D), intent(in) :: lmesh2D + class(ElementBase2D), intent(in) :: elem2D + real(RP), intent(out) :: del_flux(elem%NfpTot,3,lmesh%Ne) + real(RP), intent(out) :: alph_M(elem%NfpTot,lmesh%Ne) + real(RP), intent(out) :: alph_H(elem%NfpTot,lmesh%Ne) + real(RP), intent(in) :: DIFF_flux_z(elem%Np*lmesh%NeA,3) + real(RP), intent(in) :: DIFF_flux_z_broken(elem%Np*lmesh%NeA,3) + real(RP), intent(in) :: PROG_VARS(elem%Np*lmesh%NeX*lmesh%NeY*lmesh%NeZ,3) + real(RP), intent(in) :: DENS(elem%Np*lmesh%NeA) + real(RP), intent(in) :: NU(elem%Np*lmesh%NeA) + real(RP), intent(in) :: KH(elem%Np*lmesh%NeA) + real(RP), intent(in) :: C_IP + real(RP), intent(in) :: nz(elem%NfpTot,lmesh%Ne) + real(RP), intent(in) :: Fscale(elem%NfpTot,lmesh%Ne) + integer, intent(in) :: vmapM(elem%NfpTot,lmesh%Ne) + integer, intent(in) :: vmapP(elem%NfpTot,lmesh%Ne) + logical, intent(in) :: is_bound(elem%NfpTot,lmesh%Ne) + + integer :: ke, ke_z, ke_xy + integer :: iP(elem%NfpTot), iM(elem%NfpTot) + real(RP) :: DIFF_flux_z_P(elem%NfpTot,3) + + real(RP) :: coef(elem%NfpTot) + real(RP) :: RDENS_M(elem%NfpTot), RDENS_P(elem%NfpTot) + real(RP) :: numflux(elem%NfpTot,3) + !------------------------------------------------ + + !$omp parallel do collapse(2) & + !$omp private( ke, iM, iP, alph_M, alph_H, coef, RDENS_M, RDENS_P, numflux ) + do ke_z=1, lmesh%NeZ + do ke_xy=1, lmesh%NeX*lmesh%NeY + ke = ke_xy + (ke_z-1)*lmesh%Ne2D + iM(:) = vmapM(:,ke); iP(:) = vmapP(:,ke) + + alph_M(:,ke) = C_IP * (elem%Nnode_v)**2 * max( DENS(iM)*NU(iM), DENS(iP)*NU(iP) ) + alph_H(:,ke) = C_IP * (elem%Nnode_v)**2 * max( DENS(iM)*KH(iM), DENS(iP)*KH(iP) ) + + coef(:) = nz(:,ke) * Fscale(:,ke) + RDENS_M(:) = 1.0_RP / DENS(iM) + RDENS_P(:) = 1.0_RP / DENS(iP) + + where ( is_bound(:,ke) ) + numflux(:,1) = 0.0_RP + elsewhere + numflux(:,1) = 0.5_RP * ( DIFF_flux_z_broken(iP,1) + DIFF_flux_z(iM,1) ) + numflux(:,2) = 0.5_RP * ( DIFF_flux_z_broken(iP,2) + DIFF_flux_z(iM,2) ) + numflux(:,3) = 0.5_RP * ( DIFF_flux_z_broken(iP,3) + DIFF_flux_z(iM,3) ) + end where + + !- + del_flux(:,1,ke) = coef(:) * ( numflux(:,1) - DIFF_flux_z(iM,1) ) & + + alph_M(:,ke) * Fscale(:,ke) * ( PROG_VARS(iP,1) * RDENS_P(:)- PROG_VARS(iM,1) * RDENS_M(:) ) + + del_flux(:,2,ke) = coef(:) * ( numflux(:,2) - DIFF_flux_z(iM,2) ) & + + alph_M(:,ke) * Fscale(:,ke) * ( PROG_VARS(iP,2) * RDENS_P(:) - PROG_VARS(iM,2) * RDENS_M(:) ) + + del_flux(:,3,ke) = coef(:) * ( numflux(:,3) - DIFF_flux_z(iM,3) ) & + + alph_H(:,ke) * Fscale(:,ke) * ( PROG_VARS(iP,3) * RDENS_P(:) - PROG_VARS(iM,3) * RDENS_M(:) ) + end do + end do + + return + end subroutine cal_del_flux + +!OCL SERIAL + subroutine cal_grad_del_flux( del_flux, & + PROG_VARS, DENS, nz, Fscale,vmapM, vmapP, & + lmesh, elem, lmesh2D, elem2D ) + implicit none + class(LocalMesh3D), intent(in) :: lmesh + class(ElementBase3D), intent(in) :: elem + class(LocalMesh2D), intent(in) :: lmesh2D + class(ElementBase2D), intent(in) :: elem2D + real(RP), intent(out) :: del_flux(elem%NfpTot,3,lmesh%Ne) + real(RP), intent(in) :: PROG_VARS(elem%Np*lmesh%NeX*lmesh%NeY*lmesh%NeZ,3) + real(RP), intent(in) :: DENS(elem%Np*lmesh%Ne) + real(RP), intent(in) :: nz(elem%NfpTot,lmesh%Ne) + real(RP), intent(in) :: Fscale(elem%NfpTot,lmesh%Ne) + integer, intent(in) :: vmapM(elem%NfpTot,lmesh%Ne) + integer, intent(in) :: vmapP(elem%NfpTot,lmesh%Ne) + + integer :: ke, ke_z, ke_xy + integer :: iP(elem%NfpTot), iM(elem%NfpTot) + + real(RP) :: coef(elem%NfpTot) + real(RP) :: RDENS_M(elem%NfpTot), RDENS_P(elem%NfpTot) + !------------------------------------------------ + + !$omp parallel do collapse(2) & + !$omp private( ke, iM, iP, coef, RDENS_M, RDENS_P ) + do ke_z=1, lmesh%NeZ + do ke_xy=1, lmesh%NeX*lmesh%NeY + ke = ke_xy + (ke_z-1)*lmesh%Ne2D + iM(:) = vmapM(:,ke); iP(:) = vmapP(:,ke) + + coef(:) = 0.5_RP * nz(:,ke) * Fscale(:,ke) + RDENS_M(:) = 1.0_RP / DENS(iM) + RDENS_P(:) = 1.0_RP / DENS(iP) + + del_flux(:,1,ke) = coef(:) * ( PROG_VARS(iP,1) * RDENS_P(:) - PROG_VARS(iM,1) * RDENS_M(:) ) + del_flux(:,2,ke) = coef(:) * ( PROG_VARS(iP,2) * RDENS_P(:) - PROG_VARS(iM,2) * RDENS_M(:) ) + del_flux(:,3,ke) = coef(:) * ( PROG_VARS(iP,3) * RDENS_P(:) - PROG_VARS(iM,3) * RDENS_M(:) ) + end do + end do + + return + end subroutine cal_grad_del_flux +end module scale_atm_phy_bl_dgm_common \ No newline at end of file diff --git a/FElib/src/bl_turbulence/scale_atm_phy_bl_dgm_mynn_lv2.F90 b/FElib/src/bl_turbulence/scale_atm_phy_bl_dgm_mynn_lv2.F90 new file mode 100644 index 00000000..b3b7ecba --- /dev/null +++ b/FElib/src/bl_turbulence/scale_atm_phy_bl_dgm_mynn_lv2.F90 @@ -0,0 +1,395 @@ +!> module FElib / Atmosphere / Physics / boundary layer turbulence +!! +!! @par Description +!! A module to provide routines to calculate turbulent diffusion coefficients based on the Mellor-Yamada-Nakanishi-Niino (MYNN) Level 2 closure. +!! +!! This implementation uses the Level 2 diagnostic closure only. Prognostic turbulent kinetic energy, nonlocal mixing, and the MY Level 2.5 equations are not included. +!! +!! @par Reference +!! Mellor, G. L. and T. Yamada, 1982: +!! Development of a turbulence closure model for geophysical fluid problems. +!! Reviews of Geophysics and Space Physics, 20, 851-875. +!! +!! Nakanishi, M. and H. Niino, 2009: +!! Development of an improved turbulence closure model for the atmospheric boundary layer. +!! Journal of the Meteorological Society of Japan, 87, 895-912. +!! +!! Blackadar, A. K., 1962: +!! The vertical distribution of wind and turbulent exchange in a neutral atmosphere. +!! Journal of Geophysical Research, 67, 3095-3102. +!! +!! +!! @author Yuta Kawai, Team SCALE +!! +!< +!------------------------------------------------------------------------------- +#include "scaleFElib.h" +module scale_atm_phy_bl_dgm_mynn_lv2 + !----------------------------------------------------------------------------- + ! + !++ Used modules + ! + use scale_precision + use scale_io + use scale_prc + use scale_prof + use scale_const, only: & + EPS => CONST_EPS, & + GRAV => CONST_GRAV, & + Rdry => CONST_Rdry, & + CPdry => CONST_CPdry, & + CVdry => CONST_CVdry, & + PRES00 => CONST_PRE00, & + KARMAN => CONST_KARMAN, & + RPlanet => CONST_RADIUS + use scale_sparsemat + use scale_element_base, only: & + ElementBase2D, ElementBase3D + use scale_element_hexahedral, only: HexahedralElement + use scale_localmesh_2d, only: LocalMesh2D + use scale_localmesh_3d, only: LocalMesh3D + use scale_mesh_base3d, only: MeshBase3D + use scale_localmeshfield_base, only: LocalMeshField3D + use scale_meshfield_base, only: MeshField3D + + !----------------------------------------------------------------------------- + implicit none + private + !----------------------------------------------------------------------------- + ! + !++ Public procedures + ! + public :: atm_phy_bl_dgm_mynn_lv2_Init + public :: atm_phy_bl_dgm_mynn_lv2_Final + public :: atm_phy_bl_dgm_mynn_lv2_cal_VViscDiffCoef + + !----------------------------------------------------------------------------- + ! + !++ Private procedure + ! + + !----------------------------------------------------------------------------- + ! + !++ Private parameters & variables + ! + + !- MYNN closure constants ----------- + + real(RP), parameter :: A1 = 1.18_RP + real(RP), parameter :: A2 = 0.665_RP + + real(RP), parameter :: B1 = 24.0_RP + real(RP), parameter :: B2 = 15.0_RP + + real(RP), parameter :: C1 = 0.137_RP + real(RP), parameter :: C2 = 0.75_RP + real(RP), parameter :: C3 = 0.352_RP + real(RP), parameter :: C5 = 0.2_RP + + real(RP), parameter :: G1 = 0.235_RP + + ! Derived constants + real(RP) :: G2 + real(RP) :: F2 + real(RP) :: Rf2 + real(RP) :: RFc + + ! Limiter with mixing length + real(RP) :: L_INF = 100.0_RP + real(RP) :: L_MIN = 1.0E-6_RP + +contains + !> Initialize a module of MYNN Level 2 PBL turbulence parameterization +!OCL SERIAL + subroutine atm_phy_bl_dgm_mynn_lv2_Init( mesh ) + implicit none + class(MeshBase3D), intent(in) :: mesh + + namelist / PARAM_ATMOS_PHY_BL_DGM_MYNN_LV2 / & + L_INF, L_MIN + + integer :: ierr + !-------------------------------------------------------------------------------- + + LOG_NEWLINE + LOG_INFO("ATMOS_PHY_BL_DGM_MYNN_LV2_setup",*) 'Setup' + LOG_INFO("ATMOS_PHY_BL_DGM_MYNN_LV2_setup",*) 'MYNN Level 2 PBL turbulence parameterization' + + !--- read namelist + rewind(IO_FID_CONF) + read(IO_FID_CONF,nml=PARAM_ATMOS_PHY_BL_DGM_MYNN_LV2,iostat=ierr) + if( ierr < 0 ) then !--- missing + LOG_INFO("ATMOS_PHY_BL_DGM_MYNN_LV2_setup",*) 'Not found namelist. Default used.' + elseif( ierr > 0 ) then !--- fatal error + LOG_ERROR("ATMOS_PHY_BL_DGM_MYNN_LV2_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_BL_DGM_MYNN_LV2. Check!' + call PRC_abort + endif + LOG_NML(PARAM_ATMOS_PHY_BL_DGM_MYNN_LV2) + + !- + + G2 = ( 2.0_RP * A1 * ( 3.0_RP - 2.0_RP * C2 ) & + + B2 * ( 1.0_RP - C3 ) & + ) / B1 + + F2 = B1 * ( G1 + G2 ) & + - 3.0_RP * A1 * ( 1.0_RP - C2 ) + + Rf2 = B1 * G1 / F2 + + RFc = G1 / ( G1 + G2 ) + return + end subroutine atm_phy_bl_dgm_mynn_lv2_Init + + !> Finalize a module of MYNN Level 2 PBL turbulence parameterization +!OCL SERIAL + subroutine atm_phy_bl_dgm_mynn_lv2_Final() + implicit none + !-------------------------------------------------------------------------------- + return + end subroutine atm_phy_bl_dgm_mynn_lv2_Final + +!OCL SERIAL + subroutine atm_phy_bl_dgm_mynn_lv2_cal_VViscDiffCoef( & + Nu, Kh, TKE, & ! (out) + DDENS_, MOMX_, MOMY_, MOMZ_, DRHOT_, DENS_hyd, PRES_hyd, & ! (in) + PRES, PT, & ! (in) + Dz, Lift, lmesh, elem, is_bound ) ! (in) + implicit none + class(LocalMesh3D), intent(in) :: lmesh + class(ElementBase3D), intent(in) :: elem + real(RP), intent(out) :: Nu(elem%Np,lmesh%NeA) !< Vertical eddy viscosity + real(RP), intent(out) :: Kh(elem%Np,lmesh%NeA) !< Vertical eddy diffusivity + real(RP), intent(out) :: TKE(elem%Np,lmesh%NeA) !< Turbulent kinetic energy + real(RP), intent(in) :: DDENS_(elem%Np,lmesh%NeA) !< Density perturbation + real(RP), intent(in) :: MOMX_ (elem%Np,lmesh%NeA) !< Momentum in x1 direction + real(RP), intent(in) :: MOMY_ (elem%Np,lmesh%NeA) !< Momentum in x2 direction + real(RP), intent(in) :: MOMZ_ (elem%Np,lmesh%NeA) !< Momentum in x3 direction + real(RP), intent(in) :: DRHOT_(elem%Np,lmesh%NeA) !< Density x potential temperature perturbation + real(RP), intent(in) :: DENS_hyd(elem%Np,lmesh%NeA) !< Reference density in hydrostatic balance + real(RP), intent(in) :: PRES_hyd(elem%Np,lmesh%NeA) !< Reference pressure in hydrostatic balance + real(RP), intent(in) :: PRES(elem%Np,lmesh%NeA) !< Pressure + real(RP), intent(in) :: PT(elem%Np,lmesh%NeA) !< Potential temperature + type(SparseMat), intent(in) :: Dz, Lift + logical, intent(in) :: is_bound(elem%NfpTot,lmesh%Ne) !< Flag whether nodes are located at domain boundaries + + real(RP) :: Fz(elem%Np), LiftDelFlx(elem%Np) + real(RP) :: DENS(elem%Np), RDENS(elem%Np), RHOT(elem%Np), Q(elem%Np) + real(RP) :: DdensDz(elem%Np), DVelDz(elem%Np,3), DptDz(elem%Np) + + real(RP) :: del_flux_rho (elem%NfpTot,lmesh%Ne) + real(RP) :: del_flux_mom (elem%NfpTot,lmesh%Ne,3) + real(RP) :: del_flux_rhot(elem%NfpTot,lmesh%Ne) + + real(RP) :: S2(elem%Np) ! Square of vertical shear of horizontal wind + real(RP) :: Ri ! Richardson number + real(RP) :: Rf(elem%Np) ! Flux Richardson number + + real(RP) :: A2_loc + real(RP) :: F1 + real(RP) :: Rf1 + real(RP) :: AF12 + real(RP) :: discriminant + real(RP) :: denom_h + real(RP) :: denom_m + real(RP) :: S_M(elem%Np), S_H(elem%Np) ! Stability functions for momentum and heat + + real(RP) :: mixlen(elem%Np) ! Mixing length + real(RP) :: zsfc(elem%Nnode_h1D**2,lmesh%Ne2D) ! Surface height + real(RP) :: kz + + integer :: ke, ke2D + integer :: p, ph, pz + class(LocalMesh2D), pointer :: lmesh2D + class(ElementBase2D), pointer :: elem2D + + real(RP), parameter :: EPS_RF = 1.0e-10_RP + real(RP), parameter :: EPS_DISC = 1.0e-14_RP + !-------------------------------------------------------------------------------- + + ! For standard case of NN2009, these parameter are independent of Ri. + ! If we will implement the K2010 correction in the future, + ! these parameter should be calculated in the ke loop because they depend on Ri (i.e., A2_loc = A_2/(1+Ri)). + A2_loc = A2 + + F1 = B1 * ( G1 - C1 ) & + + 2.0_RP * A1 * ( 3.0_RP - 2.0_RP * C2 ) & + + 3.0_RP * A2_loc * ( 1.0_RP - C3 ) * ( 1.0_RP - C5 ) + + Rf1 = B1 * ( G1 - C1 ) / F1 + + AF12 = A1 * F1 / ( A2_loc * F2 ) + + lmesh2D => lmesh%lcmesh2D + elem2D => lmesh2D%refElem2D + + !--------------------------------- + + call cal_del_flux_grad( del_flux_rho, del_flux_mom, del_flux_rhot, & ! (out) + DDENS_, MOMX_, MOMY_, MOMZ_, DRHOT_, DENS_hyd, PRES_hyd, PT, & ! (in) + lmesh%normal_fn(:,:,3), lmesh%vmapM, lmesh%vmapP, & ! (in) + lmesh, elem, is_bound ) ! (in) + + !$omp parallel & + !$omp private( Fz, LiftDelFlx, DENS, RDENS, RHOT, Q, DdensDz, DVelDz, DptDz, & + !$omp S2, Ri, Rf, discriminant, denom_m, denom_h, S_M, S_H, mixlen, kz ) + + !$omp do + do ke2D=lmesh2D%NeS, lmesh2D%NeE + zsfc(:,ke2D) = lmesh%zlev(elem%Hslice(:,1),ke2D) + end do + !$omp end do + !$omp do + do ke=lmesh%NeS, lmesh%NeE + DENS (:) = DENS_hyd(:,ke) + DDENS_(:,ke) + RDENS(:) = 1.0_RP / DENS(:) + RHOT(:) = DENS(:) * PT(:,ke) + + ! gradient of density + call sparsemat_matmul( Dz, DENS, Fz ) + call sparsemat_matmul( Lift, lmesh%Fscale(:,ke) * del_flux_rho(:,ke), LiftDelFlx ) + DdensDz(:) = lmesh%Escale(:,ke,3,3) * Fz(:) + LiftDelFlx(:) + + ! gradient of u + Q(:) = MOMX_(:,ke) * RDENS(:) + call sparsemat_matmul( Dz, MOMX_(:,ke), Fz ) + call sparsemat_matmul( Lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,1), LiftDelFlx ) + DVelDz(:,1) = ( lmesh%Escale(:,ke,3,3) * Fz(:) + LiftDelFlx(:) - Q(:) * DdensDz(:) ) * RDENS(:) + + ! gradient of v + Q(:) = MOMY_(:,ke) * RDENS(:) + call sparsemat_matmul( Dz, MOMY_(:,ke), Fz ) + call sparsemat_matmul( Lift, lmesh%Fscale(:,ke) * del_flux_mom(:,ke,2), LiftDelFlx ) + DVelDz(:,2) = ( lmesh%Escale(:,ke,3,3) * Fz(:) + LiftDelFlx(:) - Q(:) * DdensDz(:) ) * RDENS(:) + + ! gradient of pt + Q(:) = RHOT(:) * RDENS(:) + call sparsemat_matmul( Dz, RHOT, Fz ) + call sparsemat_matmul( Lift, lmesh%Fscale(:,ke) * del_flux_rhot(:,ke), LiftDelFlx ) + DptDz(:) = ( lmesh%Escale(:,ke,3,3) * Fz(:) + LiftDelFlx(:) - Q(:) * DdensDz(:) ) * RDENS(:) + + ! Calculate flux Richardson number: Rf + do p=1, elem%Np + S2(p) = DVelDz(p,1)**2 + DVelDz(p,2)**2 + Ri = GRAV * DptDz(p) / ( PT(p,ke) * max(S2(p), EPS) ) + + discriminant = Ri * Ri & + + 2.0_RP * AF12 * ( Rf1 - 2.0_RP * Rf2 ) * Ri & + + ( AF12 * Rf1 )**2 + discriminant = max( discriminant, 0.0_RP ) + + Rf(p) = 0.5_RP / AF12 * ( Ri + AF12 * Rf1 - sqrt(discriminant) ) + Rf(p) = min( Rf(p), RFc - EPS_RF ) + end do + + ! Calculate stability function: S_M and S_H + do p=1, elem%Np + denom_h = max( 1.0_RP - Rf(p), EPS_RF ) + S_H(p) = 3.0_RP * A2_loc * ( G1 + G2 ) * ( RFc - Rf(p) ) / denom_h + + denom_m = Rf2 - Rf(p) + if ( abs(denom_m) < EPS_RF ) then + denom_m = sign(EPS_RF, denom_m) + end if + S_M(p) = S_H(p) * AF12 * ( Rf1 - Rf(p) ) / denom_m + + S_M(p) = max( S_M(p), 0.0_RP ) + S_H(p) = max( S_H(p), 0.0_RP ) + end do + + ! Calculate mixing length + + ke2D = lmesh%EMap3Dto2D(ke) + do pz=1, elem%Nnode_v + do ph=1, elem%Nnode_h1D**2 + p = ph + (pz-1)*elem%Nnode_h1D**2 + kz = KARMAN * ( lmesh%zlev(p,ke) - zsfc(ph,ke2D) ) + + mixlen(p) = kz * L_INF / max( kz + L_INF, L_MIN ) + end do + end do + + ! Calculate eddy viscosity and diffusivity + do p=1, elem%Np + ! q^2 + q(p) = B1 * mixlen(p)**2 & + * S_M(p) * max( 1.0_RP - Rf(p), 0.0_RP ) * S2(p) + TKE(p,ke) = 0.5_RP * max( q(p), 0.0_RP ) + + q(p) = sqrt( 2.0_RP * TKE(p,ke) ) + Kh(p,ke) = mixlen(p) * q(p) * S_H(p) + Nu(p,ke) = mixlen(p) * q(p) * S_M(p) + end do + end do + !$omp end do + !$omp end parallel + + return + end subroutine atm_phy_bl_dgm_mynn_lv2_cal_VViscDiffCoef + +!-- private -------------------------------------------------------- + +!OCL SERIAL + subroutine cal_del_flux_grad( del_flux_rho, del_flux_mom, del_flux_rhot, & ! (out) + DDENS_, MOMX_, MOMY_, MOMZ_, DRHOT_, DENS_hyd, PRES_hyd, PT_, & ! (in) + nz, vmapM, vmapP, lmesh, elem, is_bound ) ! (in) + + implicit none + + class(LocalMesh3D), intent(in) :: lmesh + class(ElementBase3D), intent(in) :: elem + real(RP), intent(out) :: del_flux_rho(elem%NfpTot*lmesh%Ne) + real(RP), intent(out) :: del_flux_mom(elem%NfpTot*lmesh%Ne,3) + real(RP), intent(out) :: del_flux_rhot(elem%NfpTot*lmesh%Ne) + real(RP), intent(in) :: DDENS_(elem%Np*lmesh%NeA) + real(RP), intent(in) :: MOMX_(elem%Np*lmesh%NeA) + real(RP), intent(in) :: MOMY_(elem%Np*lmesh%NeA) + real(RP), intent(in) :: MOMZ_(elem%Np*lmesh%NeA) + real(RP), intent(in) :: DRHOT_(elem%Np*lmesh%NeA) + real(RP), intent(in) :: DENS_hyd(elem%Np*lmesh%NeA) + real(RP), intent(in) :: PRES_hyd(elem%Np*lmesh%NeA) + real(RP), intent(in) :: PT_(elem%Np*lmesh%NeA) + real(RP), intent(in) :: nz(elem%NfpTot*lmesh%Ne) + integer, intent(in) :: vmapM(elem%NfpTot*lmesh%Ne) + integer, intent(in) :: vmapP(elem%NfpTot*lmesh%Ne) + logical, intent(in) :: is_bound(elem%NfpTot*lmesh%Ne) + + integer :: i, iP, iM + real(RP) :: densM, densP + real(RP) :: del + real(RP) :: facz + !------------------------------------------------------------------------ + + !$omp parallel do private ( iM, iP, densM, densP, del, facz ) + do i=1, elem%NfpTot * lmesh%Ne + iM = vmapM(i); iP = vmapP(i) + + densM = DDENS_(iM) + DENS_hyd(iM) + densP = DDENS_(iP) + DENS_hyd(iP) + + if ( is_bound(i) ) then + facz = 1.0_RP + else + ! facz = 1.0_RP - sign(1.0_RP,nz(i)) + facz = 1.0_RP + end if + + del = 0.5_RP * ( densP - densM ) + del_flux_rho(i) = facz * del * nz(i) + + del = 0.5_RP * ( MOMX_(iP) - MOMX_(iM) ) + del_flux_mom(i,1) = facz * del * nz(i) + + del = 0.5_RP * ( MOMY_(iP) - MOMY_(iM) ) + del_flux_mom(i,2) = facz * del * nz(i) + + del = 0.5_RP * ( MOMZ_(iP) - MOMZ_(iM) ) + del_flux_mom(i,3) = facz * del * nz(i) + + del = 0.5_RP * ( densP * PT_(iP) - densM * PT_(iM) ) + del_flux_rhot(i) = facz * del * nz(i) + end do + + return + end subroutine cal_del_flux_grad +end module scale_atm_phy_bl_dgm_mynn_lv2 \ No newline at end of file diff --git a/FElib/src/depend b/FElib/src/depend index 2f599fb8..d6ac5b5a 100644 --- a/FElib/src/depend +++ b/FElib/src/depend @@ -32,6 +32,8 @@ $(BUILD_DIR)/scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_numflux.o: fluid_dyn_solver/ $(BUILD_DIR)/scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform.o: fluid_dyn_solver/scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform.F90 $(DEPENDLIB) $(BUILD_DIR)/scale_atm_dyn_dgm_nonhydro3d_common.o $(BUILD_DIR)/scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common.o $(BUILD_DIR)/scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_numflux.o $(BUILD_DIR)/scale_element_base.o $(BUILD_DIR)/scale_element_hexahedral.o $(BUILD_DIR)/scale_element_modalfilter.o $(BUILD_DIR)/scale_element_operation_base.o $(BUILD_DIR)/scale_linalgebra.o $(BUILD_DIR)/scale_localmesh_2d.o $(BUILD_DIR)/scale_localmesh_3d.o $(BUILD_DIR)/scale_localmeshfield_base.o $(BUILD_DIR)/scale_mesh_base3d.o $(BUILD_DIR)/scale_meshfield_base.o $(BUILD_DIR)/scale_sparsemat.o $(BUILD_DIR)/scale_atm_dyn_dgm_spongelayer.o: fluid_dyn_solver/scale_atm_dyn_dgm_spongelayer.F90 $(DEPENDLIB) $(BUILD_DIR)/scale_element_base.o $(BUILD_DIR)/scale_localmesh_3d.o $(BUILD_DIR)/scale_localmesh_base.o $(BUILD_DIR)/scale_mesh_base3d.o $(BUILD_DIR)/scale_mesh_cubedom3d.o $(BUILD_DIR)/scale_mesh_cubedspheredom3d.o $(BUILD_DIR)/scale_atm_dyn_dgm_trcadvect3d_heve.o: fluid_dyn_solver/scale_atm_dyn_dgm_trcadvect3d_heve.F90 $(DEPENDLIB) $(BUILD_DIR)/scale_element_base.o $(BUILD_DIR)/scale_element_hexahedral.o $(BUILD_DIR)/scale_element_operation_base.o $(BUILD_DIR)/scale_localmesh_2d.o $(BUILD_DIR)/scale_localmesh_3d.o $(BUILD_DIR)/scale_localmeshfield_base.o $(BUILD_DIR)/scale_mesh_base2d.o $(BUILD_DIR)/scale_mesh_base3d.o $(BUILD_DIR)/scale_meshfield_base.o $(BUILD_DIR)/scale_polynomial.o $(BUILD_DIR)/scale_sparsemat.o +$(BUILD_DIR)/scale_atm_phy_bl_dgm_common.o: bl_turbulence/scale_atm_phy_bl_dgm_common.F90 $(DEPENDLIB) $(BUILD_DIR)/scale_atm_dyn_dgm_hevi_common_linalgebra.o $(BUILD_DIR)/scale_element_base.o $(BUILD_DIR)/scale_element_hexahedral.o $(BUILD_DIR)/scale_element_operation_base.o $(BUILD_DIR)/scale_localmesh_2d.o $(BUILD_DIR)/scale_localmesh_3d.o $(BUILD_DIR)/scale_localmeshfield_base.o $(BUILD_DIR)/scale_mesh_base3d.o $(BUILD_DIR)/scale_meshfield_base.o $(BUILD_DIR)/scale_sparsemat.o +$(BUILD_DIR)/scale_atm_phy_bl_dgm_mynn_lv2.o: bl_turbulence/scale_atm_phy_bl_dgm_mynn_lv2.F90 $(DEPENDLIB) $(BUILD_DIR)/scale_element_base.o $(BUILD_DIR)/scale_element_hexahedral.o $(BUILD_DIR)/scale_localmesh_2d.o $(BUILD_DIR)/scale_localmesh_3d.o $(BUILD_DIR)/scale_localmeshfield_base.o $(BUILD_DIR)/scale_mesh_base3d.o $(BUILD_DIR)/scale_meshfield_base.o $(BUILD_DIR)/scale_sparsemat.o $(BUILD_DIR)/scale_atm_phy_mp_dgm_common.o: microphysics/scale_atm_phy_mp_dgm_common.F90 $(DEPENDLIB) $(BUILD_DIR)/scale_element_base.o $(BUILD_DIR)/scale_localmesh_3d.o $(BUILD_DIR)/scale_localmeshfield_base.o $(BUILD_DIR)/scale_mesh_base3d.o $(BUILD_DIR)/scale_polynomial.o $(BUILD_DIR)/scale_sparsemat.o $(BUILD_DIR)/scale_atm_phy_mp_lscond.o: microphysics/scale_atm_phy_mp_lscond.F90 $(DEPENDLIB) $(BUILD_DIR)/scale_element_base.o $(BUILD_DIR)/scale_localmesh_3d.o $(BUILD_DIR)/scale_mesh_base3d.o $(BUILD_DIR)/scale_atm_phy_rd_dgm_common.o: radiation/scale_atm_phy_rd_dgm_common.F90 $(DEPENDLIB) $(BUILD_DIR)/scale_element_base.o $(BUILD_DIR)/scale_element_hexahedral.o $(BUILD_DIR)/scale_element_operation_base.o $(BUILD_DIR)/scale_localmesh_2d.o $(BUILD_DIR)/scale_localmesh_3d.o $(BUILD_DIR)/scale_localmeshfield_base.o $(BUILD_DIR)/scale_mesh_base3d.o $(BUILD_DIR)/scale_meshfield_base.o diff --git a/FElib/src/file/scale_file_base_meshfield.F90 b/FElib/src/file/scale_file_base_meshfield.F90 index d5e7ffec..40e1bfb9 100644 --- a/FElib/src/file/scale_file_base_meshfield.F90 +++ b/FElib/src/file/scale_file_base_meshfield.F90 @@ -1200,32 +1200,39 @@ subroutine write_axes( this, & ! (in) real(RP), allocatable :: x(:) real(RP), allocatable :: y(:) real(RP), allocatable :: z(:) - - logical :: force_uniform_grid !------------ - if ( this%mesh_type_id == MESHTYPE_2D_CUBEDSPHERE .or. this%mesh_type_id == MESHTYPE_3D_CUBEDSPHERE ) then - force_uniform_grid = .false. - else - force_uniform_grid = this%force_uniform_grid - end if - select case ( this%mesh_type_id ) case ( MESHTYPE_1D ) ! 1D mesh allocate( x(this%dimsinfo(1)%size) ) - call File_common_meshfield_get_axis( this%mesh1D, this%dimsinfo, x(:), force_uniform_grid ) - - call FILE_Write_Axis( fid, this%dimsinfo(1)%name, x(:), start(1:1) ) case ( MESHTYPE_2D_RECTDOM, MESHTYPE_2D_CUBEDSPHERE ) ! 2D mesh allocate( x(this%dimsinfo(1)%size), y(this%dimsinfo(2)%size) ) - call File_common_meshfield_get_axis( this%mesh2D, this%dimsinfo, x(:), y(:), force_uniform_grid ) + case ( MESHTYPE_3D_CUBEDOM, MESHTYPE_3D_CUBEDSPHERE ) ! 3D mesh + allocate( x(this%dimsinfo(1)%size), y(this%dimsinfo(2)%size), z(this%dimsinfo(3)%size) ) + end select + !- + select case ( this%mesh_type_id ) + case ( MESHTYPE_1D ) ! 1D mesh + call File_common_meshfield_get_axis( this%mesh1D, this%dimsinfo, x(:), this%force_uniform_grid ) + case ( MESHTYPE_2D_RECTDOM ) ! 2D mesh + call File_common_meshfield_get_axis( this%mesh2D, this%dimsinfo, x(:), y(:), this%force_uniform_grid ) + case ( MESHTYPE_2D_CUBEDSPHERE ) ! 2D mesh + call File_common_meshfield_get_axis( this%meshCS2D, this%dimsinfo, x(:), y(:) ) + case ( MESHTYPE_3D_CUBEDOM ) ! 3D mesh + call File_common_meshfield_get_axis( this%mesh3D, this%dimsinfo, x(:), y(:), z(:), this%force_uniform_grid ) + case ( MESHTYPE_3D_CUBEDSPHERE ) ! 3D mesh + call File_common_meshfield_get_axis( this%meshCS3D, this%dimsinfo, x(:), y(:), z(:) ) + end select + + !- + select case ( this%mesh_type_id ) + case ( MESHTYPE_1D ) ! 1D mesh + call FILE_Write_Axis( fid, this%dimsinfo(1)%name, x(:), start(1:1) ) + case ( MESHTYPE_2D_RECTDOM, MESHTYPE_2D_CUBEDSPHERE ) ! 2D mesh call FILE_Write_Axis( fid, this%dimsinfo(1)%name, x(:), start(1:1) ) call FILE_Write_Axis( fid, this%dimsinfo(2)%name, y(:), start(2:2) ) case ( MESHTYPE_3D_CUBEDOM, MESHTYPE_3D_CUBEDSPHERE ) ! 3D mesh - allocate( x(this%dimsinfo(1)%size), y(this%dimsinfo(2)%size), z(this%dimsinfo(3)%size) ) - call File_common_meshfield_get_axis( this%mesh3D, this%dimsinfo, x(:), y(:), z(:), force_uniform_grid ) - call FILE_Write_Axis( fid, this%dimsinfo(1)%name, x(:), start(1:1) ) call FILE_Write_Axis( fid, this%dimsinfo(2)%name, y(:), start(2:2) ) call FILE_Write_Axis( fid, this%dimsinfo(3)%name, z(:), start(3:3) ) diff --git a/FElib/src/fluid_dyn_solver/scale_atm_dyn_dgm_hydrostatic.F90 b/FElib/src/fluid_dyn_solver/scale_atm_dyn_dgm_hydrostatic.F90 index f8cdcc19..a24b3bfe 100644 --- a/FElib/src/fluid_dyn_solver/scale_atm_dyn_dgm_hydrostatic.F90 +++ b/FElib/src/fluid_dyn_solver/scale_atm_dyn_dgm_hydrostatic.F90 @@ -320,17 +320,23 @@ subroutine hydrostatic_calc_basicstate_constPTLAPS( & CPovR = CpDry / Rdry exner_sfc = (PRES_sfc / PRES00)**RovCP - !$omp parallel do private(PT, exner) - do ke=lcmesh3D%NeS, lcmesh3D%NeE - ! d exner / dz = - g / ( Cp * PT0 ) / (1 + PTLAPS/PT0 * z) - ! exner = exner(zs) - g / (Cp * PTLAPS ) * log[ 1 + PTLAPS/PT0 * z ] - PT(:) = PotTemp0 + PTLAPS * z(:,ke) - exner(:) = exner_sfc - Grav / ( CpDry * PTLAPS ) * log( 1.0_RP + PTLAPS / PotTemp0 * z(:,ke) ) - - PRES_hyd(:,ke) = PRES00 * exner(:)**CPovR - DENS_hyd(:,ke) = PRES_hyd(:,ke) / ( Rdry * exner(:) * PT(:) ) - end do - + if ( PTLAPS == 0.0_RP ) then + call hydrostatic_calc_basicstate_constPT( & + DENS_hyd, PRES_hyd, & + PotTemp0, PRES_sfc, x, y, z, lcmesh3D, elem ) + else + !$omp parallel do private(PT, exner) + do ke=lcmesh3D%NeS, lcmesh3D%NeE + ! d exner / dz = - g / ( Cp * PT0 ) / (1 + PTLAPS/PT0 * z) + ! exner = exner(zs) - g / (Cp * PTLAPS ) * log[ 1 + PTLAPS/PT0 * z ] + PT(:) = PotTemp0 + PTLAPS * z(:,ke) + exner(:) = exner_sfc - Grav / ( CpDry * PTLAPS ) * log( 1.0_RP + PTLAPS / PotTemp0 * z(:,ke) ) + + PRES_hyd(:,ke) = PRES00 * exner(:)**CPovR + DENS_hyd(:,ke) = PRES_hyd(:,ke) / ( Rdry * exner(:) * PT(:) ) + end do + end if + return end subroutine hydrostatic_calc_basicstate_constPTLAPS diff --git a/model/atm_nonhydro3d/src/admin/mod_dg_driver.F90 b/model/atm_nonhydro3d/src/admin/mod_dg_driver.F90 index ca8e06c1..22ef3c06 100644 --- a/model/atm_nonhydro3d/src/admin/mod_dg_driver.F90 +++ b/model/atm_nonhydro3d/src/admin/mod_dg_driver.F90 @@ -175,6 +175,7 @@ subroutine dg_driver( & if ( atmos%phy_mp_proc%IsActivated() ) call atmos%phy_mp_proc%vars%History() if ( atmos%phy_sfc_proc%IsActivated() ) call atmos%phy_sfc_proc%vars%History() if ( atmos%phy_rd_proc%IsActivated() ) call atmos%phy_rd_proc%vars%History() + if ( atmos%phy_bl_proc%IsActivated() ) call atmos%phy_bl_proc%vars%History() if ( ocean%IsActivated() ) call ocean%vars%History() @@ -359,14 +360,11 @@ subroutine restart_read() if ( atmos%isActivated() ) then call atmos%vars%History() - if ( atmos%phy_sfc_proc%IsActivated() ) & - call atmos%phy_sfc_proc%vars%History() - if ( atmos%phy_tb_proc%IsActivated() ) & - call atmos%phy_tb_proc%vars%History() - if ( atmos%phy_mp_proc%IsActivated() ) & - call atmos%phy_mp_proc%vars%History() - if ( atmos%phy_rd_proc%IsActivated() ) & - call atmos%phy_rd_proc%vars%History() + if ( atmos%phy_sfc_proc%IsActivated() ) call atmos%phy_sfc_proc%vars%History() + if ( atmos%phy_tb_proc%IsActivated() ) call atmos%phy_tb_proc%vars%History() + if ( atmos%phy_mp_proc%IsActivated() ) call atmos%phy_mp_proc%vars%History() + if ( atmos%phy_rd_proc%IsActivated() ) call atmos%phy_rd_proc%vars%History() + if ( atmos%phy_bl_proc%IsActivated() ) call atmos%phy_bl_proc%vars%History() call atmos%vars%Monitor() end if diff --git a/model/atm_nonhydro3d/src/atmos/mod_atmos_component.F90 b/model/atm_nonhydro3d/src/atmos/mod_atmos_component.F90 index df90274c..b199c443 100644 --- a/model/atm_nonhydro3d/src/atmos/mod_atmos_component.F90 +++ b/model/atm_nonhydro3d/src/atmos/mod_atmos_component.F90 @@ -250,6 +250,7 @@ subroutine Atmos_setup( this ) !- Setup the module for atmosphere / physics / PBL turbulence parameterization call this%phy_bl_proc%ModelComponentProc_Init( 'AtmosPhysBl', ATMOS_PHY_BL_DO ) call this%phy_bl_proc%setup( this%mesh, this%time_manager ) + call this%phy_bl_proc%SetDynBC( this%dyn_proc%dyncore_driver%boundary_cond ) !- Setup diff --git a/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_bl.F90 b/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_bl.F90 index 6e90e273..c711acb6 100644 --- a/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_bl.F90 +++ b/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_bl.F90 @@ -21,6 +21,7 @@ module mod_atmos_phy_bl use scale_const, only: & UNDEF8 => CONST_UNDEF8 + use scale_element_line, only: LineElement use scale_mesh_base, only: MeshBase use scale_mesh_base2d, only: MeshBase2D use scale_mesh_base3d, only: MeshBase3D @@ -40,10 +41,12 @@ module mod_atmos_phy_bl use scale_model_var_manager, only: ModelVarManager use scale_model_component_proc, only: ModelComponentProc - use mod_atmos_phy_bl_vars, only: AtmosPhyBlVars + use scale_atm_dyn_dgm_bnd, only: AtmDynBnd + use mod_atmos_phy_bl_vars, only: AtmosPhyBlVars use mod_atmos_vars_container, only: & AtmosVarsContainer + !----------------------------------------------------------------------------- implicit none @@ -61,11 +64,18 @@ module mod_atmos_phy_bl integer :: atm_var_container_typeid !< Type ID of variable container for PBL turbulence parameterization + real(RP) :: dtsec !< Timestep for PBL turbulence parameterization + + type(LineElement) :: v_elem1D + type(AtmDynBnd), pointer :: dyn_bnd !< Pointer to object for treating boundary conditions with atmospheric dynamics + real(RP) :: C_IP !< Parameter for symmetric interior penalty method in DGM + logical :: use_delta_form !< Flag to use delta form in the vertical implicit time integration of PBL scheme contains procedure :: setup => AtmosPhyBl_setup procedure :: calc_tendency => AtmosPhyBl_calc_tendency procedure :: update => AtmosPhyBl_update procedure :: finalize => AtmosPhyBl_finalize + procedure, public :: SetDynBC => AtmosPhyBl_SetDynBC end type AtmosPhyBl !----------------------------------------------------------------------------- @@ -79,6 +89,8 @@ module mod_atmos_phy_bl ! !++ Private parameters & variables ! + + integer, parameter :: BL_TYPEID_MYNN_LEVEL2 = 1 !< Type ID of MYNN Level 2 PBL scheme contains !> Setup a component of planetary boundary layer (PBL) turbulence parameterization in atmospheric model @@ -87,8 +99,12 @@ module mod_atmos_phy_bl !! @param tm_parent_comp Object to mange a temporal scheme in a parent component !! subroutine AtmosPhyBl_setup( this, model_mesh, tm_parent_comp ) + use scale_tracer, only: QA use mod_atmos_mesh, only: AtmosMesh use scale_time_manager, only: TIME_manager_component + + use scale_atm_phy_bl_dgm_mynn_lv2, only: & + atm_phy_bl_dgm_mynn_lv2_Init use mod_atmos_vars, only: ATM_VARS_CONTAINER_PRIMARY_ID implicit none class(AtmosPhyBl), intent(inout) :: this @@ -98,14 +114,18 @@ subroutine AtmosPhyBl_setup( this, model_mesh, tm_parent_comp ) real(DP) :: TIME_DT = UNDEF8 !< Timestep for PBL turbulence parameterization character(len=H_SHORT) :: TIME_DT_UNIT = 'SEC' !< Unit of timestep - character(len=H_MID) :: BL_TYPE = 'NONE' !< Type of a PBL turbulence parameterization scheme - integer :: atm_var_container_typeid + character(len=H_MID) :: BL_TYPE = 'NONE' !< Type of a PBL turbulence parameterization scheme + integer :: atm_var_container_typeid + real(RP) :: C_IP = 1.0_RP !< Parameter for symmetric interior penalty method in DGM + logical :: use_delta_form = .false. !< Flag to use delta form in the vertical implicit time integration of PBL scheme namelist /PARAM_ATMOS_PHY_BL/ & - TIME_DT, & - TIME_DT_UNIT, & - BL_TYPE, & - atm_var_container_typeid + TIME_DT, & + TIME_DT_UNIT, & + BL_TYPE, & + atm_var_container_typeid, & + C_IP, & + use_delta_form class(AtmosMesh), pointer :: atm_mesh class(MeshBase), pointer :: ptr_mesh @@ -136,6 +156,8 @@ subroutine AtmosPhyBl_setup( this, model_mesh, tm_parent_comp ) LOG_NML(PARAM_ATMOS_PHY_BL) this%atm_var_container_typeid = atm_var_container_typeid + this%C_IP = C_IP + this%use_delta_form = use_delta_form !- Get atmospheric mesh -------------------------------------------------- @@ -149,10 +171,18 @@ subroutine AtmosPhyBl_setup( this, model_mesh, tm_parent_comp ) call tm_parent_comp%Regist_process( 'ATMOS_PHY_BL', TIME_DT, TIME_DT_UNIT, & ! (in) this%tm_process_id ) ! (out) + + this%dtsec = tm_parent_comp%process_list(this%tm_process_id)%dtsec !--- Set the type of PBL turbulence parameterization select case( BL_TYPE ) + case( 'MYNN_LEVEL2' ) + this%BL_TYPEID = BL_TYPEID_MYNN_LEVEL2 + QS_BL = QA + QA_BL = 0 + + call atm_phy_bl_dgm_mynn_lv2_Init( atm_mesh%ptr_mesh ) case default LOG_ERROR("ATMOS_PHY_BL_setup",*) 'Not appropriate PBL turbulence parameterization type. Check!' call PRC_abort @@ -163,6 +193,9 @@ subroutine AtmosPhyBl_setup( this, model_mesh, tm_parent_comp ) !- Initialize the variables call this%vars%Init( model_mesh, QS_BL, QE_BL, QA_BL ) + !- + call this%v_elem1D%Init( atm_mesh%ptr_mesh%refElem3D%PolyOrder_v, .false. ) + return end subroutine AtmosPhyBl_setup @@ -179,6 +212,26 @@ end subroutine AtmosPhyBl_setup subroutine AtmosPhyBl_calc_tendency( & this, model_mesh, prgvars_list, trcvars_list, & auxvars_list, forcing_list, is_update ) + use scale_tracer, only: QA + use scale_atm_phy_bl_dgm_mynn_lv2, only: & + atm_phy_bl_dgm_mynn_lv2_cal_VViscDiffCoef + use scale_atm_phy_bl_dgm_common, only: & + atm_phy_bl_dgm_common_calc_tendency + + use mod_atmos_vars, only: & + AtmosVars_GetLocalMeshPrgVars, & + AtmosVars_GetLocalMeshPhyAuxVars, & + AtmosVars_GetLocalMeshQTRCVarList, & + AtmosVars_GetLocalMeshPhyTends + use mod_atmos_phy_bl_vars, only: & + AtmosPhyBLVars_GetLocalMeshFields_tend, & + RHOU_tp_ID => ATMOS_PHY_BL_RHOU_t_ID, & + RHOV_tp_ID => ATMOS_PHY_BL_RHOV_t_ID, & + RHOT_tp_ID => ATMOS_PHY_BL_RHOT_t_ID, & + TKE_ID => ATMOS_PHY_BL_DIAG_TKE_ID, & + NU_ID => ATMOS_PHY_BL_DIAG_NU_ID, & + KH_ID => ATMOS_PHY_BL_DIAG_KH_ID + implicit none class(AtmosPhyBl), intent(inout) :: this class(ModelMeshBase), intent(in) :: model_mesh @@ -187,7 +240,121 @@ subroutine AtmosPhyBl_calc_tendency( & class(ModelVarManager), intent(inout) :: auxvars_list class(ModelVarManager), intent(inout) :: forcing_list logical, intent(in) :: is_update + + class(MeshBase), pointer :: mesh + class(MeshBase3D), pointer :: mesh3D + class(LocalMesh3D), pointer :: lcmesh + + integer :: n + integer :: ke + integer :: iq + + class(LocalMeshFieldBase), pointer :: DDENS, MOMX, MOMY, MOMZ, DRHOT + class(LocalMeshFieldBase), pointer :: DENS_hyd, PRES_hyd, Rtot, CVtot, CPtot + class(LocalMeshFieldBase), pointer :: PRES, PT + + class(LocalMeshFieldBase), pointer :: DENS_tp, MOMX_tp, MOMY_tp, MOMZ_tp, RHOT_tp, RHOH_P + type(LocalMeshFieldBaseList) :: RHOQ_tp(QA) + class(LocalMeshFieldBase), pointer :: bl_RHOU_t, bl_RHOV_t, bl_RHOT_t + + type DYN_BNDInfo + logical, allocatable :: is_bound(:,:) + end type + type(DYN_BNDInfo), allocatable :: bnd_info(:) !------------------------------------------------------------------------ + + if (.not. this%IsActivated()) return + + LOG_PROGRESS(*) 'atmosphere / physics / planetary boundary layer' + + call model_mesh%GetModelMesh( mesh ) + select type(mesh) + class is (MeshBase3D) + mesh3D => mesh + end select + + !- + if ( is_update ) then + call PROF_rapstart( 'ATM_BL_tendency', 2) + + allocate( bnd_info(mesh3D%LOCAL_MESH_NUM) ) + + do n=1, mesh3D%LOCAL_MESH_NUM + call AtmosVars_GetLocalMeshPrgVars( n, & + mesh, prgvars_list, auxvars_list, & + DDENS, MOMX, MOMY, MOMZ, DRHOT, & + DENS_hyd, PRES_hyd, Rtot, CVtot, CPtot, & + lcmesh ) + + call AtmosVars_GetLocalMeshPhyAuxVars( n, & + mesh, auxvars_list, & + PRES, PT ) + + call AtmosPhyBLVars_GetLocalMeshFields_tend( n, & + mesh, this%vars%tends_manager, & + bl_RHOU_t, bl_RHOV_t, bl_RHOT_t ) + + !- + allocate( bnd_info(n)%is_bound(lcmesh%refElem3D%NfpTot,lcmesh%Ne) ) + call this%dyn_bnd%Inquire_bound_flag( bnd_info(n)%is_bound, & ! (out) + n, lcmesh%VMapM, lcmesh%VMapP, lcmesh%VMapB, & ! (in) + lcmesh, lcmesh%refElem3D ) ! (in) + + select case( this%BL_TYPEID ) + case( BL_TYPEID_MYNN_LEVEL2 ) + call atm_phy_bl_dgm_mynn_lv2_cal_VViscDiffCoef( & + this%vars%diagvars(NU_ID)%local(n)%val, & ! (out) + this%vars%diagvars(KH_ID)%local(n)%val, & ! (out) + this%vars%diagvars(TKE_ID)%local(n)%val, & ! (out) + DDENS%val, MOMX%val, MOMY%val, MOMZ%val, DRHOT%val, & ! (in) + DENS_hyd%val, PRES_hyd%val, PRES%val, PT%val, & ! (in) + model_mesh%DOptrMat(3), model_mesh%LiftOptrMat, & ! (in) + lcmesh, lcmesh%refElem3D, bnd_info(n)%is_bound ) ! (in) + end select + + call atm_phy_bl_dgm_common_calc_tendency( & + bl_RHOU_t%val, bl_RHOV_t%val, bl_RHOT_t%val, & ! (out) + DDENS%val, MOMX%val, MOMY%val, DRHOT%val, & ! (in) + PT%val, DENS_hyd%val, PRES_hyd%val, & ! (in) + this%vars%diagvars(NU_ID)%local(n)%val, & ! (in) + this%vars%diagvars(KH_ID)%local(n)%val, & ! (in) + model_mesh%element3D_operation, & ! (in) + this%C_IP, this%dtsec, & ! (in) + lcmesh, lcmesh%refElem3D, this%v_elem1D, & ! (in) + bnd_info(n)%is_bound, this%use_delta_form ) ! (in) + + end do + + do n=1, mesh3D%LOCAL_MESH_NUM + deallocate( bnd_info(n)%is_bound ) + end do + call PROF_rapend( 'ATM_BL_tendency', 2) + end if + + call PROF_rapstart('ATM_PHY_BL_add_tend', 2) + do n=1, mesh%LOCAL_MESH_NUM + call AtmosVars_GetLocalMeshPhyTends( n, & + mesh, forcing_list, & + DENS_tp, MOMX_tp, MOMY_tp, MOMZ_tp, RHOT_tp, & + RHOH_p, RHOQ_tp ) + + call AtmosPhyBLVars_GetLocalMeshFields_tend( n, & + mesh, this%vars%tends_manager, & + bl_RHOU_t, bl_RHOV_t, bl_RHOT_t, & + lcmesh ) + + !$omp parallel private(ke, iq) + !$omp do + do ke=lcmesh%NeS, lcmesh%NeE + MOMX_tp%val(:,ke) = MOMX_tp%val(:,ke) + bl_RHOU_t%val(:,ke) + MOMY_tp%val(:,ke) = MOMY_tp%val(:,ke) + bl_RHOV_t%val(:,ke) + RHOT_tp%val(:,ke) = RHOT_tp%val(:,ke) + bl_RHOT_t%val(:,ke) + end do + !$omp end parallel + end do + call PROF_rapend('ATM_PHY_BL_add_tend', 2) + + return end subroutine AtmosPhyBl_calc_tendency @@ -220,20 +387,39 @@ end subroutine AtmosPhyBl_update !! !OCL SERIAL subroutine AtmosPhyBl_finalize( this ) + use scale_atm_phy_bl_dgm_mynn_lv2, only: & + atm_phy_bl_dgm_mynn_lv2_Final implicit none class(AtmosPhyBl), intent(inout) :: this !-------------------------------------------------- if (.not. this%IsActivated()) return - ! select case ( this%BL_TYPEID ) - ! case( BL_TYPEID_LSCOND ) - ! end select + select case ( this%BL_TYPEID ) + case( BL_TYPEID_MYNN_LEVEL2 ) + call atm_phy_bl_dgm_mynn_lv2_Final() + end select call this%vars%Final() + call this%v_elem1D%Final() return end subroutine AtmosPhyBl_finalize +!> Set boundary conditions to PBL component in atmospheric model +!! +!! @param dyn_bnd Object to manage boundary conditions of dynamical core +!! +!OCL SERIAL + subroutine AtmosPhyBl_setDynBC( this, dyn_bnd ) + implicit none + class(AtmosPhyBl), intent(inout) :: this + type(AtmDynBnd), intent(in), target :: dyn_bnd + !-------------------------------------------------- + + this%dyn_bnd => dyn_bnd + + return + end subroutine AtmosPhyBl_setDynBC !- private ------------------------------------------------ end module mod_atmos_phy_bl diff --git a/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_bl_vars.F90 b/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_bl_vars.F90 index e55b869b..e7e97be2 100644 --- a/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_bl_vars.F90 +++ b/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_bl_vars.F90 @@ -56,6 +56,9 @@ module mod_atmos_phy_bl_vars type(MeshField3D), allocatable :: tends(:) !< Array of tendency variables type(ModelVarManager) :: tends_manager !< Object to manage tendencies + type(MeshField3D), allocatable :: diagvars(:) + type(ModelVarManager) :: diagvars_manager + integer :: QS !< Start index of tracer variables with PBL turbulence parameterization integer :: QE !< End index of tracer variables with PBL turbulence parameterization integer :: QA !< Number of tracer variables with PBL turbulence parameterization @@ -67,6 +70,8 @@ module mod_atmos_phy_bl_vars procedure :: History => AtmosPhyBlVars_history end type AtmosPhyBlVars + public :: AtmosPhyBlVars_GetLocalMeshFields_tend + !----------------------------------------------------------------------------- ! !++ Public variables @@ -85,6 +90,20 @@ module mod_atmos_phy_bl_vars VariableInfo( ATMOS_PHY_BL_RHOT_t_ID, 'BL_RHOT_t', 'tendency of rho*PT in BL process', & 'kg/m3.K/s', 3, 'XYZ', '' ) / + integer, public, parameter :: ATMOS_PHY_BL_DIAG_TKE_ID = 1 + integer, public, parameter :: ATMOS_PHY_BL_DIAG_NU_ID = 2 + integer, public, parameter :: ATMOS_PHY_BL_DIAG_KH_ID = 3 + integer, public, parameter :: ATMOS_PHY_BL_DIAG_NUM = 3 + + type(VariableInfo) :: ATMOS_PHY_BL_DIAG_VINFO(ATMOS_PHY_BL_DIAG_NUM) + DATA ATMOS_PHY_BL_DIAG_VINFO / & + VariableInfo( ATMOS_PHY_BL_DIAG_TKE_ID, 'TKE', 'SGS turbulence kinetic energy', & + 'm2/s2', 3, 'XYZ', '' ), & + VariableInfo( ATMOS_PHY_BL_DIAG_NU_ID, 'NU', 'eddy viscosity', & + 'm2/s', 3, 'XYZ', '' ), & + VariableInfo( ATMOS_PHY_BL_DIAG_KH_ID, 'KH', 'eddy diffusion', & + 'm2/s', 3, 'XYZ', '' ) / + !----------------------------------------------------------------------------- ! !++ Private procedures @@ -138,20 +157,16 @@ subroutine AtmosPhyBlVars_Init( this, model_mesh, & call mesh3D%GetMesh2D( mesh2D ) - !---- + !- Initialize tendency variables call this%tends_manager%Init() allocate( this%tends(this%TENDS_NUM_TOT) ) reg_file_hist = .true. - do iv = 1, ATMOS_PHY_BL_TENDS_NUM1 - call this%tends_manager%Regist( & - ATMOS_PHY_BL_TEND_VINFO(iv), mesh3D, & - this%tends(iv), reg_file_hist ) - - do n = 1, mesh3D%LOCAL_MESH_NUM - this%tends(iv)%local(n)%val(:,:) = 0.0_RP - end do + do iv=1, ATMOS_PHY_BL_TENDS_NUM1 + call this%tends_manager%Regist( & + ATMOS_PHY_BL_TEND_VINFO(iv), mesh3D, & + this%tends(iv), reg_file_hist, fill_zero=.true. ) end do qtrc_tp_vinfo_tmp%ndims = 3 @@ -167,13 +182,21 @@ subroutine AtmosPhyBlVars_Init( this, model_mesh, & reg_file_hist = .true. call this%tends_manager%Regist( & - qtrc_tp_vinfo_tmp, mesh3D, & - this%tends(iv), reg_file_hist ) - - do n = 1, mesh3D%LOCAL_MESH_NUM - this%tends(iv)%local(n)%val(:,:) = 0.0_RP - end do - end do + qtrc_tp_vinfo_tmp, mesh3D, & + this%tends(iv), reg_file_hist, fill_zero=.true. ) + end do + + !- Initialize diagnostic variables + + call this%diagvars_manager%Init() + allocate( this%diagvars(ATMOS_PHY_BL_DIAG_NUM) ) + + reg_file_hist = .true. + do iv=1, ATMOS_PHY_BL_DIAG_NUM + call this%diagvars_manager%Regist( & + ATMOS_PHY_BL_DIAG_VINFO(iv), mesh3D, & + this%diagvars(iv), reg_file_hist, fill_zero=.true. ) + end do return end subroutine AtmosPhyBlVars_Init @@ -184,15 +207,85 @@ subroutine AtmosPhyBlVars_Final( this ) implicit none class(AtmosPhyBlVars), intent(inout) :: this !---------------------------------------------------- + + LOG_INFO('AtmosPhyBlVars_Final',*) + + call this%tends_manager%Final() + deallocate( this%tends ) + + call this%diagvars_manager%Final() + deallocate( this%diagvars ) return end subroutine AtmosPhyBlVars_Final +!OCL SERIAL + subroutine AtmosPhyBlVars_GetLocalMeshFields_tend( domID, mesh, bl_tends_list, & + bl_RHOU_t, bl_RHOV_t, bl_RHOT_t, & + lcmesh3D & + ) + + use scale_mesh_base, only: MeshBase + use scale_meshfield_base, only: MeshFieldBase + implicit none + + integer, intent(in) :: domID + class(MeshBase), intent(in) :: mesh + class(ModelVarManager), intent(inout) :: bl_tends_list + class(LocalMeshFieldBase), pointer, intent(out) :: bl_RHOU_t + class(LocalMeshFieldBase), pointer, intent(out) :: bl_RHOV_t + class(LocalMeshFieldBase), pointer, intent(out) :: bl_RHOT_t + class(LocalMesh3D), pointer, intent(out), optional :: lcmesh3D + + class(MeshFieldBase), pointer :: field + class(LocalMeshBase), pointer :: lcmesh + + integer :: iq + !------------------------------------------------------- + + !-- + call bl_tends_list%Get(ATMOS_PHY_BL_RHOU_t_ID, field) + call field%GetLocalMeshField(domID, bl_RHOU_t) + + call bl_tends_list%Get(ATMOS_PHY_BL_RHOV_t_ID, field) + call field%GetLocalMeshField(domID, bl_RHOV_t) + + call bl_tends_list%Get(ATMOS_PHY_BL_RHOT_t_ID, field) + call field%GetLocalMeshField(domID, bl_RHOT_t) + + if (present(lcmesh3D)) then + call mesh%GetLocalMesh( domID, lcmesh ) + nullify( lcmesh3D ) + + select type(lcmesh) + type is (LocalMesh3D) + if (present(lcmesh3D)) lcmesh3D => lcmesh + end select + end if + + return + end subroutine AtmosPhyBlVars_GetLocalMeshFields_tend + + !> Put data with BL variables to history file !OCL SERIAL subroutine AtmosPhyBlVars_history( this ) use scale_file_history_meshfield, only: FILE_HISTORY_meshfield_put implicit none class(AtmosPhyBlVars), intent(inout) :: this + + integer :: v + integer :: hst_id !---------------------------------------------------- + + do v=1, this%TENDS_NUM_TOT + hst_id = this%tends(v)%hist_id + if ( hst_id > 0 ) call FILE_HISTORY_meshfield_put( hst_id, this%tends(v) ) + end do + + do v=1, ATMOS_PHY_BL_DIAG_NUM + hst_id = this%diagvars(v)%hist_id + if ( hst_id > 0 ) call FILE_HISTORY_meshfield_put( hst_id, this%diagvars(v) ) + end do + return end subroutine AtmosPhyBlVars_history diff --git a/model/atm_nonhydro3d/test/case/boundary_layer/Makefile b/model/atm_nonhydro3d/test/case/boundary_layer/Makefile new file mode 100644 index 00000000..a02773a1 --- /dev/null +++ b/model/atm_nonhydro3d/test/case/boundary_layer/Makefile @@ -0,0 +1,28 @@ +################################################################################ +# +# Makefile for each test program +# +################################################################################ + +PWD = $(shell pwd) +TOPDIR = $(abspath ../../../../..) +TESTDIR = ../.. + + +# user-defined source files +CODE_DIR = . +ORG_SRCS = mod_user.F90 + +# parameters for run +INITCONF = init.conf +RUNCONF = run.conf +TPROC = 1 +OMP_NUM_THREADS = 4 + + +# required data (parameters,distributed files) +DATPARAM = +DATDISTS = + +# build, makedir, run, jobshell, allclean, clean is inside of common Makefile +include $(TESTDIR)/Makefile.common diff --git a/model/atm_nonhydro3d/test/case/boundary_layer/init.conf b/model/atm_nonhydro3d/test/case/boundary_layer/init.conf new file mode 100644 index 00000000..ca90c4b3 --- /dev/null +++ b/model/atm_nonhydro3d/test/case/boundary_layer/init.conf @@ -0,0 +1,57 @@ +#--- Configuration file for a test case of boundary layer ------- +&PARAM_IO + IO_LOG_BASENAME = 'init_LOG', +/ +&PARAM_MKINIT + initname = 'boundary_layer', +/ +&PARAM_RESTART + OUTPUT_FLAG = .true., + OUT_BASENAME = 'init' +/ +&PARAM_TIME + TIME_STARTDATE = 0000, 1, 1, 0, 0, 0, + TIME_STARTMS = 0.D0, +/ +&PARAM_CONST +/ +&PARAM_EXP + U0 = 10.0D0, + TEMP0 = 300.0D0, + ENV_RH = 0D0, + NITER_RH = 3, +/ +#** ATMOS ****************************************************** +&PARAM_ATMOS + ACTIVATE_FLAG = .true., + ATMOS_DYN_DO = .true., + ATMOS_PHY_BL_DO = .true., +! ATMOS_PHY_SF_DO = .true., + ATMOS_USE_QV = .true., +/ +&PARAM_ATMOS_MESH + dom_xmin = -5.0D3, + dom_xmax = 5.0D3, + isPeriodicX = .true., + dom_ymin = -5.0D3, + dom_ymax = 5.0D3, + isPeriodicY = .true., + dom_zmin = 0.0D0, + dom_zmax = 1.0D3, + NprcX = 1, + NeX = 1, + NprcY = 1, + NeY = 1, + NeZ = 8, + PolyOrder_h = 7, + PolyOrder_v = 7, +! LumpedMassMatFlag = .true., +/ +#** ATMOS / DYN ****************************************************** +&PARAM_ATMOS_DYN + EQS_TYPE = "NONHYDRO3D_HEVE", + TINTEG_TYPE = 'ERK_SSP_3s3o', +/ +&PARAM_ATMOS_PHY_BL + BL_TYPE = 'MYNN_LEVEL2', +/ \ No newline at end of file diff --git a/model/atm_nonhydro3d/test/case/boundary_layer/mod_user.F90 b/model/atm_nonhydro3d/test/case/boundary_layer/mod_user.F90 new file mode 100644 index 00000000..ae8b864f --- /dev/null +++ b/model/atm_nonhydro3d/test/case/boundary_layer/mod_user.F90 @@ -0,0 +1,321 @@ +!------------------------------------------------------------------------------- +!> module USER +!! +!! @par Description +!! User defined module for a test case of boundary layer scheme +!! +!! @author Yuta Kawai, Team SCALE +!! +!< +!------------------------------------------------------------------------------- +#include "scaleFElib.h" +module mod_user + + !----------------------------------------------------------------------------- + ! + !++ used modules + ! + use scale_precision + use scale_io + use scale_prof + use scale_prc, only: PRC_abort + + use mod_atmos_component, only: & + AtmosComponent + + use scale_element_base, only: ElementBase3D + use scale_element_hexahedral, only: HexahedralElement + use scale_localmesh_3d, only: LocalMesh3D + use scale_meshfield_base, only: MeshField3D + + use mod_user_base, only: UserBase + use mod_experiment, only: Experiment + + !----------------------------------------------------------------------------- + implicit none + private + !----------------------------------------------------------------------------- + ! + !++ Public type & procedure + ! + type, public, extends(UserBase) :: User + contains + procedure :: mkinit_ => USER_mkinit + generic :: mkinit => mkinit_ + procedure :: setup_ => USER_setup + generic :: setup => setup_ + procedure :: calc_tendency => USER_calc_tendency + end type User + + !----------------------------------------------------------------------------- + ! + !++ Public parameters & variables + ! + !----------------------------------------------------------------------------- + ! + !++ Private procedure + ! + !----------------------------------------------------------------------------- + ! + !++ Private parameters & variables + ! + + type(MeshField3D), private :: PRES_diff + + !----------------------------------------------------------------------------- +contains +!OCL SERIAL + subroutine USER_mkinit ( this, atm ) + implicit none + class(User), intent(inout) :: this + class(AtmosComponent), intent(inout) :: atm + + type(Experiment) :: exp_manager + !------------------------------------------ + + call exp_manager%Init('boundary_layer') + call exp_manager%Regist_SetInitCond( exp_SetInitCond_boundary_layer ) + call this%UserBase%mkinit( atm, exp_manager ) + call exp_manager%Final() + + return + end subroutine USER_mkinit + +!OCL SERIAL + subroutine USER_setup( this, atm ) + implicit none + class(User), intent(inout) :: this + class(AtmosComponent), intent(inout) :: atm + + logical :: USER_do = .false. !< do user + namelist / PARAM_USER / & + USER_do + + integer :: ierr + !------------------------------------------ + + + LOG_NEWLINE + LOG_INFO("USER_setup",*) 'Setup' + + !--- read namelist + rewind(IO_FID_CONF) + read(IO_FID_CONF,nml=PARAM_USER,iostat=ierr) + if( ierr < 0 ) then !--- missing + LOG_INFO("USER_setup",*) 'Not found namelist. Default used.' + elseif( ierr > 0 ) then !--- fatal error + LOG_ERROR("USER_setup",*) 'Not appropriate names in namelist PARAM_USER. Check!' + call PRC_abort + endif + LOG_NML(PARAM_USER) + + call this%UserBase%Setup( atm, USER_do ) + + !- + if ( USER_do ) call PRES_diff%Init( 'PRES_diff', 'Pa', atm%mesh%ptr_mesh ) + + return + end subroutine USER_setup + +!OCL SERIAL + subroutine USER_calc_tendency( this, atm ) + use scale_file_history_meshfield, only: & + FILE_HISTORY_meshfield_in + implicit none + + class(User), intent(inout) :: this + class(AtmosComponent), intent(inout) :: atm + !------------------------------------------ + + if ( this%USER_do ) then + call atm%vars%Calc_diagVar( 'PRES_diff', PRES_diff ) + call FILE_HISTORY_meshfield_in( PRES_diff, "perturbation of PRES" ) + end if + + return + end subroutine USER_calc_tendency + + !------ + +!OCL SERIAL + subroutine exp_SetInitCond_boundary_layer( this, & + DENS_hyd, PRES_hyd, DDENS, MOMX, MOMY, MOMZ, DRHOT, tracer_field_list, & + x, y, z, dom_xmin, dom_xmax, dom_ymin, dom_ymax, dom_zmin, dom_zmax, & + lcmesh, elem ) + use scale_tracer, only: & + TRACER_inq_id + use scale_const, only: & + PI => CONST_PI, & + GRAV => CONST_GRAV, & + Rdry => CONST_Rdry, & + Rvap => CONST_Rvap, & + CPdry => CONST_CPdry, & + CVdry => CONST_CVdry, & + PRES00 => CONST_PRE00 + use scale_atmos_saturation, only: & + ATMOS_SATURATION_psat_all, & + ATMOS_SATURATION_pres2qsat_all + use scale_atmos_hydrometeor, only: & + CV_VAPOR, & + CV_WATER, & + CP_VAPOR, & + CP_WATER + use scale_atm_dyn_dgm_hydrostatic, only: & + hydrostatic_calc_basicstate_constPTLAPS, & + hydrostatic_build_rho_XYZ + use mod_experiment, only: & + TracerLocalMeshField_ptr + implicit none + + class(Experiment), intent(inout) :: this + type(LocalMesh3D), intent(in) :: lcmesh + class(ElementBase3D), intent(in) :: elem + real(RP), intent(out) :: DENS_hyd(elem%Np,lcmesh%NeA) + real(RP), intent(out) :: PRES_hyd(elem%Np,lcmesh%NeA) + real(RP), intent(out) :: DDENS(elem%Np,lcmesh%NeA) + real(RP), intent(out) :: MOMX(elem%Np,lcmesh%NeA) + real(RP), intent(out) :: MOMY(elem%Np,lcmesh%NeA) + real(RP), intent(out) :: MOMZ(elem%Np,lcmesh%NeA) + real(RP), intent(out) :: DRHOT(elem%Np,lcmesh%NeA) + type(TracerLocalMeshField_ptr), intent(inout) :: tracer_field_list(:) + real(RP), intent(in) :: x(elem%Np,lcmesh%Ne) + real(RP), intent(in) :: y(elem%Np,lcmesh%Ne) + real(RP), intent(in) :: z(elem%Np,lcmesh%Ne) + real(RP), intent(in) :: dom_xmin, dom_xmax + real(RP), intent(in) :: dom_ymin, dom_ymax + real(RP), intent(in) :: dom_zmin, dom_zmax + + real(RP) :: U0 = 0.0_RP !< Initial horizontal velocity [m/s] + real(RP) :: TEMP0 = 250.0_RP !< Initial temperature [K] + real(RP) :: SFC_PRES = 1.0E5_RP !< Surface pressure [Pa] + real(RP) :: PTLAPS = 0.0_RP !< Potential temperature lapse rate [K/m] + real(RP) :: VSHEAR_HVEL = 0.01_RP !< Vertical shear of horizontal velocity [m/s/m] + + real(RP) :: ENV_RH = 0.0_RP !< Relative Humidity of environment [%] + integer :: NITER_RH = 3 + namelist /PARAM_EXP/ & + TEMP0, & + PTLAPS, & + U0, & + VSHEAR_HVEL, & + ENV_RH, & + NITER_RH + integer :: ierr + + integer :: iq_QV + + integer :: ke, ke2D + integer :: ke_x, ke_y, ke_z + integer :: p, p3, p2D + + real(RP) :: RHOT_hyd(elem%Np,lcmesh%Ne) + real(RP) :: PT_tmp(elem%Np,lcmesh%NeZ,lcmesh%NeX,lcmesh%NeY) + real(RP) :: Rtot (elem%Np,lcmesh%NeZ,lcmesh%NeX,lcmesh%NeY) + real(RP) :: CPtot(elem%Np,lcmesh%NeZ,lcmesh%NeX,lcmesh%NeY) + real(RP) :: CPtot_ov_CVtot(elem%Np,lcmesh%NeZ,lcmesh%NeX,lcmesh%NeY) + real(RP) :: bnd_SFC_PRES(elem%Nnode_h1D**2,lcmesh%lcmesh2D%NeA) + real(RP) :: QV(elem%Np) + real(RP) :: PRES(elem%Np) + real(RP) :: psat0 + + integer :: itr + !----------------------------------------------------------------------------- + + rewind(IO_FID_CONF) + read(IO_FID_CONF,nml=PARAM_EXP,iostat=ierr) + if( ierr < 0 ) then !--- missing + LOG_INFO("BOUNDARY_LAYER_setup",*) 'Not found namelist. Default used.' + elseif( ierr > 0 ) then !--- fatal error + LOG_ERROR("BOUNDARY_LAYER_setup",*) 'Not appropriate names in namelist PARAM_EXP. Check!' + call PRC_abort + endif + LOG_NML(PARAM_EXP) + + !--- + + call hydrostatic_calc_basicstate_constPTLAPS( & + DENS_hyd, PRES_hyd, & ! (out) + PTLAPS, TEMP0, SFC_PRES, x, y, z, lcmesh, elem ) ! (in) + + !$omp parallel do collapse(3) private(ke_z,ke_x,ke_y,ke,ke2D) + do ke_y=1, lcmesh%NeY + do ke_x=1, lcmesh%NeX + do ke_z=1, lcmesh%NeZ + ke2D = ke_x + (ke_y-1)*lcmesh%NeX + ke = ke2D + (ke_z-1)*lcmesh%NeX*lcmesh%NeY + PT_tmp(:,ke_z,ke_x,ke_y) = TEMP0 + PTLAPS * z(:,ke) + RHOT_hyd(:,ke) = DENS_hyd(:,ke) * PT_tmp(:,ke_z,ke_x,ke_y) + end do + end do + end do + + + if ( ENV_RH > 0.0_RP ) then + LOG_INFO("BOUNDARY_LAYER_setup",*) 'Calculate QV from RH' + LOG_INFO("BOUNDARY_LAYER_setup",*) 'ENV_RH = ', ENV_RH, ' [%]' + LOG_INFO("BOUNDARY_LAYER_setup",*) 'NITER_RH = ', NITER_RH + + call TRACER_inq_id( "QV", iq_QV ) + call ATMOS_SATURATION_psat_all( TEMP0, psat0 ) + + do itr=1, NITER_RH + LOG_INFO("BOUNDARY_LAYER_setup",*) 'RH iteration: ', itr + + !$omp parallel do collapse(3) private(ke_z,ke_x,ke_y,ke,ke2D,p3,p2D,p, QV,PRES) + do ke_y=1, lcmesh%NeY + do ke_x=1, lcmesh%NeX + do ke_z=1, lcmesh%NeZ + ke2D = ke_x + (ke_y-1)*lcmesh%NeX + ke = ke2D + (ke_z-1)*lcmesh%NeX*lcmesh%NeY + + do p3=1, elem%Nnode_v + do p2D=1, elem%Nnode_h1D**2 + p = p2D + (p3 - 1)*elem%Nnode_h1D**2 + QV(p) = ENV_RH * 1.0E-2_RP * psat0 / ( ( DENS_hyd(p,ke) + DDENS(p,ke) ) * Rvap * TEMP0 ) + end do + end do + tracer_field_list(iq_QV)%ptr%val(:,ke) = QV(:) + + Rtot (:,ke_z,ke_x,ke_y) = Rdry * ( 1.0_RP - QV(:) ) + Rvap * QV(:) + CPtot (:,ke_z,ke_x,ke_y) = CPdry * ( 1.0_RP - QV(:) ) + CP_VAPOR * QV(:) + CPtot_ov_CVtot(:,ke_z,ke_x,ke_y) = CPtot(:,ke_z,ke_x,ke_y) & + / ( CVdry * ( 1.0_RP - QV(:) ) + CV_VAPOR * QV(:) ) + + PRES(:) = ( DENS_hyd(:,ke) + DDENS(:,ke) ) * Rtot(:,ke_z,ke_x,ke_y) * TEMP0 + PT_tmp(:,ke_z,ke_x,ke_y) = TEMP0 * ( PRES00 / PRES(:) )**( Rtot(:,ke_z,ke_x,ke_y) / CPtot(:,ke_z,ke_x,ke_y) ) + if ( ke_z == 1 ) then + bnd_SFC_PRES(:,ke2D) = PRES(elem%Hslice(:,1)) + end if + end do + end do + end do + + call hydrostatic_build_rho_XYZ( DDENS, & ! (out) + DENS_hyd, PRES_hyd, PT_tmp, & ! (in) + Rtot, CPtot_ov_CVtot, & ! (in) + x, y, z, lcmesh, elem, & ! (in) + bnd_SFC_PRES ) ! (in) + end do ! End of RH iteration + end if + + !$omp parallel do collapse(3) private(ke_z,ke_x,ke_y,ke,ke2D) + do ke_y=1, lcmesh%NeY + do ke_x=1, lcmesh%NeX + do ke_z=1, lcmesh%NeZ + ke2D = ke_x + (ke_y-1)*lcmesh%NeX + ke = ke2D + (ke_z-1)*lcmesh%NeX*lcmesh%NeY + DRHOT(:,ke) = ( DENS_hyd(:,ke) + DDENS(:,ke) ) * PT_tmp(:,ke_z,ke_x,ke_y) & + - RHOT_hyd(:,ke) + + ! MOMX(:,ke) = ( DENS_hyd(:,ke) + DDENS(:,ke) ) * ( U0 + VSHEAR_HVEL * z(:,ke) ) + MOMX(:,ke) = ( DENS_hyd(:,ke) + DDENS(:,ke) ) * ( U0 + 2.0_RP * cos( PI * z(:,ke) / ( dom_zmax - dom_zmin ) ) ) + MOMY(:,ke) = 0.0_RP + MOMZ(:,ke) = 0.0_RP + end do + end do + end do + + return + end subroutine exp_SetInitCond_boundary_layer + +end module mod_user diff --git a/model/atm_nonhydro3d/test/case/boundary_layer/run.conf b/model/atm_nonhydro3d/test/case/boundary_layer/run.conf new file mode 100644 index 00000000..0ee4fcfb --- /dev/null +++ b/model/atm_nonhydro3d/test/case/boundary_layer/run.conf @@ -0,0 +1,135 @@ +#--- Configuration file for a test case of boundary layer ------- +&PARAM_RESTART + IN_BASENAME = "init_00000101-000000.000", + OUTPUT_FLAG = .true., + OUT_BASENAME = 'restart' +/ +&PARAM_TIME + TIME_STARTDATE = 0000, 1, 1, 0, 0, 0, + TIME_STARTMS = 0.D0, + TIME_DURATION = 7200D0, + TIME_DURATION_UNIT = 'SEC', + TIME_DT = 180D0, + TIME_DT_UNIT = 'SEC', +/ +&PARAM_CONST +/ +&PARAM_EXP + TEMP0 = 250.0D0, +/ +&PARAM_USER + USER_do = .true., +/ +#** ATMOS ****************************************************** +&PARAM_ATMOS + ACTIVATE_FLAG = .true., + TIME_DT = 180D0, + TIME_DT_UNIT = 'SEC', + ATMOS_DYN_DO = .true., +! ATMOS_PHY_SF_DO = .true., + ATMOS_PHY_BL_DO = .true., + ATMOS_USE_QV = .true., +/ +&PARAM_ATMOS_MESH + dom_xmin = -5.0D3, + dom_xmax = 5.0D3, + isPeriodicX = .true., + dom_ymin = -5.0D3, + dom_ymax = 5.0D3, + isPeriodicY = .true., + dom_zmin = 0.0D0, + dom_zmax = 1.0D3, + NprcX = 1, + NeX = 1, + NprcY = 1, + NeY = 1, + NeZ = 8, + PolyOrder_h = 7, + PolyOrder_v = 7, + Element_Operation_Type = 'TensorProd3D', +/ +&PARAM_ATMOS_VARS + CHECK_RANGE = .true. , + CHECK_TOTAL = .false., +/ +#** ATMOS / DYN ****************************************************** +&PARAM_ATMOS_DYN + EQS_TYPE = "NONE", + TRACERADV_DISABLE_LIMITER = .true., + !- + TINTEG_TYPE = 'ERK_1s1o', + TINTEG_TYPE_TRACER = 'ERK_1s1o', + TIME_DT = 180D0, + TIME_DT_UNIT = 'SEC', + !- + MODALFILTER_FLAG = .false., + NUMDIFF_FLAG = .false., +/ +&PARAM_ATMOS_DYN_BND + btm_vel_bc = 'SLIP', + top_vel_bc = 'SLIP', + north_vel_bc = 'PERIODIC', + south_vel_bc = 'PERIODIC', + east_vel_bc = 'PERIODIC', + west_vel_bc = 'PERIODIC', + btm_thermal_bc = 'ADIABATIC', + top_thermal_bc = 'ADIABATIC', + north_thermal_bc = 'PERIODIC', + south_thermal_bc = 'PERIODIC', + west_thermal_bc = 'PERIODIC', + east_thermal_bc = 'PERIODIC', +/ + +#** ATMOS / PHYS / SFC ****************************************************** +&PARAM_ATMOS_PHY_SFC + TIME_DT = 180D0, + TIME_DT_UNIT = 'SEC', + SFCFLX_TYPE = 'CONST', + DEFAULT_SFC_TEMP = 300D0, +/ +&PARAM_ATMOS_PHY_SF_CONST + ATMOS_PHY_SF_Const_SH = 0.0D0, + ATMOS_PHY_SF_Const_LH = 0.0D0, + ATMOS_PHY_SF_Const_Cm = 0.0D0, +/ + +#** ATMOS / PHYS / BL ****************************************************** +&PARAM_ATMOS_PHY_BL + TIME_DT = 180D0, + TIME_DT_UNIT = 'SEC', + BL_TYPE = 'MYNN_LEVEL2', + C_IP = 0.01D0, +/ + +#*** OUTPUT ******************************************* +&PARAM_FILE_HISTORY + FILE_HISTORY_DEFAULT_BASENAME = "history", + FILE_HISTORY_DEFAULT_TINTERVAL = 180.0D0, + FILE_HISTORY_DEFAULT_TUNIT = "SEC", + FILE_HISTORY_DEFAULT_TAVERAGE = .false., + FILE_HISTORY_DEFAULT_DATATYPE = "REAL4", + FILE_HISTORY_OUTPUT_STEP0 = .true., +/ +!&HISTORY_ITEM name='DDENS' / +&HISTORY_ITEM name='T' / +&HISTORY_ITEM name='U' / +&HISTORY_ITEM name='BL_RHOU_t' / +&HISTORY_ITEM name='BL_RHOT_t' / +&HISTORY_ITEM name='NU' / +&HISTORY_ITEM name='KH' / +!&HISTORY_ITEM name='QV' / + +#*** Statistics ******************************************* + +&PARAM_MESHFIELD_STATISTICS + use_globalcomm = .true., +/ +&PARAM_MONITOR + MONITOR_STEP_INTERVAL = 10 +/ +&MONITOR_ITEM name='DDENS' / +&MONITOR_ITEM name='ENGT' / +&MONITOR_ITEM name='ENGK' / +&MONITOR_ITEM name='ENGI' / +&MONITOR_ITEM name='ENGP' / + diff --git a/model/atm_nonhydro3d/test/case/boundary_layer/visualize/visualize.sh b/model/atm_nonhydro3d/test/case/boundary_layer/visualize/visualize.sh new file mode 100644 index 00000000..920a3099 --- /dev/null +++ b/model/atm_nonhydro3d/test/case/boundary_layer/visualize/visualize.sh @@ -0,0 +1,27 @@ +#!/bin/bash -x + +echo "+make directory" +mkdir -p analysis + +### check error norm ### +# echo "+mkgraph monitor" +# for var in DDENS ENGT ENGP ENGK ENGI; do +# python ../common/cmd_analysis_monitor.py monitor.peall ${var} 0.25 analysis/monitor_${var}.png +# done + +### make figures ### +echo "+mkgraph" +python ../common/cmd_mkgraph.py history.pe00\*.nc@T,x=0e0,y=0e0 analysis/T.png --prc_num_xy 1 1 --figsize 8 4 +python ../common/cmd_mkgraph.py history.pe00\*.nc@BL_RHOU_t,x=0e0,y=0e0 analysis/BL_RHOU_t.png --prc_num_xy 1 1 --figsize 8 4 --range -4e-4 4e-4 +python ../common/cmd_mkgraph.py history.pe00\*.nc@BL_RHOT_t,x=0e0,y=0e0 analysis/BL_RHOT_t.png --prc_num_xy 1 1 --figsize 8 4 +python ../common/cmd_mkgraph.py history.pe00\*.nc@NU,x=0e0,y=0e0 analysis/NU.png --prc_num_xy 1 1 --figsize 8 4 --range 0 30 +python ../common/cmd_mkgraph.py history.pe00\*.nc@KH,x=0e0,y=0e0 analysis/KH.png --prc_num_xy 1 1 --figsize 8 4 --range 0 30 +python ../common/cmd_mkgraph.py history.pe00\*.nc@U,x=0e0,y=0e0 analysis/U.png --prc_num_xy 1 1 --figsize 8 4 + +### make animation ### +# echo "+make animation" +# python visualize/mkanim.py analysis/advdiff1d.mp4 -1.1 1.1 + +### check error norm ### +# echo "+check numerical errors" +# python visualize/mkgraph_numerror.py LOG_NUMERROR.peall analysis/numerror_