C $Header: /u/gcmpack/MITgcm/pkg/sphere/sphere.F,v 1.3 2007/10/09 00:11:39 jmc Exp $
C $Name: $
#include "CPP_OPTIONS.h"
c ==================================================================
c
c shpere.F: Routines that handle the projection onto sherical
c harmonics.
c
c Routines:
c
c o adfsc4dat - Adjoint routine of fsc4dat.
c o shc2grid - Evaluate a spherical harmonics model on a
c regular grid.
c o shc4grid - Evaluate a spherical harmonics model on a
c regular grid.
c
c o shcrotate - s/r used for the sph analysis ...
c o shc2zone - s/r used for the sph analysis ...
c o shc4zone - s/r used for the sph analysis ...
c o fsc2dat - s/r used for the sph analysis ...
c o frsbase - s/r used for the sph analysis ...
c o helmholtz - s/r used for the sph analysis ...
c o recur_z - s/r used for the sph analysis ...
c
c o shcError - Print error message and stop program (see below).
c
c IMSL Routines:
c
c o fftrb - Compute the real periodic sequence from its Fourier
c coefficients.
c --> replace by NAGLIB routine. C06...
c
c Platform-specific:
c
c M abort - FORTRAN intrinsic on SUN and, presumably, on CRAY to
c exit the program if an error condition is met.
c
c --> Replaced ABORT by shcError ( Christian Eckert MIT 22-Mar-2000 )
c
c Note:
c =====
c
c Where code is written in lower case letters, changes were
c introduced. The changes are mainly related to F90 syntax.
c Additionally, I replaced the call to the intrinsic *AMOD*
c by its generic name *MOD*. ( Ch.E. 25-May-1999 )
c
c
c Documentation:
c
c
c ==================================================================
SUBROUTINE FSC4DAT(N,FSC)
IMPLICIT NONE
INTEGER N
REAL FSC(N)
REAL SCALE
integer i
#ifdef USE_SPH_PROJECTION
CALL FFTRF(N,FSC,FSC)
#else
do i=1,n
fsc(i) = 0.0
enddo
#endif
SCALE = 2.0/FLOAT(N)
FSC(1) = FSC(1)/FLOAT(N)
CE FSC(2:N-1:2) = FSC(2:N-1:2)*SCALE
CE FSC(3:N :2) =-FSC(3:N: 2)*SCALE ! change sign of SINE coeffs.
do i = 2,n-1,2
fsc( i ) = fsc( i )*scale
enddo
do i = 3,n,2
fsc( i ) = -fsc( i )*scale ! change sign of sine coeffs.
enddo
IF (MOD(N,2).EQ.0) THEN
FSC(N) = FSC(N)/FLOAT(N)
ENDIF
RETURN
END
c ==================================================================
subroutine ADFSC4DAT(n,adfsc)
implicit none
integer n
real adfsc(n)
integer i
#ifdef USE_SPH_PROJECTION
call FFTRB(n,adfsc,adfsc)
#else
do i=1,n
adfsc(i) = 0.0
enddo
#endif
do i=1,n
adfsc(i) = adfsc(i)/float(n)
enddo
end
c ==================================================================
SUBROUTINE SHC2GRID( LMAX, SHC, NLON, NLAT, GRID, P )
C
C --- Evaluate a spherical harmonics model on a regular grid.
C
C SHC ! spherical harmonic coefficients in the order of
C !
C ! C00 : SHC(1)
C ! S11 C10 C11 : SHC(2) SHC(3) SHC(4)
C ! S22 S21 C20 C21 C22 : SHC(5) SHC(6) SHC(7) SHC(8) SHC(9)
C !
C ! i.e., C(L,M) = SHC(L*L+L+1+M)
C ! S(L,M) = SHC(L*L+L+1-M)
C NLAT ! number of latitude (zonal) lines.
C NLON ! number of longitude (meridional) lines.
C !
C GRID ! grid values in the order of
C ! ((west to east),south to north)
C ! ((LON=1,NLON),LAT=1,NLAT)
C ! grid values are defined on the "centers"
C ! of geographically regular equi-angular blocks.
C
C --- Coded by Dr. Myung-Chan Kim, Center for Space Research, 1997
C University of Texas at Austin, Austin, Texas, 78712
C
IMPLICIT NONE
INTEGER LMAX ! INPUT max. degree
REAL SHC((1+LMAX)*(1+LMAX)) ! INPUT spherical harmonics
INTEGER NLAT ! INPUT number of latitude lines.
INTEGER NLON ! INPUT number of longitude lines.
REAL GRID(NLON,NLAT) ! OUTPUT grid values
REAL P((LMAX+1)*(LMAX+2)/2) ! work space
INTEGER NLONLMT
PARAMETER (NLONLMT=10000)
REAL HS (NLONLMT)
REAL HN (NLONLMT)
REAL DLAT, ANGLE, XLAT1, XLAT2
INTEGER LATS, LATN
ce integer js, jn
integer i
ce IF(NLON.GT.NLONLMT) CALL ABORT('NLON.GT.NLONLMT')
ce IF(NLON.LT.LMAX*2 ) CALL ABORT('NLON.LT.LMAX*2 ')
IF(NLON.GT.NLONLMT) CALL SHCERROR('NLON.GT.NLONLMT',1)
IF(NLON.LT.LMAX*2 ) CALL SHCERROR('NLON.LT.LMAX*2 ',1)
DLAT = 180.0/FLOAT(NLAT)
ANGLE = 180.0/FLOAT(NLON)
CALL SHCROTATE(LMAX,SHC,ANGLE)
DO LATS = 1,(NLAT+1)/2
LATN = NLAT-(LATS-1)
XLAT2 =-90.0 + FLOAT(LATS)*DLAT
XLAT1 = XLAT2-DLAT
do i=1,nlon
hs(i) = 0.0
hn(i) = 0.0
enddo
CALL SHC2ZONE( LMAX, SHC, XLAT1, XLAT2, HS, HN, P )
CALL FSC2DAT( NLON, HS )
CALL FSC2DAT( NLON, HN )
do i=1,nlon
grid(i,lats) = hs(i)
grid(i,latn) = hn(i)
enddo
ENDDO
CALL SHCROTATE(LMAX,SHC,-ANGLE)
RETURN
END
c ==================================================================
SUBROUTINE SHC4GRID( LMAX, SHC, NLON, NLAT, GRID, P )
IMPLICIT NONE
INTEGER LMAX ! INPUT max. degree
REAL SHC((1+LMAX)*(1+LMAX)) ! OUTPUT spherical harmonics
INTEGER NLAT ! INPUT number of latitude lines.
INTEGER NLON ! INPUT number of longitude lines.
REAL GRID(NLON,NLAT) ! INPUT grid values
REAL P((LMAX+1)*(LMAX+2)/2) ! work space
INTEGER NLONLMT
PARAMETER (NLONLMT=10000)
REAL HS (NLONLMT)
REAL HN (NLONLMT)
REAL DLAT, ANGLE, XLAT1, XLAT2
INTEGER LATS, LATN
ce integer js, jn
integer i
IF(NLON.GT.NLONLMT) THEN
PRINT *, 'NLON = ', NLON
PRINT *, 'NLONLMT = ', NLONLMT
ce CALL ABORT('NLON.GT.NLONLMT')
CALL SHCERROR('NLON.GT.NLONLMT',1)
END
IF
IF(NLON.LT.LMAX*2 ) THEN
PRINT *, 'NLON = ', NLON
PRINT *, 'LMAX = ', LMAX
ce CALL ABORT('NLON.LT.LMAX*2')
CALL SHCERROR('NLON.LT.LMAX*2',1)
END
IF
DLAT = 180.0/FLOAT(NLAT)
ANGLE = 180.0/FLOAT(NLON)
DO LATS = 1,(NLAT+1)/2
LATN = NLAT-(LATS-1)
do i = 1,nlon
hs(i) = grid(i,lats)
hn(i) = grid(i,latn)
enddo
CALL FSC4DAT(NLON,HS)
CALL FSC4DAT(NLON,HN)
XLAT2 =-90.0 + FLOAT(LATS)*DLAT
XLAT1 = XLAT2-DLAT
CALL SHC4ZONE(LMAX,SHC,XLAT1,XLAT2,HS,HN,P)
ENDDO
CALL SHCROTATE(LMAX,SHC,-ANGLE)
RETURN
END
c ==================================================================
SUBROUTINE SHCROTATE(LMAX,SHC,ANGLE)
IMPLICIT NONE
INTEGER LMAX ! max. degree of spherical harmonics.
REAL SHC((1+LMAX)*(1+LMAX)) ! spherical harmonic coeffs.
REAL ANGLE ! in degree.
INTEGER LMAXLMT
PARAMETER (LMAXLMT=10000)
REAL H(LMAXLMT*2+1)
INTEGER K, L, M
REAL SINA, COSA, C, S
if (mod(angle,360.0) .ne. 0.0) then
CALL FRSBASE(ANGLE,H,1,1+LMAX+LMAX)
DO M = 0,LMAX
IF(M.EQ.0) THEN
SINA = 0.0
COSA = 1.0
ELSE
COSA = H(M+M)
SINA = H(M+M+1)
ENDIF
DO L = M,LMAX
K = L*L+L+1
C = SHC(K+M)
S = SHC(K-M)
SHC(K+M) = COSA*C + SINA*S
SHC(K-M) =-SINA*C + COSA*S
ENDDO
ENDDO
ENDIF
RETURN
END
c ==================================================================
SUBROUTINE SHC2ZONE(LMAX,SHC,XLAT1,XLAT2,HS,HN,P)
IMPLICIT NONE
INTEGER LMAX ! INPUT max. degree of spher. harmonics.
REAL SHC((1+LMAX)*(1+LMAX)) ! INPUT spherical harmonic coeffs.
REAL XLAT1 ! INPUT latitude.
REAL XLAT2 ! INPUT latitude.
REAL HS(1+LMAX+LMAX) ! OUTPUT
REAL HN(1+LMAX+LMAX) ! OUTPUT
REAL P((LMAX+1)*(LMAX+2)/2) ! INPUT
INTEGER LMAXLMT
PARAMETER (LMAXLMT=5000)
REAL FACT(0:LMAXLMT) !
INTEGER LMAX1, J, K, L, M
REAL DEG2RAD, SLAT
REAL A, B, E, F, Q, DA, DB, XLAT
ce IF (LMAX.GT.LMAXLMT) CALL ABORT('LMAX.GT.LMAXLMT')
cph IF (LMAX.GT.LMAXLMT) CALL shcError('LMAX.GT.LMAXLMT',1)
DEG2RAD = ACOS(-1.0)/180.0
XLAT = 0.5*(XLAT1+XLAT2)
do k = 1,lmax
fact(k) = 1.0
enddo
SLAT = SIN(XLAT*DEG2RAD)
CALL HELMHOLTZ(LMAX,SLAT,P)
LMAX1 = LMAX+1
K = 0
DO M = 0,LMAX
A = 0.0
B = 0.0
E = 0.0
F = 0.0
DO L = M,LMAX-1,2
K = K+1
J = L*L+L+1
Q = FACT(L)*P(K)
DA = Q*SHC(J+M)
DB = Q*SHC(J-M)
A = A+DA
B = B+DB
E = E+DA
F = F+DB
K = K+1
J = J+L+L+2
Q = FACT(L+1)*P(K)
DA = Q*SHC(J+M)
DB = Q*SHC(J-M)
A = A+DA
B = B+DB
E = E-DA
F = F-DB
ENDDO
IF(MOD(LMAX-L,2).EQ.0) THEN
K = K+1
J = LMAX*LMAX+LMAX+1
Q = FACT(LMAX)*P(K)
DA = Q*SHC(J+M)
DB = Q*SHC(J-M)
A = A+DA
B = B+DB
E = E+DA
F = F+DB
ENDIF
IF(M.EQ.0) THEN
HS(1) = A
HN(1) = E
ELSE
HS(M+M ) = A
HS(M+M+1) = B
HN(M+M ) = E
HN(M+M+1) = F
ENDIF
ENDDO
RETURN
END
c ==================================================================
SUBROUTINE SHC4ZONE(LMAX,SHC,XLAT1,XLAT2,HS,HN,P)
IMPLICIT NONE
INTEGER LMAX
REAL SHC((1+LMAX)*(1+LMAX)) ! spherical harmonic coeffs.
REAL XLAT1 ! latitude.
REAL XLAT2 ! latitude.
REAL P((LMAX+1)*(LMAX+2)/2)
C
C - output -
C
REAL HS(1+LMAX+LMAX)
REAL HN(1+LMAX+LMAX)
C
C - work space -
C
INTEGER LMAXLMT
PARAMETER (LMAXLMT=5000)
REAL XLAT, SIN0, SIN1, SIN2, SCALE
REAL DEG2RAD, CE, CO, SE, SO
INTEGER I, J, K, L, M
INTEGER JP, I1, I2
ce IF(LMAX.GT.LMAXLMT) CALL ABORT('LMAX.GT.LMAXLMT')
cph IF(LMAX.GT.LMAXLMT) CALL shcError('LMAX.GT.LMAXLMT',1)
IF(XLAT1 .LT. -90.0+1.E-10) THEN
do i = 1,(1+lmax)*(1+lmax)
shc(i) = 0.0
enddo
ENDIF
XLAT = 0.5*(XLAT1+XLAT2)
DEG2RAD = ACOS(-1.0)/180.0
SIN0 = SIN(XLAT *DEG2RAD)
SIN1 = SIN(AMAX1(-90.0,XLAT1)*DEG2RAD)
SIN2 = SIN(AMIN1( 0.0,XLAT2)*DEG2RAD)
SCALE = (SIN2 - SIN1)*0.25
do i = 1,lmax+lmax+1
hs(i) = hs(i)*scale
hn(i) = hn(i)*scale
enddo
CALL HELMHOLTZ(LMAX,SIN0,P)
I = 0
cadj loop = parallel
DO M = 0,LMAX
IF (M .EQ. 0) THEN
J = 1
JP = 1
ELSE
J = 2*M
JP = J+1
ENDIF
CE = HS(J)+HN(J)
CO = HS(J)-HN(J)
SE = HS(J)+HN(JP)
SO = HS(J)-HN(JP)
I1 = I+1
I2 = I+2
cadj loop = parallel
DO L = M,LMAX-1,2
K = L*L+L+1
SHC(K+M) = SHC(K+M) + P(I1) * CE
SHC(K-M) = SHC(K-M) + P(I1) * SE
K = K+L+L+2
SHC(K+M) = SHC(K+M) + P(I2) * CO
SHC(K-M) = SHC(K-M) + P(I2) * SO
ENDDO
I = I + 2
IF(MOD(LMAX-M,2).NE.1) THEN
I = I+1
K = LMAX*LMAX+LMAX+1
SHC(K+M) = SHC(K+M) + P(I) * CE
SHC(K-M) = SHC(K-M) + P(I) * SE
ENDIF
ENDDO
RETURN
END
c ==================================================================
subroutine FSC2DAT(n,fsc)
c
c --- this routine are coded to clarify the imsl's fftrf/fftrb.
c
implicit none
integer n
real fsc(n)
integer i
do i = 2,n-1,2
fsc(i ) = fsc(i )*0.5
fsc(i+1) = -fsc(i+1)*0.5 ! change sign of sine coeffs.
enddo
#ifdef USE_SPH_PROJECTION
call FFTRB(n,fsc,fsc)
#else
do i = 1,n
fsc( i ) = 0.0
enddo
#endif
return
end
c ==================================================================
SUBROUTINE FRSBASE(A,H,I,J)
C
C --- Given A (in degree), return H(I:J) = 1,COS(A),SIN(A).....
C
IMPLICIT NONE
REAL A
REAL H(1)
INTEGER I, J
REAL DEG2RAD, T, ARG
INTEGER N, L
DEG2RAD = ACOS(-1.0)/180.0
N = J-I+1
IF(N.LE.0) RETURN
H(I) = 1.0
IF(N.EQ.1) RETURN
ARG = A * DEG2RAD
H(I+1) = COS(ARG)
IF(N.EQ.2) RETURN
H(I+2) = SIN(ARG)
IF(N.EQ.3) RETURN
T = H(I+1)+H(I+1)
H(I+3) = T*H(I+1) - 1.0
IF(N.EQ.4) RETURN
H(I+4) = T*H(I+2)
IF(N.EQ.5) RETURN
DO L = I+5,J
H(L) = T*H(L-2)-H(L-4)
END
DO
RETURN
END
c ==================================================================
SUBROUTINE HELMHOLTZ(LMAX,S,P)
C
C --- compute the fully normalized associated legendre polynomials.
C
IMPLICIT NONE
INTEGER LMAX ! INPUT max. degree.
REAL S ! INPUT sin(latitude).
REAL P(1) ! OUTPUT assoc. legendre polynomials.
INTEGER LMAXLMT
PARAMETER(LMAXLMT=10800)
REAL A(0:LMAXLMT)
REAL ISECT(0:LMAXLMT)
REAL X(LMAXLMT)
REAL Y(LMAXLMT+LMAXLMT)
REAL Z(LMAXLMT)
INTEGER LMAXOLD
DATA LMAXOLD/-1/
SAVE LMAXOLD, ISECT, X, Y, Z
REAL C, CM, AK
INTEGER K, L, N, M
IF(LMAX.NE.LMAXOLD) THEN
ISECT(0) = 1
DO L = 1,LMAX
X(L) = SNGL(1.0D0/DSQRT(DBLE(FLOAT(4*L*L-1))))
Y(L) = SNGL(DSQRT(DBLE(FLOAT(L))))
Y(L+LMAX) = SNGL(DSQRT(DBLE(FLOAT(L+LMAX))))
ISECT(L) = L*LMAX-L*(L-3)/2+1
ENDDO
CALL RECUR_Z(Z,LMAX)
LMAXOLD = LMAX
ENDIF
C = SNGL(DSQRT(1.0D0-DBLE(S)*DBLE(S)))
CM = 1.0
P(1) = 1.0
DO M = 1,LMAX
K = ISECT(M)
CM = C * CM
P(K) = Z(M) * CM
ENDDO
N = 1
DO M = 0,LMAX-N
K = ISECT(M)+N
L = N+M
AK = X(L)*Y(L+M)*Y(N)
P(K) = S*P(K-1)/AK
A(M) = AK
ENDDO
DO N = 2,LMAX
DO M = 0,LMAX-N
K = ISECT(M)+N
L = N+M
AK = X(L)*Y(L+M)*Y(N)
P(K) =(S*P(K-1)-P(K-2)*A(M))/AK
A(M) = AK
ENDDO
ENDDO
RETURN
END
c ==================================================================
SUBROUTINE RECUR_Z(Z,LMAX)
C
C --- Coefficients for fully normalized sectorial Legendre polynomials.
C
IMPLICIT NONE
INTEGER LMAX
REAL Z(LMAX)
DOUBLE PRECISION ZZ
INTEGER L
ZZ = 2.0D0
DO L = 1,LMAX
ZZ = ZZ * DBLE(FLOAT(L+L+1))/DBLE(FLOAT(L+L))
Z(L) = SNGL(DSQRT(ZZ))
ENDDO
RETURN
END
c ==================================================================
subroutine SHCERROR( errstring, ierr )
c
c --- Print an error message and exit. Written in order to replace
c the machine specific routine ABORT.
c ( Christian Eckert MIT 22-Mar-2000 )
c
implicit none
integer ierr
character*(*) errstring
if ( ierr .ne. 0 ) then
print*
print*,' sphere: ',errstring
print*
stop ' ... program stopped.'
endif
return
end
c ==================================================================