206 lines
6.1 KiB
Fortran
206 lines
6.1 KiB
Fortran
SUBROUTINE NEWDM
|
|
C ================
|
|
C
|
|
C New m-scale, calculated as that corresponding to a new
|
|
C tau(Ross)-scale, which is logarithmically equidistant, with
|
|
C a generally different step in six different regions:
|
|
C
|
|
C 1. region between the original tau(1) and T0 (taken as 0.01);
|
|
C 2. region between T0 and TC0 (taken as 0.01 and 0.1) - denser
|
|
C mesh (with number of points N0)
|
|
C 3. region between TC0 and TC1 (taken as 0.1 and 10.) - the
|
|
C central region with densest mesh (N1 points)
|
|
C 4. region between TC1 and T1 (taken as 10. and 100.) - as dens
|
|
C as the second region (N0 points)
|
|
C 5. the remaining region between T1 and the original last tau
|
|
C (more precisely, the last-but-one point TAUROS(ND-1).
|
|
C
|
|
C If T1 is greater than the original last tau, the new tau-scale
|
|
C is equidistant between T0 and the last tau.
|
|
C
|
|
C The procedure also calulates all the necessary state parameters
|
|
C for the new depth scale (density, z, pressure, opacities, and
|
|
C temperature)
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
PARAMETER (TEN=1.D1)
|
|
DIMENSION TAU(MDEPTH),TAUL(MDEPTH),DM0(MDEPTH),DENS0(MDEPTH),
|
|
* ABRS0(MDEPTH),ABPL0(MDEPTH)
|
|
COMMON/PRSAUX/VSND2(MDEPTH),HG1,HR1,RR1
|
|
COMMON/FACTRS/GAMJ(MDEPTH),GAMH,FAK0
|
|
COMMON/FLXAUX/T4,PGAS,PRAD,PGM,PRADM,ITGMAX,ITGMX0
|
|
C
|
|
DATA N0,NC0 /8,24/
|
|
DATA T0,TC0,TC1,T1 /-2.D0,-1.D0,1.D0,2.D0/
|
|
C
|
|
DO ID=1,ND
|
|
DM0(ID)=DM(ID)
|
|
DENS0(ID)=DENS(ID)
|
|
ABRS0(ID)=ABROSD(ID)
|
|
ABPL0(ID)=ABPLAD(ID)
|
|
TAUL(ID)=LOG10(TAUROS(ID))
|
|
IF(TAUL(ID).LT.T0) IMIN=ID
|
|
IF(TAUL(ID).LT.T1) IMAX=ID
|
|
END DO
|
|
ND1=ND-1
|
|
NC=2*N0+NC0
|
|
NB=ND1-NC
|
|
IF(IMAX.GE.ND1) THEN
|
|
IC=0
|
|
ELSE
|
|
X=(TAUL(IMIN)-TAUL(1))/(TAUL(ND1)-TAUL(IMAX))
|
|
X1=FLOAT(NB)/(X+UN)
|
|
IC=int(X1)
|
|
END IF
|
|
NB0=NB-IC
|
|
C
|
|
C New tau-scale
|
|
C
|
|
C First, logarithmically equidistant tau-points between tau(1)
|
|
C (which is the previous TAUROS(1)), and log tau = T0.
|
|
C Their number is NB0.
|
|
c
|
|
DT=(T0-TAUL(1))/FLOAT(NB0-1)
|
|
TAU(1)=TAUL(1)
|
|
DO ID=2,NB0
|
|
TAU(ID)=TAU(ID-1)+DT
|
|
END DO
|
|
C
|
|
IF(IC.GT.0) THEN
|
|
C
|
|
C 2.region - between log tau = T0 and TC0
|
|
C
|
|
DT=(TC0-T0)/FLOAT(N0)
|
|
DO I=1,N0
|
|
TAU(NB0+I)=TAU(NB0+I-1)+DT
|
|
END DO
|
|
C
|
|
C 3. The most dense region between TC0 and TC1 (central part)
|
|
C
|
|
NB1=NB0+N0
|
|
DT=(TC1-TC0)/FLOAT(NC0)
|
|
DO I=1,NC0
|
|
TAU(NB1+I)=TAU(NB1+I-1)+DT
|
|
END DO
|
|
C
|
|
C 4. The part similar to that between T0 and TC0, this time
|
|
C betwen TC1 and T1
|
|
C
|
|
NB2=NB1+NC0
|
|
DT=(T1-TC1)/FLOAT(N0)
|
|
DO I=1,N0
|
|
TAU(NB2+I)=TAU(NB2+I-1)+DT
|
|
END DO
|
|
C
|
|
C 5. The remainig part between T1 and the last-but-one tau
|
|
C
|
|
NB3=NB2+N0
|
|
DT=(TAUL(ND1)-T1)/FLOAT(IC)
|
|
DO I=1,IC
|
|
TAU(NB3+I)=TAU(NB3+I-1)+DT
|
|
END DO
|
|
TAU(ND)=TAUL(ND)
|
|
C
|
|
C The case where the last tau is smaller than T1; in this case
|
|
C the points are logarithmically equidistant between T0 and the
|
|
C last-but-one tau
|
|
C
|
|
ELSE
|
|
DT=(TAUL(ND1)-T0)/FLOAT(NC)
|
|
DO I=1,NC
|
|
TAU(NB0+I)=TAU(NB0+I-1)+DT
|
|
END DO
|
|
TAU(ND)=TAUL(ND)
|
|
C
|
|
END IF
|
|
C
|
|
C ---------------------------------------
|
|
C
|
|
C Final new Rosseland optical depth scale
|
|
C
|
|
DO ID=1,ND
|
|
TAU(ID)=TEN**TAU(ID)
|
|
END DO
|
|
C
|
|
C Interpolation from the old to the new tau(Ross) scale to
|
|
C get the new m-scale, density and Planck mean opacity
|
|
C
|
|
CALL INTERP(TAUROS,DM0,TAU,DM,ND,ND,2,1,1)
|
|
CALL INTERP(DM0,DENS0,DM,DENS,ND,ND,2,1,1)
|
|
CALL INTERP(DM0,ABRS0,DM,ABROSD,ND,ND,2,1,1)
|
|
CALL INTERP(DM0,ABPL0,DM,ABPLAD,ND,ND,2,1,1)
|
|
C
|
|
C New Rosseland opacity and functions theta and tauthe
|
|
C
|
|
AMUV0=DMVISC**(ZETA0+UN)
|
|
AMUV1=UN-AMUV0
|
|
DO ID=1,ND
|
|
IF(DM(ID).LE.DMVISC*DM(ND)) THEN
|
|
VISCD(ID)=(UN-FRACTV)*(ZETA1+UN)/
|
|
* DMVISC**(ZETA1+UN)*(DM(ID)/DM(ND))**ZETA1
|
|
THETA(ID)=(UN-FRACTV)*(DM(ID)/DMVISC/DM(ND))**(ZETA1+UN)
|
|
ELSE
|
|
VISCD(ID)=FRACTV*(ZETA0+UN)/AMUV1*
|
|
* (DM(ID)/DM(ND))**ZETA0
|
|
THETA(ID)=(UN-FRACTV)+FRACTV*((DM(ID)/DM(ND))**(ZETA0+UN)-
|
|
* AMUV0)/AMUV1
|
|
END IF
|
|
GAMJ(ID)=UN
|
|
IF(ID.EQ.1) THEN
|
|
TAUROS(ID)=DM(ID)*ABROSD(ID)
|
|
TAUTHE(ID)=TAUROS(ID)*THETA(ID)/(ZETA1+TWO)
|
|
ANEREL=ELEC(ID)/(DENS(ID)/WMM(ID)+ELEC(ID))
|
|
ELSE
|
|
DDM=DM(ID)-DM(ID-1)
|
|
TAUROS(ID)=TAUROS(ID-1)+DDM*HALF*(ABROSD(ID-1)+ABROSD(ID))
|
|
ZETAD=ZETA0
|
|
IF(DM(ID).LE.DMVISC*DM(ND)) ZETAD=ZETA1
|
|
A0=(ABROSD(ID-1)*DM(ID)-ABROSD(ID)*DM(ID-1))/DDM/
|
|
* (ZETAD+TWO)
|
|
A1=(ABROSD(ID)-ABROSD(ID-1))/DDM/(ZETAD+3.D0)
|
|
TAUTHE(ID)=TAUTHE(ID-1)+
|
|
* A0*(THETA(ID)*DM(ID)-THETA(ID-1)*DM(ID-1))+
|
|
* A1*(THETA(ID)*DM(ID)**2-THETA(ID-1)*DM(ID-1)**2)
|
|
END IF
|
|
TAUR=TAUROS(ID)
|
|
CALL TEMPER(ID,TAUR,1)
|
|
END DO
|
|
C
|
|
C Next step - simultaneous solution of the hydrostatic
|
|
C equilibrium and the z-m relation
|
|
C
|
|
if(nconit.ge.0) CALL HESOLV
|
|
C
|
|
C New temperature and mean opacities for the current density
|
|
C and pressure
|
|
C
|
|
DO ID=1,ND
|
|
TAUR=TAUROS(ID)
|
|
CALL TEMPER(ID,TAUR,1)
|
|
END DO
|
|
C
|
|
C Once again - simultaneous solution of the hydrostatic
|
|
C equilibrium and the z-m relation
|
|
C
|
|
if(nconit.ge.0) CALL HESOLV
|
|
C
|
|
IF(IPRING.GE.1) THEN
|
|
WRITE(6,601)
|
|
DO ID=1,ND
|
|
WRITE(6,602) ID,DM(ID),TAUROS(ID),
|
|
* TEMP(ID),ELEC(ID),PTOTAL(ID),
|
|
* ZD(ID),ABROSD(ID),ABPLAD(ID)
|
|
END DO
|
|
END IF
|
|
C
|
|
601 FORMAT(1H1,' NEW DEPTH GRID ESTABLISHED, NEW MODEL:'/
|
|
* ' --------------------------------------'/
|
|
* ' ID DM TAUROSS TEMP NE P',
|
|
* 8X,'ZD ROSS.MEAN PLANCK'/)
|
|
602 FORMAT(1H ,I3,1P2D9.2,0PF8.0,1P3D9.2,2X,2D9.2)
|
|
C
|
|
RETURN
|
|
END
|