!$Author: wjr $
!$Date: 2012-03-30 09:46:26 -0500 (Fri, 30 Mar 2012) $
!$Revision: 12162 $
!$HeadURL: https://svn.weru.ksu.edu/weru/weps1/trunk/weps.src/src/sweep/erodin.for $

      subroutine erodin (inp_unit,out_unit,cmdebugflag,already_read_inputs)

!     +++ PURPOSE +++
!     Utility to read initial conditions and variables from
!     input file (stdin or erod.in) for the standalone erosion submodel
!
!     If "out_unit" == stdout (6) then input not echo'd
    
      use constants_def
      use p1werm_def
      use p1erode_def
      use c1gen_def
      use m1sim_def
      use m1subr_def
      use m1geo_def
      use m1flag_def
      use anemometer_def
      use soilmod

!     +++ ARGUMENT DECLARATIONS +++
!
      integer inp_unit, out_unit, cmdebugflag, already_read_inputs
!
!     +++ ARGUMENT DEFINITIONS +++
!
!
!     +++ PARAMETERS +++
!

      integer mrcl
      parameter (mrcl = 512)
!
!     +++ ARGUMENT DECLARATIONS & +++
!     + + + GLOBAL COMMON BLOCKS + + +
      include 'b1glob.inc'
      include 'c1glob.inc'
      include 'd1glob.inc'
      include 's1sgeo.inc'
      include 'h1db1.inc'
      include 'w1wind.inc'
      include 'w1pavg.inc'
!
!     + + + LOCAL COMMON BLOCKS + + +
!
!      include 'erosion/p1erode.inc'
      integer debugflg
      common /flags/ debugflg
!
      integer xplot,xflag
      character*12 xcharin(30)
      real xin(30)
      common /plot/ xplot, xcharin, xin
!
!     +++ LOCAL VARIABLES +++
      integer i,j,k
      integer x,y,s,b,a,l,h
      integer wflg
      real f(maxTimeStepDay), wfcalm, wuc, w0k, step, wu(maxTimeStepDay)
!
!     + + + LOCAL VARIABLE DEFINITIONS + + +
!     i, j, k			= do-loop indices
!     x,y,s,b,a,l,h		= do-loop indices
!     wflg				= flag to determine format of wind speed data (0 = Weibull, 1 = real)
!     debugflg			= flag to output debug data (0 = none, 1 = input, 2 = more, etc.)
!     xplot				= flag to put plot data in arrays
!						  (value>0 = no. indep input variable, 0= none)
!     f(mntime)			= cumulative frequency of wind at speeds <avgWindSped(i)
!     wfcalm			= wind fraction intercept (+calm, - no calm in period)
!     wuc				= Weibull wind speed distribution scale factor (m/s)
!     w0k				= Weibull wind speed distribution shape factor
!     step				= tmp real variable for ntstep
!     xcharin(i)		= indep. variable name(s) used in plot
!     xin(i)			= indep. variable value(s) used in plot
!
!     +++ FUNCTIONS CALLED +++
!     getline
!
      character*(mrcl) getline
      character*(mrcl) line

!     +++ END SPECIFICATIONS +++
!
!      read  +++ MAIN  COMMON BLOCKS +++
!
!     +++ INIT STUFF +++
      if (already_read_inputs .gt. 0) goto 999  !Only write output 

      debugflg = 0 !needs to be initialized when using full debugging compiles
!     Read the debug flag to specify level of debug support
!     (currently only input file level debug supported)

      line = getline(inp_unit)
      if (cmdebugflag .lt. 0) then  !commandline option not set - use input file setting
          read (line,*) debugflg
!         read ((line=getline(inp_unit)),*) debugflg
       else
          debugflg = cmdebugflag  !use commandline setting
      endif

!     m1flag.inc

!     EROSION initialization flag (logical)
      line = getline(inp_unit)
      read (line,*) am0eif
!     read (getline(inp_unit),*) am0eif

!     EROSION "print" flag (integer)
      line = getline(inp_unit)
      read (line,*) am0efl
!     read (getline(inp_unit),*) am0efl

!     +++ SIMULATION REGION +++

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

!     m1geo.inc

!     Simulation region diagonal corners (x1,y1) and (x2,y2)
!     Note: form of amxsim(x,pt) where x=1, y=2, pt=1 or 2.
      line = getline(inp_unit)
      read (line,*) ((amxsim (x,y), x=1,2),y=1,2)
