inletgennotemp Subroutine

  subroutine inletgennotemp
    use modglobal,   only : ib,ie,jb,je,jb,jtot,kb,ke,zf,zh,dzf,dzhi,timee,btime,totavtime,rk3step,dt,numol,iplane,lles,iinletgen,inletav,runavtime,Uinf,lwallfunc,linletRA,totinletav,lstoreplane,nstore,lfixinlet,lfixutauin,luvolflowr
    use modfields,   only : u0,v0,w0,wm,uprof
    use modsave,     only : writerestartfiles
    use modmpi,      only : slabsum,myid
    implicit none

    real,dimension(ib:ib,jb:je,kb:ke)   :: uinletbc2   ! dummy variable
    real,dimension(jb:je,kb:ke)   :: uprec             ! velocity fluctuation (up_rec = u0 - Urec)
    real,dimension(jb:je,kb:ke)   :: vprec             ! velocity fluctuation (vp_rec = v0 - 0)
    real,dimension(jb:je,kb:ke+1) :: wprec             ! velocity fluctuation (wp_rec = w0 - Wrec)
    real,dimension(jb:je,kb:ke)   :: upinli,vpinli     ! = gamma * (uprec,v interpolated to zii grid)
    real,dimension(jb:je,kb:ke)   :: upinlo,vpinlo     ! = gamma * (uprec,v interpolated to zoi grid)
    real,dimension(jb:je,kb:ke+1) :: wpinli            ! = gamma * (wprec   interpolated to zii grid)
    real,dimension(jb:je,kb:ke+1) :: wpinlo            ! = gamma * (wprec   interpolated to zoi grid)
    real,dimension(kb:ke)   :: udiff                   ! difference between Uinl and Urec
!    real,dimension(kb:ke)   :: Urecdiff                ! difference between Urec new and old
    real,dimension(kb:ke)   :: urav                    ! j-averaged u-velocity (not time-averaged)
    real,dimension(kb:ke)   :: uravdzf                 ! j-averaged u-velocity (not time-averaged) times dzf
    real,dimension(kb:ke)   :: uinldzf                 ! j-averaged u-velocity (not time-averaged) times dzf
    real,dimension(kb:ke)   :: Urecdzf                 ! Urec times dzf
    real,dimension(kb:ke+1) :: wrav                    ! j-averaged w-velocity (not time-averaged)
    real,dimension(kb:ke)   :: Uinli                   ! = gamma * (Urec interpolated to ziif grid points)
    real,dimension(kb:ke+1) :: Winli                   ! = gamma * (Wrec interpolated to ziih grid points)
    real,dimension(kb:ke)   :: Uinlo                   ! = gamma * (Urec interpolated to zioif grid points)
    real,dimension(kb:ke+1) :: Winlo                   ! = gamma * (Wrec interpolated to zoih grid points)
    real,dimension(kb:ke)   :: wfuncf                  ! weight function at full level
    real,dimension(kb:ke+1) :: wfunch                  ! weight function at half level
    real                    :: utaur2,utaui2           ! (utau)^2 at recycle station and inlet
    real                    :: gamm                    ! utaui / utaur
    real                    :: avint,avinti            ! avering interval
    real                    :: alpha,beta              ! factors used in the Weight function
!    real                    :: totalu                  ! total u-velocity at outlet
    real                    :: Urectot                  ! total u-velocity at recycle plane
    real                    :: rk3coef
