PROGRAM BCs_tool C fixed-form Fortran 90 source code by Peter Bird, UCLA, 7/98 CHARACTER*80 :: INFILE, OUTFIL, LINE, LINE2 INTEGER :: BCTYPE REAL :: LAT, LON C C INTRODUCTION C 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: the Africa-fixed version of NUVEL-1.)' +//' Press [Enter] to continue...'\) READ (*, "(A)", IOSTAT = IOS) LINE C C INPUT AND OUTPUT FILE NAMES C 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 .NE. 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 .NE. 0) THEN WRITE (*,"(/' ERROR: Supply a NEW filename.')") GO TO 4 END IF C C GET POLE POSITION AND RATE C WRITE (*,20) 20 FORMAT (/' Enter radius of planet, in kilometers: '\) READ (*,*) RADIUS C C INTERNALLY, COMPUTATIONS ARE DONE IN SI C 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.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) 32 FORMAT (/' NOTE: -16' + /' To enter very small numbers, like 5.362x10 ' + ,' radians/second' + /' on one line, use this format: 5.362E-16'/) ENDIF C C ENTER POLE C ("FACTOR" WILL CONVERT ROTATION RATES TO RADIANS/SECOND = SI) C 42 WRITE (*,43) 43 FORMAT (/' Absolute rotation pole LATITUDE in ' + ,'degrees (N = +): '\) READ (*,*) POLLAT IF ((POLLAT.GT.90.).OR.(POLLAT.LT.-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.EQ.1) THEN WRITE (*,61) 61 FORMAT (' in units of degrees per million years.') ELSEIF (ITYPE.EQ.2) THEN WRITE (*,62) 62 FORMAT (' in units of radians per second.') ENDIF C C SCAN FILE, COMPUTING VELOCITIES FOR TYPE-3 NODES C READ (1, "(A)") LINE WRITE (*,"(/' Title line from input .BCS file follows:')") WRITE (*,"(' ',A)") TRIM(LINE) WRITE (*,"(/' Enter new title line, or just [Enter] top 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 .LT. 0) EXIT READ (LINE, *) I, NODE, BCTYPE IF (BCTYPE .EQ. -1) THEN READ (LINE, *, IOSTAT = IOS) I,NODE,BCTYPE, + VELLAT,VELLON IF (IOS .NE. 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 .EQ. 0) THEN READ (LINE, *, IOSTAT = IOS) I,NODE,BCTYPE, + VELLAT,VELLON IF (IOS .NE. 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 .EQ. 1) THEN READ (LINE, *, IOSTAT = IOS) I,NODE,BCTYPE,VEL,AZIM, + VELLAT,VELLON IF (IOS .NE. 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 .EQ. 2) THEN READ (LINE, *, IOSTAT = IOS) I,NODE,BCTYPE,VEL,AZIM, + VELLAT,VELLON IF (IOS .NE. 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 .EQ. 3) THEN READ (LINE, *, IOSTAT = IOS) I,NODE,BCTYPE,VELLAT,VELLON IF (IOS .NE. 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) C 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.LT.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