Skip to content

Cherry-pick run_pre fix and related improvements from PR #10 #23

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

Merged
DarriEy merged 3 commits into from
Feb 4, 2026
Merged

Cherry-pick run_pre fix and related improvements from PR #10 #23

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
94293dd
fix run_pre option
Nov 29, 2025
a84ccec
fix uninitialized error code in get_mbands
Dec 1, 2025
e82a777
add random seed to obtain same results for different optimization trials
Feb 4, 2026
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
4 changes: 4 additions & 0 deletions build/FUSE_SRC/driver/functn.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@ FUNCTION FUNCTN(NOPT,A)
USE fuse_metric_module ! run model and compute the metric chosen as objective function
USE multiforce, only: ncid_forc ! NetCDF forcing file ID
USE fuse_fileManager,only:METRIC, TRANSFO ! metric and transformation requested in the filemanager
USE globaldata, only: nFUSE_eval ! # fuse evaluations

IMPLICIT NONE
! input
Expand All @@ -31,6 +32,9 @@ FUNCTION FUNCTN(NOPT,A)
REAL(MSP) :: FUNCTN ! objective function value

! ---------------------------------------------------------------------------------------

nFUSE_eval = nFUSE_eval + 1

! get SCE parameter set
ALLOCATE(SCE_PAR(NOPT), STAT=IERR); IF (IERR.NE.0) STOP ' problem allocating space '
SCE_PAR(1:NOPT) = A(1:NOPT) ! convert from MSP used in SCE to SP used in FUSE
Expand Down
94 changes: 37 additions & 57 deletions build/FUSE_SRC/driver/fuse_driver.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,8 @@ PROGRAM DISTRIBUTED_DRIVER
! data modules
USE model_defn,nstateFUSE=>nstate ! model definition structures
USE model_defnames ! defines the integer model options
USE globaldata, ONLY: isPrint ! flag for printing progress to screen
USE globaldata, only: nFUSE_eval ! number of fuse evaluations
USE multiforce, ONLY: forcefile,vname_aprecip ! model forcing structures
USE multiforce, ONLY: AFORCE, aValid ! time series of lumped forcing/response data
USE multiforce, ONLY: nspat1, nspat2 ! grid dimensions
Expand All @@ -40,7 +42,7 @@ PROGRAM DISTRIBUTED_DRIVER
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 multiforce, only: NUMPSET,name_psets ! number of parameter set and their names
USE multiforce, only: NUMPSET ! number of parameter sets

USE multiforce, only: ncid_forc ! NetCDF forcing file ID
USE multiforce, only: ncid_var ! NetCDF forcing variable ID
Expand Down Expand Up @@ -84,14 +86,15 @@ PROGRAM DISTRIBUTED_DRIVER
CHARACTER(LEN=256) :: DatString ! file manager
CHARACTER(LEN=256) :: dom_id ! ID of the domain
CHARACTER(LEN=10) :: fuse_mode=' ' ! fuse execution mode (run_def, run_best, run_pre, calib_sce)
CHARACTER(LEN=256) :: file_para_list ! txt file containing list of parameter sets
CHARACTER(LEN=256) :: file_param ! name of parameter file
CHARACTER(LEN=10) :: index_param ! index of desired parameter set

! ---------------------------------------------------------------------------------------
! SETUP MODELS FOR SIMULATION -- POPULATE DATA STRUCTURES
! ---------------------------------------------------------------------------------------
! fuse_file_manager
CHARACTER(LEN=1024) :: FFMFILE ! name of fuse_file_manager file
CHARACTER(LEN=1024) :: ELEV_BANDS_NC ! name of NetCDF file for elevation bands
CHARACTER(LEN=1024) :: FFMFILE ! name of fuse_file_manager file
CHARACTER(LEN=1024) :: ELEV_BANDS_NC ! name of NetCDF file for elevation bands
! get model forcing data
INTEGER(I4B) :: NTIM ! number of time steps - still needed ?
INTEGER(I4B) :: INFERN_START ! start of inference period - still needed?
Expand Down Expand Up @@ -119,15 +122,14 @@ PROGRAM DISTRIBUTED_DRIVER
! ---------------------------------------------------------------------------------------
INTEGER(I4B) :: ITIM ! loop thru time steps
INTEGER(I4B) :: IPAR ! loop thru model parameters
INTEGER(I4B) :: IPSET ! loop thru model parameter sets
INTEGER(I4B) :: IPSET ! index of desired model parameter set
TYPE(PARATT) :: PARAM_META ! parameter metadata (model parameters)
REAL(SP), DIMENSION(:), ALLOCATABLE :: BL ! vector of lower parameter bounds
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
REAL(KIND=4),DIMENSION(:), ALLOCATABLE :: URAND ! vector of quasi-random numbers U[0,1]
REAL(SP) :: METRIC_VAL ! error from the simulation