!    real                    :: di_test                 ! BL thickness as measured from Uinl
    real                    :: utop                    ! j-averaged top velocity
    real                    :: interval
    real                    :: dtinrk                  ! RK time step in inlet data
    real                    :: rk3coefin               ! Cumulative RK time step in inlet data
    real                    :: dr_old
    real                    :: scalef                      ! scale factor to scale instantaneous velocity profile with to get constant mass flux
    real                    :: totaluinl                   ! bulk velocity at the inlet
    integer i,j,k,kk

   if (iinletgen == 1) then

   u0inletbcold = u0inletbc
   v0inletbcold = v0inletbc
   w0inletbcold = w0inletbc
   totaluold    = totalu
   displold     = displ
   ddispdxold   = ddispdx

   ! compute time-average velocities
    rk3coef = dt / (4. - dble(rk3step))
    if (rk3step==1) then
      deltat = rk3coef
    elseif (rk3step==2) then
      deltat = rk3coef - (dt/3.)
    elseif (rk3step==3) then
      deltat = rk3coef - (dt/2.)
    end if

    if (linletRA) then  ! this is a switch to use 'running average'
      avint = totinletav + timee-btime ! runav interval = averaging interval previuous sim  + current elapsed sim time
    else
      avint  = inletav
    end if
    avinti = 1./avint
    uaver=0.
    do i=ib,ie
      call slabsum(uaver(i,:),kb,ke,u0(i:i,jb:je,kb:ke),i,i,jb,je,kb,ke,i,i,jb,je,kb,ke)
    end do

    wrav=0.
    call slabsum(wrav(kb:ke+1),kb,ke,w0(irecy-1:irecy-1,jb:je,kb:ke+1),irecy-1,irecy-1,jb,je,kb,ke+1,irecy-1,irecy-1,jb,je,kb,ke+1)

    uaver = uaver / jtot                ! average over j-direction
    urav = uaver(irecy,:)
    wrav = wrav / jtot                 ! average over j-direction

    do k=kb,ke
      Urec(k) =  urav(k)*deltat*avinti + (1.-deltat*avinti)*Urec(k)
    end do
    do k=kb,ke+1
      Wrec(k) =  wrav(k)*deltat*avinti + (1.-deltat*avinti)*Wrec(k)
    end do
    do k=kb,ke
    do i=ib,ie
      Utav(i,k) =  uaver(i,k)*deltat*avinti + (1.-deltat*avinti)*Utav(i,k)
    end do
    end do



   ! compute velocity fluctuation at recycle station
    do k=kb,ke
      do j=jb,je
        uprec(j,k) = u0(irecy,j,k) - Urec(k)
        vprec(j,k) = v0(irecy-1,j,k)             ! mean v is zero
      end do
    end do

    do k=kb,ke+1
      do j=jb,je
        wprec(j,k) = w0(irecy-1,j,k) - Wrec(k)   ! note that w-velocity is taken at i=irecy-1 !!
      end do
    end do


    if (lwallfunc) then
      call wallawinlet(Urec(kb),dzf(kb),numol,utaur2)    ! compute wall shear stress at recycle station
    else
      utaur2 = 2.*numol*Urec(kb)/dzf(kb)
    end if
    utaur = sqrt(abs(utaur2))                          ! compute utau at recycle station
   ! compute momentum thickness at inlet and recycle plane

   dr_old = dr
!   call blthickness(dr,utaur)                     ! also needed for thetar
   call blthicknesst(dr,Urec,0.99)
!   call momentumthickness(thetai,utaui,di)        ! di is kept fixed
   call momentumthicknessexp(thetai,Uinl)
!   call momentumthickness(thetar,utaur,dr)
   call momentumthicknessexp(thetar,Urec)
!   call blthickness(dr,utaur)

    if (.not.lfixutauin) then
      utaui  = utaur* abs(thetar/thetai)**(1./8.)   ! See Lund (1998): 'Similar to Ludwig-Tillmann correlation'
    end if
    gamm = utaui / utaur                          ! Gamma in Lund (1998)

   ! compute inner scaling coordinates
    zirf = utaur*zf / numol                       ! inner scaling zf-coordinate at recycle station
    zirh = utaur*zh / numol                       ! inner scaling zh-coordinate at recycle station
    ziif = utaui*zf / numol                       ! inner scaling zf-coordinate at inlet station
    ziih = utaui*zh / numol                       ! inner scaling zh-coordinate at inlet station

   ! compute outer scaling coordinates
    zorf = zf / dr                                ! outer scaling zf-coor as measured from Uinldinate at recycle station
    zorh = zh / dr                                ! outer scaling zh-coordinate at recycle station
    zoif = zf / di                                ! outer scaling zf-coordinate at inlet station  (could be done once, actually..)
    zoih = zh / di                                ! outer scaling zf-coordinate at inlet station  (could be done once, actually..)

