PROGRAM HFgrid14_2_GRD
    !Utility program to reformat file HFgrid14.csv of Lucazeau [2019, G**3]
    !into Peter Bird's .GRD format.
    !Note that I had to convert every ';' in the original file to ','
    !to get it to open in Excel, so that altered version HFgrid14-commas.csv
    !is the file used here.
    !The output is actually TWO .GRD files, corresponding to the two (full)
    !columns in the input file: HFgrid14.grd corresponds to column "HF_pred"
    !while sHFgrid14.grd corresponds to column "sHF_pred".
    !By Peter Bird, UCLA, 2022.12.03

    IMPLICIT NONE
    CHARACTER*80 :: inputFile = "HFgrid14-commas.csv"
    REAL*8, PARAMETER :: dElon = 0.5D0, dNLat = 0.5D0
    REAL*8, PARAMETER :: ELon_min = -179.75D0, ELon_max = +179.75D0, NLat_min = -89.75D0, NLat_max = +89.75D0
    INTEGER :: i, ios, j, lines, nX, nY
    REAL*8 :: ELon, NLat, HF_pred, sHF_pred
    REAL*8, DIMENSION(:,:), ALLOCATABLE :: HF_GRD, sHF_GRD
    
    nX = NINT(1.0D0 + (ELon_max - ELon_min) / dElon)
    nY = NINT(1.0D0 + (NLat_max - NLat_min) / dNLat)
    ALLOCATE ( HF_GRD(nY, nX) )
    ALLOCATE ( sHF_GRD(nY, nX) )
    HF_GRD  = 0.0D0
    sHF_GRD = 0.0D0

    OPEN (UNIT = 1, FILE = TRIM(inputFile), STATUS = "OLD", IOSTAT = ios)
    IF (ios /= 0) THEN
        WRITE (*, "(' ERROR: Input file ', A, ' not found in current folder.')")
        CALL Pause()
        STOP
    END IF
    READ (1, *) ! skip header line
    lines = 1
    reading: DO
        lines = lines + 1
        READ(1, *, IOSTAT = ios) ELon, NLat, HF_pred, sHF_pred
        IF (ios /= 0) EXIT reading ! at EOF, I hope
        i = NINT(1.0D0 + (NLat_max - NLat) / dNLat) ! row (starting at N)
        j = NINT(1.0D0 + (ELon - ELon_min) / dELon) ! column (starting in W)
        HF_pred = MAX(HF_pred, 30.0D0) ! limit low values in old seafloor
        sHF_pred = MAX(sHF_pred, 30.0D0) ! limit low values in old seafloor
        HF_pred = MIN(HF_pred, 300.0D0) ! limit high values along spreading ridges
        sHF_pred = MIN(sHF_pred, 300.0D0) ! limit high values along spreading ridges
        HF_GRD(i, j) = HF_pred * 1.0D-3  ! convert to SI units = W/m**2
        sHF_GRD(i, j) = sHF_pred * 1.0D-3 ! convert to SI units = W/m**2
    END DO reading
    CLOSE (1)

    OPEN (UNIT = 2, FILE = "HFgrid14.grd") ! unconditional OPEN; overwrites any old version
    WRITE (2, "(3F10.2, ' HFgrid14.grd, reformatted from HFgrid14.csv of Lucazeau, 2019')") ELon_min, dELon, ELon_max
    WRITE (2, "(3F10.2)") NLat_min, dNLat, NLat_max
    WRITE (2, "(10F8.4)") ((HF_GRD(i, j), j = 1, nX), i = 1, nY)
    CLOSE (2)

    OPEN (UNIT = 3, FILE = "sHFgrid14.grd") ! unconditional OPEN; overwrites any old version
    WRITE (3, "(3F10.2, ' sHFgrid14.grd, reformatted from HFgrid14.csv of Lucazeau, 2019')") ELon_min, dELon, ELon_max
    WRITE (3, "(3F10.2)") NLat_min, dNLat, NLat_max
    WRITE (3, "(10F8.4)") ((sHF_GRD(i, j), j = 1, nX), i = 1, nY)
    CLOSE (3)

    WRITE (*, *)
    WRITE (*, "(' Job completed.')")
    CALL Pause()

CONTAINS

    SUBROUTINE Pause()
       IMPLICIT NONE
       WRITE (*,"(' Press [Enter]...'\)")
       READ  (*,*)
    END SUBROUTINE Pause
END PROGRAM HFgrid14_2_GRD