PROGRAM laplacian

	integer nc,nr
	integer lshsec,ldwork,lshaec,lwork

	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 laplacian1(nc,nr,lwork,lshsec,ldwork,lshaec)

	end

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

	SUBROUTINE laplacian1(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)

	do i=1,nr
	do j=1,nc
  	  read(unit=52,fmt=*) sf(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, SF , NLAT , NLON ,a,b,mdab,ndab,
     &        wshaec,lshaec,work,lwork,ierror)
	print *,'shaec err: ',ierror
	
	call slapec(nlat,nlon,isym,nt,RHS,ids,jds,a,b,mdab,ndab,
     &                    wshsec,lshsec,work,lwork,ierror)

	print *,'slapec err code: ',ierror

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