Feature/removeprint

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

build/FUSE_SRC/FUSE_DMSL/fuse_driver.f90

@@ -23,6 +23,7 @@ PROGRAM DISTRIBUTED_DRIVER
 ! data modules
 USE model_defn,nstateFUSE=>nstate                         ! model definition structures
 USE model_defnames                                        ! defines the integer model options
+USE globaldata, ONLY: isPrint                             ! flag for printing progress to screen
 USE multiforce, ONLY: forcefile,vname_aprecip             ! model forcing structures
 USE multiforce, ONLY: AFORCE, aValid                      ! time series of lumped forcing/response data
 USE multiforce, ONLY: nspat1, nspat2                      ! grid dimensions
@@ -43,10 +44,11 @@ PROGRAM DISTRIBUTED_DRIVER
 
 USE multiforce, only: ncid_forc                           ! NetCDF forcing file ID
 USE multiforce, only: ncid_var                            ! NetCDF forcing variable ID
-USE multistate, only: ncid_out                            ! NetCDF output file ID
+USE globaldata, only: ncid_out                            ! NetCDF output file ID
 
 USE multibands                                            ! basin band stuctures
-USE multiparam, ONLY: LPARAM, PARATT, NUMPAR              ! parameter metadata structures
+USE data_types, ONLY: PARATT                              ! data type for metadata
+USE multiparam, ONLY: LPARAM, NUMPAR                      ! parameter metadata structures
 USE multistate, only: gState                              ! gridded state variables
 USE multistate, only: gState_3d                           ! gridded state variables with a time dimension
 USE multiroute, ONLY: AROUTE                              ! model routing structures
@@ -345,9 +347,9 @@ PROGRAM DISTRIBUTED_DRIVER
 
   ! assign algorithmic control parameters for SCE
   ! convert characters to interger/MSP
-  READ (MAXN_STR,*) MAXN		 ! maximum number of trials before optimization is terminated
+  READ (MAXN_STR,*) MAXN     ! maximum number of trials before optimization is terminated
   READ (KSTOP_STR,*) KSTOP   ! number of shuffling loops the value must change by PCENTO (MAX=9)
-  READ (PCENTO_STR,*) PCENTO    ! the percentage
+  READ (PCENTO_STR,*) PCENTO ! the percentage
 
   PRINT *, 'SCE parameters read from file manager:'
   PRINT *, 'Maximum number of trials before SCE optimization is stopped (MAXN) = ', MAXN_STR
@@ -438,11 +440,13 @@ PROGRAM DISTRIBUTED_DRIVER
 
   FNAME_ASCII = TRIM(OUTPUT_PATH)//TRIM(dom_id)//'_'//TRIM(FMODEL_ID)//'_sce_output.txt'
 
+  ! turn off printing to screen
+  isPrint = .false.
+
   ! convert from SP used in FUSE to MSP used in SCE
   ALLOCATE(APAR_MSP(NUMPAR),BL_MSP(NUMPAR),BU_MSP(NUMPAR),URAND_MSP(NUMPAR))
 
   APAR_MSP=APAR
-  PRINT *, 'BL=',BL
   BL_MSP=BL
   BU_MSP=BU
   URAND_MSP=URAND

build/FUSE_SRC/FUSE_DMSL/fuse_rmse.f90

@@ -22,6 +22,9 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG
     ! data modules
     USE model_defn, ONLY:NSTATE,SMODL                        ! number of state variables
     USE model_defnames                                       ! integer model definitions
+    USE globaldata, ONLY: isPrint                            ! flag for printing progress to screen
+    USE globaldata, ONLY: fracstate0                         ! fraction of initial state (used for initialization)
+    USE globaldata, ONLY: NA_VALUE, NA_VALUE_SP              ! NA_VALUE for the forcing
     USE multiparam, ONLY: LPARAM,NUMPAR,MPARAM               ! list of model parameters
     USE multiforce, ONLY: MFORCE,AFORCE,DELTIM,ISTART        ! model forcing data
     USE multiforce, ONLY: numtim_in, itim_in                 ! length of input time series and associated index
