#define sigmav
#define term7
#define sig2
      PROGRAM RELAX
      USE MOD_ZA  ! HYCOM array I/O interface
      IMPLICIT NONE
C
C     WIDTH OF ZONAL AVERAGE (0 FOR NO ZONAL AVERAGE).
C
      INTEGER    JDW
      PARAMETER (JDW=0)
C
C     DIAGNOSTIC/DEBUGGING VARIABLES.
C
      INTEGER        ITEST,JTEST
      COMMON/DEBUGI/ ITEST,JTEST
      SAVE  /DEBUGI/
C
C     BLKDAT VARIABLES.
C
      CHARACTER*79 CTITLE
      CHARACTER*26 SIGFMT
      INTEGER      IVERSN,IEXPT,YRFLAG,KDM,LEVTOP,MONTH,
     +             NHYBRD,NSIGMA,THFLAG
      REAL*4       DP00S,DP00,DP00X,DP00F,DS00,DS00X,DS00F,ISOTOP,
     +             SIGMA(99),THBASE,THKMIN,BLK
      LOGICAL      VSIGMA,ISOPYC
C
C     I/O VARIABLES.
C
      CHARACTER PREAMBL(5)*79,CLINE*80
      REAL*4    HMINA,HMINB,HMAXA,HMAXB
C
C     INPUT CLIM ARRAYS.
C
      REAL*4               :: ZLEV(999)
      REAL*4,  ALLOCATABLE :: TZ(:,:,:),RZ(:,:,:)
C
C     OUTPUT ARRAYS
C
      REAL*4,  ALLOCATABLE :: PMIX(:,:),PM(:,:),TM(:,:),SM(:,:),RM(:,:)
C
C     OTHER ARRAYS.
C
      INTEGER, ALLOCATABLE :: MSK(:,:)
      REAL*4,  ALLOCATABLE :: DEPTH(:,:),PLAT(:,:),
     +                        PKM1(:,:),RKM1(:,:),PKM2(:,:),WORK(:,:),
     +                        SIG3D(:,:,:)
      REAL*4               :: DP0K(99),DS0K(99),SIGSCL(99),
     +                        DPCK(0:99),DSCK(0:99),SIGCUM(0:99),
     +                        ZI(0:999),ZIL(0:999)
C
C**********
C*
C 1)  FROM A Z-LEVEL CLIMATOLOGY ON THE HYCOM REGION GRID,
C      CREATE A 'ISOPYCNAL' CLIMATOLOGY SUITABLE FOR INPUT TO HYCOM.
C
C      ONLY FOR USE WITH HYCOM 2.0.00 OR LATER.
C
C 2)  PARAMETERS:
C
C     MODEL GRID SPECIFICATION (W.R.T. PRESSURE GRID):
C
C        IDM    = 1ST DIMENSION OF MAJOR (HYCOM) MODEL ARRAYS
C        JDM    = 2ND DIMENSION OF MAJOR (HYCOM) MODEL ARRAYS
C        KZ     = NUMBER OF Z-LEVELS IN INPUT CLIMATOLOGY
C
C 3)  INPUT:
C        ON UNIT 51:  BATHYMETRY FILE
C        ON UNIT 51A: BATHYMETRY FILE
C        ON UNIT 52:  SPACIALLY VARYING ISOPYCNAL TARGET DENSITY FILE
C        ON UNIT 52A: SPACIALLY VARYING ISOPYCNAL TARGET DENSITY FILE
C        ON UNIT 71:  DENS. Z-LEVEL CLIM FILE
C        ON UNIT 72:  TEMP. Z-LEVEL CLIM FILE
C        ON UNIT 99:  A SUBSET OF blkdat.input FROM TARGET SIMULATION
C     OUTPUT:
C        ON UNIT 10:  TEMP. LAYER CLIM FILE
C        ON UNIT 10A: TEMP. LAYER CLIM FILE
C        ON UNIT 11:  SALN. LAYER CLIM FILE
C        ON UNIT 11A: SALN. LAYER CLIM FILE
C        ON UNIT 12:  INTF. LAYER CLIM FILE
C        ON UNIT 12A: INTF. LAYER CLIM FILE
C        ON UNIT 21:  DUMMY HYCOM ARCHIVE FILE
C        ON UNIT 21A: DUMMY HYCOM ARCHIVE FILE
C
C 4)  ALAN J. WALLCRAFT,  NAVAL RESEARCH LABORATORY, MARCH 2001.
C*
C**********
C
      REAL*4     ZERO,ONE,RADIAN
      PARAMETER (ZERO=0.0, ONE=1.0, RADIAN=57.2957795)
C
      INTEGER I,J,K,K400,KK,KZ,KZL
      REAL*4  DP0KF,DPMS,DS0KF,DSMS,SZK,RZK,TZK,TIME,THK,THIKMN,
     +        DMIN,
     +        PAVE,XAVE,XMAX,XMIN,TMIN,SMIN,PMIN,ZJ,ZZ,Q,
     +        PINTEG,SIGMAA,SIGMAB,THKA,THKB,ZBOT,ZTOP
#if defined(term7)
C
C     STATEMENT FUNCTIONS.
C
C --- sigma-theta as a function of temp (deg c) and salinity (psu)
C --- (7-term polynomial fit to T:[-2:30],S:[30:38])
C
      REAL*8     C1,C2,C3,C4,C5,C6,C7
#if defined(sig0)
C --- coefficients for sigma-0 (based on Brydon & Sun fit)
      PARAMETER (C1=-1.36471E-01, C2= 4.68181E-02, C3= 8.07004E-01,
     .           C4=-7.45353E-03, C5=-2.94418E-03,
     .           C6= 3.43570E-05, C7= 3.48658E-05)
#elif defined(sig2)
C --- coefficients for sigma-2 (based on Brydon & Sun fit)
      PARAMETER (C1= 9.77093E+00, C2=-2.26493E-02, C3= 7.89879E-01,
     .           C4=-6.43205E-03, C5=-2.62983E-03,
     .           C6= 2.75835E-05, C7= 3.15235E-05)
#endif
C
      REAL*8     AZERO,AHALF,ATHIRD,A1P5
      PARAMETER (AZERO  =0.D0)
      PARAMETER (AHALF  =1.D0/2.D0)
      PARAMETER (ATHIRD =1.D0/3.D0)
      PARAMETER (A1P5   =3.0D0/2.D0)
C
      REAL*4  R4
      REAL*8  R8
      REAL*8  R,T,S
      REAL*8  TOFSIG,SOFSIG,SIG
C
      REAL*8 A0,A1,A2,CUBR,CUBQ,CUBAN,CUBRL,CUBIM
C
C --- auxiliary statement for real*4 to real*8 conversion
      R8(R4)=R4
C
C --- auxiliary statements for finding root of 3rd degree polynomial
      A0(S)=(C1+C3*S)/C6
      A1(S)=(C2+C5*S)/C6
      A2(S)=(C4+C7*S)/C6
      CUBQ(S)=ATHIRD*A1(S)-(ATHIRD*A2(S))**2
      CUBR(R,S)=ATHIRD*(AHALF*A1(S)*A2(S)-A1P5*(A0(S)-R/C6))
     .   -(ATHIRD*A2(S))**3
C --- if q**3+r**2>0, water is too dense to yield real root at given
C --- salinitiy. setting q**3+r**2=0 in that case is equivalent to
C --- lowering sigma until a double real root is obtained.
      CUBAN(R,S)=ATHIRD*ATAN2(SQRT(MAX(AZERO,
     .   -(CUBQ(S)**3+CUBR(R,S)**2))),CUBR(R,S))
      CUBRL(R,S)=SQRT(-CUBQ(S))*COS(CUBAN(R,S))
      CUBIM(R,S)=SQRT(-CUBQ(S))*SIN(CUBAN(R,S))
C
C --- temp (deg c) as a function of sigma and salinity (mil)
      TOFSIG(R,S)=-CUBRL(R,S)+SQRT(3.D0)*CUBIM(R,S)-ATHIRD*A2(S)
C
C --- salinity (mil) as a function of sigma and temperature (deg c)
      SOFSIG(R,T)=(R-C1-T*(C2+T*(C4+C6*T)))/(C3+T*(C5+C7*T))
C
C --- sigma-theta as a function of temp (deg c) and salinity (mil)
C --- (polynomial fit that is cubic in T and linear in S)
      SIG(T,S)=(C1+C3*S+T*(C2+C5*S+T*(C4+C7*S+C6*T)))
