Skip to content

Enable hourly timestep #31

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
CyrilThebault wants to merge 1 commit into from
Closed

Enable hourly timestep #31

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 16 additions & 1 deletion build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ PROGRAM DISTRIBUTED_DRIVER
FMODEL_ID,&
suffix_forcing,suffix_elev_bands,&
numtim_sub_str,&
KSTOP_str, MAXN_str, PCENTO_str
KSTOP_str, MAXN_str, PCENTO_str, SEED_str

! data modules
USE model_defn,nstateFUSE=>nstate ! model definition structures
Expand Down Expand Up @@ -125,6 +125,7 @@ PROGRAM DISTRIBUTED_DRIVER
REAL(SP), DIMENSION(:), ALLOCATABLE :: BU ! vector of upper parameter bounds
REAL(SP), DIMENSION(:), ALLOCATABLE :: APAR ! model parameter set
INTEGER(KIND=4) :: ISEED ! seed for the random sequence
INTEGER :: count ! used to create a random value based on clock time
REAL(KIND=4),DIMENSION(:), ALLOCATABLE :: URAND ! vector of quasi-random numbers U[0,1]
REAL(SP) :: METRIC_VAL ! error from the simulation

Expand Down Expand Up @@ -349,6 +350,20 @@ PROGRAM DISTRIBUTED_DRIVER
READ (MAXN_STR,*) MAXN ! maximum number of trials before optimization is terminated
READ (KSTOP_STR,*) KSTOP ! number of shuffling loops the value must change by PCENTO (MAX=9)
READ (PCENTO_STR,*) PCENTO ! the percentage


! SCE random seed.
! SEED_STR = -9999 means use random seed.
! SEED_STR > 0 means use fixed user-defined seed.
READ (SEED_STR,*) ISEED

IF (ISEED <= 0) THEN
CALL SYSTEM_CLOCK(count)
ISEED = ABS(count)
ENDIF

PRINT *, 'SCE seed = ', ISEED


PRINT *, 'SCE parameters read from file manager:'
PRINT *, 'Maximum number of trials before SCE optimization is stopped (MAXN) = ', MAXN_STR
Expand Down
4 changes: 2 additions & 2 deletions build/FUSE_SRC/FUSE_ENGINE/force_info.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@ SUBROUTINE force_info(fuse_mode,ierr,message)
USE multiforce,only:vname_dtime ! name of time variable (time since reference time)
USE multiforce,only:deltim ! model timestep (days)
USE multiforce,only:istart,numtim_sim ! index for start of inference, and number steps in the reduced array
USE multiforce,only:amult_ppt,amult_pet,amult_q ! used to convert fluxes to mm/day
USE multiforce,only:amult_ppt,amult_pet,amult_q ! used to convert input flux units to FUSE internal units
USE multiforce,only:numtim_sub ! number of time steps of subperiod (will be kept in memory)

IMPLICIT NONE
Expand Down Expand Up @@ -229,7 +229,7 @@ subroutine get_multiplier(cunits, amult, ierr, message)
real(sp),parameter :: hrprday=24._sp ! number of hours per day
! initialize error control
ierr=0; message='get_multiplier/'
! if units are undefined, assume mm/day and have an early return
! if units are undefined, assume no unit conversion and have an early return
if(trim(cunits)=='undefined')then; amult=1._sp; return; endif
! find the position of the "/" character
ipos = index(trim(cunits),'/')
Expand Down
17 changes: 16 additions & 1 deletion build/FUSE_SRC/FUSE_ENGINE/get_time_indices.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,12 +19,16 @@ SUBROUTINE GET_TIME_INDICES
USE multiforce, only: sim_beg,sim_end ! timestep indices
USE multiforce, only: eval_beg,eval_end ! timestep indices
USE multiforce, only: SUB_PERIODS_FLAG ! .true. if subperiods are used to run FUSE

USE time_io, only: time_units_to_days, INPUT_DT_DAYS ! Add by Cyril Thebault

USE fuse_fileManager,only:date_start_sim,date_end_sim,&
date_start_eval,date_end_eval,&
numtim_sub_str

real(sp) :: jdate_ref_netcdf
real(sp) :: time_factor_days ! Conversion factor from NetCDF time units to days.

INTEGER(I4B) :: ERR ! error code
CHARACTER(LEN=1024) :: MESSAGE ! error message

Expand All @@ -43,7 +47,18 @@ SUBROUTINE GET_TIME_INDICES
call juldayss(iy,im,id,ih, & ! convert it to julian day
jdate_ref_netcdf,err,message)

