SUBROUTINE WRITE_OUTPUT(OUTPUT_STEP)

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

  ! WRITES VARS FOR 1 TIMESTEP TO THE OUTPUT FILES

  ! driverMod contains definitions of all global driver variables
  USE driverMod

  IMPLICIT NONE

  ! Define local variables
  INTEGER OUTPUT_STEP
  INTEGER varid
  INTEGER I,J,K,L
  INTEGER start,count
  REAL TIME_tmp
  LOGICAL LBAND
  CHARACTER(len=200) :: OUTASCIIFILE

  LBAND = .FALSE.
  start = OUTPUT_STEP
  count = 1

  ! Loop over output files

  DO K = 1,6

    ! Write Time,Timestep
    TIME_tmp = (OUTPUT_STEP-1)*OUTPUT_DT
    status = NF_INQ_VARID(OUT_NCIDS(K),'time',varid)
    status = NF_PUT_VARA_REAL(OUT_NCIDS(K),varid,start,count,TIME_tmp)
    status = NF_INQ_VARID(OUT_NCIDS(K),'timestp',varid)
    status = NF_PUT_VARA_INT(OUT_NCIDS(K),varid,start,count,OUTPUT_STEP-1)

    ! Initialize var index L
    L = 1

    IF (K == 1) THEN

!     ! Q.1) General energy balance components
!
!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,SWnet,2,COMP_OUTPUT,LBAND)
!
!      L = L + 1
!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,LWnet,2,COMP_OUTPUT,LBAND)
!
!      L = L + 1
!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,Qle,2,COMP_OUTPUT,LBAND)
!
!      L = L + 1
!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,Qh,2,COMP_OUTPUT,LBAND)
!
!      L = L + 1
!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,Qg,2,COMP_OUTPUT,LBAND)
!
!      L = L + 1
!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,Qf,2,COMP_OUTPUT,LBAND)
!
!      L = L + 1
!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,Qv,2,COMP_OUTPUT,LBAND)
!
!      L = L + 1
!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,Qa,2,COMP_OUTPUT,LBAND)
!
!      L = L + 1
!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,DelSurfHeat,2,COMP_OUTPUT,LBAND)
!
!      L = L + 1
!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,DelColdCont,2,COMP_OUTPUT,LBAND)
!
!
!    ELSE IF (K == 2) THEN

      ! Q.2) General water balance components

      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,Snowf,2,COMP_OUTPUT,LBAND)

      L = L + 1
      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,Rainf,2,COMP_OUTPUT,LBAND)

      L = L + 1
      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,Evap,2,COMP_OUTPUT,LBAND)

      L = L + 1
      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,Qs,2,COMP_OUTPUT,LBAND)

      L = L + 1
      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,Qsb,2,COMP_OUTPUT,LBAND)

      L = L + 1
      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,Qsm,2,COMP_OUTPUT,LBAND)

      L = L + 1
      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,Qfz,2,COMP_OUTPUT,LBAND)

      L = L + 1
      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,Qst,2,COMP_OUTPUT,LBAND)

      L = L + 1
      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,DelSoilMoist,2,COMP_OUTPUT,LBAND)

      L = L + 1
      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,DelSWE,2,COMP_OUTPUT,LBAND)

      L = L + 1
      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,DelSurfstor,2,COMP_OUTPUT,LBAND)

      L = L + 1
      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,DelIntercept,2,COMP_OUTPUT,LBAND)

!    ELSE IF (K == 3) THEN

      ! Q.3) Surface state variables

!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,SnowT,2,COMP_OUTPUT,LBAND)
!
!      L = L + 1
!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,VegT,2,COMP_OUTPUT,LBAND)
!
!      L = L + 1
!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,BaresoilT,2,COMP_OUTPUT,LBAND)
!
      L = L + 1
      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,AvgSurfT,2,COMP_OUTPUT,LBAND)

      L = L + 1
      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,RadT,2,COMP_OUTPUT,LBAND)

