readzincoord Subroutine

  subroutine readzincoord
  use modglobal,   only :  kb,ke,kh,ifinput,zf,zh,ylen,jb,je,dy
  use modmpi,      only :  myid, mpi_integer,comm3d,mpierr,my_real,nprocs
  implicit none
  character(72) chmess
  character(20) namezinlet
  character(20) namezinfo
  integer ierr,k,kk,kmaxin,j,jj
  real ysizeproc

      namelist/INFO/nprocsinl,jgtotinl,kmaxin,dtin,wtop,totalreadu

  namezinlet = 'zgrid.inl'
  namezinfo  = 'zgrid.inf'

    if (myid==0) then


      open(ifinput,file=namezinfo,status='old',iostat=ierr)
        if (ierr /= 0) then
          write(0, *) 'ERROR: zgrid.inf does not exist'
          stop 1
        end if
        read (ifinput,INFO,iostat=ierr)
        if (ierr > 0) then
          write(0, *) 'Problem in zgrid.inf INFO'
          write(0, *) 'iostat error: ', ierr
          stop 1
        endif
        write(6,INFO)
      close(ifinput)
    end if
    kbin = 0
    kein = kmaxin-1
    call MPI_BCAST(nprocsinl,1,MPI_INTEGER,0,comm3d,mpierr)
    call MPI_BCAST(jgtotinl ,1,MPI_INTEGER,0,comm3d,mpierr)
    call MPI_BCAST(kbin     ,1,MPI_INTEGER,0,comm3d,mpierr)
    call MPI_BCAST(kein     ,1,MPI_INTEGER,0,comm3d,mpierr)
    call MPI_BCAST(dtin     ,1,MY_REAL,0,comm3d,mpierr)
    call MPI_BCAST(wtop     ,1,MY_REAL,0,comm3d,mpierr)
    call MPI_BCAST(totalreadu ,1,MY_REAL,0,comm3d,mpierr)
    allocate(zhin(kbin   :kein+1))
    allocate(zfin(kbin   :kein+1))
    allocate(dzfin(kbin-1:kein+1))
    allocate(dzhin(kbin  :kein+1))

    if (myid==0) then
      write(6,*) 'loading ',namezinlet
      open (ifinput,file=namezinlet)
      read(ifinput,'(a72)') chmess
      read(ifinput,'(a72)') chmess

      do k=kbin,kein
        read(ifinput,*) zfin(k)
      end do
      close(ifinput)

      zhin(kbin) = 0.0
      do k=kbin,kein
        zhin(k+1) = zhin(k) + 2.0*(zfin(k)-zhin(k))
      end do
      zfin(kein+kh)  = zfin(kein)+ 2.0*(zhin(kein+kh)-zfin(kein))

      do  k=kbin,kein
        dzfin(k) = zhin(k+1) - zhin(k)
      end do
      dzfin(kein+1) = dzfin(kein)
      dzfin(kbin-1) = dzfin(kbin)

      dzhin(kbin) = 2*zfin(kbin)
      do k=kbin+1,kein+kh
        dzhin(k) = zfin(k) - zfin(k-1)
      end do

    ! check if the inlet mesh and the simulation mesh differ
    do k=kb,min(ke,kein)
      if(abs(zfin(k)-zf(k)) > 1e-7) then
        lzinzsim = .false.
      end if
    enddo

    end if ! myid==0

  ! MPI broadcast kmax elements from zf
    call MPI_BCAST( zfin,kein-kbin+2,MY_REAL   ,0,comm3d,mpierr)
    call MPI_BCAST( zhin,kein-kbin+2,MY_REAL   ,0,comm3d,mpierr)
    call MPI_BCAST(dzfin,kein-kbin+3,MY_REAL   ,0,comm3d,mpierr)
    call MPI_BCAST(dzhin,kein-kbin+2,MY_REAL   ,0,comm3d,mpierr)
    call MPI_BCAST(lzinzsim,1,MPI_INTEGER      ,0,comm3d,mpierr)


    if (.not.lzinzsim) then
      if (myid==0) then
        write(6,*) 'zgrid.inl does not equal zgrid.inp: Inlet will be interpolated in z'
      end if
!      allocate(linlf(kbin:kein))
!      allocate(linuf(kbin:kein))
!      allocate(linlh(kbin:kein+1))
!      allocate(linuh(kbin:kein+1))
      allocate(linlf(kb:ke))
      allocate(linuf(kb:ke))
      allocate(linlh(kb:ke+1))
      allocate(linuh(kb:ke+1))
   ! zf
      do k=kb,ke
        do kk=kbin,kein
          if (zfin(kk) >= zf(k)) then
            linuf(k)  = kk
            linlf(k)  = kk-1
            exit
          elseif (kk==kein) then
            linuf(k)  = kein+1            ! this means extrapolation!
            linlf(k)  = kein-1            ! waarom niet ke? of wordt dit niet gebruikt?
          end if
        end do
      end do
   ! for w-components (zh)
      do k=kb,ke+1
        do kk=kbin,kein+1
          if (zhin(kk) >= zh(k)) then
            linuh(k)  = kk
            linlh(k) = kk-1
            exit
          elseif (kk==kein+1) then
            linuh(k)  = kein+2            ! this means extrapolation!
            linlh(k)  = kein
          end if
        end do
      end do

     else  ! lzinzsim  -> grids are equal
       if (myid==0) then
         write(6,*) 'zgrid.inl equals zgrid.inp: Inlet will not be interpolated in z'
       end if
     end if

