MODULE init_unstructured_mod
  USE prec, ONLY : rstd
  USE math_const, ONLY : radian_to_degree
  USE mpipara, ONLY : is_mpi_master
  USE data_unstructured_mod
  USE geometry, ONLY : swap_geometry
  IMPLICIT NONE
  SAVE
  PRIVATE

  ! these buffers are used when reading from the grid file
  REAL(8), ALLOCATABLE :: Ddata_read1(:),Ddata_read2(:,:),Ddata_read3(:,:,:)
  INTEGER, ALLOCATABLE :: Idata_read1(:),Idata_read2(:,:),Idata_read3(:,:,:)

  CHARACTER(LEN=*),PARAMETER :: meshfile="input/mesh_information.nc"
  INTEGER :: id_nc ! NetCDF id of mesh file open by open_local_mesh_file

  PUBLIC :: open_local_mesh_file, read_local_mesh

CONTAINS

  SUBROUTINE free_data_read()
    IF(ALLOCATED(Idata_read1)) DEALLOCATE(Idata_read1) 
    IF(ALLOCATED(Idata_read2)) DEALLOCATE(Idata_read2) 
    IF(ALLOCATED(Idata_read3)) DEALLOCATE(Idata_read3) 
    IF(ALLOCATED(Ddata_read1)) DEALLOCATE(Ddata_read1) 
    IF(ALLOCATED(Ddata_read2)) DEALLOCATE(Ddata_read2) 
    IF(ALLOCATED(Ddata_read3)) DEALLOCATE(Ddata_read3) 
  END SUBROUTINE free_data_read

  SUBROUTINE read_from_gridfile(id_nc, data_type, name)
    use netcdf_mod
    CHARACTER(*)    :: data_type, name
    INTEGER :: id_nc, id_var, status
    INTEGER :: dim_1, dim_2, dim_3
    INTEGER :: numDims, dimIDs(3) !max_var_dims

    !-------------------Reading variable-------------------------------
    status = nf90_inq_varid(id_nc, name,id_var)
    IF(status /= 0) THEN
        print *, "Not able to read variable from gridfile :", name
        STOP "Exit"
    ENDIF
    !inquire dimension
    status = nf90_inquire_variable(id_nc,id_var,dimids=dimIDs,ndims=numDims)
    status = nf90_inquire_dimension(id_nc, dimIDs(1), len = dim_1)
    status = nf90_inquire_dimension(id_nc, dimIDs(2), len = dim_2)
    status = nf90_inquire_dimension(id_nc, dimIDs(3), len = dim_3)
    SELECT CASE(numDims)
    CASE(3)
    print *,"Size of array, ",name,":", dim_1,dim_2,dim_3
    CASE(2)
    print *,"Size of array, ",name,":", dim_1,dim_2
    CASE DEFAULT
    print *,"Size of array, ",name,":", dim_1
    END SELECT

    CALL free_data_read
    SELECT CASE(data_type)
    CASE('integer')
       SELECT CASE(numDims)
       CASE(3)
          allocate(Idata_read3(dim_1,dim_2,dim_3))
          status = nf90_get_var(id_nc, id_var,Idata_read3)
          print *,"First value of array, ",name,":", Idata_read3(1,1,1)
       CASE(2)
          allocate(Idata_read2(dim_1,dim_2))
          status = nf90_get_var(id_nc, id_var,Idata_read2)
          print *,"First value of array, ",name,":", Idata_read2(1,1)
       CASE DEFAULT
          allocate(Idata_read1(dim_1))
          status = nf90_get_var(id_nc, id_var,Idata_read1)
          print *,"First value of array, ",name,":", Idata_read1(1)
       END SELECT
    CASE DEFAULT
       SELECT CASE(numDims)
       CASE(3)
          allocate(Ddata_read3(dim_1,dim_2,dim_3))
          status = nf90_get_var(id_nc, id_var,Ddata_read3)
          print *,"First value of array, ",name,":", Ddata_read3(1,1,1)
       CASE(2)
          allocate(Ddata_read2(dim_1,dim_2))
          status = nf90_get_var(id_nc, id_var,Ddata_read2)
          print *,"First value of array, ",name,":", Ddata_read2(1,1)
       CASE DEFAULT
          allocate(Ddata_read1(dim_1))
          status = nf90_get_var(id_nc, id_var,Ddata_read1)
          print *,"First value of array, ",name,":", Ddata_read1(1)
       END SELECT
    END SELECT
    
    IF(status /= nf90_NoErr) THEN
        print *, "Error when reading from grid file : ", name
        STOP "Exit"
    ENDIF

  END SUBROUTINE read_from_gridfile

  SUBROUTINE open_local_mesh_file
    USE netcdf_mod
    CHARACTER(LEN= 80) :: description
    INTEGER :: ierr, status, descriptionLength

    PRINT *,"------------------ READING FILE " , meshfile, "----------------------- "
    !open and read the input file
    ierr = nf90_open(meshfile, NF90_NOWRITE, id_nc)
    if (ierr /= nf90_noerr) then
      print *, trim(nf90_strerror(ierr))
      stop "Error reading file"
    end if

    status = nf90_inquire_attribute(id_nc, nf90_global, "description", len=descriptionLength)
    IF(status /= 0 .or. len(description) < descriptionLength) THEN
        print *, "Not enough space to put NetCDF attribute values."
        STOP "Error reading file"
    ENDIF

    !-------------------Reading global attributes-----------------------
    status = nf90_get_att(id_nc, nf90_global, "description", description)
    print *,"Data file description :",description

    CALL read_from_gridfile(id_nc, 'integer', 'primal_deg')
    primal_num = SIZE(Idata_read1)
    CALL read_from_gridfile(id_nc, 'integer', 'dual_deg')
    dual_num = SIZE(Idata_read1)
    CALL read_from_gridfile(id_nc, 'integer', 'trisk_deg')
    edge_num = SIZE(Idata_read1)
  END SUBROUTINE open_local_mesh_file


  SUBROUTINE read_field(id_nc, field)
    USE field_mod, ONLY : t_field, type_integer, type_real
    INTEGER :: id_nc
    TYPE(t_field), POINTER :: field(:), fld
    fld=>field(1)
    SELECT CASE(fld%data_type)
    CASE(type_integer)
       CALL read_from_gridfile(id_nc, 'integer', TRIM(fld%name))
       SELECT CASE(field(1)%ndim)
       CASE(2)
          fld%ival2d = Idata_read1
       CASE(3)
          fld%ival3d = Idata_read2
       END SELECT
    CASE(type_real)
       CALL read_from_gridfile(id_nc, 'float', TRIM(fld%name))
       SELECT CASE(field(1)%ndim)
       CASE(2)
          fld%rval2d = Ddata_read1
       CASE(3)
          fld%rval3d = Ddata_read2
       END SELECT
    END SELECT
  END SUBROUTINE read_field


  SUBROUTINE setup_cellset(set, lon, lat)
    USE domain_mod, ONLY : t_cellset
    TYPE(t_cellset) :: set
    REAL(rstd) :: lon(:), lat(:)
    INTEGER :: ij
    set%ncell = SIZE(lon)
    ALLOCATE(set%lon, SOURCE=lon*radian_to_degree)
    ALLOCATE(set%lat, SOURCE=lat*radian_to_degree)
    ALLOCATE(set%ij, SOURCE = (/(ij,ij=1,set%ncell)/))
    ALLOCATE(set%ind_glo, SOURCE = set%ij(:)-1)
  END SUBROUTINE setup_cellset

  SUBROUTINE read_local_mesh_bounds(dom)
    USE domain_mod, ONLY : t_domain
    USE geometry, ONLY : lon_i, lat_i, lon_e, lat_e
    TYPE(t_domain), INTENT(IN) :: dom
    REAL(rstd) :: bounds_e(2,edge_num)
    REAL(rstd), ALLOCATABLE :: lon_v(:), lat_v(:)
    INTEGER :: ij, vup, vdown

    ! read dual cell positions, not stored in geom
    CALL read_from_gridfile(id_nc, 'float', 'lon_v')
    ALLOCATE(lon_v, SOURCE=Ddata_read1)
    CALL read_from_gridfile(id_nc, 'float', 'lat_v')
    ALLOCATE(lat_v, SOURCE=Ddata_read1)

    CALL setup_cellset(dom%primal_own, lon_i, lat_i)
    CALL setup_cellset(dom%dual_own, lon_v, lat_v)
    CALL setup_cellset(dom%edge_own, lon_e, lat_e)

    PRINT *, 'read_local_mesh :', primal_num, SHAPE(dom%primal_own%ij)

    ! read primal cell bounds
    CALL read_from_gridfile(id_nc, 'float', 'primal_bounds_lon')
    ALLOCATE(dom%primal_own%bnds_lon,     source=Ddata_read2*radian_to_degree)
    CALL read_from_gridfile(id_nc, 'float', 'primal_bounds_lat')
    ALLOCATE(dom%primal_own%bnds_lat,     source=Ddata_read2*radian_to_degree)
    ! read dual cell bounds
    CALL read_from_gridfile(id_nc, 'float', 'dual_bounds_lon')
    ALLOCATE(dom%dual_own%bnds_lon,     source=Ddata_read2*radian_to_degree)
    CALL read_from_gridfile(id_nc, 'float', 'dual_bounds_lat')
    ALLOCATE(dom%dual_own%bnds_lat,     source=Ddata_read2*radian_to_degree)

    CALL collect_lonlat_e(lat_v)
    ALLOCATE(dom%edge_own%bnds_lat, SOURCE=bounds_e)
    CALL collect_lonlat_e(lon_v)
    ALLOCATE(dom%edge_own%bnds_lon, SOURCE=bounds_e)

    CONTAINS
      SUBROUTINE collect_lonlat_e(lonlat)
        REAL(rstd) :: lonlat(dual_num)
        DO ij=1, edge_num
           vdown = down(ij)
           vup = up(ij)
           bounds_e(1,ij) = lonlat(vdown)*radian_to_degree
           bounds_e(2,ij) = lonlat(vup)*radian_to_degree
        END DO
      END SUBROUTINE collect_lonlat_e

  END SUBROUTINE read_local_mesh_bounds

  SUBROUTINE read_local_mesh
    USE field_mod
    USE domain_mod, ONLY : swap_needed, domain, domain_glo
    USE geometry, ONLY : geom, lon_i, lat_i, lon_e, lat_e, ep_e
    USE netcdf_mod, ONLY : nf90_close
    IMPLICIT NONE
    INTEGER :: ij, ierr
    INTEGER, ALLOCATABLE :: cell_ij(:)
    REAL(rstd), ALLOCATABLE :: angle_e(:)
    REAL(rstd) :: clon, slon, clat, slat, & ! COS/SIN of lon/lat
         elon(3), elat(3) ! lon/lat unit vectors
    
    ! if possible data is read into fields allocated in geom(1) using read_field
    CALL read_field(id_nc, geom%Ai)
    CALL read_field(id_nc, geom%Av)
    CALL read_field(id_nc, geom%fv)
    CALL read_field(id_nc, geom%le_de)
    CALL read_field(id_nc, geom%le)
    CALL read_field(id_nc, geom%de)
    CALL read_field(id_nc, geom%lon_i)
    CALL read_field(id_nc, geom%lat_i)
    CALL read_field(id_nc, geom%lon_e)
    CALL read_field(id_nc, geom%lat_e)

    swap_needed = .TRUE.
    CALL swap_geometry(1)
    
    PRINT *, 'read_local_mesh : primal_num =', primal_num, domain_glo(1)%primal_own%ncell

    ! other fields are defined only in the case of an unstructured mesh
    ! those fields are read into simple arrays
    CALL read_from_gridfile(id_nc, 'integer', 'primal_deg')
    ALLOCATE(primal_deg, source = Idata_read1)
    CALL read_from_gridfile(id_nc, 'integer', 'primal_edge')
    ALLOCATE(primal_edge, source = Idata_read2)
    CALL read_from_gridfile(id_nc, 'integer', 'primal_ne')
    ALLOCATE(primal_ne, source = Idata_read2)
    CALL read_from_gridfile(id_nc, 'integer', 'primal_vertex') 
    ALLOCATE(primal_vertex, source = Idata_read2)

    CALL read_from_gridfile(id_nc, 'integer', 'dual_deg')
    ALLOCATE(dual_deg, source = Idata_read1)
    CALL read_from_gridfile(id_nc, 'integer', 'dual_edge')
    ALLOCATE(dual_edge, source = Idata_read2)
    CALL read_from_gridfile(id_nc, 'integer', 'dual_ne')
    ALLOCATE(dual_ne, source = Idata_read2)
    CALL read_from_gridfile(id_nc, 'integer', 'dual_vertex') 
    ALLOCATE(dual_vertex, source = Idata_read2)

    CALL read_from_gridfile(id_nc, 'integer', 'left')
    ALLOCATE(left, source = Idata_read1)
    CALL read_from_gridfile(id_nc, 'integer', 'right')
    ALLOCATE(right, source = Idata_read1)
    CALL read_from_gridfile(id_nc, 'integer', 'up')
    ALLOCATE(up, source = Idata_read1)
    CALL read_from_gridfile(id_nc, 'integer', 'down')
    ALLOCATE(down, source = Idata_read1)
    CALL read_from_gridfile(id_nc, 'float', 'angle_e')
    ALLOCATE(angle_e, source = Ddata_read1)
    CALL read_from_gridfile(id_nc, 'integer', 'trisk_deg')
    ALLOCATE(trisk_deg, source = Idata_read1)
    CALL read_from_gridfile(id_nc, 'integer', 'trisk')
    ALLOCATE(trisk, source = Idata_read2)
    CALL read_from_gridfile(id_nc, 'float', 'Riv2')
    ALLOCATE(Riv2, source = Ddata_read2)
    CALL read_from_gridfile(id_nc, 'float', 'wee')
    ALLOCATE(wee(SIZE(Ddata_read2,1), SIZE(Ddata_read2,2), 1))
    wee(:,:,1) = Ddata_read2(:,:)

    ! read cell centers and bounds
    CALL read_local_mesh_bounds(domain(1))
    CALL read_local_mesh_bounds(domain_glo(1))

    CALL free_data_read ! free buffers after reading all data from grid file
    ierr = nf90_close(id_nc)

    ! now process some of the data we just read

    max_primal_deg = SIZE(primal_edge,1)
    max_dual_deg = SIZE(dual_edge,1)
    max_trisk_deg = SIZE(trisk,1)

    DO ij=1,edge_num
       clon = COS(lon_e(ij))
       slon = SIN(lon_e(ij))
       clat = COS(lat_e(ij))
       slat = SIN(lat_e(ij))
       ! x,y,z = clat*clon, clat*slon, slat
       elon = (/ -slon, clon, 0. /)
       elat = (/ -slat*clon, -slat*slon, clat /)
       ep_e(:,ij) = COS(angle_e(ij))*elon + SIN(angle_e(ij))*elat
    END DO
    
    DEALLOCATE(angle_e)

    CALL swap_geometry(1)
    swap_needed = .FALSE.

  END SUBROUTINE read_local_mesh

END MODULE init_unstructured_mod