79900 IF(ISTART.EO.IEND) 1= I END-I I+IBEGIN 



80000 IFdSTART.EO.IEND AND. I.EO.IEND) GO TO 510 



80100 C*ADX EQUALS ACTUAL DELTA X ACROSS SPACE STEP. 



80200 C*ONLY ON BOUNDARIES WHERE FORWARD OR BACKWARD DIFFERENCING. 



80300 IF(I .NE. IBEGIN) GO TO 503 



80400 ADX=DX 



80500 IP=I+1 



80600 IM=I 



80700 GO TO 505 



80800 503 IF(I .NE . lEND) GO TO 504 



80900 ADX=DX 



81000 IP=I 



81100 IM=I-1 



81200 GO TO 505 



81300 504 ADX = 2.0*DX 



81400 IP=I+1 



81500 IM=I-1 



81600 505 CONTINUE 



81700 DO 502 J = JBEGIN( I ) , dEND( I ) - 1 



81800 Jd=JEND(I )-1-J+JBEGIN(I) 



81900 HOLDd , Jd)=H(I , JJ) 



82000 YBAR=0.25*(Y(I ,dJ)+2.0*Y(I , dd+ 1 )+Y(I,dd+2)) 



82100 CALL LOC( IM.dd.dOIM.dSIM, YBAR, IMINUS) 



82200 CALL LOC( IP , dd , dOIP , dSIP , YBAR , IPLUS ) 



82300 PART13 = (H{I , dd+ 1 )**2. )*CG(I ,dd+1 ) *COS( THETA( I , dd+ 1 )) 



82400 PART2=DY(I ,dd)/ADX 



82500 IF(dSIM.NE.O) GO TO 311 



82600 PART4B-0.0 



82700 GO TO 312 



82800 31 1 TOPIMH=(H(IM.dOIM-1 )**2. ) *CG( IM, dOIM- 1 ) *( SIN( THETA(IM , dOIM- 1 ) ) ) 



82900 B0TIMH=(H(IM,dSIM)**2. ) *CG( IM, dSIM) *SIN( THETA(IM . dSIM) ) 



83000 TOTALB = 0.5*(Y(IM,dOIM) + Y(IM,dOIM-1 ) ) -0 . 5*( Y ( IM, dSIM+ 1) + Y ( IM , dSIM ) ) 



83100 DUMB=0.5*(Y(IM.dOIM)+Y(IM,dOIM-1 ))-YBAR 



83200 PART4B=( (TOTALB-DUMB)*(TOPIMH-BOTIMH)/TOTALB)+BOTIMH 



83300 312 IF(dSIP.NE .0) GO TO 313 



83400 PART4A=0.0 



83500 GO TO 314 



83600 313 TOPIPH=(H(IP,dOIP-1)**2. ) ♦CGC IP , dOI P- 1 ) *SIN( THETA( IP , dOIP- 1 ) ) 



83700 B0TIPH=(H(IP.d3IP)**2. ) *CG( IP , dSIP ) *SIN( THETA(IP , dSIP ) ) 



83800 T0TALA = O.5*(Y(IP,d0IP) + Y(IP,d0IP-1 ) ) -O. 5* ( Y( IP , dSIP+ 1 ) + Y( IP , dSIP ) ) 



83900 DUMA=O.5*(Y(IP,d0IP)+Y(IP,d0IP-1 ) )-YBAR 



84000 PART4A=( ( TOTALA-DUMA ) * ( TOPIPH-BOTIPH)/TOTALA )+BOTIPH 



84100 314 PART4=PART4A-PART4B 



84200 YBARP=0.25*(Y(I , Jd+1)+2. *Y(I ,dd+2)+Y(I .dd+3)) 



84300 CALL LOC( IM , dd+ 1 , dPOIM , dPSIM , YBARP , IMINUS ) 



84400 CALL LOC( IP , dd+ 1 , dPOIP , dPSIP , YBARP , IPLUS ) 



84500 IF(dPSIM.NE .0) GO TO 315 



84600 PART12=0.0 



84700 GO TO 316 



84800 315 T0PMH=(H(IM,dP0IM-1 ) * +2 ) *CG( IM , dPOIM- 1 ) *COS( THETA( IM, dPOIM- 1 )) 



84900 B0TMH=(H(IM,dPSIM)**2)*CG(IM,dPSIM)*C0S(THETA(IM,dPSIM)) 



85000 TOTB= . 5* ( Y ( IM, dPOIM)+Y ( IM,dP0IM-1 ) )- . 5* ( Y ( IM, dPSIM+ 1 )+Y( IM , dPSIM) ) 



85100 DUMPB=0.5*( Y( IM.dPOIM)+Y(IM,dPOIM- 1 ) )-YBARP 



85200 PART 1 2 = ( ( TOTB-DUMPB ) * ( TOPMH-BOTMH ) /TOTB ) + BOTMH 



85300 316 IF(dPSIP.NE .0) GO TO 317 



85400 PARTI 1=0.0 



85500 GO TO 318 



85600 317 T0PPH=(H(IP,dP0IP-1 ) **2 ) *CG( IP , dPOIP- 1 ) *COS( THETA( IP , dPOIP- 1 )) 



85700 B0TPH=(H(IP,dPSIP)**2)*CG(IP,dPSIP)*C0S(THETA(IP,dPSIP)) 



85800 T0TA=.5*(Y(IP,dP0IP)+Y(IP.dP0IP-1 ))- . 5* ( Y( IP . dPSIP+ 1 )+Y( IP , dPSIP ) ) 



85900 DUMPA=O.5*(Y(IP,dP0IP)+Y(IP.dP0IP-1 )) -YBARP 



86000 PARTI 1 = ( (TOTA-DUMPA)*(TOPPH-BOTPH)/TOTA)+BOTPH 



86100 318 PART1H=TAU*PART12+( 1 . -2. *TAU)*PART13 + TAU*PART1 1 



86200 IF(dPSIM. E0.O)PART1H=( 1 . -TAU ) *PART 13 + TAU*PART 1 1 



86300 IF(dPSIP.E0.O)PART1H=TAU*PART12+( 1 . -TAU)*PART13 



86400 ARG=((PART1H+PART2*PART4)/(CG(I , dd ) *COS( THETA( I , dd ) ) ) ) 



86500 C*IF THERE IS TO BE AN INVALID SORT. USE LINEAR SHOALING. 



86600 IF(ARG.GE .0. ) GO TO 44 



86700 ARG=(CG(I , dd+ 1 ) *COS( THETA( I . dd+ 1 )))/(CG(I , dd) *COS(THETA(I , dd ) ) ) 



86800 IF(ARG.LT.O.O) ARG=0.0 



86900 H( I ,dd)=H( I ,dd+1 )*SORT(ARG) 



87000 GO TO 45 



77 



