subroutine inletgennotemp
use modglobal, only : ib,ie,jb,je,jgb,jge,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 / (jge-jgb +1) ! average over j-direction
urav = uaver(irecy,:)
wrav = wrav / (jge-jgb +1) ! 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 / (jge-jgb +1) ! 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 / (jge-jgb +1) ! 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