subroutine kernel(block_size, lmax, nj, ni, Gkernel, Xlm) c This subroutine builds kernel relating to a point on a sphere c to spherical harmonic coefficients. c f(theta,phi) = sum_{l,m} c_l^m Y_l^m(theta,phi) c Using real spherical harmonics to avoid -l terms. cccccccccc c Input: c block_size = size of blocks in degrees c lmax = maximum angular degrees c nj = number of data (for declaring Gkernel) c ni = number of model coefficients (for declaring Gkernel) c Output: c Gkernel = kernel (also saved in a file) c Matrices requiring proper dimensions: c Xlm = Legendre and normalisation for given l,m at a given point cccccccccc cccccccccc c Note: all required functions are included in this file cccccccccc c July 28, 1999. c Last Modified: September 11, 1999. real block_size c block_size = size of the blocks in degrees integer lmax c maximum angular degrees. real Xlm(0:lmax, 0:lmax) c Legendre function and factors for a given l,m at a given point. integer nj, ni c nj = number of data points c ni = number of spherical coefficients real Gkernel(nj, ni) c Gkernel is the matrix containing the kernel. c Gkernel elements correpond to spherical harmonics at a given point. real Pi parameter (Pi = 3.14159265358979) c-- value of pi. real rlat, rlon c latitude and longitude as used in do loops. real theta, phi c theta = colatitude c phi = longitude real XP c XP is a function which calculates Xlm. integer lmax c maximum angular degrees. integer i, j c i = counter which keeps track of columns of kernel matrix. c j = counter which keeps track of rows of kernel matrix. character file_kernel*45 integer l, m c variables for do loops c========================PROGRAM============================ do i = 1, ni do j = 1, nj Gkernel(j,i) = 0.0 ! cleaning end do ! end i loop end do ! end j loop j = 1 ! initial value do rlat = block_size/2.0 - 90.0, 90.0, block_size theta = 0.0 ! cleaning up if (rlat.gt.0.0) then theta = (rlat-90)*Pi/180.0 elseif (rlat.lt.0.0) then theta = (90-rlat)*Pi/180.0 elseif (rlat.eq.0) then theta = Pi/2.0 end if ! if (rlat.eq.0.0)... ! Pi/180.0 terms are for converting from degrees to radians. do l = 0, lmax do m = 0, lmax Xlm(l,m) = 0.0 ! cleaning end do ! end m loop end do ! end l loop do l = 0, lmax do m = 0, l Xlm(l,m) = XP(l,m,theta) ! getting Xlm values c print*, 'l = ', l, ' m = ', m, ' and Xlm(l,m) = ', Xlm(l,m) end do ! end m loop end do ! end l loop do rlon = block_size/2.0, 360.0, block_size i = 1 ! initial value phi = 0.0 ! cleanin phi = rlon*Pi/180.0 do l = 0, lmax do m = 0, l if (i.gt.ni) then print*, '!!Error: Gkernel not finished and yet '// + 'i is greater than ni (sub_kernel.f)' print*, 'Quitting...' go to 15 end if ! if (i.gt.ni) if (m.eq.0) then ! zonal term Gkernel(j,i) = Xlm(l,m) i = i + 1 elseif (m.ne.0) then ! non-zonal terms Gkernel(j,i) = ((-1.0)**m)*Sqrt(2.0)* + Xlm(l,m)*Cos(Float(m*phi)) ! cosine part Gkernel(j,i+1) = ((-1.0)**m)*Sqrt(2.0)* + Xlm(l,m)*Sin(Float(m*phi)) ! sine part i = i + 2 end if ! if (m.ne.0)... end do ! end m loop end do ! end l loop if (i-1.ne.ni) then print*, '!!Error: i-1 = ', i-1, ' but ni = ', ni, + ' (sub_kernel.f)' print*, 'Quitting sub_kernel.f ...' go to 15 end if ! if (i-1.ne.ni) if (j.gt.nj) then print*, '!!Error: Gkernel not finished and yet j is '// + 'greater than nj (sub_kernel.f)' print*, 'Quitting...' go to 15 end if ! if (j.gt.nj) j = j + 1 end do ! end rlon loop end do ! end rlat loop Finished making Gkernel c checking dimensions if (i-1.ne.ni) then print*, '!!Error: i-1 = ', i-1, ' but ni = ', ni, ' (sub_kernel.f)' print*, 'Quitting sub_kernel.f ...' go to 15 end if ! if (i-1.ne.ni) if (j-1.ne.nj) then print*, '!!Error: j-1 = ', j-1, ' but nj = ', nj, ' (sub_kernel.f)' print*, 'Quitting sub_kernel.f ...' go to 15 end if ! if (j-1.ne.nj) file_kernel = 'kernel.dat' print*, 'Saving kernel in file: ', file_kernel(1:index(file_kernel,' ')), + '(sub_kernel.f)' open(120, file = file_kernel) rewind(120) write(120,*), block_size do j = 1, nj write(120,*), (Gkernel(j,i), i = 1, ni) ! saving kernel end do ! end j loop close(120) 15 return END cc***************************************************** real function XP(l,m,theta) c This function calculates the colatitude dependent part of the c spherical harmonics. integer l, m c l = spherical harmonics angular degree c m = spherical harmonics angular order real theta c colatitude real Pi parameter (Pi = 3.14159265358979) real factrl, plgndr XP = ((-1.0)**m)*Sqrt((2.0*float(l)+1.0)/(4.0*Pi))* + Sqrt(FACTRL(l-m)/FACTRL(l+m))* + PLGNDR(l,m,Cos(theta)) return end c******************************************* real FUNCTION FACTRL(N) c from the Numerical Recipes... integer N real A(33) integer NTOP DATA NTOP,A(1)/0,1./ integer j real gammln IF (N.LT.0) THEN PAUSE 'negative factorial' ELSE IF (N.LE.NTOP) THEN FACTRL=A(N+1) ELSE IF (N.LE.32) THEN DO 11 J=NTOP+1,N A(J+1)=J*A(J) 11 CONTINUE NTOP=N FACTRL=A(N+1) ELSE FACTRL=EXP(GAMMLN(N+1.)) ENDIF RETURN END c******************************************* real FUNCTION GAMMLN(XX) c from the Numerical Recipes... REAL*8 COF(6),STP,HALF,ONE,FPF,X,TMP,SER DATA COF,STP/76.18009173D0,-86.50532033D0,24.01409822D0, * -1.231739516D0,.120858003D-2,-.536382D-5,2.50662827465D0/ DATA HALF,ONE,FPF/0.5D0,1.0D0,5.5D0/ real xx real*8 x, tmp integer j X=Dble(XX)-ONE TMP=X+FPF TMP=(X+HALF)*LOG(TMP)-TMP SER=ONE DO 11 J=1,6 X=X+ONE SER=SER+COF(J)/X 11 CONTINUE GAMMLN=Sngl(TMP+LOG(STP*SER)) RETURN END c******************************************* real FUNCTION PLGNDR(L,M,X) c This is from Numerical Recipes... integer l, m real x real pmm, somx2, fact, pmmp1, pll integer i, ll IF(M.LT.0.OR.M.GT.L.OR.ABS(X).GT.1.)PAUSE 'bad arguments' PMM=1. IF(M.GT.0) THEN SOMX2=SQRT((1.-X)*(1.+X)) FACT=1. DO 11 I=1,M PMM=-PMM*FACT*SOMX2 FACT=FACT+2. 11 CONTINUE ENDIF IF(L.EQ.M) THEN PLGNDR=PMM ELSE PMMP1=X*(2*M+1)*PMM IF(L.EQ.M+1) THEN PLGNDR=PMMP1 ELSE DO 12 LL=M+2,L PLL=(X*(2*LL-1)*PMMP1-(LL+M-1)*PMM)/(LL-M) PMM=PMMP1 PMMP1=PLL 12 CONTINUE PLGNDR=PLL ENDIF ENDIF RETURN END