!     read (getline(inp_unit),*) ((amxsim (x,y), x=1,2),y=1,2)
!     write (out_unit,*) 'sim x1,y1 x1,y2 ',  amxsim(1,1), amxsim(1,2)
!     write (out_unit,*) 'sim x2,y1 x2,y2 ',  amxsim(2,1), amxsim(2,2)

!     Simulation region orientation angle
      line = getline(inp_unit)
      read (line,*) amasim
!     read (getline(inp_unit),*) amasim

!     +++ ACCOUNTING REGIONS +++

!     m1geo.inc

!     Number of accounting regions
      line = getline(inp_unit)
      read (line,*) nacctr
!     read (getline(inp_unit),*) nacctr

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

!     +++ BARRIERS +++

!     m1geo.inc

!     Number of barriers
      line = getline(inp_unit)
      read (line,*) nbr
!     read (getline(inp_unit),*) nbr

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

      do 30 b = 1, nbr
!       Barrier linear endpoints (x1,y1) and (x2,y2)
        line = getline(inp_unit)
        read (line,*) ((amxbr(x,y,b), x=1,2), y=1,2)
!       read (getline(inp_unit),*) ((amxbr(x,y,b), x=1,2), y=1,2)
!       Barrier ht, porosity, and width
        line = getline(inp_unit)
        read (line,*) amzbr(b), ampbr(b), amxbrw(b)
!       read (getline(inp_unit),*)  amzbr(b), ampbr(b), amxbrw(b)
   30 continue

!     +++ SUBREGION REGIONS +++

!     m1subr.inc

!     Number of subregions
      line = getline(inp_unit)
      read (line,*) numSubr
!     read (getline(inp_unit),*) numSubr

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

!     +++ BIOMASS +++

!     b1glob.inc & b1geom.inc

!       Biomass height
!        line = getline(inp_unit)
!        read (line,*)  abioMassHght(s)
!!       read (getline(inp_unit),*) abioMassHght(s)
!       Now reads in average "residue height" instead of "biomass height'
!       LEW - 1/26/06
        line = getline(inp_unit)
        read (line,*)  adzht_ave(s)

!     c1glob.inc

!       Crop height
        line = getline(inp_unit)
        read (line,*)  aczht(s)

!       Crop stem area index and leaf area index
        line = getline(inp_unit)
        read (line,*) acrsai(s), acrlai(s)

!     d1glob.inc

!       Residue stem area index and leaf area index
        line = getline(inp_unit)
        read (line,*) adrsaitot(s), adrlaitot(s)

!       use crop and residue values to find the total value
!       sum the stem area index and leaf area index values
        abrsai(s) = acrsai(s) + adrsaitot(s)
        abrlai(s) = acrlai(s) + adrlaitot(s)
      ! Compute the weighted average "biomass height" (residues and crop)
      ! which is used internally by the erosion code - LEW 1/26/06
        if (abrsai(s) .le. 0.0) then
            abioMassHght(s) = 0.0
        else
            abioMassHght(s) = (adzht_ave(s)*adrsaitot(s)+aczht(s)*acrsai(s))/  &
     &              abrsai(s)
        endif

!       addition to code for biodrag
!       crop row spacing and seed location
        line = getline(inp_unit)
        read (line,*) cropRowSpac(s), cropFurrFlg(s)
!
!       These aren't used in EROSION yet
!       Biomass stem area index by height
!       read (getline(inp_unit),*) (abrsaz(h,s), h=1,mncz)
!       Biomass leaf area index by height
!       read (getline(inp_unit),*) (abrlaz(h,s), h=1,mncz)

!       Biomass flat fraction cover, standing cover, and fraction total cover
!       read (getline(inp_unit),*) abioFracFlatCovr(s), abfscv(s), abftcv(s)
!       Only flat fraction cover used yet
        line = getline(inp_unit)
        read (line,*) abioFracFlatCovr(s)
!       read (getline(inp_unit),*) abioFracFlatCovr(s)

!     +++ SOIL +++

!       s1layr.inc, s1phys.inc & s1agg.inc

!       Number of soil layers (in this subregion)
        line = getline(inp_unit)
        read (line,*) SoilRegionData(s)%nslay
!       read (getline(inp_unit),*) SoilRegionData(s)%nslay

!       Soil layer thickness
        line = getline(inp_unit)
        read (line,*) (SoilRegionData(s)%SoilLayerData(l)%aszlyt,l=1,SoilRegionData(s)%nslay)
!       read (getline(inp_unit),*) (aszlyt(l,s),l=1,SoilRegionData(s)%nslay)

!       Soil layer bulk density
        line = getline(inp_unit)
        read (line,*) (SoilRegionData(s)%SoilLayerData(l)%asdblk, l=1,SoilRegionData(s)%nslay)
