MODULE input
    USE definitions
    IMPLICIT NONE
    
CONTAINS    
    SUBROUTINE read_data(alpha, northerness, p, x, obs)
      ! This subroutine reads the data.
      IMPLICIT NONE
      INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(8, 0)
      TYPE(variables_real), INTENT(INOUT) :: alpha
      TYPE(variables_int), INTENT(INOUT) :: northerness      
      TYPE(variables_real), INTENT(INOUT) :: p
      TYPE(variables_int), INTENT(INOUT) :: x
      TYPE(variables_int), INTENT(INOUT) :: obs    
      
      INTEGER :: i, j, temp, next, iostatus

      ! this is used for reading the whole datafile
      ! when there's an one-dimensional variable to read
      ! (to get the correct column)
      ! for this, we need to know the dimensions of data files (hard-coded)
      REAL (KIND=dp), ALLOCATABLE, DIMENSION(:) :: datafile
      TYPE(variable_int) :: current

      ! HARD-CODED: the data file name(s)
      ! these have been generated: the observations_1 file has zeros on the
      ! places where the original had *-signs. This makes reading it easier.
      ! the observations_missing has the places of the missing variables
      ! in a matrix of boolean values (0, 1).
    
      CHARACTER(LEN=*), PARAMETER :: obs_filename='observations_1.txt'
      CHARACTER(LEN=*), PARAMETER :: obs_missing='observations_missing.txt'

      CHARACTER(LEN=*), PARAMETER :: northerness_filename='squares_1.txt'

      CHARACTER(LEN=*), PARAMETER :: northerness_missing='squares_missing.txt'

      OPEN(UNIT=16, FILE=obs_filename, STATUS='old', IOSTAT=iostatus)

      IF (iostatus /= 0) THEN
          WRITE (*, *) 'Could not open ', obs_filename, '. Exiting program.'
          STOP
      END IF

      OPEN(UNIT=17, FILE=obs_missing, STATUS='old', IOSTAT=iostatus)

      IF (iostatus /= 0) THEN
          WRITE (*, *) 'Could not open ', obs_missing, '. Exiting.'
          STOP
      END IF

      next = 1
    
    ! these loop lengths come from the generator. 
    ! These are the dimensions of the data array (obs), and should match
    ! the dimensions of the matrix in the file.

    ! HARD-CODED: the loop length
    DO i = 1, 21

       !READ(16, *, ERR=13, END=14) obs % two_dim(i, :) % value
       
