utils/knudsen2/unesco.f

a37a13034c Mart*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
                0013       PARAMETER (NumLevels=15)
                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       DOUBLE PRECISION DUNESCO
                0027       real POTEM
                0028 C
                0029 C     Arrays for least squares routine
                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 C     ENTER BOUNDS FOR POLYNOMIAL FIT:
                0049 C           TMin(K) = LOWER BND OF T AT Level K  (Insitu Temperature)
                0050 C           TMax(K) = UPPER BND OF T          "  (Insitu Temperature)
                0051 C           SMin(K) = LOWER BND OF S          "
                0052 C           SMax(K) = UPPER BND OF S          "
                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 ! Local
                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 C  ENTER LEVEL THICKNESSES IN CENTIMETERS
                0074 C     DATA dz/ 4*50.E2, 2*100.E2, 200.E2, 400.E2, 7*500.E2, 6*10.E2 /
                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/ ! ,6.90e+04/
                0079
                0080 C     CALC DEPTHS OF LEVELS FROM DZ (IN METERS)
                0081 C     THE MAXIMUM ALLOWABLE DEPTH IS 8000 METERS
                0082       Z(1) = 0.             ! for level at box edges
                0083       Z(1)=.5*DZ(1)/100.    ! for level at box center
                0084       DO K=2,NumLevels
                0085         Z(K)=Z(K-1)+.5*(DZ(K)+DZ(K-1))/100.
                0086       ENDDO
                0087
                0088
                0089 C     Break the levels up into 250m bands
                0090       DO  I=1,NumLevels
                0091         IDX( I ) = I
                0092 C       Comment out the next line to use input bands as polynomial levels
                0093         IDX( I ) = IFIX(Z(I)/250.)+1
                0094       ENDDO
                0095
                0096
                0097 C     Write some diagnostics
                0098       PRINT 419
                0099       PRINT 422,(Z(I),
                0100      &    Tmin( IDX(I)),Tmax( IDX(I)),
                0101      &    SMin(IDX(I)),Smax( IDX(I)),
                0102 CCC     &    TMin(I),TMax(I),SMin(I),SMax(I),
                0103      &    (   IDX(I) ),I=1,NumLevels )
                0104
                0105
                0106 C     Loop over each level and calculate the polynomial coefficients
                0107
                0108       DO MP=1,NumLevels
                0109
                0110 C     Choose callocation points
                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 C  For each callocation point, convert insitu temperature to
                0123 C  potential temperature and calculate the corresponding density.
                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 C          DenPick(K)=DN(T,S,D)
                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 C     Let (Tbar,Sbar) = the average of (T,S) used in the set of calloc pts
                0142 C     NOTE:  Tbar is still an insitu temperature
                0143 C     Also, calculate the average density of the set of callocation points
                0144
                0145         Tbar=Tsum / FLOAT( NumCallocPt )
                0146         Sbar=Ssum / FLOAT( NumCallocPt )
                0147         DensityBar=DensitySum / FLOAT( NumCallocPt )
                0148
                0149 C     Calculate the average potential temperature of the callocation points
                0150         ThetaBar=ThetaSum/ FLOAT( NumCallocPt )
                0151
                0152 C     Set the reference temperature, salinity and density
                0153 C        DensityRef=DN(Tbar,Sbar,D)
                0154         DensityRef=DUNESCO(Tbar,Sbar,D)
                0155
                0156         SalRef = Sbar
                0157
                0158         TempRef=TBar
                0159         TempRef=ThetaBar   ! DELETE THIS LINE IF USING IN SITU TEMPERATURES
                0160
                0161 C$$$        TempRef=POTEM(TBar,Sbar,D)
                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)   ! DELETE THIS LINE IF USING IN SITU TEMPERATURES
                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 C     SET THE ARGUMENTS IN CALL TO LSQSL2
                0185 C     FIRST DIMENSION OF ARRAY A
                0186         NDIM=NumCallocPt
                0187 C
                0188 C     NUMBER OF ROWS OF A
                0189         M=NumCallocPt
                0190 C
                0191 C     NUMBER OF COLUMNS OF A
                0192         N=NTerms
                0193 C     OPTION NUMBER OF LSQSL2
                0194         IN=1
                0195 C
                0196 C     ITMAX=NUMBER OF ITERATIONS
                0197         ITMAX=100
                0198 C
                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 CCC   PRINT 411
                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 ccc C************************************************************************
                0219 ccc C     WRITE TO UNIT 50
                0220 ccc C************************************************************************
                0221 ccc       OPEN(UNIT=50,FILE='KNUDSEN_COEFS.mit.h',STATUS='UNKNOWN')
                0222 ccc
                0223 ccc C***      WRITE(50,613)
                0224 ccc 613   FORMAT('      REAL SIGREF(KM)')
                0225 ccc       ISEQ=114399990
                0226 ccc       DO  J=1,NumLevels
                0227 ccc         ISEQ=ISEQ+10
                0228 ccc C$$$        WRITE(50,615)J,.01*DZ(J)
                0229 ccc       ENDDO
                0230 ccc   615 FORMAT('      DZ(',I3,')=',F7.2,'E2')
                0231 ccc C$$$      WRITE(50,601)
                0232 ccc 601   FORMAT('C  REFERENCE DENSITIES AT T-POINTS'                )
                0233 ccc       ISEQ=114500030
                0234 ccc       DO J=1,NumLevels
                0235 ccc         ISEQ=ISEQ+10
                0236 ccc C$$$        WRITE(50,501)J,AB(2,J)
                0237 ccc       ENDDO
                0238 ccc 501   FORMAT('      SIGREF(',I3,')=',F8.4)
                0239 ccc
                0240 ccc
                0241 ccc       WRITE(50,'("      DATA SIGREF/")')
                0242 ccc
                0243 ccc       N0 = 1
                0244 ccc c      DO ILINE = 1,NumLevels/5 -1
                0245 ccc 222   CONTINUE
                0246 ccc         N1 = N0 + 4
                0247 ccc         N1 = MIN(N1, NumLevels)
                0248 ccc         WRITE(50,1280) (AB(2,I),I=N0,N1)
                0249 ccc         N0 = N1 + 1
                0250 ccc       IF ( N1 .LT. NumLevels ) GOTO 222
                0251 ccc C     N1 = N0 + 4
                0252 ccc C     IF( N1 .GT. NumLevels ) N1 = NumLevels
                0253 ccc       WRITE(50,1286)
                0254 ccc
                0255
                0256 C**********************************************************************
                0257
                0258 C MITgcmUV
                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 C**********************************************************************
                0269
                0270 c     PRINT 412,((AB(I,J),I=1,NN),J=1,NumLevels)
                0271
                0272 C     WRITE DATA STATEMENTS TO UNIT 50
                0273 caja  DO  L=1,NumLevels
                0274 caja    AB(2,L)=1.E-3*AB(2,L)
                0275 caja    AB(4,L)=1.E-3*AB(4,L)-.035
                0276 caja    AB(5,L)=1.E-3*AB(5,L)
                0277 caja    AB(7,L)=1.E-3*AB(7,L)
                0278 caja    AB(10,L)=1.E-3*AB(10,L)
                0279 caja    AB( 9,L)=1.E+3*AB( 9,L)
                0280 caja    AB(11,L)=1.E+3*AB(11,L)
                0281 caja    AB(13,L)=1.E+6*AB(13,L)
                0282 caja  ENDDO
                0283
                0284 C**********************************************************************
                0285
                0286 C MITgcm (compare01)
                0287
                0288       open(99,file='polyeos.coeffs',form='formatted',status='unknown')
                0289       do J=1,NumLevels
                0290        write(99,*) ab(3,J)  ! ref temperature
                0291       enddo
                0292       do J=1,NumLevels
                0293        write(99,*) ab(4,J)  ! ref sal
                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)  ! ref sig0
                0302       enddo
                0303 CEK      write(99,200)(ab(3,J),J=11,15)  ! ref temperature
                0304 CEK      write(99,200)(ab(4,J),J= 1, 5)  ! ref sal
                0305 CEK      write(99,200)(ab(4,J),J= 6,10)  ! ref sal
                0306 CEK      write(99,200)(ab(4,J),J= 11,15) ! ref sal
                0307 CEK      do I=5,NN
                0308 CEK         write(99,200)(AB(I,J),J= 1, 5)
                0309 CEK         write(99,200)(AB(I,J),J= 6,10)
                0310 CEK         write(99,200)(AB(I,J),J=11,15)
                0311 CEK      enddo
                0312 CEK      write(99,200)(AB(2,J),J= 1, 5)  ! ref sig0
                0313 CEK      write(99,200)(AB(2,J),J= 6,10)  ! ref sig0
                0314 CEK      write(99,200)(AB(2,J),J=11,15)  ! ref sig0
                0315 CEK  200 format(5(E14.7,1X))
                0316
                0317 C**********************************************************************
                0318
                0319
                0320
                0321 ccc       NSEQ=603800000
                0322 ccc       WRITE(50,1298)
                0323 ccc       N=0
                0324 ccc  1260 IS=N+1
                0325 ccc       IE=N+5
                0326 ccc       NSEQ=NSEQ+10
                0327 ccc       IF(IE.LT.NumLevels) THEN
                0328 ccc         WRITE(50,1280) (AB(3,I),I=IS,IE)
                0329 ccc         N=IE
                0330 ccc         GO TO 1260
                0331 ccc       ELSE
                0332 ccc         IE=NumLevels
                0333 ccc         N=IE-IS+1
                0334 ccc         GO TO (1261,1262,1263,1264,1265),N
                0335 ccc  1261     WRITE(50,1281) (AB(3,I),I=IS,IE)
                0336 ccc           GO TO 1268
                0337 ccc  1262     WRITE(50,1282) (AB(3,I),I=IS,IE)
                0338 ccc           GO TO 1268
                0339 ccc  1263     WRITE(50,1283) (AB(3,I),I=IS,IE)
                0340 ccc           GO TO 1268
                0341 ccc  1264     WRITE(50,1284) (AB(3,I),I=IS,IE)
                0342 ccc           GO TO 1268
                0343 ccc  1265     WRITE(50,1285) (AB(3,I),I=IS,IE)
                0344 ccc       ENDIF
                0345 ccc  1268 CONTINUE
                0346 ccc       NSEQ=603900000
                0347 ccc       WRITE(50,1297)
                0348 ccc       N=0
                0349 ccc  1270 IS=N+1
                0350 ccc       IE=N+5
                0351 ccc       NSEQ=NSEQ+10
                0352 ccc       IF(IE.LT.NumLevels) THEN
                0353 ccc         WRITE(50,1280) (AB(4,I),I=IS,IE)
                0354 ccc         N=IE
                0355 ccc         GO TO 1270
                0356 ccc       ELSE
                0357 ccc         IE=NumLevels
                0358 ccc         N=IE-IS+1
                0359 ccc         GO TO (1271,1272,1273,1274,1275),N
                0360 ccc  1271     WRITE(50,1281) (AB(4,I),I=IS,IE)
                0361 ccc           GO TO 1278
                0362 ccc  1272     WRITE(50,1282) (AB(4,I),I=IS,IE)
                0363 ccc           GO TO 1278
                0364 ccc  1273     WRITE(50,1283) (AB(4,I),I=IS,IE)
                0365 ccc           GO TO 1278
                0366 ccc  1274     WRITE(50,1284) (AB(4,I),I=IS,IE)
                0367 ccc           GO TO 1278
                0368 ccc  1275     WRITE(50,1285) (AB(4,I),I=IS,IE)
                0369 ccc       ENDIF
                0370 ccc  1278 CONTINUE
                0371 ccc       DO 1200 L=1,NumLevels
                0372 ccc       IF(L.EQ.1) NSEQ=604000000
                0373 ccc       WRITE(50,1296) L
                0374 ccc       NSEQ=NSEQ+10
                0375 ccc       WRITE(50,1295) (AB(I,L),I=5,8)
                0376 ccc       NSEQ=NSEQ+10
                0377 ccc       WRITE(50,1295) (AB(I,L),I=9,12)
                0378 ccc       NSEQ=NSEQ+10
                0379 ccc       WRITE(50,1294) AB(13,L)
                0380 ccc       NSEQ=NSEQ+10
                0381 ccc  1200 CONTINUE
                0382 ccc  1288 CONTINUE
                0383  1298 FORMAT(18H      DATA TRef  /,67X,I9)
                0384  1297 FORMAT(18H      DATA SRef  /,67X,I9)
                0385 C 419  FORMAT(5X,'LEVEL   TMIN      TMAX      SMIN      SMAX   ',
                0386 C     &          ' TMIN()  Tmax()   Smin()    Smax()  D/250')
                0387 C  422 FORMAT (5X,F6.1,4F10.3,4F10.3,I3)
                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 C****************************************************************************
                0409 *DECK LSQSL2
                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 C     THIS ROUTINE IS A MODIFICATION OF LSQSOL. MARCH,1968. R. HANSON.
                0414 C     LINEAR LEAST SQUARES SOLUTION
                0415 C
                0416 C     THIS ROUTINE FINDS X SUCH THAT THE EUCLIDEAN LENGTH OF
                0417 C     (*) AX-B IS A MINIMUM.
                0418 C
                0419 C     HERE A HAS K ROWS AND N COLUMNS, WHILE B IS A COLUMN VECTOR WITH
                0420 C     K COMPONENTS.
                0421 C
                0422 C     AN ORTHOGONAL MATRIX Q IS FOUND SO THAT QA IS ZERO BELOW
                0423 C     THE MAIN DIAGONAL.
                0424 C     SUPPOSE THAT RANK (A)=R
                0425 C     AN ORTHOGONAL MATRIX S IS FOUND SUCH THAT
                0426 C     QAS=T IS AN R X N UPPER TRIANGULAR MATRIX WHOSE LAST N-R COLUMNS
                0427 C     ARE ZERO.
                0428 C     THE SYSTEM TZ=C (C THE FIRST R COMPONENTS OF QB) IS THEN
                0429 C     SOLVED. WITH W=SZ, THE SOLUTION MAY BE EXPRESSED
                0430 C     AS X = W + SY, WHERE W IS THE SOLUTION OF (*) OF MINIMUM EUCLID-
                0431 C     EAN LENGTH AND Y IS ANY SOLUTION TO (QAS)Y=TY=0.
                0432 C
                0433 C     ITERATIVE IMPROVEMENTS ARE CALCULATED USING RESIDUALS AND
                0434 C     THE ABOVE PROCEDURES WITH B REPLACED BY B-AX, WHERE X IS AN
                0435 C     APPROXIMATE SOLUTION.
                0436 C
                0437       integer ndim,nhdim
                0438       DOUBLE PRECISION SJ,DP,DP1,UP,BP,AJ
                0439       LOGICAL ERM
                0440       INTEGER D,W
                0441 C
                0442 C     IN=1 FOR FIRST ENTRY.
                0443 C                   A IS DECOMPOSED AND SAVED. AX-B IS SOLVED.
                0444 C     IN = 2 FOR SUBSEQUENT ENTRIES WITH A NEW VECTOR B.
                0445 C     IN=3 TO RESTORE A FROM THE PREVIOUS ENTRY.
                0446 C     IN=4 TO CONTINUE THE ITERATIVE IMPROVEMENT FOR THIS SYSTEM.
                0447 C     IN = 5 TO CALCULATE SOLUTIONS TO AX=0, THEN STORE IN THE ARRAY H.
                0448 C     IN  =  6   DO NOT STORE A  IN AA.  OBTAIN  T = QAS, WHERE T IS
                0449 C     MIN(K,N) X MIN(K,N) AND UPPER TRIANGULAR. NOW RETURN.DO NOT OBTAIN
                0450 C     A SOLUTION.
                0451 C     NO SCALING OR COLUMN INTERCHANGES ARE PERFORMED.
                0452 C     IN  =  7   SAME AS WITH  IN = 6  EXCEPT THAT SOLN. OF MIN. LENGTH
                0453 C                IS PLACED INTO X. NO ITERATIVE REFINEMENT.  NOW RETURN.
                0454 C     COLUMN INTERCHANGES ARE PERFORMED. NO SCALING IS PERFORMED.
                0455 C     IN  = 8    SET ADDRESSES. NOW RETURN.
                0456 C
                0457 C     OPTIONS FOR COMPUTING  A MATRIX PRODUCT   Y*H  OR  H*Y ARE
                0458 C     AVAILABLE WITH THE USE OF THE ENTRY POINTS  MYH AND MHY.
                0459 C     USE OF THESE OPTIONS IN THESE ENTRY POINTS ALLOW A GREAT SAVING IN
                0460 C     STORAGE REQUIRED.
                0461 C
                0462 C
                0463       real A(NDIM,NDIM),B(1),AA(D,W),S(1), X(1),H(NHDIM,NHDIM),R(1)
                0464 C     D = DEPTH OF MATRIX.
                0465 C     W = WIDTH OF MATRIX.
                0466 C----
                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 C----
                0473       K=D
                0474       N=W
                0475       ERM=.TRUE.
                0476 C
                0477 C     IF IT=0 ON ENTRY, THE POSSIBLE ERROR MESSAGE WILL BE SUPPRESSED.
                0478 C
                0479       IF (IT.EQ.0) ERM=.FALSE.
                0480 C
                0481 C     IEQ = 2      IF COLUMN SCALING BY LEAST MAX. COLUMN LENGTH IS
                0482 C     TO BE PERFORMED.
                0483 C
                0484 C     IEQ = 1       IF SCALING OF ALL COMPONENTS IS TO BE DONE WITH
                0485 C     THE SCALAR MAX(ABS(AIJ))/K*N.
                0486 C
                0487 C     IEQ = 3 IF COLUMN SCALING AS WITH IN =2 WILL BE RETAINED IN
                0488 C     RANK DEFICIENT CASES.
                0489 C
                0490 C     THE ARRAY S MUST CONTAIN AT LEAST MAX(K,N) + 4N + 4MIN(K,N) CELLS
                0491 C        THE   ARRAY R MUST CONTAIN K+4N S.P. CELLS.
                0492 C
                0493       DATA EPS2/1.E-16/
                0494 C     THE LAST CARD CONTROLS DESIRED RELATIVE ACCURACY.
                0495 C     EPS1  CONTROLS  (EPS) RANK.
                0496 C
                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 C
                0514 C     RETURN AFTER SETTING ADDRESSES WHEN IN=8.
                0515 C
                0516       GO TO (10,360,810,390,830,10,10), IN
                0517 C
                0518 C     EQUILIBRATE COLUMNS OF A (1)-(2).
                0519 C
                0520 C     (1)
                0521 C
                0522    10 CONTINUE
                0523 C
                0524 C     SAVE DATA WHEN IN = 1.
                0525 C
                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 C
                0537 C      S(M+N+1)-S(M+2N) CONTAINS SCALING FOR OUTPUT VARIABLES.
                0538 C
                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 C     COMPUTE COLUMN LENGTHS WITH D.P. SUMS FINALLY ROUNDED TO S.P.
                0559 C
                0560 C     (2)
                0561 C
                0562   100 DO 110 J=1,N
                0563       N7=J7+J
                0564       N2=J2+J
                0565   110 S(N7)=S(N2)
                0566 C
                0567 C      S(M+1)-S(M+ N) CONTAINS VARIABLE PERMUTATIONS.
                0568 C
                0569 C     SET PERMUTATION TO IDENTITY.
                0570 C
                0571       DO 120 J=1,N
                0572       N1=J1+J
                0573   120 S(N1)=J
                0574 C
                0575 C     BEGIN ELIMINATION ON THE MATRIX A WITH ORTHOGONAL MATRICES .
                0576 C
                0577 C     IP=PIVOT ROW
                0578 C
                0579       DO 250 IP=1,LM
                0580 C
                0581 C
                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 C
                0595 C     MAXIMIZE (SIGMA)**2 BY COLUMN INTERCHANGE.
                0596 C
                0597 C      SUPRESS COLUMN INTERCHANGES WHEN IN=6.
                0598 C
                0599 C
                0600 C     EXCHANGE COLUMNS IF NECESSARY.
                0601 C
                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 C
                0608 C     RECORD PERMUTATION AND EXCHANGE SQUARES OF COLUMN LENGTHS.
                0609 C
                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 C
                0629 C     TEST FOR RANK DEFICIENCY.
                0630 C
                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 C
                0639 C     (EPS1) RANK IS DEFICIENT.
                0640 C
                0641   200 SP=DSQRT(DP)
                0642 C
                0643 C     BEGIN FRONT ELIMINATION ON COLUMN IP.
                0644 C
                0645 C     SP=SQROOT(SIGMA**2).
                0646 C
                0647       BP=1.D0/(DP+SP*ABS(A(IP,IP)))
                0648 C
                0649 C     STORE BETA IN S(3N+1)-S(3N+L).
                0650 C
                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 C
                0665 C     SJ=YJ NOW
                0666 C
                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 C
                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 C
                0681 C     BEGIN BACK PROCESSING FOR RANK DEFICIENCY CASE
                0682 C      IF IRANK IS LESS THAN N.
                0683 C
                0684       DO 280 J=1,N
                0685       N2=J2+J
                0686       N7=J7+J
                0687       L=MIN0(J,IRANK)
                0688 C
                0689 C     UNSCALE COLUMNS FOR RANK DEFICIENT MATRICES WHEN IEQ.NE.3.
                0690 C
                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 C
                0704 C     IP TH ELEMENT OF U VECTOR CALCULATED.
                0705 C
                0706       BP=1.D0/(SJ+ABS(A(IP,IP))*AJ)
                0707 C
                0708 C     BP = 2/LENGTH OF U SQUARED.
                0709 C
                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 C
                0719 C     CALC. (AJ,U), WHERE AJ=JTH ROW OF A
                0720 C
                0721       A(I,IP)=A(I,IP)-UP*DP
                0722 C
                0723 C     MODIFY ARRAY A.
                0724 C
                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 C
                0729 C     CALC. MODIFIED PIVOT.
                0730 C
                0731 C
                0732 C     SAVE BETA AND IP TH ELEMENT OF U VECTOR IN R ARRAY.
                0733 C
                0734       N6=K+IP
                0735       N7=K+N+IP
                0736       R(N6)=BP
                0737       R(N7)=UP
                0738 C
                0739 C     TEST FOR END OF BACK PROCESSING.
                0740 C
                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 C
                0749 C     SET INITIAL X VECTOR TO ZERO.
                0750 C
                0751       DO 380 J=1,N
                0752   380 X(J)=0.D0
                0753       IF (IRANK.EQ.0) GO TO 690
                0754 C
                0755 C     APPLY Q TO RT. HAND SIDE.
                0756 C
                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 C
                0785 C     APPLY AT TO RT. HAND SIDE.
                0786 C     APPLY SCALING.
                0787 C
                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 C
                0802 C     ENTRY TO CONTINUE ITERATING.  SOLVES TZ = C = 1ST IRANK
                0803 C     COMPONENTS OF QB .
                0804 C
                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 C
                0815 C     Z CALCULATED.  COMPUTE X = SZ.
                0816 C
                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 C
                0831 C     INCREMENT FOR X OF MINIMAL LENGTH CALCULATED.
                0832 C
                0833   590 DO 600 I=1,N
                0834   600 X(I)=X(I)+S(I)
                0835       IF (IN.EQ.7) GO TO 750
                0836 C
                0837 C     CALC. SUP NORM OF INCREMENT AND RESIDUALS
                0838 C
                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 C
                0853 C     CALC. SUP NORM OF X.
                0854 C
                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 C
                0860 C     COMPARE RELATIVE CHANGE IN X WITH TOLERANCE EPS .
                0861 C
                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 C
                0879 C     SAVE X ARRAY.
                0880 C
                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 C
                0888 C     CHOOSE BEST SOLUTION
                0889 C
                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 C
                0897 C     NORM OF AX - B LOCATED IN THE CELL ENORM .
                0898 C
                0899 C
                0900 C     REARRANGE VARIABLES.
                0901 C
                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 C
                0916 C     SCALE VARIABLES.
                0917 C
                0918       DO 800 J=1,N
                0919       N2=J2+J
                0920   800 X(J)=X(J)*S(N2)
                0921       RETURN
                0922 C
                0923 C     RESTORE A.
                0924 C
                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 C
                0931 C     GENERATE SOLUTIONS TO THE HOMOGENEOUS EQUATION AX = 0.
                0932 C
                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 C
                0960 C     REARRANGE ROWS OF SOLUTION MATRIX.
                0961 C
                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 C
                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 C     POTENTIAL TEMPERATURE FUNCTION
                0984 C     BASED ON FOFONOFF AND FROESE (1958) AS SHOWN IN "THE SEA" VOL. 1,
                0985 C     PAGE 17, TABLE IV
                0986 C     INPUT IS TEMPERATURE, SALINITY, PRESSURE (OR DEPTH)
                0987 C     UNITS ARE DEG.C., PPT, DBARS (OR METERS)
                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 c     convert from meters to bars
                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 ***   density at one standard atmosphere pressure (p=0) in si-units
                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 ***   density at depth d[m] = p[bar] :
                1078 ***   include depth effect of secant bulk modulus k(s,t,p) if desired
                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