1 ***************************************************************************
2 * *
3 * program eof analysis for January atlantic sst anomalies
4 *
5 * parameter *
6 * *
7 * 1.variable *
8 * sst : sea surface temperature *
9 * mon : month *
10 * iyr : year *
11 *
12 * compile command f77 file.f -llapack -lblas
13 * deback command f77 file.f -o silly -g -llapack -lblas
14 *
15 ***************************************************************************
16 parameter (i2s=1,i2e=59,j2s=1,j2e=60,k2s=26,k2e=60)
17 parameter (iys=1900,iye=1996,ipn=1069,idn=(iye-iys+1-1))
18 parameter (mode=5)
19
20 double precision t0(idn,ipn)
21 double precision a (ipn,ipn)
22 integer i2g(ipn),j2g(ipn)
23
24 open(91,file=
25 &'../dataset/sst.gisst.1900-1996.1x1.monthly.n_atl.dat'
26 & ,form='unformatted',access='direct'
27 & ,recl=_________________________
28 & ,status='old')
29
30 open(92,file='../stat_data/SSTeof101v'
31 & ,form='unformatted',access='direct'
32 & ,recl=(i2e-i2s+1)*(k2e-k2s+1)*4)
33 open(93,file='../stat_data/SSTeof101t'
34 & ,form='unformatted',access='direct'
35 & ,recl=mode*4)
36
37 write(6,*) 'dread(t0,i2g,j2g)'
38 call dread(t0,i2g,j2g)
39
40 write(6,*) 'covar(t0,a)'
41 call covar(t0,a)
42
43 write(6,*) 'eigen(t0,a,i2g,j2g)'
44 call eigen(t0,a,i2g,j2g)
45
46 stop
47 end
48 c///////////////////////////////////////////////////////////////////////
49 subroutine dread(t0,i2g,j2g)
50 parameter (i2s=1,i2e=59,j2s=1,j2e=60,k2s=26,k2e=60)
51 parameter (iys=1900,iye=1996,ipn=1069,idn=(iye-iys+1-1))
52
53 double precision t0(idn,ipn)
54 integer i2g(ipn),j2g(ipn)
55
56 real sst (i2e-i2s+1,j2e-j2s+1,iys+1:iye)
57 real*4 r4 (i2e-i2s+1,j2e-j2s+1)
58 integer n (i2e-i2s+1,j2e-j2s+1)
59
60 do 1010 i2=i2s,i2e
61 do 1010 j2=j2s,j2e
62 n (i2,j2)=0
63 do 1015 iyr=iys+1,iye
64 sst(i2,j2,iyr)=0.0
65 1015 continue
66 1010 continue
67
68 do 1020 iyr=iys+1,iye
69 mon=3
70 ir=12*(iyr-iys)+mon
71 read(91,rec=ir) r4
72
73 do 1025 i2=1,i2e-i2s+1
74 do 1025 j2=1,j2e-j2s+1
75 if(r4(i2,j2).gt.-0.9e21) then
76 sst(i2,j2,iyr)=r4(i2,j2)
77 n (i2,j2) =n (i2,j2)+1
78 endif
79 1025 continue
80 1020 continue
81
82 np=0
83 do 1030 i2=1,i2e-i2s+1
84 do 1030 j2=k2s,k2e
85 c '25' is the northern limit of southern domain
86 if(n (i2,j2).eq.(iye-iys+1-1)) then
87 np =np+1
88 i2g(np)=i2
89 j2g(np)=j2
90 do 1031 iyr=iys+1,iye
91 t0(iyr-iys,np)=dble(______________)
92 1031 continue
93 endif
94 1030 continue
95 write(6,*) np
96 pause
97
98 return
99 end
100 c///////////////////////////////////////////////////////////
101 subroutine covar(t0,a)
102 parameter (iys=1900,iye=1996,ipn=1069,idn=(iye-iys+1-1))
103
104 double precision t0(idn,ipn)
105 double precision a (ipn,ipn)
106 double precision x0,y0,vxy,sumvar,sumcovar
107
108 do 10 _=1,ipn
109 do 10 _=1,ipn
110 vxy=0.0d00
111
112 do 11 k=1,idn
113 x0=t0(k,i)
114 y0=t0(k,j)
115 vxy=________
116 11 continue
117
118 a(i,j)=(vxy)/dble(idn)
119
120 10 continue
121
122 sumvar=0.0
123 nvar =0
124 do 20 i=1,ipn
125 sumvar=sumvar+a(i,i)
126 nvar =nvar+1
127 20 continue
128
129 sumcovar=0.0
130 ncovar =0
131 do 30 i=1,ipn
132 do 30 j=1,ipn
133 if(i.ne.j) then
134 sumcovar=sumcovar+a(i,j)
135 ncovar=ncovar+1
136 endif
137 30 continue
138
139 write(6,*) sumvar ,nvar
140 write(6,*) sumcovar,ncovar,ipn*ipn-ipn
141
142 return
143 end
144 c///////////////////////////////////////////////////////////
145 subroutine eigen(t0,a,i2g,j2g)
146 parameter (i2s=1,i2e=59,j2s=1,j2e=60,k2s=26,k2e=60)
147 parameter (iys=1900,iye=1996,ipn=1069,idn=(iye-iys+1-1))
148 parameter (lwork=ipn*ipn)
149 parameter (mode=5)
150
151 double precision t0(idn,ipn)
152 double precision a (ipn,ipn)
153 integer i2g(ipn),j2g(ipn)
154
155 double precision eval(ipn)
156 double precision work(lwork)
157
158 real*4 score(mode)
159 real*4 mtrx(i2e-i2s+1,k2e-k2s+1)
160
161 real sum,vxy
162
163 c------------------------------------eigen value decomposition
164 n =ipn
165 lda1=ipn
166 call dsyev(___,___,_,_,____,____,____,_____,____)
167
168 write(6,*) info
169 c-------------------------------------------------------------
170 c------------------------------------explained variance of each mode
171 sum=0.0
172 do 10 i=ipn,1,-1
173 sum=sum+eval(i)
174 10 continue
175
176 write(6,'(5f8.5)') (eval(i)/sum,i=ipn,ipn-4,-1)
177 c----------------------------------------------------------------
178 c------------------------------------temporal coefficient of each mode
179 k=0
180 do 20 iyr=iys+1,iye
181 k =k+1
182 do 21 i=ipn,ipn-4,-1
183 vxy=0.0
184 do 22 j=1,ipn
185 c dim.1: grid, dim.2: mode > a(j,i))
186 vxy=vxy+t0(k,j)*a(j,i)
187 22 continue
188 c dim.1: mode, dim.2: time > score(ipn-1+1,k)
189 score(ipn-i+1)=vxy/sqrt(eval(i))
190 21 continue
191 itn=iyr-iys
192 c write(6,*) itn,score
193 write(93,rec=itn) score
194 20 continue
195
196 c----------------------------------------------------------------
197 c-----------------------------------check for squared all eigen vectors =1.0
198 sumev=0.0
199 do 30 j=1,ipn
200 c write(6,'(5f10.5)') (a(j,i),i=ipn,ipn-4,-1)
201 sumev=sumev+a(j,i)*a(j,i)
202 30 continue
203 c write(6,*) 'sumev=', sumev
204 c----------------------------------------------------------------
205 c------------------------------------eigen vectors on lon x lat map
206 do 40 i2=1,i2e-i2s+1
207 do 40 j2=1,k2e-k2s+1
208 mtrx(i2,j2)=-1.0e21
209 40 continue
210
211 iw=0
212 do 50 i=ipn,ipn-4,-1
213 iw=iw+1
214 do 51 j=1,ipn
215 c write(6,'(3i5,f10.5)') j,ig(j),jg(j),a(j,i)
216 mtrx(i2g(j),j2g(j)-k2s+1)=a(j,i)*sqrt(eval(i))
217 51 continue
218
219 write(92,rec=iw) mtrx
220 50 continue
221
222 return
223 end