PROGRAM Complete_Type2_bNKI !Takes nodal velocities from a successful regional NeoKinema model, !and uses them to complete type-2 entries in a b*.NKI !boundary-conditions file for a local NeoKinema (or SHELLS) model. !Note that nodes are matched by location, not by number. !CAUTION: This logic will NOT work with nodal-velocity !files produced by Shells, because in any Shells model !containing faults there will be multiple co-located nodes !with different velocities along faults. !By Peter Bird, UCLA, 2010.02.13 & 2013.11.21 IMPLICIT NONE CHARACTER*80 :: old_bNKI_pathfile, new_bNKI_pathfile, old_FEG_pathfile, old_vOUT_pathfile, title1, title2, title3 CHARACTER*80, DIMENSION(:), ALLOCATABLE :: BC_lines INTEGER :: best_node, count_type2, i, ios, j, n1, n2, n3, ncond, numnod REAL*8 :: azimuth_degrees, best_Dds, Delta_degrees_squared, latitude, longitude REAL*8, DIMENSION(:), ALLOCATABLE :: node_lon, node_lat DOUBLE PRECISION :: v_in_mps DOUBLE PRECISION, DIMENSION (:,:), ALLOCATABLE :: vw WRITE (*, "(' PROGRAM Complete_Type2_bNKI:')") WRITE (*, "(' Takes nodal velocities from a successful regional NeoKinema model,')") WRITE (*, "(' and uses them to complete type-2 entries in a b*.NKI')") WRITE (*, "(' boundary-conditions file for a local NeoKinema (or Shells) model.')") WRITE (*, "(' Note that nodes are matched by location, not by number.')") WRITE (*, "(' CAUTION: This logic will NOT work with nodal-velocity')") WRITE (*, "(' files produced by Shells, because in any Shells model')") WRITE (*, "(' containing faults there will be multiple co-located nodes')") WRITE (*, "(' with different velocities along faults.')") WRITE (*, "(' By Peter Bird, UCLA, 2010.02.13 & 2013.11.21')") WRITE (*, *) 10 WRITE (*, "(' Enter [Drive:][\path\]filename for the v*.OUT file' & & /' produced by the large NeoKinema model:')") READ (*, "(A)") old_vOUT_pathfile OPEN (UNIT = 1, FILE = TRIM(old_vOUT_pathfile), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: File not found (at this location). Try again...')") CALL Pause() GO TO 10 END IF READ (1, "(A)") title1 READ (1, "(A)") title2 READ (1, "(A)") title3 WRITE (*, "(' Here are the 3 header/title lines from that file:')") WRITE (*, "(' ', A)") TRIM(title1) WRITE (*, "(' ', A)") TRIM(title2) WRITE (*, "(' ', A)") TRIM(title3) WRITE (*, *) 20 WRITE (*, "(' Enter [Drive:][\path\]filename for the .FEG file' & & /' used in this large NeoKinema model:')") READ (*, "(A)") old_FEG_pathfile OPEN (UNIT = 2, FILE = TRIM(old_FEG_pathfile), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: File not found (at this location). Try again...')") CALL Pause() GO TO 20 END IF READ (2, *) title1 WRITE (*, "(' Here is the header/title line from that file:')") WRITE (*, "(' ', A)") TRIM(title1) READ (2, *) numnod ALLOCATE ( node_lon(numnod) ); node_lon = 0.0D0 ALLOCATE ( node_lat(numnod) ); node_lat = 0.0D0 DO i = 1, numnod READ (2, *) j, node_lon(j), node_lat(j) END DO CLOSE (2) WRITE (*, *) WRITE (*, "(' Captured node locations in old, regional .FEG file.')") ALLOCATE ( vw(2, numnod) ); vw = 0.0D0 !DO i = 1, numnod ! READ (1, *) vw(1, i), vw(2, i) ! S & E components, in m/s !END DO READ (1, *)((vw(j,i), j=1,2), i=1,numnod) CLOSE (1) WRITE (*, "(' Captured nodal velocities in old, regional v.OUT file.')") WRITE (*, *) 30 WRITE (*, "(' Enter [Drive:\[\path\]filename for incomplete b*.NKI file:')") READ (*, "(A)") old_bNKI_pathfile OPEN (UNIT = 3, FILE = TRIM(old_bNKI_pathfile), STATUS = "OLD", IOSTAT = ios) IF (ios /= 0) THEN WRITE (*, "(' ERROR: File not found (at this location). Try again...')") CALL Pause() GO TO 30 END IF READ (3, *) title1 WRITE (*, "(' Here is the header/title line from that file:')") WRITE (*, "(' ', A)") TRIM(title1) CALL Pause() ncond = 0; count_type2 = 0 counting : DO READ (3, *, IOSTAT = ios) n1, n2, n3 ! essential elements: " 1 1048 2", for example IF (ios == 0) THEN ncond = ncond + 1 IF (n3 == 2) count_type2 = count_type2 + 1 ELSE EXIT counting END IF END DO counting WRITE (*, *) WRITE (*, "(' Counted ',I6,' boundary conditions in this file.')") ncond WRITE (*, "(' of which',I6,' were type-2 BCs needing completion.')") count_type2 CLOSE (3) ALLOCATE ( BC_lines(ncond) ); BC_lines = ' ' OPEN (UNIT = 3, FILE = TRIM(old_bNKI_pathfile), STATUS = "OLD", IOSTAT = ios, PAD = "YES") READ (3, "(A)") title1 DO i = 1, ncond READ (3, "(A)") BC_lines(i) END DO CLOSE (3) WRITE (*, *) 40 WRITE (*, "(' Enter [Drive:\[\path\]filename for the NEW, completed b*.NKI file:')") READ (*, "(A)") new_bNKI_pathfile OPEN (UNIT = 4, FILE = TRIM(new_bNKI_pathfile)) WRITE (*, "(' Enter a revised title line for this new b*.NKI file:')") READ (*, "(A)") title1 WRITE (4, *) TRIM(title1) DO i = 1, ncond READ(BC_lines(i), *) n1, n2, n3 IF (n3 /= 2) THEN ! output this line unchanged WRITE (4, "(A)") TRIM(BC_lines(i)) ELSE ! must complete type-2 BC: READ (BC_lines(i), *) n1, n2, n3, latitude, longitude best_Dds = 999.9D0 best_node = 0 DO j = 1, numnod Delta_degrees_squared = (longitude - node_lon(j))**2 + (latitude - node_lat(j))**2 IF (Delta_degrees_squared < best_Dds) THEN best_Dds = Delta_degrees_squared best_node = j END IF END DO IF ((best_node == 0).OR.(best_Dds > 1.0D0)) THEN WRITE (*, "(' ERROR: Could not find any node near lon ',F9.3,', lat ',F8.3)") longitude, latitude WRITE (*, "(' This line in the .bcs file will not be modified; fix it manually.')") CALL Pause() WRITE (4, "(A)") TRIM(BC_lines(i)) ELSE IF (best_Dds > 0.0001D0) THEN WRITE (*, "(' ERROR: Could not find any node at lon ',F9.3,', lat ',F8.3)") longitude, latitude WRITE (*, "(' The velocity of the nearest node (within 1 degree) was applied.')") v_in_mps = DSQRT(vw(1, best_node)**2 + vw(2, best_node)**2) azimuth_degrees = 57.2957795130D0 * DAtan2F(vw(2, best_node), -vw(1, best_node)) ! (y, x) = (E, N) IF (azimuth_degrees < 0.0D0) azimuth_degrees = azimuth_degrees + 360.0D0 WRITE (4, "(2I6, I5, ES12.4, F8.2, F9.3, F10.3)") n1, n2, n3, v_in_mps, azimuth_degrees, latitude, longitude ELSE ! normal processing (we hope): v_in_mps = DSQRT(vw(1, best_node)**2 + vw(2, best_node)**2) azimuth_degrees = 57.2957795130D0 * DAtan2F(vw(2, best_node), -vw(1, best_node)) ! (y, x) = (E, N) IF (azimuth_degrees < 0.0D0) azimuth_degrees = azimuth_degrees + 360.0D0 WRITE (4, "(2I6, I5, ES12.4, F8.2, F9.3, F10.3)") n1, n2, n3, v_in_mps, azimuth_degrees, latitude, longitude END IF END IF END DO WRITE (*, *) WRITE (*, "(' Job completed.')") CALL Pause() CONTAINS DOUBLE PRECISION FUNCTION DAtan2F (y, x) ! corrects for possible abend in case of (0., 0.) IMPLICIT NONE DOUBLE PRECISION, INTENT(IN) :: x, y IF ((x /= 0.0D0).OR.(y /= 0.0D0)) THEN DAtan2F = DATAN2(y, x) ELSE DAtan2F = 0.0D0 END IF END FUNCTION DAtan2F SUBROUTINE Pause() IMPLICIT NONE WRITE (*,"(' Press [Enter]...'\)") READ (*,*) END SUBROUTINE Pause END PROGRAM Complete_Type2_bNKI