From f371a87ddef304791d9f8016a9672dd9283f95f8 Mon Sep 17 00:00:00 2001 From: Martyn Clark Date: Wed, 26 Nov 2025 08:40:07 -0700 Subject: [PATCH 1/9] 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 1cf34cac4498a6866920c0f369d7584e9c75a8b3 Mon Sep 17 00:00:00 2001 From: Martyn Clark Date: Sat, 29 Nov 2025 08:32:08 -0700 Subject: [PATCH 2/9] 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) #======================================================================== From 849b836dd16bca3009373bc6333501f90ad088f7 Mon Sep 17 00:00:00 2001 From: Martyn Clark Date: Sat, 29 Nov 2025 11:20:02 -0700 Subject: [PATCH 3/9] add global data and remove print statements --- build/FUSE_SRC/FUSE_DMSL/functn.f90 | 7 ++-- build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 | 11 +++--- build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 | 42 +++++++++++++---------- build/FUSE_SRC/FUSE_ENGINE/fuse_solve.f90 | 4 +-- build/FUSE_SRC/FUSE_ENGINE/get_mbands.f90 | 7 ++-- build/FUSE_SRC/FUSE_ENGINE/mean_stats.f90 | 22 ++++++------ build/FUSE_SRC/FUSE_ENGINE/q_misscell.f90 | 3 +- build/FUSE_SRC/FUSE_NETCDF/def_output.f90 | 2 +- build/FUSE_SRC/FUSE_NETCDF/def_params.f90 | 2 +- build/FUSE_SRC/FUSE_NETCDF/def_sstats.f90 | 2 +- build/FUSE_SRC/FUSE_NETCDF/get_gforce.f90 | 2 +- build/FUSE_SRC/FUSE_NETCDF/put_output.f90 | 6 ++-- build/FUSE_SRC/dshare/globaldata.f90 | 21 ++++++++++++ build/FUSE_SRC/dshare/model_defn.f90 | 17 +++------ build/FUSE_SRC/dshare/multi_flux.f90 | 1 - build/FUSE_SRC/dshare/multiforce.f90 | 4 --- build/FUSE_SRC/dshare/multiparam.f90 | 1 - build/FUSE_SRC/dshare/multistate.f90 | 6 ---- build/Makefile | 1 + 19 files changed, 86 insertions(+), 75 deletions(-) create mode 100644 build/FUSE_SRC/dshare/globaldata.f90 diff --git a/build/FUSE_SRC/FUSE_DMSL/functn.f90 b/build/FUSE_SRC/FUSE_DMSL/functn.f90 index b7cfa89..a0e78fe 100644 --- a/build/FUSE_SRC/FUSE_DMSL/functn.f90 +++ b/build/FUSE_SRC/FUSE_DMSL/functn.f90 @@ -9,13 +9,13 @@ FUNCTION FUNCTN(NOPT,A) ! Wrapper for SCE (used to compute the objective function) ! --------------------------------------------------------------------------------------- USE nrtype ! variable types, etc. -USE FUSE_RMSE_MODULE ! run model and compute the root mean squared error -USE multiforce, only: ncid_forc ! NetCDF forcing file ID +USE FUSE_RMSE_MODULE ! run model and compute the root mean squared error +USE multiforce, only: ncid_forc ! NetCDF forcing file ID IMPLICIT NONE ! input INTEGER(I4B) :: NOPT ! number of parameters -REAL(MSP), DIMENSION(100), INTENT(IN) :: A ! model parameter set - can be bumped up to 100 elements +REAL(MSP), DIMENSION(100), INTENT(IN) :: A ! model parameter set - can be bumped up to 100 elements ! internal REAL(SP), DIMENSION(:), ALLOCATABLE :: SCE_PAR ! sce parameter set @@ -39,7 +39,6 @@ FUNCTION FUNCTN(NOPT,A) ! deallocate parameter set DEALLOCATE(SCE_PAR, STAT=IERR); IF (IERR.NE.0) STOP ' problem deallocating space ' -print *, 'RMSE =', RMSE ! save objective function value FUNCTN = RMSE diff --git a/build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 b/build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 index 085cdb9..582e5d2 100644 --- a/build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 +++ b/build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 @@ -23,6 +23,7 @@ PROGRAM DISTRIBUTED_DRIVER ! data modules USE model_defn,nstateFUSE=>nstate ! model definition structures USE model_defnames ! defines the integer model options +USE globaldata, ONLY: isPrint ! flag for printing progress to screen USE multiforce, ONLY: forcefile,vname_aprecip ! model forcing structures USE multiforce, ONLY: AFORCE, aValid ! time series of lumped forcing/response data USE multiforce, ONLY: nspat1, nspat2 ! grid dimensions @@ -43,7 +44,7 @@ PROGRAM DISTRIBUTED_DRIVER USE multiforce, only: ncid_forc ! NetCDF forcing file ID USE multiforce, only: ncid_var ! NetCDF forcing variable ID -USE multistate, only: ncid_out ! NetCDF output file ID +USE globaldata, only: ncid_out ! NetCDF output file ID USE multibands ! basin band stuctures USE data_types, ONLY: PARATT ! data type for metadata @@ -346,9 +347,9 @@ PROGRAM DISTRIBUTED_DRIVER ! assign algorithmic control parameters for SCE ! convert characters to interger/MSP - READ (MAXN_STR,*) MAXN ! maximum number of trials before optimization is terminated + READ (MAXN_STR,*) MAXN ! maximum number of trials before optimization is terminated READ (KSTOP_STR,*) KSTOP ! number of shuffling loops the value must change by PCENTO (MAX=9) - READ (PCENTO_STR,*) PCENTO ! the percentage + READ (PCENTO_STR,*) PCENTO ! the percentage PRINT *, 'SCE parameters read from file manager:' PRINT *, 'Maximum number of trials before SCE optimization is stopped (MAXN) = ', MAXN_STR @@ -439,11 +440,13 @@ PROGRAM DISTRIBUTED_DRIVER FNAME_ASCII = TRIM(OUTPUT_PATH)//TRIM(dom_id)//'_'//TRIM(FMODEL_ID)//'_sce_output.txt' + ! turn off printing to screen + isPrint = .false. + ! convert from SP used in FUSE to MSP used in SCE ALLOCATE(APAR_MSP(NUMPAR),BL_MSP(NUMPAR),BU_MSP(NUMPAR),URAND_MSP(NUMPAR)) APAR_MSP=APAR - PRINT *, 'BL=',BL BL_MSP=BL BU_MSP=BU URAND_MSP=URAND diff --git a/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 b/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 index c269eff..7996ebb 100644 --- a/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 +++ b/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 @@ -22,6 +22,9 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG ! data modules USE model_defn, ONLY:NSTATE,SMODL ! number of state variables USE model_defnames ! integer model definitions + USE globaldata, ONLY: isPrint ! flag for printing progress to screen + USE globaldata, ONLY: fracstate0 ! fraction of initial state (used for initialization) + USE globaldata, ONLY: NA_VALUE, NA_VALUE_SP ! NA_VALUE for the forcing USE multiparam, ONLY: LPARAM,NUMPAR,MPARAM ! list of model parameters USE multiforce, ONLY: MFORCE,AFORCE,DELTIM,ISTART ! model forcing data USE multiforce, ONLY: numtim_in, itim_in ! length of input time series and associated index @@ -34,9 +37,7 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG USE multiforce, ONLY:nspat1,nspat2 ! spatial dimensions USE multiforce, ONLY:ncid_var ! NetCDF ID for forcing variables USE multiforce, ONLY:gForce,gForce_3d ! gridded forcing data - USE multistate, ONLY:fracstate0,TSTATE,MSTATE,FSTATE,& ! model states - HSTATE ! model states (continued) - USE multiforce, ONLY:NA_VALUE, NA_VALUE_SP ! NA_VALUE for the forcing + USE multistate, ONLY:TSTATE,MSTATE,FSTATE,HSTATE ! model state variables USE multistate, ONLY:gState,gState_3d ! gridded state variables USE multiroute, ONLY:MROUTE,AROUTE,AROUTE_3d ! routed runoff USE multistats, ONLY:MSTATS,PCOUNT,MOD_IX ! access model statistics; counter for param set @@ -116,13 +117,16 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG ! add parameter set to the data structure CALL PUT_PARSET(XPAR) - PRINT *, 'Parameter set added to data structure:' - PRINT *, XPAR + if(isPrint) PRINT *, 'Parameter set added to data structure:' + if(isPrint) PRINT *, XPAR ! compute derived model parameters (bucket sizes, etc.) CALL PAR_DERIVE(ERR,MESSAGE) IF (ERR.NE.0) WRITE(*,*) TRIM(MESSAGE); IF (ERR.GT.0) STOP + if(isPrint) PRINT *, 'Writing parameter values...' + CALL PUT_PARAMS(PCOUNT) + ! initialize model states over the 2D gridded domain (1x1 domain in catchment mode) DO iSpat2=1,nSpat2 DO iSpat1=1,nSpat1 @@ -130,10 +134,10 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG gState_3d(iSpat1,iSpat2,1) = FSTATE ! put the state into first time step of 3D structure END DO END DO - PRINT *, 'Model states initialized over the 2D gridded domain' + if(isPrint) PRINT *, 'Model states initialized over the 2D gridded domain' ! initialize elevations bands if snow module is on - PRINT *, 'N_BANDS =', N_BANDS + if(isPrint) PRINT *, 'N_BANDS =', N_BANDS IF (SMODL%iSNOWM.EQ.iopt_temp_index) THEN DO iSpat2=1,nSpat2 @@ -146,7 +150,7 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG END DO END DO END DO - PRINT *, 'Snow states initiatlized over the 2D gridded domain ' + if(isPrint) PRINT *, 'Snow states initiatlized over the 2D gridded domain ' ENDIF ! allocate 3d data structure for fluxes @@ -186,10 +190,10 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG numtim_sub_cur=MIN(numtim_sub,numtim_sim-itim_sim+1) ! load forcing for desired period into gForce_3d - PRINT *, 'New subperiod: loading forcing for ',numtim_sub_cur,' time steps' + if(isPrint) PRINT *, 'New subperiod: loading forcing for ',numtim_sub_cur,' time steps' CALL get_gforce_3d(itim_in,numtim_sub_cur,ncid_forc,err,message) IF(err/=0)THEN; WRITE(*,*) 'Error while extracting 3d forcing'; STOP; ENDIF - PRINT *, 'Forcing loaded. Running FUSE...' + if(isPrint) PRINT *, 'Forcing loaded. Running FUSE...' ENDIF @@ -345,15 +349,15 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG ! if end of subperiod: write to output file and save states IF(itim_sub.EQ.numtim_sub_cur)THEN - PRINT *, 'End of subperiod reached:' + if(isPrint) PRINT *, 'End of subperiod reached:' ! write model output IF (OUTPUT_FLAG) THEN - PRINT *, 'Write output for ',numtim_sub_cur,' time steps starting at indice', itim_sim-numtim_sub_cur+1 + if(isPrint) PRINT *, 'Write output for ',numtim_sub_cur,' time steps starting at indice', itim_sim-numtim_sub_cur+1 CALL PUT_GOUTPUT_3D(itim_sim-numtim_sub_cur+1,itim_in-numtim_sub_cur+1,numtim_sub_cur,IPSET) - PRINT *, 'Done writing output' + if(isPrint) PRINT *, 'Done writing output' ELSE - PRINT *, 'OUTPUT_FLAG is set on FALSE, no output written' + if(isPrint) PRINT *, 'OUTPUT_FLAG is set on FALSE, no output written' END IF ! TODO: set gState_3d and MBANDS_VAR_4d to NA @@ -382,20 +386,20 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG ! get timing information CALL CPU_TIME(T2) - WRITE(*,*) "TIME ELAPSED = ", t2-t1 + if(isPrint) WRITE(*,*) "TIME ELAPSED = ", t2-t1 ! calculate mean summary statistics IF(.NOT.GRID_FLAG)THEN - PRINT *, 'Calculating performance metrics...' + if(isPrint) PRINT *, 'Calculating performance metrics...' CALL MEAN_STATS() RMSE = MSTATS%RAW_RMSE + print*, "NSE = ", MSTATS%NASH_SUTT + ENDIF - PRINT *, 'Writing parameter values...' - CALL PUT_PARAMS(PCOUNT) - PRINT *, 'Writing model statistics...' + if(isPrint) PRINT *, 'Writing model statistics...' CALL PUT_SSTATS(PCOUNT) ! deallocate vectors diff --git a/build/FUSE_SRC/FUSE_ENGINE/fuse_solve.f90 b/build/FUSE_SRC/FUSE_ENGINE/fuse_solve.f90 index dd4ec5b..5ae59e3 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/fuse_solve.f90 +++ b/build/FUSE_SRC/FUSE_ENGINE/fuse_solve.f90 @@ -10,8 +10,8 @@ SUBROUTINE FUSE_SOLVE(CALCDSDT,IE_SOLVE,SI_SOLVE,B_IMPOSE,AVG_FLUX,ADD_FLUX,NEWS ! (6) add fluxes from accepted sub-steps to the total timestep flux ! (7) estimate state at end of a full step, based on sum of fluxes USE nrtype ! variable definitions, etc. -USE multi_flux, ONLY: M_FLUX,FLUX_0,FLUX_1,W_FLUX,& ! model fluxes - CURRENT_DT ! model fluxes (continued) +USE globaldata, ONLY: CURRENT_DT +USE multi_flux, ONLY: M_FLUX,FLUX_0,FLUX_1,W_FLUX ! model fluxes USE multistate, ONLY: FSTATE,MSTATE,BSTATE,ESTATE,& ! model states DY_DT,DYDT_0,DYDT_1,HSTATE ! model states (continued) USE fminln, ONLY: fmin_x0p,fmin_dtp,fmin_dt2p,fmin_dseep ! variables used for residual vector in IE diff --git a/build/FUSE_SRC/FUSE_ENGINE/get_mbands.f90 b/build/FUSE_SRC/FUSE_ENGINE/get_mbands.f90 index f05a6ba..efd006c 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/get_mbands.f90 +++ b/build/FUSE_SRC/FUSE_ENGINE/get_mbands.f90 @@ -147,11 +147,12 @@ SUBROUTINE GET_MBANDS_INFO(ELEV_BANDS_NC,err,message) use nrtype,only:I4B,LGT,SP use utilities_dmsl_kit_FUSE,only:getSpareUnit,stripTrailString USE fuse_fileManager,only:INPUT_PATH,SETNGS_PATH ! defines data directory -USE fuse_fileManager,only:MBANDS_NC ! defines elevation bands +USE fuse_fileManager,only:MBANDS_NC ! defines elevation bands USE multibands,only:N_BANDS,MBANDS,MBANDS_INFO_3d,Z_FORCING,& - Z_FORCING_grid,elev_mask ! model band structures -USE multiforce,only:nspat1,nspat2,startSpat2,NA_VALUE_SP ! dimension lengths, na_value + Z_FORCING_grid,elev_mask ! model band structures +USE multiforce,only:nspat1,nspat2,startSpat2 ! dimension lengths +USE globaldata,only:NA_VALUE_SP IMPLICIT NONE ! dummies diff --git a/build/FUSE_SRC/FUSE_ENGINE/mean_stats.f90 b/build/FUSE_SRC/FUSE_ENGINE/mean_stats.f90 index 106c3e3..28fde13 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/mean_stats.f90 +++ b/build/FUSE_SRC/FUSE_ENGINE/mean_stats.f90 @@ -19,6 +19,8 @@ SUBROUTINE MEAN_STATS() USE multiroute ! routed runoff USE multi_flux ! fluxes USE multistats ! summary statistics +USE globaldata, ONLY: isPrint ! flag for printing progress to screen +USE globaldata, only: NA_VALUE_SP ! missing value USE model_numerix ! model numerix parameters and data IMPLICIT NONE @@ -51,7 +53,7 @@ SUBROUTINE MEAN_STATS() ! --------------------------------------------------------------------------------------- ! define sample size NS = eval_end-eval_beg+1 -PRINT *, 'Number of time steps in evaluation period (EP) = ', NS +if(isPrint) PRINT *, 'Number of time steps in evaluation period (EP) = ', NS ! allocate space for observed and simulated runoff ALLOCATE(QOBS(NS),QOBS_MASK(NS),QSIM(NS),STAT=IERR) @@ -63,10 +65,10 @@ SUBROUTINE MEAN_STATS() QOBS = aValid(1,1,eval_beg-sim_beg+1:eval_end-sim_beg+1)%OBSQ ! check for missing QOBS values -QOBS_MASK = QOBS.ne.REAL(NA_VALUE, KIND(SP)) ! find the time steps for which QOBS is available -NUM_AVAIL = COUNT(QOBS_MASK) ! number of time steps for which QOBS is available +QOBS_MASK = QOBS.ne.NA_VALUE_SP ! find the time steps for which QOBS is available +NUM_AVAIL = COUNT(QOBS_MASK) ! number of time steps for which QOBS is available -PRINT *, 'Number of time steps with observed streamflow in EP = ', NUM_AVAIL +if(isPrint) PRINT *, 'Number of time steps with observed streamflow in EP = ', NUM_AVAIL IF (NUM_AVAIL.EQ.0) THEN @@ -80,11 +82,11 @@ SUBROUTINE MEAN_STATS() ALLOCATE(QOBS_AVAIL(NUM_AVAIL),QSIM_AVAIL(NUM_AVAIL),DOBS(NUM_AVAIL),DSIM(NUM_AVAIL),RAWD(NUM_AVAIL),LOGD(NUM_AVAIL),STAT=IERR) QOBS_AVAIL=PACK(QOBS,QOBS_MASK,QOBS_AVAIL) ! moves QOBS time steps indicated by QOBS_MASK to QOBS_AVAIL, - ! if no values is missing (i.e. NS = NUM_AVAIL) then QOBS_AVAIL - ! should be a copy of QOBS + ! if no values is missing (i.e. NS = NUM_AVAIL) then QOBS_AVAIL + ! should be a copy of QOBS QSIM_AVAIL=PACK(QSIM,QOBS_MASK,QSIM_AVAIL) ! moves QSIM time steps indicated by QOBS_MASK to QSIM_AVAIL - ! if no values is missing (i.e. NS = NUM_AVAIL) then QSIM_AVAIL - ! should be a copy of QSIM + ! if no values is missing (i.e. NS = NUM_AVAIL) then QSIM_AVAIL + ! should be a copy of QSIM ! compute mean XB_OBS = SUM(QOBS_AVAIL(:)) / REAL(NUM_AVAIL, KIND(SP)) @@ -130,8 +132,8 @@ SUBROUTINE MEAN_STATS() END IF -PRINT *, 'NSE = ', MSTATS%NASH_SUTT -PRINT *, 'RAW_RMSE = ', MSTATS%RAW_RMSE +if(isPrint) PRINT *, 'NSE = ', MSTATS%NASH_SUTT +if(isPrint) PRINT *, 'RAW_RMSE = ', MSTATS%RAW_RMSE ! --------------------------------------------------------------------------------------- ! (3§) COMPUTE STATISTICS ON NUMERICAL ACCURACY AND EFFICIENCY diff --git a/build/FUSE_SRC/FUSE_ENGINE/q_misscell.f90 b/build/FUSE_SRC/FUSE_ENGINE/q_misscell.f90 index b40328a..b84bad5 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/q_misscell.f90 +++ b/build/FUSE_SRC/FUSE_ENGINE/q_misscell.f90 @@ -20,9 +20,10 @@ SUBROUTINE Q_MISSCELL() USE nrtype ! variable types, etc. USE model_defn ! model definition structure USE model_defnames +USE globaldata, ONLY: CURRENT_DT USE multiparam, ONLY: MPARAM,DPARAM ! model parameters USE multistate, ONLY: MSTATE,TSTATE ! model states -USE multi_flux, ONLY: M_FLUX,CURRENT_DT ! model fluxes +USE multi_flux, ONLY: M_FLUX ! model fluxes USE model_numerix ! access model numerix decisions IMPLICIT NONE REAL(SP) :: LOGISMOOTH ! FUNCTION logistic smoothing diff --git a/build/FUSE_SRC/FUSE_NETCDF/def_output.f90 b/build/FUSE_SRC/FUSE_NETCDF/def_output.f90 index f04bbb5..ebd3dc9 100644 --- a/build/FUSE_SRC/FUSE_NETCDF/def_output.f90 +++ b/build/FUSE_SRC/FUSE_NETCDF/def_output.f90 @@ -19,7 +19,7 @@ SUBROUTINE DEF_OUTPUT(nSpat1,nSpat2,NPSET,NTIM) USE multiforce, only: name_psets,time_steps ! dimension arrays USE multiforce, only: latUnits,lonUnits ! units string USE multiforce, only: timeUnits ! units string - USE multistate, only: ncid_out ! NetCDF output file ID + USE globaldata, only: ncid_out ! NetCDF output file ID IMPLICIT NONE diff --git a/build/FUSE_SRC/FUSE_NETCDF/def_params.f90 b/build/FUSE_SRC/FUSE_NETCDF/def_params.f90 index 46b2cdb..0c9ea24 100644 --- a/build/FUSE_SRC/FUSE_NETCDF/def_params.f90 +++ b/build/FUSE_SRC/FUSE_NETCDF/def_params.f90 @@ -13,7 +13,7 @@ SUBROUTINE DEF_PARAMS(NPAR) USE model_defn ! model definition (includes filename) USE metaparams ! metadata for all model parameters USE multistats, ONLY: MSTATS ! model statistics structure -USE multistate, only: ncid_out ! NetCDF output file ID +USE globaldata, only: ncid_out ! NetCDF output file ID IMPLICIT NONE ! input INTEGER(I4B), INTENT(IN) :: NPAR ! number of parameter sets diff --git a/build/FUSE_SRC/FUSE_NETCDF/def_sstats.f90 b/build/FUSE_SRC/FUSE_NETCDF/def_sstats.f90 index 2d4b4ac..f75b1e2 100644 --- a/build/FUSE_SRC/FUSE_NETCDF/def_sstats.f90 +++ b/build/FUSE_SRC/FUSE_NETCDF/def_sstats.f90 @@ -12,7 +12,7 @@ SUBROUTINE DEF_SSTATS() USE model_defn ! model definition (includes filename) USE meta_stats ! metadata for summary statistics USE model_numerix ! model numerix decisions -USE multistate, only: ncid_out ! NetCDF output file ID +USE globaldata, only: ncid_out ! NetCDF output file ID IMPLICIT NONE ! internal INTEGER(I4B) :: IERR ! error code; NetCDF ID diff --git a/build/FUSE_SRC/FUSE_NETCDF/get_gforce.f90 b/build/FUSE_SRC/FUSE_NETCDF/get_gforce.f90 index f990a04..1445c3e 100644 --- a/build/FUSE_SRC/FUSE_NETCDF/get_gforce.f90 +++ b/build/FUSE_SRC/FUSE_NETCDF/get_gforce.f90 @@ -2,6 +2,7 @@ module get_gforce_module USE nrtype USE netcdf USE time_io +USE globaldata, only: NA_VALUE, NA_VALUE_SP ! missing value implicit none private public::read_ginfo @@ -75,7 +76,6 @@ SUBROUTINE read_ginfo(ncid,ierr,message) USE multiforce,only:latUnits,lonUnits,timeUnits ! units string for time USE multiforce,only:vname_dtime ! variable name: time sice reference time USE multiforce, only: nForce, nInput ! number of parameter set and their names - USE multiforce, only: NA_VALUE ! NA_VALUE for the forcing #ifdef __MPI__ use mpi diff --git a/build/FUSE_SRC/FUSE_NETCDF/put_output.f90 b/build/FUSE_SRC/FUSE_NETCDF/put_output.f90 index 8d1b13e..99d6676 100644 --- a/build/FUSE_SRC/FUSE_NETCDF/put_output.f90 +++ b/build/FUSE_SRC/FUSE_NETCDF/put_output.f90 @@ -15,7 +15,7 @@ SUBROUTINE PUT_OUTPUT(iSpat1,iSpat2,ITIM,IMOD,IPAR) USE varextract_module ! interface for the function to extract variables USE fuse_fileManager,only: Q_ONLY ! only write streamflow to output file? USE multiforce,ONLY: timDat ! time data - USE multistate, only: ncid_out ! NetCDF output file ID + USE globaldata, only: ncid_out ! NetCDF output file ID IMPLICIT NONE ! input @@ -95,10 +95,10 @@ SUBROUTINE PUT_GOUTPUT_3D(istart_sim,istart_in,numtim,IPSET) USE fuse_fileManager,only: Q_ONLY ! only write streamflow to output file? USE multiforce, ONLY: timDat,time_steps ! time data - USE multistate, only: ncid_out ! NetCDF output file ID USE multiforce, ONLY: nspat1,nspat2,startSpat2 ! spatial dimensions USE multiforce, ONLY: gForce_3d ! test only - USE multiforce, only: GRID_FLAG ! .true. if distributed + USE multiforce, only: GRID_FLAG ! .true. if distributed + USE globaldata, only: ncid_out ! NetCDF output file ID IMPLICIT NONE diff --git a/build/FUSE_SRC/dshare/globaldata.f90 b/build/FUSE_SRC/dshare/globaldata.f90 new file mode 100644 index 0000000..0029c83 --- /dev/null +++ b/build/FUSE_SRC/dshare/globaldata.f90 @@ -0,0 +1,21 @@ +MODULE globaldata + + USE nrtype + + ! time step + REAL(SP), save :: CURRENT_DT ! current time step (days) + + ! missing values + INTEGER(I4B),PARAMETER :: NA_VALUE=-9999 ! integer designating missing values - TODO: retrieve from NetCDF file + REAL(SP),PARAMETER :: NA_VALUE_SP=-9999_sp ! integer designating missing values - TODO: retrieve from NetCDF file + + ! NetCDF + integer(i4b), save :: ncid_out=-1 ! NetCDF output file ID + + ! initial store fraction (initialization) + real(sp), parameter :: fracState0=0.25_sp + + ! print flag + logical(lgt) :: isPrint=.true. + +end MODULE globaldata diff --git a/build/FUSE_SRC/dshare/model_defn.f90 b/build/FUSE_SRC/dshare/model_defn.f90 index 9a0c80a..0dcd28b 100644 --- a/build/FUSE_SRC/dshare/model_defn.f90 +++ b/build/FUSE_SRC/dshare/model_defn.f90 @@ -27,36 +27,27 @@ MODULE model_defn TYPE UMODEL INTEGER(I4B) :: MODIX ! model index CHARACTER(LEN=256) :: MNAME ! model name -! CHARACTER(LEN=16) :: RFERR ! rainfall error INTEGER(I4B) :: iRFERR -! CHARACTER(LEN=16) :: ARCH1 ! upper-layer architecture INTEGER(I4B) :: iARCH1 -! CHARACTER(LEN=16) :: ARCH2 ! lower-layer architecture INTEGER(I4B) :: iARCH2 -! CHARACTER(LEN=16) :: QSURF ! surface runoff INTEGER(I4B) :: iQSURF -! CHARACTER(LEN=16) :: QPERC ! percolation INTEGER(I4B) :: iQPERC -! CHARACTER(LEN=16) :: ESOIL ! evaporation INTEGER(I4B) :: iESOIL -! CHARACTER(LEN=16) :: QINTF ! interflow INTEGER(I4B) :: iQINTF -! CHARACTER(LEN=16) :: Q_TDH ! time delay in runoff INTEGER(I4B) :: iQ_TDH INTEGER(I4B) :: iSNOWM ! snow - END TYPE UMODEL + END TYPE UMODEL ! structure to hold model state names TYPE SNAMES -! CHARACTER(LEN=8) :: SNAME ! state name INTEGER(I4B) :: iSNAME ! integer value of state name END TYPE SNAMES ! structure to hold model flux names TYPE FNAMES CHARACTER(LEN=16) :: FNAME ! state name END TYPE FNAMES -! max steps in routing function - INTEGER(I4B),PARAMETER::NTDH_MAX=500 -! model definitions + ! max steps in routing function + INTEGER(I4B),PARAMETER::NTDH_MAX=500 + ! model definitions CHARACTER(LEN=256) :: FNAME_NETCDF_RUNS ! NETCDF output filename for model runs CHARACTER(LEN=256) :: FNAME_NETCDF_PARA ! NETCDF output filename for model parameters CHARACTER(LEN=256) :: FNAME_NETCDF_PARA_SCE ! NETCDF output filename for model parameters produced by SCE diff --git a/build/FUSE_SRC/dshare/multi_flux.f90 b/build/FUSE_SRC/dshare/multi_flux.f90 index fa393ff..b00bb06 100644 --- a/build/FUSE_SRC/dshare/multi_flux.f90 +++ b/build/FUSE_SRC/dshare/multi_flux.f90 @@ -7,5 +7,4 @@ MODULE multi_flux 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/dshare/multiforce.f90 b/build/FUSE_SRC/dshare/multiforce.f90 index 46f205a..468f649 100644 --- a/build/FUSE_SRC/dshare/multiforce.f90 +++ b/build/FUSE_SRC/dshare/multiforce.f90 @@ -124,9 +124,5 @@ MODULE multiforce REAL(sp) :: amult_pet=-1._dp ! convert potential ET to mm/day REAL(sp) :: amult_q=-1._dp ! convert runoff to mm/day - ! missing values - INTEGER(I4B),PARAMETER :: NA_VALUE=-9999 ! integer designating missing values - TODO: retrieve from NetCDF file - REAL(SP),PARAMETER :: NA_VALUE_SP=-9999 ! integer designating missing values - TODO: retrieve from NetCDF file - ! -------------------------------------------------------------------------------------- END MODULE multiforce diff --git a/build/FUSE_SRC/dshare/multiparam.f90 b/build/FUSE_SRC/dshare/multiparam.f90 index 7f7938d..cfaa939 100644 --- a/build/FUSE_SRC/dshare/multiparam.f90 +++ b/build/FUSE_SRC/dshare/multiparam.f90 @@ -7,7 +7,6 @@ ! --------------------------------------------------------------------------------------- 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 diff --git a/build/FUSE_SRC/dshare/multistate.f90 b/build/FUSE_SRC/dshare/multistate.f90 index 3a9a3a6..f7724f0 100644 --- a/build/FUSE_SRC/dshare/multistate.f90 +++ b/build/FUSE_SRC/dshare/multistate.f90 @@ -18,10 +18,4 @@ MODULE multistate 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/Makefile b/build/Makefile index 5b960bc..9c2edfb 100644 --- a/build/Makefile +++ b/build/Makefile @@ -87,6 +87,7 @@ FUSE_DATAMS= \ model_defn.f90 \ data_types.f90 \ model_defnames.f90 \ + globaldata.f90 \ multiconst.f90 \ multiforce.f90 \ multibands.f90 \ From f359f5d1dabee8512892c339e285db8508a25bac Mon Sep 17 00:00:00 2001 From: Martyn Clark Date: Sat, 29 Nov 2025 11:36:58 -0700 Subject: [PATCH 4/9] add random seed to obtain same results for different optimization trials --- build/FUSE_SRC/FUSE_DMSL/functn.f90 | 4 ++++ build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 | 9 +++++++-- build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 | 3 ++- build/FUSE_SRC/dshare/globaldata.f90 | 3 +++ 4 files changed, 16 insertions(+), 3 deletions(-) diff --git a/build/FUSE_SRC/FUSE_DMSL/functn.f90 b/build/FUSE_SRC/FUSE_DMSL/functn.f90 index a0e78fe..1d2e813 100644 --- a/build/FUSE_SRC/FUSE_DMSL/functn.f90 +++ b/build/FUSE_SRC/FUSE_DMSL/functn.f90 @@ -11,6 +11,7 @@ FUNCTION FUNCTN(NOPT,A) USE nrtype ! variable types, etc. USE FUSE_RMSE_MODULE ! run model and compute the root mean squared error USE multiforce, only: ncid_forc ! NetCDF forcing file ID +USE globaldata, only: nFUSE_eval ! # fuse evaluations IMPLICIT NONE ! input @@ -29,6 +30,9 @@ FUNCTION FUNCTN(NOPT,A) REAL(MSP) :: FUNCTN ! objective function value ! --------------------------------------------------------------------------------------- + +nFUSE_eval = nFUSE_eval + 1 + ! get SCE parameter set ALLOCATE(SCE_PAR(NOPT), STAT=IERR); IF (IERR.NE.0) STOP ' problem allocating space ' SCE_PAR(1:NOPT) = A(1:NOPT) ! convert from MSP used in SCE to SP used in FUSE diff --git a/build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 b/build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 index 582e5d2..1060413 100644 --- a/build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 +++ b/build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 @@ -24,6 +24,7 @@ PROGRAM DISTRIBUTED_DRIVER USE model_defn,nstateFUSE=>nstate ! model definition structures USE model_defnames ! defines the integer model options USE globaldata, ONLY: isPrint ! flag for printing progress to screen +USE globaldata, only: nFUSE_eval ! number of fuse evaluations USE multiforce, ONLY: forcefile,vname_aprecip ! model forcing structures USE multiforce, ONLY: AFORCE, aValid ! time series of lumped forcing/response data USE multiforce, ONLY: nspat1, nspat2 ! grid dimensions @@ -440,8 +441,9 @@ PROGRAM DISTRIBUTED_DRIVER FNAME_ASCII = TRIM(OUTPUT_PATH)//TRIM(dom_id)//'_'//TRIM(FMODEL_ID)//'_sce_output.txt' - ! turn off printing to screen - isPrint = .false. + ! printing + isPrint = .false. ! ! turn off printing to screen + nFUSE_eval = 0 ! number of fuse eevaluations ! convert from SP used in FUSE to MSP used in SCE ALLOCATE(APAR_MSP(NUMPAR),BL_MSP(NUMPAR),BU_MSP(NUMPAR),URAND_MSP(NUMPAR)) @@ -451,6 +453,9 @@ PROGRAM DISTRIBUTED_DRIVER BU_MSP=BU URAND_MSP=URAND + ! set random seed + ISEED = 1 + ! open up ASCII output file print *, 'Creating SCE output file:', trim(FNAME_ASCII) ISCE = 96; OPEN(ISCE,FILE=TRIM(FNAME_ASCII)) diff --git a/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 b/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 index 7996ebb..76226d0 100644 --- a/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 +++ b/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 @@ -23,6 +23,7 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG USE model_defn, ONLY:NSTATE,SMODL ! number of state variables USE model_defnames ! integer model definitions USE globaldata, ONLY: isPrint ! flag for printing progress to screen + USE globaldata, only: nFUSE_eval ! number of fuse evaluations USE globaldata, ONLY: fracstate0 ! fraction of initial state (used for initialization) USE globaldata, ONLY: NA_VALUE, NA_VALUE_SP ! NA_VALUE for the forcing USE multiparam, ONLY: LPARAM,NUMPAR,MPARAM ! list of model parameters @@ -395,7 +396,7 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG CALL MEAN_STATS() RMSE = MSTATS%RAW_RMSE - print*, "NSE = ", MSTATS%NASH_SUTT + write(*,'(i6,1x,a6,1x,f12.6,1x)') nFUSE_eval, "NSE = ", MSTATS%NASH_SUTT ENDIF diff --git a/build/FUSE_SRC/dshare/globaldata.f90 b/build/FUSE_SRC/dshare/globaldata.f90 index 0029c83..d9edb1d 100644 --- a/build/FUSE_SRC/dshare/globaldata.f90 +++ b/build/FUSE_SRC/dshare/globaldata.f90 @@ -18,4 +18,7 @@ MODULE globaldata ! print flag logical(lgt) :: isPrint=.true. + ! number of fuse evaluations + integer(i4b), save :: nFUSE_eval + end MODULE globaldata From fa48f628dcc0f130c7f8a8e66ddc7f64b73bb17c Mon Sep 17 00:00:00 2001 From: Martyn Clark Date: Sat, 29 Nov 2025 12:24:55 -0700 Subject: [PATCH 5/9] fix run_pre option --- build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 | 86 ++++++++---------------- 1 file changed, 27 insertions(+), 59 deletions(-) diff --git a/build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 b/build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 index 1060413..0c54a36 100644 --- a/build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 +++ b/build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 @@ -41,7 +41,6 @@ PROGRAM DISTRIBUTED_DRIVER USE multiforce, only: sim_beg,sim_end ! timestep indices USE multiforce, only: eval_beg,eval_end ! timestep indices USE multiforce, only: SUB_PERIODS_FLAG ! .true. if subperiods are used to run FUSE -USE multiforce, only: NUMPSET,name_psets ! number of parameter set and their names USE multiforce, only: ncid_forc ! NetCDF forcing file ID USE multiforce, only: ncid_var ! NetCDF forcing variable ID @@ -86,14 +85,15 @@ PROGRAM DISTRIBUTED_DRIVER CHARACTER(LEN=256) :: DatString ! file manager CHARACTER(LEN=256) :: dom_id ! ID of the domain CHARACTER(LEN=10) :: fuse_mode=' ' ! fuse execution mode (run_def, run_best, run_pre, calib_sce) -CHARACTER(LEN=256) :: file_para_list ! txt file containing list of parameter sets +CHARACTER(LEN=256) :: file_param ! name of parameter file +CHARACTER(LEN=10) :: index_param ! index of desired parameter set ! --------------------------------------------------------------------------------------- ! SETUP MODELS FOR SIMULATION -- POPULATE DATA STRUCTURES ! --------------------------------------------------------------------------------------- ! fuse_file_manager -CHARACTER(LEN=1024) :: FFMFILE ! name of fuse_file_manager file -CHARACTER(LEN=1024) :: ELEV_BANDS_NC ! name of NetCDF file for elevation bands +CHARACTER(LEN=1024) :: FFMFILE ! name of fuse_file_manager file +CHARACTER(LEN=1024) :: ELEV_BANDS_NC ! name of NetCDF file for elevation bands ! get model forcing data INTEGER(I4B) :: NTIM ! number of time steps - still needed ? INTEGER(I4B) :: INFERN_START ! start of inference period - still needed? @@ -121,7 +121,7 @@ PROGRAM DISTRIBUTED_DRIVER ! --------------------------------------------------------------------------------------- INTEGER(I4B) :: ITIM ! loop thru time steps INTEGER(I4B) :: IPAR ! loop thru model parameters -INTEGER(I4B) :: IPSET ! loop thru model parameter sets +INTEGER(I4B) :: IPSET ! index of desired model parameter set TYPE(PARATT) :: PARAM_META ! parameter metadata (model parameters) REAL(SP), DIMENSION(:), ALLOCATABLE :: BL ! vector of lower parameter bounds REAL(SP), DIMENSION(:), ALLOCATABLE :: BU ! vector of upper parameter bounds @@ -129,7 +129,7 @@ PROGRAM DISTRIBUTED_DRIVER INTEGER(KIND=4) :: ISEED ! seed for the random sequence REAL(KIND=4),DIMENSION(:), ALLOCATABLE :: URAND ! vector of quasi-random numbers U[0,1] REAL(SP) :: RMSE ! error from the simulation - +integer(i4b) :: NUMPSET ! number of parameter sets ! --------------------------------------------------------------------------------------- ! SCE VARIABLES ! --------------------------------------------------------------------------------------- @@ -179,22 +179,28 @@ PROGRAM DISTRIBUTED_DRIVER CALL GETARG(1,DatString) ! string defining forcinginfo file CALL GETARG(2,dom_id) ! ID of the domain CALL GETARG(3,fuse_mode) ! fuse execution mode (run_def, run_best, calib_sce) -IF(TRIM(fuse_mode).EQ.'run_pre') CALL GETARG(4,file_para_list) ! fuse execution mode txt file containing list of parameter sets +IF(TRIM(fuse_mode).EQ.'run_pre')then + CALL GETARG(4,file_param) ! name of parameter file + CALL GETARG(5,index_param) ! index of desired parameter set + IF(LEN_TRIM(index_param).EQ.0) IPSET = 1 + IF(LEN_TRIM(index_param).GT.0) read(index_param,*) IPSET +ENDIF ! check command-line arguments IF (LEN_TRIM(DatString).EQ.0) STOP '1st command-line argument is missing (fileManager)' IF (LEN_TRIM(dom_id).EQ.0) STOP '2nd command-line argument is missing (dom_id)' IF (LEN_TRIM(fuse_mode).EQ.0) STOP '3rd command-line argument is missing (fuse_mode)' IF(TRIM(fuse_mode).EQ.'run_pre')THEN - IF(LEN_TRIM(file_para_list).EQ.0) STOP '4th command-line argument is missing (file_para_list) and is required in mode run_pre' + IF(LEN_TRIM(file_param).EQ.0) STOP '4th command-line argument is missing (file_param) and is required in mode run_pre' ENDIF ! print command-line arguments -print*, '1st command-line argument (fileManager) = ', trim(DatString) -print*, '2nd command-line argument (dom_id) = ', trim(dom_id) -print*, '3rd command-line argument (fuse_mode) = ', fuse_mode +print*, '1st command-line argument (fileManager) = ', trim(DatString) +print*, '2nd command-line argument (dom_id) = ', trim(dom_id) +print*, '3rd command-line argument (fuse_mode) = ', fuse_mode IF(TRIM(fuse_mode).EQ.'run_pre')THEN - print*, '4th command-line argument (file_para_list) = ', file_para_list + print*, '4th command-line argument (file_param) = ', file_param + print*, '5th command-line argument (index_param) = ', IPSET ENDIF ! --------------------------------------------------------------------------------------- @@ -301,45 +307,14 @@ PROGRAM DISTRIBUTED_DRIVER FNAME_NETCDF_PARA = TRIM(OUTPUT_PATH)//TRIM(dom_id)//'_'//TRIM(FMODEL_ID)//'_para_def.nc' #endif - NUMPSET=1 ! only the default parameter set is run - ALLOCATE(name_psets(NUMPSET)) - name_psets(1)='default_param_set' - ELSE IF(fuse_mode == 'run_pre')THEN ! run FUSE with pre-defined parameter values - ! read file_para_list twice: - ! 1st pass: determine number of parameter set and allocate name_psets accordingly - ! 2st pass: save the names of parameter sets in name_psets - - do file_pass = 1, 2 - - NUMPSET=0 ! intialize counter - - OPEN(21,FILE=TRIM(file_para_list)) - DO ! loop through parameter files - - READ(21,*,IOSTAT=ERR) dummy_string - IF (ERR.NE.0) EXIT - NUMPSET=NUMPSET+1 ! increment counter - - if (file_pass.eq.2) THEN - name_psets(NUMPSET) = dummy_string ! save file names - ENDIF - - END DO ! looping through parameter files - - CLOSE(21) - - if(file_pass.eq.1) THEN - print *, 'NUMPSET=', NUMPSET, 'based on the number of lines in ', TRIM(file_para_list) - ALLOCATE(name_psets(NUMPSET)) - END IF - end do - ! files to which model run and parameter set will be saved FNAME_NETCDF_RUNS = TRIM(OUTPUT_PATH)//TRIM(dom_id)//'_'//TRIM(FMODEL_ID)//'_runs_pre.nc' FNAME_NETCDF_PARA = TRIM(OUTPUT_PATH)//TRIM(dom_id)//'_'//TRIM(FMODEL_ID)//'_para_pre_out.nc' + NUMPSET=1 ! only the one "desired" parameter set is run + ELSE IF(fuse_mode == 'calib_sce')THEN ! calibrate FUSE using SCE ! files to which model run and parameter set will be saved @@ -417,22 +392,15 @@ PROGRAM DISTRIBUTED_DRIVER OUTPUT_FLAG=.TRUE. - do IPSET = 1, NUMPSET - - FNAME_NETCDF_PARA_PRE=TRIM(OUTPUT_PATH)//name_psets(IPSET) - PRINT *, 'Loading parameter set ',IPSET,':' - - ! load specific parameter set - ! 2nd argument is 1 because first (and only) parameter set should be loaded - CALL GET_PRE_PARAM(FNAME_NETCDF_PARA_PRE,1,ONEMOD,NUMPAR,APAR) - - print *, 'Running FUSE with pre-defined parameter set' - CALL FUSE_RMSE(APAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET) - print *, 'Done running FUSE with pre-defined parameter set' + FNAME_NETCDF_PARA_PRE=TRIM(OUTPUT_PATH)//file_param + PRINT *, 'Loading parameter set ',IPSET,':' - end do + ! load specific parameter set + CALL GET_PRE_PARAM(FNAME_NETCDF_PARA_PRE,IPSET,ONEMOD,NUMPAR,APAR) - DEALLOCATE(name_psets) + print *, 'Running FUSE with pre-defined parameter set' + CALL FUSE_RMSE(APAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,1) ! last argument IPSET=1 + print *, 'Done running FUSE with pre-defined parameter set' ELSE IF(fuse_mode == 'calib_sce')THEN ! calibrate FUSE using SCE From c00d202650903abe75e1a418ef0b92c53360e30b Mon Sep 17 00:00:00 2001 From: Martyn Clark Date: Mon, 1 Dec 2025 09:21:48 -0700 Subject: [PATCH 6/9] updates to IE solve --- build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 | 12 +- build/FUSE_SRC/FUSE_ENGINE/logismooth.f90 | 22 -- build/FUSE_SRC/FUSE_ENGINE/q_misscell.f90 | 2 +- build/FUSE_SRC/dshare/globaldata.f90 | 3 +- build/FUSE_SRC/physics/fix_ovshoot.f90 | 30 +- build/FUSE_SRC/physics/implicit_solve.f90 | 407 ++++++++++++++------- build/FUSE_SRC/physics/mstate_rhs_diff.f90 | 19 + build/FUSE_SRC/physics/q_baseflow_diff.f90 | 2 + build/FUSE_SRC/physics/q_misscell_diff.f90 | 26 +- build/FUSE_SRC/physics/smoothers.f90 | 132 +++++++ build/Makefile | 10 +- 11 files changed, 492 insertions(+), 173 deletions(-) delete mode 100644 build/FUSE_SRC/FUSE_ENGINE/logismooth.f90 create mode 100644 build/FUSE_SRC/physics/smoothers.f90 diff --git a/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 b/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 index 76226d0..1aa8567 100644 --- a/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 +++ b/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 @@ -103,7 +103,6 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG ! differentiable model type(parent) :: fuseStruct ! parent fuse data structure - ! --------------------------------------------------------------------------------------- ! allocate state vectors ALLOCATE(STATE0(NSTATE),STATE1(NSTATE),STAT=IERR) @@ -277,9 +276,16 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG ! populate parent fuse structure call get_parent(fuseStruct) - ! solve differentiable ODEs - call implicit_solve(fuseStruct, state0, state1, nState) + call implicit_solve(fuseStruct, state0, state1, nState, ierr, cmessage) + if(ierr/=0)then + print*, trim(cmessage) + print*, 'state0 = ', state0 + call implicit_solve(fuseStruct, state0, state1, nState, ierr, cmessage, isVerbose=.true.) + stop 1 + endif + + !print*, state1 !if(ITIM_IN > sim_beg+100) stop diff --git a/build/FUSE_SRC/FUSE_ENGINE/logismooth.f90 b/build/FUSE_SRC/FUSE_ENGINE/logismooth.f90 deleted file mode 100644 index b149654..0000000 --- a/build/FUSE_SRC/FUSE_ENGINE/logismooth.f90 +++ /dev/null @@ -1,22 +0,0 @@ -PURE FUNCTION LOGISMOOTH(STATE,STATE_MAX,PSMOOTH) -! --------------------------------------------------------------------------------------- -! Creator: -! -------- -! Martyn Clark, 2007 -! --------------------------------------------------------------------------------------- -! Purpose: -! -------- -! Uses a logistic function to smooth the threshold at the top of a bucket -! --------------------------------------------------------------------------------------- -USE nrtype -IMPLICIT NONE -REAL(SP), INTENT(IN) :: STATE ! model state -REAL(SP), INTENT(IN) :: STATE_MAX ! maximum model state -REAL(SP), INTENT(IN) :: PSMOOTH ! smoothing parameter (fraction of state) -REAL(SP) :: ASMOOTH ! actual smoothing -REAL(SP) :: LOGISMOOTH ! FUNCTION name -! --------------------------------------------------------------------------------------- -ASMOOTH = PSMOOTH*STATE_MAX ! actual smoothing -LOGISMOOTH = 1._SP / ( 1._SP + EXP(-(STATE - (STATE_MAX - ASMOOTH*5._SP) ) / ASMOOTH) ) -! --------------------------------------------------------------------------------------- -END FUNCTION LOGISMOOTH diff --git a/build/FUSE_SRC/FUSE_ENGINE/q_misscell.f90 b/build/FUSE_SRC/FUSE_ENGINE/q_misscell.f90 index b84bad5..2dc7257 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/q_misscell.f90 +++ b/build/FUSE_SRC/FUSE_ENGINE/q_misscell.f90 @@ -25,8 +25,8 @@ SUBROUTINE Q_MISSCELL() USE multistate, ONLY: MSTATE,TSTATE ! model states USE multi_flux, ONLY: M_FLUX ! model fluxes USE model_numerix ! access model numerix decisions +USE smoothers, only: logismooth ! logistic smoothing function IMPLICIT NONE -REAL(SP) :: LOGISMOOTH ! FUNCTION logistic smoothing REAL(SP), PARAMETER :: PSMOOTH=0.01_SP ! smoothing parameter REAL(SP) :: W_FUNC ! result from LOGISMOOTH REAL(SP) :: DT ! current time step diff --git a/build/FUSE_SRC/dshare/globaldata.f90 b/build/FUSE_SRC/dshare/globaldata.f90 index d9edb1d..9d551d3 100644 --- a/build/FUSE_SRC/dshare/globaldata.f90 +++ b/build/FUSE_SRC/dshare/globaldata.f90 @@ -16,7 +16,8 @@ MODULE globaldata real(sp), parameter :: fracState0=0.25_sp ! print flag - logical(lgt) :: isPrint=.true. + logical(lgt), save :: isPrint=.true. + logical(lgt), save :: isDebug=.false. ! number of fuse evaluations integer(i4b), save :: nFUSE_eval diff --git a/build/FUSE_SRC/physics/fix_ovshoot.f90 b/build/FUSE_SRC/physics/fix_ovshoot.f90 index 9e51da3..16feb4b 100644 --- a/build/FUSE_SRC/physics/fix_ovshoot.f90 +++ b/build/FUSE_SRC/physics/fix_ovshoot.f90 @@ -9,6 +9,7 @@ module overshoot_module private public :: get_bounds public :: fix_ovshoot + public :: sigmoid contains @@ -28,10 +29,21 @@ pure real(sp) function softplus(x, alpha) result(y) end function softplus ! --------------------------------------------------------------------------------------- ! --------------------------------------------------------------------------------------- + ! Sigmoid + pure real(sp) function sigmoid(z) result(s) + real(sp), intent(in) :: z + if (z >= 0._sp) then + s = 1._sp / (1._sp + exp(-z)) + else + s = exp(z) / (1._sp + exp(z)) + end if + end function sigmoid + ! --------------------------------------------------------------------------------------- + ! --------------------------------------------------------------------------------------- ! --------------------------------------------------------------------------------------- ! --------------------------------------------------------------------------------------- - SUBROUTINE fix_ovshoot(X_TRY, lower, upper) + SUBROUTINE fix_ovshoot(X_TRY, lower, upper, dclamp) ! --------------------------------------------------------------------------------------- ! Creator: ! -------- @@ -45,13 +57,23 @@ SUBROUTINE fix_ovshoot(X_TRY, lower, upper) 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 + real(sp), dimension(:), intent(out) :: dclamp ! derivative ! 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) + + ! hard constraints + x_try(i) = max( min(x_try(i), upper(i)), lower(i) ) + dclamp(i) = 1._sp + + ! ! apply soft constraint to model states + ! x_try(i) = lower(i) + softplus(x_try(i)-lower(i), alpha) - softplus(x_try(i)-upper(i), alpha) + ! + ! ! compute derivative in clamp + ! dclamp(i) = sigmoid( (x_try(i) - lower(i)) * alpha ) - sigmoid( (x_try(i) - upper(i)) * alpha ) + end do ! looping through model state variables end subroutine fix_ovshoot diff --git a/build/FUSE_SRC/physics/implicit_solve.f90 b/build/FUSE_SRC/physics/implicit_solve.f90 index 6ee2325..39a30c7 100644 --- a/build/FUSE_SRC/physics/implicit_solve.f90 +++ b/build/FUSE_SRC/physics/implicit_solve.f90 @@ -11,55 +11,36 @@ module implicit_solve_module ! global data use model_defn, only: nState ! number of state variables use multiforce, only: dt => deltim ! time step + use globaldata, only: isDebug ! print flag 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) + function dx_dt(fuseStruct, x_try) result(g_x) use MOD_DERIVS_DIFF_module, only: MOD_DERIVS_DIFF ! compute dx/dt implicit none ! input + type(parent) , intent(inout) :: fuseStruct ! parent fuse data structure 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) + call XTRY_2_STR(x_try, fuseStruct%state1) ! run the fuse physics - call mod_derivs_diff(ctx) + call mod_derivs_diff(fuseStruct) ! extract dx_dt from fuse structure - call STR_2_XTRY(ctx%dx_dt, g_x) + call STR_2_XTRY(fuseStruct%dx_dt, g_x) ! track the total number of function calls NUM_FUNCS = NUM_FUNCS + 1 @@ -67,72 +48,123 @@ function dx_dt(x_try) result(g_x) end function dx_dt ! ----- calculate the Jacobian of g(x) ------------------------------------------------- - SUBROUTINE jac_flux(x,g_x,Jac) + SUBROUTINE jac_flux(fuseStruct, x_try, g_x, lower, upper, Jac) IMPLICIT NONE - REAL(SP), DIMENSION(:), INTENT(IN) :: g_x - REAL(SP), DIMENSION(:), INTENT(INOUT) :: x + ! input-output + type(parent) , intent(inout) :: fuseStruct ! parent fuse data structure + REAL(SP), DIMENSION(:), INTENT(IN) :: g_x, lower, upper + REAL(SP), DIMENSION(:), INTENT(IN) :: x_try REAL(SP), DIMENSION(:,:), INTENT(OUT) :: Jac - REAL(SP), PARAMETER :: EPS=-1.0e-4_sp ! NOTE force h to be negative + ! locals + type(parent) :: ctx_sav + real(sp), parameter :: eps_rel = 1e-4_sp + real(sp), parameter :: eps_abs = 1e-6_sp ! or smaller, but NOT 1e-9 scale + real(sp), parameter :: h_min = 1e-8_sp 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 + REAL(SP), DIMENSION(size(x_try)) :: x, xsav, g_ph + real(sp) :: h_try, h_act + + ! preliminaries + n = size(x) + ctx_sav = fuseStruct + x = x_try + xsav = x + + ! loop through columns do j=1,n - x(j)=xph(j) - Jac(:,j)=(dx_dt(x)-g_x(:))/h(j) - x(j)=xsav(j) - end do + + ! safety: save full vector and data structure + fuseStruct = ctx_sav + x=xsav + + ! propose one-sided step + h_try = -max(eps_rel*abs(xsav(j)), eps_abs) + + ! flip sign if necessary + if(xsav(j) + h_try < lower(j)) h_try = -h_try + + ! compute function from the perturbed vector + x(j) = xsav(j) + h_try + g_ph = dx_dt(fuseStruct, x) + h_act = x(j) - xsav(j) + + ! compute column in the Jacobian + Jac(:,j) = (g_ph - g_x) / h_act + + end do ! looping through Jacobian columns + NUM_JACOBIAN = NUM_JACOBIAN + 1 ! keep track of the number of iterations - call XTRY_2_STR(xsav, ctx%state1) ! restores consistency after finite differencing + fuseStruct = ctx_sav ! restores consistency after finite differencing end SUBROUTINE jac_flux ! ----- simple implicit solve for differentiable model -------------------------- - subroutine implicit_solve(fuseStruct, x0, x1, nx) + subroutine implicit_solve(fuseStruct, x0, x1, nx, ierr, message, isVerbose) 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 + 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 + 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 + ! error control + integer(i4b), intent(out) :: ierr ! error code + character(*), intent(out) :: message ! error message + logical(lgt), intent(in), optional :: isVerbose ! flag for printing (subroutine argument) + logical(lgt) :: isPrint ! flag for printing (local flag) + ! internal: newton iterations + real(sp) :: x_old(nx) ! old trial state vector + 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 + type(parent) :: ctx ! save the fuse structure + real(sp) :: xnorm ! norm used in maximum step + real(sp) :: dxnorm ! norm used to evaluate step size + real(sp) :: stpmax ! the maximum step + real(sp) :: dxScale ! used to scale dx if dxnorm > stpmax + real(sp) :: gpsi(nx) ! function gradient: func = 0.5*sum(res*res) + real(sp) :: slope ! direction of decrease + real(sp) :: lambda ! backtrack length multiplier (lambda*dx) + real(sp) :: alamin ! minimum lambda + real(sp) :: lam_i ! maximum lambda for the i-th state + real(sp) :: lam_max ! maximum lambda + real(sp) :: lower(nx) ! lower bound + real(sp) :: upper(nx) ! lower bound + real(sp) :: dclamp(nx) ! derivative in the clamp + real(sp) :: x_trial(nx) ! state vector for 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 + ! algorithmic control parameters (most passed through MODULE model_numerix) + REAL(SP), PARAMETER :: TOLMIN=1.0e-10_sp ! check for spurious minima + REAL(SP), PARAMETER :: STPMX=100.0_sp ! maximum step in lnsrch + real(sp), parameter :: shrink = 0.5_sp + real(sp), parameter :: dampen = 0.1_sp + real(sp), parameter :: phi_rel_tol = 1e-5_sp ! 0.001% + real(sp), parameter :: phi_abs_tol = 1e-6_sp + real(sp), parameter :: epsb = 1.e-10_sp ! small safety margin + integer(i4b), parameter :: ls_max = 5 + ! ----- procedure starts here -------------------------------------------------------------------- + ! initialize error control + ierr=0; message='implicit_solve/' ! check dimension size if (nx /= nState) stop "implicit_solve: nx /= nState" @@ -141,89 +173,216 @@ subroutine implicit_solve(fuseStruct, x0, x1, nx) NUM_FUNCS = 0 ! number of function calls NUM_JACOBIAN = 0 ! number of times Jacobian is calculated + ! get the flag for printing + isPrint = .false.; if (present(isVerbose)) isPrint = isVerbose + + ! save the fuse structure + ctx = fuseStruct + ! 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 + ! intialize state vector (and soft clamp) + x_try = x0 + x_old = x_try + dclamp = 1._sp + + ! fix overshoot (only if necessary) + if(any(x_try < lower) .or. any(x_try > upper)) & + call fix_ovshoot(x_try, lower, upper, dclamp) + + ! define maximum step + xnorm = sqrt( sum(x_try*x_try) ) + stpmax = STPMX * max( xnorm, real(nx, sp) ) + + ! initialize flags accepted = .false. converged = .false. - ! --- F(x) and objective phi = 0.5*||F||^2 - g_x = dx_dt(x_try) + if(isPrint) isDebug = .true. + + ! --- F(x) and objective phi + g_x = dx_dt(fuseStruct, x_try) res = x_try - (x0 + g_x*dt) - phi = 0.5_sp * sum(res*res) + phi = 0.5_sp * dot_product(res, res) + + if(isPrint) isDebug = .false. ! iterate do it = 1, maxit - if (sqrt(2.0_sp*phi) < ERR_ITER_FUNC) then + ! save x + x_old = x_try + + if(isPrint) print*, '***** start of iteration *****' + + ! check + if(isPrint)then + print*, 'x_try = ', x_try + print*, 'g_x = ', g_x + print*, 'res = ', res + print*, 'phi = ', phi + print*, 'dclamp = ', dclamp + if(it > 10) stop 1 + endif + + if (phi < ERR_ITER_FUNC) then converged = .true. exit ! exit iteration loop end if + if(isPrint) print*, 'x_try 0 = ', x_try + ! --- 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 + call jac_flux(fuseStruct, x_try, g_x, lower, upper, Jg) + do i=1,nx + Jac(:,i) = -dt*Jg(:,i) !* dclamp(i) ! multiply dt and clamp derivative + Jac(i,i) = Jac(i,i) + 1.0_sp + end do + + if(isPrint) print*, 'x_try 1 = ', x_try + + ! --- function gradient: before Jac is modified in ludcmp + gpsi = matmul(transpose(Jac), res) ! assumes func = 0.5_sp * sum(res*res) ! --- 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 + + if(isPrint) print*, 'x_try 2 = ', x_try + + if(isPrint)then + print*, 'dx = ', dx + print*, 'Jg = ', Jg + endif - ! initialize flag to check if line search is accepted - accepted = .false. + ! --- Modify dx + + ! modify dx if norm > stpmax + dxnorm = sqrt( sum(dx*dx) ) + if (dxnorm > stpmax) then + dxScale = stpmax / dxnorm + dx = dxScale * dx + end if - ! ---- backtracking line search w/ overshoot reject ---- - lambda = 1.0_sp + ! implement active-set methods + do i=1,nx + if (x_try(i) <= lower(i)+epsb .and. dx(i) < 0._sp) dx(i)=0._sp + if (x_try(i) >= upper(i)-epsb .and. dx(i) > 0._sp) dx(i)=0._sp + end do + + ! modify dx if Newton step not descending for psi + slope = dot_product(gpsi, dx) + if (slope >= 0._sp) dx = -gpsi ! fallback + + if(isPrint) print*, 'x_try 3 = ', x_try + ! ---- backtracking line search -------------- + + ! save the fuse structure to re-use in subsequent linesearch calls + ctx = fuseStruct ! <--- snapshot *now*, for this Newton step + + ! line search control + accepted = .false. ! flag to check if line search is accepted + alamin = ERR_ITER_DX / maxval( abs(dx) / max(abs(x_try), 1.0_sp) ) + + ! compute maximum lambda + lam_max = 1.0_sp + ! do i=1,nx + ! if (dx(i) > 0._sp) then + ! lam_i = (upper(i) - x_try(i) - epsb) / dx(i) + ! if (lam_i < lam_max) lam_max = max(0._sp, lam_i) + ! else if (dx(i) < 0._sp) then + ! lam_i = (lower(i) - x_try(i) + epsb) / dx(i) ! dx<0 so this is positive + ! if (lam_i < lam_max) lam_max = max(0._sp, lam_i) + ! end if + ! end do + + ! check + if(isPrint)then + print*, 'alamin = ', alamin + print*, 'slope = ', slope + print*, 'gpsi = ', gpsi + endif + + if(isPrint) isDebug = .true. + + lambda = lam_max 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 + + if(isPrint)then + print*, '***** new linesearch *****', ls_it + print*, 'dx = ', dx + endif + + ! update x + x_trial = x_try + lambda*dx + + if(isPrint)then + print*, 'x_try = ', x_try + print*, 'x_trial = ', x_trial + print *, "delta = ", x_trial - x_try + print *, "lambda*dx = ", lambda*dx + endif + + ! exit if violate bounds: line search direction is not valid + if(any(x_trial < lower) .or. any(x_trial > upper))then + accepted = .false.; exit + !lambda = lambda * shrink + !cycle + endif + + ! compute function and function eval + fuseStruct = ctx + g_trial = dx_dt(fuseStruct, x_trial) + res_trial = x_trial - (x0 + dt*g_trial) + phi_new = 0.5_sp * dot_product(res_trial, res_trial) + + if(isPrint)then + print*, 'ls_it, lambda, phi, phi_new', ls_it, lambda, phi, phi_new + print*, 'phi, phi_new, slope=', phi, phi_new, slope + print*, 'x_trial = ', x_trial + print*, 'g_trial = ', g_trial + print*, 'res _trial= ', res_trial + endif + + if (phi_new <= phi + phi_abs_tol) then + accepted = .true.; exit + endif + + ! update lambda lambda = lambda * shrink + if (lambda < alamin) exit ! give up shrinking + end do ! line search + if(isPrint) isDebug = .false. + 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) + ! ----- fallback: try a small step along the direction of steepest descent ----- + !dx = -gpsi ! use steepest descent + x_trial = x_try + dampen*dx + if(any(x_trial < lower) .or. any(x_trial > upper)) & + call fix_ovshoot(x_trial, lower, upper, dclamp) ! 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) + x_try = x_trial + fuseStruct = ctx + g_x = dx_dt(fuseStruct, x_try) + res = x_try - (x0 + g_x*dt) + phi = 0.5_sp * dot_product(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 + if (maxval( abs(x_try - x_old) / max(abs(x_try), 1._sp) ) < ERR_ITER_DX) then converged = .true. exit ! exit iteration loop end if @@ -233,11 +392,11 @@ subroutine implicit_solve(fuseStruct, x0, x1, nx) ! 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" + if( .not. converged)then + message=trim(message)//"failed to converge" + ierr=10; return + endif end subroutine implicit_solve diff --git a/build/FUSE_SRC/physics/mstate_rhs_diff.f90 b/build/FUSE_SRC/physics/mstate_rhs_diff.f90 index 1ab0107..b34eb90 100644 --- a/build/FUSE_SRC/physics/mstate_rhs_diff.f90 +++ b/build/FUSE_SRC/physics/mstate_rhs_diff.f90 @@ -1,5 +1,7 @@ module MSTATE_RHS_DIFF_module + use globaldata, only: isDebug ! print flag + implicit none private @@ -32,6 +34,8 @@ SUBROUTINE MSTATE_RHS_DIFF(fuseStruct) DX_DT => fuseStruct%dx_dt & ! time derivative in states ) ! (associate) + if(isDebug) print*, 'M_FLUX%QPERC_12 = ', M_FLUX%QPERC_12 + ! --------------------------------------------------------------------------------------- ! (1) COMPUTE TIME DERIVATIVES FOR STATES IN THE UPPER LAYER ! --------------------------------------------------------------------------------------- @@ -71,6 +75,21 @@ SUBROUTINE MSTATE_RHS_DIFF(fuseStruct) END SELECT ! --------------------------------------------------------------------------------------- + if(isDebug) print*, '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 = ', & + 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 + + if(isDebug) print*, 'DX_DT%WATR_2, M_FLUX%QPERC_12, M_FLUX%EVAP_2, M_FLUX%QBASE_2, M_FLUX%OFLOW_2 = ', & + DX_DT%WATR_2, M_FLUX%QPERC_12, M_FLUX%EVAP_2, M_FLUX%QBASE_2, M_FLUX%OFLOW_2 + + ! if(isDebug) print*, 'DX_DT%TENS_1, M_FLUX%EFF_PPT, M_FLUX%QSURF, M_FLUX%EVAP_1, M_FLUX%TENS2FREE_1 = ', & + ! DX_DT%TENS_1, M_FLUX%EFF_PPT, M_FLUX%QSURF, M_FLUX%EVAP_1, M_FLUX%TENS2FREE_1 + ! + ! if(isDebug) print*, 'DX_DT%TENS_2, M_FLUX%QPERC_12*(1._SP-MPARAM%PERCFRAC), M_FLUX%EVAP_2, M_FLUX%TENS2FREE_2 = ', & + ! DX_DT%TENS_2, M_FLUX%QPERC_12*(1._SP-MPARAM%PERCFRAC), M_FLUX%EVAP_2, M_FLUX%TENS2FREE_2 + ! + ! if(isDebug) print*, '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 = ', & + ! 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 + end associate ! end association with variables in the data structures END SUBROUTINE MSTATE_RHS_DIFF diff --git a/build/FUSE_SRC/physics/q_baseflow_diff.f90 b/build/FUSE_SRC/physics/q_baseflow_diff.f90 index 29386f4..d63adb6 100644 --- a/build/FUSE_SRC/physics/q_baseflow_diff.f90 +++ b/build/FUSE_SRC/physics/q_baseflow_diff.f90 @@ -1,5 +1,7 @@ module Q_BASEFLOW_DIFF_module + USE globaldata, only: isDebug + implicit none private diff --git a/build/FUSE_SRC/physics/q_misscell_diff.f90 b/build/FUSE_SRC/physics/q_misscell_diff.f90 index 1801c7b..c347455 100644 --- a/build/FUSE_SRC/physics/q_misscell_diff.f90 +++ b/build/FUSE_SRC/physics/q_misscell_diff.f90 @@ -27,13 +27,13 @@ SUBROUTINE Q_MISSCELL_DIFF(fuseStruct) USE data_types, only: parent ! fuse parent data type USE model_defn ! model definition structure USE model_defnames + USE smoothers, only: smoother ! smoothing function 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 + REAL(SP), PARAMETER :: PSMOOTH=0.05_SP ! smoothing parameter + REAL(SP) :: W_FUNC ! result from smoother ! ------------------------------------------------------------------------------------------------- ! associate variables with elements of data structure associate(& @@ -48,29 +48,29 @@ SUBROUTINE Q_MISSCELL_DIFF(fuseStruct) 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) + W_FUNC = SMOOTHER(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) + W_FUNC = SMOOTHER(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) + W_FUNC = SMOOTHER(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) + W_FUNC = SMOOTHER(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) + W_FUNC = SMOOTHER(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) + W_FUNC = SMOOTHER(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" @@ -81,13 +81,13 @@ SUBROUTINE Q_MISSCELL_DIFF(fuseStruct) 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) + W_FUNC = SMOOTHER(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) + W_FUNC = SMOOTHER(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) + W_FUNC = SMOOTHER(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 @@ -97,7 +97,7 @@ SUBROUTINE Q_MISSCELL_DIFF(fuseStruct) 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) + W_FUNC = SMOOTHER(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 diff --git a/build/FUSE_SRC/physics/smoothers.f90 b/build/FUSE_SRC/physics/smoothers.f90 new file mode 100644 index 0000000..71f5277 --- /dev/null +++ b/build/FUSE_SRC/physics/smoothers.f90 @@ -0,0 +1,132 @@ +module smoothers + + implicit none + + private + public:: LOGISMOOTH + public:: smoother + +contains + + ! --------------------------------------------------------------------------------------- + ! --------------------------------------------------------------------------------------- + + PURE FUNCTION smoother(STATE,STATE_MAX,PSMOOTH) result(w_func) + ! --------------------------------------------------------------------------------------- + ! Creator: + ! -------- + ! Martyn Clark, 2025 + ! --------------------------------------------------------------------------------------- + ! Purpose: + ! -------- + ! Provides the option of different smoothers + ! --------------------------------------------------------------------------------------- + USE nrtype + IMPLICIT NONE + REAL(SP), INTENT(IN) :: STATE ! model state + REAL(SP), INTENT(IN) :: STATE_MAX ! maximum model state + REAL(SP), INTENT(IN) :: PSMOOTH ! smoothing parameter (fraction of state) + real(sp) :: w_func ! smoothed threshold + real(sp) :: delta ! scale factor + + ! logistic smoothing (original) + w_func = LOGISMOOTH(STATE,STATE_MAX,PSMOOTH) + + ! qintic smoother (plays better with Newton) + !delta = MAX(PSMOOTH*STATE_MAX, 1.0e-6_SP*STATE_MAX) + !w_func = SMOOTHSTEP5_W(STATE,STATE_MAX,delta) + + end function smoother + + ! --------------------------------------------------------------------------------------- + ! --------------------------------------------------------------------------------------- + + PURE FUNCTION LOGISMOOTH(STATE,STATE_MAX,PSMOOTH) + ! --------------------------------------------------------------------------------------- + ! Creator: + ! -------- + ! Martyn Clark, 2007 + ! --------------------------------------------------------------------------------------- + ! Purpose: + ! -------- + ! Uses a logistic function to smooth the threshold at the top of a bucket + ! --------------------------------------------------------------------------------------- + USE nrtype + IMPLICIT NONE + REAL(SP), INTENT(IN) :: STATE ! model state + REAL(SP), INTENT(IN) :: STATE_MAX ! maximum model state + REAL(SP), INTENT(IN) :: PSMOOTH ! smoothing parameter (fraction of state) + real(sp) :: arg ! clamp argument + REAL(SP) :: ASMOOTH ! actual smoothing + REAL(SP) :: LOGISMOOTH ! FUNCTION name + ! --------------------------------------------------------------------------------------- + ASMOOTH = PSMOOTH*STATE_MAX ! actual smoothing + arg = -(STATE - (STATE_MAX - 5*ASMOOTH))/ASMOOTH ! argument + !arg = max(min(arg, 50._SP), -50._SP) ! clamp + LOGISMOOTH = 1._SP / ( 1._SP + EXP(arg) ) + ! --------------------------------------------------------------------------------------- + END FUNCTION LOGISMOOTH + + ! --------------------------------------------------------------------------------------- + ! --------------------------------------------------------------------------------------- + + PURE FUNCTION SMOOTHSTEP5_W(STATE, STATE_MAX, DELTA) RESULT(W) + ! --------------------------------------------------------------------------------------- + ! Creator: + ! -------- + ! Martyn Clark, 2025 + ! --------------------------------------------------------------------------------------- + ! Purpose: + ! -------- + ! Uses a qintic function to smooth the threshold at the top of a bucket + ! --------------------------------------------------------------------------------------- + USE nrtype + IMPLICIT NONE + REAL(SP), INTENT(IN) :: STATE, STATE_MAX, DELTA + REAL(SP) :: W, x + + x = (STATE - (STATE_MAX - DELTA)) / DELTA + IF (x <= 0._SP) THEN + W = 0._SP + ELSEIF (x >= 1._SP) THEN + W = 1._SP + ELSE + W = x*x*x*(10._SP + x*(-15._SP + 6._SP*x)) + END IF + END FUNCTION + + ! --------------------------------------------------------------------------------------- + ! --------------------------------------------------------------------------------------- + + PURE FUNCTION SMOOTHSTEP5_DWDS(STATE, STATE_MAX, DELTA) RESULT(DWDS) + ! --------------------------------------------------------------------------------------- + ! Creator: + ! -------- + ! Martyn Clark, 2025 + ! --------------------------------------------------------------------------------------- + ! Purpose: + ! -------- + ! Compute the derivative of the qintic function + ! --------------------------------------------------------------------------------------- + USE nrtype + IMPLICIT NONE + REAL(SP), INTENT(IN) :: STATE, STATE_MAX, DELTA + REAL(SP) :: DWDS, x + + IF (DELTA <= 0._SP) THEN + DWDS = 0._SP + RETURN + END IF + + x = (STATE - (STATE_MAX - DELTA)) / DELTA + IF (x <= 0._SP .OR. x >= 1._SP) THEN + DWDS = 0._SP + ELSE + DWDS = (30._SP * x*x * (1._SP - x)*(1._SP - x)) / DELTA + END IF + END FUNCTION + + ! --------------------------------------------------------------------------------------- + ! --------------------------------------------------------------------------------------- + +end module smoothers diff --git a/build/Makefile b/build/Makefile index 9c2edfb..b0c04f6 100644 --- a/build/Makefile +++ b/build/Makefile @@ -128,6 +128,7 @@ NR_SUB = $(patsubst %, $(NUMREC_DIR)/%, $(FUSE_NR_SUB)) # FUSE physics FUSE_PHYSICS= \ + smoothers.f90 \ get_parent.f90 \ qsatexcess_diff.f90 \ evap_upper_diff.f90 \ @@ -154,7 +155,6 @@ FUSE_MODGUT=\ qpercolate.f90 \ q_baseflow.f90 \ q_misscell.f90 \ - logismooth.f90 \ mstate_eqn.f90 \ fix_states.f90 \ meanfluxes.f90 \ @@ -248,15 +248,15 @@ endif ifeq "$(FC)" "gfortran" FLAGS_NORMA = -O3 -ffree-line-length-none -fmax-errors=0 -cpp - FLAGS_DEBUG = -p -g -Wall -ffree-line-length-none -fmax-errors=0 -fbacktrace -fcheck=bounds -cpp - FLAGS_DEBUG = -p -Og -ffree-line-length-none -fmax-errors=0 -fcheck=bounds -cpp # without -warn all + #FLAGS_DEBUG = -Wall -ffree-line-length-none -fmax-errors=0 -fbacktrace -fcheck=bounds -cpp + FLAGS_DEBUG = -O0 -ffree-line-length-none -fmax-errors=0 -fcheck=bounds -cpp # without -warn all FLAGS_FIXED = -O2 -c -ffixed-form endif # select a set of flags -FLAGS = $(FLAGS_NORMA) -#FLAGS = $(FLAGS_DEBUG) +#FLAGS = $(FLAGS_NORMA) +FLAGS = $(FLAGS_DEBUG) # MPI: FUSE with MPI has been compiled successfully with mpif90 and mpiifort. # Note: override must be specifed to enable FC passed as argument to be overridden From ee8c182184fc6b72db03b5b0484f2c0b500e6979 Mon Sep 17 00:00:00 2001 From: Martyn Clark Date: Mon, 1 Dec 2025 12:00:18 -0700 Subject: [PATCH 7/9] fix uninitialized error code in get_mbands --- build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 | 1 + build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 | 3 -- build/FUSE_SRC/FUSE_ENGINE/get_mbands.f90 | 54 +++++++++++------------ 3 files changed, 28 insertions(+), 30 deletions(-) diff --git a/build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 b/build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 index 0c54a36..8f1d81b 100644 --- a/build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 +++ b/build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90 @@ -264,6 +264,7 @@ PROGRAM DISTRIBUTED_DRIVER ! get elevation band info, in particular N_BANDS CALL GET_MBANDS_INFO(ELEV_BANDS_NC,err,message) ! read band data from NetCDF file +if(err/=0)then; write(*,*) trim(message); stop; endif ! allocate space for elevation bands allocate(MBANDS_VAR_4d(nspat1,nspat2,N_BANDS,numtim_sub+1),stat=err) diff --git a/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 b/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 index 1aa8567..e1d30e2 100644 --- a/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 +++ b/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 @@ -268,9 +268,6 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG 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) diff --git a/build/FUSE_SRC/FUSE_ENGINE/get_mbands.f90 b/build/FUSE_SRC/FUSE_ENGINE/get_mbands.f90 index efd006c..6ec4fcb 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/get_mbands.f90 +++ b/build/FUSE_SRC/FUSE_ENGINE/get_mbands.f90 @@ -161,7 +161,7 @@ SUBROUTINE GET_MBANDS_INFO(ELEV_BANDS_NC,err,message) character(*), intent(out) :: message ! internal integer(i4b),parameter::lenPath=1024 ! DK/2008/10/21: allows longer file paths -INTEGER(I4B),DIMENSION(2) :: IERR ! error codes +INTEGER(I4B) :: IERR ! error code INTEGER(I4B) :: IUNIT ! input file unit CHARACTER(LEN=lenPath) :: CFILE ! name of control file CHARACTER(LEN=lenPath) :: BFILE ! name of band file @@ -181,21 +181,21 @@ SUBROUTINE GET_MBANDS_INFO(ELEV_BANDS_NC,err,message) ! internal: NetCDF read integer(i4b) :: ivarid_af,ivarid_me ! NetCDF variable ID for area_frac and mean_area integer(i4b),parameter :: ndims=3 ! number of dimensions for frac_area -integer(i4b) :: dimid_eb ! ID elevation bands +integer(i4b) :: dimid_eb ! ID elevation bands integer(i4b) :: iDimID ! dimension ID integer(i4b) :: dimLen ! dimension length ! --------------------------------------------------------------------------------------- ! read in NetCDF file defining the elevation bands -err=0 +err=0; ierr=0 CFILE = TRIM(INPUT_PATH)//ELEV_BANDS_NC ! control file info shared in MODULE directory print *, 'Loading elevation bands from:',TRIM(CFILE) INQUIRE(FILE=CFILE,EXIST=LEXIST) ! check that control file exists IF (.NOT.LEXIST) THEN - print *, 'f-GET_MBANDS_GRID/NetCDF file ',TRIM(CFILE),' for elevation bands does not exist ' - STOP + print *, 'f-GET_MBANDS_GRID/NetCDF file ',TRIM(CFILE),' for elevation bands does not exist ' + STOP ENDIF !open netcdf file @@ -221,14 +221,14 @@ SUBROUTINE GET_MBANDS_INFO(ELEV_BANDS_NC,err,message) if(err/=0)then; message=trim(message)//trim(nf90_strerror(err)); return; endif ! allocate 1 data stucture -ALLOCATE(MBANDS(N_BANDS),STAT=IERR(1)) +ALLOCATE(MBANDS(N_BANDS),STAT=IERR) ! allocate data structures ALLOCATE(Z_FORCING_grid(nspat1,nspat2),MBANDS_INFO_3d(nspat1,nspat2,n_bands),& - AF_TEMP(nspat1,nspat2,n_bands),ME_TEMP(nspat1,nspat2,n_bands),& - elev_mask(nspat1,nspat2),STAT=IERR(1)) + AF_TEMP(nspat1,nspat2,n_bands),ME_TEMP(nspat1,nspat2,n_bands),& + elev_mask(nspat1,nspat2),STAT=IERR) -IF (ANY(IERR.NE.0)) THEN +IF (IERR.NE.0) THEN message="f-GET_MBANDS/problem allocating elevation band data structures" err=100; return ENDIF @@ -243,26 +243,26 @@ SUBROUTINE GET_MBANDS_INFO(ELEV_BANDS_NC,err,message) ! populate MBANDS_INFO_3d, Z_FORCING_grid and elev_mask DO iSpat2=1,nSpat2 - DO iSpat1=1,nSpat1 + DO iSpat1=1,nSpat1 - MBANDS_INFO_3d(iSpat1,iSpat2,:)%Z_MID = me_TEMP(iSpat1,iSpat2,:) - MBANDS_INFO_3d(iSpat1,iSpat2,:)%AF = af_TEMP(iSpat1,iSpat2,:) - Z_FORCING_grid(iSpat1,iSpat2) = sum(me_TEMP(iSpat1,iSpat2,:)*af_TEMP(iSpat1,iSpat2,:)) ! estimate mean elevation of forcing using weighted mean of EB elevation - elev_mask(iSpat1,iSpat2) = me_TEMP(iSpat1,iSpat2,1) .EQ. NA_VALUE_SP ! if mean elevation first band is NA_VALUE, mask this grid cell + MBANDS_INFO_3d(iSpat1,iSpat2,:)%Z_MID = me_TEMP(iSpat1,iSpat2,:) + MBANDS_INFO_3d(iSpat1,iSpat2,:)%AF = af_TEMP(iSpat1,iSpat2,:) + Z_FORCING_grid(iSpat1,iSpat2) = sum(me_TEMP(iSpat1,iSpat2,:)*af_TEMP(iSpat1,iSpat2,:)) ! estimate mean elevation of forcing using weighted mean of EB elevation + elev_mask(iSpat1,iSpat2) = me_TEMP(iSpat1,iSpat2,1) .EQ. NA_VALUE_SP ! if mean elevation first band is NA_VALUE, mask this grid cell + + if(.NOT.elev_mask(iSpat1,iSpat2)) THEN ! only check area fraction sum to 1 if not NA_VALUE + + if (abs(sum(MBANDS_INFO_3d(iSpat1,iSpat2,:)%AF)-1).GT.1E-2) then ! check that area fraction sum to 1 + + print *, "The area fraction of all the elevation bands do not add up to 1" + !print *, 'Difference with 1 = ', abs(sum(MBANDS_INFO_3d(iSpat1,iSpat2,:)%AF)-1) + print *, 'AF', MBANDS_INFO_3d(iSpat1,iSpat2,:)%AF + stop + + end if + end if - if(.NOT.elev_mask(iSpat1,iSpat2)) THEN ! only check area fraction sum to 1 if not NA_VALUE - - if (abs(sum(MBANDS_INFO_3d(iSpat1,iSpat2,:)%AF)-1).GT.1E-2) then ! check that area fraction sum to 1 - - print *, "The area fraction of all the elevation bands do not add up to 1" - !print *, 'Difference with 1 = ', abs(sum(MBANDS_INFO_3d(iSpat1,iSpat2,:)%AF)-1) - print *, 'AF', MBANDS_INFO_3d(iSpat1,iSpat2,:)%AF - stop - - end if - end if - - END DO + END DO END DO err = nf90_close(ncid_eb) From 40f4015bc43505817eb6002922f8aca37b0d3d20 Mon Sep 17 00:00:00 2001 From: Martyn Clark Date: Tue, 2 Dec 2025 05:07:04 -0700 Subject: [PATCH 8/9] further improvements to implicit solve --- build/FUSE_SRC/FUSE_ENGINE/fix_states.f90 | 4 +- build/FUSE_SRC/physics/conserve_clamp.f90 | 303 ++++++++++++++++++++++ build/FUSE_SRC/physics/implicit_solve.f90 | 114 ++++---- build/Makefile | 1 + 4 files changed, 366 insertions(+), 56 deletions(-) create mode 100644 build/FUSE_SRC/physics/conserve_clamp.f90 diff --git a/build/FUSE_SRC/FUSE_ENGINE/fix_states.f90 b/build/FUSE_SRC/FUSE_ENGINE/fix_states.f90 index 01a3bb2..5ef6ce0 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/fix_states.f90 +++ b/build/FUSE_SRC/FUSE_ENGINE/fix_states.f90 @@ -255,14 +255,14 @@ SUBROUTINE FIX_STATES(DT,ERROR_FLAG) ! ------------------------------------------------------------------------------------- CASE (iopt_WATR_2) IF (ESTATE%WATR_2.LT.XMIN*MPARAM%MAXWATR_2) THEN ! too much drainage - ERROR_LOSS = (ESTATE%WATR_2 - XMIN*MPARAM%MAXWATR_1)/DT ! error (L/T) + ERROR_LOSS = (ESTATE%WATR_2 - XMIN*MPARAM%MAXWATR_2)/DT ! error (L/T) TOTAL_LOSS = M_FLUX%EVAP_2 + M_FLUX%QBASE_2 M_FLUX%EVAP_2 = M_FLUX%EVAP_2 + (M_FLUX%EVAP_2 /TOTAL_LOSS)*ERROR_LOSS M_FLUX%QBASE_2 = M_FLUX%QBASE_2 + (M_FLUX%QBASE_2/TOTAL_LOSS)*ERROR_LOSS ESTATE%WATR_2 = XMIN*MPARAM%MAXWATR_2 ! (correct state) ERROR_FLAG = .TRUE. ENDIF - IF (ESTATE%FREE_2B.GT.DPARAM%MAXFREE_2B) THEN + IF (ESTATE%WATR_2.GT.MPARAM%MAXWATR_2) THEN ERROR_LOSS = (ESTATE%WATR_2 - MPARAM%MAXWATR_2)/DT M_FLUX%OFLOW_2 = M_FLUX%OFLOW_2 + ERROR_LOSS ESTATE%WATR_2 = MPARAM%MAXWATR_2 ! (correct state) diff --git a/build/FUSE_SRC/physics/conserve_clamp.f90 b/build/FUSE_SRC/physics/conserve_clamp.f90 new file mode 100644 index 0000000..03ec0d1 --- /dev/null +++ b/build/FUSE_SRC/physics/conserve_clamp.f90 @@ -0,0 +1,303 @@ +module conserve_clamp_module + + ! data types + use nrtype ! variable types, etc. + use data_types, only: parent ! parent fuse data structure + USE model_defn ! model definition structure + USE model_defnames + USE model_numerix + + implicit none + + private + public :: conserve_clamp + + contains + + SUBROUTINE conserve_clamp(fuseStruct,DT,ERROR_FLAG) + ! --------------------------------------------------------------------------------------- + ! Creator: + ! -------- + ! Martyn Clark, 2009 + ! Modified by Martyn Clark to pass fuse parent data structure, 12/2025 + ! --------------------------------------------------------------------------------------- + ! Purpose: + ! -------- + ! Ensure states are within bounds, and disaggregate fluxes if necessary + ! - This routine handles the very rare case (less than one-in-a-million) where + ! the implicit Euler solver fails to converge + ! --------------------------------------------------------------------------------------- + IMPLICIT NONE + ! input/output + type(parent) , intent(inout) :: fuseStruct ! parent fuse data structure + REAL(SP), INTENT(IN) :: DT ! time step + LOGICAL(LGT), INTENT(OUT) :: ERROR_FLAG ! .TRUE. if extrapolation error + ! internal + REAL(SP) :: XMIN ! very small number + INTEGER(I4B) :: ISTT ! loop through model states + REAL(SP) :: ERROR_LOSS ! error (L/T) + REAL(SP) :: TOTAL_LOSS ! total loss (L/T) + ! --------------------------------------------------------------------------------------- + ! associate variables with elements of data structure + associate(& + M_FLUX => fuseStruct%flux , & ! fluxes + BSTATE => fuseStruct%state0 , & ! state variables (start of step) + ESTATE => fuseStruct%state1 , & ! trial state variables (end of step) + MPARAM => fuseStruct%param_adjust , & ! adjustable model parameters + DPARAM => fuseStruct%param_derive & ! derived model parameters + ) ! (associate) + ! --------------------------------------------------------------------------------------- + ERROR_FLAG=.FALSE. ! initialize with no extrapolation error + ! --------------------------------------------------------------------------------------- + XMIN = FRACSTATE_MIN ! used to avoid zero derivatives + ! --------------------------------------------------------------------------------------- + DO ISTT=1,NSTATE + if (M_FLUX%QSURF.LT.0._sp) print *, 'start ', desc_int2str(cstate(istt)%isname), M_FLUX%QSURF + ERROR_LOSS = 0._SP ! initialize state error + SELECT CASE(CSTATE(ISTT)%iSNAME) + ! --------------------------------------------------------------------------------------- + ! (1) FIX STATES IN THE UPPER LAYER + ! ------------------------------------------------------------------------------------- + CASE (iopt_TENS1A) + IF (ESTATE%TENS_1A.LT.XMIN*DPARAM%MAXTENS_1A) THEN ! too much drainage + ERROR_LOSS = (ESTATE%TENS_1A - XMIN*DPARAM%MAXTENS_1A)/DT ! error (L/T) + TOTAL_LOSS = M_FLUX%QSURF + M_FLUX%EVAP_1A ! total loss (L/T) + M_FLUX%QSURF = M_FLUX%QSURF + (M_FLUX%QSURF /TOTAL_LOSS)*ERROR_LOSS + M_FLUX%EVAP_1A = M_FLUX%EVAP_1A + (M_FLUX%EVAP_1A/TOTAL_LOSS)*ERROR_LOSS + ESTATE%TENS_1A = XMIN*DPARAM%MAXTENS_1A ! (correct state) + ERROR_FLAG = .TRUE. + ENDIF + IF (ESTATE%TENS_1A.GT.DPARAM%MAXTENS_1A) THEN ! too much input + ERROR_LOSS = (ESTATE%TENS_1A - DPARAM%MAXTENS_1A)/DT + M_FLUX%RCHR2EXCS = M_FLUX%RCHR2EXCS + ERROR_LOSS + ESTATE%TENS_1A = DPARAM%MAXTENS_1A ! (correct state) + ESTATE%TENS_1B = BSTATE%TENS_1B + & ! (correct subsequent states) + (M_FLUX%RCHR2EXCS - M_FLUX%EVAP_1B - M_FLUX%TENS2FREE_1)*DT + ERROR_FLAG = .TRUE. + ENDIF + M_FLUX%ERR_TENS_1A = ERROR_LOSS + ! ------------------------------------------------------------------------------------- + CASE (iopt_TENS1B) + IF (ESTATE%TENS_1B.LT.XMIN*DPARAM%MAXTENS_1B) THEN ! too much drainage + ERROR_LOSS = (ESTATE%TENS_1B - XMIN*DPARAM%MAXTENS_1B)/DT + M_FLUX%EVAP_1B = M_FLUX%EVAP_1B + ERROR_LOSS + ESTATE%TENS_1B = XMIN*DPARAM%MAXTENS_1B ! (correct state) + ERROR_FLAG = .TRUE. + ENDIF + IF (ESTATE%TENS_1B.GT.DPARAM%MAXTENS_1B) THEN ! too much input + ERROR_LOSS = (ESTATE%TENS_1B - DPARAM%MAXTENS_1B)/DT + M_FLUX%TENS2FREE_1 = M_FLUX%TENS2FREE_1 + ERROR_LOSS + ESTATE%TENS_1B = DPARAM%MAXTENS_1B ! (correct state) + ESTATE%FREE_1 = BSTATE%FREE_1 + & ! (correct subsequent states) + (M_FLUX%TENS2FREE_1 - M_FLUX%QPERC_12 - M_FLUX%QINTF_1 - M_FLUX%OFLOW_1)*DT + ERROR_FLAG = .TRUE. + ENDIF + M_FLUX%ERR_TENS_1B = ERROR_LOSS + ! ------------------------------------------------------------------------------------- + CASE (iopt_TENS_1) + IF (ESTATE%TENS_1.LT.XMIN*DPARAM%MAXTENS_1) THEN ! too much drainage + ERROR_LOSS = (ESTATE%TENS_1 - XMIN*DPARAM%MAXTENS_1)/DT ! error (L/T) + TOTAL_LOSS = M_FLUX%QSURF + M_FLUX%EVAP_1 ! total loss (L/T) + M_FLUX%QSURF = M_FLUX%QSURF + (M_FLUX%QSURF /TOTAL_LOSS)*ERROR_LOSS + M_FLUX%EVAP_1 = M_FLUX%EVAP_1 + (M_FLUX%EVAP_1/TOTAL_LOSS)*ERROR_LOSS + ESTATE%TENS_1 = XMIN*DPARAM%MAXTENS_1 ! (correct state) + ERROR_FLAG = .TRUE. + ENDIF + IF (ESTATE%TENS_1.GT.DPARAM%MAXTENS_1) THEN ! too much input + ERROR_LOSS = (ESTATE%TENS_1 - DPARAM%MAXTENS_1)/DT + M_FLUX%TENS2FREE_1 = M_FLUX%TENS2FREE_1 + (ESTATE%TENS_1 - DPARAM%MAXTENS_1)/DT + ESTATE%TENS_1 = DPARAM%MAXTENS_1 ! (correct state) + ESTATE%FREE_1 = BSTATE%FREE_1 + & ! (correct subsequent states) + (M_FLUX%TENS2FREE_1 - M_FLUX%QPERC_12 - M_FLUX%QINTF_1 - M_FLUX%OFLOW_1)*DT + ERROR_FLAG = .TRUE. + ENDIF + M_FLUX%ERR_TENS_1 = ERROR_LOSS + ! ------------------------------------------------------------------------------------- + CASE (iopt_FREE_1) + IF (ESTATE%FREE_1.LT.XMIN*DPARAM%MAXFREE_1) THEN ! too much drainage + ERROR_LOSS = (ESTATE%FREE_1 - XMIN*DPARAM%MAXFREE_1)/DT ! error (L/T) + TOTAL_LOSS = M_FLUX%QPERC_12 + M_FLUX%QINTF_1 ! total loss (L/T) + M_FLUX%QPERC_12 = M_FLUX%QPERC_12 + (M_FLUX%QPERC_12/TOTAL_LOSS)*ERROR_LOSS + M_FLUX%QINTF_1 = M_FLUX%QINTF_1 + (M_FLUX%QINTF_1 /TOTAL_LOSS)*ERROR_LOSS + ESTATE%FREE_1 = XMIN*DPARAM%MAXFREE_1 ! (correct state) + ! correct subsequent states (deal appropriately with percolation) + ! NOTE: do this here because only necessary to make corrections if M_FLUX%QPERC_12 changes + SELECT CASE(SMODL%iARCH2) + CASE(iopt_tens2pll_2) ! tension reservoir plus two parallel tanks + ! fix overflow fluxes + M_FLUX%TENS2FREE_2 = MAX(0._SP, M_FLUX%QPERC_12*(1._SP-MPARAM%PERCFRAC) - (DPARAM%MAXTENS_2 - BSTATE%TENS_2 )/DT) + M_FLUX%OFLOW_2A = MAX(0._SP, (M_FLUX%QPERC_12*(MPARAM%PERCFRAC/2._SP) + M_FLUX%TENS2FREE_2/2._SP) & + - (DPARAM%MAXFREE_2A - BSTATE%FREE_2A)/DT) + M_FLUX%OFLOW_2B = MAX(0._SP, (M_FLUX%QPERC_12*(MPARAM%PERCFRAC/2._SP) + M_FLUX%TENS2FREE_2/2._SP) & + - (DPARAM%MAXFREE_2B - BSTATE%FREE_2B)/DT) + M_FLUX%OFLOW_2 = M_FLUX%OFLOW_2A + M_FLUX%OFLOW_2B + ! fix states + ESTATE%TENS_2 = BSTATE%TENS_2 + & + (M_FLUX%QPERC_12*(1._SP-MPARAM%PERCFRAC) - M_FLUX%EVAP_2 - M_FLUX%TENS2FREE_2)*DT + ESTATE%FREE_2A = BSTATE%FREE_2A + & + (M_FLUX%QPERC_12*(MPARAM%PERCFRAC/2._SP) + M_FLUX%TENS2FREE_2/2._SP - M_FLUX%QBASE_2A & + - M_FLUX%OFLOW_2A)*DT + ESTATE%FREE_2B = BSTATE%FREE_2B + & + (M_FLUX%QPERC_12*(MPARAM%PERCFRAC/2._SP) + M_FLUX%TENS2FREE_2/2._SP - M_FLUX%QBASE_2B & + - M_FLUX%OFLOW_2B)*DT + CASE(iopt_unlimfrc_2,iopt_unlimpow_2,iopt_fixedsiz_2) ! single state + ! NOTE: M_FLUX%OFLOW_2 and M_FLUX%EVAP_2 only calculated for 'fixedsiz_2' + ! fix overflow + IF (SMODL%iARCH2.EQ.iopt_fixedsiz_2) & + M_FLUX%OFLOW_2 = MAX(0._SP, M_FLUX%QPERC_12 - (MPARAM%MAXWATR_2 - BSTATE%WATR_2)/DT) + ! fix states + ESTATE%WATR_2 = BSTATE%WATR_2 + & + (M_FLUX%QPERC_12 - M_FLUX%EVAP_2 - M_FLUX%QBASE_2 - M_FLUX%OFLOW_2)*DT + CASE DEFAULT; stop ' SMODL%iARCH2 must be iopt_tens2pll_2, iopt_unlimfrc_2, iopt_unlimpow_2 or iopt_fixedsiz_2 ' + END SELECT ! deal with modified percolation of water to the lower layer + ERROR_FLAG = .TRUE. + ENDIF + IF (ESTATE%FREE_1.GT.DPARAM%MAXFREE_1) THEN ! too much input + ERROR_LOSS = (ESTATE%FREE_1 - DPARAM%MAXFREE_1)/DT + M_FLUX%OFLOW_1 = M_FLUX%OFLOW_1 + ERROR_LOSS + ESTATE%FREE_1 = DPARAM%MAXFREE_1 ! (correct state) + ERROR_FLAG = .TRUE. + ENDIF + M_FLUX%ERR_FREE_1 = ERROR_LOSS + ! ------------------------------------------------------------------------------------- + CASE (iopt_WATR_1) + IF (ESTATE%WATR_1.LT.XMIN*MPARAM%MAXWATR_1) THEN ! too much drainage + ERROR_LOSS = (ESTATE%WATR_1 - XMIN*MPARAM%MAXWATR_1)/DT ! error (L/T) + TOTAL_LOSS = M_FLUX%QSURF + M_FLUX%EVAP_1 + M_FLUX%QPERC_12 + M_FLUX%QINTF_1 + M_FLUX%QSURF = M_FLUX%QSURF + (M_FLUX%QSURF /TOTAL_LOSS)*ERROR_LOSS + M_FLUX%EVAP_1 = M_FLUX%EVAP_1 + (M_FLUX%EVAP_1 /TOTAL_LOSS)*ERROR_LOSS + M_FLUX%QINTF_1 = M_FLUX%QINTF_1 + (M_FLUX%QINTF_1 /TOTAL_LOSS)*ERROR_LOSS + M_FLUX%QPERC_12 = M_FLUX%QPERC_12 + (M_FLUX%QPERC_12/TOTAL_LOSS)*ERROR_LOSS + ESTATE%WATR_1 = XMIN*MPARAM%MAXWATR_1 ! (correct state) + ! correct subsequent states (deal appropriately with percolation) + ! NOTE: do this here because only necessary to make corrections if M_FLUX%QPERC_12 changes + SELECT CASE(SMODL%iARCH2) + CASE(iopt_tens2pll_2) ! tension reservoir plus two parallel tanks + ! fix overflow fluxes + M_FLUX%TENS2FREE_2 = MAX(0._SP, M_FLUX%QPERC_12*(1._SP-MPARAM%PERCFRAC) - (DPARAM%MAXTENS_2 - BSTATE%TENS_2 )/DT) + M_FLUX%OFLOW_2A = MAX(0._SP, (M_FLUX%QPERC_12*(MPARAM%PERCFRAC/2._SP) + M_FLUX%TENS2FREE_2/2._SP) & + - (DPARAM%MAXFREE_2A - BSTATE%FREE_2A)/DT) + M_FLUX%OFLOW_2B = MAX(0._SP, (M_FLUX%QPERC_12*(MPARAM%PERCFRAC/2._SP) + M_FLUX%TENS2FREE_2/2._SP) & + - (DPARAM%MAXFREE_2B - BSTATE%FREE_2B)/DT) + M_FLUX%OFLOW_2 = M_FLUX%OFLOW_2A + M_FLUX%OFLOW_2B + ! fix states + ESTATE%TENS_2 = BSTATE%TENS_2 + & + (M_FLUX%QPERC_12*(1._SP-MPARAM%PERCFRAC) - M_FLUX%EVAP_2 - M_FLUX%TENS2FREE_2)*DT + ESTATE%FREE_2A = BSTATE%FREE_2A + & + (M_FLUX%QPERC_12*(MPARAM%PERCFRAC/2._SP) + M_FLUX%TENS2FREE_2/2._SP - M_FLUX%QBASE_2A & + - M_FLUX%OFLOW_2A)*DT + ESTATE%FREE_2B = BSTATE%FREE_2B + & + (M_FLUX%QPERC_12*(MPARAM%PERCFRAC/2._SP) + M_FLUX%TENS2FREE_2/2._SP - M_FLUX%QBASE_2B & + - M_FLUX%OFLOW_2B)*DT + CASE(iopt_unlimfrc_2,iopt_unlimpow_2,iopt_fixedsiz_2) ! single state + ! NOTE: M_FLUX%OFLOW_2 and M_FLUX%EVAP_2 only calculated for 'fixedsiz_2' + ! fix overflow + IF (SMODL%iARCH2.EQ.iopt_fixedsiz_2) & + M_FLUX%OFLOW_2 = MAX(0._SP, M_FLUX%QPERC_12 - (MPARAM%MAXWATR_2 - BSTATE%WATR_2)/DT) + ! fix states + ESTATE%WATR_2 = BSTATE%WATR_2 + & + (M_FLUX%QPERC_12 - M_FLUX%EVAP_2 - M_FLUX%QBASE_2 - M_FLUX%OFLOW_2)*DT + CASE DEFAULT; stop ' SMODL%iARCH2 must be iopt_tens2pll_2, iopt_unlimfrc_2, iopt_unlimpow_2 or iopt_fixedsiz_2 ' + END SELECT ! deal with modified percolation of water to the lower layer + ERROR_FLAG = .TRUE. + ENDIF + IF (ESTATE%WATR_1.GT.MPARAM%MAXWATR_1) THEN ! too much input + ERROR_LOSS = (ESTATE%WATR_1 - MPARAM%MAXWATR_1)/DT + M_FLUX%OFLOW_1 = M_FLUX%OFLOW_1 + ERROR_LOSS + ESTATE%WATR_1 = MPARAM%MAXWATR_1 ! (correct state) + ERROR_FLAG = .TRUE. + ENDIF + M_FLUX%ERR_WATR_1 = ERROR_LOSS + ! ------------------------------------------------------------------------------------- + ! (2) FIX STATES IN THE LOWER LAYER + ! ------------------------------------------------------------------------------------- + CASE (iopt_TENS_2) + IF (ESTATE%TENS_2.LT.XMIN*DPARAM%MAXTENS_2) THEN ! too much drainage + ERROR_LOSS = (ESTATE%TENS_2 - XMIN*DPARAM%MAXTENS_2)/DT + M_FLUX%EVAP_2 = M_FLUX%EVAP_2 + ERROR_LOSS + ESTATE%TENS_2 = XMIN*DPARAM%MAXTENS_2 ! (correct state) + ERROR_FLAG = .TRUE. + ENDIF + IF (ESTATE%TENS_2.GT.DPARAM%MAXTENS_2) THEN ! too much input + ERROR_LOSS = (ESTATE%TENS_2 - DPARAM%MAXTENS_2)/DT + M_FLUX%TENS2FREE_2 = M_FLUX%TENS2FREE_2 + ERROR_LOSS + ESTATE%TENS_2 = DPARAM%MAXTENS_2 ! (correct state) + ! ** correct subsequent states (NOTE: 2 parallel tanks always coupled with a tension store) + ! fix overflow fluxes + M_FLUX%OFLOW_2A = MAX(0._SP, (M_FLUX%QPERC_12*(MPARAM%PERCFRAC/2._SP) + M_FLUX%TENS2FREE_2/2._SP) & + - (DPARAM%MAXFREE_2A - BSTATE%FREE_2A)/DT) + M_FLUX%OFLOW_2B = MAX(0._SP, (M_FLUX%QPERC_12*(MPARAM%PERCFRAC/2._SP) + M_FLUX%TENS2FREE_2/2._SP) & + - (DPARAM%MAXFREE_2B - BSTATE%FREE_2B)/DT) + M_FLUX%OFLOW_2 = M_FLUX%OFLOW_2A + M_FLUX%OFLOW_2B + ! fix states + ESTATE%FREE_2A = BSTATE%FREE_2A + & + (M_FLUX%QPERC_12*(MPARAM%PERCFRAC/2._SP) + M_FLUX%TENS2FREE_2/2._SP & + - M_FLUX%QBASE_2A - M_FLUX%OFLOW_2A)*DT + ESTATE%FREE_2B = BSTATE%FREE_2B + & + (M_FLUX%QPERC_12*(MPARAM%PERCFRAC/2._SP) + M_FLUX%TENS2FREE_2/2._SP & + - M_FLUX%QBASE_2B - M_FLUX%OFLOW_2B)*DT + ERROR_FLAG = .TRUE. + ENDIF + M_FLUX%ERR_TENS_2 = ERROR_LOSS + ! ------------------------------------------------------------------------------------- + CASE (iopt_FREE2A) + IF (ESTATE%FREE_2A.LT.XMIN*DPARAM%MAXFREE_2A) THEN ! too much drainage + ERROR_LOSS = (ESTATE%FREE_2A - XMIN*DPARAM%MAXFREE_2A)/DT + M_FLUX%QBASE_2A = M_FLUX%QBASE_2A + ERROR_LOSS + ESTATE%FREE_2A = XMIN*DPARAM%MAXFREE_2A ! (correct state) + ERROR_FLAG = .TRUE. + ENDIF + IF (ESTATE%FREE_2A.GT.DPARAM%MAXFREE_2A) THEN ! too much input + ERROR_LOSS = (ESTATE%FREE_2A - DPARAM%MAXFREE_2A)/DT + M_FLUX%OFLOW_2A = M_FLUX%OFLOW_2A + ERROR_LOSS + ESTATE%FREE_2A = DPARAM%MAXFREE_2A ! (correct state) + ERROR_FLAG = .TRUE. + ENDIF + M_FLUX%ERR_FREE_2A = ERROR_LOSS + ! ------------------------------------------------------------------------------------- + CASE (iopt_FREE2B) + IF (ESTATE%FREE_2B.LT.XMIN*DPARAM%MAXFREE_2B) THEN ! too much drainage + ERROR_LOSS = (ESTATE%FREE_2B - XMIN*DPARAM%MAXFREE_2B)/DT + M_FLUX%QBASE_2B = M_FLUX%QBASE_2B + ERROR_LOSS + ESTATE%FREE_2B = XMIN*DPARAM%MAXFREE_2B ! (correct state) + ERROR_FLAG = .TRUE. + ENDIF + IF (ESTATE%FREE_2B.GT.DPARAM%MAXFREE_2B) THEN ! too much input + ERROR_LOSS = (ESTATE%FREE_2B - DPARAM%MAXFREE_2B)/DT + M_FLUX%OFLOW_2B = M_FLUX%OFLOW_2B + ERROR_LOSS + ESTATE%FREE_2B = DPARAM%MAXFREE_2B ! (correct state) + ERROR_FLAG = .TRUE. + ENDIF + M_FLUX%ERR_FREE_2B = ERROR_LOSS + ! ------------------------------------------------------------------------------------- + CASE (iopt_WATR_2) + IF (ESTATE%WATR_2.LT.XMIN*MPARAM%MAXWATR_2) THEN ! too much drainage + ERROR_LOSS = (ESTATE%WATR_2 - XMIN*MPARAM%MAXWATR_2)/DT ! error (L/T) + TOTAL_LOSS = M_FLUX%EVAP_2 + M_FLUX%QBASE_2 + M_FLUX%EVAP_2 = M_FLUX%EVAP_2 + (M_FLUX%EVAP_2 /TOTAL_LOSS)*ERROR_LOSS + M_FLUX%QBASE_2 = M_FLUX%QBASE_2 + (M_FLUX%QBASE_2/TOTAL_LOSS)*ERROR_LOSS + ESTATE%WATR_2 = XMIN*MPARAM%MAXWATR_2 ! (correct state) + ERROR_FLAG = .TRUE. + ENDIF + IF (ESTATE%WATR_2.GT.MPARAM%MAXWATR_2) THEN + ERROR_LOSS = (ESTATE%WATR_2 - MPARAM%MAXWATR_2)/DT + M_FLUX%OFLOW_2 = M_FLUX%OFLOW_2 + ERROR_LOSS + ESTATE%WATR_2 = MPARAM%MAXWATR_2 ! (correct state) + ERROR_FLAG = .TRUE. + ENDIF + M_FLUX%ERR_WATR_2 = ERROR_LOSS + CASE DEFAULT; STOP ' cannot find state in fix_states() ' + END SELECT ! select state variable for processing + if (M_FLUX%QSURF.LT.0._sp) print *, 'end ', desc_int2str(cstate(istt)%isname), M_FLUX%QSURF + END DO ! loop through state variables + ! --------------------------------------------------------------------------------------- + ! compute derived fluxes, if necessary + IF (SMODL%iARCH2.EQ.iopt_tens2pll_2) THEN ! tension reservoir plus two parallel tanks + M_FLUX%QBASE_2 = M_FLUX%QBASE_2A + M_FLUX%QBASE_2B + M_FLUX%OFLOW_2 = M_FLUX%OFLOW_2A + M_FLUX%OFLOW_2B + ENDIF + ! --------------------------------------------------------------------------------------- + end associate ! end association with variables in the data structures + END SUBROUTINE conserve_clamp + +end module conserve_clamp_module diff --git a/build/FUSE_SRC/physics/implicit_solve.f90 b/build/FUSE_SRC/physics/implicit_solve.f90 index 39a30c7..db92b58 100644 --- a/build/FUSE_SRC/physics/implicit_solve.f90 +++ b/build/FUSE_SRC/physics/implicit_solve.f90 @@ -51,12 +51,12 @@ end function dx_dt SUBROUTINE jac_flux(fuseStruct, x_try, g_x, lower, upper, Jac) IMPLICIT NONE ! input-output - type(parent) , intent(inout) :: fuseStruct ! parent fuse data structure + type(parent) , intent(in) :: fuseStruct ! parent fuse data structure REAL(SP), DIMENSION(:), INTENT(IN) :: g_x, lower, upper REAL(SP), DIMENSION(:), INTENT(IN) :: x_try REAL(SP), DIMENSION(:,:), INTENT(OUT) :: Jac ! locals - type(parent) :: ctx_sav + type(parent) :: fuseStruct_local real(sp), parameter :: eps_rel = 1e-4_sp real(sp), parameter :: eps_abs = 1e-6_sp ! or smaller, but NOT 1e-9 scale real(sp), parameter :: h_min = 1e-8_sp @@ -66,18 +66,14 @@ SUBROUTINE jac_flux(fuseStruct, x_try, g_x, lower, upper, Jac) ! preliminaries n = size(x) - ctx_sav = fuseStruct + fuseStruct_local = fuseStruct x = x_try xsav = x ! loop through columns do j=1,n - ! safety: save full vector and data structure - fuseStruct = ctx_sav - x=xsav - - ! propose one-sided step + ! propose one-sided step (NOTE: negative) h_try = -max(eps_rel*abs(xsav(j)), eps_abs) ! flip sign if necessary @@ -85,16 +81,19 @@ SUBROUTINE jac_flux(fuseStruct, x_try, g_x, lower, upper, Jac) ! compute function from the perturbed vector x(j) = xsav(j) + h_try - g_ph = dx_dt(fuseStruct, x) + g_ph = dx_dt(fuseStruct_local, x) h_act = x(j) - xsav(j) ! compute column in the Jacobian Jac(:,j) = (g_ph - g_x) / h_act + ! safety: save full vector and data structure + fuseStruct_local = fuseStruct ! restores consistency after finite differencing + x = xsav + end do ! looping through Jacobian columns NUM_JACOBIAN = NUM_JACOBIAN + 1 ! keep track of the number of iterations - fuseStruct = ctx_sav ! restores consistency after finite differencing end SUBROUTINE jac_flux ! ----- simple implicit solve for differentiable model -------------------------- @@ -103,6 +102,7 @@ subroutine implicit_solve(fuseStruct, x0, x1, nx, ierr, message, isVerbose) 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 conserve_clamp_module, only: conserve_clamp ! fix overshoot and disaggregate fluxes to conserve mass 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 @@ -132,7 +132,6 @@ subroutine implicit_solve(fuseStruct, x0, x1, nx, ierr, message, isVerbose) integer(i4b), parameter :: maxit=100 ! maximum number of iterations logical(lgt) :: converged ! flag for convergence ! internal: backtracking line search w/ overshoot reject - type(parent) :: ctx ! save the fuse structure real(sp) :: xnorm ! norm used in maximum step real(sp) :: dxnorm ! norm used to evaluate step size real(sp) :: stpmax ! the maximum step @@ -153,6 +152,11 @@ subroutine implicit_solve(fuseStruct, x0, x1, nx, ierr, message, isVerbose) integer(i4b) :: ls_it ! index of line search iteration logical(lgt) :: ovshoot ! flag for overshoot logical(lgt) :: accepted ! flag for accepting newton step + real(sp) :: phi_best ! best function evaluation + real(sp) :: x_best(nx) ! best state vector + real(sp) :: g_best(nx) ! dx/dt = g(x_best) + logical(lgt) :: have_best ! check if found a state vector + logical(lgt) :: isClamped ! check if fallback is clamped ! algorithmic control parameters (most passed through MODULE model_numerix) REAL(SP), PARAMETER :: TOLMIN=1.0e-10_sp ! check for spurious minima REAL(SP), PARAMETER :: STPMX=100.0_sp ! maximum step in lnsrch @@ -169,6 +173,9 @@ subroutine implicit_solve(fuseStruct, x0, x1, nx, ierr, message, isVerbose) ! check dimension size if (nx /= nState) stop "implicit_solve: nx /= nState" + ! initialize check for best function evaluation + phi_best = huge(1._sp); have_best=.false. + ! initialize number of calls NUM_FUNCS = 0 ! number of function calls NUM_JACOBIAN = 0 ! number of times Jacobian is calculated @@ -176,9 +183,6 @@ subroutine implicit_solve(fuseStruct, x0, x1, nx, ierr, message, isVerbose) ! get the flag for printing isPrint = .false.; if (present(isVerbose)) isPrint = isVerbose - ! save the fuse structure - ctx = fuseStruct - ! 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) @@ -235,8 +239,6 @@ subroutine implicit_solve(fuseStruct, x0, x1, nx, ierr, message, isVerbose) exit ! exit iteration loop end if - if(isPrint) print*, 'x_try 0 = ', x_try - ! --- J(x) call jac_flux(fuseStruct, x_try, g_x, lower, upper, Jg) do i=1,nx @@ -244,8 +246,6 @@ subroutine implicit_solve(fuseStruct, x0, x1, nx, ierr, message, isVerbose) Jac(i,i) = Jac(i,i) + 1.0_sp end do - if(isPrint) print*, 'x_try 1 = ', x_try - ! --- function gradient: before Jac is modified in ludcmp gpsi = matmul(transpose(Jac), res) ! assumes func = 0.5_sp * sum(res*res) @@ -254,8 +254,6 @@ subroutine implicit_solve(fuseStruct, x0, x1, nx, ierr, message, isVerbose) call ludcmp(Jac, indx, d) ! J overwritten with LU call lubksb(Jac, indx, dx) ! dx becomes solution - if(isPrint) print*, 'x_try 2 = ', x_try - if(isPrint)then print*, 'dx = ', dx print*, 'Jg = ', Jg @@ -270,38 +268,22 @@ subroutine implicit_solve(fuseStruct, x0, x1, nx, ierr, message, isVerbose) dx = dxScale * dx end if + ! modify dx if Newton step not descending for psi + slope = dot_product(gpsi, dx) + if (slope >= 0._sp) dx = -gpsi ! fallback + ! implement active-set methods do i=1,nx if (x_try(i) <= lower(i)+epsb .and. dx(i) < 0._sp) dx(i)=0._sp if (x_try(i) >= upper(i)-epsb .and. dx(i) > 0._sp) dx(i)=0._sp end do - ! modify dx if Newton step not descending for psi - slope = dot_product(gpsi, dx) - if (slope >= 0._sp) dx = -gpsi ! fallback - - if(isPrint) print*, 'x_try 3 = ', x_try ! ---- backtracking line search -------------- - ! save the fuse structure to re-use in subsequent linesearch calls - ctx = fuseStruct ! <--- snapshot *now*, for this Newton step - ! line search control accepted = .false. ! flag to check if line search is accepted alamin = ERR_ITER_DX / maxval( abs(dx) / max(abs(x_try), 1.0_sp) ) - ! compute maximum lambda - lam_max = 1.0_sp - ! do i=1,nx - ! if (dx(i) > 0._sp) then - ! lam_i = (upper(i) - x_try(i) - epsb) / dx(i) - ! if (lam_i < lam_max) lam_max = max(0._sp, lam_i) - ! else if (dx(i) < 0._sp) then - ! lam_i = (lower(i) - x_try(i) + epsb) / dx(i) ! dx<0 so this is positive - ! if (lam_i < lam_max) lam_max = max(0._sp, lam_i) - ! end if - ! end do - ! check if(isPrint)then print*, 'alamin = ', alamin @@ -311,7 +293,7 @@ subroutine implicit_solve(fuseStruct, x0, x1, nx, ierr, message, isVerbose) if(isPrint) isDebug = .true. - lambda = lam_max + lambda = 1.0_sp do ls_it = 1, ls_max if(isPrint)then @@ -325,19 +307,19 @@ subroutine implicit_solve(fuseStruct, x0, x1, nx, ierr, message, isVerbose) if(isPrint)then print*, 'x_try = ', x_try print*, 'x_trial = ', x_trial + print*, 'lower = ', lower + print*, 'upper = ', upper print *, "delta = ", x_trial - x_try print *, "lambda*dx = ", lambda*dx endif - ! exit if violate bounds: line search direction is not valid + ! shrink lambda until find a value in the feasible space if(any(x_trial < lower) .or. any(x_trial > upper))then - accepted = .false.; exit - !lambda = lambda * shrink - !cycle + lambda = lambda * shrink + cycle endif ! compute function and function eval - fuseStruct = ctx g_trial = dx_dt(fuseStruct, x_trial) res_trial = x_trial - (x0 + dt*g_trial) phi_new = 0.5_sp * dot_product(res_trial, res_trial) @@ -349,7 +331,15 @@ subroutine implicit_solve(fuseStruct, x0, x1, nx, ierr, message, isVerbose) print*, 'g_trial = ', g_trial print*, 'res _trial= ', res_trial endif - + + ! save best function evaluation + if (phi_new < phi_best) then + phi_best = phi_new + x_best = x_trial + g_best = g_trial + have_best = .true. + endif + if (phi_new <= phi + phi_abs_tol) then accepted = .true.; exit endif @@ -375,10 +365,16 @@ subroutine implicit_solve(fuseStruct, x0, x1, nx, ierr, message, isVerbose) call fix_ovshoot(x_trial, lower, upper, dclamp) ! get new function evaluation x_try = x_trial - fuseStruct = ctx g_x = dx_dt(fuseStruct, x_try) res = x_try - (x0 + g_x*dt) phi = 0.5_sp * dot_product(res, res) + ! save best function evaluation + if (phi < phi_best) then + phi_best = phi + x_best = x_try + g_best = g_x + have_best = .true. + endif end if ! tiny-step convergence @@ -389,15 +385,25 @@ subroutine implicit_solve(fuseStruct, x0, x1, nx, ierr, message, isVerbose) end do ! loop through iterations + ! ----- handle the extremely rare case of non-convergence ----- + if( .not. converged)then + + ! use explicit Euler if did not find anything + if( .not. have_best) g_best = dx_dt(fuseStruct, x0) + + ! use dx/dt = g(x_best) + x_try = x0 + dt*g_best + + ! test bounds violations: if bounds exceeded, then clamp and disaggregate fluxes (conserve mass) + call XTRY_2_STR(x_try, fuseStruct%state1) + call conserve_clamp(fuseStruct, dt, isClamped) + print*, 'WARNING: '//trim(message)//"failed to converge: use best function evaluation. Clamp = ", isClamped + + endif ! if not converged + ! save final state x1 = x_try - ! check convergence - if( .not. converged)then - message=trim(message)//"failed to converge" - ierr=10; return - endif - end subroutine implicit_solve end module implicit_solve_module diff --git a/build/Makefile b/build/Makefile index b0c04f6..8b9f6f7 100644 --- a/build/Makefile +++ b/build/Makefile @@ -139,6 +139,7 @@ FUSE_PHYSICS= \ q_misscell_diff.f90 \ mstate_rhs_diff.f90 \ mod_derivs_diff.f90 \ + conserve_clamp.f90 \ fix_ovshoot.f90 \ implicit_solve.f90 PHYSICS = $(patsubst %, $(PHYSICS_DIR)/%, $(FUSE_PHYSICS)) From 26b355169259dd79c0f9c1a5d51fc1d20fe909cd Mon Sep 17 00:00:00 2001 From: Martyn Clark Date: Sun, 14 Dec 2025 08:47:51 -0700 Subject: [PATCH 9/9] initial implementation of analytical Jacobian --- build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 | 17 ++- build/FUSE_SRC/FUSE_ENGINE/str_2_xtry.f90 | 2 +- build/FUSE_SRC/FUSE_ENGINE/xtry_2_str.f90 | 6 +- build/FUSE_SRC/dshare/data_types.f90 | 1 + build/FUSE_SRC/dshare/globaldata.f90 | 3 + build/FUSE_SRC/physics/evap_lower_diff.f90 | 8 +- build/FUSE_SRC/physics/evap_upper_diff.f90 | 113 ++++++++++++----- build/FUSE_SRC/physics/get_parent.f90 | 13 +- build/FUSE_SRC/physics/implicit_solve.f90 | 141 +++++++-------------- build/FUSE_SRC/physics/mod_derivs_diff.f90 | 36 ++++-- build/FUSE_SRC/physics/mstate_rhs_diff.f90 | 65 ++++++---- build/FUSE_SRC/physics/q_baseflow_diff.f90 | 46 ++++++- build/FUSE_SRC/physics/q_misscell_diff.f90 | 17 ++- build/FUSE_SRC/physics/qinterflow_diff.f90 | 9 +- build/FUSE_SRC/physics/qpercolate_diff.f90 | 79 ++++++++++-- build/FUSE_SRC/physics/qsatexcess_diff.f90 | 62 ++++++++- build/FUSE_SRC/physics/smoothers.f90 | 107 ++++++++++++++++ 17 files changed, 525 insertions(+), 200 deletions(-) diff --git a/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 b/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 index e1d30e2..3f4c8a6 100644 --- a/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 +++ b/build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90 @@ -106,8 +106,12 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG ! --------------------------------------------------------------------------------------- ! allocate state vectors ALLOCATE(STATE0(NSTATE),STATE1(NSTATE),STAT=IERR) - IF (IERR.NE.0) STOP ' problem allocating space for state vectors in fuse_rmse ' + IF (IERR.NE.0) STOP ' problem allocating space for state vectors in fuse_rmse' + ! allocate flux derivative vector + allocate(fuseStruct%df_dS(nState), stat=ierr) + if(ierr/=0) STOP ' problem allocating space for the flux derivative vector' + ! increment parameter counter for model output IF (.NOT.PRESENT(MPARAM_FLAG)) THEN PCOUNT = PCOUNT + 1 @@ -131,6 +135,8 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG DO iSpat2=1,nSpat2 DO iSpat1=1,nSpat1 CALL INIT_STATE(fracState0) ! define FSTATE using fracState0 + CALL STR_2_XTRY(FSTATE,STATE0) ! set state at the start of the time step (STATE0) using FSTATE + CALL XTRY_2_STR(STATE0,FSTATE) ! update structure, including derived state variables gState_3d(iSpat1,iSpat2,1) = FSTATE ! put the state into first time step of 3D structure END DO END DO @@ -231,7 +237,7 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG FSTATE = gState_3d(iSpat1,iSpat2,itim_sub) MSTATE = FSTATE ! refresh model states CALL STR_2_XTRY(FSTATE,STATE0) ! set state at the start of the time step (STATE0) using FSTATE - + ! initialize model fluxes CALL INITFLUXES() ! set weighted sum of fluxes to zero @@ -273,6 +279,7 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG ! populate parent fuse structure call get_parent(fuseStruct) + ! solve differentiable ODEs call implicit_solve(fuseStruct, state0, state1, nState, ierr, cmessage) if(ierr/=0)then @@ -399,7 +406,8 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG CALL MEAN_STATS() RMSE = MSTATS%RAW_RMSE - write(*,'(i6,1x,a6,1x,f12.6,1x)') nFUSE_eval, "NSE = ", MSTATS%NASH_SUTT + write(*,'(i6,1x,a6,1x,f12.6,1x,a20,1x,f12.6)') nFUSE_eval, "NSE = ", MSTATS%NASH_SUTT, "; TIME ELAPSED = ", t2-t1 + !if(nFUSE_eval > 10) stop "checking results" ENDIF @@ -408,7 +416,8 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG ! deallocate vectors DEALLOCATE(W_FLUX_3d); IF (IERR.NE.0) STOP ' problem deallocating W_FLUX_3d in fuse_rmse ' - DEALLOCATE(STATE0,STATE1,STAT=IERR); IF (IERR.NE.0) STOP ' problem deallocating state vectors in fuse_rmse ' + DEALLOCATE(STATE0,STATE1,STAT=IERR); IF (IERR.NE.0) STOP ' problem deallocating state vectors in fuse_rmse' + deallocate(fuseStruct%df_dS, stat=ierr); if(ierr/=0) STOP ' problem deallocating space for the flux derivative vector' END SUBROUTINE FUSE_RMSE END MODULE FUSE_RMSE_MODULE diff --git a/build/FUSE_SRC/FUSE_ENGINE/str_2_xtry.f90 b/build/FUSE_SRC/FUSE_ENGINE/str_2_xtry.f90 index 71875e9..0da2581 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/str_2_xtry.f90 +++ b/build/FUSE_SRC/FUSE_ENGINE/str_2_xtry.f90 @@ -12,9 +12,9 @@ SUBROUTINE STR_2_XTRY(TMPSTR,X_TRY) ! Copy model states into vector X_TRY ! --------------------------------------------------------------------------------------- USE nrtype ! Numerical Recipes data types +USE data_types, ONLY: STATEV ! model state structure USE model_defn, ONLY: CSTATE,NSTATE ! model definitions USE model_defnames -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/xtry_2_str.f90 b/build/FUSE_SRC/FUSE_ENGINE/xtry_2_str.f90 index b68b583..33ed995 100644 --- a/build/FUSE_SRC/FUSE_ENGINE/xtry_2_str.f90 +++ b/build/FUSE_SRC/FUSE_ENGINE/xtry_2_str.f90 @@ -16,10 +16,10 @@ SUBROUTINE XTRY_2_STR(X_TRY,TMPSTR) ! Temporary model states updated in MODULE multistate ! --------------------------------------------------------------------------------------- USE nrtype ! Numerical Recipes data types +USE data_types, ONLY: STATEV ! model state definitions USE model_defn, ONLY: CSTATE,NSTATE,SMODL ! model definitions USE model_defnames -USE multistate, ONLY: STATEV ! model states -USE multiparam, ONLY: DPARAM ! model parameters +USE multiparam, ONLY: DPARAM ! derived model parameters IMPLICIT NONE ! input REAL(SP), DIMENSION(:), INTENT(IN) :: X_TRY ! vector of model states @@ -73,7 +73,7 @@ SUBROUTINE XTRY_2_STR(X_TRY,TMPSTR) TMPSTR%TENS_2 = MIN(TMPSTR%WATR_2, DPARAM%MAXTENS_2) ! tension storage TMPSTR%FREE_2 = MAX(0._sp, TMPSTR%WATR_2 - DPARAM%MAXTENS_2) ! free storage TMPSTR%FREE_2A = missingValue ! primary reservoir (undefined) - TMPSTR%FREE_2A = missingValue ! secondary reservoir (undefined) + TMPSTR%FREE_2B = missingValue ! secondary reservoir (undefined) CASE DEFAULT ! (error check) print *, "MDEFN(IMOD)%ARCH2 must be 'tens2pll_2', 'unlimfrc_2', 'unlimpow_2'" print *, " 'topmdexp_2', or 'fixedsiz_2'" diff --git a/build/FUSE_SRC/dshare/data_types.f90 b/build/FUSE_SRC/dshare/data_types.f90 index b27a0f7..afbc701 100644 --- a/build/FUSE_SRC/dshare/data_types.f90 +++ b/build/FUSE_SRC/dshare/data_types.f90 @@ -310,6 +310,7 @@ module data_types type(statev) :: state1 ! state variables (end of step) type(statev) :: dx_dt ! time derivative in state variables type(fluxes) :: flux ! fluxes + type(fluxes), allocatable :: df_dS(:) ! derivative in fluxes w.r.t. states type(runoff) :: route ! hillslope routing type(par_id) :: param_name ! parameter names type(parinfo) :: param_meta ! metadata on model parameters diff --git a/build/FUSE_SRC/dshare/globaldata.f90 b/build/FUSE_SRC/dshare/globaldata.f90 index 9d551d3..bd95d54 100644 --- a/build/FUSE_SRC/dshare/globaldata.f90 +++ b/build/FUSE_SRC/dshare/globaldata.f90 @@ -15,6 +15,9 @@ MODULE globaldata ! initial store fraction (initialization) real(sp), parameter :: fracState0=0.25_sp + ! original code + logical(lgt), save :: isOriginal=.true. + ! print flag logical(lgt), save :: isPrint=.true. logical(lgt), save :: isDebug=.false. diff --git a/build/FUSE_SRC/physics/evap_lower_diff.f90 b/build/FUSE_SRC/physics/evap_lower_diff.f90 index a07a76a..96a1f2d 100644 --- a/build/FUSE_SRC/physics/evap_lower_diff.f90 +++ b/build/FUSE_SRC/physics/evap_lower_diff.f90 @@ -7,7 +7,7 @@ module EVAP_LOWER_DIFF_MODULE contains - SUBROUTINE EVAP_LOWER_DIFF(fuseStruct) + SUBROUTINE EVAP_LOWER_DIFF(fuseStruct, want_dflux) ! ------------------------------------------------------------------------------------------------- ! Creator: ! -------- @@ -25,6 +25,9 @@ SUBROUTINE EVAP_LOWER_DIFF(fuseStruct) IMPLICIT NONE ! input-output type(parent), intent(inout) :: fuseStruct ! parent fuse data structure + logical(lgt), intent(in), optional :: want_dflux ! if we want flux derivatives + ! internal + logical(lgt) :: comp_dflux ! flag to compute flux derivatives ! ------------------------------------------------------------------------------------------------- ! associate variables with elements of data structure associate(& @@ -36,6 +39,9 @@ SUBROUTINE EVAP_LOWER_DIFF(fuseStruct) ) ! (associate) ! ------------------------------------------------------------------------------------------------- + ! check the need to compute flux derivatives + comp_dflux = .false.; if(present(want_dflux)) comp_dflux = want_dflux + ! --------------------------------------------------------------------------------------- SELECT CASE(SMODL%iARCH2) ! lower layer architecture CASE(iopt_tens2pll_2,iopt_fixedsiz_2) diff --git a/build/FUSE_SRC/physics/evap_upper_diff.f90 b/build/FUSE_SRC/physics/evap_upper_diff.f90 index fa0fd2d..030fa60 100644 --- a/build/FUSE_SRC/physics/evap_upper_diff.f90 +++ b/build/FUSE_SRC/physics/evap_upper_diff.f90 @@ -7,7 +7,7 @@ module EVAP_UPPER_DIFF_module contains - SUBROUTINE EVAP_UPPER_DIFF(fuseStruct) + SUBROUTINE EVAP_UPPER_DIFF(fuseStruct, want_dflux) ! ------------------------------------------------------------------------------------------------- ! Creator: ! -------- @@ -21,21 +21,38 @@ SUBROUTINE EVAP_UPPER_DIFF(fuseStruct) USE nrtype ! variable types, etc. USE data_types, only: parent ! fuse parent data type USE model_defn ! model definition structure - USE model_defnames + USE model_defnames ! model definition names + use smoothers, only : sfrac, dsfrac ! smoothed fraction, derivative IMPLICIT NONE ! input-output type(parent), intent(inout) :: fuseStruct ! parent fuse data structure + logical(lgt), intent(in), optional :: want_dflux ! if we want flux derivatives + ! local variables + logical(lgt) :: comp_dflux ! flag to compute flux derivatives + integer(i4b) :: iState ! state index + real(sp) :: phi ! smoothed fraction of total tension storage (0,1] + real(sp) :: phi_1a ! smoothed fraction of primary tension storage (0,1] + real(sp) :: phi_1b ! smoothed fraction of secondary tension storage (0,1] + real(sp) :: maxRate ! maximum forcing + real(sp) :: maxRate_1a ! maximum forcing for the primary tension tank + real(sp) :: maxRate_1b ! maximum forcing for the secondary tension tank + real(sp) :: dphi_dx ! derivative in fraction w.r.t. storage + real(sp) :: devap_dx ! derivative in evaporation w.r.t. storage ! ------------------------------------------------------------------------------------------------- ! associate variables with elements of data structure associate(& MFORCE => fuseStruct%force , & ! model forcing data M_FLUX => fuseStruct%flux , & ! fluxes + dfx_dS => fuseStruct%df_dS , & ! deriv in fluxes w.r.t. states TSTATE => fuseStruct%state1 , & ! trial state variables (end of step) MPARAM => fuseStruct%param_adjust , & ! adjustable model parameters DPARAM => fuseStruct%param_derive & ! derived model parameters ) ! (associate) ! ------------------------------------------------------------------------------------------------- + ! check the need to compute flux derivatives + comp_dflux = .false.; if(present(want_dflux)) comp_dflux = want_dflux + ! --------------------------------------------------------------------------------------- SELECT CASE(SMODL%iARCH1) ! upper layer architecture @@ -43,47 +60,77 @@ SUBROUTINE EVAP_UPPER_DIFF(fuseStruct) CASE(iopt_tension2_1) ! tension storage sub-divided into recharge and excess ! -------------------------------------------------------------------------------------- - ! use different evaporation schemes for the upper layer - ! ----------------------------------------------------- + ! calculate the smoothed fraction of tension storage (NOTE: use WATR_1) + phi_1a = sfrac(TSTATE%TENS_1A, DPARAM%MAXTENS_1A) + phi_1b = sfrac(TSTATE%TENS_1B, DPARAM%MAXTENS_1B) + + ! calculate the maximum evap rate for the storage 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 + maxrate_1a = MFORCE%PET + maxrate_1b = MFORCE%PET - MFORCE%PET*phi_1a 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 + maxrate_1a = MFORCE%PET * MPARAM%RTFRAC1 + maxrate_1b = MFORCE%PET * DPARAM%RTFRAC2 + CASE DEFAULT; stop "evap_upper: SMODL%iESOIL must be either iopt_sequential or iopt_rootweight" END SELECT + + ! ----- compute flux ---------------------------------------------------------------- + M_FLUX%EVAP_1A = maxrate_1a*phi_1a + M_FLUX%EVAP_1B = maxrate_1b*phi_1b + M_FLUX%EVAP_1 = M_FLUX%EVAP_1A + M_FLUX%EVAP_1B + + ! ----- compute derivatives --------------------------------------------------------------------- + if(comp_dflux) stop "evap_upper: derivatives for iopt_tension2_1 not implemented yet" + ! -------------------------------------------------------------------------------------- CASE(iopt_tension1_1,iopt_onestate_1) ! single tension store or single state ! -------------------------------------------------------------------------------------- - ! use different evaporation schemes for the upper layer - ! ----------------------------------------------------- + ! zero fluxes not used + M_FLUX%EVAP_1A = 0._sp + M_FLUX%EVAP_1B = 0._sp + + select case(SMODL%iARCH1) + case(iopt_tension1_1); phi = sfrac(TSTATE%TENS_1, DPARAM%MAXTENS_1) + case(iopt_onestate_1); phi = sfrac(TSTATE%WATR_1, DPARAM%MAXTENS_1) ! NOTE: use WATR_1 + end select ! no need for default because checked above + + ! calculate the maximum evap rate 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" + CASE(iopt_sequential); maxRate = MFORCE%PET + CASE(iopt_rootweight); maxRate = MFORCE%PET*MPARAM%RTFRAC1 + CASE DEFAULT; stop "evap_upper: 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) + + ! ----- compute flux ---------------------------------------------------------------- + M_FLUX%EVAP_1 = maxRate*phi + + ! ----- compute derivatives --------------------------------------------------------- + if(comp_dflux)then + + ! calculate the derivative in the smoothed fraction of tension storage + select case(SMODL%iARCH1) + case(iopt_tension1_1); dphi_dx = dsfrac(TSTATE%TENS_1, DPARAM%MAXTENS_1) + case(iopt_onestate_1); dphi_dx = dsfrac(TSTATE%WATR_1, DPARAM%MAXTENS_1) ! NOTE: use WATR_1 + end select ! no need for default because checked above + + ! calculate the derivative in the maximum rate + devap_dx = maxRate*dphi_dx + + ! populate derivative vector + do iState=1,nState + select case(cState(iState)%iSNAME) + case (iopt_TENS_1); dfx_dS(iState)%EVAP_1 = devap_dx ! exists if one tension tank + case (iopt_WATR_1); dfx_dS(iState)%EVAP_1 = devap_dx ! exists if one state in the upper layer + end select ! no default needed + end do ! looping through states + + endif ! if computing derivatives + + CASE DEFAULT; stop "evap_upper: SMODL%iARCH1 must be iopt_tension2_1, iopt_tension1_1, or iopt_onestate_1" + END SELECT ! (upper-layer architecture) + end associate ! end association with variables in the data structures END SUBROUTINE EVAP_UPPER_DIFF diff --git a/build/FUSE_SRC/physics/get_parent.f90 b/build/FUSE_SRC/physics/get_parent.f90 index 1a79e0d..c142bb2 100644 --- a/build/FUSE_SRC/physics/get_parent.f90 +++ b/build/FUSE_SRC/physics/get_parent.f90 @@ -1,5 +1,7 @@ module get_parent_module + use nrtype use data_types, only: parent + USE model_defn, ONLY:NSTATE implicit none contains @@ -10,16 +12,23 @@ subroutine get_parent(fuseStruct) use multi_flux, only: m_flux use multiparam, only: parMeta,mParam,dParam implicit none - type(parent), intent(out) :: fuseStruct + type(parent), intent(inout) :: fuseStruct + integer(i4b) :: iState + ! populate parent fuse structures fuseStruct%force = mForce fuseStruct%state0 = mState fuseStruct%state1 = mState - fuseStruct%flux = m_flux + fuseStruct%flux = m_flux ! initialized at zero fuseStruct%param_meta = parMeta fuseStruct%param_adjust = mParam fuseStruct%param_derive = dParam + ! initialize derivatives + do iState=1,nState + fuseStruct%df_dS(iState) = m_flux ! initialized at zero + end do + end subroutine get_parent diff --git a/build/FUSE_SRC/physics/implicit_solve.f90 b/build/FUSE_SRC/physics/implicit_solve.f90 index db92b58..f008e29 100644 --- a/build/FUSE_SRC/physics/implicit_solve.f90 +++ b/build/FUSE_SRC/physics/implicit_solve.f90 @@ -24,28 +24,35 @@ module implicit_solve_module contains ! ----- calculate dx/dt=g(x) ----------------------------------------------------------- - function dx_dt(fuseStruct, x_try) result(g_x) - use MOD_DERIVS_DIFF_module, only: MOD_DERIVS_DIFF ! compute dx/dt + subroutine dx_dt(fuseStruct, x_try, g_x, J_g) + use MOD_DERIVS_DIFF_module, only: MOD_DERIVS_DIFF ! compute dx/dt implicit none ! input - type(parent) , intent(inout) :: fuseStruct ! parent fuse data structure - real(sp) , intent(in) :: x_try(:) ! trial state vector + type(parent) , intent(inout) :: fuseStruct ! parent fuse data structure + real(sp) , intent(in) :: x_try(:) ! trial state vector ! output - real(sp) :: g_x(size(x_try)) ! dx/dt=g(x) + real(sp) , intent(out) :: g_x(:) ! dx/dt=g(x) + real(sp) , intent(out) , optional :: J_g(:,:) ! flux Jacobian matrix + ! internal + logical(lgt) :: comp_dflux ! flag to compute flux derivatives + ! -------------------------------------------------------------------------------------- + + comp_dflux = present(J_g) ! put data in structure call XTRY_2_STR(x_try, fuseStruct%state1) ! run the fuse physics - call mod_derivs_diff(fuseStruct) + if (present(J_g)) then + call mod_derivs_diff(fuseStruct, g_x, J_g) + else + call mod_derivs_diff(fuseStruct, g_x) + end if - ! extract dx_dt from fuse structure - call STR_2_XTRY(fuseStruct%dx_dt, g_x) - ! track the total number of function calls NUM_FUNCS = NUM_FUNCS + 1 - end function dx_dt + end subroutine dx_dt ! ----- calculate the Jacobian of g(x) ------------------------------------------------- SUBROUTINE jac_flux(fuseStruct, x_try, g_x, lower, upper, Jac) @@ -81,7 +88,7 @@ SUBROUTINE jac_flux(fuseStruct, x_try, g_x, lower, upper, Jac) ! compute function from the perturbed vector x(j) = xsav(j) + h_try - g_ph = dx_dt(fuseStruct_local, x) + call dx_dt(fuseStruct_local, x, g_ph) h_act = x(j) - xsav(j) ! compute column in the Jacobian @@ -121,6 +128,7 @@ subroutine implicit_solve(fuseStruct, x0, x1, nx, ierr, message, isVerbose) 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) :: Ja(nx,nx) ! Jacobian matrix (flux) real(sp) :: Jg(nx,nx) ! Jacobian matrix (flux) real(sp) :: Jac(nx,nx) ! Jacobian matrix (full) real(sp) :: dx(nx) ! state update @@ -207,42 +215,27 @@ subroutine implicit_solve(fuseStruct, x0, x1, nx, ierr, message, isVerbose) accepted = .false. converged = .false. - if(isPrint) isDebug = .true. - - ! --- F(x) and objective phi - g_x = dx_dt(fuseStruct, x_try) + ! --- F(x), J(x), and objective phi + call dx_dt(fuseStruct, x_try, g_x, Jg) ! compute analytical Jacobian res = x_try - (x0 + g_x*dt) phi = 0.5_sp * dot_product(res, res) - if(isPrint) isDebug = .false. - ! iterate do it = 1, maxit ! save x x_old = x_try - if(isPrint) print*, '***** start of iteration *****' - - ! check - if(isPrint)then - print*, 'x_try = ', x_try - print*, 'g_x = ', g_x - print*, 'res = ', res - print*, 'phi = ', phi - print*, 'dclamp = ', dclamp - if(it > 10) stop 1 - endif - + ! check convergence if (phi < ERR_ITER_FUNC) then converged = .true. exit ! exit iteration loop end if - ! --- J(x) - call jac_flux(fuseStruct, x_try, g_x, lower, upper, Jg) + ! --- compute residual Jacobian J(x) from flux Jacobian Jg(x) ---- + !call jac_flux(fuseStruct, x_try, g_x, lower, upper, Jg) do i=1,nx - Jac(:,i) = -dt*Jg(:,i) !* dclamp(i) ! multiply dt and clamp derivative + Jac(:,i) = -dt*Jg(:,i) Jac(i,i) = Jac(i,i) + 1.0_sp end do @@ -254,11 +247,6 @@ subroutine implicit_solve(fuseStruct, x0, x1, nx, ierr, message, isVerbose) call ludcmp(Jac, indx, d) ! J overwritten with LU call lubksb(Jac, indx, dx) ! dx becomes solution - if(isPrint)then - print*, 'dx = ', dx - print*, 'Jg = ', Jg - endif - ! --- Modify dx ! modify dx if norm > stpmax @@ -284,54 +272,23 @@ subroutine implicit_solve(fuseStruct, x0, x1, nx, ierr, message, isVerbose) accepted = .false. ! flag to check if line search is accepted alamin = ERR_ITER_DX / maxval( abs(dx) / max(abs(x_try), 1.0_sp) ) - ! check - if(isPrint)then - print*, 'alamin = ', alamin - print*, 'slope = ', slope - print*, 'gpsi = ', gpsi - endif - - if(isPrint) isDebug = .true. - lambda = 1.0_sp do ls_it = 1, ls_max - if(isPrint)then - print*, '***** new linesearch *****', ls_it - print*, 'dx = ', dx - endif - ! update x x_trial = x_try + lambda*dx - if(isPrint)then - print*, 'x_try = ', x_try - print*, 'x_trial = ', x_trial - print*, 'lower = ', lower - print*, 'upper = ', upper - print *, "delta = ", x_trial - x_try - print *, "lambda*dx = ", lambda*dx - endif - ! shrink lambda until find a value in the feasible space if(any(x_trial < lower) .or. any(x_trial > upper))then lambda = lambda * shrink cycle endif - ! compute function and function eval - g_trial = dx_dt(fuseStruct, x_trial) + ! compute function and function eval -- no need for the Jacobian here + call dx_dt(fuseStruct, x_trial, g_trial) res_trial = x_trial - (x0 + dt*g_trial) phi_new = 0.5_sp * dot_product(res_trial, res_trial) - if(isPrint)then - print*, 'ls_it, lambda, phi, phi_new', ls_it, lambda, phi, phi_new - print*, 'phi, phi_new, slope=', phi, phi_new, slope - print*, 'x_trial = ', x_trial - print*, 'g_trial = ', g_trial - print*, 'res _trial= ', res_trial - endif - ! save best function evaluation if (phi_new < phi_best) then phi_best = phi_new @@ -350,32 +307,26 @@ subroutine implicit_solve(fuseStruct, x0, x1, nx, ierr, message, isVerbose) end do ! line search - if(isPrint) isDebug = .false. - - if (accepted) then - x_try = x_trial - g_x = g_trial - res = res_trial - phi = phi_new - else - ! ----- fallback: try a small step along the direction of steepest descent ----- - !dx = -gpsi ! use steepest descent + ! ----- fallback: try a small step ----- + if(.not. accepted)then x_trial = x_try + dampen*dx if(any(x_trial < lower) .or. any(x_trial > upper)) & call fix_ovshoot(x_trial, lower, upper, dclamp) - ! get new function evaluation - x_try = x_trial - g_x = dx_dt(fuseStruct, x_try) - res = x_try - (x0 + g_x*dt) - phi = 0.5_sp * dot_product(res, res) - ! save best function evaluation - if (phi < phi_best) then - phi_best = phi - x_best = x_try - g_best = g_x - have_best = .true. - endif - end if + end if ! (if accepted) + + ! recompute dx_dt because we need the Jacobian + x_try = x_trial + call dx_dt(fuseStruct, x_try, g_x, Jg) ! compute analytical Jacobian + res = x_try - (x0 + g_x*dt) + phi = 0.5_sp * dot_product(res, res) + + ! save best function evaluation + if (phi < phi_best) then + phi_best = phi + x_best = x_try + g_best = g_x + have_best = .true. + endif ! tiny-step convergence if (maxval( abs(x_try - x_old) / max(abs(x_try), 1._sp) ) < ERR_ITER_DX) then @@ -389,7 +340,7 @@ subroutine implicit_solve(fuseStruct, x0, x1, nx, ierr, message, isVerbose) if( .not. converged)then ! use explicit Euler if did not find anything - if( .not. have_best) g_best = dx_dt(fuseStruct, x0) + if( .not. have_best) call dx_dt(fuseStruct, x0, g_best) ! use dx/dt = g(x_best) x_try = x0 + dt*g_best @@ -398,7 +349,7 @@ subroutine implicit_solve(fuseStruct, x0, x1, nx, ierr, message, isVerbose) call XTRY_2_STR(x_try, fuseStruct%state1) call conserve_clamp(fuseStruct, dt, isClamped) print*, 'WARNING: '//trim(message)//"failed to converge: use best function evaluation. Clamp = ", isClamped - + endif ! if not converged ! save final state diff --git a/build/FUSE_SRC/physics/mod_derivs_diff.f90 b/build/FUSE_SRC/physics/mod_derivs_diff.f90 index 5dc1752..5bdd3ce 100644 --- a/build/FUSE_SRC/physics/mod_derivs_diff.f90 +++ b/build/FUSE_SRC/physics/mod_derivs_diff.f90 @@ -18,7 +18,7 @@ module MOD_DERIVS_DIFF_module contains - SUBROUTINE MOD_DERIVS_DIFF(fuseStruct) + SUBROUTINE MOD_DERIVS_DIFF(fuseStruct, g_x, J_g) ! --------------------------------------------------------------------------------------- ! Creator: ! -------- @@ -29,21 +29,35 @@ SUBROUTINE MOD_DERIVS_DIFF(fuseStruct) ! Purpose: ! -------- ! compute the time derivative (dx/dt) of all model states (x) - ! --------------------------------------------------------------------------------------- + ! -------------------------------------------------------------------------------------- implicit none - type(parent), intent(inout) :: fuseStruct ! parent fuse data structure + ! input + type(parent) , intent(inout) :: fuseStruct ! parent fuse data structure + ! output + real(sp) , intent(out) :: g_x(:) ! dx/dt=g(x) + real(sp) , intent(out) , optional :: J_g(:,:) ! flux Jacobian matrix + ! internal + logical(lgt) :: comp_dflux ! flag to compute flux derivatives + ! -------------------------------------------------------------------------------------- + + ! check if Jacobian is desired + comp_dflux = present(J_g) ! 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) + call qsatexcess_diff(fuseStruct, comp_dflux) ! compute the saturated area and surface runoff + call evap_upper_diff(fuseStruct, comp_dflux) ! compute evaporation from the upper layer + call evap_lower_diff(fuseStruct, comp_dflux) ! compute evaporation from the lower layer + call qinterflow_diff(fuseStruct, comp_dflux) ! compute interflow from free water in the upper layer + call qpercolate_diff(fuseStruct, comp_dflux) ! compute percolation from the upper to lower soil layers + call q_baseflow_diff(fuseStruct, comp_dflux) ! compute baseflow from the lower soil layer + call q_misscell_diff(fuseStruct, comp_dflux) ! 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) + if(comp_dflux)then + call mstate_rhs_diff(fuseStruct, g_x, J_g) + else + call mstate_rhs_diff(fuseStruct, g_x) + endif END SUBROUTINE MOD_DERIVS_DIFF diff --git a/build/FUSE_SRC/physics/mstate_rhs_diff.f90 b/build/FUSE_SRC/physics/mstate_rhs_diff.f90 index b34eb90..c73502a 100644 --- a/build/FUSE_SRC/physics/mstate_rhs_diff.f90 +++ b/build/FUSE_SRC/physics/mstate_rhs_diff.f90 @@ -9,7 +9,7 @@ module MSTATE_RHS_DIFF_module contains - SUBROUTINE MSTATE_RHS_DIFF(fuseStruct) + SUBROUTINE MSTATE_RHS_DIFF(fuseStruct, g_x, J_g) ! --------------------------------------------------------------------------------------- ! Creator: ! -------- @@ -20,12 +20,18 @@ SUBROUTINE MSTATE_RHS_DIFF(fuseStruct) ! -------- ! 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 + USE nrtype ! variable types, etc. + USE data_types, only: parent ! fuse parent data type + USE model_defn ! model definition structure + USE model_defnames ! model names + use str_2_xtry_module ! puts FUSE state structure into state vector ! input-output type(parent), intent(inout) :: fuseStruct ! parent fuse data structure + ! output + real(sp) , intent(out) :: g_x(:) ! dx/dt=g(x) + real(sp) , intent(out) , optional :: J_g(:,:) ! flux Jacobian matrix + ! internal + logical(lgt) :: comp_dflux ! flag to compute flux derivatives ! ------------------------------------------------------------------------------------------------- ! associate variables with elements of data structure associate(& @@ -33,12 +39,16 @@ SUBROUTINE MSTATE_RHS_DIFF(fuseStruct) MPARAM => fuseStruct%param_adjust , & ! adjustable model parameters DX_DT => fuseStruct%dx_dt & ! time derivative in states ) ! (associate) + ! ------------------------------------------------------------------------------------------------- - if(isDebug) print*, 'M_FLUX%QPERC_12 = ', M_FLUX%QPERC_12 + ! check if Jacobian is desired + comp_dflux = present(J_g) ! --------------------------------------------------------------------------------------- - ! (1) COMPUTE TIME DERIVATIVES FOR STATES IN THE UPPER LAYER + ! (1) UPPER LAYER ! --------------------------------------------------------------------------------------- + + ! compute time derivatives 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 @@ -53,11 +63,19 @@ SUBROUTINE MSTATE_RHS_DIFF(fuseStruct) CASE DEFAULT print *, "SMODL%iARCH1 must be iopt_tension2_1, iopt_tension1_1, or iopt_onestate_1" STOP - END SELECT ! (upper layer architechure) + END SELECT ! (upper layer architecture) + + ! compute Jacobian + if(comp_dflux)then + if(SMODL%iARCH1 /= iopt_onestate_1) stop "mstate_rhs: only iopt_onestate_1 currently implemented" + J_g(1,:) = -M_FLUX%EFF_PPT*fuseStruct%df_dS%SATAREA - fuseStruct%df_dS%EVAP_1 - fuseStruct%df_dS%QPERC_12 + endif ! --------------------------------------------------------------------------------------- - ! (2) COMPUTE TIME DERIVATIVES FOR STATES IN THE LOWER LAYER + ! (2) LOWER LAYER ! --------------------------------------------------------------------------------------- + + ! compute time derivatives 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 @@ -73,22 +91,23 @@ SUBROUTINE MSTATE_RHS_DIFF(fuseStruct) print *, " iopt_topmdexp_2, or iopt_fixedsiz_2" STOP END SELECT + + ! compute Jacobian + ! NOTE: assume M_FLUX%EVAP_2=0 and M_FLUX%OFLOW_2=0 + if(comp_dflux)then + if(SMODL%iARCH2 == iopt_tens2pll_2) stop "mstate_rhs: iopt_tens2pll_2 not currently implemented" + J_g(2,:) = fuseStruct%df_dS%QPERC_12 - fuseStruct%df_dS%QBASE_2 + endif + ! --------------------------------------------------------------------------------------- - if(isDebug) print*, '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 = ', & - 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 - - if(isDebug) print*, 'DX_DT%WATR_2, M_FLUX%QPERC_12, M_FLUX%EVAP_2, M_FLUX%QBASE_2, M_FLUX%OFLOW_2 = ', & - DX_DT%WATR_2, M_FLUX%QPERC_12, M_FLUX%EVAP_2, M_FLUX%QBASE_2, M_FLUX%OFLOW_2 - - ! if(isDebug) print*, 'DX_DT%TENS_1, M_FLUX%EFF_PPT, M_FLUX%QSURF, M_FLUX%EVAP_1, M_FLUX%TENS2FREE_1 = ', & - ! DX_DT%TENS_1, M_FLUX%EFF_PPT, M_FLUX%QSURF, M_FLUX%EVAP_1, M_FLUX%TENS2FREE_1 - ! - ! if(isDebug) print*, 'DX_DT%TENS_2, M_FLUX%QPERC_12*(1._SP-MPARAM%PERCFRAC), M_FLUX%EVAP_2, M_FLUX%TENS2FREE_2 = ', & - ! DX_DT%TENS_2, M_FLUX%QPERC_12*(1._SP-MPARAM%PERCFRAC), M_FLUX%EVAP_2, M_FLUX%TENS2FREE_2 - ! - ! if(isDebug) print*, '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 = ', & - ! 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 + ! --------------------------------------------------------------------------------------- + ! (3) FINALIZE + ! --------------------------------------------------------------------------------------- + + ! extract dx_dt from fuse structure + call STR_2_XTRY(fuseStruct%dx_dt, g_x) + ! --------------------------------------------------------------------------------------- end associate ! end association with variables in the data structures END SUBROUTINE MSTATE_RHS_DIFF diff --git a/build/FUSE_SRC/physics/q_baseflow_diff.f90 b/build/FUSE_SRC/physics/q_baseflow_diff.f90 index d63adb6..805f42b 100644 --- a/build/FUSE_SRC/physics/q_baseflow_diff.f90 +++ b/build/FUSE_SRC/physics/q_baseflow_diff.f90 @@ -1,7 +1,5 @@ module Q_BASEFLOW_DIFF_module - USE globaldata, only: isDebug - implicit none private @@ -10,7 +8,7 @@ module Q_BASEFLOW_DIFF_module contains - SUBROUTINE Q_BASEFLOW_DIFF(fuseStruct) + SUBROUTINE Q_BASEFLOW_DIFF(fuseStruct, want_dflux) ! --------------------------------------------------------------------------------------- ! Creator: ! -------- @@ -28,43 +26,83 @@ SUBROUTINE Q_BASEFLOW_DIFF(fuseStruct) IMPLICIT NONE ! input-output type(parent), intent(inout) :: fuseStruct ! parent fuse data structure + logical(lgt), intent(in), optional :: want_dflux ! if we want flux derivatives + ! derivatives + logical(lgt) :: comp_dflux ! flag to compute flux derivatives + integer(i4b) :: iState ! state index + real(sp) :: phi ! scaled water storage, phi=w/ws + real(sp) :: dqb_dw ! derivative in baseflow flux w.r.t. water store ! ------------------------------------------------------------------------------------------------- ! associate variables with elements of data structure associate(& M_FLUX => fuseStruct%flux , & ! fluxes + dfx_dS => fuseStruct%df_dS , & ! deriv in fluxes w.r.t. states TSTATE => fuseStruct%state1 , & ! trial state variables (end of step) MPARAM => fuseStruct%param_adjust , & ! adjustable model parameters DPARAM => fuseStruct%param_derive & ! derived model parameters ) ! (associate) + ! check the need to compute flux derivatives + comp_dflux = .false.; if(present(want_dflux)) comp_dflux = want_dflux + ! --------------------------------------------------------------------------------------- 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 + if(comp_dflux) stop "q_baseflow: derivative not implemented yet for iopt_tens2pll_2" + ! -------------------------------------------------------------------------------------- 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) + if(comp_dflux) stop "q_baseflow: derivative not implemented yet for iopt_unlimfrc_2" + ! -------------------------------------------------------------------------------------- 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 + + associate(qbsat=>DPARAM%QBSAT, w=>TSTATE%WATR_2, ws=>MPARAM%MAXWATR_2, p=>MPARAM%QB_POWR) + + ! ----- compute flux ------------------------------------------------------------------ + phi = w/ws + M_FLUX%QBASE_2 = qbsat*phi**p + + ! ----- compute derivative ------------------------------------------------------------ + if(comp_dflux) dqb_dw = (qbsat*p/ws)*phi**(p - 1._sp) + + end associate + ! -------------------------------------------------------------------------------------- CASE(iopt_topmdexp_2) ! topmodel exponential reservoir (-HUGE to HUGE) M_FLUX%QBASE_2 = DPARAM%QBSAT * EXP( -(1. - TSTATE%WATR_2/MPARAM%MAXWATR_2) ) + if(comp_dflux) stop "q_baseflow: derivative not implemented yet for iopt_topmdexp_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 + if(comp_dflux) stop "q_baseflow: derivative not implemented yet for iopt_fixedsiz_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 ! --------------------------------------------------------------------------------------- + ! populate derivative vector + if(comp_dflux)then + do iState=1,nState + select case(cState(iState)%iSNAME) + case (iopt_WATR_2); dfx_dS(iState)%QBASE_2 = dqb_dw ! exists if one state in the upper layer + end select ! no default needed + end do ! looping through states + endif + end associate ! end association with variables in the data structures END SUBROUTINE Q_BASEFLOW_DIFF diff --git a/build/FUSE_SRC/physics/q_misscell_diff.f90 b/build/FUSE_SRC/physics/q_misscell_diff.f90 index c347455..8282832 100644 --- a/build/FUSE_SRC/physics/q_misscell_diff.f90 +++ b/build/FUSE_SRC/physics/q_misscell_diff.f90 @@ -1,5 +1,5 @@ module Q_MISSCELL_DIFF_module - + implicit none private @@ -7,7 +7,7 @@ module Q_MISSCELL_DIFF_module contains - SUBROUTINE Q_MISSCELL_DIFF(fuseStruct) + SUBROUTINE Q_MISSCELL_DIFF(fuseStruct, want_dflux) ! --------------------------------------------------------------------------------------- ! Creator: ! -------- @@ -31,7 +31,9 @@ SUBROUTINE Q_MISSCELL_DIFF(fuseStruct) IMPLICIT NONE ! input-output type(parent), intent(inout) :: fuseStruct ! parent fuse data structure + logical(lgt), intent(in), optional :: want_dflux ! if we want flux derivatives ! internal + logical(lgt) :: comp_dflux ! flag to compute flux derivatives REAL(SP), PARAMETER :: PSMOOTH=0.05_SP ! smoothing parameter REAL(SP) :: W_FUNC ! result from smoother ! ------------------------------------------------------------------------------------------------- @@ -44,6 +46,9 @@ SUBROUTINE Q_MISSCELL_DIFF(fuseStruct) ) ! (associate) ! --------------------------------------------------------------------------------------- + ! check the need to compute flux derivatives + comp_dflux = .false.; if(present(want_dflux)) comp_dflux = want_dflux + ! --------------------------------------------------------------------------------------- SELECT CASE(SMODL%iARCH1) CASE(iopt_tension2_1) ! tension storage sub-divided into recharge and excess @@ -70,8 +75,12 @@ SUBROUTINE Q_MISSCELL_DIFF(fuseStruct) M_FLUX%RCHR2EXCS = 0._SP M_FLUX%TENS2FREE_1 = 0._SP ! compute over-flow of free water - W_FUNC = SMOOTHER(TSTATE%WATR_1,MPARAM%MAXWATR_1,PSMOOTH) - M_FLUX%OFLOW_1 = W_FUNC * (M_FLUX%EFF_PPT - M_FLUX%QSURF) + if(SMODL%iQSURF == iopt_arno_x_vic)then + M_FLUX%OFLOW_1 = 0._sp ! no need for overflow since the vic parmaeterization is smoothed now + else + W_FUNC = SMOOTHER(TSTATE%WATR_1,MPARAM%MAXWATR_1,PSMOOTH) + M_FLUX%OFLOW_1 = W_FUNC * (M_FLUX%EFF_PPT - M_FLUX%QSURF) + endif CASE DEFAULT print *, "SMODL%iARCH1 must be iopt_tension2_1, iopt_tension1_1, or iopt_onestate_1" STOP diff --git a/build/FUSE_SRC/physics/qinterflow_diff.f90 b/build/FUSE_SRC/physics/qinterflow_diff.f90 index 9b1ed32..221c2c5 100644 --- a/build/FUSE_SRC/physics/qinterflow_diff.f90 +++ b/build/FUSE_SRC/physics/qinterflow_diff.f90 @@ -7,7 +7,7 @@ module QINTERFLOW_DIFF_module contains - SUBROUTINE QINTERFLOW_DIFF(fuseStruct) + SUBROUTINE QINTERFLOW_DIFF(fuseStruct, want_dflux) ! --------------------------------------------------------------------------------------- ! Creator: ! -------- @@ -25,6 +25,9 @@ SUBROUTINE QINTERFLOW_DIFF(fuseStruct) IMPLICIT NONE ! input-output type(parent), intent(inout) :: fuseStruct ! parent fuse data structure + logical(lgt), intent(in), optional :: want_dflux ! if we want flux derivatives + ! internal + logical(lgt) :: comp_dflux ! flag to compute flux derivatives ! ------------------------------------------------------------------------------------------------- ! associate variables with elements of data structure associate(& @@ -33,6 +36,10 @@ SUBROUTINE QINTERFLOW_DIFF(fuseStruct) MPARAM => fuseStruct%param_adjust , & ! adjustable model parameters DPARAM => fuseStruct%param_derive & ! derived model parameters ) ! (associate) + ! ------------------------------------------------------------------------------------------------- + + ! check the need to compute flux derivatives + comp_dflux = .false.; if(present(want_dflux)) comp_dflux = want_dflux ! --------------------------------------------------------------------------------------- SELECT CASE(SMODL%iQINTF) diff --git a/build/FUSE_SRC/physics/qpercolate_diff.f90 b/build/FUSE_SRC/physics/qpercolate_diff.f90 index 9ff599c..8e19db5 100644 --- a/build/FUSE_SRC/physics/qpercolate_diff.f90 +++ b/build/FUSE_SRC/physics/qpercolate_diff.f90 @@ -7,7 +7,7 @@ module QPERCOLATE_DIFF_module contains - SUBROUTINE QPERCOLATE_DIFF(fuseStruct) + SUBROUTINE QPERCOLATE_DIFF(fuseStruct, want_dflux) ! --------------------------------------------------------------------------------------- ! Creator: ! -------- @@ -21,37 +21,92 @@ SUBROUTINE QPERCOLATE_DIFF(fuseStruct) USE nrtype ! variable types, etc. USE data_types, only: parent ! fuse parent data type USE model_defn ! model definition structure - USE model_defnames + USE model_defnames ! model definition names + use smoothers, only : sfrac, dsfrac ! smoothed fraction, derivative IMPLICIT NONE ! input-output type(parent), intent(inout) :: fuseStruct ! parent fuse data structure + logical(lgt), intent(in), optional :: want_dflux ! if we want flux derivatives ! internal + logical(lgt) :: comp_dflux ! flag to compute flux derivatives + integer(i4b) :: iState ! state index + real(sp) :: phi ! smoothed fraction of free water + real(sp) :: dphi_dx ! derivative in smoothed fraction of free water + real(sp) :: df_dpsi ! derivative of flux w.r.t. fraction + real(sp) :: dqperc_dx ! derivative of percolation fux w.r.t. water state REAL(SP) :: LZ_PD ! lower zone percolation demand ! --------------------------------------------------------------------------------------- ! associate variables with elements of data structure associate(& M_FLUX => fuseStruct%flux , & ! fluxes + dfx_dS => fuseStruct%df_dS , & ! deriv in fluxes w.r.t. states TSTATE => fuseStruct%state1 , & ! trial state variables (end of step) MPARAM => fuseStruct%param_adjust , & ! adjustable model parameters DPARAM => fuseStruct%param_derive & ! derived model parameters ) ! (associate) ! --------------------------------------------------------------------------------------- + ! check the need to compute flux derivatives + comp_dflux = .false.; if(present(want_dflux)) comp_dflux = want_dflux + ! --------------------------------------------------------------------------------------- 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 + + ! -------------------------------------------------------------------------------------- + ! upper zone control + ! -------------------------------------------------------------------------------------- + CASE(iopt_perc_w2sat, iopt_perc_f2sat) + + ! short-cuts + associate(k=>MPARAM%PERCRTE, c=>MPARAM%PERCEXP) + + ! compute fractions + select case(SMODL%iQPERC) + case(iopt_perc_w2sat); phi = sfrac(TSTATE%WATR_1, MPARAM%MAXWATR_1) + case(iopt_perc_f2sat); phi = sfrac(TSTATE%FREE_1, DPARAM%MAXFREE_1) + end select ! no need for default since already in block + + ! ----- compute flux ---------------------------------------------------------------- + M_FLUX%QPERC_12 = k*phi**c + + ! ----- compute derivative ---------------------------------------------------------- + if(comp_dflux)then + + ! compute derivative in the fractions + select case(SMODL%iQPERC) + case(iopt_perc_w2sat); dphi_dx = dsfrac(TSTATE%WATR_1, MPARAM%MAXWATR_1) + case(iopt_perc_f2sat); dphi_dx = dsfrac(TSTATE%FREE_1, DPARAM%MAXFREE_1) + end select ! no need for default since already in block + + ! compute derivatives in the percolation flux + df_dpsi = k*c*phi**(c - 1._sp) ! derivative of flux w.r.t. fraction + dqperc_dx = df_dpsi*dphi_dx + + ! populate derivative vector + do iState=1,nState + select case(cState(iState)%iSNAME) + case (iopt_FREE_1); dfx_dS(iState)%QPERC_12 = dqperc_dx ! exists if separate free tank + case (iopt_WATR_1); dfx_dS(iState)%QPERC_12 = dqperc_dx ! exists if one state in the upper layer + end select ! no default needed + end do ! looping through states + + endif ! if computing derivatives + + end associate + + ! -------------------------------------------------------------------------------------- + ! lower zone control + ! -------------------------------------------------------------------------------------- CASE(iopt_perc_lower) ! perc defined by moisture content in lower layer (SAC) - ! (compute lower-zone percolation demand -- multiplier on maximum percolation, then percolation) + + ! ----- compute flux ---------------------------------------------------------------- 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 + + ! ----- compute derivatives --------------------------------------------------------------------- + if(comp_dflux) stop "qpercolate: derivatives for iopt_perc_lower not implemented yet" + + CASE DEFAULT; stop "qpercolate: SMODL%iQPERC must be iopt_perc_f2sat, iopt_perc_w2sat, or iopt_perc_lower" END SELECT ! -------------------------------------------------------------------------------------- diff --git a/build/FUSE_SRC/physics/qsatexcess_diff.f90 b/build/FUSE_SRC/physics/qsatexcess_diff.f90 index fa454a2..753dee2 100644 --- a/build/FUSE_SRC/physics/qsatexcess_diff.f90 +++ b/build/FUSE_SRC/physics/qsatexcess_diff.f90 @@ -7,7 +7,7 @@ module QSATEXCESS_DIFF_MODULE contains - SUBROUTINE QSATEXCESS_DIFF(fuseStruct) + SUBROUTINE QSATEXCESS_DIFF(fuseStruct, want_dflux) ! ------------------------------------------------------------------------------------------------- ! Creator: ! -------- @@ -20,13 +20,21 @@ SUBROUTINE QSATEXCESS_DIFF(fuseStruct) ! ------------------------------------------------------------------------------------------------- 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 + USE nr, ONLY : gammp ! interface for the incomplete gamma function + USE smoothers, only : smax,dsmax ! smoothed max function, derivative IMPLICIT NONE ! input-output type(parent), intent(inout) :: fuseStruct ! parent fuse data structure - ! internal variables + logical(lgt), intent(in), optional :: want_dflux ! if we want flux derivatives + ! internal variables -- vic + real(sp) :: u,xp ! temporary variables + real(sp) :: ds_dx ! derivative of saturated area w.r.t. x + real(sp) :: dx_du ! derivative of smooth max(u,0) w.r.t. u + real(sp) :: du_dw ! derivative of u w.r.t. w + real(sp) :: ds_dw ! derivative of saturated area w.r.t. w + ! internal variables -- topmodel 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 @@ -34,16 +42,23 @@ SUBROUTINE QSATEXCESS_DIFF(fuseStruct) 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 + ! derivatives + logical(lgt) :: comp_dflux ! flag to compute flux derivatives + integer(i4b) :: iState ! state index ! ------------------------------------------------------------------------------------------------- ! associate variables with elements of data structure associate(& M_FLUX => fuseStruct%flux , & ! fluxes + dfx_dS => fuseStruct%df_dS , & ! deriv in fluxes w.r.t. states TSTATE => fuseStruct%state1 , & ! trial state variables (end of step) MPARAM => fuseStruct%param_adjust , & ! adjustable model parameters DPARAM => fuseStruct%param_derive & ! derived model parameters ) ! (associate) ! ------------------------------------------------------------------------------------------------- + ! check the need to compute flux derivatives + comp_dflux = .false.; if(present(want_dflux)) comp_dflux = want_dflux + ! saturated area method SELECT CASE(SMODL%iQSURF) @@ -52,9 +67,38 @@ SUBROUTINE QSATEXCESS_DIFF(fuseStruct) ! ------------------------------------------------------------------------------------------------ CASE(iopt_arno_x_vic) + ! define variables + associate(w=>TSTATE%WATR_1, wmax=>MPARAM%MAXWATR_1, b=>MPARAM%AXV_BEXP) + ! ----- compute flux ---------------------------------------------------------------------------- - M_FLUX%SATAREA = 1._sp - ( 1._sp - MIN(TSTATE%WATR_1/MPARAM%MAXWATR_1, 1._sp) )**MPARAM%AXV_BEXP - + u = 1._sp - w/wmax + xp = smax(u, 0._sp) ! smooth version of max(u,0) + M_FLUX%SATAREA = 1._sp - xp**b + + ! ----- compute derivatives --------------------------------------------------------------------- + if(comp_dflux)then + + ! compute derivative w.r.t. saturated area + ds_dx = -b*xp**(b - 1._sp) ! derivative of saturated area w.r.t. xp + dx_du = dsmax(u, 0._sp) ! derivative of smooth max(u,0) w.r.t. u + du_dw = -1._sp/wmax ! derivative of u w.r.t. w + ds_dw = du_dw*dx_du*ds_dx ! derivative of saturated area w.r.t. w + + ! since WATR_1 is the sum of individual state variables (e.g., WATR_1=TENS_1+FREE_1) simply copy derivative + do iState=1,nState + select case(cState(iState)%iSNAME) + case (iopt_TENS1A); dfx_dS(iState)%SATAREA = ds_dw ! exists if two tension tanks + case (iopt_TENS1B); dfx_dS(iState)%SATAREA = ds_dw ! exists if two tension tanks + case (iopt_TENS_1); dfx_dS(iState)%SATAREA = ds_dw ! exists if one tension tank + case (iopt_FREE_1); dfx_dS(iState)%SATAREA = ds_dw ! exists if separate free storage + case (iopt_WATR_1); dfx_dS(iState)%SATAREA = ds_dw ! exists if one state in the upper layer + end select ! no default needed + end do ! looping through states + + endif ! if want derivatives + + end associate + ! ------------------------------------------------------------------------------------------------ ! ----- PRMS variant (fraction of upper tension storage) ----------------------------------------- ! ------------------------------------------------------------------------------------------------ @@ -63,6 +107,9 @@ SUBROUTINE QSATEXCESS_DIFF(fuseStruct) ! ----- compute flux ---------------------------------------------------------------------------- M_FLUX%SATAREA = MIN(TSTATE%TENS_1/DPARAM%MAXTENS_1, 1._sp) * MPARAM%SAREAMAX + ! ----- compute derivatives --------------------------------------------------------------------- + if(comp_dflux) stop "qsatexcess: derivatives for iopt_prms_varnt not implemented yet" + ! ------------------------------------------------------------------------------------------------ ! ----- TOPMODEL parameterization (only valid for TOPMODEL qb) ----------------------------------- ! ------------------------------------------------------------------------------------------------ @@ -87,6 +134,9 @@ SUBROUTINE QSATEXCESS_DIFF(fuseStruct) M_FLUX%SATAREA = 1._sp - GAMMP(TI_SHP, TI_ARG) ! GAMMP is the incomplete Gamma function ENDIF + ! ----- compute derivatives --------------------------------------------------------------------- + if(comp_dflux) stop "qsatexcess: derivatives for iopt_tmdl_param not implemented yet" + ! ------------------------------------------------------------------------------------------------ ! ------------------------------------------------------------------------------------------------ ! check processed surface runoff selection @@ -95,7 +145,7 @@ SUBROUTINE QSATEXCESS_DIFF(fuseStruct) STOP END SELECT ! (different surface runoff options) - + ! ...and, compute surface runoff ! ------------------------------ M_FLUX%QSURF = M_FLUX%EFF_PPT * M_FLUX%SATAREA diff --git a/build/FUSE_SRC/physics/smoothers.f90 b/build/FUSE_SRC/physics/smoothers.f90 index 71f5277..246ed59 100644 --- a/build/FUSE_SRC/physics/smoothers.f90 +++ b/build/FUSE_SRC/physics/smoothers.f90 @@ -5,9 +5,116 @@ module smoothers private public:: LOGISMOOTH public:: smoother + public:: smax,dsmax + public:: sfrac,dsfrac contains + ! --------------------------------------------------------------------------------------- + ! --------------------------------------------------------------------------------------- + + PURE FUNCTION sfrac(x,xmax) result(xf) + ! --------------------------------------------------------------------------------------- + ! Creator: + ! -------- + ! Martyn Clark, 2025 + ! --------------------------------------------------------------------------------------- + ! Purpose: + ! -------- + ! Use smoothed min function to compute smooth fraction + ! --------------------------------------------------------------------------------------- + USE nrtype + implicit none + real(sp), intent(in) :: x ! x value + real(sp), intent(in) :: xmax ! maximum value + real(sp) :: xp ! smooth min(x,xmax) + real(sp) :: xf ! smooth fraction x/xmax + xp = xmax - smax(xmax - x, 0._sp) ! smooth version of min(x, xmax) + xf = max(0._sp, xp) / xmax ! use max(0._sp, xp) to account for small neg values at zero + end function sfrac + + ! --------------------------------------------------------------------------------------- + ! --------------------------------------------------------------------------------------- + + PURE FUNCTION dsfrac(x,xmax) result(dxf_dx) + ! --------------------------------------------------------------------------------------- + ! Creator: + ! -------- + ! Martyn Clark, 2025 + ! --------------------------------------------------------------------------------------- + ! Purpose: + ! -------- + ! Get derivative of the smooth fraction + ! --------------------------------------------------------------------------------------- + USE nrtype + implicit none + real(sp), intent(in) :: x ! x value + real(sp), intent(in) :: xmax ! maximum value + real(sp) :: dxp_dx ! derivative of the max smoother + real(sp) :: dxf_dx ! derivative of the smoothed fraction + ! NOTE: ignore the hard clamp at zero (very small differences and not worth the extra expense) + dxp_dx = dsmax(xmax - x, 0._sp) ! note signs cancel out + dxf_dx = dxp_dx / xmax + end function dsfrac + + ! --------------------------------------------------------------------------------------- + ! --------------------------------------------------------------------------------------- + + PURE FUNCTION smax(x,xmin) result(xp) + ! --------------------------------------------------------------------------------------- + ! Creator: + ! -------- + ! Martyn Clark, 2025 + ! --------------------------------------------------------------------------------------- + ! Purpose: + ! -------- + ! Compute smoothed max function following Kavetski and Kuczera (2007) + ! + ! Kavetski, D. and Kuczera, G., 2007. Model smoothing strategies to remove microscale + ! discontinuities and spurious secondary optima in objective functions in hydrological + ! calibration. Water Resources Research, 43(3). + ! --------------------------------------------------------------------------------------- + USE nrtype + implicit none + real(sp), intent(in) :: x ! x value + real(sp), intent(in) :: xmin ! minimum value + real(sp), parameter :: ms=1.e-4_sp ! smoothing parameter + real(sp) :: srt ! sqrt(x*x + ms) + real(sp) :: xp ! smooth max(x,xmin) + srt = sqrt((x-xmin)**2 + ms) + xp = 0.5_sp*(x + xmin + srt) ! smooth max(x,xmin) + end function smax + + ! --------------------------------------------------------------------------------------- + ! --------------------------------------------------------------------------------------- + + PURE FUNCTION dsmax(x,xmin) result(dxp) + ! --------------------------------------------------------------------------------------- + ! Creator: + ! -------- + ! Martyn Clark, 2025 + ! --------------------------------------------------------------------------------------- + ! Purpose: + ! -------- + ! Compute derivative of smoothed max function of Kavetski and Kuczera (2007) + ! + ! Kavetski, D. and Kuczera, G., 2007. Model smoothing strategies to remove microscale + ! discontinuities and spurious secondary optima in objective functions in hydrological + ! calibration. Water Resources Research, 43(3). + ! --------------------------------------------------------------------------------------- + USE nrtype + implicit none + real(sp), intent(in) :: x ! x value + real(sp), intent(in) :: xmin ! minimum value + real(sp), parameter :: ms=1.e-4_sp ! smoothing parameter + real(sp) :: u ! x-xmin + real(sp) :: srt ! sqrt(x*x + ms) + real(sp) :: dxp ! derivative of smooth max(x,xmin) + u = x-xmin + srt = sqrt(u*u + ms) + dxp = 0.5_sp*(1._sp + u/srt) ! derivative of smooth max(x,xmin) + end function dsmax + ! --------------------------------------------------------------------------------------- ! ---------------------------------------------------------------------------------------