@@ -34,9 +37,7 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG
     USE multiforce, ONLY:nspat1,nspat2                       ! spatial dimensions
     USE multiforce, ONLY:ncid_var                            ! NetCDF ID for forcing variables
     USE multiforce, ONLY:gForce,gForce_3d                    ! gridded forcing data
-    USE multistate, ONLY:fracstate0,TSTATE,MSTATE,FSTATE,&   ! model states
-         HSTATE                              ! model states (continued)
-    USE multiforce, ONLY:NA_VALUE, NA_VALUE_SP              ! NA_VALUE for the forcing
+    USE multistate, ONLY:TSTATE,MSTATE,FSTATE,HSTATE         ! model state variables 
     USE multistate, ONLY:gState,gState_3d                    ! gridded state variables
     USE multiroute, ONLY:MROUTE,AROUTE,AROUTE_3d             ! routed runoff
     USE multistats, ONLY:MSTATS,PCOUNT,MOD_IX                ! access model statistics; counter for param set
@@ -52,6 +53,11 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG
     USE str_2_xtry_module                                    ! provide access to the routine str_2_xtry
     USE xtry_2_str_module                                    ! provide access to the routine xtry_2_str
 
+    ! differentiable model
+    use data_types, only: parent                             ! fuse parent data type
+    use get_parent_module, only: get_parent                  ! populate the parent data structure
+    use implicit_solve_module, only:implicit_solve           ! simple implicit solve for differnetiable ODE
+
     ! interface blocks
     USE interfaceb, ONLY:ode_int,fuse_solve                  ! provide access to FUSE_SOLVE through ODE_INT
 
@@ -93,6 +99,10 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG
     CHARACTER(LEN=CLEN)                    :: CMESSAGE       ! error message of downwind routine
     INTEGER(I4B),PARAMETER::UNT=6  !1701 ! 6
 
+    ! differentiable model
+    type(parent)                           :: fuseStruct     ! parent fuse data structure
+
+
     ! ---------------------------------------------------------------------------------------
     ! allocate state vectors
     ALLOCATE(STATE0(NSTATE),STATE1(NSTATE),STAT=IERR)
@@ -107,24 +117,27 @@ 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
           CALL INIT_STATE(fracState0)             ! define FSTATE using fracState0
           gState_3d(iSpat1,iSpat2,1) = FSTATE     ! put the state into first time step of 3D structure
        END DO
     END DO
-    PRINT *, 'Model states initialized over the 2D gridded domain'
+    if(isPrint) PRINT *, 'Model states initialized over the 2D gridded domain'
 
     ! initialize elevations bands if snow module is on
-    PRINT *, 'N_BANDS =', N_BANDS
+    if(isPrint) PRINT *, 'N_BANDS =', N_BANDS
 
     IF (SMODL%iSNOWM.EQ.iopt_temp_index) THEN
       DO iSpat2=1,nSpat2
@@ -137,7 +150,7 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG
           END DO
         END DO
       END DO
-      PRINT *, 'Snow states initiatlized over the 2D gridded domain '
+      if(isPrint) PRINT *, 'Snow states initiatlized over the 2D gridded domain '
     ENDIF
 
     ! allocate 3d data structure for fluxes
@@ -177,10 +190,10 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG
         numtim_sub_cur=MIN(numtim_sub,numtim_sim-itim_sim+1)
 
         ! load forcing for desired period into gForce_3d
-        PRINT *, 'New subperiod: loading forcing for ',numtim_sub_cur,' time steps'
+        if(isPrint) PRINT *, 'New subperiod: loading forcing for ',numtim_sub_cur,' time steps'
         CALL get_gforce_3d(itim_in,numtim_sub_cur,ncid_forc,err,message)
         IF(err/=0)THEN; WRITE(*,*) 'Error while extracting 3d forcing'; STOP; ENDIF
-        PRINT *, 'Forcing loaded. Running FUSE...'
+        if(isPrint) PRINT *, 'Forcing loaded. Running FUSE...'
 
       ENDIF
 
