C OMEGAXR.FOR: C C A convenient utility program to compute present velocity C vectors from a plate-tectonic model. C C Peter Bird, Earth & Space Sciences, UCLA C 25 April 1995 C 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: ',F10.2) C C Internally, computations are done in SI: C RADIUS=RADIUS*1.E3 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.LT.1).OR.(ITYPE.GT.2)) GO TO 25 IF (ITYPE.EQ.1) THEN FACTOR=(3.14159/180.)/(1.E6*3.15576E7) ELSEIF (ITYPE.EQ.2) THEN FACTOR=1. 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' + /' ------------------------------------------------' + ) ENDIF C C Enter pole, or list of poles to be summed. C ("FACTOR" will convert rotation rates to radians/second = SI) C 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. SUMY=0. SUMZ=0. 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 = +): ',F10.3) IF ((POLLAT.GT.90.).OR.(POLLAT.LT.-90.)) GO TO 42 POLLAT=POLLAT*(3.14159/180.) 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 = +): ',F10.3) POLLON=POLLON*(3.14159/180.) 469 WRITE (*,47) I 47 FORMAT (' For pole #',I1, + ', ROTATION RATE (counterclockwise = +): '\) READ (*,*,IOSTAT=IOS) ROTRAT IF (IOS.NE.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,E12.4) ROTRAT=ROTRAT*FACTOR SUMX=SUMX+ROTRAT*COS(POLLAT)*COS(POLLON) SUMY=SUMY+ROTRAT*COS(POLLAT)*SIN(POLLON) SUMZ=SUMZ+ROTRAT*SIN(POLLAT) 50 CONTINUE 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 WRITE (1,60) POLLAT, POLLON, POLLN2, ROTRAT 60 FORMAT (' ---------------------------------------' + /' Net 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.EQ.1) THEN WRITE (*,61) WRITE (1,61) 61 FORMAT (' in units of degrees per million years.') ELSEIF (ITYPE.EQ.2) THEN WRITE (*,62) WRITE (1,62) 62 FORMAT (' in units of radians per second.') ENDIF C C Velocity computations C 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.GT.90.).OR.(VELLAT.LT.-90.)) GO TO 101 VELLAT=VELLAT*(3.14159/180.) 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 = +): ', + F10.3) VELLON=VELLON*(3.14159/180.) RX=RADIUS*COS(VELLAT)*COS(VELLON) RY=RADIUS*COS(VELLAT)*SIN(VELLON) RZ=RADIUS*SIN(VELLAT) VX=SUMY*RZ-SUMZ*RY VY=SUMZ*RX-SUMX*RZ VZ=SUMX*RY-SUMY*RX VEL=SQRT(VX**2+VY**2+VZ**2) VEL2=VEL*1000.*3.15576E7 C Find unit vectors pointing North and East to extract components UNX= -SIN(VELLAT)*COS(VELLON) UNY= -SIN(VELLAT)*SIN(VELLON) UNZ=COS(VELLAT) UEX= -SIN(VELLON) UEY=COS(VELLON) 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.LT.0.0) AZIM=AZIM+360. VE2=VE*1000.*3.15576E7 VN2=VN*1000.*3.15576E7 WRITE (*,150) VEL, VEL2, AZIM, VE2, VN2 WRITE (1,150) VEL, VEL2, AZIM, VE2, VN2 150 FORMAT (' Computed velocity is ',1P,E10.3' meters per second,' + /' or ',0P,F10.3,' millimeters per year,' + /' at an azimuth of ',F5.1,' degrees clockwise' + ,' from North,' + /' which gives components of' + /' ',0P,F10.3,' mm/a Eastward and ',F10.3,' mm/a' + ,' Northward.') C C LOOP or STOP? C 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.LE.0) STOP ' See OMEGAXR.LOG for record of your work.' IF (NOW.EQ.1) GO TO 101 IF (NOW.GE.2) GO TO 35 END