#elif defined(term9)
C
C     STATEMENT FUNCTIONS.
C
C --- sigma-theta as a function of temp (deg c) and salinity (psu)
C --- (9-term polynomial fit to T:[-2:30],S:[18:38])
C
#if defined(sig0)
C --- coefficients for sigma-0.
      REAL*8, PARAMETER ::
     &   C1=-4.311829E-02,  !const. coefficent
     &   C2= 5.429948E-02,  !T      coefficent
     &   C3= 8.011774E-01,  !   S   coefficent
     &   C4=-7.641336E-03,  !T^2    coefficent
     &   C5=-3.258442E-03,  !T  S   coefficent
     &   C6= 3.757643E-05,  !T^3    coefficent
     &  RC6=1.0/C6,
     &   C7= 3.630361E-05,  !T^2S   coefficent
     &   C8= 8.675546E-05,  !   S^2 coefficent
     &   C9= 3.995086E-06   !T  S^2 coefficent
#elif defined(sig2)
C --- coefficients for sigma-2.
      REAL*8, PARAMETER ::
     &   C1= 9.903308E+00,  !const. coefficent
     &   C2=-1.618075E-02,  !T      coefficent
     &   C3= 7.819166E-01,  !   S   coefficent
     &   C4=-6.593939E-03,  !T^2    coefficent
     &   C5=-2.896464E-03,  !T  S   coefficent
     &   C6= 3.038697E-05,  !T^3    coefficent
     &  RC6= 1.0/C6,
     &   C7= 3.266933E-05,  !T^2S   coefficent
     &   C8= 1.180109E-04,  !   S^2 coefficent
     &   C9= 3.399511E-06   !T  S^2 coefficent
#endif
C
      REAL*8     AZERO,A3RD
      PARAMETER (AZERO=0.D0,A3RD=1.D0/3.D0)
C
      REAL*4  R4
      REAL*8  R8
      REAL*8  R,S,T
      REAL*8  A0,A1,A2,CUBQ,CUBR,CUBAN,CUBRL,CUBIM,TOFSIG,SOFSIG,SIG
      REAL*8  SOFSIG_A,SOFSIG_B,SOFSIG_C
C
C --- auxiliary statement for real*4 to real*8 conversion
      R8(R4)=R4
C
C --- auxiliary statements for finding root of 3rd degree polynomial
      A0(S,R)=(C1+S*(C3+S*C8)-R)*RC6  !constant  coefficient
      A1(S)  =(C2+S*(C5+S*C9)  )*RC6  !linear    coefficient
      A2(S)  =(C4+S* C7        )*RC6  !quadratic coefficient
                                      !cubic     coefficient is c6*rc6=1.0
      CUBQ(S)  =A3RD*     A1(S)                   -(A3RD*A2(S))**2
      CUBR(R,S)=A3RD*(0.5*A1(S)*A2(S)-1.5*A0(S,R))-(A3RD*A2(S))**3
C --- if q**3+r**2>0, water is too dense to yield real root at given
C --- salinitiy. setting q**3+r**2=0 in that case is equivalent to
C --- lowering sigma until a double real root is obtained.
      CUBAN(R,S)=A3RD*ATAN2(SQRT(MAX(0.0,-(CUBQ(S)**3+CUBR(R,S)**2))),
     &                      CUBR(R,S))
      CUBRL(R,S)=SQRT(-CUBQ(S))*COS(CUBAN(R,S))
      CUBIM(R,S)=SQRT(-CUBQ(S))*SIN(CUBAN(R,S))
C
C --- temp (deg c) as a function of sigma and salinity (psu)
C --- find a cubic polynominal root of t**3+a2*t**2+a1*t+a0=0
      TOFSIG(R,S)=-CUBRL(R,S)+SQRT(3.0)*CUBIM(R,S)-A3RD*A2(S)
C
C --- salinity (psu) as a function of sigma and temperature (deg c)
C --- find a quadratic polynominal root of a*s**2+b*s+c=0
      SOFSIG_A(R,T)=(C8+T* C9)                !quadratic coefficient
      SOFSIG_B(R,T)=(C3+T*(C5+T* C7))         !linear    coefficient
      SOFSIG_C(R,T)=(C1+T*(C2+T*(C4+T*C6))-R) !constant  coefficient
      SOFSIG(R,T)=(2.0*SOFSIG_C(R,T))/
     &            (-SOFSIG_B(R,T)
     &             -SIGN(SQRT(SOFSIG_B(R,T)**2-
     &                        4.0*SOFSIG_A(R,T)*SOFSIG_C(R,T)),
     &                   SOFSIG_B(R,T)))
C
C --- sigma-theta as a function of temp (deg c) and salinity (psu)
C --- (polynomial fit that is cubic in T and quadratic in S)
C
      SIG(T,S)=(C1+S*(C3+S* C8)+
     &             T*(C2+S*(C5+S*C9)+T*(C4+S*C7+T*C6)))
#endif
C
      CALL XCSPMD
C
      CALL ZHOPEN(71, 'FORMATTED', 'OLD', 0)
      READ(71,*) !5-line header
      READ(71,*)
      READ(71,*)
      READ(71,*)
      READ(71,*)
      DO K= 1,999
        READ(71,*,END=100) !one line per level
      ENDDO
  100 CONTINUE
      KZ = K-1
      write(6,*) 'kz = ',kz
      CLOSE(71)
C
      ALLOCATE( TZ(KZ+1,IDM,JDM) )
      ALLOCATE( RZ(KZ+1,IDM,JDM) )
      ALLOCATE(    PMIX(IDM,JDM) )
      ALLOCATE(      PM(IDM,JDM) )
      ALLOCATE(      TM(IDM,JDM) )
      ALLOCATE(      SM(IDM,JDM) )
      ALLOCATE(      RM(IDM,JDM) )
      ALLOCATE(   DEPTH(IDM,JDM) )
      ALLOCATE(    PLAT(IDM,JDM) )
      ALLOCATE(    PKM1(IDM,JDM) )
      ALLOCATE(    RKM1(IDM,JDM) )
      ALLOCATE(    PKM2(IDM,JDM) )
      ALLOCATE(    WORK(IDM,JDM) )
      ALLOCATE(     MSK(IDM,JDM) )
C
C     BLKDAT INPUT.
C
      CALL ZHOPEN(99, 'FORMATTED', 'OLD', 0)
C
C --- 'ctitle' = climatology title
C
      WRITE(6,*)
      READ(99,'(A79)') CTITLE
      WRITE(6,'(A79)') CTITLE
C
C --- 'month'  = month (1 to 12)
C --- 'levtop' = top level of input clim. to use (optional, default 1)
C --- 'iversn' = hycom version number x10
C --- 'iexpt'  = experiment number x10
C --- 'yrflag' = days in year flag (0=360,1=366,2=366Jan1)
C --- 'mapflg' = map flag (0=mercator,1=rotated,2=uniform,3=beta-plane)
C --- 'idm   ' = longitudinal array size
C --- 'jdm   ' = latitudinal  array size
C --- 'kdm   ' = longitudinal array size
C --- 'nhybrd' = number of hybrid levels (0=all isopycnal)
C --- 'nsigma' = number of sigma  levels (nhybrd-nsigma z-levels)
C
      WRITE(6,*)
      CALL BLKINI(MONTH, 'month ')
      CALL BLKINI2(I,J,  'levtop','iversn')
      IF     (J.EQ.1) THEN
        LEVTOP = I
        CALL BLKINI(IVERSN,'iversn')
      ELSE
        IVERSN = I
        LEVTOP = 1
        write(6,'(a6," =",i6)') 'levtop',LEVTOP
      ENDIF
      WRITE(6,*)
      CALL BLKINI(IEXPT, 'iexpt ')
      CALL BLKINI(YRFLAG,'yrflag')
      WRITE(6,*)
      CALL BLKINI(I,     'idm   ')
      CALL BLKINI(J,     'jdm   ')
C
      IF     (I.NE.IDM) THEN
        WRITE(6,*)
        WRITE(6,*) 'ERROR - wrong IDM'
        WRITE(6,*)
        CALL ZHFLSH(6)
        STOP
      ELSEIF (J.NE.JDM) THEN
        WRITE(6,*)
        WRITE(6,*) 'ERROR - wrong JDM'
        WRITE(6,*)
        CALL ZHFLSH(6)
        STOP
      ENDIF
C
      CALL BLKINI(ITEST, 'itest ')
      CALL BLKINI(JTEST, 'jtest ')
      CALL BLKINI(KDM,   'kdm   ')
      CALL BLKINI(NHYBRD,'nhybrd')
      CALL BLKINI(NSIGMA,'nsigma')
C
      ISOPYC = NHYBRD.EQ.0
C
      IF     (IVERSN.LE.20) THEN
