C DATA SET MEDLINE AT LEVEL 002 AS OF 03/19/81 SUBROUTINE CENTER(ALAT,ALONG, CLAT,CLONG, DISTMN) 00001 C*** .GIVEN. 3 GEO. PTS .. (ALAT(I), ALONG(I)), I=1,3) 00002 C .RETURNS. CENTER OF GEODESIC CIRCLE ON SPHEROID .. (CLAT,CLONG) 00003 C AND RADIUS, DISTMN 00004 C .VIA. INITIAL APPROX"N, AND ITERATION 00005 IMPLICIT REAL*8(A-H,O-Z), INTEGER(I-N) 00006**2 COMMON/IOFILS/ IN, LP 00007 COMMON/SPHDTA/ A,B, ESQ 00008 DIMENSION ALAT(3),ALONG(3),DIST(3),AZ(3) 00009 LOGICAL EQPT 00010 EQPT(X1,Y1, X2,Y2) = X1.EQ.X2 .AND. Y1.EQ.Y2 00011 OFFZRO(U)=DSIGN(DMAX1(DABS(U),1.D-9),U) 00012**2 TRIM(U,ELOW,HIGH) = DMAX1(ELOW,DMIN1(U,HIGH)) 00013**2 SIND(U)=DSIN(U*RADANS) 00014**2 COSD(U) = DCOS(U*RADANS) 00015**2 SQRT(U)=DSQRT(U) 00016**2 DATA RADANS/.017453292519943/, ERRMAX/.01/,KNTMAX/20/ 00017 C 00018 C** TEST FOR 2 OF 3 PTS THE SAME 00019 I = 0 00020 IF( EQPT(ALAT(1),ALONG(1), ALAT(3),ALONG(3)) ) I = 2 00021 IF( EQPT(ALAT(2),ALONG(2), ALAT(3),ALONG(3)) ) I = 1 00022 IF( I.NE.0 ) GO TO 80 00023 C*** COMPUTE INITIAL ESTIMATE = CENTER OF TRIANGLE ON SPHERE 00024 CALL CSPHTR(ALAT,ALONG, CLAT,CLONG) 00025 C*** ITERATE ... 00026 C*** .GIVEN. APPROXIMATE CENTER (CLAT,CLONG), 00027 C AND BEARINGS, AZ(I), FROM CENTER-PT TO 3 PTS, 00028 C MOVE DS FROM CENTER, AWAY FROM PT ( DX EAST/WEST, DY NORTH/SOUTH )00029 C ALONG EACH GEODESIC RADIUS, SO THAT 3 RADII BECOME EQUAL ... 00030 C I. E. 00031 C (S + DS)(I) = CONSTANT, I = 1, 2, 3 00032 C .WHERE. -DS = COS(AZ)*DX + SIN(AZ)*DY 00033 C 00034 C SUBTRACTING EQN. FOR PT 3 FROM EQN. FOR PT 1 ... 00035 C S1 - S3 = (COS(AZ1) - COS(AZ3))*DX + (SIN(AZ1) - SIN(AZ3))*DY 00036 C SUBTRACTING EQN. FOR PT 3 FROM EQN. FOR PT 2 ... 00037 C S2 - S3 = (COS(AZ2) - COS(AZ3))*DX + (SIN(AZ2) - SIN(AZ3))*DY 00038 C 00039 C I.E., 2 EQUATIONS IN 2 UNKNOWNS: DX AND DY 00040 C AND SOLVE VIA CRAMERS" RULE 00041 C 00042 C THEN DLAT = DX/RHO, AND DLONG = -DY/(NU * COS(CLAT)) 00043 C WHERE THE RADII OF CURVATURE, RHO AND NU (ACROSS / ALONG MERIDIAN)00044 C ARE APPROXIMATED BY NU AT LATITUDE CLAT, FROM SOLUTION ON SPHERE00045 C**** 00046 RADIUS = A / SQRT(1. - ESQ*SIND(CLAT)**2) 00047 DO 39 KOUNT = 1, KNTMAX 00048 C** FIND DISTANCES AND AZIMUTHS. 00049 DO 31 I = 1, 3 00050 31 CALL GEOINV(CLAT,CLONG, ALAT(I),ALONG(I), DIST(I),AZ(I),TEMP) 00051 C** FIND ERRORS IN DISTANCES. IF SMALL ENOUGH, QUIT. 00052 C1 = DIST(1) - DIST(3) 00053 C2 = DIST(2) - DIST(3) 00054 DISTMN=DMIN1(DIST(1),DIST(2),DIST(3)) 00055**2 IF(DMAX1(DABS(C1),DABS(C2)).LE.ERRMAX) RETURN 00056**2 C** ERRORS ARE TOO BIG. ITERATE. 00057 SINAZ3 = SIND(AZ(3)) 00058 COSAZ3 = COSD(AZ(3)) 00059 A1 = COSD(AZ(1)) - COSAZ3 00060 A2 = COSD(AZ(2)) - COSAZ3 00061 B1 = SIND(AZ(1)) - SINAZ3 00062 B2 = SIND(AZ(2)) - SINAZ3 00063 U=A1*B2-A2*B1 00064**2 DENOM=OFFZRO(U)*RADIUS*RADANS 00065**2 DX = (C1*B2 - C2*B1) / DENOM 00066 DY = (A1*C2 - A2*C1) / DENOM 00067 U=CLAT+DX 00068**2 ELOW=-89.999 00069**2 HIGH=89.999 00070**2 CLAT=TRIM(U,ELOW,HIGH) 00071**2 39 CLONG = CLONG - DY/COSD(CLAT) 00072 WRITE(LP,101) KNTMAX 00073 101 FORMAT ('0ERROR IN SUBROUTINE CENTER. MORE THAN',I4, 00074**2 ' ' ITERATIONS REQUIRED.') 00075**2 CLAT=0. 00076**2 CLONG=0. 00077**2 DISTMN = 1.E50 00078 RETURN 00079 C 00080 80 CALL GEOINV(ALAT(3),ALONG(3), ALAT(I),ALONG(I), DISTMN,AZ,AZ2) 00081 DISTMN = DISTMN/2. 00082 CALL GEODIR(ALAT(3),ALONG(3), DISTMN,AZ, CLAT,CLONG) 00083 RETURN 00084 END 00085 SUBROUTINE CSPHTR(LATI,LONGI, LAT,LONG) 00086 C*** .GIVEN. 3 GEO. PTS .. (ALAT(I), ALONG(I)), I=1,3) 00087 C .RETURNS. (LAT,LONG) =APPROX. CENTER OF 3-PT CIRCLE ON UNIT SPHERE00088 C .AS. PT WHERE NORMAL FROM ORIGIN TO 3-PT PLANE CUTS THE SPHERE 00089 IMPLICIT REAL*8(A-H, J-Z) 00090**2 REAL*8 LATI(3),LONGI(3), X(3),Y(3),Z(3) 00091**2 LOGICAL EQPT 00092 EQPT(X1,Y1, X2,Y2) = X1.EQ.X2 .AND. Y1.EQ.Y2 00093 AREA(X1,Y1, X2,Y2, X3,Y3) = 00094 " X2*Y3 - X3*Y2 + X3*Y1 - X1*Y3 + X1*Y2 - X2*Y1 00095 DET(X1,Y1,Z1, X2,Y2,Z2, X3,Y3,Z3) = 00096 = X1*Y2*Z3 - X1*Y3*Z2 + X2*Y3*Z1 - X2*Y1*Z3 + X3*Y1*Z2 - X3*Y2*Z1 00097 SIND(U)=DSIN(U*RADANS) 00098**2 COSD(U)=DCOS(U*RADANS) 00099**2 SQRT(U)=DSQRT(U) 00100**2 ATAN2(U,U1)=DATAN2(U,U1) 00101**2 AMOD(U1,U2)=DMOD(U1,U2) 00102**2 DATA RADANS/.017453292519943/ 00103 C 00104 DO 21 I = 1, 3 00105 Z(I) = SIND(LATI(I)) 00106 CLA = SQRT(1. - Z(I)**2) 00107 X(I) = CLA * SIND(LONGI(I)) 00108 21 Y(I) = CLA * COSD(LONGI(I)) 00109 C 00110 LAMBDA = AREA(Y(1),Z(1), Y(2),Z(2), Y(3),Z(3)) 00111 MU = AREA(Z(1),X(1), Z(2),X(2), Z(3),X(3)) 00112 NU = AREA(X(1),Y(1), X(2),Y(2), X(3),Y(3)) 00113 C.. TEST FOR 2 COINCIDENT PTS. 00114 I = 0 00115 IF( .NOT.(LAMBDA.EQ.0 .AND. MU.EQ.0 .AND. NU.EQ.0) ) GO TO 25 00116 I = 1 00117 IF( EQPT(LATI(I),LONGI(I), LATI(3),LONGI(3)) ) I = 2 00118 LAMBDA = (X(I) + X(3))/2. 00119 MU = (Y(I) + Y(3))/2. 00120 NU = (Z(I) + Z(3))/2. 00121 25 RR = LAMBDA**2 + MU**2 00122 LAT = ATAN2(NU, SQRT(RR)) / RADANS 00123 LONG = ATAN2(LAMBDA, MU) / RADANS 00124 C.. IF CENTER-PT AND CENTER-OF-SPHERE ON SAME SIDE OF PLANE THRU 3 PTS00125 C CHOODE ANTIPODAL CENTER. 00126 IF(I.NE.0 .OR. DET(X(1),Y(1),Z(1), X(2),Y(2),Z(2), X(3),Y(3),Z(3))00127 . .GT. RR + NU**2 ) GO TO 29 00128 LAT = -LAT 00129 LONG = LONG + 180. 00130 29 U1=LONG+360. 00131**2 U2=360. 00132**2 LONG=AMOD(U1,U2) 00133**2 C 00134 RETURN 00135 END 00136 SUBROUTINE CTRLIN(I1,I2,N1,N2) 00137**2 C*** COMPUTE MEDIAN-LINE 00138 IMPLICIT REAL*8(A-H,O-Z), INTEGER(I-N) 00139**2 COMMON/DATA/LWRK,I1Z(300),I2Z(300),IZZ(300), ISHORE(300),RAD(300) 00140 COMMON /MIDLIN/ J,LSHORC,CLAT(300),CLONG(300) 00141**2 C 00142 J = 0 00143 I1NEXT = N1 + 1 - I1 00144 I2NEXT = N2 + 1 - I2 00145 CALL MEDLIN(I1,I2) 00146 JOLD = J 00147 CALL FLIP 00148 CALL MEDLIN(I1NEXT,I2NEXT) 00149 C 00150 JOLD = JOLD + 1 00151 IF( JOLD.GT.J ) GO TO 64 00152 DO 62 I = JOLD, J 00153 I1Z(I) = N1 + 1 - I1Z(I) 00154 I2Z(I) = N2 + 1 - I2Z(I) 00155 IF( ISHORE(I).EQ.1 ) IZZ(I) = N1 + 1 - IZZ(I) 00156 IF( ISHORE(I).NE.1 ) IZZ(I) = N2 + 1 - IZZ(I) 00157 62 CONTINUE 00158 64 CALL FLIP 00159 C 00160 RETURN 00161**2 END 00162**2 SUBROUTINE DMS2D (DMS,DEG) 00163 C**THIS ROUTINE CONVERTS A (DEGREEE,MINUTE,SECOND) STRING TO (DEGREES). 00164 IMPLICIT REAL*8(A-H,O-Z), INTEGER(I-N) 00165**2 DIMENSION DMS(3) 00166 DATA C/1.E-10/ 00167 DEG = (DMS(3)/60.0+DMS(2))/60.0+DMS(1) 00168 RETURN 00169 ENTRY D2DMS(DMS,DEG) 00170**2 DMS(1)=DINT(DEG+C) 00171**2 DMS(2)=(DEG-DMS(1))*60.0+C 00172**2 I=DMS(2) 00173**2 Z=DMS(2) 00174**2 DMS(3) = (Z-DMS(2))*60.0 00175 RETURN 00176 END 00177 SUBROUTINE FLIP 00178 C*** FLIPS END-FOR-END .. CENTER- + SHORE-LINE IN /MEDLIN/, /SHORE1/ 2/00179 IMPLICIT REAL*8(A-H,O-Z), INTEGER(I-N) 00180**2 COMMON/SHORE1/ N1, LSHOR1,ALAT1(200),ALONG1(200) 00181 COMMON/SHORE2/ N2, LSHOR2,ALAT2(200),ALONG2(200) 00182 COMMON/MIDLIN/ J , LSHORC,CLAT (300),CLONG (300) 00183 COMMON/DATA/ I1Z(300),I2Z(300),IZZ(300), ISHORE(300),RAD(300) 00184 C 00185 C AND ADJUSTS PTRS IN /DATA/ 00186 JH = N1/2 00187 IF( JH.EQ.0 ) GO TO 15 00188 DO 10 K = 1,JH 00189 KC = N1+1-K 00190 T = ALAT1(KC) 00191 ALAT1(KC) = ALAT1(K) 00192 ALAT1(K) = T 00193 T = ALONG1(KC) 00194 ALONG1(KC) = ALONG1(K) 00195 10 ALONG1(K) = T 00196 15 JH = N2/2 00197 IF( JH.EQ.0 ) GO TO 25 00198 DO 20 K=1,JH 00199 KC = N2+1-K 00200 T = ALAT2(KC) 00201 ALAT2(KC) = ALAT2(K) 00202 ALAT2(K) = T 00203 T = ALONG2(KC) 00204 ALONG2(KC) = ALONG2(K) 00205 20 ALONG2(K) = T 00206 25 JH = J/2 00207 IF( JH.EQ.0 ) RETURN 00208 DO 30 K=1,JH 00209 KC = J+1-K 00210 IT = I1Z(KC) 00211 I1Z(KC) = I1Z(K) 00212 I1Z(K) = IT 00213 IT = I2Z(KC) 00214 I2Z(KC) = I2Z(K) 00215 I2Z(K) = IT 00216 IT = IZZ(KC) 00217 IZZ(KC) = IZZ(K) 00218 IZZ(K) = IT 00219 IT = ISHORE(KC) 00220 ISHORE(KC) = ISHORE(K) 00221 ISHORE(K) = IT 00222 T = RAD(KC) 00223 RAD(KC) = RAD(K) 00224 RAD(K) = T 00225 T = CLAT(KC) 00226 CLAT(KC) = CLAT(K) 00227 CLAT(K) = T 00228 T = CLONG(KC) 00229 CLONG(KC) = CLONG(K) 00230 30 CLONG(K) = T 00231 RETURN 00232 END 00233 SUBROUTINE GEODIR (LAT1,LONG1,DIST,AZ1,LAT2,LONG2,AZ2) 00234 C** THIS SUBROUTINE SOLVES THE DIRECT GEODESIC PROBLEM VIA 00235 C** THE CLARKE SOLUTION WITH BRUNAVS MODIFICATIONS. NOTE THAT 00236 C** THE PARAMETER AZ2 IS NOT COMPUTED. 00237 IMPLICIT REAL*8(A-H,O-Z), INTEGER(I-N) 00238**2 REAL*8 LAT1,LONG1,LAT2,LONG2 00239**2 COMMON/SPHDTA/ A,B, ESQ 00240 DATA RADANS/.017453292519/ 00241**2 SIN(U)=DSIN(U) 00242**2 COS(U)=DCOS(U) 00243**2 SQRT(U)=DSQRT(U) 00244**2 ATAN(U)=DATAN(U) 00245**2 ASIN(U)=DARSIN(U) 00246**2 C 00247 E1SQ = ESQ / (1. - ESQ) 00248 ALAT1 = LAT1*RADANS 00249 SINALF = SIN(AZ1*RADANS) 00250 COSALF = COS(AZ1*RADANS) 00251 SNLAT1 = SIN(ALAT1) 00252 CSLAT1 = COS(ALAT1) 00253 ANU = A / SQRT(1. - ESQ*SNLAT1**2) 00254 Q1 = DIST*COSALF/ANU 00255 Q2 = DIST*SINALF/ANU 00256 P = (E1SQ*CSLAT1**2)/6.0 00257 ONEPQ2 = 1.0+P*Q1**2 00258 THET = DIST*ONEPQ2/ANU 00259 SNTHET = SIN(THET) 00260 SINA = SINALF*ONEPQ2 00261 COSA = COSALF*(1.0-P*Q2**2) 00262 SINPHI = SNLAT1*COS(THET) + CSLAT1*SNTHET*COSA 00263 TNLAT2 = ((1.0+E1SQ)*SINPHI-E1SQ*SNLAT1)/(SQRT(1.0-SINPHI**2)) 00264 LAT2 = ATAN(TNLAT2)/RADANS 00265 SNDPHI = -SNTHET*SINA/SQRT(1.0-SINPHI**2) 00266 DPHI = ASIN(SNDPHI) 00267 LONG2 = LONG1+DPHI/RADANS 00268 AZ2 = 0.0 00269 RETURN 00270 END 00271 SUBROUTINE GEOINV (LAT1,LONG1,LAT2,LONG2,DIST,AZ,AZ2) 00272 C*** GIVEN LAT/LONG OF TWO POINTS ON CLARKE 1866 SPHEROID, RETURNS 00273 C*** GEODESIC DISTANCE, AZIMUTH, AND BACK AZIMUTH BY P.D.THOMAS METHOD. 00274 IMPLICIT REAL*8(A-H,O-Z), INTEGER(I-N) 00275**2 COMMON/SPHDTA/ A,B, ESQ 00276 REAL*8 LAT1,LONG1,LAT2,LONG2,KL,KK,L 00277**2 RCTLAT(T) = DATAN2(B0A*2.*T, 1. - T**2) 00278**2 ATAN2D(Y,X)=DATAN2(Y,X)/RADANS 00279**2 DATA RADANS/.017453292519/ 00280 TAN(U)=DTAN(U) 00281**2 SIN(U)=DSIN(U) 00282**2 COS(U)=DCOS(U) 00283**2 ACOS(U)=DARCOS(U) 00284**2 AMOD(U1,U2)=DMOD(U1,U2) 00285**2 ABS(U)=DABS(U) 00286**2 C 00287 DIST=0. 00288**2 AZ=0. 00289**2 AZ2=0. 00290**2 IF( LAT1.EQ.LAT2 .AND. LONG1.EQ.LONG2 ) RETURN 00291 B0A = B/A 00292 F = 1.0-B0A 00293 DLR = ABS(LONG2 - LONG1)*RADANS 00294 W=LAT1/2.*RADANS 00295**2 T1R=RCTLAT(TAN(W)) 00296**2 W=LAT2/2.*RADANS 00297**2 T2R=RCTLAT(TAN(W)) 00298**2 TM = (T1R+T2R)/2.0 00299 DTM = (T2R-T1R)/2.0 00300 IF( LONG2.LT.LONG1 ) DTM = -DTM 00301 STM = SIN(TM) 00302 CTM = COS(TM) 00303 SDTM = SIN(DTM) 00304 CDTM = COS(DTM) 00305 KL = STM*CDTM 00306 KK = SDTM*CTM 00307 SDLMR = SIN(DLR/2.0) 00308 L = SDTM**2 + SDLMR**2*(CDTM**2 - STM**2) 00309 CD = 1. - 2.*L 00310 DL = ACOS(CD) 00311 SD = SIN(DL) 00312 T = DL/SD 00313 U = KL**2/(1. - L) 00314 V = KK**2/L 00315 X = 2.*(U + V) 00316 Y = 2.*(U - V) 00317 D = (2.*T)**2 00318 E = -2.*CD 00319 A1 = -D*E 00320 FF64 = F*F/64.0 00321 DIST = A*SD*(T-(F/4.0)*(T*X-Y)+FF64*(X*(A1+(T-(A1+E)/2.0)*X) 00322 * +Y*(-2.0*D+E*Y)+D*X*Y)) 00323 TDLPM = TAN((DLR+(-((E*(4.0-X)+2.0*Y)*((F/2.0)*T+FF64* 00324 * (32.0*T+(A1-20.0*T)*X-2.0*(D+2.0)*Y))/4.0)* 00325 * TAN(DLR)))/2.0) 00326 W=CTM*TDLPM 00327**2 HAPBR=ATAN2D(SDTM,W) 00328**2 W=STM*TDLPM 00329**2 HAMBR=ATAN2D(CDTM,W) 00330**2 IF( LONG2.LT.LONG1 ) HAMBR = -HAMBR 00331 W=HAPBR-HAMBR+360. 00332**2 Z=360. 00333**2 AZ=AMOD(W,Z) 00334**2 W=HAPBR+HAMBR+360. 00335**2 AZ2=AMOD(W,Z) 00336**2 RETURN 00337 END 00338 SUBROUTINE MEDBGN(I1,I2) 00339 C*** FIND SHORTEST-CHORD JOINING SHORELINES 1 AND 2 00340 C (WILL START HERE, AND GENERATE MIDDLE-LINE IN 2 DIRECTIONS ..) 00341 C ( I.E. .GO FORWARD., FLIP AND CONCATENATE WITH .GO BACKWARD. ) 00342 C .RETURNS. START-PAIR OF POINTS .. I1"TH ON SHORELINE 1, I2"TH ON 200343 IMPLICIT REAL*8(A-H,O-Z), INTEGER(I-N) 00344**2 COMMON/IOFILS/ IN, LP 00345 DATA KMAX/5/ 00346 C 00347 WRITE(LP,149) 00348 149 FORMAT('1') 00349**2 C 00350 DISTMN = 1.E50 00351 I1OLD=1 00352**2 I1=1 00353**2 I2OLD=1 00354**2 I2=1 00355**2 DO 15 KOUNT = 1, KMAX 00356 CALL MEDMIN(2,N2, I1,I2, DISTMN) 00357 CALL MEDMIN(1,N1, I2,I1, DISTMN) 00358 IF( I1.EQ.I1OLD .AND. I2.EQ.I2OLD ) RETURN 00359 I1OLD = I1 00360 15 I2OLD = I2 00361 C 00362 20 WRITE(LP,120) I1,I2 00363 120 FORMAT('1CANNOT FIND START-POINTS. STOPPING AT (I1,I2)=(',I4,',',00364**2 +I4,')') 00365**2 STOP 00366**2 C 00367 END 00368 SUBROUTINE MEDCPR 00369 C*** PRINT COMPUTED MEDIAN-LINE 00370 IMPLICIT REAL*8(A-H,O-Z), INTEGER(I-N) 00371**2 COMMON/IOFILS/ IN, LP 00372 COMMON/MIDLIN/ J , LSHORC,CLAT (300),CLONG (300) 00373 COMMON/DATA/LWRK,I1Z(300),I2Z(300),IZZ(300), ISHORE(300),RAD(300) 00374 REAL*8 ADMS(3),BDMS(3) 00375**2 C 00376 WRITE(LP,140) 00377 140 FORMAT('1POINT'5X'LATITUDE'9X'LONGITUDE'9X'POINTS USED'12X 00378**2 *'RADIUS'/' NO. DEG MIN SEC DEG MIN SEC'6X'S1 S2', 00379**2 *9X'KILOMETRES'3X'NAUT. MILES'/) 00380**2 DO 50 I = 1, J 00381 CALL D2DMS(ADMS, CLAT (I)) 00382 CALL D2DMS(BDMS, CLONG(I)) 00383 RAD2=RAD(I)*5.3995677E-4 00384**2 WRITE(LP,151) I,ADMS,BDMS, I1Z(I),I2Z(I),IZZ(I), ISHORE(I), RAD(I)00385 *,RAD2 00386**2 151 FORMAT (I4,2(F7.0,F5.0,F6.2),I7,2I4,'-S',I1,-3PF12.3,2X,0PF9.2) 00387**2 50 CONTINUE 00388 C 00389 RETURN 00390 END 00391 C PROGRAM MEDIAN(INPUT,OUTPUT,TAPE8,TAPE6=INPUT,TAPE5=OUTPUT) 00392**2 C*** 00393 C*** .INPUTS. TURNING PTS ON 2 SHORELINES 00394 C.. (AND, OPTIONALLY, MODIFICATIONS .TO. INTERPOLATE MORE TURNING PTS)00395 C.. .COMPUTES. A MIDDLE-LINE, EQUIDISTANT FROM THE 2 SHORELINES 00396 C .VIA. BRUNAVS" SEARCH ALGORITHM, MIMICKING THE IHB MANUAL METHOD00397 C.. .LISTS. (GIVEN, MODIFIED) SHORELINES, AND COMPUTED MEDIAN-LINE 00398 C.. (ALSO .SAVES. THESE 3 LINES AS PT-SETS ON FILE "SAVFIL") 00399 C.. ( FOR LATER PROCESSING / PLOTTING ) 00400 C.. .NOTE. ALL LINES ARE INPUT, COMPUTED, AND LISTED AS GEOGRAPHIC 00401 C 00402 IMPLICIT REAL*8(A-H,O-Z), INTEGER(I-N) 00403**2 C*** COMMON/IOFILS/ STORES LOGICAL UNIT-NUMBERS .FOR. INPUT, OUTPUT 00404 COMMON/IOFILS/ IN, LP 00405 C*** COMMON/SHORE1/ STORES SHORELINE 1 (N1 = ACTUAL NO. OF PTS) 00406 COMMON/SHORE1/ N1, LSHOR1,ALAT1(200),ALONG1(200) 00407 C*** COMMON/SHORE2/ STORES SHORELINE 2 (N2 = ACTUAL NO. OF PTS) 00408 COMMON/SHORE2/ N2, LSHOR2,ALAT2(200),ALONG2(200) 00409 C*** COMMON/MIDLIN/ STORES MEDIAN-LINE (J = ACTUAL NO. OF PTS) 00410 COMMON/MIDLIN/ J , LSHORC,CLAT (300),CLONG (300) 00411 C*** COMMON/SHOREM/ WORKSPACE FOR MODIFYING SHORELINES 00412 COMMON/SHOREM/NMOD, LSHORM,MSHORE(100),MPT(100), EMRAD(100) 00413 C*** COMMON/DATA / WORKSPACE FOR COMPUTING MEDIAN-LINE 00414 COMMON/DATA/LWRK,I1Z(300),I2Z(300),IZZ(300), ISHORE(300),RAD(300) 00415 COMMON/SPHDTA/ A,B, ESQ 00416 INTEGER SAVFIL 00417 DATA SAVFIL/8/ 00418 C*** .INPUT. SHORELINES 1 AND 2 00419 CALL SHORIN(ALAT1,ALONG1, 1,N1) 00420 CALL SHORIN(ALAT2,ALONG2, 2,N2) 00421 C** MAKE MODIFICATIONS TO INPUT DATA. 00422 CALL SHORIM 00423 C** FIND STARTING POINTS. 00424 CALL MEDBGN(I1,I2) 00425 C*** COMPUTE MEDIAN-LINE 00426 CALL CTRLIN(I1,I2,N1,N2) 00427**2 C*** LIST MEDIAN-LINE 00428 CALL MEDCPR 00429 C*** SAVE SHORELINES (1, 2, AND MID-) 00430 CALL SHOREX(SAVFIL) 00431 C 00432 STOP 00433**2 END 00434 BLOCK DATA 00435**2 IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N) 00436**2 COMMON /IOFILS/IN,LP 00437**2 COMMON/SHORE1/N1,LSHOR1,ALAT1(200),ALONG1(200) 00438**2 COMMON /SHORE2/N2,LSHOR2,ALAT2(200),ALONG2(200) 00439**2 COMMON /MIDLIN/J,LSHORC,CLAT(300),CLONG(300) 00440**2 COMMON /SHOREM/NMOD,LSHORM,MSHORE(100),MPT(100),EMRAD(100) 00441**2 COMMON /DATA/LWRK,I1Z(300),I2Z(300),IZZ(300),ISHORE(300),RAD(300) 00442**2 COMMON /SPHDTA/A,B,ESQ 00443**2 DATA IN,LP/5,6/ 00444**2 DATA LSHOR1,LSHOR2,LSHORC/200,200,300/ 00445**2 DATA LSHORM,LWRK/100,300/ 00446**2 DATA A,B/6378206.4,6356583.8/,ESQ/.006768657997291/ 00447**2 END 00448**2 SUBROUTINE MEDLIN (II1,II2) 00449 C*** COMPUTE NEXT PT OF MEDIAN-LINE 00450 C** START-PTS .VIA. (II1, II2) 00451 IMPLICIT REAL*8(A-H,O-Z), INTEGER(I-N) 00452**2 COMMON/SHORE1/ N1, LSHOR1,ALAT1(200),ALONG1(200) 00453 COMMON/SHORE2/ N2, LSHOR2,ALAT2(200),ALONG2(200) 00454 COMMON/MIDLIN/ J , LSHORC,CLAT (300),CLONG (300) 00455 COMMON/DATA/LWRK,I1Z(300),I2Z(300),IZZ(300), ISHORE(300),RAD(300) 00456 COMMON/MEDTRI/ A(3),B(3) 00457 LOGICAL ZGOOD, ZEQPT,EQPT 00458 EQPT(X1,Y1, X2,Y2) = X1.EQ.X2 .AND. Y1.EQ.Y2 00459 C 00460 I1 = II1 00461 I2 = II2 00462 ZEQPT = EQPT(ALAT1(I1),ALONG1(I1), ALAT2(I2),ALONG2(I2)) 00463 99 IF (I1.GE.N1.AND.I2.GE.N2) RETURN 00464 A(1) = ALAT1(I1) 00465 B(1) = ALONG1(I1) 00466 A(2) = ALAT2(I2) 00467 B(2) = ALONG2(I2) 00468 CALL MEDLPT(ALAT1,ALONG1, N1,I1,I1NEW, I1,I2, DIST1,CLA1,CLO1) 00469 CALL MEDLPT(ALAT2,ALONG2, N2,I2,I2NEW, I1,I2, DIST2,CLA2,CLO2) 00470 50 IF (DIST1.EQ.1.E50.AND.DIST2.EQ.1.E50) RETURN 00471 CALL MEDLPR(DIST1,I1,I1NEW, DIST2,I2,I2NEW) 00472 IF (DIST1.GT.DIST2) GO TO 60 00473 C.. USE SHORE 1 DATA ... 00474 IF( ZEQPT .AND. I1.NE.II1 .AND. I2.EQ.II2 ) GO TO 61 00475 J = J+1 00476 I1Z(J) = I1 00477 I2Z(J) = I2 00478 IZZ(J) = I1NEW 00479 ISHORE(J) = 1 00480 RAD(J) = DIST1 00481 CLAT(J) = CLA1 00482 CLONG(J) = CLO1 00483 I1 = I1NEW 00484 GO TO 99 00485 C.. USE SHORE 2 DATA ... 00486 60 IF( ZEQPT .AND. I1.EQ.II1 .AND. I2.NE.II2 ) GO TO 62 00487 J = J + 1 00488 I1Z(J) = I1 00489 I2Z(J) = I2 00490 IZZ(J) = I2NEW 00491 ISHORE(J) = 2 00492 RAD(J) = DIST2 00493 CLAT(J) = CLA2 00494 CLONG(J) = CLO2 00495 I2 = I2NEW 00496 GO TO 99 00497 61 I1 = I1 + 1 00498 GO TO 99 00499 62 I2 = I2 + 1 00500 GO TO 99 00501 END 00502 SUBROUTINE MEDLPR(DIST1,I1,I1NEW, DIST2,I2,I2NEW) 00503 C*** PRINT CLOSEST-PT TRIPLES ON SHORELINES, AND ASSOCIATED DISTANCE 00504 IMPLICIT REAL*8(A-H,O-Z), INTEGER(I-N) 00505**2 COMMON/IOFILS/ IN, LP 00506 C 00507 DTEXT=1E50 00508**2 IF(DIST1.NE.1E50)DTEXT=DIST1 00509**2 WRITE(LP,151) DTEXT, I1,I1NEW, I2 00510 DTEXT=1E50 00511**2 IF(DIST2.NE.1E50)DTEXT=DIST2 00512**2 WRITE(LP,152) DTEXT, I1, I2,I2NEW 00513 151 FORMAT(' DISTANCE= ',E10.5,'KM. FROM SHORE 1, PT',I4,' TO SHORE 2,00514**2 + PTS',I4,',',I3) 00515**2 152 FORMAT('+',T70,' =',E10.5,' KM FROM SHORE 1 PT',I4,' TO SHORE 2, 00516**2 +PTS',I4,',',I3) 00517**2 C 00518 RETURN 00519 END 00520 SUBROUTINE MEDLPT(LAT,LONG, N,I,INEW, I1,I2, DIST,CLA,CLO) 00521 C*** 00522 C*** .GIVEN. I1"TH PT ON SHORELINE 1, I2"TH ON SHORELINE 2, 00523 C FIND ASSOCIATED NEAREST PT (INEW) ON GIVEN SHORELINE 00524 C .RETURNS. BEST CENTER (CLA,CLO) AND ASSOCIATED DISTANCE, DIST 00525 C ( DIST = INFINITE .IMPLIES. FAILURE ) 00526 IMPLICIT REAL*8(A-H,O-Z), INTEGER(I-N) 00527**2 COMMON/MEDTRI/ A(3),B(3) 00528 REAL*8 LAT(1),LONG(1) 00529**2 LOGICAL ZGOOD 00530 C 00531 DIST = 1.E50 00532 IF( I.GE.N ) RETURN 00533 C 00534 10 ISTART = I + 1 00535 DO 20 INEW = ISTART, N 00536 A(3) = LAT (INEW) 00537 B(3) = LONG(INEW) 00538 CALL CENTER(A,B, CLA,CLO, DIST) 00539 IF( DIST.EQ.1.E50 ) GO TO 20 00540 CALL MEDTST(CLA,CLO, I1,I2, ZGOOD) 00541 IF( ZGOOD ) RETURN 00542 20 CONTINUE 00543 C 00544 RETURN 00545 END 00546 SUBROUTINE MEDMIN(K,N, II,IMIN,DISTMN) 00547 C*** .GIVEN. II"TH PT OF K"TH SHORELINE, 00548 C SEARCH FOR NEAREST PT ON OTHER SHORELINE. 00549 C .RETURNS. IMIN"TH PT, AND ASSOCIATED DISTANCE DISTMN 00550 IMPLICIT REAL*8(A-H,O-Z), INTEGER(I-N) 00551**2 COMMON/SHORE1/ N1, LSHOR1,ALAT1(200),ALONG1(200) 00552 COMMON/SHORE2/ N2, LSHOR2,ALAT2(200),ALONG2(200) 00553 C 00554 DO 15 I = 1, N 00555 IF( K.EQ.1 ) DIST = S(ALAT1(II),ALONG1(II), ALAT2(I),ALONG2(I)) 00556 IF( K.NE.1 ) DIST = S(ALAT2(II),ALONG2(II), ALAT1(I),ALONG1(I)) 00557 IF( DIST.GE.DISTMN ) GO TO 15 00558 DISTMN = DIST 00559 IMIN = I 00560 15 CONTINUE 00561 C 00562 RETURN 00563**2 END 00564**2 SUBROUTINE MEDTST(CLAT,CLONG, I1,I2, ZGOOD) 00565 C*** TEST FOR SATISFACTORY CENTER-POINT 00566 IMPLICIT REAL*8(A-H,O-Z), INTEGER(I-N) 00567**2 LOGICAL ZGOOD 00568 COMMON/SHORE1/ N1, LSHOR1,ALAT1(200),ALONG1(200) 00569 COMMON/SHORE2/ N2, LSHOR2,ALAT2(200),ALONG2(200) 00570 DATA TOL/1.0/ 00571 C 00572 DIST = S (CLAT,CLONG,ALAT1(I1),ALONG1(I1)) - TOL 00573 DO 10 I=1,N1 00574 ZGOOD = S(CLAT,CLONG, ALAT1(I),ALONG1(I)).GE.DIST 00575 IF( .NOT.ZGOOD ) RETURN 00576 10 CONTINUE 00577 DO 20 I=1,N2 00578 ZGOOD = S(CLAT,CLONG, ALAT2(I),ALONG2(I)).GE.DIST 00579 IF( .NOT.ZGOOD ) RETURN 00580 20 CONTINUE 00581 C 00582 RETURN 00583 END 00584 FUNCTION S( LAT1,LONG1, LAT2,LONG2 ) 00585 C** ANDOYER-LAMBERT FORMULA USING PARAMETRIC LATS, ANGLES IN DEGREES 00586 C** ACCURACY IS 1 M IN 800 KM AND 7 M IN 10,000 KM. S.ORAAS, JAN / 73. 00587 IMPLICIT REAL*8(A-H,O-Z), INTEGER(I-N) 00588**2 REAL*8 LAT1,LAT2,LONG1,LONG2 00589**2 COMMON/SPHDTA/ A,B, ESQ 00590 LOGICAL EQPT 00591 EQPT(X1,Y1, X2,Y2) = DABS(X2 - X1) + DABS(Y2 - Y1) .LE. 1E-6 00592**2 DATA RADANS/.017453292519943/ 00593 SIN(U)=DSIN(U) 00594**2 COS(U)=DCOS(U) 00595**2 SQRT(U)=DSQRT(U) 00596**2 ATAN2(U1,U2)=DATAN2(U1,U2) 00597**2 C 00598 S = 0. 00599 IF( EQPT(LAT1,LONG1, LAT2,LONG2) ) RETURN 00600 C 00601 BA = B/A 00602 C = 1.0 - BA*BA 00603 F4 = (1.0-BA)/4.0 00604 COSDL = COS((LONG1 - LONG2) * RADANS ) 00605 SINPHI = SIN(LAT1*RADANS) 00606 SIN1 = BA*SINPHI/SQRT(1.-C*SINPHI*SINPHI) 00607 SINPHI = SIN(LAT2*RADANS) 00608 SIN2 = BA*SINPHI/SQRT(1.-C*SINPHI*SINPHI) 00609 COSD = SIN1*SIN2 + COSDL*SQRT((1.-SIN1*SIN1)*(1.-SIN2*SIN2)) 00610 SIND = SQRT(1.- COSD*COSD) 00611 D = ATAN2(SIND,COSD) 00612 TQ = SIN1 - SIN2 00613 TP = SIN1 + SIN2 00614 S = A*(D-F4*(TP*TP*(D-SIND)/(1.+COSD)+TQ*TQ*(D+SIND)/(1.-COSD))) 00615 RETURN 00616 END 00617 SUBROUTINE SHOREX( SAVFIL ) 00618 C*** SAVE SHORELINES (1, 2, AND MID-) ON FILE "SAVFIL" 00619 IMPLICIT REAL*8(A-H,O-Z), INTEGER(I-N) 00620**2 INTEGER SAVFIL 00621 COMMON/SHORE1/ N1, LSHOR1,ALAT1(200),ALONG1(200) 00622 COMMON/SHORE2/ N2, LSHOR2,ALAT2(200),ALONG2(200) 00623 COMMON/MIDLIN/ J , LSHORC,CLAT (300),CLONG (300) 00624 C 00625 CALL SHORSV(ALAT1,ALONG1, N1,1, SAVFIL) 00626 CALL SHORSV(ALAT2,ALONG2, N2,2, SAVFIL) 00627 CALL SHORSV(CLAT ,CLONG , J ,0, SAVFIL) 00628 C 00629 ENDFILE SAVFIL 00630 REWIND SAVFIL 00631 C 00632 RETURN 00633 END 00634 SUBROUTINE SHORIM 00635 C*** IF REQ"D, INPUT MODS, AND MODIFY SHORELINES 1 AND 2 00636 IMPLICIT REAL*8(A-H,O-Z), INTEGER(I-N) 00637**2 COMMON/IOFILS/ IN, LP 00638 COMMON/SHOREM/NMOD, LSHORM,MSHORE(100),MPT(100), EMRAD(100) 00639 COMMON/SHORE1/ N1, LSHOR1,ALAT1(200),ALONG1(200) 00640 COMMON/SHORE2/ N2, LSHOR2,ALAT2(200),ALONG2(200) 00641 INTEGER TEXT(30),BLANKS,BUFF(30) 00642**2 DATA BLANKS/' '/ 00643**2 C 00644 SQRT(U)=DSQRT(U) 00645**2 10 READ(IN,100,END=24)(TEXT(I),I=1,30) 00646**2 100 FORMAT(30A1) 00647**2 IF(TEXT(1).NE.BLANKS) GO TO 10 00648**2 CALL CORE(BUFF,30) 00649**2 WRITE(10,100)(TEXT(I),I=1,30) 00650**2 CALL CORE(BUFF,30) 00651**2 READ(10,114)NMOD,SCALE,ERRMAX 00652**2 C 00653 114 FORMAT(I10,2F10.0) 00654 IF( NMOD.LE.0 ) GO TO 24 00655 WRITE(LP,118) SCALE, ERRMAX 00656 118 FORMAT('1MODIFICATIONS REQUIRED. SCALE =',F10.0,' MAX ERROR= ', 00657**2 +F5.3,' INCHES.',/,'0NO SHORE POINTS RADIUS'/) 00658**2 DO 20 I = 1, NMOD 00659 READ (IN,113) MSHORE(I),MPT(I), EMRAD(I) 00660 113 FORMAT(I1,I5,3X,F10.0) 00661 ITEMP = MPT(I)+1 00662 20 WRITE(LP,112) I,MSHORE(I),MPT(I),ITEMP,EMRAD(I) 00663 112 FORMAT(I3,I6,I7,',',I3,F12.0) 00664**2 FACTOR = ERRMAX*8.0*.0254*SCALE 00665 DO 21 IZ = 1, NMOD 00666 I = NMOD+1-IZ 00667 STEP = SQRT (FACTOR*EMRAD(I)) 00668 II = MPT(I)+1 00669 IF( MSHORE(I).EQ.1 ) CALL SHORMD(ALAT1,ALONG1,1, N1,II, DIST,STEP)00670 IF( MSHORE(I).EQ.2 ) CALL SHORMD(ALAT2,ALONG2,2, N2,II, DIST,STEP)00671 21 CONTINUE 00672 C 00673 CALL SHORMP(ALAT1,ALONG1,1, N1) 00674 CALL SHORMP(ALAT2,ALONG2,2, N2) 00675 RETURN 00676 C 00677 24 WRITE(LP,117) 00678 117 FORMAT ('1NO MODIFICATIONS TO SHORELINE DATA REQUIRED.') 00679**2 C 00680 RETURN 00681 END 00682 SUBROUTINE SHORIN(LAT,LONG,K, N) 00683 C*** INPUT SHORELINES 1 AND 2 00684 IMPLICIT REAL*8(A-H,O-Z), INTEGER(I-N) 00685**2 COMMON/IOFILS/ IN, LP 00686 REAL*8 LAT(1),LONG(1), ADMS(3),BDMS(3) 00687**2 INTEGER TEXT(80),BLANKS,BUFF(80) 00688**2 DATA BLANKS/' '/ 00689**2 C 00690 N = 0 00691 GO TO 11 00692 C 00693 10 WRITE(LP,101) TEXT 00694 C... LIST, AND IGNORE, ANY COMMENT-CARDS 00695 11 READ(IN,100,END=99)(TEXT(I),I=1,80) 00696**2 C... IF NULL-FILE, .RETURNS. NO POINTS 00697 100 FORMAT(80A1) 00698**2 101 FORMAT(80A1) 00699**2 IF(TEXT(1).NE.BLANKS) GO TO 10 00700**2 CALL CORE(BUFF,80) 00701**2 WRITE(10,100)(TEXT(I),I=1,80) 00702**2 CALL CORE(BUFF,80) 00703**2 READ(10,110)N 00704**2 110 FORMAT(I5) 00705 WRITE(LP,120) K, N 00706 120 FORMAT('1SHORELINE NO.'I2,5X'NO. POINTS ='I5/) 00707**2 IF( N.LE.0 ) RETURN 00708 C 00709 DO 15 I = 1, N 00710 READ (IN,111) ADMS,BDMS 00711 111 FORMAT(2(F3.0,F4.1,F3.0)) 00712 WRITE(LP,121) I, ADMS,BDMS 00713 121 FORMAT(I6,2(F10.0,2F7.1,2X)) 00714 CALL DMS2D(ADMS, LAT (I)) 00715 15 CALL DMS2D(BDMS, LONG(I)) 00716 C 00717 99 RETURN 00718**2 END 00719 SUBROUTINE SHORMD(LAT,LONG,K, M,II, DIST,STEP) 00720 C** MAKE MODIFICATIONS TO K"TH SHORELINE 00721 IMPLICIT REAL*8(A-H,O-Z), INTEGER(I-N) 00722**2 COMMON/IOFILS/ IN, LP 00723 REAL*8 LAT(1),LONG(1), ADMS(3),BDMS(3) 00724**2 C 00725 IP = II - 1 00726 CALL GEOINV(LAT(IP),LONG(IP), LAT(II),LONG(II), DIST,AZ1,AZ2) 00727 C 00728 N = DIST/STEP 00729 IF( N.LE.0 ) RETURN 00730 STEP = DIST/(N+1) 00731 C 00732 DO 13 ITEMP = II, M 00733 IZ = M + II - ITEMP 00734 LAT (IZ+N) = LAT (IZ) 00735 13 LONG(IZ+N) = LONG(IZ) 00736 DO 14 ITEMP = 1, M 00737 IZ = IP + ITEMP 00738 14 CALL GEODIR(LAT(IP),LONG(IP),ITEMP*STEP,AZ1, LAT(IZ),LONG(IZ),AZ2)00739 M = M + N 00740 RETURN 00741 C 00742 ENTRY SHORMP(LAT,LONG,K,M) 00743**2 C*** PRINT MODIFIED SHORELINE ... 00744 C 00745 WRITE(LP,120) K, M 00746 120 FORMAT('1MODIFIED SHORELINE NO.'I2,5X'NO. POINTS ='I5/) 00747**2 DO 21 I = 1, M 00748 CALL D2DMS(ADMS, LAT (I)) 00749 CALL D2DMS(BDMS, LONG(I)) 00750 21 WRITE(LP,121) I, ADMS,BDMS 00751 121 FORMAT(I6,2(F10.0,2F7.1,2X)) 00752 C 00753 RETURN 00754 END 00755 SUBROUTINE SHORSV(LAT,LONG, N,KSEG, SAVFIL) 00756 C*** SAVE (ON SAVFIL) K"TH SHORELINE (0TH = MEDIAN-LINE) 00757 IMPLICIT REAL*8(A-H,O-Z), INTEGER(I-N) 00758**2 REAL*8 LAT(1),LONG(1) 00759**2 INTEGER SAVFIL 00760 C 00761 IF( N.LE.0 ) RETURN 00762 DO 15 I = 1, N 00763 15 WRITE(SAVFIL) LAT(I),LONG(I), KSEG 00764 C 00765 RETURN 00766 END 00767