@@ -245,9 +258,38 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG
                   RETURN
                END SELECT
 
-               ! temporally integrate the ordinary differential equations
-               CALL ODE_INT(FUSE_SOLVE,STATE0,STATE1,DT_SUB,DT_FULL,IERR,MESSAGE)
-               IF (IERR.NE.0) THEN; PRINT *, TRIM(MESSAGE); PAUSE; ENDIF
+              ! ----- start of soil physics code ------------------------------------------------------------
+
+              ! temporally integrate the ordinary differential equations
+              select case(diff_mode)
+
+                ! original code               
+                case(original)
+                 CALL ODE_INT(FUSE_SOLVE,STATE0,STATE1,DT_SUB,DT_FULL,IERR,MESSAGE)
+                 IF (IERR.NE.0) THEN; PRINT *, TRIM(MESSAGE); STOP 1; ENDIF
+
+                 !print*, state1
+                 !if(ITIM_IN > sim_beg+100) stop
+
+                ! differentiable code   
+                case(differentiable)
+
+                 ! populate parent fuse structure
+                 call get_parent(fuseStruct)
+
+                 ! solve differentiable ODEs
+                 call implicit_solve(fuseStruct, state0, state1, nState)
+                 !print*, state1
+                 !if(ITIM_IN > sim_beg+100) stop
+
+                 ! save fluxes
+                 W_FLUX = fuseStruct%flux
+
+                ! check options
+                case default; print*, "Cannot identify diff_mode"; stop 1
+              end select
+
+              ! ----- end of soil physics code --------------------------------------------------------------
 
               ! perform overland flow routing
               CALL Q_OVERLAND()
@@ -307,15 +349,15 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG
       ! if end of subperiod: write to output file and save states
       IF(itim_sub.EQ.numtim_sub_cur)THEN
 
-        PRINT *, 'End of subperiod reached:'
+        if(isPrint) PRINT *, 'End of subperiod reached:'
 
         ! write model output
         IF (OUTPUT_FLAG) THEN
-          PRINT *, 'Write output for ',numtim_sub_cur,' time steps starting at indice', itim_sim-numtim_sub_cur+1
+          if(isPrint) PRINT *, 'Write output for ',numtim_sub_cur,' time steps starting at indice', itim_sim-numtim_sub_cur+1
           CALL PUT_GOUTPUT_3D(itim_sim-numtim_sub_cur+1,itim_in-numtim_sub_cur+1,numtim_sub_cur,IPSET)
-          PRINT *, 'Done writing output'
+          if(isPrint) PRINT *, 'Done writing output'
         ELSE
-          PRINT *, 'OUTPUT_FLAG is set on FALSE, no output written'
+          if(isPrint) PRINT *, 'OUTPUT_FLAG is set on FALSE, no output written'
         END IF
 
         ! TODO: set gState_3d and MBANDS_VAR_4d to NA
@@ -344,20 +386,20 @@ SUBROUTINE FUSE_RMSE(XPAR,GRID_FLAG,NCID_FORC,RMSE,OUTPUT_FLAG,IPSET,MPARAM_FLAG
 
     ! get timing information
     CALL CPU_TIME(T2)
-    WRITE(*,*) "TIME ELAPSED = ", t2-t1
+    if(isPrint) WRITE(*,*) "TIME ELAPSED = ", t2-t1
 
     ! calculate mean summary statistics
     IF(.NOT.GRID_FLAG)THEN
 
-      PRINT *, 'Calculating performance metrics...'
+      if(isPrint) PRINT *, 'Calculating performance metrics...'
       CALL MEAN_STATS()
       RMSE = MSTATS%RAW_RMSE
 
+      print*, "NSE = ", MSTATS%NASH_SUTT
+
     ENDIF
 
-    PRINT *, 'Writing parameter values...'
-    CALL PUT_PARAMS(PCOUNT)
-    PRINT *, 'Writing model statistics...'
+    if(isPrint) PRINT *, 'Writing model statistics...'
     CALL PUT_SSTATS(PCOUNT)
 
     ! deallocate vectors

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

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 

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

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

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

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