julian_day_input=jdate_ref_netcdf+time_steps ! julian day of each time step of the input file
time_factor_days = time_units_to_days(timeUnits) ! Add by Cyril Thebault for the hourly time step
julian_day_input = jdate_ref_netcdf + time_steps * time_factor_days ! julian day of each time step of the input file

! Store the physical duration of one forcing timestep in days. Add by Cyril Thebault
! This assumes a regular input timestep.
IF (SIZE(time_steps) > 1) THEN
INPUT_DT_DAYS = (time_steps(2) - time_steps(1)) * time_factor_days
ELSE
INPUT_DT_DAYS = time_factor_days
ENDIF

PRINT *, 'INPUT_DT_DAYS = ', INPUT_DT_DAYS

call caldatss(julian_day_input(1),iy,im,id,ih,imin,isec)
print *, 'Start date input file=',iy,im,id
Expand Down
2 changes: 1 addition & 1 deletion build/FUSE_SRC/FUSE_ENGINE/mean_stats.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@ SUBROUTINE MEAN_STATS()

! extract OBS and SIM for evaluation period, note that sim_beg, eval_beg, sim_end,
! eval_end are all with respect to julian_day_input
QSIM = AROUTE_3d(1,1,eval_beg-sim_beg+1:eval_end-sim_beg+1)%Q_ROUTED
QSIM = AROUTE_3d(1,1,eval_beg-sim_beg+1:eval_end-sim_beg+1)%Q_ROUTED
QOBS = aValid(1,1,eval_beg-sim_beg+1:eval_end-sim_beg+1)%OBSQ

! check for missing QOBS values
Expand Down
4 changes: 2 additions & 2 deletions build/FUSE_SRC/FUSE_ENGINE/metaoutput.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,8 +27,8 @@ MODULE metaoutput
SUBROUTINE VARDESCRIBE()
I=0 ! initialize counter
! model forcing
I=I+1; VNAME(I)='ppt '; LNAME(I)='precipitation rate '; VUNIT(I)='mm day-1 '
I=I+1; VNAME(I)='pet '; LNAME(I)='potential evapotranspiration rate '; VUNIT(I)='mm day-1 '
I=I+1; VNAME(I)='ppt '; LNAME(I)='precipitation rate '; VUNIT(I)='mm timestep-1 '
I=I+1; VNAME(I)='pet '; LNAME(I)='potential evapotranspiration rate '; VUNIT(I)='mm timestep-1 '
I=I+1; VNAME(I)='temp '; LNAME(I)='mean air temperature '; VUNIT(I)='deg.C '
I=I+1; VNAME(I)='obsq '; LNAME(I)='observed runoff '; VUNIT(I)='mm timestep-1'
! model states
Expand Down
14 changes: 7 additions & 7 deletions build/FUSE_SRC/FUSE_ENGINE/multiforce.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ MODULE multiforce
ENDTYPE TDATA
! the response structure (will not have a spatial dimension)
TYPE VDATA
REAL(SP) :: OBSQ ! observed runoff (mm day-1)
REAL(SP) :: OBSQ ! observed runoff (mm per timestep)
END TYPE VDATA
! ancillary forcing variables used to compute ET (will have a spatial dimension)
TYPE ADATA
Expand All @@ -34,9 +34,9 @@ MODULE multiforce
END TYPE ADATA
! the forcing data structure (will have a spatial dimension)
TYPE FDATA
REAL(SP) :: PPT ! water input: rain + melt (mm day-1)
REAL(SP) :: PPT ! water input: rain + melt (mm per timestep)
REAL(SP) :: TEMP ! temperature for snow model (deg.C)
REAL(SP) :: PET ! energy input: potential ET (mm day-1)
REAL(SP) :: PET ! energy input: potential ET (mm per timestep)
ENDTYPE FDATA
! --------------------------------------------------------------------------------------
! general
Expand Down Expand Up @@ -74,7 +74,7 @@ MODULE multiforce

INTEGER(i4b) :: istart=-1 ! index for start of inference period (in reduced array)
REAL(sp) :: jdayRef ! reference time (days)
REAL(sp) :: deltim=-1._dp ! length of time step (days)
REAL(sp) :: deltim=-1._dp ! length of time step

LOGICAL(LGT) :: SUB_PERIODS_FLAG ! .true. if subperiods are used to run FUSE

Expand Down Expand Up @@ -148,9 +148,9 @@ MODULE multiforce
INTEGER(i4b) :: ivarid_q=-1 ! variable ID for runoff

