PROGRAM Bcs_Tool ! Fixed-format Fortran 90 source code (by Peter Bird, UCLA, July 1998) ! was converted to free-format Fortran 90 source code in January 2016. CHARACTER*80 :: infile, outfil, line, line2 INTEGER :: bctype REAL :: lat, lon ! INTRODUCTION WRITE (*, 1) 1 FORMAT ( & &' -= BCs_tool =-' & &/' A utility program for use with SHELLS, to help create the' & &/' Boundary ConditionS (.BCS) files required by SHELLS.' & &/' Wherever the .BCS file contains a boundary condition type 3,' & &/' this program computes an absolute velocity from an absolute' & &/' rotation pole and rate entered at the keyboard, and revises' & &/' the type-3 boundary condition to a type-2 boundary condition ' & &/' with the velocity and azimuth filled in.' & &/' (Note: SHELLS includes a similar function, but only provides' & &/' one model (NUVEL-1A) and converts the bc to type-1.)' & &//' Press [Enter] to continue...'\) READ (*, "(A)", iostat = ios) line ! INPUT AND OUTPUT FILE NAMES 2 WRITE (*, 3) 3 FORMAT ( & &//' REQUIRED INPUT FILE:' & &/' -------------------------------------------------------' & &,'------------' & &/' Title line, containing description of boundary conditions?' & &/' 1 1024 3 34.957 ' & &,' -102.455' & &/' 2 1079 0 35.123 ' & &,' -102.455' & &/' 3 1134 -1 35.765 ' & &,' -102.901' & &/' 4 1250 1 0.0 90.0 36.773 ' & &,' -103.111' & &/' 5 1298 2 3.12E-9 0.0 37.223 ' & &,' -105.234' & &/' .........................................................' & &,'..........' & &/' ---------------------------------------------------------' & &,'----------' & &/' (ordinal) (node) (bc-type) (velocity) (azimuth) (latitute' & &,') (longitude)' & &//' The easiest way to get this file is to run SHELLS with n' & &,'o .BCS file,' & &/' capture the output, and edit the list of boundary nodes' & &/' (with latitudes and longitudes) that it supplies.' & &/' When editing, insert bc-type 0, 1, 2, or 3 for each node.' & &/' This program only modifies type-3 lines, converting to type-2.' & &//' Enter [drive][\path\]filename of rough .BCS (input) file:') READ (*, '(A)') infile OPEN (UNIT = 1, FILE = infile, STATUS = 'OLD', iostat = ios, & & pad = 'YES') IF (iostat /= 0) THEN WRITE (*, "(/' ERROR: File not found (in this directory)!')") GO TO 2 END IF 4 WRITE (*, 5) 5 FORMAT ( & &//' Enter [drive][\path\]filename for NEW (output) .BCS file: ') READ (*, "(A)") outfil OPEN (UNIT = 2, FILE = outfil, STATUS = 'NEW', iostat = ios) IF (iostat /= 0) THEN WRITE (*, "(/' ERROR: Supply a NEW filename.')") GO TO 4 END IF ! GET POLE POSITION AND RATE WRITE (*, 20) 20 FORMAT (/' Enter radius of planet, in kilometers: '\) READ (*, * ) radius ! INTERNALLY, COMPUTATIONS ARE DONE IN SI radius = radius * 1.E3 25 WRITE (*, 30) 30 FORMAT (/' Select units option for rotation rate:' & & /' 1 = degrees per million years;' & & /' 2 = radians per second.' & & /' Enter integer (1-2): '\) READ (*, * ) itype IF ((itype < 1).OR.(itype > 2)) GO TO 25 IF (itype == 1) THEN factor = (3.14159 / 180.) / (1.E6 * 3.15576E7) elseif (itype == 2) THEN factor = 1. WRITE (*, 32) 32 FORMAT (/' NOTE: -16' & & /' To enter very small numbers, like 5.362x10 ' & & ,' radians/second' & & /' on one line, use this format: 5.362E-16'/) END IF ! ENTER POLE ! ("FACTOR" WILL CONVERT ROTATION RATES TO RADIANS/SECOND = SI) 42 WRITE (*, 43) 43 FORMAT (/' Absolute rotation pole LATITUDE in ' & & ,'degrees (N = +): '\) READ (*, * ) pollat IF ((pollat > 90.).OR.(pollat < -90.)) GO TO 42 pollat = pollat * (3.14159 / 180.) WRITE (*, 45) 45 FORMAT (' Absolute rotation pole LONGITUDE in ' & & ,'degrees (E = +): '\) READ (*, * ) pollon pollon = pollon * (3.14159 / 180.) WRITE (*, 47) 47 FORMAT (' Absolute rotation', & & ' RATE (counterclockwise = +): '\) READ (*, * ) rotrat rotrat = rotrat * factor sumx = rotrat * COS(pollat) * COS(pollon) sumy = rotrat * COS(pollat) * SIN(pollon) sumz = rotrat * SIN(pollat) rotrat = SQRT(sumx**2 + sumy**2 + sumz**2) upx = sumx / rotrat upy = sumy / rotrat upz = sumz / rotrat rotrat = rotrat / factor pollon = ATAN2(upy, upx) equate = SQRT(upx**2 + upy**2) pollat = ATAN2(upz, equate) pollon = pollon * (180. / 3.14159) pollat = pollat * (180. / 3.14159) pollon = AMOD(pollon + 720., 360.) polln2 = 360. - pollon WRITE (*, 60) pollat, pollon, polln2, rotrat 60 FORMAT (' ---------------------------------------' & & /' Rotation pole is located at ',F7.3,' degrees North' & & /' and ',F7.3,' degrees east (',F7.3,' degrees west)' & & /' with rotation rate of ',1P,E10.3) IF (itype == 1) THEN WRITE (*, 61) 61 FORMAT (' in units of degrees per million years.') elseif (itype == 2) THEN WRITE (*, 62) 62 FORMAT (' in units of radians per second.') END IF ! SCAN FILE, COMPUTING VELOCITIES FOR TYPE-3 NODES READ (1, "(A)") line WRITE (*, "(/' Title line from input .BCS file follows:')") WRITE (*, "(' ',A)") trim(line) WRITE (*, "(/' Enter new title line, or just [Enter] to copy:')") READ (*, "(A)") line2 IF (LEN_trim(line2) == 0) THEN WRITE (2, "(A)") trim(line) ELSE WRITE (2, "(A)") trim(line2) END IF n = 1 DO READ (1, "(A)", iostat = ios) line n = n + 1 IF (ios < 0) exit READ (line, * ) i, node, bctype IF (bctype == - 1) THEN READ (line, * , iostat = ios) i, node, bctype, & & vellat, vellon IF (ios /= 0) THEN WRITE (*, "(/' ERROR in input file, line ',I6)") n WRITE (*, "(' Did not find required data:')") WRITE (*, "(' I NODE BC-TYPE=-1 ', & & 'LATITUDE LONGITUDE')") STOP ' ' END IF WRITE (2, 89) i, node, vellat, vellon 89 FORMAT (I4,2X,I6,' -1 ',19X, & & 2X,F8.4,2X,F9.4) ELSE IF (bctype == 0) THEN READ (line, * , iostat = ios) i, node, bctype, & & vellat, vellon IF (ios /= 0) THEN WRITE (*, "(/' ERROR in input file, line ',I6)") n WRITE (*, "(' Did not find required data:')") WRITE (*, "(' I NODE BC-TYPE=0 ', & & 'LATITUDE LONGITUDE')") STOP ' ' END IF WRITE (2, 90) i, node, vellat, vellon 90 FORMAT (I4,2X,I6,' 0 ',19X, & & 2X,F8.4,2X,F9.4) ELSE IF (bctype == 1) THEN READ (line, * , iostat = ios) i, node, bctype, vel, azim, & & vellat, vellon IF (ios /= 0) THEN WRITE (*, "(/' ERROR in input file, line ',I6)") n WRITE (*, "(' Did not find required data:')") WRITE (*, "(' I NODE BC-TYPE=1 VELOCITY AZIMUTH ', & & 'LATITUDE LONGITUDE')") STOP ' ' END IF WRITE (2, 91) i, node, vel, azim, vellat, vellon 91 FORMAT (I4,2X,I6,' 1 ',1P,E10.3,1X,0P,F8.3, & & 2X,F8.4,2X,F9.4) ELSE IF (bctype == 2) THEN READ (line, * , iostat = ios) i, node, bctype, vel, azim, & & vellat, vellon IF (ios /= 0) THEN WRITE (*, "(/' ERROR in input file, line ',I6)") n WRITE (*, "(' Did not find required data:')") WRITE (*, "(' I NODE BC-TYPE=2 VELOCITY AZIMUTH ', & & 'LATITUDE LONGITUDE')") STOP ' ' END IF WRITE (2, 92) i, node, vel, azim, vellat, vellon 92 FORMAT (I4,2X,I6,' 2 ',1P,E10.3,1X,0P,F8.3, & & 2X,F8.4,2X,F9.4) ELSE IF (bctype == 3) THEN READ (line, * , iostat = ios) i, node, bctype, vellat, vellon IF (ios /= 0) THEN WRITE (*, "(/' ERROR in input file, line ',I6)") n WRITE (*, "(' Did not find required data:')") WRITE (*, "(' I NODE BC-TYPE=3 LATITUDE LONGITUDE')") STOP ' ' END IF lat = vellat * (3.14159 / 180.) lon = vellon * (3.14159 / 180.) rx = radius * COS(lat) * COS(lon) ry = radius * COS(lat) * SIN(lon) rz = radius * SIN(lat) vx = sumy * rz - sumz * ry vy = sumz * rx - sumx * rz vz = sumx * ry - sumy * rx vel = SQRT(vx**2 + vy**2 + vz**2) ! FIND UNIT VECTORS POINTING NORTH AND EAST TO EXTRACT COMPONENTS unx = -SIN(lat) * COS(lon) uny = -SIN(lat) * SIN(lon) unz = COS(lat) uex = -SIN(lon) uey = COS(lon) uez = 0. vn = vx * unx + vy * uny + vz * unz ve = vx * uex + vy * uey + vz * uez azim = 90. - (180. / 3.14159) * ATAN2(vn, ve) IF (azim < 0.0) azim = azim + 360. WRITE (2, 93) i, node, vel, azim, vellat, vellon 93 FORMAT (I4,2X,I6,' 2 ',1P,E10.3,1X,0P,F8.3, & & 2X,F8.4,2X,F9.4) ELSE WRITE (*, 99) bctype 99 FORMAT (/' Unexpected BC-type ',I2,' encountered. ???') STOP ' ' END IF END DO CLOSE (1) CLOSE (2) STOP 'Job completed.' END PROGRAM Bcs_Tool