SUBROUTINE READ_INITIAL()

  ! UW Land Surface Hydrology Group implementation of SAC/SNOW17 model
  ! modified from NLDAS implementation
  ! author: Kostas Andreadis, kostas@hydro.washington.edu

  ! READS INITIAL CONDITIONS

  ! driverMod contains definitions of all global driver variables
  USE driverMod

  IMPLICIT NONE

  ! Define local variables
  INTEGER n
  REAL tc, qa, ts, mc
  INTEGER ncid,ndims,nvars,ngatts,unlimited
  INTEGER xdimid,ydimid,zsoildimid,zsnowdimid,zsoillen_initial
  INTEGER xlen_initial,ylen_initial,NBANDS_initial,zsnowlen_initial
  INTEGER N_snow_initial,varid
  INTEGER J,I,K,L,start2d(2),count2d(2),start3d(3),count3d(3),start,count
  INTEGER landmask_initial(xlen,ylen)
  INTEGER land_idx
  REAL    TEMP(xlen,ylen)

  IF (INITIAL == '') THEN

     write(*,*)'No initial conditions file given; using default values'

     ! Set to default values
     tc = 280. 
     qa = .002 
     ts = 280. 
     mc = .3

     DO K=1,landlen
        CAT_PROGN(K)%tc1 = tc
        CAT_PROGN(K)%tc2 = tc
        CAT_PROGN(K)%tc4 = tc
        CAT_PROGN(K)%qa1 = qa
        CAT_PROGN(K)%qa2 = qa
        CAT_PROGN(K)%qa4 = qa
        CAT_PROGN(K)%capac = 0.0
        DO n=1,N_gndtmp
           IF (ts .lt. 273.16) THEN ! soil frozen
              CAT_PROGN(K)%ght(n) =                              &
                   (ts-273.16)*                                  &
                   ( 2.4e6*0.55*CP(K)%dzgt(n)             &
                   +   0.5*0.45*CP(K)%dzgt(n)*2.06e6 )    &
                   - 0.5*0.45*CP(K)%dzgt(n)*3.34e8
           ELSE ! soil thawed
              CAT_PROGN(K)%ght(n) =                              &
                   (ts-273.16)*                                  &
                   ( 2.4e6*0.55*CP(K)%dzgt(n)             &
                   +   0.5*0.45*CP(K)%dzgt(n)*4.185e6)
           END IF
        END DO
        DO n=1,N_snow
           CAT_PROGN(K)%wesn(n) = 0.0
           CAT_PROGN(K)%htsn(n) = 0.0
           CAT_PROGN(K)%sndz(n) = 0.0
        END DO
        CAT_PROGN(K)%catdef = max( (CP(K)%poros - mc), 0.05) * CP(K)%dzpr
        CAT_PROGN(K)%rzexc = 0.0
        CAT_PROGN(K)%srfexc = 0.0
     END DO

  ELSE

     ! Read initial conditions from INITIAL file

    WRITE (*,*) 'Reading initial conditions from ',trim(INITIAL)

    status = NF_OPEN(INITIAL, 0, ncid)
    IF (status .ne. NF_NOERR) THEN
      WRITE(*,*)'ERROR: cannot open initial condition file',INITIAL
      STOP
    END IF
    status = NF_INQ(ncid, ndims, nvars, ngatts, unlimited)
!    status = NF_INQ_DIMID(ncid,'z',zdimid)
!    status = NF_INQ_DIMID(ncid,'band',banddimid)
    status = NF_INQ_DIMID(ncid,'y',ydimid)
    status = NF_INQ_DIMID(ncid,'x',xdimid)
    status = NF_INQ_DIMID(ncid,'zsoil',zsoildimid)
    status = NF_INQ_DIMID(ncid,'zsnow',zsnowdimid)
!    status = NF_INQ_DIMLEN(ncid,banddimid,NBANDS_initial)
    status = NF_INQ_DIMLEN(ncid,ydimid,ylen_initial)
    status = NF_INQ_DIMLEN(ncid,xdimid,xlen_initial)
    status = NF_INQ_DIMLEN(ncid,zsoildimid,zsoillen_initial)
    status = NF_INQ_DIMLEN(ncid,zsnowdimid,zsnowlen_initial)