!!!!! Interpolation starts here

 !!! First inner coordinates
   ! determine which elements are needed when recycle velocity profile is interpolated on inlet plane
   ! for u(,v)-components (zf)
    do k=kb,ke
      do kk=kb,ke
        if (zirf(kk) >= ziif(k)) then
          locupif(k)  = kk
          loclowif(k) = kk-1
          exit
        elseif (kk==ke) then
          locupif(k)  = ke+1            ! this means extrapolation!
          loclowif(k) = ke-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=kb,ke+1
        if (zirh(kk) >= ziih(k)) then
          locupih(k)  = kk
          loclowih(k) = kk-1
          exit
        elseif (kk==ke+1) then
          locupih(k)  = ke+2            ! this means extrapolation!
          loclowih(k) = ke
        end if
      end do
    end do
 !!! Finished with inner coordinates

 !!! Do the same trick for outer coordinates
   ! determine which elements are needed when recycle velocity profile is interpolated on inlet plane
   ! for u(,v)-components (zf)
    do k=kb,ke
      do kk=kb,ke
        if (zorf(kk) >= zoif(k)) then
          locupof(k)  = kk
          loclowof(k) = kk-1
          exit
        elseif (kk==ke) then
          locupof(k)  = ke+1            ! this means extrapolation!
          loclowof(k) = ke-1
        end if
      end do
    end do
   ! for w-components (zh)
    do k=kb,ke+1
      do kk=kb,ke+1
        if (zorh(kk) >= zoih(k)) then
          locupoh(k)  = kk
          loclowoh(k) = kk-1
          exit
        elseif (kk==ke+1) then
          locupoh(k)  = ke+2            ! this means extrapolation!
          loclowoh(k) = ke
        end if
      end do
    end do

  !!! Finished with outer coordinates

!!! Now really interpolate
  !!! First inner coordinates
   ! compute Urec on zii grid
    do k=kb,ke
      if (locupif(k) == ke+1) then      ! indicator for extrapolation!
!        Uinli(k) = Urec(ke) + (Urec(ke) - Urec(ke-1)) / (zirf(ke)-zirf(ke-1)) * (ziif(k)-zirf(ke))
        Uinli(k) = Urec(ke)
      elseif (loclowif(k) == kb-1) then ! interprets this as extrapolation to bottom (use u=0 at z+=0)
        Uinli(k) = Urec(kb)/zirf(kb) * ziif(k)
      else                            ! normal interpolation
        Uinli(k) = Urec(loclowif(k)) + (Urec(locupif(k)) - Urec(loclowif(k))) / (zirf(locupif(k)) - zirf(loclowif(k)))* (ziif(k) - zirf(loclowif(k)))
        if ((ziif(k) .gt. zirf(locupif(k))) .or. (ziif(k) .lt. zirf(loclowif(k)))) then
          write(6,*) '!!!Mistake in Interpolation !!!!'
        end if
      end if
    end do

   ! compute Wrec on zii grid
    Winli(kb) = 0.0                      ! corresponds to ground level
    do k=kb+1,ke+1
      if (locupih(k) == ke+2) then     ! indicator for extrapolation!
!        Winli(k) = Wrec(ke+1) + (Wrec(ke+1) - Wrec(ke)) / (zirh(ke+1)-zirh(ke)) * (ziih(k)-zirh(ke+1))
        Winli(k) = Wrec(ke+1)
      else                            ! normal interpolation
        Winli(k) = Wrec(loclowih(k)) + (Wrec(locupih(k)) - Wrec(loclowih(k))) / (zirh(locupih(k)) - zirh(loclowih(k))) * (ziih(k) - zirh(loclowih(k)))
      end if
    end do

   ! compute u- and v- and t-fluctuation on zii grid
    do k=kb,ke
      if (locupif(k) == ke+1) then      ! indicator for extrapolation!
