      PROGRAM TSZINT
      USE MOD_ZA  ! HYCOM array I/O interface
      USE netcdf  ! NetCDF fortran 90 interface
      IMPLICIT NONE
C
C     DEFINE INPUT CLIMATOLOGY GRID.
C
C     SETUP FOR 0.25 DEGREE GDEM 4.2 NEAR-GLOBAL CLIMATOLOGY.
C
      INTEGER    IWI,JWI,KWI
      REAL*4     XFIN,YFIN,DXIN,DYIN
      PARAMETER (IWI=1440, JWI=697, KWI=80)
      PARAMETER (XFIN=0.0, YFIN=-84.0, DXIN=0.25, DYIN=0.25)
C
C     CLIM ARRAYS.
C
      INTEGER, ALLOCATABLE :: MSK(:,:)
      REAL*4,  ALLOCATABLE :: PLON(:,:),PLAT(:,:),XAF(:,:),YAF(:,:)
      REAL*4,  ALLOCATABLE :: TM(:,:),SM(:,:),RM(:,:),RT(:,:)
C
      INTEGER   KSIGMA
      REAL*4    XAMAX,XAMIN,YAMAX,YAMIN
      INTEGER*2 TSEAI2(IWI,JWI),SSEAI2(IWI,JWI)
      REAL*4    ADD_OFF,SCALE_F,DBAR,DEPTH,TINSIT
      REAL*8    ZLEV(KWI)
      CHARACTER PREAMBL(5)*79
C
C     NETCDF I/O VARIABLES.
C
      CHARACTER*(256) CFILED,CFILET,CFILES
      INTEGER         ncFIDd,ncVIDd,ncFIDt,ncVIDt,ncFIDs,ncVIDs
C
C     INTERPOLATION ARRAYS.
C
      INTEGER IBD(4)
      REAL*4  TSEAI(IWI+4,JWI+4),SSEAI(IWI+4,JWI+4)
      REAL*4  FXI(IWI+4),FYI(JWI+4),
     +        WQSEA3(IWI+4,JWI+4,3),WK(3*(IWI+JWI+8)+1)
C
C     NAMELIST.
C
      INTEGER          JPR
      COMMON/NPROCS/   JPR
      SAVE  /NPROCS/
C
      CHARACTER*40     CTITLE
      NAMELIST/AFTITL/ CTITLE
      INTEGER          ICTYPE,SIGVER,INTERP,MONTH,ITEST,JTEST
      NAMELIST/AFFLAG/ ICTYPE,SIGVER,INTERP,MONTH,ITEST,JTEST,JPR
C
C**********
C*
C 1)  FROM UNFORMATTED Sig0-T OR Sig0-S OR T-S DATA ON ITS NATIVE GRID,
C      CREATE A FORMATTED MODEL GRID CLIM FILE SUITABLE FOR INPUT TO THE
C      HYCOM ISOPYCNAL CLIMATOLOGY GENERATOR OVER THE GIVEN REGION.
C
C     INTERPOLATION IS EITHER PIECEWISE BILINEAR OR CUBIC SPLINE.
C
C 2)  PARAMETERS:
C
C     NATIVE CLIM GRID SPECIFICATION, SEE (4):
C
C        IWI    = 1ST DIMENSION OF CLIM GRID
C        JWI    = 2ND DIMENSION OF CLIM GRID
C        JWI    = 3RD DIMENSION OF CLIM GRID (NUMBER OF Z-LEVELS)
C        XFIN   = LONGITUDE OF 1ST CLIM GRID POINT
C        YFIN   = LATITUDE  OF 1ST CLIM GRID POINT
C        DXIN   = CLIM LONGITUDINAL GRID SPACING
C        DYIN   = CLIM LATITUDINAL  GRID SPACING
C
C 3)  NAMELIST INPUT:
C
C     /AFTITL/
C        CTITLE - ONE (40-CHARACTER) LINE TITLE.
C
C     /AFFLAG/
C        ICTYPE -  INPUT FILE TYPE
C                   =5; IN-SITU TEMPERATURE AND SALINITY
C                   =6; IN-SITU TEMPERATURE AND SALINITY, UNSTABLE
C        SIGVER - EQUATION OF STATE TYPE
C                   =1-2; 7-term
C                   =3-4; 9-term
C                   =5-6; 17-term
C                   =7-8; 12-term
C                   =odd;  sigma-0
C                   =even; sigma-2
C        INTERP - INTERPOLATION FLAG.
C                   =0; PIECEWISE LINEAR
C                   =1; CUBIC SPLINE (DEFAULT)
C        MONTH  - MONTH OF CLIMATOLOGY (1 TO 12)
C        ITEST  - 1ST ARRAY INDEX FOR DIAGNOSTIC PRINTOUT
C                   =0; NO DIAGNOSTIC PRINTOUT
C        JTEST  - 2ND ARRAY INDEX FOR DIAGNOSTIC PRINTOUT
C                   =0; NO DIAGNOSTIC PRINTOUT
C
C 4)  INPUT:
C        ON UNIT  5: NAMELIST /AFTITL/, /AFTIME/
C        netCDF NATIVE T CLIM FILE (ENV.VAR. CDF_TEMP), SEE (5).
C        netCDF NATIVE S CLIM FILE (ENV.VAR. CDF_SALN), SEE (5).
C     OUTPUT:
C        ON UNIT 11:    UNFORMATTED MODEL    CLIM FILE, SEE (6).
C        ON UNIT 12:    UNFORMATTED MODEL    CLIM FILE, SEE (6).
C        ON UNIT 13:    UNFORMATTED MODEL    CLIM FILE, SEE (6).
C
C 5)  THE INPUT CLIM FIELDS, VIA NetCDF, ARE ON THE 'NATIVE' LAT-LON
C      GRID, STARTING AT THE POINT 'XFIN' EAST AND 'YFIN' NORTH WITH 
C      'YFIN' NORTH WITH GRIDSIZE 'DXIN' BY 'DYIN' DEGREES.  THE
C      INPUT ARRAY SIZE IS 'IWI' BY 'JWI', AND THERE ARE NO REPEATED
C      NODES (EVEN FOR GLOBAL DATA SETS).  IN-SITU TEMERATURE IS INPUT
C      AND CONVERTED HERE TO POTENTIAL TEMPERATURE FOR OUTPUT.
C
C     ALL CLIMATOLOGY FIELDS MUST BE DEFINED AT EVERY GRID POINT,
C      INCLUDING LAND AND BELOW THE OCEAN FLOOR.
C
C     IF ICTYPE=6, DON'T INFORCE STABLE DENSITY STRATIFICATION AND
C      USE A SALINITY DEPENDENT FREEZING POINT.
C
C 6)  THE OUTPUT CLIMS ARE AT EVERY GRID POINT OF THE MODEL'S 'P' GRID.
C     ARRAY SIZE IS 'IDM' BY 'JDM'.
C
C 7)  SEVERAL STATISTICS ARE WRITTEN OUT IN ORDER TO CHECK THE 
C      INTERPOLATION BETWEEN VARIOUS MACHINES.  MIN, MAX, MEAN AND RMS 
C      OF THE ENTIRE BASIN ARE OUTPUT FOR (TM, SM, RM) FOR EACH 
C      RECORD.  NOTE HOWEVER THAT THESE VALUES MAY NOT REPRESENT THE
C      STATISTICS OF THE CLIMS AS SEEN BY THE MODEL, IF THE INPUT CLIM 
C      DATA HAS NON-REALISTIC VALUES OVER LAND.  IT IS UP TO THE USER 
C      TO CHECK THE LOG FILES FOR CONSISTENCY BETWEEN MACHINES.
C
C 8)  ALAN J. WALLCRAFT,  PLANNING SYSTEMS INC.,  OCTOBER 1995.
C      BASED ON EARILER VERSIONS BY SEVERAL AUTHORS.
C*
C**********
C
      EXTERNAL THETA,P80, LINEAR, AVERMS,MINMAX
      REAL     THETA,P80
