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

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

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


Revision 1.1 - (show annotations) (download)
Wed Aug 7 16:55:52 2002 UTC (20 years, 1 month ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint46n_post, checkpoint51k_post, checkpoint62v, checkpoint47e_post, checkpoint57m_post, checkpoint52l_pre, checkpoint62u, hrcube4, hrcube5, checkpoint46l_post, checkpoint57g_pre, checkpoint46g_pre, checkpoint47c_post, checkpoint62t, checkpoint50c_post, checkpoint57s_post, checkpoint58b_post, checkpoint57b_post, checkpoint46f_post, checkpoint52d_pre, checkpoint57g_post, checkpoint48e_post, checkpoint56b_post, checkpoint50c_pre, checkpoint57y_post, checkpoint46b_post, checkpoint52j_pre, checkpoint51o_pre, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint54d_post, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint54e_post, checkpoint65l, checkpoint65m, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, checkpoint62c, checkpoint51l_post, checkpoint48i_post, checkpoint57r_post, checkpoint46l_pre, checkpoint57d_post, checkpoint57i_post, checkpoint50d_pre, checkpoint52k_post, checkpoint59, checkpoint58, checkpoint55, checkpoint54, checkpoint57, checkpoint56, checkpoint51, checkpoint50, checkpoint53, checkpoint52, checkpoint50d_post, checkpoint58f_post, checkpoint52f_post, checkpoint57n_post, checkpoint58d_post, checkpoint62s, checkpoint58a_post, checkpoint62r, checkpoint62q, checkpoint50b_pre, checkpoint62p, checkpoint57z_post, checkpoint54f_post, checkpoint51f_post, checkpoint62a, checkpoint62g, checkpoint62f, checkpoint62e, checkpoint62d, checkpoint62k, checkpoint62j, checkpoint62i, checkpoint62h, checkpoint62o, checkpoint62n, checkpoint62m, checkpoint62l, checkpoint62w, checkpoint62z, checkpoint62y, checkpoint62x, checkpoint58y_post, checkpoint48b_post, checkpoint51d_post, checkpoint48c_pre, checkpoint47d_pre, checkpoint51t_post, checkpoint58t_post, checkpoint51n_post, checkpoint55i_post, checkpoint58m_post, checkpoint57l_post, checkpoint52i_pre, hrcube_1, hrcube_2, hrcube_3, checkpoint51s_post, checkpoint47a_post, checkpoint57t_post, checkpoint55c_post, checkpoint48d_pre, checkpoint51j_post, checkpoint47i_post, checkpoint63g, checkpoint52e_pre, checkpoint57v_post, checkpoint57f_post, checkpoint52e_post, checkpoint51n_pre, checkpoint47d_post, checkpoint53d_post, checkpoint46d_pre, checkpoint64, checkpoint65, checkpoint60, checkpoint61, checkpoint62, checkpoint63, checkpoint57a_post, checkpoint48d_post, checkpoint57h_pre, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint48f_post, checkpoint52b_pre, checkpoint54b_post, checkpoint46j_pre, checkpoint58w_post, checkpoint57h_post, checkpoint51l_pre, checkpoint52m_post, checkpoint57y_pre, checkpoint55g_post, checkpoint48h_post, checkpoint51q_post, checkpoint51b_pre, checkpoint47g_post, checkpoint52b_post, checkpoint52c_post, checkpoint46j_post, checkpoint51h_pre, checkpoint46k_post, checkpoint58o_post, checkpoint48a_post, checkpoint57c_post, checkpoint50f_post, checkpoint50a_post, checkpoint50f_pre, checkpoint58p_post, checkpoint58q_post, checkpoint52f_pre, checkpoint55d_post, checkpoint58e_post, checkpoint47j_post, checkpoint54a_pre, checkpoint63p, checkpoint63q, checkpoint63r, checkpoint63s, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint63h, checkpoint63i, checkpoint63j, checkpoint63k, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint53c_post, checkpoint63a, checkpoint63b, checkpoint63c, checkpoint55d_pre, checkpoint57c_pre, checkpoint58r_post, checkpoint55j_post, branch-exfmods-tag, branchpoint-genmake2, checkpoint54a_post, checkpoint46e_pre, checkpoint55h_post, checkpoint58n_post, checkpoint51r_post, checkpoint48c_post, checkpoint51i_post, checkpoint57e_post, checkpoint55b_post, checkpoint51b_post, checkpoint51c_post, checkpoint46c_pre, checkpoint53a_post, checkpoint65o, checkpoint47b_post, checkpoint59q, checkpoint59p, checkpoint59r, checkpoint59e, checkpoint59d, checkpoint59g, checkpoint59f, checkpoint59a, checkpoint55f_post, checkpoint59c, checkpoint59b, checkpoint59m, checkpoint59l, checkpoint59o, checkpoint59n, checkpoint59i, checkpoint59h, checkpoint59k, checkpoint46h_pre, checkpoint52d_post, checkpoint53g_post, checkpoint46m_post, checkpoint57p_post, checkpint57u_post, checkpoint50g_post, checkpoint57q_post, eckpoint57e_pre, checkpoint46g_post, checkpoint58k_post, checkpoint52a_pre, checkpoint62b, checkpoint58v_post, checkpoint50h_post, checkpoint52i_post, checkpoint50e_pre, checkpoint50i_post, checkpoint51i_pre, checkpoint52h_pre, checkpoint56a_post, checkpoint64y, checkpoint64x, checkpoint58l_post, checkpoint64z, checkpoint53f_post, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64w, checkpoint64v, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n, checkpoint64a, checkpoint64c, checkpoint64b, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f, checkpoint57h_done, checkpoint52j_post, checkpoint47f_post, checkpoint50e_post, checkpoint46i_post, checkpoint57j_post, checkpoint57f_pre, checkpoint61f, checkpoint46c_post, checkpoint58g_post, branch-netcdf, checkpoint52l_post, checkpoint58x_post, checkpoint61n, checkpoint52n_post, checkpoint53b_pre, checkpoint46e_post, checkpoint59j, checkpoint58h_post, checkpoint56c_post, checkpoint58j_post, checkpoint51e_post, checkpoint57a_pre, checkpoint55a_post, checkpoint47, checkpoint48, checkpoint49, checkpoint57o_post, checkpoint46h_post, checkpoint51o_post, checkpoint61q, checkpoint57k_post, checkpoint51f_pre, checkpoint48g_post, checkpoint53b_post, checkpoint47h_post, checkpoint52a_post, checkpoint57w_post, checkpoint61e, checkpoint58i_post, checkpoint51g_post, ecco_c52_e35, checkpoint57x_post, checkpoint46d_post, checkpoint50b_post, checkpoint58c_post, checkpoint58u_post, checkpoint51m_post, checkpoint53d_pre, checkpoint58s_post, checkpoint55e_post, checkpoint61g, checkpoint61d, checkpoint54c_post, checkpoint61b, checkpoint61c, checkpoint61a, checkpoint51a_post, checkpoint61o, checkpoint61l, checkpoint61m, checkpoint61j, checkpoint61k, checkpoint61h, checkpoint61i, checkpoint61v, checkpoint61w, checkpoint61t, checkpoint61u, checkpoint61r, checkpoint61s, checkpoint61p, checkpoint51p_post, checkpoint61z, checkpoint61x, checkpoint61y, checkpoint51u_post, HEAD
Branch point for: branch-exfmods-curt, branch-genmake2, branch-nonh, tg2-branch, ecco-branch, netcdf-sm0, checkpoint51n_branch
o Added new equation of state -> JMD95Z and JMD95P
  - EOS of Jackett and McDougall, 1995, JPO
  - moved all EOS parameters into EOS.h
  - new routines ini_eos.F, store_pressure.F
o Added UNESCO EOS, but not recommended because it requires
  in-situ temperature (see JMD95)
o Modified formatting for knudsen2.f in utils/knudsen2 and added
  unesco.f to be used with POLY3

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

  ViewVC Help
Powered by ViewVC 1.1.22