pkg/sphere/sphere.F

0bed5b371d Patr*0001 #include "CPP_OPTIONS.h"
                0002
                0003 c     ==================================================================
                0004 c
                0005 c     shpere.F: Routines that handle the projection onto sherical
                0006 c               harmonics.
                0007 c
                0008 c     Routines:
                0009 c
                0010 c     o adfsc4dat  - Adjoint routine of fsc4dat.
                0011 c     o shc2grid   - Evaluate a spherical harmonics model on a
                0012 c                    regular grid.
                0013 c     o shc4grid   - Evaluate a spherical harmonics model on a
                0014 c                    regular grid.
                0015 c
                0016 c     o shcrotate  - s/r used for the sph analysis ...
                0017 c     o shc2zone   - s/r used for the sph analysis ...
                0018 c     o shc4zone   - s/r used for the sph analysis ...
                0019 c     o fsc2dat    - s/r used for the sph analysis ...
                0020 c     o frsbase    - s/r used for the sph analysis ...
                0021 c     o helmholtz  - s/r used for the sph analysis ...
                0022 c     o recur_z    - s/r used for the sph analysis ...
                0023 c
                0024 c     o shcError   - Print error message and stop program (see below).
                0025 c
                0026 c     IMSL Routines:
                0027 c
                0028 c     o fftrb      - Compute the real periodic sequence from its Fourier
                0029 c                    coefficients.
                0030 c                    --> replace by NAGLIB routine. C06...
                0031 c
                0032 c     Platform-specific:
                0033 c
                0034 c     M abort      - FORTRAN intrinsic on SUN and, presumably, on CRAY to
                0035 c                    exit the program if an error condition is met.
                0036 c
                0037 c     -->  Replaced ABORT by shcError ( Christian Eckert  MIT  22-Mar-2000 )
                0038 c
                0039 c     Note:
                0040 c     =====
                0041 c
                0042 c     Where code is written in lower case letters, changes were
                0043 c     introduced. The changes are mainly related to F90 syntax.
                0044 c     Additionally, I replaced the call to the intrinsic *AMOD*
                0045 c     by its generic name *MOD*. ( Ch.E. 25-May-1999 )
                0046 c
                0047 c
                0048 c     Documentation:
                0049 c
                0050 c
                0051 c     ==================================================================
                0052
                0053       SUBROUTINE FSC4DAT(N,FSC)
                0054
                0055       IMPLICIT NONE
                0056
                0057       INTEGER N
                0058       REAL FSC(N)
                0059       REAL SCALE
                0060       integer i
                0061
                0062 #ifdef USE_SPH_PROJECTION
                0063       CALL FFTRF(N,FSC,FSC)
                0064 #else
                0065       do i=1,n
                0066         fsc(i) = 0.0
                0067       enddo
                0068 #endif
                0069
                0070       SCALE        = 2.0/FLOAT(N)
                0071       FSC(1)       = FSC(1)/FLOAT(N)
                0072
                0073 CE      FSC(2:N-1:2) = FSC(2:N-1:2)*SCALE
                0074 CE      FSC(3:N  :2) =-FSC(3:N:  2)*SCALE ! change sign of SINE coeffs.
                0075
                0076       do i = 2,n-1,2
                0077         fsc( i ) =  fsc( i )*scale
                0078       enddo
                0079       do i = 3,n,2
                0080         fsc( i ) = -fsc( i )*scale ! change sign of sine coeffs.
                0081       enddo
                0082
                0083       IF (MOD(N,2).EQ.0) THEN
                0084          FSC(N) = FSC(N)/FLOAT(N)
                0085       ENDIF
                0086
                0087       RETURN
                0088       END
                0089
                0090 c     ==================================================================
                0091
                0092       subroutine adfsc4dat(n,adfsc)
                0093
                0094       implicit none
                0095
                0096       integer n
                0097       real adfsc(n)
                0098
                0099       integer i
                0100
                0101 #ifdef USE_SPH_PROJECTION
                0102       call fftrb(n,adfsc,adfsc)
                0103 #else
                0104       do i=1,n
                0105         adfsc(i) = 0.0
                0106       enddo
                0107 #endif
                0108
                0109       do i=1,n
                0110         adfsc(i) = adfsc(i)/float(n)
                0111       enddo
                0112
                0113       end
                0114
                0115 c     ==================================================================
                0116
                0117       SUBROUTINE SHC2GRID( LMAX, SHC, NLON, NLAT, GRID, P )
                0118 C
                0119 C --- Evaluate a spherical harmonics model on a regular grid.
                0120 C
                0121 C      SHC  ! spherical harmonic coefficients in the order of
                0122 C           !
                0123 C           !         C00         :               SHC(1)
                0124 C           !     S11 C10 C11     :        SHC(2) SHC(3) SHC(4)
                0125 C           ! S22 S21 C20 C21 C22 : SHC(5) SHC(6) SHC(7) SHC(8) SHC(9)
                0126 C           !
                0127 C           ! i.e.,  C(L,M) = SHC(L*L+L+1+M)
                0128 C           !        S(L,M) = SHC(L*L+L+1-M)
                0129 C     NLAT  ! number of latitude (zonal) lines.
                0130 C     NLON  ! number of longitude (meridional) lines.
                0131 C           !
                0132 C     GRID  ! grid values in the order of
                0133 C           ! ((west to east),south to north)
                0134 C           ! ((LON=1,NLON),LAT=1,NLAT)
                0135 C           ! grid values are defined on the "centers"
                0136 C           ! of geographically regular equi-angular blocks.
                0137 C
                0138 C --- Coded by Dr. Myung-Chan Kim, Center for Space Research, 1997
                0139 C              University of Texas at Austin, Austin, Texas, 78712
                0140 C
                0141       IMPLICIT NONE
                0142
                0143       INTEGER LMAX                       ! INPUT  max. degree
                0144       REAL    SHC((1+LMAX)*(1+LMAX))     ! INPUT  spherical harmonics
                0145       INTEGER NLAT                       ! INPUT  number of latitude  lines.
                0146       INTEGER NLON                       ! INPUT  number of longitude lines.
                0147       REAL    GRID(NLON,NLAT)            ! OUTPUT grid values
                0148       REAL    P((LMAX+1)*(LMAX+2)/2)     ! work space
                0149
                0150       INTEGER NLONLMT
                0151       PARAMETER (NLONLMT=10000)
                0152       REAL HS (NLONLMT)
                0153       REAL HN (NLONLMT)
                0154
                0155       REAL    DLAT, ANGLE, XLAT1, XLAT2
                0156       INTEGER LATS, LATN
                0157
                0158 ce    integer js, jn
                0159
                0160       integer i
                0161
                0162 ce      IF(NLON.GT.NLONLMT) CALL ABORT('NLON.GT.NLONLMT')
                0163 ce      IF(NLON.LT.LMAX*2 ) CALL ABORT('NLON.LT.LMAX*2 ')
                0164       IF(NLON.GT.NLONLMT) CALL shcError('NLON.GT.NLONLMT',1)
                0165       IF(NLON.LT.LMAX*2 ) CALL shcError('NLON.LT.LMAX*2 ',1)
                0166
                0167       DLAT  = 180.0/FLOAT(NLAT)
                0168       ANGLE = 180.0/FLOAT(NLON)
                0169
                0170       CALL SHCROTATE(LMAX,SHC,ANGLE)
                0171       DO LATS = 1,(NLAT+1)/2
                0172          LATN = NLAT-(LATS-1)
                0173          XLAT2 =-90.0 + FLOAT(LATS)*DLAT
                0174          XLAT1 = XLAT2-DLAT
                0175          do i=1,nlon
                0176            hs(i) = 0.0
                0177            hn(i) = 0.0
                0178          enddo
                0179          CALL SHC2ZONE( LMAX, SHC, XLAT1, XLAT2, HS, HN, P )
                0180          CALL FSC2DAT( NLON, HS )
                0181          CALL FSC2DAT( NLON, HN )
                0182          do i=1,nlon
                0183            grid(i,lats) = hs(i)
                0184            grid(i,latn) = hn(i)
                0185          enddo
                0186       ENDDO
                0187       CALL SHCROTATE(LMAX,SHC,-ANGLE)
                0188       RETURN
                0189       END
                0190
                0191 c     ==================================================================
                0192
                0193       SUBROUTINE SHC4GRID( LMAX, SHC, NLON, NLAT, GRID, P )
                0194
                0195       IMPLICIT NONE
                0196
                0197       INTEGER LMAX                       ! INPUT  max. degree
                0198       REAL    SHC((1+LMAX)*(1+LMAX))     ! OUTPUT spherical harmonics
                0199       INTEGER NLAT                       ! INPUT  number of latitude  lines.
                0200       INTEGER NLON                       ! INPUT  number of longitude lines.
                0201       REAL    GRID(NLON,NLAT)            ! INPUT  grid values
                0202       REAL    P((LMAX+1)*(LMAX+2)/2)     ! work space
                0203
                0204       INTEGER NLONLMT
                0205       PARAMETER (NLONLMT=10000)
                0206       REAL HS (NLONLMT)
                0207       REAL HN (NLONLMT)
                0208
                0209       REAL    DLAT, ANGLE, XLAT1, XLAT2
                0210       INTEGER LATS, LATN
                0211
                0212 ce    integer js, jn
                0213
                0214       integer i
                0215
                0216       IF(NLON.GT.NLONLMT) THEN
                0217          PRINT *, 'NLON = ', NLON
                0218          PRINT *, 'NLONLMT = ', NLONLMT
                0219 ce         CALL ABORT('NLON.GT.NLONLMT')
                0220          CALL shcError('NLON.GT.NLONLMT',1)
                0221       END IF
                0222
                0223       IF(NLON.LT.LMAX*2 ) THEN
                0224          PRINT *, 'NLON = ', NLON
                0225          PRINT *, 'LMAX = ', LMAX
                0226 ce         CALL ABORT('NLON.LT.LMAX*2')
                0227          CALL shcError('NLON.LT.LMAX*2',1)
                0228       END IF
                0229
                0230       DLAT  = 180.0/FLOAT(NLAT)
                0231       ANGLE = 180.0/FLOAT(NLON)
                0232       DO LATS = 1,(NLAT+1)/2
                0233          LATN = NLAT-(LATS-1)
                0234          do i = 1,nlon
                0235            hs(i) = grid(i,lats)
                0236            hn(i) = grid(i,latn)
                0237          enddo
                0238          CALL FSC4DAT(NLON,HS)
                0239          CALL FSC4DAT(NLON,HN)
                0240          XLAT2 =-90.0 + FLOAT(LATS)*DLAT
                0241          XLAT1 = XLAT2-DLAT
                0242          CALL SHC4ZONE(LMAX,SHC,XLAT1,XLAT2,HS,HN,P)
                0243       ENDDO
                0244       CALL SHCROTATE(LMAX,SHC,-ANGLE)
                0245
                0246       RETURN
                0247       END
                0248
                0249 c     ==================================================================
                0250
                0251       SUBROUTINE SHCROTATE(LMAX,SHC,ANGLE)
                0252
                0253       IMPLICIT NONE
                0254
                0255       INTEGER  LMAX                   ! max. degree of spherical harmonics.
                0256       REAL     SHC((1+LMAX)*(1+LMAX)) ! spherical harmonic coeffs.
                0257       REAL     ANGLE                  ! in degree.
                0258
                0259       INTEGER LMAXLMT
                0260       PARAMETER (LMAXLMT=10000)
                0261
                0262       REAL    H(LMAXLMT*2+1)
                0263       INTEGER K, L, M
                0264       REAL    SINA, COSA, C, S
                0265
                0266       if (mod(angle,360.0) .ne. 0.0) then
                0267
                0268          CALL FRSBASE(ANGLE,H,1,1+LMAX+LMAX)
                0269
                0270          DO M = 0,LMAX
                0271             IF(M.EQ.0) THEN
                0272                SINA = 0.0
                0273                COSA = 1.0
                0274             ELSE
                0275                COSA = H(M+M)
                0276                SINA = H(M+M+1)
                0277             ENDIF
                0278             DO L = M,LMAX
                0279                K = L*L+L+1
                0280                C = SHC(K+M)
                0281                S = SHC(K-M)
                0282                SHC(K+M) = COSA*C + SINA*S
                0283                SHC(K-M) =-SINA*C + COSA*S
                0284             ENDDO
                0285          ENDDO
                0286
                0287       ENDIF
                0288
                0289       RETURN
                0290       END
                0291
                0292 c     ==================================================================
                0293
                0294       SUBROUTINE SHC2ZONE(LMAX,SHC,XLAT1,XLAT2,HS,HN,P)
                0295
                0296       IMPLICIT NONE
                0297
                0298       INTEGER LMAX                   ! INPUT  max. degree of spher. harmonics.
                0299       REAL    SHC((1+LMAX)*(1+LMAX)) ! INPUT  spherical harmonic coeffs.
                0300       REAL    XLAT1                  ! INPUT  latitude.
                0301       REAL    XLAT2                  ! INPUT  latitude.
                0302       REAL    HS(1+LMAX+LMAX)        ! OUTPUT
                0303       REAL    HN(1+LMAX+LMAX)        ! OUTPUT
                0304       REAL    P((LMAX+1)*(LMAX+2)/2) ! INPUT
                0305
                0306       INTEGER LMAXLMT
                0307       PARAMETER (LMAXLMT=5000)
                0308
                0309       REAL    FACT(0:LMAXLMT) !
                0310       INTEGER LMAX1, J, K, L, M
                0311       REAL    DEG2RAD, SLAT
                0312       REAL    A, B, E, F, Q, DA, DB, XLAT
                0313
                0314 ce      IF (LMAX.GT.LMAXLMT) CALL ABORT('LMAX.GT.LMAXLMT')
