anchor.diff:      end do
       end do
   
  -!    if ( iret == 1 ) then
  -      call cpm_Isend(cf_x(1,1), cf_sz(1), CPM_REAL, nID(X_MINUS), 0, key(X_MINUS,1), ierr)
  -    endif
  -!  endif
  +    call cpm_Isend(cf_x(1,1), cf_sz(1), CPM_REAL, nID(X_MINUS), 0, key(X_MINUS,1), ierr)
  +  endif
   
   ! recv
     if( nID(X_PLUS).ge.0 ) then
  -!    if ( iret == 1 ) then
  -      call cpm_Irecv(cf_x(1,3), cf_sz(1), CPM_REAL, nID(X_PLUS), 0, key(X_MINUS,2), ierr)
  -!    end if
  +    call cpm_Irecv(cf_x(1,3), cf_sz(1), CPM_REAL, nID(X_PLUS), 0, key(X_MINUS,2), ierr)
     endif
   
   ! X_PLUS
  
positive.code:     end do
      end do
  
  !    if ( iret == 1 ) then
        call cpm_Isend(cf_x(1,2), cf_sz(1), CPM_REAL, nID(X_PLUS), 0, key(X_PLUS,1), ierr)
  !    end if
    endif
  
  ! recv
    if( nID(X_MINUS).ge.0 ) then
  !    if ( iret == 1 ) then
        call cpm_Irecv(cf_x(1,4), cf_sz(1), CPM_REAL, nID(X_MINUS), 0, key(X_PLUS,2), ierr)
  !    endif
    endif
  
  ! Y_MINUS
  
