1 |
adcroft |
1.2 |
C $Header: /u/gcmpack/models/MITgcmUV/utils/knudsen2/knudsen2.f,v 1.1 1998/05/28 16:26:56 adcroft Exp $ |
2 |
adcroft |
1.1 |
|
3 |
|
|
PROGRAM KNUDSEN2 |
4 |
|
|
implicit none |
5 |
|
|
C |
6 |
|
|
C COEFFICIENTS FOR DENSITY COMPUTATION |
7 |
|
|
C |
8 |
|
|
C THIS PROGRAM CALCULATES THE COEFFICIENTS OF THE POLYNOMIAL |
9 |
|
|
C APPROXIMATION TO THE KNUDSEN FORMULA. THE PROGRAM IS SET UP |
10 |
|
|
C TO YIELD COEFFICIENTS THAT WILL COMPUTE DENSITY AS A FUNCTION |
11 |
|
|
C OF POTENTIAL TEMPERATURE, SALINITY, AND DEPTH. |
12 |
|
|
C |
13 |
|
|
C Number of levels |
14 |
|
|
integer NumLevels |
15 |
adcroft |
1.2 |
PARAMETER (NumLevels=15) |
16 |
adcroft |
1.1 |
|
17 |
|
|
C Number of Temperature and Salinity callocation points |
18 |
|
|
integer NumTempPt,NumSalPt,NumCallocPt |
19 |
|
|
PARAMETER (NumTempPt = 10, NumSalPt = 5 ) |
20 |
|
|
PARAMETER (NumCallocPt =NumSalPt*NumTempPt ) |
21 |
|
|
|
22 |
|
|
C Number of terms in fit polynomial |
23 |
|
|
integer NTerms |
24 |
|
|
PARAMETER (NTerms = 9 ) |
25 |
|
|
|
26 |
|
|
C Functions (Density and Potential temperature) |
27 |
|
|
DOUBLE PRECISION DN |
28 |
|
|
real POTEM |
29 |
|
|
C |
30 |
|
|
C Arrays for least squares routine |
31 |
|
|
real A(NumCallocPt,NumCallocPt), |
32 |
|
|
& B(NumCallocPt), |
33 |
|
|
& X(NTerms), |
34 |
|
|
& AB(13,NumLevels), |
35 |
|
|
& C(NumCallocPt,NTerms), |
36 |
|
|
& H(NTerms,NTerms), |
37 |
|
|
& R(NumCallocPt*2),SB(NumCallocPt*4) |
38 |
|
|
|
39 |
|
|
real |
40 |
|
|
& Z(NumLevels), |
41 |
|
|
& DZ(NumLevels) |
42 |
|
|
integer |
43 |
|
|
& IDX(NumLevels) |
44 |
|
|
|
45 |
|
|
real |
46 |
|
|
& Theta(NumCallocPt), |
47 |
|
|
& SPick(NumCallocPt),TPick(NumCallocPt),DenPick(NumCallocPt) |
48 |
|
|
|
49 |
|
|
C ENTER BOUNDS FOR POLYNOMIAL FIT: |
50 |
|
|
C TMin(K) = LOWER BND OF T AT Level K (Insitu Temperature) |
51 |
|
|
C TMax(K) = UPPER BND OF T " (Insitu Temperature) |
52 |
|
|
C SMin(K) = LOWER BND OF S " |
53 |
|
|
C SMax(K) = UPPER BND OF S " |
54 |
|
|
|
55 |
|
|
real TMin(25), TMax(25) |
56 |
|
|
& ,SMin(25), SMax(25) |
57 |
|
|
|
58 |
|
|
|
59 |
|
|
DATA TMin / 4*-2., 15*-1., 6*0. / |
60 |
|
|
DATA TMax / 29., 19., 14., 11., 9., 7., 5., 3*4., 5*3., 10*2. / |
61 |
|
|
DATA SMin / 28.5, 33.7, 34., 34.1, 34.2, 34.4, 2*34.5, |
62 |
|
|
& 15*34.6, 2*34.7 / |
63 |
|
|
DATA SMax / 36.7, 36.6, 35.8, 35.7, 35.3, 2*35.1, 7*35., |
64 |
|
|
& 9*34.9, 2*34.8 / |
65 |
|
|
|
66 |
|
|
! Local |
67 |
|
|
integer I,J,K,IT,IEQ,IRANK,NDIM,NHDIM,N,M,IN,ITMAX,MP,L,NN |
68 |
|
|
real T,S,D,EPS,DeltaT,DeltaS,ENORM,DeltaDen,DensityRef |
69 |
|
|
real Tbar,Sbar,Tsum,SSum,TempInc,SalnInc |
70 |
|
|
real DensitySum,ThetaSum |
71 |
|
|
real DensityBar,thetabar,TempRef,SAlRef |
72 |
|
|
|
73 |
|
|
|
74 |
|
|
C ENTER LEVEL THICKNESSES IN CENTIMETERS |
75 |
|
|
C DATA dz/ 4*50.E2, 2*100.E2, 200.E2, 400.E2, 7*500.E2, 6*10.E2 / |
76 |
adcroft |
1.2 |
DATA dz/ 5.00e+03,7.00e+03,1.00e+04,1.40e+04, |
77 |
|
|
&1.90e+04,2.40e+04,2.90e+04,3.40e+04,3.90e+04, |
78 |
|
|
&4.40e+04,4.90e+04,5.40e+04,5.90e+04,6.40e+04, |
79 |
|
|
&6.90e+04/ ! ,6.90e+04/ |
80 |
adcroft |
1.1 |
|
81 |
|
|
C CALC DEPTHS OF LEVELS FROM DZ (IN METERS) |
82 |
|
|
C THE MAXIMUM ALLOWABLE DEPTH IS 8000 METERS |
83 |
|
|
Z(1) = 0. ! for level at box edges |
84 |
|
|
Z(1)=.5*DZ(1)/100. ! for level at box center |
85 |
|
|
DO K=2,NumLevels |
86 |
|
|
Z(K)=Z(K-1)+.5*(DZ(K)+DZ(K-1))/100. |
87 |
|
|
ENDDO |
88 |
|
|
|
89 |
|
|
|
90 |
|
|
C Break the levels up into 250m bands |
91 |
|
|
DO I=1,NumLevels |
92 |
|
|
IDX( I ) = I |
93 |
|
|
C Comment out the next line to use input bands as polynomial levels |
94 |
|
|
IDX( I ) = IFIX(Z(I)/250.)+1 |
95 |
|
|
ENDDO |
96 |
|
|
|
97 |
|
|
|
98 |
|
|
C Write some diagnostics |
99 |
|
|
PRINT 419 |
100 |
|
|
PRINT 422,(Z(I), |
101 |
|
|
& Tmin( IDX(I)),Tmax( IDX(I)), |
102 |
|
|
& SMin(IDX(I)),Smax( IDX(I)), |
103 |
|
|
CCC & TMin(I),TMax(I),SMin(I),SMax(I), |
104 |
|
|
& ( IDX(I) ),I=1,NumLevels ) |
105 |
|
|
|
106 |
|
|
|
107 |
|
|
C Loop over each level and calculate the polynomial coefficients |
108 |
|
|
|
109 |
|
|
DO MP=1,NumLevels |
110 |
|
|
|
111 |
|
|
C Choose callocation points |
112 |
|
|
TempInc =(Tmax( IDX(MP))-TMin(IDX(MP)))/(2.*FLOAT(NumSalPt)-1.0) |
113 |
|
|
SalnInc =(Smax( IDX(MP))-SMin(IDX(MP)))/(FLOAT(NumSalPt)-1.0) |
114 |
|
|
DO I=1,NumTempPt |
115 |
|
|
DO J=1,NumSalPt |
116 |
|
|
K=NumSalPt*I+J-NumSalPt |
117 |
|
|
TPick(K)=TMin(IDX(MP))+(FLOAT(I)-1.0)*TempInc |
118 |
|
|
SPick(K)=SMin(IDX(MP))+(FLOAT(J)-1.0)*SalnInc |
119 |
|
|
ENDDO |
120 |
|
|
ENDDO |
121 |
|
|
|
122 |
|
|
|
123 |
|
|
C For each callocation point, convert insitu temperature to |
124 |
|
|
C potential temperature and calculate the corresponding density. |
125 |
|
|
Tsum=0.0 |
126 |
|
|
Ssum=0.0 |
127 |
|
|
DensitySum=0.0 |
128 |
|
|
ThetaSum=0.0 |
129 |
|
|
DO K=1,NumCallocPt |
130 |
|
|
D=Z(MP) |
131 |
|
|
S=SPick(K) |
132 |
|
|
T=TPick(K) |
133 |
|
|
DenPick(K)=DN(T,S,D) |
134 |
|
|
Theta(K)=POTEM(T,S,D) |
135 |
|
|
Tsum=Tsum+TPick(K) |
136 |
|
|
Ssum=Ssum+SPick(K) |
137 |
|
|
DensitySum = DensitySum+DenPick(K) |
138 |
|
|
ThetaSum = ThetaSum+Theta(K) |
139 |
|
|
ENDDO |
140 |
|
|
|
141 |
|
|
C Let (Tbar,Sbar) = the average of (T,S) used in the set of calloc pts |
142 |
|
|
C NOTE: Tbar is still an insitu temperature |
143 |
|
|
C Also, calculate the average density of the set of callocation points |
144 |
|
|
|
145 |
|
|
Tbar=Tsum / FLOAT( NumCallocPt ) |
146 |
|
|
Sbar=Ssum / FLOAT( NumCallocPt ) |
147 |
|
|
DensityBar=DensitySum / FLOAT( NumCallocPt ) |
148 |
|
|
|
149 |
|
|
C Calculate the average potential temperature of the callocation points |
150 |
|
|
ThetaBar=ThetaSum/ FLOAT( NumCallocPt ) |
151 |
|
|
|
152 |
|
|
C Set the reference temperature, salinity and density |
153 |
|
|
DensityRef=DN(Tbar,Sbar,D) |
154 |
|
|
|
155 |
|
|
SalRef = Sbar |
156 |
|
|
|
157 |
|
|
TempRef=TBar |
158 |
|
|
TempRef=ThetaBar ! DELETE THIS LINE IF USING IN SITU TEMPERATURES |
159 |
|
|
|
160 |
|
|
C$$$ TempRef=POTEM(TBar,Sbar,D) |
161 |
|
|
|
162 |
|
|
|
163 |
|
|
AB(1,MP)=Z(MP) |
164 |
|
|
AB(2,MP)=DensityRef |
165 |
|
|
AB(3,MP)=TempRef |
166 |
|
|
AB(4,MP)=SalRef |
167 |
|
|
DO K=1,NumCallocPt |
168 |
|
|
TPick(K)=Theta(K) ! DELETE THIS LINE IF USING IN SITU TEMPERATURES |
169 |
|
|
DeltaT = TPick(K) - TempRef |
170 |
|
|
DeltaS = SPick(K) - SalRef |
171 |
|
|
DeltaDen = DenPick(K) - DensityRef |
172 |
|
|
B(K)= DeltaDen |
173 |
|
|
A(K,1)=DeltaT |
174 |
|
|
A(K,2)=DeltaS |
175 |
|
|
A(K,3)=DeltaT*DeltaT |
176 |
|
|
A(K,4)=DeltaT*DeltaS |
177 |
|
|
A(K,5)=DeltaS*DeltaS |
178 |
|
|
A(K,6)=A(K,3)*DeltaT |
179 |
|
|
A(K,7)=A(K,4)*DeltaT |
180 |
|
|
A(K,8)=A(K,4)*DeltaS |
181 |
|
|
A(K,9)=A(K,5)*DeltaS |
182 |
|
|
ENDDO |
183 |
|
|
C SET THE ARGUMENTS IN CALL TO LSQSL2 |
184 |
|
|
C FIRST DIMENSION OF ARRAY A |
185 |
|
|
NDIM=NumCallocPt |
186 |
|
|
C |
187 |
|
|
C NUMBER OF ROWS OF A |
188 |
|
|
M=NumCallocPt |
189 |
|
|
C |
190 |
|
|
C NUMBER OF COLUMNS OF A |
191 |
|
|
N=NTerms |
192 |
|
|
C OPTION NUMBER OF LSQSL2 |
193 |
|
|
IN=1 |
194 |
|
|
C |
195 |
|
|
C ITMAX=NUMBER OF ITERATIONS |
196 |
|
|
ITMAX=100 |
197 |
|
|
C |
198 |
|
|
IT=0 |
199 |
|
|
IEQ=2 |
200 |
|
|
IRANK=0 |
201 |
|
|
EPS=1.0E-7 |
202 |
|
|
EPS=1.0E-11 |
203 |
|
|
NHDIM=NTerms |
204 |
|
|
CALL LSQSL2(NDIM,A,M,N,B,X,IRANK,IN,ITMAX,IT,IEQ,ENORM,EPS, |
205 |
|
|
& NHDIM,H,C,R,SB) |
206 |
|
|
|
207 |
|
|
CCC PRINT 411 |
208 |
|
|
DO I=1,N |
209 |
|
|
AB(I+4,MP)=X(I) |
210 |
|
|
ENDDO |
211 |
|
|
ENDDO |
212 |
|
|
PRINT 430 |
213 |
|
|
430 FORMAT(1X,' Z SIG0 T S X1 X2 |
214 |
|
|
1 X3 X4 X5 X6 X7 X8 |
215 |
|
|
2 X9 ',/) |
216 |
|
|
NN=N+4 |
217 |
|
|
ccc C************************************************************************ |
218 |
|
|
ccc C WRITE TO UNIT 50 |
219 |
|
|
ccc C************************************************************************ |
220 |
|
|
ccc OPEN(UNIT=50,FILE='KNUDSEN_COEFS.mit.h',STATUS='UNKNOWN') |
221 |
|
|
ccc |
222 |
|
|
ccc C*** WRITE(50,613) |
223 |
|
|
ccc 613 FORMAT(' REAL SIGREF(KM)') |
224 |
|
|
ccc ISEQ=114399990 |
225 |
|
|
ccc DO J=1,NumLevels |
226 |
|
|
ccc ISEQ=ISEQ+10 |
227 |
|
|
ccc C$$$ WRITE(50,615)J,.01*DZ(J) |
228 |
|
|
ccc ENDDO |
229 |
|
|
ccc 615 FORMAT(' DZ(',I3,')=',F7.2,'E2') |
230 |
|
|
ccc C$$$ WRITE(50,601) |
231 |
|
|
ccc 601 FORMAT('C REFERENCE DENSITIES AT T-POINTS' ) |
232 |
|
|
ccc ISEQ=114500030 |
233 |
|
|
ccc DO J=1,NumLevels |
234 |
|
|
ccc ISEQ=ISEQ+10 |
235 |
|
|
ccc C$$$ WRITE(50,501)J,AB(2,J) |
236 |
|
|
ccc ENDDO |
237 |
|
|
ccc 501 FORMAT(' SIGREF(',I3,')=',F8.4) |
238 |
|
|
ccc |
239 |
|
|
ccc |
240 |
|
|
ccc WRITE(50,'(" DATA SIGREF/")') |
241 |
|
|
ccc |
242 |
|
|
ccc N0 = 1 |
243 |
|
|
ccc c DO ILINE = 1,NumLevels/5 -1 |
244 |
|
|
ccc 222 CONTINUE |
245 |
|
|
ccc N1 = N0 + 4 |
246 |
|
|
ccc N1 = MIN(N1, NumLevels) |
247 |
|
|
ccc WRITE(50,1280) (AB(2,I),I=N0,N1) |
248 |
|
|
ccc N0 = N1 + 1 |
249 |
|
|
ccc IF ( N1 .LT. NumLevels ) GOTO 222 |
250 |
|
|
ccc C N1 = N0 + 4 |
251 |
|
|
ccc C IF( N1 .GT. NumLevels ) N1 = NumLevels |
252 |
|
|
ccc WRITE(50,1286) |
253 |
|
|
ccc |
254 |
|
|
|
255 |
|
|
C********************************************************************** |
256 |
|
|
|
257 |
|
|
C MITgcmUV |
258 |
|
|
|
259 |
|
|
open(99,file='POLY3.COEFFS',form='formatted',status='unknown') |
260 |
|
|
write(99,*) NumLevels |
261 |
|
|
write(99,*) (ab(3,J),ab(4,J),ab(2,J),J=1 , NumLevels) |
262 |
|
|
do J=1,NumLevels |
263 |
|
|
write(99,*)(AB(I,J),I=5,13) |
264 |
|
|
enddo |
265 |
|
|
close(99) |
266 |
|
|
|
267 |
|
|
C********************************************************************** |
268 |
|
|
|
269 |
|
|
c PRINT 412,((AB(I,J),I=1,NN),J=1,NumLevels) |
270 |
|
|
|
271 |
|
|
C WRITE DATA STATEMENTS TO UNIT 50 |
272 |
|
|
caja DO L=1,NumLevels |
273 |
|
|
caja AB(2,L)=1.E-3*AB(2,L) |
274 |
|
|
caja AB(4,L)=1.E-3*AB(4,L)-.035 |
275 |
|
|
caja AB(5,L)=1.E-3*AB(5,L) |
276 |
|
|
caja AB(7,L)=1.E-3*AB(7,L) |
277 |
|
|
caja AB(10,L)=1.E-3*AB(10,L) |
278 |
|
|
caja AB( 9,L)=1.E+3*AB( 9,L) |
279 |
|
|
caja AB(11,L)=1.E+3*AB(11,L) |
280 |
|
|
caja AB(13,L)=1.E+6*AB(13,L) |
281 |
|
|
caja ENDDO |
282 |
|
|
|
283 |
|
|
C********************************************************************** |
284 |
|
|
|
285 |
|
|
C MITgcm (compare01) |
286 |
|
|
|
287 |
|
|
open(99,file='polyeos.coeffs',form='formatted',status='unknown') |
288 |
|
|
do J=1,NumLevels |
289 |
|
|
write(99,*) ab(3,J) ! ref temperature |
290 |
|
|
enddo |
291 |
|
|
do J=1,NumLevels |
292 |
|
|
write(99,*) ab(4,J) ! ref sal |
293 |
|
|
enddo |
294 |
|
|
do J=1,NumLevels |
295 |
|
|
do I=5,13 |
296 |
|
|
write(99,*) AB(I,J) |
297 |
|
|
enddo |
298 |
|
|
enddo |
299 |
|
|
do J=1,NumLevels |
300 |
|
|
write(99,*) ab(2,J) ! ref sig0 |
301 |
|
|
enddo |
302 |
|
|
CEK write(99,200)(ab(3,J),J=11,15) ! ref temperature |
303 |
|
|
CEK write(99,200)(ab(4,J),J= 1, 5) ! ref sal |
304 |
|
|
CEK write(99,200)(ab(4,J),J= 6,10) ! ref sal |
305 |
|
|
CEK write(99,200)(ab(4,J),J= 11,15) ! ref sal |
306 |
|
|
CEK do I=5,NN |
307 |
|
|
CEK write(99,200)(AB(I,J),J= 1, 5) |
308 |
|
|
CEK write(99,200)(AB(I,J),J= 6,10) |
309 |
|
|
CEK write(99,200)(AB(I,J),J=11,15) |
310 |
|
|
CEK enddo |
311 |
|
|
CEK write(99,200)(AB(2,J),J= 1, 5) ! ref sig0 |
312 |
|
|
CEK write(99,200)(AB(2,J),J= 6,10) ! ref sig0 |
313 |
|
|
CEK write(99,200)(AB(2,J),J=11,15) ! ref sig0 |
314 |
|
|
CEK 200 format(5(E14.7,1X)) |
315 |
|
|
|
316 |
|
|
C********************************************************************** |
317 |
|
|
|
318 |
|
|
|
319 |
|
|
|
320 |
|
|
ccc NSEQ=603800000 |
321 |
|
|
ccc WRITE(50,1298) |
322 |
|
|
ccc N=0 |
323 |
|
|
ccc 1260 IS=N+1 |
324 |
|
|
ccc IE=N+5 |
325 |
|
|
ccc NSEQ=NSEQ+10 |
326 |
|
|
ccc IF(IE.LT.NumLevels) THEN |
327 |
|
|
ccc WRITE(50,1280) (AB(3,I),I=IS,IE) |
328 |
|
|
ccc N=IE |
329 |
|
|
ccc GO TO 1260 |
330 |
|
|
ccc ELSE |
331 |
|
|
ccc IE=NumLevels |
332 |
|
|
ccc N=IE-IS+1 |
333 |
|
|
ccc GO TO (1261,1262,1263,1264,1265),N |
334 |
|
|
ccc 1261 WRITE(50,1281) (AB(3,I),I=IS,IE) |
335 |
|
|
ccc GO TO 1268 |
336 |
|
|
ccc 1262 WRITE(50,1282) (AB(3,I),I=IS,IE) |
337 |
|
|
ccc GO TO 1268 |
338 |
|
|
ccc 1263 WRITE(50,1283) (AB(3,I),I=IS,IE) |
339 |
|
|
ccc GO TO 1268 |
340 |
|
|
ccc 1264 WRITE(50,1284) (AB(3,I),I=IS,IE) |
341 |
|
|
ccc GO TO 1268 |
342 |
|
|
ccc 1265 WRITE(50,1285) (AB(3,I),I=IS,IE) |
343 |
|
|
ccc ENDIF |
344 |
|
|
ccc 1268 CONTINUE |
345 |
|
|
ccc NSEQ=603900000 |
346 |
|
|
ccc WRITE(50,1297) |
347 |
|
|
ccc N=0 |
348 |
|
|
ccc 1270 IS=N+1 |
349 |
|
|
ccc IE=N+5 |
350 |
|
|
ccc NSEQ=NSEQ+10 |
351 |
|
|
ccc IF(IE.LT.NumLevels) THEN |
352 |
|
|
ccc WRITE(50,1280) (AB(4,I),I=IS,IE) |
353 |
|
|
ccc N=IE |
354 |
|
|
ccc GO TO 1270 |
355 |
|
|
ccc ELSE |
356 |
|
|
ccc IE=NumLevels |
357 |
|
|
ccc N=IE-IS+1 |
358 |
|
|
ccc GO TO (1271,1272,1273,1274,1275),N |
359 |
|
|
ccc 1271 WRITE(50,1281) (AB(4,I),I=IS,IE) |
360 |
|
|
ccc GO TO 1278 |
361 |
|
|
ccc 1272 WRITE(50,1282) (AB(4,I),I=IS,IE) |
362 |
|
|
ccc GO TO 1278 |
363 |
|
|
ccc 1273 WRITE(50,1283) (AB(4,I),I=IS,IE) |
364 |
|
|
ccc GO TO 1278 |
365 |
|
|
ccc 1274 WRITE(50,1284) (AB(4,I),I=IS,IE) |
366 |
|
|
ccc GO TO 1278 |
367 |
|
|
ccc 1275 WRITE(50,1285) (AB(4,I),I=IS,IE) |
368 |
|
|
ccc ENDIF |
369 |
|
|
ccc 1278 CONTINUE |
370 |
|
|
ccc DO 1200 L=1,NumLevels |
371 |
|
|
ccc IF(L.EQ.1) NSEQ=604000000 |
372 |
|
|
ccc WRITE(50,1296) L |
373 |
|
|
ccc NSEQ=NSEQ+10 |
374 |
|
|
ccc WRITE(50,1295) (AB(I,L),I=5,8) |
375 |
|
|
ccc NSEQ=NSEQ+10 |
376 |
|
|
ccc WRITE(50,1295) (AB(I,L),I=9,12) |
377 |
|
|
ccc NSEQ=NSEQ+10 |
378 |
|
|
ccc WRITE(50,1294) AB(13,L) |
379 |
|
|
ccc NSEQ=NSEQ+10 |
380 |
|
|
ccc 1200 CONTINUE |
381 |
|
|
ccc 1288 CONTINUE |
382 |
|
|
1298 FORMAT(18H DATA TRef /,67X,I9) |
383 |
|
|
1297 FORMAT(18H DATA SRef /,67X,I9) |
384 |
|
|
C 419 FORMAT(5X,'LEVEL TMIN TMAX SMIN SMAX ', |
385 |
|
|
C & ' TMIN() Tmax() Smin() Smax() D/250') |
386 |
|
|
C 422 FORMAT (5X,F6.1,4F10.3,4F10.3,I3) |
387 |
|
|
419 FORMAT(5X,'LEVEL TMIN TMAX SMIN SMAX D/250') |
388 |
|
|
422 FORMAT (5X,F6.1,4F10.3,I3) |
389 |
|
|
412 FORMAT(1X,F6.1,F8.4,F7.3,F6.2,9E12.5) |
390 |
|
|
1296 FORMAT(6X,'DATA (C(',I2,',N),N=1,9)/',54X,I9) |
391 |
|
|
1295 FORMAT(5X,1H*,9X,4(E13.7,1H,),10X,I9) |
392 |
|
|
1294 FORMAT(5X,1H*,9X,E13.7,1H/,52X,I9) |
393 |
|
|
1280 FORMAT(5X,1H*,8X,5(F10.7,1H,),12X,I9) |
394 |
|
|
1281 FORMAT(5X,1H*,8X,F10.7,1H/,56X,I9) |
395 |
|
|
1282 FORMAT(5X,1H*,8X,F10.7,1H,,F10.7,1H/,45X,I9) |
396 |
|
|
1283 FORMAT(5X,1H*,8X,2(F10.7,1H,),F10.7,1H/,34X,I9) |
397 |
|
|
1284 FORMAT(5X,1H*,8X,3(F10.7,1H,),F10.7,1H/,23X,I9) |
398 |
|
|
1285 FORMAT(5X,1H*,8X,4(F10.7,1H,),F10.7,1H/,12X,I9) |
399 |
|
|
1286 FORMAT(5X,1H*,1H/) |
400 |
|
|
350 FORMAT(1X,E14.7) |
401 |
|
|
351 FORMAT(1H+,T86,E14.7/) |
402 |
|
|
400 FORMAT(1X,9E14.7) |
403 |
|
|
410 FORMAT(1X,5E14.7) |
404 |
|
|
411 FORMAT(///) |
405 |
|
|
STOP |
406 |
|
|
END |
407 |
|
|
C**************************************************************************** |
408 |
|
|
*DECK LSQSL2 |
409 |
|
|
SUBROUTINE LSQSL2 (NDIM,A,D,W,B,X,IRANK,IN,ITMAX,IT,IEQ,ENORM,EPS1 |
410 |
|
|
1,NHDIM,H,AA,R,S) |
411 |
|
|
implicit none |
412 |
|
|
C THIS ROUTINE IS A MODIFICATION OF LSQSOL. MARCH,1968. R. HANSON. |
413 |
|
|
C LINEAR LEAST SQUARES SOLUTION |
414 |
|
|
C |
415 |
|
|
C THIS ROUTINE FINDS X SUCH THAT THE EUCLIDEAN LENGTH OF |
416 |
|
|
C (*) AX-B IS A MINIMUM. |
417 |
|
|
C |
418 |
|
|
C HERE A HAS K ROWS AND N COLUMNS, WHILE B IS A COLUMN VECTOR WITH |
419 |
|
|
C K COMPONENTS. |
420 |
|
|
C |
421 |
|
|
C AN ORTHOGONAL MATRIX Q IS FOUND SO THAT QA IS ZERO BELOW |
422 |
|
|
C THE MAIN DIAGONAL. |
423 |
|
|
C SUPPOSE THAT RANK (A)=R |
424 |
|
|
C AN ORTHOGONAL MATRIX S IS FOUND SUCH THAT |
425 |
|
|
C QAS=T IS AN R X N UPPER TRIANGULAR MATRIX WHOSE LAST N-R COLUMNS |
426 |
|
|
C ARE ZERO. |
427 |
|
|
C THE SYSTEM TZ=C (C THE FIRST R COMPONENTS OF QB) IS THEN |
428 |
|
|
C SOLVED. WITH W=SZ, THE SOLUTION MAY BE EXPRESSED |
429 |
|
|
C AS X = W + SY, WHERE W IS THE SOLUTION OF (*) OF MINIMUM EUCLID- |
430 |
|
|
C EAN LENGTH AND Y IS ANY SOLUTION TO (QAS)Y=TY=0. |
431 |
|
|
C |
432 |
|
|
C ITERATIVE IMPROVEMENTS ARE CALCULATED USING RESIDUALS AND |
433 |
|
|
C THE ABOVE PROCEDURES WITH B REPLACED BY B-AX, WHERE X IS AN |
434 |
|
|
C APPROXIMATE SOLUTION. |
435 |
|
|
C |
436 |
|
|
integer ndim,nhdim |
437 |
|
|
DOUBLE PRECISION SJ,DP,DP1,UP,BP,AJ |
438 |
|
|
LOGICAL ERM |
439 |
|
|
INTEGER D,W |
440 |
|
|
C |
441 |
|
|
C IN=1 FOR FIRST ENTRY. |
442 |
|
|
C A IS DECOMPOSED AND SAVED. AX-B IS SOLVED. |
443 |
|
|
C IN = 2 FOR SUBSEQUENT ENTRIES WITH A NEW VECTOR B. |
444 |
|
|
C IN=3 TO RESTORE A FROM THE PREVIOUS ENTRY. |
445 |
|
|
C IN=4 TO CONTINUE THE ITERATIVE IMPROVEMENT FOR THIS SYSTEM. |
446 |
|
|
C IN = 5 TO CALCULATE SOLUTIONS TO AX=0, THEN STORE IN THE ARRAY H. |
447 |
|
|
C IN = 6 DO NOT STORE A IN AA. OBTAIN T = QAS, WHERE T IS |
448 |
|
|
C MIN(K,N) X MIN(K,N) AND UPPER TRIANGULAR. NOW RETURN.DO NOT OBTAIN |
449 |
|
|
C A SOLUTION. |
450 |
|
|
C NO SCALING OR COLUMN INTERCHANGES ARE PERFORMED. |
451 |
|
|
C IN = 7 SAME AS WITH IN = 6 EXCEPT THAT SOLN. OF MIN. LENGTH |
452 |
|
|
C IS PLACED INTO X. NO ITERATIVE REFINEMENT. NOW RETURN. |
453 |
|
|
C COLUMN INTERCHANGES ARE PERFORMED. NO SCALING IS PERFORMED. |
454 |
|
|
C IN = 8 SET ADDRESSES. NOW RETURN. |
455 |
|
|
C |
456 |
|
|
C OPTIONS FOR COMPUTING A MATRIX PRODUCT Y*H OR H*Y ARE |
457 |
|
|
C AVAILABLE WITH THE USE OF THE ENTRY POINTS MYH AND MHY. |
458 |
|
|
C USE OF THESE OPTIONS IN THESE ENTRY POINTS ALLOW A GREAT SAVING IN |
459 |
|
|
C STORAGE REQUIRED. |
460 |
|
|
C |
461 |
|
|
C |
462 |
|
|
real A(NDIM,NDIM),B(1),AA(D,W),S(1), X(1),H(NHDIM,NHDIM),R(1) |
463 |
|
|
C D = DEPTH OF MATRIX. |
464 |
|
|
C W = WIDTH OF MATRIX. |
465 |
|
|
C---- |
466 |
|
|
integer K,N,IT,ISW,L,M,IRANK,IEQ,IN,K1 |
467 |
|
|
integer J1,J2,J3,J4,J5,J6,J7,J8,J9 |
468 |
|
|
integer N1,N2,N3,N4,N5,N6,N7,N8,NS |
469 |
|
|
integer I,ITMAX,IPM1,II,LM,J,IP,KM,IPP1,IRP1,IRM1 |
470 |
|
|
real SP,ENORM,TOP,TOP1,ENM1,TOP2,EPS1,EPS2,A1,A2,AM |
471 |
|
|
C---- |
472 |
|
|
K=D |
473 |
|
|
N=W |
474 |
|
|
ERM=.TRUE. |
475 |
|
|
C |
476 |
|
|
C IF IT=0 ON ENTRY, THE POSSIBLE ERROR MESSAGE WILL BE SUPPRESSED. |
477 |
|
|
C |
478 |
|
|
IF (IT.EQ.0) ERM=.FALSE. |
479 |
|
|
C |
480 |
|
|
C IEQ = 2 IF COLUMN SCALING BY LEAST MAX. COLUMN LENGTH IS |
481 |
|
|
C TO BE PERFORMED. |
482 |
|
|
C |
483 |
|
|
C IEQ = 1 IF SCALING OF ALL COMPONENTS IS TO BE DONE WITH |
484 |
|
|
C THE SCALAR MAX(ABS(AIJ))/K*N. |
485 |
|
|
C |
486 |
|
|
C IEQ = 3 IF COLUMN SCALING AS WITH IN =2 WILL BE RETAINED IN |
487 |
|
|
C RANK DEFICIENT CASES. |
488 |
|
|
C |
489 |
|
|
C THE ARRAY S MUST CONTAIN AT LEAST MAX(K,N) + 4N + 4MIN(K,N) CELLS |
490 |
|
|
C THE ARRAY R MUST CONTAIN K+4N S.P. CELLS. |
491 |
|
|
C |
492 |
|
|
DATA EPS2/1.E-16/ |
493 |
|
|
C THE LAST CARD CONTROLS DESIRED RELATIVE ACCURACY. |
494 |
|
|
C EPS1 CONTROLS (EPS) RANK. |
495 |
|
|
C |
496 |
|
|
ISW=1 |
497 |
|
|
L=MIN0(K,N) |
498 |
|
|
M=MAX0(K,N) |
499 |
|
|
J1=M |
500 |
|
|
J2=N+J1 |
501 |
|
|
J3=J2+N |
502 |
|
|
J4=J3+L |
503 |
|
|
J5=J4+L |
504 |
|
|
J6=J5+L |
505 |
|
|
J7=J6+L |
506 |
|
|
J8=J7+N |
507 |
|
|
J9=J8+N |
508 |
|
|
LM=L |
509 |
|
|
IF (IRANK.GE.1.AND.IRANK.LE.L) LM=IRANK |
510 |
|
|
IF (IN.EQ.6) LM=L |
511 |
|
|
IF (IN.EQ.8) RETURN |
512 |
|
|
C |
513 |
|
|
C RETURN AFTER SETTING ADDRESSES WHEN IN=8. |
514 |
|
|
C |
515 |
|
|
GO TO (10,360,810,390,830,10,10), IN |
516 |
|
|
C |
517 |
|
|
C EQUILIBRATE COLUMNS OF A (1)-(2). |
518 |
|
|
C |
519 |
|
|
C (1) |
520 |
|
|
C |
521 |
|
|
10 CONTINUE |
522 |
|
|
C |
523 |
|
|
C SAVE DATA WHEN IN = 1. |
524 |
|
|
C |
525 |
|
|
IF (IN.GT.5) GO TO 30 |
526 |
|
|
DO 20 J=1,N |
527 |
|
|
DO 20 I=1,K |
528 |
|
|
20 AA(I,J)=A(I,J) |
529 |
|
|
30 CONTINUE |
530 |
|
|
IF (IEQ.EQ.1) GO TO 60 |
531 |
|
|
DO 50 J=1,N |
532 |
|
|
AM=0.E0 |
533 |
|
|
DO 40 I=1,K |
534 |
|
|
40 AM=AMAX1(AM,ABS(A(I,J))) |
535 |
|
|
C |
536 |
|
|
C S(M+N+1)-S(M+2N) CONTAINS SCALING FOR OUTPUT VARIABLES. |
537 |
|
|
C |
538 |
|
|
N2=J2+J |
539 |
|
|
IF (IN.EQ.6) AM=1.E0 |
540 |
|
|
S(N2)=1.E0/AM |
541 |
|
|
DO 50 I=1,K |
542 |
|
|
50 A(I,J)=A(I,J)*S(N2) |
543 |
|
|
GO TO 100 |
544 |
|
|
60 AM=0.E0 |
545 |
|
|
DO 70 J=1,N |
546 |
|
|
DO 70 I=1,K |
547 |
|
|
70 AM=AMAX1(AM,ABS(A(I,J))) |
548 |
|
|
AM=AM/FLOAT(K*N) |
549 |
|
|
IF (IN.EQ.6) AM=1.E0 |
550 |
|
|
DO 80 J=1,N |
551 |
|
|
N2=J2+J |
552 |
|
|
80 S(N2)=1.E0/AM |
553 |
|
|
DO 90 J=1,N |
554 |
|
|
N2=J2+J |
555 |
|
|
DO 90 I=1,K |
556 |
|
|
90 A(I,J)=A(I,J)*S(N2) |
557 |
|
|
C COMPUTE COLUMN LENGTHS WITH D.P. SUMS FINALLY ROUNDED TO S.P. |
558 |
|
|
C |
559 |
|
|
C (2) |
560 |
|
|
C |
561 |
|
|
100 DO 110 J=1,N |
562 |
|
|
N7=J7+J |
563 |
|
|
N2=J2+J |
564 |
|
|
110 S(N7)=S(N2) |
565 |
|
|
C |
566 |
|
|
C S(M+1)-S(M+ N) CONTAINS VARIABLE PERMUTATIONS. |
567 |
|
|
C |
568 |
|
|
C SET PERMUTATION TO IDENTITY. |
569 |
|
|
C |
570 |
|
|
DO 120 J=1,N |
571 |
|
|
N1=J1+J |
572 |
|
|
120 S(N1)=J |
573 |
|
|
C |
574 |
|
|
C BEGIN ELIMINATION ON THE MATRIX A WITH ORTHOGONAL MATRICES . |
575 |
|
|
C |
576 |
|
|
C IP=PIVOT ROW |
577 |
|
|
C |
578 |
|
|
DO 250 IP=1,LM |
579 |
|
|
C |
580 |
|
|
C |
581 |
|
|
DP=0.D0 |
582 |
|
|
KM=IP |
583 |
|
|
DO 140 J=IP,N |
584 |
|
|
SJ=0.D0 |
585 |
|
|
DO 130 I=IP,K |
586 |
|
|
SJ=SJ+A(I,J)**2 |
587 |
|
|
130 CONTINUE |
588 |
|
|
IF (DP.GT.SJ) GO TO 140 |
589 |
|
|
DP=SJ |
590 |
|
|
KM=J |
591 |
|
|
IF (IN.EQ.6) GO TO 160 |
592 |
|
|
140 CONTINUE |
593 |
|
|
C |
594 |
|
|
C MAXIMIZE (SIGMA)**2 BY COLUMN INTERCHANGE. |
595 |
|
|
C |
596 |
|
|
C SUPRESS COLUMN INTERCHANGES WHEN IN=6. |
597 |
|
|
C |
598 |
|
|
C |
599 |
|
|
C EXCHANGE COLUMNS IF NECESSARY. |
600 |
|
|
C |
601 |
|
|
IF (KM.EQ.IP) GO TO 160 |
602 |
|
|
DO 150 I=1,K |
603 |
|
|
A1=A(I,IP) |
604 |
|
|
A(I,IP)=A(I,KM) |
605 |
|
|
150 A(I,KM)=A1 |
606 |
|
|
C |
607 |
|
|
C RECORD PERMUTATION AND EXCHANGE SQUARES OF COLUMN LENGTHS. |
608 |
|
|
C |
609 |
|
|
N1=J1+KM |
610 |
|
|
A1=S(N1) |
611 |
|
|
N2=J1+IP |
612 |
|
|
S(N1)=S(N2) |
613 |
|
|
S(N2)=A1 |
614 |
|
|
N7=J7+KM |
615 |
|
|
N8=J7+IP |
616 |
|
|
A1=S(N7) |
617 |
|
|
S(N7)=S(N8) |
618 |
|
|
S(N8)=A1 |
619 |
|
|
160 IF (IP.EQ.1) GO TO 180 |
620 |
|
|
A1=0.E0 |
621 |
|
|
IPM1=IP-1 |
622 |
|
|
DO 170 I=1,IPM1 |
623 |
|
|
A1=A1+A(I,IP)**2 |
624 |
|
|
170 CONTINUE |
625 |
|
|
IF (A1.GT.0.E0) GO TO 190 |
626 |
|
|
180 IF (DP.GT.0.D0) GO TO 200 |
627 |
|
|
C |
628 |
|
|
C TEST FOR RANK DEFICIENCY. |
629 |
|
|
C |
630 |
|
|
190 IF (DSQRT(DP/A1).GT.EPS1) GO TO 200 |
631 |
|
|
IF (IN.EQ.6) GO TO 200 |
632 |
|
|
II=IP-1 |
633 |
|
|
IF (ERM) WRITE (6,1140) IRANK,EPS1,II,II |
634 |
|
|
IRANK=IP-1 |
635 |
|
|
ERM=.FALSE. |
636 |
|
|
GO TO 260 |
637 |
|
|
C |
638 |
|
|
C (EPS1) RANK IS DEFICIENT. |
639 |
|
|
C |
640 |
|
|
200 SP=DSQRT(DP) |
641 |
|
|
C |
642 |
|
|
C BEGIN FRONT ELIMINATION ON COLUMN IP. |
643 |
|
|
C |
644 |
|
|
C SP=SQROOT(SIGMA**2). |
645 |
|
|
C |
646 |
|
|
BP=1.D0/(DP+SP*ABS(A(IP,IP))) |
647 |
|
|
C |
648 |
|
|
C STORE BETA IN S(3N+1)-S(3N+L). |
649 |
|
|
C |
650 |
|
|
IF (IP.EQ.K) BP=0.D0 |
651 |
|
|
N3=K+2*N+IP |
652 |
|
|
R(N3)=BP |
653 |
|
|
UP=DSIGN(DBLE(SP)+ABS(A(IP,IP)),DBLE(A(IP,IP))) |
654 |
|
|
IF (IP.GE.K) GO TO 250 |
655 |
|
|
IPP1=IP+1 |
656 |
|
|
IF (IP.GE.N) GO TO 240 |
657 |
|
|
DO 230 J=IPP1,N |
658 |
|
|
SJ=0.D0 |
659 |
|
|
DO 210 I=IPP1,K |
660 |
|
|
210 SJ=SJ+A(I,J)*A(I,IP) |
661 |
|
|
SJ=SJ+UP*A(IP,J) |
662 |
|
|
SJ=BP*SJ |
663 |
|
|
C |
664 |
|
|
C SJ=YJ NOW |
665 |
|
|
C |
666 |
|
|
DO 220 I=IPP1,K |
667 |
|
|
220 A(I,J)=A(I,J)-A(I,IP)*SJ |
668 |
|
|
230 A(IP,J)=A(IP,J)-SJ*UP |
669 |
|
|
240 A(IP,IP)=-SIGN(SP,A(IP,IP)) |
670 |
|
|
C |
671 |
|
|
N4=K+3*N+IP |
672 |
|
|
R(N4)=UP |
673 |
|
|
250 CONTINUE |
674 |
|
|
IRANK=LM |
675 |
|
|
260 IRP1=IRANK+1 |
676 |
|
|
IRM1=IRANK-1 |
677 |
|
|
IF (IRANK.EQ.0.OR.IRANK.EQ.N) GO TO 360 |
678 |
|
|
IF (IEQ.EQ.3) GO TO 290 |
679 |
|
|
C |
680 |
|
|
C BEGIN BACK PROCESSING FOR RANK DEFICIENCY CASE |
681 |
|
|
C IF IRANK IS LESS THAN N. |
682 |
|
|
C |
683 |
|
|
DO 280 J=1,N |
684 |
|
|
N2=J2+J |
685 |
|
|
N7=J7+J |
686 |
|
|
L=MIN0(J,IRANK) |
687 |
|
|
C |
688 |
|
|
C UNSCALE COLUMNS FOR RANK DEFICIENT MATRICES WHEN IEQ.NE.3. |
689 |
|
|
C |
690 |
|
|
DO 270 I=1,L |
691 |
|
|
270 A(I,J)=A(I,J)/S(N7) |
692 |
|
|
S(N7)=1.E0 |
693 |
|
|
280 S(N2)=1.E0 |
694 |
|
|
290 IP=IRANK |
695 |
|
|
300 SJ=0.D0 |
696 |
|
|
DO 310 J=IRP1,N |
697 |
|
|
SJ=SJ+A(IP,J)**2 |
698 |
|
|
310 CONTINUE |
699 |
|
|
SJ=SJ+A(IP,IP)**2 |
700 |
|
|
AJ=DSQRT(SJ) |
701 |
|
|
UP=DSIGN(AJ+ABS(A(IP,IP)),DBLE(A(IP,IP))) |
702 |
|
|
C |
703 |
|
|
C IP TH ELEMENT OF U VECTOR CALCULATED. |
704 |
|
|
C |
705 |
|
|
BP=1.D0/(SJ+ABS(A(IP,IP))*AJ) |
706 |
|
|
C |
707 |
|
|
C BP = 2/LENGTH OF U SQUARED. |
708 |
|
|
C |
709 |
|
|
IPM1=IP-1 |
710 |
|
|
IF (IPM1.LE.0) GO TO 340 |
711 |
|
|
DO 330 I=1,IPM1 |
712 |
|
|
DP=A(I,IP)*UP |
713 |
|
|
DO 320 J=IRP1,N |
714 |
|
|
DP=DP+A(I,J)*A(IP,J) |
715 |
|
|
320 CONTINUE |
716 |
|
|
DP=DP/(SJ+ABS(A(IP,IP))*AJ) |
717 |
|
|
C |
718 |
|
|
C CALC. (AJ,U), WHERE AJ=JTH ROW OF A |
719 |
|
|
C |
720 |
|
|
A(I,IP)=A(I,IP)-UP*DP |
721 |
|
|
C |
722 |
|
|
C MODIFY ARRAY A. |
723 |
|
|
C |
724 |
|
|
DO 330 J=IRP1,N |
725 |
|
|
330 A(I,J)=A(I,J)-A(IP,J)*DP |
726 |
|
|
340 A(IP,IP)=-DSIGN(AJ,DBLE(A(IP,IP))) |
727 |
|
|
C |
728 |
|
|
C CALC. MODIFIED PIVOT. |
729 |
|
|
C |
730 |
|
|
C |
731 |
|
|
C SAVE BETA AND IP TH ELEMENT OF U VECTOR IN R ARRAY. |
732 |
|
|
C |
733 |
|
|
N6=K+IP |
734 |
|
|
N7=K+N+IP |
735 |
|
|
R(N6)=BP |
736 |
|
|
R(N7)=UP |
737 |
|
|
C |
738 |
|
|
C TEST FOR END OF BACK PROCESSING. |
739 |
|
|
C |
740 |
|
|
IF (IP-1) 360,360,350 |
741 |
|
|
350 IP=IP-1 |
742 |
|
|
GO TO 300 |
743 |
|
|
360 IF (IN.EQ.6) RETURN |
744 |
|
|
DO 370 J=1,K |
745 |
|
|
370 R(J)=B(J) |
746 |
|
|
IT=0 |
747 |
|
|
C |
748 |
|
|
C SET INITIAL X VECTOR TO ZERO. |
749 |
|
|
C |
750 |
|
|
DO 380 J=1,N |
751 |
|
|
380 X(J)=0.D0 |
752 |
|
|
IF (IRANK.EQ.0) GO TO 690 |
753 |
|
|
C |
754 |
|
|
C APPLY Q TO RT. HAND SIDE. |
755 |
|
|
C |
756 |
|
|
390 DO 430 IP=1,IRANK |
757 |
|
|
N4=K+3*N+IP |
758 |
|
|
SJ=R(N4)*R(IP) |
759 |
|
|
IPP1=IP+1 |
760 |
|
|
IF (IPP1.GT.K) GO TO 410 |
761 |
|
|
DO 400 I=IPP1,K |
762 |
|
|
400 SJ=SJ+A(I,IP)*R(I) |
763 |
|
|
410 N3=K+2*N+IP |
764 |
|
|
BP=R(N3) |
765 |
|
|
IF (IPP1.GT.K) GO TO 430 |
766 |
|
|
DO 420 I=IPP1,K |
767 |
|
|
420 R(I)=R(I)-BP*A(I,IP)*SJ |
768 |
|
|
430 R(IP)=R(IP)-BP*R(N4)*SJ |
769 |
|
|
DO 440 J=1,IRANK |
770 |
|
|
440 S(J)=R(J) |
771 |
|
|
ENORM=0.E0 |
772 |
|
|
IF (IRP1.GT.K) GO TO 510 |
773 |
|
|
DO 450 J=IRP1,K |
774 |
|
|
450 ENORM=ENORM+R(J)**2 |
775 |
|
|
ENORM=SQRT(ENORM) |
776 |
|
|
GO TO 510 |
777 |
|
|
460 DO 480 J=1,N |
778 |
|
|
SJ=0.D0 |
779 |
|
|
N1=J1+J |
780 |
|
|
IP=S(N1) |
781 |
|
|
DO 470 I=1,K |
782 |
|
|
470 SJ=SJ+R(I)*AA(I,IP) |
783 |
|
|
C |
784 |
|
|
C APPLY AT TO RT. HAND SIDE. |
785 |
|
|
C APPLY SCALING. |
786 |
|
|
C |
787 |
|
|
N7=J2+IP |
788 |
|
|
N1=K+N+J |
789 |
|
|
480 R(N1)=SJ*S(N7) |
790 |
|
|
N1=K+N |
791 |
|
|
S(1)=R(N1+1)/A(1,1) |
792 |
|
|
IF (N.EQ.1) GO TO 510 |
793 |
|
|
DO 500 J=2,N |
794 |
|
|
N1=J-1 |
795 |
|
|
SJ=0.D0 |
796 |
|
|
DO 490 I=1,N1 |
797 |
|
|
490 SJ=SJ+A(I,J)*S(I) |
798 |
|
|
N2=K+J+N |
799 |
|
|
500 S(J)=(R(N2)-SJ)/A(J,J) |
800 |
|
|
C |
801 |
|
|
C ENTRY TO CONTINUE ITERATING. SOLVES TZ = C = 1ST IRANK |
802 |
|
|
C COMPONENTS OF QB . |
803 |
|
|
C |
804 |
|
|
510 S(IRANK)=S(IRANK)/A(IRANK,IRANK) |
805 |
|
|
IF (IRM1.EQ.0) GO TO 540 |
806 |
|
|
DO 530 J=1,IRM1 |
807 |
|
|
N1=IRANK-J |
808 |
|
|
N2=N1+1 |
809 |
|
|
SJ=0. |
810 |
|
|
DO 520 I=N2,IRANK |
811 |
|
|
520 SJ=SJ+A(N1,I)*S(I) |
812 |
|
|
530 S(N1)=(S(N1)-SJ)/A(N1,N1) |
813 |
|
|
C |
814 |
|
|
C Z CALCULATED. COMPUTE X = SZ. |
815 |
|
|
C |
816 |
|
|
540 IF (IRANK.EQ.N) GO TO 590 |
817 |
|
|
DO 550 J=IRP1,N |
818 |
|
|
550 S(J)=0.E0 |
819 |
|
|
DO 580 I=1,IRANK |
820 |
|
|
N7=K+N+I |
821 |
|
|
SJ=R(N7)*S(I) |
822 |
|
|
DO 560 J=IRP1,N |
823 |
|
|
SJ=SJ+A(I,J)*S(J) |
824 |
|
|
560 CONTINUE |
825 |
|
|
N6=K+I |
826 |
|
|
DO 570 J=IRP1,N |
827 |
|
|
570 S(J)=S(J)-A(I,J)*R(N6)*SJ |
828 |
|
|
580 S(I)=S(I)-R(N6)*R(N7)*SJ |
829 |
|
|
C |
830 |
|
|
C INCREMENT FOR X OF MINIMAL LENGTH CALCULATED. |
831 |
|
|
C |
832 |
|
|
590 DO 600 I=1,N |
833 |
|
|
600 X(I)=X(I)+S(I) |
834 |
|
|
IF (IN.EQ.7) GO TO 750 |
835 |
|
|
C |
836 |
|
|
C CALC. SUP NORM OF INCREMENT AND RESIDUALS |
837 |
|
|
C |
838 |
|
|
TOP1=0.E0 |
839 |
|
|
DO 610 J=1,N |
840 |
|
|
N2=J7+J |
841 |
|
|
610 TOP1=AMAX1(TOP1,ABS(S(J))*S(N2)) |
842 |
|
|
DO 630 I=1,K |
843 |
|
|
SJ=0.D0 |
844 |
|
|
DO 620 J=1,N |
845 |
|
|
N1=J1+J |
846 |
|
|
IP=S(N1) |
847 |
|
|
N7=J2+IP |
848 |
|
|
620 SJ=SJ+AA(I,IP)*X(J)*S(N7) |
849 |
|
|
630 R(I)=B(I)-SJ |
850 |
|
|
IF (ITMAX.LE.0) GO TO 750 |
851 |
|
|
C |
852 |
|
|
C CALC. SUP NORM OF X. |
853 |
|
|
C |
854 |
|
|
TOP=0.E0 |
855 |
|
|
DO 640 J=1,N |
856 |
|
|
N2=J7+J |
857 |
|
|
640 TOP=AMAX1(TOP,ABS(X(J))*S(N2)) |
858 |
|
|
C |
859 |
|
|
C COMPARE RELATIVE CHANGE IN X WITH TOLERANCE EPS . |
860 |
|
|
C |
861 |
|
|
IF (TOP1-TOP*EPS2) 690,650,650 |
862 |
|
|
650 IF (IT-ITMAX) 660,680,680 |
863 |
|
|
660 IT=IT+1 |
864 |
|
|
IF (IT.EQ.1) GO TO 670 |
865 |
|
|
IF (TOP1.GT..25*TOP2) GO TO 690 |
866 |
|
|
670 TOP2=TOP1 |
867 |
|
|
GO TO (390,460), ISW |
868 |
|
|
680 IT=0 |
869 |
|
|
690 SJ=0.D0 |
870 |
|
|
DO 700 J=1,K |
871 |
|
|
SJ=SJ+R(J)**2 |
872 |
|
|
700 CONTINUE |
873 |
|
|
ENORM=DSQRT(SJ) |
874 |
|
|
IF (IRANK.EQ.N.AND.ISW.EQ.1) GO TO 710 |
875 |
|
|
GO TO 730 |
876 |
|
|
710 ENM1=ENORM |
877 |
|
|
C |
878 |
|
|
C SAVE X ARRAY. |
879 |
|
|
C |
880 |
|
|
DO 720 J=1,N |
881 |
|
|
N1=K+J |
882 |
|
|
720 R(N1)=X(J) |
883 |
|
|
ISW=2 |
884 |
|
|
IT=0 |
885 |
|
|
GO TO 460 |
886 |
|
|
C |
887 |
|
|
C CHOOSE BEST SOLUTION |
888 |
|
|
C |
889 |
|
|
730 IF (IRANK.LT.N) GO TO 750 |
890 |
|
|
IF (ENORM.LE.ENM1) GO TO 750 |
891 |
|
|
DO 740 J=1,N |
892 |
|
|
N1=K+J |
893 |
|
|
740 X(J)=R(N1) |
894 |
|
|
ENORM=ENM1 |
895 |
|
|
C |
896 |
|
|
C NORM OF AX - B LOCATED IN THE CELL ENORM . |
897 |
|
|
C |
898 |
|
|
C |
899 |
|
|
C REARRANGE VARIABLES. |
900 |
|
|
C |
901 |
|
|
750 DO 760 J=1,N |
902 |
|
|
N1=J1+J |
903 |
|
|
760 S(J)=S(N1) |
904 |
|
|
DO 790 J=1,N |
905 |
|
|
DO 770 I=J,N |
906 |
|
|
IP=S(I) |
907 |
|
|
IF (J.EQ.IP) GO TO 780 |
908 |
|
|
770 CONTINUE |
909 |
|
|
780 S(I)=S(J) |
910 |
|
|
S(J)=J |
911 |
|
|
SJ=X(J) |
912 |
|
|
X(J)=X(I) |
913 |
|
|
790 X(I)=SJ |
914 |
|
|
C |
915 |
|
|
C SCALE VARIABLES. |
916 |
|
|
C |
917 |
|
|
DO 800 J=1,N |
918 |
|
|
N2=J2+J |
919 |
|
|
800 X(J)=X(J)*S(N2) |
920 |
|
|
RETURN |
921 |
|
|
C |
922 |
|
|
C RESTORE A. |
923 |
|
|
C |
924 |
|
|
810 DO 820 J=1,N |
925 |
|
|
N2=J2+J |
926 |
|
|
DO 820 I=1,K |
927 |
|
|
820 A(I,J)=AA(I,J) |
928 |
|
|
RETURN |
929 |
|
|
C |
930 |
|
|
C GENERATE SOLUTIONS TO THE HOMOGENEOUS EQUATION AX = 0. |
931 |
|
|
C |
932 |
|
|
830 IF (IRANK.EQ.N) RETURN |
933 |
|
|
NS=N-IRANK |
934 |
|
|
DO 840 I=1,N |
935 |
|
|
DO 840 J=1,NS |
936 |
|
|
840 H(I,J)=0.E0 |
937 |
|
|
DO 850 J=1,NS |
938 |
|
|
N2=IRANK+J |
939 |
|
|
850 H(N2,J)=1.E0 |
940 |
|
|
IF (IRANK.EQ.0) RETURN |
941 |
|
|
DO 870 J=1,IRANK |
942 |
|
|
DO 870 I=1,NS |
943 |
|
|
N7=K+N+J |
944 |
|
|
SJ=R(N7)*H(J,I) |
945 |
|
|
DO 860 K1=IRP1,N |
946 |
|
|
860 SJ=SJ+H(K1,I)*A(J,K1) |
947 |
|
|
N6=K+J |
948 |
|
|
BP=R(N6) |
949 |
|
|
DP=BP*R(N7)*SJ |
950 |
|
|
A1=DP |
951 |
|
|
A2=DP-A1 |
952 |
|
|
H(J,I)=H(J,I)-(A1+2.*A2) |
953 |
|
|
DO 870 K1=IRP1,N |
954 |
|
|
DP=BP*A(J,K1)*SJ |
955 |
|
|
A1=DP |
956 |
|
|
A2=DP-A1 |
957 |
|
|
870 H(K1,I)=H(K1,I)-(A1+2.*A2) |
958 |
|
|
C |
959 |
|
|
C REARRANGE ROWS OF SOLUTION MATRIX. |
960 |
|
|
C |
961 |
|
|
DO 880 J=1,N |
962 |
|
|
N1=J1+J |
963 |
|
|
880 S(J)=S(N1) |
964 |
|
|
DO 910 J=1,N |
965 |
|
|
DO 890 I=J,N |
966 |
|
|
IP=S(I) |
967 |
|
|
IF (J.EQ.IP) GO TO 900 |
968 |
|
|
890 CONTINUE |
969 |
|
|
900 S(I)=S(J) |
970 |
|
|
S(J)=J |
971 |
|
|
DO 910 K1=1,NS |
972 |
|
|
A1=H(J,K1) |
973 |
|
|
H(J,K1)=H(I,K1) |
974 |
|
|
910 H(I,K1)=A1 |
975 |
|
|
RETURN |
976 |
|
|
C |
977 |
|
|
1140 FORMAT (31H0WARNING. IRANK HAS BEEN SET TO,I4,6H BUT(,1PE10.3,9H) |
978 |
|
|
1 RANK IS,I4,25H. IRANK IS NOW TAKEN AS ,I4,1H.) |
979 |
|
|
END |
980 |
|
|
FUNCTION POTEM(T,S,P) |
981 |
|
|
implicit none |
982 |
|
|
C POTENTIAL TEMPERATURE FUNCTION |
983 |
|
|
C BASED ON FOFONOFF AND FROESE (1958) AS SHOWN IN "THE SEA" VOL. 1, |
984 |
|
|
C PAGE 17, TABLE IV |
985 |
|
|
C INPUT IS TEMPERATURE, SALINITY, PRESSURE (OR DEPTH) |
986 |
|
|
C UNITS ARE DEG.C., PPT, DBARS (OR METERS) |
987 |
|
|
real POTEM,T,S,P |
988 |
|
|
real B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11 |
989 |
|
|
real T2,T3,S2,P2 |
990 |
|
|
B1=-1.60E-5*P |
991 |
|
|
B2=1.014E-5*P*T |
992 |
|
|
T2=T*T |
993 |
|
|
T3=T2*T |
994 |
|
|
B3=-1.27E-7*P*T2 |
995 |
|
|
B4=2.7E-9*P*T3 |
996 |
|
|
B5=1.322E-6*P*S |
997 |
|
|
B6=-2.62E-8*P*S*T |
998 |
|
|
S2=S*S |
999 |
|
|
P2=P*P |
1000 |
|
|
B7=4.1E-9*P*S2 |
1001 |
|
|
B8=9.14E-9*P2 |
1002 |
|
|
B9=-2.77E-10*P2*T |
1003 |
|
|
B10=9.5E-13*P2*T2 |
1004 |
|
|
B11=-1.557E-13*P2*P |
1005 |
|
|
POTEM=B1+B2+B3+B4+B5+B6+B7+B8+B9+B10+B11 |
1006 |
|
|
POTEM=T-POTEM |
1007 |
|
|
RETURN |
1008 |
|
|
END |
1009 |
|
|
FUNCTION DN(T,S,D) |
1010 |
|
|
implicit none |
1011 |
|
|
real T,S,D |
1012 |
|
|
DOUBLE PRECISION DN,T3,S2,T2,S3,F1,F2,F3,FS,SIGMA,A,B1,B2,B,CO, |
1013 |
|
|
1ALPHA,ALPSTD |
1014 |
|
|
T2 = T*T |
1015 |
|
|
T3= T2* T |
1016 |
|
|
S2 = S*S |
1017 |
|
|
S3 = S2 * S |
1018 |
|
|
F1 = -(T-3.98)**2 * (T+283.)/(503.57*(T+67.26)) |
1019 |
|
|
F2 = T3*1.0843E-6 - T2*9.8185E-5 + T*4.786E-3 |
1020 |
|
|
F3 = T3*1.667E-8 - T2*8.164E-7 + T*1.803E-5 |
1021 |
|
|
FS= S3*6.76786136D-6 - S2*4.8249614D-4 + S*8.14876577D-1 |
1022 |
|
|
SIGMA= F1 + (FS+3.895414D-2)*(1.-F2+F3*(FS-.22584586D0)) |
1023 |
|
|
A= D*1.0E-4*(105.5+ T*9.50 - T2*0.158 - D*T*1.5E-4) - |
1024 |
|
|
1(227. + T*28.33 - T2*0.551 + T3* 0.004) |
1025 |
|
|
B1 = (FS-28.1324)/10.0 |
1026 |
|
|
B2 = B1 * B1 |
1027 |
|
|
B= -B1* (147.3-T*2.72 + T2*0.04 - D*1.0E-4*(32.4- 0.87*T+.02*T2)) |
1028 |
|
|
B= B+ B2*(4.5-0.1*T - D*1.0E-4*(1.8-0.06*T)) |
1029 |
|
|
CO = 4886./(1. + 1.83E-5*D) |
1030 |
|
|
ALPHA= D*1.0E-6* (CO+A+B) |
1031 |
|
|
DN=(SIGMA+ALPHA)/(1.-1.E-3*ALPHA) |
1032 |
|
|
RETURN |
1033 |
|
|
END |