-- version for calling ronfict directly
_DEFINE
@SAVE_MODEL SAVEM
@TAUPH 0.64D0
@HOLD_BAL HOLDBL -- OLD STEP FOR BALANCE
@ReadData "READ(Nc,'(A)'); Read(Nc,*)"
@U (CCL/3.65D+03)
@B (CCL/4.4D+03)
@V (CCL/5.5D+03)
@R (CCL/7.D+03)
@I (CCL/9.D+03)
@Lyman (CCL/912.D0) ;
_LABEL @F1,@F2;
SUBROUTINE BEGIN;
use mLbalsw, only: lossen; --, pribal, primag;
-- _TRACE "@wterm' Nzon=',Nzon,' Mzon=',Mzon,"
-- _TRACE "@wterm' Ncnd=',Ncnd,' Krad=',Krad,"
-- _TRACE "@wterm' Natom=',Natom,' Zn(Natom)=',Zn(Natom),"
IMPLICIT REAL*8(A-H,O-Z);
_INCLUDE snrad;
_Include opacity;
_INCLUDE abo;
-- real*8 hptrbab(Nfreq,@Ntab,@Ntab,Mzon/@skip),
-- hptrbsc(Nfreq,@Ntab,@Ntab,Mzon/@skip);
<*B: vars, PARAMETERS, data *>;
<*R: READ INPUT DATA FOR BURST *>;
RETURN;
9:WRITE(@Wres,*)' BEGIN: error in read file ',Lunit;
stop;
END;
<*S: SAVE CURRENT MODEL *>;
<*I: BLOCK DATA *>;
%B:O
CHARACTER*80 STRING;
-- Dimension dumFreq(Nfreq+1),dumFreqmn(Nfreq),dumwavel(Nfreq);
-- now in opacity.inc
Dimension indfr(6), indop(6); -- keys for freq and opacit
-- depending on Knadap
-- PARAMETER(UFRC=2.417965D+14/1.16048D+04);
CHARACTER*3 status,NFILE*80;
Logical petread; -- is Peters's file already read?
Logical fexist;
Data petread/.false./;
Data indfr/1,2,2,3,3,3/, -- 1 geom., 2 peter, 3 ronfict
indop/1,1,2,3,3,3/; -- 1 opazr, 2 opahoef,hap2int, 3 hapintron
-- DIMENSION WAVEL(NFREQ);-- WAVELENGTHS IN ANGSTREMS
-- Now WAVEL in Include Opacity
-- DATA WLMAX/5.0D+04/,WLMIN/1.D+01/;
DATA WLMAX/5.0D+04/,WLMIN/1.D+00/;
-- wavelength range in Angstrems for SN1987A
-- DATA WLMAX/5.0D+04/,WLMIN/9.D+01/; -- for SN2L
%R:
<*test* READ START_PACK DATA FOR RUN , CHANNEL - 8 *>;
<*D: READ START_PACK DATA FOR RUN , CHANNEL - 8 *>;
<*M: READ MODEL CHANNEL - 9 -- begin run
CHANNEL -10 -- continue run
READ Lcurve CHANNEL -11 *>;
<*U: UNITS & CONSTANTS FOR DIMENSIONLESS EQUATIONS *>;
inquire(file='lines1.out',exist=fexist);
if(fexist)then;
open(file='lines1.out ', unit=7, status='old');
else;
open(file='lines1.out ', unit=7, status='new');
endif;
If (ABS(Knadap) .ne. 1) then; -- see begrad for Knadap
<*haptab: read tables of happa if needed *>;
endif;
AS=ACARB; ZS=ZCARB;
N=NZON*NVARS+2*KRAD;
HMIN=HMIN/UTIME; HMAX=HMAX/UTIME;
METH=METHN;
Haplim=1.D0/(3.D0*FitTau);
Uplim=1.D0+Haplim;
<*H: PRINT HEADER OF THE RUN *>;
-- For enhanced quadratic viscosity:
-- Initially DRT in STRAD DATA defines the mass "distance"
-- in Solar mass units
-- from AM(NCND) with enhanced pseudoviscosity determined by BQ.
-- Then it is used in DFRAD & LBALRAD in the following form:
--
-- DRT=(UM/DRT)**2; -- initial DRT is destroyed here, be careful!!
-- for R-T do not change DRT !
--
-- For linear R-T viscosity (Hot or Cold)
-- DRT was the optical thickness (assuming Hydrogen & full
-- ionization where QRT begins to act)
-- later the value 30. was fixed & now
-- DRT determines the enhancement of EpsUq
-- NRT may be noninteger in RTphi - see STRAD
IF(NSTEP==0)THEN;
H=1.D+05*HMIN;
@HOLD_BAL=H; -- INITIAL CONDITIONS FOR BALANCES
YENTOT(1)=ELOST;YENTOT(2)=H*ELTOT;
<*F: PREPARE FREQUENCY ARRAYS etc. for Ronfict *>;
NFRUS=NFREQ;
_Do L=1,NFRUS;
NTHICK(L)=ncnd;
_od;
_Do KM=NCND+1,NZON;
_Do L=1,NFRUS;
LTHICK(L,KM)=.FALSE.;
_od;
_od;
ENDIF;
<*Freqob: PREPARE FREQUENCY ARRAY freqob from freqmn *>;
toldg=-1.d0;
WRITE(@Wres, '('' FREQ:'',1P,10E12.5)')FREQ;
WRITE(@Wres, '('' FREQMN:'',1P,10E12.5)')FREQMN;
WRITE(@Wres, '('' WAVE BOUNDS:'',1P,11E10.3)')
(CCL/FREQ(LW),LW=1,NFREQ+1);
WRITE(@Wres, '('' WAVES:'',1P,10E12.5)')(CCL/FREQMN(LW),LW=1,NFREQ);
WRITE(@Wres, '('' WEIGHT:'',1P,10E12.5)')WEIGHT;
<*BANDS: find numbers of freq. groups for UBV *>;
LST=2;KENTR=0;
-- RGASA=-.831434D0*(1.D0-ACARB/ASI);
-- IBURNT=0; -- INDICATES PRINTING IN BURNT ZONE
RHO(1)=3.D0*DM(1)/(Ry(1)**3-Rce**3);
_DO I=2,NZON;
-- write(@term,'(a,i3,2(a,1p,e12.4))')' i=',i,' dm=',dm(i),
-- ' ry=',Ry(i);
RHO(I)=3.D0*DM(I)/(Ry(I)**3-Ry(I-1)**3)
_OD;
IF(NSTEP==0)THEN;
CALL LOSSEN;
@HOLD_BAL=H; -- INITIAL CONDITIONS FOR BALANCES
YENTOT(1)=ELOST;YENTOT(2)=H*ELTOT;
ENDIF;
chem=0.5d0; -- to avoid too low Xe in saha
<*testhappa* *>;
%R_testhappa:
pl=1.d-8;
K=48;
_do itp=1,1;
Tp=1.d+01**(4d0+(itp-1))/UTp;
-- write(@term,'(a,i5,1p2e10.2)') ' Test happal for zone, pl*Urho, Tp*UTp',
-- K,pl*Urho,Tp*UTp;
-- _Do i=1, Natom; Yat(i)=YABUN(i,K)_od;
_Do i=1, Natom;
Yat(i)=0.d0;
If(i==14)Yat(i)=1.d0/AZ(14);
_od;
write(@term,'(a,1p3e10.2)') ' Test happal for Fe, pl*Urho, Tp*UTp',
Yat(14)*56.,pl*Urho,Tp*UTp;
RADP=.FALSE.;
CALL URSOS;
kmhap=K;
call Happa;
write(@term,'(a)') ' Test table happal:';
write(@term,'(1p10e10.2)') happal;
call HappaDirect;
write(@term,'(a)') ' Test direct happal:';
write(@term,'(1p10e10.2)') happal;
_od
%R_test:
_Repeat
READ(8,'(A)',End=110,Iostat=kio) STRING;
WRITE(@Wres,'(A)') string;
_until kio^=0;
110: stop;
%RD:
Nc=8;
@ReadData EPS,HMIN,HMAX;-- STIFF ACC.,MIN.STEP(SEC),MAX.STEP(SEC)
@wterm ' 1';
@ReadData METHN,JSTART,MAXORD,KNadap;
@wterm ' 1';
-- ADAMS(1)-GEAR(2)-BGH(3), 0-START , max order,
-- KNadap<0 for Ni from a foreign model
-- abs(KNadap): 1 freqmn - geom. progression direcr opacity
-- 2 freqmn - Peter Hoeflich, opazr
-- 3 freqmn - Peter Hoeflich, opahoef, hap2int
-- 4 freqmn - Ronfict, opahoef, hapintron
-- 5 freqmn - as 4, but homolog. vel.
-- 6 freqmn - as 4, but istim==0 always => hapron0
@ReadData NSTA,NSTB,TcurA,TcurB;-- Steps & time in days
@wterm ' 1';
-- For read:
-- if NSTA<0 then files 10 & 11, else files 12 & 13
-- used by STRAD for test runs, beginning
-- with first saved model having NSTEP>=abs(NSTA)
-- For Write:
-- if NSTB<0 then files 10 & 11, else files 12 & 13
-- To begin test runs with zero step (file 9) use
-- NSTA<0, Begrun==T, & NSTB==1
-- To begin test runs with some step in file 10 use
-- NSTA<0, Begrun==T, NSTB>0 & NSTB^=1 & NSTB>abs(NSTA)
-- To begin test runs with last step in file 10 use
-- NSTA<0, Begrun==T, NSTB>0 & NSTB^=1 & abs(NSTB)<abs(NSTA)
-- To continue test runs in files 12 & 13 use
-- NSTA>0, Begrun==F, NSTB>0 & NSTB<abs(NSTA), e.g. NSTB=1
-- For STINFO: ????
-- use always Begrun=T & NSTB^=1 (if NSTB==1 then initial model)
-- ????
@ReadData NSTMAX,NDebug,NOUT,IOUT,MBATCH;
-- Max STEP number, Debug step number, INTERVAL OUTPUT, LINE PRINT
@wterm ' 1';
Mbatch=Min(Mbatch,Lcurdm);
-- MASS OF NI CORE, NI 56 ABUNDANCE
@ReadData AMNI,XMNI; -- contaminated and true nickel initial mass
@wterm ' 1';
@ReadData AMHT,EBurst,tBurst,tbeght;-- Mass of Heated Core, Energy & tm
@wterm ' 1';
@ReadData EKOB,AI1,AI2,AI3,US; -- KINETIC ENERGY(E+50)
@wterm ' 1';
-- Mass fractions TRIANGLE VEL.PROFILE, OUT/INWARD +1/-1
-- THRESH.JAC., CK RAPID, CONVECTION, EDD TM DEPEND, CHANGE NCND
@ReadData THRMAT,CRAP,CONV,EDTM,CHNCND,Givdtl;
@wterm ' 1';
-- FLOOR FOR R V TP YCARBINV FJ - FH
@ReadData FLOOR(1),FLOOR(2),FLOOR(3),FLOOR(4);
@wterm ' 1';
-- Wacc(R) Wacc(V) Wacc(T) Wacc(RADIAT):
@ReadData Wacc(1), Wacc(2), Wacc(3), Wacc(4);
@wterm ' 1';
-- FitTau TauTol Rvis AQ BQ DRT NRT SCAT
@ReadData FitTau,TauTol,Rvis,AQ,BQ,DRT,NRT,SCAT;
@wterm ' 1';
@ReadData NnTO; -- number of outputs
_do ito=1,NnTO;
Read(Nc,*) TO(ito);
_od;
If (LSystem) then; -- for IBM the read BEGRUN is ignored and redefined here
READ(5,*) IRC; -- for IBM channel 5 is from stack
If (IRC==0) then;
BEGRUN=.FALSE.;
else;
BEGRUN=.TRUE.;
endif;
endif;
%RM:
-- READ MODEL
IF(BEGRUN & IABS(NSTB)==1)THEN;
<*Start: read channel - 9 -- begin run *>;
ELSE;
<*Conr: read Sumprf - channel 10/12 *>;
<*Curv: read Lcurve - channel 11/13 *>;
ENDIF;
Z1=0.d0; Z2=0.d0;
do K=1,Nzon;
Ry(K)=Y(K,1);
-- write(@term,'(a,i3,2(a,1p,e12.4))')' k=',k,' dm=',dm(k),
-- ' Ry=',Ry(i);
Uy(K)=Y(Nzon+K,1); -- standard, change for relativism!
If(abs(KNadap)==5)then;
If(K<NZON)Then;
DM2=DM(K+1);
else;
DM2=DMOUT;
endif;
Z1=Z1+Uy(K)**2*(DM(K)+DM2); -- kin.energy
Z2=Z2+Ry(K)**2*(DM(K)+DM2);
endif;
Ty(K)=Y(2*Nzon+K,1);
enddo;
If(abs(KNadap)==5)then;
URout=sqrt(Z1/Z2);
do K=1,Nzon;
Uy(K)=URout*Ry(K); -- const. homolog. vel.
Y(Nzon+K,1)=Uy(K);
enddo;
endif;
<*nickel: adjust AMNI, XNI for KNadap >0
or read file *.xni for KNadap<0 *>;
IF(CONV)THEN;
do I=1,NZON;
Y(I+(NVARS-2)*NZON,1)=UC(I);
Y(I+(NVARS-1)*NZON,1)=YAINV(I);
enddo;
ENDIF;
NSTMAX=NSTEP+NSTMAX;
IF(MOD(NSTMAX,MBATCH)^=0) THEN;
NSTMAX=(NSTMAX/MBATCH+1)*MBATCH;
-- CALL CMSCOM(IERR,'EXEC ST$MSG WARNING');
WRITE(@Wres,'( '' NSTMAX CHANGED:'',I6)') NSTMAX;
ENDIF;
-- IF(JSTART > MAXORD)THEN;
JSTART=0; --always
-- ELSE; -- seb changed 28 Mar 2001
-- JSTART=-IABS(JSTART); -- IN ORDER TO CHANGE METH, EPS, ETC.
-- Jstart < 0 needed to initialize EPSJ in STIFFBGH
-- ENDIF;
%RM_Start:
Lunit=9;
NFILE=Model;
@wterm '>>',Nfile,'<<';
call StradIO('rm',Lunit,NFILE);
-- Define initial values for rad:
@wterm ' Eko, Ekob: ',Eko,Ekob;
EKO=EKO+EKOB; -- EKO comes from Eve in @IOSTART, but usually = 0
If(EKO>0.)then;
<*Vel: define Uy() for nonzero EKO *>;
endif;
NCND=NZON; -- STRAD defines real NCND and NFRUS
-- any other value for NCND is not allowed here
NFRUS=NFREQ;
KRAD=(NZON-NCND)*NFRUS;
TAUOLD=0.D0;
%RM_Conr:
if(NSTA<0)then;
Lunit=10;
NFILE=Sumprf;
else;
Lunit=12;
NFILE='test.prf';
endif;
call StradIO('cm',Lunit,NFILE);
WRITE(@term,*)' Begrad READ STEP=',NSTEP;
%RM_Curv:
if(NSTA<0)then;
Lunit=11;
NFILE=Sumcur;
else;
Lunit=13;
NFILE='test.crv';
endif;
call StradIO('rc',Lunit,NFILE);
%RM_nickel:
km=1;
print *, 'AMINI=',Amini,UM;
_While (AM(km)-AMini)*UM <= AMNi*1.0000000001d0 & km<nzon _Do
print *, km,AM(km), AMNi*1.0000000001d0;
km=km+1;
_Od; -- (AM(km)-AMini)*UM > AMNi
If ( km >= nzon ) then;
write(@term,*) ' Begrad: AMNI in error! km=',km,'>=nzon =',nzon;
stop;
Endif;
kmnick=km-1;
AMNi=(AM(kmnick)-AMini)*UM;
XNI=XMNI/AMNi; --used in VOLEN
If ( XNI >= 1.d0 & KNadap >0 ) then;
write(@term,*) ' Begrad: AMNI too small ! XNI=',XNI;
stop;
Endif;
If(KNadap>0)then;
<*Fe: adjust Ni to Fe for eq.state *>;
else;
<*Nifor: Ni from a foreign model *>;
endif;
%RM_nickel_Fe:
iferrum=1;
_While iferrum <= Natom & Zn(Iferrum) <> 26 _Do
iferrum=iferrum+1;
_Od;
If (iferrum > Natom ) Then;
write(@term,*) ' Begrad: Ferrum not found !';
stop;
Endif;
_Do km=1,kmnick;
YABUN(iferrum,km)=XNI/AZ(iferrum);
sum=0.d0;
_Do j=1,Natom;
if(j<>iferrum) sum=sum+YABUN(j,km)*AZ(j);
_Od;
If(sum<1.d-5)Then;
write(@term,*) ' Begrad: Sum too small=',sum;
stop;
Endif;
_Do j=1,Natom;
if(j<>iferrum) YABUN(j,km)=YABUN(j,km)*(1.d0-XNI)/sum;
_Od;
_Od;
%RM_nickel_Nifor:
@Wterm 'Read .XNI';
read(28) XNIfor;
/*
read(28)
IOlen1, -- for WATCOM
(XNifor(i),i=1,Nzon)
,IOlen2; -- for WATCOM
*/
Xmni=0.d0;
_Do km=1,Nzon;
Xmni=Xmni+Xnifor(km)*dm(km)*UM;
_od;
write(@term,*)' Ni mass=',Xmni;
%RU:
UR=10.D0**ULGR;
CLIGHT=CS*UTIME/UR;
IF(CRAOLD^=CRAP)THEN;
T=T*(CRAOLD/CRAP);H=H*(CRAOLD/CRAP);
CK1=CK1*(CRAOLD/CRAP);CK2=CK2*(CRAOLD/CRAP);
ENDIF;
UFREQ=BOLTZK*UTP/(2.d0*PI*HPLANC);
-- WRITE(@Wres,*)' UFRC UFREQ:',UFRC,UFREQ;
-- UFRP=BOLTZK/(2.D0*PI*hPlanc);
-- WRITE(@Wres,*)' UFRP UFREQ:',UFRP,UFRP*UTP;
CKRAD=6.D+01/PI**4*CSIGM*UTP**4*UTIME**3/(URHO*UR**3);
CCL=CS*1.D+08/UFREQ; -- TO TRANSFORM FROM ANGSTREMS & V.V.
CFLUX=60.D0*CSIGM*(UTP/PI)**4;
CLUM=32.D0*PI/3.D0*(CSIGM*UR*UTP**4/URHO); -- FL0 INTO LUMINOSITY
CLUMF=4.D0*PI*UR**2*CFLUX; -- FH INTO LUMINOSITY
CIMP =CFLUX*UTIME**2/(CS*UR**2*URHO);
-- Cimp=0.d0;
-- @wterm '!!!!!****>>>> Cimp=0 for TEST with ZERO RaD.FORCE !';
%R_haptab:
-- KNadap<0 for Ni from a foreign model
-- abs(KNadap): 1 freqmn - geom. progression, direct opacity
-- 2 freqmn - Peter Hoeflich, opazr
-- 3 freqmn - Peter Hoeflich, opahoef, hap2int
-- 4,5 freqmn - Ronfict, opahoef, hapintron
-- 6 freqmn - Ronfict, opahoef, hapron0 istim==0)
If (ABS(Knadap) <=3 .and. ABS(Knadap).ne. 1) then; -- see begrad for Knadap
Open(unit=2,file=Opafile,form='unformatted');
Else;
_do nf=1,@ns;
nu=30+nf;
app=char(ichar('0')+nf); -- nf must be <=9
Opafile=Opafile(1:Length(Opafile)-1)//app;
Open(unit=nu,file=Opafile,form='unformatted');
@wterm ' opened unit=',nu,' file=',Opafile;
_od;
Open(unit=29,file=Opafile(1:Length(Opafile)-1)//'ab',form='unformatted');
endif;
_case indop(abs(KNadap)) _of
_1
_2 If(^petread)then;
read(2) Nfreq0,Msta,Nrho,Ntp,YATab,(Wavel(iif),iif=1,@Nfreq),
TpTab,RhoTab,STab,EpsBand,EppBand;
petread=.true.;
endif
_3 _do ihp=1,@ns;
read(30+ihp) nw,Stime,Nfreq0,Msta,Nrho,NTp,
dumWavel,dumFreq,dumFreqmn,
-- read into dummy files
--OLD TpTab,RhoTab,hpbanab2,hpbansc2;
TpTab,RhoTab,hpbansc2; -- actual
_do im=1,Nzon/@skip;
_do iro=1,Nrho;
_do itp=1,NTp;
_do L=1,Nfreq0;
hpsavsc(L,itp,iro,im,ihp)=hpbansc2(L,itp,iro,im);
_od_od_od_od
/* write(@wres,*)' nw,Stime,Nfreq0,Msta,Nrho,NTp', nw,Stime,Nfreq0,Msta,Nrho,NTp;
write(@wres,*)' dumWavel';
write(@wres,*)dumWavel;
write(@wres,*)' dumFreq';
write(@wres,*)dumFreq;
write(@wres,*)' dumFreqmn',dumFreqmn;
stop ' test read ' */;
@wterm ' read unit=',30+ihp;
close(30+ihp);
If(ihp==nw ! ABS(knadap)==6)then;
stmlog(ihp)=log(Stime); -- in days; m.b. better to save in Ronfict?
@wterm ' ihp stmlog=',stmlog(ihp),ihp;
else;
@wterm ' warning error reading haptab, ihp=',ihp,' nw=',nw;
-- here only faked tables are used
-- stop 324;
endif;
_od
read(29) hpbanab2; -- read absorption separately without expansion effect
close(29);
_esac;
@wterm ' Nfreq: ',Nfreq0,Nfreq,' Msta,Nrho,Ntp: ',Msta,Nrho,Ntp;
-- only two sets of tables are in the ROM for
-- two values of dvdr (given by Stime)
-- istim is the flag for the regime of saving the tables at given t
-- for istim==0 , t<t0=Stime(1)/(half of a factor of Stime step),
-- and istim==@ns+1, t>Stime(@ns), use only Stime(@ns) -- the
-- smallest dvdr without interpolation;
--
-- for istim==1 , t0<t<Stime(1), log interp
-- between Stime(@ns) and Stime(1) using t in (0.5, 1)*Stime(1)
--
-- for istim==2, 3, ... @ns, Stime(istim-1)<t<Stime(istim),
-- log interp between boundaries
-- actually Stime is a scalar and the vector is
-- stmlog(.)==log(Stime for a respective moment of time) saved in Begrad
--
tdlog = log(max(t*Utime/8.6d+04,hmin)); -- for Nstep=0, t=0 defined in EVE and
-- saved in *.mod
_select
_ tdlog <= stmlog(1)-(stmlog(2)-stmlog(1))/2.d0
.or. tdlog > stmlog(@ns) .or. abs(Knadap)==6 [istim=0]
_ tdlog > stmlog(1)-(stmlog(2)-stmlog(1))/2.d0 & tdlog <= stmlog(1)
[istim=1]
_other -- here: stmlog (1) < tdlog <= stmlog(@ns)
[ istim=2;
_while stmlog(istim) < tdlog _do istim=istim+1 _od ]
-- stmlog(istim) >= tdlog
_end
@wterm ' istim=',istim;
istold=istim; -- for Strad
_case istim+1 _of
_0 _do im=1,Nzon/@skip;
_do iro=1,@Ntab;
_do itp=1,@Ntab;
_do L=1,Nfreq;
hpbanab1(L,itp,iro,im)=hpbanab2(L,itp,iro,im);
hpbansc1(L,itp,iro,im)=hpbansc2(L,itp,iro,im);
_od_od_od_od
thaplog1=stmlog(1); thaplog2=stmlog(@ns);
-- here the values of thaplog are not important
-- (if they are not equal), since haptab is the same
_1 _do im=1,Nzon/@skip;
_do iro=1,@Ntab;
_do itp=1,@Ntab;
_do L=1,Nfreq;
hpbanab1(L,itp,iro,im)=hpbanab2(L,itp,iro,im);
hpbansc1(L,itp,iro,im)=hpbansc2(L,itp,iro,im);
_od_od_od_od
thaplog1=stmlog(1)-(stmlog(2)-stmlog(1))/2.d0;
thaplog2=stmlog(1)
_2 thaplog1=stmlog(1); thaplog2=stmlog(2)
_3 thaplog1=stmlog(2); thaplog2=stmlog(3)
_4 thaplog1=stmlog(3); thaplog2=stmlog(4)
_5 thaplog1=stmlog(4); thaplog2=stmlog(5)
_esac
If(istim^=0)then; -- for istim==0 already read
app=char(ichar('0')+istim); -- istim must be <=9
Opafile=Opafile(1:Length(Opafile)-1)//app;
Open(unit=30,file=Opafile,form='unformatted');
read(30) nw,Stime,Nfreq0,Msta,Nrho,NTp,
dumWavel,dumFreq,dumFreqmn,
-- read into dummy files
-- TpTab,RhoTab,hpbanab2,hpbansc2;
TpTab,RhoTab,hpbansc2;
close(30);
endif;
If(istim^=0 & istim^=1)then; -- for istim==0 & 1 already read
app=char(ichar('0')+istim-1); -- istim must be <=9
Opafile=Opafile(1:Length(Opafile)-1)//app;
Open(unit=30,file=Opafile,form='unformatted');
read(30) nw,Stime,Nfreq0,Msta,Nrho,NTp,
dumWavel,dumFreq,dumFreqmn,
-- read into dummy files
-- TpTab,RhoTab,hpbanab1,hpbansc1;
TpTab,RhoTab,hpbansc1;
close(30);
close(29);
Open(unit=29,file=Opafile(1:Length(Opafile)-1)//'ab',form='unformatted');
read(29) hpbanab1;
close(29);
endif;
<*check: *>;
%R_haptab_check:
@wterm 'check nw, Stime, Nfreq0, Msta, Nrho, NTp:';
@wterm nw,Stime,Nfreq0,Msta,Nrho,NTp;
_do im=1,Nzon/@skip;
_do iro=1,@Ntab;
_do itp=1,@Ntab;
_do L=1,Nfreq;
hbab=hpbanab2(L,itp,iro,im);
hbal=hpbansc2(L,itp,iro,im);
if (hbal >1.d50 .or. hbab > 1.d50)then;
write(@term,'(4(a,i4),1p,2(a,e12.4))')
' L=',L,' itp=',itp,' iro=',iro,' im=',im,
' hbal=',hbal,' hbab=',hbab;
-- stop ' Begrad: bad happa !!!';
endif;
_od_od_od_od
%RF:
_case indfr(abs(KNadap)) _of
_1 -- our old definition:
FREQ(1)=CS*1.D+08/(WLMAX*UFREQ);
FREQ(NFREQ+1)=CS*1.D+08/(WLMIN*UFREQ);
dumFreq(1) = Cs * 1.d+08 / wlmax;
dumFreq(nfreq+1) = Cs * 1.d+08 / wlmin;
Basis=(WLMAX/WLMIN)**(1.D0/(DBLE(NFREQ)));
-- Geometric PROGRESSION
_Do i=2,NFREQ;
FREQ(i)=FREQ(I-1)*BASIS;
dumFREQ(i)=FREQ(i)*Ufreq;
_OD;
_DO L=1,NFREQ;
FREQMN(L)=SQRT(FREQ(L)*FREQ(L+1));
dumFREQMN(L)=FREQMN(L)*Ufreq;
dumWavel(L) = Cs * 1.d+08 / dumFreqmn(L);
_OD;
--
_2 -- from P.Hoeflich:
If(^petread)then;
read(2) Nfreq0,Msta,Nrho,Ntp,YATab,(Wavel(iif),iif=1,@Nfreq),
TpTab,RhoTab,STab,EpsBand,EppBand;
petread=.true.;
endif;
If(Nfreq>@Nfreq)then;
_Do L=@Nfreq+1,Nfreq;
Wavel(L)=Wavel(L-1)/2.d0;
_od;
endif;
_Do L=1,NFREQ;
FREQMN(L)=CS*1.D+08/(WaveL(L)*UFREQ);
_od;
FREQ(1)=0.5d0*Freqmn(1);
FREQ(Nfreq+1)=2.d0*Freqmn(Nfreq);
_Do L=2,NFREQ;
FREQ(L)=0.5d0*(Freqmn(L-1)+Freqmn(L))
_od;
_3 -- from Ronfict (now the same as _1 but CAUTION here!)
_Do i=1,NFREQ+1;
FREQ(i)=dumFREQ(i)/Ufreq;
_od;
_Do L=1,NFREQ;
FREQMN(L)=dumFREQMN(L)/Ufreq;
_od;
_esac
write(@term,'(2(a,i5))') ' KNadap=',KNadap,
' indfr(abs(KNadap))=',indfr(abs(KNadap));
write(@term,'(a,1pe12.4)') ' Lam max,AA=',CS*1.D+08/(FREQ(1)*UFREQ);
write(@term,'(a,1pe12.4)') ' Lam min,AA=',CS*1.D+08/(FREQ(NFREQ+1)*UFREQ);
-- pause;
-- for all cases:
_DO L=1,NFREQ;
WEIGHT(L)=(FREQ(L+1)-FREQ(L))*FREQMN(L)**3;
IF(L<NFREQ)THEN;
DLOGNU(L)=1.D0/LOG(FREQMN(L)/FREQMN(L+1));
ELSE;
DLOGNU(L)=1.D0/LOG(FREQMN(L)/FREQ(L+1));
ENDIF;
_OD;
%R_freqob:
do L=1,Nfreq;
freqob(L)=freqmn(L);
enddo;
If(Mfreq>Nfreq)then;
do L=Nfreq+1,Mfreq;
freqob(L)=freqmn(Nfreq)*exp(-dble(L-Nfreq)/dlognu(1));
enddo;
endif;
%R_BANDS:
LubvU=NFRUS;
_WHILE LubvU>1 & FREQMN(LubvU) > @U _DO
LubvU=LubvU-1;
_OD; -- FREQMN(LubvU) <= @U OR LubvU==1
LubvB=NFRUS;
_WHILE LubvB>1 & FREQMN(LubvB) > @B _DO
LubvB=LubvB-1;
_OD; -- FREQMN(LubvB) <= @B OR LubvB==1
LubvV=NFRUS;
_WHILE LubvV>1 & FREQMN(LubvV) > @V _DO
LubvV=LubvV-1;
_OD; -- FREQMN(LubvV) <= @V OR LubvV==1
LubvR=NFRUS;
_WHILE LubvR>1 & FREQMN(LubvR) > @R _DO
LubvR=LubvR-1;
_OD; -- FREQMN(LubvR) <= @R OR LubvR==1
LubvI=NFRUS;
_WHILE LubvI>1 & FREQMN(LubvI) > @I _DO
LubvI=LubvI-1;
_OD; -- FREQMN(LubvI) <= @I OR LubvI==1
Lyman=NFRUS;
_WHILE Lyman>1 & FREQMN(Lyman) > @Lyman _DO
Lyman=Lyman-1;
_OD; -- FREQMN(LubvV) <= @V OR LubvV==1
WRITE(@Wres,*)' L UBVRI Lyman:',LubvU,LubvB,LubvV,LubvR,LubvI,Lyman;
%RM_Start_Vel:
<*tri: define initial velocity triangle profile i1, i2, i3 *>;
Z1=US/(AM(I2)-AM(I1));
IF(I2^=I3) Z2=US/(AM(I2)-AM(I3));
_Do K=I1,I3;
If(K<=I2)then;
Uy(K)=Z1*(AM(K)-AM(I1));
else;
Uy(K)=Z2*(AM(K)-AM(I3));
endif
_od;
Z1=0.D0;
K1=I1+1;
_Do K=K1,I3;
if(K<Nzon)then;
Z1=Z1+Uy(K)**2*(DM(K)+DM(K+1));
else;
Z1=Z1+Uy(K)**2*(DM(K)+DMout);
endif;
_od;
_Do K=I1,I3;
Y(Nzon+K,1)=Uy(K)*2.D0*SQRT((EKO/UEPRI)/Z1)
_od;
%RM_Start_Vel_tri:
i=0;
_repeat
i=i+1;
_until (am(i)-amini)/(amout-amini) >=ai1;
i1=i;
write(@term,'(a,i5,1p,e10.3)') ' i1, ai1:', i1, ai1;
_repeat
i=i+1;
_until (am(i)-amini) /(amout-amini) >=ai2;
i2=i;
write(@term,'(a,i5,1p,e10.3)') ' i2, ai2:', i2, ai2;
_repeat
i=i+1;
_until (am(i)-amini) /(amout-amini) >=ai3 .or. i==nzon;
i3=i;
write(@term,'(a,i5,1p,e10.3)') ' i3, ai3:', i3, ai3;
%RH:
WRITE(@Wres,'(''%RUN:'')');
WRITE(@Wres,'(//30X,''<===== HYDRODYNAMIC RUN OF MODEL '',A,
''=====>'',/)') Sumprf;
WRITE(@Wres,'(30X,''MASS(SOLAR)='',F6.3,7X,''RADIUS(SOLAR)='',F9.3/,
30X,''EXPLOSION ENERGY(10**50 ERG)='',1P,E12.5,/)')
AMOUT*UM,RADBEG,EKO+Eburst;
WRITE(@Wres,'(30X,''<====='',33X,''=====>'',//)');
WRITE(@Wres,'('' INPUT PARAMETERS '')');
WRITE(@Wres,'('' EPS = '',F15.5,9X,'' Rce = '',1P,E15.5,A)')
EPS, RCE*UR/RSOL,' SOL.Rad.';
WRITE(@Wres,'('' HMIN = '',1P,E15.5,A,7X,'' AMht = '',1P,E15.5,A)')
HMIN*UTIME,' S', AMht,' SOL.MASS';
WRITE(@Wres,'('' HMAX = '',1P,E15.5,A,7X,'' Tburst= '',1P,E15.5,A)')
HMAX*UTIME,' S', TBurst,' S';
WRITE(@Wres,'('' THRMAT= '',1P,E15.5,9X,'' Ebstht= '',1P,E15.5,A)')
THRMAT,EBurst,' 1e50 ergs';
WRITE(@Wres,'('' METH = '',I15,9X,'' CONV = '',L15)')METH,CONV;
WRITE(@Wres,'('' JSTART= '',I15,9X,'' EDTM = '',L15)')JSTART,EDTM;
WRITE(@Wres,'('' MAXORD= '',I15,9X,'' CHNCND= '',L15)')MAXORD,CHNCND;
WRITE(@Wres,'('' NSTA = '',I15,9X,'' FitTau= '',1P,E15.5)')NSTA,FitTau;
WRITE(@Wres,'('' NSTB = '',I15,9X,'' TauTol= '',1P,E15.5)')NSTB,TauTol;
WRITE(@Wres,'('' NOUT = '',I15,9X,'' IOUT = '',I15)')NOUT,IOUT;
WRITE(@Wres,'('' TcurA = '',F15.5,9X,'' Rvis ='',F15.5)')
TcurA,Rvis;
WRITE(@Wres,'('' TcurB = '',F15.5,9X,'' BQ = '',1P,E15.5)')TcurB,BQ;
WRITE(@Wres,'('' NSTMAX= '',I15,9X,'' DRT = '',1P,E15.5)')NSTMAX,DRT;
WRITE(@Wres,'('' XMNI = '',1P,E15.5,A,'' NRT = '',I15)')
XMNI,' SOL.MASS',NRT;
If(KNadap>0)then;
WRITE(@Wres,'('' XNI = '',1P,E15.5)')XNI;
WRITE(@Wres,'('' CONTM.= '',1P,E15.5,A,'' SCAT = '',L15)')
AMNI,' SOL.MASS',SCAT;
else;
WRITE(@Wres,'('' XNifor= '',1P,E15.5)')XNifor(1);
WRITE(@Wres,'('' MNicor= '',1P,E15.5,A,'' SCAT = '',L15)')
AMNI,' SOL.MASS',SCAT;
endif;
PRINT @F1,UTP,UTIME,URHO,UFREQ;
PRINT @F2,CK1,CK2,CFR,CKRAD,CFLUX,CLUM,CLUMF;
@F2: FORMAT(' CK1 =',1P,E12.5,' CK2=',E12.5,' CFR=',E12.5,
' CKRAD=',E12.5/
' CFLUX=',E12.5,' CLUM=',E12.5,' CLUMF=',E12.5);
@F1: FORMAT(' UTP=',1P,E12.5,' UTIME=',E12.5,' URHO=',E12.5,
' UFREQ=',E12.5);
PRINT'('' *****CK1, CK2 INCREASED'',G8.1,'' TIMES*****'')',CRAP;
WRITE(@Wres,'('' FLOOR :''/1X,1P,10E10.2)') FLOOR;
PRINT *,' N DIFF. EQS=',N,' N RAD. EQS=',2*KRAD;
%S:O
SUBROUTINE @SAVE_MODEL;
IMPLICIT REAL*8 (A-H,O-Z);
INTEGER LUNIT;
Character*80 NFILE*80;
_INCLUDEN snrad;
_INCLUDEN abo;
if(NSTB<0)then;
Lunit=10;
NFILE=Sumprf;
else;
Lunit=12;
NFILE='test.prf';
endif;
call StradIO('wm',Lunit,NFILE);
--close(Lunit);
if(NSTB<0)then;
Lunit=11;
NFILE=Sumcur;
else;
Lunit=13;
NFILE='test.crv';
endif;
call StradIO('wc',Lunit,NFILE);
--close(Lunit);
RETURN;END;
%I:O
BLOCK DATA STDATA;
IMPLICIT REAL*8(A-H,O-Z);
_include zone;
COMMON/QNRGYE/QNUC,RGASA,YELECT;
COMMON/HNUSED/HUSED,NQUSED,NFUN,NJAC,NITER,NFAIL;
-- COMMON/AQ/AQ,BQ,DRT,NRT; -- BQ pseudo-viscosity for R-T
COMMON/AZNUC/ACARB,ZCARB,ASI,ZSI,ANI,ZNI,QCSI,QSINI;
COMMON/TOO/TOO,KO,KNTO,TO(KOMAX),STO(KOMAX),NTO(KOMAX);
-- DATA AQ/2.D0/; -- standard
-- DATA AQ/8.D0/; -- enhanced
<*A: DATA FOR TOO *>;
<*P* DATA FOR Peter's prescribed moments *>;
DATA QNUC,YELECT/8.8861D+6,0.5D0/;
DATA ACARB,ZCARB,ASI,ZSI,ANI,ZNI,QCSI,QSINI/12.D0,6.D0,28.D0,
14.D0,56.D0,28.D0,88.861D0,1.884D0/;
DATA NFUN,NJAC,NITER,NFAIL/4*0/;
END;
%IA:O
DATA TO/komax*1.D+19/,STO/komax*1.D-5/,NTO/komax*0/;
-- no prescribed moments
/* DATA TO/
8.0000d-02, 8.0116d-02, 8.0231d-02,
8.0289d-02, 8.0347d-02, 8.0463d-02,
8.0694d-02, 8.0926d-02, 8.1157d-02,
8.1620d-02, 8.2083d-02, 8.2315d-02,
8.3000d-02, 8.5000d-02, 8.8000d-02,
9.2000d-02, 9.6000d-02, 1.0000d-01,
1.0500d-01, 1.1000d-01, 1.2000d-01,
1.4000d-01, 1.6000d-01, 2.0000d-01,
2.5000d-01, 3.0000d-01, 4.0000d-01,
5.0000d-01, 7.0000d-01, 1.0000d+00,
1.5d+00, 2.d0, 3.d+00, 4.d0, 5.d+00,
6.d0, 8.d0, 10.d0, 30.d0, 90.d0 -- in days
/, STO/komax*1.D-5/, NTO/komax*0/;
*/ -- STO always in secs
%IP:O
DATA TO/
8.0000d-02, 8.0116d-02, 8.0231d-02,
8.0289d-02, 8.0347d-02, 8.0463d-02,
8.0694d-02, 8.0926d-02, 8.1157d-02,
8.1620d-02, 8.2083d-02, 8.2315d-02,
8.3000d-02, 8.5000d-02, 8.8000d-02,
9.2000d-02, 9.6000d-02, 1.0000d-01,
1.0500d-01, 1.1000d-01, 1.2000d-01,
1.4000d-01, 1.6000d-01, 2.0000d-01,
2.5000d-01, 3.0000d-01, 4.0000d-01,
5.0000d-01, 7.0000d-01, 1.0000d+00,
1.5000d+00, 3.0000d+00, 5.0000d+00/, -- zdesj momenty v
-- v sutkah ili sekundah
STO/komax*0.D0/,
NTO/komax*0/;
English
Русский