PROGRAM poisson_solver

        integer nr,nc
        integer lshsec,ldwork,lshaec

        open(unit=52,status='old')
        read(unit=52,fmt=*) nc,nr

        lwork=nr*nc*10
        lshsec=nr*nc*10
        ldwork=nr+1
        lshaec=nr*nc*10

	call poisson_solver1(nc,nr,lwork,lshsec,ldwork,lshaec)

	end

C -------------------------------------------------------------

	subroutine poisson_solver1(nc,nr,lwork,lshsec,ldwork,lshaec)

	IMPLICIT NONE

	integer i,j,nr,nc,nlat,nlon,isym,nt,ids,jds,lwork,ierror
	integer lshsec,ldwork,lshaec,ndab,mdab

	real rhs(nr,nc),xlmbda,sf(nr,nc),work(lwork),pertrb
	real wshsec(lshsec),wshaec(lshaec),a(nc,nr),b(nc,nr)
	double precision dwork(ldwork)

C Read in RHS
C Note that SPHEREPACK assumes (rows,cols) whereas PP uses (cols,rows)
	open(unit=52,status='old')
	do i=1,nr
C Inner loop is along the row, because that natural for PP,
	do j=1,nc
  	  read(unit=52,fmt=*) rhs(i,j)
	enddo
	enddo

	nlat=nr
	nlon=nc
	isym=0
	nt=1
	xlmbda=0
	ids=nlat
	jds=nlon
	ndab=nlat
	mdab=nlon

	call shaeci(nlat,nlon,wshaec,lshaec,dwork,ldwork,ierror)
	print *,'shaeci err: ',ierror

	call shseci(nlat,nlon,wshsec,lshsec,dwork,ldwork,ierror)
	print *,'shseci err: ',ierror

	call shaec(nlat,nlon,isym,nt, RHS , NLAT , NLON ,a,b,mdab,ndab,
     &        wshaec,lshaec,work,lwork,ierror)
	print *,'shaec err: ',ierror
	
	call islapec(nlat,nlon,isym,nt,xlmbda,sf,ids,jds,a,b,
     &               mdab,ndab,wshsec,lshsec,work,lwork,pertrb,ierror)

	print *,'islapec err code: (pert) ',ierror,'(',pertrb,')'

	open(unit=51)
	do i=1,nr
	do j=1,nc
  	  write(unit=51,fmt=*) sf(i,j)
	enddo
	enddo
	
	end