C
C ---   'dp00s'  = sigma   spacing minimum thickness (m)
C ---   'dp00'   = z-level spacing minimum thickness (m)
C ---   'dp00x'  = z-level spacing maximum thickness (m)
C ---   'dp00f'  = z-level spacing stretching factor (1.0=const.spacing)
C
        CALL BLKINR(DP00S, 'dp00s ','(a6," =",f10.4," m")')
        CALL BLKINR(DP00,  'dp00  ','(a6," =",f10.4," m")')
        CALL BLKINR(DP00X, 'dp00x ','(a6," =",f10.4," m")')
        CALL BLKINR(DP00F, 'dp00f ','(a6," =",f10.4," ")')
        DS00   = MAX(DP00S,0.01)
        DS00X  = DS00
        DS00F  = 1.0
        IF     (ISOPYC) THEN
          ISOTOP = -1.0  !turn off isotop
        ELSE
          ISOTOP = 0.01  !first layer always fixed
        ENDIF
      ELSE
C
C ---   'isotop' = shallowest depth for isopycnal layers (m), optional
C ---   'dp00'   = deep    z-level spacing minimum thickness (m)
C ---   'dp00x'  = deep    z-level spacing maximum thickness (m)
C ---   'dp00f'  = deep    z-level spacing stretching factor (1.0=const.z)
C ---   'ds00'   = shallow z-level spacing minimum thickness (m)
C ---   'ds00x'  = shallow z-level spacing maximum thickness (m)
C ---   'ds00f'  = shallow z-level spacing stretching factor (1.0=const.z)
C
C ---   the above describe a system that is isopycnal or:
C ---       z in    deep water, based on dp00,dp00x,dp00f
C ---       z in shallow water, based on ds00,ds00x,ds00f and nsigma
C ---       sigma between them, based on ds00,ds00x,ds00f and nsigma
C
C ---   for z-only set nsigma=0 (and ds00,ds00x,ds00f=dp00,dp00x,dp00f)
C ---   for sigma-z (shallow-deep) use a very small ds00
C ---    (pure sigma-z also has ds00f=dp00f and ds00x=dp00x*ds00/dp00)
C ---   for z-sigma (shallow-deep) use a very large dp00 (not recommended)
C ---   for sigma-only set nsigma=kdm, dp00 large, and ds00 small
C
        CALL BLKINR2(BLK,J,'isotop','(a6," =",f10.4," m")',
     &                     'dp00  ','(a6," =",f10.4," m")')
        IF     (J.EQ.1) THEN
          ISOTOP = BLK
          CALL BLKINR(DP00,  'dp00  ','(a6," =",f10.4," m")')
        ELSE
          DP00   = BLK
          IF     (ISOPYC) THEN
            ISOTOP = -1.0  !turn off isotop
          ELSE
            ISOTOP = 0.01  !first layer always fixed
          ENDIF
          write(6,'(a6," =",f10.4," m")') 'isotop',ISOTOP
        ENDIF
        CALL BLKINR(DP00X, 'dp00x ','(a6," =",f10.4," m")')
        CALL BLKINR(DP00F, 'dp00f ','(a6," =",f10.4," ")')
        CALL BLKINR(DS00,  'ds00  ','(a6," =",f10.4," m")')
        CALL BLKINR(DS00X, 'ds00x ','(a6," =",f10.4," m")')
        CALL BLKINR(DS00F, 'ds00f ','(a6," =",f10.4," ")')
      ENDIF
C
      IF (NSIGMA.LE.1) THEN
        NSIGMA=1
        DS00  =DP00
        DS00X =DP00X
        DS00F =DP00F
      ENDIF
C
C --- 'thflag' = reference pressure flag (0=Sigma-0, 2=Sigma-2)
      WRITE(6,*)
      CALL BLKINI(THFLAG,'thflag')
      IF     (THFLAG.EQ.0) THEN
        IF     (ABS(SOFSIG(AZERO,AZERO)-
     +              0.136471/0.807004   ).GT.0.001) THEN
          WRITE(6,*)
          WRITE(6,*) 'ERROR - thflag not consistent with sofsig()'
          WRITE(6,*) SOFSIG(AZERO,AZERO)
          WRITE(6,*)
          CALL ZHFLSH(6)
          STOP
        ENDIF
      ELSEIF (THFLAG.EQ.2) THEN
        IF     (ABS(SOFSIG(AZERO,AZERO)+
     +              9.77093 /0.789879   ).GT.0.001) THEN
          WRITE(6,*)
          WRITE(6,*) 'ERROR - thflag not consistent with sofsig()'
          WRITE(6,*) SOFSIG(AZERO,AZERO)
          WRITE(6,*)
          CALL ZHFLSH(6)
          STOP
        ENDIF
      ELSEIF (THFLAG.EQ.4) THEN
        IF     (ABS(SOFSIG(AZERO,AZERO)+
     +              19.2362 /0.773552  ).GT.0.001) THEN
          WRITE(6,*)
          WRITE(6,*) 'ERROR - thflag not consistent with sofsig()'
          WRITE(6,*) SOFSIG(AZERO,AZERO)
          WRITE(6,*)
          CALL ZHFLSH(6)
          STOP
        ENDIF
      ELSE
        WRITE(6,*)
        WRITE(6,*) 'ERROR - thflag must be 0 or 2 or 4'
        WRITE(6,*)
        CALL ZHFLSH(6)
        STOP
      ENDIF
      IF     (THFLAG.EQ.0) THEN
        SIGFMT = '(a6," =",f10.4," sigma-0")'
      ELSEIF (THFLAG.EQ.2) THEN
        SIGFMT = '(a6," =",f10.4," sigma-2")'
      ELSEIF (THFLAG.EQ.4) THEN
        SIGFMT = '(a6," =",f10.4," sigma-4")'
      ENDIF
C
C --- 'thbase' = reference density (sigma units)
      WRITE(6,*)
      CALL BLKINR(THBASE,'thbase',SIGFMT)
#if defined(sigmav)
C
C --- 'vsigma' = spacially varying isopycnal layer target densities (0=F,1=T)
      WRITE(6,*)
      CALL BLKINL(VSIGMA,'vsigma')
#else
      VSIGMA = .FALSE.
#endif
C --- 'sigma ' = isopycnal layer target densities (sigma units)
      DO K=1,KDM
        CALL BLKINR(SIGMA(K),'sigma ',SIGFMT)
        IF     (K.GT.1) THEN
          IF     (SIGMA(K).LT.SIGMA(K-1)) THEN
            WRITE(6,*)
            WRITE(6,*) 'ERROR - sigma(k) must be > sigma(k-1)'
            WRITE(6,*)
            CALL ZHFLSH(6)
            STOP
          ENDIF
        ENDIF
      ENDDO
C
C --- 'thkmin' = minimum mixed-layer thickness (m)
      WRITE(6,*)
      CALL BLKINR(THKMIN,'thkmin','(a6," =",f10.4," m")')
      WRITE(6,*)
      CLOSE(UNIT=99)
C
C     MINIMUM (DEEP) LAYER THICKNESSES.
C
      IF     (ISOPYC) THEN
        DP0K(1)=THKMIN
      ELSE
        DP0K(1)=DP00
      ENDIF
      DPMS=DP0K(1)
      DPCK(0)=0.0
      DPCK(1)=DPMS
      IF     (THFLAG.EQ.0) THEN
        WRITE(6,6000) 1,DP0K(1),DPMS,SIGMA(1),' Sig-th'
      ELSEIF (THFLAG.EQ.0) THEN
        WRITE(6,6000) 1,DP0K(1),DPMS,SIGMA(1),' Sigma2'
      ELSEIF (THFLAG.EQ.0) THEN
        WRITE(6,6000) 1,DP0K(1),DPMS,SIGMA(1),' Sigma4'
      ENDIF
 6000 FORMAT(       'k =',I3,
     +       '   thkns =',F6.1,' m',
     +       '   depth =',F7.1,' m',
     +       '   density =',F7.3,A)
C
      DP0KF=ONE
      DO K=2,KDM
        DP0KF=DP0KF*DP00F
        IF     (K.LE.NHYBRD) THEN
          DP0K(K)=MIN(DP00*DP0KF,DP00X)
        ELSE
          DP0K(K)=ZERO
        ENDIF
        DPMS=DPMS+DP0K(K)
        DPCK(K)=DPMS
        IF     (THFLAG.EQ.0) THEN
          WRITE(6,6000) K,DP0K(K),DPMS,SIGMA(K),' Sig-th'
        ELSEIF (THFLAG.EQ.0) THEN
          WRITE(6,6000) K,DP0K(K),DPMS,SIGMA(K),' Sigma2'
        ENDIF
      ENDDO
      WRITE(6,*)
      CALL ZHFLSH(6)