C
      REAL*4     ZERO,RADIAN
      PARAMETER (ZERO=0.0, RADIAN=57.2957795)
C
      CHARACTER*80 CLINE
      CHARACTER*40 CCNAME
      INTEGER IUNIT,IWIT,JWIT,KWIT
      REAL*4  XFINT,YFINT,DXINT,DYINT
      REAL*4  HMINA,HMINB,HMAXA,HMAXB
C
      INTEGER I,IWIX,J,KREC
      REAL*4  XFD,YFD,DXD,DYD
      REAL*4  XLIN,XFDX,XOV,YOU,
     +        XMIN,XMAX,XAVE,XRMS
C
      CHARACTER*11 CMONTH(12)
      DATA CMONTH /  ', January',
     +               ', February',
     +               ', March',
     +               ', April',
     +               ', May',
     +               ', June',
     +               ', July',
     +               ', August',
     +               ', September',
     +               ', October',
     +               ', November',
     +               ', December'  /
C
C --- MODEL ARRAYS.
C
      CALL XCSPMD  !define idm,jdm
      ALLOCATE(  MSK(IDM,JDM) )
      ALLOCATE( PLON(IDM,JDM) )
      ALLOCATE( PLAT(IDM,JDM) )
      ALLOCATE(  XAF(IDM,JDM) )
      ALLOCATE(  YAF(IDM,JDM) )
      ALLOCATE(   TM(IDM,JDM) )
      ALLOCATE(   SM(IDM,JDM) )
      ALLOCATE(   RM(IDM,JDM) )
      ALLOCATE(   RT(IDM,JDM) )
C
C     NAMELIST INPUT.
C
      CALL ZHOPEN(6, 'FORMATTED', 'UNKNOWN', 0)
C
      CTITLE = ' '
      WRITE(6,*) 'READING /AFTITL/'
      CALL ZHFLSH(6)
      READ( 5,AFTITL)
      WRITE(6,AFTITL)
C
      ICTYPE = 5
      SIGVER = 0
      INTERP = 1
      MONTH  = 1
      ITEST  = 0
      JTEST  = 0
      JPR    = 8
      WRITE(6,*) 'READING /AFFLAG/'
      CALL ZHFLSH(6)
      READ( 5,AFFLAG)
      WRITE(6,AFFLAG)
      WRITE(6,*)
      CALL ZHFLSH(6)
C
      IF     (SIGVER.LT.1 .OR. SIGVER.GT.8) THEN
        WRITE(6,*)
        WRITE(6,*) 'ERROR - SIGVER MUST BE BETWEEN 0 AND 8'
        WRITE(6,*)
        STOP
      ENDIF
      IF     (SIGVER.GE.40) THEN 
        KSIGMA = 4
      ELSEIF (MOD(SIGVER,2).EQ.1) THEN 
        KSIGMA = 0  !odd  sigver
      ELSE
        KSIGMA = 2  !even sigver
      ENDIF
C
C     CONVERT ICTYPE FROM 5 OR 6 TO 3 OR 4, FOR LEGACY CODE
C
      ICTYPE = ICTYPE - 2
C
C     GRID INPUT.
C
      CALL ZAIOST
C
      CALL ZHOPNC(21, 'regional.grid.b', 'FORMATTED', 'OLD', 0)
      CALL ZAIOPF('regional.grid.a', 'OLD', 21)
C
      READ(21,*) ! skip idm
      READ(21,*) ! skip jdm
      READ(21,*) ! skip mapflg
      READ(21,'(A)') CLINE
      I = INDEX(CLINE,'=')
      READ (CLINE(I+1:),*)   HMINB,HMAXB
      CALL ZAIORD(PLON,MSK,.FALSE., HMINA,HMAXA, 21)
      IF     (ABS(HMINA-HMINB).GT.ABS(HMINB)*1.E-4 .OR.
     &        ABS(HMAXA-HMAXB).GT.ABS(HMAXB)*1.E-4     ) THEN
        WRITE(6,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)')
     &    'error - .a and .b grid files not consistent (plon):',
     &    '.a,.b min = ',HMINA,HMINB,HMINA-HMINB,
     &    '.a,.b max = ',HMAXA,HMAXB,HMAXA-HMAXB
        CALL ZHFLSH(6)
        STOP
      ENDIF
C
      READ(21,'(A)') CLINE
      I = INDEX(CLINE,'=')
      READ (CLINE(I+1:),*)   HMINB,HMAXB
      CALL ZAIORD(PLAT,MSK,.FALSE., HMINA,HMAXA, 21)
      IF     (ABS(HMINA-HMINB).GT.ABS(HMINB)*1.E-4 .OR.
     &        ABS(HMAXA-HMAXB).GT.ABS(HMAXB)*1.E-4     ) THEN
        WRITE(6,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)')
     &    'error - .a and .b grid files not consistent (plat):',
     &    '.a,.b min = ',HMINA,HMINB,HMINA-HMINB,
     &    '.a,.b max = ',HMAXA,HMAXB,HMAXA-HMAXB
        CALL ZHFLSH(6)
        STOP
      ENDIF
C
      CLOSE(UNIT=21)
      CALL ZAIOCL(21)
C
C     INITIALIZE OUTPUT.
C
      CTITLE = trim(CTITLE) // CMONTH(MONTH)
      WRITE(6,6000) 'OUTPUT:',CTITLE
      CALL ZHFLSH(6)
C
      CALL ZAIOPN('NEW', 10)
      CALL ZAIOPN('NEW', 11)
      CALL ZAIOPN('NEW', 12)
C
      CALL ZHOPEN(10, 'FORMATTED', 'NEW', 0)
      CALL ZHOPEN(11, 'FORMATTED', 'NEW', 0)
      CALL ZHOPEN(12, 'FORMATTED', 'NEW', 0)
C
      PREAMBL(1) = CTITLE
      PREAMBL(2) = ' '
      PREAMBL(3) = ' '
      PREAMBL(4) = ' '
      WRITE(PREAMBL(5),'(A,2I5)')
     +        'i/jdm =',
     +       IDM,JDM
