readinletfile Subroutine

  subroutine readinletfile
    use modglobal, only : ib,jb,je,jmax,kb,ke,cexpnr,ifinput,nstore,ltempeq,ntrun,zh,jtot,jh
    use modmpi,    only : cmyid,myid,nprocs,slabsum,excjs
!    use modinletdata, only : storeu0inletbc,storev0inletbc,storew0inletbc,nfile

    implicit none
    real, dimension(ib:ib,jb:jb+inlfactor*jmax-1,kbin:kein)        :: udummy
    real, dimension(kbin:kein)                                     :: uread
    real, dimension(kbin:kein)                                     :: ureaddzfin
    real, dimension(jb:jb+jtotin-1,kbin:kein,1:nstore)     :: storeu0inold
    real, dimension(jb:jb+jtotin-1,kbin:kein,1:nstore)     :: storev0inold
    real, dimension(jb:jb+jtotin-1,kbin:kein+1,1:nstore)   :: storew0inold
    real, dimension(jb:jb+jtotin-1,kbin:kein,1:nstore)     :: storet0inold
    real, dimension(jb:jb+jtotdum-1,kbin:kein,1:nstore)     :: storeu0indum
    real, dimension(jb:jb+jtotdum-1,kbin:kein,1:nstore)     :: storev0indum
    real, dimension(jb:jb+jtotdum-1,kbin:kein+1,1:nstore)   :: storew0indum
    real, dimension(jb:jb+jtotdum-1,kbin:kein,1:nstore)     :: storet0indum
    real, dimension(jb:je,kbin:kein,1:nstore)     :: storeu0innew
    real, dimension(jb:je,kbin:kein,1:nstore)     :: storev0innew
    real, dimension(jb:je,kbin:kein+1,1:nstore)   :: storew0innew
    real, dimension(jb:je,kbin:kein,1:nstore)     :: storet0innew
    integer filen,filee
    integer fileid
    integer j,k,n,js,jf,jfdum,jsdum
    character(24) name
      filen = -5   ! should be checked carefully if in use ever
      jfdum = jbdum-1  ! initial value
      do fileid = filenumstart, filenumstart+(filestoread-1)
        if (filen == -1) then
          filen = nprocsinl-1                      ! -1 means the last proc (periodic)
        else
          filen = fileid - floor(real(fileid)/real(nprocsinl))*nprocsinl  ! loop over proc's
        end if
       write(6,*) '!!!!! filen = ', filen
        name = 'inlet/inlet_    k   .'
        write (name(13:16)  ,'(i4.4)') nfile
        write (name(18:20)  ,'(i3.3)') filen
        name(22:24)= cexpnr
        write(6,*) 'Reading Inlet velocity: ', name
        open(unit=ifinput,file=name,form='unformatted')
        read(ifinput)  (((storeu0inold (j,k,n),j=jbin,jein),k=kbin,kein),  n=1,nstore)
        read(ifinput)  (((storev0inold (j,k,n),j=jbin,jein),k=kbin,kein),  n=1,nstore)
        read(ifinput)  (((storew0inold (j,k,n),j=jbin,jein),k=kbin,kein+1),n=1,nstore)
        close (ifinput)

        if (ltempeq ) then
          name = 'inlet/itemp_    k   .'
          write (name(13:16)  ,'(i4.4)') nfile
          write (name(18:20)  ,'(i3.3)') filen
          name(22:24)= cexpnr
          write(6,*) 'Reading Inlet temperature: ', name
          open(unit=ifinput,file=name,form='unformatted')
          read(ifinput)  (((storet0inold (j,k,n),j=jbin,jein),k=kbin,kein),  n=1,nstore)
          close (ifinput)
        end if


        ! determine start and end indices
        if (filen == procinlo) then
          js = jbeg
        else
          js = jbin
        end if
        if (filen == procinup) then
          jf = jend
        else
          jf = jein
        end if
        jsdum = jfdum + 1
        jfdum = jsdum + (jf-js)
!        if (jsdum >= 3) write(6,*) 'myid, jsdum = ',myid, jsdum
!        if (jfdum >= 3) write(6,*) 'myid, jfdum = ',myid, jfdum

 !!! put values from original in dummy variable
        storeu0indum(jsdum:jfdum,:,:)    = storeu0inold(js:jf,:,:)      ! s: start  f: final
        storev0indum(jsdum:jfdum,:,:)    = storev0inold(js:jf,:,:)      ! s: start  f: final
        storew0indum(jsdum:jfdum,:,:)    = storew0inold(js:jf,:,:)      ! s: start  f: final
        if (ltempeq ) then
          storet0indum(jsdum:jfdum,:,:)    = storet0inold(js:jf,:,:)      ! s: start  f: final
        end if
      end do  ! loop over original inlet files
  ! now interpolate in y
      call yinterpolate (storeu0indum,storeu0innew,kbin,kein)
      call yinterpolate (storev0indum,storev0innew,kbin,kein)
      call yinterpolate (storew0indum,storew0innew,kbin,kein+1)
      call yinterpolate (storet0indum,storet0innew,kbin,kein)

      if (.not.lzinzsim) then ! interpolate when zin =/ zsim
        call zinterpolate (storeu0innew(:,:,:),storeu0inletbc)   ! interpolate inlet profile to zgrid
        call zinterpolate (storev0innew(:,:,:),storev0inletbc)   ! interpolate inlet profile to zgrid
        call zinterpolatew(storew0innew(:,:,:),storew0inletbc)   ! interpolate inlet profile to zgrid
        if (ltempeq) then
          call zinterpolatet(storet0innew(:,:,:),storet0inletbc)   ! interpolate inlet profile to zgrid
        end if
      else
        storeu0inletbc(:,:,:) = storeu0inold(:,:,:)
        storev0inletbc(:,:,:) = storev0inold(:,:,:)
        storew0inletbc(:,:,:) = storew0inold(:,:,:)
        if (ltempeq) then
          storet0inletbc(:,:,:) = storet0inold(:,:,:)
        end if
      end if

      if (iangle/=0.0) then   ! modify for inflow angle
         do n=1,nstore
         do k=kb,ke
         do j=jb,je
           u0rot(n,j,k) = storeu0inletbc(j,k,n)         ! swap indices in order
           v0rot(n,j,k) = storev0inletbc(j,k,n)         ! to use excjs
         end do
         end do
         end do
         call excjs(u0rot, 1,nstore,jb,je,kb,ke,0,jh)
         call excjs(v0rot, 1,nstore,jb,je,kb,ke,0,jh)
!         write(6,*) 'v0rot(1,je+1,30) = ',v0rot(1,je+1,30)
         do n=1,nstore
         do k=kb,ke
         do j=jb,je
         ! apply horizontal rotation (neglecting the delta_x difference)
           storeu0inletbc(j,k,n) = u0rot(n,j,k)*cos(iangle) - 0.5*sin(iangle)*(v0rot(n,j,k)+v0rot(n,j+1,k))
           storev0inletbc(j,k,n) = v0rot(n,j,k)*cos(iangle) + 0.5*sin(iangle)*(u0rot(n,j,k)+u0rot(n,j-1,k))
         end do
         end do
         end do
      end if ! iangle =/0.0


  end subroutine readinletfile