c
c --- profout first (before PROGRAM) to fix an AIX xlf95 bug.
c
      subroutine profout(array,zz, plonj,platj,distj,jn,
     &                   artype,yrflag,time3,nstep,iexpt,
     &                   name,namec,names,units, kz, ncfile)
      use netcdf   ! NetCDF Interface
      
      implicit none
c
      character*(*)    name,namec,names,units,ncfile
      integer          jn,artype,yrflag,nstep,iexpt,kz
      double precision time3(3)
      real             array(kz,jn),zz(kz),
     &                 plonj(jn),platj(jn),distj(jn)
c
c     write out a 2-d z-level array to a netCDF file.
c
c     2-d array size      must be identical in all calls.
c     artype,yrflag,time3 must be identical in all calls.
c
c     the NetCDF filename is taken from ncfile.
c     the NetCDF title, institution and section_name (if any) are taken
c     from environment variables CDF_TITLE, CDF_INST and CDF_SECTION.
c     If environment variables MERSEA_B_DATE and MERSEA_B_TYPE exist,
c     then the output format is in MERSEA format (i.e. no time axis).
c
c     This version needs version 3.5 of the NetCDF library, from: 
c     http://www.unidata.ucar.edu/packages/netcdf/
c
      real             :: array_t(jn,kz)  !transposed array
      real             :: array_0(jn)     !surface array
c
      integer          :: ncfileID, status, varID, latVarID,lonVarID
      integer          :: StaDimID,StaVarID,
     &                    DepDimID,DepVarID
      integer          :: MTDimID,MTVarID,datVarID
      character        :: ncenv*240
      character        :: Ename*6
c
      logical          :: lopen,lexist,lmersea,ltrack,lstation,llon,llat
      integer          :: i,j,k,l,iyear,month,iday,ihour,
     &                            iyrms,monms,idms,ihrms
      real             :: hmin(999),hmax(999)
      double precision :: dt,dt0,year
c
      character*81,     save :: label   = ' '
      integer,          save :: ioinit  = -1
      double precision, save :: date    = 0.d0
      double precision, save :: cell    = 0.d0
      real,        parameter :: fill_value = 2.0**100
c
      character cmonth(12)*3
      data      cmonth/'jan','feb','mar','apr','may','jun',
     &                 'jul','aug','sep','oct','nov','dec'/
c
      ncenv = ' '
      call getenv('MERSEA_B_TYPE',ncenv)
      lmersea  = ncenv.ne.' '
      lstation = distj(jn).eq.0.0
      ltrack   = distj(jn).gt.0.0
      llon     = distj(jn).lt.0.0 .and. distj(jn).gt.-1.5
      llat     = distj(jn).lt.0.0 .and. distj(jn).le.-1.5
*     write(6,*) 
*     write(6,*) 'MERSEA_B_TYPE: ',trim(ncenv)
*     write(6,*) 'lmersea  = ',lmersea
*     write(6,*) 'lstation = ',lstation
*     write(6,*) 'ltrack   = ',ltrack
*     write(6,*) 'llon     = ',llon
*     write(6,*) 'llat     = ',llat
*     write(6,*) 'jn,kz    = ',jn,kz
*     write(6,*) 
c
      if     (ioinit.eq.-1) then
        ioinit = 0
c
c       initialize label.
c
        if     (yrflag.eq.0) then
          year  = 360.0d0
        elseif (yrflag.lt.3) then
          year  = 366.0d0
        else
          year  = 365.25d0
        endif
        call fordate(time3(3),yrflag, iyear,month,iday,ihour)
        date    = (iday + 100 * month + 10000 * iyear) +
     &            (time3(3)-int(time3(3)))
        if     (artype.eq.1) then
          write (label(51:72),123) cmonth(month),iday,iyear,ihour
        elseif (artype.eq.2 .and. time3(2)-time3(1).lt.1.1) then  !daily mean
          write (label(51:72),223) cmonth(month),iday,iyear
        else  ! mean or sdev archive
          write(6,*) 'time3 = ',time3
          dt = 0.5*(time3(2)-time3(1))/(nstep-1)
          if     (yrflag.eq.0) then
            dt0 = 15.0
          elseif (yrflag.eq.1) then
            dt0 = 15.25
          elseif (yrflag.eq.2) then
            dt0 = 0.0
          else
            dt0 = 0.0
          endif
          cell = (time3(2)+dt+dt0) - (time3(1)-dt+dt0)
          if     (artype.eq.2) then
            write(label(51:72),114) ' mean: ',(time3(1)-dt+dt0)/year,
     &                                        (time3(2)+dt+dt0)/year
          else
            write(label(51:72),114) ' sdev: ',(time3(1)-dt+dt0)/year,
     &                                        (time3(2)+dt+dt0)/year
          endif
        endif
        write (label(73:81),115) iexpt/10,mod(iexpt,10),'H'
 123    format ('   ',a3,i3.2,',',i5.4,i3.2,'Z   ')
 223    format ('   ',a3,i3.2,',',i5.4,' MEAN  ')
 114    format (a7,f7.2,'-',f7.2)
 115    format (' [',i2.2,'.',i1.1,a1,']')
      endif  !initialization
c
c       NetCDF I/O
c
        call ncrange(array, 1,jn,kz, fill_value, hmin(1),hmax(1))