!       read (getline(inp_unit),*) (asdblk(l,s), l=1,SoilRegionData(s)%nslay)

!       Sand, silt, and clay fractions
        line = getline(inp_unit)
        read (line,*) (SoilRegionData(s)%SoilLayerData(l)%fracSand, l=1,SoilRegionData(s)%nslay)
!       read (getline(inp_unit),*) (SoilRegionData(s)%SoilLayerData(l)%fracSand, l=1,SoilRegionData(s)%nslay)

!       read very fine sand content edit 6-9-01 LH
        line = getline(inp_unit)
        read (line,*) (SoilRegionData(s)%SoilLayerData(l)%fracVeryFineSand, l=1,SoilRegionData(s)%nslay)

        line = getline(inp_unit)
        read (line,*) (SoilRegionData(s)%SoilLayerData(l)%fracSilt, l=1,SoilRegionData(s)%nslay)
!       read (getline(inp_unit),*) (SoilRegionData(s)%SoilLayerData(l)%fracSilt, l=1,SoilRegionData(s)%nslay)

        line = getline(inp_unit)
        read (line,*) (SoilRegionData(s)%SoilLayerData(l)%fracClay, l=1,SoilRegionData(s)%nslay)
!       read (getline(inp_unit),*) (SoilRegionData(s)%SoilLayerData(l)%fracClay, l=1,SoilRegionData(s)%nslay)

!       Volume of rock fraction
        line = getline(inp_unit)
        read (line,*) (SoilRegionData(s)%SoilLayerData(l)%fracRock, l=1,SoilRegionData(s)%nslay)

!       s1agg.inc
!       Soil layer aggregate density
        line = getline(inp_unit)
        read (line,*) (SoilRegionData(s)%SoilLayerData(l)%AggDen, l=1,SoilRegionData(s)%nslay)

!       Soil layer aggregate stability
        line = getline(inp_unit)
        read (line,*) (SoilRegionData(s)%SoilLayerData(l)%aggSta, l=1,SoilRegionData(s)%nslay)

! Check these variables with ASD inc files and Hagen's EROSION inc files - LEW
!       Soil layer ASD parms (gmd, min, max, gsd)
        line = getline(inp_unit)
        !read (line,*) (asoilLayrAggGMD(l,s), l=1,SoilRegionData(s)%nslay)
        read (line,*) (SoilRegionData(s)%SoilLayerData(l)%aggGMD, l=1,SoilRegionData(s)%nslay)

        line = getline(inp_unit)
        read (line,*) (SoilRegionData(s)%SoilLayerData(l)%aggMinSiz, l=1,SoilRegionData(s)%nslay)

        line = getline(inp_unit)
        read (line,*) (SoilRegionData(s)%SoilLayerData(l)%aggMaxSiz, l=1,SoilRegionData(s)%nslay)

        line = getline(inp_unit)
        read (line,*) (SoilRegionData(s)%SoilLayerData(l)%aggGSD, l=1,SoilRegionData(s)%nslay)

!       s1surf.inc & s1sgeo.inc

!       Crust parms (fraction, thickness)
        line = getline(inp_unit)
        read (line,*) SoilRegionData(s)%asoilCrstFrac, SoilRegionData(s)%asoilCrstThck,                               &
            SoilRegionData(s)%asoilLoosCovFrac, SoilRegionData(s)%asoilLoosMass,                                      &
            SoilRegionData(s)%asoilCrstDens, SoilRegionData(s)%asoilCrstStab
!       read (getline(inp_unit),*) SoilRegionData(s)%asoilCrstFrac, SoilRegionData(s)%asoilCrstThck,
!       Crust parms (fraction cover of loose material, mass loose material)
!       Crust parms (crust density and stability)

!       Random Roughness
        line = getline(inp_unit)
        read (line,*) SoilRegionData(s)%orgRandRoug
