SUBROUTINE OPADD(MODE,ICALL,IJ,ID) 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/OPCKEY - 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 OPACF1 C = 1 - the same, but also derivatives wrt temperature, C electron density, and level populations C IJ - frequency index C ID - depth index C C Output: transmitted by COMMON/OPACAD: C C ABAD - absorption coefficient (at frequency point IJ and depth ID) C EMAD - emission coefficient (at frequency point IJ and depth ID) C SCAD - scattering coefficient (at frequency point IJ and depth ID) C Dxy - derivatives of x (=A for Absorption, =E for Emission, C =S for scattering) coefficient wrt y (=T for temperature, C =N for electron density) C DDN(I) - quantity proportional to a derivative of absorption (and C emission) coefficient wrt population of I-th level; ie. C d(abs)/d(pop) = DDN * [1-exp(-h*nu/kT)] C d(em)/d(pop) = DDN * (2h*nu**3/c**2)*exp(-h*nu/kT) C C INCLUDE 'IMPLIC.FOR' INCLUDE 'BASICS.FOR' INCLUDE 'ATOMIC.FOR' INCLUDE 'MODELQ.FOR' common/eospar/anmol(600,mdepth), * anato(100,mdepth), * anion(100,mdepth) c PARAMETER (FRRAY = 2.463D15, * FRAYHe= 5.150E15, * FRAYH2= 2.922E15, * CLS = 2.997925e18, * C18 = 2.997925D18, * CR0 = 5.799D-13, * CR1 = 1.422D-6, * CR2 = 2.784D0, * TENM4 = 1.0D-4, * THM0 = 8.7629D3, * SBHM = 1.0353D-16, * TRHA = 1.5D0, * SFF0 = 1.3727D-25, * SFF1 = 4.3748D-10, * SFFM2 = -2.5993D-7, * F0HE1 = 3.29D15, * F0HE2 = 1.316D16, * SBH0 = 4.1412D-16, * SG01 = 2.815D-16, * SG02 = 4.504D-15) SAVE T,DELTAT,ANE,HKT,T32,XHM,POPI,SB00 C AB0=0. AB1=0. DAB0=0. DAB1=0. ABAD=0. EMAD=0. SCAD=0. DBF=0. DFF=0. DAT=0. DET=0. DAN=0. DEN=0. DST=0. DSN=0. DO I=1,NLEVEL DDN(I)=0. END DO C FR=FREQ(IJ) al=2.997925e18/fr lpri=al.gt.1579.0.and.al.lt.1579.5 if(ielh.gt.0) then N0HN=NFIRST(IELH) NKH=NKA(IATH) end if c IF(ICALL.GT.0) THEN T=TEMP(ID) DELTAT=TENM4*T ANE=ELEC(ID) HKT=HK/T T32=UN/T/SQRT(T) XHM=THM0/T if(ielh.gt.0) then ah=popul(n0hn,id) ahp=popul(nkh,id) else ah=anato(1,id) ahp=anion(1,id) end if popi=ah SB00=SBHM*T32*EXP(XHM)*POPI*ANE c if(iathe.gt.0) then ahe=popul(n0a(iathe),id) else ahe=anato(2,id) end if END IF C IT=NCON C C ----------------------- C HI Rayleigh scattering C ----------------------- C IF(IRSCT.NE.0) THEN IT=IT+1 SCAD=AH*CROSS(IT,IJ) END IF C C ----------------------- C He I Rayleigh scattering C ----------------------- C IF(IRSCHE.NE.0.AND.MODE.GE.0) THEN IT=IT+1 scad=scad+ahe*cross(it,ij) END IF C C ----------------------- C H2 Rayleigh scattering C ----------------------- C IF(IRSCH2.NE.0.AND.MODE.GE.0.AND.IFMOL.GT.0) THEN IT=IT+1 sg=cross(it,ij)*anmol(2,id) if(t.lt.tmolim) scad=scad+sg END IF C IF(IOPHMI.GT.0) THEN C C ---------------------------- C H- bound-free and free-free C ---------------------------- C IT=IT+1 if(t.lt.20000.) then SB=SB00*CROSS(IT,IJ) SF=SFFHMI(AH,FR,T)*ANE AB0=SB+SF end if END IF c END IF C IF(IOPH2P.GT.0) THEN C C ----------------------------- C H2+ bound-free and free-free C ----------------------------- C IT=IT+1 X2=-CROSS(IT,IJ)/T+CROSS(IT+1,IJ) IT=IT+1 SB=0. IF(X2.GT.-150..and.fr.lt.3.28e15.and.t.le.9000.) * SB=AH*EXP(X2)*AHP AB0=AB0+SB DAB1=SB*CROSS(IT-1,IJ)/T/T IF(N0HN.GT.0) DDN(N0HN)=DDN(N0HN)+SB/POPI IF(NKH.GT.0) DDN(NKH)=DDN(NKH)+SB/POPUL(NKH,ID) END IF C C ----------------------------- C He- free-free C ----------------------------- C IF(IOPHEM.GT.0) THEN IT=IT+1 sg=cross(it,ij)*t+cross(it+1,ij)+cross(it+2,ij)/t ab0=ab0+sg*ane*ahe 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,anmol(2,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)*anmol(5,id) if(iopoh.gt.0) ab0=ab0+sbfoh(fr,t)*anmol(4,id) C C --------------------------- C CIA H2-H2 opacity C --------------------------- C if(ioh2h2.gt.0) then call cia_h2h2(t,anmol(2,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,anmol(2,id),ahe,fr,oph2) ab1=ab1+oph2 end if C C --------------------------- C CIA H2-H opacity C --------------------------- C if(ioh2h.gt.0) then call cia_h2h(t,anmol(2,id),ah,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,ah,ahe,fr,oph2) ab1=ab1+oph2 end if end if C C C ---------------------------------------------- C The user may supply more opacity sources here: C ---------------------------------------------- C C Finally, actual absorption and emission coefficients and C their derivatives IF(MODE.LT.0) RETURN X=EXP(-HKT*FR) X1=UN-X FR15=FR*1.D-15 BNX=BN*FR15*FR15*FR15*X AB1=AB1/X1 ABAD=ABAD+AB0+AB1 EMAD=EMAD+AB0+AB1 IF(MODE.EQ.1) THEN HKFT=HKT*FR/T DB=HKFT*AB0 DAT=DAB1 DET=DAB1 DAN=AB0/ANE DEN=AB0/ANE END IF RETURN END