C
C     MINIMUM (SHALLOW) LAYER THICKNESSES.
C
      IF     (ISOPYC) THEN
        DS0K(1)=THKMIN
      ELSE
        DS0K(1)=DS00
      ENDIF
      DSMS=DS0K(1)
      DSCK(0)=0.0
      DSCK(1)=DSMS
      DS0KF=ONE
      DO K=2,NSIGMA
        DS0KF=DS0KF*DS00F
        DS0K(K)=MIN(DS00*DS0KF,DS00X)
        DSMS=DSMS+DS0K(K)
        DSCK(K)=DSMS
      ENDDO
      SIGCUM(0)=0.0
      DO K=1,NSIGMA
        SIGSCL(K)=DS0K(K)/DSMS
        SIGCUM(K)=DSCK(K)/DSMS
      ENDDO
      DO K= NSIGMA+1,KDM
        DS0K(  K)=DP0K(K)
        SIGSCL(K)=0.0
        SIGCUM(K)=1.0
      ENDDO
C
C     TOPOGRAPHY INPUT.
C
      CALL ZAIOST
C
      CALL ZHOPEN(51, 'FORMATTED', 'OLD', 0)
      READ (51,'(A79)') PREAMBL
      READ (51,'(A)')   CLINE
      CLOSE(UNIT=51)
      WRITE(6,'(/(1X,A79))') PREAMBL,CLINE
C
      I = INDEX(CLINE,'=')
      READ (CLINE(I+1:),*)   HMINB,HMAXB
C
      CALL ZAIOPN('OLD', 51)
      CALL ZAIORD(DEPTH,MSK,.FALSE., HMINA,HMAXA, 51)
      CALL ZAIOCL(51)
C
      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 topography files not consistent:',
     &    '.a,.b min = ',HMINA,HMINB,HMINA-HMINB,
     &    '.a,.b max = ',HMAXA,HMAXB,HMAXA-HMAXB
        CALL ZHFLSH(6)
        STOP
      ENDIF
C
C     INITIALIZE LAND MASK.
C
      DO J= 1,JDM
        DO I= 1,IDM
          IF     (DEPTH(I,J).LT.2.0**99) THEN
            MSK(  I,J) = 1
          ELSE
            MSK(  I,J) = 0
            DEPTH(I,J) = ZERO
          ENDIF
        ENDDO
      ENDDO
C
C     TARGET DENSITIES.
C
      ALLOCATE( SIG3D(IDM,JDM,KDM) )
C
      IF     (.NOT.VSIGMA) THEN
        DO K= 1,KDM
          SIG3D(:,:,K) = SIGMA(K)
        ENDDO
      ELSE
        CALL ZAIOPN('OLD', 52)
        DO K= 1,KDM
          CALL ZAIORD(SIG3D(1,1,K),MSK,.FALSE., HMINA,HMAXA, 52)
          IF     (HMINA.GT.SIGMA(K)+0.005 .OR.
     &            HMAXA.LT.SIGMA(K)-0.005     ) THEN
            WRITE(6,'(/ a,i3,a /)')
     &        'ERROR - VARIABLE TARGET DENSITY FOR LAYER',
     &        K,' IS NOT CONSISTENT WITH SIGMA(K)'
            WRITE(6,*) 'SIGMA(K)    = ',SIGMA(K)
            WRITE(6,*) 'HMINA,HMAXA = ',HMINA,HMAXA
            STOP
          ENDIF
        ENDDO
        CALL ZAIOCL(52)
      ENDIF
C
C     CHECK ITEST,JTEST.
C
      IF     (MIN(ITEST,JTEST).GT.0) THEN
        IF     (ITEST.GT.IDM) THEN
          WRITE(6,'(/ a /)') 'error - itest > idm'
          CALL ZHFLSH(6)
          STOP
        ELSEIF (JTEST.GT.JDM) THEN
          WRITE(6,'(/ a /)') 'error - jtest > jdm'
          CALL ZHFLSH(6)
          STOP
        ELSEIF(MSK(ITEST,JTEST).EQ.0) THEN
          WRITE(6,'(/ a /)') 'error - itest,jtest is a land point'
          CALL ZHFLSH(6)
          STOP
        ENDIF
      ENDIF
C
C     LATITUDE GRID INPUT.
C
      CALL ZHOPNC(31, 'regional.grid.b', 'FORMATTED', 'OLD', 0)
      CALL ZAIOPF('regional.grid.a', 'OLD', 31)
C
      READ(31,*) ! skip idm
      READ(31,*) ! skip jdm
      READ(31,*) ! skip mapflg
      READ(31,*) ! skip plon
      CALL ZAIOSK(31)
      READ(31,'(A)') CLINE
      I = INDEX(CLINE,'=')
      READ (CLINE(I+1:),*)   HMINB,HMAXB
      CALL ZAIORD(PLAT,MSK,.FALSE., HMINA,HMAXA, 31)
      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
C     Z-LEVEL CLIMATOLOGY INPUT.
C
      CALL ZAIOPN('OLD', 71)
      CALL ZHOPEN(71, 'FORMATTED', 'OLD', 0)
      READ (71,'(A79)') PREAMBL
      WRITE(6, '(/(1X,A79))') PREAMBL
      DO K= 1,KZ
        READ (71,'(A)') CLINE
        I = INDEX(CLINE,'=')
        READ (CLINE(I+1:),*) ZLEV(K),HMINB,HMAXB
        IF     (HMINB.EQ.HMAXB) THEN  ! constant field
          DO J= 1,JDM
            DO I= 1,IDM
              RZ(K,I,J) = HMINB
            ENDDO
          ENDDO
          CALL ZAIOSK(71)
        ELSE
          CALL ZAIORD(WORK,MSK,.FALSE., HMINA,HMAXA, 71)
          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 density files not consistent:',
     &       '.a,.b min = ',HMINA,HMINB,HMINA-HMINB,
     &       '.a,.b max = ',HMAXA,HMAXB,HMAXA-HMAXB
            CALL ZHFLSH(6)
            STOP
          ENDIF
          DO J= 1,JDM
            DO I= 1,IDM
              RZ(K,I,J) = WORK(I,J)
            ENDDO
          ENDDO
        ENDIF
*       if     (min(itest,jtest).gt.0) then
*         write(6,'(a,2i5,i3,a,f7.3)')
*    +     'i,j,k =',itest,jtest,k,
*    +     '   r =',rz(k,itest,jtest)
*       endif
C
        IF     (K.GT.1) THEN
          IF     (ZLEV(K).GT.399.0 .AND. ZLEV(K-1).LT.399) THEN
            K400 = K
          ENDIF
        ENDIF
      ENDDO
      CLOSE(UNIT=71)
      CALL ZAIOCL(71)
C
C     CONVERT Z-LEVELS TO FINITE VOLUMES.
C
      ZI(0) = ZERO
      DO K= 1,KZ-1
        ZI(K) = 0.5*(ZLEV(K) + ZLEV(K+1))
      ENDDO
      ZI(KZ) = 9999.0
C
      CALL ZAIOPN('OLD', 72)
      CALL ZHOPEN(72, 'FORMATTED', 'OLD', 0)
      READ (72,'(A79)') PREAMBL
      WRITE(6, '(/(1X,A79))') PREAMBL
      DO K= 1,KZ
        READ (72,'(A)') CLINE
        I = INDEX(CLINE,'=')
        READ (CLINE(I+1:),*) ZLEV(K),HMINB,HMAXB
        IF     (HMINB.EQ.HMAXB) THEN  ! constant field
          DO J= 1,JDM
            DO I= 1,IDM
              WORK(I,J) = HMINB
            ENDDO
          ENDDO
          CALL ZAIOSK(72)
        ELSE
          CALL ZAIORD(WORK,MSK,.FALSE., HMINA,HMAXA, 72)
          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 temperature files not consistent:',
     &        '.a,.b min = ',HMINA,HMINB,HMINA-HMINB,
     &        '.a,.b max = ',HMAXA,HMAXB,HMAXA-HMAXB
            CALL ZHFLSH(6)
            STOP
          ENDIF
        ENDIF
C
        IF     (.NOT.ISOPYC) THEN
          DO J= 1,JDM
            DO I= 1,IDM
              TZ(K,I,J) = WORK(I,J)
            ENDDO
          ENDDO
        ELSE
