diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSroute_GridComp/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSroute_GridComp/CMakeLists.txt index 4c19b2cb6..0644d09cb 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSroute_GridComp/CMakeLists.txt +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSroute_GridComp/CMakeLists.txt @@ -6,8 +6,12 @@ set (srcs reservoir.F90 ) -if (CMAKE_Fortran_COMPILER_ID MATCHES GNU AND CMAKE_BUILD_TYPE MATCHES Release) - set_source_files_properties(routing_model.F90 PROPERTIES COMPILE_OPTIONS ${FOPT2}) -endif () + +# routing_model.F90-specific optimization override no longer needed with gcc15. +# - reichle + borescan, 15 July 2026 +# +#if (CMAKE_Fortran_COMPILER_ID MATCHES GNU AND CMAKE_BUILD_TYPE MATCHES Release) +# set_source_files_properties(routing_model.F90 PROPERTIES COMPILE_OPTIONS ${FOPT2}) +#endif () esma_add_library (${this} SRCS ${srcs} DEPENDENCIES MAPL GEOS_LandShared ESMF::ESMF NetCDF::NetCDF_Fortran) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSroute_GridComp/routing_model.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSroute_GridComp/routing_model.F90 index 3e44ace02..87980c72a 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSroute_GridComp/routing_model.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSroute_GridComp/routing_model.F90 @@ -79,7 +79,6 @@ MODULE routing_model ! ------------------------- !**** QS = TRANSFER OF MOISTURE FROM STREAM VARIABLE TO RIVER VARIABLE [m^3/s] !**** QOUT = TRANSFER OF RIVER WATER TO THE DOWNSTREAM (DOWNRIVER) CATCHMENT [m^3/s] - SUBROUTINE RIVER_ROUTING_HYD ( & NCAT,ROUTE_DT, & Qrunf0, RRM_ALPHA_RIV, RRM_ALPHA_STR, & @@ -87,42 +86,64 @@ SUBROUTINE RIVER_ROUTING_HYD ( & Qs,Qout) IMPLICIT NONE - + INTEGER, INTENT(IN) :: NCAT,ROUTE_DT REAL, INTENT(IN), DIMENSION (NCAT) :: Qrunf0 REAL, INTENT(IN), DIMENSION (NCAT) :: RRM_ALPHA_RIV, RRM_ALPHA_STR REAL, INTENT(INOUT),DIMENSION (NCAT) :: Ws0,Wr0 REAL, INTENT(OUT), DIMENSION (NCAT) :: Qs,Qout - real, parameter :: small = 1.e-20 + real, parameter :: small = 1.e-20 real, dimension(NCAT) :: Qrunf,Ws,Wr real, dimension(NCAT) :: Qs0,ks,Ws_last - + real :: dt ! convert volume units to mass - Qrunf = Qrunf0 * rho ! m3/s -> kg/s - Ws = Ws0 * rho ! m3 -> kg - Wr = Wr0 * rho ! m3 -> kg + Qrunf = Qrunf0 * rho ! m3/s -> kg/s + Ws = Ws0 * rho ! m3 -> kg + Wr = Wr0 * rho ! m3 -> kg - dt = ROUTE_DT ! integer -> real ! If river input is too small, set alp_r to 0 + dt = ROUTE_DT ! integer -> real + + ! Update state variables: ks, Ws, and Qs - ! Update state variables: ks, Ws, and Qs - where(Qrunf<=small)Qrunf=0. ! Set runoff to zero if it's too small - Qs0=max(0.,RRM_ALPHA_STR * Ws**(1./(1.-RRM_mm))) ! Initial flow from local stream storage (kg/s) - ks =max(0.,(RRM_ALPHA_STR/(1.-RRM_mm)) * Ws**(RRM_mm/(1.-RRM_mm))) ! Flow coefficient (s^-1) - Ws_last=Ws ! Store the current water storage + where( Qrunf<=small ) Qrunf=0. ! Set runoff to zero if it is too small + + where( Ws>0.0 ) ! Avoid 0 to the power of something, does not work with gcc15 + + Qs0=max(0.,RRM_ALPHA_STR * Ws**(1./(1.-RRM_mm))) ! Initial flow from local stream storage (kg/s) + ks =max(0.,(RRM_ALPHA_STR/(1.-RRM_mm)) * Ws**(RRM_mm/(1.-RRM_mm))) ! Flow coefficient (s^-1) + + elsewhere + + Qs0 = 0. + ks = 0. + + end where + + Ws_last=Ws ! Store the current water storage where(ks>small) Ws=Ws + (Qrunf-Qs0)/ks*(1.-exp(-ks*dt)) ! Update storage (kg) where(ks<=small) Ws=Ws + (Qrunf-Qs0)*dt ! Simplified update if ks is small Ws=max(0.,Ws) ! Ensure storage is non-negative - Qs=max(0.,Qrunf-(Ws-Ws_last)/dt) ! Calculate the local stream flow (kg/s) + Qs=max(0.,Qrunf-(Ws-Ws_last)/dt) ! Calculate local stream flow (kg/s) - ! Calculate variables related to river routing: Qr0, kr + ! Update river storage and calculate river outflow Wr=Wr+Qs*dt - Qout=max(0.,RRM_ALPHA_RIV * Wr**(1./(1.-RRM_mm))) ! River flow based on water storage (kg/s) - Qout=min(Qout,Wr/dt) - Wr=max(0.,Wr-Qout*dt) + + where( Wr>0.0 ) ! Avoid 0 to the power of something, does not work with gcc15 + + Qout=max(0.,RRM_ALPHA_RIV * Wr**(1./(1.-RRM_mm))) ! River flow based on water storage (kg/s) + + elsewhere + + Qout=0.0 + + end where + + Qout=min(Qout,Wr/dt) ! Limit outflow to available river storage + Wr=max(0.,Wr-Qout*dt) ! Update river storage and keep it non-negative ! convert mass units back to volume Ws0 = Ws /rho ! kg -> m3 @@ -132,7 +153,7 @@ SUBROUTINE RIVER_ROUTING_HYD ( & RETURN - END SUBROUTINE RIVER_ROUTING_HYD + END SUBROUTINE RIVER_ROUTING_HYD ! -------------------------------------------------------------------------------------------------------