write(*,*)'N_gndtmp_initial',N_gndtmp
!write(*,*)'NBANDS_initial',NBANDS_initial
write(*,*)'ylen_initial',ylen_initial
write(*,*)'xlen_initial',xlen_initial

    ! Validate data dimensions
    IF (zsoillen_initial .gt. MAXNSOIL) THEN
      WRITE(*,*)'ERROR: MAXNSOIL',MAXNSOIL, &
        ' (from lsc file ',TRIM(LSC),') less than zsoillen_initial ',zsoillen_initial, &
        ' (from initial condition file ',TRIM(INITIAL),')'
      STOP
    END IF
    IF (zsnowlen_initial .gt. MAXNSNOW) THEN
      WRITE(*,*)'ERROR: MAXNSNOW',MAXNSNOW, &
        ' (from lsc file ',TRIM(LSC),') less than zsnowlen_initial ',zsnowlen_initial, &
        ' (from initial condition file ',TRIM(INITIAL),')'
      STOP
    END IF
    IF (ylen_initial /= ylen) THEN
      WRITE(*,*)'ERROR: ylen',ylen, &
        ' (from lsc file ',TRIM(LSC),') /= ylen ',ylen_initial, &
        ' (from initial condition file ',TRIM(INITIAL),')'
      STOP
    END IF
    IF (xlen_initial /= xlen) THEN
      WRITE(*,*)'ERROR: xlen',xlen, &
        ' (from lsc file ',TRIM(LSC),') /= xlen ',xlen_initial, &
        ' (from initial condition file ',TRIM(INITIAL),')'
      STOP
    END IF

    ! Get landmask and compare to landmask from LSC file
    status = NF_INQ_VARID(ncid,'land',varid)
    status = NF_GET_VAR_INT(ncid,varid,landmask_initial)
    DO I=1,ylen
      DO J=1,xlen
        IF (landmask_initial(J,I) /= LANDMASK(J,I)) THEN
          WRITE(*,*)'ERROR: landmask from lsc file ',TRIM(LSC), &
          ' /= landmask from initial condition file ',TRIM(INITIAL)
          STOP
        END IF
      END DO
    END DO

    ! Get Ground Heat Content -- removed bands from everything
    status = NF_INQ_VARID(ncid,'ght',varid)
    DO L = 1, MAXNSOIL
       start3d(1) = 1
       start3d(2) = 1
       start3d(3) = L
       count3d(1) = xlen
       count3d(2) = ylen
       count3d(3) = 1
       status = NF_GET_VARA_REAL(ncid,varid,start3d,count3d,TEMP)
       land_idx = 1
       DO I = 1, ylen
          DO J = 1, xlen
             IF (LANDMASK(J,I) == 1) THEN
                CAT_PROGN(land_idx)%ght(L) = TEMP(J,I)
                land_idx = land_idx + 1
             END IF
          END DO
       END DO
    END DO

   ! Get SWE -- removed bands from everything
    status = NF_INQ_VARID(ncid,'wesn',varid)
    DO L = 1, N_snow
       start3d(1) = 1
       start3d(2) = 1
       start3d(3) = L
       count3d(1) = xlen
       count3d(2) = ylen
       count3d(3) = 1
       status = NF_GET_VARA_REAL(ncid,varid,start3d,count3d,TEMP)
       land_idx = 1
       DO I = 1, ylen
          DO J = 1, xlen
             IF (LANDMASK(J,I) == 1) THEN
                CAT_PROGN(land_idx)%wesn(L) = TEMP(J,I)
                land_idx = land_idx + 1
             END IF
          END DO
       END DO
    END DO
 
   ! Get snow heat content -- removed bands from everything
    status = NF_INQ_VARID(ncid,'htsn',varid)
    DO L = 1, N_snow
       start3d(1) = 1
       start3d(2) = 1
       start3d(3) = L
       count3d(1) = xlen
       count3d(2) = ylen
       count3d(3) = 1
       status = NF_GET_VARA_REAL(ncid,varid,start3d,count3d,TEMP)
       land_idx = 1
       DO I = 1, ylen
          DO J = 1, xlen
             IF (LANDMASK(J,I) == 1) THEN
                CAT_PROGN(land_idx)%htsn(L) = TEMP(J,I)
                land_idx = land_idx + 1
             END IF
          END DO
       END DO
    END DO

   ! Get snow denpth -- removed bands from everything
    status = NF_INQ_VARID(ncid,'sndz',varid)
    DO L = 1, N_snow
       start3d(1) = 1
       start3d(2) = 1
       start3d(3) = L
       count3d(1) = xlen
       count3d(2) = ylen
       count3d(3) = 1
       status = NF_GET_VARA_REAL(ncid,varid,start3d,count3d,TEMP)
       land_idx = 1
       DO I = 1, ylen
          DO J = 1, xlen
             IF (LANDMASK(J,I) == 1) THEN
                CAT_PROGN(land_idx)%sndz(L) = TEMP(J,I)
                land_idx = land_idx + 1
             END IF
          END DO
       END DO
    END DO

   ! Get surface/canopy temperature -- removed bands from everything
    status = NF_INQ_VARID(ncid,'tc1',varid)
    start2d(1) = 1
    start2d(2) = 1
    count2d(1) = xlen
    count2d(2) = ylen
    status = NF_GET_VARA_REAL(ncid,varid,start2d,count2d,TEMP)
    land_idx = 1
    DO I = 1, ylen
       DO J = 1, xlen
          IF (LANDMASK(J,I) == 1) THEN
             CAT_PROGN(land_idx)%tc1 = TEMP(J,I)
             land_idx = land_idx + 1
          END IF
       END DO
    END DO
     WRITE(*,*) 'read tc1'
    
   ! Get tc2 surface/canopy temperature -- removed bands from everything
    status = NF_INQ_VARID(ncid,'tc2',varid)
    start2d(1) = 1
    start2d(2) = 1
    count2d(1) = xlen
    count2d(2) = ylen
    status = NF_GET_VARA_REAL(ncid,varid,start2d,count2d,TEMP)
    land_idx = 1
    DO I = 1, ylen
       DO J = 1, xlen
          IF (LANDMASK(J,I) == 1) THEN
             CAT_PROGN(land_idx)%tc2 = TEMP(J,I)
             land_idx = land_idx + 1
          END IF
       END DO
    END DO

   ! Get tc4 surface/canopy temperature -- removed bands from everything
    status = NF_INQ_VARID(ncid,'tc4',varid)
    start2d(1) = 1
    start2d(2) = 1
    count2d(1) = xlen
    count2d(2) = ylen
    status = NF_GET_VARA_REAL(ncid,varid,start2d,count2d,TEMP)
    land_idx = 1
    DO I = 1, ylen
       DO J = 1, xlen
          IF (LANDMASK(J,I) == 1) THEN
             CAT_PROGN(land_idx)%tc4 = TEMP(J,I)
             land_idx = land_idx + 1
          END IF
       END DO
    END DO

   ! Get qa1 specific humidity in canopy air -- removed bands from everything
    status = NF_INQ_VARID(ncid,'qa1',varid)
    start2d(1) = 1
    start2d(2) = 1
    count2d(1) = xlen
    count2d(2) = ylen
    status = NF_GET_VARA_REAL(ncid,varid,start2d,count2d,TEMP)
    land_idx = 1
    DO I = 1, ylen
       DO J = 1, xlen
          IF (LANDMASK(J,I) == 1) THEN
             CAT_PROGN(land_idx)%qa1 = TEMP(J,I)
             land_idx = land_idx + 1
          END IF
       END DO
    END DO
     WRITE(*,*) 'read qa1'
    
   ! Get qa2 specific humidity in canopy air -- removed bands from everything
    status = NF_INQ_VARID(ncid,'qa2',varid)
    start2d(1) = 1
    start2d(2) = 1
    count2d(1) = xlen
    count2d(2) = ylen
    status = NF_GET_VARA_REAL(ncid,varid,start2d,count2d,TEMP)
    land_idx = 1
    DO I = 1, ylen
       DO J = 1, xlen
          IF (LANDMASK(J,I) == 1) THEN
             CAT_PROGN(land_idx)%qa2 = TEMP(J,I)
             land_idx = land_idx + 1
          END IF
       END DO
    END DO

   ! Get qa4 specific humidity in canopy air -- removed bands from everything
    status = NF_INQ_VARID(ncid,'qa4',varid)
    start2d(1) = 1
    start2d(2) = 1
    count2d(1) = xlen
    count2d(2) = ylen
    status = NF_GET_VARA_REAL(ncid,varid,start2d,count2d,TEMP)
    land_idx = 1
    DO I = 1, ylen
       DO J = 1, xlen
          IF (LANDMASK(J,I) == 1) THEN
             CAT_PROGN(land_idx)%qa4 = TEMP(J,I)
             land_idx = land_idx + 1
          END IF
       END DO
    END DO

   ! Get capac canopy interception water -- removed bands from everything
    status = NF_INQ_VARID(ncid,'capac',varid)
    start2d(1) = 1
    start2d(2) = 1
    count2d(1) = xlen
    count2d(2) = ylen
    status = NF_GET_VARA_REAL(ncid,varid,start2d,count2d,TEMP)
    land_idx = 1
    DO I = 1, ylen
       DO J = 1, xlen
          IF (LANDMASK(J,I) == 1) THEN
             CAT_PROGN(land_idx)%capac = TEMP(J,I)
             land_idx = land_idx + 1
          END IF
       END DO
    END DO

   ! Get catdef catchment deficit -- removed bands from everything
    status = NF_INQ_VARID(ncid,'catdef',varid)
    start2d(1) = 1
    start2d(2) = 1
    count2d(1) = xlen
    count2d(2) = ylen
    status = NF_GET_VARA_REAL(ncid,varid,start2d,count2d,TEMP)
    land_idx = 1
    DO I = 1, ylen
       DO J = 1, xlen
          IF (LANDMASK(J,I) == 1) THEN
             CAT_PROGN(land_idx)%catdef = TEMP(J,I)
             land_idx = land_idx + 1
          END IF
       END DO
    END DO
   ! Get rzexc root zone excess -- removed bands from everything
    status = NF_INQ_VARID(ncid,'rzexc',varid)
    start2d(1) = 1
    start2d(2) = 1
    count2d(1) = xlen
    count2d(2) = ylen
    status = NF_GET_VARA_REAL(ncid,varid,start2d,count2d,TEMP)
    land_idx = 1
    DO I = 1, ylen
       DO J = 1, xlen
          IF (LANDMASK(J,I) == 1) THEN
             CAT_PROGN(land_idx)%rzexc = TEMP(J,I)
             land_idx = land_idx + 1
          END IF
       END DO
    END DO

   ! Get srfexc surface excess -- removed bands from everything
    status = NF_INQ_VARID(ncid,'srfexc',varid)
    start2d(1) = 1
    start2d(2) = 1
    count2d(1) = xlen
    count2d(2) = ylen
    status = NF_GET_VARA_REAL(ncid,varid,start2d,count2d,TEMP)
    land_idx = 1
    DO I = 1, ylen
       DO J = 1, xlen
          IF (LANDMASK(J,I) == 1) THEN
             CAT_PROGN(land_idx)%srfexc = TEMP(J,I)
             land_idx = land_idx + 1
          END IF
       END DO
    END DO

