c**********************************************************************
c     subroutine erodin
c**********************************************************************
      subroutine erodin (
     i i_unit, o_unit)

c
c     +++ PURPOSE +++
c     Utility to read initial conditions and variables from
c     input file (stdin or eros.ipt) for the standalone erosion submodel
c
c     +++ ARGUMENT DECLARATIONS +++
c
      integer i_unit, o_unit
c
c     +++ ARGUMENT DEFINITIONS +++
c
c
c     +++ PARAMETERS +++
c
      integer mrcl
      parameter (mrcl = 132)
c
c     +++ ARGUMENT DECLARATIONS & +++
c     + + + GLOBAL COMMON BLOCKS + + +
*$noreference
      include 'p1werm.inc'
C     include 'p1unconv.inc'
      include 'p1const.inc'
      include 'm1sim.inc'
      include 'm1flag.inc'
      include 'm1subr.inc'
      include 'm1geo.inc'
      include 'b1glob.inc'
      include 's1layr.inc'
      include 's1phys.inc'
      include 's1dbh.inc'
      include 's1agg.inc'
      include 's1surf.inc'
      include 's1sgeo.inc'
      include 'h1db1.inc'
      include 'w1wind.inc'
      include 'w1pavg.inc'
*$reference
c
c     + + + LOCAL COMMON BLOCKS + + +
c
      integer debugflg
      common /flags/ debugflg
c
c     +++ LOCAL VARIABLES +++
      integer i,j,k
      integer x,y,s,b,a,l,h
      integer wflg
      real f(mntime), wfcalm, wuc, w0k, step
c
c     + + + LOCAL VARIABLE DEFINITIONS + + +
c     i, j, k = do-loop indices
c     x,y,s,b,a,l,h = do-loop indices
c     wflg = flag to determine format of wind speed data (0 = Weibull, 1 = real)
c     debugflg = flag to output debug data (0 = none, 1 = input, 2 = more, etc.)
c     f(mntime) = cumulative frequency of wind at speeds <awu(i)
c     wfcalm    = wind fraction intercept (+calm, - no calm in period)
c     wuc       = Weibull wind speed distribution scale factor (m/s)
c     w0k       = Weibull wind speed distribution shape factor
c     step      = tmp real variable for ntstep
c
c     +++ FUNCTIONS CALLED +++
c     getline
c
      character*(mrcl) getline
      character*(mrcl) line

c     +++ END SPECIFICATIONS +++
c
c      read  +++ MAIN  COMMON BLOCKS +++
c

c     +++ INIT STUFF +++

c     Read the debug flag to specify level of debug support
c     (currently only input file level debug supported)
      line = getline(i_unit)
      read (line,*) debugflg
C     read ((line=getline(i_unit)),*) debugflg

c     m1flag.inc

c     EROSION initialization flag (logical)
      line = getline(i_unit)
      read (line,*) am0eif
C     read (getline(i_unit),*) am0eif

c     EROSION "print" flag (integer)
      line = getline(i_unit)
      read (line,*) am0efl
C     read (getline(i_unit),*) am0efl

c     +++ SIMULATION REGION +++

C     Need to make sure there is no confusion on what is x1, x2, y1, and y2
C     Not clearly defined in include files

c     m1geo.inc

c     Simulation region diagonal corners (x1,y1) and (x2,y2)
      line = getline(i_unit)
      read (line,*) ((amxsim (x,y), x=1,2),y=1,2)
C     read (getline(i_unit),*) ((amxsim (x,y), x=1,2),y=1,2)
C     write (o_unit,*) 'sim x1,y1 x1,y2 ',  amxsim(1,1), amxsim(1,2)
C     write (o_unit,*) 'sim x2,y1 x2,y2 ',  amxsim(2,1), amxsim(2,2)

c     Simulation region orientation angle
      line = getline(i_unit)
      read (line,*) amasim
C     read (getline(i_unit),*) amasim


c     +++ ACCOUNTING REGIONS +++

c     m1geo.inc

c     Number of accounting regions
      line = getline(i_unit)
      read (line,*) nacctr
C     read (getline(i_unit),*) nacctr

c     Accounting Region diagonal corners (x1,y1) and (x2,y2)
      do 20 a = 1, nacctr
        line = getline(i_unit)
        read (line,*) ((amxar (x,y,a), x=1,2), y=1,2)
