/[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.2 - (show annotations) (download)
Mon Jun 22 15:26:26 1998 UTC (25 years, 11 months ago) by adcroft
Branch: MAIN
CVS Tags: checkpoint44e_post, release1_p13_pre, checkpoint44f_post, checkpoint43a-release1mods, release1_p13, checkpoint40pre3, checkpoint40pre1, checkpoint40pre7, checkpoint40pre6, checkpoint40pre9, checkpoint40pre8, chkpt44d_post, release1_p8, release1_p9, release1_p1, release1_p2, release1_p3, release1_p4, release1_p5, release1_p6, release1_p7, checkpoint44e_pre, release1_b1, checkpoint43, checkpoint38, release1_chkpt44d_post, checkpoint11, checkpoint10, checkpoint13, checkpoint15, checkpoint14, checkpoint17, checkpoint19, checkpoint18, branch-atmos-merge-shapiro, icebear5, icebear4, icebear3, icebear2, checkpoint40pre2, release1-branch_tutorials, checkpoint45d_post, chkpt44a_post, checkpoint44h_pre, checkpoint40pre4, checkpoint46a_post, checkpoint28, checkpoint29, checkpoint20, checkpoint21, checkpoint22, checkpoint23, checkpoint24, checkpoint25, checkpoint27, chkpt44c_pre, checkpoint45a_post, branch-atmos-merge-freeze, branch-atmos-merge-start, ecco_c44_e19, ecco_c44_e18, ecco_c44_e17, ecco_c44_e16, checkpoint9, checkpoint8, release1_p12, release1_p10, release1_p11, release1_p16, release1_p17, release1_p14, release1_p15, pre38tag1, checkpoint44g_post, checkpoint26, checkpoint45b_post, checkpoint46b_pre, release1-branch-end, checkpoint12, c37_adj, release1_final_v1, checkpoint16, checkpoint46, checkpoint44b_post, checkpoint46a_pre, checkpoint45c_post, ecco_ice2, ecco_ice1, checkpoint44h_post, pre38-close, release1_p12_pre, checkpoint39, checkpoint33, checkpoint32, checkpoint31, checkpoint30, checkpoint37, checkpoint36, checkpoint35, checkpoint34, ecco_c44_e22, ecco_c44_e25, checkpoint40pre5, chkpt44a_pre, ecco_c44_e23, ecco_c44_e20, ecco_c44_e21, ecco_c44_e26, ecco_c44_e27, ecco_c44_e24, ecco-branch-mod1, ecco-branch-mod2, ecco-branch-mod3, ecco-branch-mod4, ecco-branch-mod5, branch-atmos-merge-zonalfilt, release1_beta1, checkpoint44b_pre, checkpoint42, checkpoint40, checkpoint41, checkpoint44, checkpoint45, chkpt44c_post, branch-point-rdot, checkpoint44f_pre, branch-atmos-merge-phase5, branch-atmos-merge-phase4, branch-atmos-merge-phase7, branch-atmos-merge-phase6, branch-atmos-merge-phase1, branch-atmos-merge-phase3, branch-atmos-merge-phase2, release1-branch_branchpoint
Branch point for: c24_e25_ice, release1_final, release1-branch, release1, ecco-branch, branch-rdot, release1_50yr, icebear, branch-atmos-merge, pre38, release1_coupled
Changes since 1.1: +6 -7 lines 
Various changes including time-dependant forcing:
 o logic for controlling external forcing fields now allows
   for time-dependant forcing: load_external_fields.F
 o genmake.dec needed a special line for the above file.
 o theta* and salt* time-stepping algorithm were re-implemented.
The 4x4 global configuration has been "double-checked" against
CNH's version. However, we do not assume any responsibility for
the correctness of this code ...  8-)
1 C $Header: /u/gcmpack/models/MITgcmUV/utils/knudsen2/knudsen2.f,v 1.1 1998/05/28 16:26:56 adcroft Exp $
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=15)
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 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
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
  ViewVC Help
Powered by ViewVC 1.1.22