SUBROUTINE READ_LSC_NETCDF()

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

  ! READS LAND SURFACE CHARACTERISTICS

  ! driverMod contains definitions for all global driver variables
  USE driverMod

  IMPLICIT NONE

  ! Define local variables
  INTEGER ndims,natts
  INTEGER varid
  CHARACTER*20 name
  INTEGER xtype, dimids(2)
  INTEGER J,I,K,NT,start2d(2),count2d(2),start3d(3),count3d(3),start,count
  REAL    TEMP(xlen,ylen)
  INTEGER TEMP_INT(xlen,ylen)
  INTEGER land_idx
  REAL    depth_save(landlen)
  REAL    z_in_m, term1, term2

  ! Get Col
  status = NF_INQ_VARID(LSC_NCID,'col',varid)
  status = NF_GET_VAR_INT(LSC_NCID,varid,COL)

  ! Get Row
  status = NF_INQ_VARID(LSC_NCID,'row',varid)
  status = NF_GET_VAR_INT(LSC_NCID,varid,ROW)

!!$  ! Get Lon
!!$  status = NF_INQ_VARID(LSC_NCID,'lon',varid)
!!$  status = NF_GET_VAR_REAL(LSC_NCID,varid,LON)
!!$
!!$  ! Get Lat
!!$  status = NF_INQ_VARID(LSC_NCID,'lat',varid)
!!$  status = NF_GET_VAR_REAL(LSC_NCID,varid,LAT)

  ! Get Cellid
  status = NF_INQ_VARID(LSC_NCID,'CellID',varid)
  status = NF_GET_VAR_INT(LSC_NCID,varid,CELLID)

  ! Get Catchment soil parameters

  status = NF_INQ_VARID(LSC_NCID,'bee',varid)
  status = NF_GET_VAR_REAL(LSC_NCID,varid,TEMP)
  land_idx = 1
  DO I = 1,ylen
     DO J = 1,xlen
        IF (LANDMASK(J,I) == 1) THEN
           CP(land_idx)%bee = TEMP(J,I)
           land_idx = land_idx + 1
        END IF
     END DO
  END DO

  status = NF_INQ_VARID(LSC_NCID,'psis',varid)
  status = NF_GET_VAR_REAL(LSC_NCID,varid,TEMP)
  land_idx = 1
  DO I = 1,ylen
     DO J = 1,xlen
        IF (LANDMASK(J,I) == 1) THEN
           CP(land_idx)%psis = TEMP(J,I)
           land_idx = land_idx + 1
        END IF
     END DO
  END DO

  status = NF_INQ_VARID(LSC_NCID,'poros',varid)
  status = NF_GET_VAR_REAL(LSC_NCID,varid,TEMP)
  land_idx = 1
  DO I = 1,ylen
     DO J = 1,xlen
        IF (LANDMASK(J,I) == 1) THEN
           CP(land_idx)%poros = TEMP(J,I)
           land_idx = land_idx + 1
        END IF
     END DO
  END DO

  status = NF_INQ_VARID(LSC_NCID,'cond',varid)
  status = NF_GET_VAR_REAL(LSC_NCID,varid,TEMP)
  land_idx = 1
  DO I = 1,ylen
     DO J = 1,xlen
        IF (LANDMASK(J,I) == 1) THEN
           CP(land_idx)%cond = TEMP(J,I)
           land_idx = land_idx + 1
        END IF
     END DO
  END DO

  status = NF_INQ_VARID(LSC_NCID,'wpwet',varid)
  status = NF_GET_VAR_REAL(LSC_NCID,varid,TEMP)
  land_idx = 1
  DO I = 1,ylen
     DO J = 1,xlen
        IF (LANDMASK(J,I) == 1) THEN
           CP(land_idx)%wpwet = TEMP(J,I)
           land_idx = land_idx + 1
        END IF
     END DO
  END DO

  status = NF_INQ_VARID(LSC_NCID,'dpth',varid)
  status = NF_GET_VAR_REAL(LSC_NCID,varid,TEMP)
  land_idx = 1
  DO I = 1,ylen
     DO J = 1,xlen
        IF (LANDMASK(J,I) == 1) THEN
           CP(land_idx)%dpth = TEMP(J,I)
           land_idx = land_idx + 1
        END IF
     END DO
  END DO

  status = NF_INQ_VARID(LSC_NCID,'atau',varid)
  status = NF_GET_VAR_REAL(LSC_NCID,varid,TEMP)
  land_idx = 1
  DO I = 1,ylen
     DO J = 1,xlen
        IF (LANDMASK(J,I) == 1) THEN
           CP(land_idx)%atau = TEMP(J,I)
           land_idx = land_idx + 1
        END IF
     END DO
  END DO

  status = NF_INQ_VARID(LSC_NCID,'btau',varid)
  status = NF_GET_VAR_REAL(LSC_NCID,varid,TEMP)
  land_idx = 1
  DO I = 1,ylen
     DO J = 1,xlen
        IF (LANDMASK(J,I) == 1) THEN
           CP(land_idx)%btau = TEMP(J,I)
           land_idx = land_idx + 1
        END IF
     END DO
  END DO

  status = NF_INQ_VARID(LSC_NCID,'gnu',varid)
  status = NF_GET_VAR_REAL(LSC_NCID,varid,TEMP)
  land_idx = 1
  DO I = 1,ylen
     DO J = 1,xlen
        IF (LANDMASK(J,I) == 1) THEN
           CP(land_idx)%gnu = TEMP(J,I)
           land_idx = land_idx + 1
        END IF
     END DO
  END DO

  status = NF_INQ_VARID(LSC_NCID,'ars1',varid)
  status = NF_GET_VAR_REAL(LSC_NCID,varid,TEMP)
  land_idx = 1
  DO I = 1,ylen
     DO J = 1,xlen
        IF (LANDMASK(J,I) == 1) THEN
           CP(land_idx)%ars1 = TEMP(J,I)
           land_idx = land_idx + 1
        END IF
     END DO
  END DO

  status = NF_INQ_VARID(LSC_NCID,'ars2',varid)
  status = NF_GET_VAR_REAL(LSC_NCID,varid,TEMP)
  land_idx = 1
  DO I = 1,ylen
     DO J = 1,xlen
        IF (LANDMASK(J,I) == 1) THEN
           CP(land_idx)%ars2 = TEMP(J,I)
           land_idx = land_idx + 1
        END IF
     END DO
  END DO

  status = NF_INQ_VARID(LSC_NCID,'ars3',varid)
  status = NF_GET_VAR_REAL(LSC_NCID,varid,TEMP)
  land_idx = 1
  DO I = 1,ylen
     DO J = 1,xlen
        IF (LANDMASK(J,I) == 1) THEN
           CP(land_idx)%ars3 = TEMP(J,I)
           land_idx = land_idx + 1
        END IF
     END DO
  END DO

  status = NF_INQ_VARID(LSC_NCID,'ara1',varid)
  status = NF_GET_VAR_REAL(LSC_NCID,varid,TEMP)
  land_idx = 1
  DO I = 1,ylen
     DO J = 1,xlen
        IF (LANDMASK(J,I) == 1) THEN
           CP(land_idx)%ara1 = TEMP(J,I)
           land_idx = land_idx + 1
        END IF
     END DO
  END DO

  status = NF_INQ_VARID(LSC_NCID,'ara2',varid)
  status = NF_GET_VAR_REAL(LSC_NCID,varid,TEMP)
  land_idx = 1
  DO I = 1,ylen
     DO J = 1,xlen
        IF (LANDMASK(J,I) == 1) THEN
           CP(land_idx)%ara2 = TEMP(J,I)
           land_idx = land_idx + 1
        END IF
     END DO
  END DO

  status = NF_INQ_VARID(LSC_NCID,'ara3',varid)
  status = NF_GET_VAR_REAL(LSC_NCID,varid,TEMP)
  land_idx = 1
  DO I = 1,ylen
     DO J = 1,xlen
        IF (LANDMASK(J,I) == 1) THEN
           CP(land_idx)%ara3 = TEMP(J,I)
           land_idx = land_idx + 1
        END IF
     END DO
  END DO

  status = NF_INQ_VARID(LSC_NCID,'ara4',varid)
  status = NF_GET_VAR_REAL(LSC_NCID,varid,TEMP)
  land_idx = 1
  DO I = 1,ylen
     DO J = 1,xlen
        IF (LANDMASK(J,I) == 1) THEN
           CP(land_idx)%ara4 = TEMP(J,I)
           land_idx = land_idx + 1
        END IF
     END DO
  END DO

  status = NF_INQ_VARID(LSC_NCID,'arw1',varid)
  status = NF_GET_VAR_REAL(LSC_NCID,varid,TEMP)
  land_idx = 1
  DO I = 1,ylen
     DO J = 1,xlen
        IF (LANDMASK(J,I) == 1) THEN
           CP(land_idx)%arw1 = TEMP(J,I)
           land_idx = land_idx + 1
        END IF
     END DO
  END DO

  status = NF_INQ_VARID(LSC_NCID,'arw2',varid)
  status = NF_GET_VAR_REAL(LSC_NCID,varid,TEMP)
  land_idx = 1
  DO I = 1,ylen
     DO J = 1,xlen
        IF (LANDMASK(J,I) == 1) THEN
           CP(land_idx)%arw2 = TEMP(J,I)
           land_idx = land_idx + 1
        END IF
     END DO
  END DO

  status = NF_INQ_VARID(LSC_NCID,'arw3',varid)
  status = NF_GET_VAR_REAL(LSC_NCID,varid,TEMP)
  land_idx = 1
  DO I = 1,ylen
     DO J = 1,xlen
        IF (LANDMASK(J,I) == 1) THEN
           CP(land_idx)%arw3 = TEMP(J,I)
           land_idx = land_idx + 1
        END IF
     END DO
  END DO

  status = NF_INQ_VARID(LSC_NCID,'arw4',varid)
  status = NF_GET_VAR_REAL(LSC_NCID,varid,TEMP)
  land_idx = 1
  DO I = 1,ylen
     DO J = 1,xlen
        IF (LANDMASK(J,I) == 1) THEN
           CP(land_idx)%arw4 = TEMP(J,I)
           land_idx = land_idx + 1
        END IF
     END DO
  END DO

  status = NF_INQ_VARID(LSC_NCID,'bf1',varid)
  status = NF_GET_VAR_REAL(LSC_NCID,varid,TEMP)
  land_idx = 1
  DO I = 1,ylen
     DO J = 1,xlen
        IF (LANDMASK(J,I) == 1) THEN
           CP(land_idx)%bf1 = TEMP(J,I)
           land_idx = land_idx + 1
        END IF
     END DO
  END DO

  status = NF_INQ_VARID(LSC_NCID,'bf2',varid)
  status = NF_GET_VAR_REAL(LSC_NCID,varid,TEMP)
  land_idx = 1
  DO I = 1,ylen
     DO J = 1,xlen
        IF (LANDMASK(J,I) == 1) THEN
           CP(land_idx)%bf2 = TEMP(J,I)
           land_idx = land_idx + 1
        END IF
     END DO
  END DO

  status = NF_INQ_VARID(LSC_NCID,'bf3',varid)
  status = NF_GET_VAR_REAL(LSC_NCID,varid,TEMP)
  land_idx = 1
  DO I = 1,ylen
     DO J = 1,xlen
        IF (LANDMASK(J,I) == 1) THEN
           CP(land_idx)%bf3 = TEMP(J,I)
           land_idx = land_idx + 1
        END IF
     END DO
  END DO

  status = NF_INQ_VARID(LSC_NCID,'tsa1',varid)
  status = NF_GET_VAR_REAL(LSC_NCID,varid,TEMP)
  land_idx = 1
  DO I = 1,ylen
     DO J = 1,xlen
        IF (LANDMASK(J,I) == 1) THEN
           CP(land_idx)%tsa1 = TEMP(J,I)
           land_idx = land_idx + 1
        END IF
     END DO
  END DO

  status = NF_INQ_VARID(LSC_NCID,'tsa2',varid)
  status = NF_GET_VAR_REAL(LSC_NCID,varid,TEMP)
  land_idx = 1
  DO I = 1,ylen
     DO J = 1,xlen
        IF (LANDMASK(J,I) == 1) THEN
           CP(land_idx)%tsa2 = TEMP(J,I)
           land_idx = land_idx + 1
        END IF
     END DO
  END DO

  status = NF_INQ_VARID(LSC_NCID,'tsb1',varid)
  status = NF_GET_VAR_REAL(LSC_NCID,varid,TEMP)
  land_idx = 1
  DO I = 1,ylen
     DO J = 1,xlen
        IF (LANDMASK(J,I) == 1) THEN
           CP(land_idx)%tsb1 = TEMP(J,I)
           land_idx = land_idx + 1
        END IF
     END DO
  END DO

  status = NF_INQ_VARID(LSC_NCID,'tsb2',varid)
  status = NF_GET_VAR_REAL(LSC_NCID,varid,TEMP)
  land_idx = 1
  DO I = 1,ylen
     DO J = 1,xlen
        IF (LANDMASK(J,I) == 1) THEN
           CP(land_idx)%tsb2 = TEMP(J,I)
           land_idx = land_idx + 1
        END IF
     END DO
  END DO

  ! Three soil depths for soil moisture model: 
  !
  ! dzsf: surface layer
  ! dzrz: root zone        -> water capacity of the root zone
  ! dzpr: unsaturated zone -> approx depth-to-bedrock
  !
  ! NOTE: Units of dz** are [mm] while excess/deficits from catchment()
  !       are in SI units (ie kg/m^2) or loosely speaking, in mm of water.
  !       In other words, density of water (1000 kg/m^3) is built 
  !       into dz** (reichle, 5 Feb 04).

  DO K =1,land_idx

     CP(K)%dzsf = 20.
     CP(K)%dzrz = 1000.

     ! changed re-setting of dzrz back to earlier value because
     ! Sarith's parameters are in fact consistent that the earlier version
     ! reichle, 12 Sep 2007
     !
     ! cp(k)%dzpr = max(1500., cp(k)%dpth)
     !
     ! previously, root zone depth ranged from .75m to 1m, which
     ! is inconsistent with subroutine catchment(), where root
     ! zone depth is hard-wired to 1m, and with the time scale
     ! parameters, that have been derived for 1m root zone depth
     ! (THE LATTER IS IN FACT *NOT* TRUE - reichle, 12 Sep 2007)
     ! - reichle, 30 May 2003

     CP(K)%dzpr = max(1000., CP(K)%dpth)
     IF (CP(K)%dzrz > 0.75*CP(K)%dzpr) CP(K)%dzrz = 0.75*CP(K)%dzpr

     ! soil storages
     CP(K)%vgwmax = CP(K)%poros*CP(K)%dzrz

     z_in_m = CP(K)%dzpr/1000.

     term1 = -1.+((CP(K)%psis-z_in_m)/CP(K)%psis)**((CP(K)%bee-1.)/CP(K)%bee)

     term2 = CP(K)%psis*CP(K)%bee/(CP(K)%bee-1)

     CP(K)%cdcr1 = 1000.*CP(K)%poros*(z_in_m-(-term2*term1))

     CP(K)%cdcr2 = (1.-CP(K)%wpwet)*CP(K)%poros*CP(K)%dzpr

     ! soil depths for ground temperature model

!!$     if (N_gndtmp/=6) then
!!$        write (*,*) 'read_land_parameters: using N_gndtmp = ', N_gndtmp
!!$        write (*,*) 'but only 6 layer depths are specified, STOPPING!'
!!$        stop
!!$     end if

     CP(K)%dzgt(1) =  0.0988  
     CP(K)%dzgt(2) =  0.1952
     CP(K)%dzgt(3) =  0.3859
     CP(K)%dzgt(4) =  0.7626
     CP(K)%dzgt(5) =  1.5071
     CP(K)%dzgt(6) = 10.0000
  END DO

    
END SUBROUTINE READ_LSC_NETCDF