!*==erodin.spg  processed by SPAG 6.70Rc at 15:33 on 10 Dec 2012
!*------------------   SPAG Configuration Options   --------------------
!*--0323,76 000101,-1 000000102011332010100002000000210211210,136 10 --
!*--1100000011112111000000000000,10,10,10,10,10,10,900,100 200000000 --
!*--000000010000000000000,72,72 73,42,38,33 00011112110000100000000  --
!*----------------------------------------------------------------------
!$Author: joelevin $
!$Date: 2012-03-30 10:14:38 -0500 (Fri, 30 Mar 2012) $
!$Revision: 12163 $
!$HeadURL: https://eweru-dev1.eweru.ksu.edu/svn/code/weps1/branches/WEPS_F90_update/weps.src/src/sweep/erodin.for $
 
      subroutine erodin(i_unit,o_unit,cmdebugflag,already_read_inputs)
      use i_p1werm
      use i_c1gen
      use i_p1const
      use i_m1sim
      use i_m1flag
      use i_m1subr
      use i_m1geo
      use i_b1glob
      use i_c1glob
      use i_d1glob
      use i_s1layr
      use i_s1phys
      use i_s1dbh
      use i_s1agg
      use i_s1surf
      use i_s1sgeo
      use i_h1db1
      use i_w1wind
      use i_w1pavg
      use i_p1erode
      use c_flags
      use c_plot
      use s_getline
      implicit none
!*--ERODIN37
!
!*** Start of declarations rewritten by SPAG
!
! PARAMETER definitions
!
      integer,parameter :: mrcl = 512
!
! Dummy arguments
!
      integer :: already_read_inputs,cmdebugflag,i_unit,o_unit
      intent (in) already_read_inputs,cmdebugflag,o_unit
!
! Local variables
!
      integer :: a,b,h,i,j,k,l,s,wflg,x,xflag,y
      real,dimension(mntime) :: f,wu
      character(mrcl) :: line
      real :: step,w0k,wfcalm,wuc
!
!*** End of declarations rewritten by SPAG
!
 
