!$Author: joelevin $
!$Date: 2011-03-24 11:33:26 -0500 (Thu, 24 Mar 2011) $
!$Revision: 11724 $
!$HeadURL: https://svn.weru.ksu.edu/weru/weps1/trunk/weps.src/src/lib_erosion/daily_erodout.for $

      subroutine daily_erodout (o_unit, o_E_unit)

!     +++  PURPOSE +++

!     To print output desired from standalone EROSION submodel

!     +++ ARGUMENT DECLARATIONS +++

      integer o_unit, o_E_unit

      integer i, j
      real aegt, aegtss, aegt10
      real tt, lx, ly
      real topt,topss, top10, bott, botss, bot10, ritt, ritss, rit10
      real lftt, lftss, lft10, tot, totbnd

      integer initflag
      save    initflag
      integer yr, mo, da

!     + + + GLOBAL COMMON BLOCKS + + +
      include 'p1werm.inc'
      include 'm1geo.inc'
      include 'm1flag.inc'

!     + + + LOCAL COMMON BLOCKS
      include 'erosion/e2erod.inc'
      include 'erosion/m2geo.inc'
!
      integer x, y


      integer :: dt(8)
      character(len=3) :: mstring
      common / datetime / dt, mstring

!     ++++ ARGUMENT DEFINITIONS +++

!     +++ SUBROUTINES CALLED+++
!     plotout.for

!     ++++ LOCAL VARIABLES +++

!     +++ END SPECIFICATIONS +++


!     Calculate Averages Crossing Borders
!      top border
       aegt   = 0.0
       aegtss = 0.0
       aegt10 = 0.0
       j = jmax
       do 1 i = 1, imax-1
         aegt    = aegt   + egt(i,j)
         aegtss  = aegtss + egtss(i,j)
         aegt10  = aegt10 + egt10(i,j)
    1  continue
!      calc. average at top border
       topt  = aegt/(imax-1)
       topss = aegtss/(imax-1)
       top10 = aegt10/(imax-1)

!      bottom border
       aegt   = 0.0
       aegtss = 0.0
       aegt10 = 0.0
       j = 0
       do 2 i = 1, imax-1
         aegt    = aegt   + egt(i,j)
         aegtss  = aegtss + egtss(i,j)
         aegt10  = aegt10 + egt10(i,j)
    2  continue
!      calc. average at bottom border
        bott  = aegt/(imax-1)
        botss = aegtss/(imax-1)
        bot10 = aegt10/(imax-1)

!     right border
       aegt   = 0.0
       aegtss = 0.0
       aegt10 = 0.0
       i = imax
       do 3 j = 1, jmax-1
         aegt    = aegt   + egt(i,j)
         aegtss  = aegtss + egtss(i,j)
         aegt10  = aegt10 + egt10(i,j)
    3  continue
!      calc. average at right border
        ritt  = aegt/(jmax-1)
        ritss = aegtss/(jmax-1)
        rit10 = aegt10/(jmax-1)
!
!     left border
       aegt   = 0.0
       aegtss = 0.0
       aegt10 = 0.0
       i = 0
       do 4 j = 1, jmax-1
         aegt    = aegt   + egt(i,j)
         aegtss  = aegtss + egtss(i,j)
         aegt10  = aegt10 + egt10(i,j)
    4  continue
!      calc. average at left border
        lftt   = aegt/(jmax-1)
        lftss  = aegtss/(jmax-1)
        lft10  = aegt10/(jmax-1)

!     calculate averages of inner grid points
      aegt   = 0.0
      aegtss = 0.0
      aegt10 = 0.0
      do 5 j=1,jmax-1
       do 5 i= 1, imax-1
        aegt= aegt + egt(i,j)
        aegtss = aegtss + egtss(i,j)
        aegt10 = aegt10 + egt10(i,j)
    5 continue
      tt     = (imax-1)*(jmax-1)
      aegt   = aegt/tt
      aegtss = aegtss/tt
      aegt10 = aegt10/tt

!    calculate comparision of boundary and interior losses
      lx = amxsim(1,2) - amxsim(1,1)
      ly = amxsim(2,2) - amxsim(2,1)
      tot = aegt*lx*ly
      totbnd = (topt + bott + topss + botss)*lx +                       &
     &         (ritt + lftt + ritss + lftss)*ly


      if (btest(am0efl,1)) then
        if (initflag .eq. 0) then  !write header to files

            write (o_unit,*) 'Grid cell output from daily_erodout.for'
            write (o_unit,*)

           ! Print date of Run
