diff --git a/build/FUSE_SRC/driver/functn.f90 b/build/FUSE_SRC/driver/functn.f90 index 91c0f4f..fa5acf7 100644 --- a/build/FUSE_SRC/driver/functn.f90 +++ b/build/FUSE_SRC/driver/functn.f90 @@ -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 @@ -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 diff --git a/build/FUSE_SRC/driver/fuse_driver.f90 b/build/FUSE_SRC/driver/fuse_driver.f90 index 6be2133..796b6e9 100644 --- a/build/FUSE_SRC/driver/fuse_driver.f90 +++ b/build/FUSE_SRC/driver/fuse_driver.f90 @@ -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 @@ -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 @@ -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? @@ -119,7 +122,7 @@ 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 @@ -127,7 +130,6 @@ PROGRAM DISTRIBUTED_DRIVER 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 ! --------------------------------------------------------------------------------------- @@ -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 ! --------------------------------------------------------------------------------------- @@ -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) @@ -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 @@ -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 @@ -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)) @@ -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)) diff --git a/build/FUSE_SRC/driver/fuse_metric.f90 b/build/FUSE_SRC/driver/fuse_metric.f90 index c22a43b..fe1fdd1 100644 --- a/build/FUSE_SRC/driver/fuse_metric.f90 +++ b/build/FUSE_SRC/driver/fuse_metric.f90 @@ -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 @@ -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...' diff --git a/build/FUSE_SRC/dshare/globaldata.f90 b/build/FUSE_SRC/dshare/globaldata.f90 index be2049c..070d824 100644 --- a/build/FUSE_SRC/dshare/globaldata.f90 +++ b/build/FUSE_SRC/dshare/globaldata.f90 @@ -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 diff --git a/build/FUSE_SRC/netcdf/get_mbands.f90 b/build/FUSE_SRC/netcdf/get_mbands.f90 index f05a6ba..6a4e4c8 100644 --- a/build/FUSE_SRC/netcdf/get_mbands.f90 +++ b/build/FUSE_SRC/netcdf/get_mbands.f90 @@ -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 @@ -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 @@ -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 @@ -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)