SUBROUTINE RECEIV(XSOUR,YSOUR,ZSOUR,TSOUR,DT,AC,ITMAX, 1ASTART,ASTEP,AFIN,BMIN,BSTEP,BMAX,MOUT,LU1,LU2,METHOD,ITPR,NCD) C C TWO-POINT RAY TRACING C COMPLEX AMPX1,AMPX2,AMPY1,AMPY2,AMPZ1,AMPZ2,APX,APY,APZ DIMENSION JC(50,2),YDD(2),DEP(6) DIMENSION TIME(500),XCOOR(500),ZCOOR(500),INDI(500),AMPX1(500), 1AMPY1(500),AMPZ1(500),AMPX2(500),AMPY2(500),AMPZ2(500), 2p(500,3),pol(500,3),pol1(500,3),APX(2),APY(2),APZ(2) COMMON /AUXI/ IANI(20),INTR,INT1,IOUT,KRE,IREFR,LAY,NDER,IPRINT, 1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOU, 2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori INTEGER CODE COMMON /COD/ CODE(50,2),KREF,KC,ITYPE COMMON /DIST/ DST(200),NDST,REPS,PROF(2),NDSTP,PREPS,LNDST, 1xprf,yprf COMPLEX PS,CKMAH COMMON /RAY/ AY(28,400),DS(20,20),KINT(20),HHH(3,3),tmax, 1 PS(3,7,20),IS(8,20),DINC,DTOLD,N,IREF,IND,IND1 COMMON /RAY2/ DRY(3,400) COMMON/VSP/XVSP,YVSP,XNRM,YNRM,ICOD,IVSP COMMON/DYN/XDYN(3,3),ydyn(3,3) COMMON/KM/KMAH,SPROLD,TSTOLD C ITER=0 iprof=0 indr1=0 JNULL=0 II=0 DD=0. C AA=ASTART-ASTEP BMIN1=BMIN REPS1=.05*REPS INDEX=0 INUM=0 IK1=0 ICLS=0 ICR=0 INDS=ISOUR C C LOOP IN DECLINATION, FROM VALUE ASTART TO AFIN, WITH THE STEP C ASTEP C 1 AA=AA+ASTEP PNEW=AA IF(ASTEP.GT.0..AND.AA.GT.AFIN)GO TO 99 IF(ASTEP.LT.0..AND.AA.LT.AFIN)GO TO 99 4 IND=INDS CALL PROFIL(XSOUR,YSOUR,ZSOUR,TSOUR,AA,PAZM,r,XX,YY,Zz,T, 1DT,AC,BMIN,BSTEP,BMAX,ITMAX,MOUT,NCD,METHOD,ITPR,indr1) IF(IND.EQ.14.OR.IND.EQ.20)RETURN IF(ITPR.EQ.43)r=Zz IF(NDST.EQ.0.AND.LNDST.EQ.0)THEN BMIN=BMIN1 GO TO 1 END IF IF(NDST.EQ.0)GO TO 65 IF(IND.EQ.ITPR)rAX=r IF(IND.EQ.ITPR)PNEW=AA IF(IND.EQ.9)XCR=r IF(MOUT.GE.2)WRITE(LOU,100)IND,IND1,KMAH,r,T,AA IF(INUM.NE.0)GO TO 2 dOLD=AA aold=pazm IOLD=IND IOLD1=IND1 rOLD=r xold=xx yold=yy zold=zz TOLD=T INUM=1 GO TO 1 C C PARAMETERS FOR THE PRECEDING RAY: C DOLD (AA), AOLD (PAZM), ROLD (R), IOLD (IND) C PARAMETERS FOR THE NEW RAY: C DNEW (AA), ANEW (PAZM), RNEW (R), INEW (IND) C 2 INEW=IND INEW1=IND1 dNEW=AA anew=pazm rNEW=r xnew=xx ynew=yy znew=zz TNEW=T IF(INEW.EQ.ITPR.AND.IOLD.EQ.ITPR)GO TO 21 IF(INEW.EQ.ITPR)GO TO 50 IF(IOLD.EQ.ITPR)GO TO 55 IF(INEW.EQ.9.AND.IOLD.NE.9.AND.IOLD.NE.8)GO TO 30 IF(INEW.NE.9.AND.INEW.NE.8.AND.IOLD.EQ.9)GO TO 35 indr1=ind1 GO TO 3 21 IF(IOLD1.NE.INEW1)IK1=2 IF(IK1.EQ.2)GO TO 55 indr1=ind1 GO TO 40 C C NO ITERATIONS, TAKE A NEW RAY IN THE LOOP C 3 CONTINUE IOLD=INEW IOLD1=INEW1 rOLD=rNEW xold=xnew yold=ynew zold=znew dOLD=dNEW aold=anew TOLD=TNEW GO TO 1 C C REGULAR CASE: IOLD=3, INEW=3 C 40 rrNEW=rNEW xxnew=xnew yynew=ynew zznew=znew rrOLD=rOLD xxold=xold yyold=yold zzold=zold ddNEW=dNEW aanew=anew ddOLD=dOLD aaold=aold TTNEW=TNEW TTOLD=TOLD IBNEW=0 IBOLD=0 41 IF(rrNEW.GT.rrOLD.AND.DST(2).GT.DST(1))GO TO 46 IF(rrNEW.LT.rrOLD.AND.DST(2).LT.DST(1))GO TO 46 C C REGULAR CASE: RRNEW.LE.RROLD, ITREND=-1 (REVERSE BRANCH) C ITREND=-1 DO 42 I=1,NDST NDD=NDST-I+1 DD=DST(NDD) dr=rrOLD IF(IBOLD.EQ.1)dr=dr+REPS IF(DD.GE.dr)GO TO 42 dr=rrNEW IF(IBNEW.EQ.1)dr=dr-REPS IF(DD.LT.dr.AND.IK1.NE.0)GO TO 6 IF(DD.LT.dr)GO TO 3 II=NDD GO TO 43 42 CONTINUE IF(IK1.NE.0)GO TO 6 GO TO 3 C C REGULAR CASE: RRNEW.GT.RROLD, ITREND=1 (REGULAR BRANCH) C 46 ITREND=1 DO 44 I=1,NDST DD=DST(I) dr=rrOLD IF(IBOLD.EQ.1)dr=dr-REPS IF(DD.LE.dr)GO TO 44 dr=rrNEW IF(IBNEW.EQ.1)dr=dr+REPS IF(DD.GT.dr.AND.IK1.NE.0)GO TO 6 IF(DD.GT.dr)GO TO 3 II=I GO TO 43 44 CONTINUE IF(IK1.NE.0)GO TO 6 GO TO 3 43 d1=ddOLD a1=aaold d2=ddNEW a2=aanew x1=xxold y1=yyold z1=zzold x2=xxnew y2=yynew z2=zznew T1=TTOLD T2=TTNEW C C I T E R A T I O N S C 60 continue C C TRANSFORMATION OF COORDINATES OF A RECEIVER FROM CYLINDRICAL C TO CARTESIAN COORDINATES C IF(ITPR.NE.43)THEN XD=Xprf+DD*COS(PROF(1)) YD=Yprf+DD*SIN(PROF(1)) YDD(1)=XD YDD(2)=YD INTR=1 IF(ITPR.GT.100)INTR=ITPR-100 CALL DISC(YDD,DEP) ZD=DEP(1) END IF IF(ITPR.EQ.43)THEN XD=XVSP YD=YVSP ZD=DD END IF DELX=XD-X1 DELY=YD-Y1 DELZ=ZD-Z1 dr1=sqrt(delx*delx+dely*dely+delz*delz) DELX=XD-X2 DELY=YD-Y2 DELZ=ZD-Z2 dr2=sqrt(delx*delx+dely*dely+delz*delz) c ITER=0 ISIGN=1 IPR1=0 IPR2=0 GO TO 61 68 rAX=r PAX=PNEW 61 ITER=ITER+1 IF(ITER.GT.ITMAX)GO TO 80 ISIGN=-ISIGN C C PREPARATION FOR ITERATIVE SOLUTION OF TWO-POINT RAY TRACING C daux=0.5*(d1+d2) azaux=0.5*(a1+a2) if(method.eq.2.or.iter.eq.1.or.itpr.ne.ind)go to 69 C C PARAXIAL RAY APPROXIMATION C AUX1=XDYN(2,1)*XDYN(3,2)-XDYN(3,1)*XDYN(2,2) AUX2=XDYN(3,1)*XDYN(1,2)-XDYN(1,1)*XDYN(3,2) AUX3=XDYN(1,1)*XDYN(2,2)-XDYN(2,1)*XDYN(1,2) DET=AUX1*XDYN(1,3)+AUX2*XDYN(2,3)+AUX3*XDYN(3,3) IF(ABS(DET).LT..00001)GO TO 69 AUX1=DELY*XDYN(3,2)-DELZ*XDYN(2,2) AUX2=DELZ*XDYN(1,2)-DELX*XDYN(3,2) AUX3=DELX*XDYN(2,2)-DELY*XDYN(1,2) aux11=(AUX1*XDYN(1,3)+AUX2*XDYN(2,3)+AUX3*XDYN(3,3))/DET AUX1=DELZ*XDYN(2,1)-DELY*XDYN(3,1) AUX2=DELX*XDYN(3,1)-DELZ*XDYN(1,1) AUX3=DELY*XDYN(1,1)-DELX*XDYN(2,1) aux22=(AUX1*XDYN(1,3)+AUX2*XDYN(2,3)+AUX3*XDYN(3,3))/DET if(mori.eq.0)pazm=pazm+aux11 if(mori.ne.0)pazm=pazm+aux22 if(mori.eq.0)pnew=pnew+aux22 if(mori.ne.0)pnew=pnew+aux11 go to 72 C C HALVING OF INTERVAL C 69 pnew=0.5*(d1+d2) pazm=0.5*(a1+a2) C C INITIAL ANGLES FOR A NEW RAY WERE DETERMINED C 72 ind=inds NOLD=0 rOLD=0. SPROLD=0. IF(MDIM.GE.1)NDER=2 CALL RAYA(XSOUR,YSOUR,ZSOUR,TSOUR,PNEW,PAZM,PX,PY,PZ,XX,YY,Zz, 1T,DT,AC) NDER=1 XE=XX-Xprf YE=YY-Yprf r=SQRT(XE*XE+YE*YE) delx=xd-xx dely=yd-yy delz=zd-zz drs=sqrt(delx*delx+dely*dely+delz*delz) IF(IND.EQ.20)RETURN IF(ITPR.EQ.43)r=Zz IF(MOUT.GE.2)WRITE(LOU,101)IND,IND1,ITER,KMAH,DD,r,T,PNEW,PAZM if(mout.eq.4)write(lou,120)xd,yd,zd,xx,yy,zz 120 format(1x,'(x,y,z) receiver',3F15.8/1x,'(xx,yy,zz) ray',3F15.8) IF(ICLS.NE.0)GO TO 70 IF(IND.NE.ITPR)then d2=PNEW a2=pazm x2=xx y2=yy z2=zz dr2=drs GO TO 61 end if IF(drs.LT.REPS)GO TO 65 IF(dr2.LT.dr1)GO TO 63 d2=PNEW a2=pazm x2=xx y2=yy z2=zz dr2=drs T2=T IPR2=1 GO TO 68 63 continue d1=PNEW a1=pazm x1=xx y1=yy z1=zz dr1=drs T1=T IPR1=1 GO TO 68 70 ICLS=0 C C SUCCESSFUL ITERATIONS C 65 INDEX=INDEX+1 IF(LU1.EQ.0.OR.NDST.EQ.0)GO TO 800 TIME(INDEX)=T XCOOR(INDEX)=r ZCOOR(INDEX)=Zz INDI(INDEX)=II 800 CONTINUE rAX=r IF(MOUT.GE.1)WRITE(LOU,113)DD,r,XX,YY,Zz,T,PNEW,PAZM, 1IND,IND1,ITER,II,INDEX IF(LU1.NE.0.AND.NDST.NE.0)WRITE(LU1,100)N,II IF(LU1.NE.0.AND.NDST.NE.0) 1WRITE(LU1,108)(AY(2,I),AY(3,I),AY(4,I),I=1,N) IF(MDIM.EQ.0)GO TO 80 CALL AMPL(APX,APY,APZ,UU) UU=1. SPR=1. CKMAH=(1.,0.) IF(MDIM.EQ.2)THEN IF(KMAH.NE.0)THEN PH=-1.57079632*KMAH CSKMAH=COS(PH) SNKMAH=SIN(PH) CKMAH=CMPLX(CSKMAH,SNKMAH) END IF SPR1=XDYN(2,1)*XDYN(3,2)-XDYN(3,1)*XDYN(2,2) SPR2=XDYN(3,1)*XDYN(1,2)-XDYN(1,1)*XDYN(3,2) SPR3=XDYN(1,1)*XDYN(2,2)-XDYN(2,1)*XDYN(1,2) SPRv=XDYN(1,3)*SPR1+XDYN(2,3)*SPR2+XDYN(3,3)*SPR3 IF(NDST.EQ.0.and.mori.eq.0)CSF=COS(AA) IF(NDST.NE.0.and.mori.eq.0)CSF=COS(PNEW) if(mori.ne.0)csf=cos(pazm) if(mori.eq.0)sprv=-sprv VV=ay(5,n)*ay(5,n)+ay(6,n)*ay(6,n)+ay(7,n)*ay(7,n) SPR=SPRv*SQRT(VV) SPR=SQRT(ABS(SPR/csf)) IF(MOUT.GE.2)WRITE(LOU,110)XDYN IF(MOUT.GE.2)WRITE(LOU,112)yDYN END IF DO 802 I=1,2 APX(I)=APX(I)*UU*CKMAH/SPR APY(I)=APY(I)*UU*CKMAH/SPR APZ(I)=APZ(I)*UU*CKMAH/SPR 802 CONTINUE IF(MOUT.GE.1.AND.MDIM.EQ.2) 1WRITE(LOU,'(2X,F8.5,6(2X,E11.5)/10X,6(2X,E11.5),F10.5,I5)') 2UU,(APX(I),APY(I),APZ(I),I=1,2),SPR,KMAH TAUX=T TAST=0. NCC=code(1,2) ncod=ncd IF(iani(isour).eq.0.and.ncc.eq.1)NCOD=-NCD call polar(1,1,1,1) IF(LU2.NE.0.AND.NDST.NE.0)then WRITE(LU2,116)ncod,II,T,APX(1),APY(1),APZ(1),TAST if(ncc.eq.1.and.ncod.lt.0)WRITE(LU2,115)APX(2),APY(2),APZ(2) WRITE(LU2,114)ay(5,1),ay(6,1),ay(7,1) if(ncc.eq.1)WRITE(LU2,114)(hhh(1,i),i=1,3) if(ncc.eq.1.and.ncod.lt.0)WRITE(LU2,114)(hhh(2,i),i=1,3) if(ncc.eq.2)WRITE(LU2,114)(hhh(2,i),i=1,3) if(ncc.eq.3)WRITE(LU2,114)(hhh(3,i),i=1,3) end if IF(LU1.EQ.0.OR.NDST.EQ.0)GO TO 801 AMPX1(INDEX)=APX(1) AMPY1(INDEX)=APY(1) AMPZ1(INDEX)=APZ(1) if(ncc.eq.1.and.ncod.lt.0)then AMPX2(INDEX)=APX(2) AMPY2(INDEX)=APY(2) AMPZ2(INDEX)=APZ(2) end if p(index,1)=ay(5,1) p(index,2)=ay(6,1) p(index,3)=ay(7,1) if(ncc.eq.1)then pol(index,1)=hhh(1,1) pol(index,2)=hhh(1,2) pol(index,3)=hhh(1,3) end if if(ncc.eq.1.and.ncod.lt.0)then pol1(index,1)=hhh(2,1) pol1(index,2)=hhh(2,2) pol1(index,3)=hhh(2,3) end if if(ncc.eq.2)then pol(index,1)=hhh(2,1) pol(index,2)=hhh(2,2) pol(index,3)=hhh(2,3) end if if(ncc.eq.3)then pol(index,1)=hhh(3,1) pol(index,2)=hhh(3,2) pol(index,3)=hhh(3,3) end if 801 CONTINUE C C 80 IF(NDST.EQ.0.AND.LNDST.EQ.1)THEN BMIN=BMIN+BSTEP GO TO 4 END IF d1=PNEW a1=pazm dr1=drs x1=xx y1=yy z1=zz T1=TAUX d2=ddNEW a2=aanew x2=xxnew y2=yynew z2=zznew T2=TTNEW IF(ITREND.EQ.(-1))GO TO 85 II=II+1 IF(II.GT.NDST.AND.IK1.NE.0)GO TO 6 IF(method.eq.1)then aa=pnew inum=0 end if IF(II.GT.NDST)GO TO 3 DD=DST(II) if(method.eq.1)go to 60 dr=rrNEW IF(IBNEW.EQ.1)dr=dr+REPS IF(DD.GT.dr.AND.IK1.NE.0)GO TO 6 IF(DD.GT.dr)GO TO 3 GO TO 60 85 II=II-1 IF(II.LT.1.AND.IK1.NE.0)GO TO 6 IF(method.eq.1)then aa=pnew inum=0 end if IF(II.LT.1)GO TO 3 DD=DST(II) if(method.eq.1)go to 60 dr=rrNEW IF(IBNEW.EQ.1)dr=dr-REPS IF(DD.LT.dr.AND.IK1.NE.0)GO TO 6 IF(DD.LT.dr)GO TO 3 GO TO 60 C C 6 CONTINUE IF(IK1.EQ.1)GO TO 7 IK1=1 d1=dNEW a1=anew d2=ddNEW a2=aanew IOLD1=INEW1 indr1=iold1 GO TO 59 7 IK1=0 rrOLD=rrNEW xxold=xxnew yyold=yynew zzold=zznew ddOLD=ddNEW aaold=aanew TTOLD=TTNEW IBOLD=IBNEW rrNEW=rNEW xxnew=xnew yynew=ynew zznew=znew ddNEW=dNEW aanew=anew TTNEW=TNEW IBNEW=0 GO TO 41 C C E N D O F I T E R A T I O N S C C C BOUNDARY RAYS. CASE IOLD.NE.3, INEW=3 C 50 rrNEW=rNEW xxnew=xnew yynew=ynew zznew=znew TTNEW=TNEW ddNEW=dNEW aanew=anew IBNEW=0 d1=dOLD d2=dNEW a1=dOLD a2=dNEW 54 IRES=0 ITER=0 51 PNEW=0.5*(d1+d2) ITER=ITER+1 IND=INDS CALL PROFIL(XSOUR,YSOUR,ZSOUR,TSOUR,PNEW,PAZM,r,XX,YY,Zz,T, 1DT,AC,BMIN,BSTEP,BMAX,ITMAX,MOUT,NCD,METHOD,ITPR,indr1) IF(IND.EQ.20)RETURN IF(ITPR.EQ.43)r=Zz IF(MOUT.GE.2)WRITE(LOU,103)IND,IND1,ITER,r,T,PNEW IF(IND.EQ.ITPR.AND.LNDST.EQ.1)GO TO 52 d1=PNEW a1=pazm GO TO 53 52 rrOLD=r xxold=xx yyold=yy zzold=zz ddOLD=PNEW aaold=pazm TTOLD=T IBOLD=1 IF(ABS(r-rAX).LT.REPS1.AND.IRES.EQ.1)ITER=ITMAX IRES=1 rAX=r ZAX=Zz IAX=IND IAX1=IND1 T1=T d2=PNEW a2=pazm 53 IF(ITER.LT.ITMAX)GO TO 51 IF(MOUT.GE.1.AND.IRES.EQ.1) 1WRITE(LOU,107)rAX,ZAX,T1,d2,IAX,IAX1,IRES IF(MOUT.GE.1.AND.IRES.EQ.0) 1WRITE(LOU,107)r,Zz,T,PNEW,IND,IND1,IRES IF(IRES.EQ.1)GO TO 41 GO TO 3 C C BOUNDARY RAYS. CASE IOLD=3, INEW.NE.3 C OR CASE IOLD=3, INEW=3 BUT IOLD1.NE.INEW1 C 55 rrOLD=rOLD xxold=xold yyold=yold zzold=zold ddOLD=dOLD aaold=aold TTOLD=TOLD IBOLD=0 d1=dOLD d2=dNEW a1=aOLD a2=aNEW 59 IRES=0 ITER=0 56 PNEW=0.5*(d1+d2) ITER=ITER+1 IND=INDS CALL PROFIL(XSOUR,YSOUR,ZSOUR,TSOUR,PNEW,PAZM,r,XX,YY,Zz,T, 1DT,AC,BMIN,BSTEP,BMAX,ITMAX,MOUT,NCD,METHOD,ITPR,indr1) IF(IND.EQ.20)RETURN IF(ITPR.EQ.43)r=Zz IF(MOUT.GE.2)WRITE(LOU,103)IND,IND1,ITER,r,T,PNEW IF(IND.EQ.ITPR.AND.IBOLD.EQ.1.AND.LNDST.EQ.1)GO TO 57 IF(IND.EQ.ITPR.AND.IND1.EQ.IOLD1.AND.LNDST.EQ.1)GO TO 57 d2=PNEW a2=pazm GO TO 58 57 rrNEW=r xxnew=xx yynew=yy zznew=zz ddNEW=PNEW aanew=pazm TTNEW=T IBNEW=1 IF(ABS(r-rAX).LT.REPS1.AND.IRES.EQ.1)ITER=ITMAX IRES=1 rAX=r ZAX=Zz IAX=IND IAX1=IND1 T2=T d1=PNEW a1=pazm 58 IF(ITER.LT.ITMAX)GO TO 56 IF(MOUT.GE.1.AND.IRES.EQ.1) 1WRITE(LOU,107)rAX,ZAX,T2,d1,IAX,IAX1,IRES IF(MOUT.GE.1.AND.IRES.EQ.0) 1WRITE(LOU,107)r,Zz,T,PNEW,IND,IND1,IRES IF(IRES.EQ.1.AND.IK1.EQ.1)GO TO 7 IF(IRES.EQ.1)GO TO 41 GO TO 3 C C CRITICAL RAYS. CASE IOLD.NE.9, IOLD.NE.3, INEW=9 C 30 ITER=0 XCR=rNEW d1=dOLD d2=dNEW IRES=0 31 PNEW=0.5*(d1+d2) ITER=ITER+1 IND=INDS CALL PROFIL(XSOUR,YSOUR,ZSOUR,TSOUR,PNEW,PAZM,r,XX,YY,Zz,T, 1DT,AC,BMIN,BSTEP,BMAX,ITMAX,MOUT,NCD,METHOD,ITPR,indr1) IF(IND.EQ.20)RETURN IF(ITPR.EQ.43)r=Zz IF(MOUT.GE.2)WRITE(LOU,104)IND,IND1,ITER,r,T,PNEW IF(IND.EQ.9)GO TO 32 IF(IND.EQ.ITPR)GO TO 33 d1=PNEW GO TO 34 32 IF(IND1.NE.INEW1)d1=PNEW IF(IND1.NE.INEW1)GO TO 34 d2=PNEW IF(ITER.EQ.ITMAX.AND.KC.NE.0.AND.IRES.EQ.1)GO TO 89 XCR=r GO TO 34 89 DO 86 K=1,KREF DO 86 L=1,2 86 JC(K,L)=CODE(K,L) IRF3=IREF+3 DO 87 K=IRF3,KREF DO 87 L=1,2 87 CODE(K-2,L)=JC(K,L) KREF1=KREF KREF=KREF-2 ICR=1 ITER=ITMAX-1 GO TO 31 33 IF(ABS(r-rAX).LT.REPS1.AND.IRES.EQ.1)ITER=ITMAX IRES=1 rAX=r rrnew=r xxnew=xx yynew=yy zznew=zz ZAX=Zz IAX=IND IAX1=IND1 T2=T d1=PNEW PAP=PNEW 34 IF(ITER.LT.ITMAX)GO TO 31 IF(MOUT.GE.1.AND.IRES.EQ.1) 1WRITE(LOU,111)rAX,ZAX,T2,PAP,IAX,IAX1,IRES IF(MOUT.GE.1.AND.IRES.EQ.0) 1WRITE(LOU,111)r,Zz,T,PNEW,IND,IND1,IRES IF(IRES.EQ.0)GO TO 3 d2=PAP rrNEW=rAX ddNEW=d2 aanew=anew TTNEW=T2 IBNEW=1 d1=dOLD IF(ICR.EQ.0)GO TO 54 ICR=0 KREF=KREF1 DO 88 K=1,KREF DO 88 L=1,2 88 CODE(K,L)=JC(K,L) GO TO 54 C C CRITICAL RAYS. CASE IOLD=9, INEW.NE.9, INEW.NE.3. C 35 ITER=0 XCR=rOLD d1=dOLD d2=dNEW IRES=0 36 PNEW=0.5*(d1+d2) ITER=ITER+1 IND=INDS CALL PROFIL(XSOUR,YSOUR,ZSOUR,TSOUR,PNEW,PAZM,r,XX,YY,Zz,T, 1DT,AC,BMIN,BSTEP,BMAX,ITMAX,MOUT,NCD,METHOD,ITPR,indr1) IF(IND.EQ.20)RETURN IF(ITPR.EQ.43)r=Zz IF(MOUT.GE.2)WRITE(LOU,104)IND,IND1,ITER,r,T,PNEW IF(IND.EQ.9)GO TO 37 IF(IND.EQ.ITPR)GO TO 38 d2=PNEW GO TO 39 37 IF(IND1.NE.INEW1)d2=PNEW IF(IND1.NE.INEW1)GO TO 39 d1=PNEW IF(ITER.EQ.ITMAX.AND.KC.NE.0.AND.IRES.EQ.1)GO TO 94 XCR=r GO TO 39 94 DO 91 K=1,KREF DO 91 L=1,2 91 JC(K,L)=CODE(K,L) IRF3=IREF+3 DO 92 K=IRF3,KREF DO 92 L=1,2 92 CODE(K-2,L)=JC(K,L) KREF1=KREF KREF=KREF-2 ICR=1 ITER=ITMAX-1 GO TO 36 38 IF(ABS(r-rAX).LT.REPS1.AND.IRES.EQ.1)ITER=ITMAX IRES=1 rAX=r rrold=r xxold=xx yyold=yy zzold=zz ZAX=Zz IAX=IND IAX1=IND1 d2=PNEW PAP=PNEW T1=T 39 IF(ITER.LT.ITMAX)GO TO 36 IF(MOUT.GE.1.AND.IRES.EQ.1) 1WRITE(LOU,111)rAX,ZAX,T1,PAP,IAX,IAX1,IRES IF(MOUT.GE.1.AND.IRES.EQ.0)WRITE(LOU,111)r,Zz,T,PNEW,IND,IND1,IRES IF(IRES.EQ.0)GO TO 3 d1=PAP rrOLD=rAX ddOLD=d1 aaold=aold TTOLD=T1 IBOLD=1 d2=dNEW IF(ICR.EQ.0)GO TO 59 ICR=0 KREF=KREF1 DO 93 K=1,KREF DO 93 L=1,2 93 CODE(K,L)=JC(K,L) GO TO 59 C C 100 FORMAT(3I3,2F10.5,F15.10) 101 FORMAT(1X,'ITERATIONS',5X,4I3,3F10.5,2F15.10) 102 FORMAT(5X,6I3,3F10.5,F15.10) 103 FORMAT(2X,'BOUNDARY RAY',5X,3I3,2F10.5,F15.10) 104 FORMAT(2X,'CRITICAL RAY',5X,3I3,2F10.5,F15.10) 107 FORMAT(10X,3F10.5,F15.10,3I3,5X,'BOUNDARY RAY') 108 FORMAT(3E15.5) 109 FORMAT(I5,9E15.5) 110 FORMAT(1X,'XDYN',3F10.5) 111 FORMAT(10X,3F10.5,F15.10,3I3,5X,'CRITICAL RAY') 112 FORMAT(1X,'YDYN',3F10.5) 113 FORMAT(7F10.5,F15.10,5I3) 114 FORMAT(3F10.5) 115 format(6e12.6) 116 FORMAT(2I3,F10.5,12E12.6,F10.6) 117 FORMAT(6E15.5) C C 99 N=0 NAUX=0 IF(LU1.NE.0.AND.NDST.NE.0)WRITE(LU1,100)N,NAUX IF(ncc.eq.1.and.ncod.lt.0)INDEX=-INDEX IF(LU1.NE.0.AND.NDST.NE.0)WRITE(LU1,100)INDEX IF(INDEX.EQ.0)RETURN IF(INDEX.LT.0)INDEX=-INDEX IF(LU1.NE.0.AND.NDST.NE.0)then do 200 i=1,index WRITE(LU1,109)INDI(I),XCOOR(I),ZCOOR(I),TIME(I), 1 AMPX1(I),AMPY1(I),AMPZ1(I) if(ncc.eq.1.and.ncod.lt.0) 1 write(lu1,117)AMPX2(I),AMPY2(I),AMPZ2(I) write(lu1,108)(p(i,k),k=1,3) write(lu1,108)(pol(i,k),k=1,3) if(ncc.eq.1.and.ncod.lt.0) 1 write(lu1,108)(pol1(i,k),k=1,3) 200 continue end if RETURN END