program laplace
code:  LaPlace solver sample program 
code:    From "getting Started in CM Fortran", revised for C90 CF77
code:    Solution of LaPlace's equation on unit square.
code:    BC: f(x,1)=1, f(0,y)=f(x,0)=f(1,y)=2, 0<{x,y}<1.
code:    Initial interior starting iterate: f(x,y) = 0.
      parameter(maxx=32,maxy=maxx)
      real f(maxx,maxy),df(maxx,maxy),df2(maxx,maxy),dfmax(maxy)
      real rms_error,max_error
      logical cmask(maxx,maxy)
      integer iteration
      intrinsic sum,maxval,cshift
comment:  Initialize the mask for the interior points
      cmask = .false.
      cmask(2:maxx-1,2:maxy-1) = .true.
comment:  Initialize f
      f = 2.
      f(:,maxy) = 1.
      where (cmask) f = 0.
      max_error  = 1
      iteration = 0
comment:  Iterate until max_error < 1.e-3 stopping criterion.
      do while (max_error.gt.1.e-3)
         iteration = iteration + 1
comment:  Compute df, the change at each iteration, and udate.
         df = 0.
         where (cmask)
	    df = 0.25*(cshift(f,1,1) + cshift(f,-1,1)
     &         + cshift(f,1,2) + cshift(f,-1,2)) - f
Caution: CM has cshift args in opposite order: (f,dim,shift), not (f,s,d)
C	    df = 0.25*(cshift(f,1,1) + cshift(f,1,-1)
C     &         + cshift(f,2,1) + cshift(f,2,-1)) - f
	    f = f + df
         endwhere
comment:  Compute RMS and Maximum errors.
         df2= df*df
         rms_error = sqrt(sum(df2)/((maxx-2)*(maxy-2)))
         max_error = maxval(df)
comment:  Output errors every 20th iteration
         if(mod(iteration,20).eq.0) then
	    write(6,*) 'every 20th: ',iteration,rms_error,max_error
         endif
      enddo
comment:  Output final iteration count and errors (stopping criteria).
      write(6,*) 'final: ',iteration,rms_error,max_error
      stop
      end