subroutine inletgen
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,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 / (jge-jgb +1) ! average over j-direction
taver = taver / (jge-jgb +1) ! average over j-direction
urav = uaver(irecy,:)
wrav = wrav / (jge-jgb +1) ! average over j-direction
trav = trav / (jge-jgb +1) ! 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 / (jge-jgb +1) ! average over j-direction
trav = trav / (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
! 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 / (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
tminletbc = t0inletbc
end if
end if ! iinletgen
end subroutine inletgen