SUBROUTINE OPADD(MODE,ID,FR,ABAD,EMAD,SCAD)
C     ===========================================
C
C     Additional opacities
C     This is basically user-supplied procedure; here are some more
C     important non-standard opacity sources, namely
C     Rayleigh scattering, H- opacity, H2+ opacity, and additional
C     opacity of He I and He II.
C     Inclusion of these opacities is contolled by switches transmitted
C     by COMMON/OPCPAR - see description in START.
C
C     Input parameters:
C     MODE  - controls the nature and the amount of calculations
C           = -1 - (OPADD called from START) evaluation of relevant
C                  depth-dependent quantities (usually photoionization
C                  cross-sections, but also possibly other), which are
C                  stored in array CROS
C           = 0  - evaluation of an additional opacity, emissivity, and
C                  scattering - for procedure OPAC0
C     ID    - depth index
C     FR    - frequency 
C
C     Output: 
C
C     ABAD  - absorption coefficient (at frequency FR and depth ID)
C     EMAD  - emission coefficient (at frequency FR and depth ID)
C     SCAD  - scattering coefficient (at frequency FR and depth ID)
C
C
      INCLUDE 'PARAMS.FOR'
      INCLUDE 'MODELP.FOR'
      PARAMETER (FRAYH  =  2.463E15,
     *           FRAYHe =  5.150E15,
     *           FRAYH2 =  2.922E15,
     *           CLS    =  2.997925e18)
C
      AB0=0.
      AB1=0.
      ABAD=0.
      EMAD=0.
      SCAD=0.
C
      if(iath.gt.0) then
      N0HN=NFIRST(IELH)
      NKH=NKA(IATH)
C
      IF(MODE.GE.0) THEN
         T=TEMP(ID)
         ANE=ELEC(ID)
         HKT=HK/T
         T32=1./T/SQRT(T)
      END IF
      anh=dens(id)/(wmm(id)*ytot(id))
      anhe=rrr(id,1,2)
C
      IT=NLEVEL
C
C   -----------------------
C   HI  Rayleigh scattering
C   -----------------------
C
      IF(IRSCT.NE.0.AND.IOPHLI.NE.1.AND.IOPHLI.NE.2) THEN
         X=1.D0/(CLS/MIN(FR,FRAYH))**2
         SG=(5.799E-13+(1.422E-6+2.784*X)*X)*X*X
c        ABAD=POPUL(N0HN,ID)*SG 
         SCAD=POPUL(N0HN,ID)*SG
         scad=anh*sg
      END IF
      IF(IOPHMI.NE.0) THEN
C
C   ----------------------------
C   H-  bound-free and free-free
C   ----------------------------
C     Note: IOPHMI must not by taken non-zero if H- is considered
C           explicitly, because H- opacity would be taken twice
C
          SG=SBFHMI(FR)
          XHM=8762.9/T
          SB=1.0353E-16*T32*EXP(XHM)*POPUL(N0HN,ID)*ANE*SG
          SF=SFFHMI(POPUL(N0HN,ID),FR,T)*ANE
          AB0=SB+SF
      END IF
C
C   -----------------------
C   He I  Rayleigh scattering
C   -----------------------
C
      IF(IRSCHE.NE.0.AND.MODE.GE.0) THEN
         X=(CLS/MIN(FR,FRAYHe))**2
         CS=5.484E-14/X/X*(1.+(2.44E5+5.94E10/(X-2.90E5))/X)**2
         sg=anhe*cs
c        abad=abad+sg
         scad=scad+sg
      END IF
C
C   -----------------------
C   H2  Rayleigh scattering
C   -----------------------
C
      IF(IRSCH2.NE.0.AND.MODE.GE.0.AND.IFMOL.GT.0) THEN
         X=(CLS/MIN(FR,FRAYH2))**2
           X2=1./X/X
           CS=(8.14E-13+1.28E-6/X+1.61*X2)*X2
           sg=cs*anh2(id)
c          abad=abad+sg
           scad=scad+sg
        END IF
C
      IF(IOPH2P.GT.0.AND.IFMOL.GT.0.and.
     *   t.lt.tmolim.and.fr.lt.3.28e15) THEN
C
C   -----------------------------
C   H2+  bound-free and free-free
C   -----------------------------
C
         X=FR*1.E-15
         SG1=(-7.342E-3+(-2.409+(1.028+(-4.23E-1+
     *       (1.224E-1-1.351E-2*X)*X)*X)*X)*X)*1.602E-12/BOLK
         IT=IT+1
         X=LOG(FR)
         SG2=-3.0233E3+(3.7797E2+(-1.82496E1+(3.9207E-1-
     *        3.1672E-3*X)*X)*X)*X
         X2=-SG1/T+SG2
         SB=0.
         IF(X2.GT.-150.) SB=POPUL(N0HN,ID)*POPUL(NKH,ID)*EXP(X2)
         AB0=AB0+SB
      END IF
      end if
C
C   -----------------------------
C   He-  free-free
C   -----------------------------
C
      if(mode.ge.0.and.iophem.gt.0) then
         A=3.397D-46+(-5.216D-31+7.039D-15/FR)/FR
         B=-4.116D-42+(1.067D-26+8.135D-11/FR)/FR
         C=5.081D-37+(-8.724D-23-5.659D-8/FR)/FR
         cs=a*t+b+c/t
         sg=anhe*ane*cs
         ab0=ab0+sg
      end if
C
C   -----------------------------
C   H2-  free-free
C   -----------------------------
C
      IF(IOPH2M.NE.0.AND.MODE.GE.0.AND.IFMOL.GT.0.AND.T.LT.TMOLIM) THEN
         call h2minus(t,anh2(id),ane,fr,oph2)
         ab1=ab1+oph2
      END IF
C
C   -----------------------------
C     CH and OH continuuum opacity
C   -----------------------------
C 
      if(mode.ge.0.and.ifmol.gt.0.and.t.lt.tmolim) then
         if(iopch.gt.0) ab0=ab0+sbfch(fr,t)*anch(id)
         if(iopoh.gt.0) ab0=ab0+sbfoh(fr,t)*anoh(id)
C
C     ---------------------------
C     CIA H2-H2 opacity
C     ---------------------------
C
         if(ioh2h2.gt.0) then
            call cia_h2h2(t,anh2(id),fr,oph2)
            ab1=ab1+oph2
         end if
C
C     ---------------------------
C     CIA H2-He opacity
C     ---------------------------
C
         if(ioh2he.gt.0) then
            call cia_h2he(t,anh2(id),anhe,fr,oph2)
            ab1=ab1+oph2
         end if
C
C     ---------------------------
C     CIA H2-H opacity
C     ---------------------------
C
         if(ioh2h1.gt.0) then
            call cia_h2h(t,anh2(id),anh,fr,oph2)
            ab1=ab1+oph2
         end if
C
C     ---------------------------
C     CIA H-He opacity
C     ---------------------------
C  
         if(iohhe.gt.0) then
            call cia_hhe(t,anh,anhe,fr,oph2)
            ab1=ab1+oph2
         end if
      end if
C
C     ----------------------------------------------
C     The user may supply more opacity sources here:
C     ----------------------------------------------
C
C     Finally, actual absorption and emission coefficients

      IF(MODE.LT.0) RETURN
      X=EXP(-HKT*FR)
      X1=1.-X
      BNX=BN*(FR*1.E-15)**3*X
      ABAD=ABAD+X1*AB0+AB1
      EMAD=EMAD+BNX*(AB0+AB1/X1)
      RETURN
      END