SUBROUTINE ReadATI3D( FileRoot, TopATI, DepthT, PRTFile )
CHARACTER (LEN= 1), INTENT( IN ) :: TopATI ! Set to '~' if altimetry is not flat
INTEGER, INTENT( IN ) :: PRTFile ! unit number for print file
REAL (KIND=8), INTENT( IN ) :: DepthT ! Nominal top depth
CHARACTER (LEN=80), INTENT( IN ) :: FileRoot
REAL (KIND=8), ALLOCATABLE :: Temp( : )
REAL (KIND=8), ALLOCATABLE :: TopGlobalx( : ), TopGlobaly( : )
SELECT CASE ( TopATI )
CASE ( '~', '*' )
WRITE( PRTFile, * ) '__________________________________________________________________________'
WRITE( PRTFile, * )
WRITE( PRTFile, * ) 'Using top-altimetry file'
OPEN( UNIT = ATIFile, FILE = TRIM( FileRoot ) // '.ati', STATUS = 'OLD', IOSTAT = IOStat, ACTION = 'READ' )
IF ( IOsTAT /= 0 ) THEN
WRITE( PRTFile, * ) 'ATIFile = ', TRIM( FileRoot ) // '.ati'
CALL ERROUT( 'ReadATI', 'Unable to open altimetry file' )
END IF
READ( ATIFile, * ) atiType
SELECT CASE ( atiType )
CASE ( 'R' )
WRITE( PRTFile, * ) 'Regular grid for a 3D run'
CASE ( 'C' )
WRITE( PRTFile, * ) 'Regular grid for a 3D run (curvilinear)'
CASE DEFAULT
CALL ERROUT( 'ReadATI3D', 'Unknown option for selecting altimetry interpolation' )
END SELECT
! x values
READ( ATIFile, * ) NatiPts( 1 )
WRITE( PRTFile, * )
WRITE( PRTFile, * ) 'Number of altimetry points in x', NatiPts( 1 )
ALLOCATE( TopGlobalx( MAX( NatiPts( 1 ), 3 ) ), Stat = IAllocStat )
IF ( IAllocStat /= 0 ) &
CALL ERROUT( 'BELLHOP3D:ReadATI3D', 'Insufficient memory for altimetry data: reduce # ati points' )
TopGlobalx( 3 ) = -999.9
READ( ATIFile, * ) TopGlobalx( 1 : NatiPts( 1 ) )
CALL SubTab( TopGlobalx, NatiPts( 1 ) )
WRITE( PRTFile, "( 5G14.6 )" ) ( TopGlobalx( ix ), ix = 1, MIN( NatiPts( 1 ), Number_to_Echo ) )
IF ( NatiPts( 1 ) > Number_to_Echo ) WRITE( PRTFile, "( G14.6 )" ) ' ... ', TopGlobalx( NatiPts( 1 ) )
IF ( .NOT. monotonic( TopGlobalx, NatiPts( 1 ) ) ) THEN
CALL ERROUT( 'BELLHOP3D:ReadATI3D', 'Altimetry x-coordinates are not monotonically increasing' )
END IF
! y values
READ( ATIFile, * ) NatiPts( 2 )
WRITE( PRTFile, * )
WRITE( PRTFile, * ) 'Number of altimetry points in y', NatiPts( 2 )
ALLOCATE( TopGlobaly( MAX( NatiPts( 2 ), 3 ) ), Stat = IAllocStat )
IF ( IAllocStat /= 0 ) &
CALL ERROUT( 'BELLHOP3D:ReadATI3D', 'Insufficient memory for altimetry data: reduce # ati points' )
TopGlobaly( 3 ) = -999.9
READ( ATIFile, * ) TopGlobaly( 1 : NatiPts( 2 ) )
CALL SubTab( TopGlobaly, NatiPts( 2 ) )
WRITE( PRTFile, "( 5G14.6 )" ) ( TopGlobaly( iy ), iy = 1, MIN( NatiPts( 2 ), Number_to_Echo ) )
IF ( NatiPts( 2 ) > Number_to_Echo ) WRITE( PRTFile, "( G14.6 )" ) ' ... ', TopGlobaly( NatiPts( 2 ) )
IF ( .NOT. monotonic( TopGlobaly, NatiPts( 2 ) ) ) THEN
CALL ERROUT( 'BELLHOP3D:ReadATI3D', 'Altimetry y-coordinates are not monotonically increasing' )
END IF
TopGlobalx = 1000. * TopGlobalx ! convert km to m
TopGlobaly = 1000. * TopGlobaly
! z values
ALLOCATE( Top( NatiPts( 1 ), NatiPts( 2 ) ), Temp( NatiPts( 1 ) ), Stat = IAllocStat )
IF ( IAllocStat /= 0 ) &
CALL ERROUT( 'BELLHOP3D:ReadATI3D', 'Insufficient memory for altimetry data: reduce # ati points' )
WRITE( PRTFile, * )
DO iy = 1, NatiPts( 2 )
READ( ATIFile, * ) Top( :, iy )%x( 3 ) ! read a row of depths
! IF ( iy < Number_to_Echo .OR. iy == NatiPts( 2 ) ) THEN ! echo some values
! WRITE( PRTFile, FMT = "(G11.3)" ) Top( :, iy )%x( 3 )
! END IF
! IF ( ANY( Top( :, iy )%x( 3 ) > DepthB ) ) THEN
! CALL ERROUT( 'BELLHOP3D:ReadATI3D', 'Altimetry drops below lowest point in the sound speed profile' )
! END IF
END DO
CLOSE( ATIFile )
IF ( ANY( ISNAN( Top( :, : )%x( 3 ) ) ) ) THEN
WRITE( PRTFile, * ) 'Warning in BELLHOP3D - ReadATI3D : The altimetry file contains a NaN'
END IF
DO ix = 1, NatiPts( 1 )
DO iy = 1, NatiPts( 2 )
Top( ix, iy )%x( 1 ) = TopGlobalx( ix )
Top( ix, iy )%x( 2 ) = TopGlobaly( iy )
END DO
END DO
DEALLOCATE( TopGlobalx )
DEALLOCATE( TopGlobaly )
CALL ComputeBdryTangentNormal( Top, 'Top' )
CASE DEFAULT ! no altimetry given, use SSP depth for flat top
atiType = 'R'
NatiPts = [ 2, 2 ]
ALLOCATE( Top( 2, 2 ), Stat = IAllocStat )
IF ( IAllocStat /= 0 ) CALL ERROUT( 'BELLHOP3D:ReadATI3D', 'Insufficient memory' )
Top( 1, 1 )%x = [ -big, -big, DepthT ]
Top( 1, 2 )%x = [ -big, big, DepthT ]
Top( 2, 1 )%x = [ big, -big, DepthT ]
Top( 2, 2 )%x = [ big, big, DepthT ]
Top( 1, 1 )%t = [ 1.0, 0.0, 0.0 ] ! tangent to top
Top( 1, 1 )%n1 = [ 0.0, 0.0, -1.0 ] ! outward-pointing normal
Top( 1, 1 )%n2 = [ 0.0, 0.0, -1.0 ] ! outward-pointing normal
Top( 1, 2 )%t = [ 1.0, 0.0, 0.0 ] ! tangent to top
Top( 1, 1 )%n1 = [ 0.0, 0.0, -1.0 ] ! outward-pointing normal
Top( 1, 1 )%n2 = [ 0.0, 0.0, -1.0 ] ! outward-pointing normal
Top( 2, 1 )%t = [ 1.0, 0.0, 0.0 ] ! tangent to top
Top( 2, 1 )%n1 = [ 0.0, 0.0, -1.0 ] ! outward-pointing normal
Top( 2, 1 )%n2 = [ 0.0, 0.0, -1.0 ] ! outward-pointing normal
Top( 2, 2 )%t = [ 1.0, 0.0, 0.0 ] ! tangent to top
Top( 2, 2 )%n1 = [ 0.0, 0.0, -1.0 ] ! outward-pointing normal
Top( 2, 2 )%n2 = [ 0.0, 0.0, -1.0 ] ! outward-pointing normal
END SELECT
! dummy TopSeg info to force GetTopSeg to search for the active segment on first call
xTopSeg = [ +big, -big ]
yTopSeg = [ +big, -big ]
END SUBROUTINE ReadATI3D