! Now prepare everything for interpolation in y-direction
     jgbin = 1
     jgein = jgbin + jgtotinl -1
     jtotin = jgtotinl / nprocsinl
     jbin = 1
     jein = 1 + jtotin -1
     ysizeproc = ylen/nprocs
     dyin = ylen / jgtotinl
     jbdum = 1
     jtotdum = ceiling(ysizeproc/real(dyin))+1                    ! dummy indices
     jedum = jbdum + jtotdum-1

!     allocate(yf   (jb    :je))
     allocate(yf   (jb    :je+1))
     allocate(yh   (jb    :je+1))
     allocate(yfin (jgbin :jgein+1))
     allocate(yhin (jgbin-1 :jgein+1))
!     allocate(yfdum(jbdum :jedum))
     allocate(yfdum(jbdum :jedum+1))
     allocate(yhdum(jbdum :jedum+1))
     allocate(yloclowf(jb:je+1))
     allocate(ylocupf (jb:je+1))
     allocate(yloclowh(jb:je+1))
     allocate(ylocuph (jb:je+1))

    ! make global y-grid (equidistant) for inlet data

     do j = jgbin-1,jgein+1
       yhin(j) = (j-jgbin)*dyin
     end do
     do j = jgbin,jgein+1
       yfin(j) = yhin(j) + 0.5*dyin
     end do
    ! make new y-grid (equidistant)
     do j = jb,je+1
       yh(j) = myid*(ylen/nprocs) + (j-jb)*dy
       yf(j) = yh(j) + 0.5*dy
     end do

   ! check which original cells are needed for interpolation
      do j=jgein+1,jgbin,-1
        if (yhin(j)<= yh(jb)) then
          if (yfin(j)<= yf(jb)) then
            procinlo = floor(real(j-jgbin)/real(jtotin))      ! this is the first cell to consider
            filenumstart = procinlo
            jgbeg = j
            jbeg = j-(procinlo*jtotin)
!            jend = jbeg+jtotdum-1
            jj = j+jtotdum-1
!            procinup = floor((j-jgbin)/real(jtotin))
!            procinup = floor((j-jgbin+1)/real(jtotin))
            procinup = floor(real(jj-jgbin)/real(jtotin))
            filenumend = procinup
            jend = jj-(procinup*jtotin)
            procinup = procinup-floor(real(procinup)/real(nprocsinl))*nprocsinl  ! continue on first procinl again
          else
            if (j == jgbin) then
              jgbeg = j-1
              jbeg = jein
              procinlo = nprocsinl-1
              filenumstart = -1
              jj = j+jtotdum-2
              procinup = floor(real(jj-jgbin)/real(jtotin))
              filenumend = procinup
              jend = jj-(procinup*jtotin)
              procinup = procinup-floor(real(procinup)/real(nprocsinl))*nprocsinl  !continue on first procinl again
            else
              procinlo = floor(real(j-jgbin-1)/real(jtotin))    ! One cell lower is needed
              filenumstart = procinlo
              jgbeg = j-1
              jbeg = j-(procinlo*jtotin)-1
              jj = j+jtotdum-2
              procinup = floor(real(jj-jgbin)/real(jtotin))
              filenumend = procinup
              jend = jj-(procinup*jtotin)
              procinup = procinup-floor(real(procinup)/real(nprocsinl))*nprocsinl  ! continue on first procinl again
            end if ! j=jgbin
          end if
          exit
        end if
      end do

      write(6,*) '!! myid,procinlo,jbeg,procinup,jend,jgbeg = ',myid,procinlo,jbeg,procinup,jend,jgbeg


   ! make dummy y-grid (equidistant)
      do j = jbdum,jedum+1
        yhdum(j) = yhin(jgbeg) + (j-jbdum+1)*dyin
        yfdum(j) = yhdum(j) + 0.5*dyin
      end do
!      if (procoldup /= procoldlo) then
!        write(6,*) '!!! Start-cell and end-cell are not in the same file!!!'
!      end if
      filestoread = filenumend - filenumstart + 1 ! no. of files to be read
!      write(6,*) '!! procinlo,procinup = ',procinlo,procinup
!      write(6,*) '!! jbin,jein,jbeg,jend= ',jbin,jein,jbeg,jend

   ! for components defined on yf
    do j=jb,je
      do jj=jbdum+1,jedum+1
        if (yfdum(jj) >= yf(j)) then
          ylocupf(j)  = jj
          yloclowf(j) = jj-1
          exit
        end if
      end do
    end do
   ! for components defined on yh
    do j=jb,je+1
      do jj=jbdum+1,jedum+1
        if (yhdum(jj) >= yh(j)) then
          ylocuph(j)  = jj
          yloclowh(j) = jj-1
          exit
        end if
      end do
    end do

  end subroutine readzincoord