!        upinli(:,k) = uprec(:,ke) + (uprec(:,ke) - uprec(:,ke-1)) / (zirf(ke)-zirf(ke-1)) * (ziif(k)-zirf(ke))
        upinli(:,k) = 0.
        vpinli(:,k) = 0.
      elseif (loclowif(k) == kb-1) then ! interprets this as extrapolation to bottom (use u=0 at z+=0)
        upinli(:,k) = uprec(:,kb)/zirf(kb) * ziif(k)
        vpinli(:,k) = vprec(:,kb)/zirf(kb) * ziif(k)
      else                            ! normal interpolation
        upinli(:,k) = uprec(:,loclowif(k)) + (uprec(:,locupif(k)) - uprec(:,loclowif(k))) / (zirf(locupif(k)) - zirf(loclowif(k))) * (ziif(k) - zirf(loclowif(k)))
        vpinli(:,k) = vprec(:,loclowif(k)) + (vprec(:,locupif(k)) - vprec(:,loclowif(k))) / (zirf(locupif(k)) - zirf(loclowif(k))) * (ziif(k) - zirf(loclowif(k)))
      end if
    end do

   ! compute w-fluctuation on zii grid
    do k=kb+1,ke+1
!      if (locupih(k) == ke+1) then      ! indicator for extrapolation!
      if (locupih(k) == ke+2) then      ! indicator for extrapolation!
!        wpinli(:,k) = wprec(:,ke+1) + (wprec(:,ke+1) - wprec(:,ke)) / (zirh(ke+1)-zirh(ke)) * (ziih(k)-zirh(ke+1))
        wpinli(:,k) = 0.
      else                            ! normal interpolation
        wpinli(:,k) = wprec(:,loclowih(k)) + (wprec(:,locupih(k)) - wprec(:,loclowih(k))) / (zirh(locupih(k)) - zirh(loclowih(k))) * (ziih(k) - zirh(loclowih(k)))
      end if
    end do

 !! Finished with interpolating inner variables
 !! Continue with interpolating outer variables
   ! compute Urec on zoi grid
    do k=kb,ke
      if (locupof(k) == ke+1) then      ! indicator for extrapolation!
!        Uinlo(k) = Urec(ke) + (Urec(ke) - Urec(ke-1)) / (zorf(ke)-zorf(ke-1)) * (zoif(k)-zorf(ke))
!        Uinlo(k) = Urec(ke)
        Uinlo(k) = Uinf
      elseif (loclowof(k) == kb-1) then ! interprets this as extrapolation to bottom (use u=0 at z+=0)
        Uinlo(k) = Urec(kb)/zorf(kb) * zoif(k)
      else                            ! normal interpolation
        Uinlo(k) = Urec(loclowof(k)) + (Urec(locupof(k)) - Urec(loclowof(k))) / (zorf(locupof(k)) - zorf(loclowof(k))) * (zoif(k) - zorf(loclowof(k)))
      end if
    end do

   ! compute Wrec on zii grid
    Winlo(kb) = 0.0                      ! corresponds to ground level
    do k=kb+1,ke+1
      if (locupoh(k) == ke+2) then     ! indicator for extrapolation!
!        Winlo(k) = Wrec(ke+1) + (Wrec(ke+1) - Wrec(ke)) / (zorh(ke+1)-zorh(ke)) * (zoih(k)-zorh(ke+1))
        Winlo(k) = Wrec(ke+1)
      else                            ! normal interpolation
        Winlo(k) = Wrec(loclowoh(k)) + (Wrec(locupoh(k)) - Wrec(loclowoh(k))) / (zorh(locupoh(k)) - zorh(loclowoh(k))) * (zoih(k) - zorh(loclowoh(k)))
      end if
    end do

   ! compute u- and v-fluctuation on zoi grid
    do k=kb,ke
      if (locupof(k) == ke+1) then      ! indicator for extrapolation!
