diff --git a/route/build/Makefile b/route/build/Makefile index 4d928ee1..791b2dbd 100644 --- a/route/build/Makefile +++ b/route/build/Makefile @@ -148,8 +148,9 @@ DATATYPES = \ datetime_data.f90 \ csv_data.f90 \ globalData.f90 \ - historyFile.f90 \ popMetadat.f90 \ + histVars_data.f90 \ + historyFile.f90 \ allocation.f90 # define utilities UTILS = \ @@ -177,8 +178,8 @@ IO = \ standalone/read_remap.f90 \ io_rpointfile.f90 \ read_restart.f90 \ - write_restart_pio.f90 \ - write_simoutput_pio.f90 + write_simoutput_pio.f90 \ + write_restart_pio.f90 # CORE CORE = \ accum_runoff.f90 \ diff --git a/route/build/cpl/RtmTimeManager.F90 b/route/build/cpl/RtmTimeManager.F90 index 02a06e9e..028d8fe2 100644 --- a/route/build/cpl/RtmTimeManager.F90 +++ b/route/build/cpl/RtmTimeManager.F90 @@ -120,6 +120,7 @@ SUBROUTINE init_time(ierr, message) iTime = 1 simDatetime(0) = datetime(integerMissing, integerMissing, integerMissing, integerMissing, integerMissing, realMissing) simDatetime(1) = begDatetime + simDatetime(2) = simDatetime(1)%add_sec(dt, calendar, ierr, cmessage) if (masterproc .and. debug_write) then write(iulog,*) 'simStart datetime = ', trim(simStart) diff --git a/route/build/src/globalData.f90 b/route/build/src/globalData.f90 index 198498c8..54322969 100644 --- a/route/build/src/globalData.f90 +++ b/route/build/src/globalData.f90 @@ -89,7 +89,7 @@ MODULE globalData integer(i4b), public :: iTime ! time index at simulation time step real(dp), public :: timeVar ! time variables (unit given by time variable) real(dp), public :: TSEC(0:1) ! begning and end of time step since simulation started (sec) - type(datetime), public :: simDatetime(0:1) ! previous and current simulation time (yyyy:mm:dd:hh:mm:ss) + type(datetime), public :: simDatetime(0:2) ! previous, current and next simulation time (yyyy:mm:dd:hh:mm:ss) type(datetime), public :: begDatetime ! simulation start date/time (yyyy:mm:dd:hh:mm:ss) type(datetime), public :: endDatetime ! simulation end date/time (yyyy:mm:dd:hh:mm:ss) type(datetime), public :: restDatetime ! desired restart date/time (yyyy:mm:dd:hh:mm:ss) @@ -104,13 +104,13 @@ MODULE globalData ! ---------- Misc. data ------------------------------------------------------------------------- character(len=strLen), public :: runMode='standalone' ! run options: standalone or cesm-coupling - logical(lgt), public :: isHistFileOpen=.false. ! flag to indicate history output netcdf is open integer(i4b), public :: ixPrint(1:2)=integerMissing ! index of desired reach to be on-screen print integer(i4b), public :: nEns=1 ! number of ensemble + integer(i4b), public :: maxtdh ! maximum unit-hydrograph future time steps type(cMolecule), public :: nMolecule ! number of computational molecule (used for KW, MC, DW) - character(300), public :: hfileout=charMissing ! name of the history output file - character(300), public :: hfileout_gage=charMissing ! name of the gage-only history output file - character(300), public :: rfileout=charMissing ! name of the restart output file + character(300), public :: hfileout=charMissing ! history output file name + character(300), public :: hfileout_gage=charMissing ! gage-only history output file name + character(300), public :: rfileout=charMissing ! restart output file name ! ---------- MPI/OMP/PIO variables ---------------------------------------------------------------- diff --git a/route/build/src/histVars_data.f90 b/route/build/src/histVars_data.f90 new file mode 100644 index 00000000..d28b99bc --- /dev/null +++ b/route/build/src/histVars_data.f90 @@ -0,0 +1,263 @@ +MODULE histVars_data + + ! history file variable data + + USE nrtype + USE dataTypes, ONLY: STRFLX + USE var_lookup, ONLY: ixRFLX, nVarsRFLX + USE var_lookup, ONLY: ixHFLX, nVarsHFLX + USE public_var, ONLY: accumRunoff ! accum runoff ID = 0 + USE public_var, ONLY: impulseResponseFunc ! IRF routing ID = 1 + USE public_var, ONLY: kinematicWaveTracking ! KWT routing ID = 2 + USE public_var, ONLY: kinematicWave ! KW routing ID = 3 + USE public_var, ONLY: muskingumCunge ! MC routing ID = 4 + USE public_var, ONLY: diffusiveWave ! DW routing ID = 5 + USE globalData, ONLY: nRoutes + USE globalData, ONLY: routeMethods + USE globalData, ONLY: meta_rflx, meta_hflx + USE globalData, ONLY: idxSUM,idxIRF,idxKWT, & + idxKW,idxMC,idxDW + + implicit none + + public:: histVars + + type :: histVars + ! -------- + ! output varialbes: basRunoff, instRunoff, dlayRunoff, discharge, elev, volume + ! if new variables need to be written, need to add here AND each procedure + ! -------- + integer(i4b) :: nt ! number of time steps over variables are aggregated + integer(i4b) :: nHru ! number of + integer(i4b) :: nRch ! number of + real(dp) :: timeVar ! time variable [unit] at nt = 1 + real(sp), allocatable :: basRunoff(:) ! catchment average runoff [m/s] [nHru] + real(sp), allocatable :: instRunoff(:) ! instantaneous lateral inflow into each river/lake [m3/s] [nRch] + real(sp), allocatable :: dlayRunoff(:) ! lateral inflow into river or lake [m3/s] for each reach [nRch] + real(sp), allocatable :: discharge(:,:) ! river/lake discharge [m3/s] for each reach/lake and routing method [nRch,nMethod] + real(sp), allocatable :: elev(:,:) ! river/lake surface water elevation [m] for each reach/lake and routing method [nRch,nMethod] + real(sp), allocatable :: volume(:,:) ! river/lake volume [m3] for each reach/lake and routing method [nRch,nMethod] + + CONTAINS + + procedure, public :: aggregate ! Accumulating output variables + procedure, public :: finalize ! Compute aggregated values (currently only mean) to be written in netCDF + procedure, public :: refresh ! Reset output arrays to zero + procedure, public :: clean ! Deallocate all the output array memories + + end type histVars + + INTERFACE histVars + module procedure constructor + END INTERFACE histVars + + CONTAINS + + ! ----------------------------------------------------- + ! constructor - instantiate history output data structure + ! ----------------------------------------------------- + FUNCTION constructor(nHru_local, nRch_local, ierr, message) RESULT(instHistVar) + + implicit none + ! argument variables + integer(i4b), intent(in) :: nHru_local ! number of hru in each proc + integer(i4b), intent(in) :: nRch_local ! number of Rch in each proc + integer(i4b), intent(out) :: ierr ! error code + character(*), intent(out) :: message ! error message + type(histVars) :: instHistVar + ! local variables + character(strLen) :: cmessage ! error message from subroutine + + ierr=0; message='initHistVars/' + + instHistVar%nt = 0 + instHistVar%nHru = nHRU_local + instHistVar%nRch = nRch_local + + if (meta_hflx(ixHFLX%basRunoff)%varFile) then + allocate(instHistVar%basRunoff(nHRU_local), stat=ierr, errmsg=cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [instHistVar%basRunoff]'; return; endif + instHistVar%basRunoff(1:nHRU_local) = 0._sp + end if + + if (meta_rflx(ixRFLX%instRunoff)%varFile) then + allocate(instHistVar%instRunoff(nRch_local), stat=ierr, errmsg=cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [instHistVar%instRunoff]'; return; endif + instHistVar%instRunoff(1:nRch_local) = 0._sp + end if + + if (meta_rflx(ixRFLX%dlayRunoff)%varFile) then + allocate(instHistVar%dlayRunoff(nRch_local), stat=ierr, errmsg=cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [instHistVar%instRunoff]'; return; endif + instHistVar%dlayRunoff(1:nRch_local) = 0._sp + end if + + if (nRoutes>0) then ! this should be number of methods that ouput + allocate(instHistVar%discharge(nRch_local, nRoutes), stat=ierr, errmsg=cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [instHistVar%discharge]'; return; endif + instHistVar%discharge(1:nRch_local, 1:nRoutes) = 0._sp + + allocate(instHistVar%volume(nRch_local, nRoutes), stat=ierr, errmsg=cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [instHistVar%volume]'; return; endif + instHistVar%volume(1:nRch_local, 1:nRoutes) = 0._sp + end if + + END FUNCTION constructor + + ! --------------------------------------------------------------- + ! accumulate data + ! --------------------------------------------------------------- + SUBROUTINE aggregate(this, & ! inout: + timeVar_local, & ! input: time variables current + basRunoff_local, & ! input: + RCHFLX_local, & ! input: + ierr, message) ! output: error handling + + implicit none + ! argument variables + class(histVars), intent(inout) :: this + real(dp), intent(in) :: timeVar_local + real(dp), intent(in) :: basRunoff_local(:) + type(STRFLX), intent(in) :: RCHFLX_local(:,:) + integer(i4b), intent(out) :: ierr ! error code + character(*), intent(out) :: message ! error message + ! local variables + integer(i4b) :: nHRU_input + integer(i4b) :: nRch_input + integer(i4b) :: ix, iRoute ! loop indices + integer(i4b) :: idxMethod ! temporal method index + + ierr=0; message='aggregate/' + + ! -- increment number of sample + this%nt = this%nt + 1 + + if (this%nt == 1) then + this%timeVar = timeVar_local + end if + + ! -- array size checks - input data vs history output buffer + ! hru and reach size in input data + nHru_input = size(basRunoff_local) + nRch_input = size(RCHFLX_local(1,:)) + if (nHru_input/=this%nHru) then + write(message,'(2A,I,A,I)') trim(message),'history buffer hru size:',this%nHru,'/= input data hru size:',nHru_input + ierr=81; return + end if + if (nRch_input/=this%nRch) then + write(message,'(2A,I,A,I)') trim(message),'history buffer reach size:',this%nRch,'/= input data reach size:',nRch_input + ierr=81; return + end if + + ! ---- aggregate + ! 1. basin runoff + if (meta_hflx(ixHFLX%basRunoff)%varFile) then + this%basRunoff(1:this%nHRU) = this%basRunoff(1:this%nHRU) + basRunoff_local(1:this%nHru) + end if + + ! 2. instantaneous runoff into reach + if (meta_rflx(ixRFLX%instRunoff)%varFile) then + this%instRunoff(1:this%nRch) = this%instRunoff(1:this%nRch) + RCHFLX_local(1,1:this%nRch)%BASIN_QI + end if + + ! 3. delayed runoff into reach + if (meta_rflx(ixRFLX%dlayRunoff)%varFile) then + this%dlayRunoff(1:this%nRch) = this%dlayRunoff(1:this%nRch) + RCHFLX_local(1,1:this%nRch)%BASIN_QR(1) + end if + + ! 4. discharge and volume + do iRoute = 1, nRoutes + select case(routeMethods(iRoute)) + case(accumRunoff); idxMethod=idxSUM + case(impulseResponseFunc); idxMethod=idxIRF + case(kinematicWaveTracking); idxMethod=idxKWT + case(kinematicWave); idxMethod=idxKW + case(muskingumCunge); idxMethod=idxMC + case(diffusiveWave); idxMethod=idxDW + case default + write(message,'(2A,X,I,X,A)') trim(message), 'routing method index:',routeMethods(iRoute), 'must be 0-5' + ierr=81; return + end select + + do ix=1,this%nRch + this%discharge(ix,iRoute) = this%discharge(ix,iRoute) + RCHFLX_local(1,ix)%ROUTE(idxMethod)%REACH_Q + this%volume(ix,iRoute) = this%volume(ix,iRoute) + RCHFLX_local(1,ix)%ROUTE(idxMethod)%REACH_VOL(1) + end do + end do + + END SUBROUTINE aggregate + + ! --------------------------------------------------------------- + ! finalize output variable (compute mean) + ! --------------------------------------------------------------- + SUBROUTINE finalize(this) + + implicit none + ! Argument variables + class(histVars), intent(inout) :: this + + ! ---- aggregate + ! 1. basin runoff + if (allocated(this%basRunoff)) then + this%basRunoff = this%basRunoff/real(this%nt, kind=sp) + end if + + ! 2. instantaneous runoff into reach + if (allocated(this%instRunoff)) then + this%instRunoff = this%instRunoff/real(this%nt, kind=sp) + end if + + ! 3. delayed runoff into reach + if (allocated(this%dlayRunoff)) then + this%dlayRunoff = this%dlayRunoff/real(this%nt, kind=sp) + end if + + ! 4. discharge + if (allocated(this%discharge)) then + this%discharge = this%discharge/real(this%nt, kind=sp) + end if + + ! 5. volume + if (allocated(this%volume)) then + this%volume = this%volume/real(this%nt, kind=sp) + end if + + END SUBROUTINE finalize + + ! --------------------------------------------------------------- + ! re-initialzae intantiated data structure + ! --------------------------------------------------------------- + SUBROUTINE refresh(this) + + implicit none + ! Argument variables + class(histVars), intent(inout) :: this + + this%nt = 0 + + if (allocated(this%basRunoff)) this%basRunoff = 0._sp + if (allocated(this%instRunoff)) this%instRunoff = 0._sp + if (allocated(this%dlayRunoff)) this%dlayRunoff = 0._sp + if (allocated(this%discharge)) this%discharge = 0._sp + if (allocated(this%volume)) this%volume = 0._sp + + END SUBROUTINE refresh + + ! --------------------------------------------------------------- + ! Release memory (deallocate array) + ! --------------------------------------------------------------- + SUBROUTINE clean(this) + + implicit none + ! Argument variables + class(histVars), intent(inout) :: this + + if (allocated(this%basRunoff)) deallocate(this%basRunoff) + if (allocated(this%instRunoff)) deallocate(this%instRunoff) + if (allocated(this%dlayRunoff)) deallocate(this%dlayRunoff) + if (allocated(this%discharge)) deallocate(this%discharge) + if (allocated(this%volume)) deallocate(this%volume) + + END SUBROUTINE clean + +END MODULE histVars_data diff --git a/route/build/src/historyFile.f90 b/route/build/src/historyFile.f90 index 600ef934..c0a0544b 100644 --- a/route/build/src/historyFile.f90 +++ b/route/build/src/historyFile.f90 @@ -1,7 +1,6 @@ MODULE historyFile USE nrtype - USE dataTypes, ONLY: STRFLX USE var_lookup, ONLY: ixRFLX, nVarsRFLX USE var_lookup, ONLY: ixHFLX, nVarsHFLX USE globalData, ONLY: idxSUM,idxIRF,idxKWT, & @@ -17,6 +16,7 @@ MODULE historyFile USE globalData, ONLY: pio_rearranger USE globalData, ONLY: pio_root USE globalData, ONLY: pio_stride + USE histVars_data, ONLY: histVars USE pio_utils USE perf_mod, ONLY: t_startf,t_stopf ! timing start/stop (GPTL library) @@ -106,7 +106,7 @@ SUBROUTINE set_compdof_rch(this, compdof_rch, nRch_in) ! initialze iodescRchFlux call pio_decomp(this%pioSystem, & ! input: pio system descriptor - ncd_double, & ! input: data type (pio_int, pio_real, pio_double, pio_char) + ncd_float, & ! input: data type (pio_int, pio_real, pio_double, pio_char) [nRch_in], & ! input: dimension length == global array size compdof_rch, & ! input: local->global mapping this%ioDescRchFlux) @@ -127,13 +127,13 @@ SUBROUTINE set_compdof_rch_hru(this, compdof_rch, compdof_hru, nRch_in, nHRU_in) ! initialze iodescRchFlux and ioDescHruFlux call pio_decomp(this%pioSystem, & ! input: pio system descriptor - ncd_double, & ! input: data type (pio_int, pio_real, pio_double, pio_char) + ncd_float, & ! input: data type (pio_int, pio_real, pio_double, pio_char) [nRch_in], & ! input: dimension length == global array size compdof_rch, & ! input: local->global mapping this%ioDescRchFlux) call pio_decomp(this%pioSystem, & ! input: pio system descriptor - ncd_double, & ! input: data type (pio_int, pio_real, pio_double, pio_char) + ncd_float, & ! input: data type (pio_int, pio_real, pio_double, pio_char) [nHRU_in], & ! input: dimension length == global array size compdof_hru, & ! input: local->global mapping this%ioDescHruFlux) @@ -450,13 +450,12 @@ SUBROUTINE write_loc_rch(this, reach_id, ierr, message) ! --------------------------------- ! writing flux (hru and/or rch) ! --------------------------------- - SUBROUTINE write_flux(this, time, RCHFLX_local, index_write, ierr, message) + SUBROUTINE write_flux(this, hVars_local, index_write, ierr, message) implicit none ! Argument variables class(histFile), intent(inout) :: this - real(dp), intent(in) :: time - type(STRFLX), intent(in) :: RCHFLX_local(:,:) + type(histVars), intent(in) :: hVars_local integer(i4b), allocatable, intent(in) :: index_write(:) integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message @@ -468,17 +467,17 @@ SUBROUTINE write_flux(this, time, RCHFLX_local, index_write, ierr, message) this%iTime = this%iTime + 1 ! this is only line to increment time step index ! write time -- note time is just carried across from the input - call write_netcdf(this%pioFileDesc, 'time', [time], [this%iTime], [1], ierr, cmessage) + call write_netcdf(this%pioFileDesc, 'time', [hVars_local%timeVar], [this%iTime], [1], ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif if (.not.this%gageOutput) then if (meta_hflx(ixHFLX%basRunoff)%varFile) then - call this%write_flux_hru(ierr, cmessage) + call this%write_flux_hru(hVars_local, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif endif end if - call this%write_flux_rch(RCHFLX_local, index_write, ierr, cmessage) + call this%write_flux_rch(hVars_local, index_write, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif call sync_file(this%pioFileDesc, ierr, cmessage) @@ -487,42 +486,41 @@ SUBROUTINE write_flux(this, time, RCHFLX_local, index_write, ierr, message) END SUBROUTINE write_flux - SUBROUTINE write_flux_hru(this, ierr, message) - ! currently only basin Runoff to be output + SUBROUTINE write_flux_hru(this, hVars_local, ierr, message) + + ! Write HRU flux variables + ! currently only basin hru Runoff [m/s] + implicit none ! Argument variables class(histFile), intent(inout) :: this + type(histVars), intent(in) :: hVars_local integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! local variables - real(dp), allocatable :: basinRunoff(:) character(len=strLen) :: cmessage ! error message of downwind routine ierr=0; message='write_flux_hru/' if (meta_hflx(ixHFLX%basRunoff)%varFile) then - call get_proc_flux(ierr, cmessage, basinRunoff=basinRunoff) - - ! write the basin runoff at HRU (unit: the same as runoff input) - call write_pnetcdf_recdim(this%pioFileDesc, 'basRunoff',basinRunoff, this%ioDescHruFlux, this%iTime, ierr, cmessage) + call write_pnetcdf_recdim(this%pioFileDesc, 'basRunoff', hVars_local%basRunoff, this%ioDescHruFlux, this%iTime, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif endif END SUBROUTINE write_flux_hru - SUBROUTINE write_flux_rch(this, RCHFLX_local, index_write, ierr, message) + SUBROUTINE write_flux_rch(this, hVars_local, index_write, ierr, message) implicit none ! Argument variables class(histFile), intent(inout) :: this - type(STRFLX), intent(in) :: RCHFLX_local(:,:) + type(histVars), intent(in) :: hVars_local integer(i4b), allocatable, intent(in) :: index_write(:) integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! local variables - real(dp), allocatable :: array_temp(:) - integer(i4b) :: ix + real(sp), allocatable :: array_temp(:) integer(i4b) :: nRch_write character(len=strLen) :: cmessage ! error message of downwind routine @@ -544,7 +542,7 @@ SUBROUTINE write_flux_rch(this, RCHFLX_local, index_write, ierr, message) ! write instataneous local runoff in each stream segment (m3/s) if (meta_rflx(ixRFLX%instRunoff)%varFile) then if (nRch_write>0) then - array_temp(1:nRch_write) = RCHFLX_local(1,index_write)%BASIN_QI + array_temp(1:nRch_write) = hVars_local%instRunoff(index_write) end if call write_pnetcdf_recdim(this%pioFileDesc, 'instRunoff', array_temp, this%ioDescRchFlux, this%iTime, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif @@ -553,7 +551,7 @@ SUBROUTINE write_flux_rch(this, RCHFLX_local, index_write, ierr, message) ! write routed local runoff in each stream segment (m3/s) if (meta_rflx(ixRFLX%dlayRunoff)%varFile) then if (nRch_write>0) then - array_temp(1:nRch_write) = RCHFLX_local(1,index_write)%BASIN_QR(1) + array_temp(1:nRch_write) = hVars_local%dlayRunoff(index_write) end if call write_pnetcdf_recdim(this%pioFileDesc, 'dlayRunoff', array_temp, this%ioDescRchFlux, this%iTime, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif @@ -561,9 +559,7 @@ SUBROUTINE write_flux_rch(this, RCHFLX_local, index_write, ierr, message) if (meta_rflx(ixRFLX%sumUpstreamRunoff)%varFile) then if (nRch_write>0) then - do ix=1,nRch_write - array_temp(ix) = RCHFLX_local(1,index_write(ix))%ROUTE(idxSUM)%REACH_Q - end do + array_temp(1:nRch_write) = hVars_local%discharge(index_write, idxSUM) end if call write_pnetcdf_recdim(this%pioFileDesc, 'sumUpstreamRunoff', array_temp, this%ioDescRchFlux, this%iTime, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif @@ -571,9 +567,7 @@ SUBROUTINE write_flux_rch(this, RCHFLX_local, index_write, ierr, message) if (meta_rflx(ixRFLX%KWTroutedRunoff)%varFile) then if (nRch_write>0) then - do ix=1,nRch_write - array_temp(ix) = RCHFLX_local(1,index_write(ix))%ROUTE(idxKWT)%REACH_Q - end do + array_temp(1:nRch_write) = hVars_local%discharge(index_write, idxKWT) end if call write_pnetcdf_recdim(this%pioFileDesc, 'KWTroutedRunoff', array_temp, this%ioDescRchFlux, this%iTime, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif @@ -582,9 +576,7 @@ SUBROUTINE write_flux_rch(this, RCHFLX_local, index_write, ierr, message) call t_startf ('output/write_flux/irf') if (meta_rflx(ixRFLX%IRFroutedRunoff)%varFile) then if (nRch_write>0) then - do ix=1,nRch_write - array_temp(ix) = RCHFLX_local(1,index_write(ix))%ROUTE(idxIRF)%REACH_Q - end do + array_temp(1:nRch_write) = hVars_local%discharge(index_write, idxIRF) end if call write_pnetcdf_recdim(this%pioFileDesc, 'IRFroutedRunoff', array_temp, this%ioDescRchFlux, this%iTime, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif @@ -593,9 +585,7 @@ SUBROUTINE write_flux_rch(this, RCHFLX_local, index_write, ierr, message) if (meta_rflx(ixRFLX%KWroutedRunoff)%varFile) then if (nRch_write>0) then - do ix=1,nRch_write - array_temp(ix) = RCHFLX_local(1,index_write(ix))%ROUTE(idxKW)%REACH_Q - end do + array_temp(1:nRch_write) = hVars_local%discharge(index_write, idxKW) end if call write_pnetcdf_recdim(this%pioFileDesc, 'KWroutedRunoff', array_temp, this%ioDescRchFlux, this%iTime, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif @@ -603,9 +593,7 @@ SUBROUTINE write_flux_rch(this, RCHFLX_local, index_write, ierr, message) if (meta_rflx(ixRFLX%MCroutedRunoff)%varFile) then if (nRch_write>0) then - do ix=1,nRch_write - array_temp(ix) = RCHFLX_local(1,index_write(ix))%ROUTE(idxMC)%REACH_Q - end do + array_temp(1:nRch_write) = hVars_local%discharge(index_write, idxMC) end if call write_pnetcdf_recdim(this%pioFileDesc, 'MCroutedRunoff', array_temp, this%ioDescRchFlux, this%iTime, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif @@ -613,9 +601,7 @@ SUBROUTINE write_flux_rch(this, RCHFLX_local, index_write, ierr, message) if (meta_rflx(ixRFLX%DWroutedRunoff)%varFile) then if (nRch_write>0) then - do ix=1,nRch_write - array_temp(ix) = RCHFLX_local(1,index_write(ix))%ROUTE(idxDW)%REACH_Q - end do + array_temp(1:nRch_write) = hVars_local%discharge(index_write, idxDW) end if call write_pnetcdf_recdim(this%pioFileDesc, 'DWroutedRunoff', array_temp, this%ioDescRchFlux, this%iTime, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif @@ -623,9 +609,7 @@ SUBROUTINE write_flux_rch(this, RCHFLX_local, index_write, ierr, message) if (meta_rflx(ixRFLX%IRFvolume)%varFile) then if (nRch_write>0) then - do ix=1,nRch_write - array_temp(ix) = RCHFLX_local(1,index_write(ix))%ROUTE(idxIRF)%REACH_VOL(1) - end do + array_temp(1:nRch_write) = hVars_local%volume(index_write, idxIRF) end if call write_pnetcdf_recdim(this%pioFileDesc, 'IRFvolume', array_temp, this%ioDescRchFlux, this%iTime, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif @@ -633,40 +617,4 @@ SUBROUTINE write_flux_rch(this, RCHFLX_local, index_write, ierr, message) END SUBROUTINE write_flux_rch - ! --------------------------------- - ! non-type-bound subroutine - organize flux data array or structure into per processor - ! --------------------------------- - SUBROUTINE get_proc_flux(ierr, message, basinRunoff) - - USE globalData, ONLY: nHRU_mainstem ! number of mainstem HRUs - USE globalData, ONLY: basinRunoff_main ! mainstem only HRU runoff - USE globalData, ONLY: basinRunoff_trib ! tributary only HRU runoff - USE globalData, ONLY: hru_per_proc ! number of hrus assigned to each proc (size = num of procs+1) - - implicit none - ! Argument variables - real(dp), allocatable, optional, intent(out) :: basinRunoff(:) - integer(i4b), intent(out) :: ierr ! error code - character(*), intent(out) :: message ! error message - ! local variables - integer(i4b) :: nHRU_local - character(strLen) :: cmessage ! error message of downwind routine - - if (present(basinRunoff)) then - if (masterproc) then - nHRU_local = nHRU_mainstem + hru_per_proc(0) - allocate(basinRunoff(nHRU_local), stat=ierr, errmsg=cmessage) - if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [basinRunoff]'; return; endif - if (nHRU_mainstem>0) basinRunoff(1:nHRU_mainstem) = basinRunoff_main(1:nHRU_mainstem) - if (hru_per_proc(0)>0) basinRunoff(nHRU_mainstem+1:nHRU_local) = basinRunoff_trib(:) - else - nHRU_local = hru_per_proc(pid) - allocate(basinRunoff(nHRU_local), stat=ierr, errmsg=cmessage) - if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [basinRunoff]'; return; endif - basinRunoff = basinRunoff_trib - endif - end if - - END SUBROUTINE get_proc_flux - END MODULE historyFile diff --git a/route/build/src/init_model_data.f90 b/route/build/src/init_model_data.f90 index bd6c0422..20fe6681 100644 --- a/route/build/src/init_model_data.f90 +++ b/route/build/src/init_model_data.f90 @@ -312,6 +312,7 @@ SUBROUTINE update_time(finished, ierr, message) ! increment model datetime simDatetime(0) = simDatetime(1) simDatetime(1) = simDatetime(1)%add_sec(dt, calendar, ierr, cmessage) + simDatetime(2) = simDatetime(1)%add_sec(dt, calendar, ierr, cmessage) ! model time stamp variable for output - dt is in second t_unit = trim( time_units(1:index(time_units,' ')) ) diff --git a/route/build/src/mpi_process.f90 b/route/build/src/mpi_process.f90 index e9852eea..ee9d21dd 100644 --- a/route/build/src/mpi_process.f90 +++ b/route/build/src/mpi_process.f90 @@ -2951,6 +2951,7 @@ END SUBROUTINE subrch_struc2array SUBROUTINE pass_global_data(comm, ierr, message) ! output: error control USE globalData, ONLY: nRch,nHRU ! number of reaches and hrus in whole network USE globalData, ONLY: iTime ! time index + USE globalData, ONLY: maxtdh ! implicit none ! argument variables integer(i4b), intent(in) :: comm ! communicator @@ -2966,6 +2967,8 @@ SUBROUTINE pass_global_data(comm, ierr, message) ! output: error control call MPI_BCAST(calendar, strLen, MPI_CHARACTER, root, comm, ierr) call MPI_BCAST(time_units,strLen, MPI_CHARACTER, root, comm, ierr) + CALL MPI_ALLREDUCE(maxtdh, maxtdh, 1, MPI_INTEGER, MPI_MAX, comm, ierr) + END SUBROUTINE pass_global_data END MODULE mpi_process diff --git a/route/build/src/pio_utils.f90 b/route/build/src/pio_utils.f90 index d7b8251a..03e25ba7 100644 --- a/route/build/src/pio_utils.f90 +++ b/route/build/src/pio_utils.f90 @@ -60,6 +60,7 @@ MODULE pio_utils module procedure write_int_darray1D module procedure write_int_darray2D module procedure write_int_darray3D + module procedure write_float_darray1D module procedure write_double_darray1D module procedure write_double_darray2D module procedure write_double_darray3D @@ -784,6 +785,34 @@ SUBROUTINE write_int_darray3D(pioFileDesc, & END SUBROUTINE write_int_darray3D + ! --------------------------------------------------------------- + ! write distributed double vector into 1D data + SUBROUTINE write_float_darray1D(pioFileDesc, & + vname, & ! input: variable name + array, & ! input: variable data + iodesc, & ! input: ??? it is from initdecomp routine + ierr, message) ! output: error control + implicit none + ! Argument variables: + type(file_desc_t), intent(inout) :: pioFileDesc ! pio file handle + character(*), intent(in) :: vname ! variable name + real(sp), intent(in) :: array(:) ! variable data + type(io_desc_t), intent(inout) :: iodesc ! io descriptor handle that is generated in PIO_initdecomp + integer(i4b), intent(out) :: ierr ! error code + character(*), intent(out) :: message ! error message + ! local variables + type(var_desc_t) :: pioVarId ! netCDF variable ID + + ierr=0; message='write_float_darray1D/' + + ierr = pio_inq_varid(pioFileDesc, trim(vname), pioVarId) + if(ierr/=0)then; message=trim(message)//'ERROR: getting variable id'; return; endif + + call pio_write_darray(pioFileDesc, pioVarId, iodesc, array, ierr) + if(ierr/=pio_noerr)then; message=trim(message)//'cannot write data'; return; endif + + END SUBROUTINE write_float_darray1D + ! --------------------------------------------------------------- ! write distributed double vector into 1D data SUBROUTINE write_double_darray1D(pioFileDesc, & diff --git a/route/build/src/popMetadat.f90 b/route/build/src/popMetadat.f90 index c6f21e97..133f85b0 100644 --- a/route/build/src/popMetadat.f90 +++ b/route/build/src/popMetadat.f90 @@ -94,6 +94,7 @@ subroutine popMetadat(err,message) meta_dims (ixDims%pfaf ) = dim_info('pfaf' , integerMissing, integerMissing) ! max pfafstetter code length meta_stateDims(ixStateDims%seg ) = dim_info('seg', integerMissing, integerMissing) ! stream segment vector + meta_stateDims(ixStateDims%hru ) = dim_info('hru', integerMissing, integerMissing) ! catchment hru vector meta_stateDims(ixStateDims%time ) = dim_info('time', integerMissing, integerMissing) ! time meta_stateDims(ixStateDims%tbound ) = dim_info('tbound', integerMissing, 2) ! time bound (alway 2 - start and end) meta_stateDims(ixStateDims%ens ) = dim_info('ens', integerMissing, integerMissing) ! runoff ensemble diff --git a/route/build/src/process_param.f90 b/route/build/src/process_param.f90 index 6b658e4a..84032985 100644 --- a/route/build/src/process_param.f90 +++ b/route/build/src/process_param.f90 @@ -119,6 +119,7 @@ SUBROUTINE make_uh(& ! global variables USE public_var, ONLY: pi ! pi USE dataTypes, ONLY: dlength + USE globalData, ONLY: maxtdh ! maximum unit-hydrogrph future time implicit none ! Argument variables real(dp), intent(in) :: length(:) ! river segment length @@ -143,6 +144,7 @@ SUBROUTINE make_uh(& integer(i4b) :: iHrStrt ! index of UH time step where rising limb of UH start integer(i4b) :: iHrLast ! index of UH time step where recession limb of UH become zero integer(i4b) :: nTSub ! number of time steps where 1/nTsub [m] of runoff is inserted + integer(i4b) :: ntdh ! number of values on the time delay histogram integer(i4b) :: iSeg ! Loop index integer(i4b) :: iHr,jHr ! Loop index of hour integer(i4b) :: iTagg ! index for aggregated (i.e. simulation) time step @@ -157,6 +159,8 @@ SUBROUTINE make_uh(& !nTsub=floor(dt/dTUH) nSeg = size(length) + maxtdh = 0 + ! Memory allocation allocate(seg_uh(nSeg), stat=ierr, errmsg=cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage)//': seg_uh'; return; endif @@ -249,8 +253,10 @@ SUBROUTINE make_uh(& ! Re-normalize the UHQ by its sum UHQ = UHQ/INTE + ntdh = (iHrLast+nTsub-1)/nTsub + maxtdh = max(maxtdh, ntdh) !Aggregate hourly unit hydrograph to simulation time step - allocate(seg_uh(iSeg)%dat((iHrLast+nTsub-1)/nTsub),stat=ierr,errmsg=cmessage) + allocate(seg_uh(iSeg)%dat(ntdh),stat=ierr,errmsg=cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage)//': seg_uh%dat'; return; endif seg_uh(iSeg)%dat(:)=0._dp diff --git a/route/build/src/public_var.f90 b/route/build/src/public_var.f90 index 9faa9007..cc773107 100644 --- a/route/build/src/public_var.f90 +++ b/route/build/src/public_var.f90 @@ -94,7 +94,8 @@ MODULE public_var 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 character(len=strLen),public :: newFileFrequency = 'yearly' ! frequency for new output files (daily, monthly, yearly, single) - character(len=strLen),public :: timeAggregation = 'none' ! time aggregation for output: 'nh','nd','nm','ny' where n is number or 'none' no aggregation + character(len=strLen),public :: outputFrequency = '1' ! output frequency (integer for multiple of simulation time step or daily, monthly or yearly) + integer(i4b) ,public :: nOutFreq = integerMissing ! integer output frequency character(len=10) ,public :: routOpt = '0' ! routing scheme options 0: accum runoff, 1:IRF, 2:KWT, 3:KW, 4:MC, 5:DW integer(i4b) ,public :: doesBasinRoute = 1 ! basin routing options 0-> no, 1->IRF, otherwise error logical(lgt) ,public :: is_lake_sim = .false. ! logical if lakes are activated in simulation @@ -106,7 +107,7 @@ MODULE public_var logical(lgt) ,public :: is_vol_wm_jumpstart = .false. ! logical if true the volume is reset to target volume for the first time step of modeling logical(lgt) ,public :: suppress_runoff = .false. ! logical to suppress the read runoff to zero(0) logical(lgt) ,public :: suppress_P_Ep = .false. ! logical to suppress evaporation and precipitation to zero(0) - logical(lgt) ,public :: compWB = .true. ! logical if entire domain water balance is computed + logical(lgt) ,public :: compWB = .false. ! logical if entire domain water balance is computed real(dp) ,public :: dt = realMissing ! simulation time step (seconds) ! RIVER NETWORK TOPOLOGY character(len=strLen),public :: fname_ntopOld = '' ! old filename containing stream network topology information diff --git a/route/build/src/read_control.f90 b/route/build/src/read_control.f90 index b140d584..f0ea25af 100644 --- a/route/build/src/read_control.f90 +++ b/route/build/src/read_control.f90 @@ -127,7 +127,7 @@ SUBROUTINE read_control(ctl_fname, err, message) case(''); read(cData,*,iostat=io_error) is_vol_wm_jumpstart ! logical; jump to the first time step target volume is set to true case(''); read(cData,*,iostat=io_error) suppress_runoff ! logical; suppress the read runoff to zero (0) no host model case(''); read(cData,*,iostat=io_error) suppress_P_Ep ! logical; suppress the precipitation/evaporation to zero (0) no host model - case(''); read(cData,*,iostat=io_error) dt ! time interval of the simulation (To-do: change dt to dt_sim) + case(''); read(cData,*,iostat=io_error) dt ! time interval of the simulation [sec] (To-do: change dt to dt_sim) ! RIVER NETWORK TOPOLOGY case(''); fname_ntopOld = trim(cData) ! name of file containing stream network topology information case(''); read(cData,*,iostat=io_error) ntopAugmentMode ! option for river network augmentation mode. terminate the program after writing augmented ntopo. @@ -201,8 +201,8 @@ SUBROUTINE read_control(ctl_fname, err, message) case(''); read(cData,*,iostat=io_error) maxPfafLen ! maximum digit of pfafstetter code (default 32) case(''); pfafMissing = trim(cData) ! missing pfafcode (e.g., reach without any upstream area) ! OUTPUT OPTIONS - case(''); newFileFrequency = trim(cData) ! frequency for new output files (day, month, annual, single) - case(''); timeAggregation = trim(cData) ! time aggregation for output: 'nh','nd','nm','ny' where n is number + case(''); newFileFrequency = trim(cData) ! frequency for new history files (daily, monthly, yearly, single) + case(''); outputFrequency = trim(cData) ! output frequency (integer for multiple of simulation time step or daily, monthly or yearly) case(''); read(cData,*,iostat=io_error) meta_hflx(ixHFLX%basRunoff )%varFile ! default: true case(''); read(cData,*,iostat=io_error) meta_rflx(ixRFLX%instRunoff )%varFile ! default: false case(''); read(cData,*,iostat=io_error) meta_rflx(ixRFLX%dlayRunoff )%varFile ! default: false @@ -334,16 +334,17 @@ SUBROUTINE read_control(ctl_fname, err, message) end do ! looping through lines in the control file + ! ---------- Perform minor processing and checking control variables ---------------------------------------- + ! ---------- 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 - ! Case2- river network augment mode: Write full network variables over the entire network + ! option 1- river network subset mode (idSegOut>0): Write the network variables read from file over only upstream network specified idSegOut + ! option 2- river network augment mode: Write full network variables over the entire network ! River network subset mode turnes off augmentation mode. - ! Turned off ntopAugmentMode if (idSegOut>0) then ntopAugmentMode = .false. endif @@ -398,7 +399,61 @@ SUBROUTINE read_control(ctl_fname, err, message) err=81; return end select - ! ---------- output options -------------------------------------------------------------------------------------------- + ! ---------- simulation time step, output frequency, file frequency ------- + if (masterproc) then + write(iulog,'(2a)') new_line('a'), '---- output/simulation time steps --- ' + write(iulog,'(A,F10.1)') ' simulation frequency: ', dt + write(iulog,'(2A)') ' history file freqeuncy: ', trim(newFileFrequency) + write(iulog,'(2A)') ' history output freqeuncy: ', trim(outputFrequency) + end if + + ! 1. Process history output frequency + select case(trim(outputFrequency)) + case('daily', 'monthly', 'yearly') ! do nothing + case default + read(outputFrequency,'(5I)',iostat=err) nOutFreq + if (err/=0) then + message=trim(message)//' is invalid: must be daily, monthly, yearly or integer (number of time steps)'; return + end if + if (nOutFreq<0) then + message=trim(message)//' is invalid: must be positive integer AND outputFrequency x simulation step [sec] must be 86400 [sec] (one day)'; return + end if + end select + + ! 2. Check simulation time step + ! 2.1 must be less than one day + if (dt>86400._dp) then + write(message, '(2A)') trim(message), 'simulation frequency must be less than one-day (86400 sec)' + err=81; return + end if + ! 2.2. multiple of simulation time step must be one day + if (mod(86400._dp, dt)>0._dp) then + write(message, '(2A)') trim(message), 'multiple of simulation step [sec] must end up in one-day (86400 sec)' + err=81; return + end if + + ! 3. Check output frequency against simulation time step if outputFrequency is numeric + if (nOutFreq/=integerMissing) then + if (mod(86400._dp, real(nOutFreq,kind=dp)*dt)>0._dp) then + write(message, '(2A)') trim(message), 'outputFrequency x simulation step [sec] must be 86400 [sec] (one day)' + err=81; return + end if + end if + + ! 4. Check new history frequency against output frequency + if (trim(newFileFrequency)=='daily') then + if (trim(outputFrequency)=='monthly' .or. trim(outputFrequency)=='yearly') then + write(message, '(2A)') trim(message), 'you cannot output monthly or yearly output in daily file' + err=81; return + end if + else if (trim(newFileFrequency)=='monthly') then + if (trim(outputFrequency)=='yearly') then + write(message, '(2A)') trim(message), 'you cannot output yearly output in monthly file' + err=81; return + end if + end if + + ! ---- routing methods ! Assign index for each active routing method ! Make sure to turn off write option for routines not used if (trim(routOpt)=='0')then; write(iulog,'(a)') 'WARNING: routOpt=0 is accumRunoff option now. 12 is previous 0 now'; endif @@ -418,6 +473,7 @@ SUBROUTINE read_control(ctl_fname, err, message) end select end do + ! ---- history Output variables do iRoute = 0, nRouteMethods-1 select case(iRoute) case(accumRunoff) diff --git a/route/build/src/standalone/model_setup.f90 b/route/build/src/standalone/model_setup.f90 index 2844d5ff..18310f42 100644 --- a/route/build/src/standalone/model_setup.f90 +++ b/route/build/src/standalone/model_setup.f90 @@ -94,7 +94,7 @@ SUBROUTINE init_data(pid, & ! input: proc id call pass_global_data(comm, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - ! channel state initialization + ! restart initialization call init_state_data(pid, nNodes, comm, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif @@ -111,7 +111,7 @@ END SUBROUTINE init_data ! or abstraction or injection ! ********************************************************************* - SUBROUTINE init_inFile_pop (ierr, message) ! output + SUBROUTINE init_inFile_pop(ierr, message) USE public_var, ONLY: input_dir ! directory containing the text files of fname_qsim and fname_wm USE public_var, ONLY: fname_qsim ! simulated runoff txt file that includes the NetCDF file names @@ -570,8 +570,12 @@ SUBROUTINE init_time(ierr, message) endif - ! set initial model simulation time (beginning of simulation time step) + ! set initial model simulation time (beginning of simulation time step and next time step) simDatetime(1) = begDatetime + simDatetime(2) = simDatetime(1)%add_sec(dt, calendar, ierr, cmessage) + if (continue_run) then + simDatetime(0) = simDatetime(1)%add_sec(-dt, calendar, ierr, cmessage) + end if ! set simulation time step index (should be one to start) iTime = 1 @@ -590,10 +594,6 @@ SUBROUTINE init_time(ierr, message) ierr=20; message=trim(message)//'= '//trim(t_unit)//': must be seconds, minutes, hours or days.'; return end select - if (continue_run) then - simDatetime(0) = simDatetime(1)%add_sec(-dt, calendar, ierr, cmessage) - end if - ! Set restart calendar date/time and dropoff calendar date/time and ! -- For periodic restart options --------------------------------------------------------------------- ! Ensure that user-input restart month, day are valid. diff --git a/route/build/src/standalone/read_runoff.f90 b/route/build/src/standalone/read_runoff.f90 index af8e79bf..20729706 100644 --- a/route/build/src/standalone/read_runoff.f90 +++ b/route/build/src/standalone/read_runoff.f90 @@ -80,7 +80,7 @@ SUBROUTINE read_runoff_metadata(fname , & ! input: filename including else if (input_fillvalue==realMissing) then write(iulog, '(3A)') 'WARNING:', trim(var_name), '. No _FillValue exist in attributes nor provided by user in input_fillvalue' - write(iulog, '(A,x,F5.2)') ' _FillValue is set to default:', realMissing + write(iulog, '(A,x,G15.4)') ' _FillValue is set to default:', realMissing end if fillvalue = input_fillvalue end if diff --git a/route/build/src/var_lookup.f90 b/route/build/src/var_lookup.f90 index 624a9c26..6b413af3 100644 --- a/route/build/src/var_lookup.f90 +++ b/route/build/src/var_lookup.f90 @@ -34,15 +34,16 @@ MODULE var_lookup ! For routing state variables type, public :: iLook_stateDims integer(i4b) :: seg = integerMissing ! 1. stream segment vector - integer(i4b) :: time = integerMissing ! 2. time - integer(i4b) :: tbound = integerMissing ! 3. 2 elelment time bound vector - integer(i4b) :: ens = integerMissing ! 4. runoff ensemble - integer(i4b) :: wave = integerMissing ! 5. waves in a channel - integer(i4b) :: mol_kw = integerMissing ! 6. kw finite difference computational molecule - integer(i4b) :: mol_mc = integerMissing ! 7. mc finite difference computational molecule - integer(i4b) :: mol_dw = integerMissing ! 8. kw finite difference computational molecule - integer(i4b) :: tdh_irf = integerMissing ! 9. irf routed future channel flow in a segment - integer(i4b) :: tdh = integerMissing ! 10. uh routed future overland flow + integer(i4b) :: hru = integerMissing ! 2. catchment hru vector + integer(i4b) :: time = integerMissing ! 3. time + integer(i4b) :: tbound = integerMissing ! 4. 2 elelment time bound vector + integer(i4b) :: ens = integerMissing ! 5. runoff ensemble + integer(i4b) :: wave = integerMissing ! 6. waves in a channel + integer(i4b) :: mol_kw = integerMissing ! 7. kw finite difference computational molecule + integer(i4b) :: mol_mc = integerMissing ! 8. mc finite difference computational molecule + integer(i4b) :: mol_dw = integerMissing ! 9. kw finite difference computational molecule + integer(i4b) :: tdh_irf = integerMissing ! 10. irf routed future channel flow in a segment + integer(i4b) :: tdh = integerMissing ! 11. uh routed future overland flow endtype iLook_stateDims ! For river discharge variables type, public :: iLook_qDims @@ -257,7 +258,8 @@ MODULE var_lookup ! *********************************************************************************************************** type(iLook_struct) ,public,parameter :: ixStruct = iLook_struct ( 1, 2, 3, 4, 5) type(iLook_dims) ,public,parameter :: ixDims = iLook_dims ( 1, 2, 3, 4, 5, 6, 7) - type(iLook_stateDims),public,parameter :: ixStateDims = iLook_stateDims( 1, 2, 3, 4, 5, 6, 7, 8, 9,10) + type(iLook_stateDims),public,parameter :: ixStateDims = iLook_stateDims( 1, 2, 3, 4, 5, 6, 7, 8, 9,10, & + 11) 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) diff --git a/route/build/src/write_restart_pio.f90 b/route/build/src/write_restart_pio.f90 index 17b7f3cb..b91b0f7d 100644 --- a/route/build/src/write_restart_pio.f90 +++ b/route/build/src/write_restart_pio.f90 @@ -3,6 +3,8 @@ MODULE write_restart_pio ! Moudle wide shared data/external routines USE nrtype +USE var_lookup, ONLY: ixRFLX, nVarsRFLX +USE var_lookup, ONLY: ixHFLX, nVarsHFLX USE var_lookup, ONLY: ixStateDims, nStateDims USE var_lookup, ONLY: ixIRFbas, nVarsIRFbas USE var_lookup, ONLY: ixIRF, nVarsIRF @@ -12,19 +14,18 @@ MODULE write_restart_pio USE var_lookup, ONLY: ixDW, nVarsDW USE var_lookup, ONLY: ixIRFbas, nVarsIRFbas USE var_lookup, ONLY: ixBasinQ, nVarsBasinQ - +! data daype USE dataTypes, ONLY: STRFLX ! fluxes in each reach USE dataTypes, ONLY: STRSTA ! state in each reach USE dataTypes, ONLY: RCHTOPO ! Network topology -USE dataTypes, ONLY: states USE datetime_data, ONLY: datetime - +! public variables USE public_var, ONLY: iulog ! i/o logical unit number USE public_var, ONLY: integerMissing USE public_var, ONLY: realMissing -USE public_var, ONLY: verySmall - +! meta data USE globalData, ONLY: meta_stateDims ! states dimension meta +USE globalData, ONLY: meta_qDims USE globalData, ONLY: meta_irf_bas USE globalData, ONLY: meta_basinQ USE globalData, ONLY: meta_irf @@ -32,7 +33,9 @@ MODULE write_restart_pio USE globalData, ONLY: meta_kw USE globalData, ONLY: meta_mc USE globalData, ONLY: meta_dw - +USE globalData, ONLY: meta_rflx +USE globalData, ONLY: meta_hflx +! pio stuff USE globalData, ONLY: pid, nNodes USE globalData, ONLY: masterproc USE globalData, ONLY: mpicom_route @@ -43,9 +46,12 @@ MODULE write_restart_pio USE globalData, ONLY: pio_root USE globalData, ONLY: pio_stride USE globalData, ONLY: pioSystem + USE globalData, ONLY: runMode USE globalData, ONLY: rfileout - +USE globalData, ONLY: idxSUM,idxIRF,idxKWT, & + idxKW,idxMC,idxDW +! external procedures USE nr_utils, ONLY: arth USE pio_utils @@ -53,8 +59,10 @@ MODULE write_restart_pio ! The following variables used only in this module type(file_desc_t) :: pioFileDescState ! contains data identifying the file -type(io_desc_t) :: iodesc_state_int -type(io_desc_t) :: iodesc_state_double +type(io_desc_t) :: iodesc_rch_int +type(io_desc_t) :: iodesc_rch_double +type(io_desc_t) :: ioDesc_rch_float +type(io_desc_t) :: ioDesc_hru_float type(io_desc_t) :: iodesc_wave_int type(io_desc_t) :: iodesc_wave_double type(io_desc_t) :: iodesc_mesh_kw_double @@ -198,9 +206,7 @@ SUBROUTINE restart_fname(fname, timeStamp, ierr, message) USE public_var, ONLY: restart_dir USE public_var, ONLY: case_name ! simulation name ==> output filename head - USE public_var, ONLY: calendar USE public_var, ONLY: secprday - USE public_var, ONLY: dt USE globalData, ONLY: simDatetime ! current model datetime implicit none @@ -213,7 +219,6 @@ SUBROUTINE restart_fname(fname, timeStamp, ierr, message) character(*), intent(out) :: message ! error message ! local variables type(datetime) :: restartTimeStamp ! datetime corresponding to file name time stamp - character(len=strLen) :: cmessage ! error message of downwind routine integer(i4b) :: sec_in_day ! second within day 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=50),parameter :: fmtYMDS='(a,I0.4,a,I0.2,a,I0.2,a,I0.5,a)' @@ -222,7 +227,7 @@ SUBROUTINE restart_fname(fname, timeStamp, ierr, message) select case(timeStamp) case(currTimeStep); restartTimeStamp = simDatetime(1) - case(nextTimeStep); restartTimeStamp = simDatetime(1)%add_sec(dt, calendar, ierr, cmessage) + case(nextTimeStep); restartTimeStamp = simDatetime(2) 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 @@ -262,8 +267,7 @@ SUBROUTINE define_state_nc(fname, & ! input: filename USE public_var, ONLY: muskingumCunge USE public_var, ONLY: diffusiveWave USE globalData, ONLY: onRoute ! logical to indicate which routing method(s) is on - USE globalData, ONLY: rch_per_proc ! number of reaches assigned to each proc (size = num of procs+1) - USE globalData, ONLY: nRch ! number of ensembles and river reaches + USE write_simoutput_pio, ONLY:get_compdof_all_network implicit none ! argument variables @@ -271,10 +275,10 @@ SUBROUTINE define_state_nc(fname, & ! input: filename integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! local variables - integer(i4b) :: ix1, ix2 ! frst and last indices of global array for local array chunk - integer(i4b) :: ixRch(nRch) ! global reach index used for output + integer(i4b), allocatable :: compdof_rch(:) ! + integer(i4b), allocatable :: compdof_hru(:) ! integer(i4b) :: jDim ! loop index for dimension - integer(i4b) :: ixDim_common(3) ! custom dimension ID array + integer(i4b) :: ixDim_common(4) ! custom dimension ID array character(len=strLen) :: cmessage ! error message of downwind routine ! Initialize error control @@ -302,7 +306,7 @@ SUBROUTINE define_state_nc(fname, & ! input: filename if(ierr/=0)then; message=trim(cmessage)//'cannot create state netCDF'; return; endif ! For common dimension/variables - seg id, time, time-bound ----------- - ixDim_common = [ixStateDims%seg, ixStateDims%ens, ixStateDims%tbound] + ixDim_common = [ixStateDims%seg, ixStateDims%hru, ixStateDims%ens, ixStateDims%tbound] ! ---------------------------------- ! Define dimensions @@ -323,6 +327,9 @@ SUBROUTINE define_state_nc(fname, & ! input: filename call def_var(pioFileDescState, 'nNodes', ncd_int, ierr, cmessage, vdesc='Number of MPI tasks', vunit='-' ) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + call def_var(pioFileDescState, 'nt', ncd_int, ierr, cmessage, vdesc='Number of current acculated time steps in history variable', vunit='-' ) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + call def_var(pioFileDescState, 'reachID', ncd_int, ierr, cmessage, pioDimId=[dim_seg], vdesc='reach ID', vunit='-' ) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif @@ -370,10 +377,15 @@ SUBROUTINE define_state_nc(fname, & ! input: filename if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif end if + ! accumulated history variables + call define_history_state(ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + ! ---------------------------------- ! pio initialization of decomposition ! ---------------------------------- associate(nSeg => meta_stateDims(ixStateDims%seg)%dimLength, & + nHru => meta_stateDims(ixStateDims%hru)%dimLength, & nEns => meta_stateDims(ixStateDims%ens)%dimLength, & ntdh => meta_stateDims(ixStateDims%tdh)%dimLength, & ! maximum future q time steps among basins ntdh_irf => meta_stateDims(ixStateDims%tdh_irf)%dimLength, & ! maximum future q time steps among reaches @@ -383,34 +395,44 @@ SUBROUTINE define_state_nc(fname, & ! input: filename nMesh_dw => meta_stateDims(ixStateDims%mol_dw)%dimLength, & ! dw_finite difference mesh points nWave => meta_stateDims(ixStateDims%wave)%dimLength) ! maximum waves allowed in a reach - if (masterproc) then - ix1 = 1 - else - ix1 = sum(rch_per_proc(-1:pid-1))+1 - endif - ix2 = sum(rch_per_proc(-1:pid)) - ixRch = arth(1,1,nSeg) + call get_compdof_all_network(compdof_rch, compdof_hru) ! type: float dim: [dim_seg, dim_ens] -- channel runoff coming from hru call pio_decomp(pioSystem, & ! input: pio system descriptor ncd_double, & ! input: data type (pio_int, pio_real, pio_double, pio_char) [nSeg,nEns], & ! input: dimension length == global array size - ixRch(ix1:ix2), & ! input: - iodesc_state_double) + compdof_rch, & ! input: + iodesc_rch_double) ! type: int dim: [dim_seg, dim_ens] -- number of wave or uh future time steps call pio_decomp(pioSystem, & ! input: pio system descriptor ncd_int, & ! input: data type (pio_int, pio_real, pio_double, pio_char) [nSeg,nEns], & ! input: dimension length == global array size - ixRch(ix1:ix2), & ! input: - iodesc_state_int) + compdof_rch, & ! input: + iodesc_rch_int) + + ! type: single precision float dim: [dim_seg, dim_ens] -- + call pio_decomp(pioSystem, & ! input: pio system descriptor + ncd_float, & ! input: data type (pio_int, pio_real, pio_double, pio_char) + [nSeg], & ! input: dimension length == global array size + compdof_rch, & ! input: + ioDesc_rch_float) + + if (meta_hflx(ixHFLX%basRunoff)%varFile) then + ! type: single precision float dim: [dim_hru, dim_ens] -- + call pio_decomp(pioSystem, & ! input: pio system descriptor + ncd_float, & ! input: data type (pio_int, pio_real, pio_double, pio_char) + [nHru], & ! input: dimension length == global array size + compdof_hru, & ! input: + ioDesc_hru_float) + end if if (doesBasinRoute==1) then ! type: float dim: [dim_seg, dim_tdh_irf, dim_ens] call pio_decomp(pioSystem, & ! input: pio system descriptor ncd_double, & ! input: data type (pio_int, pio_real, pio_double, pio_char) [nSeg,ntdh,nEns], & ! input: dimension length == global array size - ixRch(ix1:ix2), & ! input: + compdof_rch, & ! input: iodesc_irf_bas_double) end if @@ -419,14 +441,14 @@ SUBROUTINE define_state_nc(fname, & ! input: filename call pio_decomp(pioSystem, & ! input: pio system descriptor ncd_double, & ! input: data type (pio_int, pio_real, pio_double, pio_char) [nSeg,ntdh_irf,nEns], & ! input: dimension length == global array size - ixRch(ix1:ix2), & ! input: + compdof_rch, & ! input: iodesc_irf_double) ! type: float dim: [dim_seg, dim_ens] call pio_decomp(pioSystem, & ! input: pio system descriptor ncd_double, & ! input: data type (pio_int, pio_real, pio_double, pio_char) [nSeg,nTbound,nEns], & ! input: dimension length == global array size - ixRch(ix1:ix2), & ! input: + compdof_rch, & ! input: iodesc_vol_double) end if @@ -435,14 +457,14 @@ SUBROUTINE define_state_nc(fname, & ! input: filename call pio_decomp(pioSystem, & ! input: pio system descriptor ncd_int, & ! input: data type (pio_int, pio_real, pio_double, pio_char) [nSeg,nWave,nEns], & ! input: dimension length == global array size - ixRch(ix1:ix2), & ! input: + compdof_rch, & ! input: iodesc_wave_int) ! type: float, dim: [dim_seg, dim_wave, dim_ens] call pio_decomp(pioSystem, & ! input: pio system descriptor ncd_double, & ! input: data type (pio_int, pio_real, pio_double, pio_char) [nSeg,nWave,nEns], & ! input: dimension length == global array size - ixRch(ix1:ix2), & ! input: + compdof_rch, & ! input: iodesc_wave_double) end if @@ -451,7 +473,7 @@ SUBROUTINE define_state_nc(fname, & ! input: filename call pio_decomp(pioSystem, & ! input: pio system descriptor ncd_double, & ! input: data type (pio_int, pio_real, pio_double, pio_char) [nSeg,nMesh_kw,nEns], & ! input: dimension length == global array size - ixRch(ix1:ix2), & ! input: + compdof_rch, & ! input: iodesc_mesh_kw_double) end if @@ -460,7 +482,7 @@ SUBROUTINE define_state_nc(fname, & ! input: filename call pio_decomp(pioSystem, & ! input: pio system descriptor ncd_double, & ! input: data type (pio_int, pio_real, pio_double, pio_char) [nSeg,nMesh_mc,nEns], & ! input: dimension length == global array size - ixRch(ix1:ix2), & ! input: + compdof_rch, & ! input: iodesc_mesh_mc_double) end if @@ -469,7 +491,7 @@ SUBROUTINE define_state_nc(fname, & ! input: filename call pio_decomp(pioSystem, & ! input: pio system descriptor ncd_double, & ! input: data type (pio_int, pio_real, pio_double, pio_char) [nSeg,nMesh_dw,nEns], & ! input: dimension length == global array size - ixRch(ix1:ix2), & ! input: + compdof_rch, & ! input: iodesc_mesh_dw_double) end if @@ -485,8 +507,9 @@ SUBROUTINE set_dim_len(ixDim, ierr, message1) ! populate state netCDF dimension size USE public_var, ONLY: MAXQPAR USE globalData, ONLY: nMolecule - USE globalData, ONLY: FRAC_FUTURE ! To get size of q future for basin IRF - USE globalData, ONLY: nEns, nRch ! number of ensembles and river reaches + USE globalData, ONLY: maxtdh ! maximum unit-hydrogrph future time + USE globalData, ONLY: FRAC_FUTURE ! To get size of q future for basin IRF + USE globalData, ONLY: nEns, nRch, nHru ! number of ensembles and river reaches implicit none ! input integer(i4b), intent(in) :: ixDim ! ixDim @@ -498,10 +521,11 @@ SUBROUTINE set_dim_len(ixDim, ierr, message1) select case(ixDim) case(ixStateDims%seg); meta_stateDims(ixStateDims%seg)%dimLength = nRch + case(ixStateDims%hru); meta_stateDims(ixStateDims%hru)%dimLength = nHru case(ixStateDims%ens); meta_stateDims(ixStateDims%ens)%dimLength = nEns 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 = maxtdh case(ixStateDims%mol_kw); meta_stateDims(ixStateDims%mol_kw)%dimLength = nMolecule%KW_ROUTE case(ixStateDims%mol_mc); meta_stateDims(ixStateDims%mol_mc)%dimLength = nMolecule%MC_ROUTE case(ixStateDims%mol_dw); meta_stateDims(ixStateDims%mol_dw)%dimLength = nMolecule%DW_ROUTE @@ -786,6 +810,35 @@ SUBROUTINE define_DW_state(ierr, message1) END SUBROUTINE define_DW_state + SUBROUTINE define_history_state(ierr, message1) + implicit none + ! output + integer(i4b), intent(out) :: ierr ! error code + character(*), intent(out) :: message1 ! error message + ! local + integer(i4b) :: iVar ! index loop for variables + integer(i4b) :: dim_set(1) ! dimension Id array + + ierr=0; message1='define_history_state/' + + do iVar=1,nVarsRFLX + if (.not. meta_rflx(iVar)%varFile) cycle + dim_set(1) = meta_stateDims(ixStateDims%seg)%dimId ! this should be seg dimension + call def_var(pioFileDescState, meta_rflx(iVar)%varName, meta_rflx(iVar)%varType, ierr, cmessage, & + pioDimId=dim_set, vdesc=meta_rflx(iVar)%varDesc, vunit=meta_rflx(iVar)%varUnit) + if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif + end do + + do iVar=1,nVarsHFLX + if (.not. meta_hflx(iVar)%varFile) cycle + dim_set(1) = meta_stateDims(ixStateDims%hru)%dimId ! this should be hru dimension + call def_var(pioFileDescState, meta_hflx(iVar)%varName, meta_hflx(iVar)%varType, ierr, cmessage, & + pioDimId=dim_set, vdesc=meta_hflx(iVar)%varDesc, vunit=meta_hflx(iVar)%varUnit) + if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif + end do + + END SUBROUTINE define_history_state + END SUBROUTINE define_state_nc @@ -816,6 +869,7 @@ SUBROUTINE write_state_nc(fname, & ! Input: state netcdf name USE globalData, ONLY: nRch ! number of reaches in network USE globalData, ONLY: TSEC ! beginning/ending of simulation time step [sec] USE globalData, ONLY: timeVar ! time variables (unit given by runoff data) + USE write_simoutput_pio, ONLY:hVars ! current history variable data implicit none @@ -827,7 +881,6 @@ SUBROUTINE write_state_nc(fname, & ! Input: state netcdf name ! local variables real(dp) :: secPerTime ! number of sec per time-unit. time-unit is from t_unit real(dp) :: restartTimeVar ! restart timeVar [time_units] - real(dp) :: tb_array(2) ! restart timeVar [time_units] integer(i4b) :: iens ! temporal integer(i4b) :: nRch_local ! number of reach in each processors integer(i4b) :: nRch_root ! number of reaches in root processors (including halo reaches) @@ -881,11 +934,10 @@ SUBROUTINE write_state_nc(fname, & ! Input: state netcdf name restartTimeVar = timeVar + dt/secPerTime - ! -- Write out to netCDF - call openFile(pioSystem, pioFileDescState, trim(fname),pio_typename, ncd_write, restartOpen, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + ! -- Write out to netCDF call write_scalar_netcdf(pioFileDescState, 'nNodes', nNodes, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif @@ -895,8 +947,7 @@ SUBROUTINE write_state_nc(fname, & ! Input: state netcdf name call write_scalar_netcdf(pioFileDescState, 'restart_time', restartTimeVar, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - tb_array = [TSEC(0),TSEC(1)] - call write_netcdf(pioFileDescState, 'time_bound', tb_array, [1], [2], ierr, cmessage) + call write_netcdf(pioFileDescState, 'time_bound', TSEC, [1], [2], ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif call write_basinQ_state(ierr, cmessage) @@ -932,9 +983,16 @@ SUBROUTINE write_state_nc(fname, & ! Input: state netcdf name if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif end if + call write_history_state(ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + ! clean decomposition data - call freeDecomp(pioFileDescState, iodesc_state_double) - call freeDecomp(pioFileDescState, iodesc_state_int) + call freeDecomp(pioFileDescState, iodesc_rch_double) + call freeDecomp(pioFileDescState, iodesc_rch_int) + call freeDecomp(pioFileDescState, iodesc_rch_float) + if (meta_hflx(ixHFLX%basRunoff)%varFile) then + call freeDecomp(pioFileDescState, iodesc_hru_float) + end if if (doesBasinRoute==1) then call freeDecomp(pioFileDescState, iodesc_irf_bas_double) end if @@ -956,6 +1014,7 @@ SUBROUTINE write_state_nc(fname, & ! Input: state netcdf name call freeDecomp(pioFileDescState, iodesc_mesh_dw_double) end if + ! close netCDF call closeFile(pioFileDescState, restartOpen) CONTAINS @@ -966,7 +1025,7 @@ SUBROUTINE write_basinQ_state(ierr, message1) integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message1 ! error message ! local variables - type(states) :: state ! temporal state data structures -currently 2 river routing scheme + basin IRF routing + real(dp), allocatable :: array_2d_dp(:,:) integer(i4b) :: iVar,iens,iSeg ! index loops for variables, ensembles and segments respectively ! initialize error control @@ -975,35 +1034,20 @@ SUBROUTINE write_basinQ_state(ierr, message1) associate(nSeg => size(RCHFLX_local), & nens => meta_stateDims(ixStateDims%ens)%dimLength) - allocate(state%var(nVarsBasinQ), stat=ierr, errmsg=cmessage) - if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif - do iVar=1,nVarsBasinQ select case(iVar) - case(ixBasinQ%q); allocate(state%var(iVar)%array_2d_dp(nSeg, nEns), stat=ierr) - case default; ierr=20; message1=trim(message1)//'unable to identify basin routing variable index'; return + case(ixBasinQ%q) + allocate(array_2d_dp(nSeg, nEns), stat=ierr, errmsg=cmessage) + if(ierr/=0)then; message1=trim(message1)//trim(cmessage)//':basin runoff:'//trim(meta_mc(iVar)%varName); return; endif + do iens=1,nens + do iSeg=1,nSeg + array_2d_dp(iSeg,iens) = RCHFLX_local(iSeg)%BASIN_QR(1) + enddo + enddo + call write_pnetcdf(pioFileDescState, meta_basinQ(iVar)%varName, array_2d_dp, iodesc_rch_double, ierr, cmessage) + deallocate(array_2d_dp) + case default; ierr=20; message1=trim(message1)//'unable to identify basin runoff variable index'; return end select - if(ierr/=0)then; message1=trim(message1)//'problem allocating space for basin IRF routing state '//trim(meta_basinQ(iVar)%varName); return; endif - end do - - ! --Convert data structures to arrays - do iens=1,nens - do iSeg=1,nSeg - do iVar=1,nVarsBasinQ - select case(iVar) - case(ixBasinQ%q); state%var(iVar)%array_2d_dp(iSeg,iens) = RCHFLX_local(iSeg)%BASIN_QR(1) - case default; ierr=20; message1=trim(message1)//'unable to identify basin IRF state variable index'; return - end select - enddo - enddo - enddo - - do iVar=1,nVarsBasinQ - select case(iVar) - case(ixBasinQ%q); call write_pnetcdf(pioFileDescState, meta_basinQ(iVar)%varName, state%var(iVar)%array_2d_dp, iodesc_state_double, ierr, cmessage) - case default; ierr=20; message1=trim(message1)//'unable to identify reach inflow variable index for nc writing'; return - end select - if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif end do end associate @@ -1016,7 +1060,7 @@ SUBROUTINE write_IRFbas_state(ierr, message1) integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message1 ! error message ! local variables - type(states) :: state ! temporal state data structures -currently 2 river routing scheme + basin IRF routing + real(dp), allocatable :: array_3d_dp(:,:,:) integer(i4b) :: iVar,iens,iSeg ! index loops for variables, ensembles and segments respectively ierr=0; message1='write_IRFbas_state/' @@ -1025,38 +1069,21 @@ SUBROUTINE write_IRFbas_state(ierr, message1) nEns => meta_stateDims(ixStateDims%ens)%dimLength, & ntdh => meta_stateDims(ixStateDims%tdh)%dimLength) ! maximum future q time steps among basins - allocate(state%var(nVarsIRFbas), stat=ierr, errmsg=cmessage) - if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif - - do iVar=1,nVarsIRFbas - - select case(iVar) - case(ixIRFbas%qfuture); allocate(state%var(iVar)%array_3d_dp(nSeg, ntdh, nEns), stat=ierr) - case default; ierr=20; message1=trim(message1)//'unable to identify basin routing variable index'; return - end select - if(ierr/=0)then; message1=trim(message1)//'problem allocating space for basin IRF routing state '//trim(meta_irf_bas(iVar)%varName); return; endif - - end do - - ! --Convert data structures to arrays - do iens=1,nEns - do iSeg=1,nSeg - do iVar=1,nVarsIRFbas - select case(iVar) - case(ixIRFbas%qfuture); state%var(iVar)%array_3d_dp(iSeg,:,iens) = RCHFLX_local(iSeg)%QFUTURE - case default; ierr=20; message1=trim(message1)//'unable to identify basin IRF state variable index'; return - end select - enddo - enddo - enddo - do iVar=1,nVarsIRFbas select case(iVar) - case(ixIRFbas%qfuture); call write_pnetcdf(pioFileDescState, meta_irf_bas(iVar)%varName, state%var(iVar)%array_3d_dp, iodesc_irf_bas_double, ierr, cmessage) - case default; ierr=20; message1=trim(message1)//'unable to identify basin IRF variable index for nc writing'; return + case(ixIRFbas%qfuture) + allocate(array_3d_dp(nSeg, ntdh, nEns), stat=ierr, errmsg=cmessage) + if(ierr/=0)then; message1=trim(message1)//trim(cmessage)//':basin IRF routing state:'//trim(meta_mc(iVar)%varName); return; endif + do iens=1,nEns + do iSeg=1,nSeg + array_3d_dp(iSeg,:,iens) = RCHFLX_local(iSeg)%QFUTURE + end do + end do + call write_pnetcdf(pioFileDescState, meta_irf_bas(iVar)%varName, array_3d_dp, iodesc_irf_bas_double, ierr, cmessage) + deallocate(array_3d_dp) + case default; ierr=20; message1=trim(message1)//'unable to identify basin IRF state variable index'; return end select - if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif - enddo + end do end associate @@ -1069,9 +1096,9 @@ SUBROUTINE write_IRF_state(ierr, message1) integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message1 ! error message ! local variables - type(states) :: state ! temporal state data structures -currently 2 river routing scheme + basin IRF routing - integer(i4b) :: iVar,iens,iSeg ! index loops for variables, ensembles and segments respectively - integer(i4b), allocatable :: numQF(:,:) ! number of future Q time steps for each ensemble and segment + real(dp), allocatable :: array_3d_dp(:,:,:) + integer(i4b) :: iVar,iens,iSeg ! index loops for variables, ensembles and segments respectively + integer(i4b), allocatable :: numQF(:,:) ! number of future Q time steps for each ensemble and segment ierr=0; message1='write_IRF_state/' @@ -1080,51 +1107,43 @@ SUBROUTINE write_IRF_state(ierr, message1) nTbound => meta_stateDims(ixStateDims%tbound)%dimLength, & ntdh_irf => meta_stateDims(ixStateDims%tdh_irf)%dimLength) ! maximum future q time steps among reaches - allocate(state%var(nVarsIRF), stat=ierr, errmsg=cmessage) - if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif - ! array to store number of wave per segment and ensemble - allocate(numQF(nEns,nSeg), stat=ierr) - if(ierr/=0)then; message1=trim(message1)//'problem allocating space for numQF'; return; endif - - do iVar=1,nVarsIRF - select case(iVar) - case(ixIRF%qfuture); allocate(state%var(iVar)%array_3d_dp(nSeg, ntdh_irf, nEns), stat=ierr) - case(ixIRF%vol); allocate(state%var(iVar)%array_3d_dp(nSeg, nTbound, 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 + allocate(numQF(nEns,nSeg), stat=ierr, errmsg=cmessage) + if(ierr/=0)then; message1=trim(message1)//trim(cmessage)//':numQF'; return; endif - ! --Convert data structures to arrays do iens=1,nEns do iSeg=1,nSeg numQF(iens,iseg) = size(NETOPO_local(iSeg)%UH) - do iVar=1,nVarsIRF - select case(iVar) - case(ixIRF%qfuture) - state%var(iVar)%array_3d_dp(iSeg,1:numQF(iens,iSeg),iens) = RCHFLX_local(iSeg)%QFUTURE_IRF - state%var(iVar)%array_3d_dp(iSeg,numQF(iens,iSeg)+1:ntdh_irf,iens) = realMissing - case(ixIRF%vol) - state%var(iVar)%array_3d_dp(iSeg,1:nTbound,iens) = RCHFLX_local(iSeg)%ROUTE(idxIRF)%REACH_VOL(0:1) - case default; ierr=20; message1=trim(message1)//'unable to identify variable index'; return - end select - enddo ! variable loop - enddo ! seg loop - enddo ! ensemble loop - - ! writing netcdf - call write_pnetcdf(pioFileDescState, 'numQF', numQF, iodesc_state_int, ierr, cmessage) + end do + end do + + call write_pnetcdf(pioFileDescState, 'numQF', numQF, iodesc_rch_int, ierr, cmessage) if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif do iVar=1,nVarsIRF select case(iVar) case(ixIRF%qfuture) - call write_pnetcdf(pioFileDescState, meta_irf(iVar)%varName, state%var(iVar)%array_3d_dp, iodesc_irf_double, ierr, cmessage) + allocate(array_3d_dp(nSeg,ntdh_irf,nEns), stat=ierr, errmsg=cmessage) + if(ierr/=0)then; message1=trim(message1)//trim(cmessage)//':IRF routing state:'//trim(meta_mc(iVar)%varName); return; endif + do iens=1,nEns + do iSeg=1,nSeg + array_3d_dp(iSeg,1:numQF(iens,iSeg),iens) = RCHFLX_local(iSeg)%QFUTURE_IRF + array_3d_dp(iSeg,numQF(iens,iSeg)+1:ntdh_irf,iens) = realMissing + end do + end do + call write_pnetcdf(pioFileDescState, meta_irf(iVar)%varName, array_3d_dp, iodesc_irf_double, ierr, cmessage) + deallocate(array_3d_dp) case(ixIRF%vol) - call write_pnetcdf(pioFileDescState, meta_irf(iVar)%varName, state%var(iVar)%array_3d_dp, iodesc_vol_double, ierr, cmessage) + allocate(array_3d_dp(nSeg,nTbound,nEns), stat=ierr, errmsg=cmessage) + if(ierr/=0)then; message1=trim(message1)//trim(cmessage)//':IRF routing state:'//trim(meta_mc(iVar)%varName); return; endif + do iens=1,nEns + do iSeg=1,nSeg + array_3d_dp(iSeg,1:nTbound,iens) = RCHFLX_local(iSeg)%ROUTE(idxIRF)%REACH_VOL(0:1) + end do + end do + call write_pnetcdf(pioFileDescState, meta_irf(iVar)%varName, array_3d_dp, iodesc_vol_double, ierr, cmessage) + deallocate(array_3d_dp) 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 end do @@ -1133,90 +1152,114 @@ SUBROUTINE write_IRF_state(ierr, message1) END SUBROUTINE write_IRF_state SUBROUTINE write_KWT_state(ierr, message1) - USE globalData, ONLY: idxKWT - implicit none - ! output - integer(i4b), intent(out) :: ierr ! error code - character(*), intent(out) :: message1 ! error message - ! local variables - type(states) :: state ! temporal state data structures -currently 2 river routing scheme + basin IRF routing - integer(i4b) :: iVar,iens,iSeg ! index loops for variables, ensembles and segments respectively - integer(i4b), allocatable :: RFvec(:) ! temporal vector - integer(i4b), allocatable :: numWaves(:,:) ! number of waves for each ensemble and segment - - ierr=0; message1='write_KWT_state/' - - associate(nSeg => size(RCHFLX_local), & - nEns => meta_stateDims(ixStateDims%ens)%dimLength, & - nWave => meta_stateDims(ixStateDims%wave)%dimLength) ! maximum waves allowed in a reach - - allocate(state%var(nVarsKWT), stat=ierr, errmsg=cmessage) - if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif - - ! array to store number of wave per segment and ensemble - allocate(numWaves(nEns,nSeg), stat=ierr, errmsg=cmessage) - if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif - - do iVar=1,nVarsKWT - select case(iVar) - case(ixKWT%routed); allocate(state%var(iVar)%array_3d_int(nSeg, nWave, nEns), stat=ierr) - case(ixKWT%tentry, ixKWT%texit, ixKWT%qwave, ixKWT%qwave_mod) - allocate(state%var(iVar)%array_3d_dp(nSeg, nWave, nEns), stat=ierr) - case(ixKWT%vol); allocate(state%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 KWT routing state '//trim(meta_kwt(iVar)%varName); return; endif - end do - - ! --Convert data structures to arrays - do iens=1,nEns - do iSeg=1,nSeg - numWaves(iens,iseg) = size(RCHSTA_local(iseg)%LKW_ROUTE%KWAVE) - do iVar=1,nVarsKWT - select case(iVar) - case(ixKWT%tentry) - state%var(iVar)%array_3d_dp(iSeg,1:numWaves(iens,iSeg),iens) = RCHSTA_local(iSeg)%LKW_ROUTE%KWAVE(:)%TI - state%var(iVar)%array_3d_dp(iSeg,numWaves(iens,iSeg)+1:,iens) = realMissing - case(ixKWT%texit) - state%var(iVar)%array_3d_dp(iSeg,1:numWaves(iens,iSeg),iens) = RCHSTA_local(iSeg)%LKW_ROUTE%KWAVE(:)%TR - state%var(iVar)%array_3d_dp(iSeg,numWaves(iens,iSeg)+1:,iens) = realMissing - case(ixKWT%qwave) - state%var(iVar)%array_3d_dp(iSeg,1:numWaves(iens,iSeg),iens) = RCHSTA_local(iSeg)%LKW_ROUTE%KWAVE(:)%QF - state%var(iVar)%array_3d_dp(iSeg,numWaves(iens,iSeg)+1:,iens) = realMissing - case(ixKWT%qwave_mod) - state%var(iVar)%array_3d_dp(iSeg,1:numWaves(iens,iSeg),iens) = RCHSTA_local(iSeg)%LKW_ROUTE%KWAVE(:)%QM - state%var(iVar)%array_3d_dp(iSeg,numWaves(iens,iSeg)+1:,iens) = realMissing - case(ixKWT%routed) ! this is suppposed to be logical variable, but put it as 0 or 1 in double now - if (allocated(RFvec)) deallocate(RFvec, stat=ierr) - allocate(RFvec(numWaves(iens,iSeg)),stat=ierr); RFvec=0_i4b - where (RCHSTA_local(iSeg)%LKW_ROUTE%KWAVE(:)%RF) RFvec=1_i4b - state%var(iVar)%array_3d_int(iSeg,1:numWaves(iens,iSeg),iens) = RFvec - state%var(iVar)%array_3d_int(iSeg,numWaves(iens,iSeg)+1:,iens) = integerMissing - case(ixKWT%vol); state%var(iVar)%array_2d_dp(iSeg,iens) = RCHFLX_local(iSeg)%ROUTE(idxKWT)%REACH_VOL(1) - case default; ierr=20; message1=trim(message1)//'unable to identify KWT routing state variable index'; return - end select - enddo ! variable loop - enddo ! seg loop - enddo ! ensemble loop - - ! Writing netCDF - call write_pnetcdf(pioFileDescState,'numWaves', numWaves, iodesc_state_int, ierr, cmessage) - if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif - - do iVar=1,nVarsKWT - select case(iVar) - case(ixKWT%routed) - call write_pnetcdf(pioFileDescState, meta_kwt(iVar)%varName, state%var(iVar)%array_3d_int, iodesc_wave_int, ierr, cmessage) - case(ixKWT%tentry, ixKWT%texit, ixKWT%qwave, ixKWT%qwave_mod) - call write_pnetcdf(pioFileDescState, meta_kwt(iVar)%varName, state%var(iVar)%array_3d_dp, iodesc_wave_double, ierr, cmessage) - case(ixKWT%vol) - call write_pnetcdf(pioFileDescState, meta_kwt(iVar)%varName, state%var(iVar)%array_2d_dp, iodesc_state_double, ierr, cmessage) - case default; ierr=20; message1=trim(message1)//'unable to identify KWT variable index for nc writing'; return - end select + USE globalData, ONLY: idxKWT + implicit none + ! output + integer(i4b), intent(out) :: ierr ! error code + character(*), intent(out) :: message1 ! error message + ! local variables + real(dp), allocatable :: array_2d_dp(:,:) + real(dp), allocatable :: array_3d_dp(:,:,:) + integer(i4b), allocatable :: array_3d_int(:,:,:) + integer(i4b) :: iVar,iens,iSeg ! index loops for variables, ensembles and segments respectively + integer(i4b), allocatable :: RFvec(:) ! temporal vector + integer(i4b), allocatable :: numWaves(:,:) ! number of waves for each ensemble and segment + + ierr=0; message1='write_KWT_state/' + + associate(nSeg => size(RCHFLX_local), & + nEns => meta_stateDims(ixStateDims%ens)%dimLength, & + nWave => meta_stateDims(ixStateDims%wave)%dimLength) ! maximum waves allowed in a reach + + ! array to store number of wave per segment and ensemble + allocate(numWaves(nEns,nSeg), stat=ierr, errmsg=cmessage) if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif - end do - end associate + do iens=1,nEns + do iSeg=1,nSeg + numWaves(iens,iseg) = size(RCHSTA_local(iseg)%LKW_ROUTE%KWAVE) + end do + end do + + call write_pnetcdf(pioFileDescState,'numWaves', numWaves, iodesc_rch_int, ierr, cmessage) + if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif + + do iVar=1,nVarsKWT + select case(iVar) + case(ixKWT%routed) + allocate(array_3d_int(nSeg, nWave, nEns), stat=ierr, errmsg=cmessage) + if(ierr/=0)then; message1=trim(message1)//trim(cmessage)//':KWT routing state:'//trim(meta_mc(iVar)%varName); return; endif + do iens=1,nEns + do iSeg=1,nSeg + if (allocated(RFvec)) deallocate(RFvec, stat=ierr) + allocate(RFvec(numWaves(iens,iSeg)),stat=ierr); RFvec=0_i4b + where (RCHSTA_local(iSeg)%LKW_ROUTE%KWAVE(:)%RF) RFvec=1_i4b + array_3d_int(iSeg,1:numWaves(iens,iSeg),iens) = RFvec + array_3d_int(iSeg,numWaves(iens,iSeg)+1:,iens) = integerMissing + end do + end do + call write_pnetcdf(pioFileDescState, meta_kwt(iVar)%varName, array_3d_int, iodesc_wave_int, ierr, cmessage) + deallocate(array_3d_int) + case(ixKWT%tentry) + allocate(array_3d_dp(nSeg, nWave, nEns), stat=ierr, errmsg=cmessage) + if(ierr/=0)then; message1=trim(message1)//trim(cmessage)//':KWT routing state:'//trim(meta_mc(iVar)%varName); return; endif + do iens=1,nEns + do iSeg=1,nSeg + array_3d_dp(iSeg,1:numWaves(iens,iSeg),iens) = RCHSTA_local(iSeg)%LKW_ROUTE%KWAVE(:)%TI + array_3d_dp(iSeg,numWaves(iens,iSeg)+1:,iens) = realMissing + end do + end do + call write_pnetcdf(pioFileDescState, meta_kwt(iVar)%varName, array_3d_dp, iodesc_wave_double, ierr, cmessage) + deallocate(array_3d_dp) + case(ixKWT%texit) + allocate(array_3d_dp(nSeg, nWave, nEns), stat=ierr, errmsg=cmessage) + if(ierr/=0)then; message1=trim(message1)//trim(cmessage)//':KWT routing state:'//trim(meta_mc(iVar)%varName); return; endif + do iens=1,nEns + do iSeg=1,nSeg + array_3d_dp(iSeg,1:numWaves(iens,iSeg),iens) = RCHSTA_local(iSeg)%LKW_ROUTE%KWAVE(:)%TR + array_3d_dp(iSeg,numWaves(iens,iSeg)+1:,iens) = realMissing + end do + end do + call write_pnetcdf(pioFileDescState, meta_kwt(iVar)%varName, array_3d_dp, iodesc_wave_double, ierr, cmessage) + deallocate(array_3d_dp) + case(ixKWT%qwave) + allocate(array_3d_dp(nSeg, nWave, nEns), stat=ierr, errmsg=cmessage) + if(ierr/=0)then; message1=trim(message1)//trim(cmessage)//':KWT routing state:'//trim(meta_mc(iVar)%varName); return; endif + do iens=1,nEns + do iSeg=1,nSeg + array_3d_dp(iSeg,1:numWaves(iens,iSeg),iens) = RCHSTA_local(iSeg)%LKW_ROUTE%KWAVE(:)%QF + array_3d_dp(iSeg,numWaves(iens,iSeg)+1:,iens) = realMissing + end do + end do + call write_pnetcdf(pioFileDescState, meta_kwt(iVar)%varName, array_3d_dp, iodesc_wave_double, ierr, cmessage) + deallocate(array_3d_dp) + case(ixKWT%qwave_mod) + allocate(array_3d_dp(nSeg, nWave, nEns), stat=ierr, errmsg=cmessage) + if(ierr/=0)then; message1=trim(message1)//trim(cmessage)//':KWT routing state:'//trim(meta_mc(iVar)%varName); return; endif + do iens=1,nEns + do iSeg=1,nSeg + array_3d_dp(iSeg,1:numWaves(iens,iSeg),iens) = RCHSTA_local(iSeg)%LKW_ROUTE%KWAVE(:)%QM + array_3d_dp(iSeg,numWaves(iens,iSeg)+1:,iens) = realMissing + end do + end do + call write_pnetcdf(pioFileDescState, meta_kwt(iVar)%varName, array_3d_dp, iodesc_wave_double, ierr, cmessage) + deallocate(array_3d_dp) + case(ixKWT%vol) + allocate(array_2d_dp(nSeg, nEns), stat=ierr, errmsg=cmessage) + if(ierr/=0)then; message1=trim(message1)//trim(cmessage)//':KWT routing state:'//trim(meta_mc(iVar)%varName); return; endif + do iens=1,nEns + do iSeg=1,nSeg + array_2d_dp(iSeg,iens) = RCHFLX_local(iSeg)%ROUTE(idxKWT)%REACH_VOL(1) + end do + end do + call write_pnetcdf(pioFileDescState, meta_kwt(iVar)%varName, array_2d_dp, iodesc_rch_double, ierr, cmessage) + deallocate(array_2d_dp) + case default; ierr=20; message1=trim(message1)//'unable to identify KWT routing state variable index'; return + end select + if(ierr/=0)then; message1=trim(message1)//'problem allocating space for KWT routing state '//trim(meta_kwt(iVar)%varName); return; endif + end do + + end associate END SUBROUTINE write_KWT_state @@ -1227,7 +1270,8 @@ SUBROUTINE write_KW_state(ierr, message1) integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message1 ! error message ! local variables - type(states) :: state ! temporal state data structures -currently 2 river routing scheme + basin IRF routing + real(dp), allocatable :: array_2d_dp(:,:) + real(dp), allocatable :: array_3d_dp(:,:,:) integer(i4b) :: iVar,iens,iSeg ! index loops for variables, ensembles and segments respectively ierr=0; message1='write_KW_state/' @@ -1236,41 +1280,30 @@ SUBROUTINE write_KW_state(ierr, message1) nEns => meta_stateDims(ixStateDims%ens)%dimLength, & nMesh => meta_stateDims(ixStateDims%mol_kw)%dimLength) ! maximum waves allowed in a reach - allocate(state%var(nVarsKW), stat=ierr, errmsg=cmessage) - if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif - - do iVar=1,nVarsKW - select case(iVar) - case(ixKW%qsub); allocate(state%var(iVar)%array_3d_dp(nSeg, nMesh, nEns), stat=ierr) - case(ixKW%vol); allocate(state%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 KW routing state '//trim(meta_kw(iVar)%varName); return; endif - end do - - ! --Convert data structures to arrays - do iens=1,nEns - do iSeg=1,nSeg - do iVar=1,nVarsKW - select case(iVar) - case(ixKW%qsub); state%var(iVar)%array_3d_dp(iSeg,1:nMesh,iens) = RCHSTA_local(iSeg)%KW_ROUTE%molecule%Q(1:nMesh) - case(ixKW%vol); state%var(iVar)%array_2d_dp(iSeg,iens) = RCHFLX_local(iSeg)%ROUTE(idxKW)%REACH_VOL(1) - case default; ierr=20; message1=trim(message1)//'unable to identify KW routing state variable index'; return - end select - enddo ! variable loop - enddo ! seg loop - enddo ! ensemble loop - - ! Writing netCDF do iVar=1,nVarsKW select case(iVar) - case(ixKW%qsub) - call write_pnetcdf(pioFileDescState, meta_kw(iVar)%varName, state%var(iVar)%array_3d_dp, iodesc_mesh_kw_double, ierr, cmessage) - case(ixKW%vol) - call write_pnetcdf(pioFileDescState, meta_kw(iVar)%varName, state%var(iVar)%array_2d_dp, iodesc_state_double, ierr, cmessage) - case default; ierr=20; message1=trim(message1)//'unable to identify KW variable index for nc writing'; return + case(ixKW%qsub) + allocate(array_3d_dp(nSeg, nMesh, nEns), stat=ierr, errmsg=cmessage) + if(ierr/=0)then; message1=trim(message1)//trim(cmessage)//':KW routing state:'//trim(meta_mc(iVar)%varName); return; endif + do iens=1,nEns + do iSeg=1,nSeg + array_3d_dp(iSeg,1:nMesh,iens) = RCHSTA_local(iSeg)%KW_ROUTE%molecule%Q(1:nMesh) + end do + end do + call write_pnetcdf(pioFileDescState, meta_kw(iVar)%varName, array_3d_dp, iodesc_mesh_kw_double, ierr, cmessage) + deallocate(array_3d_dp) + case(ixKW%vol) + allocate(array_2d_dp(nSeg, nEns), stat=ierr, errmsg=cmessage) + if(ierr/=0)then; message1=trim(message1)//trim(cmessage)//':KW routing state:'//trim(meta_mc(iVar)%varName); return; endif + do iens=1,nEns + do iSeg=1,nSeg + array_2d_dp(iSeg,iens) = RCHFLX_local(iSeg)%ROUTE(idxKW)%REACH_VOL(1) + end do + end do + call write_pnetcdf(pioFileDescState, meta_kw(iVar)%varName, array_2d_dp, iodesc_rch_double, ierr, cmessage) + deallocate(array_2d_dp) + case default; ierr=20; message1=trim(message1)//'unable to identify KW routing state variable index'; return end select - if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif end do end associate @@ -1281,11 +1314,12 @@ SUBROUTINE write_MC_state(ierr, message1) USE globalData, ONLY: idxMC implicit none ! output - integer(i4b), intent(out) :: ierr ! error code - character(*), intent(out) :: message1 ! error message + integer(i4b), intent(out) :: ierr ! error code + character(*), intent(out) :: message1 ! error message ! local variables - type(states) :: state ! temporal state data structures -currently 2 river routing scheme + basin IRF routing - integer(i4b) :: iVar,iens,iSeg ! index loops for variables, ensembles and segments respectively + real(dp), allocatable :: array_2d_dp(:,:) + real(dp), allocatable :: array_3d_dp(:,:,:) + integer(i4b) :: iVar,iens,iSeg ! index loops for variables, ensembles and segments respectively ierr=0; message1='write_MC_state/' @@ -1293,41 +1327,30 @@ SUBROUTINE write_MC_state(ierr, message1) nEns => meta_stateDims(ixStateDims%ens)%dimLength, & nMesh => meta_stateDims(ixStateDims%mol_mc)%dimLength) ! maximum waves allowed in a reach - allocate(state%var(nVarsMC), stat=ierr, errmsg=cmessage) - if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif - do iVar=1,nVarsMC select case(iVar) - case(ixMC%qsub); allocate(state%var(iVar)%array_3d_dp(nSeg, nMesh, nEns), stat=ierr) - case(ixMC%vol); allocate(state%var(iVar)%array_2d_dp(nSeg, nEns), stat=ierr) - case default; ierr=20; message1=trim(message1)//'unable to identify variable index'; return + case(ixMC%qsub) + allocate(array_3d_dp(nSeg,nMesh,nEns), stat=ierr, errmsg=cmessage) + if(ierr/=0)then; message1=trim(message1)//trim(cmessage)//':MC routing state:'//trim(meta_mc(iVar)%varName); return; endif + do iens=1,nEns + do iSeg=1,nSeg + array_3d_dp(iSeg,1:nMesh,iens) = RCHSTA_local(iSeg)%MC_ROUTE%molecule%Q(1:nMesh) + end do + end do + call write_pnetcdf(pioFileDescState, meta_mc(iVar)%varName, array_3d_dp, iodesc_mesh_mc_double, ierr, cmessage) + deallocate(array_3d_dp) + case(ixMC%vol) + allocate(array_2d_dp(nSeg, nEns),stat=ierr,errmsg=cmessage) + if(ierr/=0)then; message1=trim(message1)//trim(cmessage)//':MC routing state:'//trim(meta_mc(iVar)%varName); return; endif + do iens=1,nEns + do iSeg=1,nSeg + array_2d_dp(iSeg,iens) = RCHFLX_local(iSeg)%ROUTE(idxMC)%REACH_VOL(1) + end do + end do + call write_pnetcdf(pioFileDescState, meta_mc(iVar)%varName, array_2d_dp, iodesc_rch_double, ierr, cmessage) + deallocate(array_2d_dp, stat=ierr) + case default; ierr=20; message1=trim(message1)//'unable to identify MC routing state variable index'; return end select - if(ierr/=0)then; message1=trim(message1)//'problem allocating space for MC routing state '//trim(meta_mc(iVar)%varName); return; endif - end do - - ! --Convert data structures to arrays - do iens=1,nEns - do iSeg=1,nSeg - do iVar=1,nVarsMC - select case(iVar) - case(ixMC%qsub); state%var(iVar)%array_3d_dp(iSeg,1:nMesh,iens) = RCHSTA_local(iSeg)%MC_ROUTE%molecule%Q(1:nMesh) - case(ixMC%vol); state%var(iVar)%array_2d_dp(iSeg,iens) = RCHFLX_local(iSeg)%ROUTE(idxMC)%REACH_VOL(1) - case default; ierr=20; message1=trim(message1)//'unable to identify MC routing state variable index'; return - end select - enddo ! variable loop - enddo ! seg loop - enddo ! ensemble loop - - ! Writing netCDF - do iVar=1,nVarsMC - select case(iVar) - case(ixMC%qsub) - call write_pnetcdf(pioFileDescState, meta_mc(iVar)%varName, state%var(iVar)%array_3d_dp, iodesc_mesh_mc_double, ierr, cmessage) - case(ixMC%vol) - call write_pnetcdf(pioFileDescState, meta_mc(iVar)%varName, state%var(iVar)%array_2d_dp, iodesc_state_double, ierr, cmessage) - case default; ierr=20; message1=trim(message1)//'unable to identify MC variable index for nc writing'; return - end select - if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif end do end associate @@ -1341,8 +1364,9 @@ SUBROUTINE write_DW_state(ierr, message1) integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message1 ! error message ! local variables - type(states) :: state ! temporal state data structures -currently 2 river routing scheme + basin IRF routing - integer(i4b) :: iVar,iens,iSeg ! index loops for variables, ensembles and segments respectively + real(dp), allocatable :: array_2d_dp(:,:) + real(dp), allocatable :: array_3d_dp(:,:,:) + integer(i4b) :: iVar,iens,iSeg ! index loops for variables, ensembles and segments respectively ierr=0; message1='write_DW_state/' @@ -1350,47 +1374,109 @@ SUBROUTINE write_DW_state(ierr, message1) nEns => meta_stateDims(ixStateDims%ens)%dimLength, & nMesh => meta_stateDims(ixStateDims%mol_dw)%dimLength) ! maximum waves allowed in a reach - allocate(state%var(nVarsDW), stat=ierr, errmsg=cmessage) - if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif - do iVar=1,nVarsDW select case(iVar) - case(ixDW%qsub); allocate(state%var(iVar)%array_3d_dp(nSeg, nMesh, nEns), stat=ierr) - case(ixDW%vol); allocate(state%var(iVar)%array_2d_dp(nSeg, nEns), stat=ierr) - case default; ierr=20; message1=trim(message1)//'unable to identify variable index'; return + case(ixDW%qsub) + allocate(array_3d_dp(nSeg,nMesh,nEns),stat=ierr,errmsg=cmessage) + if(ierr/=0)then; message1=trim(message1)//trim(cmessage)//':DW routing state:'//trim(meta_mc(iVar)%varName); return; endif + do iens=1,nEns + do iSeg=1,nSeg + array_3d_dp(iSeg,1:nMesh,iens) = RCHSTA_local(iSeg)%DW_ROUTE%molecule%Q(1:nMesh) + end do + end do + call write_pnetcdf(pioFileDescState, meta_dw(iVar)%varName, array_3d_dp, iodesc_mesh_dw_double, ierr, cmessage) + deallocate(array_3d_dp) + case(ixDW%vol) + allocate(array_2d_dp(nSeg, nEns),stat=ierr,errmsg=cmessage) + if(ierr/=0)then; message1=trim(message1)//trim(cmessage)//':DW routing state:'//trim(meta_mc(iVar)%varName); return; endif + do iens=1,nEns + do iSeg=1,nSeg + array_2d_dp(iSeg,iens) = RCHFLX_local(iSeg)%ROUTE(idxDW)%REACH_VOL(1) + end do + end do + call write_pnetcdf(pioFileDescState, meta_dw(iVar)%varName, array_2d_dp, iodesc_rch_double, ierr, cmessage) + deallocate(array_2d_dp) + case default; ierr=20; message1=trim(message1)//'unable to identify DW routing state variable index'; return end select - if(ierr/=0)then; message1=trim(message1)//'problem allocating space for DW routing state '//trim(meta_dw(iVar)%varName); return; endif - end do - - ! --Convert data structures to arrays - do iens=1,nEns - do iSeg=1,nSeg - do iVar=1,nVarsDW - select case(iVar) - case(ixDW%qsub); state%var(iVar)%array_3d_dp(iSeg,1:nMesh,iens) = RCHSTA_local(iSeg)%DW_ROUTE%molecule%Q(1:nMesh) - case(ixDW%vol); state%var(iVar)%array_2d_dp(iSeg,iens) = RCHFLX_local(iSeg)%ROUTE(idxDW)%REACH_VOL(1) - case default; ierr=20; message1=trim(message1)//'unable to identify DW routing state variable index'; return - end select - enddo ! variable loop - enddo ! seg loop - enddo ! ensemble loop - - ! Writing netCDF - do iVar=1,nVarsDW - select case(iVar) - case(ixDW%qsub) - call write_pnetcdf(pioFileDescState, meta_dw(iVar)%varName, state%var(iVar)%array_3d_dp, iodesc_mesh_dw_double, ierr, cmessage) - case(ixDW%vol) - call write_pnetcdf(pioFileDescState, meta_dw(iVar)%varName, state%var(iVar)%array_2d_dp, iodesc_state_double, ierr, cmessage) - case default; ierr=20; message1=trim(message1)//'unable to identify DW variable index for nc writing'; return - end select - if(ierr/=0)then; message1=trim(message1)//trim(cmessage); return; endif - end do + enddo ! variable loop end associate END SUBROUTINE write_DW_state + SUBROUTINE write_history_state(ierr, message1) + implicit none + ! output + integer(i4b), intent(out) :: ierr ! error code + character(*), intent(out) :: message1 ! error message + ! local variables + real(sp), allocatable :: array_sp(:) + + ierr=0; message1='write_history_state/' + + allocate(array_sp(hVars%nRch),stat=ierr, errmsg=cmessage) + + call write_scalar_netcdf(pioFileDescState, 'nt', hVars%nt, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + + if (meta_hflx(ixHFLX%basRunoff)%varFile) then + call write_pnetcdf(pioFileDescState, 'basRunoff', hVars%basRunoff, ioDesc_hru_float, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + end if + + if (meta_rflx(ixRFLX%instRunoff)%varFile) then + call write_pnetcdf(pioFileDescState, 'instRunoff', hVars%instRunoff, ioDesc_rch_float, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + endif + + if (meta_rflx(ixRFLX%dlayRunoff)%varFile) then + call write_pnetcdf(pioFileDescState, 'dlayRunoff', hVars%dlayRunoff, ioDesc_rch_float, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + endif + + if (meta_rflx(ixRFLX%sumUpstreamRunoff)%varFile) then + array_sp = hVars%discharge(:, idxSUM) + call write_pnetcdf(pioFileDescState, 'sumUpstreamRunoff', array_sp, ioDesc_rch_float, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + endif + + if (meta_rflx(ixRFLX%KWTroutedRunoff)%varFile) then + array_sp = hVars%discharge(:, idxKWT) + call write_pnetcdf(pioFileDescState, 'KWTroutedRunoff', array_sp, ioDesc_rch_float, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + endif + + if (meta_rflx(ixRFLX%IRFroutedRunoff)%varFile) then + array_sp = hVars%discharge(:, idxIRF) + call write_pnetcdf(pioFileDescState, 'IRFroutedRunoff', array_sp, ioDesc_rch_float, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + endif + + if (meta_rflx(ixRFLX%KWroutedRunoff)%varFile) then + array_sp = hVars%discharge(:, idxKW) + call write_pnetcdf(pioFileDescState, 'KWroutedRunoff', array_sp, ioDesc_rch_float, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + endif + + if (meta_rflx(ixRFLX%MCroutedRunoff)%varFile) then + array_sp = hVars%discharge(:, idxMC) + call write_pnetcdf(pioFileDescState, 'MCroutedRunoff', array_sp, ioDesc_rch_float, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + endif + + if (meta_rflx(ixRFLX%DWroutedRunoff)%varFile) then + array_sp = hVars%discharge(:, idxDW) + call write_pnetcdf(pioFileDescState, 'DWroutedRunoff', array_sp, ioDesc_rch_float, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + endif + + if (meta_rflx(ixRFLX%IRFvolume)%varFile) then + array_sp = hVars%volume(:, idxIRF) + call write_pnetcdf(pioFileDescState, 'IRFvolume', array_sp, ioDesc_rch_float, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + endif + + END SUBROUTINE write_history_state END SUBROUTINE write_state_nc diff --git a/route/build/src/write_simoutput_pio.f90 b/route/build/src/write_simoutput_pio.f90 index 055d26d2..1c662bd8 100644 --- a/route/build/src/write_simoutput_pio.f90 +++ b/route/build/src/write_simoutput_pio.f90 @@ -12,20 +12,25 @@ MODULE write_simoutput_pio USE globalData, ONLY: masterproc USE globalData, ONLY: hfileout, hfileout_gage, rfileout USE historyFile, ONLY: histFile +USE histVars_data, ONLY: histVars USE io_rpointfile, ONLY: io_rpfile USE ascii_utils, ONLY: lower USE pio_utils implicit none +logical(lgt), save :: initHvars=.false. ! +type(histVars), save :: hVars ! type(histFile), save :: hist_all_network ! type(histFile), save :: hist_gage ! integer(i4b), allocatable, save :: index_write_gage(:) ! private +public::hVars public::main_new_file public::output public::close_all public::init_histFile +public::get_compdof_all_network ! used in write_restart_pio CONTAINS @@ -57,7 +62,7 @@ SUBROUTINE main_new_file(ierr, message) ierr=0; message='main_new_file/' - createNewFile = newFileAlarm(simDatetime, newFileFrequency, ierr, cmessage) + createNewFile = newFileAlarm(simDatetime(0:1), newFileFrequency, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif if (createNewFile) then @@ -111,8 +116,6 @@ SUBROUTINE main_new_file(ierr, message) end if - simDatetime(0) = simDatetime(1) - ! update history files call io_rpfile('w', ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif @@ -154,6 +157,9 @@ END FUNCTION newFileAlarm SUBROUTINE output(ierr, message) USE public_var, ONLY: outputAtGage ! ascii containing last restart and history files + USE public_var, ONLY: nOutFreq ! + USE public_var, ONLY: outputFrequency ! + USE globalData, ONLY: simDatetime ! previous,current and next model datetime USE globalData, ONLY: timeVar ! current simulation time variable USE globalData, ONLY: RCHFLX_trib ! reach flux data structure containing current flux variables USE globalData, ONLY: rch_per_proc ! number of reaches assigned to each proc (size = num of procs+1) @@ -166,39 +172,142 @@ SUBROUTINE output(ierr, message) integer(i4b), intent(out) :: ierr ! error code character(*), intent(out) :: message ! error message ! local variables: + logical(lgt) :: writeHistory ! integer(i4b) :: nRch_local ! number of reaches per processors - integer(i4b) :: nRch_root ! number of reaches in root processors (including halo reaches) + real(dp), allocatable :: basinRunoff(:) integer(i4b), allocatable :: index_write_all(:) ! indices in RCHFLX_trib to be written in netcdf character(strLen) :: cmessage ! error message of downwind routine ierr=0; message='output/' - ! compute index array for each processors (herer - if (masterproc) then - nRch_local = sum(rch_per_proc(-1:pid)) - nRch_root = nRch_mainstem+nTribOutlet+rch_per_proc(0) - allocate(index_write_all(nRch_local)) - if (nRch_mainstem>0) then - index_write_all(1:nRch_mainstem) = arth(1,1,nRch_mainstem) - end if - index_write_all(nRch_mainstem+1:nRch_local) = arth(nRch_mainstem+nTribOutlet+1, 1, rch_per_proc(0)) - else - nRch_local = rch_per_proc(pid) - allocate(index_write_all(nRch_local)) - index_write_all = arth(1,1,nRch_local) + ! initialize if histVars data - hVar is not initalized + if (.not. initHvars) then + call init_histVar_data(ierr, message) + initHvars = .true. end if - ! write out output variables in history files - call hist_all_network%write_flux(timeVar, RCHFLX_trib, index_write_all, ierr, cmessage) + ! get unified basinRunoff array + call get_proc_flux(ierr, cmessage, basinRunoff=basinRunoff) + + ! accumulate history variables + call hVars%aggregate(timeVar, basinRunoff, RCHFLX_trib, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif - if (outputAtGage) then - call hist_gage%write_flux(timeVar, RCHFLX_trib, index_write_gage, ierr, cmessage) + ! history output alarm + if (nOutFreq/=integerMissing) then + writeHistory = writeAlarmNumeric(nOutFreq) + else + writeHistory = writeAlarmLiteral(simDatetime(1:2), outputFrequency, ierr, cmessage) if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif end if + if (writeHistory) then + + ! compute final output values (compute mean) + call hVars%finalize() + + ! compute index array for each processors (herer + if (masterproc) then + nRch_local = sum(rch_per_proc(-1:pid)) + allocate(index_write_all(nRch_local)) + if (nRch_mainstem>0) then + index_write_all(1:nRch_mainstem) = arth(1,1,nRch_mainstem) + end if + index_write_all(nRch_mainstem+1:nRch_local) = arth(nRch_mainstem+nTribOutlet+1, 1, rch_per_proc(0)) + else + nRch_local = rch_per_proc(pid) + allocate(index_write_all(nRch_local)) + index_write_all = arth(1,1,nRch_local) + end if + + ! write out output variables in history files + call hist_all_network%write_flux(hVars, index_write_all, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + + if (outputAtGage) then + call hist_gage%write_flux(hVars, index_write_gage, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + end if + + ! refresh history output buffers + call hVars%refresh() + + end if + END SUBROUTINE output + ! ********************************************************************* + ! private subroutine: restart alarming + ! ********************************************************************* + logical(lgt) FUNCTION writeAlarmNumeric(alarmFrequency) + + implicit none + ! Argument variables + integer(i4b), intent(in) :: alarmFrequency ! numeric output frequency + + writeAlarmNumeric = .false. + if (hVars%nt == alarmFrequency) writeAlarmNumeric = .true. + + END FUNCTION writeAlarmNumeric + + ! ********************************************************************* + ! private subroutine: restart alarming + ! ********************************************************************* + logical(lgt) FUNCTION writeAlarmLiteral(inDatetime, alarmFrequency, ierr, message) + + implicit none + ! Argument variables + type(datetime), intent(in) :: inDatetime(1:2) ! datetime at previous and current timestep + character(*), intent(in) :: alarmFrequency ! literal output frequency + integer(i4b), intent(out) :: ierr ! error code + character(*), intent(out) :: message ! error message + + writeAlarmLiteral = .false. + select case(lower(trim(alarmFrequency))) + case('yearly'); writeAlarmLiteral=(inDatetime(1)%year() /=inDatetime(2)%year()) + case('monthly');writeAlarmLiteral=(inDatetime(1)%month()/=inDatetime(2)%month()) + case('daily'); writeAlarmLiteral=(inDatetime(1)%day() /=inDatetime(2)%day()) + case default; ierr=20; message=trim(message)//'unable to identify the option to define new output files'; return + end select + + END FUNCTION writeAlarmLiteral + + ! ********************************************************************* + ! private subroutine - organize flux data array or structure into per processor + ! ********************************************************************* + SUBROUTINE get_proc_flux(ierr, message, basinRunoff) + + USE globalData, ONLY: nHRU_mainstem ! number of mainstem HRUs + USE globalData, ONLY: basinRunoff_main ! mainstem only HRU runoff + USE globalData, ONLY: basinRunoff_trib ! tributary only HRU runoff + USE globalData, ONLY: hru_per_proc ! number of hrus assigned to each proc (size = num of procs+1) + + implicit none + ! Argument variables + real(dp), allocatable, optional, intent(out) :: basinRunoff(:) + integer(i4b), intent(out) :: ierr ! error code + character(*), intent(out) :: message ! error message + ! local variables + integer(i4b) :: nHRU_local + character(strLen) :: cmessage ! error message of downwind routine + + if (present(basinRunoff)) then + if (masterproc) then + nHRU_local = nHRU_mainstem + hru_per_proc(0) + allocate(basinRunoff(nHRU_local), stat=ierr, errmsg=cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [basinRunoff]'; return; endif + if (nHRU_mainstem>0) basinRunoff(1:nHRU_mainstem) = basinRunoff_main(1:nHRU_mainstem) + if (hru_per_proc(0)>0) basinRunoff(nHRU_mainstem+1:nHRU_local) = basinRunoff_trib(:) + else + nHRU_local = hru_per_proc(pid) + allocate(basinRunoff(nHRU_local), stat=ierr, errmsg=cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage)//' [basinRunoff]'; return; endif + basinRunoff = basinRunoff_trib + endif + end if + + END SUBROUTINE get_proc_flux + ! ********************************************************************* ! private subroutine: get pio local-global mapping data ! ********************************************************************* @@ -419,4 +528,36 @@ SUBROUTINE init_histFile(ierr, message) END SUBROUTINE init_histFile + ! ********************************************************************* + ! public subroutine: initialize history variable buffers + ! ********************************************************************* + SUBROUTINE init_histVar_data(ierr, message) ! output: error control + + USE globalData, ONLY: nRch_mainstem, nRch_trib + USE globalData, ONLY: nHRU_mainstem, nHRU_trib + USE globalData, ONLY: nTribOutlet + + implicit none + ! Argument variables: + integer(i4b), intent(out) :: ierr ! error code + character(*), intent(out) :: message ! error message + ! local variable + character(len=strLen) :: cmessage ! error message from subroutine + integer(i4b) :: nHru_local ! nHRU in each processor + integer(i4b) :: nRch_local ! nRch in each processor + + ierr=0; message='init_histVar_data/' + + if (masterproc) then + nHru_local=nHru_mainstem+nHru_trib + nRch_local=nRch_mainstem+nRch_trib+nTribOutlet + else + nHru_local=nHru_trib + nRch_local=nRch_trib + end if + hVars = histVars(nHru_local, nRch_local, ierr, cmessage) + if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif + + END SUBROUTINE init_histVar_data + END MODULE write_simoutput_pio