!
!     (c) 2019 Guide Star Engineering, LLC
!     This Software was developed for the US Nuclear Regulatory Commission (US NRC)
!     under contract "Multi-Dimensional Physics Implementation into Fuel Analysis under
!     Steady-state and Transients (FAST)", contract # NRC-HQ-60-17-C-0007
!
SUBMODULE (class_exodus) class_exodus_procedures

    USE iso_fortran_env, ONLY : INT64, REAL64
    USE class_psblas
    USE class_cell
    USE class_face
    USE class_vertex
    USE class_mesh
    USE class_connectivity
    USE class_iterating
    USE class_scalar_field
    USE class_vector_field
    USE class_vector
    USE class_output

    IMPLICIT NONE
    INCLUDE 'exodusII.inc'

    !! author: Hari Radhakrishnan, GSE
    !! date 12/12/2019
    !!
    !! This module implements the routines for writing an exodus file with mesh and results
    CONTAINS

        MODULE PROCEDURE write_exo_morfeus

            CHARACTER(len=100) :: filename
            INTEGER, DIMENSION(:), ALLOCATABLE :: cell_ids
            INTEGER, DIMENSION(:), ALLOCATABLE :: v2cconn
            INTEGER, DIMENSION(:), ALLOCATABLE :: icverts
            INTEGER, DIMENSION(:), ALLOCATABLE :: iproc
            REAL(psb_dpk_), DIMENSION(:,:), ALLOCATABLE :: points

            INTEGER(kind = int64) :: ierr, titlelen
            INTEGER(kind = int64) :: exodus_file_id, num_dims, num_nodes, num_elems,  &
              num_elem_blks, num_node_sets, num_side_sets
            INTEGER(kind = int64) :: cpu_ws, io_ws, num_props, prop_value
            REAL :: vers
            REAL(kind = real64), ALLOCATABLE :: xx(:), yy(:), zz(:)
            INTEGER :: info, err_act
            INTEGER :: icontxt, mypnum
            INTEGER :: i, ig, ic, nc2g, nv2c, nconn, ncells, ngc, ngroups
            INTEGER(kind = int64), ALLOCATABLE :: conn(:)
            CHARACTER(len = MXSTLN) :: coord_names(3), ctype
            INTEGER, POINTER :: ic2g(:) => NULL()
            INTEGER, POINTER :: iv2c(:) => NULL()
            INTEGER, ALLOCATABLE  :: igroup(:)
            INTEGER, ALLOCATABLE :: i_loc(:)
            TYPE(cell), ALLOCATABLE :: cells(:)
            TYPE(face), ALLOCATABLE :: faces(:)
            TYPE(vertex), ALLOCATABLE :: verts(:)
            TYPE(connectivity) :: v2f, v2c, f2c
            CHARACTER(len=32) :: path
            REAL(psb_dpk_), ALLOCATABLE :: scalar_local(:), scalar_global(:), f_s(:)
            TYPE(vector), ALLOCATABLE :: vector_local(:), vector_global(:)
            REAL(psb_dpk_), ALLOCATABLE :: f_x(:), f_y(:), f_z(:)
            ! Sets error handling for PSBLAS-2 routines
            CALL psb_erractionsave(err_act)
            mypnum  = mypnum_()
            icontxt = icontxt_()

            ! Sets output path
            ! IF(PRESENT(iter)) CALL out%set_output_path(iter)
            ! path = out%path_()
            filename = 'out_nemo.e' !TRIM(path)
            print *, filename
            ! Global number of cells
            ncells = psb_cd_get_global_cols(msh%desc_c)

            CALL psb_geall(i_loc,msh%desc_c,info)
            CALL psb_check_error(info,'write_mesh','psb_geall',icontxt)

            ALLOCATE(igroup(ncells),iproc(ncells),stat=info)
            IF(info /= 0) THEN
                WRITE(*,100)
                CALL abort_psblas
            END IF

            ! Gathers MESH components on P0
            CALL l2g_vertex(msh%verts,verts,msh%desc_v)
            CALL l2g_face(msh%faces,faces,msh%desc_f,msh%desc_c)
            CALL l2g_cell(msh%cells,cells,msh%desc_c)
            CALL l2g_conn(msh%v2f,v2f,msh%desc_v,msh%desc_f)
            CALL l2g_conn(msh%v2c,v2c,msh%desc_v,msh%desc_c)
            CALL l2g_conn(msh%f2c,f2c,msh%desc_f,msh%desc_c)

            ! Gathers processor IDs
            i_loc(:) = mypnum
            CALL psb_gather(iproc,i_loc,msh%desc_c,info,root=0)
            CALL psb_check_error(info,'write_mesh','psb_gather',icontxt)

            ! Gathers group IDs
            i_loc(:) = 0
            ngroups = msh%c2g%nel_()

            DO ig = 1, ngroups
                CALL msh%c2g%get_ith_conn(ic2g,ig)
                ngc = SIZE(ic2g) ! number of group cells
                DO i = 1, ngc
                    ic = ic2g(i)
                    i_loc(ic) = ig
                END DO
            END DO

            CALL psb_gather(igroup,i_loc,msh%desc_c,info,root=0)
            CALL psb_check_error(info,'write_mesh','psb_gather',icontxt)

            ! Open exodusII file for writing
            cpu_ws = 0
            io_ws = 8
            call exopts (EXVRBS, IERR) !Verbose error reporting of EXODUS routines
            exodus_file_id = excre("output.exo", EXCLOB, cpu_ws, io_ws, ierr)

            num_dims = msh%ncd
            num_elem_blks = ngroups
            num_node_sets = 0
            num_side_sets = msh%nbc
            num_nodes = size(verts)
            num_elems = size(cells)

            ! Write the database parameters to the exodus file
            CALL expini(exodus_file_id, "MORFEUS Output", num_dims, num_nodes, &
              num_elems, num_elem_blks, num_node_sets, num_side_sets, ierr)

            ! Write the mesh coordinates
            ALLOCATE(xx(num_nodes), yy(num_nodes), zz(num_nodes))

            xx(:) = verts(:)%x_()
            yy(:) = verts(:)%y_()
            zz(:) = verts(:)%z_()

            CALL expcor(exodus_file_id, xx, yy, zz, ierr)

            ! Write the mesh coordinate names
            coord_names(1) = "X"
            coord_names(2) = "Y"
            coord_names(3) = "Z"
            CALL expcon(exodus_file_id, coord_names, ierr)

            ! Write the element block parameters
            DO ig = 1, num_elem_blks !Loop over groups
              CALL msh%c2g%get_ith_conn(ic2g,ig)
              nc2g = SIZE(ic2g) !Number of cells in group
              CALL msh%v2c%get_ith_conn(iv2c,1) ! First cell
              nv2c = SIZE(iv2c) ! Number of nodes per cell
              ctype = "HEX8"
              CALL expelb(exodus_file_id, INT(ig,INT64), ctype, INT(nc2g,INT64), INT(nv2c,INT64), 1_INT64, ierr)

              nconn = nv2c*nc2g !Connectivity vector size

              print *, ig, nconn
              IF (ALLOCATED(conn)) DEALLOCATE(conn)
              ALLOCATE(conn(nconn))
              i = 0
              DO ic = 1, nc2g
                CALL msh%v2c%get_ith_conn(iv2c,ic)
                !print *, ic, i, iv2c
                conn(i+1:i+nv2c) = iv2c
                i = i + nv2c
              END DO
              call expelc(exodus_file_id, INT(ig,INT64), conn, ierr)
            END DO

            ! Loop over scalar fields and gather them on the root processor
            DO i = 1, size(sfield)
              ! Is the field cell-centered?
              IF(sfield(i)%on_faces_()) THEN
                  WRITE(*,100)
                  CYCLE
              END IF
              CALL sfield(i)%get_x(scalar_local)


              ALLOCATE(scalar_global(ncells),stat=info)
              IF(info /= 0) THEN
                  WRITE(*,200)
                  CALL abort_psblas
              END IF

              ! Gathers cell-centered values
              CALL psb_gather(scalar_global,scalar_local,msh%desc_c,info,root=0)

              IF(mypnum == 0) THEN
                DO ig = 1, num_elem_blks
                  CALL msh%c2g%get_ith_conn(ic2g,ig)
                  nc2g = SIZE(ic2g) !Number of cells in group
                  IF(ALLOCATED(f_s)) DEALLOCATE(f_s)
                  ALLOCATE(f_s(nc2g))
                  f_s(1:nc2g) = scalar_global(ic2g(1:nc2g))
                  CALL expev(exodus_file_id, 0.0, i, INT(ig, INT64), INT(nc2g, INT64), f_s, ierr)
                END DO
              END IF
            END DO

            ! Loop over vector fields and gather them on the root processor
            DO i = 1, size(vfield)
              ! Is the field cell-centered?
              IF(vfield(i)%on_faces_()) THEN
                  WRITE(*,100)
                  CYCLE
              END IF
              CALL vfield(i)%get_x(vector_local)


              ALLOCATE(vector_global(ncells),stat=info)
              IF(info /= 0) THEN
                  WRITE(*,200)
                  CALL abort_psblas
              END IF

              ! Gathers cell-centered values
              CALL l2g_vector(vector_global,vector_local,msh%desc_c)

              IF(mypnum == 0) THEN
                DO ig = 1, num_elem_blks
                  CALL msh%c2g%get_ith_conn(ic2g,ig)
                  nc2g = SIZE(ic2g) !Number of cells in group
                  IF(ALLOCATED(f_x)) DEALLOCATE(f_x, f_y, f_z)
                  ALLOCATE(f_x(nc2g), f_y(nc2g), f_z(nc2g))
                  f_x(1:nc2g) = vector_global(ic2g(1:nc2g))%x_()
                  f_y(1:nc2g) = vector_global(ic2g(1:nc2g))%y_()
                  f_z(1:nc2g) = vector_global(ic2g(1:nc2g))%z_()
                  CALL expev(exodus_file_id, 0.0, i, INT(ig, INT64), INT(nc2g, INT64), f_s, ierr)
                END DO
              END IF
            END DO


            ! DO i = 1, size(sfield) + 3*size(vfield)
            !   IF (i <= size(sfield)) THEN
            !     !write sfield to exodus file
            !   ELSE

            !   END IF
            ! END DO
            ! Close the exodus file
            CALL exclos(exodus_file_id, ierr)
            STOP

100         FORMAT(' WARNING! Face-centered field in WRITE_EXODUS. Field will not be output.')
200         FORMAT(' ERROR! Memory allocation failure in WRITE_EXODUS')

        END PROCEDURE write_exo_morfeus

END SUBMODULE class_exodus_procedures
