SUBROUTINE OPEN_OUTPUT()

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

  ! CREATES OUTPUT FILES

  ! driverMod contains definitions of all global driver variables
  USE driverMod

  IMPLICIT NONE

  ! Define local variables
  INTEGER nvars,ngatts,natts,unlimited
  INTEGER xid,yid,landid,levelid,tstepid
  INTEGER rowvarid,colvarid,latvarid,lonvarid,cellidvarid,landmaskvarid,varid
  CHARACTER*20 name
  INTEGER xtype, dimids(2),odim(2),outdim_2d(3),outdim_3d(4),ndims,ndims_2d,ndims_3d
  INTEGER I,J,K,L,NT
  INTEGER :: temp_int(landlen)
  REAL :: temp(landlen)
  INTEGER :: land_idx
  CHARACTER(len=200) :: OUTASCIIFILE

  IF (COMP_OUTPUT) THEN
    ndims = 1
    ndims_2d = 2
    ndims_3d = 3
  ELSE
    ndims = 2
    ndims_2d = 3
    ndims_3d = 4
  END IF

  DO K = 1,1

    status = NF_CREATE(OUTFILES(K),0,OUT_NCIDS(K))
    status = NF_DEF_DIM(OUT_NCIDS(K),'tstep',NF_UNLIMITED,tstepid)
    IF (K == 1) THEN
      status = NF_DEF_DIM(OUT_NCIDS(K),'level',MAXNSOIL,levelid)
    END IF
    IF (COMP_OUTPUT) THEN
      status = NF_DEF_DIM(OUT_NCIDS(K),'land',landlen,landid)
    ELSE
      status = NF_DEF_DIM(OUT_NCIDS(K),'y',ylen,yid)
      status = NF_DEF_DIM(OUT_NCIDS(K),'x',xlen,xid)
    END IF

    IF (COMP_OUTPUT) THEN
      odim(1) = landid
      outdim_2d(1) = landid
      outdim_2d(2) = tstepid
      outdim_3d(1) = landid
      outdim_3d(2) = levelid
      outdim_3d(3) = tstepid
    ELSE
      odim(1) = xid
      odim(2) = yid
      outdim_2d(1) = xid
      outdim_2d(2) = yid
      outdim_2d(3) = tstepid
      outdim_3d(1) = xid
      outdim_3d(2) = yid
      outdim_3d(3) = levelid
      outdim_3d(4) = tstepid
    END IF

    status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),NF_GLOBAL,'Conventions',7,'GDT 1.3')
    status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),NF_GLOBAL,'file_name',31,OUTFILES(K))

    status = NF_DEF_VAR(OUT_NCIDS(K),'row',NF_INT,ndims,odim,rowvarid)
    status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),rowvarid,'long_name',8,'Grid Row')
    status = NF_PUT_ATT_INT(OUT_NCIDS(K),rowvarid,'missing_value', NF_INT,1,NODATA_INT)

    status = NF_DEF_VAR(OUT_NCIDS(K),'col',NF_INT,ndims,odim,colvarid)
    status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),colvarid,'long_name',11,'Grid Column')
    status = NF_PUT_ATT_INT(OUT_NCIDS(K),colvarid,'missing_value', NF_INT,1,NODATA_INT)

    status = NF_DEF_VAR(OUT_NCIDS(K),'nav_lat',NF_REAL,ndims,odim,latvarid)
    status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),latvarid,'units',13,'degrees_north')
    status = NF_PUT_ATT_REAL(OUT_NCIDS(K),latvarid,'valid_min',NF_REAL,1,-90.)
    status = NF_PUT_ATT_REAL(OUT_NCIDS(K),latvarid,'valid_max',NF_REAL,1,90.)
    status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),latvarid,'long_name',8,'Latitude')
    status = NF_PUT_ATT_REAL(OUT_NCIDS(K),latvarid,'missing_value', NF_REAL,1,NODATA)

    status = NF_DEF_VAR(OUT_NCIDS(K),'nav_lon',NF_REAL,ndims,odim,lonvarid)
    status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),lonvarid,'units',12,'degrees_east')
    status = NF_PUT_ATT_REAL(OUT_NCIDS(K),lonvarid,'valid_min',NF_REAL,1,-180.)
    status = NF_PUT_ATT_REAL(OUT_NCIDS(K),lonvarid,'valid_max',NF_REAL,1,180.)
    status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),lonvarid,'long_name',9,'Longitude')
    status = NF_PUT_ATT_REAL(OUT_NCIDS(K),lonvarid,'missing_value', NF_REAL,1,NODATA)

    status = NF_DEF_VAR(OUT_NCIDS(K),'land',NF_INT,ndims,odim,landmaskvarid)
    status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),landmaskvarid,'units',1,'-')
    status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),landmaskvarid,'long_name',9,'Land Mask')
    status = NF_PUT_ATT_INT(OUT_NCIDS(K),landmaskvarid,'missing_value', NF_INT,1,NODATA_INT)

    status = NF_DEF_VAR(OUT_NCIDS(K),'CellID',NF_INT,ndims,odim,cellidvarid)
    status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),cellidvarid,'long_name',7,'Cell ID')
    status = NF_PUT_ATT_INT(OUT_NCIDS(K),cellidvarid,'missing_value', NF_INT,1,NODATA_INT)

    status = NF_DEF_VAR(OUT_NCIDS(K),'time',NF_REAL,1,tstepid,varid)
    status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),varid,'units',7,'seconds')
    status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),varid,'calendar',9,'gregorian')
    status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),varid,'title',4,'Time')
    status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),varid,'long_name',9,'Time axis')
    status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),varid,'time_origin',19,MODEL_START_TIME)

    status = NF_DEF_VAR(OUT_NCIDS(K),'timestp',NF_INT,1,tstepid,varid)
    status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),varid,'units',9,'timesteps')
    status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),varid,'title',10,'Time steps')
    status = NF_PUT_ATT_REAL(OUT_NCIDS(K),varid,'tstep_sec',NF_REAL,1,OUTPUT_DT_REAL)
    status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),varid,'long_name',14,'Time step axis')
    status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),varid,'time_origin',19,MODEL_START_TIME)

    ! Initialize var index L
    L = 1

    IF (K == 1) THEN
 
      status = NF_DEF_VAR(OUT_NCIDS(K),'Snowf',NF_REAL,ndims_2d,outdim_2d, VARIDS(K,L))
      status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),VARIDS(K,L),'units',7,'kg/m^2s')
      status = NF_PUT_ATT_REAL(OUT_NCIDS(K),VARIDS(K,L),'missing_value', NF_REAL,1,NODATA)

      L = L + 1
      status = NF_DEF_VAR(OUT_NCIDS(K),'Rainf',NF_REAL,ndims_2d,outdim_2d, VARIDS(K,L))
      status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),VARIDS(K,L),'units',7,'kg/m^2s')
      status = NF_PUT_ATT_REAL(OUT_NCIDS(K),VARIDS(K,L),'missing_value', NF_REAL,1,NODATA)

      L = L + 1
      status = NF_DEF_VAR(OUT_NCIDS(K),'Evap',NF_REAL,ndims_2d,outdim_2d,VARIDS(K,L))
      status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),VARIDS(K,L),'units',7,'kg/m^2s')
      status = NF_PUT_ATT_REAL(OUT_NCIDS(K),VARIDS(K,L),'missing_value', NF_REAL,1,NODATA)

      L = L + 1
      status = NF_DEF_VAR(OUT_NCIDS(K),'Qs',NF_REAL,ndims_2d,outdim_2d,VARIDS(K,L))
      status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),VARIDS(K,L),'units',7,'kg/m^2s')
      status = NF_PUT_ATT_REAL(OUT_NCIDS(K),VARIDS(K,L),'missing_value', NF_REAL,1,NODATA)

      L = L + 1
      status = NF_DEF_VAR(OUT_NCIDS(K),'Qsb',NF_REAL,ndims_2d,outdim_2d,VARIDS(K,L))
      status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),VARIDS(K,L),'units',7,'kg/m^2s')
      status = NF_PUT_ATT_REAL(OUT_NCIDS(K),VARIDS(K,L),'missing_value', NF_REAL,1,NODATA)

      L = L + 1
      status = NF_DEF_VAR(OUT_NCIDS(K),'Qsm',NF_REAL,ndims_2d,outdim_2d,VARIDS(K,L))
      status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),VARIDS(K,L),'units',7,'kg/m^2s')
      status = NF_PUT_ATT_REAL(OUT_NCIDS(K),VARIDS(K,L),'missing_value', NF_REAL,1,NODATA)

      L = L + 1
      status = NF_DEF_VAR(OUT_NCIDS(K),'Qfz',NF_REAL,ndims_2d,outdim_2d,VARIDS(K,L))
      status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),VARIDS(K,L),'units',7,'kg/m^2s')
      status = NF_PUT_ATT_REAL(OUT_NCIDS(K),VARIDS(K,L),'missing_value', NF_REAL,1,NODATA)

      L = L + 1
      status = NF_DEF_VAR(OUT_NCIDS(K),'Qst',NF_REAL,ndims_2d,outdim_2d,VARIDS(K,L))
      status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),VARIDS(K,L),'units',7,'kg/m^2s')
      status = NF_PUT_ATT_REAL(OUT_NCIDS(K),VARIDS(K,L),'missing_value', NF_REAL,1,NODATA)

      L = L + 1
      status = NF_DEF_VAR(OUT_NCIDS(K),'DelSoilMoist',NF_REAL,ndims_2d,outdim_2d, VARIDS(K,L))
      status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),VARIDS(K,L),'units',6, 'kg/m^2')
      status = NF_PUT_ATT_REAL(OUT_NCIDS(K),VARIDS(K,L),'missing_value', NF_REAL,1,NODATA)

      L = L + 1
      status = NF_DEF_VAR(OUT_NCIDS(K),'DelSWE',NF_REAL,ndims_2d,outdim_2d,VARIDS(K,L))
      status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),VARIDS(K,L),'units',6,'kg/m^2')
      status = NF_PUT_ATT_REAL(OUT_NCIDS(K),VARIDS(K,L),'missing_value', NF_REAL,1,NODATA)

      L = L + 1
      status = NF_DEF_VAR(OUT_NCIDS(K),'DelSurfStor',NF_REAL,ndims_2d,outdim_2d, VARIDS(K,L))
      status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),VARIDS(K,L),'units',6, 'kg/m^2')
      status = NF_PUT_ATT_REAL(OUT_NCIDS(K),VARIDS(K,L),'missing_value', NF_REAL,1,NODATA)

      L = L + 1
      status = NF_DEF_VAR(OUT_NCIDS(K),'DelIntercept',NF_REAL,ndims_2d,outdim_2d, VARIDS(K,L))
      status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),VARIDS(K,L),'units',6,'kg/m^2')
      status = NF_PUT_ATT_REAL(OUT_NCIDS(K),VARIDS(K,L),'missing_value', NF_REAL,1,NODATA)

      L = L + 1
      status = NF_DEF_VAR(OUT_NCIDS(K),'AvgSurfT',NF_REAL,ndims_2d,outdim_2d, VARIDS(K,L))
      status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),VARIDS(K,L),'units',1,'K')
      status = NF_PUT_ATT_REAL(OUT_NCIDS(K),VARIDS(K,L),'missing_value', NF_REAL,1,NODATA)

      L = L + 1
      status = NF_DEF_VAR(OUT_NCIDS(K),'RadT',NF_REAL,ndims_2d,outdim_2d, VARIDS(K,L))
      status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),VARIDS(K,L),'units',1,'K')
      status = NF_PUT_ATT_REAL(OUT_NCIDS(K),VARIDS(K,L),'missing_value', NF_REAL,1,NODATA)

      L = L + 1
      status = NF_DEF_VAR(OUT_NCIDS(K),'SWE',NF_REAL,ndims_2d,outdim_2d, VARIDS(K,L))
      status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),VARIDS(K,L),'units',6,'kg/m^2')
      status = NF_PUT_ATT_REAL(OUT_NCIDS(K),VARIDS(K,L),'missing_value', NF_REAL,1,NODATA)

      L = L + 1
      status = NF_DEF_VAR(OUT_NCIDS(K),'SWEVeg',NF_REAL,ndims_2d,outdim_2d, VARIDS(K,L))
      status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),VARIDS(K,L),'units',6,'kg/m^2')
      status = NF_PUT_ATT_REAL(OUT_NCIDS(K),VARIDS(K,L),'missing_value', NF_REAL,1,NODATA)

      L = L + 1
      status = NF_DEF_VAR(OUT_NCIDS(K),'SurfStor',NF_REAL,ndims_2d,outdim_2d, VARIDS(K,L))
      status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),VARIDS(K,L),'units',6,'kg/m^2')
      status = NF_PUT_ATT_REAL(OUT_NCIDS(K),VARIDS(K,L),'missing_value', NF_REAL,1,NODATA)

