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
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'
c     include 'c1gen.inc'
      include 'b1glob.inc'
      include 'b1geom.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'
c     include 'w1pavg.inc'
c
c     + + + LOCAL COMMON BLOCKS + + +
     
c
*$reference
c
c     +++ LOCAL VARIABLES +++
      integer i,j,k
      real f(mntime), wfcalm, wuc, w0k, step
      character*80 line
c
c     + + + LOCAL VARIABLE DEFINITIONS + + +
c     i, j, k = do-loop indices
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     +++ SUBROUTINES CALLED +++
c     sbr2
c
c     +++ END SPECIFICATIONS +++
c
c      read  +++ MAIN  COMMON BLOCKS +++
c
      call sbr2(i_unit, 18)

c     m1sim.inc  +++ SIMULATION REGION +++
      read (i_unit,*) am0jd, ntstep, amalat, amalon, amzele

      call sbr2(i_unit, 3)

c     m1flag.inc
      read (i_unit,*) am0ifl, am0eif, am0cgf, am0cif, am0dif
      call sbr2(i_unit, 1)

      read (i_unit,*) am0hfl, am0sfl, am0tfl, am0cfl, am0dfl, am0efl
      call sbr2(i_unit, 3)
c     m1subr.inc

      read (i_unit,*) nsubr, am0csr, (am0irr(i),i=1,nsubr),
     &  (amrslp(i), i=1,nsubr)
      call sbr2(i_unit, 3)
c      
c     m1geo.inc
      read (i_unit,*) nacctr, nbr, amasim
            call sbr2(i_unit, 1)
    
      read (i_unit,*) (((amxar (i,j,k), i=1,2), j=1,2), k=1,nacctr)
      
c
      if (nbr .gt. 0) then
         call sbr2(i_unit, 1)
         do 3  k = 1, nbr
    3    read (i_unit,*) ((amxbr(i,j,k), i=1,2), j=1,2)
         call sbr2(i_unit, 1)
c      
  
         do 4 k = 1,nbr
    4    read (i_unit,*)  amzbr(k), ampbr(k), amxbrw(k)
      else
         call sbr2(i_unit, 4)
      endif
c         
      call sbr2(i_unit, 1)
      read (i_unit,*) ((amxsim (i,j), i=1,2),j=1,2)
      call sbr2(i_unit, 1)
      read (i_unit,*) (((amxsr (i,j,k), i=1,2),j=1,2),k=1,nsubr)
c ^^^ tmp
c      go to 500 
c     read   +++ BIOMASS +++
      call sbr2(i_unit, 3)

c    not used
c     c1gen.inc
c      do 3  i=1, nsubr
c    3 read (i_unit,*) ac0rg(i), acxrow(i), acdpop(i)
c      call sbr2(i_unit, 3)

c     b1glob.inc
      do 5   i=1,nsubr
    5 read (i_unit,*) abzht(i), abz(i), abzrtd(i), abmst(i), abmrt(i),
     & abrahi(i)
      call sbr2(i_unit, 3)

c     b1geom.inc
      do 10  i=1, nsubr
   10 read (i_unit,*) abrsai(i), abrlai(i),
     & (abrsaz(j,i), j=1,1), (abrlaz(j,i), j=1,1)
c
      call sbr2(i_unit, 1)
      do 11 i=1, nsubr
   11  read (i_unit,*) abffcv(i), abfscv(i), abftcv(i)
    
c     read   +++ SOIL +++
      call sbr2(i_unit, 3)
      
c     s1layr.inc
      read (i_unit,*) (nslay(i), i=1,nsubr)
      call sbr2(i_unit, 1)
      do 20 j=1,nsubr
   20 read (i_unit,*) (aszlyt(i,j),i=1,nslay(j))
      call sbr2(i_unit, 3)
c     s1phys.inc
      do 30 j=1,nsubr
   30 read (i_unit,*) (asdblk(i,j), i=1,nslay(j))
      call sbr2(i_unit, 3)
     
c     s1dbh.inc
      do 40  j=1, nsubr
   40 read (i_unit,*) (asfsan(i,j), i=1,3 )
c     i=1,nslay(j))
      call sbr2(i_unit, 1)

      do 50 j=1,nsubr
   50 read (i_unit,*) (asfsil(i,j), i=1,nslay(j))
      call sbr2(i_unit, 1)
      do 60 j=1,nsubr
   60 read (i_unit,*)(asfcla(i,j), i=1, nslay(j))
      call sbr2(i_unit, 1)
      do 70 j=1,nsubr
   70 read (i_unit,*)(asvroc(i,j), i=1, nslay(j))
      call sbr2(i_unit, 3)
     
