inletgen Subroutine

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

    real,dimension(ib:ib,jb:je,kb:ke)   :: uinletbc2   ! dummy variable
    real,dimension(ib:ib,jb:je,kb:ke)   :: tinletbc2   ! 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)   :: tprec             ! temperature fluctuation (tp_rec = t0 - Trec)
    real,dimension(jb:je,kb:ke)   :: upinli,vpinli     ! = gamma * (uprec,v interpolated to zii grid)
    real,dimension(jb:je,kb:ke)   :: tpinli            ! = lambda  * (tprec   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)   :: tpinlo            ! = lambda  * (tprec   interpolated to zoti 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)   :: trav                    ! j-averaged temperature (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)   :: Tinli                   ! = lambda  * (Trec interpolated to ziif 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)   :: Tinlo                   ! = lambda  * (Trec interpolated to zoti grid points)
    real,dimension(kb:ke)   :: wfuncf                  ! weight function at full level
    real,dimension(kb:ke+1) :: wfunch                  ! weight function at half level
    real,dimension(kb:ke)   :: wfunct                  ! weight function at full level
    real                    :: utaur2,utaui2           ! (utau)^2 at recycle station and inlet
    real                    :: gamm                    ! utaui / utaur
    real                    :: lamb                    ! ttaui / ttaur
    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
!    real                    :: q0                      ! heat flux
    integer i,j,k,kk,kdamp


   if (iinletgen == 1) then

   u0inletbcold = u0inletbc
   v0inletbcold = v0inletbc
   w0inletbcold = w0inletbc
   t0inletbcold = t0inletbc    ! temperature
   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.
    taver=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)
      call slabsum(taver(i,:),kb,ke,thl0(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)
    trav=0.
    call slabsum(trav(kb:ke),  kb,ke,thl0(irecy-1:irecy-1,jb:je,kb:ke), irecy-1,irecy-1,jb,je,kb,ke,  irecy-1,irecy-1,jb,je,kb,ke)

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

    do k=kb,ke
      Urec(k) =  urav(k)*deltat*avinti + (1.-deltat*avinti)*Urec(k)
      Trec(k) =  trav(k)*deltat*avinti + (1.-deltat*avinti)*Trec(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)
      Ttav(i,k) =  taver(i,k)*deltat*avinti + (1.-deltat*avinti)*Ttav(i,k)
    end do
    end do

!    Urec = Urec +(Uinf-Urec(ke))     ! make sure at the recycle plane the top velocity equals Uinf


!    Urecdiff = Urecdiff - Urec
!    if (myid==0) then
!      write(6,*) 'Urec_old - Urec_new (kb+40)=',Urecdiff(kb+40)
!    end if

!! check if Urec contains NaN
!    if (myid==0) then
!      write(6,*) 'Checking Urec for NaN'
!      do k=kb,ke
!        if (ISNAN(Urec(k))) then
!          write(6,*) 'Urec(k)=NaN at k=kb+', k-kb
!        end if
!      end do
!      write(6,*) 'Finished checking Urec for NaN'
!    end if


!    if (myid==0) then
!      write(6,*) 'myid, Urec(ke)=',myid, Urec(ke)
!      write(6,*) 'wrav(ke), Wrec(ke)=',wrav(ke), Wrec(ke)
!      write(6,*) 'wrav(ke-1), Wrec(ke-1)=',wrav(ke-1), Wrec(ke-1)
!      write(6,*) 'wrav(ke-10), Wrec(ke-10)=',wrav(ke-10), Wrec(ke-10)
!      write(6,*) 'wrav(ke-30), Wrec(ke-30)=',wrav(ke-30), Wrec(ke-30)
!      write(6,*) 'wrav(kb+10), Wrec(kb+10)=',wrav(kb+10), Wrec(kb+10)
!      write(6,*) 'wrav(kb+11), Wrec(kb+11)=',wrav(kb+11), Wrec(kb+11)
!    end if


   ! 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
        tprec(j,k) = thl0(irecy-1,j,k) - Trec(k)
      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
    ! heat flux at recycle station (isothermal wall) q = alpha * dT/dz = (nu/prandtl) * dT/dz
!    q0 = numol*prandtlmoli*(Trec(kb) - Trec(kb-1)) * dzhi(kb)
    q0 = numol*prandtlmoli*2*(Trec(kb) - thls) / dzf(kb)
    ttaur = q0/utaur                ! ttau = q/(rho*cp*utau) =  (alpha *dT/dz) / utau
   ! compute momentum thickness at inlet and recycle plane

   if(lbuoyancy) then
     lmor = (thls* utaur**2 )/ (0.41 * grav * ttaur)      ! L = -T0*utau^3 / kappa*g*<w'T'> =
!     write(6,*) 'Initial dr,myid, utaur, ttaur, Lmor =', dr,myid,utaur,ttaur,lmor
!     lmor = 0.3;
     lmoi = (thls* utaui**2 )/ (0.41 * grav * ttaui)      ! L = -T0*utau^3 / kappa*g*<w'T'> =
!     lmoi = 0.3;
!     write(6,*) 'Initial di_test,myid, utaui, ttaui, Lmoi =', di_test,myid,utaui,ttaui,lmoi
     dr_old = dr

!     call blthicknessmo(dr,utaur,lmor)                    ! Also needed for momentumthickness
     call blthicknesst(dr,Urec,0.99)            ! changed back to this one (instead of the above)

!     call momentumthicknessmo(thetai,utaui,di,lmoi)
!     call momentumthicknessmo(thetar,utaur,dr,lmor)
     call momentumthicknessexp(thetai,Uinl)
     call momentumthicknessexp(thetar,Urec)
   else
!     call blthickness(dr,utaur)                           ! Also needed for momentumthickness
     call blthicknesst(dr,Urec,0.99)
!     call momentumthickness(thetai,utaui,di)
!     call momentumthickness(thetar,utaur,dr)
     call momentumthicknessexp(thetai,Uinl)
     call momentumthicknessexp(thetar,Urec)
   end if
   call enthalpythickness(thetati,Tinl,Uinl)
   call enthalpythickness(thetatr,Trec,Urec)
!   call blthickness(dr,utaur)
   call blthicknesst(dtr,Trec-thls,0.99)
   ! compute utau at inlet from interior field
!    if (thetai == 0.) then
!      write(6,*) '!!! thetai = 0, myid=',myid
!    else if (thetar == 0.) then
!      write(6,*) '!!! thetar = 0, myid=',myid
!      thetar=0.00001
!    else
!      utaui = utaur* (thetar/thetai)**(1./8.)    ! See Lund (1998): 'Similar to Ludwig-Tillmann correlation'
    if (.not.lfixutauin) then
      utaui  = utaur* abs(thetar/thetai)**(1./8.)   ! See Lund (1998): 'Similar to Ludwig-Tillmann correlation'
    end if
    if (thetati == 0.) then
      thetati = 0.0000001
    end if
    ttaui = ttaur*abs(thetatr/thetati)**(1./8.)   ! See Kong (2000):
!    end if
    gamm = utaui / utaur                          ! Gamma in Lund (1998)
    if (ttaur == 0.) then
      ttaur = 0.0000001
    end if
    lamb = ttaui / ttaur                          ! Lambda in Kong (2000)

   ! 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..)
    zotr = zf / dtr                               ! temperature outer scaling zf-coordinate at recycle station
    zoti = zf / dti                               ! temperature outer scaling zf-coordinate at inlet station

