SUBROUTINE AERATE
  USE GLOBAL;USE MAIN; USE KINETIC; USE TRANS; USE SCREENC
  
  Real, Allocatable, Dimension (:) ::  dzmulta,smass,atimon,atimoff,cumdomass
  Real, Allocatable, Dimension (:,:) ::  dzmult
  Integer, Allocatable, Dimension (:) ::   iaseg,ktopa,kbota,iprb,kprb
  Real, Allocatable, Dimension (:) :: dooff,doon
  Real, Allocatable, Dimension (:) :: actual_mass
  Logical, Allocatable, Dimension (:) :: aerateo2
  integer naer, nlayers, ktop,kbot
  character(16) :: conaer
  SAVE
  
  allocate(dzmult(kmx,imx))
   dzmult=1.0
   open(9540,file='w2_aerate.npt',status='old')
   READ (9540,'(//i8,a16)')NAER,conaer

     if(naer.eq.0)naer=1    
     allocate(dzmulta(naer),iaseg(naer),ktopa(naer),kbota(naer),smass(naer),atimon(naer),atimoff(naer),dooff(naer),doon(naer),iprb(naer),kprb(naer),actual_mass(naer),aerateo2(naer),cumdomass(naer))
     dzmulta=0.0; actual_mass=0.0; cumdomass=0.0
     read(9540,1013)
1013 FORMAT(/)
      do i=1,naer
        Read(9540,'(i8,i8,i8,f8.0,f8.0,f8.0,3f8.0,2i8)')iaseg(i),ktopa(i),kbota(i),smass(i),atimon(i),atimoff(i),dzmulta(i),dooff(i),doon(i),iprb(i),kprb(i)
      end do

	close(9540)

if(restart_in)then
  open(9540,file=conaer,POSITION='APPEND')
        JDAY1=0.0
        REWIND (9540)
        READ   (9540,'(//)')
        DO WHILE (JDAY1 < JDAY)
          READ (9540,'(F9.0)',END=106) JDAY1
        END DO
        BACKSPACE (9540)
        106     JDAY1=0.0
else
  open(9540,file=conaer,status='unknown')
  write(9540,'(a,i3,a)')'Output file for Aeration Input with',naer,' input(s).'
  write(9540,'(a114)')'JDAY,InstMassRate#1(kgO2/d),CumMass#1(kgO2),DOProbe#1(mg/l),InstMassRate#2(kgO2/d),CumMass#2(kgO2),DOProbe#2(mg/l)'
endif
  dzmult=1.0  ! always reset mixing coefficient for aeration SW IPC 2/01/01
 RETURN
 
 ENTRY DZAERATE   ! from w2_ main code
    DO I=1,NAER
    DZ(KTOPA(I):KBOTA(I),IASEG(I))=DZ(KTOPA(I):KBOTA(I),IASEG(I))*DZMULT(KTOPA(I):KBOTA(I),IASEG(I))
    ENDDO   
 RETURN
    
 ENTRY AERATEMASS   ! from wq_constituents

! section for hypolimnetic aeration 
 DZMULT=1.0   ! ALWAYS RESET TO 1.0 IN CASE NO AERATION
              do ii=1,naer
                if(jday.ge.atimon(ii).and.jday.le.atimoff(ii))then

           ! Find branch and waterbody for iseg
                  do jjb=1,nbr
                    if(iaseg(ii).ge.us(jjb).and.iaseg(ii).le.ds(jjb))then
                      exit
                    end if
                  end do

                  if(jjb.eq.jb)then   ! If this branch doesn't have aeration skip it

                    ktop=max(ktwb(jw),ktopa(ii))
                    kbot=min(kb(iaseg(ii)),kbota(ii))

                    nlayers=kbot-ktop+1

                    aerateo2(ii)=.true.

                    if(o2(kprb(ii),iprb(ii)).gt.DOON(ii).and.o2(kprb(ii),iprb(ii)).gt.DOOFF(ii))then
                      aerateo2(ii)=.false.
                      actual_mass(ii)=0.0
                    end if

                    if(aerateo2(ii))then
                      actual_mass(ii)=smass(ii)  ! SW 12/28/01
                      cumdomass(ii)=cumdomass(ii)+smass(ii)*DLT/86400.
                      do k=ktop,kbot
            
                        dzmult(k,iaseg(ii))=dzmulta(ii)   ! this means the increase in DZ is lagged one time step
                        cssb(k,iaseg(ii),NDO)=cssb(k,iaseg(ii),NDO)+smass(ii)/(86.4*real(nlayers))
! Units of smass are in kg/day
! Typical units of CSSB: (mg/l)*(m3/s) to oonvert multiply by (1kg/10^6mg)*(1000l/1m3)*(86400s/day)==86.4
                        
                      end do          
                    end if
                  end if
                else
                aerateo2(ii)=.false.
                actual_mass(ii)=0.0
                end if
                
              end do

RETURN

ENTRY AERATEOUTPUT     ! from w2_main code

        write(9540,'(f9.3,<naer>(1x,e12.3,1x,e12.3,1x,f8.3))')jday,(actual_mass(ii),cumdomass(ii),o2(kprb(ii),iprb(ii)),ii=1,naer)

RETURN
ENTRY DEALLOCATE_AERATE
DEALLOCATE(dzmulta,iaseg,ktopa,kbota,smass,atimon,atimoff,dooff,doon,iprb,kprb,actual_mass,aerateo2,cumdomass)
DEALLOCATE(DZMULT)
close(9540)
RETURN
END