Stella

radiative hydro
code stella

stradsep5tt.trf

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

begradron.trf

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

mStiffbgh.trf

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;