C
      PREAMBL(2) = 'Potential Temperature'
      WRITE(10,4101) PREAMBL
C
      PREAMBL(2) = 'Salinity'
      WRITE(11,4101) PREAMBL
C
      IF     (KSIGMA.EQ.0) THEN
        WRITE(PREAMBL(2),"(a,i3)")
     &    'Potential Density (Sigma-0), sigver =',sigver
        WRITE(12,4101) PREAMBL
      ELSEIF (KSIGMA.EQ.2) THEN
        WRITE(PREAMBL(2),"(a,i3)")
     &    'Potential Density (Sigma-2), sigver =',sigver
        WRITE(12,4101) PREAMBL
      ELSEIF (KSIGMA.EQ.4) THEN
        WRITE(PREAMBL(2),"(a,i3)")
     &    'Potential Density (Sigma-4), sigver =',sigver
        WRITE(12,4101) PREAMBL
      ELSE
        WRITE(6,*)
        WRITE(6,*) 'ERROR - KSIGMA MUST BE 0 OR 2 OR 4'
        WRITE(6,*)
        STOP
      ENDIF
C
      WRITE(6,*)
      WRITE(6, 4101) PREAMBL
      WRITE(6,*)
C
C     INITIALIZE CLIMS.
C
      CALL GETENV('CDF_TEMP',CFILET)
      WRITE(6,*)
      WRITE(6,*) 'CDF_TEMP = ',trim(CFILET)
      CALL ZHFLSH(6)
      ! open NetCDF file
      call ncheck(nf90_open(trim(CFILET), nf90_nowrite, ncFIDt))
      ! get ZLEV
      call ncheck(nf90_inq_varid(ncFIDt,'depth',ncVIDt))
      call ncheck(nf90_get_var(  ncFIDt,        ncVIDt,ZLEV(:)))
      ! inquire variable ID
      call ncheck(nf90_inq_varid(ncFIDt,
     &                           'water_temp',
     &                           ncVIDt))
C
      CALL GETENV('CDF_SALN',CFILES)
      WRITE(6,*)
      WRITE(6,*) 'CDF_SALN = ',trim(CFILES)
      CALL ZHFLSH(6)
      ! open NetCDF file
      call ncheck(nf90_open(trim(CFILES), nf90_nowrite, ncFIDs))
      ! inquire variable ID
      call ncheck(nf90_inq_varid(ncFIDs,
     &                           'salinity',
     &                           ncVIDs))
C
C     DEFINE THE GRID COORDINATES.
C
      IF     (IWI*DXIN.GE.359.9) THEN
        IF     (ABS(IWI * DXIN - 360.0) .GT. 0.01) THEN
          WRITE(6,9050)
          CALL ZHFLSH(6)
          STOP
        ENDIF
        IWIX = IWI + 1
        IBD(1) = 3
        IBD(2) = 3
        IBD(3) = 2
        IBD(4) = 2
      ELSE
        IWIX = IWI
        IBD(1) = 2
        IBD(2) = 2
        IBD(3) = 2
        IBD(4) = 2
      ENDIF
C
C     CONVERT HYCOM LON,LAT TO CLIMATOLOGY ARRAY COORDS.
C
      XLIN  = XFIN + (IWIX-1)*DXIN
      XAMIN = 2*IWI
      XAMAX = 0
      DO J= 1,JDM
        DO I= 1,IDM
          XOV = PLON(I,J)
          IF     (XOV.LT.XFIN) THEN
            XOV = XOV + 360.0
          ELSEIF (XOV.GE.XLIN) THEN
            XOV = XOV - 360.0
          ENDIF
C
          XAF(I,J) = 3.0 + (XOV - XFIN)/DXIN
C
          IF     (MOD(J,100).EQ.1 .OR. J.EQ.JDM) THEN
            IF     (MOD(I,10).EQ.1 .OR. I.EQ.IDM) THEN
              WRITE(6,'("I,J,LONV,XAF =",2I5,2F10.3)') I,J,XOV,XAF(I,J)
            ENDIF
          ENDIF
          XAMIN  = MIN( XAMIN, XAF(I,J) )
          XAMAX  = MAX( XAMAX, XAF(I,J) )
        ENDDO
      ENDDO
C
      YAMIN = 2*JWI
      YAMAX = 0
      DO I= 1,IDM
        DO J= 1,JDM
          YOU = PLAT(I,J)
C
          YAF(I,J) = 3.0 + (YOU - YFIN)/DYIN
C
          IF     (MOD(I,100).EQ.1 .OR. I.EQ.IDM) THEN
            IF     (MOD(J,10).EQ.1 .OR. J.EQ.JDM) THEN
              WRITE(6,'("I,J,LATU,YAF =",2I5,2F10.3)') I,J,YOU,YAF(I,J)
            ENDIF
          ENDIF
          YAMIN  = MIN( YAMIN, YAF(I,J) )
          YAMAX  = MAX( YAMAX, YAF(I,J) )
        ENDDO
      ENDDO
C
      WRITE(6,6200) XAMIN,XAMAX,YAMIN,YAMAX
      CALL ZHFLSH(6)
C
C     CHECK THAT THE INTERPOLATION IS 'SAFE',
C
      IF     (INT(XAMIN).LT.3 .OR. INT(XAMAX).GT.IWI+2 .OR.
     +        INT(YAMIN).LT.3 .OR. INT(YAMAX).GT.JWI+2     ) THEN
        WRITE(6,9150)
        CALL ZHFLSH(6)
        STOP
      ENDIF
C
C     PROCESS ALL THE CLIM RECORDS.
C
      DO J= 1,JDM
        DO I= 1,IDM
          RT(I,J) =  ZERO
        ENDDO
      ENDDO
C
C     CONVERSION FACTORS ARE THE SAME FOR T AND S
C
      ADD_OFF = 15.0
      SCALE_F = 0.001
C
      DO 810 KREC= 1,KWI
C
C       READ THE INPUT CLIMS.
C
        call ncheck(nf90_get_var(ncFIDt,ncVIDt,
     &                           TSEAI2(:,:),
     &                           (/ 1,1,KREC /) ))
        call ncheck(nf90_get_var(ncFIDs,ncVIDs,
     &                           SSEAI2(:,:),
     &                           (/ 1,1,KREC /) ))
C
C       COPY INTO THE (LARGER) INTERPOLATION ARRAYS.
C       CONVERT FROM IN-SITU TO POTENTIAL TEMPERATURE.
C
        DEPTH = ZLEV(KREC)
        DO 310 J= 1,JWI
          DO 311 I= 1,IWI
            SSEAI(I+2,J+2) = SSEAI2(I,J)*SCALE_F + ADD_OFF
            TINSIT         = TSEAI2(I,J)*SCALE_F + ADD_OFF
            DBAR           = P80(DEPTH,PLAT(I,J))
            TSEAI(I+2,J+2) = THETA(SSEAI(I+2,J+2),TINSIT, DBAR,0.0)
  311     CONTINUE
  310   CONTINUE