!        upinlo(:,k) = uprec(:,ke) + (uprec(:,ke) - uprec(:,ke-1)) / (zorf(ke)-zorf(ke-1)) * (zoif(k)-zorf(ke))
        upinlo(:,k) = 0.
        vpinlo(:,k) = 0.
      elseif (loclowof(k) == kb-1) then ! interprets this as extrapolation to bottom (use u=0 at z+=0)
        upinlo(:,k) = uprec(:,kb)/zorf(kb) * zoif(k)
        vpinlo(:,k) = vprec(:,kb)/zorf(kb) * zoif(k)
      else                            ! normal interpolation
        upinlo(:,k) = uprec(:,loclowof(k)) + (uprec(:,locupof(k)) - uprec(:,loclowof(k))) / (zorf(locupof(k)) - zorf(loclowof(k))) * (zoif(k) - zorf(loclowof(k)))
        vpinlo(:,k) = vprec(:,loclowof(k)) + (vprec(:,locupof(k)) - vprec(:,loclowof(k))) / (zorf(locupof(k)) - zorf(loclowof(k))) * (zoif(k) - zorf(loclowof(k)))
      end if
    end do

   ! compute w-fluctuation on zoi grid
    do k=kb+1,ke+1
      if (locupoh(k) == ke+2) then      ! indicator for extrapolation!
!        wpinlo(:,k) = wprec(:,ke+1) + (wprec(:,ke+1) - wprec(:,ke)) / (zorh(ke+1)-zorh(ke)) * (zoih(k)-zorh(ke+1))
        wpinlo(:,k) = 0.
      else                            ! normal interpolation
        wpinlo(:,k) = wprec(:,loclowoh(k)) + (wprec(:,locupoh(k)) - wprec(:,loclowoh(k))) / (zorh(locupoh(k)) - zorh(loclowoh(k))) * (zoih(k) - zorh(loclowoh(k)))
      end if
    end do
 !! Finished interpolating outer velocity variables
!!!!! Finished Interpolation! !!!!!


   ! compute rescaled inner variables ! Winli = Winli (interpolation is enough)
    Uinli = gamm* Uinli
    upinli = gamm* upinli
    vpinli = gamm* vpinli
    wpinli = gamm* wpinli
   ! compute rescaled outer variables ! Winlo = Winlo (interpolation is enough)
    Uinlo = gamm* Uinlo  + (1.- gamm)*Uinf
!    Uinlo = gamm* Uinlo  + (1.- gamm)*Urec(ke)
    upinlo = gamm* upinlo
    vpinlo = gamm* vpinlo
    wpinlo = gamm* wpinlo


   ! Compute weight function (alpha=4, b=0.2)
    alpha = 4.
    beta = 0.2
    wfuncf = 0.5*(1. + tanh( alpha*(zoif-beta)/((1.-2.*beta)*zoif+beta) )/tanh(alpha) ) ! for full level height
    wfunch = 0.5*(1. + tanh( alpha*(zoih-beta)/((1.-2.*beta)*zoih+beta) )/tanh(alpha) ) ! for half level height
    do k=kb,ke
      if (wfuncf(k) .gt. 1.) then
        wfuncf(k) = 1.
      end if
    end do
    do k=kb,ke+1
      if (wfunch(k) .gt. 1.) then
        wfunch(k) = 1.
      end if
    end do


   ! Compute the velocity components for the inlet BC
    do k=kb,ke
    do j=jb,je
      u0inletbc(j,k) = (Uinli(k)+ upinli(j,k))*(1.-wfuncf(k)) +  (Uinlo(k) + upinlo(j,k))* wfuncf(k)
      v0inletbc(j,k) =            vpinli(j,k) *(1.-wfuncf(k)) +              vpinlo(j,k) * wfuncf(k)
    end do
    end do


    do k=kb,ke+1
    do j=jb,je
      w0inletbc(j,k) = (Winli(k)+ wpinli(j,k))*(1-wfunch(k)) +  (Winlo(k) + wpinlo(j,k))* wfunch(k)
    end do
    end do
    w0inletbc(:,kb) = 0.
    w0inletbc(:,ke+1) = 0.


   ! Compute j-averaged inlet U  (used for compute thetai)
    uinletbc2(ib,jb:je,kb:ke) = u0inletbc(jb:je,kb:ke)  ! this is just a dummy variable to give uninletbc the right dimension in slabsum
    urav = 0.
    call slabsum(urav  ,kb,ke,uinletbc2 ,ib,ib,jb,je,kb,ke,ib,ib,jb,je,kb,ke)
    urav = urav / jtot                   ! average over j-direction

