116 lines
3.7 KiB
Fortran
116 lines
3.7 KiB
Fortran
SUBROUTINE NEWDMT
|
|
C =================
|
|
C
|
|
C New m-scale, calculated as that corresponding to the new
|
|
C grid better representing temperature variations
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
DIMENSION DM0(MDEPTH),DM11(MDEPTH),DENS0(MDEPTH),ZD0(MDEPTH),
|
|
* T0(MDEPTH),T1(MDEPTH),ELEC0(MDEPTH),PT0(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
|
|
DO ID=1,ND
|
|
DM0(ID)=LOG10(DM(ID))
|
|
T0(ID)=LOG10(TAUROS(ID))
|
|
ELEC0(ID)=ELEC(ID)
|
|
DENS0(ID)=DENS(ID)
|
|
PT0(ID)=PTOTAL(ID)
|
|
ZD0(ID)=ZD(ID)
|
|
ABRS0(ID)=ABROSD(ID)
|
|
ABPL0(ID)=ABPLAD(ID)
|
|
END DO
|
|
ND1=ND-1
|
|
CALL GRIDP(DM0,T0,DM11,T1,ND1)
|
|
DM11(ND)=DM0(ND)
|
|
T1(ND)=T0(ND)
|
|
DO ID=1,ND
|
|
DM(ID)=EXP(2.3025851*DM11(ID))
|
|
TAUROS(ID)=EXP(2.3025851*T1(ID))
|
|
END DO
|
|
CALL INTERP(DM0,ELEC0,DM11,ELEC,ND,ND,2,0,1)
|
|
CALL INTERP(DM0,DENS0,DM11,DENS,ND,ND,2,0,1)
|
|
CALL INTERP(DM0,PT0,DM11,PTOTAL,ND,ND,2,0,1)
|
|
CALL INTERP(DM0,ZD0,DM11,ZD,ND,ND,2,0,0)
|
|
CALL INTERP(DM0,ABRS0,DM11,ABROSD,ND,ND,2,0,1)
|
|
CALL INTERP(DM0,ABPL0,DM11,ABPLAD,ND,ND,2,0,1)
|
|
DO ID=1,ND
|
|
VSND2(ID)=PTOTAL(ID)/DENS(ID)
|
|
END DO
|
|
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
|