!!!!! 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
  !!! Outer coordinates for temperature
    do k=kb,ke
      do kk=kb,ke
        if (zotr(kk) >= zoti(k)) then
          locupot(k)  = kk
          loclowot(k) = kk-1
          exit
        elseif (kk==ke) then
          locupot(k)  = ke+1            ! this means extrapolation!
          loclowot(k) = ke-1
        end if
      end do
    end do
  !!! Finished with outer coordinates temperature

!!! 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)
        Tinli(k) = Trec(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)
!        Tinli(k) = thls + Trec(kb)/zirf(kb)*ziif(k)
!        Tinli(k) = (Trec(kb)-thls)/zirf(kb)*ziif(k)
        Tinli(k) = thls + (Trec(kb)-thls)/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)))
        Tinli(k) = Trec(loclowif(k)) + (Trec(locupif(k)) - Trec(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.
        tpinli(:,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)
        tpinli(:,k) = tprec(:,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)))
        tpinli(:,k) = tprec(:,loclowif(k)) + (tprec(:,locupif(k)) - tprec(:,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
 !! Interpolating outer temperature
   ! mean temperature
    do k=kb,ke
      if (locupot(k) == ke+1) then      ! indicator for extrapolation!
        Tinlo(k) = thl_top
      elseif (loclowot(k) == kb-1) then ! interprets this as extrapolation to bottom (use Tinlo=thls at z+=0)
!        Tinlo(k) = Trec(kb)/zotr(kb) * zoti(k)
!        Tinlo(k) = (Trec(kb)-thls)/zotr(kb) * zoti(k)
        Tinlo(k) = thls + (Trec(kb)-thls)/zotr(kb) * zoti(k)
      else                            ! normal interpolation
        Tinlo(k) = Trec(loclowot(k)) + (Trec(locupot(k)) - Trec(loclowot(k))) / (zotr(locupot(k)) - zotr(loclowot(k))) * (zoti(k) - zotr(loclowot(k)))
      end if
    end do

   ! fluctuating temperature
    do k=kb,ke
      if (locupot(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))
        tpinlo(:,k) = 0.
      elseif (loclowot(k) == kb-1) then ! interprets this as extrapolation to bottom (use t=0 at z+=0)
        tpinlo(:,k) = tprec(:,kb)/zotr(kb) * zoti(k)
      else                            ! normal interpolation
        tpinlo(:,k) = tprec(:,loclowot(k)) + (tprec(:,locupot(k)) - tprec(:,loclowot(k))) / (zotr(locupot(k)) - zotr(loclowot(k))) * (zoti(k) - zotr(loclowot(k)))
      end if
    end do

 !! Finished interpolating out temperature
!!!!! Finished Interpolation! !!!!!


   ! compute rescaled inner variables ! Winli = Winli (interpolation is enough)
    Uinli = gamm* Uinli
    Tinli = lamb* Tinli + (1.-lamb)*thls     ! this is different for isoflux wall!
    upinli = gamm* upinli
    vpinli = gamm* vpinli
    wpinli = gamm* wpinli
    tpinli = lamb* tpinli         ! See Kong (2000)
   ! compute rescaled outer variables ! Winlo = Winlo (interpolation is enough)
    Uinlo = gamm* Uinlo  + (1.- gamm)*Uinf
    Tinlo = lamb* Tinlo  + (1.- lamb)*thl_top
!    Uinlo = gamm* Uinlo  + (1.- gamm)*Urec(ke)
    upinlo = gamm* upinlo
    vpinlo = gamm* vpinlo
    wpinlo = gamm* wpinlo
    tpinlo = lamb* tpinlo         ! See Kong (2000)

!    utop = Uinlo(ke)
!    Uinlo = Uinlo +(Uinf-utop)     ! make sure at the inlet the mean top velocity equals Uinf
!! add defect velocity to make sure the j-averaged velocity at the top equals Uinf
!    utop = Uinlo(ke)
!    do k=kb,ke
!        Uinlo(k) = Uinlo(k)*Uinf/utop
!    end do


   ! 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
    wfunct = 0.5*(1. + tanh( alpha*(zoti-beta)/((1.-2.*beta)*zoti+beta) )/tanh(alpha) ) ! for temperature (full level height)
    do k=kb,ke
      if (wfuncf(k) .gt. 1.) then
        wfuncf(k) = 1.
      end if
      if (wfunct(k) .gt. 1.) then
        wfunct(k) = 1.
      end if
    end do
    do k=kb,ke+1
      if (wfunch(k) .gt. 1.) then
        wfunch(k) = 1.
      end if
    end do



!    write(6,*) 'maxval(wfuncf)=', maxval(wfuncf)
!    write(6,*) 'maxval(wfunch)=', maxval(wfunch)


   ! 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)
!      t0inletbc(j,k) = (Tinli(k)+ tpinli(j,k))*(1.-wfunct(k)) +  (Tinlo(k) + tpinlo(j,k))* wfunct(k)
      u0inletbc(j,k) = (Uinli(k)+ upinli(j,k)*heavif(k))*(1.-wfuncf(k)) +  (Uinlo(k) + upinlo(j,k)*heavif(k))* wfuncf(k)
      v0inletbc(j,k) =            vpinli(j,k)*heavif(k) *(1.-wfuncf(k)) +              vpinlo(j,k)*heavif(k) * wfuncf(k)
      t0inletbc(j,k) = (Tinli(k)+ tpinli(j,k)*heavit(k))*(1.-wfunct(k)) +  (Tinlo(k) + tpinlo(j,k)*heavit(k))* wfunct(k)
    end do
    end do

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

!!    kdamp = kb + floor(0.75*(ke-kb+1))
!    kdamp = kb + 144  ! => zf = 2.24
!    do k=kdamp,ke
!    do j=jb,je
!      if (u0inletbc(j,k) > Uinf) then
!        u0inletbc(j,k) = Uinf
!      end if
!    end do
!    end do



   ! 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
    tinletbc2(ib,jb:je,kb:ke) = t0inletbc(jb:je,kb:ke)  ! this is just a dummy variable to give tninletbc the right dimension in slabsum
    urav = 0.
    trav = 0.
    call slabsum(urav  ,kb,ke,uinletbc2 ,ib,ib,jb,je,kb,ke,ib,ib,jb,je,kb,ke)
    call slabsum(trav  ,kb,ke,tinletbc2 ,ib,ib,jb,je,kb,ke,ib,ib,jb,je,kb,ke)
!    call slabsum(urav  ,kb,ke,u0 ,ib-1,ie+1,jb-1,je+1,kb-1,ke+1,ib,ib,jb,je,kb,ke)
    urav = urav / jtot                    ! average over j-direction
    trav = trav / 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

! rescale the instantaneous profile to keep mass flux constant (tot avoid pressure fluctuations)
    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
      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

!! add defect velocity to make sure the mass flow is the same as the initial mass flow
 !   u0inletbc = u0inletbc + (ubulk-totalu)
 !   urav      = urav      + (ubulk-totalu)

!! add defect velocity to make sure the j-averaged velocity at the top equals Uinf
!    utop = urav(ke)
!    do k=kb,ke
!      do j=jb,je
!        u0inletbc(j,k) = u0inletbc(j,k)*Uinf/utop
!      end do
!      urav(k) = urav(k)*Uinf/utop
!    end do

!    u0inletbc = u0inletbc + (Uinf-utop)
!    urav      = urav      + (Uinf-utop)


!    if (myid==0) then
!    write(6,*) 'u0inletbc(jb+2,ke)', u0inletbc
!    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
    do k=kb,ke
      Tinl(k) =  trav(k)*deltat*avinti + (1.-deltat*avinti)*Tinl(k)
    end do

!    utop = Uinl(ke)
!    Uinl = Uinl +(Uinf-utop)     ! make sure at the inlet the mean top velocity equals Uinf
!    uminletbc = uminletbc + (Uinf-utop)

! write inletplane to array (and to file after 1000 time steps)
    if (lstoreplane ) then
      storeu0inletbc(:,:,nstepread) = u0inletbc(:,:)
      storev0inletbc(:,:,nstepread) = v0inletbc(:,:)
      storew0inletbc(:,:,nstepread) = w0inletbc(:,:)
      storet0inletbc(:,:,nstepread) = t0inletbc(:,:)
      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
      tminletbc = t0inletbc
    end if

   if (lbuoyancy ) then
!     call blthicknessmo(di_test,utaui,lmoi)
     call blthicknesst(di_test,Uinl,0.99)
!     call dispthicknessmo(displ)  ! needed in top BC
     call dispthicknessexp(displ)
   else
!     call blthickness(di_test,utaui)
     call blthicknesst(di_test,Uinl,0.99)
!     call dispthickness(displ)  ! needed in top BC
     call dispthicknessexp(displ)
   end if
   call blthicknesst(dti_test,Tinl-thls,0.99)
   if ((myid==0) .and. (rk3step==3)) then
     write(6,*) 'Inlet Gen: gamma,lambda=',gamm,lamb
     write(6,*) 'Inlet Gen: Uinl(ke),Tinl(ke)=',Uinl(ke),Tinl(ke)
     write(6,*) 'Inlet Gen: utaui,utaur =',utaui,utaur
     write(6,*) 'Inlet Gen: ttaui,ttaur =',ttaui,ttaur
     write(6,*) 'Inlet Gen: Lmoi,Lmor =',lmoi,lmor
     write(6,*) 'Inlet Gen: deltar, deltai_test', dr, di_test
     write(6,*) 'Inlet Gen: deltatr, deltati_test', dtr, dti_test
     write(6,*) 'Inlet Gen: d*i, d*r=',displ(ib),displ(irecy)
     write(6,*) 'Inlet Gen: thetai,thetar',thetai,thetar
     write(6,*) 'Inlet Gen: thetati,thetatr',thetati,thetatr
     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
      t0inletbcold = t0inletbc

  ! 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)
        t0inletbc(:,:) = storet0inletbc(:,:,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)
      t0inletbc(:,:) = (1. - deltat/interval)*t0inletbc(:,:) + (deltat/interval)*storet0inletbc(:,:,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
        tminletbc = t0inletbc
      end if
    end if  ! iinletgen

  end subroutine inletgen