From 483ef4557341501697f2d1d48f42484501c01db6 Mon Sep 17 00:00:00 2001 From: Martyn Clark Date: Wed, 26 Nov 2025 08:40:07 -0700 Subject: [PATCH 1/2] minor changes to get working on a mac --- build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 | 2 +- build/FUSE_SRC/FUSE_NETCDF/def_output.f90 | 14 +++++++------- build/FUSE_SRC/FUSE_NETCDF/put_output.f90 | 4 ++-- build/Makefile | 22 ++++++++-------------- 4 files changed, 18 insertions(+), 24 deletions(-) diff --git a/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 b/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 index b64355c..bca594e 100644 --- a/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 +++ b/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 @@ -247,7 +247,7 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG ! 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 + IF (IERR.NE.0) THEN; PRINT *, TRIM(MESSAGE); STOP 1; ENDIF ! perform overland flow routing CALL Q_OVERLAND() diff --git a/build/FUSE_SRC/FUSE_NETCDF/def_output.f90 b/build/FUSE_SRC/FUSE_NETCDF/def_output.f90 index 6e71c93..f04bbb5 100644 --- a/build/FUSE_SRC/FUSE_NETCDF/def_output.f90 +++ b/build/FUSE_SRC/FUSE_NETCDF/def_output.f90 @@ -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/put_output.f90 b/build/FUSE_SRC/FUSE_NETCDF/put_output.f90 index bc2e361..8d1b13e 100644 --- a/build/FUSE_SRC/FUSE_NETCDF/put_output.f90 +++ b/build/FUSE_SRC/FUSE_NETCDF/put_output.f90 @@ -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) @@ -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/Makefile b/build/Makefile index 4019dfc..01cb576 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 = \ From f9e3fb4d06886ea0cd1daff143b7a050e991e120 Mon Sep 17 00:00:00 2001 From: Martyn Clark Date: Sat, 29 Nov 2025 08:32:08 -0700 Subject: [PATCH 2/2] new scaffolding for differentiable model --- build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 | 3 +- build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 | 44 +++- build/FUSE_SRC/FUSE_ENGINE/assign_par.f90 | 3 +- build/FUSE_SRC/FUSE_ENGINE/getnumerix.f90 | 10 +- build/FUSE_SRC/FUSE_ENGINE/getpar_str.f90 | 3 +- build/FUSE_SRC/FUSE_ENGINE/getparmeta.f90 | 2 +- build/FUSE_SRC/FUSE_ENGINE/multi_flux.f90 | 42 --- build/FUSE_SRC/FUSE_ENGINE/multistate.f90 | 53 ---- build/FUSE_SRC/FUSE_ENGINE/multistats.f90 | 35 --- build/FUSE_SRC/FUSE_ENGINE/putpar_str.f90 | 3 +- build/FUSE_SRC/FUSE_ENGINE/str_2_xtry.f90 | 2 +- .../multiparam.f90 => dshare/data_types.f90} | 201 +++++++++++++-- .../{FUSE_ENGINE => dshare}/model_defn.f90 | 0 .../model_defnames.f90 | 0 .../{FUSE_ENGINE => dshare}/model_numerix.f90 | 3 + build/FUSE_SRC/dshare/multi_flux.f90 | 11 + .../{FUSE_ENGINE => dshare}/multibands.f90 | 0 .../{FUSE_ENGINE => dshare}/multiconst.f90 | 0 .../{FUSE_ENGINE => dshare}/multiforce.f90 | 31 +-- build/FUSE_SRC/dshare/multiparam.f90 | 22 ++ .../{FUSE_ENGINE => dshare}/multiroute.f90 | 0 build/FUSE_SRC/dshare/multistate.f90 | 27 ++ build/FUSE_SRC/dshare/multistats.f90 | 9 + build/FUSE_SRC/physics/evap_lower_diff.f90 | 88 +++++++ build/FUSE_SRC/physics/evap_upper_diff.f90 | 91 +++++++ build/FUSE_SRC/physics/fix_ovshoot.f90 | 139 ++++++++++ build/FUSE_SRC/physics/get_parent.f90 | 26 ++ build/FUSE_SRC/physics/implicit_solve.f90 | 244 ++++++++++++++++++ build/FUSE_SRC/physics/mod_derivs_diff.f90 | 50 ++++ build/FUSE_SRC/physics/mstate_rhs_diff.f90 | 77 ++++++ build/FUSE_SRC/physics/q_baseflow_diff.f90 | 69 +++++ build/FUSE_SRC/physics/q_misscell_diff.f90 | 116 +++++++++ build/FUSE_SRC/physics/qinterflow_diff.f90 | 52 ++++ build/FUSE_SRC/physics/qpercolate_diff.f90 | 61 +++++ build/FUSE_SRC/physics/qsatexcess_diff.f90 | 106 ++++++++ build/Makefile | 37 ++- 36 files changed, 1455 insertions(+), 205 deletions(-) delete mode 100644 build/FUSE_SRC/FUSE_ENGINE/multi_flux.f90 delete mode 100644 build/FUSE_SRC/FUSE_ENGINE/multistate.f90 delete mode 100644 build/FUSE_SRC/FUSE_ENGINE/multistats.f90 rename build/FUSE_SRC/{FUSE_ENGINE/multiparam.f90 => dshare/data_types.f90} (52%) rename build/FUSE_SRC/{FUSE_ENGINE => dshare}/model_defn.f90 (100%) rename build/FUSE_SRC/{FUSE_ENGINE => dshare}/model_defnames.f90 (100%) rename build/FUSE_SRC/{FUSE_ENGINE => dshare}/model_numerix.f90 (96%) create mode 100644 build/FUSE_SRC/dshare/multi_flux.f90 rename build/FUSE_SRC/{FUSE_ENGINE => dshare}/multibands.f90 (100%) rename build/FUSE_SRC/{FUSE_ENGINE => dshare}/multiconst.f90 (100%) rename build/FUSE_SRC/{FUSE_ENGINE => dshare}/multiforce.f90 (84%) create mode 100644 build/FUSE_SRC/dshare/multiparam.f90 rename build/FUSE_SRC/{FUSE_ENGINE => dshare}/multiroute.f90 (100%) create mode 100644 build/FUSE_SRC/dshare/multistate.f90 create mode 100644 build/FUSE_SRC/dshare/multistats.f90 create mode 100644 build/FUSE_SRC/physics/evap_lower_diff.f90 create mode 100644 build/FUSE_SRC/physics/evap_upper_diff.f90 create mode 100644 build/FUSE_SRC/physics/fix_ovshoot.f90 create mode 100644 build/FUSE_SRC/physics/get_parent.f90 create mode 100644 build/FUSE_SRC/physics/implicit_solve.f90 create mode 100644 build/FUSE_SRC/physics/mod_derivs_diff.f90 create mode 100644 build/FUSE_SRC/physics/mstate_rhs_diff.f90 create mode 100644 build/FUSE_SRC/physics/q_baseflow_diff.f90 create mode 100644 build/FUSE_SRC/physics/q_misscell_diff.f90 create mode 100644 build/FUSE_SRC/physics/qinterflow_diff.f90 create mode 100644 build/FUSE_SRC/physics/qpercolate_diff.f90 create mode 100644 build/FUSE_SRC/physics/qsatexcess_diff.f90 diff --git a/build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 b/build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 index 4d2c7d0..085cdb9 100644 --- a/build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 +++ b/build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 @@ -46,7 +46,8 @@ PROGRAM DISTRIBUTED_DRIVER USE multistate, 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 diff --git a/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 b/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 index bca594e..c269eff 100644 --- a/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 +++ b/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 @@ -52,6 +52,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 +98,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) @@ -245,9 +254,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); STOP 1; 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() 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/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/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/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_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/FUSE_ENGINE/model_defn.f90 b/build/FUSE_SRC/dshare/model_defn.f90 similarity index 100% rename from build/FUSE_SRC/FUSE_ENGINE/model_defn.f90 rename to build/FUSE_SRC/dshare/model_defn.f90 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..fa393ff --- /dev/null +++ b/build/FUSE_SRC/dshare/multi_flux.f90 @@ -0,0 +1,11 @@ +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 + REAL(SP) :: CURRENT_DT ! current time step (days) +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 84% rename from build/FUSE_SRC/FUSE_ENGINE/multiforce.f90 rename to build/FUSE_SRC/dshare/multiforce.f90 index 900befd..46f205a 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 diff --git a/build/FUSE_SRC/dshare/multiparam.f90 b/build/FUSE_SRC/dshare/multiparam.f90 new file mode 100644 index 0000000..7f7938d --- /dev/null +++ b/build/FUSE_SRC/dshare/multiparam.f90 @@ -0,0 +1,22 @@ +! --------------------------------------------------------------------------------------- +! 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 model_defn,ONLY:NTDH_MAX + 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..3a9a3a6 --- /dev/null +++ b/build/FUSE_SRC/dshare/multistate.f90 @@ -0,0 +1,27 @@ +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) + + ! 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/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 01cb576..5b960bc 100644 --- a/build/Makefile +++ b/build/Makefile @@ -59,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= \ @@ -83,6 +85,7 @@ NRUTIL = $(patsubst %, $(NUMREC_DIR)/%, $(FUSE_NRUTIL)) # Data modules FUSE_DATAMS= \ model_defn.f90 \ + data_types.f90 \ model_defnames.f90 \ multiconst.f90 \ multiforce.f90 \ @@ -93,7 +96,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= \ @@ -122,6 +125,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 \ @@ -208,7 +227,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) #========================================================================