      PROGRAM SPPOT                                                                                                                 
C     PROGRAM SPPOT (SPIN,       SPOUT,       OUTPUT,                   )       E2                                                  
C    1         TAPE1=SPIN, TAPE2=SPOUT, TAPE6=OUTPUT)                                                                               
C     -----------------------------------------------                           E2                                                  
C                                                                               E2                                                  
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)       E2                                                  
C     FEB 15/94 - F.MAJAESS (REVISE FOR TRUNCATIONS > 99)                                                                           
C     JAN 29/92 - E. CHAN  (CONVERT HOLLERITH LITERALS TO ASCII)                                                                    
C     SEP   /86 - G.J.BOER                                                                                                          
C                                                                               E2                                                  
CSPPOT   - PRODUCES THE SPHERICAL HARMONIC "POTENTIAL FUNCTION"         1  1    E1                                                  
C                                                                               E3                                                  
CAUTHOR  - G.J.BOER                                                             E3                                                  
C                                                                               E3                                                  
CPURPOSE - PRODUCES THE SPHERICAL HARMONIC "POTENTIAL FUNCTION" FROM THE        E3                                                  
C          (M,N) SPACE ANALOG TO THE POISSON EQUATION. THE METHOD USED          E3                                                  
C          IS TO COPY THE TRIANGLE TO A SQUARE AND TO SOLVE THE POISSON         E3                                                  
C          EQUATION USING COSINE SERIES.                                        E3                                                  
C          NOTE - CONSIDERABLE CARE IS NEEDED IN THE FORM OF THE                E3                                                  
C                 QUADRATURE.                                                   E3                                                  
C                 IT WORKS ONLY FOR TRIANGULAR TRUNCATION .LT. $MP1$            E3                                                  
C                                                                               E3                                                  
CINPUT FILE...                                                                  E3                                                  
C                                                                               E3                                                  
C      SPIN  = CONTAINS SERIES OF THE ONE-SIDED "INTERACTION TERM"  J,          E3                                                  
C              I, OR M WHICH IS THE OUTPUT OF DECKS SPJ, SPI, SPM...            E3                                                  
C                                                                               E3                                                  
COUTPUT FILE...                                                                 E3                                                  
C                                                                               E3                                                  
C      SPOUT = SPHERICAL HARMONIC "POTENTIAL FUNCTION" OF SPIN                  E3                                                  
C----------------------------------------------------------------------------                                                       
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      LOGICAL OK                                                                                                                    
      INTEGER P,Q                                                                                                                   
      DIMENSION LSR(2,66),G(66,66),GT(66,66),                                                                                       
     1          COSN(66,66)                                                                                                         
      COMMON/LCM/F                                                                                                                  
      COMPLEX F(2210)                                                                                                               
      COMMON/BUFCOM/IBUF(8),IDAT(4420)                                                                                              
      DATA MAXX/4420/                                                                                                               
C-----------------------------------------------------------------------                                                            
      NFF=3                                                                                                                         
      CALL JCLPNT(NFF,1,2,6)                                                                                                        
      NRECS=0                                                                                                                       
      REWIND 1                                                                                                                      
      REWIND 2                                                                                                                      
                                                                                                                                    
  100 CALL GETFLD2(1,F,NC4TO8("SPEC"),0,0,0,IBUF,MAXX,OK)                                                                           
      IF(.NOT.OK)THEN                                                                                                               
        IF(NRECS.EQ.0) CALL                        XIT('SPPOT',-1)                                                                  
        WRITE(6,6020) NRECS                                                                                                         
        CALL                                       XIT('SPPOT',0)                                                                   
      ENDIF                                                                                                                         
      IF(NRECS.EQ.0) CALL PRTLAB (IBUF)                                                                                             
      CALL DIMGT(LSR,LA,LR,LM,KTR,IBUF(7))                                                                                          
C                                                                                                                                   
C     * GET THE INPUT DATA INTO THE PROPER SQUARE                                                                                   
C     * DO NOT DIVIDE M>0 VALUES BY 2 - KEEP ONESIDED INTERACTION TERM                                                              
C                                                                                                                                   
      DO 210 M=1,LM                                                                                                                 
         DO 210 N=M,LM                                                                                                              
         G(M,N)=REAL( F(LSR(1,M)+N-M) )                                                                                             
         G(N,M)=G(M,N)                                                                                                              
  210 CONTINUE                                                                                                                      
