readrestartfiles Subroutine

public subroutine readrestartfiles()

Uses

  • proc~~readrestartfiles~~UsesGraph proc~readrestartfiles readrestartfiles module~initfac initfac proc~readrestartfiles->module~initfac module~modfields modfields proc~readrestartfiles->module~modfields module~modglobal modglobal proc~readrestartfiles->module~modglobal module~modinlet modinlet proc~readrestartfiles->module~modinlet module~modmpi modmpi proc~readrestartfiles->module~modmpi module~modsubgriddata modsubgriddata proc~readrestartfiles->module~modsubgriddata module~modsurfdata modsurfdata proc~readrestartfiles->module~modsurfdata module~initfac->module~modglobal module~initfac->module~modmpi mpi mpi module~initfac->mpi netcdf netcdf module~initfac->netcdf decomp_2d decomp_2d module~modfields->decomp_2d module~modinletdata modinletdata module~modinlet->module~modinletdata module~modinlet->mpi module~modmpi->mpi

Arguments

None

Calls

proc~~readrestartfiles~~CallsGraph proc~readrestartfiles readrestartfiles mpi_abort mpi_abort proc~readrestartfiles->mpi_abort mpi_bcast mpi_bcast proc~readrestartfiles->mpi_bcast proc~zinterpolate1d zinterpolate1d proc~readrestartfiles->proc~zinterpolate1d proc~zinterpolate2d zinterpolate2d proc~readrestartfiles->proc~zinterpolate2d proc~zinterpolatet1d zinterpolatet1d proc~readrestartfiles->proc~zinterpolatet1d proc~zinterpolatew1d zinterpolatew1d proc~readrestartfiles->proc~zinterpolatew1d

Called by

proc~~readrestartfiles~~CalledByGraph proc~readrestartfiles readrestartfiles proc~readinitfiles readinitfiles proc~readinitfiles->proc~readrestartfiles program~udales uDALES program~udales->proc~readinitfiles

