C
C   Copyright (c) 1998 Silvano Bonazzola
C
C    This file is part of LORENE.
C
C    LORENE is free software; you can redistribute it and/or modify
C    it under the terms of the GNU General Public License as published by
C    the Free Software Foundation; either version 2 of the License, or
C    (at your option) any later version.
C
C    LORENE is distributed in the hope that it will be useful,
C    but WITHOUT ANY WARRANTY; without even the implied warranty of
C    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C    GNU General Public License for more details.
C
C    You should have received a copy of the GNU General Public License
C    along with LORENE; if not, write to the Free Software
C    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
C
C
C
C $Id: grgp2s.f,v 1.2 2012/03/30 12:12:43 j_novak Exp $
C $Log: grgp2s.f,v $
C Revision 1.2  2012/03/30 12:12:43  j_novak
C Cleaning of fortran files
C
C Revision 1.1.1.1  2001/11/20 15:19:30  e_gourgoulhon
C LORENE
C
c Revision 1.1  1998/06/22  10:38:29  eric
c Initial revision
c
C
C $Header: /cvsroot/Lorene/F77/Source/Poisson2d/grgp2s.f,v 1.2 2012/03/30 12:12:43 j_novak Exp $
C
C

	SUBROUTINE GRGP2S(NDEG,NDR,NDT,IPAR,CC,DEN1,BB,SOLH,CDEN)
C
C## Routine modifiee le 20.02.1995 : suppression des variables inutiles
C
	IMPLICIT NONE

