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;