C
C         LIMIT MAXIMUM DENSITY TO SIGMA(KDM)
C
          DO J= 1,JDM
            DO I= 1,IDM
              TZ(K,I,J) = WORK(I,J)
              IF     (RZ(K,I,J).GT.SIG3D(I,J,KDM)) THEN
                IF     (RZ(MAX(K-1,1),I,J).EQ.SIG3D(I,J,KDM)) THEN
                  RZ(K,I,J) = SIG3D(I,J,KDM)
                  TZ(K,I,J) = TZ(K-1,I,J)
                ELSE
                  SZK       = SOFSIG( R8(RZ(K,I,J)),  R8(TZ(K,I,J)) )
                  RZ(K,I,J) =            SIG3D(I,J,KDM)
                  TZ(K,I,J) = TOFSIG( R8(SIG3D(I,J,KDM)), R8(SZK) )
                ENDIF
              ENDIF
            ENDDO
          ENDDO
        ENDIF
*       if     (min(itest,jtest).gt.0) then
*         write(6,'(a,2i5,i3,a,f7.3)')
*    +     'i,j,k =',itest,jtest,k,
*    +     '   r =',rz(k,itest,jtest)
*         write(6,'(a,2i5,i3,a,f7.3)')
*    +     'i,j,k =',itest,jtest,k,
*    +     '   t =',tz(k,itest,jtest)
*         write(6,'(a,2i5,i3,a,f7.3)')
*    +     'i,j,k =',itest,jtest,k,
*    +     '   w =',work(itest,jtest)
*       endif
      ENDDO
      CLOSE (UNIT=72)
      CALL ZAIOCL(72)
C
C     MAKE LEVELS 1:LEVTOP IDENTICAL
C
      IF     (LEVTOP.GT.1) THEN
        DO J= 1,JDM
          DO I= 1,IDM
            DO K= 1,LEVTOP-1
              RZ(K,I,J) = RZ(LEVTOP,I,J)
              TZ(K,I,J) = TZ(LEVTOP,I,J)
            ENDDO
          ENDDO
        ENDDO
      ENDIF
C
      WRITE(6,'(/A,10F7.1/(7X,10F7.1))') 'ZI = ',ZI
C
C     DIAGNOSTIC PRINTOUT.
C
      IF     (MIN(ITEST,JTEST).GT.0) THEN
        WRITE(6,*)
        DO KZL= 0,KZ
          ZIL(KZL) = ZI(KZL)
          IF     (ZI(KZL).GT.DEPTH(ITEST,JTEST)) THEN
            ZIL(KZL) = DEPTH(ITEST,JTEST)
            EXIT
          ENDIF
        ENDDO
        DO K= 1,MAX(KZL-1,2)
          SZK = SOFSIG( R8(RZ(K,ITEST,JTEST)), R8(TZ(K,ITEST,JTEST)) )
          WRITE(6,'(A,2I5,I3,A,F8.2,A,3F7.3)')
     +     'I,J,K =',ITEST,JTEST,K,
     +     '  ZIL =',ZIL(K),
     +     '   R,T,S =',RZ(K,ITEST,JTEST),
     +                  TZ(K,ITEST,JTEST),SZK
        ENDDO
        IF     (KZL.GE.3) THEN
          ZZ  = 0.5*(ZIL(KZL-1)+ZIL(KZL))
          Q   = (ZZ - ZI(KZL-2))/(ZI(KZL-1) - ZI(KZL-2))
          RZK = (1.0-Q)*RZ(KZL-2,ITEST,JTEST) +
     +               Q *RZ(KZL-1,ITEST,JTEST)
          TZK = (1.0-Q)*TZ(KZL-2,ITEST,JTEST) +
     +               Q *TZ(KZL-1,ITEST,JTEST)
          SZK = SOFSIG( R8(RZK), R8(TZK) )
          WRITE(6,'(A,2I5,I3,A,F8.2,A,3F7.3 /)')
     +     'I,J,K =',ITEST,JTEST,KZL,
     +     '  ZIL =',ZIL(K),
     +     '   R,T,S =',RZK,TZK,SZK
          WRITE(6,'(A,5F10.3 / A,5F10.3 /)')
     +      'BATHY - ZZ,Q,RA,RB,RZ =',
     +      ZZ,Q,RZ(KZL-2,ITEST,JTEST),RZ(KZL-1,ITEST,JTEST),RZK,
     +      'BATHY - ZZ,Q,TA,TB,TZ =',
     +      ZZ,Q,TZ(KZL-2,ITEST,JTEST),TZ(KZL-1,ITEST,JTEST),TZK
          CALL ZHFLSH(6)
        ENDIF
      ENDIF
C
C     INITIALIZE CLIMATOLOGY OUTPUT.
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
C
      IF     (NSIGMA.EQ.1) THEN
        WRITE(PREAMBL(2),4000) IEXPT/10,MOD(IEXPT,10),
     +                         NHYBRD,0,
     +                         DP00,DP00,DP00X,DP00F
      ELSE
        WRITE(PREAMBL(2),4000) IEXPT/10,MOD(IEXPT,10),
     +                         NHYBRD,NSIGMA,
     +                         DS00,DP00,DP00X,DP00F
      ENDIF
C
      WRITE(PREAMBL(3),'(A,I1,A,I1)')
     +        'Layered averages w.r.t. Sigma-',THFLAG,
     +        ',  levtop=',LEVTOP
C
      WRITE(PREAMBL(5),'(A,2I5,I3,F9.3,F9.2,2F6.3)')
     +        'i/jdm =',
     +       IDM,JDM
C
      PREAMBL(4) = 'Potential Temperature'
      WRITE(10,4101) PREAMBL
C
      PREAMBL(4) = 'Salinity'
      WRITE(11,4101) PREAMBL
C
      PREAMBL(4) = 'Interface Depths'
      WRITE(12,4101) PREAMBL
C
      WRITE(6,*)
      PREAMBL(4) = 'Potential Temperature, or ' //
     +             'Salinity, or ' //
     +             'Interface Depths'
      WRITE(6, 4101) PREAMBL
      WRITE(6,*)
C
C     INITIALIZE DUMMY ARCHIVE OUTPUT.
C
      CALL ZAIOPN('NEW', 21)
      CALL ZHOPEN(21, 'FORMATTED', 'NEW', 0)
      YRFLAG = MAX(0,MIN(2,YRFLAG))
      WRITE(21,4200) MONTH,PREAMBL(1),PREAMBL(2),
     &               IVERSN,IEXPT,YRFLAG,IDM,JDM

      IF     (YRFLAG.EQ.0) THEN  ! 360 days, starting Jan 16
        TIME = (MONTH-1)*30.0
      ELSEIF (YRFLAG.EQ.1) THEN  ! 366 days, starting Jan 16
        TIME = (MONTH-1)*30.5
      ELSEIF (YRFLAG.EQ.2) THEN  ! 366 days, starting Jan 1
        TIME = (MONTH-1)*30.5 + 15.0
      ENDIF
C
C     FIRST INTERFACE IS AT SURFACE.
C     MIXED LAYER AVERAGE TEMPERATURE IS NOMINALLY SST-0.25 DEG.
C     THIS IS EQUIVALENT TO A 0.5 DEG CHANGE ACROSS THE MIXED
C     LAYER, WHICH IS LARGER THAN TYPICALLY USED FOR ACTUAL
C     PROFILES BECAUSE THIS IS A SMOOTH CLIMATOLOGICAL PROFILE.
C
      DO J= 1,JDM
        DO I= 1,IDM
          IF     (DEPTH(I,J).GT.ZERO) THEN
            PKM2(I,J) = ZERO
            PKM1(I,J) = ZERO
            TM(  I,J) = TZ(1,I,J)-0.25
            SM(  I,J) = SOFSIG( R8(RZ(1,I,J)),
     +                          R8(TZ(1,I,J)) )
            SIGMAA    = SIG(    R8(TZ(1,I,J)-0.5),
     +                          R8(SM(  I,J)    ) )
            PMIX(I,J) = 400.0
            DO K= 2,K400
              IF     (RZ(K,I,J).GE.SIGMAA) THEN
                SIGMAB = RZ(K,I,J)-RZ(K-1,I,J)
                IF     (SIGMAB.GT.1.E-4) THEN
                  Q = (RZ(K,I,J)-SIGMAA)/SIGMAB
                  PMIX(I,J) = ZLEV(K-1) + (1.0-Q)*(ZLEV(K)-ZLEV(K-1))
                ELSE
                  PMIX(I,J) = ZLEV(K)
                ENDIF
                EXIT
              ENDIF
            ENDDO
            PMIX(I,J) = MAX( THKMIN, MIN( PMIX(I,J), DEPTH(I,J) ) )
          ELSE
            PKM2(I,J) = ZERO
            PKM1(I,J) = ZERO
            PMIX(I,J) = THKMIN
            TM(  I,J) = ZERO
            SM(  I,J) = ZERO
            RM(  I,J) = ZERO
