From 038621b89a53dc0db491fd722e40caf42034c095 Mon Sep 17 00:00:00 2001 From: ykawai Date: Wed, 5 Aug 2026 15:10:09 +0900 Subject: [PATCH 1/6] Fix bugs --- .../scale_atm_phy_bl_dgm_common.F90 | 30 +++++++++---------- .../scale_atm_phy_bl_dgm_mynn_lv2.F90 | 4 +-- .../radiation/scale_atm_phy_rd_dgm_simple.F90 | 3 +- .../src/atmos/mod_atmos_phy_bl.F90 | 10 +++---- .../src/atmos/mod_atmos_phy_bl_vars.F90 | 4 +-- .../src/atmos/mod_atmos_phy_sfc.F90 | 14 ++++----- .../src/ocean/mod_ocean_mesh_gm.F90 | 4 ++- 7 files changed, 36 insertions(+), 33 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 6729cbef..4abd36a3 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 @@ -55,11 +55,11 @@ module scale_atm_phy_bl_dgm_common !OCL SERIAL subroutine atm_phy_bl_dgm_common_calc_tendency( & 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) - use_delta_form ) ! (in) + DDENS_, MOMX_, MOMY_, DRHOT_, QTRC_list, & ! (in) + PT_, DENS_hyd, PRES_hyd, NU, KH, & ! (in) + element3D_operation, C_IP, dtsec, & ! (in) + lmesh, elem, elem1D, is_bound, & ! (in) + use_delta_form ) ! (in) use scale_atm_dyn_dgm_hevi_common_linalgebra, only: & atm_dyn_dgm_hevi_common_linalgebra_get_param implicit none @@ -74,7 +74,7 @@ subroutine atm_phy_bl_dgm_common_calc_tendency( & 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) + type(LocalMeshFieldBaseList), intent(in) :: QTRC_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) @@ -90,7 +90,7 @@ subroutine atm_phy_bl_dgm_common_calc_tendency( & class(ElementBase2D), pointer :: elem2D integer :: iq - real(RP) :: RHOQ00_(elem%Np,QA,lmesh%Ne) + real(RP) :: QTRC00_(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) @@ -130,7 +130,7 @@ subroutine atm_phy_bl_dgm_common_calc_tendency( & !$omp do collapse(2) do iq = 1, QA do ke=lmesh%NeS, lmesh%NeE - RHOQ00_(:,iq,ke) = RHOQ_list(iq)%ptr%val(:,ke) + QTRC00_(:,iq,ke) = QTRC_list(iq)%ptr%val(:,ke) end do end do !$omp end do @@ -147,7 +147,7 @@ subroutine atm_phy_bl_dgm_common_calc_tendency( & 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) + PROG_VARS(:,ke_xy,ke_z,3+iq) = DENS(:,ke) * QTRC00_(:,iq,ke) end do do p=1, elem%Np @@ -159,7 +159,7 @@ subroutine atm_phy_bl_dgm_common_calc_tendency( & !$omp end parallel 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) + PROG_VARS, MOMX_, MOMY_, PT_, QTRC00_, 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) @@ -180,7 +180,7 @@ subroutine atm_phy_bl_dgm_common_calc_tendency( & 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 + RHOQ_tp_list(iq)%ptr%val(:,ke) = ( PROG_VARS(:,ke_xy,ke_z,3+iq) - DENS(:,ke) * QTRC00_(:,iq,ke) ) * r_impl_fac end do end do end do @@ -788,7 +788,7 @@ end subroutine construct_sip_face_blocks_lgl !OCL SERIAL 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, & + PROG_VARS, MOMX00, MOMY00, PT00, QTRC00, 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 @@ -804,7 +804,7 @@ subroutine eval_Ax( MOMX_t, MOMY_t, DRHOT_t, RHOQ_t_list, alph_M, alph_H, & 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) :: QTRC00(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) @@ -964,7 +964,7 @@ subroutine eval_Ax( MOMX_t, MOMY_t, DRHOT_t, RHOQ_t_list, alph_M, alph_H, & 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) + + DENS(p,ke) * QTRC00(p,iq,ke) end do end do end do @@ -990,7 +990,7 @@ subroutine eval_Ax( MOMX_t, MOMY_t, DRHOT_t, RHOQ_t_list, alph_M, alph_H, & 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) + b(i,pv,iv,j,ke) = DENS(p,ke) * QTRC00(p,iq,ke) 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 e4c418ad..01ecf82c 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 @@ -234,8 +234,8 @@ subroutine atm_phy_bl_dgm_mynn_lv2_cal_VViscDiffCoef( & lmesh, elem, is_bound ) ! (in) !$omp parallel & - !$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 private( ke2D, 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 diff --git a/FElib/src/radiation/scale_atm_phy_rd_dgm_simple.F90 b/FElib/src/radiation/scale_atm_phy_rd_dgm_simple.F90 index 034f83b6..11a88586 100644 --- a/FElib/src/radiation/scale_atm_phy_rd_dgm_simple.F90 +++ b/FElib/src/radiation/scale_atm_phy_rd_dgm_simple.F90 @@ -433,7 +433,8 @@ subroutine calc_optical_thick( this, dtau_lw, dtau_sw, dtau_sw_kadd, & A = this%OPTDEP_VALLIS_EQ5_PARAMS(1) * this%OPTDEP_VALLIS_EQ5_PARAMS(2) B = this%OPTDEP_VALLIS_EQ5_PARAMS(3) C = this%OPTDEP_VALLIS_EQ5_PARAMS(4) - + dtau_sw_kadd = 0.0_RP + do ke_z=1, NeZ do p_z=1, Nnode_v-1 qv_tmp = 0.5_RP * ( qv(p_z,ke_z) + qv(p_z+1,ke_z) ) 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 43fbc5cb..cb82eb6c 100644 --- a/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_bl.F90 +++ b/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_bl.F90 @@ -253,7 +253,7 @@ 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) + type(LocalMeshFieldBaseList) :: QTRC_list(QA) class(LocalMeshFieldBase), pointer :: DENS_tp, MOMX_tp, MOMY_tp, MOMZ_tp, RHOT_tp, RHOH_P type(LocalMeshFieldBaseList) :: RHOQ_tp(QA) @@ -296,7 +296,7 @@ subroutine AtmosPhyBl_calc_tendency( & call AtmosVars_GetLocalMeshQTRCVarList( n, & mesh, trcvars_list, & - 1, RHOQ_list ) + 1, QTRC_list ) call AtmosPhyBLVars_GetLocalMeshFields_tend( n, & mesh, this%vars%tends_manager, & @@ -324,7 +324,7 @@ subroutine AtmosPhyBl_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) + QTRC_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) @@ -365,8 +365,8 @@ subroutine AtmosPhyBl_calc_tendency( & 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) + 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 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 4fc5a8cc..ff9ad0ee 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 @@ -188,8 +188,8 @@ subroutine AtmosPhyBlVars_Setup( this, model_mesh ) 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' - qtrc_tp_vinfo_tmp%DESC = 'tendency of rho*'//trim(TRACER_NAME(this%QS+iq-1))//' in BL process' + qtrc_tp_vinfo_tmp%NAME = 'BL_'//trim(TRACER_NAME(iq))//'_t' + qtrc_tp_vinfo_tmp%DESC = 'tendency of rho*'//trim(TRACER_NAME(iq))//' in BL process' qtrc_tp_vinfo_tmp%UNIT = 'kg/m3/s' reg_file_hist = .true. 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 d4bd8ab0..ccca7155 100644 --- a/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_sfc.F90 +++ b/model/atm_nonhydro3d/src/atmos/mod_atmos_phy_sfc.F90 @@ -489,7 +489,7 @@ subroutine cal_tend_from_sfcflx( this, is_update_sflx, & call convert_LocalOrth2UVVec( & this%mesh%ptr_mesh, lcmesh2D%pos_en(:,:,1), lcmesh2D%pos_en(:,:,2), Z1(:,:), elem2D%Np*lcmesh2D%Ne, & - SFLX_MU, SFLX_MV ) ! (inout) + SFLX_MU(:,lcmesh2D%NeS:lcmesh2D%NeE), SFLX_MV(:,lcmesh2D%NeS:lcmesh2D%NeE) ) ! (inout) end if @@ -551,12 +551,12 @@ subroutine cal_del_flux( del_flux, & class(LocalMesh2D), intent(in) :: lmesh2D class(ElementBase2D), intent(in) :: elem2D real(RP), intent(out) :: del_flux(elem%NfpTot*lmesh%Ne,5) - real(RP), intent(in) :: sflx_mu(elem2D%Np,lmesh2D%Ne) - real(RP), intent(in) :: sflx_mv(elem2D%Np,lmesh2D%Ne) - real(RP), intent(in) :: sflx_mw(elem2D%Np,lmesh2D%Ne) - real(RP), intent(in) :: sflx_sh(elem2D%Np,lmesh2D%Ne) - real(RP), intent(in) :: sflx_qv(elem2D%Np,lmesh2D%Ne) - real(RP), intent(in) :: SFC_TEMP(elem2D%Np,lmesh2D%Ne) + real(RP), intent(in) :: sflx_mu(elem2D%Np,lmesh2D%NeA) + real(RP), intent(in) :: sflx_mv(elem2D%Np,lmesh2D%NeA) + real(RP), intent(in) :: sflx_mw(elem2D%Np,lmesh2D%NeA) + real(RP), intent(in) :: sflx_sh(elem2D%Np,lmesh2D%NeA) + real(RP), intent(in) :: sflx_qv(elem2D%Np,lmesh2D%NeA) + real(RP), intent(in) :: SFC_TEMP(elem2D%Np,lmesh2D%NeA) real(RP), intent(in) :: nz(elem%NfpTot*lmesh%Ne) integer :: ke2D, p diff --git a/model/atm_nonhydro3d/src/ocean/mod_ocean_mesh_gm.F90 b/model/atm_nonhydro3d/src/ocean/mod_ocean_mesh_gm.F90 index b1a3c9e6..e3fa7d7e 100644 --- a/model/atm_nonhydro3d/src/ocean/mod_ocean_mesh_gm.F90 +++ b/model/atm_nonhydro3d/src/ocean/mod_ocean_mesh_gm.F90 @@ -195,7 +195,9 @@ subroutine OceanMeshGM_Init( this ) this%mesh%mesh2D, dom_zmin ) this%vcoord_type_id = MeshUtil_get_VCoord_TypeID( VERTICAL_COORD_NAME ) - call this%Setup_vcoordinate() + ! We tentatively disable the setup of vertical coordinate in ocean component. + ! For future, we will implement a more flexible setup of vertical coordinate. + ! call this%Setup_vcoordinate() !- this%comm_use_mpi_pc = COMM_USE_MPI_PC From a3f10679b519692cf2499a36d142c6a61f15e343 Mon Sep 17 00:00:00 2001 From: ykawai Date: Wed, 5 Aug 2026 15:15:38 +0900 Subject: [PATCH 2/6] Add water vapor flux variable to ocean component --- .../src/coupler/mod_cpl_vars.F90 | 42 ++++--- .../src/ocean/mod_ocean_component.F90 | 114 +++++++++++++++--- .../src/ocean/mod_ocean_vars.F90 | 3 +- 3 files changed, 124 insertions(+), 35 deletions(-) diff --git a/model/atm_nonhydro3d/src/coupler/mod_cpl_vars.F90 b/model/atm_nonhydro3d/src/coupler/mod_cpl_vars.F90 index 424e0eb4..54724094 100644 --- a/model/atm_nonhydro3d/src/coupler/mod_cpl_vars.F90 +++ b/model/atm_nonhydro3d/src/coupler/mod_cpl_vars.F90 @@ -82,8 +82,9 @@ module mod_cpl_vars integer, parameter, public :: OCN_SFLX_MV_ID = 4 integer, parameter, public :: OCN_SFLX_SH_ID = 5 integer, parameter, public :: OCN_SFLX_LH_ID = 6 - integer, parameter, public :: OCN_SFC_ALBEDO_ID = 7 - integer, parameter, public :: OCN_VAR_NUM = 7 + integer, parameter, public :: OCN_SFLX_QV_ID = 7 + integer, parameter, public :: OCN_SFC_ALBEDO_ID = 8 + integer, parameter, public :: OCN_VAR_NUM = 8 ! Output to ocean component integer, parameter, public :: OCN_ATM_SFC_DENS_ID = 1 @@ -226,7 +227,8 @@ subroutine CouplerVars_putOCN( this, ocn_vars, countup ) SFLX_MU_ID => OCN_SFLX_MU_ID, & SFLX_MV_ID => OCN_SFLX_MV_ID, & SFLX_SH_ID => OCN_SFLX_SH_ID, & - SFLX_LH_ID => OCN_SFLX_LH_ID + SFLX_LH_ID => OCN_SFLX_LH_ID, & + SFLX_QV_ID => OCN_SFLX_QV_ID implicit none class(CouplerVars), intent(inout) :: this type(OceanVars), intent(in) :: ocn_vars @@ -246,13 +248,15 @@ subroutine CouplerVars_putOCN( this, ocn_vars, countup ) ocn_vars%OCN_SFLX(SFLX_MV_ID)%local(ldomID)%val, & ocn_vars%OCN_SFLX(SFLX_SH_ID)%local(ldomID)%val, & ocn_vars%OCN_SFLX(SFLX_LH_ID)%local(ldomID)%val, & + ocn_vars%OCN_SFLX(SFLX_QV_ID)%local(ldomID)%val, & !- this%OCN_vars(OCN_SFC_TEMP_ID)%local(ldomID)%val, & this%OCN_vars(OCN_SFLX_MW_ID)%local(ldomID)%val, & this%OCN_vars(OCN_SFLX_MU_ID)%local(ldomID)%val, & this%OCN_vars(OCN_SFLX_MV_ID)%local(ldomID)%val, & this%OCN_vars(OCN_SFLX_SH_ID)%local(ldomID)%val, & - this%OCN_vars(OCN_SFLX_LH_ID)%local(ldomID)%val ) + this%OCN_vars(OCN_SFLX_LH_ID)%local(ldomID)%val, & + this%OCN_vars(OCN_SFLX_QV_ID)%local(ldomID)%val ) end do ! Update counter @@ -265,9 +269,10 @@ end subroutine CouplerVars_putOCN !> Get SFC variables from CPL buffer for ATM component !OCL SERIAL - subroutine CouplerVars_get_SFC_ATM( this, & - SFC_TEMP, SFC_ALBEDO, & - SFLX_MW, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_LH ) + subroutine CouplerVars_get_SFC_ATM( this, & + SFC_TEMP, SFC_ALBEDO, & + SFLX_MW, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_LH, & + SFLX_QV ) implicit none class(CouplerVars), intent(inout) :: this type(MeshField2D), intent(inout) :: SFC_TEMP @@ -277,6 +282,7 @@ subroutine CouplerVars_get_SFC_ATM( this, & type(MeshField2D), intent(inout) :: SFLX_MV type(MeshField2D), intent(inout) :: SFLX_SH type(MeshField2D), intent(inout) :: SFLX_LH + type(MeshField2D), intent(inout) :: SFLX_QV integer :: ldomID class(LocalMesh2D), pointer :: lmesh, lmesh_o @@ -295,11 +301,12 @@ subroutine CouplerVars_get_SFC_ATM( this, & this%OCN_vars(OCN_SFLX_MV_ID)%local(ldomID)%val, & this%OCN_vars(OCN_SFLX_SH_ID)%local(ldomID)%val, & this%OCN_vars(OCN_SFLX_LH_ID)%local(ldomID)%val, & + this%OCN_vars(OCN_SFLX_QV_ID)%local(ldomID)%val, & !- SFC_TEMP%local(ldomID)%val, SFC_ALBEDO%local(ldomID)%val, & SFLX_MW%local(ldomID)%val, SFLX_MU%local(ldomID)%val, & SFLX_MV%local(ldomID)%val, SFLX_SH%local(ldomID)%val, & - SFLX_LH%local(ldomID)%val ) + SFLX_LH%local(ldomID)%val, SFLX_QV%local(ldomID)%val ) end do return end subroutine CouplerVars_get_SFC_ATM @@ -395,7 +402,7 @@ subroutine putATM_local( this, & real(RP), intent(out) :: O_ATM_SFLX_ENGI(elem2D_a%Np,lmesh2D_a%NeA) integer :: ke2D - integer :: coef1, coef2 + real(RP) :: coef1, coef2 integer :: hSlice0(elem2D_a%Np), hSlice1(elem2D_a%Np) real(RP) :: dens_(elem2D_a%Np) @@ -419,7 +426,7 @@ subroutine putATM_local( this, & O_ATM_TEMP(:,ke2D) = O_ATM_TEMP(:,ke2D) * coef1 + temp_(:) O_ATM_PRES(:,ke2D) = O_ATM_PRES(:,ke2D) * coef1 + ATM_PRES(hSlice0(:),ke2D) - O_ATM_W(:,ke2D) = O_ATM_W(:,ke2D) * coef1 + ATM_MOMZ(hSlice0(:),ke2D) / dens_(:) + O_ATM_W(:,ke2D) = O_ATM_W(:,ke2D) * coef1 + ATM_MOMZ(hSlice1(:),ke2D) / dens_(:) O_ATM_U(:,ke2D) = O_ATM_U(:,ke2D) * coef1 + ATM_MOMX(hSlice0(:),ke2D) / dens_(:) O_ATM_V(:,ke2D) = O_ATM_V(:,ke2D) * coef1 + ATM_MOMY(hSlice0(:),ke2D) / dens_(:) O_ATM_QV(:,ke2D) = O_ATM_QV(:,ke2D) * coef1 + ATM_QV(hSlice0(:),ke2D) @@ -451,8 +458,8 @@ end subroutine putATM_local !OCL SERIAL subroutine putOCN_local( this, & lmesh_o, elem_o, & - SFC_TEMP, SFLX_MW, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_LH, & - SFC_TEMP_out, SFLX_MW_out, SFLX_MU_out, SFLX_MV_out, SFLX_SH_out, SFLX_LH_out ) + SFC_TEMP, SFLX_MW, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_LH, SFLX_QV, & + SFC_TEMP_out, SFLX_MW_out, SFLX_MU_out, SFLX_MV_out, SFLX_SH_out, SFLX_LH_out, SFLX_QV_out ) implicit none class(CouplerVars), intent(inout) :: this class(LocalMesh2D), intent(in) :: lmesh_o @@ -463,12 +470,14 @@ subroutine putOCN_local( this, & real(RP), intent(in) :: SFLX_MV(elem_o%Np,lmesh_o%NeA) real(RP), intent(in) :: SFLX_SH(elem_o%Np,lmesh_o%NeA) real(RP), intent(in) :: SFLX_LH(elem_o%Np,lmesh_o%NeA) + real(RP), intent(in) :: SFLX_QV(elem_o%Np,lmesh_o%NeA) real(RP), intent(out) :: SFC_TEMP_out(elem_o%Np,lmesh_o%NeA) real(RP), intent(out) :: SFLX_MW_out(elem_o%Np,lmesh_o%NeA) real(RP), intent(out) :: SFLX_MU_out(elem_o%Np,lmesh_o%NeA) real(RP), intent(out) :: SFLX_MV_out(elem_o%Np,lmesh_o%NeA) real(RP), intent(out) :: SFLX_SH_out(elem_o%Np,lmesh_o%NeA) real(RP), intent(out) :: SFLX_LH_out(elem_o%Np,lmesh_o%NeA) + real(RP), intent(out) :: SFLX_QV_out(elem_o%Np,lmesh_o%NeA) integer :: ke real(RP) :: coef1, coef2 @@ -486,6 +495,7 @@ subroutine putOCN_local( this, & SFLX_MV_out(:,ke) = SFLX_MV_out(:,ke) * coef1 + SFLX_MV(:,ke) SFLX_SH_out(:,ke) = SFLX_SH_out(:,ke) * coef1 + SFLX_SH(:,ke) SFLX_LH_out(:,ke) = SFLX_LH_out(:,ke) * coef1 + SFLX_LH(:,ke) + SFLX_QV_out(:,ke) = SFLX_QV_out(:,ke) * coef1 + SFLX_QV(:,ke) end do !$omp do do ke=lmesh_o%NeS, lmesh_o%NeE @@ -495,6 +505,7 @@ subroutine putOCN_local( this, & SFLX_MV_out(:,ke) = SFLX_MV_out(:,ke) * coef2 SFLX_SH_out(:,ke) = SFLX_SH_out(:,ke) * coef2 SFLX_LH_out(:,ke) = SFLX_LH_out(:,ke) * coef2 + SFLX_QV_out(:,ke) = SFLX_QV_out(:,ke) * coef2 end do !$omp end parallel return @@ -503,8 +514,8 @@ end subroutine putOCN_local !OCL SERIAL subroutine get_SFC_ATM_local( this, & lmesh_a, elem_a, lmesh_o, elem_o, & - O_SFC_TEMP, O_SFC_ALBEDO, O_SFLX_MW, O_SFLX_MU, O_SFLX_MV, O_SFLX_SH, O_SFLX_LH, & - SFC_TEMP, SFC_ALBEDO, SFLX_MW, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_LH ) + O_SFC_TEMP, O_SFC_ALBEDO, O_SFLX_MW, O_SFLX_MU, O_SFLX_MV, O_SFLX_SH, O_SFLX_LH, O_SFLX_QV, & + SFC_TEMP, SFC_ALBEDO, SFLX_MW, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_LH, SFLX_QV ) implicit none class(CouplerVars), intent(inout) :: this class(LocalMesh2D), intent(in) :: lmesh_a @@ -518,6 +529,7 @@ subroutine get_SFC_ATM_local( this, & real(RP), intent(in) :: O_SFLX_MV(elem_o%Np,lmesh_o%NeA) real(RP), intent(in) :: O_SFLX_SH(elem_o%Np,lmesh_o%NeA) real(RP), intent(in) :: O_SFLX_LH(elem_o%Np,lmesh_o%NeA) + real(RP), intent(in) :: O_SFLX_QV(elem_o%Np,lmesh_o%NeA) real(RP), intent(out) :: SFC_TEMP(elem_a%Np,lmesh_a%NeA) real(RP), intent(out) :: SFC_ALBEDO(elem_a%Np,lmesh_a%NeA) real(RP), intent(out) :: SFLX_MW(elem_a%Np,lmesh_a%NeA) @@ -525,6 +537,7 @@ subroutine get_SFC_ATM_local( this, & real(RP), intent(out) :: SFLX_MV(elem_a%Np,lmesh_a%NeA) real(RP), intent(out) :: SFLX_SH(elem_a%Np,lmesh_a%NeA) real(RP), intent(out) :: SFLX_LH(elem_a%Np,lmesh_a%NeA) + real(RP), intent(out) :: SFLX_QV(elem_a%Np,lmesh_a%NeA) integer :: ke !---------------------------------------------------------------- @@ -538,6 +551,7 @@ subroutine get_SFC_ATM_local( this, & SFLX_MV(:,ke) = O_SFLX_MV(:,ke) SFLX_SH(:,ke) = O_SFLX_SH(:,ke) SFLX_LH(:,ke) = O_SFLX_LH(:,ke) + SFLX_QV(:,ke) = O_SFLX_QV(:,ke) end do return end subroutine get_SFC_ATM_local diff --git a/model/atm_nonhydro3d/src/ocean/mod_ocean_component.F90 b/model/atm_nonhydro3d/src/ocean/mod_ocean_component.F90 index b580dd0a..6c590865 100644 --- a/model/atm_nonhydro3d/src/ocean/mod_ocean_component.F90 +++ b/model/atm_nonhydro3d/src/ocean/mod_ocean_component.F90 @@ -242,6 +242,7 @@ subroutine Ocean_calc_tendency( this, force ) SFLX_MV_ID => OCN_SFLX_MV_ID, & SFLX_SH_ID => OCN_SFLX_SH_ID, & SFLX_LH_ID => OCN_SFLX_LH_ID, & + SFLX_QV_ID => OCN_SFLX_QV_ID, & SFLX_RD_SW_DIR_ID => ATMVAR2D_SFLX_RD_SW_DIR_ID, & SFLX_RD_LW_DIF_ID => ATMVAR2D_SFLX_RD_LW_DIF_ID, & RHOH_ID => PHYTEND_RHOH_ID @@ -272,12 +273,13 @@ subroutine Ocean_calc_tendency( this, force ) !- Calculate surface fluxes - call calculate_surface_flux( & + call calculate_surface_flux( this, & this%vars%OCN_SFLX(SFLX_MW_ID)%local(idom)%val, & this%vars%OCN_SFLX(SFLX_MU_ID)%local(idom)%val, & this%vars%OCN_SFLX(SFLX_MV_ID)%local(idom)%val, & this%vars%OCN_SFLX(SFLX_SH_ID)%local(idom)%val, & this%vars%OCN_SFLX(SFLX_LH_ID)%local(idom)%val, & + this%vars%OCN_SFLX(SFLX_QV_ID)%local(idom)%val, & !- this%vars%AUX_VARS2D(SFC_TEMP_ID)%local(idom)%val, & this%vars%ATM_VARS2D(ATM_SFC_DENS_ID)%local(idom)%val, & @@ -292,13 +294,15 @@ subroutine Ocean_calc_tendency( this, force ) this%vars%ATM_VARS2D(DZ_A_ID)%local(idom)%val, & lmesh2D, lmesh2D%refElem2D ) - !- Calculate ground heat flux + !- Calculate ocean surface heat flux - call calculate_ground_heat_flux( SFLX_GH, & ! (out) + call calculate_sfc_heat_flux( SFLX_GH, & ! (out) this%vars%AUX_VARS2D(SFC_TEMP_ID)%local(idom)%val, & ! (in) this%vars%ATM_VARS2D(SFLX_RD_SW_DIR_ID)%local(idom)%val, & ! (in) this%vars%ATM_VARS2D(SFLX_RD_LW_DIF_ID)%local(idom)%val, & ! (in) this%vars%AUX_VARS2D(ALB_VIS_DIR_ID)%local(idom)%val, & ! (in) + this%vars%OCN_SFLX(SFLX_SH_ID)%local(idom)%val, & ! (in) + this%vars%OCN_SFLX(SFLX_LH_ID)%local(idom)%val, & ! (in) lmesh2D, lmesh2D%refElem2D ) ! (in) !- Calculate tendencies with physics @@ -475,7 +479,8 @@ end subroutine Ocean_finalize !> Calculate momentum and heat flux at the surface (tentative) !OCL SERIAL - subroutine calculate_surface_flux( SFLX_MW, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_LH, & + subroutine calculate_surface_flux( this, & + SFLX_MW, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_LH, SFLX_QV, & SFC_TEMP, SFC_DENS, SFC_PRES, & ATM_TEMP, ATM_DENS, ATM_PRES, ATM_W, ATM_U, ATM_V, ATM_QV, zlev_a, & lmesh, elem ) @@ -484,6 +489,7 @@ subroutine calculate_surface_flux( SFLX_MW, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_LH, use scale_atm_phy_sf_bulk_simple, only: & ATMOS_PHY_SF_simple_flux implicit none + class(OceanComponent), intent(in) :: this class(LocalMesh2D), intent(in) :: lmesh class(ElementBase2D), intent(in) :: elem real(RP), intent(out) :: SFLX_MW(elem%Np,lmesh%NeA) @@ -491,6 +497,7 @@ subroutine calculate_surface_flux( SFLX_MW, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_LH, real(RP), intent(out) :: SFLX_MV(elem%Np,lmesh%NeA) real(RP), intent(out) :: SFLX_SH(elem%Np,lmesh%NeA) real(RP), intent(out) :: SFLX_LH(elem%Np,lmesh%NeA) + real(RP), intent(out) :: SFLX_QV(elem%Np,lmesh%NeA) real(RP), intent(in) :: SFC_TEMP(elem%Np,lmesh%NeA) real(RP), intent(in) :: SFC_DENS(elem%Np,lmesh%NeA) real(RP), intent(in) :: SFC_PRES(elem%Np,lmesh%NeA) @@ -503,8 +510,6 @@ subroutine calculate_surface_flux( SFLX_MW, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_LH, real(RP), intent(in) :: ATM_QV(elem%Np,lmesh%NeA) real(RP), intent(in) :: zlev_a(elem%Np,lmesh%NeA) - real(RP) :: SFLX_QV(elem%Np,lmesh%NeA) - ! dummy real(RP) :: U10(elem%Np,lmesh%NeA) real(RP) :: V10(elem%Np,lmesh%NeA) @@ -514,38 +519,54 @@ subroutine calculate_surface_flux( SFLX_MW, SFLX_MU, SFLX_MV, SFLX_SH, SFLX_LH, real(RP) :: DZ1 (elem%Np,lmesh%NeA) real(RP) :: Z1 (elem%Np,lmesh%NeA) + real(RP) :: ATM_W_lc(elem%Np,lmesh%Ne) + real(RP) :: ATM_U_lc(elem%Np,lmesh%Ne) + real(RP) :: ATM_V_lc(elem%Np,lmesh%Ne) !-------------------------------------------------- !$omp parallel do collapse(2) do ke2D=lmesh%NeS, lmesh%NeE do ij=1, elem%Np + ATM_W_lc(ij,ke2D) = ATM_W(ij,ke2D) + ATM_U_lc(ij,ke2D) = ATM_U(ij,ke2D) + ATM_V_lc(ij,ke2D) = ATM_V(ij,ke2D) + Z1(ij,ke2D) = zlev_a(ij,ke2D) DZ1(ij,ke2D) = Z1(ij,ke2D) - zlev_a(ij,ke2D) Z1(ij,ke2D) = Z1(ij,ke2D) + RPlanet end do end do - call ATMOS_PHY_SF_simple_flux( & - elem%Np, 1, elem%Np, lmesh%NeA, 1, lmesh%Ne, & ! (in) - ATM_W(:,:), ATM_U(:,:), ATM_V(:,:), ATM_TEMP(:,:), ATM_PRES(:,:), & ! (in) Note: ATM_PRES is not used - ATM_QV(:,:), & ! (in) - SFC_DENS(:,:), SFC_TEMP(:,:), SFC_PRES(:,:), & ! (in) - DZ1(:,:), & ! (in) - SFLX_MW(:,:), SFLX_MU(:,:), SFLX_MV(:,:), & ! (out) - SFLX_SH(:,:), SFLX_LH(:,:), SFLX_QV(:,:), & ! (out) - U10(:,:), V10(:,:) ) ! (out) + call convert_UV2LocalOrthVec( & + this%mesh%ptr_mesh, lmesh%pos_en(:,:,1), lmesh%pos_en(:,:,2), Z1(:,:), elem%Np*lmesh%Ne, & + ATM_U_lc, ATM_V_lc ) ! (inout) + call ATMOS_PHY_SF_simple_flux( & + elem%Np, 1, elem%Np, lmesh%NeA, 1, lmesh%Ne, & ! (in) + ATM_W_lc(:,:), ATM_U_lc(:,:), ATM_V_lc(:,:), ATM_TEMP(:,:), ATM_PRES(:,:), & ! (in) Note: ATM_PRES is not used + ATM_QV(:,:), & ! (in) + SFC_DENS(:,:), SFC_TEMP(:,:), SFC_PRES(:,:), & ! (in) + DZ1(:,:), & ! (in) + SFLX_MW(:,:), SFLX_MU(:,:), SFLX_MV(:,:), & ! (out) + SFLX_SH(:,:), SFLX_LH(:,:), SFLX_QV(:,:), & ! (out) + U10(:,:), V10(:,:) ) ! (out) + + call convert_LocalOrth2UVVec( & + this%mesh%ptr_mesh, lmesh%pos_en(:,:,1), lmesh%pos_en(:,:,2), Z1(:,:), elem%Np*lmesh%Ne, & + SFLX_MU(:,lmesh%NeS:lmesh%NeE), SFLX_MV(:,lmesh%NeS:lmesh%NeE) ) ! (inout) + return end subroutine calculate_surface_flux !> Calculate ground heat flux (tentative) !! !OCL SERIAL - subroutine calculate_ground_heat_flux( sflx_GH, & - SFC_TEMP, & + subroutine calculate_sfc_heat_flux( sflx_GH, & + SFC_TEMP, & SFLX_RD_SW_dn_dir, SFLX_RD_LW_dn_dif, & - SFC_ALB_dir_vis, & + SFC_ALB_dir_vis, & + SFLX_SH, SFLX_LH, & lmesh, elem ) use scale_const, only: & STB => CONST_STB @@ -557,6 +578,8 @@ subroutine calculate_ground_heat_flux( sflx_GH, & real(RP), intent(in) :: SFLX_RD_SW_dn_dir(elem%Np,lmesh%NeA) real(RP), intent(in) :: SFLX_RD_LW_dn_dif(elem%Np,lmesh%NeA) real(RP), intent(in) :: SFC_ALB_dir_vis(elem%Np,lmesh%NeA) + real(RP), intent(in) :: SFLX_SH(elem%Np,lmesh%NeA) + real(RP), intent(in) :: SFLX_LH(elem%Np,lmesh%NeA) integer :: ke real(RP) :: emis(elem%Np) @@ -572,10 +595,10 @@ subroutine calculate_ground_heat_flux( sflx_GH, & SWD(:) = SFLX_RD_SW_dn_dir(:,ke) SWU(:) = SFC_ALB_dir_vis(:,ke) * SWD(:) - sflx_GH(:,ke) = SWD(:) - SWU(:) + LWD(:) - LWU(:) + sflx_GH(:,ke) = SWD(:) - SWU(:) + LWD(:) - LWU(:) - SFLX_SH(:,ke) - SFLX_LH(:,ke) end do return - end subroutine calculate_ground_heat_flux + end subroutine calculate_sfc_heat_flux !OCL SERIAL subroutine calculate_phys_tendency( RHOH, & @@ -618,4 +641,55 @@ subroutine calculate_phys_tendency( RHOH, & !$omp end parallel return end subroutine calculate_phys_tendency + +!OCL SERIAL + subroutine convert_UV2LocalOrthVec( mesh, x1, x2, r, N, & ! (in) + U, V ) ! (inout) + use scale_cubedsphere_coord_cnv, only: & + CubedSphereCoordCnv_CS2LocalOrthVec_alpha + use scale_mesh_cubedspheredom3d, only: MeshCubedSphereDom3D + + implicit none + integer, intent(in) :: N + class(MeshBase3D), intent(in) :: mesh + real(RP), intent(in) :: x1(N) + real(RP), intent(in) :: x2(N) + real(RP), intent(in) :: r(N) + real(RP), intent(inout) :: U(N) + real(RP), intent(inout) :: V(N) + + select type(mesh) + type is (MeshCubedSphereDom3D) + call CubedSphereCoordCnv_CS2LocalOrthVec_alpha( & + x1, x2, r, N, & ! (in) + U, V ) ! (inout) + end select + return + end subroutine convert_UV2LocalOrthVec + +!OCL SERIAL + subroutine convert_LocalOrth2UVVec( mesh, x1, x2, r, N, & ! (in) + U, V ) ! (inout) + use scale_cubedsphere_coord_cnv, only: & + CubedSphereCoordCnv_LocalOrth2CSVec_alpha + use scale_mesh_cubedspheredom3d, only: MeshCubedSphereDom3D + + implicit none + integer, intent(in) :: N + class(MeshBase3D), intent(in) :: mesh + real(RP), intent(in) :: x1(N) + real(RP), intent(in) :: x2(N) + real(RP), intent(in) :: r(N) + real(RP), intent(inout) :: U(N) + real(RP), intent(inout) :: V(N) + + select type(mesh) + type is (MeshCubedSphereDom3D) + call CubedSphereCoordCnv_LocalOrth2CSVec_alpha( & + x1, x2, r, N, & ! (in) + U, V ) ! (inout) + end select + return + end subroutine convert_LocalOrth2UVVec + end module mod_ocean_component \ No newline at end of file diff --git a/model/atm_nonhydro3d/src/ocean/mod_ocean_vars.F90 b/model/atm_nonhydro3d/src/ocean/mod_ocean_vars.F90 index d977604a..ee816e85 100644 --- a/model/atm_nonhydro3d/src/ocean/mod_ocean_vars.F90 +++ b/model/atm_nonhydro3d/src/ocean/mod_ocean_vars.F90 @@ -156,7 +156,8 @@ module mod_ocean_vars integer, public, parameter :: OCN_SFLX_MV_ID = 3 !< v-momentum flux at the ocean surface [kg/m/s2] integer, public, parameter :: OCN_SFLX_SH_ID = 4 !< sensible heat flux at the ocean surface [W/m2] integer, public, parameter :: OCN_SFLX_LH_ID = 5 !< latent heat flux at the ocean surface [W/m2] - integer, public, parameter :: OCN_SFLX_NUM = 5 !< Number of ocean surface fluxes + integer, public, parameter :: OCN_SFLX_QV_ID = 6 !< water vapor flux at the ocean surface [kg/m2/s] + integer, public, parameter :: OCN_SFLX_NUM = 6 !< Number of ocean surface fluxes ! Diagnostic variables From f6a34ae8b46de01d4990a0f029f79d26bb6186f7 Mon Sep 17 00:00:00 2001 From: ykawai Date: Wed, 5 Aug 2026 15:17:25 +0900 Subject: [PATCH 3/6] Consider CP component when precipitation fluxes is output --- .../src/atmos/mod_atmos_component.F90 | 26 ++++- .../src/atmos/mod_atmos_phy_cp_vars.F90 | 31 ++++++ .../src/atmos/mod_atmos_vars.F90 | 99 ++++++++++++++++--- 3 files changed, 136 insertions(+), 20 deletions(-) diff --git a/model/atm_nonhydro3d/src/atmos/mod_atmos_component.F90 b/model/atm_nonhydro3d/src/atmos/mod_atmos_component.F90 index 1c28ecf1..fa3f6e6e 100644 --- a/model/atm_nonhydro3d/src/atmos/mod_atmos_component.F90 +++ b/model/atm_nonhydro3d/src/atmos/mod_atmos_component.F90 @@ -283,7 +283,7 @@ subroutine Atmos_setup_vars( this ) 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%Regist_physvar_manager( mp_AUXVARS2D_manager=this%phy_mp_proc%vars%auxvars2D_manager ) call this%vars%Setup_container( this%phy_mp_proc%atm_var_container_typeid, this%mesh ) end if !- Surface component @@ -309,6 +309,7 @@ subroutine Atmos_setup_vars( this ) !- Cumulus parameterization component if ( this%phy_cp_proc%IsActivated() ) then call this%phy_cp_proc%vars%Setup( this%mesh ) + call this%vars%Regist_physvar_manager( cp_AUXVARS2D_manager=this%phy_cp_proc%vars%auxvars2D_manager ) end if !- PBL component if ( this%phy_bl_proc%IsActivated() ) then @@ -687,6 +688,8 @@ subroutine Atmos_set_surface( this, countup ) PREC_ENGI_ID => ATMOS_AUXVARS2D_PREC_ENGI_ID use mod_atmos_phy_mp_vars, only: & AtmosPhyMpVars_GetLocalMeshFields_sfcflx + use mod_atmos_phy_cp_vars, only: & + AtmosPhyCpVars_GetLocalMeshFields_sfcflx use mod_atmos_phy_rd_vars, only: & RD_SFLX_LW_dif_ID => ATMOS_PHY_RD_AUX2D_SFLX_LW_dn_ID, & RD_SFLX_SW_dir_ID => ATMOS_PHY_RD_AUX2D_SFLX_SW_dn_ID @@ -703,6 +706,7 @@ subroutine Atmos_set_surface( this, countup ) class(LocalMeshFieldBase), pointer :: PREC, PREC_ENGI class(LocalMeshFieldBase), pointer :: SFLX_rain_MP, SFLX_snow_MP, SFLX_ENGI_MP + class(LocalMeshFieldBase), pointer :: SFLX_rain_CP, SFLX_snow_CP, SFLX_ENGI_CP integer :: iq !-------------------------------------------------- @@ -715,7 +719,7 @@ subroutine Atmos_set_surface( this, countup ) call mesh%GetMesh2D( mesh2D ) end select - !- sum of rainfall from mp and cp + !- Sum up precipitation and energy fluxes from cloud microphysics and cumulus parameterization components do n=1, mesh2D%LOCAL_MESH_NUM call AtmosVars_GetLocalMeshSfcVar( n, & @@ -739,6 +743,18 @@ subroutine Atmos_set_surface( this, countup ) PREC_ENGI%val(:,ke) = PREC_ENGI%val(:,ke) + SFLX_ENGI_MP%val(:,ke) end do end if + + if ( this%phy_cp_proc%IsActivated() ) then + call AtmosPhyCpVars_GetLocalMeshFields_sfcflx( n, & + mesh2D, this%phy_cp_proc%vars%auxvars2D_manager, & ! (in) + SFLX_rain_CP, SFLX_snow_CP, SFLX_ENGI_CP ) ! (out) + + !$omp parallel do private(ke) + do ke=lcmesh%NeS, lcmesh%NeE + PREC %val(:,ke) = PREC %val(:,ke) + SFLX_rain_CP%val(:,ke) + SFLX_snow_CP%val(:,ke) + PREC_ENGI%val(:,ke) = PREC_ENGI%val(:,ke) + SFLX_ENGI_CP%val(:,ke) + end do + end if end do if ( this%coupler_ptr%IsActivated() ) then @@ -779,7 +795,8 @@ subroutine Atmos_get_surface( this ) SFLX_MU_ID => ATMOS_PHY_SF_SFLX_MU_ID, & SFLX_MV_ID => ATMOS_PHY_SF_SFLX_MV_ID, & SFLX_SH_ID => ATMOS_PHY_SF_SFLX_SH_ID, & - SFLX_LH_ID => ATMOS_PHY_SF_SFLX_LH_ID + SFLX_LH_ID => ATMOS_PHY_SF_SFLX_LH_ID, & + SFLX_QV_ID => ATMOS_PHY_SF_SFLX_QV_ID use mod_cpl_component, only: CouplerComponent implicit none class(AtmosComponent), intent(inout), target :: this @@ -795,7 +812,8 @@ subroutine Atmos_get_surface( this ) this%phy_sfc_proc%vars%SFC_FLX(SFLX_MU_ID), & this%phy_sfc_proc%vars%SFC_FLX(SFLX_MV_ID), & this%phy_sfc_proc%vars%SFC_FLX(SFLX_SH_ID), & - this%phy_sfc_proc%vars%SFC_FLX(SFLX_LH_ID) ) + this%phy_sfc_proc%vars%SFC_FLX(SFLX_LH_ID), & + this%phy_sfc_proc%vars%SFC_FLX(SFLX_QV_ID) ) call PROF_rapend( 'ATM_sfc_exch', 1) return end subroutine Atmos_get_surface 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 4ce986a1..55cfcd3e 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 @@ -68,6 +68,7 @@ module mod_atmos_phy_cp_vars end type AtmosPhyCpVars public :: AtmosPhyCpVars_GetLocalMeshFields_tend + public :: AtmosPhyCPVars_GetLocalMeshFields_sfcflx !----------------------------------------------------------------------------- ! @@ -282,6 +283,36 @@ subroutine AtmosPhyCpVars_GetLocalMeshFields_tend( domID, mesh, bl_tends_list, & return end subroutine AtmosPhyCpVars_GetLocalMeshFields_tend +!OCL SERIAL + subroutine AtmosPhyCPVars_GetLocalMeshFields_sfcflx( domID, mesh, sfcflx_list, & + SFLX_rain, SFLX_snow, SFLX_engi ) + + 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) :: sfcflx_list + class(LocalMeshFieldBase), pointer, intent(out) :: SFLX_rain + class(LocalMeshFieldBase), pointer, intent(out) :: SFLX_snow + class(LocalMeshFieldBase), pointer, intent(out) :: SFLX_engi + + class(MeshFieldBase), pointer :: field + !------------------------------------------------------- + + call sfcflx_list%Get(ATMOS_PHY_CP_AUX2D_SFLX_RAIN_ID, field) + call field%GetLocalMeshField(domID, SFLX_rain) + + call sfcflx_list%Get(ATMOS_PHY_CP_AUX2D_SFLX_SNOW_ID, field) + call field%GetLocalMeshField(domID, SFLX_snow) + + call sfcflx_list%Get(ATMOS_PHY_CP_AUX2D_SFLX_engi_ID, field) + call field%GetLocalMeshField(domID, SFLX_engi) + + return + end subroutine AtmosPhyCPVars_GetLocalMeshFields_sfcflx + !OCL SERIAL subroutine AtmosPhyCpVars_history( this ) use scale_file_history_meshfield, only: FILE_HISTORY_meshfield_put diff --git a/model/atm_nonhydro3d/src/atmos/mod_atmos_vars.F90 b/model/atm_nonhydro3d/src/atmos/mod_atmos_vars.F90 index 847f3a1c..dc306f08 100644 --- a/model/atm_nonhydro3d/src/atmos/mod_atmos_vars.F90 +++ b/model/atm_nonhydro3d/src/atmos/mod_atmos_vars.F90 @@ -93,6 +93,7 @@ module mod_atmos_vars !- type(ModelVarManager), pointer :: ptr_MP_AUXVARS2D_manager + type(ModelVarManager), pointer :: ptr_CP_AUXVARS2D_manager logical :: moist type(MeshField3D), pointer :: QV @@ -458,17 +459,24 @@ subroutine AtmosVars_Final( this ) return end subroutine AtmosVars_Final + !> Set the pointer of 2D ModelVarManager with cloud microphysics or cumlus parameterization + !! to output surface variables with precipitation fluxes !OCL SERIAL subroutine AtmosVars_Regist_physvar_manager( this, & - mp_AUXVARS2D_manager ) + mp_AUXVARS2D_manager, cp_AUXVARS2D_manager ) implicit none class(AtmosVars), target, intent(inout) :: this - type(ModelVarManager), intent(in), target :: mp_AUXVARS2D_manager + type(ModelVarManager), intent(in), target, optional:: mp_AUXVARS2D_manager + type(ModelVarManager), intent(in), target, optional:: cp_AUXVARS2D_manager !---------------------------------------------- - this%ptr_MP_AUXVARS2D_manager => mp_AUXVARS2D_manager - + if ( present(mp_AUXVARS2D_manager) ) then + this%ptr_MP_AUXVARS2D_manager => mp_AUXVARS2D_manager + end if + if ( present(cp_AUXVARS2D_manager) ) then + this%ptr_CP_AUXVARS2D_manager => cp_AUXVARS2D_manager + end if return end subroutine AtmosVars_Regist_physvar_manager @@ -950,7 +958,7 @@ subroutine AtmosVars_CalcDiagvar2D( this, field_name, field_work ) do n=1, field_work%mesh%LOCAL_MESH_NUM lcmesh2D => field_work%mesh%lcmesh_list(n) call vars_calc_diagnoseVar2D_lc( field_name, field_work%local(n)%val, & - this%ptr_MP_AUXVARS2D_manager, & + this%ptr_MP_AUXVARS2D_manager, this%ptr_CP_AUXVARS2D_manager, & field_work%mesh, lcmesh2D, lcmesh2D%refElem2D ) end do !$acc wait(1) @@ -962,48 +970,107 @@ end subroutine AtmosVars_CalcDiagvar2D !OCL SERIAL subroutine vars_calc_diagnoseVar2D_lc( field_name, & ! (in) var_out, & ! (out) - MP_auxvars2D, mesh2D, lcmesh, elem ) ! (in) + MP_auxvars2D, CP_auxvars2D, mesh2D, lcmesh, elem ) ! (in) use mod_atmos_phy_mp_vars, only: & AtmosPhyMpVars_GetLocalMeshFields_sfcflx + use mod_atmos_phy_cp_vars, only: & + AtmosPhyCpVars_GetLocalMeshFields_sfcflx implicit none class(LocalMesh2D), intent(in) :: lcmesh class(ElementBase2D), intent(in) :: elem character(*), intent(in) :: field_name real(RP), intent(out) :: var_out(elem%Np,lcmesh%NeA) - class(ModelVarManager), intent(inout) :: MP_auxvars2D + type(ModelVarManager), intent(inout), pointer :: MP_auxvars2D + type(ModelVarManager), intent(inout), pointer :: CP_auxvars2D class(MeshBase2D), intent(in) :: mesh2D integer :: ke, p + logical :: sw_MP, sw_CP class(LocalMeshFieldBase), pointer :: SFLX_rain_MP, SFLX_snow_MP, SFLX_ENGI_MP + class(LocalMeshFieldBase), pointer :: SFLX_rain_CP, SFLX_snow_CP, SFLX_ENGI_CP !------------------------------------------------------------------------- + sw_MP = associated(MP_auxvars2D) + sw_CP = associated(CP_auxvars2D) + select case(trim(field_name)) case('RAIN', 'SNOW') - call AtmosPhyMpVars_GetLocalMeshFields_sfcflx( & - lcmesh%lcdomID, mesh2D, MP_auxvars2D, & - SFLX_rain_MP, SFLX_snow_MP, SFLX_ENGI_MP ) + if ( sw_MP ) then + call AtmosPhyMpVars_GetLocalMeshFields_sfcflx( & + lcmesh%lcdomID, mesh2D, MP_auxvars2D, & + SFLX_rain_MP, SFLX_snow_MP, SFLX_ENGI_MP ) + end if + if ( sw_CP ) then + call AtmosPhyCpVars_GetLocalMeshFields_sfcflx( & + lcmesh%lcdomID, mesh2D, CP_auxvars2D, & + SFLX_rain_CP, SFLX_snow_CP, SFLX_ENGI_CP ) + end if end select select case(trim(field_name)) case('RAIN') - !$omp parallel do - !$acc parallel loop collapse(2) present(SFLX_rain_MP%val, var_out) async(1) + !$omp parallel + !$acc parallel present(SFLX_rain_MP%val, var_out) async(1) + !$omp do + !$acc loop collapse(2) do ke=lcmesh%NeS, lcmesh%NeE do p=1, elem%Np - var_out(p,ke) = SFLX_rain_MP%val(p,ke) + var_out(p,ke) = 0.0_RP end do end do + if ( sw_MP ) then + !$omp do + !$acc loop collapse(2) + do ke=lcmesh%NeS, lcmesh%NeE + do p=1, elem%Np + var_out(p,ke) = var_out(p,ke) + SFLX_rain_MP%val(p,ke) + end do + end do + end if + if ( sw_CP ) then + !$omp do + !$acc loop collapse(2) + do ke=lcmesh%NeS, lcmesh%NeE + do p=1, elem%Np + var_out(p,ke) = var_out(p,ke) + SFLX_rain_CP%val(p,ke) + end do + end do + end if + !$omp end parallel + !$acc end parallel case('SNOW') - !$omp parallel do - !$acc parallel loop collapse(2) present(var_out, SFLX_snow_MP) async(1) + !$omp parallel + !$acc parallel present(SFLX_snow_MP%val, var_out) async(1) + !$omp do + !$acc loop collapse(2) do ke=lcmesh%NeS, lcmesh%NeE do p=1, elem%Np - var_out(p,ke) = SFLX_snow_MP%val(p,ke) + var_out(p,ke) = 0.0_RP end do end do + if ( sw_MP ) then + !$omp do + !$acc loop collapse(2) + do ke=lcmesh%NeS, lcmesh%NeE + do p=1, elem%Np + var_out(p,ke) = var_out(p,ke) + SFLX_snow_MP%val(p,ke) + end do + end do + end if + if ( sw_CP ) then + !$omp do + !$acc loop collapse(2) + do ke=lcmesh%NeS, lcmesh%NeE + do p=1, elem%Np + var_out(p,ke) = var_out(p,ke) + SFLX_snow_CP%val(p,ke) + end do + end do + end if + !$omp end parallel + !$acc end parallel case default LOG_ERROR("AtmosVars_calc_diagnoseVar2D_lc",*) 'The name of diagnostic variable is not suported. Check!', field_name call PRC_abort From 5bf0a86443332bee4379230f356d0b006057ee68 Mon Sep 17 00:00:00 2001 From: ykawai Date: Wed, 5 Aug 2026 16:35:06 +0900 Subject: [PATCH 4/6] Improve the treatmensts of no energy root issues in the MCA scheme with available-water constraints. --- .../scale_atm_phy_cp_dgm_mconv_adjustment.F90 | 1314 ++++++++++++----- 1 file changed, 925 insertions(+), 389 deletions(-) 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 index 1901b3b2..59d8762b 100644 --- a/FElib/src/cumulus/scale_atm_phy_cp_dgm_mconv_adjustment.F90 +++ b/FElib/src/cumulus/scale_atm_phy_cp_dgm_mconv_adjustment.F90 @@ -1,12 +1,20 @@ !> module FElib / Atmosphere / Physics cumulus parameterization !! !! @par Description -!! A module to provide a moist convective adjustment scheme for cumulus parameterization in atmospheric model -!! +!! A module providing a moist convective adjustment scheme. +!! +!! Saturated convective instability is diagnosed from the vertical gradient of saturated moist static energy in nearly saturated layers. +!! The diagnosed layer is adjusted toward a moist-neutral temperature profile subject to column-integrated moist-static-energy and available-water constraints. +!! Condensed water is immediately removed as precipitation. +!! !! @author Yuta Kawai, Team SCALE !! !! @par Reference +!! Manabe, S., J. Smagorinsky, and R. F. Strickler (1965): +!! Simulated Climatology of a General Circulation Model with a +!! Hydrologic Cycle. Monthly Weather Review, 93, 769-798. !! +!! !------------------------------------------------------------------------------- #include "scaleFElib.h" module scale_atm_phy_cp_dgm_mconv_adjustment @@ -24,6 +32,7 @@ module scale_atm_phy_cp_dgm_mconv_adjustment CPdry => CONST_CPDRY, & CVdry => CONST_CVdry, & Rdry => CONST_Rdry, & + Rvap => CONST_Rvap, & PRES0 => CONST_PRE00, & EPS => CONST_EPS use scale_atmos_hydrometeor, only: & @@ -61,32 +70,60 @@ module scale_atm_phy_cp_dgm_mconv_adjustment !++ 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 + + !* MCA paramters + + ! Iteration limits + integer, parameter :: MCA_MAX_ADJUST_ITER = 200 + integer, parameter :: MCA_MAX_ENERGY_ITER = 50 + integer, parameter :: MCA_MAX_WATER_ITER = 60 + + ! Column-energy root solve + real(RP), parameter :: MCA_COLUMN_ENERGY_RTOL = 5.0E-8_RP + real(RP), parameter :: MCA_COLUMN_ENERGY_ATOL = 1.0_RP + real(RP), parameter :: MCA_BOUNDARY_ENERGY_RTOL = 1.0E-7_RP + + ! Local saturated-MSE solve for moist-adiabat integration + real(RP), parameter :: MCA_LOCAL_HMSE_RTOL = 1.0E-11_RP + real(RP), parameter :: MCA_LOCAL_HMSE_ATOL = 1.0E-8_RP + real(RP), parameter :: MCA_TEMP_TOL = 1.0E-8_RP + + ! Water-feasibility constraint + real(RP), parameter :: MCA_WATER_RTOL = 1.0E-12_RP + real(RP), parameter :: MCA_WATER_ATOL = 1.0E-12_RP + + ! Saturation and instability thresholds + real(RP), parameter :: MCA_RH_FORCED_SATURATION = 0.999_RP + real(RP), parameter :: MCA_RH_TRIGGER = 0.999_RP + + ! Thresholds for suppressing numerical-noise detection + real(RP), parameter :: MCA_HMSE_GRAD_TOL = 1.0E-3_RP + real(RP), parameter :: MCA_HMSE_DIFF_TOL = 1.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) + real(RP), parameter :: MCA_Z_TOL = 100.0_RP * epsilon(1.0_RP) + + ! Moist-neutral-profile solver status + integer, parameter :: MCA_PROFILE_SUCCESS = 0 + integer, parameter :: MCA_PROFILE_NO_ENERGY_ROOT = 1 + integer, parameter :: MCA_PROFILE_ADIABAT_FAILURE = 2 + integer, parameter :: MCA_PROFILE_ENERGY_MAXITER = 3 + integer, parameter :: MCA_PROFILE_NO_ADJUSTED_NODE = 4 + integer, parameter :: MCA_PROFILE_NO_WATER_FEASIBLE_STATE = 5 contains + !> Setup a module for moist convective adjustment scheme subroutine atm_phy_cp_dgm_mconv_adjustment_setup() implicit none !---------------------------------------------------- return end subroutine atm_phy_cp_dgm_mconv_adjustment_setup + !> Calculate tendencies with moist convective adjustment scheme !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 ) + DENS_t, RHOT_t, RHOQV_t, SFLX_RAIN, & ! (out) + DDENS, DRHOT, QV, PT, PRES, DENS_hyd, Rtot, CPtot, & ! (in) + dtsec, lmesh, elem, elem1D ) ! (in) implicit none class(LocalMesh3D), intent(in) :: lmesh class(ElementBase3D), intent(in) :: elem @@ -113,8 +150,9 @@ subroutine atm_phy_cp_dgm_mconv_adjustment_calc_tendency( & 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) :: qsat + real(RP) :: rh_z + real(RP) :: rtot_, cptot_, qdry real(RP) :: elem_width_z real(RP) :: int_weight(elem%Nnode_v,lmesh%NeZ) @@ -124,59 +162,68 @@ subroutine atm_phy_cp_dgm_mconv_adjustment_calc_tendency( & 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_expand + integer :: MCA_MAX_MASK_EXPAND + integer :: lbase, ltop, lbase_try, ltop_try + real(RP) :: zbase_try, ztop_try + logical :: mask_expanded + logical :: is_unstable + logical :: unstable_core_mask(elem%Nnode_v,lmesh%NeZ) + logical :: forced_saturation_mask(elem%Nnode_v,lmesh%NeZ) + logical :: adjustment_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) + real(RP) :: temp_adj(elem%Nnode_v,lmesh%NeZ) + real(RP) :: qvap_adj(elem%Nnode_v,lmesh%NeZ) + real(RP) :: rhoqvap_adj(elem%Nnode_v,lmesh%NeZ) + real(RP) :: rhoprecip_adj(elem%Nnode_v,lmesh%NeZ) ! Water-vapor mass density diagnosed to be removed as precipitation during the current adjustment iteration [kg m-3] + real(RP) :: rhoprecip_accum(elem%Nnode_v,lmesh%NeZ) ! Accumulated vapor mass density removed as precipitation [kg m-3] + real(RP) :: rhoqdry_adj(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 :: profile_converged + integer :: profile_status logical :: do_adjustment logical :: adjustment_converged integer :: Hslice_b(elem%Nnode_h1D**2), Hslice_t(elem%Nnode_h1D**2) - real(RP) :: Rvap + logical :: debug_flag - integer, parameter :: MAX_ITER_ADJ = 10 + integer :: lactive_base + integer :: lactive_top + logical :: active_region_initialized !---------------------------------------------------- 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) ) + allocate( zlev_diag(nlev_diag) ) Hslice_b(:) = elem%Hslice(:,1) Hslice_t(:) = elem%Hslice(:,elem%Nnode_v) - Rvap = CP_VAPOR - CV_VAPOR + MCA_MAX_MASK_EXPAND = nlev_diag !$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 ) + !$omp dens_z, pres_z, temp_z, zlev_z, pott_z, qvap_z, qsat, rh_z, rtot_,cptot_, & + !$omp dens_ini, pott_ini, qvap_ini, rhoqvap_ini, & + !$omp zlev_diag, & + !$omp lbase, ltop, lbase_try, ltop_try, zbase_try, ztop_try, iter_expand, mask_expanded, & + !$omp is_unstable, unstable_core_mask, forced_saturation_mask, adjustment_mask, & + !$omp temp_adj, qdry, qvap_adj, rhoqdry_adj, rhoqvap_adj, rhoprecip_adj, rhoprecip_accum, & + !$omp profile_converged, profile_status, do_adjustment, adjustment_converged, & + !$omp int_weight,elem_width_z, debug_flag, lactive_base, lactive_top, active_region_initialized ) do ke_xy=1, lmesh%Ne2D do ph=1, elem%Nnode_h1D**2 + ! if ( lmesh%PRC_myrank == 1 .and. ph == 8 .and. ke_xy == 8 ) then + ! debug_flag = .true. + ! else + debug_flag = .false. + ! end if + !- Extract vertical 1D DG column do ke_z=1, lmesh%NeZ @@ -204,89 +251,181 @@ subroutine atm_phy_cp_dgm_mconv_adjustment_calc_tendency( & qvap_ini(:,ke_z) = qvap_z(:,ke_z) rhoqvap_ini(:,ke_z) = dens_ini(:,ke_z) * qvap_ini(:,ke_z) + temp_adj(:,ke_z) = temp_z(:,ke_z) + rhoqvap_adj(:,ke_z) = rhoqvap_ini(:,ke_z) rhoprecip_accum(:,ke_z) = 0.0_RP end do - !** Loop for iteration of the moist convective adjustment + !* Iteratively diagnose and remove remaining saturated convective instability. *************************************** + ! A single profile adjustment may generate a new unstable layer adjacent to or outside the previously adjusted region. 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 + active_region_initialized = .false. + lactive_base = 0 + lactive_top = 0 + forced_saturation_mask(:,:) = .false. - 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 + do iter_adj=1, MCA_MAX_ADJUST_ITER + + qvap_adj(:,:) = rhoqvap_adj(:,:) / dens_z(:,:) + + ! Reconstruct the saturation constraint from the current state. + forced_saturation_mask(:,:) = .false. !- 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) + call diagnose_convective_layer( lbase, ltop, is_unstable, zlev_diag, & ! (out) + temp_adj, qvap_adj, pres_z, zlev_z, lmesh, elem, nlev_diag, debug_flag ) ! (in) - ! write(*,*) "-- ph=", ph, "ke_xy=", ke_xy, "iter_adj=", iter_adj - ! write(*,*) " lbase, ltop, is_unstable = ", lbase, ltop, is_unstable + if ( debug_flag ) then + write(*,*) "ke_xy=", ke_xy, "ph=", ph, "iter_adj=", iter_adj + write(*,*) " lbase, ltop, is_unstable = ", lbase, ltop, is_unstable + end if 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 + ! Maintain the smallest diagnostic-level interval containing all unstable or expanded regions encountered during this column adjustment. + if ( .not. active_region_initialized ) then + lactive_base = lbase + lactive_top = ltop + active_region_initialized = .true. + else if ( lbase <= lactive_top + 1 .and. & + ltop >= lactive_base - 1 ) then + ! The new unstable layer overlaps or directly touches the previously active connected layer. + lactive_base = min(lactive_base, lbase) + lactive_top = max(lactive_top, ltop) + else + ! Union of the previously adjusted region and the newly diagnosed unstable region. + lactive_base = min(lactive_base, lbase) + lactive_top = max(lactive_top, ltop) end if - - !- Immediate precipitation removal - rhoprecip_accum(:,:) = rhoprecip_accum(:,:) + rhoprecip_aj(:,:) + if ( debug_flag ) then + write(*,*) " persistent active region=", lactive_base, lactive_top + end if + + ! Mark nearly saturated DG nodes inside the diagnosed unstable core. + ! These nodes will be constrained to saturation in each trial moist-neutral profile. + call make_dg_vertical_range_mask( & + unstable_core_mask, & ! (out) + zlev_z, zlev_diag(lbase), zlev_diag(ltop), & ! (in) + elem%Nnode_v, lmesh%NeZ ) ! (in) + + do ke_z = 1, lmesh%NeZ + do pz = 1, elem%Nnode_v + if ( unstable_core_mask(pz,ke_z) ) then + call ATMOS_SATURATION_pres2qsat_liq( temp_adj(pz,ke_z), pres_z(pz,ke_z), & + qsat ) ! (out) + + rh_z = qvap_adj(pz,ke_z) / max(qsat, EPS) + forced_saturation_mask(pz,ke_z) = rh_z >= MCA_RH_FORCED_SATURATION + end if + end do + end do - rhoqdry_aj(:,:) = dens_z(:,:) - rhoqvap_aj(:,:) - rhoprecip_aj(:,:) - dens_z(:,:) = dens_z(:,:) - rhoprecip_aj(:,:) + !- Construct a moist-neutral profile. + ! If no water-feasible and energy-conserving solution exists, expand the adjustment layer and retry. - qvap_aj(:,:) = rhoqvap_aj(:,:) / dens_z(:,:) - qdry_aj(:,:) = rhoqdry_aj(:,:) / dens_z(:,:) + lbase_try = lactive_base + ltop_try = lactive_top - !- Use the post-precipitation state in the next iteration + profile_converged = .false. + profile_status = MCA_PROFILE_ENERGY_MAXITER - 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) + do iter_expand=0, MCA_MAX_MASK_EXPAND + zbase_try = zlev_diag(lbase_try) + ztop_try = zlev_diag(ltop_try) + + !- Construct a mask for the adjustment layer + call make_dg_vertical_range_mask( adjustment_mask, & ! (out) + zlev_z, zbase_try, ztop_try, elem%Nnode_v, lmesh%NeZ ) ! (in) + + if ( debug_flag ) then + write(*,*) "* Adjustment-profile attempt: iter_expand=", iter_expand + write(*,*) " lactive_base, lactive_top=", lactive_base, lactive_top + write(*,*) " lbase_try, ltop_try=", lbase_try, ltop_try + write(*,*) " zbase_try, ztop_try=", zbase_try, ztop_try + write(*,*) " saturation, mixing nodes=", count(forced_saturation_mask), count(adjustment_mask) + end if + + !- Solve for a water-feasible and energy-conserving saturated profile in the adjustment layer + + call solve_moist_neutral_adjustment( & + temp_adj, qvap_adj, rhoqvap_adj, rhoprecip_adj, profile_converged, profile_status, & ! (out) + qvap_z, pres_z, zlev_z, dens_z, temp_z, & ! (in) + int_weight, adjustment_mask, forced_saturation_mask, elem%Nnode_v, lmesh%NeZ, debug_flag ) ! (in) + + if ( profile_converged ) then + ! The region added by water-feasibility expansion also becomes part of the persistent active region. + lactive_base = min(lactive_base, lbase_try) + lactive_top = max(lactive_top, ltop_try) + + if ( debug_flag ) then + write(*,*) "Moist-neutral profile converged: iter_expand=", iter_expand, "final: lbase, ltop=", lbase_try, ltop_try, ", zbase, ztop=", zbase_try, ztop_try + + ! write(*,*) " check:" + ! call diagnose_convective_layer( lbase, ltop, is_unstable, zlev_diag, & ! (out) + ! temp_aj, qvap_aj, pres_z, zlev_z, lmesh, elem, nlev_diag, debug_flag, .true., lbase_try, ltop_try ) ! (in) + end if + exit + end if + + select case ( profile_status ) + case ( MCA_PROFILE_NO_ENERGY_ROOT ) + call expand_adjustment_layer( lbase_try, ltop_try, mask_expanded, nlev_diag ) + if ( .not. mask_expanded ) then + exit + end if + case default + exit + end select + + end do ! end loop for iter_expand + + if ( .not. profile_converged ) then + LOG_INFO("atm_phy_cp_dgm_mconv_adjustment_calc_tendency",*) "iter_adj=", iter_adj, "ph=", ph, "ke_xy=", ke_xy, "iter_expand=", iter_expand + select case ( profile_status ) + case ( MCA_PROFILE_NO_ENERGY_ROOT ) + LOG_INFO("atm_phy_cp_dgm_mconv_adjustment_calc_tendency",*) "No energy-conserving root exists in the water-feasible temperature interval, even after adjustment-layer expansion." + LOG_INFO("atm_phy_cp_dgm_mconv_adjustment_calc_tendency",*) "original lbase,ltop=", lbase, ltop, "final lbase,ltop=", lbase_try, ltop_try + do_adjustment = .false. + exit + ! call PRC_abort + case ( MCA_PROFILE_NO_WATER_FEASIBLE_STATE ) + LOG_INFO("atm_phy_cp_dgm_mconv_adjustment_calc_tendency",*) "No water-feasible trial moist-neutral profile exists." + call PRC_abort + case ( MCA_PROFILE_ADIABAT_FAILURE ) + LOG_INFO("atm_phy_cp_dgm_mconv_adjustment_calc_tendency",*) "Moist-adiabat profile integration failed." + call PRC_abort + case ( MCA_PROFILE_ENERGY_MAXITER ) + LOG_INFO("atm_phy_cp_dgm_mconv_adjustment_calc_tendency",*) "Energy root was bracketed, but bisection did not converge." + call PRC_abort + case ( MCA_PROFILE_NO_ADJUSTED_NODE ) + LOG_INFO("atm_phy_cp_dgm_mconv_adjustment_calc_tendency",*) "No adjusted DG node was found." + call PRC_abort + case default + LOG_INFO("atm_phy_cp_dgm_mconv_adjustment_calc_tendency",*) "Unknown moist-neutral-profile construction error." + call PRC_abort + end select + end if + + ! Accumulate the precipitation density from the current iteration + rhoprecip_accum(:,:) = rhoprecip_accum(:,:) + rhoprecip_adj(:,:) - 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 + ! Update the state for the next iteration + temp_z(:,:) = temp_adj(:,:) + qvap_z(:,:) = qvap_adj(:,:) + + 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 + if ( do_adjustment .and. ( .not. adjustment_converged ) ) then + LOG_INFO("atm_phy_cp_dgm_mconv_adjustment_calc_tendency",*) "Moist convective adjustment did not converge: ph=", ph, "ke_xy=", ke_xy call PRC_abort end if @@ -296,6 +435,19 @@ subroutine atm_phy_cp_dgm_mconv_adjustment_calc_tendency( & if ( do_adjustment ) then do ke_z=1, lmesh%NeZ ke = ke_xy + (ke_z-1)*lmesh%Ne2D + + dens_z(:,ke_z) = dens_ini(:,ke_z) - rhoprecip_accum(:,ke_z) + qvap_z(:,ke_z) = rhoqvap_adj(:,ke_z) / dens_z(:,ke_z) + + do pz=1, elem%Nnode_v + qdry = 1.0_RP - qvap_z(pz,ke_z) + rtot_ = Rdry * qdry + Rvap * qvap_z(pz,ke_z) + cptot_ = CPdry * qdry + CP_VAPOR * qvap_z(pz,ke_z) + + pres_z(pz,ke_z) = dens_z(pz,ke_z) * rtot_ * temp_z(pz,ke_z) + pott_z(pz,ke_z) = temp_adj(pz,ke_z) * ( PRES0 / pres_z(pz,ke_z) )**( rtot_ / cptot_ ) + end do + 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 @@ -303,6 +455,9 @@ subroutine atm_phy_cp_dgm_mconv_adjustment_calc_tendency( & RHOQV_t(p,ke) = ( dens_z(pz,ke_z) * qvap_z(pz,ke_z) - rhoqvap_ini(pz,ke_z) ) / dtsec end do end do + + SFLX_RAIN(ph,ke_xy) = sum( rhoprecip_accum(:,:) * int_weight(:,:) ) / dtsec + else do ke_z=1, lmesh%NeZ ke = ke_xy + (ke_z-1)*lmesh%Ne2D @@ -313,16 +468,16 @@ subroutine atm_phy_cp_dgm_mconv_adjustment_calc_tendency( & RHOQV_t(p,ke) = 0.0_RP end do end do + SFLX_RAIN(ph,ke_xy) = 0.0_RP 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 + !> Finalize a module for moist convective adjustment scheme subroutine atm_phy_cp_dgm_mconv_adjustment_finalize() implicit none !---------------------------------------------------- @@ -330,23 +485,21 @@ subroutine atm_phy_cp_dgm_mconv_adjustment_finalize() end subroutine atm_phy_cp_dgm_mconv_adjustment_finalize -!- private subroutines +!- 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 ) + subroutine build_diagnostic_profile_from_dg( & + hmse_diag, rh_diag, zlev_diag, & ! (out) + hmse_z, rh_z, zlev_z, lmesh, elem, nlev ) ! (in) 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) :: hmse_diag(nlev) + real(RP), intent(out) :: rh_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) :: hmse_z(elem%Nnode_v,lmesh%NeZ) + real(RP), intent(in) :: rh_z(elem%Nnode_v,lmesh%NeZ) real(RP), intent(in) :: zlev_z(elem%Nnode_v,lmesh%NeZ) integer :: pz, ke_z @@ -357,86 +510,104 @@ subroutine build_diag_zprofile_from_dg( & 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) + hmse_diag(l) = hmse_z(1,1) + rh_diag(l) = rh_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) + hmse_diag(l) = hmse_z(pz,ke_z) + rh_diag(l) = rh_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) ) + hmse_diag(l) = 0.5_RP * ( hmse_z(Nnode_v,ke_z) + hmse_z(1,ke_z+1) ) + rh_diag(l) = 0.5_RP * ( rh_z(Nnode_v,ke_z) + rh_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) + hmse_diag(l) = hmse_z(Nnode_v,ke_z) + rh_diag(l) = rh_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 + end subroutine build_diagnostic_profile_from_dg !OCL SERIAL - subroutine diagnose_convective_layer( lbase, ltop, is_unstable, & - temp, pott, qvap, pres, zlev, nlev ) + subroutine diagnose_convective_layer( lbase, ltop, is_unstable, zlev_diag, & ! (out) + temp_z, qvap_z, pres_z, zlev_z, lmesh, elem, nlev, debug_flag ) ! (in) implicit none + class(LocalMesh3D), intent(in) :: lmesh + class(ElementBase3D), intent(in) :: elem + integer, intent(in) :: nlev 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) + real(RP), intent(out) :: zlev_diag(nlev) + real(RP), intent(in) :: temp_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) + logical, intent(in) :: debug_flag + integer :: pz, ke_z integer :: l - real(RP) :: qsat(nlev) + + real(RP) :: hmse_sat_z(elem%Nnode_v,lmesh%NeZ) + real(RP) :: rh_z(elem%Nnode_v,lmesh%NeZ) + real(RP) :: qsat, cptot + real(RP) :: rh(nlev) real(RP) :: hmse_sat(nlev) - real(RP) :: cptot real(RP) :: dz - real(RP) :: dhmse_sat_dz + real(RP) :: dhmse_sat, dhmse_sat_dz real(RP) :: rh_lyr logical :: unstable_pair logical :: is_inside + + integer :: Nnode_v !------------------------------------ - call ATMOS_SATURATION_pres2qsat_liq( & - nlev, 1, nlev, & ! (in) - temp, pres, & ! (in) - qsat ) ! (out) + Nnode_v = elem%Nnode_v + + do ke_z=1, lmesh%NeZ + do pz=1, Nnode_v - do l=1, nlev - rh(l) = qvap(l) / max(qsat(l), EPS) + call ATMOS_SATURATION_pres2qsat_liq( temp_z(pz,ke_z), pres_z(pz,ke_z), & + qsat ) ! (out) - cptot = CPdry * (1.0_RP - qsat(l)) + CP_VAPOR * qsat(l) - hmse_sat(l) = cptot * temp(l) + GRAV * zlev(l) + LHV0 * qsat(l) + rh_z(pz,ke_z) = qvap_z(pz,ke_z) / max(qsat, EPS) + + cptot = CPdry * (1.0_RP - qsat) + CP_VAPOR * qsat + hmse_sat_z(pz,ke_z) = cptot * temp_z(pz,ke_z) + GRAV * zlev_z(pz,ke_z) + LHV0 * qsat end do - + end do + + call build_diagnostic_profile_from_dg( hmse_sat, rh, zlev_diag, & + hmse_sat_z, rh_z, zlev_z, & + lmesh, elem, nlev ) + + !- Diagnose unstable layers based on the saturated MSE profile and relative humidity + 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 + dz = zlev_diag(l+1) - zlev_diag(l) + dhmse_sat = hmse_sat(l+1) - hmse_sat(l) + dhmse_sat_dz = dhmse_sat / 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 ) - + unstable_pair = & + ( min(rh(l), rh(l+1)) >= MCA_RH_TRIGGER ) & + .and. ( dhmse_sat < -MCA_HMSE_DIFF_TOL ) & + .and. ( dhmse_sat_dz < -MCA_HMSE_GRAD_TOL ) + + if ( unstable_pair ) then if ( .not. is_inside ) then lbase = l @@ -449,18 +620,27 @@ subroutine diagnose_convective_layer( lbase, ltop, is_unstable, & end do if ( is_inside ) then - if ( zlev(ltop) - zlev(lbase) >= MCA_MIN_UNSTABLE_DEPTH ) is_unstable = .true. + if ( zlev_diag(ltop) - zlev_diag(lbase) >= MCA_MIN_UNSTABLE_DEPTH ) is_unstable = .true. + end if + + if ( debug_flag .and. is_unstable ) then + write(*,*) "Diagnose convective layer:" + write(*,*) "hmse_sat=", hmse_sat(lbase:ltop) + write(*,*) "rh=", rh(lbase:ltop) + write(*,*) "hmse_sat_dg=", hmse_sat_z(:,:) + write(*,*) "rh_dg=", rh_z(:,:) end if + return end subroutine diagnose_convective_layer !OCL SERIAL - subroutine make_dg_adjustment_mask( adjust_mask, & ! (out) - zlev, zbase, ztop, npz, nez ) ! (in) + subroutine make_dg_vertical_range_mask( adjustment_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) + logical, intent(out) :: adjustment_mask(npz,nez) real(RP), intent(in) :: zlev(npz,nez) real(RP), intent(in) :: zbase real(RP), intent(in) :: ztop @@ -468,91 +648,161 @@ subroutine make_dg_adjustment_mask( adjust_mask, & ! (out) integer :: ke_z, pz !---------------------------------------------------- - adjust_mask(:,:) = .false. + adjustment_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. + adjustment_mask(pz,ke_z) = .true. end if end do - end do + end do return - end subroutine make_dg_adjustment_mask + end subroutine make_dg_vertical_range_mask + !> Expand the adjustment layer by one diagnostic level above and below, if possible + !! Numerical fallback: enlarge the adjustment region by one diagnostic level on each available side. This is not a physically based entrainment closure. +!OCL SERIAL + subroutine expand_adjustment_layer( & + lbase, ltop, expanded, & ! (inout,out) + nlev ) ! (in) + implicit none + + integer, intent(inout) :: lbase + integer, intent(inout) :: ltop + logical, intent(out) :: expanded + integer, intent(in) :: nlev + + logical :: can_expand_below + logical :: can_expand_above + !------------------------------------------------------------ + + can_expand_below = lbase > 1 + can_expand_above = ltop < nlev + + expanded = .false. + + ! Symmetrically expand by one diagnostic level whenever possible. + if ( can_expand_below ) then + lbase = lbase - 1 + expanded = .true. + end if + + if ( can_expand_above ) then + ltop = ltop + 1 + expanded = .true. + end if + return + end subroutine expand_adjustment_layer + + !> Solve for a trial moist-neutral profile satisfying: 1. available-water constraint, 2. column-integrated moist-static-energy conservation. + !! + !! The pressure and density profiles are held fixed during the solve. + !! Nodes selected by forced_saturation_mask are set to saturation. + !! At other adjusted nodes, the original vapor mixing ratio is retained unless it exceeds saturation. + !! !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 ) + subroutine solve_moist_neutral_adjustment( & + temp_adj, qvap_adj, rhoqvap_aj, rhoprecip_aj, converged, profile_status, & ! (out) + qvap, pres, zlev, dens, temp, & ! (in) + int_weight, mix_mask, saturation_mask, npz, nez, debug_flag ) ! (in) 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) :: temp_adj(npz,nez) + real(RP), intent(out) :: qvap_adj(npz,nez) real(RP), intent(out) :: rhoqvap_aj(npz,nez) real(RP), intent(out) :: rhoprecip_aj(npz,nez) logical, intent(out) :: converged + integer, intent(out) :: profile_status 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) + logical, intent(in) :: mix_mask(npz,nez) + logical, intent(in) :: saturation_mask(npz,nez) + logical, intent(in) :: debug_flag + ! Initial column-integrated quantities real(RP) :: energy_target - real(RP) :: energy_aj + real(RP) :: water_mass_ini + + ! Trial-profile quantities + real(RP) :: energy_trial + real(RP) :: water_trial real(RP) :: residual - real(RP) :: tbase_lo - real(RP) :: tbase_hi - real(RP) :: tbase_mid + ! Energy root bracket + real(RP) :: tbase_lo, tbase_mid, tbase_hi + real(RP) :: residual_lo, residual_hi + real(RP) :: energy_lo, energy_hi + real(RP) :: water_lo, water_hi - integer :: ke_z, pz - integer :: iter + ! Water-feasible upper-bound search + real(RP) :: twater_lo, twater_mid, twater_hi + real(RP) :: tbase_water_max + real(RP) :: water_mid + real(RP) :: energy_dummy - real(RP) :: qsat_aj(npz,nez) + ! Work arrays + real(RP) :: temp_trial(npz,nez) + real(RP) :: qvap_trial(npz,nez) - real(RP) :: water_mass_ini - real(RP) :: water_mass_aj - real(RP) :: precip_mass_total + ! Precipitation diagnostics + real(RP) :: water_mass_adj real(RP) :: precip_mass_local + real(RP) :: precip_mass_total real(RP) :: positive_cond_mass real(RP) :: precip_scale + ! Tolerances + real(RP) :: energy_tol + real(RP) :: boundary_energy_tol + real(RP) :: water_tol + + integer :: ke_z, pz + integer :: iter + integer :: iter_water + logical :: found_base - real(RP), parameter :: MCA_TBASE_RANGE = 40.0_RP + logical :: adiabat_profile_ok - logical :: is_converged_ref_profile - logical :: water_feasible - logical :: residual_evaluated - !---------------------------------------------------- + logical :: water_lo_feasible + logical :: water_hi_feasible + + real(RP), parameter :: MCA_TBASE_RANGE = 40.0_RP + !------------------------------------------------------------ converged = .false. + profile_status = MCA_PROFILE_ENERGY_MAXITER ! Initialize the adjusted profile with the original state - temp_aj(:,:) = temp(:,:) - qvap_aj(:,:) = qvap(:,:) + temp_adj(:,:) = temp(:,:) + qvap_adj(:,:) = qvap(:,:) rhoqvap_aj(:,:) = dens(:,:) * qvap(:,:) rhoprecip_aj(:,:) = 0.0_RP !- Calculate initial column-integrated MSE - call column_mse( energy_target, & ! (out) + call integ_masked_column_mse( energy_target, & ! (out) temp, qvap, zlev, dens, & ! (in) - int_weight, adjust_mask, npz, nez ) ! (in) + int_weight, mix_mask, npz, nez ) ! (in) + + energy_tol = MCA_COLUMN_ENERGY_ATOL + MCA_COLUMN_ENERGY_RTOL * abs(energy_target) + boundary_energy_tol = MCA_COLUMN_ENERGY_ATOL + MCA_BOUNDARY_ENERGY_RTOL * abs(energy_target) !- 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 + if ( mix_mask(pz,ke_z) ) then tbase_mid = temp(pz,ke_z) found_base = .true. exit @@ -560,110 +810,325 @@ subroutine construct_moist_neutral_profile( & end do if ( found_base ) exit end do + if ( .not. found_base ) then - write(*,*) "Error: No adjusted nodes found in the column. Aborting adjustment." + profile_status = MCA_PROFILE_NO_ADJUSTED_NODE + if ( debug_flag ) then + write(*,*) "construct_moist_neutral_profile: No adjusted DG node was found." + end if 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 + if ( mix_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 + water_tol = MCA_WATER_ATOL + MCA_WATER_RTOL * max(abs(water_mass_ini),1.0_RP) - tbase_lo = tbase_mid - MCA_TBASE_RANGE - tbase_hi = tbase_mid + MCA_TBASE_RANGE + !- Step 1: + ! Determine the upper limit of the water-feasible base temperature. + ! A trial moist-neutral profile is water-feasible when its vapor mass does not exceed the initial vapor mass in the adjustment region. - water_feasible = .false. - residual_evaluated = .false. + twater_lo = tbase_mid - MCA_TBASE_RANGE + twater_hi = tbase_mid + MCA_TBASE_RANGE - !- - do iter=1, MCA_MAX_ENERGY_ITER + ! Evaluate the water mass at the lower bound of the temperature range - tbase_mid = 0.5_RP * ( tbase_lo + tbase_hi ) + call evaluate_trial_moist_neutral_profile( temp_trial, qvap_trial, water_lo, energy_dummy, adiabat_profile_ok, & ! (out) + twater_lo, temp, qvap, dens, pres, zlev, int_weight, mix_mask, saturation_mask, nez, npz, debug_flag ) ! (in) - ! Evaluate the moist adiabat temperature profile for the current base temperature + if ( adiabat_profile_ok ) then + water_lo_feasible = ( water_lo <= water_mass_ini + water_tol ) + else + profile_status = MCA_PROFILE_ADIABAT_FAILURE + if ( debug_flag ) then + write(*,*) "construct_moist_neutral_profile: Moist-adiabat integration failed. tbase_lo=", twater_lo + end if + return + end if - 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. water_lo_feasible ) then + profile_status = MCA_PROFILE_NO_WATER_FEASIBLE_STATE - if ( .not. is_converged_ref_profile ) then - converged = .false. - exit + if ( debug_flag ) then + write(*,*) "No water-feasible saturated profile exists." + write(*,*) " tbase_lo=", twater_lo, "water_mass_lo=", water_lo, "water_mass_ini =", water_mass_ini end if + return + end if + + ! Evaluate the water mass at the upper bound of the temperature range - ! Set the adjusted vapor to stauration + call evaluate_trial_moist_neutral_profile( temp_trial, qvap_trial, water_hi, energy_dummy, adiabat_profile_ok, & ! (out) + twater_hi, temp, qvap, dens, pres, zlev, int_weight, mix_mask, saturation_mask, nez, npz, debug_flag ) ! (in) - qvap_aj(:,:) = qvap(:,:) - where ( adjust_mask(:,:) ) - qvap_aj(:,:) = qsat_aj(:,:) - end where + if ( adiabat_profile_ok ) then + water_hi_feasible = ( water_hi <= water_mass_ini + water_tol ) + else + water_hi_feasible = .false. + water_hi = huge(1.0_RP) - 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) + if ( debug_flag ) then + write(*,*) "High-temperature trial profile is not admissible." + write(*,*) " Treating it as the infeasible upper bracket: twater_hi=", twater_hi + end if + end if + + !- + + if ( debug_flag ) then + write(*,*) "Water feasibility at initial bounds:" + write(*,*) " t_lo, water_lo =", twater_lo, water_lo + write(*,*) " t_hi, water_hi =", twater_hi, water_hi + write(*,*) " lower feasible =", water_lo_feasible + write(*,*) " upper feasible =", water_hi_feasible + write(*,*) " water_mass_ini =", water_mass_ini + write(*,*) " water_tol =", water_tol + end if + + + if ( water_hi_feasible ) then + ! If the high side is still feasible, the chosen initial temperature range does not yet contain the water or thermodynamic upper boundary. + tbase_water_max = twater_hi + else + + ! Search for the maximum base temperature for which the moist-adiabat construction succeeds + ! and the resulting vapor mass remains feasible. + ! The lower endpoint is always: + ! - thermodynamically admissible + ! - water feasible + ! + ! The upper endpoint is either: + ! - water infeasible, or + ! - thermodynamically inadmissible + + do iter_water = 1, MCA_MAX_WATER_ITER + twater_mid = 0.5_RP * (twater_lo + twater_hi) + + call evaluate_trial_moist_neutral_profile( temp_trial, qvap_trial, water_mid, energy_dummy, adiabat_profile_ok, & ! (out) + twater_mid, temp, qvap, dens, pres, zlev, int_weight, mix_mask, saturation_mask, nez, npz, debug_flag ) ! (in) + + if ( .not. adiabat_profile_ok ) then + ! Thermodynamically inadmissible midpoint: + ! move the upper boundary downward. + twater_hi = twater_mid + water_hi = huge(1.0_RP) + + if ( debug_flag ) then + write(*,*) "water iter=", iter_water + write(*,*) " inadmissible midpoint=", twater_mid, " new bracket=", twater_lo, twater_hi + end if + + if ( abs(twater_hi - twater_lo) <= MCA_TEMP_TOL ) exit + + cycle end if + + if ( water_mid <= water_mass_ini + water_tol ) then + ! Admissible and water feasible. + twater_lo = twater_mid + water_lo = water_mid + else + ! Admissible but water infeasible. + twater_hi = twater_mid + water_hi = water_mid + end if + + if ( debug_flag ) then + write(*,*) "water iter=", iter_water, ", tbase=", twater_lo, twater_mid, twater_hi, ", water=", water_lo, water_mid, water_hi, ", target=", water_mass_ini + end if + + if ( abs(twater_hi - twater_lo) <= MCA_TEMP_TOL ) exit + if ( abs(water_mid - water_mass_ini) <= water_tol ) exit 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 + ! The lower endpoint is maintained on the valid and water-feasible side. + tbase_water_max = twater_lo + end if + + + !- Step 2: + ! Evaluate the energy residual at both ends of the water-feasible interval. + + tbase_lo = tbase_mid - MCA_TBASE_RANGE + tbase_hi = tbase_water_max + + !- + call evaluate_trial_moist_neutral_profile( temp_trial, qvap_trial, water_lo, energy_lo, adiabat_profile_ok, & ! (out) + tbase_lo, temp, qvap, dens, pres, zlev, int_weight, mix_mask, saturation_mask, nez, npz, debug_flag ) ! (in) + + if ( .not. adiabat_profile_ok ) then + profile_status = MCA_PROFILE_ADIABAT_FAILURE + return + end if + + residual_lo = energy_lo - energy_target + + !- + call evaluate_trial_moist_neutral_profile( temp_trial, qvap_trial, water_hi, energy_hi, adiabat_profile_ok, & ! (out) + tbase_hi, temp, qvap, dens, pres, zlev, int_weight, mix_mask, saturation_mask, nez, npz, debug_flag ) ! (in) + + if ( .not. adiabat_profile_ok ) then + profile_status = MCA_PROFILE_ADIABAT_FAILURE + + if ( debug_flag ) then + write(*,*) "Internal inconsistency:" + write(*,*) " tbase_water_max should be admissible but evaluation failed." + write(*,*) " tbase_hi=", tbase_hi end if - water_feasible = .true. + return + end if - call column_mse( energy_aj, & ! (out) - temp_aj, qvap_aj, zlev, dens, & ! (in) - int_weight, adjust_mask, npz, nez ) ! (in) + residual_hi = energy_hi - energy_target + + if ( debug_flag ) then + write(*,*) "Energy-root feasibility check:" + write(*,*) " energy_target =", energy_target + write(*,*) " energy_tol =", energy_tol + write(*,*) " tbase_lo =", tbase_lo + write(*,*) " energy_lo =", energy_lo + write(*,*) " residual_lo =", residual_lo + write(*,*) " tbase_water_max =", tbase_hi + write(*,*) " energy_water_max =", energy_hi + write(*,*) " residual_water_max =", residual_hi + write(*,*) " water_at_upper =", water_hi + write(*,*) " water_mass_ini =", water_mass_ini + end if + + ! Check whether either endpoint is already an energy root. + if ( abs(residual_lo) <= boundary_energy_tol ) then + + call evaluate_trial_moist_neutral_profile( temp_adj, qvap_adj, water_mass_adj, energy_trial, adiabat_profile_ok, & ! (out) + tbase_lo, temp, qvap, dens, pres, zlev, int_weight, mix_mask, saturation_mask, nez, npz, debug_flag ) ! (in) + + if ( .not. adiabat_profile_ok ) then + profile_status = MCA_PROFILE_ADIABAT_FAILURE + return + end if + converged = .true. + profile_status = MCA_PROFILE_SUCCESS - 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 + else if ( abs(residual_hi) <= boundary_energy_tol ) then - if ( abs(residual) <= MCA_ENERGY_RTOL * max(abs(energy_target),1.0_RP) ) then - converged = .true. - exit + call evaluate_trial_moist_neutral_profile( temp_adj, qvap_adj, water_mass_adj, energy_trial, adiabat_profile_ok, & ! (out) + tbase_hi, temp, qvap, dens, pres, zlev, int_weight, mix_mask, saturation_mask, nez, npz, debug_flag ) ! (in) + + if ( .not. adiabat_profile_ok ) then + profile_status = MCA_PROFILE_ADIABAT_FAILURE + return end if - if ( residual > 0.0_RP ) then - tbase_hi = tbase_mid - else - tbase_lo = tbase_mid + converged = .true. + profile_status = MCA_PROFILE_SUCCESS + + else if ( residual_lo * residual_hi > 0.0_RP ) then + profile_status = MCA_PROFILE_NO_ENERGY_ROOT + if ( debug_flag ) then + write(*,*) "No energy root exists in the water-feasible interval." + write(*,*) " tbase_lo =", tbase_lo + write(*,*) " tbase_water_max =", tbase_hi + write(*,*) " residual_lo =", residual_lo + write(*,*) " residual_upper =", residual_hi + write(*,*) " water_mass_ini =", water_mass_ini + write(*,*) " water_mass_upper =", water_hi end if - end do ! end loop for iteration + return + end if + !- Step 3: + ! If the endpoint check found a valid bracket, solve the column-energy constraint by bisection. + 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." + do iter=1, MCA_MAX_ENERGY_ITER + + tbase_mid = 0.5_RP * ( tbase_lo + tbase_hi ) + + call evaluate_trial_moist_neutral_profile( temp_adj, qvap_adj, water_trial, energy_trial, adiabat_profile_ok, & ! (out) + tbase_mid, temp, qvap, dens, pres, zlev, int_weight, mix_mask, saturation_mask, nez, npz, debug_flag ) ! (in) + + if ( .not. adiabat_profile_ok ) then + profile_status = MCA_PROFILE_ADIABAT_FAILURE + + if ( debug_flag ) then + write(*,*) "Moist-adiabat integration failed during" + write(*,*) "energy-root search: iteration=", iter, ", tbase=", tbase_mid + end if + return + end if + + ! The energy search must remain inside the previously determined water-feasible interval. + + if ( water_trial > water_mass_ini + water_tol ) then + profile_status = MCA_PROFILE_NO_ENERGY_ROOT + if ( debug_flag ) then + write(*,*) "Energy bisection entered a water-infeasible state." + write(*,*) " tbase =", tbase_mid + write(*,*) " water_trial =", water_trial + write(*,*) " water_mass_ini =", water_mass_ini + end if + return + end if + + residual = energy_trial - energy_target + + if ( debug_flag ) then + write(*,*) "energy iter=", iter + write(*,*) " tbase=", tbase_lo, tbase_mid, tbase_hi + write(*,*) " energy_trial=", energy_trial, ", energy_target=", energy_target + write(*,*) " residual=", residual + write(*,*) " water_trial=", water_trial + end if + + if ( abs(residual) <= energy_tol ) then + water_mass_adj = water_trial + converged = .true. + profile_status = MCA_PROFILE_SUCCESS + exit + end if + + ! General bisection update based on the residual signs. + + if ( residual_lo * residual <= 0.0_RP ) then + tbase_hi = tbase_mid + residual_hi = residual + else + tbase_lo = tbase_mid + residual_lo = residual + end if + + end do ! end loop for iteration + end if + + if ( .not. converged ) then + profile_status = MCA_PROFILE_ENERGY_MAXITER + + if ( debug_flag ) then + write(*,*) "Energy bisection did not converge." + write(*,*) " tbase_lo=", tbase_lo, ", tbase_hi=", tbase_hi + write(*,*) " residual_lo=", residual_lo, ", residual_hi=", residual_hi + write(*,*) " energy_tol=", energy_tol end if return end if - ! Convert the converged state to density variables + ! Step 4: + ! Convert the converged vapor mixing ratio to density form and diagnose precipitation. - rhoqvap_aj(:,:) = dens(:,:) * qvap_aj(:,:) + rhoqvap_aj(:,:) = dens(:,:) * qvap_adj(:,:) - precip_mass_total = max( 0.0_RP, water_mass_ini - water_mass_aj ) + precip_mass_total = max( 0.0_RP, water_mass_ini - water_mass_adj ) positive_cond_mass = 0.0_RP do ke_z=1, nez do pz=1, npz - if ( adjust_mask(pz,ke_z) ) then + if ( mix_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) @@ -671,13 +1136,13 @@ subroutine construct_moist_neutral_profile( & end do end do - ! Distribute precipitation mass over locally condensing nodes + ! Distribute the column-integrated removed water mass in proportion to positive local vapor-density reductions. 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 + if ( mix_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 @@ -686,64 +1151,154 @@ subroutine construct_moist_neutral_profile( & 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_ - !------------------------------------------------ + end subroutine solve_moist_neutral_adjustment - energy = 0.0_RP +!OCL SERIAL + subroutine evaluate_trial_moist_neutral_profile( & + temp_work, qvap_work, water_mass, energy, profile_ok, & ! (out) + tbase, temp, qvap, dens, pres, zlev, int_weight, & ! (in) + adj_mask, forced_saturation_mask, nez, npz, debug_flag ) ! (in) + implicit none - 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_) + integer, intent(in) :: nez + integer, intent(in) :: npz + real(RP), intent(out) :: temp_work(npz,nez) + real(RP), intent(out) :: qvap_work(npz,nez) + real(RP), intent(out) :: water_mass + real(RP), intent(out) :: energy + logical, intent(out) :: profile_ok + real(RP), intent(in) :: tbase + real(RP), intent(in) :: temp(npz,nez) + 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) :: int_weight(npz,nez) + logical, intent(in) :: adj_mask(npz,nez) + logical, intent(in) :: forced_saturation_mask(npz,nez) + logical, intent(in) :: debug_flag + + real(RP) :: qsat_work(npz,nez) - 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 + integer :: ke_z, pz + !------------------------------------------------ + + temp_work(:,:) = temp(:,:) + qvap_work(:,:) = qvap(:,:) + qsat_work(:,:) = 0.0_RP + + !- Evaluate the reference temperature profile over mix_mask. + + call build_moist_adiabat_temperature( & + temp_work, qsat_work, profile_ok, & ! (inout,out) + tbase, pres, zlev, adj_mask, npz, nez, debug_flag ) ! (in) + + if ( .not. profile_ok ) then + water_mass = huge(1.0_RP) + energy = huge(1.0_RP) return - end subroutine column_mse - end subroutine construct_moist_neutral_profile + end if + + !- Humidity construction: + ! 1. Nearly saturated nodes in the diagnosed unstable core: + ! qv = qsat. + ! 2. Other nodes inside the adjustment region: + ! retain the input qv unless the trial temperature produces supersaturation. + ! 3. Outside the adjustment region: + ! retain the input temperature and qv. + do ke_z = 1, nez + do pz = 1, npz + if ( adj_mask(pz,ke_z) ) then + if ( forced_saturation_mask(pz,ke_z) ) then + qvap_work(pz,ke_z) = qsat_work(pz,ke_z) + else + qvap_work(pz,ke_z) = min(qvap(pz,ke_z), qsat_work(pz,ke_z)) + end if + end if + qvap_work(pz,ke_z) = max(0.0_RP, qvap_work(pz,ke_z)) + end do + end do + + !- Column water over the full mixing region. + + water_mass = 0.0_RP + do ke_z = 1, nez + do pz = 1, npz + if ( adj_mask(pz,ke_z) ) then + water_mass = water_mass + int_weight(pz,ke_z) * dens(pz,ke_z) * qvap_work(pz,ke_z) + end if + end do + end do + + !- Column MSE over the full mixing region. + + call integ_masked_column_mse( energy, & ! (out) + temp_work, qvap_work, zlev, dens, & ! (in) + int_weight, adj_mask, npz, nez ) ! (in) + + return + end subroutine evaluate_trial_moist_neutral_profile + + ! The column-MSE constraint is imposed while holding the density and pressure profiles fixed. + ! Thermodynamic composition includes dry air and water vapor only. + ! Condensate sensible heat and the energy carried away by precipitation are neglected. +!OCL SERIAL + subroutine integ_masked_column_mse( energy, & ! (out) + temp, qv, zlev, dens, & ! (in) + int_weight, integ_mask, npz, nez ) ! (in) + 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) :: integ_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 ( integ_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 integ_masked_column_mse + + !---- !OCL SERIAL - subroutine evaluate_moist_adiabat_temperature( & - temp_ref, qsat_ref, is_converged_ref_profile, & ! (out) - tbase, pres, zlev, mask, npz, nez ) ! (in) + subroutine build_moist_adiabat_temperature( & + temp_profile, qsat_profile, profile_converged, & ! (out) + tbase, pres, zlev, mask, npz, nez, debug_flag ) ! (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(inout) :: temp_profile(npz,nez) + real(RP), intent(inout) :: qsat_profile(npz,nez) + logical, intent(out) :: profile_converged 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) + logical, intent(in) :: debug_flag integer :: ke_z, pz @@ -755,7 +1310,7 @@ subroutine evaluate_moist_adiabat_temperature( & !---------------------------------------------------- started = .false. - is_converged_ref_profile = .true. + profile_converged = .true. do ke_z = 1, nez do pz = 1, npz @@ -767,22 +1322,21 @@ subroutine evaluate_moist_adiabat_temperature( & 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 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) + call solve_next_moist_adiabat_temp( 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. + profile_converged = .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) + temp_profile(pz,ke_z) = t_now + call ATMOS_SATURATION_pres2qsat_liq( t_now, p_now, & ! (in) + qsat_profile(pz,ke_z) ) ! (out) t_prev = t_now p_prev = p_now @@ -791,11 +1345,12 @@ subroutine evaluate_moist_adiabat_temperature( & end do return - end subroutine evaluate_moist_adiabat_temperature + end subroutine build_moist_adiabat_temperature + !> Solve the temperature at (pres1,zlev1) such that the saturated MSE equals that at (temp0,pres0,zlev0). !OCL SERIAL - subroutine advance_moist_adiabat( temp1, converged, & ! (out) - temp0, pres0, zlev0, pres1, zlev1 ) ! (in) + subroutine solve_next_moist_adiabat_temp( temp1, converged, & ! (out) + temp0, pres0, zlev0, pres1, zlev1 ) ! (in) implicit none real(RP), intent(out) :: temp1 logical, intent(out) :: converged @@ -804,108 +1359,89 @@ subroutine advance_moist_adiabat( temp1, converged, & ! (out) real(RP), intent(in) :: zlev0 real(RP), intent(in) :: pres1 real(RP), intent(in) :: zlev1 + + real(RP) :: temp_lo, temp_mid, temp_hi + real(RP) :: f_lo, f_mid, f_hi + real(RP) :: hs_target, hs_mid + real(RP) :: temp_tol, energy_tol + 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 + integer, parameter :: MAX_ITER = 60 + real(RP), parameter :: TEMP_RANGE = 40.0_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 + call saturated_mse_point( temp0, pres0, zlev0, & + hs_target ) - 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 + temp_lo = temp0 - TEMP_RANGE + temp_hi = temp0 + TEMP_RANGE - !- 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) + call saturated_mse_point( temp_lo, pres1, zlev1, & + hs_mid ) + f_lo = hs_mid - hs_target - tguess = temp0 - Gamma_m * dz - tnew = tguess + call saturated_mse_point( temp_hi, pres1, zlev1, & + hs_mid ) + f_hi = hs_mid - hs_target - !---------------------------------------------- - ! Iteration to solve the implicit equation for moist adiabatic lapse rate - ! T1^* = T0 - Gamma_m ( (T0 + T1^*)/2, pmid ) * dz - !---------------------------------------------- + if ( f_lo * f_hi > 0.0_RP ) then + return + end if - do iter = 1, MCA_MAX_LOCAL_ITER - - t_mid = 0.5_RP * ( temp0 + tguess ) + energy_tol = MCA_LOCAL_HMSE_ATOL + MCA_LOCAL_HMSE_RTOL * abs(hs_target) + temp_tol = MCA_TEMP_TOL - call ATMOS_SATURATION_pres2qsat_liq( t_mid, p_mid, & ! (in) - qsat_mid ) ! (out) + do iter = 1, MAX_ITER + temp_mid = 0.5_RP * ( temp_lo + temp_hi ) - call moist_adiabatic_lapse_rate( Gamma_m, & ! (out) - t_mid, qsat_mid ) ! (in) + call saturated_mse_point( temp_mid, pres1, zlev1, & + hs_mid ) - tnew = temp0 - Gamma_m * dz + f_mid = hs_mid - hs_target - 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 + if ( f_lo * f_mid <= 0.0_RP ) then + temp_hi = temp_mid + f_hi = f_mid + else + temp_lo = temp_mid + f_lo = f_mid + 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 + if ( ( abs(temp_hi - temp_lo) <= temp_tol ) & + .or. ( abs(f_mid) <= energy_tol ) ) then + temp1 = temp_mid + converged = .true. + return + end if + end do - temp1 = tnew + temp1 = temp_mid return contains - subroutine moist_adiabatic_lapse_rate( GammaM, & - temp, qsat ) + subroutine saturated_mse_point( & + temp, pres, zlev, hmse_sat ) 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 + real(RP), intent(in) :: temp + real(RP), intent(in) :: pres + real(RP), intent(in) :: zlev + real(RP), intent(out) :: hmse_sat + + real(RP) :: qsat + real(RP) :: cptot !------------------------------------------------ - 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 ) ) + call ATMOS_SATURATION_pres2qsat_liq( & + temp, pres, qsat ) + + cptot = CPdry * (1.0_RP-qsat) + CP_VAPOR * qsat + + hmse_sat = cptot * temp + GRAV * zlev + LHV0 * qsat return - end subroutine moist_adiabatic_lapse_rate - end subroutine advance_moist_adiabat + end subroutine saturated_mse_point + end subroutine solve_next_moist_adiabat_temp + end module scale_atm_phy_cp_dgm_mconv_adjustment From 9a3fe8fe799b550654de2245fe4efe33f47ecaba Mon Sep 17 00:00:00 2001 From: ykawai Date: Wed, 5 Aug 2026 20:24:10 +0900 Subject: [PATCH 5/6] Update auxiliary variable management and fix array dimensions --- model/atm_nonhydro3d/src/atmos/mod_atmos_vars.F90 | 5 +++++ model/atm_nonhydro3d/src/ocean/mod_ocean_component.F90 | 10 +++++----- 2 files changed, 10 insertions(+), 5 deletions(-) diff --git a/model/atm_nonhydro3d/src/atmos/mod_atmos_vars.F90 b/model/atm_nonhydro3d/src/atmos/mod_atmos_vars.F90 index dc306f08..8e58100e 100644 --- a/model/atm_nonhydro3d/src/atmos/mod_atmos_vars.F90 +++ b/model/atm_nonhydro3d/src/atmos/mod_atmos_vars.F90 @@ -375,6 +375,11 @@ subroutine AtmosVars_Init( this, atm_mesh ) LOG_INFO("ATMOS_vars_setup",*) 'Check value range of variables? : ', CHECK_RANGE LOG_INFO("ATMOS_vars_setup",*) 'Check total value of variables? : ', CHECK_TOTAL + !-- Set the pointer of 2D auxiliary variable manager with MP and CP components to output precipitation fluxes + + nullify( this%ptr_MP_AUXVARS2D_manager ) + nullify( this%ptr_CP_AUXVARS2D_manager ) + return end subroutine AtmosVars_Init diff --git a/model/atm_nonhydro3d/src/ocean/mod_ocean_component.F90 b/model/atm_nonhydro3d/src/ocean/mod_ocean_component.F90 index 6c590865..e2de1c40 100644 --- a/model/atm_nonhydro3d/src/ocean/mod_ocean_component.F90 +++ b/model/atm_nonhydro3d/src/ocean/mod_ocean_component.F90 @@ -518,10 +518,10 @@ subroutine calculate_surface_flux( this, & integer :: ke2D, ij real(RP) :: DZ1 (elem%Np,lmesh%NeA) - real(RP) :: Z1 (elem%Np,lmesh%NeA) - real(RP) :: ATM_W_lc(elem%Np,lmesh%Ne) - real(RP) :: ATM_U_lc(elem%Np,lmesh%Ne) - real(RP) :: ATM_V_lc(elem%Np,lmesh%Ne) + real(RP) :: Z1 (elem%Np,lmesh%Ne) + real(RP) :: ATM_W_lc(elem%Np,lmesh%NeA) + real(RP) :: ATM_U_lc(elem%Np,lmesh%NeA) + real(RP) :: ATM_V_lc(elem%Np,lmesh%NeA) !-------------------------------------------------- !$omp parallel do collapse(2) @@ -540,7 +540,7 @@ subroutine calculate_surface_flux( this, & call convert_UV2LocalOrthVec( & this%mesh%ptr_mesh, lmesh%pos_en(:,:,1), lmesh%pos_en(:,:,2), Z1(:,:), elem%Np*lmesh%Ne, & - ATM_U_lc, ATM_V_lc ) ! (inout) + ATM_U_lc(:,lmesh%NeS:lmesh%NeE), ATM_V_lc(:,lmesh%NeS:lmesh%NeE) ) ! (inout) call ATMOS_PHY_SF_simple_flux( & elem%Np, 1, elem%Np, lmesh%NeA, 1, lmesh%Ne, & ! (in) From 1fdb342e9872d4e7a81540a9891d4c41d4d37004 Mon Sep 17 00:00:00 2001 From: ykawai Date: Wed, 5 Aug 2026 20:27:52 +0900 Subject: [PATCH 6/6] Fix OpenACC directives --- model/atm_nonhydro3d/src/atmos/mod_atmos_vars.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/model/atm_nonhydro3d/src/atmos/mod_atmos_vars.F90 b/model/atm_nonhydro3d/src/atmos/mod_atmos_vars.F90 index 8e58100e..8ba07fd1 100644 --- a/model/atm_nonhydro3d/src/atmos/mod_atmos_vars.F90 +++ b/model/atm_nonhydro3d/src/atmos/mod_atmos_vars.F90 @@ -1018,7 +1018,7 @@ subroutine vars_calc_diagnoseVar2D_lc( field_name, & ! (in) select case(trim(field_name)) case('RAIN') !$omp parallel - !$acc parallel present(SFLX_rain_MP%val, var_out) async(1) + !$acc parallel present(var_out) async(1) !$omp do !$acc loop collapse(2) do ke=lcmesh%NeS, lcmesh%NeE @@ -1048,7 +1048,7 @@ subroutine vars_calc_diagnoseVar2D_lc( field_name, & ! (in) !$acc end parallel case('SNOW') !$omp parallel - !$acc parallel present(SFLX_snow_MP%val, var_out) async(1) + !$acc parallel present(var_out) async(1) !$omp do !$acc loop collapse(2) do ke=lcmesh%NeS, lcmesh%NeE