radiative hydro
code stella
-- _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;
-- version for calling ronfict directly _DEFINE @SAVE_MODEL SAVEM @TAUPH 0.64D0 @HOLD_BAL HOLDBL -- OLD STEP FOR BALANCE @ReadData "READ(Nc,'(A)'); Read(Nc,*)" @U (CCL/3.65D+03) @B (CCL/4.4D+03) @V (CCL/5.5D+03) @R (CCL/7.D+03) @I (CCL/9.D+03) @Lyman (CCL/912.D0) ; _LABEL @F1,@F2; SUBROUTINE BEGIN; use mLbalsw, only: lossen; --, pribal, primag; -- _TRACE "@wterm' Nzon=',Nzon,' Mzon=',Mzon," -- _TRACE "@wterm' Ncnd=',Ncnd,' Krad=',Krad," -- _TRACE "@wterm' Natom=',Natom,' Zn(Natom)=',Zn(Natom)," IMPLICIT REAL*8(A-H,O-Z); _INCLUDE snrad; _Include opacity; _INCLUDE abo; -- real*8 hptrbab(Nfreq,@Ntab,@Ntab,Mzon/@skip), -- hptrbsc(Nfreq,@Ntab,@Ntab,Mzon/@skip); <*B: vars, PARAMETERS, data *>; <*R: READ INPUT DATA FOR BURST *>; RETURN; 9:WRITE(@Wres,*)' BEGIN: error in read file ',Lunit; stop; END; <*S: SAVE CURRENT MODEL *>; <*I: BLOCK DATA *>; %B:O CHARACTER*80 STRING; -- Dimension dumFreq(Nfreq+1),dumFreqmn(Nfreq),dumwavel(Nfreq); -- now in opacity.inc Dimension indfr(6), indop(6); -- keys for freq and opacit -- depending on Knadap -- PARAMETER(UFRC=2.417965D+14/1.16048D+04); CHARACTER*3 status,NFILE*80; Logical petread; -- is Peters's file already read? Logical fexist; Data petread/.false./; Data indfr/1,2,2,3,3,3/, -- 1 geom., 2 peter, 3 ronfict indop/1,1,2,3,3,3/; -- 1 opazr, 2 opahoef,hap2int, 3 hapintron -- DIMENSION WAVEL(NFREQ);-- WAVELENGTHS IN ANGSTREMS -- Now WAVEL in Include Opacity -- DATA WLMAX/5.0D+04/,WLMIN/1.D+01/; DATA WLMAX/5.0D+04/,WLMIN/1.D+00/; -- wavelength range in Angstrems for SN1987A -- DATA WLMAX/5.0D+04/,WLMIN/9.D+01/; -- for SN2L %R: <*test* READ START_PACK DATA FOR RUN , CHANNEL - 8 *>; <*D: READ START_PACK DATA FOR RUN , CHANNEL - 8 *>; <*M: READ MODEL CHANNEL - 9 -- begin run CHANNEL -10 -- continue run READ Lcurve CHANNEL -11 *>; <*U: UNITS & CONSTANTS FOR DIMENSIONLESS EQUATIONS *>; inquire(file='lines1.out',exist=fexist); if(fexist)then; open(file='lines1.out ', unit=7, status='old'); else; open(file='lines1.out ', unit=7, status='new'); endif; If (ABS(Knadap) .ne. 1) then; -- see begrad for Knadap <*haptab: read tables of happa if needed *>; endif; AS=ACARB; ZS=ZCARB; N=NZON*NVARS+2*KRAD; HMIN=HMIN/UTIME; HMAX=HMAX/UTIME; METH=METHN; Haplim=1.D0/(3.D0*FitTau); Uplim=1.D0+Haplim; <*H: PRINT HEADER OF THE RUN *>; -- For enhanced quadratic viscosity: -- Initially DRT in STRAD DATA defines the mass "distance" -- in Solar mass units -- from AM(NCND) with enhanced pseudoviscosity determined by BQ. -- Then it is used in DFRAD & LBALRAD in the following form: -- -- DRT=(UM/DRT)**2; -- initial DRT is destroyed here, be careful!! -- for R-T do not change DRT ! -- -- For linear R-T viscosity (Hot or Cold) -- DRT was the optical thickness (assuming Hydrogen & full -- ionization where QRT begins to act) -- later the value 30. was fixed & now -- DRT determines the enhancement of EpsUq -- NRT may be noninteger in RTphi - see STRAD IF(NSTEP==0)THEN; H=1.D+05*HMIN; @HOLD_BAL=H; -- INITIAL CONDITIONS FOR BALANCES YENTOT(1)=ELOST;YENTOT(2)=H*ELTOT; <*F: PREPARE FREQUENCY ARRAYS etc. for Ronfict *>; NFRUS=NFREQ; _Do L=1,NFRUS; NTHICK(L)=ncnd; _od; _Do KM=NCND+1,NZON; _Do L=1,NFRUS; LTHICK(L,KM)=.FALSE.; _od; _od; ENDIF; <*Freqob: PREPARE FREQUENCY ARRAY freqob from freqmn *>; toldg=-1.d0; WRITE(@Wres, '('' FREQ:'',1P,10E12.5)')FREQ; WRITE(@Wres, '('' FREQMN:'',1P,10E12.5)')FREQMN; WRITE(@Wres, '('' WAVE BOUNDS:'',1P,11E10.3)') (CCL/FREQ(LW),LW=1,NFREQ+1); WRITE(@Wres, '('' WAVES:'',1P,10E12.5)')(CCL/FREQMN(LW),LW=1,NFREQ); WRITE(@Wres, '('' WEIGHT:'',1P,10E12.5)')WEIGHT; <*BANDS: find numbers of freq. groups for UBV *>; LST=2;KENTR=0; -- RGASA=-.831434D0*(1.D0-ACARB/ASI); -- IBURNT=0; -- INDICATES PRINTING IN BURNT ZONE RHO(1)=3.D0*DM(1)/(Ry(1)**3-Rce**3); _DO I=2,NZON; -- write(@term,'(a,i3,2(a,1p,e12.4))')' i=',i,' dm=',dm(i), -- ' ry=',Ry(i); RHO(I)=3.D0*DM(I)/(Ry(I)**3-Ry(I-1)**3) _OD; IF(NSTEP==0)THEN; CALL LOSSEN; @HOLD_BAL=H; -- INITIAL CONDITIONS FOR BALANCES YENTOT(1)=ELOST;YENTOT(2)=H*ELTOT; ENDIF; chem=0.5d0; -- to avoid too low Xe in saha <*testhappa* *>; %R_testhappa: pl=1.d-8; K=48; _do itp=1,1; Tp=1.d+01**(4d0+(itp-1))/UTp; -- write(@term,'(a,i5,1p2e10.2)') ' Test happal for zone, pl*Urho, Tp*UTp', -- K,pl*Urho,Tp*UTp; -- _Do i=1, Natom; Yat(i)=YABUN(i,K)_od; _Do i=1, Natom; Yat(i)=0.d0; If(i==14)Yat(i)=1.d0/AZ(14); _od; write(@term,'(a,1p3e10.2)') ' Test happal for Fe, pl*Urho, Tp*UTp', Yat(14)*56.,pl*Urho,Tp*UTp; RADP=.FALSE.; CALL URSOS; kmhap=K; call Happa; write(@term,'(a)') ' Test table happal:'; write(@term,'(1p10e10.2)') happal; call HappaDirect; write(@term,'(a)') ' Test direct happal:'; write(@term,'(1p10e10.2)') happal; _od %R_test: _Repeat READ(8,'(A)',End=110,Iostat=kio) STRING; WRITE(@Wres,'(A)') string; _until kio^=0; 110: stop; %RD: Nc=8; @ReadData EPS,HMIN,HMAX;-- STIFF ACC.,MIN.STEP(SEC),MAX.STEP(SEC) @wterm ' 1'; @ReadData METHN,JSTART,MAXORD,KNadap; @wterm ' 1'; -- ADAMS(1)-GEAR(2)-BGH(3), 0-START , max order, -- KNadap<0 for Ni from a foreign model -- abs(KNadap): 1 freqmn - geom. progression direcr opacity -- 2 freqmn - Peter Hoeflich, opazr -- 3 freqmn - Peter Hoeflich, opahoef, hap2int -- 4 freqmn - Ronfict, opahoef, hapintron -- 5 freqmn - as 4, but homolog. vel. -- 6 freqmn - as 4, but istim==0 always => hapron0 @ReadData NSTA,NSTB,TcurA,TcurB;-- Steps & time in days @wterm ' 1'; -- For read: -- if NSTA<0 then files 10 & 11, else files 12 & 13 -- used by STRAD for test runs, beginning -- with first saved model having NSTEP>=abs(NSTA) -- For Write: -- if NSTB<0 then files 10 & 11, else files 12 & 13 -- To begin test runs with zero step (file 9) use -- NSTA<0, Begrun==T, & NSTB==1 -- To begin test runs with some step in file 10 use -- NSTA<0, Begrun==T, NSTB>0 & NSTB^=1 & NSTB>abs(NSTA) -- To begin test runs with last step in file 10 use -- NSTA<0, Begrun==T, NSTB>0 & NSTB^=1 & abs(NSTB)<abs(NSTA) -- To continue test runs in files 12 & 13 use -- NSTA>0, Begrun==F, NSTB>0 & NSTB<abs(NSTA), e.g. NSTB=1 -- For STINFO: ???? -- use always Begrun=T & NSTB^=1 (if NSTB==1 then initial model) -- ???? @ReadData NSTMAX,NDebug,NOUT,IOUT,MBATCH; -- Max STEP number, Debug step number, INTERVAL OUTPUT, LINE PRINT @wterm ' 1'; Mbatch=Min(Mbatch,Lcurdm); -- MASS OF NI CORE, NI 56 ABUNDANCE @ReadData AMNI,XMNI; -- contaminated and true nickel initial mass @wterm ' 1'; @ReadData AMHT,EBurst,tBurst,tbeght;-- Mass of Heated Core, Energy & tm @wterm ' 1'; @ReadData EKOB,AI1,AI2,AI3,US; -- KINETIC ENERGY(E+50) @wterm ' 1'; -- Mass fractions TRIANGLE VEL.PROFILE, OUT/INWARD +1/-1 -- THRESH.JAC., CK RAPID, CONVECTION, EDD TM DEPEND, CHANGE NCND @ReadData THRMAT,CRAP,CONV,EDTM,CHNCND,Givdtl; @wterm ' 1'; -- FLOOR FOR R V TP YCARBINV FJ - FH @ReadData FLOOR(1),FLOOR(2),FLOOR(3),FLOOR(4); @wterm ' 1'; -- Wacc(R) Wacc(V) Wacc(T) Wacc(RADIAT): @ReadData Wacc(1), Wacc(2), Wacc(3), Wacc(4); @wterm ' 1'; -- FitTau TauTol Rvis AQ BQ DRT NRT SCAT @ReadData FitTau,TauTol,Rvis,AQ,BQ,DRT,NRT,SCAT; @wterm ' 1'; @ReadData NnTO; -- number of outputs _do ito=1,NnTO; Read(Nc,*) TO(ito); _od; If (LSystem) then; -- for IBM the read BEGRUN is ignored and redefined here READ(5,*) IRC; -- for IBM channel 5 is from stack If (IRC==0) then; BEGRUN=.FALSE.; else; BEGRUN=.TRUE.; endif; endif; %RM: -- READ MODEL IF(BEGRUN & IABS(NSTB)==1)THEN; <*Start: read channel - 9 -- begin run *>; ELSE; <*Conr: read Sumprf - channel 10/12 *>; <*Curv: read Lcurve - channel 11/13 *>; ENDIF; Z1=0.d0; Z2=0.d0; do K=1,Nzon; Ry(K)=Y(K,1); -- write(@term,'(a,i3,2(a,1p,e12.4))')' k=',k,' dm=',dm(k), -- ' Ry=',Ry(i); Uy(K)=Y(Nzon+K,1); -- standard, change for relativism! If(abs(KNadap)==5)then; If(K<NZON)Then; DM2=DM(K+1); else; DM2=DMOUT; endif; Z1=Z1+Uy(K)**2*(DM(K)+DM2); -- kin.energy Z2=Z2+Ry(K)**2*(DM(K)+DM2); endif; Ty(K)=Y(2*Nzon+K,1); enddo; If(abs(KNadap)==5)then; URout=sqrt(Z1/Z2); do K=1,Nzon; Uy(K)=URout*Ry(K); -- const. homolog. vel. Y(Nzon+K,1)=Uy(K); enddo; endif; <*nickel: adjust AMNI, XNI for KNadap >0 or read file *.xni for KNadap<0 *>; IF(CONV)THEN; do I=1,NZON; Y(I+(NVARS-2)*NZON,1)=UC(I); Y(I+(NVARS-1)*NZON,1)=YAINV(I); enddo; ENDIF; NSTMAX=NSTEP+NSTMAX; IF(MOD(NSTMAX,MBATCH)^=0) THEN; NSTMAX=(NSTMAX/MBATCH+1)*MBATCH; -- CALL CMSCOM(IERR,'EXEC ST$MSG WARNING'); WRITE(@Wres,'( '' NSTMAX CHANGED:'',I6)') NSTMAX; ENDIF; -- IF(JSTART > MAXORD)THEN; JSTART=0; --always -- ELSE; -- seb changed 28 Mar 2001 -- JSTART=-IABS(JSTART); -- IN ORDER TO CHANGE METH, EPS, ETC. -- Jstart < 0 needed to initialize EPSJ in STIFFBGH -- ENDIF; %RM_Start: Lunit=9; NFILE=Model; @wterm '>>',Nfile,'<<'; call StradIO('rm',Lunit,NFILE); -- Define initial values for rad: @wterm ' Eko, Ekob: ',Eko,Ekob; EKO=EKO+EKOB; -- EKO comes from Eve in @IOSTART, but usually = 0 If(EKO>0.)then; <*Vel: define Uy() for nonzero EKO *>; endif; NCND=NZON; -- STRAD defines real NCND and NFRUS -- any other value for NCND is not allowed here NFRUS=NFREQ; KRAD=(NZON-NCND)*NFRUS; TAUOLD=0.D0; %RM_Conr: if(NSTA<0)then; Lunit=10; NFILE=Sumprf; else; Lunit=12; NFILE='test.prf'; endif; call StradIO('cm',Lunit,NFILE); WRITE(@term,*)' Begrad READ STEP=',NSTEP; %RM_Curv: if(NSTA<0)then; Lunit=11; NFILE=Sumcur; else; Lunit=13; NFILE='test.crv'; endif; call StradIO('rc',Lunit,NFILE); %RM_nickel: km=1; print *, 'AMINI=',Amini,UM; _While (AM(km)-AMini)*UM <= AMNi*1.0000000001d0 & km<nzon _Do print *, km,AM(km), AMNi*1.0000000001d0; km=km+1; _Od; -- (AM(km)-AMini)*UM > AMNi If ( km >= nzon ) then; write(@term,*) ' Begrad: AMNI in error! km=',km,'>=nzon =',nzon; stop; Endif; kmnick=km-1; AMNi=(AM(kmnick)-AMini)*UM; XNI=XMNI/AMNi; --used in VOLEN If ( XNI >= 1.d0 & KNadap >0 ) then; write(@term,*) ' Begrad: AMNI too small ! XNI=',XNI; stop; Endif; If(KNadap>0)then; <*Fe: adjust Ni to Fe for eq.state *>; else; <*Nifor: Ni from a foreign model *>; endif; %RM_nickel_Fe: iferrum=1; _While iferrum <= Natom & Zn(Iferrum) <> 26 _Do iferrum=iferrum+1; _Od; If (iferrum > Natom ) Then; write(@term,*) ' Begrad: Ferrum not found !'; stop; Endif; _Do km=1,kmnick; YABUN(iferrum,km)=XNI/AZ(iferrum); sum=0.d0; _Do j=1,Natom; if(j<>iferrum) sum=sum+YABUN(j,km)*AZ(j); _Od; If(sum<1.d-5)Then; write(@term,*) ' Begrad: Sum too small=',sum; stop; Endif; _Do j=1,Natom; if(j<>iferrum) YABUN(j,km)=YABUN(j,km)*(1.d0-XNI)/sum; _Od; _Od; %RM_nickel_Nifor: @Wterm 'Read .XNI'; read(28) XNIfor; /* read(28) IOlen1, -- for WATCOM (XNifor(i),i=1,Nzon) ,IOlen2; -- for WATCOM */ Xmni=0.d0; _Do km=1,Nzon; Xmni=Xmni+Xnifor(km)*dm(km)*UM; _od; write(@term,*)' Ni mass=',Xmni; %RU: UR=10.D0**ULGR; CLIGHT=CS*UTIME/UR; IF(CRAOLD^=CRAP)THEN; T=T*(CRAOLD/CRAP);H=H*(CRAOLD/CRAP); CK1=CK1*(CRAOLD/CRAP);CK2=CK2*(CRAOLD/CRAP); ENDIF; UFREQ=BOLTZK*UTP/(2.d0*PI*HPLANC); -- WRITE(@Wres,*)' UFRC UFREQ:',UFRC,UFREQ; -- UFRP=BOLTZK/(2.D0*PI*hPlanc); -- WRITE(@Wres,*)' UFRP UFREQ:',UFRP,UFRP*UTP; CKRAD=6.D+01/PI**4*CSIGM*UTP**4*UTIME**3/(URHO*UR**3); CCL=CS*1.D+08/UFREQ; -- TO TRANSFORM FROM ANGSTREMS & V.V. CFLUX=60.D0*CSIGM*(UTP/PI)**4; CLUM=32.D0*PI/3.D0*(CSIGM*UR*UTP**4/URHO); -- FL0 INTO LUMINOSITY CLUMF=4.D0*PI*UR**2*CFLUX; -- FH INTO LUMINOSITY CIMP =CFLUX*UTIME**2/(CS*UR**2*URHO); -- Cimp=0.d0; -- @wterm '!!!!!****>>>> Cimp=0 for TEST with ZERO RaD.FORCE !'; %R_haptab: -- KNadap<0 for Ni from a foreign model -- abs(KNadap): 1 freqmn - geom. progression, direct opacity -- 2 freqmn - Peter Hoeflich, opazr -- 3 freqmn - Peter Hoeflich, opahoef, hap2int -- 4,5 freqmn - Ronfict, opahoef, hapintron -- 6 freqmn - Ronfict, opahoef, hapron0 istim==0) If (ABS(Knadap) <=3 .and. ABS(Knadap).ne. 1) then; -- see begrad for Knadap Open(unit=2,file=Opafile,form='unformatted'); Else; _do nf=1,@ns; nu=30+nf; app=char(ichar('0')+nf); -- nf must be <=9 Opafile=Opafile(1:Length(Opafile)-1)//app; Open(unit=nu,file=Opafile,form='unformatted'); @wterm ' opened unit=',nu,' file=',Opafile; _od; Open(unit=29,file=Opafile(1:Length(Opafile)-1)//'ab',form='unformatted'); endif; _case indop(abs(KNadap)) _of _1 _2 If(^petread)then; read(2) Nfreq0,Msta,Nrho,Ntp,YATab,(Wavel(iif),iif=1,@Nfreq), TpTab,RhoTab,STab,EpsBand,EppBand; petread=.true.; endif _3 _do ihp=1,@ns; read(30+ihp) nw,Stime,Nfreq0,Msta,Nrho,NTp, dumWavel,dumFreq,dumFreqmn, -- read into dummy files --OLD TpTab,RhoTab,hpbanab2,hpbansc2; TpTab,RhoTab,hpbansc2; -- actual _do im=1,Nzon/@skip; _do iro=1,Nrho; _do itp=1,NTp; _do L=1,Nfreq0; hpsavsc(L,itp,iro,im,ihp)=hpbansc2(L,itp,iro,im); _od_od_od_od /* write(@wres,*)' nw,Stime,Nfreq0,Msta,Nrho,NTp', nw,Stime,Nfreq0,Msta,Nrho,NTp; write(@wres,*)' dumWavel'; write(@wres,*)dumWavel; write(@wres,*)' dumFreq'; write(@wres,*)dumFreq; write(@wres,*)' dumFreqmn',dumFreqmn; stop ' test read ' */; @wterm ' read unit=',30+ihp; close(30+ihp); If(ihp==nw ! ABS(knadap)==6)then; stmlog(ihp)=log(Stime); -- in days; m.b. better to save in Ronfict? @wterm ' ihp stmlog=',stmlog(ihp),ihp; else; @wterm ' warning error reading haptab, ihp=',ihp,' nw=',nw; -- here only faked tables are used -- stop 324; endif; _od read(29) hpbanab2; -- read absorption separately without expansion effect close(29); _esac; @wterm ' Nfreq: ',Nfreq0,Nfreq,' Msta,Nrho,Ntp: ',Msta,Nrho,Ntp; -- only two sets of tables are in the ROM for -- two values of dvdr (given by Stime) -- istim is the flag for the regime of saving the tables at given t -- for istim==0 , t<t0=Stime(1)/(half of a factor of Stime step), -- and istim==@ns+1, t>Stime(@ns), use only Stime(@ns) -- the -- smallest dvdr without interpolation; -- -- for istim==1 , t0<t<Stime(1), log interp -- between Stime(@ns) and Stime(1) using t in (0.5, 1)*Stime(1) -- -- for istim==2, 3, ... @ns, Stime(istim-1)<t<Stime(istim), -- log interp between boundaries -- actually Stime is a scalar and the vector is -- stmlog(.)==log(Stime for a respective moment of time) saved in Begrad -- tdlog = log(max(t*Utime/8.6d+04,hmin)); -- for Nstep=0, t=0 defined in EVE and -- saved in *.mod _select _ tdlog <= stmlog(1)-(stmlog(2)-stmlog(1))/2.d0 .or. tdlog > stmlog(@ns) .or. abs(Knadap)==6 [istim=0] _ tdlog > stmlog(1)-(stmlog(2)-stmlog(1))/2.d0 & tdlog <= stmlog(1) [istim=1] _other -- here: stmlog (1) < tdlog <= stmlog(@ns) [ istim=2; _while stmlog(istim) < tdlog _do istim=istim+1 _od ] -- stmlog(istim) >= tdlog _end @wterm ' istim=',istim; istold=istim; -- for Strad _case istim+1 _of _0 _do im=1,Nzon/@skip; _do iro=1,@Ntab; _do itp=1,@Ntab; _do L=1,Nfreq; hpbanab1(L,itp,iro,im)=hpbanab2(L,itp,iro,im); hpbansc1(L,itp,iro,im)=hpbansc2(L,itp,iro,im); _od_od_od_od thaplog1=stmlog(1); thaplog2=stmlog(@ns); -- here the values of thaplog are not important -- (if they are not equal), since haptab is the same _1 _do im=1,Nzon/@skip; _do iro=1,@Ntab; _do itp=1,@Ntab; _do L=1,Nfreq; hpbanab1(L,itp,iro,im)=hpbanab2(L,itp,iro,im); hpbansc1(L,itp,iro,im)=hpbansc2(L,itp,iro,im); _od_od_od_od thaplog1=stmlog(1)-(stmlog(2)-stmlog(1))/2.d0; thaplog2=stmlog(1) _2 thaplog1=stmlog(1); thaplog2=stmlog(2) _3 thaplog1=stmlog(2); thaplog2=stmlog(3) _4 thaplog1=stmlog(3); thaplog2=stmlog(4) _5 thaplog1=stmlog(4); thaplog2=stmlog(5) _esac If(istim^=0)then; -- for istim==0 already read app=char(ichar('0')+istim); -- istim must be <=9 Opafile=Opafile(1:Length(Opafile)-1)//app; Open(unit=30,file=Opafile,form='unformatted'); read(30) nw,Stime,Nfreq0,Msta,Nrho,NTp, dumWavel,dumFreq,dumFreqmn, -- read into dummy files -- TpTab,RhoTab,hpbanab2,hpbansc2; TpTab,RhoTab,hpbansc2; close(30); endif; If(istim^=0 & istim^=1)then; -- for istim==0 & 1 already read app=char(ichar('0')+istim-1); -- istim must be <=9 Opafile=Opafile(1:Length(Opafile)-1)//app; Open(unit=30,file=Opafile,form='unformatted'); read(30) nw,Stime,Nfreq0,Msta,Nrho,NTp, dumWavel,dumFreq,dumFreqmn, -- read into dummy files -- TpTab,RhoTab,hpbanab1,hpbansc1; TpTab,RhoTab,hpbansc1; close(30); close(29); Open(unit=29,file=Opafile(1:Length(Opafile)-1)//'ab',form='unformatted'); read(29) hpbanab1; close(29); endif; <*check: *>; %R_haptab_check: @wterm 'check nw, Stime, Nfreq0, Msta, Nrho, NTp:'; @wterm nw,Stime,Nfreq0,Msta,Nrho,NTp; _do im=1,Nzon/@skip; _do iro=1,@Ntab; _do itp=1,@Ntab; _do L=1,Nfreq; hbab=hpbanab2(L,itp,iro,im); hbal=hpbansc2(L,itp,iro,im); if (hbal >1.d50 .or. hbab > 1.d50)then; write(@term,'(4(a,i4),1p,2(a,e12.4))') ' L=',L,' itp=',itp,' iro=',iro,' im=',im, ' hbal=',hbal,' hbab=',hbab; -- stop ' Begrad: bad happa !!!'; endif; _od_od_od_od %RF: _case indfr(abs(KNadap)) _of _1 -- our old definition: FREQ(1)=CS*1.D+08/(WLMAX*UFREQ); FREQ(NFREQ+1)=CS*1.D+08/(WLMIN*UFREQ); dumFreq(1) = Cs * 1.d+08 / wlmax; dumFreq(nfreq+1) = Cs * 1.d+08 / wlmin; Basis=(WLMAX/WLMIN)**(1.D0/(DBLE(NFREQ))); -- Geometric PROGRESSION _Do i=2,NFREQ; FREQ(i)=FREQ(I-1)*BASIS; dumFREQ(i)=FREQ(i)*Ufreq; _OD; _DO L=1,NFREQ; FREQMN(L)=SQRT(FREQ(L)*FREQ(L+1)); dumFREQMN(L)=FREQMN(L)*Ufreq; dumWavel(L) = Cs * 1.d+08 / dumFreqmn(L); _OD; -- _2 -- from P.Hoeflich: If(^petread)then; read(2) Nfreq0,Msta,Nrho,Ntp,YATab,(Wavel(iif),iif=1,@Nfreq), TpTab,RhoTab,STab,EpsBand,EppBand; petread=.true.; endif; If(Nfreq>@Nfreq)then; _Do L=@Nfreq+1,Nfreq; Wavel(L)=Wavel(L-1)/2.d0; _od; endif; _Do L=1,NFREQ; FREQMN(L)=CS*1.D+08/(WaveL(L)*UFREQ); _od; FREQ(1)=0.5d0*Freqmn(1); FREQ(Nfreq+1)=2.d0*Freqmn(Nfreq); _Do L=2,NFREQ; FREQ(L)=0.5d0*(Freqmn(L-1)+Freqmn(L)) _od; _3 -- from Ronfict (now the same as _1 but CAUTION here!) _Do i=1,NFREQ+1; FREQ(i)=dumFREQ(i)/Ufreq; _od; _Do L=1,NFREQ; FREQMN(L)=dumFREQMN(L)/Ufreq; _od; _esac write(@term,'(2(a,i5))') ' KNadap=',KNadap, ' indfr(abs(KNadap))=',indfr(abs(KNadap)); write(@term,'(a,1pe12.4)') ' Lam max,AA=',CS*1.D+08/(FREQ(1)*UFREQ); write(@term,'(a,1pe12.4)') ' Lam min,AA=',CS*1.D+08/(FREQ(NFREQ+1)*UFREQ); -- pause; -- for all cases: _DO L=1,NFREQ; WEIGHT(L)=(FREQ(L+1)-FREQ(L))*FREQMN(L)**3; IF(L<NFREQ)THEN; DLOGNU(L)=1.D0/LOG(FREQMN(L)/FREQMN(L+1)); ELSE; DLOGNU(L)=1.D0/LOG(FREQMN(L)/FREQ(L+1)); ENDIF; _OD; %R_freqob: do L=1,Nfreq; freqob(L)=freqmn(L); enddo; If(Mfreq>Nfreq)then; do L=Nfreq+1,Mfreq; freqob(L)=freqmn(Nfreq)*exp(-dble(L-Nfreq)/dlognu(1)); enddo; endif; %R_BANDS: LubvU=NFRUS; _WHILE LubvU>1 & FREQMN(LubvU) > @U _DO LubvU=LubvU-1; _OD; -- FREQMN(LubvU) <= @U OR LubvU==1 LubvB=NFRUS; _WHILE LubvB>1 & FREQMN(LubvB) > @B _DO LubvB=LubvB-1; _OD; -- FREQMN(LubvB) <= @B OR LubvB==1 LubvV=NFRUS; _WHILE LubvV>1 & FREQMN(LubvV) > @V _DO LubvV=LubvV-1; _OD; -- FREQMN(LubvV) <= @V OR LubvV==1 LubvR=NFRUS; _WHILE LubvR>1 & FREQMN(LubvR) > @R _DO LubvR=LubvR-1; _OD; -- FREQMN(LubvR) <= @R OR LubvR==1 LubvI=NFRUS; _WHILE LubvI>1 & FREQMN(LubvI) > @I _DO LubvI=LubvI-1; _OD; -- FREQMN(LubvI) <= @I OR LubvI==1 Lyman=NFRUS; _WHILE Lyman>1 & FREQMN(Lyman) > @Lyman _DO Lyman=Lyman-1; _OD; -- FREQMN(LubvV) <= @V OR LubvV==1 WRITE(@Wres,*)' L UBVRI Lyman:',LubvU,LubvB,LubvV,LubvR,LubvI,Lyman; %RM_Start_Vel: <*tri: define initial velocity triangle profile i1, i2, i3 *>; Z1=US/(AM(I2)-AM(I1)); IF(I2^=I3) Z2=US/(AM(I2)-AM(I3)); _Do K=I1,I3; If(K<=I2)then; Uy(K)=Z1*(AM(K)-AM(I1)); else; Uy(K)=Z2*(AM(K)-AM(I3)); endif _od; Z1=0.D0; K1=I1+1; _Do K=K1,I3; if(K<Nzon)then; Z1=Z1+Uy(K)**2*(DM(K)+DM(K+1)); else; Z1=Z1+Uy(K)**2*(DM(K)+DMout); endif; _od; _Do K=I1,I3; Y(Nzon+K,1)=Uy(K)*2.D0*SQRT((EKO/UEPRI)/Z1) _od; %RM_Start_Vel_tri: i=0; _repeat i=i+1; _until (am(i)-amini)/(amout-amini) >=ai1; i1=i; write(@term,'(a,i5,1p,e10.3)') ' i1, ai1:', i1, ai1; _repeat i=i+1; _until (am(i)-amini) /(amout-amini) >=ai2; i2=i; write(@term,'(a,i5,1p,e10.3)') ' i2, ai2:', i2, ai2; _repeat i=i+1; _until (am(i)-amini) /(amout-amini) >=ai3 .or. i==nzon; i3=i; write(@term,'(a,i5,1p,e10.3)') ' i3, ai3:', i3, ai3; %RH: WRITE(@Wres,'(''%RUN:'')'); WRITE(@Wres,'(//30X,''<===== HYDRODYNAMIC RUN OF MODEL '',A, ''=====>'',/)') Sumprf; WRITE(@Wres,'(30X,''MASS(SOLAR)='',F6.3,7X,''RADIUS(SOLAR)='',F9.3/, 30X,''EXPLOSION ENERGY(10**50 ERG)='',1P,E12.5,/)') AMOUT*UM,RADBEG,EKO+Eburst; WRITE(@Wres,'(30X,''<====='',33X,''=====>'',//)'); WRITE(@Wres,'('' INPUT PARAMETERS '')'); WRITE(@Wres,'('' EPS = '',F15.5,9X,'' Rce = '',1P,E15.5,A)') EPS, RCE*UR/RSOL,' SOL.Rad.'; WRITE(@Wres,'('' HMIN = '',1P,E15.5,A,7X,'' AMht = '',1P,E15.5,A)') HMIN*UTIME,' S', AMht,' SOL.MASS'; WRITE(@Wres,'('' HMAX = '',1P,E15.5,A,7X,'' Tburst= '',1P,E15.5,A)') HMAX*UTIME,' S', TBurst,' S'; WRITE(@Wres,'('' THRMAT= '',1P,E15.5,9X,'' Ebstht= '',1P,E15.5,A)') THRMAT,EBurst,' 1e50 ergs'; WRITE(@Wres,'('' METH = '',I15,9X,'' CONV = '',L15)')METH,CONV; WRITE(@Wres,'('' JSTART= '',I15,9X,'' EDTM = '',L15)')JSTART,EDTM; WRITE(@Wres,'('' MAXORD= '',I15,9X,'' CHNCND= '',L15)')MAXORD,CHNCND; WRITE(@Wres,'('' NSTA = '',I15,9X,'' FitTau= '',1P,E15.5)')NSTA,FitTau; WRITE(@Wres,'('' NSTB = '',I15,9X,'' TauTol= '',1P,E15.5)')NSTB,TauTol; WRITE(@Wres,'('' NOUT = '',I15,9X,'' IOUT = '',I15)')NOUT,IOUT; WRITE(@Wres,'('' TcurA = '',F15.5,9X,'' Rvis ='',F15.5)') TcurA,Rvis; WRITE(@Wres,'('' TcurB = '',F15.5,9X,'' BQ = '',1P,E15.5)')TcurB,BQ; WRITE(@Wres,'('' NSTMAX= '',I15,9X,'' DRT = '',1P,E15.5)')NSTMAX,DRT; WRITE(@Wres,'('' XMNI = '',1P,E15.5,A,'' NRT = '',I15)') XMNI,' SOL.MASS',NRT; If(KNadap>0)then; WRITE(@Wres,'('' XNI = '',1P,E15.5)')XNI; WRITE(@Wres,'('' CONTM.= '',1P,E15.5,A,'' SCAT = '',L15)') AMNI,' SOL.MASS',SCAT; else; WRITE(@Wres,'('' XNifor= '',1P,E15.5)')XNifor(1); WRITE(@Wres,'('' MNicor= '',1P,E15.5,A,'' SCAT = '',L15)') AMNI,' SOL.MASS',SCAT; endif; PRINT @F1,UTP,UTIME,URHO,UFREQ; PRINT @F2,CK1,CK2,CFR,CKRAD,CFLUX,CLUM,CLUMF; @F2: FORMAT(' CK1 =',1P,E12.5,' CK2=',E12.5,' CFR=',E12.5, ' CKRAD=',E12.5/ ' CFLUX=',E12.5,' CLUM=',E12.5,' CLUMF=',E12.5); @F1: FORMAT(' UTP=',1P,E12.5,' UTIME=',E12.5,' URHO=',E12.5, ' UFREQ=',E12.5); PRINT'('' *****CK1, CK2 INCREASED'',G8.1,'' TIMES*****'')',CRAP; WRITE(@Wres,'('' FLOOR :''/1X,1P,10E10.2)') FLOOR; PRINT *,' N DIFF. EQS=',N,' N RAD. EQS=',2*KRAD; %S:O SUBROUTINE @SAVE_MODEL; IMPLICIT REAL*8 (A-H,O-Z); INTEGER LUNIT; Character*80 NFILE*80; _INCLUDEN snrad; _INCLUDEN abo; if(NSTB<0)then; Lunit=10; NFILE=Sumprf; else; Lunit=12; NFILE='test.prf'; endif; call StradIO('wm',Lunit,NFILE); --close(Lunit); if(NSTB<0)then; Lunit=11; NFILE=Sumcur; else; Lunit=13; NFILE='test.crv'; endif; call StradIO('wc',Lunit,NFILE); --close(Lunit); RETURN;END; %I:O BLOCK DATA STDATA; IMPLICIT REAL*8(A-H,O-Z); _include zone; COMMON/QNRGYE/QNUC,RGASA,YELECT; COMMON/HNUSED/HUSED,NQUSED,NFUN,NJAC,NITER,NFAIL; -- COMMON/AQ/AQ,BQ,DRT,NRT; -- BQ pseudo-viscosity for R-T COMMON/AZNUC/ACARB,ZCARB,ASI,ZSI,ANI,ZNI,QCSI,QSINI; COMMON/TOO/TOO,KO,KNTO,TO(KOMAX),STO(KOMAX),NTO(KOMAX); -- DATA AQ/2.D0/; -- standard -- DATA AQ/8.D0/; -- enhanced <*A: DATA FOR TOO *>; <*P* DATA FOR Peter's prescribed moments *>; DATA QNUC,YELECT/8.8861D+6,0.5D0/; DATA ACARB,ZCARB,ASI,ZSI,ANI,ZNI,QCSI,QSINI/12.D0,6.D0,28.D0, 14.D0,56.D0,28.D0,88.861D0,1.884D0/; DATA NFUN,NJAC,NITER,NFAIL/4*0/; END; %IA:O DATA TO/komax*1.D+19/,STO/komax*1.D-5/,NTO/komax*0/; -- no prescribed moments /* DATA TO/ 8.0000d-02, 8.0116d-02, 8.0231d-02, 8.0289d-02, 8.0347d-02, 8.0463d-02, 8.0694d-02, 8.0926d-02, 8.1157d-02, 8.1620d-02, 8.2083d-02, 8.2315d-02, 8.3000d-02, 8.5000d-02, 8.8000d-02, 9.2000d-02, 9.6000d-02, 1.0000d-01, 1.0500d-01, 1.1000d-01, 1.2000d-01, 1.4000d-01, 1.6000d-01, 2.0000d-01, 2.5000d-01, 3.0000d-01, 4.0000d-01, 5.0000d-01, 7.0000d-01, 1.0000d+00, 1.5d+00, 2.d0, 3.d+00, 4.d0, 5.d+00, 6.d0, 8.d0, 10.d0, 30.d0, 90.d0 -- in days /, STO/komax*1.D-5/, NTO/komax*0/; */ -- STO always in secs %IP:O DATA TO/ 8.0000d-02, 8.0116d-02, 8.0231d-02, 8.0289d-02, 8.0347d-02, 8.0463d-02, 8.0694d-02, 8.0926d-02, 8.1157d-02, 8.1620d-02, 8.2083d-02, 8.2315d-02, 8.3000d-02, 8.5000d-02, 8.8000d-02, 9.2000d-02, 9.6000d-02, 1.0000d-01, 1.0500d-01, 1.1000d-01, 1.2000d-01, 1.4000d-01, 1.6000d-01, 2.0000d-01, 2.5000d-01, 3.0000d-01, 4.0000d-01, 5.0000d-01, 7.0000d-01, 1.0000d+00, 1.5000d+00, 3.0000d+00, 5.0000d+00/, -- zdesj momenty v -- v sutkah ili sekundah STO/komax*0.D0/, NTO/komax*0/;
module mStiffbgh; use zone_environment_types, only: zone_environment_set_type; use cp_error_handling, only: cp_error_type; use mHcdfnrad, only: difmat; implicit real*8 (A-H,O-Z); -- implicit none; private; character(len=*), parameter, private :: mdl_name = 'mStiffbgh'; public vtstif; contains; -- GEAR integrator for ODE with sparse Jacobian and option METH=3 -- for BRAYTON, GUSTAVSON & HACHTEL method _Define @DNAMEDVAR DNDVAR -- work array of derivatives @STRING STR020 -- for debugger @LINE LIN130 -- for debugger @THR 1.D-12 -- Zlatev threshold -- @THR 1.D-15 -- Zlatev threshold --_TRACE "WRITE(@Wres,*)' ifail =',ifail,' NumSing=',NUMSING," --_TRACE "WRITE(@Wres,*)' Y11 =',Y(1,1),' Ytp11 =',Y(301,1)," --_TRACE "IF(NSTEP>8000) write(*,'(a,L3,1p,2e12.3,a)')'EVALJA H Hused=',EVALJA,H,Hused," --_TRACE "@wterm' lmax=',lmax,' maxord=', maxord," --_TRACE "@wterm' EVALJA=',EVALJA,' maxer=', maxer," --_TRACE "write(@term,'(a,4L3,i4,a)')' EVALJA needbr corrco badste ifail=', -- EVALJA, needbr, corrco, badste, ifail," --_TRACE "@wterm' EVALJA=',EVALJA,' ifail =',ifail,' H=',H," --_TRACE "WRITE(*,*)' N=',N,' Nzmod=',Nzmod," subroutine vtstif(envs, err); IMPLICIT REAL*8(A-H,O-Z); /* VTSTIF performs one step of the integration of an initial value problem for a system of ordinary differential equations. Communication with VTSTIF is done with the following variables: Y an NYDIM by LMAX array containing the dependent variables and their scaled derivatives. LMAX is currently 13 for the Adams methods and 6 for the Gear methods (though LMAX=5 is recommended due to stability considerations). LMAX-1 is MAXORD, the maximum order used. See subroutine COSET. Y(I,J+1) contains the J-th derivative of Y(I), scaled by H**J/Factorial(J). Only Y(I), 1 <= I <= N, need be set by the calling program on the first entry. If it is desired to interpolate to non-mesh points, the Y array can be used. If the current step size is H and the value at t+E is needed, form S = E/H, and then compute NQ Y(I)(t+E) = SUM Y(I,J+1)*S**J . J=0 The Y array should not be altered by the calling program. When referencing Y as a 2-dimensional array, use a column length of NYDIM, as this is the value used in STIFF. N The number of first order differential equations. N may be decreased on later calls if the number of active equations reduces, but it must not be increased without calling with JSTART=0. NYDIM A constant integer >= N, used for dimensioning purposes. NYDIM must not be changed without setting JSTART=0. t The independent variable. t is updated on each step taken. H The step size to be attempted on the next step. H may be adjusted up or down by the routine in order to achieve an economical integration. However, if the H provided by the user does not cause a larger error than requested, it will be used. To save computer time, the user is advised to use a fairly small step for the first call. It will be automatically increased later. H can be either positive or negative, but its sign must remain constant throughout the problem. HMIN The minimum absolute value of the step size that will be used for the problem. On starting this must be much smaller than the average Abs(H) expected, since a first order method is used initially. HMAX The maximum absolute value of the step size that will be used for the problem. EPSJ The relative error test constant. Single step error estimates divided by YMAX(I) must be less than this in the Euclidean norm. The step and/or order is adjusted to achieve this. METH The method flag. METH=1 means the implicit ADAMS methods METH=2 means the GEAR method for stiff problems METH=3 means the GEAR method, corrected by R.K.Brayton, F.G.Gustavson & G.D.Hachtel in the Proceedings of the IEEE v.60, No.1, January 1972, p.98-108 YMAX An array of N locations which contains the maximum absolute value of each Y seen so far. ERROR An array of N elements proportional to the estimated one step error in each component. KFLAG A completion code with the following meanings: 0 the step was succesful -1 the requested error could not be achieved with abs(H)=HMIN -2 corrector convergence could not be achieved for abs(H)>HMIN. On a return with KFLAG negative, the values of t and the Y array are as of the beginning of the last step, and H is the last step size attempted. JSTART An integer used on input and output. On input it has the following values and meanings: 0 perform the first step. This value enables the subroutine to initialize itself. >0 take a new step continuing from the last. Assumes the last step was succesful and user has not changed any parameters. <0 repeat the last step with a new value of H and/or EPS and/or METH. This may be either in redoing a step that failed, or in continuing from a succesful step. On exit, JSTART is set to NQ, the current order of the method. This is also the order of the maximum derivative available in the Y array. After a succesful step, JSTART need not be reset for the next call. AJAC A block used for partial derivatives. It keeps the matrix AJAC=1-EL(1)*H*Jacobian & is computed by DIFMAT. HL = H*EL(1) is communicated to DIFMAT to compute AJAC. FSAVE A block of at least 2*NYDIM locations for temporary storage. The parameters which must be input by the user are: N, NYDIM, t, Y, H, HMIN, HMAX, EPS, METH, JSTART, MAXORD. Additional Subroutines required are: DIFMAT computes DY/Dt given t and Y, and if EVALJA computes AJAC, stored in the form appropriate for M28Y12, COSET, RESCAL, M28Y12, M30Y12, M28BYS - sparse matrix solvers. The calling program must contain the common declarations given in node %C: */ type(zone_environment_set_type), pointer :: envs; type(cp_error_type), intent(inout), optional :: err; character(len=*), parameter :: subrtn_name = 'vtstif' , fullPathSubrtn = mdl_name//'.'//subrtn_name; _INCLUDE snrad; <*C: COMMONS & VAR STIFF *>; <*B: DATA FOR STIFF *>; <*D: DATA FOR SPARSE MATRIX SubroutineS M28Y12 *>; <*F: PREPARATIONS DEPENDING ON JSTART *>; HNT(1)=H; CALL COSET(METH,NQ,EL,TQ,MAXORD,IDOUB);-- needed for continuation --runs RC=RC*EL(1)/OLDL0;OLDL0=EL(1); -- RC is the ratio of new to old values -- of the coefficient H*EL(1). -- When RC differs from 1 by more than RCTEST, -- EVALJA is set to .TRUE. to force -- jacobian to be calculated in DIFMAT. <*A: the constants E,EUP,EDN,BND *>; IF(MOD(NSTEP,MBATCH)==0) NEEDBR=.TRUE.; RETCOD=0; <*G: GENERAL STEP OF STIFF WHILE RETCOD==0 *>; <*X: ALL RETURNS ARE MADE THROUGH THIS SECTION DEPENDING ON RETCOD. H IS USED IN HOLD TO ALLOW THE CALLER TO CHANGE H ON THE NEXT STEP. *>; RETURN; end subroutine vtstif; end module mStiffbgh; %C:O _include nstep; -- COMMON/NSTEP/NSTEP,NDebug,MAXER,IOUT,NOUT; -- COMMON/CREAD/TAUOLD,NSTMAX,MBATCH,MAXORD; _include stsave; -- COMMON/STSAVE/RMAX,TREND,OLDL0,RC,HOLD,EDN,E,EUP,BND,EPSOLD,TOLD, -- MEO,NOLD,NQ,LNQ,IDOUB; REAL*8 ERROR(NYDIM),EL(13),TQ(4); REAL*8 XMOD,ERREST,ERROLD; PARAMETER (NCHANL=6, ETA=1.D-04); -- PIVOT THRESHOLD IN M28BYS PARAMETER (DeltaDeb=1.d-10); -- FOR DEBUGGER CHARACTER @STRING*20,@LINE*130; -- FOR DEBUGGER CHARACTER*20 CharVarName(5); -- FOR DEBUGGER Logical FOUNDI; -- FOR DEBUGGER REAL*8 @DNAMEDVAR(NYDIM),DSAVE(NYDIM);-- FOR DEBUGGER Integer RETCOD; -- RETURN CODE Real Tm1,Tm2; -- Timer for CMS -- CRAY REAL LOGICAL CORRCO; -- CORRECTOR ITERATIONS HAVE CONVERGED LOGICAL NUMSING; -- Numerical singularity in decomposition <*L: COMMON & VARIABLES FOR LINEAR EQS SOLVER F01BRF or M28Y12 *>; %CL:O Parameter(MAXIT=15); LOGICAL LBLOCK,GROW,ABORT(4); COMMON/JSAVE/ASAVE(NZ),IRS(NZ),ICS(NZ); Dimension DB(NYDIM), XSAVE(NYDIM); %B:O DATA ANOISE/1.D-12/;--should be set to the noise level of the machine DATA DFLTZR/1.D-75/, MAXITE/3/, MAXFAI/3/; DATA RMXINI/1.D+4/, RMXNOR /10.D0/, RMXFAI/2.0D0/, IDELAY /10/; DATA RHCORR/.25D0/, RHERR3/.1D0/, RCTEST/.3D0/; DATA BIAS1/1.3D0/, BIAS2/1.2D0/, BIAS3/1.4D0/; DATA CharVarName/'radius','velocity','temperature','FJ','FH'/; %D:O DATA LBLOCK /.True./,GROW/.TRUE./, ABORT/.FALSE.,.TRUE.,.FALSE.,.TRUE./, PIVOT/1.D-01/; %F: KFLAG=0;tOLD=t; NUMSING=.FALSE.; -- WRITE(@term,*)'Entering Stiff NSTEP=', NSTEP; /* On the first call, the order is set to 1 and the initial derivatives are calculated. RMAX is the maximum ratio by which H can be increased in a single step. It is initially RMXINI to compensate for the small initial H, but then is normally equal to RMXNOR. If a failure occurs (in corrector convergence or error test), RMAX is set at RMXFAI for the next increase. */ _SELECT _ JSTART==0 [<*B: INITIAL *> ] _ JSTART<0 [ <*C: reset Y & other variables for changed METH, EPS, N, H *> ] _END; %FB: HUSED=0.D0; NQUSED=0; NOLD=N; EVALJA=.FALSE.; NQ=1; _DO I=1,N; FSAVE(I)=Y(I,1)_OD; HL=H*EL(1); CALL DIFMAT(envs, err); _DO I=1,N; Y(I,2)=FSAVE(NYDIM+I)*H _OD; LNQ=2; IDOUB=LNQ+1; RMAX=RMXINI; EPSJ=EPS*SQRT(Real(N)); TREND=1.D0; OLDL0=1.D0; RC=0.D0; -- RC is the ratio of new to old values -- of the coefficient H*EL(1). -- When RC differs from 1 by more than RCTEST, -- EVALJA is set to .TRUE. to force -- Jacobian to be calculated in DIFMAT. HOLD=H; MEO=METH; EVALJA=.TRUE.; OLDJAC=.FALSE.; LMAX=MAXORD+1; EPSOLD=EPSJ; %A: EDN=(TQ(1)*EPSJ)**2;-- is to test for decreasing the order, E=(TQ(2)*EPSJ)**2; -- comparison for errors of current order NQ, EUP=(TQ(3)*EPSJ)**2;-- is to test for increasing the order, BND=(TQ(4)*EPSJ)**2;-- is used to test for convergence of the -- correction iterates %FC: -- If the caller change EPSJ or METH , the constants E,EUP,EDN,BND -- must be reseted. EPSJ=EPS*SQRT(Real(N)); IF(METH==MEO) THEN; IF(EPSJ^=EPSOLD)THEN; EPSOLD=EPSJ; ENDIF; ELSE; IDOUB=LNQ+1; ENDIF; <*C: TEST N==NOLD, H==HOLD ? *>;EPSOLD=EPSJ; %FCC: IF(N^=NOLD)THEN; IDOUB=LNQ+1;EVALJA=.TRUE.;OLDJAC=.FALSE.;NOLD=N; ENDIF; IF(H^=HOLD) THEN; RH=H/HOLD;H=HOLD; CALL RESCAL(Y,NYDIM,RH,RMAX,RC,LNQ); IDOUB=LNQ+1; ENDIF; %G: _WHILE RETCOD==0 _DO IF(ABS(RC-1.D0)>RCTEST .or. MOD(NSTEP,MBATCH)==0) EVALJA=.TRUE.; If(ifail==13)Then; Needbr=.true.; Evalja=.true.; Endif; t=t+h; <*E: Predictor. This section computes the predicted values by effectively multiplying the Y array by the Pascal triangle matrix *> _Repeat -- loop until .not.EVALJA ITER=0; _DO I=1,N; ERROR(I)=0.D0 _OD; _DO I=1,N; FSAVE(I)=Y(I,1) _OD; CORRCO=.FALSE.; -- Corrector convergence. It is tested -- by requiring changes relative to YMAX(I) -- to be less, in Euclidian norm, than BND, -- which is dependent on EPSJ. _WHILE ^CORRCO & ITER<MAXITE _DO <*Cor: Up to MAXITE corrector iterations are taken. The sum of the corrections is accumulated in the vector ERROR(i). It is approximately equal to the Lnq-th derivative of Y multiplied by H**Lnq/(factorial(Lnq-1)*EL(Lnq)), and is thus proportional to the actual errors to the lowest power of H present (H**Lnq). The Y array is not altered in the correction loop. The updated Y vector is stored temporarily in FSAVE. The norm of the iterate difference is stored in D. *> _OD; -- CORRCO ! ITER==MAXITE _IF CORRCO _THEN <*X: The corrector has converged. OLDJAC is set to .true. if partial derivatives were used, to signal that they may need updating on subsequent steps. The error test is made and control passes to node %GXA if it succeeds. *>; _ELSE -- ITER==MAXITE -- The corrector iterations failed to converge in MAXITE -- tries. NITER=NITER+1; IF(OLDJAC)THEN; -- If partials are not up to date , they are -- reevaluated for the next try. EVALJA=.TRUE.; -- Refresh JACOBIAN ELSE; <*Y: Otherwise, the Y array is retracted to its values before prediction, and H is reduced, if possible. If not, a no convergence exit is taken. *>; ENDIF; _FI; _Until ^EVALJA .or. RETCOD^=0; _OD; -- while RETCOD==0 %GE: do J1=1,NQ; do J2=J1,NQ; J=NQ-J2+J1; do I=1,N; Y(I,J)=Y(I,J)+Y(I,J+1); enddo; enddo; enddo; --_OD_OD_OD; %G_Cor: BADSTE=.FALSE.; -- Logical variable BADSTE may be changed by DIFMAT HL=H*EL(1); -- if(Nstep>1750 & EVALJA)then; -- if(Nstep>1586 & EVALJA)then; if(Nstep>=NDebug & EVALJA)then; <*DEB: full debugger of Jacobian *>; end if; if(ifail==13)then; @wterm ' G_Cor: ifail 13'; endif; CALL DIFMAT(envs, err); -- if necessary (i.e. EVALJA=.true.) the partials -- are reevaluated by DIFMAT <*MatrOut : output Jacobian matrix *>; -- stop ' Jacobian is output'; _IF BADSTE _THEN -- Something is wrong in DIFMAT (i.e. Y(I) is out of physical -- range for some I) WRITE(*,*)' G_Cor: BADSTE=',BADSTE; ITER=MAXITE; EVALJA=.FALSE.; OLDJAC=.FALSE.; Ifail=13; -- Variables are given the values such as to guarantee that -- control passes to the node %GY - to try a new smaller step _ELSE _DO Isw=1,N; XSAVE(Isw)=0.D0; _OD; IF(EVALJA) THEN; <*L: L-U Decomposition of matrix AJAC=1-El*H*Jacobian *>; EVALJA=.FALSE.;OLDJAC=.FALSE.;RC=1.D0; NJAC=NJAC+1; ENDIF; <*U: Find in FSAVE(1:N) the solution X of equation AJAC!X>=H*F(Y)-Y(I,2)-ERROR(I), where AJAC is the matrix AJAC=1-El*H*Jacobian, then obtain in ERROR the full correction, in D - the norm of the iterate difference & in FSAVE a new approximation to Y *> If(Ifail^=13)Then; IF(ITER^=0) TREND=MAX(.9D0*TREND,D/D1); -- If ITER > 0 an estimate of the convergence rate constant -- is stored in TREND, and this is used in the convergence test CORRCO=D*MIN(1.D0,2.D0*TREND)<=BND; D1=D; ITER=ITER+1; Else; --WRITE(@Wres,*)'Ifail=',Ifail; ITER=MAXITE; EVALJA=.FALSE.; OLDJAC=.FALSE.; -- Variables are given the values such as to guarantee that -- control passes to the node %GY - to try a new smaller step Endif; _FI; -- FOR BADSTE %G_Cor.L: /*_Do I=1,NZMOD; If(IRS(I)==201 & Nstep>=160 & Nstep<=180) then; WRITE(NCHANL,'(2(A,I5),1P,A,E15.8,A,L3)') ' ROW=',IRS(I),' COL=',ICS(I),' AJAC=',AJAC(I), ' Needbr=',Needbr; endif _od; */ write(Nchanl,*)' Nzmod=',Nzmod; _Do Isw=1,NZMOD; ASAVE(Isw)=AJAC(Isw); _od; --_repeat If(Needbr)then; <*SOLY12: *> IRW=Idisp(11); IF (IFAIL==2)then; NUMSING=.true.; NEEDBR=.true.; ITER=MAXITE; EVALJA=.FALSE.; OLDJAC=.FALSE.; Ifail=13; -- Variables are given the values such as to guarantee that -- control passes to the node %GY - to try a new smaller step else; NUMSING=.false.; endif; else; <*SOLBSF: *> IRW=IW(1); IF (IFAIL==6)then; NUMSING=.true.; NEEDBR=.true.; ITER=MAXITE; EVALJA=.FALSE.; OLDJAC=.FALSE.; Ifail=13; -- Variables are given the values such as to guarantee that -- control passes to the node %GY - to try a new smaller step else; NUMSING=.false.; NEEDBR=.false.; endif; endif; IF (NUMSING) THEN; <*NUMSING: process the numerical singularity in row IR *>; -- stop 26; ENDIF; --_until ^NUMSING ; %G_Cor.L_SOLY12: -- IFAIL=110; -- hard failure IFAIL=111; -- soft failure /* if(NUMSING)then; NZMOD=N; do isw=1,N; IRN(isw)=isw; ICN(isw)=isw; AJAC(isw)=1.d0; enddo; else; */ _DO Isw=1,NZMOD; IRN(Isw)=IRS(Isw); ICN(Isw)=ICS(Isw); _OD; /* endif; */ -- Call VTIME(Tm1); -- CALL F01BRF(N , NZMOD, AJAC, LICN, IRN, LIRN, ICN, PIVOT, -- IKEEP, IW, WJAC, LBLOCK, GROW, ABORT, IDISP, IFAIL); -- write(*,*) ' In stiffbgh.trf NZMOD=',NZMOD; -- pause; CALL M28Y12(N , NZMOD, AJAC, LICN, IRN, LIRN, ICN, PIVOT, IKEEP, IW, WJAC, LBLOCK, GROW, ABORT, IDISP, IFAIL, @THR); NEEDBR=.FALSE.; -- F01BRF or M28Y12 NOT NEEDED -- Call VTIME(Tm2); IF (GROW) WRITE (NCHANL,89996) WJAC(1); 89996:FORMAT (' ON EXIT FROM M28Y12: W(1) = ',1P,E12.4); WRITE (NCHANL,99996)IDISP(2), IDISP(6), IDISP(7), IDISP(3),IDISP(4); IF (LBLOCK) WRITE (NCHANL,99994) (IDISP(I),I=8,10); 99999:FORMAT (6A4, A3); 99998:FORMAT (4(1X/), 1X , 5A4, A3, 'RESULTS'/1X); 99997:FORMAT (5(2X, G8.0, 2X, I1, 2X, I1)); 99996:FORMAT (' NUMBER OF NON-ZEROS IN DECOMPOSITION', 9X, '=', I8 /' MINIMUM SIZE OF ARRAY IRN', 20X, '=', I8 /' MINIMUM SIZE OF ARRAYS A AND ICN', 13X, '=', I8 /' NUMBER OF COMPRESSES ON IRN (IDISP(3))', 7X, '=', I8 /' NUMBER OF COMPRESSES ON A AND ICN (IDISP(4)) =', I8); 99994:FORMAT (' STRUCTURAL RANK', 16X, '=', I7 /' NUMBER OF DIAGONAL BLOCKS', 6X, '=', I7 /' SIZE OF LARGEST DIAGONAL BLOCK =', I7); -- WRITE(NCHANL, -- '('' M28Y12-1 done at step :'',I6,'' Time(mS):'',I6)') -- CRAY '('' M28Y12 done at step :'',I6,'' Time(S):'',1P,G11.5)') -- NSTEP,Tm2-Tm1; %G_Cor.L_SOLBSF: IFAIL=11; -- Call VTIME(Tm1); CALL M28BYS(N , NZMOD, AJAC, LICN, IRS, ICS, ICN, IKEEP, IW, WJAC, GROW, ETA, RPMIN, ABORT(4), IDISP, IFAIL); -- Call VTIME(Tm2); -- WRITE(NCHANL,'('' M28BYS done at step :'',I6,'' Time(ms):'',I6)') -- NSTEP,Tm2-Tm1; /* IF (GROW) WRITE (NCHANL,89995) WJAC(1); WRITE (NCHANL,89994) RPMIN; 89995:FORMAT (' ON EXIT FROM M28BYS: W(1) = ',1P,E12.4); 89994:FORMAT (' VALUE OF RPMIN = ', G12.4); */ IF(IFAIL^=0 .OR. WJAC(1)>1.D50)THEN; IF(IFAIL==6)THEN; <*NUMSING: process the numerical singularity in row IR *>; ENDIF; -- IFAIL=110; -- hard failure IFAIL=111; -- soft failure _Do Isw=1,NZMOD; IRN(Isw)=IRS(Isw); ICN(Isw)=ICS(Isw); AJAC(Isw)=ASAVE(Isw); _od; -- Call VTIME(Tm1); CALL M28Y12(N , NZMOD, AJAC, LICN, IRN, LIRN, ICN, PIVOT, IKEEP, IW, WJAC, LBLOCK, GROW, ABORT, IDISP, IFAIL, @THR); NEEDBR=.FALSE.; -- F01BRF or M28Y12 NOT NEEDED -- Call VTIME(Tm2); /* IF (GROW) WRITE (NCHANL,89996) WJAC(1); WRITE (NCHANL,99996)IDISP(2), IDISP(6), IDISP(7), IDISP(3),IDISP(4); IF (LBLOCK) WRITE (NCHANL,99994) (IDISP(I),I=8,10); */ -- WRITE(NCHANL, -- '('' M28Y12 done again at step :'',I6,'' Time(ms):'',I6)') -- CRAY '('' M28Y12 done at step :'',I6,'' Time(S):'',1P,G11.5)') -- NSTEP,Tm2-Tm1; ENDIF; %G_Cor.L_NUMSING: WRITE(@term,'(3(A,1P,E12.4))') 'Y(IR)=',Y(IRW,1),' Fsave(IR)=',Fsave(IRW), ' Ydot(IR)=',FSAVE(NYDIM+IRW); _DO I=1,NZMOD; IF(IRS(I)==IRW)THEN; WRITE(@term,'(2(A,I5),A,1P,E12.4)') 'IR=',IRS(I), ' IC=',ICS(I),' AJAC=',ASAVE(I); ICW=ICS(I); If(ICW<=NZON*NVARS)then; Izon=MOD(ICW-1,NZON)+1; Ivar=(ICW-1)/NZON +1; WRITE(@term,'(2(A,I5))') 'Ivar=',Ivar,' Izon=',Izon; else; ICW=ICW-NZON*NVARS; If(ICW<=KRAD)Then; L=(ICW-1)/(NZON-NCND)+1; Izon=NCND+MOD(ICW-1,(NZON-NCND))+1; WRITE(@term,'(2(A,I5))') 'FJ 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,I5))') 'FH in L=',L,' Izon=',Izon; endif; endif; ENDIF; _OD; WRITE(@term,'(A)') 'For Ydot of: '; If(IRW<=NZON*NVARS)then; Izon=MOD(IRW-1,NZON)+1; Ivar=(IRW-1)/NZON +1; WRITE(@term,'(2(A,I5))') 'Ivar=',Ivar,' Izon=',Izon; else; IRW=IRW-NZON*NVARS; If(IRW<=KRAD)Then; L=(IRW-1)/(NZON-NCND)+1; Izon=NCND+MOD(IRW-1,(NZON-NCND))+1; WRITE(@term,'(2(A,I5))') 'FJ in L=',L,' Izon=',Izon; else; IRW=IRW-KRAD; L=(IRW-1)/(NZON-NCND)+1; Izon=NCND+MOD(IRW-1,(NZON-NCND))+1; WRITE(@term,'(2(A,I5))') 'FH in L=',L,' Izon=',Izon; endif; endif; %G_Cor.L_SOLBSF_NUMSING:=G_Cor.L_NUMSING: %G_Cor.U: <*K: Find solution by M28CYN *>; If (Ifail^=13) Then; <*V: Find ERROR, D and put new Y in FSAVE *>; Endif; %G_Cor.UK: _DO I=1,N; FSAVE(I+NYDIM)=FSAVE(I+NYDIM)*H-Y(I,2)-ERROR(I); -- R-H SIDE IN FSAVE(NYDIM+1:NYDIM+N) DB(I)=FSAVE(I+NYDIM); -- TO SAVE R-H SIDE _OD; -- ITERATIVE REFINEMENT: ITQ = 0; If(Ifail^=13)Then; _Repeat CALL M28CYN(N,AJAC,LICN,ICN,IKEEP,FSAVE(NYDIM+1), WJAC,IW,1,IDISP,RESID,Ifail); -- SOLUTION IN FSAVE(NYDIM+1:NYDIM+N) If(Ifail^=13)Then; ERROLD=ERREST; ERREST=0.D0; XMOD=0.D0; _DO ID=1,N; XSAVE(ID)=XSAVE(ID)+FSAVE(ID+NYDIM); XMOD=XMOD+ABS(XSAVE(ID)); -- NORMA ERREST=ERREST+ABS(FSAVE(ID+NYDIM)); FSAVE(ID)=0.D0; _OD; _DO K=1,NZMOD; FSAVE(IRS(K))=FSAVE(IRS(K)) -- Polish strategy +XSAVE(ICS(K))*ASAVE(K); _OD; _DO IL=1,N; FSAVE(IL+NYDIM)=DB(IL)-FSAVE(IL); _OD; -- WRITE(NCHANL,'(2(A,1P,G12.3))')' ERROLD=',ERROLD,' ERREST=',ERREST; ITQ = ITQ+1; Else; _DO Isw=1,NZMOD; IRN(Isw)=IRS(Isw); ICN(Isw)=ICS(Isw); AJAC(Isw)=ASAVE(Isw); _OD; NEEDBR=.TRUE.; EVALJA=.TRUE.; @wterm ' G_Cor.UK: ifail 13'; Endif; -- _Until ERREST<=ANOISE*XMOD ! ITQ > MAXIT; _Until ERREST<=ANOISE*XMOD .or. Ifail==13 .or. (ITQ>2 & ERREST>ERROLD) .or. ITQ > MAXIT; IF(ERREST > XMOD*.1d0*EPS) NEEDBR=.TRUE.; -- Esaulov correction -- IF(XMOD>EPS**3)WRITE(NCHANL,*)' RELEST=',ERREST/XMOD; Endif; %G_Cor.UV: D=0.D0; _DO I=1,N; ERROR(I)=ERROR(I)+XSAVE(I); D=D+(XSAVE(I)/YMAX(I))**2; FSAVE(I)=Y(I,1)+EL(1)*ERROR(I) _OD; %GX: D=0.D0; ERMAX=0.D0; MAXER=0; do I=1,N; IF((ERROR(I)/YMAX(I))**2>ERMAX)THEN; ERMAX=(ERROR(I)/YMAX(I))**2; MAXER=I; ENDIF; D=D+(ERROR(I)/YMAX(I))**2; enddo; OLDJAC=.TRUE.; if( D<=E) then; -- After a succesful step, update the Y array and YMAX. -- Consider changing H if IDOUB==1. Otherwise decrease -- IDOUB by 1. -- If a change in H is considered, an increase -- or decrease in order by one is considered also. HUSED=H;NQUSED=NQ;KFLAG=0; IF(METH==3)THEN; -- Save time intervals in HNT do J=7,3,-1; HNT(J)=HNT(J-1) + H; -- HNT(LNQ+1) saves old H for possible order increase -- Initially HNT(2:6) is undefined here! enddo; HNT(2)=H; ENDIF; do J=1,LNQ; do I=1,N; Y(I,J)=Y(I,J)+EL(J)*ERROR(I); enddo; enddo; if ( IDOUB==1 ) then; <*A: Factors PR1, PR2 and PR3 are computed, by which H could be multiplied at order NQ-1, order NQ, or order NQ+1, respectively. The largest of these is determined and the new order is chosen accordingly. If the order is to be increased, we compute one additional scaled derivative. *> RETCOD=3; -- EXIT else; IDOUB=IDOUB-1; IF((IDOUB==1) & (NQ^=MAXORD)) THEN; -- Error is saved for use in a possible -- order increase on the next step. -- _DO I=1,N; Y(I,LMAX)=ERROR(I) _OD; Y(1:N,LMAX)=ERROR(1:N); ENDIF; IF(METH==3)THEN; <*TWO: RESET PR2 *>; RH=PR2; CALL RESCAL(Y,NYDIM,RH,RMAX,RC,LNQ); -- NO IDOUB HERE!!! ENDIF; RETCOD=4; -- EXIT endif; else; -- The error test failed. KFLAG keeps track -- of multiple failures. Restore t and the Y array to their -- previous values, and prepare to try the step again. -- Compute the optimum step size for this or one lower order. NFAIL=NFAIL+1;KFLAG=KFLAG-1;t=tOLD; do J1=1,NQ; do J2=J1,NQ; J=NQ-J2+J1; do I=1,N; Y(I,J)=Y(I,J)-Y(I,J+1); enddo;enddo;enddo; RMAX=RMXFAI; if ( ABS(H)<=(HMIN*1.00001D0) ) then; RETCOD=1; else; if( KFLAG<=-MAXFAI) then; <*Z: Control reaches this section if MAXFAI or more failures have occured. It is assumed that the derivatives that have accumulated in the Y array have errors of the wrong order. Hence the first derivative is recomputed, and the order is set to one. Then H is reduced by factor of RHERR3, and the step is retried. *> else; <*2: RESET PR2 *>; if(NQ==1) then; RH=PR2; CALL RESCAL (Y,NYDIM,RH,RMAX,RC,LNQ); IDOUB=LNQ+1; else; <*1: RESET PR1 *>; if(PR1<=PR2) then; RH=PR2; if(meth==3)then; HNT(1)=H*RH; CALL COSET(METH,NQ,EL,TQ,MAXORD,IDOUB); RC=RC*EL(1)/OLDL0;OLDL0=EL(1); <*E: the constants E,EUP,EDN,BND *>; CALL RESCAL (Y,NYDIM,RH,RMAX,RC,LNQ); else; CALL RESCAL (Y,NYDIM,RH,RMAX,RC,LNQ); IDOUB=LNQ+1; endif; else; <*R: CHANGE NQ,RESET CONST. *>; endif; endif; endif; endif; endif; %GXE=A: %GX1: SUM=0.D0; _DO I=1,N;SUM=SUM+(Y(I,LNQ)/YMAX(I))**2 _OD; PR1=1.D0/(((SUM/EDN)**(.5D0/DBLE(NQ)))*BIAS1+DFLTZR); %GX2: PR2=1.D0/(((D/E)**(.5D0/DBLE(LNQ)))*BIAS2+DFLTZR); %GX_TWO=GX2: %GXR: LNQ=NQ;NQ=NQ-1;RH=PR1; HNT(1)=H*RH; CALL COSET(METH,NQ,EL,TQ,MAXORD,IDOUB); RC=RC*EL(1)/OLDL0;OLDL0=EL(1); <*A: EDN,E,EUP,BND *>; CALL RESCAL(Y,NYDIM,RH,RMAX,RC,LNQ); IDOUB=LNQ+1; %GXRA=A: %GXZ: RH=RHERR3;RH=MAX(HMIN/ABS(H),RH);H=H*RH; HL=H*EL(1); _DO I=1,N; FSAVE(I)=Y(I,1) _OD; CALL DIFMAT(envs, err); _DO I=1,N; Y(I,2)=FSAVE(NYDIM+I)*H _OD; -- EVALJA=.TRUE.; EVALJA=.False.; -- to get out the _repeat loop RC=0.d0; -- to force a new Jacobian IDOUB=IDELAY; IF(NQ^=1) THEN; NQ=1;LNQ=2; CALL COSET(METH,NQ,EL,TQ,MAXORD,IDOUB); OLDL0=EL(1); <*A: EDN,E,EUP,BND *>; ENDIF; %GXZA=A: %GXA: PR3=DFLTZR; IF(NQ^=MAXORD) THEN; <*3: RESET PR3 *>; ENDIF; <*2: RESET PR2 *>; PR1=DFLTZR; IF(NQ^=1) THEN;<*1: RESET PR1 *>;ENDIF; IF((PR2>PR1) & (PR2>PR3))THEN; RH=PR2; ELSE; /* IF(PR1>PR3) THEN; NEWQ=NQ-1;RH=PR1; ELSE; NEWQ=LNQ;RH=PR3; -- PR3 IS GREATER _DO I=1,N;Y(I,NEWQ+1)=ERROR(I)*EL(LNQ)/DBLE(LNQ)_OD; ENDIF; */ -- bug (for METH==3) in above lines corrected: IF(PR1>PR3) THEN; NEWQ=NQ-1;RH=PR1; ELSEIF(NQ^=MAXORD)THEN; NEWQ=LNQ;RH=PR3; -- PR3 IS GREATER _DO I=1,N;Y(I,NEWQ+1)=ERROR(I)*EL(LNQ)/DBLE(LNQ) _OD; ELSE; -- NQ==MAXORD NEWQ=NQ;RH=PR3; -- PR3 IS GREATER ENDIF; -- If there is a change of order, reset NQ and LNQ. -- In any case H is reset according to -- RH and the Y array and the coefficients EL() -- are rescaled. Then exit. NQ=NEWQ;LNQ=NQ+1; ENDIF; CALL RESCAL(Y,NYDIM,RH,RMAX,RC,LNQ); IDOUB=LNQ+1; %GXA3: SUM=0.D0; IF(METH^=3)THEN; _DO I=1,N;SUM=SUM+((ERROR(I)-Y(I,LMAX))/YMAX(I))**2 _OD; PR3=1.D0/(((SUM/EUP)**(.5D0/DBLE(LNQ+1)))*BIAS3+DFLTZR); ELSE; _DO I=1,N; SUM=SUM+ ((ERROR(I)*EL(1)/DBLE(NQ+2)-Y(I,LMAX)/TQ(3))/YMAX(I))**2 _OD; PR3=1.D0/(((SUM/EPSJ**2)**(.5D0/DBLE(LNQ+1)))*BIAS3+DFLTZR); ENDIF; %GXA2=GX2: %GXA1=GX1: %GY: t=tOLD;RMAX=RMXFAI; _DO J1=1,NQ; _DO J2=J1,NQ; J=NQ-J2+J1; _DO I=1,N; Y(I,J)=Y(I,J)-Y(I,J+1) _OD_OD_OD; IF(ABS(H)<=(HMIN*1.00001D0)) THEN; RETCOD=2; ELSE; write(*,'(a,i7,i3,1p,2e12.3)') ' GY: maxer, kflag, H, Hu:', maxer, kflag, H, Hu; RH=RHCORR; IF(METH==3)THEN; HNT(1)=H*RH; CALL COSET(METH,NQ,EL,TQ,MAXORD,IDOUB); <*A: the constants E,EUP,EDN,BND *>; RC=RC*EL(1)/OLDL0;OLDL0=EL(1); CALL RESCAL(Y,NYDIM,RH,RMAX,RC,LNQ); ELSE; CALL RESCAL(Y,NYDIM,RH,RMAX,RC,LNQ); IDOUB=LNQ+1; ENDIF; ENDIF; if(ifail==13)then; @wterm ' GY: ifail 13'; endif; %GYA=A: %G_COR_MatrOut: _DO JTST=1,N; --@VAR <*ROW: put in IRN(NZ+1:ILAST) all indeces I such that ICN(I)==JTST *>: _Do KTST=1,N; --@NAMELEFT,@NAMERIGHT; <*COMPR: IF in IRN(NZ+1:ILAST) there is I such that IRN(I)==KTST THEN output relevant vars *> _od; _OD %G_COR_MatrOut_ROW: ILAST=0; do I=1,NZMOD; IF(ICS(I)==JTST)THEN; ILAST=ILAST+1; IRN(NZ+ILAST)=I; ENDIF; enddo; -- WRITE(5,*)'ILAST=',ILAST; %G_COR_MatrOut_COMPR: FOUNDI=.FALSE.; IL=1; _WHILE IL<=ILAST & ^FOUNDI _DO IF(IRS(IRN(NZ+IL))==KTST)THEN; FOUNDI=.TRUE.; ELSE; IL=IL+1; ENDIF; _OD; IF(FOUNDI)THEN; If(KTST^=JTST)then; if(AJAC(IRN(NZ+IL))>1.d-6) WRITE(@Wres,'(A,2I7,1P,3E12.3)') ' K,J,AJAC : ', KTST,JTST,AJAC(IRN(NZ+IL)); else; WRITE(@Wres,'(A,2I7,1P,3E12.3)') ' K,J,AJAC,1-AJAC : ', KTST,JTST,AJAC(IRN(NZ+IL)),1.d0-AJAC(IRN(NZ+IL)); endif; ENDIF; <*PhysInfo : give info on variable number KTST *>; %G_COR_MatrOut_COMPR_PhysInfo: do I=1,NZMOD; if(IRS(I)==KTST .and. ICS(I)==JTST)then; WRITE(@Wres,'(3(A,1P,E12.4))') ' Y(K)=',Y(KTST,1),' Fsave(K)=',Fsave(KTST), ' Ydot(K)=',FSAVE(NYDIM+KTST); WRITE(@Wres,'(2(A,I5),A,1P,E12.4)') ' Control IR=',KTST, ' IC=',JTST,' AJAC=',AJAC(I); endif; enddo; -- this is for control of found AJAC; ICW=JTST; If(ICW<=NZON*NVARS)then; Izon=MOD(ICW-1,NZON)+1; Ivar=(ICW-1)/NZON +1; WRITE(@Wres,'(3A,I3)') ' Derivative over ',CharVarName(Ivar)(1:length(CharVarName(Ivar))), ' in zon=',Izon; else; ICW=ICW-NZON*NVARS; If(ICW<=KRAD)Then; L=(ICW-1)/(NZON-NCND)+1; Izon=NCND+MOD(ICW-1,(NZON-NCND))+1; WRITE(@Wres,'(2(A,I3))') ' over FJ in L=',L,' Izon=',Izon; else; ICW=ICW-KRAD; L=(ICW-1)/(NZON-NCND)+1; Izon=NCND+MOD(ICW-1,(NZON-NCND))+1; WRITE(@Wres,'(2(A,I3))') ' over FH in L=',L,' Izon=',Izon; endif; endif; If(KTST<=NZON*NVARS)then; Izon=MOD(KTST-1,NZON)+1; Ivar=(KTST-1)/NZON +1; WRITE(@Wres,'(3A,I3)') ' for Ydot of ',CharVarName(Ivar)(1:length(CharVarName(Ivar))), ' in zon=',Izon; else; KTST1=KTST-NZON*NVARS; If(KTST1<=KRAD)Then; L=(KTST1-1)/(NZON-NCND)+1; Izon=NCND+MOD(KTST1-1,(NZON-NCND))+1; WRITE(@Wres,'(2(A,I5))') ' Der. of dot(FJ) in L=',L,' Izon=',Izon; else; KTST1=KTST1-KRAD; L=(KTST1-1)/(NZON-NCND)+1; Izon=NCND+MOD(KTST1-1,(NZON-NCND))+1; WRITE(@Wres,'(2(A,I5))') ' Der. of dot(FH) in L=',L,' Izon=',Izon; endif; endif; %G_COR_DEB: _DEFINE @NOISE 1.D-4 -- for IBM Double Precision 3.D-8 -- 1.D-7 -- for Cray -- input variables: -- @NAME - zero starting position in FSAVE-array for testing -- @VAR - zero starting position in FSAVE-array for indep.var -- @NOISE- ignore derivatives less than @NOISE If(^BADSTE)Then; t=tOLD; _DO J1=1,NQ; _DO J2=J1,NQ; J=NQ-J2+J1; _DO I=1,N; Y(I,J)=Y(I,J)-Y(I,J+1) _OD_OD_OD; _Do Ich=1,Mzon; Chem0(Ich)=0.1; -- check if we need this _od; _Do I=1,N; FSAVE(I)=Y(I,1)_od; --************** ************** ************** ************** -- EVALJA=.FALSE.; -- never put it here: some errors may be lost --************** ************** ************** ************** Nperturb=0; -- index of perturbation, now NO perturbation CALL DIFMAT(envs, err); _Do KTST=1,N; -- save old YDOT DSAVE(KTST)=FSAVE(NYDIM+KTST); _od; _DO JTST=1,N; --@VAR -- _DO JTST=Ncnd+1,Ncnd+2; -- only d over dr is tested @wterm' JTST=',JTST; EVALJA=.FALSE.; OLD=FSAVE(JTST); If(JTST<Nzon)then; -- radius If(JTST==1)then; DeltY=DeltaDeb*min(FSAVE(JTST+1)-OLD,OLD-Rce); else; DeltY=DeltaDeb*min(FSAVE(JTST+1)-OLD,OLD-FSAVE(JTST-1)); endif; elseIf(Nzon<JTST & JTST<2*Nzon)then; -- velocity If(JTST==Nzon+1)then; DeltY=DeltaDeb*min(abs(FSAVE(JTST+1)-OLD),abs(OLD)); else; DeltY=DeltaDeb*min(abs(FSAVE(JTST+1)-OLD), abs(OLD-FSAVE(JTST-1))); endif; else; DeltY=OLD*DeltaDeb; endif; FSAVE(JTST)=OLD+DeltY; _Do Ich=1,Mzon; Chem0(Ich)=0.1 _od; Nperturb=1; -- index of perturbation, now IS perturbation CALL DIFMAT(envs, err); FSAVE(JTST)=OLD; <*ROW: put in IRN(NZ+1:ILAST) all indeces I such that ICN(I)==JTST *>; _Do KTST=1,N; --@NAMELEFT,@NAMERIGHT; @DNAMEDVAR(KTST) = (FSAVE(NYDIM+KTST)-DSAVE(KTST)); <*COMPR: IF in IRN(NZ+1:ILAST) there is I such that IRN(I)==KTST THEN IF AJAC(I) differs strongly from -HL*@DNAMEDVAR for KTST^=JTST or from 1.-HL*@DNAMEDVAR for KTST==JTST THEN output relevant vars ELSE IF @DNAMEDVAR is HIGHER than @NOISE THEN output relevant vars *> _od; _OD; STOP; else; WRITE(@Wres,*)' Testing not possible: BADSTE=',BADSTE; endif; %G_COR_DEB_ROW: ILAST=0; _DO I=1,NZMOD; IF(ICS(I)==JTST)THEN; ILAST=ILAST+1; IRN(NZ+ILAST)=I; ENDIF; _OD; -- WRITE(5,*)'ILAST=',ILAST; %G_COR_DEB_COMPR: FOUNDI=.FALSE.; IL=1; _WHILE IL<=ILAST & ^FOUNDI _DO IF(IRS(IRN(NZ+IL))==KTST)THEN; FOUNDI=.TRUE.; ELSE; IL=IL+1; ENDIF; _OD; IF(FOUNDI)THEN; If(KTST^=JTST)THEN; If( ABS(AJAC(IRN(NZ+IL))*DeltY +HL*@DNAMEDVAR(KTST)) > 1.D-1*ABS(AJAC(IRN(NZ+IL))*DeltY) & ABS(AJAC(IRN(NZ+IL))/HL) > @NOISE ) Then; If(abs(DeltY)>DFLTZR*abs(@DNAMEDVAR(KTST)))then; WRITE(@Wres,'(/A/2I7,1P,3G12.3)') ' K,J,AJAC,Stella df/dY,Num DN(K)/DV(J): ', KTST,JTST,AJAC(IRN(NZ+IL)), -AJAC(IRN(NZ+IL))/HL,@DNAMEDVAR(KTST)/DeltY; else; WRITE(@Wres,'(/A/2I7,1P,3G12.3)') ' K,J,AJAC,(df/dY,DN/DV)*DeltaDeb: ', KTST,JTST,AJAC(IRN(NZ+IL)), -AJAC(IRN(NZ+IL))*DeltY/HL,@DNAMEDVAR(KTST); endif; <*PhysInfo: give info on variable number KTST *>; endif; else; If( ABS((AJAC(IRN(NZ+IL))-1.D0)*DeltY+HL*@DNAMEDVAR(KTST)) > 1.D-1*ABS((AJAC(IRN(NZ+IL))-1.D0)*DeltY) & ABS((AJAC(IRN(NZ+IL))-1.D0)/HL) > @NOISE ) Then; If(abs(DeltY)>DFLTZR*abs(@DNAMEDVAR(KTST)))then; WRITE(@Wres,'(/A/2I7,1P,3G12.3)') ' K,J,AJAC,Stella df/dY,Num DN(K)/DV(J): ', KTST,JTST,AJAC(IRN(NZ+IL)),(1.D0 -AJAC(IRN(NZ+IL)) )/HL, @DNAMEDVAR(KTST)/DeltY ; else; WRITE(@Wres,'(/A/2I7,1P,3G12.3)') ' K,J,AJAC,(df/dY,DN/DV)*DeltaDeb: ', KTST,JTST,AJAC(IRN(NZ+IL)),(1.D0-AJAC(IRN(NZ+IL)))*DeltY/HL, @DNAMEDVAR(KTST); endif; <*PhysInfo2: give info on variable number KTST *>; endif; endif; ELSE; IF(ABS(@DNAMEDVAR(KTST)) > @NOISE) then; If(abs(DeltY)>DFLTZR*abs(@DNAMEDVAR(KTST)))then; WRITE(@Wres,'(/A,2I7,1P,2G12.3)')' NO JAC: K,J,DN(K)/DV(J) ', KTST,JTST, @DNAMEDVAR(KTST)/DeltY; else; WRITE(@Wres,'(/A,2I7,1P,2G12.3)')' NO JAC: K,J,DN(K)/DV(J)*DeltaDeb ', KTST,JTST, @DNAMEDVAR(KTST); endif; endif; ENDIF; %G_COR_DEB_COMPR_PhysInfo: do I=1,NZMOD; if(IRS(I)==KTST .and. ICS(I)==JTST)then; WRITE(@Wres,'(3(A,1P,E12.4))') ' Y(K)=',Y(KTST,1),' Fsave(K)=',Fsave(KTST), ' Ydot(K)=',FSAVE(NYDIM+KTST); WRITE(@Wres,'(2(A,I5),A,1P,E12.4)') ' Control IR=',KTST, ' IC=',JTST,' AJAC=',AJAC(I); endif; enddo; -- this is for control of found AJAC; ICW=JTST; If(ICW<=NZON*NVARS)then; Izon=MOD(ICW-1,NZON)+1; Ivar=(ICW-1)/NZON +1; WRITE(@Wres,'(3A,I3)') ' Derivative over ',CharVarName(Ivar)(1:length(CharVarName(Ivar))), ' in zon=',Izon; else; ICW=ICW-NZON*NVARS; If(ICW<=KRAD)Then; L=(ICW-1)/(NZON-NCND)+1; Izon=NCND+MOD(ICW-1,(NZON-NCND))+1; WRITE(@Wres,'(2(A,I3))') ' over FJ in L=',L,' Izon=',Izon; else; ICW=ICW-KRAD; L=(ICW-1)/(NZON-NCND)+1; Izon=NCND+MOD(ICW-1,(NZON-NCND))+1; WRITE(@Wres,'(2(A,I3))') ' over FH in L=',L,' Izon=',Izon; endif; endif; If(KTST<=NZON*NVARS)then; Izon=MOD(KTST-1,NZON)+1; Ivar=(KTST-1)/NZON +1; WRITE(@Wres,'(3A,I3)') ' for Ydot of ',CharVarName(Ivar)(1:length(CharVarName(Ivar))), ' in zon=',Izon; else; KTST1=KTST-NZON*NVARS; If(KTST1<=KRAD)Then; L=(KTST1-1)/(NZON-NCND)+1; Izon=NCND+MOD(KTST1-1,(NZON-NCND))+1; WRITE(@Wres,'(2(A,I5))') ' Der. of dot(FJ) in L=',L,' Izon=',Izon; else; KTST1=KTST1-KRAD; L=(KTST1-1)/(NZON-NCND)+1; Izon=NCND+MOD(KTST1-1,(NZON-NCND))+1; WRITE(@Wres,'(2(A,I5))') ' Der. of dot(FH) in L=',L,' Izon=',Izon; endif; endif; %G_COR_DEB_COMPR_PhysInfo2=G_COR_DEB_COMPR_PhysInfo: %X: -- _Case RETCOD _Of -- -- Control passes to _Esac for RETCOD==4 -- _1 KFLAG=-1 -- _2 KFLAG=-2 -- _3 RMAX=RMXNOR -- _Esac; select case(RETCOD); case(1); KFLAG=-1; case(2); KFLAG=-2; case(3); RMAX=RMXNOR; end select; HOLD=H; JSTART=NQ;