91 lines
2.8 KiB
Fortran
91 lines
2.8 KiB
Fortran
SUBROUTINE TAUFR1(IJ)
|
|
C =====================
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
INCLUDE 'ALIPAR.FOR'
|
|
INCLUDE 'ITERAT.FOR'
|
|
COMMON/OPTDPT/DT(MDEPTH)
|
|
dimension ST0(MDEPTH),SS0(MDEPTH),ab0(mdepth),
|
|
* tau(mdepth),taus(mdepth)
|
|
c PARAMETER (TAUREF = 0.6666666666667)
|
|
PARAMETER (TAUREF = 1.)
|
|
PARAMETER (XCON=8.0935D-21,YCON=1.68638E-10)
|
|
C
|
|
FR=FREQ(IJ)
|
|
DO ID=1,ND
|
|
AB0(ID)=ABSO1(ID)
|
|
C put in a floor to avoid division by zero:
|
|
IF(AB0(ID)-SCAT1(ID).GT.1.D-100) THEN
|
|
ST0(ID)=EMIS1(ID)/(AB0(ID)-scat1(id))
|
|
ELSE
|
|
ST0(ID)=0.
|
|
ENDIF
|
|
if(st0(id).eq.0) st0(id)=1.d-20*scat1(id)
|
|
SS0(ID)=-SCAT1(ID)/AB0(ID)
|
|
END DO
|
|
C
|
|
id=1
|
|
TAUMIN=ABSO1(1)*DEDM1
|
|
tau(1)=taumin
|
|
C to avoid a negative square root:
|
|
taus(1)=sqrt(3.*ab0(id)*max(ab0(id)-scat1(id),0.d0))*DEDM1
|
|
C
|
|
IREF=1
|
|
IREFs=1
|
|
DO ID=1,ND-1
|
|
DT(ID)=DELDMZ(ID)*(ABSOT(ID+1)+ABSOT(ID))
|
|
tau(id+1)=tau(id)+dt(id)
|
|
C to avoid negative square root:
|
|
eps0=sqrt(ab0(id)*max(ab0(id)-scat1(id),0.d0))
|
|
eps1=sqrt(ab0(id+1)*max(ab0(id+1)-scat1(id+1),0.d0))
|
|
dts=deldm(id)*(eps0*dens1(id)+eps1*dens1(id+1))*sqrt(3.)
|
|
taus(id+1)=taus(id)+dts
|
|
IF(TAU(Id).LE.TAUREF.AND.TAU(Id+1).GT.TAUREF) IREF=Id
|
|
IF(TAUs(Id).LE.TAUREF.AND.TAUs(Id+1).GT.TAUREF) IREFs=Id
|
|
END DO
|
|
if(iref.eq.1.and.tau(nd).le.tauref) iref=nd
|
|
if(irefs.eq.1.and.taus(nd).le.tauref) irefs=nd
|
|
C
|
|
t0=1.
|
|
iref0=iref
|
|
iref=irefs
|
|
if(irefs.lt.nd) then
|
|
T0=LOG(TAUs(IREF+1)/TAUs(IREF))
|
|
X0=LOG(TAUs(IREF+1)/TAUREF)/T0
|
|
X1=LOG(TAUREF/TAUs(IREF))/T0
|
|
DMREF=EXP(LOG(DM(IREF))*X0+LOG(DM(IREF+1))*X1)
|
|
TREF=EXP(LOG(TEMP(IREF))*X0+LOG(TEMP(IREF+1))*X1)
|
|
abREF=EXP(LOG(ab0(IREF))*X0+LOG(ab0(IREF+1))*X1)
|
|
scREF=EXP(LOG(scat1(IREF))*X0+LOG(scat1(IREF+1))*X1)
|
|
STREF=EXP(LOG(ST0(IREF))*X0+LOG(ST0(IREF+1))*X1)
|
|
tauef=EXP(LOG(TAU(IREF))*X0+LOG(TAU(IREF+1))*X1)
|
|
else
|
|
x0=1.
|
|
x1=0.
|
|
dmref=dm(nd)
|
|
tref=temp(nd)
|
|
abref=ab0(nd)
|
|
scref=scat1(nd)
|
|
stref=st0(nd)
|
|
tauef=tau(nd)
|
|
end if
|
|
epref=(abref-scref)/abref
|
|
CX add if statement to avoid overflow:
|
|
IF(hk*fr/tref.lt.200.) then
|
|
bref=1.4743e-2*(fr*1.e-15)**3/(exp(hk*fr/tref)-1.)
|
|
ELSE
|
|
bref=0.
|
|
END IF
|
|
taur=tauef*tauef
|
|
if(tauef.gt.taur) taur=tauef
|
|
yref=4.*ycon*tref*taur
|
|
ALM=2.997925E18/FREQ(IJ)
|
|
r1=rad(nfreq-kij(ij)+1,1)
|
|
if(epref.ge.0) rb1=sqrt(epref)*bref
|
|
if(epref.ge.0) rs1=sqrt(epref)*stref
|
|
C
|
|
return
|
|
end
|