C
            PM(  I,J) = ZERO
C
            DEPTH(I,J) = ZERO  ! should already be zero
          ENDIF
        ENDDO
      ENDDO
C
C     PUT "SURFACE" FIELDS IN THE THE DUMMY ARCHIVE.
C
      DO J= 1,JDM
        DO I= 1,IDM
          WORK(I,J) = ZERO
        ENDDO
      ENDDO
      CALL ZAIOWR(WORK,MSK,.TRUE.,  XMIN,XMAX, 21, .FALSE.)
      WRITE(21,4201) 'montg1  ',MONTH,TIME,0,ZERO,XMIN,XMAX
      CALL ZAIOWR(WORK,MSK,.TRUE.,  XMIN,XMAX, 21, .FALSE.)
      WRITE(21,4201) 'srfhgt  ',MONTH,TIME,0,ZERO,XMIN,XMAX
      CALL ZAIOWR(WORK,MSK,.TRUE.,  XMIN,XMAX, 21, .FALSE.)
      WRITE(21,4201) 'surflx  ',MONTH,TIME,0,ZERO,XMIN,XMAX
      CALL ZAIOWR(WORK,MSK,.TRUE.,  XMIN,XMAX, 21, .FALSE.)
      WRITE(21,4201) 'salflx  ',MONTH,TIME,0,ZERO,XMIN,XMAX
C
      DO J= 1,JDM
        DO I= 1,IDM
          IF     (DEPTH(I,J).GT.ZERO) THEN
            WORK(I,J) = PMIX(I,J)*9806.0
          ENDIF
        ENDDO
      ENDDO
      CALL ZAIOWR(WORK,MSK,.TRUE.,  XMIN,XMAX, 21, .FALSE.)
      WRITE(21,4201) 'bl_dpth ',MONTH,TIME,0,ZERO,XMIN,XMAX
      CALL ZAIOWR(WORK,MSK,.TRUE.,  XMIN,XMAX, 21, .FALSE.)
      WRITE(21,4201) 'mix_dpth',MONTH,TIME,0,ZERO,XMIN,XMAX
      DO J= 1,JDM
        DO I= 1,IDM
          IF     (DEPTH(I,J).GT.ZERO) THEN
            TM(  I,J) = TZ(1,I,J)-0.25
            SM(  I,J) = SOFSIG( R8(RZ(1,I,J)),
     +                          R8(TZ(1,I,J)) )
            WORK(I,J) = SIG(    R8(TM(  I,J)),
     +                          R8(SM(  I,J)) ) - THBASE
          ELSE
            TM(  I,J) = ZERO
            SM(  I,J) = ZERO
            WORK(I,J) = ZERO
          ENDIF
        ENDDO
      ENDDO
      CALL ZAIOWR(TM,  MSK,.TRUE.,  XMIN,XMAX, 21, .FALSE.)
      WRITE(21,4201) 'tmix    ',MONTH,TIME,0,  ZERO,XMIN,XMAX
      CALL ZAIOWR(SM,  MSK,.TRUE.,  XMIN,XMAX, 21, .FALSE.)
      WRITE(21,4201) 'smix    ',MONTH,TIME,0,  ZERO,XMIN,XMAX
      CALL ZAIOWR(WORK,MSK,.TRUE.,  XMIN,XMAX, 21, .FALSE.)
      WRITE(21,4201) 'thmix   ',MONTH,TIME,0,THBASE,XMIN,XMAX
C
      DO J= 1,JDM
        DO I= 1,IDM
          WORK(I,J) = ZERO
        ENDDO
      ENDDO
      CALL ZAIOWR(WORK,MSK,.FALSE., XMIN,XMAX, 21, .FALSE.)
      WRITE(21,4201) 'umix    ',MONTH,TIME,0,ZERO,XMIN,XMAX
      CALL ZAIOWR(WORK,MSK,.FALSE., XMIN,XMAX, 21, .FALSE.)
      WRITE(21,4201) 'vmix    ',MONTH,TIME,0,ZERO,XMIN,XMAX
      CALL ZAIOWR(WORK,MSK,.FALSE., XMIN,XMAX, 21, .FALSE.)
      WRITE(21,4201) 'u_btrop ',MONTH,TIME,0,ZERO,XMIN,XMAX
      CALL ZAIOWR(WORK,MSK,.FALSE., XMIN,XMAX, 21, .FALSE.)
      WRITE(21,4201) 'v_btrop ',MONTH,TIME,0,ZERO,XMIN,XMAX
C
C     PROCESS ALL THE LAYERS.
C
      DO J= 1,JDM
        DO I= 1,IDM
          IF     (DEPTH(I,J).GT.ZERO) THEN
            RKM1(I,J) = ZERO
            PKM1(I,J) = ZERO
            PKM2(I,J) = ZERO
          ENDIF
        ENDDO
      ENDDO
      DO 810 K= 1,KDM
        DO J= 1,JDM
          DO I= 1,IDM
            IF     (DEPTH(I,J).GT.ZERO) THEN
C
C             ALLOW FOR BATHYMETRY.
C
              DO KZL= 0,KZ
                ZIL(KZL) = ZI(KZL)
                IF     (ZI(KZL).GT.DEPTH(I,J)) THEN
                  ZIL(KZL) = DEPTH(I,J)
                  EXIT
                ENDIF
              ENDDO
              IF     (K.EQ.1 .AND. KZL.GE.3) THEN
C
C               MODIFY LOWEST FINITE VOLUME USING EXTRAPOLATION
C               FROM ABOVE, BECAUSE LOWEST VALUES LIKELY TO BE
C               BELOW THE BOTTOM AND THEREFORE INACCURATE.
C
                ZZ = 0.5*(ZIL(KZL-1)+ZIL(KZL))
                Q  = (ZZ - ZI(KZL-2))/(ZI(KZL-1) - ZI(KZL-2))
                RZ(KZL,I,J) = (1.0-Q)*RZ(KZL-2,I,J) +
     +                             Q *RZ(KZL-1,I,J)
                TZ(KZL,I,J) = (1.0-Q)*TZ(KZL-2,I,J) +
     +                             Q *TZ(KZL-1,I,J)
              ENDIF
C
              IF     (K.GT.1) THEN
                RM(I,J) =     SIG3D(I,J,K)
              ELSE
                RM(I,J) = MIN(SIGMA(1),RZ(1,I,J))
              ENDIF
C
C             FIND RM AND PM.
C
              IF     (NSIGMA.GT.1 .AND.
     +                DEPTH(I,J).LT.DPCK(NSIGMA)) THEN
                DMIN = SIGCUM(K-1)*DEPTH(I,J)
              ELSE
                DMIN = DPCK(K-1)
              ENDIF
*             if (i.eq.itest .and. j.eq.jtest) then
*               WRITE(6,'(A,I3,F10.3)')
*    +            'ISOTOP - K,DMIN =',K,DMIN
*             endif
              IF     (ISOPYC .AND. K.EQ.1) THEN
C
C               UPPER LAYER IS THE MIXED LAYER.
C
                SIGMAA = ZERO
                DO KK= 1,KZL
                  ZTOP   =              ZIL(KK-1)
                  ZBOT   = MAX(ZTOP,MIN(ZIL(KK  ),PMIX(I,J)))
                  SIGMAA = SIGMAA + RZ(KK,I,J)*(ZBOT-ZTOP)
                ENDDO
                RM(I,J) = SIGMAA/PMIX(I,J)
                PM(I,J) = PMIX(I,J)
              ELSEIF (K.EQ.KDM) THEN
C
C               LOWEST LAYER.
C
                PM(I,J) = DEPTH(I,J)
              ELSEIF (DMIN.LE.ISOTOP) THEN
C
C               FIXED GRID NEAR THE SURFACE.
C
                IF     (NSIGMA.GT.1 .AND.
     +                  DEPTH(I,J).LT.DPCK(NSIGMA)) THEN
                  PM(I,J) = MAX( PKM1(I,J), SIGCUM(K)*DEPTH(I,J) )
                ELSE
                  PM(I,J) = MAX( PKM1(I,J), DPCK(K) )
                ENDIF
*               if (i.eq.itest .and. j.eq.jtest) then
*                 WRITE(6,'(A,I3,F10.3)')
*    +              'FIXED - K,PM =',K,PM(I,J)
*               endif
              ELSEIF (RM(I,J).LT.RKM1(I,J)) THEN
