utils/knudsen2/knudsen2.f

725bce7a44 Alis*0001       PROGRAM KNUDSEN2
                0002       implicit none
                0003 C
                0004 C     COEFFICIENTS FOR DENSITY COMPUTATION
                0005 C
                0006 C     THIS PROGRAM CALCULATES THE COEFFICIENTS OF THE POLYNOMIAL
                0007 C     APPROXIMATION TO THE KNUDSEN FORMULA.  THE PROGRAM IS SET UP
                0008 C     TO YIELD COEFFICIENTS THAT WILL COMPUTE DENSITY AS A FUNCTION
                0009 C     OF POTENTIAL TEMPERATURE, SALINITY, AND DEPTH.
                0010 C
                0011 C     Number of levels
                0012       integer NumLevels
5e81fe432e Alis*0013       PARAMETER (NumLevels=15)
725bce7a44 Alis*0014
                0015 C     Number of Temperature and Salinity callocation points
                0016       integer NumTempPt,NumSalPt,NumCallocPt
                0017       PARAMETER (NumTempPt = 10, NumSalPt = 5 )
                0018       PARAMETER (NumCallocPt =NumSalPt*NumTempPt )
                0019
                0020 C     Number of terms in fit polynomial
                0021       integer NTerms
                0022       PARAMETER (NTerms = 9 )
                0023
                0024 C     Functions (Density and Potential temperature)
                0025       DOUBLE PRECISION DN
                0026       real POTEM
                0027 C
                0028 C     Arrays for least squares routine
                0029       real A(NumCallocPt,NumCallocPt),
                0030      &    B(NumCallocPt),
                0031      &    X(NTerms),
                0032      &    AB(13,NumLevels),
                0033      &    C(NumCallocPt,NTerms),
                0034      &    H(NTerms,NTerms),
                0035      &    R(NumCallocPt*2),SB(NumCallocPt*4)
                0036
                0037       real
                0038      &    Z(NumLevels),
                0039      &    DZ(NumLevels)
                0040       integer
                0041      &    IDX(NumLevels)
                0042
                0043       real
                0044      &    Theta(NumCallocPt),
                0045      &    SPick(NumCallocPt),TPick(NumCallocPt),DenPick(NumCallocPt)
                0046
                0047 C     ENTER BOUNDS FOR POLYNOMIAL FIT:
                0048 C           TMin(K) = LOWER BND OF T AT Level K  (Insitu Temperature)
                0049 C           TMax(K) = UPPER BND OF T          "  (Insitu Temperature)
                0050 C           SMin(K) = LOWER BND OF S          "
                0051 C           SMax(K) = UPPER BND OF S          "
                0052
                0053       real TMin(25), TMax(25)
                0054      &    ,SMin(25), SMax(25)
                0055
                0056
                0057       DATA TMin / 4*-2., 15*-1., 6*0.     /
                0058       DATA TMax / 29., 19., 14., 11., 9., 7., 5., 3*4., 5*3., 10*2. /
                0059       DATA SMin / 28.5, 33.7, 34., 34.1, 34.2, 34.4, 2*34.5,
                0060      &    15*34.6, 2*34.7 /
                0061       DATA SMax / 36.7, 36.6, 35.8, 35.7, 35.3, 2*35.1, 7*35.,
                0062      &    9*34.9, 2*34.8 /
                0063
                0064 ! Local
                0065       integer I,J,K,IT,IEQ,IRANK,NDIM,NHDIM,N,M,IN,ITMAX,MP,L,NN
                0066       real T,S,D,EPS,DeltaT,DeltaS,ENORM,DeltaDen,DensityRef
                0067       real Tbar,Sbar,Tsum,SSum,TempInc,SalnInc
                0068       real DensitySum,ThetaSum
                0069       real DensityBar,thetabar,TempRef,SAlRef
                0070
                0071
                0072 C  ENTER LEVEL THICKNESSES IN CENTIMETERS
                0073 C     DATA dz/ 4*50.E2, 2*100.E2, 200.E2, 400.E2, 7*500.E2, 6*10.E2 /
5e81fe432e Alis*0074        DATA dz/   5.00e+03,7.00e+03,1.00e+04,1.40e+04,
                0075      &1.90e+04,2.40e+04,2.90e+04,3.40e+04,3.90e+04,
                0076      &4.40e+04,4.90e+04,5.40e+04,5.90e+04,6.40e+04,
                0077      &6.90e+04/ ! ,6.90e+04/
