Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F120681396
sphere.F
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sun, Jul 6, 06:31
Size
16 KB
Mime Type
text/x-c
Expires
Tue, Jul 8, 06:31 (2 d)
Engine
blob
Format
Raw Data
Handle
27215086
Attached To
R2795 mitgcm_lac_leman_abirani
sphere.F
View Options
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 ==================================================================
Event Timeline
Log In to Comment