program gcg2016_interp
    implicit none

    ! Dynamische Variablen für Rasterdimensionen
    integer :: n_lat, n_lon, total_lines
    real(8) :: lat_min, lat_max, lon_min, lon_max
    real(8) :: dlat, dlon

    ! Dynamic Allocatable Grid (Spalte = lon, Zeile = lat)
    real(8), allocatable :: grid(:,:)

    ! Eingabe- und Hilfsvariablen
    real(8) :: lat_in, lon_in, h_ellip
    real(8) :: lat_read, lon_read, val_read
    real(8) :: zeta, h_dhhn2016
    real(8) :: t, u
    integer :: i, j, row, col, ios

    ! 1. Eingabewerte festlegen (München-Koordinaten)
    lat_in  = 48.16989345_8
    lon_in  = 11.55152710_8
    h_ellip = 610.77_8

    ! 2. Erster Durchlauf: Header & Geometrie der Datei ermitteln
    open(unit=10, file='GCG2016.txt', status='old', action='read', iostat=ios)
    if (ios /= 0) then
        print *, 'Fehler: GCG2016.txt konnte nicht geoffnet werden!'
        stop
    end if

    ! Erste Zeile lesen -> Startwerte setzen
    read(10, *, iostat=ios) lat_read, lon_read, val_read
    lat_max = lat_read
    lon_min = lon_read
    lat_min = lat_read
    lon_max = lon_read

    total_lines = 1
    n_lon = 1

    ! Spaltenanzahl (n_lon) ermitteln: solange Breite konstant bleibt
    do
        read(10, *, iostat=ios) lat_read, lon_read, val_read
        if (ios /= 0) exit
        total_lines = total_lines + 1

        if (abs(lat_read - lat_max) < 1.0e-5_8) then
            n_lon = n_lon + 1
            if (lon_read > lon_max) lon_max = lon_read
        else
            ! Erste Zeile beendet, Rest der Datei durchzählen
            exit
        end if
    end do

    ! Restliche Zeilen zählen, um lat_min zu finden
    do
        read(10, *, iostat=ios) lat_read, lon_read, val_read
        if (ios /= 0) exit
        total_lines = total_lines + 1
        if (lat_read < lat_min) lat_min = lat_read
        if (lon_read > lon_max) lon_max = lon_read
    end do
    close(10)

    n_lat = total_lines / n_lon
    dlat = (lat_max - lat_min) / dble(n_lat - 1)
    dlon = (lon_max - lon_min) / dble(n_lon - 1)

    print '(A, I7, A, I5, A, I5)', 'Grid eingelesen: ', total_lines, ' Punkte (', n_lon, ' x ', n_lat, ')'

    allocate(grid(n_lon, n_lat))

    ! 3. Zweiter Durchlauf: Matrix im RAM befüllen (Nord -> Süd, West -> Ost)
    open(unit=10, file='GCG2016.txt', status='old', action='read')
    do row = n_lat, 1, -1
        do col = 1, n_lon
            read(10, *) lat_read, lon_read, val_read
            grid(col, row) = val_read
        end do
    end do
    close(10)

    ! 4. Grid-Index für die gesuchte Koordinate bestimmen
    i = int((lon_in - lon_min) / dlon) + 1
    j = int((lat_in - lat_min) / dlat) + 1

    if (i < 1 .or. i >= n_lon .or. j < 1 .or. j >= n_lat) then
        print *, 'Fehler: Koordinaten liegen ausserhalb des GCG2016-Rasters!'
        deallocate(grid)
        stop
    end if

    ! Normierte Faktoren im Quadrat [0, 1]
    t = (lon_in - (lon_min + (i - 1) * dlon)) / dlon
    u = (lat_in - (lat_min + (j - 1) * dlat)) / dlat

    ! 5. Bilineare Interpolation für Zeta
    zeta = (1.0_8 - t) * (1.0_8 - u) * grid(i, j)     + &
           t           * (1.0_8 - u) * grid(i+1, j)   + &
           (1.0_8 - t) * u           * grid(i, j+1)   + &
           t           * u           * grid(i+1, j+1)

    if (zeta > 90000.0_8) then
        print *, 'Fehler: Punkt liegt im Bereich ohne Geoiddaten (999999 / NoData)!'
        deallocate(grid)
        stop
    end if

    ! 6. DHHN2016 Normalhöhe berechnen: H = h - zeta
    h_dhhn2016 = h_ellip - zeta

    ! 7. Ausgabe
    print '(A)',                  '=========================================='
    print '(A)',                  ' GCG2016 Transformation -> DHHN2016'
    print '(A)',                  '=========================================='
    print '(A, F12.8, A, F12.8)', 'Koordinate:            ', lat_in, ' N, ', lon_in, ' E'
    print '(A, F10.3, A)',        'Ellipsoidische Hohe h: ', h_ellip, ' m'
    print '(A, F10.4, A)',        'Quasigeoidhohe Zeta:   ', zeta, ' m'
    print '(A, F10.3, A)',        'Normalhohe H DHHN2016: ', h_dhhn2016, ' m'
    print '(A)',                  '=========================================='

    deallocate(grid)
end program gcg2016_interp