!!$      L = L + 1
!!$      status = NF_DEF_VAR(OUT_NCIDS(K),'SoilMoist',NF_REAL,ndims_3d,outdim_3d, VARIDS(K,L))
!!$      status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),VARIDS(K,L),'units',6,'kg/m^2')
!!$      status = NF_PUT_ATT_REAL(OUT_NCIDS(K),VARIDS(K,L),'missing_value', NF_REAL,1,NODATA)
      L = L + 1
      status = NF_DEF_VAR(OUT_NCIDS(K),'SoilMoistTotal',NF_REAL,ndims_2d,outdim_2d, VARIDS(K,L))
      status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),VARIDS(K,L),'units',6,'kg/m^2')
      status = NF_PUT_ATT_REAL(OUT_NCIDS(K),VARIDS(K,L),'missing_value', NF_REAL,1,NODATA)

      L = L + 1
      status = NF_DEF_VAR(OUT_NCIDS(K),'SoilWet',NF_REAL,ndims_2d,outdim_2d, VARIDS(K,L))
      status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),VARIDS(K,L),'units',1,'-')
      status = NF_PUT_ATT_REAL(OUT_NCIDS(K),VARIDS(K,L),'missing_value', NF_REAL,1,NODATA)

      L = L + 1
      status = NF_DEF_VAR(OUT_NCIDS(K),'PotEvap',NF_REAL,ndims_2d,outdim_2d, VARIDS(K,L))
      status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),VARIDS(K,L),'units',7,'kg/m^2s')
      status = NF_PUT_ATT_REAL(OUT_NCIDS(K),VARIDS(K,L),'missing_value', NF_REAL,1,NODATA)

      L = L + 1
      status = NF_DEF_VAR(OUT_NCIDS(K),'CanopInt',NF_REAL,ndims_2d,outdim_2d, VARIDS(K,L))
      status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),VARIDS(K,L),'units',6,'kg/m^2')
      status = NF_PUT_ATT_REAL(OUT_NCIDS(K),VARIDS(K,L),'missing_value', NF_REAL,1,NODATA)