c     s1agg.inc
      do 75 j=1,nsubr
      k = j
   75 read (i_unit,*) (asdagd(i,j), i=1,nslay(k))
      call sbr2(i_unit, 1)
      do 80 j=1, nsubr
   80 read (i_unit,*) (aseags(i,j),i=1,nslay(j))
      call sbr2(i_unit, 1)
     
      do 90 j=1,nsubr
   90 read (i_unit,*)(aslagm(i,j),i=1,nslay(j))
      call sbr2(i_unit, 1)

      do 100 j=1,nsubr
  100 read (i_unit,*)(aslagn(i,j), i=1,nslay(j))
      call sbr2(i_unit, 1)

      do 110 j=1,nsubr
  110 read (i_unit,*)(aslagx(i,j), i=1,nslay(j))
      call sbr2(i_unit, 1)
      do 120 j=1,nsubr
  120 read (i_unit,*)(as0ags(i,j), i=1,nslay(j))
      call sbr2(i_unit, 3)
    
c     s1surf.inc
      do 130 i=1,nsubr
  130 read (i_unit,*)aszcr(i), asfcr(i), asmlos(i), asflos(i),
     &  asdcr(i), asecr(i)
      call sbr2(i_unit, 1)
      do 135 i = 1,nsubr
  135 read (i_unit,*) asfald(i), asfalw(i)
      call sbr2(i_unit, 3)

c     s1sgeo.inc
      do 140 i=1,nsubr
  140    read(i_unit,*) aszrgh(i), asxrgw(i), asxrgs(i), asargo(i),
     &          asxdks(i), aslrr(i)
     
      call sbr2(i_unit, 3)
     
c     read   +++ HYDROLOGY +++

c     h1db1.inc
      do 160 j = 1,nsubr
  160 read (i_unit,*) (ahrwcw(i,j), i=1,nslay(j))
      call sbr2(i_unit, 1)

      do 170 j = 1,nsubr
  170 read (i_unit,*) (ahrwcf(i,j), i=1,nslay(j))
      call sbr2(i_unit, 1)

      do 180 j = 1,nsubr
  180 read (i_unit,*) (ahrwca(i,j), i=1,nslay(j))
      call sbr2(i_unit, 1)
      do 190 j = 1,nsubr
  190 read (i_unit,*)(ahrwcs(i,j), i=1,nslay(j))
       call sbr2(i_unit, 1)

      do 195 j = 1,nsubr
  195 read (i_unit,*)(ahrwc0(i,j), i=1,24)
      call sbr2(i_unit, 2)
      
      read (i_unit,*) (ahzsnd(i), i = 1,nsubr)
      call sbr2(i_unit, 3)
c     read   +++ WEATHER +++
      
c     w1pavg.inc  skip
c     read (i_unit,*) awtmmx,awtmmn,awadir,awztpt,awtpav,awepir
      call sbr2(i_unit, 1)
c     read (i_unit,*) awupav, awnurt, aweuet
      call sbr2(i_unit, 2)
      call sbr2(i_unit, 3)

c     w1wind.inc
c
c     wind data inputs as the Weibull paramters
c     (wfcalm, wuc, w0k) is indicated by code ntstep = 99
      
      if (ntstep .eq. 99) then
         read (i_unit,*) awadir
         call sbr2(i_unit, 1)
         read (i_unit,*) wfcalm, wuc, w0k, ntstep
c ^^^ tmp
c         write(*,*) 'wfcalm=', wfcalm, wuc, w0k, ntstep 
c        
c     calculate daily max wind speed
         awudmx = wuc*(-log((1-0.99)/(1-wfcalm)))**(1/w0k)
c ^^^ tmp
c         write(*,*) 'awudmx=', awudmx
c         
c     calculate period wind speeds
         step = ntstep
         do 198 i= 1, ntstep
c          find center of each step and add empirical last term from 
c                  file ntstep.mcd         
           f(i) = (1./(2.0*step)) + ((i-1)/step) +0.3/(step*w0k)
           
c ^^^ tmp
c        write(*,*) 'i=', i, '  f(i)=', f(i), 'ntstep=', ntstep
c          
c            to prevent out-of-range           
             if (f(i) .lt. wfcalm) then
                 f(i) = wfcalm
             endif    
           awu(i) = wuc*(-log((1-f(i))/(1-wfcalm)))**(1/w0k)
c ^^^tmp
c         write(*,*) 'f(i)=', f(i), '  awu(i)=', awu(i)           
  198 continue
          awudmx = awu(ntstep)
c  
c     when ntstep .ne. 99 read input wind period data directly
      else
         read (i_unit,*) awadir, awrmxn, awhrmx, awudmx, awudmn, awudav
          call sbr2(i_unit, 3)
          read (i_unit,*) (awu(i),i=1,ntstep)
      endif 
      
c     + + + OUTPUT SECTION + + + 
c
c     can skip output, if desired
c     go to 900
c
  200 format (1x, f6.2, 6f8.2)
  350 format (1x, i3, 6f6.0)
  300 format (1x, 8i5)
  400 format (1x, i3, 1x, l6, 4x, i3, 4x, i3, 4x, i3, 4x, f6.0)
  
