op_obj_Init Subroutine

    impure subroutine op_obj_Init(this,block,parallel,order)
      !> Initializes object.
      ! The operator order can be specified by the optional parameter.
      ! Otherwise, the object initializes 2nd order operators by default.
      implicit none
      class(op_obj),      intent(inout) :: this                                !! Differential operators utility
      type(block_obj),    intent(in),    &
                                 target :: block                               !! A block object
      type(parallel_obj), intent(in),    &
                                 target :: parallel                            !! Parallel structure to link with
      integer,            intent(in),    &
                               optional :: Order                               !! Order of interpolation/differentiation schemes
      ! Work variables
      integer :: i,j,k,n

      ! Point to the master objects
      this%parallel => parallel
      this%block    => block

      associate (                                 &
        lo => this%block%lo, hi => this%block%hi, &
        xm => this%block%xm, x  => this%block%x,  &
        ym => this%block%ym, y  => this%block%y,  &
        zm => this%block%zm, z  => this%block%z,  &
        st=>this%st, scheme_order=>this%scheme_order)

        ! Set scheme order, if explicitely stated,
        ! or use default 2nd order
        if (present(order)) scheme_order = order

        ! Set stencil size
        st = scheme_order/2
        if (this%block%ngc.lt.2*st-1) call this%parallel%Stop("Insufficient number of ghost cells for the requested scheme")

        ! Basic interpolation and differentiation
        allocate(this%morinishi_int(2*st, st)); this%morinishi_int = 0.0_wp
        allocate(this%morinishi_ddx(2*st, st)); this%morinishi_ddx = 0.0_wp
        allocate(this%alpha(st))

        ! Eq. 139 and 140
        ! Fourth order: Eq. 32 and 35 in paper
        ! Second order: Eq. 87
        do n=1,st
          this%morinishi_int(st - (n-1),n) =  0.5_wp
          this%morinishi_int(st +  n   ,n) =  0.5_wp
          this%morinishi_ddx(st - (n-1),n) = -1.0_wp/real(2*n-1,wp)
          this%morinishi_ddx(st +  n   ,n) =  1.0_wp/real(2*n-1,wp)
        end do

        ! Solutions to Eq. 12 in Desjardins et al, 2008.
        select case (scheme_order)
        case (2)
          this%alpha=[1.0_wp]
        case (4)
          this%alpha=[9.0_wp/8.0_wp, -1.0_wp/8.0_wp]
        case (6)
          this%alpha=[75.0_wp/64.0_wp, -25.0_wp/128.0_wp, 3.0_wp/128.0_wp]
        case default
          error stop "Only second, fourth, and sixth order schemes for now"
        end select

        ! Coefficients for differentiation and interpolation
        ! -----------------------------------------
        allocate(this%c_d1dx1  (-st+1:st,lo(1)-st:hi(1)+st-1)); this%c_d1dx1   = 0.0_wp
        allocate(this%c_d1dx2  (-st+1:st,lo(2)-st:hi(2)+st-1)); this%c_d1dx2   = 0.0_wp
        allocate(this%c_d1dx3  (-st+1:st,lo(3)-st:hi(3)+st-1)); this%c_d1dx3   = 0.0_wp
        allocate(this%c_intrp1 (-st+1:st,lo(1)-st:hi(1)+st-1)); this%c_intrp1  = 0.0_wp
        allocate(this%c_intrp2 (-st+1:st,lo(2)-st:hi(2)+st-1)); this%c_intrp2  = 0.0_wp
        allocate(this%c_intrp3 (-st+1:st,lo(3)-st:hi(3)+st-1)); this%c_intrp3  = 0.0_wp

        allocate(this%c_d1dx1m (-st:st-1,lo(1)-st+1:hi(1)+st)); this%c_d1dx1m  = 0.0_wp
        allocate(this%c_d1dx2m (-st:st-1,lo(2)-st+1:hi(2)+st)); this%c_d1dx2m  = 0.0_wp
        allocate(this%c_d1dx3m (-st:st-1,lo(3)-st+1:hi(3)+st)); this%c_d1dx3m  = 0.0_wp
        allocate(this%c_intrp1m(-st:st-1,lo(1)-st+1:hi(1)+st)); this%c_intrp1m = 0.0_wp
        allocate(this%c_intrp2m(-st:st-1,lo(2)-st+1:hi(2)+st)); this%c_intrp2m = 0.0_wp
        allocate(this%c_intrp3m(-st:st-1,lo(3)-st+1:hi(3)+st)); this%c_intrp3m = 0.0_wp

        ! 1-dir
        do i=lo(1)-st,hi(1)+st-1
          do n=1,st
            this%c_d1dx1  (-st+1:st,i) = this%c_d1dx1 (-st+1:st,i)  + this%alpha(n)*this%morinishi_ddx(:,n)/(x (i+1)-x (i  ))
            this%c_intrp1 (-st+1:st,i) = this%c_intrp1 (-st+1:st,i) + this%alpha(n)*this%morinishi_int(:,n)
          end do
        end do
        do i=lo(1)-st+1,hi(1)+st
          do n=1,st
            this%c_d1dx1m (-st:st-1,i) = this%c_d1dx1m(-st:st-1,i)  + this%alpha(n)*this%morinishi_ddx(:,n)/(xm(i  )-xm(i-1))
            this%c_intrp1m(-st:st-1,i) = this%c_intrp1m(-st:st-1,i) + this%alpha(n)*this%morinishi_int(:,n)
          end do
        end do
        ! 2-dir
        do j=lo(2)-st,hi(2)+st-1
          do n=1,st
            this%c_d1dx2  (-st+1:st,j) = this%c_d1dx2 (-st+1:st,j)  + this%alpha(n)*this%morinishi_ddx(:,n)/(y (j+1)-y (j  ))
            this%c_intrp2 (-st+1:st,j) = this%c_intrp2 (-st+1:st,j) + this%alpha(n)*this%morinishi_int(:,n)
          end do
        end do
        do j=lo(2)-st+1,hi(2)+st
          do n=1,st
            this%c_d1dx2m (-st:st-1,j) = this%c_d1dx2m(-st:st-1,j)  + this%alpha(n)*this%morinishi_ddx(:,n)/(ym(j  )-ym(j-1))
            this%c_intrp2m(-st:st-1,j) = this%c_intrp2m(-st:st-1,j) + this%alpha(n)*this%morinishi_int(:,n)
          end do
        end do
        ! 3-dir
        do k=lo(3)-st,hi(3)+st-1
          do n=1,st
            this%c_d1dx3  (-st+1:st,k) = this%c_d1dx3 (-st+1:st,k)  + this%alpha(n)*this%morinishi_ddx(:,n)/(z (k+1)-z (k  ))
            this%c_intrp3 (-st+1:st,k) = this%c_intrp3 (-st+1:st,k) + this%alpha(n)*this%morinishi_int(:,n)
          end do
        end do
        do k=lo(3)-st+1,hi(3)+st
          do n=1,st
            this%c_d1dx3m (-st:st-1,k) = this%c_d1dx3m(-st:st-1,k)  + this%alpha(n)*this%morinishi_ddx(:,n)/(zm(k  )-zm(k-1))
            this%c_intrp3m(-st:st-1,k) = this%c_intrp3m(-st:st-1,k) + this%alpha(n)*this%morinishi_int(:,n)
          end do
        end do
      end associate

      return
    end subroutine op_obj_Init