negative.code: return
  end subroutine set_neum8
  
  
  !> ********************************************************************
  !! @brief Residula-Cut SOR
  !! @param[in,out] p 圧力
  !! @param sz 配列長
  !! @param g ガイドセル長
  !! @param omg 加速係数
  !! @param[out] res 絶対残差と相対残差
  !! @param src0 固定ソース項
  !! @param src1 反復毎に変化するソース項
  !! @param bp BCindex P
  !! @param[out] flop
  !<
    subroutine psor_m1 (para_key, p, sz, g, omg, e, res, wrk_m, bnd, cm_mode, do_time, cm_time, dtype)
    implicit none
    include 'ffv_f_params.h'
    integer                                                   ::  para_key, g, dtype
    integer, dimension(3)                                     ::  sz
    real                                                      ::  omg, e, res
    real, dimension(1-g:sz(1)+g, 1-g:sz(2)+g, 1-g:sz(3)+g)    ::  wrk_m, p, bnd
    integer                                                   ::  cm_mode ! 0:CommBndCell, 1:CommBndCell2
    integer                                                   ::  do_time
    double precision                                          ::  cm_time
  
    integer                                                   ::  ierr, wait_num, req(12)
    integer                                                   ::  i, j, k, ix, jx, kx, iret
    real                                                      ::  cpd, s0, ee
    real                                                      ::  f1, f2, f3, f4, f5, f6, f0
    real                                                      ::  tmp, r1, dp
    double precision                                          ::  prev, TuneGetTime
  
    ix = sz(1)
    jx = sz(2)
    kx = sz(3)
    ee= 1.0+e
    r1 = 0.0
  
  !$OMP PARALLEL SHARED  (ix,jx,kx,p,bnd,wrk,omg,e,ee) &
  !              PRIVATE (k,j,i,s0,ss,f0,f1,f2,f3,f4,f5,f6,cpd)
  !$OMP BARRIER
  !$OMP DO
    do k=1,kx
    do j=1,jx
    do i=1,ix
      f0=bnd(i,j,k)
      f1=bnd(i+1,j  ,k  )
      f2=bnd(i-1,j  ,k  )
      f3=bnd(i  ,j+1,k  )
      f4=bnd(i  ,j-1,k  )
      f5=bnd(i  ,j  ,k+1)
      f6=bnd(i  ,j  ,k-1)
      cpd=ee/(e+(f1+f2+f3+f4+f5+f6)*f0)
  
      s0=  (f1*p(i+1,j  ,k  )+f2*p(i-1,j  ,k  )  &
          + f3*p(i  ,j+1,k  )+f4*p(i  ,j-1,k  )  &
          + f5*p(i  ,j  ,k+1)+f6*p(i  ,j  ,k-1))*f0
      dp = wrk_m(i,j,k) + (cpd*s0-p(i,j,k))*f0
      p(i,j,k)=p(i,j,k) + omg*dp
      r1 = r1 + dp*dp
    end do
    end do
    end do
  !$OMP END DO NOWAIT
  !$OMP END PARALLEL
  
    res = r1
    call SklIsParallel(iret)
    if ( iret .eq. 1 ) then
      tmp = r1
      if ( dtype == 4 ) then
        call SklAllreduce(tmp, res, 1, SKL_REAL4, SKL_SUM, SKL_DEFAULT_GROUP, ierr)
      else
        call SklAllreduce(tmp, res, 1, SKL_REAL8, SKL_SUM, SKL_DEFAULT_GROUP, ierr)
      endif
    end if
  
    if( do_time.eq.1 ) then
      call SklBarrier(SKL_DEFAULT_GROUP, ierr)
      prev = TuneGetTime()
    endif
  
    if( cm_mode.eq.0 ) then    ! $BF14|DL?.(B
      call SklCommBndCell(p, 1, ierr)
    else                        ! $BHsF14|DL?.(B
      call SklCommBndCell2(p, 1, wait_num, req, ierr)
      call SklWaitAll(wait_num, req, ierr);
    endif
  
    if( do_time.eq.1 ) then
      cm_time = cm_time + (TuneGetTime()-prev)
    endif
  
    return
    end subroutine psor_m1
  
  !   ***************************************************
    subroutine cbs3d_gmressor (para_key, p, sz, g, omg, res, wrk, bnd, commmode, do_time, commtime, dtype, iparam, rparam)
    implicit none
    include 'sklparaf.h'
  
    integer                                                   ::  para_key, g, dtype
    integer                                                   ::  iparam(*)
    integer                                                   ::  commmode ! 0:CommBndCell, 1:Commface(NoHide), 2:Commface(Hide)
    integer                                                   ::  do_time, commtime
    integer, dimension(3)                                     ::  sz
    real                                                      ::  omg, res
    real*8                                                    ::  rparam(*)
    real, dimension(1-g:sz(1)+g, 1-g:sz(2)+g, 1-g:sz(3)+g)    ::  wrk, p, bnd
  
    real,dimension(:,:,:),allocatable                         ::  xt, yt, rest, res_fct
    real,dimension(:,:,:,:),allocatable                       ::  vm, zm
    integer                                                   ::  i, j, k, ix, jx, kx, it, m
    integer                                                   ::  ig, jg, kg, lg, im, jm, km, lm
    integer                                                   ::  i_iter, n_iter, tneum
    integer                                                   ::  iter, ierr, nrm, isfin
    integer, parameter                                        ::  m_max = 15, iter_max = 100
    integer, parameter                                        ::  isneum = 0
    integer                                                   ::  i_inner, n_inner, oki, step
    real                                                      ::  e, ee, fc, eps_abs, res_abs
    real*8                                                    ::  rgm(m_max, m_max)
    real*8                                                    ::  hgm(m_max+1, m_max)
    real*8                                                    ::  cgm(m_max), sgm(m_max)
    real*8                                                    ::  bgm(m_max), ygm(m_max)
    real*8                                                    ::  prev, TuneGetTime
    real*8                                                    ::  al, t_eps, f_v, beta, beta_1
    real*8, parameter                                         ::  fct1 = 1.0d-4, fct2 = 6.0d0
    integer, dimension(2)                                     ::  err
    real*4, dimension(5)                                      ::  r4
  
    if(iparam(8).lt.0) return
    err(1) = 0
    isfin = 0
  !    n_iter = min(iparam(8), iter_max)
  !    nrc    = min(iparam(8), nrc_max)
    oki = 4; step = 6
    m       = m_max
    n_iter  = iter_max
  
    if (dtype == 4) then
      e = 2.4e-7
      t_eps = 0.995*rparam(1)*rparam(1)
    else
      e = 4.4e-16
      t_eps = 0.999*rparam(1)*rparam(1)
    endif
    ee = 1.0 + e
  
    ix = sz(1)
    jx = sz(2)
    kx = sz(3)
    lg = 1 - g
    ig = ix + g
    jg = jx + g
    kg = kx + g
    tneum = 0
    if((isneum.ne.0).and.(g.gt.0)) then
      tneum = 1
    endif
  
  ! Distributed like "wrk", "p", "bnd"
    allocate(xt(lg:ig,lg:jg,lg:kg),      yt(lg:ig,lg:jg,lg:kg), rest(lg:ig,lg:jg,lg:kg), res_fct(lg:ig,lg:jg,lg:kg))
    allocate(vm(ix,jx,kx,m_max), zm(lg:ig,lg:jg,lg:kg,m_max))
  
  !$OMP PARALLEL SHARED  (ix,jx,kx,p,bnd,wrk,omg,e) &
  !              PRIVATE (pl,pc,pr,ytmp,k,j,i,it,fc)
    f_v = 0.0
  !$OMP BARRIER
  !$OMP DO
    do k=1,kx
    do j=1,jx
    do i=1,ix
      fc = e + (bnd(i+1,j  ,k  ) + bnd(i-1,j  ,k  ) &
              + bnd(i  ,j+1,k  ) + bnd(i  ,j-1,k  ) &
              + bnd(i  ,j  ,k+1) + bnd(i  ,j  ,k-1))*bnd(i,j,k)
      res_fct(i,j,k) = bnd(i,j,k)*ee/fc
      al = res_fct(i,j,k)*wrk(i,j,k)
      f_v = f_v + al*al
    end do
    end do
    end do
  !$OMP END DO
  !$OMP BARRIER
  
    if (f_v.lt.1.0e-30) return
    t_eps   = f_v*t_eps
  
    call cbs3d_psor (para_key, p, sz, g, omg, res, wrk, bnd, commmode, do_time, commtime, dtype)
    r4(1) = sqrt(res/f_v)
  
    if(res.lt.t_eps) then
      print*,'Final',r4(1)
      isfin = 1
      goto 100
    endif
  
    if(tneum.eq.1) then
      call cbs3d_set_neum (p, sz, g)
    endif
  
  !    if((res.gt.(fct1*f_v)).and.(res.gt.(fct2*t_eps))) then
  
    if(res.gt.(fct2*t_eps)) then
      t_eps = res/fct2
  !    else
  !      t_eps = sqrt(res*t_eps)
    endif
  
    eps_abs = 0.15*t_eps
    iter = 0
    nrm = 0
  
    do i_iter=1,n_iter
      if(tneum.eq.1) then
        call cbs3d_psor_neumann (para_key, p, sz, g, omg, res, wrk, bnd, commmode, do_time, commtime, dtype)
        call cbs3d_set_neum (p, sz, g)
      else
        call cbs3d_psor (para_key, p, sz, g, omg, res, wrk, bnd, commmode, do_time, commtime, dtype)
      endif
  
      iparam(2) = iparam(2) + 1
  
      if(res.lt.t_eps) goto 70
  !      if(res.lt.t_eps.and.i_iter.gt.1) goto 70
  !    print*,'Intermid',i_iter,res/f_v
  
      do k=1,kx
      do j=1,jx
      do i=1,ix
        rest(i,j,k) = res_fct(i,j,k)*wrk(i,j,k)
      end do
      end do
      end do
  
      call cbs3d_mvprod (p, sz, g, e, bnd, yt, res_fct)
      iparam(3) = iparam(3) + 1
  
      call cbs3d_ressmpl (sz, g, res_abs, rest, yt)
      beta = sqrt(res_abs)
  
      if(beta.lt.1.0e-30) goto 60
      beta_1 = 1.0/beta
  
  !$OMP DO
      do k=1,kx
      do j=1,jx
      do i=1,ix
        vm(i,j,k,1) = beta_1*rest(i,j,k)
      end do
      end do
      end do
  
      bgm(1) = beta
  
      do lm=2,m
        bgm(lm) = 0.0;
      enddo
  
      if(m.gt.m_max) m = m_max
  
  !$OMP END DO
      do im=1,m
        iter = iter + 1
  !$OMP BARRIER
  
        do k=1,kx
        do j=1,jx
        do i=1,ix
          rest(i,j,k) = vm(i,j,k,im)
        end do
        end do
        end do
  
  !$OMP DO
        do k=lg,kg
        do j=lg,jg
        do i=lg,ig
          xt(i,j,k) = 0.0
        end do
        end do
        end do
  !$OMP END DO
  
        if(    iter.le.  oki) then; n_inner =   step
        elseif(iter.le.2*oki) then; n_inner = 2*step
        elseif(iter.le.3*oki) then; n_inner = 3*step
        elseif(iter.le.4*oki) then; n_inner = 4*step
        else;                       n_inner = 5*step
        endif
  
        if(n_inner.lt.1) n_inner = 1
  
        do i_inner = 1, n_inner
          if(tneum.eq.1) then
            call cbs3d_psor_neumann_m1(para_key, xt, sz, g, omg, e, res, rest, bnd, commmode, do_time, commtime, dtype)
            call cbs3d_set_neum (xt, sz, g)
          else
            call cbs3d_psor_m1 (para_key, xt, sz, g, omg, e, res, rest, bnd, commmode, do_time, commtime, dtype)
          endif
  
          iparam(2) = iparam(2) + 1
        enddo
  
  
        do k=lg,kg
        do j=lg,jg
        do i=lg,ig
          zm(i,j,k,im) = xt(i,j,k)
        end do
        end do
        end do
  
        call cbs3d_mvprod (xt, sz, g, e, bnd, yt, res_fct)
        iparam(3) = iparam(3) + 1
  
        do km=1,im
          al = 0.0
  
          do k=1,kx
          do j=1,jx
          do i=1,ix
            al = al + vm(i,j,k,km)*yt(i,j,k)
          end do
          end do
          end do
  
          hgm(km,im) = al
  
          do k=1,kx
          do j=1,jx
          do i=1,ix
            yt(i,j,k) = yt(i,j,k) - al*vm(i,j,k,km)
          end do
          end do
          end do
        end do
  
        al = 0.0
    
        do k=1,kx
        do j=1,jx
        do i=1,ix
          al = al + yt(i,j,k)*yt(i,j,k)
        end do
        end do
        end do
  
        hgm(im+1,im) = sqrt(al)
  
        if(hgm(im+1,im).lt.1.0e-30) then
          nrm = im-1
          goto 50
        endif
  
        if(im.lt.m) then
          al = 1.0/hgm(im+1,im)
  
          do k=1,kx
          do j=1,jx
          do i=1,ix
            vm(i,j,k,im+1) = al*yt(i,j,k)
          end do
          end do
          end do
        endif
  
        rgm(1,im) = hgm(1,im)
  
        do km=2,im
          rgm(km  ,im) = cgm(km-1)*hgm(km,im) - sgm(km-1)*rgm(km-1,im)
          rgm(km-1,im) = sgm(km-1)*hgm(km,im) + cgm(km-1)*rgm(km-1,im)
        end do
  
        al = sqrt(rgm(im,im)*rgm(im,im) + hgm(im+1,im)*hgm(im+1,im))
  
        if(al.lt.1.0e-30) then
          nrm = im-1
          goto 50
        endif
  
        cgm(im) = rgm(im,im)/al
        sgm(im) = hgm(im+1,im)/al
        rgm(im,im) = cgm(im)*rgm(im,im) + sgm(im)*hgm(im+1,im)
        bgm(im+1) = - sgm(im)*bgm(im)
        bgm(im) = cgm(im)*bgm(im)
        al = bgm(im+1)*bgm(im+1)
        iparam(1) = iparam(1) + 1
  
        if(al.lt.eps_abs) then
          nrm = im
          goto 50
        endif
  !      end "im"
      end do
  
      nrm = m
  
  50  ygm(nrm) = bgm(nrm)/rgm(nrm,nrm)
  
      do im = nrm-1, 1, -1
        al = bgm(im)
  
        do jm = im+1, nrm
          al = al - rgm(im,jm)*ygm(jm)
        enddo
  
        ygm(im) = al/rgm(im,im)
      end do
  
      do im = 1, nrm
        al = ygm(im)
  
        do k=lg,kg
        do j=lg,jg
        do i=lg,ig
          p(i,j,k) = p(i,j,k) + al*zm(i,j,k,im)
        end do
        end do
        end do
      end do
  !      end "i_iter"
    end do
  
  60  continue
  
    if(tneum.eq.1) then
      call cbs3d_psor_neumann (para_key, p, sz, g, omg, res, wrk, bnd, commmode, do_time, commtime, dtype)
      call cbs3d_set_neum (p, sz, g)
    else
      call cbs3d_psor (para_key, p, sz, g, omg, res, wrk, bnd, commmode, do_time, commtime, dtype)
    endif
  
  70 continue
  !    call cbs3d_mvprod (p, sz, g, e, bnd, yt, res_fct)
  !    do k=1,kx; do j=1,jx; do i=1,ix
  !      rest(i,j,k) = res_fct(i,j,k)*wrk(i,j,k)
  !    end do; end do; end do
  !    call cbs3d_ressmpl (sz, g, res_abs, rest, yt)
  !    r4(2) = sqrt(res/f_v); r4(3) = sqrt(res_abs/f_v)
  !    print*,'Intermid',i_iter-1,nrm,r4(1),r4(2),r4(3)
  !    r4(2) = sqrt(res/f_v)
  !    print*,'Intermid',i_iter-1,nrm,r4(1),r4(2)
  !$OMP END PARALLEL
  100 continue
  
    if((isfin.eq.0).and.(res.lt.(f_v*rparam(1)*rparam(1)))) res = 100.0*f_v*rparam(1)*rparam(1)
  
    if(allocated(xt))	deallocate(xt)
    if(allocated(yt))	deallocate(yt)
    if(allocated(rest))	deallocate(rest)
    if(allocated(res_fct))	deallocate(res_fct)
    if(allocated(vm))	deallocate(vm)
    if(allocated(zm))	deallocate(zm)
  
    return
    end subroutine cbs3d_gmressor
  
  
  !   ***************************************************
  subroutine cbs3d_gmresadi (para_key, p, sz, g, res, wrk, bnd,         &
  commmode, do_time, commtime, dtype, iparam, rparam)
  implicit none
  include 'sklparaf.h'
  
  integer                                                   ::  para_key, g, dtype
  integer                                                   ::  iparam(*)
  integer                                                   ::  commmode ! 0:CommBndCell, 1:Commface(NoHide), 2:Commface(Hide)
  integer                                                   ::  do_time, commtime
  integer, dimension(3)                                     ::  sz
  real                                                      ::  res
  real*8                                                    ::  rparam(*)
  real, dimension(1-g:sz(1)+g, 1-g:sz(2)+g, 1-g:sz(3)+g)    ::  wrk, p, bnd
  
  real,dimension(:),allocatable                             ::  ytmp, pl, pc, pr
  real,dimension(:,:,:),allocatable                         ::  xt, yt, rest, res_fct
  real,dimension(:,:,:,:),allocatable                       ::  vm, zm
  integer                                                   ::  i, j, k, ix, jx, kx, it, m
  integer                                                   ::  ig, jg, kg, lg, im, jm, km, lm
  integer                                                   ::  i_iter, n_iter, tneum
  integer                                                   ::  iter, ierr, nrm, dir, la
  integer, parameter                                        ::  m_max = 15, iter_max = 100
  integer, parameter                                        ::  isneum = 0, is_fast = 1
  integer                                                   ::  i_inner, n_inner, oki, step
  real                                                      ::  e, ee, fc, tomg
  real                                                      ::  eps_abs, res_abs
  real*8                                                    ::  rgm(m_max, m_max)
  real*8                                                    ::  hgm(m_max+1, m_max)
  real*8                                                    ::  cgm(m_max), sgm(m_max)
  real*8                                                    ::  bgm(m_max), ygm(m_max)
  real*8                                                    ::  prev, TuneGetTime
  real*8                                                    ::  al, t_eps, f_v, t_tiny
  real*8                                                    ::  beta, beta_1
  
  integer, dimension(2)                                     ::  err
  
  if(iparam(8).lt.0) return
  err(1) = 0
  iparam(1) = 0
  iparam(2) = 0
  iparam(3) = 0
  !    n_iter = min(iparam(8), iter_max)
  !    nrc    = min(iparam(8), nrc_max)
  oki = 3; step = 3
  m       = m_max
  n_iter  = iter_max
  if (dtype == 4) then
  e = 2.4e-7
  t_eps = 0.995*rparam(1)*rparam(1)
  else
  e = 4.4e-16
  t_eps = 0.999*rparam(1)*rparam(1)
  endif
  eps_abs = 0.15*rparam(1)*rparam(1)
  tomg    = 1.2
  ee = 1.0 + e
  t_tiny = 1.0e-36
  
  ix = sz(1)
  jx = sz(2)
  kx = sz(3)
  lg = 1 - g; ig = ix + g; jg = jx + g; kg = kx + g
  la = max(ix, jx, kx); tneum = 0
  if((isneum.ne.0).and.(g.gt.0)) then; la = la + 2; tneum = 1; endif
  ! Distributed like "wrk", "p", "bnd"
  allocate(xt(lg:ig,lg:jg,lg:kg),      yt(lg:ig,lg:jg,lg:kg), &
  rest(lg:ig,lg:jg,lg:kg), res_fct(lg:ig,lg:jg,lg:kg))
  allocate(vm(ix,jx,kx,m_max), zm(lg:ig,lg:jg,lg:kg,m_max))
  ! Allocate on each processor
  allocate(ytmp(la), pl(la), pc(la), pr(la))
  
  !$OMP PARALLEL SHARED  (ix,jx,kx,p,bnd,wrk,tomg,e) &
  !              PRIVATE (pl,pc,pr,ytmp,k,j,i,it,fc)
  f_v = 0.0
  !$OMP BARRIER
  !$OMP DO
  do k=1,kx; do j=1,jx; do i=1,ix
  fc = e + (bnd(i+1,j  ,k  ) + bnd(i-1,j  ,k  ) &
  + bnd(i  ,j+1,k  ) + bnd(i  ,j-1,k  ) &
  + bnd(i  ,j  ,k+1) + bnd(i  ,j  ,k-1))*bnd(i,j,k)
  res_fct(i,j,k) = bnd(i,j,k)*ee/fc
  al = res_fct(i,j,k)*wrk(i,j,k)
  f_v = f_v + al*al
  end do; end do; end do
  !$OMP END DO
  !$OMP BARRIER
  if(f_v.lt.1.0e-30) return
  t_eps   = f_v*t_eps
  eps_abs = f_v*eps_abs
  iter = 0
  dir = 1
  nrm = 0
  
  do i_iter=1,n_iter
  call cbs3d_psor (para_key, p, sz, g, tomg, res, wrk, bnd, commmode, do_time, commtime, dtype)
  iparam(2) = iparam(2) + 1
  if(res.lt.t_eps) goto 70
  !    print*,'Intermid',i_iter,res
  do k=1,kx; do j=1,jx; do i=1,ix
  rest(i,j,k) = res_fct(i,j,k)*wrk(i,j,k)
  end do; end do; end do
  call cbs3d_mvprod (p, sz, g, e, bnd, yt, res_fct)
  call cbs3d_ressmpl (sz, g, res_abs, rest, yt)
  iparam(3) = iparam(3) + 1
  beta = sqrt(res_abs)
  if(beta.lt.1.0e-30) goto 70
  beta_1 = 1.0/beta
  !$OMP DO
  do k=1,kx; do j=1,jx; do i=1,ix
  vm(i,j,k,1) = beta_1*rest(i,j,k)
  end do; end do; end do
  !$OMP END DO
  bgm(1) = beta; do lm=2,m; bgm(lm) = 0.0; enddo
  if(m.gt.m_max) m = m_max
  
  do im=1,m
  iter = iter + 1
  !$OMP BARRIER
  !$OMP DO
  do k=1,kx; do j=1,jx; do i=1,ix
  if(abs(res_fct(i,j,k)).gt.e) then
  rest(i,j,k) = vm(i,j,k,im)/res_fct(i,j,k)
  else
  rest(i,j,k) = vm(i,j,k,im)/e
  endif
  end do; end do; end do
  !$OMP END DO
  !$OMP DO
  do k=lg,kg; do j=lg,jg; do i=lg,ig
  xt(i,j,k) = 0.0
  end do; end do; end do
  !$OMP END DO
  
  if(    iter.le.  oki) then; n_inner =   step
  elseif(iter.le.2*oki) then; n_inner = 2*step
  elseif(iter.le.3*oki) then; n_inner = 3*step
  elseif(iter.le.4*oki) then; n_inner = 4*step
  else;                       n_inner = 5*step
  endif
  
  if(n_inner.lt.1) n_inner = 1
  do i_inner = 1, n_inner
  if(is_fast.eq.1) then
  call cbs3d_adiiter_x1 (xt, sz, g, rest, bnd, tneum, t_tiny, e,    &
  ytmp, pl, pc, pr, err,     &
  commmode, do_time, commtime)
  else
  if(dir.eq.1) then
  call cbs3d_adiiter_x (xt, sz, g, rest, bnd, tneum, t_tiny, e,    &
  ytmp, pl, pc, pr, err,     &
  commmode, do_time, commtime)
  if(err(1).ne.0) return
  dir = 2
  elseif(dir.eq.2) then
  call cbs3d_adiiter_y (xt, sz, g, rest, bnd, tneum, t_tiny, e,    &
  ytmp, pl, pc, pr, err,     &
  commmode, do_time, commtime)
  if(err(1).ne.0) return
  dir = 3
  else
  call cbs3d_adiiter_z (xt, sz, g, rest, bnd, tneum, t_tiny, e,    &
  ytmp, pl, pc, pr, err,     &
  commmode, do_time, commtime)
  if(err(1).ne.0) return
  dir = 1
  endif
  endif
  iparam(2) = iparam(2) + 1
  enddo
  if(tneum.eq.1) then
  call cbs3d_set_neum (xt, sz, g)
  endif
  
  do k=lg,kg; do j=lg,jg; do i=lg,ig
  zm(i,j,k,im) = xt(i,j,k)
  end do; end do; end do
  call cbs3d_mvprod (xt, sz, g, e, bnd, yt, res_fct)
  iparam(3) = iparam(3) + 1
  do km=1,im
  al = 0.0
  do k=1,kx; do j=1,jx; do i=1,ix
  al = al + vm(i,j,k,km)*yt(i,j,k)
  end do; end do; end do
  hgm(km,im) = al
  do k=1,kx; do j=1,jx; do i=1,ix
  yt(i,j,k) = yt(i,j,k) - al*vm(i,j,k,km)
  end do; end do; end do
  end do
  al = 0.0
  do k=1,kx; do j=1,jx; do i=1,ix
  al = al + yt(i,j,k)*yt(i,j,k)
  end do; end do; end do
  hgm(im+1,im) = sqrt(al)
  if(hgm(im+1,im).lt.1.0e-30) then; nrm = im-1; goto 50; endif
  if(im.lt.m) then
  al = 1.0/hgm(im+1,im)
  do k=1,kx; do j=1,jx; do i=1,ix
  vm(i,j,k,im+1) = al*yt(i,j,k)
  end do; end do; end do
  endif
  rgm(1,im) = hgm(1,im)
  do km=2,im
  rgm(km  ,im) = cgm(km-1)*hgm(km,im) - sgm(km-1)*rgm(km-1,im)
  rgm(km-1,im) = sgm(km-1)*hgm(km,im) + cgm(km-1)*rgm(km-1,im)
  end do
  al = sqrt(rgm(im,im)*rgm(im,im) + hgm(im+1,im)*hgm(im+1,im))
  if(al.lt.1.0e-30) then; nrm = im-1; goto 50; endif
  cgm(im) = rgm(im,im)/al; sgm(im) = hgm(im+1,im)/al
  rgm(im,im) = cgm(im)*rgm(im,im) + sgm(im)*hgm(im+1,im)
  bgm(im+1) = - sgm(im)*bgm(im); bgm(im) = cgm(im)*bgm(im)
  al = bgm(im+1)*bgm(im+1)
  iparam(1) = iparam(1) + 1
  if(al.lt.eps_abs) then; nrm = im; goto 50; endif
  !      end "im"
  end do
  nrm = m
  50  ygm(nrm) = bgm(nrm)/rgm(nrm,nrm)
  do im = nrm-1, 1, -1
  al = bgm(im)
  do jm = im+1, nrm; al = al - rgm(im,jm)*ygm(jm); enddo
  ygm(im) = al/rgm(im,im)
  end do
  do im = 1, nrm
  al = ygm(im)
  do k=lg,kg; do j=lg,jg; do i=lg,ig
  p(i,j,k) = p(i,j,k) + al*zm(i,j,k,im)
  end do; end do; end do
  end do
  !      end "i_iter"
  end do
  
  60  continue
  call cbs3d_psor (para_key, p, sz, g, tomg, res, wrk, bnd, commmode, do_time, commtime, dtype)
  70  continue
  !    call cbs3d_mvprod (p, sz, g, e, bnd, yt, res_fct)
  !    do k=1,kx; do j=1,jx; do i=1,ix
  !      rest(i,j,k) = res_fct(i,j,k)*wrk(i,j,k)
  !    end do; end do; end do
  !    call cbs3d_ressmpl (sz, g, res_abs, rest, yt)
  !    print*,i_iter-1,nrm,res,res_abs,sngl(sqrt(res_abs/f_v))
  !$OMP END PARALLEL
  100 continue
  
  if(allocated(ytmp))	deallocate(ytmp)
  if(allocated(pl))	deallocate(pl)
  if(allocated(pc))	deallocate(pc)
  if(allocated(pr))	deallocate(pr)
  
  if(allocated(xt))	deallocate(xt)
  if(allocated(yt))	deallocate(yt)
  if(allocated(rest))	deallocate(rest)
  if(allocated(res_fct))	deallocate(res_fct)
  if(allocated(vm))	deallocate(vm)
  if(allocated(zm))	deallocate(zm)
  
  return
  end subroutine cbs3d_gmresadi
  