!       read (getline(inp_unit),*) SoilRegionData(s)%orgRandRoug

        !Lower and upper limits of grid cell RR allowed by erosion submodel
        if (SoilRegionData(s)%orgRandRoug < SLRR_MIN) then
           write(0,*) 'randRoug: ', SoilRegionData(s)%orgRandRoug,' < ', SLRR_MIN
        end if
        if (SoilRegionData(s)%orgRandRoug > SLRR_MAX) then
           write(0,*) 'randRoug: ', SoilRegionData(s)%orgRandRoug,' < ', SLRR_MIN
        end if

        !Lower and upper limits of grid cell aerodynamic roughness allowed
        !by erosion submodel (currently determined by equation used here)
        if (SoilRegionData(s)%orgRandRoug < (WZZO_MIN/0.3)) then
           write(0,*) 'randRoug: ', SoilRegionData(s)%orgRandRoug
           write(0,*) 'wzzo < WZZO_MIN: ', SoilRegionData(s)%orgRandRoug*0.3,' < ', WZZO_MIN
        else if(SoilRegionData(s)%orgRandRoug > (WZZO_MAX/0.3)) then
           write(0,*) 'randRoug: ', SoilRegionData(s)%orgRandRoug
           write(0,*) 'wzzo > WZZO_MAX: ', SoilRegionData(s)%orgRandRoug*0.3,' > ', WZZO_MAX
        end if

!       Oriented Roughness (ridge ht, spacing, width, orientation)
        line = getline(inp_unit)
        read (line,*) SoilRegionData(s)%orgRidgHght, SoilRegionData(s)%orgRidgSpac,                             &
!       read (getline(inp_unit),*) SoilRegionData(s)%orgRidgHght, SoilRegionData(s)%orgRidgSpac,
     &       SoilRegionData(s)%orgRidgWid, SoilRegionData(s)%orgRidgOrient

!       Oriented Roughness ( spacing)
        line = getline(inp_unit)
        read (line,*) SoilRegionData(s)%orgDikeSpac
!       read (getline(inp_unit),*) SoilRegionData(s)%orgDikeHght, SoilRegionData(s)%orgDikeSpac

!     +++ HYDROLOGY +++

!       h1db1.inc

!       Snow depth
        line = getline(inp_unit)
        read (line,*) asoilSurfSnowDept(s)
!       read (getline(inp_unit),*) asoilSurfSnowDept(s)

!       Soil layer wilting point
        line = getline(inp_unit)
        read (line,*) (SoilRegionData(s)%SoilLayerData(l)%asoilLayrWiltPt, l=1,SoilRegionData(s)%nslay)
!       read (getline(inp_unit),*) (asoilLayrWiltPt(l,s), l=1,SoilRegionData(s)%nslay)
!       Soil layer water content
        line = getline(inp_unit)
        read (line,*) (SoilRegionData(s)%SoilLayerData(l)%asoilSurfAvalWatrCont, l=1,SoilRegionData(s)%nslay)
!       read (getline(inp_unit),*) (asoilSurfAvalWatrCont(l,s), l=1,SoilRegionData(s)%nslay)
!       Soil surface hourly water content
        line = getline(inp_unit)
        read (line,*) (asoilSurfWatrCont(h,s), h=1,12)
!       read (getline(inp_unit),*) (asoilSurfWatrCont(h,s), h=1,12)
        line = getline(inp_unit)
        read (line,*) (asoilSurfWatrCont(h,s), h=13,24)
!       read (getline(inp_unit),*) (asoilSurfWatrCont(h,s), h=13,24)
  100 continue

!     +++ WEATHER +++

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

!     Wind Direction
      line = getline(inp_unit)
      read (line,*) windDir
!     read (getline(inp_unit),*) windDir

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

!     anemometer height, zo at anemom, and location (station or field)
!     note if flag=1, at field, awwzo will be changed to field value
      line = getline(inp_unit)
      read (line,*) anemHght, aeroRougStatHght, aeroAnemFlg

!     Weibull wind flag (0 - read Weibull parms, 1 - read wind speeds)
      line = getline(inp_unit)
      read (line,*) wflg
!     read (getline(inp_unit),*) wflg

!     w1wind.inc

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

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

!       calculate period wind speeds
        step = ntstep
        do 198 i= 1, ntstep
!         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)
!         to prevent out-of-range
          if (f(i) .lt. wfcalm) then
            f(i) = wfcalm
          endif
          avgWindSped(i) = wuc*(-log((1.0-f(i))/(1.0-wfcalm)))**(1.0/w0k)
  198   end do
!       Use greatest interval wind speed rather than 99% speed above
        dailyMaxWindSped = avgWindSped(ntstep)
!
!       change weibull wind speed dist. to a symmetric shape similar
!       to the daily distribution from wind gen
!
!       insure that ntstep is an even no.
        ntstep = (ntstep/2)*2
!
!       store wind speed in temp array
        do 110 i = 1, ntstep
          wu(i) = avgWindSped(i)
  110   end do
!
!       generate the symmetric distribution
        i = -1
        do 115 j = 1, ntstep/2
           i = i+2
           avgWindSped(j) = wu(i)
  115   continue
        i = ntstep+2
        do 125 j = (ntstep/2+1),ntstep
           i = i-2
           avgWindSped(j) = wu(i)
  125   continue
