diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/CMakeLists.txt index 1e283ab8f..b44f9f5f2 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/CMakeLists.txt +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/CMakeLists.txt @@ -4,6 +4,7 @@ set (alldirs GEOSlandice_GridComp GEOSlake_GridComp GEOSland_GridComp + GEOSroute_GridComp GEOSsaltwater_GridComp ) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/CMakeLists.txt index acaec2953..fd4506218 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/CMakeLists.txt +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/CMakeLists.txt @@ -5,7 +5,6 @@ set (alldirs GEOScatch_GridComp GEOScatchCN_GridComp GEOSlana_GridComp - GEOSroute_GridComp GEOSigni_GridComp ) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOS_LandGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOS_LandGridComp.F90 index 69d73008e..3fa6c86e3 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOS_LandGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOS_LandGridComp.F90 @@ -31,7 +31,6 @@ module GEOS_LandGridCompMod use GEOS_CatchGridCompMod, only : CatchSetServices => SetServices use GEOS_CatchCNGridCompMod, only : CatchCNSetServices => SetServices use GEOS_IgniGridCompMod, only : IgniSetServices => SetServices -! use GEOS_RouteGridCompMod, only : RouteSetServices => SetServices implicit none private @@ -45,8 +44,8 @@ module GEOS_LandGridCompMod integer :: VEGDYN - integer, allocatable :: CATCH(:), ROUTE (:), CATCHCN (:) - integer :: LSM_CHOICE, RUN_ROUTE, DO_GOSWIM + integer, allocatable :: CATCH(:), CATCHCN(:) + integer :: LSM_CHOICE, DO_GOSWIM integer :: IGNI logical :: DO_FIRE_DANGER @@ -68,7 +67,7 @@ subroutine SetServices ( GC, RC ) ! !DESCRIPTION: The SetServices for the Physics GC needs to register its ! Initialize and Run. It uses the MAPL\_Generic construct for defining ! state specs and couplings among its children. In addition, it creates the -! children GCs (VegDyn, Catch, CatchCN, Route) and runs their respective SetServices. +! children GCs (VegDyn, Catch, CatchCN) and runs their respective SetServices. !EOP @@ -84,7 +83,7 @@ subroutine SetServices ( GC, RC ) character(len=ESMF_MAXSTR) :: GCName type(ESMF_Config) :: CF, SCF - integer :: NUM_CATCH + integer :: NUM_CATCH_ENS integer :: I character(len=ESMF_MAXSTR) :: TMP type(MAPL_MetaComp),pointer :: MAPL=>null() @@ -134,7 +133,7 @@ subroutine SetServices ( GC, RC ) call MAPL_GridCompSetEntryPoint ( GC, ESMF_METHOD_RUN, Run2, RC=STATUS ) VERIFY_(STATUS) - call ESMF_ConfigGetAttribute ( CF, NUM_CATCH, Label="NUM_CATCH_ENSEMBLES:", default=1, RC=STATUS) + call ESMF_ConfigGetAttribute ( CF, NUM_CATCH_ENS, Label="NUM_CATCH_ENSEMBLES:", default=1, RC=STATUS) VERIFY_(STATUS) !------------------------------------------------------------ @@ -154,7 +153,6 @@ subroutine SetServices ( GC, RC ) call MAPL_GetResource (MAPL, SURFRC, label = 'SURFRC:', default = 'GEOS_SurfaceGridComp.rc', RC=STATUS) ; VERIFY_(STATUS) SCF = ESMF_ConfigCreate(rc=status) ; VERIFY_(STATUS) call ESMF_ConfigLoadFile(SCF,SURFRC,rc=status) ; VERIFY_(STATUS) - call MAPL_GetResource (SCF, RUN_ROUTE, label='RUN_ROUTE:', DEFAULT=0, __RC__ ) call MAPL_GetResource (SCF, DO_GOSWIM, label='N_CONST_LAND4SNWALB:', DEFAULT=0, __RC__ ) call MAPL_GetResource (SCF, DO_FIRE_DANGER, label='FIRE_DANGER:', DEFAULT=.false., __RC__ ) call ESMF_ConfigDestroy (SCF, __RC__) @@ -163,13 +161,13 @@ subroutine SetServices ( GC, RC ) CASE (1) - allocate (CATCH(NUM_CATCH), stat=status) + allocate (CATCH(NUM_CATCH_ENS), stat=status) VERIFY_(STATUS) - if (NUM_CATCH == 1) then + if (NUM_CATCH_ENS == 1) then CATCH(1) = MAPL_AddChild(GC, NAME='CATCH'//trim(tmp), SS=CatchSetServices, RC=STATUS) VERIFY_(STATUS) else - do I = 1, NUM_CATCH + do I = 1, NUM_CATCH_ENS WRITE(TMP,'(I3.3)') I GCName = 'ens' // trim(TMP) // ':CATCH' CATCH(I) = MAPL_AddChild(GC, NAME=GCName, SS=CatchSetServices, RC=STATUS) @@ -179,13 +177,13 @@ subroutine SetServices ( GC, RC ) CASE (2,3) - allocate (CATCHCN(NUM_CATCH), stat=status) + allocate (CATCHCN(NUM_CATCH_ENS), stat=status) VERIFY_(STATUS) - if (NUM_CATCH == 1) then + if (NUM_CATCH_ENS == 1) then CATCHCN(1) = MAPL_AddChild(GC, NAME='CATCHCN'//trim(tmp), SS=CatchCNSetServices, RC=STATUS) VERIFY_(STATUS) else - do I = 1, NUM_CATCH + do I = 1, NUM_CATCH_ENS WRITE(TMP,'(I3.3)') I GCName = 'ens' // trim(TMP) // ':CATCHCN' CATCHCN(I) = MAPL_AddChild(GC, NAME=GCName, SS=CatchCNSetServices, RC=STATUS) @@ -195,20 +193,6 @@ subroutine SetServices ( GC, RC ) END SELECT -! IF(RUN_ROUTE == 1) THEN -! if (NUM_CATCH == 1) then -! ROUTE(1) = MAPL_AddChild(GC, NAME='ROUTE', SS=RouteSetServices, RC=STATUS) -! VERIFY_(STATUS) -! else -! do I = 1, NUM_CATCH -! WRITE(TMP,'(I3.3)') I -! GCName = 'ens' // trim(TMP) // ':ROUTE' -! ROUTE(I) = MAPL_AddChild(GC, NAME=GCName, SS=RouteSetServices, RC=STATUS) -! VERIFY_(STATUS) -! end do -! end if -! ENDIF - if (DO_FIRE_DANGER) then IGNI = MAPL_AddChild(GC, NAME='IGNI'//trim(tmp), SS=IgniSetServices, RC=STATUS) VERIFY_(STATUS) @@ -1377,14 +1361,6 @@ subroutine SetServices ( GC, RC ) CHILD_ID = VEGDYN,& RC=STATUS ) VERIFY_(STATUS) -! IF(RUN_ROUTE == 1) THEN -! call MAPL_AddExportSpec ( GC, & -! SHORT_NAME = 'QOUTFLOW', & -! CHILD_ID = ROUTE(1), & -! RC=STATUS ) -! VERIFY_(STATUS) -! ENDIF - if (DO_FIRE_DANGER) then call MAPL_AddExportSpec ( GC, SHORT_NAME = 'FFMC', CHILD_ID = IGNI, __RC__ ) @@ -1426,7 +1402,7 @@ subroutine SetServices ( GC, RC ) ! !CONNECTIONS: - DO I = 1, NUM_CATCH + DO I = 1, NUM_CATCH_ENS SELECT CASE (LSM_CHOICE) @@ -1453,17 +1429,6 @@ subroutine SetServices ( GC, RC ) VERIFY_(STATUS) end if -! IF(RUN_ROUTE == 1) THEN -! call MAPL_AddConnectivity ( & -! GC ,& -! SHORT_NAME = (/'RUNOFF '/) ,& -! SRC_ID = CATCH(I) ,& -! DST_ID = ROUTE(I) ,& -! -! RC=STATUS ) -! VERIFY_(STATUS) -! ENDIF - CASE (2,3) call MAPL_AddConnectivity ( & GC , & @@ -1486,19 +1451,9 @@ subroutine SetServices ( GC, RC ) VERIFY_(STATUS) end if -! IF(RUN_ROUTE == 1) THEN -! call MAPL_AddConnectivity ( & -! GC ,& -! SHORT_NAME = (/'RUNOFF '/) ,& -! SRC_ID = CATCHCN(I) ,& -! DST_ID = ROUTE(I) ,& -! -! RC=STATUS ) -! VERIFY_(STATUS) -! ENDIF END SELECT END DO - + call MAPL_TimerAdd(GC, name="INITIALIZE" ,RC=STATUS) VERIFY_(STATUS) @@ -1714,7 +1669,7 @@ subroutine Run2(GC, IMPORT, EXPORT, CLOCK, RC ) type (ESMF_State), pointer :: GEX(:) character(len=ESMF_MAXSTR),pointer :: GCnames(:) - integer :: I + integer :: I, phase !============================================================================= diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/CMakeLists.txt deleted file mode 100644 index 8a502e3e7..000000000 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/CMakeLists.txt +++ /dev/null @@ -1,10 +0,0 @@ -esma_set_this () - -set (srcs - #GEOS_RouteGridComp.F90 - routing_model.F90 - ) - -esma_add_library (${this} SRCS ${srcs} DEPENDENCIES MAPL ESMF::ESMF NetCDF::NetCDF_Fortran) - -install(PROGRAMS build_rivernetwork.py DESTINATION bin) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/GEOS_RouteGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/GEOS_RouteGridComp.F90 deleted file mode 100644 index 4a28c0da2..000000000 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/GEOS_RouteGridComp.F90 +++ /dev/null @@ -1,1116 +0,0 @@ - -#include "MAPL_Generic.h" - -#ifndef RUN_FOR_REAL -#define MAPL_DimsCatchOnly MAPL_DimsTileOnly -#endif - -!============================================================================= -module GEOS_RouteGridCompMod - -!BOP -! !MODULE: GEOS_Route -- child to the "Land" gridded component. - -! !DESCRIPTION: -! {\tt GEOS\_Route} is a gridded component to route total runoff produced in -! {\tt GEOS\_Catch} (RUNOFF in {\tt GEOScatch\_GridComp} or {\tt GEOScatchCN\_GridComp}) through a 290,188 -! watersheds globally (excluding submerged catchments (watersheds) from the global list of 291,284. -! All of its calculations are done on Pfafstetter watershed space. {\tt GEOS\_Route} has no children. \\ -! -! IMPORTS : RUNOFF \\ -! INTERNALS : AREACAT, LENGSC2, DNSTR, WSTREAM, WRIVER, LRIVERMOUTH, ORIVERMOUTH \\ -! EXPORTS : QSFLOW, QINFLOW, QOUTFLOW \\ - -! !USES: - - use ESMF - use MAPL_Mod - use MAPL_ConstantsMod - use ROUTING_MODEL, ONLY: & - river_routing, ROUTE_DT -#if 0 - USE catch_constants, ONLY: & - N_CatG => N_Pfaf_Catchs -#endif - - implicit none - integer, parameter :: N_CatG = 291284 - private - - type T_RROUTE_STATE - private - type (ESMF_RouteHandle) :: routeHandle - type (ESMF_Field) :: field - integer :: nTiles - integer :: comm - integer :: nDes - integer :: myPe - integer :: minCatch - integer :: maxCatch - integer, pointer :: pfaf(:) => NULL() - real, pointer :: tile_area(:) => NULL() - end type T_RROUTE_STATE - - ! Wrapper for extracting internal state - ! ------------------------------------- - type RROUTE_WRAP - type (T_RROUTE_STATE), pointer :: PTR => null() - end type RROUTE_WRAP - - include "mpif.h" - - -! !PUBLIC MEMBER FUNCTIONS: - - public SetServices - -!EOP - -contains - -!BOP - -! !IROUTINE: SetServices -- Sets ESMF services for this component - -! !INTERFACE: - - subroutine SetServices ( GC, RC ) - -! !ARGUMENTS: - - type(ESMF_GridComp), intent(INOUT) :: GC ! gridded component - integer, optional :: RC ! return code - -! !DESCRIPTION: This version uses the MAPL\_GenericSetServices. This function sets -! the Initialize and Finalize services, as well as allocating -! our instance of a generic state and putting it in the -! gridded component (GC). Here we only need to set the run method and -! add the state variable specifications (also generic) to our instance -! of the generic state. This is the way our true state variables get into -! the ESMF\_State INTERNAL, which is in the MAPL\_MetaComp. - -!EOP - -!============================================================================= -! -! ErrLog Variables - - - character(len=ESMF_MAXSTR) :: IAm - integer :: STATUS - character(len=ESMF_MAXSTR) :: COMP_NAME - -! Local derived type aliases - - type (ESMF_Config ) :: CF - - type (T_RROUTE_STATE), pointer :: route_internal_state => null() - type (RROUTE_wrap) :: wrap - - integer :: RUN_DT - real :: DT - -!============================================================================= - -! Begin... - -!------------------------------------------------------------ -! Get my name and set-up traceback handle -!------------------------------------------------------------ - - call ESMF_GridCompGet(GC ,& - NAME=COMP_NAME ,& - RC=STATUS ) - - VERIFY_(STATUS) - - Iam = trim(COMP_NAME) // 'SetServices' - -! ----------------------------------------------------------- -! Set the Initialize and Run entry points -! ----------------------------------------------------------- - - call MAPL_GridCompSetEntryPoint ( GC, ESMF_METHOD_INITIALIZE, Initialize, RC=STATUS ) - VERIFY_(STATUS) - call MAPL_GridCompSetEntryPoint (GC, ESMF_METHOD_RUN, Run, RC=STATUS) - VERIFY_(STATUS) - -!------------------------------------------------------------ -! Set generic final method -!------------------------------------------------------------ - - -! ----------------------------------------------------------- -! Get the configuration -! ----------------------------------------------------------- -! - call ESMF_GridCompGet( GC, CONFIG = CF, RC=STATUS ) - VERIFY_(STATUS) -! -! ----------------------------------------------------------- -! Get the intervals -! ----------------------------------------------------------- -! - call ESMF_ConfigGetAttribute (CF, DT ,& - Label="RUN_DT:" ,& - RC=STATUS) - - VERIFY_(STATUS) - - RUN_DT = nint(DT) - -! ----------------------------------------------------------- -! At the moment, this will refresh when the land parent -! needs to refresh. - - call ESMF_ConfigGetAttribute ( CF, DT, Label=trim(COMP_NAME)//"_DT:", & - default=DT, RC=STATUS) - VERIFY_(STATUS) - -! ----------------------------------------------------------- -! Set the state variable specs. -! ----------------------------------------------------------- - -!BOS - -! ----------------------------------------------------------- -! Import States -! ----------------------------------------------------------- - - call MAPL_AddImportSpec(GC, & - LONG_NAME = 'runoff_total_flux' ,& - UNITS = 'kg m-2 s-1' ,& - SHORT_NAME = 'RUNOFF' ,& - DIMS = MAPL_DimsTileOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RC=STATUS ) - VERIFY_(STATUS) - -! ----------------------------------------------------------- -! INTERNAL STATE -! ----------------------------------------------------------- - - call MAPL_AddInternalSpec(GC ,& - LONG_NAME = 'area_of_catchment' ,& - UNITS = 'km+2' ,& - SHORT_NAME = 'AREACAT' ,& - DIMS = MAPL_DimsCatchOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RESTART = MAPL_RestartRequired ,& - RC=STATUS ) - - call MAPL_AddInternalSpec(GC ,& - LONG_NAME = 'length_of_channel_segment',& - UNITS = 'km+2' ,& - SHORT_NAME = 'LENGSC' ,& - DIMS = MAPL_DimsCatchOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RESTART = MAPL_RestartRequired ,& - RC=STATUS ) - - call MAPL_AddInternalSpec(GC ,& - LONG_NAME = 'index_of_downtream_catchment',& - UNITS = '1' ,& - SHORT_NAME = 'DNSTR' ,& - DIMS = MAPL_DimsCatchOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RESTART = MAPL_RestartRequired ,& - RC=STATUS ) - - call MAPL_AddInternalSpec(GC ,& - LONG_NAME = 'volume_of_water_in_local_stream',& - UNITS = 'm+3' ,& - SHORT_NAME = 'WSTREAM' ,& - DIMS = MAPL_DimsCatchOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RESTART = MAPL_RestartRequired ,& - RC=STATUS ) - - call MAPL_AddInternalSpec(GC ,& - LONG_NAME = 'volume_of_water_in_river' ,& - UNITS = 'm+3' ,& - SHORT_NAME = 'WRIVER' ,& - DIMS = MAPL_DimsCatchOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RESTART = MAPL_RestartRequired ,& - RC=STATUS ) - - call MAPL_AddInternalSpec(GC ,& - LONG_NAME = 'TileID_of_the_lake_tile_at_the_river_mouth' ,& - UNITS = '1' ,& - SHORT_NAME = 'LRIVERMOUTH' ,& - DIMS = MAPL_DimsCatchOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RESTART = MAPL_RestartRequired ,& - RC=STATUS ) - - call MAPL_AddInternalSpec(GC ,& - LONG_NAME = 'TileID_of_the_ocean_tile_at_the_river_mouth' ,& - UNITS = '1' ,& - SHORT_NAME = 'ORIVERMOUTH' ,& - DIMS = MAPL_DimsCatchOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RESTART = MAPL_RestartRequired ,& - RC=STATUS ) -! ----------------------------------------------------------- -! EXPORT STATE: -! ----------------------------------------------------------- - - call MAPL_AddExportSpec(GC, & - LONG_NAME = 'transfer_of_moisture_from_stream_variable_to_river_variable' ,& - UNITS = 'm+3 s-1' ,& - SHORT_NAME = 'QSFLOW' ,& - DIMS = MAPL_DimsCatchOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RC=STATUS ) - - call MAPL_AddExportSpec(GC, & - LONG_NAME = 'transfer_of_river_water_from_upstream_catchments' ,& - UNITS = 'm+3 s-1' ,& - SHORT_NAME = 'QINFLOW' ,& - DIMS = MAPL_DimsCatchOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RC=STATUS ) - - call MAPL_AddExportSpec(GC, & - LONG_NAME = 'transfer_of_river_water_to_downstream_catchments' ,& - UNITS = 'm+3 s-1' ,& - SHORT_NAME = 'QOUTFLOW' ,& - DIMS = MAPL_DimsCatchOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RC=STATUS ) - -!EOS - - call MAPL_TimerAdd(GC, name="RUN" ,RC=STATUS) - VERIFY_(STATUS) - call MAPL_TimerAdd(GC, name="-RRM" ,RC=STATUS) - VERIFY_(STATUS) - - - -! Allocate this instance of the internal state and put it in wrapper. -! ------------------------------------------------------------------- - - allocate( route_internal_state, stat=status ) - VERIFY_(STATUS) - wrap%ptr => route_internal_state - -! Save pointer to the wrapped internal state in the GC -! ---------------------------------------------------- - - call ESMF_UserCompSetInternalState ( GC, 'RiverRoute_state',wrap,status ) - VERIFY_(STATUS) - -! Clocks -!------- - - call MAPL_TimerAdd(GC, name="INITIALIZE" ,RC=STATUS) - VERIFY_(STATUS) - call MAPL_TimerAdd(GC, name="RUN" ,RC=STATUS) - VERIFY_(STATUS) - -! All done -!--------- - - call MAPL_GenericSetServices(GC, RC=STATUS) - VERIFY_(STATUS) - - RETURN_(ESMF_SUCCESS) - - end subroutine SetServices - -! ----------------------------------------------------------- -! INITIALIZE -- Initialize method for the route component -! ----------------------------------------------------------- - - subroutine INITIALIZE (GC,IMPORT, EXPORT, CLOCK, RC ) - -! ----------------------------------------------------------- -! !ARGUMENTS: -! ----------------------------------------------------------- - - type(ESMF_GridComp), intent(inout) :: GC - type(ESMF_State), intent(inout) :: IMPORT - type(ESMF_State), intent(inout) :: EXPORT - type(ESMF_Clock), intent(inout) :: CLOCK - integer, optional, intent( out) :: RC - -! ----------------------------------------------------------- -! ErrLog Variables -! ----------------------------------------------------------- - - character(len=ESMF_MAXSTR) :: IAm="Initialize" - integer :: STATUS - character(len=ESMF_MAXSTR) :: COMP_NAME - -! ----------------------------------------------------------- -! Locals -! ----------------------------------------------------------- - - type (ESMF_VM) :: VM - integer :: comm - integer :: nDEs - integer :: myPE - integer :: beforeMe, minCatch, maxCatch, pf, i - integer :: ntiles, nt_global - - type(ESMF_Grid) :: tileGrid - type(ESMF_Grid) :: newTileGrid - type(ESMF_Grid) :: catchGrid - type(ESMF_DistGrid) :: distGrid - type(ESMF_Field) :: field, field0 - - type(MAPL_MetaComp), pointer :: MAPL - type(MAPL_LocStream) :: locstream - - integer, pointer :: ims(:) => NULL() - integer, pointer :: pfaf(:) => NULL() - integer, pointer :: arbSeq(:) => NULL() - integer, allocatable :: arbIndex(:,:) - real, pointer :: tile_area_src(:) => NULL() - real, pointer :: tile_area(:) => NULL() - real, pointer :: ptr2(:) => NULL() - - type (T_RROUTE_STATE), pointer :: route => null() - type (RROUTE_wrap) :: wrap - - ! ------------------ - ! begin - - call ESMF_UserCompGetInternalState ( GC, 'RiverRoute_state',wrap,status ) - VERIFY_(STATUS) - - route => wrap%ptr - - ! get vm - ! extract comm - call ESMF_VMGetCurrent(VM, RC=STATUS) - VERIFY_(STATUS) - call ESMF_VMGet (VM, mpiCommunicator =comm, RC=STATUS) - VERIFY_(STATUS) - call ESMF_VMGet (VM, localpet=MYPE, petcount=nDEs, RC=STATUS) - VERIFY_(STATUS) - - call MAPL_GetObjectFromGC ( GC, MAPL, RC=STATUS) - VERIFY_(STATUS) - - route%comm = comm - route%ndes = ndes - route%mype = mype - - ! define minCatch, maxCatch - call MAPL_DecomposeDim ( n_catg,ims,ndes ) ! ims(mype+1) gives the size of my partition - ! myPE is 0-based! - beforeMe = sum(ims(1:mype)) - minCatch = beforeMe + 1 - maxCatch = beforeMe + ims(myPe+1) - - ! get LocStream - call MAPL_Get(MAPL, LocStream = locstream, RC=status) - VERIFY_(STATUS) - ! extract Pfaf (TILEI on the "other" grid) - call MAPL_LocStreamGet(locstream, tilei=pfaf, OnAttachedGrid=.false., & - tileGrid=tilegrid, nt_global=nt_global, RC=status) - VERIFY_(STATUS) - - ! exchange Pfaf across PEs - - ntiles = 0 - !loop over total_n_tiles - do i = 1, nt_global - pf = pfaf(i) - if (pf >= minCatch .and. pf <= maxCatch) then ! I want this! - ntiles = ntiles+1 - !realloc if needed - arbSeq(ntiles) = i - end if - end do ! global tile loop - - distgrid = ESMF_DistGridCreate(arbSeqIndexList=arbSeq, rc=status) - VERIFY_(STATUS) - - newTileGRID = ESMF_GridEmptyCreate(rc=status) - VERIFY_(STATUS) - - allocate(arbIndex(nTiles,1), stat=status) - VERIFY_(STATUS) - - arbIndex(:,1) = arbSeq - - call ESMF_GridSet(newTileGrid, & - name='redist_tile_grid_for_'//trim(COMP_NAME), & - distgrid=distgrid, & - gridMemLBound=(/1/), & - indexFlag=ESMF_INDEX_USER, & - distDim = (/1/), & - localArbIndexCount=ntiles, & - localArbIndex=arbIndex, & - minIndex=(/1/), & - maxIndex=(/NT_GLOBAL/), & - rc=status) - VERIFY_(STATUS) - - deallocate(arbIndex) - - call ESMF_GridCommit(newTileGrid, rc=status) - VERIFY_(STATUS) - - - ! now create a "catch" grid to be the "native" grid for this component - distgrid = ESMF_DistGridCreate(arbSeqIndexList=(/minCatch:maxCatch/), & - rc=status) - VERIFY_(STATUS) - - catchGRID = ESMF_GridEmptyCreate(rc=status) - VERIFY_(STATUS) - - allocate(arbIndex(ims(myPE+1),1), stat=status) - VERIFY_(STATUS) - - arbIndex(:,1) = (/minCatch:maxCatch/) - - call ESMF_GridSet(catchGrid, & - name='catch_grid_for_'//trim(COMP_NAME), & - distgrid=distgrid, & - gridMemLBound=(/1/), & - indexFlag=ESMF_INDEX_USER, & - distDim = (/1/), & - localArbIndexCount=ims(myPE+1), & - localArbIndex=arbIndex, & - minIndex=(/1/), & - maxIndex=(/N_CatG/), & - rc=status) - VERIFY_(STATUS) - - deallocate(arbIndex) - - call ESMF_GridCommit(catchGrid, rc=status) - VERIFY_(STATUS) - - call ESMF_GridCompSet(gc, grid=catchGrid, RC=status) - VERIFY_(STATUS) - - call MAPL_LocStreamGet(locstream, TILEAREA = tile_area_src, RC=status) - VERIFY_(STATUS) - - field0 = ESMF_FieldCreate(grid=tilegrid, datacopyflag=ESMF_DATACOPY_VALUE, & - farrayPtr=tile_area_src, name='TILE_AREA_SRC', RC=STATUS) - VERIFY_(STATUS) - ! create field on the "new" tile grid - allocate(tile_area(ntiles), stat=status) - VERIFY_(STATUS) - field = ESMF_FieldCreate(grid=newtilegrid, datacopyflag=ESMF_DATACOPY_VALUE, & - farrayPtr=tile_area, name='TILE_AREA', RC=STATUS) - VERIFY_(STATUS) - - ! create routehandle - call ESMF_FieldRedistStore(srcField=field0, dstField=field, & - routehandle=route%routehandle, rc=status) - VERIFY_(STATUS) - - ! redist tile_area - call ESMF_FieldRedist(srcField=FIELD0, dstField=FIELD, & - routehandle=route%routehandle, rc=status) - VERIFY_(STATUS) - - call ESMF_FieldDestroy(field, rc=status) - VERIFY_(STATUS) - call ESMF_FieldDestroy(field0, rc=status) - VERIFY_(STATUS) - - ! redist pfaf (NOTE: me might need a second routehandle for integers) - - route%pfaf => arbSeq - route%ntiles = ntiles - route%minCatch = minCatch - route%maxCatch = maxCatch - - allocate(ptr2(ntiles), stat=status) - VERIFY_(STATUS) - route%field = ESMF_FieldCreate(grid=newtilegrid, datacopyflag=ESMF_DATACOPY_VALUE, & - farrayPtr=ptr2, name='RUNOFF', RC=STATUS) - VERIFY_(STATUS) - - deallocate(ims) - call MAPL_GenericInitialize ( GC, import, export, clock, rc=status ) - VERIFY_(STATUS) - - RETURN_(ESMF_SUCCESS) - end subroutine INITIALIZE - -! ----------------------------------------------------------- -! RUN -- Run method for the route component -! ----------------------------------------------------------- - - subroutine RUN (GC,IMPORT, EXPORT, CLOCK, RC ) - -! ----------------------------------------------------------- -! !ARGUMENTS: -! ----------------------------------------------------------- - - type(ESMF_GridComp), intent(inout) :: GC - type(ESMF_State), intent(inout) :: IMPORT - type(ESMF_State), intent(inout) :: EXPORT - type(ESMF_Clock), intent(inout) :: CLOCK - integer, optional, intent( out) :: RC - -! ----------------------------------------------------------- -! ErrLog Variables -! ----------------------------------------------------------- - - character(len=ESMF_MAXSTR) :: IAm="Run" - integer :: STATUS - character(len=ESMF_MAXSTR) :: COMP_NAME - -! ----------------------------------------------------------- -! Locals -! ----------------------------------------------------------- - - type (MAPL_MetaComp), pointer :: MAPL - type (ESMF_State ) :: INTERNAL -! type(ESMF_Alarm) :: ALARM - type (ESMF_Config ) :: CF - type(ESMF_VM) :: VM - -! ----------------------------------------------------- -! IMPORT pointers -! ----------------------------------------------------- - - real, dimension(:), pointer :: RUNOFF - -! ----------------------------------------------------- -! INTERNAL pointers -! ----------------------------------------------------- - - real, dimension(:), pointer :: AREACAT - real, dimension(:), pointer :: LENGSC - real, dimension(:), pointer :: DNSTR - real, dimension(:), pointer :: WSTREAM - real, dimension(:), pointer :: WRIVER - real, dimension(:), pointer :: LRIVERMOUTH - real, dimension(:), pointer :: ORIVERMOUTH - -! ----------------------------------------------------- -! EXPORT pointers -! ----------------------------------------------------- - - real, dimension(:), pointer :: QSFLOW - real, dimension(:), pointer :: QINFLOW - real, dimension(:), pointer :: QOUTFLOW - -! Time attributes and placeholders - -! type(ESMF_Time) :: CURRENT_TIME - -! Others - - type(ESMF_Grid) :: TILEGRID - type (MAPL_LocStream) :: LOCSTREAM - integer :: NTILES, N_CatL, N_CYC - logical, save :: FirstTime=.true. - real, pointer, dimension(:) :: tile_area - integer, pointer, dimension(:) :: pfaf_code - - INTEGER, DIMENSION(:,:), POINTER, SAVE :: AllActive,DstCatchID - INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: srcProcsID, LocDstCatchID - integer, dimension (:),allocatable, SAVE :: GlbActive - INTEGER, SAVE :: N_Active, ThisCycle - INTEGER :: Local_Min, Local_Max - integer :: K, N, I, req - REAL :: mm2m3, rbuff, HEARTBEAT - REAL, ALLOCATABLE, DIMENSION(:) :: RUNOFF_CATCH, RUNOFF_ACT,AREACAT_ACT,& - LENGSC_ACT, WSTREAM_ACT,WRIVER_ACT, QSFLOW_ACT,QOUTFLOW_ACT, runoff_save - INTEGER, ALLOCATABLE, DIMENSION(:) :: tmp_index - type(ESMF_Field) :: runoff_src - - integer :: ndes, mype - type (T_RROUTE_STATE), pointer :: route => null() - type (RROUTE_wrap) :: wrap - - ! ------------------ - ! begin - - call ESMF_UserCompGetInternalState ( GC, 'RiverRoute_state',wrap,status ) - VERIFY_(STATUS) - - route => wrap%ptr - -! Get the target components name and set-up traceback handle. -! ----------------------------------------------------------- - - call ESMF_GridCompGet(GC, name=COMP_NAME, CONFIG=CF, RC=STATUS ) - VERIFY_(STATUS) - - Iam = trim(COMP_NAME) // "RUN" - -! Get my internal MAPL_Generic state -! ----------------------------------------------------------- - - call MAPL_GetObjectFromGC(GC, MAPL, STATUS) - VERIFY_(STATUS) - - call MAPL_Get(MAPL, HEARTBEAT = HEARTBEAT, RC=STATUS) - VERIFY_(STATUS) - -! Start timers -! ------------ - - call MAPL_TimerOn(MAPL,"RUN") - -! Get parameters from generic state -! --------------------------------- - - call MAPL_Get(MAPL, INTERNAL_ESMF_STATE=INTERNAL, RC=STATUS) - VERIFY_(STATUS) - -! get pointers to inputs variables -! ---------------------------------- - - ! get the field from IMPORT - call ESMF_StateGet(IMPORT, 'RUNOFF', field=runoff_src, RC=STATUS) - VERIFY_(STATUS) - - ! redist RunOff - call ESMF_FieldRedist(srcField=runoff_src, dstField=route%field, & - routehandle=route%routehandle, rc=status) - VERIFY_(STATUS) - - call ESMF_FieldGet(route%field, farrayPtr=RUNOFF, rc=status) - VERIFY_(STATUS) - - pfaf_code => route%pfaf - tile_area => route%tile_area - -! get pointers to internal variables -! ---------------------------------- - - call MAPL_GetPointer(INTERNAL, AREACAT , 'AREACAT', RC=STATUS) - VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL, LENGSC , 'LENGSC', RC=STATUS) - VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL, DNSTR , 'DNSTR' , RC=STATUS) - VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL, WSTREAM , 'WSTREAM', RC=STATUS) - VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL, WRIVER , 'WRIVER' , RC=STATUS) - VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL, LRIVERMOUTH, 'LRIVERMOUTH' , RC=STATUS) - VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL, ORIVERMOUTH, 'ORIVERMOUTH' , RC=STATUS) - VERIFY_(STATUS) - -! get pointers to EXPORTS -! ----------------------- - - call MAPL_GetPointer(EXPORT, QSFLOW, 'QSFLOW' , RC=STATUS) - VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT, QINFLOW, 'QINFLOW' , RC=STATUS) - VERIFY_(STATUS) - call MAPL_GetPointer(EXPORT, QOUTFLOW, 'QOUTFLOW', RC=STATUS) - VERIFY_(STATUS) - - call MAPL_Get(MAPL, LocStream=LOCSTREAM, RC=STATUS) - VERIFY_(STATUS) - call MAPL_LocStreamGet(LOCSTREAM, TILEGRID=TILEGRID, RC=STATUS) - VERIFY_(STATUS) - - call MAPL_TimerOn ( MAPL, "-RRM" ) - - call MAPL_LocStreamGet(LocStream, NT_LOCAL=NTILES, RC=STATUS ) - N_CatL = size(AREACAT) - -!@@ ALLOCATE (pfaf_code (1:NTILES)) ! 9th_coulumn_in_TILFILE - - ! NOTES : - !Need below area and pfaf_index from the .til file (Maybe, they are already in LocStream) - ! - ! TILFILE: /discover/nobackup/smahanam/bcs/Heracles-4_3/Heracles-4_3_MERRA-3/CF0090x6C_DE1440xPE0720/CF0090x6C_DE1440xPE0720-Pfafstetter.til - ! The 8-line header is followed by 1061481 number of rows. - ! do n = 1,475330 - ! read (10,*)type,area, longitude, latitude, ig, jg, cell_frac, integer, & - ! pfaf_code, pfaf_index, pfaf_frac - ! end do - ! - ! where for each tile: - ! (1) type [-] tile type (100-land; 19-lakes; 20-ice) - ! (2) area [x EarthRadius^2 km2] tile area - ! (3) longitude [degree] longitude at the centroid of the tile - ! (4) latitude [degree] latitude at the centroid of the tile - ! (5) ig [-] i-index of the AGCM grid cell where the tile is located - ! (6) jg [-] j-index of the AGCM grid cell where the tile is located - ! (7) cell_frac [-] fraction of the AGCM grid cell - ! (8) integer some integer that matters only for OGCM tiles, I suppose. - ! (9) pfaf_code [-] catchment index (1-291284) after sorting Pfafstetter codes in ascending order - ! (10) pfaf_index[-] catchment index (1-290188) after sorting Pfafstetter codes - ! and removing submerged in ascending order - ! (11) pfaf_frac [-] fraction of the pfafstetter catchment - - !call MAPL_LocStreamGet(LocStream, 9th_coulumn_in_TILFILE=pfaf_code, RC=STATUS ) - - Local_Min = route%minCatch - Local_Max = route%maxCatch - - FIRST_TIME : IF (FirstTime) THEN - - ! Pfafstetter catchment Domain Decomposition : - ! -------------------------------------------- - - ! AllActive : Processor(s) where the catchment is active (identical in any processor). - ! srcProcsID : For all active catchments anywhere tells which processor is the principal owner of the catchment (identical in any processor). - ! DstCatchID : 2-D array contains downstream catchID and downstream processor (identical in any processor) - ! LocDstCatchID : Downstream catchID when for catchments that are local to the processor. - - ndes = route%ndes - mype = route%mype - allocate (AllActive (1:N_CatG, 1: nDEs)) - allocate (DstCatchID(1:N_CatG, 1: nDEs)) - allocate (srcProcsID (1:N_CatG )) - allocate (LocDstCatchID(1:N_CatG )) - - AllActive = -9999 - srcProcsID = -9999 - DstCatchID = -9999 - LocDstCatchID = NINT(DNSTR) - - call InitializeRiverRouting(MYPE, nDEs, MAPL_am_I_root(vm),pfaf_code, & - AllActive, DstCatchID, srcProcsID, LocDstCatchID, rc=STATUS) - - VERIFY_(STATUS) - - N_Active = count (srcProcsID == MYPE) - - allocate (GlbActive(1 : N_Active)) - allocate (tmp_index(1 : N_CatG )) - - forall (N=1:N_CatG) tmp_index(N) = N - - GlbActive = pack (tmp_index, mask = (srcProcsID == MYPE)) - - ! Initialize the cycle counter and sum (runoff) - - allocate (runoff_save (1:NTILES)) - - runoff_save = 0. - ThisCycle = 1 - - FirstTime = .false. - - deallocate (tmp_index) - - ENDIF FIRST_TIME - - ! For efficiency, the time step to call the river routing model is set at ROUTE_DT - - N_CYC = ROUTE_DT/HEARTBEAT - - RUN_MODEL : if (ThisCycle == N_CYC) then - - runoff_save = runoff_save + runoff/real (N_CYC) - - ! Here we aggreagate GEOS_Catch/GEOS_CatchCN produced RUNOFF from TILES to CATCHMENTS - ! Everything is local to the parallel block. Units: RUNOFF [kg m-2 s-1], - ! RUNOFF_CATCH [m3 s-1] - ! ----------------------------------------------------------------------------------- - - ! Unit conversion - - mm2m3 = MAPL_RADIUS * MAPL_RADIUS / 1000. - - ALLOCATE (RUNOFF_CATCH(1:N_CatG)) - - RUNOFF_CATCH = 0. - - DO N = 1, NTILES - RUNOFF_CATCH (pfaf_code(n)) = RUNOFF_CATCH (pfaf_code(n)) + mm2m3 * RUNOFF_SAVE (N) * TILE_AREA (N) - END DO - - ! Inter-processor communication 1 - ! For catchment-tiles that contribute to the main catchment in some other processor, - ! send runoff to the corresponding srcProcsID(N) - ! ----------------------------------------------------------------------------------- - - do N = Local_Min, Local_Max - - if ((AllActive (N,MYPE+1) > 0).and.(srcProcsID(N) /= MYPE)) then - - rbuff = RUNOFF_CATCH (N) - - call MPI_ISend(rbuff,1,MPI_real,srcProcsID(N),999,MPI_COMM_WORLD,req,status) - call MPI_WAIT (req ,MPI_STATUS_IGNORE,status) - - RUNOFF_CATCH (N) = 0. - - else - - if(srcProcsID(N) == MYPE) then - - do i = 1,nDEs - if((i-1 /= MYPE).and.(AllActive (N,i) > 0)) then - - call MPI_RECV(rbuff,1,MPI_real,i-1,999,MPI_COMM_WORLD,MPI_STATUS_IGNORE,status) - RUNOFF_CATCH (N) = RUNOFF_CATCH (N) + rbuff - - endif - end do - endif - endif - end do - - ! Now compress and create subsets of arrays that only contain active catchments - ! in the local processor - ! ----------------------------------------------------------------------------- - - if(allocated (LENGSC_ACT ) .eqv. .false.) allocate (LENGSC_ACT (1:N_Active)) - if(allocated (AREACAT_ACT ) .eqv. .false.) allocate (AREACAT_ACT (1:N_Active)) - if(allocated (WSTREAM_ACT ) .eqv. .false.) allocate (WSTREAM_ACT (1:N_Active)) - if(allocated (WRIVER_ACT ) .eqv. .false.) allocate (WRIVER_ACT (1:N_Active)) - if(allocated (QSFLOW_ACT ) .eqv. .false.) allocate (QSFLOW_ACT (1:N_Active)) - if(allocated (QOUTFLOW_ACT) .eqv. .false.) allocate (QOUTFLOW_ACT(1:N_Active)) - if(allocated (RUNOFF_ACT ) .eqv. .false.) allocate (RUNOFF_ACT (1:N_Active)) - - DO N = 1, size (GlbActive) - - I = GlbActive (N) - RUNOFF_ACT (N) = RUNOFF_CATCH (I) - - I = GlbActive (N) - Local_Min + 1 - WSTREAM_ACT (N) = WSTREAM (I) - WRIVER_ACT (N) = WRIVER (I) - LENGSC_ACT (N) = LENGSC (I) - AREACAT_ACT (N) = AREACAT (I) - - END DO - - QSFLOW_ACT = 0. - QOUTFLOW_ACT = 0. - QSFLOW = 0. - QOUTFLOW = 0. - QINFLOW = 0. - - ! Call river_routing_model - ! ------------------------ - - CALL RIVER_ROUTING (N_Active, RUNOFF_ACT,AREACAT_ACT,LENGSC_ACT, & - WSTREAM_ACT,WRIVER_ACT, QSFLOW_ACT,QOUTFLOW_ACT) - - DO N = 1, size (GlbActive) - - I = GlbActive (N) - Local_Min + 1 - - WSTREAM (I) = WSTREAM_ACT (N) - WRIVER (I) = WRIVER_ACT (N) - QSFLOW (I) = QSFLOW_ACT (N) - QOUTFLOW(I) = QOUTFLOW_ACT(N) - - if (LocDstCatchID (GlbActive (N)) == GlbActive (N)) then - - ! This catchment drains to the ocean, lake or a sink - ! if(ORIVERMOUTH(... ) > 0) send QOUTFLOW(I) [m3/s] to ORIVERMOUTH(N) th ocean tile - ! if(LRIVERMOUTH(... ) > 0) send QOUTFLOW(I) [m3/s] to LRIVERMOUTH(N) th lake tile - - endif - END DO - - ! Inter-processor communication-2 - ! Update down stream catchments - ! ------------------------------- - - do N = 1,N_CatG - - if ((srcProcsID (N) == MYPE).and.(srcProcsID (LocDstCatchID (N)) == MYPE)) then ! destination is local - - I = LocDstCatchID (N) - Local_Min + 1 ! Downstream index in the local processor - K = N - Local_Min + 1 ! Source index in the local processor - - if(LocDstCatchID (N) /= N) then ! ensure not to refill the reservoir by itself - - QINFLOW(I) = QINFLOW(I) + QOUTFLOW (K) - WRIVER (I) = WRIVER (I) + QOUTFLOW (K) * real(route_dt) - - endif - - elseif ((srcProcsID (N) == MYPE).and.(srcProcsID (LocDstCatchID (N)) /= MYPE)) then - - if(srcProcsID (LocDstCatchID (N)) >= 0) then - - ! Send to downstream processor - - K = N - Local_Min + 1 ! Source index in the local processor - - call MPI_ISend(QOUTFLOW(K),1,MPI_real,srcProcsID (LocDstCatchID (N)),999,MPI_COMM_WORLD,req,status) - call MPI_WAIT(req,MPI_STATUS_IGNORE,status) - - endif - - elseif ((srcProcsID (N) /= MYPE).and.(srcProcsID (N) >= 0)) then - - K = srcProcsID (dstCatchID(N,srcProcsID (N)+1)) - - if (k == MYPE) then - - do i = 1,nDEs - - if(MYPE /= i-1) then - - if((srcProcsID (n) == i-1).and.(srcProcsID (dstCatchID(N, i)) == MYPE))then - call MPI_RECV(rbuff,1,MPI_real, srcProcsID (N),999,MPI_COMM_WORLD,MPI_STATUS_IGNORE,status) - K = dstCatchID(N,i) - Local_Min + 1 - QINFLOW (K) = QINFLOW (K) + rbuff - WRIVER (K) = WRIVER (K) + rbuff * real(route_dt) - - endif - endif - end do - endif - - endif - - end do - - ! initialize the cycle counter and sum (runoff_tile) - - runoff_save = 0. - ThisCycle = 1 - - - else - - runoff_save = runoff_save + runoff/real (N_CYC) - - ThisCycle = ThisCycle + 1 - - endif RUN_MODEL - - call MAPL_TimerOff ( MAPL, "-RRM" ) - -! All done -! -------- - - call MAPL_TimerOff(MAPL,"RUN") - - RETURN_(ESMF_SUCCESS) - - end subroutine RUN - -! --------------------------------------------------------------------------- - - subroutine InitializeRiverRouting(MYPE, numprocs, root_proc, & - pfaf_code, AllActive, AlldstCatchID, srcProcsID, LocDstCatchID, rc) - - implicit none - INTEGER, INTENT (IN) :: MYPE, numprocs - LOGICAL, INTENT (IN) :: root_proc - INTEGER, DIMENSION (:), INTENT (IN) :: pfaf_code - INTEGER, DIMENSION (N_CatG), INTENT (INOUT) :: srcProcsID, LocDstCatchID - INTEGER, DIMENSION (N_CatG,numprocs), INTENT (INOUT) :: Allactive, AlldstCatchID - - INTEGER, DIMENSION(:) ,ALLOCATABLE :: global_buff, scounts, rdispls, rcounts, LocalActive - INTEGER :: N_active, I,J,K,N,i1,i2,NProcs, Local_Min, Local_Max - - integer, optional, intent(OUT):: rc - - integer :: mpierr - character(len=ESMF_MAXSTR), parameter :: Iam='InitializeRiverRouting' - - ! STEP 1: Identify active catchments within the local processor. If the catchment is active in - ! more than 1 processor, choose an owner. - ! -------------------------------------------------------------------------------------------- - - allocate (LocalActive (1:N_CatG)) - LocalActive = -9999 - - Local_Min = minval (pfaf_code) - Local_Max = maxval (pfaf_code) - - do N = 1, size (pfaf_code) - LocalActive(pfaf_code(n)) = pfaf_code(n) - end do - - allocate (global_buff (N_CatG * numprocs)) - allocate (scounts(numprocs),rdispls(numprocs),rcounts(numprocs)) - - scounts = N_CatG - rcounts = N_CatG - - rdispls(1) = 0 - global_buff= 0 - - do i=2,numprocs - rdispls(i)=rdispls(i-1)+rcounts(i-1) - enddo - - call MPI_allgatherv ( & - LocalActive, scounts ,MPI_INTEGER, & - global_buff, rcounts, rdispls,MPI_INTEGER, & - MPI_COMM_WORLD, mpierr) - - do i=1,numprocs - Allactive (:,i) = global_buff((i-1)*N_CatG+1:i*N_CatG) - enddo - - if (root_proc) then - - DO N = 1, N_CatG - NPROCS = count(Allactive(N,:) >= 1) - if(NPROCS > 0)then - if (NPROCS == 1) then - srcProcsID (N) = maxloc(Allactive(N,:),dim=1) - 1 - else - i1 = MAX(N - 5,1) - i2 = MIN(N + 5, N_CatG) - N_active = 0 - do I = 1,numprocs - if(Allactive (N,I) >= 1) then - if(count (Allactive(I1:I2,I) > 0) > N_active) then - N_active = count (Allactive(I1:I2,I) > 0) - J = I - endif - endif - end do - srcProcsID (N) = J - 1 - endif - endif - END DO - - endif - - call MPI_BCAST (srcProcsID, N_CatG, MPI_INTEGER, 0,MPI_COMM_WORLD,mpierr) - - ! STEP 2: reset downstream catchment indeces (from -1 OR 1:291284) of catchments that are - ! in the local processor to full domain indeces. - ! ------------------------------------------------------------------------------------------ - - do N = Local_Min, Local_Max - - if(LocalActive (N) >=1) then - - if (LocDstCatchID (N) == -1) then - ! (a) DNST Catch is a sink hole, ocean or lake so water drains to self - LocDstCatchID (N) = N - - endif - - else - - LocDstCatchID (N) = -9999 ! is inactive - - endif - end do - - global_buff= 0 - - call MPI_allgatherv ( & - LocDstCatchID, scounts ,MPI_INTEGER, & - global_buff, rcounts, rdispls,MPI_INTEGER, & - MPI_COMM_WORLD, mpierr) - - do i=1,numprocs - AlldstCatchID (:,i) = global_buff((i-1)*N_CatG+1:i*N_CatG) - enddo - - deallocate (global_buff, scounts, rdispls, rcounts, LocalActive) - - RETURN_(ESMF_SUCCESS) - end subroutine InitializeRiverRouting -end module GEOS_RouteGridCompMod diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/build_rivernetwork.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/build_rivernetwork.py deleted file mode 100644 index a22817312..000000000 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/build_rivernetwork.py +++ /dev/null @@ -1,403 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Sun May 12 16:42:19 2019 - -@author: smahanam -""" -############### -# Load modules -############### - -from netCDF4 import Dataset, stringtochar -import numpy as np -import csv -import os -import datetime -import matplotlib.pyplot as plt -from matplotlib.backends.backend_pdf import PdfPages -from mpl_toolkits.basemap import Basemap -from statistics import mode -from scipy.io import FortranFile -import math -import sys -import gmplot - -###################################### -# Read streamflow station information -###################################### - -file = input("Enter station lat/lon information file name with full path \ - \n \t (Note : The .csv file format is as follows:\ - \n \t One line header is followed by a seperate line\ - \n \t for each station, please see the below example.)\ - \n \ - \n STA_LAT, STA_LON, STA_NAME\ - \n 38.709372, -91.407608, Missouri Hermann\ - \n 43.065813, -98.563029, Missouri Ft. Randall Dam\ - \n \ - \n ") -sta_info = [] -with open(file) as fh: - rd = csv.DictReader(fh, delimiter=',') - for row in rd: - sta_info.append(row) -N_STA = len(sta_info) -lats = np.zeros(shape=(N_STA)); lons = np.zeros(shape=(N_STA)); snames = ["" for x in range(N_STA)] -for x in range(0,N_STA) : lats[x], lons[x], snames[x] = sta_info[x].values() - -################################ -# Create basin information file -################################ - -OUTDIR = '../river_basin_infor/' - -if not os.path.exists(OUTDIR): - os.mkdir(OUTDIR) - -MAX_CAT_BASIN = 12500 -ncFidOut = Dataset(OUTDIR + '/Riverflow_Station_Information.nc4', mode='w', clobber='True',format='NETCDF4') -StaDim = ncFidOut.createDimension('n_sta' , N_STA) -UniDim = ncFidOut.createDimension('n_catb', MAX_CAT_BASIN) -StrDim = ncFidOut.createDimension('strl' , 40) - -# Check whether AGCM or GEOSldas -# ------------------------------ - -til_file = None -if os.path.isfile('LDAS.rc'): - til_file = os.popen('ls -1 ../output/*/rc_out/*_tilecoord.bin').read().strip() - trn_file = os.popen("grep BCS_PATH * | cut -d':' -f3").read().strip() + "/clsm/Grid2Catch_TransferData.nc" - - with FortranFile(til_file, 'r') as tf: - NTILES = tf.read_ints(np.int32) - tile_id = np.array(tf.read_ints(np.int32)) - pfaf_domain = np.array(tf.read_ints(np.int32)) - pfaf_domain = np.array(tf.read_ints(np.int32)) - com_lon = np.array(tf.read_reals(float)) - com_lat = np.array(tf.read_reals(float)) - min_lon = np.array(tf.read_reals(float)) - max_lon = np.array(tf.read_reals(float)) - min_lat = np.array(tf.read_reals(float)) - max_lat = np.array(tf.read_reals(float)) - i_indg = np.array(tf.read_ints(np.int32)) - j_indg = np.array(tf.read_ints(np.int32)) - - cat_index_domain = np.array(np.unique(np.sort(pfaf_domain))) - IM = i_indg.max() - i_indg.min() + 1 - JM = j_indg.max() - j_indg.min() + 1 - - # read Grid2Catch_TransferData.nc - Grid2Catch = Dataset(trn_file,mode='r') - ncats_in_grid = np.array (Grid2Catch.variables['NCats_in_GRID'][:]) - pfaf_index = np.array (Grid2Catch.variables['Pfaf_Index'][:]) - Grid2Catch.close() - sys.exit() -else: - print ('Enter output AGCM grid resolution 288x181, 576x361 .. etc : ' ) - imjm = input ('IMxJM : ') - im, jm = imjm.split('x') - IM = int(im) - JM = int(jm) - - DX = 360. / IM - if (JM % 2) == 1: - DY = 180. / (JM - 1) - elif(JM % 2) == 0 : - DY = 180. / JM - xdim = ncFidOut.createDimension('lon' , IM) - ydim = ncFidOut.createDimension('lat' , JM) - -ncFidOut.description = "Pfafstetter catchment indentifications for river basins" -ncFidOut.history = "Created " + datetime.datetime.now().strftime('%Y-%m-%d %H:%M:%S') - -latsOut = ncFidOut.createVariable('sta_lat','f4',('n_sta',)) -#latsOut.long_name, 'Station Latitude' -latsOut [:] = lats -lonsOut = ncFidOut.createVariable('sta_lon','f4',('n_sta',)) -#lonsOut.long_name, 'Station Longitude' -lonsOut [:] = lons -NameOut = ncFidOut.createVariable('sta_name','S1',('n_sta','strl')) -#NameOut.long_name, 'Station Name' -#for nrec in range(nrecs): -data = [] -data = np.empty((N_STA,),'S'+repr(40)) -for n in range(N_STA): - tstr = snames[n] - data[n] = tstr[0:len(tstr)] -datac = stringtochar(data) -NameOut[:] = datac -AreaOut = ncFidOut.createVariable('CompBasinArea','f4',('n_sta',)) -#AreaOut.long_name, 'Computed Basin Area' -#AreaOut.units, 'km^2' -NcatOut = ncFidOut.createVariable('NCatB','i4',('n_sta',)) -#NcatOut.long_name, 'Number of catchments in the basin' -GCIDOut = ncFidOut.createVariable('GlobalID','i4',('n_catb','n_sta')) -#GCIDOut.long_name, 'Catchment Index numbers of catchments within the basin' - -if os.path.isfile('LDAS.rc'): - NcellsBOut = ncFidOut.createVariable('N_cells_basin','i4',('n_sta',)) - SMAPIDOut = ncFidOut.createVariable('SMAPID','i4',('n_catb','n_sta')) -else: - maskOut = ncFidOut.createVariable('basin_mask','i4',('lat','lon','n_sta')) - lonout = ncFidOut.createVariable('lon','f4',('lon')) - varr = np.full (IM, 0.) - for i in range (IM): - varr [i] = -180. + DX*i - lonout [:] = varr - latout = ncFidOut.createVariable('lat','f4',('lat')) - varr = np.full (JM,0.) - varr [0]= -90. + DY/4. - varr [JM-1] = 90. - DY/4. - for i in range (1, JM -1): - varr [i] = DY*i + -90. - latout [:] = varr - -####################### -# Opening PfafCode file -####################### - -# (1) Pfafstetter Index raster -SRTMPfaf = Dataset('SRTM_PfafData.nc',mode='r') -glat = np.array (SRTMPfaf.variables['latitude'][:]) -glon = np.array (SRTMPfaf.variables['longitude'][:]) -catid = np.array (SRTMPfaf.variables['CatchIndex'][:]) -DXY = 360./21600. # Pfafstetter index raster resolution -SRTMPfaf.close() - -# (2) River Network information -RNetWork = Dataset ('RiverNetwork_information.nc4',mode='r') -num_catchs = np.array (RNetWork.variables['NUM_CATCHS'][:]) -downstream_lon = np.array (RNetWork.variables['DownStream_lon'][:]) -downstream_lat = np.array (RNetWork.variables['DownStream_lat'][:]) -upstream_lon = np.array (RNetWork.variables['UpStream_lon'][:]) -upstream_lat = np.array (RNetWork.variables['UpStream_lat'][:]) -catchment_index = np.array (RNetWork.variables['CatchmentIndex'][:]) -upstream_index = np.array (RNetWork.variables['UPSTRIndex'][:]) -downstream_index = np.array (RNetWork.variables['DNSTRIndex'][:]) -catch_area = np.array (RNetWork.variables['CATCH_AREA'][:]) - -RNetWork.close() - -################# -# Define classes -################# - -class RiverBasin(object): - - def find_pfid_at_station (self,sta_lat, sta_lon): - global catid, DXY, glat, glon - j_sta = np.asscalar(np.where(np.logical_and(glat >= sta_lat, glat < sta_lat + DXY))[0]) - i_sta = np.asscalar(np.where(np.logical_and(glon >= sta_lon, glon < sta_lon + DXY))[0]) - loc_index = np.asscalar(catid [j_sta, i_sta]) - if loc_index < 1: - loc_index = np.max(np.array(catid [j_sta-2:j_sta+3, i_sta-2:i_sta+3])) - return loc_index - - def get_continent_id (self,loc_index): - if (1 <= loc_index <= 75368): - cid = 0 # Asia - elif (75369 <= loc_index <= 140751): - cid = 1 # Africa - elif (140752 <= loc_index <= 189105): - cid = 2 # North America - elif (189106 <= loc_index <= 229074): - cid = 3 # Europe - elif (229075 <= loc_index <= 267083): - cid = 4 # South America - elif (267084 <= loc_index <= 291284): - cid = 5 # Australia - return cid - - def update_array (self, a, n, v): - if (a.ndim == 1): - a[n] = v - elif(a.ndim > 1): - a[:,n] = v - return a - - def link_upstream_catchments (self, catch_in, DOM_MIN, upst, mouth_coord, cum_area, cat_area, catch_down): - if (upst[0,catch_in] < 0): - if (catch_in != catch_down): - cum_area[catch_down] = cum_area[catch_down] + cat_area [catch_in] - return mouth_coord,cum_area - for i in range(30): - if (upst[i,catch_in] > 0): - this_catch = upst[i,catch_in] - DOM_MIN - mouth_coord[:,this_catch] = mouth_coord[:,catch_in] - if (catch_in != catch_down): - cum_area[this_catch] = cum_area[this_catch] + cat_area [this_catch] - cum_area[catch_in] = cum_area[catch_in] + cat_area [this_catch] - cum_area[catch_down] = cum_area[catch_down] + cat_area [this_catch] - if(upst[0,this_catch] > 0): - mouth_coord,cum_area = self.link_upstream_catchments (this_catch, DOM_MIN, upst, mouth_coord, - cum_area, cat_area, catch_down) - return mouth_coord,cum_area - - def derive_basin_mask (catid_mask): - global catid, DXY, glat, glon, IM, JM, DX, DY - basin_mask = np.full ((JM, IM), 0) - dx2 = math.ceil(DX / DXY /2.) - dy2 = math.ceil(DY / DXY /2.) - for j in range (1, JM - 1): - sta_lat = DY*j + -90. - j_sta = (np.where(np.logical_and(glat >= sta_lat, glat < sta_lat + 2*DXY))[0]).min() - for i in range (1, IM): - sta_lon = DX*i -180. - i_sta = (np.where(np.logical_and(glon >= sta_lon, glon < sta_lon + 2*DXY))[0]).min() - catid_domain = catid_mask [j_sta - dy2: j_sta + dy2+1, i_sta - dx2: i_sta + dx2+1] - if any(np.reshape(catid_domain,((2*dx2+1)*(2*dy2+1)))): - basin_mask [j,i] = 1 - return basin_mask - - def __init__(self, sta_lat, sta_lon): - global num_catchs, catchment_index, catch_area, upstream_index - - pfid_at_station = self.find_pfid_at_station (sta_lat, sta_lon) - self.continent_id = self.get_continent_id (pfid_at_station) - - ncat_continent = num_catchs[self.continent_id] - domain_index = np.array (catchment_index[self.continent_id,0:ncat_continent]) - catch_area_continent = np.array (catch_area [self.continent_id,0:ncat_continent]) - upst = np.array (upstream_index [self.continent_id,:,0:ncat_continent]) - - mouth_coord = np.full ((2,ncat_continent),-9999.) - cum_area = np.full (ncat_continent,0.) - this_catch = pfid_at_station - domain_index.min() - - DOM_MIN = domain_index.min() - mouth_coord = self.update_array (mouth_coord, this_catch, [sta_lat, sta_lon]) - cum_area = self.update_array (cum_area, this_catch, catch_area_continent [this_catch]) - - mouth_coord,cum_area = self.link_upstream_catchments (this_catch, DOM_MIN, upst, mouth_coord, \ - cum_area, catch_area_continent, this_catch) - - self.upst_catchs = domain_index[(np.where(mouth_coord[0,:] != -9999.))] - self.comp_basin_area = cum_area[this_catch] - -####################### -# Loop through stations -####################### - -basin_areas = [] -catch_in_basin = np.full ((MAX_CAT_BASIN,N_STA),-9999.) -ncatch_in_basin= [] -continent_id = [] -basin_masks = np.full ((JM, IM, N_STA), 0) - -for n in range(N_STA): - this_basin = RiverBasin(lats[n], lons[n]) - basin_areas.append(this_basin.comp_basin_area) - ncatch_in_basin.append(np.size(this_basin.upst_catchs)) - continent_id.append(this_basin.continent_id) - catch_in_basin [0:ncatch_in_basin[n], n] = this_basin.upst_catchs - if os.path.isfile('LDAS.rc'): - - catid_mask = np.in1d(catid,[catch_in_basin[0:ncatch_in_basin[n], n]]).reshape(catid.shape) - else: - catid_mask = np.in1d(catid,[catch_in_basin[0:ncatch_in_basin[n], n]]).reshape(catid.shape) - basin_masks [:,:,n] = RiverBasin.derive_basin_mask (catid_mask) - del this_basin - -AreaOut [:] = basin_areas -NcatOut [:] = ncatch_in_basin -GCIDOut [:] = catch_in_basin -maskOut [:] = basin_masks - -ncFidOut.close() - -############################################## -# PLOTTING -############################################## - -class BasinMaps(object): - def draw_basinmaps_html (name, lats_down,lats_up,lons_down,lons_up): - global OUTDIR - gmap = gmplot.GoogleMapPlotter(np.array([lats_down,lats_up]).mean(),np.array([lons_down,lons_up]).mean(), 5) - for i in range(lats_down.size): - gmap.scatter([lats_down[i], lats_up[i]], [lons_down[i], lons_up [i]], '# FF0000', size = 40, marker = False) - gmap.plot([lats_down[i], lats_up[i]], [lons_down[i]],'cornflowerblue', edge_width = 2.5) - gmap.draw( OUTDIR + name + ".html" ) - - def plot_basinmaps (sta_name, boundary,lats_down,lats_up,lons_down,lons_up, mask = None): - global IM, JM, DX, DY - m = Basemap(projection = 'cyl', llcrnrlat= boundary[0] - 0.25,urcrnrlat=boundary[2]+ 0.25, llcrnrlon=boundary[1]- 0.25,urcrnrlon=boundary[3] + 0.25, resolution ='c') - m.drawcoastlines() - m.fillcontinents(color='white') - m.drawmapboundary(fill_color='white') - m.drawstates(color='black') - m.drawcountries(color='black') - for i in range(lats_down.size): - x, y = m([lons_down[i],lons_up[i]], [lats_down[i], lats_up[i]]) - m.plot(x, y, linewidth=1, color = 'cornflowerblue') - if mask.any(): - for j in range (1, JM): - sta_lat = DY*j + -90. - for i in range (1, IM): - sta_lon = DX*i -180. - if mask [j,i] == 1: - xx = sta_lon + np.array([-1,-1, 1, 1,-1]) * DX/2. - yy = sta_lat + np.array([-1, 1, 1,-1,-1]) * DY/2. -# m.plot(xx, yy, linewidth=1, color = 'r') - m.plot(xx.mean(), yy.mean(), 'ro', markersize=1) - plt.title (sta_name) - -# 1) html table with google maps -################################ - -outfile = open(OUTDIR + "index.html", "w") -outfile.write(""""" - - Maps of River Basins .. - - -""") - -outfile.write ('') - -with PdfPages(OUTDIR + 'basin_maps.pdf') as pdf: - for n in range(N_STA): - sname_strip = snames[n].replace(' ','') - domain_index = np.array (catchment_index[continent_id[n],0:num_catchs[continent_id[n]]]) - DOM_MIN = domain_index.min() - upst_catid = np.array(catch_in_basin [0:ncatch_in_basin[n], n] - DOM_MIN).astype(int) - lats_up = np.array(upstream_lat[continent_id[n],upst_catid]) - lats_down = np.array(downstream_lat[continent_id[n],upst_catid]) - lons_up = np.array(upstream_lon[continent_id[n],upst_catid]) - lons_down = np.array(downstream_lon[continent_id[n],upst_catid]) - boundary = [] - boundary.append (np.array([lats_up,lats_down]).min()) - boundary.append (np.array([lons_up,lons_down]).min()) - boundary.append (np.array([lats_up,lats_down]).max()) - boundary.append (np.array([lons_up,lons_down]).max()) - - # 1) draw google maps - BasinMaps.draw_basinmaps_html (sname_strip,lats_down,lats_up,lons_down,lons_up) - outfile.write ("" % (snames[n], '' + sname_strip + '')) - - # 2) plot maps - BasinMaps.plot_basinmaps(snames[n], boundary,lats_down,lats_up,lons_down,lons_up, mask = basin_masks[:,:,n]) - pdf.savefig() - plt.close() - -outfile.write( """
NameMap
%s%s
-""") -outfile.close() - - - - - - - - - - - - - - - - - - diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/routing_model.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/routing_model.F90 deleted file mode 100644 index 4d0f6a2da..000000000 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/routing_model.F90 +++ /dev/null @@ -1,139 +0,0 @@ -MODULE routing_model - - IMPLICIT NONE - - private - public :: river_routing, SEARCH_DNST, ROUTE_DT - integer , parameter :: ROUTE_DT = 3600 - - CONTAINS - - - ! ------------------------------------------------------------------------ - - SUBROUTINE RIVER_ROUTING ( & - NCAT, & - RUNCATCH,AREACAT,LENGSC, & - WSTREAM,WRIVER, & - QSFLOW,QOUTFLOW) - - IMPLICIT NONE - INTEGER, INTENT(IN) :: NCAT - REAL, INTENT(IN), DIMENSION (NCAT) :: RUNCATCH,AREACAT,LENGSC - REAL, INTENT(INOUT),DIMENSION (NCAT) :: WSTREAM, WRIVER - REAL, INTENT(OUT), DIMENSION (NCAT) :: QSFLOW,QOUTFLOW - - REAL, PARAMETER :: K_SIMPLE = 0.111902, K_RES_MAX = 0.8 ! m1_r2com_c1 - REAL, PARAMETER :: P1 = 0.010611, P2 = 0.188556, P3 = 0.096864, & - P4 = 0.691310, P5 = 0.365747, P6 = 0.009831 ! m5_calib_240 - - INTEGER :: N,I,J - REAL :: COEFF, LS, COEFF1, COEFF2,ROFF - - ! Routing Model Input Parameters - ! ------------------------------ - !**** NCAT = NUMBER OF CATCHMENTS IN THE STUDY DOMAIN - !**** RUNCATCH = RUNOFF PRODUCED BY LAND SURFACE MODEL IN THE CATCHMENT [m3/s] - !**** AREACAT = AREA OF CATCHMENT [km^2] - !**** LENGSC = LENGTHSCALE OF CATCHMENT FOR RIVER CALCULATION [km] - ! Note: We assume LENGSC for stream to river calculation as AREACAT/LENGSC - - ! Routing Model Prognostics - ! ------------------------- - !**** WSTREAM = AMOUNT OF WATER IN "LOCAL STREAM" [m^3] - !**** WRIVER = AMOUNT OF WATER IN RIVER [m^3] - - ! Routing Model Diagnostics - ! ------------------------- - !**** QSFLOW = TRANSFER OF MOISTURE FROM STREAM VARIABLE TO RIVER VARIABLE [m^3/s] - !**** QINFLOW = TRANSFER OF RIVER WATER FROM UPSTREAM CATCHMENTS [m^3/s] - i.e. sum of - ! QOUTFLOWs from all upstream catchments. This is computed outside this subroutine - !**** QOUTFLOW = TRANSFER OF RIVER WATER TO THE DOWNSTREAM CATCHMENT [m^3/s] - - QSFLOW = 0. - QOUTFLOW = 0. - - DO N=1,NCAT - - ! Updating WSTREAM - - WSTREAM(N) = WSTREAM(N) + RUNCATCH(N) * REAL (ROUTE_DT) - LS = AREACAT(N) / (AMAX1(1.,LENGSC (N))) - ROFF = RUNCATCH(N) * AREACAT(N) - IF(ROFF < 2. ) THEN - COEFF = RESCONST (LS, P1, P2) - ELSEIF(ROFF > 10.) THEN - COEFF = RESCONST (LS, P3, P4) - ELSE - COEFF1 = RESCONST (LS, P1, P2) - COEFF2 = RESCONST (LS, P3, P4) - COEFF = COEFF1 + (ROFF - 2.)*(COEFF2 - COEFF1)/8. - ENDIF - - IF(COEFF > K_RES_MAX) COEFF = K_SIMPLE - - QSFLOW(N) = COEFF * WSTREAM(N) - WSTREAM(N) = WSTREAM(N) - QSFLOW(N) - WRIVER(N) = WRIVER(N) + QSFLOW(N) - QSFLOW(N) = QSFLOW(N) / REAL (ROUTE_DT) - - ! Updating WRIVER - - LS = AMAX1(1.,LENGSC (N)) - COEFF = RESCONST (LS, P5, P6) - IF(COEFF > K_RES_MAX) COEFF = K_SIMPLE - - QOUTFLOW(N) = COEFF * WRIVER(N) - WRIVER(N) = WRIVER(N) - QOUTFLOW(N) - QOUTFLOW(N) = QOUTFLOW(N) / REAL (ROUTE_DT) - - ENDDO - - RETURN - - END SUBROUTINE RIVER_ROUTING - -! ------------------------------------------------------------------------------------------------------- - - REAL FUNCTION RESCONST (LS, P1, P2) - - IMPLICIT NONE - - REAL, INTENT (IN) :: LS, P1, P2 - - RESCONST = P1 * ((1./LS)**P2) - - END FUNCTION RESCONST - -! ------------------------------------------------------------------------------------------------------- - - RECURSIVE SUBROUTINE SEARCH_DNST (K, NCAT_G, DNST, Pfaf_all, DNST_OUT) - - implicit none - - integer, intent (in) :: NCAT_G, K - integer, intent (in), dimension (NCAT_G) :: Pfaf_all, DNST - integer, intent (inout) :: DNST_OUT - - if (DNST(K) == -1) then - DNST_OUT = -1 - else - - if(Pfaf_all(DNST(K)) >= 1) then - DNST_OUT = Pfaf_all(DNST(K)) - else - if(DNST(DNST(K)) == -1) then - DNST_OUT = -1 - else - call SEARCH_DNST (DNST(DNST(K)), NCAT_G, DNST, Pfaf_all, DNST_OUT) - endif - endif - endif - - RETURN - - END SUBROUTINE SEARCH_DNST - -! ------------------------------------------------------------------------------------------------------- - -END MODULE routing_model diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSroute_GridComp/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSroute_GridComp/CMakeLists.txt new file mode 100644 index 000000000..c2a067edf --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSroute_GridComp/CMakeLists.txt @@ -0,0 +1,9 @@ +esma_set_this () + +set (srcs + GEOS_RouteGridComp.F90 + routing_model.F90 + reservoir.F90 + ) + +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/GEOS_RouteGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSroute_GridComp/GEOS_RouteGridComp.F90 new file mode 100644 index 000000000..8075b6b76 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSroute_GridComp/GEOS_RouteGridComp.F90 @@ -0,0 +1,1091 @@ +#include "MAPL_Generic.h" + +module GEOS_RouteGridCompMod + +!BOP +! !MODULE: GEOS_Route -- child to the "Land" gridded component. + +! !DESCRIPTION: +! {\tt GEOS\_Route} is a gridded component to route total runoff produced in +! {\tt GEOS\_Catch} (RUNOFF in {\tt GEOScatch\_GridComp} or {\tt GEOScatchCN\_GridComp}) through a 290,188 +! watersheds globally (excluding submerged catchments (watersheds) from the global list of 291,284. +! All of its calculations are done on Pfafstetter watershed space. {\tt GEOS\_Route} has no children. \\ +! +! IMPORTS : RUNOFF \\ + +! !USES: + + use ESMF + use MAPL_Mod + use MAPL_ConstantsMod + use ROUTING_MODEL, ONLY: river_routing_hyd + use reservoirMod, ONLY: RES_STATE, Reservoir + use catch_constants, ONLY: N_pfaf_g => CATCH_N_PFAFS + + implicit none + + private + + type T_RROUTE_STATE !routing related variables + private + type (ESMF_RouteHandle) :: routeHandle + type (ESMF_Field) :: field_src + type (ESMF_Field) :: field + type (RES_STATE) :: reservoir + + integer :: n_pfaf_local + integer :: nt_global + integer :: nt_local + integer :: comm + integer :: nDes + integer :: myPe + integer :: minCatch + integer :: maxCatch + integer :: route_dt + + real, allocatable :: areacat(:) ! m2 + real, allocatable :: lengsc(:) ! m + integer, allocatable :: downid(:) + + real, allocatable :: runoff_acc(:) + + real, allocatable :: lstr(:) ! m + real, allocatable :: qri_clmt(:) ! m3/s + real, allocatable :: qin_clmt(:) ! m3/s + real, allocatable :: qstr_clmt(:) ! m3/s + real, allocatable :: K(:) + real, allocatable :: Kstr(:) + + integer, allocatable :: re_order(:), to_down(:) + integer, allocatable :: send_count(:), displ_send(:) + integer, allocatable :: recv_count(:), displ_recv(:) + integer :: total_send, total_recv + end type T_RROUTE_STATE + + ! Wrapper for extracting internal state + ! ------------------------------------- + type RROUTE_WRAP + type (T_RROUTE_STATE), pointer :: PTR => null() + end type RROUTE_WRAP + + include "mpif.h" + + +! !PUBLIC MEMBER FUNCTIONS: + + public SetServices + +!EOP + +contains + +!BOP + +! !IROUTINE: SetServices -- Sets ESMF services for this component + +! !INTERFACE: + + subroutine SetServices ( GC, RC ) + +! !ARGUMENTS: + + type(ESMF_GridComp), intent(INOUT) :: GC ! gridded component + integer, optional :: RC ! return code + +! !DESCRIPTION: This version uses the MAPL\_GenericSetServices. This function sets +! the Initialize and Finalize services, as well as allocating +! our instance of a generic state and putting it in the +! gridded component (GC). Here we only need to set the run method and +! add the state variable specifications (also generic) to our instance +! of the generic state. This is the way our true state variables get into +! the ESMF\_State INTERNAL, which is in the MAPL\_MetaComp. + +!EOP + +!============================================================================= +! +! ErrLog Variables + + + character(len=ESMF_MAXSTR) :: IAm + integer :: STATUS + character(len=ESMF_MAXSTR) :: COMP_NAME + +! Local derived type aliases + + type (ESMF_Config ) :: CF + + type (T_RROUTE_STATE), pointer :: route_internal_state => null() + type (RROUTE_wrap) :: wrap + +!============================================================================= + +! Begin... + +!------------------------------------------------------------ +! Get my name and set-up traceback handle +!------------------------------------------------------------ + + call ESMF_GridCompGet(GC ,& + NAME=COMP_NAME ,& + RC=STATUS ) + + VERIFY_(STATUS) + + Iam = trim(COMP_NAME) // 'SetServices' + +! ----------------------------------------------------------- +! Set the Initialize and Run entry points +! ----------------------------------------------------------- + + call MAPL_GridCompSetEntryPoint ( GC, ESMF_METHOD_INITIALIZE, Initialize, RC=STATUS) + VERIFY_(STATUS) + call MAPL_GridCompSetEntryPoint ( GC, ESMF_METHOD_RUN, Run, RC=STATUS) + VERIFY_(STATUS) + call MAPL_GridCompSetEntryPoint ( GC, ESMF_METHOD_FINALIZE, Finalize, RC=status) + VERIFY_(status) + +! ----------------------------------------------------------- +! Get the configuration +! ----------------------------------------------------------- +! + call ESMF_GridCompGet( GC, CONFIG = CF, RC=STATUS ) + VERIFY_(STATUS) +! +! ----------------------------------------------------------- +! Set the state variable specs. +! ----------------------------------------------------------- + +!BOS + +! ----------------------------------------------------------- +! Import States +! ----------------------------------------------------------- +! WY note: Here TileOnly is on tile space + + call MAPL_AddImportSpec(GC, & + LONG_NAME = 'runoff_total_flux' ,& + UNITS = 'kg m-2 s-1' ,& + SHORT_NAME = 'RUNOFF' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + _RC ) + +!!!!!!!!!!!!!!!! +! Internal +!!!!!!!!!!!!!! +! WY note: Here TileOnly is on *Pfafstetter* catchment space + + call MAPL_AddInternalSpec(GC ,& + LONG_NAME = 'volume_of_water_in_local_stream',& + UNITS = 'm+3' ,& + SHORT_NAME = 'WSTREAM' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RESTART = MAPL_RestartRequired ,& + _RC ) + + call MAPL_AddInternalSpec(GC ,& + LONG_NAME = 'volume_of_water_in_river' ,& + UNITS = 'm+3' ,& + SHORT_NAME = 'WRIVER' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RESTART = MAPL_RestartRequired ,& + _RC ) + + + call MAPL_AddInternalSpec(GC ,& + LONG_NAME = 'reservoir_storage' ,& + UNITS = 'm+3' ,& + SHORT_NAME = 'WRES' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RESTART = MAPL_RestartRequired ,& + _RC ) + +!!!!!!!!!!!!!!!! +! Export +!!!!!!!!!!!!!!! +! WY note: Here TileOnly is on *Pfafstetter* catchment space + + call MAPL_AddExportSpec(GC, & + LONG_NAME = 'transfer_of_moisture_from_stream_variable_to_river_variable' ,& + UNITS = 'm+3 s-1' ,& + SHORT_NAME = 'QSFLOW' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + _RC ) + + call MAPL_AddExportSpec(GC, & + LONG_NAME = 'transfer_of_river_water_from_upstream_catchments' ,& + UNITS = 'm+3 s-1' ,& + SHORT_NAME = 'QINFLOW' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + _RC ) + + call MAPL_AddExportSpec(GC, & + LONG_NAME = 'transfer_of_river_water_to_downstream_catchments' ,& + UNITS = 'm+3 s-1' ,& + SHORT_NAME = 'QOUTFLOW' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + _RC ) + + call MAPL_AddExportSpec(GC, & + LONG_NAME = 'reservoir_discharge' ,& + UNITS = 'm+3 s-1' ,& + SHORT_NAME = 'QRES' ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + _RC ) +!EOS + + call MAPL_TimerAdd(GC, name="-RRM" ,RC=STATUS) + VERIFY_(STATUS) + + + +! Allocate this instance of the internal state and put it in wrapper. +! ------------------------------------------------------------------- + + allocate( route_internal_state, stat=status ) + VERIFY_(STATUS) + wrap%ptr => route_internal_state + +! Save pointer to the wrapped internal state in the GC +! ---------------------------------------------------- + + call ESMF_UserCompSetInternalState ( GC, 'RiverRoute_state',wrap,status ) + VERIFY_(STATUS) + +! Clocks +!------- + + call MAPL_TimerAdd(GC, name="INITIALIZE" ,RC=STATUS) + VERIFY_(STATUS) + call MAPL_TimerAdd(GC, name="RUN" ,RC=STATUS) + VERIFY_(STATUS) + +! All done +!--------- + + call MAPL_GenericSetServices(GC, RC=STATUS) + VERIFY_(STATUS) + + RETURN_(ESMF_SUCCESS) + + end subroutine SetServices + +! ----------------------------------------------------------- +! INITIALIZE -- Initialize method for the route component +! ----------------------------------------------------------- + + subroutine INITIALIZE (GC,IMPORT, EXPORT, CLOCK, RC ) + +! ----------------------------------------------------------- +! !ARGUMENTS: +! ----------------------------------------------------------- + + type(ESMF_GridComp), intent(inout) :: GC + type(ESMF_State), intent(inout) :: IMPORT + type(ESMF_State), intent(inout) :: EXPORT + type(ESMF_Clock), intent(inout) :: CLOCK + integer, optional, intent( out) :: RC + +! ----------------------------------------------------------- +! ErrLog Variables +! ----------------------------------------------------------- + + character(len=ESMF_MAXSTR) :: IAm="Initialize" + integer :: STATUS + character(len=ESMF_MAXSTR) :: COMP_NAME + +! ----------------------------------------------------------- +! Locals +! ----------------------------------------------------------- + type (ESMF_VM) :: VM + integer :: comm + integer :: nDEs + integer :: myPE + integer :: beforeMe, minCatch, maxCatch, i + integer :: n_pfaf_local, nt_global + integer :: ROUTE_DT, route_flag + + type(ESMF_Grid) :: tileGrid + type(ESMF_Grid) :: pfaf_tilegrid, pfaf_grid + character(len=ESMF_MAXSTR) :: SURFRC + type(ESMF_Config) :: SCF + + type(MAPL_MetaComp), pointer :: MAPL + type(MAPL_LocStream) :: locstream, pfaf_LocStream + + character(len=ESMF_MAXSTR) :: RIVER_INPUT_FILE + character(len=ESMF_MAXSTR) :: TILE_PFAF_FILE + + type(ESMF_Grid) :: agrid + type(ESMF_DELayout) :: layout + type(ESMF_DistGrid) :: dist_grid + + integer, pointer :: ims(:) => NULL() + + type (T_RROUTE_STATE), pointer :: route => null() + type (RROUTE_wrap) :: wrap + real, allocatable :: tmp_real(:) + integer, allocatable :: tmp_int(:) + + type(ESMF_Time) :: CurrentTime + type(ESMF_Alarm) :: CollectWaterAlarm + type(ESMF_TimeInterval) :: CollectWater_DT, ModelTimeStep + character(len=3) :: resname + type(Netcdf4_Fileformatter) :: formatter + integer :: j,nt_local, mpierr + logical :: use_res + + ! ------------------ + ! begin + + call MAPL_GetObjectFromGC ( GC, MAPL, RC=STATUS) + VERIFY_(STATUS) + ! get LocStream + call MAPL_Get(MAPL, LocStream = locstream, RC=status) + VERIFY_(STATUS) + + call ESMF_UserCompGetInternalState ( GC, 'RiverRoute_state',wrap,status ) + VERIFY_(STATUS) + + route => wrap%ptr + ! get vm + ! extract comm + call ESMF_VMGetCurrent(VM, RC=STATUS) + VERIFY_(STATUS) + call ESMF_VMGet (VM, mpiCommunicator =comm, RC=STATUS) + VERIFY_(STATUS) + call ESMF_VMGet (VM, localpet=MYPE, petcount=nDEs, RC=STATUS) + VERIFY_(STATUS) + + + route%comm = comm + route%ndes = ndes + route%mype = mype + + allocate(ims(1:ndes)) + ! define catchment space for this processor + call MAPL_DecomposeDim ( N_pfaf_g,ims,ndes ) ! ims(mype+1) gives the size of my partition + ! myPE is 0-based! + beforeMe = sum(ims(1:mype)) + minCatch = beforeMe + 1 + maxCatch = beforeMe + ims(myPe+1) + + ! Get grid info from the gridcomp + call ESMF_GridCompGet(gc, grid=agrid, rc=status) + VERIFY_(status) + call ESMF_GridGet(agrid, distGrid=dist_grid, _RC) + call ESMF_DistGridGet(dist_grid, delayout=layout, _RC) + + call MAPL_LocStreamGet(locstream, & + tileGrid=tilegrid, nt_local=nt_local, nt_global=nt_global, & + RC=status) + VERIFY_(STATUS) + + route%nt_global = nt_global + route%nt_local = nt_local + n_pfaf_local = maxCatch-minCatch+1 + route%n_pfaf_local= n_pfaf_local + route%minCatch = minCatch + route%maxCatch = maxCatch + + call MAPL_GetResource (MAPL, SURFRC, label = 'SURFRC:', default = 'GEOS_SurfaceGridComp.rc', RC=STATUS) ; VERIFY_(STATUS) + SCF = ESMF_ConfigCreate(rc=status) ; VERIFY_(STATUS) + call ESMF_ConfigLoadFile(SCF,SURFRC,rc=status) ; VERIFY_(STATUS) + call MAPL_GetResource (SCF, route_flag, label='RUN_ROUTE:', DEFAULT=1, RC=STATUS ) + if(route_flag==2)then + use_res=.True. + else + use_res=.False. + endif + + call MAPL_GetResource (SCF, ROUTE_DT, label='RRM_DT:', DEFAULT=3600, RC=STATUS ) + route%route_dt = ROUTE_DT + + call MAPL_GetResource (MAPL, RIVER_INPUT_FILE, label = 'RIVER_INPUT_FILE:', default = '../input/river_input.nc', RC=STATUS ) + if (MAPL_AM_I_Root()) then + call formatter%open(RIVER_INPUT_FILE, PFIO_READ, _RC) + endif + allocate(tmp_real(n_pfaf_g)) + allocate(tmp_int(n_pfaf_g)) + + ! read areacat + if (MAPL_AM_I_Root()) then + call formatter%get_var('area_catch', tmp_real(:), _RC) + endif + call MAPL_CommsBcast(layout, tmp_real, n_pfaf_g, MAPL_Root, status) + allocate(route%areacat(n_pfaf_local), source = tmp_real(minCatch:maxCatch)) + route%areacat=route%areacat*1.e6 + + !read lengsc + if (MAPL_AM_I_Root()) then + call formatter%get_var('lengsc', tmp_real(:), _RC) + endif + call MAPL_CommsBcast(layout, tmp_real, n_pfaf_g, MAPL_Root, status) + allocate(route%lengsc(n_pfaf_local), source = tmp_real(minCatch:maxCatch)) + route%lengsc=route%lengsc*1.e3 + + !read downid + if (MAPL_AM_I_Root()) then + call formatter%get_var('downid', tmp_int(:), _RC) + endif + call MAPL_CommsBcast(layout, tmp_int, n_pfaf_g, MAPL_Root, status) + allocate(route%downid(n_pfaf_local), source = tmp_int(minCatch:maxCatch)) + + !read lstr + if (MAPL_AM_I_Root()) then + call formatter%get_var('lstr', tmp_real(:), _RC) + endif + call MAPL_CommsBcast(layout, tmp_real, n_pfaf_g, MAPL_Root, status) + allocate(route%lstr(n_pfaf_local), source = tmp_real(minCatch:maxCatch)) + route%lstr = route%lstr*1.e3 + + !read K + if (MAPL_AM_I_Root()) then + call formatter%get_var('K', tmp_real(:), _RC) + endif + call MAPL_CommsBcast(layout, tmp_real, n_pfaf_g, MAPL_Root, status) + allocate(route%K(n_pfaf_local), source = tmp_real(minCatch:maxCatch)) + + !read Kstr + if (MAPL_AM_I_Root()) then + call formatter%get_var('Kstr', tmp_real(:), _RC) + endif + call MAPL_CommsBcast(layout, tmp_real, n_pfaf_g, MAPL_Root, status) + allocate(route%Kstr(n_pfaf_local), source = tmp_real(minCatch:maxCatch)) + + !read qri_clmt + if (MAPL_AM_I_Root()) then + call formatter%get_var('qri_clmt', tmp_real(:), _RC) + endif + call MAPL_CommsBcast(layout, tmp_real, n_pfaf_g, MAPL_Root, status) + allocate(route%qri_clmt(n_pfaf_local), source = tmp_real(minCatch:maxCatch)) + + !read qin_clmt + if (MAPL_AM_I_Root()) then + call formatter%get_var('qin_clmt', tmp_real(:), _RC) + endif + call MAPL_CommsBcast(layout, tmp_real, n_pfaf_g, MAPL_Root, status) + allocate(route%qin_clmt(n_pfaf_local), source = tmp_real(minCatch:maxCatch)) + + !read qstr_clmt + if (MAPL_AM_I_Root()) then + call formatter%get_var('qstr_clmt', tmp_real(:), _RC) + endif + call MAPL_CommsBcast(layout, tmp_real, n_pfaf_g, MAPL_Root, status) + allocate(route%qstr_clmt(n_pfaf_local), source = tmp_real(minCatch:maxCatch)) + + deallocate(tmp_real, tmp_int) + + if (MAPL_AM_I_Root()) then + call formatter%close() + endif + + allocate(route%runoff_acc(nt_local), source = 0.) + + + !Initial reservoir module + route%reservoir = Reservoir(layout, trim(RIVER_INPUT_FILE), N_pfaf_g,n_pfaf_local, minCatch,maxCatch,use_res, _RC) + + if(mapl_am_I_root()) print *,"reservoir init success" + + pfaf_grid = create_pfaf_grid(_RC) + call MAPL_LocstreamCreate(pfaf_Locstream, pfaf_Grid, _RC) + call MAPL_LocStreamGet(pfaf_LocStream, TILEGRID=pfaf_tilegrid, _RC) + + call MAPL_GetResource (MAPL, TILE_PFAF_File, label = 'TILE_PFAF_FILE:', default = '../input/tile_pfaf.nc4', RC=STATUS ) + + call create_mapping_handler(trim(TILE_PFAF_File), tilegrid, pfaf_tilegrid, _RC) + + call MAPL%grid%set(pfaf_grid, _RC) + + call ESMF_GridCompSet(gc, grid=pfaf_grid, RC=status) + VERIFY_(STATUS) + + call MAPL_set(MAPL, locstream = pfaf_locstream, rc=status) + VERIFY_(STATUS) + + call setup_exchange_water(pfaf_tilegrid, _RC) + + call ESMF_TimeIntervalSet(CollectWater_DT, s=ROUTE_DT, rc=status) + VERIFY_(status) + call ESMF_ClockGet(clock, currTime=CurrentTime, timeStep=ModelTimeStep, rc=status) + CollectWaterAlarm = ESMF_AlarmCreate( & + clock, & + name='CollectWater', & + ringTime=CurrentTime - ModelTimeStep+CollectWater_DT, & + ringInterval=CollectWater_DT, & + ringTimeStepCount=1, & + sticky=.false., & + rc=status & + ) + VERIFY_(status) + + + deallocate(ims) + call MAPL_GenericInitialize ( GC, import, export, clock, rc=status ) + VERIFY_(STATUS) + RETURN_(ESMF_SUCCESS) + + contains + + function create_pfaf_grid(rc) result(pfaf_grid) + type (ESMF_Grid) :: pfaf_grid + integer, optional, intent(out) :: rc + integer :: status + real(kind=8), pointer :: centers(:,:) + ! create catchment grid and it is tile space + pfaf_Grid = ESMF_GridCreate( & + name='CATCHMENT_GRID', & + countsPerDEDim1=IMs, & + countsPerDEDim2=[1], & + indexFlag=ESMF_INDEX_DELOCAL, & + coordDep1 = (/1,2/), & + coordDep2 = (/1,2/), & + gridEdgeLWidth = (/0,0/), & + gridEdgeUWidth = (/0,0/), & + _RC) + ! coord and centers are required for a valid grid, + ! even if their values don't make sense; + ! later on, the coord will be set to catchment's lat lon. + call ESMF_GridAddCoord(pfaf_Grid, _RC) + _VERIFY(status) + + call ESMF_GridGetCoord(pfaf_Grid, coordDim=1, localDE=0, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + farrayPtr=centers, _RC) + centers = 0 ! ?? just assign + call ESMF_GridGetCoord(pfaf_Grid, coordDim=2, localDE=0, & + staggerloc=ESMF_STAGGERLOC_CENTER, & + farrayPtr=centers, _RC) + centers = 0 ! + + _RETURN(_SUCCESS) + end function create_pfaf_grid + + ! ---------------------------------------------------- + + subroutine create_mapping_handler(tile_pfaf_file, tilegrid, pfaf_tilegrid, rc) + + character(*), intent(in) :: tile_pfaf_file + type(ESMF_Grid), intent(in) :: tilegrid + type(ESMF_Grid), intent(in) :: pfaf_tilegrid + integer, optional, intent(out) :: rc + + integer :: status, nWeights, nLocal_weights + integer, allocatable :: global_id(:) + logical, allocatable :: mask(:) + integer, allocatable :: srcIndices(:), positions(:), factorIndexList(:,:) + real, allocatable :: weights(:), global_frac(:), global_area(:) + integer, allocatable :: local_src(:), local_dst(:), global_src(:), global_dst(:) + integer :: unit + integer, pointer :: pfaf_index(:), local_id(:) + real , pointer :: tilearea(:) + type(Netcdf4_Fileformatter) :: formatter + type(Filemetadata) :: meta + character(len=MAPL_TileNameLength), pointer :: GNAMES(:) + + ! create source for orignal tile space + route%field_src = ESMF_FieldCreate(grid=tilegrid, typekind=ESMF_TYPEKIND_R4, _RC) + + ! create destination for pfaf tile space + route%field = ESMF_FieldCreate(grid=pfaf_tilegrid, typekind=ESMF_TYPEKIND_R4,_RC) + + call MAPL_LocstreamGet(LOCSTREAM, GRIDNAMES=GNAMES, pfaf_index=pfaf_index, tilearea=tilearea, local_id=local_id, _RC) + ! ESMF use global indices increasing with mpi_rank, no mask here for tile grid + allocate(global_id(route%nt_global)) + call ESMFL_Fcollect(tilegrid, global_id, local_id, _RC) + + nweights = route%nt_global + + if (index(GNAMES(1), 'EASEv') /=0) then + + if (MAPL_AM_I_ROOT()) then + call formatter%open(tile_pfaf_file, PFIO_READ, _RC) + meta = formatter%read(rc=status) + nweights = meta%get_dimension('tile') + endif + call MAPL_CommsBcast(layout, nWeights, 1, MAPL_Root, status) + allocate(global_src(nWeights), global_dst(nWeights), global_frac(nWeights)) + if (MAPL_AM_I_ROOT()) then + call formatter%get_var('tile_id', global_src, _RC) + call formatter%get_var('pfaf_index', global_dst, _RC) + call formatter%get_var('pfaf_frac', global_frac, _RC) + endif + call MAPL_CommsBcast(layout, global_src, nWeights, MAPL_Root, status) + call MAPL_CommsBcast(layout, global_dst, nWeights, MAPL_Root, status) + call MAPL_CommsBcast(layout, global_frac,nWeights, MAPL_Root, status) + else + allocate(global_src(nWeights), global_dst(nWeights), global_frac(nWeights)) + allocate(global_area(nWeights)) + global_src = global_id + call ESMFL_Fcollect(tilegrid, global_dst, pfaf_index, _RC) + call ESMFL_Fcollect(tilegrid, global_area, tilearea, _RC) + global_frac = global_area/route%areacat(global_dst) + deallocate(global_area) + endif + + allocate(mask(nWeights)) + mask = minCatch <= global_dst(:) .and. global_dst(:) <=maxCatch + local_src = pack(global_src(:), mask) + local_dst = pack(global_dst(:), mask) + weights = pack(global_frac(:), mask) + deallocate(global_src, global_dst, global_frac) + + ! mapping form local to global index + nLocal_weights = count(mask) + allocate(srcIndices(nLocal_weights)) + do i =1, nLocal_weights + positions = pack([(j, j=1, nt_global)], global_id == local_src(i)) + srcIndices(i) = positions(1) + enddo + + allocate(factorIndexList(2, nlocal_weights)) + factorIndexList(1,:) = srcIndices + factorIndexList(2,:) = local_dst + call ESMF_FieldSMMStore(route%field_src, route%field, & + routeHandle=route%routeHandle, & + factorList=weights, & + factorIndexList= factorIndexList, & + _RC) + + ! testing + !call ESMF_FieldGet(route%field_src, farrayPtr=ptr, rc=status) + !VERIFY_(STATUS) + !ptr(:) = 1. + + !call ESMF_FieldSMM(route%field_src, route%Field, & + ! route%routeHandle, rc=status) + !VERIFY_(STATUS) + + !call ESMF_FieldGet(route%Field, farrayPtr=ptr, rc=status) + !VERIFY_(STATUS) + !! after remapping, all values should be 1 + !if (route%mype == 20) then + ! do i = 1, n_pfaf_local + ! print* , 'ptr(i): ',i, ptr(i) + ! enddo + !endif + _RETURN(_SUCCESS) + + end subroutine create_mapping_handler + + ! ------------------------------------------------------------------------------------------- + + subroutine setup_exchange_water(pfaf_tilegrid, rc) + type(ESMF_Grid), intent(in) :: pfaf_tilegrid + integer, optional, intent(out) :: rc + integer :: pf, down_id, rank + integer, allocatable :: cat_to_ranks_global(:) + integer, allocatable :: cat_to_ranks_local(:) + integer :: status, mype,mpierr, k + integer, allocatable :: to_downstream(:), positions(:) + + mype = route%mype + allocate(cat_to_ranks_local(route%n_pfaf_local)) + allocate(cat_to_ranks_global(n_pfaf_g)) + cat_to_ranks_local = mype + + call ESMFL_FCollect(pfaf_tilegrid, cat_to_ranks_global, cat_to_ranks_local, _RC) + + allocate(route%send_count(route%ndes), source=0) + allocate(route%recv_count(route%ndes), source=0) + + do pf = 1, route%n_pfaf_local + down_id = route%downid(pf) + if (down_id == -1) cycle + ! down stream is not in the process + if (down_id < minCatch .or. maxCatch < down_id) then + rank = cat_to_ranks_global(down_id) + route%send_count(rank+1) = route%send_count(rank+1) + 1 + endif + enddo + + call MPI_AlltoALL(route%send_count, 1, MPI_INTEGER, & + route%recv_count, 1, MPI_INTEGER, route%comm, status) + VERIFY_(STATUS) + + allocate(route%displ_send(ndes), source = 0) + allocate(route%displ_recv(ndes), source = 0) + do rank = 1, ndes-1 + route%displ_send(rank+1) = route%displ_send(rank) + route%send_count(rank) + route%displ_recv(rank+1) = route%displ_recv(rank) + route%recv_count(rank) + enddo + k = sum(route%send_count) + route%total_send = k + allocate(to_downstream(k)) + allocate(route%re_order(k)) + allocate(positions(ndes), source = route%displ_send ) + positions = positions + 1 + k = 0 + do pf = 1, route%n_pfaf_local + down_id = route%downid(pf) + if (down_id == -1) cycle + if (down_id < minCatch .or. maxCatch < down_id) then + rank = cat_to_ranks_global(down_id) + k = k + 1 + to_downstream(k) = down_id + route%re_order(k)= positions(rank+1) + positions(rank+1)= positions(rank+1) + 1 + endif + enddo + if (route%total_send > 0) then + to_downstream = to_downstream(route%re_order) + endif + k = sum(route%recv_count) + route%total_recv = k + allocate(route%to_down(k)) + call MPI_AllToAllv(to_downstream, route%send_count, route%displ_send, MPI_Integer, & + route%to_down, route%recv_count, route%displ_recv, MPI_INTEGER, & + route%comm, status) + + ! testing + !do i = 1, route%total_recv + ! down_id = route%to_down(i) + ! if (down_id < minCatch .or. maxCatch < down_id) then + ! _ASSERT(.false., "Got the down_id that does not belong to me") + ! endif + !enddo + + _VERIFY(STATUS) + _RETURN(_SUCCESS) + end subroutine setup_exchange_water + + end subroutine INITIALIZE + + ! -------------------------------------------------------------------------------- + + subroutine RUN (GC,IMPORT, EXPORT, CLOCK, RC ) + +! ----------------------------------------------------------- +! !ARGUMENTS: +! ----------------------------------------------------------- + + type(ESMF_GridComp), intent(inout) :: GC + type(ESMF_State), intent(inout) :: IMPORT + type(ESMF_State), intent(inout) :: EXPORT + type(ESMF_Clock), intent(inout) :: CLOCK + integer, optional, intent( out) :: RC + +! ----------------------------------------------------------- +! ErrLog Variables +! ----------------------------------------------------------- + + character(len=ESMF_MAXSTR) :: IAm='RUN' + integer :: STATUS + character(len=ESMF_MAXSTR) :: COMP_NAME + +! ----------------------------------------------------------- +! Locals +! ----------------------------------------------------------- + type (ESMF_State ) :: INTERNAL + type (MAPL_MetaComp), pointer :: MAPL +! type(ESMF_Alarm) :: ALARM + type (ESMF_Config ) :: CF + +! ----------------------------------------------------- +! IMPORT pointers +! ----------------------------------------------------- + + real, dimension(:), pointer :: RUNOFF_SRC0 + +! ----------------------------------------------------- +! INTERNAL pointers +! ----------------------------------------------------- + + real, dimension(:), pointer :: AREACAT + real, dimension(:), pointer :: LENGSC + real, dimension(:), pointer :: DNSTR + real, dimension(:), pointer :: WSTREAM + real, dimension(:), pointer :: WRIVER + real, dimension(:), pointer :: WRES + real, dimension(:), pointer :: LRIVERMOUTH + real, dimension(:), pointer :: ORIVERMOUTH + +! ----------------------------------------------------- +! EXPORT pointers +! ----------------------------------------------------- + + real, dimension(:), pointer :: QSFLOW + real, dimension(:), pointer :: QINFLOW + real, dimension(:), pointer :: QOUTFLOW + real, dimension(:), pointer :: QRES + +! Time attributes and placeholders + +! type(ESMF_Time) :: CURRENT_TIME + +! Others + + type(ESMF_Grid) :: TILEGRID + type (MAPL_LocStream) :: LOCSTREAM + + integer :: n_pfaf_local, N_CYC + integer :: ROUTE_DT + + integer :: I + REAL :: HEARTBEAT + REAL, ALLOCATABLE, DIMENSION(:) :: RUNOFF_IN,AREACAT_ACT,& + LENGSC_ACT, QSFLOW_OUT,QOUTFLOW_OUT,QRES_OUT,QOUT_CAT + type(ESMF_Field) :: runoff_src + + integer :: ndes, mype + type (T_RROUTE_STATE), pointer :: route + type (RROUTE_wrap) :: wrap + + integer :: nt_global,nt_local + integer,save :: nstep_per_day + + type(ESMF_Time) :: CurrentTime, nextTime + integer :: YY,MM,DD,HH,MMM,SS,YY_next,MM_next,DD_next + character(len=4) :: yr_s + character(len=2) :: mon_s,day_s + + real, pointer :: arrayPtr(:) + type (RES_STATE), pointer :: res + + real, allocatable :: WTOT_BEFORE(:),QINFLOW_LOCAL(:) + + type(ESMF_Alarm) :: CollectWaterAlarm + + ! ------------------ + ! begin + call ESMF_UserCompGetInternalState ( GC, 'RiverRoute_state',wrap,status ) + VERIFY_(STATUS) + route => wrap%ptr + +! Get the target components name and set-up traceback handle. +! ----------------------------------------------------------- + + call ESMF_GridCompGet(GC, name=COMP_NAME, CONFIG=CF, RC=STATUS ) + VERIFY_(STATUS) + Iam = trim(COMP_NAME) // "::RUN" + +! Get my internal MAPL_Generic state +! ----------------------------------------------------------- + + call MAPL_GetObjectFromGC(GC, MAPL, STATUS) + VERIFY_(STATUS) + call MAPL_Get(MAPL, HEARTBEAT = HEARTBEAT, RC=STATUS) + VERIFY_(STATUS) + + call ESMF_ClockGetAlarm(clock, 'CollectWater', CollectWaterAlarm, _RC) + !if (mapl_am_I_root()) print *, "HEARTBEAT=",HEARTBEAT + +! internal + call MAPL_Get(MAPL, INTERNAL_ESMF_STATE=INTERNAL, _RC) + call MAPL_GetPointer(INTERNAL, WRIVER, 'WRIVER', _RC ) + call MAPL_GetPointer(INTERNAL, WSTREAM,'WSTREAM', _RC) + call MAPL_GetPointer(INTERNAL, WRES, 'WRES', _RC) + +! export + call MAPL_GetPointer(EXPORT, QSFLOW, 'QSFLOW', _RC) + call MAPL_GetPointer(EXPORT, QINFLOW, 'QINFLOW' , _RC) + call MAPL_GetPointer(EXPORT, QOUTFLOW, 'QOUTFLOW', _RC) + call MAPL_GetPointer(EXPORT, QRES, 'QRES', _RC) + +! Start timers +! ------------ + call MAPL_TimerOn(MAPL,"RUN") +! Get parameters from generic state +! --------------------------------- + + ! call MAPL_Get(MAPL, INTERNAL_ESMF_STATE=INTERNAL, RC=STATUS) + ! VERIFY_(STATUS) + +! get pointers to inputs variables +! ---------------------------------- + + ndes = route%ndes + mype = route%mype + n_pfaf_local = route%n_pfaf_local + nt_global = route%nt_global + nt_local = route%nt_local + res => route%reservoir + ROUTE_DT = route%route_dt + + ! get the field from IMPORT + call ESMF_StateGet(IMPORT, 'RUNOFF', field=runoff_src, RC=STATUS) + VERIFY_(STATUS) + call ESMF_FieldGet(runoff_src, farrayPtr=RUNOFF_SRC0, rc=status) + VERIFY_(STATUS) + call MAPL_Get(MAPL, LocStream=LOCSTREAM, RC=STATUS) + VERIFY_(STATUS) + call MAPL_LocStreamGet(LOCSTREAM, TILEGRID=TILEGRID, RC=STATUS) + VERIFY_(STATUS) + call MAPL_TimerOn ( MAPL, "-RRM" ) + + ! For efficiency, the time step to call the river routing model is set at ROUTE_DT + N_CYC = ROUTE_DT/HEARTBEAT + if (ESMF_AlarmIsRinging(CollectWaterAlarm)) then + + !accumulates runoff + route%runoff_acc = (route%runoff_acc + RUNOFF_SRC0)/real (N_CYC) + + !Gets time used for output and restart + call ESMF_ClockGet(clock, currTime=CurrentTime, rc=status) + call ESMF_TimeGet(CurrentTime, yy=YY, mm=MM, dd=DD, h=HH, m=MMM, s=SS, rc=status) + call ESMF_ClockGetNextTime(clock, nextTime=nextTime, rc=status) + call ESMF_TimeGet(nextTime, yy=YY_next, mm=MM_next, dd=DD_next, rc=status) + write(yr_s, '(I4.4)')YY + write(mon_s,'(I2.2)')MM + write(day_s,'(I2.2)')DD + + ! redistribute runoff from tile space to catchment space + call ESMF_FieldGet(route%field_src, farrayPtr=arrayPtr, rc=status) + VERIFY_(STATUS) + ArrayPtr = route%runoff_acc(:) + call ESMF_FieldSMM(srcField=route%field_src, dstField=route%Field, & + routeHandle=route%routeHandle, rc=rc) + call ESMF_FieldGet(route%field, farrayPtr=arrayPtr, rc=status) + VERIFY_(STATUS) + RUNOFF_IN = arrayPtr * route%areacat/1000. + + + ! Prepares to conduct routing model + allocate (AREACAT_ACT (n_pfaf_local)) + allocate (LENGSC_ACT (n_pfaf_local)) + allocate (QSFLOW_OUT (n_pfaf_local)) + allocate (QOUTFLOW_OUT(n_pfaf_local),QRES_OUT(n_pfaf_local),QOUT_CAT(n_pfaf_local)) + + QRES_OUT=0. + LENGSC_ACT =route%lengsc/1.e3 !m->km + AREACAT_ACT=route%areacat/1.e6 !m2->km2 + + allocate(WTOT_BEFORE(n_pfaf_local)) + WTOT_BEFORE=WSTREAM + WRIVER + WRES + + ! Call river_routing_model + ! ------------------------ + CALL RIVER_ROUTING_HYD (n_pfaf_local, ROUTE_DT,& + RUNOFF_IN, route%lengsc, route%lstr, & + route%qstr_clmt, route%qri_clmt, route%qin_clmt, & + route%K, route%Kstr, & + WSTREAM,WRIVER, & + QSFLOW_OUT,QOUTFLOW_OUT) + + ! Call reservoir module + call res%calc( QOUTFLOW_OUT, QRES_OUT, WRES, real(route_dt), _RC) + QOUT_CAT = QOUTFLOW_OUT + where(res%active_res==1) QOUT_CAT=QRES_OUT + allocate(QINFLOW_LOCAL(n_pfaf_local)) + call exchange_water(QOUT_CAT, QINFLOW_LOCAL, _RC) + WRIVER = WRIVER + QINFLOW_LOCAL*real(route_dt) + + ! Check balance if needed + !call check_balance(route,n_pfaf_local,nt_local,runoff_acc,WRIVER_ACT,WSTREAM_ACT,WTOT_BEFORE,RUNOFF_ACT,QINFLOW_LOCAL,QOUT_CAT,FirstTime,yr_s,mon_s) + + if (associated(QSFLOW)) QSFLOW = QSFLOW_OUT + if (associated(QOUTFLOW)) QOUTFLOW = QOUTFLOW_OUT + if (associated(QRES)) QRES = QRES_OUT + deallocate(RUNOFF_IN,AREACAT_ACT,LENGSC_ACT,QOUTFLOW_OUT,QINFLOW_LOCAL,QSFLOW_OUT,WTOT_BEFORE,QRES_OUT,QOUT_CAT) + + route%runoff_acc = 0. + else + + route%runoff_acc = route%runoff_acc + RUNOFF_SRC0 + + endif ! alarm + + ! All done + ! -------- + call MAPL_TimerOff ( MAPL, "-RRM" ) + call MAPL_TimerOff(MAPL,"RUN") + !call MPI_Barrier(route%comm, mpierr) + + RETURN_(ESMF_SUCCESS) + + contains + + subroutine exchange_water(cat_out, cat_in, rc) + real, intent(in) :: cat_out(:) + real, intent(out) :: cat_in(:) + integer, optional, intent(out) :: rc + + integer :: status + real, allocatable :: send_data(:), recv_data(:) + integer :: pf, down_id, k, pos, total_send, total_recv, mpierr + + cat_in = 0.0 + total_send = route%total_send + allocate(send_data(total_send)) + pos = 0 + do pf = 1, route%n_pfaf_local + down_id = route%downid(pf) + if (down_id == -1) cycle + if (route%minCatch <= down_id .and. down_id <= route%maxCatch) then + ! from local + k = down_id - route%minCatch + 1 + cat_in(k) = cat_in(k) + cat_out(pf) + else ! from the other process + pos = pos + 1 + send_data(pos) = cat_out(pf) + endif + enddo + if (total_send > 0) then + send_data = send_data(route%re_order) + endif + total_recv = route%total_recv + allocate(recv_data(total_recv)) + call MPI_AllToAllv(send_data, route%send_count, route%displ_send, MPI_REAL, & + recv_data, route%recv_count, route%displ_recv, MPI_REAL, & + route%comm, mpierr) + do i = 1, total_recv + down_id = route%to_down(i) + k = down_id - route%minCatch + 1 + cat_in(k) = cat_in(k) + recv_data(i) + enddo + RETURN_(ESMF_SUCCESS) + + end subroutine exchange_water + + end subroutine RUN + + ! ---------------------------------------------------------------------------------------- + + subroutine Finalize(gc, import, export, clock, rc) + type(ESMF_GridComp), intent(inout) :: gc ! Gridded component + type(ESMF_State), intent(inout) :: import ! Import state + type(ESMF_State), intent(inout) :: export ! Export state + type(ESMF_Clock), intent(inout) :: clock ! The clock + integer, optional, intent( out) :: rc ! Error code + + ! !DESCRIPTION: + ! Clean-up. + !EOP + ! ErrLog variables + integer :: status + character(len=ESMF_MAXSTR) :: Iam + character(len=ESMF_MAXSTR) :: comp_name + type (T_RROUTE_STATE), pointer :: route => null() + type (RROUTE_wrap) :: wrap + + ! Begin... + ! Get component's name and setup traceback handle + call ESMF_GridCompget(gc, name=comp_name, rc=status) + VERIFY_(status) + Iam = trim(comp_name) // "::Finalize" + + call ESMF_UserCompGetInternalState ( GC, 'RiverRoute_state',wrap, _RC) + route => wrap%ptr + + CALL ESMF_FieldSMMRelease(routeHandle=route%routeHandle, _RC) + + ! Call Finalize for every child + call MAPL_GenericFinalize(gc, import, export, clock, _RC) + ! End + RETURN_(ESMF_SUCCESS) + + end subroutine Finalize + +end module GEOS_RouteGridCompMod + +! ======================= EOF ========================================================= + diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSroute_GridComp/reservoir.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSroute_GridComp/reservoir.F90 new file mode 100644 index 000000000..493891f47 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSroute_GridComp/reservoir.F90 @@ -0,0 +1,384 @@ +#include "MAPL_Generic.h" + +module reservoirMod + + use MAPL + + implicit none + private + public :: res_state, Reservoir + + !----Reservoir module constants---------- + integer,parameter :: nres=7250 + integer,parameter :: nlake=3917 + + real, parameter :: fac_elec_a = 0.30 ! Coefficient for hydropower calculation + real, parameter :: fac_elec_b = 2.00 ! Exponent for hydropower calculation + real, parameter :: fac_irr_a = 0.01 ! Coefficient for irrigation calculation (arid areas) + real, parameter :: fac_irr_b = 3.00 ! Scaling factor for irrigation (arid areas) + real, parameter :: fac_sup_a = 0.03 ! Coefficient for water supply calculation + real, parameter :: fac_sup_b = 2.00 ! Exponent for water supply calculation + real, parameter :: fac_other_a = 0.20 ! Coefficient for other reservoir types + real, parameter :: fac_other_b = 2.00 ! Exponent for other reservoir types + integer, parameter :: fac_fld = 1 ! Flood control parameter + + real, parameter :: fac_a_slake = 0.003 ! Factor for small lakes + real, parameter :: fac_b_slake = 0.40 ! Exponent for small lakes + real, parameter :: fac_a_llake = 0.01 ! Factor for large lakes + real, parameter :: fac_b_llake = 0.60 ! Exponent for large lakes + real, parameter :: thr_wid_lake = 1.e5 ! Threshold lake width (in m) + + + type RES_STATE !reserver related variables + integer, allocatable :: active_res(:) + integer, allocatable :: active_up( :,:) + integer, allocatable :: type_res( :) + real, allocatable :: cap_res( :) !m3 + real, allocatable :: wid_res( :) !m + integer, allocatable :: fld_res( :) + real, allocatable :: Qfld_thres(:) !m3/s + integer, allocatable :: cat2res( :) + real, allocatable :: qres_acc( :) + contains + procedure :: calc + end type RES_STATE + + interface Reservoir + module procedure new_Reservoir + end interface Reservoir + + !----------------------------------------- + +contains + + !--------------------------------------------------------------------------------- + ! Initialization subroutine for reservoirs + + function new_Reservoir(layout, river_input_file, nall, nc, minCatch, maxCatch, use_res, rc) result (res) + + type(RES_STATE) :: res + type (ESMF_DELayout), intent(in) :: layout + character(len=*), intent(in) :: river_input_file + + ! Define the number of reservoirs (nres) and the number of catchments (nc) + integer, intent(in) :: nall,nc,minCatch,maxCatch + + ! Logical variable to check if reservoirs are used + logical, intent(in) :: use_res + + integer, optional, intent(out) :: rc + + ! --------------- + + integer, allocatable :: active_res(:),type_res(:),fld_res(:),cat2res(:) + real, allocatable :: cap_res(:),Qfld_thres(:) + real, allocatable :: wid_res(:) + + ! Internal arrays for various reservoir-related data + integer,allocatable,dimension(:) :: flag_grand,catid_grand,elec_grand,fld_grand,supply_grand,irr_grand,realuse_grand + integer,allocatable,dimension(:) :: nav_grand,rec_grand,other_grand + integer,allocatable,dimension(:) :: type_res_all,cat2res_all + real,allocatable,dimension(:) :: cap_grand,area_max_res,area_grand,area_res + real,pointer :: buff_global(:)=>NULL(),area_all(:)=>NULL() + integer,pointer :: fld_all(:)=>NULL() !buff_global_int(:)=>NULL() + real :: value_max + + integer,allocatable,dimension(:) :: flag_lake,catid_lake + real,allocatable,dimension(:) :: area_lake + + ! Define the flood threshold variable and a counter variable + character(len=2) :: fld_thres + integer :: i,cid, status + type(Netcdf4_fileformatter) :: formatter + + !----------reservoir module-------------- + ! Allocate memory for each array + allocate(flag_grand(nres),catid_grand(nres),active_res(nc),Qfld_thres(nc)) + allocate(elec_grand(nres),type_res(nc),type_res_all(nall),cap_grand(nres),cap_res(nc),area_grand(nres)) + allocate(area_res(nc),area_max_res(nall)) + allocate(fld_grand(nres),fld_res(nc),supply_grand(nres)) + allocate(irr_grand(nres)) + allocate(cat2res(nc),cat2res_all(nall)) + allocate(nav_grand(nres),rec_grand(nres)) + allocate(other_grand(nres)) + allocate(wid_res(nc)) + allocate(realuse_grand(nres)) + + allocate(flag_lake(nlake),catid_lake(nlake),area_lake(nlake)) + + if (MAPL_AM_I_ROOT()) then + call formatter%open(river_input_file, PFIO_READ) + call formatter%get_var('catid_grand', catid_grand, _RC) + call formatter%get_var('flag_grand', flag_grand, _RC) + call formatter%get_var('cap_grand', cap_grand, _RC) + call formatter%get_var('elec_grand', elec_grand, _RC) + call formatter%get_var('fld_grand', fld_grand, _RC) + call formatter%get_var('supply_grand', supply_grand, _RC) + call formatter%get_var('irr_grand', irr_grand, _RC) + call formatter%get_var('nav_grand', nav_grand, _RC) + call formatter%get_var('rec_grand', rec_grand, _RC) + call formatter%get_var('other_grand', other_grand, _RC) + call formatter%get_var('area_grand', area_grand, _RC) + endif + call MAPL_CommsBcast(layout, catid_grand, nres, MAPL_Root, _RC) + call MAPL_CommsBcast(layout, flag_grand, nres, MAPL_Root, _RC) + call MAPL_CommsBcast(layout, cap_grand, nres, MAPL_Root, _RC) + call MAPL_CommsBcast(layout, elec_grand, nres, MAPL_Root, _RC) + call MAPL_CommsBcast(layout, fld_grand, nres, MAPL_Root, _RC) + call MAPL_CommsBcast(layout, supply_grand, nres, MAPL_Root, _RC) + call MAPL_CommsBcast(layout, irr_grand, nres, MAPL_Root, _RC) + call MAPL_CommsBcast(layout, nav_grand, nres, MAPL_Root, _RC) + call MAPL_CommsBcast(layout, rec_grand, nres, MAPL_Root, _RC) + call MAPL_CommsBcast(layout, other_grand, nres, MAPL_Root, _RC) + call MAPL_CommsBcast(layout, area_grand, nres, MAPL_Root, _RC) + + cap_grand=cap_grand*1.e6! Convert capacity from million cubic meters (MCM) to m3 + write(fld_thres,'(I2.2)')fac_fld + area_grand=area_grand*1.e6 ! Convert area from square kilometers (km2) to square meters (m2) + + Qfld_thres=0.!buff_global(minCatch:maxCatch) + + !lake input + if (MAPL_AM_I_ROOT()) then + call formatter%get_var('catid_lake', catid_lake, _RC) + call formatter%get_var('flag_lake', flag_lake, _RC) + call formatter%get_var('area_lake', area_lake, _RC) + call formatter%close(_RC) + endif + call MAPL_CommsBcast(layout, catid_lake, nlake, MAPL_Root, _RC) + call MAPL_CommsBcast(layout, flag_lake, nlake, MAPL_Root, _RC) + call MAPL_CommsBcast(layout, area_lake, nlake, MAPL_Root, _RC) + area_lake=area_lake*1.e6 + + ! Set initial reservoir ID mapping + cat2res_all=0 + do i=1,nres + if(flag_grand(i)==1)then + cid=catid_grand(i) + cat2res_all(cid)=i ! Link reservoirs with catchments: multiple reservoirs in a catchment share attributes that can be accessed via cat2res + endif + enddo + + ! Initialize reservoir properties + cap_res = 0.0 ! Set reservoir capacity to zero + area_res = 0.0 ! Set reservoir area to zero + area_max_res = 0.0 ! Set max reservoir area to zero + type_res_all = 0 ! Set reservoir type to zero + fld_res = 0 ! Set flood status to zero + active_res = 0 ! Set active reservoirs to zero + realuse_grand = 0 ! Initialize real use for each reservoir to zero + + ! Loop over all reservoirs + allocate(buff_global(nall),fld_all(nall),area_all(nall)) + buff_global=0. + area_all=0. + fld_all=0 + do i = 1, nres + if(flag_grand(i) == 1) then ! If the reservoir is flagged as active + cid = catid_grand(i) ! Get the catchment ID for the reservoir + buff_global(cid) = buff_global(cid) + cap_grand(i) ! Sum up the capacities for reservoirs in the same catchment + area_all(cid) = area_all(cid) + area_grand(i) ! Sum up the areas for reservoirs in the same catchment + !Qavg_res(cid) = Qavg_grand(i) ! Assign average flow rate to the catchment + if(fld_grand(i) == 1) fld_all(cid) = 1 ! Mark the catchment if it has flood control + endif + enddo + cap_res=buff_global(minCatch:maxCatch) + value_max=huge(value_max) + where(cap_res==0.) cap_res=value_max + !area_res=buff_global2(minCatch:maxCatch) + fld_res=fld_all(minCatch:maxCatch) + deallocate(buff_global) + + ! Assign reservoir type 6 (Other use) to the largest reservoir in a catchment + do i = 1, nres + if(flag_grand(i) == 1) then + cid = catid_grand(i) + if(other_grand(i) == 1 .and. area_grand(i) >= area_max_res(cid)) then + type_res_all(cid) = 6 + cat2res_all(cid) = i ! Map the catchment to the reservoir + area_max_res(cid) = area_grand(i) ! Update the maximum area for the catchment + endif + endif + enddo + + ! Assign reservoir type 5 (Recreational use) to the largest reservoir in a catchment + do i = 1, nres + if(flag_grand(i) == 1) then + cid = catid_grand(i) + if(rec_grand(i) == 1 .and. area_grand(i) >= area_max_res(cid)) then + type_res_all(cid) = 5 + cat2res_all(cid) = i + area_max_res(cid) = area_grand(i) + endif + endif + enddo + + ! Assign reservoir type 4 (Navigational use) to the largest reservoir in a catchment + do i = 1, nres + if(flag_grand(i) == 1) then + cid = catid_grand(i) + if(nav_grand(i) == 1 .and. area_grand(i) >= area_max_res(cid)) then + type_res_all(cid) = 4 + cat2res_all(cid) = i + area_max_res(cid) = area_grand(i) + endif + endif + enddo + + ! Assign reservoir type 3 (Water supply) to the largest reservoir in a catchment + do i = 1, nres + if(flag_grand(i) == 1) then + cid = catid_grand(i) + if(supply_grand(i) == 1 .and. area_grand(i) >= area_max_res(cid)) then + type_res_all(cid) = 3 + cat2res_all(cid) = i + area_max_res(cid) = area_grand(i) + endif + endif + enddo + + ! Assign reservoir type 2 (Electricity generation) to the largest reservoir in a catchment + do i = 1, nres + if(flag_grand(i) == 1) then + cid = catid_grand(i) + if(elec_grand(i) == 1 .and. area_grand(i) >= area_max_res(cid)) then + type_res_all(cid) = 2 + cat2res_all(cid) = i + area_max_res(cid) = area_grand(i) + endif + endif + enddo + + ! Assign reservoir type 1 (Irrigation) to the largest reservoir in a catchment + do i = 1, nres + if(flag_grand(i) == 1) then + cid = catid_grand(i) + if(irr_grand(i) == 1 .and. area_grand(i) >= area_max_res(cid)) then + type_res_all(cid) = 1 + cat2res_all(cid) = i + area_max_res(cid) = area_grand(i) + endif + endif + enddo + + ! Set up natural lakes + do i = 1, nlake + if(flag_lake(i) == 1 .and. catid_lake(i) > 0) then + cid = catid_lake(i) + if(type_res_all(cid)==0.and.fld_all(cid)==0)then + type_res_all(cid) = -1 !for lake + cat2res_all(cid) = i + area_all(cid) = area_lake(i) + endif + endif + enddo + + type_res=type_res_all(minCatch:maxCatch) + cat2res=cat2res_all(minCatch:maxCatch) + area_res=area_all(minCatch:maxCatch) + ! Compute reservoir width from area (square root of the area) + wid_res = sqrt(area_res)!m + + ! Mark active reservoirs based on type or flood control status + do i = 1, nc + if(type_res(i) /= 0 .or. fld_res(i) == 1) then + active_res(i) = 1 + endif + enddo + + ! Deactivate reservoirs if the use_res flag is set to False + if(use_res .eqv. .False.) active_res = 0 + + allocate(res%active_res, source=active_res) + allocate(res%type_res, source=type_res) + allocate(res%cap_res, source=cap_res) + allocate(res%wid_res, source=wid_res) + allocate(res%fld_res, source=fld_res) + allocate(res%Qfld_thres, source=Qfld_thres) + allocate(res%cat2res, source=cat2res) + + deallocate(flag_grand,catid_grand,elec_grand,type_res_all,cap_grand,area_grand) + deallocate(area_res,area_max_res,fld_grand,supply_grand,irr_grand) + deallocate(cat2res_all,nav_grand,rec_grand,other_grand,realuse_grand) + deallocate(flag_lake,catid_lake,area_lake,area_all,fld_all) + + _RETURN(_SUCCESS) + + end function new_Reservoir + + !--------------------------------------------------------------------------------- + ! Reservoir calculation subroutine + + subroutine calc( this, Qout, Q_res, Wr_res, dt, rc) + + class(RES_STATE), intent(in) :: this + real, intent(in) :: Qout(:) + real, intent(inout) :: Q_res(:), Wr_res(:) + real, intent(in) :: dt + integer, optional, intent(out) :: rc + + ! ------------------------ + + real :: Qin_res, alp_res ! Variables for inflow, coefficients, and factors + integer :: i, status + + ! If the reservoir is active + do i = 1, size(this%active_res) + + if (this%active_res(i) == 1) then + + ! Determine the inflow to the reservoir + Qin_res = Qout(i) ! Inflow from river + + ! Irrigation reservoir + if (this%type_res(i) == 1) then + alp_res = fac_irr_a * ((1.0 / (this%wid_res(i) / 1.e3)) ** fac_irr_b) / 3600.0 ! irrigation coefficient + + ! Hydropower reservoir + else if (this%type_res(i) == 2) then + alp_res = fac_elec_a * ((1.0 / (this%wid_res(i) / 1.e3)) ** fac_elec_b) / 3600.0 ! Hydropower coefficient + + ! Water supply reservoir + else if (this%type_res(i) == 3) then + alp_res = fac_sup_a * ((1.0 / (this%wid_res(i) / 1.e3)) ** fac_sup_b) / 3600.0 ! Supply coefficient + + ! Other reservoir types + else if (this%type_res(i) == 4 .or. this%type_res(i) == 5 .or. this%type_res(i) == 6 .or. this%type_res(i) == 0) then + alp_res = fac_other_a * ((1.0 / (this%wid_res(i) / 1.e3)) ** fac_other_b) / 3600.0 ! Generic reservoir coefficient + + ! Natural lake + else if (this%type_res(i) == -1) then + ! Determine lake type based on area and calculate alpha + if (this%wid_res(i) >= thr_wid_lake) then + alp_res = fac_a_llake * ( (1. / (this%wid_res(i) / 1.e3)) ** fac_b_llake ) / 3600. + else + alp_res = fac_a_slake * ( (1./ (this%wid_res(i) / 1.e3)) ** fac_b_slake ) / 3600. + endif + + endif + + Q_res(i) = alp_res * Wr_res(i) + + ! Ensure outflow is within reasonable bounds + Q_res(i) = max(0.0, Q_res(i)) ! Ensure non-negative outflow + Q_res(i) = min(Q_res(i), Wr_res(i) / dt + Qin_res) ! Limit outflow to prevent exceeding inflow and storage + !if (fld_res == 1) Q_res = min(Q_res, Qfld_thres) ! Limit outflow for flood control + Wr_res(i) = Wr_res(i) + dt * (Qin_res - Q_res(i)) ! Update water storage in the reservoir + Wr_res(i) = max(0.0, Wr_res(i)) ! Ensure non-negative storage + + ! If the storage exceeds capacity, adjust outflow and storage + if (Wr_res(i) > this%cap_res(i)) then + Q_res(i) = Q_res(i) + (Wr_res(i) - this%cap_res(i)) / dt ! Adjust outflow for overflow + Wr_res(i) = this%cap_res(i) ! Limit storage to reservoir capacity + endif + + endif + enddo + + _RETURN(_SUCCESS) + + end subroutine calc + +end module reservoirMod + +! =============== EOF =================================================================== 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 new file mode 100644 index 000000000..e15cdd55a --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSroute_GridComp/routing_model.F90 @@ -0,0 +1,168 @@ +MODULE routing_model + + use MAPL + + IMPLICIT NONE + + private + + real, parameter :: M = 0.45 ! Parameter in hydraulic geometry formula + real, parameter :: mm = 0.35 ! Parameter in hydraulic geometry formula + + public :: river_routing_hyd, SEARCH_DNST, M, mm + +CONTAINS + +! ------------------------------------------------------------------------------------------------------- + + REAL FUNCTION RESCONST (LS, P1, P2) + + IMPLICIT NONE + + REAL, INTENT (IN) :: LS, P1, P2 + + RESCONST = P1 * ((1./LS)**P2) + + END FUNCTION RESCONST + +! ------------------------------------------------------------------------------------------------------- + + RECURSIVE SUBROUTINE SEARCH_DNST (K, NCAT_G, DNST, Pfaf_all, DNST_OUT) + + implicit none + + integer, intent (in) :: NCAT_G, K + integer, intent (in), dimension (NCAT_G) :: Pfaf_all, DNST + integer, intent (inout) :: DNST_OUT + + if (DNST(K) == -1) then + DNST_OUT = -1 + else + + if(Pfaf_all(DNST(K)) >= 1) then + DNST_OUT = Pfaf_all(DNST(K)) + else + if(DNST(DNST(K)) == -1) then + DNST_OUT = -1 + else + call SEARCH_DNST (DNST(DNST(K)), NCAT_G, DNST, Pfaf_all, DNST_OUT) + endif + endif + endif + + RETURN + + END SUBROUTINE SEARCH_DNST + +! ------------------------------------------------------------------------------------------------------- + + ! ====================================================================================== + ! + ! HYDRAULIC GEOMETRY ROUTING MODEL + ! + ! -------------------------------------------------------------------------------------- + ! Routing Model Input Parameters + ! ------------------------------ + !**** NCAT = NUMBER OF CATCHMENTS IN THE STUDY DOMAIN + !**** Qrunf0 = RUNOFF PRODUCED BY LAND SURFACE MODEL IN THE CATCHMENT [m^3/s] + !**** llc_ori = MAIN RIVER LENGTH SCALE [m] + !**** lstr = LOCAL STREAMS LENGTH SCALE [m] + !**** qstr_clmt0= CLIMATOLOGY RUNOFF [m^3/s] + !**** qri_clmt0 = CLIMATOLOGY DISCHAR [m^3/s] + !**** qin_clmt0 = CLIMATOLOGY INFLOW [m^3/s] + !**** K = K PARAMETER FOR MAIN RIVER + !**** Kstr0 = K PARAMETER FOR LOCAL STREAM [m^3/s] + + ! Routing Model Prognostics + ! ------------------------- + !**** Ws0 = AMOUNT OF WATER IN "LOCAL STREAM" [m^3] + !**** Wr0 = AMOUNT OF WATER IN RIVER [m^3] + + ! Routing Model Diagnostics + ! ------------------------- + !**** QS = TRANSFER OF MOISTURE FROM STREAM VARIABLE TO RIVER VARIABLE [m^3/s] + !**** QOUT = TRANSFER OF RIVER WATER TO THE DOWNSTREAM CATCHMENT [m^3/s] + + SUBROUTINE RIVER_ROUTING_HYD ( & + NCAT,ROUTE_DT, & + Qrunf0,llc_ori,lstr, & + qstr_clmt0, qri_clmt0, qin_clmt0, & + K, Kstr0, & + Ws0,Wr0, & + Qs,Qout) + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: NCAT,ROUTE_DT + REAL, INTENT(IN), DIMENSION (NCAT) :: Qrunf0,llc_ori,lstr + REAL, INTENT(IN), DIMENSION (NCAT) :: qstr_clmt0,qri_clmt0,qin_clmt0 + REAL, INTENT(IN), DIMENSION (NCAT) :: K, Kstr0 + REAL, INTENT(INOUT),DIMENSION (NCAT) :: Ws0,Wr0 + REAL, INTENT(OUT), DIMENSION (NCAT) :: Qs,Qout + + + + real, parameter :: small = 1.e-20 + real, parameter :: fac_kstr = 0.01 ! Factor for local stream scaling + + real,dimension(NCAT) :: Qrunf,qstr_clmt,qri_clmt,qin_clmt,Ws,Wr,Kstr + real,dimension(NCAT) :: nume,deno,llc,alp_s,alp_r,Qs0,ks,Ws_last + + real :: dt, rho + + integer :: i,j + + rho = MAPL_RHOWTR + Qrunf = Qrunf0 * rho !m3/s -> kg/s + !llc_ori = llc_ori0 * 1.e3 !km -> m + !lstr = lstr0 * 1.e3 !km -> m + qstr_clmt = qstr_clmt0 * rho !m3/s -> kg/s + qri_clmt = qri_clmt0 * rho !m3/s -> kg/s + qin_clmt = qin_clmt0 * rho !m3/s -> kg/s + Ws = Ws0 * rho !m3 -> kg + Wr = Wr0 * rho !m3 -> kg + Kstr = fac_kstr * Kstr0 + dt = ROUTE_DT + + ! Adjust llc (length of river channel) + nume = qri_clmt**(2.-M) - qin_clmt**(2.-M) ! Numerator for the llc calculation + deno = (2.-M) * (qri_clmt - qin_clmt) * (qri_clmt**(1.-M)) ! Denominator for the llc calculation + where(abs(deno) > small) llc = llc_ori * (nume / deno) ! Compute llc where denominator is not too small + where(abs(deno) <= small) llc = llc_ori * 0.5 ! Set llc to half of original value if denominator is small + + ! Calculate alp_s (stream coefficient) and alp_r (river coefficient) + where(qstr_clmt > small) alp_s = (rho**(-M) * qstr_clmt**(M-mm) * Kstr * (0.5*lstr)**(-1.))**(1./(1.-mm)) ! For non-zero streamflow + where(qstr_clmt <= small) alp_s = 0. ! If streamflow is too small, set alp_s to 0 + where(qri_clmt > small) alp_r = (rho**(-M) * qri_clmt**(M-mm) * K * llc**(-1.))**(1./(1.-mm)) ! For non-zero river input + where(qri_clmt <= small) alp_r = 0. ! If river input is too small, set alp_r to 0 + + ! Update state variables: ks, Ws, and Qs + where(Qrunf<=small)Qrunf=0. ! Set runoff to zero if it's too small + Qs0=max(0.,alp_s * Ws**(1./(1.-mm))) ! Initial flow from stream storage (kg/s) + ks=max(0.,(alp_s/(1.-mm)) * Ws**(mm/(1.D0-mm))) ! Flow coefficient (s^-1) + 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 stream flow (kg/s) + + ! Calculate variables related to river routing: Qr0, kr + Wr=Wr+Qs*dt + Qout=max(0.,alp_r * Wr**(1./(1.-mm))) ! River flow based on water storage (kg/s) + Qout=min(Qout,Wr/dt) + Wr=max(0.,Wr-Qout*dt) + + Ws0 = Ws/rho !kg -> m3 + Wr0 = Wr/rho !kg -> m3 + Qs = Qs/rho !kg/s -> m3/s + Qout = Qout/rho !kg/s -> m3/s + + RETURN + + END SUBROUTINE RIVER_ROUTING_HYD + + ! ------------------------------------------------------------------------------------------------------- + +END MODULE routing_model + +! ======================== EOF ========================================================= diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc index 1ba28b70b..da6252580 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc @@ -187,7 +187,8 @@ # ---- Run river routing model interactively # # 0 : Default - NO -# 1 : YES +# 1 : YES without reservoir model +# 2 : YES with reservoir model # # GEOSagcm=>RUN_ROUTE: 0 # GEOSldas=>RUN_ROUTE: 0 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CMakeLists.txt index 895b59f9d..07bf914eb 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CMakeLists.txt +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CMakeLists.txt @@ -7,6 +7,7 @@ rasterize.F90 read_riveroutlet.F90 CubedSphere_GridMod.F90 rmTinyCatchParaMod.F90 +pfaf_frac.F90 zip.c util.c ) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_shared.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_shared.py index b5e651deb..d252e8044 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_shared.py +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_shared.py @@ -92,6 +92,7 @@ def get_script_mv(grid_type): /bin/mv {GRIDNAME}.j geometry/{GRIDNAME}/. /bin/cp til/{GRIDNAME}{RS}.til geometry/{GRIDNAME}/. /bin/cp til/{GRIDNAME}{RS}.nc4 geometry/{GRIDNAME}/. +/bin/cp til/{GRIDNAME}_tile_pfaf.nc4 geometry/{GRIDNAME}/. if( {TRIPOL_OCEAN} == True ) /bin/cp til/{GRIDNAME}{RS}.TRN geometry/{GRIDNAME}/. /bin/mv rst til geometry/{GRIDNAME}/. @@ -172,6 +173,13 @@ def get_script_mv(grid_type): echo "Successfully copied CO2_MonthlyMean_DiurnalCycle.nc4 to bcs dir." endif +if(-f land/shared/river_input.nc ) then + echo "river_input.nc already present in bcs dir." +else + /bin/cp -p {MAKE_BCS_INPUT_DIR}/../test/stuff/route_model/v2/river_input.nc land/shared/river_input.nc + echo "Successfully copied river_input.nc to bcs dir." +endif + # adjust permissions chmod +rX -R geometry land logs diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 index 5856811b8..8e19a35ae 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 @@ -39,6 +39,7 @@ PROGRAM mkEASETilesParam use MAPL_ConstantsMod, only : MAPL_PI_r8, MAPL_RADIUS use MAPL_ExceptionHandling use MAPL, only : MAPL_ease_extent, MAPL_ease_convert, MAPL_ease_inverse, MAPL_WriteTilingNC4 + use pfaf_fracMod, only : get_pfaf_frac use netcdf implicit none @@ -979,6 +980,8 @@ PROGRAM mkEASETilesParam deallocate( tileid_index, catid_index,veg ) deallocate( tile_area, ease_grid_area, tile_elev, my_land, all_id ) + + call get_pfaf_frac('til/'//trim(EASELabel)//'_tile_pfaf.nc4', MAKE_BCS_INPUT_DIR, EASELabel) ! Commented out "empty" if-block. -rreichle, 15 Jun 2023 ! diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/pfaf_frac.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/pfaf_frac.F90 new file mode 100755 index 000000000..0ccc8f3e5 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/pfaf_frac.F90 @@ -0,0 +1,281 @@ +#include "MAPL_ErrLog.h" + +module pfaf_fracMod +!Main purpose: Assigns a catchment‐tile index from catchment definition files to each model grid cell for M36 grid. +use MAPL +use LogRectRasterizeMod, only: nc=>SRTM_maxcat + +implicit none +private + +public :: get_Pfaf_frac + +contains + + subroutine get_Pfaf_frac(file_out, BCS_PATH, GridName) + character(*), intent(in) :: file_out, BCS_PATH + character(*), intent(in) :: GridName + integer,parameter :: nlon=21600 ! Number of longitude grid points in the original grid + integer,parameter :: nlat=10800 ! Number of latitude grid points in the original grid + ! Variable declarations: + integer :: id, xi, yi, i,j, flag, subi,nmax,nmax2,ntot,it,ns_tot + integer,allocatable :: nsub(:) ! Array to store the number of sub-areas for each catchment + integer, allocatable :: xsub(:,:), ysub(:,:), isub(:,:) + ! xsub and ysub: 2D arrays to store mapped x and y indices for sub-catchments (not using subi_global in this code) + real, allocatable :: asub(:,:) ! 2D array to store aggregated area for each sub-catchment + + real*8, allocatable :: lon(:), lat(:) ! Arrays to hold longitude and latitude values from the NetCDF file + real*8, allocatable :: lons(:), lats(:) ! Arrays to hold longitude and latitude values from the NetCDF file + integer, allocatable :: loni(:), lati(:) + integer, allocatable :: loni2(:), lati2(:) + ! loni and lati: Arrays holding mapping indices from 1-minute resolution data files + integer, allocatable :: catchind(:,:) ! 2D array holding catchment indices for each grid cell + real, allocatable :: cellarea(:,:) ! 2D array containing the area of each grid cell + + ! Declare allocatable arrays for catchment ID, parent catchment ID, and boundary coordinates + integer, allocatable, dimension(:) :: tid, catid, lati_tile, loni_tile + real*8, allocatable, dimension(:) :: lon_left, lon_right, lat_bottom, lat_top, latc,lonc + integer,allocatable,dimension(:,:) :: map_tile + + real,allocatable,dimension(:) :: area_cat, pfaf_frac(:) + real,allocatable,dimension(:,:) :: frac,frac_tile + integer,allocatable,dimension(:) :: nsub_tile + integer,allocatable ::csub(:,:),flag2(:,:), tile_id(:), pfaf_index(:) + type(NetCDF4_FileFormatter) :: formatter + type (FileMetadata) :: meta + type(Variable) :: v + integer :: nc_ease, nr_ease + real :: tmp_lat, tmp_lon + + ! Define file path for input routing data: + character(len=256) :: pfafData_file + character(len=256) :: cellarea_file + + pfafData_file = trim(BCS_PATH)//"/land/topo/v1/SRTM-TopoData/SRTM_PfafData.nc" + cellarea_file = trim(BCS_PATH)//"../test/stuff/route_model/v2/cellarea.nc" + call MAPL_ease_extent( trim(GridName), nc_ease, nr_ease) + + ! Allocate arrays with the specified dimensions: + allocate(catchind(nlon, nlat), cellarea(nlon, nlat)) + allocate(lon(nlon), lat(nlat)) + allocate(loni(nlon), lati(nlat),loni2(nlon), lati2(nlat)) + allocate(lons(nc_ease),lats(nr_ease)) + allocate(nsub(nc)) + + call formatter%open(trim(pfafData_file), PFIO_READ) + call formatter%get_var("latitude", lat) + call formatter%get_var("longitude", lon) + call formatter%get_var("CatchIndex", catchind) + call formatter%close() + + do i = 1, nc_ease + call MAPL_ease_inverse( trim(GridName), real(i-1), 0.0, tmp_lat, tmp_lon) + lons(i) = tmp_lon + enddo + do j = 1, nr_ease + call MAPL_ease_inverse( trim(GridName), 0.0, real(j-1), tmp_lat, tmp_lon) + lats(j) = tmp_lat + enddo + call nearest_index_vector(lat, lats, lati) + call nearest_index_vector(lon, lons, loni) + + call formatter%open(trim(cellarea_file), PFIO_READ) + call formatter%get_var("data", cellarea) + call formatter%close() + cellarea = cellarea / 1.e6 ! Convert cell area (e.g., from m^2 to km^2) + ! Initialize aggregation arrays to zero: + allocate(xsub(9999, nc), ysub(9999, nc)) + nsub = 0 + ! Loop over all grid cells to aggregate cell areas by catchment and sub-area: + do xi = 1, nlon + do yi = 1, nlat + if (catchind(xi, yi) >= 1) then + ! The grid cell belongs to a catchment + id = catchind(xi, yi) ! Get the catchment id for the current cell + flag = 0 ! Reset flag to indicate whether a matching sub-area is found + ! If the catchment already has one or more sub-areas, check for a matching sub-area: + if (nsub(id) >= 1) then + do i = 1, nsub(id) + if (loni(xi) == xsub(i, id) .and. lati(yi) == ysub(i, id)) then + flag = 1 + exit ! Exit the inner loop since a matching sub-area has been found + endif + end do + endif + ! If no matching sub-area was found, create a new sub-area: + if (flag == 0) then + nsub(id) = nsub(id) + 1 + xsub(nsub(id), id) = loni(xi) + ysub(nsub(id), id) = lati(yi) + endif + endif + end do + end do + nmax = maxval(nsub) + !print *,nmax + deallocate(xsub,ysub) + allocate(xsub(nmax, nc), ysub(nmax, nc), asub(nmax, nc)) + ! Initialize aggregation arrays to zero: + nsub = 0;xsub = 0;ysub = 0;asub = 0. + ! Loop over all grid cells to aggregate cell areas by catchment and sub-area: + do xi = 1, nlon + do yi = 1, nlat + if (catchind(xi, yi) >= 1) then + ! The grid cell belongs to a catchment + id = catchind(xi, yi) ! Get the catchment id for the current cell + flag = 0 ! Reset flag to indicate whether a matching sub-area is found + ! If the catchment already has one or more sub-areas, check for a matching sub-area: + if (nsub(id) >= 1) then + do i = 1, nsub(id) + if (loni(xi) == xsub(i, id) .and. lati(yi) == ysub(i, id)) then + flag = 1 + ! If a match is found, accumulate the cell area into the existing sub-area: + asub(i, id) = asub(i, id) + cellarea(xi, yi) + exit ! Exit the inner loop since a matching sub-area has been found + endif + end do + endif + ! If no matching sub-area was found, create a new sub-area: + if (flag == 0) then + nsub(id) = nsub(id) + 1 + xsub(nsub(id), id) = loni(xi) + ysub(nsub(id), id) = lati(yi) + asub(nsub(id), id) = cellarea(xi, yi) + endif + endif + end do + end do + + ! Open the catchment definition file for the M36 grid and read the total number of catchments (header) + open(77, file="clsm/catchment.def");read(77, *) ntot + ! Allocate arrays with size nt + allocate(tid(ntot), catid(ntot), lon_left(ntot), lon_right(ntot), lat_bottom(ntot), lat_top(ntot),latc(ntot),lonc(ntot)) + allocate(lati_tile(ntot),loni_tile(ntot),map_tile(nc_ease,nr_ease)) + ! Loop over each catchment and read: id, catchment id, left/right longitudes, bottom/top latitudes + do i = 1, ntot + read(77, *) tid(i), catid(i), lon_left(i), lon_right(i), lat_bottom(i), lat_top(i) + end do + latc = (lat_bottom + lat_top) / 2. + lonc = (lon_left + lon_right) / 2. + call nearest_index_vector(latc,lats,lati_tile) + call nearest_index_vector(lonc,lons,loni_tile) + map_tile=-9999 + do i =1, ntot + map_tile(loni_tile(i),lati_tile(i)) = i + enddo + + allocate(isub(nmax, nc)) + isub=0 + ! Loop over each catchment + do i = 1, nc + ! Loop over each potential sub-area within the current catchment + do j = 1, nmax + xi = xsub(j, i) ! Retrieve the x-coordinate for the sub-area + yi = ysub(j, i) ! Retrieve the y-coordinate for the sub-area + if (xi /= 0) then ! Check if a valid sub-area exists (non-zero x-coordinate) + isub(j, i) = map_tile(xi, yi) ! Map the sub-area coordinates to a global tile index using map_tile + endif + enddo + enddo + + allocate(area_cat(nc),frac(nmax,nc),nsub_tile(ntot),flag2(nmax,nc)) + where(isub==-9999) asub=0. + area_cat=0. + do i=1,nc + area_cat(i)=sum(asub(:,i)) + enddo + nsub_tile=0 + frac=0. + do i=1,nc + do j=1,nmax + if(isub(j,i)>0)then + it=isub(j,i) + nsub_tile(it)=nsub_tile(it)+1 + frac(j,i)=asub(j,i)/area_cat(i) + endif + if(isub(j,i)==0)exit + enddo + enddo + nmax2=maxval(nsub_tile) + allocate(csub(nmax2,ntot),frac_tile(nmax2,ntot)) + csub=0 + nsub_tile=0 + frac_tile=0. + do i=1,nc + do j=1,nmax + if(isub(j,i)>0)then + it=isub(j,i) + nsub_tile(it)=nsub_tile(it)+1 + csub(nsub_tile(it),it)=i + frac_tile(nsub_tile(it),it)=frac(j,i) + endif + if(isub(j,i)==0)exit + enddo + enddo + ns_tot=sum(nsub_tile) + allocate(tile_id(ns_tot), pfaf_index(ns_tot), pfaf_frac(ns_tot)) + it = 0 + do i=1,ntot + do j=1,nmax2 + if(csub(j,i)>0)then + it = it+1 + tile_id(it) = i + pfaf_index(it) = csub(j,i) + pfaf_frac(it) = frac_tile(j,i) + endif + if(csub(j,i)==0)exit + enddo + enddo + + call meta%add_dimension('tile', ns_tot) + + v = Variable(type=PFIO_INT32, dimensions='tile') + call v%add_attribute('units', '1') + call v%add_attribute('long_name', 'tile_id') + call meta%add_variable('tile_id', v) + + v = Variable(type=pFio_INT32, dimensions='tile') + call v%add_attribute('units', '1') + call v%add_attribute('long_name', 'pfaf_index') + call meta%add_variable('pfaf_index', v) + + v = Variable(type=pFio_REAL32, dimensions='tile') + call v%add_attribute('units', '1') + call v%add_attribute('long_name', 'area_fraction_of_catchment') + call meta%add_variable('pfaf_frac', v) + + call formatter%create(file_out, mode=PFIO_CLOBBER) + call formatter%write(meta) + call formatter%put_var('tile_id', tile_id) + call formatter%put_var('pfaf_index', pfaf_index) + call formatter%put_var('pfaf_frac', pfaf_frac) + call formatter%close() + + end subroutine get_Pfaf_frac + + subroutine nearest_index_vector(candidates, targets, idx) + ! For each targets(k), find argmin_j |candidates(j) - targets(k)| + real*8, intent(in) :: targets(:) + real*8, intent(in) :: candidates(:) + integer, intent(out) :: idx(size(candidates)) ! 1-based indices + integer :: k, j, nT, nC, best_j + real*8 :: best_d, d + + nT = size(targets) + nC = size(candidates) + + do k = 1, nC + best_d = huge(1.0d0) + best_j = 1 + do j = 1, nT + d = abs(candidates(k) - targets(j)) + if (d < best_d) then + best_d = d + best_j = j + end if + end do + idx(k) = best_j ! already 1-based + end do + end subroutine nearest_index_vector + +end module pfaf_fracMod diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/CMakeLists.txt index d4eca3cfc..013d625b5 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/CMakeLists.txt +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/CMakeLists.txt @@ -1 +1 @@ -esma_add_subdirectories (soil routing) +esma_add_subdirectories (soil routing routing_model) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/CMakeLists.txt index 46bdc0dc2..da311583c 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/CMakeLists.txt +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/CMakeLists.txt @@ -16,6 +16,7 @@ set (exe_srcs esma_add_library (${this} SRCS ${srcs} + DEPENDENCIES MAPL GEOS_SurfaceShared GEOS_LandShared ESMF::ESMF NetCDF::NetCDF_Fortran GEOS_CatchCNShared ) foreach (src ${exe_srcs}) string (REGEX REPLACE ".f90" ".x" exe ${src}) @@ -23,6 +24,8 @@ foreach (src ${exe_srcs}) TARGET ${exe} SOURCES ${src} NOINSTALL + LIBS MAPL GFTL_SHARED::gftl-shared GEOS_SurfaceShared GEOSroute_GridComp GEOS_LandShared GEOS_CatchCNShared + ${this} LIBS ${this} NetCDF::NetCDF_Fortran MPI::MPI_Fortran) endforeach () diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/routing_constant.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/routing_constant.f90 index c4fbfa3cb..ff47945c4 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/routing_constant.f90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/routing_constant.f90 @@ -1,5 +1,6 @@ module routing_constant - + + use catch_constants, ONLY: nc => CATCH_N_PFAFS ! hardwired constants for GEOS river routing scheme implicit none @@ -15,7 +16,6 @@ module routing_constant integer,parameter :: loni_max = 20400 ! = loni_min + nlon_G - 1, ind of lon end of Greenland grid in world grid (30 sec resolution) integer,parameter :: lati_min = 16801 ! ind of lat start of Greenland grid in world grid (30 sec resolution) integer,parameter :: lati_max = 21600 ! = lati_min + nlat_G - 1, ind of lat end of Greenland grid in world grid (30 sec resolution) - integer,parameter :: nc = 291284 ! number of catchments in land integer,parameter :: ng = 525 ! number of catchments in Greenland integer,parameter :: nl = 22116 ! number of outlets to ocean in land (not including Greenland) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/CMakeLists.txt new file mode 100644 index 000000000..0956f57f7 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/CMakeLists.txt @@ -0,0 +1,48 @@ +esma_set_this () + +set(srcs + constant.f90 + river_read.f90 + k_module_cali.f90 +) + +set (exe_srcs + get_Pfaf_file.f90 + get_num_sub_catchment.f90 + get_lonlat_bond.f90 + get_isub.f90 + get_area.f90 + get_Qr_clmt.f90 + get_river_length.f90 + get_K_model_calik.f90 + read_input_TopoCat.f90 +) + +esma_add_library (${this} + SRCS ${srcs} + DEPENDENCIES MAPL GEOS_SurfaceShared GEOS_LandShared ESMF::ESMF NetCDF::NetCDF_Fortran GEOS_CatchCNShared) + +foreach (src ${exe_srcs}) + string (REGEX REPLACE ".f90" ".x" exe ${src}) + ecbuild_add_executable ( + TARGET ${exe} + SOURCES ${src} + NOINSTALL + LIBS MAPL GFTL_SHARED::gftl-shared GEOS_SurfaceShared GEOSroute_GridComp GEOS_LandShared GEOS_CatchCNShared + ${this} + LIBS ${this} NetCDF::NetCDF_Fortran MPI::MPI_Fortran) +endforeach () + +# copy to the build directory +set(SCRIPT_FILES + create_river_input.py + get_latloni_cellarea.py + process_lake_data.py + get_dam_data.py + get_lonlati_maptile.py + run_routing_preproc.py +) + +foreach (SCRIPT ${SCRIPT_FILES}) + configure_file(${SCRIPT} ${SCRIPT}) +endforeach() diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/constant.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/constant.f90 new file mode 100644 index 000000000..e1658f09a --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/constant.f90 @@ -0,0 +1,36 @@ +module constant +!module for constants used in the river routing pre-processing package + use catch_constants, ONLY: nc => CATCH_N_PFAFS + + implicit none + public + + ! Define constant parameters + integer,parameter :: nmax09=458 ! Maximum number of sub-catchments per catchment for M09 + integer,parameter :: nlon=21600 ! Number of longitude grid points in the original grid + integer,parameter :: nlat=10800 ! Number of latitude grid points in the original grid + integer,parameter :: nlat09=1624, nlon09=3856 ! Dimensions for the M09 grid + integer,parameter :: nt_global09=1684725 ! Total number of global tiles for area mapping for M09 + ! Define grid dimensions for 15-second resolution data (HydroSHEDS high-res grid) + integer,parameter :: nlonh = 86400 + integer,parameter :: nlath = 33600 + + integer,parameter :: nl_USGS = 3352492 ! Total number of USGS data records + integer,parameter :: nt09=1684725 !Total number of catchment gridcell in M09 + integer,parameter :: nupmax = 34 ! Maximum number of upstream catchments to record + + !river curve parameters + real,parameter :: cur_avg = 1.4 + real,parameter :: cur_min = 0.5 + real,parameter :: cur_max = 5.0 + + integer,parameter :: nga = 9067 ! Number of GAGE-II records + + !lake input parameters: + integer,parameter :: no = 1459201 ! Total number of outlet records in the outlet files + integer,parameter :: nvl = 3409 ! Number of lakes that pass the filtering criteria (area >= 50) + integer,parameter :: nvo = 3917 ! Number of outlet records after matching with lakes + integer,parameter :: nl_lake = 1426967 ! Total number of lake records in the input files + + +end module constant \ No newline at end of file diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/create_river_input.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/create_river_input.py new file mode 100644 index 000000000..aac54a333 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/create_river_input.py @@ -0,0 +1,89 @@ +import numpy as np +from netCDF4 import Dataset + +# Dimensions +nc = 291284 +nr = 7250 +nl = 3917 + +def read_ascii(filename, count, dtype=float): + return np.loadtxt(filename, dtype=dtype, max_rows=count) + +# ---------- Catchment datasets ---------- +Kstr_catchment = read_ascii("temp/Pfaf_Kstr_PR_fac1_0p35_0p45_0p2_n0p2.txt", nc) +Kv_catchment = read_ascii("temp/Pfaf_Kv_PR_0p35_0p45_0p2_n0p2.txt", nc) +area_catchment = read_ascii("temp/Pfaf_area.txt", nc) +lriv_catchment = read_ascii("temp/Pfaf_lriv_PR.txt", nc) +lstr_catchment = read_ascii("temp/Pfaf_lstr_PR.txt", nc) +Qin_catchment = read_ascii("temp/Pfaf_qin.txt", nc) +Qri_catchment = read_ascii("temp/Pfaf_qri.txt", nc) +Qstr_catchment = read_ascii("temp/Pfaf_qstr.txt", nc) +tosink_catchment = read_ascii("temp/Pfaf_tosink.txt", nc, dtype=int) +dnid_catchment = read_ascii("temp/downstream_1D_new_noadj.txt", nc, dtype=int) + +# ---------- Reservoir datasets ---------- +area_reservoir = read_ascii("temp/area_skm_grand.txt", nr) +capmax_reservoir = read_ascii("temp/cap_max_grand.txt", nr) +catid_reservoir = read_ascii("temp/catid_dam_corr_aca_grand5000.txt", nr, dtype=int) +flag_reservoir = read_ascii("temp/flag_all_res.txt", nr, dtype=int) +fldmainsec_reservoir = read_ascii("temp/fldmainsec_grand.txt", nr, dtype=int) +elec_reservoir = read_ascii("temp/hydroelec_grand.txt", nr, dtype=int) +irr_reservoir = read_ascii("temp/irr_grand.txt", nr, dtype=int) +nav_reservoir = read_ascii("temp/nav_grand.txt", nr, dtype=int) +other_reservoir = read_ascii("temp/other_grand.txt", nr, dtype=int) +rec_reservoir = read_ascii("temp/rec_grand.txt", nr, dtype=int) +supply_reservoir = read_ascii("temp/watersupply_grand.txt", nr, dtype=int) + +# ---------- Lake datasets ---------- +catid_lake = read_ascii("temp/lake_outlet_catid.txt", nl, dtype=int) +flag_lake = read_ascii("temp/lake_outlet_flag_valid_2097.txt", nl, dtype=int) +area_lake = read_ascii("temp/lake_outlet_lakearea.txt", nl) + +# ---------- Create NetCDF ---------- +fout = Dataset("output/river_input.nc", "w", format="NETCDF4") + +# define dimensions +fout.createDimension("pfaf", nc) +fout.createDimension("reservoir", nr) +fout.createDimension("lake", nl) + +# ---------- Create variables ---------- +def create_var(name, data, dims, units, long_name): + var = fout.createVariable(name, data.dtype, dims) + var[:] = data + var.units = units + var.long_name = long_name + return var + +# Catchment variables +create_var("Kstr", Kstr_catchment, ("pfaf",), "1", "K parameter for local streams") +create_var("K", Kv_catchment, ("pfaf",), "1", "K parameter for main rivers") +create_var("lengsc", lriv_catchment, ("pfaf",), "km", "main river length scale") +create_var("lstr", lstr_catchment, ("pfaf",), "km", "local streams length scale") +create_var("qin_clmt", Qin_catchment, ("pfaf",), "m3/s", "climatology of catchment inflow") +create_var("qri_clmt", Qri_catchment, ("pfaf",), "m3/s", "climatology of catchment outflow") +create_var("qstr_clmt", Qstr_catchment, ("pfaf",), "m3/s", "climatology of catchment stream flow") +create_var("downid", dnid_catchment, ("pfaf",), "1", "downstream catchment id") +create_var("area_catch", area_catchment, ("pfaf",), "km2", "catchment area") + +# Reservoir variables +create_var("area_grand", area_reservoir, ("reservoir",), "km2", "reservoir area") +create_var("cap_grand", capmax_reservoir, ("reservoir",), "million m3", "reservoir maximum capacity") +create_var("catid_grand", catid_reservoir, ("reservoir",), "1", "catchment id for reservoir") +create_var("flag_grand", flag_reservoir, ("reservoir",), "1", "active flag for reservoir") +create_var("fld_grand", fldmainsec_reservoir, ("reservoir",), "1", "flood-control flag") +create_var("elec_grand", elec_reservoir, ("reservoir",), "1", "hydro-power flag") +create_var("irr_grand", irr_reservoir, ("reservoir",), "1", "irrigation flag") +create_var("nav_grand", nav_reservoir, ("reservoir",), "1", "navigation flag") +create_var("other_grand", other_reservoir, ("reservoir",), "1", "other use flag") +create_var("rec_grand", rec_reservoir, ("reservoir",), "1", "recreation flag") +create_var("supply_grand", supply_reservoir, ("reservoir",), "1", "water supply flag") + +# Lake variables +create_var("catid_lake", catid_lake, ("lake",), "1", "catchment id for lake outlet") +create_var("flag_lake", flag_lake, ("lake",), "1", "active flag for lake") +create_var("area_lake", area_lake, ("lake",), "km2", "lake area") + +# Close file +fout.close() +print("river_input.nc created successfully!") diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_K_model_calik.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_K_model_calik.f90 new file mode 100755 index 000000000..301c88d01 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_K_model_calik.f90 @@ -0,0 +1,124 @@ +program main +!Main purpose: Calculates the K parameter used in the river routing model. + + use k_module ! Import custom module "k_module" which contains necessary subroutines and functions + use constant, only: nl=>nl_USGS,nlat,nlon,nc + use routing_model, only: MU=>M, mm=>mm + implicit none + + ! Declare variables and allocatable arrays + integer, allocatable :: lati(:), loni(:) ! Arrays to store grid indices (latitude and longitude) for stations + real, allocatable :: data(:, :) ! 2D data array to store USGS data + integer, allocatable :: catid_full(:), catid(:) ! Arrays to store full and filtered catchment ids for stations + real, allocatable, dimension(:) :: vel, dis ! Arrays to store velocity and distance data from the USGS dataset + integer, allocatable :: nv(:), flag_gageii(:) ! Arrays for number of values per station and GAGE-II validation flags + real, allocatable :: Qclmt_full(:), slp_full(:), KKobs_full(:), KImodel_full(:) + ! climatology discharge (Qclmt), slope (slp), observed K factor (KKobs), and initial appromaxtion of modeled K factor (KImodel) + real, allocatable :: Qclmt(:), slp(:), KKobs(:), KImodel(:) + ! Arrays for filtered parameters after station selection + real, allocatable :: KImodel_all(:) ! Array to store all modeled K values from the dataset + real, allocatable :: lats_full(:), lons_full(:) ! Arrays to store the full latitude and longitude values for the grid + real*8, allocatable :: MU_axis(:), slp_axis(:), clmt_axis(:), p_axis(:) + ! Arrays for parameter axes: M factor (MU), slope exponent (slp_axis), runoff exponent (clmt_axis), and an additional parameter axis (p_axis) + + ! Set model parameters and scaling factors + real :: exp_slp = 0.2, exp_clmt = -0.2 ! Base model parameters (mm, MU, and exponents for slope and climatology discharge) + !real :: mm=0.4, MU=0.1, exp_slp=0.5, exp_clmt=0.2 ! Alternative parameter set (commented out) + + ! Declare additional integer and real variables for looping and statistical calculations + integer :: nt, ns, np, i, j, k, p, count + real :: ccr(10,10,10), rms(10,10,10) ! 3D arrays to hold correlation coefficients and RMS errors over a parameter space + !real :: ccr(20,10), rms(20,10) ! Alternative array dimensions (commented out) + real :: ccrp, rmsp ! Variables to store computed correlation coefficient and RMS error for a given parameter set + + character(len=900) :: file_vel + character(len=900) :: file_dis + character(len=900) :: file_usid + character(len=900) :: file_lats + character(len=900) :: file_lons + character(len=900) :: file_lat1m + character(len=900) :: file_lon1m + character(len=900) :: file_pfafmap + character(len=900) :: file_gage_id + character(len=900) :: file_gage_acar + + if (command_argument_count() /= 10) then + print *, "no appropriate files found" + stop + endif + call get_command_argument(1, file_vel) + call get_command_argument(2, file_dis) + call get_command_argument(3, file_usid) + call get_command_argument(4, file_lats) + call get_command_argument(5, file_lons) + call get_command_argument(6, file_lat1m) + call get_command_argument(7, file_lon1m) + call get_command_argument(8, file_pfafmap) + call get_command_argument(9, file_gage_id) + call get_command_argument(10, file_gage_acar) + + ! Read USGS data and process it + call read_usgs_data(file_vel, file_dis, nl, data) ! Read USGS data (nl records) into the 2D array "data" + call process_usgs_data(file_usid, nl, ns, data, nv, nt, vel, dis) + ! Process the USGS data to extract the number of stations (ns), velocity, distance + + ! Determine the nearest grid coordinates for each station based on the full grid latitude and longitude arrays + call find_nearest_coords(file_lats, file_lons, file_lat1m, file_lon1m, ns, nlat, nlon, lats_full, lons_full, lati, loni) + + ! Allocate arrays for parameter axes (each with 10 discrete values) + allocate(MU_axis(10), slp_axis(10), clmt_axis(10)) + ! Initialize the correlation and RMS error arrays with a default invalid value + ccr = -9999. + rms = -9999. + count = 0 + ! Set up the parameter axis for MU (M factor): values from 0 to 0.45 in increments of 0.05 + do k = 1, 10 + MU_axis(k) = (k - 1) * 0.05 + enddo + ! Set up the parameter axis for slope exponent: values from 0 to 0.9 in increments of 0.1 + do i = 1, 10 + slp_axis(i) = (i - 1) * 0.1 + enddo + ! Set up the parameter axis for climate exponent: values from -0.8 to 1.2 in increments of 0.2 + do j = 1, 10 + clmt_axis(j) = (j - 1) * 0.2 - 0.8 + enddo + + !do k=1,10 + !do i=1,10 + !do j=1,10 + + ! count = count + 1 ! Increment the count of parameter combinations (currently only one iteration) + + !MU = MU_axis(k) + !exp_slp = slp_axis(i) + !exp_clmt = clmt_axis(j) + +! print *, "count=", count + print *, "M=", MU, ", exp_slp=", exp_slp, ", exp_clmt=", exp_clmt + + ! Retrieve station information and associated parameter data based on grid indices and model parameters + call get_station_inf(file_pfafmap, ns, nc, nlat, nlon, lati, loni, catid_full, Qclmt_full, slp_full, KImodel_all, exp_slp, exp_clmt) + ! filtering stations using the GAGE-II dataset criteria + call get_valide_stations_gageii(file_gage_id, file_gage_acar, ns, nc, catid_full, flag_gageii) + ! Perform regression analysis using the USGS data + call regression(nt, vel, dis, nv, ns, Qclmt_full, slp_full, KKobs_full, KImodel_full, exp_slp, exp_clmt, mm, MU) + ! Filter stations based on predefined criteria + call filter_station(nc, ns, np, lats_full, lons_full, Qclmt_full, slp_full, catid_full, KKobs_full, KImodel_full, Qclmt, slp, catid, KKobs, KImodel, flag_gageii) + ! Calculate the modeled K parameter for each station + !call cal_Kmodel(ns, np, nc, MU, exp_slp, exp_clmt, Qclmt, slp, KKobs, KImodel, KImodel_all, catid, catid_full, ccr(k,i,j), rms(k,i,j)) + call cal_Kmodel(ns, np, nc, MU, exp_slp, exp_clmt, Qclmt, slp, KKobs, KImodel, KImodel_all, catid, catid_full, ccrp, rmsp) + + ! Print the computed correlation coefficient and RMS error + print *, "ccr=", ccrp + print *, "rms=", rmsp + + !enddo + !enddo + !enddo + + ! The following calls would write the 3D parameter space results to NetCDF files (currently commented out) + !call create_ncfile_real3d("ccr_clmtxslpxMU_10x10x10_mm0p35.nc", "data", ccr, MU_axis, slp_axis, clmt_axis, 10, 10, 10) + !call create_ncfile_real3d("rms_clmtxslpxMU_10x10x10_mm0p35.nc", "data", rms, MU_axis, slp_axis, clmt_axis, 10, 10, 10) + +end program main \ No newline at end of file diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_Pfaf_file.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_Pfaf_file.f90 new file mode 100644 index 000000000..1b28e2f58 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_Pfaf_file.f90 @@ -0,0 +1,206 @@ +program main +!main purpose: Reads the Pfafstetter code dataset and generates files for the connectivity of catchments in the routing network. + + use constant, only: nc, nupmax + + implicit none + + ! Declare allocatable arrays for routing and Pfafstetter information: + integer, allocatable, dimension(:) :: downid, finalid + real*8, allocatable, dimension(:) :: pfaf ! Pfafstetter number for each catchment + integer, allocatable, dimension(:,:) :: pfaf_digit, upstream + integer*8, allocatable, dimension(:) :: res ! Temporary storage for digit extraction + integer, allocatable, dimension(:) :: pfaf_last, pfaf_msk, code, behind + integer, allocatable, dimension(:) :: first, last, nup, nts, nts_old + real, allocatable, dimension(:) :: pfaf_area, pfaf_acar, pfaf_acar_old + + ! Declare loop and temporary variables: + integer :: i, j, jj, k, p, down, cur, idx, num, ok, samed, did, nmax + integer :: fulli(12), fullj(12) + real :: val(5) + + ! Define file path for input routing data: + character(len=900) :: file_path !"input/Pfafcatch-routing.dat" + + if (command_argument_count() /= 1) then + print *, "no found" + stop + endif + call get_command_argument(1, file_path) + + !--------------------------------------------------------------------------- + ! Read routing data from the input file: + open(77, file=trim(file_path), form="formatted", status="old") + read(77, *) num + + ! Allocate arrays based on the total number of catchments (nc): + allocate(downid(nc), finalid(nc), pfaf(nc), pfaf_digit(nc,12), res(nc), & + pfaf_last(nc), pfaf_msk(nc), pfaf_area(nc)) + allocate(first(nc), last(nc)) + + do i = 1, nc + read(77, *) idx, pfaf(i), val(1:5), pfaf_area(i) + end do + + !--------------------------------------------------------------------------- + ! Separate the Pfafstetter number into its 12 individual digits. + res = int8(pfaf) ! Convert Pfafstetter numbers to 64-bit integers + pfaf_digit(:,1) = res / (int8(10) ** int8(11)) + do i = 2, 12 + res = res - int8(10) ** int8(13-i) * int8(pfaf_digit(:, i-1)) + pfaf_digit(:, i) = res / (int8(10) ** int8(12-i)) + end do + + !--------------------------------------------------------------------------- + ! Determine the positions of the last nonzero digit (pfaf_last) + ! and the position of the last digit that is not 1 (stored in 'last'). + first = 2 ! Initialize 'first' to 2 by default + last = 2 ! Initialize 'last' to 2 by default + do i = 1, nc + do j = 12, 1, -1 + if (pfaf_digit(i, j) /= 0) then + pfaf_last(i) = j + do k = 0, j-1 + if (pfaf_digit(i, j-k) /= 1) then + last(i) = j - k + exit + endif + end do + exit + endif + end do + end do + do i = 1, nc + if (last(i) <= 1) last(i) = 2 + end do + + !--------------------------------------------------------------------------- + ! Determine the position of the final zero that has nonzero digits after it. + do i = 1, nc + do j = last(i), 2, -1 + if (pfaf_digit(i, j) == 0) then + first(i) = j + exit + endif + end do + end do + + !--------------------------------------------------------------------------- + ! Determine the immediate downstream catchment for each catchment. + do i = 1, nc + + if (first(i) > last(i) - 1) then + ! No valid downstream digit exists; mark as terminal (sink) + downid(i) = -1 + else + + allocate(code(1 : last(i) - first(i))) + code = pfaf_digit(i, first(i) : last(i)-1) + if (any(code == 2) .or. any(code == 4) .or. any(code == 6) .or. any(code == 8)) then + ! If any digit in the extracted part is even, then the catchment is non-coastal. + fulli = pfaf_digit(i, :) + do j = i-1, 1, -1 ! Loop backward to find a catchment just downstream of catchment i + ok = 1 + fullj = pfaf_digit(j, :) + samed = 0 + do k = 1, min(pfaf_last(i), pfaf_last(j)) + if (fulli(k) == fullj(k)) then + samed = samed + 1 + else + exit + endif + end do ! End of k loop: number of matching leading digits stored in 'samed' + if (samed + 1 <= pfaf_last(j)) then + ! Check that none of catchment j's remaining digits (after the common part) + ! are even, which would indicate a branching downstream. + allocate(behind(1 : pfaf_last(j) - samed)) + behind = fullj(samed+1 : pfaf_last(j)) + if (any(mod(behind, 2) == 0)) ok = 0 + deallocate(behind) + else + ok = 0 + endif + if (ok == 1) then + downid(i) = j ! Found the immediate downstream catchment for catchment i + exit + endif + end do ! End of j loop + else + downid(i) = -1 ! If extracted digits are not even, mark as sink (or coastal) + endif + deallocate(code) + + endif ! End if for determining downstream catchment for catchment i + + end do + + !--------------------------------------------------------------------------- + ! Write the downstream catchment IDs to an output file: + open(88, file="temp/downstream_1D_new_noadj.txt") + do i = 1, nc + write(88, *) downid(i) + end do + + ! Write catchment areas to an output file: + open(88, file="temp/Pfaf_area.txt") + do i = 1, nc + write(88, *) pfaf_area(i) + end do + + !--------------------------------------------------------------------------- + ! Build an upstream connectivity matrix: + allocate(upstream(nupmax, nc), nup(nc)) + nup = 0 + upstream = -1 + do i = 1, nc + did = downid(i) + if (did >= 1) then + nup(did) = nup(did) + 1 + upstream(nup(did), did) = i + end if + end do + open(88, file="temp/upstream_1D.txt") + do i = 1, nc + write(88, '(34(I8))') upstream(:, i) + end do + open(88, file="temp/Pfaf_upnum.txt") + do i = 1, nc + write(88, *) nup(i) + end do + + !--------------------------------------------------------------------------- + ! Calculate the number of steps (nts) from each catchment to the sink: + allocate(nts(nc), pfaf_acar(nc)) + nts = -9999 + do i = 1, nc + k = 0 + cur = i + do while (downid(cur) /= -1) + k = k + 1 + cur = downid(cur) + end do + nts(i) = k + end do + open(88, file="temp/Pfaf_tosink.txt") + do i = 1, nc + write(88, *) nts(i) + end do + + !--------------------------------------------------------------------------- + ! Aggregate catchment areas along the flow network: + nmax = maxval(nts) + pfaf_acar = pfaf_area + do j = nmax, 1, -1 + do i = 1, nc + if (nts(i) == j) then + did = downid(i) + pfaf_acar(did) = pfaf_acar(did) + pfaf_acar(i) + endif + end do + end do + open(88, file="temp/Pfaf_acar.txt") + do i = 1, nc + write(88, *) pfaf_acar(i) + end do + +end program main diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_Qr_clmt.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_Qr_clmt.f90 new file mode 100644 index 000000000..ccf1a78c9 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_Qr_clmt.f90 @@ -0,0 +1,183 @@ +program main +!Main purpose: Reads SMAP L4 runoff data (2016–2023) from a NetCDF file and computes the climatological mean discharge for each catchment. + + use omp_lib + use river_read + use constant, only : nlat=>nlat09, nlon=>nlon09, nc, nupmax + + implicit none + + ! Define variables: + real, allocatable :: runoff(:,:), qrunf(:), temp(:,:), qri(:), qin(:) + integer, allocatable :: nts(:), downid(:), upstream(:,:) + integer :: i, j, nmax, did + + character(len=900) :: file_path !"input/SMAPL4_OL7000_runoff_mean_2016_2023.nc" + + if (command_argument_count() /= 1) then + print *, "no found" + stop + endif + call get_command_argument(1, file_path) + + ! Allocate arrays for runoff (grid), catchment runoff, and a temporary grid array: + allocate(runoff(nlon, nlat), qrunf(nc), temp(nlon, nlat)) + + ! Read the "mean_runoff_flux" variable from the NetCDF file: + call read_ncfile_real2d(trim(file_path), "mean_runoff_flux", runoff, nlon, nlat) + ! Replace missing values (-9999) with 0: + where(runoff == -9999.) runoff = 0. + + ! Flip the grid vertically (reverse the latitude order) and assign back to runoff: + temp = runoff(:, nlat:1:-1) + runoff = temp + + ! Convert runoff from [mm/s] to [mm/d] + runoff = runoff * 86400. + + ! Map runoff from the M09 grid to catchments using the function M09_to_cat. + ! The result is in kg/s; then convert to m^3/s by dividing by 1.e3. + qrunf = M09_to_cat(runoff, nlon, nlat, nc) ! kg/s + qrunf = qrunf / 1.e3 ! m^3/s + + ! Write catchment runoff (qrunf): + open(88, file="temp/Pfaf_qstr.txt") + do i = 1, nc + write(88, *) qrunf(i) + end do + + ! Allocate arrays for "steps to sink" (nts), downstream id (downid) and aggregated runoff (qri): + allocate(nts(nc), downid(nc), qri(nc)) + ! Read the number of steps to sink for each catchment from file: + open(77, file="temp/Pfaf_tosink.txt") + read(77, *) nts + ! Read the downstream connectivity (immediate downstream catchment id) from file: + open(77, file="temp/downstream_1D_new_noadj.txt") + read(77, *) downid + + ! Get the maximum number of steps among all catchments: + nmax = maxval(nts) + ! Initialize qri with the catchment runoff values: + qri = qrunf + ! Aggregate runoff upstream: For each catchment with a given number of steps j, + ! add its runoff to its downstream catchment. + do j = nmax, 1, -1 + do i = 1, nc + if (nts(i) == j) then + did = downid(i) + qri(did) = qri(did) + qri(i) + endif + end do + end do + + ! Write the aggregated runoff (qri) to file "Pfaf_qri.txt": + open(88, file="temp/Pfaf_qri.txt") + do i = 1, nc + write(88, *) qri(i) + end do + + ! Allocate arrays for upstream connectivity and inlet discharge (qin): + allocate(upstream(nupmax, nc), qin(nc)) + ! Read upstream connectivity information from file "upstream_1D.txt": + open(77, file="temp/upstream_1D.txt") + read(77, *) upstream + ! Initialize qin to -9999: + qin = -9999. + ! For catchments that have upstream connectivity (upstream(1,:) /= -1), + ! set qin as the difference between outlet discharge (qri) and runoff (qrunf); + ! for catchments with no upstream (upstream(1,:) == -1), set qin to half of direct runoff. + where(upstream(1,:) /= -1) qin = qri - qrunf + where(upstream(1,:) == -1) qin = qrunf / 2. + + ! Write the inlet discharge (qin): + open(88, file="temp/Pfaf_qin.txt") + do i = 1, nc + write(88, *) qin(i) + end do + +contains + !------------------------------------------------------------------------------ + ! Function: M09_to_cat + ! Purpose : Maps runoff data from the M09 grid resolution to catchments using + ! sub-area information. It aggregates runoff from sub-areas weighted by + ! their area fractions. + ! + ! Input: + ! runoff - Runoff array of size (nlon, nlat) [in mm/d] + ! nlon - Number of longitude grid cells. + ! nlat - Number of latitude grid cells. + ! ncat - Number of catchments. + ! + ! Output: + ! Qrunf - Runoff mapped to catchments (in kg/s, then converted to m^3/s). + !------------------------------------------------------------------------------ + function M09_to_cat(runoff, nlon, nlat, ncat) result(Qrunf) + integer, intent(in) :: nlon, nlat, ncat ! Grid dimensions and number of catchments + real, intent(in) :: runoff(nlon, nlat) ! Input runoff array at grid resolution + real :: Qrunf(ncat) ! Output catchment runoff array + + real, parameter :: small = 1.e-12 + + ! Define sub-area parameters (same as in the M09 dataset) + integer, parameter :: nmax = 458 ! Maximum number of sub-areas per catchment + + ! Declare allocatable arrays to hold sub-area data: + real, allocatable, dimension(:,:) :: subarea, frac ! subarea: area of each sub-area, frac: fraction of total + integer, allocatable, dimension(:,:) :: subx, suby ! Coordinates of sub-areas in the grid + real, allocatable, dimension(:) :: tot, runfC, fracA ! tot: total catchment area; runfC: aggregated runoff; fracA: fraction sum + integer, allocatable, dimension(:) :: nsub ! nsub: number of sub-areas per catchment + + integer :: i, j, sx, sy ! Loop counters and sub-area grid coordinates + + ! Allocate arrays for sub-area information and total area: + allocate(nsub(nc), subarea(nmax, nc), subx(nmax, nc), suby(nmax, nc), tot(nc)) + + ! Read sub-area data from text files: + open(77, file="temp/Pfaf_nsub_M09.txt"); read(77, *) nsub + open(77, file="temp/Pfaf_asub_M09.txt"); read(77, *) subarea + open(77, file="temp/Pfaf_xsub_M09.txt"); read(77, *) subx + open(77, file="temp/Pfaf_ysub_M09.txt"); read(77, *) suby + open(77, file="temp/Pfaf_area.txt"); read(77, *) tot + + ! Allocate fraction array (fraction of sub-area relative to total catchment area) + allocate(frac(nmax, nc)) + + ! Compute the fraction for each sub-area: + do i = 1, nc + frac(:, i) = subarea(:, i) / tot(i) + end do + + ! Allocate arrays to accumulate runoff and fraction sums per catchment: + allocate(runfC(nc), fracA(nc)) + runfC = 0. ! Initialize aggregated runoff for each catchment to zero + fracA = 0. ! Initialize fraction accumulation to zero + + !$OMP PARALLEL default(shared) private(i,j,sx,sy) + !$OMP DO + ! Loop over all catchments and their sub-areas: + do i = 1, nc + do j = 1, nsub(i) + sy = suby(j, i) ! Get y-coordinate of the sub-area + sx = subx(j, i) ! Get x-coordinate of the sub-area + ! Only consider valid sub-areas (non-zero fraction and valid runoff values) + if (frac(j, i) > 0. .and. runoff(sx, sy) < 1.e14 .and. runoff(sx, sy) >= 0.) then + runfC(i) = runfC(i) + frac(j, i) * runoff(sx, sy) + fracA(i) = fracA(i) + frac(j, i) + endif + end do + end do + !$OMP END DO + !$OMP END PARALLEL + + ! Convert aggregated runoff to kg/s by multiplying by total catchment area (in m²) + ! and dividing by the number of seconds per day (86400): + Qrunf = runfC * (tot * 1.e6) / 86400. + + ! Deallocate allocated arrays to free memory: + deallocate(subarea, subx, suby, tot, frac, & + runfC, fracA, nsub) + + end function M09_to_cat + !------------------------------------------------------------------------------ + +end program main diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_area.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_area.f90 new file mode 100755 index 000000000..c70a876e8 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_area.f90 @@ -0,0 +1,103 @@ +program main +!Main purpose: Gets the area for each catchment-tile for M09 grid. + +use river_read +use constant, only: nmax=>nmax09,nc,nlon,nlat,nlat09,nlon09,nt_global=>nt_global09 + +implicit none +! Require explicit declaration of all variables + +! Declare variables for indices, flags, and temporary storage +integer :: id, xi, yi, i, j, flag, subi, x_m09, y_m09, it +! Allocatable arrays to hold sub-catchment coordinate indices and global sub-catchment information +integer,allocatable :: xsub(:,:), ysub(:,:), subi_global(:,:) + +! Allocatable array to store sub-catchment area data +real,allocatable :: asub(:,:) + +! Allocatable double precision arrays for storing longitude and latitude values from file +real*8,allocatable :: lon(:), lat(:) +! Allocatable integer arrays for mapping longitude and latitude indices +integer,allocatable :: loni(:), lati(:) +! 2D arrays: catchind holds catchment index for each grid cell; map_tile maps M09 grid cells to global indices +integer,allocatable :: catchind(:,:), map_tile(:,:) +! Arrays for cell areas from the original grid, aggregated area on the M grid, and area per global tile +real,allocatable :: cellarea(:,:), area_m09(:,:), area_tile(:) + +! Define file path for input routing data: +character(len=900) :: file_path !"input/CatchIndex.nc" + +if (command_argument_count() /= 1) then + print *, "no found" + stop +endif +call get_command_argument(1, file_path) + +! Allocate arrays for sub-catchment data +! Allocate 2D arrays with dimensions (nmax, nc) for sub-catchment coordinate indices and areas +allocate(xsub(nmax,nc), ysub(nmax,nc), asub(nmax,nc)) + +! Allocate arrays for the catchment index grid and cell area data +allocate(catchind(nlon,nlat), cellarea(nlon,nlat)) +! Allocate 1D arrays for longitude and latitude values +allocate(lon(nlon), lat(nlat)) +! Allocate arrays for integer mappings of longitude and latitude indices +allocate(loni(nlon), lati(nlat)) + +! Read longitude and latitude data from the NetCDF file +call read_ncfile_double1d(trim(file_path), "longitude", lon, nlon) +call read_ncfile_double1d(trim(file_path), "latitude", lat, nlat) +! Read 2D catchment index data from the same file +call read_ncfile_int2d(trim(file_path), "CatchIndex", catchind, nlon, nlat) +! Read cell area data +call read_ncfile_real2d("output/cellarea.nc", "data", cellarea, nlon, nlat) +! Scale cell area values (from m^2 to km^2) +cellarea = cellarea/1.e6 + +! Read mapping indices from text files for the M09 grid conversion +! Read integer latitude indices mapping for each original grid +open(10, file="temp/lati_1m_M09.txt") +read(10, *) lati +! Read integer longitude indices mapping for each original grid +open(11, file="temp/loni_1m_M09.txt") +read(11, *) loni + +! Allocate and initialize the aggregated area array for the M09 grid +allocate(area_m09(nlon09, nlat09)) +area_m09 = 0. +! Loop over the original grid and accumulate cell areas into the M09 grid using mapping indices +do xi = 1, nlon + do yi = 1, nlat + if (catchind(xi,yi) >= 1) then + x_m09 = loni(xi) + y_m09 = lati(yi) + ! For grid cells with a valid catchment index, add their cell area to the corresponding M09 grid cell + area_m09(x_m09, y_m09) = area_m09(x_m09, y_m09) + cellarea(xi,yi) + endif + enddo +enddo + +! Allocate the map_tile array and read its data from a NetCDF file +allocate(map_tile(nlon09, nlat09)) +call read_ncfile_int2d("temp/map_tile_M09.nc", "data", map_tile, nlon09, nlat09) +! Allocate the global area array to hold area data for each tile +allocate(area_tile(nt_global)) +area_tile = -9999. + +! Map the aggregated M09 grid areas to the global tile indices using the map_tile array +do i = 1, nlon09 + do j = 1, nlat09 + it = map_tile(i, j) + if (it > 0) then + area_tile(it) = area_m09(i, j) + endif + enddo +enddo + +! Write the global tile area data to an output text file +open(88, file="temp/area_M09_1d.txt") +do i = 1, nt_global + write(88, *) area_tile(i) +enddo + +end diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_dam_data.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_dam_data.py new file mode 100755 index 000000000..506e12416 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_dam_data.py @@ -0,0 +1,235 @@ +import sys +import numpy as np +from netCDF4 import Dataset +import os +import glob + +#Main purpose: Processes reservoir (dam) data: reads dam locations and usage information from GRanD database. + +# Function to find the nearest index in a coordinate array +def ind_nearest_coord(coord_array1, coord_array2): + """ + Find the index of the nearest value in coord_array2 for each value in coord_array1. + """ + indices = [] + for coord in coord_array1: + index = np.argmin(np.abs(coord_array2 - coord)) + indices.append(index) + return np.array(indices) + +if __name__ == '__main__': + + file_latdam, file_londam, file_lat1m, file_lon1m, file_catmap, file_acadam, file_damcat_manfix, file_dam_manflag, file_dam_use, file_flood = sys.argv[1:11] + # Parameter settings + ns = 7250 #nr + nlat = 10800 + nlon = 21600 + nc = 291284 + thres = 5000.0 + +#----get dam lat lon ind---- + + # Read data from ASCII files + lats = np.loadtxt(file_latdam, dtype=np.float64, max_rows=ns) + lons = np.loadtxt(file_londam, dtype=np.float64, max_rows=ns) + lat1m = np.loadtxt(file_lat1m, dtype=np.float64, max_rows=nlat) + lon1m = np.loadtxt(file_lon1m, dtype=np.float64, max_rows=nlon) + + # For each target coordinate, find the nearest index in the reference array. + lati = ind_nearest_coord(lats, lat1m) + loni = ind_nearest_coord(lons, lon1m) + +#----get dam catchemnt indices---- + # Read the catchment indices file + nc_file = file_catmap + with Dataset(nc_file, 'r') as ncf: + catchind = ncf.variables["CatchIndex"][:] + if np.ma.is_masked(catchind): + catchind = catchind.filled(-9999) + + # Initialize the array to store the output catchment ID values + catid = np.empty(ns, dtype=int) + + # Loop over each index + for i in range(ns): + # For each index, retrieve the value from catchind. + catid[i] = catchind[ lati[i], loni[i] ] + +#----get dam drainage area-- + # Read full dataset for acar(drainage area) and catchment area from ASCII files + acar_all = np.loadtxt("temp/Pfaf_acar.txt", dtype=float, max_rows=nc) + area_all = np.loadtxt("temp/Pfaf_area.txt", dtype=float, max_rows=nc) + + # Initialize arrays to store the selected values for each catchment + acar = np.empty(ns, dtype=float) + area = np.empty(ns, dtype=float) + + # Loop over each catchment index and assign values based on catid + for i in range(ns): + cid = catid[i] + if cid != -9999: + # Subtract 1 from cid to convert 1-based index to 0-based index for Python + acar[i] = acar_all[cid - 1] + area[i] = area_all[cid - 1] + else: + acar[i] = -9999.0 + area[i] = -9999.0 + +#----look for station: model drainage area is too small------ + # Read drainage area from GRAND database + grand = np.loadtxt(file_acadam, dtype=float, max_rows=ns) + + # Initialize lists to store error information + id_error = [] + + # Loop over each catchment index + for i in range(ns): + # At here only care about large-scale dams + if grand[i] > thres: + if acar[i] < 0.8 * grand[i]: + # Append error information; add 1 to i for 1-based indexing + id_error.append(i + 1) + + # Get the number of errors found + ne = len(id_error) + +#----get corrected catid for above station-------------------- + # Read manually corrected catid array from ASCII files. + catid_error = np.loadtxt(file_damcat_manfix, dtype=int, max_rows=ne) + + # Loop over each error index and update catid. + # Note: We subtract 1 from resid_error values to convert from 1-based to 0-based indexing. + for i in range(ne): + rid = id_error[i] + catid[rid - 1] = catid_error[i] + +#----get dam drainage area with corrected catid-------------------- + # Initialize arrays to store the selected acar and area values for each catchment + acar = np.empty(ns, dtype=float) + area = np.empty(ns, dtype=float) + + # Loop over each catchment index + for i in range(ns): + cid = catid[i] + if cid != -9999: + # Adjust for 1-based indexing: subtract 1 when accessing the full dataset arrays + acar[i] = acar_all[cid - 1] + area[i] = area_all[cid - 1] + else: + acar[i] = -9999.0 + area[i] = -9999.0 + +#----look for station: model drainage area is too large------ + model = acar + # we use a list to collect error indices. + id_error = [] + + # Loop over each catchment index + for i in range(ns): + # Check if the model value is greater than the threshold and grand is less than 80% of model + if model[i] > thres: + if grand[i] < 0.8 * model[i]: + # Append 1-based index (i+1) to the error list + id_error.append(i + 1) + + ne = len(id_error) + +#----create flag for all dams------ + + # Three more manual adjustment dams + resid_man = np.array([5179, 289, 7070], dtype=int) + catid_man = np.array([46616, 142851, 199281], dtype=int) + nman = resid_man.size + + # Update specific indices in catid_all with manual adjustments. + for i in range(nman): + catid[resid_man[i] - 1] = catid_man[i] + + # Write the updated catid_all to an ASCII file + np.savetxt("temp/catid_dam_corr_aca_grand5000.txt", catid, fmt='%d') + + # Read dams flag (whether we still need it in the model) for the above uncorrect dams from a manually checked flag file + flag_error = np.loadtxt(file_dam_manflag, dtype=int, max_rows=ne) + + # Initialize flag_all array with ones (default flag value) + flag_all = np.ones(ns, dtype=int) + + # For each error entry, update flag_all at the specified index. + # Adjust id from 1-based to 0-based indexing. + for i in range(ne): + id_val = id_error[i] + flag_all[id_val - 1] = flag_error[i] + + # If drainage area in GRAND is small and also less than 0.5 times model drainage area, we do not need the dam. + for i in range(ns): + if grand[i] < 1.e3: + if grand[i] < 0.5 * acar[i]: + flag_all[i] = 0 + + # If model drainage area is negative, set flag_all to 0 for that dam. + for i in range(ns): + if acar[i] < 0.: + flag_all[i] = 0 + + # Write the final flag_all array to an ASCII file. + np.savetxt("temp/flag_all_res.txt", flag_all, fmt='%d') + +#----get dam main use--------------- + # Define category strings and corresponding output tags + use_string = ["Irrigation", "Hydroelectricity", "Water supply", "Navigation", "Recreation"] + use_out = ["irr", "hydroelec", "watersupply", "nav", "rec"] + nu = len(use_string) + + # Read the main use data as strings from the GRAND file + with open(file_dam_use, "r") as f: + use = [line.strip() for line in f] + if len(use) != ns: + print(f"Warning: expected {ns} lines, but got {len(use)} lines.") + + # For each category in use_string, create a flag array and output the result + for j in range(nu): + flag = np.zeros(ns, dtype=int) + # Set flag to 1 where the use value matches the current category + for i in range(ns): + if use[i] == use_string[j]: + flag[i] = 1 + + # Write the flag array + out_filename = os.path.join("temp", use_out[j] + "_grand.txt") + np.savetxt(out_filename, flag, fmt='%d') +#----flood use-------------------- + # Read the use_irr strings from the GRAND file + with open(file_flood, "r") as f: + use_irr = [line.strip() for line in f] + + # Initialize the flag array with zeros + flag = np.zeros(ns, dtype=int) + + # Loop over each entry and set flag to 1 if the entry is not "NA" + for i in range(ns): + if use_irr[i] != "NA": + flag[i] = 1 + + # Write the flag array + np.savetxt("temp/fldmainsec_grand.txt", flag, fmt='%d') + +#----other use-------------------- + use_out = "other" + + # Read the main use data from the GRAND file + with open(file_dam_use, "r") as f: + use = [line.strip() for line in f] + if len(use) != ns: + print(f"Warning: expected {ns} entries, but got {len(use)} entries.") + + # Initialize the flag array with zeros + flag = np.zeros(ns, dtype=int) + + # Loop over each entry and set flag to 1 if the entry matches the specified categories + for i in range(ns): + if use[i] == "Fisheries" or use[i] == "NA" or use[i] == "Other": + flag[i] = 1 + + # Write the flag array + np.savetxt("temp/" + use_out + "_grand.txt", flag, fmt='%d') + diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_isub.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_isub.f90 new file mode 100755 index 000000000..f9c314c67 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_isub.f90 @@ -0,0 +1,55 @@ +program main +!Main purpose: Assigns a catchment‐tile index from maptile files to each sub-catchment for M09 grid. + +use river_read ! Use custom module for reading NetCDF files +use constant, only: nlat=>nlat09, nlon=>nlon09, nmax=>nmax09, nc +implicit none + +! Declare allocatable arrays for grid mapping and sub-catchment indices +integer,allocatable :: map_tile(:,:), subx(:,:), suby(:,:), subi(:,:) + +! Declare integer variables for looping and temporary storage +integer :: i, x, y, j, it + +! Allocate the map_tile array for the aggregated grid (M09) with dimensions (nlon, nlat) +allocate(map_tile(nlon,nlat)) +! Read the mapping data from a NetCDF file into the map_tile array +call read_ncfile_int2d("temp/map_tile_M09.nc", "data", map_tile, nlon, nlat) + +! Allocate subx, suby, and subi arrays to store sub-catchment coordinate data and indices +allocate(subx(nmax,nc), suby(nmax,nc), subi(nmax,nc)) + +! Open and read the x-coordinates of sub-catchments from a text file into subx +open(77, file="temp/Pfaf_xsub_M09.txt") +read(77, *) subx + +! Open and read the y-coordinates of sub-catchments from a text file into suby +open(77, file="temp/Pfaf_ysub_M09.txt") +read(77, *) suby + +! Initialize the subi array to zero +subi = 0 + +! Loop over each catchment +do i = 1, nc + ! Loop over each possible sub-area within a catchment + do j = 1, nmax + x = subx(j, i) + y = suby(j, i) + ! If the x-coordinate is non-zero, then the sub-area exists + if (x /= 0) then + ! If x exists but y is zero, then there is an error and the program stops + if (y == 0) stop + ! Map the sub-area indices from the aggregated grid using the map_tile array + subi(j, i) = map_tile(x, y) + endif + enddo +enddo + +! Open an output file to write the computed sub-catchment tile indices +open(88, file="temp/Pfaf_isub_M09.txt") +do i = 1, nc + write(88, '(150(i8))') subi(:, i) +enddo + +end program main diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_latloni_cellarea.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_latloni_cellarea.py new file mode 100755 index 000000000..7c5c6caf5 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_latloni_cellarea.py @@ -0,0 +1,130 @@ +import sys +import numpy as np +import os +from netCDF4 import Dataset + +#Main purpose: Computes grid-cell index arrays and per-cell areas for 1-m high-res grid. + +lat09_file, lon09_file, lat1m_file, lon1m_file = sys.argv[1:5] + +lati09_output_file = "temp/lati_1m_M09.txt" +loni09_output_file = "temp/loni_1m_M09.txt" +cellarea_output_file = "output/cellarea.nc" + +# Grid dimensions +nlat1m, nlon1m = 10800, 21600 +nlat09, nlon09 = 1624, 3856 + +# Read data +lat1m = np.loadtxt(lat1m_file, dtype=float, max_rows=nlat1m) +lon1m = np.loadtxt(lon1m_file, dtype=float, max_rows=nlon1m) +lat09 = np.loadtxt(lat09_file, dtype=float, max_rows=nlat09) +lon09 = np.loadtxt(lon09_file, dtype=float, max_rows=nlon09) + +# Define nearest coordinate function +def ind_nearest_coord(coord_array1, coord_array2): + """ + Find the index of the nearest value in coord_array2 for each value in coord_array1. + """ + indices = [] + for coord in coord_array1: + index = np.argmin(np.abs(coord_array2 - coord)) + indices.append(index) + return np.array(indices) + +# Find nearest coordinates +lati09 = ind_nearest_coord(lat1m, lat09) +loni09 = ind_nearest_coord(lon1m, lon09) + +# Save indices to files (1-based index) +np.savetxt(lati09_output_file, lati09 + 1, fmt='%d') +np.savetxt(loni09_output_file, loni09 + 1, fmt='%d') + +# Compute global grid cell area +def area_global_rectilinear_grid(lat, lon, rearth=6371.22): + """ + Calculate the approximate area of each grid cell on a global rectilinear grid. + + Parameters: + lat : numpy.ndarray + Array of latitude values (degrees). + lon : numpy.ndarray + Array of longitude values (degrees). + rearth : float + Earth radius in kilometers. Default is 6371.22 km. + + Returns: + area_grid : numpy.ndarray + 2D array representing the area of each grid cell (km^2). + """ + # Convert degrees to radians + rad = np.pi / 180.0 + rr = rearth * rad + + # Longitude spacing (constant across latitudes) + dlon = rr * (lon[1] - lon[0]) # Assuming uniform spacing in longitude + + # Compute longitude spacing at each latitude (dx) + dx = dlon * np.cos(lat * rad) + + # Handle rounding issues at poles + dx[0] = 0.0 if lat[0] < -89.9999 else dx[0] + dx[-1] = 0.0 if lat[-1] > 89.9999 else dx[-1] + + # Latitude spacing (dy), can be variable + dy = np.zeros_like(lat) + dy[0] = (lat[1] - lat[0]) * rr + dy[1:-1] = (lat[2:] - lat[:-2]) * rr / 2.0 + dy[-1] = (lat[-1] - lat[-2]) * rr + + # Area per latitude band + area_lat = dx * dy + + # Extend latitude areas to all longitudes + area_grid = np.outer(area_lat, np.ones(len(lon))) + + # Total area of all grid cells + area_total = np.sum(area_grid) + + # Total surface area of the sphere + area_sphere = 4.0 * np.pi * (rearth ** 2) + + # Add metadata as a dictionary + metadata = { + "long_name": "area of each grid cell", + "units": "km^2", + "area_total": area_total, + "area_lat": area_lat, + "rearth": rearth, + "area_sphere": area_sphere, + "area_ratio": area_total / area_sphere + } + + return area_grid, metadata + +# Calculate and save cell area +area, metadata = area_global_rectilinear_grid(lat1m, lon1m) +area *= 1.e6 # Convert to m² + +# Remove existing file and write new cell area to NetCDF +if os.path.exists(cellarea_output_file): + os.remove(cellarea_output_file) + +with Dataset(cellarea_output_file, "w", format="NETCDF4") as fout: + # Create dimensions + fout.createDimension("lat", nlat1m) + fout.createDimension("lon", nlon1m) + + # Create variables for lat and lon + lat_var = fout.createVariable("lat", "f4", ("lat",)) + lon_var = fout.createVariable("lon", "f4", ("lon",)) + lat_var[:] = lat1m + lon_var[:] = lon1m + # Assign units attribute to lat and lon + lat_var.units = "degrees_north" + lon_var.units = "degrees_east" + + # Create the area variable + area_var = fout.createVariable("data", "f8", ("lat", "lon")) + area_var[:] = area + area_var.units = "m2" diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_lonlat_bond.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_lonlat_bond.f90 new file mode 100755 index 000000000..bb0dc5c74 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_lonlat_bond.f90 @@ -0,0 +1,62 @@ +program main +!Main purpose: Extracts the latitude/longitude boundaries of each catchment-tile from catchment definition files for M09 grid. + +use constant, only : nt=>nt09 + +implicit none + +! Declare allocatable arrays for catchment ID, parent catchment ID, +! and geographical boundaries (longitude and latitude extents) +integer, allocatable, dimension(:) :: id, catid +real, allocatable, dimension(:) :: lon_left, lon_right, lat_bottom, lat_top + +integer :: i, ntot ! Loop counter and total number of catchments read from file + +! Define file path for input routing data: +character(len=900) :: file_path !"input/catchment_M09.def" + +if (command_argument_count() /= 1) then + print *, "no found" + stop +endif +call get_command_argument(1, file_path) + +! Allocate arrays with size nt +allocate(id(nt), catid(nt), lon_left(nt), lon_right(nt), lat_bottom(nt), lat_top(nt)) + +! Open input file that contains catchment definitions +open(77, file=trim(file_path)) +! Read total number of catchments (ntot) from the file header +read(77, *) ntot + +! Loop over each catchment and read the definitions: +! id, catchment id, left and right longitudes, bottom and top latitudes +do i = 1, nt + read(77, *) id(i), catid(i), lon_left(i), lon_right(i), lat_bottom(i), lat_top(i) +end do + +! Write the left longitude values to a temporary output file +open(88, file="temp/lon_left_M09.txt") +do i = 1, nt + write(88, *) lon_left(i) +end do + +! Write the right longitude values to a temporary output file +open(88, file="temp/lon_right_M09.txt") +do i = 1, nt + write(88, *) lon_right(i) +end do + +! Write the bottom latitude values to a temporary output file +open(88, file="temp/lat_bottom_M09.txt") +do i = 1, nt + write(88, *) lat_bottom(i) +end do + +! Write the upper (top) latitude values to a temporary output file +open(88, file="temp/lat_upper_M09.txt") +do i = 1, nt + write(88, *) lat_top(i) +end do + +end \ No newline at end of file diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_lonlati_maptile.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_lonlati_maptile.py new file mode 100755 index 000000000..9ecf8396d --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_lonlati_maptile.py @@ -0,0 +1,70 @@ +import sys +import numpy as np +from netCDF4 import Dataset +import os +#Main purpose: Assigns a catchment‐tile index from catchment definition files to each model grid cell for M09 grid. + +# Load data +nt = 1684725 +nlat = 1624 +nlon = 3856 +nc = 291284 + +# Read input data from text files +lat_bot = np.loadtxt("temp/lat_bottom_M09.txt", dtype=float) +lat_up = np.loadtxt("temp/lat_upper_M09.txt", dtype=float) +lon_left = np.loadtxt("temp/lon_left_M09.txt", dtype=float) +lon_right = np.loadtxt("temp/lon_right_M09.txt", dtype=float) + +# Calculate the center latitudes and longitudes +latc = (lat_bot + lat_up) / 2.0 +lonc = (lon_left + lon_right) / 2.0 + +# Read latitudes and longitudes for the grid +lat09_file, lon09_file = sys.argv[1:3] + +lat09m = np.loadtxt(lat09_file, dtype=float) +lon09m = np.loadtxt(lon09_file, dtype=float) + +# Find the nearest coordinates +def ind_nearest_coord(coord_array1, coord_array2): + """ + Find the index of the nearest value in coord_array2 for each value in coord_array1. + """ + indices = [] + for coord in coord_array1: + index = np.argmin(np.abs(coord_array2 - coord)) + indices.append(index) + return np.array(indices) + +lati = ind_nearest_coord(latc, lat09m) +loni = ind_nearest_coord(lonc, lon09m) + +# Save the indices to files (1-based index) +np.savetxt("temp/lati_tile_M09.txt", lati + 1, fmt='%d') +np.savetxt("temp/loni_tile_M09.txt", loni + 1, fmt='%d') + +# Initialize the map_tile array +map_tile = np.full((nlat, nlon), -9999, dtype=int) + +# Fill the map_tile with data +for i in range(nt): + map_tile[lati[i], loni[i]] = i + 1 + +# Remove the existing file if it exists +if os.path.exists("temp/map_tile_M09.nc"): + os.remove("temp/map_tile_M09.nc") + +# Create a NetCDF file and write the data +with Dataset("temp/map_tile_M09.nc", "w", format="NETCDF4") as fout: + # Create dimensions + fout.createDimension("lat", nlat) + fout.createDimension("lon", nlon) + + # Create variable to store the map_tile data with fill_value set during creation + map_tile_var = fout.createVariable("data", "i4", ("lat", "lon"), fill_value=-9999) + map_tile_var[:] = map_tile + +# Print a sample of the map_tile data +#print(map_tile[62, 10]) # Corresponds to map_tile(63-1, 11-1) in NCL (1-based to 0-based) + diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_num_sub_catchment.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_num_sub_catchment.f90 new file mode 100755 index 000000000..649a2131f --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_num_sub_catchment.f90 @@ -0,0 +1,113 @@ +program main +!Main purpose: Assigns a catchment‐tile index from catchment definition files to each model grid cell for M09 grid. + +use river_read +use constant, only : nmax=>nmax09, nc, nlon, nlat + +implicit none + +! Variable declarations: +integer :: id, xi, yi, i, flag, subi +integer :: nsub(nc) ! Array storing the number of sub-areas for each catchment + +! Allocatable arrays for sub-catchment information: +integer, allocatable :: xsub(:,:), ysub(:,:) +real, allocatable :: asub(:,:) ! Aggregated area for each sub-catchment + +! Arrays for grid data and mapping: +real*8, allocatable :: lon(:), lat(:) ! Longitude and latitude arrays from NetCDF file +integer, allocatable :: loni(:), lati(:) ! Mapped integer indices from 1-minute resolution files +integer, allocatable :: catchind(:,:) ! 2D array of catchment indices for each grid cell +real, allocatable :: cellarea(:,:) ! 2D array of cell areas + + +! Define file path for input routing data: +character(len=900) :: file_path !"input/CatchIndex.nc" + +if (command_argument_count() /= 1) then + print *, "no found" + stop +endif +call get_command_argument(1, file_path) + +! Allocate arrays based on the defined dimensions: +allocate(xsub(nmax, nc), ysub(nmax, nc), asub(nmax, nc)) +allocate(catchind(nlon, nlat), cellarea(nlon, nlat)) +allocate(lon(nlon), lat(nlat)) +allocate(loni(nlon), lati(nlat)) + +! Read grid longitude, latitude, catchment index, and cell area data from NetCDF files: +call read_ncfile_double1d(trim(file_path), "longitude", lon, nlon) +call read_ncfile_double1d(trim(file_path), "latitude", lat, nlat) +call read_ncfile_int2d(trim(file_path), "CatchIndex", catchind, nlon, nlat) +call read_ncfile_real2d("output/cellarea.nc", "data", cellarea, nlon, nlat) +cellarea = cellarea / 1.e6 ! Convert cell area units (from m^2 to km^2) + +! Read mapped grid indices for 1-minute resolution from text files: +open(10, file="temp/lati_1m_M09.txt") +read(10, *) lati +open(11, file="temp/loni_1m_M09.txt") +read(11, *) loni + +! Initialize aggregation arrays: +nsub = 0 ! Set number of sub-areas per catchment to zero +xsub = 0 ! Initialize x-coordinate array for sub-catchments to zero +ysub = 0 ! Initialize y-coordinate array for sub-catchments to zero +asub = 0. ! Initialize aggregated area values to zero + +! Loop over each 1m grid cell to accumulate cell areas into sub-catchments: +do xi = 1, nlon + do yi = 1, nlat + if (catchind(xi, yi) >= 1) then + ! The cell belongs to a catchment: + id = catchind(xi, yi) ! Retrieve the catchment id for the current cell + flag = 0 ! Reset flag; will be set to 1 if a matching sub-area is found + + ! Check if this catchment already has at least one sub-area: + if (nsub(id) >= 1) then + do i = 1, nsub(id) + ! If the mapped indices of the current cell match an existing sub-area: + if (loni(xi) == xsub(i, id) .and. lati(yi) == ysub(i, id)) then + flag = 1 + ! Accumulate the cell area into the existing sub-area: + asub(i, id) = asub(i, id) + cellarea(xi, yi) + exit ! Exit the loop once the match is found + endif + end do + endif + + ! If no matching sub-area was found, create a new sub-area for this catchment: + if (flag == 0) then + nsub(id) = nsub(id) + 1 + xsub(nsub(id), id) = loni(xi) + ysub(nsub(id), id) = lati(yi) + asub(nsub(id), id) = cellarea(xi, yi) + endif + + endif + end do +end do + +! Open output files to write the aggregated sub-catchment information: +open(50, file="temp/Pfaf_nsub_M09.txt") +open(51, file="temp/Pfaf_xsub_M09.txt") +open(52, file="temp/Pfaf_ysub_M09.txt") +open(53, file="temp/Pfaf_asub_M09.txt") + +! For each catchment, write: +! - Number of sub-areas +! - X indices of sub-areas (formatted in groups of 458 integers) +! - Y indices of sub-areas (formatted similarly) +! - Aggregated area values of sub-areas (formatted as floating-point numbers) +do i = 1, nc + write(50, *) nsub(i) + write(51, '(458(1x,i4))') xsub(:, i) + write(52, '(458(1x,i4))') ysub(:, i) + write(53, '(458(1x,f10.4))') asub(:, i) +end do + +! Print the maximum number of sub-areas found for any catchment and its location: +print *, maxval(nsub) +print *, maxloc(nsub) + +end \ No newline at end of file diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_river_length.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_river_length.f90 new file mode 100755 index 000000000..df2cbbe56 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/get_river_length.f90 @@ -0,0 +1,250 @@ +program main +!Main purpose: Determines main river channel lengths for each catchment by using HydroSHEDS data of distance to sink. + +use river_read +use constant, only : nc, nlon, nlat, nlonh, nlath, cur_avg, cur_min, cur_max + +implicit none + +real*8, allocatable :: lon(:), lat(:) +real, allocatable :: ldn1m(:,:), elev1m(:,:) +integer, allocatable :: catid(:,:), flag_slp(:) + +real*8, allocatable :: lonh(:), lath(:) +real, allocatable :: ldnh(:,:), elev_15s(:,:) + +! Declare arrays to hold routing and catchment characteristics: +real, allocatable, dimension(:) :: lon_dn, lat_dn, lon_up, lat_up, dist_ref, dist_ref2, ldn_min, ldn_max, riv_len, str_len, slp +real, allocatable, dimension(:) :: lon_min, lon_max, lat_min, lat_max, area, elevdiff_ref, elevdiff +integer, allocatable, dimension(:) :: xi_min, yi_min, xi_max, yi_max +integer, allocatable, dimension(:) :: downid + +! Loop indices and temporary variables +integer xi, yi +integer :: num, i, j, cid, did, k +integer :: data1, data12 +real*8 :: data2 +real :: data7, data9, data10 +real :: elev_temp + +character(len=900) :: file_pfafmap !input/SRTM_PfafData.nc +character(len=900) :: file_ldn !input/hyd_glo_ldn_15s.nc +character(len=900) :: file_hyelev !input/hyd_glo_dem_15s.nc +character(len=900) :: file_pfafrout !input/Pfafcatch-routing.dat + + if (command_argument_count() /= 4) then + print *, "no found" + stop + endif + call get_command_argument(1, file_pfafmap) + call get_command_argument(2, file_ldn) + call get_command_argument(3, file_hyelev) + call get_command_argument(4, file_pfafrout) + +!----------------------------------------------------------------------- +! Regrid LDN (length to sink) from HydroSHEDS data + +allocate(ldn1m(nlon, nlat), catid(nlon, nlat)) +allocate(lon(nlon), lat(nlat)) +! Read longitude, latitude, and catchment index data from SRTM Pfaf data +call read_ncfile_double1d(trim(file_pfafmap), "longitude", lon, nlon) +call read_ncfile_double1d(trim(file_pfafmap), "latitude", lat, nlat) +call read_ncfile_int2d(trim(file_pfafmap), "CatchIndex", catid, nlon, nlat) +ldn1m = -1. +where(catid == -9999) ldn1m = -9999. + +! Allocate high-resolution LDN array and read data from HydroSHEDS 15s file +allocate(ldnh(nlonh, nlath)) +call read_ncfile_real2d(trim(file_ldn), "Band1", ldnh, nlonh, nlath) +where(ldnh.lt.4.e9) ldnh = ldnh / 1.e3 ! Convert from meters to kilometers + +! Regrid: For each grid cell in the M09 grid, assign the minimum LDN value from the corresponding high-res block. +do xi = 1, nlon + do yi = 2041, 10440 + if (ldn1m(xi, yi) .ne. -9999.) then + ldn1m(xi, yi) = minval(ldnh(4*xi-3:4*xi, 4*yi-3-8160:4*yi-8160)) + if (ldn1m(xi, yi) .gt. 4.e9) ldn1m(xi, yi) = -1. + end if + end do +end do +print *, maxval(ldn1m) + +! Allocate arrays to store minimum and maximum LDN for each catchment and their corresponding grid indices +allocate(ldn_min(nc), ldn_max(nc), xi_min(nc), yi_min(nc), xi_max(nc), yi_max(nc)) +ldn_min = 1.e20 +ldn_max = -9999. +xi_min = -9999; yi_min = -9999; xi_max = -9999; yi_max = -9999 +do i = 1, nlon + do j = 1, nlat + if (catid(i, j) >= 1) then + cid = catid(i, j) + if (ldn1m(i, j) > 0. .and. ldn1m(i, j) < ldn_min(cid)) then + ldn_min(cid) = ldn1m(i, j) + xi_min(cid) = i + yi_min(cid) = j + endif + if (ldn1m(i, j) > 0. .and. ldn1m(i, j) > ldn_max(cid)) then + ldn_max(cid) = ldn1m(i, j) + xi_max(cid) = i + yi_max(cid) = j + endif + endif + end do +end do +where(ldn_min == 1.e20) ldn_min = -9999 + +!----------------------------------------------------------------------- +! Compute elevation at 1-minute resolution from high-resolution DEM (15s) +allocate(elev_15s(nlonh, nlath), elev1m(nlon, nlat)) +call read_ncfile_real2d(trim(file_hyelev), "Band1", elev_15s, nlonh, nlath) +where(elev_15s > 30000.) elev_15s = 0. +elev1m = 0. +do xi = 1, nlon + do yi = 2041, 10440 + elev1m(xi, yi) = sum(elev_15s(4*xi-3:4*xi, 4*yi-3-8160:4*yi-8160)) / 16. + end do +end do + + +deallocate(ldnh, elev_15s) +!----------------------------------------------------------------------- +! Get reference distances using routing data + +open(77, file=trim(file_pfafrout), form="formatted", status="old") +read(77, *) num +allocate(lon_dn(nc), lat_dn(nc), lon_up(nc), lat_up(nc), dist_ref(nc), dist_ref2(nc)) +allocate(lon_min(nc), lon_max(nc), lat_min(nc), lat_max(nc), area(nc), elevdiff_ref(nc), elevdiff(nc)) + +! Read routing and catchment geometry data from the Pfafcatch routing file +do i = 1, nc + read(77, *) data1, data2, lon_min(i), lon_max(i), lat_min(i), lat_max(i), data7, area(i), data9, data10, elevdiff_ref(i), data12, lon_dn(i), lat_dn(i), lon_up(i), lat_up(i) +end do + +! Compute spherical distances reference +do i = 1, nc + dist_ref(i) = spherical_distance(lon_dn(i), lat_dn(i), lon_up(i), lat_up(i)) + dist_ref2(i) = spherical_distance(lon_min(i), lat_min(i), lon_max(i), lat_max(i)) +end do +where(dist_ref > dist_ref2 .or. dist_ref == 0.) dist_ref = 0.5 * dist_ref2 + +!-------------------------------------------------------------------- +! Get initial guess of river length (riv_len) based on LDN differences and elevation differences + +allocate(riv_len(nc), downid(nc), flag_slp(nc)) +open(77, file="temp/downstream_1D_new_noadj.txt") +read(77, *) downid + +flag_slp = 1 + +riv_len = -9999. +elevdiff = -9999. +do i = 1, nc + if (downid(i) >= 1) then + did = downid(i) + if (.not. (riv_len(did) >= cur_min * dist_ref(did) .and. riv_len(did) <= cur_max * dist_ref(did))) then + riv_len(did) = ldn_min(i) - ldn_min(did) + if (xi_min(i) > 0 .and. xi_min(did) > 0) then + elevdiff(did) = max(0., elev1m(xi_min(i), yi_min(i)) - elev1m(xi_min(did), yi_min(did))) + flag_slp(did) = 1 + else + elevdiff(did) = elevdiff_ref(did) + flag_slp(did) = 0 + endif + else if (flag_slp(did) == 0 .or. elevdiff(did) == 0.) then + riv_len(did) = ldn_min(i) - ldn_min(did) + if (xi_min(i) > 0 .and. xi_min(did) > 0) then + elevdiff(did) = max(0., elev1m(xi_min(i), yi_min(i)) - elev1m(xi_min(did), yi_min(did))) + flag_slp(did) = 1 + else + elevdiff(did) = elevdiff_ref(did) + flag_slp(did) = 0 + endif + endif + endif +end do + +do i = 1, nc + if (riv_len(i) == -9999.) then + riv_len(i) = (ldn_max(i) - ldn_min(i)) * 0.5 + if (xi_min(i) > 0) then + elevdiff(i) = max(0., 0.5 * elev1m(xi_max(i), yi_max(i)) - 0.5 * elev1m(xi_min(i), yi_min(i))) + else + elevdiff(i) = elevdiff_ref(i) + flag_slp(i) = 0 + endif + endif +end do + +k = 0 +do i = 1, nc + if (.not. (riv_len(i) >= cur_min * dist_ref(i) .and. riv_len(i) <= cur_max * dist_ref(i))) then + riv_len(i) = cur_avg * dist_ref(i) + elevdiff(i) = elevdiff_ref(i) + flag_slp(i) = 0 + k = k + 1 + endif +end do +open(88, file="temp/Pfaf_lriv_PR.txt") +do i = 1, nc + write(88, *) riv_len(i) +end do + +!-------------------------------------------------------------------- +! Calculate the length scale of local streams based on catchment area and river length. +allocate(str_len(nc)) +str_len = area / riv_len / 4. * cur_avg +open(88, file="temp/Pfaf_lstr_PR.txt") +do i = 1, nc + write(88, *) str_len(i) +end do +!-------------------------------------------------------------------- +! Calculate the catchment slope from elevation difference and river length. +allocate(slp(nc)) +slp = elevdiff * 1.e-3 / riv_len +where(slp.lt.1.e-5) flag_slp = 0 +where(slp.lt.1.e-5) slp = 1.e-5 +print *, sum(flag_slp) +open(88, file="temp/Pfaf_slope.txt") +do i = 1, nc + write(88, *) slp(i) +end do +print *, minval(slp) +open(88, file="temp/Pfaf_slope_flag.txt") +do i = 1, nc + write(88, *) flag_slp(i) +end do + +!-------------------------------------------------------------------- +contains + +function spherical_distance(lon_dn, lat_dn, lon_up, lat_up) result(distance) + implicit none + !------------------------------------------------------------ + ! Function: spherical_distance + ! Purpose : Calculates the great-circle distance between two geographic + ! points using the Haversine formula. + ! + ! Input: + ! lon_dn, lat_dn - Longitude and latitude of the first point (degrees) + ! lon_up, lat_up - Longitude and latitude of the second point (degrees) + ! + ! Output: + ! distance - Great-circle distance between the two points (kilometers) + !------------------------------------------------------------ + real, intent(in) :: lon_dn, lat_dn ! Coordinates of downstream point + real, intent(in) :: lon_up, lat_up ! Coordinates of upstream point + real :: distance ! Computed distance (km) + real :: R, dlon, dlat, a, c ! Intermediate variables + + R = 6371.0 ! Earth's radius in kilometers + dlon = (lon_up - lon_dn) * (acos(-1.0) / 180.0) ! Delta longitude (radians) + dlat = (lat_up - lat_dn) * (acos(-1.0) / 180.0) ! Delta latitude (radians) + + a = sin(dlat / 2.0)**2 + cos(lat_dn * (acos(-1.0) / 180.0)) * & + cos(lat_up * (acos(-1.0) / 180.0)) * sin(dlon / 2.0)**2 + c = 2.0 * atan2(sqrt(a), sqrt(1.0 - a)) + distance = R * c + +end function spherical_distance + +end program main \ No newline at end of file diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/k_module_cali.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/k_module_cali.f90 new file mode 100755 index 000000000..50daa6538 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/k_module_cali.f90 @@ -0,0 +1,1029 @@ +module k_module +!module for K parameter calculations. + + use river_read + use constant, only: nga + + implicit none + private + public :: read_usgs_data, process_usgs_data, find_nearest_coords, get_station_inf, regression + public :: filter_station, cal_Kmodel, get_valide_stations_gageii + +contains +!------------------------------------------------------------ + subroutine read_usgs_data(file_vel, file_dis, nl, data) + !------------------------------------------------------------ + ! Subroutine: read_usgs_data + ! Purpose : Reads USGS velocity and discharge data from text files + ! and stores the data in a 2D array. + ! + ! Input: + ! nl - Total number of records (lines) to read. + ! + ! Output: + ! data - 2D array (nl x 2) where column 1 contains velocity and + ! column 2 contains discharge. + !------------------------------------------------------------ + character(len=900),intent(in) :: file_vel, file_dis + integer, intent(in) :: nl + real, allocatable, intent(out) :: data(:,:) + + character(len=900) :: filename + character(len=100) :: line + character(len=100) :: x(100) + + integer :: i, j, l, io, k + integer, allocatable :: nv(:) + + + ! Allocate the data array with nl rows and 2 columns + allocate(data(nl, 2)) + + ! Loop over both "velocity" and "discharge" files + do l = 1, 2 + ! Construct the file name + if(l==1)filename = file_vel + if(l==2)filename = file_dis + open(unit=77, file=trim(filename), status='old') + + ! Allocate a temporary array to count the number of valid tokens per line + allocate(nv(nl)) + + ! Read each line from the file + do k = 1, nl + read(77, '(A)', iostat=io) line + if (io /= 0) then + print *, "Error reading line ", k, " from file: ", trim(filename) + exit + endif + + ! Read tokens from the line and store them in array x + do i = 1, 100 + read(line, *, iostat=io) (x(j), j = 1, i) + if (io == -1) then + exit + endif + end do + + ! Record the number of valid tokens read from this line + nv(k) = i - 1 + + ! If valid data is found, extract the first value for the data array; + ! otherwise assign a missing value (-9999) + if (nv(k) >= 1) then + read(x(1), *, iostat=io) data(k, l) + else + data(k, l) = -9999 + end if + end do + + ! Deallocate the temporary token count array + deallocate(nv) + + ! Close the current file + close(77) + end do + end subroutine read_usgs_data + +!------------------------------------------------------------ + subroutine process_usgs_data(file_usid, nl, ns, data, nv, nt, vel, dis) + !------------------------------------------------------------ + ! Subroutine: process_usgs_data + ! Purpose : Processes the raw USGS data by reading unique station IDs, + ! counting valid records per station. + ! + ! Input: + ! nl - Total number of records. + ! data - 2D array with raw velocity and discharge data. + ! + ! Output: + ! ns - Number of unique stations. + ! nv - Array with the count of valid records per station. + ! nt - Total number of valid records. + ! vel - Array of velocity values. + ! dis - Array of discharge values. + !------------------------------------------------------------ + character(len=900),intent(in) :: file_usid + integer, intent(in) :: nl + integer, intent(out) :: ns + real, intent(inout), allocatable :: data(:,:) + real, allocatable, intent(out) :: vel(:), dis(:) + integer, allocatable, intent(out) :: nv(:) + integer, intent(out) :: nt + + character(len=20), allocatable :: id(:) + integer, allocatable :: nu(:) + character(len=20), allocatable :: idu(:) + integer :: i, k, ii + + ! Allocate array to hold station IDs for each record + allocate(id(nl)) + + ! Read station IDs from file "input/USGSID.txt" + open(unit=11, file=trim(file_usid), status="old") + read(11, *) id + close(11) + + ! Initialize station count and count unique IDs + k = 1 + do i = 2, nl + if (.not.(trim(id(i)) == trim(id(i-1)))) then + k = k + 1 + end if + end do + ns = k + allocate(nu(ns), nv(ns)) + allocate(idu(ns)) + + nu(1) = 1 + idu(1) = id(1) + k = 1 + do i = 2, nl + if (trim(id(i)) == trim(id(i-1))) then + nu(k) = nu(k) + 1 + else + k = k + 1 + nu(k) = 1 + idu(k) = id(i) + end if + end do + + ! Write unique station IDs to files (with and without commas) + open(unit=13, file="temp/id_for_site.txt") + do i = 1, ns + write(13, '(A)') trim(idu(i)) // "," + end do + close(13) + open(unit=13, file="temp/id_for_site_nocomma.txt") + do i = 1, ns + write(13, '(A)') trim(idu(i)) + end do + close(13) + + ! Read record + nv = 0 + nv(1) = 1 + k = 1 + ii = 0 + k = k - 1 + do i = 2, nl + if (id(i) == id(i - 1)) then + k = k + 1 + if (data(i,1) <= 0.0) then + k = k - 1 + else if (data(i,2) <= 0.0) then + k = k - 1 + end if + else + nv(ii + 1) = k + k = 1 + ii = ii + 1 + if (data(i,1) <= 0.0) then + k = k - 1 + else if (data(i,2) <= 0.0) then + k = k - 1 + end if + end if + end do + nv(ii + 1) = k + nt = sum(nv) + allocate(vel(nt), dis(nt)) + k = 0 + do i = 1, nl + if (data(i,1) > 0.0 .and. data(i,2) > 0.0) then + k = k + 1 + vel(k) = data(i,1) + dis(k) = data(i,2) + endif + enddo + + ! Deallocate temporary arrays + deallocate(id) + deallocate(nu) + deallocate(idu) + deallocate(data) + + end subroutine process_usgs_data + +!------------------------------------------------------------ + subroutine find_nearest_coords(file_lats, file_lons, file_lat1m, file_lon1m, ns, nlat, nlon, lats, lons, lati, loni) + !------------------------------------------------------------ + ! Subroutine: find_nearest_coords + ! Purpose : For each station, finds the nearest grid point in a 1-minute + ! resolution grid and returns the corresponding indices. + ! + ! Input: + ! ns - Number of stations. + ! nlat - Number of latitude grid points in the high-resolution grid. + ! nlon - Number of longitude grid points in the high-resolution grid. + ! + ! In/Out: + ! lats, lons - Arrays to store the station latitude and longitude values. + ! + ! Output: + ! lati, loni - Arrays of indices corresponding to the nearest grid points. + !------------------------------------------------------------ + character(len=900),intent(in) :: file_lats, file_lons, file_lat1m, file_lon1m + integer, intent(in) :: ns, nlat, nlon + real, allocatable, intent(inout) :: lats(:), lons(:) + integer, allocatable, intent(out) :: lati(:), loni(:) + real, allocatable :: lat1m(:), lon1m(:) + real :: min_dist_lat, min_dist_lon, dist + integer :: i, j, idx_min_lat, idx_min_lon + + ! Allocate output arrays for grid indices for each station + allocate(lati(ns), loni(ns)) + ! Allocate arrays for station coordinates + allocate(lats(ns), lons(ns)) + ! Allocate arrays for 1-minute grid coordinates + allocate(lat1m(nlat), lon1m(nlon)) + + ! Read station latitudes from file "input/lat_for_site_200.txt" + open(unit=10, file=trim(file_lats), status='old') + do i = 1, ns + read(10, *) lats(i) + end do + close(10) + + ! Read station longitudes from file "input/lon_for_site_200.txt" + open(unit=11, file=trim(file_lons), status='old') + do i = 1, ns + read(11, *) lons(i) + end do + close(11) + + ! Read high-resolution latitude grid from file "input/lat_1m.txt" + open(unit=12, file=trim(file_lat1m), status='old') + do i = 1, nlat + read(12, *) lat1m(i) + end do + close(12) + + ! Read high-resolution longitude grid from file "input/lon_1m.txt" + open(unit=13, file=trim(file_lon1m), status='old') + do i = 1, nlon + read(13, *) lon1m(i) + end do + close(13) + + ! For each station, determine the nearest latitude and longitude indices + do i = 1, ns + min_dist_lat = 1.0e20 + min_dist_lon = 1.0e20 + idx_min_lat = -1 + idx_min_lon = -1 + + ! Find nearest latitude index + do j = 1, nlat + dist = abs(lats(i) - lat1m(j)) + if (dist < min_dist_lat) then + min_dist_lat = dist + idx_min_lat = j + end if + end do + lati(i) = idx_min_lat + + ! Find nearest longitude index + do j = 1, nlon + dist = abs(lons(i) - lon1m(j)) + if (dist < min_dist_lon) then + min_dist_lon = dist + idx_min_lon = j + end if + end do + loni(i) = idx_min_lon + end do + + ! Deallocate high-resolution grid arrays + deallocate(lat1m) + deallocate(lon1m) + end subroutine find_nearest_coords +!------------------------------------------------------------ + subroutine get_station_inf(file_pfafmap, ns, nc, nlat, nlon, lati, loni, catid, Qclmt, slp, KImodel_all, exp_slp, exp_clmt) + !------------------------------------------------------------ + ! Subroutine: get_station_inf + ! Purpose : Retrieves station catchment information from a NetCDF file, + ! assigns climate runoff (Qclmt) and slope values for each station, + ! and computes modeled K values for all catchments. + ! + ! Input: + ! ns - Number of stations. + ! nc - Total number of catchments. + ! nlat, nlon - Dimensions of the grid. + ! lati, loni - Grid indices for each station. + ! exp_slp, exp_clmt - Exponents for slope and climatology discharge. + ! + ! Output: + ! catid - Array of catchment IDs for each station. + ! Qclmt - Array of climatology discharge values for stations. + ! slp - Array of slope values for stations. + ! KImodel_all - Array of modeled K values for all catchments. + !------------------------------------------------------------ + character(len=900),intent(in) :: file_pfafmap + integer, intent(in) :: ns, nc, nlat, nlon + integer, intent(in) :: lati(nlon), loni(nlon) + integer, allocatable, intent(out) :: catid(:) + real, allocatable, intent(out) :: Qclmt(:), slp(:) + real, allocatable, intent(out) :: KImodel_all(:) + real, intent(in) :: exp_slp, exp_clmt + + integer, allocatable :: catchind(:,:) + real, allocatable, dimension(:) :: Qclmt_all, slp_all, Kstr_all, Qstr_all + integer :: i + + ! Allocate arrays for the catchment index and station outputs + allocate(catchind(nlon, nlat), catid(ns)) + allocate(Qclmt_all(nc), slp_all(nc)) + allocate(Qclmt(ns), slp(ns)) + allocate(KImodel_all(nc), Kstr_all(nc), Qstr_all(nc)) + + ! Read catchment index data from the NetCDF file "input/SRTM_PfafData.nc" + call read_ncfile_int2d(trim(file_pfafmap), "CatchIndex", catchind, nlon, nlat) + + ! For each station, assign the catchment ID based on its grid location + do i = 1, ns + catid(i) = catchind(loni(i), lati(i)) + end do + + ! Write station catchment IDs to a temporary file + open(88, file="temp/catid_for_site_200.txt") + do i = 1, ns + write(88, *) catid(i) + end do + close(88) + + ! Read climate runoff data from file "output/Pfaf_qri.txt" + open(77, file="temp/Pfaf_qri.txt") + read(77, *) Qclmt_all + where(Qclmt_all < 1.e-8) Qclmt_all = 1.e-8 + ! Read slope data from file "temp/Pfaf_slope.txt" + open(77, file="temp/Pfaf_slope.txt") + read(77, *) slp_all + ! Read clmt discharge data from file "output/Pfaf_qstr.txt" + open(77, file="temp/Pfaf_qstr.txt") + read(77, *) Qstr_all + where(Qstr_all < 1.e-8) Qstr_all = 1.e-8 + + ! For each station, assign Qclmt and slope using the catchment ID + do i = 1, ns + if (catid(i) /= -9999) then + Qclmt(i) = Qclmt_all(catid(i)) + slp(i) = slp_all(catid(i)) + else + Qclmt(i) = -9999 + slp(i) = -9999 + endif + end do + + ! Calculate modeled K values for all catchments + KImodel_all = (Qclmt_all**(exp_clmt)) * (slp_all**(exp_slp)) + + ! Calculate stream K values using the scaling factor + Kstr_all = (Qstr_all**(exp_clmt)) * (slp_all**(exp_slp)) + + ! Write stream K values to an output file + open(88, file="temp/Pfaf_Kstr_PR_fac1_0p35_0p45_0p2_n0p2.txt") + do i = 1, nc + write(88, *) Kstr_all(i) + end do + close(88) + + ! Deallocate temporary arrays + deallocate(catchind, Qclmt_all, slp_all, Kstr_all, Qstr_all) + + end subroutine get_station_inf +!------------------------------------------------------------ + subroutine get_valide_stations_gageii(file_gage_id, file_gage_acar, ns, nc, catid_sta, flag_thres) + !------------------------------------------------------------ + ! Subroutine: get_valide_stations_gageii + ! Purpose : Compares station drainage area with GAGE-II dataset and applies an + ! area ratio threshold to determine valid stations. + ! + ! Input: + ! ns - Number of stations. + ! nc - Total number of catchments. + ! catid_sta - Array of catchment IDs for stations. + ! + ! Output: + ! flag_thres - Array indicating valid stations (1 for valid, 0 otherwise). + !------------------------------------------------------------ + character(len=900),intent(in) :: file_gage_id, file_gage_acar + integer, intent(in) :: ns, nc + integer, intent(in) :: catid_sta(ns) + integer, allocatable, intent(out) :: flag_thres(:) + + real :: thr_sel = 0.3 ! Threshold selection factor + + real, dimension(:), allocatable :: acar_pfaf + integer :: i, j, k, cid + character(len=20) :: id_gages(nga) + character(len=20) :: id_sta(ns) + integer :: flag_gageii(ns) + real :: acar_gages(nga) + real :: acar_gages_sta(ns), acar_sta(ns) + character(len=20) :: line + integer :: ios + + allocate(flag_thres(ns)) + + ! Initialize station area ratios with a missing value + acar_sta = -9999.0 + k = 0 + + ! Read GAGE-II station IDs from "input/id_gagesii.txt" + open(unit=10, file=trim(file_gage_id), status="old", action="read") + do j = 1, nga + read(10, '(A)', iostat=ios) id_gages(j) + if (ios /= 0) then + print *, "Error reading id_gagesii.txt" + stop + end if + end do + close(10) + + ! Read station IDs for the sites from "temp/id_for_site_nocomma.txt" + open(unit=11, file="temp/id_for_site_nocomma.txt", status="old", action="read") + do i = 1, ns + read(11, '(A)', iostat=ios) id_sta(i) + if (ios /= 0) then + print *, "Error reading id_for_site_nocomma.txt" + stop + end if + end do + close(11) + + ! Read area ratios for GAGE-II stations from "input/acar_gagesii.txt" + open(unit=12, file=trim(file_gage_acar), status="old", action="read") + do j = 1, nga + read(12, *, iostat=ios) acar_gages(j) + if (ios /= 0) then + print *, "Error reading acar_gagesii.txt" + stop + end if + end do + close(12) + + ! Initialize the GAGE-II flag array to zero (no match) + flag_gageii = 0 + ! Compare station IDs with GAGE-II IDs and mark matches + do i = 1, ns + do j = 1, nga + if (trim(id_gages(j)) == trim(id_sta(i))) then + acar_gages_sta(i) = acar_gages(j) + flag_gageii(i) = 1 + k = k + 1 + exit ! Exit loop after a match is found + end if + end do + end do + + print *, "Number of matches:", sum(flag_gageii) + + allocate(acar_pfaf(nc)) + open(77, file="temp/Pfaf_acar.txt") + read(77, *) acar_pfaf + close(77) + + ! For each station, assign the area ratio based on its catchment ID + do i = 1, ns + if (catid_sta(i) /= -9999) then + cid = catid_sta(i) + acar_sta(i) = acar_pfaf(cid) + else + acar_sta(i) = -9999. + end if + end do + + ! Apply threshold criteria to flag valid stations + flag_thres = 0 + do i = 1, ns + if (flag_gageii(i) == 1 .and. catid_sta(i) /= -9999) then + if (acar_sta(i) .ge. (1. - thr_sel) * acar_gages_sta(i) .and. & + acar_sta(i) .le. (1. + thr_sel) * acar_gages_sta(i)) then + flag_thres(i) = 1 + endif + endif + end do + + print *, "Number of valid:", sum(flag_thres) + + deallocate(acar_pfaf) + end subroutine get_valide_stations_gageii +!------------------------------------------------------------ + subroutine regression(nt, vel_ori, dis_ori, nv, ns, Qclmt, slp, KKobs, KImodel, exp_slp, exp_clmt, mm, MU) + !------------------------------------------------------------ + ! Subroutine: regression + ! Purpose : For each station with sufficient valid records, performs a + ! regression between discharge and velocity to obtain a calibration + ! factor, and then computes the observed K value (KKobs) for that station. + ! + ! Input: + ! nt - Total number of valid records. + ! ns - Number of stations. + ! nv - Array containing the count of valid records per station. + ! vel_ori - Original velocity data (in ft/s, will be converted). + ! dis_ori - Original discharge data (in ft^3/s, will be converted). + ! Qclmt - Climatology discharge data for each station. + ! slp - Slope data for each station. + ! exp_slp, exp_clmt - Exponents for slope and climatology discharge. + ! mm, MU - Model parameters. + ! + ! Output: + ! KKobs - Array of observed K values for each station. + ! KImodel - Array of modeled K values (init guess) for each station. + !------------------------------------------------------------ + integer, intent(in) :: nt, ns + real, intent(inout), allocatable :: vel_ori(:), dis_ori(:) + integer, intent(in) :: nv(ns) + real, intent(inout), allocatable :: Qclmt(:), slp(:) + real, intent(out), allocatable :: KKobs(:), KImodel(:) + real, intent(in) :: exp_slp, exp_clmt, mm, MU + + real, allocatable, dimension(:) :: x, y, yest + integer :: thres = 100 + integer :: i, j + real :: k(ns), cdtm(ns), med + integer :: acc(ns) + real, allocatable :: vel(:), dis(:) + + ! Convert velocity from ft/s to m/s and discharge from ft^3/s to m^3/s + allocate(vel(nt), dis(nt)) + vel = vel_ori * 0.3048 + dis = dis_ori * 0.0283168 + + ! Calculate cumulative counts to index into the valid records for each station + acc(1) = nv(1) + do i = 2, ns + acc(i) = acc(i - 1) + nv(i) + end do + + ! For each station with enough valid records, perform regression + do i = 1, ns + if (nv(i) >= thres) then + allocate(x(nv(i)), y(nv(i)), yest(nv(i))) + x = dis(acc(i) - nv(i) + 1 : acc(i))**mm + y = vel(acc(i) - nv(i) + 1 : acc(i)) + k(i) = sum(x * y) / sum(x * x) + yest = k(i) * x + cdtm(i) = cal_cdtm(y, yest) + deallocate(x, y, yest) + else + k(i) = -9999. + cdtm(i) = -9999. + endif + end do + med = median(cdtm) + + ! Invalidate calibration factors for stations with low determination coefficient + where(cdtm < 0.5) k = -9999. + + allocate(KKobs(ns)) + do i = 1, ns + if (k(i) /= -9999. .and. Qclmt(i) /= -9999.) then + KKobs(i) = k(i) / (Qclmt(i)**(MU - mm)) + else + KKobs(i) = -9999. + endif + end do + + ! Calculate modeled K values (init guess) using the provided exponents + allocate(KImodel(ns)) + KImodel = (Qclmt**(exp_clmt)) * (slp**(exp_slp)) + + deallocate(vel, dis) + end subroutine regression +!------------------------------------------------------------ + subroutine filter_station(nc, ns, np, lats_full, lons_full, Qclmt_full, slp_full, catid_full, KKobs_full, KImodel_full, Qclmt, slp, catid, KKobs, KImodel, flag_gageii) + !------------------------------------------------------------ + ! Subroutine: filter_station + ! Purpose : Filters out stations that do not meet several criteria: + ! valid catchment ID, valid K values, minimum slope threshold, + ! and a positive GAGE-II flag. It then outputs the filtered data. + ! + ! Input: + ! nc - Total number of catchments. + ! ns - Number of stations. + ! lats_full, lons_full - Full arrays of station latitudes and longitudes. + ! Qclmt_full, slp_full - Full climatology discharge and slope data for stations. + ! catid_full - Full catchment ID array for stations. + ! KKobs_full, KImodel_full - Full observed and modeled K values (initial guess). + ! flag_gageii - GAGE-II validation flags. + ! + ! Output: + ! np - Number of stations that passed the filter. + ! Qclmt, slp, KKobs, KImodel - Filtered arrays for clmt discharge, slope, observed and modeled K (init guess). + ! catid - Filtered catchment IDs for the valid stations. + !------------------------------------------------------------ + integer, intent(in) :: ns, nc + integer, intent(out) :: np + real, intent(inout), allocatable :: lats_full(:), lons_full(:), Qclmt_full(:), slp_full(:), KKobs_full(:), KImodel_full(:) + real, intent(out), allocatable :: Qclmt(:), slp(:), KKobs(:), KImodel(:) + integer, intent(inout), allocatable :: catid_full(:) + integer, intent(out), allocatable :: catid(:) + integer, intent(inout), allocatable :: flag_gageii(:) + + integer, allocatable :: flag_slp(:) + real, allocatable :: lats(:), lons(:) + integer :: i, k + integer, allocatable :: flag_7065(:) + + ! Allocate and read slope flag data from file "temp/Pfaf_slope_flag.txt" + allocate(flag_slp(nc)) + open(77, file="temp/Pfaf_slope_flag.txt") + read(77, *) flag_slp + + allocate(flag_7065(ns)) + flag_7065 = 0 + + k = 0 + ! Count stations that meet all filtering criteria + do i = 1, ns + if (catid_full(i) .ne. -9999 .and. KKobs_full(i) /= -9999. .and. & + slp_full(i) > 1.e-5 .and. flag_slp(catid_full(i)) == 1 .and. & + flag_gageii(i) == 1) then + k = k + 1 + endif + end do + np = k + print *, "number of valid stations: ", np + + ! Allocate filtered output arrays + allocate(Qclmt(np), slp(np), catid(np), KKobs(np), KImodel(np)) + allocate(lats(np), lons(np)) + k = 0 + do i = 1, ns + if (catid_full(i) .ne. -9999 .and. KKobs_full(i) /= -9999. .and. & + slp_full(i) > 1.e-5 .and. flag_slp(catid_full(i)) == 1 .and. & + flag_gageii(i) == 1) then + k = k + 1 + Qclmt(k) = Qclmt_full(i) + slp(k) = slp_full(i) + KKobs(k) = KKobs_full(i) + KImodel(k) = KImodel_full(i) + catid(k) = catid_full(i) + lats(k) = lats_full(i) + lons(k) = lons_full(i) + flag_7065(i) = 1 + endif + end do + + ! Deallocate temporary full arrays that are no longer needed + deallocate(Qclmt_full, slp_full, KKobs_full, KImodel_full, flag_slp, flag_gageii, lats, lons) + + end subroutine filter_station +!------------------------------------------------------------ + subroutine cal_Kmodel(ns, np, nc, MU, exp_slp, exp_clmt, Qclmt, slp, KKobs, KImodel, KImodel_all, catid, catid_full, ccr, rms) + !------------------------------------------------------------ + ! Subroutine: cal_Kmodel + ! Purpose : Calibrates the model by adjusting catchment K values with a scaling + ! factor computed from the percentiles of observed and modeled K values. + ! It then computes the correlation coefficient (ccr) and RMS error. + ! + ! Input/Output: + ! ns - Number of stations. + ! np - Number of valid stations. + ! nc - Total number of catchments. + ! MU, exp_slp, exp_clmt - Model parameters. + ! Qclmt, slp, KKobs, KImodel - Arrays for station data. + ! KImodel_all - Modeled K values for all catchments. + ! catid, catid_full - Filtered and full catchment ID arrays. + ! + ! Output: + ! ccr - Correlation coefficient between observed and calibrated K. + ! rms - RMS error between observed and calibrated K. + !------------------------------------------------------------ + integer, intent(in) :: ns, np, nc + real, intent(in) :: MU, exp_slp, exp_clmt + real, intent(inout), allocatable :: Qclmt(:), slp(:), KKobs(:), KImodel(:) + real, intent(inout), allocatable :: KImodel_all(:) + integer, intent(inout), allocatable :: catid(:), catid_full(:) + real, intent(inout) :: ccr, rms + + real, allocatable :: KKobs_sort(:), KImodel_sort(:), KKmodel_full(:) + real, allocatable, dimension(:) :: dis, sca, Kv, KKmodel + integer, allocatable, dimension(:) :: gear + + character(len=50) :: MU_s, exp_slp_s, exp_clmt_s + + integer :: bulk, i, lev + real :: Kper(11), KMper(11), rat(11), dis_full(11) + + ! Format model parameters into strings for output naming purposes + write(MU_s, '(f4.2)') MU + write(exp_slp_s, '(f4.2)') exp_slp + if (exp_clmt >= 0.) then + write(exp_clmt_s, '(f4.2)') exp_clmt + else + write(exp_clmt_s, '(f4.2)') -1.*exp_clmt + exp_clmt_s = "n" // trim(exp_clmt_s) + endif + + ! Allocate arrays for sorted K values + allocate(KKobs_sort(np), KImodel_sort(np)) + call sort(np, KKobs, KKobs_sort) + call sort(np, KImodel, KImodel_sort) + + ! Compute percentile thresholds by dividing sorted arrays into 10 equal parts + bulk = np / 10 + Kper(1) = KKobs_sort(1) + KMper(1) = KImodel_sort(1) + do i = 2, 10 + Kper(i) = KKobs_sort(bulk * (i - 1)) + KMper(i) = KImodel_sort(bulk * (i - 1)) + end do + Kper(11) = KKobs_sort(np) + KMper(11) = KImodel_sort(np) + rat = Kper / KMper + + ! Allocate arrays for scaling calculations over all catchments + allocate(gear(nc), dis(nc), sca(nc), Kv(nc)) + + ! Initialize gear to default (12) and compute distance to percentile thresholds + gear = 12 + dis = -9999. + do i = 1, nc + do lev = 1, 11 + if (KImodel_all(i) <= KMper(lev)) then + gear(i) = lev + dis(i) = KMper(lev) - KImodel_all(i) + exit + endif + end do + end do + + ! Calculate differences between consecutive percentile thresholds + dis_full(1) = KMper(1) + do i = 2, 11 + dis_full(i) = KMper(i) - KMper(i - 1) + end do + + ! Compute scaling factors for each catchment based on its percentile position + do i = 1, nc + if (gear(i) == 1) then + sca(i) = rat(1) + elseif (gear(i) == 12) then + sca(i) = rat(11) + else + sca(i) = ( rat(gear(i)-1) * dis(i) + rat(gear(i)) * (dis_full(gear(i)) - dis(i)) ) / dis_full(gear(i)) + endif + Kv(i) = KImodel_all(i) * sca(i) + end do + + ! Write scaled K values for each catchment to an output file + open(88, file="temp/Pfaf_Kv_PR_0p35_0p45_0p2_n0p2.txt") + do i = 1, nc + write(88, *) Kv(i) + end do + close(88) + + ! For each station, assign the corresponding scaled K value from its catchment + allocate(KKmodel_full(ns)) + do i = 1, ns + if (catid_full(i) /= -9999) then + KKmodel_full(i) = Kv(catid_full(i)) + else + KKmodel_full(i) = -9999. + endif + end do + + ! For filtered stations, extract the modeled K values + allocate(KKmodel(np)) + do i = 1, np + KKmodel(i) = Kv(catid(i)) + end do + + ! Compute correlation coefficient and RMS error between observed and modeled K values + ccr = cal_ccr(KKobs, KKmodel) + rms = cal_rms(KKobs, KKmodel, np) + + ! Deallocate temporary arrays and full data arrays + deallocate(KKobs_sort, KImodel_sort) + deallocate(KImodel_all, gear, dis, sca, Kv) + deallocate(Qclmt, slp, KKobs, KImodel, catid, KKmodel, catid_full, KKmodel_full) + end subroutine cal_Kmodel + + subroutine sort(np, data, data_sort) + !------------------------------------------------------------ + ! Subroutine: sort + ! Purpose : Sorts an array of real numbers in ascending order using bubble sort. + ! + ! Input: + ! np - Number of elements in the array. + ! data - Input array to be sorted. + ! + ! Output: + ! data_sort - Sorted array. + !------------------------------------------------------------ + integer, intent(in) :: np ! Size of the array + real, intent(in) :: data(np) ! Input array + real, intent(out) :: data_sort(np) ! Output sorted array + integer :: i, j + real :: temp + + ! Copy the input array to the output array + data_sort = data + + ! Bubble sort algorithm + do i = 1, np - 1 + do j = 1, np - i + if (data_sort(j) > data_sort(j + 1)) then + temp = data_sort(j) + data_sort(j) = data_sort(j + 1) + data_sort(j + 1) = temp + end if + end do + end do + end subroutine sort + + function cal_ccr(y, yest) result(ccr) + !------------------------------------------------------------ + ! Function: cal_ccr + ! Purpose : Calculates the correlation coefficient between observed and + ! estimated arrays. + ! + ! Input: + ! y - Observed data array. + ! yest - Estimated (modeled) data array. + ! + ! Output: + ! ccr - Correlation coefficient. + !------------------------------------------------------------ + real, intent(in) :: y(:) + real, intent(in) :: yest(:) + real :: ccr + real :: mean_y, mean_yest + real :: sum_y, sum_yest + real :: sum_num, sum_den_y, sum_den_yest + integer :: n + integer :: i + + n = size(y) + if (n /= size(yest)) then + print *, "Error: Arrays must have the same length" + ccr = 0.0 + return + endif + + ! Compute means + sum_y = sum(y) + sum_yest = sum(yest) + mean_y = sum_y / n + mean_yest = sum_yest / n + + ! Compute numerator and denominators for the correlation coefficient + sum_num = 0.0 + sum_den_y = 0.0 + sum_den_yest = 0.0 + do i = 1, n + sum_num = sum_num + (y(i) - mean_y) * (yest(i) - mean_yest) + sum_den_y = sum_den_y + (y(i) - mean_y)**2 + sum_den_yest = sum_den_yest + (yest(i) - mean_yest)**2 + end do + + if (sum_den_y == 0.0 .or. sum_den_yest == 0.0) then + print *, "Error: Zero variance in input arrays" + ccr = 0.0 + else + ccr = sum_num / sqrt(sum_den_y * sum_den_yest) + end if + + end function cal_ccr + + function cal_rms(k_obs, k_model, n) result(rms) + !------------------------------------------------------------ + ! Function: cal_rms + ! Purpose : Calculates the relative root mean square error between observed + ! and modeled K values. + ! + ! Input: + ! k_obs - Observed K values array. + ! k_model - Modeled K values array. + ! n - Number of elements. + ! + ! Output: + ! rms - Relative RMS error. + !------------------------------------------------------------ + implicit none + integer, intent(in) :: n + real, intent(in) :: k_obs(n), k_model(n) + real :: rms + real :: sum_sq_diff + integer :: i + + sum_sq_diff = 0.0 + + do i = 1, n + sum_sq_diff = sum_sq_diff + ((k_model(i) - k_obs(i)) / k_obs(i))**2 + end do + + rms = sqrt(sum_sq_diff / n) + end function cal_rms + + function cal_cdtm(y, yest) result(dtmc) + !------------------------------------------------------------ + ! Function: cal_cdtm + ! Purpose : Computes the coefficient of determination (R^2) between observed + ! and estimated data. + ! + ! Input: + ! y - Observed data array. + ! yest - Estimated data array. + ! + ! Output: + ! dtmc - Coefficient of determination (R^2). + !------------------------------------------------------------ + real, intent(in) :: y(:) + real, intent(in) :: yest(:) + real :: dtmc + real :: ss_tot, ss_res + real :: mean_y + integer :: n, i + + n = size(y) + if (n /= size(yest)) then + print *, "Error: Arrays must have the same length" + dtmc = 0.0 + return + endif + + mean_y = sum(y) / n + ss_tot = sum((y - mean_y)**2) + ss_res = sum((y - yest)**2) + + if (ss_tot == 0.0) then + print *, "Error: Zero total sum of squares" + dtmc = 0.0 + else + dtmc = 1.0 - (ss_res / ss_tot) + endif + + end function cal_cdtm + + function median(data) result(med) + !------------------------------------------------------------ + ! Function: median + ! Purpose : Computes the median of an array, ignoring values equal to -9999.0. + ! + ! Input: + ! data - Array of real numbers. + ! + ! Output: + ! med - Median value. + !------------------------------------------------------------ + implicit none + real, intent(in) :: data(:) + real :: med + real :: sorted_data(size(data)) + integer :: n_valid + integer :: i + + n_valid = 0 + do i = 1, size(data) + if (data(i) /= -9999.0) then + n_valid = n_valid + 1 + sorted_data(n_valid) = data(i) + end if + end do + + if (n_valid == 0) then + med = -9999.0 + return + end if + + call sort2(sorted_data(1:n_valid)) + + if (mod(n_valid, 2) == 0) then + med = (sorted_data(n_valid/2) + sorted_data(n_valid/2 + 1)) / 2.0 + else + med = sorted_data((n_valid + 1) / 2) + end if + + end function median + + subroutine sort2(arr) + !------------------------------------------------------------ + ! Subroutine: sort2 + ! Purpose : Sorts an array of real numbers in ascending order using + ! insertion sort. + ! + ! Input/Output: + ! arr - Array to be sorted. + !------------------------------------------------------------ + implicit none + real, intent(inout) :: arr(:) + integer :: i, j + real :: temp + + do i = 2, size(arr) + temp = arr(i) + j = i - 1 + do while (j >= 1 .and. arr(j) > temp) + arr(j + 1) = arr(j) + j = j - 1 + end do + arr(j + 1) = temp + end do + end subroutine sort2 + +!------------------------------------------------------------ +end module k_module \ No newline at end of file diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/process_lake_data.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/process_lake_data.py new file mode 100755 index 000000000..86633f55c --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/process_lake_data.py @@ -0,0 +1,173 @@ +import sys +import numpy as np +from netCDF4 import Dataset +#Main purpose: Processes lake data to be used in the river routing model. + +file_lat1m, file_lon1m, file_catmap, file_lake_mantag, file_lakecat_manfix = sys.argv[1:6] + +# Define constants +nlat = 10800 +nlon = 21600 + +# Read data from files +lats = np.loadtxt("temp/outlet_lat.txt", dtype=float) # Latitude of outlets +lons = np.loadtxt("temp/outlet_lon.txt", dtype=float) # Longitude of outlets +lat1m = np.loadtxt(file_lat1m, dtype=float) # Latitude grid +lon1m = np.loadtxt(file_lon1m, dtype=float) # Longitude grid + +# Function to find the nearest index in a coordinate array +def ind_nearest_coord(coord_array1, coord_array2): + """ + Find the index of the nearest value in coord_array2 for each value in coord_array1. + """ + indices = [] + for coord in coord_array1: + index = np.argmin(np.abs(coord_array2 - coord)) + indices.append(index) + return np.array(indices) + +# Find nearest indices for latitudes and longitudes +lati = ind_nearest_coord(lats, lat1m)+1 +loni = ind_nearest_coord(lons, lon1m)+1 + +#------------------------------------------------------------------------------------------------------ +ns = 3917 + +# Allocate array +catchind = np.zeros((nlat, nlon), dtype=int) + +# Read NetCDF file +def read_ncfile_int2d(filepath, varname, shape): + # Open the NetCDF file and read the specified variable + with Dataset(filepath, 'r') as nc: + data = nc.variables[varname][:].reshape(shape) + # Check for missing values and replace them with a default value + fill_value = nc.variables[varname]._FillValue if hasattr(nc.variables[varname], '_FillValue') else None + if fill_value is not None: + data = np.where(data == fill_value, -9999, data) # Replace missing values with 0 + return data + +catchind = read_ncfile_int2d(file_catmap, "CatchIndex", (nlat, nlon)) + +# Calculate catid +catid = np.zeros(ns, dtype=int) +for i in range(ns): + # Ensure indices are within bounds + if 0 < loni[i] <= nlon and 0 < lati[i] <= nlat: + catid[i] = catchind[lati[i] - 1, loni[i] - 1] # Adjust for 0-based indexing in Python + else: + catid[i] = -1 # Assign a default value for out-of-bounds indices + +#------------------------------------------------------------------------------------------------------ +# Constants +nall = 291284 +nv = 1782 +nv3 = 2097 + +# Read input data +aca_all = np.loadtxt("temp/Pfaf_acar.txt") + +# Initialize aca_model array +aca_model = np.full(ns, -9999.0) + +# Map aca_model using catid +for i in range(ns): + if catid[i] != -9999: + cid = catid[i] + aca_model[i] = aca_all[cid - 1] + +# Read observation data +aca_obs = np.loadtxt("temp/outlet_lakeacaOBS.txt") +outid_INCON = np.zeros(nv, dtype=int) + +# Filter inconsistent data +k = 0 +for i in range(ns): + if not (0.7 * aca_model[i] <= aca_obs[i] <= 1.3 * aca_model[i]): + outid_INCON[k] = i + 1 + k += 1 + +#print(k) + +#------------------------------------------------------------------------------------------------------ +# Read tags +tag_INCON = np.loadtxt(file_lake_mantag, dtype=int) + +# Update catid and aca_model based on tags +for i in range(nv): + oid = outid_INCON[i] + tag = tag_INCON[i] + if tag >= 1: + cid = tag + catid[oid - 1] = cid + aca_model[oid - 1] = aca_all[cid - 1] + else: + catid[oid - 1] = -9999 + aca_model[oid - 1] = -9999.0 + +# Compute flag_out +flag_out = np.where(aca_model != -9999, 1, 0) +#print(np.sum(flag_out)) + +#------------------------------------------------------------------------------------------------------ +# Read lakeid_out and compute absolute differences +lakeid_out = np.loadtxt("temp/outlet_lakeid.txt", dtype=int) +acaABSDIF_out = np.abs(aca_model - aca_obs) + +# Initialize collections +lakeid_collect = np.zeros(nv3, dtype=int) +outletid_collect = np.zeros(nv3, dtype=int) +acaABSDIF_collect = np.full(nv3, 1e10) +flag_2097_out = np.zeros(ns, dtype=int) +k = 0 + +# Collect valid outlets +for i in range(ns): + if flag_out[i] == 1: + lid = lakeid_out[i] + flag = 1 + if k >= 1: + for j in range(k): + if lid == lakeid_collect[j]: + flag = 0 + if acaABSDIF_out[i] < acaABSDIF_collect[j]: + flag_2097_out[outletid_collect[j]] = 0 + outletid_collect[j] = i + acaABSDIF_collect[j] = acaABSDIF_out[i] + flag_2097_out[i] = 1 + if flag == 1: + lakeid_collect[k] = lid + outletid_collect[k] = i + acaABSDIF_collect[k] = acaABSDIF_out[i] + flag_2097_out[i] = 1 + k += 1 + +#print(np.sum(flag_2097_out)) +np.savetxt("temp/lake_outlet_flag_valid_2097.txt", flag_2097_out, fmt="%d") + +#------------------------------------------------------------------------------------------------------ +# Update catid with valid flags +catid = np.where(flag_2097_out == 0, -9999, catid) + +# Collect valid outlet IDs +outidV = np.zeros(nv3, dtype=int) +k = 0 +for i in range(ns): + if flag_2097_out[i] == 1: + outidV[k] = i + k += 1 + +outidV += 1 + +# Fix multiple outlets in same catchment +catid_outfix_2097 = np.loadtxt(file_lakecat_manfix, dtype=int) +catid_outfix_out = np.full(ns, -9999, dtype=int) + +for i in range(nv3): + oid = outidV[i] + catid_outfix_out[oid - 1] = catid_outfix_2097[i] + +catid = np.where((catid_outfix_out != 0) & (catid_outfix_out != -9999), catid_outfix_out, catid) +np.savetxt("temp/lake_outlet_catid.txt", catid, fmt="%d") + + diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/read_input_TopoCat.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/read_input_TopoCat.f90 new file mode 100755 index 000000000..f6fc8897d --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/read_input_TopoCat.f90 @@ -0,0 +1,129 @@ +program main +!Main purpose: Reads lake and lake outlets information from Lake-TopoCat database. + +use constant, only : no, nvl, nvo, nl=>nl_lake + + implicit none + + ! Declare arrays for outlet and lake data: + integer, allocatable, dimension(:) :: lakeid_out, outid_out, lakeid_lake, lakeid_outV, outid_outV + real, allocatable, dimension(:) :: lat_out, lon_out, lat_outV, lon_outV, & + lakeaca_lake, lakearea_lake, lakeaca_outV, lakearea_outV + + ! Arrays for raw lake data: + real, allocatable, dimension(:) :: area_lake, aca_lake + integer, allocatable, dimension(:) :: id_lake + character(len=900) :: area_lake_file, id_lake_file, aca_lake_file + integer :: i, j, k + + character(len=900) :: file_lake_area + character(len=900) :: file_lake_id + character(len=900) :: file_lake_aca + character(len=900) :: file_lakeo_lakeid + character(len=900) :: file_lakeo_id + character(len=900) :: file_lakeo_lat + character(len=900) :: file_lakeo_lon + + if (command_argument_count() /= 7) then + print *, "no appropriate files found" + stop + endif + call get_command_argument(1, file_lake_area) + call get_command_argument(2, file_lake_id) + call get_command_argument(3, file_lake_aca) + call get_command_argument(4, file_lakeo_lakeid) + call get_command_argument(5, file_lakeo_id) + call get_command_argument(6, file_lakeo_lat) + call get_command_argument(7, file_lakeo_lon) + + ! Allocate arrays for raw lake data (size nl): + allocate(area_lake(nl), aca_lake(nl), id_lake(nl)) + + ! Initialize file names for input lake data: + area_lake_file = file_lake_area + id_lake_file = file_lake_id + aca_lake_file = file_lake_aca + + ! Read lake area, lake ID, and lake "aca" data from the input CSV files: + open(77, file=trim(area_lake_file), status="old") + read(77, *) area_lake + open(77, file=trim(id_lake_file), status="old") + read(77, *) id_lake + open(77, file=trim(aca_lake_file), status="old") + read(77, *) aca_lake + + ! Allocate arrays for filtered lake data (size nvl): + allocate(lakearea_lake(nvl)) + allocate(lakeid_lake(nvl)) + allocate(lakeaca_lake(nvl)) + + k = 0 + ! Filter lakes: select only those with an area greater than or equal to 50. + do i = 1, nl + if (area_lake(i) .ge. 50.0) then + k = k + 1 + lakearea_lake(k) = area_lake(i) + lakeid_lake(k) = id_lake(i) + lakeaca_lake(k) = aca_lake(i) + end if + end do +!------------------------------------------------------------------------------------- + + ! Allocate arrays for outlet data (raw arrays with size 'no'): + allocate(lakeid_out(no), outid_out(no), lat_out(no), lon_out(no)) + ! Allocate arrays for matched outlet data (size 'nvo'): + allocate(lakeid_outV(nvo), outid_outV(nvo), lat_outV(nvo), lon_outV(nvo), lakeaca_outV(nvo), lakearea_outV(nvo)) + + ! Read outlet data from CSV files: + open(77, file=trim(file_lakeo_lakeid)) + read(77, *) lakeid_out + open(77, file=trim(file_lakeo_id)) + read(77, *) outid_out + open(77, file=trim(file_lakeo_lat)) + read(77, *) lat_out + open(77, file=trim(file_lakeo_lon)) + read(77, *) lon_out + + ! Match outlet records to filtered lakes: + k = 0 + do i = 1, no + do j = 1, nvl + if (lakeid_out(i) == lakeid_lake(j)) then + k = k + 1 + outid_outV(k) = outid_out(i) + lat_outV(k) = lat_out(i) + lon_outV(k) = lon_out(i) + lakeid_outV(k) = lakeid_lake(j) + lakeaca_outV(k) = lakeaca_lake(j) + lakearea_outV(k) = lakearea_lake(j) + end if + end do + end do + + ! Write matched outlet data to output text files: + open(88, file="temp/outlet_lat.txt") + do i = 1, nvo + write(88, *) lat_outV(i) + end do + + open(88, file="temp/outlet_lon.txt") + do i = 1, nvo + write(88, *) lon_outV(i) + end do + + open(88, file="temp/outlet_lakeid.txt") + do i = 1, nvo + write(88, *) lakeid_outV(i) + end do + + open(88, file="temp/outlet_lakeacaOBS.txt") + do i = 1, nvo + write(88, *) lakeaca_outV(i) + end do + + open(88, file="temp/lake_outlet_lakearea.txt") + do i = 1, nvo + write(88, *) lakearea_outV(i) + end do + +end \ No newline at end of file diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/readme.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/readme.txt new file mode 100644 index 000000000..d6144964c --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/readme.txt @@ -0,0 +1,76 @@ +v1 12/05/2025, Yujin Zeng + +The "preproc/routing_model" package is used for creating input data related to the GEOS routing model that would be used in the make_bcs. + +Usage: + • On NCCS/Discover: + 1. source g5_modules.sh + 2. build GEOSldas code + 3. python3 run_routing_preproc.py under the build/src directory + • Off Discover: contact yujin.zeng@nasa.gov + +The tasks completed by each F90 or Python program are briefly described as follows: + +1. get_Pfaf_file.f90 + Reads the Pfafstetter code dataset and generates + files for the connectivity of catchments in the routing network. + +2. get_latloni_cellarea.py + Computes grid-cell index arrays and per-cell areas for 1-m high-res grid. + +3. get_num_sub_catchment.f90 + Parses high-res map of catchment index to get the area and + coordinates (of model grid) of each sub-catchments within each main catchment. + +4. get_lonlat_bond.f90 + Extracts the latitude/longitude boundaries of each catchment-tile from + catchment definition files. + +5. get_lonlati_maptile.py + Assigns a catchment‐tile index from catchment definition files to each model grid cell. + +6. get_isub.f90 + Assigns a catchment‐tile index from maptile files to each sub-catchment. + +7. get_area.f90 + Gets the area for each catchment-tile. + +8. get_Qr_clmt.f90 + Reads SMAP L4 runoff data (2016–2023) from a NetCDF file and computes the climatological + mean discharge for each catchment. + +9. get_river_length.f90 + Determines main river channel lengths for each catchment by using HydroSHEDS + data of distance to sink. + +10. get_K_model_calik.f90 + Calculates the K parameter used in the river routing model. + +11. get_dam_data.py + Processes reservoir (dam) data: reads dam locations and usage information from GRanD database. + +12. read_input_TopoCat.f90 + Reads lake and lake outlets information from Lake-TopoCa database. + +13. process_lake_data.py + Processes lake data to be used in the river routing model. + +14. create_river_input.py + Combines all inputs created above to the one nc file river_input.nc + +The explanations for the input files of this package can be found in the input directory. + +The outputs from this package (cellarea.nc and river_input.nc located in the created output folder) will be used as input to the make_bcs. They are listed as follows: + + cellarea.nc: cell area [m^2] of the 1-min grid. + + river_input.nc: River parameters on the catchment space. More detail can be seen in the attributes of each variable in the netcdf file. + + + + + + + + + diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/river_read.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/river_read.f90 new file mode 100755 index 000000000..1ed66ce0b --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/river_read.f90 @@ -0,0 +1,193 @@ +module river_read +!module for reading river routing-related netcdf data + + implicit none + include 'netcdf.inc' + + public :: read_ncfile_int1d + public :: read_ncfile_real1d + public :: read_ncfile_double1d + + public :: read_ncfile_int2d + public :: read_ncfile_int3d + public :: read_ncfile_real2d + public :: read_ncfile_real3d + public :: read_ncfile_double2d + public :: read_ncfile_double3d + + contains +!------------------------------------------------------------------------------------------ + subroutine read_ncfile_int1d(filename,varname,var,n) + character(len=*), intent(in) :: filename + character(len=*), intent(in) :: varname + integer, intent(in) :: n + integer, intent(inout) :: var(n) + + character(len=4) :: subname="read" + integer :: ncid, varid + + call check_ret(nf_open(filename,0,ncid),subname) + call check_ret(nf_inq_varid(ncid,varname,varid),subname) + call check_ret(nf_get_var_int(ncid,varid,var),subname) + call check_ret(nf_close(ncid), subname) + + end subroutine read_ncfile_int1d +!------------------------------------------------------------------------------------------ + subroutine read_ncfile_real1d(filename,varname,var,n) + character(len=*), intent(in) :: filename + character(len=*), intent(in) :: varname + integer, intent(in) :: n + real, intent(inout) :: var(n) + + character(len=4) :: subname="read" + integer :: ncid, varid + + call check_ret(nf_open(filename,0,ncid),subname) + call check_ret(nf_inq_varid(ncid,varname,varid),subname) + call check_ret(nf_get_var_real(ncid,varid,var),subname) + call check_ret(nf_close(ncid), subname) + + end subroutine read_ncfile_real1d +!------------------------------------------------------------------------------------------ + subroutine read_ncfile_double1d(filename,varname,var,n) + character(len=*), intent(in) :: filename + character(len=*), intent(in) :: varname + integer, intent(in) :: n + real*8, intent(inout) :: var(n) + + character(len=4) :: subname="read" + integer :: ncid, varid + + call check_ret(nf_open(filename,0,ncid),subname) + call check_ret(nf_inq_varid(ncid,varname,varid),subname) + call check_ret(nf_get_var_double(ncid,varid,var),subname) + call check_ret(nf_close(ncid), subname) + + end subroutine read_ncfile_double1d +!------------------------------------------------------------------------------------------ + subroutine read_ncfile_int2d(filename,varname,var,nlon,nlat) + character(len=*), intent(in) :: filename + character(len=*), intent(in) :: varname + integer, intent(in) :: nlon, nlat + integer, intent(inout) :: var(nlon,nlat) + + character(len=4) :: subname="read" + integer :: ncid, varid + + call check_ret(nf_open(filename,0,ncid),subname) + call check_ret(nf_inq_varid(ncid,varname,varid),subname) + call check_ret(nf_get_var_int(ncid,varid,var),subname) + call check_ret(nf_close(ncid), subname) + + end subroutine read_ncfile_int2d +!------------------------------------------------------------------------------------------ + subroutine read_ncfile_int3d(filename,varname,var,nlon,nlat,nlev) + character(len=*), intent(in) :: filename + character(len=*), intent(in) :: varname + integer, intent(in) :: nlon, nlat, nlev + integer, intent(inout) :: var(nlon,nlat,nlev) + + character(len=4) :: subname="read" + integer :: ncid, varid + + call check_ret(nf_open(filename,0,ncid),subname) + call check_ret(nf_inq_varid(ncid,varname,varid),subname) + call check_ret(nf_get_var_int(ncid,varid,var),subname) + call check_ret(nf_close(ncid), subname) + + end subroutine read_ncfile_int3d +!------------------------------------------------------------------------------------------ + subroutine read_ncfile_real2d(filename,varname,var,nlon,nlat) + character(len=*), intent(in) :: filename + character(len=*), intent(in) :: varname + integer, intent(in) :: nlon, nlat + real, intent(inout) :: var(nlon,nlat) + + character(len=4) :: subname="read" + integer :: ncid, varid + + call check_ret(nf_open(filename,0,ncid),subname) + call check_ret(nf_inq_varid(ncid,varname,varid),subname) + call check_ret(nf_get_var_real(ncid,varid,var),subname) + call check_ret(nf_close(ncid), subname) + + end subroutine read_ncfile_real2d +!------------------------------------------------------------------------------------------ + subroutine read_ncfile_real3d(filename,varname,var,nlon,nlat,nlev) + character(len=*), intent(in) :: filename + character(len=*), intent(in) :: varname + integer, intent(in) :: nlon, nlat, nlev + real, intent(inout) :: var(nlon,nlat,nlev) + + character(len=4) :: subname="read" + integer :: ncid, varid + + call check_ret(nf_open(filename,0,ncid),subname) + call check_ret(nf_inq_varid(ncid,varname,varid),subname) + call check_ret(nf_get_var_real(ncid,varid,var),subname) + call check_ret(nf_close(ncid), subname) + + end subroutine read_ncfile_real3d +!------------------------------------------------------------------------------------------ + subroutine read_ncfile_double2d(filename,varname,var,nlon,nlat) + character(len=*), intent(in) :: filename + character(len=*), intent(in) :: varname + integer, intent(in) :: nlon, nlat + real*8, intent(inout) :: var(nlon,nlat) + + character(len=4) :: subname="read" + integer :: ncid, varid + + call check_ret(nf_open(filename,0,ncid),subname) + call check_ret(nf_inq_varid(ncid,varname,varid),subname) + call check_ret(nf_get_var_double(ncid,varid,var),subname) + call check_ret(nf_close(ncid), subname) + + end subroutine read_ncfile_double2d +!------------------------------------------------------------------------------------------ + subroutine read_ncfile_double3d(filename,varname,var,nlon,nlat,nlev) + character(len=*), intent(in) :: filename + character(len=*), intent(in) :: varname + integer, intent(in) :: nlon, nlat, nlev + real*8, intent(inout) :: var(nlon,nlat,nlev) + + character(len=4) :: subname="read" + integer :: ncid, varid + + call check_ret(nf_open(filename,0,ncid),subname) + call check_ret(nf_inq_varid(ncid,varname,varid),subname) + call check_ret(nf_get_var_double(ncid,varid,var),subname) + call check_ret(nf_close(ncid), subname) + + end subroutine read_ncfile_double3d +!------------------------------------------------------------------------------------------ + subroutine check_ret(ret, calling) + integer, intent(in) :: ret + character(len=*) :: calling + + if (ret /= NF_NOERR) then + write(6,*)'netcdf error from ',trim(calling) + call endrun(nf_strerror(ret)) + end if + end subroutine check_ret +!----------------------------------------------------------------------- + subroutine endrun(msg,subname) + character(len=*), intent(in), optional :: msg + character(len=*), intent(in), optional :: subname + + if (present (subname)) then + write(6,*) 'ERROR in subroutine :', trim(subname) + end if + + if (present (msg)) then + write(6,*)'ENDRUN:', msg + else + write(6,*) 'ENDRUN: called without a message string' + end if + + stop + end subroutine endrun +!----------------------------------------------------------------------- + +end module river_read + diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/run_routing_preproc.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/run_routing_preproc.py new file mode 100644 index 000000000..e42ec8a48 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing_model/run_routing_preproc.py @@ -0,0 +1,160 @@ +#!/usr/bin/env python3 + +#source g5_modules before run to get the necessary env + +import os +import shutil +import subprocess +from pathlib import Path + +def run(cmd): + """Run a command and exit on failure.""" + print(f">>> {' '.join(cmd)}") + subprocess.run(cmd, check=True) + +def main(): + # ----------------------------- + # Define input and output paths + # ----------------------------- + +# Path to "bcs_shared" directory: + file_path1 = "/discover/nobackup/projects/gmao/bcs_shared" + file_path2 = "/discover/nobackup/projects/gmao/bcs_shared/preprocessing_bcs_inputs/land/routing/routing_model" + + + file_pfafrout = file_path1 + "/make_bcs_inputs/land/topo/v1/SRTM-TopoData/Pfafcatch-routing.dat" + file_pfafmap = file_path1 + "/make_bcs_inputs/land/topo/v1/SRTM-TopoData/SRTM_PfafData.nc" + file_catdef = file_path1 + "/fvInput/ExtData/esm/tiles/v14/land/EASEv2_M09/clsm/catchment.def" + + file_lat1m = file_path2 + "/lat_1m.txt" + file_lon1m = file_path2 + "/lon_1m.txt" + file_lat = file_path2 + "/lat_SMAPL4.txt" + file_lon = file_path2 + "/lon_SMAPL4.txt" + file_clmtrunf = file_path2 + "/SMAPL4_OL7000_runoff_mean_2016_2023.nc" + file_ldn = file_path2 + "/hyd_glo_ldn_15s.nc" + file_hyelev = file_path2 + "/hyd_glo_dem_15s.nc" + file_vel = file_path2 + "/velocity.txt" + file_dis = file_path2 + "/discharge.txt" + file_usid = file_path2 + "/USGSID.txt" + file_lats = file_path2 + "/lat_for_site.txt" + file_lons = file_path2 + "/lon_for_site.txt" + file_gage_id = file_path2 + "/id_gagesii.txt" + file_gage_acar = file_path2 + "/acar_gagesii.txt" + + file_latdam = file_path2 + "/lat_dam_grand.txt" + file_londam = file_path2 + "/lon_dam_grand.txt" + file_acadam = file_path2 + "/catch_acar_grand.txt" + file_damcat_manfix = file_path2 + "/catid_dam_manfix.txt" + file_dam_manflag = file_path2 + "/flag_dam_manfix.txt" + file_dam_use = file_path2 + "/main_use_grand.txt" + file_damflood = file_path2 + "/flood_use_grand.txt" + file_damarea = file_path2 + "/area_skm_grand.txt" + file_damcap = file_path2 + "/cap_max_grand.txt" + + file_lake_area = file_path2 + "/Lake_area.csv" + file_lake_id = file_path2 + "/Hylak_id_lake.csv" + file_lake_aca = file_path2 + "/acar_lake.csv" + file_lakeo_lakeid = file_path2 + "/Hylak_id_lakeout.csv" + file_lakeo_id = file_path2 + "/Outlet_id_lakeout.csv" + file_lakeo_lat = file_path2 + "/Outlet_lat_lakeout.csv" + file_lakeo_lon = file_path2 + "/Outlet_lon_lakeout.csv" + file_lake_mantag = file_path2 + "/catid_lake_manfix.txt" + file_lakecat_manfix = file_path2 + "/catid_lake_multout_manfix.txt" + + # ----------------------------- + # Ensure output and temp directory exists + # ----------------------------- + subprocess.run(["mkdir", "-p", "output"], check=True) + subprocess.run(["mkdir", "-p", "temp"], check=True) + + # ----------------------------- + # Copy dam area and capacity files + # ----------------------------- + shutil.copy(file_damarea, "temp/") + shutil.copy(file_damcap, "temp/") + + # ----------------------------- + # River processing section + # ----------------------------- + # Compile and run Pfafcatch routing generator + run(["./get_Pfaf_file.x", file_pfafrout]) + + # Generate latitude/longitude indices and cell areas + run([ + "python3", "get_latloni_cellarea.py", + file_lat, file_lon, + file_lat1m, file_lon1m, + ]) + + # Compute number of sub-catchments + run([f"./get_num_sub_catchment.x", file_pfafmap]) + + # longitude-latitude boundary + run([f"./get_lonlat_bond.x", file_catdef]) + + # Map tile longitude/latitude + run(["python3", f"get_lonlati_maptile.py", file_lat, file_lon]) + # isub calculators + run([f"./get_isub.x"]) + + # Calculate area of each catchment + run([f"./get_area.x", file_pfafmap]) + + # Compute climatological runoff + run(["./get_Qr_clmt.x", file_clmtrunf]) + + # Determine river lengths + run([ + "./get_river_length.x", + file_pfafmap, file_ldn, + file_hyelev, file_pfafrout + ]) + + # Calibrate K model using velocity and discharge data + run([ + "./get_K_model_calik.x", + file_vel, file_dis, file_usid, + file_lats, file_lons, + file_lat1m, file_lon1m, + file_pfafmap, + file_gage_id, file_gage_acar + ]) + + # ----------------------------- + # Reservoir (dam) processing + # ----------------------------- + run([ + "python3", "get_dam_data.py", + file_latdam, file_londam, + file_lat1m, file_lon1m, + file_pfafmap, file_acadam, + file_damcat_manfix, file_dam_manflag, + file_dam_use, file_damflood + ]) + + # ----------------------------- + # Lake processing section + # ----------------------------- + run([ + "./read_input_TopoCat.x", + file_lake_area, file_lake_id, + file_lake_aca, + file_lakeo_lakeid, file_lakeo_id, + file_lakeo_lat, file_lakeo_lon + ]) + + run([ + "python3", "process_lake_data.py", + file_lat1m, file_lon1m, + file_pfafmap, + file_lake_mantag, file_lakecat_manfix + ]) + + run([ + "python3", "create_river_input.py", + ]) + + subprocess.run(["rm", "-rf", "temp"], check=True) + +if __name__ == "__main__": + main()