C
C       FILL IN THE PADDING AREA AS NECESSARY.
C
        IF     (INT(XAMAX).GE.IWI+1) THEN
          IF     (IWIX.GT.IWI) THEN
            DO 320 J= 3,JWI+2
              TSEAI(IWI+3,J) = TSEAI(3,J)
              TSEAI(IWI+4,J) = TSEAI(4,J)
              SSEAI(IWI+3,J) = SSEAI(3,J)
              SSEAI(IWI+4,J) = SSEAI(4,J)
  320       CONTINUE
          ELSE
            DO 325 J= 3,JWI+2
              TSEAI(IWI+3,J) = 2.0*TSEAI(IWI+2,J) -     TSEAI(IWI+1,J)
              TSEAI(IWI+4,J) = 3.0*TSEAI(IWI+2,J) - 2.0*TSEAI(IWI+1,J)
              SSEAI(IWI+3,J) = 2.0*SSEAI(IWI+2,J) -     SSEAI(IWI+1,J)
              SSEAI(IWI+4,J) = 3.0*SSEAI(IWI+2,J) - 2.0*SSEAI(IWI+1,J)
  325       CONTINUE
          ENDIF
        ENDIF
        IF     (INT(XAMIN).LE.3) THEN
          IF     (IWIX.GT.IWI) THEN
            DO 330 J= 3,JWI+2
              TSEAI(1,J) = TSEAI(IWI+1,J)
              TSEAI(2,J) = TSEAI(IWI+2,J)
              SSEAI(1,J) = SSEAI(IWI+1,J)
              SSEAI(2,J) = SSEAI(IWI+2,J)
  330       CONTINUE
          ELSE
            DO 335 J= 3,JWI+2
              TSEAI(1,J) = 3.0*TSEAI(3,J) - 2.0*TSEAI(4,J)
              TSEAI(2,J) = 2.0*TSEAI(3,J) -     TSEAI(4,J)
              SSEAI(1,J) = 3.0*SSEAI(3,J) - 2.0*SSEAI(4,J)
              SSEAI(2,J) = 2.0*SSEAI(3,J) -     SSEAI(4,J)
  335       CONTINUE
          ENDIF
        ENDIF
        IF     (INT(YAMAX).GE.JWI+1) THEN
          DO 340 I= 1,IWI+4
            TSEAI(I,JWI+3) = 2.0*TSEAI(I,JWI+2) -     TSEAI(I,JWI+1)
            TSEAI(I,JWI+4) = 3.0*TSEAI(I,JWI+2) - 2.0*TSEAI(I,JWI+1)
            SSEAI(I,JWI+3) = 2.0*SSEAI(I,JWI+2) -     SSEAI(I,JWI+1)
            SSEAI(I,JWI+4) = 3.0*SSEAI(I,JWI+2) - 2.0*SSEAI(I,JWI+1)
  340     CONTINUE
        ENDIF
        IF     (INT(YAMIN).LE.3) THEN
          DO 350 I= 1,IWI+4
            TSEAI(I,1) = 3.0*TSEAI(I,3) - 2.0*TSEAI(I,4)
            TSEAI(I,2) = 2.0*TSEAI(I,3) -     TSEAI(I,4)
            SSEAI(I,1) = 3.0*SSEAI(I,3) - 2.0*SSEAI(I,4)
            SSEAI(I,2) = 2.0*SSEAI(I,3) -     SSEAI(I,4)
  350     CONTINUE
        ENDIF
C
C       INTERPOLATE FROM NATIVE TO MODEL CLIM GRIDS.
C       ALSO INFORCE A STABLE DENSITY PROFILE.
C
        IF     (ICTYPE.EQ.3) THEN
          IF     (INTERP.EQ.0) THEN
            CALL LINEAR(TM,XAF,YAF,IDM,IDM,JDM,
     +                  TSEAI,IWI+4,IWI+4,JWI+4)
            CALL LINEAR(SM,XAF,YAF,IDM,IDM,JDM,
     +                  SSEAI,IWI+4,IWI+4,JWI+4)
          ELSE
            CALL CUBSPL(TM,XAF,YAF,IDM,IDM,JDM,
     +                  TSEAI,IWI+4,IWI+4,JWI+4, IBD, FXI,FYI,WQSEA3,WK)
            CALL CUBSPL(SM,XAF,YAF,IDM,IDM,JDM,
     +                  SSEAI,IWI+4,IWI+4,JWI+4, IBD, FXI,FYI,WQSEA3,WK)
          ENDIF
C
C         ASSUME ICE FORMS (I.E. MIN SST) AT -1.8 DEGC.
C         USE SURFACE FREEZING POINT FOR ALL DEPTHS.
C
          DO J= 1,JDM
            DO I= 1,IDM
              IF     (TM(I,J).LT.-1.8) THEN
                TM(I,J) = -1.8
              ENDIF
            ENDDO
          ENDDO
        ELSEIF (ICTYPE.EQ.4) THEN
          IF     (INTERP.EQ.0) THEN
            CALL LINEAR(TM,XAF,YAF,IDM,IDM,JDM,
     +                  TSEAI,IWI+4,IWI+4,JWI+4)
            CALL LINEAR(SM,XAF,YAF,IDM,IDM,JDM,
     +                  SSEAI,IWI+4,IWI+4,JWI+4)
          ELSE
            CALL CUBSPL(TM,XAF,YAF,IDM,IDM,JDM,
     +                  TSEAI,IWI+4,IWI+4,JWI+4, IBD, FXI,FYI,WQSEA3,WK)
            CALL CUBSPL(SM,XAF,YAF,IDM,IDM,JDM,
     +                  SSEAI,IWI+4,IWI+4,JWI+4, IBD, FXI,FYI,WQSEA3,WK)
          ENDIF
C
C         USE SURFACE FREEZING POINT AS MINUMUM FOR ALL DEPTHS.
C
          DO J= 1,JDM
            DO I= 1,IDM
              TM(I,J) = MAX( TM(I,J), -0.054*SM(I,J) )
            ENDDO
          ENDDO
        ENDIF
        IF     (KSIGMA.EQ.0) THEN
          IF     (SIGVER.EQ.1) THEN
            CALL SIGMA_1(RM,TM,SM,RT,IDM,JDM, ICTYPE)
          ELSEIF (SIGVER.EQ.3) THEN
            CALL SIGMA_3(RM,TM,SM,RT,IDM,JDM, ICTYPE)
          ELSEIF (SIGVER.EQ.5) THEN
            CALL SIGMA_5(RM,TM,SM,RT,IDM,JDM, ICTYPE)
          ELSEIF (SIGVER.EQ.7) THEN
            CALL SIGMA_7(RM,TM,SM,RT,IDM,JDM, ICTYPE)
          ENDIF
        ELSEIF (KSIGMA.EQ.2) THEN
          IF     (SIGVER.EQ.2) THEN
            CALL SIGMA_2(RM,TM,SM,RT,IDM,JDM, ICTYPE)
          ELSEIF (SIGVER.EQ.4) THEN
            CALL SIGMA_4(RM,TM,SM,RT,IDM,JDM, ICTYPE)
          ELSEIF (SIGVER.EQ.6) THEN
            CALL SIGMA_6(RM,TM,SM,RT,IDM,JDM, ICTYPE)
          ELSEIF (SIGVER.EQ.8) THEN
            CALL SIGMA_8(RM,TM,SM,RT,IDM,JDM, ICTYPE)
          ENDIF
        ELSE
          IF     (SIGVER.EQ.46) THEN
            CALL SIGMA_46(RM,TM,SM,RT,IDM,JDM, ICTYPE)
          ELSEIF (SIGVER.EQ.48) THEN
            CALL SIGMA_48(RM,TM,SM,RT,IDM,JDM, ICTYPE)
          ENDIF
        ENDIF