!
      else     ! when (wflg .eq. 1) input wind period data directly
        do 191 j = 1, ntstep/6
          line = getline(inp_unit)
          read (line,*) (avgWindSped(i),i=(j-1)*6+1,(j-1)*6+6)
!         read (getline(inp_unit),*) (avgWindSped(i),i=(j-1)*6+1,(j-1)*6+6)
191     end do
!       If not divisible evenly by 6, then get the remaining values
        if (mod(ntstep,6) .ne. 0) then
          line = getline(inp_unit)
          read (line,*)                                                 &
!         read (getline(inp_unit),*)
     &      (avgWindSped(i),i=(j-1)*6+1,(j-1)*6+mod(ntstep,6))
        endif

!     Determine the maximum wind speed during the day
        dailyMaxWindSped = 0.0
        do 193 i = 1, ntstep
           if (dailyMaxWindSped .lt. avgWindSped(i)) then
              dailyMaxWindSped = avgWindSped(i)
           endif
  193   end do

      endif

!     + + + PLOT SECTION + + +
!     selectively reads succesive indep variable names and values
!     if xplot is set to zero and the xflag is set to
!      1 for a given variable in the erodin file
!
!     test if plotout file to be created
      line = getline(inp_unit)
      read (line,*) xplot
      if (xplot .eq. 0) then
!     intialize xin array
      do 194 i=1,30
         xin(i) = 0.0
  194 continue
!
!     field length (good for wind parallel x-axis)
        line = getline(inp_unit)
        read (line,*) xflag
        if (xflag .eq. 1) then
          xplot = xplot + 1
          line = getline(inp_unit)
          xcharin(xplot) = line(1:30)
          xin(xplot) = amxsim(1,2) - amxsim(1,1)
        else
          line = getline(inp_unit)
        endif
! ^^^ tmp out
!      write (*,*) 'tmp out from erodin.for line 452'
!     write (*,*)  'xplot =', xplot
!      write (*,*)  'amxsim(2,1)=', amxsim(1,2)
!      write (*,*) 'xin(xplot)=', xin(xplot)
! ^^^ end tmp
!
!     biomass height
        line = getline(inp_unit)
        read (line,*) xflag
        if (xflag .eq. 1) then
          xplot = xplot + 1
          line = getline(inp_unit)
          xcharin(xplot) = line(1:30)
          xin(xplot) = abioMassHght(1)
        else
          line = getline(inp_unit)
        endif
! ^^^ tmp out
!     write (*,*) 'tmp out from erod.for line 474'
!     write (*,*)  'xplot =', xplot
!     write (*,*)  'abioMassHght(1)=', abioMassHght(1)
!     write (*,*) 'xin(xplot=', xin(xplot)
! ^^^ end tmp
!     biomass stem area index
        line = getline(inp_unit)
        read (line,*) xflag
        if (xflag .eq. 1) then
          xplot = xplot + 1
          line = getline(inp_unit)
          xcharin(xplot) = line(1:30)
          xin(xplot) = abrsai(1)
        else
          line = getline(inp_unit)
        endif
!     biomass leaf area index
        line = getline(inp_unit)
        read (line,*) xflag
        if (xflag .eq. 1) then
          xplot = xplot + 1
          line = getline(inp_unit)
          xcharin(xplot) = line(1:30)
          xin(xplot) = abrlai(1)
        else
          line = getline(inp_unit)
        endif
!     biomass flat cover
        line = getline(inp_unit)
        read (line,*) xflag
        if (xflag .eq. 1) then
          xplot = xplot + 1
          line = getline(inp_unit)
          xcharin(xplot) = line(1:30)
          xin(xplot) = abioFracFlatCovr(1)
        else
          line = getline(inp_unit)
        endif
!     very fine sand
        line = getline(inp_unit)
        read (line,*) xflag
        if (xflag .eq. 1) then
          xplot = xplot + 1
          line = getline(inp_unit)
          xcharin(xplot) = line(1:30)
          xin(xplot) = SoilRegionData(1)%SoilLayerData(1)%fracVeryFineSand
        else
          line = getline(inp_unit)
        endif
!     sand
        line = getline(inp_unit)
        read (line,*) xflag
        if (xflag .eq. 1) then
          xplot = xplot + 1
          line = getline(inp_unit)
          xcharin(xplot) = line(1:30)
          xin(xplot) =  SoilRegionData(1)%SoilLayerData(1)%fracSand
        else
          line = getline(inp_unit)
        endif
