   program proc_oml

!  Program to read ocean temp data from HYCOM data and estimate mixed layer depth 
!     and write it out in intermediate file format
!
!  WW: Feb, 2009
!
!===========================================================================
!
!  To make executable:
!
!   Linux flags for PGI compiler:
!      pgf90 proc_oml.f -L/usr/local/netcdf/lib -lnetcdf \
!      -I/usr/local/netcdf/include  -Mfree -byteswapio -o proc_oml
!
!   Linux flags for ifort compiler:
!      ifort proc_oml.f -L/usr/local/netcdf/lib -lnetcdf \
!      -I/usr/local/netcdf/include  -FR -convert big_endian -o proc_oml
!
!   SGI flags
!      f90 proc_oml.f -L/usr/local/netcdf/lib -lnetcdf \
!      -I/usr/local/netcdf/include  -freeform  -o proc_oml
!
!   IBM flags (NCAR bluefire)
!      xlf proc_oml.f -L/usr/local/netcdf/lib -lnetcdf \
!      -I/usr/local/netcdf/include  -qfree=f90  -o proc_oml
!      xlf proc_oml.f -L/usr/local/netcdf/lib -lnetcdf -I/usr/local/netcdf/include -qfree=f90  -o proc_oml
!
!   Mac flags (with g95 compiler)
!	g95 proc_oml.f -L/usr/local/netcdf/lib -lnetcdf \
!	-I/usr/local/netcdf/include -ffree-form -fendian=big -o proc_oml
!
!===========================================================================
!
!   To run:
!      proc_oml file.nc yyyymmddhh
!
      implicit none
      INCLUDE 'netcdf.inc'
      integer ncid, status
      integer  :: jdim
      parameter (jdim=6)
      integer ishape(jdim)
      integer ishape2(jdim)
      character cval*50
      character name*31
      character (len=31), allocatable, dimension(:)       :: dname, dname2
      integer,            allocatable, dimension(:)       :: dval, dval2
      real,               allocatable, dimension(:,:,:,:) :: data
      real,               allocatable, dimension(:,:,:  ) :: temp
      real,               allocatable, dimension(:,:    ) :: temp1
      real,               allocatable, dimension(:,:)     :: data2d
      real,               allocatable, dimension(:)       :: data1d
      double precision,   allocatable, dimension(:,:)     :: ddata
      double precision,   allocatable, dimension(:)       :: ddata1
      integer,            allocatable, dimension(:)       :: idata
      real,               allocatable, dimension(:,:)     :: mld, h, tmean, tot_depth
      real,               allocatable, dimension(:,:)     :: mld_out, tm_out, t0_out
      real,               allocatable, dimension(:)       :: depth
      real,               allocatable, dimension(:,:)     :: xlat, xlon
      real,               allocatable, dimension(:)       :: xlat_out, xlon_out
      character,          allocatable, dimension(:,:)     :: text
      character omit(10)*80
      integer             :: start_dims(2)
      integer             :: dims_in(4), dims_out(4)
      integer             :: idm, ndims, nvars, natt, ngatts, nunlimdimid
      integer             :: i, ii, i1, j, j1, ix, iy, it, iz, k, iy_in, kbot, ix_out, iy_out
      integer             :: ilen, itype, ival, na, ndif
      integer             :: mcid
      integer             :: varid
      real                :: dx, rval, xlat_tmp, xlat_dif, hdif, gamma, delta_lat, delta_last, delta_lon
      real                :: f1, f2, msgval, msgval2
      character (len=80)  :: input_file
      character (len=10)  :: option
      logical             :: debug=.TRUE.
      logical             :: bit64=.FALSE.
      character (len=10)  :: input_date
      character (len=19)  :: date
      integer :: year,month,day,hour

      integer, external   :: iargc
      integer             :: numarg