C
C               ZERO THICKNESS LAYER REQUIRED FOR STABILITY.
C
                PM(I,J) = PKM1(I,J)
              ELSE  !K.LT.KDM
C
C               REMAP TO EXACTLY ISOPYCNAL LAYER
C
                PINTEG  = ZERO
                SIGMAA = ZERO
                SIGMAB = ZERO
                  THKA = ZERO
                  THKB = ZERO
                DO KK= 1,KZL
                  ZTOP   =          MAX(ZIL(KK-1),PKM1(I,J))
                  ZBOT   = MAX(ZTOP,MIN(ZIL(KK  ),PM(  I,J)))
                  IF     (ZBOT.NE.ZTOP) THEN
                    SIGMAA = SIGMAA + (ZBOT-ZTOP)*RZ(KK,I,J)
                    THKA   = THKA   + (ZBOT-ZTOP)
*                   if (i.eq.itest .and. j.eq.jtest) then
*                     WRITE(6,'(A,I3,3F10.3)')
*    +                  'REMAP - K,SA,SB,THK =',
*    +                  KK,SIGMAA/THKA,SIGMAB/MAX(0.001,THKB),THKA
*                   endif
                    IF     (SIGMAA/THKA.GT.SIG3D(I,J,K))
C
C                     LAYER INTERFACE IS SOMEWHERE INSIDE LEVEL KK.
C
                      THK    = THKB * (SIG3D(I,J,K)-SIGMAB/THKB)/
     +                                (RZ(KK,I,J)-SIG3D(I,J,K))
                      SIGMAB = SIGMAB + THK*RZ(KK,I,J)
                        THKB =   THKB + THK
                      EXIT
                    ELSEIF (SIGMAA/THKA.EQ.SIG3D(I,J,K))
C
C                     LAYER INTERFACE IS AT BOTTOM OF LEVEL KK
C
                      SIGMAB = SIGMAA
                      THKB   = THKA
                      EXIT
                    ENDIF
                    IF     (THKA.NE.THKB) THEN
                      THKB   = THKA
                      SIGMAB = SIGMAA  !layer so far's volume density
                    ENDIF
                  ENDIF
                ENDDO !kk
                IF     (THKB.GT.ZERO) THEN
C
C                 RM, NORMALLY SIG3D.K
C
                  PM(I,J) = MIN(DEPTH(I,J), PKM1(I,J)+THKB )
                  RM(I,J) = MAX( RKM1(I,J), SIGMAB/THKB )
                ELSEIF (RM(I,J).LT.RZ(KZL,I,J)) THEN
C
C                 AT THE TOP.
C
                  PM(I,J) = ZERO
                ELSE
C
C                 AT THE BOTTOM.
C
                  PM(I,J) = DEPTH(I,J)
                ENDIF
              ENDIF
C
C             MODIFY RM AND PM (IF NECESSARY).
C
              IF     (DMIN.LE.ISOTOP) THEN
                DPMS = PM(I,J)
              ELSE
                DPMS = PKM1(I,J) + MIN(DP0K(K),
     +                                 MAX(DS0K(K),
     +                                     SIGSCL(K)*DEPTH(I,J)))
              ENDIF
              IF     (.NOT. ISOPYC .AND.
     +                PM(I,J).EQ.DEPTH(I,J)) THEN
                IF     (PM(I,J).GT.PKM1(I,J)) THEN
C
C                 MAKE LOWEST NON-ZERO LAYER EXACTLY ISOPYNCAL,
C                 UNLESS THE LAYER ABOVE IS NOT ISOPYCNAL.
C
                  IF     (K.GT.1) THEN
                    IF     (RKM1(I,J).NE.SIG3D(I,J,K-1)) THEN
                      THK    = PKM1(I,J)-PKM2(I,J)
                      THIKMN = MIN(DP0K(K-1),
     +                             MAX(DS0K(K-1),
     +                                 SIGSCL(K-1)*DEPTH(I,J)))
                    ELSE  ! LAYER ABOVE IS DEFINATELY ISOPYCNAL
                      THK    = 1.0
                      THIKMN = 0.0
                    ENDIF
                  ELSE
                    THK    = 1.0
                    THIKMN = 1.0
                  ENDIF
                  IF     (        THIKMN .EQ.ZERO .OR.
     +                    ABS(THK-THIKMN).GT.0.01     ) THEN
                    RM(I,J) = MAX(RKM1(I,J),SIG3D(I,J,K))
                  ELSE
C
C                   LOWEST NON-ZERO LAYER IS NOT-ISOPYCNAL.
C
                    SIGMAA = ZERO
                    DO KK= 1,KZL
                      ZTOP   =          MAX(ZIL(KK-1),PKM1(I,J))
                      ZBOT   = MAX(ZTOP,MIN(ZIL(KK  ),PM(  I,J)))
                      SIGMAA = SIGMAA + RZ(KK,I,J)*(ZBOT-ZTOP)
                    ENDDO
                    RM(I,J) = MAX(RKM1(I,J),SIGMAA/(PM(I,J)-PKM1(I,J)))
                  ENDIF
                ELSE
                  RM(I,J) =     RKM1(I,J)
                ENDIF
              ELSEIF (.NOT. ISOPYC .AND. 
     +                PM(I,J).LE.MIN(DEPTH(I,J),DPMS)) THEN
C
C               HYBRID (MINIMUM THICKNESS) LAYER.
C
                PM(I,J) = MIN(DPMS,DEPTH(I,J))
*               if (i.eq.itest .and. j.eq.jtest) then
*                 WRITE(6,'(A,I3,2F10.3)')
*    +              'HYBRID - K,PKM1,PM =',K,PKM1(I,J),PM(I,J)
*               endif
                IF     (PM(I,J).GT.PKM1(I,J)) THEN
                  SIGMAA = ZERO
                  DO KK= 1,KZL
                    ZTOP   =          MAX(ZIL(KK-1),PKM1(I,J))
                    ZBOT   = MAX(ZTOP,MIN(ZIL(KK  ),PM(  I,J)))
                    SIGMAA = SIGMAA + RZ(KK,I,J)*(ZBOT-ZTOP)
                  ENDDO
                  RM(I,J) = MAX(RKM1(I,J),SIGMAA/(PM(I,J)-PKM1(I,J)))
C               ELSE
C                 RM UNCHANGED
                ENDIF
              ENDIF
C
C             FIND T & S.
C
              IF     (PM(I,J).GT.PKM1(I,J)) THEN
                SIGMAA = ZERO
                DO KK= 1,KZL
                  ZTOP   =          MAX(ZIL(KK-1),PKM1(I,J))
                  ZBOT   = MAX(ZTOP,MIN(ZIL(KK  ),PM(  I,J)))
                  SIGMAA = SIGMAA + TZ(KK,I,J)*(ZBOT-ZTOP)
                ENDDO
                TM(I,J) = SIGMAA/(PM(I,J)-PKM1(I,J))
C             ELSE
C               TM UNCHANGED
              ENDIF
              SM(I,J) = SOFSIG(R8(RM(I,J)),R8(TM(I,J)))
            ENDIF  !DEPTH>0
          ENDDO  !J=1,JDM
        ENDDO  !I=1,IDM
C
C       STATISTICS
C
        DO J= 1,JDM
          DO I= 1,IDM
            IF     (DEPTH(I,J).LE.ZERO) THEN
                RM(I,J) = SIG3D(I,J,K)
              WORK(I,J) = 0.0
            ELSE
              WORK(I,J) = PM(I,J) - PKM1(I,J)
            ENDIF
          ENDDO
        ENDDO
        CALL LAYSTAT(TM,  WORK,IDM,JDM, TMIN,XAVE,XMAX)
        WRITE(6,8100) ' tem', TMIN,XAVE,XMAX, K,SIGMA(K)
        CALL LAYSTAT(SM,  WORK,IDM,JDM, SMIN,XAVE,XMAX)
        WRITE(6,8100) ' sal', SMIN,XAVE,XMAX, K,SIGMA(K)
        CALL LAYSTAT(RM,  WORK,IDM,JDM, XMIN,XAVE,XMAX)
        WRITE(6,8100) ' den', XMIN,XAVE,XMAX, K,SIGMA(K)
        CALL LAYSTAT(PKM1, WORK,IDM,JDM, PMIN,XAVE,XMAX)
        WRITE(6,8100) ' inf', PMIN,XAVE,XMAX, K,SIGMA(K)
        CALL LAYSTAT(WORK, WORK,IDM,JDM, XMIN,XAVE,XMAX)
        WRITE(6,8100) ' thk', XMIN,XAVE,XMAX, K,SIGMA(K)
