yowpdlibmain.F90

Go to the documentation of this file.

!PDLIB Software License
!
!Software, as understood herein, shall be broadly interpreted as being inclusive of algorithms,
!source code, object code, data bases and related documentation, all of which shall be furnished
!free of charge to the Licensee. Corrections, upgrades or enhancements may be furnished and, if
!furnished, shall also be furnished to the Licensee without charge. NOAA, however, is not
!required to develop or furnish such corrections, upgrades or enhancements.
!Roland & Partner software, whether that initially furnished or corrections or upgrades,
!are furnished "as is". Roland & Partner furnishes its software without any warranty
!whatsoever and is not responsible for any direct, indirect or consequential damages
!that may be incurred by the Licensee. Warranties of merchantability, fitness for any
!particular purpose, title, and non-infringement, are specifically negated.
!The Licensee is not required to develop any software related to the licensed software.
!However, in the event that the Licensee does so, the Licensee is required to offer same
!to Roland & Partner for inclusion under the instant licensing terms with Roland & Partner
!licensed software along with documentation regarding its principles, use and its advantages.
!This includes changes to the wave model proper including numerical and physical approaches
!to wave modeling, and boundary layer parameterizations embedded in the wave model
!A Licensee may reproduce sufficient software to satisfy its needs.
!All copies shall bear the name of the software with any version number
!as well as replicas of any applied copyright notice, trademark notice,
!other notices and credit lines. Additionally, if the copies have been modified,
!e.g. with deletions or additions, this shall be so stated and identified.
!All of Licensee's employees who have a need to use the software may have access
!to the software but only after reading the instant license and stating, in writing,
!that they have read and understood the license and have agreed to its terms.
!Licensee is responsible for employing reasonable efforts to assure
!that only those of its employees that should have access to the software, in fact, have access.
!The Licensee may use the software for any purpose relating to sea state prediction.
!No disclosure of any portion of the software, whether by means of a media or verbally,
!may be made to any third party by the Licensee or the Licensee's employees
!The Licensee is responsible for compliance with any applicable export or
!import control laws of the United States, the European Union and Germany.
!
!© 2009 Roland&Partner, Georgenstr.32, 64297 Germany. All rights reserved.
!PDLIB is a trademark of Roland & Partner. No unauthorized use without permission.
!
module yowpdlibmain
  use yowerr
  use yowdatapool, only : rkind
  use w3servmd,    only : print_memcheck
 
  !module default
  implicit none
 
  private
  public :: initfromgriddim, finalizepd
 