! multipliers for variables to convert fluxes to mm/day
REAL(sp) :: amult_ppt=-1._dp ! convert precipitation to mm/day
REAL(sp) :: amult_pet=-1._dp ! convert potential ET to mm/day
REAL(sp) :: amult_q=-1._dp ! convert runoff to mm/day
REAL(sp) :: amult_ppt=-1._dp ! precipitation unit multiplier
REAL(sp) :: amult_pet=-1._dp ! potential ET unit multiplier
REAL(sp) :: amult_q=-1._dp ! runoff unit multiplier

! missing values
INTEGER(I4B),PARAMETER :: NA_VALUE=-9999 ! integer designating missing values - TODO: retrieve from NetCDF file
Expand Down
2 changes: 1 addition & 1 deletion build/FUSE_SRC/FUSE_ENGINE/multiroute.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ MODULE multiroute
TYPE RUNOFF
REAL(SP) :: Q_INSTNT ! instantaneous runoff
REAL(SP) :: Q_ROUTED ! routed runoff
REAL(SP) :: Q_ACCURATE ! "accurate" runoff estimate (mm day-1)
REAL(SP) :: Q_ACCURATE ! "accurate" runoff estimate (same units as Q_ROUTED)
END TYPE RUNOFF
REAL(SP), DIMENSION(NTDH_MAX) :: FUTURE ! runoff placed in future time steps
TYPE(RUNOFF), DIMENSION(:), POINTER :: AROUTE ! runoff for all time steps
Expand Down
10 changes: 6 additions & 4 deletions build/FUSE_SRC/FUSE_ENGINE/varextract.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@ PURE FUNCTION VAREXTRACT(VARNAME)
USE multibands ! model snow bands
USE multiroute ! routed runoff
USE model_numerix ! model numerix parameters
USE time_io, ONLY: INPUT_DT_DAYS ! physical duration of one input timestep in days. Add by Cyril Thebault
IMPLICIT NONE
! input
CHARACTER(*), INTENT(IN) :: VARNAME ! variable name
Expand Down Expand Up @@ -234,8 +235,8 @@ PURE FUNCTION VAREXTRACT(VARNAME)
! time check
CASE ('chk_time') ; XVAR = W_FLUX%CHK_TIME
! extract model runoff
CASE ('q_instnt') ; XVAR = MROUTE%Q_INSTNT
CASE ('q_routed') ; XVAR = MROUTE%Q_ROUTED
CASE ('q_instnt') ; XVAR = MROUTE%Q_INSTNT * INPUT_DT_DAYS ! Modified by Cyril Thebault
CASE ('q_routed') ; XVAR = MROUTE%Q_ROUTED * INPUT_DT_DAYS ! Modified by Cyril Thebault
! extract information on numerical solution (shared in MODULE model_numerix)
CASE ('num_funcs') ; XVAR = NUM_FUNCS
CASE ('numjacobian'); XVAR = NUM_JACOBIAN
Expand Down Expand Up @@ -269,6 +270,7 @@ PURE FUNCTION VAREXTRACT_3d(VARNAME,numtim)
USE multibands ! model snow bands
USE multiroute ! routed runoff
USE model_numerix ! model numerix parameters
USE time_io, ONLY: INPUT_DT_DAYS ! physical duration of one input timestep in days. Add by Cyril Thebault
IMPLICIT NONE
! input
CHARACTER(*), INTENT(IN) :: VARNAME ! variable name
Expand Down Expand Up @@ -488,8 +490,8 @@ PURE FUNCTION VAREXTRACT_3d(VARNAME,numtim)
! time check
CASE ('chk_time') ; XVAR_3d = W_FLUX_3d%CHK_TIME
! extract model runoff
CASE ('q_instnt') ; XVAR_3d = AROUTE_3d%Q_INSTNT
CASE ('q_routed') ; XVAR_3d = AROUTE_3d%Q_ROUTED
CASE ('q_instnt') ; XVAR_3d = AROUTE_3d%Q_INSTNT * INPUT_DT_DAYS ! Modified by Cyril Thebault
CASE ('q_routed') ; XVAR_3d = AROUTE_3d%Q_ROUTED * INPUT_DT_DAYS ! Modified by Cyril Thebault
! extract information on numerical solution (shared in MODULE model_numerix)
CASE ('num_funcs') ; XVAR_3d = NUM_FUNCS
CASE ('numjacobian'); XVAR_3d = NUM_JACOBIAN
Expand Down
2 changes: 2 additions & 0 deletions build/FUSE_SRC/FUSE_HOOK/fuse_fileManager.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,7 @@ MODULE fuse_filemanager
CHARACTER(len=20) :: KSTOP_str ! number of shuffling loops the value must change by PCENTO
CHARACTER(len=20) :: MAXN_str ! maximum number of trials before optimization is terminated
CHARACTER(len=20) :: PCENTO_str ! the percentage
CHARACTER(len=20) :: SEED_str ! SCE seed; -9999 keeps legacy seed loop. Add by Cyril Thebault

