diff --git a/build/FUSE_SRC/FUSE_DMSL/functn.f90 b/build/FUSE_SRC/FUSE_DMSL/functn.f90 index b7cfa89..a0e78fe 100644 --- a/build/FUSE_SRC/FUSE_DMSL/functn.f90 +++ b/build/FUSE_SRC/FUSE_DMSL/functn.f90 @@ -9,13 +9,13 @@ FUNCTION FUNCTN(NOPT,A) ! Wrapper for SCE (used to compute the objective function) ! --------------------------------------------------------------------------------------- USE nrtype ! variable types, etc. -USE FUSE_RMSE_MODULE ! run model and compute the root mean squared error -USE multiforce, only: ncid_forc ! NetCDF forcing file ID +USE FUSE_RMSE_MODULE ! run model and compute the root mean squared error +USE multiforce, only: ncid_forc ! NetCDF forcing file ID IMPLICIT NONE ! input INTEGER(I4B) :: NOPT ! number of parameters -REAL(MSP), DIMENSION(100), INTENT(IN) :: A ! model parameter set - can be bumped up to 100 elements +REAL(MSP), DIMENSION(100), INTENT(IN) :: A ! model parameter set - can be bumped up to 100 elements ! internal REAL(SP), DIMENSION(:), ALLOCATABLE :: SCE_PAR ! sce parameter set @@ -39,7 +39,6 @@ FUNCTION FUNCTN(NOPT,A) ! deallocate parameter set DEALLOCATE(SCE_PAR, STAT=IERR); IF (IERR.NE.0) STOP ' problem deallocating space ' -print *, 'RMSE =', RMSE ! save objective function value FUNCTN = RMSE diff --git a/build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 b/build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 index 4d2c7d0..582e5d2 100644 --- a/build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 +++ b/build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 @@ -23,6 +23,7 @@ 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 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 @@ -43,10 +44,11 @@ PROGRAM DISTRIBUTED_DRIVER USE multiforce, only: ncid_forc ! NetCDF forcing file ID USE multiforce, only: ncid_var ! NetCDF forcing variable ID -USE multistate, only: ncid_out ! NetCDF output file ID +USE globaldata, only: ncid_out ! NetCDF output file ID USE multibands ! basin band stuctures -USE multiparam, ONLY: LPARAM, PARATT, NUMPAR ! parameter metadata structures +USE data_types, ONLY: PARATT ! data type for metadata +USE multiparam, ONLY: LPARAM, NUMPAR ! parameter metadata structures USE multistate, only: gState ! gridded state variables USE multistate, only: gState_3d ! gridded state variables with a time dimension USE multiroute, ONLY: AROUTE ! model routing structures @@ -345,9 +347,9 @@ PROGRAM DISTRIBUTED_DRIVER ! assign algorithmic control parameters for SCE ! convert characters to interger/MSP - READ (MAXN_STR,*) MAXN ! maximum number of trials before optimization is terminated + 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 + READ (PCENTO_STR,*) PCENTO ! the percentage PRINT *, 'SCE parameters read from file manager:' PRINT *, 'Maximum number of trials before SCE optimization is stopped (MAXN) = ', MAXN_STR @@ -438,11 +440,13 @@ PROGRAM DISTRIBUTED_DRIVER FNAME_ASCII = TRIM(OUTPUT_PATH)//TRIM(dom_id)//'_'//TRIM(FMODEL_ID)//'_sce_output.txt' + ! turn off printing to screen + isPrint = .false. + ! convert from SP used in FUSE to MSP used in SCE ALLOCATE(APAR_MSP(NUMPAR),BL_MSP(NUMPAR),BU_MSP(NUMPAR),URAND_MSP(NUMPAR)) APAR_MSP=APAR - PRINT *, 'BL=',BL BL_MSP=BL BU_MSP=BU URAND_MSP=URAND diff --git a/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 b/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 index b64355c..7996ebb 100644 --- a/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 +++ b/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 @@ -22,6 +22,9 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG ! 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: fracstate0 ! fraction of initial state (used for initialization) + USE globaldata, ONLY: NA_VALUE, NA_VALUE_SP ! NA_VALUE for the forcing 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 @@ -34,9 +37,7 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG USE multiforce, ONLY:nspat1,nspat2 ! spatial dimensions USE multiforce, ONLY:ncid_var ! NetCDF ID for forcing variables USE multiforce, ONLY:gForce,gForce_3d ! gridded forcing data - USE multistate, ONLY:fracstate0,TSTATE,MSTATE,FSTATE,& ! model states - HSTATE ! model states (continued) - USE multiforce, ONLY:NA_VALUE, NA_VALUE_SP ! NA_VALUE for the forcing + USE multistate, ONLY:TSTATE,MSTATE,FSTATE,HSTATE ! model state variables USE multistate, ONLY:gState,gState_3d ! gridded state variables USE multiroute, ONLY:MROUTE,AROUTE,AROUTE_3d ! routed runoff USE multistats, ONLY:MSTATS,PCOUNT,MOD_IX ! access model statistics; counter for param set @@ -52,6 +53,11 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG USE str_2_xtry_module ! provide access to the routine str_2_xtry USE xtry_2_str_module ! provide access to the routine xtry_2_str + ! differentiable model + use data_types, only: parent ! fuse parent data type + use get_parent_module, only: get_parent ! populate the parent data structure + use implicit_solve_module, only:implicit_solve ! simple implicit solve for differnetiable ODE + ! interface blocks USE interfaceb, ONLY:ode_int,fuse_solve ! provide access to FUSE_SOLVE through ODE_INT @@ -93,6 +99,10 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG CHARACTER(LEN=CLEN) :: CMESSAGE ! error message of downwind routine INTEGER(I4B),PARAMETER::UNT=6 !1701 ! 6 + ! differentiable model + type(parent) :: fuseStruct ! parent fuse data structure + + ! --------------------------------------------------------------------------------------- ! allocate state vectors ALLOCATE(STATE0(NSTATE),STATE1(NSTATE),STAT=IERR) @@ -107,13 +117,16 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG ! add parameter set to the data structure CALL PUT_PARSET(XPAR) - PRINT *, 'Parameter set added to data structure:' - PRINT *, XPAR + if(isPrint) PRINT *, 'Parameter set added to data structure:' + if(isPrint) PRINT *, XPAR ! compute derived model parameters (bucket sizes, etc.) CALL PAR_DERIVE(ERR,MESSAGE) IF (ERR.NE.0) WRITE(*,*) TRIM(MESSAGE); IF (ERR.GT.0) STOP + if(isPrint) PRINT *, 'Writing parameter values...' + CALL PUT_PARAMS(PCOUNT) + ! initialize model states over the 2D gridded domain (1x1 domain in catchment mode) DO iSpat2=1,nSpat2 DO iSpat1=1,nSpat1 @@ -121,10 +134,10 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG gState_3d(iSpat1,iSpat2,1) = FSTATE ! put the state into first time step of 3D structure END DO END DO - PRINT *, 'Model states initialized over the 2D gridded domain' + if(isPrint) PRINT *, 'Model states initialized over the 2D gridded domain' ! initialize elevations bands if snow module is on - PRINT *, 'N_BANDS =', N_BANDS + if(isPrint) PRINT *, 'N_BANDS =', N_BANDS IF (SMODL%iSNOWM.EQ.iopt_temp_index) THEN DO iSpat2=1,nSpat2 @@ -137,7 +150,7 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG END DO END DO END DO - PRINT *, 'Snow states initiatlized over the 2D gridded domain ' + if(isPrint) PRINT *, 'Snow states initiatlized over the 2D gridded domain ' ENDIF ! allocate 3d data structure for fluxes @@ -177,10 +190,10 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG numtim_sub_cur=MIN(numtim_sub,numtim_sim-itim_sim+1) ! load forcing for desired period into gForce_3d - PRINT *, 'New subperiod: loading forcing for ',numtim_sub_cur,' time steps' + if(isPrint) PRINT *, 'New subperiod: loading forcing for ',numtim_sub_cur,' time steps' CALL get_gforce_3d(itim_in,numtim_sub_cur,ncid_forc,err,message) IF(err/=0)THEN; WRITE(*,*) 'Error while extracting 3d forcing'; STOP; ENDIF - PRINT *, 'Forcing loaded. Running FUSE...' + if(isPrint) PRINT *, 'Forcing loaded. Running FUSE...' ENDIF @@ -245,9 +258,38 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG RETURN END SELECT - ! temporally integrate the ordinary differential equations - CALL ODE_INT(FUSE_SOLVE,STATE0,STATE1,DT_SUB,DT_FULL,IERR,MESSAGE) - IF (IERR.NE.0) THEN; PRINT *, TRIM(MESSAGE); PAUSE; ENDIF + ! ----- start of soil physics code ------------------------------------------------------------ + + ! temporally integrate the ordinary differential equations + select case(diff_mode) + + ! original code + case(original) + CALL ODE_INT(FUSE_SOLVE,STATE0,STATE1,DT_SUB,DT_FULL,IERR,MESSAGE) + IF (IERR.NE.0) THEN; PRINT *, TRIM(MESSAGE); STOP 1; ENDIF + + !print*, state1 + !if(ITIM_IN > sim_beg+100) stop + + ! differentiable code + case(differentiable) + + ! populate parent fuse structure + call get_parent(fuseStruct) + + ! solve differentiable ODEs + call implicit_solve(fuseStruct, state0, state1, nState) + !print*, state1 + !if(ITIM_IN > sim_beg+100) stop + + ! save fluxes + W_FLUX = fuseStruct%flux + + ! check options + case default; print*, "Cannot identify diff_mode"; stop 1 + end select + + ! ----- end of soil physics code -------------------------------------------------------------- ! perform overland flow routing CALL Q_OVERLAND() @@ -307,15 +349,15 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG ! if end of subperiod: write to output file and save states IF(itim_sub.EQ.numtim_sub_cur)THEN - PRINT *, 'End of subperiod reached:' + if(isPrint) PRINT *, 'End of subperiod reached:' ! write model output IF (OUTPUT_FLAG) THEN - PRINT *, 'Write output for ',numtim_sub_cur,' time steps starting at indice', itim_sim-numtim_sub_cur+1 + if(isPrint) PRINT *, 'Write output for ',numtim_sub_cur,' time steps starting at indice', itim_sim-numtim_sub_cur+1 CALL PUT_GOUTPUT_3D(itim_sim-numtim_sub_cur+1,itim_in-numtim_sub_cur+1,numtim_sub_cur,IPSET) - PRINT *, 'Done writing output' + if(isPrint) PRINT *, 'Done writing output' ELSE - PRINT *, 'OUTPUT_FLAG is set on FALSE, no output written' + if(isPrint) PRINT *, 'OUTPUT_FLAG is set on FALSE, no output written' END IF ! TODO: set gState_3d and MBANDS_VAR_4d to NA @@ -344,20 +386,20 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG ! get timing information CALL CPU_TIME(T2) - WRITE(*,*) "TIME ELAPSED = ", t2-t1 + if(isPrint) WRITE(*,*) "TIME ELAPSED = ", t2-t1 ! calculate mean summary statistics IF(.NOT.GRID_FLAG)THEN - PRINT *, 'Calculating performance metrics...' + if(isPrint) PRINT *, 'Calculating performance metrics...' CALL MEAN_STATS() RMSE = MSTATS%RAW_RMSE + print*, "NSE = ", MSTATS%NASH_SUTT + ENDIF - PRINT *, 'Writing parameter values...' - CALL PUT_PARAMS(PCOUNT) - PRINT *, 'Writing model statistics...' + if(isPrint) PRINT *, 'Writing model statistics...' CALL PUT_SSTATS(PCOUNT) ! deallocate vectors diff --git a/build/FUSE_SRC/FUSE_ENGINE/assign_par.f90 b/build/FUSE_SRC/FUSE_ENGINE/assign_par.f90 index 5558eae..3bf82e9 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/assign_par.f90 +++ b/build/FUSE_SRC/FUSE_ENGINE/assign_par.f90 @@ -16,7 +16,8 @@ SUBROUTINE ASSIGN_PAR() USE nrtype ! variable types, etc. USE model_defn ! model definition structure USE model_defnames -USE multiparam, ONLY : lparam, paratt, numpar ! model parameter structures +USE data_types, ONLY : paratt ! data type for metadata +USE multiparam, ONLY : lparam, numpar ! model parameter structures USE getpar_str_module ! access to SUBROUTINE get_par_str IMPLICIT NONE INTEGER(I4B) :: MPAR ! counter for number of parameters diff --git a/build/FUSE_SRC/FUSE_ENGINE/fuse_solve.f90 b/build/FUSE_SRC/FUSE_ENGINE/fuse_solve.f90 index dd4ec5b..5ae59e3 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/fuse_solve.f90 +++ b/build/FUSE_SRC/FUSE_ENGINE/fuse_solve.f90 @@ -10,8 +10,8 @@ SUBROUTINE FUSE_SOLVE(CALCDSDT,IE_SOLVE,SI_SOLVE,B_IMPOSE,AVG_FLUX,ADD_FLUX,NEWS ! (6) add fluxes from accepted sub-steps to the total timestep flux ! (7) estimate state at end of a full step, based on sum of fluxes USE nrtype ! variable definitions, etc. -USE multi_flux, ONLY: M_FLUX,FLUX_0,FLUX_1,W_FLUX,& ! model fluxes - CURRENT_DT ! model fluxes (continued) +USE globaldata, ONLY: CURRENT_DT +USE multi_flux, ONLY: M_FLUX,FLUX_0,FLUX_1,W_FLUX ! model fluxes USE multistate, ONLY: FSTATE,MSTATE,BSTATE,ESTATE,& ! model states DY_DT,DYDT_0,DYDT_1,HSTATE ! model states (continued) USE fminln, ONLY: fmin_x0p,fmin_dtp,fmin_dt2p,fmin_dseep ! variables used for residual vector in IE diff --git a/build/FUSE_SRC/FUSE_ENGINE/get_mbands.f90 b/build/FUSE_SRC/FUSE_ENGINE/get_mbands.f90 index f05a6ba..efd006c 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/get_mbands.f90 +++ b/build/FUSE_SRC/FUSE_ENGINE/get_mbands.f90 @@ -147,11 +147,12 @@ SUBROUTINE GET_MBANDS_INFO(ELEV_BANDS_NC,err,message) use nrtype,only:I4B,LGT,SP use utilities_dmsl_kit_FUSE,only:getSpareUnit,stripTrailString USE fuse_fileManager,only:INPUT_PATH,SETNGS_PATH ! defines data directory -USE fuse_fileManager,only:MBANDS_NC ! defines elevation bands +USE fuse_fileManager,only:MBANDS_NC ! defines elevation bands USE multibands,only:N_BANDS,MBANDS,MBANDS_INFO_3d,Z_FORCING,& - Z_FORCING_grid,elev_mask ! model band structures -USE multiforce,only:nspat1,nspat2,startSpat2,NA_VALUE_SP ! dimension lengths, na_value + Z_FORCING_grid,elev_mask ! model band structures +USE multiforce,only:nspat1,nspat2,startSpat2 ! dimension lengths +USE globaldata,only:NA_VALUE_SP IMPLICIT NONE ! dummies diff --git a/build/FUSE_SRC/FUSE_ENGINE/getnumerix.f90 b/build/FUSE_SRC/FUSE_ENGINE/getnumerix.f90 index 9a4d3c5..43b15e7 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/getnumerix.f90 +++ b/build/FUSE_SRC/FUSE_ENGINE/getnumerix.f90 @@ -18,13 +18,15 @@ SUBROUTINE GETNUMERIX(err,message) USE model_numerix,only:SOLUTION_METHOD,& ! defines numerix decisions TEMPORAL_ERROR_CONTROL,INITIAL_NEWTON,JAC_RECOMPUTE,CHECK_OVERSHOOT,SMALL_ENDSTEP,& ERR_TRUNC_ABS,ERR_TRUNC_REL,ERR_ITER_FUNC,ERR_ITER_DX,THRESH_FRZE,FRACSTATE_MIN,& - SAFETY,RMIN,RMAX,NITER_TOTAL,MIN_TSTEP,MAX_TSTEP + SAFETY,RMIN,RMAX,NITER_TOTAL,MIN_TSTEP,MAX_TSTEP,diff_mode +USE model_numerix,only:original,differentiable ! named variables for diff_mode IMPLICIT NONE ! dummies integer(I4B),intent(out) :: err character(*),intent(out) :: message ! locals INTEGER(I4B) :: IUNIT ! file unit +integer(i4b) :: ios ! io status flag integer(i4b),parameter::lenPath=1024 !DK/2008/10/21: allows longer file paths CHARACTER(LEN=lenPath) :: CFILE ! name of constraints file LOGICAL(LGT) :: LEXIST ! .TRUE. if file exists @@ -65,6 +67,12 @@ SUBROUTINE GETNUMERIX(err,message) READ(IUNIT,*) NITER_TOTAL ! Total number of iterations used in the implicit scheme READ(IUNIT,*) MIN_TSTEP ! Minimum time step length (minutes) READ(IUNIT,*) MAX_TSTEP ! Maximum time step length (minutes) +! new option -- ensure backwards compatible +read(iunit,*, iostat=ios) diff_mode ! Mode for differentiable models (non-differentiable; differentiable) +if(ios/=0)then + diff_mode = original + print*, "WARNING: diff_mode is not specified; setting option to original. Continuing" +endif ! if problem reading CLOSE(IUNIT) MIN_TSTEP = MIN_TSTEP/(24._SP*60._SP) ! Convert from minutes to days MAX_TSTEP = MAX_TSTEP/(24._SP*60._SP) ! Convert from minutes to days diff --git a/build/FUSE_SRC/FUSE_ENGINE/getpar_str.f90 b/build/FUSE_SRC/FUSE_ENGINE/getpar_str.f90 index 4b51b3e..bf3fd77 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/getpar_str.f90 +++ b/build/FUSE_SRC/FUSE_ENGINE/getpar_str.f90 @@ -13,7 +13,8 @@ SUBROUTINE GETPAR_STR(PARNAME,METADAT) ! Inserts parameter metadata into data structures ! --------------------------------------------------------------------------------------- USE nrtype ! variable types, etc. -USE multiparam, ONLY : PARATT, PARMETA ! derived type for parameter metadata +USE data_types, ONLY : PARATT ! derived type for parameter metadata +USE multiparam, ONLY : PARMETA ! parameter metadata IMPLICIT NONE ! input CHARACTER(*), INTENT(IN) :: PARNAME ! parameter name diff --git a/build/FUSE_SRC/FUSE_ENGINE/getparmeta.f90 b/build/FUSE_SRC/FUSE_ENGINE/getparmeta.f90 index 8c6313d..774fd7b 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/getparmeta.f90 +++ b/build/FUSE_SRC/FUSE_ENGINE/getparmeta.f90 @@ -14,7 +14,7 @@ SUBROUTINE GETPARMETA(err,message) ! --------------------------------------------------------------------------------------- USE nrtype ! variable types, etc. USE fuse_fileManager,only:SETNGS_PATH,CONSTRAINTS ! defines data directory -USE multiparam, ONLY: PARATT ! parameter attribute structure +USE data_types, ONLY: PARATT ! parameter attribute structure USE putpar_str_module ! provide access to SUBROUTINE putpar_str USE par_insert_module ! provide access to SUBROUTINE par_insert IMPLICIT NONE diff --git a/build/FUSE_SRC/FUSE_ENGINE/mean_stats.f90 b/build/FUSE_SRC/FUSE_ENGINE/mean_stats.f90 index 106c3e3..28fde13 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/mean_stats.f90 +++ b/build/FUSE_SRC/FUSE_ENGINE/mean_stats.f90 @@ -19,6 +19,8 @@ SUBROUTINE MEAN_STATS() USE multiroute ! routed runoff USE multi_flux ! fluxes USE multistats ! summary statistics +USE globaldata, ONLY: isPrint ! flag for printing progress to screen +USE globaldata, only: NA_VALUE_SP ! missing value USE model_numerix ! model numerix parameters and data IMPLICIT NONE @@ -51,7 +53,7 @@ SUBROUTINE MEAN_STATS() ! --------------------------------------------------------------------------------------- ! define sample size NS = eval_end-eval_beg+1 -PRINT *, 'Number of time steps in evaluation period (EP) = ', NS +if(isPrint) PRINT *, 'Number of time steps in evaluation period (EP) = ', NS ! allocate space for observed and simulated runoff ALLOCATE(QOBS(NS),QOBS_MASK(NS),QSIM(NS),STAT=IERR) @@ -63,10 +65,10 @@ SUBROUTINE MEAN_STATS() QOBS = aValid(1,1,eval_beg-sim_beg+1:eval_end-sim_beg+1)%OBSQ ! check for missing QOBS values -QOBS_MASK = QOBS.ne.REAL(NA_VALUE, KIND(SP)) ! find the time steps for which QOBS is available -NUM_AVAIL = COUNT(QOBS_MASK) ! number of time steps for which QOBS is available +QOBS_MASK = QOBS.ne.NA_VALUE_SP ! find the time steps for which QOBS is available +NUM_AVAIL = COUNT(QOBS_MASK) ! number of time steps for which QOBS is available -PRINT *, 'Number of time steps with observed streamflow in EP = ', NUM_AVAIL +if(isPrint) PRINT *, 'Number of time steps with observed streamflow in EP = ', NUM_AVAIL IF (NUM_AVAIL.EQ.0) THEN @@ -80,11 +82,11 @@ SUBROUTINE MEAN_STATS() ALLOCATE(QOBS_AVAIL(NUM_AVAIL),QSIM_AVAIL(NUM_AVAIL),DOBS(NUM_AVAIL),DSIM(NUM_AVAIL),RAWD(NUM_AVAIL),LOGD(NUM_AVAIL),STAT=IERR) QOBS_AVAIL=PACK(QOBS,QOBS_MASK,QOBS_AVAIL) ! moves QOBS time steps indicated by QOBS_MASK to QOBS_AVAIL, - ! if no values is missing (i.e. NS = NUM_AVAIL) then QOBS_AVAIL - ! should be a copy of QOBS + ! if no values is missing (i.e. NS = NUM_AVAIL) then QOBS_AVAIL + ! should be a copy of QOBS QSIM_AVAIL=PACK(QSIM,QOBS_MASK,QSIM_AVAIL) ! moves QSIM time steps indicated by QOBS_MASK to QSIM_AVAIL - ! if no values is missing (i.e. NS = NUM_AVAIL) then QSIM_AVAIL - ! should be a copy of QSIM + ! if no values is missing (i.e. NS = NUM_AVAIL) then QSIM_AVAIL + ! should be a copy of QSIM ! compute mean XB_OBS = SUM(QOBS_AVAIL(:)) / REAL(NUM_AVAIL, KIND(SP)) @@ -130,8 +132,8 @@ SUBROUTINE MEAN_STATS() END IF -PRINT *, 'NSE = ', MSTATS%NASH_SUTT -PRINT *, 'RAW_RMSE = ', MSTATS%RAW_RMSE +if(isPrint) PRINT *, 'NSE = ', MSTATS%NASH_SUTT +if(isPrint) PRINT *, 'RAW_RMSE = ', MSTATS%RAW_RMSE ! --------------------------------------------------------------------------------------- ! (3§) COMPUTE STATISTICS ON NUMERICAL ACCURACY AND EFFICIENCY diff --git a/build/FUSE_SRC/FUSE_ENGINE/multi_flux.f90 b/build/FUSE_SRC/FUSE_ENGINE/multi_flux.f90 deleted file mode 100644 index b3c884f..0000000 --- a/build/FUSE_SRC/FUSE_ENGINE/multi_flux.f90 +++ /dev/null @@ -1,42 +0,0 @@ -MODULE multi_flux - USE nrtype - TYPE FLUXES - REAL(SP) :: EFF_PPT ! effective precipitation (mm day-1) - REAL(SP) :: SATAREA ! saturated area (-) - REAL(SP) :: QSURF ! surface runoff (mm day-1) - REAL(SP) :: EVAP_1A ! evaporation from soil excess zone (mm day-1) - REAL(SP) :: EVAP_1B ! evaporation from soil recharge zone (mm day-1) - REAL(SP) :: EVAP_1 ! evaporation from upper soil layer (mm day-1) - REAL(SP) :: EVAP_2 ! evaporation from lower soil layer (mm day-1) - REAL(SP) :: RCHR2EXCS ! flow from recharge to excess (mm day-1) - REAL(SP) :: TENS2FREE_1 ! flow from tension storage to free storage (mm day-1) - REAL(SP) :: TENS2FREE_2 ! flow from tension storage to free storage (mm day-1) - REAL(SP) :: QINTF_1 ! interflow from free water (mm day-1) - REAL(SP) :: QPERC_12 ! percolation from upper to lower soil layers (mm day-1) - REAL(SP) :: QBASE_2 ! baseflow (mm day-1) - REAL(SP) :: QBASE_2A ! baseflow from primary linear resvr (mm day-1) - REAL(SP) :: QBASE_2B ! baseflow from secondary linear resvr (mm day-1) - REAL(SP) :: OFLOW_1 ! bucket overflow (mm day-1) - REAL(SP) :: OFLOW_2 ! bucket overflow (mm day-1) - REAL(SP) :: OFLOW_2A ! bucket overflow (mm day-1) - REAL(SP) :: OFLOW_2B ! bucket overflow (mm day-1) - REAL(SP) :: ERR_WATR_1 ! excessive extrapolation: total storage in layer1 (mm day-1) - REAL(SP) :: ERR_TENS_1 ! excessive extrapolation: tension storage in layer1 (mm day-1) - REAL(SP) :: ERR_FREE_1 ! excessive extrapolation: free storage in layer 1 (mm day-1) - REAL(SP) :: ERR_TENS_1A ! excessive extrapolation: storage in the recharge zone (mm day-1) - REAL(SP) :: ERR_TENS_1B ! excessive extrapolation: storage in the lower zone (mm day-1) - REAL(SP) :: ERR_WATR_2 ! excessive extrapolation: total storage in layer2 (mm day-1) - REAL(SP) :: ERR_TENS_2 ! excessive extrapolation: tension storage in layer2 (mm day-1) - REAL(SP) :: ERR_FREE_2 ! excessive extrapolation: free storage in layer2 (mm day-1) - REAL(SP) :: ERR_FREE_2A ! excessive extrapolation: storage in the primary resvr (mm day-1) - REAL(SP) :: ERR_FREE_2B ! excessive extrapolation: storage in the secondary resvr (mm day-1) - REAL(SP) :: CHK_TIME ! time elapsed during time step (days) - ENDTYPE FLUXES - TYPE(FLUXES) :: M_FLUX ! model fluxes - TYPE(FLUXES) :: FLUX_0 ! model fluxes at start of step - TYPE(FLUXES) :: FLUX_1 ! model fluxes at end of step - TYPE(FLUXES), DIMENSION(:), POINTER :: FDFLUX=>NULL() ! finite difference fluxes - TYPE(FLUXES) :: W_FLUX ! weighted sum of model fluxes over a time step - TYPE(FLUXES), dimension(:,:,:), allocatable :: W_FLUX_3d ! weighted sum of model fluxes over a time step for several time steps - REAL(SP) :: CURRENT_DT ! current time step (days) -END MODULE multi_flux diff --git a/build/FUSE_SRC/FUSE_ENGINE/multistate.f90 b/build/FUSE_SRC/FUSE_ENGINE/multistate.f90 deleted file mode 100644 index 51c563c..0000000 --- a/build/FUSE_SRC/FUSE_ENGINE/multistate.f90 +++ /dev/null @@ -1,53 +0,0 @@ -MODULE multistate - USE nrtype - ! -------------------------------------------------------------------------------------- - ! model state structure - ! -------------------------------------------------------------------------------------- - TYPE STATEV - ! snow layer - REAL(SP) :: SWE_TOT ! total storage as snow (mm) - ! upper layer - REAL(SP) :: WATR_1 ! total storage in layer1 (mm) - REAL(SP) :: TENS_1 ! tension storage in layer1 (mm) - REAL(SP) :: FREE_1 ! free storage in layer 1 (mm) - REAL(SP) :: TENS_1A ! storage in the recharge zone (mm) - REAL(SP) :: TENS_1B ! storage in the lower zone (mm) - ! lower layer - REAL(SP) :: WATR_2 ! total storage in layer2 (mm) - REAL(SP) :: TENS_2 ! tension storage in layer2 (mm) - REAL(SP) :: FREE_2 ! free storage in layer2 (mm) - REAL(SP) :: FREE_2A ! storage in the primary resvr (mm) - REAL(SP) :: FREE_2B ! storage in the secondary resvr (mm) - END TYPE STATEV - ! -------------------------------------------------------------------------------------- - ! model time structure - ! -------------------------------------------------------------------------------------- - TYPE M_TIME - REAL(SP) :: STEP ! (time interval to advance model states) - END TYPE M_TIME - ! -------------------------------------------------------------------------------------- - ! variable definitions - ! -------------------------------------------------------------------------------------- - type(statev),dimension(:,:),pointer :: gState ! (grid of model states) - type(statev),dimension(:,:,:),pointer :: gState_3d ! (grid of model states with a time dimension) - TYPE(STATEV) :: ASTATE ! (model states at the start of full timestep) - TYPE(STATEV) :: FSTATE ! (model states at start of sub-timestep) - TYPE(STATEV) :: MSTATE ! (model states at start/middle of sub-timestep) - TYPE(STATEV) :: TSTATE ! (temporary copy of model states) - TYPE(STATEV) :: BSTATE ! (temporary copy of model states) - TYPE(STATEV) :: ESTATE ! (temporary copy of model states) - TYPE(STATEV) :: DSTATE ! (default model states) - TYPE(STATEV) :: DYDT_0 ! (derivative of model states at start of sub-step) - TYPE(STATEV) :: DYDT_1 ! (derivative of model states at end of sub-step) - TYPE(STATEV) :: DY_DT ! (derivative of model states) - TYPE(STATEV) :: DYDT_OLD ! (derivative of model states for final solution) - TYPE(M_TIME) :: HSTATE ! (time interval to advance model states) - ! -------------------------------------------------------------------------------------- - - ! NetCDF - integer(i4b) :: ncid_out=-1 ! NetCDF output file ID - - ! initial store fraction (initialization) - real(sp),parameter::fracState0=0.25_sp - -END MODULE multistate diff --git a/build/FUSE_SRC/FUSE_ENGINE/multistats.f90 b/build/FUSE_SRC/FUSE_ENGINE/multistats.f90 deleted file mode 100644 index 74096ca..0000000 --- a/build/FUSE_SRC/FUSE_ENGINE/multistats.f90 +++ /dev/null @@ -1,35 +0,0 @@ -MODULE multistats - USE nrtype - TYPE SUMMARY - ! DMSL diagnostix - REAL(SP) :: VAR_RESIDUL ! variance of the model residuals - REAL(SP) :: LOGP_SIMULN ! log density of the model simulation - REAL(SP) :: JUMP_TAKEN ! defines a jump in the MCMC production run - ! comparisons between model output and observations - REAL(SP) :: QOBS_MEAN ! mean observed runoff (mm day-1) - REAL(SP) :: QSIM_MEAN ! mean simulated runoff (mm day-1) - REAL(SP) :: QOBS_CVAR ! coefficient of variation of observed runoff (-) - REAL(SP) :: QSIM_CVAR ! coefficient of variation of simulated runoff (-) - REAL(SP) :: QOBS_LAG1 ! lag-1 correlation of observed runoff (-) - REAL(SP) :: QSIM_LAG1 ! lag-1 correlation of simulated runoff (-) - REAL(SP) :: RAW_RMSE ! root-mean-squared-error of flow (mm day-1) - REAL(SP) :: LOG_RMSE ! root-mean-squared-error of LOG flow (mm day-1) - REAL(SP) :: NASH_SUTT ! Nash-Sutcliffe score - ! attributes of model output - REAL(SP) :: NUM_RMSE ! error of the approximate solution - REAL(SP) :: NUM_FUNCS ! number of function calls - REAL(SP) :: NUM_JACOBIAN ! number of times Jacobian is calculated - REAL(SP) :: NUMSUB_ACCEPT ! number of sub-steps taken - REAL(SP) :: NUMSUB_REJECT ! number of sub-steps taken - REAL(SP) :: NUMSUB_NOCONV ! number of sub-steps tried that did not converge - INTEGER(I4B) :: MAXNUM_ITERNS ! maximum number of iterations in implicit scheme - REAL(SP), DIMENSION(20) :: NUMSUB_PROB ! probability distribution for number of sub-steps - ! error checking - CHARACTER(LEN=1024) :: ERR_MESSAGE ! error message - ENDTYPE SUMMARY - ! final data structures - TYPE(SUMMARY) :: MSTATS ! (model summary statistics) - INTEGER(I4B) :: MOD_IX=1 ! (model index) - INTEGER(I4B) :: PCOUNT ! (number of parameter sets in model output files) - INTEGER(I4B) :: FCOUNT ! (number of model simulations) -END MODULE multistats diff --git a/build/FUSE_SRC/FUSE_ENGINE/putpar_str.f90 b/build/FUSE_SRC/FUSE_ENGINE/putpar_str.f90 index 139aea3..3481b65 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/putpar_str.f90 +++ b/build/FUSE_SRC/FUSE_ENGINE/putpar_str.f90 @@ -13,7 +13,8 @@ SUBROUTINE PUTPAR_STR(METADAT,PARNAME) ! Inserts parameter metadata into data structures ! --------------------------------------------------------------------------------------- USE nrtype ! variable types, etc. -USE multiparam, ONLY : PARATT, PARMETA ! derived type for parameter metadata +USE data_types, ONLY : PARATT ! derived type for parameter metadata +USE multiparam, ONLY : PARMETA ! parameter metadata IMPLICIT NONE ! input TYPE(PARATT), INTENT(IN) :: METADAT ! parameter metadata diff --git a/build/FUSE_SRC/FUSE_ENGINE/q_misscell.f90 b/build/FUSE_SRC/FUSE_ENGINE/q_misscell.f90 index b40328a..b84bad5 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/q_misscell.f90 +++ b/build/FUSE_SRC/FUSE_ENGINE/q_misscell.f90 @@ -20,9 +20,10 @@ SUBROUTINE Q_MISSCELL() USE nrtype ! variable types, etc. USE model_defn ! model definition structure USE model_defnames +USE globaldata, ONLY: CURRENT_DT USE multiparam, ONLY: MPARAM,DPARAM ! model parameters USE multistate, ONLY: MSTATE,TSTATE ! model states -USE multi_flux, ONLY: M_FLUX,CURRENT_DT ! model fluxes +USE multi_flux, ONLY: M_FLUX ! model fluxes USE model_numerix ! access model numerix decisions IMPLICIT NONE REAL(SP) :: LOGISMOOTH ! FUNCTION logistic smoothing diff --git a/build/FUSE_SRC/FUSE_ENGINE/str_2_xtry.f90 b/build/FUSE_SRC/FUSE_ENGINE/str_2_xtry.f90 index cb0ac71..71875e9 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/str_2_xtry.f90 +++ b/build/FUSE_SRC/FUSE_ENGINE/str_2_xtry.f90 @@ -14,7 +14,7 @@ SUBROUTINE STR_2_XTRY(TMPSTR,X_TRY) USE nrtype ! Numerical Recipes data types USE model_defn, ONLY: CSTATE,NSTATE ! model definitions USE model_defnames -USE multistate, ONLY: STATEV ! model state structure +USE data_types, ONLY: STATEV ! model state structure IMPLICIT NONE ! input TYPE(STATEV), INTENT(IN) :: TMPSTR ! temporary state structure diff --git a/build/FUSE_SRC/FUSE_NETCDF/def_output.f90 b/build/FUSE_SRC/FUSE_NETCDF/def_output.f90 index 6e71c93..ebd3dc9 100644 --- a/build/FUSE_SRC/FUSE_NETCDF/def_output.f90 +++ b/build/FUSE_SRC/FUSE_NETCDF/def_output.f90 @@ -19,7 +19,7 @@ SUBROUTINE DEF_OUTPUT(nSpat1,nSpat2,NPSET,NTIM) USE multiforce, only: name_psets,time_steps ! dimension arrays USE multiforce, only: latUnits,lonUnits ! units string USE multiforce, only: timeUnits ! units string - USE multistate, only: ncid_out ! NetCDF output file ID + USE globaldata, only: ncid_out ! NetCDF output file ID IMPLICIT NONE @@ -64,9 +64,9 @@ SUBROUTINE DEF_OUTPUT(nSpat1,nSpat2,NPSET,NTIM) !IERR = NF_REDEF(ncid_out); CALL HANDLE_ERR(IERR) ! define dimensions - IERR = NF_DEF_DIM(ncid_out,'time',NF_UNLIMITED,NTIM_DIM); CALL HANDLE_ERR(IERR) !record dimension (unlimited length) - IERR = NF_DEF_DIM(ncid_out,'longitude',nSpat1,lon_dim); CALL HANDLE_ERR(IERR) - IERR = NF_DEF_DIM(ncid_out,'latitude',nSpat2,lat_dim); CALL HANDLE_ERR(IERR) + IERR = NF_DEF_DIM(ncid_out,'time',NF_UNLIMITED, NTIM_DIM); CALL HANDLE_ERR(IERR) !record dimension (unlimited length) + IERR = NF_DEF_DIM(ncid_out,'longitude',nSpat1, lon_dim); CALL HANDLE_ERR(IERR) + IERR = NF_DEF_DIM(ncid_out,'latitude',nSpat2, lat_dim); CALL HANDLE_ERR(IERR) IF(.NOT.GRID_FLAG)THEN IERR = NF_DEF_DIM(ncid_out,'param_set',NPSET,param_dim); CALL HANDLE_ERR(IERR) ENDIF @@ -141,23 +141,23 @@ SUBROUTINE DEF_OUTPUT(nSpat1,nSpat2,NPSET,NTIM) END DO ! ivar ! define the time variable - ierr = nf_def_var(ncid_out,'time',nf_real,1,ntim_dim,ivar_id); call handle_err(ierr) + ierr = nf_def_var(ncid_out,'time',nf_real,1,(/ntim_dim/),ivar_id); call handle_err(ierr) ierr = nf_put_att_text(ncid_out,ivar_id,'units',len_trim(timeUnits),trim(timeUnits)) call handle_err(ierr) ! define the latitude variable - ierr = nf_def_var(ncid_out,'latitude',nf_real,1,lat_dim,ivar_id); call handle_err(ierr) + ierr = nf_def_var(ncid_out,'latitude',nf_real,1,(/lat_dim/),ivar_id); call handle_err(ierr) ierr = nf_put_att_text(ncid_out,ivar_id,'units',8,'degreesN'); call handle_err(ierr) ierr = nf_put_att_text(ncid_out,ivar_id,'axis',1,'Y'); call handle_err(ierr) ! define the longitude variable - ierr = nf_def_var(ncid_out,'longitude',nf_real,1,lon_dim,ivar_id); call handle_err(ierr) + ierr = nf_def_var(ncid_out,'longitude',nf_real,1,(/lon_dim/),ivar_id); call handle_err(ierr) ierr = nf_put_att_text(ncid_out,ivar_id,'units',8,'degreesE'); call handle_err(ierr) ierr = nf_put_att_text(ncid_out,ivar_id,'axis',1,'X'); call handle_err(ierr) IF(.NOT.GRID_FLAG)THEN ! define the param_set variable - ierr = nf_def_var(ncid_out,'param_set',nf_char,1,param_dim,ivar_id); call handle_err(ierr) + ierr = nf_def_var(ncid_out,'param_set',nf_char,1,(/param_dim/),ivar_id); call handle_err(ierr) ierr = nf_put_att_text(ncid_out,ivar_id,'units',1,'-'); call handle_err(ierr) ENDIF diff --git a/build/FUSE_SRC/FUSE_NETCDF/def_params.f90 b/build/FUSE_SRC/FUSE_NETCDF/def_params.f90 index 46b2cdb..0c9ea24 100644 --- a/build/FUSE_SRC/FUSE_NETCDF/def_params.f90 +++ b/build/FUSE_SRC/FUSE_NETCDF/def_params.f90 @@ -13,7 +13,7 @@ SUBROUTINE DEF_PARAMS(NPAR) USE model_defn ! model definition (includes filename) USE metaparams ! metadata for all model parameters USE multistats, ONLY: MSTATS ! model statistics structure -USE multistate, only: ncid_out ! NetCDF output file ID +USE globaldata, only: ncid_out ! NetCDF output file ID IMPLICIT NONE ! input INTEGER(I4B), INTENT(IN) :: NPAR ! number of parameter sets diff --git a/build/FUSE_SRC/FUSE_NETCDF/def_sstats.f90 b/build/FUSE_SRC/FUSE_NETCDF/def_sstats.f90 index 2d4b4ac..f75b1e2 100644 --- a/build/FUSE_SRC/FUSE_NETCDF/def_sstats.f90 +++ b/build/FUSE_SRC/FUSE_NETCDF/def_sstats.f90 @@ -12,7 +12,7 @@ SUBROUTINE DEF_SSTATS() USE model_defn ! model definition (includes filename) USE meta_stats ! metadata for summary statistics USE model_numerix ! model numerix decisions -USE multistate, only: ncid_out ! NetCDF output file ID +USE globaldata, only: ncid_out ! NetCDF output file ID IMPLICIT NONE ! internal INTEGER(I4B) :: IERR ! error code; NetCDF ID diff --git a/build/FUSE_SRC/FUSE_NETCDF/get_gforce.f90 b/build/FUSE_SRC/FUSE_NETCDF/get_gforce.f90 index f990a04..1445c3e 100644 --- a/build/FUSE_SRC/FUSE_NETCDF/get_gforce.f90 +++ b/build/FUSE_SRC/FUSE_NETCDF/get_gforce.f90 @@ -2,6 +2,7 @@ module get_gforce_module USE nrtype USE netcdf USE time_io +USE globaldata, only: NA_VALUE, NA_VALUE_SP ! missing value implicit none private public::read_ginfo @@ -75,7 +76,6 @@ SUBROUTINE read_ginfo(ncid,ierr,message) USE multiforce,only:latUnits,lonUnits,timeUnits ! units string for time USE multiforce,only:vname_dtime ! variable name: time sice reference time USE multiforce, only: nForce, nInput ! number of parameter set and their names - USE multiforce, only: NA_VALUE ! NA_VALUE for the forcing #ifdef __MPI__ use mpi diff --git a/build/FUSE_SRC/FUSE_NETCDF/put_output.f90 b/build/FUSE_SRC/FUSE_NETCDF/put_output.f90 index bc2e361..99d6676 100644 --- a/build/FUSE_SRC/FUSE_NETCDF/put_output.f90 +++ b/build/FUSE_SRC/FUSE_NETCDF/put_output.f90 @@ -15,7 +15,7 @@ SUBROUTINE PUT_OUTPUT(iSpat1,iSpat2,ITIM,IMOD,IPAR) USE varextract_module ! interface for the function to extract variables USE fuse_fileManager,only: Q_ONLY ! only write streamflow to output file? USE multiforce,ONLY: timDat ! time data - USE multistate, only: ncid_out ! NetCDF output file ID + USE globaldata, only: ncid_out ! NetCDF output file ID IMPLICIT NONE ! input @@ -71,7 +71,7 @@ SUBROUTINE PUT_OUTPUT(iSpat1,iSpat2,ITIM,IMOD,IPAR) ! write the time tDat = timDat%dtime ! convert to actual single precision ierr = nf_inq_varid(ncid_out,'time',ivar_id); CALL handle_err(ierr) ! get variable ID for time - ierr = nf_put_var1_real(ncid_out,ivar_id,itim,tDat); CALL handle_err(ierr) ! write time variable + ierr = nf_put_var1_real(ncid_out,ivar_id,(/itim/),tDat); CALL handle_err(ierr) ! write time variable ! close NetCDF file IERR = NF_CLOSE(ncid_out) @@ -95,10 +95,10 @@ SUBROUTINE PUT_GOUTPUT_3D(istart_sim,istart_in,numtim,IPSET) USE fuse_fileManager,only: Q_ONLY ! only write streamflow to output file? USE multiforce, ONLY: timDat,time_steps ! time data - USE multistate, only: ncid_out ! NetCDF output file ID USE multiforce, ONLY: nspat1,nspat2,startSpat2 ! spatial dimensions USE multiforce, ONLY: gForce_3d ! test only - USE multiforce, only: GRID_FLAG ! .true. if distributed + USE multiforce, only: GRID_FLAG ! .true. if distributed + USE globaldata, only: ncid_out ! NetCDF output file ID IMPLICIT NONE @@ -180,7 +180,7 @@ SUBROUTINE PUT_GOUTPUT_3D(istart_sim,istart_in,numtim,IPSET) time_steps_sub = time_steps(istart_in:(istart_in+numtim-1)) ! extract time for subperiod tDat = time_steps_sub ! convert to actual single precision ierr = nf_inq_varid(ncid_out,'time',ivar_id); CALL handle_err(ierr) ! get variable ID for time - ierr = nf_put_vara_real(ncid_out,ivar_id,istart_sim,numtim,tDat); CALL handle_err(ierr) ! write time variable + ierr = nf_put_vara_real(ncid_out,ivar_id,(/istart_sim/),(/numtim/),tDat); CALL handle_err(ierr) ! write time variable ! close NetCDF file IERR = NF_CLOSE(ncid_out) diff --git a/build/FUSE_SRC/FUSE_ENGINE/multiparam.f90 b/build/FUSE_SRC/dshare/data_types.f90 similarity index 52% rename from build/FUSE_SRC/FUSE_ENGINE/multiparam.f90 rename to build/FUSE_SRC/dshare/data_types.f90 index dd1188e..b27a0f7 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/multiparam.f90 +++ b/build/FUSE_SRC/dshare/data_types.f90 @@ -1,15 +1,120 @@ -! --------------------------------------------------------------------------------------- -! Creator: -! -------- -! Martyn Clark -! Modified by Brian Henn to include snow model, 6/2013 -! --------------------------------------------------------------------------------------- -MODULE multiparam - USE nrtype - USE model_defn,ONLY:NTDH_MAX - ! -------------------------------------------------------------------------------------- - ! (1) PARAMETER METADATA +module data_types + + use nrtype + use model_defn, only:NTDH_MAX + ! -------------------------------------------------------------------------------------- + ! model time structure + ! -------------------------------------------------------------------------------------- + TYPE M_TIME + REAL(SP) :: STEP ! (time interval to advance model states) + END TYPE M_TIME + + ! -------------------------------------------------------------------------------------- + ! model forcing structures + ! -------------------------------------------------------------------------------------- + + ! the time data structure (will have no spatial dimension) + TYPE TDATA + INTEGER(I4B) :: IY ! year + INTEGER(I4B) :: IM ! month + INTEGER(I4B) :: ID ! day + INTEGER(I4B) :: IH ! hour + INTEGER(I4B) :: IMIN ! minute + REAL(SP) :: DSEC ! second + REAL(SP) :: DTIME ! time in seconds since year dot + ENDTYPE TDATA + + ! the response structure (will not have a spatial dimension) + TYPE VDATA + REAL(SP) :: OBSQ ! observed runoff (mm day-1) + END TYPE VDATA + + ! ancillary forcing variables used to compute ET (will have a spatial dimension) + TYPE ADATA + REAL(SP) :: AIRTEMP ! air temperature (K) + REAL(SP) :: SPECHUM ! specific humidity (g/g) + REAL(SP) :: AIRPRES ! air pressure (Pa) + REAL(SP) :: SWDOWN ! downward sw radiation (W m-2) + REAL(SP) :: NETRAD ! net radiation (W m-2) + 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) :: TEMP ! temperature for snow model (deg.C) + REAL(SP) :: PET ! energy input: potential ET (mm day-1) + ENDTYPE FDATA + + ! -------------------------------------------------------------------------------------- + ! model state structure + ! -------------------------------------------------------------------------------------- + TYPE STATEV + ! snow layer + REAL(SP) :: SWE_TOT ! total storage as snow (mm) + ! upper layer + REAL(SP) :: WATR_1 ! total storage in layer1 (mm) + REAL(SP) :: TENS_1 ! tension storage in layer1 (mm) + REAL(SP) :: FREE_1 ! free storage in layer 1 (mm) + REAL(SP) :: TENS_1A ! storage in the recharge zone (mm) + REAL(SP) :: TENS_1B ! storage in the lower zone (mm) + ! lower layer + REAL(SP) :: WATR_2 ! total storage in layer2 (mm) + REAL(SP) :: TENS_2 ! tension storage in layer2 (mm) + REAL(SP) :: FREE_2 ! free storage in layer2 (mm) + REAL(SP) :: FREE_2A ! storage in the primary resvr (mm) + REAL(SP) :: FREE_2B ! storage in the secondary resvr (mm) + END TYPE STATEV + + ! -------------------------------------------------------------------------------------- + ! model flux structure + ! -------------------------------------------------------------------------------------- + TYPE FLUXES + REAL(SP) :: EFF_PPT ! effective precipitation (mm day-1) + REAL(SP) :: SATAREA ! saturated area (-) + REAL(SP) :: QSURF ! surface runoff (mm day-1) + REAL(SP) :: EVAP_1A ! evaporation from soil excess zone (mm day-1) + REAL(SP) :: EVAP_1B ! evaporation from soil recharge zone (mm day-1) + REAL(SP) :: EVAP_1 ! evaporation from upper soil layer (mm day-1) + REAL(SP) :: EVAP_2 ! evaporation from lower soil layer (mm day-1) + REAL(SP) :: RCHR2EXCS ! flow from recharge to excess (mm day-1) + REAL(SP) :: TENS2FREE_1 ! flow from tension storage to free storage (mm day-1) + REAL(SP) :: TENS2FREE_2 ! flow from tension storage to free storage (mm day-1) + REAL(SP) :: QINTF_1 ! interflow from free water (mm day-1) + REAL(SP) :: QPERC_12 ! percolation from upper to lower soil layers (mm day-1) + REAL(SP) :: QBASE_2 ! baseflow (mm day-1) + REAL(SP) :: QBASE_2A ! baseflow from primary linear resvr (mm day-1) + REAL(SP) :: QBASE_2B ! baseflow from secondary linear resvr (mm day-1) + REAL(SP) :: OFLOW_1 ! bucket overflow (mm day-1) + REAL(SP) :: OFLOW_2 ! bucket overflow (mm day-1) + REAL(SP) :: OFLOW_2A ! bucket overflow (mm day-1) + REAL(SP) :: OFLOW_2B ! bucket overflow (mm day-1) + REAL(SP) :: ERR_WATR_1 ! excessive extrapolation: total storage in layer1 (mm day-1) + REAL(SP) :: ERR_TENS_1 ! excessive extrapolation: tension storage in layer1 (mm day-1) + REAL(SP) :: ERR_FREE_1 ! excessive extrapolation: free storage in layer 1 (mm day-1) + REAL(SP) :: ERR_TENS_1A ! excessive extrapolation: storage in the recharge zone (mm day-1) + REAL(SP) :: ERR_TENS_1B ! excessive extrapolation: storage in the lower zone (mm day-1) + REAL(SP) :: ERR_WATR_2 ! excessive extrapolation: total storage in layer2 (mm day-1) + REAL(SP) :: ERR_TENS_2 ! excessive extrapolation: tension storage in layer2 (mm day-1) + REAL(SP) :: ERR_FREE_2 ! excessive extrapolation: free storage in layer2 (mm day-1) + REAL(SP) :: ERR_FREE_2A ! excessive extrapolation: storage in the primary resvr (mm day-1) + REAL(SP) :: ERR_FREE_2B ! excessive extrapolation: storage in the secondary resvr (mm day-1) + REAL(SP) :: CHK_TIME ! time elapsed during time step (days) + ENDTYPE FLUXES + + ! -------------------------------------------------------------------------------------- + ! model runoff structure + ! -------------------------------------------------------------------------------------- + TYPE RUNOFF + REAL(SP) :: Q_INSTNT ! instantaneous runoff + REAL(SP) :: Q_ROUTED ! routed runoff + REAL(SP) :: Q_ACCURATE ! "accurate" runoff estimate (mm day-1) + END TYPE RUNOFF + + ! -------------------------------------------------------------------------------------- + ! parameter metadata + ! -------------------------------------------------------------------------------------- + ! data structure to hold metadata for adjustable model parameters TYPE PARATT LOGICAL(LGT) :: PARFIT ! flag to determine if parameter is fitted @@ -29,6 +134,7 @@ MODULE multiparam CHARACTER(LEN=256) :: CHILD1 ! name of 1st parameter child CHARACTER(LEN=256) :: CHILD2 ! name of 2nd parameter child END TYPE PARATT + ! data structure to hold metadata for each parameter TYPE PARINFO ! rainfall error parameters (adjustable) @@ -48,7 +154,7 @@ MODULE multiparam TYPE(PARATT) :: FPRIMQB ! SAC: fraction of baseflow in primary resvr (-) ! evaporation (adjustable) TYPE(PARATT) :: RTFRAC1 ! fraction of roots in the upper layer (-) - ! percolation (adjustable) + ! percolation (adjustable) TYPE(PARATT) :: PERCRTE ! percolation rate (mm day-1) TYPE(PARATT) :: PERCEXP ! percolation exponent (-) TYPE(PARATT) :: SACPMLT ! multiplier in the SAC model for dry lower layer (-) @@ -75,11 +181,12 @@ MODULE multiparam TYPE(PARATT) :: MFMAX ! maximum melt factor (mm melt deg C.-1 6hrs-1) TYPE(PARATT) :: MFMIN ! minimum melt factor (mm melt deg C.-1 6hrs-1) TYPE(PARATT) :: PXTEMP ! rain-snow partition temperature (deg. C) - TYPE(PARATT) :: OPG ! precipitation gradient (-) + TYPE(PARATT) :: OPG ! precipitation gradient (-) TYPE(PARATT) :: LAPSE ! temperature gradient (deg. C) ENDTYPE PARINFO + ! -------------------------------------------------------------------------------------- - ! (2) ADJUSTABLE PARAMETERS + ! adjustable parameters ! -------------------------------------------------------------------------------------- TYPE PARADJ ! rainfall error parameters (adjustable) @@ -126,11 +233,12 @@ MODULE multiparam REAL(SP) :: MFMAX ! maximum melt factor (mm melt deg C.-1 6hrs-1) REAL(SP) :: MFMIN ! minimum melt factor (mm melt deg C.-1 6hrs-1) REAL(SP) :: PXTEMP ! rain-snow partition temperature (deg. C) - REAL(SP) :: OPG ! precipitation gradient (-) + REAL(SP) :: OPG ! precipitation gradient (-) REAL(SP) :: LAPSE ! temperature gradient (deg. C) END TYPE PARADJ + ! -------------------------------------------------------------------------------------- - ! (3) DERIVED PARAMETERS + ! derived parameters ! -------------------------------------------------------------------------------------- TYPE PARDVD ! bucket sizes (derived) @@ -153,22 +261,61 @@ MODULE multiparam REAL(SP), DIMENSION(NTDH_MAX) :: FRAC_FUTURE ! fraction of runoff in future time steps INTEGER(I4B) :: NTDH_NEED ! number of time-steps with non-zero routing contribution END TYPE PARDVD + ! -------------------------------------------------------------------------------------- - ! (4) LIST OF PARAMETERS FOR A GIVEN MODEL + ! list of parameters for a given model ! -------------------------------------------------------------------------------------- TYPE PAR_ID CHARACTER(LEN=9) :: PARNAME ! list of parameter names ENDTYPE PAR_ID + + ! -------------------------------------------------------------------------------------- + ! model statistics structure ! -------------------------------------------------------------------------------------- - ! (5) FINAL DATA STRUCTURES + TYPE SUMMARY + ! DMSL diagnostix + REAL(SP) :: VAR_RESIDUL ! variance of the model residuals + REAL(SP) :: LOGP_SIMULN ! log density of the model simulation + REAL(SP) :: JUMP_TAKEN ! defines a jump in the MCMC production run + ! comparisons between model output and observations + REAL(SP) :: QOBS_MEAN ! mean observed runoff (mm day-1) + REAL(SP) :: QSIM_MEAN ! mean simulated runoff (mm day-1) + REAL(SP) :: QOBS_CVAR ! coefficient of variation of observed runoff (-) + REAL(SP) :: QSIM_CVAR ! coefficient of variation of simulated runoff (-) + REAL(SP) :: QOBS_LAG1 ! lag-1 correlation of observed runoff (-) + REAL(SP) :: QSIM_LAG1 ! lag-1 correlation of simulated runoff (-) + REAL(SP) :: RAW_RMSE ! root-mean-squared-error of flow (mm day-1) + REAL(SP) :: LOG_RMSE ! root-mean-squared-error of LOG flow (mm day-1) + REAL(SP) :: NASH_SUTT ! Nash-Sutcliffe score + ! attributes of model output + REAL(SP) :: NUM_RMSE ! error of the approximate solution + REAL(SP) :: NUM_FUNCS ! number of function calls + REAL(SP) :: NUM_JACOBIAN ! number of times Jacobian is calculated + REAL(SP) :: NUMSUB_ACCEPT ! number of sub-steps taken + REAL(SP) :: NUMSUB_REJECT ! number of sub-steps taken + REAL(SP) :: NUMSUB_NOCONV ! number of sub-steps tried that did not converge + INTEGER(I4B) :: MAXNUM_ITERNS ! maximum number of iterations in implicit scheme + REAL(SP), DIMENSION(20) :: NUMSUB_PROB ! probability distribution for number of sub-steps + ! error checking + CHARACTER(LEN=1024) :: ERR_MESSAGE ! error message + ENDTYPE SUMMARY + ! -------------------------------------------------------------------------------------- - INTEGER(I4B), PARAMETER :: MAXPAR=50 ! maximum number of parameters for a single model - TYPE(PARADJ), DIMENSION(:), POINTER :: APARAM=>null() ! all model parameter sets; DK/2008/10/21: explicit null - TYPE(PARADJ) :: MPARAM ! single model parameter set - TYPE(PARDVD) :: DPARAM ! derived model parameters - TYPE(PARINFO) :: PARMETA ! parameter metadata (all parameters) - TYPE(PAR_ID), DIMENSION(MAXPAR) :: LPARAM ! list of model parameter names (need to modify to 16 for SCE) - INTEGER(I4B) :: NUMPAR ! number of model parameters for current model - INTEGER(I4B) :: SOBOL_INDX ! code to re-assemble Sobol parameters + ! parent FUSE structure ! -------------------------------------------------------------------------------------- -END MODULE multiparam + type parent + type(m_time) :: time ! time step + type(fdata) :: force ! model forcing data + type(statev) :: state0 ! state variables (start of step) + type(statev) :: state1 ! state variables (end of step) + type(statev) :: dx_dt ! time derivative in state variables + type(fluxes) :: flux ! fluxes + type(runoff) :: route ! hillslope routing + type(par_id) :: param_name ! parameter names + type(parinfo) :: param_meta ! metadata on model parameters + type(paradj) :: param_adjust ! adjustable model parametrs + type(pardvd) :: param_derive ! derived model parameters + type(summary) :: sim_stats ! simulation statistics + end type parent + +end module data_types diff --git a/build/FUSE_SRC/dshare/globaldata.f90 b/build/FUSE_SRC/dshare/globaldata.f90 new file mode 100644 index 0000000..0029c83 --- /dev/null +++ b/build/FUSE_SRC/dshare/globaldata.f90 @@ -0,0 +1,21 @@ +MODULE globaldata + + USE nrtype + + ! time step + REAL(SP), save :: CURRENT_DT ! current time step (days) + + ! missing values + INTEGER(I4B),PARAMETER :: NA_VALUE=-9999 ! integer designating missing values - TODO: retrieve from NetCDF file + REAL(SP),PARAMETER :: NA_VALUE_SP=-9999_sp ! integer designating missing values - TODO: retrieve from NetCDF file + + ! NetCDF + integer(i4b), save :: ncid_out=-1 ! NetCDF output file ID + + ! initial store fraction (initialization) + real(sp), parameter :: fracState0=0.25_sp + + ! print flag + logical(lgt) :: isPrint=.true. + +end MODULE globaldata diff --git a/build/FUSE_SRC/FUSE_ENGINE/model_defn.f90 b/build/FUSE_SRC/dshare/model_defn.f90 similarity index 81% rename from build/FUSE_SRC/FUSE_ENGINE/model_defn.f90 rename to build/FUSE_SRC/dshare/model_defn.f90 index 9a0c80a..0dcd28b 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/model_defn.f90 +++ b/build/FUSE_SRC/dshare/model_defn.f90 @@ -27,36 +27,27 @@ MODULE model_defn TYPE UMODEL INTEGER(I4B) :: MODIX ! model index CHARACTER(LEN=256) :: MNAME ! model name -! CHARACTER(LEN=16) :: RFERR ! rainfall error INTEGER(I4B) :: iRFERR -! CHARACTER(LEN=16) :: ARCH1 ! upper-layer architecture INTEGER(I4B) :: iARCH1 -! CHARACTER(LEN=16) :: ARCH2 ! lower-layer architecture INTEGER(I4B) :: iARCH2 -! CHARACTER(LEN=16) :: QSURF ! surface runoff INTEGER(I4B) :: iQSURF -! CHARACTER(LEN=16) :: QPERC ! percolation INTEGER(I4B) :: iQPERC -! CHARACTER(LEN=16) :: ESOIL ! evaporation INTEGER(I4B) :: iESOIL -! CHARACTER(LEN=16) :: QINTF ! interflow INTEGER(I4B) :: iQINTF -! CHARACTER(LEN=16) :: Q_TDH ! time delay in runoff INTEGER(I4B) :: iQ_TDH INTEGER(I4B) :: iSNOWM ! snow - END TYPE UMODEL + END TYPE UMODEL ! structure to hold model state names TYPE SNAMES -! CHARACTER(LEN=8) :: SNAME ! state name INTEGER(I4B) :: iSNAME ! integer value of state name END TYPE SNAMES ! structure to hold model flux names TYPE FNAMES CHARACTER(LEN=16) :: FNAME ! state name END TYPE FNAMES -! max steps in routing function - INTEGER(I4B),PARAMETER::NTDH_MAX=500 -! model definitions + ! max steps in routing function + INTEGER(I4B),PARAMETER::NTDH_MAX=500 + ! model definitions CHARACTER(LEN=256) :: FNAME_NETCDF_RUNS ! NETCDF output filename for model runs CHARACTER(LEN=256) :: FNAME_NETCDF_PARA ! NETCDF output filename for model parameters CHARACTER(LEN=256) :: FNAME_NETCDF_PARA_SCE ! NETCDF output filename for model parameters produced by SCE diff --git a/build/FUSE_SRC/FUSE_ENGINE/model_defnames.f90 b/build/FUSE_SRC/dshare/model_defnames.f90 similarity index 100% rename from build/FUSE_SRC/FUSE_ENGINE/model_defnames.f90 rename to build/FUSE_SRC/dshare/model_defnames.f90 diff --git a/build/FUSE_SRC/FUSE_ENGINE/model_numerix.f90 b/build/FUSE_SRC/dshare/model_numerix.f90 similarity index 96% rename from build/FUSE_SRC/FUSE_ENGINE/model_numerix.f90 rename to build/FUSE_SRC/dshare/model_numerix.f90 index 8aefa42..030073e 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/model_numerix.f90 +++ b/build/FUSE_SRC/dshare/model_numerix.f90 @@ -30,6 +30,9 @@ MODULE model_numerix ! 6. Method used to process the small interval at the end of a time step INTEGER(I4B), PARAMETER :: STEP_TRUNC=0, LOOK_AHEAD=1, STEP_ABSORB=2 INTEGER(I4B) :: SMALL_ENDSTEP +! 7. Flag for differentiable model +integer(i4b), parameter :: original=0, differentiable=1 +integer(i4b) :: diff_mode ! --------------------------------------------------------------------------------------- ! (B) PARAMETERS ! --------------------------------------------------------------------------------------- diff --git a/build/FUSE_SRC/dshare/multi_flux.f90 b/build/FUSE_SRC/dshare/multi_flux.f90 new file mode 100644 index 0000000..b00bb06 --- /dev/null +++ b/build/FUSE_SRC/dshare/multi_flux.f90 @@ -0,0 +1,10 @@ +MODULE multi_flux + USE nrtype + use data_types, only: fluxes + TYPE(FLUXES) :: M_FLUX ! model fluxes + TYPE(FLUXES) :: FLUX_0 ! model fluxes at start of step + TYPE(FLUXES) :: FLUX_1 ! model fluxes at end of step + TYPE(FLUXES), DIMENSION(:), POINTER :: FDFLUX=>NULL() ! finite difference fluxes + TYPE(FLUXES) :: W_FLUX ! weighted sum of model fluxes over a time step + TYPE(FLUXES), dimension(:,:,:), allocatable :: W_FLUX_3d ! weighted sum of model fluxes over a time step for several time steps +END MODULE multi_flux diff --git a/build/FUSE_SRC/FUSE_ENGINE/multibands.f90 b/build/FUSE_SRC/dshare/multibands.f90 similarity index 100% rename from build/FUSE_SRC/FUSE_ENGINE/multibands.f90 rename to build/FUSE_SRC/dshare/multibands.f90 diff --git a/build/FUSE_SRC/FUSE_ENGINE/multiconst.f90 b/build/FUSE_SRC/dshare/multiconst.f90 similarity index 100% rename from build/FUSE_SRC/FUSE_ENGINE/multiconst.f90 rename to build/FUSE_SRC/dshare/multiconst.f90 diff --git a/build/FUSE_SRC/FUSE_ENGINE/multiforce.f90 b/build/FUSE_SRC/dshare/multiforce.f90 similarity index 82% rename from build/FUSE_SRC/FUSE_ENGINE/multiforce.f90 rename to build/FUSE_SRC/dshare/multiforce.f90 index 900befd..468f649 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/multiforce.f90 +++ b/build/FUSE_SRC/dshare/multiforce.f90 @@ -4,40 +4,13 @@ ! Martyn Clark ! Modified by Brian Henn to include snow model, 6/2013 ! Modified by Nans Addor to enable distributed modeling, 9/2016 +! Modified by Martyn Clark to separate derived types from shard data, 12/2025 ! --------------------------------------------------------------------------------------- MODULE multiforce USE nrtype + USE data_types, only: tdata, vdata, adata, fdata SAVE ! -------------------------------------------------------------------------------------- - ! the time data structure (will have no spatial dimension) - TYPE TDATA - INTEGER(I4B) :: IY ! year - INTEGER(I4B) :: IM ! month - INTEGER(I4B) :: ID ! day - INTEGER(I4B) :: IH ! hour - INTEGER(I4B) :: IMIN ! minute - REAL(SP) :: DSEC ! second - REAL(SP) :: DTIME ! time in seconds since year dot - ENDTYPE TDATA - ! the response structure (will not have a spatial dimension) - TYPE VDATA - REAL(SP) :: OBSQ ! observed runoff (mm day-1) - END TYPE VDATA - ! ancillary forcing variables used to compute ET (will have a spatial dimension) - TYPE ADATA - REAL(SP) :: AIRTEMP ! air temperature (K) - REAL(SP) :: SPECHUM ! specific humidity (g/g) - REAL(SP) :: AIRPRES ! air pressure (Pa) - REAL(SP) :: SWDOWN ! downward sw radiation (W m-2) - REAL(SP) :: NETRAD ! net radiation (W m-2) - 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) :: TEMP ! temperature for snow model (deg.C) - REAL(SP) :: PET ! energy input: potential ET (mm day-1) - ENDTYPE FDATA - ! -------------------------------------------------------------------------------------- ! general INTEGER(I4B),PARAMETER :: STRLEN=256 ! length of the character string ! time data structures @@ -151,9 +124,5 @@ MODULE multiforce REAL(sp) :: amult_pet=-1._dp ! convert potential ET to mm/day REAL(sp) :: amult_q=-1._dp ! convert runoff to mm/day - ! missing values - INTEGER(I4B),PARAMETER :: NA_VALUE=-9999 ! integer designating missing values - TODO: retrieve from NetCDF file - REAL(SP),PARAMETER :: NA_VALUE_SP=-9999 ! integer designating missing values - TODO: retrieve from NetCDF file - ! -------------------------------------------------------------------------------------- END MODULE multiforce diff --git a/build/FUSE_SRC/dshare/multiparam.f90 b/build/FUSE_SRC/dshare/multiparam.f90 new file mode 100644 index 0000000..cfaa939 --- /dev/null +++ b/build/FUSE_SRC/dshare/multiparam.f90 @@ -0,0 +1,21 @@ +! --------------------------------------------------------------------------------------- +! Creator: +! -------- +! Martyn Clark +! Modified by Brian Henn to include snow model, 6/2013 +! Modified by Martyn Clark to separate derived types from shard data, 12/2025 +! --------------------------------------------------------------------------------------- +MODULE multiparam + USE nrtype + USE data_types,ONLY:par_id,parinfo,paradj,pardvd + ! -------------------------------------------------------------------------------------- + INTEGER(I4B), PARAMETER :: MAXPAR=50 ! maximum number of parameters for a single model + TYPE(PARADJ), DIMENSION(:), POINTER :: APARAM=>null() ! all model parameter sets; DK/2008/10/21: explicit null + TYPE(PARADJ) :: MPARAM ! single model parameter set + TYPE(PARDVD) :: DPARAM ! derived model parameters + TYPE(PARINFO) :: PARMETA ! parameter metadata (all parameters) + TYPE(PAR_ID), DIMENSION(MAXPAR) :: LPARAM ! list of model parameter names (need to modify to 16 for SCE) + INTEGER(I4B) :: NUMPAR ! number of model parameters for current model + INTEGER(I4B) :: SOBOL_INDX ! code to re-assemble Sobol parameters + ! -------------------------------------------------------------------------------------- +END MODULE multiparam diff --git a/build/FUSE_SRC/FUSE_ENGINE/multiroute.f90 b/build/FUSE_SRC/dshare/multiroute.f90 similarity index 100% rename from build/FUSE_SRC/FUSE_ENGINE/multiroute.f90 rename to build/FUSE_SRC/dshare/multiroute.f90 diff --git a/build/FUSE_SRC/dshare/multistate.f90 b/build/FUSE_SRC/dshare/multistate.f90 new file mode 100644 index 0000000..f7724f0 --- /dev/null +++ b/build/FUSE_SRC/dshare/multistate.f90 @@ -0,0 +1,21 @@ +MODULE multistate + USE nrtype + use data_types, only: statev, m_time ! <— import canonical types + + ! variable definitions + type(statev),dimension(:,:),pointer :: gState ! (grid of model states) + type(statev),dimension(:,:,:),pointer :: gState_3d ! (grid of model states with a time dimension) + TYPE(STATEV) :: ASTATE ! (model states at the start of full timestep) + TYPE(STATEV) :: FSTATE ! (model states at start of sub-timestep) + TYPE(STATEV) :: MSTATE ! (model states at start/middle of sub-timestep) + TYPE(STATEV) :: TSTATE ! (temporary copy of model states) + TYPE(STATEV) :: BSTATE ! (temporary copy of model states) + TYPE(STATEV) :: ESTATE ! (temporary copy of model states) + TYPE(STATEV) :: DSTATE ! (default model states) + TYPE(STATEV) :: DYDT_0 ! (derivative of model states at start of sub-step) + TYPE(STATEV) :: DYDT_1 ! (derivative of model states at end of sub-step) + TYPE(STATEV) :: DY_DT ! (derivative of model states) + TYPE(STATEV) :: DYDT_OLD ! (derivative of model states for final solution) + TYPE(M_TIME) :: HSTATE ! (time interval to advance model states) + +END MODULE multistate diff --git a/build/FUSE_SRC/dshare/multistats.f90 b/build/FUSE_SRC/dshare/multistats.f90 new file mode 100644 index 0000000..70907f7 --- /dev/null +++ b/build/FUSE_SRC/dshare/multistats.f90 @@ -0,0 +1,9 @@ +MODULE multistats + USE nrtype + Use data_types, only: summary + ! final data structures + TYPE(SUMMARY) :: MSTATS ! (model summary statistics) + INTEGER(I4B) :: MOD_IX=1 ! (model index) + INTEGER(I4B) :: PCOUNT ! (number of parameter sets in model output files) + INTEGER(I4B) :: FCOUNT ! (number of model simulations) +END MODULE multistats diff --git a/build/FUSE_SRC/physics/evap_lower_diff.f90 b/build/FUSE_SRC/physics/evap_lower_diff.f90 new file mode 100644 index 0000000..a07a76a --- /dev/null +++ b/build/FUSE_SRC/physics/evap_lower_diff.f90 @@ -0,0 +1,88 @@ +module EVAP_LOWER_DIFF_MODULE + + implicit none + + private + public :: EVAP_LOWER_DIFF + +contains + + SUBROUTINE EVAP_LOWER_DIFF(fuseStruct) + ! ------------------------------------------------------------------------------------------------- + ! Creator: + ! -------- + ! Martyn Clark, 2007 + ! Modified by Martyn Clark to create a differentiable model, 12/25 + ! ------------------------------------------------------------------------------------------------- + ! Purpose: + ! -------- + ! Computes evaporation from the lower soil layer + ! ------------------------------------------------------------------------------------------------- + USE nrtype ! variable types, etc. + USE data_types, only: parent ! fuse parent data type + USE model_defn ! model definition structure + USE model_defnames + IMPLICIT NONE + ! input-output + type(parent), intent(inout) :: fuseStruct ! parent fuse data structure + ! ------------------------------------------------------------------------------------------------- + ! associate variables with elements of data structure + associate(& + MFORCE => fuseStruct%force , & ! model forcing data + M_FLUX => fuseStruct%flux , & ! fluxes + TSTATE => fuseStruct%state1 , & ! trial state variables (end of step) + MPARAM => fuseStruct%param_adjust , & ! adjustable model parameters + DPARAM => fuseStruct%param_derive & ! derived model parameters + ) ! (associate) + ! ------------------------------------------------------------------------------------------------- + + ! --------------------------------------------------------------------------------------- + SELECT CASE(SMODL%iARCH2) ! lower layer architecture + CASE(iopt_tens2pll_2,iopt_fixedsiz_2) + + ! ------------------------------------------------------------------------------------- + SELECT CASE(SMODL%iARCH1) + ! ------------------------------------------------------------------------------------ + CASE(iopt_tension1_1,iopt_onestate_1) ! lower-layer evap is valid + + ! ------------------------------------------------------------------------------------ + ! use different evaporation schemes for the lower layer + ! ----------------------------------------------------- + SELECT CASE(SMODL%iESOIL) + CASE(iopt_sequential) + M_FLUX%EVAP_2 = (MFORCE%PET-M_FLUX%EVAP_1) * (TSTATE%TENS_2/DPARAM%MAXTENS_2) + CASE(iopt_rootweight) + M_FLUX%EVAP_2 = MFORCE%PET * DPARAM%RTFRAC2 * (TSTATE%TENS_2/DPARAM%MAXTENS_2) + CASE DEFAULT + print *, "SMODL%iESOIL must be either iopt_sequential or iopt_rootweight" + END SELECT ! (evaporation schemes) + + ! ------------------------------------------------------------------------------------ + CASE(iopt_tension2_1) ! lower-layer evap is zero + M_FLUX%EVAP_2 = 0._sp + + ! ------------------------------------------------------------------------------------ + CASE DEFAULT + print *, "SMODL%iARCH1 must be iopt_tension2_1, iopt_tension1_1, or iopt_onestate_1" + STOP + + ! ------------------------------------------------------------------------------------ + END SELECT ! (upper-layer architechure) + + ! -------------------------------------------------------------------------------------- + CASE(iopt_unlimfrc_2,iopt_unlimpow_2,iopt_topmdexp_2) + M_FLUX%EVAP_2 = 0._sp + + ! -------------------------------------------------------------------------------------- + CASE DEFAULT + print *, "SMODL%iARCH2 must be iopt_tens2pll_2, iopt_unlimfrc_2, iopt_unlimpow_2" + print *, " iopt_topmdexp_2, or iopt_fixedsiz_2" + STOP + + END SELECT + ! --------------------------------------------------------------------------------------- + + end associate ! end association with variables in the data structures + END SUBROUTINE EVAP_LOWER_DIFF + +end module EVAP_LOWER_DIFF_module diff --git a/build/FUSE_SRC/physics/evap_upper_diff.f90 b/build/FUSE_SRC/physics/evap_upper_diff.f90 new file mode 100644 index 0000000..fa0fd2d --- /dev/null +++ b/build/FUSE_SRC/physics/evap_upper_diff.f90 @@ -0,0 +1,91 @@ +module EVAP_UPPER_DIFF_module + + implicit none + + private + public :: EVAP_UPPER_DIFF + +contains + + SUBROUTINE EVAP_UPPER_DIFF(fuseStruct) + ! ------------------------------------------------------------------------------------------------- + ! Creator: + ! -------- + ! Martyn Clark, 2007 + ! Modified by Martyn Clark to create a differentiable model, 12/25 + ! ------------------------------------------------------------------------------------------------- + ! Purpose: + ! -------- + ! Computes evaporation from the upper soil layer + ! ------------------------------------------------------------------------------------------------- + USE nrtype ! variable types, etc. + USE data_types, only: parent ! fuse parent data type + USE model_defn ! model definition structure + USE model_defnames + IMPLICIT NONE + ! input-output + type(parent), intent(inout) :: fuseStruct ! parent fuse data structure + ! ------------------------------------------------------------------------------------------------- + ! associate variables with elements of data structure + associate(& + MFORCE => fuseStruct%force , & ! model forcing data + M_FLUX => fuseStruct%flux , & ! fluxes + TSTATE => fuseStruct%state1 , & ! trial state variables (end of step) + MPARAM => fuseStruct%param_adjust , & ! adjustable model parameters + DPARAM => fuseStruct%param_derive & ! derived model parameters + ) ! (associate) + ! ------------------------------------------------------------------------------------------------- + + ! --------------------------------------------------------------------------------------- + SELECT CASE(SMODL%iARCH1) ! upper layer architecture + + ! -------------------------------------------------------------------------------------- + CASE(iopt_tension2_1) ! tension storage sub-divided into recharge and excess + ! -------------------------------------------------------------------------------------- + + ! use different evaporation schemes for the upper layer + ! ----------------------------------------------------- + SELECT CASE(SMODL%iESOIL) + CASE(iopt_sequential) + M_FLUX%EVAP_1A = MFORCE%PET * TSTATE%TENS_1A/DPARAM%MAXTENS_1A + M_FLUX%EVAP_1B = (MFORCE%PET - M_FLUX%EVAP_1A) * TSTATE%TENS_1B/DPARAM%MAXTENS_1B + M_FLUX%EVAP_1 = M_FLUX%EVAP_1A + M_FLUX%EVAP_1B + CASE(iopt_rootweight) + M_FLUX%EVAP_1A = MFORCE%PET * MPARAM%RTFRAC1 * TSTATE%TENS_1A/DPARAM%MAXTENS_1A + M_FLUX%EVAP_1B = MFORCE%PET * DPARAM%RTFRAC2 * TSTATE%TENS_1B/DPARAM%MAXTENS_1B + M_FLUX%EVAP_1 = M_FLUX%EVAP_1A + M_FLUX%EVAP_1B + CASE DEFAULT + print *, "SMODL%iESOIL must be either iopt_sequential or iopt_rootweight" + STOP + END SELECT + ! -------------------------------------------------------------------------------------- + CASE(iopt_tension1_1,iopt_onestate_1) ! single tension store or single state + ! -------------------------------------------------------------------------------------- + + ! use different evaporation schemes for the upper layer + ! ----------------------------------------------------- + SELECT CASE(SMODL%iESOIL) + CASE(iopt_sequential) + M_FLUX%EVAP_1A = 0._sp + M_FLUX%EVAP_1B = 0._sp + M_FLUX%EVAP_1 = MFORCE%PET * TSTATE%TENS_1/DPARAM%MAXTENS_1 + CASE(iopt_rootweight) + M_FLUX%EVAP_1A = 0._sp + M_FLUX%EVAP_1B = 0._sp + M_FLUX%EVAP_1 = MFORCE%PET * MPARAM%RTFRAC1 * TSTATE%TENS_1/DPARAM%MAXTENS_1 + CASE DEFAULT + print *, "SMODL%iESOIL must be either iopt_sequential or iopt_rootweight" + END SELECT ! (evaporation schemes) + + ! -------------------------------------------------------------------------------------- + CASE DEFAULT + print *, "SMODL%iARCH1 must be iopt_tension2_1, iopt_tension1_1, or iopt_onestate_1" + STOP + ! -------------------------------------------------------------------------------------- + + END SELECT ! (upper-layer architechure) + + end associate ! end association with variables in the data structures + END SUBROUTINE EVAP_UPPER_DIFF + +end module EVAP_UPPER_DIFF_module diff --git a/build/FUSE_SRC/physics/fix_ovshoot.f90 b/build/FUSE_SRC/physics/fix_ovshoot.f90 new file mode 100644 index 0000000..9e51da3 --- /dev/null +++ b/build/FUSE_SRC/physics/fix_ovshoot.f90 @@ -0,0 +1,139 @@ +module overshoot_module + + USE nrtype ! variable types, etc. + USE data_types, only: parent ! fuse parent data type + USE model_defn, only: CSTATE,NSTATE,SMODL ! model definition structures + USE model_defnames + implicit none + + private + public :: get_bounds + public :: fix_ovshoot + +contains + + ! --------------------------------------------------------------------------------------- + ! --------------------------------------------------------------------------------------- + ! Numerically-stable softplus with sharpness alpha + pure real(sp) function softplus(x, alpha) result(y) + implicit none + real(sp), intent(in) :: x, alpha + real(sp) :: ax + ax = alpha * x + if (ax > 0.0_sp) then + y = (ax + log(1.0_sp + exp(-ax))) / alpha + else + y = log(1.0_sp + exp(ax)) / alpha + end if + end function softplus + ! --------------------------------------------------------------------------------------- + ! --------------------------------------------------------------------------------------- + + ! --------------------------------------------------------------------------------------- + ! --------------------------------------------------------------------------------------- + SUBROUTINE fix_ovshoot(X_TRY, lower, upper) + ! --------------------------------------------------------------------------------------- + ! Creator: + ! -------- + ! Martyn Clark, 2025 + ! --------------------------------------------------------------------------------------- + ! Purpose: + ! -------- + ! Apply soft constraints to model state variables + ! --------------------------------------------------------------------------------------- + ! input/output + REAL(SP), DIMENSION(:), INTENT(INOUT) :: X_TRY ! vector of model states + real(sp), dimension(:), intent(in) :: lower ! lower bound + real(sp), dimension(:), intent(in) :: upper ! upper bound + ! internal + integer(i4b) :: i ! index of model state variable + real(sp), parameter :: alpha=10_sp ! controls sharpness in smoothing + + ! apply soft constraint to model states + do i=1,NSTATE + x_try(i) = lower(i) + softplus(x_try(i)-lower(i), alpha) - softplus(x_try(i)-upper(i), alpha) + end do ! looping through model state variables + + end subroutine fix_ovshoot + + ! --------------------------------------------------------------------------------------- + ! --------------------------------------------------------------------------------------- + SUBROUTINE get_bounds(fuseStruct, lower, upper) + ! --------------------------------------------------------------------------------------- + ! Creator: + ! -------- + ! Martyn Clark, 2007 + ! Modified to return lower and upper bounds by Martyn Clark, 12/2025 + ! --------------------------------------------------------------------------------------- + ! Purpose: + ! -------- + ! Identify lower and upper bounds for the vector of model states + ! --------------------------------------------------------------------------------------- + USE model_numerix ! model numerix + IMPLICIT NONE + ! input/output + type(parent), intent(in) :: fuseStruct ! parent fuse data structure + real(sp), dimension(:), intent(out) :: lower ! lower bound for states + real(sp), dimension(:), intent(out) :: upper ! upper bound for states + ! internal + REAL(SP) :: XMIN ! very small number + INTEGER(I4B) :: ISTT ! loop through model states + ! --------------------------------------------------------------------------------------- + associate(MPARAM => fuseStruct%param_adjust, & ! adjuustable model parameters + DPARAM => fuseStruct%param_derive) ! derived model parameters + ! --------------------------------------------------------------------------------------- + XMIN=FRACSTATE_MIN ! used to avoid zero derivatives + ! --------------------------------------------------------------------------------------- + ! loop through model states + DO ISTT=1,NSTATE + SELECT CASE(CSTATE(ISTT)%iSNAME) + ! upper tanks + CASE (iopt_TENS1A) + lower(ISTT) = XMIN*DPARAM%MAXTENS_1A + upper(ISTT) = DPARAM%MAXTENS_1A + CASE (iopt_TENS1B) + lower(ISTT) = XMIN*DPARAM%MAXTENS_1B + upper(ISTT) = DPARAM%MAXTENS_1B + CASE (iopt_TENS_1) + lower(ISTT) = XMIN*DPARAM%MAXTENS_1 + upper(ISTT) = DPARAM%MAXTENS_1 + CASE (iopt_FREE_1) + lower(ISTT) = XMIN*DPARAM%MAXFREE_1 + upper(ISTT) = DPARAM%MAXFREE_1 + CASE (iopt_WATR_1) + lower(ISTT) = XMIN*MPARAM%MAXWATR_1 + upper(ISTT) = MPARAM%MAXWATR_1 + ! lower tanks + CASE (iopt_TENS_2) + lower(ISTT) = XMIN*DPARAM%MAXTENS_2 + upper(ISTT) = DPARAM%MAXTENS_2 + CASE (iopt_FREE2A) + lower(ISTT) = XMIN*DPARAM%MAXFREE_2A + upper(ISTT) = DPARAM%MAXFREE_2A + CASE (iopt_FREE2B) + lower(ISTT) = XMIN*DPARAM%MAXFREE_2B + upper(ISTT) = DPARAM%MAXFREE_2B + CASE (iopt_WATR_2) + ! *** SET LOWER LIMITS *** + IF (SMODL%iARCH2.NE.iopt_topmdexp_2) THEN + ! enforce lower limit + lower(ISTT) = XMIN*MPARAM%MAXWATR_2 + ELSE + ! MPARAM%MAXWATR_2 is just a scaling parameter, but don't allow stupid values + lower(ISTT) = -MPARAM%MAXWATR_2*10._sp + ENDIF + ! *** SET UPPER LIMITS *** + IF (SMODL%iARCH2.EQ.iopt_tens2pll_2 .OR. SMODL%iARCH2.EQ.iopt_fixedsiz_2) THEN + ! cannot exceed capacity + upper(ISTT) = MPARAM%MAXWATR_2 + ELSE + ! unlimited storage, but make sure the values are still sensible + upper(ISTT) = MPARAM%MAXWATR_2*1000._sp + ENDIF + END SELECT + END DO ! (loop through states) + end associate ! end association with variables in the data structures + ! --------------------------------------------------------------------------------------- + END SUBROUTINE get_bounds + +END MODULE overshoot_module diff --git a/build/FUSE_SRC/physics/get_parent.f90 b/build/FUSE_SRC/physics/get_parent.f90 new file mode 100644 index 0000000..1a79e0d --- /dev/null +++ b/build/FUSE_SRC/physics/get_parent.f90 @@ -0,0 +1,26 @@ +module get_parent_module + use data_types, only: parent + implicit none + +contains + + subroutine get_parent(fuseStruct) + use multiforce, only: mForce + use multistate, only: mState + use multi_flux, only: m_flux + use multiparam, only: parMeta,mParam,dParam + implicit none + type(parent), intent(out) :: fuseStruct + ! populate parent fuse structures + fuseStruct%force = mForce + fuseStruct%state0 = mState + fuseStruct%state1 = mState + fuseStruct%flux = m_flux + fuseStruct%param_meta = parMeta + fuseStruct%param_adjust = mParam + fuseStruct%param_derive = dParam + + end subroutine get_parent + + +end module get_parent_module diff --git a/build/FUSE_SRC/physics/implicit_solve.f90 b/build/FUSE_SRC/physics/implicit_solve.f90 new file mode 100644 index 0000000..6ee2325 --- /dev/null +++ b/build/FUSE_SRC/physics/implicit_solve.f90 @@ -0,0 +1,244 @@ +module implicit_solve_module + + ! data types + use nrtype ! variable types, etc. + use data_types, only: parent ! parent fuse data structure + + ! modules + use xtry_2_str_module ! puts state vector into FUSE state structure + use str_2_xtry_module ! puts FUSE state structure into state vector + + ! global data + use model_defn, only: nState ! number of state variables + use multiforce, only: dt => deltim ! time step + + use model_numerix, only: NUM_FUNCS ! number of function calls + use model_numerix, only: NUM_JACOBIAN ! number of times Jacobian is calculated + + implicit none + + ! provide access to the fuse parent structure + type(parent), pointer, save :: ctx => null() + + private + public :: implicit_solve + + contains + + ! ----- point to the fuse parent structure --------------------------------------------- + + subroutine set_dxdt_context(fuseStruct) + type(parent), target, intent(inout) :: fuseStruct + ctx => fuseStruct + end subroutine set_dxdt_context + + subroutine clear_dxdt_context() + nullify(ctx) + end subroutine clear_dxdt_context + + ! -------------------------------------------------------------------------------------- + ! -------------------------------------------------------------------------------------- + ! -------------------------------------------------------------------------------------- + + ! ----- calculate dx/dt=g(x) ----------------------------------------------------------- + function dx_dt(x_try) result(g_x) + use MOD_DERIVS_DIFF_module, only: MOD_DERIVS_DIFF ! compute dx/dt + implicit none + ! input + real(sp) , intent(in) :: x_try(:) ! trial state vector + ! output + real(sp) :: g_x(size(x_try)) ! dx/dt=g(x) + + ! check made the association to ctx (ctx=>fuseStruct) + if (.not. associated(ctx)) stop "dx_dt: context not set" + + ! put data in structure + call XTRY_2_STR(x_try, ctx%state1) + + ! run the fuse physics + call mod_derivs_diff(ctx) + + ! extract dx_dt from fuse structure + call STR_2_XTRY(ctx%dx_dt, g_x) + + ! track the total number of function calls + NUM_FUNCS = NUM_FUNCS + 1 + + end function dx_dt + + ! ----- calculate the Jacobian of g(x) ------------------------------------------------- + SUBROUTINE jac_flux(x,g_x,Jac) + IMPLICIT NONE + REAL(SP), DIMENSION(:), INTENT(IN) :: g_x + REAL(SP), DIMENSION(:), INTENT(INOUT) :: x + REAL(SP), DIMENSION(:,:), INTENT(OUT) :: Jac + REAL(SP), PARAMETER :: EPS=-1.0e-4_sp ! NOTE force h to be negative + INTEGER(I4B) :: j,n + REAL(SP), DIMENSION(size(x)) :: xsav,xph,h + xsav=x + n=size(x) + h=EPS*abs(xsav) + where (h == 0.0) h=EPS + xph=xsav+h + h=xph-xsav + do j=1,n + x(j)=xph(j) + Jac(:,j)=(dx_dt(x)-g_x(:))/h(j) + x(j)=xsav(j) + end do + NUM_JACOBIAN = NUM_JACOBIAN + 1 ! keep track of the number of iterations + call XTRY_2_STR(xsav, ctx%state1) ! restores consistency after finite differencing + end SUBROUTINE jac_flux + + ! ----- simple implicit solve for differentiable model -------------------------- + + subroutine implicit_solve(fuseStruct, x0, x1, nx) + USE nr, ONLY : lubksb,ludcmp + USE overshoot_module, only : get_bounds ! get state bounds + USE overshoot_module, only : fix_ovshoot ! fix overshoot (soft clamp) + USE model_numerix, only: ERR_ITER_FUNC ! Iteration convergence tolerance for function values + USE model_numerix, only: ERR_ITER_DX ! Iteration convergence tolerance for dx + implicit none + ! input-output + type(parent), intent(inout) :: fuseStruct ! parent fuse data structure + real(sp) , intent(in) :: x0(:) ! state vector at start of step + real(sp) , intent(out) :: x1(:) ! state vector at end of step + integer(i4b), intent(in) :: nx ! number of state variables + ! internal: newton iterations + real(sp) :: x_try(nx) ! trial state vector + real(sp) :: g_x(nx) ! dx/dt=g(x) + real(sp) :: res(nx) ! residual vector + real(sp) :: Jg(nx,nx) ! Jacobian matrix (flux) + real(sp) :: Jac(nx,nx) ! Jacobian matrix (full) + real(sp) :: dx(nx) ! state update + real(sp) :: phi ! half squared residual norm + real(sp) :: d ! determinant sign tracker + integer(i4b) :: indx(nx) ! LU pivot indices (row-swap bookkeeping) + integer(i4b) :: i ! index of state + integer(i4b) :: it ! index of newton iteration + integer(i4b), parameter :: maxit=100 ! maximum number of iterations + logical(lgt) :: converged ! flag for convergence + ! internal: backtracking line search w/ overshoot reject + real(sp) :: lambda ! backtrack length multiplier (lambda*dx) + real(sp) :: lower(nx) ! lower bound + real(sp) :: upper(nx) ! lower bound + real(sp) :: x_trial(nx) ! state vectorfor backtrack + real(sp) :: g_trial(nx) ! dx/dt=g(x) for backtrack + real(sp) :: res_trial(nx) ! residual for backtrack + real(sp) :: phi_new ! half squared residual norm + integer(i4b) :: ls_it ! index of line search iteration + logical(lgt) :: ovshoot ! flag for overshoot + logical(lgt) :: accepted ! flag for accepting newton step + ! line search params + real(sp), parameter :: shrink = 0.5_sp + real(sp), parameter :: c_armijo = 1e-4_sp + integer(i4b), parameter :: ls_max = 5 + + ! check dimension size + if (nx /= nState) stop "implicit_solve: nx /= nState" + + ! initialize number of calls + NUM_FUNCS = 0 ! number of function calls + NUM_JACOBIAN = 0 ! number of times Jacobian is calculated + + ! get the bounds for the state variables + ! NOTE: This can be done outside of the time and iteration loops (keeping here for now) + call get_bounds(fuseStruct, lower, upper) + + ! point to the fuse parent structure so that it is available in other routines + call set_dxdt_context(fuseStruct) + + ! put state vector into the fuse data structure + call XTRY_2_STR(x0, fuseStruct%state0) + + ! intialize state vector and convergence flag + x_try = x0 + accepted = .false. + converged = .false. + + ! --- F(x) and objective phi = 0.5*||F||^2 + g_x = dx_dt(x_try) + res = x_try - (x0 + g_x*dt) + phi = 0.5_sp * sum(res*res) + + ! iterate + do it = 1, maxit + + if (sqrt(2.0_sp*phi) < ERR_ITER_FUNC) then + converged = .true. + exit ! exit iteration loop + end if + + ! --- J(x) + call jac_flux(x_try, g_x, Jg) + Jac = -dt*Jg ! multiply dt + do i=1,nx; Jac(i,i) = Jac(i,i) + 1.0_sp; end do ! add identity matrix + + ! --- Solve J dx = -F (Newton step) + dx = -res + call ludcmp(Jac, indx, d) ! J overwritten with LU + call lubksb(Jac, indx, dx) ! dx becomes solution + + ! initialize flag to check if line search is accepted + accepted = .false. + + ! ---- backtracking line search w/ overshoot reject ---- + lambda = 1.0_sp + do ls_it = 1, ls_max + x_trial = x_try + lambda*dx + + ! check overshoot + ovshoot = any(x_trial < lower) .or. any(x_trial > upper) + if (.not. ovshoot) then + ! new function and residual + g_trial = dx_dt(x_trial) + res_trial = x_trial - (x0 + dt*g_trial) + phi_new = 0.5_sp * sum(res_trial*res_trial) + ! check for sufficient decrease (Armijo-lite) + if (phi_new <= (1.0_sp - c_armijo*lambda) * phi)then + accepted = .true. + exit + endif + end if + lambda = lambda * shrink + end do ! line search + + if (accepted) then + x_try = x_trial + g_x = g_trial + res = res_trial + phi = phi_new + else + ! ----- fallback: soft clamp a very small Newton step ----- + x_trial = x_try + lambda*dx + call fix_ovshoot(x_trial, lower, upper) + ! get new function evaluation + x_try = x_trial + g_x = dx_dt(x_try) + res = x_try - (x0 + g_x*dt) + phi = 0.5_sp * sum(res*res) + end if + + ! re-populate fuse data structure + call XTRY_2_STR(x_try, fuseStruct%state1) + + ! tiny-step convergence + if (maxval(abs(lambda*dx)) < ERR_ITER_DX) then + converged = .true. + exit ! exit iteration loop + end if + + end do ! loop through iterations + + ! save final state + x1 = x_try + + ! nullify pointer to the fuse structure + call clear_dxdt_context() + + ! check convergence + if( .not. converged) STOP "failed to converge in implicit_solve" + + end subroutine implicit_solve + +end module implicit_solve_module diff --git a/build/FUSE_SRC/physics/mod_derivs_diff.f90 b/build/FUSE_SRC/physics/mod_derivs_diff.f90 new file mode 100644 index 0000000..5dc1752 --- /dev/null +++ b/build/FUSE_SRC/physics/mod_derivs_diff.f90 @@ -0,0 +1,50 @@ +module MOD_DERIVS_DIFF_module + + USE nrtype + USE data_types, only: parent, statev + USE qsatexcess_diff_module, only: qsatexcess_diff + USE evap_upper_diff_module, only: evap_upper_diff + USE evap_lower_diff_module, only: evap_lower_diff + USE qinterflow_diff_module, only: qinterflow_diff + USE qpercolate_diff_module, only: qpercolate_diff + USE q_baseflow_diff_module, only: q_baseflow_diff + USE q_misscell_diff_module, only: q_misscell_diff + USE mstate_rhs_diff_module, only: mstate_rhs_diff + + implicit none + + private + public :: MOD_DERIVS_DIFF + +contains + + SUBROUTINE MOD_DERIVS_DIFF(fuseStruct) + ! --------------------------------------------------------------------------------------- + ! Creator: + ! -------- + ! Martyn Clark, 2007 + ! Modified to include snow model by Brian Henn, 6/2013 + ! Modified to include analytical derivatives by Martyn Clark, 12/2025 + ! --------------------------------------------------------------------------------------- + ! Purpose: + ! -------- + ! compute the time derivative (dx/dt) of all model states (x) + ! --------------------------------------------------------------------------------------- + implicit none + type(parent), intent(inout) :: fuseStruct ! parent fuse data structure + + ! compute fluxes + call qsatexcess_diff(fuseStruct) ! compute the saturated area and surface runoff + call evap_upper_diff(fuseStruct) ! compute evaporation from the upper layer + call evap_lower_diff(fuseStruct) ! compute evaporation from the lower layer + call qinterflow_diff(fuseStruct) ! compute interflow from free water in the upper layer + call qpercolate_diff(fuseStruct) ! compute percolation from the upper to lower soil layers + call q_baseflow_diff(fuseStruct) ! compute baseflow from the lower soil layer + call q_misscell_diff(fuseStruct) ! compute miscellaneous fluxes (NOTE: need sat area, evap, and perc) + + ! compute the time derivative (dx/dt) of all model states (x) + call mstate_rhs_diff(fuseStruct) + + END SUBROUTINE MOD_DERIVS_DIFF + +end module MOD_DERIVS_DIFF_module diff --git a/build/FUSE_SRC/physics/mstate_rhs_diff.f90 b/build/FUSE_SRC/physics/mstate_rhs_diff.f90 new file mode 100644 index 0000000..1ab0107 --- /dev/null +++ b/build/FUSE_SRC/physics/mstate_rhs_diff.f90 @@ -0,0 +1,77 @@ +module MSTATE_RHS_DIFF_module + + implicit none + + private + public :: MSTATE_RHS_DIFF + +contains + + SUBROUTINE MSTATE_RHS_DIFF(fuseStruct) + ! --------------------------------------------------------------------------------------- + ! Creator: + ! -------- + ! Martyn Clark, 2007 + ! Modified by Martyn Clark to create a differentiable model, 12/25 + ! --------------------------------------------------------------------------------------- + ! Purpose: + ! -------- + ! Computes time derivatives of all states for all model combinations + ! --------------------------------------------------------------------------------------- + USE nrtype ! variable types, etc. + USE data_types, only: parent ! fuse parent data type + USE model_defn ! model definition structure + USE model_defnames + ! input-output + type(parent), intent(inout) :: fuseStruct ! parent fuse data structure + ! ------------------------------------------------------------------------------------------------- + ! associate variables with elements of data structure + associate(& + M_FLUX => fuseStruct%flux , & ! fluxes + MPARAM => fuseStruct%param_adjust , & ! adjustable model parameters + DX_DT => fuseStruct%dx_dt & ! time derivative in states + ) ! (associate) + + ! --------------------------------------------------------------------------------------- + ! (1) COMPUTE TIME DERIVATIVES FOR STATES IN THE UPPER LAYER + ! --------------------------------------------------------------------------------------- + SELECT CASE(SMODL%iARCH1) + CASE(iopt_tension2_1) ! tension storage sub-divided into recharge and excess + DX_DT%TENS_1A = M_FLUX%EFF_PPT - M_FLUX%QSURF - M_FLUX%EVAP_1A - M_FLUX%RCHR2EXCS + DX_DT%TENS_1B = M_FLUX%RCHR2EXCS - M_FLUX%EVAP_1B - M_FLUX%TENS2FREE_1 + DX_DT%FREE_1 = M_FLUX%TENS2FREE_1 - M_FLUX%QPERC_12 - M_FLUX%QINTF_1 - M_FLUX%OFLOW_1 + CASE(iopt_tension1_1) ! upper layer broken up into tension and free storage + DX_DT%TENS_1 = M_FLUX%EFF_PPT - M_FLUX%QSURF - M_FLUX%EVAP_1 - M_FLUX%TENS2FREE_1 + DX_DT%FREE_1 = M_FLUX%TENS2FREE_1 - M_FLUX%QPERC_12 - M_FLUX%QINTF_1 - M_FLUX%OFLOW_1 + CASE(iopt_onestate_1) ! upper layer defined by a single state variable + DX_DT%WATR_1 = M_FLUX%EFF_PPT - M_FLUX%QSURF - M_FLUX%EVAP_1 - M_FLUX%QPERC_12 - M_FLUX%QINTF_1 & + - M_FLUX%OFLOW_1 + CASE DEFAULT + print *, "SMODL%iARCH1 must be iopt_tension2_1, iopt_tension1_1, or iopt_onestate_1" + STOP + END SELECT ! (upper layer architechure) + + ! --------------------------------------------------------------------------------------- + ! (2) COMPUTE TIME DERIVATIVES FOR STATES IN THE LOWER LAYER + ! --------------------------------------------------------------------------------------- + SELECT CASE(SMODL%iARCH2) + CASE(iopt_tens2pll_2) ! tension reservoir plus two parallel tanks + DX_DT%TENS_2 = M_FLUX%QPERC_12*(1._SP-MPARAM%PERCFRAC) - M_FLUX%EVAP_2 - M_FLUX%TENS2FREE_2 + DX_DT%FREE_2A = M_FLUX%QPERC_12*(MPARAM%PERCFRAC/2._SP) + M_FLUX%TENS2FREE_2/2._SP - M_FLUX%QBASE_2A & + - M_FLUX%OFLOW_2A + DX_DT%FREE_2B = M_FLUX%QPERC_12*(MPARAM%PERCFRAC/2._SP) + M_FLUX%TENS2FREE_2/2._SP - M_FLUX%QBASE_2B & + - M_FLUX%OFLOW_2B + CASE(iopt_unlimfrc_2,iopt_unlimpow_2,iopt_topmdexp_2,iopt_fixedsiz_2) ! single state + ! (NOTE: M_FLUX%OFLOW_2=0 for 'unlimfrc_2','unlimpow_2','topmdexp_2') + DX_DT%WATR_2 = M_FLUX%QPERC_12 - M_FLUX%EVAP_2 - M_FLUX%QBASE_2 - M_FLUX%OFLOW_2 + CASE DEFAULT + print *, "SMODL%iARCH2 must be iopt_tens2pll_2, iopt_unlimfrc_2, iopt_unlimpow_2" + print *, " iopt_topmdexp_2, or iopt_fixedsiz_2" + STOP + END SELECT + ! --------------------------------------------------------------------------------------- + + end associate ! end association with variables in the data structures + END SUBROUTINE MSTATE_RHS_DIFF + +end module MSTATE_RHS_DIFF_module diff --git a/build/FUSE_SRC/physics/q_baseflow_diff.f90 b/build/FUSE_SRC/physics/q_baseflow_diff.f90 new file mode 100644 index 0000000..29386f4 --- /dev/null +++ b/build/FUSE_SRC/physics/q_baseflow_diff.f90 @@ -0,0 +1,69 @@ +module Q_BASEFLOW_DIFF_module + + implicit none + + private + public :: Q_BASEFLOW_DIFF + +contains + + + SUBROUTINE Q_BASEFLOW_DIFF(fuseStruct) + ! --------------------------------------------------------------------------------------- + ! Creator: + ! -------- + ! Martyn Clark, 2007 + ! Modified by Martyn Clark to create a differentiable model, 12/25 + ! --------------------------------------------------------------------------------------- + ! Purpose: + ! -------- + ! Computes the baseflow from the lower soil layer + ! --------------------------------------------------------------------------------------- + USE nrtype ! variable types, etc. + USE data_types, only: parent ! fuse parent data type + USE model_defn ! model definition structure + USE model_defnames + IMPLICIT NONE + ! input-output + type(parent), intent(inout) :: fuseStruct ! parent fuse data structure + ! ------------------------------------------------------------------------------------------------- + ! associate variables with elements of data structure + associate(& + M_FLUX => fuseStruct%flux , & ! fluxes + TSTATE => fuseStruct%state1 , & ! trial state variables (end of step) + MPARAM => fuseStruct%param_adjust , & ! adjustable model parameters + DPARAM => fuseStruct%param_derive & ! derived model parameters + ) ! (associate) + + ! --------------------------------------------------------------------------------------- + SELECT CASE(SMODL%iARCH2) + ! -------------------------------------------------------------------------------------- + CASE(iopt_tens2pll_2) ! tension reservoir plus two parallel tanks + M_FLUX%QBASE_2A = MPARAM%QBRATE_2A * TSTATE%FREE_2A ! qbrate_2a is a fraction (T-1) + M_FLUX%QBASE_2B = MPARAM%QBRATE_2B * TSTATE%FREE_2B ! qbrate_2b is a fraction (T-1) + M_FLUX%QBASE_2 = M_FLUX%QBASE_2A + M_FLUX%QBASE_2B ! total baseflow + ! -------------------------------------------------------------------------------------- + CASE(iopt_unlimfrc_2) ! baseflow resvr of unlimited size (0-HUGE), frac rate + M_FLUX%QBASE_2 = MPARAM%QB_PRMS * TSTATE%WATR_2 ! qb_prms is a fraction (T-1) + ! -------------------------------------------------------------------------------------- + CASE(iopt_unlimpow_2) ! baseflow resvr of unlimited size (0-HUGE), power recession + M_FLUX%QBASE_2 = DPARAM%QBSAT * (TSTATE%WATR_2/MPARAM%MAXWATR_2)**MPARAM%QB_POWR + ! -------------------------------------------------------------------------------------- + CASE(iopt_topmdexp_2) ! topmodel exponential reservoir (-HUGE to HUGE) + M_FLUX%QBASE_2 = DPARAM%QBSAT * EXP( -(1. - TSTATE%WATR_2/MPARAM%MAXWATR_2) ) + ! -------------------------------------------------------------------------------------- + CASE(iopt_fixedsiz_2) ! baseflow reservoir of fixed size + M_FLUX%QBASE_2 = MPARAM%BASERTE * (TSTATE%WATR_2/MPARAM%MAXWATR_2)**MPARAM%QB_POWR + ! -------------------------------------------------------------------------------------- + CASE DEFAULT + print *, "SMODL%iARCH2 must be iopt_tens2pll_2, iopt_unlimfrc_2, iopt_unlimpow_2" + print *, " iopt_topmdexp_2, or iopt_fixedsiz_2" + STOP + ! -------------------------------------------------------------------------------------- + END SELECT + ! --------------------------------------------------------------------------------------- + + end associate ! end association with variables in the data structures + END SUBROUTINE Q_BASEFLOW_DIFF + +end module Q_BASEFLOW_DIFF_module diff --git a/build/FUSE_SRC/physics/q_misscell_diff.f90 b/build/FUSE_SRC/physics/q_misscell_diff.f90 new file mode 100644 index 0000000..1801c7b --- /dev/null +++ b/build/FUSE_SRC/physics/q_misscell_diff.f90 @@ -0,0 +1,116 @@ +module Q_MISSCELL_DIFF_module + + implicit none + + private + public :: Q_MISSCELL_DIFF + +contains + + SUBROUTINE Q_MISSCELL_DIFF(fuseStruct) + ! --------------------------------------------------------------------------------------- + ! Creator: + ! -------- + ! Martyn Clark, 2007 + ! Modified by Martyn Clark to create a differentiable model, 12/25 + ! --------------------------------------------------------------------------------------- + ! Purpose: + ! -------- + ! Computes miscellaneous fluxes: + ! RCHR2EXCS = flow from recharge to excess (mm day-1) + ! TENS2FREE_1 = flow from tension storage to free storage in the upper layer (mm day-1) + ! TENS2FREE_2 = flow from tension storage to free storage in the lower layer (mm day-1) + ! OFLOW_1 = overflow from the upper soil layer (mm day-1) + ! OFLOW_2 = overflow from the lower soil layer (mm day-1) + ! --------------------------------------------------------------------------------------- + USE nrtype ! variable types, etc. + USE data_types, only: parent ! fuse parent data type + USE model_defn ! model definition structure + USE model_defnames + IMPLICIT NONE + ! input-output + type(parent), intent(inout) :: fuseStruct ! parent fuse data structure + ! internal + REAL(SP) :: LOGISMOOTH ! FUNCTION logistic smoothing + REAL(SP), PARAMETER :: PSMOOTH=0.01_SP ! smoothing parameter + REAL(SP) :: W_FUNC ! result from LOGISMOOTH + ! ------------------------------------------------------------------------------------------------- + ! associate variables with elements of data structure + associate(& + M_FLUX => fuseStruct%flux , & ! fluxes + TSTATE => fuseStruct%state1 , & ! trial state variables (end of step) + MPARAM => fuseStruct%param_adjust , & ! adjustable model parameters + DPARAM => fuseStruct%param_derive & ! derived model parameters + ) ! (associate) + ! --------------------------------------------------------------------------------------- + + ! --------------------------------------------------------------------------------------- + SELECT CASE(SMODL%iARCH1) + CASE(iopt_tension2_1) ! tension storage sub-divided into recharge and excess + ! compute flow from recharge to excess (mm s-1) + W_FUNC = LOGISMOOTH(TSTATE%TENS_1A,DPARAM%MAXTENS_1A,PSMOOTH) + M_FLUX%RCHR2EXCS = W_FUNC * (M_FLUX%EFF_PPT - M_FLUX%QSURF) + ! compute flow from tension storage to free storage (mm s-1) + W_FUNC = LOGISMOOTH(TSTATE%TENS_1B,DPARAM%MAXTENS_1B,PSMOOTH) + M_FLUX%TENS2FREE_1 = W_FUNC * M_FLUX%RCHR2EXCS + ! compute over-flow of free water + W_FUNC = LOGISMOOTH(TSTATE%FREE_1,DPARAM%MAXFREE_1,PSMOOTH) + M_FLUX%OFLOW_1 = W_FUNC * M_FLUX%TENS2FREE_1 + CASE(iopt_tension1_1) ! upper layer broken up into tension and free storage + ! no separate recharge zone (flux should never be used) + M_FLUX%RCHR2EXCS = 0._SP + ! compute flow from tension storage to free storage (mm s-1) + W_FUNC = LOGISMOOTH(TSTATE%TENS_1,DPARAM%MAXTENS_1,PSMOOTH) + M_FLUX%TENS2FREE_1 = W_FUNC * (M_FLUX%EFF_PPT - M_FLUX%QSURF) + ! compute over-flow of free water + W_FUNC = LOGISMOOTH(TSTATE%FREE_1,DPARAM%MAXFREE_1,PSMOOTH) + M_FLUX%OFLOW_1 = W_FUNC * M_FLUX%TENS2FREE_1 + CASE(iopt_onestate_1) ! upper layer defined by a single state variable + ! no tension stores + M_FLUX%RCHR2EXCS = 0._SP + M_FLUX%TENS2FREE_1 = 0._SP + ! compute over-flow of free water + W_FUNC = LOGISMOOTH(TSTATE%WATR_1,MPARAM%MAXWATR_1,PSMOOTH) + M_FLUX%OFLOW_1 = W_FUNC * (M_FLUX%EFF_PPT - M_FLUX%QSURF) + CASE DEFAULT + print *, "SMODL%iARCH1 must be iopt_tension2_1, iopt_tension1_1, or iopt_onestate_1" + STOP + END SELECT + + ! --------------------------------------------------------------------------------------- + SELECT CASE(SMODL%iARCH2) + CASE(iopt_tens2pll_2) ! tension reservoir plus two parallel tanks + ! compute flow from tension storage to free storage (mm s-1) + W_FUNC = LOGISMOOTH(TSTATE%TENS_2,DPARAM%MAXTENS_2,PSMOOTH) + M_FLUX%TENS2FREE_2 = W_FUNC * M_FLUX%QPERC_12*(1._SP-MPARAM%PERCFRAC) + ! compute over-flow of free water in the primary reservoir + W_FUNC = LOGISMOOTH(TSTATE%FREE_2A,DPARAM%MAXFREE_2A,PSMOOTH) + M_FLUX%OFLOW_2A = W_FUNC * (M_FLUX%QPERC_12*(MPARAM%PERCFRAC/2._SP) + M_FLUX%TENS2FREE_2/2._SP) + ! compute over-flow of free water in the secondary reservoir + W_FUNC = LOGISMOOTH(TSTATE%FREE_2B,DPARAM%MAXFREE_2B,PSMOOTH) + M_FLUX%OFLOW_2B = W_FUNC * (M_FLUX%QPERC_12*(MPARAM%PERCFRAC/2._SP) + M_FLUX%TENS2FREE_2/2._SP) + ! compute total overflow + M_FLUX%OFLOW_2 = M_FLUX%OFLOW_2A + M_FLUX%OFLOW_2B + CASE(iopt_fixedsiz_2) + ! no tension store + M_FLUX%TENS2FREE_2 = 0._SP + M_FLUX%OFLOW_2A = 0._SP + M_FLUX%OFLOW_2B = 0._SP + ! compute over-flow of free water + W_FUNC = LOGISMOOTH(TSTATE%WATR_2,MPARAM%MAXWATR_2,PSMOOTH) + M_FLUX%OFLOW_2 = W_FUNC * M_FLUX%QPERC_12 + CASE(iopt_unlimfrc_2,iopt_unlimpow_2,iopt_topmdexp_2) ! unlimited size + M_FLUX%TENS2FREE_2 = 0._SP + M_FLUX%OFLOW_2 = 0._SP + M_FLUX%OFLOW_2A = 0._SP + M_FLUX%OFLOW_2B = 0._SP + CASE DEFAULT + print *, "SMODL%iARCH2 must be iopt_tens2pll_2, iopt_unlimfrc_2, iopt_unlimpow_2" + print *, " iopt_topmdexp_2, or iopt_fixedsiz_2" + STOP + END SELECT + + end associate ! end association with variables in the data structures + END SUBROUTINE Q_MISSCELL_DIFF + +end module Q_MISSCELL_DIFF_module diff --git a/build/FUSE_SRC/physics/qinterflow_diff.f90 b/build/FUSE_SRC/physics/qinterflow_diff.f90 new file mode 100644 index 0000000..9b1ed32 --- /dev/null +++ b/build/FUSE_SRC/physics/qinterflow_diff.f90 @@ -0,0 +1,52 @@ +module QINTERFLOW_DIFF_module + + implicit none + + private + public :: QINTERFLOW_DIFF + +contains + + SUBROUTINE QINTERFLOW_DIFF(fuseStruct) + ! --------------------------------------------------------------------------------------- + ! Creator: + ! -------- + ! Martyn Clark, 2007 + ! Modified by Martyn Clark to create a differentiable model, 12/25 + ! --------------------------------------------------------------------------------------- + ! Purpose: + ! -------- + ! Computes the interflow from free water in the upper soil layer + ! --------------------------------------------------------------------------------------- + USE nrtype ! variable types, etc. + USE data_types, only: parent ! fuse parent data type + USE model_defn ! model definition structure + USE model_defnames + IMPLICIT NONE + ! input-output + type(parent), intent(inout) :: fuseStruct ! parent fuse data structure + ! ------------------------------------------------------------------------------------------------- + ! associate variables with elements of data structure + associate(& + M_FLUX => fuseStruct%flux , & ! fluxes + TSTATE => fuseStruct%state1 , & ! trial state variables (end of step) + MPARAM => fuseStruct%param_adjust , & ! adjustable model parameters + DPARAM => fuseStruct%param_derive & ! derived model parameters + ) ! (associate) + + ! --------------------------------------------------------------------------------------- + SELECT CASE(SMODL%iQINTF) + CASE(iopt_intflwsome) ! interflow + M_FLUX%QINTF_1 = MPARAM%IFLWRTE * (TSTATE%FREE_1/DPARAM%MAXFREE_1) + CASE(iopt_intflwnone) ! no interflow + M_FLUX%QINTF_1 = 0. + CASE DEFAULT ! check for errors + print *, "SMODL%iQINTF must be either iopt_intflwsome or iopt_intflwnone" + STOP + END SELECT + ! --------------------------------------------------------------------------------------- + + end associate ! end association with variables in the data structures + END SUBROUTINE QINTERFLOW_DIFF + +end module QINTERFLOW_DIFF_module diff --git a/build/FUSE_SRC/physics/qpercolate_diff.f90 b/build/FUSE_SRC/physics/qpercolate_diff.f90 new file mode 100644 index 0000000..9ff599c --- /dev/null +++ b/build/FUSE_SRC/physics/qpercolate_diff.f90 @@ -0,0 +1,61 @@ +module QPERCOLATE_DIFF_module + + implicit none + + private + public :: QPERCOLATE_DIFF + +contains + + SUBROUTINE QPERCOLATE_DIFF(fuseStruct) + ! --------------------------------------------------------------------------------------- + ! Creator: + ! -------- + ! Martyn Clark, 2007 + ! Modified by Martyn Clark to create a differentiable model, 12/25 + ! --------------------------------------------------------------------------------------- + ! Purpose: + ! -------- + ! Computes the percolation from the upper soil layer to the lower soil layer + ! --------------------------------------------------------------------------------------- + USE nrtype ! variable types, etc. + USE data_types, only: parent ! fuse parent data type + USE model_defn ! model definition structure + USE model_defnames + IMPLICIT NONE + ! input-output + type(parent), intent(inout) :: fuseStruct ! parent fuse data structure + ! internal + REAL(SP) :: LZ_PD ! lower zone percolation demand + ! --------------------------------------------------------------------------------------- + ! associate variables with elements of data structure + associate(& + M_FLUX => fuseStruct%flux , & ! fluxes + TSTATE => fuseStruct%state1 , & ! trial state variables (end of step) + MPARAM => fuseStruct%param_adjust , & ! adjustable model parameters + DPARAM => fuseStruct%param_derive & ! derived model parameters + ) ! (associate) + ! --------------------------------------------------------------------------------------- + + ! --------------------------------------------------------------------------------------- + SELECT CASE(SMODL%iQPERC) + CASE(iopt_perc_f2sat) ! water from (field cap to sat) avail for percolation + M_FLUX%QPERC_12 = MPARAM%PERCRTE * (TSTATE%FREE_1/DPARAM%MAXFREE_1)**MPARAM%PERCEXP + CASE(iopt_perc_w2sat) ! water from (wilt pt to sat) avail for percolation + M_FLUX%QPERC_12 = MPARAM%PERCRTE * (TSTATE%WATR_1/MPARAM%MAXWATR_1)**MPARAM%PERCEXP + CASE(iopt_perc_lower) ! perc defined by moisture content in lower layer (SAC) + ! (compute lower-zone percolation demand -- multiplier on maximum percolation, then percolation) + LZ_PD = 1._SP + MPARAM%SACPMLT*(1._SP - TSTATE%WATR_2/MPARAM%MAXWATR_2)**MPARAM%SACPEXP + M_FLUX%QPERC_12 = DPARAM%QBSAT*LZ_PD * (TSTATE%FREE_1/DPARAM%MAXFREE_1) + !print *, 'lz_pd = ', LZ_PD, MPARAM%SACPMLT, TSTATE%WATR_2/MPARAM%MAXWATR_2, MPARAM%SACPEXP + !print *, 'qperc_12 = ', M_FLUX%QPERC_12, DPARAM%QBSAT, LZ_PD, TSTATE%FREE_1/DPARAM%MAXFREE_1 + CASE DEFAULT ! check for errors + print *, "SMODL%iQPERC must be iopt_perc_f2sat, iopt_perc_w2sat, or iopt_perc_lower" + STOP + END SELECT + ! -------------------------------------------------------------------------------------- + + end associate ! end association with variables in the data structures- + END SUBROUTINE QPERCOLATE_DIFF + +end module QPERCOLATE_DIFF_module diff --git a/build/FUSE_SRC/physics/qsatexcess_diff.f90 b/build/FUSE_SRC/physics/qsatexcess_diff.f90 new file mode 100644 index 0000000..fa454a2 --- /dev/null +++ b/build/FUSE_SRC/physics/qsatexcess_diff.f90 @@ -0,0 +1,106 @@ +module QSATEXCESS_DIFF_MODULE + + implicit none + + private + public :: QSATEXCESS_DIFF + +contains + + SUBROUTINE QSATEXCESS_DIFF(fuseStruct) + ! ------------------------------------------------------------------------------------------------- + ! Creator: + ! -------- + ! Martyn Clark, 2007 + ! Modified by Martyn Clark to create a differentiable model, 12/25 + ! ------------------------------------------------------------------------------------------------- + ! Purpose: + ! -------- + ! Computes the saturated area and surface runoff + ! ------------------------------------------------------------------------------------------------- + USE nrtype ! variable types, etc. + USE data_types, only: parent ! fuse parent data type + USE nr, ONLY : gammp ! interface for the incomplete gamma function + USE model_defn ! model definition structure + USE model_defnames + IMPLICIT NONE + ! input-output + type(parent), intent(inout) :: fuseStruct ! parent fuse data structure + ! internal variables + REAL(SP) :: TI_SAT ! topographic index where saturated + REAL(SP) :: TI_LOG ! critical value of topo index in log space + REAL(SP) :: TI_OFF ! offset in the Gamma distribution + REAL(SP) :: TI_SHP ! shape of the Gamma distribution + REAL(SP) :: TI_CHI ! CHI, see Sivapalan et al., 1987 + REAL(SP) :: TI_ARG ! argument of the Gamma function + REAL(SP) :: NO_ZERO=1.E-8 ! avoid divide by zero + ! ------------------------------------------------------------------------------------------------- + ! associate variables with elements of data structure + associate(& + M_FLUX => fuseStruct%flux , & ! fluxes + TSTATE => fuseStruct%state1 , & ! trial state variables (end of step) + MPARAM => fuseStruct%param_adjust , & ! adjustable model parameters + DPARAM => fuseStruct%param_derive & ! derived model parameters + ) ! (associate) + ! ------------------------------------------------------------------------------------------------- + + ! saturated area method + SELECT CASE(SMODL%iQSURF) + + ! ------------------------------------------------------------------------------------------------ + ! ----- ARNO/Xzang/VIC parameterization (upper zone control) ------------------------------------- + ! ------------------------------------------------------------------------------------------------ + CASE(iopt_arno_x_vic) + + ! ----- compute flux ---------------------------------------------------------------------------- + M_FLUX%SATAREA = 1._sp - ( 1._sp - MIN(TSTATE%WATR_1/MPARAM%MAXWATR_1, 1._sp) )**MPARAM%AXV_BEXP + + ! ------------------------------------------------------------------------------------------------ + ! ----- PRMS variant (fraction of upper tension storage) ----------------------------------------- + ! ------------------------------------------------------------------------------------------------ + CASE(iopt_prms_varnt) + + ! ----- compute flux ---------------------------------------------------------------------------- + M_FLUX%SATAREA = MIN(TSTATE%TENS_1/DPARAM%MAXTENS_1, 1._sp) * MPARAM%SAREAMAX + + ! ------------------------------------------------------------------------------------------------ + ! ----- TOPMODEL parameterization (only valid for TOPMODEL qb) ----------------------------------- + ! ------------------------------------------------------------------------------------------------ + CASE(iopt_tmdl_param) + + ! ----- compute flux ---------------------------------------------------------------------------- + + ! compute the minimum value of the topographic index where the basin is saturated + ! (this is correct, as MPARAM%MAXWATR_2 is m*n -- units are meters**(1/n) + TI_SAT = DPARAM%POWLAMB / (TSTATE%WATR_2/MPARAM%MAXWATR_2 + NO_ZERO) + ! compute the saturated area + IF (TI_SAT.GT.DPARAM%MAXPOW) THEN + M_FLUX%SATAREA = 0. + ELSE + ! convert the topographic index to log space + TI_LOG = LOG( TI_SAT**MPARAM%QB_POWR ) + ! compute the saturated area (NOTE: critical value of the topographic index is in log space) + TI_OFF = 3._sp ! offset in the Gamma distribution (the "3rd" parameter) + TI_SHP = MPARAM%TISHAPE ! shape of the Gamma distribution (the "2nd" parameter) + TI_CHI = (MPARAM%LOGLAMB - TI_OFF) / MPARAM%TISHAPE ! Chi -- loglamb is the first parameter (mean) + TI_ARG = MAX(0._sp, TI_LOG - TI_OFF) / TI_CHI ! argument to the incomplete Gamma function + M_FLUX%SATAREA = 1._sp - GAMMP(TI_SHP, TI_ARG) ! GAMMP is the incomplete Gamma function + ENDIF + + ! ------------------------------------------------------------------------------------------------ + ! ------------------------------------------------------------------------------------------------ + ! check processed surface runoff selection + CASE DEFAULT + print *, "SMODL%iQSURF must be iopt_arno_x_vic, iopt_prms_varnt, or iopt_tmdl_param" + STOP + + END SELECT ! (different surface runoff options) + + ! ...and, compute surface runoff + ! ------------------------------ + M_FLUX%QSURF = M_FLUX%EFF_PPT * M_FLUX%SATAREA + + end associate ! end association with variables in the data structures + END SUBROUTINE QSATEXCESS_DIFF + +end module QSATEXCESS_DIFF_MODULE diff --git a/build/Makefile b/build/Makefile index 4019dfc..9c2edfb 100644 --- a/build/Makefile +++ b/build/Makefile @@ -7,7 +7,7 @@ #======================================================================== # Define core directory below which everything resides -F_MASTER = ${HOME}/fuse/ +F_MASTER = ${HOME}/Documents/analysis/diffModel/FUSE/source/fuse/ # Core directory that contains FUSE source code F_KORE_DIR = $(F_MASTER)build/FUSE_SRC/ @@ -22,14 +22,9 @@ EXE_PATH = $(F_MASTER)bin/ # PART 1: Define the libraries, driver programs, and executables #======================================================================== -# Define the fortran compiler. You have a few options: i) set the FC -# variable in your environment, ii) set it when you compile this -# Make file (e.g. make FC=ifort), iii) or if don't define it, the compiler -# specified below is used -ifndef FC - #FC = ifort - FC = gfortran -endif +# Define the fortran compiler. +#FC = ifort +FC = gfortran # Define the NetCDF and HDF5 libraries. Use the libraries associated with # the compiler you selected above. Note that these paths are machine-dependent @@ -41,13 +36,12 @@ ifeq "$(FC)" "ifort" endif ifeq "$(FC)" "gfortran" - NCDF_LIB_PATH = /usr/lib/x86_64-linux-gnu#${NCDF_PATH} - HDF_LIB_PATH = /usr/lib/x86_64-linux-gnu/hdf5/serial#${HDF_PATH} - INCLUDE_PATH = /usr#${IN_PATH} + INC_NETCDF := $(shell nf-config --fflags) + LIB_NETCDF := $(shell nf-config --flibs) $(shell nc-config --libs) endif -LIBRARIES = -L$(NCDF_LIB_PATH)/lib -lnetcdff -lnetcdf -L$(HDF_LIB_PATH)/lib -lhdf5_hl -lhdf5 -INCLUDE = -I$(INCLUDE_PATH)/include -I$(INCLUDE_PATH)/include +LIBRARIES = $(LIB_NETCDF) +INCLUDE = $(INC_NETCDF) # Define the driver program and associated subroutines for the fidelity test FUSE_DRIVER = \ @@ -65,13 +59,15 @@ DRIVER_EX = fuse.exe #======================================================================== # Define directories -NUMREC_DIR = $(F_KORE_DIR)FUSE_NR -HOOKUP_DIR = $(F_KORE_DIR)FUSE_HOOK -DRIVER_DIR = $(F_KORE_DIR)FUSE_DMSL -NETCDF_DIR = $(F_KORE_DIR)FUSE_NETCDF -ENGINE_DIR = $(F_KORE_DIR)FUSE_ENGINE -SCE_DIR = $(F_KORE_DIR)FUSE_SCE -TIME_DIR = $(F_KORE_DIR)FUSE_TIME +NUMREC_DIR = $(F_KORE_DIR)FUSE_NR +HOOKUP_DIR = $(F_KORE_DIR)FUSE_HOOK +DRIVER_DIR = $(F_KORE_DIR)FUSE_DMSL +NETCDF_DIR = $(F_KORE_DIR)FUSE_NETCDF +DSHARE_DIR = $(F_KORE_DIR)dshare +PHYSICS_DIR = $(F_KORE_DIR)physics +ENGINE_DIR = $(F_KORE_DIR)FUSE_ENGINE +SCE_DIR = $(F_KORE_DIR)FUSE_SCE +TIME_DIR = $(F_KORE_DIR)FUSE_TIME # Utility modules FUSE_UTILMS= \ @@ -89,7 +85,9 @@ NRUTIL = $(patsubst %, $(NUMREC_DIR)/%, $(FUSE_NRUTIL)) # Data modules FUSE_DATAMS= \ model_defn.f90 \ + data_types.f90 \ model_defnames.f90 \ + globaldata.f90 \ multiconst.f90 \ multiforce.f90 \ multibands.f90 \ @@ -99,7 +97,7 @@ FUSE_DATAMS= \ multiroute.f90 \ multistats.f90 \ model_numerix.f90 -DATAMS = $(patsubst %, $(ENGINE_DIR)/%, $(FUSE_DATAMS)) +DATAMS = $(patsubst %, $(DSHARE_DIR)/%, $(FUSE_DATAMS)) # Time I/O modules FUSE_TIMEMS= \ @@ -128,6 +126,22 @@ FUSE_NR_SUB= \ gammln.f90 gammp.f90 gcf.f90 gser.f90 NR_SUB = $(patsubst %, $(NUMREC_DIR)/%, $(FUSE_NR_SUB)) +# FUSE physics +FUSE_PHYSICS= \ + get_parent.f90 \ + qsatexcess_diff.f90 \ + evap_upper_diff.f90 \ + evap_lower_diff.f90 \ + qinterflow_diff.f90 \ + qpercolate_diff.f90 \ + q_baseflow_diff.f90 \ + q_misscell_diff.f90 \ + mstate_rhs_diff.f90 \ + mod_derivs_diff.f90 \ + fix_ovshoot.f90 \ + implicit_solve.f90 +PHYSICS = $(patsubst %, $(PHYSICS_DIR)/%, $(FUSE_PHYSICS)) + # Model guts FUSE_MODGUT=\ mod_derivs.f90 \ @@ -214,7 +228,7 @@ SCE = \ # ... and stitch it all together... FUSE_ALL = $(UTILMS) $(NRUTIL) $(DATAMS) $(TIMUTILS) $(INFOMS) \ - $(NR_SUB) $(MODGUT) $(SOLVER) $(PRELIM) $(MODRUN) \ + $(NR_SUB) $(PHYSICS) $(MODGUT) $(SOLVER) $(PRELIM) $(MODRUN) \ $(NETCDF) $(SCE) #========================================================================