! determine bulk velocity of new profile
      do k=kb,ke
        uravdzf(k) = urav(k)*dzf(k)
      end do
      totalu = sum(uravdzf(kb:ke))/(zh(ke+1)-zh(kb))      ! Area-averaged outflow velocity

! correct instantaneous inflow velocity for constant mass-flux
      if (luvolflowr ) then
        do k=kb,ke
          uinldzf(k) = Uinl(k)*dzf(k)
        end do
        totaluinl = sum(uinldzf(kb:ke))/(zh(ke+1)-zh(kb))      ! Area-averaged inflow velocity that should be kept
        scalef = totaluinl/totalu ! compute factor to scale the velocity profile with
        u0inletbc(:,:) = u0inletbc(:,:)*scalef  ! rescale the velocity profile to have constant mass-flux
        urav(:) = urav(:)*scalef                ! also rescale the part that is added to the mean
      end if

   ! Compute j- and time-averaged  inlet U  (used for compute thetai)
    if (.not.lfixinlet) then  ! only update the average inlet profiles when lfixinlet .eqv..false.
      do k=kb,ke
        Uinl(k) =  urav(k)*deltat*avinti + (1.-deltat*avinti)*Uinl(k)
      end do
    end if

! write inletplane to array (and to file after 1000 time steps)
    if (lstoreplane ) then
      storeu0inletbc(:,:,nstepread) = u0inletbc(:,:)
      storev0inletbc(:,:,nstepread) = v0inletbc(:,:)
      storew0inletbc(:,:,nstepread) = w0inletbc(:,:)
      nstepread = nstepread +1
      if (nstepread == nstore+1) then
        nfile = nfile +1       ! next file number
        call writeinletfile    ! write 1000 time steps to file
        call writerestartfiles
        nstepread = 1          ! reset counter
      end if   ! nstepread == 1001
    end if ! lstoreplane


    if (rk3step==1) then
      uminletbc = u0inletbc
      vminletbc = v0inletbc
      wminletbc = w0inletbc
    end if

!   call blthickness(di_test,utaui)
   call blthicknesst(di_test,Uinl,0.99)

!   call dispthickness(displ)  ! needed in top BC
   call dispthicknessexp(displ)  ! needed in top BC

   if ((myid==0) .and. (rk3step==3)) then
     write(6,*) 'Inlet Gen: gamma=',gamm
     write(6,*) 'Inlet Gen: Uinl(ke)=',Uinl(ke)
     write(6,*) 'Inlet Gen: utaui,utaur =',utaui,utaur
     write(6,*) 'Inlet Gen: deltar, deltai_test', dr, di_test
     write(6,*) 'Inlet Gen: d*i, d*r=',displ(ib),displ(irecy)
     write(6,*) 'Inlet Gen: thetai,thetar',thetai,thetar
     if (luvolflowr ) then
       write(6,*) 'Inlet Gen: mass flux correction factor = ',scalef
