subroutine molini
c     =================
c
c     Initialization of the molecular equilibrium
c
      INCLUDE 'PARAMS.FOR'
      INCLUDE 'MODELP.FOR'
      common/moltst/pfmol(600,mdepth),anmol(600,mdepth),
     *              pfato(100,mdepth),anato(100,mdepth),
     *              pfion(100,mdepth),anion(100,mdepth)
      dimension hpo(mdepth)
c
      aeinit=1.0
c
      do 10 id=1,nd
         t=temp(id)
         tln=log(t)*1.5
         thl=11605./t
         t32=exp(tln)

         do i=1,MMOLEC
            rrmol(i,id)=0.
         end do
 
         hpo(id)=DENS(ID)/WMM(ID)/YTOT(ID)

         if(t.gt.tmolim) go to 10
         HPOP=DENS(ID)/WMM(ID)/YTOT(ID)

         an=dens(id)/wmm(id)+elec(id)
         aeinit=0.1*an
         if(t.lt.4000.) aeinit=0.01*an
         call moleq(id,t,an,aeinit,ane,0)
         

c     next initial guess will be the last ane determined for
c     previous depth point

         aeinit=ane
c
         if (id.eq.idstd) then
            write(6,600)
            nmol=nmolec
            if(id.eq.1) nmol=32
            do i=1,nmol
               write(6,601) i, cmol(i), rrmol(i,id), rrmol(i,id)/hpop
            end do
         end if
 600     format(/ 'Molecular number densities at the standard depth'/)
 601     format(i4,1x,A8,1x,1pe12.2,1x,e12.2)

   10 continue


c     update atomic populations once molecular densities are calculated

      if(imode.lt.-4) then
      do i=1,nlevel
         iat=numat(iatm(i))
         ion=iz(iel(i))
         ii=nfirst(iel(i))
         ener=(enion(ii)-enion(i))/bolk
         if((enion(i).eq.0).and.(ilk(i).gt.0)) then
            ener=0.
            ion=ion+1
         end if
         if(ifwop(i).ge.0) then
            do id=1,nd
               popul(i,id)=rrr(id,ion,iat)*g(i)
     *              *exp(-ener/temp(id))
               if(iat.eq.1.and.ion.eq.0) popul(i,id)=anhm(id)
            end do
         endif
      end do
      end if
c
      return
      end