diff --git a/build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 b/build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 index 32b4cee..af418eb 100644 --- a/build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 +++ b/build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 @@ -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 @@ -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 @@ -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 diff --git a/build/FUSE_SRC/FUSE_ENGINE/force_info.f90 b/build/FUSE_SRC/FUSE_ENGINE/force_info.f90 index 02fcbcb..5356f80 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/force_info.f90 +++ b/build/FUSE_SRC/FUSE_ENGINE/force_info.f90 @@ -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 @@ -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),'/') diff --git a/build/FUSE_SRC/FUSE_ENGINE/get_time_indices.f90 b/build/FUSE_SRC/FUSE_ENGINE/get_time_indices.f90 index e71333e..82e22d1 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/get_time_indices.f90 +++ b/build/FUSE_SRC/FUSE_ENGINE/get_time_indices.f90 @@ -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 @@ -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 diff --git a/build/FUSE_SRC/FUSE_ENGINE/mean_stats.f90 b/build/FUSE_SRC/FUSE_ENGINE/mean_stats.f90 index e80b641..e509573 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/mean_stats.f90 +++ b/build/FUSE_SRC/FUSE_ENGINE/mean_stats.f90 @@ -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 diff --git a/build/FUSE_SRC/FUSE_ENGINE/metaoutput.f90 b/build/FUSE_SRC/FUSE_ENGINE/metaoutput.f90 index 66765a1..e6e354a 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/metaoutput.f90 +++ b/build/FUSE_SRC/FUSE_ENGINE/metaoutput.f90 @@ -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 diff --git a/build/FUSE_SRC/FUSE_ENGINE/multiforce.f90 b/build/FUSE_SRC/FUSE_ENGINE/multiforce.f90 index 90d6ec6..6c648fd 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/multiforce.f90 +++ b/build/FUSE_SRC/FUSE_ENGINE/multiforce.f90 @@ -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 @@ -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 @@ -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 @@ -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 diff --git a/build/FUSE_SRC/FUSE_ENGINE/multiroute.f90 b/build/FUSE_SRC/FUSE_ENGINE/multiroute.f90 index f9d046b..c09fcc6 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/multiroute.f90 +++ b/build/FUSE_SRC/FUSE_ENGINE/multiroute.f90 @@ -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 diff --git a/build/FUSE_SRC/FUSE_ENGINE/varextract.f90 b/build/FUSE_SRC/FUSE_ENGINE/varextract.f90 index f73f766..8f2f346 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/varextract.f90 +++ b/build/FUSE_SRC/FUSE_ENGINE/varextract.f90 @@ -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 @@ -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 @@ -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 @@ -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 diff --git a/build/FUSE_SRC/FUSE_HOOK/fuse_fileManager.f90 b/build/FUSE_SRC/FUSE_HOOK/fuse_fileManager.f90 index 2352ec5..a5f741d 100644 --- a/build/FUSE_SRC/FUSE_HOOK/fuse_fileManager.f90 +++ b/build/FUSE_SRC/FUSE_HOOK/fuse_fileManager.f90 @@ -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 @@ -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 diff --git a/build/FUSE_SRC/FUSE_NETCDF/get_gforce.f90 b/build/FUSE_SRC/FUSE_NETCDF/get_gforce.f90 index f990a04..1f7b782 100644 --- a/build/FUSE_SRC/FUSE_NETCDF/get_gforce.f90 +++ b/build/FUSE_SRC/FUSE_NETCDF/get_gforce.f90 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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. @@ -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 diff --git a/build/FUSE_SRC/FUSE_TIME/time_io.f90 b/build/FUSE_SRC/FUSE_TIME/time_io.f90 index 9345aac..174e88b 100644 --- a/build/FUSE_SRC/FUSE_TIME/time_io.f90 +++ b/build/FUSE_SRC/FUSE_TIME/time_io.f90 @@ -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) ! --------------------------------------------------------------------------------------- @@ -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 @@ -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)