32 subroutine splat(IDRT,JMAX,SLAT,WLAT)
37 integer,
intent(in ) :: idrt
38 integer,
intent(in ) :: jmax
39 real(kind=
kind_real),
intent(out) :: slat(jmax)
40 real(kind=
kind_real),
intent(out) :: wlat(jmax)
43 real(kind=
kind_real):: pk(jmax/2),pkm1(jmax/2),pkm2(jmax/2)
44 real(kind=
kind_real):: slatd(jmax/2),sp,spmax,eps
45 integer,
parameter:: jz=50
48 real(kind=
kind_real):: awork((jmax+1)/2,((jmax+1)/2))
49 real(kind=
kind_real):: bwork(((jmax+1)/2))
50 integer:: jh,jhe,jho,j0
51 integer:: ipvt((jmax+1)/2)
55 data bz / 2.4048255577_kind_real, 5.5200781103_kind_real, &
56 & 8.6537279129_kind_real, 11.7915344391_kind_real, &
57 & 14.9309177086_kind_real, 18.0710639679_kind_real, &
58 & 21.2116366299_kind_real, 24.3524715308_kind_real, &
59 & 27.4934791320_kind_real, 30.6346064684_kind_real, &
60 & 33.7758202136_kind_real, 36.9170983537_kind_real, &
61 & 40.0584257646_kind_real, 43.1997917132_kind_real, &
62 & 46.3411883717_kind_real, 49.4826098974_kind_real, &
63 & 52.6240518411_kind_real, 55.7655107550_kind_real, &
64 & 58.9069839261_kind_real, 62.0484691902_kind_real, &
65 & 65.1899648002_kind_real, 68.3314693299_kind_real, &
66 & 71.4729816036_kind_real, 74.6145006437_kind_real, &
67 & 77.7560256304_kind_real, 80.8975558711_kind_real, &
68 & 84.0390907769_kind_real, 87.1806298436_kind_real, &
69 & 90.3221726372_kind_real, 93.4637187819_kind_real, &
70 & 96.6052679510_kind_real, 99.7468198587_kind_real, &
71 & 102.888374254_kind_real, 106.029930916_kind_real, &
72 & 109.171489649_kind_real, 112.313050280_kind_real, &
73 & 115.454612653_kind_real, 118.596176630_kind_real, &
74 & 121.737742088_kind_real, 124.879308913_kind_real, &
75 & 128.020877005_kind_real, 131.162446275_kind_real, &
76 & 134.304016638_kind_real, 137.445588020_kind_real, &
77 & 140.587160352_kind_real, 143.728733573_kind_real, &
78 & 146.870307625_kind_real, 150.011882457_kind_real, &
79 & 153.153458019_kind_real, 156.295034268_kind_real /
95 slatd(j)=cos((bz(jz)+(j-jz)*
pi)*r)
112 sp=pk(j)*(
one-slatd(j)**2)/(jmax*(pkm1(j)-slatd(j)*pk(j)))
114 spmax=max(spmax,abs(sp))
119 wlat(j)=(
two*(
one-slatd(j)**2))/(jmax*pkm1(j))**2
120 slat(jmax+1-j)=-slat(j)
121 wlat(jmax+1-j)=wlat(j)
125 wlat(jhe)=
two/jmax**2
127 wlat(jhe)=wlat(jhe)*n**2/(n-1)**2
130 else if(idrt==0)
then
142 awork(js,j)=cos(real(2*(js-1)*j,
kind_real)*dlt)
146 bwork(js)=-d1/real((4*(js-1)**2-1),
kind_real)
148 call ludcmp(awork,jho,jhe,ipvt)
149 call lubksb(awork,jho,jhe,ipvt,bwork)
155 slat(jmax+1-j)=-slat(j)
156 wlat(jmax+1-j)=wlat(j)
160 wlat(jhe)=
two*wlat(jhe)
162 else if(idrt==256)
then
178 bwork(js)=-d1/real(4*(js-1)**2-1,
kind_real)
180 call ludcmp(awork,jho,jhe,ipvt,d)
181 call lubksb(awork,jho,jhe,ipvt,bwork)
187 slat(jmax+1-j)=-slat(j)
188 wlat(jmax+1-j)=wlat(j)
192 wlat(jhe)=
two*wlat(jhe)
207 integer,
intent(in):: NP
208 integer,
intent(in):: N
209 real(kind=
kind_real),
intent(in):: a(np,np)
210 real(kind=
kind_real),
intent(inout):: b(n)
211 integer,
intent(in):: INDX(N)
226 else if(abs(sum)>
zero)
then
252 integer,
intent(in):: N
253 integer,
intent(in):: NP
254 real(kind=
kind_real),
intent(inout):: a(np,np)
255 integer,
intent(out):: INDX(N)
256 real(kind=
kind_real),
intent(out),
optional:: d
259 real(kind=
kind_real),
parameter:: tiny=1.0e-20
269 if(abs(a(i,j))>aamax) aamax=abs(a(i,j))
271 if(.NOT.(abs(aamax)>
zero)) print *,
'SINGULAR MATRIX.'
280 sum=sum-a(i,k)*a(k,j)
291 sum=sum-a(i,k)*a(k,j)
312 if(.NOT.(abs(a(j,j))>
zero)) a(j,j)=tiny
319 if(.NOT.(abs(a(n,n))>
zero)) a(n,n)=tiny