c
        inquire(file= ncfile, exist=lexist)
        if (.not.lexist) then
          ! create a new NetCDF and write data to it
          ! netcdf-4 classic model, netcdf version 4.3 and later
          call nchek('nf90_create',
     &                nf90_create(trim(ncfile),
     &                            or(nf90_clobber,
     &                               or(nf90_hdf5,
     &                                  nf90_classic_model)),
     &                            ncfileID))
          ! define the dimensions
          if     (.not.lmersea) then
            call nchek("nf90_def_dim",
     &                  nf90_def_dim(ncfileID,
     &                               "MT", nf90_unlimited,MTDimID))
          endif
          if     (llon)   then
            call nchek("nf90_def_dim",
     &                  nf90_def_dim(ncfileID,"longitude",jn,StaDimID))
          elseif (llat)   then
            call nchek("nf90_def_dim",
     &                  nf90_def_dim(ncfileID,"latitude", jn,StaDimID))
          elseif (ltrack) then
            call nchek("nf90_def_dim",
     &                  nf90_def_dim(ncfileID,"track",    jn,StaDimID))
          else !station
            call nchek("nf90_def_dim",
     &                  nf90_def_dim(ncfileID,"station",  jn,StaDimID))
          endif
          call nchek("nf90_def_dim",
     &                nf90_def_dim(ncfileID,"depth",  kz,DepDimID))
          ! create the global attributes
          call nchek("nf90_put_att (global)",
     &                nf90_put_att(ncfileID,nf90_global,
     &                             "Conventions",
     &                             "CF-1.0"))
            ncenv = ' '
            call getenv('CDF_TITLE',ncenv)
            if     (ncenv.eq.' ') then
              ncenv = "HYCOM"
            endif
            call nchek("nf90_put_att (global)",
     &                  nf90_put_att(ncfileID,nf90_global,
     &                               "title",
     &                               trim(ncenv)))
            ncenv = ' '
            call getenv('CDF_INST',ncenv)
            if     (ncenv.ne.' ') then
              call nchek("nf90_put_att (global)",
     &                    nf90_put_att(ncfileID,nf90_global,
     &                                 "institution",
     &                                 trim(ncenv)))
            endif
            ncenv = ' '
            call getenv('CDF_SECTION',ncenv)
            if     (ncenv.ne.' ') then
              call nchek("nf90_put_att (global)",
     &                    nf90_put_att(ncfileID,nf90_global,
     &                                 "section_name",
     &                                 trim(ncenv)))
            endif
            if     (artype.eq.1) then
              call nchek("nf90_put_att (global)",
     &                    nf90_put_att(ncfileID,nf90_global,
     &                                 "source",
     &                                 "HYCOM archive file"))
            elseif (artype.eq.2) then
              call nchek("nf90_put_att (global)",
     &                    nf90_put_att(ncfileID,nf90_global,
     &                                 "source",
     &                                 "HYCOM mean archive file"))
            else
              call nchek("nf90_put_att (global)",
     &                    nf90_put_att(ncfileID,nf90_global,
     &                                 "source",
     &                                 "HYCOM std. archive file"))
            endif
            call nchek("nf90_put_att (global)",
     &                  nf90_put_att(ncfileID,nf90_global,
     &                               "experiment",
     &                               label(75:78)))
            call nchek("nf90_put_att (global)",
     &                  nf90_put_att(ncfileID,nf90_global,
     &                               "history",
     &                               "hycom_profile2z_nc"))
            if     (lmersea) then
              if     (artype.eq.2) then
                call nchek("nf90_put_att (global)",
     &                      nf90_put_att(ncfileID,nf90_global,
     &                                   "field_type",
     &                                   "daily average"))
              else
                call nchek("nf90_put_att (global)",
     &                      nf90_put_att(ncfileID,nf90_global,
     &                                   "field_type",
     &                                   "instantaneous"))
              endif
              write(ncenv,
     &          '(i4.4,"-",i2.2,"-",i2.2," ",i2.2,":00:00")')
     &          iyear,month,iday,ihour
              call nchek("nf90_put_att (global)",
     &                    nf90_put_att(ncfileID,nf90_global,
     &                                 "field_date",
     &                                 trim(ncenv)))
              ncenv = ' '
              call getenv('MERSEA_B_DATE',ncenv)
              if     (ncenv.eq.'TODAY') then
                write(ncenv,
     &            '(i4.4,"-",i2.2,"-",i2.2," ",i2.2,":00:00")')
     &            iyear,month,iday,ihour
              endif
              if     (ncenv.ne.' ') then
                read(ncenv,'(i4,1x,i2,1x,i2,1x,i2)')
     &            iyrms,monms,idms,ihrms
                if     (iyrms.lt.iyear) then
                  call nchek("nf90_put_att (global)",
     &                        nf90_put_att(ncfileID,nf90_global,
     &                                     "forecast_type",
     &                                     "forecast"))
                elseif (iyrms.gt.iyear) then
                  call nchek("nf90_put_att (global)",
     &                        nf90_put_att(ncfileID,nf90_global,
     &                                     "forecast_type",
     &                                     "hindcast"))
                else   !iyrms.eq.iyear
                  if     (monms.lt.month) then
                    call nchek("nf90_put_att (global)",
     &                          nf90_put_att(ncfileID,nf90_global,
     &                                       "forecast_type",
     &                                       "forecast"))
                  elseif (monms.gt.month) then
                    call nchek("nf90_put_att (global)",
     &                          nf90_put_att(ncfileID,nf90_global,
     &                                       "forecast_type",
     &                                       "hindcast"))
                  else   !monms.eq.month
                    if     (idms.lt.iday) then
                      call nchek("nf90_put_att (global)",
     &                            nf90_put_att(ncfileID,nf90_global,
     &                                         "forecast_type",
     &                                         "forecast"))
                    elseif (idms.gt.iday) then
                      call nchek("nf90_put_att (global)",
     &                            nf90_put_att(ncfileID,nf90_global,
     &                                         "forecast_type",
     &                                         "hindcast"))
                    else   !idms.eq.iday
                      if     (ihrms.lt.ihour) then
                        call nchek("nf90_put_att (global)",
     &                              nf90_put_att(ncfileID,nf90_global,
     &                                           "forecast_type",
     &                                           "forecast"))
                      elseif (ihrms.gt.ihour) then
                        call nchek("nf90_put_att (global)",
     &                              nf90_put_att(ncfileID,nf90_global,
     &                                           "forecast_type",
     &                                           "hindcast"))
                      else   !ihrms.eq.ihour
                        call nchek("nf90_put_att (global)",
     &                              nf90_put_att(ncfileID,nf90_global,
     &                                           "forecast_type",
     &                                           "nowcast"))
                      endif  !ihrms
                    endif !idms
                  endif !monms
                endif !iyrms
                call nchek("nf90_put_att (global)",
     &                      nf90_put_att(ncfileID,nf90_global,
     &                                   "bulletin_date",
     &                                   trim(ncenv)))
              endif
              ncenv = ' '
              call getenv('MERSEA_B_TYPE',ncenv)
              if     (ncenv.ne.' ') then
                call nchek("nf90_put_att (global)",
     &                      nf90_put_att(ncfileID,nf90_global,
     &                                   "bulletin_type",
     &                                   trim(ncenv)))
              endif
            endif !lmersea
          ! create the variables and attributes
          if     (.not.lmersea) then
            call nchek("nf90_def_var (MT)",
     &                  nf90_def_var(ncfileID,"MT",  nf90_double,
     &                               MTDimID,MTVarID))
            if     (yrflag.eq.0) then
              call nchek("nf90_put_att",
     &                    nf90_put_att(ncfileID,MTVarID,
     &                                 "long_name",
     &                                 "model time"))
              call nchek("nf90_put_att",
     &                    nf90_put_att(ncfileID,MTVarID,
     &                                 "units",
     &                            "days since 0001-01-01 00:00:00"))
              call nchek("nf90_put_att",
     &                    nf90_put_att(ncfileID,MTVarID,
     &                                 "calendar",
     &                                 "360_day"))
            elseif (yrflag.eq.1) then
              call nchek("nf90_put_att",
     &                    nf90_put_att(ncfileID,MTVarID,
     &                                 "long_name",
     &                                 "model time"))
              call nchek("nf90_put_att",
     &                    nf90_put_att(ncfileID,MTVarID,
     &                                 "units",
     &                            "days since 0001-01-16 00:00:00"))
              call nchek("nf90_put_att",
     &                    nf90_put_att(ncfileID,MTVarID,
     &                                 "calendar",
     &                                 "366_day"))
            elseif (yrflag.eq.2) then
              call nchek("nf90_put_att",
     &                    nf90_put_att(ncfileID,MTVarID,
     &                                 "long_name",
     &                                 "model time"))
              call nchek("nf90_put_att",
     &                    nf90_put_att(ncfileID,MTVarID,
     &                                 "units",
     &                            "days since 0001-01-01 00:00:00"))
              call nchek("nf90_put_att",
     &                    nf90_put_att(ncfileID,MTVarID,
     &                                 "calendar",
     &                                 "366_day"))
            else
              call nchek("nf90_put_att",
     &                    nf90_put_att(ncfileID,MTVarID,
     &                                 "long_name",
     &                                 "time"))
              call nchek("nf90_put_att",
     &                    nf90_put_att(ncfileID,MTVarID,
     &                                 "units",
     &                            "days since 1900-12-31 00:00:00"))
              call nchek("nf90_put_att",
     &                    nf90_put_att(ncfileID,MTVarID,
     &                                 "calendar",
     &                                 "standard"))
            endif
            call nchek("nf90_put_att",
     &                  nf90_put_att(ncfileID,MTVarID,
     &                               "axis","T"))
            call nchek("nf90_def_var (Date)",
     &                  nf90_def_var(ncfileID,"Date", nf90_double,
     &                               MTDimID,datVarID))
            call nchek("nf90_put_att",
     &                  nf90_put_att(ncfileID,datVarID,
     &                               "long_name",
     &                               "date"))
            call nchek("nf90_put_att",
     &                  nf90_put_att(ncfileID,datVarID,
     &                               "units",
     &                               "day as %Y%m%d.%f"))
            call nchek("nf90_put_att",
     &                  nf90_put_att(ncfileID,datVarID,
     &                               "C_format",
     &                               "%13.4f"))
            call nchek("nf90_put_att",
     &                  nf90_put_att(ncfileID,datVarID,
     &                               "FORTRAN_format",
     &                               "(f13.4)"))
            if     (artype.eq.2) then
              call nchek("nf90_put_att",
     &                    nf90_put_att(ncfileID,MTVarID,
     &                                 "cell_methods",
     &                                 "mean"))
              call nchek("nf90_put_att",
     &                    nf90_put_att(ncfileID,MTVarID,
     &                                 "cell_extent",
     &                                 cell))
              call nchek("nf90_put_att",
     &                    nf90_put_att(ncfileID,datVarID,
     &                                 "cell_methods",
     &                                 "mean"))
              call nchek("nf90_put_att",
     &                    nf90_put_att(ncfileID,datVarID,
     &                                 "cell_extent",
     &                                 cell))
            elseif (artype.eq.3) then
              call nchek("nf90_put_att",
     &                    nf90_put_att(ncfileID,MTVarID,
     &                                 "cell_methods",
     &                                 "standard_deviation"))
              call nchek("nf90_put_att",
     &                    nf90_put_att(ncfileID,MTVarID,
     &                                 "cell_extent",
     &                                 cell))
              call nchek("nf90_put_att",
     &                    nf90_put_att(ncfileID,datVarID,
     &                                 "cell_methods",
     &                                 "standard_deviation"))
              call nchek("nf90_put_att",
     &                    nf90_put_att(ncfileID,datVarID,
     &                                 "cell_extent",
     &                                 cell))
            endif
          endif !.not.lmersea
          call nchek("nf90_def_var (depth)",
     &                nf90_def_var(ncfileID,"depth", nf90_float,
     &                             DepDimID,DepVarID))
          call nchek("nf90_put_att",
     &                nf90_put_att(ncfileID,DepVarID,
     &                             "units","m"))
          call nchek("nf90_put_att",
     &                nf90_put_att(ncfileID,DepVarID,
     &                             "positive","down"))
          call nchek("nf90_put_att",
     &                nf90_put_att(ncfileID,DepVarID,
     &                             "axis","Z"))
          if     (ltrack) then
            call nchek("nf90_def_var (track)",
     &                  nf90_def_var(ncfileID,"track",    nf90_float,
     &                               StaDimID,StaVarID))
            call nchek("nf90_put_att",
     &                  nf90_put_att(ncfileID,StaVarID,
     &                               "units","km"))
            if     (jn.ne.1 .and.
     &              abs( (distj(jn)-distj(1))-
     &                   (distj( 2)-distj(1))*(jn-1) ).lt.0.01) then
              call nchek("nf90_put_att",
     &                    nf90_put_att(ncfileID,StaVarID,
     &                                 "point_spacing","even"))  !ferret
            endif
          elseif (lstation) then
            call nchek("nf90_def_var (station)",
     &                  nf90_def_var(ncfileID,"station",  nf90_int,
     &                               StaDimID,StaVarID))
          endif
          ! longitude and latitude
          call nchek("nf90_def_var (latitude)",
     &                nf90_def_var(ncfileID,"latitude",  nf90_float,
     &                             StaDimID,latVarID))
          call nchek("nf90_put_att",
     &                nf90_put_att(ncfileID,latVarID,
     &                             "standard_name","latitude"))
          call nchek("nf90_put_att",
     &                nf90_put_att(ncfileID,latVarID,
     &                             "units","degrees_north"))
          if     (llat) then
            call nchek("nf90_put_att",
     &                  nf90_put_att(ncfileID,latVarID,
     &                               "axis","Y"))
          endif
          call nchek("nf90_put_att",
     &                nf90_put_att(ncfileID,latVarID,
     &                             "C_format","%8.3f"))
          call nchek("nf90_put_att",
     &                nf90_put_att(ncfileID,latVarID,
     &                             "FORTRAN_format","(f8.3)"))
          call nchek("nf90_def_var (longitude)",
     &                nf90_def_var(ncfileID,"longitude", nf90_float,
     &                             StaDimID,lonVarID))
          call nchek("nf90_put_att",
     &                nf90_put_att(ncfileID,lonVarID,
     &                             "standard_name","longitude"))
          call nchek("nf90_put_att",
     &                nf90_put_att(ncfileID,lonVarID,
     &                             "units","degrees_east"))
          if     (llon) then
            call nchek("nf90_put_att",
     &                  nf90_put_att(ncfileID,lonVarID,
     &                               "axis","X"))
          endif
          call nchek("nf90_put_att",
     &                nf90_put_att(ncfileID,lonVarID,
     &                             "C_format","%8.3f"))
          call nchek("nf90_put_att",
     &                nf90_put_att(ncfileID,lonVarID,
     &                             "FORTRAN_format","(f8.3)"))
          ! model 2Z variable
          if     (kz.gt.1) then
            if     (lmersea) then
*             write(6,*) 'StaDimID = ',StaDimID
*             write(6,*) 'DepDimID = ',DepDimID
              call nchek("nf90_def_var (namec)",
     &                    nf90_def_var(ncfileID,trim(namec),nf90_float,
     &                                 (/StaDimID, DepDimID/),
     &                                 varID))
              if     (llon) then
                call nchek("nf90_put_att",
     &                      nf90_put_att(ncfileID,varID,
     &                                   "coordinates",
     &                                   "latitude"))
              elseif (llat) then
                call nchek("nf90_put_att",
     &                      nf90_put_att(ncfileID,varID,
     &                                   "coordinates",
     &                                   "longitude"))
              else !track/station
                call nchek("nf90_put_att",
     &                      nf90_put_att(ncfileID,varID,
     &                                   "coordinates",
     &                                   "longitude latitude"))
              endif !lon:lat:else
            else !.not.lmersea
              call nchek("nf90_def_var (namec)",
     &                    nf90_def_var(ncfileID,trim(namec),nf90_float,
     &                                 (/StaDimID,
     &                                   DepDimID, MTDimID/),
     &                                 varID))
              if     (llon) then
                call nchek("nf90_put_att",
     &                      nf90_put_att(ncfileID,varID,
     &                                   "coordinates",
     &                                   "latitude Date"))
              elseif (llat) then
                call nchek("nf90_put_att",
     &                      nf90_put_att(ncfileID,varID,
     &                                   "coordinates",
     &                                   "longitude Date"))
              else !track/station
                call nchek("nf90_put_att",
     &                      nf90_put_att(ncfileID,varID,
     &                                   "coordinates",
     &                                   "longitude latitude Date"))
              endif !lon:lat:else
            endif !mersea:else
          else !kz.eq.1
            if     (lmersea) then