C
C       WRITE OUT STATISTICS.
C
        CALL MINMAX(TM,IDM,JDM, XMIN,XMAX)
        CALL AVERMS(TM,IDM,JDM, XAVE,XRMS)
        WRITE(6,8100) 'TSEA', XMIN,XMAX,XAVE,XRMS
C
        CALL MINMAX(SM,IDM,JDM, XMIN,XMAX)
        CALL AVERMS(SM,IDM,JDM, XAVE,XRMS)
        WRITE(6,8100) 'SSEA', XMIN,XMAX,XAVE,XRMS
C
        CALL MINMAX(RM,IDM,JDM, XMIN,XMAX)
        CALL AVERMS(RM,IDM,JDM, XAVE,XRMS)
        WRITE(6,8100) 'RSEA', XMIN,XMAX,XAVE,XRMS
C
C       DIAGNOSTIC PRINTOUT.
C
        IF     (MIN(ITEST,JTEST).GT.0) THEN
          WRITE(6,*)
          WRITE(6,'(A,2I5,I3,A,F8.2,A,3F6.2)')
     +      'I,J,K =',ITEST,JTEST,KREC,
     +      '   ZLEV =',ZLEV(KREC),
     +      '   R,T,S =',RM(ITEST,JTEST),
     +                   TM(ITEST,JTEST),
     +                   SM(ITEST,JTEST)
          WRITE(6,*)
          CALL ZHFLSH(6)
        ENDIF
C
C       WRITE OUT HYCOM CLIMS.
C
        CALL ZAIOWR(TM,MSK,.FALSE., XMIN,XMAX, 10, .FALSE.)
        WRITE(10,4102) 'potential temperature',ZLEV(KREC),XMIN,XMAX
C
        CALL ZAIOWR(SM,MSK,.FALSE., XMIN,XMAX, 11, .FALSE.)
        WRITE(11,4102) '             salinity',ZLEV(KREC),XMIN,XMAX
C
        CALL ZAIOWR(RM,MSK,.FALSE., XMIN,XMAX, 12, .FALSE.)
        IF     (KSIGMA.EQ.0) THEN
          WRITE(12,4102) '              sigma-0',ZLEV(KREC),XMIN,XMAX
        ELSEIF (KSIGMA.EQ.2) THEN
          WRITE(12,4102) '              sigma-2',ZLEV(KREC),XMIN,XMAX
        ELSE
          WRITE(12,4102) '              sigma-4',ZLEV(KREC),XMIN,XMAX
        ENDIF
C
        WRITE(6,6300) KREC,ZLEV(KREC)
        CALL ZHFLSH(6)
  810 CONTINUE
C
      CALL ZAIOCL(10)
      CLOSE( UNIT=10)
      CALL ZAIOCL(11)
      CLOSE( UNIT=11)
      CALL ZAIOCL(12)
      CLOSE( UNIT=12)
C
C     SUMMARY.
C
      WRITE(6,6400) KWI
      CALL ZHFLSH(6)
      STOP
