SUBROUTINE INTERP_FORCING(day_of_year,tstep_of_day,MODEL_STEP,I)

  ! UW Land Surface Hydrology Group implementation of SAC/SNOW17 model
  ! modified from NLDAS implementation
  ! author: Ted Bohn, tbohn@hydro.washington.edu

  ! INTERPOLATES FORCING DATA TO MODEL TIME STEP

  ! driverMod contains definitions for all global driver variables
  USE driverMod

  IMPLICIT NONE
      
  ! Define local variables
  INTEGER, INTENT(IN) :: day_of_year,tstep_of_day,MODEL_STEP,I
  REAL    :: GMT, LHOUR, BTIME, ETIME, MTIME, WEIGHT1, WEIGHT2
  INTEGER :: ZONE,J

  GMT=MODEL_STEP*MODEL_DT/3600.+ (tstep_of_day-1)*FORCING_DT/3600
  IF(GMT.GE.24)GMT=GMT-24
  CALL LOCALTIME (GMT,LON(X(I),Y(I)),LHOUR,ZONE)
  CALL COSZENITH(LON(X(I),Y(I)),LAT(X(I),Y(I)),LHOUR,ZONE,day_of_year,CZMODEL(I))

  CALL FINTERP(0,'I',0,SWdownnow(I),SWdownplus(I),0,MADTT,MODEL_STEP,DT_SOLDNB(I))
  CALL FINTERP(0,'I',0,LWdownnow(I),LWdownplus(I),0,MADTT,MODEL_STEP,DT_LWDN(I))
  CALL FINTERP(0,'I',0,Tairnow(I),Tairplus(I),0,MADTT,MODEL_STEP,DT_TAIR(I))
  CALL FINTERP(0,'I',0,Qairnow(I),Qairplus(I),0,MADTT,MODEL_STEP,DT_SPFH(I))
  CALL FINTERP(0,'I',0,Snowfnow(I),Snowfplus(I),0,MADTT,MODEL_STEP,DT_SNOW(I))
  CALL FINTERP(0,'I',0,Rainfnow(I),Rainfplus(I),0,MADTT,MODEL_STEP,DT_RAIN(I))
  CALL FINTERP(0,'I',0,PSurfnow(I),PSurfplus(I),0,MADTT,MODEL_STEP,DT_PSFC(I))
  CALL FINTERP(0,'I',0,Windnow(I),Windplus(I),0,MADTT,MODEL_STEP,DT_WIND(I))
  DO J = 1,NBANDS
    CALL FINTERP(0,'I',0,PotEvapnow(I,J),PotEvapplus(I,J),0,MADTT,MODEL_STEP,DT_EDMND(I,J))
  END DO

  BTIME=(tstep_of_day-1)*FORCING_DT/3600
  ETIME=tstep_of_day*FORCING_DT/3600
  IF(ETIME.GE.24)ETIME=ETIME-24
  MTIME=GMT

!  CALL ZTERP(1,LAT(X(I),Y(I)),LON(X(I),Y(I)),BTIME,ETIME,MTIME,day_of_year,WEIGHT1,WEIGHT2)
!  DT_SOLDN(I) = SWdownnow(I)*WEIGHT1 + SWdownplus(I)*WEIGHT2

!  ! Arbitrary cutoff scheme from Brian Cosgrove
!  IF (CZMODEL(I).LE.0.1.and.DT_SOLDN(I).GE.400) THEN
!    DT_SOLDN(I) = DT_SOLDNB(I)
!  ENDIF

  DT_SOLDN(I) = DT_SOLDNB(I)

  DT_PRCP(I)  = DT_RAIN(I)+DT_SNOW(I)

! This has been moved to inside the elevation band loop, since each band can
! have a different amount of snow cover
!  ! determine precip type
!  FFROZP(I) = 0.0
!  IF (DT_PRCP(I) .GT. 0.0) THEN
!    IF (DT_TAIR(I) .LE. TFREEZ) THEN
!      FFROZP(I) = 1.0
!    ENDIF
!  ENDIF
!
!  ! determine net incoming solar
!  IF (SNEQV(I) .EQ. 0.0) THEN
!    Albedo(I) = ALBEDO_D(I)
!  ELSE
!    Albedo(I) = ALBEDO_D(I) + SNCOVR(I)*(SNOALB(I)-ALBEDO_D(I))
!    IF (Albedo(I) .GT. SNOALB(I)) Albedo(I)=SNOALB(I)
!  ENDIF
!  DT_SOLNET(I) = DT_SOLDN(I)*(1.0-Albedo(I))

END