/[MITgcm]/MITgcm/utils/knudsen2/knudsen2.f
ViewVC logotype

Contents of /MITgcm/utils/knudsen2/knudsen2.f

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (show annotations) (download)
Thu May 28 16:26:56 1998 UTC (24 years, 8 months ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint5, checkpoint4, checkpoint7, checkpoint6
Branch point for: checkpoint7-4degree-ref
knudsen2.f is an auxilliary program for calculating the polynomial
coefficients for the POL3 appriximation to the EOS.
Currently creates two files:
 o POLY3.COEFFS is read by MITgcmUV
 o polyeos.coeffs is read by compare01

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

  ViewVC Help
Powered by ViewVC 1.1.22