C
 4101 FORMAT(A79)
 4102 FORMAT(A,': depth,range = ',F7.1,1P2E16.7)
 6000 FORMAT(1X,A,2X,A40 //)
 6200 FORMAT(/ 1X,'MIN,MAX I COORDS = ',F8.2,',',F8.2 
     +       / 1X,'MIN,MAX J COORDS = ',F8.2,',',F8.2 /)
 6300 FORMAT(10X,'WRITING CLIM RECORD',I3,'     ZLEV =',F7.1 /)
 6400 FORMAT(I5,' LEVEL CLIMATOLOGY COMPLETED.')
 8100 FORMAT(1X,A,': MIN=',F13.8,' MAX=',F13.8,
     +             ' AVE=',F13.8,' RMS=',F13.8)
 9000 FORMAT(// 20X,'*****  ERROR ON UNIT -',I3,
     +   ' INPUT DOES NOT AGREE WITH PARAMETERS  *****' //
     +   1X,'IWI,JWI,KWI = ',I5,  I10,  I4 /
     +   1X,'XFIN,YFIN   = ',F8.2,F10.2    /
     +   1X,'DXIN,DYIN   = ',F9.3, F9.3    //)
 9050 FORMAT(// 20X,'**********  ERROR - ',
     +   'INPUT IS GLOBAL AND IWI*DXIN IS NOT 360 DEGREES  ****' //)
 9150 FORMAT(// 20X,'**********  ERROR - ',
     +   'THE OUTPUT GRID IS NOT INSIDE THE INPUT GRID  **********' //)
C     END OF PROGRAM WNDINT.
      END
      SUBROUTINE SIGMA_1(RSEA,TSEA,SSEA,RTOP,IW,JW, ICTYPE)
      IMPLICIT NONE
C
      INTEGER IW,JW,ICTYPE
      REAL*4  TSEA(IW,JW),SSEA(IW,JW),RSEA(IW,JW),RTOP(IW,JW)
C
C     CALCULATE THE DIAGNOSTIC THERMODYNAMEIC FIELD.
C
      INTEGER I,J
C
      INCLUDE '../../include/stmt_fns_SIGMA0_7term.h'
C
      IF     (ICTYPE.EQ.3) THEN
        DO J= 1,JW
          DO I= 1,IW
            RSEA(I,J) =    SIG( R8(TSEA(I,J)), R8(SSEA(I,J)) )
            RSEA(I,J) = MAX(       RSEA(I,J),     RTOP(I,J) )
*           TSEA(I,J) = TOFSIG( R8(RSEA(I,J)), R8(SSEA(I,J)) )
            SSEA(I,J) = SOFSIG( R8(RSEA(I,J)), R8(TSEA(I,J)) )
            RTOP(I,J) =            RSEA(I,J)
          ENDDO
        ENDDO
      ELSEIF (ICTYPE.EQ.4) THEN
        DO J= 1,JW
          DO I= 1,IW
            RSEA(I,J) = SIG( R8(TSEA(I,J)), R8(SSEA(I,J)) )
            RTOP(I,J) = RSEA(I,J)
          ENDDO
        ENDDO
      ENDIF
      RETURN
      END
      SUBROUTINE SIGMA_2(RSEA,TSEA,SSEA,RTOP,IW,JW, ICTYPE)
      IMPLICIT NONE
C
      INTEGER IW,JW,ICTYPE
      REAL*4  TSEA(IW,JW),SSEA(IW,JW),RSEA(IW,JW),RTOP(IW,JW)
C
C     CALCULATE THE DIAGNOSTIC THERMODYNAMEIC FIELD.
C
      INTEGER I,J
C
      INCLUDE '../../include/stmt_fns_SIGMA2_7term.h'
C
      IF     (ICTYPE.EQ.3) THEN
        DO J= 1,JW
          DO I= 1,IW
            RSEA(I,J) =    SIG( R8(TSEA(I,J)), R8(SSEA(I,J)) )
            RSEA(I,J) = MAX(       RSEA(I,J),     RTOP(I,J) )
*           TSEA(I,J) = TOFSIG( R8(RSEA(I,J)), R8(SSEA(I,J)) )
            SSEA(I,J) = SOFSIG( R8(RSEA(I,J)), R8(TSEA(I,J)) )
            RTOP(I,J) =            RSEA(I,J)
          ENDDO
        ENDDO
      ELSEIF (ICTYPE.EQ.4) THEN
        DO J= 1,JW
          DO I= 1,IW
            RSEA(I,J) = SIG( R8(TSEA(I,J)), R8(SSEA(I,J)) )
            RTOP(I,J) = RSEA(I,J)
          ENDDO
        ENDDO
      ENDIF
      RETURN
      END
      SUBROUTINE SIGMA_3(RSEA,TSEA,SSEA,RTOP,IW,JW, ICTYPE)
      IMPLICIT NONE
C
      INTEGER IW,JW,ICTYPE
      REAL*4  TSEA(IW,JW),SSEA(IW,JW),RSEA(IW,JW),RTOP(IW,JW)
C
C     CALCULATE THE DIAGNOSTIC THERMODYNAMEIC FIELD.
C
      INTEGER I,J
C
      INCLUDE '../../include/stmt_fns_SIGMA0_9term.h'
C
      IF     (ICTYPE.EQ.3) THEN
        DO J= 1,JW
          DO I= 1,IW
            RSEA(I,J) =    SIG( R8(TSEA(I,J)), R8(SSEA(I,J)) )
            RSEA(I,J) = MAX(       RSEA(I,J),     RTOP(I,J) )
*           TSEA(I,J) = TOFSIG( R8(RSEA(I,J)), R8(SSEA(I,J)) )
            SSEA(I,J) = SOFSIG( R8(RSEA(I,J)), R8(TSEA(I,J)) )
            RTOP(I,J) =            RSEA(I,J)
          ENDDO
        ENDDO
      ELSEIF (ICTYPE.EQ.4) THEN
        DO J= 1,JW
          DO I= 1,IW
            RSEA(I,J) = SIG( R8(TSEA(I,J)), R8(SSEA(I,J)) )
            RTOP(I,J) = RSEA(I,J)
          ENDDO
        ENDDO
      ENDIF
      RETURN
      END
      SUBROUTINE SIGMA_4(RSEA,TSEA,SSEA,RTOP,IW,JW, ICTYPE)
      IMPLICIT NONE
C
      INTEGER IW,JW,ICTYPE
      REAL*4  TSEA(IW,JW),SSEA(IW,JW),RSEA(IW,JW),RTOP(IW,JW)
C
C     CALCULATE THE DIAGNOSTIC THERMODYNAMEIC FIELD.
C
      INTEGER I,J
C
      INCLUDE '../../include/stmt_fns_SIGMA2_9term.h'
C
      IF     (ICTYPE.EQ.3) THEN
        DO J= 1,JW
          DO I= 1,IW
            RSEA(I,J) =    SIG( R8(TSEA(I,J)), R8(SSEA(I,J)) )
            RSEA(I,J) = MAX(       RSEA(I,J),     RTOP(I,J) )
*           TSEA(I,J) = TOFSIG( R8(RSEA(I,J)), R8(SSEA(I,J)) )
            SSEA(I,J) = SOFSIG( R8(RSEA(I,J)), R8(TSEA(I,J)) )
            RTOP(I,J) =            RSEA(I,J)
          ENDDO
        ENDDO
      ELSEIF (ICTYPE.EQ.4) THEN
        DO J= 1,JW
          DO I= 1,IW
            RSEA(I,J) = SIG( R8(TSEA(I,J)), R8(SSEA(I,J)) )
            RTOP(I,J) = RSEA(I,J)
          ENDDO
        ENDDO
      ENDIF
      RETURN
      END
      SUBROUTINE SIGMA_5(RSEA,TSEA,SSEA,RTOP,IW,JW, ICTYPE)
      IMPLICIT NONE
C
      INTEGER IW,JW,ICTYPE
      REAL*4  TSEA(IW,JW),SSEA(IW,JW),RSEA(IW,JW),RTOP(IW,JW)
C
C     CALCULATE THE DIAGNOSTIC THERMODYNAMEIC FIELD.
C
      INTEGER I,J,NN
      REAL*8  RR,SN,SO,TT
C
      REAL*8, PARAMETER :: TOL=1.D-6
C
      INCLUDE '../../include/stmt_fns_SIGMA0_17term.h'
C
      IF     (ICTYPE.EQ.3) THEN
        DO J= 1,JW
          DO I= 1,IW
            RSEA(I,J) =    SIG( R8(TSEA(I,J)), R8(SSEA(I,J)) )
            RSEA(I,J) = MAX(       RSEA(I,J),     RTOP(I,J) )
            RTOP(I,J) =            RSEA(I,J)
          ENDDO
        ENDDO
C       sofsig via Newton iteration from a 12-term 1st guess
        CALL SOFSIG_7(RSEA,TSEA,SSEA,IW,JW)
        DO J= 1,JW
          DO I= 1,IW
            SN = R8(SSEA(I,J))  !non-negative
            TT = R8(TSEA(I,J))
            RR = R8(RSEA(I,J))
            DO NN= 1,10
              SO = SN
              SN = SO - (SIG(TT,SO)-RR)/DSIGDS(TT,SO)
              IF     (NN.EQ.10 .OR. ABS(SN-SO).LT.TOL) THEN
                EXIT
              ENDIF
            ENDDO !nn
            SSEA(I,J) = SN
          ENDDO
        ENDDO
      ELSEIF (ICTYPE.EQ.4) THEN
        DO J= 1,JW
          DO I= 1,IW
            RSEA(I,J) = SIG( R8(TSEA(I,J)), R8(SSEA(I,J)) )
            RTOP(I,J) = RSEA(I,J)
          ENDDO
        ENDDO
      ENDIF
      RETURN
      END
      SUBROUTINE SIGMA_6(RSEA,TSEA,SSEA,RTOP,IW,JW, ICTYPE)
      IMPLICIT NONE
C
      INTEGER IW,JW,ICTYPE
      REAL*4  TSEA(IW,JW),SSEA(IW,JW),RSEA(IW,JW),RTOP(IW,JW)
C
C     CALCULATE THE DIAGNOSTIC THERMODYNAMEIC FIELD.
C
      INTEGER I,J,NN
      REAL*8  RR,SN,SO,TT
C
      REAL*8, PARAMETER :: TOL=1.D-6
C
      INCLUDE '../../include/stmt_fns_SIGMA2_17term.h'
C
      IF     (ICTYPE.EQ.3) THEN
        DO J= 1,JW
          DO I= 1,IW
            RSEA(I,J) =    SIG( R8(TSEA(I,J)), R8(SSEA(I,J)) )
            RSEA(I,J) = MAX(       RSEA(I,J),     RTOP(I,J) )
            RTOP(I,J) =            RSEA(I,J)
          ENDDO
        ENDDO
C       sofsig via Newton iteration from a 12-term 1st guess
        CALL SOFSIG_8(RSEA,TSEA,SSEA,IW,JW)
        DO J= 1,JW
          DO I= 1,IW
            SN = R8(SSEA(I,J))  !non-negative
            TT = R8(TSEA(I,J))
            RR = R8(RSEA(I,J))
            DO NN= 1,10
              SO = SN
              SN = SO - (SIG(TT,SO)-RR)/DSIGDS(TT,SO)
              IF     (NN.EQ.10 .OR. ABS(SN-SO).LT.TOL) THEN
                EXIT
              ENDIF
            ENDDO !nn
            SSEA(I,J) = SN
          ENDDO
        ENDDO
      ELSEIF (ICTYPE.EQ.4) THEN
        DO J= 1,JW
          DO I= 1,IW
            RSEA(I,J) = SIG( R8(TSEA(I,J)), R8(SSEA(I,J)) )
            RTOP(I,J) = RSEA(I,J)
          ENDDO
        ENDDO
      ENDIF
      RETURN
      END
      SUBROUTINE SIGMA_46(RSEA,TSEA,SSEA,RTOP,IW,JW, ICTYPE)
      IMPLICIT NONE
C
      INTEGER IW,JW,ICTYPE
      REAL*4  TSEA(IW,JW),SSEA(IW,JW),RSEA(IW,JW),RTOP(IW,JW)
C
C     CALCULATE THE DIAGNOSTIC THERMODYNAMEIC FIELD.
C
      INTEGER I,J,NN
      REAL*8  RR,SN,SO,TT
C
      REAL*8, PARAMETER :: TOL=1.D-6
C
      INCLUDE '../../include/stmt_fns_SIGMA4_17term.h'
C
      IF     (ICTYPE.EQ.3) THEN
        DO J= 1,JW
          DO I= 1,IW
            RSEA(I,J) =    SIG( R8(TSEA(I,J)), R8(SSEA(I,J)) )
            RSEA(I,J) = MAX(       RSEA(I,J),     RTOP(I,J) )
            RTOP(I,J) =            RSEA(I,J)
          ENDDO
        ENDDO
C       sofsig via Newton iteration from a 12-term 1st guess
        CALL SOFSIG_48(RSEA,TSEA,SSEA,IW,JW)
        DO J= 1,JW
          DO I= 1,IW
            SN = R8(SSEA(I,J))  !non-negative
            TT = R8(TSEA(I,J))
            RR = R8(RSEA(I,J))
            DO NN= 1,10
              SO = SN
              SN = SO - (SIG(TT,SO)-RR)/DSIGDS(TT,SO)
              IF     (NN.EQ.10 .OR. ABS(SN-SO).LT.TOL) THEN
                EXIT
              ENDIF
            ENDDO !nn
            SSEA(I,J) = SN
          ENDDO
        ENDDO
      ELSEIF (ICTYPE.EQ.4) THEN
        DO J= 1,JW
          DO I= 1,IW
            RSEA(I,J) = SIG( R8(TSEA(I,J)), R8(SSEA(I,J)) )
            RTOP(I,J) = RSEA(I,J)
          ENDDO
        ENDDO
      ENDIF
      RETURN
      END
      SUBROUTINE SIGMA_7(RSEA,TSEA,SSEA,RTOP,IW,JW, ICTYPE)
      IMPLICIT NONE
C
      INTEGER IW,JW,ICTYPE
      REAL*4  TSEA(IW,JW),SSEA(IW,JW),RSEA(IW,JW),RTOP(IW,JW)
C
C     CALCULATE THE DIAGNOSTIC THERMODYNAMEIC FIELD.
C
      INTEGER I,J
C
      INCLUDE '../../include/stmt_fns_SIGMA0_12term.h'
C
      IF     (ICTYPE.EQ.3) THEN
        DO J= 1,JW
          DO I= 1,IW
            RSEA(I,J) =    SIG( R8(TSEA(I,J)), R8(SSEA(I,J)) )
            RSEA(I,J) = MAX(       RSEA(I,J),     RTOP(I,J) )
*           TSEA(I,J) = TOFSIG( R8(RSEA(I,J)), R8(SSEA(I,J)) )
            SSEA(I,J) = SOFSIG( R8(RSEA(I,J)), R8(TSEA(I,J)) )
            RTOP(I,J) =            RSEA(I,J)
          ENDDO
        ENDDO
      ELSEIF (ICTYPE.EQ.4) THEN
        DO J= 1,JW
          DO I= 1,IW
            RSEA(I,J) = SIG( R8(TSEA(I,J)), R8(SSEA(I,J)) )
            RTOP(I,J) = RSEA(I,J)
          ENDDO
        ENDDO
      ENDIF
      RETURN
      END
      SUBROUTINE SIGMA_8(RSEA,TSEA,SSEA,RTOP,IW,JW, ICTYPE)
      IMPLICIT NONE
C
      INTEGER IW,JW,ICTYPE
      REAL*4  TSEA(IW,JW),SSEA(IW,JW),RSEA(IW,JW),RTOP(IW,JW)
C
C     CALCULATE THE DIAGNOSTIC THERMODYNAMEIC FIELD.
C
      INTEGER I,J
C
      INCLUDE '../../include/stmt_fns_SIGMA2_12term.h'
C
      IF     (ICTYPE.EQ.3) THEN
        DO J= 1,JW
          DO I= 1,IW
            RSEA(I,J) =    SIG( R8(TSEA(I,J)), R8(SSEA(I,J)) )
            RSEA(I,J) = MAX(       RSEA(I,J),     RTOP(I,J) )
*           TSEA(I,J) = TOFSIG( R8(RSEA(I,J)), R8(SSEA(I,J)) )
            SSEA(I,J) = SOFSIG( R8(RSEA(I,J)), R8(TSEA(I,J)) )
            RTOP(I,J) =            RSEA(I,J)
          ENDDO
        ENDDO
      ELSEIF (ICTYPE.EQ.4) THEN
        DO J= 1,JW
          DO I= 1,IW
            RSEA(I,J) = SIG( R8(TSEA(I,J)), R8(SSEA(I,J)) )
            RTOP(I,J) = RSEA(I,J)
          ENDDO
        ENDDO
      ENDIF
      RETURN
      END
      SUBROUTINE SIGMA_48(RSEA,TSEA,SSEA,RTOP,IW,JW, ICTYPE)
      IMPLICIT NONE
C
      INTEGER IW,JW,ICTYPE
      REAL*4  TSEA(IW,JW),SSEA(IW,JW),RSEA(IW,JW),RTOP(IW,JW)
C
C     CALCULATE THE DIAGNOSTIC THERMODYNAMEIC FIELD.
C
      INTEGER I,J
C
      INCLUDE '../../include/stmt_fns_SIGMA4_12term.h'
C
      IF     (ICTYPE.EQ.3) THEN
        DO J= 1,JW
          DO I= 1,IW
            RSEA(I,J) =    SIG( R8(TSEA(I,J)), R8(SSEA(I,J)) )
            RSEA(I,J) = MAX(       RSEA(I,J),     RTOP(I,J) )
*           TSEA(I,J) = TOFSIG( R8(RSEA(I,J)), R8(SSEA(I,J)) )
            SSEA(I,J) = SOFSIG( R8(RSEA(I,J)), R8(TSEA(I,J)) )
            RTOP(I,J) =            RSEA(I,J)
          ENDDO
        ENDDO
      ELSEIF (ICTYPE.EQ.4) THEN
        DO J= 1,JW
          DO I= 1,IW
            RSEA(I,J) = SIG( R8(TSEA(I,J)), R8(SSEA(I,J)) )
            RTOP(I,J) = RSEA(I,J)
          ENDDO
        ENDDO
      ENDIF
      RETURN
      END
      SUBROUTINE SOFSIG_7(RSEA,TSEA,SSEA,IW,JW)
      IMPLICIT NONE
C
      INTEGER IW,JW
      REAL*4  TSEA(IW,JW),SSEA(IW,JW),RSEA(IW,JW)
C
C     CALCULATE SALINITY FROM POT.TEMPERATURE AND POT.DENSITY.
C
      INTEGER I,J
C
      INCLUDE '../../include/stmt_fns_SIGMA0_12term.h'
C
      DO J= 1,JW
        DO I= 1,IW
          SSEA(I,J) = SOFSIG( R8(RSEA(I,J)), R8(TSEA(I,J)) )
        ENDDO
      ENDDO
      RETURN
      END
      SUBROUTINE SOFSIG_8(RSEA,TSEA,SSEA,IW,JW)
      IMPLICIT NONE
C
      INTEGER IW,JW
      REAL*4  TSEA(IW,JW),SSEA(IW,JW),RSEA(IW,JW)
C
C     CALCULATE SALINITY FROM POT.TEMPERATURE AND POT.DENSITY.
C
      INTEGER I,J
C
      INCLUDE '../../include/stmt_fns_SIGMA2_12term.h'
C
      DO J= 1,JW
        DO I= 1,IW
          SSEA(I,J) = SOFSIG( R8(RSEA(I,J)), R8(TSEA(I,J)) )
        ENDDO
      ENDDO
      RETURN
      END
      SUBROUTINE SOFSIG_48(RSEA,TSEA,SSEA,IW,JW)
      IMPLICIT NONE
C
      INTEGER IW,JW
      REAL*4  TSEA(IW,JW),SSEA(IW,JW),RSEA(IW,JW)
C
C     CALCULATE SALINITY FROM POT.TEMPERATURE AND POT.DENSITY.
C
      INTEGER I,J
C
      INCLUDE '../../include/stmt_fns_SIGMA4_12term.h'
C
      DO J= 1,JW
        DO I= 1,IW
          SSEA(I,J) = SOFSIG( R8(RSEA(I,J)), R8(TSEA(I,J)) )
        ENDDO
      ENDDO
      RETURN
      END

      subroutine ncheck(status)
      use netcdf   ! NetCDF fortran 90 interface
      implicit none
c
      integer, intent(in) :: status
c
c     subroutine to handle NetCDF errors
c
      if (status /= nf90_noerr) then
        write(6,'(/a)')   'error from NetCDF library'
        write(6,'(a/)')   trim(nf90_strerror(status))
        call zhflsh(6)
        stop
      end if
      end subroutine ncheck

c separate set of subroutines, adapted from WHOI CTD group
c     real function atg(s,t,p)
c     real function theta(s,t0,p0,pr)
c      function p80(dpth,xlat)
c
c n fofonoff & r millard
c
      real function atg(s,t,p)
c ****************************
c adiabatic temperature gradient deg c per decibar
c ref: bryden,h.,1973,deep-sea res.,20,401-408
c units:
c       pressure        p        decibars
c       temperature     t        deg celsius (ipts-68)
c       salinity        s        (ipss-78)
c       adiabatic       atg      deg. c/decibar
c checkvalue: atg=3.255976e-4 c/dbar for s=40 (ipss-78),
c t=40 deg c,p0=10000 decibars
      ds = s - 35.0
      atg = (((-2.1687e-16*t+1.8676e-14)*t-4.6206e-13)*p
     x+((2.7759e-12*t-1.1351e-10)*ds+((-5.4481e-14*t
     x+8.733e-12)*t-6.7795e-10)*t+1.8741e-8))*p
     x+(-4.2393e-8*t+1.8932e-6)*ds
     x+((6.6228e-10*t-6.836e-8)*t+8.5258e-6)*t+3.5803e-5
      return
      end

      real function theta(s,t0,p0,pr)
c ***********************************
c to compute local potential temperature at pr
c using bryden 1973 polynomial for adiabatic lapse rate
c and runge-kutta 4-th order integration algorithm.
c ref: bryden,h.,1973,deep-sea res.,20,401-408
c fofonoff,n.,1977,deep-sea res.,24,489-491
c units:
c       pressure        p0       decibars
c       temperature     t0       deg celsius (ipts-68)
c       salinity        s        (ipss-78)
c       reference prs   pr       decibars
c       potential tmp.  theta    deg celsius
c checkvalue: theta= 36.89073 c,s=40 (ipss-78),t0=40 deg c,
c p0=10000 decibars,pr=0 decibars
c
c      set-up intermediate temperature and pressure variables
      p=p0
      t=t0
c**************
      h = pr - p
      xk = h*atg(s,t,p)
      t = t + 0.5*xk
      q = xk
      p = p + 0.5*h
      xk = h*atg(s,t,p)
      t = t + 0.29289322*(xk-q)
      q = 0.58578644*xk + 0.121320344*q
      xk = h*atg(s,t,p)
      t = t + 1.707106781*(xk-q)
      q = 3.414213562*xk - 4.121320344*q
      p = p + 0.5*h
      xk = h*atg(s,t,p)
      theta = t + (xk-2.0*q)/6.0
      return
      end
c
c pressure from depth from saunder's formula with eos80.
c reference: saunders,peter m., practical conversion of pressure
c            to depth., j.p.o. , april 1981.
c r millard
c march 9, 1983
c check value: p80=7500.004 dbars;for lat=30 deg., depth=7321.45 meters
      function p80(dpth,xlat)
      parameter pi=3.141592654
      plat=abs(xlat*pi/180.)
      d=sin(plat)
      c1=5.92e-3+d**2*5.25e-3
      p80=((1-c1)-sqrt(((1-c1)**2)-(8.84e-6*dpth)))/4.42e-6
      return
      end
