module solver
    use const
    use var
    use comm
    use io
implicit real*8 (a-h,o-z)

contains

subroutine init_solver
    ! -------------------------
    ! evp elliptic solver data
    ! used for surface pressure
    ! -------------------------
!    open(99,file='prepglo/data/evp'//grd,form='unformatted')
    if(lsolver.eq.1) then
	if(myid.eq.0) write(*,*)'evp solver'
!	read(99) al,ar,ab,at,ac,rinv,rinv1,ie
!	close(99)
	elseif(lsolver.eq.2) then
	if(myid.eq.0) write(*,*)'bi-cgstab solver'
        call rdnc_evp
!	read(99) al,ar,ab,at,ac,cl,cb,cc,cr,ct
!	close(99)
    else
	if(myid.eq.0)write(*,*)'lsolver=',lsolver,' (incorrect solver)'
	call mpi_finalize(ierr)
	stop
    endif
end subroutine



! *****************
! elliptic solver *
! *****************
! ----------------------------------------------------------------------
subroutine p_bicgstab_le(ab,al,ac,ar,at,b,x,cb,cl,cc,cr,ct,r,rh,p,v,s,t,ph,sh)
    parameter(i0j0=i0*j0)
    real*8 x,r,rh,p,v,s,t,ph,sh
    real*8 rho,rhn,alpha,beta,w,res,tmp,tp,tol
    dimension ab(i2,j2),al(i2,j2),ac(i2,j2),ar(i2,j2),at(i2,0:j2), &
    b(i2,j2),x(i0,j0),r(i0,j0),rh(i0,j0),p(i0,j0),v(i0,j0),  &
    s(i0,j0),t(i0,j0),ph(i0,j0),sh(i0,j0), &
    cb(i2,j2),cl(i2,j2),cc(i2,j2),cr(i2,j2),ct(i2,j2)

    iter=0
    10   rho=1.d0
    alpha=1.d0
    w=1.d0
    call p_sqprod(ab,al,ac,ar,at,x,r)

    do 50 j=2,j1
    do 50 i=2,i1
    50   r(i,j)=b(i-1,j-1)-r(i,j)

    call dcopy(i0j0,r,1,rh,1)
    call p_inprod(r,r,tp)

    res=sqrt(tp)
    if (res .lt. 1e-7) then
	call mpi_sendrecv(x(2,1),1,m_v8lon,m_w,1,x(i0,1),1,m_v8lon,m_e,1,m_cart,istat,ierr)
	call mpi_sendrecv(x(i1,1),1,m_v8lon,m_e,1,x(1,1),1,m_v8lon,m_w,1,m_cart,istat,ierr)
	call mpi_sendrecv(x(1,2),1,m_v8lat,m_s,1,x(1,j0),1,m_v8lat,m_n,1,m_cart,istat,ierr)
	call mpi_sendrecv(x(1,j1),1,m_v8lat,m_n,1,x(1,1),1,m_v8lat,m_s,1,m_cart,istat,ierr)
	return
    endif

    100  iter=iter+1
    call p_inprod(r,rh,rhn)
    beta=(rhn/rho)*(alpha/w)

    ! print*,iter,'beta=',beta

    call daxpy(i0j0,-1.d0*w,v,1,p,1)
    call dscal(i0j0,beta,p,1)
    call daxpy(i0j0,1.d0,r,1,p,1)

    ! do 200 j=2,j1
    ! do 200 i=2,i1
    ! 200 ph(i,j)=p(i,j)

    ! call p_sqelim(cb,cl,cc,cr,ct,p,ph)
    ! call p_lusqelim(cb,cl,cc,cr,ct,p,ph)
    call p_sipsqelim(cb,cl,cc,cr,ct,p,ph)
    call p_sqprod(ab,al,ac,ar,at,ph,v)
    call p_inprod(rh,v,tmp)
    alpha=rhn/tmp

    ! print*,iter,'alpha=', alpha

    call dcopy(i0j0,r,1,s,1)
    call daxpy(i0j0,-1.d0*alpha,v,1,s,1)
    ! call dcopy(i0j0,s,1,sh,1)

    ! do 300 j=2,j1
    ! do 300 i=2,i1
    ! 300 sh(i,j)=s(i,j)

    call p_sipsqelim(cb,cl,cc,cr,ct,s,sh)
    ! call p_lusqelim(cb,cl,cc,cr,ct,s,sh)
    ! call p_sqelim(cb,cl,cc,cr,ct,s,sh)
    call p_sqprod(ab,al,ac,ar,at,sh,t)
    call p_inprod(t,t,tmp)
    call p_inprod(t,s,w)
    w=w/tmp

    ! print*,iter,'wn=',w

    tmp=0.d0

    call daxpy(i0j0,alpha,ph,1,x,1)
    call daxpy(i0j0,w,sh,1,x,1)
    call dcopy(i0j0,s,1,r,1)
    call daxpy(i0j0,-1.d0*w,t,1,r,1)

    call p_inprod(r,r,tp)

    rho=rhn
    ! if (myid .eq. 0) print*,iter,sqrt(tp)
    ! if (sqrt(tp)/res .gt. 1e-3 .and. mod(iter,100) .ne. 0) goto 100
    ! if (myid .eq. 0) print*, iter,sqrt(tp)
    !if (sqrt(tp) .gt. 1e-7 .and. iter .lt. 10000) goto 10
    if (sqrt(tp).gt.1.e-9.or.sqrt(tp)/res .gt. 1e-3.and. mod(iter,1000) .ne. 0) goto 100
    call mpi_barrier(m_cart,ierr)
    call mpi_sendrecv(x(2,1),1,m_v8lon,m_w,1,x(i0,1),1,m_v8lon,m_e,1,m_cart,istat,ierr)
    call mpi_sendrecv(x(i1,1),1,m_v8lon,m_e,1,x(1,1),1,m_v8lon,m_w,1,m_cart,istat,ierr)
    call mpi_sendrecv(x(1,2),1,m_v8lat,m_s,1,x(1,j0),1,m_v8lat,m_n,1,m_cart,istat,ierr)
    call mpi_sendrecv(x(1,j1),1,m_v8lat,m_n,1,x(1,1),1,m_v8lat,m_s,1,m_cart,istat,ierr)

    continue
end subroutine


! ----------------------------------------------------------------------
subroutine p_sqprod(ab,al,ac,ar,at,x,b)
    ! real*8 ab,al,ac,ar,at,x,b
    real*8 x,b
    dimension ab(i2,j2),al(i2,j2),ac(i2,j2),ar(i2,j2),at(i2,0:j2),x(i0,j0),b(i0,j0)

    call mpi_barrier(m_cart,ierr)
    call mpi_sendrecv(x(2,1),1,m_v8lon,m_w,1,x(i0,1),1,m_v8lon,m_e,1,m_cart,istat,ierr)
    call mpi_sendrecv(x(i1,1),1,m_v8lon,m_e,1,x(1,1),1,m_v8lon,m_w,1,m_cart,istat,ierr)
    call mpi_sendrecv(x(1,2),1,m_v8lat,m_s,1,x(1,j0),1,m_v8lat,m_n,1,m_cart,istat,ierr)
    call mpi_sendrecv(x(1,j1),1,m_v8lat,m_n,1,x(1,1),1,m_v8lat,m_s,1,m_cart,istat,ierr)

    do 100 j=2,j1
    do 100 i=2,i1
    it=i-1
    jt=j-1
    100  b(i,j)=ab(it,jt)*x(i,jt)+al(it,jt)*x(it,j)+ac(it,jt)*x(i,j)+ar(it,jt)*x(i+1,j)+at(it,jt)*x(i,j+1)
end subroutine


! ----------------------------------------------------------------------
subroutine p_inprod(a1,a2,pd)
    real*8 a1,a2,tmp,pd
    dimension a1(i0,j0),a2(i0,j0)

    tmp=0.d0
    do 100 j=2,j1
    do 100 i=2,i1
    100  tmp=tmp+a1(i,j)*a2(i,j)
    call mpi_barrier(m_cart,ierr)
    call mpi_allreduce(tmp,pd,1,mpi_real8,mpi_sum,m_cart,ierr)

    return
end subroutine


! --------------------------
subroutine p_sqelim(cb,cl,cc,cr,ct,b,bh)
    parameter(i3=i2-1,j3=j2-1,i0j0=i0*j0)
    ! real*8 cb,cl,cc,cr,ct,
    real*8 b,bh
    dimension cb(i2,j2),cl(i2,j2),cc(i2,j2),cr(i2,j2),ct(i2,j2),b(i0,j0),bh(i0,j0)

    call dcopy(i0j0,b,1,bh,1)
    do 100 j=2,j2
    jt=j-1
    do 150 i=2,i2
    it=i-1
    bh(i,j)=bh(i,j)*cc(it,jt)
    bh(i+1,j)=bh(i+1,j)-cr(it,jt)*bh(i,j)
    150  bh(i,j+1)=bh(i,j+1)-ct(it,jt)*bh(i,j)
    bh(i1,j)=bh(i1,j)*cc(i2,jt)
    100  bh(i1,j+1)=bh(i1,j+1)-ct(i2,j)*bh(i1,j)

    do 200 i=2,i2
    it=i-1
    bh(i,j1)=bh(i,j1)*cc(it,j2)
    200  bh(i+1,j1)=bh(i+1,j1)-cr(it,j2)*bh(i,j1)
    bh(i1,j1)=bh(i1,j1)*cc(i2,j2)

    do 300 j=j1,3,-1
    jt=j-1
    do 350 i=i1,3,-1
    it=i-1
    bh(i,j)=bh(i,j)*cc(it,jt)
    bh(i-1,j)=bh(i-1,j)-cr(it,jt)*bh(i,j)
    350  bh(i,j-1)=bh(i,j-1)-ct(it,jt)*bh(i,j)
    bh(2,j)=bh(2,j)*cc(1,jt)
    300  bh(2,j-1)=bh(2,j-1)-ct(1,jt)*bh(2,j)

    do 400 i=i1,3,-1
    it=i-1
    bh(i,2)=bh(i,2)*cc(it,1)
    400  bh(i-1,2)=bh(i-1,2)-cr(it,1)*bh(i,2)
    bh(2,2)=bh(2,2)*cc(1,1)
end subroutine


! ---------------------------------------------------------------------
subroutine p_lusqelim(cb,cl,cc,cr,ct,b,bh)
    parameter(i3=i2-1,j3=j2-1,i0j0=i0*j0)
    ! real*8 cb,cl,cc,cr,ct,b,bh
    real*8 b,bh
    dimension cb(i2,j2),cl(i2,j2),cc(i2,j2),cr(i2,j2),ct(i2,j2),b(i0,j0),bh(i0,j0)

    call dcopy(i0j0,b,1,bh,1)
    do 100 j=2,j2
    jt=j-1
    do 150 i=2,i2
    it=i-1
    bh(i+1,j)=bh(i+1,j)-cl(it,jt)*bh(i,j)
    150  bh(i,j+1)=bh(i,j+1)-cb(it,jt)*bh(i,j)
    100  bh(i1,j+1)=bh(i1,j+1)-cb(i2,j)*bh(i1,j)

    do 200 i=2,i2
    it=i-1
    200  bh(i+1,j1)=bh(i+1,j1)-cl(it,j2)*bh(i,j1)
    bh(i1,j1)=bh(i1,j1)

    do 300 j=j1,3,-1
    jt=j-1
    do 350 i=i1,3,-1
    it=i-1
    bh(i,j)=bh(i,j)*cc(it,jt)
    bh(i-1,j)=bh(i-1,j)-cr(it,jt)*bh(i,j)
    350  bh(i,j-1)=bh(i,j-1)-ct(it,jt)*bh(i,j)
    bh(2,j)=bh(2,j)*cc(1,jt)
    300  bh(2,j-1)=bh(2,j-1)-ct(1,jt)*bh(2,j)

    do 400 i=i1,3,-1
    it=i-1
    bh(i,2)=bh(i,2)*cc(it,1)
    400  bh(i-1,2)=bh(i-1,2)-cr(it,1)*bh(i,2)
    bh(2,2)=bh(2,2)*cc(1,1)
end subroutine

! ---------------------------------------------------------------------
subroutine p_sipsqelim(cb,cl,cc,cr,ct,b,bh)
    parameter(i3=i2-1,j3=j2-1,i0j0=i0*j0)
    ! real*8 cb,cl,cc,cr,ct,b,bh
    real*8 b,bh
    dimension cb(i2,j2),cl(i2,j2),cc(i2,j2),cr(i2,j2),ct(i2,j2),b(i0,j0),bh(i0,j0)

    call dcopy(i0j0,b,1,bh,1)
    do 100 j=2,j2
    jt=j-1
    do 150 i=2,i2
    it=i-1
    bh(i,j)=bh(i,j)*cc(it,jt)
    bh(i+1,j)=bh(i+1,j)-cl(it,jt)*bh(i,j)
    150  bh(i,j+1)=bh(i,j+1)-cb(it,jt)*bh(i,j)
    bh(i1,j)=bh(i1,j)*cc(i2,jt)
    100  bh(i1,j+1)=bh(i1,j+1)-cb(i2,jt)*bh(i1,j)

    do 200 i=2,i2
    it=i-1
    bh(i,j1)=bh(i,j1)*cc(it,j2)
    200  bh(i+1,j1)=bh(i+1,j1)-cl(it,j2)*bh(i,j1)
    bh(i1,j1)=bh(i1,j1)*cc(i2,j2)

    do 300 j=j1,3,-1
    jt=j-1
    do 350 i=i1,3,-1
    it=i-1
    bh(i-1,j)=bh(i-1,j)-cr(it,jt)*bh(i,j)
    350  bh(i,j-1)=bh(i,j-1)-ct(it,jt)*bh(i,j)
    300  bh(2,j-1)=bh(2,j-1)-ct(1,jt)*bh(2,j)

    do 400 i=i1,3,-1
    400  bh(i-1,2)=bh(i-1,2)-cr(i-1,1)*bh(i,2)
end subroutine

! ----------------------------------------------------------------------
subroutine rep2(ax,ay,bb,cx,cy,rinv,rinv1,dum0,dum1,dum2,dumg,f,x,ie)
    ! ----------------------------------------------------------------------
    !ibm  real*8 rinv,rinv1,dum0,dum1,dum2,x,h,row
    parameter(nbg0=nb0*ngy0,nbg1=nbg0-1,ibir=i2t/ngy0)
    real*8 x,rinv,rinv1,dum1,dum2,dum0,dumg,dd
    dimension ax(i2,j2),ay(i2,j2),bb(i2,j2),cx(i2,j2),cy(i2,0:j2),f(i2,j2),rinv(ibir,i0,nbg0),rinv1(ibir,i0,nbg1),ie(nb0),dum0(i0,nb0),dum1(ibir),dum2(i0),x(i0,j0),dumg(i2t)

    character grd*7
    common/grd/grd


    do 160 ng=0,ngy0-1
    jsp=1
    do 150 nb=1,nb0
    nbg=ng*nb0+nb
    if (mylat .eq. ng) then
	jfp=ie(nb)-2
	do 106 j=jsp,jfp
	do 105 i=1,i2
	105  x(i+1,j+2)=(f(i,j)-ax(i,j)*x(i,j+1)-ay(i,j)*x(i+1,j)-bb(i,j)*x(i+1,j+1)-cx(i,j)*x(i+2,j+1))/cy(i,j)
	call mpi_sendrecv(x(2,j+2),1,mpi_real8,m_w,1,x(i0,j+2),1,mpi_real8,m_e,1,m_cart,istat,ierr)
	call mpi_sendrecv(x(i1,j+2),1,mpi_real8,m_e,1,x(1,j+2),1,mpi_real8,m_w,1,m_cart,istat,ierr)
	106  continue
	if (nbg.eq.nbg0) go to 160
	j=ie(nb)-1
	do 115 i=1,i2
	115  dum2(i)=f(i,j)-ax(i,j)*x(i,j+1)-ay(i,j)*x(i+1,j)-bb(i,j)*x(i+1,j+1)-cx(i,j)*x(i+2,j+1)-cy(i,j)*x(i+1,j+2)
	j=ie(nb)
	call mpi_allgather(dum2,i2,mpi_real8,dumg,i2,mpi_real8,m_comm_lon,ierr)
	do 116 n=1,i0
	dum2(n)=x(n,j)
	116  dum0(n,nb)=x(n,j)
	dd=1.d0
    else
	if (nbg .eq. nbg0) go to 160
	dd=0.d0
    endif
    call mpi_scatter(dumg,ibir,mpi_real8,dum1,ibir,mpi_real8,ng,m_comm_lat,ierr)
    call dgemv('t',ibir,i0,-1.d0,rinv1(1,1,nbg),ibir,dum1,1,dd,dum2,1)
    if (nb .eq. nb0) go to 150
    call mpi_reduce(dum2,x(1,ie(nb)),i0,mpi_real8,mpi_sum,ng,m_comm_lat,ierr)
    jsp=ie(nb)

    150  continue
    call mpi_reduce(dum2,x(1,1),i0,mpi_real8,mpi_sum,ng+1,m_comm_lat,ierr)
    160  continue

    !  wenien
    !       open(160,file='opt/rep160'//grd,form='unformatted')
    !       write(160) x,f
    !       close(160)


    do 300 ng=ngy0-1,0,-1
    do 250 nbs=1,nb0
    nb=nb0-nbs+1
    nbg=nb0*ng+nb

    jsp=1
    if (nb.ne.1) jsp=ie(nb-1)
    if (mylat .eq. ng) then
	if (nbg.eq.nbg0) go to 201
	j=ie(nb)
	do 200 n=1,i0
	200  x(n,j)=dum0(n,nb)
	201  j=ie(nb)-1
	do 210 i=1,i2
	210  dum2(i)=f(i,j)-ax(i,j)*x(i,j+1)-ay(i,j)*x(i+1,j)-bb(i,j)*x(i+1,j+1)-cx(i,j)*x(i+2,j+1)-cy(i,j)*x(i+1,j+2)
	call mpi_allgather(dum2,i2,mpi_real8,dumg,i2,mpi_real8,m_comm_lon,ierr)
    endif

    call mpi_scatter(dumg,ibir,mpi_real8,dum1,ibir,mpi_real8,ng,m_comm_lat,ierr)
    call dgemv('t',ibir,i0,1.d0,rinv(1,1,nbg),ibir,dum1,1,0.d0,dum2,1)
    call mpi_reduce(dum2,dum0(1,nb),i0,mpi_real8,mpi_sum,ng,m_comm_lat,ierr)

    if (mylat .eq. ng) then
	do 220 n=1,i0
	220  x(n,jsp+1)=x(n,jsp+1)+dum0(n,nb)
    endif
    if (nbg.eq.1) go to 300
    if (mylat .eq. ng) then
	do 230 m=1,i2
	230  dum2(m)=dum0(m+1,nb)*cy(m,jsp-1)
	call mpi_allgather(dum2,i2,mpi_real8,dumg,i2,mpi_real8,m_comm_lon,ierr)
    endif
    call mpi_scatter(dumg,ibir,mpi_real8,dum1,ibir,mpi_real8,ng,m_comm_lat,ierr)
    call dgemv('t',ibir,i0,1.d0,rinv1(1,1,nbg-1),ibir,dum1,1,0.d0,dum2,1)
    call mpi_reduce(dum2,dum0(1,nb),i0,mpi_real8,mpi_sum,ng,m_comm_lat,ierr)
    if (mylat .eq. ng) then
	do 245 n=1,i0
	245  dum0(n,nb)=x(n,jsp)+dum0(n,nb)
    endif
    250  continue
    if (mylat .eq. ng) call mpi_send(x(1,2),i0,mpi_real8,m_s,1,m_cart,ierr)
    if (mylat .eq. ng-1) call mpi_recv(x(1,j0),i0,mpi_real8,m_n,1,m_cart,istat,ierr)
    300  continue
    if (mylat .eq. 0) then
	do 350 i=1,i0
	350  dum0(i,1)=0.d0
    endif

    jsp=1
    do 400 nb=1,nb0
    jfp=ie(nb)-2
    if (nb.eq.nb0) jfp=ie(nb)-2
    do 410 i=1,i0
    410  x(i,jsp)=dum0(i,nb)
    do 420 j=jsp,jfp
    do 430 i=1,i2
    430  x(i+1,j+2)=(f(i,j)-ax(i,j)*x(i,j+1)-ay(i,j)*x(i+1,j)-bb(i,j)*x(i+1,j+1)-cx(i,j)*x(i+2,j+1))/cy(i,j)
    call mpi_sendrecv(x(2,j+2),1,mpi_real8,m_w,1,x(i0,j+2),1,mpi_real8,m_e,1,m_cart,istat,ierr)
    call mpi_sendrecv(x(i1,j+2),1,mpi_real8,m_e,1,x(1,j+2),1,mpi_real8,m_w,1,m_cart,istat,ierr)
    420  continue
    400  jsp=ie(nb)
end subroutine


! ----------------------------------------------------------------------
subroutine rep(ax,ay,bb,cx,cy,rinv,rinv1,dum0,dum1,dum2,dumg,f,x,ie)
    parameter(nbg0=nb0*ngy0,nbg1=nbg0-1,ig0=i2*ngx0+2,ig2=ig0-2,ib=ig2/ngy0)
    real*8 x,rinv,rinv1,dum1,dum2,dum0,dumg,dd
    ! real*8 ax,ay,bb,cx,cy,f
    dimension ax(i2,j2),ay(i2,j2),bb(i2,j2),cx(i2,j2),cy(i2,0:j2),f(i2,j2),rinv(ib,i0,nbg0),rinv1(ib,i0,nbg1),ie(nb0),dum0(i0,nb0),dum1(ib),dum2(i0),x(i0,j0),dumg(i2t)

    real*8 timer
    common/timer/timer(10)

    timer(7)=mpi_wtime()

    do 160 ng=0,ngy0-1
    jsp=1
    do 150 nb=1,nb0
    nbg=ng*nb0+nb
    if (mylat .eq. ng) then
	jfp=ie(nb)-2
	do 106 j=jsp,jfp
	do 105 i=1,i2
	105  x(i+1,j+2)=(f(i,j)-ax(i,j)*x(i,j+1)-ay(i,j)*x(i+1,j)-bb(i,j)*x(i+1,j+1)-cx(i,j)*x(i+2,j+1))/cy(i,j)
	timer(8)=mpi_wtime()
	call mpi_sendrecv(x(2,j+2),1,mpi_real8,m_w,1,x(i0,j+2),1,mpi_real8,m_e,1,m_cart,istat,ierr)
	call mpi_sendrecv(x(i1,j+2),1,mpi_real8,m_e,1,x(1,j+2),1,mpi_real8,m_w,1,m_cart,istat,ierr)

	106  timer(5)=timer(5)+mpi_wtime()-timer(8)
	if (nbg.eq.nbg0) go to 160
	j=ie(nb)-1
	do 115 i=1,i2
	115  dum2(i)=f(i,j)-ax(i,j)*x(i,j+1)-ay(i,j)*x(i+1,j)-bb(i,j)*x(i+1,j+1)-cx(i,j)*x(i+2,j+1)-cy(i,j)*x(i+1,j+2)
	j=ie(nb)
	timer(8)=mpi_wtime()
	call mpi_allgather(dum2,i2,mpi_real8,dumg,i2,mpi_real8,m_comm_lon,ierr)
	timer(5)=timer(5)+mpi_wtime()-timer(8)
	do 116 n=1,i0
	dum2(n)=x(n,j)
	116  dum0(n,nb)=x(n,j)
	dd=1.d0
    else
	if (nbg .eq. nbg0) go to 160
	dd=0.d0
    endif
    timer(8)=mpi_wtime()
    call mpi_scatter(dumg,ib,mpi_real8,dum1,ib,mpi_real8,ng,m_comm_lat,ierr)
    timer(5)=timer(5)+mpi_wtime()-timer(8)
    call dgemv('t',ib,i0,-1.d0,rinv1(1,1,nbg),ib,dum1,1,dd,dum2,1)
    if (nb .eq. nb0) go to 150
    timer(8)=mpi_wtime()
    call mpi_reduce(dum2,x(1,ie(nb)),i0,mpi_real8,mpi_sum,ng,m_comm_lat,ierr)
    timer(5)=timer(5)+mpi_wtime()-timer(8)
    jsp=ie(nb)

    150  continue
    timer(8)=mpi_wtime()
    call mpi_reduce(dum2,x(1,1),i0,mpi_real8,mpi_sum,ng+1,m_comm_lat,ierr)
    timer(5)=timer(5)+mpi_wtime()-timer(8)
    160  continue

    do 300 ng=ngy0-1,0,-1
    do 250 nbs=1,nb0
    nb=nb0-nbs+1
    nbg=nb0*ng+nb

    jsp=1
    if (nb.ne.1) jsp=ie(nb-1)
    if (mylat .eq. ng) then
	if (nbg.eq.nbg0) go to 201
	j=ie(nb)
	do 200 n=1,i0
	200  x(n,j)=dum0(n,nb)
	201  j=ie(nb)-1
	do 210 i=1,i2
	210  dum2(i)=f(i,j)-ax(i,j)*x(i,j+1)-ay(i,j)*x(i+1,j)-bb(i,j)*x(i+1,j+1)-cx(i,j)*x(i+2,j+1)-cy(i,j)*x(i+1,j+2)
	timer(8)=mpi_wtime()
	call mpi_allgather(dum2,i2,mpi_real8,dumg,i2,mpi_real8,m_comm_lon,ierr)
	timer(5)=timer(5)+mpi_wtime()-timer(8)
    endif

    timer(8)=mpi_wtime()
    call mpi_scatter(dumg,ib,mpi_real8,dum1,ib,mpi_real8,ng,m_comm_lat,ierr)
    timer(5)=timer(5)+mpi_wtime()-timer(8)
    call dgemv('t',ib,i0,1.d0,rinv(1,1,nbg),ib,dum1,1,0.d0,dum2,1)
    timer(8)=mpi_wtime()
    call mpi_reduce(dum2,dum0(1,nb),i0,mpi_real8,mpi_sum,ng,m_comm_lat,ierr)
    timer(5)=timer(5)+mpi_wtime()-timer(8)

    if (mylat .eq. ng) then
	do 220 n=1,i0
	220  x(n,jsp+1)=x(n,jsp+1)+dum0(n,nb)
    endif
    if (nbg.eq.1) go to 300
    if (mylat .eq. ng) then
	do 230 m=1,i2
	230  dum2(m)=dum0(m+1,nb)*cy(m,jsp-1)
	timer(8)=mpi_wtime()
	call mpi_allgather(dum2,i2,mpi_real8,dumg,i2,mpi_real8,m_comm_lon,ierr)
	timer(5)=timer(5)+mpi_wtime()-timer(8)
    endif
    timer(8)=mpi_wtime()
    call mpi_scatter(dumg,ib,mpi_real8,dum1,ib,mpi_real8,ng,m_comm_lat,ierr)
    timer(5)=timer(5)+mpi_wtime()-timer(8)
    call dgemv('t',ib,i0,1.d0,rinv1(1,1,nbg-1),ib,dum1,1,0.d0,dum2,1)
    timer(8)=mpi_wtime()
    call mpi_reduce(dum2,dum0(1,nb),i0,mpi_real8,mpi_sum,ng,m_comm_lat,ierr)
    timer(5)=timer(5)+mpi_wtime()-timer(8)
    if (mylat .eq. ng) then
	do 245 n=1,i0
	245  dum0(n,nb)=x(n,jsp)+dum0(n,nb)
    endif
    250  continue
    if (mylat .eq. ng) call mpi_send(x(1,2),i0,mpi_real8,m_s,1,m_cart,ierr)
    if (mylat .eq. ng-1) call mpi_recv(x(1,j0),i0,mpi_real8,m_n,1,m_cart,istat,ierr)
    300  continue
    if (mylat .eq. 0) then
	do 350 i=1,i0
	350  dum0(i,1)=0.d0
    endif

    jsp=1
    do 400 nb=1,nb0
    jfp=ie(nb)-2
    if (nb.eq.nb0) jfp=ie(nb)-2
    do 410 i=1,i0
    410  x(i,jsp)=dum0(i,nb)
    do 420 j=jsp,jfp
    do 430 i=1,i2
    430  x(i+1,j+2)=(f(i,j)-ax(i,j)*x(i,j+1)-ay(i,j)*x(i+1,j)-bb(i,j)*x(i+1,j+1)-cx(i,j)*x(i+2,j+1))/cy(i,j)
    timer(8)=mpi_wtime()
    call mpi_sendrecv(x(2,j+2),1,mpi_real8,m_w,1,x(i0,j+2),1,mpi_real8,m_e,1,m_cart,istat,ierr)
    call mpi_sendrecv(x(i1,j+2),1,mpi_real8,m_e,1,x(1,j+2),1,mpi_real8,m_w,1,m_cart,istat,ierr)
    420  timer(5)=timer(5)+mpi_wtime()-timer(8)
    400  jsp=ie(nb)
end subroutine


end module