!     silt
        line = getline(inp_unit)
        read (line,*) xflag
        if (xflag .eq. 1) then
          xplot = xplot + 1
          line = getline(inp_unit)
          xcharin(xplot) = line(1:30)
          xin(xplot) = SoilRegionData(1)%SoilLayerData(1)%fracSilt
        else
          line = getline(inp_unit)
        endif
!     clay
        line = getline(inp_unit)
        read (line,*) xflag
        if (xflag .eq. 1) then
          xplot = xplot + 1
          line = getline(inp_unit)
          xcharin(xplot) = line(1:30)
          xin(xplot) =  SoilRegionData(1)%SoilLayerData(1)%fracClay
        else
          line = getline(inp_unit)
        endif
!     rock vol.
        line = getline(inp_unit)
        read (line,*) xflag
        if (xflag .eq. 1) then
          xplot = xplot + 1
          line = getline(inp_unit)
          xcharin(xplot) = line(1:30)
          xin(xplot) =  SoilRegionData(1)%SoilLayerData(1)%fracRock
        else
          line = getline(inp_unit)
        endif
!     aggregate density
        line = getline(inp_unit)
        read (line,*) xflag
        if (xflag .eq. 1) then
          xplot = xplot + 1
          line = getline(inp_unit)
          xcharin(xplot) = line(1:30)
          xin(xplot) = SoilRegionData(1)%SoilLayerData(1)%AggDen
        else
          line = getline(inp_unit)
        endif
!     aggregate stability
        line = getline(inp_unit)
        read (line,*) xflag
        if (xflag .eq. 1) then
          xplot = xplot + 1
          line = getline(inp_unit)
          xcharin(xplot) = line(1:30)
          xin(xplot) =  SoilRegionData(1)%SoilLayerData(1)%aggSta
        else
          line = getline(inp_unit)
        endif
!     agregate geometric mean diameter
        line = getline(inp_unit)
        read (line,*) xflag
        if (xflag .eq. 1) then
          xplot = xplot + 1
          line = getline(inp_unit)
          xcharin(xplot) = line(1:30)
          xin(xplot) =  SoilRegionData(1)%SoilLayerData(1)%aggGMD
        else
          line = getline(inp_unit)
        endif
!     aggreate minimum diameter
        line = getline(inp_unit)
        read (line,*) xflag
        if (xflag .eq. 1) then
          xplot = xplot + 1
          line = getline(inp_unit)
          xcharin(xplot) = line(1:30)
          xin(xplot) =  SoilRegionData(1)%SoilLayerData(1)%aggMinSiz
        else
          line = getline(inp_unit)
        endif
!     aggregate maximum diameter
        line = getline(inp_unit)
        read (line,*) xflag
        if (xflag .eq. 1) then
          xplot = xplot + 1
          line = getline(inp_unit)
          xcharin(xplot) = line(1:30)
          xin(xplot) =  SoilRegionData(1)%SoilLayerData(1)%aggMaxSiz
        else
          line = getline(inp_unit)
        endif
!     aggregate geometric std dev
        line = getline(inp_unit)
        read (line,*) xflag
        if (xflag .eq. 1) then
          xplot = xplot + 1
          line = getline(inp_unit)
          xcharin(xplot) = line(1:30)
          xin(xplot) =  SoilRegionData(1)%SoilLayerData(1)%aggGSD
        else
          line = getline(inp_unit)
        endif
!     soil fraction crust cover
        line = getline(inp_unit)
        read (line,*) xflag
        if (xflag .eq. 1) then
          xplot = xplot + 1
          line = getline(inp_unit)
          xcharin(xplot) = line(1:30)
          xin(xplot) = SoilRegionData(1)%asoilCrstFrac
        else
          line = getline(inp_unit)
        endif
!
!     surface crust thickness
        line = getline(inp_unit)
        read (line,*) xflag
        if (xflag .eq. 1) then
          xplot = xplot + 1
          line = getline(inp_unit)
          xcharin(xplot) = line(1:30)
          xin(xplot) = SoilRegionData(1)%asoilCrstThck
        else
          line = getline(inp_unit)
        endif
!     fraction loose material on crust
        line = getline(inp_unit)
        read (line,*) xflag
        if (xflag .eq. 1) then
          xplot = xplot + 1
          line = getline(inp_unit)
          xcharin(xplot) = line(1:30)
          xin(xplot) =  SoilRegionData(1)%asoilLoosCovFrac
        else
          line = getline(inp_unit)
        endif
