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