C       read (getline(i_unit),*) ((amxar (x,y,a), x=1,2), y=1,2)
   20 continue

c     +++ BARRIERS +++

c     m1geo.inc

c     Number of barriers
      line = getline(i_unit)
      read (line,*) nbr
C     read (getline(i_unit),*) nbr

C     NOTE: Barrier data must not be in the input file if "nbr = 0"

      do 30 b = 1, nbr
c       Barrier diagonal corners (x1,y1) and (x2,y2)
        line = getline(i_unit)
        read (line,*) ((amxbr(x,y,b), x=1,2), y=1,2)
C       read (getline(i_unit),*) ((amxbr(x,y,b), x=1,2), y=1,2)
c       Barrier ht, porosity, and width
        line = getline(i_unit)
        read (line,*) amzbr(b), ampbr(b), amxbrw(b)
C       read (getline(i_unit),*)  amzbr(b), ampbr(b), amxbrw(b)
   30 continue


c     +++ SUBREGION REGIONS +++

c     m1subr.inc

c     Number of subregions
      line = getline(i_unit)
      read (line,*) nsubr
C     read (getline(i_unit),*) nsubr

c     Dimensions, Biomass, Soil, and Hydrology (by subregion)
      do 100  s=1, nsubr
c       Subregion diagonal corners (x1,y1) and (x2,y2)
        line = getline(i_unit)
        read (line,*) ((amxsr (x,y,s), x=1,2),y=1,2)
C       read (getline(i_unit),*) ((amxsr (x,y,s), x=1,2),y=1,2)

c     +++ BIOMASS +++

c     b1glob.inc & b1geom.inc

c       Biomass height
        line = getline(i_unit)
        read (line,*)  abzht(s)
C       read (getline(i_unit),*) abzht(s)
c       Biomass stem area index and leaf area index
        line = getline(i_unit)
        read (line,*) abrsai(s), abrlai(s)
C       read (getline(i_unit),*) abrsai(s), abrlai(s)

C These aren't used in EROSION yet
*       Biomass stem area index by height
C       read (getline(i_unit),*) (abrsaz(h,s), h=1,mncz)
c       Biomass leaf area index by height
*       read (getline(i_unit),*) (abrlaz(h,s), h=1,mncz)

c       Biomass flat fraction cover, standing cover, and fraction total cover
*       read (getline(i_unit),*) abffcv(s), abfscv(s), abftcv(s)
C Only flat fraction cover used yet
        line = getline(i_unit)
        read (line,*) abffcv(s)
C       read (getline(i_unit),*) abffcv(s)

c     +++ SOIL +++

c       s1layr.inc, s1phys.inc & s1agg.inc

c       Number of soil layers (in this subregion)
        line = getline(i_unit)
        read (line,*) nslay(s)
C       read (getline(i_unit),*) nslay(s)

c       Soil layer thickness
        line = getline(i_unit)
        read (line,*) (aszlyt(l,s),l=1,nslay(s))
C       read (getline(i_unit),*) (aszlyt(l,s),l=1,nslay(s))
c       Soil layer bulk density
        line = getline(i_unit)
        read (line,*) (asdblk(l,s), l=1,nslay(s))
C       read (getline(i_unit),*) (asdblk(l,s), l=1,nslay(s))
c       Sand, silt, and clay fractions
        line = getline(i_unit)
        read (line,*) (asfsan(l,s), l=1,nslay(s))
C       read (getline(i_unit),*) (asfsan(l,s), l=1,nslay(s))
        line = getline(i_unit)
        read (line,*) (asfsil(l,s), l=1,nslay(s))
C       read (getline(i_unit),*) (asfsil(l,s), l=1,nslay(s))
        line = getline(i_unit)
        read (line,*) (asfcla(l,s), l=1,nslay(s))
C       read (getline(i_unit),*) (asfcla(l,s), l=1,nslay(s))
c       Volume of rock fraction
        line = getline(i_unit)
        read (line,*) (asvroc(l,s), l=1,nslay(s))
C
c       s1agg.inc
c       Soil layer aggregate density
        line = getline(i_unit)
        read (line,*) (asdagd(l,s), l=1,nslay(s))
C
c       Soil layer aggregate stability
        line = getline(i_unit)
        read (line,*) (aseags(l,s), l=1,nslay(s))