725bce7a44 Alis*0078
                0079 C     CALC DEPTHS OF LEVELS FROM DZ (IN METERS)
                0080 C     THE MAXIMUM ALLOWABLE DEPTH IS 8000 METERS
                0081       Z(1) = 0.             ! for level at box edges
                0082       Z(1)=.5*DZ(1)/100.    ! for level at box center
                0083       DO K=2,NumLevels
                0084         Z(K)=Z(K-1)+.5*(DZ(K)+DZ(K-1))/100.
                0085       ENDDO
                0086
                0087
                0088 C     Break the levels up into 250m bands
                0089       DO  I=1,NumLevels
                0090         IDX( I ) = I
                0091 C       Comment out the next line to use input bands as polynomial levels
                0092         IDX( I ) = IFIX(Z(I)/250.)+1
                0093       ENDDO
                0094
                0095
                0096 C     Write some diagnostics
                0097       PRINT 419
                0098       PRINT 422,(Z(I),
                0099      &    Tmin( IDX(I)),Tmax( IDX(I)),
                0100      &    SMin(IDX(I)),Smax( IDX(I)),
                0101 CCC     &    TMin(I),TMax(I),SMin(I),SMax(I),
                0102      &    (   IDX(I) ),I=1,NumLevels )
                0103
                0104
                0105 C     Loop over each level and calculate the polynomial coefficients
                0106
                0107       DO MP=1,NumLevels
                0108
                0109 C     Choose callocation points
                0110         TempInc =(Tmax( IDX(MP))-TMin(IDX(MP)))/(2.*FLOAT(NumSalPt)-1.0)
                0111         SalnInc =(Smax( IDX(MP))-SMin(IDX(MP)))/(FLOAT(NumSalPt)-1.0)
                0112         DO I=1,NumTempPt
                0113           DO J=1,NumSalPt
                0114             K=NumSalPt*I+J-NumSalPt
                0115             TPick(K)=TMin(IDX(MP))+(FLOAT(I)-1.0)*TempInc
                0116             SPick(K)=SMin(IDX(MP))+(FLOAT(J)-1.0)*SalnInc
                0117           ENDDO
                0118         ENDDO
                0119
                0120
                0121 C  For each callocation point, convert insitu temperature to
                0122 C  potential temperature and calculate the corresponding density.
                0123         Tsum=0.0
                0124         Ssum=0.0
                0125         DensitySum=0.0
                0126         ThetaSum=0.0
                0127         DO K=1,NumCallocPt
                0128           D=Z(MP)
                0129           S=SPick(K)
                0130           T=TPick(K)
                0131           DenPick(K)=DN(T,S,D)
                0132           Theta(K)=POTEM(T,S,D)
                0133           Tsum=Tsum+TPick(K)
                0134           Ssum=Ssum+SPick(K)
                0135           DensitySum = DensitySum+DenPick(K)
                0136           ThetaSum   = ThetaSum+Theta(K)
                0137         ENDDO
                0138
                0139 C     Let (Tbar,Sbar) = the average of (T,S) used in the set of calloc pts
                0140 C     NOTE:  Tbar is still an insitu temperature
                0141 C     Also, calculate the average density of the set of callocation points
                0142
                0143         Tbar=Tsum / FLOAT( NumCallocPt )
                0144         Sbar=Ssum / FLOAT( NumCallocPt )
                0145         DensityBar=DensitySum / FLOAT( NumCallocPt )
                0146
                0147 C     Calculate the average potential temperature of the callocation points
                0148         ThetaBar=ThetaSum/ FLOAT( NumCallocPt )
                0149
                0150 C     Set the reference temperature, salinity and density
                0151         DensityRef=DN(Tbar,Sbar,D)
                0152
                0153         SalRef = Sbar
                0154
                0155         TempRef=TBar
                0156         TempRef=ThetaBar   ! DELETE THIS LINE IF USING IN SITU TEMPERATURES
                0157
                0158 C$$$        TempRef=POTEM(TBar,Sbar,D)
                0159
                0160
                0161         AB(1,MP)=Z(MP)
                0162         AB(2,MP)=DensityRef
                0163         AB(3,MP)=TempRef
                0164         AB(4,MP)=SalRef
                0165         DO K=1,NumCallocPt
                0166           TPick(K)=Theta(K)   ! DELETE THIS LINE IF USING IN SITU TEMPERATURES
                0167           DeltaT   = TPick(K)   - TempRef
                0168           DeltaS   = SPick(K)   - SalRef
                0169           DeltaDen = DenPick(K) - DensityRef
                0170           B(K)= DeltaDen
                0171           A(K,1)=DeltaT
                0172           A(K,2)=DeltaS
                0173           A(K,3)=DeltaT*DeltaT
                0174           A(K,4)=DeltaT*DeltaS
                0175           A(K,5)=DeltaS*DeltaS
                0176           A(K,6)=A(K,3)*DeltaT
                0177           A(K,7)=A(K,4)*DeltaT
                0178           A(K,8)=A(K,4)*DeltaS
                0179           A(K,9)=A(K,5)*DeltaS
                0180         ENDDO
                0181 C     SET THE ARGUMENTS IN CALL TO LSQSL2
                0182 C     FIRST DIMENSION OF ARRAY A
                0183         NDIM=NumCallocPt
                0184 C
                0185 C     NUMBER OF ROWS OF A
                0186         M=NumCallocPt
                0187 C
                0188 C     NUMBER OF COLUMNS OF A
                0189         N=NTerms
                0190 C     OPTION NUMBER OF LSQSL2
                0191         IN=1
                0192 C
                0193 C     ITMAX=NUMBER OF ITERATIONS
                0194         ITMAX=100
                0195 C
                0196         IT=0
                0197         IEQ=2
                0198         IRANK=0
                0199         EPS=1.0E-7
                0200         EPS=1.0E-11
                0201         NHDIM=NTerms
                0202         CALL LSQSL2(NDIM,A,M,N,B,X,IRANK,IN,ITMAX,IT,IEQ,ENORM,EPS,
                0203      &      NHDIM,H,C,R,SB)
                0204
                0205 CCC   PRINT 411
                0206         DO I=1,N
                0207           AB(I+4,MP)=X(I)
                0208         ENDDO
                0209       ENDDO
                0210       PRINT 430
                0211   430 FORMAT(1X,'  Z    SIG0    T    S       X1          X2
                0212      1    X3          X4          X5          X6          X7          X8
                0213      2    X9     ',/)
                0214       NN=N+4
                0215 ccc C************************************************************************
                0216 ccc C     WRITE TO UNIT 50
                0217 ccc C************************************************************************
                0218 ccc       OPEN(UNIT=50,FILE='KNUDSEN_COEFS.mit.h',STATUS='UNKNOWN')
                0219 ccc
                0220 ccc C***      WRITE(50,613)
                0221 ccc 613   FORMAT('      REAL SIGREF(KM)')
                0222 ccc       ISEQ=114399990
                0223 ccc       DO  J=1,NumLevels
                0224 ccc         ISEQ=ISEQ+10
                0225 ccc C$$$        WRITE(50,615)J,.01*DZ(J)
                0226 ccc       ENDDO
                0227 ccc   615 FORMAT('      DZ(',I3,')=',F7.2,'E2')
                0228 ccc C$$$      WRITE(50,601)
                0229 ccc 601   FORMAT('C  REFERENCE DENSITIES AT T-POINTS'                )
                0230 ccc       ISEQ=114500030
                0231 ccc       DO J=1,NumLevels
                0232 ccc         ISEQ=ISEQ+10
                0233 ccc C$$$        WRITE(50,501)J,AB(2,J)
                0234 ccc       ENDDO
                0235 ccc 501   FORMAT('      SIGREF(',I3,')=',F8.4)
                0236 ccc
                0237 ccc
                0238 ccc       WRITE(50,'("      DATA SIGREF/")')
                0239 ccc
                0240 ccc       N0 = 1
                0241 ccc c      DO ILINE = 1,NumLevels/5 -1
                0242 ccc 222   CONTINUE
                0243 ccc         N1 = N0 + 4
                0244 ccc         N1 = MIN(N1, NumLevels)
                0245 ccc         WRITE(50,1280) (AB(2,I),I=N0,N1)
                0246 ccc         N0 = N1 + 1
                0247 ccc       IF ( N1 .LT. NumLevels ) GOTO 222
                0248 ccc C     N1 = N0 + 4
                0249 ccc C     IF( N1 .GT. NumLevels ) N1 = NumLevels
                0250 ccc       WRITE(50,1286)
                0251 ccc
                0252
                0253 C**********************************************************************
                0254
                0255 C MITgcmUV
                0256
                0257       open(99,file='POLY3.COEFFS',form='formatted',status='unknown')
                0258       write(99,*) NumLevels
a37a13034c Mart*0259       write(99,'(3(G20.13))') (ab(3,J),ab(4,J),ab(2,J),J=1 , NumLevels)
725bce7a44 Alis*0260       do J=1,NumLevels
a37a13034c Mart*0261          write(99,'(3(G20.13))')(AB(I,J),I=5,13)
725bce7a44 Alis*0262       enddo
                0263       close(99)
                0264
                0265 C**********************************************************************
                0266
                0267 c     PRINT 412,((AB(I,J),I=1,NN),J=1,NumLevels)
                0268
                0269 C     WRITE DATA STATEMENTS TO UNIT 50
                0270 caja  DO  L=1,NumLevels
                0271 caja    AB(2,L)=1.E-3*AB(2,L)
                0272 caja    AB(4,L)=1.E-3*AB(4,L)-.035
                0273 caja    AB(5,L)=1.E-3*AB(5,L)
                0274 caja    AB(7,L)=1.E-3*AB(7,L)
                0275 caja    AB(10,L)=1.E-3*AB(10,L)
                0276 caja    AB( 9,L)=1.E+3*AB( 9,L)
                0277 caja    AB(11,L)=1.E+3*AB(11,L)
                0278 caja    AB(13,L)=1.E+6*AB(13,L)
                0279 caja  ENDDO
                0280
                0281 C**********************************************************************
                0282
                0283 C MITgcm (compare01)
                0284
                0285       open(99,file='polyeos.coeffs',form='formatted',status='unknown')
                0286       do J=1,NumLevels
                0287        write(99,*) ab(3,J)  ! ref temperature
                0288       enddo
                0289       do J=1,NumLevels
                0290        write(99,*) ab(4,J)  ! ref sal
                0291       enddo
                0292       do J=1,NumLevels
                0293        do I=5,13
                0294         write(99,*) AB(I,J)
                0295        enddo
                0296       enddo
                0297       do J=1,NumLevels
                0298        write(99,*) ab(2,J)  ! ref sig0
                0299       enddo
                0300 CEK      write(99,200)(ab(3,J),J=11,15)  ! ref temperature
                0301 CEK      write(99,200)(ab(4,J),J= 1, 5)  ! ref sal
                0302 CEK      write(99,200)(ab(4,J),J= 6,10)  ! ref sal
                0303 CEK      write(99,200)(ab(4,J),J= 11,15) ! ref sal
                0304 CEK      do I=5,NN
                0305 CEK         write(99,200)(AB(I,J),J= 1, 5)
                0306 CEK         write(99,200)(AB(I,J),J= 6,10)
                0307 CEK         write(99,200)(AB(I,J),J=11,15)
                0308 CEK      enddo
                0309 CEK      write(99,200)(AB(2,J),J= 1, 5)  ! ref sig0
                0310 CEK      write(99,200)(AB(2,J),J= 6,10)  ! ref sig0
                0311 CEK      write(99,200)(AB(2,J),J=11,15)  ! ref sig0
                0312 CEK  200 format(5(E14.7,1X))
                0313
                0314 C**********************************************************************
                0315
                0316
                0317
                0318 ccc       NSEQ=603800000
                0319 ccc       WRITE(50,1298)
                0320 ccc       N=0
                0321 ccc  1260 IS=N+1
                0322 ccc       IE=N+5
                0323 ccc       NSEQ=NSEQ+10
                0324 ccc       IF(IE.LT.NumLevels) THEN
                0325 ccc         WRITE(50,1280) (AB(3,I),I=IS,IE)
                0326 ccc         N=IE
                0327 ccc         GO TO 1260
                0328 ccc       ELSE
                0329 ccc         IE=NumLevels
                0330 ccc         N=IE-IS+1
                0331 ccc         GO TO (1261,1262,1263,1264,1265),N
                0332 ccc  1261     WRITE(50,1281) (AB(3,I),I=IS,IE)
                0333 ccc           GO TO 1268
                0334 ccc  1262     WRITE(50,1282) (AB(3,I),I=IS,IE)
                0335 ccc           GO TO 1268
                0336 ccc  1263     WRITE(50,1283) (AB(3,I),I=IS,IE)
                0337 ccc           GO TO 1268
                0338 ccc  1264     WRITE(50,1284) (AB(3,I),I=IS,IE)
                0339 ccc           GO TO 1268
                0340 ccc  1265     WRITE(50,1285) (AB(3,I),I=IS,IE)
                0341 ccc       ENDIF
                0342 ccc  1268 CONTINUE
                0343 ccc       NSEQ=603900000
                0344 ccc       WRITE(50,1297)
                0345 ccc       N=0
                0346 ccc  1270 IS=N+1
                0347 ccc       IE=N+5
                0348 ccc       NSEQ=NSEQ+10
                0349 ccc       IF(IE.LT.NumLevels) THEN
                0350 ccc         WRITE(50,1280) (AB(4,I),I=IS,IE)
                0351 ccc         N=IE
                0352 ccc         GO TO 1270
                0353 ccc       ELSE
                0354 ccc         IE=NumLevels
                0355 ccc         N=IE-IS+1
                0356 ccc         GO TO (1271,1272,1273,1274,1275),N
                0357 ccc  1271     WRITE(50,1281) (AB(4,I),I=IS,IE)
                0358 ccc           GO TO 1278
                0359 ccc  1272     WRITE(50,1282) (AB(4,I),I=IS,IE)
                0360 ccc           GO TO 1278
                0361 ccc  1273     WRITE(50,1283) (AB(4,I),I=IS,IE)
                0362 ccc           GO TO 1278
                0363 ccc  1274     WRITE(50,1284) (AB(4,I),I=IS,IE)
                0364 ccc           GO TO 1278
                0365 ccc  1275     WRITE(50,1285) (AB(4,I),I=IS,IE)
                0366 ccc       ENDIF
                0367 ccc  1278 CONTINUE
                0368 ccc       DO 1200 L=1,NumLevels
                0369 ccc       IF(L.EQ.1) NSEQ=604000000
                0370 ccc       WRITE(50,1296) L
                0371 ccc       NSEQ=NSEQ+10
                0372 ccc       WRITE(50,1295) (AB(I,L),I=5,8)
                0373 ccc       NSEQ=NSEQ+10
                0374 ccc       WRITE(50,1295) (AB(I,L),I=9,12)
                0375 ccc       NSEQ=NSEQ+10
                0376 ccc       WRITE(50,1294) AB(13,L)
                0377 ccc       NSEQ=NSEQ+10
                0378 ccc  1200 CONTINUE
                0379 ccc  1288 CONTINUE
                0380  1298 FORMAT(18H      DATA TRef  /,67X,I9)
                0381  1297 FORMAT(18H      DATA SRef  /,67X,I9)
                0382 C 419  FORMAT(5X,'LEVEL   TMIN      TMAX      SMIN      SMAX   ',
                0383 C     &          ' TMIN()  Tmax()   Smin()    Smax()  D/250')
                0384 C  422 FORMAT (5X,F6.1,4F10.3,4F10.3,I3)
                0385   419 FORMAT(5X,'LEVEL   TMIN      TMAX      SMIN      SMAX  D/250')
                0386   422 FORMAT (5X,F6.1,4F10.3,I3)
                0387   412 FORMAT(1X,F6.1,F8.4,F7.3,F6.2,9E12.5)
                0388  1296 FORMAT(6X,'DATA (C(',I2,',N),N=1,9)/',54X,I9)
                0389  1295 FORMAT(5X,1H*,9X,4(E13.7,1H,),10X,I9)
                0390  1294 FORMAT(5X,1H*,9X,E13.7,1H/,52X,I9)
                0391  1280 FORMAT(5X,1H*,8X,5(F10.7,1H,),12X,I9)
                0392  1281 FORMAT(5X,1H*,8X,F10.7,1H/,56X,I9)
                0393  1282 FORMAT(5X,1H*,8X,F10.7,1H,,F10.7,1H/,45X,I9)
                0394  1283 FORMAT(5X,1H*,8X,2(F10.7,1H,),F10.7,1H/,34X,I9)
                0395  1284 FORMAT(5X,1H*,8X,3(F10.7,1H,),F10.7,1H/,23X,I9)
                0396  1285 FORMAT(5X,1H*,8X,4(F10.7,1H,),F10.7,1H/,12X,I9)
                0397  1286 FORMAT(5X,1H*,1H/)
                0398   350 FORMAT(1X,E14.7)
                0399   351 FORMAT(1H+,T86,E14.7/)
                0400   400 FORMAT(1X,9E14.7)
                0401   410 FORMAT(1X,5E14.7)
                0402   411 FORMAT(///)
                0403       STOP
                0404       END
                0405 C****************************************************************************
                0406 *DECK LSQSL2
                0407       SUBROUTINE LSQSL2 (NDIM,A,D,W,B,X,IRANK,IN,ITMAX,IT,IEQ,ENORM,EPS1
                0408      1,NHDIM,H,AA,R,S)
                0409       implicit none
                0410 C     THIS ROUTINE IS A MODIFICATION OF LSQSOL. MARCH,1968. R. HANSON.
                0411 C     LINEAR LEAST SQUARES SOLUTION
                0412 C
                0413 C     THIS ROUTINE FINDS X SUCH THAT THE EUCLIDEAN LENGTH OF
                0414 C     (*) AX-B IS A MINIMUM.
                0415 C
                0416 C     HERE A HAS K ROWS AND N COLUMNS, WHILE B IS A COLUMN VECTOR WITH
                0417 C     K COMPONENTS.
                0418 C
                0419 C     AN ORTHOGONAL MATRIX Q IS FOUND SO THAT QA IS ZERO BELOW
                0420 C     THE MAIN DIAGONAL.
                0421 C     SUPPOSE THAT RANK (A)=R
                0422 C     AN ORTHOGONAL MATRIX S IS FOUND SUCH THAT
                0423 C     QAS=T IS AN R X N UPPER TRIANGULAR MATRIX WHOSE LAST N-R COLUMNS
                0424 C     ARE ZERO.
                0425 C     THE SYSTEM TZ=C (C THE FIRST R COMPONENTS OF QB) IS THEN
                0426 C     SOLVED. WITH W=SZ, THE SOLUTION MAY BE EXPRESSED
                0427 C     AS X = W + SY, WHERE W IS THE SOLUTION OF (*) OF MINIMUM EUCLID-
                0428 C     EAN LENGTH AND Y IS ANY SOLUTION TO (QAS)Y=TY=0.
                0429 C
                0430 C     ITERATIVE IMPROVEMENTS ARE CALCULATED USING RESIDUALS AND
                0431 C     THE ABOVE PROCEDURES WITH B REPLACED BY B-AX, WHERE X IS AN
                0432 C     APPROXIMATE SOLUTION.
                0433 C
                0434       integer ndim,nhdim
                0435       DOUBLE PRECISION SJ,DP,DP1,UP,BP,AJ
                0436       LOGICAL ERM
                0437       INTEGER D,W
                0438 C
                0439 C     IN=1 FOR FIRST ENTRY.
                0440 C                   A IS DECOMPOSED AND SAVED. AX-B IS SOLVED.
                0441 C     IN = 2 FOR SUBSEQUENT ENTRIES WITH A NEW VECTOR B.
                0442 C     IN=3 TO RESTORE A FROM THE PREVIOUS ENTRY.
                0443 C     IN=4 TO CONTINUE THE ITERATIVE IMPROVEMENT FOR THIS SYSTEM.
                0444 C     IN = 5 TO CALCULATE SOLUTIONS TO AX=0, THEN STORE IN THE ARRAY H.
                0445 C     IN  =  6   DO NOT STORE A  IN AA.  OBTAIN  T = QAS, WHERE T IS
                0446 C     MIN(K,N) X MIN(K,N) AND UPPER TRIANGULAR. NOW RETURN.DO NOT OBTAIN
                0447 C     A SOLUTION.
                0448 C     NO SCALING OR COLUMN INTERCHANGES ARE PERFORMED.
                0449 C     IN  =  7   SAME AS WITH  IN = 6  EXCEPT THAT SOLN. OF MIN. LENGTH
                0450 C                IS PLACED INTO X. NO ITERATIVE REFINEMENT.  NOW RETURN.
                0451 C     COLUMN INTERCHANGES ARE PERFORMED. NO SCALING IS PERFORMED.
                0452 C     IN  = 8    SET ADDRESSES. NOW RETURN.
                0453 C
                0454 C     OPTIONS FOR COMPUTING  A MATRIX PRODUCT   Y*H  OR  H*Y ARE
                0455 C     AVAILABLE WITH THE USE OF THE ENTRY POINTS  MYH AND MHY.
                0456 C     USE OF THESE OPTIONS IN THESE ENTRY POINTS ALLOW A GREAT SAVING IN
                0457 C     STORAGE REQUIRED.
                0458 C
                0459 C
                0460       real A(NDIM,NDIM),B(1),AA(D,W),S(1), X(1),H(NHDIM,NHDIM),R(1)
                0461 C     D = DEPTH OF MATRIX.
                0462 C     W = WIDTH OF MATRIX.
                0463 C----
                0464       integer K,N,IT,ISW,L,M,IRANK,IEQ,IN,K1
                0465       integer J1,J2,J3,J4,J5,J6,J7,J8,J9
                0466       integer N1,N2,N3,N4,N5,N6,N7,N8,NS
                0467       integer I,ITMAX,IPM1,II,LM,J,IP,KM,IPP1,IRP1,IRM1
                0468       real SP,ENORM,TOP,TOP1,ENM1,TOP2,EPS1,EPS2,A1,A2,AM
                0469 C----
                0470       K=D
                0471       N=W
                0472       ERM=.TRUE.
                0473 C
                0474 C     IF IT=0 ON ENTRY, THE POSSIBLE ERROR MESSAGE WILL BE SUPPRESSED.
                0475 C
                0476       IF (IT.EQ.0) ERM=.FALSE.
                0477 C
                0478 C     IEQ = 2      IF COLUMN SCALING BY LEAST MAX. COLUMN LENGTH IS
                0479 C     TO BE PERFORMED.
                0480 C
                0481 C     IEQ = 1       IF SCALING OF ALL COMPONENTS IS TO BE DONE WITH
                0482 C     THE SCALAR MAX(ABS(AIJ))/K*N.
                0483 C
                0484 C     IEQ = 3 IF COLUMN SCALING AS WITH IN =2 WILL BE RETAINED IN
                0485 C     RANK DEFICIENT CASES.
                0486 C
                0487 C     THE ARRAY S MUST CONTAIN AT LEAST MAX(K,N) + 4N + 4MIN(K,N) CELLS
                0488 C        THE   ARRAY R MUST CONTAIN K+4N S.P. CELLS.
                0489 C
                0490       DATA EPS2/1.E-16/
                0491 C     THE LAST CARD CONTROLS DESIRED RELATIVE ACCURACY.
                0492 C     EPS1  CONTROLS  (EPS) RANK.
                0493 C
                0494       ISW=1
                0495       L=MIN0(K,N)
                0496       M=MAX0(K,N)
                0497       J1=M
                0498       J2=N+J1
                0499       J3=J2+N
                0500       J4=J3+L
                0501       J5=J4+L
                0502       J6=J5+L
                0503       J7=J6+L
                0504       J8=J7+N
                0505       J9=J8+N
                0506       LM=L
                0507       IF (IRANK.GE.1.AND.IRANK.LE.L) LM=IRANK
                0508       IF (IN.EQ.6) LM=L
                0509       IF (IN.EQ.8) RETURN
                0510 C
                0511 C     RETURN AFTER SETTING ADDRESSES WHEN IN=8.
                0512 C
                0513       GO TO (10,360,810,390,830,10,10), IN
                0514 C
                0515 C     EQUILIBRATE COLUMNS OF A (1)-(2).
                0516 C
                0517 C     (1)
                0518 C
                0519    10 CONTINUE
                0520 C
                0521 C     SAVE DATA WHEN IN = 1.
                0522 C
                0523       IF (IN.GT.5) GO TO 30
                0524       DO 20 J=1,N
                0525       DO 20 I=1,K
                0526    20 AA(I,J)=A(I,J)
                0527    30 CONTINUE
                0528       IF (IEQ.EQ.1) GO TO 60
                0529       DO 50 J=1,N
                0530       AM=0.E0
                0531       DO 40 I=1,K
                0532    40 AM=AMAX1(AM,ABS(A(I,J)))
                0533 C
                0534 C      S(M+N+1)-S(M+2N) CONTAINS SCALING FOR OUTPUT VARIABLES.
                0535 C
                0536       N2=J2+J
                0537       IF (IN.EQ.6) AM=1.E0
                0538       S(N2)=1.E0/AM
                0539       DO 50 I=1,K
                0540    50 A(I,J)=A(I,J)*S(N2)
                0541       GO TO 100
                0542    60 AM=0.E0
                0543       DO 70 J=1,N
                0544       DO 70 I=1,K
                0545    70 AM=AMAX1(AM,ABS(A(I,J)))
                0546       AM=AM/FLOAT(K*N)
                0547       IF (IN.EQ.6) AM=1.E0
                0548       DO 80 J=1,N
                0549       N2=J2+J
                0550    80 S(N2)=1.E0/AM
                0551       DO 90 J=1,N
                0552       N2=J2+J
                0553       DO 90 I=1,K
                0554    90 A(I,J)=A(I,J)*S(N2)
                0555 C     COMPUTE COLUMN LENGTHS WITH D.P. SUMS FINALLY ROUNDED TO S.P.
                0556 C
                0557 C     (2)
                0558 C
                0559   100 DO 110 J=1,N
                0560       N7=J7+J
                0561       N2=J2+J
                0562   110 S(N7)=S(N2)
                0563 C
                0564 C      S(M+1)-S(M+ N) CONTAINS VARIABLE PERMUTATIONS.
                0565 C
                0566 C     SET PERMUTATION TO IDENTITY.
                0567 C
                0568       DO 120 J=1,N
                0569       N1=J1+J
                0570   120 S(N1)=J
                0571 C
                0572 C     BEGIN ELIMINATION ON THE MATRIX A WITH ORTHOGONAL MATRICES .
                0573 C
                0574 C     IP=PIVOT ROW
                0575 C
                0576       DO 250 IP=1,LM
                0577 C
                0578 C
                0579       DP=0.D0
                0580       KM=IP
                0581       DO 140 J=IP,N
                0582       SJ=0.D0
                0583       DO 130 I=IP,K
                0584       SJ=SJ+A(I,J)**2
                0585   130 CONTINUE
                0586       IF (DP.GT.SJ) GO TO 140
                0587       DP=SJ
                0588       KM=J
                0589       IF (IN.EQ.6) GO TO 160
                0590   140 CONTINUE
                0591 C
                0592 C     MAXIMIZE (SIGMA)**2 BY COLUMN INTERCHANGE.
                0593 C
                0594 C      SUPRESS COLUMN INTERCHANGES WHEN IN=6.
                0595 C
                0596 C
                0597 C     EXCHANGE COLUMNS IF NECESSARY.
                0598 C
                0599       IF (KM.EQ.IP) GO TO 160
                0600       DO 150 I=1,K
                0601       A1=A(I,IP)
                0602       A(I,IP)=A(I,KM)
                0603   150 A(I,KM)=A1
                0604 C
                0605 C     RECORD PERMUTATION AND EXCHANGE SQUARES OF COLUMN LENGTHS.
                0606 C
                0607       N1=J1+KM
                0608       A1=S(N1)
                0609       N2=J1+IP
                0610       S(N1)=S(N2)
                0611       S(N2)=A1
                0612       N7=J7+KM
                0613       N8=J7+IP
                0614       A1=S(N7)
                0615       S(N7)=S(N8)
                0616       S(N8)=A1
                0617   160 IF (IP.EQ.1) GO TO 180
                0618       A1=0.E0
                0619       IPM1=IP-1
                0620       DO 170 I=1,IPM1
                0621       A1=A1+A(I,IP)**2
                0622   170 CONTINUE
                0623       IF (A1.GT.0.E0) GO TO 190
                0624   180 IF (DP.GT.0.D0) GO TO 200
                0625 C
                0626 C     TEST FOR RANK DEFICIENCY.
                0627 C
                0628   190 IF (DSQRT(DP/A1).GT.EPS1) GO TO 200
                0629       IF (IN.EQ.6) GO TO 200
                0630       II=IP-1
                0631       IF (ERM) WRITE (6,1140) IRANK,EPS1,II,II
                0632       IRANK=IP-1
                0633       ERM=.FALSE.
                0634       GO TO 260
                0635 C
                0636 C     (EPS1) RANK IS DEFICIENT.
                0637 C
                0638   200 SP=DSQRT(DP)
                0639 C
                0640 C     BEGIN FRONT ELIMINATION ON COLUMN IP.
                0641 C
                0642 C     SP=SQROOT(SIGMA**2).
                0643 C
                0644       BP=1.D0/(DP+SP*ABS(A(IP,IP)))
                0645 C
                0646 C     STORE BETA IN S(3N+1)-S(3N+L).
                0647 C
                0648       IF (IP.EQ.K) BP=0.D0
                0649       N3=K+2*N+IP
                0650       R(N3)=BP
                0651       UP=DSIGN(DBLE(SP)+ABS(A(IP,IP)),DBLE(A(IP,IP)))
                0652       IF (IP.GE.K) GO TO 250
                0653       IPP1=IP+1
                0654       IF (IP.GE.N) GO TO 240
                0655       DO 230 J=IPP1,N
                0656       SJ=0.D0
                0657       DO 210 I=IPP1,K
                0658   210 SJ=SJ+A(I,J)*A(I,IP)
                0659       SJ=SJ+UP*A(IP,J)
                0660       SJ=BP*SJ
                0661 C
                0662 C     SJ=YJ NOW
                0663 C
                0664       DO 220 I=IPP1,K
                0665   220 A(I,J)=A(I,J)-A(I,IP)*SJ
                0666   230 A(IP,J)=A(IP,J)-SJ*UP
                0667   240 A(IP,IP)=-SIGN(SP,A(IP,IP))
                0668 C
                0669       N4=K+3*N+IP
                0670       R(N4)=UP
                0671   250 CONTINUE
                0672       IRANK=LM
                0673   260 IRP1=IRANK+1
                0674       IRM1=IRANK-1
                0675       IF (IRANK.EQ.0.OR.IRANK.EQ.N) GO TO 360
                0676       IF (IEQ.EQ.3) GO TO 290
                0677 C
                0678 C     BEGIN BACK PROCESSING FOR RANK DEFICIENCY CASE
                0679 C      IF IRANK IS LESS THAN N.
                0680 C
                0681       DO 280 J=1,N
                0682       N2=J2+J
                0683       N7=J7+J
                0684       L=MIN0(J,IRANK)
                0685 C
                0686 C     UNSCALE COLUMNS FOR RANK DEFICIENT MATRICES WHEN IEQ.NE.3.
                0687 C
                0688       DO 270 I=1,L
                0689   270 A(I,J)=A(I,J)/S(N7)
                0690       S(N7)=1.E0
                0691   280 S(N2)=1.E0
                0692   290 IP=IRANK
                0693   300 SJ=0.D0
                0694       DO 310 J=IRP1,N
                0695       SJ=SJ+A(IP,J)**2
                0696   310 CONTINUE
                0697       SJ=SJ+A(IP,IP)**2
                0698       AJ=DSQRT(SJ)
                0699       UP=DSIGN(AJ+ABS(A(IP,IP)),DBLE(A(IP,IP)))
                0700 C
                0701 C     IP TH ELEMENT OF U VECTOR CALCULATED.
                0702 C
                0703       BP=1.D0/(SJ+ABS(A(IP,IP))*AJ)
                0704 C
                0705 C     BP = 2/LENGTH OF U SQUARED.
                0706 C
                0707       IPM1=IP-1
                0708       IF (IPM1.LE.0) GO TO 340
                0709       DO 330 I=1,IPM1
                0710       DP=A(I,IP)*UP
                0711       DO 320 J=IRP1,N
                0712       DP=DP+A(I,J)*A(IP,J)
                0713   320 CONTINUE
                0714       DP=DP/(SJ+ABS(A(IP,IP))*AJ)
                0715 C
                0716 C     CALC. (AJ,U), WHERE AJ=JTH ROW OF A
                0717 C
                0718       A(I,IP)=A(I,IP)-UP*DP
                0719 C
                0720 C     MODIFY ARRAY A.
                0721 C
                0722       DO 330 J=IRP1,N
                0723   330 A(I,J)=A(I,J)-A(IP,J)*DP
                0724   340 A(IP,IP)=-DSIGN(AJ,DBLE(A(IP,IP)))
                0725 C
                0726 C     CALC. MODIFIED PIVOT.
                0727 C
                0728 C
                0729 C     SAVE BETA AND IP TH ELEMENT OF U VECTOR IN R ARRAY.
                0730 C
                0731       N6=K+IP
                0732       N7=K+N+IP
                0733       R(N6)=BP
                0734       R(N7)=UP
                0735 C
                0736 C     TEST FOR END OF BACK PROCESSING.
                0737 C
                0738       IF (IP-1) 360,360,350
                0739   350 IP=IP-1
                0740       GO TO 300
                0741   360 IF (IN.EQ.6) RETURN
                0742       DO 370 J=1,K
                0743   370 R(J)=B(J)
                0744       IT=0
                0745 C
                0746 C     SET INITIAL X VECTOR TO ZERO.
                0747 C
                0748       DO 380 J=1,N
                0749   380 X(J)=0.D0
                0750       IF (IRANK.EQ.0) GO TO 690
                0751 C
                0752 C     APPLY Q TO RT. HAND SIDE.
                0753 C
                0754   390 DO 430 IP=1,IRANK
                0755       N4=K+3*N+IP
                0756       SJ=R(N4)*R(IP)
                0757       IPP1=IP+1
                0758       IF (IPP1.GT.K) GO TO 410
                0759       DO 400 I=IPP1,K
                0760   400 SJ=SJ+A(I,IP)*R(I)
                0761   410 N3=K+2*N+IP
                0762       BP=R(N3)
                0763       IF (IPP1.GT.K) GO TO 430
                0764       DO 420 I=IPP1,K
                0765   420 R(I)=R(I)-BP*A(I,IP)*SJ
                0766   430 R(IP)=R(IP)-BP*R(N4)*SJ
                0767       DO 440 J=1,IRANK
                0768   440 S(J)=R(J)
                0769       ENORM=0.E0
                0770       IF (IRP1.GT.K) GO TO 510
                0771       DO 450 J=IRP1,K
                0772   450 ENORM=ENORM+R(J)**2
                0773       ENORM=SQRT(ENORM)
                0774       GO TO 510
                0775   460 DO 480 J=1,N
                0776       SJ=0.D0
                0777       N1=J1+J
                0778       IP=S(N1)
                0779       DO 470 I=1,K
                0780   470 SJ=SJ+R(I)*AA(I,IP)
                0781 C
                0782 C     APPLY AT TO RT. HAND SIDE.
                0783 C     APPLY SCALING.
                0784 C
                0785       N7=J2+IP
                0786       N1=K+N+J
                0787   480 R(N1)=SJ*S(N7)
                0788       N1=K+N
                0789       S(1)=R(N1+1)/A(1,1)
                0790       IF (N.EQ.1) GO TO 510
                0791       DO 500 J=2,N
                0792       N1=J-1
                0793       SJ=0.D0
                0794       DO 490 I=1,N1
                0795   490 SJ=SJ+A(I,J)*S(I)
                0796       N2=K+J+N
                0797   500 S(J)=(R(N2)-SJ)/A(J,J)
                0798 C
                0799 C     ENTRY TO CONTINUE ITERATING.  SOLVES TZ = C = 1ST IRANK
                0800 C     COMPONENTS OF QB .
                0801 C
                0802   510 S(IRANK)=S(IRANK)/A(IRANK,IRANK)
                0803       IF (IRM1.EQ.0) GO TO 540
                0804       DO 530 J=1,IRM1
                0805       N1=IRANK-J
                0806       N2=N1+1
                0807       SJ=0.
                0808       DO 520 I=N2,IRANK
                0809   520 SJ=SJ+A(N1,I)*S(I)
                0810   530 S(N1)=(S(N1)-SJ)/A(N1,N1)
                0811 C
                0812 C     Z CALCULATED.  COMPUTE X = SZ.
                0813 C
                0814   540 IF (IRANK.EQ.N) GO TO 590
                0815       DO 550 J=IRP1,N
                0816   550 S(J)=0.E0
                0817       DO 580 I=1,IRANK
                0818       N7=K+N+I
                0819       SJ=R(N7)*S(I)
                0820       DO 560 J=IRP1,N
                0821       SJ=SJ+A(I,J)*S(J)
                0822   560 CONTINUE
                0823       N6=K+I
                0824       DO 570 J=IRP1,N
                0825   570 S(J)=S(J)-A(I,J)*R(N6)*SJ
                0826   580 S(I)=S(I)-R(N6)*R(N7)*SJ
                0827 C
                0828 C     INCREMENT FOR X OF MINIMAL LENGTH CALCULATED.
                0829 C
                0830   590 DO 600 I=1,N
                0831   600 X(I)=X(I)+S(I)
                0832       IF (IN.EQ.7) GO TO 750
                0833 C
                0834 C     CALC. SUP NORM OF INCREMENT AND RESIDUALS
                0835 C
                0836       TOP1=0.E0
                0837       DO 610 J=1,N
                0838       N2=J7+J
                0839   610 TOP1=AMAX1(TOP1,ABS(S(J))*S(N2))
                0840       DO 630 I=1,K
                0841       SJ=0.D0
                0842       DO 620 J=1,N
                0843       N1=J1+J
                0844       IP=S(N1)
                0845       N7=J2+IP
                0846   620 SJ=SJ+AA(I,IP)*X(J)*S(N7)
                0847   630 R(I)=B(I)-SJ
                0848       IF (ITMAX.LE.0) GO TO 750
                0849 C
                0850 C     CALC. SUP NORM OF X.
                0851 C
                0852       TOP=0.E0
                0853       DO 640 J=1,N
                0854       N2=J7+J
                0855   640 TOP=AMAX1(TOP,ABS(X(J))*S(N2))
                0856 C
                0857 C     COMPARE RELATIVE CHANGE IN X WITH TOLERANCE EPS .
                0858 C
                0859       IF (TOP1-TOP*EPS2) 690,650,650
                0860   650 IF (IT-ITMAX) 660,680,680
                0861   660 IT=IT+1
                0862       IF (IT.EQ.1) GO TO 670
                0863       IF (TOP1.GT..25*TOP2) GO TO 690
                0864   670 TOP2=TOP1
                0865       GO TO (390,460), ISW
                0866   680 IT=0
                0867   690 SJ=0.D0
                0868       DO 700 J=1,K
                0869       SJ=SJ+R(J)**2
                0870   700 CONTINUE
                0871       ENORM=DSQRT(SJ)
                0872       IF (IRANK.EQ.N.AND.ISW.EQ.1) GO TO 710
                0873       GO TO 730
                0874   710 ENM1=ENORM
                0875 C
                0876 C     SAVE X ARRAY.
                0877 C
                0878       DO 720 J=1,N
                0879       N1=K+J
                0880   720 R(N1)=X(J)
                0881       ISW=2
                0882       IT=0
                0883       GO TO 460
                0884 C
                0885 C     CHOOSE BEST SOLUTION
                0886 C
                0887   730 IF (IRANK.LT.N) GO TO 750
                0888       IF (ENORM.LE.ENM1) GO TO 750
                0889       DO 740 J=1,N
                0890       N1=K+J
                0891   740 X(J)=R(N1)
                0892       ENORM=ENM1
                0893 C
                0894 C     NORM OF AX - B LOCATED IN THE CELL ENORM .
                0895 C
                0896 C
                0897 C     REARRANGE VARIABLES.
                0898 C
                0899   750 DO 760 J=1,N
                0900       N1=J1+J
                0901   760 S(J)=S(N1)
                0902       DO 790 J=1,N
                0903       DO 770 I=J,N
                0904       IP=S(I)
                0905       IF (J.EQ.IP) GO TO 780
                0906   770 CONTINUE
                0907   780 S(I)=S(J)
                0908       S(J)=J
                0909       SJ=X(J)
                0910       X(J)=X(I)
                0911   790 X(I)=SJ
                0912 C
                0913 C     SCALE VARIABLES.
                0914 C
                0915       DO 800 J=1,N
                0916       N2=J2+J
                0917   800 X(J)=X(J)*S(N2)
                0918       RETURN
                0919 C
                0920 C     RESTORE A.
                0921 C
                0922   810 DO 820 J=1,N
                0923       N2=J2+J
                0924       DO 820 I=1,K
                0925   820 A(I,J)=AA(I,J)
                0926       RETURN
                0927 C
                0928 C     GENERATE SOLUTIONS TO THE HOMOGENEOUS EQUATION AX = 0.
                0929 C
                0930   830 IF (IRANK.EQ.N) RETURN
                0931       NS=N-IRANK
                0932       DO 840 I=1,N
                0933       DO 840 J=1,NS
                0934   840 H(I,J)=0.E0
                0935       DO 850 J=1,NS
                0936       N2=IRANK+J
                0937   850 H(N2,J)=1.E0
                0938       IF (IRANK.EQ.0) RETURN
                0939       DO 870 J=1,IRANK
                0940       DO 870 I=1,NS
                0941       N7=K+N+J
                0942       SJ=R(N7)*H(J,I)
                0943       DO 860 K1=IRP1,N
                0944   860 SJ=SJ+H(K1,I)*A(J,K1)
                0945       N6=K+J
                0946       BP=R(N6)
                0947       DP=BP*R(N7)*SJ
                0948       A1=DP
                0949       A2=DP-A1
                0950       H(J,I)=H(J,I)-(A1+2.*A2)
                0951       DO 870 K1=IRP1,N
                0952       DP=BP*A(J,K1)*SJ
                0953       A1=DP
                0954       A2=DP-A1
                0955   870 H(K1,I)=H(K1,I)-(A1+2.*A2)
                0956 C
                0957 C     REARRANGE ROWS OF SOLUTION MATRIX.
                0958 C
                0959       DO 880 J=1,N
                0960       N1=J1+J
                0961   880 S(J)=S(N1)
                0962       DO 910 J=1,N
                0963       DO 890 I=J,N
                0964       IP=S(I)
                0965       IF (J.EQ.IP) GO TO 900
                0966   890 CONTINUE
                0967   900 S(I)=S(J)
                0968       S(J)=J
                0969       DO 910 K1=1,NS
                0970       A1=H(J,K1)
                0971       H(J,K1)=H(I,K1)
                0972   910 H(I,K1)=A1
                0973       RETURN
                0974 C
                0975  1140 FORMAT (31H0WARNING. IRANK HAS BEEN SET TO,I4,6H  BUT(,1PE10.3,9H)
                0976      1 RANK IS,I4,25H.  IRANK IS NOW TAKEN AS ,I4,1H.)
                0977       END
                0978       FUNCTION POTEM(T,S,P)
                0979       implicit none
                0980 C     POTENTIAL TEMPERATURE FUNCTION
                0981 C     BASED ON FOFONOFF AND FROESE (1958) AS SHOWN IN "THE SEA" VOL. 1,
                0982 C     PAGE 17, TABLE IV
                0983 C     INPUT IS TEMPERATURE, SALINITY, PRESSURE (OR DEPTH)
                0984 C     UNITS ARE DEG.C., PPT, DBARS (OR METERS)
                0985       real POTEM,T,S,P
                0986       real B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11
                0987       real T2,T3,S2,P2
                0988       B1=-1.60E-5*P
                0989       B2=1.014E-5*P*T
                0990       T2=T*T
                0991       T3=T2*T
                0992       B3=-1.27E-7*P*T2
                0993       B4=2.7E-9*P*T3
                0994       B5=1.322E-6*P*S
                0995       B6=-2.62E-8*P*S*T
                0996       S2=S*S
                0997       P2=P*P
                0998       B7=4.1E-9*P*S2
                0999       B8=9.14E-9*P2
                1000       B9=-2.77E-10*P2*T
                1001       B10=9.5E-13*P2*T2
                1002       B11=-1.557E-13*P2*P
                1003       POTEM=B1+B2+B3+B4+B5+B6+B7+B8+B9+B10+B11
                1004       POTEM=T-POTEM
                1005       RETURN
                1006       END
                1007       FUNCTION DN(T,S,D)
                1008       implicit none
                1009       real T,S,D
                1010       DOUBLE PRECISION DN,T3,S2,T2,S3,F1,F2,F3,FS,SIGMA,A,B1,B2,B,CO,
                1011      1ALPHA,ALPSTD
                1012       T2 = T*T
                1013       T3= T2* T
                1014       S2 = S*S
                1015       S3 = S2 * S
                1016       F1 = -(T-3.98)**2 * (T+283.)/(503.57*(T+67.26))
                1017       F2 = T3*1.0843E-6 - T2*9.8185E-5 + T*4.786E-3
                1018       F3 = T3*1.667E-8 - T2*8.164E-7 + T*1.803E-5
                1019       FS= S3*6.76786136D-6 - S2*4.8249614D-4 + S*8.14876577D-1
                1020       SIGMA= F1 + (FS+3.895414D-2)*(1.-F2+F3*(FS-.22584586D0))
                1021       A= D*1.0E-4*(105.5+ T*9.50 - T2*0.158 - D*T*1.5E-4)  -
                1022      1(227. + T*28.33 - T2*0.551 + T3* 0.004)
                1023       B1 = (FS-28.1324)/10.0
                1024       B2 = B1 * B1
                1025       B= -B1* (147.3-T*2.72 + T2*0.04 - D*1.0E-4*(32.4- 0.87*T+.02*T2))
                1026       B= B+ B2*(4.5-0.1*T - D*1.0E-4*(1.8-0.06*T))
                1027       CO = 4886./(1. + 1.83E-5*D)
                1028       ALPHA=     D*1.0E-6* (CO+A+B)
                1029       DN=(SIGMA+ALPHA)/(1.-1.E-3*ALPHA)
                1030       RETURN
                1031       END