538ab314d4 Patr*0315 cph      IF (LMAX.GT.LMAXLMT) CALL shcError('LMAX.GT.LMAXLMT',1)
0bed5b371d Patr*0316
                0317       DEG2RAD = ACOS(-1.0)/180.0
                0318       XLAT = 0.5*(XLAT1+XLAT2)
                0319       do k = 1,lmax
                0320         fact(k) = 1.0
                0321       enddo
                0322       SLAT = SIN(XLAT*DEG2RAD)
                0323
                0324       CALL HELMHOLTZ(LMAX,SLAT,P)
                0325       LMAX1 = LMAX+1
                0326       K    = 0
                0327       DO M = 0,LMAX
                0328          A    = 0.0
                0329          B    = 0.0
                0330          E    = 0.0
                0331          F    = 0.0
                0332          DO L = M,LMAX-1,2
                0333             K    = K+1
                0334             J    = L*L+L+1
                0335             Q    = FACT(L)*P(K)
                0336             DA   = Q*SHC(J+M)
                0337             DB   = Q*SHC(J-M)
                0338             A    = A+DA
                0339             B    = B+DB
                0340             E    = E+DA
                0341             F    = F+DB
                0342             K    = K+1
                0343             J    = J+L+L+2
                0344             Q    = FACT(L+1)*P(K)
                0345             DA   = Q*SHC(J+M)
                0346             DB   = Q*SHC(J-M)
                0347             A    = A+DA
                0348             B    = B+DB
                0349             E    = E-DA
                0350             F    = F-DB
                0351          ENDDO
                0352          IF(MOD(LMAX-L,2).EQ.0) THEN
                0353             K    = K+1
                0354             J    = LMAX*LMAX+LMAX+1
                0355             Q    = FACT(LMAX)*P(K)
                0356             DA   = Q*SHC(J+M)
                0357             DB   = Q*SHC(J-M)
                0358             A    = A+DA
                0359             B    = B+DB
                0360             E    = E+DA
                0361             F    = F+DB
                0362          ENDIF
                0363          IF(M.EQ.0) THEN
                0364             HS(1) = A
                0365             HN(1) = E
                0366          ELSE
                0367             HS(M+M  ) = A
                0368             HS(M+M+1) = B
                0369             HN(M+M  ) = E
                0370             HN(M+M+1) = F
                0371          ENDIF
                0372       ENDDO
                0373
                0374       RETURN
                0375       END
                0376
                0377 c     ==================================================================
                0378
                0379       SUBROUTINE SHC4ZONE(LMAX,SHC,XLAT1,XLAT2,HS,HN,P)
                0380
                0381       IMPLICIT NONE
                0382
                0383       INTEGER LMAX
                0384       REAL    SHC((1+LMAX)*(1+LMAX))    ! spherical harmonic coeffs.
                0385       REAL    XLAT1                     ! latitude.
                0386       REAL    XLAT2                     ! latitude.
                0387       REAL    P((LMAX+1)*(LMAX+2)/2)
                0388 C
                0389 C   - output -
                0390 C
                0391       REAL HS(1+LMAX+LMAX)
                0392       REAL HN(1+LMAX+LMAX)
                0393 C
                0394 C   - work space -
                0395 C
                0396       INTEGER LMAXLMT
                0397       PARAMETER (LMAXLMT=5000)
                0398
                0399       REAL    XLAT, SIN0, SIN1, SIN2, SCALE
                0400       REAL    DEG2RAD, CE, CO, SE, SO
                0401       INTEGER I, J, K, L, M
                0402       INTEGER JP, I1, I2
                0403
                0404 ce      IF(LMAX.GT.LMAXLMT) CALL ABORT('LMAX.GT.LMAXLMT')