!      L = L + 1
!      status = NF_DEF_VAR(OUT_NCIDS(K),'SubSnow',NF_REAL,ndims_2d,outdim_2d,VARIDS(K,L))
!      status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),VARIDS(K,L),'units',7,'kg/m^2s')
!      status = NF_PUT_ATT_REAL(OUT_NCIDS(K),VARIDS(K,L),'missing_value', NF_REAL,1,NODATA)

      L = L + 1
      status = NF_DEF_VAR(OUT_NCIDS(K),'SnowDepth',NF_REAL,ndims_2d,outdim_2d, VARIDS(K,L))
      status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),VARIDS(K,L),'units',1,'m')
      status = NF_PUT_ATT_REAL(OUT_NCIDS(K),VARIDS(K,L),'missing_value', NF_REAL,1,NODATA)

      L = L + 1
      status = NF_DEF_VAR(OUT_NCIDS(K),'CanEvap',NF_REAL,ndims_2d,outdim_2d, VARIDS(K,L))
      status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),VARIDS(K,L),'units',7,'kg/m^2s')
      status = NF_PUT_ATT_REAL(OUT_NCIDS(K),VARIDS(K,L),'missing_value', NF_REAL,1,NODATA)

      L = L + 1
      status = NF_DEF_VAR(OUT_NCIDS(K),'Transp',NF_REAL,ndims_2d,outdim_2d, VARIDS(K,L))
      status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),VARIDS(K,L),'units',7,'kg/m^2s')
      status = NF_PUT_ATT_REAL(OUT_NCIDS(K),VARIDS(K,L),'missing_value', NF_REAL,1,NODATA)

      L = L + 1
      status = NF_DEF_VAR(OUT_NCIDS(K),'BareEvap',NF_REAL,ndims_2d,outdim_2d, VARIDS(K,L))
      status = NF_PUT_ATT_TEXT(OUT_NCIDS(K),VARIDS(K,L),'units',7,'kg/m^2s')
      status = NF_PUT_ATT_REAL(OUT_NCIDS(K),VARIDS(K,L),'missing_value', NF_REAL,1,NODATA)

    END IF

    status = NF_ENDDEF(OUT_NCIDS(K))

    ! Write row, col, etc.

    IF (COMP_OUTPUT) THEN
      land_idx = 0
      DO I = 1,ylen
        DO J = 1,xlen
          IF (LANDMASK(J,I) == 1) THEN
            land_idx = land_idx + 1
            temp_int(land_idx) = ROW(J,I)
          END IF
        END DO
      END DO
      status = NF_PUT_VAR_INT(OUT_NCIDS(K),rowvarid,temp_int)
      land_idx = 0
      DO I = 1,ylen
        DO J = 1,xlen
          IF (LANDMASK(J,I) == 1) THEN
            land_idx = land_idx + 1
            temp_int(land_idx) = COL(J,I)
          END IF
        END DO
      END DO
      status = NF_PUT_VAR_INT(OUT_NCIDS(K),colvarid,temp_int)
      land_idx = 0
      DO I = 1,ylen
        DO J = 1,xlen
          IF (LANDMASK(J,I) == 1) THEN
            land_idx = land_idx + 1
            temp(land_idx) = LON(J,I)
          END IF
        END DO
      END DO
      status = NF_PUT_VAR_REAL(OUT_NCIDS(K),lonvarid,temp)
      land_idx = 0
      DO I = 1,ylen
        DO J = 1,xlen
          IF (LANDMASK(J,I) == 1) THEN
            land_idx = land_idx + 1
            temp(land_idx) = LAT(J,I)
          END IF
        END DO
      END DO
      status = NF_PUT_VAR_REAL(OUT_NCIDS(K),latvarid,temp)
      land_idx = 0
      DO I = 1,ylen
        DO J = 1,xlen
          IF (LANDMASK(J,I) == 1) THEN
            land_idx = land_idx + 1
            temp_int(land_idx) = LANDMASK(J,I)
          END IF
        END DO
      END DO
      status = NF_PUT_VAR_INT(OUT_NCIDS(K),landmaskvarid,temp_int)
      land_idx = 0
      DO I = 1,ylen
        DO J = 1,xlen
          IF (LANDMASK(J,I) == 1) THEN
            land_idx = land_idx + 1
            temp_int(land_idx) = CELLID(J,I)
          END IF
        END DO
      END DO
      status = NF_PUT_VAR_INT(OUT_NCIDS(K),cellidvarid,temp_int)
    ELSE
      status = NF_PUT_VAR_INT(OUT_NCIDS(K),rowvarid,ROW)
      status = NF_PUT_VAR_INT(OUT_NCIDS(K),colvarid,COL)
      status = NF_PUT_VAR_REAL(OUT_NCIDS(K),lonvarid,LON)
      status = NF_PUT_VAR_REAL(OUT_NCIDS(K),latvarid,LAT)
      status = NF_PUT_VAR_INT(OUT_NCIDS(K),landmaskvarid,LANDMASK)
      status = NF_PUT_VAR_INT(OUT_NCIDS(K),cellidvarid,CELLID)
    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='NEW')
!!$     CLOSE(10)
!!$  END DO

END SUBROUTINE OPEN_OUTPUT