! data for intermediate file
      character(len=24) :: hdate
      character(len=9)  :: field
      character(len=25) :: units
      character(len=46) :: desc
      character(len=32) :: map_source
      character(len=8)  :: startloc
      integer :: nx, ny
      real :: xlvl, xfcst, earth_radius
      real :: startlat, startlon
      real :: deltalat, deltalon
      character(len=80) :: flnm, flnm2
      character(len=2)  :: verbose = " "
      logical           :: is_wind_grid_rel

      integer :: ifv, iproj, ilev, ierr, iferr, iop
      real    :: dy, truelat1, truelat2, lon2

      msgval  = 10.
      msgval2 = 200.
      numarg = iargc()

      if (numarg .lt. 1) call help_info

      call getarg(1,input_file)
      call getarg(2,input_date)

      print *, 'input_date ', input_date
      read (input_date(1:4), '(I4)') year
      read (input_date(5:6), '(I2)') month
      read (input_date(7:8), '(I2)') day
      read (input_date(9:10), '(I2)') hour

      write (date, '(I4,"-",I2.2,"-",I2.2,"_",i2.2,":00:00")') year,month,day,hour
      print *,'Date : ', date

! input and output files:

      status = nf_open(input_file, 0, ncid)
      if (status .ne. nf_noerr) call handle_err(status)
!     if (bit64) then
!       status = nf_create(output_file, NF_64BIT_OFFSET, mcid)
!     else
!       status = nf_create(output_file, 0, mcid)
!     endif
!     if (status .ne. nf_noerr) call handle_err(status)

! GET BASIC INFORMTION ABOUT THE FILE
! most important 
!   ndims:  number of dimensions
!   nvars:  number of variables
!   ngatts: number of global attributes

      status = nf_inq(ncid, ndims, nvars, ngatts, nunlimdimid)
      if (status .ne. nf_noerr) call handle_err(status)
      IF (debug) THEN
        write(6,*) 'ndims = ',ndims,' nvars = ',nvars,' ngatts = ',ngatts, &
                   ' nunlimdimid =',nunlimdimid
      ENDIF

! ALLOCATE SOME VARIABLES
      allocate  (dval(ndims))
      allocate (dname(ndims))
      allocate  (dval2(ndims))
      allocate (dname2(ndims))

! READ ALL DIMS FROM INPUT FILE AND CREATE DIMS FOR OUTPUT FILE
      j = 1
      do i = 1, ndims
        status = nf_inq_dim(ncid, i, dname(i), dval(i))
        if (debug) write (6,*) dname(i), trim(dname(i))
        if (dname(i) .eq. 'X' .or. dname(i) .eq. 'Y') then
           dval2(j) = dval(i)
           dname2(j) = dname(i)
!       if ( dname(i) == "Time" ) then
!         status = nf_def_dim(mcid, dname(i), NF_UNLIMITED, i)
!       else
!         status = nf_def_dim(mcid, dname(i), dval2(i), i)
!       end if
        IF (debug) THEN
          write(6,'(" i = ",i2," : ",A," input = ",i4," (output = ",i4,")")') &
                i,dname(i),dval(i),dval2(j)
        ENDIF
        j = j+1
        endif
      enddo

      IF (debug) THEN
        write(6,*)
        write(6,*) "FILE variables:"
      ENDIF
! some preliminary stuff for new netcdf file
      j1 = 1
      do i = 1, nvars
        if (i.ne.6) cycle
        status = nf_inq_var(ncid, i, cval, itype, idm, ishape, natt)
        if (idm .gt. jdim) stop
        status = nf_def_var(mcid, cval, itype, idm, ishape, j)
        print *, 'j = ', j
        do na = 1, natt
          status = nf_inq_attname(ncid, i, na, name)
          status = nf_copy_att(ncid, i, name, mcid, j)
        enddo
        print *, 'j1 = ', j1
        j1 = j1+1
      enddo
      status = nf_enddef(mcid)

! ########## LOOP THROUGH THE DATA 
      do i = 1, nvars
        status = nf_inq_var(ncid, i, cval, itype, idm, ishape, natt)
        IF (debug) THEN
          write(6,*) 
	  write(6,*) 'VARIABLE: ',trim(cval), '  DIM: ', idm
        ENDIF

! GET THE DIMS FOR INPUT AND OUTPUT FROM THE SHAPE
        dims_in = 1
        do ii = 1,idm
           dims_in(ii)  = dval(ishape(ii))
        enddo
        IF (debug) THEN
          write(6,*) '   DIMS IN: ', dims_in
        ENDIF