C
C       DIAGNOSTIC PRINTOUT.
C
        IF     (MIN(ITEST,JTEST).GT.0) THEN
          WRITE(6,*)
          WRITE(6,'(A,2I5,I3,A,F9.2,F8.2,A,3F7.3)')
     +     'I,J,K =',ITEST,JTEST,K,
     +     '   INF,THK =',PKM1(ITEST,JTEST),
     +                    WORK(ITEST,JTEST),
     +     '   R,T,S =',    RM(ITEST,JTEST),
     +                      TM(ITEST,JTEST),
     +                      SM(ITEST,JTEST)
          WRITE(6,*)
          CALL ZHFLSH(6)
        ENDIF
C
C       ZONAL INTERFACE DEPTHS.
C
        IF     (K.GT.1 .AND. JDW.NE.0) THEN
          XAVE = ZERO
          DO J= 1,JDM-1
            CALL LAYSTAT(PLAT(1,J), WORK(1,J),IDM,MIN(JDW,JDM-J),
     +                   XMIN,PAVE,XMAX)
            CALL LAYSTAT(PKM1(1,J), WORK(1,J),IDM,MIN(JDW,JDM-J),
     +                   XMIN,XAVE,XMAX)
            WRITE(6,8200) K-1,J,MIN(J+JDW-1,JDM-1),
     +                    SIGMA(K-1),PAVE,XAVE
          ENDDO
        ENDIF
C
C       WRITE OUT CLIMS.
C
        CALL ZAIOWR(TM,MSK,.TRUE., XMIN,XMAX, 10, .FALSE.)
        WRITE(10,4102) ' tem',MONTH,K,SIGMA(K),XMIN,XMAX
C
        CALL ZAIOWR(SM,MSK,.TRUE., XMIN,XMAX, 11, .FALSE.)
        WRITE(11,4102) ' sal',MONTH,K,SIGMA(K),XMIN,XMAX
C
C       ONEM IN HYCOM MKS PRESSURE UNITS IS 9806.
C
        DO J= 1,JDM
          DO I= 1,IDM
            IF     (DEPTH(I,J).GT.ZERO) THEN
              WORK(I,J) = PKM1(I,J)*9806.0
            ENDIF
          ENDDO
        ENDDO
        CALL ZAIOWR(WORK,MSK,.TRUE., XMIN,XMAX, 12, .FALSE.)
        WRITE(12,4102) ' int',MONTH,K,SIGMA(K),XMIN,XMAX
C
C       WRITE OUT DUMMY ARCHIVE.
C
        DO J= 1,JDM
          DO I= 1,IDM
            WORK(I,J) = ZERO
          ENDDO
        ENDDO
        CALL ZAIOWR(WORK,MSK,.FALSE., XMIN,XMAX, 21, .FALSE.)
        WRITE(21,4201) 'u-vel.  ',MONTH,TIME,K,SIGMA(K),XMIN,XMAX
        CALL ZAIOWR(WORK,MSK,.FALSE., XMIN,XMAX, 21, .FALSE.)
        WRITE(21,4201) 'v-vel.  ',MONTH,TIME,K,SIGMA(K),XMIN,XMAX
        DO J= 1,JDM
          DO I= 1,IDM
            WORK(I,J) = (PM(I,J) - PKM1(I,J))*9806.0
          ENDDO
        ENDDO
        CALL ZAIOWR(WORK,MSK,.TRUE.,  XMIN,XMAX, 21, .FALSE.)
        WRITE(21,4201) 'thknss  ',MONTH,TIME,K,SIGMA(K),XMIN,XMAX
        CALL ZAIOWR(TM,  MSK,.TRUE.,  XMIN,XMAX, 21, .FALSE.)
        WRITE(21,4201) 'temp    ',MONTH,TIME,K,SIGMA(K),XMIN,XMAX
        CALL ZAIOWR(SM,  MSK,.TRUE.,  XMIN,XMAX, 21, .FALSE.)
        WRITE(21,4201) 'salin   ',MONTH,TIME,K,SIGMA(K),XMIN,XMAX
        DO J= 1,JDM
          DO I= 1,IDM
            WORK(I,J) = RM(I,J) - THBASE
          ENDDO
        ENDDO
        CALL ZAIOWR(WORK,MSK,.TRUE.,  XMIN,XMAX, 21, .FALSE.)
        WRITE(21,4201) 'density ',MONTH,TIME,K,SIGMA(K),XMIN,XMAX
C
        WRITE(6,6300) K,SIGMA(K)
        CALL ZHFLSH(6)
C
C       PREPARE FOR NEXT K-LOOP
C
        DO J= 1,JDM
          DO I= 1,IDM
            IF     (DEPTH(I,J).GT.ZERO) THEN
              PKM2(I,J) = PKM1(I,J)
              PKM1(I,J) = PM(  I,J)
              RKM1(I,J) = RM(  I,J)
            ENDIF
          ENDDO
        ENDDO
C
 810  CONTINUE  !K=1,KDM
C
      CLOSE (UNIT=10)
      CLOSE (UNIT=11)
      CLOSE (UNIT=12)
      CLOSE (UNIT=21)
      CALL ZAIOCL(10)
      CALL ZAIOCL(11)
      CALL ZAIOCL(12)
      CALL ZAIOCL(21)
      STOP
C
 4000 FORMAT('Expt ',I2.2,'.',I1.1,
     +       '  nhybrd=',I2,
     +        ' nsigma=',I2,
     +          ' ds00=',F5.2,
     +          ' dp00=',F5.2,
     +         ' dp00x=',F6.1,
     +         ' dp00f=',F5.3)
 4101 FORMAT(A79)
 4102 FORMAT(A,': month,layer,dens,range = ',I2.2,I4.2,F7.3,1P2E16.7)
 4200 FORMAT(
     + 'Dummy HYCOM archive from climatology for month ',I2.2,'.' /
     + A80/A80/
     + '1234567890123456789012345678901234567890',
     + '1234567890123456789012345678901234567890'/
     & I5,4X,'''iversn'' = hycom version number x10'/
     & I5,4X,'''iexpt '' = experiment number x10'/
     & I5,4x,'''yrflag'' = days in year flag'/
     & I5,4x,'''idm   '' = longitudinal array size'/
     & I5,4x,'''jdm   '' = latitudinal  array size'/
     & 'field       time step  model day',
     & '  k  dens        min              max')
 4201 FORMAT(a8,' =',i11,f11.2,i3,f7.3,1p2e16.7)
 5000 FORMAT(A40)
 5500 FORMAT(6E13.6)
 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 LAYER',I3,'     SIGMA =',F7.3 /)
 8100 FORMAT(1X,A,': min=',F9.2,' ave=',F9.2,' max=',F9.2,
     +   '   (k,sigma =',i3,F7.2,')')
 8200 FORMAT(' k =',I3.2,' j = ',I4.4,' to ',I4.4,
     +       ' sig =',F7.3,' lat =',F6.1,' inf =',F8.2)
C     END OF PROGRAM WNDINT.
      END
      SUBROUTINE LAYSTAT(PM, THICK, IDM,JDM, PMIN,PAVE,PMAX)
      IMPLICIT NONE
C
      INTEGER IDM,JDM
      REAL*4  PM(IDM,JDM),THICK(IDM,JDM), PMIN,PAVE,PMAX
C
C --- CALCULATE STATISTICS FOR PM.
C --- ONLY WHERE LAYER THICKNESS IS AT LEAST 10 CM.
C --- AVERAGE DOES NOT ALLOW FOR VARIATION IN GRID CELL SIZE.
C
      REAL*4     TENCM
      PARAMETER (TENCM=0.1)
C
      INTEGER I,J,IJSUM
      REAL*4  PPMIN,PPMAX
      REAL*8  PPSUM
C
      PPMIN =  1.E10
      PPMAX = -1.E10
      IJSUM =  0
      PPSUM =  0.0D0
      DO J= 1,JDM
        DO I= 1,IDM
          IF     (THICK(I,J).GT.TENCM) THEN
            PPMIN = MIN( PPMIN, PM(I,J) )
            PPMAX = MAX( PPMAX, PM(I,J) )
            IJSUM = IJSUM + 1
            PPSUM = PPSUM + PM(I,J)
          ENDIF
        ENDDO
      ENDDO
      IF     (IJSUM.NE.0) THEN
        PMIN = PPMIN
        PMAX = PPMAX
        PAVE = PPSUM / IJSUM
      ELSE
        PMIN = 99.9
        PMAX = 99.9
        PAVE = 99.9
      ENDIF
      RETURN
C     END OF LAYSTAT.
      END
