269 lines
6.8 KiB
Fortran
269 lines
6.8 KiB
Fortran
SUBROUTINE LUCY
|
|
C ===============
|
|
C
|
|
C NLTE Lucy-Unsold temperature correction scheme, following
|
|
C Werner & Dreizler
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'ATOMIC.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
INCLUDE 'ODFPAR.FOR'
|
|
INCLUDE 'ITERAT.FOR'
|
|
INCLUDE 'ALIPAR.FOR'
|
|
INCLUDE 'ARRAY1.FOR'
|
|
DIMENSION HEAT(MDEPTH),ABSZ(MDEPTH),ABSP(MDEPTH),ABSH(MDEPTH),
|
|
* TOTJ(MDEPTH),TOTH(MDEPTH),TOTB(MDEPTH),DELH(MDEPTH),
|
|
* EDDF(MDEPTH),DELTAT(MDEPTH),dt1(mdepth),dt2(mdepth),
|
|
* ANTC(MDEPTH),XE(MDEPTH),heab(mdepth),tau(mdepth),
|
|
* TEM0(MDEPTH),TEM1(MDEPTH),TEM2(MDEPTH),TEM3(MDEPTH)
|
|
DIMENSION COL(MTRANS),CLOC(MTRANS)
|
|
PARAMETER (THIRD=1./3.)
|
|
DATA ILINIT /0/
|
|
C
|
|
c
|
|
IF(ILINIT.EQ.0) THEN
|
|
ILINIT=1
|
|
do itr=1,ntrans
|
|
iluctr(itr)=0
|
|
end do
|
|
ntrl=0
|
|
read(2,*,end=5,err=5) ntrl
|
|
5 continue
|
|
if(ntrl.gt.0) then
|
|
do i=1,ntrl
|
|
read(2,*,end=6,err=6) itrl
|
|
iluctr(itrl)=1
|
|
end do
|
|
end if
|
|
6 continue
|
|
end if
|
|
c
|
|
c IF(MOD(ILAM,4).NE.3) RETURN
|
|
if(itlucy.le.0) return
|
|
LAC2T=.FALSE.
|
|
iacc0t=IACLT-3
|
|
C
|
|
ilucy=1
|
|
10 continue
|
|
CALL OPAINI(0)
|
|
C
|
|
DO ID=1,ND
|
|
HEAT(ID)=0.
|
|
heab(id)=0.
|
|
ABSZ(ID)=0.
|
|
ABSP(ID)=0.
|
|
ABSH(ID)=0.
|
|
TOTJ(ID)=0.
|
|
TOTH(ID)=0.
|
|
TOTB(ID)=0.
|
|
EDDF(ID)=0.
|
|
END DO
|
|
EDDH=0.
|
|
C
|
|
DO IJ=1,NFREQ
|
|
W0=W(IJ)
|
|
FR=FREQ(IJ)
|
|
FR15=FR*1.D-15
|
|
BNU=BN*FR15*FR15*FR15
|
|
CALL OPACFL(IJ)
|
|
CALL RTEFR1(IJ)
|
|
tau(1)=abso1(1)/dens(1)*dm(1)
|
|
DO ID=1,ND
|
|
if(id.ge.2) tau(id)=tau(id-1)+half*(abso1(id)/dens(id)+
|
|
* abso1(id-1)/dens(id-1))*(dm(id)-dm(id-1))
|
|
PLAND=BNU/(EXP(HK*FR/TEMP(ID))-UN)
|
|
ABSOT0=ABSO1L(ID)-SCAT1(ID)
|
|
HEAT(ID)=HEAT(ID)+W0*(ABSOT0*RAD1(ID)-EMIS1L(ID))
|
|
HEAB(ID)=HEAB(ID)+W0*(ABSOT0*RAD1(ID)-EMIS1L(ID))/
|
|
* (ABSOT0*PLAND)
|
|
ABSP(ID)=ABSP(ID)+W0*ABSOT0*PLAND
|
|
TOTJ(ID)=TOTJ(ID)+W0*RAD1(ID)
|
|
TOTB(ID)=TOTB(ID)+W0*PLAND
|
|
EDDF(ID)=EDDF(ID)+W0*FAK1(ID)*RAD1(ID)
|
|
END DO
|
|
C
|
|
ID=1
|
|
EDDH=EDDH+W0*RAD1(ID)*FH(IJ)
|
|
TOTH(ID)=TOTH(ID)+W0*(RAD1(ID)*FH(IJ)-HEXTRD(IJ))
|
|
ABSOT0=ABSO1(ID)
|
|
ABSH(ID)=ABSH(ID)+W0*ABSOT0*(RAD1(ID)*FH(IJ)-HEXTRD(IJ))
|
|
DO ID=2,ND
|
|
ABSOT0=HALF*(ABSO1(ID)+ABSO1(ID-1))
|
|
DTM=DELDM(ID-1)*(ABSO1(ID)*DENS1(ID)+ABSO1(ID-1)*DENS1(ID-1))
|
|
FLUZ=(RAD1(ID)*FAK1(ID)-RAD1(ID-1)*FAK1(ID-1))/DTM
|
|
TOTH(ID)=TOTH(ID)+W0*FLUZ
|
|
ABSH(ID)=ABSH(ID)+W0*ABSOT0*FLUZ
|
|
END DO
|
|
END DO
|
|
C
|
|
EDDH=EDDH/TOTJ(1)
|
|
TEF4=SIG4P*TEFF**4
|
|
TOTH(ND)=TEF4
|
|
HEAT(ND)=0.
|
|
HEAT(ND-1)=0.
|
|
DO ID=1,ND
|
|
ABSZ(ID)=(HEAT(ID)+ABSP(ID))/TOTJ(ID)
|
|
ABSP(ID)=ABSP(ID)/TOTB(ID)
|
|
ABSH(ID)=ABSH(ID)/TOTH(ID)
|
|
EDDF(ID)=EDDF(ID)/TOTJ(ID)
|
|
DELH(ID)=-TOTH(ID)+TEF4
|
|
END DO
|
|
DHHMX1=0.
|
|
DO ID=1,ND-1
|
|
DHH=ABS(DELH(ID)/TOTH(ID))
|
|
IF(DHH.GT.DHHMX1) DHHMX1=DHH
|
|
END DO
|
|
C
|
|
ID=1
|
|
TP3=TEMP1(ID)*TEMP1(ID)*TEMP1(ID)
|
|
XX=EDDF(ID)/EDDH*DELH(ID)
|
|
XX1=XX
|
|
DT1(id)=HEAT(ID)/16./SIG4P*TP3/ABSP(ID)
|
|
DT2(id)=ABSZ(ID)/EDDF(ID)*XX/16./SIG4P*TP3/ABSP(ID)
|
|
DELTAT(ID)=dt1(id)+dt2(id)
|
|
DO ID=2,ND
|
|
TP3=TEMP1(ID)*TEMP1(ID)*TEMP1(ID)
|
|
XX=XX+DELDM(ID)*(ABSH(ID)*DENS1(ID)*DELH(ID)+
|
|
* ABSH(ID+1)*DENS1(ID+1)*DELH(ID+1))
|
|
DT1(id)=HEAT(ID)/16./SIG4P*TP3/ABSP(ID)
|
|
DT2(id)=ABSZ(ID)/EDDF(ID)*XX/16./SIG4P*TP3/ABSP(ID)
|
|
DELTAT(ID)=dt1(id)+dt2(id)
|
|
END DO
|
|
C
|
|
C New temperature
|
|
C
|
|
DO ID=1,ND
|
|
TEMP(ID)=TEMP(ID)+DELTAT(ID)
|
|
tem0(id)=temp(id)
|
|
AOLD=DENS(ID)/WMM(ID)+ELEC(ID)
|
|
XE(ID)=UN-ELEC(ID)/AOLD
|
|
END DO
|
|
c
|
|
c acceleration
|
|
c
|
|
if(itlucy.lt.IACLT .or. ilucy.lt.iacc0t) go to 20
|
|
ipng=1
|
|
if(IACLDT.gt.0) ipng=mod((ilucy-IACLT),IACLDT)
|
|
if(.not.lac2t) then
|
|
IPT=MOD(ILUCY,3)
|
|
IPT0=MOD(IACLT,3)
|
|
IPT1=MOD((IACLT+1),3)
|
|
IPT2=MOD((IACLT+2),3)
|
|
IF(ILUCY.EQ.IACC0T) THEN
|
|
DO ID=1,ND
|
|
TEM3(ID)=TEM0(ID)
|
|
END DO
|
|
ELSE IF(IPT.EQ.IPT1) THEN
|
|
DO ID=1,ND
|
|
TEM2(ID)=TEM0(ID)
|
|
END DO
|
|
ELSE IF(IPT.EQ.IPT2) THEN
|
|
DO ID=1,ND
|
|
TEM1(ID)=TEM0(ID)
|
|
END DO
|
|
ENDIF
|
|
else if (ipng.ne.0) then
|
|
DO ID=1,ND
|
|
TEM3(ID)=TEM2(ID)
|
|
END DO
|
|
DO ID=1,ND
|
|
TEM2(ID)=TEM1(ID)
|
|
END DO
|
|
DO ID=1,ND
|
|
TEM1(ID)=TEM0(ID)
|
|
END DO
|
|
GO TO 20
|
|
end if
|
|
|
|
IF(ILUCY.LT.IACLT) go to 20
|
|
C
|
|
A1=0.
|
|
B1=0.
|
|
B2=0.
|
|
C1=0.
|
|
C2=0.
|
|
DO ID=1,ND
|
|
WT=0.
|
|
IF(TEM0(ID).NE.0.) WT=1./ABS(TEM0(ID))
|
|
D0=TEM0(ID)-TEM1(ID)
|
|
D1=D0-TEM1(ID)+TEM2(ID)
|
|
D2=D0-TEM2(ID)+TEM3(ID)
|
|
A1=A1+WT*D1*D1
|
|
B1=B1+WT*D1*D2
|
|
B2=B2+WT*D2*D2
|
|
C1=C1+WT*D0*D1
|
|
C2=C2+WT*D0*D2
|
|
END DO
|
|
AB=B2*A1-B1*B1
|
|
IF(AB.EQ.0.) THEN
|
|
WRITE(6,601) ILUCY,AB
|
|
IACLT=IACLT+IACLDT
|
|
IACC0T=IACLT-3
|
|
GO TO 20
|
|
ENDIF
|
|
A0=(B2*C1-B1*C2)/AB
|
|
B0=(A1*C2-B1*C1)/AB
|
|
DO ID=1,ND
|
|
TEM0(ID)=(1.-A0-B0)*TEM0(ID)+A0*TEM1(ID)+
|
|
* B0*TEM2(ID)
|
|
TEMP(ID)=TEM0(ID)
|
|
END DO
|
|
LAC2T=.TRUE.
|
|
601 FORMAT(/,' **** ACCELT, ITER=',I4,' AB = ',F7.3,/)
|
|
20 CONTINUE
|
|
|
|
CALL TDPINI
|
|
C
|
|
C Integrate hydrostatic equilibrium
|
|
C
|
|
IF(IHECOR.GE.1) THEN
|
|
ID=1
|
|
PTUR=HALF*VTURB(ID)*VTURB(ID)*DENS(ID)
|
|
ANTC(ID)=(DM(ID)*GRAV-PRD0-PTUR)/BOLK/TEMP(ID)
|
|
IF(ANTC(ID).LE.0) ANTC(ID)=DENS(ID)/WMM(ID)+ELEC(ID)
|
|
DO ID=2,ND
|
|
PTUR=HALF*VTURB(ID)*VTURB(ID)*DENS(ID)
|
|
PTURM=HALF*VTURB(ID-1)*VTURB(ID-1)*DENS(ID-1)
|
|
ANTC(ID)=(GRAV*(DM(ID)-DM(ID-1))+
|
|
* BOLK*TEMP(ID-1)*ANTC(ID-1)-
|
|
* PRADT(ID)+PRADT(ID-1)-PTUR+PTURM)/
|
|
* BOLK/TEMP(ID)
|
|
END DO
|
|
DO ID=1,ND
|
|
ELEC(ID)=(UN-XE(ID))*ANTC(ID)
|
|
DENS(ID)=WMM(ID)*(ANTC(ID)-ELEC(ID))
|
|
dens1(id)=un/dens(id)
|
|
END DO
|
|
END IF
|
|
C
|
|
C Other depth-dependent quantities
|
|
C
|
|
DO ID=1,ND
|
|
PGS(ID)=(DENS(ID)/WMM(ID)+ELEC(ID))*BOLK*TEMP(ID)
|
|
CALL WNSTOR(ID)
|
|
IF(LTE) GO TO 60
|
|
CALL SABOLF(ID)
|
|
CALL COLIS(ID,TEMP(ID),COL,CLOC)
|
|
DO I=1,NTRANS
|
|
COLRAT(I,ID)=COL(I)
|
|
COLTAR(I,ID)=CLOC(I)
|
|
END DO
|
|
60 CONTINUE
|
|
END DO
|
|
C
|
|
C new populations
|
|
C
|
|
DO ID=1,ND
|
|
CALL STEQEQ(ID,POP,1)
|
|
IF(.NOT.LCHC.and.iter.lt.ielcor) CALL ELCOR(ID)
|
|
END DO
|
|
C
|
|
CALL CONCOR
|
|
CALL ODFMER
|
|
ilucy=ilucy+1
|
|
if(ilucy.le.itlucy) go to 10
|
|
C
|
|
RETURN
|
|
END
|