From 755430ba686194511c13379049fe7362529b2b2a Mon Sep 17 00:00:00 2001 From: Tian Zhou Date: Sat, 17 May 2025 00:01:34 -0700 Subject: [PATCH 01/17] Add scheme to redirect total Qgwl to top 100 river outlets --- components/mosart/bld/build-namelist | 1 + .../namelist_defaults_mosart.xml | 1 + .../namelist_definition_mosart.xml | 8 + components/mosart/src/cpl/rof_comp_mct.F90 | 2 +- components/mosart/src/riverroute/RtmMod.F90 | 305 +++++++++++++++++- components/mosart/src/riverroute/RtmVar.F90 | 1 + 6 files changed, 311 insertions(+), 7 deletions(-) diff --git a/components/mosart/bld/build-namelist b/components/mosart/bld/build-namelist index 2ce54f39a2c6..9175d3d12136 100755 --- a/components/mosart/bld/build-namelist +++ b/components/mosart/bld/build-namelist @@ -330,6 +330,7 @@ else { add_default($nl, 'RoutingMethod'); add_default($nl, 'DLevelH2R'); add_default($nl, 'DLevelR'); + add_default($nl, 'redirect_total_qgwl'); # add_default($nl, 'finidat_rtm' , 'rof_grid'=>$ROF_GRID, 'simyr'=>$opts{'simyr'}, 'nofail'=>1 ); my $val = 0; if ($NCPL_BASE_PERIOD eq "year") { diff --git a/components/mosart/bld/namelist_files/namelist_defaults_mosart.xml b/components/mosart/bld/namelist_files/namelist_defaults_mosart.xml index b3e2ad8c9269..e8e0b9f0a1e2 100644 --- a/components/mosart/bld/namelist_files/namelist_defaults_mosart.xml +++ b/components/mosart/bld/namelist_files/namelist_defaults_mosart.xml @@ -24,6 +24,7 @@ for the CLM data in the CESM distribution .false. 0 .false. +.false. .false. 283.15 50 diff --git a/components/mosart/bld/namelist_files/namelist_definition_mosart.xml b/components/mosart/bld/namelist_files/namelist_definition_mosart.xml index c21879908ed6..1d3b54fe068a 100644 --- a/components/mosart/bld/namelist_files/namelist_definition_mosart.xml +++ b/components/mosart/bld/namelist_files/namelist_definition_mosart.xml @@ -126,6 +126,14 @@ Default: .false. If .true., mosart will read in bgc nutrients and send to ocean. + +Default: .false. +If .true., mosart will redirect total qgwl to the ocean. + + 0) then + allocate(all_outlet_gindices(num_all_outlets_global)) + allocate(all_outlet_discharges(num_all_outlets_global)) + else + ! Handle cases where no outlets are found if necessary, or let current_top_n_count be 0 + num_all_outlets_global = 0 ! Ensure it's zero if no outlets + endif + else + num_all_outlets_global = 0 ! Non-master tasks initialize to 0 + endif + call MPI_Bcast(num_all_outlets_global, 1, MPI_INTEGER, mastertask, mpicom_rof, ier) + + if(num_all_outlets_global > 0) then + call MPI_Gatherv(outlet_gindices_local, local_outlet_count, MPI_INTEGER, & + all_outlet_gindices, recvcounts, displs, MPI_INTEGER, & + mastertask, mpicom_rof, ier) + call MPI_Gatherv(outlet_discharges_local, local_outlet_count, MPI_REAL8, & + all_outlet_discharges, recvcounts, displs, MPI_REAL8, & + mastertask, mpicom_rof, ier) + endif + + if (allocated(outlet_gindices_local)) deallocate(outlet_gindices_local) + if (allocated(outlet_discharges_local)) deallocate(outlet_discharges_local) + + ! Rank on master, calculate proportions, and prepare corrections + current_top_n_count = 0 ! Initialize for all tasks + if (masterproc) then + if (num_all_outlets_global > 0) then + allocate(sorted_outlets(num_all_outlets_global)) + do i = 1, num_all_outlets_global + sorted_outlets(i)%gidx = all_outlet_gindices(i) + sorted_outlets(i)%discharge = all_outlet_discharges(i) + enddo + if (allocated(all_outlet_gindices)) deallocate(all_outlet_gindices) + if (allocated(all_outlet_discharges)) deallocate(all_outlet_discharges) + + call sort_outlets_by_discharge_desc(sorted_outlets, num_all_outlets_global) + + current_top_n_count = min(num_top_outlets_for_qgwl, num_all_outlets_global) + sum_discharge_top_n = 0.0_r8 + do i = 1, current_top_n_count + sum_discharge_top_n = sum_discharge_top_n + sorted_outlets(i)%discharge + enddo + if(current_top_n_count > 0) then + if (abs(sum_discharge_top_n) > 1.0e-9_r8) then ! Avoid division by zero + block + real(r8) :: qgwl_to_discharge_ratio_percent + + qgwl_to_discharge_ratio_percent = (global_net_qgwl / sum_discharge_top_n) * 100.0_r8 + + ! Print the ratio + write(iulog, *) trim(subname), & + 'Debug QGWL Ratio: (GlobalNetQGWL / SumTopNDischarge) = ', & + qgwl_to_discharge_ratio_percent, '%', & + ' (N_used = ', current_top_n_count, & + ', N_param = ', num_top_outlets_for_qgwl, ')' ! Using your var name + + ! Check and print warning if magnitude is > 5% + if (abs(qgwl_to_discharge_ratio_percent) > 5.0_r8) then + call shr_sys_flush(iulog) ! Flush previous message + write(iulog, *) trim(subname), & + 'WARNING: QGWL_Redist_Ratio magnitude > 5% (', & + qgwl_to_discharge_ratio_percent, '%). ', & + 'N_used = ', current_top_n_count, & + ', N_param = ', num_top_outlets_for_qgwl, & + '. Consider increasing N_param.' + call shr_sys_flush(iulog) ! Flush warning + endif + end block + else + write(iulog, *) trim(subname), & + 'Debug QGWL Ratio: Sum of Top ', current_top_n_count, & + ' Outlet Discharges is near zero. Ratio not calculated.' + endif + else + write(iulog, *) trim(subname), & + 'Debug QGWL Ratio: No top outlets selected (N_used = 0). Redistribution skipped.' + endif + if (current_top_n_count > 0) then + allocate(qgwl_correction_gindices(current_top_n_count)) + allocate(qgwl_correction_values(current_top_n_count)) + qgwl_correction_values = 0.0_r8 + do i = 1, current_top_n_count + qgwl_correction_gindices(i) = sorted_outlets(i)%gidx + qgwl_correction_values(i) = (sorted_outlets(i)%discharge / sum_discharge_top_n) * global_net_qgwl + enddo + else + write(iulog, '(A)') trim(subname), & + 'Debug: No top outlets selected for QGWL redistribution (current_top_n_count is 0).' + allocate(qgwl_correction_gindices(current_top_n_count)) + allocate(qgwl_correction_values(current_top_n_count)) + qgwl_correction_values = 0.0_r8 + if (sum_discharge_top_n > 1.0e-9_r8) then ! Avoid division by zero + do i = 1, current_top_n_count + qgwl_correction_gindices(i) = sorted_outlets(i)%gidx + qgwl_correction_values(i) = (sorted_outlets(i)%discharge / sum_discharge_top_n) * global_net_qgwl + enddo + else ! If sum of top N discharge is zero (or tiny), distribute equally (or handle as error/no distribution) + do i = 1, current_top_n_count + qgwl_correction_gindices(i) = sorted_outlets(i)%gidx + qgwl_correction_values(i) = global_net_qgwl / real(current_top_n_count, kind=r8) + enddo + endif + endif + if(allocated(sorted_outlets)) deallocate(sorted_outlets) + else + current_top_n_count = 0 ! No outlets globally + endif ! num_all_outlets_global > 0 + endif ! masterproc + + ! Broadcast correction info + call MPI_Bcast(current_top_n_count, 1, MPI_INTEGER, mastertask, mpicom_rof, ier) + if (.not. masterproc .and. current_top_n_count > 0) then + if (allocated(qgwl_correction_gindices)) deallocate(qgwl_correction_gindices) + if (allocated(qgwl_correction_values)) deallocate(qgwl_correction_values) + allocate(qgwl_correction_gindices(current_top_n_count)) + allocate(qgwl_correction_values(current_top_n_count)) + endif + if (current_top_n_count > 0) then + call MPI_Bcast(qgwl_correction_gindices, current_top_n_count, MPI_INTEGER, mastertask, mpicom_rof, ier) + call MPI_Bcast(qgwl_correction_values, current_top_n_count, MPI_REAL8, mastertask, mpicom_rof, ier) + + ! Populate local correction array + do k = 1, current_top_n_count + target_gidx = qgwl_correction_gindices(k) + do nr_local = rtmCTL%begr, rtmCTL%endr + if (rtmCTL%gindex(nr_local) == target_gidx) then + qgwl_correction_local(nr_local) = qgwl_correction_local(nr_local) + qgwl_correction_values(k) + exit + endif + enddo + enddo + endif ! current_top_n_count > 0 for broadcast and local application + + if (allocated(qgwl_correction_gindices)) deallocate(qgwl_correction_gindices) + if (allocated(qgwl_correction_values)) deallocate(qgwl_correction_values) + if (allocated(recvcounts)) deallocate(recvcounts) + if (allocated(displs)) deallocate(displs) + + call t_stopf('mosartr_qgwl_redir_dist') + endif ! global_net_qgwl /= 0.0_r8 (or just global_net_qgwl condition removed as per your last request) + endif ! redirect_total_qgwl_flag + + + ! This loop should now use the qgwl_correction_local array + do nt = 1,nt_rtm + do nr = rtmCTL%begr,rtmCTL%endr + volr_init = rtmCTL%volr(nr,nt) + + rtmCTL%runoff(nr,nt) = flow(nr,nt) ! This is main channel outflow from routing + rtmCTL%runofftot(nr,nt) = rtmCTL%direct(nr,nt) ! This is ice, qdto (if flag off), etc. + + if (rtmCTL%mask(nr) == 1) then ! Land + rtmCTL%runofflnd(nr,nt) = rtmCTL%runoff(nr,nt) + rtmCTL%dvolrdtlnd(nr,nt)= rtmCTL%dvolrdt(nr,nt) + elseif (rtmCTL%mask(nr) >= 2) then ! Ocean or Outlet + rtmCTL%runoffocn(nr,nt) = rtmCTL%runoff(nr,nt) ! Routed component + + ! Apply the distributed global qgwl correction ONLY to liquid tracer at target outlets + if (redirect_total_qgwl_flag .and. nt == nt_nliq) then + rtmCTL%runoffocn(nr,nt) = rtmCTL%runoffocn(nr,nt) + qgwl_correction_local(nr) + endif + rtmCTL%runofftot(nr,nt) = rtmCTL%runofftot(nr,nt) + rtmCTL%runoffocn(nr,nt) ! Add potentially corrected routed flow + rtmCTL%dvolrdtocn(nr,nt)= rtmCTL%dvolrdt(nr,nt) + endif + enddo ! nr + enddo ! nt + if (allocated(qgwl_correction_local)) deallocate(qgwl_correction_local) + + call t_stopf('mosartr_subcycling') !----------------------------------- @@ -5073,5 +5348,23 @@ end subroutine SubTimestep !----------------------------------------------------------------------- -end module RtmMod +subroutine sort_outlets_by_discharge_desc(outlets_array, count) !TZ negative runoff + implicit none + integer, intent(in) :: count + type(outlet_discharge_info_type), intent(inout) :: outlets_array(:) ! Assumed-shape + integer :: i, j + type(outlet_discharge_info_type) :: temp_outlet_info + ! Bubble sort + if (count < 2) return + do i = 1, count - 1 + do j = 1, count - i + if (outlets_array(j)%discharge < outlets_array(j+1)%discharge) then + temp_outlet_info = outlets_array(j) + outlets_array(j) = outlets_array(j+1) + outlets_array(j+1) = temp_outlet_info + endif + enddo + enddo +end subroutine sort_outlets_by_discharge_desc +end module RtmMod \ No newline at end of file diff --git a/components/mosart/src/riverroute/RtmVar.F90 b/components/mosart/src/riverroute/RtmVar.F90 index 7818d3b4bdf6..7cb72100732b 100644 --- a/components/mosart/src/riverroute/RtmVar.F90 +++ b/components/mosart/src/riverroute/RtmVar.F90 @@ -40,6 +40,7 @@ module RtmVar integer, public :: nlayers = 30 ! Maximum number of reservoir layers for stratification logical, public :: noland = .false. ! true => no valid land points -- do NOT run logical, public :: data_bgc_fluxes_to_ocean_flag = .false.! read in and send BGC fluxes to ocean flag + logical, public :: redirect_total_qgwl = .false. ! redirect total qgwl flag character(len=32) , public :: decomp_option ! decomp option character(len=32) , public :: smat_option ! smatrix multiply option (opt, Xonly, Yonly) ! opt = XandY in MCT From 900af6794671e680e4c683aa9eb97547f8ca1326 Mon Sep 17 00:00:00 2001 From: Tian Zhou Date: Mon, 2 Jun 2025 17:53:08 -0700 Subject: [PATCH 02/17] Implementation to only redistribute negative Qgwl --- components/mosart/bld/build-namelist | 2 +- .../namelist_defaults_mosart.xml | 2 +- .../namelist_definition_mosart.xml | 4 +- components/mosart/src/cpl/rof_comp_mct.F90 | 2 +- components/mosart/src/riverroute/RtmMod.F90 | 72 +++++++++++-------- components/mosart/src/riverroute/RtmVar.F90 | 2 +- 6 files changed, 47 insertions(+), 37 deletions(-) diff --git a/components/mosart/bld/build-namelist b/components/mosart/bld/build-namelist index 9175d3d12136..74ffc6bc0abe 100755 --- a/components/mosart/bld/build-namelist +++ b/components/mosart/bld/build-namelist @@ -330,7 +330,7 @@ else { add_default($nl, 'RoutingMethod'); add_default($nl, 'DLevelH2R'); add_default($nl, 'DLevelR'); - add_default($nl, 'redirect_total_qgwl'); + add_default($nl, 'redirect_negative_qgwl'); # add_default($nl, 'finidat_rtm' , 'rof_grid'=>$ROF_GRID, 'simyr'=>$opts{'simyr'}, 'nofail'=>1 ); my $val = 0; if ($NCPL_BASE_PERIOD eq "year") { diff --git a/components/mosart/bld/namelist_files/namelist_defaults_mosart.xml b/components/mosart/bld/namelist_files/namelist_defaults_mosart.xml index e8e0b9f0a1e2..c68944769b2a 100644 --- a/components/mosart/bld/namelist_files/namelist_defaults_mosart.xml +++ b/components/mosart/bld/namelist_files/namelist_defaults_mosart.xml @@ -24,7 +24,7 @@ for the CLM data in the CESM distribution .false. 0 .false. -.false. +.false. .false. 283.15 50 diff --git a/components/mosart/bld/namelist_files/namelist_definition_mosart.xml b/components/mosart/bld/namelist_files/namelist_definition_mosart.xml index 1d3b54fe068a..41ce96eac36e 100644 --- a/components/mosart/bld/namelist_files/namelist_definition_mosart.xml +++ b/components/mosart/bld/namelist_files/namelist_definition_mosart.xml @@ -126,12 +126,12 @@ Default: .false. If .true., mosart will read in bgc nutrients and send to ocean. - Default: .false. -If .true., mosart will redirect total qgwl to the ocean. +If .true., mosart will redirect negative qgwl to the ocean. 1.0e-9_r8) then ! Avoid division by zero do i = 1, current_top_n_count qgwl_correction_gindices(i) = sorted_outlets(i)%gidx - qgwl_correction_values(i) = (sorted_outlets(i)%discharge / sum_discharge_top_n) * global_net_qgwl + qgwl_correction_values(i) = (sorted_outlets(i)%discharge / sum_discharge_top_n) * global_net_negative_qgwl enddo else ! If sum of top N discharge is zero (or tiny), distribute equally (or handle as error/no distribution) do i = 1, current_top_n_count qgwl_correction_gindices(i) = sorted_outlets(i)%gidx - qgwl_correction_values(i) = global_net_qgwl / real(current_top_n_count, kind=r8) + qgwl_correction_values(i) = global_net_negative_qgwl / real(current_top_n_count, kind=r8) enddo endif endif @@ -3086,8 +3096,8 @@ subroutine Rtmrun(rstwr,nlend,rdate) if (allocated(displs)) deallocate(displs) call t_stopf('mosartr_qgwl_redir_dist') - endif ! global_net_qgwl /= 0.0_r8 (or just global_net_qgwl condition removed as per your last request) - endif ! redirect_total_qgwl_flag + endif ! global_net_negative_qgwl /= 0.0_r8 (or just global_net_negative_qgwl condition removed as per your last request) + endif ! redirect_negative_qgwl_flag ! This loop should now use the qgwl_correction_local array @@ -3105,7 +3115,7 @@ subroutine Rtmrun(rstwr,nlend,rdate) rtmCTL%runoffocn(nr,nt) = rtmCTL%runoff(nr,nt) ! Routed component ! Apply the distributed global qgwl correction ONLY to liquid tracer at target outlets - if (redirect_total_qgwl_flag .and. nt == nt_nliq) then + if (redirect_negative_qgwl_flag .and. nt == nt_nliq) then rtmCTL%runoffocn(nr,nt) = rtmCTL%runoffocn(nr,nt) + qgwl_correction_local(nr) endif rtmCTL%runofftot(nr,nt) = rtmCTL%runofftot(nr,nt) + rtmCTL%runoffocn(nr,nt) ! Add potentially corrected routed flow diff --git a/components/mosart/src/riverroute/RtmVar.F90 b/components/mosart/src/riverroute/RtmVar.F90 index 7cb72100732b..64c8e6307159 100644 --- a/components/mosart/src/riverroute/RtmVar.F90 +++ b/components/mosart/src/riverroute/RtmVar.F90 @@ -40,7 +40,7 @@ module RtmVar integer, public :: nlayers = 30 ! Maximum number of reservoir layers for stratification logical, public :: noland = .false. ! true => no valid land points -- do NOT run logical, public :: data_bgc_fluxes_to_ocean_flag = .false.! read in and send BGC fluxes to ocean flag - logical, public :: redirect_total_qgwl = .false. ! redirect total qgwl flag + logical, public :: redirect_negative_qgwl = .false. ! redirect negative qgwl flag character(len=32) , public :: decomp_option ! decomp option character(len=32) , public :: smat_option ! smatrix multiply option (opt, Xonly, Yonly) ! opt = XandY in MCT From b96bd31374afe4ae679274d240307cf88e585602 Mon Sep 17 00:00:00 2001 From: Tian Zhou Date: Wed, 6 Aug 2025 16:40:04 -0700 Subject: [PATCH 03/17] implement offsetting scheme --- components/mosart/src/riverroute/RtmMod.F90 | 109 +++++++++++++------- 1 file changed, 71 insertions(+), 38 deletions(-) diff --git a/components/mosart/src/riverroute/RtmMod.F90 b/components/mosart/src/riverroute/RtmMod.F90 index 32985aed296a..677fb45bfd49 100644 --- a/components/mosart/src/riverroute/RtmMod.F90 +++ b/components/mosart/src/riverroute/RtmMod.F90 @@ -142,7 +142,8 @@ module RtmMod ! Variables for TNR redirection integer, parameter :: num_top_outlets_for_qgwl = 500 ! Number of top outlets to use - real(r8), save :: global_net_negative_qgwl = 0.0_r8 ! Global Total negative qgwl + real(r8), save :: global_positive_qgwl_sum = 0.0_r8 ! Sum of all positive qgwl + real(r8), save :: global_negative_qgwl_sum = 0.0_r8 ! Sum of all negative qgwl real(r8), save :: delt_save ! previous delt @@ -2238,7 +2239,10 @@ subroutine Rtmrun(rstwr,nlend,rdate) !scs ! Variables for negative runoff redirection - real(r8) :: local_qgwl_sum + + real(r8) :: local_positive_qgwl_sum, local_negative_qgwl_sum + real(r8) :: net_global_qgwl, original_cell_qgwl, reduction + integer, allocatable :: outlet_gindices_local(:) ! Local array of global indices of outlets on this task real(r8), allocatable :: outlet_discharges_local(:) ! Local array of discharges for these outlets integer :: local_outlet_count @@ -2387,42 +2391,68 @@ subroutine Rtmrun(rstwr,nlend,rdate) ! data for euler solver, in m3/s here ! Aggregate global Qgwl (TNR) if flag is true --- - global_net_negative_qgwl = 0.0_r8 ! Reset global sum each step - if (redirect_negative_qgwl_flag) then - local_qgwl_sum = 0.0_r8 ! This will sum local NEGATIVE qgwl - do nr = rtmCTL%begr, rtmCTL%endr - if (rtmCTL%qgwl(nr, nt_nliq) < 0.0_r8) then - local_qgwl_sum = local_qgwl_sum + rtmCTL%qgwl(nr, nt_nliq) - TRunoff%qgwl(nr, nt_nliq) = 0.0_r8 ! Negative part is removed for global redistribution - else - TRunoff%qgwl(nr, nt_nliq) = rtmCTL%qgwl(nr, nt_nliq) ! Positive part goes to local routing + global_negative_qgwl_sum = 0.0_r8 ! Reset global sum each step + + if (redirect_negative_qgwl_flag) then + local_negative_qgwl_sum = 0.0_r8 ! This will sum local NEGATIVE qgwl + local_positive_qgwl_sum = 0.0_r8 ! This will sum local POSITIVE qgwl + + do nr = rtmCTL%begr, rtmCTL%endr + if (rtmCTL%qgwl(nr, nt_nliq) > 0.0_r8) then + local_positive_qgwl_sum = local_positive_qgwl_sum + rtmCTL%qgwl(nr, nt_nliq) + elseif (rtmCTL%qgwl(nr, nt_nliq) < 0.0_r8) then + local_negative_qgwl_sum = local_negative_qgwl_sum + rtmCTL%qgwl(nr, nt_nliq) endif - ! Handle other tracers for qgwl if necessary - assuming qgwl for other tracers - ! should follow the original rtmCTL input if not the liquid one being modified. - do nt = 1, nt_rtm - if (nt /= nt_nliq) then - TRunoff%qgwl(nr, nt) = rtmCTL%qgwl(nr, nt) ! Other tracers pass through - endif - enddo - ! If qgwl for nt_nliq was positive, it's already set above. - ! If qgwl for nt_nliq was negative, it was zeroed out for nt_nliq. - enddo - - call MPI_Allreduce(local_qgwl_sum, global_net_negative_qgwl, 1, MPI_REAL8, MPI_SUM, mpicom_rof, ier) - if (ier /= MPI_SUCCESS) then - if (masterproc) write(iulog,*) trim(subname),' ERROR in MPI_Allreduce for global_net_negative_qgwl, ier=',ier - call shr_sys_abort(trim(subname)//' MPI_Allreduce error for global_net_negative_qgwl') - endif - + enddo + + ! Use direct MPI_Allreduce to get the global sums A and B on all processes + call MPI_Allreduce(local_positive_qgwl_sum, global_positive_qgwl_sum, 1, MPI_REAL8, MPI_SUM, mpicom_rof, ier) + call MPI_Allreduce(local_negative_qgwl_sum, global_negative_qgwl_sum, 1, MPI_REAL8, MPI_SUM, mpicom_rof, ier) + + net_global_qgwl = global_positive_qgwl_sum + global_negative_qgwl_sum + if (masterproc) then - write(iulog,'(A,ES15.6)') trim(subname)//' Global negative sum of Qgwl this step: ', global_net_negative_qgwl + write(iulog, *) trim(subname), 'Debug QGWL Balancing: Global Positive Sum =', global_positive_qgwl_sum + write(iulog, *) trim(subname), 'Debug QGWL Balancing: Global Negative Sum =', global_negative_qgwl_sum + write(iulog, *) trim(subname), 'Debug QGWL Balancing: Global Net Sum =', net_global_qgwl endif - else + + ! Decide how to set TRunoff%qgwl for local routing + if (net_global_qgwl >= 0.0_r8) then + ! --- SCENARIO A: Net global QGWL is non-negative. Balance the negative by proportionally reducing positive. --- + do nr = rtmCTL%begr, rtmCTL%endr + original_cell_qgwl = rtmCTL%qgwl(nr, nt_nliq) + if (original_cell_qgwl > 0.0_r8) then ! This cell has positive qgwl + if (global_positive_qgwl_sum > 1.0e-9_r8) then ! Avoid division by zero + ! Calculate how much this cell should contribute to offsetting the negative sum + reduction = (original_cell_qgwl / global_positive_qgwl_sum) * abs(global_negative_qgwl_sum) + TRunoff%qgwl(nr, nt_nliq) = max(0.0_r8, original_cell_qgwl - reduction) + else ! No positive qgwl anywhere, so this branch shouldn't be hit if original_cell_qgwl > 0 + TRunoff%qgwl(nr, nt_nliq) = original_cell_qgwl + endif + else ! This cell has negative or zero qgwl + TRunoff%qgwl(nr, nt_nliq) = 0.0_r8 ! Its negativity is balanced by the global positives + endif + ! Pass through other tracers unmodified + do nt = 1, nt_rtm + if (nt /= nt_nliq) then + TRunoff%qgwl(nr, nt) = rtmCTL%qgwl(nr, nt) + endif + enddo + enddo + else + ! --- SCENARIO B: Net global QGWL is negative. No local qgwl input. --- + ! The deficit (net_global_qgwl) will be redistributed to top N outlets later. + do nr = rtmCTL%begr, rtmCTL%endr + TRunoff%qgwl(nr, :) = 0.0_r8 + enddo + endif + else ! If flag is false, ensure TRunoff gets the qgwl from rtmCTL for routing do nr = rtmCTL%begr, rtmCTL%endr TRunoff%qgwl(nr, :) = rtmCTL%qgwl(nr, :) enddo - endif + endif do nr = rtmCTL%begr,rtmCTL%endr do nt = 1,nt_rtm @@ -2912,12 +2942,12 @@ subroutine Rtmrun(rstwr,nlend,rdate) enddo enddo - ! Collect outlet discharge, Rank, Redistribute global_net_negative_qgwl, Update discharge map --- + ! Collect outlet discharge, Rank, Redistribute net_global_qgwl, Update discharge map --- allocate(qgwl_correction_local(rtmCTL%begr:rtmCTL%endr)) ! Moved allocation here qgwl_correction_local = 0.0_r8 ! Initialize if (redirect_negative_qgwl_flag) then ! Check the main flag first - if (global_net_negative_qgwl < 0.0_r8) then ! Only proceed if there's some net qgwl to redistribute (can be pos or neg) + if (net_global_qgwl < 0.0_r8) then ! Only proceed if there's negative net qgwl to redistribute call t_startf('mosartr_qgwl_redir_dist') ! Identify local outlets and their current discharges @@ -3004,7 +3034,7 @@ subroutine Rtmrun(rstwr,nlend,rdate) block real(r8) :: qgwl_to_discharge_ratio_percent - qgwl_to_discharge_ratio_percent = (global_net_negative_qgwl / sum_discharge_top_n) * 100.0_r8 + qgwl_to_discharge_ratio_percent = (net_global_qgwl / sum_discharge_top_n) * 100.0_r8 ! Print the ratio write(iulog, *) trim(subname), & @@ -3040,7 +3070,7 @@ subroutine Rtmrun(rstwr,nlend,rdate) qgwl_correction_values = 0.0_r8 do i = 1, current_top_n_count qgwl_correction_gindices(i) = sorted_outlets(i)%gidx - qgwl_correction_values(i) = (sorted_outlets(i)%discharge / sum_discharge_top_n) * global_net_negative_qgwl + qgwl_correction_values(i) = (sorted_outlets(i)%discharge / sum_discharge_top_n) * net_global_qgwl enddo else write(iulog, '(A)') trim(subname), & @@ -3051,12 +3081,12 @@ subroutine Rtmrun(rstwr,nlend,rdate) if (sum_discharge_top_n > 1.0e-9_r8) then ! Avoid division by zero do i = 1, current_top_n_count qgwl_correction_gindices(i) = sorted_outlets(i)%gidx - qgwl_correction_values(i) = (sorted_outlets(i)%discharge / sum_discharge_top_n) * global_net_negative_qgwl + qgwl_correction_values(i) = (sorted_outlets(i)%discharge / sum_discharge_top_n) * net_global_qgwl enddo else ! If sum of top N discharge is zero (or tiny), distribute equally (or handle as error/no distribution) do i = 1, current_top_n_count qgwl_correction_gindices(i) = sorted_outlets(i)%gidx - qgwl_correction_values(i) = global_net_negative_qgwl / real(current_top_n_count, kind=r8) + qgwl_correction_values(i) = net_global_qgwl / real(current_top_n_count, kind=r8) enddo endif endif @@ -3096,7 +3126,10 @@ subroutine Rtmrun(rstwr,nlend,rdate) if (allocated(displs)) deallocate(displs) call t_stopf('mosartr_qgwl_redir_dist') - endif ! global_net_negative_qgwl /= 0.0_r8 (or just global_net_negative_qgwl condition removed as per your last request) + else + write(iulog, *) trim(subname), & + 'Global net QGWL is positive, offsetting the negative values.' + endif ! net_global_qgwl < 0.0_r8 endif ! redirect_negative_qgwl_flag From 8cb9508ead1a241180c5aeaf38a6b5734caac3e6 Mon Sep 17 00:00:00 2001 From: Tian Zhou Date: Sat, 18 Oct 2025 00:38:45 -0700 Subject: [PATCH 04/17] Reprosum instead of MPI_Allreduce --- components/mosart/src/riverroute/RtmMod.F90 | 174 ++++++++++++++++---- 1 file changed, 140 insertions(+), 34 deletions(-) diff --git a/components/mosart/src/riverroute/RtmMod.F90 b/components/mosart/src/riverroute/RtmMod.F90 index 677fb45bfd49..2758e406d097 100644 --- a/components/mosart/src/riverroute/RtmMod.F90 +++ b/components/mosart/src/riverroute/RtmMod.F90 @@ -15,6 +15,7 @@ module RtmMod use rof_cpl_indices , only : nt_rtm, rtm_tracers, KW, DW use seq_flds_mod , only : rof_sed use RtmSpmd , only : masterproc, npes, iam, mpicom_rof, ROFID, mastertask + use shr_reprosum_mod, only : shr_reprosum_calc use mpi use RtmVar , only : re, spval, rtmlon, rtmlat, iulog, ice_runoff, & frivinp_rtm, frivinp_mesh, finidat_rtm, nrevsn_rtm,rstraflag,ngeom,nlayers,rinittemp, & @@ -2255,7 +2256,7 @@ subroutine Rtmrun(rstwr,nlend,rdate) real(r8), allocatable :: qgwl_correction_values(:) ! Values to apply integer, allocatable :: qgwl_correction_gindices(:) ! Global indices for these corrections real(r8), allocatable :: qgwl_correction_local(:) ! Local portion of correction array - + real(r8) :: qgwl_to_redistribute ! Amount to redistribute (Scenario A or B) character(len=*),parameter :: subname = '(Rtmrun) ' !----------------------------------------------------------------------- @@ -2405,9 +2406,30 @@ subroutine Rtmrun(rstwr,nlend,rdate) endif enddo - ! Use direct MPI_Allreduce to get the global sums A and B on all processes - call MPI_Allreduce(local_positive_qgwl_sum, global_positive_qgwl_sum, 1, MPI_REAL8, MPI_SUM, mpicom_rof, ier) - call MPI_Allreduce(local_negative_qgwl_sum, global_negative_qgwl_sum, 1, MPI_REAL8, MPI_SUM, mpicom_rof, ier) + ! Use reproducible sum for bit-for-bit reproducibility across PE layouts + block + real(r8) :: pos_local(1,1), pos_global(1) + real(r8) :: neg_local(1,1), neg_global(1) + + ! Reproducible sum for positive qgwl + pos_local(1,1) = local_positive_qgwl_sum + call shr_reprosum_calc(pos_local, pos_global, 1, 1, 1, & + commid=mpicom_rof) + global_positive_qgwl_sum = pos_global(1) + + ! Reproducible sum for negative qgwl + neg_local(1,1) = local_negative_qgwl_sum + call shr_reprosum_calc(neg_local, neg_global, 1, 1, 1, & + commid=mpicom_rof) + global_negative_qgwl_sum = neg_global(1) + + ! Diagnostic output for debugging + if (masterproc) then + write(iulog, '(A,ES24.16)') trim(subname)//' REPROSUM Positive Sum = ', global_positive_qgwl_sum + write(iulog, '(A,ES24.16)') trim(subname)//' REPROSUM Negative Sum = ', global_negative_qgwl_sum + write(iulog, '(A,ES24.16)') trim(subname)//' REPROSUM Net Sum = ', global_positive_qgwl_sum + global_negative_qgwl_sum + endif + end block net_global_qgwl = global_positive_qgwl_sum + global_negative_qgwl_sum @@ -2418,31 +2440,18 @@ subroutine Rtmrun(rstwr,nlend,rdate) endif ! Decide how to set TRunoff%qgwl for local routing + ! For BOTH scenarios, we let original qgwl flow to direct discharge, + ! and then add corrections at outlets to balance negative values if (net_global_qgwl >= 0.0_r8) then - ! --- SCENARIO A: Net global QGWL is non-negative. Balance the negative by proportionally reducing positive. --- + ! --- SCENARIO A: Net global QGWL is non-negative. --- + ! Let original qgwl flow to rtmCTL%direct (no modification here) + ! The negative qgwl will be balanced by adding to outlets later do nr = rtmCTL%begr, rtmCTL%endr - original_cell_qgwl = rtmCTL%qgwl(nr, nt_nliq) - if (original_cell_qgwl > 0.0_r8) then ! This cell has positive qgwl - if (global_positive_qgwl_sum > 1.0e-9_r8) then ! Avoid division by zero - ! Calculate how much this cell should contribute to offsetting the negative sum - reduction = (original_cell_qgwl / global_positive_qgwl_sum) * abs(global_negative_qgwl_sum) - TRunoff%qgwl(nr, nt_nliq) = max(0.0_r8, original_cell_qgwl - reduction) - else ! No positive qgwl anywhere, so this branch shouldn't be hit if original_cell_qgwl > 0 - TRunoff%qgwl(nr, nt_nliq) = original_cell_qgwl - endif - else ! This cell has negative or zero qgwl - TRunoff%qgwl(nr, nt_nliq) = 0.0_r8 ! Its negativity is balanced by the global positives - endif - ! Pass through other tracers unmodified - do nt = 1, nt_rtm - if (nt /= nt_nliq) then - TRunoff%qgwl(nr, nt) = rtmCTL%qgwl(nr, nt) - endif - enddo + TRunoff%qgwl(nr, :) = rtmCTL%qgwl(nr, :) enddo else - ! --- SCENARIO B: Net global QGWL is negative. No local qgwl input. --- - ! The deficit (net_global_qgwl) will be redistributed to top N outlets later. + ! --- SCENARIO B: Net global QGWL is negative. Redistribute ALL qgwl mass to outlets.--- + ! Set all local qgwl to zero - entire net amount will be added at outlets do nr = rtmCTL%begr, rtmCTL%endr TRunoff%qgwl(nr, :) = 0.0_r8 enddo @@ -2454,6 +2463,8 @@ subroutine Rtmrun(rstwr,nlend,rdate) enddo endif + ! Remove the incorrect "before" diagnostic from here - will be added at the right location + do nr = rtmCTL%begr,rtmCTL%endr do nt = 1,nt_rtm TRunoff%qsur(nr,nt) = rtmCTL%qsur(nr,nt) @@ -2942,12 +2953,29 @@ subroutine Rtmrun(rstwr,nlend,rdate) enddo enddo - ! Collect outlet discharge, Rank, Redistribute net_global_qgwl, Update discharge map --- + ! Collect outlet discharge, Rank, Redistribute qgwl, Update discharge map --- allocate(qgwl_correction_local(rtmCTL%begr:rtmCTL%endr)) ! Moved allocation here qgwl_correction_local = 0.0_r8 ! Initialize if (redirect_negative_qgwl_flag) then ! Check the main flag first - if (net_global_qgwl < 0.0_r8) then ! Only proceed if there's negative net qgwl to redistribute + ! Determine amount to redistribute: + ! Scenario A (net >= 0): redistribute abs(negative_sum) to balance negatives + ! Scenario B (net < 0): redistribute net_global_qgwl (entire net amount) + if (net_global_qgwl >= 0.0_r8) then + qgwl_to_redistribute = abs(global_negative_qgwl_sum) ! Scenario A + if (masterproc) then + write(iulog, *) trim(subname), 'Using Scenario A: redistributing negative qgwl (', & + qgwl_to_redistribute, ' m3/s) to outlets' + endif + else + qgwl_to_redistribute = net_global_qgwl ! Scenario B (negative value) + if (masterproc) then + write(iulog, *) trim(subname), 'Using Scenario B: redistributing entire net qgwl (', & + qgwl_to_redistribute, ' m3/s) to outlets' + endif + endif + + if (abs(qgwl_to_redistribute) > 1.0e-15_r8) then ! Only proceed if there's something to redistribute call t_startf('mosartr_qgwl_redir_dist') ! Identify local outlets and their current discharges @@ -3033,12 +3061,12 @@ subroutine Rtmrun(rstwr,nlend,rdate) if (abs(sum_discharge_top_n) > 1.0e-9_r8) then ! Avoid division by zero block real(r8) :: qgwl_to_discharge_ratio_percent - - qgwl_to_discharge_ratio_percent = (net_global_qgwl / sum_discharge_top_n) * 100.0_r8 - + + qgwl_to_discharge_ratio_percent = (qgwl_to_redistribute / sum_discharge_top_n) * 100.0_r8 + ! Print the ratio write(iulog, *) trim(subname), & - 'Debug QGWL Ratio: (GlobalNetQGWL / SumTopNDischarge) = ', & + 'Debug QGWL Ratio: (QGWL_to_redistribute / SumTopNDischarge) = ', & qgwl_to_discharge_ratio_percent, '%', & ' (N_used = ', current_top_n_count, & ', N_param = ', num_top_outlets_for_qgwl, ')' ! Using your var name @@ -3070,7 +3098,7 @@ subroutine Rtmrun(rstwr,nlend,rdate) qgwl_correction_values = 0.0_r8 do i = 1, current_top_n_count qgwl_correction_gindices(i) = sorted_outlets(i)%gidx - qgwl_correction_values(i) = (sorted_outlets(i)%discharge / sum_discharge_top_n) * net_global_qgwl + qgwl_correction_values(i) = (sorted_outlets(i)%discharge / sum_discharge_top_n) * qgwl_to_redistribute enddo else write(iulog, '(A)') trim(subname), & @@ -3081,12 +3109,12 @@ subroutine Rtmrun(rstwr,nlend,rdate) if (sum_discharge_top_n > 1.0e-9_r8) then ! Avoid division by zero do i = 1, current_top_n_count qgwl_correction_gindices(i) = sorted_outlets(i)%gidx - qgwl_correction_values(i) = (sorted_outlets(i)%discharge / sum_discharge_top_n) * net_global_qgwl + qgwl_correction_values(i) = (sorted_outlets(i)%discharge / sum_discharge_top_n) * qgwl_to_redistribute enddo else ! If sum of top N discharge is zero (or tiny), distribute equally (or handle as error/no distribution) do i = 1, current_top_n_count qgwl_correction_gindices(i) = sorted_outlets(i)%gidx - qgwl_correction_values(i) = net_global_qgwl / real(current_top_n_count, kind=r8) + qgwl_correction_values(i) = qgwl_to_redistribute / real(current_top_n_count, kind=r8) enddo endif endif @@ -3147,6 +3175,46 @@ subroutine Rtmrun(rstwr,nlend,rdate) elseif (rtmCTL%mask(nr) >= 2) then ! Ocean or Outlet rtmCTL%runoffocn(nr,nt) = rtmCTL%runoff(nr,nt) ! Routed component + ! ===================================================== + ! DIAGNOSTIC: Global sum BEFORE outlet corrections (once per timestep) + ! ===================================================== + if (nt == nt_nliq .and. nr == rtmCTL%begr) then ! Only do this once per timestep + block + real(r8) :: local_total_liquid_before, global_total_liquid_before + integer :: nr_temp + local_total_liquid_before = 0.0_r8 + do nr_temp = rtmCTL%begr, rtmCTL%endr + ! Sum what will be sent to coupler: rtmCTL%direct + rtmCTL%runoff + local_total_liquid_before = local_total_liquid_before + & + rtmCTL%direct(nr_temp, nt_nliq) + rtmCTL%runoff(nr_temp, nt_nliq) + enddo + + ! Use reproducible sum for bit-for-bit reproducibility across PE layouts + block + real(r8) :: before_local(1,1), before_global(1) + + ! Reproducible sum for total liquid before + before_local(1,1) = local_total_liquid_before + call shr_reprosum_calc(before_local, before_global, 1, 1, 1, & + commid=mpicom_rof) + global_total_liquid_before = before_global(1) + + ! Diagnostic output for debugging + if (masterproc) then + write(iulog, '(A,ES24.16)') trim(subname)//' REPROSUM Total Liquid Before = ', global_total_liquid_before + endif + end block + + if (masterproc) then + if (redirect_negative_qgwl_flag) then + write(iulog, '(A,ES24.16)') trim(subname)//' CONSERVATION CHECK: Global liquid runoff BEFORE outlet corrections (redirect=ON) = ', global_total_liquid_before + else + write(iulog, '(A,ES24.16)') trim(subname)//' CONSERVATION CHECK: Global liquid runoff BEFORE outlet corrections (baseline) = ', global_total_liquid_before + endif + endif + end block + endif + ! Apply the distributed global qgwl correction ONLY to liquid tracer at target outlets if (redirect_negative_qgwl_flag .and. nt == nt_nliq) then rtmCTL%runoffocn(nr,nt) = rtmCTL%runoffocn(nr,nt) + qgwl_correction_local(nr) @@ -3156,6 +3224,44 @@ subroutine Rtmrun(rstwr,nlend,rdate) endif enddo ! nr enddo ! nt + + ! ===================================================== + ! DIAGNOSTIC: Global sum of total liquid runoff AFTER redirection + ! This should equal the "before" sum for conservation (only when redirect is ON) + ! ===================================================== + if (redirect_negative_qgwl_flag) then + block + real(r8) :: local_total_liquid_after, global_total_liquid_after + local_total_liquid_after = 0.0_r8 + do nr = rtmCTL%begr, rtmCTL%endr + ! Sum what actually gets sent to coupler: rtmCTL%direct + rtmCTL%runoff + ! This matches exactly what rof_comp_mct.F90 sends to ocean + local_total_liquid_after = local_total_liquid_after + & + rtmCTL%direct(nr, nt_nliq) + rtmCTL%runoff(nr, nt_nliq) + enddo + + ! Use reproducible sum for bit-for-bit reproducibility across PE layouts + block + real(r8) :: after_local(1,1), after_global(1) + + ! Reproducible sum for total liquid after + after_local(1,1) = local_total_liquid_after + call shr_reprosum_calc(after_local, after_global, 1, 1, 1, & + commid=mpicom_rof) + global_total_liquid_after = after_global(1) + + ! Diagnostic output for debugging + if (masterproc) then + write(iulog, '(A,ES24.16)') trim(subname)//' REPROSUM Total Liquid After = ', global_total_liquid_after + endif + end block + + if (masterproc) then + write(iulog, '(A,ES24.16)') trim(subname)//' CONSERVATION CHECK: Global liquid runoff AFTER redirection = ', global_total_liquid_after + endif + end block + endif + if (allocated(qgwl_correction_local)) deallocate(qgwl_correction_local) From 330622278d51b82e6d04d0de1467073c49ad28f5 Mon Sep 17 00:00:00 2001 From: Tian Zhou Date: Sat, 18 Oct 2025 02:33:44 -0700 Subject: [PATCH 05/17] Fixing error for Scenario B --- components/mosart/src/riverroute/RtmMod.F90 | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/components/mosart/src/riverroute/RtmMod.F90 b/components/mosart/src/riverroute/RtmMod.F90 index 2758e406d097..28c7967323d6 100644 --- a/components/mosart/src/riverroute/RtmMod.F90 +++ b/components/mosart/src/riverroute/RtmMod.F90 @@ -2959,13 +2959,12 @@ subroutine Rtmrun(rstwr,nlend,rdate) if (redirect_negative_qgwl_flag) then ! Check the main flag first ! Determine amount to redistribute: - ! Scenario A (net >= 0): redistribute abs(negative_sum) to balance negatives - ! Scenario B (net < 0): redistribute net_global_qgwl (entire net amount) + ! Scenario A (net >= 0): Original qgwl flows to direct (no outlet correction needed) + ! Scenario B (net < 0): All qgwl zeroed, redistribute net amount to outlets if (net_global_qgwl >= 0.0_r8) then - qgwl_to_redistribute = abs(global_negative_qgwl_sum) ! Scenario A + qgwl_to_redistribute = 0.0_r8 ! Scenario A: negatives already in direct discharge if (masterproc) then - write(iulog, *) trim(subname), 'Using Scenario A: redistributing negative qgwl (', & - qgwl_to_redistribute, ' m3/s) to outlets' + write(iulog, *) trim(subname), 'Using Scenario A: original qgwl flows to direct (no outlet correction)' endif else qgwl_to_redistribute = net_global_qgwl ! Scenario B (negative value) @@ -3216,10 +3215,12 @@ subroutine Rtmrun(rstwr,nlend,rdate) endif ! Apply the distributed global qgwl correction ONLY to liquid tracer at target outlets + ! Apply to rtmCTL%runoff (not runoffocn) since runoff is what gets sent to coupler if (redirect_negative_qgwl_flag .and. nt == nt_nliq) then - rtmCTL%runoffocn(nr,nt) = rtmCTL%runoffocn(nr,nt) + qgwl_correction_local(nr) + rtmCTL%runoff(nr,nt) = rtmCTL%runoff(nr,nt) + qgwl_correction_local(nr) + rtmCTL%runoffocn(nr,nt) = rtmCTL%runoff(nr,nt) ! Keep runoffocn in sync endif - rtmCTL%runofftot(nr,nt) = rtmCTL%runofftot(nr,nt) + rtmCTL%runoffocn(nr,nt) ! Add potentially corrected routed flow + rtmCTL%runofftot(nr,nt) = rtmCTL%runofftot(nr,nt) + rtmCTL%runoffocn(nr,nt) ! Add corrected routed flow rtmCTL%dvolrdtocn(nr,nt)= rtmCTL%dvolrdt(nr,nt) endif enddo ! nr @@ -3235,7 +3236,7 @@ subroutine Rtmrun(rstwr,nlend,rdate) local_total_liquid_after = 0.0_r8 do nr = rtmCTL%begr, rtmCTL%endr ! Sum what actually gets sent to coupler: rtmCTL%direct + rtmCTL%runoff - ! This matches exactly what rof_comp_mct.F90 sends to ocean + ! rtmCTL%runoff includes the qgwl correction applied at line 3221 local_total_liquid_after = local_total_liquid_after + & rtmCTL%direct(nr, nt_nliq) + rtmCTL%runoff(nr, nt_nliq) enddo @@ -3301,7 +3302,7 @@ subroutine Rtmrun(rstwr,nlend,rdate) if (rtmCTL%mask(nr) >= 2) then ! (2 -- Ocean; 3 -- Outlet. --Inund.) - budget_terms(br_ocnout,nt) = budget_terms(br_ocnout,nt) + rtmCTL%runoff(nr,nt)*delt_coupling ! (Volume of outflows to oceans. Note that rtmCTL%runoff is averge value of sub-steps used by MOSART. --Inund.) + budget_terms(br_ocnout,nt) = budget_terms(br_ocnout,nt) + rtmCTL%runoff(nr,nt)*delt_coupling ! (Volume of outflows to oceans. runoff includes qgwl corrections from line 3221. --Inund.) budget_terms(br_erolpo,nt) = budget_terms(br_erolpo,nt) + eroup_lagi(nr,nt)*delt_coupling ! (Volume of outflows to oceans in the previous MOSART sub-step. ! Also eroup_lagi is averge value of several previous MOSART sub-steps. --Inund.) budget_terms(br_erolco,nt) = budget_terms(br_erolco,nt) + eroup_lagf(nr,nt)*delt_coupling ! (Volume of outflows to oceans. Actually eroup_lagf is same as the above rtmCTL%runoff. --Inund.) From 9078cbad38b3590357f06c8abb5d61a5fdfd0517 Mon Sep 17 00:00:00 2001 From: Tian Zhou Date: Sat, 18 Oct 2025 03:22:13 -0700 Subject: [PATCH 06/17] Fixing error for Scenario A --- components/mosart/src/riverroute/RtmMod.F90 | 26 ++++++++++++++++----- 1 file changed, 20 insertions(+), 6 deletions(-) diff --git a/components/mosart/src/riverroute/RtmMod.F90 b/components/mosart/src/riverroute/RtmMod.F90 index 28c7967323d6..57a3c6cdc874 100644 --- a/components/mosart/src/riverroute/RtmMod.F90 +++ b/components/mosart/src/riverroute/RtmMod.F90 @@ -2242,7 +2242,7 @@ subroutine Rtmrun(rstwr,nlend,rdate) ! Variables for negative runoff redirection real(r8) :: local_positive_qgwl_sum, local_negative_qgwl_sum - real(r8) :: net_global_qgwl, original_cell_qgwl, reduction + real(r8) :: net_global_qgwl, original_cell_qgwl, reduction, scaling_factor integer, allocatable :: outlet_gindices_local(:) ! Local array of global indices of outlets on this task real(r8), allocatable :: outlet_discharges_local(:) ! Local array of discharges for these outlets @@ -2440,14 +2440,28 @@ subroutine Rtmrun(rstwr,nlend,rdate) endif ! Decide how to set TRunoff%qgwl for local routing - ! For BOTH scenarios, we let original qgwl flow to direct discharge, - ! and then add corrections at outlets to balance negative values if (net_global_qgwl >= 0.0_r8) then ! --- SCENARIO A: Net global QGWL is non-negative. --- - ! Let original qgwl flow to rtmCTL%direct (no modification here) - ! The negative qgwl will be balanced by adding to outlets later + ! Proportionally reduce positive qgwl to offset negatives, zero out negative cells + ! This ensures NO negative qgwl enters the ocean at any grid cell + + if (global_positive_qgwl_sum > 0.0_r8) then + ! Calculate scaling factor: (positive - |negative|) / positive = net / positive + scaling_factor = net_global_qgwl / global_positive_qgwl_sum + else + scaling_factor = 1.0_r8 ! No positive qgwl, keep as is + endif + do nr = rtmCTL%begr, rtmCTL%endr - TRunoff%qgwl(nr, :) = rtmCTL%qgwl(nr, :) + do nt = 1, nt_rtm + if (rtmCTL%qgwl(nr, nt) > 0.0_r8) then + ! Positive cell: scale down proportionally + TRunoff%qgwl(nr, nt) = rtmCTL%qgwl(nr, nt) * scaling_factor + else + ! Negative or zero cell: set to zero + TRunoff%qgwl(nr, nt) = 0.0_r8 + endif + enddo enddo else ! --- SCENARIO B: Net global QGWL is negative. Redistribute ALL qgwl mass to outlets.--- From 2fc7068e79704e256d7f43ba6e2bb58bab11eb45 Mon Sep 17 00:00:00 2001 From: Tian Zhou Date: Sun, 19 Oct 2025 01:49:42 -0700 Subject: [PATCH 07/17] optimize diagnostic outputs --- components/mosart/src/riverroute/RtmMod.F90 | 332 +++++--------------- 1 file changed, 84 insertions(+), 248 deletions(-) diff --git a/components/mosart/src/riverroute/RtmMod.F90 b/components/mosart/src/riverroute/RtmMod.F90 index 57a3c6cdc874..0dc179c42c24 100644 --- a/components/mosart/src/riverroute/RtmMod.F90 +++ b/components/mosart/src/riverroute/RtmMod.F90 @@ -2257,6 +2257,13 @@ subroutine Rtmrun(rstwr,nlend,rdate) integer, allocatable :: qgwl_correction_gindices(:) ! Global indices for these corrections real(r8), allocatable :: qgwl_correction_local(:) ! Local portion of correction array real(r8) :: qgwl_to_redistribute ! Amount to redistribute (Scenario A or B) + real(r8) :: local_total_outlet_discharge, global_total_outlet_discharge + real(r8) :: outlet_discharge_local(1,1), outlet_discharge_global(1) + real(r8) :: qgwl_to_discharge_ratio_percent + real(r8) :: conservation_error + real(r8) :: global_total_liquid_before, global_total_liquid_after + real(r8) :: local_correction_sum, global_correction_sum + real(r8) :: correction_local(1,1), correction_global(1) character(len=*),parameter :: subname = '(Rtmrun) ' !----------------------------------------------------------------------- @@ -2423,8 +2430,8 @@ subroutine Rtmrun(rstwr,nlend,rdate) commid=mpicom_rof) global_negative_qgwl_sum = neg_global(1) - ! Diagnostic output for debugging - if (masterproc) then + ! Diagnostic output only when do_budget == 3 + if (masterproc .and. do_budget == 3) then write(iulog, '(A,ES24.16)') trim(subname)//' REPROSUM Positive Sum = ', global_positive_qgwl_sum write(iulog, '(A,ES24.16)') trim(subname)//' REPROSUM Negative Sum = ', global_negative_qgwl_sum write(iulog, '(A,ES24.16)') trim(subname)//' REPROSUM Net Sum = ', global_positive_qgwl_sum + global_negative_qgwl_sum @@ -2433,7 +2440,7 @@ subroutine Rtmrun(rstwr,nlend,rdate) net_global_qgwl = global_positive_qgwl_sum + global_negative_qgwl_sum - if (masterproc) then + if (masterproc .and. do_budget == 3) then write(iulog, *) trim(subname), 'Debug QGWL Balancing: Global Positive Sum =', global_positive_qgwl_sum write(iulog, *) trim(subname), 'Debug QGWL Balancing: Global Negative Sum =', global_negative_qgwl_sum write(iulog, *) trim(subname), 'Debug QGWL Balancing: Global Net Sum =', net_global_qgwl @@ -2977,12 +2984,12 @@ subroutine Rtmrun(rstwr,nlend,rdate) ! Scenario B (net < 0): All qgwl zeroed, redistribute net amount to outlets if (net_global_qgwl >= 0.0_r8) then qgwl_to_redistribute = 0.0_r8 ! Scenario A: negatives already in direct discharge - if (masterproc) then + if (masterproc .and. do_budget == 3) then write(iulog, *) trim(subname), 'Using Scenario A: original qgwl flows to direct (no outlet correction)' endif else qgwl_to_redistribute = net_global_qgwl ! Scenario B (negative value) - if (masterproc) then + if (masterproc .and. do_budget == 3) then write(iulog, *) trim(subname), 'Using Scenario B: redistributing entire net qgwl (', & qgwl_to_redistribute, ' m3/s) to outlets' endif @@ -2991,190 +2998,68 @@ subroutine Rtmrun(rstwr,nlend,rdate) if (abs(qgwl_to_redistribute) > 1.0e-15_r8) then ! Only proceed if there's something to redistribute call t_startf('mosartr_qgwl_redir_dist') - ! Identify local outlets and their current discharges - local_outlet_count = 0 - do nr = rtmCTL%begr, rtmCTL%endr - if (rtmCTL%mask(nr) == 3) then ! Is it an outlet? - local_outlet_count = local_outlet_count + 1 - endif - enddo - - if (allocated(outlet_gindices_local)) deallocate(outlet_gindices_local) - if (allocated(outlet_discharges_local)) deallocate(outlet_discharges_local) - allocate(outlet_gindices_local(local_outlet_count)) - allocate(outlet_discharges_local(local_outlet_count)) - local_outlet_count = 0 ! Reset for population + ! Calculate total discharge from all outlets globally using reprosum + local_total_outlet_discharge = 0.0_r8 do nr = rtmCTL%begr, rtmCTL%endr - if (rtmCTL%mask(nr) == 3) then - local_outlet_count = local_outlet_count + 1 - outlet_gindices_local(local_outlet_count) = rtmCTL%gindex(nr) - outlet_discharges_local(local_outlet_count) = rtmCTL%runoffocn(nr, nt_nliq) ! Current discharge - endif + if (rtmCTL%mask(nr) == 3) then ! Outlet cell + local_total_outlet_discharge = local_total_outlet_discharge + rtmCTL%runoffocn(nr, nt_nliq) + endif enddo - ! Gather all outlet discharges to master for ranking - if (allocated(recvcounts)) deallocate(recvcounts) - if (allocated(displs)) deallocate(displs) - allocate(recvcounts(npes), displs(npes)) - call MPI_Gather(local_outlet_count, 1, MPI_INTEGER, recvcounts, 1, MPI_INTEGER, mastertask, mpicom_rof, ier) + ! Use reproducible sum for bit-for-bit reproducibility across PE layouts + outlet_discharge_local(1,1) = local_total_outlet_discharge + call shr_reprosum_calc(outlet_discharge_local, outlet_discharge_global, 1, 1, 1, commid=mpicom_rof) + global_total_outlet_discharge = outlet_discharge_global(1) + ! Print diagnostic on master if (masterproc) then - displs(1) = 0 - num_all_outlets_global = recvcounts(1) - do i = 2, npes - displs(i) = displs(i-1) + recvcounts(i-1) - num_all_outlets_global = num_all_outlets_global + recvcounts(i) - enddo - if (allocated(all_outlet_gindices)) deallocate(all_outlet_gindices) - if (allocated(all_outlet_discharges)) deallocate(all_outlet_discharges) - if (num_all_outlets_global > 0) then - allocate(all_outlet_gindices(num_all_outlets_global)) - allocate(all_outlet_discharges(num_all_outlets_global)) - else - ! Handle cases where no outlets are found if necessary, or let current_top_n_count be 0 - num_all_outlets_global = 0 ! Ensure it's zero if no outlets - endif - else - num_all_outlets_global = 0 ! Non-master tasks initialize to 0 - endif - call MPI_Bcast(num_all_outlets_global, 1, MPI_INTEGER, mastertask, mpicom_rof, ier) - - if(num_all_outlets_global > 0) then - call MPI_Gatherv(outlet_gindices_local, local_outlet_count, MPI_INTEGER, & - all_outlet_gindices, recvcounts, displs, MPI_INTEGER, & - mastertask, mpicom_rof, ier) - call MPI_Gatherv(outlet_discharges_local, local_outlet_count, MPI_REAL8, & - all_outlet_discharges, recvcounts, displs, MPI_REAL8, & - mastertask, mpicom_rof, ier) - endif + if (abs(global_total_outlet_discharge) > 1.0e-9_r8) then + qgwl_to_discharge_ratio_percent = (qgwl_to_redistribute / global_total_outlet_discharge) * 100.0_r8 - if (allocated(outlet_gindices_local)) deallocate(outlet_gindices_local) - if (allocated(outlet_discharges_local)) deallocate(outlet_discharges_local) + if (do_budget == 3) then + write(iulog, *) trim(subname), & + 'Debug QGWL Ratio: (QGWL_to_redistribute / TotalOutletDischarge) = ', & + qgwl_to_discharge_ratio_percent, '% (', qgwl_to_redistribute, ' m3/s)' + endif - ! Rank on master, calculate proportions, and prepare corrections - current_top_n_count = 0 ! Initialize for all tasks - if (masterproc) then - if (num_all_outlets_global > 0) then - allocate(sorted_outlets(num_all_outlets_global)) - do i = 1, num_all_outlets_global - sorted_outlets(i)%gidx = all_outlet_gindices(i) - sorted_outlets(i)%discharge = all_outlet_discharges(i) - enddo - if (allocated(all_outlet_gindices)) deallocate(all_outlet_gindices) - if (allocated(all_outlet_discharges)) deallocate(all_outlet_discharges) - - call sort_outlets_by_discharge_desc(sorted_outlets, num_all_outlets_global) - - current_top_n_count = min(num_top_outlets_for_qgwl, num_all_outlets_global) - sum_discharge_top_n = 0.0_r8 - do i = 1, current_top_n_count - sum_discharge_top_n = sum_discharge_top_n + sorted_outlets(i)%discharge - enddo - if(current_top_n_count > 0) then - if (abs(sum_discharge_top_n) > 1.0e-9_r8) then ! Avoid division by zero - block - real(r8) :: qgwl_to_discharge_ratio_percent - - qgwl_to_discharge_ratio_percent = (qgwl_to_redistribute / sum_discharge_top_n) * 100.0_r8 - - ! Print the ratio - write(iulog, *) trim(subname), & - 'Debug QGWL Ratio: (QGWL_to_redistribute / SumTopNDischarge) = ', & - qgwl_to_discharge_ratio_percent, '%', & - ' (N_used = ', current_top_n_count, & - ', N_param = ', num_top_outlets_for_qgwl, ')' ! Using your var name - - ! Check and print warning if magnitude is > 5% - if (abs(qgwl_to_discharge_ratio_percent) > 5.0_r8) then - call shr_sys_flush(iulog) ! Flush previous message - write(iulog, *) trim(subname), & - 'WARNING: QGWL_Redist_Ratio magnitude > 5% (', & - qgwl_to_discharge_ratio_percent, '%). ', & - 'N_used = ', current_top_n_count, & - ', N_param = ', num_top_outlets_for_qgwl, & - '. Consider increasing N_param.' - call shr_sys_flush(iulog) ! Flush warning - endif - end block - else - write(iulog, *) trim(subname), & - 'Debug QGWL Ratio: Sum of Top ', current_top_n_count, & - ' Outlet Discharges is near zero. Ratio not calculated.' - endif - else + ! Always check and warn if magnitude is > 5% (regardless of do_budget) + if (abs(qgwl_to_discharge_ratio_percent) > 5.0_r8) then + call shr_sys_flush(iulog) write(iulog, *) trim(subname), & - 'Debug QGWL Ratio: No top outlets selected (N_used = 0). Redistribution skipped.' - endif - if (current_top_n_count > 0) then - allocate(qgwl_correction_gindices(current_top_n_count)) - allocate(qgwl_correction_values(current_top_n_count)) - qgwl_correction_values = 0.0_r8 - do i = 1, current_top_n_count - qgwl_correction_gindices(i) = sorted_outlets(i)%gidx - qgwl_correction_values(i) = (sorted_outlets(i)%discharge / sum_discharge_top_n) * qgwl_to_redistribute - enddo - else - write(iulog, '(A)') trim(subname), & - 'Debug: No top outlets selected for QGWL redistribution (current_top_n_count is 0).' - allocate(qgwl_correction_gindices(current_top_n_count)) - allocate(qgwl_correction_values(current_top_n_count)) - qgwl_correction_values = 0.0_r8 - if (sum_discharge_top_n > 1.0e-9_r8) then ! Avoid division by zero - do i = 1, current_top_n_count - qgwl_correction_gindices(i) = sorted_outlets(i)%gidx - qgwl_correction_values(i) = (sorted_outlets(i)%discharge / sum_discharge_top_n) * qgwl_to_redistribute - enddo - else ! If sum of top N discharge is zero (or tiny), distribute equally (or handle as error/no distribution) - do i = 1, current_top_n_count - qgwl_correction_gindices(i) = sorted_outlets(i)%gidx - qgwl_correction_values(i) = qgwl_to_redistribute / real(current_top_n_count, kind=r8) - enddo - endif - endif - if(allocated(sorted_outlets)) deallocate(sorted_outlets) - else - current_top_n_count = 0 ! No outlets globally - endif ! num_all_outlets_global > 0 - endif ! masterproc - - ! Broadcast correction info - call MPI_Bcast(current_top_n_count, 1, MPI_INTEGER, mastertask, mpicom_rof, ier) - if (.not. masterproc .and. current_top_n_count > 0) then - if (allocated(qgwl_correction_gindices)) deallocate(qgwl_correction_gindices) - if (allocated(qgwl_correction_values)) deallocate(qgwl_correction_values) - allocate(qgwl_correction_gindices(current_top_n_count)) - allocate(qgwl_correction_values(current_top_n_count)) + 'WARNING: QGWL_Redist_Ratio magnitude > 5% (', & + qgwl_to_discharge_ratio_percent, '%). This may impact global hydrology.' + call shr_sys_flush(iulog) + endif + else + if (do_budget == 3) then + write(iulog, *) trim(subname), & + 'Debug QGWL Ratio: Total outlet discharge is near zero. Ratio not calculated.' + endif + endif endif - if (current_top_n_count > 0) then - call MPI_Bcast(qgwl_correction_gindices, current_top_n_count, MPI_INTEGER, mastertask, mpicom_rof, ier) - call MPI_Bcast(qgwl_correction_values, current_top_n_count, MPI_REAL8, mastertask, mpicom_rof, ier) - - ! Populate local correction array - do k = 1, current_top_n_count - target_gidx = qgwl_correction_gindices(k) - do nr_local = rtmCTL%begr, rtmCTL%endr - if (rtmCTL%gindex(nr_local) == target_gidx) then - qgwl_correction_local(nr_local) = qgwl_correction_local(nr_local) + qgwl_correction_values(k) - exit - endif - enddo - enddo - endif ! current_top_n_count > 0 for broadcast and local application - if (allocated(qgwl_correction_gindices)) deallocate(qgwl_correction_gindices) - if (allocated(qgwl_correction_values)) deallocate(qgwl_correction_values) - if (allocated(recvcounts)) deallocate(recvcounts) - if (allocated(displs)) deallocate(displs) + ! Each PE calculates corrections for its local outlets proportionally + if (abs(global_total_outlet_discharge) > 1.0e-9_r8) then + do nr = rtmCTL%begr, rtmCTL%endr + if (rtmCTL%mask(nr) == 3) then ! Outlet cell + ! Proportional correction based on this outlet's discharge + qgwl_correction_local(nr) = (rtmCTL%runoffocn(nr, nt_nliq) / global_total_outlet_discharge) * qgwl_to_redistribute + endif + enddo + else + ! If total discharge is zero, no redistribution (corrections remain 0) + if (masterproc) then + write(iulog, *) trim(subname), & + 'Warning: Cannot redistribute qgwl - total outlet discharge is zero' + endif + endif call t_stopf('mosartr_qgwl_redir_dist') - else - write(iulog, *) trim(subname), & - 'Global net QGWL is positive, offsetting the negative values.' - endif ! net_global_qgwl < 0.0_r8 - endif ! redirect_negative_qgwl_flag + endif ! abs(qgwl_to_redistribute) > 1.0e-15 + endif ! redirect_negative_qgwl_flag - ! This loop should now use the qgwl_correction_local array + ! Apply outlet corrections to runoff do nt = 1,nt_rtm do nr = rtmCTL%begr,rtmCTL%endr volr_init = rtmCTL%volr(nr,nt) @@ -3188,46 +3073,6 @@ subroutine Rtmrun(rstwr,nlend,rdate) elseif (rtmCTL%mask(nr) >= 2) then ! Ocean or Outlet rtmCTL%runoffocn(nr,nt) = rtmCTL%runoff(nr,nt) ! Routed component - ! ===================================================== - ! DIAGNOSTIC: Global sum BEFORE outlet corrections (once per timestep) - ! ===================================================== - if (nt == nt_nliq .and. nr == rtmCTL%begr) then ! Only do this once per timestep - block - real(r8) :: local_total_liquid_before, global_total_liquid_before - integer :: nr_temp - local_total_liquid_before = 0.0_r8 - do nr_temp = rtmCTL%begr, rtmCTL%endr - ! Sum what will be sent to coupler: rtmCTL%direct + rtmCTL%runoff - local_total_liquid_before = local_total_liquid_before + & - rtmCTL%direct(nr_temp, nt_nliq) + rtmCTL%runoff(nr_temp, nt_nliq) - enddo - - ! Use reproducible sum for bit-for-bit reproducibility across PE layouts - block - real(r8) :: before_local(1,1), before_global(1) - - ! Reproducible sum for total liquid before - before_local(1,1) = local_total_liquid_before - call shr_reprosum_calc(before_local, before_global, 1, 1, 1, & - commid=mpicom_rof) - global_total_liquid_before = before_global(1) - - ! Diagnostic output for debugging - if (masterproc) then - write(iulog, '(A,ES24.16)') trim(subname)//' REPROSUM Total Liquid Before = ', global_total_liquid_before - endif - end block - - if (masterproc) then - if (redirect_negative_qgwl_flag) then - write(iulog, '(A,ES24.16)') trim(subname)//' CONSERVATION CHECK: Global liquid runoff BEFORE outlet corrections (redirect=ON) = ', global_total_liquid_before - else - write(iulog, '(A,ES24.16)') trim(subname)//' CONSERVATION CHECK: Global liquid runoff BEFORE outlet corrections (baseline) = ', global_total_liquid_before - endif - endif - end block - endif - ! Apply the distributed global qgwl correction ONLY to liquid tracer at target outlets ! Apply to rtmCTL%runoff (not runoffocn) since runoff is what gets sent to coupler if (redirect_negative_qgwl_flag .and. nt == nt_nliq) then @@ -3240,43 +3085,34 @@ subroutine Rtmrun(rstwr,nlend,rdate) enddo ! nr enddo ! nt - ! ===================================================== - ! DIAGNOSTIC: Global sum of total liquid runoff AFTER redirection - ! This should equal the "before" sum for conservation (only when redirect is ON) - ! ===================================================== - if (redirect_negative_qgwl_flag) then - block - real(r8) :: local_total_liquid_after, global_total_liquid_after - local_total_liquid_after = 0.0_r8 - do nr = rtmCTL%begr, rtmCTL%endr - ! Sum what actually gets sent to coupler: rtmCTL%direct + rtmCTL%runoff - ! rtmCTL%runoff includes the qgwl correction applied at line 3221 - local_total_liquid_after = local_total_liquid_after + & - rtmCTL%direct(nr, nt_nliq) + rtmCTL%runoff(nr, nt_nliq) - enddo - - ! Use reproducible sum for bit-for-bit reproducibility across PE layouts - block - real(r8) :: after_local(1,1), after_global(1) - - ! Reproducible sum for total liquid after - after_local(1,1) = local_total_liquid_after - call shr_reprosum_calc(after_local, after_global, 1, 1, 1, & - commid=mpicom_rof) - global_total_liquid_after = after_global(1) - - ! Diagnostic output for debugging - if (masterproc) then - write(iulog, '(A,ES24.16)') trim(subname)//' REPROSUM Total Liquid After = ', global_total_liquid_after - endif - end block + ! Conservation check: verify that corrections were applied correctly + if (redirect_negative_qgwl_flag .and. abs(qgwl_to_redistribute) > 1.0e-15_r8) then + ! Sum up all the corrections that were applied + local_correction_sum = 0.0_r8 + do nr = rtmCTL%begr, rtmCTL%endr + local_correction_sum = local_correction_sum + qgwl_correction_local(nr) + enddo + + ! Use reprosum for bit-for-bit reproducibility + correction_local(1,1) = local_correction_sum + call shr_reprosum_calc(correction_local, correction_global, 1, 1, 1, commid=mpicom_rof) + global_correction_sum = correction_global(1) - if (masterproc) then - write(iulog, '(A,ES24.16)') trim(subname)//' CONSERVATION CHECK: Global liquid runoff AFTER redirection = ', global_total_liquid_after + ! Check if total corrections match what we intended to redistribute + if (masterproc) then + conservation_error = abs(global_correction_sum - qgwl_to_redistribute) + if (conservation_error > 1.0e-10_r8) then + call shr_sys_flush(iulog) + write(iulog, '(A)') trim(subname)//' WARNING: QGWL redistribution error detected!' + write(iulog, '(A,ES24.16)') trim(subname)//' Intended to redistribute = ', qgwl_to_redistribute + write(iulog, '(A,ES24.16)') trim(subname)//' Actually redistributed = ', global_correction_sum + write(iulog, '(A,ES24.16)') trim(subname)//' Difference = ', conservation_error + call shr_sys_flush(iulog) endif - end block + endif endif + ! Deallocate qgwl_correction_local after use if (allocated(qgwl_correction_local)) deallocate(qgwl_correction_local) From df6114fb6932b03acc1781e971e7786877477371 Mon Sep 17 00:00:00 2001 From: Tian Zhou Date: Mon, 27 Oct 2025 15:46:12 -0700 Subject: [PATCH 08/17] try to fix PEM test fail --- components/mosart/src/riverroute/RtmMod.F90 | 98 +++++++++------------ 1 file changed, 40 insertions(+), 58 deletions(-) diff --git a/components/mosart/src/riverroute/RtmMod.F90 b/components/mosart/src/riverroute/RtmMod.F90 index 0dc179c42c24..3e3fd68a9291 100644 --- a/components/mosart/src/riverroute/RtmMod.F90 +++ b/components/mosart/src/riverroute/RtmMod.F90 @@ -132,11 +132,6 @@ module RtmMod logical :: do_rtmflood logical :: do_rtm - type :: outlet_discharge_info_type - integer :: gidx - real(r8) :: discharge - end type outlet_discharge_info_type - ! Namelist variable and flag for TNR redirection logical, public :: redirect_negative_qgwl = .false. ! Namelist control logical, save :: redirect_negative_qgwl_flag = .false. ! Module flag @@ -2244,17 +2239,6 @@ subroutine Rtmrun(rstwr,nlend,rdate) real(r8) :: local_positive_qgwl_sum, local_negative_qgwl_sum real(r8) :: net_global_qgwl, original_cell_qgwl, reduction, scaling_factor - integer, allocatable :: outlet_gindices_local(:) ! Local array of global indices of outlets on this task - real(r8), allocatable :: outlet_discharges_local(:) ! Local array of discharges for these outlets - integer :: local_outlet_count - integer, allocatable :: all_outlet_gindices(:) ! Gathered on master - real(r8), allocatable :: all_outlet_discharges(:) ! Gathered on master - integer, allocatable :: recvcounts(:), displs(:) ! For MPI_Gatherv - type(outlet_discharge_info_type), allocatable :: sorted_outlets(:) - integer :: num_all_outlets_global, current_top_n_count, target_gidx - real(r8) :: sum_discharge_top_n - real(r8), allocatable :: qgwl_correction_values(:) ! Values to apply - integer, allocatable :: qgwl_correction_gindices(:) ! Global indices for these corrections real(r8), allocatable :: qgwl_correction_local(:) ! Local portion of correction array real(r8) :: qgwl_to_redistribute ! Amount to redistribute (Scenario A or B) real(r8) :: local_total_outlet_discharge, global_total_outlet_discharge @@ -2264,6 +2248,7 @@ subroutine Rtmrun(rstwr,nlend,rdate) real(r8) :: global_total_liquid_before, global_total_liquid_after real(r8) :: local_correction_sum, global_correction_sum real(r8) :: correction_local(1,1), correction_global(1) + real(r8) :: correction_ratio ! Ratio to apply: qgwl_to_redistribute / global_total_outlet_discharge character(len=*),parameter :: subname = '(Rtmrun) ' !----------------------------------------------------------------------- @@ -2414,21 +2399,21 @@ subroutine Rtmrun(rstwr,nlend,rdate) enddo ! Use reproducible sum for bit-for-bit reproducibility across PE layouts + ! Combine both positive and negative sums in one call for efficiency block - real(r8) :: pos_local(1,1), pos_global(1) - real(r8) :: neg_local(1,1), neg_global(1) + real(r8) :: local_sums(2,1), global_sums(2) - ! Reproducible sum for positive qgwl - pos_local(1,1) = local_positive_qgwl_sum - call shr_reprosum_calc(pos_local, pos_global, 1, 1, 1, & - commid=mpicom_rof) - global_positive_qgwl_sum = pos_global(1) + ! Pack both positive and negative sums + local_sums(1,1) = local_positive_qgwl_sum + local_sums(2,1) = local_negative_qgwl_sum - ! Reproducible sum for negative qgwl - neg_local(1,1) = local_negative_qgwl_sum - call shr_reprosum_calc(neg_local, neg_global, 1, 1, 1, & + ! Single reproducible sum call for both fields + call shr_reprosum_calc(local_sums, global_sums, 2, 1, 1, & commid=mpicom_rof) - global_negative_qgwl_sum = neg_global(1) + + ! Unpack results + global_positive_qgwl_sum = global_sums(1) + global_negative_qgwl_sum = global_sums(2) ! Diagnostic output only when do_budget == 3 if (masterproc .and. do_budget == 3) then @@ -2452,13 +2437,19 @@ subroutine Rtmrun(rstwr,nlend,rdate) ! Proportionally reduce positive qgwl to offset negatives, zero out negative cells ! This ensures NO negative qgwl enters the ocean at any grid cell - if (global_positive_qgwl_sum > 0.0_r8) then - ! Calculate scaling factor: (positive - |negative|) / positive = net / positive - scaling_factor = net_global_qgwl / global_positive_qgwl_sum - else - scaling_factor = 1.0_r8 ! No positive qgwl, keep as is + ! Master PE calculates scaling factor, then broadcasts for bit-for-bit reproducibility + if (masterproc) then + if (global_positive_qgwl_sum > 0.0_r8) then + ! Calculate scaling factor: (positive - |negative|) / positive = net / positive + scaling_factor = net_global_qgwl / global_positive_qgwl_sum + else + scaling_factor = 1.0_r8 ! No positive qgwl, keep as is + endif endif + ! Broadcast scaling factor from master to all PEs + call MPI_Bcast(scaling_factor, 1, MPI_REAL8, 0, mpicom_rof, ier) + do nr = rtmCTL%begr, rtmCTL%endr do nt = 1, nt_rtm if (rtmCTL%qgwl(nr, nt) > 0.0_r8) then @@ -2998,10 +2989,11 @@ subroutine Rtmrun(rstwr,nlend,rdate) if (abs(qgwl_to_redistribute) > 1.0e-15_r8) then ! Only proceed if there's something to redistribute call t_startf('mosartr_qgwl_redir_dist') - ! Calculate total discharge from all outlets globally using reprosum + ! Calculate total discharge from outlets with discharge > 1.0 m3/s + ! This threshold filters out small outlets and ensures numerical stability local_total_outlet_discharge = 0.0_r8 do nr = rtmCTL%begr, rtmCTL%endr - if (rtmCTL%mask(nr) == 3) then ! Outlet cell + if (rtmCTL%mask(nr) == 3 .and. rtmCTL%runoffocn(nr, nt_nliq) > 1.0_r8) then local_total_outlet_discharge = local_total_outlet_discharge + rtmCTL%runoffocn(nr, nt_nliq) endif enddo @@ -3038,19 +3030,28 @@ subroutine Rtmrun(rstwr,nlend,rdate) endif endif - ! Each PE calculates corrections for its local outlets proportionally + ! Master PE calculates the correction ratio, then broadcast to all PEs + ! This ensures bit-for-bit reproducibility across different PE layouts (PEM test) if (abs(global_total_outlet_discharge) > 1.0e-9_r8) then + if (masterproc) then + correction_ratio = qgwl_to_redistribute / global_total_outlet_discharge + endif + + ! Broadcast the correction ratio from master to all PEs + call MPI_Bcast(correction_ratio, 1, MPI_REAL8, 0, mpicom_rof, ier) + + ! All PEs apply the same ratio to their local outlets > 1.0 m3/s do nr = rtmCTL%begr, rtmCTL%endr - if (rtmCTL%mask(nr) == 3) then ! Outlet cell - ! Proportional correction based on this outlet's discharge - qgwl_correction_local(nr) = (rtmCTL%runoffocn(nr, nt_nliq) / global_total_outlet_discharge) * qgwl_to_redistribute + if (rtmCTL%mask(nr) == 3 .and. rtmCTL%runoffocn(nr, nt_nliq) > 1.0_r8) then + ! Apply the correction ratio + qgwl_correction_local(nr) = rtmCTL%runoffocn(nr, nt_nliq) * correction_ratio endif enddo else ! If total discharge is zero, no redistribution (corrections remain 0) if (masterproc) then write(iulog, *) trim(subname), & - 'Warning: Cannot redistribute qgwl - total outlet discharge is zero' + 'Warning: Cannot redistribute qgwl - total outlet discharge (>1 m3/s) is zero' endif endif @@ -5348,23 +5349,4 @@ end subroutine SubTimestep !----------------------------------------------------------------------- -subroutine sort_outlets_by_discharge_desc(outlets_array, count) !TZ negative runoff - implicit none - integer, intent(in) :: count - type(outlet_discharge_info_type), intent(inout) :: outlets_array(:) ! Assumed-shape - integer :: i, j - type(outlet_discharge_info_type) :: temp_outlet_info - ! Bubble sort - if (count < 2) return - do i = 1, count - 1 - do j = 1, count - i - if (outlets_array(j)%discharge < outlets_array(j+1)%discharge) then - temp_outlet_info = outlets_array(j) - outlets_array(j) = outlets_array(j+1) - outlets_array(j+1) = temp_outlet_info - endif - enddo - enddo -end subroutine sort_outlets_by_discharge_desc - end module RtmMod \ No newline at end of file From 8948e32f838bffff5337f5978c5ad7ab3c86bfda Mon Sep 17 00:00:00 2001 From: Tian Zhou Date: Thu, 30 Oct 2025 11:58:14 -0700 Subject: [PATCH 09/17] Revert "try to fix PEM test fail" This reverts commit df6114fb6932b03acc1781e971e7786877477371. --- components/mosart/src/riverroute/RtmMod.F90 | 98 ++++++++++++--------- 1 file changed, 58 insertions(+), 40 deletions(-) diff --git a/components/mosart/src/riverroute/RtmMod.F90 b/components/mosart/src/riverroute/RtmMod.F90 index 3e3fd68a9291..0dc179c42c24 100644 --- a/components/mosart/src/riverroute/RtmMod.F90 +++ b/components/mosart/src/riverroute/RtmMod.F90 @@ -132,6 +132,11 @@ module RtmMod logical :: do_rtmflood logical :: do_rtm + type :: outlet_discharge_info_type + integer :: gidx + real(r8) :: discharge + end type outlet_discharge_info_type + ! Namelist variable and flag for TNR redirection logical, public :: redirect_negative_qgwl = .false. ! Namelist control logical, save :: redirect_negative_qgwl_flag = .false. ! Module flag @@ -2239,6 +2244,17 @@ subroutine Rtmrun(rstwr,nlend,rdate) real(r8) :: local_positive_qgwl_sum, local_negative_qgwl_sum real(r8) :: net_global_qgwl, original_cell_qgwl, reduction, scaling_factor + integer, allocatable :: outlet_gindices_local(:) ! Local array of global indices of outlets on this task + real(r8), allocatable :: outlet_discharges_local(:) ! Local array of discharges for these outlets + integer :: local_outlet_count + integer, allocatable :: all_outlet_gindices(:) ! Gathered on master + real(r8), allocatable :: all_outlet_discharges(:) ! Gathered on master + integer, allocatable :: recvcounts(:), displs(:) ! For MPI_Gatherv + type(outlet_discharge_info_type), allocatable :: sorted_outlets(:) + integer :: num_all_outlets_global, current_top_n_count, target_gidx + real(r8) :: sum_discharge_top_n + real(r8), allocatable :: qgwl_correction_values(:) ! Values to apply + integer, allocatable :: qgwl_correction_gindices(:) ! Global indices for these corrections real(r8), allocatable :: qgwl_correction_local(:) ! Local portion of correction array real(r8) :: qgwl_to_redistribute ! Amount to redistribute (Scenario A or B) real(r8) :: local_total_outlet_discharge, global_total_outlet_discharge @@ -2248,7 +2264,6 @@ subroutine Rtmrun(rstwr,nlend,rdate) real(r8) :: global_total_liquid_before, global_total_liquid_after real(r8) :: local_correction_sum, global_correction_sum real(r8) :: correction_local(1,1), correction_global(1) - real(r8) :: correction_ratio ! Ratio to apply: qgwl_to_redistribute / global_total_outlet_discharge character(len=*),parameter :: subname = '(Rtmrun) ' !----------------------------------------------------------------------- @@ -2399,21 +2414,21 @@ subroutine Rtmrun(rstwr,nlend,rdate) enddo ! Use reproducible sum for bit-for-bit reproducibility across PE layouts - ! Combine both positive and negative sums in one call for efficiency block - real(r8) :: local_sums(2,1), global_sums(2) - - ! Pack both positive and negative sums - local_sums(1,1) = local_positive_qgwl_sum - local_sums(2,1) = local_negative_qgwl_sum + real(r8) :: pos_local(1,1), pos_global(1) + real(r8) :: neg_local(1,1), neg_global(1) - ! Single reproducible sum call for both fields - call shr_reprosum_calc(local_sums, global_sums, 2, 1, 1, & + ! Reproducible sum for positive qgwl + pos_local(1,1) = local_positive_qgwl_sum + call shr_reprosum_calc(pos_local, pos_global, 1, 1, 1, & commid=mpicom_rof) + global_positive_qgwl_sum = pos_global(1) - ! Unpack results - global_positive_qgwl_sum = global_sums(1) - global_negative_qgwl_sum = global_sums(2) + ! Reproducible sum for negative qgwl + neg_local(1,1) = local_negative_qgwl_sum + call shr_reprosum_calc(neg_local, neg_global, 1, 1, 1, & + commid=mpicom_rof) + global_negative_qgwl_sum = neg_global(1) ! Diagnostic output only when do_budget == 3 if (masterproc .and. do_budget == 3) then @@ -2437,19 +2452,13 @@ subroutine Rtmrun(rstwr,nlend,rdate) ! Proportionally reduce positive qgwl to offset negatives, zero out negative cells ! This ensures NO negative qgwl enters the ocean at any grid cell - ! Master PE calculates scaling factor, then broadcasts for bit-for-bit reproducibility - if (masterproc) then - if (global_positive_qgwl_sum > 0.0_r8) then - ! Calculate scaling factor: (positive - |negative|) / positive = net / positive - scaling_factor = net_global_qgwl / global_positive_qgwl_sum - else - scaling_factor = 1.0_r8 ! No positive qgwl, keep as is - endif + if (global_positive_qgwl_sum > 0.0_r8) then + ! Calculate scaling factor: (positive - |negative|) / positive = net / positive + scaling_factor = net_global_qgwl / global_positive_qgwl_sum + else + scaling_factor = 1.0_r8 ! No positive qgwl, keep as is endif - ! Broadcast scaling factor from master to all PEs - call MPI_Bcast(scaling_factor, 1, MPI_REAL8, 0, mpicom_rof, ier) - do nr = rtmCTL%begr, rtmCTL%endr do nt = 1, nt_rtm if (rtmCTL%qgwl(nr, nt) > 0.0_r8) then @@ -2989,11 +2998,10 @@ subroutine Rtmrun(rstwr,nlend,rdate) if (abs(qgwl_to_redistribute) > 1.0e-15_r8) then ! Only proceed if there's something to redistribute call t_startf('mosartr_qgwl_redir_dist') - ! Calculate total discharge from outlets with discharge > 1.0 m3/s - ! This threshold filters out small outlets and ensures numerical stability + ! Calculate total discharge from all outlets globally using reprosum local_total_outlet_discharge = 0.0_r8 do nr = rtmCTL%begr, rtmCTL%endr - if (rtmCTL%mask(nr) == 3 .and. rtmCTL%runoffocn(nr, nt_nliq) > 1.0_r8) then + if (rtmCTL%mask(nr) == 3) then ! Outlet cell local_total_outlet_discharge = local_total_outlet_discharge + rtmCTL%runoffocn(nr, nt_nliq) endif enddo @@ -3030,28 +3038,19 @@ subroutine Rtmrun(rstwr,nlend,rdate) endif endif - ! Master PE calculates the correction ratio, then broadcast to all PEs - ! This ensures bit-for-bit reproducibility across different PE layouts (PEM test) + ! Each PE calculates corrections for its local outlets proportionally if (abs(global_total_outlet_discharge) > 1.0e-9_r8) then - if (masterproc) then - correction_ratio = qgwl_to_redistribute / global_total_outlet_discharge - endif - - ! Broadcast the correction ratio from master to all PEs - call MPI_Bcast(correction_ratio, 1, MPI_REAL8, 0, mpicom_rof, ier) - - ! All PEs apply the same ratio to their local outlets > 1.0 m3/s do nr = rtmCTL%begr, rtmCTL%endr - if (rtmCTL%mask(nr) == 3 .and. rtmCTL%runoffocn(nr, nt_nliq) > 1.0_r8) then - ! Apply the correction ratio - qgwl_correction_local(nr) = rtmCTL%runoffocn(nr, nt_nliq) * correction_ratio + if (rtmCTL%mask(nr) == 3) then ! Outlet cell + ! Proportional correction based on this outlet's discharge + qgwl_correction_local(nr) = (rtmCTL%runoffocn(nr, nt_nliq) / global_total_outlet_discharge) * qgwl_to_redistribute endif enddo else ! If total discharge is zero, no redistribution (corrections remain 0) if (masterproc) then write(iulog, *) trim(subname), & - 'Warning: Cannot redistribute qgwl - total outlet discharge (>1 m3/s) is zero' + 'Warning: Cannot redistribute qgwl - total outlet discharge is zero' endif endif @@ -5349,4 +5348,23 @@ end subroutine SubTimestep !----------------------------------------------------------------------- +subroutine sort_outlets_by_discharge_desc(outlets_array, count) !TZ negative runoff + implicit none + integer, intent(in) :: count + type(outlet_discharge_info_type), intent(inout) :: outlets_array(:) ! Assumed-shape + integer :: i, j + type(outlet_discharge_info_type) :: temp_outlet_info + ! Bubble sort + if (count < 2) return + do i = 1, count - 1 + do j = 1, count - i + if (outlets_array(j)%discharge < outlets_array(j+1)%discharge) then + temp_outlet_info = outlets_array(j) + outlets_array(j) = outlets_array(j+1) + outlets_array(j+1) = temp_outlet_info + endif + enddo + enddo +end subroutine sort_outlets_by_discharge_desc + end module RtmMod \ No newline at end of file From c2ca073bc0bc004ec606f57fdb0c03889cf33260 Mon Sep 17 00:00:00 2001 From: Tian Zhou Date: Thu, 30 Oct 2025 14:55:10 -0700 Subject: [PATCH 10/17] fix reprosum error --- components/mosart/src/riverroute/RtmMod.F90 | 87 +++++++++++++-------- 1 file changed, 56 insertions(+), 31 deletions(-) diff --git a/components/mosart/src/riverroute/RtmMod.F90 b/components/mosart/src/riverroute/RtmMod.F90 index 0dc179c42c24..9121593e3aec 100644 --- a/components/mosart/src/riverroute/RtmMod.F90 +++ b/components/mosart/src/riverroute/RtmMod.F90 @@ -2240,9 +2240,10 @@ subroutine Rtmrun(rstwr,nlend,rdate) !scs ! Variables for negative runoff redirection - + real(r8) :: local_positive_qgwl_sum, local_negative_qgwl_sum real(r8) :: net_global_qgwl, original_cell_qgwl, reduction, scaling_factor + integer :: local_neg_count, global_neg_count integer, allocatable :: outlet_gindices_local(:) ! Local array of global indices of outlets on this task real(r8), allocatable :: outlet_discharges_local(:) ! Local array of discharges for these outlets @@ -2264,6 +2265,7 @@ subroutine Rtmrun(rstwr,nlend,rdate) real(r8) :: global_total_liquid_before, global_total_liquid_after real(r8) :: local_correction_sum, global_correction_sum real(r8) :: correction_local(1,1), correction_global(1) + real(r8) :: correction_ratio ! Ratio to apply: qgwl_to_redistribute / global_total_outlet_discharge character(len=*),parameter :: subname = '(Rtmrun) ' !----------------------------------------------------------------------- @@ -2404,61 +2406,63 @@ subroutine Rtmrun(rstwr,nlend,rdate) if (redirect_negative_qgwl_flag) then local_negative_qgwl_sum = 0.0_r8 ! This will sum local NEGATIVE qgwl local_positive_qgwl_sum = 0.0_r8 ! This will sum local POSITIVE qgwl + local_neg_count = 0 do nr = rtmCTL%begr, rtmCTL%endr if (rtmCTL%qgwl(nr, nt_nliq) > 0.0_r8) then local_positive_qgwl_sum = local_positive_qgwl_sum + rtmCTL%qgwl(nr, nt_nliq) elseif (rtmCTL%qgwl(nr, nt_nliq) < 0.0_r8) then local_negative_qgwl_sum = local_negative_qgwl_sum + rtmCTL%qgwl(nr, nt_nliq) + local_neg_count = local_neg_count + 1 endif enddo - ! Use reproducible sum for bit-for-bit reproducibility across PE layouts - block - real(r8) :: pos_local(1,1), pos_global(1) - real(r8) :: neg_local(1,1), neg_global(1) + ! Count total negative cells globally + call MPI_Reduce(local_neg_count, global_neg_count, 1, MPI_INTEGER, MPI_SUM, 0, mpicom_rof, ier) - ! Reproducible sum for positive qgwl - pos_local(1,1) = local_positive_qgwl_sum - call shr_reprosum_calc(pos_local, pos_global, 1, 1, 1, & - commid=mpicom_rof) - global_positive_qgwl_sum = pos_global(1) + ! Use combined reproducible sum for bit-for-bit reproducibility across PE layouts + ! This sums both positive and negative qgwl in a single call for efficiency + block + real(r8) :: local_sums(1,2), global_sums(2) - ! Reproducible sum for negative qgwl - neg_local(1,1) = local_negative_qgwl_sum - call shr_reprosum_calc(neg_local, neg_global, 1, 1, 1, & + ! Combined reprosum: arr(dsummands, nflds) where nflds=2 + local_sums(1,1) = local_positive_qgwl_sum + local_sums(1,2) = local_negative_qgwl_sum + call shr_reprosum_calc(local_sums, global_sums, 1, 1, 2, & commid=mpicom_rof) - global_negative_qgwl_sum = neg_global(1) + global_positive_qgwl_sum = global_sums(1) + global_negative_qgwl_sum = global_sums(2) ! Diagnostic output only when do_budget == 3 if (masterproc .and. do_budget == 3) then - write(iulog, '(A,ES24.16)') trim(subname)//' REPROSUM Positive Sum = ', global_positive_qgwl_sum - write(iulog, '(A,ES24.16)') trim(subname)//' REPROSUM Negative Sum = ', global_negative_qgwl_sum - write(iulog, '(A,ES24.16)') trim(subname)//' REPROSUM Net Sum = ', global_positive_qgwl_sum + global_negative_qgwl_sum + write(iulog, '(A,I10)') trim(subname)//' Global count of negative qgwl cells = ', global_neg_count + write(iulog, '(A,ES24.16)') trim(subname)//' Global Positive qgwl Sum = ', global_positive_qgwl_sum + write(iulog, '(A,ES24.16)') trim(subname)//' Global Negative qgwl Sum = ', global_negative_qgwl_sum + write(iulog, '(A,ES24.16)') trim(subname)//' Global Net qgwl Sum = ', global_positive_qgwl_sum + global_negative_qgwl_sum endif end block net_global_qgwl = global_positive_qgwl_sum + global_negative_qgwl_sum - if (masterproc .and. do_budget == 3) then - write(iulog, *) trim(subname), 'Debug QGWL Balancing: Global Positive Sum =', global_positive_qgwl_sum - write(iulog, *) trim(subname), 'Debug QGWL Balancing: Global Negative Sum =', global_negative_qgwl_sum - write(iulog, *) trim(subname), 'Debug QGWL Balancing: Global Net Sum =', net_global_qgwl - endif - ! Decide how to set TRunoff%qgwl for local routing if (net_global_qgwl >= 0.0_r8) then ! --- SCENARIO A: Net global QGWL is non-negative. --- ! Proportionally reduce positive qgwl to offset negatives, zero out negative cells ! This ensures NO negative qgwl enters the ocean at any grid cell - if (global_positive_qgwl_sum > 0.0_r8) then - ! Calculate scaling factor: (positive - |negative|) / positive = net / positive - scaling_factor = net_global_qgwl / global_positive_qgwl_sum - else - scaling_factor = 1.0_r8 ! No positive qgwl, keep as is + ! Master PE calculates scaling factor, then broadcasts for bit-for-bit reproducibility + if (masterproc) then + if (global_positive_qgwl_sum > 0.0_r8) then + ! Calculate scaling factor: (positive - |negative|) / positive = net / positive + scaling_factor = net_global_qgwl / global_positive_qgwl_sum + else + scaling_factor = 1.0_r8 ! No positive qgwl, keep as is + endif endif + ! Broadcast scaling factor from master to all PEs + call MPI_Bcast(scaling_factor, 1, MPI_REAL8, 0, mpicom_rof, ier) + do nr = rtmCTL%begr, rtmCTL%endr do nt = 1, nt_rtm if (rtmCTL%qgwl(nr, nt) > 0.0_r8) then @@ -3038,12 +3042,33 @@ subroutine Rtmrun(rstwr,nlend,rdate) endif endif - ! Each PE calculates corrections for its local outlets proportionally + ! Master PE calculates the correction ratio, then broadcast to all PEs + ! This ensures bit-for-bit reproducibility across different PE layouts (PEM test) if (abs(global_total_outlet_discharge) > 1.0e-9_r8) then + if (masterproc) then + correction_ratio = qgwl_to_redistribute / global_total_outlet_discharge + + ! Check if correction ratio magnitude exceeds 100% + if (correction_ratio < -1.0_r8) then + call shr_sys_flush(iulog) + write(iulog, *) trim(subname), & + 'WARNING: Correction ratio < -100% (', correction_ratio, ').' + write(iulog, *) trim(subname), & + 'Negative runoff to ocean is unavoidable as negative qgwl magnitude (', & + abs(qgwl_to_redistribute), ' m3/s) exceeds total outlet discharge (', & + global_total_outlet_discharge, ' m3/s).' + call shr_sys_flush(iulog) + endif + endif + + ! Broadcast the correction ratio from master to all PEs + call MPI_Bcast(correction_ratio, 1, MPI_REAL8, 0, mpicom_rof, ier) + + ! All PEs apply the same ratio to their local outlets do nr = rtmCTL%begr, rtmCTL%endr if (rtmCTL%mask(nr) == 3) then ! Outlet cell - ! Proportional correction based on this outlet's discharge - qgwl_correction_local(nr) = (rtmCTL%runoffocn(nr, nt_nliq) / global_total_outlet_discharge) * qgwl_to_redistribute + ! Apply the correction ratio + qgwl_correction_local(nr) = rtmCTL%runoffocn(nr, nt_nliq) * correction_ratio endif enddo else From 1bf9a82089a5056732011afc9b2d83290e553687 Mon Sep 17 00:00:00 2001 From: Tian Zhou Date: Thu, 30 Oct 2025 15:37:15 -0700 Subject: [PATCH 11/17] remove some diag calculations --- components/mosart/src/riverroute/RtmMod.F90 | 7 ------- 1 file changed, 7 deletions(-) diff --git a/components/mosart/src/riverroute/RtmMod.F90 b/components/mosart/src/riverroute/RtmMod.F90 index 9121593e3aec..739ffac84019 100644 --- a/components/mosart/src/riverroute/RtmMod.F90 +++ b/components/mosart/src/riverroute/RtmMod.F90 @@ -2243,7 +2243,6 @@ subroutine Rtmrun(rstwr,nlend,rdate) real(r8) :: local_positive_qgwl_sum, local_negative_qgwl_sum real(r8) :: net_global_qgwl, original_cell_qgwl, reduction, scaling_factor - integer :: local_neg_count, global_neg_count integer, allocatable :: outlet_gindices_local(:) ! Local array of global indices of outlets on this task real(r8), allocatable :: outlet_discharges_local(:) ! Local array of discharges for these outlets @@ -2406,20 +2405,15 @@ subroutine Rtmrun(rstwr,nlend,rdate) if (redirect_negative_qgwl_flag) then local_negative_qgwl_sum = 0.0_r8 ! This will sum local NEGATIVE qgwl local_positive_qgwl_sum = 0.0_r8 ! This will sum local POSITIVE qgwl - local_neg_count = 0 do nr = rtmCTL%begr, rtmCTL%endr if (rtmCTL%qgwl(nr, nt_nliq) > 0.0_r8) then local_positive_qgwl_sum = local_positive_qgwl_sum + rtmCTL%qgwl(nr, nt_nliq) elseif (rtmCTL%qgwl(nr, nt_nliq) < 0.0_r8) then local_negative_qgwl_sum = local_negative_qgwl_sum + rtmCTL%qgwl(nr, nt_nliq) - local_neg_count = local_neg_count + 1 endif enddo - ! Count total negative cells globally - call MPI_Reduce(local_neg_count, global_neg_count, 1, MPI_INTEGER, MPI_SUM, 0, mpicom_rof, ier) - ! Use combined reproducible sum for bit-for-bit reproducibility across PE layouts ! This sums both positive and negative qgwl in a single call for efficiency block @@ -2435,7 +2429,6 @@ subroutine Rtmrun(rstwr,nlend,rdate) ! Diagnostic output only when do_budget == 3 if (masterproc .and. do_budget == 3) then - write(iulog, '(A,I10)') trim(subname)//' Global count of negative qgwl cells = ', global_neg_count write(iulog, '(A,ES24.16)') trim(subname)//' Global Positive qgwl Sum = ', global_positive_qgwl_sum write(iulog, '(A,ES24.16)') trim(subname)//' Global Negative qgwl Sum = ', global_negative_qgwl_sum write(iulog, '(A,ES24.16)') trim(subname)//' Global Net qgwl Sum = ', global_positive_qgwl_sum + global_negative_qgwl_sum From 39fa31cade4c3ec5f9c8c2d51edb3b48e2c9e94f Mon Sep 17 00:00:00 2001 From: Tian Zhou Date: Thu, 30 Oct 2025 23:26:40 -0700 Subject: [PATCH 12/17] detect reprosum problem --- components/mosart/src/riverroute/RtmMod.F90 | 43 ++++++++++++++++++++- 1 file changed, 41 insertions(+), 2 deletions(-) diff --git a/components/mosart/src/riverroute/RtmMod.F90 b/components/mosart/src/riverroute/RtmMod.F90 index 739ffac84019..693aef3b325f 100644 --- a/components/mosart/src/riverroute/RtmMod.F90 +++ b/components/mosart/src/riverroute/RtmMod.F90 @@ -2418,19 +2418,35 @@ subroutine Rtmrun(rstwr,nlend,rdate) ! This sums both positive and negative qgwl in a single call for efficiency block real(r8) :: local_sums(1,2), global_sums(2) + real(r8) :: check_pos_mpi, check_neg_mpi ! For debugging PEM ! Combined reprosum: arr(dsummands, nflds) where nflds=2 local_sums(1,1) = local_positive_qgwl_sum local_sums(1,2) = local_negative_qgwl_sum + + ! DEBUG: Check what reprosum receives on master PE + if (masterproc .and. do_budget == 3) then + write(iulog, '(A,ES24.16)') trim(subname)//' [PE0] LOCAL Positive input = ', local_sums(1,1) + write(iulog, '(A,ES24.16)') trim(subname)//' [PE0] LOCAL Negative input = ', local_sums(1,2) + endif + call shr_reprosum_calc(local_sums, global_sums, 1, 1, 2, & commid=mpicom_rof) global_positive_qgwl_sum = global_sums(1) global_negative_qgwl_sum = global_sums(2) + ! DEBUG: Compare reprosum with MPI_Allreduce for PEM debugging + call MPI_Allreduce(local_positive_qgwl_sum, check_pos_mpi, 1, MPI_REAL8, MPI_SUM, mpicom_rof, ier) + call MPI_Allreduce(local_negative_qgwl_sum, check_neg_mpi, 1, MPI_REAL8, MPI_SUM, mpicom_rof, ier) + ! Diagnostic output only when do_budget == 3 if (masterproc .and. do_budget == 3) then - write(iulog, '(A,ES24.16)') trim(subname)//' Global Positive qgwl Sum = ', global_positive_qgwl_sum - write(iulog, '(A,ES24.16)') trim(subname)//' Global Negative qgwl Sum = ', global_negative_qgwl_sum + write(iulog, '(A,ES24.16)') trim(subname)//' REPROSUM Positive Sum = ', global_positive_qgwl_sum + write(iulog, '(A,ES24.16)') trim(subname)//' REPROSUM Negative Sum = ', global_negative_qgwl_sum + write(iulog, '(A,ES24.16)') trim(subname)//' MPI_ALLREDUCE Positive Sum = ', check_pos_mpi + write(iulog, '(A,ES24.16)') trim(subname)//' MPI_ALLREDUCE Negative Sum = ', check_neg_mpi + write(iulog, '(A,ES24.16)') trim(subname)//' Diff (REPROSUM - MPI) Positive = ', global_positive_qgwl_sum - check_pos_mpi + write(iulog, '(A,ES24.16)') trim(subname)//' Diff (REPROSUM - MPI) Negative = ', global_negative_qgwl_sum - check_neg_mpi write(iulog, '(A,ES24.16)') trim(subname)//' Global Net qgwl Sum = ', global_positive_qgwl_sum + global_negative_qgwl_sum endif end block @@ -2451,11 +2467,22 @@ subroutine Rtmrun(rstwr,nlend,rdate) else scaling_factor = 1.0_r8 ! No positive qgwl, keep as is endif + ! DEBUG: Show scaling factor calculation + if (do_budget == 3) then + write(iulog, '(A,ES24.16)') trim(subname)//' [Scenario A] net_global_qgwl = ', net_global_qgwl + write(iulog, '(A,ES24.16)') trim(subname)//' [Scenario A] global_positive_qgwl_sum = ', global_positive_qgwl_sum + write(iulog, '(A,ES24.16)') trim(subname)//' [Scenario A] scaling_factor (before bcast) = ', scaling_factor + endif endif ! Broadcast scaling factor from master to all PEs call MPI_Bcast(scaling_factor, 1, MPI_REAL8, 0, mpicom_rof, ier) + ! DEBUG: Verify broadcast on all PEs + if (do_budget == 3) then + write(iulog, '(A,I5,A,ES24.16)') trim(subname)//' [PE ', iam, '] scaling_factor (after bcast) = ', scaling_factor + endif + do nr = rtmCTL%begr, rtmCTL%endr do nt = 1, nt_rtm if (rtmCTL%qgwl(nr, nt) > 0.0_r8) then @@ -3041,6 +3068,13 @@ subroutine Rtmrun(rstwr,nlend,rdate) if (masterproc) then correction_ratio = qgwl_to_redistribute / global_total_outlet_discharge + ! DEBUG: Show correction ratio calculation + if (do_budget == 3) then + write(iulog, '(A,ES24.16)') trim(subname)//' [Scenario B] qgwl_to_redistribute = ', qgwl_to_redistribute + write(iulog, '(A,ES24.16)') trim(subname)//' [Scenario B] global_total_outlet_discharge = ', global_total_outlet_discharge + write(iulog, '(A,ES24.16)') trim(subname)//' [Scenario B] correction_ratio (before bcast) = ', correction_ratio + endif + ! Check if correction ratio magnitude exceeds 100% if (correction_ratio < -1.0_r8) then call shr_sys_flush(iulog) @@ -3057,6 +3091,11 @@ subroutine Rtmrun(rstwr,nlend,rdate) ! Broadcast the correction ratio from master to all PEs call MPI_Bcast(correction_ratio, 1, MPI_REAL8, 0, mpicom_rof, ier) + ! DEBUG: Verify broadcast on all PEs + if (do_budget == 3) then + write(iulog, '(A,I5,A,ES24.16)') trim(subname)//' [PE ', iam, '] correction_ratio (after bcast) = ', correction_ratio + endif + ! All PEs apply the same ratio to their local outlets do nr = rtmCTL%begr, rtmCTL%endr if (rtmCTL%mask(nr) == 3) then ! Outlet cell From f50555635123a7449bff6e53f64016c78bf4180f Mon Sep 17 00:00:00 2001 From: Tian Zhou Date: Fri, 31 Oct 2025 00:39:32 -0700 Subject: [PATCH 13/17] revert to seperate reprosum calls --- components/mosart/src/riverroute/RtmMod.F90 | 38 ++++++++------------- 1 file changed, 14 insertions(+), 24 deletions(-) diff --git a/components/mosart/src/riverroute/RtmMod.F90 b/components/mosart/src/riverroute/RtmMod.F90 index 693aef3b325f..cc344518ca42 100644 --- a/components/mosart/src/riverroute/RtmMod.F90 +++ b/components/mosart/src/riverroute/RtmMod.F90 @@ -2414,39 +2414,29 @@ subroutine Rtmrun(rstwr,nlend,rdate) endif enddo - ! Use combined reproducible sum for bit-for-bit reproducibility across PE layouts - ! This sums both positive and negative qgwl in a single call for efficiency + ! Use separate reproducible sum calls for bit-for-bit reproducibility across PE layouts + ! NOTE: Combined reprosum (nflds=2) was found to be non-reproducible for positive field + ! while negative field was reproducible, suggesting a bug in multi-field reprosum block - real(r8) :: local_sums(1,2), global_sums(2) - real(r8) :: check_pos_mpi, check_neg_mpi ! For debugging PEM + real(r8) :: pos_local(1,1), pos_global(1) + real(r8) :: neg_local(1,1), neg_global(1) - ! Combined reprosum: arr(dsummands, nflds) where nflds=2 - local_sums(1,1) = local_positive_qgwl_sum - local_sums(1,2) = local_negative_qgwl_sum - - ! DEBUG: Check what reprosum receives on master PE - if (masterproc .and. do_budget == 3) then - write(iulog, '(A,ES24.16)') trim(subname)//' [PE0] LOCAL Positive input = ', local_sums(1,1) - write(iulog, '(A,ES24.16)') trim(subname)//' [PE0] LOCAL Negative input = ', local_sums(1,2) - endif - - call shr_reprosum_calc(local_sums, global_sums, 1, 1, 2, & + ! Separate reprosum for positive qgwl + pos_local(1,1) = local_positive_qgwl_sum + call shr_reprosum_calc(pos_local, pos_global, 1, 1, 1, & commid=mpicom_rof) - global_positive_qgwl_sum = global_sums(1) - global_negative_qgwl_sum = global_sums(2) + global_positive_qgwl_sum = pos_global(1) - ! DEBUG: Compare reprosum with MPI_Allreduce for PEM debugging - call MPI_Allreduce(local_positive_qgwl_sum, check_pos_mpi, 1, MPI_REAL8, MPI_SUM, mpicom_rof, ier) - call MPI_Allreduce(local_negative_qgwl_sum, check_neg_mpi, 1, MPI_REAL8, MPI_SUM, mpicom_rof, ier) + ! Separate reprosum for negative qgwl + neg_local(1,1) = local_negative_qgwl_sum + call shr_reprosum_calc(neg_local, neg_global, 1, 1, 1, & + commid=mpicom_rof) + global_negative_qgwl_sum = neg_global(1) ! Diagnostic output only when do_budget == 3 if (masterproc .and. do_budget == 3) then write(iulog, '(A,ES24.16)') trim(subname)//' REPROSUM Positive Sum = ', global_positive_qgwl_sum write(iulog, '(A,ES24.16)') trim(subname)//' REPROSUM Negative Sum = ', global_negative_qgwl_sum - write(iulog, '(A,ES24.16)') trim(subname)//' MPI_ALLREDUCE Positive Sum = ', check_pos_mpi - write(iulog, '(A,ES24.16)') trim(subname)//' MPI_ALLREDUCE Negative Sum = ', check_neg_mpi - write(iulog, '(A,ES24.16)') trim(subname)//' Diff (REPROSUM - MPI) Positive = ', global_positive_qgwl_sum - check_pos_mpi - write(iulog, '(A,ES24.16)') trim(subname)//' Diff (REPROSUM - MPI) Negative = ', global_negative_qgwl_sum - check_neg_mpi write(iulog, '(A,ES24.16)') trim(subname)//' Global Net qgwl Sum = ', global_positive_qgwl_sum + global_negative_qgwl_sum endif end block From a91ff9e1bb871fda05c000a942c2748cb3dc182f Mon Sep 17 00:00:00 2001 From: Tian Zhou Date: Fri, 31 Oct 2025 02:01:20 -0700 Subject: [PATCH 14/17] another attempt --- components/mosart/src/riverroute/RtmMod.F90 | 47 ++++++++++++--------- 1 file changed, 26 insertions(+), 21 deletions(-) diff --git a/components/mosart/src/riverroute/RtmMod.F90 b/components/mosart/src/riverroute/RtmMod.F90 index cc344518ca42..9cc08f3c209a 100644 --- a/components/mosart/src/riverroute/RtmMod.F90 +++ b/components/mosart/src/riverroute/RtmMod.F90 @@ -44,7 +44,7 @@ module RtmMod use MOSART_physics_mod, only : Euler use MOSART_physics_mod, only : updatestate_hillslope, updatestate_subnetwork, & updatestate_mainchannel - use MOSART_BGC_type, only : TSedi, TSedi_para, MOSART_sediment_init + use MOSART_BGC_type, only : TSedi, TSedi_para, MOSART_sediment_init, TINYVALUE_s use MOSART_RES_type, only : Tres, MOSART_reservoir_sed_init, Tres_para !#ifdef INCLUDE_WRM use WRM_type_mod , only : ctlSubwWRM, WRMUnit, StorWater @@ -2241,7 +2241,7 @@ subroutine Rtmrun(rstwr,nlend,rdate) ! Variables for negative runoff redirection - real(r8) :: local_positive_qgwl_sum, local_negative_qgwl_sum + real(r8) :: local_positive_qgwl_sum, local_negative_qgwl_sum, local_total_qgwl_sum real(r8) :: net_global_qgwl, original_cell_qgwl, reduction, scaling_factor integer, allocatable :: outlet_gindices_local(:) ! Local array of global indices of outlets on this task @@ -2405,39 +2405,44 @@ subroutine Rtmrun(rstwr,nlend,rdate) if (redirect_negative_qgwl_flag) then local_negative_qgwl_sum = 0.0_r8 ! This will sum local NEGATIVE qgwl local_positive_qgwl_sum = 0.0_r8 ! This will sum local POSITIVE qgwl + local_total_qgwl_sum = 0.0_r8 ! This will sum ALL qgwl (for verification) + ! Use TINYVALUE_s threshold to avoid non-reproducibility from near-zero values + ! that can flip sign across different PE layouts due to floating-point noise do nr = rtmCTL%begr, rtmCTL%endr - if (rtmCTL%qgwl(nr, nt_nliq) > 0.0_r8) then + if (rtmCTL%qgwl(nr, nt_nliq) > TINYVALUE_s) then local_positive_qgwl_sum = local_positive_qgwl_sum + rtmCTL%qgwl(nr, nt_nliq) - elseif (rtmCTL%qgwl(nr, nt_nliq) < 0.0_r8) then + elseif (rtmCTL%qgwl(nr, nt_nliq) < -TINYVALUE_s) then local_negative_qgwl_sum = local_negative_qgwl_sum + rtmCTL%qgwl(nr, nt_nliq) endif + ! Sum ALL qgwl values without threshold for verification + local_total_qgwl_sum = local_total_qgwl_sum + rtmCTL%qgwl(nr, nt_nliq) enddo - ! Use separate reproducible sum calls for bit-for-bit reproducibility across PE layouts - ! NOTE: Combined reprosum (nflds=2) was found to be non-reproducible for positive field - ! while negative field was reproducible, suggesting a bug in multi-field reprosum + ! Use combined reproducible sum for efficiency (3 fields: positive, negative, total) block - real(r8) :: pos_local(1,1), pos_global(1) - real(r8) :: neg_local(1,1), neg_global(1) + real(r8) :: local_sums(1,3), global_sums(3) + real(r8) :: global_total_qgwl_sum - ! Separate reprosum for positive qgwl - pos_local(1,1) = local_positive_qgwl_sum - call shr_reprosum_calc(pos_local, pos_global, 1, 1, 1, & - commid=mpicom_rof) - global_positive_qgwl_sum = pos_global(1) + ! Pack all three sums into combined array + local_sums(1,1) = local_positive_qgwl_sum + local_sums(1,2) = local_negative_qgwl_sum + local_sums(1,3) = local_total_qgwl_sum - ! Separate reprosum for negative qgwl - neg_local(1,1) = local_negative_qgwl_sum - call shr_reprosum_calc(neg_local, neg_global, 1, 1, 1, & + ! Single combined reprosum call for all three fields + call shr_reprosum_calc(local_sums, global_sums, 1, 1, 3, & commid=mpicom_rof) - global_negative_qgwl_sum = neg_global(1) + global_positive_qgwl_sum = global_sums(1) + global_negative_qgwl_sum = global_sums(2) + global_total_qgwl_sum = global_sums(3) ! Diagnostic output only when do_budget == 3 if (masterproc .and. do_budget == 3) then - write(iulog, '(A,ES24.16)') trim(subname)//' REPROSUM Positive Sum = ', global_positive_qgwl_sum - write(iulog, '(A,ES24.16)') trim(subname)//' REPROSUM Negative Sum = ', global_negative_qgwl_sum - write(iulog, '(A,ES24.16)') trim(subname)//' Global Net qgwl Sum = ', global_positive_qgwl_sum + global_negative_qgwl_sum + write(iulog, '(A,ES24.16)') trim(subname)//' REPROSUM Positive Sum (with threshold) = ', global_positive_qgwl_sum + write(iulog, '(A,ES24.16)') trim(subname)//' REPROSUM Negative Sum (with threshold) = ', global_negative_qgwl_sum + write(iulog, '(A,ES24.16)') trim(subname)//' REPROSUM Total Sum (all values) = ', global_total_qgwl_sum + write(iulog, '(A,ES24.16)') trim(subname)//' Net from pos+neg = ', global_positive_qgwl_sum + global_negative_qgwl_sum + write(iulog, '(A,ES24.16)') trim(subname)//' Difference (Total - Net) = ', global_total_qgwl_sum - (global_positive_qgwl_sum + global_negative_qgwl_sum) endif end block From a28ca9cf9e29912899c6fdfeb774479ae6f5a806 Mon Sep 17 00:00:00 2001 From: Tian Zhou Date: Fri, 31 Oct 2025 23:54:46 -0700 Subject: [PATCH 15/17] PEM test passed --- components/mosart/src/riverroute/RtmMod.F90 | 97 +++++++++++---------- 1 file changed, 52 insertions(+), 45 deletions(-) diff --git a/components/mosart/src/riverroute/RtmMod.F90 b/components/mosart/src/riverroute/RtmMod.F90 index 9cc08f3c209a..664492b5d64e 100644 --- a/components/mosart/src/riverroute/RtmMod.F90 +++ b/components/mosart/src/riverroute/RtmMod.F90 @@ -44,7 +44,7 @@ module RtmMod use MOSART_physics_mod, only : Euler use MOSART_physics_mod, only : updatestate_hillslope, updatestate_subnetwork, & updatestate_mainchannel - use MOSART_BGC_type, only : TSedi, TSedi_para, MOSART_sediment_init, TINYVALUE_s + use MOSART_BGC_type, only : TSedi, TSedi_para, MOSART_sediment_init use MOSART_RES_type, only : Tres, MOSART_reservoir_sed_init, Tres_para !#ifdef INCLUDE_WRM use WRM_type_mod , only : ctlSubwWRM, WRMUnit, StorWater @@ -143,7 +143,8 @@ module RtmMod ! Variables for TNR redirection integer, parameter :: num_top_outlets_for_qgwl = 500 ! Number of top outlets to use - real(r8), save :: global_positive_qgwl_sum = 0.0_r8 ! Sum of all positive qgwl + real(r8), parameter :: TINYVALUE_s = 1.0e-14_r8 ! Threshold for near-zero qgwl values + real(r8), save :: global_positive_qgwl_sum = 0.0_r8 ! Sum of all positive qgwl real(r8), save :: global_negative_qgwl_sum = 0.0_r8 ! Sum of all negative qgwl real(r8), save :: delt_save ! previous delt @@ -2403,47 +2404,62 @@ subroutine Rtmrun(rstwr,nlend,rdate) global_negative_qgwl_sum = 0.0_r8 ! Reset global sum each step if (redirect_negative_qgwl_flag) then - local_negative_qgwl_sum = 0.0_r8 ! This will sum local NEGATIVE qgwl - local_positive_qgwl_sum = 0.0_r8 ! This will sum local POSITIVE qgwl - local_total_qgwl_sum = 0.0_r8 ! This will sum ALL qgwl (for verification) - - ! Use TINYVALUE_s threshold to avoid non-reproducibility from near-zero values - ! that can flip sign across different PE layouts due to floating-point noise - do nr = rtmCTL%begr, rtmCTL%endr - if (rtmCTL%qgwl(nr, nt_nliq) > TINYVALUE_s) then - local_positive_qgwl_sum = local_positive_qgwl_sum + rtmCTL%qgwl(nr, nt_nliq) - elseif (rtmCTL%qgwl(nr, nt_nliq) < -TINYVALUE_s) then - local_negative_qgwl_sum = local_negative_qgwl_sum + rtmCTL%qgwl(nr, nt_nliq) - endif - ! Sum ALL qgwl values without threshold for verification - local_total_qgwl_sum = local_total_qgwl_sum + rtmCTL%qgwl(nr, nt_nliq) - enddo - - ! Use combined reproducible sum for efficiency (3 fields: positive, negative, total) + ! Use sparse packing approach: only include values exceeding threshold in reprosum + ! This avoids near-zero values that can flip sign across PE layouts block - real(r8) :: local_sums(1,3), global_sums(3) + real(r8), allocatable :: local_qgwl_array(:,:) + real(r8) :: global_sums(2) real(r8) :: global_total_qgwl_sum + integer :: num_positive, num_negative, num_total + integer :: idx_pos, idx_neg, idx_tot + integer :: local_count + integer :: nstep, yr, mon, day, tod + + ! First pass: count how many values exceed threshold + num_positive = 0 + num_negative = 0 + do nr = rtmCTL%begr, rtmCTL%endr + if (rtmCTL%qgwl(nr, nt_nliq) > TINYVALUE_s) then + num_positive = num_positive + 1 + elseif (rtmCTL%qgwl(nr, nt_nliq) < -TINYVALUE_s) then + num_negative = num_negative + 1 + endif + enddo + + ! Total entries to pack + local_count = num_positive + num_negative - ! Pack all three sums into combined array - local_sums(1,1) = local_positive_qgwl_sum - local_sums(1,2) = local_negative_qgwl_sum - local_sums(1,3) = local_total_qgwl_sum + ! Allocate sparse array for reprosum (2 fields: positive and negative) + allocate(local_qgwl_array(local_count, 2)) + local_qgwl_array(:,:) = 0.0_r8 - ! Single combined reprosum call for all three fields - call shr_reprosum_calc(local_sums, global_sums, 1, 1, 3, & + ! Second pass: pack only values exceeding threshold + idx_pos = 0 + idx_neg = 0 + do nr = rtmCTL%begr, rtmCTL%endr + if (rtmCTL%qgwl(nr, nt_nliq) > TINYVALUE_s) then + idx_pos = idx_pos + 1 + local_qgwl_array(idx_pos, 1) = rtmCTL%qgwl(nr, nt_nliq) + elseif (rtmCTL%qgwl(nr, nt_nliq) < -TINYVALUE_s) then + idx_neg = idx_neg + 1 + local_qgwl_array(num_positive + idx_neg, 2) = rtmCTL%qgwl(nr, nt_nliq) + endif + enddo + + ! Reproducible sum with sparse packing + call shr_reprosum_calc(local_qgwl_array, global_sums, local_count, local_count, 2, & commid=mpicom_rof) global_positive_qgwl_sum = global_sums(1) global_negative_qgwl_sum = global_sums(2) - global_total_qgwl_sum = global_sums(3) - ! Diagnostic output only when do_budget == 3 + ! Diagnostic output - total positive, negative, and net sums if (masterproc .and. do_budget == 3) then - write(iulog, '(A,ES24.16)') trim(subname)//' REPROSUM Positive Sum (with threshold) = ', global_positive_qgwl_sum - write(iulog, '(A,ES24.16)') trim(subname)//' REPROSUM Negative Sum (with threshold) = ', global_negative_qgwl_sum - write(iulog, '(A,ES24.16)') trim(subname)//' REPROSUM Total Sum (all values) = ', global_total_qgwl_sum - write(iulog, '(A,ES24.16)') trim(subname)//' Net from pos+neg = ', global_positive_qgwl_sum + global_negative_qgwl_sum - write(iulog, '(A,ES24.16)') trim(subname)//' Difference (Total - Net) = ', global_total_qgwl_sum - (global_positive_qgwl_sum + global_negative_qgwl_sum) + write(iulog, '(A,ES24.16)') trim(subname)//' Global positive qgwl (threshold)=', global_positive_qgwl_sum + write(iulog, '(A,ES24.16)') trim(subname)//' Global negative qgwl (threshold)=', global_negative_qgwl_sum + write(iulog, '(A,ES24.16)') trim(subname)//' Net global qgwl (pos+neg) =', global_positive_qgwl_sum + global_negative_qgwl_sum endif + + deallocate(local_qgwl_array) end block net_global_qgwl = global_positive_qgwl_sum + global_negative_qgwl_sum @@ -2473,18 +2489,14 @@ subroutine Rtmrun(rstwr,nlend,rdate) ! Broadcast scaling factor from master to all PEs call MPI_Bcast(scaling_factor, 1, MPI_REAL8, 0, mpicom_rof, ier) - ! DEBUG: Verify broadcast on all PEs - if (do_budget == 3) then - write(iulog, '(A,I5,A,ES24.16)') trim(subname)//' [PE ', iam, '] scaling_factor (after bcast) = ', scaling_factor - endif - + ! Apply scaling with TINYVALUE_s threshold to ensure reproducibility do nr = rtmCTL%begr, rtmCTL%endr do nt = 1, nt_rtm - if (rtmCTL%qgwl(nr, nt) > 0.0_r8) then + if (rtmCTL%qgwl(nr, nt) > TINYVALUE_s) then ! Positive cell: scale down proportionally TRunoff%qgwl(nr, nt) = rtmCTL%qgwl(nr, nt) * scaling_factor else - ! Negative or zero cell: set to zero + ! Negative, zero, or near-zero cell: set to zero TRunoff%qgwl(nr, nt) = 0.0_r8 endif enddo @@ -3086,11 +3098,6 @@ subroutine Rtmrun(rstwr,nlend,rdate) ! Broadcast the correction ratio from master to all PEs call MPI_Bcast(correction_ratio, 1, MPI_REAL8, 0, mpicom_rof, ier) - ! DEBUG: Verify broadcast on all PEs - if (do_budget == 3) then - write(iulog, '(A,I5,A,ES24.16)') trim(subname)//' [PE ', iam, '] correction_ratio (after bcast) = ', correction_ratio - endif - ! All PEs apply the same ratio to their local outlets do nr = rtmCTL%begr, rtmCTL%endr if (rtmCTL%mask(nr) == 3) then ! Outlet cell From ec68102cf4ab4144e2a0aaed6d2168b7a0cb3bb1 Mon Sep 17 00:00:00 2001 From: Tian Zhou Date: Sat, 1 Nov 2025 01:07:14 -0700 Subject: [PATCH 16/17] remove unused stuff --- components/mosart/src/riverroute/RtmMod.F90 | 55 ++------------------- 1 file changed, 5 insertions(+), 50 deletions(-) diff --git a/components/mosart/src/riverroute/RtmMod.F90 b/components/mosart/src/riverroute/RtmMod.F90 index 664492b5d64e..047b4c5e2a3a 100644 --- a/components/mosart/src/riverroute/RtmMod.F90 +++ b/components/mosart/src/riverroute/RtmMod.F90 @@ -132,17 +132,9 @@ module RtmMod logical :: do_rtmflood logical :: do_rtm - type :: outlet_discharge_info_type - integer :: gidx - real(r8) :: discharge - end type outlet_discharge_info_type - - ! Namelist variable and flag for TNR redirection + ! Namelist variable and flag for negative qgwl redirection logical, public :: redirect_negative_qgwl = .false. ! Namelist control logical, save :: redirect_negative_qgwl_flag = .false. ! Module flag - - ! Variables for TNR redirection - integer, parameter :: num_top_outlets_for_qgwl = 500 ! Number of top outlets to use real(r8), parameter :: TINYVALUE_s = 1.0e-14_r8 ! Threshold for near-zero qgwl values real(r8), save :: global_positive_qgwl_sum = 0.0_r8 ! Sum of all positive qgwl real(r8), save :: global_negative_qgwl_sum = 0.0_r8 ! Sum of all negative qgwl @@ -2241,31 +2233,16 @@ subroutine Rtmrun(rstwr,nlend,rdate) !scs ! Variables for negative runoff redirection - - real(r8) :: local_positive_qgwl_sum, local_negative_qgwl_sum, local_total_qgwl_sum - real(r8) :: net_global_qgwl, original_cell_qgwl, reduction, scaling_factor - - integer, allocatable :: outlet_gindices_local(:) ! Local array of global indices of outlets on this task - real(r8), allocatable :: outlet_discharges_local(:) ! Local array of discharges for these outlets - integer :: local_outlet_count - integer, allocatable :: all_outlet_gindices(:) ! Gathered on master - real(r8), allocatable :: all_outlet_discharges(:) ! Gathered on master - integer, allocatable :: recvcounts(:), displs(:) ! For MPI_Gatherv - type(outlet_discharge_info_type), allocatable :: sorted_outlets(:) - integer :: num_all_outlets_global, current_top_n_count, target_gidx - real(r8) :: sum_discharge_top_n - real(r8), allocatable :: qgwl_correction_values(:) ! Values to apply - integer, allocatable :: qgwl_correction_gindices(:) ! Global indices for these corrections - real(r8), allocatable :: qgwl_correction_local(:) ! Local portion of correction array - real(r8) :: qgwl_to_redistribute ! Amount to redistribute (Scenario A or B) + real(r8) :: net_global_qgwl, scaling_factor + real(r8), allocatable :: qgwl_correction_local(:) + real(r8) :: qgwl_to_redistribute real(r8) :: local_total_outlet_discharge, global_total_outlet_discharge real(r8) :: outlet_discharge_local(1,1), outlet_discharge_global(1) real(r8) :: qgwl_to_discharge_ratio_percent real(r8) :: conservation_error - real(r8) :: global_total_liquid_before, global_total_liquid_after real(r8) :: local_correction_sum, global_correction_sum real(r8) :: correction_local(1,1), correction_global(1) - real(r8) :: correction_ratio ! Ratio to apply: qgwl_to_redistribute / global_total_outlet_discharge + real(r8) :: correction_ratio character(len=*),parameter :: subname = '(Rtmrun) ' !----------------------------------------------------------------------- @@ -2405,7 +2382,6 @@ subroutine Rtmrun(rstwr,nlend,rdate) if (redirect_negative_qgwl_flag) then ! Use sparse packing approach: only include values exceeding threshold in reprosum - ! This avoids near-zero values that can flip sign across PE layouts block real(r8), allocatable :: local_qgwl_array(:,:) real(r8) :: global_sums(2) @@ -5405,25 +5381,4 @@ subroutine SubTimestep end do end subroutine SubTimestep -!----------------------------------------------------------------------- - -subroutine sort_outlets_by_discharge_desc(outlets_array, count) !TZ negative runoff - implicit none - integer, intent(in) :: count - type(outlet_discharge_info_type), intent(inout) :: outlets_array(:) ! Assumed-shape - integer :: i, j - type(outlet_discharge_info_type) :: temp_outlet_info - ! Bubble sort - if (count < 2) return - do i = 1, count - 1 - do j = 1, count - i - if (outlets_array(j)%discharge < outlets_array(j+1)%discharge) then - temp_outlet_info = outlets_array(j) - outlets_array(j) = outlets_array(j+1) - outlets_array(j+1) = temp_outlet_info - endif - enddo - enddo -end subroutine sort_outlets_by_discharge_desc - end module RtmMod \ No newline at end of file From 89d781c9a37e8fde0ca4d792e8cf4d6c9636fc36 Mon Sep 17 00:00:00 2001 From: Tian Zhou Date: Wed, 5 Nov 2025 11:07:32 -0800 Subject: [PATCH 17/17] consolidate diagnostic outputs --- components/mosart/src/riverroute/RtmMod.F90 | 28 ++++++++++----------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/components/mosart/src/riverroute/RtmMod.F90 b/components/mosart/src/riverroute/RtmMod.F90 index 047b4c5e2a3a..20adb886705c 100644 --- a/components/mosart/src/riverroute/RtmMod.F90 +++ b/components/mosart/src/riverroute/RtmMod.F90 @@ -2458,7 +2458,7 @@ subroutine Rtmrun(rstwr,nlend,rdate) if (do_budget == 3) then write(iulog, '(A,ES24.16)') trim(subname)//' [Scenario A] net_global_qgwl = ', net_global_qgwl write(iulog, '(A,ES24.16)') trim(subname)//' [Scenario A] global_positive_qgwl_sum = ', global_positive_qgwl_sum - write(iulog, '(A,ES24.16)') trim(subname)//' [Scenario A] scaling_factor (before bcast) = ', scaling_factor + write(iulog, '(A,ES24.16)') trim(subname)//' [Scenario A] scaling_factor = ', scaling_factor endif endif @@ -3029,7 +3029,19 @@ subroutine Rtmrun(rstwr,nlend,rdate) qgwl_to_discharge_ratio_percent, '% (', qgwl_to_redistribute, ' m3/s)' endif - ! Always check and warn if magnitude is > 5% (regardless of do_budget) + ! Check for concerning redistribution ratios + if (abs(qgwl_to_discharge_ratio_percent) > 100.0_r8) then + call shr_sys_flush(iulog) + write(iulog, *) trim(subname), & + 'WARNING: QGWL_Redist_Ratio magnitude > 100% (', & + qgwl_to_discharge_ratio_percent, '%).' + write(iulog, *) trim(subname), & + 'Negative runoff to ocean is unavoidable as negative qgwl magnitude (', & + abs(qgwl_to_redistribute), ' m3/s) exceeds total outlet discharge (', & + global_total_outlet_discharge, ' m3/s).' + call shr_sys_flush(iulog) + endif + if (abs(qgwl_to_discharge_ratio_percent) > 5.0_r8) then call shr_sys_flush(iulog) write(iulog, *) trim(subname), & @@ -3057,18 +3069,6 @@ subroutine Rtmrun(rstwr,nlend,rdate) write(iulog, '(A,ES24.16)') trim(subname)//' [Scenario B] global_total_outlet_discharge = ', global_total_outlet_discharge write(iulog, '(A,ES24.16)') trim(subname)//' [Scenario B] correction_ratio (before bcast) = ', correction_ratio endif - - ! Check if correction ratio magnitude exceeds 100% - if (correction_ratio < -1.0_r8) then - call shr_sys_flush(iulog) - write(iulog, *) trim(subname), & - 'WARNING: Correction ratio < -100% (', correction_ratio, ').' - write(iulog, *) trim(subname), & - 'Negative runoff to ocean is unavoidable as negative qgwl magnitude (', & - abs(qgwl_to_redistribute), ' m3/s) exceeds total outlet discharge (', & - global_total_outlet_discharge, ' m3/s).' - call shr_sys_flush(iulog) - endif endif ! Broadcast the correction ratio from master to all PEs