*             write(6,*) 'StaDimID = ',StaDimID
              call nchek("nf90_def_var (namec)",
     &                    nf90_def_var(ncfileID,trim(namec),nf90_float,
     &                                 StaDimID, varID))
              if     (llon) then
                call nchek("nf90_put_att",
     &                      nf90_put_att(ncfileID,varID,
     &                                   "coordinates",
     &                                   "latitude"))
              elseif (llat) then
                call nchek("nf90_put_att",
     &                      nf90_put_att(ncfileID,varID,
     &                                   "coordinates",
     &                                   "longitude"))
              else !track/station
                call nchek("nf90_put_att",
     &                      nf90_put_att(ncfileID,varID,
     &                                   "coordinates",
     &                                   "longitude latitude"))
              endif !lon:lat:else
            else !.not.lmersea
              call nchek("nf90_def_var (namec)",
     &                    nf90_def_var(ncfileID,trim(namec),nf90_float,
     &                                 (/StaDimID, MTDimID/),
     &                                 varID))
              if     (llon) then
                call nchek("nf90_put_att",
     &                      nf90_put_att(ncfileID,varID,
     &                                   "coordinates",
     &                                   "latitude Date"))
              elseif (llat) then
                call nchek("nf90_put_att",
     &                      nf90_put_att(ncfileID,varID,
     &                                   "coordinates",
     &                                   "longitude Date"))
              else !track/station
                call nchek("nf90_put_att",
     &                      nf90_put_att(ncfileID,varID,
     &                                   "coordinates",
     &                                   "longitude latitude Date"))
              endif !lon:lat:else
            endif !mersea:else
          endif  !kz.gt.1:else
          if     (names.ne." ") then
            call nchek("nf90_put_att",
     &                  nf90_put_att(ncfileID,varID,
     &                               "standard_name",trim(names)))
          endif
          call nchek("nf90_put_att",
     &                nf90_put_att(ncfileID,varID,"units",trim(units)))
          call nchek("nf90_put_att",
     &                nf90_put_att(ncfileID,varID,
     &                             "_FillValue",fill_value))
          call nchek("nf90_put_att",
     &                nf90_put_att(ncfileID,varID,
     &                             "valid_range",
     &                             (/hmin(1), hmax(1)/)))
          if     (artype.eq.1) then
            call nchek("nf90_put_att",
     &                  nf90_put_att(ncfileID,varID,
     &                               "long_name",
     &                               trim(name)//label(73:81)))
          else
            call nchek("nf90_put_att",
     &                  nf90_put_att(ncfileID,varID,
     &                               "long_name",
     &                 trim(name)//label(51:55)//label(73:81)))
          endif
          ! leave def mode
          call nchek("nf90_enddef",
     &                nf90_enddef(ncfileID))
          ! write data into coordinate variables
          if     (.not.lmersea) then
            call nchek("nf90_put_var",
     &                  nf90_put_var(ncfileID,MTVarID, time3(3)))
            call nchek("nf90_put_var",
     &                  nf90_put_var(ncfileID,datVarID,date    ))
          endif
          if     (ltrack) then
            call nchek("nf90_put_var",
     &                  nf90_put_var(ncfileID,StaVarID,distj(:)))
          elseif (lstation) then
            call nchek("nf90_put_var",
     &                  nf90_put_var(ncfileID,StaVarID,
     &                                         (/ (j,j=1,jn) /)))
          endif
          call nchek("nf90_put_var",
     &                nf90_put_var(ncfileID,DepVarID,zz(:)   ))
          call nchek("nf90_put_var",
     &                nf90_put_var(ncfileID,latVarID,platj(:)))
          call nchek("nf90_put_var",
     &                nf90_put_var(ncfileID,lonVarID,plonj(:)))
          ! write to model variable
          if     (kz.gt.1) then
            do j= 1,jn
              do k= 1,kz
                array_t(j,k) = array(k,j)
              enddo
            enddo
            call nchek("nf90_put_var",
     &                  nf90_put_var(ncfileID,varID,array_t(:,:)))
          else
            do j= 1,jn
              array_0(j) = array(1,j)
            enddo
            call nchek("nf90_put_var",
     &                  nf90_put_var(ncfileID,varID,array_0(:)))
          endif
          ! close NetCDF file
          call nchek("nf90_close",
     &                nf90_close(ncfileID))
        else
c
c        Write data to existing NetCDF file
c
          ! open NetCDF file
          call nchek("nf90_open",
     &                nf90_open(trim(ncfile), nf90_write, ncfileID))
          ! get dimension ID's
          if     (.not.lmersea) then
            call nchek("nf90_inq_dimid",
     &                  nf90_inq_dimid(ncfileID,"MT",MTDimID))
          endif
          call nchek("nf90_inq_dimid",
     &                nf90_inq_dimid(ncfileID,"depth",  DepDimID))
          if     (llon) then
            call nchek("nf90_inq_dimid",
     &                  nf90_inq_dimid(ncfileID,"longitude",StaDimID))
          elseif (llat) then
            call nchek("nf90_inq_dimid",
     &                  nf90_inq_dimid(ncfileID,"latitude", StaDimID))
          elseif (ltrack) then
            call nchek("nf90_inq_dimid",
     &                  nf90_inq_dimid(ncfileID,"track",    StaDimID))
          elseif (lstation) then
            call nchek("nf90_inq_dimid",
     &                  nf90_inq_dimid(ncfileID,"station",  StaDimID))
          endif
          !  switch to define mode
          call nchek("nf90_redef",
     &                nf90_redef(ncfileID))
          ! define new variable
          if     (kz.gt.1) then
            if     (lmersea) then
              call nchek("nf90_def_var (namec)",
     &                    nf90_def_var(ncfileID,trim(namec),nf90_float,
     &                                 (/StaDimID, DepDimID/),
     &                                 varID))
              call nchek("nf90_put_att",
     &                    nf90_put_att(ncfileID,varID,
     &                                 "coordinates",
     &                                 "longitude latitude"))
            else !.not.lmersea
              call nchek("nf90_def_var (namec)",
     &                    nf90_def_var(ncfileID,trim(namec),nf90_float,
     &                                 (/StaDimID, DepDimID, MTDimID/),
     &                                 varID))
              call nchek("nf90_put_att",
     &                    nf90_put_att(ncfileID,varID,
     &                                 "coordinates",
     &                                 "longitude latitude Date"))
            endif !lmersea:else
          else !kz.eq.1
            if     (lmersea) then
              call nchek("nf90_def_var (namec)",
     &                    nf90_def_var(ncfileID,trim(namec),nf90_float,
     &                                 StaDimID,
     &                                 varID))
              call nchek("nf90_put_att",
     &                    nf90_put_att(ncfileID,varID,
     &                                 "coordinates",
     &                                 "longitude latitude"))
            else !.not.lmersea
              call nchek("nf90_def_var (namec)",
     &                    nf90_def_var(ncfileID,trim(namec),nf90_float,
     &                                 (/StaDimID, MTDimID/),
     &                                 varID))
              call nchek("nf90_put_att",
     &                    nf90_put_att(ncfileID,varID,
     &                                 "coordinates",
     &                                 "longitude latitude Date"))
            endif !lmersea:else
          endif !kz
          if     (names.ne." ") then
            call nchek("nf90_put_att",
     &                  nf90_put_att(ncfileID,varID,
     &                               "standard_name",trim(names)))
          endif
          call nchek("nf90_put_att",
     &                nf90_put_att(ncfileID,varID,"units",trim(units)))
          call nchek("nf90_put_att",
     &                nf90_put_att(ncfileID,varID,
     &                             "_FillValue",fill_value))
          call nchek("nf90_put_att",
     &                nf90_put_att(ncfileID,varID,
     &                             "valid_range",
     &                             (/hmin(1), hmax(1)/)))
          if     (artype.eq.1) then
            call nchek("nf90_put_att",
     &                  nf90_put_att(ncfileID,varID,
     &                               "long_name",
     &                               trim(name)//label(73:81)))
          else
            call nchek("nf90_put_att",
     &                  nf90_put_att(ncfileID,varID,
     &                               "long_name",
     &                 trim(name)//label(51:55)//label(73:81)))
          endif
          ! leave define mode
          call nchek("nf90_enddef",
     &                nf90_enddef(ncfileID))
          ! inquire variable ID
          call nchek("nf90_inq_varid",
     &                nf90_inq_varid(ncfileID,trim(namec),varID))
          !write values into array
          if     (kz.gt.1) then
            do j= 1,jn
              do k= 1,kz
                array_t(j,k) = array(k,j)
              enddo
            enddo
            call nchek("nf90_put_var",
     &                  nf90_put_var(ncfileID,varID,array_t(:,:)))
          else
            do j= 1,jn
              array_0(j) = array(1,j)
            enddo
            call nchek("nf90_put_var",
     &                  nf90_put_var(ncfileID,varID,array_0(:)))
          endif
          ! close file
          call nchek("nf90_close",
     &                nf90_close(ncfileID))
        endif
        write(6,'(a49,a,2g15.6)') 
     &    trim(name)//label(51:81),':',hmin(1),hmax(1)
*       call flush(6)
      return

      contains

      subroutine nchek(cnf90,status)
      implicit none
c
      character*(*), intent(in) :: cnf90
      integer,       intent(in) :: status
c
c     subroutine to handle NetCDF errors
c
*     write(6,'(a)') trim(cnf90)
*     call flush(6)
*
      if (status /= nf90_noerr) then
        write(6,'(/a)')   'error in profout - from NetCDF library'
        write(6,'(a/)')   trim(cnf90)
        write(6,'(a/)')   trim(nf90_strerror(status))
*       call flush(6)
        stop
      end if
      end subroutine nchek

      end subroutine profout

      subroutine ncrange(h,ii,jj,kk, fill_value, hmin,hmax)
      implicit none
c
      integer, intent(in ) :: ii,jj,kk
      real,    intent(in ) :: h(ii,jj,kk),fill_value
      real,    intent(out) :: hmin,hmax
c
c     return range of array, ignoring fill_value
c
      integer i,j,k
      real    hhmin,hhmax
c
      hhmin =  abs(fill_value)
      hhmax = -abs(fill_value)
      do k= 1,kk
        do j= 1,jj
          do i= 1,ii
            if     (h(i,j,k).ne.fill_value) then
              hhmin = min(hhmin,h(i,j,k))
              hhmax = max(hhmax,h(i,j,k))
            endif
          enddo
        enddo
      enddo
      hmin = hhmin
      hmax = hhmax
      end subroutine ncrange

      PROGRAM PROFILE
      IMPLICIT NONE
C
C  hycom_profile2z_nc - Usage:  hycom_profile2z_nc archv.a list.txt z.txt archz.nc [itype [bot]]
C
C                 produce a netCDF z-level profile file at a list of points
C
C   archv.a   is assumed to be an HYCOM archive data file, with companion
C             header file archv.b.
C   list.txt  contains a list of points,   one per line
C   z.txt     contains a list of z-depths, one per line
C   archz.nc  will be the output z-level NetCDF profile file
C   itype     is the interpolation type (default 1)
C                =-2; piecewise quadratic cross cell
C                =-1; piecewise linear   across cell
C                = 0; piecewise constant across cell
C                = 1; linear interpolation between cell centers
C   bot       ignore layers within bot of the bottom (default 0.0)
C
C  Each line of list.txt should contain three values: x,y,dist.
C  The x,y values are the p-grid location location of the point,
C  note that hycom/bin/hycom_lonlat2xy will convert lon,lat to x,y.
C  The dist value is the distance from the first point in km,
C  or 0.0 if the points do not represent a track or transect.
C  To use lon (lat) as the standard axis set dist to -1.0 (-2.0).
C
C  The files regional.grid.[ab] for this domain must be in the current
C  directory.  They are used to calculate lon,lat values for each point.
C
C  The following environment variables control the netCDF output:
C
C     CDF_TITLE     - NetCDF title
C     CDF_INST      - NetCDF institution
C     CDF_SECTION   - NetCDF section_name
C     MERSEA_B_DATE - MERSEA Bulletin date
C     MERSEA_B_TYPE - MERSEA Bulletin type
C
C  If environment variables MERSEA_B_DATE and MERSEA_B_TYPE exist,
C  then the output format is in MERSEA format (i.e. no time axis).
C
C  Primarily for MERSEA, so only T,S,u,v (and fsd,sst,sss) are output.
C
C  this version for "serial" Unix systems.
C
C  Alan J. Wallcraft,  Naval Research Laboratory,  July 2003.
C
      REAL*4     ONEM,SPVAL
      PARAMETER (ONEM=9806.0, SPVAL=2.0**100)
C
      REAL*4, ALLOCATABLE :: A(:,:),AA(:,:)
      REAL*4, ALLOCATABLE :: XP(:),YP(:),DP(:)
      REAL*4, ALLOCATABLE :: PLON(:),PLAT(:),PANG(:)
      REAL*4, ALLOCATABLE :: P(:,:),TK(:,:),SK(:,:),
     +                       UK(:,:),VK(:,:),
     +                       UB(:),VB(:),FSD(:),SST(:),SSS(:)
      REAL*4, ALLOCATABLE :: ZZ(:),TZ(:,:),SZ(:,:),
     +                       UZ(:,:),VZ(:,:)
      REAL*4              :: PAD(4096)
C
      INTEGER       IARGC
      INTEGER       NARG
      CHARACTER*240 CARG
C
      LOGICAL       LSTERIC,LSEAICE,LFATAL
      INTEGER       IDM,JDM,KDM,KTR,NSURF,NLAY,NDIF,NTRC,
     +              ARTYPE,IEXPT,YRFLAG,SIGVER,I_ARCH(17)
      INTEGER       ITYPE,IDEBUG,KZ,NPAD,N,NX
      REAL*4        DUMZ,FLAG,THBASE,BOT
      INTEGER       NSTEP
      REAL*8        TIME3(3)
      CHARACTER*240 CFILEA,CFILEB,CFILEL,CFILEP,CFILEZ,CFILEVS,CFORMAT
C
      CHARACTER*18 CASN
      LOGICAL      LDONE
      INTEGER      I,II,IP,IPP1,IPP2,J,JJ,JP,JPP1,JPP2,K,LANDF,LANDL,
     +             KREC,KREC0,IOS,NRECL
      REAL*4       A00,A01,A10,A11,AXY,DX,DY
#ifdef CRAY
      INTEGER*8    IU8,IOS8
#endif
C
C     READ ARGUMENTS.
C
      NARG = IARGC()
C
      IF     (NARG.EQ.4) THEN
        CALL GETARG(1,CFILEA)
        CALL GETARG(2,CFILEL)
        CALL GETARG(3,CFILEZ)
        CALL GETARG(4,CFILEP)
        ITYPE  = 1
        IDEBUG = 0
        BOT    = 0.0
      ELSEIF (NARG.EQ.5) THEN
        CALL GETARG(1,CFILEA)
        CALL GETARG(2,CFILEL)
        CALL GETARG(3,CFILEZ)
        CALL GETARG(4,CFILEP)
        CALL GETARG(5,CARG)
        READ(CARG,*) ITYPE
        IDEBUG = 0
        BOT    = 0.0
      ELSEIF (NARG.EQ.6) THEN 
        CALL GETARG(1,CFILEA)
        CALL GETARG(2,CFILEL)
        CALL GETARG(3,CFILEZ)
        CALL GETARG(4,CFILEP)
        CALL GETARG(5,CARG)
        READ(CARG,*) ITYPE
        CALL GETARG(6,CARG)
        READ(CARG,*) BOT
        IDEBUG = 0
      ELSEIF (NARG.EQ.7) THEN !undocumented debug option
        CALL GETARG(1,CFILEA)
        CALL GETARG(2,CFILEL)
        CALL GETARG(3,CFILEZ)
        CALL GETARG(4,CFILEP)
        CALL GETARG(5,CARG)
        READ(CARG,*) ITYPE
        CALL GETARG(6,CARG)
        READ(CARG,*) BOT
        CALL GETARG(7,CARG)
        READ(CARG,*) IDEBUG  !>0 print this station
      ELSE
        WRITE(6,"(2a)") 
     +    'Usage: hycom_profile2z_nc',
     +    ' archv.a list.txt z.txt archz.nc [itype [bot]]'
        CALL EXIT(1)
        STOP
      ENDIF
C
C     EXTRACT MODEL PARAMETERS FROM ".b" FILE.
C
      cfilevs = ' '
      call getenv('ARCHVS',cfilevs)
      if     (cfilevs.eq.' ') then
        CFILEB = CFILEA(1:LEN_TRIM(CFILEA)-1) // 'b'
        CALL READ_BS(CFILEB,
     +               IEXPT,YRFLAG,IDM,JDM,KDM,NSURF,NLAY,NDIF,NTRC,
     +               LSTERIC,ARTYPE,SIGVER,THBASE,TIME3,NSTEP)
        I_ARCH( :) = -1  !turn off i_arch
      else
        CFILEB = CFILEA(1:LEN_TRIM(CFILEA)-1) // 'b'
        CALL READ_BSCS(CFILEB,CFILEVS,I_ARCH,
     +                 IEXPT,YRFLAG,IDM,JDM,KDM,NSURF,NLAY,NDIF,NTRC,
     +                 LSEAICE,LSTERIC,ARTYPE,SIGVER,THBASE,TIME3,NSTEP)
      endif
*     write(6,*) 'IDM = ',idm
*     write(6,*) 'JDM = ',jdm
*     write(6,*) 'KDM = ',kdm
C
      CALL SIG_I(SIGVER)
C
C     FIND THE NUMBER OF PROFILES AND THE NUMBER OF Z-LEVELS.
C
      OPEN(UNIT=31, FILE=CFILEL, FORM='FORMATTED', STATUS='OLD',
     +         IOSTAT=IOS)
      IF     (IOS.NE.0) THEN
        WRITE(6,*) 'Error: can''t open ',TRIM(CFILEL)
        WRITE(6,*) 'ios   = ',ios
        CALL EXIT(3)
        STOP
      ENDIF
C
      DO N= 1,99999
        READ(31,*,IOSTAT=IOS) DX,DY,DUMZ
        IF     (IOS.NE.0) THEN
          EXIT
        ENDIF
      ENDDO !n
      N = N-1
      REWIND(31)
*     write(6,*) 'N = ',n
C
      OPEN(UNIT=32, FILE=CFILEZ, FORM='FORMATTED', STATUS='OLD',
     +         IOSTAT=IOS)
      IF     (IOS.NE.0) THEN
        WRITE(6,*) 'Error: can''t open ',TRIM(CFILEZ)
        WRITE(6,*) 'ios   = ',ios
        CALL EXIT(3)
        STOP
      ENDIF
C
      DO KZ= 1,99999
        READ(32,*,IOSTAT=IOS) DUMZ
        IF     (IOS.NE.0) THEN
          EXIT
        ENDIF
      ENDDO !kz
      KZ = KZ-1
      REWIND(32)
*     write(6,*) 'KZ  = ',kz
C
C     ALLOCATE ARRAYS.
C
      NPAD = 4096 - MOD(IDM*JDM,4096)
      IF     (NPAD.EQ.4096) THEN
        NPAD = 0
      ENDIF
C
      ALLOCATE( A(IDM,JDM), STAT=IOS )
      IF     (IOS.NE.0) THEN
        WRITE(6,*) 'Error in hycom_profile2z_nc: could not allocate ',
     +             IDM*JDM,' words'
        CALL EXIT(2)
        STOP
      ENDIF
      ALLOCATE( AA(IDM+1,JDM), STAT=IOS )
      IF     (IOS.NE.0) THEN
        WRITE(6,*) 'Error in hycom_profile2z_nc: could not allocate ',
     +             (IDM+1)*JDM,' words'
        CALL EXIT(2)
        STOP
      ENDIF
      ALLOCATE(    P(0:KDM,N),
     +            TK(1:KDM,N),
     +            SK(1:KDM,N),
     +            UK(1:KDM,N),
     +            VK(1:KDM,N), 
     +            ZZ(1:KZ   ),
     +            TZ(1:KZ, N),
     +            SZ(1:KZ, N),
     +            UZ(1:KZ, N),
     +            VZ(1:KZ, N), STAT=IOS )
      IF     (IOS.NE.0) THEN
        WRITE(6,*) 'Error in hycom_profile2z_nc: could not allocate ',
     +             '(1:kz,1:n) arrays'
        CALL EXIT(2)
        STOP
      ENDIF
      ALLOCATE(   UB(      N),
     +            VB(      N),
     +           FSD(      N),
     +           SST(      N),
     +           SSS(      N),
     +          PLON(      N),
     +          PLAT(      N),
     +          PANG(      N),
     +            XP(      N),
     +            YP(      N),
     +            DP(      N), STAT=IOS )
      IF     (IOS.NE.0) THEN
        WRITE(6,*) 'Error in hycom_profile2z_nc: could not allocate ',
     +             '(1:n) arrays'
        CALL EXIT(2)
        STOP
      ENDIF
C
C     READ THE Z FILE.
C
      DO K= 1,KZ
        READ(32,*) ZZ(K)
        IF     (IOS.NE.0) THEN
        ENDIF
      ENDDO !k
      CLOSE(32)
C
      IF     (KZ.EQ.1 .AND. ZZ(1).EQ.0.0) THEN
C
C       SURFACE ONLY.
C
        KDM = 1
      ENDIF
C
C     READ PROFILE LOCATIONS
C
      DO I= 1,N
        READ(31,*,IOSTAT=IOS) XP(I),YP(I),DP(I)
        IF     (IOS.NE.0) THEN
          WRITE(6,*) 'Error: can''t read ',TRIM(CFILEL)
          WRITE(6,*) 'ios   = ',ios
          WRITE(6,*) 'i,n   = ',i,n
          CALL EXIT(3)
          STOP
        ENDIF
      ENDDO !i
      CLOSE(31)
C
C     OPEN "regional.grid.a" FILE.
C
      IF     (NPAD.EQ.0) THEN
        INQUIRE( IOLENGTH=NRECL) A
      ELSE
        INQUIRE( IOLENGTH=NRECL) A,PAD(1:NPAD)
      ENDIF
*     write(6,*) 'nrecl = ',nrecl
#ifdef CRAY
#ifdef t3e
      IF     (MOD(NRECL,4096).EQ.0) THEN
        WRITE(CASN,8000) NRECL/4096
 8000   FORMAT('-F cachea:',I4.4,':1:0')
        IU8 = 11
        CALL ASNUNIT(IU8,CASN,IOS8)
        IF     (IOS8.NE.0) THEN
          WRITE(6,*) 'Error: can''t asnunit 11'
          WRITE(6,*) 'ios  = ',ios8
          WRITE(6,*) 'casn = ',casn
          CALL EXIT(5)
          STOP
        ENDIF
      ENDIF
#else
      CALL ASNUNIT(11,'-F syscall -N ieee',IOS)
      IF     (IOS.NE.0) THEN
        WRITE(6,*) 'Error: can''t asnunit 11'
        WRITE(6,*) 'ios = ',ios
        CALL EXIT(5)
        STOP
      ENDIF
#endif
#endif
      OPEN(UNIT=11, FILE="regional.grid.a",
     +         FORM='UNFORMATTED', STATUS='OLD',
     +         ACCESS='DIRECT', RECL=NRECL, IOSTAT=IOS)
      IF     (IOS.NE.0) THEN
        WRITE(6,*) 'Error: can''t open regional.grid.a'
        WRITE(6,*) 'ios   = ',ios
        WRITE(6,*) 'nrecl = ',nrecl
        CALL EXIT(3)
        STOP
      ENDIF
C
      READ(11,REC=1,IOSTAT=IOS) A  !plon
#ifdef ENDIAN_IO
      CALL ENDIAN_SWAP(A,IDM*JDM)
#endif
      IF     (IOS.NE.0) THEN
        WRITE(6,*) 'can''t read regional.grid.a'
        WRITE(6,*) 'ios  = ',ios
        CALL EXIT(5)
        STOP
      ENDIF
      AA(1:IDM,:) = A(:,:)
      DO J= 1,JDM
        A00 = MOD(AA(1,J) - AA(IDM,J),360.0)
        IF     (A00.LT.-180.0) THEN
          A00 = A00 + 360.0
        ELSEIF (A00.GT. 180.0) THEN
          A00 = A00 - 360.0
        ENDIF
        AA(IDM+1,J) = AA(IDM,J) + A00
      ENDDO !j
      DO I= 1,N
        IP = XP(I)
        JP = YP(I)
        DX = XP(I) - IP
        DY = YP(I) - JP
        IF     (DX.EQ.0.0) THEN
          IPP1 = IP
        ELSEIF (DX.EQ.1.0) THEN
          IPP1 = IP+1
          IP   = IPP1
        ELSE
          IPP1 = IP+1
        ENDIF
        IF     (DY.EQ.0.0) THEN
          JPP1 = JP
        ELSEIF (DY.EQ.1.0) THEN
          JPP1 = JP+1
          JP   = JPP1
        ELSE
          JPP1 = JP+1
        ENDIF
        AXY = (1.0-DX)*(1.0-DY)*AA(IP,  JP)   +
     +        (1.0-DX)*     DY *AA(IP,  JPP1) +
     +             DX *(1.0-DY)*AA(IPP1,JP)   +
     +             DX *     DY *AA(IPP1,JPP1)
        PLON(I) = AXY
        PLON(I) = MOD( PLON(I) + 1080.0, 360.0 )
        IF     (PLON(I).GT.180.0) THEN
          PLON(I) = PLON(I) - 360.0
        ENDIF
        IF     (IDEBUG.EQ.I) THEN
          WRITE(6,*) 'debug: i,plon = ',I,PLON(I)
        ENDIF
      ENDDO !i
C
      READ(11,REC=2,IOSTAT=IOS) A  !plat
#ifdef ENDIAN_IO
      CALL ENDIAN_SWAP(A,IDM*JDM)
#endif
      IF     (IOS.NE.0) THEN
        WRITE(6,*) 'can''t read regional.grid.a'
        WRITE(6,*) 'ios  = ',ios
        CALL EXIT(5)
        STOP
      ENDIF
      AA(1:IDM,:) = A(:,:)
      DO J= 1,JDM
        AA(IDM+1,J) = AA(IDM,J)
      ENDDO !j
      DO I= 1,N
        IP = XP(I)
        JP = YP(I)
        DX = XP(I) - IP
        DY = YP(I) - JP
        IF     (DX.EQ.0.0) THEN
          IPP1 = IP
        ELSEIF (DX.EQ.1.0) THEN
          IPP1 = IP+1
          IP   = IPP1
        ELSEIF (IP.EQ.IDM) THEN
          IPP1 = 1
        ELSE
          IPP1 = IP+1
        ENDIF
        IF     (DY.EQ.0.0) THEN
          JPP1 = JP
        ELSEIF (DY.EQ.1.0) THEN
          JPP1 = JP+1
          JP   = JPP1
        ELSE
          JPP1 = JP+1
        ENDIF
        AXY = (1.0-DX)*(1.0-DY)*AA(IP,  JP)   +
     +        (1.0-DX)*     DY *AA(IP,  JPP1) +
     +             DX *(1.0-DY)*AA(IPP1,JP)   +
     +             DX *     DY *AA(IPP1,JPP1)
        PLAT(I) = AXY
        IF     (IDEBUG.EQ.I) THEN
          WRITE(6,*) 'debug: i,plat = ',I,PLAT(I)
        ENDIF
      ENDDO !i
C
      READ(11,REC=9,IOSTAT=IOS) A  !pang
#ifdef ENDIAN_IO
      CALL ENDIAN_SWAP(A,IDM*JDM)
#endif
      IF     (IOS.NE.0) THEN
        WRITE(6,*) 'can''t read regional.grid.a'
        WRITE(6,*) 'ios  = ',ios
        CALL EXIT(5)
        STOP
      ENDIF
      DO I= 1,N
        IP = NINT(XP(I))
        JP = NINT(YP(I))
        IF     (IP.GT.IDM) THEN
          IP = 1 !periodic wrap
        ENDIF
        IF     (JP.GT.JDM) THEN
          JP = JDM
        ENDIF
        PANG(I) = A(IP,JP)  !no interpolation of angle, use nearest point
        IF     (IDEBUG.EQ.I) THEN
          WRITE(6,*) 'debug: i,pang = ',I,PANG(I)
        ENDIF
      ENDDO !i
      CLOSE(UNIT=11)
      DEALLOCATE( AA )
C
C     OPEN ARCHIVE ".a" FILE.
C
      OPEN(UNIT=11, FILE=CFILEA, FORM='UNFORMATTED', STATUS='OLD',
     +         ACCESS='DIRECT', RECL=NRECL, IOSTAT=IOS)
      IF     (IOS.NE.0) THEN
        WRITE(6,*) 'Error: can''t open ',TRIM(CFILEA)
        WRITE(6,*) 'ios   = ',ios
        WRITE(6,*) 'nrecl = ',nrecl
        CALL EXIT(3)
        STOP
      ENDIF
C
C     SURFACE FIELDS (JUST FSD,UB,VB).
C
      IOS = 0
C
      IF     (NLAY.LE.6) THEN
        LFATAL = .FALSE.
        IF     (I_ARCH(2).LT.0) THEN
          READ(11,REC=2,IOSTAT=IOS) A
#ifdef ENDIAN_IO
          CALL ENDIAN_SWAP(A,IDM*JDM)
#endif
        ELSEIF (I_ARCH(2).GT.0) THEN
          READ(11,REC=I_ARCH(2),IOSTAT=IOS) A
#ifdef ENDIAN_IO
          CALL ENDIAN_SWAP(A,IDM*JDM)
#endif
        ELSE
          A(:,:) = 0.0
        ENDIF
        IF     (IOS.NE.0) THEN
          WRITE(6,*) 'can''t read ',TRIM(CFILEA)
          CALL EXIT(4)
          STOP
        ENDIF
        LANDF = 0
        LANDL = 0
        DO I= 1,N
          IP = XP(I)
          JP = YP(I)
          DX = XP(I) - IP
          DY = YP(I) - JP
          IF     (DX.EQ.0.0) THEN !1-d interpolation
            IPP1 = IP
          ELSEIF (DX.EQ.1.0) THEN !1-d interpolation
            IPP1 = IP+1
            IP   = IPP1
          ELSEIF (IP.EQ.IDM) THEN
            IPP1 = 1
          ELSE
            IPP1 = IP+1
          ENDIF
          IF     (DY.EQ.0.0) THEN !1-d interpolation
            JPP1 = JP
          ELSEIF (DY.EQ.1.0) THEN !1-d interpolation
            JPP1 = JP+1
            JP   = JPP1
          ELSE
            JPP1 = JP+1
          ENDIF
          IF     (MAX( A(IP,  JP),
     +                 A(IP,  JPP1),
     +                 A(IPP1,JP),
     +                 A(IPP1,JPP1) ).GE.SPVAL) THEN
C
C           NEAR LAND, IS NEAREST POINT ACTUALLY LAND?
C
            IF     (DX.LE.0.5 .AND. DY.LE.0.5) THEN
              IP = IP
              JP = JP
            ELSEIF (DX.LE.0.5 .AND. DY.GT.0.5) THEN
              IP = IP
              JP = JPP1
            ELSEIF (DX.GT.0.5 .AND. DY.LE.0.5) THEN
              IP = IPP1
              JP = JP
            ELSE ! (DX.GT.0.5 .AND. DY.GT.0.5) THEN
              IP = IPP1
              JP = JPP1
            ENDIF
*           IF     (.TRUE.) THEN !require 4 good points
            IF     (A(IP,JP).GE.SPVAL) THEN
C
C             NEAREST POINT IS LAND.
C
              FSD(I) = SPVAL
              IF     (LANDF.EQ.0) THEN
                LANDF = I
              ENDIF
              LANDL = I
              IF     (IDEBUG.LT.0) THEN
                WRITE(6,'(a,2i6)') 'nearest point is land: ',IP,JP
                LDONE = .FALSE.
                DO JJ= MAX(1,JP-1),MIN(JDM,JP+1)
                  DO II= MAX(1,IP-1),MIN(IDM,IP+1)
                    IF     (A(II,JJ).LT.SPVAL) THEN
                      LDONE = .TRUE.
                      WRITE(6,'(a,2i6)') '  but sea at: ',II,JJ
                    ENDIF
                  ENDDO
                ENDDO
                IF     (.NOT.LDONE) THEN !25-pt search
                  DO JJ= MAX(1,JP-2),MIN(JDM,JP+2)
                    DO II= MAX(1,IP-2),MIN(IDM,IP+2)
                      IF     (A(II,JJ).LT.SPVAL) THEN
                        LDONE = .TRUE.
                        WRITE(6,'(a,2i6)') '  but sea at: ',II,JJ
                      ENDIF
                    ENDDO
                  ENDDO
                ENDIF
              ENDIF !IDEBUG
            ELSE
C
C             USE NEAREST POINT AS THE EXACT VALUE.
C
              IF     (LANDF.NE.0) THEN
                WRITE(6,'(a,i6,a,i6,a,i6,a)')
     +            'WARNING - STATIONS',LANDF,
     +                           ' TO',LANDL,
     +                           ' OF',N,    ' ARE TREATED AS LAND'
                LANDF = 0
              ENDIF
              FSD(I) = A(IP,JP)/9.806
               XP(I) = IP  !for all subsequent interpolations
               YP(I) = JP  !for all subsequent interpolations
            ENDIF
          ELSE  !all 4 points ok
            IF     (LANDF.NE.0) THEN
              WRITE(6,'(a,i6,a,i6,a,i6,a)')
     +          'WARNING - STATIONS',LANDF,
     +                         ' TO',LANDL,
     +                         ' OF',N,    ' ARE TREATED AS LAND'
              LANDF = 0
            ENDIF
            AXY = (1.0-DX)*(1.0-DY)*A(IP,  JP)   +
     +            (1.0-DX)*     DY *A(IP,  JPP1) +
     +                 DX *(1.0-DY)*A(IPP1,JP)   +
     +                 DX *     DY *A(IPP1,JPP1)
            FSD(I) = AXY/9.806
          ENDIF
          IF     (IDEBUG.EQ.I) THEN
            WRITE(6,*) 'debug: i,k,fsd= ',I,K,FSD(I)
          ENDIF
        ENDDO !i
        IF     (LANDF.NE.0) THEN
          WRITE(6,'(a,i6,a,i6,a,i6,a)')
     +      'WARNING - STATIONS',LANDF,
     +                     ' TO',LANDL,
     +                     ' OF',N,    ' ARE TREATED AS LAND'
        ENDIF
        IF     (I_ARCH(11).LT.0) THEN
          READ(11,REC=NSURF-1,IOSTAT=IOS) A
#ifdef ENDIAN_IO
          CALL ENDIAN_SWAP(A,IDM*JDM)
#endif
        ELSEIF (I_ARCH(11).GT.0) THEN
          READ(11,REC=I_ARCH(11),IOSTAT=IOS) A
#ifdef ENDIAN_IO
          CALL ENDIAN_SWAP(A,IDM*JDM)
#endif
        ELSE
          A(:,:) = 0.0
        ENDIF
        IF     (IOS.NE.0) THEN
          WRITE(6,*) 'can''t read ',TRIM(CFILEA)
          CALL EXIT(4)
          STOP
        ENDIF
        DO I= 1,N
          IF     (FSD(I).EQ.SPVAL) THEN
            CYCLE
          ENDIF
          IP = XP(I)
          JP = YP(I)
          DX = XP(I) - IP
          DY = YP(I) - JP
          IF     (DX.EQ.0.0) THEN
            IPP1 = IP
          ELSEIF (DX.EQ.1.0) THEN
            IPP1 = IP+1
            IP   = IPP1
          ELSEIF (IP.EQ.IDM) THEN
            IPP1 = 1
          ELSE
            IPP1 = IP+1
          ENDIF
          IF     (IPP1.EQ.IDM) THEN
            IPP2 = 1
          ELSE
            IPP2 = IPP1+1
          ENDIF
          IF     (DY.EQ.0.0) THEN
            JPP1 = JP
          ELSEIF (DY.EQ.1.0) THEN
            JPP1 = JP+1
            JP   = JPP1
          ELSE
            JPP1 = JP+1
          ENDIF
          IF (A(IP,  JP)  .EQ.SPVAL) A(IP,  JP)   = 0.0
          IF (A(IPP1,JP)  .EQ.SPVAL) A(IPP1,JP)   = 0.0
          IF (A(IPP2,JP)  .EQ.SPVAL) A(IPP2,JP)   = 0.0
          IF (A(IP,  JPP1).EQ.SPVAL) A(IP,  JPP1) = 0.0
          IF (A(IPP1,JPP1).EQ.SPVAL) A(IPP1,JPP1) = 0.0
          IF (A(IPP2,JPP1).EQ.SPVAL) A(IPP2,JPP1) = 0.0
          IF     (ARTYPE.GT.0) THEN
            A00 = 0.5*(A(IP,  JP)   + A(IP+1,JP)   )
            A01 = 0.5*(A(IP,  JPP1) + A(IP+1,JPP1) )
            A10 = 0.5*(A(IPP1,JP)   + A(IPP2,JP)   )
            A11 = 0.5*(A(IPP1,JPP1) + A(IPP2,JPP1) )
          ELSE !p-grid
            A00 =      A(IP,  JP)
            A01 =      A(IP,  JPP1)
            A10 =      A(IPP1,JP)
            A11 =      A(IPP1,JPP1)
          ENDIF !artype
          AXY = (1.0-DX)*(1.0-DY)*A00 +
     +          (1.0-DX)*     DY *A01 +
     +               DX *(1.0-DY)*A10 +
     +               DX *     DY *A11
          UB( I) = AXY
          IF     (A00.EQ.0.0 .OR. A01.EQ.0.0 .OR.
     +            A10.EQ.0.0 .OR. A11.EQ.0.0     ) THEN
*           LFATAL = .TRUE.
*           WRITE(6,*) 'ERROR - STATION ',I,' HAS BAD U',
*    +                 IP,IPP1,IPP2,JP,JPP1
            WRITE(6,*) 'WARNING - STATION ',I,' HAS BAD U',
     +                 IP,IPP1,IPP2,JP,JPP1
            WRITE(6,*) 
     +        'DX,DY = ',DX,DY
            WRITE(6,*) 
     +        'U.0 = ',A(IP,  JP),  A(IPP1,JP),  A(IPP2,JP)
            WRITE(6,*) 
     +        'U.1 = ',A(IP,  JPP1),A(IPP1,JPP1),A(IPP2,JPP1)
          ENDIF
        ENDDO !i
C
        IF     (I_ARCH(12).LT.0) THEN
          READ(11,REC=NSURF,  IOSTAT=IOS) A
#ifdef ENDIAN_IO
          CALL ENDIAN_SWAP(A,IDM*JDM)
#endif
        ELSEIF (I_ARCH(12).GT.0) THEN
          READ(11,REC=I_ARCH(12),IOSTAT=IOS) A
#ifdef ENDIAN_IO
          CALL ENDIAN_SWAP(A,IDM*JDM)
#endif
        ELSE
          A(:,:) = 0.0
        ENDIF
        IF     (IOS.NE.0) THEN
          WRITE(6,*) 'can''t read ',TRIM(CFILEA)
          CALL EXIT(4)
          STOP
        ENDIF
        DO I= 1,N
          IF     (FSD(I).EQ.SPVAL) THEN
            CYCLE
          ENDIF
          IP = XP(I)
          JP = YP(I)
          DX = XP(I) - IP
          DY = YP(I) - JP
          IF     (DX.EQ.0.0) THEN
            IPP1 = IP
          ELSEIF (DX.EQ.1.0) THEN
            IPP1 = IP+1
            IP   = IPP1
          ELSEIF (IP.EQ.IDM) THEN
            IPP1 = 1
          ELSE
            IPP1 = IP+1
          ENDIF
          IF     (DY.EQ.0.0) THEN
            JPP1 = JP
          ELSEIF (DY.EQ.1.0) THEN
            JPP1 = JP+1
            JP   = JPP1
          ELSE
            JPP1 = JP+1
          ENDIF
          IF     (JPP1.NE.JDM) THEN
            JPP2 = JPP1 + 1
          ELSE
            JPP2 = JPP1
          ENDIF
          IF (A(IP,  JP)  .EQ.SPVAL) A(IP,  JP)   = 0.0
          IF (A(IP,  JP+1).EQ.SPVAL) A(IP,  JP+1) = 0.0
          IF (A(IP,  JPP2).EQ.SPVAL) A(IP,  JPP2) = 0.0
          IF (A(IPP1,JP)  .EQ.SPVAL) A(IPP1,JP)   = 0.0
          IF (A(IPP1,JP+1).EQ.SPVAL) A(IPP1,JP+1) = 0.0
          IF (A(IPP1,JPP2).EQ.SPVAL) A(IPP1,JPP2) = 0.0
          IF     (ARTYPE.GT.0) THEN
            A00 = 0.5*(A(IP,  JP)   + A(IP,  JP+1) )
            A01 = 0.5*(A(IP,  JPP1) + A(IP,  JPP2) )
            A10 = 0.5*(A(IPP1,JP)   + A(IPP1,JP+1) )
            A11 = 0.5*(A(IPP1,JPP1) + A(IPP1,JPP2) )
          ELSE !p-grid
            A00 =      A(IP,  JP)
            A01 =      A(IP,  JPP1)
            A10 =      A(IPP1,JP)
            A11 =      A(IPP1,JPP1)
          ENDIF !artype
          AXY = (1.0-DX)*(1.0-DY)*A00 +
     +          (1.0-DX)*     DY *A01 +
     +               DX *(1.0-DY)*A10 +
     +               DX *     DY *A11
          VB( I) = AXY
          P(0,I) = 0.0
          IF     (A00.EQ.0.0 .OR. A01.EQ.0.0 .OR.
     +            A10.EQ.0.0 .OR. A11.EQ.0.0     ) THEN
*           LFATAL = .TRUE.
*           WRITE(6,*) 'ERROR - STATION ',I,' HAS BAD V',
*    +                 IP,IPP1,JP,JPP1,JPP2
            WRITE(6,*) 'WARNING - STATION ',I,' HAS BAD V',
     +                 IP,IPP1,JP,JPP1,JPP2
            WRITE(6,*) 
     +        'DX,DY = ',DX,DY
            WRITE(6,*) 
     +        'V.0 = ',A(IP,  JP),A(IP,  JP+1),A(IP,  JPP2)
            WRITE(6,*) 
     +        'V.1 = ',A(IPP1,JP),A(IPP1,JP+1),A(IPP1,JPP2)
          ENDIF
        ENDDO !i
        IF     (LFATAL) THEN
          CALL EXIT(9)
          STOP
        ENDIF
      ELSE  !mean archive
        DO I= 1,N
          UB( I) = 0.0
          VB( I) = 0.0
          P(0,I) = 0.0
        ENDDO !i
      ENDIF
C
C     LAYER FIELDS.
C
      DO K= 1,KDM
        KREC0 = NSURF+(NLAY+NDIF+NTRC)*(K-1)
c
c ---   u
        IF     (I_ARCH(13).LT.0) THEN
          READ(11,REC=KREC0+1,IOSTAT=IOS) A
#ifdef ENDIAN_IO
          CALL ENDIAN_SWAP(A,IDM*JDM)
#endif
        ELSEIF (I_ARCH(13).GT.0) THEN
          READ(11,REC=I_ARCH(13),IOSTAT=IOS) A
#ifdef ENDIAN_IO
          CALL ENDIAN_SWAP(A,IDM*JDM)
#endif
        ELSE
          A(:,:) = 0.0
        ENDIF
        IF     (IOS.NE.0) THEN
          WRITE(6,*) 'can''t read ',TRIM(CFILEA)
          CALL EXIT(4)
          STOP
        ENDIF
        DO I= 1,N
          IF     (FSD(I).EQ.SPVAL) THEN
            CYCLE
          ENDIF
          IP = XP(I)
          JP = YP(I)
          DX = XP(I) - IP
          DY = YP(I) - JP
          IF     (DX.EQ.0.0) THEN
            IPP1 = IP
          ELSEIF (DX.EQ.1.0) THEN
            IPP1 = IP+1
            IP   = IPP1
          ELSEIF (IP.EQ.IDM) THEN
            IPP1 = 1
          ELSE
            IPP1 = IP+1
          ENDIF
          IF     (IPP1.EQ.IDM) THEN
            IPP2 = 1
          ELSE
            IPP2 = IPP1+1
          ENDIF
          IF     (DY.EQ.0.0) THEN
            JPP1 = JP
          ELSEIF (DY.EQ.1.0) THEN
            JPP1 = JP+1
            JP   = JPP1
          ELSE
            JPP1 = JP+1
          ENDIF
          IF (A(IP,  JP)  .EQ.SPVAL) A(IP,  JP)   = 0.0
          IF (A(IPP1,JP)  .EQ.SPVAL) A(IPP1,JP)   = 0.0
          IF (A(IPP2,JP)  .EQ.SPVAL) A(IPP2,JP)   = 0.0
          IF (A(IP,  JPP1).EQ.SPVAL) A(IP,  JPP1) = 0.0
          IF (A(IPP1,JPP1).EQ.SPVAL) A(IPP1,JPP1) = 0.0
          IF (A(IPP2,JPP1).EQ.SPVAL) A(IPP2,JPP1) = 0.0
          IF     (ARTYPE.GT.0) THEN
            A00 = 0.5*(A(IP,  JP)   + A(IP+1,  JP)   )
            A01 = 0.5*(A(IP,  JPP1) + A(IP+1,  JPP1) )
            A10 = 0.5*(A(IPP1,JP)   + A(IPP2,  JP)   )
            A11 = 0.5*(A(IPP1,JPP1) + A(IPP2,  JPP1) )
          ELSE !p-grid
            A00 =      A(IP,  JP)
            A01 =      A(IP,  JPP1)
            A10 =      A(IPP1,JP)
            A11 =      A(IPP1,JPP1)
          ENDIF !artype
          AXY = (1.0-DX)*(1.0-DY)*A00 +
     +          (1.0-DX)*     DY *A01 +
     +               DX *(1.0-DY)*A10 +
     +               DX *     DY *A11
          UK(K,I) = UB(I) + AXY
        ENDDO !i
c
c ---   v
        IF     (I_ARCH(14).LT.0) THEN
          READ(11,REC=KREC0+2,IOSTAT=IOS) A
#ifdef ENDIAN_IO
          CALL ENDIAN_SWAP(A,IDM*JDM)
#endif
        ELSEIF (I_ARCH(14).GT.0) THEN
          READ(11,REC=I_ARCH(14),IOSTAT=IOS) A
#ifdef ENDIAN_IO
          CALL ENDIAN_SWAP(A,IDM*JDM)
#endif
        ELSE
          A(:,:) = 0.0
        ENDIF
        IF     (IOS.NE.0) THEN
          WRITE(6,*) 'can''t read ',TRIM(CFILEA)
          CALL EXIT(4)
          STOP
        ENDIF
        DO I= 1,N
          IF     (FSD(I).EQ.SPVAL) THEN
            CYCLE
          ENDIF
          IP = XP(I)
          JP = YP(I)
          DX = XP(I) - IP
          DY = YP(I) - JP
          IF     (DX.EQ.0.0) THEN
            IPP1 = IP
          ELSEIF (DX.EQ.1.0) THEN
            IPP1 = IP+1
            IP   = IPP1
          ELSEIF (IP.EQ.IDM) THEN
            IPP1 = 1
          ELSE
            IPP1 = IP+1
          ENDIF
          IF     (DY.EQ.0.0) THEN
            JPP1 = JP
          ELSEIF (DY.EQ.1.0) THEN
            JPP1 = JP+1
            JP   = JPP1
          ELSE
            JPP1 = JP+1
          ENDIF
          IF     (JPP1.NE.JDM) THEN
            JPP2 = JPP1 + 1
          ELSE
            JPP2 = JPP1
          ENDIF
          IF (A(IP,  JP)  .EQ.SPVAL) A(IP,  JP)   = 0.0
          IF (A(IP,  JP+1).EQ.SPVAL) A(IP,  JP+1) = 0.0
          IF (A(IP,  JPP2).EQ.SPVAL) A(IP,  JPP2) = 0.0
          IF (A(IPP1,JP)  .EQ.SPVAL) A(IPP1,JP)   = 0.0
          IF (A(IPP1,JP+1).EQ.SPVAL) A(IPP1,JP+1) = 0.0
          IF (A(IPP1,JPP2).EQ.SPVAL) A(IPP1,JPP2) = 0.0
          IF     (ARTYPE.GT.0) THEN
            A00 = 0.5*(A(IP,  JP)   + A(IP,  JP+1) )
            A01 = 0.5*(A(IP,  JPP1) + A(IP,  JPP2) )
            A10 = 0.5*(A(IPP1,JP)   + A(IPP1,JP+1) )
            A11 = 0.5*(A(IPP1,JPP1) + A(IPP1,JPP2) )
          ELSE !p-grid
            A00 =      A(IP,  JP)
            A01 =      A(IP,  JPP1)
            A10 =      A(IPP1,JP)
            A11 =      A(IPP1,JPP1)
          ENDIF !artype
          AXY = (1.0-DX)*(1.0-DY)*A00 +
     +          (1.0-DX)*     DY *A01 +
     +               DX *(1.0-DY)*A10 +
     +               DX *     DY *A11
          VK(K,I) = VB(I) + AXY
        ENDDO !i
c
c ---   p
        IF     (ABS(ARTYPE).EQ.1) THEN
          IF     (I_ARCH(15).LT.0) THEN
            READ(11,REC=KREC0+3,IOSTAT=IOS) A
#ifdef ENDIAN_IO
            CALL ENDIAN_SWAP(A,IDM*JDM)
#endif
          ELSEIF (I_ARCH(15).GT.0) THEN
            READ(11,REC=I_ARCH(15),IOSTAT=IOS) A
#ifdef ENDIAN_IO
            CALL ENDIAN_SWAP(A,IDM*JDM)
#endif
          ELSE
            A(:,:) = 0.0
          ENDIF
        ELSE  !mean archive
          READ(11,REC=KREC0+4,IOSTAT=IOS) A
#ifdef ENDIAN_IO
          CALL ENDIAN_SWAP(A,IDM*JDM)
#endif
        ENDIF
        IF     (IOS.NE.0) THEN
          WRITE(6,*) 'can''t read ',TRIM(CFILEA)
          CALL EXIT(4)
          STOP
        ENDIF
        DO I= 1,N
          IF     (FSD(I).EQ.SPVAL) THEN
            CYCLE
          ENDIF
          IP = XP(I)
          JP = YP(I)
          DX = XP(I) - IP
          DY = YP(I) - JP
          IF     (DX.EQ.0.0) THEN
            IPP1 = IP
          ELSEIF (DX.EQ.1.0) THEN
            IPP1 = IP+1
            IP   = IPP1
          ELSEIF (IP.EQ.IDM) THEN
            IPP1 = 1
          ELSE
            IPP1 = IP+1
          ENDIF
          IF     (DY.EQ.0.0) THEN
            JPP1 = JP
          ELSEIF (DY.EQ.1.0) THEN
            JPP1 = JP+1
            JP   = JPP1
          ELSE
            JPP1 = JP+1
          ENDIF
          AXY = (1.0-DX)*(1.0-DY)*A(IP,  JP)   +
     +          (1.0-DX)*     DY *A(IP,  JPP1) +
     +               DX *(1.0-DY)*A(IPP1,JP)   +
     +               DX *     DY *A(IPP1,JPP1)
          P(K,I) = P(K-1,I) + AXY/ONEM  !interface depth in m
          IF     (IDEBUG.EQ.I) THEN
            WRITE(6,*) 'debug: i,k,p  = ',I,K,P(K,I)
          ENDIF
        ENDDO !i
c
c ---   t
        IF     (ABS(ARTYPE).EQ.1) THEN
          IF     (I_ARCH(16).LT.0) THEN
            READ(11,REC=KREC0+4,IOSTAT=IOS) A
#ifdef ENDIAN_IO
            CALL ENDIAN_SWAP(A,IDM*JDM)
#endif
          ELSEIF (I_ARCH(16).GT.0) THEN
            READ(11,REC=I_ARCH(16),IOSTAT=IOS) A
#ifdef ENDIAN_IO
            CALL ENDIAN_SWAP(A,IDM*JDM)
#endif
          ELSE
            A(:,:) = 0.0
          ENDIF
        ELSE  !mean archive
          READ(11,REC=KREC0+5,IOSTAT=IOS) A
#ifdef ENDIAN_IO
          CALL ENDIAN_SWAP(A,IDM*JDM)
#endif
        ENDIF
        IF     (IOS.NE.0) THEN
          WRITE(6,*) 'can''t read ',TRIM(CFILEA)
          CALL EXIT(4)
          STOP
        ENDIF
        DO I= 1,N
          IF     (FSD(I).EQ.SPVAL) THEN
            CYCLE
          ENDIF
          IP = XP(I)
          JP = YP(I)
          DX = XP(I) - IP
          DY = YP(I) - JP
          IF     (DX.EQ.0.0) THEN
            IPP1 = IP
          ELSEIF (DX.EQ.1.0) THEN
            IPP1 = IP+1
            IP   = IPP1
          ELSEIF (IP.EQ.IDM) THEN
            IPP1 = 1
          ELSE
            IPP1 = IP+1
          ENDIF
          IF     (DY.EQ.0.0) THEN
            JPP1 = JP
          ELSEIF (DY.EQ.1.0) THEN
            JPP1 = JP+1
            JP   = JPP1
          ELSE
            JPP1 = JP+1
          ENDIF
          AXY = (1.0-DX)*(1.0-DY)*A(IP,  JP)   +
     +          (1.0-DX)*     DY *A(IP,  JPP1) +
     +               DX *(1.0-DY)*A(IPP1,JP)   +
     +               DX *     DY *A(IPP1,JPP1)
          TK(K,I) = AXY
          IF     (IDEBUG.EQ.I) THEN
            WRITE(6,*) 'debug: i,k,tk = ',I,K,TK(K,I)
          ENDIF
        ENDDO !i
c
c ---   s
        IF     (ABS(ARTYPE).EQ.1) THEN
          IF     (I_ARCH(17).LT.0) THEN
            READ(11,REC=KREC0+5,IOSTAT=IOS) A
#ifdef ENDIAN_IO
            CALL ENDIAN_SWAP(A,IDM*JDM)
#endif
          ELSEIF (I_ARCH(17).GT.0) THEN
            READ(11,REC=I_ARCH(17),IOSTAT=IOS) A
#ifdef ENDIAN_IO
            CALL ENDIAN_SWAP(A,IDM*JDM)
#endif
          ELSE
            A(:,:) = 0.0
          ENDIF
        ELSE  !mean archive
          READ(11,REC=KREC0+6,IOSTAT=IOS) A
#ifdef ENDIAN_IO
          CALL ENDIAN_SWAP(A,IDM*JDM)
#endif
        ENDIF
        IF     (IOS.NE.0) THEN
          WRITE(6,*) 'can''t read ',TRIM(CFILEA)
          CALL EXIT(4)
          STOP
        ENDIF
        DO I= 1,N
          IF     (FSD(I).EQ.SPVAL) THEN
            CYCLE
          ENDIF
          IP = XP(I)
          JP = YP(I)
          DX = XP(I) - IP
          DY = YP(I) - JP
          IF     (DX.EQ.0.0) THEN
            IPP1 = IP
          ELSEIF (DX.EQ.1.0) THEN
            IPP1 = IP+1
            IP   = IPP1
          ELSEIF (IP.EQ.IDM) THEN
            IPP1 = 1
          ELSE
            IPP1 = IP+1
          ENDIF
          IF     (DY.EQ.0.0) THEN
            JPP1 = JP
          ELSEIF (DY.EQ.1.0) THEN
            JPP1 = JP+1
            JP   = JPP1
          ELSE
            JPP1 = JP+1
          ENDIF
          AXY = (1.0-DX)*(1.0-DY)*A(IP,  JP)   +
     +          (1.0-DX)*     DY *A(IP,  JPP1) +
     +               DX *(1.0-DY)*A(IPP1,JP)   +
     +               DX *     DY *A(IPP1,JPP1)
          SK(K,I) = AXY
        ENDDO !i
      ENDDO !k
C
C     INTERPOLATE ALL PROFILES TO Z
C
      DO I= 1,N
        IF     (FSD(I).EQ.SPVAL) THEN
          SST( I) = SPVAL
          SSS( I) = SPVAL
          DO K= 1,KZ
            UZ(K,I) = SPVAL
            VZ(K,I) = SPVAL
            TZ(K,I) = SPVAL
            SZ(K,I) = SPVAL
          ENDDO !k
          CYCLE
        ENDIF
C
        IF     (BOT.NE.0.0) THEN
C
C         IGNORE LAYERS WITHIN BOT OF THE BOTTOM
C
          DO K= 2,KDM
            IF      (P(K-1,I).GE.P(KDM,I)-BOT) THEN
              UK(K,I) = UK(K-1,I)
              VK(K,I) = VK(K-1,I)
              TK(K,I) = TK(K-1,I)
              SK(K,I) = SK(K-1,I)
            ENDIF
          ENDDO !k
        ENDIF
C
        SST(I) = TK(1,I)
        SSS(I) = SK(1,I)
C
        FLAG = SPVAL
        IF     (ITYPE.EQ.-2) THEN
          CALL LAYER2Z_PPM(UK(1,I),VK(1,I),TK(1,I),SK(1,I),P(0,I), KDM,
     &                     UZ(1,I),VZ(1,I),TZ(1,I),SZ(1,I),ZZ,KZ, FLAG)
        ELSEIF (ITYPE.EQ.-1) THEN
          CALL LAYER2Z_PLM(UK(1,I),VK(1,I),TK(1,I),SK(1,I),P(0,I), KDM,
     &                     UZ(1,I),VZ(1,I),TZ(1,I),SZ(1,I),ZZ,KZ, FLAG)
        ELSEIF (ITYPE.EQ. 0) THEN
          CALL LAYER2Z_PCM(UK(1,I),VK(1,I),TK(1,I),SK(1,I),P(0,I), KDM,
     &                     UZ(1,I),VZ(1,I),TZ(1,I),SZ(1,I),ZZ,KZ, FLAG)
        ELSEIF (ITYPE.EQ. 1 .AND. IDEBUG.EQ.I) THEN
          CALL LAYER2Z_LIN_DEBUG(
     &                     UK(1,I),VK(1,I),TK(1,I),SK(1,I),P(0,I), KDM,
     &                     UZ(1,I),VZ(1,I),TZ(1,I),SZ(1,I),ZZ,KZ, FLAG)
        ELSEIF (ITYPE.EQ. 1) THEN
          CALL LAYER2Z_LIN(UK(1,I),VK(1,I),TK(1,I),SK(1,I),P(0,I), KDM,
     &                     UZ(1,I),VZ(1,I),TZ(1,I),SZ(1,I),ZZ,KZ, FLAG)
        ELSE
          WRITE(6,"(2a)") 
     +      'Usage: hycom_profile2z_nc',
     +      ' archv.a list.txt z.txt archz.nc [itype [bot]]'
          WRITE(6,*) 'unsupported itype value'
          CALL EXIT(1)
          STOP
        ENDIF
        IF     (PANG(I).NE.0.0) THEN
C
C         ROTATE TO EASTWARDS,NORTHWARDS
C
          DO K= 1,KZ
            IF     (UZ(K,I).NE.SPVAL) THEN
              DX = UZ(K,I)
              DY = VZ(K,I)
              UZ(K,I) = COS(PANG(I))*DX + SIN(-PANG(I))*DY
              VZ(K,I) = COS(PANG(I))*DY - SIN(-PANG(I))*DX
            ENDIF
          ENDDO !k
        ENDIF
        IF     (IDEBUG.EQ.I) THEN
          DO K= 1,KZ
            WRITE(6,*) 'debug: i,k,tz = ',I,K,TZ(K,I)
          ENDDO !k
        ENDIF
      ENDDO !i
C
C     NETCDF OUTPUT.
C
      ARTYPE = ABS(ARTYPE)
      CALL PROFOUT(UZ,ZZ, PLON,PLAT,DP,N,
     &             ARTYPE,YRFLAG,TIME3,NSTEP,IEXPT,
     &             ' u-veloc.',                       ! plot name
     &             'u',                               ! ncdf name (mersea)
     &             'eastward_sea_water_velocity',     ! ncdf standard_name
     &             'm/s',                             ! units
     &             KZ, CFILEP)
      CALL PROFOUT(VZ,ZZ, PLON,PLAT,DP,N,
     &             ARTYPE,YRFLAG,TIME3,NSTEP,IEXPT,
     &             ' v-veloc.',                       ! plot name
     &             'v',                               ! ncdf name (mersea)
     &             'northward_sea_water_velocity',    ! ncdf standard_name
     &             'm/s',                             ! units
     &             KZ, CFILEP)
      IF     (KZ.NE.1 .OR. ZZ(1).NE.0.0) THEN
      CALL PROFOUT(TZ,ZZ, PLON,PLAT,DP,N,
     &             ARTYPE,YRFLAG,TIME3,NSTEP,IEXPT,
     &             '  temp   ',                       ! plot name
     &             'temperature',                     ! ncdf name
     &             'sea_water_potential_temperature', ! ncdf standard_name
     &             'degC',                            ! units
     &             KZ, CFILEP)
      CALL PROFOUT(SZ,ZZ, PLON,PLAT,DP,N,
     &             ARTYPE,YRFLAG,TIME3,NSTEP,IEXPT,
     &             ' salinity',                       ! plot name
     &             'salinity',                        ! ncdf name
     &             'sea_water_salinity',              ! ncdf standard_name
     &             'psu',                             ! units
     &             KZ, CFILEP)
      ENDIF
      CALL PROFOUT(FSD,ZZ, PLON,PLAT,DP,N,
     &             ARTYPE,YRFLAG,TIME3,NSTEP,IEXPT,
     &             '  SSH    ',                       ! plot name
     &             'ssh',                             ! ncdf name (mersea)
     &             'sea_surface_elevation',           ! ncdf standard_name
     &             'm',                               ! units
     &             1, CFILEP)
      CALL PROFOUT(SST,ZZ, PLON,PLAT,DP,N,
     &             ARTYPE,YRFLAG,TIME3,NSTEP,IEXPT,
     &             '  SST    ',                       ! plot name
     &             'sst',                             ! ncdf name
     &             'sea_surface_temperature',         ! ncdf standard_name
     &             'degC',                            ! units
     &             1, CFILEP)
      CALL PROFOUT(SSS,ZZ, PLON,PLAT,DP,N,
     &             ARTYPE,YRFLAG,TIME3,NSTEP,IEXPT,
     &             ' SSS     ',                       ! plot name
     &             'sss',                             ! ncdf name
     &             'sea_surface_salinity',            ! ncdf standard_name
     &             'psu',                             ! units
     &             1, CFILEP)
      END
