program cnv_eps2p

c  Reads secondary sensor data with 3 opticals (no PAR).  Modified for GP806-KA.
c  1/4/07.

c  Reads ASCII CTD profile data (.cnv) created by SEASOFT processing modules,
c  applies post-cruise calibrations to temp and cond as given in the .coef 
c  file, derives standard variables, and writes 1 db data into NetCDF format.
c  Program determines first good salinity value and copies back raw data
c  (P, T, C) to 0 db.  Program skips over initial records containing -9.990e-29 
c  and checks for potential density inversions.  Program removes data where 
c  package slowdowns and reversals contaminate the profile (SBE module LOOPEDIT 
c  does not handle this well enough) and linearly interpolates any gaps such 
c  that a record exists every 1 db.

c  Modified from /home/nansen2/ctd/src/cnv_eps.f for Ka'imimoana upgrade from
c  VAX/VMS to Windows NT with Cygnus UNIX interface.  K. McTaggart, 5/1/01.

c  Using Cygnus/UNIX compile and link with
c  g77 cnv_eps2p.f -o cnv_eps2p.x -L/usr/local/epic/lib -ldmc -lphyprops -leps -lnetcdf -lepic

c  Modified according to nansen2/ctd/src/cnv_eps2.f to run on PC.  KEM, 4/5/05.

c  Modified to read files with raw optical voltages for Pete Strutton.  KEM, 10/26/05.
c  Voltages are not copied back to the surface, but they are interpolated along with
c  the CTD data over missing dbars and density inversions.

       implicit none

c  Include attributes required for use with EPIC programs:

       include '/home/ctd/src/epsystemk.fh'

       integer nx, ny, nz, nt, nvar
       parameter(nx=1, ny=1, nz=6000, nt=1, nvar=11)
       integer itime(2,nt), eptime(2), epfout
       real lat(ny), lon(nx), pres(nz)

       real rawp,rawt,rawc,raws,rawth,rawn,pfill,sal,tht
       real rawv0,rawv1,rawv2,rawv4
       real pr(1,1,nz,1),t068(1,1,nz,1),c0(1,1,nz,1),qual(1,1,nz,1)
       real sa(1,1,nz,1),ths(1,1,nz,1),nbin(1,1,nz,1),topp(10),tops(10)
       real v0(1,1,nz,1),v1(1,1,nz,1),v2(1,1,nz,1),v4(1,1,nz,1)
       real fullp(1,1,nz,1),fullt(1,1,nz,1),fullc(1,1,nz,1)
       real fulln(1,1,nz,1),fullq(1,1,nz,1),sigth(nz)
       real fullv0(1,1,nz,1),fullv1(1,1,nz,1),fullv2(1,1,nz,1)
       real fullv4(1,1,nz,1)
       real t90(1,1,nz,1),scal(1,1,nz,1),sig(1,1,nz,1),dh(1,1,nz,1)
       real ccal(1,1,nz,1),th(1,1,nz,1),sth(1,1,nz,1) 
       real dp,ds,sv,svl,du,sal78,theta,svan,xlat,xlon,dry,wet
       real drift,beta1,bias1,csta(100),slope(100),slopes
       real beta2,bias2,offset   

       integer i, j, jj, k, kk, ndata, nskip, nfill, ifill, nhdr 
       integer idy, imo, iyr, ihr, imin, ilat, ilon
       integer bar, dep, wdir, wspd, ista

       character*1 hlat,hlon,test,dum,state,vis,amt,cloud,weat
       character mon*3, stacst*4, cruise*5, crs_str*10, outfile*13
       character fname*11, station*4

c -----------------------------------------------------------------------

	call getenv('CTD_CRUISE',cruise)
	call getenv('CTD_STRING',crs_str)