C		ROUTINE POUR LA SOLUTION DE L'EQUATION DE POISSON EN
C		COORDONNEES SPHERIQUES EN GENERAL ET DANS LE CAS OU IL 
C		EXISTERAIT UNE SYMETRIE PAR RAPPORT AU PLAN EQUATORIAL
C		z=0, OU DANS LE CAS SUPERSYMETRIQUE. LES CONDITIONS AU
C		CONTOUR SONT CELLES IMPOSEES PAR UN POTENTIEL GRAVITA-
C		TIONNEL OU ELECTROSTATIQUE. (CORDONNEES RAREFIEES A
C		L'ORIGINE)
C
C	N.B.	AUCUN DESALIASAGE EST EFFECTUE. LE CAS IPAR=3 N'A PAS ETE TESTE
C       ----
C
C	ARGUMENTS DE LA ROUTINE:
C
C		NDEG	=TABLEAU CONTENANT LE NOMBRE DE DEGREES DE LIBERTE'
C			 POUR LES COORDONNES r, POUR TETA, ET FI
C		NDR	=DIMENSION DU PREMIER INDICE DES TABLEAUX DEN ET
C			 DEN1.
C		NDT  = COMME NDR MAIS POUR LE 2ME INDICE. NDR ET NDT
C			 SONT DEFINIS DANS LE PROGRAMME APPELANT LA ROUTINE.
C		IPAR	=DRAPEAU: 
C
C		         IPAR=0 AUCUNE SYMETRIE
C
C			 IPAR=2, LA SOLUTION CHERCHEE EST SCALAIRE AVEC
C                        UNE SYMETRIE PAR RAPPORT AU PLAN EQUATORIALE z=0.
C			 IPAR=3, LA SOLUTION CHERCHEE EST ANTISYMETRIQUE
C			 PAR RAPPORT LE PLAN z=0
C
C			 IPAR=4, CAS SUPERSYMETRIQUE (SYMETRIE PAR RAPORT
C		         LE PLAN z=0 E PAR RAPORT LA TRANSFORMATION x,y -> -x,y
C			 IPAR=5, CAS SUPERSYMETRIQUE ANTISYMETRIQUE PAR
C			 RAPPORT LE PLAN z=0
C
C		DEN1,SOLH =TABLEAUX DE TRAVAIL. LES DIMENSIONS DE CES TABLEAUX
C			     DOIVENT ETRE SUPERIEUR A NR1*NT1*MAX(2,NF)
C			     LES DEGREES DE LIBERTE' DE LA VARIABLE FI, ET
C			     NR1 CEUX DE LA VARIABLE r.
C		BB	=TABLEAU DE TRAVAIL DE DIMENSION .GE.(NR1+1)*12
C		CDEN	=TABLEAU CONTENANT EN IMPUT LES COEFFICIENTS DE LA 
C			 DENSITE' (FOURIER EN FI (TROISIEME INNDICE,LEGENDRE
C		         EN TETA(2me INDICE) ET TCHEBITCHEF EN r (1ere INDICE)
C			 DEN(LR,LY,LF) ET EN OUTPUT LA SOLUTION.
C			 ETRE AU MOINS DEN(NR1+1,NT1,NF)
C
C			 LE STOCKAGE DES COEFFICIENTS EST LE MEME QUE DANS
C			 FCIR3S POUR IDR=7. C'EST A DIRE: POUR LE 3me
C			 INDICE( PARTIE EN PHI) DANS LF=1 IL-Y-A
C			 LE COEFFICIENT CORRESPONDENT A LA FREQUENCE ZERO,
C			 DANS LF =2,3 LES COEFFICIENTS CORRESPONDANT A LA
C			 FREQUENCE 1 (COS ET SIN) ET AINSI DE SUITE.
C			 LE STOCAGE DE FONCTIONS ASSOCIEES DE LEGENDRE
C			 EST LE SUIVANT: LE COEFFICIENT CORRESPONDENT A LA
C			 FONCTION ASSOCIEE DE LEGENDRE
C                         m
C	                P (theta,fi) SE TROUVE DANS DEN(j+1,LF).
C                         m+j	    
C
C		ROUTINE testee le 29/10/94
C
	character*120 header
	data header/'$Header: /cvsroot/Lorene/F77/Source/Poisson2d/grgp2s.f,v 1.2 2012/03/30 12:12:43 j_novak Exp $'/

	double PRECISION
     1	DEN1,CDEN,BB,CC,SOLH,VH1,VA1,DE1,RAP	
C
	INTEGER NDR,NDT,LR,LY,LF,NR1,NR,NF,NDEG
     1	,N257,NT1,LT,IPA2,LF2,LF21
     1	,LEL,LELLE,J,MM,ILT,IPAR,JF,LF1,IPA,IPP
C
	PARAMETER (N257=260)
C
	DIMENSION CC(*),CDEN(NDR,NDT,*),ILT(N257)
	DIMENSION DEN1(NDR,*),BB(NDR,*),SOLH(NDR,NDT,*),NDEG(3)
C
	NR1=NDEG(1)
	NT1=NDEG(2)
	NF= NDEG(3)
	IF(NR1+2.GT.N257.OR.NF.GT.N257) THEN
	PRINT*,'DIMENSIONS INSUFF. DANS LA SUB. GRGP2S'
	PRINT*,'NR1=,NT1=,NF=',NR1,NT1,NF
	STOP
	ENDIF
C
	IF(NR1+2.GT.NDR.OR.NT1.GT.NDT) THEN
	PRINT*,'DIMENSIONS INSUFF. DANS LA SUB. GRGP2S'
	PRINT*,'NR1=,NT1=,NF=',NR1,NT1,NF
	PRINT*,'NDR,NDT=',NDR,NDT
	STOP
	ENDIF
C
	IF(IPAR.EQ.3) THEN
	PRINT*,'ROUTINE GRGP2S: CE CAS EST IMPLEMENTE MAIS PAS TESTE'
	PAUSE
	ENDIF
C
	NR=NR1-1
C
C.............................................................................
C
C			LELLE REPRESENTE LE NOMBRE QUANTIQUE l.
	IPA=1
	IPA2=IPAR/2
	IF((IPAR/2)*2.EQ.IPAR) IPA=0
C
C		INVERSION DE L'OPERATEUR
C
C		CAS SUPERSYMETRIQUE, IPAR=4,5
C
	IF(IPA2.EQ.2) THEN

	DO 3 LY=1,NT1
	LEL=LY+LY-1
	LELLE=LEL-1+IPA
	LF1=LEL
	IF(LF1.GT.NF) LF1=NF
C
C		RELATION ENTRE LES VALEURS DE LF ET LE NOMBRE QUANTIQUE l
C
	DO 2 LF=1,LF1
	MM=(LF/2)*2
	LT=(LEL-MM)/2+1
C
	ILT(LF)=LT
  2	CONTINUE
	CC(2)=1
	CC(3)=-LELLE**2
C
	DO 4 LF=1,LF1
	DO 5 LR=1,NR1
	DEN1(LR,LF)=CDEN(LR,ILT(LF),LF)
   5	CONTINUE
   4	CONTINUE
C
	CALL ILGPRS(NR,NDR,LF1,IPA,CC,BB,DEN1,SOLH)
C
	DO 9 LF=1,LF1
	DO 10 LR=1,NR1
	CDEN(LR,ILT(LF),LF)=SOLH(LR,LF,1)
   10	CONTINUE
    9	CONTINUE
    3	CONTINUE
C
	DO LF=2,NF
	MM=(LF/2)-1
	DO LY=NT1-MM,NT1
	DO LR=1,NR1
	CDEN(LR,LY,LF)=0
	ENDDO
	ENDDO
	ENDDO
C
C	PREPARATION SOLUTIONS HOMGENES
C
	DO LR=1,NR1
	CC(LR)=0
	ENDDO
	CC(1)=2
	IF(IPA.EQ.1) CC(1)=-1
C
	DO LY=1,NT1
	DO LR=1,NR1
	SOLH(LR,LY,1)=CC(LR)
	ENDDO
	CALL DIRA2S(NR,IPA,1,CC,BB)
	DO LR=1,NR1
	CC(LR)=BB(LR,1)
	ENDDO
	ENDDO
C
C		RACCORDEMENT AVEC LE VIDE CAS SUPERSYMETRIQUE,IPAR=4,5
C
	DO LF=1,NF
	MM=(LF/2)
	DO LY=1,NT1-MM
	LEL=LY+LY-1+IPA+MM*2
	LELLE=LEL-1
C
C@@@
	IF(LELLE.GT.0.OR.IPA.EQ.1) THEN
	LT=LELLE/2+1
C	VH1=LELLE+LEL
	VH1=LELLE+LELLE
	DO LR=1,NR1
	CC(LR)=CDEN(LR,LY,LF)	
	ENDDO
C
	CALL EXTR1S(NR,1,0,IPA,CC,VA1)
	CALL EXTR1S(NR,1,1,IPA,CC,DE1)
C
	RAP=(VA1*LELLE+DE1)/VH1
	DO LR=1,NR1
	CDEN(LR,LY,LF)=CDEN(LR,LY,LF)-RAP*SOLH(LR,LT,1)
	ENDDO
	ENDIF
C
	IF(IPA.EQ.1) CDEN(NR1,LY,LF)=0	
	ENDDO
	ENDDO
	RETURN
	ENDIF
C
C		CAS AVEC SYMETRIES PAR RAPPORT LE PLAN z=0, MAIS SENS
C		SYMETRIES EN phi, IPAR=2,3
C
	IF(IPA2.EQ.1) THEN
C
	DO 777 LY=1,NT1
	LEL=LY+LY-2
	LELLE=LEL+IPA
C
	DO LR=1,NR1
	DEN1(LR,1)=CDEN(LR,LY,1)
	ENDDO
C
	CC(2)=1
	CC(3)=-LELLE**2
	CALL ILGPRS(NR,NDR,1,IPA,CC,BB,DEN1,SOLH)
C
	DO LR=1,NR1
	CDEN(LR,LY,1)=SOLH(LR,1,1)
	ENDDO
C
	IF(NF.LT.4) GO TO 776
C
	ILT(1)=1
	LF2=1
	LF1=MIN0(LEL+LEL+1,NF)
	IF(LF1.GE.4) THEN
	DO LF=4,LF1,4	
	MM=LF/2
	LT=LY-LF/4
	DO JF=LF,MIN0(LF+1,LF1)
	LF2=LF2+1
	ILT(LF2)=LT
	DO LR=1,NR1
	DEN1(LR,LF2)=CDEN(LR,LT,JF)
	ENDDO
	ENDDO
	ENDDO
C
	CC(2)=1
	CC(3)=-LELLE**2
C
	CALL ILGPRS(NR,NDR,LF2,IPA,CC,BB,DEN1,SOLH)
C
	LF2=1
	DO LF=4,LF1,4
	DO JF=LF,MIN0(LF+1,LF1)
	LF2=LF2+1
	DO LR=1,NR1
	CDEN(LR,ILT(LF2),JF)=SOLH(LR,LF2,1)
	ENDDO
	ENDDO
	ENDDO
	ENDIF
	LF21=LF2+1
C
  776	CONTINUE
  777 	CONTINUE
C
	IF(NF.GE.2) THEN
C
C		CAS l-ipa impaire
C			
	DO 778 LT=1,NT1	
	LEL=LT+LT-1
	LELLE=LEL+IPA
	LF1=MIN0(LEL+LEL+1,NF-1)
	LF2=0
	DO LF=2,LF1,4
	MM=LF/2
	LY=LT-LF/4
	DO JF=LF,MIN0(LF+1,LF1)
	LF2=LF2+1
	ILT(LF2)=LY
	DO LR=1,NR1
	DEN1(LR,LF2)=CDEN(LR,LY,JF)
	ENDDO
	ENDDO
	ENDDO
	CC(2)=1
	CC(3)=-LELLE**2
	CALL ILGPRS(NR,NDR,LF2,(1-IPA),CC,BB,DEN1,SOLH)
	LF2=0
	DO LF=2,LF1,4
	DO JF=LF,MIN0(LF+1,LF1)
	LF2=LF2+1
	DO LR=1,NR1
	CDEN(LR,ILT(LF2),JF)=SOLH(LR,LF2,1)
	ENDDO
	ENDDO
	ENDDO
  778	CONTINUE
	ENDIF
C
C		CALCUL DE LASOLUTION HOMOGENE (POUR l-ipa paire)
C
	DO JF=1,2
C
	DO LR=1,NR1
	CC(LR)=0
	ENDDO
C
	LF2=IPA-JF+1
	IF((LF2/2)*LF2.EQ.LF2) THEN
	LF2=0
	CC(1)=2
	ELSE
	LF2=1
	CC(1)=-1
	ENDIF
C
	DO LY=1,NT1
	DO LR=1,NR1
	SOLH(LR,LY,JF)=CC(LR)
	ENDDO
	CALL DIRA2S(NR,LF2,1,CC,BB)
	DO LR=1,NR1
	CC(LR)=BB(LR,1)
	ENDDO	
	ENDDO
	ENDDO
C
C		RACCORDEMENT AVEC LE VIDE CAS SYMETIQUE PAR RAPPORT z=0
C		IPAR=2,3
C
C@@
	DO LT=1,NT1
	LEL=LT+LT-1+IPA
	VH1=LEL+LEL-1
C
	DO LR=1,NR1
	CC(LR)=CDEN(LR,LT,1)
	ENDDO
	CALL EXTR1S(NR,1,0,IPA,CC,VA1)
	CALL EXTR1S(NR,1,1,IPA,CC,DE1)
	RAP=(VA1*LELLE+DE1)/VH1	
	DO LR=1,NR1
	CDEN(LR,LT,1)=CDEN(LR,LT,1)-RAP*SOLH(LR,LT,1+IPA)
	ENDDO
	ENDDO
C
	IF(NF.EQ.1) RETURN
C
	IPP=1
	IF(IPA.EQ.1) IPP=0
	DO LF=2,NF,4
	MM=LF/2
	DO JF=LF,MIN0(LF+1,NF)
	DO LT=1,NT1-MM/2
	LELLE=LT+LT-2+IPA+MM	
	LEL=LELLE/2+1
	VH1=LELLE+LELLE+1
C
	DO LR=1,NR1
	CC(LR)=CDEN(LR,LT,JF)
	ENDDO
	CALL EXTR1S(NR,1,0,IPP,CC,VA1)
	CALL EXTR1S(NR,1,1,IPP,CC,DE1)
	RAP=(VA1*LELLE+DE1)/VH1	
	DO LR=1,NR1
	CDEN(LR,LT,JF)=CDEN(LR,LT,JF)-RAP*SOLH(LR,LEL,2-IPA)
	ENDDO
	ENDDO
	ENDDO
	ENDDO
C
	IF(NF.LT.4) RETURN
C
	DO LF=4,NF,4
	MM=LF/2
	DO JF=LF,MIN0(LF+1,NF)
	DO LT=1,NT1-MM/2
	LELLE=LT+LT-2+IPA+MM	
	VH1=LELLE+LELLE+1
	LEL=LELLE/2+1
C
	DO LR=1,NR1
	CC(LR)=CDEN(LR,LT,JF)
	ENDDO
	CALL EXTR1S(NR,1,0,IPA,CC,VA1)
	CALL EXTR1S(NR,1,1,IPA,CC,DE1)
	RAP=(VA1*(LELLE+1)+DE1)/VH1	
	DO LR=1,NR1
	CDEN(LR,LT,JF)=CDEN(LR,LT,JF)-RAP*SOLH(LR,LEL,1)
	ENDDO
	ENDDO
	ENDDO
	ENDDO
C
	RETURN
	ENDIF
C
C		CAS IPAR=0: AUCUNE SYMETRIE
C
	IF(IPAR.EQ.0) THEN
	DO LEL=1,NT1
	LELLE=LEL-1
	LF1=MIN0(NF,LEL*2-1)
C
C		RELATION ENTRE LES VALEURS DE LF ET LE NOMBRE QUANTIQUE l
C
	DO LF=1,LF1
	LF2=LF/2
	LT=LEL-LF2
	IF((LF2/2)*2.NE.LF2) LT=LT+1
	ILT(LF)=LT
	ENDDO
C
	IPA=1
	IF((LELLE/2)*2.EQ.LELLE) IPA=0
C
	CC(1)=1
	CC(2)=1
	CC(3)=-LELLE**2
C
	DO LF=1,LF1
	DO LR=1,NR1
	DEN1(LR,LF)=CDEN(LR,ILT(LF),LF)
	ENDDO
	ENDDO
C
	CALL ILGPRS(NR,NDR,LF1,IPA,CC,BB,DEN1,SOLH)
C		
	DO  LF=1,LF1
	DO  LR=1,NR1
	CDEN(LR,ILT(LF),LF)=SOLH(LR,LF,1)
	ENDDO
	IF(IPA.EQ.1) CDEN(NR1,ILT(LF),LF)=0
	ENDDO
	ENDDO
C
	
	DO LF=2,NF
	MM=(LF/2)-1
	DO LY=NT1-MM,NT1
	DO LR=1,NR1
	CDEN(LR,LY,LF)=0
	ENDDO
	ENDDO
	ENDDO
C
C	PREPARATION SOLUTIONS HOMGENES
C
	DO LR=1,NR1
	CC(LR)=0
	ENDDO
	CC(1)=2
	IF(IPA.EQ.1) CC(1)=-1
C
	DO LY=1,NT1,2
	DO LR=1,NR1
	SOLH(LR,LY,1)=CC(LR)
	ENDDO
	CALL DIRA2S(NR,0,1,CC,BB)
	DO LR=1,NR1
	CC(LR)=BB(LR,1)
	ENDDO
	ENDDO
C
	DO LR=1,NR1
	CC(LR)=0
	ENDDO
C
	CC(1)=1
C
	DO LY=2,NT1,2
	DO LR=1,NR1
	SOLH(LR,LY,1)=CC(LR)
	ENDDO
	CALL DIRA2S(NR,1,1,CC,BB)
	DO LR=1,NR1
	CC(LR)=BB(LR,1)
	ENDDO
	ENDDO	
C
C		RACCORDEMENT AVEC LE VIDE CAS SENS SYMETRIES , IPAR=0
C
	DO LY=2,NT1,2
	DO LR=1,NR1
	SOLH(LR,LY,1)=-SOLH(LR,LY,1)
	ENDDO
	ENDDO
C
	DO LY=2,NT1
	LELLE=LY-1
	VH1=2*LY-2
	DO LR=1,NR1
	CC(LR)=CDEN(LR,LY,1)	
	ENDDO
C
	IF((LELLE/2)*2.EQ.LELLE) THEN
	CALL EXTR1S(NR,1,0,0,CC,VA1)
	CALL EXTR1S(NR,1,1,0,CC,DE1)
	ELSE
	CALL EXTR1S(NR,1,0,1,CC,VA1)
	CALL EXTR1S(NR,1,1,1,CC,DE1)
	ENDIF
C
	RAP=(VA1*LELLE+DE1)/VH1
	DO LR=1,NR1
	CDEN(LR,LY,1)=CDEN(LR,LY,1)-RAP*SOLH(LR,LY,1)
	ENDDO
	ENDDO
C
	DO LF=2,NF,4
	MM=(LF/2)
	DO J=LF,MIN0(LF+1,NF)
	DO LY=1,NT1-MM
	LEL=LY+MM-1
	LELLE=LEL-1
	VH1=2*LEL-1
	DO LR=1,NR1
	CC(LR)=CDEN(LR,LY,J)	
	ENDDO
C
	IF((LELLE/2)*2.EQ.LELLE) THEN
	CALL EXTR1S(NR,1,0,0,CC,VA1)
	CALL EXTR1S(NR,1,1,0,CC,DE1)
	ELSE
	CALL EXTR1S(NR,1,0,1,CC,VA1)
	CALL EXTR1S(NR,1,1,1,CC,DE1)
	ENDIF
C
	RAP=(VA1*LELLE+DE1)/VH1
	DO LR=1,NR1
	CDEN(LR,LY,J)=CDEN(LR,LY,J)-RAP*SOLH(LR,LEL,1)
	ENDDO
	ENDDO
	ENDDO
	ENDDO
C
	DO LF=4,NF,4
	MM=(LF/2)
	DO J=LF,MIN0(LF+1,NF)
	DO LY=1,NT1-MM
	LEL=LY+MM
	LELLE=LEL-1
	VH1=2*LEL-1
	DO LR=1,NR1
	CC(LR)=CDEN(LR,LY,J)	
	ENDDO
C
	IF((LELLE/2)*2.EQ.LELLE) THEN
	CALL EXTR1S(NR,1,0,0,CC,VA1)
	CALL EXTR1S(NR,1,1,0,CC,DE1)
	ELSE
	CALL EXTR1S(NR,1,0,1,CC,VA1)
	CALL EXTR1S(NR,1,1,1,CC,DE1)
	ENDIF
C
	RAP=(VA1*LELLE+DE1)/VH1
	DO LR=1,NR1
	CDEN(LR,LY,J)=CDEN(LR,LY,J)-RAP*SOLH(LR,LEL,1)
	ENDDO
	ENDDO
	ENDDO
	ENDDO
	RETURN
	ENDIF
C 
C  100	FORMAT(1X,10E10.3)
C  101	FORMAT(1X,' ')
C  200	FORMAT(10X,'LT=',I3)
C  202	FORMAT(10X,'LT=',I3,'  N64=',I3)
C  203	FORMAT(10X,'LT=',I3,' N64=',I3,' N257=',I3,' LF1=',I3,' LF0=',I3)
C  204   FORMAT(1X,20I3)
	RETURN
	END
C