-- _outcom ;
-- _trace '--------';
-- calling gdepos (6) each step after some time
-- calling tt4strad and Lbol in one run
/* PACKAGE STELLA: IMPLICIT HYDRODYNAMICS, KINETICS & VARIOUS TYPES
OF TRANSFER.
Version R4 continuum PHOTON TRANSFER without REZONING
but with dropping of outer zones for tau<@tau_floor
and limiting the flux for zones with jumps in tau > taujmp
12 Mar 1989 at ITEP, Moscow
27 May 1990 at MPA, Garching
18 Dec 1991 at MPA, Garching -- scattering
24 Apr 2001 at UCO/Lick, Santa Cruz -- call beams -- now in eddirel.trf
subroutine references:
**** BEGRAD ****
BEGIN - initialization of run.
**** LBALRAD ****
LOSSEN - losses and gains of energy
PRIBAL - print model & balances
**** STIFFBGH ****
VTSTIF - integrates one time step
**** DFRAD ****
DIFMAT for HydroDynamics & Transfer EQUATIONS--
**** HAPBLA ****
HAPPA -- opacity in frequency bands
**** SAHAZ ****
URSOS
**** OPAZR ****
OPACIT -- Rosseland mean (Zeldovich-Rayzer approximation)
**** NTH ****
Subroutine NTH: finds in NTHICK(L) the zone numbers of
optically thick zones for all NFRUS frequency groups *>;
**** commons are included by trefor _include ****
fundrad include a
SNRAD INCLUDE A
ABO INCLUDE A
BLACK INCLUDE A;
BLACK, BLACKD -- black-body radiation occ. number & its derivative
*/
-- _TRACE "WRITE(@term,*)' jstart=',jstart,"
-- _TRACE "WRITE(@term,*)' Ncnd,Nfrus=',Ncnd,Nfrus,"
-- _TRACE "WRITE(@term,*)' Y11=',Y(1,1),"
-- _TRACE "WRITE(@term,'(a,1p,e12.4)')' chem=',chem,"
-- _TRACE "@wterm' Strad Ltau,Nthick=',Ltau,Nthick(Ltau),"
-- _TRACE "@wterm' Strad Ltau,=',Ltau,"
_Define
@Lfn ' ' -- for Unix
-- @Lfn ',80' -- length of the Filein, Fileout
_DEFINE
@INTERRUPTION INTRP
@FREQ_TOL 1.6D0
@FJ2 Y(NVARS*NZON+(NZON-Ncnd)*(L-1)-Ncnd+KM+1,1) -- FJ(Km+1,L)
@FJ1 Y(NVARS*NZON+(NZON-Ncnd)*(L-1)-Ncnd+KM,1) -- FJ(Km,L)
@FH1 Y(NVARS*NZON+(NZON-Ncnd)*(L-1)-Ncnd+KM+KRAD,1)
@NCOND_MAX NZON-3
@TIME CT001 -- CURRENT TIME
@DATE D0001 -- CURRENT DATE
@BEG_TIME B0001
@SAVE_MODEL SAVEM
@TAU_FLOOR 3.D-4
@Nrezone 8
@N_EDD 50 -- N of steps for updating Eddington factors
@Nrlx 25.d0 -- N of steps for Edd relaxation (real)
@LEFT 6.D0 -- COEFF. FOR LEFT BOUNDARY AT R=RCORE
-- @LEFT==6*HEDD AT CORE WHERE H==.5*INTEGRAL(F*MU*DMU)==FH*HEDD
@DateTime "If(Lsystem) then;
CALL DATTIM(@DATE,@TIME);
else;
call DateTime(@Date,@Time);
endif;"
--@DateTime ";" ;
-- DEBUG SUBCHK,UNIT(6);ENDDEBUG;
PROGRAM STELLA; -- CALCULATE STELLAR MODEL
<*usemodules: include use stantments*>;
IMPLICIT REAL*8(A-H,O-Z);
type(cp_para_env_type), pointer :: para_env;
type(cp_error_type) error;
integer mpi_comm_default;
integer :: istat, ierr = 0;
integer, parameter :: default_output_unit = 6;
logical :: init_mpi = .true.;
integer handle;
character(len=*), parameter :: subrtn_name = 'prog STELLA';
REAL LTIM,MTIM; -- TIMER
REAL @BEG_TIME,TIMEND; -- TIMER IBM, F90 & Ndp
--REAL*8 @BEG_TIME,TIMEND; -- TIMER RISC_IBM
--CHARACTER*8 @DATE, @TIME*8; -- TIMER IBM
CHARACTER*10 @DATE;
CHARACTER*8 @TIME; -- TIMER F90
_INCLUDE snrad;
_INCLUDE opacity;
_INCLUDE abo;
-- real*8 hptrbab(Nfreq,@Ntab,@Ntab,Mzon/@skip),
-- hptrbsc(Nfreq,@Ntab,@Ntab,Mzon/@skip);
-- Parameter(Taukli=0.1d0); -- tau for flux limiter
-- Parameter(Taujmp=0.5d0); -- tau jump for flux limiter
-- Common/Lim/Klim(Mzon);
Common/Volanm/Ryzer,Kmcor; -- emulate Ni for analyt.model
PARAMETER(NzMIN=Mzon*5/6); -- minimum Nzon - number of mass zones
PARAMETER(NCNMIN=5); -- minimum Ncnd - number of conduction zones
Parameter(TREZON=3.D18); -- allow rezoning after Trezon (sec)
-- Parameter(TREZON=3.D1); -- allow rezoning after Trezon (sec)
Parameter(tsmoth=4.D6); -- smoothing/forced Ncnd=0 (sec, proper tim)
Parameter(tret=1.d0/3.d0);
Logical HOLDFR; -- Decreasing NFRUS is forbidden
Logical LEXIST1,LEXIST2;
Character*80 Dfile,Rfile,bmfile,Nidist,swdet, wMaratfile,trdet;
character*160 wordm(10),runname,filestr;
integer lwordm(10),lrun,nwordm,errnum;
Character*8 STRBAT;
-- Dimension dumFreq(Nfreq+1),dumFreqmn(Nfreq),dumwavel(Nfreq);
-- now in opacity.inc
Common/FRAD/ FRADJ(Mzon,Nfreq),FRADH(Mzon,Nfreq);
Common/out/ tretard,tout,Yout(NYdim); -- for interpolation
-- Dimension Wacc(3); -- for Wacc(IFL)<1. accuracy is higher
-- Data Wacc /1., 1., 1./;
-- DATA freqfind/.false./;
<*C* CONVECTIVE VARIABLES *>;
_Include black;
-- BLACK, BLACKD -- black-body radiation occ. number & its derivative
-- Theta(Z)=0.5D0*(SIGN(1.D0,Z)+1.D0);
-- *************** UNIT ************************
nullify(para_env);
call init_stl(init_mpi, ierr);
call cp_para_env_create(para_env, group=mpi_comm_default, owns_group=.false.,error=error);
-- *************** END UNIT ************************
CALL timeset(subrtn_name,'I','',handle);
<*openfiles: unit=1 strad.1
unit=8 dat file
unit=28 Nidist etc. *>;
call Vtime(@BEG_TIME);
If(Lsystem) then;
-- CALL SPM001(12); -- underflow correction
-- LTIM=LTIME(0); -- IBM
else;
Call Tremain(LTIM,@BEG_TIME); -- Run time remaining
endif;
<*I: INITIALIZATION: call BEGIN *>;
If ( Givdtl ) Then; -- for test runs
Open(unit=18,file='details.res');
Endif;
<*A: COMPUTE NEW MODEL CALLING VTSTIF, ANALYSE THE STEP,
CHECK BALANCES & PRINT CALLING PRIBAL, WRITE WHERE NEEDED.
For successful runs KFLAG==0 in the end *>;
<*closecrvfiles: check if crv flx are closed by savemodel - stradio *>;
<*closefiles: final closing all files *>;
if (KFLAG==0)then;
call tt4strad(runname);
call Lbol(runname);
endif;
<*packing* using command SYSTEM() send
.tt, .ph, .swd, .lbol to ~/prg/stella/res/
gzip .res, rm .bm, .avg etc *>;
CALL timestop(0.0_dp,handle);
CALL timeprint (6,1.e-1_dp);
STOP ' Normal stop of stella';
END;
%_openfiles:
Open(unit=1,file='strad.1',status='old');
READ(1,'(A)'); -- header
read(1,'(A)') filestr;
close(1);
<*filenames: filestr contains filenames for input and output *>;
@wterm Dfile;@wterm Rfile;@wterm Model;
@wterm Nidist;@wterm Opafile;@wterm Sumprf;@wterm Sumcur;
INQUIRE(FILE=Sumprf,EXIST=LEXIST1);
INQUIRE(FILE=sumcur,EXIST=LEXIST2);
IF(LEXIST1 & LEXIST2)THEN;
BEGRUN=.FALSE.;
ELSEIF(^LEXIST1 & ^LEXIST2) THEN;
BEGRUN=.TRUE.;
ELSE;
Write(@Term,'(A)') ' File Sumprf or SUMCUR is missing';
stop 28;
ENDIF;
write(@term,*) Lexist1,Lexist2,BegRun;
Open(unit=8,file=Dfile,status='old'); -- Data file for start
Open(unit=@WRES,file=Rfile,status='new'); -- Results
Open(unit=@WMarat,file=wMaratfile,status='new'); -- Results
Open(unit=60,file=bmfile,status='new'); -- file for beams
write(@term,*) ' opened unit=60 file=',bmfile;
-- read(*,*);
Open(unit=81,file=swdet,status='new'); -- SW details
-- Open(unit=82,file=trdet,status='unknown'); -- for A.Tolstov
Open(unit=28,file=Nidist,form='unformatted',status='unknown');
-- open(450,file='findTradByMNKLum.dat');
-- write(450,*)'Hello'; stop;
-- Open(unit=7,file='delzon.res');
%_openfiles_filenames:
call words(filestr, wordm,lwordm, nwordm);
@wterm filestr;
@wterm nwordm;
if (nwordm < 5 ) then; --
@wterm ' check your start file, something is missing !';
stop 12;
endif;
lrun=lwordm(1);
runname(1:)=wordm(1)(1:lrun);
rfile(1:lwordm(2))=wordm(2)(1:lwordm(2));
Narg=Iargc();
_Case Narg+1 _Of
write(*,*) ' Extra arguments.'; Stop 16;
_0 write(*,*) ' using runname from strad.1 file';
_1 write(*,*) ' using runname from argument';
call GetArg(1,wordm(1) @Lfn);
runname=wordm(1)(1:len_trim(wordm(1)));
rfile=runname(1:len_trim(runname))//'.res';
lrun=len_trim(runname);
lwordm(3)=8; -- for c******m models
wordm(3)=runname(1:lwordm(3))//'_2000h';
lwordm(3)=lwordm(3)+6; -- for '_2000h' models
lwordm(4)=lwordm(3);
wordm(4)=wordm(3);
_Esac;
dfile(1:)=runname(1:lrun)//'.dat';
wMaratfile(1:)=runname(1:lrun)//'.mrt';
bmfile(1:)=runname(1:lrun)//'.bm';
model(1:)='../eve/'//wordm(3)(1:lwordm(3))//'.mod';
Nidist(1:)='../eve/'//wordm(3)(1:lwordm(3))//'.xni';
-- Nidist(1:)='../eve/'//wordm(4)(1:lwordm(4))//'.xni';
Opafile(1:)='../vladsf/'//wordm(5)(1:lwordm(5)); -- standard path
Sumprf(1:)=runname(1:lrun)//'.prf';
SumCur(1:)=runname(1:lrun)//'.crv'; -- opened in stradio unit=11
Depfile(1:)=runname(1:lrun)//'.dep';
Flxfile(1:)=runname(1:lrun)//'.flx'; -- opened in stradio unit=23
swdet(1:)=runname(1:lrun)//'.swd';
%_closefiles:
Close(unit=2);
Close(unit=8);
Close(unit=@WRES);
Close(unit=@WMarat);
Close(unit=9);
Close(unit=10);
Close(unit=11); -- crv
Close(unit=22); -- deposition
close(unit=23); -- flx
close(unit=60); -- bm
close(unit=81); -- swdet
-- close(unit=82); -- trdet
-- close(450); -- findTradByMNKLum
If ( Givdtl ) Then;
Close(unit=18);
Endif;
-- _do nf=1,@ns;
-- nu=30+nf;
-- close(nu); -- closed immediately
-- _od;
-- Close(unit=7); -- delzone.res
%_packing:
/*
Isys = SYSTEM('gzip '//runname(1:lrun)//'.res');
write(*,*) ' gzip res is done Isys=',Isys, ' Csys=',Csys;
If (Isys == -1) then;
errnum = ierrno( );
write(*,'(a,i2)') ' Error gzip res= ', errnum;
endif;
Isys = SYSTEM('mv '//runname(1:lrun)//'.{ph,tt,swd,lbol,res.gz}'
*/
Isys = SYSTEM('mv '//runname(1:lrun)//'.{ph,tt,swd,lbol,res}'
//' ../../res/');
write(*,*) ' mv to res is done Isys=',Isys, ' Csys=',Csys;
If (Isys == -1) then;
errnum = ierrno( );
write(*,'(a,i2)') ' Error mv to res= ', errnum;
endif;
Isys = SYSTEM('rm -f '//runname(1:lrun)//'.{bm,avg}');
write(*,*) ' rm bm,avg is done Isys=',Isys, ' Csys=',Csys;
If (Isys == -1) then;
errnum = ierrno( );
write(*,'(a,i2)') ' Error rm bm,avg = ', errnum;
endif;
%C:O
EQUIVALENCE (Y((NVARS-1)*Mzon+1,1),YAINV(1));
EQUIVALENCE (Y((NVARS-2)*Mzon+1,1),UC(1));
COMMON/CONV/ELMX(Mzon),FLCON(Mzon);
%I:
EIT=1.D-5;<*ACCURACY OF ITERATIONS*> DST=1.D-4;
CALL BEGIN;
close(8); -- this is needed to detach dat file for use by ttfit
Close(unit=28); -- to detach xni
FLnois=Floor(NVARS+1); -- for radiation
If(Lsystem) CALL BASTAT('BEGIN RD');
PRINT *,' JSTART=',JSTART;
NEEDBR=.TRUE.; -- F01BRF/M28Y12 NEEDED
NSTEP0=NSTEP;
Call LOSSEN;
<*T: CHECK Ncnd & KRAD *>;
RTphi(0)=0.; RTphi(Mzon+1)=0.;
Uscale=1.d9*Utime/UR;
<*RT: R-T coeff. in RTphi *>;
<*Tpcomp: find in Tpcomp the highest Tp(i) for i>Ncnd
for use in approximate treatment of Compton *>
Ycomp=(Boltzk*UTp*Tpcomp)/(Amelec*Cs**2); -- Comptonization par.
Write(@term,*)' Tpcomp Ycomp: ',Tpcomp,Ycomp;
<*N: CHECK NFRUS *>;
_do L=1,NFREQ;
_DO IK=1,NZON; -- for EDTM==T for all later steps
Eddn(IK,L)=tret;
EddO(IK,L)=Eddn(IK,L);
EDDJ(IK,L)=Eddn(IK,L);
_OD;
HEdN(L)=0.5d0;
HEdO(L)=HEdN(L);
HEDD(L)=HEDN(L);
write(*,'(a,i5,1p,e12.4)')' 0 step L Hedd:', L,Hedd(L);
_od;
%IT:
tday=t*UTIME/8.64d+04; -- CRAP=1.
if(NSTEP==0)t_eve=tday; -- needed by happa (hapdirect)
/*Ncnew=NTHICK(1); LTAU=1;
_DO L=2, NFRUS;
IF( NCNEW > NTHICK(L)) THEN;
NCNEW=NTHICK(L);
LTAU=L; -- find the most transparent energy group
ENDIF;
_OD; */
@wterm' Ltau=',Ltau;
I=NZON; TAU(I)=0.D0;
-- _WHILE I>1 & (TAU(I)<=FitTau*TauTol ! I>=Ncnd ) _DO
_WHILE I>1 _DO
TP=Ty(I); Pl=3.D0*DM(I)/(Ry(I)**3-Ry(I-1)**3);
_Do ii=1, Natom [ Yat(ii)=YABUN(ii,I) ];
RADP=.FALSE.;
CALL URSOS;
Rho(I)=Pl;
Press(I)=P*UPI;
kmhap=I; CALL HAPPA;
-- @wterm' after CALL HAPPA kmhap=',kmhap;
Hapmin=Happal(1);
Ltau=1;
_Do L=2,Nfrus;
If(Happal(L)<Hapmin)then;
Hapmin=Happal(L);
Ltau=L;
endif;
_od;
-- @wterm' Ltau=',Ltau;
TAUste=(Ry(I)-Ry(I-1))*HAPPAL(LTAU)*PL;
If(Scat)then;
TAUabs=(Ry(I)-Ry(I-1))*HAPabs(LTAU)*PL;
TAUsca=Tauste-Tauabs;
---If(Nstep>7250)then;
-- @wterm' Ltau=',Ltau;
-- Write(@term,'(A,1P,3E11.3)')' tau tauabs tausca',tauste,tauabs,tausca;
---endif;
If(Tausca>1.d0)then;
Tauste=sqrt(Tauabs*Tausca);
else;
Tauste=Tauabs;
endif;
endif;
TAU(I-1)=TAU(I)+Tauste;
I=I-1
_OD; -- I=0 --old I<Ncnd & TAU(I)>FitTau*TauTol
Ncnew=Nthick(Ltau);
-- _Do K=Ncnew-1,I,-1; -- old wrong for small diff. at ncnew & I
_Do K=Ncnew-1,1,-1;
tau(K)=tau(K)-tau(Ncnew); -- to try to catch transparent zones < Ncnew
_od;
-- Ncnew=NZON;
_WHILE TAU(max(Ncnew,1)) <= FitTau & Ncnew>=NCNMIN _DO
-- & TAU(Ncnew-1) <= FitTau*TauTol _DO
Ncnew=Ncnew-1
_OD; -- TAU(Ncnew) > FitTau or Ncnew<NCNMIN or
-- TAU(Ncnew-1) > FitTau*TauTol
-- WRITE(@term,*)' Ncnew=',Ncnew,' Ncnd=',Ncnd,
-- ' TAU(Ncnd or 1)=',TAU(max(Ncnd,1));
/* WRITE(@term,'(2I6,I3,4(A,1P,G10.3))')
NSTEP,MAXER,NQUSED,' H=',H,' Hused=',HUSED,
' Tp=',Ty(1),' Rho=',Rho(1);
WRITE(@term,'(1X,A,1P,G11.3,A,I3,A,G11.2)')
' Time(days)=',tday,' Ncnd=',Ncnd,' TAU(Ncnd or 1)=',TAU(max(Ncnd,1));
*/
WRITE(@term,'(2I6,I3,3(A,1P,G10.3),A,I3,A,G11.2)')
NSTEP,MAXER,NQUSED,' Hu=',HUSED,' Ro=',Rho(1),
' td=',tday,' Ncnd=',Ncnd,' Tau=',TAU(max(Ncnd,1));
IF(Ncnew<NCNMIN) Ncnew=0;
--If(Ncnew<NCNMIN .or. t*Utime>tsmoth) Ncnew=0; -- force Ncnd=0
Ncnew=MIN(Ncnew,@NCOND_MAX);
HOLDFR=.TRUE.; -- forbids the decrease of NFRUS if there is no
-- change of Ncnd or if Ncnd grows (i.e. hot wave
-- moves out)
_If (TAU(max(Ncnd,1))>TauTol*FitTau ! -- Ncnew==0 ! -- for forced Ncnd
TAU(max(Ncnd,1))*TauTol<FitTau) & Ncnew^=Ncnd & (CHNCND ! Nstep==0)_THEN
-- TAU(Ncnd)*TauTol<FitTau) & Ncnew< Ncnd & CHNCND _THEN
_Select
_ Ncnew<Ncnd
[ <*O: ADD FJ,FH FOR ZONES Ncnew+1:Ncnd *>;
HOLDFR=.FALSE. ] -- for Ncnd moving inward
-- the decrease of Nfrus is allowed
_ Ncnew>Ncnd & Ncnew<=@NCOND_MAX
[ <*U: DELETE FJ,FH FOR ZONES Ncnd+1:Ncnew *>]
_ Ncnew>@NCOND_MAX [ KFLAG=-3 ] -- not active now
_end;
If(KFLAG==0)Then;
<*LIM* limiting FH *>;
<*D* test printing*>;
-- test as coll:
JSTART=0;
-- Hmax=max(Hmax,t*1.d-3);
@wterm ' 2nd Jstart=',Jstart;
NEEDBR=.TRUE.; -- F01BRF/M28Y12 NEEDED
-- Evalja=.TRUE.; -- not needed when Jstart=0;
<*H* CHECK FH & PUT SMALLER THAN FJ WITH JSTART=0 *>;
-- XJPH=RPH; -- for limiting DFJ in Dfrad %Derata
endif;
_fi;
-- If(Ncnd>0) Nclast=Ncnd;
%A:
@DateTime;
WRITE(@Wres,'(1X,'' ====> STEP:'',I6,2X,A,2X,'' ON '',A)')
NSTEP,@TIME,@DATE;
-- IF(Ncnd==0)THEN;
-- IF(t*UTIME > tsmoth) CALL SMOOTH;
<*eddtau: eddington factors and optical depth to start with
call pribal to print the first model of the run*>;
-- Klim=min(I+1,Nzon);
<*HLIM* *>;
_do kko=1,komax;
TO(kko)=TO(kko)*8.64d+04; -- if TO in days
_od;
KFLAG=0;
-- KO=1; -- wrong for continuation runs
KNTO=1;
-- TOO=TO(1); -- here TO in secs
nmom=1;
TOO=TO(1);
/*
_while t*Utime > TO(nmom) _do -- this loop must be corrected for R/c
-- if below we output interpolation for t_observed
-- maybe better return to interpolation outputs in comoving time
@wterm' t*Utime > TO(n)',t*Utime,TO(nmom),nmom;
-- pause;
nmom=nmom+1;
TOO=TO(nmom);
_od;
*/ --
KO=nmom;
_WHILE KFLAG==0 & NSTEP<NSTMAX _DO
<*ymax: define for all y(.) from 1 to N *>;
tday=t*UTIME/8.64d+04; -- CRAP=1.
if(NSTEP==0)t_eve=tday;
_If abs(KNadap)>=4 _then
<*haptab: read two relevant tables of hptab{ab/sc}, depending
on the interval in which log(tday) falls in stmlog(.)
array and save in hpban{ab/sc}1 & hpban{ab/sc}2
respectively *>;
_fi
CALL VTSTIF;
NSTEP=NSTEP+1;
NZM=NZMOD;
_Do K=1,Nzon;
Ry(K)=Y(K,1);
Uy(K)=Y(Nzon+K,1); -- NR
-- Uy(K)=Y(Nzon+K,1)/sqrt(1.d0+(Y(Nzon+K,1)/clight)**2);
Ty(K)=Y(2*Nzon+K,1);
_Od;
if(tday > 2.d0) CALL GDEPOS; -- each step after 2 days
-- CALL GDEPOS; -- each step
<*L: analyse the step: find new regime in NREG(1:NZON):
NREG(I)==1 Carbon burning
2 Silicon burning
3 no burning
find integral energy losses,
print & write model when needed*>;
--WRITE(@term,*)' NSTEP=', NSTEP;
IF (KFLAG==0) THEN ;
-- WRITE(STRBAT,'(''step'',I5)') NSTEP;
-- CALL BASTAT(STRBAT);
If(t*UTIME>=8.64D4*tcurB) NSTMAX=NSTEP;
If((MOD(NSTEP,MBATCH)==0 ! NSTEP==NSTMAX) & ^Givdtl ) then;
CALL @SAVE_MODEL; -- File OUT RES
@DateTime;
WRITE(@Wres,'(1X,'' ====> SAVED:'',I6,2X,A,2X,'' ON '',A)')
NSTEP,@TIME,@DATE;
call Tremain(MTIM,@BEG_TIME);
-- IF(MTIM*10<(LTIM-MTIM)*12) KFLAG=-4; -- KFLAG = -4 for time limit
LTIM=MTIM;
call timeprint (@Wres,1.e-1_dp);
endif;
ENDIF;
_OD;
_Case -KFLAG _of
Write (6,*)' finish with KFLAG=',KFLAG;
_1 Write (6,*)' the error test fails for H=HMIN: KFLAG=',KFLAG;
CALL PRIBAL(0)
_2 Write (6,*)' corrector iterations fail to converge for',
' H=HMIN: KFLAG=',KFLAG;
Write (@term,*)' corrector iterations fail to converge for',
' H=HMIN: KFLAG=',KFLAG;
CALL PRIBAL(0)
_3 WRITE(@Wres,*)' Ncnew=',Ncnew,'>=@NCOND_MAX'
_4 WRITE(@Wres,*)' time limit: KFLAG=',KFLAG
_esac;
IF(KFLAG<0) CALL PRIBAL(1);
CALL @SAVE_MODEL; -- File OUT RES
@DateTime;
WRITE(@Wres,'(1X,'' ====> SAVED:'',I6,2X,A,2X,'' ON '',A)')
NSTEP,@TIME,@DATE;
CALL VTIME(TIMEND);
TIMEND=(TIMEND-@BEG_TIME);
WRITE(@Wres,'('' RUN_TIME(s):'',1P,G12.3)')TIMEND;
@DateTime;
WRITE(@Wres,'(1X,'' ====> STEP:'',I6,2X,A,2X,'' ON '',A)')
NSTEP,@TIME,@DATE;
%A_eddtau:
CALL @FEAUTRIER(Ncnd,Ncnd);
trlx=@Nrlx*h; @wterm trlx;
_do L=1,Nfrus;
_DO IK=Ncnd+1,NZON; -- for EDTM==T
EddO(IK,L)=EddN(IK,L);
EddJ(IK,L)=EddN(IK,L);
if(L==25) write(*,'(a,i5,1p,3e12.4)')' L Edd{O,J,N}:', L, EddO(IK,L),
EddJ(IK,L),EddN(IK,L);
_OD;
HEdO(L)=HEdN(L);
HEdd(L)=HEdN(L);
_od;
-- pause;
_do L=Nfrus+1,Nfreq;
HEdO(L)=HEdN(Nfrus);
HEdN(L)=HEdN(Nfrus);
HEdd(L)=HEdN(Nfrus); -- to define for primag
write(*,'(a,i5,1p,e12.4)')' added L Hedd:', L,Hedd(L);
_od;
CALL GDEPOS;
NEEDBR=.TRUE.;
-- ENDIF;
IF(NSTEP==0 ! NOUT>=200) THEN;
<*K: PRINT & WRITE FIRST MODEL OF THE RUN *>;
ENDIF;
I=NZON; TAU(I)=0.D0;
_WHILE I>1 & I>=Ncnd _DO
TP=Ty(I); PL=3.D0*DM(I)/(Ry(I)**3-Ry(I-1)**3);
_Do ii=1, Natom [ Yat(ii)=YABUN(ii,I) ];
RADP=.FALSE.;
CALL URSOS; kmhap=I; CALL HAPPA;
TAUste=(Ry(I)-Ry(I-1))*HAPPAL(@LM)*PL;
If(Scat)then;
TAUabs=(Ry(I)-Ry(I-1))*HAPabs(@LM)*PL;
TAUsca=Tauste-Tauabs;
-- Write(@wterm,'(A,1P,3E11.3)')
-- ' tau tauabs tausca',tauste,tauabs,tausca;
If(Tausca>1.d0)then;
Tauste=sqrt(Tauabs*Tausca);
else;
Tauste=Tauabs;
endif;
endif;
TAU(I-1)=TAU(I)+Tauste;
/* If(TAU(I-1)>=Taukli & TAU(I)/TAU(I-1)<Taujmp)Then;
Klim(I)=1;
Else;
Klim(I)=0;
Endif; */
I=I-1
_OD; -- I<Ncnd --& TAU(I)>Taukli
%A_eddtau.K:
<*A: FIND DENSITY IN Rho*>;
tout=t;
_Do I=1,N;
Yout(I)=Y(I,1);
_od
tretard=Ry(NZON-1)/CLIGHT;
t_ob=(t-tretard)*UTIME; -- observer time in sec
CALL PRIBAL(IOUT);
%A_eddtau.KA:
Rho(1)=3.D0*DM(1)/(Ry(1)**3-Rce**3);
_DO I=1,NZON;
-- RADIUS(I)=Ry(I); VELOC(I)=Uy(I); TMPR(I)=Ty(I);
IF(I>1)Rho(I)=3.D0*DM(I)/(Ry(I)**3-Ry(I-1)**3)
_OD;
ELOST=YENTOT(1);
%A_ymax:
_DO IFL=1,NVARS;
_DO I=1,NZON;
YMAX(I+(IFL-1)*NZON)=
MAX(ABS(Y(I+(IFL-1)*NZON,1)),FLOOR(IFL))*Wacc(IFL)
_OD _OD;
--_DO I=NZON*NVARS+1,N; -- for FJ & FH
-- YMAX(I)=MAX(ABS(Y(I,1)),FLOOR(NVARS+1));
--_OD;
_DO IFL=1,Nfrus; -- for FJ
_DO I=1,NZON-Ncnd;
If(Freqmn(IFL)<1.d0)Then;
YMAX(I+Nvars*Nzon+(IFL-1)*(NZON-Ncnd))=
MAX(
ABS(Y(I+Nvars*Nzon+(IFL-1)*(NZON-Ncnd),1)),FLOOR(Nvars+1));
else;
YMAX(I+Nvars*Nzon+(IFL-1)*(NZON-Ncnd))=
MAX(
ABS(Y(I+Nvars*Nzon+(IFL-1)*(NZON-Ncnd),1)),
FLOOR(Nvars+1)/(5.d0*Freqmn(IFL)**3));
endif;
_OD_OD;
_DO I=NZON*NVARS+Krad+1,N; -- for FH
YMAX(I)=MAX(ABS(Y(I,1)),FLOOR(NVARS+1));
_OD;
%A_haptab:
-- only two sets of tables were saved in the RAM for
-- two values of dvdr (given by Stime)
-- now all tables are in hpsavsc
-- 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(tday,hmin)); -- for Nstep=0, t=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; -- , ' tday=',tday;
if(istim^=istold)then;
istold=istim;
_case istim+1 _of
if(istim>1 .and. istim <= @ns)then;
-- Opafile=Opafile(1:len_trim(Opafile)-1)//'2';
ihp=istim;
thaplog1=stmlog(istim-1); thaplog2=stmlog(istim);
else;
write(*,*)' in strad istim=',istim;
stop ' wrong istim in strad';
endif;
_0 ihp=@ns; -- Opafile=Opafile(1:len_trim(Opafile)-1)//'5';
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 ihp=1; -- Opafile=Opafile(1:len_trim(Opafile)-1)//'1';
thaplog1=stmlog(1)-(stmlog(2)-stmlog(1))/2.d0;
thaplog2=stmlog(1)
/*
_2 Opafile=Opafile(1:len_trim(Opafile)-1)//'2';
thaplog1=stmlog(1); thaplog2=stmlog(2)
_3 Opafile=Opafile(1:len_trim(Opafile)-1)//'3';
thaplog1=stmlog(2); thaplog2=stmlog(3)
_4 Opafile=Opafile(1:len_trim(Opafile)-1)//'4';
thaplog1=stmlog(3); thaplog2=stmlog(4)
_5 Opafile=Opafile(1:len_trim(Opafile)-1)//'5';
thaplog1=stmlog(4); thaplog2=stmlog(5)
*/
_esac;
_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);
hpbansc2(L,itp,iro,im)=hpsavsc(L,itp,iro,im,ihp);
_od _od _od _od
-- we do not need read(30) anymore - use hpsavsc
/* 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;-- the first haptab (@ns)
TpTab,RhoTab,hpbansc2;-- the first haptab (@ns)
close(30); */
endif;
%AL:O
-- KM=1;
-- WRITE(@Wres,'(1P,8G10.3)')(@FJ1,L=1,NFRUS);
--
Rho(1)=3.D0*DM(1)/(Ry(1)**3-Rce**3);
_DO I=2,NZON; Rho(I)=3.D0*DM(I)/(Ry(I)**3-Ry(I-1)**3)_OD;
--
<*MaxEr* information on variables, causing max error in STIFF *>;
If ( Givdtl ) Then; -- for test runs
write(18,'(a)')' nstep, Ncnd, Ty(26), Ty(Ncnd+1)';
write(18,'(2I5,1P,4E12.4)')
nstep, Ncnd, Ty(26), Ty(Ncnd+1);
-- nstep,rho(10),rho(11), TY(10), TY(11);
-- Y(NVARS*NZON+(NZON-Ncnd)*(L-1)-Ncnd+KM,1) -- FJ(Km,L)
write(18,'(a,1P,4E12.4)')' k=26 L=40 FJ Bb: ',Y(NVARS*NZON+(NZON-Ncnd)*(40-1)-Ncnd+26,1),
Black(40,Ty(26));
write(18,'(a,1P,4E12.4)')' k=Ncnd+1 L=40 FJ Bb: ',Y(NVARS*NZON+(NZON-Ncnd)*(40-1)+1,1),
Black(40,Ty(Ncnd+1));
Endif;
--
<*A: find integral energy losses: ELvol--volume luminosity,
ELsurf--surface luminosity, Tpsurf & ELtot, calling LOSSEN.
Find lost energy in YEntot(1)<0 for losses, >0 for gains.
Find optical depth from surface TAU & put Ncnd:
TAU(Ncnd)>Fitting_tau. If Ncnd changed then Jstart=0, N-new.
For the version with rezoning:
when number of transparent zones (with Tau < @Tau_floor)
exceeds @Nrezone delete out @Nrezone zones & insert
@Nrezone zones for I>Ncnd *>;
--
<*E* for burning put *E: find ITpmax - number of burning zone
of highest Tp, find in Tpmax(1) predicted Tp,
then: limit step H, checking the rate for Tpmax.
find NREG, add energy DEnuc to YEntot for burnt zones and
put Jstart =0 *>;
<*CORCN* correct convection where flcon < 0 *>;
--
<*O: print model for every NOUT step or in prescribed moments *>;
%AL_MaxEr:
ICW=MaxEr;
If(ICW<=NZON*NVARS)then;
Izon=MOD(ICW-1,NZON)+1;
Ivar=(ICW-1)/NZON +1;
_Case Ivar _of
_1 WRITE(@term,'(2(A,1P,E10.3),2(A,I3))')
' Ry=',Y(MaxEr,1),' Rydot=',Y(MaxEr,2)/H,
' Ivar=',Ivar,' Izon=',Izon;
_2 WRITE(@term,'(2(A,1P,E10.3),2(A,I3))')
' Uy=',Y(MaxEr,1),' Uydot=',Y(MaxEr,2)/H,
' Ivar=',Ivar,' Izon=',Izon;
_3 Pl=Rho(Izon); Tp=Ty(Izon);
_Do ii=1, Natom [ Yat(ii)=YABUN(ii,Izon) ];
RADP=.FALSE.;
CALL URSOS; kmhap=Izon; CALL HAPPA;
WRITE(@term,'(4(A,1P,E10.3),2(A,I3))')
' Ty=',Y(MaxEr,1),' Tydot=',Y(MaxEr,2)/H,
' Rho=',Pl,' Haab6=',Hapabs(6),
' Ivar=',Ivar,' Izon=',Izon;
_esac
else;
ICW=ICW-NZON*NVARS;
If(ICW<=KRAD)Then;
L=(ICW-1)/(NZON-Ncnd)+1;
Izon=Ncnd+MOD(ICW-1,(NZON-Ncnd))+1;
Pl=Rho(Izon); Tp=Ty(Izon);
_Do ii=1, Natom [ Yat(ii)=YABUN(ii,Izon) ];
RADP=.FALSE.;
CALL URSOS; kmhap=Izon; CALL HAPPA;
WRITE(@term,'(4(A,1P,E10.3),2(A,I3))')
' FJ=',Y(MaxEr,1),' FJdot=',Y(MaxEr,2)/H,
' Rho=',Pl,' Haab6=',Hapabs(6),
' in L=',L,' Izon=',Izon;
else;
ICW=ICW-KRAD;
L=(ICW-1)/(NZON-Ncnd)+1;
Izon=Ncnd+MOD(ICW-1,(NZON-Ncnd))+1;
WRITE(@term,'(2(A,1P,E10.3),2(A,I3))')
' FH=',Y(MaxEr,1),' FHdot=',Y(MaxEr,2)/H,
' in L=',L,' Izon=',Izon;
endif;
endif;
km=Izon;
If(Km>1)Then;
U0=Uy(Km-1); RR0=Ry(Km-1);
else;
U0=0.D0; RR0=Rce;
endif;
%ALA:S
<*Loss: find ERLUM & YENTOT *>;
Call NTH(Ncnd); -- find new NTHICK(L), LTHICK(L,KM)
-- _IF Ncnd >= NCNMIN _THEN -- remove this for collapse
<*T: FIND TAU,Ncnd,N *>;
If(KFLAG==0)then;
<*R: REZONE WHEN NUMBER OF TRANSPARENT ZONES
(WITH TAU < @TAU_FLOOR) EXCEEDS @NREZONE *>;
<*F* REGROUPING FREQUENCY TO EMBRACE MAXIMUM FLUX *>;
endif;
-- _FI;
<*Tpcomp: find in Tpcomp the highest Tp(i) for i>Ncnd
for use in approximate treatment of Compton *>
<*N: change NFRUS, KRAD, N if Tpcomp changes
& put JSTART=-JSTART*>;
If(N>NYDIM)then;
WRITE(@Wres,*)'N=',N,' > NYDIM';
STOP 390;
endif;
<*CORRFJ: correct negative values of FJ *>;
-- If(^EDTM &
If((NSTEP-KTAIL>=@N_EDD ! MOD(NSTEP,MBATCH)==0
! MOD(NSTEP,@N_EDD)==0 ) )then;
-- If(t*UTIME > tsmoth) CALL SMOOTH;
<*RT: R-T coeff. in RTphi *>;
if( NSTEP>0)then;
_do L=1,NFRUSO;
_DO IK=Ncnd+1,NZON; -- for EDTM==T
EddO(IK,L)=EddJ(IK,L);
_OD;
HEdO(L)=HEdd(L);
_od;
endif;
Write(@term,*)' Call 1 Feau ';
CALL @FEAUTRIER(Ncnd,Ncnd);
_do L=NFRUSO+1,NFRUS;
_DO IK=Ncnd+1,NZON; -- for EDTM==T
EddO(IK,L)=Eddn(IK,L);
EDDJ(IK,L)=Eddn(IK,L);
_OD;
HEdO(L)=HEdn(L);
HEDD(L)=HEDN(L);
_od;
trlx=@Nrlx*h; @wterm trlx;
CALL GDEPOS;
KTAIL=NSTEP; -- TO TRACK LAST FEAUTRIER
Ycomp=(Boltzk*UTp*Tpcomp)/(Amelec*Cs**2); -- Comptonization par.
Write(@term,*)' Tpcomp Ycomp after Feau: ',Tpcomp,Ycomp;
-- H=1.d-4*H; -- test as coll
-- JSTART=0;
-- Hmax=max(Hmax,t*1.d-3);
@wterm ' 1st Jstart=',Jstart;
NEEDBR=.TRUE.; -- F01BRF/M28Y12 NEEDED
Evalja=.TRUE.;
<*Lim* *>;
endif;
If(JSTART==0)then;
/* Hold=H;
H=MAX(H*1.D-3,1.D1*HMIN);
_Do iY=Nzon+1,2*Nzon;
Y(iY,2)=Y(iY,2)*(H/Hold); -- corrects the acceleration
_od; */
CALL LOSSEN;
YENTOT(2)=H*ELTOT;
@HOLD_BAL=H;
endif;
%ALA_Loss:
RX=1.D0;
IF(HUSED^=@HOLD_BAL)THEN;
_Do J=2,NQused+1;
RX=RX*(HUSED/@HOLD_BAL);
YENTOT(J)=YENTOT(J)*RX
_od;
endif;
_Do J1=1,NQUSED;
_Do J2=J1,NQUSED;
J=NQUSED-J2+J1;
YENTOT(J)=YENTOT(J)+YENTOT(J+1)
_od
_od;
CALL LOSSEN;
ERLUM=HUSED*ELTOT-YENTOT(2);
CALL COSET(METH,max(NQUSED,1),EL,TQ,MAXDER,IDOUB);
_Do J=1,NQused+1; YENTOT(J)=YENTOT(J)+EL(J)*ERLUM _od;
IF(IABS(JSTART)>NQUSED)YENTOT(NQused+1+1)=
ERLUM*EL(NQused+1)/FLOAT(NQused+1);
@HOLD_BAL=HUSED;
%ALAT:
/*Ncnew=NTHICK(1); LTAU=1;
_DO L=2, NFRUS;
IF( NCNEW > NTHICK(L)) THEN;
NCNEW=NTHICK(L);
LTAU=L; -- find the most transparent energy group
ENDIF;
_OD; */
@wterm' Ltau=',Ltau;
I=NZON; TAU(I)=0.D0;
-- _WHILE I>1 & (TAU(I)<=FitTau*TauTol ! I>=Ncnd ) _DO
_WHILE I>1 _DO
TP=Ty(I); Pl=3.D0*DM(I)/(Ry(I)**3-Ry(I-1)**3);
_Do ii=1, Natom [ Yat(ii)=YABUN(ii,I) ];
RADP=.FALSE.;
CALL URSOS;
Rho(I)=Pl;
Press(I)=P*UPI;
kmhap=I; CALL HAPPA;
-- @wterm' after CALL HAPPA kmhap=',kmhap;
Hapmin=Happal(1);
Ltau=1;
_Do L=2,Nfrus;
If(Happal(L)<Hapmin)then;
Hapmin=Happal(L);
Ltau=L;
endif;
_od;
--@wterm' Ltau=',Ltau;
TAUste=(Ry(I)-Ry(I-1))*HAPPAL(LTAU)*PL;
If(Scat)then;
TAUabs=(Ry(I)-Ry(I-1))*HAPabs(LTAU)*PL;
TAUsca=Tauste-Tauabs;
---If(Nstep>7250)then;
----@wterm' Ltau=',Ltau;
----Write(0,'(A,1P,3E11.3)')' tau tauabs tausca',tauste,tauabs,tausca;
---endif;
If(Tausca>1.d0)then;
Tauste=sqrt(Tauabs*Tausca);
else;
Tauste=Tauabs;
endif;
endif;
TAU(I-1)=TAU(I)+Tauste;
I=I-1
_OD; -- I=0 --old I<Ncnd & TAU(I)>FitTau*TauTol
Ncnew=Nthick(Ltau);
-- _Do K=Ncnew-1,I,-1; -- old wrong for small diff. at ncnew & I
_Do K=Ncnew-1,1,-1;
tau(K)=tau(K)-tau(Ncnew); -- to try to catch transparent zones < Ncnew
_od;
-- Ncnew=NZON;
_WHILE TAU(max(Ncnew,1)) <= FitTau & Ncnew>=NCNMIN _DO
-- & TAU(Ncnew-1) <= FitTau*TauTol _DO
Ncnew=Ncnew-1
_OD; -- TAU(Ncnew) > FitTau or Ncnew<NCNMIN or
-- TAU(Ncnew-1) > FitTau*TauTol
-- WRITE(@term,*)' Ncnew=',Ncnew,' Ncnd=',Ncnd,
-- ' TAU(Ncnd or 1)=',TAU(max(Ncnd,1));
/* WRITE(@term,'(2I6,I3,4(A,1P,G10.3))')
NSTEP,MAXER,NQUSED,' H=',H,' Hused=',HUSED,
' Tp=',Ty(1),' Rho=',Rho(1);
WRITE(@term,'(1X,A,1P,G11.3,A,I3,A,G11.2)')
' Time(days)=',tday,' Ncnd=',Ncnd,' TAU(Ncnd or 1)=',TAU(max(Ncnd,1));
*/
WRITE(@term,'(2I6,I3,3(A,1P,G10.3),A,I3,A,G11.2)')
NSTEP,MAXER,NQUSED,' Hu=',HUSED,' Ro=',Rho(1),
' td=',tday,' Ncnd=',Ncnd,' Tau=',TAU(max(Ncnd,1));
IF(Ncnew<NCNMIN) Ncnew=0;
--If(Ncnew<NCNMIN .or. t*Utime>tsmoth) Ncnew=0; -- force Ncnd=0
Ncnew=MIN(Ncnew,@NCOND_MAX);
HOLDFR=.TRUE.; -- forbids the decrease of NFRUS if there is no
-- change of Ncnd or if Ncnd grows (i.e. hot wave
-- moves out)
_If (TAU(max(Ncnd,1))>TauTol*FitTau ! -- Ncnew==0 ! -- for forced Ncnd
TAU(max(Ncnd,1))*TauTol<FitTau) & Ncnew^=Ncnd & (CHNCND ! Nstep==0)_THEN
-- TAU(Ncnd)*TauTol<FitTau) & Ncnew< Ncnd & CHNCND _THEN
_Select
_ Ncnew<Ncnd
[ <*O: ADD FJ,FH FOR ZONES Ncnew+1:Ncnd *>;
HOLDFR=.FALSE. ] -- for Ncnd moving inward
-- the decrease of Nfrus is allowed
_ Ncnew>Ncnd & Ncnew<=@NCOND_MAX
[ <*U: DELETE FJ,FH FOR ZONES Ncnd+1:Ncnew *>]
_ Ncnew>@NCOND_MAX [ KFLAG=-3 ] -- not active now
_end;
If(KFLAG==0)Then;
<*LIM* limiting FH *>;
<*D* test printing*>;
-- test as coll:
JSTART=0;
-- Hmax=max(Hmax,t*1.d-3);
@wterm ' 2nd Jstart=',Jstart;
NEEDBR=.TRUE.; -- F01BRF/M28Y12 NEEDED
-- Evalja=.TRUE.; -- not needed when Jstart=0;
<*H* CHECK FH & PUT SMALLER THAN FJ WITH JSTART=0 *>;
-- XJPH=RPH; -- for limiting DFJ in Dfrad %Derata
endif;
_fi;
-- If(Ncnd>0) Nclast=Ncnd;
%ALATO:
WRITE(@Wres,'(2(A,I4),2(A,1P,E10.2))')
' Ncnew=',Ncnew,' Ncnd=',Ncnd,' TAU(Ncnd)=',TAU(Ncnd),
' TAU(Ncnew)=',TAU(max(Ncnew,1));
WRITE(@term,'(2(A,I4),2(A,1P,E10.2))')
' Ncnew=',Ncnew,' Ncnd=',Ncnd,' TAU(Ncnd)=',TAU(Ncnd),
' TAU(Ncnew)=',TAU(max(Ncnew,1));
_do L=1,Nfrus;
_DO IK=Ncnd+1,NZON; -- for EDTM==T
EddO(IK,L)=EddJ(IK,L);
_OD;
HEdO(L)=HEdd(L);
_od;
Write(@term,*)' Call 2 Feau ';
CALL @FEAUTRIER(Ncnew,Ncnd); -- new position of this call,
-- important for eddi, since it uses old Y(...ncnd), nfrus
LFR=NFRUS;
<*freq: increase NFRUS if needed for Ncnew
new NFRUS is temporarily saved in LFR*>;
KRAD1=(NZON-Ncnew)*LFR;
<*savclean: save FH and FJ in all old zones for all
old frequencies and clean FJ and FH for all
new frequencies for zones K > Ncnd (old) *>;
/* _do L=1,Nfrus;
_DO IK=Ncnd+1,NZON; -- for EDTM==T
EddO(IK,L)=EddJ(IK,L);
_OD;
HEdO(L)=HEdd(L);
_od;
*/
NFRUSO=NFRUS;
NFRUS=LFR; -- will not change if not changed by ALATO_freq
Call NTH(Ncnew);
-- Write(@term,*)' Call 2 Feau ';
-- CALL @FEAUTRIER(Ncnew,Ncnd); -- old position, wrong for eddi
_do L=NFRUSO+1,NFRUS;
_DO IK=Ncnd+1,NZON; -- for EDTM==T
EddO(IK,L)=Eddn(IK,L);
EDDJ(IK,L)=Eddn(IK,L);
_OD;
-- HEdO(L)=HEdn(L); -- nor defined yet
-- HEDD(L)=HEDN(L);
HEdO(L)= HEDD(L);
write(*,'(a,i5,1p,e12.4)')' 2 L Hedd:', L,Hedd(L);
_od;
trlx=@Nrlx*h; @wterm trlx;
@wterm' calling GDepos';
CALL GDEPOS;
KTAIL=NSTEP; -- TO TRACK LAST FEAUTRIER
--
Tp=Ty(Ncnew+1);
If(Ncnew>0)Then;
PL=3.D0*DM(Ncnew+1)/(Ry(Ncnew+1)**3-Ry(Ncnew)**3);
Else;
PL=3.D0*DM(Ncnew+1)/(Ry(Ncnew+1)**3-Rce**3);
Endif;
RADP=.FALSE.; -- NO rad pressure
_Do ii=1, Natom; Yat(ii)=YABUN(ii,Ncnew+1) _od;
CALL URSOS; kmhap=Ncnew+1; CALL OPACIT; CAP2=CAPPA;
--
_DO K=Ncnew+1,Ncnd; -- ADD NEW FJ ,FH
If(K<Nzon)Then;
TP1=Tp;
CAP1=CAP2;
Tp=Ty(K+1);
PL=3.D0*DM(K+1)/(Ry(K+1)**3-Ry(K)**3);
_Do ii=1, Natom; Yat(ii)=YABUN(ii,K+1) _od;
CALL URSOS; kmhap=K+1; CALL OPACIT; -- called by HAPPA ( sometimes called )
call HAPPA;
CAP2=CAPPA;
FLLF1=((Ry(K)*TP1)**4-(Ry(K)*TP)**4)*TP1**4/
(CAP1*(DM(K)+DM(K+1) )*(TP1**4+TP**4));
FLRT1=((Ry(K)*TP1)**4-(Ry(K)*TP)**4)*TP**4/
(CAP2*(DM(K)+DM(K+1) )*(TP1**4+TP**4));
FLUM=CLUM*(FLLF1+FLRT1);
Endif;
If(NSTEP==0 & EddJ(K,@LM)>0.25D0)Then;
dillm=0.75d0-0.25d0*sqrt(12.d0*EddJ(K,@LM)-3.D0);
-- actually another root of square equation for edd<1/3
Tpm=Freqmn(@LM)/LOG(1.D0+dillm/FradJ(K,@LM));
else;
dillm=0.75d0;
TPM=MAX(Ty(K),@TP_FLOOR);
endif;
FLUN=0.D0;
_DO L=1,NFRUS;
IX=NVARS*NZON+(NZON-Ncnew)*(L-1);
-- If(tau(K)>1.d0)then;
Y(IX+(K-Ncnew),1)=min(FradJ(K,L),Black(L,Ty(k))); -- FJ
/*
If(K<=Nthick(L))then;
Y(IX+(K-Ncnew),1)=Black(L,Ty(k));
-- FJ
-- WRITE(@term,'(2(A,1P,E12.5),A,2I3)')
-- ' thick YJ=',Y(IX+(K-Ncnew),1),' Bl=',Black(L,Ty(k)),
-- ' Km L:',K,L
else;
Y(IX+(K-Ncnew),1)=FradJ(K,L);
-- FJ
-- WRITE(@term,'(2(A,1P,E12.5),A,2I3)')
-- ' YJ=',Y(IX+(K-Ncnew),1),' Bl=',Black(L,Ty(k)),
-- ' Km L:',K,L;
endif;
*/
/* IF(Nstep>0)Then;
Y(IX+(K-Ncnew),1)=BLACK(L,TPM);-- FJ
Else;
If(EddJ(K,L)>0.25D0)Then;
dill=0.75d0-0.25d0*sqrt(12.d0*EddJ(K,L)-3.D0);
-- actually another root of square equation for edd<1/3
else;
dill=0.75d0;
endif;
Y(IX+(K-Ncnew),1)=FradJ(K,@LM)*BLACK(L,TPM)*dill
/(BLACK(@LM,TPM)*dillm);
Endif; */
-- IF(K==1 & L==1)
-- WRITE(@Wres,*)'T BL @FJ FRAD=',Ty(K),BLACK(L,TPM),
-- Y(IX+(K-Ncnew),1),FRADJ(K,@LM);
-- insert here *CAPPA/HAPPAL(L)
If(K<Nzon)Then;
If(Nstep>1 & K>Nthick(L))then;
Y(IX+(K-Ncnew)+KRAD1,1)=FradH(K,L); -- FH
elseif(Nstep>1)then;
Y(IX+(K-Ncnew)+KRAD1,1)= -- FH
(FLLF1+FLRT1)*Blackd(L,0.5D0*(Tp+Tp1))*CAP2
/(@Left*(0.5D0*(Tp+Tp1))**4*Ry(K)**2*HAPPAL(L));
else;
-- Y(IX+(K-Ncnew)+KRAD1,1)=0.1d0*FradH(K,L); -- FH
Y(IX+(K-Ncnew)+KRAD1,1)=0.d0; -- let Stella find the flux
endif;
/* If(Nstep>0)Then;
Y(IX+(K-Ncnew)+KRAD1,1)= -- FH
(FLLF1+FLRT1)*Blackd(L,0.5D0*(Tp+Tp1))*CAP2
/(@Left*(0.5D0*(Tp+Tp1))**4*Ry(K)**2*HAPPAL(L));
-- WRITE(@Wres,'(A,1P,E11.3,A,I2)')
-- ' HAP=',HAPPAL(L)/(URHO*UR),' L=',L;
Else; -- Nstep==0
IF( K>NTHICK(L) & NTHICK(L)<Nzon-3 ) THEN;
Y(IX+(K-Ncnew)+KRAD1,1)=
0.1d0*FRADH(K,@LM)*BLACKD(L,Tpm)*dill
/(BLACKD(@LM,Tpm)*dillm); -- FH
-- print *,'K,L,Fradh(k,@LM)',K,L,Fradh(k,@LM);
ELSE;
Y(IX+(K-Ncnew)+KRAD1,1)= -- FH
0.1d0*
(FLLF1+FLRT1)*Blackd(L,0.5D0*(Tp+Tp1))*CAP2
/(@Left*(0.5D0*(Tp+Tp1))**4*Ry(K)**2*HAPPAL(L));
ENDIF;
Endif; */
-- Y(IX+(K-Ncnew)+KRAD1,1)=
-- MIN(Y(IX+(K-Ncnew)+KRAD1,1),Y(IX+(K-Ncnew),1)*HEDD(K));
FLUN=FLUN+Y(IX+(K-Ncnew)+KRAD1,1)*WEIGHT(L);
-- print *,K,L,FLUN;
-- WRITE(@Wres,'(2(A,1P,E12.5),A,2I3)')
--' YH=',Y(IX+(K-Ncnew)+Krad1,1),' FradH=',FradH(K,L),' Km L:',K,L;
Endif;
-- @wterm'K L Nth Bld @FH FRADH=',K,L,
-- Nthick(L),BLACKD(L,Tpm),Y(IX+(K-Ncnew)+Krad1,1),FRADH(K,L);
_OD;
If(K<Nzon)Then;
FLUN=CLUMF*FLUN*Ry(K)**2;
WRITE(@Wres,'(2(A,1P,E12.5),A,I3)')
' Cond Lum=',Flum,' New Lum=',Flun,' Km=',k;
<*CORR* correct new luminosity *>;
Endif;
_OD;
Ncnd=Ncnew; KRAD=KRAD1; N=NZON*NVARS+2*KRAD;
@wterm' NEW Ncnd KRAD N: ',Ncnd,KRAD,N,' AT STEP=',NSTEP;
WRITE(@Wres,*)' NEW Ncnd KRAD N: ',Ncnd,KRAD,N,' AT STEP=',NSTEP;
%ALATO_freq:
-- find Tpcomp for for Ncnew:
Tpcomp=Ty(NZON);
_Do I=Ncnew+1,NZON;
If(Ty(I)>Tpcomp) Tpcomp=Ty(I)
_od;
-- Tpcomp=Ty(Ncnew+Ismax(Nzon-Ncnew,Ty(Ncnew+1),1)); -- CRAY
_WHILE BLACK(MIN(LFR+1,NFREQ),Tpcomp)*FREQMN(MIN(LFR+1,NFREQ))**3 >
FLNOIS _DO
LFR=LFR+1;
_UNTIL LFR>=NFREQ;
LFR=min(LFR,NFREQ);
IF(LFR ^= NFRUS)THEN;
WRITE(@Wres,*)'Tpcomp=',Tpcomp;
WRITE(@Wres,*)' NEW NFRUS -->',LFR,' AT STEP=',NSTEP;
WRITE(@term,*)' NEW NFRUS -->',LFR,' AT STEP=',NSTEP;
ENDIF;
%ALATO_savclean:
_DO L=NFRUS,1,-1; -- SAVE OLD FH
IX=NVARS*NZON+(NZON-Ncnd)*(L-1);
_DO K=NZON-1,Ncnd+1,-1;
Y(IX+(K-Ncnd)+(Ncnd-Ncnew)*L+KRAD1,1)=Y(IX+(K-Ncnd)+KRAD,1);
_OD_OD;
_DO L=NFRUS,1,-1; -- SAVE OLD FJ
IX=NVARS*NZON+(NZON-Ncnd)*(L-1);
_DO K=NZON,Ncnd+1,-1;
Y(IX+(K-Ncnd)+(Ncnd-Ncnew)*L,1)=Y(IX+(K-Ncnd),1);
_OD_OD;
_DO K=Ncnd+1,NZON; -- clean new FJ ,FH
_DO L=NFRUS+1,LFR;
IX=NVARS*NZON+(NZON-Ncnew)*(L-1);
Y(IX+(K-Ncnew),1)=0.D0; -- FJ
Y(IX+(K-Ncnew)+KRAD1,1)= 0.D0; -- FH
_OD;
_OD;
%ALATO_CORR:
If(Nstep>0)Then;
Flum=Flum/Flun;
Flun=0.;
_Do L=1,NFRUS;
IX=NVARS*NZON+(NZON-Ncnew)*(L-1);
Y(IX+(K-Ncnew)+KRAD1,1)=Flum*Y(IX+(K-Ncnew)+KRAD1,1);
FLUN=FLUN+Y(IX+(K-Ncnew)+KRAD1,1)*WEIGHT(L);
_Od;
FLUN=CLUMF*FLUN*Ry(K)**2;
WRITE(@Wres,*)' Corrected New Lum=',Flun;
Endif;
%ALATU:
WRITE(@Wres,'(2(A,I4),2(A,1P,E10.2))')
' Ncnew=',Ncnew,' Ncnd=',Ncnd,' TAU(Ncnd)=',TAU(max(Ncnd,1)),
' TAU(Ncnew)=',TAU(max(Ncnew,1));
WRITE(0,'(2(A,I4),2(A,1P,E10.2))')
' Ncnew=',Ncnew,' Ncnd=',Ncnd,' TAU(Ncnd)=',TAU(max(Ncnd,1)),
' TAU(Ncnew)=',TAU(max(Ncnew,1));
_do L=1,Nfrus;
_DO IK=Ncnd+1,NZON; -- for EDTM==T
EddO(IK,L)=EddJ(IK,L);
_OD;
HEdO(L)=HEdd(L);
_od;
Write(@term,*)' Call 3 Feau ';
CALL @FEAUTRIER(Ncnew,Ncnd); -- new pos.
KRAD1=(NZON-Ncnew)*NFRUS;
_DO L=1,NFRUS; -- DEL FJ
IX=NVARS*NZON+(NZON-Ncnd)*(L-1);
_DO K=Ncnew+1,NZON;
Y(NVARS*NZON+(NZON-Ncnew)*(L-1)+(K-Ncnew),1)=
Y(IX+(K-Ncnd),1);
_OD;
_OD;
_DO L=1,NFRUS; -- DEL FH
-- CHECK FOR K==NZON & FH!!!
IX=NVARS*NZON+(NZON-Ncnd)*(L-1);
_DO K=Ncnew+1,NZON-1;
Y(NVARS*NZON+(NZON-Ncnew)*(L-1)+(K-Ncnew)+KRAD1,1)=
Y(IX+(K-Ncnd)+KRAD,1);
_OD
_OD;
/* old _do L=1,Nfrus;
_DO IK=Ncnd+1,NZON; -- for EDTM==T
EddO(IK,L)=EddJ(IK,L);
_OD;
HEdO(L)=HEdd(L);
_od;
Write(@term,*)' Call 3 Feau ';
CALL @FEAUTRIER(Ncnew,Ncnd); */
trlx=@Nrlx*h; @wterm trlx;
CALL GDEPOS;
KTAIL=NSTEP; -- TO TRACK LAST FEAUTRIER
Ncnd=Ncnew; KRAD=KRAD1; N=NZON*NVARS+2*KRAD;
WRITE(@Wres,*)' NEW Ncnd KRAD N: ',Ncnd,KRAD,N,' AT STEP=',NSTEP;
WRITE(@Wres,*)'Ty(Ncnd+1)=',Ty(Ncnd+1);
%ALAT_LIM:
I=NZON; TAU(I)=0.D0;
_WHILE I>1 & (TAU(I)<=FitTau ! I>=Ncnd) _DO
TP=Ty(I); PL=3.D0*DM(I)/(Ry(I)**3-Ry(I-1)**3);
_Do ii=1, Natom [ Yat(ii)=YABUN(ii,I) ];
RADP=.FALSE.;
CALL URSOS; kmhap=I; CALL HAPPA;
TAU(I-1)=TAU(I)+(Ry(I)-Ry(I-1))*HAPPAL(@LM)*PL;
If(TAU(I-1)>=Taukli & TAU(I)/TAU(I-1)<Taujmp)Then;
Klim(I)=1;
WRITE(@Wres,*)' New Ncnd, Klim=1 for Km=',I;
Else;
Klim(I)=0;
Endif;
I=I-1
_OD; -- I<Ncnd & TAU(I)>FitTau
_DO KT=Ncnd+1,NZON;
TAUZON(KT)=TAU(KT-1)-TAU(KT);
_OD;
%ALAR:
-------_INCLUDE alar;
<*A: FIND IN NTRANS THE NUMBER OF ZONES WITH TAU < @TAU_FLOOR *>;
_IF Nzon>Nzmin & NTRANS >= @NREZONE & T*Utime>TREZON _THEN
--NDELZN=MIN(@NREZONE,(NZON-Ncnd)/2);
NDELZN=NTrans;
WRITE(@Wres,'('' ==> Dropping '',I3,'' zones at STEP:'',I5)')
Ndelzn, Nstep;
WRITE(@term,'('' ==> Dropping '',I3,'' zones at STEP:'',I5)')
Ndelzn, Nstep;
Nznew=max(Nzmin,Nzon-Ndelzn);
<*S* SAVE DELETED ZONES INTO FILE 'DELZON RES',CHANNEL 7
& CHANGE ELOST *>;
<*E: DELETE NDELZN OUTER ZONES *>;
<*D* DEBUGGER FJ,FH *>;
Nzon=Nznew; KRAD=KRAD1; N=NZON*NVARS+2*KRAD;
Call NTH(Ncnd);
PRINT *,' @@@>>> New Nzon Krad N: ',Nzon,KRAD,N;
JSTART=0;
Hmax=max(Hmax,t*1.d-3);
_do L=1,Nfrus;
_DO IK=Ncnd+1,NZON; -- for EDTM==T
EddO(IK,L)=EddJ(IK,L);
_OD;
HEdO(L)=HEdd(L);
_od;
CALL @FEAUTRIER(Ncnd,Ncnd); -- CHANGE EDDJ & HEDD -- will not work for
-- eddi!
trlx=@Nrlx*h; @wterm trlx;
CALL GDEPOS;
NEEDBR=.TRUE.; -- F01BRF/M28Y12 NEEDED
_FI
%ALARD=AS:
%ALARA:
NTRANS=NZON;
_WHILE TAU(NTRANS) < @TAU_FLOOR & NTRANS > Ncnd+1_DO
NTRANS=NTRANS-1
_OD; -- TAU(NTRANS) >= @TAU_FLOOR
NTRANS=NZON-NTRANS;
%ALARS:
WRITE(7,@F0) NSTEP,NDELZN,T*UTIME;
WRITE(7,'(5X,''R'',12X,''V'',12X,''AM'',11X,''DM'',11X,''Rho'')');
WRITE(7,@F1) (Ry(K)*UR,Uy(K)*UR/(UTIME*CRAP),AM(K)*UM,DM(K)*UM,
Rho(K)*URHO,K=NZON,NZON-NDELZN+1,-1);
@F0: Format(' NSTEP:',I5,' NREZONE:',I3,
' STRIPPING TIME(S):',F16.2);
@F1: Format(1X,1P,5E13.5);
ETERM=0.D0;EKIN=0.D0;EGRAV=0.D0;ERADD=0.D0;
ELOST=YENTOT(1);
_DO KM=NZON-NDELZN+1,NZON;
<*E: FIND PL,TP,P,Q,EGAS,ERAD,ACCEL,ACCRAD *>;
IF(KM^=NZON)THEN;
DM2=DM(KM+1);
ELSE;
DM2=DMOUT;
ENDIF;
RX=(DM(KM)+DM2)/2.D0;
EKIN=EKIN+RX*Uy(KM)**2;
EGRAV=EGRAV-RX*AM(KM)/Ry(KM);
ETERM=ETERM+EGAS*DM(KM)*UEI;
ERADD=ERADD+ERAD*DM(KM)*UEI;
_OD;
EKIN=EKIN/2.D0;
ELOST=ELOST+(ETERM+ERADD+EGRAV+EKIN);
YENTOT(1)=ELOST;
%ALARSE:
PL=Rho(KM);TP=Ty(KM);CHEM=CHEM0(KM);
_Do ii=1, Natom [ Yat(ii)=YABUN(ii,KM) ];
-- _IF PL>0.D0 & TP>0.D0_THEN
-- NZR=KM; CALL URSOS;
-- _FI;
IF(KM<=Ncnd)THEN;
ERAD=0.D0; RADP=.TRUE.; CALL URSOS;
ELSE;
RADP=.FALSE.; CALL URSOS;
SUMJ=0.D0;
_DO L=1,NFRUS;
SUMJ=SUMJ+@FJ1*WEIGHT(L);
_OD;
SUMJ=(1.5D1/PI**4)*SUMJ; --TPRAD=SUMJ**.25D0;
ERAD=CSIGM*4.D+13/(PL*CS)*SUMJ;
ENDIF;
%ALARE:
<*R: Drop Ry, Uy, Ty in outer zones *>;
<*F: Compress FH & FJ on new grid *>;
%ALARER:
-- Move Uy,Ty TO NEW positions in Y array
_DO K=1,Nznew; --SAVE OLD Uy,Ty
Y(Nznew+K,1)=Uy(K); -- NR
-- Y(Nznew+K,1)=Uy(K)/sqrt(1.d0-(Uy(K)/clight)**2);
Y(2*Nznew+K,1)=Ty(K);
_OD;
%ALAREF:
KRAD1=(Nznew-Ncnd)*NFRUS;
_DO L=2,NFRUS; -- DEL FJ
IX=NVARS*NZON+(NZON-Ncnd)*(L-1);
_DO K=Ncnd+1,Nznew;
Y(NVARS*Nznew+(Nznew-Ncnd)*(L-1)+(K-Ncnd),1)=
Y(IX+(K-Ncnd),1);
_OD;
_OD;
_DO L=2,NFRUS; -- DEL FH
-- CHECK FOR K==NZON & FH!!!
IX=NVARS*NZON+(NZON-Ncnd)*(L-1);
_DO K=Ncnd+1,Nznew-1;
Y(NVARS*Nznew+(Nznew-Ncnd)*(L-1)+(K-Ncnd)+KRAD1,1)=
Y(IX+(K-Ncnd)+KRAD,1);
_OD
_OD;
%ALA_Tpcomp:
Tpcomp=Ty(NZON);
_Do I=Ncnd+1,NZON;
If(Ty(I)>Tpcomp) Tpcomp=Ty(I)
_od;
-- Tpcomp=Ty(Ncnd+Ismax(Nzon-Ncnd,Ty(Ncnd+1),1)); -- CRAY
%ALAN:
LFR=NFRUS;
/* -- for Watcom
_repeat [
Buff1=FREQMN(LFR)**3; --=
Buff2=BLACK(LFR,Tpcomp); --=
Buff3=Buff1*Buff2; --=
_IF (Buff3 >= FLNOIS/@FREQ_TOL .OR. HOLDFR) _Then; _leave; _FI;
LFR=LFR-1;
]
*/
_WHILE BLACK(LFR,Tpcomp)*FREQMN(LFR)**3 < FLNOIS/@FREQ_TOL
& ^HOLDFR & Lfr>NFRMIN _DO
LFR=LFR-1
_OD;
IF(LFR < NFRMIN) LFR=NFRMIN;
IF(LFR ^= NFRUS)THEN;
-- WRITE(@Wres,*)'Tpcomp=',Tpcomp;
NFRUS=LFR;
KRAD1=(NZON-Ncnd)*NFRUS;
_DO L=1,NFRUS;
IX=NVARS*NZON+(NZON-Ncnd)*(L-1);
_DO K=Ncnd+1,NZON-1; -- SAVE FH
Y(IX+(K-Ncnd)+KRAD1,1)=Y(IX+(K-Ncnd)+KRAD,1);
_OD
_OD;
KRAD=KRAD1;
N=NZON*NVARS+2*KRAD;
JSTART=0;
Hmax=max(Hmax,t*1.d-3);
NEEDBR=.TRUE.; -- F01BRF/M28Y12 NEEDED
WRITE(@Wres,*)' NEW NFRUS KRAD N-->',NFRUS,KRAD,N,' AT STEP=',NSTEP;
Call NTH(Ncnd);
KTAIL=-@N_EDD; -- to force call FEAU
<*B* TEST PRINTING *>;
ELSEIF(LFR<NFREQ)THEN; -- LFR==NFRUS
/* -- for Watcom
_repeat [
Buff1=FREQMN(LFR+1)**3; --=
Buff2=BLACK(LFR+1,Tpcomp); --=
Buff3=Buff1*Buff2; --=
_IF Buff3 <= FLNOIS _Then; _leave; _FI;
LFR=LFR+1;
_IF LFR==NFREQ _THEN _leave; _FI;
]
*/
_WHILE BLACK(LFR+1,Tpcomp)*FREQMN(LFR+1)**3 > FLNOIS _DO
LFR=LFR+1;
_UNTIL LFR==NFREQ;
IF(LFR ^= NFRUS)THEN;
WRITE(@Wres,*)'Tpcomp=',Tpcomp;
KRAD1=(NZON-Ncnd)*LFR;
<*A: DEFINE Y FOR NEW FREQUENCIES *>;
NFRUSO=NFRUS;
NFRUS=LFR;
KRAD=KRAD1;
N=NZON*NVARS+2*KRAD;
JSTART=0;
Hmax=max(Hmax,t*1.d-3);
NEEDBR=.TRUE.; -- F01BRF/M28Y12 NEEDED
WRITE(@Wres,*)' NEW NFRUS KRAD N-->',NFRUS,KRAD,N,' AT STEP=',NSTEP;
WRITE(@term,*)' NEW NFRUS KRAD N-->',NFRUS,KRAD,N,' AT STEP=',NSTEP;
Call NTH(Ncnd);
KTAIL=-@N_EDD; -- to force call FEAU
<*D* TEST PRINTING *>;
ENDIF;
ENDIF;
%ALANA:
_DO L=NFRUS,1,-1; -- SAVE OLD FH
IX=NVARS*NZON+(NZON-Ncnd)*(L-1);
_DO K=NZON-1,Ncnd+1,-1;
Y(IX+(K-Ncnd)+KRAD1,1)=Y(IX+(K-Ncnd)+KRAD,1);
_OD_OD;
_DO K=Ncnd+1,NZON; -- ADD NEW FJ ,FH
/* CM1=CM2; TT1=TP;
If (K<NZON) Then;
TP=Ty(K+1);PL=3.D0*DM(K+1)/(Ry(K+1)**3-Ry(K)**3);
RADP=.FALSE.;
_Do i=1, Natom [ Yat(i)=YABUN(I,K+1) ];
CALL URSOS; kmhap=K+1;CALL OPACIT;
CM2=CAPPA*DM(K+1);
FL1=((Ry(K)*TT1)**4-(Ry(K)*Tp)**4)/(CM1+CM2);
Endif; */
_DO L=NFRUS+1,LFR;
IX=NVARS*NZON+(NZON-Ncnd)*(L-1);
-- Y(IX+(K-Ncnd),1)=BLACK(L,Tpcomp);
Y(IX+(K-Ncnd),1)=Y(IX+(K-Ncnd)-(NZON-Ncnd),1) -- FJ
*EXP((FREQMN(L-1)-FREQMN(L))/Tpcomp);
-- *EXP((FREQMN(L-1)-FREQMN(L))/Ty(K));
Y(IX+(K-Ncnd)+KRAD1,1)= 0.D0; -- FH
-- IF(K<NZON)
-- Y(IX+K+KRAD1,1)=Y(IX+K+KRAD1,1) -- FH
-- *EXP((FREQMN(NFREQ-1)-FREQMN(NFREQ))/Ty(Ncnd))
_OD;
_OD;
%ALA_CORRFJ:
_Do K=NZON*NVARS+1,NZON*NVARS+KRAD;
Y(K,1)=max(0.D0,Y(K,1)) -- FJ>=0
_od;
%ALA_RT:
Ryzer=Ry(Kmcor);
-- RTphi(Nzon)=0.d0; -- try to use non-zero for the last zone also
RTphi(Nzon+1)=0.d0;
-- must be done for every change of Nzon
EpsUq=0.d0;
-- _Do K=1,NZON-1;
_Do K=1,NZON; -- try to use non-zero for the last zone also
EpsUq=max(EpsUq,DRT*Eps*abs(Uy(K)));
-- we continue to apply nonzero QRT_cold even
-- for positive UU2 since we have no infinite
-- accuracy: Eps is the measure of the Stiff accuracy
TauRT=(Amout-Am(K))*(0.4d0*Urho*UR)/Ry(K)**2; -- appr. tau
RTphi(K)=BQ*dM(K)*Ry(K)**(NRT-1)/(1.d0+0.033d0*TauRT); -- cold
_od;
%ALAF:
_IF Ty(Ncnd+1)*@FREQ_TOL > FREQ(NFREQ+1) _THEN
@wterm ' >>>>> FREQUENCIES INCREASED AT STEP',NSTEP;
<*A: REDEFINE FREQ. GROUPS INCREASING FREQ(I) *>;
JSTART=0;
Hmax=max(Hmax,t*1.d-3);
NEEDBR=.TRUE.; -- F01BRF/M28Y12 NEEDED
_ELSE _IF Ty(Ncnd+1)*@FREQ_TOL < FREQ(NFREQ-1) _THEN
@wterm ' <<<<< FREQUENCIES DECREASED AT STEP',NSTEP;
<*E: REDEFINE FREQ. GROUPS DECREASING FREQ(I) *>;
JSTART=0;
Hmax=max(Hmax,t*1.d-3);
NEEDBR=.TRUE.; -- F01BRF/M28Y12 NEEDED
_FI_FI;
%ALAFA:
_DO L=1,NFREQ;
FREQ(L)=FREQ(L+1)
_OD;
FREQ(NFREQ+1)=FREQ(NFREQ)*FREQ(2)/FREQ(1);
W1OLD=WEIGHT(1); -- SAVE OLD W1
<*B: DEFINE FREQMN,DLOGNU,WEIGHT *>;
IX=NVARS*NZON+(NZON-Ncnd)*(L-1)-Ncnd;
_DO K=Ncnd+1,NZON;
FJ1OLD(K)=Y(NVARS*NZON-Ncnd+K,1); -- SAVE FJOLD
FH1OLD(K)=Y(NVARS*NZON-Ncnd+K+KRAD,1); -- SAVE FHOLD
_OD;
_DO L=1,NFREQ-1;
IX=NVARS*NZON+(NZON-Ncnd)*(L-1)-Ncnd;
_DO K=Ncnd+1,NZON;
Y(IX+K,1)=Y(IX+(NZON-Ncnd)+K,1); -- FJ
IF(K<NZON)
Y(IX+K+KRAD,1)=Y(IX+(NZON-Ncnd)+K+KRAD,1); -- FH
_OD;
_OD;
IX=NVARS*NZON+(NZON-Ncnd)*(NFREQ-1)-Ncnd;
_DO K=Ncnd+1,NZON;
Y(IX+K,1)=FJ1OLD(K)*W1OLD/WEIGHT(NFREQ); -- FJ
-- Y(IX+K,1)=Y(IX+K,1) -- FJ
-- *EXP((FREQMN(NFREQ-1)-FREQMN(NFREQ))/Ty(Ncnd));
IF(K<NZON)
-- Y(IX+K+KRAD,1)=Y(IX+K+KRAD,1) -- FH
-- *EXP((FREQMN(NFREQ-1)-FREQMN(NFREQ))/Ty(Ncnd))
-- Y(IX+K+KRAD,1)=0.D0 -- FH
Y(IX+K+KRAD,1)=FH1OLD(K)*W1OLD/WEIGHT(NFREQ); -- FH
_OD;
<*D: TEST PRINTING *>;
%ALAFAB:
_DO L=1,NFREQ;
FREQMN(L)=SQRT(FREQ(L)*FREQ(L+1));
@wterm cs/(freqmn(l)*UFreq*1.d8);
_OD;
_DO L=1,NFREQ;
-- IF(L==1)THEN;
-- WL=(FREQ(L+1)-FREQ(L))*FREQMN(L)**3;
-- DLOGNU(L)=1.D0/LOG(FREQMN(L)/FREQ(L));
-- ELSE;
WEIGHT(L)=(FREQ(L+1)-FREQ(L))*FREQMN(L)**3;
DLOGNU(L)=1.D0/LOG(FREQMN(L)/FREQMN(L-1));
-- ENDIF;
_OD;
TSS=0.D0;
_DO L=1,NFREQ-1;
TSS=TSS+BLACKD(L,.5D0*(Ty(Ncnd)+Ty(Ncnd+1)))*WEIGHT(L);
_OD;
WL=WEIGHT(NFREQ);
-- insert here *CAPPA/HAPPAL(L)
WEIGHT(NFREQ)=(4.D0*(PI*0.5D0*(Ty(Ncnd)+Ty(Ncnd+1)))**4/15.D0-TSS)
/BLACKD(NFREQ,.5D0*(Ty(Ncnd)+Ty(Ncnd+1)));
TSS=TSS+BLACKD(NFREQ,.5D0*(Ty(Ncnd)+Ty(Ncnd+1)))*WEIGHT(NFREQ);
PRINT *,' TSS=',TSS,
4.D0*(PI*0.5D0*(Ty(Ncnd)+Ty(Ncnd+1)))**4/15.D0;
PRINT *,' WL=',WL,' WEIGHT(NFREQ)=',WEIGHT(NFREQ);
@wterm ' FREQ FREQMN';
PRINT'(1X,1P,10G13.3)',FREQ;
PRINT'(1X,1P,10G13.3)',FREQMN;
%ALAFE:
_DO L=NFREQ,1,-1;
FREQ(L+1)=FREQ(L)
_OD;
FREQ(1)=FREQ(1)*FREQ(2)/FREQ(3);
<*B: DEFINE FREQMN,DLOGNU,WEIGHT *>;
_DO L=NFREQ-1,1,-1;
IX=NVARS*NZON+(NZON-Ncnd)*(L-1)-Ncnd;
_DO K=Ncnd+1,NZON;
Y(IX+(NZON-Ncnd)+K,1)=Y(IX+K,1); -- FJ
IF(K<NZON) Y(IX+(NZON-Ncnd)+K+KRAD,1)=Y(IX+K+KRAD,1); -- FH
_OD;
_OD;
IX=NVARS*NZON-Ncnd; -- L==1
_DO K=Ncnd+1,NZON;
Y(IX+K,1)=Y(IX+K,1)*FREQMN(2)/FREQMN(1); -- FJ
IF(K<NZON) Y(IX+K+KRAD,1)=Y(IX+K+KRAD,1)*FREQMN(2)/FREQMN(1);
_OD;
<*D: TEST PRINTING *>;
%ALAFEB=ALAFAB:
%ALAFAD=AS:
%ALATD=AS:
%ALANB=AS:
%ALAND=AS:
%ALAFED=AS:
/* now separate node: %IT=ALAT:
*/
%I_RT=ALA_RT:
%I_Tpcomp=ALA_Tpcomp:
%IN=ALAN:
%ALA_Lim:
I=NZON; TAU(I)=0.D0;
_WHILE I>1 & I>=Ncnd _DO
TP=Ty(I); PL=3.D0*DM(I)/(Ry(I)**3-Ry(I-1)**3);
_Do ii=1, Natom; Yat(ii)=YABUN(ii,I) _od;
RADP=.FALSE.;
CALL URSOS; kmhap=I; CALL HAPPA;
TAU(I-1)=TAU(I)+(Ry(I)-Ry(I-1))*HAPPAL(@LM)*PL;
If(TAU(I-1)>=Taukli & TAU(I)/TAU(I-1)<Taujmp)Then;
Klim(I)=1;
WRITE(@Wres,*)' New Feau, Klim=1 for Km=',I;
Else;
Klim(I)=0;
Endif;
I=I-1
_OD; -- I<Ncnd --& TAU(I)>Taukli
%ALE:N
<*F: FIND ITPMAX *>;
_IF JSTART>0_THEN
_DO J1=1,JSTART+1;
TPMAX(J1)=Y(2*NZON+ITPMAX,J1)
_OD;
_DO J1=1,JSTART;
_DO J2=J1,JSTART;
J=JSTART-J2+J1;
TPMAX(J)=TPMAX(J)+TPMAX(J+1)
_OD
_OD;
_ELSE
TPMAX(1)=Y(2*NZON+ITPMAX,1)
_FI;
<*R: CONTROL ZONE ITPMAX & ADD ENERGY DENUC FOR BURNT ZONES *>;
191:FORMAT(' TP IN ZONE',I3,' FAILS TO CONVERGE IN',I3,' ITERATIONS');
%ALEF:
ITPMAX=NZON;
_DO I=1,NZON;
IF(Ty(I)>Ty(ITPMAX) & NREG(I)==1) ITPMAX=I
_OD;
%ALER:
_DO I=1,NZON;
GOTO(101,102,103),NREG(I);
101: Continue;
_IF I==ITPMAX_THEN
<*B: FOR I==ITPMAX FIND ET, TMGROW==TP/TPDOT,
FOR TMGROW*EPS <H: H=EPS*TMGROW *>;
IF(TMGROW>=TMCRIT&1.D0>=YCARB*YDOT*TMCRIT*UTIME)GOTO 103;
JSTART=0; NPRIN=NOUT;
IBURNT=I;
XCARB=ACARB*YCARB;
PRINT'('' ZONE'',I4,'' BURNT AT NSTEP='',I5,'' TIME='',G15.10,
'' SEC HUSED='',E8.3)',I,NSTEP,T*UTIME,HUSED*UTIME;
YCAR1=YCARB;
_DO K=MAX0(1,I-1),MIN0(I+1,NZON);
-- RADIUS(K)=Ry(K); VELOC(K)=Uy(K); TMPR(K)=Ty(K)
_OD;
CALL PRIBUR(IBURNT); NREG(I)=2; TP=Ty(I); PL=Rho(I); YCARB=YCAR1;
<*E: FIND NEW Ty(I) & DENUC - ENERGY RELEASED BY 12C - 28SI*>;
PRINT'('' TP 28SI ='',E11.5,'' Q='',E11.5)',Ty(I),QCSI*YCARB;
_ELSE GOTO 103_FI;
YENTOT(1)=YENTOT(1)+DENUC;
102: IF(Ty(I)<TPNSE)GOTO 103;
NREG(I)=3; JSTART=0; NPRIN=NOUT;
IBURNT=I;
<*I: =\ASEDE, BUT FOR 28SI - 56NI *>;
YENTOT(1)=YENTOT(1)+DENUC;
PRINT'('' TP 56NI ='',E11.5,'' Q='',E11.5)',Ty(I),QSINI;
103:TMGROW=1.D0<* TMGROW GT TMCRIT IS NEEDED *>;
_OD;
%ALERB:
TP=TPMAX(1); PL=Rho(I); CALL BURNC;
YCARB=1.D0/YAINV(I);
AS=ASI/(1.D0+(ASI-ACARB)*YCARB); ZS=AS*YELECT;
-- _Do i=1, Natom [ Yat(i)=YABUN(I,KM+1) ];
CALL URSOS; <* HENCE ET *>
TMGROW=(TP*ET*UEI)/(CK2*QNUC*YDOT*YCARB**2);
_IF H>(EPGROW*TMGROW)_THEN
H=(EPGROW*TMGROW);JSTART=-JSTART
_FI;
%ALERE:
AS=ASI/(1.D0+(ASI-ACARB)*YCARB);
ZS=YELECT*AS;
-- Xurs=XArr(Km+1); Yurs=YArr(Km+1); Zurs=ZArr(Km+1);
CALL URSOS;
EGAS0=EGAS+QCSI*YCARB;
AS=ASI; ZS=ZSI;
<*T: FIND NEW TP BY NEWTON ITERATIONS *>;
DENUC=QCSI*YCARB*DM(I)*UEI;
%ALERET:
DELTP=(EGAS0-EGAS)/ET; ITER=0; TP=TP+DELTP;
_WHILE ABS(DELTP)>EIT&ITER<=10_DO
-- Xurs=XArr(Km+1); Yurs=YArr(Km+1); Zurs=ZArr(Km+1);
CALL URSOS;
DELTP=(EGAS0-EGAS)/ET; TP=TP+DELTP;
IF(TP<0.D0) TP=(TP-DELTP)/2.D0;
ITER=ITER+1
_OD;
_IF ITER>10_THEN
PRINT 191,I,ITER;
_ELSE Ty(I)=TP_FI;
%ALERI:
AS=ASI; ZS=ZSI;
-- Xurs=XArr(Km+1); Yurs=YArr(Km+1); Zurs=ZArr(Km+1);
CALL URSOS;
EGAS0=EGAS+QSINI; AS=ANI; ZS=ZNI;
<*T: FIND NEW TP BY NEWTON ITERATIONS *>;
DENUC=QSINI*DM(I)*UEI
%ALERIT=ALERET:
%AL_CorCn:
_DO IK=1,NZON;
IF(FLCON(IK)<0.D0 ! UC(IK)<0.D0)THEN;
_DO IKK=1,MAXDER+1;Y(IK+(NVARS-2)*NZON,IKK)=0.D0_OD;
ENDIF
_OD
%ALO:
/* TOO - moment of prescribed output in secs,
KO - number of prescribed output,
TO(1:KOMAX) - beginning of prescribed output series,
STO(1:KOMAX) - steps after TO,
NTO(1:KOMAX) - how much outputs with step STO, if not zero,
KNTO - number of output done with step STO.
initially TOO=TO(1); KO=1; KNTO=1 */
-- CORRECT below for continuation runs!!! First too's must be
-- ignored
tretard=Ry(NZON-1)/CLIGHT;
t_ob=(t-tretard)*UTIME; -- observer time in sec
_WHILE t_ob >= TOO & KO<=KOMAX & Nstep > 10 _DO
<*F: INTERPOLATE FOR TIME
t = TOO/UTIME + Ry(NZON-1)/CLIGHT
AND OUTPUT MODEL *>;
if( NTO(KO)^=0 )then;
if( KNTO<=NTO(KO) )then;
TOO=TOO+STO(KO);
KNTO=KNTO+1;
else;
KNTO=1;
KO=KO+1;
IF(KO<=KOMAX) TOO=TO(KO); -- here TO in secs
endif;
else;
KO=KO+1;
IF(KO<=KOMAX) TOO=TO(KO); -- here TO in secs
endif;
_OD;
-- call beams; -- now in eddirel.trf
ELOST=YENTOT(1);
if( MOD(NSTEP,NOUT)==0 )then;
@DateTime;
WRITE(@Wres,'(1X,'' ====> STEP:'',I6,2X,A,2X,'' ON '',A)')
NSTEP,@TIME,@DATE;
tout=t;
_Do I=1,N;
Yout(I)=Y(I,1);
_od
CALL PRIBAL(IOUT); -- PRINT & BALANCES
--if(Nstep>0)
CALL PRIMAG;
-- PRIMAG IS not CALLED BY PRIBAL now
else;
-- @wterm' calling primag';
if(Nstep>1)CALL PRIMAG;
-- @wterm' after primag';
-- if(mod(Nstep,Nout/10)==0) then;
-- call pribal(+100); -- for swdet
-- endif;
-- If(Givdtl)then;
-- CALL PRIBAL(50);
-- <*Details: *>;
-- endif;
endif;
--IF(JSTART==0 & ((ITPMAX-1)/10)*10==(ITPMAX-1)) CALL PRIBAL(IOUT);
<*H* CHECK FH & PUT SMALLER THAN FJ WITH JSTART=0 *>;
%ALOF:
--@wterm ' too/UTIME tretard t ', too/UTIME, tretard, t ;
--@wterm ' TOO/UTIME+tretard-t H ', TOO/UTIME+tretard-t, H;
tout=TOO/UTIME+tretard;
_DO I=1,N;
Yout(I)=0.D0;
RX=1.D0;
_DO K=1,JSTART+1;
Yout(I)=Y(I,K)*RX+Yout(I); --interpolation using predictor
-- RX=RX*((TOO/UTIME-t)/H)
RX=RX*((TOO/UTIME+tretard-t)/H)
_OD
_OD;
-- NOW WE HAVE TO EXTRACT Ry, TEMPERATURE AND VELOCITY FROM THE BUFFER Yout
_DO I=1,NZON;
Ry(i)=Yout(i);
Uy(I)=Yout(NZON+I); -- NR
-- Uy(I)=Yout(NZON+I)/sqrt(1.d0+(Yout(Nzon+K,1)/clight)**2
Ty(I)=Yout(2*NZON+I);
-- @wterm ' I Ty =', Ty(I);
_OD;
RX=1.D0; ELOST=0.D0;
_DO K=1,JSTART+1;
ELOST=YENTOT(K)*RX+ELOST;
RX=RX*((TOO/UTIME-T)/H)
_OD;
<*CALCULATE DENSITY IN Rho*>
Rho(1)=3.D0*DM(1)/(Ry(1)**3-Rce**3);
_DO I=2,NZON;
Rho(I)=3.D0*DM(I)/(Ry(I)**3-Ry(I-1)**3);
@wterm ' I rho =', I, Rho(I);
_OD;
@wterm '('' ***TIME='',G9.3,'' SEC'')',TOO;
CALL PRIBAL(IOUT);
CALL PRIBAL(+100); -- for swdet
%ALO_DETAILS:
WRITE(18,*) ' NSTEP tobs Mass Rbeg E';
WRITE(18,*)
NSTEP, (T-Ry(NZON-1)/CLIGHT)*(UTIME*CRAP)/86400.D0,
AMOUT*UM,RADBEG,EKO+Eburst;
WRITE(18,*)' Km Lum accsel QRT UU ';
WRITE(18,*)Nzon;
_Do Km=1,NZON;
WRITE(18,*) Real(Km),ArrLum(Km),Acc(Km),QRTarr(Km),UUarr(Km);
_od;
%ALOH:
-- WRITE(5,*) ' OLD JSTART=',JSTART;
IF(Ncnd==0)THEN;
_DO KM=Ncnd+1, NZON-1;
_DO L=1,NFRUS;
IF(@FH1<0.D0)THEN;
IF(ABS(@FH1)>ABS(@FJ2)*@TOLF & ABS(@FH1)>FLOOR(NVARS+1))
THEN;
-- WRITE(5,'(1X,A,1P,2E11.3,3I6)')
-- ' FH1>FJ2 NSTEP KM L ',@FH1,@FJ2,NSTEP,KM,L;
@FH1=-ABS(@FJ2)*EDDH(KM) ;
-- JSTART=0;
ENDIF;
ELSE;
IF(ABS(@FH1)>ABS(@FJ1)*@TOLF & ABS(@FH1)>FLOOR(NVARS+1))
THEN;
-- WRITE(5,'(1X,A,1P,2E11.3,3I6)')
-- ' FH1>FJ1 NSTEP KM L ',@FH1,@FJ1,NSTEP,KM,L;
@FH1=ABS(@FJ1)*EDDH(KM) ;
-- JSTART=0;
ENDIF;
ENDIF;
_OD;
_OD;
ENDIF;
%ALATH=ALOH:
%A_HLIM=ALOH:
%_usemodules:
/* use chem_composition_types, only: chem_composition_type, chem_composition_create, chem_composition_retain
, chem_composition_release, write_chem_composition, chem_composition_set
, p_method_saha, p_method_lucy;
use data_chemical, only: p_sampleZ;
*/
/* use rad_intensity_types, only: rad_intensity_type, rad_intensity_create, rad_intensity_set
, rad_intensity_release, write_rad_intensity;
use stl_ursos_lucy, only: ursos_bridge;
use stl_units, only: sub_unitsStella2World, sub_unitsWorld2Stella, p_ulgpl, p_ulgr ;
use td_types, only: td_type, td_type_create
, td_type_set, td_type_release
, write_td_type;
*/
use f77_interface, only: init_stl, finalize_stl;
use kinds, only: dp;
use timings, only: trace_debug, timeset, timestop,timeprint;
use cp_para_types, only: cp_para_env_type;
use cp_para_env, only: cp_para_env_create, cp_para_env_release,
cp_para_env_retain, cp_para_env_write;
use message_passing, only: mp_sync,
mpi_comm_world, add_mp_perf_env, get_mp_perf_env, mp_max,
mp_perf_env_release, mp_perf_env_retain, mp_perf_env_type,
mp_world_finalize, mp_world_init, rm_mp_perf_env;
/*_C
#include "cp_common_uses.h"
_trefor*/
use cp_log_handling , only: cp_failure_level;
use cp_error_handling, only: cp_error_type, cp_a_l, cp_failure_level;
English
Русский