! ALLOCATE THE INPUT AND OUTPUT ARRAYS
! READ THE DATA FROM INPUT FILE

	IF     (itype .eq. NF_CHAR ) THEN          ! character
          allocate (text(dims_in(1), dims_in(2)))
	  status = nf_get_var_text(ncid, i, text)
          IF (debug) write(6,*) '   SAMPLE VALUE = ',text(:,1)
          deallocate (text)

	ELSEIF (itype .eq. NF_INT ) THEN          ! integer (X and Y)
          allocate (idata(dims_in(1)))
	  status = nf_get_var_int(ncid, i, idata)
          IF (debug) write(6,*) '   SAMPLE VALUE = ',idata(100)
          IF (trim(cval) .eq. 'X') THEN
              ix = dims_in(1)
          ELSE IF (trim(cval) .eq. 'Y') THEN
              iy = dims_in(1)
          ENDIF
          deallocate (idata)

	ELSEIF (itype .eq. NF_FLOAT) THEN          ! real
          IF (idm .eq. 1) THEN                     ! depth
             allocate (data1d(dims_in(1)))
             status = nf_get_var_real(ncid, i, data1d)
             iz = dims_in(1)
             IF (debug) write(6,*) '   SAMPLE VALUE = ',data1d(1:10), iz
             IF (trim(cval) .eq. 'Depth') THEN
                allocate (depth(dims_in(1)))
                depth = data1d
             ENDIF
             deallocate (data1d)
          ENDIF
          IF (idm .eq. 2) THEN                     ! lat and lon
             allocate (data2d(dims_in(1), dims_in(2)))
             status = nf_get_var_real(ncid, i, data2d)
             IF (debug) write(6,*) '   SAMPLE VALUE = ',data2d(1:12,1:12)
             IF (trim(cval) .eq. 'Latitude') THEN
                allocate (xlat(dims_in(1), dims_in(2)))
                xlat = data2d
             ELSE IF (trim(cval) .eq. 'Longitude') THEN
                allocate (xlon(dims_in(1), dims_in(2)))
                xlon = data2d
             ENDIF
             deallocate (data2d)
          ENDIF
          IF (idm .eq. 4) THEN                     ! temp (C)
             allocate (data(dims_in(1), dims_in(2), dims_in(3), dims_in(4)))
             status = nf_get_var_real(ncid, i, data)
             IF (debug) write(6,*) '   SAMPLE VALUE = ',data(1:12,1:12,1,1)
             print *, dims_in(1), dims_in(2), dims_in(3)
!            allocate (temp(dims_in(1), dims_in(2), dims_in(3), dims_in(4)))
! allocate 10 layers only
             allocate (temp(dims_in(1), dims_in(2), 10) )
             temp = data(:,:,1:10,1)
             IF (debug) write(6,*) '   SAMPLE VALUE = ',temp(1000,1000,1)
             deallocate (data)
          ENDIF

        ELSEIF (itype .eq. NF_DOUBLE) THEN          ! double
          IF (debug) write(6,*) '   READING DBL PRECISION DATA'
          IF (idm .eq. 1) THEN
            allocate (ddata1(dims_in(1)))
            status = nf_get_var_double(ncid, i, ddata1)
            IF (debug) write(6,*) '   SAMPLE VALUE = ',ddata1(1)
            deallocate (ddata1)
          ELSE IF (idm .eq. 2) THEN
            allocate (ddata(dims_in(1), dims_in(2)))
            status = nf_get_var_double(ncid, i, ddata)
            IF (debug) write(6,*) '   SAMPLE VALUE = ',ddata(1,1)
            deallocate (ddata)
          ENDIF

        ENDIF

      ENDDO     ! END OF VARIABLE LOOP

      print *, 'Have read netcdf data, xdim, ydim, zdim are ', ix, iy, iz
      print *, 'lat/lon(1,1), lat/lon(4500,3298) ',xlat(1,1),xlon(1,1),xlat(ix,iy),xlon(ix,iy)