! ---------------------------------------------------------------------------------------
! SCE VARIABLES
! ---------------------------------------------------------------------------------------
Expand Down Expand Up @@ -177,22 +179,28 @@ PROGRAM DISTRIBUTED_DRIVER
CALL GETARG(1,DatString) ! string defining forcinginfo file
CALL GETARG(2,dom_id) ! ID of the domain
CALL GETARG(3,fuse_mode) ! fuse execution mode (run_def, run_best, calib_sce)
IF(TRIM(fuse_mode).EQ.'run_pre') CALL GETARG(4,file_para_list) ! fuse execution mode txt file containing list of parameter sets
IF(TRIM(fuse_mode).EQ.'run_pre')then
CALL GETARG(4,file_param) ! name of parameter file
CALL GETARG(5,index_param) ! index of desired parameter set
IF(LEN_TRIM(index_param).EQ.0) IPSET = 1
IF(LEN_TRIM(index_param).GT.0) read(index_param,*) IPSET
ENDIF

! check command-line arguments
IF (LEN_TRIM(DatString).EQ.0) STOP '1st command-line argument is missing (fileManager)'
IF (LEN_TRIM(dom_id).EQ.0) STOP '2nd command-line argument is missing (dom_id)'
IF (LEN_TRIM(fuse_mode).EQ.0) STOP '3rd command-line argument is missing (fuse_mode)'
IF(TRIM(fuse_mode).EQ.'run_pre')THEN
IF(LEN_TRIM(file_para_list).EQ.0) STOP '4th command-line argument is missing (file_para_list) and is required in mode run_pre'
IF(LEN_TRIM(file_param).EQ.0) STOP '4th command-line argument is missing (file_param) and is required in mode run_pre'
ENDIF

! print command-line arguments
print*, '1st command-line argument (fileManager) = ', trim(DatString)
print*, '2nd command-line argument (dom_id) = ', trim(dom_id)
print*, '3rd command-line argument (fuse_mode) = ', fuse_mode
print*, '1st command-line argument (fileManager) = ', trim(DatString)
print*, '2nd command-line argument (dom_id) = ', trim(dom_id)
print*, '3rd command-line argument (fuse_mode) = ', fuse_mode
IF(TRIM(fuse_mode).EQ.'run_pre')THEN
print*, '4th command-line argument (file_para_list) = ', file_para_list
print*, '4th command-line argument (file_param) = ', file_param
print*, '5th command-line argument (index_param) = ', IPSET
ENDIF

! ---------------------------------------------------------------------------------------
Expand Down Expand Up @@ -256,6 +264,7 @@ PROGRAM DISTRIBUTED_DRIVER

! get elevation band info, in particular N_BANDS
CALL GET_MBANDS_INFO(ELEV_BANDS_NC,err,message) ! read band data from NetCDF file
if(err/=0)then; write(*,*) trim(message); stop; endif

! allocate space for elevation bands
allocate(MBANDS_VAR_4d(nspat1,nspat2,N_BANDS,numtim_sub+1),stat=err)
Expand Down Expand Up @@ -300,44 +309,15 @@ PROGRAM DISTRIBUTED_DRIVER
#endif

NUMPSET=1 ! only the default parameter set is run
ALLOCATE(name_psets(NUMPSET))
name_psets(1)='default_param_set'

ELSE IF(fuse_mode == 'run_pre')THEN ! run FUSE with pre-defined parameter values

