PROGRAM OmegaXR ! A convenient utility program to compute present velocity ! vectors from a plate-tectonic model. ! Peter Bird, Dept. of Earth, Planetary, & Space Sciences, UCLA ! 25 April 1995; updated to .f90 on 2015.07.19. IMPLICIT DOUBLE PRECISION (A-H, O-Z) IMPLICIT INTEGER (I-N) REAL*8, PARAMETER :: Two_Pi = 6.283185307179586D0 REAL*8, PARAMETER :: Pi = 3.141592653589793D0 REAL*8, PARAMETER :: Pi_over_2 = 1.5707963267948965D0 REAL*8, PARAMETER :: Pi_over_4 = 0.78539816339744825D0 REAL*8, PARAMETER :: degrees_per_radian = 57.2957795130D0 REAL*8, PARAMETER :: radians_per_degree = 0.01745329252D0 REAL*8, PARAMETER :: s_per_year = 365.25 * 24.0D0 * 60.0D0 * 60.0D0 OPEN (UNIT = 1, FILE = 'OmegaXR.log') WRITE (*, 10) WRITE (1, 11) 10 FORMAT ( & & ' --== OmegaXR ==--:' & & //' A convenient utility program to compute present velocity' & & ,' vectors' & & /' from a plate-tectonic model.'/) 11 FORMAT (' OmegaXr:' & & /' A convenient utility program to compute present velocity' & & ,' vectors' & & /' from a plate-tectonic model.'/) WRITE (*, 20) 20 FORMAT (' Enter RADIUS of planet in kilometers: '\) READ (*, * ) radius WRITE (1, 21) radius 21 FORMAT (' Enter RADIUS of planet in kilometers: ',F13.5) ! Internally, computations are done in SI: radius = radius * 1.0D3 25 WRITE (*, 30) 30 FORMAT (/' Select unit identifier for rotation rates:' & & /' 1 = degrees per million years;' & & /' 2 = radians per second.' & & /' Enter INTEGER (1-2): '\) READ (*, * ) itype WRITE (1, 31) itype 31 FORMAT (/' Select unit identifier for rotation rates:' & & /' 1 = degrees per million years;' & & /' 2 = radians per second.' & & /' Enter INTEGER (1-2): ',I1) IF ((itype < 1).OR.(itype > 2)) GO TO 25 IF (itype == 1) THEN factor = radians_per_degree / (1.0D6 * s_per_year) elseif (itype == 2) THEN factor = 1.0D0 WRITE (*, 32) WRITE (1, 32) 32 FORMAT (/' ------------------------------------------------' & & /' NOTE: -16' & & /' To enter very small numbers, like 1.234x10 ' & & ,' radians/second' & & /' on one line, use exponential format: 1.234E-16' & & /' ------------------------------------------------' & & ) END IF ! Enter pole, or list of poles to be summed. ! ("FACTOR" will convert rotation rates to radians/second = SI) 35 WRITE (*, 40) 40 FORMAT (/' Enter NUMBER of rotation poles to be summed (1-?): '\) READ (*, * ) npole WRITE (1, 41) npole 41 FORMAT (/' Enter NUMBER of rotation poles to be summed (1-?): ', & & I2) sumx = 0.0D0 sumy = 0.0D0 sumz = 0.0D0 DO 50 i = 1, npole 42 WRITE (*, 43) i 43 FORMAT (/' For pole #',I1,', LATITUDE in ' & & ,'degrees (N = +): '\) READ (*, * ) pollat WRITE (1, 44) i, pollat 44 FORMAT (/' For pole #',I1,', LATITUDE in ' & & ,'degrees (N = +): ',F12.5) IF ((pollat > 90.0D0).OR.(pollat < -90.0D0)) GO TO 42 pollat = pollat * radians_per_degree WRITE (*, 45) i 45 FORMAT (' For pole #',I1,', LONGITUDE in ' & & ,'degrees (E = +): '\) READ (*, * ) pollon WRITE (1, 46) i, pollon 46 FORMAT (' For pole #',I1,', LONGITUDE in ' & & ,'degrees (E = +): ',F12.5) pollon = pollon * radians_per_degree 469 WRITE (*, 47) i 47 FORMAT (' For pole #',I1, & & ', ROTATION RATE (counterclockwise = +): '\) READ (*, * , iostat = ios) rotrat IF (ios /= 0) THEN WRITE (*, "(' ERROR: Number format not understood.' & & /' Press Enter to try again...')") READ (*, * ) GO TO 469 END IF WRITE (1, 48) i, rotrat 48 FORMAT (' For pole #',I1, & & ', ROTATION RATE (counterclockwise = +): ', & & 1P,E13.5) rotrat = rotrat * factor sumx = sumx + rotrat * DCOS(pollat) * DCOS(pollon) sumy = sumy + rotrat * DCOS(pollat) * DSIN(pollon) sumz = sumz + rotrat * DSIN(pollat) 50 CONTINUE rotrat = DSQRT(sumx**2 + sumy**2 + sumz**2) upx = sumx / rotrat upy = sumy / rotrat upz = sumz / rotrat rotrat = rotrat / factor pollon = DATAN2(upy, upx) equate = DSQRT(upx**2 + upy**2) pollat = DATAN2(upz, equate) pollon = pollon * degrees_per_radian pollat = pollat * degrees_per_radian pollon = DMOD(pollon + 720.0D0, 360.0D0) polln2 = 360.0D0 - pollon WRITE (*, 60) pollat, pollon, polln2, rotrat WRITE (1, 60) pollat, pollon, polln2, rotrat 60 FORMAT (' ---------------------------------------' & & /' Net rotation pole is located at ',F9.5,' degrees North' & & /' and ',F9.5,' degrees East (',F9.5,' degrees West)' & & /' with rotation rate of ',1P,E12.5) IF (itype == 1) THEN WRITE (*, 61) WRITE (1, 61) 61 FORMAT (' in units of degrees per million years.') elseif (itype == 2) THEN WRITE (*, 62) WRITE (1, 62) 62 FORMAT (' in units of radians per second.') END IF ! Velocity computations 101 WRITE (*, 110) 110 FORMAT (/' Enter LATITUDE for velocity, in degrees (N = +): '\) READ (*, * ) vellat WRITE (1, 111) vellat 111 FORMAT (/' Enter LATITUDE for velocity, in degrees (N = +): ', & & F10.3) IF ((vellat > 90.0D0).OR.(vellat < -90.0D0)) GO TO 101 vellat = vellat * radians_per_degree WRITE (*, 120) 120 FORMAT (' Enter LONGITUDE for velocity, in degrees (E = +): '\) READ (*, * ) vellon WRITE (1, 121) vellon 121 FORMAT (' Enter LONGITUDE for velocity, in degrees (E = +): ', & & F12.5) vellon = vellon * radians_per_degree rx = radius * DCOS(vellat) * DCOS(vellon) ry = radius * DCOS(vellat) * DSIN(vellon) rz = radius * DSIN(vellat) vx = sumy * rz - sumz * ry vy = sumz * rx - sumx * rz vz = sumx * ry - sumy * rx vel = DSQRT(vx**2 + vy**2 + vz**2) vel2 = vel * 1000.0D0 * s_per_year ! Find unit vectors pointing North and East to extract components unx = -DSIN(vellat) * DCOS(vellon) uny = -DSIN(vellat) * DSIN(vellon) unz = DCOS(vellat) uex = -DSIN(vellon) uey = DCOS(vellon) uez = 0.0D0 vn = vx * unx + vy * uny + vz * unz ve = vx * uex + vy * uey + vz * uez azim = 90.0D0 - degrees_per_radian * DATAN2(vn, ve) IF (azim < 0.0D0) azim = azim + 360.0D0 ve2 = ve * 1000.0D0 * s_per_year vn2 = vn * 1000.0D0 * s_per_year WRITE (*, 150) vel, vel2, azim, ve2, vn2 WRITE (1, 150) vel, vel2, azim, ve2, vn2 150 FORMAT (' Computed velocity is ',1P,E12.5' meters per second,' & & /' or ',0P,F12.5,' millimeters per year,' & & /' at an azimuth of ',F9.5,' degrees clockwise' & & ,' from North,' & & /' which gives components of' & & /' ',0P,F12.5,' mm/a Eastward and ',F12.5,' mm/a' & & ,' Northward.') ! LOOP or STOP? WRITE (*, 160) 160 FORMAT (/' ENTER 0 TO QUIT, 1 FOR NEW VELOCITY, 2 FOR' & & , ' NEW POLE(S): '\) READ (*, * ) now WRITE (1, 161) now 161 FORMAT (/' Enter 0 to quit, 1 for new velocity, 2 for' & & ,' new pole(s): ',I1) IF (now <= 0) STOP ' See OmegaXR.log for record of your work.' IF (now == 1) GO TO 101 IF (now >= 2) GO TO 35 END PROGRAM OmegaXR