File indexing completed on 2018-03-02 18:45:14 UTC
view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
a37a13034c Mart*0001 PROGRAM KNUDSEN2
0002 implicit none
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012 integer NumLevels
0013 PARAMETER (NumLevels=15)
0014
0015
0016 integer NumTempPt,NumSalPt,NumCallocPt
0017 PARAMETER (NumTempPt = 10, NumSalPt = 5 )
0018 PARAMETER (NumCallocPt =NumSalPt*NumTempPt )
0019
0020
0021 integer NTerms
0022 PARAMETER (NTerms = 9 )
0023
0024
0025 DOUBLE PRECISION DN
0026 DOUBLE PRECISION DUNESCO
0027 real POTEM
0028
0029
0030 real A(NumCallocPt,NumCallocPt),
0031 & B(NumCallocPt),
0032 & X(NTerms),
0033 & AB(13,NumLevels),
0034 & C(NumCallocPt,NTerms),
0035 & H(NTerms,NTerms),
0036 & R(NumCallocPt*2),SB(NumCallocPt*4)
0037
0038 real
0039 & Z(NumLevels),
0040 & DZ(NumLevels)
0041 integer
0042 & IDX(NumLevels)
0043
0044 real
0045 & Theta(NumCallocPt),
0046 & SPick(NumCallocPt),TPick(NumCallocPt),DenPick(NumCallocPt)
0047
0048
0049
0050
0051
0052
0053
0054 real TMin(25), TMax(25)
0055 & ,SMin(25), SMax(25)
0056
0057
0058 DATA TMin / 4*-2., 15*-1., 6*0. /
0059 DATA TMax / 29., 19., 14., 11., 9., 7., 5., 3*4., 5*3., 10*2. /
0060 DATA SMin / 28.5, 33.7, 34., 34.1, 34.2, 34.4, 2*34.5,
0061 & 15*34.6, 2*34.7 /
0062 DATA SMax / 36.7, 36.6, 35.8, 35.7, 35.3, 2*35.1, 7*35.,
0063 & 9*34.9, 2*34.8 /
0064
0065
0066 integer I,J,K,IT,IEQ,IRANK,NDIM,NHDIM,N,M,IN,ITMAX,MP,L,NN
0067 real T,S,D,EPS,DeltaT,DeltaS,ENORM,DeltaDen,DensityRef
0068 real Tbar,Sbar,Tsum,SSum,TempInc,SalnInc
0069 real DensitySum,ThetaSum
0070 real DensityBar,thetabar,TempRef,SAlRef
0071
0072
0073
0074
0075 DATA dz/ 5.00e+03,7.00e+03,1.00e+04,1.40e+04,
0076 &1.90e+04,2.40e+04,2.90e+04,3.40e+04,3.90e+04,
0077 &4.40e+04,4.90e+04,5.40e+04,5.90e+04,6.40e+04,
0078 &6.90e+04/
0079
0080
0081
0082 Z(1) = 0.
0083 Z(1)=.5*DZ(1)/100.
0084 DO K=2,NumLevels
0085 Z(K)=Z(K-1)+.5*(DZ(K)+DZ(K-1))/100.
0086 ENDDO
0087
0088
0089
0090 DO I=1,NumLevels
0091 IDX( I ) = I
0092
0093 IDX( I ) = IFIX(Z(I)/250.)+1
0094 ENDDO
0095
0096
0097
0098 PRINT 419
0099 PRINT 422,(Z(I),
0100 & Tmin( IDX(I)),Tmax( IDX(I)),
0101 & SMin(IDX(I)),Smax( IDX(I)),
0102
0103 & ( IDX(I) ),I=1,NumLevels )
0104
0105
0106
0107
0108 DO MP=1,NumLevels
0109
0110
0111 TempInc =(Tmax( IDX(MP))-TMin(IDX(MP)))/(2.*FLOAT(NumSalPt)-1.0)
0112 SalnInc =(Smax( IDX(MP))-SMin(IDX(MP)))/(FLOAT(NumSalPt)-1.0)
0113 DO I=1,NumTempPt
0114 DO J=1,NumSalPt
0115 K=NumSalPt*I+J-NumSalPt
0116 TPick(K)=TMin(IDX(MP))+(FLOAT(I)-1.0)*TempInc
0117 SPick(K)=SMin(IDX(MP))+(FLOAT(J)-1.0)*SalnInc
0118 ENDDO
0119 ENDDO
0120
0121
0122
0123
0124 Tsum=0.0
0125 Ssum=0.0
0126 DensitySum=0.0
0127 ThetaSum=0.0
0128 DO K=1,NumCallocPt
0129 D=Z(MP)
0130 S=SPick(K)
0131 T=TPick(K)
0132
0133 DenPick(K)=DUNESCO(T,S,D)
0134 Theta(K)=POTEM(T,S,D)
0135 Tsum=Tsum+TPick(K)
0136 Ssum=Ssum+SPick(K)
0137 DensitySum = DensitySum+DenPick(K)
0138 ThetaSum = ThetaSum+Theta(K)
0139 ENDDO
0140
0141
0142
0143
0144
0145 Tbar=Tsum / FLOAT( NumCallocPt )
0146 Sbar=Ssum / FLOAT( NumCallocPt )
0147 DensityBar=DensitySum / FLOAT( NumCallocPt )
0148
0149
0150 ThetaBar=ThetaSum/ FLOAT( NumCallocPt )
0151
0152
0153
0154 DensityRef=DUNESCO(Tbar,Sbar,D)
0155
0156 SalRef = Sbar
0157
0158 TempRef=TBar
0159 TempRef=ThetaBar
0160
0161
0162
0163
0164 AB(1,MP)=Z(MP)
0165 AB(2,MP)=DensityRef
0166 AB(3,MP)=TempRef
0167 AB(4,MP)=SalRef
0168 DO K=1,NumCallocPt
0169 TPick(K)=Theta(K)
0170 DeltaT = TPick(K) - TempRef
0171 DeltaS = SPick(K) - SalRef
0172 DeltaDen = DenPick(K) - DensityRef
0173 B(K)= DeltaDen
0174 A(K,1)=DeltaT
0175 A(K,2)=DeltaS
0176 A(K,3)=DeltaT*DeltaT
0177 A(K,4)=DeltaT*DeltaS
0178 A(K,5)=DeltaS*DeltaS
0179 A(K,6)=A(K,3)*DeltaT
0180 A(K,7)=A(K,4)*DeltaT
0181 A(K,8)=A(K,4)*DeltaS
0182 A(K,9)=A(K,5)*DeltaS
0183 ENDDO
0184
0185
0186 NDIM=NumCallocPt
0187
0188
0189 M=NumCallocPt
0190
0191
0192 N=NTerms
0193
0194 IN=1
0195
0196
0197 ITMAX=100
0198
0199 IT=0
0200 IEQ=2
0201 IRANK=0
0202 EPS=1.0E-7
0203 EPS=1.0E-11
0204 NHDIM=NTerms
0205 CALL LSQSL2(NDIM,A,M,N,B,X,IRANK,IN,ITMAX,IT,IEQ,ENORM,EPS,
0206 & NHDIM,H,C,R,SB)
0207
0208
0209 DO I=1,N
0210 AB(I+4,MP)=X(I)
0211 ENDDO
0212 ENDDO
0213 PRINT 430
0214 430 FORMAT(1X,' Z SIG0 T S X1 X2
0215 1 X3 X4 X5 X6 X7 X8
0216 2 X9 ',/)
0217 NN=N+4
0218
0219
0220
0221
0222
0223
0224
0225
0226
0227
0228
0229
0230
0231
0232
0233
0234
0235
0236
0237
0238
0239
0240
0241
0242
0243
0244
0245
0246
0247
0248
0249
0250
0251
0252
0253
0254
0255
0256
0257
0258
0259
0260 open(99,file='POLY3.COEFFS',form='formatted',status='unknown')
0261 write(99,*) NumLevels
0262 write(99,'(3(G20.13))') (ab(3,J),ab(4,J),ab(2,J),J=1 , NumLevels)
0263 do J=1,NumLevels
0264 write(99,'(3(G20.13))')(AB(I,J),I=5,13)
0265 enddo
0266 close(99)
0267
0268
0269
0270
0271
0272
0273
0274
0275
0276
0277
0278
0279
0280
0281
0282
0283
0284
0285
0286
0287
0288 open(99,file='polyeos.coeffs',form='formatted',status='unknown')
0289 do J=1,NumLevels
0290 write(99,*) ab(3,J)
0291 enddo
0292 do J=1,NumLevels
0293 write(99,*) ab(4,J)
0294 enddo
0295 do J=1,NumLevels
0296 do I=5,13
0297 write(99,*) AB(I,J)
0298 enddo
0299 enddo
0300 do J=1,NumLevels
0301 write(99,*) ab(2,J)
0302 enddo
0303
0304
0305
0306
0307
0308
0309
0310
0311
0312
0313
0314
0315
0316
0317
0318
0319
0320
0321
0322
0323
0324
0325
0326
0327
0328
0329
0330
0331
0332
0333
0334
0335
0336
0337
0338
0339
0340
0341
0342
0343
0344
0345
0346
0347
0348
0349
0350
0351
0352
0353
0354
0355
0356
0357
0358
0359
0360
0361
0362
0363
0364
0365
0366
0367
0368
0369
0370
0371
0372
0373
0374
0375
0376
0377
0378
0379
0380
0381
0382
0383 1298 FORMAT(18H DATA TRef /,67X,I9)
0384 1297 FORMAT(18H DATA SRef /,67X,I9)
0385
0386
0387
0388 419 FORMAT(5X,'LEVEL TMIN TMAX SMIN SMAX D/250')
0389 422 FORMAT (5X,F6.1,4F10.3,I3)
0390 412 FORMAT(1X,F6.1,F8.4,F7.3,F6.2,9E12.5)
0391 1296 FORMAT(6X,'DATA (C(',I2,',N),N=1,9)/',54X,I9)
0392 1295 FORMAT(5X,1H*,9X,4(E13.7,1H,),10X,I9)
0393 1294 FORMAT(5X,1H*,9X,E13.7,1H/,52X,I9)
0394 1280 FORMAT(5X,1H*,8X,5(F10.7,1H,),12X,I9)
0395 1281 FORMAT(5X,1H*,8X,F10.7,1H/,56X,I9)
0396 1282 FORMAT(5X,1H*,8X,F10.7,1H,,F10.7,1H/,45X,I9)
0397 1283 FORMAT(5X,1H*,8X,2(F10.7,1H,),F10.7,1H/,34X,I9)
0398 1284 FORMAT(5X,1H*,8X,3(F10.7,1H,),F10.7,1H/,23X,I9)
0399 1285 FORMAT(5X,1H*,8X,4(F10.7,1H,),F10.7,1H/,12X,I9)
0400 1286 FORMAT(5X,1H*,1H/)
0401 350 FORMAT(1X,E14.7)
0402 351 FORMAT(1H+,T86,E14.7/)
0403 400 FORMAT(1X,9E14.7)
0404 410 FORMAT(1X,5E14.7)
0405 411 FORMAT(///)
0406 STOP
0407 END
0408
0409
0410 SUBROUTINE LSQSL2 (NDIM,A,D,W,B,X,IRANK,IN,ITMAX,IT,IEQ,ENORM,EPS1
0411 1,NHDIM,H,AA,R,S)
0412 implicit none
0413
0414
0415
0416
0417
0418
0419
0420
0421
0422
0423
0424
0425
0426
0427
0428
0429
0430
0431
0432
0433
0434
0435
0436
0437 integer ndim,nhdim
0438 DOUBLE PRECISION SJ,DP,DP1,UP,BP,AJ
0439 LOGICAL ERM
0440 INTEGER D,W
0441
0442
0443
0444
0445
0446
0447
0448
0449
0450
0451
0452
0453
0454
0455
0456
0457
0458
0459
0460
0461
0462
0463 real A(NDIM,NDIM),B(1),AA(D,W),S(1), X(1),H(NHDIM,NHDIM),R(1)
0464
0465
0466
0467 integer K,N,IT,ISW,L,M,IRANK,IEQ,IN,K1
0468 integer J1,J2,J3,J4,J5,J6,J7,J8,J9
0469 integer N1,N2,N3,N4,N5,N6,N7,N8,NS
0470 integer I,ITMAX,IPM1,II,LM,J,IP,KM,IPP1,IRP1,IRM1
0471 real SP,ENORM,TOP,TOP1,ENM1,TOP2,EPS1,EPS2,A1,A2,AM
0472
0473 K=D
0474 N=W
0475 ERM=.TRUE.
0476
0477
0478
0479 IF (IT.EQ.0) ERM=.FALSE.
0480
0481
0482
0483
0484
0485
0486
0487
0488
0489
0490
0491
0492
0493 DATA EPS2/1.E-16/
0494
0495
0496
0497 ISW=1
0498 L=MIN0(K,N)
0499 M=MAX0(K,N)
0500 J1=M
0501 J2=N+J1
0502 J3=J2+N
0503 J4=J3+L
0504 J5=J4+L
0505 J6=J5+L
0506 J7=J6+L
0507 J8=J7+N
0508 J9=J8+N
0509 LM=L
0510 IF (IRANK.GE.1.AND.IRANK.LE.L) LM=IRANK
0511 IF (IN.EQ.6) LM=L
0512 IF (IN.EQ.8) RETURN
0513
0514
0515
0516 GO TO (10,360,810,390,830,10,10), IN
0517
0518
0519
0520
0521
0522 10 CONTINUE
0523
0524
0525
0526 IF (IN.GT.5) GO TO 30
0527 DO 20 J=1,N
0528 DO 20 I=1,K
0529 20 AA(I,J)=A(I,J)
0530 30 CONTINUE
0531 IF (IEQ.EQ.1) GO TO 60
0532 DO 50 J=1,N
0533 AM=0.E0
0534 DO 40 I=1,K
0535 40 AM=AMAX1(AM,ABS(A(I,J)))
0536
0537
0538
0539 N2=J2+J
0540 IF (IN.EQ.6) AM=1.E0
0541 S(N2)=1.E0/AM
0542 DO 50 I=1,K
0543 50 A(I,J)=A(I,J)*S(N2)
0544 GO TO 100
0545 60 AM=0.E0
0546 DO 70 J=1,N
0547 DO 70 I=1,K
0548 70 AM=AMAX1(AM,ABS(A(I,J)))
0549 AM=AM/FLOAT(K*N)
0550 IF (IN.EQ.6) AM=1.E0
0551 DO 80 J=1,N
0552 N2=J2+J
0553 80 S(N2)=1.E0/AM
0554 DO 90 J=1,N
0555 N2=J2+J
0556 DO 90 I=1,K
0557 90 A(I,J)=A(I,J)*S(N2)
0558
0559
0560
0561
0562 100 DO 110 J=1,N
0563 N7=J7+J
0564 N2=J2+J
0565 110 S(N7)=S(N2)
0566
0567
0568
0569
0570
0571 DO 120 J=1,N
0572 N1=J1+J
0573 120 S(N1)=J
0574
0575
0576
0577
0578
0579 DO 250 IP=1,LM
0580
0581
0582 DP=0.D0
0583 KM=IP
0584 DO 140 J=IP,N
0585 SJ=0.D0
0586 DO 130 I=IP,K
0587 SJ=SJ+A(I,J)**2
0588 130 CONTINUE
0589 IF (DP.GT.SJ) GO TO 140
0590 DP=SJ
0591 KM=J
0592 IF (IN.EQ.6) GO TO 160
0593 140 CONTINUE
0594
0595
0596
0597
0598
0599
0600
0601
0602 IF (KM.EQ.IP) GO TO 160
0603 DO 150 I=1,K
0604 A1=A(I,IP)
0605 A(I,IP)=A(I,KM)
0606 150 A(I,KM)=A1
0607
0608
0609
0610 N1=J1+KM
0611 A1=S(N1)
0612 N2=J1+IP
0613 S(N1)=S(N2)
0614 S(N2)=A1
0615 N7=J7+KM
0616 N8=J7+IP
0617 A1=S(N7)
0618 S(N7)=S(N8)
0619 S(N8)=A1
0620 160 IF (IP.EQ.1) GO TO 180
0621 A1=0.E0
0622 IPM1=IP-1
0623 DO 170 I=1,IPM1
0624 A1=A1+A(I,IP)**2
0625 170 CONTINUE
0626 IF (A1.GT.0.E0) GO TO 190
0627 180 IF (DP.GT.0.D0) GO TO 200
0628
0629
0630
0631 190 IF (DSQRT(DP/A1).GT.EPS1) GO TO 200
0632 IF (IN.EQ.6) GO TO 200
0633 II=IP-1
0634 IF (ERM) WRITE (6,1140) IRANK,EPS1,II,II
0635 IRANK=IP-1
0636 ERM=.FALSE.
0637 GO TO 260
0638
0639
0640
0641 200 SP=DSQRT(DP)
0642
0643
0644
0645
0646
0647 BP=1.D0/(DP+SP*ABS(A(IP,IP)))
0648
0649
0650
0651 IF (IP.EQ.K) BP=0.D0
0652 N3=K+2*N+IP
0653 R(N3)=BP
0654 UP=DSIGN(DBLE(SP)+ABS(A(IP,IP)),DBLE(A(IP,IP)))
0655 IF (IP.GE.K) GO TO 250
0656 IPP1=IP+1
0657 IF (IP.GE.N) GO TO 240
0658 DO 230 J=IPP1,N
0659 SJ=0.D0
0660 DO 210 I=IPP1,K
0661 210 SJ=SJ+A(I,J)*A(I,IP)
0662 SJ=SJ+UP*A(IP,J)
0663 SJ=BP*SJ
0664
0665
0666
0667 DO 220 I=IPP1,K
0668 220 A(I,J)=A(I,J)-A(I,IP)*SJ
0669 230 A(IP,J)=A(IP,J)-SJ*UP
0670 240 A(IP,IP)=-SIGN(SP,A(IP,IP))
0671
0672 N4=K+3*N+IP
0673 R(N4)=UP
0674 250 CONTINUE
0675 IRANK=LM
0676 260 IRP1=IRANK+1
0677 IRM1=IRANK-1
0678 IF (IRANK.EQ.0.OR.IRANK.EQ.N) GO TO 360
0679 IF (IEQ.EQ.3) GO TO 290
0680
0681
0682
0683
0684 DO 280 J=1,N
0685 N2=J2+J
0686 N7=J7+J
0687 L=MIN0(J,IRANK)
0688
0689
0690
0691 DO 270 I=1,L
0692 270 A(I,J)=A(I,J)/S(N7)
0693 S(N7)=1.E0
0694 280 S(N2)=1.E0
0695 290 IP=IRANK
0696 300 SJ=0.D0
0697 DO 310 J=IRP1,N
0698 SJ=SJ+A(IP,J)**2
0699 310 CONTINUE
0700 SJ=SJ+A(IP,IP)**2
0701 AJ=DSQRT(SJ)
0702 UP=DSIGN(AJ+ABS(A(IP,IP)),DBLE(A(IP,IP)))
0703
0704
0705
0706 BP=1.D0/(SJ+ABS(A(IP,IP))*AJ)
0707
0708
0709
0710 IPM1=IP-1
0711 IF (IPM1.LE.0) GO TO 340
0712 DO 330 I=1,IPM1
0713 DP=A(I,IP)*UP
0714 DO 320 J=IRP1,N
0715 DP=DP+A(I,J)*A(IP,J)
0716 320 CONTINUE
0717 DP=DP/(SJ+ABS(A(IP,IP))*AJ)
0718
0719
0720
0721 A(I,IP)=A(I,IP)-UP*DP
0722
0723
0724
0725 DO 330 J=IRP1,N
0726 330 A(I,J)=A(I,J)-A(IP,J)*DP
0727 340 A(IP,IP)=-DSIGN(AJ,DBLE(A(IP,IP)))
0728
0729
0730
0731
0732
0733
0734 N6=K+IP
0735 N7=K+N+IP
0736 R(N6)=BP
0737 R(N7)=UP
0738
0739
0740
0741 IF (IP-1) 360,360,350
0742 350 IP=IP-1
0743 GO TO 300
0744 360 IF (IN.EQ.6) RETURN
0745 DO 370 J=1,K
0746 370 R(J)=B(J)
0747 IT=0
0748
0749
0750
0751 DO 380 J=1,N
0752 380 X(J)=0.D0
0753 IF (IRANK.EQ.0) GO TO 690
0754
0755
0756
0757 390 DO 430 IP=1,IRANK
0758 N4=K+3*N+IP
0759 SJ=R(N4)*R(IP)
0760 IPP1=IP+1
0761 IF (IPP1.GT.K) GO TO 410
0762 DO 400 I=IPP1,K
0763 400 SJ=SJ+A(I,IP)*R(I)
0764 410 N3=K+2*N+IP
0765 BP=R(N3)
0766 IF (IPP1.GT.K) GO TO 430
0767 DO 420 I=IPP1,K
0768 420 R(I)=R(I)-BP*A(I,IP)*SJ
0769 430 R(IP)=R(IP)-BP*R(N4)*SJ
0770 DO 440 J=1,IRANK
0771 440 S(J)=R(J)
0772 ENORM=0.E0
0773 IF (IRP1.GT.K) GO TO 510
0774 DO 450 J=IRP1,K
0775 450 ENORM=ENORM+R(J)**2
0776 ENORM=SQRT(ENORM)
0777 GO TO 510
0778 460 DO 480 J=1,N
0779 SJ=0.D0
0780 N1=J1+J
0781 IP=S(N1)
0782 DO 470 I=1,K
0783 470 SJ=SJ+R(I)*AA(I,IP)
0784
0785
0786
0787
0788 N7=J2+IP
0789 N1=K+N+J
0790 480 R(N1)=SJ*S(N7)
0791 N1=K+N
0792 S(1)=R(N1+1)/A(1,1)
0793 IF (N.EQ.1) GO TO 510
0794 DO 500 J=2,N
0795 N1=J-1
0796 SJ=0.D0
0797 DO 490 I=1,N1
0798 490 SJ=SJ+A(I,J)*S(I)
0799 N2=K+J+N
0800 500 S(J)=(R(N2)-SJ)/A(J,J)
0801
0802
0803
0804
0805 510 S(IRANK)=S(IRANK)/A(IRANK,IRANK)
0806 IF (IRM1.EQ.0) GO TO 540
0807 DO 530 J=1,IRM1
0808 N1=IRANK-J
0809 N2=N1+1
0810 SJ=0.
0811 DO 520 I=N2,IRANK
0812 520 SJ=SJ+A(N1,I)*S(I)
0813 530 S(N1)=(S(N1)-SJ)/A(N1,N1)
0814
0815
0816
0817 540 IF (IRANK.EQ.N) GO TO 590
0818 DO 550 J=IRP1,N
0819 550 S(J)=0.E0
0820 DO 580 I=1,IRANK
0821 N7=K+N+I
0822 SJ=R(N7)*S(I)
0823 DO 560 J=IRP1,N
0824 SJ=SJ+A(I,J)*S(J)
0825 560 CONTINUE
0826 N6=K+I
0827 DO 570 J=IRP1,N
0828 570 S(J)=S(J)-A(I,J)*R(N6)*SJ
0829 580 S(I)=S(I)-R(N6)*R(N7)*SJ
0830
0831
0832
0833 590 DO 600 I=1,N
0834 600 X(I)=X(I)+S(I)
0835 IF (IN.EQ.7) GO TO 750
0836
0837
0838
0839 TOP1=0.E0
0840 DO 610 J=1,N
0841 N2=J7+J
0842 610 TOP1=AMAX1(TOP1,ABS(S(J))*S(N2))
0843 DO 630 I=1,K
0844 SJ=0.D0
0845 DO 620 J=1,N
0846 N1=J1+J
0847 IP=S(N1)
0848 N7=J2+IP
0849 620 SJ=SJ+AA(I,IP)*X(J)*S(N7)
0850 630 R(I)=B(I)-SJ
0851 IF (ITMAX.LE.0) GO TO 750
0852
0853
0854
0855 TOP=0.E0
0856 DO 640 J=1,N
0857 N2=J7+J
0858 640 TOP=AMAX1(TOP,ABS(X(J))*S(N2))
0859
0860
0861
0862 IF (TOP1-TOP*EPS2) 690,650,650
0863 650 IF (IT-ITMAX) 660,680,680
0864 660 IT=IT+1
0865 IF (IT.EQ.1) GO TO 670
0866 IF (TOP1.GT..25*TOP2) GO TO 690
0867 670 TOP2=TOP1
0868 GO TO (390,460), ISW
0869 680 IT=0
0870 690 SJ=0.D0
0871 DO 700 J=1,K
0872 SJ=SJ+R(J)**2
0873 700 CONTINUE
0874 ENORM=DSQRT(SJ)
0875 IF (IRANK.EQ.N.AND.ISW.EQ.1) GO TO 710
0876 GO TO 730
0877 710 ENM1=ENORM
0878
0879
0880
0881 DO 720 J=1,N
0882 N1=K+J
0883 720 R(N1)=X(J)
0884 ISW=2
0885 IT=0
0886 GO TO 460
0887
0888
0889
0890 730 IF (IRANK.LT.N) GO TO 750
0891 IF (ENORM.LE.ENM1) GO TO 750
0892 DO 740 J=1,N
0893 N1=K+J
0894 740 X(J)=R(N1)
0895 ENORM=ENM1
0896
0897
0898
0899
0900
0901
0902 750 DO 760 J=1,N
0903 N1=J1+J
0904 760 S(J)=S(N1)
0905 DO 790 J=1,N
0906 DO 770 I=J,N
0907 IP=S(I)
0908 IF (J.EQ.IP) GO TO 780
0909 770 CONTINUE
0910 780 S(I)=S(J)
0911 S(J)=J
0912 SJ=X(J)
0913 X(J)=X(I)
0914 790 X(I)=SJ
0915
0916
0917
0918 DO 800 J=1,N
0919 N2=J2+J
0920 800 X(J)=X(J)*S(N2)
0921 RETURN
0922
0923
0924
0925 810 DO 820 J=1,N
0926 N2=J2+J
0927 DO 820 I=1,K
0928 820 A(I,J)=AA(I,J)
0929 RETURN
0930
0931
0932
0933 830 IF (IRANK.EQ.N) RETURN
0934 NS=N-IRANK
0935 DO 840 I=1,N
0936 DO 840 J=1,NS
0937 840 H(I,J)=0.E0
0938 DO 850 J=1,NS
0939 N2=IRANK+J
0940 850 H(N2,J)=1.E0
0941 IF (IRANK.EQ.0) RETURN
0942 DO 870 J=1,IRANK
0943 DO 870 I=1,NS
0944 N7=K+N+J
0945 SJ=R(N7)*H(J,I)
0946 DO 860 K1=IRP1,N
0947 860 SJ=SJ+H(K1,I)*A(J,K1)
0948 N6=K+J
0949 BP=R(N6)
0950 DP=BP*R(N7)*SJ
0951 A1=DP
0952 A2=DP-A1
0953 H(J,I)=H(J,I)-(A1+2.*A2)
0954 DO 870 K1=IRP1,N
0955 DP=BP*A(J,K1)*SJ
0956 A1=DP
0957 A2=DP-A1
0958 870 H(K1,I)=H(K1,I)-(A1+2.*A2)
0959
0960
0961
0962 DO 880 J=1,N
0963 N1=J1+J
0964 880 S(J)=S(N1)
0965 DO 910 J=1,N
0966 DO 890 I=J,N
0967 IP=S(I)
0968 IF (J.EQ.IP) GO TO 900
0969 890 CONTINUE
0970 900 S(I)=S(J)
0971 S(J)=J
0972 DO 910 K1=1,NS
0973 A1=H(J,K1)
0974 H(J,K1)=H(I,K1)
0975 910 H(I,K1)=A1
0976 RETURN
0977
0978 1140 FORMAT (31H0WARNING. IRANK HAS BEEN SET TO,I4,6H BUT(,1PE10.3,9H)
0979 1 RANK IS,I4,25H. IRANK IS NOW TAKEN AS ,I4,1H.)
0980 END
0981 FUNCTION POTEM(T,S,P)
0982 implicit none
0983
0984
0985
0986
0987
0988 real POTEM,T,S,P
0989 real B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11
0990 real T2,T3,S2,P2
0991 B1=-1.60E-5*P
0992 B2=1.014E-5*P*T
0993 T2=T*T
0994 T3=T2*T
0995 B3=-1.27E-7*P*T2
0996 B4=2.7E-9*P*T3
0997 B5=1.322E-6*P*S
0998 B6=-2.62E-8*P*S*T
0999 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
1035
1036 double precision function dunesco(t,s,d)
1037
1038 implicit none
1039 real t,s,d
1040 double precision t2,t3,t4,s2,s3o2,p,p2,r,comp
1041
1042
1043 p = d*0.1D0
1044
1045 t2 = t*t
1046 t3 = t2*t
1047 t4 = t3*t
1048
1049 s2 = s*s
1050 if (s .lt. 0.) then
1051 write(*,*) '****************************'
1052 write(*,*) '* subroutine density_field *'
1053 write(*,*) '****************************'
1054 write(*,*) ' WARNING: salinity takes on'
1055 write(*,*)
1056 & ' negative values '
1057 write(*,*) ' S = ', s
1058 write(*,*) ' This is not supposed to ',
1059 & 'happen. Have a look!'
1060 stop
1061 end if
1062 s3o2 = sqrt(s2*s)
1063
1064
1065
1066
1067 r = 999.842594 + 6.793952e-02 * t
1068 & - 9.095290e-3 * t2 + 1.001685e-04 * t3
1069 & - 1.120083e-06 * t4 + 6.536332e-9 * t4*t
1070 & + s * ( 8.24493e-01 - 4.0899e-3 * t
1071 & + 7.6438e-05 * t2 - 8.2467e-07 * t3
1072 & + 5.3875e-09 * t4 )
1073 & + s3o2 * ( - 5.72466e-03
1074 & + 1.0227e-04 * t - 1.6546e-06 * t2 )
1075 & + 4.8314e-04 * s2
1076
1077
1078
1079
1080 p2 = p*p
1081
1082 comp = 19652.21
1083 & + 148.4206 * t - 2.327105 * t2
1084 & + 1.360477e-02 * t3 - 5.155288e-05 * t4
1085 & + p * ( 3.239908 + 1.43713e-03 * t
1086 & + 1.16092e-04 * t2 - 5.77905e-07 * t3 )
1087 & + p2 * ( 8.50935e-05 - 6.12293e-06 * t
1088 & + 5.2787e-08 * t2 )
1089 & + s * (54.6746 - 0.603459 * t
1090 & + 1.09987e-02 * t2 - 6.1670e-5 * t3 )
1091 & + s3o2 * ( 7.944e-02 + 1.6483e-02 * t
1092 & - 5.3009e-04*t2)
1093 & + p*s * ( 2.2838e-03
1094 & - 1.0981e-05 * t - 1.6078e-06 * t2 )
1095 & + 1.91075e-04 * p*s3o2
1096 & + p2*s * ( - 9.9348e-07
1097 & + 2.0816e-08 * t + 9.1697e-10 * t2 )
1098
1099 dunesco = r/(1.D0-(p/comp))-1000.D0
1100
1101 return
1102 end