c ^^^ tmp
c  500 continue
c   
      write (o_unit,*)  '      REPORT OF INPUTS (read/write by sbinpt) '
      write (o_unit,*)  '+++ SIMULATION REGION +++'
c 
      write (o_unit,*) 'ntstep  am0eif  nsubr  nacctr  nbr  amasim'
      write (o_unit,400) ntstep, am0eif, nsubr, nacctr, nbr, amasim
c
      if (nbr .gt. 0) then
         write (o_unit,*)  
         write (o_unit,*) 'amxbr    (start and end of each barrier)  '
         write (o_unit,200) (((amxbr(i,j,k), i=1,2), j=1,2), k=1,nbr)
       
         write (o_unit,*) 'amzbr     ampbr   amxbrw '
         do 750 k = 1, nbr
  750    write (o_unit,200) amzbr(k), ampbr(k), amxbrw(k)
       else 
         write (o_unit,*) ' no barriers'
       endif  
c      write (o_unit,*) 
      write (o_unit,*) 'amasim  amxsim  (angle & corners of sim region)'
      write (o_unit,200) amasim,  ((amxsim(i,j), i=1,2), j=1,2)
c
c ^^^ tmp
c      go to 900 
c
      write (o_unit,*)  
      write(o_unit,*) '+++ BIOMASS +++'
      write (o_unit,*) 'abzht     abrsai  abrlai  abffcv  '
      do 800 i=1,nsubr
  800 write (o_unit,200) abzht(i), abrsai(i), abrlai(i), abffcv(i)
c
      write (o_unit,*) 
      write (o_unit,*) '+++ SOIL +++ '
      write (o_unit,*) 'nslay      aszlyt    '
      do 810 j=1,nsubr
  810 write (o_unit,350) nslay(j), (aszlyt(i,j), i=1,nslay(j))
c
c      write (o_unit,*) 
      write (o_unit,*) 'asdblk '
      do 820 j=1,nsubr
  820 write (o_unit,200) (asdblk(i,j), i=1,nslay(j)) 
c
c      write (o_unit,*)  
      write (o_unit,*) 'For top layer and subregion 1 '
      write (o_unit,*) 'asfsan    asfsil  asfcla   asvroc '
      write (o_unit,200) asfsan(1,1),asfsil(1,1),asfcla(1,1),asvroc(1,1)
c
      write (o_unit,*) 
      write (o_unit,*) 'asdagd   aseags  aslagm  aslagn  aslagx  as0ags'
      write (o_unit,200) asdagd(1,1), aseags(1,1), aslagm(1,1),
     & aslagn(1,1), aslagx(1,1), as0ags(1,1) 
c 
      write (o_unit,*)     
      write (o_unit,*) 'asfcr     asmlos  asflos  asdcr  asecr '
      do 830 i= 1,nsubr
  830 write (o_unit,200) asfcr(i),asmlos(i),asflos(i),asdcr(i),asecr(i)
c 
      write (o_unit,*)  
      write (o_unit,*) 'aszrgh    asxrgw  asxrgs asargo  asxdks  aslrr' 
      do 840 i=1,nsubr    
  840 write (o_unit,200) aszrgh(i), asxrgw(i), asxrgs(i), asargo(i),
     & asxdks(i), aslrr(i)    
c
      write (o_unit,*) 
      write (o_unit,*) '+++ HYDROLOGY +++ '    
      write (o_unit,*) 'ahrwcw    ahrwca  ahrwc0  ahzsnd '
      do 850 i=1,nsubr
  850 write (o_unit,200) ahrwcw(1,i), ahrwca(1,i),ahrwc0(24,i),ahzsnd(i)
c
      write (o_unit,*)  
      write(o_unit,*) '+++ WEATHER +++'
      write (o_unit,*) ' awadir   awudmx'
      write(o_unit,200)  awadir, awudmx
c      
      write (o_unit,*) 'awu (ntstep) '
      k = 6
      j = 1
  860 if((k) .lt. ntstep) then
        write (o_unit,200) (awu(i), i=j,k)
        j = k+1
        k = k+6
        go to 860
      else
        k = ntstep
        write (o_unit,200) (awu(i), i=j,k)
      endif      
      write (o_unit,*) 'END OF INPUTS'
c
c     use continue if skip report section      
  900 continue
c
      close (i_unit)
      return
      end
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

c**********************************************************************
c**********************************************************************
      subroutine sbr2(i_unit, n)
c**********************************************************************
c
c     +++ PURPOSE +++
c     skips 'n' comment lines in the input file
c
	  integer i_unit
      integer n,i
      do 10 i = 1,n 
      read (i_unit,*)
   10 continue
      return
      end
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++