!      L = L + 1
!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,Albedo,2,COMP_OUTPUT,LBAND)

      L = L + 1
      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,SWE,2,COMP_OUTPUT,LBAND)

      L = L + 1
      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,SWEVeg,2,COMP_OUTPUT,LBAND)

      L = L + 1
      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,Surfstor,2,COMP_OUTPUT,LBAND)

!    ELSE IF (K == 4) THEN

      ! Q.4) Subsurface state variables

      L = L + 1
!!$      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,SoilMoist,3,COMP_OUTPUT,LBAND)
      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,SoilMoistTotal,2,COMP_OUTPUT,LBAND)

!      L = L + 1
!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,SoilTemp,3,COMP_OUTPUT,LBAND)
!
!      L = L + 1
!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,SMLiqFrac,3,COMP_OUTPUT,LBAND)
!
!      L = L + 1
!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,SMFrozFrac,3,COMP_OUTPUT,LBAND)

      L = L + 1
      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,SoilWet,2,COMP_OUTPUT,LBAND)

!    ELSE IF (K == 5) THEN

      ! Q.5) Evaporation components

      L = L + 1
      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,PotEvap,2,COMP_OUTPUT,LBAND)

!      L = L + 1
!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,Ecanop,2,COMP_OUTPUT,LBAND)
!
!      L = L + 1
!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,TVeg,2,COMP_OUTPUT,LBAND)
!
!      L = L + 1
!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,ESoil,2,COMP_OUTPUT,LBAND)
!
!      L = L + 1
!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,EWater,2,COMP_OUTPUT,LBAND)
!
!      L = L + 1
!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,RootMoist,2,COMP_OUTPUT,LBAND)
!
      L = L + 1
      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,CanopInt,2,COMP_OUTPUT,LBAND)
!
!      L = L + 1
!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,EvapSnow,2,COMP_OUTPUT,LBAND)
!
!      L = L + 1
!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,SubSnow,2,COMP_OUTPUT,LBAND)
!
!      L = L + 1
!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,SubSurf,2,COMP_OUTPUT,LBAND)
!
!      L = L + 1
!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,ACond,2,COMP_OUTPUT,LBAND)

!    ELSE IF (K == 6) THEN

      ! Q.7) Cold season processes

!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,SnowFrac,2,COMP_OUTPUT,LBAND)
!
!      L = L + 1
!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,Fdepth,2,COMP_OUTPUT,LBAND)
!
!      L = L + 1
!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,Tdepth,2,COMP_OUTPUT,LBAND)
!
!      L = L + 1
!      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,SAlbedo,2,COMP_OUTPUT,LBAND)
!
      L = L + 1
      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,SnowDepth,2,COMP_OUTPUT,LBAND)

      L = L + 1
      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,CanEvap,2,COMP_OUTPUT,LBAND)

      L = L + 1
      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,Transp,2,COMP_OUTPUT,LBAND)

      L = L + 1
      CALL WRITE_OUT_VAR(OUT_NCIDS(K),VARIDS(K,L),OUTPUT_STEP,BareEvap,2,COMP_OUTPUT,LBAND)

    END IF

  END DO

!!$DO I=1,landlen
!!$   IF (LONL(I) <= -100.0) THEN
!!$      WRITE(OUTASCIIFILE,'("/wb_",F5.2,"_",F7.2)')LATL(I),LONL(I)
!!$   ELSE
!!$       WRITE(OUTASCIIFILE,'("/wb_",F5.2,"_",F6.2)')LATL(I),LONL(I)
!!$    END IF
!!$   OUTASCIIFILE = TRIM(RESULT)//TRIM(OUTASCIIFILE)
!!$   OPEN(10,FILE=OUTASCIIFILE,STATUS='OLD',POSITION='APPEND')
!!$   WRITE(10,*) SoilMoistTotal(I),Evap(I),CanEvap(I),Transp(I),BareEvap(I)
!!$   CLOSE(10)
!!$END DO

END SUBROUTINE WRITE_OUTPUT