! prepare data

! Allocate spaces
      allocate ( h(ix, iy) )
      allocate ( mld(ix, iy) )
      allocate ( tmean(ix, iy) )
      allocate ( tot_depth(ix, iy) )

! 1. compute ocean heat content in top 100 m -> top 7 layers: 0,10,20,30,50,75,100
      print *, 'Compute OHC'
      kbot = 7
      do i = 1, ix
      do j = 1, iy
         h(i,j) = 0.
         do k = 2, kbot
         if ( temp(i,j,kbot) .lt. 100. ) then
         h(i,j)=h(i,j)+(273.15+(0.5*(temp(i,j,k-1)+temp(i,j,k)))) *(depth(k)-depth(k-1))
         endif
         enddo
      enddo
      enddo

!     do j = 1, iy
!     do i = 1, ix
!     if ((xlat(i,j).gt.24. .and. xlat(i,j).lt.27.).and.(xlon(i,j).gt.267..and.xlon(i,j).lt.270.))then
!     if (i.gt.2245 .and. i.lt.3380 .and. j.gt.2170 .and. j.lt.2180) then
!        print *, 'check ', i,j,xlat(i,j),xlon(i,j),temp(i,j,1),h(i,j)
!     endif
!     enddo
!     enddo

! 2. assume gamma below 100 m is about 0.2, compute estimated mld

      print *, 'Compute estimated MLD'
      gamma = 0.20
      ndif = 0
      do i = 1, ix
      do j = 1, iy
         mld(i,j) = msgval
         if ( temp(i,j,kbot) .lt. 100. ) then
             hdif = (temp(i,j,1)+273.15)*depth(kbot)-h(i,j)
             if (hdif.le.0.) then
!               if(hdif.lt.-30.)print *, 'negative hdif, reset to 0 ', i,j,hdif
                hdif=0.
                ndif=ndif+1
             endif
             mld(i,j)=depth(kbot)-sqrt(2.*hdif/gamma)
      if(i.eq.1324.and.(j.eq.1190.or.j.eq.1192))print*,hdif,mld(i,j),depth(kbot),h(i,j),temp(i,j,kbot)
             if(mld(i,j).lt.0.) then
      if(xlat(i,j).gt.23..and.xlat(i,j).lt.30..and.xlon(i,j).gt.265..and.xlon(i,j).lt.320.) &
                print *, 'negative mld ', i,j,mld(i,j)
                mld(i,j) = 10.
             endif
!     if((i.gt.2415.and.i.lt.2420).and.j.eq.1840) then
!        print *, 'check2 ',hdif, sqrt(2.*hdif/gamma), depth(kbot)
!     endif
         else if ( temp(i,j,1) .lt. 100. ) then
             mld(i,j) = 10.
         endif
!     if ((xlat(i,j).gt.24. .and. xlat(i,j).lt.27.).and.(xlon(i,j).gt.267..and.xlon(i,j).lt.270.))then
!        print *, 'check ', i,j,h(i,j),xlat(i,j),xlon(i,j),mld(i,j)
!     endif
      enddo
      enddo

      print *, 'check for bad values ', ndif, ' points have negative hdif'
      print *, (mld(1324,j),j=1190,1197)

! 2b. compute mean temp in the ocean
!     kbot = 10 corresponding top 10 levels and 200 m depth

      print *, 'Compute mean temp'
      kbot = 10
      ndif = 0
      do i = 1, ix
      do j = 1, iy
         tmean(i,j) = msgval2
         do k = 2, kbot
         if ( temp(i,j,k) .lt. 100. ) then
           if ( k.eq.2) tmean(i,j) = 0.
           tmean(i,j)=tmean(i,j)+(0.5*(temp(i,j,k-1)+temp(i,j,k))) *(depth(k)-depth(k-1))
           tot_depth(i,j)=depth(k)
      if ((xlat(i,j).gt.29. .and. xlat(i,j).lt.29.3).and.(xlon(i,j).gt.271.6.and.xlon(i,j).lt.271.9))then
        print *, 'tmean check ', i,j,tmean(i,j), tot_depth(i,j)
      end if
         endif
         enddo
         if (tmean(i,j).ne.msgval2) tmean(i,j)=tmean(i,j)/tot_depth(i,j)+273.15
      enddo
      enddo