538ab314d4 Patr*0405 cph      IF(LMAX.GT.LMAXLMT) CALL shcError('LMAX.GT.LMAXLMT',1)
0bed5b371d Patr*0406
                0407       IF(XLAT1 .LT. -90.0+1.E-10) THEN
                0408         do i = 1,(1+lmax)*(1+lmax)
                0409           shc(i) = 0.0
                0410         enddo
                0411       ENDIF
                0412
                0413       XLAT  = 0.5*(XLAT1+XLAT2)
                0414       DEG2RAD  = ACOS(-1.0)/180.0
                0415          SIN0  = SIN(XLAT *DEG2RAD)
                0416          SIN1  = SIN(AMAX1(-90.0,XLAT1)*DEG2RAD)
                0417          SIN2  = SIN(AMIN1(  0.0,XLAT2)*DEG2RAD)
                0418          SCALE = (SIN2 - SIN1)*0.25
                0419       do i = 1,lmax+lmax+1
                0420         hs(i) = hs(i)*scale
                0421         hn(i) = hn(i)*scale
                0422       enddo
                0423
                0424       CALL HELMHOLTZ(LMAX,SIN0,P)
                0425       I    = 0
                0426 cadj loop = parallel
                0427       DO M = 0,LMAX
                0428          IF (M .EQ. 0) THEN
                0429             J  = 1
                0430             JP = 1
                0431          ELSE
                0432             J  = 2*M
                0433             JP = J+1
                0434          ENDIF
                0435          CE = HS(J)+HN(J)
                0436          CO = HS(J)-HN(J)
                0437          SE = HS(J)+HN(JP)
                0438          SO = HS(J)-HN(JP)
                0439          I1 = I+1
                0440          I2 = I+2
                0441 cadj loop = parallel
                0442          DO L = M,LMAX-1,2
                0443             K = L*L+L+1
                0444             SHC(K+M) = SHC(K+M) + P(I1) * CE
                0445             SHC(K-M) = SHC(K-M) + P(I1) * SE
                0446             K = K+L+L+2
