diff --git a/README.md b/README.md index b8993628..b89e1ace 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.3953155.svg)](https://doi.org/10.5281/zenodo.3953155) +[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.4395155.svg)](https://doi.org/10.5281/zenodo.4395155) [![Documentation Status](https://readthedocs.org/projects/mizuroute/badge/?version=master)](https://mizuroute.readthedocs.io/en/master/?badge=master) # mizuRoute diff --git a/docs/source/Control_file.rst b/docs/source/Control_file.rst index e8985dc5..0603a918 100644 --- a/docs/source/Control_file.rst +++ b/docs/source/Control_file.rst @@ -1,10 +1,10 @@ Control file ============ -Control file is a simple text file, mainly defining model control such as simulation time, file name and locations, routing option etc. +Control file is a simple text file, defining various model controls such as simulation time, file names and locations, routing options etc. Variables in control file are read in the beginning of the code (see ``./build/src/read_control.f90``) and saved in fortran variable specified by tag (inside <> in table) and as public variables (see ``./build/src/public_var.f90``) . -Some of such public varialbes have some default values, but most of them are not defined. +Some of control varialbes have their default values, but most of them are not defined. Those undefined variables need to be defined in control file. Other variables in supplement table have their default values but can be also included in control file to overwrite the values. The order of variables in the control file does not matter. However, grouping variables into similar themes is recommended for readibility. @@ -19,6 +19,7 @@ Some of rules * Format: variable ! comments * tag is Fortran variable name and cannot be changed and have to be enclosed by <> * need ! after variable, otherwise getting error. +* Do not leave any lines empty in control file The following variables (not pre-defined in the code) need to be defined in control file. @@ -86,12 +87,6 @@ The following variables (not pre-defined in the code) need to be defined in cont +--------+------------------------+-------------------------------------------------------------------------------------------+ | 2,3 | | dimension name for data | +--------+------------------------+-------------------------------------------------------------------------------------------+ -| 1,2,3 | | restart ouput timing options. N[n]ever, L[l]ast, S[s]pecified. | -+--------+------------------------+-------------------------------------------------------------------------------------------+ -| 1,2,3 | | specified restart date in yyyy-mm-dd (hh:mm:ss) if = "Specified" | -+--------+------------------------+-------------------------------------------------------------------------------------------+ -| 1,2,3 | | input restart netCDF name. If not specified, simulation start with cold start | -+--------+------------------------+-------------------------------------------------------------------------------------------+ | 1,2,3 | | option for routing schemes 0-> both, 1->IRF, 2->KWT, otherwise error | +--------+------------------------+-------------------------------------------------------------------------------------------+ @@ -120,6 +115,8 @@ Variables that have default values but can be overwritten +------------------------+------------------------+--------------------------------------------------------------------------+ | | From runoff input | specified time units since yyyy-mm-dd (hh:mm:ss). See note 4 | +------------------------+------------------------+--------------------------------------------------------------------------+ +| | netcdf4 | netcdf format for output netcdf. other options: classic, 64bit_offset. | ++------------------------+------------------------+--------------------------------------------------------------------------+ 1. River network subset mode. @@ -138,10 +135,42 @@ Often case, river network data has different variable names than defaults. In th See :doc:`River parameters `. +Restart options +--------------------- + +mizuRoute does not write restart netCDF as default. The following control variables are used to control restart dropoff timing and use restart file for continuous run from the previous simulations. +The restart file is written at previous time step to the specified time. In other words, if ``Specified`` is used for and ``1981-01-01-00000`` is specified in , mizuRoute writes restart file +at ``1980-12-31 00:00:00`` for daily time step. The restart file name uses the time stamp at user specified timing. ``Annual``, ``Monthly``, ``Daily`` options also follow This convention. + +The restart file name convension: .r.yyyy-mm-dd-sssss.nc + + ++---------------------+---------------------------------------------------------------------------------------------------------+ +| tag | Description | ++=====================+=========================================================================================================+ +| | directory for restart files. defualt is | ++---------------------+---------------------------------------------------------------------------------------------------------+ +| | restart ouput options. N[n]ever (default), L[l]ast, S[s]pecified, Annual, M[m]onthly, D[d]aily. | ++---------------------+---------------------------------------------------------------------------------------------------------+ +| | restart time in yyyy-mm-dd (hh:mm:ss). required if = "Specified" | ++---------------------+---------------------------------------------------------------------------------------------------------+ +| | periodic restart month (default 1). Effective if ="Annual" | ++---------------------+---------------------------------------------------------------------------------------------------------+ +| | periodic restart day (default 1). Effective if ="Annual" or "Monthly" | ++---------------------+---------------------------------------------------------------------------------------------------------+ +| | periodic restart hour (default 0). Effective if ="Annual", "Monthly", or "Daily" | ++---------------------+---------------------------------------------------------------------------------------------------------+ +| | input restart netCDF name. If not specified, simulation start with cold start | ++---------------------+---------------------------------------------------------------------------------------------------------+ + + Output variables --------------------- The following variables, besides time, basinID (RN_hru ID) and reachID can be output in netCDF. Users can control which variables are output by setting to T or F in control file. All the variables are set to T by default. +The output file name includes a timie stamp at the first time step. + +The output file name convension: .h.yyyy-mm-dd-sssss.nc +------------------------+------------------------------------------------------------------------------------------------+ diff --git a/route/build/Makefile b/route/build/Makefile index 33d87e02..0423ebf5 100644 --- a/route/build/Makefile +++ b/route/build/Makefile @@ -120,6 +120,8 @@ DATATYPES = \ public_var.f90 \ dataTypes.f90 \ var_lookup.f90 \ + time_utils.f90 \ + datetime_data.f90 \ globalData.f90 \ popMetadat.f90 \ allocation.f90 @@ -127,9 +129,15 @@ DATATYPES = \ UTILS = \ nr_utility.f90 \ ascii_util.f90 \ - time_utils.f90 \ ncio_utils.f90 \ gamma_func.f90 +# initialization +INIT = \ + network_topo.f90 \ + process_param.f90 \ + process_ntopo.f90 \ + pfafstetter.f90 \ + domain_decomposition.f90 # read/write files IO = \ remap.f90 \ @@ -143,16 +151,9 @@ IO = \ read_restart.f90 \ write_restart.f90 \ write_simoutput.f90 -# initialization -INIT = \ - process_time.f90 \ - network_topo.f90 \ - process_param.f90 \ - process_ntopo.f90 \ - pfafstetter.f90 \ - domain_decomposition.f90 # CORE CORE = \ + model_finalize.f90 \ accum_runoff.f90 \ basinUH.f90 \ irf_route.f90 \ @@ -161,7 +162,7 @@ CORE = \ model_setup.f90 # concatanate model subroutines -TEMP_MODSUB = $(DATATYPES) $(UTILS) $(IO) $(INIT) $(CORE) +TEMP_MODSUB = $(DATATYPES) $(UTILS) $(INIT) $(IO) $(CORE) # insert appropriate directory name MODSUB = $(patsubst %, $(F_KORE_DIR)%, $(TEMP_MODSUB)) diff --git a/route/build/src/accum_runoff.f90 b/route/build/src/accum_runoff.f90 index d19fb747..8efa546f 100644 --- a/route/build/src/accum_runoff.f90 +++ b/route/build/src/accum_runoff.f90 @@ -33,7 +33,8 @@ SUBROUTINE accum_runoff(iEns, & ! input: index of runoff ensemble to be ! ! ---------------------------------------------------------------------------------------- - USE dataTypes, ONLY: subbasin_omp ! mainstem+tributary data structures + USE dataTypes, ONLY: subbasin_omp ! mainstem+tributary data structures + USE model_finalize, ONLY : handle_err implicit none ! input @@ -56,14 +57,9 @@ SUBROUTINE accum_runoff(iEns, & ! input: index of runoff ensemble to be integer(i4b) :: iTrib, ix ! loop indices logical(lgt), allocatable :: doRoute(:) ! logical to indicate which reaches are processed character(len=strLen) :: cmessage ! error message from subroutines - integer*8 :: cr ! rate - integer*8 :: startTime,endTime ! date/time for the start and end of the initialization - real(dp) :: elapsedTime ! elapsed time for the process ierr=0; message='accum_runoff/' - call system_clock(count_rate=cr) - ! check if (size(NETOPO_in)/=size(RCHFLX_out(iens,:))) then ierr=20; message=trim(message)//'sizes of NETOPO and RCHFLX mismatch'; return endif @@ -84,8 +80,6 @@ SUBROUTINE accum_runoff(iEns, & ! input: index of runoff ensemble to be nDom = size(river_basin) - call system_clock(startTime) - do ix = 1,nDom ! 1. Route tributary reaches (parallel) ! compute the sum of all upstream runoff at each point in the river network @@ -107,17 +101,13 @@ SUBROUTINE accum_runoff(iEns, & ! input: index of runoff ensemble to be if (.not. doRoute(jSeg)) cycle call accum_qupstream(iens, jSeg, ixDesire, NETOPO_in, RCHFLX_out, ierr, cmessage) - !if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + if(ierr/=0) call handle_err(ierr, trim(message)//trim(cmessage)) end do end do !$OMP END PARALLEL DO - end do ! looping through stream segments - - call system_clock(endTime) - elapsedTime = real(endTime-startTime, kind(dp))/real(cr) - !write(*,"(A,1PG15.7,A)") ' elapsed-time [routing/accum] = ', elapsedTime, ' s' + end do END SUBROUTINE accum_runoff @@ -142,13 +132,13 @@ subroutine accum_qupstream(iEns, & ! input: index of runoff ensemble to integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! Local variables to - real(dp), allocatable :: uprflux(:) ! upstream Reach fluxes + real(dp) :: q_upstream ! upstream Reach fluxes integer(i4b) :: nUps ! number of upstream segment integer(i4b) :: iUps ! upstream reach index integer(i4b) :: iRch_ups ! index of upstream reach in NETOPO + character(len=strLen) :: fmt1,fmt2 ! format string character(len=strLen) :: cmessage ! error message from subroutine - ! initialize error control ierr=0; message='accum_qupstream/' ! identify number of upstream segments of the reach being processed @@ -156,26 +146,30 @@ subroutine accum_qupstream(iEns, & ! input: index of runoff ensemble to RCHFLX_out(iEns,segIndex)%UPSTREAM_QI = RCHFLX_out(iEns,segIndex)%BASIN_QR(1) + q_upstream = 0._dp if (nUps>0) then - allocate(uprflux(nUps), stat=ierr, errmsg=cmessage) - if(ierr/=0)then; message=trim(message)//trim(cmessage)//': uprflux'; return; endif - do iUps = 1,nUps iRch_ups = NETOPO_in(segIndex)%UREACHI(iUps) ! index of upstream of segIndex-th reach - uprflux(iUps) = RCHFLX_out(iens,iRch_ups)%UPSTREAM_QI + q_upstream = q_upstream + RCHFLX_out(iens,iRch_ups)%UPSTREAM_QI end do - RCHFLX_out(iEns,segIndex)%UPSTREAM_QI = RCHFLX_out(iEns,segIndex)%UPSTREAM_QI + sum(uprflux) + RCHFLX_out(iEns,segIndex)%UPSTREAM_QI = RCHFLX_out(iEns,segIndex)%UPSTREAM_QI + q_upstream endif ! check - if(NETOPO_in(segIndex)%REACHIX == ixDesire)then - print*, 'CHECK ACCUM_RUNOFF' - print*, ' UREACHK, uprflux = ', (NETOPO_in(segIndex)%UREACHK(iUps), uprflux(iUps), iUps=1,nUps) - print*, ' RCHFLX_out(iEns,segIndex)%BASIN_QR(1) = ', RCHFLX_out(iEns,segIndex)%BASIN_QR(1) - print*, ' RCHFLX_out%UPSTREAM_QI = ', RCHFLX_out(iens,segIndex)%UPSTREAM_QI + if(segIndex == ixDesire)then + write(fmt1,'(A,I5,A)') '(A,1X',nUps,'(1X,I10))' + write(fmt2,'(A,I5,A)') '(A,1X',nUps,'(1X,F20.7))' + write(*,'(2a)') new_line('a'),'** Check upstream discharge accumulation **' + write(*,'(a,x,I10,x,I10)') ' Reach index & ID =', segIndex, NETOPO_in(segIndex)%REACHID + write(*,'(a)') ' * upstream reach index (NETOPO_in%UREACH) and discharge (uprflux) [m3/s] :' + write(*,fmt1) ' UREACHK =', (NETOPO_in(segIndex)%UREACHK(iUps), iUps=1,nUps) + write(*,fmt2) ' prflux =', (RCHFLX_out(iens,NETOPO_in(segIndex)%UREACHI(iUps))%UPSTREAM_QI, iUps=1,nUps) + write(*,'(a)') ' * local area discharge (RCHFLX_out%BASIN_QR(1)) and final discharge (RCHFLX_out%UPSTREAM_QI) [m3/s] :' + write(*,'(a,x,F15.7)') ' RCHFLX_out%BASIN_QR(1) =', RCHFLX_out(iEns,segIndex)%BASIN_QR(1) + write(*,'(a,x,F15.7)') ' RCHFLX_out%UPSTREAM_QI =', RCHFLX_out(iens,segIndex)%UPSTREAM_QI endif end subroutine accum_qupstream diff --git a/route/build/src/basinUH.f90 b/route/build/src/basinUH.f90 index 51e6d717..571d4a83 100644 --- a/route/build/src/basinUH.f90 +++ b/route/build/src/basinUH.f90 @@ -36,12 +36,8 @@ SUBROUTINE IRF_route_basin(iens, & ! input: ensemble index integer(i4b) :: iSeg ! reach loop indix logical(lgt), allocatable :: doRoute(:) ! logical to indicate which reaches are processed character(len=strLen) :: cmessage ! error message from subroutines - integer*8 :: cr ! rate - integer*8 :: startTime,endTime ! date/time for the start and end of the initialization - real(dp) :: elapsedTime ! elapsed time for the process ierr=0; message='IRF_route_basin/' - call system_clock(count_rate=cr) nSeg = size(RCHFLX_out(iens,:)) @@ -56,7 +52,6 @@ SUBROUTINE IRF_route_basin(iens, & ! input: ensemble index doRoute(:) = .true. endif - call system_clock(startTime) !$OMP PARALLEL DO schedule(dynamic,1) & !$OMP private(iSeg) & ! loop index !$OMP private(ierr, cmessage) & ! private for a given thread @@ -75,10 +70,6 @@ SUBROUTINE IRF_route_basin(iens, & ! input: ensemble index end do !$OMP END PARALLEL DO - call system_clock(endTime) - elapsedTime = real(endTime-startTime, kind(dp))/real(cr) -! write(*,"(A,1PG15.7,A)") ' elapsed-time [routing/irf_hru] = ', elapsedTime, ' s' - END SUBROUTINE IRF_route_basin diff --git a/route/build/src/dataTypes.f90 b/route/build/src/dataTypes.f90 index 5f794d7e..9fa5964e 100644 --- a/route/build/src/dataTypes.f90 +++ b/route/build/src/dataTypes.f90 @@ -2,10 +2,11 @@ module dataTypes ! used to create/save specific data types -USE nrtype, only: i4b,dp,lgt -USE nrtype, only: strLen ! string length +USE nrtype, only: i4b,i8b,dp,lgt +USE nrtype, only: strLen USE public_var, only: realMissing USE public_var, only: integerMissing +USE public_var, only: charMissing implicit none @@ -41,16 +42,6 @@ module dataTypes logical(lgt) :: varFile = .true. ! .true. if the variable should be read from a file end type var_info - ! ---------- time structures ------------------------------------------------------------------------------ - type,public :: time - integer(i4b) :: iy = integerMissing ! year - integer(i4b) :: im = integerMissing ! month - integer(i4b) :: id = integerMissing ! day - integer(i4b) :: ih = integerMissing ! hour - integer(i4b) :: imin = integerMissing ! minute - real(dp) :: dsec = realMissing ! second - endtype time - ! ---------- states structure -------------------------------------------------------------------------- ! type,public :: var @@ -69,7 +60,7 @@ module dataTypes ! ---------- output netcdf structure -------------------------------------------------------------------------- ! type,public :: nc - character(len=strLen) :: ncname = 'empty' ! netcdf name + character(len=strLen) :: ncname = charMissing ! netcdf name integer(i4b) :: ncid = integerMissing ! netcdf id integer(i4b) :: status = integerMissing ! status: 1=defined, 2=open, 3=closed end type nc @@ -144,8 +135,8 @@ module dataTypes ! data to remap runoff hru to river network hrus type, public :: remap ! information in the mapping file - integer(i4b) , allocatable :: hru_id(:) ! Id of hrus associated with river network (="hru") - integer(i4b) , allocatable :: qhru_id(:) ! Id of hrus associated with runoff simulation (="qhru") + integer(i8b) , allocatable :: hru_id(:) ! Id of hrus associated with river network (="hru") + integer(i8b) , allocatable :: qhru_id(:) ! Id of hrus associated with runoff simulation (="qhru") integer(i4b) , allocatable :: num_qhru(:) ! number of "qhru" within "hru" integer(i4b) , allocatable :: i_index(:) ! Index in the y dimension of the runoff grid (starting with 1,1 in LL corner) integer(i4b) , allocatable :: j_index(:) ! Index in the x dimension of the runoff grid (starting with 1,1 in LL corner) @@ -162,7 +153,7 @@ module dataTypes real(dp) :: time ! time variable at one time step real(dp) , allocatable :: qsim(:) ! runoff(HM_HRU) at one time step (size: nSpace(1)) real(dp) , allocatable :: qsim2D(:,:) ! runoff(x,y) at one time step (size: /nSpace(1),nSpace(2)/) - integer(i4b) , allocatable :: hru_id(:) ! id of HM_HRUs or RN_HRUs at which runoff is stored (size: nSpace(1)) + integer(i8b) , allocatable :: hru_id(:) ! id of HM_HRUs or RN_HRUs at which runoff is stored (size: nSpace(1)) integer(i4b) , allocatable :: hru_ix(:) ! Index of RN_HRUs associated with river network (used only if HM_HRUs = RN_HRUs) real(dp) , allocatable :: basinRunoff(:)! remapped river network catchment runoff (size: number of nHRU) end type runoff @@ -178,6 +169,7 @@ module dataTypes real(DP) :: UPSAREA ! upstream area (zero if headwater basin) real(DP) :: BASAREA ! local basin area real(DP) :: TOTAREA ! UPSAREA + BASAREA + real(DP) :: QTAKE ! target abstraction/injection [m3/s] real(DP) :: MINFLOW ! minimum environmental flow end type RCHPRP @@ -242,11 +234,11 @@ module dataTypes REAL(DP), allocatable :: QFUTURE_IRF(:) ! runoff volume in future time steps for IRF routing (m3/s) REAL(DP) :: BASIN_QI ! instantaneous runoff volume from the local basin (m3/s) REAL(DP) :: BASIN_QR(0:1) ! routed runoff volume from the local basin (m3/s) - REAL(DP) :: UPSBASIN_QR ! routed runoff depth from the upstream basins (m/s) REAL(DP) :: BASIN_QR_IRF(0:1) ! routed runoff volume from all the upstream basin (m3/s) REAL(DP) :: REACH_Q ! time-step average streamflow (m3/s) REAL(DP) :: REACH_Q_IRF ! time-step average streamflow (m3/s) from IRF routing REAL(DP) :: UPSTREAM_QI ! sum of upstream streamflow (m3/s) + REAL(DP) :: REACH_VOL(0:1) ! volume of water at a reach [m3] REAL(DP) :: TAKE ! average take logical(lgt) :: CHECK_IRF ! .true. if the reach is routed ENDTYPE STRFLX @@ -296,28 +288,29 @@ END MODULE dataTypes MODULE objTypes - USE nrtype, only: i4b,dp,lgt - USE nrtype, only: strLen ! string length - USE public_var, only: realMissing - USE public_var, only: integerMissing + USE nrtype, ONLY: i4b,dp,lgt + USE nrtype, ONLY: strLen + USE public_var, ONLY: realMissing + USE public_var, ONLY: integerMissing + USE public_var, ONLY: charMissing ! define derived type for model variables, including name, description, and units - type, public :: var_info_new - character(len=strLen) :: varName = 'empty' ! variable name - character(len=strLen) :: varDesc = 'empty' ! variable description - character(len=strLen) :: varUnit = 'empty' ! variable units + type, public :: meta_var + character(len=strLen) :: varName = charMissing ! variable name + character(len=strLen) :: varDesc = charMissing ! variable description + character(len=strLen) :: varUnit = charMissing ! variable units integer(i4b) :: varType = integerMissing ! variable type integer(i4b),allocatable :: varDim(:) ! dimension ID associated with variable logical(lgt) :: varFile = .true. ! .true. if the variable should be read from a file CONTAINS procedure, pass :: init - end type var_info_new + end type meta_var CONTAINS SUBROUTINE init(this, vName, vDesc, vUnit, vType, vDim, vFile) implicit none - class(var_info_new) :: this + class(meta_var) :: this character(*), intent(in) :: vName ! variable name character(*), intent(in) :: vDesc ! variable description character(*), intent(in) :: vUnit ! variable units @@ -332,7 +325,7 @@ SUBROUTINE init(this, vName, vDesc, vUnit, vType, vDim, vFile) this%varDesc = vDesc this%varUnit = vUnit this%varType = vType - this%varDim(1:n) = vDim(1:n) + this%varDim(1:n) = vDim(1:n) this%varFile = vFile END SUBROUTINE init diff --git a/route/build/src/datetime_data.f90 b/route/build/src/datetime_data.f90 new file mode 100644 index 00000000..990b9db3 --- /dev/null +++ b/route/build/src/datetime_data.f90 @@ -0,0 +1,438 @@ +MODULE date_time + +! datetime class + +USE nrtype +USE public_var, ONLY: realMissing, integerMissing +USE public_var, ONLY: secprday, hr_per_day +USE time_utils_module, ONLY: extractTime ! +USE time_utils_module, ONLY: compJulday,& ! compute julian day + compJulday_noleap ! compute julian day for noleap calendar +USE time_utils_module, ONLY: compCalday,& ! compute calendar date and time + compCalday_noleap ! compute calendar date and time for noleap calendar + +implicit none + +type, public :: datetime + + private + integer(i4b) :: iy = integerMissing + integer(i4b) :: im = integerMissing + integer(i4b) :: id = integerMissing + integer(i4b) :: ih = integerMissing + integer(i4b) :: imin = integerMissing + real(dp) :: dsec = realMissing + +CONTAINS + + procedure, public :: set_datetime => sub_set_datetime + procedure, public :: jul2datetime => sub_jul2datetime + procedure, public :: str2datetime => sub_str2datetime + procedure, public :: year => fn_get_year + procedure, public :: month => fn_get_month + procedure, public :: day => fn_get_day + procedure, public :: hour => fn_get_hour + procedure, public :: minute => fn_get_min + procedure, public :: sec => fn_get_sec + procedure, public :: is_leap_year => fn_is_leap_year + procedure, public :: ndays_month => fn_ndays_month + procedure, public :: julianday => sub_julian_day + procedure, public :: add_mon => fn_add_months + procedure, public :: add_day => fn_add_days + procedure, public :: add_hr => fn_add_hours + procedure, public :: add_sec => fn_add_sec + procedure, public :: is_equal_mon => fn_is_equal_month + procedure, public :: is_equal_day => fn_is_equal_day + procedure, public :: is_equal_time => fn_is_equal_time + + procedure, private :: fn_is_equal + procedure, private :: fn_is_gt + procedure, private :: fn_is_ge + procedure, private :: fn_is_lt + procedure, private :: fn_is_le + procedure, private :: sub_assign + generic :: operator(==) => fn_is_equal + generic :: operator(>) => fn_is_gt + generic :: operator(>=) => fn_is_ge + generic :: operator(<) => fn_is_lt + generic :: operator(<=) => fn_is_le + generic :: assignment(=) => sub_assign + +end type datetime + +private :: sub_set_datetime, sub_jul2datetime, sub_str2datetime +private :: fn_get_year, fn_get_month, fn_get_day, fn_get_hour, fn_get_min, fn_get_sec +private :: fn_is_leap_year, sub_julian_day, fn_ndays_month +private :: fn_add_months, fn_add_days, fn_add_hours, fn_add_sec +private :: fn_is_equal_month, fn_is_equal_day, fn_is_equal_time + +CONTAINS + + SUBROUTINE sub_set_datetime(this, iy, im, id, ih, imin, dsec) + + implicit none + class(datetime) :: this + integer(i4b), intent(in) :: iy + integer(i4b), intent(in) :: im + integer(i4b), intent(in) :: id + integer(i4b), intent(in) :: ih + integer(i4b), intent(in) :: imin + real(dp), intent(in) :: dsec + + this%iy = iy + this%im = im + this%id = id + this%ih = ih + this%imin = imin + this%dsec = dsec + + END SUBROUTINE sub_set_datetime + + SUBROUTINE sub_jul2datetime(this, julday, calendar, ierr, message) + + implicit none + class(datetime) :: this + real(dp), intent(in) :: julday + character(*), intent(in) :: calendar + integer(i4b), intent(out) :: ierr + character(len=strLen),intent(out) :: message + ! local variable + character(len=strLen) :: cmessage + + ierr=0; message='sub_jul2datetime/' + + select case(trim(calendar)) + case ('noleap','365_day') + call compCalday_noleap(julday,this%iy,this%im,this%id,this%ih,this%imin,this%dsec, ierr, cmessage) + case ('standard','gregorian','proleptic_gregorian') + call compCalday(julday,this%iy,this%im,this%id,this%ih,this%imin,this%dsec, ierr, cmessage) + case default; ierr=20; message=trim(message)//'calendar name: '//trim(calendar)//' invalid'; return + end select + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + + END SUBROUTINE sub_jul2datetime + + SUBROUTINE sub_str2datetime(this, str, ierr, message) + + implicit none + class(datetime) :: this + character(*), intent(in) :: str + integer(i4b), intent(out) :: ierr + character(len=strLen),intent(out) :: message + ! local variable + character(len=strLen) :: cmessage + + ierr=0; message='sub_str2datetime/' + + call extractTime(str,this%iy,this%im,this%id,this%ih,this%imin,this%dsec, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + + END SUBROUTINE sub_str2datetime + + SUBROUTINE sub_assign(this, that) + implicit none + class(datetime),intent(out) :: this + type(datetime), intent(in) :: that + this%iy = that%iy + this%im = that%im + this%id = that%id + this%ih = that%ih + this%imin = that%imin + this%dsec = that%dsec + END SUBROUTINE sub_assign + + + integer(i4b) FUNCTION fn_get_year(this) + implicit none + class(datetime), intent(in) :: this + fn_get_year = this%iy + END FUNCTION fn_get_year + + integer(i4b) FUNCTION fn_get_month(this) + implicit none + class(datetime), intent(in) :: this + fn_get_month = this%im + END FUNCTION fn_get_month + + integer(i4b) FUNCTION fn_get_day(this) + implicit none + class(datetime), intent(in) :: this + fn_get_day = this%id + END FUNCTION fn_get_day + + integer(i4b) FUNCTION fn_get_hour(this) + implicit none + class(datetime), intent(in) :: this + fn_get_hour = this%ih + END FUNCTION fn_get_hour + + integer(i4b) FUNCTION fn_get_min(this) + implicit none + class(datetime), intent(in) :: this + fn_get_min = this%imin + END FUNCTION fn_get_min + + real(dp) FUNCTION fn_get_sec(this) + implicit none + class(datetime), intent(in) :: this + fn_get_sec = this%dsec + END FUNCTION fn_get_sec + + + logical(lgt) FUNCTION fn_is_leap_year(this) + implicit none + class(datetime), intent(in) :: this + if (mod(this%iy, 4) == 0) then + if (mod(this%iy, 100) == 0) then + if (mod(this%iy, 400) == 0) then + fn_is_leap_year = .True. + else + fn_is_leap_year = .False. + end if + else + fn_is_leap_year = .True. + end if + else + fn_is_leap_year = .False. + end if + END FUNCTION fn_is_leap_year + + integer(i4b) FUNCTION fn_ndays_month(this, calendar, ierr, message) + implicit none + class(datetime), intent(in) :: this + character(*), intent(in) :: calendar + integer(i4b), intent(out) :: ierr + character(len=strLen),intent(out) :: message + ! local variables + integer(i4b) :: nmonths(12) + + ierr=0; message="ndays_month/" + + fn_ndays_month = integerMissing + select case(trim(calendar)) + case ('standard','gregorian','proleptic_gregorian') + if ( fn_is_leap_year(this)) then + nmonths = [31,29,31,30,31,30,31,31,30,31,30,31] + else + nmonths = [31,28,31,30,31,30,31,31,30,31,30,31] + end if + case('noleap') + nmonths = [31,28,31,30,31,30,31,31,30,31,30,31] + case default; ierr=20; message=trim(message)//'calendar name: '//trim(calendar)//' invalid'; return + end select + fn_ndays_month = nmonths(this%im) + + END FUNCTION fn_ndays_month + + + SUBROUTINE sub_julian_day(this, calendar, julianDay, ierr, message) + implicit none + class(datetime), intent(in) :: this + character(*), intent(in) :: calendar + real(dp), intent(out) :: julianDay + integer(i4b), intent(out) :: ierr + character(len=strLen),intent(out) :: message + ! local variables + character(len=strLen) :: cmessage + + ierr=0; message='sub_julian_day/' + select case(trim(calendar)) + case ('standard','gregorian','proleptic_gregorian') + call compJulday(this%iy, this%im, this%id, this%ih, this%imin, this%dsec, julianDay, ierr, cmessage) + case('noleap') + call compJulday_noleap(this%iy, this%im, this%id, this%ih, this%imin, this%dsec, julianDay, ierr, cmessage) + case default; ierr=20; message=trim(message)//'calendar name: '//trim(calendar)//' invalid'; return + end select + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + END SUBROUTINE sub_julian_day + + ! ----- + ! increment(+)/decrement(+) datetime + ! ----- + + type(datetime) FUNCTION fn_add_months(this, nmonths) + implicit none + class(datetime), intent(inout) :: this + integer(i4b), intent(in) :: nmonths + fn_add_months%iy = this%iy + (nmonths/12_i4b) + fn_add_months%im = this%im + mod(nmonths,12) + fn_add_months%id = this%id + fn_add_months%ih = this%ih + fn_add_months%imin = this%imin + fn_add_months%dsec = this%dsec + END FUNCTION fn_add_months + + type(datetime) FUNCTION fn_add_days(this, days, calendar, ierr, message) + implicit none + class(datetime), intent(inout) :: this + integer(i4b), intent(in) :: days + character(*), intent(in) :: calendar + integer(i4b), intent(out) :: ierr + character(len=strLen),intent(out) :: message + ! local variables + real(dp) :: julday + character(len=strLen) :: cmessage + + ierr=0; message='fn_add_days/' + + call sub_julian_day(this, calendar, julday, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); endif + + julday = julday + real(days,dp) + call sub_jul2datetime(fn_add_days, julday, calendar, ierr, message) + if(ierr/=0)then; message=trim(message)//trim(cmessage); endif + + END FUNCTION fn_add_days + + type(datetime) FUNCTION fn_add_hours(this, hrs, calendar, ierr, message) + implicit none + class(datetime), intent(in) :: this + integer(i4b), intent(in) :: hrs + character(*), intent(in) :: calendar + integer(i4b), intent(out) :: ierr + character(len=strLen),intent(out) :: message + ! local variables + real(dp) :: julday + character(len=strLen) :: cmessage + + ierr=0; message='fn_add_hours/' + call sub_julian_day(this, calendar, julday, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); endif + + julday = julday + real(hrs,dp)/hr_per_day + call sub_jul2datetime(fn_add_hours, julday, calendar, ierr, message) + if(ierr/=0)then; message=trim(message)//trim(cmessage); endif + + END FUNCTION fn_add_hours + + type(datetime) FUNCTION fn_add_sec(this, sec, calendar, ierr, message) + implicit none + class(datetime), intent(in) :: this + real(dp), intent(in) :: sec + character(*), intent(in) :: calendar + integer(i4b), intent(out) :: ierr + character(len=strLen),intent(out) :: message + ! local variables + real(dp) :: julday + character(len=strLen) :: cmessage + + ierr=0; message='fn_add_sec/' + call sub_julian_day(this, calendar, julday, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); endif + + julday = julday + sec/secprday + call sub_jul2datetime(fn_add_sec, julday, calendar, ierr, message) + if(ierr/=0)then; message=trim(message)//trim(cmessage); endif + + END FUNCTION fn_add_sec + + ! ----- + ! check relational logic between two datetime + ! ----- + + logical(lgt) FUNCTION fn_is_equal_month(this, that) + implicit none + class(datetime), intent(in) :: this + class(datetime), intent(in) :: that + fn_is_equal_month = (this%im==that%im) + END FUNCTION fn_is_equal_month + + logical(lgt) FUNCTION fn_is_equal_day(this, that) + implicit none + class(datetime), intent(in) :: this + class(datetime), intent(in) :: that + fn_is_equal_day = (this%id==that%id) + END FUNCTION fn_is_equal_day + + logical(lgt) FUNCTION fn_is_equal_time(this, that) + implicit none + class(datetime), intent(in) :: this + class(datetime), intent(in) :: that + fn_is_equal_time = (this%ih==that%ih .and. this%imin==that%imin .and. abs(this%dsec-that%dsec)<=epsilon(this%dsec)) + END FUNCTION fn_is_equal_time + + logical(lgt) FUNCTION fn_is_equal(this, that) + ! if this == that T, otherwise F + implicit none + class(datetime), intent(in) :: this + class(datetime), intent(in) :: that + fn_is_equal = (this%iy==that%iy .and. this%im==that%im .and. this%id==that%id .and. & + this%ih==that%ih .and. this%imin==that%imin .and. abs(this%dsec-that%dsec)<=epsilon(this%dsec)) + END FUNCTION fn_is_equal + + logical(lgt) FUNCTION fn_is_gt(this, that) + ! if this >= that T, otherwise F + implicit none + class(datetime), intent(in) :: this + class(datetime), intent(in) :: that + if (this%iy > that%iy) then + fn_is_gt = .true. + else if (this%iy < that%iy) then + fn_is_gt = .false. + else + if (this%im > that%im) then + fn_is_gt = .true. + else if (this%im < that%im) then + fn_is_gt = .false. + else + if (this%id > that%id) then + fn_is_gt = .true. + else if (this%id < that%id) then + fn_is_gt = .false. + else + if (this%ih > that%ih) then + fn_is_gt = .true. + else if (this%ih < that%ih) then + fn_is_gt = .false. + else + if (this%imin > that%imin) then + fn_is_gt = .true. + else if (this%imin < that%imin) then + fn_is_gt = .false. + else + if (this%dsec > that%dsec) then + fn_is_gt = .true. + else + fn_is_gt = .false. + endif ! dsec + endif ! min + endif ! hr + endif ! day + endif ! mon + endif ! yr + END FUNCTION fn_is_gt + + logical(lgt) FUNCTION fn_is_ge(this, that) + ! if this >= that T, otherwise F + implicit none + class(datetime), intent(in) :: this + class(datetime), intent(in) :: that + fn_is_ge = .false. + if ( fn_is_gt(this, that) .or. fn_is_equal(this, that) ) then + fn_is_ge = .true. + end if + END FUNCTION fn_is_ge + + logical(lgt) FUNCTION fn_is_lt(this, that) + ! if this < that T, otherwise F + implicit none + class(datetime), intent(in) :: this + class(datetime), intent(in) :: that + fn_is_lt = .false. + if (.not.( fn_is_ge(this, that))) then + fn_is_lt = .true. + end if + END FUNCTION fn_is_lt + + logical(lgt) FUNCTION fn_is_le(this, that) + ! if this <= that T, otherwise F + implicit none + class(datetime), intent(in) :: this + class(datetime), intent(in) :: that + fn_is_le = .false. + if (.not.(fn_is_gt(this, that))) then + fn_is_le = .true. + end if + END FUNCTION fn_is_le + +END MODULE date_time diff --git a/route/build/src/globalData.f90 b/route/build/src/globalData.f90 index d6fdd055..73c2390d 100644 --- a/route/build/src/globalData.f90 +++ b/route/build/src/globalData.f90 @@ -10,7 +10,7 @@ module globalData use dataTypes, only : struct_info ! metadata type use dataTypes, only : dim_info ! metadata type use dataTypes, only : var_info ! metadata type - use objTypes, only : var_info_new ! metadata type + use objTypes, only : meta_var ! metadata type ! parameter structures USE dataTypes, only : RCHPRP ! Reach parameters (properties) @@ -33,7 +33,7 @@ module globalData use dataTypes, only : subbasin_omp ! mainstem+tributary data structures ! time data structure - use dataTypes, only : time ! time data + use date_time, only : datetime ! time data ! time data structure use dataTypes, only : nc ! netCDF data @@ -62,7 +62,6 @@ module globalData ! ---------- conversion factors ------------------------------------------------------------------- - real(dp) , public :: convTime2Days ! conversion factor to convert time to units of days real(dp) , public :: time_conv ! time conversion factor -- used to convert to mm/s real(dp) , public :: length_conv ! length conversion factor -- used to convert to mm/s @@ -90,10 +89,10 @@ module globalData type(var_info) , public :: meta_SEG (nVarsSEG ) ! stream segment properties type(var_info) , public :: meta_NTOPO (nVarsNTOPO ) ! network topology type(var_info) , public :: meta_PFAF (nVarsPFAF ) ! pfafstetter code - type(var_info_new) , public :: meta_rflx (nVarsRFLX ) ! reach flux variables - type(var_info_new) , public :: meta_irf_bas(nVarsIRFbas ) ! basin IRF routing fluxes/states - type(var_info_new) , public :: meta_kwt (nVarsKWT ) ! KWT routing fluxes/states - type(var_info_new) , public :: meta_irf (nVarsIRF ) ! IRF routing fluxes/states + type(meta_var) , public :: meta_rflx (nVarsRFLX ) ! reach flux variables + type(meta_var) , public :: meta_irf_bas(nVarsIRFbas ) ! basin IRF routing fluxes/states + type(meta_var) , public :: meta_kwt (nVarsKWT ) ! KWT routing fluxes/states + type(meta_var) , public :: meta_irf (nVarsIRF ) ! IRF routing fluxes/states ! ---------- data structures ---------------------------------------------------------------------- integer(i4b) , public :: nEns=1 ! number of ensemble @@ -102,20 +101,18 @@ module globalData integer(i4b) , public :: nRch ! number of reaches in the whole river network ! basin and reach IDs (to be removed) - integer(i4b) , allocatable , public :: basinID(:) ! HRU id + integer(i8b) , allocatable , public :: basinID(:) ! HRU id integer(i4b) , allocatable , public :: reachID(:) ! reach id ! DataTime data/variables integer(i4b) , public :: iTime ! time index at simulation time step - real(dp) , public :: startJulday ! julian day: start of routing simulation - real(dp) , public :: endJulday ! julian day: end of routing simulation - real(dp) , public :: refJulday ! julian day: reference - real(dp) , public :: modJulday ! julian day: simulation time step - real(dp) , public :: restartJulday ! julian day: restart drop off - real(dp) , allocatable , public :: roJulday(:) ! julian day: runoff input time real(dp) , allocatable , public :: timeVar(:) ! time variables (unit given by time variable) real(dp) , public :: TSEC(0:1) ! begning and end of time step (sec) - type(time) , public :: modTime(0:1) ! previous and current model time (yyyy:mm:dd:hh:mm:ss) + type(datetime) , public :: modTime(0:1) ! previous and current model time (yyyy:mm:dd:hh:mm:ss) + type(datetime) , public :: endCal ! simulation end date/time (yyyy:mm:dd:hh:mm:ss) + type(datetime) , public :: restCal ! desired restart date/time (yyyy:mm:dd:hh:mm:ss) + type(datetime) , public :: dropCal ! restart dropoff date/time (yyyy:mm:dd:hh:mm:ss) + logical(lgt) , public :: restartAlarm ! alarm to triger restart file writing ! simulation output netcdf type(nc) , public :: simout_nc diff --git a/route/build/src/irf_route.f90 b/route/build/src/irf_route.f90 index 1bd6c4f4..1002dcb2 100644 --- a/route/build/src/irf_route.f90 +++ b/route/build/src/irf_route.f90 @@ -5,13 +5,14 @@ module irf_route_module ! data type USE dataTypes, only : STRFLX ! fluxes in each reach USE dataTypes, only : RCHTOPO ! Network topology +USE dataTypes, only : RCHPRP ! Reach parameter ! global parameters USE public_var, only : realMissing ! missing value for real number USE public_var, only : integerMissing ! missing value for integer number USE globalData, only : nThreads ! number of threads used for openMP -! privary implicit none + private public::irf_route @@ -25,12 +26,13 @@ subroutine irf_route(iEns, & ! input: index of runoff ensemble to be p river_basin, & ! input: river basin information (mainstem, tributary outlet etc.) ixDesire, & ! input: reachID to be checked by on-screen pringing NETOPO_in, & ! input: reach topology data structure + RPARAM_in, & ! input: reach parameter data structure RCHFLX_out, & ! inout: reach flux data structure ierr, message, & ! output: error control ixSubRch) ! optional input: subset of reach indices to be processed - ! global routing data - USE dataTypes, only : subbasin_omp ! mainstem+tributary data structures + USE dataTypes, ONLY : subbasin_omp ! mainstem+tributary data structures + USE model_finalize, ONLY : handle_err implicit none ! Input @@ -38,6 +40,7 @@ subroutine irf_route(iEns, & ! input: index of runoff ensemble to be p type(subbasin_omp), intent(in), allocatable :: river_basin(:) ! river basin information (mainstem, tributary outlet etc.) integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output ! Output type(RCHTOPO), intent(in), allocatable :: NETOPO_in(:) ! River Network topology + type(RCHPRP), intent(in), allocatable :: RPARAM_in(:) ! River reach parameter ! inout TYPE(STRFLX), intent(inout), allocatable :: RCHFLX_out(:,:) ! Reach fluxes (ensembles, space [reaches]) for decomposed domains ! output variables @@ -55,18 +58,13 @@ subroutine irf_route(iEns, & ! input: index of runoff ensemble to be p integer(i4b) :: iTrib ! loop indices - branch integer(i4b) :: ix ! loop indices stream order ! variables needed for timing - integer*8 :: cr ! rate - integer*8 :: startTime,endTime ! date/time for the start and end of the initialization - real(dp) :: elapsedTime ! elapsed time for the process !integer(i4b) :: omp_get_thread_num !integer(i4b), allocatable :: ixThread(:) ! thread id !integer*8, allocatable :: openMPend(:) ! time for the start of the parallelization section !integer*8, allocatable :: timeTribStart(:) ! time Tributaries start !real(dp), allocatable :: timeTrib(:) ! time spent on each Tributary - ! initialize error control ierr=0; message='irf_route/' - call system_clock(count_rate=cr) ! number of reach check if (size(NETOPO_in)/=size(RCHFLX_out(iens,:))) then @@ -89,7 +87,6 @@ subroutine irf_route(iEns, & ! input: index of runoff ensemble to be p nOrder = size(river_basin) - call system_clock(startTime) do ix = 1,nOrder @@ -107,6 +104,7 @@ subroutine irf_route(iEns, & ! input: index of runoff ensemble to be p !$OMP shared(river_basin) & ! data structure shared !$OMP shared(doRoute) & ! data array shared !$OMP shared(NETOPO_in) & ! data structure shared +!$OMP shared(RPARAM_in) & ! data structure shared !$OMP shared(RCHFLX_out) & ! data structure shared !$OMP shared(ix, iEns, ixDesire) & ! indices shared !$OMP firstprivate(nTrib) @@ -120,8 +118,8 @@ subroutine irf_route(iEns, & ! input: index of runoff ensemble to be p seg:do iSeg=1,river_basin(ix)%branch(iTrib)%nRch jSeg = river_basin(ix)%branch(iTrib)%segIndex(iSeg) if (.not. doRoute(jSeg)) cycle - call segment_irf(iEns, jSeg, ixDesire, NETOPO_IN, RCHFLX_out, ierr, cmessage) -! if(ierr/=0)then; ixmessage(iTrib)=trim(message)//trim(cmessage); exit; endif + call segment_irf(iEns, jSeg, ixDesire, NETOPO_IN, RPARAM_in, RCHFLX_out, ierr, cmessage) + if(ierr/=0) call handle_err(ierr, trim(message)//trim(cmessage)) end do seg ! call system_clock(openMPend(iTrib)) ! timeTrib(iTrib) = real(openMPend(iTrib)-timeTribStart(iTrib), kind(dp)) @@ -137,10 +135,6 @@ subroutine irf_route(iEns, & ! input: index of runoff ensemble to be p end do ! basin loop - call system_clock(endTime) - elapsedTime = real(endTime-startTime, kind(dp))/real(cr) - !write(*,"(A,1PG15.7,A)") ' elapsed-time [routing/irf] = ', elapsedTime, ' s' - end subroutine irf_route @@ -153,6 +147,7 @@ subroutine segment_irf(& segIndex, & ! input: index of runoff ensemble to be processed ixDesire, & ! input: reachID to be checked by on-screen pringing NETOPO_in, & ! input: reach topology data structure + RPARAM_in, & ! input: reach parameter data structure ! inout RCHFLX_out, & ! inout: reach flux data structure ! output @@ -160,27 +155,29 @@ subroutine segment_irf(& implicit none ! Input - INTEGER(I4B), intent(IN) :: iEns ! runoff ensemble to be routed - INTEGER(I4B), intent(IN) :: segIndex ! segment where routing is performed - INTEGER(I4B), intent(IN) :: ixDesire ! index of the reach for verbose output + integer(i4b), intent(in) :: iEns ! runoff ensemble to be routed + integer(i4b), intent(in) :: segIndex ! segment where routing is performed + integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output type(RCHTOPO),intent(in), allocatable :: NETOPO_in(:) ! River Network topology + type(RCHPRP), intent(in), allocatable :: RPARAM_in(:) ! River reach parameter ! inout - TYPE(STRFLX), intent(inout), allocatable :: RCHFLX_out(:,:) ! Reach fluxes (ensembles, space [reaches]) for decomposed domains + type(STRFLX), intent(inout), allocatable :: RCHFLX_out(:,:) ! Reach fluxes (ensembles, space [reaches]) for decomposed domains ! Output integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! Local variables to - type(STRFLX), allocatable :: uprflux(:) ! upstream Reach fluxes - INTEGER(I4B) :: nUps ! number of upstream segment - INTEGER(I4B) :: iUps ! upstream reach index - INTEGER(I4B) :: iRch_ups ! index of upstream reach in NETOPO - INTEGER(I4B) :: ntdh ! number of time steps in IRF + real(dp) :: q_upstream ! total discharge at top of the reach being processed + integer(i4b) :: nUps ! number of upstream segment + integer(i4b) :: iUps ! upstream reach index + integer(i4b) :: iRch_ups ! index of upstream reach in NETOPO + integer(i4b) :: ntdh ! number of time steps in IRF + integer(i4b) :: itdh ! loop index for unit hydrograph + character(len=strLen) :: fmt1 ! format string character(len=strLen) :: cmessage ! error message from subroutine - ! initialize error control ierr=0; message='segment_irf/' - ! route streamflow through the river network + ! initialize future discharge array at first time if (.not.allocated(RCHFLX_out(iens,segIndex)%QFUTURE_IRF))then ntdh = size(NETOPO_in(segIndex)%UH) @@ -192,23 +189,22 @@ subroutine segment_irf(& end if - ! identify number of upstream segments of the reach being processed + ! get discharge coming from upstream nUps = size(NETOPO_in(segIndex)%UREACHI) - - allocate(uprflux(nUps), stat=ierr, errmsg=cmessage) - if(ierr/=0)then; message=trim(message)//trim(cmessage)//': uprflux'; return; endif - + q_upstream = 0.0_dp if (nUps>0) then do iUps = 1,nUps iRch_ups = NETOPO_in(segIndex)%UREACHI(iUps) ! index of upstream of segIndex-th reach - uprflux(iUps) = RCHFLX_out(iens,iRch_ups) + q_upstream = q_upstream + RCHFLX_out(iens, iRch_ups)%REACH_Q_IRF end do endif - ! perform river network UH routing + ! perform UH convolution call conv_upsbas_qr(NETOPO_in(segIndex)%UH, & ! input: reach unit hydrograph - uprflux, & ! input: upstream reach fluxes + q_upstream, & ! input: total discharge at top of the reach being processed RCHFLX_out(iens,segIndex), & ! inout: updated fluxes at reach + RPARAM_in(segIndex)%QTAKE, & ! input: abstraction(-)/injection(+) [m3/s] + RPARAM_in(segIndex)%MINFLOW, & ! input: minimum environmental flow [m3/s] ierr, message) ! output: error control if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif @@ -216,9 +212,16 @@ subroutine segment_irf(& RCHFLX_out(iEns,segIndex)%CHECK_IRF=.True. ! check - if(NETOPO_in(segIndex)%REACHIX == ixDesire)then - print*, 'RCHFLX_out(iens,segIndex)%BASIN_QR(1),RCHFLX_out(iens,segIndex)%REACH_Q_IRF = ', & - RCHFLX_out(iens,segIndex)%BASIN_QR(1),RCHFLX_out(iens,segIndex)%REACH_Q_IRF + if(segIndex==ixDesire)then + ntdh = size(NETOPO_in(segIndex)%UH) + write(fmt1,'(A,I5,A)') '(A, 1X',ntdh,'(1X,F20.7))' + write(*,'(2a)') new_line('a'),'** Check Impulse Response Function routing **' + write(*,'(a,x,I10,x,I10)') ' Reach index & ID =', segIndex, NETOPO_in(segIndex)%REACHID + write(*,fmt1) ' Unit-Hydrograph =', (NETOPO_in(segIndex)%UH(itdh), itdh=1,ntdh) + write(*,'(a)') ' * total discharge from upstream(q_upstream) [m3/s], local area discharge [m3/s], and Final discharge [m3/s]:' + write(*,'(a,x,F15.7)') ' q_upstream =', q_upstream + write(*,'(a,x,F15.7)') ' RCHFLX_out%BASIN_QR(1) =', RCHFLX_out(iens,segIndex)%BASIN_QR(1) + write(*,'(a,x,F15.7)') ' RCHFLX_out%REACH_Q_IRF =', RCHFLX_out(iens,segIndex)%REACH_Q_IRF endif end subroutine segment_irf @@ -228,63 +231,73 @@ end subroutine segment_irf ! subroutine: Compute delayed runoff from the upstream segments ! ********************************************************************* subroutine conv_upsbas_qr(reach_uh, & ! input: reach unit hydrograph - rflux_ups, & ! input: upstream reach fluxes + q_upstream, & ! input: rflux, & ! input: input flux at reach + Qtake, & ! input: abstraction(-)/injection(+) [m3/s] + Qmin, & ! input: minimum environmental flow [m3/s] ierr, message) ! output: error control ! ---------------------------------------------------------------------------------------- ! Details: Convolute runoff volume of upstream at one reach at one time step ! ---------------------------------------------------------------------------------------- + USE public_var, ONLY: dt + implicit none ! Input real(dp), intent(in) :: reach_uh(:) ! reach unit hydrograph - type(STRFLX), intent(in) :: rflux_ups(:) ! upstream Reach fluxes + real(dp), intent(in) :: q_upstream ! total discharge at top of the reach being processed + real(dp), intent(in) :: Qtake ! abstraction(-)/injection(+) [m3/s] + real(dp), intent(in) :: Qmin ! minimum environmental flow [m3/s] ! inout type(STRFLX), intent(inout) :: rflux ! current Reach fluxes ! Output integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! Local variables to - real(dp) :: q_upstream ! total discharge at top of the reach being processed - INTEGER(I4B) :: ntdh ! number of UH data - INTEGER(I4B) :: itdh ! index of UH data (i.e.,future time step) - INTEGER(I4B) :: nUps ! number of all upstream segment - INTEGER(I4B) :: iUps ! loop indices for u/s reaches + real(dp) :: QupMod ! modified total discharge at top of the reach being processed + real(dp) :: Qabs ! maximum allowable water abstraction rate [m3/s] + real(dp) :: Qmod ! abstraction rate to be taken from outlet discharge [m3/s] + integer(i4b) :: ntdh ! number of UH data (i.e., number of future time step + integer(i4b) :: itdh ! index of UH data - ! initialize error control ierr=0; message='conv_upsbas_qr/' - ! identify number of upstream segments of the reach being processed - nUps = size(rflux_ups) - - ! Find out total q at top of a segment - q_upstream = 0.0_dp - if(nUps>0)then - do iUps = 1,nUps - q_upstream = q_upstream + rflux_ups(iUps)%REACH_Q_IRF - end do - endif + ! Q injection, add at top of reach + QupMod = q_upstream + if (Qtake>0) then + QupMod = QupMod+ Qtake + end if ! place a fraction of runoff in future time steps - ntdh = size(reach_uh) ! identify the number of future time steps of UH for a given segment + ntdh = size(reach_uh) ! number of future time steps of UH for a given segment do itdh=1,ntdh - rflux%QFUTURE_IRF(itdh) = rflux%QFUTURE_IRF(itdh) & - + reach_uh(itdh)*q_upstream + rflux%QFUTURE_IRF(itdh) = rflux%QFUTURE_IRF(itdh)+ reach_uh(itdh)*QupMod enddo - ! Add local routed flow + ! compute volume in reach + rflux%REACH_VOL(0) = rflux%REACH_VOL(1) + rflux%REACH_VOL(1) = rflux%REACH_VOL(0) + (QupMod - rflux%QFUTURE_IRF(1))/dt + + ! Add local routed flow at the bottom of reach rflux%REACH_Q_IRF = rflux%QFUTURE_IRF(1) + rflux%BASIN_QR(1) - ! move array back use eoshift - !rflux%QFUTURE_IRF=eoshift(rflux%QFUTURE_IRF,shift=1) + ! Q abstraction + ! Compute actual abstraction (Qabs) m3/s - values should be negative + ! Compute abstraction (Qmod) m3 taken from outlet discharge (REACH_Q_IRF) + ! Compute REACH_Q_IRF subtracted from Qmod abstraction + ! Compute REACH_VOL subtracted from total abstraction minus abstraction from outlet discharge + if (Qtake<0) then + Qabs = max(-(rflux%REACH_VOL(1)/dt+rflux%REACH_Q_IRF), Qtake) + Qmod = min(rflux%REACH_VOL(1) + Qabs*dt, 0._dp) + rflux%REACH_Q_IRF = max(rflux%REACH_Q_IRF + Qmod/dt, Qmin) + rflux%REACH_VOL(1) = rflux%REACH_VOL(1) + (Qabs*dt - Qmod) + end if - do itdh=2,ntdh - rflux%QFUTURE_IRF(itdh-1) = rflux%QFUTURE_IRF(itdh) - enddo + ! move array back use eoshift + rflux%QFUTURE_IRF=eoshift(rflux%QFUTURE_IRF,shift=1) rflux%QFUTURE_IRF(ntdh) = 0._dp end subroutine conv_upsbas_qr end module irf_route_module - diff --git a/route/build/src/kwt_route.f90 b/route/build/src/kwt_route.f90 index d55e645b..97829e99 100644 --- a/route/build/src/kwt_route.f90 +++ b/route/build/src/kwt_route.f90 @@ -1,28 +1,28 @@ -module kwt_route_module +MODULE kwt_route_module !numeric type -use nrtype +USE nrtype ! data types -USE dataTypes, only : FPOINT ! particle -USE dataTypes, only : KREACH ! collection of particles in a given reach -USE dataTypes, only : STRFLX ! fluxes in each reach -USE dataTypes, only : RCHTOPO ! Network topology -USE dataTypes, only : RCHPRP ! Reach parameter +USE dataTypes, ONLY : FPOINT ! particle +USE dataTypes, ONLY : KREACH ! collection of particles in a given reach +USE dataTypes, ONLY : STRFLX ! fluxes in each reach +USE dataTypes, ONLY : RCHTOPO ! Network topology +USE dataTypes, ONLY : RCHPRP ! Reach parameter ! global data -USE public_var, only : verySmall ! a very small value -USE public_var, only : realMissing ! missing value for real number -USE public_var, only : integerMissing ! missing value for integer number +USE public_var, ONLY : runoffMin ! minimum runoff +USE public_var, ONLY : verySmall ! a very small value +USE public_var, ONLY : realMissing ! missing value for real number +USE public_var, ONLY : integerMissing ! missing value for integer number ! utilities -use nr_utility_module, only : arth ! Num. Recipies utilities -USE time_utils_module, only : elapsedSec ! calculate the elapsed time +USE nr_utility_module, ONLY : arth ! Num. Recipies utilities -! privary implicit none + private public::kwt_route -contains +CONTAINS ! ********************************************************************* ! subroutine: route kinematic waves through the river network @@ -38,7 +38,9 @@ SUBROUTINE kwt_route(iens, & ! input: ensemble index ierr,message, & ! output: error control ixSubRch) ! optional input: subset of reach indices to be processed - USE dataTypes, only : subbasin_omp ! mainstem+tributary data strucuture + USE dataTypes, ONLY : subbasin_omp ! mainstem+tributary data strucuture + USE model_finalize, ONLY : handle_err + implicit none ! Input integer(i4b), intent(in) :: iEns ! ensemble member @@ -66,18 +68,13 @@ SUBROUTINE kwt_route(iens, & ! input: ensemble index integer(i4b) :: iTrib ! loop indices - branch integer(i4b) :: ix ! loop indices stream order ! variables needed for timing - integer*8 :: cr ! rate - integer*8 :: startTime,endTime ! start and end time stamps - real(dp) :: elapsedTime ! elapsed time for the process ! integer(i4b) :: omp_get_thread_num ! integer(i4b), allocatable :: ixThread(:) ! thread id ! integer*8, allocatable :: openMPend(:) ! time for the start of the parallelization section ! integer*8, allocatable :: timeTribStart(:) ! time Tributaries start ! real(dp), allocatable :: timeTrib(:) ! time spent on each Tributary - ! initialize error control ierr=0; message='kwt_route/' - call system_clock(count_rate=cr) ! number of reach check if (size(NETOPO_in)/=size(RCHFLX_out(iens,:))) then @@ -98,8 +95,6 @@ SUBROUTINE kwt_route(iens, & ! input: ensemble index nOrder = size(river_basin) - call system_clock(startTime) - do ix = 1, nOrder nTrib=size(river_basin(ix)%branch) @@ -134,7 +129,7 @@ SUBROUTINE kwt_route(iens, & ! input: ensemble index jSeg = river_basin(ix)%branch(iTrib)%segIndex(iSeg) if (.not. doRoute(jSeg)) cycle ! route kinematic waves through the river network - call QROUTE_RCH(iEns,jSeg, & ! input: array indices + call qroute_rch(iEns,jSeg, & ! input: array indices ixDesire, & ! input: index of the desired reach T0,T1, & ! input: start and end of the time step LAKEFLAG, & ! input: flag if lakes are to be processed @@ -143,7 +138,7 @@ SUBROUTINE kwt_route(iens, & ! input: ensemble index KROUTE_out, & ! inout: reach state data structure RCHFLX_out, & ! inout: reach flux data structure ierr,cmessage) ! output: error control - !if (ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + if(ierr/=0) call handle_err(ierr, trim(message)//trim(cmessage)) end do seg ! call system_clock(openMPend(iTrib)) ! timeTrib(iTrib) = real(openMPend(iTrib)-timeTribStart(iTrib), kind(dp)) @@ -159,17 +154,13 @@ SUBROUTINE kwt_route(iens, & ! input: ensemble index end do ! basin loop - call system_clock(endTime) - elapsedTime = real(endTime-startTime, kind(dp))/real(cr) -! write(*,"(A,1PG15.7,A)") ' elapsed-time [routing/kwt] = ', elapsedTime, ' s' - END SUBROUTINE kwt_route ! ********************************************************************* ! subroutine: route kinematic waves at one segment ! ********************************************************************* - subroutine QROUTE_RCH(IENS,JRCH, & ! input: array indices + SUBROUTINE qroute_rch(IENS,JRCH, & ! input: array indices ixDesire, & ! input: index of the reach for verbose output T0,T1, & ! input: start and end of the time step LAKEFLAG, & ! input: flag if lakes are to be processed @@ -179,8 +170,7 @@ subroutine QROUTE_RCH(IENS,JRCH, & ! input: array indices RCHFLX_out, & ! inout: reach flux data structure ierr,message, & ! output: error control RSTEP) ! optional input: retrospective time step offset - ! public data - USE public_var, only : MAXQPAR ! maximum number of waves per reach + USE public_var, ONLY : MAXQPAR ! maximum number of waves per reach ! ---------------------------------------------------------------------------------------- ! Creator(s): ! Ross Woods, 1997 (original code) @@ -223,99 +213,94 @@ subroutine QROUTE_RCH(IENS,JRCH, & ! input: array indices ! ! * added additional comments ! - ! * all variables are defined (IMPLICIT NONE) and described (comments) + ! * all variables are defined (implicit none) and described (comments) ! ! * use of a new data structure (KROUTE_out) to hold and update the flow particles ! ! * upgrade to F90 (especially structured variables and dynamic memory allocation) ! ! ---------------------------------------------------------------------------------------- - ! Future revisions: - ! - ! (none planned) - ! - ! ---------------------------------------------------------------------------------------- - implicit none - ! Input - integer(i4b), intent(in) :: IENS ! ensemble member - integer(i4b), intent(in) :: JRCH ! reach to process - integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output - real(dp), intent(in) :: T0,T1 ! start and end of the time step (seconds) - integer(i4b), intent(in) :: LAKEFLAG ! >0 if processing lakes - type(RCHTOPO),intent(in), allocatable :: NETOPO_in(:) ! River Network topology - type(RCHPRP), intent(in), allocatable :: RPARAM_in(:) ! River reach parameter - integer(i4b), intent(in), optional :: RSTEP ! retrospective time step offset - ! inout - type(KREACH), intent(inout), allocatable :: KROUTE_out(:,:) ! reach state data - TYPE(STRFLX), intent(inout), allocatable :: RCHFLX_out(:,:) ! Reach fluxes (ensembles, space [reaches]) for decomposed domains - ! output variables - integer(i4b), intent(out) :: ierr ! error code - character(*), intent(out) :: message ! error message - ! (1) extract flow from upstream reaches and append to the non-routed flow in JRCH - INTEGER(I4B) :: NUPS ! number of upstream reaches - REAL(DP),DIMENSION(:),allocatable :: Q_JRCH ! flow in downstream reach JRCH - REAL(DP),DIMENSION(:),allocatable :: TENTRY ! entry time to JRCH (exit time u/s) - INTEGER(I4B) :: NQ1 ! # flow particles - ! (2) route flow within the current [JRCH] river segment - INTEGER(I4B) :: ROFFSET ! retrospective offset due to rstep - REAL(DP) :: T_START ! start of time step - REAL(DP) :: T_END ! end of time step - REAL(DP),DIMENSION(:),allocatable :: T_EXIT ! time particle expected exit JRCH - LOGICAL(LGT),DIMENSION(:),allocatable :: FROUTE ! routing flag .T. if particle exits - INTEGER(I4B) :: NQ2 ! # flow particles (<=NQ1 b/c merge) - ! (3) calculate time-step averages - INTEGER(I4B) :: NR ! # routed particles - INTEGER(I4B) :: NN ! # non-routed particles - REAL(DP),DIMENSION(2) :: TNEW ! start/end of time step - REAL(DP),DIMENSION(1) :: QNEW ! interpolated flow - ! (4) housekeeping - REAL(DP) :: Q_END ! flow at the end of the timestep - REAL(DP) :: TIMEI ! entry time at the end of the timestep - TYPE(FPOINT),allocatable,DIMENSION(:) :: NEW_WAVE ! temporary wave - LOGICAL(LGT) :: INIT=.TRUE. ! used to initialize pointers - ! random stuff - integer(i4b) :: IWV ! rech index - character(len=strLen) :: fmt1,fmt2 ! format string - character(len=strLen) :: CMESSAGE ! error message for downwind routine - - ! initialize error control - ierr=0; message='QROUTE_RCH/' - ! ---------------------------------------------------------------------------------------- - ! (0) INITIALIZE POINTERS - ! ---------------------------------------------------------------------------------------- - if(INIT) then - INIT=.false. - endif + implicit none + ! Input + integer(i4b), intent(in) :: IENS ! ensemble member + integer(i4b), intent(in) :: JRCH ! reach to process + integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output + real(dp), intent(in) :: T0,T1 ! start and end of the time step (seconds) + integer(i4b), intent(in) :: LAKEFLAG ! >0 if processing lakes + type(RCHTOPO),intent(in), allocatable :: NETOPO_in(:) ! River Network topology + type(RCHPRP), intent(in), allocatable :: RPARAM_in(:) ! River reach parameter + integer(i4b), intent(in), optional :: RSTEP ! retrospective time step offset + ! inout + type(KREACH), intent(inout), allocatable :: KROUTE_out(:,:) ! reach state data + type(STRFLX), intent(inout), allocatable :: RCHFLX_out(:,:) ! Reach fluxes (ensembles, space [reaches]) for decomposed domains + ! output variables + integer(i4b), intent(out) :: ierr ! error code + character(*), intent(out) :: message ! error message + ! (1) extract flow from upstream reaches and append to the non-routed flow in JRCH + integer(i4b) :: NUPS ! number of upstream reaches + real(dp),dimension(:),allocatable :: Q_JRCH ! flow in downstream reach JRCH + real(dp),dimension(:),allocatable :: TENTRY ! entry time to JRCH (exit time u/s) + integer(i4b) :: NQ1 ! # flow particles + ! (2) route flow within the current [JRCH] river segment + integer(I4B) :: ROFFSET ! retrospective offset due to rstep + real(dp) :: T_START ! start of time step + real(dp) :: T_END ! end of time step + real(dp),dimension(:),allocatable :: T_EXIT ! time particle expected exit JRCH + logical(LGT),dimension(:),allocatable :: FROUTE ! routing flag .T. if particle exits + integer(I4B) :: NQ2 ! # flow particles (<=NQ1 b/c merge) + ! (3) calculate time-step averages + integer(I4B) :: NR ! # routed particles + integer(I4B) :: NN ! # non-routed particles + real(dp),dimension(2) :: TNEW ! start/end of time step + real(dp),dimension(1) :: QNEW ! interpolated flow + ! (4) housekeeping + real(dp) :: Q_END ! flow at the end of the timestep + real(dp) :: TIMEI ! entry time at the end of the timestep + TYPE(FPOINT),allocatable,dimension(:) :: NEW_WAVE ! temporary wave + ! random stuff + integer(i4b) :: IWV ! rech index + character(len=strLen) :: fmt1,fmt2 ! format string + character(len=strLen) :: CMESSAGE ! error message for downwind routine + + ierr=0; message='qroute_rch/' - if(JRCH==ixDesire) write(*,"('JRCH=',I10)") JRCH - if(JRCH==ixDesire) write(*,"('T0-T1=',F20.7,1x,F20.7)") T0, T1 + if(JRCH==ixDesire) then + write(*,'(2a)') new_line('a'),'** Check kinematic wave tracking routing **' + write(*,"(a,x,I10,x,I10)") ' Reach index & ID =', JRCH, NETOPO_in(JRCH)%REACHID + write(*,"(a,x,F20.7,1x,F20.7)") ' time step(T0,T1) =', T0, T1 + write(*,'(a,x,F15.7)') ' RPARAM_in%R_SLOPE =', RPARAM_in(JRCH)%R_SLOPE + write(*,'(a,x,F15.7)') ' RPARAM_in%R_MAN_N =', RPARAM_in(JRCH)%R_MAN_N + write(*,'(a,x,F15.7)') ' RPARAM_in%R_WIDTH =', RPARAM_in(JRCH)%R_WIDTH + end if + + RCHFLX_out(IENS,JRCH)%TAKE=0.0_dp ! initialize take from this reach - RCHFLX_out(IENS,JRCH)%TAKE=0.0_dp ! initialize take from this reach ! ---------------------------------------------------------------------------------------- ! (1) EXTRACT FLOW FROM UPSTREAM REACHES & APPEND TO THE NON-ROUTED FLOW PARTICLES IN JRCH ! ---------------------------------------------------------------------------------------- NUPS = count(NETOPO_in(JRCH)%goodBas) ! number of desired upstream reaches !NUPS = size(NETOPO_in(JRCH)%UREACHI) ! number of upstream reaches - IF (NUPS.GT.0) THEN - call GETUSQ_RCH(IENS,JRCH,LAKEFLAG,T0,T1,ixDesire, & ! input + if (NUPS.GT.0) then + call getusq_rch(IENS,JRCH,LAKEFLAG,T0,T1,ixDesire, & ! input NETOPO_in,RPARAM_in,RCHFLX_out, & ! input KROUTE_out, & ! inout Q_JRCH,TENTRY,T_EXIT,ierr,cmessage,& ! output RSTEP) ! optional input if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - ! check for negative flow - if (MINVAL(Q_JRCH).lt.0.0_dp) then + + if (minval(Q_JRCH).lt.0.0_dp) then ierr=20; message=trim(message)//'negative flow extracted from upstream reach'; return endif - ! check + if(JRCH==ixDesire)then - write(fmt1,'(A,I5,A)') '(A,1X',size(Q_JRCH),'(1X,F20.7))' - write(*,fmt1) 'Initial_Q_JRCH=', (Q_JRCH(IWV), IWV=0,size(Q_JRCH)-1) + write(fmt1,'(A,I5,A)') '(A, 1X',size(Q_JRCH),'(1X,F20.7))' + write(*,'(a)') ' * Wave discharge from upstream reaches (Q_JRCH) [m2/s]:' + write(*,fmt1) ' Q_JRCH=',(Q_JRCH(IWV), IWV=0,size(Q_JRCH)-1) endif else ! set flow in headwater reaches to modelled streamflow from time delay histogram RCHFLX_out(IENS,JRCH)%REACH_Q = RCHFLX_out(IENS,JRCH)%BASIN_QR(1) - if (allocated(KROUTE_out(IENS,JRCH)%KWAVE)) THEN + if (allocated(KROUTE_out(IENS,JRCH)%KWAVE)) then deallocate(KROUTE_out(IENS,JRCH)%KWAVE,STAT=IERR) if(ierr/=0)then; message=trim(message)//'problem deallocating space for KROUTE_out'; return; endif endif @@ -326,20 +311,25 @@ subroutine QROUTE_RCH(IENS,JRCH, & ! input: array indices KROUTE_out(IENS,JRCH)%KWAVE(0)%TR=-9999 KROUTE_out(IENS,JRCH)%KWAVE(0)%RF=.False. KROUTE_out(IENS,JRCH)%KWAVE(0)%QM=-9999 - ! check - if(JRCH==ixDesire) print*, 'JRCH, RCHFLX_out(IENS,JRCH)%REACH_Q = ', JRCH, RCHFLX_out(IENS,JRCH)%REACH_Q + + if(JRCH==ixDesire) then + write(*,'(a)') ' * Final discharge (RCHFLX_out(IENS,JRCH)%REACH_Q) [m3/s]:' + write(*,'(x,F20.7)') RCHFLX_out(IENS,JRCH)%REACH_Q + end if return ! no upstream reaches (routing for sub-basins done using time-delay histogram) endif + ! ---------------------------------------------------------------------------------------- ! (2) REMOVE FLOW PARTICLES (REDUCE MEMORY USAGE AND PROCESSING TIME) ! ---------------------------------------------------------------------------------------- if (size(Q_JRCH).GT.MAXQPAR) then - call REMOVE_RCH(MAXQPAR,Q_JRCH,TENTRY,T_EXIT,ierr,cmessage) + call remove_rch(MAXQPAR,Q_JRCH,TENTRY,T_EXIT,ierr,cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif endif - NQ1 = SIZE(Q_JRCH)-1 ! -1 because of the zero element + NQ1 = size(Q_JRCH)-1 ! -1 because of the zero element + ! ---------------------------------------------------------------------------------------- - ! (3) ROUTE FLOW WITHIN THE CURRENT [JRCH] RIVER SEGMENT + ! (x) Water use - take out (Qtake is negative) ! ---------------------------------------------------------------------------------------- ! set the retrospective offset if (.not.present(RSTEP)) then @@ -350,35 +340,63 @@ subroutine QROUTE_RCH(IENS,JRCH, & ! input: array indices ! set time boundaries T_START = T0 - (T1 - T0)*ROFFSET T_END = T1 - (T1 - T0)*ROFFSET + + if (RPARAM_in(jrch)%QTAKE < 0._dp) then + call extract_from_rch(iens, jrch, & ! input: ensemble and reach indices + T_START, T_END, & ! input: time [sec] of current time step bounds + RPARAM_in, & ! input: river reach parameters + RPARAM_in(jrch)%QTAKE, & ! input: target Qtake (minus) + ixDesire, & ! input: + Q_JRCH, T_EXIT, TENTRY, & ! inout: discharge and exit time for particle + ierr,cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + end if + + ! ---------------------------------------------------------------------------------------- + ! (3) ROUTE FLOW WITHIN THE CURRENT [JRCH] RIVER SEGMENT + ! ---------------------------------------------------------------------------------------- allocate(FROUTE(0:NQ1),STAT=IERR) if(ierr/=0)then; message=trim(message)//'problem allocating space for FROUTE'; return; endif - FROUTE(0) = .TRUE.; FROUTE(1:NQ1)=.FALSE. ! init. routing flags + FROUTE(0) = .true.; FROUTE(1:NQ1)=.false. ! init. routing flags + ! route flow through the current [JRCH] river segment (Q_JRCH in units of m2/s) - call KINWAV_RCH(JRCH,T_START,T_END,ixDesire, & ! input: location and time + call kinwav_rch(JRCH,T_START,T_END,ixDesire, & ! input: location and time NETOPO_in, RPARAM_in, & ! input: river data structure Q_JRCH(1:NQ1),TENTRY(1:NQ1),T_EXIT(1:NQ1),FROUTE(1:NQ1), & ! inout: kwt states NQ2,ierr,cmessage) ! output: if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + if(JRCH == ixDesire)then write(fmt1,'(A,I5,A)') '(A,1X',NQ1+1,'(1X,F20.7))' write(fmt2,'(A,I5,A)') '(A,1X',NQ1+1,'(1X,L))' - write(*,fmt1) 'Q_JRCH=',(Q_JRCH(IWV), IWV=0,NQ1) - write(*,fmt2) 'FROUTE=',(FROUTE(IWV), IWV=0,NQ1) - write(*,fmt1) 'TENTRY=',(TENTRY(IWV), IWV=0,NQ1) - write(*,fmt1) 'T_EXIT=',(T_EXIT(IWV), IWV=0,NQ1) + write(*,'(a)') ' * After routed: wave discharge (Q_JRCH) [m2/s], isExit(FROUTE), entry time (TENTRY) [s], and exit time (T_EXIT) [s]:' + write(*,fmt1) ' Q_JRCH=',(Q_JRCH(IWV), IWV=0,NQ1) + write(*,fmt1) ' TENTRY=',(TENTRY(IWV), IWV=0,NQ1) + write(*,fmt1) ' T_EXIT=',(T_EXIT(IWV), IWV=0,NQ1) + write(*,fmt2) ' FROUTE=',(FROUTE(IWV), IWV=0,NQ1) endif + ! ---------------------------------------------------------------------------------------- ! (4) COMPUTE TIME-STEP AVERAGES ! ---------------------------------------------------------------------------------------- - NR = COUNT(FROUTE)-1 ! -1 because of the zero element (last routed) + NR = count(FROUTE)-1 ! -1 because of the zero element (last routed) NN = NQ2-NR ! number of non-routed points TNEW = (/T_START,T_END/) + ! (zero position last routed; use of NR+1 instead of NR keeps next expected routed point) - call INTERP_RCH(T_EXIT(0:NR+1),Q_JRCH(0:NR+1),TNEW,QNEW,IERR,CMESSAGE) + call interp_rch(T_EXIT(0:NR+1),Q_JRCH(0:NR+1),TNEW,QNEW,IERR,CMESSAGE) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - if(JRCH == ixDesire) write(*,"('QNEW(1)=',1x,F10.7)") QNEW(1) + ! m2/s --> m3/s + instantaneous runoff from basin RCHFLX_out(IENS,JRCH)%REACH_Q = QNEW(1)*RPARAM_in(JRCH)%R_WIDTH + RCHFLX_out(IENS,JRCH)%BASIN_QR(1) + + if(JRCH == ixDesire)then + write(*,'(a)') ' * Time-ave. wave discharge that exit (QNEW(1)) [m2/s], local-area discharge (RCHFLX_out%BASIN_QR(1)) [m3/s] and Final discharge (RCHFLX_out%REACH_Q) [m3/s]:' + write(*,"(A,1x,F15.7)") ' QNEW(1) =', QNEW(1) + write(*,"(A,1x,F15.7)") ' RCHFLX_out%BASIN_QR(1) =', RCHFLX_out(IENS,JRCH)%BASIN_QR(1) + write(*,"(A,1x,F15.7)") ' RCHFLX_out%REACH_Q =', RCHFLX_out(IENS,JRCH)%REACH_Q + endif + ! ---------------------------------------------------------------------------------------- ! (5) HOUSEKEEPING ! ---------------------------------------------------------------------------------------- @@ -400,7 +418,7 @@ subroutine QROUTE_RCH(IENS,JRCH, & ! input: array indices endif ! insert the interpolated point (TI is irrelevant, as the point is "routed") KROUTE_out(IENS,JRCH)%KWAVE(NR+1)%QF=Q_END; KROUTE_out(IENS,JRCH)%KWAVE(NR+1)%TI=TIMEI - KROUTE_out(IENS,JRCH)%KWAVE(NR+1)%TR=T_END; KROUTE_out(IENS,JRCH)%KWAVE(NR+1)%RF=.TRUE. + KROUTE_out(IENS,JRCH)%KWAVE(NR+1)%TR=T_END; KROUTE_out(IENS,JRCH)%KWAVE(NR+1)%RF=.true. ! add the output from kinwave... - skip NR+1 ! (when JRCH becomes IR routed points will be stripped out & the structures updated again) KROUTE_out(IENS,JRCH)%KWAVE(0:NR)%QF=Q_JRCH(0:NR); KROUTE_out(IENS,JRCH)%KWAVE(NR+2:NQ2+1)%QF=Q_JRCH(NR+1:NQ2) @@ -408,25 +426,27 @@ subroutine QROUTE_RCH(IENS,JRCH, & ! input: array indices KROUTE_out(IENS,JRCH)%KWAVE(0:NR)%TR=T_EXIT(0:NR); KROUTE_out(IENS,JRCH)%KWAVE(NR+2:NQ2+1)%TR=T_EXIT(NR+1:NQ2) KROUTE_out(IENS,JRCH)%KWAVE(0:NR)%RF=FROUTE(0:NR); KROUTE_out(IENS,JRCH)%KWAVE(NR+2:NQ2+1)%RF=FROUTE(NR+1:NQ2) KROUTE_out(IENS,JRCH)%KWAVE(0:NQ2+1)%QM=-9999 + ! implement water use !IF (NUSER.GT.0.AND.UCFFLAG.GE.1) THEN !CALL EXTRACT_FROM_RCH(IENS,JRCH,NR,Q_JRCH,T_EXIT,T_END,TNEW) !ENDIF + ! free up space for the next reach deallocate(Q_JRCH,TENTRY,T_EXIT,FROUTE,STAT=IERR) ! FROUTE defined in this sub-routine if(ierr/=0)then; message=trim(message)//'problem deallocating space for [Q_JRCH, TENTRY, T_EXIT, FROUTE]'; return; endif + ! *** ! remove flow particles from the most downstream reach ! if the last reach or lake inlet (and lakes are enabled), remove routed elements from memory -! IF ((NETOPO_in(JRCH)%DREACHI.LT.0 .and. basinType==2).OR. & ! if the last reach, then there is no downstream reach IF ((NETOPO_in(JRCH)%DREACHK<=0 ).OR. & ! if the last reach (down reach ID:DREACHK is negative), then there is no downstream reach (LAKEFLAG.EQ.1.AND.NETOPO_in(JRCH)%LAKINLT)) THEN ! if lake inlet ! copy data to a temporary wave if (allocated(NEW_WAVE)) THEN - DEALLOCATE(NEW_WAVE,STAT=IERR) + deallocate(NEW_WAVE,STAT=IERR) if(ierr/=0)then; message=trim(message)//'problem deallocating space for NEW_WAVE'; return; endif endif - ALLOCATE(NEW_WAVE(0:NN),STAT=IERR) ! NN = number non-routed (the zero element is the last routed point) + allocate(NEW_WAVE(0:NN),STAT=IERR) ! NN = number non-routed (the zero element is the last routed point) if(ierr/=0)then; message=trim(message)//'problem allocating space for NEW_WAVE'; return; endif NEW_WAVE(0:NN) = KROUTE_out(IENS,JRCH)%KWAVE(NR+1:NQ2+1) ! +1 because of the interpolated point ! re-size wave structure @@ -438,16 +458,120 @@ subroutine QROUTE_RCH(IENS,JRCH, & ! input: array indices if(ierr/=0)then; message=trim(message)//'problem allocating space for KROUTE_out'; return; endif ! copy data back to the wave structure and deallocate space for the temporary wave KROUTE_out(IENS,JRCH)%KWAVE(0:NN) = NEW_WAVE(0:NN) - DEALLOCATE(NEW_WAVE,STAT=IERR) - if(ierr/=0)then; message=trim(message)//'problem deallocating space for NEW_WAVE'; return; endif endif ! (if JRCH is the last reach) - end subroutine QROUTE_RCH + END SUBROUTINE qroute_rch + + ! ********************************************************************* + ! subroutine: wave discharge mod to extract water from the JRCH reach + ! ********************************************************************* + SUBROUTINE extract_from_rch(iens, jrch, & ! input: ensemble and reach indices + T_START, T_END, & ! input: start and end time [sec] for this time step + RPARAM_in, & ! input: river reach parameters + Qtake, & ! input: target Qtake (minus) + ixDesire, & ! input: + Q_JRCH, T_EXIT, TENTRY, & ! inout: discharge and exit time for particle + ierr,message) + implicit none + ! Input + integer(i4b), intent(in) :: iens ! ensemble member + integer(i4b), intent(in) :: jRch ! index of reach to process + real(dp), intent(in) :: T_START ! start time [s] + real(dp), intent(in) :: T_END ! end time [s] + type(RCHPRP), allocatable, intent(in) :: RPARAM_in(:) ! River reach parameter + real(dp), intent(in) :: Qtake ! target Q abstraction [m3/s] + integer(i4b), intent(in) :: ixDesire ! index of the reach for verbose output + ! inout + real(dp), allocatable, intent(inout) :: Q_JRCH(:) ! discharge of particle [m2/s] -- discharge for unit channel width + real(dp), allocatable, intent(inout) :: T_EXIT(:) ! time flow is expected to exit JR + real(dp), allocatable, intent(inout) :: TENTRY(:) ! time flow entered JR + integer(i4b), intent(out) :: ierr ! error code + character(*), intent(out) :: message ! error message + ! Local variables + real(dp) :: totQ ! total available flow + real(dp) :: Qabs ! actual abstracted water + real(dp) :: Qfrac ! fraction of target abstraction to total available flow + real(dp) :: alfa ! constant = 5/3 + real(dp) :: K ! sqrt(slope)/mannings N + real(dp) :: TP(2) ! start/end of time step + real(dp) :: Qavg(1) ! time average flow [m2/s] + real(dp), allocatable :: wc(:) ! wave celerity [m/s] + real(dp), allocatable :: q_jrch_mod(:) ! wave flow remaining after abstract [m2/s] + real(dp), allocatable :: q_jrch_abs(:) ! wave abstracted flow [m2/s] + real(dp), allocatable :: t_exit_mod(:) ! wave expected exit time [sec] + integer(i4b) :: iw ! loop index for wave + integer(i4b) :: NR ! number of particle (not including zero index) + character(len=strLen) :: fmt1 ! format string + character(len=strLen) :: cmessage ! error message for downwind routine + + ierr=0; message='extract_from_rch/' + + ! uniform flow parameters + alfa = 5._dp/3._dp + K = sqrt(RPARAM_in(JRCH)%R_SLOPE)/RPARAM_in(JRCH)%R_MAN_N + + ! number of waves + NR = size(Q_JRCH) + + allocate(q_jrch_mod(0:NR-1), t_exit_mod(1:NR-1), wc(1:NR-1)) + + ! total "available" discharge in current time step + ! (zero position last routed; use of NR+1 instead of NR keeps next expected routed point) + TP = [T_START,T_END] + call interp_rch(TENTRY(0:NR-1),Q_JRCH(0:NR-1), TP, Qavg, ierr,cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + totQ = Qavg(1)*RPARAM_in(JRCH)%R_WIDTH + + if (abs(Qtake) < totQ) then + + ! modified wave Q [m3/s] + Qfrac = abs(Qtake)/TotQ ! fraction of total wave Q to target abstracted Q + Q_jrch_mod(0) = Q_JRCH(0) + Q_jrch_mod(1:NR-1) = Q_JRCH(1:NR-1)*(1._dp-Qfrac) ! remaining wave Q after abstraction + + else ! everything taken.... + + q_jrch_mod = runoffMin ! remaining wave Q after abstraction + + end if + + if(JRCH == ixDesire)then + ! compute actual abstracted water + allocate(Q_jrch_abs(0:NR-1)) + Q_jrch_abs = Q_JRCH - Q_jrch_mod + call interp_rch(TENTRY(0:NR-1),Q_jrch_abs(0:NR-1), TP, Qavg, ierr,cmessage) + Qabs = Qavg(1)*RPARAM_in(JRCH)%R_WIDTH + write(*,'(a)') ' * Target abstraction (Qtake) [m3/s], Available discharge (totQ) [m3/s], Actual abstraction (Qabs) [m3/s] ' + write(*,'(a,x,F15.7)') ' Qtake =', Qtake + write(*,'(a,x,F15.7)') ' totQ =', totQ + write(*,'(a,x,F15.7)') ' Qabs =', Qabs + end if + + ! modify wave speed at modified wave discharge and re-compute exit time + wc = alfa*K**(1._dp/alfa)*Q_jrch_mod**((alfa-1._dp)/alfa) ! modified wave celerity [m/s] + do iw = 1,NR-1 + t_exit_mod(iw) = min(RPARAM_in(JRCH)%RLENGTH/wc(iw)+TENTRY(iw) , huge(TENTRY(iw))) + end do + + ! modified final q_jrch and t_exit + Q_JRCH(0:NR-1) = q_jrch_mod(0:NR-1) + T_EXIT(1:NR-1) = t_exit_mod(1:NR-1) + + if(JRCH == ixDesire)then + write(fmt1,'(A,I5,A)') '(A,1X',NR,'(1X,E15.7))' + write(*,'(a)') ' * After abstracted: wave discharge (Q_JRCH) [m2/s], entry time (TENTRY) [s], and exit time (T_EXIT) [s]:' + write(*,fmt1) ' Q_JRCH=',(Q_JRCH(iw), iw=0,NR-1) + write(*,fmt1) ' TENTRY=',(TENTRY(iw), iw=0,NR-1) + write(*,fmt1) ' T_EXIT=',(T_EXIT(iw), iw=0,NR-1) + endif + + END SUBROUTINE extract_from_rch + ! ********************************************************************* ! subroutine: extract flow from the reaches upstream of JRCH ! ********************************************************************* - subroutine GETUSQ_RCH(IENS,JRCH,LAKEFLAG,T0,T1,ixDesire, & ! input + subroutine getusq_rch(IENS,JRCH,LAKEFLAG,T0,T1,ixDesire, & ! input NETOPO_in,RPARAM_in,RCHFLX_in, & ! input KROUTE_out, & ! inout Q_JRCH,TENTRY,T_EXIT,ierr,message, & ! output @@ -484,14 +608,9 @@ subroutine GETUSQ_RCH(IENS,JRCH,LAKEFLAG,T0,T1,ixDesire, & ! input ! T_EXIT(:): Vector of times flow particles are expected to exit reach JRCH ! ! ---------------------------------------------------------------------------------------- - ! Future revisions: - ! - ! (none planned) - ! - ! ---------------------------------------------------------------------------------------- USE globalData, only : LKTOPO ! Lake topology USE globalData, only : LAKFLX ! Lake fluxes - IMPLICIT NONE + implicit none ! Input integer(I4B), intent(in) :: IENS ! ensemble member integer(I4B), intent(in) :: JRCH ! reach to process @@ -505,46 +624,49 @@ subroutine GETUSQ_RCH(IENS,JRCH,LAKEFLAG,T0,T1,ixDesire, & ! input ! inout type(KREACH), intent(inout), allocatable :: KROUTE_out(:,:) ! reach state data ! Output - REAL(DP),allocatable, intent(out) :: Q_JRCH(:) ! merged (non-routed) flow in JRCH - REAL(DP),allocatable, intent(out) :: TENTRY(:) ! time flow particles entered JRCH - REAL(DP),allocatable, intent(out) :: T_EXIT(:) ! time flow is expected to exit JR + real(dp),allocatable, intent(out) :: Q_JRCH(:) ! merged (non-routed) flow in JRCH + real(dp),allocatable, intent(out) :: TENTRY(:) ! time flow particles entered JRCH + real(dp),allocatable, intent(out) :: T_EXIT(:) ! time flow is expected to exit JR integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! Local variables to hold the merged inputs to the downstream reach - INTEGER(I4B) :: ROFFSET ! retrospective offset due to rstep - REAL(DP) :: DT ! model time step - REAL(DP), allocatable :: QD(:) ! merged downstream flow - REAL(DP), allocatable :: TD(:) ! merged downstream time - INTEGER(I4B) :: ND ! # points shifted downstream - INTEGER(I4B) :: NJ ! # points in the JRCH reach - INTEGER(I4B) :: NK ! # points for routing (NJ+ND) - INTEGER(I4B) :: ILAK ! lake index + integer(i4b) :: ROFFSET ! retrospective offset due to rstep + real(dp) :: DT ! model time step + real(dp), allocatable :: QD(:) ! merged downstream flow + real(dp), allocatable :: TD(:) ! merged downstream time + integer(i4b) :: iw ! wave loop index + integer(i4b) :: ND ! # of waves routed from upstreams + integer(i4b) :: NJ ! # of waves in the JRCH reach + integer(i4b) :: NK ! # of waves for routing (NJ+ND) + integer(i4b) :: ILAK ! lake index + character(len=strLen) :: fmt1 ! format string character(len=strLen) :: cmessage ! error message for downwind routine - ! initialize error control - ierr=0; message='GETUSQ_RCH/' - ! ---------------------------------------------------------------------------------------- - ! (1) EXTRACT (AND MERGE) FLOW FROM UPSTREAM REACHES OR LAKE - ! ---------------------------------------------------------------------------------------- - ! define dt + + ierr=0; message='getusq_rch/' + + ! set the retrospective offset and model time step [sec] DT = (T1 - T0) - ! set the retrospective offset - IF (.NOT.PRESENT(RSTEP)) THEN + if (.not.present(RSTEP)) then ROFFSET = 0 - ELSE + else ROFFSET = RSTEP - END IF + end if + + ! ---------------------------------------------------------------------------------------- + ! (1) EXTRACT (AND MERGE) FLOW FROM UPSTREAM REACHES OR LAKE + ! ---------------------------------------------------------------------------------------- if (LAKEFLAG.EQ.1) then ! lakes are enabled ! get lake outflow and only lake outflow if reach is a lake outlet reach, else do as normal ILAK = NETOPO_in(JRCH)%LAKE_IX ! lake index if (ILAK.GT.0) then ! part of reach is in lake if (NETOPO_in(JRCH)%REACHIX.eq.LKTOPO(ILAK)%DREACHI) then ! we are in a lake outlet reach ND = 1 - ALLOCATE(QD(1),TD(1),STAT=IERR) + allocate(QD(1),TD(1),STAT=IERR) if(ierr/=0)then; message=trim(message)//'problem allocating array for QD and TD'; return; endif QD(1) = LAKFLX(IENS,ILAK)%LAKE_Q / RPARAM_in(JRCH)%R_WIDTH ! lake outflow per unit reach width TD(1) = T1 - DT*ROFFSET else - call QEXMUL_RCH(IENS,JRCH,T0,T1,ixDesire, & ! input + call qexmul_rch(IENS,JRCH,T0,T1,ixDesire, & ! input NETOPO_in,RPARAM_in,RCHFLX_in, & ! input KROUTE_out, & ! inout ND,QD,TD,ierr,cmessage, & ! output @@ -552,7 +674,7 @@ subroutine GETUSQ_RCH(IENS,JRCH,LAKEFLAG,T0,T1,ixDesire, & ! input if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif endif else - call QEXMUL_RCH(IENS,JRCH,T0,T1,ixDesire, & ! input + call qexmul_rch(IENS,JRCH,T0,T1,ixDesire, & ! input NETOPO_in,RPARAM_in,RCHFLX_in, & ! input KROUTE_out, & ! inout ND,QD,TD,ierr,cmessage, & ! output @@ -560,23 +682,30 @@ subroutine GETUSQ_RCH(IENS,JRCH,LAKEFLAG,T0,T1,ixDesire, & ! input if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif endif else ! lakes disabled - call QEXMUL_RCH(IENS,JRCH,T0,T1,ixDesire, & ! input + call qexmul_rch(IENS,JRCH,T0,T1,ixDesire, & ! input NETOPO_in,RPARAM_in,RCHFLX_in, & ! input KROUTE_out, & ! inout ND,QD,TD,ierr,cmessage, & ! output RSTEP) ! optional input if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - if(JRCH == ixDesire) print*, 'after QEXMUL_RCH: JRCH, ND, QD = ', JRCH, ND, QD + if(JRCH == ixDesire) then + write(fmt1,'(A,I5,A)') '(A,1X',ND,'(1X,F15.7))' + write(*,'(a)') ' * After qexmul_rch: # of routed wave from upstreams (ND) and wave discharge (QD) [m2/s]:' + write(*,'(A,x,I5)') ' ND=', ND + write(*,fmt1) ' QD=', (QD(iw), iw=1,ND) + end if endif + ! ---------------------------------------------------------------------------------------- ! (2) EXTRACT NON-ROUTED FLOW FROM THE REACH JRCH & APPEND TO THE FLOW JUST ROUTED D/S ! ---------------------------------------------------------------------------------------- ! check that the routing structure is associated - if(allocated(KROUTE_out).eqv..FALSE.)THEN + if(allocated(KROUTE_out).eqv..false.)THEN ierr=20; message='routing structure KROUTE_out is not associated'; return endif + ! check that the wave has been initialized - if (allocated(KROUTE_out(IENS,JRCH)%KWAVE).eqv..FALSE.) THEN + if (allocated(KROUTE_out(IENS,JRCH)%KWAVE).eqv..false.) THEN ! if not initialized, then set initial flow to first flow ! (this will only occur for a cold start in the case of no streamflow observations) allocate(KROUTE_out(IENS,JRCH)%KWAVE(0:0),STAT=IERR) @@ -584,31 +713,32 @@ subroutine GETUSQ_RCH(IENS,JRCH,LAKEFLAG,T0,T1,ixDesire, & ! input KROUTE_out(IENS,JRCH)%KWAVE(0)%QF = QD(1) KROUTE_out(IENS,JRCH)%KWAVE(0)%TI = T0 - DT - DT*ROFFSET KROUTE_out(IENS,JRCH)%KWAVE(0)%TR = T0 - DT*ROFFSET - KROUTE_out(IENS,JRCH)%KWAVE(0)%RF = .TRUE. + KROUTE_out(IENS,JRCH)%KWAVE(0)%RF = .true. endif + ! now extract the non-routed flow ! NB: routed flows were stripped out in the previous timestep when JRCH was index of u/s reach ! {only non-routed flows remain in the routing structure [ + zero element (last routed)]} - NJ = SIZE(KROUTE_out(IENS,JRCH)%KWAVE) - 1 ! number of elements not routed (-1 for 0) + NJ = size(KROUTE_out(IENS,JRCH)%KWAVE) - 1 ! number of elements not routed (-1 for 0) NK = NJ + ND ! pts still in reach + u/s pts just routed - ALLOCATE(Q_JRCH(0:NK),TENTRY(0:NK),T_EXIT(0:NK),STAT=IERR) ! include zero element for INTERP later + + allocate(Q_JRCH(0:NK),TENTRY(0:NK),T_EXIT(0:NK),STAT=IERR) ! include zero element for INTERP later if(ierr/=0)then; message=trim(message)//'problem allocating array for [Q_JRCH, TENTRY, T_EXIT]'; return; endif + Q_JRCH(0:NJ) = KROUTE_out(IENS,JRCH)%KWAVE(0:NJ)%QF ! extract the non-routed flow from reach JR TENTRY(0:NJ) = KROUTE_out(IENS,JRCH)%KWAVE(0:NJ)%TI ! extract the non-routed time from reach JR T_EXIT(0:NJ) = KROUTE_out(IENS,JRCH)%KWAVE(0:NJ)%TR ! extract the expected exit time Q_JRCH(NJ+1:NJ+ND) = QD(1:ND) ! append u/s flow just routed downstream TENTRY(NJ+1:NJ+ND) = TD(1:ND) ! append u/s time just routed downstream T_EXIT(NJ+1:NJ+ND) = -9999.0D0 ! set un-used T_EXIT to missing - deallocate(QD,TD,STAT=IERR) ! routed flow appended, no longer needed - if(ierr/=0)then; message=trim(message)//'problem deallocating array for QD and TD'; return; endif - end subroutine GETUSQ_RCH + END SUBROUTINE getusq_rch ! ********************************************************************* ! subroutine: extract flow from multiple reaches and merge into ! a single series ! ********************************************************************* - subroutine QEXMUL_RCH(IENS,JRCH,T0,T1,ixDesire, & ! input + SUBROUTINE qexmul_rch(IENS,JRCH,T0,T1,ixDesire, & ! input NETOPO_in,RPARAM_in,RCHFLX_in, & ! input KROUTE_out, & ! inout ND,QD,TD,ierr,message, & ! output @@ -646,12 +776,7 @@ subroutine QEXMUL_RCH(IENS,JRCH,T0,T1,ixDesire, & ! input ! TD(:): Vector of times flow particles entered reach JRCH (exited upstream reaches) ! ! ---------------------------------------------------------------------------------------- - ! Future revisions: - ! - ! (none planned) - ! - ! ---------------------------------------------------------------------------------------- - IMPLICIT NONE + implicit none ! Input INTEGER(i4b), intent(in) :: IENS ! ensemble member INTEGER(i4b), intent(in) :: JRCH ! reach to process @@ -670,59 +795,52 @@ subroutine QEXMUL_RCH(IENS,JRCH,T0,T1,ixDesire, & ! input integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! Local variables to hold flow/time from upstream reaches - REAL(DP) :: DT ! model time step - INTEGER(I4B) :: ROFFSET ! retrospective offset due to rstep - INTEGER(I4B) :: IUPS ! loop through u/s reaches - INTEGER(I4B) :: NUPB ! number of upstream basins - INTEGER(I4B) :: NUPR ! number of upstream reaches - INTEGER(I4B) :: INDX ! index of the IUPS u/s reach - INTEGER(I4B) :: MUPR ! # reaches u/s of IUPS u/s reach - INTEGER(I4B) :: NUPS ! number of upstream elements + real(dp) :: DT ! model time step + integer(i4b) :: ROFFSET ! retrospective offset due to rstep + integer(i4b) :: IUPS ! loop through u/s reaches + integer(i4b) :: NUPB ! number of upstream basins + integer(i4b) :: NUPR ! number of upstream reaches + integer(i4b) :: INDX ! index of the IUPS u/s reach + integer(i4b) :: MUPR ! # reaches u/s of IUPS u/s reach + integer(i4b) :: NUPS ! number of upstream elements TYPE(KREACH), allocatable :: USFLOW(:) ! waves for all upstream segments - REAL(DP), allocatable :: UWIDTH(:) ! width of all upstream segments - INTEGER(I4B) :: IMAX ! max number of upstream particles - INTEGER(I4B) :: IUPR ! counter for reaches with particles - INTEGER(I4B) :: IR ! index of the upstream reach - INTEGER(I4B) :: NS ! size of the wave - INTEGER(I4B) :: NR ! # routed particles in u/s reach - INTEGER(I4B) :: NQ ! NR+1, if non-routed particle exists + real(dp), allocatable :: UWIDTH(:) ! width of all upstream segments + integer(i4b) :: IMAX ! max number of upstream particles + integer(i4b) :: IUPR ! counter for reaches with particles + integer(i4b) :: IR ! index of the upstream reach + integer(i4b) :: NS ! size of the wave + integer(i4b) :: NR ! # routed particles in u/s reach + integer(i4b) :: NQ ! NR+1, if non-routed particle exists TYPE(FPOINT), allocatable :: NEW_WAVE(:) ! temporary wave - LOGICAL(LGT) :: INIT=.TRUE. ! used to initialize pointers ! Local variables to merge flow - LOGICAL(LGT), DIMENSION(:), ALLOCATABLE :: MFLG ! T = all particles processed - INTEGER(I4B), DIMENSION(:), ALLOCATABLE :: ITIM ! processing point for all u/s segments - REAL(DP), DIMENSION(:), ALLOCATABLE :: CTIME ! central time for each u/s segment - INTEGER(I4B) :: JUPS ! index of reach with the earliest time - REAL(DP) :: Q_AGG ! aggregarted flow at a given time - INTEGER(I4B) :: IWAV ! index of particle in the IUPS reach - REAL(DP) :: SCFAC ! scale to conform to d/s reach width - REAL(DP) :: SFLOW ! scaled flow at CTIME(JUPS) - INTEGER(I4B) :: IBEG,IEND ! indices for particles that bracket time - REAL(DP) :: SLOPE ! slope for the interpolation - REAL(DP) :: PREDV ! value predicted by the interpolation - INTEGER(I4B) :: IPRT ! counter for flow particles - INTEGER(I4B) :: JUPS_OLD ! check that we don't get stuck in do-forever - INTEGER(I4B) :: ITIM_OLD ! check that we don't get stuck in do-forever - REAL(DP) :: TIME_OLD ! previous time -- used to check for duplicates - REAL(DP), allocatable :: QD_TEMP(:)! flow particles just enetered JRCH - REAL(DP), allocatable :: TD_TEMP(:)! time flow particles entered JRCH - ! initialize error control - ierr=0; message='QEXMUL_RCH/' - ! ---------------------------------------------------------------------------------------- - ! (0) INITIALIZE POINTERS - ! ---------------------------------------------------------------------------------------- - IF(INIT) THEN - INIT=.FALSE. - !deallocate(USFLOW,NEW_WAVE,QD_TEMP,TD_TEMP) - ENDIF - ! set the retrospective offset - IF (.NOT.PRESENT(RSTEP)) THEN + logical(lgt), dimension(:), allocatable :: MFLG ! T = all particles processed + integer(i4b), dimension(:), allocatable :: ITIM ! processing point for all u/s segments + real(dp), dimension(:), allocatable :: CTIME ! central time for each u/s segment + integer(i4b) :: JUPS ! index of reach with the earliest time + real(dp) :: Q_AGG ! aggregarted flow at a given time + integer(i4b) :: IWAV ! index of particle in the IUPS reach + real(dp) :: SCFAC ! scale to conform to d/s reach width + real(dp) :: SFLOW ! scaled flow at CTIME(JUPS) + integer(i4b) :: IBEG,IEND ! indices for particles that bracket time + real(dp) :: SLOPE ! slope for the interpolation + real(dp) :: PREDV ! value predicted by the interpolation + integer(i4b) :: IPRT ! counter for flow particles + integer(i4b) :: JUPS_OLD ! check that we don't get stuck in do-forever + integer(i4b) :: ITIM_OLD ! check that we don't get stuck in do-forever + real(dp) :: TIME_OLD ! previous time -- used to check for duplicates + real(dp), allocatable :: QD_TEMP(:)! flow particles just enetered JRCH + real(dp), allocatable :: TD_TEMP(:)! time flow particles entered JRCH + + ierr=0; message='qexmul_rch/' + + ! set the retrospective offset and model time step [sec] + if (.not.PRESENT(RSTEP)) then ROFFSET = 0 - ELSE + else ROFFSET = RSTEP - END IF - ! define dt + end if DT = (T1 - T0) + ! ---------------------------------------------------------------------------------------- ! (1) DETERMINE THE NUMBER OF UPSTREAM REACHES ! ---------------------------------------------------------------------------------------- @@ -733,80 +851,89 @@ subroutine QEXMUL_RCH(IENS,JRCH,T0,T1,ixDesire, & ! input ! the number of series merged from upstream reaches is the number of upstream basins + ! the number of upstream reaches that are not headwater basins. NUPR = 0 ! number of upstream reaches - NUPB = SIZE(NETOPO_in(JRCH)%UREACHI) ! number of upstream basins + NUPB = size(NETOPO_in(JRCH)%UREACHI) ! number of upstream basins !NUPB = count(NETOPO_in(JRCH)%goodBas) ! number of upstream basins - DO IUPS=1,NUPB + do IUPS=1,NUPB INDX = NETOPO_in(JRCH)%UREACHI(IUPS) ! index of the IUPS upstream reach !MUPR = SIZE(NETOPO_in(INDX)%UREACHI) ! # reaches upstream of the IUPS upstream reach MUPR = count(NETOPO_in(INDX)%goodBas) ! # reaches upstream of the IUPS upstream reach - IF (MUPR.GT.0) NUPR = NUPR + 1 ! reach has streamflow in it, so add that as well - END DO ! iups + if (MUPR.GT.0) NUPR = NUPR + 1 ! reach has streamflow in it, so add that as well + end do ! iups NUPS = NUPB + NUPR ! number of upstream elements (basins + reaches) !print*, 'NUPB, NUPR, NUPS', NUPB, NUPR, NUPS !print*, 'NETOPO_in(JRCH)%UREACHK = ', NETOPO_in(JRCH)%UREACHK !print*, 'NETOPO_in(JRCH)%goodBas = ', NETOPO_in(JRCH)%goodBas - ! if nups eq 1, then ** SPECIAL CASE ** of just one upstream basin that is a headwater - IF (NUPS.EQ.1) THEN + + ! ** SPECIAL CASE ** of just one upstream basin that is a headwater + if (NUPS.EQ.1) then ND = 1 - ALLOCATE(QD(1),TD(1),STAT=IERR) + allocate(QD(1),TD(1),STAT=IERR) if(ierr/=0)then; message=trim(message)//'problem allocating array QD and TD'; return; endif ! get reach index IR = NETOPO_in(JRCH)%UREACHI(1) ! get flow in m2/s (scaled by with of downstream reach) QD(1) = RCHFLX_in(IENS,IR)%BASIN_QR(1)/RPARAM_in(JRCH)%R_WIDTH TD(1) = T1 - if(JRCH == ixDesire) print*, 'special case: JRCH, IR, NETOPO_in(IR)%REACHID = ', JRCH, IR, NETOPO_in(IR)%REACHID - RETURN - ENDIF + + if(JRCH == ixDesire) then + write(*,'(A,x,I8,x,I8)') ' * Special case - This reach has one headwater upstream: IR, NETOPO_in(IR)%REACHID = ', IR, NETOPO_in(IR)%REACHID + end if + + return + endif + ! allocate space for the upstream flow, time, and flags - ALLOCATE(USFLOW(NUPS),UWIDTH(NUPS),CTIME(NUPS),STAT=IERR) + allocate(USFLOW(NUPS),UWIDTH(NUPS),CTIME(NUPS),STAT=IERR) if(ierr/=0)then; message=trim(message)//'problem allocating arrays [USFLOW, UWIDTH, CTIME]'; return; endif ! define the minimum size of the routed data structure (number of flow particles) ! (IMAX is increased when looping through the reaches -- section 3 below) IMAX = NUPB ! flow from basins (one particle / timestep) + ! ---------------------------------------------------------------------------------------- ! (2) EXTRACT FLOW FROM UPSTREAM BASINS ! ---------------------------------------------------------------------------------------- - DO IUPS=1,NUPB + do IUPS=1,NUPB ! identify the index for the IUPS upstream segment IR = NETOPO_in(JRCH)%UREACHI(IUPS) ! allocate space for the IUPS stream segment (flow, time, and flags) - ALLOCATE(USFLOW(IUPS)%KWAVE(0:1),STAT=IERR) ! basin, has flow @start and @end of the time step + allocate(USFLOW(IUPS)%KWAVE(0:1),STAT=IERR) ! basin, has flow @start and @end of the time step if(ierr>0)then; message=trim(message)//'problem allocating array USFLOW(IUPS)%KWAVE'; return; endif ! place flow and time in the KWAVE array (routing done with time-delay histogram in TIMDEL_BAS.F90) USFLOW(IUPS)%KWAVE(0:1)%QF = RCHFLX_in(IENS,IR)%BASIN_QR(0:1) ! flow USFLOW(IUPS)%KWAVE(0:1)%TI = (/T0,T1/) - DT*ROFFSET ! entry time (not used) USFLOW(IUPS)%KWAVE(0:1)%TR = (/T0,T1/) - DT*ROFFSET ! exit time - USFLOW(IUPS)%KWAVE(0:1)%RF = .TRUE. ! routing flag + USFLOW(IUPS)%KWAVE(0:1)%RF = .true. ! routing flag !write(*,'(a,i4,1x,2(e20.10,1x))') 'IR, USFLOW(IUPS)%KWAVE(0:1)%QF = ', IR, USFLOW(IUPS)%KWAVE(0:1)%QF ! save the upstream width UWIDTH(IUPS) = 1.0D0 ! basin = unit width ! save the the time for the first particle in each reach CTIME(IUPS) = USFLOW(IUPS)%KWAVE(1)%TR ! central time - END DO ! (loop through upstream basins) + end do ! (loop through upstream basins) + ! ---------------------------------------------------------------------------------------- ! (3) EXTRACT FLOW FROM UPSTREAM REACHES ! ---------------------------------------------------------------------------------------- IUPR = 0 - DO IUPS=1,NUPB + do IUPS=1,NUPB INDX = NETOPO_in(JRCH)%UREACHI(IUPS) ! index of the IUPS upstream reach !MUPR = SIZE(NETOPO_in(INDX)%UREACHI) ! # reaches upstream of the IUPS upstream reach MUPR = count(NETOPO_in(INDX)%goodBas) ! # reaches upstream of the IUPS upstream reach - IF (MUPR.GT.0) THEN ! reach has streamflow in it, so add that as well + if (MUPR.GT.0) then ! reach has streamflow in it, so add that as well IUPR = IUPR + 1 ! identify the index for the IUPS upstream segment IR = NETOPO_in(JRCH)%UREACHI(IUPS) ! identify the size of the wave - NS = SIZE(KROUTE_out(IENS,IR)%KWAVE) + NS = size(KROUTE_out(IENS,IR)%KWAVE) ! identify number of routed flow elements in the IUPS upstream segment - NR = COUNT(KROUTE_out(IENS,IR)%KWAVE(:)%RF) + NR = count(KROUTE_out(IENS,IR)%KWAVE(:)%RF) ! include a non-routed point, if it exists NQ = MIN(NR+1,NS) ! allocate space for the IUPS stream segment (flow, time, and flags) - ALLOCATE(USFLOW(NUPB+IUPR)%KWAVE(0:NQ-1),STAT=IERR) ! (zero position = last routed) + allocate(USFLOW(NUPB+IUPR)%KWAVE(0:NQ-1),STAT=IERR) ! (zero position = last routed) if(ierr/=0)then; message=trim(message)//'problem allocating array USFLOW(NUPB+IUPR)%KWAVE(0:NQ-1)'; return; endif ! place data in the new arrays USFLOW(NUPB+IUPR)%KWAVE(0:NQ-1) = KROUTE_out(IENS,IR)%KWAVE(0:NQ-1) + ! here a statement where we check for a modification in the upstream reach; ! if flow upstream is modified, then copy KROUTE_out(:,:)%KWAVE(:)%QM to USFLOW(..)%KWAVE%QF !IF (NUSER.GT.0.AND.SIMDAT%UCFFLAG.GE.1) THEN !if the irrigation module is active and there are users @@ -815,30 +942,31 @@ subroutine QEXMUL_RCH(IENS,JRCH,T0,T1,ixDesire, & ! input ! USFLOW(NUPB+IUPR)%KWAVE(0:NQ-1)%QF = KROUTE_out(IENS,IR)%KWAVE(0:NQ-1)%QM ! ENDIF !ENDIF + ! ...and REMOVE the routed particles from the upstream reach ! (copy the wave to a temporary wave) - IF (allocated(NEW_WAVE)) THEN - DEALLOCATE(NEW_WAVE,STAT=IERR) ! (so we can allocate) + if (allocated(NEW_WAVE)) then + deallocate(NEW_WAVE,STAT=IERR) ! (so we can allocate) if(ierr/=0)then; message=trim(message)//'problem deallocating array NEW_WAVE'; return; endif - END IF - ALLOCATE(NEW_WAVE(0:NS-1),STAT=IERR) ! get new wave + end if + allocate(NEW_WAVE(0:NS-1),STAT=IERR) ! get new wave if(ierr/=0)then; message=trim(message)//'problem allocating array NEW_WAVE'; return; endif NEW_WAVE(0:NS-1) = KROUTE_out(IENS,IR)%KWAVE(0:NS-1) ! copy ! (re-size wave structure) - IF (.NOT.allocated(KROUTE_out(IENS,IR)%KWAVE))then; print*,' not allocated. in qex ';return; endif - IF (allocated(KROUTE_out(IENS,IR)%KWAVE)) THEN + if (.not.allocated(KROUTE_out(IENS,IR)%KWAVE))then; print*,' not allocated. in qex ';return; endif + if (allocated(KROUTE_out(IENS,IR)%KWAVE)) then deallocate(KROUTE_out(IENS,IR)%KWAVE,STAT=IERR) if(ierr/=0)then; message=trim(message)//'problem deallocating array KROUTE_out'; return; endif - END IF - ALLOCATE(KROUTE_out(IENS,IR)%KWAVE(0:NS-NR),STAT=IERR) ! reduced size + end if + allocate(KROUTE_out(IENS,IR)%KWAVE(0:NS-NR),STAT=IERR) ! reduced size if(ierr/=0)then; message=trim(message)//'problem allocating array KROUTE_out'; return; endif ! (copy "last routed" and "non-routed" elements) KROUTE_out(IENS,IR)%KWAVE(0:NS-NR) = NEW_WAVE(NR-1:NS-1) ! (de-allocate temporary wave) - DEALLOCATE(NEW_WAVE,STAT=IERR) + deallocate(NEW_WAVE,STAT=IERR) if(ierr/=0)then; message=trim(message)//'problem deallocating array NEW_WAVE'; return; endif ! save the upstream width @@ -848,8 +976,9 @@ subroutine QEXMUL_RCH(IENS,JRCH,T0,T1,ixDesire, & ! input ! keep track of the total number of points that must be routed downstream IMAX = IMAX + (NR-1) ! exclude zero point for the last routed - ENDIF ! if reach has particles in it - END DO ! iups + endif ! if reach has particles in it + end do ! iups + ! ---------------------------------------------------------------------------------------- ! (4) MERGE FLOW FROM MULTIPLE UPSTREAM REACHES ! ---------------------------------------------------------------------------------------- @@ -872,123 +1001,124 @@ subroutine QEXMUL_RCH(IENS,JRCH,T0,T1,ixDesire, & ! input ! ---------------------------------------------------------------------------------------- IPRT = 0 ! initialize counter for flow particles in the output array ! allocate space for the merged flow at the downstream reach - ALLOCATE(QD_TEMP(IMAX),TD_TEMP(IMAX),STAT=IERR) + allocate(QD_TEMP(IMAX),TD_TEMP(IMAX),STAT=IERR) if(ierr/=0)then; message=trim(message)//'problem allocating arrays [QD_TEMP, TD_TEMP]'; return; endif ! allocate positional arrays - ALLOCATE(MFLG(NUPS),ITIM(NUPS),STAT=IERR) + allocate(MFLG(NUPS),ITIM(NUPS),STAT=IERR) if(ierr/=0)then; message=trim(message)//'problem allocating arrays [MFLG, ITIM]'; return; endif + ! initalize the flag that defines whether all particles in a given reach are processed - MFLG(1:NUPS) = .FALSE. ! false until all particles are processed + MFLG(1:NUPS) = .false. ! false until all particles are processed ! initialize the search vector ITIM(1:NUPS) = 1 ! start with the first element of the wave ! initialize jups_old and itim_old (used to check we don't get stuck in the do-forever loop) JUPS_OLD = HUGE(JUPS_OLD) ITIM_OLD = HUGE(ITIM_OLD) - DO ! loop through all the times in the upstream reaches until no more routed flows + do ! loop through all the times in the upstream reaches until no more routed flows ! find the reach with the earliest time in all upstream reaches ! (NB: the time at the start of the timestep is the earliest possible time and ! the time at the end of the timestep is the latest possible time) JUPS = MINLOC(CTIME,DIM=1) ! JUPS = reach w/ earliest time ! check that we're not stuck in a continuous do loop - IF (JUPS.EQ.JUPS_OLD .AND. ITIM(JUPS).EQ.ITIM_OLD) THEN + if (JUPS.EQ.JUPS_OLD .and. ITIM(JUPS).EQ.ITIM_OLD) then ierr=20; message=trim(message)//'stuck in the continuous do-loop'; return - ENDIF + endif ! save jups and itim(jups) to check that we don't get stuck in a continuous do-loop JUPS_OLD = JUPS ITIM_OLD = ITIM(JUPS) ! check that there are still particles in the given reach that require processing - IF (.NOT.MFLG(JUPS)) THEN + if (.not.MFLG(JUPS)) then ! check that the particle in question is a particle routed (if not, then don't process) - IF (USFLOW(JUPS)%KWAVE(ITIM(JUPS))%RF.EQV..FALSE.) THEN - MFLG(JUPS) = .TRUE. ! if routing flag is false, then have already processed all particles + if (USFLOW(JUPS)%KWAVE(ITIM(JUPS))%RF.EQV..false.) then + MFLG(JUPS) = .true. ! if routing flag is false, then have already processed all particles CTIME(JUPS) = HUGE(SFLOW) ! largest possible number = ensure reach is not selected again ! the particle is in need of processing - ELSE + else ! define previous time - IF (IPRT.GE.1) THEN + if (IPRT.GE.1) then TIME_OLD = TD_TEMP(IPRT) - ELSE ! (if no particles, set to largest possible negative number) + else ! (if no particles, set to largest possible negative number) TIME_OLD = -HUGE(SFLOW) - END IF + end if ! check that the particles are being processed in the correct order - IF (CTIME(JUPS).LT.TIME_OLD) THEN + IF (CTIME(JUPS).LT.TIME_OLD) then ierr=30; message=trim(message)//'expect process in order of time'; return - ENDIF + endif ! don't process if time already exists - IF (CTIME(JUPS).NE.TIME_OLD) THEN + if (CTIME(JUPS).NE.TIME_OLD) then ! ------------------------------------------------------------------------------------- ! compute sum of scaled flow for all reaches Q_AGG = 0.0D0 - DO IUPS=1,NUPS + do IUPS=1,NUPS ! identify the element of the wave for the IUPS upstream reach IWAV = ITIM(IUPS) ! compute scale factor (scale upstream flow by width of downstream reach) SCFAC = UWIDTH(IUPS) / RPARAM_in(JRCH)%R_WIDTH ! case of the upstream reach with the minimum time (no interpolation required) - IF (IUPS.EQ.JUPS) THEN + if (IUPS.EQ.JUPS) then SFLOW = USFLOW(IUPS)%KWAVE(IWAV)%QF * SCFAC ! scaled flow ! case of all other upstream reaches (*** now, interpolate ***) - ELSE + else ! identify the elements that bracket the flow particle in the reach JUPS ! why .GE.? Why not .GT.?? IBEG = IWAV; IF (USFLOW(IUPS)%KWAVE(IBEG)%TR.GE.CTIME(JUPS)) IBEG=IWAV-1 IEND = IBEG+1 ! *** check the elements are ordered as we think *** ! test if we have bracketed properly - IF (USFLOW(IUPS)%KWAVE(IEND)%TR.LT.CTIME(JUPS) .OR. & - USFLOW(IUPS)%KWAVE(IBEG)%TR.GT.CTIME(JUPS)) THEN + if (USFLOW(IUPS)%KWAVE(IEND)%TR.LT.CTIME(JUPS) .or. & + USFLOW(IUPS)%KWAVE(IBEG)%TR.GT.CTIME(JUPS)) then ierr=40; message=trim(message)//'the times are not ordered as we assume'; return - ENDIF ! test for bracketing + endif ! test for bracketing ! estimate flow for the IUPS upstream reach at time CTIME(JUPS) SLOPE = (USFLOW(IUPS)%KWAVE(IEND)%QF - USFLOW(IUPS)%KWAVE(IBEG)%QF) / & (USFLOW(IUPS)%KWAVE(IEND)%TR - USFLOW(IUPS)%KWAVE(IBEG)%TR) PREDV = USFLOW(IUPS)%KWAVE(IBEG)%QF + SLOPE*(CTIME(JUPS)-USFLOW(IUPS)%KWAVE(IBEG)%TR) SFLOW = PREDV * SCFAC ! scaled flow - ENDIF ! (if interpolating) + endif ! (if interpolating) ! aggregate flow Q_AGG = Q_AGG + SFLOW - END DO ! looping through upstream elements + end do ! looping through upstream elements ! ------------------------------------------------------------------------------------- ! place Q_AGG and CTIME(JUPS) in the output arrays IPRT = IPRT + 1 QD_TEMP(IPRT) = Q_AGG TD_TEMP(IPRT) = CTIME(JUPS) - ENDIF ! (check that time doesn't already exist) + endif ! (check that time doesn't already exist) ! check if the particle just processed is the last element - IF (ITIM(JUPS).EQ.SIZE(USFLOW(JUPS)%KWAVE)-1) THEN ! -1 because of the zero element - MFLG(JUPS) = .TRUE. ! have processed all particles in a given u/s reach - CTIME(JUPS) = HUGE(SFLOW) ! largest possible number = ensure reach is not selected again - ELSE + if (ITIM(JUPS).EQ.size(USFLOW(JUPS)%KWAVE)-1) then ! -1 because of the zero element + MFLG(JUPS) = .true. ! have processed all particles in a given u/s reach + CTIME(JUPS) = huge(SFLOW) ! largest possible number = ensure reach is not selected again + else ITIM(JUPS) = ITIM(JUPS) + 1 ! move on to the next flow element CTIME(JUPS) = USFLOW(JUPS)%KWAVE(ITIM(JUPS))%TR ! save the time - ENDIF ! (check if particle is the last element) - ENDIF ! (check if the particle is a routed element) - ENDIF ! (check that there are still particles to process) + endif ! (check if particle is the last element) + endif ! (check if the particle is a routed element) + endif ! (check that there are still particles to process) ! if processed all particles in all upstream reaches, then EXIT - IF (COUNT(MFLG).EQ.NUPS) EXIT - END DO ! do-forever + IF (count(MFLG).EQ.NUPS) exit + end do ! do-forever + ! free up memory - DO IUPS=1,NUPS ! de-allocate each element of USFLOW - DEALLOCATE(USFLOW(IUPS)%KWAVE,STAT=IERR) + do IUPS=1,NUPS ! de-allocate each element of USFLOW + deallocate(USFLOW(IUPS)%KWAVE,STAT=IERR) if(ierr/=0)then; message=trim(message)//'problem deallocating array USFLOW(IUPS)%KWAVE'; return; endif - END DO ! looping thru elements of USFLOW - DEALLOCATE(USFLOW,UWIDTH,CTIME,ITIM,MFLG,STAT=IERR) + end do ! looping thru elements of USFLOW + deallocate(USFLOW,UWIDTH,CTIME,ITIM,MFLG,STAT=IERR) if(ierr/=0)then; message=trim(message)//'problem deallocating arrays [USFLOW, UWIDTH, CTIME, ITIM, MFLG]'; return; endif + ! ...and, save reduced arrays in QD and TD ND = IPRT - ALLOCATE(QD(ND),TD(ND),STAT=IERR) + allocate(QD(ND),TD(ND),STAT=IERR) if(ierr/=0)then; message=trim(message)//'problem allocating arrays [QD, TD]'; return; endif QD(1:ND) = QD_TEMP(1:ND) TD(1:ND) = TD_TEMP(1:ND) - DEALLOCATE(QD_TEMP,TD_TEMP,STAT=IERR) - if(ierr/=0)then; message=trim(message)//'problem deallocating arrays [QD_TEMP, TD_TEMP]'; return; endif - end subroutine QEXMUL_RCH + END SUBROUTINE qexmul_rch ! ********************************************************************* ! subroutine: removes flow particles from the routing structure, ! to reduce memory usage and processing time ! ********************************************************************* - subroutine REMOVE_RCH(MAXQPAR,& ! input + SUBROUTINE remove_rch(MAXQPAR,& ! input Q_JRCH,TENTRY,T_EXIT,ierr,message) ! output ! ---------------------------------------------------------------------------------------- ! Creator(s): @@ -1008,125 +1138,119 @@ subroutine REMOVE_RCH(MAXQPAR,& ! input ! T EXIT(:): Vector of times flow particles are EXPECTED to exit reach JRCH ! ! ---------------------------------------------------------------------------------------- - ! Future revisions: - ! - ! (none planned) - ! - ! ---------------------------------------------------------------------------------------- - IMPLICIT NONE + implicit none ! Input - INTEGER(I4B), INTENT(IN) :: MAXQPAR ! maximum number of flow particles allowed + integer(i4b), intent(in) :: MAXQPAR ! maximum number of flow particles allowed ! output - REAL(DP), allocatable, intent(inout) :: Q_JRCH(:)! merged (non-routed) flow in JRCH - REAL(DP), allocatable, intent(inout) :: TENTRY(:)! time flow particles entered JRCH - REAL(DP), allocatable, intent(inout) :: T_EXIT(:)! time flow particles exited JRCH + real(dp), allocatable, intent(inout) :: Q_JRCH(:)! merged (non-routed) flow in JRCH + real(dp), allocatable, intent(inout) :: TENTRY(:)! time flow particles entered JRCH + real(dp), allocatable, intent(inout) :: T_EXIT(:)! time flow particles exited JRCH integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! Local variables - INTEGER(I4B) :: NPRT ! number of flow particles - INTEGER(I4B) :: IPRT ! loop through flow particles - REAL(DP), DIMENSION(:), ALLOCATABLE :: Q,T,Z ! copies of Q_JRCH and T_JRCH - LOGICAL(LGT), DIMENSION(:), ALLOCATABLE :: PARFLG ! .FALSE. if particle removed - INTEGER(I4B), DIMENSION(:), ALLOCATABLE :: INDEX0 ! indices of original vectors - REAL(DP), DIMENSION(:), ALLOCATABLE :: ABSERR ! absolute error btw interp and orig - REAL(DP) :: Q_INTP ! interpolated particle - INTEGER(I4B) :: MPRT ! local number of flow particles - INTEGER(I4B), DIMENSION(:), ALLOCATABLE :: INDEX1 ! indices of particles retained - REAL(DP), DIMENSION(:), ALLOCATABLE :: E_TEMP ! temp abs error btw interp and orig - INTEGER(I4B), DIMENSION(1) :: ITMP ! result of minloc function - INTEGER(I4B) :: ISEL ! index of local minimum value - INTEGER(I4B) :: INEG ! lower boundary for interpolation - INTEGER(I4B) :: IMID ! desired point for interpolation - INTEGER(I4B) :: IPOS ! upper boundary for interpolation - ! initialize error control - ierr=0; message='REMOVE_RCH/' + integer(i4b) :: NPRT ! number of flow particles + integer(i4b) :: IPRT ! loop through flow particles + real(dp), dimension(:), allocatable :: Q,T,Z ! copies of Q_JRCH and T_JRCH + logical(lgt), dimension(:), allocatable :: PARFLG ! .FALSE. if particle removed + integer(i4b), dimension(:), allocatable :: INDEX0 ! indices of original vectors + real(dp), dimension(:), allocatable :: ABSERR ! absolute error btw interp and orig + real(dp) :: Q_INTP ! interpolated particle + integer(i4b) :: MPRT ! local number of flow particles + integer(i4b), dimension(:), allocatable :: INDEX1 ! indices of particles retained + real(dp), dimension(:), allocatable :: E_TEMP ! temp abs error btw interp and orig + integer(i4b), dimension(1) :: ITMP ! result of minloc function + integer(i4b) :: ISEL ! index of local minimum value + integer(i4b) :: INEG ! lower boundary for interpolation + integer(i4b) :: IMID ! desired point for interpolation + integer(i4b) :: IPOS ! upper boundary for interpolation + + ierr=0; message='remove_rch/' + ! ---------------------------------------------------------------------------------------- ! (1) INITIALIZATION ! ---------------------------------------------------------------------------------------- ! get the number of particles - NPRT = SIZE(Q_JRCH)-1 ! -1 because of zero element + NPRT = size(Q_JRCH)-1 ! -1 because of zero element ! allocate and initialize arrays - ALLOCATE(Q(0:NPRT),T(0:NPRT),Z(0:NPRT),PARFLG(0:NPRT),INDEX0(0:NPRT),ABSERR(0:NPRT),STAT=IERR) + allocate(Q(0:NPRT),T(0:NPRT),Z(0:NPRT),PARFLG(0:NPRT),INDEX0(0:NPRT),ABSERR(0:NPRT),STAT=IERR) if(ierr/=0)then; message=trim(message)//'problem allocating arrays [Q, T, Z, PARFLG, INDEX0, ABSERR]'; return; endif Q = Q_JRCH; T = TENTRY ! get copies of Q_JRCH and TENTRY Z = T_EXIT ! (not used in the interp, but include for consistency) - PARFLG = .TRUE. ! particle flag = start with all points + PARFLG = .true. ! particle flag = start with all points INDEX0 = arth(0,1,NPRT+1) ! index = (0,1,2,...,NPRT) ABSERR = HUGE(Q) ! largest possible double-precision number ! get the absolte difference between actual points and interpolated points - DO IPRT=1,NPRT-1 + do IPRT=1,NPRT-1 ! interpolate at point (iprt) Q_INTP = INTERP(T(IPRT),Q(IPRT-1),Q(IPRT+1),T(IPRT-1),T(IPRT+1)) ! save the absolute difference between the actual value and the interpolated value - ABSERR(IPRT) = ABS(Q_INTP-Q(IPRT)) - END DO + ABSERR(IPRT) = abs(Q_INTP-Q(IPRT)) + end do + ! ---------------------------------------------------------------------------------------- ! (2) REMOVAL ! ---------------------------------------------------------------------------------------- - DO ! continue looping until the number of particles is below the limit + do ! continue looping until the number of particles is below the limit ! get the number of particles still in the structure - MPRT = COUNT(PARFLG)-1 ! -1 because of the zero element + MPRT = count(PARFLG)-1 ! -1 because of the zero element ! get a copy of (1) indices of selected points, and (2) the interpolation errors - ALLOCATE(INDEX1(0:MPRT),E_TEMP(0:MPRT),STAT=IERR) + allocate(INDEX1(0:MPRT),E_TEMP(0:MPRT),STAT=IERR) if(ierr/=0)then; message=trim(message)//'problem allocating arrays [INDEX1, E_TEMP]'; return; endif - INDEX1 = PACK(INDEX0,PARFLG) ! (restrict attention to the elements still present) - E_TEMP = PACK(ABSERR,PARFLG) + INDEX1 = pack(INDEX0,PARFLG) ! (restrict attention to the elements still present) + E_TEMP = pack(ABSERR,PARFLG) ! check for exit condition (exit after "pack" b.c. indices used to construct final vectors) - IF (MPRT.LT.MAXQPAR) EXIT + if (MPRT.LT.MAXQPAR) exit ! get the index of the minimum value - ITMP = MINLOC(E_TEMP) - ISEL = LBOUND(E_TEMP,DIM=1) + ITMP(1) - 1 ! MINLOC assumes count from 1, here (0,1,2,...NPRT) + ITMP = minloc(E_TEMP) + ISEL = lbound(E_TEMP,dim=1) + ITMP(1) - 1 ! MINLOC assumes count from 1, here (0,1,2,...NPRT) ! re-interpolate the point immediately before the point flagged for removal - IF (INDEX1(ISEL-1).GT.0) THEN + if (INDEX1(ISEL-1).GT.0) then INEG=INDEX1(ISEL-2); IMID=INDEX1(ISEL-1); IPOS=INDEX1(ISEL+1) Q_INTP = INTERP(T(IMID),Q(INEG),Q(IPOS),T(INEG),T(IPOS)) - ABSERR(IMID) = ABS(Q_INTP-Q(IMID)) - ENDIF + ABSERR(IMID) = abs(Q_INTP-Q(IMID)) + endif ! re-interpolate the point immediately after the point flagged for removal - IF (INDEX1(ISEL+1).LT.NPRT) THEN + if (INDEX1(ISEL+1).LT.NPRT) then INEG=INDEX1(ISEL-1); IMID=INDEX1(ISEL+1); IPOS=INDEX1(ISEL+2) Q_INTP = INTERP(T(IMID),Q(INEG),Q(IPOS),T(INEG),T(IPOS)) - ABSERR(IMID) = ABS(Q_INTP-Q(IMID)) - ENDIF + ABSERR(IMID) = abs(Q_INTP-Q(IMID)) + endif ! flag the point as "removed" - PARFLG(INDEX1(ISEL)) = .FALSE. + PARFLG(INDEX1(ISEL)) = .false. ! de-allocate arrays - DEALLOCATE(INDEX1,E_TEMP,STAT=IERR) + deallocate(INDEX1,E_TEMP,STAT=IERR) if(ierr/=0)then; message=trim(message)//'problem deallocating arrays [INDEX1, E_TEMP]'; return; endif - END DO ! keep looping until a sufficient number of points are removed + end do ! keep looping until a sufficient number of points are removed + ! ---------------------------------------------------------------------------------------- ! (3) RE-SIZE DATA STRUCTURES ! ---------------------------------------------------------------------------------------- - DEALLOCATE(Q_JRCH,TENTRY,T_EXIT,STAT=IERR) + deallocate(Q_JRCH,TENTRY,T_EXIT,STAT=IERR) if(ierr/=0)then; message=trim(message)//'problem deallocating arrays [Q_JRCH, TENTRY, T_EXIT]'; return; endif - ALLOCATE(Q_JRCH(0:MPRT),TENTRY(0:MPRT),T_EXIT(0:MPRT),STAT=IERR) + allocate(Q_JRCH(0:MPRT),TENTRY(0:MPRT),T_EXIT(0:MPRT),STAT=IERR) if(ierr/=0)then; message=trim(message)//'problem allocating arrays [Q_JRCH, TENTRY, T_EXIT]'; return; endif Q_JRCH = Q(INDEX1) TENTRY = T(INDEX1) T_EXIT = Z(INDEX1) - DEALLOCATE(INDEX1,E_TEMP,STAT=IERR) - if(ierr/=0)then; message=trim(message)//'problem deallocating arrays [INDEX1, E_TEMP]'; return; endif - DEALLOCATE(Q,T,Z,PARFLG,ABSERR,INDEX0,STAT=IERR) - if(ierr/=0)then; message=trim(message)//'problem deallocating arrays [Q, T, Z, PARFLG, INDEX0, ABSERR]'; return; endif - contains + CONTAINS function INTERP(T0,Q1,Q2,T1,T2) - REAL(DP),INTENT(IN) :: Q1,Q2 ! flow at neighbouring times - REAL(DP),INTENT(IN) :: T1,T2 ! neighbouring times - REAL(DP),INTENT(IN) :: T0 ! desired time - REAL(DP) :: INTERP ! function name + real(dp),intent(in) :: Q1,Q2 ! flow at neighbouring times + real(dp),intent(in) :: T1,T2 ! neighbouring times + real(dp),intent(in) :: T0 ! desired time + real(dp) :: INTERP ! function name INTERP = Q1 + ( (Q2-Q1) / (T2-T1) ) * (T0-T1) end function INTERP - end subroutine + END SUBROUTINE ! ********************************************************************* ! new subroutine: calculate the propagation of kinematic waves in a ! single stream segment, including the formation and ! propagation of a kinematic shock ! ********************************************************************* - subroutine KINWAV_RCH(JRCH,T_START,T_END,ixDesire, & ! input: location and time + SUBROUTINE kinwav_rch(JRCH,T_START,T_END,ixDesire, & ! input: location and time NETOPO_in, RPARAM_in, & ! input: river data structure Q_JRCH,TENTRY,T_EXIT,FROUTE, & ! inout: kwt states NQ2, & ! output: @@ -1164,9 +1288,9 @@ subroutine KINWAV_RCH(JRCH,T_START,T_END,ixDesire, & ! input: loca ! Outputs: ! -------- ! FROUTE: array of routing flags -- All inputs are .FALSE., but flags change to .TRUE. - ! if element is routed INTENT(OUT) + ! if element is routed intent(out) ! T_EXIT: array of time elements -- identify the time each element is EXPECTED to exit - ! the stream segment, INTENT(OUT). Used in INTERPTS + ! the stream segment, intent(out). Used in INTERPTS ! NQ2: number of particles -- <= input becuase multiple particles may merge ! ! ---------------------------------------------------------------------------------------- @@ -1193,7 +1317,7 @@ subroutine KINWAV_RCH(JRCH,T_START,T_END,ixDesire, & ! input: loca ! ---------------------------------------------------------------------------------------- ! Modifications to source (mclark@ucar.edu): ! - ! * All variables are now defined (IMPLICIT NONE) and described (comments) + ! * All variables are now defined (implicit none) and described (comments) ! ! * Parameters are defined within the subroutine (for ease of readibility) ! @@ -1204,12 +1328,7 @@ subroutine KINWAV_RCH(JRCH,T_START,T_END,ixDesire, & ! input: loca ! and use F90 dynamic memory features ! ! ---------------------------------------------------------------------------------------- - ! Future revisions: - ! - ! (none planned) - ! - ! ---------------------------------------------------------------------------------------- - IMPLICIT NONE + implicit none ! Input integer(i4b), intent(in) :: JRCH ! Reach to process real(dp), intent(in) :: T_START ! start of the time step @@ -1227,30 +1346,31 @@ subroutine KINWAV_RCH(JRCH,T_START,T_END,ixDesire, & ! input: loca integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! Internal - REAL(DP) :: ALFA ! constant, 5/3 - REAL(DP) :: K ! sqrt(slope)/mannings N - REAL(DP) :: XMX ! length of the stream segment - INTEGER(I4B) :: NN ! number of input points - INTEGER(I4B) :: NI ! original size of the input - INTEGER(I4B) :: NM ! mumber of merged elements - INTEGER(I4B), DIMENSION(SIZE(Q_JRCH)) :: IX ! minimum index of each merged element - INTEGER(I4B), DIMENSION(SIZE(Q_JRCH)) :: MF ! index for input element merged - REAL(DP), DIMENSION(SIZE(Q_JRCH)) :: T0,T1,T2 ! copy of input time - REAL(DP), DIMENSION(SIZE(Q_JRCH)) :: Q0,Q1,Q2 ! flow series - REAL(DP), DIMENSION(SIZE(Q_JRCH)) :: WC ! wave celerity - INTEGER(I4B) :: IW,JW ! looping variables, break check - REAL(DP) :: X,XB ! define smallest, biggest shock - REAL(DP) :: WDIFF ! difference in wave celerity-1 - REAL(DP) :: XXB ! wave break - INTEGER(I4B) :: IXB,JXB ! define position of wave break - REAL(DP) :: A1,A2 ! stage - different sides of break - REAL(DP) :: CM ! merged celerity - REAL(DP) :: TEXIT ! expected exit time of "current" particle - REAL(DP) :: TNEXT ! expected exit time of "next" particle - REAL(DP) :: TEXIT2 ! exit time of "bottom" of merged element - INTEGER(I4B) :: IROUTE ! looping variable for routing - INTEGER(I4B) :: JROUTE ! looping variable for routing - INTEGER(I4B) :: ICOUNT ! used to account for merged pts + real(dp) :: ALFA ! constant, 5/3 + real(dp) :: K ! sqrt(slope)/mannings N + real(dp) :: XMX ! length of the stream segment + integer(i4b) :: NN ! number of input points + integer(i4b) :: NI ! original size of the input + integer(i4b) :: NM ! mumber of merged elements + integer(i4b), dimension(size(Q_JRCH)) :: IX ! minimum index of each merged element + integer(i4b), dimension(size(Q_JRCH)) :: MF ! index for input element merged + real(dp), dimension(size(Q_JRCH)) :: T0,T1,T2 ! copy of input time + real(dp), dimension(size(Q_JRCH)) :: Q0,Q1,Q2 ! flow series + real(dp), dimension(size(Q_JRCH)) :: WC ! wave celerity + integer(i4b) :: IW,JW ! looping variables, break check + real(dp) :: X,XB ! define smallest, biggest shock + real(dp) :: WDIFF ! difference in wave celerity-1 + real(dp) :: XXB ! wave break + integer(i4b) :: IXB,JXB ! define position of wave break + real(dp) :: A1,A2 ! stage - different sides of break + real(dp) :: CM ! merged celerity + real(dp) :: TEXIT ! expected exit time of "current" particle + real(dp) :: TNEXT ! expected exit time of "next" particle + real(dp) :: TEXIT2 ! exit time of "bottom" of merged element + integer(i4b) :: IROUTE ! looping variable for routing + integer(i4b) :: JROUTE ! looping variable for routing + integer(i4b) :: ICOUNT ! used to account for merged pts + character(len=strLen) :: fmt1 ! format string character(len=strLen) :: cmessage ! error message of downwind routine ! ---------------------------------------------------------------------------------------- ! NOTE: If merged particles DO NOT exit the reach in the current time step, they are @@ -1271,16 +1391,19 @@ subroutine KINWAV_RCH(JRCH,T_START,T_END,ixDesire, & ! input: loca ! the 2nd, 3rd, and 4th elements of MF = 2), populate the output vector with the ! selected elements (2,3,4) of the input vector. ! ---------------------------------------------------------------------------------------- - ! initialize error control - ierr=0; message='KINWAV_RCH/' + + ierr=0; message='kinwav_rch/' + ! Get the reach parameters ALFA = 5._dp/3._dp ! should this be initialized here or in a parameter file? - K = SQRT(RPARAM_in(JRCH)%R_SLOPE)/RPARAM_in(JRCH)%R_MAN_N + K = sqrt(RPARAM_in(JRCH)%R_SLOPE)/RPARAM_in(JRCH)%R_MAN_N XMX = RPARAM_in(JRCH)%RLENGTH + ! Identify the number of points to route - NN = SIZE(Q1) ! modified when elements are merged + NN = size(Q1) ! modified when elements are merged NI = NN ! original size of the input - IF(NN.EQ.0) RETURN ! don't do anything if no points in the reach + if(NN.EQ.0) return ! don't do anything if no points in the reach + ! Initialize the vector that indicates which output element the input elements are merged MF = arth(1,1,NI) ! Num. Rec. intrinsic: see MODULE nrutil.f90 ! Initialize the vector that indicates the minumum index of each merged element @@ -1290,32 +1413,37 @@ subroutine KINWAV_RCH(JRCH,T_START,T_END,ixDesire, & ! input: loca T0=TENTRY; T1=TENTRY; T2=TENTRY ! compute wave celerity for all flow points (array operation) WC(1:NN) = ALFA*K**(1./ALFA)*Q1(1:NN)**((ALFA-1.)/ALFA) - ! check - if(jRch==ixDesire) print*, 'q1(1:nn), wc(1:nn), RPARAM_in(JRCH)%R_SLOPE, nn = ', & - q1(1:nn), wc(1:nn), RPARAM_in(JRCH)%R_SLOPE, nn + + if(jRch==ixDesire) then + write(fmt1,'(A,I5,A)') '(A,1X',NN,'(1X,F15.7))' + write(*,'(a)') ' * Wave discharge (q1) [m2/s] and wave celertiy (wc) [m/s]:' + write(*,'(a,x,I3)') ' Number of wave =', NN + write(*,fmt1) ' q1=', (q1(iw), iw=1,NN) + write(*,fmt1) ' wc=', (wc(iw), iw=1,NN) + end if ! handle breaking waves - GT_ONE: IF(NN.GT.1) THEN ! no breaking if just one point + GT_ONE: if(NN.GT.1) then ! no breaking if just one point X = 0. ! altered later to describe "closest" shock - GOTALL: DO ! keep going until all shocks are merged + GOTALL: do ! keep going until all shocks are merged XB = XMX ! initialized to length of the stream segment ! -------------------------------------------------------------------------------------- ! check for breaking ! -------------------------------------------------------------------------------------- - WCHECK: DO IW=2,NN + WCHECK: do IW=2,NN JW=IW-1 - IF(WC(IW).EQ.0. .OR. WC(JW).EQ.0.) CYCLE ! waves not moving + if(WC(IW).EQ.0. .or. WC(JW).EQ.0.) cycle ! waves not moving WDIFF = 1./WC(JW) - 1./WC(IW) ! difference in wave celerity - IF(WDIFF.EQ.0.) CYCLE ! waves moving at the same speed - IF(WC(IW).EQ.WC(JW)) CYCLE ! identical statement to the above? + if(WDIFF.EQ.0.) cycle ! waves moving at the same speed + if(WC(IW).EQ.WC(JW)) cycle ! identical statement to the above? XXB = (T1(IW)-T1(JW)) / WDIFF ! XXB is point of breaking in x direction - IF(XXB.LT.X .OR. XXB.GT.XB) CYCLE ! XB init at LENGTH, so > XB do in next reach + if(XXB.LT.X .or. XXB.GT.XB) cycle ! XB init at LENGTH, so > XB do in next reach ! if get to here, the wave is breaking XB = XXB ! identify break "closest to upstream" first IXB = IW - END DO WCHECK + end do WCHECK ! -------------------------------------------------------------------------------------- - IF(XB.EQ.XMX) EXIT ! got all breaking waves, exit gotall + if (XB.EQ.XMX) exit ! got all breaking waves, exit gotall ! -------------------------------------------------------------------------------------- ! combine waves ! -------------------------------------------------------------------------------------- @@ -1323,8 +1451,8 @@ subroutine KINWAV_RCH(JRCH,T_START,T_END,ixDesire, & ! input: loca JXB = IXB-1 ! indices for the point of breaking NM = NI-NN ! number of merged elements ! calculate merged shockwave celerity (CM) using finite-difference approximation - Q2(JXB) =MAX(Q2(JXB),Q2(IXB)) ! flow of largest merged point - Q1(JXB) =MIN(Q1(JXB),Q1(IXB)) ! flow of smallest merged point + Q2(JXB) =max(Q2(JXB),Q2(IXB)) ! flow of largest merged point + Q1(JXB) =min(Q1(JXB),Q1(IXB)) ! flow of smallest merged point A2 = (Q2(JXB)/K)**(1./ALFA) ! Q = (1./MAN_N) H**(ALFA) sqrt(SLOPE) A1 = (Q1(JXB)/K)**(1./ALFA) ! H = (Q/K)**(1./ALFA) (K=sqrt(SLOPE)/MAN_N) CM = (Q2(JXB)-Q1(JXB))/(A2-A1) ! NB: A1,A2 are river stage @@ -1342,68 +1470,73 @@ subroutine KINWAV_RCH(JRCH,T_START,T_END,ixDesire, & ! input: loca ! update X - already got the "closest shock to start", see if there are any other shocks X = XB ! -------------------------------------------------------------------------------------- - END DO GOTALL - ENDIF GT_ONE + end do GOTALL + endif GT_ONE + + ! check + if(jRch==ixDesire) then + write(fmt1,'(A,I5,A)') '(A,1X',NN,'(1X,F15.7))' + write(*,'(a)') ' * After wave merge: wave celertiy (wc) [m/s]:' + write(*,'(a,x,I3)') ' Number of wave =', NN + write(*,fmt1) ' wc=', (wc(iw), iw=1,NN) + end if ICOUNT=0 ! ---------------------------------------------------------------------------------------- ! perform the routing ! ---------------------------------------------------------------------------------------- - DO IROUTE = 1,NN ! loop through the remaining particles (shocks,waves) (NM=NI-NN have been merged) - ! check - if(jRch==ixDesire) print*, 'wc(iRoute), nn = ', wc(iRoute), nn + do IROUTE = 1,NN ! loop through the remaining particles (shocks,waves) (NM=NI-NN have been merged) ! check that we have non-zero flow if(WC(IROUTE) < verySmall)then write(message,'(a,i0)') trim(message)//'zero flow for reach id ', NETOPO_in(jRch)%REACHID ierr=20; return endif ! compute the time the shock will exit the reach - TEXIT = MIN(XMX/WC(IROUTE) + T1(IROUTE), HUGE(T1)) + TEXIT = min(XMX/WC(IROUTE) + T1(IROUTE), huge(T1)) ! compute the time the next shock will exit the reach - IF (IROUTE.LT.NN) TNEXT = MIN(XMX/WC(IROUTE+1) + T1(IROUTE+1), HUGE(T1)) - IF (IROUTE.EQ.NN) TNEXT = HUGE(T1) + if (IROUTE.LT.NN) TNEXT = min(XMX/WC(IROUTE+1) + T1(IROUTE+1), huge(T1)) + if (IROUTE.EQ.NN) TNEXT = huge(T1) ! check if element is merged - MERGED: IF(Q1(IROUTE).NE.Q2(IROUTE)) THEN + MERGED: if (Q1(IROUTE).NE.Q2(IROUTE)) then ! check if merged element has exited - IF(TEXIT.LT.T_END) THEN + if (TEXIT.LT.T_END) then ! when a merged element exits, save just the top and the bottom of the shock ! (identify the exit time for the "slower" particle) - TEXIT2 = MIN(TEXIT+1.0D0, TEXIT + 0.5D0*(MIN(TNEXT,T_END)-TEXIT)) + TEXIT2 = min(TEXIT+1.0D0, TEXIT + 0.5D0*(min(TNEXT,T_END)-TEXIT)) ! unsure what will happen in the rare case if TEXIT and TEXIT2 are the same - IF (TEXIT2.EQ.TEXIT) THEN + if (TEXIT2.EQ.TEXIT) then ierr=30; message=trim(message)//'TEXIT equals TEXIT2 in kinwav'; return - ENDIF + end if ! fill output arrays - CALL RUPDATE(Q1(IROUTE),T1(IROUTE),TEXIT,ierr,cmessage) ! fill arrays w/ Q1, T1, + run checks + call RUPDATE(Q1(IROUTE),T1(IROUTE),TEXIT,ierr,cmessage) ! fill arrays w/ Q1, T1, + run checks if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - CALL RUPDATE(Q2(IROUTE),T1(IROUTE),TEXIT2,ierr,cmessage) ! fill arrays w/ Q2, T1, + run checks + call RUPDATE(Q2(IROUTE),T1(IROUTE),TEXIT2,ierr,cmessage) ! fill arrays w/ Q2, T1, + run checks if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - ELSE ! merged elements have not exited + else ! merged elements have not exited ! when a merged element does not exit, need to disaggregate into original particles - DO JROUTE=1,NI ! loop thru # original inputs - IF(MF(JROUTE).EQ.IROUTE) & - CALL RUPDATE(Q0(JROUTE),T0(JROUTE),TEXIT,ierr,cmessage) ! fill arrays w/ Q0, T0, + run checks + do JROUTE=1,NI ! loop thru # original inputs + if (MF(JROUTE).EQ.IROUTE) & + call RUPDATE(Q0(JROUTE),T0(JROUTE),TEXIT,ierr,cmessage) ! fill arrays w/ Q0, T0, + run checks if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - END DO ! JROUTE - ENDIF ! TEXIT + end do ! JROUTE + end if ! TEXIT ! now process un-merged particles - ELSE MERGED ! (i.e., not merged) - CALL RUPDATE(Q1(IROUTE),T1(IROUTE),TEXIT,ierr,cmessage) ! fill arrays w/ Q1, T1, + run checks + else MERGED ! (i.e., not merged) + call RUPDATE(Q1(IROUTE),T1(IROUTE),TEXIT,ierr,cmessage) ! fill arrays w/ Q1, T1, + run checks if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - ENDIF MERGED - END DO + end if MERGED + end do ! update arrays NQ2 = ICOUNT - RETURN - contains + CONTAINS - subroutine RUPDATE(QNEW,TOLD,TNEW,ierr,message) - REAL(DP),INTENT(IN) :: QNEW ! Q0,Q1, or Q2 - REAL(DP),INTENT(IN) :: TOLD,TNEW ! entry/exit times + SUBROUTINE RUPDATE(QNEW,TOLD,TNEW,ierr,message) + real(dp),intent(in) :: QNEW ! Q0,Q1, or Q2 + real(dp),intent(in) :: TOLD,TNEW ! entry/exit times integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message - ! initialize error control + ierr=0; message='RUPDATE/' ! --------------------------------------------------------------------------------------- ! Used to compute the time each element will exit stream segment & update routing flag @@ -1411,29 +1544,29 @@ subroutine RUPDATE(QNEW,TOLD,TNEW,ierr,message) ! --------------------------------------------------------------------------------------- ICOUNT=ICOUNT+1 ! check for array bounds exceeded - IF (ICOUNT.GT.SIZE(Q_JRCH)) THEN + if (ICOUNT.GT.size(Q_JRCH)) then ierr=60; message=trim(message)//'array bounds exceeded'; return - ENDIF + endif ! fill output arrays Q_JRCH(ICOUNT) = QNEW ! flow (Q1 always smaller than Q2) TENTRY(ICOUNT) = TOLD ! time - note, T1 altered if element merged T_EXIT(ICOUNT) = TNEW ! time check -- occurs when disaggregating merged elements - IF (ICOUNT.GT.1) THEN - IF (T_EXIT(ICOUNT).LE.T_EXIT(ICOUNT-1)) T_EXIT(ICOUNT)=T_EXIT(ICOUNT-1)+1. - ENDIF + if (ICOUNT.GT.1) then + if (T_EXIT(ICOUNT).LE.T_EXIT(ICOUNT-1)) T_EXIT(ICOUNT)=T_EXIT(ICOUNT-1)+1. + end if ! another time check -- rare problem when the shock can get the same time as tstart - IF(ICOUNT.EQ.1.AND.T_EXIT(ICOUNT).LE.T_START) T_EXIT(ICOUNT)=T_START+1. + if (ICOUNT.EQ.1.and.T_EXIT(ICOUNT).LE.T_START) T_EXIT(ICOUNT)=T_START+1. ! update flag for routed elements - IF(T_EXIT(ICOUNT).LT.T_END) FROUTE(ICOUNT) =.TRUE. - end subroutine RUPDATE + if (T_EXIT(ICOUNT).LT.T_END) FROUTE(ICOUNT) =.true. + END SUBROUTINE RUPDATE - end subroutine KINWAV_RCH + END SUBROUTINE kinwav_rch ! ********************************************************************* ! new subroutine: calculate time-step averages from irregular values ! ********************************************************************* - subroutine INTERP_RCH(TOLD,QOLD,TNEW,QNEW,IERR,MESSAGE) + SUBROUTINE interp_rch(TOLD,QOLD,TNEW,QNEW,IERR,MESSAGE) ! ---------------------------------------------------------------------------------------- ! Creator(s): ! Unknown (original Tideda routine?), fairly old @@ -1460,11 +1593,6 @@ subroutine INTERP_RCH(TOLD,QOLD,TNEW,QNEW,IERR,MESSAGE) ! IERR: error code (1 = bad bounds) ! ! ---------------------------------------------------------------------------------------- - ! Structures Used: - ! - ! NONE - ! - ! ---------------------------------------------------------------------------------------- ! Method: ! ! Loop through all output times (can be just 2, start and end of the time step) @@ -1503,125 +1631,120 @@ subroutine INTERP_RCH(TOLD,QOLD,TNEW,QNEW,IERR,MESSAGE) ! ---------------------------------------------------------------------------------------- ! Modifications to source (mclark@ucar.edu): ! - ! * All variables are now defined (IMPLICIT NONE) and described (comments) + ! * All variables are now defined (implicit none) and described (comments) ! ! * Added extra comments ! ! * Replaced GOTO statements with DO loops and IF statements ! - ! ---------------------------------------------------------------------------------------- - ! Future revisions: - ! - ! (none planned) - ! ! -------------------------------------------------------------------------------------------- - IMPLICIT NONE + implicit none ! Input - REAL(DP), DIMENSION(:), INTENT(IN) :: TOLD ! input time array - REAL(DP), DIMENSION(:), INTENT(IN) :: QOLD ! input flow array - REAL(DP), DIMENSION(:), INTENT(IN) :: TNEW ! desired output times + real(dp), dimension(:), intent(in) :: TOLD ! input time array + real(dp), dimension(:), intent(in) :: QOLD ! input flow array + real(dp), dimension(:), intent(in) :: TNEW ! desired output times ! Output - REAL(DP), DIMENSION(:), INTENT(OUT) :: QNEW ! flow averaged for desired times - INTEGER(I4B), INTENT(OUT) :: IERR ! error, 1= bad bounds + real(dp), dimension(:), intent(out) :: QNEW ! flow averaged for desired times + integer(i4b), intent(out) :: IERR ! error, 1= bad bounds character(*), intent(out) :: MESSAGE ! error message ! Internal - INTEGER(I4B) :: NOLD ! number of elements in input array - INTEGER(I4B) :: NNEW ! number of desired new times - INTEGER(I4B) :: IOLDLOOP ! loop through input times - INTEGER(I4B) :: INEWLOOP ! loop through desired times - REAL(DP) :: T0,T1 ! time at start/end of the time step - INTEGER(I4B) :: IBEG ! identify input times spanning T0 - INTEGER(I4B) :: IEND ! identify input times spanning T1 - INTEGER(I4B) :: IMID ! input times in middle of the curve - REAL(DP) :: AREAB ! area at the start of the time step - REAL(DP) :: AREAE ! area at the end of the time step - REAL(DP) :: AREAM ! area at the middle of the time step - REAL(DP) :: AREAS ! sum of all areas - REAL(DP) :: SLOPE ! slope between two input data values - REAL(DP) :: QEST0 ! flow estimate at point T0 - REAL(DP) :: QEST1 ! flow estimate at point T1 - ! -------------------------------------------------------------------------------------------- - IERR=0; message='INTERP_RCH/' + integer(i4b) :: NOLD ! number of elements in input array + integer(i4b) :: NNEW ! number of desired new times + integer(i4b) :: IOLDLOOP ! loop through input times + integer(i4b) :: INEWLOOP ! loop through desired times + real(dp) :: T0,T1 ! time at start/end of the time step + integer(i4b) :: IBEG ! identify input times spanning T0 + integer(i4b) :: IEND ! identify input times spanning T1 + integer(i4b) :: IMID ! input times in middle of the curve + real(dp) :: AREAB ! area at the start of the time step + real(dp) :: AREAE ! area at the end of the time step + real(dp) :: AREAM ! area at the middle of the time step + real(dp) :: AREAS ! sum of all areas + real(dp) :: SLOPE ! slope between two input data values + real(dp) :: QEST0 ! flow estimate at point T0 + real(dp) :: QEST1 ! flow estimate at point T1 + + IERR=0; message='interp_rch/' ! get array size - NOLD = SIZE(TOLD); NNEW = SIZE(TNEW) + NOLD = size(TOLD); NNEW = size(TNEW) ! check that the input time series starts before the first required output time ! and ends after the last required output time - IF( (TOLD(1).GT.TNEW(1)) .OR. (TOLD(NOLD).LT.TNEW(NNEW)) ) THEN - IERR=1; message=trim(message)//'bad bounds'; RETURN - ENDIF + if( (TOLD(1).GT.TNEW(1)) .OR. (TOLD(NOLD).LT.TNEW(NNEW)) ) then + IERR=1; message=trim(message)//'bad bounds'; return + end if ! loop through the output times - DO INEWLOOP=2,NNEW + do INEWLOOP=2,NNEW T0 = TNEW(INEWLOOP-1) ! start of the time step T1 = TNEW(INEWLOOP) ! end of the time step IBEG=1 ! identify the index values that span the start of the time step - BEG_ID: DO IOLDLOOP=2,NOLD - IF(T0.LE.TOLD(IOLDLOOP)) THEN + BEG_ID: do IOLDLOOP=2,NOLD + if (T0.LE.TOLD(IOLDLOOP)) then IBEG = IOLDLOOP - EXIT - ENDIF - END DO BEG_ID + exit + end if + end do BEG_ID IEND=1 ! identify the index values that span the end of the time step - END_ID: DO IOLDLOOP=1,NOLD - IF(T1.LE.TOLD(IOLDLOOP)) THEN + END_ID: do IOLDLOOP=1,NOLD + if (T1.LE.TOLD(IOLDLOOP)) then IEND = IOLDLOOP - EXIT - ENDIF - END DO END_ID + exit + end if + end do END_ID ! initialize the areas AREAB=0D0; AREAE=0D0; AREAM=0D0 ! special case: both TNEW(INEWLOOP-1) and TNEW(INEWLOOP) are within two original values ! (implies IBEG=IEND) -- estimate values at both end-points and average - IF(T1.LT.TOLD(IBEG)) THEN + if (T1.LT.TOLD(IBEG)) then SLOPE = (QOLD(IBEG)-QOLD(IBEG-1))/(TOLD(IBEG)-TOLD(IBEG-1)) QEST0 = SLOPE*(T0-TOLD(IBEG-1)) + QOLD(IBEG-1) QEST1 = SLOPE*(T1-TOLD(IBEG-1)) + QOLD(IBEG-1) QNEW(INEWLOOP-1) = 0.5*(QEST0 + QEST1) - CYCLE ! loop back to the next desired time - ENDIF + cycle ! loop back to the next desired time + end if ! estimate the area under the curve at the start of the time step - IF(T0.LT.TOLD(IBEG)) THEN ! if equal process as AREAM + if (T0.LT.TOLD(IBEG)) then ! if equal process as AREAM SLOPE = (QOLD(IBEG)-QOLD(IBEG-1))/(TOLD(IBEG)-TOLD(IBEG-1)) QEST0 = SLOPE*(T0-TOLD(IBEG-1)) + QOLD(IBEG-1) AREAB = (TOLD(IBEG)-T0) * 0.5*(QEST0 + QOLD(IBEG)) - ENDIF + end if ! estimate the area under the curve at the end of the time step - IF(T1.LT.TOLD(IEND)) THEN ! if equal process as AREAM + if (T1.LT.TOLD(IEND)) then ! if equal process as AREAM SLOPE = (QOLD(IEND)-QOLD(IEND-1))/(TOLD(IEND)-TOLD(IEND-1)) QEST1 = SLOPE*(T1-TOLD(IEND-1)) + QOLD(IEND-1) AREAE = (T1-TOLD(IEND-1)) * 0.5*(QOLD(IEND-1) + QEST1) - ENDIF + endif ! check if there are extra points to process - IF(IBEG.LT.IEND) THEN + if (IBEG.LT.IEND) then ! loop through remaining points - DO IMID=IBEG+1,IEND - IF(IMID.LT.IEND .OR. & + do IMID=IBEG+1,IEND + if (IMID.LT.IEND .or. & ! process the end slice as AREAM, but only if not already AREAB - (IMID.EQ.IEND.AND.T1.EQ.TOLD(IEND).AND.T0.LT.TOLD(IEND-1)) ) THEN + (IMID.EQ.IEND.and.T1.EQ.TOLD(IEND).and.T0.LT.TOLD(IEND-1)) ) then ! compute AREAM AREAM = AREAM + (TOLD(IMID) - TOLD(IMID-1)) * 0.5*(QOLD(IMID-1) + QOLD(IMID)) - ENDIF ! if point is valid - END DO ! IMID - ENDIF ! If there is a possibility that middle points even exist + end if ! if point is valid + end do ! IMID + end if ! If there is a possibility that middle points even exist ! compute time step average AREAS = AREAB + AREAE + AREAM ! sum of all areas QNEW(INEWLOOP-1) = AREAS / (T1-T0) ! T1-T0 is the sum of all time slices - END DO + end do - end subroutine INTERP_RCH + END SUBROUTINE interp_rch -end module kwt_route_module +END MODULE kwt_route_module diff --git a/route/build/src/main_route.f90 b/route/build/src/main_route.f90 index ce100106..ee8a3b25 100644 --- a/route/build/src/main_route.f90 +++ b/route/build/src/main_route.f90 @@ -151,6 +151,7 @@ subroutine main_route(iens, & ! input: ensemble index river_basin, & ! input: river basin data type ixPrint, & ! input: index of the desired reach NETOPO, & ! input: reach topology data structure + RPARAM, & ! input: reach parameter data structure RCHFLX, & ! inout: reach flux data structure ierr,cmessage) ! output: error control if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif diff --git a/route/build/src/model_finalize.f90 b/route/build/src/model_finalize.f90 new file mode 100644 index 00000000..8fa3c7ab --- /dev/null +++ b/route/build/src/model_finalize.f90 @@ -0,0 +1,42 @@ +MODULE model_finalize + +USE nrtype, ONLY: i4b +USE public_var, ONLY: iulog ! i/o logical unit number + +implicit none + +private + +public :: finalize +public :: handle_err + +CONTAINS + + ! ********************************************************************* + ! public subroutine: finalize model + ! ********************************************************************* + SUBROUTINE finalize() + implicit none + write(iulog,'(a)') new_line('a'), '--------------------' + write(iulog,'(a)') 'Finished simulation' + write(iulog,'(a)') '--------------------' + stop + END SUBROUTINE finalize + + + ! ********************************************************************* + ! public subroutine: error handling + ! ********************************************************************* + SUBROUTINE handle_err(err,message) + implicit none + integer(i4b),intent(in)::err ! error code + character(*),intent(in)::message ! error message + if(err/=0)then + write(iulog,*) 'FATAL ERROR: '//trim(message) + call flush(6) + stop + endif + END SUBROUTINE handle_err + + +END MODULE model_finalize diff --git a/route/build/src/model_setup.f90 b/route/build/src/model_setup.f90 index 3a517402..c8f5a4fc 100644 --- a/route/build/src/model_setup.f90 +++ b/route/build/src/model_setup.f90 @@ -1,7 +1,7 @@ MODULE model_setup ! data types -USE nrtype, ONLY : i4b,dp,lgt ! variable types, etc. +USE nrtype, ONLY : i4b,i8b,dp,lgt ! variable types, etc. USE nrtype, ONLY : strLen ! length of characters USE dataTypes, ONLY : var_ilength ! integer type: var(:)%dat USE dataTypes, ONLY : var_clength ! integer type: var(:)%dat @@ -16,6 +16,7 @@ MODULE model_setup USE io_netcdf, ONLY : close_nc ! close netcdf +USE nr_utility_module, ONLY : findIndex ! get array index of matching element USE nr_utility_module, ONLY : unique ! get unique element array USE nr_utility_module, ONLY : indexx ! get rank of data value @@ -88,6 +89,7 @@ subroutine init_data(ierr, message) USE public_var, ONLY : is_remap ! logical whether or not runnoff needs to be mapped to river network HRU USE public_var, ONLY : ntopAugmentMode ! River network augmentation mode USE public_var, ONLY : idSegOut ! outlet segment ID (-9999 => no outlet segment specified) + USE public_var, ONLY : desireId ! ID of reach to be checked by on-screen printing USE var_lookup, ONLY : ixHRU2SEG ! index of variables for data structure USE var_lookup, ONLY : ixNTOPO ! index of variables for data structure USE globalData, ONLY : RCHFLX ! Reach flux data structures (entire river network) @@ -97,6 +99,7 @@ subroutine init_data(ierr, message) USE globalData, ONLY : nEns ! number of ensembles USE globalData, ONLY : basinID ! HRU id vector USE globalData, ONLY : reachID ! reach ID vector + USE globalData, ONLY : ixPrint ! reach index to be examined by on-screen printing USE globalData, ONLY : runoff_data ! runoff data structure USE globalData, ONLY : remap_data ! runoff mapping data structure @@ -148,6 +151,9 @@ subroutine init_data(ierr, message) reachID(iRch) = structNTOPO(iRch)%var(ixNTOPO%segId)%dat(1) end do + ! get reach index to be examined by on-screen printing + if (desireId/=integerMissing) ixPrint = findIndex(reachID, desireId, integerMissing) + ! runoff and remap data initialization (TO DO: split runoff and remap initialization) call init_runoff(is_remap, & ! input: logical whether or not runnoff needs to be mapped to river network HRU remap_data, & ! output: data structure to remap data @@ -172,12 +178,12 @@ end subroutine init_data SUBROUTINE update_time(finished, ierr, message) USE public_var, ONLY : dt + USE public_var, ONLY : calendar USE globalData, ONLY : TSEC ! beginning/ending of simulation time step [sec] USE globalData, ONLY : iTime ! time index at simulation time step - USE globalData, ONLY : roJulday ! julian day: runoff input time - USE globalData, ONLY : modJulday ! julian day: at model time step - USE globalData, ONLY : endJulday ! julian day: at end of simulation USE globalData, ONLY : simout_nc ! netCDF meta data + USE globalData, ONLY : endCal ! model ending datetime + USE globalData, ONLY : modTime ! model datetime implicit none ! output @@ -190,17 +196,13 @@ SUBROUTINE update_time(finished, ierr, message) ! initialize error control ierr=0; message='update_time/' - if (abs(modJulday-endJulday) both, 1->IRF, 2->KWT, otherwise error USE public_var, ONLY : fname_state_in ! name of state input file - USE public_var, ONLY : output_dir ! directory containing output data + USE public_var, ONLY : restart_dir ! directory containing output data USE globalData, ONLY : RCHFLX ! reach flux structure USE globalData, ONLY : TSEC ! begining/ending of simulation time step [sec] @@ -246,7 +249,7 @@ subroutine init_state(ierr, message) ! read restart file and initialize states if (trim(fname_state_in)/=charMissing) then - call read_state_nc(trim(output_dir)//trim(fname_state_in), routOpt, T0, T1, ierr, cmessage) + call read_state_nc(trim(restart_dir)//trim(fname_state_in), routOpt, T0, T1, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif TSEC(0)=T0; TSEC(1)=T1 @@ -258,6 +261,8 @@ subroutine init_state(ierr, message) RCHFLX(:,:)%BASIN_QI = 0._dp RCHFLX(:,:)%BASIN_QR(0) = 0._dp RCHFLX(:,:)%BASIN_QR(1) = 0._dp + RCHFLX(:,:)%REACH_VOL(0) = 0._dp + RCHFLX(:,:)%REACH_VOL(1) = 0._dp ! initialize time TSEC(0)=0._dp; TSEC(1)=dt @@ -273,11 +278,11 @@ SUBROUTINE init_time(nTime, & ! input: number of time steps ierr, message) ! output ! subroutines: - USE process_time_module, ONLY : process_time ! process time information - USE process_time_module, ONLY : process_calday! compute data and time from julian day - USE io_netcdf, ONLY : get_nc ! netcdf input + USE io_netcdf, ONLY : open_nc ! netcdf input + USE io_netcdf, ONLY : close_nc ! netcdf input + USE io_netcdf, ONLY : get_nc ! netcdf input ! derived datatype - USE dataTypes, ONLY : time ! time data type + USE date_time, ONLY : datetime ! time data type ! public data USE public_var, ONLY : input_dir ! directory containing input data USE public_var, ONLY : fname_qsim ! simulated runoff netCDF name @@ -286,18 +291,20 @@ SUBROUTINE init_time(nTime, & ! input: number of time steps USE public_var, ONLY : simStart ! date string defining the start of the simulation USE public_var, ONLY : simEnd ! date string defining the end of the simulation USE public_var, ONLY : calendar ! calendar name + USE public_var, ONLY : dt + USE public_var, ONLY : secprday USE public_var, ONLY : restart_write ! restart write option USE public_var, ONLY : restart_date ! restart date + USE public_var, ONLY : restart_month ! + USE public_var, ONLY : restart_day ! + USE public_var, ONLY : restart_hour ! ! saved time variables USE globalData, ONLY : timeVar ! time variables (unit given by runoff data) USE globalData, ONLY : iTime ! time index at runoff input time step - USE globalData, ONLY : refJulday ! julian day: reference - USE globalData, ONLY : roJulday ! julian day: runoff input time - USE globalData, ONLY : startJulday ! julian day: start of routing simulation - USE globalData, ONLY : endJulday ! julian day: end of routing simulation - USE globalData, ONLY : modJulday ! julian day: at model time step - USE globalData, ONLY : restartJulday ! julian day: at restart USE globalData, ONLY : modTime ! model time data (yyyy:mm:dd:hh:mm:sec) + USE globalData, ONLY : endCal ! simulation end time data (yyyy:mm:dd:hh:mm:sec) + USE globalData, ONLY : restCal ! restart time data (yyyy:mm:dd:hh:mm:sec) + USE globalData, ONLY : dropCal ! restart dropoff calendar date/time implicit none @@ -307,10 +314,15 @@ SUBROUTINE init_time(nTime, & ! input: number of time steps integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! local variable + integer(i4b) :: ncidRunoff integer(i4b) :: ix - type(time) :: rofCal - type(time) :: simCal - real(dp) :: convTime2Days + type(datetime) :: refCal + type(datetime) :: roCal(nTime) + type(datetime) :: startCal + type(datetime) :: dummyCal + integer(i4b) :: nDays ! number of days in a month + real(dp) :: convTime2sec + real(dp) :: sec(nTime) character(len=7) :: t_unit character(len=strLen) :: cmessage ! error message of downwind routine character(len=50) :: fmt1='(a,I4,a,I2.2,a,I2.2,x,I2.2,a,I2.2,a,F5.2)' @@ -319,102 +331,120 @@ SUBROUTINE init_time(nTime, & ! input: number of time steps ierr=0; message='init_time/' ! time initialization - allocate(timeVar(nTime), roJulday(nTime), stat=ierr) + allocate(timeVar(nTime), stat=ierr) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif ! get the time data - call get_nc(trim(input_dir)//trim(fname_qsim), vname_time, timeVar, 1, nTime, ierr, cmessage) + call open_nc(trim(input_dir)//trim(fname_qsim), 'r', ncidRunoff, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + call get_nc(ncidRunoff, vname_time, timeVar, 1, nTime, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + call close_nc(ncidRunoff, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif ! get the time multiplier needed to convert time to units of days t_unit = trim( time_units(1:index(time_units,' ')) ) select case( trim(t_unit) ) - case('seconds','second','sec','s'); convTime2Days=86400._dp - case('minutes','minute','min'); convTime2Days=1440._dp - case('hours','hour','hr','h'); convTime2Days=24._dp - case('days','day','d'); convTime2Days=1._dp + case('seconds','second','sec','s'); convTime2sec=1._dp + case('minutes','minute','min'); convTime2sec=60._dp + case('hours','hour','hr','h'); convTime2sec=3600._dp + case('days','day','d'); convTime2sec=86400._dp case default ierr=20; message=trim(message)//'= '//trim(t_unit)//': must be seconds, minutes, hours or days.'; return end select - ! extract time information from the control information - call process_time(time_units, calendar, refJulday, ierr, cmessage) - if(ierr/=0) then; message=trim(message)//trim(cmessage)//' [refJulday]'; return; endif - call process_time(trim(simStart),calendar, startJulday, ierr, cmessage) - if(ierr/=0) then; message=trim(message)//trim(cmessage)//' [startJulday]'; return; endif - call process_time(trim(simEnd), calendar, endJulday, ierr, cmessage) - if(ierr/=0) then; message=trim(message)//trim(cmessage)//' [endJulday]'; return; endif + ! extract datetime from the control information + call refCal%str2datetime(time_units, ierr, cmessage) + if(ierr/=0) then; message=trim(message)//trim(cmessage)//' [refCal]'; return; endif + call startCal%str2datetime(simStart, ierr, cmessage) + if(ierr/=0) then; message=trim(message)//trim(cmessage)//' [startCal]'; return; endif + call endCal%str2datetime(simEnd, ierr, cmessage) + if(ierr/=0) then; message=trim(message)//trim(cmessage)//' [endCal]'; return; endif - ! Julian day at first time step in runoff data - ! convert time unit in runoff netCDF to day + ! calendar in runoff data + sec(:) = timeVar(:)*convTime2sec do ix = 1, nTime - roJulday(ix) = refJulday + timeVar(ix)/convTime2Days + roCal(ix) = refCal%add_sec(sec(ix), calendar, ierr, cmessage) end do ! check that the dates are aligned - if(endJulday= ', trim(simStart), new_line('a'), '= ', trim(simEnd) ierr=20; message=trim(message)//trim(cmessage); return endif ! check sim_start is before the last time step in runoff data - if(startJulday>roJulday(nTime)) then - call process_calday(roJulday(nTime), calendar, rofCal, ierr, cmessage) - call process_calday(startJulday, calendar, simCal, ierr, cmessage) + if(startCal > roCal(nTime)) then write(iulog,'(2a)') new_line('a'),'ERROR: is after the first time step in input runoff' - write(iulog,fmt1) ' runoff_end : ', rofCal%iy,'-',rofCal%im,'-',rofCal%id, rofCal%ih,':', rofCal%imin,':',rofCal%dsec - write(iulog,fmt1) ' : ', simCal%iy,'-',simCal%im,'-',simCal%id, simCal%ih,':', simCal%imin,':',simCal%dsec + write(iulog,fmt1) ' runoff_end : ', roCal(nTime)%year(),'-',roCal(nTime)%month(),'-',roCal(nTime)%day(),roCal(nTime)%hour(),':',roCal(nTime)%minute(),':',roCal(nTime)%sec() + write(iulog,fmt1) ' : ', startCal%year(),'-',startCal%month(),'-',startCal%day(), startCal%hour(),':', startCal%minute(),':',startCal%sec() ierr=20; message=trim(message)//'check against runoff input time'; return endif ! Compare sim_start vs. time at first time step in runoff data - if (startJulday < roJulday(1)) then - call process_calday(roJulday(1), calendar, rofCal, ierr, cmessage) - call process_calday(startJulday, calendar, simCal, ierr, cmessage) + if (startCal < roCal(1)) then write(iulog,'(2a)') new_line('a'),'WARNING: is before the first time step in input runoff' - write(iulog,fmt1) ' runoff_start: ', rofCal%iy,'-',rofCal%im,'-',rofCal%id, rofCal%ih,':', rofCal%imin,':',rofCal%dsec - write(iulog,fmt1) ' : ', simCal%iy,'-',simCal%im,'-',simCal%id, simCal%ih,':', simCal%imin,':',simCal%dsec + write(iulog,fmt1) ' runoff_start : ', roCal(1)%year(),'-',roCal(1)%month(),'-',roCal(1)%day(), roCal(1)%hour(),':', roCal(1)%minute(),':',roCal(1)%sec() + write(iulog,fmt1) ' : ', startCal%year(),'-',startCal%month(),'-',startCal%day(), startCal%hour(),':', startCal%minute(),':',startCal%sec() write(iulog,'(a)') ' Reset to runoff_start' - startJulday = roJulday(1) + startCal = roCal(1) endif ! Compare sim_end vs. time at last time step in runoff data - if (endJulday > roJulday(nTime)) then - call process_calday(roJulday(nTime), calendar, rofCal, ierr, cmessage) - call process_calday(endJulday, calendar, simCal, ierr, cmessage) + if (endCal > roCal(nTime)) then write(iulog,'(2a)') new_line('a'),'WARNING: is after the last time step in input runoff' - write(iulog,fmt1) ' runoff_end: ', rofCal%iy,'-',rofCal%im,'-',rofCal%id, rofCal%ih,':', rofCal%imin,':',rofCal%dsec - write(iulog,fmt1) ' : ', simCal%iy,'-',simCal%im,'-',simCal%id, simCal%ih,':', simCal%imin,':',simCal%dsec + write(iulog,fmt1) ' runoff_end : ', roCal(nTime)%year(),'-',roCal(nTime)%month(),'-',roCal(nTime)%day(),roCal(nTime)%hour(),':',roCal(nTime)%minute(),':',roCal(nTime)%sec() + write(iulog,fmt1) ' : ', endCal%year(),'-',endCal%month(),'-',endCal%day(), endCal%hour(),':', endCal%minute(),':',endCal%sec() write(iulog,'(a)') ' Reset to runoff_end' - endJulday = roJulday(nTime) + endCal = roCal(nTime) endif - ! fast forward time to time index at simStart and save iTime and modJulday + ! fast forward time to time index at simStart and save iTime and modTime(1) do ix = 1, nTime - modJulday = roJulday(ix) - if( modJulday < startJulday ) cycle + modTime(1) = roCal(ix) + if( modTime(1) < startCal ) cycle exit enddo iTime = ix - ! initialize previous model time - modTime(0) = time(integerMissing, integerMissing, integerMissing, integerMissing, integerMissing, realMissing) + ! Set restart calendar date/time and dropoff calendar date/time and + ! -- For periodic restart options --------------------------------------------------------------------- + ! Ensure that user-input restart month, day are valid. + ! "Annual" option: if user input day exceed number of days given user input month, set to last day + ! "Monthly" option: use 2000-01 as template calendar yr/month + ! "Daily" option: use 2000-01-01 as template calendar yr/month/day + select case(trim(restart_write)) + case('Annual','annual') + call dummyCal%set_datetime(2000, restart_month, 1, 0, 0, 0.0_dp) + nDays = dummyCal%ndays_month(calendar, ierr, cmessage) + if(ierr/=0) then; message=trim(message)//trim(cmessage); return; endif + if (restart_day > nDays) restart_day=nDays + case('Monthly','monthly'); restart_month = 1 + case('Daily','daily'); restart_month = 1; restart_day = 1 + end select - ! restart drop off time select case(trim(restart_write)) case('last','Last') - call process_time(trim(simEnd), calendar, restartJulday, ierr, cmessage) - if(ierr/=0) then; message=trim(message)//trim(cmessage)//' [restartDate]'; return; endif - case('never','Never') - restartJulday = 0.0_dp + dropCal = endCal + restart_month = dropCal%month(); restart_day = dropCal%day(); restart_hour = dropCal%hour() case('specified','Specified') if (trim(restart_date) == charMissing) then ierr=20; message=trim(message)//' must be provided when option is "specified"'; return end if - call process_time(trim(restart_date),calendar, restartJulday, ierr, cmessage) - if(ierr/=0) then; message=trim(message)//trim(cmessage)//' [restartDate]'; return; endif + call restCal%str2datetime(restart_date, ierr, cmessage) + if(ierr/=0) then; message=trim(message)//trim(cmessage)//' [restart_date]'; return; endif + dropCal = restCal%add_sec(-dt, calendar, ierr, cmessage) + if(ierr/=0) then; message=trim(message)//trim(cmessage)//' [restCal->dropCal]'; return; endif + restart_month = dropCal%month(); restart_day = dropCal%day(); restart_hour = dropCal%hour() + case('Annual','Monthly','Daily','annual','monthly','daily') + call restCal%set_datetime(2000, restart_month, restart_day, restart_hour, 0, 0._dp) + dropCal = restCal%add_sec(-dt, calendar, ierr, cmessage) + if(ierr/=0) then; message=trim(message)//trim(cmessage)//' [ dropCal for periodical restart]'; return; endif + restart_month = dropCal%month(); restart_day = dropCal%day(); restart_hour = dropCal%hour() + case('never','Never') + call dropCal%set_datetime(integerMissing, integerMissing, integerMissing, integerMissing, integerMissing, realMissing) case default - ierr=20; message=trim(message)//'Current accepted options: last[Last], never[Never], specified[Specified]'; return + ierr=20; message=trim(message)//'Current accepted options: L[l]ast, N[n]ever, S[s]pecified, A[a]nnual, M[m]onthly, D[d]aily'; return end select END SUBROUTINE init_time @@ -616,7 +646,7 @@ SUBROUTINE init_runoff(remap_flag, & ! input: logical whether or not runno integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! local variables - integer(i4b), allocatable :: unq_qhru_id(:) + integer(i8b), allocatable :: unq_qhru_id(:) integer(i4b), allocatable :: unq_idx(:) character(len=strLen) :: cmessage ! error message from subroutine @@ -711,8 +741,8 @@ SUBROUTINE get_qix(qid,qidMaster,qix,ierr,message) implicit none ! input - integer(i4b), intent(in) :: qid(:) ! ID of input vector - integer(i4b), intent(in) :: qidMaster(:) ! ID of master vector + integer(i8b), intent(in) :: qid(:) ! ID of input vector + integer(i8b), intent(in) :: qidMaster(:) ! ID of master vector ! output integer(i4b), intent(out) :: qix(:) ! index within master vector integer(i4b), intent(out) :: ierr ! error code @@ -783,4 +813,3 @@ SUBROUTINE get_qix(qid,qidMaster,qix,ierr,message) END SUBROUTINE get_qix END MODULE model_setup - diff --git a/route/build/src/ncio_utils.f90 b/route/build/src/ncio_utils.f90 index f81af79c..ef86d603 100644 --- a/route/build/src/ncio_utils.f90 +++ b/route/build/src/ncio_utils.f90 @@ -1,4 +1,4 @@ -module io_netcdf +MODULE io_netcdf USE nrtype USE netcdf @@ -25,6 +25,7 @@ module io_netcdf module procedure get_iscalar module procedure get_dscalar module procedure get_ivec + module procedure get_ivec_long module procedure get_dvec module procedure get_2d_iarray module procedure get_2d_darray @@ -176,83 +177,67 @@ end subroutine get_var_dims ! ********************************************************************* ! subroutine: get vector dimension from netCDF ! ********************************************************************* - subroutine get_nc_dim_len(fname, & ! input: filename + subroutine get_nc_dim_len(ncid, & ! input: netcdf ID dname, & ! input: variable name nDim, & ! output: Size of dimension ierr, message) ! output: error control implicit none ! input variables - character(*), intent(in) :: fname ! filename + integer(i4b), intent(in) :: ncid ! NetCDF file ID character(*), intent(in) :: dname ! dimension name ! output variables integer(i4b), intent(out) :: nDim ! size of dimension integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! local variables - integer(i4b) :: ncid ! NetCDF file ID integer(i4b) :: iDimID ! NetCDF dimension ID - ! initialize error control - ierr=0; message='get_nc_dim_len/' - ! open NetCDF file - ierr = nf90_open(trim(fname),nf90_nowrite,ncid) - if(ierr/=0)then; message=trim(message)//'['//trim(nf90_strerror(ierr))//'; file='//trim(fname)//']'; return; endif + ierr=0; message='get_nc_dim_len/' ! get the ID of the dimension - ierr = nf90_inq_dimid(ncid, dname, iDimID) + ierr = nf90_inq_dimid(ncid, trim(dname), iDimID) if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr))//'; name='//trim(dname); return; endif ! get the length of the dimension ierr = nf90_inquire_dimension(ncid, iDimID, len=nDim) if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif - ! close output file - ierr = nf90_close(ncid) - if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif - - end subroutine + end subroutine get_nc_dim_len ! ********************************************************************* ! subroutine: get attribute values for a variable ! ********************************************************************* - FUNCTION check_attr(fname, vname, attr_name) + FUNCTION check_attr(ncid, vname, attr_name) implicit none ! input - character(*), intent(in) :: fname ! filename + integer(i4b), intent(in) :: ncid ! NetCDF file ID character(*), intent(in) :: vname ! variable name character(*), intent(in) :: attr_name ! attribute name logical(lgt) :: check_attr ! local integer(i4b) :: ierr ! error code - integer(i4b) :: ncid ! NetCDF file ID integer(i4b) :: iVarID ! variable ID - ! open file for reading - ierr = nf90_open(fname, nf90_nowrite, ncid) - ! get the ID of the variable ierr = nf90_inq_varid(ncid, trim(vname), iVarID) ierr = nf90_inquire_attribute(ncid, iVarID, attr_name) check_attr = (ierr == nf90_noerr) - ! close output file - ierr = nf90_close(ncid) - END FUNCTION check_attr ! ********************************************************************* ! subroutine: get attribute values for a variable ! ********************************************************************* - subroutine get_var_attr_char(fname, & ! input: filename + subroutine get_var_attr_char(ncid, & ! input: netcdf id vname, & ! input: variable name attr_name, & ! inpu: attribute name attr_value, & ! output: attribute value ierr, message) ! output: error control implicit none ! input variables - character(*), intent(in) :: fname ! filename + integer(i4b), intent(in) :: ncid ! NetCDF file ID character(*), intent(in) :: vname ! variable name character(*), intent(in) :: attr_name ! attribute name ! output variables @@ -261,14 +246,9 @@ subroutine get_var_attr_char(fname, & ! input: filename character(*), intent(out) :: message ! error message ! local variables integer(i4b) :: var_type ! attribute variable type - integer(i4b) :: ncid ! NetCDF file ID integer(i4b) :: iVarID ! variable ID - ! initialize error control - ierr=0; message='get_var_attr_char/' - ! open file for reading - ierr = nf90_open(fname, nf90_nowrite, ncid) - if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr))//'; name='//trim(fname); return; endif + ierr=0; message='get_var_attr_char/' ! get the ID of the variable ierr = nf90_inq_varid(ncid, trim(vname), ivarID) @@ -276,7 +256,7 @@ subroutine get_var_attr_char(fname, & ! input: filename ! Inquire attribute type, NF90_CHAR(=2) ierr = nf90_inquire_attribute(ncid, ivarID, attr_name, xtype=var_type) - if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr))//'; nc='//trim(fname)//'; attr='//trim(attr_name); return; endif + if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr))//'; attr='//trim(attr_name); return; endif if (var_type /= nf90_char)then; ierr=20; message=trim(message)//'attribute type must be character'; return; endif @@ -284,23 +264,19 @@ subroutine get_var_attr_char(fname, & ! input: filename ierr = nf90_get_att(ncid, ivarID, attr_name, attr_value) if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif - ! close the NetCDF file - ierr = nf90_close(ncid) - if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif - end subroutine get_var_attr_char ! ********************************************************************* ! subroutine: get attribute values for a real variable ! ********************************************************************* - subroutine get_var_attr_real(fname, & ! input: filename + subroutine get_var_attr_real(ncid, & ! input: netcdf id vname, & ! input: variable name attr_name, & ! inpu: attribute name attr_value, & ! output: attribute value in real ierr, message) ! output: error control implicit none ! input variables - character(*), intent(in) :: fname ! filename + integer(i4b), intent(in) :: ncid ! NetCDF file ID character(*), intent(in) :: vname ! variable name character(*), intent(in) :: attr_name ! attribute name ! output variables @@ -309,14 +285,9 @@ subroutine get_var_attr_real(fname, & ! input: filename character(*), intent(out) :: message ! error message ! local variables integer(i4b) :: var_type ! attribute variable type - integer(i4b) :: ncid ! NetCDF file ID integer(i4b) :: iVarID ! variable ID - ! initialize error control - ierr=0; message='get_var_attr_real/' - ! open file for reading - ierr = nf90_open(fname, nf90_nowrite, ncid) - if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr))//'; name='//trim(fname); return; endif + ierr=0; message='get_var_attr_real/' ! get the ID of the variable ierr = nf90_inq_varid(ncid, trim(vname), ivarID) @@ -324,7 +295,7 @@ subroutine get_var_attr_real(fname, & ! input: filename ! Inquire attribute type, NF90_CHAR(=2) ierr = nf90_inquire_attribute(ncid, ivarID, attr_name, xtype=var_type) - if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr))//'; nc='//trim(fname)//'; attr='//trim(attr_name); return; endif + if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr))//'; attr='//trim(attr_name); return; endif if (var_type /= nf90_float .and. var_type /= nf90_double)then; ierr=20; message=trim(message)//'attribute type must be real'; return; endif @@ -332,24 +303,20 @@ subroutine get_var_attr_real(fname, & ! input: filename ierr = nf90_get_att(ncid, ivarID, attr_name, attr_value) if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif - ! close the NetCDF file - ierr = nf90_close(ncid) - if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif - end subroutine get_var_attr_real ! ********************************************************************* ! subroutine: get integer scalar value from netCDF ! ********************************************************************* - subroutine get_iscalar(fname, & ! input: filename + subroutine get_iscalar(ncid, & ! input: netcdf id vname, & ! input: variable name array, & ! output: variable data iStart, & ! input: start index ierr, message) ! output: error control implicit none ! input variables - character(*), intent(in) :: fname ! filename + integer(i4b), intent(in) :: ncid ! NetCDF file ID character(*), intent(in) :: vname ! variable name integer(i4b), intent(in) :: iStart ! start index ! output variables @@ -360,10 +327,9 @@ subroutine get_iscalar(fname, & ! input: filename character(len=strLen) :: cmessage ! error message of downwind routine integer(i4b) :: array_vec(1) ! output variable data - ! initialize error control ierr=0; message='get_iscalar/' - call get_ivec(fname, vname, array_vec, iStart, 1, ierr, cmessage) + call get_ivec(ncid, vname, array_vec, iStart, 1, ierr, cmessage) array = array_vec(1) end subroutine @@ -371,14 +337,14 @@ subroutine get_iscalar(fname, & ! input: filename ! ********************************************************************* ! subroutine: double precision scalar value from netCDF ! ********************************************************************* - subroutine get_dscalar(fname, & ! input: filename + subroutine get_dscalar(ncid, & ! input: netcdf id vname, & ! input: variable name array, & ! output: variable data iStart, & ! input: start index ierr, message) ! output: error control implicit none ! input variables - character(*), intent(in) :: fname ! filename + integer(i4b), intent(in) :: ncid ! NetCDF file ID character(*), intent(in) :: vname ! variable name integer(i4b), intent(in) :: iStart ! start index ! output variables @@ -389,10 +355,9 @@ subroutine get_dscalar(fname, & ! input: filename character(len=strLen) :: cmessage ! error message of downwind routine real(dp) :: array_vec(1) ! output variable data - ! initialize error control ierr=0; message='get_dscalar/' - call get_dvec(fname, vname, array_vec, iStart, 1, ierr, cmessage) + call get_dvec(ncid, vname, array_vec, iStart, 1, ierr, cmessage) array = array_vec(1) end subroutine @@ -400,7 +365,7 @@ subroutine get_dscalar(fname, & ! input: filename ! ********************************************************************* ! subroutine: get integer vector value from netCDF ! ********************************************************************* - subroutine get_ivec(fname, & ! input: filename + subroutine get_ivec(ncid, & ! input: netcdf id vname, & ! input: variable name array, & ! output: variable data iStart, & ! input: start index @@ -408,7 +373,7 @@ subroutine get_ivec(fname, & ! input: filename ierr, message) ! output: error control implicit none ! input variables - character(*), intent(in) :: fname ! filename + integer(i4b), intent(in) :: ncid ! NetCDF file ID character(*), intent(in) :: vname ! variable name integer(i4b), intent(in) :: iStart ! start index integer(i4b), intent(in) :: iCount ! length of vector to be read in @@ -417,16 +382,10 @@ subroutine get_ivec(fname, & ! input: filename integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! local variables - integer(i4b) :: ncid ! NetCDF file ID integer(i4b) :: iVarID ! NetCDF variable ID - ! initialize error control ierr=0; message='get_ivec/' - ! open NetCDF file - ierr = nf90_open(trim(fname),nf90_nowrite,ncid) - if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif - ! get variable ID ierr = nf90_inq_varid(ncid,trim(vname),iVarId) if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif @@ -435,16 +394,46 @@ subroutine get_ivec(fname, & ! input: filename ierr = nf90_get_var(ncid, iVarID, array, start=(/iStart/), count=(/iCount/)) if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif - ! close output file - ierr = nf90_close(ncid) - if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif + end subroutine + + ! ********************************************************************* + ! subroutine: get integer vector value from netCDF + ! ********************************************************************* + subroutine get_ivec_long(ncid, & ! input: netcdf id + vname, & ! input: variable name + array, & ! output: variable data + iStart, & ! input: start index + iCount, & ! input: length of vector + ierr, message) ! output: error control + implicit none + ! input variables + integer(i4b), intent(in) :: ncid ! NetCDF file ID + character(*), intent(in) :: vname ! variable name + integer(i4b), intent(in) :: iStart ! start index + integer(i4b), intent(in) :: iCount ! length of vector to be read in + ! output variables + integer(i8b), intent(out) :: array(:) ! output variable data + integer(i4b), intent(out) :: ierr ! error code + character(*), intent(out) :: message ! error message + ! local variables + integer(i4b) :: iVarID ! NetCDF variable ID + + ierr=0; message='get_ivec_long/' + + ! get variable ID + ierr = nf90_inq_varid(ncid,trim(vname),iVarId) + if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif + + ! get the data + ierr = nf90_get_var(ncid, iVarID, array, start=(/iStart/), count=(/iCount/)) + if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif end subroutine ! ********************************************************************* ! subroutine: read a double precision vector ! ********************************************************************* - subroutine get_dvec(fname, & ! input: filename + subroutine get_dvec(ncid, & ! input: netcdf id vname, & ! input: variable name array, & ! output: variable data iStart, & ! input: start index @@ -452,7 +441,7 @@ subroutine get_dvec(fname, & ! input: filename ierr, message) ! output: error control implicit none ! input variables - character(*), intent(in) :: fname ! filename + integer(i4b), intent(in) :: ncid ! NetCDF file ID character(*), intent(in) :: vname ! variable name integer(i4b), intent(in) :: iStart ! start index integer(i4b), intent(in) :: iCount ! length of vector to be read in @@ -461,16 +450,10 @@ subroutine get_dvec(fname, & ! input: filename integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! local variables - integer(i4b) :: ncid ! NetCDF file ID integer(i4b) :: iVarID ! NetCDF variable ID - ! initialize error control ierr=0; message='get_dvec/' - ! open NetCDF file - ierr = nf90_open(trim(fname),nf90_nowrite,ncid) - if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif - ! get variable ID ierr = nf90_inq_varid(ncid,trim(vname),iVarId) if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif @@ -479,16 +462,12 @@ subroutine get_dvec(fname, & ! input: filename ierr = nf90_get_var(ncid, iVarID, array, start=(/iStart/), count=(/iCount/)) if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif - ! close output file - ierr = nf90_close(ncid) - if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif - end subroutine ! ********************************************************************* ! subroutine: read a integer 2D array ! ********************************************************************* - subroutine get_2d_iarray(fname, & ! input: filename + subroutine get_2d_iarray(ncid, & ! input: netcdf id vname, & ! input: variable name array, & ! output: variable data iStart, & ! input: start index @@ -496,7 +475,7 @@ subroutine get_2d_iarray(fname, & ! input: filename ierr, message) ! output: error control implicit none ! input variables - character(*), intent(in) :: fname ! filename + integer(i4b), intent(in) :: ncid ! NetCDF file ID character(*), intent(in) :: vname ! variable name integer(i4b), intent(in) :: iStart(1:2) ! start indices integer(i4b), intent(in) :: iCount(1:2) ! length of vector @@ -505,16 +484,10 @@ subroutine get_2d_iarray(fname, & ! input: filename integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! local variables - integer(i4b) :: ncid ! NetCDF file ID integer(i4b) :: iVarId ! NetCDF variable ID - ! initialize error control ierr=0; message='get_2d_iarray/' - ! open NetCDF file - ierr = nf90_open(trim(fname),nf90_nowrite,ncid) - if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif - ! get variable ID ierr = nf90_inq_varid(ncid,trim(vname),iVarId) if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif @@ -523,16 +496,12 @@ subroutine get_2d_iarray(fname, & ! input: filename ierr = nf90_get_var(ncid,iVarId,array,start=iStart,count=iCount) if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif - ! close output file - ierr = nf90_close(ncid) - if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif - end subroutine ! ********************************************************************* ! subroutine: read a integer 3D array ! ********************************************************************* - subroutine get_3d_iarray(fname, & ! input: filename + subroutine get_3d_iarray(ncid, & ! input: netcdf id vname, & ! input: variable name array, & ! output: variable data iStart, & ! input: start index @@ -540,7 +509,7 @@ subroutine get_3d_iarray(fname, & ! input: filename ierr, message) ! output: error control implicit none ! input variables - character(*), intent(in) :: fname ! filename + integer(i4b), intent(in) :: ncid ! NetCDF file ID character(*), intent(in) :: vname ! variable name integer(i4b), intent(in) :: iStart(1:3) ! start indices integer(i4b), intent(in) :: iCount(1:3) ! length of vector @@ -549,14 +518,9 @@ subroutine get_3d_iarray(fname, & ! input: filename integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! local variables - integer(i4b) :: ncid ! NetCDF file ID integer(i4b) :: iVarId ! NetCDF variable ID - ! initialize error control - ierr=0; message='get_3d_iarray/' - ! open NetCDF file - ierr = nf90_open(trim(fname),nf90_nowrite,ncid) - if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif + ierr=0; message='get_3d_iarray/' ! get variable ID ierr = nf90_inq_varid(ncid,trim(vname),iVarId) @@ -566,16 +530,12 @@ subroutine get_3d_iarray(fname, & ! input: filename ierr = nf90_get_var(ncid,iVarId,array,start=iStart,count=iCount) if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif - ! close output file - ierr = nf90_close(ncid) - if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif - end subroutine ! ********************************************************************* ! subroutine: read a integer 4D array ! ********************************************************************* - subroutine get_4d_iarray(fname, & ! input: filename + subroutine get_4d_iarray(ncid, & ! input: netcdf id vname, & ! input: variable name array, & ! output: variable data iStart, & ! input: start index @@ -583,7 +543,7 @@ subroutine get_4d_iarray(fname, & ! input: filename ierr, message) ! output: error control implicit none ! input variables - character(*), intent(in) :: fname ! filename + integer(i4b), intent(in) :: ncid ! NetCDF file ID character(*), intent(in) :: vname ! variable name integer(i4b), intent(in) :: iStart(1:4) ! start indices integer(i4b), intent(in) :: iCount(1:4) ! length of vector @@ -592,14 +552,9 @@ subroutine get_4d_iarray(fname, & ! input: filename integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! local variables - integer(i4b) :: ncid ! NetCDF file ID integer(i4b) :: iVarId ! NetCDF variable ID - ! initialize error control - ierr=0; message='get_4d_iarray/' - ! open NetCDF file - ierr = nf90_open(trim(fname),nf90_nowrite,ncid) - if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif + ierr=0; message='get_4d_iarray/' ! get variable ID ierr = nf90_inq_varid(ncid,trim(vname),iVarId) @@ -609,17 +564,13 @@ subroutine get_4d_iarray(fname, & ! input: filename ierr = nf90_get_var(ncid,iVarId,array,start=iStart,count=iCount) if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif - ! close output file - ierr = nf90_close(ncid) - if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif - end subroutine ! ********************************************************************* ! subroutine: read a double precision 2D array ! ********************************************************************* - subroutine get_2d_darray(fname, & ! input: filename + subroutine get_2d_darray(ncid, & ! input: netcdf id vname, & ! input: variable name array, & ! output: variable data iStart, & ! input: start index @@ -627,8 +578,8 @@ subroutine get_2d_darray(fname, & ! input: filename ierr, message) ! output: error control implicit none ! input variables - character(*), intent(in) :: fname ! filename - character(*), intent(in) :: vname ! variable name + integer(i4b), intent(in) :: ncid ! NetCDF file ID + character(*), intent(in) :: vname ! variable name integer(i4b), intent(in) :: iStart(1:2) ! start indices integer(i4b), intent(in) :: iCount(1:2) ! length of vector ! output variables @@ -636,14 +587,9 @@ subroutine get_2d_darray(fname, & ! input: filename integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! local variables - integer(i4b) :: ncid ! NetCDF file ID integer(i4b) :: iVarId ! NetCDF variable ID - ! initialize error control - ierr=0; message='get_2d_darray/' - ! open NetCDF file - ierr = nf90_open(trim(fname),nf90_nowrite,ncid) - if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif + ierr=0; message='get_2d_darray/' ! get variable ID ierr = nf90_inq_varid(ncid,trim(vname),iVarId) @@ -653,16 +599,12 @@ subroutine get_2d_darray(fname, & ! input: filename ierr = nf90_get_var(ncid,iVarId,array,start=iStart,count=iCount) if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif - ! close output file - ierr = nf90_close(ncid) - if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif - end subroutine ! ********************************************************************* ! subroutine: read a double precision 3D array ! ********************************************************************* - subroutine get_3d_darray(fname, & ! input: filename + subroutine get_3d_darray(ncid, & ! input: netcdf id vname, & ! input: variable name array, & ! output: variable data iStart, & ! input: start index @@ -670,7 +612,7 @@ subroutine get_3d_darray(fname, & ! input: filename ierr, message) ! output: error control implicit none ! input variables - character(*), intent(in) :: fname ! filename + integer(i4b), intent(in) :: ncid ! NetCDF file ID character(*), intent(in) :: vname ! variable name integer(i4b), intent(in) :: iStart(1:3) ! start indices integer(i4b), intent(in) :: iCount(1:3) ! length of vector @@ -679,14 +621,9 @@ subroutine get_3d_darray(fname, & ! input: filename integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! local variables - integer(i4b) :: ncid ! NetCDF file ID integer(i4b) :: iVarId ! NetCDF variable ID - ! initialize error control - ierr=0; message='get_3d_darray/' - ! open NetCDF file - ierr = nf90_open(trim(fname),nf90_nowrite,ncid) - if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif + ierr=0; message='get_3d_darray/' ! get variable ID ierr = nf90_inq_varid(ncid,trim(vname),iVarId) @@ -696,16 +633,12 @@ subroutine get_3d_darray(fname, & ! input: filename ierr = nf90_get_var(ncid,iVarId,array,start=iStart,count=iCount) if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif - ! close output file - ierr = nf90_close(ncid) - if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif - end subroutine ! ********************************************************************* ! subroutine: read a double precision 4D array ! ********************************************************************* - subroutine get_4d_darray(fname, & ! input: filename + subroutine get_4d_darray(ncid, & ! input: netcdf id vname, & ! input: variable name array, & ! output: variable data iStart, & ! input: start index @@ -713,7 +646,7 @@ subroutine get_4d_darray(fname, & ! input: filename ierr, message) ! output: error control implicit none ! input variables - character(*), intent(in) :: fname ! filename + integer(i4b), intent(in) :: ncid ! NetCDF file ID character(*), intent(in) :: vname ! variable name integer(i4b), intent(in) :: iStart(1:4) ! start indices integer(i4b), intent(in) :: iCount(1:4) ! length of vector @@ -722,14 +655,9 @@ subroutine get_4d_darray(fname, & ! input: filename integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! local variables - integer(i4b) :: ncid ! NetCDF file ID integer(i4b) :: iVarId ! NetCDF variable ID - ! initialize error control - ierr=0; message='get_4d_darray/' - ! open NetCDF file - ierr = nf90_open(trim(fname),nf90_nowrite,ncid) - if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif + ierr=0; message='get_4d_darray/' ! get variable ID ierr = nf90_inq_varid(ncid,trim(vname),iVarId) @@ -739,10 +667,6 @@ subroutine get_4d_darray(fname, & ! input: filename ierr = nf90_get_var(ncid,iVarId,array,start=iStart,count=iCount) if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif - ! close output file - ierr = nf90_close(ncid) - if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif - end subroutine ! ********************************************************************* @@ -1144,7 +1068,7 @@ SUBROUTINE close_nc(ncid, ierr, message) implicit none ! input - integer(i4b), intent(in) :: ncid ! Input: netcdf fine ID + integer(i4b), intent(in) :: ncid ! Input: netcdf fine ID ! output integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message diff --git a/route/build/src/network_topo.f90 b/route/build/src/network_topo.f90 index a1ad21d6..7d69bf46 100644 --- a/route/build/src/network_topo.f90 +++ b/route/build/src/network_topo.f90 @@ -783,9 +783,9 @@ END SUBROUTINE REACH_LIST ! ********************************************************************* ! new subroutine: identify all reaches above a given reach ! ********************************************************************* - SUBROUTINE REACH_MASK(& + SUBROUTINE reach_mask(& ! input - desireId, & ! input: reach index + outletId, & ! input: outlet reach id structNTOPO, & ! input: network topology structures structSeg, & ! input: Reach property structures nHRU, & ! input: number of HRUs @@ -810,7 +810,7 @@ SUBROUTINE REACH_MASK(& USE nr_utility_module, ONLY : arth ! Num. Recipies utilities IMPLICIT NONE ! input variables - integer(i4b) , intent(in) :: desireId ! id of the desired reach + integer(i4b) , intent(in) :: outletId ! id of the outlet reach type(var_ilength) , intent(inout) :: structNTOPO(:) ! network topology structure type(var_dlength) , intent(in) :: structSeg(:) ! stream segment properties integer(i4b) , intent(in) :: nHRU ! number of HRUs @@ -838,10 +838,10 @@ SUBROUTINE REACH_MASK(& integer(i4b) :: jxDesire ! index of desired reach ! ---------------------------------------------------------------------------------------- ! initialize error control - ierr=0; message='REACH_MASK/' + ierr=0; message='reach_mask/' ! check if we actually want the mask - if(desireId<0)then + if(outletId<0)then ! ---------- case 1: no mask desired --------------------------------------------------------------------------------------------------- @@ -866,7 +866,7 @@ SUBROUTINE REACH_MASK(& do iRch = 1,nRch idSeg_vec(iRch) = structNTOPO(iRch)%var(ixNTOPO%segId)%dat(1) end do - ixDesire = findIndex(idSeg_vec,desireId,integerMissing) + ixDesire = findIndex(idSeg_vec,outletId,integerMissing) if(ixDesire==integerMissing)then message=trim(message)//'unable to find index of desired reach id' ierr=20; return diff --git a/route/build/src/nr_utility.f90 b/route/build/src/nr_utility.f90 index f662e4ea..5f4987bc 100644 --- a/route/build/src/nr_utility.f90 +++ b/route/build/src/nr_utility.f90 @@ -3,8 +3,24 @@ module nr_utility_module ! contains functions that should really be part of the fortran standard, but are not implicit none INTERFACE arth - MODULE PROCEDURE arth_r, arth_d, arth_i + MODULE PROCEDURE arth_r, arth_d, arth_i4b, arth_i8b END INTERFACE + +interface indexx +module procedure indexx_i4b +module procedure indexx_i8b +end interface + +interface swap +module procedure swap_i4b +module procedure swap_i8b +end interface + +interface unique +module procedure unique_i4b +module procedure unique_i8b +end interface + ! (everything private unless otherwise specifed) private public::arth @@ -45,23 +61,36 @@ FUNCTION arth_d(first,increment,n) end if END FUNCTION arth_d ! ------------------------------------------------------------------------------------------------ - FUNCTION arth_i(first,increment,n) + FUNCTION arth_i4b(first,increment,n) implicit none INTEGER(I4B), INTENT(IN) :: first,increment,n - INTEGER(I4B), DIMENSION(n) :: arth_i + INTEGER(I4B), DIMENSION(n) :: arth_i4b INTEGER(I4B) :: k - arth_i(1)=first + arth_i4b(1)=first if(n>1)then do k=2,n - arth_i(k) = arth_i(k-1) + increment + arth_i4b(k) = arth_i4b(k-1) + increment end do end if - END FUNCTION arth_i + END FUNCTION arth_i4b + ! ------------------------------------------------------------------------------------------------ + FUNCTION arth_i8b(first,increment,n) + implicit none + INTEGER(I8B), INTENT(IN) :: first,increment,n + INTEGER(I8B), DIMENSION(n) :: arth_i8b + INTEGER(I8B) :: k + arth_i8b(1)=first + if(n>1)then + do k=2,n + arth_i8b(k) = arth_i8b(k-1) + increment + end do + end if + END FUNCTION arth_i8b ! ************************************************************************************************* ! * sort function, used to sort numbers in ascending order ! ************************************************************************************************* - SUBROUTINE indexx(arr,index) + SUBROUTINE indexx_i4b(arr,index) IMPLICIT NONE INTEGER(I4B), DIMENSION(:), INTENT(IN) :: arr INTEGER(I4B), DIMENSION(:), INTENT(OUT) :: index @@ -136,16 +165,103 @@ SUBROUTINE icomp_xchg(i,j) j=swp end if END SUBROUTINE icomp_xchg - END SUBROUTINE indexx + END SUBROUTINE indexx_i4b + ! ------------------------------------------------------------------------------------------------ + SUBROUTINE indexx_i8b(arr,index) + IMPLICIT NONE + INTEGER(I8B), DIMENSION(:), INTENT(IN) :: arr + INTEGER(I4B), DIMENSION(:), INTENT(OUT) :: index + INTEGER(I4B), PARAMETER :: NN=15, NSTACK=50 + INTEGER(I8B) :: a + INTEGER(I4B) :: n,k,i,j,indext,jstack,l,r + INTEGER(I4B), DIMENSION(NSTACK) :: istack + n=size(arr) + index=arth(1,1,n) + jstack=0 + l=1 + r=n + do + if (r-l < NN) then + do j=l+1,r + indext=index(j) + a=arr(indext) + do i=j-1,1,-1 + if (arr(index(i)) <= a) exit + index(i+1)=index(i) + end do + index(i+1)=indext + end do + if (jstack == 0) RETURN + r=istack(jstack) + l=istack(jstack-1) + jstack=jstack-2 + else + k=(l+r)/2 + call swap(index(k),index(l+1)) + call icomp_xchg(index(l),index(r)) + call icomp_xchg(index(l+1),index(r)) + call icomp_xchg(index(l),index(l+1)) + i=l+1 + j=r + indext=index(l+1) + a=arr(indext) + do + do + i=i+1 + if (arr(index(i)) >= a) exit + end do + do + j=j-1 + if (arr(index(j)) <= a) exit + end do + if (j < i) exit + call swap(index(i),index(j)) + end do + index(l+1)=index(j) + index(j)=indext + jstack=jstack+2 + if (r-i+1 >= j-l) then + istack(jstack)=r + istack(jstack-1)=i + r=j-1 + else + istack(jstack)=j-1 + istack(jstack-1)=l + l=i + end if + end if + end do + CONTAINS + ! internal subroutine + SUBROUTINE icomp_xchg(i,j) + INTEGER(I4B), INTENT(INOUT) :: i,j + INTEGER(I4B) :: swp + if (arr(j) < arr(i)) then + swp=i + i=j + j=swp + end if + END SUBROUTINE icomp_xchg + END SUBROUTINE indexx_i8b + ! ************************************************************************************************ ! private subroutine - SUBROUTINE swap(a,b) + ! ************************************************************************************************ + SUBROUTINE swap_i4b(a,b) INTEGER(I4B), INTENT(INOUT) :: a,b INTEGER(I4B) :: dum dum=a a=b b=dum - END SUBROUTINE swap + END SUBROUTINE swap_i4b + ! ------------------------------------------------------------------------------------------------ + SUBROUTINE swap_i8b(a,b) + INTEGER(I8B), INTENT(INOUT) :: a,b + INTEGER(I8B) :: dum + dum=a + a=b + b=dum + END SUBROUTINE swap_i8b ! ************************************************************************************************ ! * findIndex: find the first index within a vector @@ -202,7 +318,7 @@ subroutine indexTrue(TF,pos) end subroutine indexTrue - SUBROUTINE unique(array, unq, idx) + SUBROUTINE unique_i4b(array, unq, idx) implicit none ! Input variables integer(i4b), intent(in) :: array(:) ! integer array including duplicated elements @@ -234,6 +350,40 @@ SUBROUTINE unique(array, unq, idx) idx = pack(arth(1,1,size(array)), flg_tmp) unq = unq_tmp(idx) - END SUBROUTINE unique + END SUBROUTINE unique_i4b + + SUBROUTINE unique_i8b(array, unq, idx) + implicit none + ! Input variables + integer(i8b), intent(in) :: array(:) ! integer array including duplicated elements + ! outpu variables + integer(i8b),allocatable,intent(out) :: unq(:) ! integer array including unique elements + integer(i4b),allocatable,intent(out) :: idx(:) ! integer array including unique element index + ! local + integer(i4b) :: ranked(size(array)) ! + integer(i8b) :: unq_tmp(size(array)) ! + logical(lgt) :: flg_tmp(size(array)) ! + integer(i4b) :: ix ! loop index, counter + integer(i8b) :: last_unique ! last unique element + + flg_tmp = .false. + call indexx(array, ranked) + + unq_tmp(ranked(1)) = array(ranked(1)) + flg_tmp(ranked(1)) = .true. + last_unique = array(ranked(1)) + do ix = 2,size(ranked) + if (last_unique==array(ranked(ix))) cycle + flg_tmp(ranked(ix)) = .true. + unq_tmp(ranked(ix)) = array(ranked(ix)) + last_unique = array(ranked(ix)) + end do + + allocate(unq(count(flg_tmp)),idx(count(flg_tmp))) + + idx = pack(arth(1,1,size(array)), flg_tmp) + unq = unq_tmp(idx) + + END SUBROUTINE unique_i8b end module nr_utility_module diff --git a/route/build/src/nrtype.f90 b/route/build/src/nrtype.f90 index 6f76067a..e478b3fc 100644 --- a/route/build/src/nrtype.f90 +++ b/route/build/src/nrtype.f90 @@ -1,5 +1,6 @@ MODULE nrtype ! variable types + INTEGER, PARAMETER :: I8B = SELECTED_INT_KIND(15) INTEGER, PARAMETER :: I4B = SELECTED_INT_KIND(9) INTEGER, PARAMETER :: I2B = SELECTED_INT_KIND(4) INTEGER, PARAMETER :: I1B = SELECTED_INT_KIND(2) diff --git a/route/build/src/popMetadat.f90 b/route/build/src/popMetadat.f90 index 2d4e49a9..8d3b00be 100644 --- a/route/build/src/popMetadat.f90 +++ b/route/build/src/popMetadat.f90 @@ -120,6 +120,7 @@ subroutine popMetadat(err,message) meta_SEG (ixSEG%totalArea ) = var_info('totalArea' , 'area above the bottom of the reach -- bas + ups' ,'m2' ,ixDims%seg , .false.) meta_SEG (ixSEG%basUnderLake ) = var_info('basUnderLake' , 'Area of basin under lake' ,'m2' ,ixDims%seg , .false.) meta_SEG (ixSEG%rchUnderLake ) = var_info('rchUnderLake' , 'Length of reach under lake' ,'m' ,ixDims%seg , .false.) + meta_SEG (ixSEG%Qtake ) = var_info('Qtake' , 'target abstraction(-)/injection(+)' ,'m3 s-1',ixDims%seg , .false.) meta_SEG (ixSEG%minFlow ) = var_info('minFlow' , 'minimum environmental flow' ,'m s-1' ,ixDims%seg , .false.) ! NTOPO varName varDesc varUnit, varType, varFile @@ -162,6 +163,7 @@ subroutine popMetadat(err,message) ! Impulse Response Function varName varDesc unit, varType, varDim, writeOut call meta_irf(ixIRF%qfuture)%init('irf_qfuture', 'future flow series', 'm3/sec' ,nf90_double, [ixStateDims%seg,ixStateDims%tdh_irf,ixStateDims%ens,ixStateDims%time] , .true.) + call meta_irf(ixIRF%irfVol) %init('irf_volume' , 'IRF reach volume' , 'm3' ,nf90_double, [ixStateDims%seg,ixStateDims%ens,ixStateDims%time] , .true.) ! Basin Impulse Response Function varName varDesc unit, varType, varDim, writeOut call meta_irf_bas(ixIRFbas%qfuture)%init('qfuture', 'future flow series', 'm3/sec' ,nf90_double, [ixStateDims%seg,ixStateDims%tdh,ixStateDims%ens,ixStateDims%time], .true.) diff --git a/route/build/src/process_ntopo.f90 b/route/build/src/process_ntopo.f90 index d5a719d7..98161a4b 100644 --- a/route/build/src/process_ntopo.f90 +++ b/route/build/src/process_ntopo.f90 @@ -11,6 +11,7 @@ module process_ntopo USE public_var, only : idSegOut ! ID for stream segment at the bottom of the subset ! options +USE public_var, only : qtakeOption ! option to compute network topology USE public_var, only : topoNetworkOption ! option to compute network topology USE public_var, only : computeReachList ! option to compute reach list USE public_var, only : hydGeometryOption ! option to obtain routing parameters @@ -230,6 +231,15 @@ subroutine augment_ntopo(& !print*, trim(message)//'PAUSE : '; read(*,*) ! ---------- Compute routing parameters -------------------------------------------------------------------- + do iSeg=1,nSeg + structSEG(iSeg)%var(ixSEG%minFlow)%dat(1) = 1.e-15_dp ! Minimum environmental flow + end do + + if(.not.qtakeOption)then + do iSeg=1,nSeg + structSEG(iSeg)%var(ixSEG%Qtake)%dat(1) = 0._dp ! no abstraction/injection + end do + end if ! compute hydraulic geometry (width and Manning's "n") if(hydGeometryOption==compute)then @@ -442,7 +452,10 @@ subroutine put_data_struct(nSeg, structSEG, structNTOPO, & RPARAM_in(iSeg)%UPSAREA = structSEG(iSeg)%var(ixSEG%upsArea)%dat(1) RPARAM_in(iSeg)%TOTAREA = structSEG(iSeg)%var(ixSEG%totalArea)%dat(1) - ! NOT USED: MINFLOW -- minimum environmental flow + ! Abstraction/Injection coefficient + RPARAM_in(iSeg)%QTAKE = structSEG(iSeg)%var(ixSEG%Qtake)%dat(1) + + ! MINFLOW -- minimum environmental flow RPARAM_in(iSeg)%MINFLOW = structSEG(iSeg)%var(ixSEG%minFlow)%dat(1) ! ----- network topology ----- diff --git a/route/build/src/process_time.f90 b/route/build/src/process_time.f90 deleted file mode 100644 index e6f353de..00000000 --- a/route/build/src/process_time.f90 +++ /dev/null @@ -1,100 +0,0 @@ -module process_time_module - -! data types -USE nrtype, only : i4b,dp ! variable types, etc. -USE nrtype, only : strLen ! length of characters -USE dataTypes, only : time ! time data type - -! subroutines: model time info -USE time_utils_module, only : extractTime ! get time from units string -USE time_utils_module, only : compJulday,& ! compute julian day - compJulday_noleap ! compute julian day for noleap calendar -USE time_utils_module, only : compcalday,& ! compute calendar date and time - compcalday_noleap ! compute calendar date and time for noleap calendar -implicit none - -! privacy -- everything private unless declared explicitly -private -public::process_time -public::process_calday - -contains - - ! ********************************************************************* - ! public subroutine: extract time information from the control information - ! ********************************************************************* - subroutine process_time(& - ! input - timeUnits, & ! time units string - calendar, & ! calendar - ! output - julianDate, & ! julian date - ierr, message) - implicit none - ! input - character(*) , intent(in) :: timeUnits ! time units string - character(*) , intent(in) :: calendar ! calendar string - ! output - real(dp) , intent(out) :: julianDate ! julian date - integer(i4b) , intent(out) :: ierr ! error code - character(*) , intent(out) :: message ! error message - ! -------------------------------------------------------------------------------------------------------------- - ! local variables - type(time) :: timeStruct ! time data structure - character(len=strLen) :: cmessage ! error message of downwind routine - ! initialize error control - ierr=0; message='process_time/' - - ! extract time from the units string - call extractTime(timeUnits,timeStruct%iy,timeStruct%im,timeStruct%id,timeStruct%ih,timeStruct%imin,timeStruct%dsec,ierr,cmessage) - if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - - ! calculate the julian day for the start time - select case(trim(calendar)) - case ('noleap','365_day') - call compjulday_noleap(timeStruct%iy,timeStruct%im,timeStruct%id,timeStruct%ih,timeStruct%imin,timeStruct%dsec,julianDate,ierr,cmessage) - case ('standard','gregorian','proleptic_gregorian') - call compjulday(timeStruct%iy,timeStruct%im,timeStruct%id,timeStruct%ih,timeStruct%imin,timeStruct%dsec,julianDate,ierr,cmessage) - case default; ierr=20; message=trim(message)//trim(calendar)//': calendar invalid; accept either noleap, 365_day, standard, gregorian, or proleptic_gregorian'; return - end select - if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - - end subroutine process_time - - ! ********************************************************************* - ! public subroutine: extract calendar date/time information from julian day - ! ********************************************************************* - subroutine process_calday(& - ! input - julianDate, & ! julian date - calendar, & ! calendar - ! output - datetime, & ! calendar - ierr, message) - implicit none - ! input - real(dp) , intent(in) :: julianDate ! julian date - character(*) , intent(in) :: calendar ! calendar string - ! output - type(time) , intent(out) :: datetime ! time data structure - integer(i4b) , intent(out) :: ierr ! error code - character(*) , intent(out) :: message ! error message - ! -------------------------------------------------------------------------------------------------------------- - ! local variables - character(len=strLen) :: cmessage ! error message of downwind routine - ! initialize error control - ierr=0; message='process_calday/' - - ! calculate the julian day for the start time - select case(trim(calendar)) - case ('noleap','365_day') - call compcalday_noleap(julianDate, datetime%iy,datetime%im,datetime%id,datetime%ih,datetime%imin,datetime%dsec, ierr, cmessage) - case ('standard','gregorian','proleptic_gregorian') - call compcalday(julianDate, datetime%iy,datetime%im,datetime%id,datetime%ih,datetime%imin,datetime%dsec, ierr, cmessage) - case default; ierr=20; message=trim(message)//trim(calendar)//': calendar invalid; accept either noleap, 365_day, standard, gregorian, or proleptic_gregorian'; return - end select - if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - - end subroutine process_calday - -end module process_time_module diff --git a/route/build/src/public_var.f90 b/route/build/src/public_var.f90 index 94951aad..401c491c 100644 --- a/route/build/src/public_var.f90 +++ b/route/build/src/public_var.f90 @@ -9,7 +9,7 @@ module public_var save ! ---------- mizuRoute version ------------------------------------------------------------------- - character(len=strLen), parameter, public :: mizuRouteVersion='v1.2' + character(len=strLen), parameter, public :: mizuRouteVersion='v1.2.1' ! ---------- common constants --------------------------------------------------------------------- @@ -69,9 +69,10 @@ module public_var ! Control file variables ! DIRECTORIES - character(len=strLen),public :: ancil_dir = '' ! directory containing ancillary data - character(len=strLen),public :: input_dir = '' ! directory containing input data - character(len=strLen),public :: output_dir = '' ! directory containing output data + character(len=strLen),public :: ancil_dir = '' ! directory containing ancillary data (network, mapping, namelist) + character(len=strLen),public :: input_dir = '' ! directory containing input runoff netCDF + character(len=strLen),public :: output_dir = '' ! directory for routed flow output (netCDF) + character(len=strLen),public :: restart_dir = charMissing ! directory for restart output (netCDF) ! SIMULATION TIME character(len=strLen),public :: simStart = '' ! date string defining the start of the simulation character(len=strLen),public :: simEnd = '' ! date string defining the end of the simulation @@ -93,6 +94,7 @@ module public_var character(len=strLen),public :: units_qsim = '' ! units of simulated runoff data real(dp) ,public :: dt = realMissing ! time step (seconds) real(dp) ,public :: ro_fillvalue = realMissing ! fillvalue used for runoff depth variable + logical(lgt) ,public :: userRunoffFillvalue = .false. ! true -> runoff depth fillvalue used in netcdf is specified here, otherwise -> false ! RUNOFF REMAPPING logical(lgt),public :: is_remap = .false. ! logical whether or not runnoff needs to be mapped to river network HRU character(len=strLen),public :: fname_remap = '' ! runoff mapping netCDF name @@ -108,12 +110,16 @@ module public_var character(len=strLen),public :: case_name = '' ! name of simulation. used as head of model output and restart file character(len=strLen),public :: newFileFrequency = 'annual' ! frequency for new output files (day, month, annual, single) ! STATES - character(len=strLen),public :: restart_write = 'never' ! restart write option: never-> N[n]ever write, L[l]ast -> write at last time step, S[s]pecified + character(len=strLen),public :: restart_write = 'never' ! restart write option: N[n]ever-> never write, L[l]ast -> write at last time step, S[s]pecified, Monthly, Daily character(len=strLen),public :: restart_date = charMissing ! specifed restart date + integer(i4b) ,public :: restart_month = 1 ! restart periodic month. Default Jan (write every January of year) + integer(i4b) ,public :: restart_day = 1 ! restart periodic day. Default 1st (write every 1st of month) + integer(i4b) ,public :: restart_hour = 0 ! restart periodic hour. Default 0hr (write every 00 hr of day) character(len=strLen),public :: fname_state_in = charMissing ! name of state file ! SPATIAL CONSTANT PARAMETERS character(len=strLen),public :: param_nml = '' ! name of the namelist file ! USER OPTIONS + logical(lgt) ,public :: qtakeOption = .false. ! option for abstraction/injection integer(i4b) ,public :: hydGeometryOption = compute ! option for hydraulic geometry calculations (0=read from file, 1=compute) integer(i4b) ,public :: topoNetworkOption = compute ! option for network topology calculations (0=read from file, 1=compute) integer(i4b) ,public :: computeReachList = compute ! option to compute list of upstream reaches (0=do not compute, 1=compute) diff --git a/route/build/src/read_control.f90 b/route/build/src/read_control.f90 index 6ea7fe19..1b3de3bc 100644 --- a/route/build/src/read_control.f90 +++ b/route/build/src/read_control.f90 @@ -97,9 +97,10 @@ subroutine read_control(ctl_fname, err, message) select case(trim(cName)) ! DIRECTORIES - case(''); ancil_dir = trim(cData) ! directory containing ancillary data - case(''); input_dir = trim(cData) ! directory containing input data - case(''); output_dir = trim(cData) ! directory containing output data + case(''); ancil_dir = trim(cData) ! directory containing ancillary data (network, mapping, namelist) + case(''); input_dir = trim(cData) ! directory containing input runoff netCDF + case(''); output_dir = trim(cData) ! directory for routed flow output (netCDF) + case(''); restart_dir = trim(cData) ! directory for restart output (netCDF) ! SIMULATION TIME case(''); simStart = trim(cData) ! date string defining the start of the simulation case(''); simEnd = trim(cData) ! date string defining the end of the simulation @@ -120,7 +121,9 @@ subroutine read_control(ctl_fname, err, message) case(''); dname_ylat = trim(cData) ! name of y (i,lat) dimension case(''); units_qsim = trim(cData) ! units of runoff case(''); read(cData,*,iostat=io_error) dt ! time interval of the gridded runoff - case(''); read(cData,*,iostat=io_error) ro_fillvalue ! fillvalue used for runoff depth variable + case('') + read(cData,*,iostat=io_error) ro_fillvalue ! fillvalue used for runoff depth variable + userRunoffFillvalue = .true. ! true -> runoff depth fillvalue used in netcdf is specified here, otherwise -> false ! RUNOFF REMAPPING case(''); read(cData,*,iostat=io_error) is_remap ! logical whether or not runnoff needs to be mapped to river network HRU case(''); fname_remap = trim(cData) ! name of runoff mapping netCDF @@ -135,13 +138,17 @@ subroutine read_control(ctl_fname, err, message) ! ROUTED FLOW OUTPUT case(''); case_name = trim(cData) ! name of simulation. used as head of model output and restart file case(''); newFileFrequency = trim(cData) ! frequency for new output files (day, month, annual, single) - ! STATES - case(''); restart_write = trim(cData) ! restart write option: N[n]ever, L[l]ast - case(''); restart_date = trim(cData) ! specified restart date, yyyy-mm-dd (hh:mm:ss) + ! RESTART + case(''); restart_write = trim(cData) ! restart write option: N[n]ever, L[l]ast, S[s]pecified, Monthly, Daily + case(''); restart_date = trim(cData) ! specified restart date, yyyy-mm-dd (hh:mm:ss) for Specified option + case(''); read(cData,*,iostat=io_error) restart_month ! restart periodic month + case(''); read(cData,*,iostat=io_error) restart_day ! restart periodic day + case(''); read(cData,*,iostat=io_error) restart_hour ! restart periodic hour case(''); fname_state_in = trim(cData) ! filename for the channel states ! SPATIAL CONSTANT PARAMETERS case(''); param_nml = trim(cData) ! name of namelist including routing parameter value ! USER OPTIONS: Define options to include/skip calculations + case(''); read(cData,*,iostat=io_error) qtakeOption ! option for abstraction/injection option case(''); read(cData,*,iostat=io_error) hydGeometryOption ! option for hydraulic geometry calculations (0=read from file, 1=compute) case(''); read(cData,*,iostat=io_error) topoNetworkOption ! option for network topology calculations (0=read from file, 1=compute) case(''); read(cData,*,iostat=io_error) computeReachList ! option to compute list of upstream reaches (0=do not compute, 1=compute) @@ -186,6 +193,7 @@ subroutine read_control(ctl_fname, err, message) case('' ); meta_SEG (ixSEG%upsArea )%varName =trim(cData) ! area above the top of the reach -- zero if headwater (m2) case('' ); meta_SEG (ixSEG%basUnderLake )%varName =trim(cData) ! Area of basin under lake (m2) case('' ); meta_SEG (ixSEG%rchUnderLake )%varName =trim(cData) ! Length of reach under lake (m) + case('' ); meta_SEG (ixSEG%Qtake )%varName =trim(cData) ! abstraction(-)/injection(+) (m3 s-1) case('' ); meta_SEG (ixSEG%minFlow )%varName =trim(cData) ! minimum environmental flow ! network topology case('' ); meta_NTOPO (ixNTOPO%hruContribIx )%varName =trim(cData) ! indices of the vector of HRUs that contribute flow to each segment @@ -221,6 +229,11 @@ subroutine read_control(ctl_fname, err, message) end do ! looping through lines in the control file + ! ---------- directory option --------------------------------------------------------------------- + if (trim(restart_dir)==charMissing) then + restart_dir = output_dir + endif + ! ---------- control river network writing option --------------------------------------------------------------------- ! Case1- river network subset mode (idSegOut>0): Write the network variables read from file over only upstream network specified idSegOut diff --git a/route/build/src/read_remap.f90 b/route/build/src/read_remap.f90 index 6665c34c..8cf15fa8 100644 --- a/route/build/src/read_remap.f90 +++ b/route/build/src/read_remap.f90 @@ -6,8 +6,10 @@ module read_remap USE public_var ! Netcdf -use io_netcdf, only:get_nc -use io_netcdf, only:get_nc_dim_len +USE io_netcdf, ONLY: open_nc +USE io_netcdf, ONLY: close_nc +USE io_netcdf, ONLY: get_nc +USE io_netcdf, ONLY: get_nc_dim_len implicit none @@ -34,20 +36,23 @@ subroutine get_remap_data(fname, & ! input: file name integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! local variables + integer(i4b) :: ncidMapping ! mapping netcdf id integer(i4b) :: iVar ! index of variables integer(i4b) :: nHRU ! number of HRU in mapping files (this should match up with river network hru) integer(i4b) :: nData ! number of data (weight, runoff hru id) in mapping files character(len=strLen) :: cmessage ! error message from subroutine - ! initialize error control ierr=0; message='get_remap_data/' + call open_nc(fname, 'r', ncidMapping, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + ! get the number of HRUs - call get_nc_dim_len(fname, dname_hru_remap, nHRU, ierr, cmessage) + call get_nc_dim_len(ncidMapping, dname_hru_remap, nHRU, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif ! get the number of spatial elements in the runoff file - call get_nc_dim_len(fname, dname_data_remap, nData, ierr, cmessage) + call get_nc_dim_len(ncidMapping, dname_data_remap, nData, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif ! allocate space for info in mapping file @@ -70,12 +75,12 @@ subroutine get_remap_data(fname, & ! input: file name if (nSpatial(2) == integerMissing .and. iVar >= 5) cycle if (nSpatial(2) /= integerMissing .and. iVar == 4) cycle select case(iVar) - case(1); call get_nc(fname, vname_hruid_in_remap, remap_data_in%hru_id, 1, nHRU, ierr, cmessage) - case(2); call get_nc(fname, vname_num_qhru, remap_data_in%num_qhru, 1, nHRU, ierr, cmessage) - case(3); call get_nc(fname, vname_weight, remap_data_in%weight, 1, nData, ierr, cmessage) - case(4); call get_nc(fname, vname_qhruid, remap_data_in%qhru_id, 1, nData, ierr, cmessage) - case(5); call get_nc(fname, vname_i_index, remap_data_in%i_index, 1, nData, ierr, cmessage) - case(6); call get_nc(fname, vname_j_index, remap_data_in%j_index, 1, nData, ierr, cmessage) + case(1); call get_nc(ncidMapping, vname_hruid_in_remap, remap_data_in%hru_id, 1, nHRU, ierr, cmessage) + case(2); call get_nc(ncidMapping, vname_num_qhru, remap_data_in%num_qhru, 1, nHRU, ierr, cmessage) + case(3); call get_nc(ncidMapping, vname_weight, remap_data_in%weight, 1, nData, ierr, cmessage) + case(4); call get_nc(ncidMapping, vname_qhruid, remap_data_in%qhru_id, 1, nData, ierr, cmessage) + case(5); call get_nc(ncidMapping, vname_i_index, remap_data_in%i_index, 1, nData, ierr, cmessage) + case(6); call get_nc(ncidMapping, vname_j_index, remap_data_in%j_index, 1, nData, ierr, cmessage) case default; ierr=20; message=trim(message)//'unable to find variable'; return end select if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif @@ -85,6 +90,9 @@ subroutine get_remap_data(fname, & ! input: file name call check_remap_data(remap_data_in, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + call close_nc(ncidMapping, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + end subroutine ! ***** @@ -109,8 +117,7 @@ subroutine check_remap_data(remap_data_in, & ! inout: data structure to remap integer(i4b) :: nZero ! number of array variables with zero logical(lgt), allocatable :: logical_array(:) ! real(dp), allocatable :: real_array(:) ! - integer(i4b), allocatable :: int_array(:) ! - character(len=strLen) :: cmessage ! error message from subroutine + integer(i8b), allocatable :: int_array(:) ! ! initialize error control ierr=0; message='check_remap_data/' diff --git a/route/build/src/read_restart.f90 b/route/build/src/read_restart.f90 index 4e897dae..261e6416 100644 --- a/route/build/src/read_restart.f90 +++ b/route/build/src/read_restart.f90 @@ -1,8 +1,12 @@ MODULE read_restart + ! Moudle wide external modules -USE nrtype, only: i4b, dp, & - strLen +USE nrtype, only: i4b, dp, strLen USE public_var +USE io_netcdf, ONLY: open_nc +USE io_netcdf, ONLY: close_nc +USE io_netcdf, ONLY: get_nc +USE io_netcdf, ONLY: get_nc_dim_len implicit none @@ -20,14 +24,11 @@ SUBROUTINE read_state_nc(& opt, & ! input: which routing options T0, T1, & ! output: start and end time [sec] ierr, message) ! Output: error control - ! External module - USE io_netcdf, ONLY: get_nc, & - get_nc_dim_len + USE dataTypes, ONLY: states - ! meta data USE globalData, ONLY: meta_stateDims ! dimension for state variables - ! Named variables USE var_lookup, ONLY: ixStateDims, nStateDims + implicit none ! input variables character(*), intent(in) :: fname ! filename @@ -38,6 +39,7 @@ SUBROUTINE read_state_nc(& integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! local variables + integer(i4b) :: ncidRestart ! restart netcdf id real(dp) :: TB(2) ! 2 element-time bound vector type(states) :: state(0:2) ! temporal state data structures -currently 2 river routing scheme + basin IRF routing integer(i4b) :: nSeg,nens ! dimenion sizes @@ -46,20 +48,22 @@ SUBROUTINE read_state_nc(& integer(i4b) :: jDim ! index loops for dimension character(len=strLen) :: cmessage ! error message of downwind routine - ! initialize error control ierr=0; message='read_state_nc/' ! get Dimension sizes ! For common dimension/variables - seg id, time, time-bound ----------- ixDim_common = (/ixStateDims%seg, ixStateDims%ens, ixStateDims%time, ixStateDims%tbound/) + call open_nc(fname, 'r', ncidRestart, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + do jDim=1,size(ixDim_common) associate (ixDim_tmp => ixDim_common(jDim)) select case(ixDim_tmp) - case(ixStateDims%seg); call get_nc_dim_len(fname, trim(meta_stateDims(ixDim_tmp)%dimName), nSeg, ierr, cmessage) - case(ixStateDims%ens); call get_nc_dim_len(fname, trim(meta_stateDims(ixDim_tmp)%dimName), nens, ierr, cmessage) - case(ixStateDims%time); call get_nc_dim_len(fname, trim(meta_stateDims(ixDim_tmp)%dimName), nTime, ierr, cmessage) - case(ixStateDims%tbound); call get_nc_dim_len(fname, trim(meta_stateDims(ixDim_tmp)%dimName), ntbound, ierr, cmessage) + case(ixStateDims%seg); call get_nc_dim_len(ncidRestart, meta_stateDims(ixDim_tmp)%dimName, nSeg, ierr, cmessage) + case(ixStateDims%ens); call get_nc_dim_len(ncidRestart, meta_stateDims(ixDim_tmp)%dimName, nens, ierr, cmessage) + case(ixStateDims%time); call get_nc_dim_len(ncidRestart, meta_stateDims(ixDim_tmp)%dimName, nTime, ierr, cmessage) + case(ixStateDims%tbound); call get_nc_dim_len(ncidRestart, meta_stateDims(ixDim_tmp)%dimName, ntbound, ierr, cmessage) case default; ierr=20; message=trim(message)//'unable to identify dimension name index'; return end select if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif @@ -68,13 +72,15 @@ SUBROUTINE read_state_nc(& ! Read variables ! time bound - call get_nc(fname,'time_bound',TB(:), 1, 2, ierr, cmessage) + call get_nc(ncidRestart,'time_bound',TB(:), 1, 2, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif T0=TB(1); T1=TB(2) ! routing specific variables - call read_IRFbas_state(ierr, cmessage) - if(ierr/=0)then; message=trim(message)//trim(cmessage); return;endif + if (doesBasinRoute == 1) then + call read_IRFbas_state(ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return;endif + endif if (opt==allRoutingMethods .or. opt==kinematicWave) then call read_KWT_state(ierr, cmessage) @@ -86,6 +92,9 @@ SUBROUTINE read_state_nc(& if(ierr/=0)then; message=trim(message)//trim(cmessage);return; endif end if + call close_nc(ncidRestart, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + CONTAINS SUBROUTINE read_IRFbas_state(ierr, message1) @@ -106,7 +115,7 @@ SUBROUTINE read_IRFbas_state(ierr, message1) ! initialize error control ierr=0; message1='read_IRFbas_state/' - call get_nc_dim_len(fname, trim(meta_stateDims(ixStateDims%tdh)%dimName), ntdh, ierr, cmessage) + call get_nc_dim_len(ncidRestart, meta_stateDims(ixStateDims%tdh)%dimName, ntdh, ierr, cmessage) if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif allocate(state(0)%var(nVarsIRFbas), stat=ierr, errmsg=cmessage) @@ -126,8 +135,8 @@ SUBROUTINE read_IRFbas_state(ierr, message1) do iVar=1,nVarsIRFbas select case(iVar) - case(ixIRFbas%q); call get_nc(fname, meta_irf_bas(iVar)%varName, state(0)%var(iVar)%array_2d_dp, (/1,1/), (/nSeg,nens/), ierr, cmessage) - case(ixIRFbas%qfuture); call get_nc(fname, meta_irf_bas(iVar)%varName, state(0)%var(iVar)%array_3d_dp, (/1,1,1/), (/nSeg,ntdh,nens/), ierr, cmessage) + case(ixIRFbas%q); call get_nc(ncidRestart, meta_irf_bas(iVar)%varName, state(0)%var(iVar)%array_2d_dp, (/1,1/), (/nSeg,nens/), ierr, cmessage) + case(ixIRFbas%qfuture); call get_nc(ncidRestart, meta_irf_bas(iVar)%varName, state(0)%var(iVar)%array_3d_dp, (/1,1,1/), (/nSeg,ntdh,nens/), ierr, cmessage) case default; ierr=20; message1=trim(message1)//'unable to identify basin IRF variable index for nc writing'; return end select if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif @@ -169,10 +178,10 @@ SUBROUTINE read_IRF_state(ierr, message1) integer(i4b) :: iVar,iens,iSeg ! index loops for variables, ensembles, reaches respectively integer(i4b), allocatable :: numQF(:,:) ! number of future Q time steps for each ensemble and segment integer(i4b) :: ntdh_irf ! dimenion sizes - ! initialize error control + ierr=0; message1='read_IRF_state/' - call get_nc_dim_len(fname, trim(meta_stateDims(ixStateDims%tdh_irf)%dimName), ntdh_irf, ierr, cmessage) + call get_nc_dim_len(ncidRestart, meta_stateDims(ixStateDims%tdh_irf)%dimName, ntdh_irf, ierr, cmessage) if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif allocate(state(impulseResponseFunc)%var(nVarsIRF), stat=ierr, errmsg=cmessage) @@ -185,19 +194,21 @@ SUBROUTINE read_IRF_state(ierr, message1) select case(iVar) case(ixIRF%qfuture); allocate(state(impulseResponseFunc)%var(iVar)%array_3d_dp(nSeg, ntdh_irf, nens), stat=ierr) + case(ixIRF%irfVol); allocate(state(impulseResponseFunc)%var(iVar)%array_2d_dp(nSeg, nens), stat=ierr) case default; ierr=20; message1=trim(message1)//'unable to identify variable index'; return end select if(ierr/=0)then; message1=trim(message1)//'problem allocating space for IRF routing state '//trim(meta_irf(iVar)%varName); return; endif end do - call get_nc(fname,'numQF',numQF,(/1,1/),(/nSeg,nens/),ierr,cmessage) + call get_nc(ncidRestart,'numQF',numQF,(/1,1/),(/nSeg,nens/),ierr,cmessage) if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif do iVar=1,nVarsIRF select case(iVar) - case(ixIRF%qfuture); call get_nc(fname, meta_irf(iVar)%varName, state(impulseResponseFunc)%var(iVar)%array_3d_dp, (/1,1,1/), (/nSeg,ntdh_irf,nens/), ierr, cmessage) + case(ixIRF%qfuture); call get_nc(ncidRestart, meta_irf(iVar)%varName, state(impulseResponseFunc)%var(iVar)%array_3d_dp, (/1,1,1/), (/nSeg,ntdh_irf,nens/), ierr, cmessage) + case(ixIRF%irfVol); call get_nc(ncidRestart, meta_irf(iVar)%varName, state(impulseResponseFunc)%var(iVar)%array_2d_dp, (/1,1/), (/nSeg, nens/), ierr, cmessage) case default; ierr=20; message1=trim(message1)//'unable to identify IRF variable index for nc reading'; return end select if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif @@ -213,7 +224,8 @@ SUBROUTINE read_IRF_state(ierr, message1) do iVar=1,nVarsIRF select case(iVar) - case(ixIRF%qfuture); RCHFLX(iens,iSeg)%QFUTURE_IRF = state(impulseResponseFunc)%var(iVar)%array_3d_dp(iSeg,1:numQF(iens,iSeg),iens) + case(ixIRF%qfuture); RCHFLX(iens,iSeg)%QFUTURE_IRF = state(impulseResponseFunc)%var(iVar)%array_3d_dp(iSeg,1:numQF(iens,iSeg),iens) + case(ixIRF%irfVol); RCHFLX(iens,iSeg)%REACH_VOL(1) = state(impulseResponseFunc)%var(iVar)%array_2d_dp(iSeg,iens) case default; ierr=20; message1=trim(message1)//'unable to identify variable index'; return end select @@ -249,7 +261,7 @@ SUBROUTINE read_KWT_state(ierr, message1) if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif ! get Dimension sizes - call get_nc_dim_len(fname, trim(meta_stateDims(ixStateDims%wave)%dimName), nwave, ierr, cmessage) + call get_nc_dim_len(ncidRestart, meta_stateDims(ixStateDims%wave)%dimName, nwave, ierr, cmessage) if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif do iVar=1,nVarsKWT @@ -263,16 +275,16 @@ SUBROUTINE read_KWT_state(ierr, message1) if(ierr/=0)then; message1=trim(message1)//'problem allocating space for KWT routing state '//trim(meta_kwt(iVar)%varName); return; endif end do - call get_nc(fname,'numWaves',numWaves, (/1,1/), (/nSeg,nens/), ierr, cmessage) + call get_nc(ncidRestart,'numWaves',numWaves, (/1,1/), (/nSeg,nens/), ierr, cmessage) if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif do iVar=1,nVarsKWT select case(iVar) case(ixKWT%routed) - call get_nc(fname,trim(meta_kwt(iVar)%varName), state(kinematicWave)%var(iVar)%array_3d_dp, (/1,1,1/), (/nSeg,nwave,nens/), ierr, cmessage) + call get_nc(ncidRestart, meta_kwt(iVar)%varName, state(kinematicWave)%var(iVar)%array_3d_dp, (/1,1,1/), (/nSeg,nwave,nens/), ierr, cmessage) case(ixKWT%tentry, ixKWT%texit, ixKWT%qwave, ixKWT%qwave_mod) - call get_nc(fname,trim(meta_kwt(iVar)%varName), state(kinematicWave)%var(iVar)%array_3d_dp, (/1,1,1/), (/nSeg,nwave,nens/), ierr, cmessage) + call get_nc(ncidRestart, meta_kwt(iVar)%varName, state(kinematicWave)%var(iVar)%array_3d_dp, (/1,1,1/), (/nSeg,nwave,nens/), ierr, cmessage) case default; ierr=20; message1=trim(message)//'unable to identify KWT variable index for nc reading'; return end select if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif diff --git a/route/build/src/read_runoff.f90 b/route/build/src/read_runoff.f90 index 5271fe35..91edb115 100644 --- a/route/build/src/read_runoff.f90 +++ b/route/build/src/read_runoff.f90 @@ -4,6 +4,7 @@ module read_runoff USE nrtype USE public_var USE io_netcdf, only:open_nc +USE io_netcdf, only:close_nc USE io_netcdf, only:get_nc USE io_netcdf, only:get_var_attr USE io_netcdf, only:check_attr @@ -41,84 +42,95 @@ subroutine read_runoff_metadata(& integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! local variables - integer(i4b) :: ncid ! netcdf id + integer(i4b) :: ncidRunoff ! netcdf id integer(i4b) :: ivarID ! variable id integer(i4b) :: nDims ! number of dimension in runoff file character(len=strLen) :: cmessage ! error message from subroutine - ! initialize error control + ierr=0; message='read_runoff_metadata/' - ! open NetCDF file - call open_nc(trim(fname), 'r', ncid, ierr, cmessage) + call open_nc(fname, 'r', ncidRunoff, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif ! get the ID of runoff variable - ierr = nf90_inq_varid(ncid, trim(vname_qsim), ivarID) + ierr = nf90_inq_varid(ncidRunoff, trim(vname_qsim), ivarID) if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif ! get the number of dimensions - must be 2D(hru, time) or 3D(y, x, time) - ierr= nf90_inquire_variable(ncid, ivarID, ndims = nDims) + ierr= nf90_inquire_variable(ncidRunoff, ivarID, ndims = nDims) if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif ! get runoff metadata select case( nDims ) - case(2); call read_1D_runoff_metadata(fname, runoff_data_in, timeUnits, calendar, ierr, cmessage) - case(3); call read_2D_runoff_metadata(fname, runoff_data_in, timeUnits, calendar, ierr, cmessage) + case(2); call read_1D_runoff_metadata(ncidRunoff, runoff_data_in, timeUnits, calendar, ierr, cmessage) + case(3); call read_2D_runoff_metadata(ncidRunoff, runoff_data_in, timeUnits, calendar, ierr, cmessage) case default; ierr=20; message=trim(message)//'runoff input must be 2-dimension (e.g, [time, hru]) or 3-dimension (e.g., [time, lat, lon]'; return end select if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + call close_nc(ncidRunoff, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + end subroutine read_runoff_metadata ! ***** ! private subroutine: get 2D runoff (hru, time) metadata... ! ****************************************** - subroutine read_1D_runoff_metadata(& - ! input - fname , & ! filename - ! output - runoff_data_in , & ! runoff data structure - timeUnits , & ! time units - calendar , & ! calendar - ! error control - ierr, message) ! output: error control + subroutine read_1D_runoff_metadata(ncidRunoff , & ! input: netcdf id + runoff_data_in , & ! output: runoff data structure + timeUnits , & ! output: time units + calendar , & ! output: calendar + ierr, message) ! output: error control implicit none ! input variables - character(*), intent(in) :: fname ! filename + integer(i4b), intent(in) :: ncidRunoff ! netcdf id ! output variables - type(runoff), intent(out) :: runoff_data_in ! runoff for one time step for all HRUs + type(runoff), intent(out) :: runoff_data_in ! runoff for one time step for all HRUs character(*), intent(out) :: timeUnits ! time units character(*), intent(out) :: calendar ! calendar ! error control integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! local variables + logical(lgt) :: existFillVal character(len=strLen) :: cmessage ! error message from subroutine - ! initialize error control + ierr=0; message='read_1D_runoff_metadata/' runoff_data_in%nSpace(2) = integerMissing ! get the number of HRUs - call get_nc_dim_len(fname, trim(dname_hruid), runoff_data_in%nSpace(1), ierr, cmessage) + call get_nc_dim_len(ncidRunoff, trim(dname_hruid), runoff_data_in%nSpace(1), ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif ! get number of time steps from the runoff file - call get_nc_dim_len(fname, trim(dname_time), runoff_data_in%nTime, ierr, cmessage) + call get_nc_dim_len(ncidRunoff, trim(dname_time), runoff_data_in%nTime, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif ! get the time units if (trim(timeUnits) == charMissing) then - call get_var_attr(fname, trim(vname_time), 'units', timeUnits, ierr, cmessage) + call get_var_attr(ncidRunoff, trim(vname_time), 'units', timeUnits, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif end if ! get the calendar if (trim(calendar) == charMissing) then - call get_var_attr(fname, trim(vname_time), 'calendar', calendar, ierr, cmessage) + call get_var_attr(ncidRunoff, trim(vname_time), 'calendar', calendar, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif end if + ! get the _fill_values for runoff variable + if (.not.userRunoffFillvalue) then + existFillVal = check_attr(ncidRunoff, vname_qsim, '_FillValue') + if (existFillVal) then + call get_var_attr(ncidRunoff, vname_qsim, '_FillValue', ro_fillvalue, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + else + write(iulog,'(a)') 'WARNING: User did not provide runoff fillvalue in control file nor runoff netcdf does not have fillvalue in attribute.' + write(iulog,'(a,x,F8.1)') ' Default missing values used is', ro_fillvalue + end if + end if + ! allocate space for hru_id allocate(runoff_data_in%hru_id(runoff_data_in%nSpace(1)), stat=ierr) if(ierr/=0)then; message=trim(message)//'problem allocating runoff_data_in%hruId'; return; endif @@ -128,7 +140,7 @@ subroutine read_1D_runoff_metadata(& if(ierr/=0)then; message=trim(message)//'problem allocating runoff_data_in%qsim'; return; endif ! get HRU ids from the runoff file - call get_nc(fname, vname_hruid, runoff_data_in%hru_id, 1, runoff_data_in%nSpace(1), ierr, cmessage) + call get_nc(ncidRunoff, vname_hruid, runoff_data_in%hru_id, 1, runoff_data_in%nSpace(1), ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif end subroutine read_1D_runoff_metadata @@ -136,52 +148,61 @@ end subroutine read_1D_runoff_metadata ! ***** ! private subroutine: get 3D runoff (lat, lon, time) metadata... ! ****************************************** - subroutine read_2D_runoff_metadata(& - ! input - fname , & ! filename - ! output - runoff_data_in , & ! runoff data structure - timeUnits , & ! time units - calendar , & ! calendar - ! error control - ierr, message) ! output: error control + subroutine read_2D_runoff_metadata(ncidRunoff , & ! input: netcdf id + runoff_data_in , & ! output: runoff data structure + timeUnits , & ! output: time units + calendar , & ! output: calendar + ierr, message) ! output: error control implicit none ! input variables - character(*), intent(in) :: fname ! filename + integer(i4b), intent(in) :: ncidRunoff ! netcdf id ! output variables - type(runoff), intent(out) :: runoff_data_in ! runoff for one time step for all HRUs + type(runoff), intent(out) :: runoff_data_in ! runoff for one time step for all HRUs character(*), intent(out) :: timeUnits ! time units character(*), intent(out) :: calendar ! calendar ! error control integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! local variables + logical(lgt) :: existFillVal character(len=strLen) :: cmessage ! error message from subroutine ! initialize error control ierr=0; message='read_2D_runoff_metadata/' ! get number of time steps from the runoff file - call get_nc_dim_len(fname, trim(dname_time), runoff_data_in%nTime, ierr, cmessage) + call get_nc_dim_len(ncidRunoff, trim(dname_time), runoff_data_in%nTime, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif ! get the time units if (trim(timeUnits) == charMissing) then - call get_var_attr(fname, trim(vname_time), 'units', timeUnits, ierr, cmessage) + call get_var_attr(ncidRunoff, trim(vname_time), 'units', timeUnits, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif end if ! get the calendar if (trim(calendar) == charMissing) then - call get_var_attr(fname, trim(vname_time), 'calendar', calendar, ierr, cmessage) + call get_var_attr(ncidRunoff, trim(vname_time), 'calendar', calendar, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif end if + ! get the _fill_values for runoff variable + if (.not.userRunoffFillvalue) then + existFillVal = check_attr(ncidRunoff, vname_qsim, '_FillValue') + if (existFillVal) then + call get_var_attr(ncidRunoff, vname_qsim, '_FillValue', ro_fillvalue, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + else + write(iulog,'(a)') 'WARNING: User did not provide runoff fillvalue in control file nor runoff netcdf does not have fillvalue attribute.' + write(iulog,'(a,x,F8.1)') ' Default missing values used is', ro_fillvalue + end if + end if + ! get size of ylat dimension - call get_nc_dim_len(fname, trim(dname_ylat), runoff_data_in%nSpace(1), ierr, cmessage) + call get_nc_dim_len(ncidRunoff, trim(dname_ylat), runoff_data_in%nSpace(1), ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif ! get size of xlon dimension - call get_nc_dim_len(fname, trim(dname_xlon), runoff_data_in%nSpace(2), ierr, cmessage) + call get_nc_dim_len(ncidRunoff, trim(dname_xlon), runoff_data_in%nSpace(2), ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif ! allocate space for simulated runoff. qSim2d = runoff(lon, lat) @@ -194,7 +215,7 @@ end subroutine read_2D_runoff_metadata ! ********************************************************************* ! public subroutine: read runoff data ! ********************************************************************* - subroutine read_runoff_data(fname, & ! input: filename + subroutine read_runoff_data(fname, & ! input: runoff netcdf name iTime, & ! input: time index runoff_data_in, & ! inout: runoff data structure ierr, message) ! output: error control @@ -208,32 +229,38 @@ subroutine read_runoff_data(fname, & ! input: filename integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! local variables + integer(i4b) :: ncidRunoff ! runoff netCDF ID character(len=strLen) :: cmessage ! error message from subroutine - ! initialize error control ierr=0; message='read_runoff_data/' + call open_nc(fname, 'r', ncidRunoff, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + if (runoff_data_in%nSpace(2) == integerMissing) then - call read_1D_runoff(fname, iTime, runoff_data_in%nSpace(1), runoff_data_in, ierr, cmessage) + call read_1D_runoff(ncidRunoff, iTime, runoff_data_in%nSpace(1), runoff_data_in, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif else - call read_2D_runoff(fname, iTime, runoff_data_in%nSpace, runoff_data_in, ierr, cmessage) + call read_2D_runoff(ncidRunoff, iTime, runoff_data_in%nSpace, runoff_data_in, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif endif + call close_nc(ncidRunoff, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + end subroutine read_runoff_data ! ********************************************************************* ! private subroutine: read 2D runoff data ! ********************************************************************* - subroutine read_1D_runoff(fname, & ! input: filename + subroutine read_1D_runoff(ncidRunoff, & ! input: runoff netcdf ID iTime, & ! input: time index nSpace, & ! input: size of HRUs runoff_data_in, & ! inout: runoff data structure ierr, message) ! output: error control implicit none ! input variables - character(*), intent(in) :: fname ! filename + integer(i4b), intent(in) :: ncidRunoff ! runoff netCDF ID integer(i4b), intent(in) :: iTime ! index of time element integer(i4b), intent(in) :: nSpace ! size of spatial dimensions ! input/output variables @@ -242,7 +269,8 @@ subroutine read_1D_runoff(fname, & ! input: filename integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! local variables - logical(lgt) :: existFillVal + integer(i4b) :: iStart(2) + integer(i4b) :: iCount(2) real(dp) :: dummy(nSpace,1) ! data read character(len=strLen) :: cmessage ! error message from subroutine @@ -250,20 +278,15 @@ subroutine read_1D_runoff(fname, & ! input: filename ierr=0; message='read_1D_runoff/' ! get the time data - call get_nc(trim(fname), vname_time, runoff_data_in%time, iTime, ierr, cmessage) + call get_nc(ncidRunoff, vname_time, runoff_data_in%time, iTime, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif ! get the simulated runoff data - call get_nc(trim(fname), vname_qsim, dummy, (/1,iTime/), (/nSpace,1/), ierr, cmessage) + iStart = [1,iTime] + iCount = [nSpace,1] + call get_nc(ncidRunoff, vname_qsim, dummy, iStart, iCount, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - ! get the _fill_values for runoff variable if exist - existFillVal = check_attr(trim(fname), vname_qsim, '_FillValue') - if (existFillval) then - call get_var_attr(trim(fname), vname_qsim, '_FillValue', ro_fillvalue, ierr, cmessage) - if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - end if - ! replace _fill_value with -999 for dummy where ( abs(dummy - ro_fillvalue) < verySmall ) dummy = realMissing @@ -275,14 +298,14 @@ end subroutine read_1D_runoff ! ********************************************************************* ! private subroutine: read 2D runoff data ! ********************************************************************* - subroutine read_2D_runoff(fname, & ! input: filename + subroutine read_2D_runoff(ncidRunoff, & ! input: runoff netcdf ID iTime, & ! input: time index nSpace, & ! input: size of HRUs runoff_data_in, & ! output: runoff data structure ierr, message) ! output: error control implicit none ! input variables - character(*), intent(in) :: fname ! filename + integer(i4b), intent(in) :: ncidRunoff ! runoff netCDF ID integer(i4b), intent(in) :: iTime ! index of time element integer(i4b), intent(in) :: nSpace(1:2) ! size of spatial dimensions ! input/output variables @@ -291,28 +314,23 @@ subroutine read_2D_runoff(fname, & ! input: filename integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! local variables - logical(lgt) :: existFillVal + integer(i4b) :: iStart(3) + integer(i4b) :: iCount(3) real(dp) :: dummy(nSpace(2),nSpace(1),1) ! data read character(len=strLen) :: cmessage ! error message from subroutine - ! initialize error control ierr=0; message='read_2D_runoff/' ! get the time data - call get_nc(trim(fname), vname_time, runoff_data_in%time, iTime, ierr, cmessage) + call get_nc(ncidRunoff, vname_time, runoff_data_in%time, iTime, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif ! get the simulated runoff data - call get_nc(trim(fname), vname_qsim, dummy, (/1,1,iTime/), (/nSpace(2), nSpace(1), 1/), ierr, cmessage) + iStart = [1,1,iTime] + iCount = [nSpace(2),nSpace(1),1] + call get_nc(ncidRunoff, vname_qsim, dummy, iStart, iCount, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - ! get the _fill_values for runoff variable - existFillVal = check_attr(trim(fname), vname_qsim, '_FillValue') - if (existFillval) then - call get_var_attr(trim(fname), vname_qsim, '_FillValue', ro_fillvalue, ierr, cmessage) - if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - end if - ! replace _fill_value with -999. for dummy where ( abs(dummy - ro_fillvalue) < verySmall ) dummy = realMissing diff --git a/route/build/src/read_streamSeg.f90 b/route/build/src/read_streamSeg.f90 index 0a42fe78..9a0e5373 100644 --- a/route/build/src/read_streamSeg.f90 +++ b/route/build/src/read_streamSeg.f90 @@ -142,6 +142,10 @@ subroutine getData(& ! ----------------------------------------------------------------------------------------------------------------- ! ---------- read in data ----------------------------------------------------------------------------------------- ! ----------------------------------------------------------------------------------------------------------------- + ! set flags if we want to turn on abstraction/injection option (require Qtake in network data) + if(qtakeOption)then + meta_SEG(ixSEG%Qtake)%varFile = .true. + endif ! set flags if we want to read hdraulic geometry from file if(hydGeometryOption==readFromFile)then diff --git a/route/build/src/remap.f90 b/route/build/src/remap.f90 index 74554eed..44764d03 100644 --- a/route/build/src/remap.f90 +++ b/route/build/src/remap.f90 @@ -206,8 +206,8 @@ subroutine remap_1D_runoff(runoff_data_in, remap_data_in, basinRunoff, ierr, mes cycle endif - !print*, 'remap_data_in%hru_id(iHRU), structHRU2seg(jHRU)%var(ixHRU2seg%hruId)%dat(1), remap_data_in%num_qhru(iHRU) = ', & - ! remap_data_in%hru_id(iHRU), structHRU2seg(jHRU)%var(ixHRU2seg%hruId)%dat(1), remap_data_in%num_qhru(iHRU) + !print*, 'remap_data_in%hru_id(iHRU), basinID(jHRU), remap_data_in%num_qhru(iHRU) = ', & + ! remap_data_in%hru_id(iHRU), basinID(jHRU), remap_data_in%num_qhru(iHRU) ! initialize the weighted average sumWeights = 0._dp @@ -296,6 +296,9 @@ subroutine sort_runoff(runoff_data_in, basinRunoff, ierr, message) ierr=0; message="sort_runoff/" + ! initialize zero at all the River network HRUs + basinRunoff(:) = 0._dp + ! loop through hrus in the runoff layer do iHRU=1,size(runoff_data_in%hru_ix) diff --git a/route/build/src/route_runoff.f90 b/route/build/src/route_runoff.f90 index e82cfa62..471ac29f 100644 --- a/route/build/src/route_runoff.f90 +++ b/route/build/src/route_runoff.f90 @@ -14,12 +14,14 @@ program route_runoff USE model_setup, only : init_data ! initialize river reach data USE model_setup, only : update_time ! Update simulation time information at each time step ! subroutines: routing -USE main_route_module, only : main_route ! +USE main_route_module, only : main_route ! main routing routine ! subroutines: model I/O USE get_runoff , only : get_hru_runoff ! USE write_simoutput, only : prep_output ! USE write_simoutput, only : output ! -USE write_restart, only : output_state ! write netcdf state output file +USE write_restart, only : main_restart ! write netcdf restart file +USE model_finalize, ONLY : finalize +USE model_finalize, ONLY : handle_err implicit none @@ -68,11 +70,9 @@ program route_runoff ! *********************************** do while (.not.finished) - ! prepare simulation output netCDF call prep_output(ierr, cmessage) if(ierr/=0) call handle_err(ierr, cmessage) - ! Get river network hru runoff at current time step call system_clock(startTime) call get_hru_runoff(ierr, cmessage) if(ierr/=0) call handle_err(ierr, cmessage) @@ -94,29 +94,14 @@ program route_runoff elapsedTime = real(endTime-startTime, kind(dp))/real(cr) write(*,"(A,1PG15.7,A)") ' elapsed-time [output] = ', elapsedTime, ' s' - ! write state netCDF - call output_state(ierr, cmessage) + call main_restart(ierr, cmessage) if(ierr/=0) call handle_err(ierr, cmessage) call update_time(finished, ierr, cmessage) if(ierr/=0) call handle_err(ierr, cmessage) -end do ! looping through time +end do -stop - -contains - - subroutine handle_err(err,message) - ! handle error codes - implicit none - integer(i4b),intent(in)::err ! error code - character(*),intent(in)::message ! error message - if(err/=0)then - print*,'FATAL ERROR: '//trim(message) - call flush(6) - stop - endif - end subroutine handle_err +call finalize() end program route_runoff diff --git a/route/build/src/time_utils.f90 b/route/build/src/time_utils.f90 index 31d58905..6d564820 100644 --- a/route/build/src/time_utils.f90 +++ b/route/build/src/time_utils.f90 @@ -26,13 +26,84 @@ module time_utils_module implicit none private public::extractTime -public::compjulday -public::compjulday_noleap -public::compcalday -public::compcalday_noleap +public::ndays_month +public::isLeapYear +public::compJulday +public::compJulday_noleap +public::compCalday +public::compCalday_noleap public::elapsedSec contains + ! ****************************************************************************************** + ! public function: check leap year or not + ! ****************************************************************************************** + logical(lgt) function isLeapYear(yr) + implicit none + integer(i4b), intent(in) :: yr + if (mod(yr, 4) == 0) then + if (mod(yr, 100) == 0) then + if (mod(yr, 400) == 0) then + isLeapYear = .True. + else + isLeapYear = .False. + end if + else + isLeapYear = .True. + end if + else + isLeapYear = .False. + end if + end function isLeapYear + + ! ****************************************************************************************** + ! public subroutine: get number of days within a month + ! ****************************************************************************************** + subroutine ndays_month(yr, mo, calendar, ndays, ierr, message) + implicit none + ! input + integer(i4b),intent(in) :: yr + integer(i4b),intent(in) :: mo + character(*),intent(in) :: calendar + ! output + integer(i4b) ,intent(out) :: ndays ! + integer(i4b), intent(out) :: ierr ! error code + character(len=strLen),intent(out) :: message ! error message + ! local variables + integer(i4b) :: yr_next, mo_next + real(dp) :: julday1, julday2 + character(len=strLen) :: cmessage ! error message of downwind routine + + ierr=0; message="ndays_month/" + + select case(trim(calendar)) + case ('standard','gregorian','proleptic_gregorian') + call compJulday(yr,mo,1,0,0,0._dp,julday1,ierr,cmessage) + case('noleap') + call compJulday_noleap(yr,mo,1,0,0,0._dp,julday1,ierr,cmessage) + case default; ierr=20; message=trim(message)//'calendar name: '//trim(calendar)//' invalid'; return + end select + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + + if (mo == 12) then + mo_next = 1 + yr_next = yr+1 + else + yr_next = yr + mo_next = mo+1 + end if + select case(trim(calendar)) + case ('standard','gregorian','proleptic_gregorian') + call compJulday(yr_next,mo_next,1,0,0,0._dp,julday2,ierr,cmessage) + case('noleap') + call compJulday_noleap(yr_next,mo_next,1,0,0,0._dp,julday2,ierr,cmessage) + case default; ierr=20; message=trim(message)//'calendar name: '//trim(calendar)//' invalid'; return + end select + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + + ndays = nint(julday2-julday1) + + end subroutine ndays_month ! ****************************************************************************************** ! public subroutine extractTime: extract year/month/day/hour/minute/second from units string @@ -137,9 +208,9 @@ end subroutine extractTime ! *************************************************************************************** - ! public subroutine compjulday: convert date to julian day (units of days) + ! public subroutine compJulday: convert date to julian day (units of days) ! *************************************************************************************** - subroutine compjulday(iyyy,mm,id,ih,imin,dsec,& ! input + subroutine compJulday(iyyy,mm,id,ih,imin,dsec,& ! input juldayss,err,message) ! output implicit none ! input variables @@ -157,7 +228,7 @@ subroutine compjulday(iyyy,mm,id,ih,imin,dsec,& ! input real(dp) :: jfrac ! fraction of julian day ! initialize errors - err=0; message="juldayss" + err=0; message="compJulday/" ! compute julian day jy=iyyy @@ -181,13 +252,13 @@ subroutine compjulday(iyyy,mm,id,ih,imin,dsec,& ! input ! and return the julian day, expressed in fraction of a day juldayss = real(julday,kind(dp)) + jfrac - end subroutine compjulday + end subroutine compJulday ! *************************************************************************************** - ! public subroutine compjulday: convert date to julian day (units of days) for noleap calendar + ! public subroutine compJulday: convert date to julian day (units of days) for noleap calendar ! reference: https://github.com/nmizukami/VIC/blob/VIC.5.0.0/vic/drivers/shared_all/src/vic_time.c ! *************************************************************************************** - subroutine compjulday_noleap(iyyy,mm,id,ih,imin,dsec,& ! input + subroutine compJulday_noleap(iyyy,mm,id,ih,imin,dsec,& ! input juldayss,err,message) ! output implicit none ! input variables @@ -204,7 +275,7 @@ subroutine compjulday_noleap(iyyy,mm,id,ih,imin,dsec,& ! input real(dp) :: dfrac ! fraction of day ! initialize errors - err=0; message="compjulday_noleap" + err=0; message="compJulday_noleap/" ! compute fraction of the day dfrac = real(id,kind(dp))+(real(ih,kind(dp))*secprhour + real(imin,kind(dp))*secprmin + dsec) / secprday @@ -220,14 +291,13 @@ subroutine compjulday_noleap(iyyy,mm,id,ih,imin,dsec,& ! input juldayss = real(days_per_yr*(iyyy_tmp + 4716)) + & real(floor(30.6001 * real(mm_tmp+1))) + dfrac - 1524.5 - end subroutine compjulday_noleap + end subroutine compJulday_noleap ! *************************************************************************************** ! public subroutine compgregcal: convert julian day (units of days) to calendar date ! source: https://en.wikipedia.org/wiki/Julian_day#Julian_or_Gregorian_calendar_from_Julian_day_number ! *************************************************************************************** - - subroutine compcalday(julday, & !input + subroutine compCalday(julday, & !input iyyy,mm,id,ih,imin,dsec,err,message) !output implicit none @@ -265,7 +335,7 @@ subroutine compcalday(julday, & !input real(dp) :: remainder ! remainder of modulus operation ! initialize errors - err=0; message="compcalday" + err=0; message="compCalday/" if(julday<=0)then;err=10;message=trim(message)//"no negative julian days/"; return; end if ! step 1 @@ -305,13 +375,13 @@ subroutine compcalday(julday, & !input remainder = remainder*min_per_hour - imin dsec = nint(remainder*secprmin) - end subroutine compcalday + end subroutine compCalday ! *************************************************************************************** ! public subroutine compgregcal_noleap: compute yy,mm,dd,hr,min,hr from a noleap julian day ! source: https://github.com/nmizukami/VIC/blob/VIC.5.0.0/vic/drivers/shared_all/src/vic_time.c ! *************************************************************************************** - subroutine compcalday_noleap(julday, & !input + subroutine compCalday_noleap(julday, & !input iyyy,mm,id,ih,imin,dsec,err,message) !output implicit none @@ -336,7 +406,7 @@ subroutine compcalday_noleap(julday, & !input real(dp) :: remainder ! remainder of modulus operation ! initialize errors - err=0; message="compcalday_noleap" + err=0; message="compCalday_noleap/" if(julday<=0)then;err=10;message=trim(message)//"no negative julian days/"; return; end if A = floor(julday+0.5_dp) @@ -383,7 +453,7 @@ subroutine compcalday_noleap(julday, & !input remainder = remainder*min_per_hour - imin dsec = nint(remainder*secprmin) - end subroutine compcalday_noleap + end subroutine compCalday_noleap ! *************************************************************************************** ! public function elapsedSec: calculate difference of two time marks obtained by date_and_time() diff --git a/route/build/src/var_lookup.f90 b/route/build/src/var_lookup.f90 index 98aafd86..e6354851 100644 --- a/route/build/src/var_lookup.f90 +++ b/route/build/src/var_lookup.f90 @@ -80,6 +80,8 @@ MODULE var_lookup integer(i4b) :: basArea = integerMissing ! area of the local HRUs contributing to each reach (m2) integer(i4b) :: upsArea = integerMissing ! area above the top of the reach -- zero if headwater (m2) integer(i4b) :: totalArea = integerMissing ! basArea + upsArea -- area at the bottom of the reach (m2) + ! abstraction/injection from reach + integer(i4b) :: QTAKE = integerMissing ! abstraction(-)/injection(+) coefficient [m3/s] ! lakes integer(i4b) :: basUnderLake = integerMissing ! Area of basin under lake (m2) integer(i4b) :: rchUnderLake = integerMissing ! Length of reach under lake (m) @@ -152,6 +154,7 @@ MODULE var_lookup !IRF state/fluxes type, public :: iLook_IRF integer(i4b) :: qfuture = integerMissing ! future routed flow + integer(i4b) :: irfVol = integerMissing ! reach volume endtype iLook_IRF ! *********************************************************************************************************** ! ** define data vectors @@ -162,12 +165,12 @@ MODULE var_lookup type(iLook_qDims) ,public,parameter :: ixqDims = iLook_qDims (1,2,3,4) type(iLook_HRU) ,public,parameter :: ixHRU = iLook_HRU (1) type(iLook_HRU2SEG) ,public,parameter :: ixHRU2SEG = iLook_HRU2SEG (1,2,3,4) - type(iLook_SEG) ,public,parameter :: ixSEG = iLook_SEG (1,2,3,4,5,6,7,8,9,10,11,12,13) + type(iLook_SEG) ,public,parameter :: ixSEG = iLook_SEG (1,2,3,4,5,6,7,8,9,10,11,12,13,14) type(iLook_NTOPO) ,public,parameter :: ixNTOPO = iLook_NTOPO (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17) type(iLook_PFAF) ,public,parameter :: ixPFAF = iLook_PFAF (1) type(iLook_RFLX) ,public,parameter :: ixRFLX = iLook_RFLX (1,2,3,4,5,6) type(iLook_KWT) ,public,parameter :: ixKWT = iLook_KWT (1,2,3,4,5) - type(iLook_IRF) ,public,parameter :: ixIRF = iLook_IRF (1) + type(iLook_IRF) ,public,parameter :: ixIRF = iLook_IRF (1,2) type(iLook_IRFbas ) ,public,parameter :: ixIRFbas = iLook_IRFbas (1,2) ! *********************************************************************************************************** ! ** define size of data vectors diff --git a/route/build/src/write_restart.f90 b/route/build/src/write_restart.f90 index 43bf9096..87da7434 100644 --- a/route/build/src/write_restart.f90 +++ b/route/build/src/write_restart.f90 @@ -1,8 +1,9 @@ MODULE write_restart ! Moudle wide external modules -USE nrtype, ONLY: i4b, dp, strLen +USE nrtype, ONLY: i4b, dp, lgt, strLen USE public_var +USE date_time, ONLY: datetime USE io_netcdf, ONLY: ncd_int USE io_netcdf, ONLY: ncd_float, ncd_double USE io_netcdf, ONLY: ncd_unlimited @@ -16,66 +17,187 @@ MODULE write_restart implicit none +integer(i4b), parameter :: currTimeStep = 1 +integer(i4b), parameter :: nextTimeStep = 2 + private -public::output_state +public::main_restart CONTAINS - SUBROUTINE output_state(ierr, message) + ! ********************************************************************* + ! public subroutine: restart write main routine + ! ********************************************************************* + SUBROUTINE main_restart(ierr, message) + + USE globalData, ONLY: restartAlarm ! logical to make alarm for restart writing + + implicit none + ! output variables + integer(i4b), intent(out) :: ierr ! error code + character(*), intent(out) :: message ! error message + ! local variables + character(len=strLen) :: cmessage ! error message of downwind routine + + ierr=0; message='main_restart/' + + call restart_alarm(ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + + if (restartAlarm) then + call restart_output(ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + end if + + END SUBROUTINE main_restart + + + ! ********************************************************************* + ! private subroutine: restart alarming + ! ********************************************************************* + SUBROUTINE restart_alarm(ierr, message) + + USE public_var, ONLY: calendar + USE public_var, ONLY: restart_write ! restart write options + USE public_var, ONLY: restart_day + USE globalData, ONLY: restartAlarm ! logical to make alarm for restart writing + USE globalData, ONLY: restCal ! restart Calendar time + USE globalData, ONLY: dropCal ! restart drop off Calendar time + USE globalData, ONLY: modTime ! previous and current model time + + implicit none + + ! output + integer(i4b), intent(out) :: ierr ! error code + character(*), intent(out) :: message ! error message + ! local variables + character(len=strLen) :: cmessage ! error message of downwind routine + integer(i4b) :: nDays ! number of days in a month + + ierr=0; message='restart_alarm/' + + ! adjust restart dropoff day if the dropoff day is outside number of days in particular month + call dropCal%set_datetime(dropCal%year(), dropCal%month(), restart_day, dropCal%hour(), dropCal%minute(), dropCal%sec()) + nDays = modTime(1)%ndays_month(calendar, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + if (dropCal%day() > nDays) then + call dropCal%set_datetime(dropCal%year(), dropCal%month(), nDays, dropCal%hour(), dropCal%minute(), dropCal%sec()) + end if + + ! adjust dropoff day further if restart day is actually outside number of days in a particular month + if (restCal%day() > nDays) then + dropCal = dropCal%add_day(-1, calendar, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + end if + + select case(trim(restart_write)) + case('Specified','specified','Last','last') + restartAlarm = (dropCal==modTime(1)) + case('Annual','annual') + restartAlarm = (dropCal%is_equal_mon(modTime(1)) .and. dropCal%is_equal_day(modTime(1)) .and. dropCal%is_equal_time(modTime(1))) + case('Monthly','monthly') + restartAlarm = (dropCal%is_equal_day(modTime(1)) .and. dropCal%is_equal_time(modTime(1))) + case('Daily','daily') + restartAlarm = dropCal%is_equal_time(modTime(1)) + case('Never','never') + restartAlarm = .false. + case default + ierr=20; message=trim(message)//'Current accepted options: L[l]ast, N[n]ever, S[s]pecified, Annual, Monthly, or Daily '; return + end select + + END SUBROUTINE restart_alarm + + + ! ********************************************************************* + ! private subroutine: write restart netCDF + ! ********************************************************************* + SUBROUTINE restart_output(ierr, message) - ! Saved Data - USE public_var, ONLY: output_dir - USE public_var, ONLY: case_name ! simulation name ==> output filename head USE public_var, ONLY: routOpt USE public_var, ONLY: time_units USE public_var, ONLY: dt - USE globalData, ONLY: runoff_data ! runoff data for one time step for LSM HRUs and River network HRUs + USE globalData, ONLY: runoff_data ! runoff data for one time step for LSM HRUs and River network HRUs USE globalData, ONLY: TSEC USE globalData, ONLY: reachID - USE globalData, ONLY: modTime ! previous and current model time - USE globalData, ONLY: modJulday ! current model Julian day - USE globalData, ONLY: restartJulday ! restart Julian day implicit none + ! output variables integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! local variables real(dp) :: TSEC1, TSEC2 character(len=strLen) :: cmessage ! error message of downwind routine - integer(i4b) :: sec_in_day ! second within day - character(len=strLen) :: fileout_state ! name of the output file - character(len=50),parameter :: fmtYMDS='(a,I0.4,a,I0.2,a,I0.2,a,I0.5,a)' - character(len=50),parameter :: fmtYMDHMS='(2a,I0.4,a,I0.2,a,I0.2,x,I0.2,a,I0.2,a,I0.2)' + character(len=strLen) :: fnameRestart ! name of the restart file name - if (abs(restartJulday-modJulday) output filename head + USE public_var, ONLY: calendar + USE public_var, ONLY: secprday + USE public_var, ONLY: dt + USE globalData, ONLY: modTime ! current model datetime + + implicit none + + ! input + integer(i4b), intent(in) :: timeStamp ! optional: + ! output + character(*), intent(out) :: fnameRestart ! name of the restart file name + integer(i4b), intent(out) :: ierr ! error code + character(*), intent(out) :: message ! error message + ! local variables + character(len=strLen) :: cmessage ! error message of downwind routine + type(datetime) :: timeStampCal ! datetime corresponding to file name time stamp + integer(i4b) :: sec_in_day ! second within day + character(len=50),parameter :: fmtYMDS='(a,I0.4,a,I0.2,a,I0.2,a,I0.5,a)' + character(len=50),parameter :: fmtYMDHMS='(2a,I0.4,a,I0.2,a,I0.2,x,I0.2,a,I0.2,a,I0.2)' + + ierr=0; message='restart_fname/' + + select case(timeStamp) + case(currTimeStep); timeStampCal = modTime(1) + case(nextTimeStep); timeStampCal = modTime(1)%add_sec(dt, calendar, ierr, cmessage) + case default; ierr=20; message=trim(message)//'time stamp option in restart filename: invalid -> 1: current time Step or 2: next time step'; return + end select + + write(iulog,fmtYMDHMS) new_line('a'),'Write restart file for ', & + timeStampCal%year(),'-',timeStampCal%month(),'-',timeStampCal%day(),timeStampCal%hour(),':',timeStampCal%minute(),':',nint(timeStampCal%sec()) + + sec_in_day = timeStampCal%hour()*60*60+timeStampCal%minute()*60+nint(timeStampCal%sec()) + + write(fnameRestart, fmtYMDS) trim(restart_dir)//trim(case_name)//'.r.', & + timeStampCal%year(), '-', timeStampCal%month(), '-', timeStampCal%day(), '-',sec_in_day,'.nc' + + END SUBROUTINE restart_fname - END SUBROUTINE output_state ! ********************************************************************* ! subroutine: define restart NetCDF file @@ -192,7 +314,7 @@ SUBROUTINE set_dim_len(ixDim, ierr, message1) case(ixStateDims%ens); meta_stateDims(ixStateDims%ens)%dimLength = 1 case(ixStateDims%tbound); meta_stateDims(ixStateDims%tbound)%dimLength = 2 case(ixStateDims%tdh); meta_stateDims(ixStateDims%tdh)%dimLength = size(FRAC_FUTURE) - case(ixStateDims%tdh_irf); meta_stateDims(ixStateDims%tdh_irf)%dimLength = 20 !just temporarily + case(ixStateDims%tdh_irf); meta_stateDims(ixStateDims%tdh_irf)%dimLength = 50 !just temporarily case(ixStateDims%wave); meta_stateDims(ixStateDims%wave)%dimLength = MAXQPAR case default; ierr=20; message1=trim(message1)//'unable to identify dimension variable index'; return end select @@ -618,6 +740,7 @@ SUBROUTINE write_IRF_state(ierr, message1) do iVar=1,nVarsIRF select case(iVar) case(ixIRF%qfuture); allocate(state(impulseResponseFunc)%var(iVar)%array_3d_dp(nSeg, ntdh_irf, nens), stat=ierr) + case(ixIRF%irfVol); allocate(state(impulseResponseFunc)%var(iVar)%array_2d_dp(nSeg, nens), stat=ierr) case default; ierr=20; message1=trim(message1)//'unable to identify variable index'; return end select if(ierr/=0)then; message1=trim(message1)//'problem allocating space for IRF routing state '//trim(meta_irf(iVar)%varName); return; endif @@ -635,6 +758,8 @@ SUBROUTINE write_IRF_state(ierr, message1) case(ixIRF%qfuture) state(impulseResponseFunc)%var(iVar)%array_3d_dp(iSeg,1:numQF(iens,iSeg),iens) = RCHFLX(iens,iSeg)%QFUTURE_IRF state(impulseResponseFunc)%var(iVar)%array_3d_dp(iSeg,numQF(iens,iSeg)+1:ntdh_irf,iens) = realMissing + case(ixIRF%irfVol) + state(impulseResponseFunc)%var(iVar)%array_2d_dp(iSeg,iens) = RCHFLX(iens,iSeg)%REACH_VOL(1) case default; ierr=20; message1=trim(message1)//'unable to identify variable index'; return end select @@ -651,6 +776,8 @@ SUBROUTINE write_IRF_state(ierr, message1) select case(iVar) case(ixIRF%qfuture) call write_nc(ncid, trim(meta_irf(iVar)%varName), state(impulseResponseFunc)%var(iVar)%array_3d_dp, (/1,1,1,iTime/), (/nSeg,ntdh_irf,nens,1/), ierr, cmessage) + case(ixIRF%irfVol) + call write_nc(ncid, trim(meta_irf(iVar)%varName), state(impulseResponseFunc)%var(iVar)%array_2d_dp, (/1,1,iTime/), (/nSeg,nens,1/), ierr, cmessage) case default; ierr=20; message1=trim(message1)//'unable to identify IRF variable index for nc writing'; return if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif end select diff --git a/route/build/src/write_simoutput.f90 b/route/build/src/write_simoutput.f90 index d2afea5c..bb2d6d1b 100644 --- a/route/build/src/write_simoutput.f90 +++ b/route/build/src/write_simoutput.f90 @@ -117,12 +117,8 @@ SUBROUTINE prep_output(ierr, message) ! out: error control USE public_var, only : time_units ! time units (seconds, hours, or days) ! saved global data USE globalData, only : basinID,reachID ! HRU and reach ID in network - USE globalData, only : modJulday ! julian day: at model time step USE globalData, only : modTime ! previous and current model time USE globalData, only : nEns, nHRU, nRch ! number of ensembles, HRUs and river reaches - ! subroutines - USE time_utils_module, only : compCalday ! compute calendar day - USE time_utils_module, only : compCalday_noleap ! compute calendar day implicit none @@ -136,37 +132,22 @@ SUBROUTINE prep_output(ierr, message) ! out: error control logical(lgt) :: defnewoutputfile ! flag to define new output file character(len=50),parameter :: fmtYMDS='(a,I0.4,a,I0.2,a,I0.2,a,I0.5,a)' - ! initialize error control ierr=0; message='prep_output/' - ! get the time - select case(trim(calendar)) - case('noleap') - call compCalday_noleap(modJulday,modTime(1)%iy,modTime(1)%im,modTime(1)%id,modTime(1)%ih,modTime(1)%imin,modTime(1)%dsec,ierr,cmessage) - case ('standard','gregorian','proleptic_gregorian') - call compCalday(modJulday,modTime(1)%iy,modTime(1)%im,modTime(1)%id,modTime(1)%ih,modTime(1)%imin,modTime(1)%dsec,ierr,cmessage) - case default; ierr=20; message=trim(message)//'calendar name: '//trim(calendar)//' invalid'; return - end select - if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - - ! print progress - write(iulog,'(a,I4,4(x,I4))') new_line('a'), modTime(1)%iy, modTime(1)%im, modTime(1)%id, modTime(1)%ih, modTime(1)%imin - - ! ***** - ! *** Define model output file... - ! ******************************* + ! print progress + write(iulog,'(a,I4,4(x,I4))') new_line('a'), modTime(1)%year(), modTime(1)%month(), modTime(1)%day(), modTime(1)%hour(), modTime(1)%minute() - ! check need for the new file - select case(trim(newFileFrequency)) - case('single'); defNewOutputFile=(modTime(0)%iy==integerMissing) - case('annual'); defNewOutputFile=(modTime(1)%iy/=modTime(0)%iy) - case('month'); defNewOutputFile=(modTime(1)%im/=modTime(0)%im) - case('day'); defNewOutputFile=(modTime(1)%id/=modTime(0)%id) + ! check need for the new file + select case(trim(newFileFrequency)) + case('single'); defNewOutputFile=(modTime(0)%year() ==integerMissing) + case('annual'); defNewOutputFile=(modTime(1)%year() /=modTime(0)%year()) + case('month'); defNewOutputFile=(modTime(1)%month()/=modTime(0)%month()) + case('day'); defNewOutputFile=(modTime(1)%day() /=modTime(0)%day()) case default; ierr=20; message=trim(message)//'unable to identify the option to define new output files'; return - end select + end select - ! define new file - if(defNewOutputFile)then + ! define new file + if(defNewOutputFile)then if (simout_nc%status == 2) then call close_nc(simout_nc%ncid, ierr, cmessage) @@ -180,11 +161,10 @@ SUBROUTINE prep_output(ierr, message) ! out: error control jTime=1 ! update filename - sec_in_day = modTime(1)%ih*60*60+modTime(1)%imin*60+nint(modTime(1)%dsec) - write(simout_nc%ncname, fmtYMDS) trim(output_dir)//trim(case_name)//'.mizuRoute.h.', & - modTime(1)%iy, '-', modTime(1)%im, '-', modTime(1)%id, '-',sec_in_day,'.nc' + sec_in_day = modTime(1)%hour()*60*60+modTime(1)%minute()*60+nint(modTime(1)%sec()) + write(simout_nc%ncname, fmtYMDS) trim(output_dir)//trim(case_name)//'.h.', & + modTime(1)%year(), '-', modTime(1)%month(), '-', modTime(1)%day(), '-',sec_in_day,'.nc' - ! define output file call defineFile(simout_nc%ncname, & ! input: file name nEns, & ! input: number of ensembles nHRU, & ! input: number of HRUs @@ -199,22 +179,18 @@ SUBROUTINE prep_output(ierr, message) ! out: error control simout_nc%status = 2 - ! define basin ID - call write_nc(simout_nc%ncid, 'basinID', basinID, (/1/), (/nHRU/), ierr, cmessage) + call write_nc(simout_nc%ncid, 'basinID', int(basinID,kind(i4b)), (/1/), (/nHRU/), ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - ! define reach ID call write_nc(simout_nc%ncid, 'reachID', reachID, (/1/), (/nRch/), ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif ! no new file requested: increment time - else + else jTime = jTime+1 - endif - - modTime(0) = modTime(1) + endif END SUBROUTINE prep_output @@ -324,15 +300,12 @@ SUBROUTINE defineFile(fname, & ! input: filename end do - ! put global attribute call put_global_attr(ncid, 'version', trim(mizuRouteVersion) ,ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - ! end definitions call end_def(ncid, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - ! close NetCDF file call close_nc(ncid, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif