subroutine readzincoord
use modglobal, only : kb,ke,kh,ifinput,zf,zh,ysize,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 = ysize/nprocs
dyin = ysize / 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*(ysize/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