7c7521a1da Jean*0447             SHC(K+M) = SHC(K+M) + P(I2) * CO
0bed5b371d Patr*0448             SHC(K-M) = SHC(K-M) + P(I2) * SO
                0449          ENDDO
                0450          I = I + 2
                0451          IF(MOD(LMAX-M,2).NE.1) THEN
                0452             I = I+1
                0453             K = LMAX*LMAX+LMAX+1
                0454             SHC(K+M) = SHC(K+M) + P(I) * CE
                0455             SHC(K-M) = SHC(K-M) + P(I) * SE
                0456          ENDIF
                0457       ENDDO
                0458       RETURN
                0459       END
                0460
                0461 c     ==================================================================
                0462
                0463       subroutine fsc2dat(n,fsc)
                0464 c
                0465 c --- this routine are coded to clarify the imsl's fftrf/fftrb.
                0466 c
                0467       implicit none
                0468
                0469       integer n
                0470       real fsc(n)
                0471
                0472       integer i
                0473
                0474       do i = 2,n-1,2
                0475         fsc(i  ) =  fsc(i  )*0.5
                0476         fsc(i+1) = -fsc(i+1)*0.5   ! change sign of sine coeffs.
                0477       enddo
                0478
                0479 #ifdef USE_SPH_PROJECTION
                0480       call fftrb(n,fsc,fsc)
                0481 #else
                0482       do i = 1,n
                0483         fsc( i ) = 0.0
                0484       enddo
                0485 #endif
                0486
                0487       return
                0488       end
                0489
                0490 c     ==================================================================
                0491
                0492       SUBROUTINE FRSBASE(A,H,I,J)
                0493 C
                0494 C --- Given A (in degree), return H(I:J) = 1,COS(A),SIN(A).....
                0495 C
                0496       IMPLICIT NONE
                0497
                0498       REAL A
                0499       REAL H(1)
                0500       INTEGER I, J
                0501
                0502       REAL DEG2RAD, T, ARG
                0503       INTEGER N, L
                0504
                0505       DEG2RAD = ACOS(-1.0)/180.0
                0506
                0507       N      = J-I+1
                0508       IF(N.LE.0) RETURN
                0509       H(I)   = 1.0
                0510       IF(N.EQ.1) RETURN
                0511       ARG    = A * DEG2RAD
                0512       H(I+1) = COS(ARG)
                0513       IF(N.EQ.2) RETURN
                0514       H(I+2) = SIN(ARG)
                0515       IF(N.EQ.3) RETURN
                0516       T      = H(I+1)+H(I+1)
                0517       H(I+3) = T*H(I+1) - 1.0
                0518       IF(N.EQ.4) RETURN
                0519       H(I+4) = T*H(I+2)
                0520       IF(N.EQ.5) RETURN
                0521       DO L = I+5,J
                0522       H(L)   = T*H(L-2)-H(L-4)
                0523       END DO
                0524
                0525       RETURN
