-- 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/;