ReadFlowField.F90

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                         ! 
!    FILE: inirea.F90                                     !
!    CONTAINS: subroutine ReadFlowField                   !
!                                                         ! 
!    PURPOSE: Initialization routine. Reads in a flow     !
!     snapshot generated by a previous simulation, and    !
!     interpolates to the current grid if necessary       !
!                                                         !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      subroutine ReadFlowField
      use mpih
      use decomp_2d
      use local_arrays
      use param
      implicit none
      character*70 :: filnam,dsetname
      integer :: nzom,nyom,nxo,nyo,nzo
      integer :: xs2og,xe2og,xs3og,xe3og
      integer :: istro3
      integer :: vxflag,vyflag,vzflag,teflag
!     integer (kind=MPI_ADDRESS_KIND) :: extent,lb
      real :: stro3
      real, allocatable, dimension(:,:,:) :: tempold,vyold,vxold,vzold
      logical :: fexist
      
!EP   Reading old grid information by master
      filnam = trim('continua_master.h5')

      inquire(file=filnam,exist=fexist)
      if (.not.fexist) then 
        write(*,*) 'Continuation files not found'
        call MpiAbort
      end if

      if (ismaster) then
       dsetname = trim('nx')
       call HdfSerialReadIntScalar(dsetname,filnam,nxo)
       dsetname = trim('ny')
       call HdfSerialReadIntScalar(dsetname,filnam,nyo)
       dsetname = trim('nz')
       call HdfSerialReadIntScalar(dsetname,filnam,nzo)
       dsetname = trim('time')
       call HdfSerialReadRealScalar(dsetname,filnam,time)
       dsetname = trim('istr3')
       call HdfSerialReadIntScalar(dsetname,filnam,istro3)
       dsetname = trim('str3')
       call HdfSerialReadRealScalar(dsetname,filnam,stro3)
      endif
      
      call MpiBarrier
      call MpiBcastInt(nxo)
      call MpiBcastInt(nyo)
      call MpiBcastInt(nzo)
      call MpiBcastInt(istro3)
      call MpiBcastReal(stro3)
      call MpiBcastReal(time)

!RO   Codes for passing to other routines which variable we are dealing with 
      vxflag = 1
      vyflag = 2
      vzflag = 3
      teflag = 4
      
!EP   Check whether grid specifications have been updated
      if(nyo.ne.ny.or.nxo.ne.nx.or.nzo.ne.nz &
       .or.istro3.ne.istr3.or.(abs(stro3-str3).gt.1e-8)) then

       if(ismaster) write(*,*) "Interpolating new grid"

!RO    Using more than double the amount of points results in crashes

       if(nz.gt.nzo*2.or.ny.gt.nyo*2.or.nx.gt.nxo*2) then

        if(ismaster) then 
          write(*,*) "New grid resolution cannot be more ",&
          "than twice the old resolution"
        end if

        call MpiAbort

       endif

!RO   Interpolate

      nyom = nyo - 1
      nzom = nzo - 1
      
      call MpiBarrier

      xs2og = floor(real(xstart(2)*nyom/nym))
      xe2og = ceiling(real(xend(2)*nyom/nym))
      xs3og = floor(real(xstart(3)*nzom/nzm))
      xe3og = ceiling(real(xend(3)*nzom/nzm))

      xs2og = max(xs2og,1)
      xe2og   = min(xe2og,nyom)
      xs3og = max(xs3og,1)
      xe3og   = min(xe3og,nzom)
      
!RO   I think this is deprecated, must check
!     lb=0
!     call MPI_TYPE_GET_EXTENT(MPI_DOUBLE_PRECISION,lb,extent,ierr)


!EP   vz
      allocate(vzold(0:nxo+1,xs2og-lvlhalo:xe2og+lvlhalo, &
        xs3og-lvlhalo:xe3og+lvlhalo))
      
      call HdfReadContinua(nzo,nyo,nxo,xs2og,xe2og, &
       xs3og,xe3og,vzflag,vzold(1:nxo,xs2og-lvlhalo:xe2og+lvlhalo, &
       xs3og-lvlhalo:xe3og+lvlhalo))

      call interp(vzold,vz(1:nx,xstart(2):xend(2),xstart(3):xend(3)) &
       ,nzo,nyo,nxo,istro3,stro3,vzflag,xs2og,xe2og,xs3og,xe3og)

      deallocate(vzold)

!EP   vy
      allocate(vyold(0:nxo+1,xs2og-lvlhalo:xe2og+lvlhalo, &
     & xs3og-lvlhalo:xe3og+lvlhalo))
      
      call HdfReadContinua(nzo,nyo,nxo,xs2og,xe2og, &
     & xs3og,xe3og,vyflag,vyold(1:nxo,xs2og-lvlhalo:xe2og+lvlhalo, &
     & xs3og-lvlhalo:xe3og+lvlhalo))

      call interp(vyold,vy(1:nx,xstart(2):xend(2),xstart(3):xend(3)) &
     & ,nzo,nyo,nxo,istro3,stro3,vyflag,xs2og,xe2og,xs3og,xe3og)

      deallocate(vyold)

!EP   vx
      allocate(vxold(0:nxo+1,xs2og-lvlhalo:xe2og+lvlhalo, &
     & xs3og-lvlhalo:xe3og+lvlhalo))
      
      call HdfReadContinua(nzo,nyo,nxo,xs2og,xe2og, &
     & xs3og,xe3og,vxflag,vxold(1:nxo,xs2og-lvlhalo:xe2og+lvlhalo,  &
     & xs3og-lvlhalo:xe3og+lvlhalo))

      call interp(vxold,vx(1:nx,xstart(2):xend(2),xstart(3):xend(3)) &
     & ,nzo,nyo,nxo,istro3,stro3,vxflag,xs2og,xe2og,xs3og,xe3og)

      deallocate(vxold)

!EP   temp
      allocate(tempold(0:nxo+1,xs2og-lvlhalo:xe2og+lvlhalo, &
     & xs3og-lvlhalo:xe3og+lvlhalo))
      
      call HdfReadContinua(nzo,nyo,nxo,xs2og,xe2og, &
       xs3og,xe3og,teflag,tempold(1:nxo,xs2og-lvlhalo:xe2og+lvlhalo, &
       xs3og-lvlhalo:xe3og+lvlhalo))

      call interp(tempold,temp(1:nx,xstart(2):xend(2),xstart(3):xend(3)) &
       ,nzo,nyo,nxo,istro3,stro3,teflag,xs2og,xe2og,xs3og,xe3og)

      deallocate(tempold)


      else

!EP   Normal one to one HDF read (no interpolation)
      call HdfReadContinua(nz,ny,nx,xstart(2),xend(2) &
     & ,xstart(3),xend(3),vzflag,vz)
      call HdfReadContinua(nz,ny,nx,xstart(2),xend(2) &
     & ,xstart(3),xend(3),vyflag,vy)
      call HdfReadContinua(nz,ny,nx,xstart(2),xend(2) &
     & ,xstart(3),xend(3),vxflag,vx)
      call HdfReadContinua(nz,ny,nx,xstart(2),xend(2) &
     & ,xstart(3),xend(3),teflag,temp)

      endif

      if (resetlogstime) time=0.

!EP   Increase the maximum simulation time by the end time in the
!continuation files

       tmax = tmax + time

      return                                                            
      end                                                               
!