! read file_para_list twice:
! 1st pass: determine number of parameter set and allocate name_psets accordingly
! 2st pass: save the names of parameter sets in name_psets

do file_pass = 1, 2

NUMPSET=0 ! intialize counter

OPEN(21,FILE=TRIM(file_para_list))
DO ! loop through parameter files

READ(21,*,IOSTAT=ERR) dummy_string
IF (ERR.NE.0) EXIT
NUMPSET=NUMPSET+1 ! increment counter

if (file_pass.eq.2) THEN
name_psets(NUMPSET) = dummy_string ! save file names
ENDIF

END DO ! looping through parameter files

CLOSE(21)

if(file_pass.eq.1) THEN
print *, 'NUMPSET=', NUMPSET, 'based on the number of lines in ', TRIM(file_para_list)
ALLOCATE(name_psets(NUMPSET))
END IF
end do

! files to which model run and parameter set will be saved
FNAME_NETCDF_RUNS = TRIM(OUTPUT_PATH)//TRIM(dom_id)//'_'//TRIM(FMODEL_ID)//'_runs_pre.nc'
FNAME_NETCDF_PARA = TRIM(OUTPUT_PATH)//TRIM(dom_id)//'_'//TRIM(FMODEL_ID)//'_para_pre_out.nc'

NUMPSET=1 ! only the one "desired" parameter set is run

ELSE IF(fuse_mode == 'calib_sce')THEN ! calibrate FUSE using SCE

! files to which model run and parameter set will be saved
Expand Down Expand Up @@ -415,22 +395,15 @@ PROGRAM DISTRIBUTED_DRIVER

OUTPUT_FLAG=.TRUE.

do IPSET = 1, NUMPSET

FNAME_NETCDF_PARA_PRE=TRIM(OUTPUT_PATH)//name_psets(IPSET)
PRINT *, 'Loading parameter set ',IPSET,':'

! load specific parameter set
! 2nd argument is 1 because first (and only) parameter set should be loaded
CALL GET_PRE_PARAM(FNAME_NETCDF_PARA_PRE,1,ONEMOD,NUMPAR,APAR)
FNAME_NETCDF_PARA_PRE=TRIM(OUTPUT_PATH)//file_param
PRINT *, 'Loading parameter set ',IPSET,':'

print *, 'Running FUSE with pre-defined parameter set'
CALL FUSE_METRIC(APAR,GRID_FLAG,NCID_FORC,METRIC_VAL,OUTPUT_FLAG,IPSET)
print *, 'Done running FUSE with pre-defined parameter set'
! load specific parameter set
CALL GET_PRE_PARAM(FNAME_NETCDF_PARA_PRE,IPSET,ONEMOD,NUMPAR,APAR)

end do

DEALLOCATE(name_psets)
print *, 'Running FUSE with pre-defined parameter set'
CALL FUSE_METRIC(APAR,GRID_FLAG,NCID_FORC,METRIC_VAL,OUTPUT_FLAG,1) ! last argument IPSET=1
print *, 'Done running FUSE with pre-defined parameter set'

ELSE IF(fuse_mode == 'calib_sce')THEN ! calibrate FUSE using SCE

Expand All @@ -439,6 +412,10 @@ PROGRAM DISTRIBUTED_DRIVER

FNAME_ASCII = TRIM(OUTPUT_PATH)//TRIM(dom_id)//'_'//TRIM(FMODEL_ID)//'_sce_output.txt'

! printing
isPrint = .false. ! turn off printing to screen
nFUSE_eval = 0 ! number of fuse evaluations

! convert from SP used in FUSE to MSP used in SCE
ALLOCATE(APAR_MSP(NUMPAR),BL_MSP(NUMPAR),BU_MSP(NUMPAR),URAND_MSP(NUMPAR))

Expand All @@ -448,6 +425,9 @@ PROGRAM DISTRIBUTED_DRIVER
BU_MSP=BU
URAND_MSP=URAND

! set random seed
ISEED = 1