7c7521a1da Jean*0526       END
0bed5b371d Patr*0527
                0528 c     ==================================================================
                0529
                0530       SUBROUTINE HELMHOLTZ(LMAX,S,P)
                0531 C
                0532 C --- compute the fully normalized associated legendre polynomials.
                0533 C
                0534       IMPLICIT NONE
                0535
                0536       INTEGER LMAX        ! INPUT  max. degree.
                0537       REAL    S           ! INPUT  sin(latitude).
                0538       REAL    P(1)        ! OUTPUT assoc. legendre polynomials.
                0539
                0540       INTEGER LMAXLMT
                0541       PARAMETER(LMAXLMT=10800)
                0542
                0543       REAL A(0:LMAXLMT)
                0544       REAL ISECT(0:LMAXLMT)
                0545       REAL X(LMAXLMT)
                0546       REAL Y(LMAXLMT+LMAXLMT)
                0547       REAL Z(LMAXLMT)
                0548
                0549       INTEGER LMAXOLD
                0550       DATA    LMAXOLD/-1/
                0551       SAVE    LMAXOLD, ISECT, X, Y, Z
                0552
                0553       REAL    C, CM, AK
                0554       INTEGER K, L, N, M
                0555
                0556       IF(LMAX.NE.LMAXOLD) THEN
                0557          ISECT(0) = 1
                0558          DO L = 1,LMAX
                0559             X(L) = SNGL(1.0D0/DSQRT(DBLE(FLOAT(4*L*L-1))))
                0560             Y(L) = SNGL(DSQRT(DBLE(FLOAT(L))))
                0561             Y(L+LMAX) = SNGL(DSQRT(DBLE(FLOAT(L+LMAX))))
                0562             ISECT(L) = L*LMAX-L*(L-3)/2+1
                0563          ENDDO
                0564          CALL RECUR_Z(Z,LMAX)
                0565          LMAXOLD = LMAX
                0566       ENDIF
                0567
                0568       C  = SNGL(DSQRT(1.0D0-DBLE(S)*DBLE(S)))
                0569       CM = 1.0
                0570          P(1) = 1.0
                0571          DO M = 1,LMAX
                0572             K = ISECT(M)
                0573             CM   =    C * CM
                0574             P(K) = Z(M) * CM
                0575          ENDDO
                0576          N    = 1
                0577          DO M = 0,LMAX-N
                0578             K = ISECT(M)+N
                0579             L = N+M
                0580             AK   = X(L)*Y(L+M)*Y(N)
                0581             P(K) = S*P(K-1)/AK
                0582             A(M) = AK
                0583          ENDDO
                0584          DO N = 2,LMAX
                0585             DO M = 0,LMAX-N
                0586                K = ISECT(M)+N
                0587                L = N+M
                0588                AK   = X(L)*Y(L+M)*Y(N)
                0589                P(K) =(S*P(K-1)-P(K-2)*A(M))/AK
                0590                A(M) = AK
                0591             ENDDO
                0592          ENDDO
                0593       RETURN
                0594       END
                0595
                0596 c     ==================================================================
                0597
                0598       SUBROUTINE RECUR_Z(Z,LMAX)
                0599 C
                0600 C --- Coefficients for fully normalized sectorial Legendre polynomials.
                0601 C
                0602       IMPLICIT NONE
                0603
                0604       INTEGER LMAX
                0605       REAL    Z(LMAX)
                0606
                0607       DOUBLE PRECISION ZZ
                0608       INTEGER L
                0609
                0610       ZZ   = 2.0D0
                0611       DO L = 1,LMAX
                0612          ZZ   = ZZ * DBLE(FLOAT(L+L+1))/DBLE(FLOAT(L+L))
                0613          Z(L) = SNGL(DSQRT(ZZ))
                0614       ENDDO
                0615
                0616       RETURN
                0617       END
                0618
                0619 c     ==================================================================
                0620
                0621       subroutine shcError( errstring, ierr )
                0622 c
                0623 c --- Print an error message and exit. Written in order to replace
                0624 c     the machine specific routine ABORT.
                0625 c     ( Christian Eckert  MIT  22-Mar-2000 )
                0626 c
                0627       implicit none
                0628
                0629       integer       ierr
                0630       character*(*) errstring
                0631
                0632       if ( ierr .ne. 0 ) then
                0633         print*
                0634         print*,' sphere: ',errstring
                0635         print*
                0636         stop   ' ... program stopped.'
                0637       endif
                0638
                0639       return
                0640       end
                0641
                0642 c     ==================================================================