! 2a. now we are mostly done with temp array (which is huge), let's try to get rid of it. 
!     but we do need temp from the first level

      allocate ( temp1(ix, iy) )
      temp1 = temp(:,:,1)
      deallocate ( temp )

! 3. check lat / long grid

      delta_last = 0.032
!     do j = 2, iy-1
!        delta_lat = xlat(1110,j)-xlat(1110,j-1)
!        if (j.eq.2) delta_last = delta_lat
!     if(abs(delta_lat).gt.abs(delta_last+0.00001))print *,'J,delta_lat,delta_last ', j,delta_lat,delta_last,xlat(1110,j)
!        delta_last = delta_lat
!     enddo
!     print *, 'delta_lat ', delta_lat

!     delta_last = 0.08
!     do i = 2, ix-1
!        delta_lon = xlon(i,2300)-xlon(i-1,2300)
!        if (i.eq.2) delta_last = delta_lon
!     if(abs(delta_lon).gt.abs(delta_last+0.00001))print *,'I,delta_lon,delta_last ',i,delta_lon,delta_last,xlon(i,2300)
!        delta_last = delta_lon
!     enddo

!     print *, 'delta_lat, delta_lon ', delta_lat, delta_lon, xlat(1,1), xlon(1,1)
      print *
      print *, 'STOP HERE'
      print *

      delta_lat = 0.08
      delta_lon = 0.08

! 4. interpolate mld data to regular lat/log grid. Note that the data grid in longitude
!       direction is already regular, and the finest resolution in latitude direction is 0.08 deg
!
!     First prepare the output lat/lon grid
!     xlat_out=(-30N, 55N), xlon_out=(200,340) at a resolution of 5 min or 1/12 degrees
!     9/8/09: extend the east boundary to 350 degree 
!     5/5/11: extend the east boundary 180 E to 10 E, and from 20 S to 65 N

!     ix_out = int((350-220)/delta_lon)+1
!     iy_out = int((55+30)/delta_lat)+1
      ix_out = int((380-180)/delta_lon)+1
      iy_out = int((65+20)/delta_lat)+1

      print *, 'Dimensions of output grid: ', ix_out, iy_out

      allocate (xlat_out(iy_out))
      allocate (xlon_out(ix_out))

      xlat_out(1) = -20.
      do j = 2, iy_out
         xlat_out(j) = xlat_out(j-1)+delta_lat
      end do
      print *, xlat_out(800), xlat_out(iy)

      xlon_out(1) = 180.
      do i = 2, ix_out
         xlon_out(i) = xlon_out(i-1)+delta_lon
      end do
      print *, xlon_out(800), xlon_out(iy)
      print *, xlon(1324,2000)

      allocate ( mld_out(ix_out, iy_out) )
      allocate ( tm_out(ix_out, iy_out) )
      allocate ( t0_out(ix_out, iy_out) )
      mld_out = msgval
      tm_out  = msgval2
      t0_out  = msgval2

!     do i = 1, ix_out
      do i1 = 1324, 1324+ix_out-1
         i = i1-1323
        do j = 1, iy_out
           xlat_tmp = -1.
           do j1 = 1, iy
              xlat_tmp = xlat(i1,j1)-xlat_out(j)
              if (xlat_tmp .le. 0.) cycle
              if(mld(i1,j1-1).ne.msgval.and.mld(i1,j1).ne.msgval)then
              xlat_dif = xlat(i1,j1)-xlat(i1,j1-1)
              f1 = (xlat(i1,j1)-xlat_out(j))/xlat_dif
              f2 = (xlat_out(j)-xlat(i1,j1-1))/xlat_dif
              mld_out(i,j) = f1*mld(i1,j1-1) + f2*mld(i1,j1)