! open up ASCII output file
print *, 'Creating SCE output file:', trim(FNAME_ASCII)
ISCE = 96; OPEN(ISCE,FILE=TRIM(FNAME_ASCII))
Expand Down
4 changes: 4 additions & 0 deletions build/FUSE_SRC/driver/fuse_metric.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,8 @@ SUBROUTINE FUSE_METRIC(XPAR,GRID_FLAG,NCID_FORC,METRIC_VAL,OUTPUT_FLAG,IPSET,MPA
! data modules
USE model_defn, ONLY:NSTATE,SMODL ! number of state variables
USE model_defnames ! integer model definitions
USE globaldata, ONLY: isPrint ! flag for printing progress to screen
USE globaldata, only: nFUSE_eval ! number of fuse evaluations
USE multiparam, ONLY: LPARAM,NUMPAR,MPARAM ! list of model parameters
USE multiforce, ONLY: MFORCE,AFORCE,DELTIM,ISTART ! model forcing data
USE multiforce, ONLY: numtim_in, itim_in ! length of input time series and associated index
Expand Down Expand Up @@ -357,6 +359,8 @@ SUBROUTINE FUSE_METRIC(XPAR,GRID_FLAG,NCID_FORC,METRIC_VAL,OUTPUT_FLAG,IPSET,MPA
CALL MEAN_STATS()
METRIC_VAL = MSTATS%METRIC_VAL

write(*,'(i6,1x,a12,1x,f12.6)') nFUSE_eval, "METRIC_VAL =", METRIC_VAL

ENDIF

PRINT *, 'Writing parameter values...'
Expand Down
2 changes: 1 addition & 1 deletion build/FUSE_SRC/dshare/globaldata.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,6 @@ MODULE globaldata
integer(i4b), parameter :: iPERR=7 ! not a snow parameter, but used here

! number of fuse evaluations
integer(i4b), save :: nFUSE_eval
integer(i4b), save :: nFUSE_eval

end MODULE globaldata
54 changes: 27 additions & 27 deletions build/FUSE_SRC/netcdf/get_mbands.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -160,7 +160,7 @@ SUBROUTINE GET_MBANDS_INFO(ELEV_BANDS_NC,err,message)
character(*), intent(out) :: message
! internal
integer(i4b),parameter::lenPath=1024 ! DK/2008/10/21: allows longer file paths
INTEGER(I4B),DIMENSION(2) :: IERR ! error codes
INTEGER(I4B) :: IERR ! error code
INTEGER(I4B) :: IUNIT ! input file unit
CHARACTER(LEN=lenPath) :: CFILE ! name of control file
CHARACTER(LEN=lenPath) :: BFILE ! name of band file
Expand All @@ -180,21 +180,21 @@ SUBROUTINE GET_MBANDS_INFO(ELEV_BANDS_NC,err,message)
! internal: NetCDF read
integer(i4b) :: ivarid_af,ivarid_me ! NetCDF variable ID for area_frac and mean_area
integer(i4b),parameter :: ndims=3 ! number of dimensions for frac_area
integer(i4b) :: dimid_eb ! ID elevation bands
integer(i4b) :: dimid_eb ! ID elevation bands
integer(i4b) :: iDimID ! dimension ID
integer(i4b) :: dimLen ! dimension length

! ---------------------------------------------------------------------------------------

! read in NetCDF file defining the elevation bands
err=0
err=0; ierr=0
CFILE = TRIM(INPUT_PATH)//ELEV_BANDS_NC ! control file info shared in MODULE directory
print *, 'Loading elevation bands from:',TRIM(CFILE)

INQUIRE(FILE=CFILE,EXIST=LEXIST) ! check that control file exists
IF (.NOT.LEXIST) THEN
print *, 'f-GET_MBANDS_GRID/NetCDF file ',TRIM(CFILE),' for elevation bands does not exist '
STOP
print *, 'f-GET_MBANDS_GRID/NetCDF file ',TRIM(CFILE),' for elevation bands does not exist '
STOP
ENDIF

!open netcdf file
Expand All @@ -220,14 +220,14 @@ SUBROUTINE GET_MBANDS_INFO(ELEV_BANDS_NC,err,message)
if(err/=0)then; message=trim(message)//trim(nf90_strerror(err)); return; endif

! allocate 1 data stucture
ALLOCATE(MBANDS(N_BANDS),STAT=IERR(1))
ALLOCATE(MBANDS(N_BANDS),STAT=IERR)

! allocate data structures
ALLOCATE(Z_FORCING_grid(nspat1,nspat2),MBANDS_INFO_3d(nspat1,nspat2,n_bands),&
AF_TEMP(nspat1,nspat2,n_bands),ME_TEMP(nspat1,nspat2,n_bands),&
elev_mask(nspat1,nspat2),STAT=IERR(1))
AF_TEMP(nspat1,nspat2,n_bands),ME_TEMP(nspat1,nspat2,n_bands),&
elev_mask(nspat1,nspat2),STAT=IERR)

IF (ANY(IERR.NE.0)) THEN
IF (IERR.NE.0) THEN
message="f-GET_MBANDS/problem allocating elevation band data structures"
err=100; return
ENDIF
Expand All @@ -242,26 +242,26 @@ SUBROUTINE GET_MBANDS_INFO(ELEV_BANDS_NC,err,message)

! populate MBANDS_INFO_3d, Z_FORCING_grid and elev_mask
DO iSpat2=1,nSpat2
DO iSpat1=1,nSpat1
DO iSpat1=1,nSpat1

MBANDS_INFO_3d(iSpat1,iSpat2,:)%Z_MID = me_TEMP(iSpat1,iSpat2,:)
MBANDS_INFO_3d(iSpat1,iSpat2,:)%AF = af_TEMP(iSpat1,iSpat2,:)
Z_FORCING_grid(iSpat1,iSpat2) = sum(me_TEMP(iSpat1,iSpat2,:)*af_TEMP(iSpat1,iSpat2,:)) ! estimate mean elevation of forcing using weighted mean of EB elevation
elev_mask(iSpat1,iSpat2) = me_TEMP(iSpat1,iSpat2,1) .EQ. NA_VALUE_SP ! if mean elevation first band is NA_VALUE, mask this grid cell
MBANDS_INFO_3d(iSpat1,iSpat2,:)%Z_MID = me_TEMP(iSpat1,iSpat2,:)
MBANDS_INFO_3d(iSpat1,iSpat2,:)%AF = af_TEMP(iSpat1,iSpat2,:)
Z_FORCING_grid(iSpat1,iSpat2) = sum(me_TEMP(iSpat1,iSpat2,:)*af_TEMP(iSpat1,iSpat2,:)) ! estimate mean elevation of forcing using weighted mean of EB elevation
elev_mask(iSpat1,iSpat2) = me_TEMP(iSpat1,iSpat2,1) .EQ. NA_VALUE_SP ! if mean elevation first band is NA_VALUE, mask this grid cell

if(.NOT.elev_mask(iSpat1,iSpat2)) THEN ! only check area fraction sum to 1 if not NA_VALUE

if (abs(sum(MBANDS_INFO_3d(iSpat1,iSpat2,:)%AF)-1).GT.1E-2) then ! check that area fraction sum to 1

print *, "The area fraction of all the elevation bands do not add up to 1"
!print *, 'Difference with 1 = ', abs(sum(MBANDS_INFO_3d(iSpat1,iSpat2,:)%AF)-1)
print *, 'AF', MBANDS_INFO_3d(iSpat1,iSpat2,:)%AF
stop

end if
end if

if(.NOT.elev_mask(iSpat1,iSpat2)) THEN ! only check area fraction sum to 1 if not NA_VALUE

if (abs(sum(MBANDS_INFO_3d(iSpat1,iSpat2,:)%AF)-1).GT.1E-2) then ! check that area fraction sum to 1

print *, "The area fraction of all the elevation bands do not add up to 1"
!print *, 'Difference with 1 = ', abs(sum(MBANDS_INFO_3d(iSpat1,iSpat2,:)%AF)-1)
print *, 'AF', MBANDS_INFO_3d(iSpat1,iSpat2,:)%AF
stop

end if
end if

END DO
END DO
END DO

err = nf90_close(ncid_eb)
Expand Down