C Check these variables with ASD inc files and Hagen's EROSION inc files - LEW
c       Soil layer ASD parms (gmd, min, max, gsd)
        line = getline(i_unit)
        read (line,*) (aslagm(l,s), l=1,nslay(s))
C
        line = getline(i_unit)
        read (line,*) (aslagn(l,s), l=1,nslay(s))
C
        line = getline(i_unit)
        read (line,*) (aslagx(l,s), l=1,nslay(s))
C
        line = getline(i_unit)
        read (line,*) (as0ags(l,s), l=1,nslay(s))

c       s1surf.inc & s1sgeo.inc

c       Crust parms (fraction, thickness)
        line = getline(i_unit)
        read (line,*) asfcr(s), aszcr(s),
C       read (getline(i_unit),*) asfcr(s), aszcr(s),
c       Crust parms (fraction cover of loose material, mass loose material)
     &       asflos(s), asmlos(s),
c       Crust parms (crust density and stability)
     &       asdcr(s), asecr(s)

c       Random Roughness
        line = getline(i_unit)
        read (line,*) aslrr(s)
C       read (getline(i_unit),*) aslrr(s)
c       Oriented Roughness (ridge ht, spacing, width, orientation)
        line = getline(i_unit)
        read (line,*) aszrgh(s), asxrgs(s),
C       read (getline(i_unit),*) aszrgh(s), asxrgs(s),
     &       asxrgw(s), asargo(s)
c       Oriented Roughness ( spacing)
        line = getline(i_unit)
        read (line,*) asxdks(s)
C       read (getline(i_unit),*) asxdkh(s), asxdks(s)


c     +++ HYDROLOGY +++

c       h1db1.inc

c       Snow depth
        line = getline(i_unit)
        read (line,*) ahzsnd(s)
C       read (getline(i_unit),*) ahzsnd(s)

c       Soil layer wilting point
        line = getline(i_unit)
        read (line,*) (ahrwcw(l,s), l=1,nslay(s))
C       read (getline(i_unit),*) (ahrwcw(l,s), l=1,nslay(s))
c       Soil layer water content
        line = getline(i_unit)
        read (line,*) (ahrwca(l,s), l=1,nslay(s))
C       read (getline(i_unit),*) (ahrwca(l,s), l=1,nslay(s))
c       Soil surface hourly water content
        line = getline(i_unit)
        read (line,*) (ahrwc0(h,s), h=1,12)
C       read (getline(i_unit),*) (ahrwc0(h,s), h=1,12)
        line = getline(i_unit)
        read (line,*) (ahrwc0(h,s), h=13,24)
C       read (getline(i_unit),*) (ahrwc0(h,s), h=13,24)
  100 continue

c     +++ WEATHER +++

C We need to check on the units for air density - w1pavg.inc says (kg/m^3)
C Also, we need to see why it currently isn't being used - LJH said it was
c     Air density
      line = getline(i_unit)
      read (line,*) awdair
C     read (getline(i_unit),*) awdair

c     Wind Direction
      line = getline(i_unit)
      read (line,*) awadir
C     read (getline(i_unit),*) awadir

c     Number of "steps" during 24 hours (96 = 15 minute intervals)
      line = getline(i_unit)
      read (line,*) ntstep
C     read (getline(i_unit),*) ntstep

c     Weibull wind flag (0 - read Weibull parms, 1 - read wind speeds)
      line = getline(i_unit)
      read (line,*) wflg
C     read (getline(i_unit),*) wflg

c     w1wind.inc

c     wind data inputs as the Weibull paramters
c     (wfcalm, wuc, w0k) is indicated by code ntstep = 99
      if (wflg .eq. 0) then
c       Weibull parms (fraction calm, c, k)
        line = getline(i_unit)
        read (line,*) wfcalm, wuc, w0k
C       read (getline(i_unit),*) wfcalm, wuc, w0k

c       calculate daily max wind speed (99% speed)
C       awudmx = wuc*(-log((1.0-0.99)/(1-wfcalm)))**(1.0/w0k)

c       calculate period wind speeds
        step = ntstep
        do 198 i= 1, ntstep
c         find center of each step and add empirical last term from file ntstep.mcd
            f(i) = (1.0/(2.0*step)) + ((i-1)/step) +0.3/(step*w0k)
