subroutine edeu(km,nl450,aw,bw,ac,bc,WtoQ,facirr,Ed,Es,H,P,
     * acdom,rmud,tirrq,ihra,avgq,deltaE)
 
!  Model of irradiance in the water column.  Accounts for three 
!  irradiance streams:
 
!  Edz = direct downwelling irradiance
!  Esz = diffuse downwelling irradiance
!  Euz = diffuse upwelling irradiance
 
!  Propagation is done in energy units, tests are done in quanta,
!  final is quanta for phytoplankton growth.
 
#include "definebio.h"
#include "comlte.h"
      real Ed(nlt),Es(nlt)
      real WtoQ(nlt)
      real H(km)
      real P(km,ntyp)
      real acdom(km,nlt)
      real tirrq(km),avgq(km)
      real deltaE(km)
      real aw(nlt),bw(nlt)
      real ac(nchl,nlt),bc(nchl,nlt)
      real Edz(nlt,km),Esz(nlt,km)
      real Euz(nlt,km)
      real Edtop(nlt),Estop(nlt)
      real fchl(nchl)
      real facirr(nh,nch,5,ncd)
      real bbc(5)
      data bbc /0.002, 0.00071, 0.0032, 0.00071, 0.0029/
      data bbw /0.5/       !backscattering to forward scattering ratio
      data Dmax /500.0/    !depth at which Ed = 0
 
!  Constants and initialize
      rmus = 1.0/0.83            !avg cosine diffuse down
      tirrq = 0.0
      deltaE = 0.0
 
       Ebot = 0.0
       do nl = 1,nlt
        Edtop(nl) = Ed(nl)
        Estop(nl) = Es(nl)
        Ebot = Ebot + (Ed(nl)+Es(nl))
       enddo
!  Convert to quanta: divide by Avos # to get moles quanta; then mult by
!  1E6 to get uM or uEin
       Ebotq = 0.0
       do nl = npst,npnd   !PAR range only 350-700nm
        Ebotq = Ebotq + (Edtop(nl)+Estop(nl))*WtoQ(nl)*1.0E6
       enddo
       do k = 1,km
        Etop = Ebot
        Etopq = Ebotq
        zd = min(Dmax,H(k))
        zirr = 0.0
        zirrq = 0.0
        do nl = 1,nlt
         Edz(nl,k) = 0.0
         Esz(nl,k) = 0.0
         Euz(nl,k) = 0.0
         actot = 0.0
         bctot = 0.0
         bbctot = 0.0
         do n = 1,nchl
          actot  = actot  + P(k,nnut+n)*ac(n,nl)
          bctot  = bctot  + P(k,nnut+n)*bc(n,nl)
          bbctot = bbctot + P(k,nnut+n)*bbc(n)*bc(n,nl)
         enddo
         a  = aw(nl) + acdom(k,nl) + actot
         bt = bw(nl) + bctot
         bb = bbw*bw(nl) + bbctot
         bb = max(bb,0.0002)
         if (Edtop(nl) .ge. 1.0E-4 .or. Estop(nl) .ge. 1.0E-4)then
          call radmod(zd,Edtop(nl),Estop(nl),rmud,a,bt,bb,
     *      Dmax,Edz(nl,k),Esz(nl,k),Euz(nl,k))
         endif
         Edtop(nl) = Edz(nl,k)
         Estop(nl) = Esz(nl,k)
         zirr = zirr + (Edz(nl,k)+Esz(nl,k)+Euz(nl,k))
        enddo
        Ebot = zirr
        deltaE(k) = Etop - Ebot
        do nl = npst,npnd
         zirrq = zirrq + (Edz(nl,k)+Esz(nl,k)+Euz(nl,k))*WtoQ(nl)*1.0E6
        enddo
        Ebotq = zirrq
        ih = nint(H(k))
        acdom450 = acdom(k,nl450)
        ih = min(ih,nh)
        icd = nint((alog10(acdom450)*100.0+300.0)/10.0) + 1
        icd = max(icd,1)
        icd = min(icd,ncd)
        chl = 0.0
        do n = nnut+1,npe
         chl = chl + P(k,n) 
        enddo
        chlm = max(chl,0.001)
        ich = nint((alog10(chlm)*100.0+300.0)/10.0) + 1
        ich = max(ich,1)
        ich = min(ich,nch)
        chlm = max(chl,0.00001)
        do n = 1,nchl
         fchl(n) = P(k,nnut+n)/chlm
        enddo
        fac = 0.0
        do n = 1,nchl
         fac = fac + fchl(n)*facirr(ih,ich,n,icd)
        enddo
 
        tirrq(k) = fac*((Etopq+Ebotq)*0.5)*rmus
       enddo
 
!  Irradiance summary loops
      ihra = ihra + 1
      do k = 1,km
       avgq(k) = avgq(k) + tirrq(k)
      enddo
 
      return
      end