C                                                                                                                                   
C     * CALCULATE COSINE SERIES WHICH HAS THE DESIRED SYMMETRY                                                                      
C                                                                                                                                   
      LM1=LM-1                                                                                                                      
C                                                                                                                                   
      D=3.1415926536E0/FLOAT(LM1)                                                                                                   
C                                                                                                                                   
      DO 220 P=1,LM                                                                                                                 
      DO 220 M=1,LM                                                                                                                 
  220 COSN(M,P)=COS(D*FLOAT((M-1)*(P-1)))                                                                                           
C                                                                                                                                   
C     * TRANSFORM THE INPUT FUNCTION G                                                                                              
C                                                                                                                                   
      FAC=4.0E0/FLOAT(LM1*LM1)                                                                                                      
C                                                                                                                                   
      DO 230 P=1,LM                                                                                                                 
      DO 230 Q=1,LM                                                                                                                 
C                                                                                                                                   
      GT(P,Q)=0.0E0                                                                                                                 
C                                                                                                                                   
         DO 240 M=1,LM                                                                                                              
         DO 240 N=1,LM                                                                                                              
C                                                                                                                                   
         AA=G(M,N)*FAC*COSN(M,P)*COSN(N,Q)                                                                                          
C                                                                                                                                   
         IF(M.EQ.1.OR.M.EQ.LM) AA=0.5E0*AA                                                                                          
         IF(N.EQ.1.OR.N.EQ.LM) AA=0.5E0*AA                                                                                          
C                                                                                                                                   
  240  GT(P,Q)=GT(P,Q)+AA                                                                                                           
C                                                                                                                                   
       IF(P.EQ.1.OR.Q.EQ.1) GT(P,Q)=0.5E0*GT(P,Q)                                                                                   
C                                                                                                                                   
  230 CONTINUE                                                                                                                      
C                                                                                                                                   
C     * DIVIDE COEFFICIENTS BY P2+Q2 TO GET LAPLACIAN                                                                               
C                                                                                                                                   
      GT(1,1)=0.0E0                                                                                                                 
C                                                                                                                                   
      DO 235 P=1,LM                                                                                                                 
      DO 235 Q=1,LM                                                                                                                 
  235 IF(P+Q.GT.2) GT(P,Q)=GT(P,Q)/FLOAT((P-1)*(P-1)+(Q-1)*(Q-1))                                                                   
C                                                                                                                                   
C                                                                                                                                   
C     * TRANSFORM THE RESULTING POTENTIAL FUNCTION BACK                                                                             
C                                                                                                                                   
       DO 250 M=1,LM                                                                                                                
       DO 250 N=1,LM                                                                                                                
C                                                                                                                                   
       G(M,N)=0.0E0                                                                                                                 
C                                                                                                                                   
              DO 260 P=1,LM                                                                                                         
              DO 260 Q=1,LM                                                                                                         
  260         G(M,N)=G(M,N)+GT(P,Q)*COSN(M,P)*COSN(N,Q)                                                                             
C                                                                                                                                   
  250 CONTINUE                                                                                                                      
C                                                                                                                                   
C     * PUT THE DATA BACK IN THE USUAL VECTOR AND CHANGE SIGN                                                                       
C                                                                                                                                   
      DO 270 M=1,LM                                                                                                                 
      DO 270 N=M,LM                                                                                                                 
      INDEX=LSR(1,M)+N-M                                                                                                            
  270 F(INDEX)=-CMPLX(G(M,N))                                                                                                       
C                                                                                                                                   
      CALL PUTFLD2(2,F,IBUF,MAXX)                                                                                                   
      IF(NRECS.EQ.0) CALL PRTLAB (IBUF)                                                                                             
      NRECS=NRECS+1                                                                                                                 
      GO TO 100                                                                                                                     
C-----------------------------------------------------------------------                                                            
6020  FORMAT('0SPPOT COMPUTED ',I5,' RECORDS')                                                                                      
      END                                                                                                                           