!     +++ PURPOSE +++
!     Utility to read initial conditions and variables from
!     input file (stdin or erod.in) for the standalone erosion submodel
!
!     If "o_unit" == stdout (6) then input not echo'd
!
!     +++ ARGUMENT DECLARATIONS +++
!
!
!     +++ ARGUMENT DEFINITIONS +++
!
!
!     +++ PARAMETERS +++
!
!
!     +++ ARGUMENT DECLARATIONS & +++
!     + + + GLOBAL COMMON BLOCKS + + +
!
!     + + + LOCAL COMMON BLOCKS + + +
!
!
!
!     +++ LOCAL VARIABLES +++
!
!     + + + 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 <awu(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
!
 
!     +++ END SPECIFICATIONS +++
!
!      read  +++ MAIN  COMMON BLOCKS +++
!
!     +++ INIT STUFF +++
      if (already_read_inputs<=0) then          !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(i_unit)
         if (cmdebugflag<0) then    !commandline option not set - use input file setting
            read (line,*) debugflg
!         read ((line=getline(i_unit)),*) debugflg
         else
            debugflg = cmdebugflag
                                  !use commandline setting
         end if
 
!     m1flag.inc
 
!     EROSION initialization flag (logical)
         line = getline(i_unit)
         read (line,*) am0eif
!     read (getline(i_unit),*) am0eif
 
!     EROSION "print" flag (integer)
         line = getline(i_unit)
         read (line,*) am0efl
!     read (getline(i_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(i_unit)
         read (line,*) ((amxsim(x,y),x=1,2),y=1,2)
!     read (getline(i_unit),*) ((amxsim (x,y), x=1,2),y=1,2)
!     write (o_unit,*) 'sim x1,y1 x1,y2 ',  amxsim(1,1), amxsim(1,2)
!     write (o_unit,*) 'sim x2,y1 x2,y2 ',  amxsim(2,1), amxsim(2,2)
 
!     Simulation region orientation angle
         line = getline(i_unit)
         read (line,*) amasim
!     read (getline(i_unit),*) amasim
 
!     +++ ACCOUNTING REGIONS +++
 
!     m1geo.inc
 
!     Number of accounting regions
         line = getline(i_unit)
         read (line,*) nacctr
!     read (getline(i_unit),*) nacctr
 
!     Accounting Region diagonal corners (x1,y1) and (x2,y2)
         do a = 1,nacctr
            line = getline(i_unit)
            read (line,*) ((amxar(x,y,a),x=1,2),y=1,2)
!       read (getline(i_unit),*) ((amxar (x,y,a), x=1,2), y=1,2)
         end do
 
!     +++ BARRIERS +++
 
!     m1geo.inc
 
!     Number of barriers
         line = getline(i_unit)
         read (line,*) nbr
!     read (getline(i_unit),*) nbr
 
!     NOTE: Barrier data must not be in the input file if "nbr = 0"
 
         do b = 1,nbr
!       Barrier linear endpoints (x1,y1) and (x2,y2)
            line = getline(i_unit)
            read (line,*) ((amxbr(x,y,b),x=1,2),y=1,2)
!       read (getline(i_unit),*) ((amxbr(x,y,b), x=1,2), y=1,2)
!       Barrier ht, porosity, and width
            line = getline(i_unit)
            read (line,*) amzbr(b),ampbr(b),amxbrw(b)
!       read (getline(i_unit),*)  amzbr(b), ampbr(b), amxbrw(b)
         end do
 
!     +++ SUBREGION REGIONS +++
 
!     m1subr.inc
 
!     Number of subregions
         line = getline(i_unit)
         read (line,*) nsubr
!     read (getline(i_unit),*) nsubr
 
!     Dimensions, Biomass, Soil, and Hydrology (by subregion)
         do s = 1,nsubr
!       Subregion diagonal corners (x1,y1) and (x2,y2)
            line = getline(i_unit)
            read (line,*) ((amxsr(x,y,s),x=1,2),y=1,2)
!       read (getline(i_unit),*) ((amxsr (x,y,s), x=1,2),y=1,2)
 
!     +++ BIOMASS +++
 
!     b1glob.inc & b1geom.inc
 
!       Biomass height
!        line = getline(i_unit)
!        read (line,*)  abzht(s)
!!       read (getline(i_unit),*) abzht(s)
!       Now reads in average "residue height" instead of "biomass height'
!       LEW - 1/26/06
            line = getline(i_unit)
            read (line,*) adzht_ave(s)
 
!     c1glob.inc
 
!       Crop height
            line = getline(i_unit)
            read (line,*) aczht(s)
 
!       Crop stem area index and leaf area index
            line = getline(i_unit)
            read (line,*) acrsai(s),acrlai(s)
 
!     d1glob.inc
 
!       Residue stem area index and leaf area index
            line = getline(i_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)<=0.0) then
               abzht(s) = 0.0
            else
               abzht(s) = (adzht_ave(s)*adrsaitot(s)+aczht(s)*acrsai(s))&
                        & /&abrsai(s)
            end if
!     c1gen.inc
 
!       addition to code for biodrag
!       crop row spacing and seed location
            line = getline(i_unit)
            read (line,*) acxrow(s),ac0rg(s)
!
!       These aren't used in EROSION yet
!       Biomass stem area index by height
!       read (getline(i_unit),*) (abrsaz(h,s), h=1,mncz)
!       Biomass leaf area index by height
!       read (getline(i_unit),*) (abrlaz(h,s), h=1,mncz)
 
!       Biomass flat fraction cover, standing cover, and fraction total cover
!       read (getline(i_unit),*) abffcv(s), abfscv(s), abftcv(s)
!       Only flat fraction cover used yet
            line = getline(i_unit)
            read (line,*) abffcv(s)
!       read (getline(i_unit),*) abffcv(s)
 
!     +++ SOIL +++
 
!       s1layr.inc, s1phys.inc & s1agg.inc
 
!       Number of soil layers (in this subregion)
            line = getline(i_unit)
            read (line,*) nslay(s)
!       read (getline(i_unit),*) nslay(s)
 
!       Soil layer thickness
            line = getline(i_unit)
            read (line,*) (aszlyt(l,s),l=1,nslay(s))
!       read (getline(i_unit),*) (aszlyt(l,s),l=1,nslay(s))
 
!       Soil layer bulk density
            line = getline(i_unit)
            read (line,*) (asdblk(l,s),l=1,nslay(s))
!       read (getline(i_unit),*) (asdblk(l,s), l=1,nslay(s))
 
!       Sand, silt, and clay fractions
            line = getline(i_unit)
            read (line,*) (asfsan(l,s),l=1,nslay(s))
!       read (getline(i_unit),*) (asfsan(l,s), l=1,nslay(s))
 
!       read very fine sand content edit 6-9-01 LH
            line = getline(i_unit)
            read (line,*) (asfvfs(l,s),l=1,nslay(s))
 
            line = getline(i_unit)
            read (line,*) (asfsil(l,s),l=1,nslay(s))
!       read (getline(i_unit),*) (asfsil(l,s), l=1,nslay(s))
 
            line = getline(i_unit)
            read (line,*) (asfcla(l,s),l=1,nslay(s))
!       read (getline(i_unit),*) (asfcla(l,s), l=1,nslay(s))
 
!       Volume of rock fraction
            line = getline(i_unit)
            read (line,*) (asvroc(l,s),l=1,nslay(s))
 
!       s1agg.inc
!       Soil layer aggregate density
            line = getline(i_unit)
            read (line,*) (asdagd(l,s),l=1,nslay(s))
 
!       Soil layer aggregate stability
            line = getline(i_unit)
            read (line,*) (aseags(l,s),l=1,nslay(s))
 
! Check these variables with ASD inc files and Hagen's EROSION inc files - LEW
!       Soil layer ASD parms (gmd, min, max, gsd)
            line = getline(i_unit)
            read (line,*) (aslagm(l,s),l=1,nslay(s))
 
            line = getline(i_unit)
            read (line,*) (aslagn(l,s),l=1,nslay(s))
 
            line = getline(i_unit)
            read (line,*) (aslagx(l,s),l=1,nslay(s))
 
            line = getline(i_unit)
            read (line,*) (as0ags(l,s),l=1,nslay(s))
 
!       s1surf.inc & s1sgeo.inc
 
!       Crust parms (fraction, thickness)
            line = getline(i_unit)
!       read (getline(i_unit),*) asfcr(s), aszcr(s),
!       Crust parms (fraction cover of loose material, mass loose material)
!       Crust parms (crust density and stability)
            read (line,*) asfcr(s),aszcr(s),&asflos(s),asmlos(s),       &
                        & &asdcr(s),asecr(s)
 
!       Random Roughness
            line = getline(i_unit)
            read (line,*) aslrr(s)
!       read (getline(i_unit),*) aslrr(s)
 
        !Lower and upper limits of grid cell RR allowed by erosion submodel
            if (aslrr(s)<slrr_min) write (0,*) 'slrr: ',aslrr(s),' < ', &
              & slrr_min
            if (aslrr(s)>slrr_max) write (0,*) 'slrr: ',aslrr(s),' < ', &
              & slrr_min
 
        !Lower and upper limits of grid cell aerodynamic roughness allowed
        !by erosion submodel (currently determined by equation used here)
            if (aslrr(s)<(wzzo_min/0.3)) then
               write (0,*) 'slrr: ',aslrr(s)
               write (0,*) 'wzzo < WZZO_MIN: ',aslrr(s)*0.3,' < ',      &
                         & wzzo_min
            else if (aslrr(s)>(wzzo_max/0.3)) then
               write (0,*) 'slrr: ',aslrr(s)
               write (0,*) 'wzzo > WZZO_MAX: ',aslrr(s)*0.3,' > ',      &
                         & wzzo_max
            end if
 
!       Oriented Roughness (ridge ht, spacing, width, orientation)
            line = getline(i_unit)
!       read (getline(i_unit),*) aszrgh(s), asxrgs(s),
            read (line,*) aszrgh(s),asxrgs(s),&asxrgw(s),asargo(s)
 
!       Oriented Roughness ( spacing)
            line = getline(i_unit)
            read (line,*) asxdks(s)
!       read (getline(i_unit),*) asxdkh(s), asxdks(s)
 
!     +++ HYDROLOGY +++
 
!       h1db1.inc
 
!       Snow depth
            line = getline(i_unit)
            read (line,*) ahzsnd(s)
!       read (getline(i_unit),*) ahzsnd(s)
 
!       Soil layer wilting point
            line = getline(i_unit)
            read (line,*) (ahrwcw(l,s),l=1,nslay(s))
!       read (getline(i_unit),*) (ahrwcw(l,s), l=1,nslay(s))
!       Soil layer water content
            line = getline(i_unit)
            read (line,*) (ahrwca(l,s),l=1,nslay(s))
!       read (getline(i_unit),*) (ahrwca(l,s), l=1,nslay(s))
!       Soil surface hourly water content
            line = getline(i_unit)
            read (line,*) (ahrwc0(h,s),h=1,12)
!       read (getline(i_unit),*) (ahrwc0(h,s), h=1,12)
            line = getline(i_unit)
            read (line,*) (ahrwc0(h,s),h=13,24)
!       read (getline(i_unit),*) (ahrwc0(h,s), h=13,24)
         end do
 
!     +++ 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(i_unit)
         read (line,*) awdair
!     read (getline(i_unit),*) awdair
 
!     Wind Direction
         line = getline(i_unit)
         read (line,*) awadir
!     read (getline(i_unit),*) awadir
 
!     Number of "steps" during 24 hours (96 = 15 minute intervals)
         line = getline(i_unit)
         read (line,*) ntstep
!     read (getline(i_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(i_unit)
         read (line,*) anemht,awzzo,wzoflg
 
!     Weibull wind flag (0 - read Weibull parms, 1 - read wind speeds)
         line = getline(i_unit)
         read (line,*) wflg
!     read (getline(i_unit),*) wflg
 
!     w1wind.inc
 
!     wind data inputs as the Weibull paramters
!     (wfcalm, wuc, w0k) is indicated by code ntstep = 99
         if (wflg==0) then
!       Weibull parms (fraction calm, c, k)
            line = getline(i_unit)
            read (line,*) wfcalm,wuc,w0k
!       read (getline(i_unit),*) wfcalm, wuc, w0k
 
!       calculate daily max wind speed (99% speed)
!       awudmx = wuc*(-log((1.0-0.99)/(1-wfcalm)))**(1.0/w0k)
 
!       calculate period wind speeds
            step = ntstep
            do 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)<wfcalm) f(i) = wfcalm
               awu(i) = wuc*(-log((1.0-f(i))/(1.0-wfcalm)))**(1.0/w0k)
            end do
!       Use greatest interval wind speed rather than 99% speed above
            awudmx = awu(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 i = 1,ntstep
               wu(i) = awu(i)
            end do
!
!       generate the symmetric distribution
            i = -1
            do j = 1,ntstep/2
               i = i + 2
               awu(j) = wu(i)
            end do
            i = ntstep + 2
            do j = (ntstep/2+1),ntstep
               i = i - 2
               awu(j) = wu(i)
            end do
!
         else  ! when (wflg .eq. 1) input wind period data directly
            do j = 1,ntstep/6
               line = getline(i_unit)
               read (line,*) (awu(i),i=(j-1)*6+1,(j-1)*6+6)
!         read (getline(i_unit),*) (awu(i),i=(j-1)*6+1,(j-1)*6+6)
            end do
!       If not divisible evenly by 6, then get the remaining values
            if (mod(ntstep,6)/=0) then
               line = getline(i_unit)
!         read (getline(i_unit),*)
               read (line,*) &(awu(i),i=(j-1)*6+1,(j-1)*6+mod(ntstep,6))
            end if
 
!     Determine the maximum wind speed during the day
            awudmx = 0.0
            do i = 1,ntstep
               if (awudmx<awu(i)) awudmx = awu(i)
            end do
 
         end if
 
!     + + + 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(i_unit)
         read (line,*) xplot
         if (xplot==0) then
!     intialize xin array
            do i = 1,30
               xin(i) = 0.0
            end do
!
!     field length (good for wind parallel x-axis)
            line = getline(i_unit)
            read (line,*) xflag
            if (xflag==1) then
               xplot = xplot + 1
               line = getline(i_unit)
               xcharin(xplot) = line
               xin(xplot) = amxsim(1,2) - amxsim(1,1)
            else
               line = getline(i_unit)
            end if
! ^^^ 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(i_unit)
            read (line,*) xflag
            if (xflag==1) then
               xplot = xplot + 1
               line = getline(i_unit)
               xcharin(xplot) = line
               xin(xplot) = abzht(1)
            else
               line = getline(i_unit)
            end if
! ^^^ tmp out
!     write (*,*) 'tmp out from erod.for line 474'
!     write (*,*)  'xplot =', xplot
!     write (*,*)  'abzht(1)=', abzht(1)
!     write (*,*) 'xin(xplot=', xin(xplot)
! ^^^ end tmp
!     biomass stem area index
            line = getline(i_unit)
            read (line,*) xflag
            if (xflag==1) then
               xplot = xplot + 1
               line = getline(i_unit)
               xcharin(xplot) = line
               xin(xplot) = abrsai(1)
            else
               line = getline(i_unit)
            end if
!     biomass leaf area index
            line = getline(i_unit)
            read (line,*) xflag
            if (xflag==1) then
               xplot = xplot + 1
               line = getline(i_unit)
               xcharin(xplot) = line
               xin(xplot) = abrlai(1)
            else
               line = getline(i_unit)
            end if
!     biomass flat cover
            line = getline(i_unit)
            read (line,*) xflag
            if (xflag==1) then
               xplot = xplot + 1
               line = getline(i_unit)
               xcharin(xplot) = line
               xin(xplot) = abffcv(1)
            else
               line = getline(i_unit)
            end if
!     very fine sand
            line = getline(i_unit)
            read (line,*) xflag
            if (xflag==1) then
               xplot = xplot + 1
               line = getline(i_unit)
               xcharin(xplot) = line
               xin(xplot) = asfvfs(1,1)
            else
               line = getline(i_unit)
            end if
!     sand
            line = getline(i_unit)
            read (line,*) xflag
            if (xflag==1) then
               xplot = xplot + 1
               line = getline(i_unit)
               xcharin(xplot) = line
               xin(xplot) = asfsan(1,1)
            else
               line = getline(i_unit)
            end if
!     silt
            line = getline(i_unit)
            read (line,*) xflag
            if (xflag==1) then
               xplot = xplot + 1
               line = getline(i_unit)
               xcharin(xplot) = line
               xin(xplot) = asfsil(1,1)
            else
               line = getline(i_unit)
            end if
!     clay
            line = getline(i_unit)
            read (line,*) xflag
            if (xflag==1) then
               xplot = xplot + 1
               line = getline(i_unit)
               xcharin(xplot) = line
               xin(xplot) = asfcla(1,1)
            else
               line = getline(i_unit)
            end if
!     rock vol.
            line = getline(i_unit)
            read (line,*) xflag
            if (xflag==1) then
               xplot = xplot + 1
               line = getline(i_unit)
               xcharin(xplot) = line
               xin(xplot) = asvroc(1,1)
            else
               line = getline(i_unit)
            end if
!     aggregate density
            line = getline(i_unit)
            read (line,*) xflag
            if (xflag==1) then
               xplot = xplot + 1
               line = getline(i_unit)
               xcharin(xplot) = line
               xin(xplot) = asdagd(1,1)
            else
               line = getline(i_unit)
            end if
!     aggregate stability
            line = getline(i_unit)
            read (line,*) xflag
            if (xflag==1) then
               xplot = xplot + 1
               line = getline(i_unit)
               xcharin(xplot) = line
               xin(xplot) = aseags(1,1)
            else
               line = getline(i_unit)
            end if
!     agregate geometric mean diameter
            line = getline(i_unit)
            read (line,*) xflag
            if (xflag==1) then
               xplot = xplot + 1
               line = getline(i_unit)
               xcharin(xplot) = line
               xin(xplot) = aslagm(1,1)
            else
               line = getline(i_unit)
            end if
!     aggreate minimum diameter
            line = getline(i_unit)
            read (line,*) xflag
            if (xflag==1) then
               xplot = xplot + 1
               line = getline(i_unit)
               xcharin(xplot) = line
               xin(xplot) = aslagn(1,1)
            else
               line = getline(i_unit)
            end if
!     aggregate maximum diameter
            line = getline(i_unit)
            read (line,*) xflag
            if (xflag==1) then
               xplot = xplot + 1
               line = getline(i_unit)
               xcharin(xplot) = line
               xin(xplot) = aslagx(1,1)
            else
               line = getline(i_unit)
            end if
!     aggregate geometric std dev
            line = getline(i_unit)
            read (line,*) xflag
            if (xflag==1) then
               xplot = xplot + 1
               line = getline(i_unit)
               xcharin(xplot) = line
               xin(xplot) = as0ags(1,1)
            else
               line = getline(i_unit)
            end if
!     soil fraction crust cover
            line = getline(i_unit)
            read (line,*) xflag
            if (xflag==1) then
               xplot = xplot + 1
               line = getline(i_unit)
               xcharin(xplot) = line
               xin(xplot) = asfcr(1)
            else
               line = getline(i_unit)
            end if
!
!     surface crust thickness
            line = getline(i_unit)
            read (line,*) xflag
            if (xflag==1) then
               xplot = xplot + 1
               line = getline(i_unit)
               xcharin(xplot) = line
               xin(xplot) = aszcr(1)
            else
               line = getline(i_unit)
            end if
!     fraction loose material on crust
            line = getline(i_unit)
            read (line,*) xflag
            if (xflag==1) then
               xplot = xplot + 1
               line = getline(i_unit)
               xcharin(xplot) = line
               xin(xplot) = asflos(1)
            else
               line = getline(i_unit)
            end if
!     mass of loose material on crust
            line = getline(i_unit)
            read (line,*) xflag
            if (xflag==1) then
               xplot = xplot + 1
               line = getline(i_unit)
               xcharin(xplot) = line
               xin(xplot) = asmlos(1)
            else
               line = getline(i_unit)
            end if
!     soil crust stability
            line = getline(i_unit)
            read (line,*) xflag
            if (xflag==1) then
               xplot = xplot + 1
               line = getline(i_unit)
               xcharin(xplot) = line
               xin(xplot) = asecr(1)
            else
               line = getline(i_unit)
            end if
!     random roughness
            line = getline(i_unit)
            read (line,*) xflag
            if (xflag==1) then
               xplot = xplot + 1
               line = getline(i_unit)
               xcharin(xplot) = line
               xin(xplot) = aslrr(1)
            else
               line = getline(i_unit)
            end if
!     ridge height
            line = getline(i_unit)
            read (line,*) xflag
            if (xflag==1) then
               xplot = xplot + 1
               line = getline(i_unit)
               xcharin(xplot) = line
               xin(xplot) = aszrgh(1)
            else
               line = getline(i_unit)
            end if
!     ridge spacing
            line = getline(i_unit)
            read (line,*) xflag
            if (xflag==1) then
               xplot = xplot + 1
               line = getline(i_unit)
               xcharin(xplot) = line
               xin(xplot) = asxrgs(1)
            else
               line = getline(i_unit)
            end if
!     ridge width
            line = getline(i_unit)
            read (line,*) xflag
            if (xflag==1) then
               xplot = xplot + 1
               line = getline(i_unit)
               xcharin(xplot) = line
               xin(xplot) = asxrgw(1)
            else
               line = getline(i_unit)
            end if
!     ridge orientation
            line = getline(i_unit)
            read (line,*) xflag
            if (xflag==1) then
               xplot = xplot + 1
               line = getline(i_unit)
               xcharin(xplot) = line
               xin(xplot) = asargo(1)
            else
               line = getline(i_unit)
            end if
         end if
!
! ^^^ 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 + + +
!
 1000    format (1x,f6.2,6F8.2)
 1100    format (1x,f8.2)
 1200    format (1x,2F8.2)
 1300    format (1x,3F8.2)
 1400    format (1x,4F8.3)
 1500    format (1x,5F8.2)
 1600    format (1x,6F8.2)
 1700    format (1x,7F8.2)
 1800    format (1x,8F8.2)
 1900    format (1x,8I5)
 2000    format (1x,i3,7F8.2)
 2100    format (1x,i3,1x,l6,4x,i3,4x,i3,4x,i3,4x,i3)
      end if
 
      if (o_unit/=6) then      !Only echo input if stdout not specified
         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 am0efl'
         write (o_unit,2100) ntstep,am0eif,nsubr,nacctr,nbr,am0efl
 
         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,1500) 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 a = 1,nacctr
            write (o_unit,1500) ((amxar(x,y,a),x=1,2),y=1,2)
         end do
 
         write (o_unit,*)
         write (o_unit,*) '+++ BARRIERS +++'
         write (o_unit,*)
         write (o_unit,*) 'nbr - number of barriers'
         write (o_unit,*) nbr
         if (nbr>0) then
            write (o_unit,*) 'barrier dim (x1,y1) (x2,y2) ',            &
                            &&'and height, porosity, and width'
            do b = 1,nbr
               write (o_unit,1700) ((amxbr(x,y,b),x=1,2),y=1,2),        &
                                 & &amzbr(b),ampbr(b),amxbrw(b)
            end do
         end if
 
         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 s = 1,nsubr
            write (o_unit,1500) ((amxsr(i,j,s),i=1,2),j=1,2)
         end do
 
         do 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 ht,  flat cover'
            write (o_unit,1200) abzht(s),abffcv(s)
 
            write (o_unit,*)
            write (o_unit,*) 'Crop height, SAI,    LAI'
            write (o_unit,1300) aczht(s),acrsai(s),acrlai(s)
 
            write (o_unit,*)
            write (o_unit,*) 'Residue height, SAI,    LAI'
            write (o_unit,1300) adzht_ave(s),adrsaitot(s),adrlaitot(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 depth b.density ',                  &
                            &&'vfsand   sand   silt   clay    rock vol'
            do l = 1,nslay(s)
               write (o_unit,2000) l,aszlyt(l,s),asdblk(l,s),           &
                                 & &asfvfs(l,s),asfsan(l,s),asfsil(l,s),&
                                 & asfcla(l,s),asvroc(l,s)
            end do
 
            write (o_unit,*)
            write (o_unit,*) 'layer    AgD     AgS ',                   &
                            &&' GMD    GMDmn     GMDmx    GSD'
            do l = 1,nslay(s)
               write (o_unit,2000) l,asdagd(l,s),aseags(l,s),           &
                                 & &aslagm(l,s),aslagn(l,s),aslagx(l,s),&
                                 & as0ags(l,s)
            end do
 
            write (o_unit,*)
            write (o_unit,*) 'Crust frac thick mass LOS frac.LOS, ',    &
                            &&'density stability'
            write (o_unit,1600) asfcr(s),aszcr(s),asmlos(s),asflos(s),  &
                              & &asdcr(s),asecr(s)
 
            write (o_unit,*)
            write (o_unit,*) '    RR,    Rg ht,  width, spacing, ',     &
                            &&'orient., dike spacing'
            write (o_unit,1700) 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 l = 1,nslay(s)
               write (o_unit,2000) l,ahrwcw(l,s),ahrwca(l,s)
            end do
            write (o_unit,*) 'Hourly water contents - ahrwc0'
            write (o_unit,1600) (ahrwc0(h,s),h=1,6)
            write (o_unit,1600) (ahrwc0(h,s),h=7,12)
            write (o_unit,1600) (ahrwc0(h,s),h=13,18)
            write (o_unit,1600) (ahrwc0(h,s),h=19,24)
 
         end do
 
         write (o_unit,*)
         write (o_unit,*) '+++ WEATHER +++'
         write (o_unit,*)
         write (o_unit,*) ' anemht    awwzo   wzoflg '
         write (o_unit,*) anemht,awzzo,wzoflg
         write (o_unit,*) ' wind dir (deg) and max wind speed (m/s)'
         write (o_unit,1200) awadir,awudmx
 
         write (o_unit,*)
         write (o_unit,*) 'Wind speeds (m/s) - ',ntstep,' intervals'
         k = 6
         j = 1
         do
            if (k<ntstep) then
               write (o_unit,1600) (awu(i),i=j,k)
               j = k + 1
               k = k + 6
               cycle
            else
               k = ntstep
               write (o_unit,1600) (awu(i),i=j,k)
            end if
            write (o_unit,*)
            write (o_unit,*) 'END OF INPUTS'
            exit
         end do
 
      end if
            !(o_unit .ne. 6)
 
      !close (i_unit)
      end subroutine erodin
!*==getline.spg  processed by SPAG 6.70Rc at 15:33 on 10 Dec 2012
!*------------------   SPAG Configuration Options   --------------------
!*--0323,76 000101,-1 000000102011332010100002000000210211210,136 10 --
!*--1100000011112111000000000000,10,10,10,10,10,10,900,100 200000000 --
!*--000000010000000000000,72,72 73,42,38,33 00011112110000100000000  --
!*----------------------------------------------------------------------
 
!**********************************************************************
      function getline(i_unit)
      use c_flags
      implicit none
!*--GETLINE1010
!
!*** Start of declarations rewritten by SPAG
!
! Dummy arguments
!
      integer :: i_unit
      character(*) :: getline
      intent (in) i_unit
!
! Local variables
!
      integer,save :: dataline,linecount
!
!*** End of declarations rewritten by SPAG
!
 
 
 
      do
 
 
         read (i_unit,'(A)') getline
         linecount = linecount + 1
         if (btest(debugflg,0)) write (6,*) linecount,': ',trim(getline)
 
         if (getline(1:1)/='#') then
 
            dataline = dataline + 1
            if (btest(debugflg,1)) write (6,*) linecount,':',dataline,  &
               &': ',trim(getline)
            exit
         end if
      end do
 
      end function getline