12          format(1x,'Date of run: ',a3,' ',i2.2,', ',i4,' ',          &
     &          i2.2,':',i2.2,':',i2.2)
            write (o_unit,12) mstring, dt(3), dt(1), dt(5), dt(6), dt(7) 
            write (o_unit,*)

            write (o_unit,fmt="(1x,a,5f10.2)") "<field dimensions>",    &
     &                   amasim,((amxsim(x,y),x=1,2),y=1,2)
            write (o_unit,*)
            write (o_unit,*)                                            &
     &                 'Total grid size: (', imax+1,',', jmax+1, ')   ',&
     &                 'Inner grid size: (', imax-1,',', jmax-1, ')'
            write (o_unit,*)
            initflag = initflag + 1
        endif

        call caldatw (da, mo, yr) !get day, month and year
        write (UNIT=o_unit,FMT="(' ',i5,i3,i3,' ')",ADVANCE="YES")      &
     &         yr, mo, da
        write (o_unit,*)

      write (o_unit,*)                                                  &
     & '  top(i=1,imax-1,j=jmax) ',                                     &
     & 'bottom(i=1,imax-1,j=0) ',                                       &
     & 'right(i=imax,j=1,jmax-1) ',                                     &
     & 'left(i=0, j=1,jmax-1) '
      write (o_unit,10)  (egt(i,jmax)+egtss(i,jmax), i = 1, imax-1)
      write (o_unit,10)  (egt(i,0)+egtss(i,0), i = 1, imax-1)
      write (o_unit,10)  (egt(imax,j)+egtss(imax,j), j = 1, jmax-1)
      write (o_unit,10)  (egt(0,j)+egtss(0,j), j = 1, jmax-1)

      write (o_unit,*)
      write (o_unit,7)
      write (o_unit,*)                                                  &
     & '  top(i=1,imax-1,j=jmax) ',                                     &
     & 'bottom(i=1,imax-1,j=0) ',                                       &
     & 'right(i=imax,j=1,jmax-1) ',                                     &
     & 'left(i=0, j=1,jmax-1) '
      write (o_unit,10)  (egt(i,jmax), i = 1, imax-1)
      write (o_unit,10)  (egt(i,0), i = 1, imax-1)
      write (o_unit,10)  (egt(imax,j), j = 1, jmax-1)
      write (o_unit,10)  (egt(0,j), j = 1, jmax-1)

      write (o_unit,*)
      write (o_unit,8)
      write (o_unit,*)                                                  &
     & '  top(i=1,imax-1,j=jmax) ',                                     &
     & 'bottom(i=1,imax-1,j=0) ',                                       &
     & 'right(i=imax,j=1,jmax-1) ',                                     &
     & 'left(i=0, j=1,jmax-1) '
      write (o_unit,10)  (egtss(i,jmax), i = 1, imax-1)
      write (o_unit,10)  (egtss(i,0), i = 1, imax-1)
      write (o_unit,10)  (egtss(imax,j), j = 1, jmax-1)
      write (o_unit,10)  (egtss(0,j), j = 1, jmax-1)

      write (o_unit,*)
      write (o_unit,9)
      write (o_unit,*)                                                  &
     & '  top(i=1,imax-1,j=jmax) ',                                     &
     & 'bottom(i=1,imax-1,j=0) ',                                       &
     & 'right(i=imax,j=1,jmax-1) ',                                     &
     & 'left(i=0,j=1,jmax-1) '
      write (o_unit,11)  (egt10(i,jmax), i = 1, imax-1)
      write (o_unit,11)  (egt10(i,0), i = 1, imax-1)
      write (o_unit,11)  (egt10(imax,j), j = 1, jmax-1)
      write (o_unit,11)  (egt10(0,j), j = 1, jmax-1)

      write (o_unit,*)
      write (o_unit,fmt="(' <grid data> | ',3('|',a))")                 &
     & 'Total Soil Loss', 'soil loss', '(kg/m^2)'
      do 19  j = jmax-1, 1, -1
      write (o_unit,10)  (egt(i,j), i = 1, imax-1)
   19 continue
      write (o_unit,fmt="(' </grid data>')")

      write (o_unit,*)
      write (o_unit,fmt="(' <grid data> | ',3('|',a))")                 &                              
     & 'Saltation/Creep Soil Loss', 'salt/creep soil loss', '(kg/m^2)'
      do 29  j = jmax-1, 1, -1
      write (o_unit,10)  (egt(i,j)-egtss(i,j), i = 1, imax-1)
   29 continue
      write (o_unit,fmt="(' </grid data>')")

      write (o_unit,*)
      write (o_unit,fmt="(' <grid data> | ',3('|',a))")                 &                              
     & 'Suspension Soil Loss', 'suspension soil loss', '(kg/m^2)'
      do 39  j = jmax-1, 1, -1
      write (o_unit,10)  (egtss(i,j), i = 1, imax-1)
   39 continue
      write (o_unit,fmt="(' </grid data>')")

      write (o_unit,*)
      write (o_unit,fmt="(' <grid data> | ',3('|',a))")                 &
     & 'PM10 Soil Loss', 'PM10 soil loss', '(kg/m^2)'
      do 49  j = jmax-1, 1, -1
      write (o_unit,11)  (egt10(i,j), i = 1, imax-1)
   49 continue
      write (o_unit,fmt="(' </grid data>')")

      write (o_unit,*)
      write (o_unit,*) '**Averages - Field'
      write (o_unit,*) '     Total    salt/creep      susp       PM10 '
      write (o_unit,*) '     egt                      egtss      egt10'
      write (o_unit,*) '   -----------------kg/m^2--------------------'
      write (o_unit,15)    aegt, aegt-aegtss, aegtss, aegt10
      write (o_unit,*)
      write (o_unit,*) '**Averages - Crossing Boundaries '
      write (o_unit,*) 'Location      Total  Salt/Creep   Susp    PM10'
      write (o_unit,*) '--------------------kg/m----------------------'
      write (o_unit,21) topt+topss, topt, topss, top10
      write (o_unit,22) bott+botss, bott, botss, bot10
      write (o_unit,23) ritt+ritss, ritt, ritss, rit10
      write (o_unit,24) lftt+lftss, lftt, lftss, lft10
      write (o_unit,*)
      write (o_unit,*) '   Comparision of interior & boundary loss'
      write (o_unit,*) '      interior       boundary    int/bnd ratio'
      if( totbnd.gt.1.0e-9 ) then
          write (o_unit,16) tot, totbnd, tot/totbnd
      else
          !Boundary loss near or equal to zero
          write (o_unit,16) tot, totbnd, 1.0e-9
      end if

