| 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 |
Parent Directory
| 
Revision Log
| 
Revision Graph