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..59910cdd6 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,7 @@ 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 + use GEOS_RouteGridCompMod, only : RouteSetServices => SetServices implicit none private @@ -195,19 +195,21 @@ 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 + allocate (ROUTE(NUM_CATCH), stat=status) + VERIFY_(STATUS) + 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) @@ -1453,16 +1455,16 @@ 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 + 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 ( & @@ -1486,16 +1488,16 @@ 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 + 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 @@ -1669,6 +1671,7 @@ subroutine Run1(GC, IMPORT, EXPORT, CLOCK, RC ) !-------------------------------- DO I = 1, size(GCS) + if (trim(GCnames(i)) == "ROUTE") cycle call MAPL_TimerOn(MAPL,trim(GCnames(i)), RC=STATUS ); VERIFY_(STATUS) call ESMF_GridCompRun(GCS(I), importState=GIM(I), exportState=GEX(I), & CLOCK=CLOCK, PHASE=1, userRC=STATUS) 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 index 8a502e3e7..7629ae8f1 100644 --- 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 @@ -1,10 +1,10 @@ esma_set_this () set (srcs - #GEOS_RouteGridComp.F90 + GEOS_RouteGridComp.F90 routing_model.F90 ) -esma_add_library (${this} SRCS ${srcs} DEPENDENCIES MAPL ESMF::ESMF NetCDF::NetCDF_Fortran) +esma_add_library (${this} SRCS ${srcs} DEPENDENCIES MAPL 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 index ad2be4db2..f7628ee38 100644 --- 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 @@ -18,8 +18,6 @@ module GEOS_RouteGridCompMod ! 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: @@ -27,14 +25,19 @@ module GEOS_RouteGridCompMod use MAPL_Mod use MAPL_ConstantsMod use ROUTING_MODEL, ONLY: & - river_routing, ROUTE_DT + river_routing_lin, river_routing_hyd, ROUTE_DT #if 0 USE catch_constants, ONLY: & N_CatG => N_Pfaf_Catchs #endif + use, intrinsic :: iso_c_binding implicit none integer, parameter :: N_CatG = 291284 + integer,parameter :: upmax=34 + character(len=500) :: inputdir="/discover/nobackup/yzeng3/data/river_input/" + integer,save :: nmax + private type T_RROUTE_STATE @@ -42,15 +45,57 @@ module GEOS_RouteGridCompMod type (ESMF_RouteHandle) :: routeHandle type (ESMF_Field) :: field integer :: nTiles + integer :: nt_global + integer :: nt_local integer :: comm integer :: nDes integer :: myPe integer :: minCatch integer :: maxCatch integer, pointer :: pfaf(:) => NULL() - real, pointer :: tile_area(:) => NULL() + real, pointer :: tile_area(:) => NULL() !m2 + integer, pointer :: nsub(:) => NULL() + integer, pointer :: subi(:,:) => NULL() + real, pointer :: subarea(:,:) => NULL() !m2 + + integer, pointer :: scounts_global(:) => NULL() + integer, pointer :: rdispls_global(:) => NULL() + integer, pointer :: scounts_cat(:) => NULL() + integer, pointer :: rdispls_cat(:) => NULL() + + real, pointer :: runoff_save(:) => NULL() + real, pointer :: areacat(:) => NULL() !m2 + real, pointer :: lengsc(:) => NULL() !m + + real, pointer :: wstream(:) => NULL() !m3 + real, pointer :: wriver(:) => NULL() !m3 + integer, pointer :: downid(:) => NULL() + integer, pointer :: upid(:,:) => NULL() + + real, pointer :: wriver_acc(:) => NULL() + real, pointer :: wstream_acc(:) => NULL() + real, pointer :: qoutflow_acc(:) => NULL() + real, pointer :: qsflow_acc(:) => NULL() + + real, pointer :: lstr(:) => NULL() !m + real, pointer :: qri_clmt(:) => NULL() !m3/s + real, pointer :: qin_clmt(:) => NULL() !m3/s + real, pointer :: qstr_clmt(:) =>NULL() !m3/s + real, pointer :: K(:) => NULL() + real, pointer :: Kstr(:) => NULL() + end type T_RROUTE_STATE + + interface + function mkdir(path,mode) bind(c,name="mkdir") + use iso_c_binding + integer(c_int) :: mkdir + character(kind=c_char,len=1) :: path(*) + integer(c_int16_t), value :: mode + end function mkdir + end interface + ! Wrapper for extracting internal state ! ------------------------------------- type RROUTE_WRAP @@ -132,7 +177,11 @@ subroutine SetServices ( GC, RC ) call MAPL_GridCompSetEntryPoint ( GC, ESMF_METHOD_INITIALIZE, Initialize, RC=STATUS ) VERIFY_(STATUS) - call MAPL_GridCompSetEntryPoint (GC, ESMF_METHOD_RUN, Run, RC=STATUS) +! call MAPL_GridCompSetEntryPoint (GC, ESMF_METHOD_RUN, Run, RC=STATUS) +! VERIFY_(STATUS) + call MAPL_GridCompSetEntryPoint ( GC, ESMF_METHOD_RUN, RUN1, RC=STATUS ) + VERIFY_(STATUS) + call MAPL_GridCompSetEntryPoint ( GC, ESMF_METHOD_RUN, RUN2, RC=STATUS ) VERIFY_(STATUS) !------------------------------------------------------------ @@ -186,104 +235,8 @@ subroutine SetServices ( GC, RC ) 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) @@ -307,7 +260,9 @@ subroutine SetServices ( GC, RC ) call MAPL_TimerAdd(GC, name="INITIALIZE" ,RC=STATUS) VERIFY_(STATUS) - call MAPL_TimerAdd(GC, name="RUN" ,RC=STATUS) + call MAPL_TimerAdd(GC, name="RUN1" ,RC=STATUS) + VERIFY_(STATUS) + call MAPL_TimerAdd(GC, name="RUN2" ,RC=STATUS) VERIFY_(STATUS) ! All done @@ -367,22 +322,50 @@ subroutine INITIALIZE (GC,IMPORT, EXPORT, CLOCK, RC ) integer, pointer :: ims(:) => NULL() integer, pointer :: pfaf(:) => NULL() integer, pointer :: arbSeq(:) => NULL() + integer, pointer :: arbSeq_pf(:) => NULL() + integer, pointer :: arbSeq_ori(:) => NULL() integer, allocatable :: arbIndex(:,:) real, pointer :: tile_area_src(:) => NULL() - real, pointer :: tile_area(:) => NULL() + integer,pointer :: local_id(:) => NULL() + real, pointer :: tile_area_local(:) => NULL(), tile_area_global(:) => NULL() + real, pointer :: tile_area(:) => NULL() real, pointer :: ptr2(:) => NULL() + + real,pointer :: subarea_global(:,:)=> NULL(),subarea(:,:)=> NULL() ! Arrays for sub-area and fractions + integer,pointer :: subi_global(:,:)=> NULL(),subi(:,:)=> NULL() + integer,pointer :: nsub_global(:)=> NULL(),nsub(:)=> NULL() + real,pointer :: area_cat_global(:)=> NULL(),area_cat(:)=> NULL() + integer,pointer :: scounts(:)=>NULL() + integer,pointer :: scounts_global(:)=>NULL(),rdispls_global(:)=>NULL() + integer,pointer :: scounts_cat(:)=>NULL(),rdispls_cat(:)=>NULL() + + real,pointer :: runoff_save(:)=>NULL(), areacat(:)=>NULL() + real,pointer :: lengsc_global(:)=>NULL(), lengsc(:)=>NULL(), buff_global(:)=>NULL() + integer,pointer :: downid_global(:)=>NULL(), downid(:)=>NULL() + integer,pointer :: upid_global(:,:)=>NULL(), upid(:,:)=>NULL() + + real,pointer :: wstream(:)=>NULL(),wriver(:)=>NULL() + real,pointer :: wstream_global(:)=>NULL(),wriver_global(:)=>NULL() type (T_RROUTE_STATE), pointer :: route => null() type (RROUTE_wrap) :: wrap + type(ESMF_Time) :: CurrentTime + integer :: YY,MM,DD,HH,MMM,SS + character(len=4) :: yr_s + character(len=2) :: mon_s,day_s + character(len=3) :: resname + + real, pointer :: dataPtr(:) + integer :: j,nt_local,mpierr,it ! ------------------ ! begin + call ESMF_UserCompGetInternalState ( GC, 'RiverRoute_state',wrap,status ) VERIFY_(STATUS) route => wrap%ptr - ! get vm ! extract comm call ESMF_VMGetCurrent(VM, RC=STATUS) @@ -398,152 +381,268 @@ subroutine INITIALIZE (GC,IMPORT, EXPORT, CLOCK, RC ) route%comm = comm route%ndes = ndes route%mype = mype - + + allocate(ims(1:ndes)) ! 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., & + VERIFY_(STATUS) + ! extract Pfaf (TILEI on the "other" grid) + call MAPL_LocStreamGet(locstream, & tileGrid=tilegrid, nt_global=nt_global, RC=status) - VERIFY_(STATUS) - + VERIFY_(STATUS) + route%nt_global = nt_global + + if(nt_global==112573)then + resname="M36" + nmax=150 + else if(nt_global==1684725)then + resname="M09" + nmax=458 + else + if(mapl_am_I_root())then + print *,"unknown grid for routing model" + stop + endif + endif ! 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) + call MAPL_LocStreamGet(locstream, TILEAREA = tile_area_src, LOCAL_ID=local_id, RC=status) + VERIFY_(STATUS) + nt_local=size(tile_area_src,1) + route%nt_local=nt_local + ntiles = maxCatch-minCatch+1 + allocate(arbSeq_pf(maxCatch-minCatch+1)) + arbSeq_pf = [(i, i = minCatch, maxCatch)] + ! redist pfaf (NOTE: me might need a second routehandle for integers) - route%pfaf => arbSeq - route%ntiles = ntiles + route%pfaf => arbSeq_pf + 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) + route%maxCatch = maxCatch + ! Read sub-area data from text files + allocate(nsub_global(N_CatG),subarea_global(nmax,N_CatG)) + open(77,file=trim(inputdir)//"/Pfaf_nsub_"//trim(resname)//".txt",status="old",action="read"); read(77,*)nsub_global; close(77) + open(77,file=trim(inputdir)//"/Pfaf_asub_"//trim(resname)//".txt",status="old",action="read"); read(77,*)subarea_global; close(77) + allocate(nsub(ntiles),subarea(nmax,ntiles)) + nsub=nsub_global(minCatch:maxCatch) + subarea=subarea_global(:,minCatch:maxCatch) + subarea=subarea*1.e6 !km2->m2 + deallocate(nsub_global,subarea_global) + + route%nsub => nsub + route%subarea => subarea + allocate(subi_global(nmax,N_CatG),subi(nmax,ntiles)) + open(77,file=trim(inputdir)//"/Pfaf_isub_"//trim(resname)//".txt",status="old",action="read");read(77,*)subi_global;close(77) + subi=subi_global(:,minCatch:maxCatch) + route%subi => subi + deallocate(subi_global) + + + allocate(scounts(ndes),scounts_global(ndes),rdispls_global(ndes)) + scounts=0 + scounts(mype+1)=nt_local + call MPI_Allgather(scounts(mype+1), 1, MPI_INTEGER, scounts_global, 1, MPI_INTEGER, MPI_COMM_WORLD, mpierr) + rdispls_global(1)=0 + do i=2,nDes + rdispls_global(i)=rdispls_global(i-1)+scounts_global(i-1) + enddo + deallocate(scounts) + route%scounts_global=>scounts_global + route%rdispls_global=>rdispls_global + + allocate(scounts(ndes),scounts_cat(ndes),rdispls_cat(ndes)) + scounts=0 + scounts(mype+1)=ntiles + call MPI_Allgather(scounts(mype+1), 1, MPI_INTEGER, scounts_cat, 1, MPI_INTEGER, MPI_COMM_WORLD, mpierr) + rdispls_cat(1)=0 + do i=2,nDes + rdispls_cat(i)=rdispls_cat(i-1)+scounts_cat(i-1) + enddo + deallocate(scounts) + route%scounts_cat=>scounts_cat + route%rdispls_cat=>rdispls_cat + + allocate(runoff_save(1:nt_local)) + route%runoff_save => runoff_save + route%runoff_save=0. + + allocate(tile_area_local(nt_local),tile_area_global(nt_global)) + open(77,file=trim(inputdir)//"/area_"//trim(resname)//"_1d.txt",status="old",action="read");read(77,*)tile_area_global;close(77) + tile_area_local=tile_area_global(rdispls_global(mype+1)+1:rdispls_global(mype+1)+nt_local)*1.e6 !km2->m2 + route%tile_area => tile_area_local + deallocate(tile_area_global) + + allocate(areacat(1:ntiles)) + areacat=0. + do i=1,ntiles + do j=1,nmax + it=route%subi(j,i) + if(it>0)then + areacat(i)=areacat(i)+route%subarea(j,i) + endif + if(it==0)exit + enddo + enddo + route%areacat=>areacat + + allocate(lengsc_global(n_catg),lengsc(ntiles)) + open(77,file=trim(inputdir)//"/Pfaf_lriv_PR.txt",status="old",action="read");read(77,*)lengsc_global;close(77) + lengsc=lengsc_global(minCatch:maxCatch)*1.e3 !km->m + route%lengsc=>lengsc + deallocate(lengsc_global) + + allocate(downid_global(n_catg),downid(ntiles)) + open(77,file=trim(inputdir)//"/downstream_1D_new_noadj.txt",status="old",action="read");read(77,*)downid_global;close(77) + downid=downid_global(minCatch:maxCatch) + route%downid=>downid + deallocate(downid_global) + + allocate(upid_global(upmax,n_catg),upid(upmax,ntiles)) + open(77,file=trim(inputdir)//"/upstream_1D.txt",status="old",action="read");read(77,*)upid_global;close(77) + upid=upid_global(:,minCatch:maxCatch) + route%upid=>upid + deallocate(upid_global) + + 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) + write(yr_s,'(I4.4)')YY + write(mon_s,'(I2.2)')MM + write(day_s,'(I2.2)')DD + if(mapl_am_I_root())print *, "init time is ", YY, "/", MM, "/", DD, " ", HH, ":", MMM, ":", SS + allocate(wriver(ntiles),wstream(ntiles)) + allocate(wriver_global(n_catg),wstream_global(n_catg)) + open(77,file="../input/restart/river_storage_rs_"//trim(yr_s)//trim(mon_s)//trim(day_s)//".txt",status="old",action="read",iostat=status) + if(status==0)then + read(77,*)wriver_global;close(77) + else + close(77) + open(78,file=trim(inputdir)//"/river_storage_rs_"//trim(yr_s)//trim(mon_s)//trim(day_s)//".txt",status="old",action="read",iostat=status) + if(status==0)then + read(78,*)wriver_global;close(78) + else + close(78) + open(79,file=trim(inputdir)//"/river_storage_rs.txt",status="old",action="read",iostat=status) + if(status==0)then + read(79,*)wriver_global;close(79) + else + close(79) + wriver_global=0. + endif + endif + endif + open(77,file="../input/restart/stream_storage_rs_"//trim(yr_s)//trim(mon_s)//trim(day_s)//".txt",status="old",action="read",iostat=status) + if(status==0)then + read(77,*)wstream_global;close(77) + else + close(77) + open(78,file=trim(inputdir)//"/stream_storage_rs_"//trim(yr_s)//trim(mon_s)//trim(day_s)//".txt",status="old",action="read",iostat=status) + if(status==0)then + read(78,*)wstream_global;close(78) + else + close(78) + open(79,file=trim(inputdir)//"/stream_storage_rs.txt",status="old",action="read",iostat=status) + if(status==0)then + read(79,*)wstream_global;close(79) + else + close(79) + wstream_global=0. + endif + endif + endif + if(mapl_am_I_root())print *, "init river storage is: ",sum(wriver_global)/1.e9 + if(mapl_am_I_root())print *, "init stream storage is: ",sum(wstream_global)/1.e9 + wriver=wriver_global(minCatch:maxCatch) + wstream=wstream_global(minCatch:maxCatch) + deallocate(wriver_global,wstream_global) + route%wstream=>wstream + route%wriver=>wriver + + allocate(route%wriver_acc(ntiles),route%wstream_acc(ntiles),route%qoutflow_acc(ntiles),route%qsflow_acc(ntiles)) + route%wriver_acc=0. + route%wstream_acc=0. + route%qoutflow_acc=0. + route%qsflow_acc=0. + + !input for geometry hydraulic + allocate(buff_global(n_catg),route%lstr(ntiles)) + open(77,file=trim(inputdir)//"/Pfaf_lstr_PR.txt",status="old",action="read");read(77,*)buff_global;close(77) + route%lstr=buff_global(minCatch:maxCatch)*1.e3 !km->m + deallocate(buff_global) + + allocate(buff_global(n_catg),route%K(ntiles)) + open(77,file=trim(inputdir)//"/Pfaf_Kv_PR_0p35_0p45_0p2_n0p2.txt",status="old",action="read");read(77,*)buff_global;close(77) + route%K=buff_global(minCatch:maxCatch) + deallocate(buff_global) + + allocate(buff_global(n_catg),route%Kstr(ntiles)) + open(77,file=trim(inputdir)//"/Pfaf_Kstr_PR_fac1_0p35_0p45_0p2_n0p2.txt",status="old",action="read");read(77,*)buff_global;close(77) + route%Kstr=buff_global(minCatch:maxCatch) + deallocate(buff_global) + + allocate(buff_global(n_catg),route%qri_clmt(ntiles)) + open(77,file=trim(inputdir)//"/Pfaf_qri.txt",status="old",action="read");read(77,*)buff_global;close(77) + route%qri_clmt=buff_global(minCatch:maxCatch) !m3/s + deallocate(buff_global) + + allocate(buff_global(n_catg),route%qin_clmt(ntiles)) + open(77,file=trim(inputdir)//"/Pfaf_qin.txt",status="old",action="read");read(77,*)buff_global;close(77) + route%qin_clmt=buff_global(minCatch:maxCatch) !m3/s + deallocate(buff_global) + + allocate(buff_global(n_catg),route%qstr_clmt(ntiles)) + open(77,file=trim(inputdir)//"/Pfaf_qstr.txt",status="old",action="read");read(77,*)buff_global;close(77) + route%qstr_clmt=buff_global(minCatch:maxCatch) !m3/s + deallocate(buff_global) + + !if (mapl_am_I_root())then + ! open(88,file="nsub.txt",action="write") + ! open(89,file="subarea.txt",action="write") + ! open(90,file="subi.txt",action="write") + ! open(91,file="tile_area.txt",action="write") + ! do i=1,nTiles + ! write(88,*)route%nsub(i) + ! write(89,'(150(1x,f10.4))')route%subarea(:,i) + ! write(90,'(150(i7))')route%subi(:,i) + ! write(91,*)route%tile_area(i) + ! enddo + ! stop + !endif + deallocate(ims) call MAPL_GenericInitialize ( GC, import, export, clock, rc=status ) VERIFY_(STATUS) - RETURN_(ESMF_SUCCESS) + RETURN_(ESMF_SUCCESS) end subroutine INITIALIZE ! ----------------------------------------------------------- ! RUN -- Run method for the route component ! ----------------------------------------------------------- + subroutine RUN1 (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 + end subroutine RUN1 - subroutine RUN (GC,IMPORT, EXPORT, CLOCK, RC ) + + subroutine RUN2 (GC,IMPORT, EXPORT, CLOCK, RC ) ! ----------------------------------------------------------- ! !ARGUMENTS: @@ -559,7 +658,7 @@ subroutine RUN (GC,IMPORT, EXPORT, CLOCK, RC ) ! ErrLog Variables ! ----------------------------------------------------------- - character(len=ESMF_MAXSTR) :: IAm="Run" + character(len=ESMF_MAXSTR) :: IAm="Run2" integer :: STATUS character(len=ESMF_MAXSTR) :: COMP_NAME @@ -578,6 +677,7 @@ subroutine RUN (GC,IMPORT, EXPORT, CLOCK, RC ) ! ----------------------------------------------------- real, dimension(:), pointer :: RUNOFF + real, dimension(:), pointer :: RUNOFF_SRC0 ! ----------------------------------------------------- ! INTERNAL pointers @@ -607,6 +707,7 @@ subroutine RUN (GC,IMPORT, EXPORT, CLOCK, RC ) type(ESMF_Grid) :: TILEGRID type (MAPL_LocStream) :: LOCSTREAM + integer :: NTILES, N_CatL, N_CYC logical, save :: FirstTime=.true. real, pointer, dimension(:) :: tile_area @@ -615,502 +716,403 @@ subroutine RUN (GC,IMPORT, EXPORT, CLOCK, RC ) INTEGER, DIMENSION(:,:), POINTER, SAVE :: AllActive,DstCatchID INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: srcProcsID, LocDstCatchID integer, dimension (:),allocatable, SAVE :: GlbActive - INTEGER, SAVE :: N_Active, ThisCycle + INTEGER, SAVE :: N_Active, ThisCycle=1 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 + LENGSC_ACT, QSFLOW_ACT,QOUTFLOW_ACT 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 + INTEGER, DIMENSION(:) ,ALLOCATABLE :: scounts, scounts_global,rdispls, rcounts + real, dimension(:), pointer :: runoff_global,runoff_local,area_local,runoff_cat_global + + integer :: mpierr, nt_global,nt_local, it, j, upid,cid,temp(1),tid,istat + 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 :: runoff_save(:)=>NULL() + real,pointer :: WSTREAM_ACT(:)=>NULL() + real,pointer :: WRIVER_ACT(:)=>NULL() + real,allocatable :: runoff_save_m3(:),runoff_global_m3(:),QOUTFLOW_GLOBAL(:) + real,allocatable :: WTOT_BEFORE(:),WTOT_AFTER(:),QINFLOW_LOCAL(:),UNBALANCE(:),UNBALANCE_GLOBAL(:),ERROR(:),ERROR_GLOBAL(:) + real,allocatable :: QFLOW_SINK(:),QFLOW_SINK_GLOBAL(:),WTOT_BEFORE_GLOBAL(:),WTOT_AFTER_GLOBAL(:) + real,allocatable :: wriver_global(:),wstream_global(:),qsflow_global(:) + ! ------------------ + ! 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" + VERIFY_(STATUS) + Iam = trim(COMP_NAME) // "RUN2" ! Get my internal MAPL_Generic state ! ----------------------------------------------------------- call MAPL_GetObjectFromGC(GC, MAPL, STATUS) VERIFY_(STATUS) - call MAPL_Get(MAPL, HEARTBEAT = HEARTBEAT, RC=STATUS) VERIFY_(STATUS) - + !if (mapl_am_I_root()) print *, "HEARTBEAT=",HEARTBEAT ! Start timers ! ------------ - call MAPL_TimerOn(MAPL,"RUN") - + call MAPL_TimerOn(MAPL,"RUN2") ! Get parameters from generic state ! --------------------------------- - call MAPL_Get(MAPL, INTERNAL_ESMF_STATE=INTERNAL, RC=STATUS) - VERIFY_(STATUS) - + ! call MAPL_Get(MAPL, INTERNAL_ESMF_STATE=INTERNAL, RC=STATUS) + ! VERIFY_(STATUS) ! get pointers to inputs variables ! ---------------------------------- + ndes = route%ndes + mype = route%mype + ntiles = route%ntiles + nt_global = route%nt_global + runoff_save => route%runoff_save + nt_local = route%nt_local + ! 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 -! ----------------------- + VERIFY_(STATUS) + call ESMF_FieldGet(runoff_src, farrayPtr=RUNOFF_SRC0, rc=status) + VERIFY_(STATUS) - 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) + 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) + ! For efficiency, the time step to call the river routing model is set at ROUTE_DT - N_Active = count (srcProcsID == MYPE) + N_CYC = ROUTE_DT/HEARTBEAT + RUN_MODEL : if (ThisCycle == N_CYC) then - allocate (GlbActive(1 : N_Active)) - allocate (tmp_index(1 : N_CatG )) + runoff_save = runoff_save + RUNOFF_SRC0/real (N_CYC) - forall (N=1:N_CatG) tmp_index(N) = N + 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 - GlbActive = pack (tmp_index, mask = (srcProcsID == MYPE)) + allocate(runoff_global(nt_global)) + call MPI_allgatherv ( & + runoff_save, route%scounts_global(mype+1) ,MPI_REAL, & + runoff_global, route%scounts_global, route%rdispls_global,MPI_REAL, & + MPI_COMM_WORLD, mpierr) - ! Initialize the cycle counter and sum (runoff) + if(FirstTime.and.mapl_am_I_root()) print *,"nmax=",nmax + allocate(RUNOFF_ACT(ntiles)) + RUNOFF_ACT=0. + do i=1,ntiles + do j=1,nmax + it=route%subi(j,i) + if(it>0)then + RUNOFF_ACT(i)=RUNOFF_ACT(i)+route%subarea(j,i)*runoff_global(it)/1000. + endif + if(it==0)exit + enddo + enddo - allocate (runoff_save (1:NTILES)) + deallocate(runoff_global) - runoff_save = 0. - ThisCycle = 1 - FirstTime = .false. + allocate (AREACAT_ACT (1:ntiles)) + allocate (LENGSC_ACT (1:ntiles)) + allocate (QSFLOW_ACT (1:ntiles)) + allocate (QOUTFLOW_ACT(1:ntiles)) - deallocate (tmp_index) - - ENDIF FIRST_TIME + LENGSC_ACT=route%lengsc/1.e3 !m->km + AREACAT_ACT=route%areacat/1.e6 !m2->km2 - ! For efficiency, the time step to call the river routing model is set at ROUTE_DT + WSTREAM_ACT => route%wstream + WRIVER_ACT => route%wriver - N_CYC = ROUTE_DT/HEARTBEAT + + allocate(WTOT_BEFORE(ntiles)) + WTOT_BEFORE=WSTREAM_ACT+WRIVER_ACT - 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) + ! ------------------------ + !CALL RIVER_ROUTING_LIN (ntiles, RUNOFF_ACT,AREACAT_ACT,LENGSC_ACT, & + ! WSTREAM_ACT,WRIVER_ACT, QSFLOW_ACT,QOUTFLOW_ACT) + + CALL RIVER_ROUTING_HYD (ntiles, & + RUNOFF_ACT, route%lengsc, route%lstr, & + route%qstr_clmt, route%qri_clmt, route%qin_clmt, & + route%K, route%Kstr, & + WSTREAM_ACT,WRIVER_ACT, & + QSFLOW_ACT,QOUTFLOW_ACT) + + allocate(QOUTFLOW_GLOBAL(n_catg)) + call MPI_allgatherv ( & + QOUTFLOW_ACT, route%scounts_cat(mype+1) ,MPI_REAL, & + QOUTFLOW_GLOBAL, route%scounts_cat, route%rdispls_cat,MPI_REAL, & + MPI_COMM_WORLD, mpierr) + + allocate(QINFLOW_LOCAL(ntiles)) + QINFLOW_LOCAL=0. + do i=1,nTiles + do j=1,upmax + if(route%upid(j,i)>0)then + upid=route%upid(j,i) + WRIVER_ACT(i)=WRIVER_ACT(i)+QOUTFLOW_GLOBAL(upid)*real(route_dt) + QINFLOW_LOCAL(i)=QINFLOW_LOCAL(i)+QOUTFLOW_GLOBAL(upid) + else + exit + endif + enddo + enddo + + !call check_balance(route,ntiles,nt_local,runoff_save,WRIVER_ACT,WSTREAM_ACT,WTOT_BEFORE,RUNOFF_ACT,QINFLOW_LOCAL,QOUTFLOW_ACT,FirstTime,yr_s,mon_s) + + if(FirstTime) nstep_per_day = 86400/route_dt + route%wriver_acc = route%wriver_acc + WRIVER_ACT/real(nstep_per_day) + route%wstream_acc = route%wstream_acc + WSTREAM_ACT/real(nstep_per_day) + route%qoutflow_acc = route%qoutflow_acc + QOUTFLOW_ACT/real(nstep_per_day) + route%qsflow_acc = route%qsflow_acc + QSFLOW_ACT/real(nstep_per_day) + + deallocate(RUNOFF_ACT,AREACAT_ACT,LENGSC_ACT,QOUTFLOW_ACT,QINFLOW_LOCAL,QOUTFLOW_GLOBAL,QSFLOW_ACT,WTOT_BEFORE) + !initialize the cycle counter and sum (runoff_tile) + WSTREAM_ACT=>NULL() + WRIVER_ACT=>NULL() runoff_save = 0. - ThisCycle = 1 + ThisCycle = 1 + ! output + !if(mapl_am_I_root())print *, "nstep_per_day=",nstep_per_day + if(mapl_am_I_root())print *, "Current time is ", YY, "/", MM, "/", DD, " ", HH, ":", MMM, ":", SS, ", next MM_next:",MM_next + if(FirstTime)then + if(mapl_am_I_root()) istat = mkdir("../river", int(o'755',c_int16_t)) + endif + if(HH==23)then + allocate(wriver_global(n_catg),wstream_global(n_catg),qoutflow_global(n_catg),qsflow_global(n_catg)) + !call MPI_allgatherv ( & + ! route%wriver_acc, route%scounts_cat(mype+1) ,MPI_REAL, & + ! wriver_global, route%scounts_cat, route%rdispls_cat,MPI_REAL, & + ! MPI_COMM_WORLD, mpierr) + !call MPI_allgatherv ( & + ! route%wstream_acc, route%scounts_cat(mype+1) ,MPI_REAL, & + ! wstream_global, route%scounts_cat, route%rdispls_cat,MPI_REAL, & + ! MPI_COMM_WORLD, mpierr) + call MPI_allgatherv ( & + route%qoutflow_acc, route%scounts_cat(mype+1) ,MPI_REAL, & + qoutflow_global, route%scounts_cat, route%rdispls_cat,MPI_REAL, & + MPI_COMM_WORLD, mpierr) + !call MPI_allgatherv ( & + ! route%qsflow_acc, route%scounts_cat(mype+1) ,MPI_REAL, & + ! qsflow_global, route%scounts_cat, route%rdispls_cat,MPI_REAL, & + ! MPI_COMM_WORLD, mpierr) + if(mapl_am_I_root())then + !open(88,file="../river/river_storage_"//trim(yr_s)//trim(mon_s)//trim(day_s)//".txt",action="write") + !open(89,file="../river/stream_storage_"//trim(yr_s)//trim(mon_s)//trim(day_s)//".txt",action="write") + open(90,file="../river/river_flow_"//trim(yr_s)//trim(mon_s)//trim(day_s)//".txt",action="write") + !open(91,file="../river/stream_flow_"//trim(yr_s)//trim(mon_s)//trim(day_s)//".txt",action="write") + do i=1,n_catg + !write(88,*)wriver_global(i) + !write(89,*)wstream_global(i) + write(90,*)qoutflow_global(i) + !write(91,*)qsflow_global(i) + enddo + close(90) + !close(88);close(89);close(90)!;close(91) + !print *, "output river storage is: ",sum(wriver_global)/1.e9 + !print *, "output stream storage is: ",sum(wstream_global)/1.e9 + endif + deallocate(wriver_global,wstream_global,qoutflow_global,qsflow_global) + route%wriver_acc = 0. + route%wstream_acc = 0. + route%qoutflow_acc = 0. + route%qsflow_acc = 0. + endif + + !restart + if(MM_next/=MM)then + allocate(wriver_global(n_catg),wstream_global(n_catg)) + call MPI_allgatherv ( & + route%wstream, route%scounts_cat(mype+1) ,MPI_REAL, & + wstream_global, route%scounts_cat, route%rdispls_cat,MPI_REAL, & + MPI_COMM_WORLD, mpierr) + call MPI_allgatherv ( & + route%wriver, route%scounts_cat(mype+1) ,MPI_REAL, & + wriver_global, route%scounts_cat, route%rdispls_cat,MPI_REAL, & + MPI_COMM_WORLD, mpierr) + if(mapl_am_I_root())then + write(yr_s,'(I4.4)')YY_next + write(mon_s,'(I2.2)')MM_next + write(day_s,'(I2.2)')DD_next + open(88,file="../input/restart/river_storage_rs_"//trim(yr_s)//trim(mon_s)//trim(day_s)//".txt",action="write") + open(89,file="../input/restart/stream_storage_rs_"//trim(yr_s)//trim(mon_s)//trim(day_s)//".txt",action="write") + do i=1,n_catg + write(88,*)wriver_global(i) + write(89,*)wstream_global(i) + enddo + close(88);close(89) + print *, "saved river storage is: ",sum(wriver_global)/1.e9 + print *, "saved stream storage is: ",sum(wstream_global)/1.e9 + endif + deallocate(wriver_global,wstream_global) + endif + + if(FirstTime) FirstTime=.False. else - runoff_save = runoff_save + runoff/real (N_CYC) + runoff_save = runoff_save + RUNOFF_SRC0/real (N_CYC) ThisCycle = ThisCycle + 1 - - endif RUN_MODEL - call MAPL_TimerOff ( MAPL, "-RRM" ) + endif RUN_MODEL + + runoff_save => NULL() ! All done ! -------- + call MAPL_TimerOff ( MAPL, "-RRM" ) + call MAPL_TimerOff(MAPL,"RUN2") + !call MPI_Barrier(MPI_COMM_WORLD, mpierr) - call MAPL_TimerOff(MAPL,"RUN") RETURN_(ESMF_SUCCESS) + end subroutine RUN2 - 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' + subroutine check_balance(route,ntiles,nt_local,runoff_save,WRIVER_ACT,WSTREAM_ACT,WTOT_BEFORE,RUNOFF_ACT,QINFLOW_LOCAL,QOUTFLOW_ACT,FirstTime,yr_s,mon_s) + + type(T_RROUTE_STATE), intent(in) :: route + integer, intent(in) :: ntiles,nt_local + real,intent(in) :: runoff_save(nt_local),WRIVER_ACT(ntiles),WSTREAM_ACT(ntiles),WTOT_BEFORE(ntiles),RUNOFF_ACT(ntiles) + real,intent(in) :: QINFLOW_LOCAL(ntiles),QOUTFLOW_ACT(ntiles) + logical,intent(in) :: FirstTime + character(len=*), intent(in) :: yr_s,mon_s + + real,allocatable :: runoff_cat_global(:) + real,allocatable :: runoff_save_m3(:),runoff_global_m3(:) + real,allocatable :: WTOT_AFTER(:),UNBALANCE(:),UNBALANCE_GLOBAL(:),ERROR(:),ERROR_GLOBAL(:) + real,allocatable :: QFLOW_SINK(:),QFLOW_SINK_GLOBAL(:),WTOT_BEFORE_GLOBAL(:),WTOT_AFTER_GLOBAL(:) - ! STEP 1: Identify active catchments within the local processor. If the catchment is active in - ! more than 1 processor, choose an owner. - ! -------------------------------------------------------------------------------------------- + integer :: i, nt_global,mype,cid,temp(1),tid,mpierr + real :: wr_error, wr_tot, runf_tot - 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 + nt_global = route%nt_global + mype = route%mype - allocate (global_buff (N_CatG * numprocs)) - allocate (scounts(numprocs),rdispls(numprocs),rcounts(numprocs)) + allocate(WTOT_AFTER(ntiles),UNBALANCE(ntiles),UNBALANCE_GLOBAL(n_catg),runoff_cat_global(n_catg)) + allocate(QFLOW_SINK(ntiles),QFLOW_SINK_GLOBAL(n_catg),WTOT_BEFORE_GLOBAL(n_catg),WTOT_AFTER_GLOBAL(n_catg)) + allocate(runoff_save_m3(nt_local),runoff_global_m3(nt_global),ERROR(ntiles),ERROR_GLOBAL(n_catg)) - 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 + WTOT_AFTER=WRIVER_ACT+WSTREAM_ACT + ERROR = WTOT_AFTER - (WTOT_BEFORE + RUNOFF_ACT*route_dt + QINFLOW_LOCAL*route_dt - QOUTFLOW_ACT*route_dt) + where(QOUTFLOW_ACT>0.) UNBALANCE = abs(ERROR)/(QOUTFLOW_ACT*route_dt) + where(QOUTFLOW_ACT<=0.) UNBALANCE = 0. + call MPI_allgatherv ( & + UNBALANCE, route%scounts_cat(mype+1) ,MPI_REAL, & + UNBALANCE_GLOBAL, route%scounts_cat, route%rdispls_cat,MPI_REAL, & + MPI_COMM_WORLD, mpierr) + QFLOW_SINK=0. + do i=1,ntiles + if(route%downid(i)==-1)then + QFLOW_SINK(i) = QOUTFLOW_ACT(i) + endif + enddo + call MPI_allgatherv ( & + QFLOW_SINK, route%scounts_cat(mype+1) ,MPI_REAL, & + QFLOW_SINK_GLOBAL, route%scounts_cat, route%rdispls_cat,MPI_REAL, & + MPI_COMM_WORLD, mpierr) + call MPI_allgatherv ( & + WTOT_BEFORE, route%scounts_cat(mype+1) ,MPI_REAL, & + WTOT_BEFORE_GLOBAL, route%scounts_cat, route%rdispls_cat,MPI_REAL, & + MPI_COMM_WORLD, mpierr) + call MPI_allgatherv ( & + WTOT_AFTER, route%scounts_cat(mype+1) ,MPI_REAL, & + WTOT_AFTER_GLOBAL, route%scounts_cat, route%rdispls_cat,MPI_REAL, & + MPI_COMM_WORLD, mpierr) + runoff_save_m3=runoff_save*route%tile_area/1000. + call MPI_allgatherv ( & + runoff_save_m3, route%scounts_global(mype+1) ,MPI_REAL, & + runoff_global_m3, route%scounts_global, route%rdispls_global,MPI_REAL, & + MPI_COMM_WORLD, mpierr) + call MPI_allgatherv ( & + RUNOFF_ACT, route%scounts_cat(mype+1) ,MPI_REAL, & + runoff_cat_global, route%scounts_cat, route%rdispls_cat,MPI_REAL, & + MPI_COMM_WORLD, mpierr) + if(mapl_am_I_root())then + open(88,file="../runoff_tile_global_"//trim(yr_s)//"_"//trim(mon_s)//".txt",status="unknown", position="append") + write(88,*)sum(runoff_global_m3) + close(88) + open(88,file="../runoff_cat_global_"//trim(yr_s)//"_"//trim(mon_s)//".txt",status="unknown", position="append") + write(88,*)sum(runoff_cat_global) + close(88) + !print *,"sum(runoff_global_m3)=",sum(runoff_global_m3) + !print *,"sum(runoff_cat_global)=",sum(runoff_cat_global) + endif + if(mapl_am_I_root())then + open(88,file="../WTOT_AFTER_"//trim(yr_s)//"_"//trim(mon_s)//".txt",status="unknown", position="append") + write(88,*)sum(WTOT_AFTER_GLOBAL) + close(88) + open(88,file="../WTOT_BEFORE_RUNOFF_QSINK_"//trim(yr_s)//"_"//trim(mon_s)//".txt",status="unknown", position="append") + write(88,*) sum(WTOT_BEFORE_GLOBAL)+sum(runoff_global_m3)*route_dt-sum(QFLOW_SINK_GLOBAL)*route_dt + close(88) + wr_error=sum(WTOT_AFTER_GLOBAL)-(sum(WTOT_BEFORE_GLOBAL)+sum(runoff_global_m3)*route_dt-sum(QFLOW_SINK_GLOBAL)*route_dt) + runf_tot=sum(runoff_global_m3)*route_dt + wr_tot=sum(WTOT_AFTER_GLOBAL) + open(88,file="../WTOT_ERROR_2_RUNOFF_"//trim(yr_s)//"_"//trim(mon_s)//".txt",status="unknown", position="append") + write(88,*) wr_error/runf_tot + close(88) + open(88,file="../WTOT_ERROR_2_WTOT_"//trim(yr_s)//"_"//trim(mon_s)//".txt",status="unknown", position="append") + write(88,*) wr_error/wr_tot + close(88) + !print *,"WTOT_ERROR_2_RUNOFF:",(sum(WTOT_AFTER_GLOBAL)-(sum(WTOT_BEFORE_GLOBAL)+sum(runoff_global_m3)*route_dt-sum(QFLOW_SINK_GLOBAL)*route_dt))/(sum(runoff_global_m3)*route_dt) + endif + + call MPI_allgatherv ( & + ERROR, route%scounts_cat(mype+1) ,MPI_REAL, & + ERROR_GLOBAL, route%scounts_cat, route%rdispls_cat,MPI_REAL, & + MPI_COMM_WORLD, mpierr) + temp = maxloc(abs(ERROR_GLOBAL)) + cid = temp(1) + if(cid>=route%minCatch.and.cid<=route%maxCatch)then + tid=cid-route%minCatch+1 + print *,"my PE is:",mype,", max abs value of ERROR=", ERROR(tid)," at pfafid: ",route%minCatch+tid-1,", W_BEFORE=",WTOT_BEFORE(tid),", RUNOFF=",RUNOFF_ACT(tid)*route_dt,", QINFLOW=",QINFLOW_LOCAL(tid)*route_dt,", QOUTFLOW=",QOUTFLOW_ACT(tid)*route_dt,", W_AFTER=",WTOT_AFTER(tid) + endif + !if(FirstTime)then + ! if(mapl_am_I_root())then + ! open(88,file="ERROR_TOTAL.txt",action="write") + ! do i=1,n_catg + ! write(88,*)ERROR_GLOBAL(i) + ! enddo + ! endif + !endif + + deallocate(WTOT_AFTER,UNBALANCE,UNBALANCE_GLOBAL,ERROR,QFLOW_SINK,QFLOW_SINK_GLOBAL,WTOT_BEFORE_GLOBAL,WTOT_AFTER_GLOBAL) + deallocate(runoff_save_m3,runoff_global_m3,ERROR_GLOBAL,runoff_cat_global) + + + end subroutine check_balance - 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/offline/interp_M36toPfaf.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/offline/interp_M36toPfaf.f90 new file mode 100644 index 000000000..e3b831370 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/offline/interp_M36toPfaf.f90 @@ -0,0 +1,157 @@ +module interp + +use omp_lib ! Use OpenMP library for parallel processing +use rwncfile ! Use custom module for reading NetCDF files +implicit none + +private +public :: M36_to_cat ! Make the M36_to_cat function public +public :: M09_to_cat ! Make the M09_to_cat function public + +contains + +!------------------------------------------------------------------------------ +! This function maps runoff data from M36 resolution to catchments (cat) +function M36_to_cat(runoff,nlon,nlat,ncat,inputdir) result(Qrunf) + + integer,intent(in) :: nlon,nlat,ncat ! Input: number of longitude, latitude, and catchments + real*8,intent(in) :: runoff(nlon,nlat) ! Input: runoff array of size (nlon, nlat) + character(len=500),intent(in) :: inputdir ! Input: directory path for input files + real*8 :: Qrunf(ncat) ! Output: runoff mapped to catchments + + real*8,parameter :: small=1.D-12 ! Small value to avoid division by zero + + integer,parameter :: nmax=150 ! Maximum number of sub-areas per catchment + integer,parameter :: nc=291284 ! Total number of catchments + + real*8,allocatable,dimension(:,:) :: subarea,frac ! Arrays for sub-area and fractions + integer,allocatable,dimension(:,:) :: subx,suby ! Arrays for x and y coordinates of sub-areas + real*8,allocatable,dimension(:) :: tot,runfC,fracA ! Arrays for total area, calculated runoff, and fraction + integer,allocatable,dimension(:) :: nsub ! Array for number of sub-areas per catchment + + integer :: i,j,sx,sy ! Loop variables and coordinates for sub-areas + + ! Allocate memory for arrays + allocate(nsub(nc),subarea(nmax,nc),subx(nmax,nc),suby(nmax,nc),tot(nc)) + + ! Read sub-area data from text files + open(77,file=trim(inputdir)//"/Pfaf_nsub_M36.txt"); read(77,*)nsub + open(77,file=trim(inputdir)//"/Pfaf_asub_M36.txt"); read(77,*)subarea + open(77,file=trim(inputdir)//"/Pfaf_xsub_M36.txt"); read(77,*)subx + open(77,file=trim(inputdir)//"/Pfaf_ysub_M36.txt"); read(77,*)suby + open(77,file=trim(inputdir)//"/Pfaf_area.txt"); read(77,*)tot + + ! Allocate memory for fraction array + allocate(frac(nmax,nc)) + + ! Compute fraction of each sub-area relative to the total catchment area + do i=1,nc + frac(:,i)=subarea(:,i)/tot(i) + enddo + + ! Allocate memory for runoff and fraction arrays + allocate(runfC(nc),fracA(nc)) + runfC=0.D0 ! Initialize runoff array to zero + fracA=0.D0 ! Initialize fraction array to zero + + !$OMP PARALLEL default(shared) private(i,j,sx,sy) ! Start OpenMP parallel region + !$OMP DO + ! Loop over all catchments and sub-areas + do i=1,nc + if(nsub(i)>=1)then + 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 + ! Check for valid fraction and runoff values + if(frac(j,i)>0.D0.and.runoff(sx,sy)<1.D14)then + runfC(i)=runfC(i)+frac(j,i)*runoff(sx,sy) ! Accumulate runoff for the catchment + fracA(i)=fracA(i)+frac(j,i) ! Accumulate fraction + endif + enddo + endif + enddo + !$OMP END DO + !$OMP END PARALLEL ! End OpenMP parallel region + + ! Convert to kg/s by multiplying by area (in m²) and dividing by time (in seconds) + Qrunf=runfC*(tot*1.D6)/86400.D0 + + ! Deallocate arrays to free memory + deallocate(subarea,subx,suby,tot,frac,& + runfC,fracA,nsub) + +end function M36_to_cat +!------------------------------------------------------------------------------ + +!------------------------------------------------------------------------------ +! This function maps runoff data from M09 resolution to catchments (cat) +function M09_to_cat(runoff,nlon,nlat,ncat,inputdir) result(Qrunf) + + integer,intent(in) :: nlon,nlat,ncat ! Input: number of longitude, latitude, and catchments + real*8,intent(in) :: runoff(nlon,nlat) ! Input: runoff array of size (nlon, nlat) + character(len=500),intent(in) :: inputdir ! Input: directory path for input files + real*8 :: Qrunf(ncat) ! Output: runoff mapped to catchments + + real*8,parameter :: small=1.D-12 ! Small value to avoid division by zero + + integer,parameter :: nmax=458 ! Maximum number of sub-areas per catchment + integer,parameter :: nc=291284 ! Total number of catchments + + real*8,allocatable,dimension(:,:) :: subarea,frac ! Arrays for sub-area and fractions + integer,allocatable,dimension(:,:) :: subx,suby ! Arrays for x and y coordinates of sub-areas + real*8,allocatable,dimension(:) :: tot,runfC,fracA ! Arrays for total area, calculated runoff, and fraction + integer,allocatable,dimension(:) :: nsub ! Array for number of sub-areas per catchment + + integer :: i,j,sx,sy ! Loop variables and coordinates for sub-areas + + ! Allocate memory for arrays + allocate(nsub(nc),subarea(nmax,nc),subx(nmax,nc),suby(nmax,nc),tot(nc)) + + ! Read sub-area data from text files + open(77,file=trim(inputdir)//"/Pfaf_nsub_M09.txt"); read(77,*)nsub + open(77,file=trim(inputdir)//"/Pfaf_asub_M09.txt"); read(77,*)subarea + open(77,file=trim(inputdir)//"/Pfaf_xsub_M09.txt"); read(77,*)subx + open(77,file=trim(inputdir)//"/Pfaf_ysub_M09.txt"); read(77,*)suby + open(77,file=trim(inputdir)//"/Pfaf_area.txt"); read(77,*)tot + + ! Allocate memory for fraction array + allocate(frac(nmax,nc)) + + ! Compute fraction of each sub-area relative to the total catchment area + do i=1,nc + frac(:,i)=subarea(:,i)/tot(i) + enddo + + ! Allocate memory for runoff and fraction arrays + allocate(runfC(nc),fracA(nc)) + runfC=0.D0 ! Initialize runoff array to zero + fracA=0.D0 ! Initialize fraction array to zero + + !$OMP PARALLEL default(shared) private(i,j,sx,sy) ! Start OpenMP parallel region + !$OMP DO + ! Loop over all catchments and 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 + ! Check for valid fraction and runoff values + if(frac(j,i)>0.D0.and.runoff(sx,sy)<1.D14.and.runoff(sx,sy)>=0.D0)then + runfC(i)=runfC(i)+frac(j,i)*runoff(sx,sy) ! Accumulate runoff for the catchment + fracA(i)=fracA(i)+frac(j,i) ! Accumulate fraction + endif + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL ! End OpenMP parallel region + + ! Convert to kg/s by multiplying by area (in m²) and dividing by time (in seconds) + Qrunf=runfC*(tot*1.D6)/86400.D0 + + ! Deallocate arrays to free memory + deallocate(subarea,subx,suby,tot,frac,& + runfC,fracA,nsub) + +end function M09_to_cat +!------------------------------------------------------------------------------ + +end module interp \ No newline at end of file diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/offline/lake_mod.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/offline/lake_mod.f90 new file mode 100644 index 000000000..0aab6d316 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/offline/lake_mod.f90 @@ -0,0 +1,109 @@ +module lake + + +implicit none +private +public :: lake_init, lake_cal + +! Define parameters for small and large lakes +real*8, parameter :: fac_a_slake = 0.003D0 ! Factor for small lakes +real*8, parameter :: fac_b_slake = 0.40D0 ! Exponent for small lakes +real*8, parameter :: fac_a_llake = 0.01D0 ! Factor for large lakes +real*8, parameter :: fac_b_llake = 0.60D0 ! Exponent for large lakes +real*8, parameter :: thr_area_lake = 1D4 ! Threshold lake area (in km^2) + +! Define constants +real*8, parameter :: dt = 86400.D0 ! Time step in seconds (1 day) +real*8, parameter :: rho = 1.D3 ! Water density in kg/m^3 + +contains + +!------------------------------ +! Initialization subroutine for lakes +subroutine lake_init(input_dir, use_lake, nc, nlake, nres, active_res, active_lake, area_lake, Wr_lake, Q_lake) + character(len=500),intent(in) :: input_dir + logical, intent(in) :: use_lake ! Flag to use lake module + integer, intent(in) :: nc, nlake, nres ! Number of catchments, lakes, reservoirs + integer, intent(in) :: active_res(nres) ! Active reservoirs + integer, allocatable, intent(inout) :: active_lake(:) ! Active lakes (output) + real*8, allocatable, intent(inout) :: area_lake(:), Wr_lake(:), Q_lake(:) ! Lake areas, water storage, outflow + + integer, allocatable :: flag_valid_laked(:), catid_laked(:) + real*8, allocatable :: area_laked(:) + + integer :: i, cid + + ! Allocate arrays for lake attributes + allocate(flag_valid_laked(nlake), catid_laked(nlake), area_laked(nlake)) + allocate(active_lake(nc), area_lake(nc)) + allocate(Wr_lake(nc), Q_lake(nc)) + + ! Read lake outlet and area data from external files + open(77, file = trim(input_dir)//"/lake_outlet_flag_valid_2097.txt") + read(77, *) flag_valid_laked + open(77, file = trim(input_dir)//"/lake_outlet_catid.txt") + read(77, *) catid_laked + open(77, file = trim(input_dir)//"/lake_outlet_lakearea.txt") + read(77, *) area_laked ! km^2 + + ! Initialize lake attributes to zero + area_lake = 0.D0 + active_lake = 0 + + ! Assign active lakes and their areas based on data + do i = 1, nlake + if (flag_valid_laked(i) == 1) then + cid = catid_laked(i) + active_lake(cid) = 1 + area_lake(cid) = area_laked(i) + endif + enddo + + ! Deactivate lakes where reservoirs are active + where (active_res == 1) active_lake = 0 + + ! If lakes are not being used, set active lakes to zero + if (use_lake == .False.) active_lake = 0 + +end subroutine lake_init + +!------------------------------ +! Calculation subroutine for lakes +subroutine lake_cal(active_lake, area_lake, Q_lake, Wr_lake, Qout, B1, B2) + integer, intent(in) :: active_lake ! Flag indicating if lake is active + real*8, intent(in) :: area_lake, Qout ! Lake area, outlet flow rate + real*8, intent(inout) :: Q_lake, Wr_lake ! Lake inflow, water storage + real*8, intent(inout) :: B1, B2 ! Output variables (Q_lake, some other parameter) + + real*8 :: alp_lake ! Alpha parameter for lake flow calculation + + ! Process only active lakes + if (active_lake == 1) then + + ! Determine lake type based on area and calculate alpha + if (area_lake >= thr_area_lake) then + alp_lake = fac_a_llake * ( (1.D0 / sqrt(area_lake)) ** fac_b_llake ) / 3600.D0 + else + alp_lake = fac_a_slake * ( (1.D0 / sqrt(area_lake)) ** fac_b_slake ) / 3600.D0 + endif + + ! Compute lake outflow based on alpha and water storage + Q_lake = alp_lake * Wr_lake + + ! Ensure that outflow is non-negative and does not exceed available water + Q_lake = max(0.D0, Q_lake) + Q_lake = min(Q_lake, Wr_lake / dt + Qout) + + ! Update water storage in lake + Wr_lake = Wr_lake + dt * (Qout - Q_lake) + Wr_lake = max(0.D0, Wr_lake) + + ! Assign output values + B1 = Q_lake + B2 = 0.D0 + + endif + +end subroutine lake_cal + +end module lake \ No newline at end of file diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/offline/ncdioMod.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/offline/ncdioMod.f90 new file mode 100644 index 000000000..94b50af1a --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/offline/ncdioMod.f90 @@ -0,0 +1,2355 @@ + +module ncdio + use netcdf +!----------------------------------------------------------------------- +!BOP +! +! !MODULE: ncdioMod +! +! !DESCRIPTION: +! Generic interfaces to write fields to netcdf files +! +! !USES: +! +! !PUBLIC TYPES: + implicit none + include 'netcdf.inc' ! + save + public :: check_ret ! checks return status of netcdf calls + public :: check_var ! determine if variable is on netcdf file + public :: check_dim ! validity check on dimension + public :: ncd_defvar +! +! !REVISION HISTORY: +! +!EOP +! +! !PRIVATE METHODS: +! + interface ncd_iolocal + module procedure ncd_iolocal_int_1d + module procedure ncd_iolocal_real_1d + module procedure ncd_iolocal_double_1d + module procedure ncd_iolocal_int_2d + module procedure ncd_iolocal_real_2d + module procedure ncd_iolocal_double_2d + end interface + + interface ncd_ioglobal + module procedure ncd_ioglobal_int_var + module procedure ncd_ioglobal_real_var + module procedure ncd_ioglobal_double_var + module procedure ncd_ioglobal_int_1d + module procedure ncd_ioglobal_real_1d + module procedure ncd_ioglobal_double_1d + module procedure ncd_ioglobal_byte_2d + module procedure ncd_ioglobal_short_2d + module procedure ncd_ioglobal_int_2d + module procedure ncd_ioglobal_long_2d + module procedure ncd_ioglobal_real_2d + module procedure ncd_ioglobal_double_2d + module procedure ncd_ioglobal_int_3d + module procedure ncd_ioglobal_short_3d + module procedure ncd_ioglobal_real_3d + module procedure ncd_ioglobal_double_3d + end interface + + private :: endrun + logical, public, parameter :: nc_masterproc = .true. ! proc 0 logical for printing msgs + +!----------------------------------------------------------------------- + +contains + +!----------------------------------------------------------------------- +!BOP +! +! !IROUTINE: check_dim +! +! !INTERFACE: + subroutine check_dim(ncid, dimname, value) +! +! !DESCRIPTION: +! Validity check on dimension +! !ARGUMENTS: + implicit none + integer, intent(in) :: ncid + character(len=*), intent(in) :: dimname + integer, intent(in) :: value +! +! !REVISION HISTORY: +! +!EOP +! +! !LOCAL VARIABLES: + integer :: dimid, dimlen ! temporaries +!----------------------------------------------------------------------- + + call check_ret(nf_inq_dimid (ncid, trim(dimname), dimid), 'check_dim') + call check_ret(nf_inq_dimlen (ncid, dimid, dimlen), 'check_dim') + if (dimlen /= value) then + write (6,*) 'CHECK_DIM error: mismatch of input dimension ',dimlen, & + ' with expected value ',value,' for variable ',trim(dimname) + call endrun() + end if + + end subroutine check_dim + +!----------------------------------------------------------------------- +!BOP +! +! !IROUTINE: check_var +! +! !INTERFACE: + subroutine check_var(ncid, varname, varid, readvar) +! !DESCRIPTION: +! Check if variable is on netcdf file +! +! !ARGUMENTS: + implicit none + integer, intent(in) :: ncid + character(len=*), intent(in) :: varname + integer, intent(out) :: varid + logical, intent(out) :: readvar +! +! !REVISION HISTORY: +! +!EOP +! +! !LOCAL VARIABLES: + integer :: ret ! return value +!----------------------------------------------------------------------- + + readvar = .true. + if (nc_masterproc) then + ret = nf_inq_varid (ncid, varname, varid) + if (ret/=NF_NOERR) then + write(6,*)'CHECK_VAR: variable ',trim(varname),' is not on initial dataset' + readvar = .false. + end if + end if + end subroutine check_var + +!----------------------------------------------------------------------- +!BOP +! +! !IROUTINE: check_ret +! +! !INTERFACE: + subroutine check_ret(ret, calling) +! !DESCRIPTION: +! Check return status from netcdf call +! +! !ARGUMENTS: + implicit none + integer, intent(in) :: ret + character(len=*) :: calling +! +! !REVISION HISTORY: +! +!EOP +!----------------------------------------------------------------------- + + if (ret /= NF_NOERR) then + write(6,*)'netcdf error from ',trim(calling) + call endrun(nf_strerror(ret)) + end if + + end subroutine check_ret + +!----------------------------------------------------------------------- +!BOP +! +! !IROUTINE: ncd_defvar +! +! !INTERFACE: + subroutine ncd_defvar(ncid, varname, xtype, & + dim1name, dim2name, dim3name, dim4name, dim5name, & + long_name, units, cell_method, missing_value, fill_value, & + imissing_value, ifill_value) +! !DESCRIPTION: +! Define a netcdf variable +! +! !ARGUMENTS: + implicit none + integer , intent(in) :: ncid ! input unit + character(len=*), intent(in) :: varname ! variable name + integer , intent(in) :: xtype ! external type + character(len=*), intent(in), optional :: dim1name ! dimension name + character(len=*), intent(in), optional :: dim2name ! dimension name + character(len=*), intent(in), optional :: dim3name ! dimension name + character(len=*), intent(in), optional :: dim4name ! dimension name + character(len=*), intent(in), optional :: dim5name ! dimension name + character(len=*), intent(in), optional :: long_name ! attribute + character(len=*), intent(in), optional :: units ! attribute + character(len=*), intent(in), optional :: cell_method ! attribute + real , intent(in), optional :: missing_value ! attribute for real + real , intent(in), optional :: fill_value ! attribute for real + integer , intent(in), optional :: imissing_value ! attribute for int + integer , intent(in), optional :: ifill_value ! attribute for int +! +! !REVISION HISTORY: +! +!EOP +! +! !LOCAL VARIABLES: + integer :: n ! indices + integer :: ndims ! dimension counter + integer :: dimid(5) ! dimension ids + integer :: varid ! variable id + integer :: itmp ! temporary + character(len=256) :: str ! temporary + character(len=32) :: subname='NCD_DEFVAR_REAL' ! subroutine name +!----------------------------------------------------------------------- + + if (.not. nc_masterproc) return + + ! Determine dimension ids for variable + + dimid(:) = 0 + ndims=0 + if (present(dim1name)) then + ndims=ndims+1 + call check_ret(nf_inq_dimid(ncid, dim1name, dimid(ndims)), subname) + end if + if (present(dim2name)) then + ndims=ndims+1 + call check_ret(nf_inq_dimid(ncid, dim2name, dimid(ndims)), subname) + end if + if (present(dim3name)) then + ndims=ndims+1 + call check_ret(nf_inq_dimid(ncid, dim3name, dimid(ndims)), subname) + end if + if (present(dim4name)) then + ndims=ndims+1 + call check_ret(nf_inq_dimid(ncid, dim4name, dimid(ndims)), subname) + end if + if (present(dim5name)) then + ndims=ndims+1 + call check_ret(nf_inq_dimid(ncid, dim5name, dimid(ndims)), subname) + end if + + + ! Define variable + + if (present(dim1name) .or. present(dim2name) .or. present(dim3name) .or. & + present(dim4name) .or. present(dim5name)) then + call check_ret(nf_def_var(ncid, trim(varname), xtype, ndims, dimid(1:ndims), varid), subname) + else + call check_ret(nf_def_var(ncid, varname, xtype, 0, 0, varid), subname) + end if + if (present(long_name)) then + call check_ret(nf_put_att_text(ncid, varid, 'long_name', len_trim(long_name), trim(long_name)), subname) + end if + if (present(units)) then + call check_ret(nf_put_att_text(ncid, varid, 'units', len_trim(units), trim(units)), subname) + end if + if (present(cell_method)) then + str = 'time: ' // trim(cell_method) + call check_ret(nf_put_att_text(ncid, varid, 'cell_method', len_trim(str), trim(str)), subname) + end if + if (present(fill_value)) then + call check_ret(nf_put_att_real(ncid, varid, '_FillValue', xtype, 1, fill_value), subname) + end if + if (present(missing_value)) then + call check_ret(nf_put_att_real(ncid, varid, 'missing_value', xtype, 1, missing_value), subname) + end if + if (present(ifill_value)) then + call check_ret(nf_put_att_int(ncid, varid, '_FillValue', xtype, 1, ifill_value), subname) + end if + if (present(imissing_value)) then + call check_ret(nf_put_att_int(ncid, varid, 'missing_value', xtype, 1, imissing_value), subname) + end if + + end subroutine ncd_defvar + +!----------------------------------------------------------------------- +!BOP +! +! !IROUTINE: ncd_iolocal_int_1d +! +! !INTERFACE: + + subroutine ncd_iolocal_int_1d(varname, data, flag, ncid, & + lb_lon, lb_lat, lb_lvl, lb_t, ub_lon, ub_lat, ub_lvl, ub_t, & + long_name, units, readvar) +! !DESCRIPTION: +! I/O for 1d int field +! +! !USES: +! +! !ARGUMENTS: + implicit none + character(len=*), intent(in) :: varname ! variable name + integer , intent(inout) :: data(:) ! local decomposition data + character(len=*), intent(in) :: flag ! 'read' or 'write' + integer , intent(in) :: ncid ! input unit + integer , optional, intent(in) :: lb_lon ! start for longitude + integer , optional, intent(in) :: lb_lat ! start for latitute sizes + integer , optional, intent(in) :: lb_lvl ! start for level size + integer , optional, intent(in) :: lb_t ! start for time size + integer , optional, intent(in) :: ub_lon ! start for longitude + integer , optional, intent(in) :: ub_lat ! start for latitute sizes + integer , optional, intent(in) :: ub_lvl ! start for level size + integer , optional, intent(in) :: ub_t ! start for time size + character(len=*), optional, intent(in) :: long_name ! variable long name + character(len=*), optional, intent(in) :: units ! variable units + logical , optional, intent(out):: readvar ! true => variable is on initial dataset (read only) +! !REVISION HISTORY: +! +!EOP +! +! !LOCAL VARIABLES: + integer :: varid ! variable id + integer :: ndim ! dimension counter + integer :: start(4) ! starting indices for netcdf field + integer :: count(4) ! count values for netcdf field + character(len=32) :: inq_name ! inquid variable name + character(len=8) :: inq_xtype ! inquid variable type + integer :: inq_ndims ! inquid variable dimention + integer :: inq_dimids(4) ! inquid variable dimention id + character(len=255) :: inq_natts ! inquid variable attachment + character(len=32) :: subname='NCD_IOLOCAL_INT_1D' ! subroutine name + logical :: varpresent ! if true, variable is on tape +!----------------------------------------------------------------------- + + ! Write field as 1d field + if (flag == 'write') then + if (nc_masterproc) then + call check_ret(nf_inq_varid(ncid, varname, varid), subname) + ! Write 1d field + ndim=0 + count=1 + if (present(lb_lon) .and. present(ub_lon)) then + ndim=ndim+1 + start(ndim)=lb_lon + count(ndim)=ub_lon-lb_lon+1 + else if(present(lb_lon) .neqv. present(ub_lon))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_lat)) then + ndim=ndim+1 + start(ndim)=lb_lat + count(ndim)=ub_lat-lb_lat+1 + else if(present(lb_lat) .neqv. present(ub_lat))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_lvl)) then + ndim=ndim+1 + start(ndim)=lb_lvl + count(ndim)=ub_lvl-lb_lvl+1 + else if(present(lb_lvl) .neqv. present(ub_lvl))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_t)) then + ndim=ndim+1 + start(ndim)=lb_t + count(ndim)=ub_t-lb_t+1 + else if(present(lb_t) .neqv. present(lb_t))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + + if ((.not. present(lb_lon)) .and. (.not. present(lb_lat)) .and. & + (.not. present(lb_lvl)) .and. (.not. present(lb_t))) then + call endrun('must specify one dimention!',subname) + endif + + call check_ret(nf_put_vara_int(ncid, varid, start(1:ndim), count(1:ndim), data), subname) + if (present(long_name)) then + call check_ret(nf_put_att_text(ncid, varid, 'long_name', len_trim(long_name), trim(long_name)), subname) + end if + if (present(units)) then + call check_ret(nf_put_att_text(ncid, varid, 'units', len_trim(units), trim(units)), subname) + end if + end if ! end of if-nc_masterproc block + ! Read field as 1d field + else if (flag == 'read') then + if (nc_masterproc) then + call check_var(ncid, varname, varid, varpresent) + if (varpresent) then + ndim=0 + count=1 + if (present(lb_lon) .and. present(ub_lon)) then + ndim=ndim+1 + start(ndim)=lb_lon + count(ndim)=ub_lon-lb_lon+1 + else if(present(lb_lon) .neqv. present(ub_lon))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_lat)) then + ndim=ndim+1 + start(ndim)=lb_lat + count(ndim)=ub_lat-lb_lat+1 + else if(present(lb_lat) .neqv. present(ub_lat))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_lvl)) then + ndim=ndim+1 + start(ndim)=lb_lvl + count(ndim)=ub_lvl-lb_lvl+1 + else if(present(lb_lvl) .neqv. present(ub_lvl))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_t)) then + ndim=ndim+1 + start(ndim)=lb_t + count(ndim)=ub_t-lb_t+1 + else if(present(lb_t) .neqv. present(lb_t))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + + if ((.not. present(lb_lon)) .and. (.not. present(lb_lat)) .and. & + (.not. present(lb_lvl)) .and. (.not. present(lb_t))) then + call endrun('must specify one dimention!',subname) + endif + + !read data + call check_ret(nf_get_vara_int(ncid, varid, start(1:ndim), count(1:ndim), data), subname) + else + call endrun('the varibal does not difined!',subname) + end if + end if + if (present(readvar)) readvar = varpresent + end if + + end subroutine ncd_iolocal_int_1d + +!----------------------------------------------------------------------- +!BOP +! +! !IROUTINE: ncd_iolocal_real_1d +! +! !INTERFACE: + subroutine ncd_iolocal_real_1d(varname, data, flag, ncid, & + lb_lon, lb_lat, lb_lvl, lb_t, ub_lon, ub_lat, ub_lvl, ub_t, & + long_name, units, readvar) +! !DESCRIPTION: +! I/O for 1d int field +! +! !USES: +! +! !ARGUMENTS: + implicit none + character(len=*), intent(in) :: varname ! variable name + real, intent(inout) :: data(:) ! local decomposition data + character(len=*), intent(in) :: flag ! 'read' or 'write' + integer , intent(in) :: ncid ! input unit + integer , optional, intent(in) :: lb_lon ! start for longitude + integer , optional, intent(in) :: lb_lat ! start for latitute sizes + integer , optional, intent(in) :: lb_lvl ! start for level size + integer , optional, intent(in) :: lb_t ! start for time size + integer , optional, intent(in) :: ub_lon ! start for longitude + integer , optional, intent(in) :: ub_lat ! start for latitute sizes + integer , optional, intent(in) :: ub_lvl ! start for level size + integer , optional, intent(in) :: ub_t ! start for time size + character(len=*), optional, intent(in) :: long_name ! variable long name + character(len=*), optional, intent(in) :: units ! variable units + logical , optional, intent(out):: readvar ! true => variable is on initial dataset (read only) +! !REVISION HISTORY: +! +!EOP +! +! !LOCAL VARIABLES: + integer :: varid ! variable id + integer :: ndim ! dimension counter + integer :: start(4) ! starting indices for netcdf field + integer :: count(4) ! count values for netcdf field + character(len=32) :: inq_name ! inquid variable name + character(len=8) :: inq_xtype ! inquid variable type + integer :: inq_ndims ! inquid variable dimention + integer :: inq_dimids(4) ! inquid variable dimention id + character(len=255) :: inq_natts ! inquid variable attachment + character(len=32) :: subname='NCD_IOLOCAL_REAL_1D' ! subroutine name + logical :: varpresent ! if true, variable is on tape +!----------------------------------------------------------------------- + + ! Write field as 1d field + if (flag == 'write') then + if (nc_masterproc) then + call check_ret(nf_inq_varid(ncid, varname, varid), subname) + ! Write 1d field + ndim=0 + count=1 + if (present(lb_lon) .and. present(ub_lon)) then + ndim=ndim+1 + start(ndim)=lb_lon + count(ndim)=ub_lon-lb_lon+1 + else if(present(lb_lon) .neqv. present(ub_lon))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_lat)) then + ndim=ndim+1 + start(ndim)=lb_lat + count(ndim)=ub_lat-lb_lat+1 + else if(present(lb_lat) .neqv. present(ub_lat))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_lvl)) then + ndim=ndim+1 + start(ndim)=lb_lvl + count(ndim)=ub_lvl-lb_lvl+1 + else if(present(lb_lvl) .neqv. present(ub_lvl))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_t)) then + ndim=ndim+1 + start(ndim)=lb_t + count(ndim)=ub_t-lb_t+1 + else if(present(lb_t) .neqv. present(lb_t))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + + if ((.not. present(lb_lon)) .and. (.not. present(lb_lat)) .and. & + (.not. present(lb_lvl)) .and. (.not. present(lb_t))) then + call endrun('must specify one dimention!',subname) + endif + + call check_ret(nf_put_vara_real(ncid, varid, start(1:ndim), count(1:ndim), data), subname) + if (present(long_name)) then + call check_ret(nf_put_att_text(ncid, varid, 'long_name', len_trim(long_name), trim(long_name)), subname) + end if + if (present(units)) then + call check_ret(nf_put_att_text(ncid, varid, 'units', len_trim(units), trim(units)), subname) + end if + end if ! end of if-nc_masterproc block + ! Read field as 1d field + else if (flag == 'read') then + if (nc_masterproc) then + call check_var(ncid, varname, varid, varpresent) + if (varpresent) then + ndim=0 + count=1 + if (present(lb_lon) .and. present(ub_lon)) then + ndim=ndim+1 + start(ndim)=lb_lon + count(ndim)=ub_lon-lb_lon+1 + else if(present(lb_lon) .neqv. present(ub_lon))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_lat)) then + ndim=ndim+1 + start(ndim)=lb_lat + count(ndim)=ub_lat-lb_lat+1 + else if(present(lb_lat) .neqv. present(ub_lat))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_lvl)) then + ndim=ndim+1 + start(ndim)=lb_lvl + count(ndim)=ub_lvl-lb_lvl+1 + else if(present(lb_lvl) .neqv. present(ub_lvl))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_t)) then + ndim=ndim+1 + start(ndim)=lb_t + count(ndim)=ub_t-lb_t+1 + else if(present(lb_t) .neqv. present(lb_t))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + + if ((.not. present(lb_lon)) .and. (.not. present(lb_lat)) .and. & + (.not. present(lb_lvl)) .and. (.not. present(lb_t))) then + call endrun('must specify one dimention!',subname) + endif + + !read data + call check_ret(nf_get_vara_real(ncid, varid, start(1:ndim), count(1:ndim), data), subname) + else + call endrun('the varibal does not difined!',subname) + end if + end if + if (present(readvar)) readvar = varpresent + end if + + end subroutine ncd_iolocal_real_1d +!----------------------------------------------------------------------- +!BOP +! +! !IROUTINE: ncd_iolocal_real_1d +! +! !INTERFACE: + subroutine ncd_iolocal_double_1d(varname, data, flag, ncid, & + lb_lon, lb_lat, lb_lvl, lb_t, ub_lon, ub_lat, ub_lvl, ub_t, & + long_name, units, readvar) +! !DESCRIPTION: +! I/O for 1d int field +! +! !USES: +! +! !ARGUMENTS: + implicit none + character(len=*), intent(in) :: varname ! variable name + real*8, intent(inout) :: data(:) ! local decomposition data + character(len=*), intent(in) :: flag ! 'read' or 'write' + integer , intent(in) :: ncid ! input unit + integer , optional, intent(in) :: lb_lon ! start for longitude + integer , optional, intent(in) :: lb_lat ! start for latitute sizes + integer , optional, intent(in) :: lb_lvl ! start for level size + integer , optional, intent(in) :: lb_t ! start for time size + integer , optional, intent(in) :: ub_lon ! start for longitude + integer , optional, intent(in) :: ub_lat ! start for latitute sizes + integer , optional, intent(in) :: ub_lvl ! start for level size + integer , optional, intent(in) :: ub_t ! start for time size + character(len=*), optional, intent(in) :: long_name ! variable long name + character(len=*), optional, intent(in) :: units ! variable units + logical , optional, intent(out):: readvar ! true => variable is on initial dataset (read only) +! !REVISION HISTORY: +! +!EOP +! +! !LOCAL VARIABLES: + integer :: varid ! variable id + integer :: ndim ! dimension counter + integer :: start(4) ! starting indices for netcdf field + integer :: count(4) ! count values for netcdf field + character(len=32) :: inq_name ! inquid variable name + character(len=8) :: inq_xtype ! inquid variable type + integer :: inq_ndims ! inquid variable dimention + integer :: inq_dimids(4) ! inquid variable dimention id + character(len=255) :: inq_natts ! inquid variable attachment + character(len=32) :: subname='NCD_IOLOCAL_REAL_1D' ! subroutine name + logical :: varpresent ! if true, variable is on tape +!----------------------------------------------------------------------- + + ! Write field as 1d field + if (flag == 'write') then + if (nc_masterproc) then + call check_ret(nf_inq_varid(ncid, varname, varid), subname) + ! Write 1d field + ndim=0 + count=1 + if (present(lb_lon) .and. present(ub_lon)) then + ndim=ndim+1 + start(ndim)=lb_lon + count(ndim)=ub_lon-lb_lon+1 + else if(present(lb_lon) .neqv. present(ub_lon))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_lat)) then + ndim=ndim+1 + start(ndim)=lb_lat + count(ndim)=ub_lat-lb_lat+1 + else if(present(lb_lat) .neqv. present(ub_lat))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_lvl)) then + ndim=ndim+1 + start(ndim)=lb_lvl + count(ndim)=ub_lvl-lb_lvl+1 + else if(present(lb_lvl) .neqv. present(ub_lvl))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_t)) then + ndim=ndim+1 + start(ndim)=lb_t + count(ndim)=ub_t-lb_t+1 + else if(present(lb_t) .neqv. present(lb_t))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + + if ((.not. present(lb_lon)) .and. (.not. present(lb_lat)) .and. & + (.not. present(lb_lvl)) .and. (.not. present(lb_t))) then + call endrun('must specify one dimention!',subname) + endif + + call check_ret(nf_put_vara_double(ncid, varid, start(1:ndim), count(1:ndim), data), subname) + if (present(long_name)) then + call check_ret(nf_put_att_text(ncid, varid, 'long_name', len_trim(long_name), trim(long_name)), subname) + end if + if (present(units)) then + call check_ret(nf_put_att_text(ncid, varid, 'units', len_trim(units), trim(units)), subname) + end if + end if ! end of if-nc_masterproc block + ! Read field as 1d field + else if (flag == 'read') then + if (nc_masterproc) then + call check_var(ncid, varname, varid, varpresent) + if (varpresent) then + ndim=0 + count=1 + if (present(lb_lon) .and. present(ub_lon)) then + ndim=ndim+1 + start(ndim)=lb_lon + count(ndim)=ub_lon-lb_lon+1 + else if(present(lb_lon) .neqv. present(ub_lon))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_lat)) then + ndim=ndim+1 + start(ndim)=lb_lat + count(ndim)=ub_lat-lb_lat+1 + else if(present(lb_lat) .neqv. present(ub_lat))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_lvl)) then + ndim=ndim+1 + start(ndim)=lb_lvl + count(ndim)=ub_lvl-lb_lvl+1 + else if(present(lb_lvl) .neqv. present(ub_lvl))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_t)) then + ndim=ndim+1 + start(ndim)=lb_t + count(ndim)=ub_t-lb_t+1 + else if(present(lb_t) .neqv. present(lb_t))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + + if ((.not. present(lb_lon)) .and. (.not. present(lb_lat)) .and. & + (.not. present(lb_lvl)) .and. (.not. present(lb_t))) then + call endrun('must specify one dimention!',subname) + endif + + !read data + call check_ret(nf_get_vara_double(ncid, varid, start(1:ndim), count(1:ndim), data), subname) + else + call endrun('the varibal does not difined!',subname) + end if + end if + if (present(readvar)) readvar = varpresent + end if + + end subroutine ncd_iolocal_double_1d +!----------------------------------------------------------------------- +!BOP +! +! !IROUTINE: ncd_iolocal_int_2d +! +! !INTERFACE: + subroutine ncd_iolocal_int_2d(varname, data, flag, ncid, & + lb_lon, lb_lat, lb_lvl, lb_t, ub_lon, ub_lat, ub_lvl, ub_t, & + long_name, units, readvar) +! !DESCRIPTION: +! I/O for 2d real field +! +! !USES: +! +! !ARGUMENTS: + implicit none + character(len=*), intent(in) :: varname ! variable name + integer , intent(inout) :: data(:,:) ! local decomposition data + character(len=*), intent(in) :: flag ! 'read' or 'write' + integer , intent(in) :: ncid ! input unit + integer , optional, intent(in) :: lb_lon ! start for longitude + integer , optional, intent(in) :: lb_lat ! start for latitute sizes + integer , optional, intent(in) :: lb_lvl ! start for level size + integer , optional, intent(in) :: lb_t ! start for time size + integer , optional, intent(in) :: ub_lon ! start for longitude + integer , optional, intent(in) :: ub_lat ! start for latitute sizes + integer , optional, intent(in) :: ub_lvl ! start for level size + integer , optional, intent(in) :: ub_t ! start for time size + character(len=*), optional, intent(in) :: long_name ! variable long name + character(len=*), optional, intent(in) :: units ! variable units + logical , optional, intent(out):: readvar ! true => variable is on initial dataset (read only) +! !REVISION HISTORY: +! +!EOP +! +! !LOCAL VARIABLES: + integer :: varid ! variable id + integer :: ndim ! dimension counter + integer :: start(4) ! starting indices for netcdf field + integer :: count(4) ! count values for netcdf field + character(len=32) :: inq_name ! inquid variable name + character(len=8) :: inq_xtype ! inquid variable type + integer :: inq_ndims ! inquid variable dimention + integer :: inq_dimids(4) ! inquid variable dimention id + character(len=255) :: inq_natts ! inquid variable attachment + character(len=32) :: subname='NCD_IOLOCAL_INT_2D' ! subroutine name + logical :: varpresent ! if true, variable is on tape +!----------------------------------------------------------------------- + + ! Write field as 2d field + if (flag == 'write') then + if (nc_masterproc) then + call check_ret(nf_inq_varid(ncid, varname, varid), subname) + ! Write 2d field + ndim=0 + count=1 + if (present(lb_lon) .and. present(ub_lon)) then + ndim=ndim+1 + start(ndim)=lb_lon + count(ndim)=ub_lon-lb_lon+1 + else if(present(lb_lon) .neqv. present(ub_lon))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_lat)) then + ndim=ndim+1 + start(ndim)=lb_lat + count(ndim)=ub_lat-lb_lat+1 + else if(present(lb_lat) .neqv. present(ub_lat))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_lvl)) then + ndim=ndim+1 + start(ndim)=lb_lvl + count(ndim)=ub_lvl-lb_lvl+1 + else if(present(lb_lvl) .neqv. present(ub_lvl))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_t)) then + ndim=ndim+1 + start(ndim)=lb_t + count(ndim)=ub_t-lb_t+1 + else if(present(lb_t) .neqv. present(lb_t))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + + if ((.not. present(lb_lon)) .and. (.not. present(lb_lat)) .and. & + (.not. present(lb_lvl)) .and. (.not. present(lb_t))) then + call endrun('must specify one dimention!',subname) + endif + + call check_ret(nf_put_vara_int(ncid, varid, start(1:ndim), count(1:ndim), data), subname) + if (present(long_name)) then + call check_ret(nf_put_att_text(ncid, varid, 'long_name', len_trim(long_name), trim(long_name)), subname) + end if + if (present(units)) then + call check_ret(nf_put_att_text(ncid, varid, 'units', len_trim(units), trim(units)), subname) + end if + end if ! end of if-nc_masterproc block + ! Read field as 1d field + else if (flag == 'read') then + if (nc_masterproc) then + call check_var(ncid, varname, varid, varpresent) + if (varpresent) then + ndim=0 + count=1 + if (present(lb_lon) .and. present(ub_lon)) then + ndim=ndim+1 + start(ndim)=lb_lon + count(ndim)=ub_lon-lb_lon+1 + else if(present(lb_lon) .neqv. present(ub_lon))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_lat)) then + ndim=ndim+1 + start(ndim)=lb_lat + count(ndim)=ub_lat-lb_lat+1 + else if(present(lb_lat) .neqv. present(ub_lat))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_lvl)) then + ndim=ndim+1 + start(ndim)=lb_lvl + count(ndim)=ub_lvl-lb_lvl+1 + else if(present(lb_lvl) .neqv. present(ub_lvl))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_t)) then + ndim=ndim+1 + start(ndim)=lb_t + count(ndim)=ub_t-lb_t+1 + else if(present(lb_t) .neqv. present(lb_t))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + + if ((.not. present(lb_lon)) .and. (.not. present(lb_lat)) .and. & + (.not. present(lb_lvl)) .and. (.not. present(lb_t))) then + call endrun('must specify one dimention!',subname) + endif + + call check_ret(nf_get_vara_int(ncid, varid, start(1:ndim), count(1:ndim), data), subname) + else + call endrun('the varibal does not difined!',subname) + end if + end if + if (present(readvar)) readvar = varpresent + end if + + end subroutine ncd_iolocal_int_2d + +!----------------------------------------------------------------------- +!BOP +! +! !IROUTINE: ncd_iolocal_real_2d +! +! !INTERFACE: + subroutine ncd_iolocal_real_2d(varname, data, flag, ncid, & + lb_lon, lb_lat, lb_lvl, lb_t, ub_lon, ub_lat, ub_lvl, ub_t, & + long_name, units, readvar) +! !DESCRIPTION: +! I/O for 2d real field +! +! !USES: +! +! !ARGUMENTS: + implicit none + character(len=*), intent(in) :: varname ! variable name + real, intent(inout) :: data(:,:) ! local decomposition data + character(len=*), intent(in) :: flag ! 'read' or 'write' + integer , intent(in) :: ncid ! input unit + integer , optional, intent(in) :: lb_lon ! start for longitude + integer , optional, intent(in) :: lb_lat ! start for latitute sizes + integer , optional, intent(in) :: lb_lvl ! start for level size + integer , optional, intent(in) :: lb_t ! start for time size + integer , optional, intent(in) :: ub_lon ! start for longitude + integer , optional, intent(in) :: ub_lat ! start for latitute sizes + integer , optional, intent(in) :: ub_lvl ! start for level size + integer , optional, intent(in) :: ub_t ! start for time size + character(len=*), optional, intent(in) :: long_name ! variable long name + character(len=*), optional, intent(in) :: units ! variable units + logical , optional, intent(out):: readvar ! true => variable is on initial dataset (read only) +! !REVISION HISTORY: +! +!EOP +! +! !LOCAL VARIABLES: + integer :: varid ! variable id + integer :: ndim ! dimension counter + integer :: start(4) ! starting indices for netcdf field + integer :: count(4) ! count values for netcdf field + character(len=32) :: inq_name ! inquid variable name + character(len=8) :: inq_xtype ! inquid variable type + integer :: inq_ndims ! inquid variable dimention + integer :: inq_dimids(4) ! inquid variable dimention id + character(len=255) :: inq_natts ! inquid variable attachment + character(len=32) :: subname='NCD_IOLOCAL_REAL_2D' ! subroutine name + logical :: varpresent ! if true, variable is on tape +!----------------------------------------------------------------------- + + ! Write field as 2d field + if (flag == 'write') then + if (nc_masterproc) then + call check_ret(nf_inq_varid(ncid, varname, varid), subname) + ! Write 2d field + ndim=0 + count=1 + if (present(lb_lon) .and. present(ub_lon)) then + ndim=ndim+1 + start(ndim)=lb_lon + count(ndim)=ub_lon-lb_lon+1 + else if(present(lb_lon) .neqv. present(ub_lon))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_lat)) then + ndim=ndim+1 + start(ndim)=lb_lat + count(ndim)=ub_lat-lb_lat+1 + else if(present(lb_lat) .neqv. present(ub_lat))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_lvl)) then + ndim=ndim+1 + start(ndim)=lb_lvl + count(ndim)=ub_lvl-lb_lvl+1 + else if(present(lb_lvl) .neqv. present(ub_lvl))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_t)) then + ndim=ndim+1 + start(ndim)=lb_t + count(ndim)=ub_t-lb_t+1 + else if(present(lb_t) .neqv. present(lb_t))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + + if ((.not. present(lb_lon)) .and. (.not. present(lb_lat)) .and. & + (.not. present(lb_lvl)) .and. (.not. present(lb_t))) then + call endrun('must specify one dimention!',subname) + endif + + call check_ret(nf_put_vara_real(ncid, varid, start(1:ndim), count(1:ndim), data), subname) + if (present(long_name)) then + call check_ret(nf_put_att_text(ncid, varid, 'long_name', len_trim(long_name), trim(long_name)), subname) + end if + if (present(units)) then + call check_ret(nf_put_att_text(ncid, varid, 'units', len_trim(units), trim(units)), subname) + end if + end if ! end of if-nc_masterproc block + ! Read field as 1d field + else if (flag == 'read') then + if (nc_masterproc) then + call check_var(ncid, varname, varid, varpresent) + if (varpresent) then + ndim=0 + count=1 + if (present(lb_lon) .and. present(ub_lon)) then + ndim=ndim+1 + start(ndim)=lb_lon + count(ndim)=ub_lon-lb_lon+1 + else if(present(lb_lon) .neqv. present(ub_lon))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_lat)) then + ndim=ndim+1 + start(ndim)=lb_lat + count(ndim)=ub_lat-lb_lat+1 + else if(present(lb_lat) .neqv. present(ub_lat))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_lvl)) then + ndim=ndim+1 + start(ndim)=lb_lvl + count(ndim)=ub_lvl-lb_lvl+1 + else if(present(lb_lvl) .neqv. present(ub_lvl))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_t)) then + ndim=ndim+1 + start(ndim)=lb_t + count(ndim)=ub_t-lb_t+1 + else if(present(lb_t) .neqv. present(lb_t))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + + if ((.not. present(lb_lon)) .and. (.not. present(lb_lat)) .and. & + (.not. present(lb_lvl)) .and. (.not. present(lb_t))) then + call endrun('must specify one dimention!',subname) + endif + + call check_ret(nf_get_vara_real(ncid, varid, start(1:ndim), count(1:ndim), data), subname) + else + call endrun('the varibal does not difined!',subname) + end if + end if + if (present(readvar)) readvar = varpresent + end if + + end subroutine ncd_iolocal_real_2d + + +!----------------------------------------------------------------------- +!BOP +! +! !IROUTINE: ncd_iolocal_real_2d +! +! !INTERFACE: + subroutine ncd_iolocal_double_2d(varname, data, flag, ncid, & + lb_lon, lb_lat, lb_lvl, lb_t, ub_lon, ub_lat, ub_lvl, ub_t, & + long_name, units, readvar) +! !DESCRIPTION: +! I/O for 2d real field +! +! !USES: +! +! !ARGUMENTS: + implicit none + character(len=*), intent(in) :: varname ! variable name + real*8, intent(inout) :: data(:,:) ! local decomposition data + character(len=*), intent(in) :: flag ! 'read' or 'write' + integer , intent(in) :: ncid ! input unit + integer , optional, intent(in) :: lb_lon ! start for longitude + integer , optional, intent(in) :: lb_lat ! start for latitute sizes + integer , optional, intent(in) :: lb_lvl ! start for level size + integer , optional, intent(in) :: lb_t ! start for time size + integer , optional, intent(in) :: ub_lon ! start for longitude + integer , optional, intent(in) :: ub_lat ! start for latitute sizes + integer , optional, intent(in) :: ub_lvl ! start for level size + integer , optional, intent(in) :: ub_t ! start for time size + character(len=*), optional, intent(in) :: long_name ! variable long name + character(len=*), optional, intent(in) :: units ! variable units + logical , optional, intent(out):: readvar ! true => variable is on initial dataset (read only) +! !REVISION HISTORY: +! +!EOP +! +! !LOCAL VARIABLES: + integer :: varid ! variable id + integer :: ndim ! dimension counter + integer :: start(4) ! starting indices for netcdf field + integer :: count(4) ! count values for netcdf field + character(len=32) :: inq_name ! inquid variable name + character(len=8) :: inq_xtype ! inquid variable type + integer :: inq_ndims ! inquid variable dimention + integer :: inq_dimids(4) ! inquid variable dimention id + character(len=255) :: inq_natts ! inquid variable attachment + character(len=32) :: subname='NCD_IOLOCAL_REAL_2D' ! subroutine name + logical :: varpresent ! if true, variable is on tape +!----------------------------------------------------------------------- + + ! Write field as 2d field + if (flag == 'write') then + if (nc_masterproc) then + call check_ret(nf_inq_varid(ncid, varname, varid), subname) + ! Write 2d field + ndim=0 + count=1 + if (present(lb_lon) .and. present(ub_lon)) then + ndim=ndim+1 + start(ndim)=lb_lon + count(ndim)=ub_lon-lb_lon+1 + else if(present(lb_lon) .neqv. present(ub_lon))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_lat)) then + ndim=ndim+1 + start(ndim)=lb_lat + count(ndim)=ub_lat-lb_lat+1 + else if(present(lb_lat) .neqv. present(ub_lat))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_lvl)) then + ndim=ndim+1 + start(ndim)=lb_lvl + count(ndim)=ub_lvl-lb_lvl+1 + else if(present(lb_lvl) .neqv. present(ub_lvl))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_t)) then + ndim=ndim+1 + start(ndim)=lb_t + count(ndim)=ub_t-lb_t+1 + else if(present(lb_t) .neqv. present(lb_t))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + + if ((.not. present(lb_lon)) .and. (.not. present(lb_lat)) .and. & + (.not. present(lb_lvl)) .and. (.not. present(lb_t))) then + call endrun('must specify one dimention!',subname) + endif + + call check_ret(nf_put_vara_double(ncid, varid, start(1:ndim), count(1:ndim), data), subname) + if (present(long_name)) then + call check_ret(nf_put_att_text(ncid, varid, 'long_name', len_trim(long_name), trim(long_name)), subname) + end if + if (present(units)) then + call check_ret(nf_put_att_text(ncid, varid, 'units', len_trim(units), trim(units)), subname) + end if + end if ! end of if-nc_masterproc block + ! Read field as 1d field + else if (flag == 'read') then + if (nc_masterproc) then + call check_var(ncid, varname, varid, varpresent) + if (varpresent) then + ndim=0 + count=1 + if (present(lb_lon) .and. present(ub_lon)) then + ndim=ndim+1 + start(ndim)=lb_lon + count(ndim)=ub_lon-lb_lon+1 + else if(present(lb_lon) .neqv. present(ub_lon))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_lat)) then + ndim=ndim+1 + start(ndim)=lb_lat + count(ndim)=ub_lat-lb_lat+1 + else if(present(lb_lat) .neqv. present(ub_lat))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_lvl)) then + ndim=ndim+1 + start(ndim)=lb_lvl + count(ndim)=ub_lvl-lb_lvl+1 + else if(present(lb_lvl) .neqv. present(ub_lvl))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + if (present(lb_t)) then + ndim=ndim+1 + start(ndim)=lb_t + count(ndim)=ub_t-lb_t+1 + else if(present(lb_t) .neqv. present(lb_t))then + call endrun('must specify the low and up boundary at the same time!',subname) + endif + + if ((.not. present(lb_lon)) .and. (.not. present(lb_lat)) .and. & + (.not. present(lb_lvl)) .and. (.not. present(lb_t))) then + call endrun('must specify one dimention!',subname) + endif + + call check_ret(nf_get_vara_double(ncid, varid, start(1:ndim), count(1:ndim), data), subname) + else + call endrun('the varibal does not difined!',subname) + end if + end if + if (present(readvar)) readvar = varpresent + end if + + end subroutine ncd_iolocal_double_2d + + +!----------------------------------------------------------------------- +!BOP +! +! !IROUTINE: ncd_ioglobal_int_var +! +! !INTERFACE: + subroutine ncd_ioglobal_int_var(varname, data, flag, ncid, long_name, units, nt, readvar) +! !DESCRIPTION: +! I/O of integer variable +! + +! !ARGUMENTS: + implicit none + character(len=*), intent(in) :: varname ! variable name + integer , intent(inout) :: data ! local decomposition data + character(len=*), intent(in) :: flag ! 'read' or 'write' + integer , intent(in) :: ncid ! input unit + character(len=*), optional, intent(in) :: long_name ! variable long name + character(len=*), optional, intent(in) :: units ! variable units + integer , optional, intent(in) :: nt ! time sample index + logical , optional, intent(out):: readvar ! true => variable is on initial dataset (read only) +! +! !REVISION HISTORY: +! +!EOP +! +! !LOCAL VARIABLES: + integer :: ier ! error status + integer :: dimid(1) ! dimension id + integer :: start(1), count(1) ! output bounds + integer :: varid ! variable id + logical :: varpresent ! if true, variable is on tape + character(len=32) :: subname='NCD_IOGLOBAL_INT_VAR' ! subroutine name +!----------------------------------------------------------------------- + + if (flag == 'write') then + + if (nc_masterproc) then + call check_ret(nf_inq_varid(ncid, varname, varid), subname) + if (present(nt)) then + start(1) = nt; count(1) = 1 + call check_ret(nf_put_vara_int(ncid, varid, start, count, data), subname) + else + call check_ret(nf_put_var_int(ncid, varid, data), subname) + end if + if (present(long_name)) then + call check_ret(nf_put_att_text(ncid, varid, 'long_name', len_trim(long_name), trim(long_name)), subname) + end if + if (present(units)) then + call check_ret(nf_put_att_text(ncid, varid, 'units', len_trim(units), trim(units)), subname) + end if + end if + + else if (flag == 'read') then + + if (nc_masterproc) then + call check_var(ncid, varname, varid, varpresent) + if (varpresent) then + call check_ret(nf_get_var_int(ncid, varid, data), subname) + else + call endrun('the varibal does not difined!',subname) + endif + end if + if (present(readvar)) readvar = varpresent + + end if + + end subroutine ncd_ioglobal_int_var + +!----------------------------------------------------------------------- +!BOP +! +! !IROUTINE: ncd_ioglobal_real_var +! +! !INTERFACE: + subroutine ncd_ioglobal_real_var(varname, data, flag, ncid, long_name, units, nt, readvar) +! !DESCRIPTION: +! I/O of real variable +! + +! !ARGUMENTS: + implicit none + character(len=*), intent(in) :: varname ! variable name + real , intent(inout) :: data ! local decomposition data + character(len=*), intent(in) :: flag ! 'read' or 'write' + integer , intent(in) :: ncid ! input unit + character(len=*), optional, intent(in) :: long_name ! variable long name + character(len=*), optional, intent(in) :: units ! variable units + logical , optional, intent(out):: readvar ! true => variable is on initial dataset (read only) + integer , optional, intent(in) :: nt ! time sample index +! +! !REVISION HISTORY: +! +!EOP +! +! !LOCAL VARIABLES: + integer :: ier ! error status + integer :: dimid(1) ! dimension id + integer :: start(1), count(1) ! output bounds + integer :: varid ! variable id + logical :: varpresent ! if true, variable is on tape + character(len=32) :: subname='NCD_IOGLOBAL_REAL_VAR' ! subroutine name +!----------------------------------------------------------------------- + + if (flag == 'write') then + + if (nc_masterproc) then + call check_ret(nf_inq_varid(ncid, varname, varid), subname) + if (present(nt)) then + start(1) = nt; count(1) = 1 + call check_ret(nf_put_vara_real(ncid, varid, start, count, data), subname) + else + call check_ret(nf_put_var_real(ncid, varid, data), subname) + end if + if (present(long_name)) then + call check_ret(nf_put_att_text(ncid, varid, 'long_name', len_trim(long_name), trim(long_name)), subname) + end if + if (present(units)) then + call check_ret(nf_put_att_text(ncid, varid, 'units', len_trim(units), trim(units)), subname) + end if + end if + + else if (flag == 'read') then + + if (nc_masterproc) then + call check_var(ncid, varname, varid, varpresent) + if (varpresent) then + call check_ret(nf_get_var_real(ncid, varid, data), subname) + else + call endrun('the varibal does not difined!',subname) + endif + end if + if (present(readvar)) readvar = varpresent + + end if + + end subroutine ncd_ioglobal_real_var + +!----------------------------------------------------------------------- +!BOP +! +! !IROUTINE: ncd_ioglobal_real_var +! +! !INTERFACE: + subroutine ncd_ioglobal_double_var(varname, data, flag, ncid, long_name, units, nt, readvar) +! !DESCRIPTION: +! I/O of real variable +! + +! !ARGUMENTS: + implicit none + character(len=*), intent(in) :: varname ! variable name + real*8 , intent(inout) :: data ! local decomposition data + character(len=*), intent(in) :: flag ! 'read' or 'write' + integer , intent(in) :: ncid ! input unit + character(len=*), optional, intent(in) :: long_name ! variable long name + character(len=*), optional, intent(in) :: units ! variable units + logical , optional, intent(out):: readvar ! true => variable is on initial dataset (read only) + integer , optional, intent(in) :: nt ! time sample index +! +! !REVISION HISTORY: +! +!EOP +! +! !LOCAL VARIABLES: + integer :: ier ! error status + integer :: dimid(1) ! dimension id + integer :: start(1), count(1) ! output bounds + integer :: varid ! variable id + logical :: varpresent ! if true, variable is on tape + character(len=32) :: subname='NCD_IOGLOBAL_REAL_VAR' ! subroutine name +!----------------------------------------------------------------------- + + if (flag == 'write') then + + if (nc_masterproc) then + call check_ret(nf_inq_varid(ncid, varname, varid), subname) + if (present(nt)) then + start(1) = nt; count(1) = 1 + call check_ret(nf_put_vara_double(ncid, varid, start, count, data), subname) + else + call check_ret(nf_put_var_double(ncid, varid, data), subname) + end if + if (present(long_name)) then + call check_ret(nf_put_att_text(ncid, varid, 'long_name', len_trim(long_name), trim(long_name)), subname) + end if + if (present(units)) then + call check_ret(nf_put_att_text(ncid, varid, 'units', len_trim(units), trim(units)), subname) + end if + end if + + else if (flag == 'read') then + + if (nc_masterproc) then + call check_var(ncid, varname, varid, varpresent) + if (varpresent) then + call check_ret(nf_get_var_double(ncid, varid, data), subname) + else + call endrun('the varibal does not difined!',subname) + endif + end if + if (present(readvar)) readvar = varpresent + + end if + + end subroutine ncd_ioglobal_double_var + +!---------------------------------------------------------------------------------- +!BOP +! +! !IROUTINE: ncd_ioglobal_int_1d +! +! !INTERFACE: + subroutine ncd_ioglobal_int_1d(varname, data, flag, ncid, long_name, units, nt, readvar) +! !DESCRIPTION: +! Master I/O for 1d integer data +! +! !ARGUMENTS: + implicit none + character(len=*), intent(in) :: flag ! 'read' or 'write' + integer , intent(in) :: ncid ! input unit + character(len=*), intent(in) :: varname ! variable name + integer , intent(inout) :: data(:) ! local decomposition data + character(len=*), optional, intent(in) :: long_name ! variable long name + character(len=*), optional, intent(in) :: units ! variable units + logical , optional, intent(out):: readvar ! true => variable is on initial dataset (read only) + integer , optional, intent(in) :: nt ! time sample index +! +! !REVISION HISTORY: +! +!EOP +! +! !LOCAL VARIABLES: + integer :: varid ! netCDF variable id + integer :: dimid(2), ndims ! dimension ids + integer :: start(2), count(2) ! output bounds + integer :: ier ! error code + logical :: varpresent ! if true, variable is on tape + character(len=32) :: subname='NCD_IOGLOBAL_INT_1D' ! subroutine name +!----------------------------------------------------------------------- + + if (flag == 'write') then + + if (nc_masterproc) then + call check_ret(nf_inq_varid(ncid, varname, varid), subname) + if (present(nt)) then + start(1) = 1; count(1) = size(data) + start(2) = nt; count(2) = 1 + call check_ret(nf_put_vara_int(ncid, varid, start, count, data), subname) + else + call check_ret(nf_put_var_int(ncid, varid, data), subname) + end if + if (present(long_name)) then + call check_ret(nf_put_att_text(ncid, varid, 'long_name', len_trim(long_name), trim(long_name)), subname) + end if + if (present(units)) then + call check_ret(nf_put_att_text(ncid, varid, 'units', len_trim(units), trim(units)), subname) + end if + end if + + else if (flag == 'read') then + + if (nc_masterproc) then + call check_var(ncid, varname, varid, varpresent) + if (varpresent) then + call check_ret(nf_get_var_int(ncid, varid, data), subname) + else + call endrun('the varibal does not difined!',subname) + endif + end if + if (present(readvar)) readvar = varpresent + + end if + + end subroutine ncd_ioglobal_int_1d + +!----------------------------------------------------------------------- +!BOP +! +! !IROUTINE: ncd_ioglobal_real_1d +! +! !INTERFACE: + subroutine ncd_ioglobal_real_1d(varname, data, flag, ncid, long_name, units, nt, readvar) +! !DESCRIPTION: +! Master I/O for 1d real data +! +! !ARGUMENTS: + implicit none + character(len=*), intent(in) :: flag ! 'read' or 'write' + integer , intent(in) :: ncid ! input unit + character(len=*), intent(in) :: varname ! variable name + character(len=*), optional, intent(in) :: long_name ! variable long name + character(len=*), optional, intent(in) :: units ! variable units + real , intent(inout) :: data(:) ! local decomposition input data + logical , optional, intent(out):: readvar ! true => variable is on initial dataset (read only) + integer , optional, intent(in) :: nt ! time sample index +! +! !REVISION HISTORY: +! +!EOP +! +! !LOCAL VARIABLES: + integer :: varid ! netCDF variable id + integer :: ier ! error code + integer :: dimid(2), ndims ! dimension ids + integer :: start(2), count(2) ! output bounds + logical :: varpresent ! if true, variable is on tape + character(len=32) :: subname='NCD_IOGLOBAL_REAL_1D' ! subroutine name +!----------------------------------------------------------------------- + + if (flag == 'write') then + + if (nc_masterproc) then + call check_ret(nf_inq_varid(ncid, varname, varid), subname) + if (present(nt)) then + start(1) = 1; count(1) = size(data) + start(2) = nt; count(2) = 1 + call check_ret(nf_put_vara_real(ncid, varid, start, count, data), subname) + else +! call check_ret(nf_put_var_real(ncid, varid, data), subname) +call check_ret(nf_put_var_real(ncid, varid, data), subname) + end if + if (present(long_name)) then + call check_ret(nf_put_att_text(ncid, varid, 'long_name', len_trim(long_name), trim(long_name)), subname) + end if + if (present(units)) then + call check_ret(nf_put_att_text(ncid, varid, 'units', len_trim(units), trim(units)), subname) + end if + end if + + else if (flag == 'read') then + + if (nc_masterproc) then + call check_var(ncid, varname, varid, varpresent) + if (varpresent) then + call check_ret(nf_get_var_real(ncid, varid, data), subname) + else + call endrun('the varibal does not difined!',subname) + endif + end if + if (present(readvar)) readvar = varpresent + + end if + + end subroutine ncd_ioglobal_real_1d + +!----------------------------------------------------------------------- +!BOP +! +! !IROUTINE: ncd_ioglobal_real_1d +! +! !INTERFACE: + subroutine ncd_ioglobal_double_1d(varname, data, flag, ncid, long_name, units, nt, readvar) +! !DESCRIPTION: +! Master I/O for 1d real data +! +! !ARGUMENTS: + implicit none + character(len=*), intent(in) :: flag ! 'read' or 'write' + integer , intent(in) :: ncid ! input unit + character(len=*), intent(in) :: varname ! variable name + character(len=*), optional, intent(in) :: long_name ! variable long name + character(len=*), optional, intent(in) :: units ! variable units + real*8 , intent(inout) :: data(:) ! local decomposition input data + logical , optional, intent(out):: readvar ! true => variable is on initial dataset (read only) + integer , optional, intent(in) :: nt ! time sample index +! +! !REVISION HISTORY: +! +!EOP +! +! !LOCAL VARIABLES: + integer :: varid ! netCDF variable id + integer :: ier ! error code + integer :: dimid(2), ndims ! dimension ids + integer :: start(2), count(2) ! output bounds + logical :: varpresent ! if true, variable is on tape + character(len=32) :: subname='NCD_IOGLOBAL_REAL_1D' ! subroutine name +!----------------------------------------------------------------------- + + if (flag == 'write') then + + if (nc_masterproc) then + call check_ret(nf_inq_varid(ncid, varname, varid), subname) + if (present(nt)) then + start(1) = 1; count(1) = size(data) + start(2) = nt; count(2) = 1 + call check_ret(nf_put_vara_double(ncid, varid, start, count, data), subname) + else +! call check_ret(nf_put_var_double(ncid, varid, data), subname) +call check_ret(nf_put_var_double(ncid, varid, data), subname) + end if + if (present(long_name)) then + call check_ret(nf_put_att_text(ncid, varid, 'long_name', len_trim(long_name), trim(long_name)), subname) + end if + if (present(units)) then + call check_ret(nf_put_att_text(ncid, varid, 'units', len_trim(units), trim(units)), subname) + end if + end if + + else if (flag == 'read') then + + if (nc_masterproc) then + call check_var(ncid, varname, varid, varpresent) + if (varpresent) then + call check_ret(nf_get_var_double(ncid, varid, data), subname) + else + call endrun('the varibal does not difined!',subname) + endif + end if + if (present(readvar)) readvar = varpresent + + end if + + end subroutine ncd_ioglobal_double_1d + +!----------------------------------------------------------------------- +!BOP +! +! !IROUTINE: ncd_ioglobal_int_2d +! +! !INTERFACE: + subroutine ncd_ioglobal_int_2d(varname, data, flag, ncid, long_name, units, nt, readvar) +! !DESCRIPTION: +! netcdf I/O of global 2d integer array +! +! !ARGUMENTS: + implicit none + character(len=*), intent(in) :: flag ! 'read' or 'write' + integer , intent(in) :: ncid ! input unit + character(len=*), intent(in) :: varname ! variable name + integer , intent(inout) :: data(:,:) ! local decomposition input data + character(len=*), optional, intent(in) :: long_name ! variable long name + character(len=*), optional, intent(in) :: units ! variable units + logical , optional, intent(out):: readvar ! true => variable is on initial dataset (read only) + integer , optional, intent(in) :: nt ! time sample index +! +! !REVISION HISTORY: +! +!EOP +! +! !LOCAL VARIABLES: + integer :: varid ! netCDF variable id + integer :: dimid(3), ndims ! dimension ids + integer :: start(3), count(3) ! output bounds + integer :: ier ! error code + logical :: varpresent ! if true, variable is on tape + character(len=32) :: subname='NCD_IOGLOBAL_2D_INT_IO' ! subroutine name +!----------------------------------------------------------------------- + + if (flag == 'write') then + + if (nc_masterproc) then + call check_ret(nf_inq_varid(ncid, varname, varid), subname) + if (present(nt)) then + start(1) = 1; count(1) = size(data, dim=1) + start(2) = 1; count(2) = size(data, dim=2) + start(3) = nt; count(3) = 1 + call check_ret(nf_put_vara_int(ncid, varid, start, count, data), subname) + else + call check_ret(nf_put_var_int(ncid, varid, data), subname) + end if + if (present(long_name)) then + call check_ret(nf_put_att_text(ncid, varid, 'long_name', len_trim(long_name), trim(long_name)), subname) + end if + if (present(units)) then + call check_ret(nf_put_att_text(ncid, varid, 'units', len_trim(units), trim(units)), subname) + end if + end if + + else if (flag == 'read') then + + if (nc_masterproc) then + call check_var(ncid, varname, varid, varpresent) + if (varpresent) then + call check_ret(nf_get_var_int(ncid, varid, data), subname) + else + call endrun('the varibal does not difined!',subname) + endif + end if + if (present(readvar)) readvar = varpresent + + end if + + end subroutine ncd_ioglobal_int_2d + +!----------------------------------------------------------------------- + +!BOP +! +! !IROUTINE: ncd_ioglobal_int_2d +! +! !INTERFACE: + subroutine ncd_ioglobal_long_2d(varname, data, flag, ncid, long_name, units, nt, readvar) +! !DESCRIPTION: +! netcdf I/O of global 2d integer array +! +! !ARGUMENTS: + implicit none + character(len=*), intent(in) :: flag ! 'read' or 'write' + integer , intent(in) :: ncid ! input unit + character(len=*), intent(in) :: varname ! variable name + integer*8 , intent(inout) :: data(:,:) ! local decomposition input data + character(len=*), optional, intent(in) :: long_name ! variable long name + character(len=*), optional, intent(in) :: units ! variable units + logical , optional, intent(out):: readvar ! true => variable is on initial dataset (read only) + integer , optional, intent(in) :: nt ! time sample index +! +! !REVISION HISTORY: +! +!EOP +! +! !LOCAL VARIABLES: + integer :: varid ! netCDF variable id + integer :: dimid(3), ndims ! dimension ids + integer :: start(3), count(3) ! output bounds + integer :: ier ! error code + logical :: varpresent ! if true, variable is on tape + character(len=32) :: subname='NCD_IOGLOBAL_2D_INT_IO' ! subroutine name +!----------------------------------------------------------------------- + + if (flag == 'write') then + + if (nc_masterproc) then + call check_ret(nf_inq_varid(ncid, varname, varid), subname) + if (present(nt)) then + start(1) = 1; count(1) = size(data, dim=1) + start(2) = 1; count(2) = size(data, dim=2) + start(3) = nt; count(3) = 1 + call check_ret(nf_put_vara_int(ncid, varid, start, count, data), subname) + else + call check_ret(nf_put_var_int(ncid, varid, data), subname) + end if + if (present(long_name)) then + call check_ret(nf_put_att_text(ncid, varid, 'long_name', len_trim(long_name), trim(long_name)), subname) + end if + if (present(units)) then + call check_ret(nf_put_att_text(ncid, varid, 'units', len_trim(units), trim(units)), subname) + end if + end if + + else if (flag == 'read') then + + if (nc_masterproc) then + call check_var(ncid, varname, varid, varpresent) + if (varpresent) then + call check_ret(nf_get_var_int(ncid, varid, data), subname) + else + call endrun('the varibal does not difined!',subname) + endif + end if + if (present(readvar)) readvar = varpresent + + end if + + end subroutine ncd_ioglobal_long_2d + +!----------------------------------------------------------------------- + +!BOP +! +! !IROUTINE: ncd_ioglobal_byte_2d +! +! !INTERFACE: + subroutine ncd_ioglobal_byte_2d(varname, data, flag, ncid, long_name, units, nt, readvar) +! !DESCRIPTION: +! netcdf I/O of global 2d integer array +! +! !ARGUMENTS: + implicit none + character(len=*), intent(in) :: flag ! 'read' or 'write' + integer , intent(in) :: ncid ! input unit + character(len=*), intent(in) :: varname ! variable name + byte, intent(inout) :: data(:,:) ! local decomposition input data + character(len=*), optional, intent(in) :: long_name ! variable long name + character(len=*), optional, intent(in) :: units ! variable units + logical , optional, intent(out):: readvar ! true => variable is on initial dataset (read only) + integer , optional, intent(in) :: nt ! time sample index +! +! !REVISION HISTORY: +! +!EOP +! +! !LOCAL VARIABLES: + integer :: varid ! netCDF variable id + integer :: dimid(3), ndims ! dimension ids + integer :: start(3), count(3) ! output bounds + integer :: ier ! error code + logical :: varpresent ! if true, variable is on tape + character(len=32) :: subname='NCD_IOGLOBAL_2D_INT1_IO' ! subroutine name +!----------------------------------------------------------------------- + + if (flag == 'write') then + + if (nc_masterproc) then + call check_ret(nf_inq_varid(ncid, varname, varid), subname) + if (present(nt)) then + start(1) = 1; count(1) = size(data, dim=1) + start(2) = 1; count(2) = size(data, dim=2) + start(3) = nt; count(3) = 1 + call check_ret(nf_put_vara_int1(ncid, varid, start, count, data), subname) + else + call check_ret(nf_put_var_int1(ncid, varid, data), subname) + end if + if (present(long_name)) then + call check_ret(nf_put_att_text(ncid, varid, 'long_name', len_trim(long_name), trim(long_name)), subname) + end if + if (present(units)) then + call check_ret(nf_put_att_text(ncid, varid, 'units', len_trim(units), trim(units)), subname) + end if + end if + + else if (flag == 'read') then + + if (nc_masterproc) then + call check_var(ncid, varname, varid, varpresent) + if (varpresent) then + call check_ret(nf_get_var_int1(ncid, varid, data), subname) + else + call endrun('the varibal does not difined!',subname) + endif + end if + if (present(readvar)) readvar = varpresent + + end if + + end subroutine ncd_ioglobal_byte_2d +!----------------------------------------------------------------------- +!BOP +! +! !IROUTINE: ncd_ioglobal_short_2d +! +! !INTERFACE: + subroutine ncd_ioglobal_short_2d(varname, data, flag, ncid, long_name, units, nt, readvar) +! !DESCRIPTION: +! netcdf I/O of global 2d integer array +! +! !ARGUMENTS: + implicit none + character(len=*), intent(in) :: flag ! 'read' or 'write' + integer , intent(in) :: ncid ! input unit + character(len=*), intent(in) :: varname ! variable name + integer*2, intent(inout) :: data(:,:) ! local decomposition input data + character(len=*), optional, intent(in) :: long_name ! variable long name + character(len=*), optional, intent(in) :: units ! variable units + logical , optional, intent(out):: readvar ! true => variable is on initial dataset (read only) + integer , optional, intent(in) :: nt ! time sample index +! +! !REVISION HISTORY: +! +!EOP +! +! !LOCAL VARIABLES: + integer :: varid ! netCDF variable id + integer :: dimid(3), ndims ! dimension ids + integer :: start(3), count(3) ! output bounds + integer :: ier ! error code + logical :: varpresent ! if true, variable is on tape + character(len=32) :: subname='NCD_IOGLOBAL_2D_INT2_IO' ! subroutine name +!----------------------------------------------------------------------- + + if (flag == 'write') then + + if (nc_masterproc) then + call check_ret(nf_inq_varid(ncid, varname, varid), subname) + if (present(nt)) then + start(1) = 1; count(1) = size(data, dim=1) + start(2) = 1; count(2) = size(data, dim=2) + start(3) = nt; count(3) = 1 + call check_ret(nf_put_vara_int2(ncid, varid, start, count, data), subname) + else + call check_ret(nf_put_var_int2(ncid, varid, data), subname) + end if + if (present(long_name)) then + call check_ret(nf_put_att_text(ncid, varid, 'long_name', len_trim(long_name), trim(long_name)), subname) + end if + if (present(units)) then + call check_ret(nf_put_att_text(ncid, varid, 'units', len_trim(units), trim(units)), subname) + end if + end if + + else if (flag == 'read') then + + if (nc_masterproc) then + call check_var(ncid, varname, varid, varpresent) + if (varpresent) then + call check_ret(nf_get_var_int2(ncid, varid, data), subname) + else + call endrun('the varibal does not difined!',subname) + endif + end if + if (present(readvar)) readvar = varpresent + + end if + + end subroutine ncd_ioglobal_short_2d +!----------------------------------------------------------------------- +!BOP +! +! !IROUTINE: ncd_ioglobal_real_2d +! +! !INTERFACE: + subroutine ncd_ioglobal_real_2d(varname, data, flag, & + ncid, long_name, units, nt, readvar) +! !DESCRIPTION: +! netcdf I/O of global 2d real array +! +! !ARGUMENTS: + implicit none + character(len=*), intent(in) :: flag ! 'read' or 'write' + integer , intent(in) :: ncid ! input unit + character(len=*), intent(in) :: varname ! variable name + real , intent(inout) :: data(:,:) ! local decomposition input data + character(len=*), optional, intent(in) :: long_name ! variable long name + character(len=*), optional, intent(in) :: units ! variable units + logical , optional, intent(out):: readvar ! true => variable is on initial dataset (read only) + integer , optional, intent(in) :: nt ! time sample index +! +! !REVISION HISTORY: +! +!EOP +! +! !LOCAL VARIABLES: + integer :: varid ! netCDF variable id + integer :: ier ! error code + integer :: dimid(3), ndims ! dimension ids + integer :: start(3), count(3) ! output bounds + logical :: varpresent ! if true, variable is on tape + character(len=32) :: subname='NCD_IOGLOBAL_REAL_2D' ! subroutine name +!----------------------------------------------------------------------- + + if (flag == 'write') then + + if (nc_masterproc) then + call check_ret(nf_inq_varid(ncid, varname, varid), subname) + if (present(nt)) then + start(1) = 1; count(1) = size(data, dim=1) + start(2) = 1; count(2) = size(data, dim=2) + start(3) = nt; count(3) = 1 +! call check_ret(nf_put_vara_real(ncid, varid, start, count, data), subname) +call check_ret(nf_put_vara_real(ncid, varid, start, count, data), subname) + else + call check_ret(nf_put_var_real(ncid, varid, data), subname) + end if + if (present(long_name)) then + call check_ret(nf_put_att_text(ncid, varid, 'long_name', len_trim(long_name), trim(long_name)), subname) + end if + if (present(units)) then + call check_ret(nf_put_att_text(ncid, varid, 'units', len_trim(units), trim(units)), subname) + end if + end if + + else if (flag == 'read') then + + if (nc_masterproc) then + call check_var(ncid, varname, varid, varpresent) + if (varpresent) then + call check_ret(nf_get_var_real(ncid, varid, data), subname) + else + call endrun('the varibal does not difined!',subname) + endif + end if + if (present(readvar)) readvar = varpresent + + end if + + end subroutine ncd_ioglobal_real_2d + +!----------------------------------------------------------------------- +!BOP +! +! !IROUTINE: ncd_ioglobal_real_2d +! +! !INTERFACE: + subroutine ncd_ioglobal_double_2d(varname, data, flag, & + ncid, long_name, units, nt, readvar) +! !DESCRIPTION: +! netcdf I/O of global 2d real array +! +! !ARGUMENTS: + implicit none + character(len=*), intent(in) :: flag ! 'read' or 'write' + integer , intent(in) :: ncid ! input unit + character(len=*), intent(in) :: varname ! variable name + real*8 , intent(inout) :: data(:,:) ! local decomposition input data + character(len=*), optional, intent(in) :: long_name ! variable long name + character(len=*), optional, intent(in) :: units ! variable units + logical , optional, intent(out):: readvar ! true => variable is on initial dataset (read only) + integer , optional, intent(in) :: nt ! time sample index +! +! !REVISION HISTORY: +! +!EOP +! +! !LOCAL VARIABLES: + integer :: varid ! netCDF variable id + integer :: ier ! error code + integer :: dimid(3), ndims ! dimension ids + integer :: start(3), count(3) ! output bounds + logical :: varpresent ! if true, variable is on tape + character(len=32) :: subname='NCD_IOGLOBAL_REAL_2D' ! subroutine name +!----------------------------------------------------------------------- + + if (flag == 'write') then + + if (nc_masterproc) then + call check_ret(nf_inq_varid(ncid, varname, varid), subname) + if (present(nt)) then + start(1) = 1; count(1) = size(data, dim=1) + start(2) = 1; count(2) = size(data, dim=2) + start(3) = nt; count(3) = 1 +! call check_ret(nf_put_vara_double(ncid, varid, start, count, data), subname) +call check_ret(nf_put_vara_double(ncid, varid, start, count, data), subname) + else + call check_ret(nf_put_var_double(ncid, varid, data), subname) + end if + if (present(long_name)) then + call check_ret(nf_put_att_text(ncid, varid, 'long_name', len_trim(long_name), trim(long_name)), subname) + end if + if (present(units)) then + call check_ret(nf_put_att_text(ncid, varid, 'units', len_trim(units), trim(units)), subname) + end if + end if + + else if (flag == 'read') then + + if (nc_masterproc) then + call check_var(ncid, varname, varid, varpresent) + if (varpresent) then + call check_ret(nf_get_var_double(ncid, varid, data), subname) + else + call endrun('the varibal does not difined!',subname) + endif + end if + if (present(readvar)) readvar = varpresent + + end if + + end subroutine ncd_ioglobal_double_2d +!----------------------------------------------------------------------- +!BOP +! +! !IROUTINE: ncd_ioglobal_short_3d +! +! !INTERFACE: + subroutine ncd_ioglobal_short_3d(varname, data, flag, & + ncid, long_name, units, nt, readvar) +! !DESCRIPTION: +! netcdf I/O of global 3d integer array +! +! !ARGUMENTS: + implicit none + character(len=*), intent(in) :: flag ! 'read' or 'write' + integer , intent(in) :: ncid ! input unit + character(len=*), intent(in) :: varname ! variable name + integer*2 , intent(inout) :: data(:,:,:) ! local decomposition input data + character(len=*), optional, intent(in) :: long_name ! variable long name + character(len=*), optional, intent(in) :: units ! variable units + logical , optional, intent(out):: readvar ! true => variable is on initial dataset (read only) + integer , optional, intent(in) :: nt ! time sample index +! +! !REVISION HISTORY: +! +!EOP +! +! !LOCAL VARIABLES: + integer :: varid ! netCDF variable id + integer :: dimid(4), ndims ! dimension ids + integer :: start(4), count(4) ! output bounds + integer :: ier ! error code + logical :: varpresent ! if true, variable is on tape + character(len=32) :: subname='NCD_IOGLOBAL_3D_INT2_IO' ! subroutine name +!----------------------------------------------------------------------- + + if (flag == 'write') then + + if (nc_masterproc) then + call check_ret(nf_inq_varid(ncid, varname, varid), subname) + if (present(nt)) then + start(1) = 1; count(1) = size(data, dim=1) + start(2) = 1; count(2) = size(data, dim=2) + start(3) = 1; count(3) = size(data, dim=3) + start(4) = nt; count(4) = 1 + call check_ret(nf_put_vara_int2(ncid, varid, start, count, data), subname) + else + call check_ret(nf_put_var_int2(ncid, varid, data), subname) + end if + if (present(long_name)) then + call check_ret(nf_put_att_text(ncid, varid, 'long_name', len_trim(long_name), trim(long_name)), subname) + end if + if (present(units)) then + call check_ret(nf_put_att_text(ncid, varid, 'units', len_trim(units), trim(units)), subname) + end if + end if + + else if (flag == 'read') then + + if (nc_masterproc) then + call check_var(ncid, varname, varid, varpresent) + if (varpresent) then + if (present(nt)) then + start(1) = 1; count(1) = size(data, dim=1) + start(2) = 1; count(2) = size(data, dim=2) + start(3) = 1; count(3) = size(data, dim=3) + start(4) = nt; count(4) = 1 + call check_ret(nf_get_vara_int2(ncid, varid, start, count, data), subname) + else + call check_ret(nf_get_var_int2(ncid, varid, data), subname) + end if + else + call endrun('the varibal does not difined!',subname) + endif + end if + if (present(readvar)) readvar = varpresent + + end if + + end subroutine ncd_ioglobal_short_3d +!----------------------------------------------------------------------- +!BOP +! +! !IROUTINE: ncd_ioglobal_int_3d +! +! !INTERFACE: + subroutine ncd_ioglobal_int_3d(varname, data, flag, & + ncid, long_name, units, nt, readvar) +! !DESCRIPTION: +! netcdf I/O of global 3d integer array +! +! !ARGUMENTS: + implicit none + character(len=*), intent(in) :: flag ! 'read' or 'write' + integer , intent(in) :: ncid ! input unit + character(len=*), intent(in) :: varname ! variable name + integer , intent(inout) :: data(:,:,:) ! local decomposition input data + character(len=*), optional, intent(in) :: long_name ! variable long name + character(len=*), optional, intent(in) :: units ! variable units + logical , optional, intent(out):: readvar ! true => variable is on initial dataset (read only) + integer , optional, intent(in) :: nt ! time sample index +! +! !REVISION HISTORY: +! +!EOP +! +! !LOCAL VARIABLES: + integer :: varid ! netCDF variable id + integer :: dimid(4), ndims ! dimension ids + integer :: start(4), count(4) ! output bounds + integer :: ier ! error code + logical :: varpresent ! if true, variable is on tape + character(len=32) :: subname='NCD_IOGLOBAL_3D_INT_IO' ! subroutine name +!----------------------------------------------------------------------- + + if (flag == 'write') then + + if (nc_masterproc) then + call check_ret(nf_inq_varid(ncid, varname, varid), subname) + if (present(nt)) then + start(1) = 1; count(1) = size(data, dim=1) + start(2) = 1; count(2) = size(data, dim=2) + start(3) = 1; count(3) = size(data, dim=3) + start(4) = nt; count(4) = 1 + call check_ret(nf_put_vara_int(ncid, varid, start, count, data), subname) + else + call check_ret(nf_put_var_int(ncid, varid, data), subname) + end if + if (present(long_name)) then + call check_ret(nf_put_att_text(ncid, varid, 'long_name', len_trim(long_name), trim(long_name)), subname) + end if + if (present(units)) then + call check_ret(nf_put_att_text(ncid, varid, 'units', len_trim(units), trim(units)), subname) + end if + end if + + else if (flag == 'read') then + + if (nc_masterproc) then + call check_var(ncid, varname, varid, varpresent) + if (varpresent) then + if (present(nt)) then + start(1) = 1; count(1) = size(data, dim=1) + start(2) = 1; count(2) = size(data, dim=2) + start(3) = 1; count(3) = size(data, dim=3) + start(4) = nt; count(4) = 1 + call check_ret(nf_get_vara_int(ncid, varid, start, count, data), subname) + else + call check_ret(nf_get_var_int(ncid, varid, data), subname) + end if + else + call endrun('the varibal does not difined!',subname) + endif + end if + if (present(readvar)) readvar = varpresent + + end if + + end subroutine ncd_ioglobal_int_3d + +!----------------------------------------------------------------------- +!BOP +! +! !IROUTINE: ncd_ioglobal_real_3d +! +! !INTERFACE: + subroutine ncd_ioglobal_real_3d(varname, data, flag, & + ncid, long_name, units, nt, readvar) +! !DESCRIPTION: +! netcdf I/O of global 3d real array +! +! !ARGUMENTS: + implicit none + character(len=*), intent(in) :: flag ! 'read' or 'write' + integer , intent(in) :: ncid ! input unit + character(len=*), intent(in) :: varname ! variable name + real, intent(inout) :: data(:,:,:) ! local decomposition input data + character(len=*), optional, intent(in) :: long_name ! variable long name + character(len=*), optional, intent(in) :: units ! variable units + logical , optional, intent(out):: readvar ! true => variable is on initial dataset (read only) + integer , optional, intent(in) :: nt ! time sample index +! +! !REVISION HISTORY: +! +!EOP +! +! !LOCAL VARIABLES: + integer :: varid ! netCDF variable id + integer :: ier ! error code + integer :: dimid(4), ndims ! dimension ids + integer :: start(4), count(4) ! output bounds + logical :: varpresent ! if true, variable is on tape + character(len=32) :: subname='NCD_IOGLOBAL_REAL_3D' ! subroutine name +!----------------------------------------------------------------------- + + if (flag == 'write') then + + if (nc_masterproc) then + call check_ret(nf_inq_varid(ncid, varname, varid), subname) + if (present(nt)) then + start(1) = 1; count(1) = size(data, dim=1) + start(2) = 1; count(2) = size(data, dim=2) + start(3) = 1; count(3) = size(data, dim=3) + start(4) = nt; count(4) = 1 + call check_ret(nf_put_vara_real(ncid, varid, start, count, data), subname) + else + call check_ret(nf_put_var_real(ncid, varid, data), subname) + end if + if (present(long_name)) then + call check_ret(nf_put_att_text(ncid, varid, 'long_name', len_trim(long_name), trim(long_name)), subname) + end if + if (present(units)) then + call check_ret(nf_put_att_text(ncid, varid, 'units', len_trim(units), trim(units)), subname) + end if + end if + + else if (flag == 'read') then + + if (nc_masterproc) then + call check_var(ncid, varname, varid, varpresent) + if (varpresent) then + if (present(nt)) then + start(1) = 1; count(1) = size(data, dim=1) + start(2) = 1; count(2) = size(data, dim=2) + start(3) = 1; count(3) = size(data, dim=3) + start(4) = nt; count(4) = 1 + call check_ret(nf_get_vara_real(ncid, varid, start, count, data), subname) + else + call check_ret(nf_get_var_real(ncid, varid, data), subname) + end if + else + call endrun('the varibal does not difined!',subname) + endif + end if + if (present(readvar)) readvar = varpresent + + end if + + end subroutine ncd_ioglobal_real_3d + +!----------------------------------------------------------------------- +!BOP +! +! !IROUTINE: ncd_ioglobal_real_3d +! +! !INTERFACE: + subroutine ncd_ioglobal_double_3d(varname, data, flag, & + ncid, long_name, units, nt, readvar) +! !DESCRIPTION: +! netcdf I/O of global 3d real array +! +! !ARGUMENTS: + implicit none + character(len=*), intent(in) :: flag ! 'read' or 'write' + integer , intent(in) :: ncid ! input unit + character(len=*), intent(in) :: varname ! variable name + real*8, intent(inout) :: data(:,:,:) ! local decomposition input data + character(len=*), optional, intent(in) :: long_name ! variable long name + character(len=*), optional, intent(in) :: units ! variable units + logical , optional, intent(out):: readvar ! true => variable is on initial dataset (read only) + integer , optional, intent(in) :: nt ! time sample index +! +! !REVISION HISTORY: +! +!EOP +! +! !LOCAL VARIABLES: + integer :: varid ! netCDF variable id + integer :: ier ! error code + integer :: dimid(4), ndims ! dimension ids + integer :: start(4), count(4) ! output bounds + logical :: varpresent ! if true, variable is on tape + character(len=32) :: subname='NCD_IOGLOBAL_REAL_3D' ! subroutine name +!----------------------------------------------------------------------- + + if (flag == 'write') then + + if (nc_masterproc) then + call check_ret(nf_inq_varid(ncid, varname, varid), subname) + if (present(nt)) then + start(1) = 1; count(1) = size(data, dim=1) + start(2) = 1; count(2) = size(data, dim=2) + start(3) = 1; count(3) = size(data, dim=3) + start(4) = nt; count(4) = 1 + call check_ret(nf_put_vara_double(ncid, varid, start, count, data), subname) + else + call check_ret(nf_put_var_double(ncid, varid, data), subname) + end if + if (present(long_name)) then + call check_ret(nf_put_att_text(ncid, varid, 'long_name', len_trim(long_name), trim(long_name)), subname) + end if + if (present(units)) then + call check_ret(nf_put_att_text(ncid, varid, 'units', len_trim(units), trim(units)), subname) + end if + end if + + else if (flag == 'read') then + + if (nc_masterproc) then + call check_var(ncid, varname, varid, varpresent) + if (varpresent) then + if (present(nt)) then + start(1) = 1; count(1) = size(data, dim=1) + start(2) = 1; count(2) = size(data, dim=2) + start(3) = 1; count(3) = size(data, dim=3) + start(4) = nt; count(4) = 1 + call check_ret(nf_get_vara_double(ncid, varid, start, count, data), subname) + else + call check_ret(nf_get_var_double(ncid, varid, data), subname) + end if + else + call endrun('the varibal does not difined!',subname) + endif + end if + if (present(readvar)) readvar = varpresent + + end if + + end subroutine ncd_ioglobal_double_3d + +!----------------------------------------------------------------------- +!BOP +! !IROUTINE: endrun +! +! !INTERFACE: +subroutine endrun(msg,subname) +! +! !DESCRIPTION: +! Abort the model for abnormal termination + implicit none +! !ARGUMENTS: + character(len=*), intent(in), optional :: msg ! string to be printed + character(len=*), intent(in), optional :: subname ! 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 ncdio + + + + + + diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/offline/readme.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/offline/readme.txt new file mode 100644 index 000000000..51eba9619 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/offline/readme.txt @@ -0,0 +1,49 @@ +README - River Routing Model Offline Version +Last Updated: 10/28/2024 +Contact: yujin.zeng@nasa.gov +Overview + +This directory contains the code required to run the offline version of the river routing model. Note that not all files in this directory pertain to the offline model. Key files include: + + run: Script for building and running the model. + ncdioMod.f90: Local NetCDF library. + rwncMod.f90: Local NetCDF I/O library. + interp_M36toPfaf.f90: Interpolation module. + river_io_mod.f90: I/O module. + res_mod.f90: Reservoir module. + lake_mod.f90: Lake module. + river_routing.f90: Main program. + +Running the Offline Model + + Set Directory Paths + In river_io_mod.f90, set: + input_dir: Path for input data, e.g., /discover/nobackup/yzeng3/work/river_routing_model_offline/input/ + runoff_dir: Path for runoff data (e.g., Catchment model 2D output in M36 or M09 resolutions). + Example for M36 resolution: /discover/nobackup/yzeng3/GEOldas_output + output_dir: Path for output data. + + Define Start and End Dates + In river_routing.f90, set step_start (start date) and step_end (end date) as days since January 1, 1990 (Day 1). Ensure these dates align with the runoff forcing period. + + Build and Run + Compile and run the model using: + ./run river_routing.f90 + +Output Format + +The output files are in .txt format, generated daily with date information in each filename. The output variables are as follows: + + Main river discharge (Pfaf_Qr_Kv_*.txt) [m³/s] + Main river storage (Pfaf_Wr_Kv_*.txt) [kg] + Local stream storage (Pfaf_Ws_Kv_*.txt) [kg] + Reservoir outflow (Pfaf_Q_res_Kv_*.txt) [m³/s] (0 for catchments without reservoirs) + Reservoir water storage (Pfaf_Wr_res_Kv_*.txt) [kg] (0 for catchments without reservoirs) + Lake outflow (Pfaf_Q_lake_Kv_*.txt) [m³/s] (0 for catchments without lakes) + Lake water storage (Pfaf_Wr_lake_Kv_*.txt) [kg] (0 for catchments without lakes) + +Each .txt file contains a list of 291,284 values corresponding to catchments indexed from 1 to 291,284. To convert these lists into spatial maps, use the catchment distribution map at 1-minute resolution in CatchIndex from SRTM_PfafData.nc: + + Path: /discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/land/topo/v1/SRTM-TopoData/SRTM_PfafData.nc + +For further assistance, please contact yujin.zeng@nasa.gov. \ No newline at end of file diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/offline/res_mod.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/offline/res_mod.f90 new file mode 100644 index 000000000..f3aed2c77 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/offline/res_mod.f90 @@ -0,0 +1,316 @@ +module reservoir + +use rwncfile + +implicit none +private +public :: res_init, res_cal + +!----Reservoir module constants---------- + +real*8, parameter :: fac_elec_a = 0.30D0 ! Coefficient for hydropower calculation +real*8, parameter :: fac_elec_b = 2.00D0 ! Exponent for hydropower calculation +real*8, parameter :: fac_irr_a = 0.01D0 ! Coefficient for irrigation calculation (arid areas) +real*8, parameter :: fac_irr_b = 3.00D0 ! Scaling factor for irrigation (arid areas) +real*8, parameter :: fac_sup_a = 0.03D0 ! Coefficient for water supply calculation +real*8, parameter :: fac_sup_b = 2.00D0 ! Exponent for water supply calculation +real*8, parameter :: fac_other_a = 0.20D0 ! Coefficient for other reservoir types +real*8, parameter :: fac_other_b = 2.00D0 ! Exponent for other reservoir types +integer, parameter :: fac_fld = 1 ! Flood control parameter + +real*8, parameter :: dt = 86400.D0 ! Time step in seconds (1 day) + +real*8, parameter :: ai_thres = 0.5D0 ! Aridity index threshold for irrigation reservoirs +real*8, parameter :: rho = 1.D3 ! Water density (kg/m^3) + +!----------------------------------------- + +contains + +!------------------------------------------ +! Initialization subroutine for reservoirs +subroutine res_init(input_dir,nres,nc,use_res,active_res,Wr_res,Q_res,type_res,cap_res,Qavg_res,ai_res,fld_res,Qfld_thres,irr_sea_frac,cat2res,wid_res) + character(len=500),intent(in) :: input_dir + ! Define the number of reservoirs (nres) and the number of catchments (nc) + integer,intent(in) :: nres,nc + ! Logical variable to check if reservoirs are used + logical,intent(in) :: use_res + ! Input/output arrays for reservoir attributes: active reservoirs, types, capacities, etc. + integer,intent(inout),allocatable :: active_res(:),type_res(:),fld_res(:),cat2res(:) + real*8,intent(inout),allocatable :: Wr_res(:),Q_res(:),cap_res(:),Qavg_res(:),ai_res(:),Qfld_thres(:),irr_sea_frac(:,:) + real*8,intent(inout),allocatable :: wid_res(:) + + ! Internal arrays for various reservoir-related data + integer,allocatable,dimension(:) :: flag_grand,catid_grand,elec_grand,irrsup_grand,fld_grand,supply_grand,irr_grand,realuse_grand + integer,allocatable,dimension(:) :: nav_grand,rec_grand,other_grand + real*8,allocatable,dimension(:) :: cap_grand,area_max_res,Qavg_grand,ai_grand,area_grand,power_grand,area_res + real*8,allocatable,dimension(:,:) :: Wres_tar + + ! Define the flood threshold variable and a counter variable + character(len=2) :: fld_thres + integer :: i,cid,rid + +!----------reservoir module-------------- + ! Allocate memory for each array + allocate(flag_grand(nres),catid_grand(nres),active_res(nc)) + allocate(Wr_res(nc),Q_res(nc)) + allocate(elec_grand(nres),type_res(nc),cap_grand(nres),cap_res(nc),area_grand(nres)) + allocate(area_res(nc),area_max_res(nc)) + !allocate(irrsup_grand(nres)) + allocate(fld_grand(nres),fld_res(nc),Qfld_thres(nc),supply_grand(nres)) + allocate(irr_grand(nres)) + allocate(cat2res(nc)) + allocate(nav_grand(nres),rec_grand(nres)) + allocate(other_grand(nres)) + allocate(wid_res(nc)) + allocate(realuse_grand(nres)) + + ! Open reservoir-related data files and read the corresponding arrays + open(77,file=trim(input_dir)//"/catid_dam_corr_aca_grand5000.txt") + read(77,*)catid_grand + open(77,file=trim(input_dir)//"/flag_all_res.txt") + read(77,*)flag_grand + open(77,file=trim(input_dir)//"/cap_max_grand.txt") + read(77,*)cap_grand + cap_grand=cap_grand*1.D6*rho ! Convert capacity from million cubic meters (MCM) to kilograms (kg) + open(77,file=trim(input_dir)//"/hydroelec_grand.txt") + read(77,*)elec_grand + !open(77,file=trim(input_dir)//"/Qavg_res_2016_2020_OL7000.txt") + !read(77,*)Qavg_grand + !Qavg_grand=Qavg_grand*rho ! Convert flow rate from cubic meters per second (m3/s) to kilograms per second (kg/s) + !open(77,file=trim(input_dir)//"/ai_grand.txt") + !read(77,*)ai_grand + !open(77,file=trim(input_dir)//"/irrmainsec_noelec_grand.txt") + !read(77,*)irrsup_grand + open(77,file=trim(input_dir)//"/fldmainsec_grand.txt") + read(77,*)fld_grand + write(fld_thres,'(I2.2)')fac_fld + open(77,file=trim(input_dir)//"/Pfaf_flood_qr_thres"//trim(fld_thres)//".txt") + read(77,*)Qfld_thres ! Read flood thresholds in cubic meters per second (m3/s) + Qfld_thres=Qfld_thres*rho ! Convert threshold from cubic meters per second to kilograms per second (kg/s) + open(77,file=trim(input_dir)//"/watersupply_grand.txt") + read(77,*)supply_grand + open(77,file=trim(input_dir)//"/irr_grand.txt") + read(77,*)irr_grand + open(77,file=trim(input_dir)//"/nav_grand.txt") + read(77,*)nav_grand + open(77,file=trim(input_dir)//"/rec_grand.txt") + read(77,*)rec_grand + open(77,file=trim(input_dir)//"/other_grand.txt") + read(77,*)other_grand + open(77,file=trim(input_dir)//"/area_skm_grand.txt") + read(77,*)area_grand + area_grand=area_grand*1.D6 ! Convert area from square kilometers (km2) to square meters (m2) + !open(77,file=trim(input_dir)//"/power_grand.txt") + !read(77,*)power_grand + + ! Set initial reservoir ID mapping + cat2res=0 + do i=1,nres + if(flag_grand(i)==1)then + cid=catid_grand(i) + cat2res(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.D0 ! Set reservoir capacity to zero + area_res = 0.D0 ! Set reservoir area to zero + area_max_res = 0.D0 ! Set max reservoir area to zero + type_res = 0 ! Set reservoir type to zero + !Qavg_res = 0.D0 ! Set average reservoir flow rate to zero + !ai_res = 0.D0 ! Set irrigation index 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 + 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 + cap_res(cid) = cap_res(cid) + cap_grand(i) ! Sum up the capacities for reservoirs in the same catchment + area_res(cid) = area_res(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_res(cid) = 1 ! Mark the catchment if it has flood control + endif + enddo + + ! Compute reservoir width from area (square root of the area) + wid_res = sqrt(area_res) + + ! 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(cid) = 6 ! Type 7 for other uses + cat2res(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(cid) = 5 + cat2res(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(cid) = 4 + cat2res(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(cid) = 3 + cat2res(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(cid) = 2 + cat2res(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(cid) = 1 + cat2res(cid) = i + area_max_res(cid) = area_grand(i) + endif + endif + enddo + + ! 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 + + ! Assign real reservoir usage based on type, with error checking + do i = 1, nres + if(flag_grand(i) == 1) then + cid = catid_grand(i) + rid = cat2res(cid) + if(rid > 0) then + if(type_res(cid) == 0 .and. fld_res(cid) == 0) then + print *, "type_res(cid) == 0" + stop + endif + if(type_res(cid) == 0) then + realuse_grand(i) = -1 ! Invalid reservoir use type + else + realuse_grand(i) = type_res(cid) ! Assign the actual use type + endif + else + print *, "rid == 0" + stop + endif + endif + enddo + + ! Read irrigation and reservoir target data from NetCDF files + ! call read_ncfile_double2d(trim(input_dir)//"/irr_grand_frac.nc", "data", irr_sea_frac, nres, 12) + ! call read_ncfile_double2d(trim(input_dir)//"/Wr_tar_Dang.nc", "data", Wres_tar, 365, nres) + + ! Wres_tar = Wres_tar * 1.D6 * rho ! Convert from million cubic meters (MCM) to kilograms (kg) + + ! Deactivate reservoirs if the use_res flag is set to False + if(use_res == .False.) active_res = 0 + +end subroutine res_init + +!----------------------- +! Reservoir calculation subroutine +subroutine res_cal(active_res,active_lake,Qout,Q_lake,type_res,cat2res,Q_res,wid_res,fld_res,Wr_res,Qfld_thres,cap_res,B1,B2) + integer, intent(in) :: active_res, type_res, active_lake, cat2res, fld_res + real*8, intent(in) :: Qout, Q_lake, wid_res, Qfld_thres, cap_res + real*8, intent(inout) :: Q_res, Wr_res, B1, B2 + + integer :: rid ! Reservoir ID + real*8 :: Qin_res, coe, irrfac, alp_res ! Variables for inflow, coefficients, and factors + + ! If the reservoir is active + if (active_res == 1) then + + ! Determine the inflow to the reservoir (from river or lake) + if (active_lake == 0) then + Qin_res = Qout ! Inflow from river + else + Qin_res = Q_lake ! Inflow from lake + endif + + ! Irrigation reservoir + if (type_res == 1) then + alp_res = fac_irr_a * ((1.D0 / (wid_res / 1.D3)) ** fac_irr_b) / 3600.D0 ! irrigation coefficient + Q_res = alp_res * Wr_res ! Outflow based on water storage + + ! Hydropower reservoir + else if (type_res == 2) then + alp_res = fac_elec_a * ((1.D0 / (wid_res / 1.D3)) ** fac_elec_b) / 3600.D0 ! Hydropower coefficient + Q_res = alp_res * Wr_res ! Outflow based on water storage + + ! Water supply reservoir + else if (type_res == 3) then + alp_res = fac_sup_a * ((1.D0 / (wid_res / 1.D3)) ** fac_sup_b) / 3600.D0 ! Supply coefficient + Q_res = alp_res * Wr_res ! Outflow based on water storage + + ! Other reservoir types + else if (type_res == 4 .or. type_res == 5 .or. type_res == 6 .or. type_res == 0) then + alp_res = fac_other_a * ((1.D0 / (wid_res / 1.D3)) ** fac_other_b) / 3600.D0 ! Generic reservoir coefficient + Q_res = alp_res * Wr_res ! Outflow based on water storage + endif + + ! Ensure outflow is within reasonable bounds + Q_res = max(0.D0, Q_res) ! Ensure non-negative outflow + Q_res = min(Q_res, Wr_res / 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 = Wr_res + dt * (Qin_res - Q_res) ! Update water storage in the reservoir + Wr_res = max(0.D0, Wr_res) ! Ensure non-negative storage + + ! If the storage exceeds capacity, adjust outflow and storage + if (Wr_res > cap_res) then + Q_res = Q_res + (Wr_res - cap_res) / dt ! Adjust outflow for overflow + Wr_res = cap_res ! Limit storage to reservoir capacity + endif + + ! Output the calculated outflow and zero out the second output variable (B2) + B1 = Q_res + B2 = 0.D0 + + endif + +end subroutine res_cal + +end module reservoir diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/offline/river_io_mod.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/offline/river_io_mod.f90 new file mode 100644 index 000000000..d8472d8fe --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/offline/river_io_mod.f90 @@ -0,0 +1,319 @@ +module river_io + +use interp +use rwncfile + +implicit none +private + +public :: read_input,read_restart,read_runoff,write_output + +real*8, parameter :: rho = 1.D3 ! Water density in kg/m^3 +character(len=500) :: input_dir="/discover/nobackup/yzeng3/work/river_routing_model_offline/input/" ! Directory for input files +character(len=500) :: output_dir="/discover/nobackup/yzeng3/river_output/" ! Directory for output files +character(len=500) :: runoff_dir="/discover/nobackup/yzeng3/GEOldas_output/" ! Directory for runoff files + +integer :: nlon=964 !for M36, change to 3856 for M09 +integer :: nlat=406 !for M36, change to 1624 for M09 + +contains + +!------------------------------ +subroutine read_input(nc,ny,upmax,days_in_year,fac_kstr,qstr_clmt,qri_clmt,nts,upID,nup,llc_ori,lstr,qin_clmt,K,Kstr,days_acc_year,days_acc_noleap,days_acc_leap,inputdir) + ! Input parameters: + integer,intent(in) :: nc, ny, upmax ! nc: number of catchments, ny: number of years, upmax: max number of upstream catchments + integer,intent(in) :: days_in_year(ny) ! Array of days in each year + real*8,intent(in) :: fac_kstr ! Scaling factor for streamflow + real*8,intent(out) :: qstr_clmt(nc), qri_clmt(nc) ! Climate streamflow (qstr_clmt) and routing inflow (qri_clmt) in kg/s + integer,intent(out) :: nts(nc), upID(upmax,nc), nup(nc) ! Number of time steps, upstream IDs, and number of upstream catchments + real*8,intent(out) :: llc_ori(nc), lstr(nc), qin_clmt(nc), K(nc), Kstr(nc) ! Original stream length (llc_ori), stream length (lstr), climate inflow (qin_clmt), and hydraulic parameters (K, Kstr) + integer,intent(out) :: days_acc_year(ny), days_acc_noleap(12), days_acc_leap(12) ! Accumulated days in regular and leap years + character(len=500),intent(out) :: inputdir + + ! Days in each month for no-leap and leap years + integer,dimension(12) :: days_in_mon_noleap=(/31,28,31,30,31,30,31,31,30,31,30,31/) + integer,dimension(12) :: days_in_mon_leap=(/31,29,31,30,31,30,31,31,30,31,30,31/) + integer :: i + + inputdir=input_dir + ! Read input data from files + open(77,file=trim(input_dir)//"/Pfaf_qstr.txt") + read(77,*)qstr_clmt ! Read streamflow climatology (m3/s) + qstr_clmt=qstr_clmt*rho ! Convert to kg/s + + open(77,file=trim(input_dir)//"/Pfaf_qri.txt") + read(77,*)qri_clmt ! Read routing inflow (m3/s) + qri_clmt=qri_clmt*rho ! Convert to kg/s + + open(77,file=trim(input_dir)//"/Pfaf_qin.txt") + read(77,*)qin_clmt ! Read climate inflow (m3/s) + qin_clmt=qin_clmt*rho ! Convert to kg/s + + open(77,file=trim(input_dir)//"/Pfaf_tosink.txt") + read(77,*)nts ! Read number of steps to endpoint + + open(77,file=trim(input_dir)//"/upstream_1D.txt") + read(77,*)upID ! Read upstream IDs + + open(77,file=trim(input_dir)//"/Pfaf_upnum.txt") + read(77,*)nup ! Read number of upstream catchments + + open(77,file=trim(input_dir)//"/Pfaf_lriv_PR.txt") + read(77,*)llc_ori ! Read original stream length (km) + llc_ori=llc_ori*1.D3 ! Convert km to meters + + open(77,file=trim(input_dir)//"/Pfaf_lstr_PR.txt") + read(77,*)lstr ! Read stream length (km) + lstr=lstr*1.D3 ! Convert km to meters + + open(77,file=trim(input_dir)//"Pfaf_Kv_PR_0p35_0p45_0p2_n0p2.txt") + read(77,*)K ! Read hydraulic parameter K + + open(77,file=trim(input_dir)//"Pfaf_Kstr_PR_fac1_0p35_0p45_0p2_n0p2.txt") + read(77,*)Kstr ! Read hydraulic parameter Kstr + Kstr=fac_kstr*Kstr ! Apply scaling factor to Kstr + + ! Calculate accumulated days for regular years + days_acc_year(1)=0 + do i=2,ny + days_acc_year(i)=days_acc_year(i-1)+days_in_year(i-1) + end do + + ! Calculate accumulated days for no-leap and leap years + days_acc_noleap(1)=0 + days_acc_leap(1)=0 + do i=2,12 + days_acc_noleap(i)=days_acc_noleap(i-1)+days_in_mon_noleap(i-1) + days_acc_leap(i)=days_acc_leap(i-1)+days_in_mon_leap(i-1) + end do + +end subroutine read_input +!------------------------------ +subroutine read_restart(iter,is_coldstart,ny,nc,days_acc_year,days_acc_noleap,days_acc_leap,Ws,Wr,Wr_res,Wr_lake) + ! Input parameters: + integer,intent(in) :: iter ! Current iteration + logical,intent(inout) :: is_coldstart ! Flag for cold start condition + integer,intent(in) :: ny, nc ! ny: number of years, nc: number of catchments + integer,intent(in) :: days_acc_year(ny), days_acc_noleap(12), days_acc_leap(12) ! Accumulated days for each year and for no-leap/leap years + real*8,intent(inout) :: Ws(nc), Wr(nc), Wr_res(nc), Wr_lake(nc) ! Water storage in soil (Ws), routing (Wr), reservoir (Wr_res), and lake (Wr_lake) + + ! Local variables: + character(len=50) :: iter_s, yr_s, mon_s, day_s ! Strings for iteration, year, month, and day + integer :: step_prev, i, yr_cur, mon_cur, day_cur, d_res ! Step count, loop index, current year, month, day, and day residual + integer :: days_acc_mon(12) ! Accumulated days per month + + ! Convert iteration number to string format + write(iter_s,'(I5.5)')iter + print *,trim(iter_s) + + ! If first iteration or cold start, read initial data + if(iter==1.or.is_coldstart)then + ! Read initial water storage data from files for cold start + open(77,file=trim(input_dir)//"/Pfaf_Ws_Kv_M0.10_mm0.40_20170330_OL7000.txt") + read(77,*)Ws ! Read soil water storage (Ws) + + open(77,file=trim(input_dir)//"/Pfaf_Wr_Kv_M0.10_mm0.40_20170330_OL7000.txt") + read(77,*)Wr ! Read routing water storage (Wr) + + !----reservoir module------- + open(77,file=trim(input_dir)//"/Pfaf_Wr_res_Kv_M0.10_mm0.40_20170330_OL7000.txt") + read(77,*)Wr_res ! Read reservoir water storage (Wr_res) + + !----lake module------------ + open(77,file=trim(input_dir)//"/Pfaf_Wr_lake_Kv_M0.10_mm0.40_20170330_OL7000.txt") + read(77,*)Wr_lake ! Read lake water storage (Wr_lake) + + ! Set cold start flag to False after initialization + is_coldstart=.False. + + else + ! For non-cold start, calculate the current year and day from the previous iteration + step_prev = iter - 1 + do i = ny, 1, -1 + if(step_prev > days_acc_year(i))then + yr_cur = 1989 + i ! Calculate the current year + d_res = step_prev - days_acc_year(i) ! Calculate residual days + exit + endif + enddo + + ! Determine whether the current year is a leap year + if(mod(yr_cur,4) == 0)then + days_acc_mon = days_acc_leap ! Use leap year days if it is a leap year + else + days_acc_mon = days_acc_noleap ! Use no-leap year days if it is not a leap year + endif + + ! Determine the current month and day from the residual days + do i = 12, 1, -1 + if(d_res > days_acc_mon(i))then + mon_cur = i ! Current month + day_cur = d_res - days_acc_mon(i) ! Current day + exit + endif + enddo + + ! Convert year, month, and day to string format + write(yr_s,'(I4)')yr_cur + write(mon_s,'(I2.2)')mon_cur + write(day_s,'(I2.2)')day_cur + + ! Read water storage data for the specific date (year, month, day) + open(77,file=trim(output_dir)//"/Pfaf_Ws_Kv_"//trim(yr_s)//trim(mon_s)//trim(day_s)//"_OL7000.txt") + read(77,*)Ws ! Read soil water storage (Ws) + + open(77,file=trim(output_dir)//"/Pfaf_Wr_Kv_"//trim(yr_s)//trim(mon_s)//trim(day_s)//"_OL7000.txt") + read(77,*)Wr ! Read routing water storage (Wr) + + !----reservoir module------- + open(77,file=trim(output_dir)//"Pfaf_Wr_res_Kv_"//trim(yr_s)//trim(mon_s)//trim(day_s)//"_OL7000.txt") + read(77,*)Wr_res ! Read reservoir water storage (Wr_res) + + !----lake module------------ + open(77,file=trim(output_dir)//"Pfaf_Wr_lake_Kv_"//trim(yr_s)//trim(mon_s)//trim(day_s)//"_OL7000.txt") + read(77,*)Wr_lake ! Read lake water storage (Wr_lake) + + ! Optionally scale the water storage values (commented out) + ! Ws = Ws * 1.D9 + ! Wr = Wr * 1.D9 + endif + +end subroutine read_restart +!------------------------------ +subroutine read_runoff(nc,ny,iter,days_acc_year,days_acc_noleap,days_acc_leap,Qrunf,yr_s,mon_s,day_s,d_res,mon_cur) + integer,intent(in) :: nc,ny,iter + integer,intent(in) :: days_acc_year(ny),days_acc_noleap(12),days_acc_leap(12) + real*8,intent(inout) :: Qrunf(nc) + character(len=50),intent(inout) :: yr_s,mon_s,day_s + integer,intent(out) :: d_res,mon_cur + + real*8,allocatable,dimension(:,:,:) :: runoff,runoffr,baseflow ! Declare 3D arrays for runoff and baseflow + + integer :: i,yr_cur,day_cur + integer :: days_acc_mon(12) ! Array to store accumulated days for current month + + + ! Determine current year based on iteration days + do i=ny,1,-1 + if(iter>days_acc_year(i))then + yr_cur=1989+i ! Set current year + d_res=iter-days_acc_year(i) ! Calculate residual days + exit + endif + enddo + + ! Set days_acc_mon based on whether the current year is a leap year + if(mod(yr_cur,4)==0)then + days_acc_mon=days_acc_leap ! Use leap year days + else + days_acc_mon=days_acc_noleap ! Use non-leap year days + endif + + ! Determine current month and day based on residual days + do i=12,1,-1 + if(d_res>days_acc_mon(i))then + mon_cur=i ! Set current month + day_cur=d_res-days_acc_mon(i) ! Set current day + exit + endif + enddo + + ! Write current year, month, and day as strings + write(yr_s,'(I4)')yr_cur + write(mon_s,'(I2.2)')mon_cur + write(day_s,'(I2.2)')day_cur + print *,trim(yr_s)," ",trim(mon_s)," ",trim(day_s) + + ! Allocate memory for runoff, runoffr, and baseflow arrays + allocate(runoff(nlon,nlat,1),runoffr(nlon,nlat,1),baseflow(nlon,nlat,1)) + + ! Read runoff and baseflow data from NetCDF files + call read_ncfile_double3d(trim(runoff_dir)//"/Y"//trim(yr_s)//"/M"//trim(mon_s)//"/SMAP_Nature_v10.0_M36.tavg24_2d_lnd_Nx."//trim(yr_s)//trim(mon_s)//trim(day_s)//"_1200z.nc4","RUNOFF",runoff,nlon,nlat,1) + call read_ncfile_double3d(trim(runoff_dir)//"/Y"//trim(yr_s)//"/M"//trim(mon_s)//"/SMAP_Nature_v10.0_M36.tavg24_2d_lnd_Nx."//trim(yr_s)//trim(mon_s)//trim(day_s)//"_1200z.nc4","BASEFLOW",baseflow,nlon,nlat,1) + + ! Combine runoff and baseflow, and convert to daily values + runoff=runoff+baseflow + runoff=runoff*86400.D0 ! Convert to mm/day + + ! Reverse the y-direction of the runoff array + do i=1,406 + runoffr(:,i,:)=runoff(:,407-i,:) + enddo + runoff=runoffr + + ! Convert from mm/day to kg/s and store in Qrunf + Qrunf=M36_to_cat(runoff(:,:,1),nlon,nlat,nc,input_dir) + + ! Deallocate the arrays to free memory + deallocate(runoff,runoffr,baseflow) + + ! The following lines are commented out, but they suggest reading runoff from a text file instead of NetCDF + !open(77,file="/Users/zsp/Desktop/work/river/OL7000_Pfaf/runoff_"//trim(yr_s)//trim(mon_s)//trim(day_s)//".txt") + !read(77,*)Qrunf + !Qrunf=Qrunf*rho !m3/s -> kg/s + +end subroutine read_runoff +!------------------------------ +subroutine write_output(nc,yr_s,mon_s,day_s,Qout,Ws,Wr,Q_res,Wr_res,Q_lake,Wr_lake) + integer,intent(in) :: nc + character(len=50),intent(in) :: yr_s,mon_s,day_s + real*8,intent(in) :: Qout(nc),Ws(nc),Wr(nc),Q_res(nc),Wr_res(nc),Q_lake(nc),Wr_lake(nc) + + integer :: i + + ! Open file to write Qout (discharge) values and write to the file + open(88,file=trim(output_dir)//"/Pfaf_Qr_Kv_"//trim(yr_s)//trim(mon_s)//trim(day_s)//"_OL7000.txt") + do i=1,nc + write(88,*)Qout(i)/1.D3 ! Convert from m^3/s to km^3/s + enddo + + ! Open file to write Ws (soil water storage) values and write to the file + open(88,file=trim(output_dir)//"/Pfaf_Ws_Kv_"//trim(yr_s)//trim(mon_s)//trim(day_s)//"_OL7000.txt") + do i=1,nc + write(88,*)Ws(i) ! Write Ws values, unit in kg + enddo + + ! Open file to write Wr (river water storage) values and write to the file + open(88,file=trim(output_dir)//"/Pfaf_Wr_Kv_"//trim(yr_s)//trim(mon_s)//trim(day_s)//"_OL7000.txt") + do i=1,nc + write(88,*)Wr(i) ! Write Wr values, unit in kg + enddo + + !-----------reservoir module---------------- + ! Open file to write Q_res (reservoir discharge) values and write to the file + open(88,file=trim(output_dir)//"/Pfaf_Q_res_Kv_"//trim(yr_s)//trim(mon_s)//trim(day_s)//"_OL7000.txt") + do i=1,nc + write(88,*)Q_res(i)/1.D3 ! Convert from m^3/s to km^3/s + enddo + + ! Open file to write Wr_res (reservoir water storage) values and write to the file + open(88,file=trim(output_dir)//"/Pfaf_Wr_res_Kv_"//trim(yr_s)//trim(mon_s)//trim(day_s)//"_OL7000.txt") + do i=1,nc + write(88,*)Wr_res(i) ! Write Wr_res values, unit in kg + enddo + !------------------------------------------- + + !-----------lake module--------------------- + ! Open file to write Q_lake (lake discharge) values and write to the file + open(88,file=trim(output_dir)//"/Pfaf_Q_lake_Kv_"//trim(yr_s)//trim(mon_s)//trim(day_s)//"_OL7000.txt") + do i=1,nc + write(88,*)Q_lake(i)/1.D3 ! Convert from m^3/s to km^3/s + enddo + + ! Open file to write Wr_lake (lake water storage) values and write to the file + open(88,file=trim(output_dir)//"/Pfaf_Wr_lake_Kv_"//trim(yr_s)//trim(mon_s)//trim(day_s)//"_OL7000.txt") + do i=1,nc + write(88,*)Wr_lake(i) ! Write Wr_lake values, unit in kg + enddo + !------------------------------------------- + + ! Print out the sum of Wr (river water storage) in petagrams (10^12 kg) + print *,"sum of Wr is ", sum(Wr)/1.D12 + ! Print out the sum of Wr_lake (lake water storage) in petagrams (10^12 kg) + print *,"sum of Wr_lake is ", sum(Wr_lake)/1.D12 + ! Print out the sum of Wr_res (reservoir water storage) in petagrams (10^12 kg) + print *,"sum of Wr_res is ", sum(Wr_res)/1.D12 + +end subroutine write_output + +end module river_io diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/offline/river_routing.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/offline/river_routing.f90 new file mode 100644 index 000000000..f17167930 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/offline/river_routing.f90 @@ -0,0 +1,248 @@ +program main + +use omp_lib ! OpenMP library for parallel computing +use reservoir ! Module for reservoir operations +use lake ! Module for lake operations +use river_io ! Module for river input/output + +implicit none + +! Define parameters and constants +real*8, parameter :: small = 1.D-48 ! A small value threshold for numerical comparisons +integer, parameter :: step_start = 9221 ! Start timestep (represents 1990-01-01) +integer, parameter :: step_end = 9226 ! End timestep (adjusted for different ranges) +logical :: is_coldstart = .True. ! Logical flag for cold start +integer, parameter :: ny = 33 ! Number of years (33 years) + +real*8, parameter :: fac_kstr = 0.01D0 ! Factor for local stream scaling +real*8, parameter :: M = 0.45D0 ! Parameter in hydraulic geometry formula +real*8, parameter :: mm = 0.35D0 ! Parameter in hydraulic geometry formula + +real*8, parameter :: dt = 86400.D0 ! Time step in seconds (1 day) + +integer, parameter :: nmax = 373 ! Maximum number of catchments in a river +integer, parameter :: upmax = 34 ! Maximum number of upstream basins +integer, parameter :: nc = 291284 ! Total number of river cells +real*8, parameter :: rho = 1.D3 ! Water density in kg/m^3 + +! Declare variables +integer :: i, j, n, iter ! Loop indices and iteration variable + +! Allocate dynamic arrays for variables +integer, allocatable, dimension(:) :: nts ! Array for timestep indices +real*8, allocatable, dimension(:) :: qstr_clmt, qri_clmt, qin_clmt, & + llc_ori, llc, lstr, & + Qrunf, nume, deno, & + alp_s, alp_r, K, Kstr +real*8, allocatable, dimension(:) :: Ws, Wr ! Water storage arrays for stream and river +real*8, allocatable, dimension(:) :: Qs0, ks, Ws_last, Qs, & + Qr0, kr, Cl, Al +real*8, allocatable, dimension(:) :: C1, C2, Qout, Qin, A1, P, B1, B2 +integer, allocatable, dimension(:) :: nup ! Number of upstream nodes +integer, allocatable, dimension(:,:) :: upID ! IDs of upstream cells +real*8 :: co1, co2, co3 ! Coefficients used in calculations +integer :: ui ! Temporary upstream index variable + +real*8, allocatable, dimension(:) :: lon, lat ! Longitude and latitude arrays + +! Reservoir module variables +logical, parameter :: use_res = .True. ! Flag to enable reservoir module +integer, parameter :: nres = 7250 ! Number of reservoirs +integer, allocatable, dimension(:) :: active_res, fld_res, cat2res ! Reservoir attributes +real*8, allocatable, dimension(:) :: Wr_res, Q_res, cap_res, Qavg_res, ai_res, Qfld_thres, wid_res +integer, allocatable, dimension(:) :: type_res ! Type of reservoir (0=inactive, 1-7=different functions) +real*8, allocatable, dimension(:,:) :: irr_sea_frac ! Irrigation and sea fraction for reservoirs + +! Lake module variables +logical, parameter :: use_lake = .True. ! Flag to enable lake module +integer, parameter :: nlake = 3917 ! Number of lakes +integer, allocatable, dimension(:) :: active_lake ! Active lake flag +real*8, allocatable, dimension(:) :: area_lake, Wr_lake, Q_lake ! Lake attributes + +! Time-related variables +integer,dimension(ny) :: days_in_year=(/365,365,366,365,& + 365,365,366,365,& + 365,365,366,365,& + 365,365,366,365,& + 365,365,366,365,& + 365,365,366,365,& + 365,365,366,365,& + 365,365,366,365,365/) ! Number of days per year from 1990 to 2020 +integer :: days_acc_year(ny), days_acc_noleap(12), days_acc_leap(12) ! Accumulated days for leap and non-leap years +integer :: yr_cur, mon_cur, day_cur, d_res, step_prev ! Current date variables and previous step +character(len=50) :: yr_s, mon_s, day_s ! Year, month, day strings +character(len=500) :: inputdir ! Input directory path + +! Allocate memory for variables +allocate(nts(nc)) +allocate(qstr_clmt(nc), qri_clmt(nc), qin_clmt(nc)) +allocate(llc_ori(nc), llc(nc), lstr(nc)) +allocate(Qrunf(nc), nume(nc), deno(nc), alp_s(nc), alp_r(nc)) +allocate(Ws(nc), Wr(nc)) +allocate(Qs0(nc), ks(nc), Ws_last(nc), Qs(nc)) +allocate(Qr0(nc), kr(nc), Cl(nc), Al(nc)) +allocate(C1(nc), C2(nc), Qout(nc), Qin(nc), A1(nc), P(nc), B1(nc), B2(nc)) +allocate(nup(nc)) +allocate(upID(upmax,nc)) +allocate(K(nc), Kstr(nc)) + +! Read input data +call read_input(nc, ny, upmax, days_in_year, fac_kstr, qstr_clmt, qri_clmt, nts, upID, nup, llc_ori, lstr, qin_clmt, K, Kstr, days_acc_year, days_acc_noleap, days_acc_leap, inputdir) + +! Initialize reservoir module +call res_init(inputdir, nres, nc, use_res, active_res, Wr_res, Q_res, type_res, cap_res, Qavg_res, ai_res, fld_res, Qfld_thres, irr_sea_frac, cat2res, wid_res) + +! Initialize lake module +call lake_init(inputdir, use_lake, nc, nlake, nres, active_res, active_lake, area_lake, Wr_lake, Q_lake) + +! Calculate llc (length of river channel) +nume = qri_clmt**(2.D0-M) - qin_clmt**(2.D0-M) ! Numerator for the llc calculation +deno = (2.D0-M) * (qri_clmt - qin_clmt) * (qri_clmt**(1.D0-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.5D0 ! Set llc to half of original value if denominator is small + +! Calculate alp_s (slope coefficient) and alp_r (river coefficient) +where(qstr_clmt > small) alp_s = (rho**(-M) * qstr_clmt**(M-mm) * Kstr * (0.5D0*lstr)**(-1.D0))**(1.D0/(1.D0-mm)) ! For non-zero streamflow +where(qstr_clmt <= small) alp_s = 0.D0 ! 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.D0))**(1.D0/(1.D0-mm)) ! For non-zero river input +where(qri_clmt <= small) alp_r = 0.D0 ! If river input is too small, set alp_r to 0 + +!temporal loop +DO iter=step_start,step_end + + ! Read the state of the system from a restart file for the current iteration + call read_restart(iter,is_coldstart,ny,nc,days_acc_year,days_acc_noleap,days_acc_leap,Ws,Wr,Wr_res,Wr_lake) + + ! Read runoff data for the current time step + call read_runoff(nc,ny,iter,days_acc_year,days_acc_noleap,days_acc_leap,Qrunf,yr_s,mon_s,day_s,d_res,mon_cur) + + !$omp parallel default(shared) + !$omp workshare + + ! Update state variables: ks, Ws, and Qs + where(Qrunf<=small)Qrunf=0.D0 ! Set runoff to zero if it's too small + Qs0=max(0.D0,alp_s * Ws**(1.D0/(1.D0-mm))) ! Initial flow from stream storage (kg/s) + ks=max(0.D0,(alp_s/(1.D0-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.D0-exp(-ks*dt)) ! Update storage (kg) + where(ks<=small) Ws=Ws + (Qrunf-Qs0)*dt ! Simplified update if ks is small + Ws=max(0.D0,Ws) ! Ensure storage is non-negative + Qs=max(0.D0,Qrunf-(Ws-Ws_last)/dt) ! Calculate the stream flow (kg/s) + + ! Calculate variables related to river routing: Qr0, kr + Qr0=max(0.D0,alp_r * Wr**(1.D0/(1.D0-mm))) ! River flow based on water storage (kg/s) + kr=max(0.D0,(alp_r/(1.D0-mm)) * Wr**(mm/(1.D0-mm))) ! Flow coefficient for river (s^-1) + + ! Update Cl and Al + where(kr>small.and.abs(kr-ks)>small) Cl=Wr + (Qrunf-Qr0)/kr*(1.D0-exp(-kr*dt)) + (Qrunf-Qs0)/(kr-ks)*(exp(-kr*dt)-exp(-ks*dt)) + where(kr>small.and.abs(kr-ks)<=small) Cl=Wr + (Qrunf-Qr0)/kr*(1.D0-exp(-kr*dt)) - (Qrunf-Qs0)*dt*exp(-kr*dt) + where(kr<=small.and.ks>small) Cl=Wr + (Qrunf-Qr0)*dt - (Qrunf-Qs0)/ks*(1.D0-exp(-ks*dt)) + where(kr<=small.and.ks<=small) Cl=Wr + (Qs0-Qr0)*dt + Al=Qs+Wr/dt-Cl/dt ! Update flow variables + + ! Initialize variables for river routing process + C1=0.D0 + C2=0.D0 + Qin=0.D0 + Qout=0.D0 + A1=0.D0 + P=0.D0 + B1=0.D0 + B2=0.D0 + + !$omp end workshare + !$omp end parallel + + ! Reservoir module: reset reservoir flow + Q_res=0.D0 + if(d_res==366)d_res=365 ! Handle leap year day adjustment + + ! Lake module: reset lake flow + Q_lake=0.D0 + + ! Process river routing by going through each node from upstream to downstream + do n=nmax,0,-1 + + !$OMP PARALLEL default(shared) private(i,j,ui,co1,co2,co3) + !$OMP DO + + ! Loop over each catchment to update the water storage and flow + do i=1,nc + if(nts(i)==n)then ! If the current node matches the iteration step + + ! Process upstream dependencies if any exist + if(nup(i)>=1)then + do j=1,nup(i) + ui=upID(j,i) + if(ui==-1)exit ! Exit loop if no more upstream IDs + + ! Calculate flow coefficients based on flow conditions + if(kr(i)>small)then + co1=max(0.D0,(1.D0-exp(-kr(i)*dt))/kr(i)) + else + co1=dt + endif + C1(i)=C1(i)+co1*B1(ui) + + if(abs(kr(i)-kr(ui))>small)then + co2=-(exp(-kr(i)*dt)-exp(-kr(ui)*dt))/(kr(i)-kr(ui)) + else + co2=dt*exp(-kr(i)*dt) + endif + C2(i)=C2(i)+co2*B2(ui) + + ! Process reservoir and lake flows, if active + if(active_res(ui)==1.and.active_lake(ui)==0)then + Qin(i)=Qin(i)+Q_res(ui) + else if(active_res(ui)==0.and.active_lake(ui)==1)then + Qin(i)=Qin(i)+Q_lake(ui) + else if(active_res(ui)==1.and.active_lake(ui)==1)then + Qin(i)=Qin(i)+Q_res(ui) + else + Qin(i)=Qin(i)+Qout(ui) + endif + enddo + endif + + ! Update water storage in the current node + Wr(i)=max(0.D0,Cl(i)+C1(i)+C2(i)) + A1(i)=Qin(i)-C1(i)/dt-C2(i)/dt + Qout(i)=max(0.D0,Al(i)+A1(i)) + + ! Calculate flow parameters based on river flow characteristics + if(kr(i)>small.and.Qin(i)+Qrunf(i)>small)then + co3=max(0.D0,(1.D0-exp(-kr(i)*dt))/kr(i)) + P(i)=(dt*Qout(i)-co3*Qr0(i))/((Qin(i)+Qrunf(i))*(dt-co3)) + if(P(i)>0.5D0.and.P(i)<1.5D0)then + B1(i)=P(i)*(Qin(i)+Qrunf(i)) + B2(i)=-P(i)*(Qin(i)+Qrunf(i))+Qr0(i) + else + B1(i)=Qout(i) + B2(i)=0.D0 + endif + else + B1(i)=Qout(i) + B2(i)=0.D0 + P(i)=-9999. + endif + + ! Call lake and reservoir calculation subroutines + call lake_cal(active_lake(i),area_lake(i),Q_lake(i),Wr_lake(i),Qout(i),B1(i),B2(i)) + call res_cal(active_res(i),active_lake(i),Qout(i),Q_lake(i),type_res(i),cat2res(i),& + Q_res(i),wid_res(i),fld_res(i),Wr_res(i),Qfld_thres(i),cap_res(i),B1(i),B2(i)) + + endif + enddo + + !$OMP END DO + !$OMP END PARALLEL + + enddo + + ! Write the output for the current time step + call write_output(nc,yr_s,mon_s,day_s,Qout,Ws,Wr,Q_res,Wr_res,Q_lake,Wr_lake) + +ENDDO + +end diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/offline/run b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/offline/run new file mode 100755 index 000000000..d8ec1089f --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/offline/run @@ -0,0 +1,22 @@ +#!/bin/bash + +module load comp/intel/2021.3.0 + +if [ $# -lt 1 ]; then + echo "no f90 specified" + exit +fi + +string=$1 +array=(${string//./ }) + +FILENAME=${array[0]} + +NETCDF_PATH=/discover/nobackup/yzeng3/apps/netcdf-4.2.1.1 +LD_LIBRARY_PATH=$NETCDF_PATH/lib:$LD_LIBRARY_PATH + +#ifort -qopenmp ncdioMod.f90 rwncMod.f90 interp_M36toPfaf.f90 river_io_mod.f90 res_mod.f90 lake_mod.f90 ${FILENAME}.f90 -I$NETCDF_PATH/include -L$NETCDF_PATH/lib -L/usr/local/intel/oneapi/2021/compiler/2021.4.0/linux/lib -lnetcdf -lnetcdff -o ${FILENAME}.out + +ifort -qopenmp ncdioMod.f90 rwncMod.f90 interp_M36toPfaf.f90 river_io_mod.f90 res_mod.f90 lake_mod.f90 ${FILENAME}.f90 -I$NETCDF_PATH/include -L$NETCDF_PATH/lib -lnetcdf -lnetcdff -o ${FILENAME}.out + +./${FILENAME}.out diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/offline/rwncMod.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/offline/rwncMod.f90 new file mode 100644 index 000000000..3b076e14a --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOSroute_GridComp/offline/rwncMod.f90 @@ -0,0 +1,516 @@ +module rwncfile + + use ncdio + implicit none + + 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 + + public :: write_ncfile_int2d + public :: write_ncfile_real2d + public :: write_ncfile_double2d + + public :: create_ncfile_byte2d + public :: create_ncfile_short2d + public :: create_ncfile_short3d + public :: create_ncfile_int3d + public :: create_ncfile_int2d + + public :: create_ncfile_long2d + public :: create_ncfile_real2d + public :: create_ncfile_real3d + public :: create_ncfile_double2d + + 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 write_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="write" + integer :: ncid, varid, omode + + call check_ret(nf_open(filename, nf_write, ncid), subname) + call check_ret(nf_set_fill(ncid, nf_nofill, omode), subname) + call ncd_ioglobal(varname=varname, data=var, ncid=ncid, flag='write') + call check_ret(nf_sync(ncid), subname) + call check_ret(nf_close(ncid), subname) + end subroutine write_ncfile_int2d +!------------------------------------------------------------------------------------------ + subroutine write_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="write" + integer :: ncid, varid, omode + + call check_ret(nf_open(filename, nf_write, ncid), subname) + call check_ret(nf_set_fill(ncid, nf_nofill, omode), subname) + call ncd_ioglobal(varname=varname, data=var, ncid=ncid, flag='write') + call check_ret(nf_sync(ncid), subname) + call check_ret(nf_close(ncid), subname) + end subroutine write_ncfile_real2d +!------------------------------------------------------------------------------------------ + subroutine write_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="write" + integer :: ncid, varid, omode + + call check_ret(nf_open(filename, nf_write, ncid), subname) + call check_ret(nf_set_fill(ncid, nf_nofill, omode), subname) + call ncd_ioglobal(varname=varname, data=var, ncid=ncid, flag='write') + call check_ret(nf_sync(ncid), subname) + call check_ret(nf_close(ncid), subname) + end subroutine write_ncfile_double2d +!------------------------------------------------------------------------------------------ + subroutine create_ncfile_int2d(filename,varname,var,lon,lat,nlon,nlat) + character(len=*), intent(in) :: filename + character(len=*), intent(in) :: varname + integer, intent(in) :: nlon, nlat + integer, intent(inout) :: var(nlon,nlat) + real*8, intent(in) :: lon(nlon),lat(nlat) + + character(len=4) :: subname="create" + integer :: ncid, varid, dimid + real*8 :: lon1(nlon), lat1(nlat) + + lon1=lon + lat1=lat + call check_ret(nf_create(trim(filename), nf_clobber, ncid), subname) + call check_ret(nf_def_dim(ncid,'lon',nlon, dimid), subname) + call check_ret(nf_def_dim(ncid,'lat',nlat, dimid), subname) + call ncd_defvar(ncid=ncid, varname='lon', xtype=nf_double, dim1name='lon', & + long_name='longtitude', units='degrees_east') + call ncd_defvar(ncid=ncid, varname='lat', xtype=nf_double, dim2name='lat', & + long_name='latitude', units='degrees_north') + call ncd_defvar(ncid=ncid, varname=varname, xtype=nf_int, dim1name='lon', & + dim2name='lat', long_name=varname, units='unitless', fill_value=-9999.) + call check_ret(nf_enddef(ncid), subname) + call ncd_ioglobal(varname='lon', data=lon1, flag='write',ncid=ncid) + call ncd_ioglobal(varname='lat', data=lat1, flag='write',ncid=ncid) + call ncd_ioglobal(varname=varname, data=var, flag='write',ncid=ncid) + call check_ret(nf_close(ncid), subname) + end subroutine create_ncfile_int2d + + subroutine create_ncfile_long2d(filename,varname,var,lon,lat,nlon,nlat) + character(len=*), intent(in) :: filename + character(len=*), intent(in) :: varname + integer, intent(in) :: nlon, nlat + integer*8, intent(inout) :: var(nlon,nlat) + real*8, intent(in) :: lon(nlon),lat(nlat) + + character(len=4) :: subname="create" + integer :: ncid, varid, dimid + real*8 :: lon1(nlon), lat1(nlat) + + lon1=lon + lat1=lat + call check_ret(nf_create(trim(filename), NF_NETCDF4, ncid), subname) + call check_ret(nf_def_dim(ncid,'lon',nlon, dimid), subname) + call check_ret(nf_def_dim(ncid,'lat',nlat, dimid), subname) + call ncd_defvar(ncid=ncid, varname='lon', xtype=nf_double, dim1name='lon',& + long_name='longtitude', units='degrees_east') + call ncd_defvar(ncid=ncid, varname='lat', xtype=nf_double, dim2name='lat',& + long_name='latitude', units='degrees_north') + call ncd_defvar(ncid=ncid, varname=varname, xtype=nf_int64, dim1name='lon',& + dim2name='lat', long_name=varname, units='unitless',fill_value=-9999.) + call check_ret(nf_enddef(ncid), subname) + call ncd_ioglobal(varname='lon', data=lon1, flag='write',ncid=ncid) + call ncd_ioglobal(varname='lat', data=lat1, flag='write',ncid=ncid) + call ncd_ioglobal(varname=varname, data=var, flag='write',ncid=ncid) + call check_ret(nf_close(ncid), subname) + end subroutine create_ncfile_long2d + +!------------------------------------------------------------------------------------------ + subroutine create_ncfile_byte2d(filename,varname,var,lon,lat,nlon,nlat) + character(len=*), intent(in) :: filename + character(len=*), intent(in) :: varname + integer, intent(in) :: nlon, nlat + byte, intent(inout) :: var(nlon,nlat) + real*8, intent(in) :: lon(nlon),lat(nlat) + + character(len=4) :: subname="create" + integer :: ncid, varid, dimid + real*8 :: lon1(nlon), lat1(nlat) + + lon1=lon + lat1=lat + call check_ret(nf_create(trim(filename), nf_clobber, ncid), subname) + call check_ret(nf_def_dim(ncid,'lon',nlon, dimid), subname) + call check_ret(nf_def_dim(ncid,'lat',nlat, dimid), subname) + call ncd_defvar(ncid=ncid, varname='lon', xtype=nf_double, dim1name='lon', & + long_name='longtitude', units='degrees_east') + call ncd_defvar(ncid=ncid, varname='lat', xtype=nf_double, dim2name='lat', & + long_name='latitude', units='degrees_north') + call ncd_defvar(ncid=ncid, varname=varname, xtype=nf_byte, dim1name='lon', & + dim2name='lat', long_name=varname, units='unitless',fill_value=-128. ) + call check_ret(nf_enddef(ncid), subname) + call ncd_ioglobal(varname='lon', data=lon1, flag='write',ncid=ncid) + call ncd_ioglobal(varname='lat', data=lat1, flag='write',ncid=ncid) + call ncd_ioglobal(varname=varname, data=var, flag='write',ncid=ncid) + call check_ret(nf_close(ncid), subname) + end subroutine create_ncfile_byte2d + +!------------------------------------------------------------------------------------------ + subroutine create_ncfile_short2d(filename,varname,var,lon,lat,nlon,nlat) + character(len=*), intent(in) :: filename + character(len=*), intent(in) :: varname + integer, intent(in) :: nlon, nlat + integer*2, intent(inout) :: var(nlon,nlat) + real*8, intent(in) :: lon(nlon),lat(nlat) + + character(len=4) :: subname="create" + integer :: ncid, varid, dimid + real*8 :: lon1(nlon), lat1(nlat) + + lon1=lon + lat1=lat + call check_ret(nf_create(trim(filename), nf_clobber, ncid), subname) + call check_ret(nf_def_dim(ncid,'lon',nlon, dimid), subname) + call check_ret(nf_def_dim(ncid,'lat',nlat, dimid), subname) + call ncd_defvar(ncid=ncid, varname='lon', xtype=nf_double, dim1name='lon', & + long_name='longtitude', units='degrees_east') + call ncd_defvar(ncid=ncid, varname='lat', xtype=nf_double, dim2name='lat', & + long_name='latitude', units='degrees_north') + call ncd_defvar(ncid=ncid, varname=varname, xtype=nf_short, dim1name='lon', & + dim2name='lat', long_name=varname, units='unitless',fill_value=-9999. ) + call check_ret(nf_enddef(ncid), subname) + call ncd_ioglobal(varname='lon', data=lon1, flag='write',ncid=ncid) + call ncd_ioglobal(varname='lat', data=lat1, flag='write',ncid=ncid) + call ncd_ioglobal(varname=varname, data=var, flag='write',ncid=ncid) + call check_ret(nf_close(ncid), subname) + end subroutine create_ncfile_short2d + + +!------------------------------------------------------------------------------------------ + subroutine create_ncfile_real2d(filename,varname,var,lon,lat,nlon,nlat) + character(len=*), intent(in) :: filename + character(len=*), intent(in) :: varname + integer, intent(in) :: nlon, nlat + real, intent(inout) :: var(nlon,nlat) + real*8, intent(in) :: lon(nlon),lat(nlat) + + character(len=4) :: subname="create" + integer :: ncid, varid, dimid + real*8 :: lon1(nlon), lat1(nlat) + + lon1=lon + lat1=lat + call check_ret(nf_create(trim(filename), nf_clobber, ncid), subname) + call check_ret(nf_def_dim(ncid,'lon',nlon, dimid), subname) + call check_ret(nf_def_dim(ncid,'lat',nlat, dimid), subname) + call ncd_defvar(ncid=ncid, varname='lon', xtype=nf_double, dim1name='lon', & + long_name='longtitude', units='degrees_east') + call ncd_defvar(ncid=ncid, varname='lat', xtype=nf_double, dim2name='lat', & + long_name='latitude', units='degrees_north') + call ncd_defvar(ncid=ncid, varname=varname, xtype=nf_float, dim1name='lon', & + dim2name='lat', long_name=varname, units='unitless', fill_value=-9999.) + call check_ret(nf_enddef(ncid), subname) + call ncd_ioglobal(varname='lon', data=lon1, flag='write',ncid=ncid) + call ncd_ioglobal(varname='lat', data=lat1, flag='write',ncid=ncid) + call ncd_ioglobal(varname=varname, data=var, flag='write',ncid=ncid) + call check_ret(nf_close(ncid), subname) + end subroutine create_ncfile_real2d + +!------------------------------------------------------------------------------------------ + subroutine create_ncfile_short3d(filename,varname,var,lon,lat,lev,nlon,nlat,nlev) + character(len=*), intent(in) :: filename + character(len=*), intent(in) :: varname + integer, intent(in) :: nlon, nlat, nlev + integer*2, intent(inout) :: var(nlon,nlat,nlev) + real*8, intent(in) :: lon(nlon),lat(nlat),lev(nlev) + + character(len=4) :: subname="create" + integer :: ncid, varid, dimid + real*8 :: lon1(nlon), lat1(nlat), lev1(nlev) + + lon1=lon + lat1=lat + lev1=lev + call check_ret(nf_create(trim(filename), nf_clobber, ncid), subname) + call check_ret(nf_def_dim(ncid,'lon',nlon, dimid), subname) + call check_ret(nf_def_dim(ncid,'lat',nlat, dimid), subname) + call check_ret(nf_def_dim(ncid,'lev',nlev, dimid), subname) + + call ncd_defvar(ncid=ncid, varname='lon', xtype=nf_double, dim1name='lon', & + long_name='longtitude', units='degrees_east') + call ncd_defvar(ncid=ncid, varname='lat', xtype=nf_double, dim2name='lat', & + long_name='latitude', units='degrees_north') + call ncd_defvar(ncid=ncid, varname='lev', xtype=nf_double, dim2name='lev', & + long_name='level', units='unitless') + + call ncd_defvar(ncid=ncid, varname=varname, xtype=nf_short, dim1name='lon', & + dim2name='lat', dim3name='lev', long_name=varname, units='unitless', fill_value=-9999.) + call check_ret(nf_enddef(ncid), subname) + call ncd_ioglobal(varname='lon', data=lon1, flag='write',ncid=ncid) + call ncd_ioglobal(varname='lat', data=lat1, flag='write',ncid=ncid) + call ncd_ioglobal(varname='lev', data=lev1, flag='write',ncid=ncid) + call ncd_ioglobal(varname=varname, data=var, flag='write',ncid=ncid) + call check_ret(nf_close(ncid), subname) + end subroutine create_ncfile_short3d +!------------------------------------------------------------------------------------------ + subroutine create_ncfile_int3d(filename,varname,var,lon,lat,lev,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) + real*8, intent(in) :: lon(nlon),lat(nlat),lev(nlev) + + character(len=4) :: subname="create" + integer :: ncid, varid, dimid + real*8 :: lon1(nlon), lat1(nlat), lev1(nlev) + + lon1=lon + lat1=lat + lev1=lev + call check_ret(nf_create(trim(filename), nf_clobber, ncid), subname) + call check_ret(nf_def_dim(ncid,'lon',nlon, dimid), subname) + call check_ret(nf_def_dim(ncid,'lat',nlat, dimid), subname) + call check_ret(nf_def_dim(ncid,'lev',nlev, dimid), subname) + + call ncd_defvar(ncid=ncid, varname='lon', xtype=nf_double, dim1name='lon', & + long_name='longtitude', units='degrees_east') + call ncd_defvar(ncid=ncid, varname='lat', xtype=nf_double, dim2name='lat', & + long_name='latitude', units='degrees_north') + call ncd_defvar(ncid=ncid, varname='lev', xtype=nf_double, dim2name='lev', & + long_name='level', units='unitless') + + call ncd_defvar(ncid=ncid, varname=varname, xtype=nf_int, dim1name='lon', & + dim2name='lat', dim3name='lev', long_name=varname, units='unitless', fill_value=-9999.) + call check_ret(nf_enddef(ncid), subname) + call ncd_ioglobal(varname='lon', data=lon1, flag='write',ncid=ncid) + call ncd_ioglobal(varname='lat', data=lat1, flag='write',ncid=ncid) + call ncd_ioglobal(varname='lev', data=lev1, flag='write',ncid=ncid) + call ncd_ioglobal(varname=varname, data=var, flag='write',ncid=ncid) + call check_ret(nf_close(ncid), subname) + end subroutine create_ncfile_int3d +!------------------------------------------------------------------------------------------ + subroutine create_ncfile_real3d(filename,varname,var,lon,lat,lev,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) + real*8, intent(in) :: lon(nlon),lat(nlat),lev(nlev) + + character(len=4) :: subname="create" + integer :: ncid, varid, dimid + real*8 :: lon1(nlon), lat1(nlat), lev1(nlev) + + lon1=lon + lat1=lat + lev1=lev + call check_ret(nf_create(trim(filename), nf_clobber, ncid), subname) + call check_ret(nf_def_dim(ncid,'lon',nlon, dimid), subname) + call check_ret(nf_def_dim(ncid,'lat',nlat, dimid), subname) + call check_ret(nf_def_dim(ncid,'lev',nlev, dimid), subname) + + call ncd_defvar(ncid=ncid, varname='lon', xtype=nf_double, dim1name='lon', & + long_name='longtitude', units='degrees_east') + call ncd_defvar(ncid=ncid, varname='lat', xtype=nf_double, dim2name='lat', & + long_name='latitude', units='degrees_north') + call ncd_defvar(ncid=ncid, varname='lev', xtype=nf_double, dim2name='lev', & + long_name='level', units='unitless') + + call ncd_defvar(ncid=ncid, varname=varname, xtype=nf_float, dim1name='lon', & + dim2name='lat', dim3name='lev', long_name=varname, units='unitless', fill_value=-9999.) + call check_ret(nf_enddef(ncid), subname) + call ncd_ioglobal(varname='lon', data=lon1, flag='write',ncid=ncid) + call ncd_ioglobal(varname='lat', data=lat1, flag='write',ncid=ncid) + call ncd_ioglobal(varname='lev', data=lev1, flag='write',ncid=ncid) + call ncd_ioglobal(varname=varname, data=var, flag='write',ncid=ncid) + call check_ret(nf_close(ncid), subname) + end subroutine create_ncfile_real3d + +!------------------------------------------------------------------------------------------ + subroutine create_ncfile_double2d(filename,varname,var,lon,lat,nlon,nlat) + character(len=*), intent(in) :: filename + character(len=*), intent(in) :: varname + integer, intent(in) :: nlon, nlat + real*8, intent(inout) :: var(nlon,nlat) + real*8, intent(in) :: lon(nlon),lat(nlat) + + character(len=4) :: subname="create" + integer :: ncid, varid, dimid + real*8 :: lon1(nlon), lat1(nlat) + + lon1=lon + lat1=lat + call check_ret(nf_create(trim(filename), nf_clobber, ncid), subname) + call check_ret(nf_def_dim(ncid,'lon',nlon, dimid), subname) + call check_ret(nf_def_dim(ncid,'lat',nlat, dimid), subname) + call ncd_defvar(ncid=ncid, varname='lon', xtype=nf_double, dim1name='lon', & + long_name='longtitude', units='degrees_east') + call ncd_defvar(ncid=ncid, varname='lat', xtype=nf_double, dim2name='lat', & + long_name='latitude', units='degrees_north') + call ncd_defvar(ncid=ncid, varname=varname, xtype=nf_double, dim1name='lon', & + dim2name='lat', long_name=varname, units='unitless', fill_value=-9999.) + call check_ret(nf_enddef(ncid), subname) + call ncd_ioglobal(varname='lon', data=lon1, flag='write',ncid=ncid) + call ncd_ioglobal(varname='lat', data=lat1, flag='write',ncid=ncid) + call ncd_ioglobal(varname=varname, data=var, flag='write',ncid=ncid) + call check_ret(nf_close(ncid), subname) + end subroutine create_ncfile_double2d +!------------------------------------------------------------------------------------------ +end module rwncfile + 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 index 4d0f6a2da..922a17e14 100644 --- 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 @@ -3,7 +3,7 @@ MODULE routing_model IMPLICIT NONE private - public :: river_routing, SEARCH_DNST, ROUTE_DT + public :: river_routing_lin, river_routing_hyd, SEARCH_DNST, ROUTE_DT integer , parameter :: ROUTE_DT = 3600 CONTAINS @@ -11,7 +11,7 @@ MODULE routing_model ! ------------------------------------------------------------------------ - SUBROUTINE RIVER_ROUTING ( & + SUBROUTINE RIVER_ROUTING_LIN ( & NCAT, & RUNCATCH,AREACAT,LENGSC, & WSTREAM,WRIVER, & @@ -24,8 +24,9 @@ SUBROUTINE RIVER_ROUTING ( & REAL, INTENT(OUT), DIMENSION (NCAT) :: QSFLOW,QOUTFLOW REAL, PARAMETER :: K_SIMPLE = 0.111902, K_RES_MAX = 0.8 ! m1_r2com_c1 + REAL, PARAMETER :: CUR_AVG = 1.4 REAL, PARAMETER :: P1 = 0.010611, P2 = 0.188556, P3 = 0.096864, & - P4 = 0.691310, P5 = 0.365747, P6 = 0.009831 ! m5_calib_240 + P4 = 0.691310, P5 = 0.1, P6 = 0.009831 ! m5_calib_240, ori P5 = 0.365747, INTEGER :: N,I,J REAL :: COEFF, LS, COEFF1, COEFF2,ROFF @@ -58,7 +59,7 @@ SUBROUTINE RIVER_ROUTING ( & ! Updating WSTREAM WSTREAM(N) = WSTREAM(N) + RUNCATCH(N) * REAL (ROUTE_DT) - LS = AREACAT(N) / (AMAX1(1.,LENGSC (N))) + LS = AREACAT(N) / (AMAX1(1.,LENGSC (N))) /4. * CUR_AVG ROFF = RUNCATCH(N) * AREACAT(N) IF(ROFF < 2. ) THEN COEFF = RESCONST (LS, P1, P2) @@ -84,6 +85,7 @@ SUBROUTINE RIVER_ROUTING ( & IF(COEFF > K_RES_MAX) COEFF = K_SIMPLE QOUTFLOW(N) = COEFF * WRIVER(N) + QOUTFLOW(N) = MIN(QOUTFLOW(N), WRIVER(N)) !make WRIVER(N) >=0. WRIVER(N) = WRIVER(N) - QOUTFLOW(N) QOUTFLOW(N) = QOUTFLOW(N) / REAL (ROUTE_DT) @@ -91,7 +93,7 @@ SUBROUTINE RIVER_ROUTING ( & RETURN - END SUBROUTINE RIVER_ROUTING + END SUBROUTINE RIVER_ROUTING_LIN ! ------------------------------------------------------------------------------------------------------- @@ -136,4 +138,91 @@ END SUBROUTINE SEARCH_DNST ! ------------------------------------------------------------------------------------------------------- + SUBROUTINE RIVER_ROUTING_HYD ( & + NCAT, & + Qrunf0,llc_ori,lstr, & + qstr_clmt0, qri_clmt0, qin_clmt0, & + K, Kstr0, & + Ws0,Wr0, & + Qs,Qout) + + IMPLICIT NONE + INTEGER, INTENT(IN) :: NCAT + 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.025 ! Factor for local stream scaling + real, parameter :: M = 0.45 ! Parameter in hydraulic geometry formula + real, parameter :: mm = 0.35 ! Parameter in hydraulic geometry formula + real, parameter :: rho = 1000. + real, parameter :: cur_avg = 1.4 + + 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 + + integer :: i,j + + + 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 + + ! Calculate 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