c         to prevent out-of-range
          if (f(i) .lt. wfcalm) then
            f(i) = wfcalm
          endif
          awu(i) = wuc*(-log((1.0-f(i))/(1.0-wfcalm)))**(1.0/w0k)
  198   continue
c      Use greatest interval wind speed rather than 99% speed above
       awudmx = awu(ntstep)

c     when (wflg .eq. 1) input wind period data directly
      else
c       if (mod(ntstep,6) .eq. 0) then
c         h = ntstep/6
c       else
c         h = ntstep/6+1
c       endif
C        write(o_unit,*) 'h,mod,ntstep/6: ',h, mod(ntstep,6), ntstep/6
        do 191 j = 1, ntstep/6
          line = getline(i_unit)
          read (line,*) (awu(i),i=(j-1)*6+1,(j-1)*6+6)
C         read (getline(i_unit),*) (awu(i),i=(j-1)*6+1,(j-1)*6+6)
191     continue
c       If not divisible evenly by 6, then get the remaining values
        if (mod(ntstep,6) .ne. 0) then
          line = getline(i_unit)
          read (line,*)
C         read (getline(i_unit),*)
     &      (awu(i),i=(j-1)*6+1,(j-1)*6+mod(ntstep,6))
        endif

c       Determine the maximum wind speed during the day
        awudmx = 0.0
        do 199 i = 1, ntstep
          if (awudmx .lt. awu(i)) then
            awudmx = awu(i)
          endif
  199   continue
      endif