! if two-dimensional and spans the whole matrix, reading is easy
       READ(16, *, IOSTAT=iostatus) obs % two_dim(i, :) % value
       
       IF (iostatus /= 0) THEN
           WRITE (*, *) 'Could not read ', obs_filename, '. Exiting program.'
           STOP
       END IF

       DO j=1, 4

          IF (j == 4) THEN
             READ(17, FMT='(I2)', IOSTAT=iostatus) temp
          ELSE
             READ(17, FMT='(I2)', ADVANCE='NO', IOSTAT=iostatus) temp
          END IF

          IF (iostatus < 0) THEN
             WRITE (*, *) 'Ran out of data when reading ', obs_missing, '. Exiting program.'
             STOP
          ELSE IF (iostatus > 0) THEN
             WRITE (*, *) 'Could not read ', obs_missing, '. Exiting program.'
             STOP
          END IF  

          IF (temp == 1) THEN
             obs % two_dim(i, j) % is_data = .FALSE.
             obs % two_dim_missing(next, 1) = i
             obs % two_dim_missing(next, 2) = j
             next = next + 1
          END IF
       END DO
    END DO
    
    CLOSE(16)
    CLOSE(17)


    OPEN(UNIT=16, FILE=northerness_filename, STATUS='old', IOSTAT=iostatus)
    IF (iostatus /= 0) THEN
        WRITE (*, *) 'Could not open ', northerness_filename, '. Exiting program.'
        STOP
    END IF


    OPEN(UNIT=17, FILE=northerness_missing, STATUS='old', IOSTAT=iostatus)
    IF (iostatus /= 0) THEN
        WRITE (*, *) 'Could not open ', northerness_missing, '. Exiting program.'
        STOP
    END IF

    ! this is one-dimensional and on column 1 (hard-coded)
    ! the general reading algorithm is:
    ! 1. read a line of the file into an array of correct dimensions (hard-coded)    
    ! 2. read the column 1 (hard-coded) of the array into the 
    !    variable group's value field.

    ! in this case, the array is 
    ALLOCATE(datafile(1:1)) 
    ! 1 columns on a line (the Entity's lineLength)

    !    loop length from data length (hard-coded)
    DO i = 1, 21
       READ(16, FMT=*, IOSTAT=iostatus) datafile(:)
       IF (iostatus < 0) THEN
          WRITE (*, *) 'Ran out of data when reading ', northerness_filename, '. Exiting program.'
          STOP
       ELSE IF (iostatus > 0) THEN
          WRITE (*, *) 'Could not read ', northerness_filename, '. Exiting program.'
          STOP
       END IF
       
       ! data column (1) hard-coded
       northerness % one_dim(i) % value = NINT(datafile(1))

       READ(17, FMT='(I2)', IOSTAT=iostatus) temp
       
       IF (iostatus < 0) THEN
          WRITE (*, *) 'Ran out of data when reading ', northerness_missing, '. Exiting program.'
          STOP
       ELSE IF (iostatus > 0) THEN
          WRITE (*, *) 'Could not read ', northerness_missing, '. Exiting program.'
          STOP
       END IF
       
       IF (temp == 1) THEN
          northerness % one_dim(i) % is_data = .FALSE.
          northerness % one_dim_missing(next) = i
          next = next + 1
       END IF
    END DO
    
    CLOSE(16)
    CLOSE(17)
   
  END SUBROUTINE read_data


  SUBROUTINE set_initial_values(alpha, northerness, p, x, obs)
    ! This subroutine reads the initial values of the parameters and sets them.
    IMPLICIT NONE

    TYPE(variables_real), INTENT(INOUT) :: alpha
    TYPE(variables_int), INTENT(INOUT) :: northerness      
    TYPE(variables_real), INTENT(INOUT) :: p
    TYPE(variables_int), INTENT(INOUT) :: x
    TYPE(variables_int), INTENT(INOUT) :: obs    
    
    CHARACTER(LEN=30) :: line, var_name
    ! HARD-CODED: the filename of the initial values file
    CHARACTER(LEN=*), PARAMETER :: initial_values_filename='initial_values.txt'
    INTEGER :: i, j, linenum
    INTEGER :: vertical_start, vertical_end, horizontal_start, horizontal_end
    INTEGER :: ind, iostatus

    ! The format of the file is specified in the Design document 
    ! (section Simulation parameters input files)

    OPEN(UNIT=16, FILE=initial_values_filename, STATUS='old', IOSTAT=iostatus)
    IF (iostatus /= 0) THEN
        WRITE (*, *) 'Could not open ', initial_values_filename, '. Exiting program.'
        STOP
    END IF

    !WRITE(*, *) 'Reading ', initial_values_filename

    READ (16, FMT='(A)', IOSTAT=iostatus) line

    IF (iostatus < 0) THEN
        WRITE (*, *) 'Ran out of data when reading ', initial_values_filename, '. Exiting program.'
        STOP
    ELSE IF (iostatus > 0) THEN
        WRITE (*, *) 'Could not read ', initial_values_filename, '. Exiting program.'
        STOP
    END IF  


    linenum = 1

    IF (.NOT. line == '?? initial values') THEN
       WRITE(*, *) 'Expected initial values legend in ' , initial_values_filename
       WRITE(*, *) 'on line 1 not found. Quitting.'
       CLOSE(1)
       STOP
    END IF
    
    DO WHILE (.TRUE.)       
       READ (16, FMT='(A)', IOSTAT=iostatus) line
       
       linenum = linenum + 1
       
       IF (iostatus < 0) THEN
           WRITE (*, *) 'Debug: reading file ', initial_values_filename, 'finished.'
           RETURN
       ELSE IF (iostatus > 0) THEN
           WRITE (*, *) 'Could not read ', initial_values_filename, &
               ', line ', linenum, '. Exiting program.'
           STOP
       END IF  

       
       line = TRIM(line)
   !    PRINT(*, *) line
       
       IF (LEN_TRIM(line) == 0) CYCLE
       IF (line(1:1) == '#') CYCLE
       
       ! now we have a "? alpha 1:something" at hand
       IF (line(1:1) /= '?')  THEN
          WRITE(*, *) 'Expected legend on line ', linenum, &
               ' not found, quitting.'
          STOP
       END IF
       line = line (3:LEN_TRIM(line))
       ind = INDEX(line, ' ')
       IF (ind <= 1) THEN 
          var_name = line
       ELSE
          var_name = TRIM(line(1:ind))
       END IF
!       WRITE(*, *) "variable name: ", var_name
       line = line(LEN_TRIM(var_name)+2:LEN_TRIM(line)) 
       ! trim the variable name out 
       ! and the trailing space
       
       vertical_start  = 0
       vertical_end    = 0
       horizontal_start = 0
       horizontal_end   = 0
       
       SELECT CASE (var_name)
          ! HARD-CODED: the name of the variable
       CASE ('alpha')
          
          WRITE (*, *) 'Debug: Reading initial values for alpha.'
          
          ! HARD-CODED: alpha is global so no indices are present
          ! and no loops are needed
          
          ! WRITE(*, *) 'Vertical indices are: ', vertical_start, vertical_end
             
          READ(16, *, IOSTAT=iostatus) alpha % one_dim (1) % value
          linenum = linenum + 1
          
          IF (iostatus < 0) THEN
              WRITE (*, *) 'Ran out of data when reading ', initial_values_filename, '. Exiting program.'
              STOP
          ELSE IF (iostatus > 0) THEN
              WRITE (*, *) 'Could not read ', initial_values_filename, &
                  ', line ', linenum, '. Exiting program.'
              STOP
          END IF  
             
          alpha % one_dim(1) % initial_value_set = .TRUE.
          ! end of case alpha
          
       ! HARD-CODED: the name of the variable
       CASE ('x')
          
          WRITE (*, *) 'Debug: Reading initial values for x.'
          ! HARD-CODED: x is two-dimensional so we consider the 
          ! format "x 1:20 5:400" for example
          DO i = 1, 2
             ind = INDEX(line, ':')
             line(ind:ind) = ' '
          END DO

          READ (line, FMT=*, IOSTAT=iostatus) vertical_start, vertical_end, horizontal_start, horizontal_end 

          IF (iostatus < 0) THEN
              WRITE (*, *) 'Ran out of data when reading line ', linenum, &
                  ' in ', initial_values_filename, '. Exiting program.'
              STOP
          ELSE IF (iostatus > 0) THEN
              WRITE (*, *) 'Could not read ', initial_values_filename, &
                  ', line ', linenum, '. Exiting program.'
              STOP
          END IF  

          DO i = vertical_start, vertical_end
             READ(16, *, IOSTAT=iostatus) x % two_dim (i, horizontal_start:horizontal_end) % value
             linenum = linenum +1
             
             IF (iostatus < 0) THEN
                 WRITE (*, *) 'Ran out of data when reading ', initial_values_filename, '. Exiting program.'
                 STOP
             ELSE IF (iostatus > 0) THEN
                 WRITE (*, *) 'Could not read ', initial_values_filename, &
                     ', line ', linenum, '. Exiting program.'
                 STOP
             END IF  
             
             x % two_dim(i, horizontal_start:horizontal_end) % initial_value_set = .TRUE.
             
          END DO
          ! end of case x

       CASE DEFAULT
       END SELECT
       
    END DO ! Ends the file-reading loop.
    CLOSE(16)

    ! test printout
    !    DO i = vertical_start, vertical_end
    !       DO j = horizontal_start, horizontal_end
    !          WRITE (*, *) x % two_dim (i, j) % value
    !       END DO
    !    END DO

  END SUBROUTINE set_initial_values

END MODULE input
