initglobal Subroutine

   subroutine initglobal
      use modmpi,   only : myid, comm3d, my_real, mpierr
      use decomp_2d
      implicit none

      integer :: advarr(4)
      real phi, colat, silat, omega, omega_gs
      integer :: i, j, k, n
      character(80) chmess

      !timestepping
      if (courant < 0) then
         select case (iadv_mom)
         case (iadv_cd2)
            courant = 1.5
         case default
            courant = 1.4
         end select
         if (any(iadv_sv(1:nsv) == iadv_kappa) .or. any((/iadv_thl, iadv_qt, iadv_tke/) == iadv_kappa)) then
            courant = min(courant, 1.1)
         elseif (any(iadv_sv(1:nsv) == iadv_upw) .or. any((/iadv_thl, iadv_qt, iadv_tke/) == iadv_upw)) then
            courant = min(courant, 1.1)
         elseif (any(iadv_sv(1:nsv) == iadv_cd2) .or. any((/iadv_thl, iadv_qt, iadv_tke/) == iadv_cd2)) then
            courant = min(courant, 1.5)
         end if
      end if

      ! phsgrid

      !jmax = jtot/nprocy ! Only in z-pencil and not true generally - uneven no.
      !imax = itot/nprocx ! Only in z-pencil
      isen = imax ! Only in z-pencil - replace eventually so it is pencil-independent (in poisson)
      jsen = jmax ! Only in z-pencil - replace eventually so it is pencil-independent (in poisson)
      !set the number of ghost cells. NB: This switch has to run in order of required ghost cells
      advarr = (/iadv_mom, iadv_tke, iadv_thl, iadv_qt/)
      if (any(advarr == iadv_kappa)) then
         ih = 2
         jh = 2
         kh = 1
         ! SO: think this is inconsistent
      elseif (any(advarr == iadv_cd2) .or. any(iadv_sv(1:nsv) == iadv_cd2)) then
         ih = 1
         jh = 1
         kh = 1
         ihc = 1
         jhc = 1
         khc = 1
      end if

      ! J. Tomas added this for using only kappa scheme for sv(:)
      if (any(iadv_sv(1:nsv) == iadv_kappa) .or. (iadv_thl == iadv_kappa)) then
         ih = 1
         jh = 1
         kh = 1
         ihc = 2
         jhc = 2
         khc = 2
      end if

      ! Eventually ib etc should be completely replaced.
      ! All arrays start at 1, like in 2DECOMP, and end at e.g. zsize(1) = imax in old terminology
      ib = 1    ! Remove eventually
      jb = 1
      jgb = jb ! global j range (starting at the same as j as the processor j range)
      jge = jtot ! global j range
      kb = 1 ! Make redundant

      !kmax = ktot

      ! Define indices in terms of 2DECOMP's. Subject to change! z=pencil 'special' for now, but could rename e.g. imax -> imax3
      imax1 = xsize(1) !=itot
      imax2 = ysize(1)
      imax = zsize(1)

      jmax1 = xsize(2)
      jmax2 = ysize(2) !=jtot
      jmax = zsize(2)

      kmax1 = xsize(3)
      kmax2 = ysize(3)
      kmax = zsize(3)

      ie = imax
      je = jmax
      ke = kmax

      decomp_main%zlevel = (/ih, jh, kh/)

      if (zstart(1) == 1) then
        ibrank = .true.
      else
        ibrank = .false.
      end if

      if (zend(1) == itot) then
        ierank = .true.
      else
        ierank = .false.
      end if

      if (zstart(2) == 1) then
        jbrank = .true.
      else
        jbrank = .false.
      end if

      if (zend(2) == jtot) then
        jerank = .true.
      else
        jerank = .false.
      end if

      !write(*,*) "myid, ibrank, ierank", myid, ibrank, ierank

      ! Global constants

      ! Select advection scheme for scalars. If not set in the options file, the momentum scheme is used
      if (iadv_tke < 0) iadv_tke = iadv_mom
      if (iadv_thl < 0) iadv_thl = iadv_mom
      if (iadv_qt < 0) iadv_qt = iadv_mom

      !CvH remove where
      !where (iadv_sv<0)  iadv_sv  = iadv_mom

      !tg3315 added - only uses kappa advection scheme...
      do n = 1, nsv
         iadv_sv(n) = iadv_kappa
      end do
      !ends here

      phi = xlat*pi/180.
      colat = cos(phi)
      silat = sin(phi)

      omega = 7.292e-5
      omega_gs = 7.292e-5
      om22 = 2.*omega*colat
      om23 = 2.*omega*silat
      om22_gs = 2.*omega_gs*colat
      om23_gs = 2.*omega_gs*silat

      ! Variables
      allocate (dsv(nsv))

      ! Create the physical grid variables
      allocate (dzf(kb - kh:ke + kh))
      allocate (dzf2(kb - kh:ke + kh))
      allocate (dzfi(kb - kh:ke + kh))
      allocate (dzfiq(kb - kh:ke + kh))
      allocate (dzfi5(kb - kh:ke + kh))
      allocate (dzh(kb:ke + kh))
      allocate (dzhi(kb:ke + kh))
      allocate (dzhiq(kb:ke + kh))
      allocate (dzh2i(kb:ke + kh))
      allocate (zh(kb:ke + kh))
      allocate (zf(kb:ke + kh))

      allocate (dxf(ib-ih:itot+ih))
      allocate (dxf2(ib-ih:itot+ih))
      allocate (dxfi(ib-ih:itot+ih))
      allocate (dxfiq(ib-ih:itot+ih))
      allocate (dxfi5(ib-ih:itot+ih))
      allocate (dxh(ib:itot+ih))
      allocate (dxhi(ib:itot+ih))
      allocate (dxhiq(ib:itot+ih))
      allocate (dxh2i(ib:itot+ih))
      allocate (xh(ib:itot+ih))
      allocate (xf(ib:itot+ih))
      allocate (yh(jb:jtot+jh))
      allocate (yf(jb:jtot+jh))
      allocate (delta(ib-ih:itot+ih, kb:ke + kh))

      rslabs = real(itot*jtot)

      dx = xlen/float(itot)
      dy = ylen/float(jtot)

      ! MPI

      ! Note, that the loop for reading zf and calculating zh
      ! has been split so that reading is only done on PE 1

      write (cexpnr, '(i3.3)') iexpnr
      if (nrank == 0) then
         open (ifinput, file='prof.inp.'//cexpnr)
         read (ifinput, '(a72)') chmess
         read (ifinput, '(a72)') chmess

         do k = kb, ke
            read (ifinput, *) zf(k)
         end do
         close (ifinput)

         ! ! J. Tomas: Read the x-coordinates of the cell centers from xgrid.inp.XXX
         ! ! SO: still reads for now, but need to remove any reference to xf, xh, etc eventually
         ! open (ifinput, file='xgrid.inp.'//cexpnr)
         ! read (ifinput, '(a72)') chmess
         ! read (ifinput, '(a72)') chmess
         !
         ! do i = ib, itot
         !    read (ifinput, *) xf(i)
         ! end do
         ! close (ifinput)

      end if ! end if nrank==0

      ! MPI broadcast ktot elements from zf
      call MPI_BCAST(zf, ktot, MY_REAL, 0, comm3d, mpierr)
      ! MPI broadcast itot elements from xf
      ! call MPI_BCAST(xf, itot, MY_REAL, 0, comm3d, mpierr)

      zh(kb) = 0.0
      do k = kb, ke
         zh(k + 1) = zh(k) + 2.0*(zf(k) - zh(k))
      end do
      zf(ke + kh) = zf(ke) + 2.0*(zh(ke + kh) - zf(ke))

      do k = kb, ke
         dzf(k) = zh(k + 1) - zh(k)
      end do
      dzf(ke + 1) = dzf(ke)
      dzf(kb - 1) = dzf(kb)

      dzh(kb) = 2*zf(kb)
      do k = kb + 1, ke + kh
         dzh(k) = zf(k) - zf(k - 1)
      end do

      ! j. tomas: same trick for x-direction...
      ! xh(ib) = 0.0
      ! do i = ib, itot
      !    xh(i + 1) = xh(i) + 2.0*(xf(i) - xh(i))
      ! end do
      ! xf(itot + ih) = xf(itot) + 2.0*(xh(itot + ih) - xf(itot))

      do i=ib,itot+ih
        xh(i) = (i-1) * dx
        xf(i) = xh(i) + dx/2
      end do

      do j=jb,jtot+jh
        yh(j) = (j-1) * dy
        yf(j) = yh(j) + dy/2
      end do

      ! These should be removed eventually
      do i = ib, itot
         dxf(i) = xh(i + 1) - xh(i)
      end do
      dxf(itot + 1) = dxf(itot)
      dxf(ib - 1) = dxf(ib)

      dxh(ib) = 2*xf(ib)
      do i = 2, itot + ih
         dxh(i) = xf(i) - xf(i - 1)
      end do

      do k = kb, ke + kh
         do i = ib - ih, itot + ih
            delta(i, k) = (dxf(i)*dy*dzf(k))**(1./3.)
         end do
      end do

      !--------------------------------------------------
      ! *** Check whether the grid is equidistant *****
      !--------------------------------------------------

      !if (myid == 0) then
      !do k=kb,ke+kh
      !if (.not.(dzf(k).eq.dzf(1)))
      !      write (6, *) &
      !      'WARNING, You are working with a non-equidistant grid!!!!'
      !end if
      !end do
      !end if ! end if myid==0

      dzhi = 1./dzh
      dzfi = 1./dzf
      dzf2 = dzf*dzf
      dxhi = 1./dxh
      dxfi = 1./dxf
      dxf2 = dxf*dxf
      dxi = 1./dx
      dx2 = dx*dx
      dyi = 1./dy
      dy2 = dy*dy

      dzhiq = 0.25*dzhi
      dzfiq = 0.25*dzfi
      dxhiq = 0.25*dxhi
      dxfiq = 0.25*dxfi
      dyiq = 0.25*dyi
      dxiq = 0.25*dxi

      dzh2i = dzhi*dzhi
      dxh2i = dxhi*dxhi
      dy2i = dyi*dyi
      dx2i = dxi*dxi

      dzfi5 = 0.5*dzfi
      dxfi5 = 0.5*dxfi
      dyi5 = 0.5*dyi
      dxi5 = 0.5*dxi

      ! Grid used in kappa scheme advection (extra ghost nodes)
      if (any(iadv_sv(1:nsv) == iadv_kappa) .or. (iadv_thl == iadv_kappa)) then
         allocate (dzfc(kb - khc:ke + khc))
         allocate (dxfc(ib - ihc:itot + ihc))
         allocate (dzfci(kb - khc:ke + khc))
         allocate (dxfci(ib - ihc:itot + ihc))
         allocate (dzhci(kb - 1:ke + khc))
         allocate (dxhci(ib - 1:itot + ihc))

         dzfc(kb - kh:ke + kh) = dzf(kb - kh:ke + kh)
         dzfc(kb - khc) = dzfc(kb - kh)
         dzfc(ke + khc) = dzfc(ke + kh)

         dxfc(ib - ih:itot + ih) = dxf(ib - ih:itot + ih)
         dxfc(ib - ihc) = dxfc(ib - ih)
         dxfc(itot + ihc) = dxfc(itot + ih)

         dzhci(kb:ke + kh) = dzhi(kb:ke + kh)
         dzhci(kb - 1) = dzhci(kb)
         dzhci(ke + khc) = dzhci(ke + kh)

         dxhci(ib:itot + ih) = dxhi(ib:itot + ih)
         dxhci(ib - 1) = dxhci(ib)
         dxhci(itot + ihc) = dxhci(itot + ih)

         dzfci = 1./dzfc
         dxfci = 1./dxfc
      end if

      tnextrestart = trestart
      tnextfielddump = tfielddump
!    tnextstatsdump = tstatsdump
      timeleft = runtime ! tg3315 previously btime + runtime

   end subroutine initglobal