!             tm_out(i,j) = f1*tmean(i1,j1-1)+ f2*tmean(i1,j1)
              t0_out(i,j) = f1*temp1(i1,j1-1)+ f2*temp1(i1,j1)
              t0_out(i,j) = t0_out(i,j)+273.15
              else if (mld(i1,j1-1).ne.msgval.and.mld(i1,j1).eq.msgval)then
                mld_out(i,j) = mld(i1,j1-1)
!               tm_out(i,j) = tmean(i1,j1-1)
                t0_out(i,j) = temp1(i1,j1-1) + 273.15
              else if (mld(i1,j1-1).eq.msgval.and.mld(i1,j1).ne.msgval)then
                mld_out(i,j) = mld(i1,j1)
!               tm_out(i,j) = tmean(i1,j1)
                t0_out(i,j) = temp1(i1,j1) + 273.15
              else
                mld_out(i,j) = msgval
!               tm_out(i,j) = msgval2
                t0_out(i,j) = msgval2
              endif
! cannot do proper interpolation beyond 45 N
!             if (xlat_out(j).ge.47.) then
!                mld_out(i,j) = 70.
!                mld_out(i,j) = mld_out(i,j-1)
!                t0_out(i,j) = 287.
!             endif
!             print *, i1,j1,j,xlat(i1,j1),xlat_out(j),xlat(i1,j1-1),mld(i1,j1),mld_out(i,j),mld(i1,j1-1)
              exit
           enddo
        enddo
      enddo

      do i1 = 1324, 1324+ix_out-1
         i = i1-1323
       do j = 1, iy_out
           xlat_tmp = -1.
           do j1 = 1, iy
              xlat_tmp = xlat(i1,j1)-xlat_out(j)
              if (xlat_tmp .le. 0.) cycle
              if(tmean(i1,j1-1).ne.msgval2.and.tmean(i1,j1).ne.msgval2)then
              xlat_dif = xlat(i1,j1)-xlat(i1,j1-1)
              f1 = (xlat(i1,j1)-xlat_out(j))/xlat_dif
              f2 = (xlat_out(j)-xlat(i1,j1-1))/xlat_dif
              tm_out(i,j) = f1*tmean(i1,j1-1)+ f2*tmean(i1,j1)
              else if (tmean(i1,j1-1).ne.msgval2.and.tmean(i1,j1).eq.msgval2)then
                tm_out(i,j) = tmean(i1,j1-1)
              else if (tmean(i1,j1-1).eq.msgval2.and.tmean(i1,j1).ne.msgval2)then
                tm_out(i,j) = tmean(i1,j1)
              else
                tm_out(i,j) = msgval2
              endif
! it is ok to do interpolation beyond 45 N for Atlantic - lat variation is reasonble
!             if (xlat_out(j).ge.47.) then
!                tm_out(i,j) = 287.
!                tm_out(i,j) = tm_out(i,j-1)
!             endif
!             print *, i1,j1,j,xlat(i1,j1),xlat_out(j),xlat(i1,j1-1),tmean(i1,j1),tm_out(i,j),tmean(i1,j1-1)
              exit
           enddo
       enddo
      enddo

!     do i = 1, ix_out
!     do j = 1, iy_out
!        if (mld_out(i,j) .eq. msgval) print *, i,j,mld_out(i,j)
!        if (mld_out(i,j) .gt. 1.e+10 .and. xlat_out(j) .lt. 50.) print *, i,j,mld_out(i,j)
!        if (mld_out(i,j) .eq. 0.) mld_out(i,j) = 1.e20
!     enddo
!     enddo

! 5. Write the fields out in intermediate format (T0 is not used in the model at the moment)

      open ( 15, file='MLD', form='unformatted' )
      ifv = 5
      hdate      = '                        '
      hdate(1:19)= date
      xfcst      = 0.
      map_source = 'HYCOM                          '
      field      = 'H0ML     '
      units      = 'm                        '
      desc       = 'Ocean Mixed Layer Depth Estimate              '
      xlvl       = 200100.
! for cylindrical or regular lat/lon grid
      iproj      = 0
      startloc   = 'SWCORNER'
      startlat   = xlat_out(1)