c  Open input data and log files

	open (unit=3,file=cruise//'s1_eps.dat',form='formatted',status='old')
	open (unit=10,file=cruise//'s1_eps.log',form='formatted',status='new')

c  Read calibration coefficients

c	open (unit=9,file=cruise//'s.coef',form='formatted',status='old')
c	read (9,14) drift
c  14	format (23x,f9.6)
c	read (9,16) bias1
c	read (9,16) bias2
c	read (9,16) beta1
c	read (9,16) beta2
c  16	format (23x,e15.7)
c	do 18 i=1,100
c	  read (9,17,end=19) csta(i),slope(i)
c  17	  format (f6.1,16x,f9.7)
c  18 	continue

c  Open input .CNV file 

  19	read (3,20,end=900) fname
  20	format (a11)
        read (fname(4:7),30) station
  30    format (a4)
	read (fname(4:6),32) ista
  32    format (i3)
        write (6,35) fname
 	write (10,35) fname
  35    format (' Converting ',a11,' to netCDF format...')

c  Open and read from .HDR file:

       open (unit=7,file=cruise//'s1.hdr',form='formatted',status='old')
  38   read (7,40,end=42) stacst,ilat,xlat,hlat,ilon,xlon,hlon,idy,mon,
     .   iyr,ihr,imin,dry,wet,bar,state,vis,wdir,wspd,amt,cloud,weat,
     .   dep
  40   format (a4,i2,f5.1,a1,i3,f5.1,a1,i2,1x,a3,1x,3i2,2f4.1,i2,2a1,i3,
     .   i2,3a1,i4)

       if (station.eq.stacst) then
         goto 48
       else
         goto 38
       endif

  42   write (6,45) station
  45   format (' Station not found in .HDR: ',a4)
       goto 900

c  Create lat/lon as real numbers:

  48   lat(ny)=ilat + (xlat/60)
       if (hlat.eq.'S') lat(ny)=lat(ny)*-1
       lon(nx)=ilon + (xlon/60)
       if (hlon.eq.'W') lon(nx)=360-lon(nx)

c  Create date/time string:

       if (mon.eq.'Jan') imo=1
       if (mon.eq.'Feb') imo=2
       if (mon.eq.'Mar') imo=3
       if (mon.eq.'Apr') imo=4
       if (mon.eq.'May') imo=5
       if (mon.eq.'Jun') imo=6
       if (mon.eq.'Jul') imo=7
       if (mon.eq.'Aug') imo=8
       if (mon.eq.'Sep') imo=9
       if (mon.eq.'Oct') imo=10
       if (mon.eq.'Nov') imo=11
       if (mon.eq.'Dec') imo=12
       iyr=iyr + 2000
       call mdyhmstoeptime(imo,idy,iyr,ihr,imin,0.0,eptime)
       itime(1,1) = eptime(1)
       itime(2,1) = eptime(2)

c  Open .CNV input data file, skip .CNV file header, and skip bad records near
c  the surface:

       nhdr=0
       open (unit=1,file=fname,form='formatted',status='old')
  54   read (1,55) test
  55   format (a1)
       if (test.eq.'*'.or.test.eq.'#') then
	   nhdr=nhdr+1
	   goto 54
       else
         close (1)
       endif

       open (unit=1,file=fname,form='formatted',status='old')
       do i=1,nhdr
	 read (1,55) dum
       enddo

       do i=1,10
         read (1,56) topp(i),tops(i)
  56     format (22x,f11.3,99x,f11.4)
	 enddo
       close (1)

       do 59 i=1,9
         if (topp(i).lt.15) then
         if (abs(tops(i)-tops(i+1)).gt..005) then
           goto 59
         else
           nskip=i-1
           write (6,58) topp(i),tops(i)
	   write (10,58) topp(i),tops(i)
  58       format (' First good record at',f6.1,f9.4)
           goto 60
         endif
         endif
  59   continue

c  Read and fill raw data arrays (P, T, C, nbin) omitting those records
c  containing -9.990e-29:

  60   open (unit=1,file=fname,form='formatted',status='old')
       do i=1,nhdr+nskip
         read (1,55) dum
	 enddo

       ndata=1
       do 70 k=1,6000
         read (1,65,end=75) rawp,rawt,rawc,rawv0,rawv1,rawv2,
     .     rawn,raws
  65     format (22x,f11.3,11x,f11.4,11x,f11.6,3f11.4,f11.0,11x,f11.4)
         if (rawp.eq.-9.990e-29.or.rawt.eq.-9.990e-29.or.rawc.eq.
     .     -9.990e-29.or.raws.eq.-9.990e-29)
     .     then
           goto 70
         else
           pr(1,1,ndata,1)=rawp
           t068(1,1,ndata,1)=rawt
           c0(1,1,ndata,1)=rawc
           v0(1,1,ndata,1)=rawv0
	   v1(1,1,ndata,1)=rawv1
	   v2(1,1,ndata,1)=rawv2
	   v4(1,1,ndata,1)=0.0
           nbin(1,1,ndata,1)=rawn
	   sa(1,1,ndata,1)=raws
           qual(1,1,ndata,1)=222.0
	   ths(1,1,ndata,1)=theta(raws,rawt,rawp,0.)
           if (pr(1,1,ndata,1).lt.pr(1,1,ndata-1,1)) goto 75
           ndata=ndata+1
         endif
  70   continue

  75   ndata=ndata-1
       write (6,76) pr(1,1,1,1),pr(1,1,ndata,1),ndata
       write (10,76) pr(1,1,1,1),pr(1,1,ndata,1),ndata
  76   format (' Raw data range:',2f8.1,i6)
       close (1)

c  Copy back data within 15 db to 0 db. WOCE flag set to 8 (undefined):
c  Do not copy back voltages; set to zero.

        nfill=0
        if (pr(1,1,1,1).gt.0.and.pr(1,1,1,1).le.15) then
          ifill=pr(1,1,1,1)
          pfill=-1.0
          do 77 i=1,ifill
            nfill=nfill+1
            pfill=pfill+1.0
            fullp(1,1,i,1)=pfill
            fullt(1,1,i,1)=theta(sa(1,1,1,1),ths(1,1,1,1),0,pfill)
            fullc(1,1,i,1)=sal78(sa(1,1,1,1),fullt(1,1,i,1),pfill,
     .        1)*42.914
	    fullv0(1,1,i,1)=0.0
	    fullv1(1,1,i,1)=0.0
	    fullv2(1,1,i,1)=0.0
	    fullv4(1,1,i,1)=0.0
            fulln(1,1,i,1)=0.0
            fullq(1,1,i,1)=888.0
  77      continue
          nfill=nfill+1
        else
          nfill=1
        endif

        fullp(1,1,nfill,1)=pr(1,1,1,1)
        fullt(1,1,nfill,1)=t068(1,1,1,1)
        fullc(1,1,nfill,1)=c0(1,1,1,1)
        fullv0(1,1,nfill,1)=v0(1,1,1,1)
        fullv1(1,1,nfill,1)=v1(1,1,1,1)
        fullv2(1,1,nfill,1)=v2(1,1,1,1)
        fullv4(1,1,nfill,1)=v4(1,1,1,1)
        fulln(1,1,nfill,1)=nbin(1,1,1,1)
        fullq(1,1,nfill,1)=qual(1,1,1,1)
        nfill=nfill+1

c  Continue filling 'full' raw data arrays, linearly interpolating over missing
c  1 db records from the SBE processing modules.  WOCE quality flags are set to
c  666.0 where values are interpolated over >2 db:

        do 98 i=2,ndata
          dp = pr(1,1,i,1)-pr(1,1,i-1,1)
          if (dp.eq.1.0) then
            fullp(1,1,nfill,1)=pr(1,1,i,1)
            fullt(1,1,nfill,1)=t068(1,1,i,1)
            fullc(1,1,nfill,1)=c0(1,1,i,1)
	    fullv0(1,1,nfill,1)=v0(1,1,i,1)
	    fullv1(1,1,nfill,1)=v1(1,1,i,1)
	    fullv2(1,1,nfill,1)=v2(1,1,i,1)
	    fullv4(1,1,nfill,1)=v4(1,1,i,1)
            fulln(1,1,nfill,1)=nbin(1,1,i,1)
            fullq(1,1,nfill,1)=qual(1,1,i,1)
            nfill=nfill+1
          elseif (dp.eq.2.0) then
            do 96 j=1,dp
              fullp(1,1,nfill+j-1,1)=pr(1,1,i-1,1)+j
              fullt(1,1,nfill+j-1,1)=t068(1,1,i-1,1)+(((t068(1,1,i,1)
     .          -t068(1,1,i-1,1))/dp)*j)
              fullc(1,1,nfill+j-1,1)=c0(1,1,i-1,1)+(((c0(1,1,i,1)-
     .          c0(1,1,i-1,1))/dp)*j)
              fullv0(1,1,nfill+j-1,1)=v0(1,1,i-1,1)+(((v0(1,1,i,1)
     .          -v0(1,1,i-1,1))/dp)*j)
              fullv1(1,1,nfill+j-1,1)=v1(1,1,i-1,1)+(((v1(1,1,i,1)
     .          -v1(1,1,i-1,1))/dp)*j)
              fullv2(1,1,nfill+j-1,1)=v2(1,1,i-1,1)+(((v2(1,1,i,1)
     .          -v2(1,1,i-1,1))/dp)*j)
              fullv4(1,1,nfill+j-1,1)=v4(1,1,i-1,1)+(((v4(1,1,i,1)
     .          -v4(1,1,i-1,1))/dp)*j)
              if (j.ne.dp) then
                fulln(1,1,nfill+j-1,1)=0.0
                fullq(1,1,nfill+j-1,1)=222.0
              else
                fulln(1,1,nfill+j-1,1)=nbin(1,1,i,1)
                fullq(1,1,nfill+j-1,1)=qual(1,1,i,1)
              endif
  96        continue
            nfill=nfill+dp
          else
            do 97 j=1,dp
              fullp(1,1,nfill+j-1,1)=pr(1,1,i-1,1)+j
              fullt(1,1,nfill+j-1,1)=t068(1,1,i-1,1)+(((t068
     .          (1,1,i,1)-t068(1,1,i-1,1))/dp)*j)
              fullc(1,1,nfill+j-1,1)=c0(1,1,i-1,1)+(((c0(1,1,i,1)
     .          -c0(1,1,i-1,1))/dp)*j)
              fullv0(1,1,nfill+j-1,1)=v0(1,1,i-1,1)+(((v0(1,1,i,1
     .          )-v0(1,1,i-1,1))/dp)*j)
              fullv1(1,1,nfill+j-1,1)=v1(1,1,i-1,1)+(((v1(1,1,i,1
     .          )-v1(1,1,i-1,1))/dp)*j)
              fullv2(1,1,nfill+j-1,1)=v2(1,1,i-1,1)+(((v2(1,1,i,1
     .          )-v2(1,1,i-1,1))/dp)*j)
              fullv4(1,1,nfill+j-1,1)=v4(1,1,i-1,1)+(((v4(1,1,i,1
     .          )-v4(1,1,i-1,1))/dp)*j)
              if (j.ne.dp) then
                fulln(1,1,nfill+j-1,1)=0.0
                fullq(1,1,nfill+j-1,1)=666.0
              else
                fulln(1,1,nfill+j-1,1)=nbin(1,1,i,1)
                fullq(1,1,nfill+j-1,1)=qual(1,1,i,1)
              endif
  97        continue
            nfill=nfill+dp
          endif
  98    continue
        nfill=nfill-1

c  Compute potential density array for inversion criteria:

        do 100 i=1,nfill
          sal=sal78(fullc(1,1,i,1)/42.914,fullt(1,1,i,1),
     .      fullp(1,1,i,1),0)
          sv=svan(sal,fullt(1,1,i,1),fullp(1,1,i,1),du)
          tht=theta(sal,fullt(1,1,i,1),fullp(1,1,i,1),0.)
          du=svan(sal,tht,0.,sigth(i))
 100    continue

c  Linearly interpolate raw data over potential density inversions:

        dp=0.0
        do 150 i=2,nfill
          ds = sigth(i)-sigth(i-1)
          if (ds.lt.-0.01) then
            dp=2.0
            do 110 jj=1,50
              if (sigth(i+jj).lt.sigth(i-1)) then
                dp=dp+1.0
              else
                goto 115
              endif
 110        continue
 115        write (6,120) fullp(1,1,i,1),ds,dp
            write (10,120) fullp(1,1,i,1),ds,dp
 120        format (' Sigth inversion at',f11.3,' db of',f11.3,' over',
     .        f5.1,' db.')

            do 130 j=1,dp-1
              fullp(1,1,i-1+j,1)=fullp(1,1,i-1,1)+(((fullp(1,1,i+dp-1,1)
     .          -fullp(1,1,i-1,1))/dp)*j)
              fullt(1,1,i-1+j,1)=fullt(1,1,i-1,1)+(((fullt(1,1,i+dp-1,1)
     .          -fullt(1,1,i-1,1))/dp)*j)
              fullc(1,1,i-1+j,1)=fullc(1,1,i-1,1)+(((fullc(1,1,i+dp-1,1)
     .          -fullc(1,1,i-1,1))/dp)*j)
              fullv0(1,1,i-1+j,1)=fullv0(1,1,i-1,1)+(((fullv0(1,1,i+dp-1
     .          ,1)-fullv0(1,1,i-1,1))/dp)*j)
              fullv1(1,1,i-1+j,1)=fullv1(1,1,i-1,1)+(((fullv1(1,1,i+dp-1
     .          ,1)-fullv1(1,1,i-1,1))/dp)*j)
              fullv2(1,1,i-1+j,1)=fullv2(1,1,i-1,1)+(((fullv2(1,1,i+dp-1
     .          ,1)-fullv2(1,1,i-1,1))/dp)*j)
              fullv4(1,1,i-1+j,1)=fullv4(1,1,i-1,1)+(((fullv4(1,1,i+dp-1
     .          ,1)-fullv4(1,1,i-1,1))/dp)*j)
              if (dp.gt.2) then
                fulln(1,1,i-1+j,1)=0.0
                fullq(1,1,i-1+j,1)=666.0
              else
                fulln(1,1,i-1+j,1)=0.0
                fullq(1,1,i-1+j,1)=222.0
              endif
 130        continue
c           i=i+dp-1
          endif
 150    continue

c  Apply temperature corrections and conductivity slope and offset (and beta);
c  compute ITS-90, final salinity, theta, sigma-t, sigma-th, and dynamic height:

	do 165 kk=1,100
	  if (INT(csta(kk)).eq.ista) slopes=slope(kk)
 165	continue

	dh(1,1,1,1)=0.0
        svl=0.0

c  Calcos0 results for GP806 secondary cond s/n 1467:

	drift = 0.00078
	bias1 = -0.007508372
	slopes = 1.0005433
	beta1 = 0.0
	offset = -0.005

        do 180 i=1,nfill
	  fullt(1,1,i,1)=fullt(1,1,i,1)+drift-0.0006
c	  if (ista.le.19) then
	    ccal(1,1,i,1)=fullc(1,1,i,1)*slopes+bias1+beta1*fullp(1,1,i,1)
c	  else
c	    ccal(1,1,i,1)=fullc(1,1,i,1)*slopes+bias2+beta2*fullp(1,1,i,1)
c	  endif
          scal(1,1,i,1)=sal78(ccal(1,1,i,1)/42.914,fullt(1,1,i,1),
     .      fullp(1,1,i,1),0)+offset
          sv=svan(scal(1,1,i,1),fullt(1,1,i,1),fullp(1,1,i,1),du)
          th(1,1,i,1)=theta(scal(1,1,i,1),fullt(1,1,i,1),
     .      fullp(1,1,i,1),0.)
          du=svan(scal(1,1,i,1),fullt(1,1,i,1),0.,sig(1,1,i,1))
          du=svan(scal(1,1,i,1),th(1,1,i,1),0.,sth(1,1,i,1))
          if (i.ne.1) dh(1,1,i,1)=dh(1,1,i-1,1)+(fullp(1,1,i,1)-
     .      fullp(1,1,i-1,1))*(sv+svl)/2.*1e-5
          t90(1,1,i,1)=fullt(1,1,i,1)*0.99976
          svl=sv
 180    continue

c  Open output file:

       write (outfile,200) cruise,stacst
 200   format (a5,'c',a4,'.nc')
	call epopen(epfout, ' ', EPCREATE +EPCDF)
	call epmesgdo(epfout, 6, EPSTOP)
	call epsetnext(epfout, outfile)
 
c  Put character attributes into header: 

       call epputattvalc(epfout, EPINSTTYPE,   20, 'TAO SBE 911+ CTD 439
     .')
       call epputattvalc(epfout, EPDATATYPE,    3, 'CTD')
       call epputattvalc(epfout, EPSUBTYPE,     3, 'CTD')
       call epputattvalc(epfout, EPDATAORIGIN, 23, 'NOAA/PMEL/OCRD/MCPHA
     .DEN')
c      call epputattvalc(epfout, EPWATERMASS,   1, 'P')
       call epputattvalc(epfout, EPCOORDSYS,   12, 'GEOGRAPHICAL')
       call epputattvalc(epfout, 'Conventions',13, 'PMEL-EPIC/CTD')

       call epputattvalc(epfout, EPCAST,        4, stacst)
       call epputattvalc(epfout, EPCLOUDAMOUNT, 1, amt)
       call epputattvalc(epfout, EPCLOUDTYPE,   1, cloud)
       call epputattvalc(epfout, EPCRUISE,     10, crs_str)
       call epputattvalc(epfout, EPSEASTATE,    1, state)
       call epputattvalc(epfout, EPVISIBILITY,  1, vis)
       call epputattvalc(epfout, EPWEATHER,     1, weat)

c  Put integer attributes into header:

       call epputattval (epfout, EPBAROMETER,  EPINT, 1, bar)
       call epputattval (epfout, EPBOTTLE,     EPINT, 1, 0)
       call epputattval (epfout, EPWATERDEPTH, EPINT, 1, dep)
       call epputattval (epfout, EPWINDDIR,    EPINT, 1, wdir)
       call epputattval (epfout, EPWINDSPEED,  EPINT, 1, wspd)

c  Put real attributes into header:

       call epputattval (epfout, EPAIRTEMP,   EPREAL, 1, dry)
       call epputattval (epfout, EPWETBULB,   EPREAL, 1, wet)

c  Write to output using subroutine writep:

       do k=1,nfill
         pres(k)=fullp(1,1,k,1)
       enddo
       call writep(20, epfout, fullt, 1, 1, nfill, 1,
     .   lon, lat, pres, itime, nx, ny, nz, nt)
       call writep(28, epfout, t90, 1, 1, nfill, 1,
     .   lon, lat, pres, itime, nx, ny, nz, nt)
       call writep(41, epfout, scal, 1, 1, nfill, 1,
     .   lon, lat, pres, itime, nx, ny, nz, nt)
       call writep(70, epfout, sig, 1, 1, nfill, 1,
     .   lon, lat, pres, itime, nx, ny, nz, nt)
       call writep(10, epfout, dh, 1, 1, nfill, 1,
     .   lon, lat, pres, itime, nx, ny, nz, nt)
       call writep(50, epfout, ccal, 1, 1, nfill, 1,
     .   lon, lat, pres, itime, nx, ny, nz, nt)
       call writep(30, epfout, th, 1, 1, nfill, 1,
     .   lon, lat, pres, itime, nx, ny, nz, nt)
       call writep(71, epfout, sth, 1, 1, nfill, 1,
     .   lon, lat, pres, itime, nx, ny, nz, nt)
       call writep(105, epfout, fulln, 1, 1, nfill, 1,
     .   lon, lat, pres, itime, nx, ny, nz, nt)
       call writep(104, epfout, fullq, 1, 1, nfill, 1,
     .   lon, lat, pres, itime, nx, ny, nz, nt)
       call writep(978, epfout, fullv0, 1, 1, nfill,1,
     .   lon, lat, pres, itime, nx, ny, nz, nt)
       call writep(971, epfout, fullv1, 1, 1, nfill,1,
     .   lon, lat, pres, itime, nx, ny, nz, nt)
       call writep(4012, epfout, fullv2, 1, 1,nfill,1,
     .   lon, lat, pres, itime, nx, ny, nz, nt)
       call writep(977, epfout, fullv4, 1, 1, nfill,1,
     .   lon, lat, pres, itime, nx, ny, nz, nt)


       goto 19
 900   close (3)
       close (10)
       close (7)
c      close (9)
       call epclose(epfout)
       stop
       end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
       subroutine writep(varid, epfout, data, nlon, nlat, ndep, ntim, 
     .                            lon, lat, depth, itime,
     .                               nx, ny, nz, nt)
c
c Given an input variable contained in the 4-D array called data
c   with EPIC key code varid, create a NetCDF file containing that 
c   variable and its axes and named by character variable outfile. 
c   nx, ny, nz, and nt specify the dimensions of data and its axes 
c   as declared in the calling program.
c
c Programmed by Dai McClurg, NOAA/PMEL/OCRD, Aug 1993
c
c ----------------------------------------------------------------------
c Arguments:
c
c I/O	name	description
c
c  I	varid		integer EPIC key code for variable in array data
c  I	outfile		character output file name
c  I	data		real array containing variable
c  I	nlon		integer length of longitude axis
c  I	nlat		integer length of latitude axis
c  I	ndep		integer length of depth axis
c  I	ntim		integer length of time axis
c  I	lon		real longitude axis
c  I	lat		real latitude axis
c  I	depth		real depth axis
c  I	itime		integer EPIC time axis
c  I	nx, ny, nz, nt	integer dimensions of data and axes
c ----------------------------------------------------------------------
c
      implicit none
c
      include '/home/ctd/src/epsystemk.fh'
c
      integer nx, ny, nz, nt
c
      real data(nx,ny,nz,nt), depth(nz), lon(nx), lat(ny)
c
      integer nlon, nlat, ndep, ntim
      integer itime(2,nt), varid
c
      integer epfout, opts, lci(4), uci(4), origin(4), dims(4)
      parameter(opts = EPCREATE + EPCDF)
c
      character axtype*6
c
      integer i, j, k
      real delta
c
c.........................................................................
c
      call epsetindx(nx,ny,nz,nt,dims)
c
      call epsetindx(1,1,1,1,lci)
      call epsetindx(1,1,1,1,origin)
      call epsetindx(nlon,nlat,ndep,ntim,uci)
c
c Longitude
c
      axtype = 'EVEN'
c
      if(nlon .gt. 2) then
        delta = lon(2) - lon(1)
        do i = 3, nlon
          if((lon(i) - lon(i-1)) .ne. delta) then
            axtype(1:6) = 'UNEVEN'
            goto 10
          endif
        enddo
      endif
c
   10 continue
c
      call epsetaxis(epfout, varid, 1, 502, 'lon', ' ', 
     .                    'degree_east', axtype)
c
      call epputaxis(epfout,varid,1,origin,lci,uci,lon)
c
c  Latitude
c
      axtype(1:4) = 'EVEN'
      axtype(5:6) = '  '
c
      if(nlat .gt. 2) then
        delta = lat(2) - lat(1)
        do j = 3, nlat
          if((lat(j) - lat(j-1)) .ne. delta) then
            axtype(1:6) = 'UNEVEN'
            goto 20
          endif
        enddo
      endif
c
   20 continue
c
      call epsetaxis(epfout, varid, 2, 500, 'lat', ' ',
     .                    'degree_north', axtype)
c
      call epputaxis(epfout,varid,2,origin,lci,uci,lat)
c
c Depth
c
      axtype(1:4) = 'EVEN'
      axtype(5:6) = '  '
c
      if(ndep .gt. 2) then
        delta = depth(2) - depth(1)
        do k = 3, ndep
          if((depth(k) - depth(k-1)) .ne. delta) then
            axtype(1:6) = 'UNEVEN'
            goto 30
          endif
        enddo
      endif
c
   30 continue
c
      call epsetaxis(epfout,varid,3,1,'pressure',' ','db',axtype)
      call epputaxis(epfout,varid,3,origin,lci,uci,depth)
c
c Time
c
      call epsettaxis(epfout,varid,'time',' ','Julian',EPEVEN)
      call epputtaxis(epfout,varid,origin,lci,uci,itime)
c
c     write(*,*) ' '
c     write(*,*) 'varid = ',  varid
c     write(*,*) 'dims = ',   dims
c     write(*,*) 'origin = ', origin
c     write(*,*) 'lci = ',    lci
c     write(*,*) 'uci = ',    uci
c
      call epputvar(epfout,varid,origin,lci,uci,data,dims)
c
      return
c
      end