anchor.diff:  return
   end subroutine set_neum8
   
  -!   *******************************************************************************************
  -subroutine psor_m1 (para_key, p, sz, g, omg, e, res, wrk_m, bnd, cm_mode, do_time, cm_time, dtype)
  -implicit none
  -include 'sklparaf.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
  -real*8                                                    ::  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
  -real*8                                                    ::  prev, TuneGetTime
  -
  -ix = sz(1)
  -jx = sz(2)
  -kx = sz(3)
  -ee= 1.0+e
  -r1 = 0.0
  +
  +!> ********************************************************************
  +!! @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
  +  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)
  +  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
  -call SklAllreduce(tmp, res, 1, SKL_REAL8, SKL_SUM, SKL_DEFAULT_GROUP, ierr)
  +e = 4.4e-16
  +t_eps = 0.999*rparam(1)*rparam(1)
   endif
  -end if
  -
  -if( do_time.eq.1 ) then
  -call SklBarrier(SKL_DEFAULT_GROUP, ierr)
  -prev = TuneGetTime()
  +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
  -
  -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);
  +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( do_time.eq.1 ) then
  -cm_time = cm_time + (TuneGetTime()-prev)
  +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 psor_m1
  -
  +end subroutine cbs3d_gmresadi
  
positive.code:   res_0 = sqrt(res_abs/f_v)
  !    if(res_abs.lt.t_eps) goto 10
  
    do i_iter=1,n_iter
  
  !$OMP BARRIER
  !$OMP DO
      do k=lg,kg
      do j=lg,jg
      do i=lg,ig
        xt(i,j,k) = 0.0
      end do
      end do
  
negative.code:     key(i,2) = -1 !recv
    end do
  
  
  ! X_MINUS
  ! send
  