c     + + + OUTPUT SECTION + + +
c
  200 format (1x, f6.2, 6f8.2)
  220 format (1x, 2f8.2)
  230 format (1x, 3f8.2)
  240 format (1x, 4f8.2)
  250 format (1x, 5f8.2)
  260 format (1x, 6f8.2)
  270 format (1x, 7f8.2)
  280 format (1x, 8f8.2)
  300 format (1x, 8i5)
  350 format (1x, i3, 6f8.2)
  400 format (1x, i3, 1x, l6, 4x, i3, 4x, i3, 4x, i3, 4x)

      write (o_unit,*)
     &  '      REPORT OF INPUTS (read by erodin.for) '
      write (o_unit,*)
      write (o_unit,*)  '+++ Control Flags, etc. +++'
      write (o_unit,*)
      write (o_unit,*) 'ntstep  am0eif  nsubr  nacctr  nbr'
      write (o_unit,400) ntstep, am0eif, nsubr, nacctr, nbr

      write (o_unit,*)
      write (o_unit,*)  '+++ SIMULATION REGION +++'
      write (o_unit,*)
      write (o_unit,*) 'orientation and dimensions of sim region'
      write (o_unit,*) 'amasim(deg)  amxsim - (x1,y1) (x2,y2)'
      write (o_unit,250) amasim,  ((amxsim(x,y), x=1,2), y=1,2)

      write (o_unit,*)
      write (o_unit,*)  '+++ ACCOUNTING REGIONS +++'
      write (o_unit,*)
      write (o_unit,*) 'nacctr - number of accounting regions'
      write (o_unit,*) nacctr
      write (o_unit,*) 'accounting region dimensions (x1,y1) (x2,y2)'
      do 1000 a = 1, nacctr
        write (o_unit,240) ((amxar (x,y,a), x=1,2), y=1,2)
 1000 continue

      write (o_unit,*)
      write (o_unit,*)  '+++ BARRIERS +++'
      write (o_unit,*)
      if (nbr .gt. 0) then
        write (o_unit,*) 'nbr - number of barriers'
        write (o_unit,*) nbr
        write (o_unit,*) 'barrier dim (x1,y1) (x2,y2) ',
     &                   'and height, porosity, and width'
        do 1010 b = 1, nbr
          write (o_unit,270) ((amxbr(x,y,b), x=1,2), y=1,2),
     &                       amzbr(b), ampbr(b), amxbrw(b)
 1010   continue
      else
         write (o_unit,*) ' no barriers'
      endif

      write (o_unit,*)
      write (o_unit,*)  '+++ SUBREGIONS +++'
      write (o_unit,*)
      write (o_unit,*) 'nsubr - number of subregions'
      write (o_unit,*) nsubr
      write (o_unit,*) 'subregion dimensions (x1,y1) (x2,y2)'
      do 1020 s = 1, nsubr
        write (o_unit,240) ((amxsr (i,j,s), i=1,2), j=1,2)
 1020 continue

      do 1100 s = 1, nsubr
        write (o_unit,*)
        write (o_unit,*)
        write(o_unit,*) '*********************** Subregion ', s,
     &                 ' ***********************'

        write (o_unit,*)
        write(o_unit,*) '+++ BIOMASS +++'
        write (o_unit,*)
        write (o_unit,*) 'Biomass height, SAI, LAI, and flat cover'
        write (o_unit,240) abzht(s), abrsai(s), abrlai(s), abffcv(s)

        write (o_unit,*)
        write (o_unit,*) '+++ SOIL +++ '
        write (o_unit,*)
        write (o_unit,*) 'nslay - number of soil layers'
        write (o_unit,*) nslay(s)

        write (o_unit,*)
        write (o_unit,*) 'layer  thickness  bulk density  ',
     &                   'sand   silt   clay    rock vol'
        do 1030 l = 1, nslay(s)
          write (o_unit,350) l, aszlyt(l,s), asdblk(l,s),
     &         asfsan(l,s), asfsil(l,s), asfcla(l,s), asvroc(l,s)
 1030   continue

        write (o_unit,*)
          write (o_unit,*) 'layer    AgD     AgS ',
     &                   ' GMD    GMDmn     GMDmx    GSD'
        do 1040 l = 1, nslay(s)
          write (o_unit,350) l, asdagd(l,s), aseags(l,s),
     &           aslagm(l,s),aslagn(l,s), aslagx(l,s), as0ags(l,s)
 1040   continue

        write (o_unit,*)
        write (o_unit,*) 'Crust fraction, mass of LOS, frac. LOS, ',
     &                   'density, and stability'
        write (o_unit,250)asfcr(s),asmlos(s),asflos(s),asdcr(s),asecr(s)

        write (o_unit,*)
        write (o_unit,*) 'RR, Ridge ht, width, spacing, ',
     &                   'orientation, &  dike spacing'
        write (o_unit,270) aslrr(s), aszrgh(s), asxrgw(s), asxrgs(s),
     &                     asargo(s),  asxdks(s)

        write (o_unit,*)
        write (o_unit,*) '+++ HYDROLOGY +++ '
        write (o_unit,*)

        write (o_unit,*) 'Snow depth (mm)'
        write (o_unit,*) ahzsnd(s)

        write (o_unit,*)
        write (o_unit,*) 'layer  wilting and actual water contents'
        do 1050 l = 1, nslay(s)
          write (o_unit,350) l, ahrwcw (l,s), ahrwca(l,s)
 1050   continue
        write (o_unit,*) 'Hourly water contents - ahrwc0'
        write (o_unit,260) (ahrwc0(h,s), h=1,6)
        write (o_unit,260) (ahrwc0(h,s), h=7,12)
        write (o_unit,260) (ahrwc0(h,s), h=13,18)
        write (o_unit,260) (ahrwc0(h,s), h=19,24)

 1100 continue

      write (o_unit,*)
      write(o_unit,*) '+++ WEATHER +++'
      write (o_unit,*)
      write (o_unit,*) ' wind dir (deg) and max wind speed (m/s)'
      write(o_unit,220)  awadir, awudmx

      write (o_unit,*)
      write (o_unit,*) 'Wind speeds (m/s) - ', ntstep,' intervals'
      k = 6
      j = 1
  860 if(k .lt. ntstep) then
        write (o_unit,260) (awu(i), i=j,k)
        j = k+1
        k = k+6
        go to 860
      else
        k = ntstep
        write (o_unit,260) (awu(i), i=j,k)
      endif
      write (o_unit,*)
      write (o_unit,*) 'END OF INPUTS'

      close (i_unit)
      return
      end

c**********************************************************************
      character*(*) function getline(i_unit)

      integer debugflg
      common /flags/ debugflg

      integer i_unit
      integer dataline, linecount

      save dataline, linecount

1     read (i_unit, '(A)') getline
      linecount = linecount + 1
      if (BTEST(debugflg,0)) then
        write (6, *) linecount, ': ', getline
      endif

      if (getline(1:1) .ne. '#') goto 2
      goto 1

2     dataline = dataline + 1
      if (BTEST(debugflg,1)) then
        write (6, *) linecount, ':', dataline, ': ', getline
      endif

      return
      end
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++