!     startlon   = xlon_out(1)
      startlon   = -180.
      deltalat   = delta_lat
      deltalon   = delta_lon
      earth_radius = 6370.
      is_wind_grid_rel = .false.

      write (15) ifv
      write (15) hdate, xfcst, map_source, field, units, desc, xlvl, ix_out, iy_out, iproj
      write (15) startloc, startlat, startlon, deltalat, deltalon, earth_radius
      write (15) is_wind_grid_rel
      write (15) mld_out

      xfcst      = 0.
      map_source = 'HYCOM                          '
      field      = 'TMOML    '
      units      = 'K                        '
      desc       = 'Mean Ocean Temperature in Top 200 m           '
      xlvl       = 200100.
! for cylindrical or regular lat/lon grid
      iproj      = 0
      startloc   = 'SWCORNER'
      startlat   = xlat_out(1)
!     startlon   = xlon_out(1)
      startlon   = -180.
      deltalat   = delta_lat
      deltalon   = delta_lon
      earth_radius = 6370.
      is_wind_grid_rel = .false.

      write (15) ifv
      write (15) hdate, xfcst, map_source, field, units, desc, xlvl, ix_out, iy_out, iproj
      write (15) startloc, startlat, startlon, deltalat, deltalon, earth_radius
      write (15) is_wind_grid_rel
      write (15) tm_out

      xfcst      = 0.
      map_source = 'HYCOM                          '
      field      = 'T0       '
      units      = 'K                        '
      desc       = 'Top Level Ocean Temp                          '
      xlvl       = 200100.
! for cylindrical or regular lat/lon grid
      iproj      = 0
      startloc   = 'SWCORNER'
      startlat   = xlat_out(1)
!     startlon   = xlon_out(1)
      startlon   = -180.
      deltalat   = delta_lat
      deltalon   = delta_lon
      earth_radius = 6370.
      is_wind_grid_rel = .false.

      write (15) ifv
      write (15) hdate, xfcst, map_source, field, units, desc, xlvl, ix_out, iy_out, iproj
      write (15) startloc, startlat, startlon, deltalat, deltalon, earth_radius
      write (15) is_wind_grid_rel
      write (15) t0_out

! write a new netcdf file for 3 variables: mld, xlat and xlon - skip

      go to 102
      start_dims = 1
      dims_out = 0
      j = 1
      do i = 1, 2
         status = nf_inq_dim(mcid, i, dname2(j), dval2(j))
!        status = nf_inq_var(ncid, i, cval, itype, idm, ishape, natt)
!        do ii = 1,idm
!          dims_in(ii)  = dval(ishape(ii))
!          dims_out(ii) = dims_in(ii)
!        enddo
         if ( i .eq. 1 ) then
            dname2(i) = 'Delta Lat'
            print *, 'j, dims_out ', j, dims_out,trim(cval)
         status = nf_put_vara_real (mcid, j, start_dims, dims_out, xlat_out)
            j = j+1
         elseif (trim(cval) .eq. 'lon') then
            print *, 'j, dims_in  ', j, dims_out,trim(cval)
            dims_out(1) = ix
            dims_out(2) = 0
            dims_out(3) = 0
            print *, 'j, dims_out ', j, dims_out,trim(cval)
         status = nf_put_vara_real (mcid, j, start_dims, dims_out, xlon_out)
            j = j+1
         elseif (trim(cval) .eq. 'mld') then
            print *, 'j, dims_in  ', j, dims_out,trim(cval)
            dims_out(2) = iy
            dims_out(3) = 0
            print *, 'j, dims_out ', j, dims_out,trim(cval)
         status = nf_put_vara_real (mcid, j, start_dims, dims_out, mld)
            j = j+1
         endif
      enddo
      status = nf_close(mcid)
102   continue

      write(6,*) 
      write(6,*) " Success"      
      write(6,*) "#########################################"

      end program proc_oml

!---------------------------------------------------------------------
      subroutine handle_err(status)
      integer status
      write(6,*) 'Error number ',status
      stop
      end subroutine

!------------------------------------------------------------------------------
  subroutine help_info

  print*," "
  print*," readnc  netcdf_input "
! print*," readnc  netcdf_input netcdf_output"
  print*," "
  STOP
  end subroutine help_info