!----------------------------------------------------
contains
Expand Down Expand Up @@ -170,6 +171,7 @@ subroutine fuse_SetDirsUndPhiles(fuseMusterDirektorIn,fuseFileManagerIn,err,mess
read(unt,*)MAXN_STR
read(unt,*)KSTOP_STR
read(unt,*)PCENTO_STR
read(unt,*)SEED_STR
close(unt)

! Convert Q_ONLY to logical
Expand Down
12 changes: 7 additions & 5 deletions build/FUSE_SRC/FUSE_NETCDF/get_gforce.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -306,6 +306,7 @@ SUBROUTINE get_gforce(itim,ncid_forc,ierr,message)
! MODULE multiforce -- populate structure GFORCE(*,*)%(*)
! ---------------------------------------------------------------------------------------
USE fuse_fileManager,only:INPUT_PATH ! defines data directory

USE multiforce,only:forcefile ! name of forcing file
USE multiforce,only:vname_aprecip ! variable name: precipitation
USE multiforce,only:vname_airtemp ! variable name: temperature
Expand All @@ -323,7 +324,7 @@ SUBROUTINE get_gforce(itim,ncid_forc,ierr,message)

USE multiforce,only:nspat1,nspat2,startSpat2 ! dimension lengths
USE multiforce,only:ncid_var ! NetCDF ID for forcing variables
USE multiforce,only:amult_ppt,amult_pet ! multipliers o convert to mm/day
USE multiforce,only:amult_ppt,amult_pet ! multipliers to convert input flux units to FUSE internal units
USE multiforce,only:gForce ! gridded forcing data
USE multiforce,only:ancilF ! ancillary forcing data
USE multiforce,only:nForce ! number of forcing variables
Expand Down Expand Up @@ -375,7 +376,7 @@ SUBROUTINE get_gforce(itim,ncid_forc,ierr,message)
ierr = nf90_get_var(ncid_forc, ncid_var(ivar), gTemp, start=(/1,startSpat2,iTim/), count=(/nSpat1,nSpat2,1/)); CALL HANDLE_ERR(IERR)
if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif

! save the data in the structure -- and convert fluxes to mm/day
! save the data in the structure -- and apply input unit multipliers
if(trim(cVec(iVar)%vname) == trim(vname_aprecip) )then; gForce(:,:)%ppt = gTemp(:,:,1)*amult_ppt; lCheck(ilook_aprecip) = .true.; endif
if(trim(cVec(iVar)%vname) == trim(vname_potevap) )then; gForce(:,:)%pet = gTemp(:,:,1)*amult_pet; lCheck(ilook_potevap) = .true.; endif
if(trim(cVec(iVar)%vname) == trim(vname_airtemp) )then; gForce(:,:)%temp = gTemp(:,:,1); lCheck(ilook_airtemp) = .true.; endif
Expand Down Expand Up @@ -411,6 +412,7 @@ SUBROUTINE get_gforce_3d(itim_start,numtim,ncid_forc,ierr,message)
! MODULE multiforce -- populate structure GFORCE_3d(*,*)%(*)
! ---------------------------------------------------------------------------------------
USE fuse_fileManager,only:INPUT_PATH ! defines data directory

USE multiforce,only:forcefile ! name of forcing file
USE multiforce,only:vname_aprecip ! variable name: precipitation
USE multiforce,only:vname_airtemp ! variable name: temperature
Expand All @@ -430,7 +432,7 @@ SUBROUTINE get_gforce_3d(itim_start,numtim,ncid_forc,ierr,message)

USE multiforce,only:nspat1,nspat2,startSpat2 ! dimension lengths
USE multiforce,only:ncid_var ! NetCDF ID for forcing variables
USE multiforce,only:amult_ppt,amult_pet ! multipliers o convert to mm/day
USE multiforce,only:amult_ppt,amult_pet, amult_q ! multipliers to convert input flux units to FUSE internal units
USE multiforce,only:gForce_3d ! gridded forcing data
USE multiforce,only:ancilF_3d ! ancillary forcing data
USE multiforce,only:nForce, nInput ! number of forcing variables
Expand Down Expand Up @@ -483,7 +485,7 @@ SUBROUTINE get_gforce_3d(itim_start,numtim,ncid_forc,ierr,message)
ierr = nf90_get_var(ncid_forc, ncid_var(ivar), gTemp, start=(/1,startSpat2,itim_start/), count=(/nSpat1,nSpat2,numtim/)); CALL HANDLE_ERR(IERR)
if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif

! save the data in the structure -- and convert fluxes to mm/day
! save the data in the structure -- and apply input unit multipliers
if(trim(cVec(iVar)%vname) == trim(vname_aprecip) )then

gForce_3d(:,:,1:numtim)%ppt = gTemp(:,:,:)*amult_ppt; lCheck(ilook_aprecip) = .true.
Expand All @@ -499,7 +501,7 @@ SUBROUTINE get_gforce_3d(itim_start,numtim,ncid_forc,ierr,message)
endif

if(trim(cVec(iVar)%vname) == trim(vname_q) )then
aValid(:,:,1:numtim)%obsq = gTemp(:,:,:); lCheck(ilook_q) = .true.
aValid(:,:,1:numtim)%obsq = gTemp(:,:,:) * amult_q; lCheck(ilook_q) = .true. ! Modified by Cyril Thebault
endif

! save the other variables required to compute PET
Expand Down
49 changes: 48 additions & 1 deletion build/FUSE_SRC/FUSE_TIME/time_io.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,51 @@ module time_io

public::get_modtim

! Physical duration of one input timestep in days. Add by Cyril Thebault.
! daily = 1.0
! hourly = 1.0 / 24.0
! Used to convert historical day-based process rates to mm/timestep.
REAL(SP), SAVE :: INPUT_DT_DAYS = 1.0_SP

contains


! ---------------------------------------------------------------------------------------
!
! Add by Cyril Thebault
!
! ---------------------------------------------------------------------------------------

REAL(SP) FUNCTION time_units_to_days(timeUnits) ! Converts NetCDF time units to days per time unit.
CHARACTER(*), INTENT(IN) :: timeUnits
CHARACTER(LEN=1024) :: units_lc
INTEGER(I4B) :: i

units_lc = ADJUSTL(timeUnits)

! lower-case simple
DO i=1,LEN_TRIM(units_lc)
IF (units_lc(i:i) >= 'A' .AND. units_lc(i:i) <= 'Z') &
units_lc(i:i) = ACHAR(IACHAR(units_lc(i:i)) + 32)
END DO

IF (INDEX(units_lc,'seconds since') > 0) THEN
time_units_to_days = 1.0_SP / 86400.0_SP
ELSEIF (INDEX(units_lc,'minutes since') > 0) THEN
time_units_to_days = 1.0_SP / 1440.0_SP
ELSEIF (INDEX(units_lc,'hours since') > 0) THEN
time_units_to_days = 1.0_SP / 24.0_SP
ELSEIF (INDEX(units_lc,'days since') > 0) THEN
time_units_to_days = 1.0_SP
ELSE
PRINT *, 'Error: unsupported NetCDF time units: ', TRIM(timeUnits)
STOP
ENDIF

END FUNCTION time_units_to_days




SUBROUTINE get_modtim(itim,ncid,ierr,message)
! ---------------------------------------------------------------------------------------
Expand All @@ -21,6 +65,9 @@ SUBROUTINE get_modtim(itim,ncid,ierr,message)
! Modules Modified:
! -----------------
! MODULE multiforce -- populate structure timDat%(*)
!
! Modified to work at the hourly time step by Cyril Thebault
!
! ---------------------------------------------------------------------------------------
USE fuse_fileManager,only:INPUT_PATH ! defines data directory
USE multiforce,only:forcefile ! name of forcing file
Expand Down Expand Up @@ -61,7 +108,7 @@ SUBROUTINE get_modtim(itim,ncid,ierr,message)
if(ierr/=0)then; message=trim(message)//trim(nf90_strerror(ierr)); return; endif

! put the time into the structure
timDat%dtime = aTime(1)
timDat%dtime = aTime(1) * time_units_to_days(timeUnits)

! compute the year, month, day, hour, minute, second
call caldatss(jdayRef+timDat%dtime,timDat%iy,timDat%im,timDat%id,timDat%ih,timDat%imin,timDat%dsec)
Expand Down