From cc0b904e9d6b20b1398bcf0e588fd9e4a1e92672 Mon Sep 17 00:00:00 2001 From: Yuta Kawai Date: Sun, 26 Jul 2026 14:46:53 +0900 Subject: [PATCH 1/9] Set var0_*D values when euler scheme is used --- FElib/src/common/scale_timeint_rk.F90 | 42 ++++------------------- FElib/src/common/scale_timeint_rk.F90.erb | 24 ++----------- 2 files changed, 8 insertions(+), 58 deletions(-) diff --git a/FElib/src/common/scale_timeint_rk.F90 b/FElib/src/common/scale_timeint_rk.F90 index 150472a8..f083c011 100644 --- a/FElib/src/common/scale_timeint_rk.F90 +++ b/FElib/src/common/scale_timeint_rk.F90 @@ -1895,10 +1895,11 @@ subroutine rk_advance_general1D( this, nowstage, q, varID, is, ie , IA, var_num, tintbuf_ind = this%tend_buf_indmap(nowstage) - if ( this%nstage == 1 ) then + if ( this%nstage == 1 .and. (.not. this%imex_flag) ) then !$omp parallel do !$acc parallel loop collapse(1) present( q, varTmp_1d ) do i=is, ie + var0_1d(i,varID) = q(i) varTmp_1d(i,varID) = q(i) end do end if @@ -1933,15 +1934,6 @@ subroutine rk_advance_general1D( this, nowstage, q, varID, is, ie , IA, var_num, !$omp private( i ) !$acc parallel present( q, var0_1D, varTmp_1D, tend_buf1D_ex, tend_buf1D_im ) - if ( nowstage == 1 .and. (.not. this%imex_flag) ) then - !$omp do - !$acc loop - do i=is, ie - var0_1D(i,varID) = q(i) - varTmp_1D(i,varID) = q(i) - end do - end if - if ( this%imex_flag ) then coef_b_ex_dt = this%coef_b_ex(nowstage) * this%dt coef_b_im_dt = this%coef_b_im(nowstage) * this%dt @@ -2227,11 +2219,12 @@ subroutine rk_advance_general2D( this, nowstage, q, varID, is, ie ,js, je , IA,J tintbuf_ind = this%tend_buf_indmap(nowstage) - if ( this%nstage == 1 ) then + if ( this%nstage == 1 .and. (.not. this%imex_flag) ) then !$omp parallel do !$acc parallel loop collapse(2) present( q, varTmp_2d ) do j=js, je do i=is, ie + var0_2d(i,j,varID) = q(i,j) varTmp_2d(i,j,varID) = q(i,j) end do end do @@ -2271,17 +2264,6 @@ subroutine rk_advance_general2D( this, nowstage, q, varID, is, ie ,js, je , IA,J !$omp private( i,j ) !$acc parallel present( q, var0_2D, varTmp_2D, tend_buf2D_ex, tend_buf2D_im ) - if ( nowstage == 1 .and. (.not. this%imex_flag) ) then - !$omp do - !$acc loop collapse(2) - do j=js, je - do i=is, ie - var0_2D(i,j,varID) = q(i,j) - varTmp_2D(i,j,varID) = q(i,j) - end do - end do - end if - if ( this%imex_flag ) then coef_b_ex_dt = this%coef_b_ex(nowstage) * this%dt coef_b_im_dt = this%coef_b_im(nowstage) * this%dt @@ -2601,12 +2583,13 @@ subroutine rk_advance_general3D( this, nowstage, q, varID, is, ie ,js, je ,ks, k tintbuf_ind = this%tend_buf_indmap(nowstage) - if ( this%nstage == 1 ) then + if ( this%nstage == 1 .and. (.not. this%imex_flag) ) then !$omp parallel do collapse(2) !$acc parallel loop collapse(3) present( q, varTmp_3d ) do k=ks, ke do j=js, je do i=is, ie + var0_3d(i,j,k,varID) = q(i,j,k) varTmp_3d(i,j,k,varID) = q(i,j,k) end do end do @@ -2651,19 +2634,6 @@ subroutine rk_advance_general3D( this, nowstage, q, varID, is, ie ,js, je ,ks, k !$omp private( i,j,k ) !$acc parallel present( q, var0_3D, varTmp_3D, tend_buf3D_ex, tend_buf3D_im ) - if ( nowstage == 1 .and. (.not. this%imex_flag) ) then - !$omp do collapse(2) - !$acc loop collapse(3) - do k=ks, ke - do j=js, je - do i=is, ie - var0_3D(i,j,k,varID) = q(i,j,k) - varTmp_3D(i,j,k,varID) = q(i,j,k) - end do - end do - end do - end if - if ( this%imex_flag ) then coef_b_ex_dt = this%coef_b_ex(nowstage) * this%dt coef_b_im_dt = this%coef_b_im(nowstage) * this%dt diff --git a/FElib/src/common/scale_timeint_rk.F90.erb b/FElib/src/common/scale_timeint_rk.F90.erb index b54e4afa..af769a26 100644 --- a/FElib/src/common/scale_timeint_rk.F90.erb +++ b/FElib/src/common/scale_timeint_rk.F90.erb @@ -998,7 +998,7 @@ contains tintbuf_ind = this%tend_buf_indmap(nowstage) - if ( this%nstage == 1 ) then + if ( this%nstage == 1 .and. (.not. this%imex_flag) ) then % if (d > 2) then !$omp parallel do collapse(<%=(d-1)%>) % else @@ -1008,6 +1008,7 @@ contains % for i in 1..d do <%=ind_name[d-i]%>=<%=ind_range_list[d-i]%> % end + var0_<%=d%>d(<%=ind_ary%>,varID) = q(<%=ind_ary%>) varTmp_<%=d%>d(<%=ind_ary%>,varID) = q(<%=ind_ary%>) % for i in 1..d end do @@ -1060,27 +1061,6 @@ contains !$omp private( <%=ind_ary%> ) !$acc parallel present( q, var0_<%=d%>D, varTmp_<%=d%>D, tend_buf<%=d%>D_ex, tend_buf<%=d%>D_im ) - if ( nowstage == 1 .and. (.not. this%imex_flag) ) then -% if (d > 2) then - !$omp do collapse(<%=(d-1)%>) -% else - !$omp do -% end -% if (d >= 2) then - !$acc loop collapse(<%=(d)%>) -% else - !$acc loop -% end -% for i in 1..d - do <%=ind_name[d-i]%>=<%=ind_range_list[d-i]%> -% end - var0_<%=d%>D(<%=ind_ary%>,varID) = q(<%=ind_ary%>) - varTmp_<%=d%>D(<%=ind_ary%>,varID) = q(<%=ind_ary%>) -% for i in 1..d - end do -% end - end if - if ( this%imex_flag ) then coef_b_ex_dt = this%coef_b_ex(nowstage) * this%dt coef_b_im_dt = this%coef_b_im(nowstage) * this%dt From 044e176123c2de9833670300ddff20e5c54cc485 Mon Sep 17 00:00:00 2001 From: Yuta Kawai Date: Sun, 26 Jul 2026 14:47:48 +0900 Subject: [PATCH 2/9] Add local mesh assignment in tracer advection loop --- FElib/src/fluid_dyn_solver/scale_atm_dyn_dgm_driver_trcadv3d.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/FElib/src/fluid_dyn_solver/scale_atm_dyn_dgm_driver_trcadv3d.F90 b/FElib/src/fluid_dyn_solver/scale_atm_dyn_dgm_driver_trcadv3d.F90 index e66a561d..1e63521a 100644 --- a/FElib/src/fluid_dyn_solver/scale_atm_dyn_dgm_driver_trcadv3d.F90 +++ b/FElib/src/fluid_dyn_solver/scale_atm_dyn_dgm_driver_trcadv3d.F90 @@ -528,6 +528,7 @@ subroutine AtmDynDGMDriver_trcadv3d_update( this, & end do ! end for RK loop do n=1, mesh3D%LOCAL_MESH_NUM + lcmesh3D => mesh3D%lcmesh_list(n) !$omp parallel do do ke=lcmesh3D%NeS, lcmesh3D%NeE QTRC%local(n)%val(:,ke) = ( DENS_hyd%local(n)%val(:,ke) + DDENS_TRC%local(n)%val(:,ke) ) & From f43d9959575a19fc7e6b9dd9101b85bb8166cc6b Mon Sep 17 00:00:00 2001 From: Yuta Kawai Date: Sun, 26 Jul 2026 14:48:51 +0900 Subject: [PATCH 3/9] Implement a moist convective adjustment scheme --- .../scale_atm_phy_cp_dgm_mconv_adjustment.F90 | 911 ++++++++++++++++++ 1 file changed, 911 insertions(+) create mode 100644 FElib/src/cumulus/scale_atm_phy_cp_dgm_mconv_adjustment.F90 diff --git a/FElib/src/cumulus/scale_atm_phy_cp_dgm_mconv_adjustment.F90 b/FElib/src/cumulus/scale_atm_phy_cp_dgm_mconv_adjustment.F90 new file mode 100644 index 00000000..1901b3b2 --- /dev/null +++ b/FElib/src/cumulus/scale_atm_phy_cp_dgm_mconv_adjustment.F90 @@ -0,0 +1,911 @@ +!> module FElib / Atmosphere / Physics cumulus parameterization +!! +!! @par Description +!! A module to provide a moist convective adjustment scheme for cumulus parameterization in atmospheric model +!! +!! @author Yuta Kawai, Team SCALE +!! +!! @par Reference +!! +!------------------------------------------------------------------------------- +#include "scaleFElib.h" +module scale_atm_phy_cp_dgm_mconv_adjustment + !----------------------------------------------------------------------------- + ! + !++ Used modules + ! + use scale_precision + use scale_io + use scale_prc + use scale_prof + use scale_const, only: & + GRAV => CONST_GRAV, & + LHV0 => CONST_LHV0, & + CPdry => CONST_CPDRY, & + CVdry => CONST_CVdry, & + Rdry => CONST_Rdry, & + PRES0 => CONST_PRE00, & + EPS => CONST_EPS + use scale_atmos_hydrometeor, only: & + CV_VAPOR, & + CP_VAPOR + use scale_atmos_saturation, only: & + ATMOS_SATURATION_psat_liq, & + ATMOS_SATURATION_pres2qsat_liq + + 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 + + !----------------------------------------------------------------------------- + implicit none + private + !----------------------------------------------------------------------------- + ! + !++ Public type & procedure + ! + public :: atm_phy_cp_dgm_mconv_adjustment_setup + public :: atm_phy_cp_dgm_mconv_adjustment_calc_tendency + public :: atm_phy_cp_dgm_mconv_adjustment_finalize + + !----------------------------------------------------------------------------- + !++ Public parameters & variables + ! + !----------------------------------------------------------------------------- + ! + !++ Private procedure + ! + !----------------------------------------------------------------------------- + ! + integer, parameter :: MCA_MAX_ADJUST_ITER = 10 + integer, parameter :: MCA_MAX_ENERGY_ITER = 40 + integer, parameter :: MCA_MAX_LOCAL_ITER = 20 + + real(RP), parameter :: MCA_RH_TRIGGER = 1.00_RP + real(RP), parameter :: MCA_ENERGY_RTOL = 1.0E-10_RP + real(RP), parameter :: MCA_TEMP_TOL = 1.0E-8_RP + + ! Parameters to avoid detection of numerical noise + real(RP), parameter :: MCA_HMSE_GRAD_TOL = 0.0_RP + real(RP), parameter :: MCA_MIN_UNSTABLE_DEPTH = 0.0_RP + real(RP), parameter :: MCA_Z_TOL = 100.0_RP * epsilon(1.0_RP) + +contains + subroutine atm_phy_cp_dgm_mconv_adjustment_setup() + implicit none + !---------------------------------------------------- + return + end subroutine atm_phy_cp_dgm_mconv_adjustment_setup + +!OCL SERIAL + subroutine atm_phy_cp_dgm_mconv_adjustment_calc_tendency( & + DENS_t, RHOT_t, RHOQV_t, SFLX_RAIN, & + DDENS, DRHOT, QV, PT, PRES, DENS_hyd, Rtot, CPtot, & + dtsec, lmesh, elem, elem1D ) + implicit none + class(LocalMesh3D), intent(in) :: lmesh + class(ElementBase3D), intent(in) :: elem + class(ElementBase1D), intent(in) :: elem1D + real(RP), intent(out) :: DENS_t(elem%Np,lmesh%NeA) + real(RP), intent(out) :: RHOT_t(elem%Np,lmesh%NeA) + real(RP), intent(out) :: RHOQV_t(elem%Np,lmesh%NeA) + real(RP), intent(out) :: SFLX_RAIN(elem%Nnode_h1D**2,lmesh%Ne2DA) + real(RP), intent(in) :: DDENS(elem%Np,lmesh%NeA) + real(RP), intent(in) :: DRHOT(elem%Np,lmesh%NeA) + real(RP), intent(in) :: QV(elem%Np,lmesh%NeA) + real(RP), intent(in) :: PT(elem%Np,lmesh%NeA) + real(RP), intent(in) :: PRES(elem%Np,lmesh%NeA) + real(RP), intent(in) :: DENS_hyd(elem%Np,lmesh%NeA) + real(RP), intent(in) :: Rtot(elem%Np,lmesh%NeA) + real(RP), intent(in) :: CPtot(elem%Np,lmesh%NeA) + real(RP), intent(in) :: dtsec + + integer :: ke, ke_xy, ke_z + integer :: ph, pz, p + real(RP) :: dens_z(elem%Nnode_v,lmesh%NeZ) + real(RP) :: pres_z(elem%Nnode_v,lmesh%NeZ) + real(RP) :: zlev_z(elem%Nnode_v,lmesh%NeZ) + real(RP) :: temp_z(elem%Nnode_v,lmesh%NeZ) + real(RP) :: pott_z(elem%Nnode_v,lmesh%NeZ) ! work array + real(RP) :: qvap_z(elem%Nnode_v,lmesh%NeZ) ! work array + real(RP) :: rtot_(elem%Nnode_v) + real(RP) :: cptot_(elem%Nnode_v) + + real(RP) :: elem_width_z + real(RP) :: int_weight(elem%Nnode_v,lmesh%NeZ) + + real(RP) :: pott_ini(elem%Nnode_v,lmesh%NeZ) + real(RP) :: qvap_ini(elem%Nnode_v,lmesh%NeZ) + real(RP) :: dens_ini(elem%Nnode_v,lmesh%NeZ) + real(RP) :: rhoqvap_ini(elem%Nnode_v,lmesh%NeZ) + + integer :: lbase, ltop + logical :: is_unstable + logical :: adjust_mask(elem%Nnode_v,lmesh%NeZ) + + integer :: iter_adj + real(RP) :: temp_aj(elem%Nnode_v,lmesh%NeZ) + real(RP) :: qdry_aj(elem%Nnode_v,lmesh%NeZ) + real(RP) :: qvap_aj(elem%Nnode_v,lmesh%NeZ) + real(RP) :: rhoqvap_aj(elem%Nnode_v,lmesh%NeZ) + real(RP) :: rhoprecip_aj(elem%Nnode_v,lmesh%NeZ) ! Condensate water density [kg/m^3] with each iteration, which is removed by precipitation + real(RP) :: rhoprecip_accum(elem%Nnode_v,lmesh%NeZ) + real(RP) :: rhoqdry_aj(elem%Nnode_v,lmesh%NeZ) + + integer :: l + real(RP), allocatable :: pott_diag(:) + real(RP), allocatable :: qvap_diag(:) + real(RP), allocatable :: pres_diag(:) + real(RP), allocatable :: temp_diag(:) + real(RP), allocatable :: zlev_diag(:) + real(RP) :: qdry_diag, Rtot_diag, CPtot_diag + integer :: nlev_diag + + logical :: is_profile_converged + logical :: do_adjustment + logical :: adjustment_converged + + integer :: Hslice_b(elem%Nnode_h1D**2), Hslice_t(elem%Nnode_h1D**2) + + real(RP) :: Rvap + + integer, parameter :: MAX_ITER_ADJ = 10 + !---------------------------------------------------- + + nlev_diag = lmesh%NeZ * (elem%Nnode_v-1) + 1 + allocate( pott_diag(nlev_diag), qvap_diag(nlev_diag), pres_diag(nlev_diag), & + temp_diag(nlev_diag), zlev_diag(nlev_diag) ) + + Hslice_b(:) = elem%Hslice(:,1) + Hslice_t(:) = elem%Hslice(:,elem%Nnode_v) + + Rvap = CP_VAPOR - CV_VAPOR + + !$omp parallel do collapse(2) private(ke,p, & + !$omp dens_z, pres_z, temp_z, zlev_z, pott_z, qvap_z, rtot_,cptot_, & + !$omp dens_ini, pott_ini, qvap_ini, rhoqvap_ini, & + !$omp pott_diag, qvap_diag, pres_diag, temp_diag, zlev_diag, qdry_diag, Rtot_diag, CPtot_diag, & + !$omp lbase, ltop, is_unstable, adjust_mask, & + !$omp temp_aj, qdry_aj, qvap_aj, rhoqdry_aj, rhoqvap_aj, rhoprecip_aj, rhoprecip_accum, & + !$omp is_profile_converged, do_adjustment, adjustment_converged, & + !$omp int_weight,elem_width_z ) + do ke_xy=1, lmesh%Ne2D + do ph=1, elem%Nnode_h1D**2 + + !- Extract vertical 1D DG column + + do ke_z=1, lmesh%NeZ + ke = ke_xy + (ke_z-1)*lmesh%Ne2D + + elem_width_z = lmesh%zlev(Hslice_t(ph),ke) - lmesh%zlev(Hslice_b(ph),ke) + do pz=1, elem%Nnode_v + p = ph + (pz-1)*elem%Nnode_h1D**2 + dens_z(pz,ke_z) = DDENS(p,ke) + DENS_hyd(p,ke) + pres_z(pz,ke_z) = PRES(p,ke) + zlev_z(pz,ke_z) = lmesh%zlev(p,ke) + int_weight(pz,ke_z) = elem1D%IntWeight_lgl(pz) * 0.5_RP * elem_width_z + + temp_z(pz,ke_z) = PRES(p,ke)/ ( dens_z(pz,ke_z) * Rtot(p,ke) ) + qvap_z(pz,ke_z) = QV(p,ke) + pott_z(pz,ke_z) = PT(p,ke) + end do + end do + + !- Set initial state + + do ke_z=1, lmesh%NeZ + dens_ini(:,ke_z) = dens_z(:,ke_z) + pott_ini(:,ke_z) = pott_z(:,ke_z) + qvap_ini(:,ke_z) = qvap_z(:,ke_z) + rhoqvap_ini(:,ke_z) = dens_ini(:,ke_z) * qvap_ini(:,ke_z) + + rhoprecip_accum(:,ke_z) = 0.0_RP + end do + + !** Loop for iteration of the moist convective adjustment + + adjustment_converged = .false. + do_adjustment = .false. + do iter_adj=1, MAX_ITER_ADJ + + !- Build a vertical profile which is continuous at vertical element interfaces from DG solution + + call build_diag_zprofile_from_dg( & + pott_diag, qvap_diag, pres_diag, zlev_diag, & ! (out) + pott_z, qvap_z, pres_z, zlev_z, & ! (in) + lmesh, elem, nlev_diag ) ! (in) + + do l=1, nlev_diag + qdry_diag = 1.0_RP - qvap_diag(l) + CPtot_diag = CPdry * qdry_diag + CP_VAPOR * qvap_diag(l) + Rtot_diag = Rdry * qdry_diag + Rvap * qvap_diag(l) + temp_diag(l) = pott_diag(l) * ( pres_diag(l) / PRES0 )**( Rtot_diag / CPtot_diag ) + end do + + !- Diagnose unstable layers + + call diagnose_convective_layer( lbase, ltop, is_unstable, & ! (out) + temp_diag, pott_diag, qvap_diag, pres_diag, zlev_diag, & ! (in) + nlev_diag ) ! (in) + + ! write(*,*) "-- ph=", ph, "ke_xy=", ke_xy, "iter_adj=", iter_adj + ! write(*,*) " lbase, ltop, is_unstable = ", lbase, ltop, is_unstable + + if ( .not. is_unstable ) then + adjustment_converged = .true. + exit + end if + + do_adjustment = .true. + + !- Generate a mask for DG vertical nodes where adjustment is needed + + call make_dg_adjustment_mask( adjust_mask, & ! (out) + zlev_z, zlev_diag(lbase), zlev_diag(ltop), elem%Nnode_v, lmesh%NeZ ) ! (in) + + !- Construct a moist neutral profile in the unstable layer + + call construct_moist_neutral_profile( & + temp_aj, qvap_aj,rhoqvap_aj, rhoprecip_aj, is_profile_converged, & ! (out) + qvap_z, pres_z, zlev_z, dens_z, temp_z, & ! (in) + int_weight, adjust_mask, elem%Nnode_v, lmesh%NeZ ) ! (in) + + if ( .not. is_profile_converged ) then + LOG_INFO("atm_phy_cp_dgm_mconv_adjustment_calc_tendency",*) "Moist neutral profile construction did not converge. Aborting adjustment." + LOG_INFO("atm_phy_cp_dgm_mconv_adjustment_calc_tendency",*) "iter_adj=", iter_adj, "ph=", ph, "ke_xy=", ke_xy + call PRC_abort + end if + + !- Immediate precipitation removal + + rhoprecip_accum(:,:) = rhoprecip_accum(:,:) + rhoprecip_aj(:,:) + + rhoqdry_aj(:,:) = dens_z(:,:) - rhoqvap_aj(:,:) - rhoprecip_aj(:,:) + dens_z(:,:) = dens_z(:,:) - rhoprecip_aj(:,:) + + qvap_aj(:,:) = rhoqvap_aj(:,:) / dens_z(:,:) + qdry_aj(:,:) = rhoqdry_aj(:,:) / dens_z(:,:) + + !- Use the post-precipitation state in the next iteration + + do ke_z=1, lmesh%NeZ + pres_z(:,ke_z) = ( Rdry * rhoqdry_aj(:,ke_z) + Rvap * rhoqvap_aj(:,ke_z) ) * temp_aj(:,ke_z) + temp_z(:,ke_z) = temp_aj(:,ke_z) + qvap_z(:,ke_z) = qvap_aj(:,ke_z) + + rtot_ (:) = Rdry * qdry_aj(:,ke_z) + Rvap * qvap_aj(:,ke_z) + cptot_(:) = CPdry * qdry_aj(:,ke_z) + CP_VAPOR * qvap_aj(:,ke_z) + pott_z(:,ke_z) = temp_aj(:,ke_z) * ( PRES0 / pres_z(:,ke_z) )**( rtot_(:) / cptot_(:) ) + end do + end do ! end loop for iteration + + if ( .not. adjustment_converged ) then + LOG_INFO("atm_phy_cp_dgm_mconv_adjustment_calc_tendency",*) "Moist convective adjustment did not converge." + LOG_INFO("atm_phy_cp_dgm_mconv_adjustment_calc_tendency",*) "ph=", ph, "ke_xy=", ke_xy + call PRC_abort + end if + + + !- Convert the adjusted state to DG tendencies + + if ( do_adjustment ) then + do ke_z=1, lmesh%NeZ + ke = ke_xy + (ke_z-1)*lmesh%Ne2D + do pz=1, elem%Nnode_v + p = ph + (pz-1)*elem%Nnode_h1D**2 + DENS_t (p,ke) = ( dens_z(pz,ke_z) - dens_ini(pz,ke_z) ) / dtsec + RHOT_t (p,ke) = ( dens_z(pz,ke_z) * pott_z(pz,ke_z) - dens_ini(pz,ke_z) * pott_ini(pz,ke_z) ) / dtsec + RHOQV_t(p,ke) = ( dens_z(pz,ke_z) * qvap_z(pz,ke_z) - rhoqvap_ini(pz,ke_z) ) / dtsec + end do + end do + else + do ke_z=1, lmesh%NeZ + ke = ke_xy + (ke_z-1)*lmesh%Ne2D + do pz=1, elem%Nnode_v + p = ph + (pz-1)*elem%Nnode_h1D**2 + DENS_t (p,ke) = 0.0_RP + RHOT_t (p,ke) = 0.0_RP + RHOQV_t(p,ke) = 0.0_RP + end do + end do + end if + + SFLX_RAIN(ph,ke_xy) = sum( rhoprecip_accum(:,:) * int_weight(:,:) ) / dtsec + + end do ! end loop for ph + end do ! end loop for ke_xy + + return + end subroutine atm_phy_cp_dgm_mconv_adjustment_calc_tendency + + subroutine atm_phy_cp_dgm_mconv_adjustment_finalize() + implicit none + !---------------------------------------------------- + return + end subroutine atm_phy_cp_dgm_mconv_adjustment_finalize + + +!- private subroutines +!OCL SERIAL + subroutine build_diag_zprofile_from_dg( & + pott_diag, qvap_diag, pres_diag, zlev_diag, & + pott_z, qvap_z, pres_z, zlev_z, & + lmesh, elem, nlev ) + implicit none + class(LocalMesh3D), intent(in) :: lmesh + class(ElementBase3D), intent(in) :: elem + integer, intent(in) :: nlev + real(RP), intent(out) :: pott_diag(nlev) + real(RP), intent(out) :: qvap_diag(nlev) + real(RP), intent(out) :: pres_diag(nlev) + real(RP), intent(out) :: zlev_diag(nlev) + real(RP), intent(in) :: pott_z(elem%Nnode_v,lmesh%NeZ) + real(RP), intent(in) :: qvap_z(elem%Nnode_v,lmesh%NeZ) + real(RP), intent(in) :: pres_z(elem%Nnode_v,lmesh%NeZ) + real(RP), intent(in) :: zlev_z(elem%Nnode_v,lmesh%NeZ) + + integer :: pz, ke_z + integer :: l + integer :: Nnode_v + !---------------------------------------------------- + + Nnode_v = elem%Nnode_v + l = 1 + + pott_diag(l) = pott_z(1,1) + qvap_diag(l) = qvap_z(1,1) + pres_diag(l) = pres_z(1,1) + zlev_diag(l) = zlev_z(1,1) + do ke_z=1, lmesh%NeZ + do pz=2, Nnode_v-1 + l = l + 1 + pott_diag(l) = pott_z(pz,ke_z) + qvap_diag(l) = qvap_z(pz,ke_z) + pres_diag(l) = pres_z(pz,ke_z) + zlev_diag(l) = zlev_z(pz,ke_z) + end do + l = l + 1 + if ( ke_z < lmesh%NeZ ) then + pott_diag(l) = 0.5_RP * ( pott_z(Nnode_v,ke_z) + pott_z(1,ke_z+1) ) + qvap_diag(l) = 0.5_RP * ( qvap_z(Nnode_v,ke_z) + qvap_z(1,ke_z+1) ) + pres_diag(l) = 0.5_RP * ( pres_z(Nnode_v,ke_z) + pres_z(1,ke_z+1) ) + zlev_diag(l) = 0.5_RP * ( zlev_z(Nnode_v,ke_z) + zlev_z(1,ke_z+1) ) + else + pott_diag(l) = pott_z(Nnode_v,ke_z) + qvap_diag(l) = qvap_z(Nnode_v,ke_z) + pres_diag(l) = pres_z(Nnode_v,ke_z) + zlev_diag(l) = zlev_z(Nnode_v,ke_z) + end if + end do + return + end subroutine build_diag_zprofile_from_dg + +!OCL SERIAL + subroutine diagnose_convective_layer( lbase, ltop, is_unstable, & + temp, pott, qvap, pres, zlev, nlev ) + implicit none + integer, intent(out) :: lbase + integer, intent(out) :: ltop + logical, intent(out) :: is_unstable + integer, intent(in) :: nlev + real(RP), intent(in) :: temp(nlev) + real(RP), intent(in) :: pott(nlev) + real(RP), intent(in) :: qvap(nlev) + real(RP), intent(in) :: pres(nlev) + real(RP), intent(in) :: zlev(nlev) + + integer :: l + real(RP) :: qsat(nlev) + real(RP) :: rh(nlev) + real(RP) :: hmse_sat(nlev) + real(RP) :: cptot + + real(RP) :: dz + real(RP) :: dhmse_sat_dz + real(RP) :: rh_lyr + logical :: unstable_pair + logical :: is_inside + !------------------------------------ + + call ATMOS_SATURATION_pres2qsat_liq( & + nlev, 1, nlev, & ! (in) + temp, pres, & ! (in) + qsat ) ! (out) + + do l=1, nlev + rh(l) = qvap(l) / max(qsat(l), EPS) + + cptot = CPdry * (1.0_RP - qsat(l)) + CP_VAPOR * qsat(l) + hmse_sat(l) = cptot * temp(l) + GRAV * zlev(l) + LHV0 * qsat(l) + end do + + is_unstable = .false. + is_inside = .false. + lbase = 0; ltop = 0 + + do l=1, nlev-1 + dz = zlev(l+1) - zlev(l) + dhmse_sat_dz = ( hmse_sat(l+1) - hmse_sat(l) ) / dz + rh_lyr = 0.5_RP * ( rh(l) + rh(l+1) ) + ! write(*,*) "l=", l, "zlev(l)=", zlev(l), "zlev(l+1)=", zlev(l+1), & + ! "rh_lyr=", rh_lyr, "dhmse_sat_dz=", dhmse_sat_dz + unstable_pair = ( rh_lyr >= MCA_RH_TRIGGER ) & + .and. ( dhmse_sat_dz < - MCA_HMSE_GRAD_TOL ) + + if ( unstable_pair ) then + if ( .not. is_inside ) then + lbase = l + is_inside = .true. + end if + ltop = l + 1 + else if ( is_inside ) then + exit + end if + end do + + if ( is_inside ) then + if ( zlev(ltop) - zlev(lbase) >= MCA_MIN_UNSTABLE_DEPTH ) is_unstable = .true. + end if + return + end subroutine diagnose_convective_layer + +!OCL SERIAL + subroutine make_dg_adjustment_mask( adjust_mask, & ! (out) + zlev, zbase, ztop, npz, nez ) ! (in) + implicit none + integer, intent(in) :: npz + integer, intent(in) :: nez + logical, intent(out) :: adjust_mask(npz,nez) + real(RP), intent(in) :: zlev(npz,nez) + real(RP), intent(in) :: zbase + real(RP), intent(in) :: ztop + + integer :: ke_z, pz + !---------------------------------------------------- + + adjust_mask(:,:) = .false. + do ke_z = 1, nez + do pz = 1, npz + if ( zlev(pz,ke_z) >= zbase - MCA_Z_TOL .and. & + zlev(pz,ke_z) <= ztop + MCA_Z_TOL ) then + adjust_mask(pz,ke_z) = .true. + end if + end do + end do + return + end subroutine make_dg_adjustment_mask + +!OCL SERIAL + subroutine construct_moist_neutral_profile( & + temp_aj, qvap_aj, rhoqvap_aj, rhoprecip_aj, & + converged, & + qvap, pres, zlev, dens, temp, & + int_weight, adjust_mask, npz, nez ) + + implicit none + integer, intent(in) :: npz + integer, intent(in) :: nez + real(RP), intent(out) :: temp_aj(npz,nez) + real(RP), intent(out) :: qvap_aj(npz,nez) + real(RP), intent(out) :: rhoqvap_aj(npz,nez) + real(RP), intent(out) :: rhoprecip_aj(npz,nez) + logical, intent(out) :: converged + real(RP), intent(in) :: qvap(npz,nez) + real(RP), intent(in) :: pres(npz,nez) + real(RP), intent(in) :: zlev(npz,nez) + real(RP), intent(in) :: dens(npz,nez) + real(RP), intent(in) :: temp(npz,nez) + real(RP), intent(in) :: int_weight(npz,nez) + logical, intent(in) :: adjust_mask(npz,nez) + + real(RP) :: energy_target + real(RP) :: energy_aj + real(RP) :: residual + + real(RP) :: tbase_lo + real(RP) :: tbase_hi + real(RP) :: tbase_mid + + integer :: ke_z, pz + integer :: iter + + real(RP) :: qsat_aj(npz,nez) + + real(RP) :: water_mass_ini + real(RP) :: water_mass_aj + real(RP) :: precip_mass_total + real(RP) :: precip_mass_local + real(RP) :: positive_cond_mass + real(RP) :: precip_scale + + logical :: found_base + real(RP), parameter :: MCA_TBASE_RANGE = 40.0_RP + + logical :: is_converged_ref_profile + logical :: water_feasible + logical :: residual_evaluated + !---------------------------------------------------- + + converged = .false. + + ! Initialize the adjusted profile with the original state + + temp_aj(:,:) = temp(:,:) + qvap_aj(:,:) = qvap(:,:) + rhoqvap_aj(:,:) = dens(:,:) * qvap(:,:) + rhoprecip_aj(:,:) = 0.0_RP + + !- Calculate initial column-integrated MSE + + call column_mse( energy_target, & ! (out) + temp, qvap, zlev, dens, & ! (in) + int_weight, adjust_mask, npz, nez ) ! (in) + + !- Find temperature at the lowest adjusted node + + tbase_mid = 0.0_RP + found_base = .false. + do ke_z = 1, nez + do pz = 1, npz + if ( adjust_mask(pz,ke_z) ) then + tbase_mid = temp(pz,ke_z) + found_base = .true. + exit + end if + end do + if ( found_base ) exit + end do + if ( .not. found_base ) then + write(*,*) "Error: No adjusted nodes found in the column. Aborting adjustment." + return + end if + + !- Calculate the total water mass in the column before adjustment + + water_mass_ini = 0.0_RP + do ke_z = 1, nez + do pz = 1, npz + if ( adjust_mask(pz,ke_z) ) then + water_mass_ini = water_mass_ini + int_weight(pz,ke_z) * dens(pz,ke_z) * qvap(pz,ke_z) + end if + end do + end do + + !- Set the initial bounds for the base temperature search + + tbase_lo = tbase_mid - MCA_TBASE_RANGE + tbase_hi = tbase_mid + MCA_TBASE_RANGE + + water_feasible = .false. + residual_evaluated = .false. + + !- + do iter=1, MCA_MAX_ENERGY_ITER + + tbase_mid = 0.5_RP * ( tbase_lo + tbase_hi ) + + ! Evaluate the moist adiabat temperature profile for the current base temperature + + temp_aj(:,:) = temp(:,:) + call evaluate_moist_adiabat_temperature( & + temp_aj, qsat_aj, is_converged_ref_profile, & ! (out) + tbase_mid, pres, zlev, adjust_mask, npz, nez ) ! (in) + + if ( .not. is_converged_ref_profile ) then + converged = .false. + exit + end if + + ! Set the adjusted vapor to stauration + + qvap_aj(:,:) = qvap(:,:) + where ( adjust_mask(:,:) ) + qvap_aj(:,:) = qsat_aj(:,:) + end where + + water_mass_aj = 0.0_RP + do ke_z=1, nez + do pz=1, npz + if ( adjust_mask(pz,ke_z) ) then + water_mass_aj = water_mass_aj + int_weight(pz,ke_z) * dens(pz,ke_z) * qvap_aj(pz,ke_z) + end if + end do + end do + + ! Reject the current profile becuase it requires more water than is available in the column. + if ( water_mass_aj > water_mass_ini ) then + ! write(*,*) "water_mass_ini=", water_mass_ini, "water_mass_aj=", water_mass_aj, "tbase_mid=", tbase_mid + tbase_hi = tbase_mid + cycle + end if + water_feasible = .true. + + call column_mse( energy_aj, & ! (out) + temp_aj, qvap_aj, zlev, dens, & ! (in) + int_weight, adjust_mask, npz, nez ) ! (in) + + residual = energy_aj - energy_target + residual_evaluated = .true. + ! write(*,*) "iter=", iter, "tbase=", tbase_lo, tbase_mid, tbase_hi,"energy_aj=", energy_aj, "energy_target=", energy_target, "residual=", residual + + if ( abs(residual) <= MCA_ENERGY_RTOL * max(abs(energy_target),1.0_RP) ) then + converged = .true. + exit + end if + + if ( residual > 0.0_RP ) then + tbase_hi = tbase_mid + else + tbase_lo = tbase_mid + end if + end do ! end loop for iteration + + if ( .not. converged ) then + write(*,*) "Error: Moist neutral profile construction did not converge." + if ( residual_evaluated ) then + write(*,*) " residual = ", residual + else if ( .not. water_feasible ) then + write(*,*) " No water-feasible saturated profile was evaluated." + end if + return + end if + + ! Convert the converged state to density variables + + rhoqvap_aj(:,:) = dens(:,:) * qvap_aj(:,:) + + precip_mass_total = max( 0.0_RP, water_mass_ini - water_mass_aj ) + positive_cond_mass = 0.0_RP + do ke_z=1, nez + do pz=1, npz + if ( adjust_mask(pz,ke_z) ) then + precip_mass_local = max( 0.0_RP, dens(pz,ke_z) * qvap(pz,ke_z) - rhoqvap_aj(pz,ke_z) ) + + positive_cond_mass = positive_cond_mass + precip_mass_local * int_weight(pz,ke_z) + end if + end do + end do + + ! Distribute precipitation mass over locally condensing nodes + + if ( positive_cond_mass > EPS ) then + precip_scale = precip_mass_total / positive_cond_mass + do ke_z=1, nez + do pz=1, npz + if ( adjust_mask(pz,ke_z) ) then + precip_mass_local = max( 0.0_RP, dens(pz,ke_z) * qvap(pz,ke_z) - rhoqvap_aj(pz,ke_z) ) + rhoprecip_aj(pz,ke_z) = precip_scale * precip_mass_local + end if + end do + end do + end if + + return + contains + ! The column MSE constraint is imposed on the fixed-density, vapor-only adjusted state. + ! Energy carried away by precipitation is neglected in the present simplified adjustment scheme. + subroutine column_mse( energy, & + temp_, qv_, zlev_, dens_, & + int_weight_, adjust_mask_, npz_, nez_ ) + implicit none + integer, intent(in) :: npz_ + integer, intent(in) :: nez_ + real(RP), intent(out) :: energy + real(RP), intent(in) :: temp_(npz_,nez_) + real(RP), intent(in) :: qv_(npz_,nez_) + real(RP), intent(in) :: zlev_(npz_,nez_) + real(RP), intent(in) :: dens_(npz_,nez_) + real(RP), intent(in) :: int_weight_(npz_,nez_) + logical, intent(in) :: adjust_mask_(npz_,nez_) + + integer :: ke_z_, pz_ + real(RP) :: qdry_ + real(RP) :: cptot_ + !------------------------------------------------ + + energy = 0.0_RP + + do ke_z_ = 1, nez_ + do pz_ = 1, npz_ + if ( adjust_mask_(pz_,ke_z_) ) then + qdry_ = 1.0_RP - qv_(pz_,ke_z_) !- qcon_(pz_,ke_z_) + + cptot_ = CPdry * qdry_ & + + CP_VAPOR * qv_(pz_,ke_z_) + ! Liquid-water sensible heat is neglected. + + energy = energy + int_weight_(pz_,ke_z_) * dens_(pz_,ke_z_) * & + ( cptot_ * temp_(pz_,ke_z_) & + + GRAV * zlev_(pz_,ke_z_) & + + LHV0 * qv_(pz_,ke_z_) ) + end if + end do + end do + return + end subroutine column_mse + end subroutine construct_moist_neutral_profile + +!OCL SERIAL + subroutine evaluate_moist_adiabat_temperature( & + temp_ref, qsat_ref, is_converged_ref_profile, & ! (out) + tbase, pres, zlev, mask, npz, nez ) ! (in) + implicit none + integer, intent(in) :: npz + integer, intent(in) :: nez + real(RP), intent(inout) :: temp_ref(npz,nez) + real(RP), intent(inout) :: qsat_ref(npz,nez) + logical, intent(out) :: is_converged_ref_profile + real(RP), intent(in) :: tbase + real(RP), intent(in) :: pres(npz,nez) + real(RP), intent(in) :: zlev(npz,nez) + logical, intent(in) :: mask(npz,nez) + + integer :: ke_z, pz + + real(RP) :: t_prev, p_prev, z_prev + real(RP) :: t_now, p_now, z_now + + logical :: started + logical :: is_converge_local + !---------------------------------------------------- + + started = .false. + is_converged_ref_profile = .true. + + do ke_z = 1, nez + do pz = 1, npz + if ( .not. mask(pz,ke_z) ) cycle + + z_now = zlev(pz,ke_z) + p_now = pres(pz,ke_z) + + if ( .not. started ) then + t_now = tbase + started = .true. + else if ( abs(z_now-z_prev) <= MCA_Z_TOL ) then + t_now = t_prev + else + call advance_moist_adiabat( t_now, is_converge_local, & ! (out) + t_prev, p_prev, z_prev, p_now, z_now ) ! (in) + + if ( .not. is_converge_local ) then + is_converged_ref_profile = .false. + return + end if + end if + + temp_ref(pz,ke_z) = t_now + call ATMOS_SATURATION_pres2qsat_liq( & + t_now, p_now, & ! (in) + qsat_ref(pz,ke_z) ) ! (out) + + t_prev = t_now + p_prev = p_now + z_prev = z_now + end do + end do + + return + end subroutine evaluate_moist_adiabat_temperature + +!OCL SERIAL + subroutine advance_moist_adiabat( temp1, converged, & ! (out) + temp0, pres0, zlev0, pres1, zlev1 ) ! (in) + implicit none + real(RP), intent(out) :: temp1 + logical, intent(out) :: converged + real(RP), intent(in) :: temp0 + real(RP), intent(in) :: pres0 + real(RP), intent(in) :: zlev0 + real(RP), intent(in) :: pres1 + real(RP), intent(in) :: zlev1 + integer :: iter + + real(RP) :: dz + real(RP) :: p_mid, t_mid + real(RP) :: tguess, tnew + real(RP) :: qsat_mid + real(RP) :: Gamma_m + real(RP) :: error + + ! Local temperature iteration tolerances + real(RP), parameter :: MCA_TEMP_ATOL = 1.0E-8_RP ! [K] + real(RP), parameter :: MCA_TEMP_RTOL = 1.0E-10_RP + real(RP), parameter :: MCA_PRES_RTOL = 1.0E-12_RP + !---------------------------------------------------- + + converged = .false. + dz = zlev1 - zlev0 + + if ( abs(dz) <= MCA_Z_TOL ) then + temp1 = temp0 + converged = .true. + return + end if + + !- Evaluate the midpoint pressure for the moist adiabatic lapse rate calculation + + if ( abs(pres1 - pres0) <= MCA_PRES_RTOL * max(pres0, pres1) ) then + p_mid = 0.5_RP * ( pres0 + pres1 ) + else + p_mid = ( pres1 - pres0 ) / log( pres1 / pres0 ) + end if + + !- Calculate an initial guess + + + call ATMOS_SATURATION_pres2qsat_liq( & + temp0, p_mid, & ! (in) + qsat_mid ) ! (out) + + call moist_adiabatic_lapse_rate( Gamma_m, & ! (out) + temp0, qsat_mid ) ! (in) + + tguess = temp0 - Gamma_m * dz + tnew = tguess + + !---------------------------------------------- + ! Iteration to solve the implicit equation for moist adiabatic lapse rate + ! T1^* = T0 - Gamma_m ( (T0 + T1^*)/2, pmid ) * dz + !---------------------------------------------- + + do iter = 1, MCA_MAX_LOCAL_ITER + + t_mid = 0.5_RP * ( temp0 + tguess ) + + call ATMOS_SATURATION_pres2qsat_liq( t_mid, p_mid, & ! (in) + qsat_mid ) ! (out) + + call moist_adiabatic_lapse_rate( Gamma_m, & ! (out) + t_mid, qsat_mid ) ! (in) + + tnew = temp0 - Gamma_m * dz + + error = abs(tnew - tguess) + if ( error <= MCA_TEMP_ATOL & + + MCA_TEMP_RTOL * max(abs(tnew), abs(tguess)) ) then + temp1 = tnew + converged = .true. + return + end if + + tguess = tnew + ! If you want to stabilize the fixed-point iteration, + ! tguess = ( 1.0_RP - MCA_LOCAL_RELAX ) * tguess + MCA_LOCAL_RELAX * tnew + end do + + temp1 = tnew + return + contains + subroutine moist_adiabatic_lapse_rate( GammaM, & + temp, qsat ) + implicit none + real(RP), intent(out) :: GammaM + real(RP), intent(in) :: temp + real(RP), intent(in) :: qsat + + real(RP) :: Rvap + real(RP) :: Rtot, Cptot + real(RP) :: EpsV + real(RP) :: fac, A + !------------------------------------------------ + Rvap = CP_VAPOR - CV_VAPOR + Rtot = Rdry * ( 1.0_RP - qsat ) + Rvap * qsat + Cptot = CPdry * ( 1.0_RP - qsat ) + CP_VAPOR * qsat + EpsV = Rdry / Rvap + fac = 1.0_RP + ( 1.0_RP - EpsV ) / EpsV * qsat + A = LHV0 + ( CP_VAPOR - CPdry ) * temp + + GammaM = GRAV & + * ( 1.0_RP + A * qsat * fac / ( Rtot * temp ) ) & + / ( Cptot + A * LHV0 * qsat * fac / ( Rvap * temp**2 ) ) + return + end subroutine moist_adiabatic_lapse_rate + end subroutine advance_moist_adiabat +end module scale_atm_phy_cp_dgm_mconv_adjustment + From e02eedad33e71b5ee387f7bbe5090232f811f4ea Mon Sep 17 00:00:00 2001 From: Yuta Kawai Date: Sun, 26 Jul 2026 14:49:22 +0900 Subject: [PATCH 4/9] Add moist convective adjustment scheme to cumulus parameterization --- FElib/src/Makefile | 5 + FElib/src/depend | 1 + .../src/admin/mod_dg_driver.F90 | 2 + .../src/atmos/mod_atmos_phy_cp.F90 | 132 +++++++++++++++++- .../src/atmos/mod_atmos_phy_cp_vars.F90 | 62 ++++++++ 5 files changed, 199 insertions(+), 3 deletions(-) diff --git a/FElib/src/Makefile b/FElib/src/Makefile index eab10590..c6a1053e 100644 --- a/FElib/src/Makefile +++ b/FElib/src/Makefile @@ -25,6 +25,7 @@ VPATH = \ turbulence: \ bl_turbulence: \ microphysics: \ + cumulus: \ radiation: \ model_framework: @@ -180,6 +181,9 @@ OBJS_NAME_MICROPHYS = \ scale_atm_phy_mp_dgm_common.o \ scale_atm_phy_mp_lscond.o +OBJS_NAME_CP = \ + scale_atm_phy_cp_dgm_mconv_adjustment.o + OBJS_NAME_RADIATION = \ scale_atm_phy_rd_solarins_simple.o \ scale_atm_phy_rd_dgm_simple.o \ @@ -204,6 +208,7 @@ OBJS_NAME = \ $(OBJS_NAME_TURBULENCE) \ $(OBJS_NAME_BL_TURBULENCE) \ $(OBJS_NAME_MICROPHYS) \ + $(OBJS_NAME_CP) \ $(OBJS_NAME_RADIATION) \ $(OBJS_NAME_MODEL_FRAMEWORK) diff --git a/FElib/src/depend b/FElib/src/depend index d6ac5b5a..8bfec80c 100644 --- a/FElib/src/depend +++ b/FElib/src/depend @@ -34,6 +34,7 @@ $(BUILD_DIR)/scale_atm_dyn_dgm_spongelayer.o: fluid_dyn_solver/scale_atm_dyn_dgm $(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_cp_dgm_mconv_adjustment.o: cumulus/scale_atm_phy_cp_dgm_mconv_adjustment.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_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/model/atm_nonhydro3d/src/admin/mod_dg_driver.F90 b/model/atm_nonhydro3d/src/admin/mod_dg_driver.F90 index 22ef3c06..f3075aa2 100644 --- a/model/atm_nonhydro3d/src/admin/mod_dg_driver.F90 +++ b/model/atm_nonhydro3d/src/admin/mod_dg_driver.F90 @@ -176,6 +176,7 @@ subroutine dg_driver( & 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 ( atmos%phy_cp_proc%IsActivated() ) call atmos%phy_cp_proc%vars%History() if ( ocean%IsActivated() ) call ocean%vars%History() @@ -365,6 +366,7 @@ subroutine restart_read() 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() + if ( atmos%phy_cp_proc%IsActivated() ) call atmos%phy_cp_proc%vars%History() call atmos%vars%Monitor() end if diff --git a/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_cp.F90 b/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_cp.F90 index e3ba96db..b1a16cf5 100644 --- a/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_cp.F90 +++ b/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_cp.F90 @@ -21,6 +21,8 @@ module mod_atmos_phy_cp 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 @@ -61,6 +63,9 @@ module mod_atmos_phy_cp integer :: atm_var_container_typeid !< Type ID of variable container for cumulus parameterization + real(RP) :: dtsec !< Timestep for cumulus parameterization + + type(LineElement) :: v_elem1D contains procedure :: setup => AtmosPhyCp_setup procedure :: calc_tendency => AtmosPhyCp_calc_tendency @@ -79,6 +84,8 @@ module mod_atmos_phy_cp ! !++ Private parameters & variables ! + integer, parameter :: CP_TYPEID_MCONV_ADJUSTMENT = 1 !< Type ID of a moist convective adjustment scheme + contains !> Setup a component of cumulus parameterization in atmospheric model @@ -92,6 +99,9 @@ subroutine AtmosPhyCp_setup( this, model_mesh, tm_parent_comp ) use mod_atmos_mesh, only: AtmosMesh use scale_time_manager, only: TIME_manager_component use mod_atmos_vars, only: ATM_VARS_CONTAINER_PRIMARY_ID + + use scale_atm_phy_cp_dgm_mconv_adjustment, only: & + atm_phy_cp_dgm_mconv_adjustment_setup implicit none class(AtmosPhyCp), intent(inout) :: this class(ModelMeshBase), target, intent(in) :: model_mesh @@ -150,9 +160,14 @@ subroutine AtmosPhyCp_setup( this, model_mesh, tm_parent_comp ) call tm_parent_comp%Regist_process( 'ATMOS_PHY_CP', 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 cumulus parameterization select case( CP_TYPE ) + case( 'MOIST_CONV_ADJUSTMENT' ) + this%CP_TYPEID = CP_TYPEID_MCONV_ADJUSTMENT + call atm_phy_cp_dgm_mconv_adjustment_setup() case default LOG_ERROR("ATMOS_PHY_CP_setup",*) 'Not appropriate cumulus parameterization type. Check!' call PRC_abort @@ -161,6 +176,9 @@ subroutine AtmosPhyCp_setup( this, model_mesh, tm_parent_comp ) !- Initialize the variables call this%vars%Init( model_mesh ) + !- + call this%v_elem1D%Init( atm_mesh%ptr_mesh%refElem3D%PolyOrder_v, .false. ) + return end subroutine AtmosPhyCp_setup @@ -177,6 +195,20 @@ end subroutine AtmosPhyCp_setup subroutine AtmosPhyCp_calc_tendency( & this, model_mesh, prgvars_list, trcvars_list, & auxvars_list, forcing_list, is_update ) + use scale_tracer, only: & + QA + use scale_atm_phy_cp_dgm_mconv_adjustment, only: & + atm_phy_cp_dgm_mconv_adjustment_calc_tendency + + use mod_atmos_vars, only: & + AtmosVars_GetLocalMeshPrgVars, & + AtmosVars_GetLocalMeshPhyAuxVars, & + AtmosVars_GetLocalMeshQTRCVarList, & + AtmosVars_GetLocalMeshPhyTends, & + AtmosVars_GetLocalMeshQTRC_Qv + use mod_atmos_phy_cp_vars, only: & + AtmosPhyCpVars_GetLocalMeshFields_tend, & + SFLX_RAIN_ID => ATMOS_PHY_CP_AUX2D_SFLX_RAIN_ID implicit none class(AtmosPhyCp), intent(inout) :: this class(ModelMeshBase), intent(in) :: model_mesh @@ -185,7 +217,97 @@ subroutine AtmosPhyCp_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, QV + 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 + class(LocalMeshFieldBase), pointer :: RHOQ_tp + class(LocalMeshFieldBase), pointer :: cp_DENS_t, cp_RHOT_t, cp_RHOQv_t !------------------------------------------------------------------------ + + if (.not. this%IsActivated()) return + + LOG_PROGRESS(*) 'atmosphere / physics / cumulus parameterization' + + call model_mesh%GetModelMesh( mesh ) + select type(mesh) + class is (MeshBase3D) + mesh3D => mesh + end select + + !- + if ( is_update ) then + call PROF_rapstart( 'ATM_CP_tendency', 2) + 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 AtmosVars_GetLocalMeshQTRC_Qv( n, & + mesh, trcvars_list, forcing_list, & + QV ) + + call AtmosPhyCpVars_GetLocalMeshFields_tend( n, & + mesh, this%vars%tends_manager, & + cp_DENS_t, cp_RHOT_t, cp_RHOQv_t ) + + select case( this%CP_TYPEID ) + case( CP_TYPEID_MCONV_ADJUSTMENT ) + call atm_phy_cp_dgm_mconv_adjustment_calc_tendency( & + cp_DENS_t%val, cp_RHOT_t%val, cp_RHOQv_t%val, & ! (out) + this%vars%auxvars2D(SFLX_RAIN_ID)%local(n)%val, & ! (out) + DDENS%val, DRHOT%val, QV%val, PT%val, PRES%val, & ! (in) + DENS_hyd%val, Rtot%val, CPtot%val, this%dtsec, & ! (in) + lcmesh, lcmesh%refElem3D, this%v_elem1D ) ! (in) + end select + end do + + call PROF_rapend( 'ATM_CP_tendency', 2) + end if + + call PROF_rapstart('ATM_PHY_CP_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 ) + + call AtmosVars_GetLocalMeshQTRC_Qv( n, & + mesh, trcvars_list, forcing_list, & + QV, RHOQ_tp ) + + call AtmosPhyCpVars_GetLocalMeshFields_tend( n, & + mesh, this%vars%tends_manager, & + cp_DENS_t, cp_RHOT_t, cp_RHOQv_t, & + lcmesh ) + + !$omp parallel private(ke) + !$omp do + do ke=lcmesh%NeS, lcmesh%NeE + DENS_tp%val(:,ke) = DENS_tp%val(:,ke) + cp_DENS_t%val(:,ke) + RHOT_tp%val(:,ke) = RHOT_tp%val(:,ke) + cp_RHOT_t%val(:,ke) + RHOQ_tp%val(:,ke) = RHOQ_tp%val(:,ke) + cp_RHOQv_t%val(:,ke) + end do + !$omp end parallel + end do + call PROF_rapend('ATM_PHY_CP_add_tend', 2) + return end subroutine AtmosPhyCp_calc_tendency @@ -218,17 +340,21 @@ end subroutine AtmosPhyCp_update !! !OCL SERIAL subroutine AtmosPhyCp_finalize( this ) + use scale_atm_phy_cp_dgm_mconv_adjustment, only: & + atm_phy_cp_dgm_mconv_adjustment_finalize implicit none class(AtmosPhyCp), intent(inout) :: this !-------------------------------------------------- if (.not. this%IsActivated()) return - ! select case ( this%CP_TYPEID ) - ! case( CP_TYPEID_LSCOND ) - ! end select + select case ( this%CP_TYPEID ) + case( CP_TYPEID_MCONV_ADJUSTMENT ) + call atm_phy_cp_dgm_mconv_adjustment_finalize() + end select call this%vars%Final() + call this%v_elem1D%Final() return end subroutine AtmosPhyCp_finalize diff --git a/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_cp_vars.F90 b/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_cp_vars.F90 index 71423294..090d7f1f 100644 --- a/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_cp_vars.F90 +++ b/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_cp_vars.F90 @@ -66,6 +66,8 @@ module mod_atmos_phy_cp_vars procedure :: History => AtmosPhyCpVars_history end type AtmosPhyCpVars + public :: AtmosPhyCpVars_GetLocalMeshFields_tend + !----------------------------------------------------------------------------- ! !++ Public variables @@ -223,12 +225,72 @@ subroutine AtmosPhyCpVars_Final( this ) return end subroutine AtmosPhyCpVars_Final +!OCL SERIAL + subroutine AtmosPhyCpVars_GetLocalMeshFields_tend( domID, mesh, bl_tends_list, & + cp_DENS_t, cp_RHOT_t, cp_RHOQv_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) :: cp_DENS_t + class(LocalMeshFieldBase), pointer, intent(out) :: cp_RHOT_t + class(LocalMeshFieldBase), pointer, intent(out) :: cp_RHOQv_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_CP_DENS_t_ID, field) + call field%GetLocalMeshField(domID, cp_DENS_t) + + call bl_tends_list%Get(ATMOS_PHY_CP_RHOT_t_ID, field) + call field%GetLocalMeshField(domID, cp_RHOT_t) + + call bl_tends_list%Get(ATMOS_PHY_CP_RHOQv_t_ID, field) + call field%GetLocalMeshField(domID, cp_RHOQv_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 AtmosPhyCpVars_GetLocalMeshFields_tend + !OCL SERIAL subroutine AtmosPhyCpVars_history( this ) use scale_file_history_meshfield, only: FILE_HISTORY_meshfield_put implicit none class(AtmosPhyCpVars), 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_CP_AUX2D_NUM + hst_id = this%auxvars2D(v)%hist_id + if ( hst_id > 0 ) call FILE_HISTORY_meshfield_put( hst_id, this%auxvars2D(v) ) + end do return end subroutine AtmosPhyCpVars_history From f8c51051b166a369ba16eb82752ee59ee7ee393d Mon Sep 17 00:00:00 2001 From: Yuta Kawai Date: Sun, 26 Jul 2026 16:48:29 +0900 Subject: [PATCH 5/9] Modify setup procedures of variables for various atmospheric components --- .../src/atmos/mod_atmos_component.F90 | 35 +++++++--- .../src/atmos/mod_atmos_phy_bl_vars.F90 | 37 ++++++---- .../src/atmos/mod_atmos_phy_cp_vars.F90 | 16 ++++- .../src/atmos/mod_atmos_phy_mp_vars.F90 | 32 +++++---- .../src/atmos/mod_atmos_phy_rd_vars.F90 | 30 +++++++- .../src/atmos/mod_atmos_phy_sfc.F90 | 5 +- .../src/atmos/mod_atmos_phy_sfc_vars.F90 | 68 +++++++------------ .../src/atmos/mod_atmos_phy_tb_vars.F90 | 31 +++++---- 8 files changed, 156 insertions(+), 98 deletions(-) diff --git a/model/atm_nonhydro3d/src/atmos/mod_atmos_component.F90 b/model/atm_nonhydro3d/src/atmos/mod_atmos_component.F90 index b199c443..1c28ecf1 100644 --- a/model/atm_nonhydro3d/src/atmos/mod_atmos_component.F90 +++ b/model/atm_nonhydro3d/src/atmos/mod_atmos_component.F90 @@ -274,27 +274,46 @@ subroutine Atmos_setup_vars( this ) call PROF_rapstart( 'ATM_setup_vars', 1) + LOG_INFO('AtmosComponent_setup_vars',*) 'Atmosphere model components ' + call this%vars%Init( this%mesh ) + !- Cloud microphysics component if ( this%phy_mp_proc%IsActivated() ) then - call this%vars%Regist_physvar_manager( & - this%phy_mp_proc%vars%auxvars2D_manager ) + call this%phy_mp_proc%vars%Setup( this%mesh ) + call this%phy_mp_proc%Set_primary_atmvars_container( this%vars%container ) + call this%vars%Regist_physvar_manager( this%phy_mp_proc%vars%auxvars2D_manager ) call this%vars%Setup_container( this%phy_mp_proc%atm_var_container_typeid, this%mesh ) - call this%phy_mp_proc%Set_primary_atmvars_container( this%vars%container ) end if + !- Surface component if ( this%phy_sfc_proc%IsActivated() ) then + call this%phy_sfc_proc%vars%Setup( this%mesh ) + call this%vars%Setup_container( this%phy_sfc_proc%atm_var_container_typeid, this%mesh ) end if + !- Turbulence component + if ( this%phy_tb_proc%IsActivated() ) then + call this%phy_tb_proc%vars%Setup( this%mesh ) + end if + !- Radiation component if ( this%phy_rd_proc%IsActivated() ) then - if ( .not. this%phy_sfc_proc%IsActivated() ) then - LOG_ERROR('ATM_setup_vars',*) 'ATMOS_PHY_RD_DO requires ATMOS_PHY_SF_DO to provide SFC_TEMP.' - call PRC_abort - end if + if ( .not. this%phy_sfc_proc%IsActivated() ) then + LOG_ERROR('ATM_setup_vars',*) 'ATMOS_PHY_RD_DO requires ATMOS_PHY_SF_DO to provide SFC_TEMP.' + call PRC_abort + end if + call this%phy_rd_proc%vars%Setup( this%mesh ) call this%phy_rd_proc%SetSfcVars( this%phy_sfc_proc%vars%SFC_VARS(SFCTEMP_ID), & this%phy_sfc_proc%vars%SFC_VARS(SFC_ALB_ID) ) end if - + !- Cumulus parameterization component + if ( this%phy_cp_proc%IsActivated() ) then + call this%phy_cp_proc%vars%Setup( this%mesh ) + end if + !- PBL component + if ( this%phy_bl_proc%IsActivated() ) then + call this%phy_bl_proc%vars%Setup( this%mesh ) + end if call PROF_rapend( 'ATM_setup_vars', 1) return 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 e7e97be2..b6c613a5 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 @@ -67,6 +67,7 @@ module mod_atmos_phy_bl_vars contains procedure :: Init => AtmosPhyBlVars_Init procedure :: Final => AtmosPhyBlVars_Final + procedure :: Setup => AtmosPhyBlVars_Setup procedure :: History => AtmosPhyBlVars_history end type AtmosPhyBlVars @@ -113,18 +114,35 @@ module mod_atmos_phy_bl_vars !OCL SERIAL subroutine AtmosPhyBlVars_Init( this, model_mesh, & QS_BL, QE_BL, QA_BL ) + implicit none + class(AtmosPhyBlVars), target, intent(inout) :: this + class(ModelMeshBase), target, intent(in) :: model_mesh + integer, intent(in) :: QS_BL + integer, intent(in) :: QE_BL + integer, intent(in) :: QA_BL + !---------------------------------------------------- + + LOG_INFO('AtmosPhyBlVars_Init',*) + + this%QS = QS_BL + this%QE = QE_BL + this%QA = QA_BL + this%TENDS_NUM_TOT = ATMOS_PHY_BL_TENDS_NUM1 + QE_BL - QS_BL + 1 + + return + end subroutine AtmosPhyBlVars_Init + !> Setup variable objects +!OCL SERIAL + subroutine AtmosPhyBlVars_Setup( this, model_mesh ) use scale_tracer, only: & - TRACER_NAME, TRACER_DESC, TRACER_UNIT + TRACER_NAME, TRACER_DESC, TRACER_UNIT, & + QA use scale_file_history, only: & FILE_HISTORY_reg - implicit none class(AtmosPhyBlVars), target, intent(inout) :: this class(ModelMeshBase), target, intent(in) :: model_mesh - integer, intent(in) :: QS_BL - integer, intent(in) :: QE_BL - integer, intent(in) :: QA_BL integer :: iv integer :: iq @@ -139,13 +157,6 @@ subroutine AtmosPhyBlVars_Init( this, model_mesh, & type(VariableInfo) :: qtrc_vterm_vinfo_tmp !---------------------------------------------------- - LOG_INFO('AtmosPhyBlVars_Init',*) - - this%QS = QS_BL - this%QE = QE_BL - this%QA = QA_BL - this%TENDS_NUM_TOT = ATMOS_PHY_BL_TENDS_NUM1 + QE_BL - QS_BL + 1 - !- Initialize auxiliary and diagnostic variables nullify( atm_mesh ) @@ -199,7 +210,7 @@ subroutine AtmosPhyBlVars_Init( this, model_mesh, & end do return - end subroutine AtmosPhyBlVars_Init + end subroutine AtmosPhyBlVars_Setup !> Finalize an object to manage variables with planetary boundary layer (PBL) turbulence parameterization component !OCL SERIAL diff --git a/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_cp_vars.F90 b/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_cp_vars.F90 index 090d7f1f..4ce986a1 100644 --- a/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_cp_vars.F90 +++ b/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_cp_vars.F90 @@ -63,6 +63,7 @@ module mod_atmos_phy_cp_vars contains procedure :: Init => AtmosPhyCpVars_Init procedure :: Final => AtmosPhyCpVars_Final + procedure :: Setup => AtmosPhyCpVars_Setup procedure :: History => AtmosPhyCpVars_history end type AtmosPhyCpVars @@ -111,6 +112,17 @@ module mod_atmos_phy_cp_vars !> Setup an object to manage variables with a cumulus parameterization component !OCL SERIAL subroutine AtmosPhyCpVars_Init( this, model_mesh ) + implicit none + class(AtmosPhyCpVars), target, intent(inout) :: this + class(ModelMeshBase), target, intent(in) :: model_mesh + !---------------------------------------------------- + + LOG_INFO('AtmosPhyCpVars_Init',*) + return + end subroutine AtmosPhyCpVars_Init + + !> Setup variable objects with cumulus parameterization component + subroutine AtmosPhyCpVars_Setup( this, model_mesh ) use scale_atmos_hydrometeor, only: & N_HYD, & HYD_NAME @@ -118,7 +130,6 @@ subroutine AtmosPhyCpVars_Init( this, model_mesh ) TRACER_NAME, TRACER_DESC, TRACER_UNIT use scale_file_history, only: & FILE_HISTORY_reg - implicit none class(AtmosPhyCpVars), target, intent(inout) :: this class(ModelMeshBase), target, intent(in) :: model_mesh @@ -136,7 +147,6 @@ subroutine AtmosPhyCpVars_Init( this, model_mesh ) type(VariableInfo) :: qtrc_vterm_vinfo_tmp !---------------------------------------------------- - LOG_INFO('AtmosPhyCpVars_Init',*) this%TENDS_NUM_TOT = ATMOS_PHY_CP_TENDS_NUM1 + N_HYD @@ -205,7 +215,7 @@ subroutine AtmosPhyCpVars_Init( this, model_mesh ) end do return - end subroutine AtmosPhyCpVars_Init + end subroutine AtmosPhyCpVars_Setup !> Finalize an object to manage variables with cumulus parameterization component !OCL SERIAL diff --git a/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_mp_vars.F90 b/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_mp_vars.F90 index 08d37933..eb0436e1 100644 --- a/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_mp_vars.F90 +++ b/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_mp_vars.F90 @@ -69,6 +69,7 @@ module mod_atmos_phy_mp_vars contains procedure :: Init => AtmosPhyMpVars_Init procedure :: Final => AtmosPhyMpVars_Final + procedure :: Setup => AtmosPhyMpVars_Setup procedure :: History => AtmosPhyMpVars_history end type AtmosPhyMpVars @@ -134,18 +135,32 @@ module mod_atmos_phy_mp_vars !OCL SERIAL subroutine AtmosPhyMpVars_Init( this, model_mesh, & QS_MP, QE_MP, QA_MP ) + implicit none + class(AtmosPhyMpVars), target, intent(inout) :: this + class(ModelMeshBase), target, intent(in) :: model_mesh + integer, intent(in) :: QS_MP + integer, intent(in) :: QE_MP + integer, intent(in) :: QA_MP + !-------------------------------------------------- + + LOG_INFO('AtmosPhyMpVars_Init',*) + + this%QS = QS_MP + this%QE = QE_MP + this%QA = QA_MP + this%TENDS_NUM_TOT = ATMOS_PHY_MP_TENDS_NUM1 + QE_MP - QS_MP + 1 + return + end subroutine AtmosPhyMpVars_Init + !> Setup variable objects with cloud microphysics component + subroutine AtmosPhyMpVars_Setup( this, model_mesh ) use scale_tracer, only: & TRACER_NAME, TRACER_DESC, TRACER_UNIT use scale_file_history, only: & FILE_HISTORY_reg - implicit none class(AtmosPhyMpVars), target, intent(inout) :: this class(ModelMeshBase), target, intent(in) :: model_mesh - integer, intent(in) :: QS_MP - integer, intent(in) :: QE_MP - integer, intent(in) :: QA_MP integer :: iv integer :: iq @@ -160,13 +175,6 @@ subroutine AtmosPhyMpVars_Init( this, model_mesh, & type(VariableInfo) :: qtrc_vterm_vinfo_tmp !-------------------------------------------------- - LOG_INFO('AtmosPhyMpVars_Init',*) - - this%QS = QS_MP - this%QE = QE_MP - this%QA = QA_MP - this%TENDS_NUM_TOT = ATMOS_PHY_MP_TENDS_NUM1 + QE_MP - QS_MP + 1 - !- Initialize auxiliary and diagnostic variables nullify( atm_mesh ) @@ -250,7 +258,7 @@ subroutine AtmosPhyMpVars_Init( this, model_mesh, & end do return - end subroutine AtmosPhyMpVars_Init + end subroutine AtmosPhyMpVars_Setup !> Finalize an object to manage variables with cloud microphysics component !OCL SERIAL diff --git a/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_rd_vars.F90 b/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_rd_vars.F90 index c94bcdab..570a5fe9 100644 --- a/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_rd_vars.F90 +++ b/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_rd_vars.F90 @@ -67,6 +67,7 @@ module mod_atmos_phy_rd_vars contains procedure :: Init => AtmosPhyRdVars_Init procedure :: Final => AtmosPhyRdVars_Final + procedure :: Setup => AtmosPhyRdVars_Setup procedure :: History => AtmosPhyRdVars_history procedure :: Read_restart_file => AtmosPhyRdVars_Read_restart_file procedure :: Write_restart_file_prep => AtmosPhyRdVars_Write_restart_file_prep @@ -150,6 +151,31 @@ subroutine AtmosPhyRdVars_Init( this, model_mesh ) logical :: reg_file_hist !---------------------------------------------------- + LOG_INFO('AtmosPhyRdVars_Init',*) + return + end subroutine AtmosPhyRdVars_Init + + !> Setup variables objects with radiation component + subroutine AtmosPhyRdVars_Setup( this, model_mesh ) + use scale_tracer, only: & + TRACER_NAME, TRACER_DESC, TRACER_UNIT + use scale_file_history, only: & + FILE_HISTORY_reg + + implicit none + class(AtmosPhyRdVars), target, intent(inout) :: this + class(ModelMeshBase), target, intent(in) :: model_mesh + + class(AtmosMesh), pointer :: atm_mesh + class(MeshBase2D), pointer :: mesh2D + class(MeshBase3D), pointer :: mesh3D + + integer :: iv + integer :: iq + integer :: n + logical :: reg_file_hist + !---------------------------------------------------- + LOG_INFO('AtmosPhyRdVars_Init',*) this%TENDS_NUM_TOT = ATMOS_PHY_RD_TENDS_NUM1 @@ -200,9 +226,9 @@ subroutine AtmosPhyRdVars_Init( this, model_mesh ) !- Setup restart file for variables with radiation component call atm_mesh%Setup_restartfile( this%restart_file, & ATMOS_PHY_RD_RESTART_VAR_NUM ) - + return - end subroutine AtmosPhyRdVars_Init + end subroutine AtmosPhyRdVars_Setup !> Finalize an object to manage variables with radiation component !OCL SERIAL diff --git a/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_sfc.F90 b/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_sfc.F90 index e76d2767..d4bd8ab0 100644 --- a/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_sfc.F90 +++ b/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_sfc.F90 @@ -161,7 +161,7 @@ subroutine AtmosPhySfc_setup( this, model_mesh, tm_parent_comp ) this%tm_process_id ) ! (out) !- initialize the variables - call this%vars%Init( model_mesh ) + call this%vars%Init( model_mesh, DEFAULT_SFC_TEMP ) !--- Set the type of surface flux scheme @@ -183,9 +183,6 @@ subroutine AtmosPhySfc_setup( this, model_mesh, tm_parent_comp ) this%mesh => model_mesh end select - !-- Set default values - call this%vars%SetDefaultVal( DEFAULT_SFC_TEMP ) - return end subroutine AtmosPhySfc_setup diff --git a/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_sfc_vars.F90 b/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_sfc_vars.F90 index 9b4a1776..d19d57f2 100644 --- a/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_sfc_vars.F90 +++ b/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_sfc_vars.F90 @@ -62,8 +62,8 @@ module mod_atmos_phy_sfc_vars contains procedure :: Init => AtmosPhySfcVars_Init procedure :: Final => AtmosPhySfcVars_Final + procedure :: Setup => AtmosPhySfcVars_Setup procedure :: History => AtmosPhySfcVars_history - procedure :: SetDefaultVal => AtmosPhySfcVars_set_default_val end type AtmosPhySfcVars ! Variable information for surface variables @@ -119,10 +119,12 @@ module mod_atmos_phy_sfc_vars !> Setup an object to manage variables with a surface component !! !! @param model_mesh Object to manage computational mesh of atmospheric model - subroutine AtmosPhySfcVars_Init( this, model_mesh ) + subroutine AtmosPhySfcVars_Init( this, model_mesh, & + DEFAULT_SFC_TEMP ) implicit none class(AtmosPhySfcVars), target, intent(inout) :: this class(ModelMeshBase), target, intent(in) :: model_mesh + real(RP), intent(in) :: DEFAULT_SFC_TEMP integer :: v integer :: n @@ -131,27 +133,31 @@ subroutine AtmosPhySfcVars_Init( this, model_mesh ) class(AtmosMesh), pointer :: atm_mesh class(MeshBase2D), pointer :: mesh2D class(MeshBase3D), pointer :: mesh3D - - namelist / PARAM_ATMOS_PHY_SFC_VARS / & - ATMOS_PHY_SFC_DEFAULT_SFC_TEMP - - integer :: ierr !-------------------------------------------------- LOG_INFO('AtmosPhySfcVars_Init',*) + ATMOS_PHY_SFC_DEFAULT_SFC_TEMP = DEFAULT_SFC_TEMP + return + end subroutine AtmosPhySfcVars_Init + + !> Setup variable objects with a surface component + subroutine AtmosPhySfcVars_Setup( this, model_mesh ) + implicit none + class(AtmosPhySfcVars), target, intent(inout) :: this + class(ModelMeshBase), target, intent(in) :: model_mesh + + integer :: v + integer :: n + logical :: reg_file_hist + + class(AtmosMesh), pointer :: atm_mesh + class(MeshBase2D), pointer :: mesh2D + class(MeshBase3D), pointer :: mesh3D + !-------------------------------------------------- + this%SFCFLX_NUM_TOT = ATMOS_PHY_SF_SFLX_NUM1 - !--- read namelist - rewind(IO_FID_CONF) - read(IO_FID_CONF,nml=PARAM_ATMOS_PHY_SFC_VARS,iostat=ierr) - if( ierr < 0 ) then !--- missing - LOG_INFO("ATMOS_PHY_SFC_vars_setup",*) 'Not found namelist. Default used.' - elseif( ierr > 0 ) then !--- fatal error - LOG_ERROR("ATMOS_PHY_SFC_vars_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_SFC_VARS. Check!' - call PRC_abort - endif - LOG_NML(PARAM_ATMOS_PHY_SFC_VARS) !- Initialize auxiliary and diagnostic variables @@ -191,7 +197,7 @@ subroutine AtmosPhySfcVars_Init( this, model_mesh ) end do return - end subroutine AtmosPhySfcVars_Init + end subroutine AtmosPhySfcVars_Setup !> Finalize an object to manage variables with a surface component !! @@ -212,32 +218,6 @@ subroutine AtmosPhySfcVars_Final( this ) return end subroutine AtmosPhySfcVars_Final - !> Set default value for surface variables -!OCL SERIAL - subroutine AtmosPhySfcVars_set_default_val( this, & - SFC_TEMP ) - implicit none - class(AtmosPhySfcVars), intent(inout), target :: this - real(RP), intent(in) :: SFC_TEMP - - integer :: ldomID - class(MeshBase2D), pointer :: mesh2D - class(LocalMesh2D), pointer :: lcmesh2D - integer :: ke - !-------------------------------------------------- - - mesh2D => this%SFC_VARS(ATMOS_PHY_SF_SVAR_TEMP_ID)%mesh - - do ldomID=1, mesh2D%LOCAL_MESH_NUM - lcmesh2D => mesh2D%lcmesh_list(ldomID) - !$omp parallel do - do ke=lcmesh2D%NeS, lcmesh2D%NeE - this%SFC_VARS(ATMOS_PHY_SF_SVAR_TEMP_ID)%local(ldomID)%val(:,ke) = SFC_TEMP - end do - end do - return - end subroutine AtmosPhySfcVars_set_default_val - !> Write history data for surface variables !OCL SERIAL subroutine AtmosPhySfcVars_history( this ) diff --git a/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_tb_vars.F90 b/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_tb_vars.F90 index fc0edd45..7b16adc0 100644 --- a/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_tb_vars.F90 +++ b/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_tb_vars.F90 @@ -77,6 +77,7 @@ module mod_atmos_phy_tb_vars contains procedure :: Init => AtmosPhyTbVars_Init procedure :: Final => AtmosPhyTbVars_Final + procedure :: Setup => AtmosPhyTbVars_Setup procedure :: History => AtmosPhyTbVars_history end type AtmosPhyTbVars @@ -93,16 +94,24 @@ module mod_atmos_phy_tb_vars !! !! @param model_mesh Object to manage computational mesh of atmospheric model !OCL SERIAL - subroutine AtmosPhyTbVars_Init( this, model_mesh ) + subroutine AtmosPhyTbVars_Init( this, model_mesh ) + implicit none + class(AtmosPhyTbVars), target, intent(inout) :: this + class(ModelMeshBase), target, intent(in) :: model_mesh + !-------------------------------------------------- - use scale_tracer, only: & - TRACER_NAME, TRACER_DESC, TRACER_UNIT + LOG_INFO('AtmosPhyTbVars_Init',*) + return + end subroutine AtmosPhyTbVars_Init + + !> Setup variables objects +!OCL SERIAL + subroutine AtmosPhyTbVars_Setup( this, model_mesh ) use scale_atm_phy_tb_dgm_common, only: & atm_phy_tb_dgm_common_setup_variables - implicit none - class(AtmosPhyTbVars), target, intent(inout) :: this - class(ModelMeshBase), target, intent(in) :: model_mesh + class(AtmosPhyTbVars), intent(inout) :: this + class(ModelMeshBase), intent(in), target :: model_mesh integer :: iv integer :: iq @@ -110,13 +119,11 @@ subroutine AtmosPhyTbVars_Init( this, model_mesh ) class(AtmosMesh), pointer :: atm_mesh class(MeshBase2D), pointer :: mesh2D - class(MeshBase3D), pointer :: mesh3D + class(MeshBase3D), pointer :: mesh3D !-------------------------------------------------- - LOG_INFO('AtmosPhyTbVars_Init',*) - this%TENDS_NUM_TOT = ATMOS_PHY_TB_TENDS_NUM1 + QA - + !- Initialize auxiliary and diagnostic variables nullify( atm_mesh ) @@ -154,9 +161,9 @@ subroutine AtmosPhyTbVars_Init( this, model_mesh ) this%auxvars_manager, & ! (inout) this%auxvars(:), & ! (in) this%auxvars_commid ) ! (out) - + return - end subroutine AtmosPhyTbVars_Init + end subroutine AtmosPhyTbVars_Setup !> Finalize an object to manage variables with atmospheric SGS turbulent component !OCL SERIAL From 1d9cc2c03873bf236faff3952b3e6af9d5e4727d Mon Sep 17 00:00:00 2001 From: Yuta Kawai Date: Mon, 27 Jul 2026 11:38:06 +0900 Subject: [PATCH 6/9] Enhance PBL turbulence model to consider the vertical diffusion of tracers --- .../scale_atm_phy_bl_dgm_common.F90 | 460 ++++++++++++------ .../scale_atm_phy_bl_dgm_mynn_lv2.F90 | 76 +-- .../src/atmos/mod_atmos_phy_bl.F90 | 51 +- .../src/atmos/mod_atmos_phy_bl_vars.F90 | 18 +- 4 files changed, 392 insertions(+), 213 deletions(-) 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 index 63ddfc04..6729cbef 100644 --- a/FElib/src/bl_turbulence/scale_atm_phy_bl_dgm_common.F90 +++ b/FElib/src/bl_turbulence/scale_atm_phy_bl_dgm_common.F90 @@ -18,6 +18,8 @@ module scale_atm_phy_bl_dgm_common use scale_prof use scale_const, only: & EPS => CONST_EPS + use scale_tracer, only: QA + use scale_sparsemat use scale_element_base, only: & ElementBase1D, ElementBase2D, ElementBase3D @@ -25,7 +27,7 @@ module scale_atm_phy_bl_dgm_common 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_localmeshfield_base, only: LocalMeshField3D, LocalMeshFieldBaseList use scale_meshfield_base, only: MeshField3D use scale_element_operation_base, only: ElementOperationBase3D @@ -49,10 +51,11 @@ module scale_atm_phy_bl_dgm_common !++ Private parameters & variables ! contains + !> Calculate tendency with PBL turbulence models !OCL SERIAL subroutine atm_phy_bl_dgm_common_calc_tendency( & - RHOU_tp, RHOV_tp, DRHOT_tp, & ! (out) - DDENS_, MOMX_, MOMY_, DRHOT_, & ! (in) + RHOU_tp, RHOV_tp, DRHOT_tp, RHOQ_tp_list, & ! (out) + DDENS_, MOMX_, MOMY_, DRHOT_, RHOQ_list, & ! (in) PT_, DENS_hyd, PRES_hyd, NU, KH, & ! (in) element3D_operation, C_IP, dtsec, & ! (in) lmesh, elem, elem1D, is_bound, & ! (in) @@ -66,10 +69,12 @@ subroutine atm_phy_bl_dgm_common_calc_tendency( & 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) + type(LocalMeshFieldBaseList), intent(inout) :: RHOQ_tp_list(QA) 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) + type(LocalMeshFieldBaseList), intent(in) :: RHOQ_list(QA) 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) @@ -84,7 +89,10 @@ subroutine atm_phy_bl_dgm_common_calc_tendency( & class(LocalMesh2D), pointer :: lmesh2D class(ElementBase2D), pointer :: elem2D - real(RP) :: PROG_VARS (elem%Np,lmesh%NeX*lmesh%NeY,lmesh%NeZ,3) + integer :: iq + real(RP) :: RHOQ00_(elem%Np,QA,lmesh%Ne) + + real(RP) :: PROG_VARS (elem%Np,lmesh%NeX*lmesh%NeY,lmesh%NeZ,3+QA) real(RP) :: alph_M(elem%NfpTot,lmesh%Ne) real(RP) :: alph_H(elem%NfpTot,lmesh%Ne) real(RP) :: GsqrtV(elem%Np,lmesh%Ne) @@ -102,7 +110,7 @@ subroutine atm_phy_bl_dgm_common_calc_tendency( & real(RP), allocatable :: BndMatD(:,:,:,:,:) real(RP), allocatable :: G(:,:,:,:,:,:) - real(RP) :: impl_fac + real(RP) :: impl_fac, r_impl_fac !------------------------------------------------------------------------ lmesh2D => lmesh%lcmesh2D @@ -113,12 +121,21 @@ subroutine atm_phy_bl_dgm_common_calc_tendency( & 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( b1D_ij(im*elem%Nnode_v,3+QA,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 ) + !$omp parallel private( ke, ke2D ) + !$omp do collapse(2) + do iq = 1, QA + do ke=lmesh%NeS, lmesh%NeE + RHOQ00_(:,iq,ke) = RHOQ_list(iq)%ptr%val(:,ke) + end do + end do + !$omp end do + + !$omp do collapse(2) do ke_z =1, lmesh%NeZ do ke_xy=1, lmesh%NeX * lmesh%NeY ke = ke_xy + (ke_z-1)*lmesh%Ne2D @@ -129,17 +146,22 @@ subroutine atm_phy_bl_dgm_common_calc_tendency( & 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 iq = 1, QA + PROG_VARS(:,ke_xy,ke_z,3+iq) = RHOQ00_(:,iq,ke) + end do do p=1, elem%Np GsqrtV(p,ke) = lmesh%Gsqrt(p,ke) / lmesh%GsqrtH(elem%IndexH2Dto3D(p),ke2D) end do end do end do + !$omp end do + !$omp end parallel - 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 eval_Ax( RHOU_tp, RHOV_tp, DRHOT_tp, RHOQ_tp_list, alph_M, alph_H, & !(out) + PROG_VARS, MOMX_, MOMY_, PT_, RHOQ00_, 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) @@ -147,13 +169,19 @@ subroutine atm_phy_bl_dgm_common_calc_tendency( & im, jm, lmesh, elem, elem1D, use_delta_form ) ! (in) !--- - !$omp parallel do collapse(2) private( ke ) + r_impl_fac = 1.0_RP / impl_fac + + !$omp parallel do collapse(2) private( ke, iq ) 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 + RHOU_tp (:,ke) = ( PROG_VARS(:,ke_xy,ke_z,1) - MOMX_ (:,ke) ) * r_impl_fac + RHOV_tp (:,ke) = ( PROG_VARS(:,ke_xy,ke_z,2) - MOMY_ (:,ke) ) * r_impl_fac + DRHOT_tp(:,ke) = ( PROG_VARS(:,ke_xy,ke_z,3) - DENS(:,ke) * PT_(:,ke) ) * r_impl_fac + + do iq=1, QA + RHOQ_tp_list(iq)%ptr%val(:,ke) = ( PROG_VARS(:,ke_xy,ke_z,3+iq) - RHOQ00_(:,iq,ke) ) * r_impl_fac + end do end do end do @@ -162,6 +190,7 @@ end subroutine atm_phy_bl_dgm_common_calc_tendency !- Private subroutines ----------------------------- + !> Solve the block tridiagonal system which is generated by the vertical diffusion equation for PBL turbulence models !OCL SERIAL subroutine vi_solve( PROG_VARS, & ! (inout) BndMatL, BndMatD, G, b1D_ij, & ! (inout) @@ -172,11 +201,11 @@ subroutine vi_solve( PROG_VARS, & ! (inout) 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) :: PROG_VARS(elem%Nnode_h1D**2,elem%Nnode_v,lmesh%Ne2D,lmesh%NeZ,3+QA) 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(inout) :: b1D_ij(im*elem%Nnode_v,3+QA,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) @@ -189,9 +218,11 @@ subroutine vi_solve( PROG_VARS, & ! (inout) integer :: ke_z, ke_xy integer :: i, j, ij integer :: pv, pv1, pv2, pp, p2 + integer :: iv real(RP) :: tmp(im*elem%Nnode_v,2) - real(RP) :: tmp_b(im*elem%Nnode_v,3) + real(RP) :: tmp_b(im*elem%Nnode_v,2) + real(RP) :: tmp_b2(im*elem%Nnode_v) !------------------------------------------------------- do ke_z=1, lmesh%NeZ @@ -206,7 +237,7 @@ subroutine vi_solve( PROG_VARS, & ! (inout) 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 ) + !$omp parallel do collapse(2) private( pp, p2, tmp, tmp_b,tmp_b2 ) do ke_xy=1, lmesh%NeX * lmesh%NeY do j=1, jm !* D_k <- D_k - L_k * G_{k-1} ------------------ @@ -233,13 +264,25 @@ subroutine vi_solve( PROG_VARS, & ! (inout) 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) + + do iv=3, 3+QA + tmp_b2(:) = 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_b2(pp) = tmp_b2(pp) + BndMatL(pp,pv,j,ke_xy,2) * b1D_ij(p2,iv,j,ke_xy,ke_z-1) + end do + end do + end do + b1D_ij(:,iv,j,ke_xy,ke_z) = b1D_ij(:,iv,j,ke_xy,ke_z) - tmp_b2(:) + end do + end do end do end if @@ -249,7 +292,7 @@ subroutine vi_solve( PROG_VARS, & ! (inout) end do do ke_z=lmesh%NeZ-1, 1, -1 - !$omp parallel do collapse(2) private( pp, p2, tmp_b ) + !$omp parallel do collapse(2) private( pp, p2, tmp_b,tmp_b2 ) do ke_xy=1, lmesh%NeX * lmesh%NeY do j=1, jm tmp_b(:,:) = 0.0_RP @@ -259,13 +302,25 @@ subroutine vi_solve( PROG_VARS, & ! (inout) 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) + + do iv=3, 3+QA + tmp_b2(:) = 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_b2(pp) = tmp_b2(pp) + G(pp,pv2,j,ke_xy,ke_z,2) * b1D_ij(p2,iv,j,ke_xy,ke_z+1) + end do + end do + end do + b1D_ij(:,iv,j,ke_xy,ke_z) = b1D_ij(:,iv,j,ke_xy,ke_z) - tmp_b2(:) + end do + end do end do end do @@ -274,24 +329,24 @@ subroutine vi_solve( PROG_VARS, & ! (inout) 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 + if ( use_delta_form ) then + do iv=1, 3+QA + 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,iv) = PROG_VARS(pp,pv,ke_xy,ke_z,iv) + b1D_ij(p2,iv,j,ke_xy,ke_z) + end do + 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 + do iv=1, 3+QA + 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,iv) = b1D_ij(p2,iv,j,ke_xy,ke_z) + end do + end do + end do end if end do end do @@ -300,8 +355,9 @@ subroutine vi_solve( PROG_VARS, & ! (inout) return end subroutine vi_solve + !> Solve the linear system associated with the diagonal block of the block tridiagonal system !OCL SERIAL - subroutine linkernel_solve_sip( b, G1, G2, & ! (inout) + subroutine linkernel_solve_sip( b, G1, G2, & ! (inout) BndMatD1, BndMatD2, nv, im, jm, Ne2D, top_flag ) ! (in) implicit none @@ -309,7 +365,7 @@ subroutine linkernel_solve_sip( b, G1, G2, & ! (inout) 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) :: b(im*nv,3+QA,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) @@ -329,7 +385,7 @@ subroutine linkernel_solve_sip( b, G1, G2, & ! (inout) real(RP) :: Amat1(nv,nv) real(RP) :: RHS1(nv,nv+2) real(RP) :: Amat2(nv,nv) - real(RP) :: RHS2(nv,nv+2) + real(RP) :: RHS2(nv,nv+1+QA) !------------------------------------------------------------ !$omp parallel do collapse(2) & @@ -351,40 +407,44 @@ subroutine linkernel_solve_sip( b, G1, G2, & ! (inout) 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 + 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 + nrhs2 = 1 + QA + do iv=1, 1+QA + do pv1=1, nv + pp = i + (pv1-1)*im + RHS2(pv1,iv) = b(pp,2+iv,j,ke_xy) + end do + 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 + nrhs1 = nv + 2 + nrhs2 = nv + 1 + QA + + 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 iv=1, 1+QA + do pv1=1, nv + pp = i + (pv1-1)*im + RHS2(pv1,nv+iv) = b(pp,2+iv,j,ke_xy) + end do + end do end if @@ -413,50 +473,55 @@ subroutine linkernel_solve_sip( b, G1, G2, & ! (inout) 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 + ! 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 iv=1, 1+QA + do pv1=1, nv + pp = i + (pv1-1)*im + b(pp,2+iv,j,ke_xy) = RHS2(pv1,iv) + end do + 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 + ! 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 + ! 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 + ! RHS1(:,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 + ! RHS2(:,nv+1:nv+1+QA) contains D^{-1} b. + do iv=1, 1+QA + do pv1=1, nv + pp = i + (pv1-1)*im + b(pp,2+iv,j,ke_xy) = RHS2(pv1,nv+iv) + end do + end do end if - end do ! i + end do ! i end do ! j end do ! ke_xy !$omp end parallel do @@ -464,6 +529,7 @@ subroutine linkernel_solve_sip( b, G1, G2, & ! (inout) return end subroutine linkernel_solve_sip + !> Construct the block tridiagonal matrix for the vertical diffusion equation !OCL SERIAL subroutine construct_matbnd_sip( BndMatL, BndMatD, BndMatU, & ! (out) RHO, KDIFF, GsqrtV, Dx1D, M1D,invM1D, & ! (in) @@ -721,8 +787,8 @@ subroutine construct_sip_face_blocks_lgl( & 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, & + subroutine eval_Ax( MOMX_t, MOMY_t, DRHOT_t, RHOQ_t_list, alph_M, alph_H, & + PROG_VARS, MOMX00, MOMY00, PT00, RHOQ00, 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 @@ -731,12 +797,14 @@ subroutine eval_Ax( MOMX_t, MOMY_t, DRHOT_t, alph_M, alph_H, & 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) + type(LocalMeshFieldBaseList), intent(inout) :: RHOQ_t_list(QA) 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) :: PROG_VARS(elem%Np,lmesh%NeX*lmesh%NeY*lmesh%NeZ,3+QA) 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) :: RHOQ00(elem%Np,QA,lmesh%Ne) 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) @@ -749,23 +817,25 @@ subroutine eval_Ax( MOMX_t, MOMY_t, DRHOT_t, alph_M, alph_H, & 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) + real(RP), intent(out), optional :: b(im,elem%Nnode_v,3+QA,jm,lmesh%Ne) logical, intent(in), optional :: use_delta_form real(RP) :: Flux(elem%Np,3), DFlux(elem%Np,2,3) + real(RP) :: Flux_q(elem%Np), DFlux_q(elem%Np,2) 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) :: DIFF_flux_z_broken(elem%Np,lmesh%NeA,3+QA) + real(RP) :: DIFF_flux_z(elem%Np,lmesh%NeA,3+QA) real(RP) :: del_flux(elem%NfpTot,3,lmesh%Ne) + real(RP) :: del_flux_q(elem%NfpTot,QA,lmesh%Ne) integer :: ke_xy, ke_z integer :: ke, ke2D integer :: p, fp - integer :: iv + integer :: iv, iq integer :: i, j integer :: pv @@ -782,12 +852,12 @@ subroutine eval_Ax( MOMX_t, MOMY_t, DRHOT_t, alph_M, alph_H, & end if if ( flag_use_delta_form ) then - call cal_grad_del_flux( del_flux, & ! (out) + call cal_grad_del_flux( del_flux, del_flux_q, & ! (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 ) + !$omp parallel do collapse(2) private( ke,ke2D,iv, DFlux,DFlux_q, 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 @@ -810,17 +880,30 @@ subroutine eval_Ax( MOMX_t, MOMY_t, DRHOT_t, alph_M, alph_H, & 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 + + do iq=1, QA + iv = 3 + iq + call element3D_operation%Dz( PROG_VARS(:,ke,iv) * RDENS(:), DFlux_q(:,1) ) + call element3D_operation%Lift( del_flux_q(:,iq,ke), DFlux_q(:,2) ) + 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,iv) = DENS(p,ke) * KH(p,ke) * E33 * DFlux_q(p,1) * RGsqrtV + DIFF_flux_z(p,ke,iv) = DENS(p,ke) * KH(p,ke) * ( E33 * DFlux_q(p,1) + DFlux_q(p,2) ) * RGsqrtV + end do + end do + end do end do !--------------------------------------------------------- - call cal_del_flux( del_flux, alph_M, alph_H, & ! (out) + call cal_del_flux( del_flux, del_flux_q, 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 ) + !$omp parallel do collapse(2) private( ke,iv, ke2D, DFlux,DFlux_q, 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 @@ -837,43 +920,81 @@ subroutine eval_Ax( MOMX_t, MOMY_t, DRHOT_t, alph_M, alph_H, & 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 + + do iq=1, QA + iv = 3 + iq + call element3D_operation%Dz( DIFF_flux_z(:,ke,iv), DFlux_q(:,1) ) + call element3D_operation%Lift( del_flux_q(:,iq,ke), DFlux_q(:,2) ) + do p=1, elem%Np + RGsqrtV = 1.0_RP / GsqrtV(p,ke) + E33 = lmesh%Escale(p,ke,3,3) + RHOQ_t_list(iq)%ptr%val(p,ke) = ( E33 * DFlux_q(p,1) + DFlux_q(p,2) ) * RGsqrtV + end do + end do + end do end do end if if ( flag_cal_b ) then if ( use_delta_form ) then - !$omp parallel do private(p) + !$omp parallel do private(p,iv) 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 + + 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 + !- + do iq=1, QA + iv = 3 + iq + do j=1, jm + do i=1, im + p = i + (j-1)*im + (pv-1)*im*jm + b(i,pv,iv,j,ke) = impl_fac * RHOQ_t_list(iq)%ptr%val(p,ke) & + - PROG_VARS (p,ke,iv) & + + RHOQ00(p,iq,ke) + end do + end do + end do + end do end do else - !$omp parallel do private(p) + !$omp parallel do private(p,iv) 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 + + 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 + !- + do iq=1, QA + iv = 3 + iq + do j=1, jm + do i=1, im + p = i + (j-1)*im + (pv-1)*im*jm + b(i,pv,iv,j,ke) = RHOQ00(p,iq,ke) + end do + end do + end do + end do end do end if @@ -882,7 +1003,7 @@ subroutine eval_Ax( MOMX_t, MOMY_t, DRHOT_t, alph_M, alph_H, & end subroutine eval_Ax !OCL SERIAL - subroutine cal_del_flux( del_flux, alph_M, alph_H, & + subroutine cal_del_flux( del_flux, del_flux_q, 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 @@ -891,11 +1012,12 @@ subroutine cal_del_flux( del_flux, alph_M, alph_H, & 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) :: del_flux_q(elem%NfpTot,QA,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) :: DIFF_flux_z(elem%Np*lmesh%NeA,3+QA) + real(RP), intent(in) :: DIFF_flux_z_broken(elem%Np*lmesh%NeA,3+QA) + real(RP), intent(in) :: PROG_VARS(elem%Np*lmesh%NeX*lmesh%NeY*lmesh%NeZ,3+QA) 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) @@ -907,16 +1029,19 @@ subroutine cal_del_flux( del_flux, alph_M, alph_H, & logical, intent(in) :: is_bound(elem%NfpTot,lmesh%Ne) integer :: ke, ke_z, ke_xy + integer :: iq + 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) + real(RP) :: numflux_q(elem%NfpTot) !------------------------------------------------ !$omp parallel do collapse(2) & - !$omp private( ke, iM, iP, alph_M, alph_H, coef, RDENS_M, RDENS_P, numflux ) + !$omp private( ke,iq, iM,iP, coef, RDENS_M, RDENS_P, numflux,numflux_q ) do ke_z=1, lmesh%NeZ do ke_xy=1, lmesh%NeX*lmesh%NeY ke = ke_xy + (ke_z-1)*lmesh%Ne2D @@ -931,6 +1056,8 @@ subroutine cal_del_flux( del_flux, alph_M, alph_H, & where ( is_bound(:,ke) ) numflux(:,1) = 0.0_RP + numflux(:,2) = 0.0_RP + numflux(:,3) = 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) ) @@ -946,6 +1073,17 @@ subroutine cal_del_flux( del_flux, alph_M, alph_H, & 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(:) ) + + do iq=1, QA + where ( is_bound(:,ke) ) + numflux_q(:) = 0.0_RP + elsewhere + numflux_q(:) = 0.5_RP * ( DIFF_flux_z_broken(iP,3+iq) + DIFF_flux_z(iM,3+iq) ) + end where + + del_flux_q(:,iq,ke) = coef(:) * ( numflux_q(:) - DIFF_flux_z(iM,3+iq) ) & + + alph_H(:,ke) * Fscale(:,ke) * ( PROG_VARS(iP,3+iq) * RDENS_P(:) - PROG_VARS(iM,3+iq) * RDENS_M(:) ) + end do end do end do @@ -953,7 +1091,7 @@ subroutine cal_del_flux( del_flux, alph_M, alph_H, & end subroutine cal_del_flux !OCL SERIAL - subroutine cal_grad_del_flux( del_flux, & + subroutine cal_grad_del_flux( del_flux, del_flux_q, & PROG_VARS, DENS, nz, Fscale,vmapM, vmapP, & lmesh, elem, lmesh2D, elem2D ) implicit none @@ -962,7 +1100,8 @@ subroutine cal_grad_del_flux( del_flux, & 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(out) :: del_flux_q(elem%NfpTot,QA,lmesh%Ne) + real(RP), intent(in) :: PROG_VARS(elem%Np*lmesh%NeX*lmesh%NeY*lmesh%NeZ,3+QA) 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) @@ -970,6 +1109,7 @@ subroutine cal_grad_del_flux( del_flux, & integer, intent(in) :: vmapP(elem%NfpTot,lmesh%Ne) integer :: ke, ke_z, ke_xy + integer :: iq integer :: iP(elem%NfpTot), iM(elem%NfpTot) real(RP) :: coef(elem%NfpTot) @@ -977,7 +1117,7 @@ subroutine cal_grad_del_flux( del_flux, & !------------------------------------------------ !$omp parallel do collapse(2) & - !$omp private( ke, iM, iP, coef, RDENS_M, RDENS_P ) + !$omp private( ke,iq, 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 @@ -990,6 +1130,10 @@ subroutine cal_grad_del_flux( del_flux, & 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(:) ) + + do iq=1, QA + del_flux_q(:,iq,ke) = coef(:) * ( PROG_VARS(iP,3+iq) * RDENS_P(:) - PROG_VARS(iM,3+iq) * RDENS_M(:) ) + end do end do end do 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 index b3b7ecba..e4c418ad 100644 --- 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 @@ -153,7 +153,7 @@ end subroutine atm_phy_bl_dgm_mynn_lv2_Final 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) + Rtot, PRES, PT, & ! (in) Dz, Lift, lmesh, elem, is_bound ) ! (in) implicit none class(LocalMesh3D), intent(in) :: lmesh @@ -168,6 +168,7 @@ subroutine atm_phy_bl_dgm_mynn_lv2_cal_VViscDiffCoef( & 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) :: Rtot(elem%Np,lmesh%NeA) !< Gas constant 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 @@ -175,13 +176,15 @@ subroutine atm_phy_bl_dgm_mynn_lv2_cal_VViscDiffCoef( & 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) :: DdensDz(elem%Np), DVelDz(elem%Np,3), DptDz(elem%Np), DrtotDz(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) :: del_flux_Rtot(elem%NfpTot,lmesh%Ne) real(RP) :: S2(elem%Np) ! Square of vertical shear of horizontal wind + real(RP) :: N2 ! Square of Brunt-Vaisala frequency real(RP) :: Ri ! Richardson number real(RP) :: Rf(elem%Np) ! Flux Richardson number @@ -225,14 +228,14 @@ subroutine atm_phy_bl_dgm_mynn_lv2_cal_VViscDiffCoef( & !--------------------------------- - 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) + call cal_del_flux_grad( del_flux_rho, del_flux_mom, del_flux_rhot, del_flux_Rtot, & ! (out) + DDENS_, MOMX_, MOMY_, MOMZ_, DRHOT_, Rtot, DENS_hyd, PRES_hyd, PT, & ! (in) + lmesh%normal_fn(:,:,3), lmesh%Fscale, 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 private( Fz, LiftDelFlx, DENS, RDENS, RHOT, Q, DdensDz, DVelDz, DptDz, DrtotDz, & + !$omp N2, S2, Ri, Rf, discriminant, denom_m, denom_h, S_M, S_H, mixlen, kz ) !$omp do do ke2D=lmesh2D%NeS, lmesh2D%NeE @@ -247,31 +250,38 @@ subroutine atm_phy_bl_dgm_mynn_lv2_cal_VViscDiffCoef( & ! gradient of density call sparsemat_matmul( Dz, DENS, Fz ) - call sparsemat_matmul( Lift, lmesh%Fscale(:,ke) * del_flux_rho(:,ke), LiftDelFlx ) + call sparsemat_matmul( Lift, 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 ) + call sparsemat_matmul( Lift, 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 ) + call sparsemat_matmul( Lift, 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 ) + call sparsemat_matmul( Dz, RHOT(:) , Fz ) + call sparsemat_matmul( Lift, del_flux_rhot(:,ke), LiftDelFlx ) DptDz(:) = ( lmesh%Escale(:,ke,3,3) * Fz(:) + LiftDelFlx(:) - Q(:) * DdensDz(:) ) * RDENS(:) + ! gradient of Rtot + call sparsemat_matmul( Dz, Rtot(:,ke), Fz ) + call sparsemat_matmul( Lift, del_flux_Rtot(:,ke), LiftDelFlx ) + DrtotDz(:) = lmesh%Escale(:,ke,3,3) * Fz(:) + LiftDelFlx(:) + ! Calculate flux Richardson number: Rf do p=1, elem%Np + N2 = GRAV * ( DptDz(p) / PT(p,ke) + DrtotDz(p) / Rtot(p,ke) ) + S2(p) = DVelDz(p,1)**2 + DVelDz(p,2)**2 - Ri = GRAV * DptDz(p) / ( PT(p,ke) * max(S2(p), EPS) ) + Ri = N2 / max(S2(p), EPS) discriminant = Ri * Ri & + 2.0_RP * AF12 * ( Rf1 - 2.0_RP * Rf2 ) * Ri & @@ -330,9 +340,9 @@ 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) + subroutine cal_del_flux_grad( del_flux_rho, del_flux_mom, del_flux_rhot, del_flux_Rtot, & ! (out) + DDENS_, MOMX_, MOMY_, MOMZ_, DRHOT_, Rtot, DENS_hyd, PRES_hyd, PT_, & ! (in) + nz, Fscale, vmapM, vmapP, lmesh, elem, is_bound ) ! (in) implicit none @@ -341,26 +351,28 @@ subroutine cal_del_flux_grad( del_flux_rho, del_flux_mom, del_flux_rhot, & ! (o 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(out) :: del_flux_Rtot(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) :: DRHOT_(elem%Np*lmesh%NeA) + real(RP), intent(in) :: Rtot(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) + 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 :: i, iP, iM real(RP) :: densM, densP - real(RP) :: del real(RP) :: facz !------------------------------------------------------------------------ - !$omp parallel do private ( iM, iP, densM, densP, del, facz ) + !$omp parallel do private ( iM, iP, densM, densP, facz ) do i=1, elem%NfpTot * lmesh%Ne iM = vmapM(i); iP = vmapP(i) @@ -368,26 +380,20 @@ subroutine cal_del_flux_grad( del_flux_rho, del_flux_mom, del_flux_rhot, & ! (o densP = DDENS_(iP) + DENS_hyd(iP) if ( is_bound(i) ) then - facz = 1.0_RP + facz = 0.5_RP * Fscale(i) * nz(i) else - ! facz = 1.0_RP - sign(1.0_RP,nz(i)) - facz = 1.0_RP + ! facz = ( 1.0_RP - sign(1.0_RP,nz(i)) ) * Fscale(i) * nz(i) + facz = 0.5_RP * Fscale(i) * nz(i) 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_flux_rho(i) = facz * ( densP - densM ) - del = 0.5_RP * ( MOMZ_(iP) - MOMZ_(iM) ) - del_flux_mom(i,3) = facz * del * nz(i) + del_flux_mom(i,1) = facz * ( MOMX_(iP) - MOMX_(iM) ) + del_flux_mom(i,2) = facz * ( MOMY_(iP) - MOMY_(iM) ) + del_flux_mom(i,3) = facz * ( MOMZ_(iP) - MOMZ_(iM) ) - del = 0.5_RP * ( densP * PT_(iP) - densM * PT_(iM) ) - del_flux_rhot(i) = facz * del * nz(i) + del_flux_rhot(i) = facz * ( densP * PT_(iP) - densM * PT_(iM) ) + del_flux_Rtot(i) = facz * ( Rtot(iP) - Rtot(iM) ) end do return 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 c711acb6..43fbc5cb 100644 --- a/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_bl.F90 +++ b/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_bl.F90 @@ -212,7 +212,8 @@ 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_tracer, only: & + TRACER_ADVC, 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: & @@ -252,10 +253,13 @@ subroutine AtmosPhyBl_calc_tendency( & class(LocalMeshFieldBase), pointer :: DDENS, MOMX, MOMY, MOMZ, DRHOT class(LocalMeshFieldBase), pointer :: DENS_hyd, PRES_hyd, Rtot, CVtot, CPtot class(LocalMeshFieldBase), pointer :: PRES, PT + type(LocalMeshFieldBaseList) :: RHOQ_list(QA) 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 + class(LocalMeshFieldBase), pointer :: bl_RHOQ_t + type(LocalMeshFieldBaseList) :: bl_RHOQ_t_list(QA) type DYN_BNDInfo logical, allocatable :: is_bound(:,:) @@ -290,10 +294,14 @@ subroutine AtmosPhyBl_calc_tendency( & mesh, auxvars_list, & PRES, PT ) - call AtmosPhyBLVars_GetLocalMeshFields_tend( n, & - mesh, this%vars%tends_manager, & - bl_RHOU_t, bl_RHOV_t, bl_RHOT_t ) - + call AtmosVars_GetLocalMeshQTRCVarList( n, & + mesh, trcvars_list, & + 1, RHOQ_list ) + + call AtmosPhyBLVars_GetLocalMeshFields_tend( n, & + mesh, this%vars%tends_manager, & + bl_RHOU_t, bl_RHOV_t, bl_RHOT_t, bl_RHOQ_t_list ) + !- allocate( bnd_info(n)%is_bound(lcmesh%refElem3D%NfpTot,lcmesh%Ne) ) call this%dyn_bnd%Inquire_bound_flag( bnd_info(n)%is_bound, & ! (out) @@ -303,18 +311,20 @@ subroutine AtmosPhyBl_calc_tendency( & 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) + 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, Rtot%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) + bl_RHOQ_t_list, & ! (out) DDENS%val, MOMX%val, MOMY%val, DRHOT%val, & ! (in) + RHOQ_list, & ! (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) @@ -339,9 +349,9 @@ subroutine AtmosPhyBl_calc_tendency( & RHOH_p, RHOQ_tp ) call AtmosPhyBLVars_GetLocalMeshFields_tend( n, & - mesh, this%vars%tends_manager, & - bl_RHOU_t, bl_RHOV_t, bl_RHOT_t, & - lcmesh ) + mesh, this%vars%tends_manager, & + bl_RHOU_t, bl_RHOV_t, bl_RHOT_t, bl_RHOQ_t_list, & + lcmesh ) !$omp parallel private(ke, iq) !$omp do @@ -350,11 +360,20 @@ subroutine AtmosPhyBl_calc_tendency( & 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 do + do iq=1, QA + if ( .not. TRACER_ADVC(iq) ) cycle + !$omp do + do ke=lcmesh%NeS, lcmesh%NeE + RHOQ_tp(iq)%ptr%val(:,ke) = RHOQ_tp(iq)%ptr%val(:,ke) & + + bl_RHOQ_t_list(iq)%ptr%val(:,ke) + end do + !$omp end do + end do !$omp end parallel end do call PROF_rapend('ATM_PHY_BL_add_tend', 2) - return end subroutine AtmosPhyBl_calc_tendency 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 b6c613a5..4fc5a8cc 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 @@ -17,6 +17,7 @@ module mod_atmos_phy_bl_vars use scale_precision use scale_io use scale_prc + use scale_tracer, only: QA use scale_element_base, only: ElementBase3D use scale_mesh_base, only: MeshBase @@ -127,8 +128,6 @@ subroutine AtmosPhyBlVars_Init( this, model_mesh, & this%QS = QS_BL this%QE = QE_BL this%QA = QA_BL - this%TENDS_NUM_TOT = ATMOS_PHY_BL_TENDS_NUM1 + QE_BL - QS_BL + 1 - return end subroutine AtmosPhyBlVars_Init @@ -157,6 +156,8 @@ subroutine AtmosPhyBlVars_Setup( this, model_mesh ) type(VariableInfo) :: qtrc_vterm_vinfo_tmp !---------------------------------------------------- + this%TENDS_NUM_TOT = ATMOS_PHY_BL_TENDS_NUM1 + QA + !- Initialize auxiliary and diagnostic variables nullify( atm_mesh ) @@ -184,7 +185,7 @@ subroutine AtmosPhyBlVars_Setup( this, model_mesh ) qtrc_tp_vinfo_tmp%dim_type = 'XYZ' qtrc_tp_vinfo_tmp%STDNAME = '' - do iq = 1, this%QA + do iq = 1, QA iv = ATMOS_PHY_BL_TENDS_NUM1 + iq qtrc_tp_vinfo_tmp%keyID = iv qtrc_tp_vinfo_tmp%NAME = 'BL_'//trim(TRACER_NAME(this%QS+iq-1))//'_t' @@ -231,7 +232,7 @@ end subroutine AtmosPhyBlVars_Final !OCL SERIAL subroutine AtmosPhyBlVars_GetLocalMeshFields_tend( domID, mesh, bl_tends_list, & - bl_RHOU_t, bl_RHOV_t, bl_RHOT_t, & + bl_RHOU_t, bl_RHOV_t, bl_RHOT_t, bl_RHOQ_t, & lcmesh3D & ) @@ -245,6 +246,7 @@ subroutine AtmosPhyBlVars_GetLocalMeshFields_tend( domID, mesh, 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 + type(LocalMeshFieldBaseList), intent(out), optional :: bl_RHOQ_t(:) class(LocalMesh3D), pointer, intent(out), optional :: lcmesh3D class(MeshFieldBase), pointer :: field @@ -263,6 +265,14 @@ subroutine AtmosPhyBlVars_GetLocalMeshFields_tend( domID, mesh, bl_tends_list, & call bl_tends_list%Get(ATMOS_PHY_BL_RHOT_t_ID, field) call field%GetLocalMeshField(domID, bl_RHOT_t) + !--- + if ( present(bl_RHOQ_t) ) then + do iq = 1, size(bl_RHOQ_t) + call bl_tends_list%Get(ATMOS_PHY_BL_TENDS_NUM1 + iq, field) + call field%GetLocalMeshField(domID, bl_RHOQ_t(iq)%ptr) + end do + end if + if (present(lcmesh3D)) then call mesh%GetLocalMesh( domID, lcmesh ) nullify( lcmesh3D ) From 122e2cbd92758e3e60f8359b7d3db97b13b6389b Mon Sep 17 00:00:00 2001 From: Yuta Kawai Date: Mon, 27 Jul 2026 11:38:50 +0900 Subject: [PATCH 7/9] Update a boundary layer test to include the moisture --- .../test/case/boundary_layer/init.conf | 6 +-- .../test/case/boundary_layer/mod_user.F90 | 54 ++++++++++++++----- .../test/case/boundary_layer/run.conf | 5 +- .../boundary_layer/visualize/visualize.sh | 8 +-- 4 files changed, 52 insertions(+), 21 deletions(-) diff --git a/model/atm_nonhydro3d/test/case/boundary_layer/init.conf b/model/atm_nonhydro3d/test/case/boundary_layer/init.conf index ca90c4b3..b442c302 100644 --- a/model/atm_nonhydro3d/test/case/boundary_layer/init.conf +++ b/model/atm_nonhydro3d/test/case/boundary_layer/init.conf @@ -16,9 +16,9 @@ &PARAM_CONST / &PARAM_EXP - U0 = 10.0D0, - TEMP0 = 300.0D0, - ENV_RH = 0D0, + U0 = 10.0D0, + SFC_POTT = 300.0D0, + ENV_RH = 1D0, NITER_RH = 3, / #** ATMOS ****************************************************** diff --git a/model/atm_nonhydro3d/test/case/boundary_layer/mod_user.F90 b/model/atm_nonhydro3d/test/case/boundary_layer/mod_user.F90 index ae8b864f..1e310a4a 100644 --- a/model/atm_nonhydro3d/test/case/boundary_layer/mod_user.F90 +++ b/model/atm_nonhydro3d/test/case/boundary_layer/mod_user.F90 @@ -186,7 +186,7 @@ subroutine exp_SetInitCond_boundary_layer( this, & 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_POTT = 250.0_RP !< Initial surface potential 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] @@ -194,8 +194,8 @@ subroutine exp_SetInitCond_boundary_layer( this, & real(RP) :: ENV_RH = 0.0_RP !< Relative Humidity of environment [%] integer :: NITER_RH = 3 namelist /PARAM_EXP/ & - TEMP0, & PTLAPS, & + SFC_POTT, & U0, & VSHEAR_HVEL, & ENV_RH, & @@ -213,10 +213,19 @@ subroutine exp_SetInitCond_boundary_layer( this, & 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) :: sfc_rhot(elem%Nfp_v) real(RP) :: bnd_SFC_PRES(elem%Nnode_h1D**2,lcmesh%lcmesh2D%NeA) + + real(RP) :: temp_z(elem%Nnode_v,lcmesh%NeZ) + real(RP) :: pres_z(elem%Nnode_v,lcmesh%NeZ) + real(RP) :: dens_z(elem%Nnode_v,lcmesh%NeZ) + real(RP) :: qv_z (elem%Nnode_v,lcmesh%NeZ) + real(RP) :: rtot_z(elem%Nnode_v,lcmesh%NeZ) + real(RP) :: cptot_z(elem%Nnode_v,lcmesh%NeZ) + real(RP) :: cvtot_z(elem%Nnode_v,lcmesh%NeZ) + real(RP) :: psat_z(elem%Nnode_v,lcmesh%NeZ) real(RP) :: QV(elem%Np) - real(RP) :: PRES(elem%Np) - real(RP) :: psat0 integer :: itr !----------------------------------------------------------------------------- @@ -234,8 +243,8 @@ subroutine exp_SetInitCond_boundary_layer( this, & !--- call hydrostatic_calc_basicstate_constPTLAPS( & - DENS_hyd, PRES_hyd, & ! (out) - PTLAPS, TEMP0, SFC_PRES, x, y, z, lcmesh, elem ) ! (in) + DENS_hyd, PRES_hyd, & ! (out) + PTLAPS, SFC_POTT, 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 @@ -243,7 +252,7 @@ subroutine exp_SetInitCond_boundary_layer( this, & 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) + PT_tmp(:,ke_z,ke_x,ke_y) = SFC_POTT + PTLAPS * z(:,ke) RHOT_hyd(:,ke) = DENS_hyd(:,ke) * PT_tmp(:,ke_z,ke_x,ke_y) end do end do @@ -256,12 +265,31 @@ subroutine exp_SetInitCond_boundary_layer( this, & 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_z=1, lcmesh%NeZ + do p3=1, elem%Nnode_v + ke = 1 + (ke_z - 1)*lcmesh%Ne2D + p = 1 + (p3 - 1)*elem%Nnode_h1D**2 + + QV(p) = tracer_field_list(iq_QV)%ptr%val(p,ke) + rtot_z (p3,ke_z) = Rdry * ( 1.0_RP - QV(p) ) + Rvap * QV(p) + cptot_z(p3,ke_z) = CPdry * ( 1.0_RP - QV(p) ) + CP_VAPOR * QV(p) + cvtot_z(p3,ke_z) = CVdry * ( 1.0_RP - QV(p) ) + CV_VAPOR * QV(p) + + dens_z(p3,ke_z) = DENS_hyd(p,ke) + DDENS(p,ke) + pres_z(p3,ke_z) = PRES00 * ( rtot_z(p3,ke_z) / PRES00 * dens_z(p3,ke_z) * PT_tmp(p,ke_z,1,1) )**(cptot_z(p3,ke_z)/cvtot_z(p3,ke_z)) + temp_z(p3,ke_z) = pres_z(p3,ke_z) / ( rtot_z(p3,ke_z) * dens_z(p3,ke_z) ) + end do + end do + call ATMOS_SATURATION_psat_all( & + elem%Nnode_v, 1, elem%Nnode_v, lcmesh%NeZ, 1, lcmesh%NeZ, & + temp_z, & + psat_z ) ! [OUT] + + !$omp parallel do collapse(3) private(ke_z,ke_x,ke_y,ke,ke2D,p3,p2D,p, QV,sfc_rhot) do ke_y=1, lcmesh%NeY do ke_x=1, lcmesh%NeX do ke_z=1, lcmesh%NeZ @@ -271,7 +299,7 @@ subroutine exp_SetInitCond_boundary_layer( this, & 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 ) + QV(p) = ENV_RH * 1.0E-2_RP * psat_z(p3,ke_z) / ( ( DENS_hyd(p,ke) + DDENS(p,ke) ) * Rvap * temp_z(p3,ke_z) ) end do end do tracer_field_list(iq_QV)%ptr%val(:,ke) = QV(:) @@ -281,10 +309,9 @@ subroutine exp_SetInitCond_boundary_layer( this, & 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)) + sfc_rhot(:) = ( DDENS(elem%Hslice(:,1),ke2D) + DENS_hyd(elem%Hslice(:,1),ke2D) ) * PT_tmp(elem%Hslice(:,1),ke_z,ke_x,ke_y) + bnd_SFC_PRES(:,ke2D) = PRES00 * ( Rtot(elem%Hslice(:,1),ke_z,ke_x,ke_y) * sfc_rhot(:) / PRES00 )**CPtot_ov_CVtot(elem%Hslice(:,1),ke_z,ke_x,ke_y) end if end do end do @@ -304,6 +331,7 @@ subroutine exp_SetInitCond_boundary_layer( this, & 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) diff --git a/model/atm_nonhydro3d/test/case/boundary_layer/run.conf b/model/atm_nonhydro3d/test/case/boundary_layer/run.conf index 0ee4fcfb..e8d9e493 100644 --- a/model/atm_nonhydro3d/test/case/boundary_layer/run.conf +++ b/model/atm_nonhydro3d/test/case/boundary_layer/run.conf @@ -9,7 +9,7 @@ TIME_STARTMS = 0.D0, TIME_DURATION = 7200D0, TIME_DURATION_UNIT = 'SEC', - TIME_DT = 180D0, + TIME_DT = 180D0, TIME_DT_UNIT = 'SEC', / &PARAM_CONST @@ -115,9 +115,10 @@ &HISTORY_ITEM name='U' / &HISTORY_ITEM name='BL_RHOU_t' / &HISTORY_ITEM name='BL_RHOT_t' / +&HISTORY_ITEM name='BL_QV_t' / &HISTORY_ITEM name='NU' / &HISTORY_ITEM name='KH' / -!&HISTORY_ITEM name='QV' / +&HISTORY_ITEM name='QV' / #*** Statistics ******************************************* diff --git a/model/atm_nonhydro3d/test/case/boundary_layer/visualize/visualize.sh b/model/atm_nonhydro3d/test/case/boundary_layer/visualize/visualize.sh index 920a3099..142dc53b 100644 --- a/model/atm_nonhydro3d/test/case/boundary_layer/visualize/visualize.sh +++ b/model/atm_nonhydro3d/test/case/boundary_layer/visualize/visualize.sh @@ -12,11 +12,13 @@ mkdir -p analysis ### 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@U,x=0e0,y=0e0 analysis/U.png --prc_num_xy 1 1 --figsize 8 4 +python ../common/cmd_mkgraph.py history.pe00\*.nc@QV,x=0e0,y=0e0 analysis/QV.png --prc_num_xy 1 1 --figsize 8 4 --range 1e-4 3e-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 +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@BL_QV_t,x=0e0,y=0e0 analysis/BL_QV_t.png --prc_num_xy 1 1 --figsize 8 4 --range 0 5e-8 ### make animation ### # echo "+make animation" From cbc64978864c61e45906a395b78a58a08194c1e2 Mon Sep 17 00:00:00 2001 From: Yuta Kawai Date: Mon, 27 Jul 2026 11:51:20 +0900 Subject: [PATCH 8/9] Fix stage initialization of RK scheme --- FElib/src/common/scale_timeint_rk.F90 | 6 +++--- FElib/src/common/scale_timeint_rk.F90.erb | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/FElib/src/common/scale_timeint_rk.F90 b/FElib/src/common/scale_timeint_rk.F90 index f083c011..708fb47a 100644 --- a/FElib/src/common/scale_timeint_rk.F90 +++ b/FElib/src/common/scale_timeint_rk.F90 @@ -1895,7 +1895,7 @@ subroutine rk_advance_general1D( this, nowstage, q, varID, is, ie , IA, var_num, tintbuf_ind = this%tend_buf_indmap(nowstage) - if ( this%nstage == 1 .and. (.not. this%imex_flag) ) then + if ( nowstage == 1 .and. (.not. this%imex_flag) ) then !$omp parallel do !$acc parallel loop collapse(1) present( q, varTmp_1d ) do i=is, ie @@ -2219,7 +2219,7 @@ subroutine rk_advance_general2D( this, nowstage, q, varID, is, ie ,js, je , IA,J tintbuf_ind = this%tend_buf_indmap(nowstage) - if ( this%nstage == 1 .and. (.not. this%imex_flag) ) then + if ( nowstage == 1 .and. (.not. this%imex_flag) ) then !$omp parallel do !$acc parallel loop collapse(2) present( q, varTmp_2d ) do j=js, je @@ -2583,7 +2583,7 @@ subroutine rk_advance_general3D( this, nowstage, q, varID, is, ie ,js, je ,ks, k tintbuf_ind = this%tend_buf_indmap(nowstage) - if ( this%nstage == 1 .and. (.not. this%imex_flag) ) then + if ( nowstage == 1 .and. (.not. this%imex_flag) ) then !$omp parallel do collapse(2) !$acc parallel loop collapse(3) present( q, varTmp_3d ) do k=ks, ke diff --git a/FElib/src/common/scale_timeint_rk.F90.erb b/FElib/src/common/scale_timeint_rk.F90.erb index af769a26..2d2d5ca5 100644 --- a/FElib/src/common/scale_timeint_rk.F90.erb +++ b/FElib/src/common/scale_timeint_rk.F90.erb @@ -998,7 +998,7 @@ contains tintbuf_ind = this%tend_buf_indmap(nowstage) - if ( this%nstage == 1 .and. (.not. this%imex_flag) ) then + if ( nowstage == 1 .and. (.not. this%imex_flag) ) then % if (d > 2) then !$omp parallel do collapse(<%=(d-1)%>) % else From d0ec172acbe151f1e70c7ea976def15f92511d27 Mon Sep 17 00:00:00 2001 From: Yuta Kawai Date: Mon, 27 Jul 2026 12:40:33 +0900 Subject: [PATCH 9/9] Fix OpenACC parallel loop --- FElib/src/common/scale_timeint_rk.F90 | 6 +++--- FElib/src/common/scale_timeint_rk.F90.erb | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/FElib/src/common/scale_timeint_rk.F90 b/FElib/src/common/scale_timeint_rk.F90 index 708fb47a..045a1436 100644 --- a/FElib/src/common/scale_timeint_rk.F90 +++ b/FElib/src/common/scale_timeint_rk.F90 @@ -1897,7 +1897,7 @@ subroutine rk_advance_general1D( this, nowstage, q, varID, is, ie , IA, var_num, if ( nowstage == 1 .and. (.not. this%imex_flag) ) then !$omp parallel do - !$acc parallel loop collapse(1) present( q, varTmp_1d ) + !$acc parallel loop collapse(1) present( q, var0_1d, varTmp_1d ) do i=is, ie var0_1d(i,varID) = q(i) varTmp_1d(i,varID) = q(i) @@ -2221,7 +2221,7 @@ subroutine rk_advance_general2D( this, nowstage, q, varID, is, ie ,js, je , IA,J if ( nowstage == 1 .and. (.not. this%imex_flag) ) then !$omp parallel do - !$acc parallel loop collapse(2) present( q, varTmp_2d ) + !$acc parallel loop collapse(2) present( q, var0_2d, varTmp_2d ) do j=js, je do i=is, ie var0_2d(i,j,varID) = q(i,j) @@ -2585,7 +2585,7 @@ subroutine rk_advance_general3D( this, nowstage, q, varID, is, ie ,js, je ,ks, k if ( nowstage == 1 .and. (.not. this%imex_flag) ) then !$omp parallel do collapse(2) - !$acc parallel loop collapse(3) present( q, varTmp_3d ) + !$acc parallel loop collapse(3) present( q, var0_3d, varTmp_3d ) do k=ks, ke do j=js, je do i=is, ie diff --git a/FElib/src/common/scale_timeint_rk.F90.erb b/FElib/src/common/scale_timeint_rk.F90.erb index 2d2d5ca5..c27cd87d 100644 --- a/FElib/src/common/scale_timeint_rk.F90.erb +++ b/FElib/src/common/scale_timeint_rk.F90.erb @@ -1004,7 +1004,7 @@ contains % else !$omp parallel do % end - !$acc parallel loop collapse(<%=(d)%>) present( q, varTmp_<%=d%>d ) + !$acc parallel loop collapse(<%=(d)%>) present( q, var0_<%=d%>d, varTmp_<%=d%>d ) % for i in 1..d do <%=ind_name[d-i]%>=<%=ind_range_list[d-i]%> % end