!     mass of loose material on crust
        line = getline(inp_unit)
        read (line,*) xflag
        if (xflag .eq. 1) then
          xplot = xplot + 1
          line = getline(inp_unit)
          xcharin(xplot) = line(1:30)
          xin(xplot) =  SoilRegionData(1)%asoilLoosMass
        else
          line = getline(inp_unit)
        endif
!     soil crust stability
        line = getline(inp_unit)
        read (line,*) xflag
        if (xflag .eq. 1) then
          xplot = xplot + 1
          line = getline(inp_unit)
          xcharin(xplot) = line(1:30)
          xin(xplot) =  SoilRegionData(1)%asoilCrstStab
        else
          line = getline(inp_unit)
        endif
!     random roughness
        line = getline(inp_unit)
        read (line,*) xflag
        if (xflag .eq. 1) then
          xplot = xplot + 1
          line = getline(inp_unit)
          xcharin(xplot) = line(1:30)
          xin(xplot) =SoilRegionData(1)%orgRandRoug
        else
          line = getline(inp_unit)
        endif
!     ridge height
        line = getline(inp_unit)
        read (line,*) xflag
        if (xflag .eq. 1) then
          xplot = xplot + 1
          line = getline(inp_unit)
          xcharin(xplot) = line(1:30)
          xin(xplot) =  SoilRegionData(1)%orgRidgHght
        else
          line = getline(inp_unit)
        endif
!     ridge spacing
        line = getline(inp_unit)
        read (line,*) xflag
        if (xflag .eq. 1) then
          xplot = xplot + 1
          line = getline(inp_unit)
          xcharin(xplot) = line(1:30)
          xin(xplot) =  SoilRegionData(1)%orgRidgSpac
        else
          line = getline(inp_unit)
        endif
!     ridge width
        line = getline(inp_unit)
        read (line,*) xflag
        if (xflag .eq. 1) then
          xplot = xplot + 1
          line = getline(inp_unit)
          xcharin(xplot) = line(1:30)
          xin(xplot) =  SoilRegionData(1)%orgRidgWid
        else
          line = getline(inp_unit)
        endif
!     ridge orientation
        line = getline(inp_unit)
        read (line,*) xflag
        if (xflag .eq. 1) then
          xplot = xplot + 1
          line = getline(inp_unit)
          xcharin(xplot) = line(1:30)
          xin(xplot) =  SoilRegionData(1)%orgRidgOrient
        else
          line = getline(inp_unit)
        endif
      endif
!
! ^^^ tmp out
!     write(*,*)'^^^'
!     write(*,*)'tmp out from erodin.for'
!     write(*,*)'xcharin(1)=',xcharin(1)
!     write(*,*)'xin(1)=',xin(1)
!     write(*,*)'^^^'
! ^^^ end temp out

!     + + + OUTPUT SECTION + + +
!
  220 format (1x, 2f8.2)
  230 format (1x, 3f8.2)
  250 format (1x, 5f8.2)
  260 format (1x, 6f8.2)
  270 format (1x, 7f8.2)
  350 format (1x, i3, 7f8.2)
  400 format (1x, i3, 1x, l6, 4x, i3, 4x, i3, 4x, i3, 4x, i3)