!     OPEN(10, file=INITIAL, form='unformatted', status='old', &
!          action='read')
!    IF (NBANDS_initial /= NBANDS) THEN
!      WRITE(*,*)'ERROR: NBANDS',NBANDS, &
!        ' (from lsc file ',TRIM(LSC),') /= NBANDS ',NBANDS_initial, &
!        ' (from initial condition file ',TRIM(INITIAL),')'
!      STOP
!    END IF
!     READ (10) (CAT_PROGN(n)%tc1,    n=1,landlen)
!     WRITE(*,*) 'read tc1'
!     READ (10) (CAT_PROGN(n)%tc2,    n=1,landlen)
!     READ (10) (CAT_PROGN(n)%tc4,    n=1,landlen)
!
!     READ (10) (CAT_PROGN(n)%qa1,    n=1,landlen)
!     READ (10) (CAT_PROGN(n)%qa2,    n=1,landlen)
!     READ (10) (CAT_PROGN(n)%qa4,    n=1,landlen)
!
!     READ (10) (CAT_PROGN(n)%capac,  n=1,landlen)
!
!     READ (10) (CAT_PROGN(n)%catdef, n=1,landlen)
!     READ (10) (CAT_PROGN(n)%rzexc,  n=1,landlen)
!     READ (10) (CAT_PROGN(n)%srfexc, n=1,landlen)
!
!     DO K=1,N_gndtmp
!        READ (10) (CAT_PROGN(n)%ght(k),  n=1,landlen)
!     end do
!
!     DO K=1,N_snow
!        READ (10) (CAT_PROGN(n)%wesn(k), n=1,landlen)
!     end do
!     DO K=1,N_snow
!        READ (10) (CAT_PROGN(n)%htsn(k), n=1,landlen)
!     end do
!     DO K=1,N_snow
!        READ (10) (CAT_PROGN(n)%sndz(k), n=1,landlen)
!     end do
!
!     CLOSE (10,status='keep')

  END IF

END SUBROUTINE READ_INITIAL