!^^^tmp out
!      write (o_unit,*) 'lx=', lx, 'ly=', ly,'tot',tot
!^^^end tmpout

!     additional output statements for easy shell script parsing
      write (o_unit,*)
!     write deposition as positive numbers
      write (o_unit,17) aegt, aegt-aegtss, aegtss, aegt10
   17 format (' repeat of total, salt/creep, susp, PM10:', 3f12.4,f12.6)

!     output formats

    6 format (1x,'  Passing Border Grid Cells - Total  egt+egtss(kg/m)')
    7 format (1x,'  Passing Border Grid Cells - Salt/Creep   egt(kg/m)')
    8 format (1x,'  Passing Border Grid Cells - Suspension egtss(kg/m)')
    9 format (1x,'  Passing Border Grid Cells - PM10       egt10(kg/m)')

   50 format (1x,'  Leaving Field Grid Cells - Total       egt(kg/m^2)')
   60 format (1x,'  Leaving Field Grid Cells - Salt/Creep egt-egtss(kg/m&
     &^2)')
   70 format (1x,'  Leaving Field Grid Cells - Suspension egtss(kg/m^2) &
     &')
   80 format (1x,'  Leaving Field Grid Cells - PM10      egt10(kg/m^2)')

   10 format (1x, 500f12.4)
   11 format (1x, 500f12.6)
   15 format (1x, 3(f12.4,2x), f12.6)
   16 format (1x, 2(f13.4,2x),2x, f13.4)
   21 format (1x, 'top   ', 1x, 4(f9.2,1x))
   22 format (1x, 'bottom', 1x, 4(f9.2,1x))
   23 format (1x, 'right ', 1x, 4(f9.2,1x))
   24 format (1x, 'left  ', 1x, 4(f9.2,1x))

      endif !if (btest(am0efl,1)) then

      !Erosion summary - total, salt/creep, susp, pm10
      !(deposition values are positive - erosion values are negative)
      if (btest(am0efl,0)) then
         call caldatw (da, mo, yr) !get day, month and year
         write (UNIT=o_E_unit,FMT="(' ',i5,i3,i3,' ')",ADVANCE="NO")    &
     &         yr, mo, da
         write (UNIT=o_E_unit,FMT="(4(f12.6),' ')",ADVANCE="YES")       &
     &          aegt, (aegt-aegtss), aegtss, aegt10
      endif

!    end of plot section

      return
      end