contains
 
  subroutine initfromgriddim(MNP, MNE, INE_global, secDim, MPIcomm)
    use yowdatapool,       only: myrank, debugprepartition, debugpostpartition
    use yownodepool,       only: np_global, np, np_perprocsum, ng, ipgl, iplg, npa
    use yowelementpool,    only: ne_global,ne
    use yowsidepool,       only: ns, ns_global
    use yowexchangemodule, only: nconndomains, setdimsize
    use yowrankmodule,     only: initrankmodule, ipgl_npa
 
    integer, intent(in) :: mnp, mne
    integer, intent(in) :: ine_global(3,mne)
    integer, intent(in) :: secdim
    integer, intent(in) :: mpicomm
    integer :: istat, memunit
 
    ! note: myrank=0 until after initMPI is called, so only rank=0 file
    ! contains the 'section 1' information
    memunit = 70000+myrank
 
    call setdimsize(secdim)
    call print_memcheck(memunit, 'memcheck_____:'//' WW3_PDLIB SECTION 1')
#ifdef W3_DEBUGINIT
    print *, '1: MPIcomm=', mpicomm
#endif
    call initmpi(mpicomm)
 
    memunit = 70000+myrank
    call print_memcheck(memunit, 'memcheck_____:'//' WW3_PDLIB SECTION 2')
#ifdef W3_DEBUGINIT
    print *, '2: After initMPI'
#endif
    call assignmesh(mnp, mne)
    call print_memcheck(memunit, 'memcheck_____:'//' WW3_PDLIB SECTION 3')
#ifdef W3_DEBUGINIT
    print *, '3: After assignMesh'
#endif
    call prepartition()
    call print_memcheck(memunit, 'memcheck_____:'//' WW3_PDLIB SECTION 4')
#ifdef W3_DEBUGINIT
    print *, '3: After prePartition'
#endif
    call findconnnodes(ine_global)
    call print_memcheck(memunit, 'memcheck_____:'//' WW3_PDLIB SECTION 5')
#ifdef W3_DEBUGINIT
    print *, '4: After findConnNodes'
#endif
    if(debugprepartition) then
      if(myrank == 0) then
        write(*,*) "pre-partition"
        write(*,*) "# Nodes: ", np_global
        write(*,*) "# Elements: ", ne_global
        write(*,*) "# Sides ", ns_global
        write(*,*) "np_perProcSum :", np_perprocsum
      end if
      write(*,*) "Thread", myrank, "# local nodes ", np
      write(*,*) "Thread", myrank, "# local sides" , ns
    endif
 
    !   call writeMesh()
#ifdef W3_DEBUGINIT
    print *, '4.1: After findConnNodes'
#endif
    !    CALL REAL_MPI_BARRIER_PDLIB(MPIcomm, "Before call to runParmetis")
    call runparmetis(mnp)
    call print_memcheck(memunit, 'memcheck_____:'//' WW3_PDLIB SECTION 6')
    !    CALL REAL_MPI_BARRIER_PDLIB(MPIcomm, "After call to runParmetis")
#ifdef W3_DEBUGINIT
    print *, '5: After runParmetis'
#endif
    call postpartition
    call print_memcheck(memunit, 'memcheck_____:'//' WW3_PDLIB SECTION 7')
#ifdef W3_DEBUGINIT
    print *, 'Before findGhostNodes'
#endif
    call findghostnodes
    call print_memcheck(memunit, 'memcheck_____:'//' WW3_PDLIB SECTION 8')
 
    call findconndomains
    call print_memcheck(memunit, 'memcheck_____:'//' WW3_PDLIB SECTION 9')
 
    call exchangeghostids
    call print_memcheck(memunit, 'memcheck_____:'//' WW3_PDLIB SECTION 10')
 
    call postpartition2(ine_global)
    call print_memcheck(memunit, 'memcheck_____:'//' WW3_PDLIB SECTION 11')
 
    call initrankmodule
    call print_memcheck(memunit, 'memcheck_____:'//' WW3_PDLIB SECTION 12')
 
    call computetria_ien_si_ccon
    call print_memcheck(memunit, 'memcheck_____:'//' WW3_PDLIB SECTION 13')
 
    call computeia_ja_posi_nnz
    call print_memcheck(memunit, 'memcheck_____:'//' WW3_PDLIB SECTION 14')
 
    if(debugpostpartition) then
      if(myrank == 0) then
        write(*,*) "New data after partition"
        write(*,*) "# Nodes: ", np_global
        write(*,*) "# Elements: ", ne_global
        write(*,*) "# Sides ", ns_global
        write(*,*) "np_perProcSum :", np_perprocsum
      end if
      write(*,*) "Thread", myrank, "# local elements ", ne
      write(*,*) "Thread", myrank, "# local nodes ", np
      write(*,*) "Thread", myrank, "# local sides" , ns
      write(*,*) "Thread", myrank, "# of ghosts", ng
      write(*,*) "Thread", myrank, "# of neighbor domains", nconndomains
    endif
  end subroutine initfromgriddim
 
 
 
  SUBROUTINE real_mpi_barrier_pdlib(TheComm, string)
 
    include "mpif.h"
    integer, intent(in) :: TheComm
    character(*), intent(in) :: string
    integer NbProc, eRank
    integer :: istatus(MPI_STATUS_SIZE)
    integer ierr, iField(1), iProc
    !      Print *, 'Start of REAL_MPI_BARRIER_PDLIB'
    CALL mpi_comm_rank(thecomm, erank, ierr)
    CALL mpi_comm_size(thecomm, nbproc, ierr)
    !      Print *, '   eRank=', eRank, ' NbProc=', NbProc
    ifield(1)=1
    IF (erank .eq. 0) THEN
      DO iproc=2,nbproc
        !          Print *, '   Before MPI_RECV 1 iProc=', iProc
        CALL mpi_recv(ifield, 1, mpi_integer, iproc-1, 711, thecomm, istatus, ierr)
        !          Print *, '   Before MPI_SEND 1'
        CALL mpi_send(ifield, 1, mpi_integer, iproc-1, 712, thecomm, ierr)
      END DO
    ELSE
      !        Print *, '   Before MPI_SEND 2 eRank=', eRank
      CALL mpi_send(ifield, 1, mpi_integer, 0, 711, thecomm, ierr)
      !        Print *, '   Before MPI_RECV 2 eRank=', eRank
      CALL mpi_recv(ifield, 1, mpi_integer, 0, 712, thecomm, istatus, ierr)
    END IF
    !      Print *, 'Passing barrier string=', string
  END SUBROUTINE real_mpi_barrier_pdlib
  !--------------------------------------------------------------------------
  ! Init MPI
  !--------------------------------------------------------------------------
 
  subroutine initmpi(MPIcomm)
    use yowdatapool, only: comm, ntasks, myrank
    use yowerr
    use mpi
 
    integer, intent(in) :: MPIcomm
    logical :: flag
    integer :: ierr
#ifdef W3_DEBUGINIT
    print *, '2: MPIcomm=', mpicomm
#endif
    if(mpicomm == mpi_comm_null) then
      CALL abort("A null communicator is not allowed")
    endif
 
    comm = mpicomm
    call mpi_initialized(flag, ierr)
    if(ierr/=mpi_success) call parallel_abort(error=ierr)
 
    if(flag .eqv. .false.) then
#ifdef W3_DEBUGINIT
      print *, 'Before MPI_INIT yowpdlibmain'
#endif
      call mpi_init(ierr)
#ifdef W3_DEBUGINIT
      print *, 'After MPI_INIT yowpdlibmain'
#endif
      if(ierr/=mpi_success) call parallel_abort(error=ierr)
    endif
 
    ! Get number of processors
    call mpi_comm_size(comm, ntasks,ierr)
    if(ierr/=mpi_success) call parallel_abort(error=ierr)
 
    ! Get rank
    call mpi_comm_rank(comm, myrank,ierr)
    if(ierr/=mpi_success) call parallel_abort(error=ierr)
  end subroutine initmpi
 
 
  subroutine assignmesh(MNP, MNE)
    use yownodepool,    only: nodes_global, np_global
    use yowelementpool, only: ne_global
    use yowerr,       only: parallel_abort
 
    integer, intent(in) :: MNP, MNE
    !integer, intent(in) :: INE(3,MNE)
    integer :: stat, i
 
    np_global=mnp
    if(allocated(nodes_global)) deallocate(nodes_global)
    allocate(nodes_global(np_global), stat=stat);
    if(stat/=0) CALL abort('nodes_global() allocate failure')
 
    do i =1, np_global
      nodes_global(i)%id_global = i
    end do
 
    ne_global=mne
 
    !if(allocated(INE_global)) deallocate(INE_global)
    !allocate(INE_global(3, ne_global), stat=stat);
    !if(stat/=0) CALL ABORT('INE_global allocate failure')
    !INE_global = INE
  end subroutine assignmesh
 
 
  !-------------------------------------------------------------------------
  ! pre-partition: divide the mesh into nTasks parts.
  !-------------------------------------------------------------------------
 
  subroutine prepartition
    use yowdatapool, only: ntasks ,myrank
    use yowerr,    only: parallel_abort
    use yownodepool, only: np_global, np, np_perproc, np_perprocsum, iplg
 
    integer :: i, stat
 
    ! determine equal number of nodes in each processor (except for the last one).
    ! and create a provisional node local to global mapping iplg
 
    ! start the arrays from 0, because the first thread id is 0
    if(allocated(np_perproc)) deallocate(np_perproc)
    allocate(np_perproc(0:ntasks-1), stat=stat)
    if(stat/=0) call parallel_abort('np_perProc allocation failure')
 
    if(allocated(np_perprocsum)) deallocate(np_perprocsum)
    allocate(np_perprocsum(0:ntasks), stat=stat)
    if(stat/=0) call parallel_abort('np_perProcSum allocation failure')
 
    np_perprocsum = 0
    np = np_global / ntasks
 
    do i = 0, ntasks-2
      np_perproc(i) = np
      np_perprocsum(i+1) = np_perprocsum(i) + np
    end do
 
    np_perproc(ntasks-1) = np_global - np_perprocsum(ntasks-1)
    np_perprocsum(ntasks) = np_perprocsum(ntasks-1) + np_perproc(ntasks-1)
    np = np_perproc(myrank)
 
    ! create a provisional node local to global mapping iplg
    ! this will override later with the data from parmetis
    if(allocated(iplg)) deallocate(iplg)
    allocate(iplg(np), stat=stat);
    if(stat/=0) call parallel_abort(' iplg allocate failure')
    do i = 1, np
      iplg(i) = i + np_perprocsum(myrank)
    end do
  end subroutine prepartition
 
  !-------------------------------------------------------------------------
  !  Create the connected Nodes array
  !-------------------------------------------------------------------------
 
  subroutine findconnnodes(INE_global)
    use yowerr,       only: parallel_abort
    use yownodepool,    only: np, np_global, nodes_global, nodes, maxconnnodes, t_node, connnodes_data
    use yowelementpool, only: ne_global
    use yowsidepool,    only: ns, ns_global
 
    integer, intent(in) :: INE_global(3,ne_global)
    integer :: i, j, stat
    type(t_node), pointer :: node
    integer JPREV, JNEXT
 
    ! Loop over all nlements
    ! look at their nodes
    ! get the node 1, insert node 2 and 3 into the connected nodes array
    ! do that for node 2 and 3 again
 
    ! implementation is some different
    ! loop over alle elements to get the # of connected nodes
    ! allocate space
    ! loop a second time to insert the connected nodes
 
    ! first loop
    do i = 1, ne_global
      do j = 1, 3
        node => nodes_global(ine_global(j,i))
        call node%insertConnNode()
        call node%insertConnNode()
      end do
    end do
 
    maxconnnodes = maxval(nodes_global(:)%nConnNodes)
    nodes_global(:)%nConnNodes = 0
 
    ! allocate space
    if(allocated(connnodes_data)) deallocate(connnodes_data)
    allocate(connnodes_data(np_global, maxconnnodes), stat = stat)
    if(stat/=0) call parallel_abort('connNodes allocation failure')
 
    ! second loop
    do i = 1, ne_global
      DO j=1,3
        IF (j .eq. 3) THEN
          jnext = 1
        ELSE
          jnext = j + 1
        END IF
        IF (j .eq. 1) THEN
          jprev = 3
        ELSE
          jprev = j - 1
        END IF
        node => nodes_global(ine_global(j,i))
        call node%insertConnNode(ine_global(jnext,i))
        call node%insertConnNode(ine_global(jprev,i))
      END DO
 
    end do
 
    ns = 0
    ! calc # sides local
    do i = 1, np
      node => nodes(i)
      ns = ns + node%nConnNodes
    end do
 
    do i = 1, np_global
      node => nodes_global(i)
      ns_global = ns_global + node%nConnNodes
    end do
  end subroutine findconnnodes
 
 
  !------------------------------------------------------------------------
  ! Collect all data for parmetis und partition the mesh
  !------------------------------------------------------------------------
 
  subroutine runparmetis(MNP)
    use yowerr,    only: parallel_abort
    use yowdatapool, only: debugparmetis,debugpartition, ntasks, myrank, itype, comm
    use yownodepool, only: np, npa, np_global, nodes, nodes_global, t_node, iplg, np_perprocsum, np_perproc
    use yowsidepool, only: ns
    use yowelementpool, only: ne, ne_global
    use w3gdatmd, only: xgrd, ygrd
    use mpi
 
    integer, intent(in) :: MNP
 
    ! Parmetis
    ! Node neighbor information
    integer :: wgtflag, numflag, ndims, nparts, edgecut, ncon
    integer, allocatable :: xadj(:), part(:), vwgt(:), adjwgt(:), vtxdist(:), options(:), adjncy(:), iweights(:)
    ! parmetis need single precision
    real(4), allocatable :: xyz(:), tpwgts(:), ubvec(:)
    integer :: IP_glob, itmp
    integer :: ref
    logical :: lexist = .false.
 
    ! Node to domain mapping.
    ! np_global long. give the domain number for die global node number
    integer, allocatable :: node2domain(:)
 
    ! Mics
    integer :: i, j, stat, ierr
    type(t_node), pointer :: node, nodeNeighbor
 
    !    CALL REAL_MPI_BARRIER_PDLIB(comm, "runParmetis, step 1")
    ! Create xadj and adjncy arrays. They holds the nodes neighbors in CSR Format
    ! Here, the adjacency structure of a graph is represented by two arrays,
    ! xadj[n+1] and adjncy[m]; n vertices and m edges. Every edge is listen twice
    allocate(adjncy(ns), stat=stat)
    if(stat/=0) call parallel_abort('adjncy allocation failure')
    allocate(xadj(np+1), stat=stat)
    if(stat/=0) call parallel_abort('xadj allocation failure')
    !    CALL REAL_MPI_BARRIER_PDLIB(comm, "runParmetis, step 2")
 
    xadj = 0
    xadj(1) = 1
    adjncy = 0
    do i=1,np
      node => nodes(i)
      xadj(i+1) = xadj(i) + node%nConnNodes
      if(debugparmetis) write(710+myrank,*) i, xadj(i), xadj(i+1)
      if(node%nConnNodes == 0) then
        write(*,*) "Thread", myrank,"global node has no conn nodes", node%id_global
      endif
      do j=1, node%nConnNodes
        nodeneighbor => node%connNodes(j)
        adjncy(j + xadj(i) - 1) = nodeneighbor%id_global
        if(debugparmetis) write(710+myrank,*) i, j, j + xadj(i) - 1, adjncy(j + xadj(i) - 1), nodeneighbor%id_global
      end do
    end do
    !    CALL REAL_MPI_BARRIER_PDLIB(comm, "runParmetis, step 3")
 
    ! Option for Parmetis
    allocate(options(3))
    options(1)=1   ! 0: default options; 1: user options
    if(debugparmetis) then
      options(2)=15
    else
      options(2)=0  ! Level of information returned: see defs.h in ParMETIS-Lib dir
    endif
    options(3)=15  ! Random number seed
    CALL real_mpi_barrier_pdlib(comm, "runParmetis, step 4")
    ! Fortran-style numbering that starts from 1
    numflag = 1
    ! # dimensions of the space in which the graph is embedded
    ndims = 2
    ! # parts the mesh is divide. Usually nTasks
    nparts = ntasks
    ! # weight per node
    ncon = 1
    ! wgtflag: 0: none (vwgt and adjwgt are NULL); 1: edges (vwgt is NULL); 2: vertices (adjwgt is NULL); 3: both vertices & edges;
 
#ifdef WEIGHTS
    wgtflag = 2
    INQUIRE ( file='weights.ww3', exist = lexist )
    IF (lexist) THEN
      OPEN(100001,file='weights.ww3',form='FORMATTED',status='unknown')
      allocate(iweights(np_global)); iweights = 0
      do i = 1, np_global
        read(100001,*) iweights(i)
      enddo
      CLOSE(100001)
    ELSE
      wgtflag = 0
    ENDIF
#else
    wgtflag = 0
#endif
 
    ! Create weights
    allocate(vwgt(np*ncon), stat=stat)
    if(stat/=0) call parallel_abort('vwgt allocation failure')
 
#ifdef WEIGHTS
    if (lexist) then
      do i = 1, np
        itmp = max(1,int(real((iweights(iplg(i))+100))))
        !vwgt(i) = max(1,int(real(iweights(iplg(i)))/maxval(iweights)) * 10)
        !vwgt(i) = max(1,iweights(iplg(i)))
        !vwgt(i) = max(1,itmp)
        vwgt(i) = itmp
      enddo
      vwgt = 1
      deallocate(iweights)
    else
      vwgt = 1
    endif
#else
    vwgt = 1
#endif
 
    allocate(adjwgt(ns), stat=stat)
    if(stat/=0) call parallel_abort('adjwgt allocation failure')
    adjwgt = 1
 
    ! Vertex weight fraction
    allocate(tpwgts(ncon*ntasks),stat=stat)
    if(stat/=0) call parallel_abort('partition: tpwgts allocation failure')
    tpwgts=1.0/real(ntasks)
 
    ! Imbalance tolerance
    allocate(ubvec(ncon),stat=stat)
    if(stat/=0) call parallel_abort('partition: ubvec allocation failure')
    ubvec=1.01
    CALL real_mpi_barrier_pdlib(comm, "runParmetis, step 5")
 
    ! Partition dual graph
    allocate(xyz(2*np),stat=stat)
    if(stat/=0) call parallel_abort('xyz: ubvec allocation failure')
    if(debugparmetis) write(710+myrank,*) 'np_global, ne_global, np, npa, ne'
    if(debugparmetis) write(710+myrank,*) np_global, ne_global, np, npa, ne
    do i = 1, np
      ip_glob = iplg(i)
      !AR: this is questionable ...
      xyz(2*(i-1)+1) = real(xgrd(1,ip_glob))
      xyz(2*(i-1)+2) = real(ygrd(1,ip_glob))
      if(debugparmetis) then
        write(710+myrank,*) i, np, xyz(2*(i-1)+1), xyz(2*(i-1)+2)
        call flush(710+myrank)
      endif
    end do
    CALL real_mpi_barrier_pdlib(comm, "runParmetis, step 6")
 
    ! Partition array returned from ParMeTiS
    allocate(part(np),stat=stat)
    if(stat/=0) call parallel_abort('part: ubvec allocation failure')
 
    ! ParMeTiS vertex distribution array (starts at 1)
    allocate(vtxdist(ntasks+1),stat=stat)
    if(stat/=0) call parallel_abort('partition: vtxdist allocation failure')
 
    call mpi_allgather(np_perprocsum(myrank)+1, 1, itype, vtxdist, 1, itype, comm, ierr)
    if(ierr/=mpi_success) call parallel_abort('partition: mpi_allgather',ierr)
    vtxdist(ntasks+1)=np_global+1
    !    CALL REAL_MPI_BARRIER_PDLIB(comm, "runParmetis, step 7")
 
    ! check vtxdist
    ! myrank stats from 0
    if((vtxdist(myrank+2) - vtxdist(myrank+1)) < 1) then
      write(*,*) "Thread", myrank, "has no nodes"
      write(*,*) "vtxdist", vtxdist
      CALL abort("Poor initial vertex distribution detected")
    endif
 
 
 
    !  My notes from manual:
    !  p: # of processors;
    !  n: total # of vertices (local) in graph sense;
    !  m: total # of neighboring vertices ("edges"); double counted between neighboring vertice u and v.
    !  ncon: # of weights for each vertex;
    !  int(in) vtxdist(p+1): Processor j stores vertices vtxdist(j):vtxdist(j+1)-1
    !  int (in) xadj(n+1), adjncy(m):
    !           locally, vertex j's neighboring vertices are adjncy(xadj(j):xadj(j+1)-1). adjncy points to global index;
    !  int(in) vwgt(ncon*n), adjwgt(m): weights at vertices and "edges". Format of adjwgt follows adjncy;
    !  int(in) wgtflag: 0: none (vwgt and adjwgt are NULL);
    !          1: edges (vwgt is NULL); 2: vertices (adjwgt is NULL); 3: both vertices & edges;
    !  int(in) numflag: 0: C-style numbering from 0; 1: FORTRAN style from 1;
    !  int(in) ndims: 2 or 3 (D);
    !  float(in) xyz(ndims*n): coordinate for vertex j is xyz(j*ndims:(j+1)*ndims-1);
    !  int(in)   nparts: # of desired sub-domains (usually nTasks);
    !  float(in) tpwgts(ncon*nparts): =1/nparts if sub-domains are to be of same size for each vertex weight;
    !  float(in) ubvec(ncon): imbalance tolerance for each weight;
    !  int(in)   options: additonal parameters for the routine (see above);
    !  int(out)  edgecut: # of edges that are cut by the partitioning;
    !  int(out)  part(): array size = # of local vertices. It stores indices of local vertices.
 
    ! write(1112+myrank,*) "Thread",myrank,"sum;vtxdist", sum(vtxdist), vtxdist
    ! write(1112+myrank,*) "Thread",myrank,"sum;xadj", sum(xadj), xadj
    ! write(1112+myrank,*) "Thread",myrank,"sum;adjncy", sum(adjncy), adjncy
 
    CALL real_mpi_barrier_pdlib(comm, "runParmetis, step 8")
 
    if(debugparmetis) then
      write(710+myrank,*) vtxdist, xadj, adjncy, &
           vwgt, & !vwgt - ignore weights
           adjwgt, & ! adjwgt - ignore weights
           wgtflag, &
           numflag,ndims,ncon,nparts,tpwgts,ubvec,options, &
           edgecut,part,comm
      call flush(710+myrank)
    endif
 
    !if(debugParmetis) write(710+myrank,*) "Run ParMETIS now..."
#ifdef W3_SCOTCH
    call scotch_parmetis_v3_partgeomkway(vtxdist, xadj, adjncy, &
         vwgt, & !vwgt - ignore weights
         adjwgt, & ! adjwgt - ignore weights
         wgtflag, &
         numflag,ndims,xyz,ncon,nparts,tpwgts,ubvec,options, &
         edgecut,part, comm,ref)
#endif
 
#ifdef W3_METIS
    call parmetis_v3_partgeomkway(vtxdist, xadj, adjncy, &
         vwgt, & !vwgt - ignore weights
         adjwgt, & ! adjwgt - ignore weights
         wgtflag, &
         numflag,ndims,xyz,ncon,nparts,tpwgts,ubvec,options, &
         edgecut,part, comm)
#endif
 
 
    CALL real_mpi_barrier_pdlib(comm, "runParmetis, step 9")
 
    if(ntasks == 1) then
      !    write(*,*) myrank, "minval part", minval(part)
      if(minval(part) == 0) then
        part(:) = part(:) + 1
      endif
    endif
 
    ! write(*,*) myrank, "edge cuted", edgecut
 
    ! Collect the parmetis data from all threads
    ! and create a global node to domain number mapping
    allocate(node2domain(np_global),stat=stat)
    if(stat/=0) call parallel_abort(' node2domain allocation failure')
    !
    call mpi_allgatherv(part, np, itype, node2domain, np_perproc, np_perprocsum, itype, comm, ierr)
    if(ierr/=mpi_success) call parallel_abort('mpi_allgatherv ',ierr)
    !
    do i = 1, np_global
      node => nodes_global(i)
      node%domainID = node2domain(node%id_global)
    end do
    !    CALL REAL_MPI_BARRIER_PDLIB(comm, "runParmetis, step 10")
 
 
    ! write out partition info for katerfempresenter
    if(debugpartition) write(600,*) node2domain
    !    Print *, 'runparmetis step 1'
 
    if(allocated(xadj))        deallocate(xadj)
    !    Print *, 'runparmetis step 2'
    if(allocated(adjncy))      deallocate(adjncy)
    !    Print *, 'runparmetis step 3'
    if(allocated(part))        deallocate(part)
    !    Print *, 'runparmetis step 4'
    if(allocated(vwgt))        deallocate(vwgt)
    !    Print *, 'runparmetis step 5'
    if(allocated(adjwgt))      deallocate(adjwgt)
    !    Print *, 'runparmetis step 6'
    if(allocated(xyz))         deallocate(xyz)
    !    Print *, 'runparmetis step 7'
    if(allocated(tpwgts))      deallocate(tpwgts)
    !    Print *, 'runparmetis step 8'
    if(allocated(ubvec))       deallocate(ubvec)
    !    Print *, 'runparmetis step 9'
    if(allocated(vtxdist))     deallocate(vtxdist)
    !    Print *, 'runparmetis step 10'
    if(allocated(node2domain)) deallocate(node2domain)
    !    Print *, 'runparmetis step 11'
  end subroutine runparmetis
 
  !------------------------------------------------------------------------
  ! with the new data from parmetis, recalculate some variables
  !------------------------------------------------------------------------
 
  ! Create for every boundary node a dummy node and for every dummy node two local elements
  ! extend x,y,z, INE
  ! rebuild  NCONE, CONE, NCONN, CONN
  ! alter nd, nde
  !subroutine dummyNodes
  ! use yowExchangeModule
  ! use yowNodepool, only: npa, nb, nd, PDLIB_NCONN, PDLIB_CONN, boundaryNodes, nodes, x, y, z, outwardNormal, XY, dummy2boundary, BNDprevnext
  ! use yowElementpool, only: ne, nde, INE, NCONE, CONE
  ! implicit none
  ! integer :: i, IP, IPprev, IPnext, newIP, newIE
  ! integer, allocatable :: INEnew(:,:)
  ! real(rkind) :: vec(2), newPoint(2), newLength
  ! real(rkind), allocatable :: xNew(:), yNew(:), zNew(:)
  !
  !    nd = nb
  !
  !    allocate(dummy2boundary(nd))
  !
  !    nde = 2*nd
  !    allocate(INEnew(3,ne+nde))
  !    INEnew(:,1:ne) = INE(:,1:ne)
  !
  !    allocate(xNew(npa+nd), yNew(npa+nd), zNew(npa+nd))
  !    xNew(1:npa) = x(1:npa)
  !    yNew(1:npa) = y(1:npa)
  !    zNew(1:npa) = z(1:npa)
  !
  !    nd = 0
  !    nde = 0
  !   do i=1, nb
  !     nd = nd + 1
  !     IP = boundaryNodes(i)
  !     dummy2boundary(nd) = IP
  !
  !     call BNDprevnext(IP, IPprev, IPnext)
  !     newLength = max(distanceTo(XY(IP), XY(IPprev)), distanceTo(XY(IP), XY(IPnext)))
  !     newPoint  = XY(IP) + outwardNormal(:,i) * newLength
  !    newIP = npa+i
  !    xNew(newIP) = newPoint(1)
  !    yNew(newIP) = newPoint(2)
  !    zNew(newIP) = z(IP)
  !    NCONN(IP) = NCONN(IP) + 1
  !    CONN(NCONN(IP),IP) = newIP
  !
  !    nde = nde +1
  !     newIE = ne+nde
  !     INEnew(1,newIE) = IP
  !     INEnew(2,newIE) = IPprev
  !     INEnew(3,newIE) = newIP
  !     NCONE(IP) = NCONE(IP) +1
  !     CONE(NCONE(IP),IP) = newIE
  !
  !     nde = nde +1
  !     newIE = ne+nde
  !     INEnew(1,newIE) = IP
  !     INEnew(2,newIE) = newIP
  !     INEnew(3,newIE) = IPnext
  !     NCONE(IP) = NCONE(IP) +1
  !     CONE(NCONE(IP),IP) = newIE
  !   end do
  !
  !   call move_alloc(INEnew, INE)
  !   call move_alloc(xNew, x)
  !   call move_alloc(yNew, y)
  !   call move_alloc(zNew, z)
  !
  !   ! Die Anzahl der angeschlossenen Knoten fuer Geisterknoten wird fuer findBoundaryNodes and Dummynodes gebraucht.
  !   call exchange(NCONN)
  !   call exchange(NCONE)
  !end subroutine
 
  subroutine postpartition
    use yowerr,    only: parallel_abort
    use yowdatapool, only: myrank, ntasks
    use yownodepool, only: np_global, np, nodes_global, nodes, np_perproc, np_perprocsum, iplg, ipgl, t_node
    use yowsidepool, only: ns
 
    integer :: i, j, stat
    type(t_node), pointer :: node
 
    ! determine how many nodes now belong to which thread
    ! and set the nodes local id
    np_perproc = 0
    do i = 1, np_global
      ! fortran counts from 1. np_perProc from 0
      np_perproc(nodes_global(i)%domainID-1) = np_perproc(nodes_global(i)%domainID-1)+1
      ! set the new local id
      nodes_global(i)%id = np_perproc(nodes_global(i)%domainID-1)
    end do
 
    np_perprocsum(0) = 0
    do i = 1, ntasks-1
      np_perprocsum(i) = np_perprocsum(i-1) + np_perproc(i-1)
    end do
 
    np = np_perproc(myrank)
 
    ! create the new node local to global mapping iplg. This is not the final one.
    ! create a iplg with ghost nodes in findGhostNodes
    if(allocated(iplg)) deallocate(iplg)
    allocate(iplg(np), stat=stat)
    if(stat/=0) call parallel_abort('iplg second allocation failure')
    iplg = 0
 
    if(allocated(ipgl)) deallocate(ipgl)
    allocate(ipgl(np_global), stat=stat)
    if(stat/=0) call parallel_abort('ipgl allocation failure')
    ipgl = 0
 
    j = 1
    do i = 1, np_global
      node => nodes_global(i)
      if(node%domainID == myrank+1) then
        iplg(j) = i
        ipgl(i) = j
        j = j + 1
      endif
    end do
 
    ! calc # sides local again, because the nodes now belongs to another domain
    ns = 0
    do i = 1, np
      node => nodes(i)
      ns = ns + node%nConnNodes
    end do
  end subroutine postpartition
 
 
  !----------------------------------------------------------------------
  ! find the ghost nodes of the local domain
  !----------------------------------------------------------------------
 
  subroutine findghostnodes
    use yowerr,    only: parallel_abort
    use yowdatapool, only: myrank
    use yownodepool, only: t_node, np, nodes, ghosts, nodes_global, ng, ghostlg, ghostgl, npa, np_global, iplg
 
    integer :: i, j, k, stat
    type(t_node), pointer :: node, nodeNeighbor, nodeGhost
    integer, save, allocatable :: ghostTemp(:)
#ifdef W3_DEBUGINIT
    print *, 'Passing in findGhostNodes'
#endif
 
    ! iterate over all local nodes and look at their neighbors
    ! has the neighbor another domain id, than it is a ghost
 
    ! implementation is some different
    ! loop over all nodes to get the # of ghost nodes
    ! allocate space
    ! loop a second time to insert the ghost nodes
 
    ! first loop. find out how many ghost nodes we have (with double entries)
    ng = 0
    do i = 1, np
      node => nodes(i)
      do j = 1, node%nConnNodes
        nodeneighbor => node%connNodes(j)
        if(nodeneighbor%domainID /= node%domainID) then
          ! yes, we found a ghost
          ng = ng + 1
        end if
      end do
    end do
 
    allocate(ghosttemp(ng), stat=stat)
    if(stat/=0) call parallel_abort('ghostTemp allocation failure')
 
#ifdef W3_DEBUGINIT
    print *, 'np_global=', np_global
#endif
    IF (allocated(ghostgl)) THEN
      print *, 'ghostgl is already allocated'
    END IF
    allocate(ghostgl(np_global), stat=stat)
    if(stat/=0) call parallel_abort('ghostgl allocation failure')
    ghostgl = 0
 
    ! second loop. fill ghostlg. ignore double entries
    ng = 0
    ! iterate over all local nodes
    do i = 1, np
      node => nodes(i)
 
      ! check their neighbors
      secondloop: do j = 1, node%nConnNodes
        nodeneighbor => node%connNodes(j)
 
        if(nodeneighbor%domainID /= node%domainID) then
          ! yes, we found a ghost
          ! check if this ghost is allready in the ghost list
          do k = 1, ng
            nodeghost => nodes_global(ghosttemp(k))
            if(nodeneighbor%id_global == nodeghost%id_global) then
              ! yes, we allready know this ghost.
              ! check the next neighbor
              cycle secondloop
            end if
          end do
 
          ! no we don't know this ghost. insert it
          ng = ng + 1
          ghosttemp(ng) = nodeneighbor%id_global
          ghostgl(nodeneighbor%id_global) = ng
        end if
      end do secondloop
    end do
 
    ! reallocate the gosttemp array becouse it is longer then the new ng
    if(allocated(ghostlg)) deallocate(ghostlg)
    allocate(ghostlg(ng), stat=stat)
    if(stat/=0) call parallel_abort('ghostlg allocation failure')
    ghostlg = ghosttemp(1:ng)
    deallocate(ghosttemp)
 
    npa = np + ng
 
    ! check if ghostlg contains only uniqe values
    do i=1, ng
      do j=i+1, ng
        if(ghostlg(i) == ghostlg(j)) then
          write(*,*) "double global ghost id in ghostlg(i,j)", i, j
          stop "double global ghost id in ghostlg(i,j)"
        endif
      end do
    end do
 
 
    ! create the new node local to global mapping iplg with ghost. final one.
    if(allocated(iplg)) deallocate(iplg)
    allocate(iplg(npa), stat=stat)
    if(stat/=0) call parallel_abort('iplg second allocation failure')
    iplg = 0
 
    j = 1
    do i = 1, np_global
      node => nodes_global(i)
      if(node%domainID == myrank+1) then
        iplg(j) = i
        j = j + 1
      endif
    end do
 
    iplg(np+1: npa) = ghostlg(1:ng)
  end subroutine findghostnodes
  !-------------------------------------------------------------------------------
  ! find the number of connected domains and their ghosts
  !-------------------------------------------------------------------------------
 
  subroutine findconndomains
    use yowerr,          only: parallel_abort
    use yownodepool,       only: ghosts, ng, t_node
    use yowdatapool,       only: ntasks, myrank
    use yowexchangemodule, only: neighbordomains, initnbrdomains
 
    integer :: i, stat, itemp
    type(t_node), pointer :: ghost
 
    ! # of ghost per neighbor domain
    integer, allocatable :: numberGhostPerNeighborDomainTemp(:)
    ! look up table. domainID to neighbor (ID)
    integer, allocatable :: domainID2NeighborTemp(:)
 
    ! Part 1) find # neighbor domains
 
    ! allocate this array with a fixed size of nTasks. even if we not have
    ! so many neighbor domains, nTasks will never be very large
    allocate(numberghostperneighbordomaintemp(ntasks), stat=stat)
    if(stat/=0) call parallel_abort('numberGhostPerNeighborDomainTemp allocation failure')
    numberghostperneighbordomaintemp = 0
 
    allocate(domainid2neighbortemp(ntasks), stat=stat)
    if(stat/=0) call parallel_abort('domainID2NeighborTemp allocation failure')
    domainid2neighbortemp = 0
 
 
    ! iterate over all ghost nodes an get their thread id
    itemp = 0
    do i = 1, ng
      ghost => ghosts(i)
 
      ! sum how many ghost belongs to the ghost domainID
      numberghostperneighbordomaintemp(ghost%domainID) = numberghostperneighbordomaintemp(ghost%domainID) + 1
 
      ! check if this ghost domainID is allready in the domains list
      if(domainid2neighbortemp(ghost%domainID) /= 0) then
        ! yes we have allready insert this domain id
        cycle
      end if
 
      ! no we dont know this domain id. insert it
      itemp = itemp + 1
      domainid2neighbortemp(ghost%domainID) = itemp
    end do
 
    ! Part 2)  assign the ghost nodes to their domains
    call initnbrdomains(itemp)
 
    do i = 1, ntasks
      if(numberghostperneighbordomaintemp(i) /= 0) then
        neighbordomains(domainid2neighbortemp(i) )%domainID = i
 
        allocate(neighbordomains(domainid2neighbortemp(i))%nodesToReceive(numberghostperneighbordomaintemp(i)), stat=stat)
        if(stat/=0) call parallel_abort('neighborDomains%ghosts allocation failure')
        neighbordomains(domainid2neighbortemp(i))%nodesToReceive = 0
      end if
    end do
 
    do i = 1, ng
      ghost => ghosts(i)
      itemp = domainid2neighbortemp(ghost%domainID)
 
      neighbordomains(itemp)%numNodesToReceive = neighbordomains(itemp)%numNodesToReceive + 1
      neighbordomains(itemp)%nodesToReceive(neighbordomains(itemp)%numNodesToReceive) = ghost%id_global
    end do
 
    if(allocated(numberghostperneighbordomaintemp)) deallocate(numberghostperneighbordomaintemp)
    if(allocated(domainid2neighbortemp)) deallocate(domainid2neighbortemp)
  end subroutine findconndomains
 
 
  !-------------------------------------------------------------------------------
  ! exchange Ghost Ids so every thread knows which nodes he has to send to the
  ! other parition
  !-------------------------------------------------------------------------------
 
  subroutine exchangeghostids
    use yowerr
    use yownodepool,       only: np, t_node, nodes
    use yowdatapool,       only: ntasks, myrank, comm
    use yowexchangemodule, only: neighbordomains, nconndomains, creatempitypes
    use mpi
 
    integer :: i, j, k
    integer :: ierr
    ! uniq tag that identify the sender and which information he sends
    integer :: tag
    ! we use non-blocking send and recv subroutines
    ! store the send status
    integer :: sendRequest(nConnDomains)
    ! store the revc status
    integer :: recvRequest(nConnDomains)
    ! status to verify if one communication fails or not
    integer :: status(MPI_STATUS_SIZE, nConnDomains);
 
 
    type(t_node), pointer :: node
 
    ! send to all domain neighbors how many ghosts nodes we want from him and which ones
    do i=1, nconndomains
      ! create a uniq tag for this domain
      tag = neighbordomains(i)%domainID*10 + 1
      ! send to the neighbor how many ghost nodes we want from him
      call mpi_isend(neighbordomains(i)%numNodesToReceive, &
           1, &
           mpi_int, &
           neighbordomains(i)%domainID-1, &
           tag, &
           comm, &
           sendrequest(i), &
           ierr);
      if(ierr/=mpi_success) then
        write(*,*) "mpi send failure"
      endif
 
      tag = neighbordomains(i)%domainID*10 + 2
      ! send to the neighbor which ghost nodes we want from him
      call mpi_isend(neighbordomains(i)%nodesToReceive, &
           neighbordomains(i)%numNodesToReceive, &
           mpi_int, &
           neighbordomains(i)%domainID-1, &
           tag, &
           comm, &
           sendrequest(i), &
           ierr);
      if(ierr/=mpi_success) then
        write(*,*) "mpi send failure"
      endif
 
      ! receive from neighbor how many ghost nodes we have to send him
      tag = (myrank+1)*10 + 1
      call mpi_irecv(neighbordomains(i)%numNodesToSend, &
           1, &
           mpi_int, &
           neighbordomains(i)%domainID-1, &
           tag, &
           comm, &
           recvrequest(i), &
           ierr)
      if(ierr/=mpi_success) then
        write(*,*) "mpi recv failure"
      endif
    end do
 
    ! wait for communication end
    call mpi_waitall(nconndomains, recvrequest, status, ierr)
 
    ! test for all neighbor domains
    do i=1, nconndomains
      ! test if the neighbor wants more ghost nodes than we have
      if(neighbordomains(i)%numNodesToSend > np) then
        write(*,'(i5, a, i5, a, i5,a, i5, a, /)', advance='no') myrank, " ERROR neighbordomain ", neighbordomains(i)%domainID, &
             " wants ", neighbordomains(i)%numNodesToSend, &
             " nodes, but we have only ", np, " nodes"
        CALL abort("")
      end if
    end do
 
    ! receive from all neighbor domains which nodes we must send him
    do i=1, nconndomains
      allocate(neighbordomains(i)%nodesToSend(neighbordomains(i)%numNodesToSend))
      neighbordomains(i)%nodesToSend = 0
 
      ! receive from neighbor which nodes we must send
      tag = (myrank+1)*10 + 2
      call mpi_irecv(neighbordomains(i)%nodesToSend, &
           neighbordomains(i)%numNodesToSend, &
           mpi_int, &
           neighbordomains(i)%domainID-1, &
           tag, &
           comm, &
           recvrequest(i), &
           ierr)
      if(ierr/=mpi_success) then
        CALL parallel_abort("mpi recv failure", ierr)
      endif
    end do
 
    ! wait for communication end
    call mpi_waitall(nconndomains, recvrequest, status, ierr)
 
    ! test for all neighbor domains
    do i=1, nconndomains
      ! test if the neighbor wants nodes that we don't own
      outerloop: do j=1, neighbordomains(i)%numNodesToSend
        ! compare with all local nodes
        do k=1, np
          node => nodes(k)
          if(node%id_global == neighbordomains(i)%nodesToSend(j)) then
            cycle outerloop
          end if
        end do
        write(*,*) myrank, "Neighbordomain", neighbordomains(i)%domainID, &
             " want Node", neighbordomains(i)%nodesToSend(j), &
             " but we don't own this node"
        stop
      end do outerloop
    end do
 
    call creatempitypes()
  end subroutine exchangeghostids
 
  subroutine postpartition2(INE_global)
    use yowelementpool, only: ne, ne_global, ine, belongto, ielg
    use yowerr,         only: parallel_abort
    use yowdatapool,    only: myrank
    use yownodepool,    only: np_global, np, nodes_global, iplg, t_node, ghostlg, ng, npa
    use yownodepool,    only: x, y, z
    use w3gdatmd,       only: xgrd, ygrd, zb
 
    integer, intent(in) :: INE_global(3,ne_global)
 
    integer :: i, j, k, stat, IP_glob
    type(t_node), pointer :: node
    logical :: assigned
 
    ! to find the local elements, iterate over all global elements and check their node domain IDs
 
    ! step 1: calc the number of local elements
    ne = 0
    do i=1, ne_global
      if (belongto(ine_global(:,i))) then
        ne = ne +1
      endif
    end do
 
    ! step 2: fill the local element index array
    if(allocated(ine)) deallocate(ine)
    allocate(ine(3, ne), stat=stat)
    if(stat/=0) call parallel_abort('INE allocation failure')
    ine = 0
 
    ne = 0
    do i=1, ne_global
      ! yes, this element belongs to this domain
      if (belongto(ine_global(:,i))) then
        ne = ne + 1
        do j=1, 3
          assigned = .false.
          node => nodes_global(ine_global(j,i))
          if(node%domainID == myrank+1) then
            ine(j, ne) = node%id
            assigned = .true.
          else
            ! the element have a ghost node
            do k=1, ng
              if(node%id_global == ghostlg(k)) then
                ! conversion: the ghost nodes are stored behind the local nodes.
                if(ine(j,ne) /= 0) then
                  write(*,*) "will write to INE(j, ne) but there is allready a value", j, ne, ine(j, ne)
                endif
                ine(j, ne) = np + k
                node%id = np+k
                assigned = .true.
                !                 write(*,*) myrank, "node to ele", node%id_global-1, i-1, np+k
                exit
              endif
            end do
          endif
          if(assigned .eqv. .false.) then
            write(*,*) "Can't assign global node to INE", node%id_global
          endif
        end do
      endif
    end do
 
    ! check if INE contains 0.
    do i=1, ne
      if(minval(abs(ine(:,i))) == 0) then
        write(*,*) "0 in INE ne=", ne
        stop "0 in INE"
      endif
    end do
 
    if(maxval(ine) /= npa) then
      write(*,*) "MAXVAL(INE) /= npa ERROR?"
    endif
 
    ! create element local to global mapping ielg
    if(allocated(ielg)) deallocate(ielg)
    allocate(ielg(ne), stat=stat)
    if(stat/=0) call parallel_abort('ielg allocation failure')
    ielg = 0
 
    j = 0
    do i=1, ne_global
      if (belongto(ine_global(:,i))) then
        j = j +1
        ielg(j) = i
      end if
    end do
 
    ! fill the local x,y arrays
    if(allocated(x)) deallocate(x)
    allocate(x(npa), stat=stat)
    if(stat/=0) call parallel_abort('x allocation failure')
 
    if(allocated(y)) deallocate(y)
    allocate(y(npa), stat=stat)
    if(stat/=0) call parallel_abort('y allocation failure')
 
    if(allocated(z)) deallocate(z)
    allocate(z(npa), stat=stat)
    if(stat/=0) call parallel_abort('z allocation failure')
 
    do i=1, np
      ip_glob = iplg(i)
      x(i) = xgrd(1,ip_glob)
      y(i) = ygrd(1,ip_glob)
      z(i) = zb(ip_glob)
    end do
 
    do i=1, ng
      ip_glob = ghostlg(i)
      x(np+i) = xgrd(1,ip_glob)
      y(np+i) = ygrd(1,ip_glob)
      z(np+i) = zb(ip_glob)
    end do
 
  end subroutine postpartition2
  !**********************************************************************
  !*                                                                    *
  !**********************************************************************
  subroutine computetria_ien_si_ccon
    use yowelementpool, only: ne, ne_global, ine, ielg
    use yowexchangemodule, only : pdlib_exchange1dreal
    use yowerr,       only: parallel_abort
    use yowdatapool,    only: myrank
    use yownodepool,    only: np_global, np, iplg, t_node, ghostlg, ng, npa
    use yownodepool,    only: x, y, z, pdlib_si, pdlib_ien, pdlib_tria, pdlib_ccon, pdlib_tria03
 
    integer I1, I2, I3, stat, IE, NI(3)
    real  :: DXP1, DXP2, DXP3, DYP1, DYP2, DYP3, DBLTMP, TRIA03
    logical :: CROSSES_DATELINE
 
    allocate(pdlib_si(npa), pdlib_ccon(npa), pdlib_ien(6,ne), pdlib_tria(ne), pdlib_tria03(ne), stat=stat)
    if(stat/=0) call parallel_abort('SI allocation failure')
 
    pdlib_si(:)   = 0.0d0 ! Median Dual Patch Area of each Node
    pdlib_ccon(:) = 0     ! Number of connected Elements
    DO ie = 1 , ne
      i1 = ine(1,ie)
      i2 = ine(2,ie)
      i3 = ine(3,ie)
      ni = ine(:,ie)
 
      dxp1=x(i2) - x(i1)
      dyp1=y(i2) - y(i1)
      dxp2=x(i3) - x(i2)
      dyp2=y(i3) - y(i2)
      dxp3=x(i1) - x(i3)
      dyp3=y(i1) - y(i3)
      CALL element_crosses_dateline(dxp1, dxp2, dxp3, crosses_dateline)
      IF (crosses_dateline) THEN
        CALL correct_dx_gt180(dxp1)
        CALL correct_dx_gt180(dxp2)
        CALL correct_dx_gt180(dxp3)
      ENDIF
 
      pdlib_ien(1,ie) = - dyp2
      pdlib_ien(2,ie) =   dxp2
      pdlib_ien(3,ie) = - dyp3
      pdlib_ien(4,ie) =   dxp3
      pdlib_ien(5,ie) = - dyp1
      pdlib_ien(6,ie) =   dxp1
      dbltmp = (dxp3*dyp1 - dyp3*dxp1)*0.5
      pdlib_tria(ie) = dbltmp
      IF (pdlib_tria(ie) .lt. tiny(1.)) THEN
        WRITE(*,*) pdlib_ien(:,ie)
        WRITE(*,*)
        WRITE(*,*) 'AREA SMALLER ZERO IN PDLIB', ie, ne, pdlib_tria(ie)
        stop
      ENDIF
 
      pdlib_ccon(i1) = pdlib_ccon(i1) + 1
      pdlib_ccon(i2) = pdlib_ccon(i2) + 1
      pdlib_ccon(i3) = pdlib_ccon(i3) + 1
      tria03         = pdlib_tria(ie)/3.d0
      pdlib_si(i1) = pdlib_si(i1) + tria03
      pdlib_si(i2) = pdlib_si(i2) + tria03
      pdlib_si(i3) = pdlib_si(i3) + tria03
      pdlib_tria03(ie) = tria03
    ENDDO
    CALL pdlib_exchange1dreal(pdlib_si)
  end subroutine computetria_ien_si_ccon
  !**********************************************************************
  !*                                                                    *
  !**********************************************************************
  subroutine element_crosses_dateline(RX1, RX2, RX3, CROSSES_DATELINE)
    !   Purpose: understanding if an element crosses the dateline.
    !   An element crossing the dateline has, e.g. a node with lon < 180
    !   and another 2 with lon > -180
 
    REAL(rkind),  INTENT(IN)  :: RX1, RX2, RX3
    LOGICAL, INTENT(OUT) :: CROSSES_DATELINE
    INTEGER :: R1GT180, R2GT180, R3GT180
    r1gt180 = merge(1, 0, abs(rx1).GT.180)
    r2gt180 = merge(1, 0, abs(rx2).GT.180)
    r3gt180 = merge(1, 0, abs(rx3).GT.180)
    ! if R1GT180+R2GT180+R3GT180 .eq. 0 the element does not cross the dateline
    ! if R1GT180+R2GT180+R3GT180 .eq. 1 the element contains the pole
    ! if R1GT180+R2GT180+R3GT180 .eq. 2 the element crosses the dateline
    crosses_dateline = r1gt180+r2gt180+r3gt180 .EQ. 2
  end subroutine element_crosses_dateline
  !**********************************************************************
  !*                                                                    *
  !**********************************************************************
  subroutine correct_dx_gt180(DXP)
    !   Purpose: the absolute zonal distance between 2 points is always <= 180
    !            This subroutine corrects the zonal distance to satifsy
    !            this requirement
 
    REAL(rkind), INTENT(INOUT) :: DXP
    IF (dxp .le. -180) THEN
      dxp=dxp + 360
    END IF
    IF (dxp .ge. 180) THEN
      dxp=dxp - 360
    END IF
  end subroutine correct_dx_gt180
  !**********************************************************************
  !*                                                                    *
  !**********************************************************************
  subroutine computeia_ja_posi_nnz
    use yowelementpool, only: ne, ne_global, ine, ielg
    use yowerr,       only: parallel_abort
    use yowdatapool,    only: myrank
    use yownodepool,    only: np_global, np, nodes_global, iplg, t_node, ghostlg, ng, npa
    use yownodepool,    only: pdlib_ccon, pdlib_ia, pdlib_ja, pdlib_ja_ie, pdlib_ia_p, pdlib_ja_p
    use yownodepool,    only: pdlib_nnz, pdlib_posi, pdlib_ie_cell, pdlib_pos_cell, pdlib_ie_cell2
    use yownodepool,    only: pdlib_pos_cell2, pdlib_i_diag
 
    integer CHILF(npa)
    integer istat
    integer MAXMNECON
    integer IE, J, I, IP, K, IP_J, IP_K, IP_I
    integer I1, I2, I3, POS, POS_J, POS_K
    integer COUNT_MAX
    integer, allocatable :: CELLVERTEX(:,:,:), PTABLE(:,:)
    integer :: ITMP(npa)
    maxmnecon  = maxval(pdlib_ccon)
    ALLOCATE(cellvertex(npa,maxmnecon,2), stat=istat)
    IF (istat/=0) CALL parallel_abort('ComputeIA_JA_POSI_NNZ, allocate error 4')
    !
    cellvertex(:,:,:) = 0
    chilf             = 0
    DO ie = 1, ne
      DO j=1,3
        i = ine(j,ie)
        chilf(i) = chilf(i)+1
        cellvertex(i,chilf(i),1) = ie
        cellvertex(i,chilf(i),2) = j
      END DO
    ENDDO
    !
    !        Emulates loop structure and counts max. entries in the different pointers that have to be designed
    !
    j = 0
    DO ip = 1, npa
      DO i = 1, pdlib_ccon(ip)
        j = j + 1
      END DO
    END DO
    count_max = j
    ALLOCATE (pdlib_ie_cell(count_max), pdlib_pos_cell(count_max), stat=istat)
    IF (istat/=0) CALL parallel_abort('wwm_fluctsplit, allocate error 5a')
    ALLOCATE (pdlib_ie_cell2(maxmnecon,npa), pdlib_pos_cell2(maxmnecon, npa), stat=istat)
    IF (istat/=0) CALL parallel_abort('wwm_fluctsplit, allocate error 5b')
    ! Just a remapping from CELLVERTEX ... Element number in the
    ! order of the occurence in the loop during runtime
    pdlib_ie_cell  = 0
    ! Just a remapping from CELLVERTEX ... Position of the node
    ! in the Element index -"-
    pdlib_pos_cell = 0
    j = 0
    DO ip = 1, npa
      DO i = 1, pdlib_ccon(ip)
        j = j + 1
        pdlib_ie_cell(j)      = cellvertex(ip,i,1)
        pdlib_pos_cell(j)     = cellvertex(ip,i,2)
        pdlib_ie_cell2(i,ip)  = cellvertex(ip,i,1)
        pdlib_pos_cell2(i,ip) = cellvertex(ip,i,2)
      END DO
    END DO
    deallocate(cellvertex)
 
    ALLOCATE(ptable(count_max,7), stat=istat)
    IF (istat/=0) CALL parallel_abort('wwm_fluctsplit, allocate error 6')
    ALLOCATE(pdlib_ja_ie(3,3,ne), stat=istat)
    IF (istat/=0) CALL parallel_abort('wwm_fluctsplit, allocate error 6.1')
 
    j = 0
    ptable(:,:) = 0.
    DO ip = 1, npa
      DO i = 1, pdlib_ccon(ip)
        j = j + 1
        ie    = pdlib_ie_cell(j)
        pos   = pdlib_pos_cell(j)
        i1 = ine(1,ie)
        i2 = ine(2,ie)
        i3 = ine(3,ie)
        IF (pos == 1) THEN
          pos_j = 2
          pos_k = 3
        ELSE IF (pos == 2) THEN
          pos_j = 3
          pos_k = 1
        ELSE
          pos_j = 1
          pos_k = 2
        END IF
        ip_i = ip
        ip_j = ine(pos_j,ie)
        ip_k = ine(pos_k,ie)
        ptable(j,1) = ip_i ! Node numbers of the connected elements
        ptable(j,2) = ip_j
        ptable(j,3) = ip_k
        ptable(j,4) = pos  ! Position of the nodes in the element index
        ptable(j,5) = pos_j
        ptable(j,6) = pos_k
        ptable(j,7) = ie   ! Element numbers same as PDLIB_IE_CELL
      END DO
    END DO
    !
    ! Count number of nonzero entries in the matrix ...
    ! Basically, each connected element may have two off-diagonal
    ! contribution and one diagonal related to the connected vertex itself ...
    !
    j = 0
    pdlib_nnz = 0
    itmp = 0
    DO ip = 1, npa
      itmp = 0
      DO i = 1, pdlib_ccon(ip)
        j = j + 1
        ip_j  = ptable(j,2)
        ip_k  = ptable(j,3)
        itmp(ip)   = 1
        itmp(ip_j) = 1
        itmp(ip_k) = 1
      END DO
      pdlib_nnz = pdlib_nnz + sum(itmp)
    END DO
    !
    ! Allocate sparse matrix pointers using the Compressed Sparse Row Format CSR ... this is now done only of npa nodes
    ! The next step is to do it for the whole Matrix npa * MSC * MDC
    ! see ...:x
    !
    ALLOCATE (pdlib_ja(pdlib_nnz), pdlib_ia(npa+1), pdlib_ja_p(pdlib_nnz), stat=istat)
    IF (istat/=0) CALL parallel_abort('wwm_fluctsplit, allocate error 6a')
    ALLOCATE (pdlib_ia_p(npa+1), pdlib_posi(3,count_max), pdlib_i_diag(npa), stat=istat)
    IF (istat/=0) CALL parallel_abort('wwm_fluctsplit, allocate error 6b')
    pdlib_ja = 0
    pdlib_ia = 0
    pdlib_ja_p = 0
    pdlib_ia_p = 0
    pdlib_posi = 0
    ! Points to the position of the matrix entry in the mass matrix
    ! according to the CSR matrix format see p. 124
    j = 0
    k = 0
    pdlib_ia(1) = 1
    pdlib_ia_p(1) = 0
    DO ip = 1, npa ! Run through all rows
      itmp=0
      DO i = 1, pdlib_ccon(ip) ! Check how many entries there are ...
        j = j + 1
        ip_j = ptable(j,2)
        ip_k = ptable(j,3)
        itmp(ip)   = 1
        itmp(ip_j) = 1
        itmp(ip_k) = 1
      END DO
      DO i = 1, npa ! Run through all columns
        IF (itmp(i) .GT. 0) THEN
          k = k + 1
          pdlib_ja(k) = i
          pdlib_ja_p(k) = i-1
        END IF
      END DO
      pdlib_ia(ip + 1) = k + 1
      pdlib_ia_p(ip + 1) = k
    END DO
    pdlib_posi = 0
    j = 0
    DO ip = 1, npa
      DO i = 1, pdlib_ccon(ip)
        j = j + 1
        ip_j  = ptable(j,2)
        ip_k  = ptable(j,3)
        DO k = pdlib_ia(ip), pdlib_ia(ip+1) - 1
          IF (ip   == pdlib_ja(k)) pdlib_posi(1,j)  = k
          IF (ip   == pdlib_ja(k)) pdlib_i_diag(ip) = k
          IF (ip_j == pdlib_ja(k)) pdlib_posi(2,j)  = k
          IF (ip_k == pdlib_ja(k)) pdlib_posi(3,j)  = k
        END DO
      END DO
    END DO
    j=0
    DO ip=1,npa
      DO i = 1, pdlib_ccon(ip)
        j = j + 1
        ie    =  pdlib_ie_cell(j)
        pos   =  pdlib_pos_cell(j)
        i1    =  pdlib_posi(1,j)
        i2    =  pdlib_posi(2,j)
        i3    =  pdlib_posi(3,j)
        pdlib_ja_ie(pos,1,ie) = i1
        pdlib_ja_ie(pos,2,ie) = i2
        pdlib_ja_ie(pos,3,ie) = i3
      END DO
    END DO
    deallocate(ptable)
  end subroutine computeia_ja_posi_nnz
  !**********************************************************************
  !*                                                                    *
  !**********************************************************************
  subroutine finalizepd()
    use yowexchangemodule, only: finalizeexchangemodule
    use yownodepool,    only: finalizenodepool
    use yowelementpool, only: finalizeelementpool
    use yowrankmodule,  only: finalizerankmodule
 
    call finalizerankmodule()
    call finalizeexchangemodule()
    call finalizeelementpool()
    call finalizenodepool()
  end subroutine finalizepd
 
end module yowpdlibmain
!**********************************************************************
!*                                                                    *
!**********************************************************************