readinitfiles Subroutine

   subroutine readinitfiles
      use modfields, only:u0, v0, w0, um, vm, wm, thlm, thl0, thl0h, qtm, qt0, qt0h, uinit, vinit, &
         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, u0h, thl0c
            use modglobal,         only : ib,ie,ih,ihc,jb,je,jh,jhc,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, imax, xh, xf, linoutflow, lper2inout, iinletgen, lreadminl, &
         uflowrate, vflowrate,ltempeq, prandtlmoli, freestreamav, &
         tnextfielddump, tfielddump, tsample, tstatsdump, startfile, lprofforc, lchem, k1, JNO2,&
         idriver,dtdriver,driverstore,tdriverstart,tdriverstart_cold,tdriverdump,lchunkread,xlen,ylen,itot,jtot,ibrank,ierank,jbrank,jerank,BCxm,BCym,lrandomize,BCxq,BCxs,BCxT, BCyq,BCys,BCyT,BCxm_driver,&
         tEB,tnextEB,dtEB,BCxs_custom,lEB,lfacTlyrs,tfac,tnextfac,dtfac
      use modsubgriddata, only:ekm, ekh, loneeqn
      use modsurfdata, only:wtsurf, wqsurf, wsvsurf, &
         thls, thvs, ps, qts, svs, sv_top
      ! use modsurface,        only : surface,dthldz
      use modboundary, only:boundary, tqaver, halos
      use modmpi, only:slabsum, myid, comm3d, mpierr, my_real, avexy_ibm, myidx, myidy
      use modthermodynamics, only:thermodynamics, calc_halflev
      use modinletdata, only:Uinl, Urec, Wrec, u0inletbc, v0inletbc, w0inletbc, ubulk, vbulk, irecy, Utav, Ttav, &
         uminletbc, vminletbc, wminletbc, u0inletbcold, v0inletbcold, w0inletbcold, &
         storeu0inletbc, storev0inletbc, storew0inletbc, nstepread, nfile, Tinl, &
         Trec, tminletbc, t0inletbcold, t0inletbc, storet0inletbc, utaui, ttaui, iangle,&
         u0driver,umdriver,v0driver,vmdriver,w0driver,e120driver,tdriver,thl0driver,qt0driver,storetdriver,&
         storeu0driver,storeumdriver,storev0driver,storew0driver,storee120driver,storethl0driver,storeqt0driver,&
         nstepreaddriver
      use modinlet, only:readinletfile
      use moddriver, only: readdriverfile,initdriver,drivergen,readdriverfile_chunk
      use decomp_2d, only : exchange_halo_z, update_halo, decomp_main

      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, dimension(kb:ke+kh)  :: u_init, v_init, thl_init, qt_init
      real tv, ran, ran1

      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)') &
               !    write (*, *) &
               !       height(k), &
               !       thlprof(k), &
               !       qtprof(k), &
               !       uprof(k), &
               !       vprof(k), &
               !       e12prof(k)
               !
               ! end do

               if (loneeqn) then
                 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 ! 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

            do k = kb, ke
               do j = jb - jhc, je + jhc
                  do i = ib - ihc, ie + ihc
                     thl0c(i, j, k) = thlprof(k)
                  end do
               end do
            end do

            ! if (ibrank) then
            ! do j=jb-1,je+1
            !   um(ib,j,kb:ke) = uprof
            !   um(ib-1,j,kb:ke) = uprof
            ! end do
            ! end if
            !
            ! if (ierank) then
            ! do j=jb-1,je+1
            !   !um(ie,j,kb:ke) = uprof
            !   um(ie+1,j,kb:ke) = uprof
            ! end do
            ! end if


            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

            if (lrandomize) then
              !! 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
            end if

            !       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

            ! SO: Manually override fields

            ! if (((BCxm == 1) .and. (BCym == 1)) .or. ((BCxm == 6) .and. (BCym == 6))) then
            !   ! TGV (assumes equidistant x grid)
            !   do i = ib-1,ie+1
            !     do j = jb-1,je+1
            !       do k = kb-1,ke+1
            !         um(i,j,k) = 1. * sin(4.*atan(1.) * 2. * (dx*((i-1)+myidx*imax)) / ylen) &
            !         * cos(4.*atan(1.) * 2. * (dy*((0.5+(j-1))+myidy*jmax)) / ylen) !&
            !         !* cos(4.*atan(1.) * 2. * zf(k) / ylen)
            !         vm(i,j,k) = 1. *-cos(4.*atan(1.) * 2. * (dx*((0.5+(i-1))+myidx*imax)) / ylen)  &
            !         * sin(4.*atan(1.) * 2. * (dy*((j-1)+myidy*jmax)) / ylen) !&
            !         !* cos(4.*atan(1.) * 2. * zf(k) / ylen)
            !         wm(i,j,k) = 0.
            !       end do
            !     end do
            !   end do
            ! end if

            ! ! For shear case
            ! um(:,:,ke+1) = um(:,:,ke)
            ! um(ib-1,:,:) = um(ib,:,:)
            ! um(ie-1,:,:) = um(ie,:,:)
            ! um(:,jb-1,:) = um(:,jb,:)
            ! um(:,je+1,:) = um(:,je,:)

            u0 = um
            v0 = vm
            w0 = wm

            call halos

            uinit = um
            vinit = vm

            ! ! zeros
            ! do i = ib,ie
            !   do j = jb,je
            !     do k = kb,ke
            !       um(i,j,k) = 0.
            !       u0(i,j,k) = um(i,j,k)
            !       uinit(i,j,k) = um(i,j,k)
            !       vm(i,j,k) = 0.
            !       v0(i,j,k) = vm(i,j,k)
            !       vinit(i,j,k) = vm(i,j,k)
            !       wm(i,j,k) = 0.
            !       w0(i,j,k) = 0.
            !     end do
            !   end do
            ! end do

            ! ! ones
            ! do i = ib,ie
            !   do j = jb,je
            !     do k = kb,ke
            !       um(i,j,k) = 1.
            !       u0(i,j,k) = um(i,j,k)
            !       uinit(i,j,k) = um(i,j,k)
            !       vm(i,j,k) = 1.
            !       v0(i,j,k) = vm(i,j,k)
            !       vinit(i,j,k) = vm(i,j,k)
            !       wm(i,j,k) = 0.
            !       w0(i,j,k) = 0.
            !     end do
            !   end do
            ! end do

            ! do i = ib,ie
            !   do j = jb,je
            !     do k = kb,ke
            !       um(i,j,k) = (i + myidx*imax) * (j + myidy*jmax)
            !       u0(i,j,k) = (i + myidx*imax) * (j + myidy*jmax)
            !       vm(i,j,k) = 0.
            !       v0(i,j,k) = 0.
            !       wm(i,j,k) = 0.
            !       w0(i,j,k) = 0.
            !     end do
            !   end do
            ! end do

            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
            !write (6, *) 'Modstartup: ubulk=', ubulk
            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

               if (ibrank) then
                  if (lchunkread) then
                      call readdriverfile_chunk
                  else
                     call readdriverfile
                  end if
                  call drivergen
               end if

               ! do k = kb, ke
               !    do j = jb-1, je+1
               !       do i = ib-1, ie+1
               !          u0(i, j, k) = u0driver(j, k)
               !          um(i, j, k) = umdriver(j, k)
               !          v0(i, j, k) = v0driver(j, k)
               !          vm(i, j, k) = vmdriver(j, k)
               !       end do
               !    end do
               ! end do

               ! 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

               if (runtime < tdriverstart) then
                  if (myid==0) write(*,*) 'Warning! No driver files will be written as runtime < tdriverstart.'
               else
                  if (trestart /= (tdriverstart + (driverstore-1)*dtdriver)) then
                     trestart = (tdriverstart + (driverstore-1)*dtdriver)
                  end if
                  if (myid==0) then
                     write(*,'(A,F15.5)') 'Warning! for driver simulation, trestart gets set as &
                                           (tdriverstart + (driverstore-1)*dtdriver), ignoring the &
                                           trestart mentioned in namoptions. Hence, trestart = ',(tdriverstart + (driverstore-1)*dtdriver)
                     if (runtime >= tdriverstart .and. runtime+1e-10 < (tdriverstart + (driverstore-1)*dtdriver)) then
                        write(*,*) 'Warning! Driver files cannot be written upto ', driverstore, ' steps. &
                                    &Consider taking runtime >= (tdriverstart + (driverstore-1)*dtdriver).'
                     end if
                  end if
               end if

              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)

            if (BCxs /= BCxs_custom) then
               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
            end if

            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.
         ! SO: can't do this yet because uses u0 and it does not include halo cells
         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

            ! average initial profiles
            call avexy_ibm(u_init(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 avexy_ibm(v_init(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 avexy_ibm(thl_init(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 avexy_ibm(qt_init(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.)

            if (myid == 0) then
               ! Read profiles from file (potentially for forcing)
               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 initial profile
               open (ifinput, file='prof_restart.'//cexpnr)
               write (ifinput, *) '# SDBL flow'
               write (ifinput, *) '# z thl qt u v e12'
               do k = kb, ke
                  write (ifinput, '(f20.15,5f12.6)') &
                     height(k), &
                     thl_init(k), &
                     qt_init(k), &
                     u_init(k), &
                     v_init(k), &
                     e12prof(k)
               end do
               close (ifinput)

               ! 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

               if (loneeqn) then
                 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 ! 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)
            btime = timee
            um = u0
            vm = v0
            wm = w0
            thlm = thl0
            qtm = qt0
            svm = sv0 ! What if nsv=0?
            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
               if (timee>=tdriverstart) then

                  tdriverstart_cold = tdriverstart
                  tdriverstart = timee
                  if (trestart /= (driverstore-1)*dtdriver) then
                     trestart = (driverstore-1)*dtdriver
                  end if

                  if (myid==0) then
                     write(*,'(A,F15.5)') "Warning! during warmstart of driver simulat ion, tdriverstart &
                                           &gets overwritten by the time instant of initd restartfile, ignoring the &
                                           &tdriverstart mentioned in namoptions. Hence, tdriverstart = ",timee
                     write(*,'(A,F15.5)') 'Warning! for this driver simulation, trestart gets set as &
                                           (driverstore-1)*dtdriver, ignoring the trestart mentioned &
                                           in namoptions. Hence, trestart = ',(driverstore-1)*dtdriver
                     if ( runtime < (driverstore-1)*dtdriver ) then
                        write(*,*) 'Warning! Driver files cannot be written upto ', driverstore, ' steps. &
                                    &Consider taking runtime >= (driverstore-1)*dtdriver).'
                     end if
                  end if

               else ! if (timee<tdriverstart)
                  if (trestart /= (tdriverstart + (driverstore-1)*dtdriver - btime)) then
                     trestart = (tdriverstart + (driverstore-1)*dtdriver) - btime
                  end if
                  if (myid==0) then
                     write(*,'(A,F15.5)') 'Warning! for this driver simulation, trestart gets set as &
                                           (tdriverstart + (driverstore-1)*dtdriver - btime), ignoring the &
                                           trestart mentioned in namoptions. Hence, trestart = ',(tdriverstart + (driverstore-1)*dtdriver) - btime
                     if ( (timee + runtime) < (tdriverstart + (driverstore-1)*dtdriver) ) then
                        write(*,*) 'Warning! Driver files cannot be written upto ', driverstore, ' steps. &
                                    &Consider taking runtime + ',timee,' >= (tdriverstart + (driverstore-1)*dtdriver).'
                     end if
                  end if
               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

               if (ibrank) then
                  if (lchunkread) then
                      call readdriverfile_chunk
                  else
                     call readdriverfile
                  end if
                  call drivergen
               end if

              !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

            if (lEB .and. (lfacTlyrs .eqv. .false.)) then
               if (myid==0) write(*,*) "Warmstarting an EB simulation - consider setting internal facet temperatures"
            end if
         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), '(i8)') ntrun
         ! ntrun = ichar(startfile(6:13))
      else
         ntrun = 0
      end if
      ntimee = nint(timee/dtmax)
      tnextrestart = btime + trestart
      tnextfielddump = btime + tfielddump
      tEB = btime
      tnextEB = btime + dtEB
      tfac = btime
      tnextfac = btime + dtfac
      deallocate (height, th0av)

      !    call boundary

   end subroutine readinitfiles