!       write(6,*) 'Inlet Gen: mass flux                   = ',totalreadu
       write(6,*) 'Inlet Gen: mass flux                   = ',totaluinl
     end if
   end if

    elseif (iinletgen == 2) then
      if (myid==0) then
        write(6,*) 'nstepread=',nstepread
      end if
      u0inletbcold = u0inletbc
      v0inletbcold = v0inletbc
      w0inletbcold = w0inletbc

  ! determine time step interval in simulation
      rk3coef   = dt / (4. - dble(rk3step))
      if (rk3step==1) then
        deltat = rk3coef
      elseif (rk3step==2) then
        deltat = rk3coef   - (dt/3.)
      elseif (rk3step==3) then
        deltat = rk3coef - (dt/2.)
      end if
  ! determine time step interval in inlet data
      rk3coefin = dtin / (4. - dble(rk3stepin))
      if (rk3stepin==1) then
        dtinrk = rk3coefin
      elseif (rk3stepin==2) then
        dtinrk = rk3coefin - (dtin/3.)
      elseif (rk3stepin==3) then
        dtinrk = rk3coefin - (dtin/2.)
      end if

      interval = dtinrk - elapstep
      elapstep = elapstep + deltat
      if (elapstep > dtinrk) then      ! use new value at next time step
        nstepread = nstepread +1
        elapstep = mod(elapstep,dtinrk)
        rk3stepin = mod(rk3stepin,3) + 1
        rk3coefin = dtin / (4. - dble(rk3stepin))
        if (rk3stepin==1) then
          dtinrk = rk3coefin
        elseif (rk3stepin==2) then
          dtinrk = rk3coefin - (dtin/3.)
        elseif (rk3stepin==3) then
          dtinrk = rk3coefin - (dtin/2.)
        end if
        u0inletbc(:,:) = storeu0inletbc(:,:,nstepread)
        v0inletbc(:,:) = storev0inletbc(:,:,nstepread)
        w0inletbc(:,:) = storew0inletbc(:,:,nstepread)
        if (nstepread == nstore) then
          nfile = nfile + 1
          call readinletfile
          call writerestartfiles
          nstepread = 0
        end if
        interval = dtinrk
        deltat = elapstep
!        write(6,*) 'dtinrk,deltat=', dtinrk,deltat
      end if
      u0inletbc(:,:) = (1. - deltat/interval)*u0inletbc(:,:) + (deltat/interval)*storeu0inletbc(:,:,nstepread+1)
      v0inletbc(:,:) = (1. - deltat/interval)*v0inletbc(:,:) + (deltat/interval)*storev0inletbc(:,:,nstepread+1)
      w0inletbc(:,:) = (1. - deltat/interval)*w0inletbc(:,:) + (deltat/interval)*storew0inletbc(:,:,nstepread+1)


!! massflow correction
      uinletbc2(ib,jb:je,kb:ke) = u0inletbc(jb:je,kb:ke)  ! this is just a dummy variable to give uninletbc the right dimension in slabsum
      urav = 0.
      call slabsum(urav  ,kb,ke,uinletbc2 ,ib,ib,jb,je,kb,ke,ib,ib,jb,je,kb,ke)
      urav = urav / jtot                    ! average over j-direction

   ! determine bulk velocity of new (interpolated) profile
      do k=kb,ke
        uravdzf(k) = urav(k)*dzf(k)
      end do
      totalu    = sum(uravdzf(kb:ke))/(zh(ke+1)-zh(kb))      ! Area-averaged outflow velocity

   ! rescale the instantaneous profile to keep mass flux constant (tot avoid pressure fluctuations)
      scalef = totalreadu/totalu ! compute factor to scale the velocity profile with
      u0inletbc(:,:) = u0inletbc(:,:)*scalef  ! rescale the velocity profile to have constant mass-flux
!! end of massflow correction of interpolated streamwise velocity



      if (rk3step==1) then
        uminletbc = u0inletbc
        vminletbc = v0inletbc
        wminletbc = w0inletbc
      end if
    end if  ! iinletgen

  end subroutine inletgennotemp