999   if (out_unit .ne. 6) then  !Only echo input if stdout not specified
      write (out_unit,*)                                                  &
     &  '      REPORT OF INPUTS (read by erodin.for) '
      write (out_unit,*)
      write (out_unit,*)  '+++ Control Flags, etc. +++'
      write (out_unit,*)
      write (out_unit,*) 'ntstep  am0eif  numSubr  nacctr  nbr am0efl'
      write (out_unit,400) ntstep, am0eif, numSubr, nacctr, nbr,am0efl

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

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

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

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

      do 1100 s = 1, numSubr
        write (out_unit,*)
        write (out_unit,*)
        write(out_unit,*) '*********************** Subregion ', s,        &
     &                 ' ***********************'

        write (out_unit,*)
        write(out_unit,*) '+++ BIOMASS +++'
        write (out_unit,*)
        write (out_unit,*) 'Biomass ht,  flat cover'
        write (out_unit,220) abioMassHght(s), abioFracFlatCovr(s)

        write (out_unit,*)
        write (out_unit,*) 'Crop height, SAI,    LAI'
        write (out_unit,230) aczht(s), acrsai(s), acrlai(s)

        write (out_unit,*)
        write (out_unit,*) 'Residue height, SAI,    LAI'
        write (out_unit,230) adzht_ave(s), adrsaitot(s), adrlaitot(s)

        write (out_unit,*)
        write (out_unit,*) '+++ SOIL +++ '
        write (out_unit,*)
        write (out_unit,*) 'nslay - number of soil layers'
        write (out_unit,*) SoilRegionData(s)%nslay

        write (out_unit,*)
        write (out_unit,*) 'layer depth b.density ',                      &
     &                   'vfsand   sand   silt   clay    rock vol'
        do 1030 l = 1, SoilRegionData(s)%nslay
          write (out_unit,350) l, SoilRegionData(s)%SoilLayerData(l)%aszlyt, &
                             SoilRegionData(s)%SoilLayerData(l)%asdblk,    &
                             SoilRegionData(s)%SoilLayerData(l)%fracVeryFineSand,    &
                             SoilRegionData(s)%SoilLayerData(l)%fracSand,  &
                             SoilRegionData(s)%SoilLayerData(l)%fracSilt,  &
                             SoilRegionData(s)%SoilLayerData(l)%fracClay,  &
                             SoilRegionData(s)%SoilLayerData(l)%fracRock
 1030   continue

        write (out_unit,*)
          write (out_unit,*) 'layer    AgD     AgS ',                     &
     &                   ' GMD    GMDmn     GMDmx    GSD'
        do 1040 l = 1, SoilRegionData(s)%nslay
          write (out_unit,350) l, SoilRegionData(s)%SoilLayerData(l)%AggDen,  &
                             SoilRegionData(s)%SoilLayerData(l)%aggSta,     &
                             SoilRegionData(s)%SoilLayerData(l)%aggGMD,     &
                             SoilRegionData(s)%SoilLayerData(l)%aggMinSiz,  &
                             SoilRegionData(s)%SoilLayerData(l)%aggMaxSiz,  &
                             SoilRegionData(s)%SoilLayerData(l)%aggGSD
 1040   continue

        write (out_unit,*)
        write (out_unit,*) 'Crust frac thick mass LOS frac.LOS, ',        &
     &                   'density stability'
        write (out_unit,260)SoilRegionData(s)%asoilCrstFrac, SoilRegionData(s)%asoilCrstThck, &
                          SoilRegionData(s)%asoilLoosMass, SoilRegionData(s)%asoilLoosCovFrac, &
                          SoilRegionData(s)%asoilCrstDens, SoilRegionData(s)%asoilCrstStab

        write (out_unit,*)
        write (out_unit,*) '    RR,    Rg ht,  width, spacing, ',         &
     &                   'orient., dike spacing'
        write (out_unit,270) SoilRegionData(s)%orgRandRoug, SoilRegionData(s)%orgRidgHght, &
                           SoilRegionData(s)%orgRidgWid, SoilRegionData(s)%orgRidgSpac,  &
                           SoilRegionData(s)%orgRidgOrient,  SoilRegionData(s)%orgDikeSpac

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

        write (out_unit,*) 'Snow depth (mm)'
        write (out_unit,*) asoilSurfSnowDept(s)

        write (out_unit,*)
        write (out_unit,*) 'layer  wilting and actual water contents'
        do 1050 l = 1, SoilRegionData(s)%nslay
          write (out_unit,350) l, SoilRegionData(s)%SoilLayerData(l)%asoilLayrWiltPt, &
                             SoilRegionData(s)%SoilLayerData(l)%asoilSurfAvalWatrCont
 1050   continue
        write (out_unit,*) 'Hourly water contents - asoilSurfWatrCont'
        write (out_unit,260) (asoilSurfWatrCont(h,s), h=1,6)
        write (out_unit,260) (asoilSurfWatrCont(h,s), h=7,12)
        write (out_unit,260) (asoilSurfWatrCont(h,s), h=13,18)
        write (out_unit,260) (asoilSurfWatrCont(h,s), h=19,24)

 1100 continue

      write (out_unit,*)
      write (out_unit,*) '+++ WEATHER +++'
      write (out_unit,*)
      write (out_unit,*) ' anemHght    awwzo   aeroAnemFlg '
      write (out_unit,*)   anemHght, aeroRougStatHght, aeroAnemFlg
      write (out_unit,*) ' wind dir (deg) and max wind speed (m/s)'
      write (out_unit,220) windDir, dailyMaxWindSped

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

      endif !(out_unit .ne. 6)

      !close (inp_unit)
      return
      end

!**********************************************************************
      character*(*) function getline(inp_unit)

      integer debugflg
      common /flags/ debugflg

      integer inp_unit
      integer dataline, linecount

      save dataline, linecount


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

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

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

      return
      end
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++