Source Code

   subroutine readrestartfiles

      use modsurfdata, only:thls
      use modfields, only:u0, v0, w0, thl0, qt0, ql0, ql0h, qtav, qlav, e120, sv0, mindist, wall, &
         uav, vav, wav, uuav, vvav, wwav, uvav, uwav, vwav, svav, thlav, thl2av, sv2av, pres0, &
         svprof, viscratioav, thluav, thlvav, thlwav, svuav, svvav, svwav, presav, &
         uusgsav, vvsgsav, wwsgsav, uwsgsav, thlusgsav, thlwsgsav, svusgsav, svwsgsav, tkesgsav, &
         strain2av, nusgsav
      use modglobal, only:ib, ie, ih, jb, je, jh, kb, ke, kh, startfile, timee, totavtime, &
         ifinput, nsv, dt, cexpnr, lreadmean, lreadminl, &
         totinletav, lreadscal, ltempeq, dzf, numol, prandtlmoli, &
         lreadfacT, nfcts, nfaclyrs, lEB, iwallmom, iwalltemp, iwallmoist, ifacrestartversion, &
         tEB, tnextEB, lwarmstart
      use modmpi, only:cmyid, cmyidx, cmyidy, myid, comm3d, mpierr, my_real, mpi_logical
      use modsubgriddata, only:ekm
      use modinlet, only:zinterpolate1d, zinterpolatet1d, zinterpolatew1d, zinterpolate2d, &
         Uinl, Urec, Wrec, Utav, Tinl, Trec, &
         kbin, kein, lzinzsim, utaui, Ttav, ttaui
      use initfac, only:facT, facwsoil, facqsat, fachurel, facf, fachfi, facefi

      real, dimension(ib:ie, jb:je, kb:ke)  ::  dummy3d
      real, dimension(ib:ie, kbin:kein)    ::  Utavin
      real, dimension(ib:ie, kbin:kein)    ::  Ttavin
      real, dimension(kbin:kein)          ::  Uinlin
      real, dimension(kbin:kein)          ::  Urecin
      real, dimension(kbin:kein)          ::  Tinlin
      real, dimension(kbin:kein)          ::  Trecin
      real, dimension(kbin:kein + 1)        ::  Wrecin
      character(len(startfile)) :: name
      character(50) :: name2, name4
      character(50) :: namefac
      logical :: lfacrestart = .false., lscalrestart = .false.
      integer :: ifacversion, nfctsfile, nfaclyrsfile, iostatfac
      real :: timeefac, tEBfac, tnextEBfac
      real dummy
      integer i, j, k, n
      !********************************************************************

      !    1.0 Read initfiles
      !-----------------------------------------------------------------
      name = startfile
      name(5:5) = 'd'
      name(15:17) = cmyidx
      name(19:21) = cmyidy
      !write (6, *) 'loading ', name
      open (unit=ifinput, file=name, form='unformatted', status='old')

      read (ifinput) (((mindist(i, j, k), i=ib, ie), j=jb, je), k=kb, ke)
      read (ifinput) ((((wall(i, j, k, n), i=ib, ie), j=jb, je), k=kb, ke), n=1, 5)
      read (ifinput) (((u0(i, j, k), i=ib - ih, ie + ih), j=jb - jh, je + jh), k=kb, ke + kh)
      read (ifinput) (((v0(i, j, k), i=ib - ih, ie + ih), j=jb - jh, je + jh), k=kb, ke + kh)
      read (ifinput) (((w0(i, j, k), i=ib - ih, ie + ih), j=jb - jh, je + jh), k=kb, ke + kh)
      read (ifinput) (((pres0(i, j, k), i=ib - ih, ie + ih), j=jb - jh, je + jh), k=kb, ke + kh)
      read (ifinput) (((thl0(i, j, k), i=ib - ih, ie + ih), j=jb - jh, je + jh), k=kb, ke + kh)
      read (ifinput) (((e120(i, j, k), i=ib - ih, ie + ih), j=jb - jh, je + jh), k=kb, ke + kh)
      read (ifinput) (((ekm(i, j, k), i=ib - ih, ie + ih), j=jb - jh, je + jh), k=kb, ke + kh)
      read (ifinput) (((qt0(i, j, k), i=ib - ih, ie + ih), j=jb - jh, je + jh), k=kb, ke + kh)
      read (ifinput) (((ql0(i, j, k), i=ib - ih, ie + ih), j=jb - jh, je + jh), k=kb, ke + kh)
      read (ifinput) (((ql0h(i, j, k), i=ib - ih, ie + ih), j=jb - jh, je + jh), k=kb, ke + kh)
      read (ifinput) timee, dt
      close (ifinput)
      !write (6, *) 'finished loading ', name

      ! Scalar fields. On a warm start they are read from the inits restart set by
      ! default (lreadscal); the file must exist, there is no fallback. With
      ! lreadscal = .false. (or for lstratstart) sv0 keeps the scalar.inp profile that
      ! readinitfiles assigned before calling this routine, i.e. the cold-start value.
      if ((nsv > 0) .and. (lwarmstart) .and. (lreadscal)) then
         name(5:5) = 's'
         inquire (file=name, exist=lscalrestart)
         if (.not. lscalrestart) then
            write (0, *) 'ERROR: scalar restart file ', trim(name), ' not found.'
            write (0, *) '       A warm start with nsv > 0 and lreadscal = .true. needs the inits files written with the initd set.'
            write (0, *) '       Set lreadscal = .false. to start the scalars from scalar.inp.'//cexpnr//' instead.'
            call MPI_ABORT(MPI_COMM_WORLD, 1, mpierr)
         end if
         open (unit=ifinput, file=name, form='unformatted', status='old')
         read (ifinput) ((((sv0(i, j, k, n), i=ib - ih, ie + ih), j=jb - jh, je + jh), k=kb, ke + kh), n=1, nsv)
         read (ifinput) timee
         close (ifinput)
      end if

      ! Read the facet state written by writerestartfiles. The facet arrays are global
      ! rather than domain-decomposed, so the file has no rank fields and is read by
      ! rank 0 and broadcast. With lreadfacT, readfacetfiles (called before readinitfiles
      ! in program.f90) has not read Tfacinit*.inp, so this is the only source of the
      ! facet state: a missing, mismatched or truncated file aborts the run.
      if ((nfcts > 0) .and. (lwarmstart) .and. (lreadfacT) .and. &
          ((lEB) .or. (iwalltemp == 2) .or. (iwallmom == 2) .or. (iwallmoist == 2))) then
         namefac = 'initf        .'
         namefac(6:13) = startfile(6:13)
         namefac(15:17) = cexpnr

         if (myid == 0) then
            inquire (file=trim(namefac), exist=lfacrestart)
         end if
         call MPI_BCAST(lfacrestart, 1, MPI_LOGICAL, 0, comm3d, mpierr)

         if (lfacrestart) then
            if (myid == 0) then
               open (unit=ifinput, file=trim(namefac), form='unformatted', status='old')
               read (ifinput, iostat=iostatfac) ifacversion, nfctsfile, nfaclyrsfile
               if (iostatfac /= 0) then
                  write (0, *) 'ERROR: cannot read the header of facet restart file ', trim(namefac)
                  call MPI_ABORT(MPI_COMM_WORLD, 1, mpierr)
               end if
               if (ifacversion /= ifacrestartversion) then
                  write (0, *) 'ERROR: facet restart file ', trim(namefac), ' has format version ', ifacversion, &
                     ' but this build reads version ', ifacrestartversion
                  call MPI_ABORT(MPI_COMM_WORLD, 1, mpierr)
               end if
               if ((nfctsfile /= nfcts) .or. (nfaclyrsfile /= nfaclyrs)) then
                  write (0, *) 'ERROR: facet restart file ', trim(namefac), ' does not match this simulation'
                  write (0, *) '       file       : nfcts = ', nfctsfile, ' nfaclyrs = ', nfaclyrsfile
                  write (0, *) '       namoptions : nfcts = ', nfcts, ' nfaclyrs = ', nfaclyrs
                  call MPI_ABORT(MPI_COMM_WORLD, 1, mpierr)
               end if
               read (ifinput, iostat=iostatfac) ((facT(n, j), n=1, nfcts), j=1, nfaclyrs + 1)
               if ((iostatfac == 0) .and. (lEB)) read (ifinput, iostat=iostatfac) (facwsoil(n), n=1, nfcts)
               if (iostatfac == 0) read (ifinput, iostat=iostatfac) (facqsat(n), n=1, nfcts)
               if (iostatfac == 0) read (ifinput, iostat=iostatfac) (fachurel(n), n=1, nfcts)
               if (iostatfac == 0) read (ifinput, iostat=iostatfac) ((facf(n, j), n=1, nfcts), j=1, 5)
               if (iostatfac == 0) read (ifinput, iostat=iostatfac) timeefac
               ! the energy-balance interval in progress at the checkpoint: time of the last
               ! and next EB update and the fluxes accumulated since the last one (rank-0
               ! arrays). Restoring it continues that interval instead of restarting it at
               ! btime, which would drop the accumulated flux.
               if ((iostatfac == 0) .and. (lEB)) then
                  read (ifinput, iostat=iostatfac) tEBfac, tnextEBfac
                  if (iostatfac == 0) read (ifinput, iostat=iostatfac) (fachfi(n), n=1, nfcts)
                  if (iostatfac == 0) read (ifinput, iostat=iostatfac) (facefi(n), n=1, nfcts)
               end if
               close (ifinput)
               if (iostatfac /= 0) then
                  write (0, *) 'ERROR: facet restart file ', trim(namefac), ' is truncated or corrupt (iostat ', &
                     iostatfac, '); the facet state cannot be restored'
                  call MPI_ABORT(MPI_COMM_WORLD, 1, mpierr)
               end if
               if (lEB) then
                  tEB = tEBfac
                  tnextEB = tnextEBfac
                  lEBclockrestored = .true.
               end if
               write (*, *) 'Restored facet state from ', trim(namefac), ' (format version ', ifacversion, &
                  ', timee ', timeefac, ')'
               if (lEBclockrestored) then
                  write (*, *) 'Energy-balance interval restored from restart: tEB = ', tEB, &
                     ' tnextEB = ', tnextEB
               end if
            end if

            ! facwsoil is allocated on rank 0 only and used only there, so it is not broadcast
            call MPI_BCAST(facT(1:nfcts, 1:nfaclyrs + 1), nfcts*(nfaclyrs + 1), MY_REAL, 0, comm3d, mpierr)
            call MPI_BCAST(facqsat(1:nfcts), nfcts, MY_REAL, 0, comm3d, mpierr)
            call MPI_BCAST(fachurel(1:nfcts), nfcts, MY_REAL, 0, comm3d, mpierr)
            call MPI_BCAST(facf(1:nfcts, 1:5), nfcts*5, MY_REAL, 0, comm3d, mpierr)
            ! tEB and the accumulators are used on rank 0 only; tnextEB is tested on every rank
            call MPI_BCAST(lEBclockrestored, 1, MPI_LOGICAL, 0, comm3d, mpierr)
            if (lEBclockrestored) call MPI_BCAST(tnextEB, 1, MY_REAL, 0, comm3d, mpierr)
         else
            ! No fallback: readfacetfiles did not read Tfacinit*.inp for this warm start,
            ! so without the file there is no facet state at all.
            if (myid == 0) then
               write (0, *) 'ERROR: facet restart file ', trim(namefac), ' not found.'
               write (0, *) '       A warm start with lreadfacT = .true. needs the initf file written with the initd set.'
               write (0, *) '       Set lreadfacT = .false. to re-seed the facet temperatures from Tfacinit.inp.'//cexpnr
            end if
            call MPI_ABORT(MPI_COMM_WORLD, 1, mpierr)
         end if
      end if

      ! read mean variables if asked for by lreadmean
      name2 = 'means   .'
      name2(6:8) = cmyid
      name2(10:12) = cexpnr
      if (lreadmean) then
         write (6, *) 'Reading meansXXX.XXX, proc = ', myid
         open (unit=ifinput, file=name2, form='unformatted')
         read (ifinput) totavtime, nsv
         read (ifinput) (((uav(i, j, k), i=ib - ih, ie + ih), j=jb - jh, je + jh), k=kb - kh, ke + kh)
         read (ifinput) (((vav(i, j, k), i=ib - ih, ie + ih), j=jb - jh, je + jh), k=kb - kh, ke + kh)
         read (ifinput) (((wav(i, j, k), i=ib - ih, ie + ih), j=jb - jh, je + jh), k=kb - kh, ke + kh)
         read (ifinput) (((thlav(i, j, k), i=ib - ih, ie + ih), j=jb - jh, je + jh), k=kb - kh, ke + kh)
         read (ifinput) (((qtav(i, j, k), i=ib - ih, ie + ih), j=jb - jh, je + jh), k=kb - kh, ke + kh)
         read (ifinput) (((qlav(i, j, k), i=ib - ih, ie + ih), j=jb - jh, je + jh), k=kb - kh, ke + kh)
         read (ifinput) (((presav(i, j, k), i=ib - ih, ie + ih), j=jb - jh, je + jh), k=kb - kh, ke + kh)
         read (ifinput) ((((svav(i, j, k, n), i=ib - ih, ie + ih), j=jb - jh, je + jh), k=kb - kh, ke + kh), n=1, nsv)
         read (ifinput) (((viscratioav(i, j, k), i=ib - ih, ie + ih), j=jb - jh, je + jh), k=kb - kh, ke + kh)
         read (ifinput) (((uuav(i, j, k), i=ib - ih, ie + ih), j=jb - jh, je + jh), k=kb - kh, ke + kh)
         read (ifinput) (((vvav(i, j, k), i=ib - ih, ie + ih), j=jb - jh, je + jh), k=kb - kh, ke + kh)
         read (ifinput) (((wwav(i, j, k), i=ib - ih, ie + ih), j=jb - jh, je + jh), k=kb - kh, ke + kh)
         read (ifinput) (((thl2av(i, j, k), i=ib - ih, ie + ih), j=jb - jh, je + jh), k=kb - kh, ke + kh)
         read (ifinput) ((((sv2av(i, j, k, n), i=ib - ih, ie + ih), j=jb - jh, je + jh), k=kb - kh, ke + kh), n=1, nsv)
         read (ifinput) (((uvav(i, j, k), i=ib, ie + ih), j=jb, je + jh), k=kb, ke)
         read (ifinput) (((uwav(i, j, k), i=ib, ie + ih), j=jb, je), k=kb, ke + kh)
         read (ifinput) (((vwav(i, j, k), i=ib, ie), j=jb, je + jh), k=kb, ke + kh)
         read (ifinput) (((thluav(i, j, k), i=ib, ie), j=jb, je), k=kb, ke)
         read (ifinput) (((thlvav(i, j, k), i=ib, ie), j=jb, je + jh), k=kb, ke)
         read (ifinput) (((thlwav(i, j, k), i=ib, ie), j=jb, je), k=kb, ke + kh)
         read (ifinput) ((((svuav(i, j, k, n), i=ib, ie), j=jb, je), k=kb, ke), n=1, nsv)
         read (ifinput) ((((svvav(i, j, k, n), i=ib, ie), j=jb, je + jh), k=kb, ke), n=1, nsv)
         read (ifinput) ((((svwav(i, j, k, n), i=ib, ie), j=jb, je), k=kb, ke + kh), n=1, nsv)
         close (ifinput)
         write (6, *) 'Total averaging time so far: ', totavtime

         ! read <x'y'>_SGS to file.
         name2 = 'SGS__   .'
         name2(6:8) = cmyid
         name2(10:12) = cexpnr
         open (unit=ifinput, file=name2, form='unformatted')
         read (ifinput) dummy, dummy
         read (ifinput) (((uusgsav(i, j, k), i=ib - ih, ie + ih), j=jb - jh, je + jh), k=kb - kh, ke + kh)
         read (ifinput) (((vvsgsav(i, j, k), i=ib - ih, ie + ih), j=jb - jh, je + jh), k=kb - kh, ke + kh)
         read (ifinput) (((wwsgsav(i, j, k), i=ib - ih, ie + ih), j=jb - jh, je + jh), k=kb - kh, ke + kh)
         read (ifinput) (((uwsgsav(i, j, k), i=ib, ie + ih), j=jb, je), k=kb, ke + kh)
         read (ifinput) (((dummy3d(i, j, k), i=ib, ie), j=jb, je), k=kb, ke) ! this is dissresav, which will be computed using other mean quantities
         read (ifinput) (((tkesgsav(i, j, k), i=ib, ie), j=jb, je), k=kb, ke)
         read (ifinput) (((dummy3d(i, j, k), i=ib, ie), j=jb, je), k=kb, ke) ! this is disssgsav, which will be computed using other mean quantities
         read (ifinput) (((strain2av(i, j, k), i=ib, ie), j=jb, je), k=kb, ke) ! <SijSij> (NOT <Sij><Sij> !!) (average over time)
         read (ifinput) (((nusgsav(i, j, k), i=ib, ie), j=jb, je), k=kb, ke) ! <nu_sgs> (average over time)
         read (ifinput) (((thlusgsav(i, j, k), i=ib, ie + ih), j=jb, je), k=kb, ke)
         read (ifinput) (((thlwsgsav(i, j, k), i=ib, ie), j=jb, je), k=kb, ke + kh)
         read (ifinput) ((((svusgsav(i, j, k, n), i=ib, ie + ih), j=jb, je), k=kb, ke), n=1, nsv)
         read (ifinput) ((((svwsgsav(i, j, k, n), i=ib, ie), j=jb, je), k=kb, ke + kh), n=1, nsv)
         close (ifinput)

      end if

      ! read mean profiles for inlet generator
      if (lreadminl) then
         if (.not. lzinzsim) then
            name4 = 'meaninlet.   '
            name4(11:13) = cexpnr
            open (unit=ifinput, file=name4, form='unformatted')
            read (ifinput) totinletav ! interval of time-average
            read (ifinput) (Uinlin(k), k=kbin, kein)
            read (ifinput) (Urecin(k), k=kbin, kein)
            read (ifinput) (Wrecin(k), k=kbin, kein + 1)
            read (ifinput) ((Utavin(i, k), i=ib, ie), k=kbin, kein)
            close (ifinput)

            call zinterpolate1d(Uinlin, Uinl) ! interpolate inlet profile to zgrid
            call zinterpolate1d(Urecin, Urec)
            call zinterpolatew1d(Wrecin, Wrec)
            call zinterpolate2d(Utavin, Utav)

            if (ltempeq) then
               name4 = 'tempinlet.   '
               name4(11:13) = cexpnr
               open (unit=ifinput, file=name4, form='unformatted')
               read (ifinput) totinletav ! interval of time-average
               read (ifinput) (Tinlin(k), k=kbin, kein)
               read (ifinput) (Trecin(k), k=kbin, kein)
               read (ifinput) ((Ttavin(i, k), i=ib, ie), k=kbin, kein)
               close (ifinput)

               call zinterpolatet1d(Tinlin, Tinl)
               call zinterpolatet1d(Trecin, Trec)
               call zinterpolate2d(Ttavin, Ttav)
            end if ! ltempeq

         else !lzinzsim=.true. -> inlet grid equals sim grid
            name4 = 'meaninlet.   '
            name4(11:13) = cexpnr
            open (unit=ifinput, file=name4, form='unformatted')
            read (ifinput) totinletav ! interval of time-average
            read (ifinput) (Uinl(k), k=kb, ke)
            read (ifinput) (Urec(k), k=kb, ke)
            read (ifinput) (Wrec(k), k=kb, ke + 1)
            read (ifinput) ((Utav(i, k), i=ib, ie), k=kb, ke)

            close (ifinput)

            if (ltempeq) then
               name4 = 'tempinlet.   '
               name4(11:13) = cexpnr
               open (unit=ifinput, file=name4, form='unformatted')
               read (ifinput) totinletav ! interval of time-average
               read (ifinput) (Tinl(k), k=kb, ke)
               read (ifinput) (Trec(k), k=kb, ke)
               read (ifinput) ((Ttav(i, k), i=ib, ie), k=kb, ke)

               close (ifinput)
            end if ! ltempeq
         end if ! lzinzsim

         utaui = sqrt(abs(2*numol*Uinl(kb)/dzf(kb))) ! average streamwise friction at inlet (need for first time step)
         if (ltempeq) then
            ttaui = numol*prandtlmoli*2.*(Tinl(kb) - thls)/(dzf(kb)*utaui)
         end if
      end if !(lreadminl)

   end subroutine readrestartfiles