subroutine readinitfiles
use modfields, only:u0, v0, w0, um, vm, wm, thlm, thl0, thl0h, qtm, qt0, qt0h, &
ql0, ql0h, thv0h, sv0, svm, e12m, e120, &
dudxls, dudyls, dvdxls, dvdyls, dthldxls, dthldyls, &
dqtdxls, dqtdyls, dqtdtls, dpdx, dpdxl, dpdyl, &
wfls, whls, ug, vg, pgx, pgy, uprof, vprof, thlprof, qtprof, e12prof, svprof, &
v0av, u0av, qt0av, ql0av, thl0av, qt0av, sv0av, exnf, exnh, presf, presh, rhof, &
thlpcar, uav, thvh, thvf, IIc, IIcs, IIu, IIus, IIv, IIvs, IIw, IIws
use modglobal, only : ib,ie,ih,jb,je,jh,kb,ke,kh,khc,kmax,dtmax,dt,runtime,timeleft,timee,ntimee,ntrun,btime,dt_lim,nsv,&
zf, zh, dzf, dzh, rv, rd, grav, cp, rlv, pref0, om23_gs, jgb, jge, Uinf, Vinf, dy, &
rslabs, e12min, dzh, dtheta, dqt, dsv, cexpnr, ifinput, lwarmstart, lstratstart, trestart, numol, &
ladaptive, tnextrestart, jmax, linoutflow, lper2inout, iinletgen, lreadminl, &
uflowrate, vflowrate,ltempeq, prandtlmoli, freestreamav, &
tnextfielddump, tfielddump, tsample, tstatsdump, startfile, lprofforc, lchem, k1, JNO2,&
idriver,dtdriver,driverstore,tdriverstart,tdriverdump
use modsubgriddata, only:ekm, ekh
use modsurfdata, only:wtsurf, wqsurf, wsvsurf, &
thls, thvs, ps, qts, svs, sv_top
! use modsurface, only : surface,dthldz
use modboundary, only:boundary, tqaver
use modmpi, only:slabsum, myid, comm3d, mpierr, my_real, avexy_ibm
use modthermodynamics, only:thermodynamics, calc_halflev
use modinletdata, only:Uinl, Urec, Wrec, u0inletbc, v0inletbc, w0inletbc, ubulk, irecy, Utav, Ttav, &
uminletbc, vminletbc, wminletbc, u0inletbcold, v0inletbcold, w0inletbcold, &
storeu0inletbc, storev0inletbc, storew0inletbc, nstepread, nfile, Tinl, &
Trec, tminletbc, t0inletbcold, t0inletbc, storet0inletbc, utaui, ttaui, iangle,&
u0driver,v0driver,w0driver,e120driver,tdriver,thl0driver,qt0driver,storetdriver,&
storeu0driver,storev0driver,storew0driver,storee120driver,storethl0driver,storeqt0driver,&
nstepreaddriver
use modinlet, only:readinletfile
use moddriver, only: readdriverfile,initdriver,drivergen
integer i, j, k, n
real, allocatable :: height(:), th0av(:)
real, dimension(ib - ih:ie + ih, jb - jh:je + jh, kb:ke + kh) :: thv0
real, dimension(kb:ke) :: uaverage ! volume averaged u-velocity
real, dimension(kb:ke) :: vaverage ! volume averaged v-velocity
real, dimension(kb:ke) :: uaverager ! recycle plane
real, dimension(kb:ke) :: uaveragei ! inlet plane
real, dimension(kb:ke) :: taverager ! recycle plane
real, dimension(kb:ke) :: taveragei ! inlet plane
real, dimension(kb:ke + 1) :: waverage
real, dimension(kb:ke + 1) :: uprofrot
real, dimension(kb:ke + 1) :: vprofrot
real tv, ran, ran1, vbulk
character(80) chmess
allocate (height(kb:ke + kh))
allocate (th0av(kb:ke + kh))
if (lstratstart) then ! Switch
! Read restart files as in lwarmstart
call readrestartfiles
um = u0
vm = v0
wm = w0
thlm = thl0 !do this before or just not?
qtm = qt0
svm = sv0
e12m = e120
! Overwrite thlm, thl0, qtm, qt0 from prof.inp.xxx
if (myid == 0) then
open (ifinput, file='prof.inp.'//cexpnr)
read (ifinput, '(a80)') chmess
write (*, '(a80)') chmess
read (ifinput, '(a80)') chmess
do k = kb, ke
read (ifinput, *) &
height(k), &
thlprof(k), &
qtprof(k), &
uprof(k), &
vprof(k), &
e12prof(k)
end do
close (ifinput)
write (*, *) 'height thl qt u v e12'
do k = ke, kb, -1
write (*, '(f7.1,2f8.1,3f7.1)') &
height(k), &
thlprof(k), &
qtprof(k), &
uprof(k), &
vprof(k), &
e12prof(k)
end do
end if !myid=0
! MPI broadcast thl and qt
call MPI_BCAST(thlprof, kmax, MY_REAL, 0, comm3d, mpierr)
call MPI_BCAST(qtprof, kmax, MY_REAL, 0, comm3d, mpierr)
do k = kb, ke
do j = jb - 1, je + 1
do i = ib - 1, ie + 1
thl0(i, j, k) = thlprof(k)
thlm(i, j, k) = thlprof(k)
qt0(i, j, k) = qtprof(k)
qtm(i, j, k) = qtprof(k)
end do
end do
end do
!ILS13 reintroduced thv !tg3315 this part may wrong, could need to use
call calc_halflev
! exnf = (presf/pref0)**(rd/cp) !exner functions not in restart files
! anymore.. or at least not read
! exnh = (presh/pref0)**(rd/cp)
! write(*,*) "exnf",enf
! write(*,*) "exnh",exnh
do k = kb, ke + kh
do j = jb, je
do i = ib, ie
!write(*,*) "thl0h",thl0h(i,j,k)
thv0h(i, j, k) = (thl0h(i, j, k) + rlv*ql0h(i, j, k)/(cp)) &
*(1 + (rv/rd - 1)*qt0h(i, j, k) - rv/rd*ql0h(i, j, k))
end do
end do
end do
do j = j, je
do i = ib, ie
do k = kb, ke + kh
thv0(i, j, k) = (thl0(i, j, k) + rlv*ql0(i, j, k)/(cp)) &
*(1 + (rv/rd - 1)*qt0(i, j, k) - rv/rd*ql0(i, j, k))
end do
end do
end do
thvh = 0.
! call slabsum(thvh,kb,ke,thv0h,ib-ih,ie+ih,jb-jh,je+jh,kb-kh,ke+kh,ib,ie,jb,je,kb,ke) ! redefine halflevel thv using calculated thv
call avexy_ibm(thvh(kb:ke+kh),thv0h(ib:ie,jb:je,kb:ke+kh),ib,ie,jb,je,kb,ke,ih,jh,kh,IIw(ib:ie,jb:je,kb:ke+kh),IIws(kb:ke+kh),.false.)
! thvh = thvh/rslabs
thvf = 0.0
call avexy_ibm(thvf(kb:ke+kh),thv0(ib:ie,jb:je,kb:ke+kh),ib,ie,jb,je,kb,ke,ih,jh,kh,IIc(ib:ie,jb:je,kb:ke+kh),IIcs(kb:ke+kh),.false.)
! call slabsum(thvf,kb,ke,thv0,ib-ih,ie+ih,jb-jh,je+jh,kb-kh,ke+kh,ib,ie,jb,je,kb,ke)
! thvf = thvf/rslabs
else !if not lstratstart
if (.not. lwarmstart) then
!********************************************************************
! 1.0 prepare initial fields from files 'prof.inp' and 'scalar.inp'
! ----------------------------------------------------------------
!--------------------------------------------------------------------
! 1.1 read fields
!-----------------------------------------------------------------
dt = dtmax/100.
timee = 0.
if (myid == 0) then
open (ifinput, file='prof.inp.'//cexpnr)
read (ifinput, '(a80)') chmess
write (*, '(a80)') chmess
read (ifinput, '(a80)') chmess
do k = kb, ke
read (ifinput, *) &
height(k), &
thlprof(k), &
qtprof(k), &
uprof(k), &
vprof(k), &
e12prof(k)
end do
! Apply rotation in horizontal
!write (6, *) 'iangle = ', iangle
!uprofrot = uprof*cos(iangle) - vprof*sin(iangle)
!vprofrot = vprof*cos(iangle) + uprof*sin(iangle)
!uprof = uprofrot
!vprof = vprofrot
close (ifinput)
write (*, *) 'height thl qt u v e12'
do k = ke, kb, -1
write (*, '(f7.1,2f8.1,3f7.1)') &
height(k), &
thlprof(k), &
qtprof(k), &
uprof(k), &
vprof(k), &
e12prof(k)
end do
if (minval(e12prof(kb:ke)) < e12min) then
write (*, *) 'e12 value is zero (or less) in prof.inp'
do k = kb, ke
e12prof(k) = max(e12prof(k), e12min)
end do
end if
end if ! end if myid==0
! MPI broadcast numbers reading
call MPI_BCAST(thlprof, kmax, MY_REAL, 0, comm3d, mpierr)
call MPI_BCAST(qtprof, kmax, MY_REAL, 0, comm3d, mpierr)
call MPI_BCAST(uprof, kmax, MY_REAL, 0, comm3d, mpierr)
call MPI_BCAST(vprof, kmax, MY_REAL, 0, comm3d, mpierr)
call MPI_BCAST(e12prof, kmax, MY_REAL, 0, comm3d, mpierr)
do k = kb, ke
do j = jb - 1, je + 1
do i = ib - 1, ie + 1
thl0(i, j, k) = thlprof(k)
thlm(i, j, k) = thlprof(k)
qt0(i, j, k) = qtprof(k)
qtm(i, j, k) = qtprof(k)
u0(i, j, k) = uprof(k)
um(i, j, k) = uprof(k)
v0(i, j, k) = vprof(k)
vm(i, j, k) = vprof(k)
w0(i, j, k) = 0.0
wm(i, j, k) = 0.0
e120(i, j, k) = e12prof(k)
e12m(i, j, k) = e12prof(k)
! ekm (i,j,k) = 0.0
! ekh (i,j,k) = 0.0
ekm(i, j, k) = numol
ekh(i, j, k) = numol
end do
end do
end do
ekh(:, :, ke + 1) = ekh(:, :, ke) ! also for start up
! ILS13 30.11.17, added, not sure if necessary
! ILS13 30.11.1, commented
do j = jb - jh, je + jh
do i = ib - ih, ie + ih
thl0(i, j, ke + 1) = thl0(i, j, ke)
thl0(i, j, kb - 1) = thl0(i, j, kb)
end do
end do
!! add random fluctuations
krand = min(krand, ke)
do k = kb, krand
call randomnize(um, k, randu, irandom, ih, jh)
end do
do k = kb, krand
call randomnize(vm, k, randu, irandom, ih, jh)
end do
do k = kb, krand
call randomnize(wm, k, randu, irandom, ih, jh)
end do
! do k=kb+1,ke-1
! do j=jb,je
! do i=ib+1,ie-1
! call random_number(ran)
! ran1 = -1. +2.*ran
! wm(i,j,k)=wm(i,j,k)+ 0.1*Uinf*ran1
! end do
! end do
! end do
!
! do k=kb+1,ke-1
! do j=jb,je
! do i=ib+2,ie-1
! call random_number(ran)
! ran1 = -1. +2.*ran
! um(i,j,k)=um(i,j,k)+ 0.1*Uinf*ran1
! end do
! end do
! end do
!
! do k=kb+1,ke-1
! do j=jb,je
! do i=ib+1,ie-1
! call random_number(ran)
! ran1 = -1. +2.*ran
! vm(i,j,k)=vm(i,j,k)+ 0.1*Uinf*ran1
! end do
! end do
! end do
u0 = um
v0 = vm
w0 = wm
uaverage = 0.
! call slabsum(uaverage, kb, ke, um, ib - 1, ie + 1, jb - 1, je + 1, kb - 1, ke + 1, ib, ie, jb, je, kb, ke)
do k = kb, ke
uaverage(k) = uprof(k)*dzf(k)
end do
ubulk = sum(uaverage(kb:ke))/(zh(ke + 1) - zh(kb)) ! averaged u-velocity inflow profile
vaverage = 0.
! call slabsum(vaverage, kb, ke, vm, ib - 1, ie + 1, jb - 1, je + 1, kb - 1, ke + 1, ib, ie, jb, je, kb, ke)
do k = kb, ke
vaverage(k) = vprof(k)*dzf(k)
end do
vbulk = sum(vaverage(kb:ke))/(zh(ke + 1) - zh(kb)) ! averaged u-velocity inflow profile
! Set average inlet profile to initial inlet profile in case of inletgenerator mode
if (iinletgen == 1) then
uaverage = 0.
call slabsum(uaverage, kb, ke, um, ib - 1, ie + 1, jb - 1, je + 1, kb - 1, ke + 1, ib, ie, jb, je, kb, ke)
do k = kb, ke
uaverage(k) = uprof(k)*dzf(k)
end do
ubulk = sum(uaverage(kb:ke))/(zh(ke + 1) - zh(kb)) ! volume-averaged u-velocity
write (6, *) 'Modstartup: ubulk=', ubulk
Utav(ib:ie, kb:ke) = um(ib:ie, jb, kb:ke)
Uinl = um(ib, jb, kb:ke) ! set the initial time-averaged inlet profile equal to um
Urec = um(ib, jb, kb:ke) ! set the initial time-averaged inlet profile equal to um
Wrec(kb:ke + 1) = wm(ib, jb, kb:ke + 1) ! set the initial time-averaged inlet profile equal to mean w-profile
u0inletbcold(jb:je, kb:ke) = um(ib, jb:je, kb:ke)
v0inletbcold(jb:je, kb:ke) = vm(ib - 1, jb:je, kb:ke)
w0inletbcold(jb:je, kb:ke + 1) = wm(ib - 1, jb:je, kb:ke + 1)
uminletbc(jb:je, kb:ke) = um(ib, jb:je, kb:ke)
vminletbc(jb:je, kb:ke) = vm(ib - 1, jb:je, kb:ke)
wminletbc(jb:je, kb:ke) = wm(ib - 1, jb:je, kb:ke)
u0inletbc(jb:je, kb:ke) = um(ib, jb:je, kb:ke)
v0inletbc(jb:je, kb:ke) = vm(ib - 1, jb:je, kb:ke)
w0inletbc(jb:je, kb:ke + 1) = wm(ib - 1, jb:je, kb:ke + 1)
utaui = sqrt(abs(2*numol*Uinl(kb)/dzf(kb))) ! average streamwise friction at inlet (need for first time step)
if (ltempeq) then
Ttav(ib:ie, kb:ke) = thlm(ib:ie, jb, kb:ke) ! set the initial time-averaged inlet profile equal to thlm
Tinl = thlm(ib, jb, kb:ke) ! set the initial time-averaged inlet profile equal to thlm
Trec = thlm(ib, jb, kb:ke) ! set the initial time-averaged inlet profile equal to thlm
t0inletbcold(jb:je, kb:ke) = thlm(ib - 1, jb:je, kb:ke)
t0inletbc(jb:je, kb:ke) = thl0(ib - 1, jb:je, kb:ke)
tminletbc(jb:je, kb:ke) = thlm(ib - 1, jb:je, kb:ke)
ttaui = numol*prandtlmoli*2.*(Tinl(kb) - thls)/(dzf(kb)*utaui) ! average friction temp. at inlet (need for first time step)
end if
! add random perturbations
if (myid == 0) then
call random_number(ran)
ran1 = -1.+2.*ran
write (6, *) 'random=', ran, ran1
call random_number(ran)
ran1 = -1.+2.*ran
write (6, *) 'random=', ran, ran1
call random_number(ran)
ran1 = -1.+2.*ran
write (6, *) 'random=', ran, ran1
call random_number(ran)
ran1 = -1.+2.*ran
write (6, *) 'random=', ran, ran1
call random_number(ran)
ran1 = -1.+2.*ran
write (6, *) 'random=', ran, ran1
call random_number(ran)
ran1 = -1.+2.*ran
write (6, *) 'random=', ran, ran1
end if
do k = kb + 1, kb + 48
do j = jb, je
do i = ib + 1, ie - 1
call random_number(ran)
ran1 = -1.+2.*ran
wm(i, j, k) = wm(i, j, k) + 0.1*Uinf*ran1
end do
end do
end do
! do k=kb+1,ke-1
do k = kb + 1, kb + 48
do j = jb, je
do i = ib + 2, ie - 1
call random_number(ran)
ran1 = -1.+2.*ran
um(i, j, k) = um(i, j, k) + 0.1*Uinf*ran1
end do
end do
end do
! do k=kb+1,ke-1
do k = kb + 1, kb + 48
do j = jb, je
do i = ib + 1, ie - 1
call random_number(ran)
ran1 = -1.+2.*ran
vm(i, j, k) = vm(i, j, k) + 0.1*Uinf*ran1
end do
end do
end do
u0 = um
v0 = vm
w0 = wm
else if (iinletgen == 2) then
nfile = nfile + 1
call readinletfile
u0inletbc(:, :) = storeu0inletbc(:, :, nstepread)
v0inletbc(:, :) = storev0inletbc(:, :, nstepread)
w0inletbc(:, :) = storew0inletbc(:, :, nstepread)
uminletbc(:, :) = storeu0inletbc(:, :, nstepread)
vminletbc(:, :) = storev0inletbc(:, :, nstepread)
wminletbc(:, :) = storew0inletbc(:, :, nstepread)
! determine bulk velocity
call slabsum(uaverage, kb, ke, u0, ib - 1, ie + 1, jb - 1, je + 1, kb - 1, ke + 1, ib, ie, jb, je, kb, ke)
uaverage = uaverage/((ie - ib + 1)*(jge - jgb + 1)) ! this gives the i-j-averaged velocity (only correct for equidistant grid?)
do k = kb, ke
uaverage(k) = uaverage(k)*dzf(k)
end do
ubulk = sum(uaverage(kb:ke))/(zh(ke + 1) - zh(kb)) ! volume-averaged u-velocity
write (6, *) 'Modstartup: ubulk=', ubulk
elseif (idriver==2) then ! idriver
call readdriverfile
! if(myid==0) then
! write(*,*) 'Driver inlet velocity'
! do n=1,driverstore
! write (*,'(f9.2,e20.12)') storetdriver(n), storeu0driver(1,32,n)
! end do
! endif
! call slabsum(uaverage,kb,ke,u0,ib-1,ie+1,jb-1,je+1,kb-1,ke+1,ib,ie,jb,je,kb,ke)
! uaverage = uaverage / ((ie-ib+1)*(jge-jgb+1)) ! this gives the i-j-averaged velocity (only correct for equidistant grid?)
call avexy_ibm(uaverage(kb:ke),u0(ib:ie,jb:je,kb:ke),ib,ie,jb,je,kb,ke,ih,jh,kh,IIu(ib:ie,jb:je,kb:ke),IIus(kb:ke),.false.)
do k=kb,ke
uaverage(k) = uaverage(k)*dzf(k)
end do
ubulk = sum(uaverage(kb:ke))/(zh(ke+1)-zh(kb)) !volume-averaged u-velocity
if (myid==0) then
write(6,*) 'Modstartup: ubulk=',ubulk
end if
elseif (idriver==1) then
call drivergen
end if
!---------------------------------------------------------------
! 1.2 randomnize fields
!---------------------------------------------------------------
! if (iinletgen /= 2 .and. iinletgen /= 1) then
! write(6,*) 'randomnizing temperature!'
! krand = min(krand,ke)
! do k = kb,ke !edited tg3315 krand --> ke
! call randomnize(thlm,k,randthl,irandom,ih,jh)
! call randomnize(thl0,k,randthl,irandom,ih,jh)
! end do
! end if
svprof = 0.
if (myid == 0) then
if (nsv > 0) then
open (ifinput, file='scalar.inp.'//cexpnr)
read (ifinput, '(a80)') chmess
read (ifinput, '(a80)') chmess
do k = kb, ke
read (ifinput, *) &
height(k), &
(svprof(k, n), n=1, nsv)
end do
open (ifinput, file='scalar.inp.'//cexpnr)
write (6, *) 'height sv(1) --------- sv(nsv) '
do k = ke, kb, -1
write (6, *) &
height(k), &
(svprof(k, n), n=1, nsv)
end do
end if
end if ! end if myid==0
call MPI_BCAST(svprof, (ke + kh - (kb - kh))*nsv, MY_REAL, 0, comm3d, mpierr)
do k = kb, ke
do j = jb - 1, je + 1
do i = ib - 1, ie + 1
do n = 1, nsv
sv0(i, j, k, n) = svprof(k, n)
svm(i, j, k, n) = svprof(k, n)
end do
end do
end do
end do
if (nsv>0) then !tg3315 set these variables here for now and repeat for warmstart
allocate(sv_top(1:nsv))
sv_top(:) = svprof(ke,1:nsv)
call MPI_BCAST(sv_top, nsv, MY_REAL, 0, comm3d, mpierr)
write(*,*) 'svprof', svprof
write(*,*) 'sv_top', sv_top
end if
!do n = 1,nsv
! do j = jb - jhc, je + jhc
! do i = ib - ihc, ie + ihc
! svm(i, j, ke + 1, n) = svm(i, j, ke)
! sv0(i, j, kb - 1, n) = sv0(i, j, kb)
! end do
! end do
!end do
!-----------------------------------------------------------------
! 2.2 Initialize surface layer
!-----------------------------------------------------------------
!ILS13 reintroduced thv !tg3315 this part may wrong, could need to use
call calc_halflev
! exnf = (presf/pref0)**(rd/cp) !exner functions not in restart files
! anymore.. or at least not read
! exnh = (presh/pref0)**(rd/cp)
call boundary ! tg3315 17.10.17 having this in startup was causing issues for running with lmoist ! turned of when pot. temp = temp.
call thermodynamics ! turned off when pot. temp = temp.
call boundary
call thermodynamics ! turned off when pot. temp = temp.
else !if lwarmstart
write (*, *) "doing warmstart"
call readrestartfiles
um = u0
vm = v0
wm = w0
thlm = thl0
qtm = qt0
svm = sv0
e12m = e120
ekm(:, :, :) = numol
ekh(:, :, :) = numol*prandtlmoli !tg3315 added because wttop using ekh in modboundary which is called in startup
ekh(:, :, ke + 1) = ekh(:, :, ke) ! also for start up
if (idriver==1) then
!driverstore = (timeleft - tdriverstart)/dtdriver + 1
!if(myid==0) then
! write(*,*) 'driverstore: ', driverstore
!end if
call drivergen
tdriverdump = tdriverstart
endif
!ILS13 reintroduced thv
call calc_halflev
! exnf = (presf/pref0)**(rd/cp) !exner functions not in restart files
! anymore.. or at least not read
! exnh = (presh/pref0)**(rd/cp)
! write(*,*) "exnf",enf
! write(*,*) "exnh",exnh
do j = jb, je
do i = ib, ie
do k = kb, ke + kh
!write(*,*) "thl0h",thl0h(i,j,k)
thv0h(i, j, k) = (thl0h(i, j, k) + rlv*ql0h(i, j, k)/(cp)) &
*(1 + (rv/rd - 1)*qt0h(i, j, k) - rv/rd*ql0h(i, j, k))
end do
end do
end do
do j = j, je
do i = ib, ie
do k = kb, ke + kh
thv0(i, j, k) = (thl0(i, j, k) + rlv*ql0(i, j, k)/(cp)) &
*(1 + (rv/rd - 1)*qt0(i, j, k) - rv/rd*ql0(i, j, k))
end do
end do
end do
thvh = 0.
! call slabsum(thvh,kb,ke,thv0h,ib-ih,ie+ih,jb-jh,je+jh,kb-kh,ke+kh,ib,ie,jb,je,kb,ke) ! redefine halflevel thv using calculated thv
call avexy_ibm(thvh(kb:ke+kh),thv0h(ib:ie,jb:je,kb:ke+kh),ib,ie,jb,je,kb,ke,ih,jh,kh,IIw(ib:ie,jb:je,kb:ke+kh),IIws(kb:ke+kh),.false.)
! thvh = thvh/rslabs
thvf = 0.0
call avexy_ibm(thvf(kb:ke+kh),thv0(ib:ie,jb:je,kb:ke+kh),ib,ie,jb,je,kb,ke,ih,jh,kh,IIc(ib:ie,jb:je,kb:ke+kh),IIcs(kb:ke+kh),.false.)
! call slabsum(thvf,kb,ke,thv0,ib-ih,ie+ih,jb-jh,je+jh,kb-kh,ke+kh,ib,ie,jb,je,kb,ke)
! thvf = thvf/rslabs
! Set average inlet profile to initial inlet profile in case of inletgenerator mode
uaverage = 0.
uaveragei = 0.
uaverager = 0.
waverage = 0.
taveragei = 0.
taverager = 0.
if (iinletgen == 1) then
call slabsum(uaveragei, kb, ke, u0, ib - 1, ie + 1, jb - 1, je + 1, kb - 1, ke + 1, ib, ib, jb, je, kb, ke)
call slabsum(uaverager, kb, ke, u0, ib - 1, ie + 1, jb - 1, je + 1, kb - 1, ke + 1, irecy, irecy, jb, je, kb, ke)
call slabsum(waverage, kb, ke + 1, w0, ib - 1, ie + 1, jb - 1, je + 1, kb - 1, ke + 1, ib, ie, jb, je, kb, ke + 1)
call slabsum(uaverage, kb, ke, u0, ib - 1, ie + 1, jb - 1, je + 1, kb - 1, ke + 1, ib, ie, jb, je, kb, ke)
uaverage = uaverage/((ie - ib + 1)*(jge - jgb + 1)) ! this gives the i-j-averaged velocity (only correct for equidistant grid?)
uaveragei = uaveragei/(jge - jgb + 1) ! this gives the j-averaged u-velocity at the inlet
uaverager = uaverager/(jge - jgb + 1) ! this gives the j-averaged u-velocity at the recycle plane
waverage = waverage/((ie - ib + 1)*(jge - jgb + 1)) ! this gives the i-j-averaged w-velocity (only correct for equidistant grid?)
if (ltempeq) then
call slabsum(taveragei, kb, ke, thl0, ib - 1, ie + 1, jb - 1, je + 1, kb - 1, ke + 1, ib, ie, jb, je, kb, ke)
call slabsum(taverager, kb, ke, thl0, ib - 1, ie + 1, jb - 1, je + 1, kb - 1, ke + 1, irecy - 1, irecy - 1, jb, je, kb, ke)
taveragei = taveragei/((ie - ib + 1)*(jge - jgb + 1)) ! this gives the j-averaged temperature at the inlet
taverager = taverager/(jge - jgb + 1) ! this gives the j-averaged temperature at the recycle plane
end if
if (.not. lreadminl) then
if (myid == 0) then
write (6, *) 'uaverage(kb)=', uaverage(kb)
write (6, *) 'uaverage(ke)=', uaverage(ke)
write (6, *) 'waverage(ke)=', waverage(ke)
write (6, *) 'waverage(ke-20)=', waverage(ke - 20)
write (6, *) 'taveragei(kb)=', taveragei(kb)
write (6, *) 'taveragei(ke)=', taveragei(ke)
end if
Utav = 0.
do i = ib, ie
Utav(i, :) = uaverage
end do
Uinl = uaverage ! set the initial time-averaged inlet profile equal to mean u-profile read from means
write (6, *) 'Uinl(kb+10)=', Uinl(kb + 10)
utaui = sqrt(abs(2*numol*Uinl(kb)/dzf(kb))) ! average streamwise friction at inlet (need for first time step)
Urec = uaverage ! set the initial time-averaged inlet profile equal to mean u-profile
Wrec(kb:ke + 1) = waverage(kb:ke + 1) ! set the initial time-averaged inlet profile equal to mean w-profile
Wrec(kb) = 0. ! set the initial time-averaged inlet profile equal to zero
if (ltempeq) then
Ttav = 0.
do i = ib, ie
Ttav(i, :) = taveragei(:)
end do
Tinl = taveragei
Trec = taveragei
ttaui = numol*prandtlmoli*2.*(Tinl(kb) - thls)/(dzf(kb)*utaui) ! friction temp. at inlet (need at first time step)
end if
else ! -> lreadminl -> Uinl, Urec, Wrec already read
call slabsum(uaverage, kb, ke, u0, ib - 1, ie + 1, jb - 1, je + 1, kb - 1, ke + 1, ib, ie, jb, je, kb, ke)
uaverage = uaverage/((ie - ib + 1)*(jge - jgb + 1)) ! this gives the i-j-averaged velocity (only correct for equidistant grid?)
end if
! determine bulk velocity
do k = kb, ke
uaverage(k) = uaverage(k)*dzf(k)
end do
ubulk = sum(uaverage(kb:ke))/(zh(ke + 1) - zh(kb)) ! volume-averaged u-velocity
write (6, *) 'Modstartup: ubulk=', ubulk
do k = kb, ke
do j = jb, je
uminletbc(j, k) = um(ib, j, k)
vminletbc(j, k) = vm(ib - 1, j, k)
u0inletbcold(j, k) = um(ib, j, k)
v0inletbcold(j, k) = vm(ib - 1, j, k)
u0inletbc(j, k) = um(ib, j, k)
v0inletbc(j, k) = vm(ib - 1, j, k)
end do
end do
do k = kb, ke + 1
do j = jb, je
wminletbc(j, k) = wm(ib - 1, j, k)
w0inletbcold(j, k) = wm(ib - 1, j, k)
w0inletbc(j, k) = wm(ib - 1, j, k)
end do
end do
if (ltempeq) then
do k = kb, ke
do j = jb, je
tminletbc(j, k) = thlm(ib - 1, j, k)
t0inletbcold(j, k) = thlm(ib - 1, j, k)
t0inletbc(j, k) = thlm(ib - 1, j, k)
end do
end do
end if
write (6, *) 'uminletbc(jb,kb),um(ib,jb,kb)=', uminletbc(jb, kb), um(ib, jb, kb)
write (6, *) 'uminletbc(jb+1,kb+10),um(ib,jb+1,kb+10)=', uminletbc(jb + 1, kb + 10), um(ib, jb + 1, kb + 10)
write (6, *) 'uminletbc(je,kb+10),um(ib,je,kb+10)=', uminletbc(je, kb + 10), um(ib, je, kb + 10)
else if (iinletgen == 2) then
nfile = nfile + 1
write (6, *) 'Loading inletfile'
call readinletfile
u0inletbc(:, :) = storeu0inletbc(:, :, nstepread)
v0inletbc(:, :) = storev0inletbc(:, :, nstepread)
w0inletbc(:, :) = storew0inletbc(:, :, nstepread)
uminletbc(:, :) = storeu0inletbc(:, :, nstepread)
vminletbc(:, :) = storev0inletbc(:, :, nstepread)
wminletbc(:, :) = storew0inletbc(:, :, nstepread)
if (ltempeq) then
t0inletbc(:, :) = storet0inletbc(:, :, nstepread)
tminletbc(:, :) = storet0inletbc(:, :, nstepread)
end if
! determine bulk velocity
call slabsum(uaverage, kb, ke, u0, ib - 1, ie + 1, jb - 1, je + 1, kb - 1, ke + 1, ib, ie, jb, je, kb, ke)
uaverage = uaverage/((ie - ib + 1)*(jge - jgb + 1)) ! this gives the i-j-averaged velocity (only correct for equidistant grid?)
do k = kb, ke
uaverage(k) = uaverage(k)*dzf(k)
end do
ubulk = sum(uaverage(kb:ke))/(zh(ke + 1) - zh(kb)) ! volume-averaged u-velocity
write (6, *) 'Modstartup: ubulk=', ubulk
elseif (idriver==2) then ! idriver
call readdriverfile
call drivergen
!call slabsum(uaverage,kb,ke,u0,ib-1,ie+1,jb-1,je+1,kb-1,ke+1,ib,ie,jb,je,kb,ke)
!uaverage = uaverage / ((ie-ib+1)*(jge-jgb+1)) ! this gives the i-j-averaged velocity (only correct for equidistant grid?)
call avexy_ibm(uaverage(kb:ke),u0(ib:ie,jb:je,kb:ke),ib,ie,jb,je,kb,ke,ih,jh,kh,IIu(ib:ie,jb:je,kb:ke),IIus(kb:ke),.false.)
do k=kb,ke
uaverage(k) = uaverage(k)*dzf(k)
end do
ubulk = sum(uaverage(kb:ke))/(zh(ke+1)-zh(kb)) !volume-averaged u-velocity
if (myid==0) then
write(6,*) 'Modstartup: ubulk=',ubulk
end if
end if ! iinletgen/idriver
if (lper2inout) then ! if the restart starts from a periodic simulation to in/outflow, lper2inout should be set to .true.
if (myid == 0) then
write (6, *) 'per2inout=.true. -> reading inlet profile from prof.inp.XXX and scalar.inp.XXX'
open (ifinput, file='prof.inp.'//cexpnr) ! read the inlet profile from prof.inp
read (ifinput, '(a80)') chmess
write (*, '(a80)') chmess
read (ifinput, '(a80)') chmess
do k = kb, ke
read (ifinput, *) &
height(k), &
thlprof(k), &
qtprof(k), &
uprof(k), &
vprof(k), &
e12prof(k)
end do
svprof = 0.
if (nsv > 0) then
open (ifinput, file='scalar.inp.'//cexpnr)
read (ifinput, '(a80)') chmess
read (ifinput, '(a80)') chmess
do k = kb, ke
read (ifinput, *) &
height(k), &
(svprof(k, n), n=1, nsv)
end do
open (ifinput, file='scalar.inp.'//cexpnr)
write (6, *) 'height sv(1) --------- sv(nsv) '
do k = ke, kb, -1
write (6, *) &
height(k), &
(svprof(k, n), n=1, nsv)
end do
end if
end if ! end if myid==0
! MPI broadcast numbers reading
call MPI_BCAST(thlprof, kmax, MY_REAL, 0, comm3d, mpierr)
call MPI_BCAST(uprof, kmax, MY_REAL, 0, comm3d, mpierr)
call MPI_BCAST(vprof, kmax, MY_REAL, 0, comm3d, mpierr)
call MPI_BCAST(e12prof, kmax, MY_REAL, 0, comm3d, mpierr)
call MPI_BCAST(qtprof, kmax, MY_REAL, 0, comm3d, mpierr)
call MPI_BCAST(svprof, (ke + kh - (kb - kh))*nsv, MY_REAL, 0, comm3d, mpierr)
else if (linoutflow) then ! restart of inoutflow simulation: reproduce inlet boundary condition from restartfile
do j = jb - 1, je + 1
do k = kb, ke + 1
uprof(k) = u0(ib, j, k)
vprof(k) = (v0(ib - 1, j, k) + v0(ib, j, k))/2
thlprof(k) = (thl0(ib - 1, j, k) + thl0(ib, j, k))/2
qtprof(k) = (qt0(ib - 1, j, k) + qt0(ib, j, k))/2
e12prof(k) = (e120(ib - 1, j, k) + e120(ib, j, k))/2
do n = 1, nsv
svprof(k, n) = (sv0(ib - 1, j, k, n) + sv0(ib, j, k, n))/2
enddo
enddo
enddo
! outlet bulk velocity
call slabsum(uaverage, kb, ke, u0, ib - 1, ie + 1, jb - 1, je + 1, kb - 1, ke + 1, ib, ie, jb, je, kb, ke)
uaverage = uaverage/((ie - ib + 1)*(jge - jgb + 1)) ! this gives the i-j-averaged velocity (only correct for equidistant grid?)
! determine bulk velocity
do k = kb, ke
uaverage(k) = uaverage(k)*dzf(k)
end do
ubulk = sum(uaverage(kb:ke))/(zh(ke + 1) - zh(kb)) ! volume-averaged u-velocity
write (6, *) 'Modstartup: ubulk=', ubulk
else ! else per2per... read svprof regardless...
! tg3315 read svprof (but do not use regardless of above...)
svprof = 0.
if (myid == 0) then
if (nsv > 0) then
open (ifinput, file='scalar.inp.'//cexpnr)
read (ifinput, '(a80)') chmess
read (ifinput, '(a80)') chmess
do k = kb, ke
read (ifinput, *) &
height(k), &
(svprof(k, n), n=1, nsv)
end do
open (ifinput, file='scalar.inp.'//cexpnr)
write (6, *) 'height sv(1) --------- sv(nsv) '
do k = ke, kb, -1
write (6, *) &
height(k), &
(svprof(k, n), n=1, nsv)
end do
end if
end if ! end if myid==0
call MPI_BCAST(svprof, (ke + kh - (kb - kh))*nsv, MY_REAL, 0, comm3d, mpierr)
if (nsv>0) then !tg3315 set these variables here for now and repeat for warmstart
allocate(sv_top(1:nsv))
sv_top(:) = svprof(ke,1:nsv)
call MPI_BCAST(sv_top, nsv, MY_REAL, 0, comm3d, mpierr)
end if
end if ! end if lper2inout
u0av = 0.0
v0av = 0.0
thl0av = 0.0
qt0av = 0.0
th0av = 0.0
sv0av = 0.
! call slabsum(u0av ,kb,ke+kh,u0(:,:,kb:ke+kh) ,ib-ih,ie+ih,jb-jh,je+jh,kb,ke+kh,ib,ie,jb,je,kb,ke+kh)
call avexy_ibm(u0av(kb:ke+kh),u0(ib:ie,jb:je,kb:ke+kh),ib,ie,jb,je,kb,ke,ih,jh,kh,IIu(ib:ie,jb:je,kb:ke+kh),IIus(kb:ke+kh),.false.)
! call slabsum(v0av ,kb,ke+kh,v0(:,:,kb:ke+kh) ,ib-ih,ie+ih,jb-jh,je+jh,kb,ke+kh,ib,ie,jb,je,kb,ke+kh)
call avexy_ibm(v0av(kb:ke+kh),v0(ib:ie,jb:je,kb:ke+kh),ib,ie,jb,je,kb,ke,ih,jh,kh,IIv(ib:ie,jb:je,kb:ke+kh),IIvs(kb:ke+kh),.false.)
! call slabsum(thl0av,kb,ke+kh,thl0(:,:,kb:ke+kh),ib-ih,ie+ih,jb-jh,je+jh,kb,ke+kh,ib,ie,jb,je,kb,ke+kh)
call avexy_ibm(thl0av(kb:ke+kh),thl0(ib:ie,jb:je,kb:ke+kh),ib,ie,jb,je,kb,ke,ih,jh,kh,IIc(ib:ie,jb:je,kb:ke+kh),IIcs(kb:ke+kh),.false.)
! call slabsum(qt0av,kb,ke+kh,qt0(:,:,kb:ke+kh),ib-ih,ie+ih,jb-jh,je+jh,kb,ke+kh,ib,ie,jb,je,kb,ke+kh)
call avexy_ibm(qt0av(kb:ke+kh),qt0(ib:ie,jb:je,kb:ke+kh),ib,ie,jb,je,kb,ke,ih,jh,kh,IIc(ib:ie,jb:je,kb:ke+kh),IIcs(kb:ke+kh),.false.)
do n = 1, nsv
! call slabsum(sv0av(kb,n),kb,ke+kh,sv0(ib-ih,jb-jh,kb,n),ib-ih,ie+ih,jb-jh,je+jh,kb,ke+kh,ib,ie,jb,je,kb,ke+kh)
call avexy_ibm(sv0av(kb:ke+khc,n),sv0(ib:ie,jb:je,kb:ke+khc,n),ib,ie,jb,je,kb,ke,ih,jh,kh,IIc(ib:ie,jb:je,kb:ke+khc),IIcs(kb:ke+khc),.false.)
end do
! CvH - only do this for fixed timestepping. In adaptive dt comes from restartfile
if (.not. ladaptive) dt = dtmax
! call boundary
end if ! lwarmstart
end if ! not lstratstart
!-----------------------------------------------------------------
! 2.1 read and initialise fields
!-----------------------------------------------------------------
if (myid == 0) then
open (ifinput, file='lscale.inp.'//cexpnr)
read (ifinput, '(a80)') chmess
read (ifinput, '(a80)') chmess
write (6, *) ' height u_geo v_geo pgx pgy subs ' &
, ' dqtdx dqtdy dqtdtls thl_rad '
do k = kb, ke
read (ifinput, *) &
height(k), &
ug(k), &
vg(k), &
pgx(k), &
pgy(k), &
wfls(k), &
dqtdxls(k), &
dqtdyls(k), &
dqtdtls(k), &
thlpcar(k)
end do
close (ifinput)
do k = ke, kb, -1
write (6, '(3f7.1,5e12.4)') &
height(k), &
ug(k), &
vg(k), &
pgx(k), &
pgy(k), &
wfls(k), &
dqtdxls(k), &
dqtdyls(k), &
dqtdtls(k), &
thlpcar(k)
end do
end if ! end myid==0
! MPI broadcast variables read in
call MPI_BCAST(ug, kmax, MY_REAL, 0, comm3d, mpierr)
call MPI_BCAST(vg, kmax, MY_REAL, 0, comm3d, mpierr)
call MPI_BCAST(pgx, kmax, MY_REAL, 0, comm3d, mpierr)
call MPI_BCAST(pgy, kmax, MY_REAL, 0, comm3d, mpierr)
call MPI_BCAST(wfls, kmax, MY_REAL, 0, comm3d, mpierr)
call MPI_BCAST(dqtdxls, kmax, MY_REAL, 0, comm3d, mpierr)
call MPI_BCAST(dqtdyls, kmax, MY_REAL, 0, comm3d, mpierr)
call MPI_BCAST(dqtdtls, kmax, MY_REAL, 0, comm3d, mpierr)
call MPI_BCAST(thlpcar, kmax, MY_REAL, 0, comm3d, mpierr)
!-----------------------------------------------------------------
! 2.3 make large-scale horizontal pressure gradient
!-----------------------------------------------------------------
!******include rho if rho = rho(z) /= 1.0 ***********
if (lprofforc) then !tg3315
do k = kb, ke
dpdxl(k) = -pgx(k) - dpdx !-om23_gs*ug(k)-pgx(k)-dpdx
dpdyl(k) = -pgy(k)
end do
else
do k = kb, ke
!dpdxl(k) = om23_gs*vg(k)
!dpdyl(k) = -om23_gs*ug(k)
!dpdxl(k) = -ug(k)
!dpdyl(k) = -vg(k)
dpdxl(k) = om23_gs*vg(k) - pgx(k) - dpdx !corriolis forcing and pressure gradient
dpdyl(k) = -om23_gs*ug(k) - pgy(k)
end do
endif
!-----------------------------------------------------------------
! 2.4 large-scale subsidence, reintroduced ILS13 05.06.2014
!-----------------------------------------------------------------
whls(kb) = 0.0
do k = kb + 1, ke
whls(k) = (wfls(k)*dzf(k - 1) + wfls(k - 1)*dzf(k))/(2*dzh(k))
end do
whls(ke + 1) = (wfls(ke) + dzf(ke)*(wfls(ke) - wfls(ke - 1))/dzh(ke)) ! tg3315 31/07/18 removed a 0.5
! idtmax = floor(dtmax/tres)
btime = timee
! timeleft=ceiling(runtime/tres)
timeleft = runtime
dt_lim = timeleft
! write(6,*) 'real(dt)*tres= ',rdt, ' dtmax/100= ',dtmax/100
if ((lwarmstart) .or. (lstratstart)) then ! tg3315 to have cumulative number on restart files
read (startfile(6:13), '(i4)') ntrun
! ntrun = ichar(startfile(6:13))
else
ntrun = 0
end if
ntimee = nint(timee/dtmax)
tnextrestart = btime + trestart
tnextfielddump = btime + tfielddump
deallocate (height, th0av)
! call boundary
end subroutine readinitfiles