#include "options.inc"

c
c=======================================================================
c     convert binary output of averaged quantities 
c     to NetCDF format
c     real variables are 8 or 4 bytes long set by
c=======================================================================
c
c#define real_8
c
c     use less memory 
#define low_memory
c
c    unfold domain
#define unfold

      module kind_mod
#ifdef real_8
      integer, parameter::  r_knd  = selected_real_kind(13)
#else
      integer, parameter::  r_knd  = selected_real_kind(6)
#endif
      end module kind_mod

      module cv_ave_module
      use kind_mod
      integer :: imt,jmt,km,nt,imtf

      real (kind=r_knd),allocatable,dimension(:,:,:,:):: mean_u,mean_t
      real (kind=r_knd),allocatable,dimension(:,:)    :: mean_psi
      real (kind=r_knd), allocatable, dimension(:,:,:):: mean_stf
      real (kind=r_knd), allocatable, dimension(:,:,:):: mean_smf
      real (kind=r_knd), allocatable, dimension(:) :: xt,yt,xu,yu,zt,zw
      real (kind=r_knd), allocatable, dimension(:) :: csu
      real (kind=r_knd), allocatable, dimension(:,:,:):: mean_ka
      real (kind=r_knd), allocatable, dimension(:,:,:,:):: mean_ugm
      real (kind=r_knd), allocatable, dimension(:,:):: mean_pco2
      real (kind=r_knd), allocatable, dimension(:,:):: mean_dco2
      real (kind=r_knd), allocatable, dimension(:,:):: mean_o2
      real (kind=r_knd) :: mean_pco2_air
      logical :: found_kappa = .false.
      logical :: found_gm    = .false.
      logical :: found_pco2  = .false.
      logical :: found_o2    = .false.

      integer :: counter
      integer, allocatable, dimension(:,:)     :: kmt,kmu
      logical :: enable_monthly_averages = .false.
      logical :: enable_annual_averages = .false.
      logical :: enable_daily_averages = .false.
      logical :: enable_seasonal_averages = .false.
      logical :: enable_freesurf = .false.
      real, parameter :: pi       = 
     &        3.14159265358979323846264338327950588
      end module cv_ave_module



      program cv_ave
      use cv_ave_module
      implicit none
      character (len=80) :: infile,outfile
      integer :: ierr,i
      character (len=9) :: c9
      character (len=1) :: c1
      character (len=2) :: c2
      character (len=4) :: c4
      integer :: year,month,day,domain

      print*,' Converting averaged quantities written by SPFLAME ',
     &       ' to NetCDF file '

      call get_commandline_arg(1,infile,ierr)
      if (ierr/=0) goto 99
      print*,' reading from file ',infile(1:len_trim(infile))

      read(infile,'(a9,i3,a2,i4,a1,i2,a1,i2,a4)') 
     &            c9,domain,c2,year,c1,month,c1,day,c4
      write(outfile,'("averages_",i3,"_y",i4,"m",i2,"d",i2,".cdf")') 
     &            domain,year,month,day
      do i=1,len_trim(outfile);
        if (outfile(i:i)==' ')outfile(i:i)='0';
      enddo
      print*,' sub domain = ',domain
      print*,' year/month/day  = ',year,month,day
      print*,' writing to file ',outfile
#ifdef unfold
      print*,' unfolding file using kmt/grid.ufl'
#endif
      call read_averages(infile)
      call write_averages(infile,outfile,year,month,day,domain)
      stop

 99   print*,' program needs command line parameter'
      print*,' 1. input binary file name (written by spflame)'
      end program cv_ave


      subroutine read_averages(infile)
      use cv_ave_module
      implicit none
      character (len=*) :: infile
      character (len=3) :: ssh
      integer :: i,j,k,io=10,n,iok=20,iog=30
      character (len=80) :: ident
      logical :: ok
      real (kind=r_knd) :: xx

      open(io,file=infile,form='unformatted')
      read(io) imtf,jmt,km,nt
      imt=imtf

#ifdef unfold
      print*,' reading from file kmt.ufl and grid.ufl'
      open(iok,file='kmt.ufl',form='unformatted')
      read (iok)
      read (iok) 
      read (iok) imt
      print*,' folded imt = ',imtf,' found in kmt.ufl imt = ',imt
      open(iog,file='grid.ufl',form='unformatted')
      read (iog)
      read (iog) imt
      print*,' folded imt = ',imtf,' found in grid.ufl imt = ',imt
#endif


      print*,' imt=',imt,' jmt=',jmt,' km=',km,' nt=',nt
#ifdef low_memory
      allocate( mean_t(imt,1,jmt,nt) )
      allocate( mean_u(imt,1,jmt,2) )
#else
      allocate( mean_t(imt,km,jmt,nt) )
      allocate( mean_u(imt,km,jmt,2) )
#endif
      allocate( mean_psi(imt,jmt) )
      allocate( mean_stf(imt,jmt,nt) )
      allocate( mean_smf(imt,jmt,2) )
      allocate( xt(imt),yt(jmt),xu(imt),yu(jmt),zt(km),zw(km),csu(jmt))
      allocate( kmt(imt,jmt),kmu(imt,jmt) ); kmu=0

      read(io) xt(1:imtf),yt,xu(1:imtf),yu,zt,zw
#ifdef unfold
c   overwrite with full grid
      read (iog) (xx,i=1,(imt+jmt)*2+km*2+1),xt, xu
      close(iog)
#endif
      read(io) counter
      print*,' counter for averages=',counter
      read(io) enable_daily_averages
      read(io) enable_monthly_averages
      read(io) enable_seasonal_averages
      read(io) enable_annual_averages
      read(io) kmt(1:imtf,:)

#ifdef unfold
c   overwrite with full topography
      read (iok) 
      read (iok) kmt
      close(iok)
#endif

      read(io) mean_psi(1:imtf,:),ssh
#ifdef unfold
      mean_psi(imtf-1:imt,:)= mean_psi(1:imt-imtf+2,:)
#endif
      if (ssh == 'ssh') then
        enable_freesurf=.true.
        print*,'found sea surface height in file'
      else
        enable_freesurf=.false.
        print*,'found streamfunction in file'
      endif
      do n=1,2
       read(io) mean_smf(1:imtf,:,n)
#ifdef unfold
       mean_smf(imtf-1:imt,:,n)= mean_smf(1:imt-imtf+2,:,n)
#endif
       do k=1,km
#ifdef low_memory
         read(io) 
#else
         read(io) mean_u(1:imtf,k,:,n)
#ifdef unfold
         mean_u(imtf-1:imt,k,:,n)= mean_u(1:imt-imtf+2,k,:,n)
#endif
#endif
       enddo
      enddo
      do n=1,nt
       read(io) mean_stf(1:imtf,:,n)
#ifdef unfold
       mean_stf(imtf-1:imt,:,n)= mean_stf(1:imt-imtf+2,:,n)
#endif
       do k=1,km
#ifdef low_memory
         read(io) 
#else
         read(io) mean_t(1:imtf,k,:,n)
#ifdef unfold
         mean_t(imtf-1:imt,k,:,n)= mean_t(1:imt-imtf+2,k,:,n)
#endif
#endif
       enddo
      enddo

 10    ident=''; ok = .true.
       read(io,end=20,err=50) ident
       goto 30
 20    ok = .false.
 30    continue

       if (ident(1:8) == 'diff_cbt') then
        print*, ' found vertical diffusivity in file'
        found_kappa = .true.
        allocate(mean_ka(imt,km,jmt) )
        do k=1,km
          read(io) mean_ka(1:imtf,k,:)
#ifdef unfold
          mean_ka(imtf-1:imt,k,:)= mean_ka(1:imt-imtf+2,k,:)
#endif
        enddo
       elseif (ident(1:4) == 'u_gm') then
        print*, ' found GM bolus velocities in file'
        found_gm = .true.
        allocate(mean_ugm(imt,km,jmt,3) )
        do k=1,km
          read(io) mean_ugm(1:imtf,k,:,1)
          read(io) mean_ugm(1:imtf,k,:,2)
          read(io) mean_ugm(1:imtf,k,:,3)
#ifdef unfold
          mean_ugm(imtf-1:imt,k,:,:)= mean_ugm(1:imt-imtf+2,k,:,:)
#endif
        enddo
       elseif (ident(1:10) == 'pco2_water' ) then
        found_pco2=.true.
        print*,' found sea surface CO2 in file'
        allocate(mean_pco2(imt,jmt) )
        allocate(mean_dco2(imt,jmt) )
        read(io) mean_pco2(1:imtf,:)
        read(io) mean_dco2(1:imtf,:)
#ifdef unfold
          mean_pco2(imtf-1:imt,:)= mean_pco2(1:imt-imtf+2,:)
          mean_dco2(imtf-1:imt,:)= mean_dco2(1:imt-imtf+2,:)
#endif
        read(io) mean_pco2_air
       elseif (ident(1:6) == 'o2_sat' ) then
        found_o2=.true.
        print*,' found sea surface O2 in file'
        allocate(mean_o2(imt,jmt) )
        read(io) mean_o2(1:imtf,:)
#ifdef unfold
          mean_o2(imtf-1:imt,:)= mean_o2(1:imt-imtf+2,:)
#endif
       else
         print*,' cannot read ',ident,' in file '
       endif
       if (ok) goto 10  ! read next item til end of file

      close(io)
c
c   construct other variables
c
      do j=1,jmt
       csu(j)   = cos(yu(j)/180*pi)
      enddo
      do j=1,jmt-1
        do i=1,imt-1
           kmu(i,j) = min (kmt(i,j), kmt(i+1,j), 
     &                     kmt(i,j+1), kmt(i+1,j+1))
        enddo
      enddo
      return
 50    print*,' unexpected end of file '
      end subroutine read_averages


      subroutine write_averages(infile,outfile,year,month,day,domain)
      use cv_ave_module
      implicit none
      character (len=80) :: infile,outfile,name,unit
      integer :: year,month,day,domain
#include "netcdf.inc"
      integer ncid,iret,i,j,k,dims(4),start(4),count(4)
      real (kind=r_knd) :: spval=-9.9e12, tt

      integer :: lon_tdim,lon_tid
      integer :: lon_udim,lon_uid
      integer :: lat_tdim,lat_tid
      integer :: lat_udim,lat_uid
      integer :: depth_tdim,depth_tid
      integer :: depth_wdim,depth_wid
      integer :: itimedim,itimeid
      integer :: uid,vid,trid(nt),n,len,diff_cbtid
      integer :: psiid,tauxid, tauyid,strid(nt)
      integer :: ugm_id,vgm_id,wgm_id,pco2_id,o2_id,dco2_id,airco2_id,kk
      integer :: 
     &days_in_month(12) = (/31,28,31,30,31,30,31,31,30,31,30,31/)

      tt=(year-1900)*365 + day
      do n=1,month-1
       tt=tt+days_in_month(n)
      enddo

      if (enable_daily_averages) then
        tt=tt-0.5
      elseif (enable_monthly_averages) then
        tt=tt-days_in_month(month)/2.
      elseif (enable_seasonal_averages) then
        tt=tt-(days_in_month(month)    ! month should be always ge 3
     &        +days_in_month(month-1)
     &        +days_in_month(month-2))/2.
      elseif (enable_annual_averages) then
        tt=tt-365./2.
      endif

      print*,' writing to file ',outfile

      ncid = nccre (outfile, NCCLOB, iret)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
c       dimensions
      lon_tdim  = ncddef(ncid, 'Longitude_t', imt, iret)
      Lon_udim  = ncddef(ncid, 'Longitude_u', imt, iret)
      Lat_tdim  = ncddef(ncid, 'Latitude_t',  jmt, iret)
      Lat_udim  = ncddef(ncid, 'Latitude_u',  jmt, iret)
      depth_wdim = ncddef(ncid, 'depth_w',  km, iret)
      depth_tdim = ncddef(ncid, 'depth_t',  km, iret)
      iTimedim  = ncddef(ncid, 'Time', 1, iret)
c     grid variables
      dims(1)  = Lon_tdim
      Lon_tid  = ncvdef (ncid,'Longitude_t',NCFLOAT,1,dims,iret)
      dims(1)  = Lon_udim
      Lon_uid  = ncvdef (ncid,'Longitude_u',NCFLOAT,1,dims,iret)
      dims(1)  = Lat_tdim
      Lat_tid  = ncvdef (ncid,'Latitude_t', NCFLOAT,1,dims,iret)
      dims(1)  = Lat_udim
      Lat_uid  = ncvdef (ncid,'Latitude_u', NCFLOAT,1,dims,iret)
      dims(1)  = iTimedim
      iTimeid   = ncvdef(ncid,'Time',       NCFLOAT,1,dims,iret)
      dims(1)  = depth_wdim
      depth_wid = ncvdef (ncid,'depth_w', NCFLOAT,1,dims,iret)
      dims(1)  = depth_tdim
      depth_tid = ncvdef (ncid,'depth_t', NCFLOAT,1,dims,iret)
c     attributes of the grid
      name = 'Longitude on T grid     '; unit = 'degrees_W       '
      call ncaptc(ncid, Lon_tid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lon_tid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Longitude on U grid     '; unit = 'degrees_W       '
      call ncaptc(ncid, Lon_uid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lon_uid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Latitude on T grid      '; unit = 'degrees_N       '
      call ncaptc(ncid, Lat_tid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lat_tid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Latitude on U grid      '; unit = 'degrees_N       '
      call ncaptc(ncid, Lat_uid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lat_uid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Depth of T Grid points  '; unit = 'm '
      call ncaptc(ncid, depth_tid, 'long_name',NCCHAR, 24, name, iret) 
      call ncaptc(ncid, depth_tid, 'units',    NCCHAR, 16, unit, iret) 
      name = 'Depth of W Grid points  '; unit = 'm '
      call ncaptc(ncid, depth_wid, 'long_name', NCCHAR,24, name, iret) 
      call ncaptc(ncid, depth_wid, 'units',     NCCHAR,16, unit, iret) 
      name = 'Time                    '; unit = 'days            '
      call ncaptc(ncid, iTimeid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, iTimeid, 'units',     NCCHAR, 16, unit, iret) 
      call ncaptc(ncid, iTimeid,'time_origin',NCCHAR, 20,
     &'31-DEC-1899 00:00:00', iret)
c     variables
      dims=(/lon_udim,lat_udim,depth_tdim,itimedim/)
      uid=ncvdef (ncid,'u', NCFLOAT,4,dims,iret)
      vid=ncvdef (ncid,'v', NCFLOAT,4,dims,iret)
      do n=1,nt
         name='temp';if (n==2) name='salt'
         if (n>2) then
          write(name,'("tracer_",i2)') n
          call replace_space_zero(name)
         endif
         dims=(/lon_tdim,lat_tdim,depth_tdim,itimedim/)
         trid(n)=ncvdef (ncid,name, NCFLOAT,4,dims,iret)
         name='heat_fl';if (n==2) name='salt_fl'
         if (n>2) then
          write(name,'("tr",i2,"_flux")') n
          call replace_space_zero(name)
         endif
         dims=(/lon_tdim,lat_tdim,itimedim,1/)
         strid(n)=ncvdef (ncid,name, NCFLOAT,3,dims,iret)
      enddo
      dims=(/lon_udim,lat_udim,itimedim,1/)
      tauxid=ncvdef (ncid,'taux', NCFLOAT,3,dims,iret)
      tauyid=ncvdef (ncid,'tauy', NCFLOAT,3,dims,iret)
      dims=(/lon_tdim,lat_tdim,itimedim,1/)
      if (enable_freesurf) then
       psiid=ncvdef (ncid,'ssh', NCFLOAT,3,dims,iret)
      else
       psiid=ncvdef (ncid,'psi', NCFLOAT,3,dims,iret)
      endif
      if (found_kappa) then
       dims=(/lon_tdim,lat_tdim,depth_wdim,itimedim/)
       diff_cbtid=ncvdef (ncid,'diff_cbt', NCFLOAT,4,dims,iret)
      endif
      if (found_gm) then
       dims=(/lon_udim,lat_tdim,depth_tdim,itimedim/)
       ugm_id=ncvdef (ncid,'u_gm', NCFLOAT,4,dims,iret)
       dims=(/lon_tdim,lat_udim,depth_tdim,itimedim/)
       vgm_id=ncvdef (ncid,'v_gm', NCFLOAT,4,dims,iret)
       dims=(/lon_tdim,lat_tdim,depth_wdim,itimedim/)
       wgm_id=ncvdef (ncid,'w_gm', NCFLOAT,4,dims,iret)
      endif
      if (found_pco2) then
         dims=(/lon_tdim,lat_tdim,itimedim,1/)
         pco2_id=ncvdef (ncid,'pco2_water', NCFLOAT,3,dims,iret)
         dims=(/lon_tdim,lat_tdim,itimedim,1/)
         dco2_id=ncvdef (ncid,'dco2', NCFLOAT,3,dims,iret)
         dims=(/itimedim,1,1,1/)
         airco2_id=ncvdef (ncid,'pco2_air', NCFLOAT,1,dims,iret)
      endif
      if (found_o2) then
       dims=(/lon_tdim,lat_tdim,itimedim,1/)
       o2_id=ncvdef (ncid,'o2_sat', NCFLOAT,3,dims,iret)
      endif
c     attributes of variables
      name = 'Zonal velocity'; unit = 'm/s'
      call dvcdf_here(ncid,uid,name,24,unit,16,spval)
      name = 'Meridional velocity'; unit = 'm/s'
      call dvcdf_here(ncid,vid,name,24,unit,16,spval)
      do n=1,nt
         name = 'Potential temperature'; unit = 'deg C'
         if (n==2) then
          name = 'Salinity'; unit = 'psu'
         elseif (n>2) then
          write(name,'("Passive tracer ",i2)') n
          unit = '?'
         endif
         call dvcdf_here(ncid,trid(n),name,len_trim(name),
     &              unit,len_trim(unit),spval)
         name = 'Surface heat flux'; unit = 'W/m^2'
         if (n==2) then
          name = 'Surface salt flux'; unit = 'kg/m^2/s'
         elseif (n>2) then
          write(name,'("Surface flux passive tracer ",i2)') n
          unit = '?'
         endif
         call dvcdf_here(ncid,strid(n),name,len_trim(name),
     &              unit,len_trim(unit),spval)
      enddo
      if (enable_freesurf) then
       name = 'Sea Surface height'; unit = 'm'
      else
       name = 'Barotropic streamfunction'; unit = 'Sv'
      endif
      call dvcdf_here(ncid,psiid,name,24,unit,16,spval)
      name = 'Zonal wind stress'; unit = 'N/m^2'
      call dvcdf_here(ncid,tauxid,name,24,unit,16,spval)
      name = 'Meridional wind stress'; unit = 'N/m^2'
      call dvcdf_here(ncid,tauyid,name,24,unit,16,spval)

      if (found_kappa) then
       name = 'Vertical diffusivity'; unit = 'm^2/s'
       call dvcdf_here(ncid,diff_cbtid,name,24,unit,16,spval)
      endif
      if (found_gm) then
       name = 'Zonal bolus velocity'; unit = 'm/s'
       call dvcdf_here(ncid,ugm_id,name,24,unit,16,spval)
       name = 'Meridional bolus velocity'; unit = 'm/s'
       call dvcdf_here(ncid,vgm_id,name,24,unit,16,spval)
       name = 'Vertical bolus velocity'; unit = 'm/s'
       call dvcdf_here(ncid,wgm_id,name,24,unit,16,spval)
      endif
      if (found_pco2) then
         name = 'Sea surface CO2 partial pressure'; unit = 'ppmv'
         call dvcdf_here(ncid,pco2_id,name,
     &            len_trim(name),unit,len_trim(unit),spval)
         name = 'Air/Sea difference in CO2 concentration'
          unit = 'mol/m^3'
         call dvcdf_here(ncid,dco2_id,name,
     &            len_trim(name),unit,len_trim(unit),spval)
         name = 'Atmospheric CO2 partial pressure'; unit = 'ppmv'
         call dvcdf_here(ncid,airco2_id,name,
     &            len_trim(name),unit,len_trim(unit),spval)
      endif
      if (found_o2) then
         name = 'Sea surface O2 saturation concentration'; 
         unit = 'mol/m^3'
         call dvcdf_here(ncid,o2_id,name,
     &            len_trim(name),unit,len_trim(unit),spval)
      endif
      iret= nf_put_att_int(ncid,nf_global,
     &   'Number_of_averaged_time_steps ',nf_int,1,counter)
        if (iret.ne.0) print*,nf_strerror(iret)

      write( name,'("SFLAME sub domain # ",i3)') domain
      if (enable_daily_averages) 
     &   name=name(1:len_trim(name))//' daily means'
      if (enable_monthly_averages) 
     &   name=name(1:len_trim(name))//' monthly means'
      if (enable_seasonal_averages) 
     &   name=name(1:len_trim(name))//' seasonal means'
      if (enable_annual_averages) 
     &   name=name(1:len_trim(name))//' annual means'
      len=len_trim(name)
      iret= nf_put_att_text(ncid,nf_global,'history',len,name)

      call ncendf(ncid, iret)

      zt=zt/100.; zw=zw/100.

#ifdef real_8
      iret= nf_put_vara_double (ncid,lon_tid,1,imt,xt)
      iret= nf_put_vara_double (ncid,lon_uid,1,imt,xu)
      iret= nf_put_vara_double (ncid,lat_tid,1,jmt,yt)
      iret= nf_put_vara_double (ncid,lat_uid,1,jmt,yu)
      iret= nf_put_vara_double (ncid,depth_tid,1,km,zt)
      iret= nf_put_vara_double (ncid,depth_wid,1,km,zw)
      iret= nf_put_vara_double (ncid,itimeid,1,1,tt)
#else
      iret= nf_put_vara_real (ncid,lon_tid,1,imt,xt)
      iret= nf_put_vara_real (ncid,lon_uid,1,imt,xu)
      iret= nf_put_vara_real (ncid,lat_tid,1,jmt,yt)
      iret= nf_put_vara_real (ncid,lat_uid,1,jmt,yu)
      iret= nf_put_vara_real (ncid,depth_tid,1,km,zt)
      iret= nf_put_vara_real (ncid,depth_wid,1,km,zw)
      iret= nf_put_vara_real (ncid,itimeid,1,1,tt)
#endif


c     convert to SI units and divide by N
#ifdef low_memory
#else
      mean_u=mean_u/100.   /counter
      mean_t = mean_t/counter
      mean_t(:,:,:,2)=mean_t(:,:,:,2)*1000.+35.
#endif
      mean_stf=mean_stf/counter
      mean_stf(:,:,1)=mean_stf(:,:,1)*41868.
      mean_stf(:,:,2)=mean_stf(:,:,2)*10. 

      mean_smf=mean_smf/10. /counter
      if (enable_freesurf) then
       mean_psi=mean_psi/100. /counter
      else
       mean_psi=mean_psi/1.e12 /counter
      endif
c     apply land mask and write to file
      do k=1,km
#ifdef low_memory
         call read_tsuv(infile,k)
         kk=1
      mean_u=mean_u/100.   /counter
      mean_t = mean_t/counter
      mean_t(:,:,:,2)=mean_t(:,:,:,2)*1000.+35.
#else
         kk=k
#endif
         do n=1,nt
          where( kmt <k ) mean_t(:,kk,:,n) = spval
         enddo
         do n=1,2
          where( kmu <k ) mean_u(:,kk,:,n) = spval
         enddo
         start=(/1,1,k,1/); count=(/imt,jmt,1,1/)
#ifdef real_8
         iret=nf_put_vara_double(ncid,uid,start,count,mean_u(:,kk,:,1))
         iret=nf_put_vara_double(ncid,vid,start,count,mean_u(:,kk,:,2))
#else
         iret= nf_put_vara_real (ncid,uid,start,count,mean_u(:,kk,:,1))
         iret= nf_put_vara_real (ncid,vid,start,count,mean_u(:,kk,:,2))
#endif
         do n=1,nt
#ifdef real_8
          iret=nf_put_vara_double(ncid,trid(n),
     &                          start,count,mean_t(:,kk,:,n))
#else
          iret=nf_put_vara_real(ncid,trid(n),
     &                          start,count,mean_t(:,kk,:,n))
#endif
         enddo
      enddo
      where( kmt == 0) mean_psi = spval
      do n=1,nt
       where( kmt == 0) mean_stf(:,:,n) = spval
      enddo
      do n=1,2
       where( kmu == 0) mean_smf(:,:,n) = spval
      enddo
      start=(/1,1,1,1/); count=(/imt,jmt,1,1/)
#ifdef real_8
      iret= nf_put_vara_double (ncid,psiid,start,count,mean_psi)
#else
      iret= nf_put_vara_real (ncid,psiid,start,count,mean_psi)
#endif
      do n=1,nt
#ifdef real_8
      iret=nf_put_vara_double(ncid,strid(n),start,count,mean_stf(:,:,n))
#else
       iret=nf_put_vara_real(ncid,strid(n),start,count,mean_stf(:,:,n))
#endif
      enddo
#ifdef real_8
      iret= nf_put_vara_double (ncid,tauxid,start,count,mean_smf(:,:,1))
      iret= nf_put_vara_double (ncid,tauyid,start,count,mean_smf(:,:,2))
#else
      iret= nf_put_vara_real (ncid,tauxid,start,count,mean_smf(:,:,1))
      iret= nf_put_vara_real (ncid,tauyid,start,count,mean_smf(:,:,2))
#endif
      if (found_kappa) then
       mean_ka=mean_ka/1.0e4 /counter
       do k=1,km
        where( kmt <k ) mean_ka(:,k,:) = spval
        start=(/1,1,k,1/); count=(/imt,jmt,1,1/)
#ifdef real_8
         iret= nf_put_vara_double (ncid,
     &                 diff_cbtid,start,count,mean_ka(:,k,:))
#else
         iret= nf_put_vara_real (ncid,
     &                 diff_cbtid,start,count,mean_ka(:,k,:))
#endif
       enddo
      endif
      if (found_gm) then
       mean_ugm=mean_ugm/100.0/counter
       do k=1,km
        where( kmt <k ) mean_ugm(:,k,:,1) = spval
        where( kmt <k ) mean_ugm(:,k,:,2) = spval
        where( kmt <k ) mean_ugm(:,k,:,3) = spval
        start=(/1,1,k,1/); count=(/imt,jmt,1,1/)
#ifdef real_8
        iret= nf_put_vara_double(ncid,
     &           ugm_id,start,count,mean_ugm(:,k,:,1))
        iret= nf_put_vara_double(ncid,
     &           vgm_id,start,count,mean_ugm(:,k,:,2))
        iret= nf_put_vara_double(ncid,
     &           wgm_id,start,count,mean_ugm(:,k,:,3))
#else
        iret= nf_put_vara_real(ncid,
     &           ugm_id,start,count,mean_ugm(:,k,:,1))
        iret= nf_put_vara_real(ncid,
     &           vgm_id,start,count,mean_ugm(:,k,:,2))
        iret= nf_put_vara_real(ncid,
     &           wgm_id,start,count,mean_ugm(:,k,:,3))
#endif
       enddo        
      endif
      if (found_pco2) then
        mean_pco2=mean_pco2/counter
        mean_dco2=mean_dco2/counter
        mean_pco2_air=mean_pco2_air/counter
        where(kmt==0.) mean_pco2 = spval
        where(kmt==0.) mean_dco2 = spval
        start=(/1,1,1,1/); count=(/imt,jmt,1,1/)
#ifdef real_8
        iret= nf_put_vara_double(ncid,pco2_id,start,count,mean_pco2)
        iret= nf_put_vara_double(ncid,dco2_id,start,count,mean_dco2)
        iret= nf_put_vara_double(ncid,airco2_id,1,1,mean_pco2_air )
#else
        iret= nf_put_vara_real(ncid,pco2_id,start,count,mean_pco2)
        iret= nf_put_vara_real(ncid,dco2_id,start,count,mean_dco2)
        iret= nf_put_vara_real(ncid,airco2_id,1,1,mean_pco2_air )
#endif
      endif
      if (found_pco2) then
        mean_o2=mean_o2/counter
        where(kmt==0.) mean_o2 = spval
        start=(/1,1,1,1/); count=(/imt,jmt,1,1/)
#ifdef real_8
        iret= nf_put_vara_double(ncid,o2_id,start,count,mean_o2)
#else
        iret= nf_put_vara_real(ncid,o2_id,start,count,mean_o2)
#endif
      endif
      call ncclos (ncid, iret)
      end subroutine write_averages







      subroutine read_tsuv(infile,kk)
      use cv_ave_module
      implicit none
      character (len=*) :: infile
      character (len=3) :: ssh
      integer :: i,j,k,io,n,kk
      character (len=80) :: ident
      logical :: ok

      io=10
      print*,'reading from file ',infile(1:len_trim(infile)),
     &      ' UVTS at level ',kk
      open(io,file=infile,form='unformatted')
      read(io) !i,j,k,n
      read(io) !xt,yt,xu,yu,zt,zw
      read(io) !counter
      read(io) !enable_daily_averages
      read(io) !enable_monthly_averages
      read(io) !enable_seasonal_averages
      read(io) !enable_annual_averages
      read(io) !kmt
      read(io) !mean_psi,ssh
      do n=1,2
       read(io) !mean_smf(:,:,n)
       do k=1,kk-1
         read(io) 
       enddo
       read(io) mean_u(1:imtf,1,:,n)
#ifdef unfold
       mean_u(imtf-1:imt,1,:,n)= mean_u(1:imt-imtf+2,1,:,n)
#endif
       do k=kk+1,km
         read(io) 
       enddo
      enddo
      do n=1,nt
       read(io) !mean_stf(:,:,n)
       do k=1,kk-1
         read(io) !mean_t(:,k,:,n)
       enddo
       read(io) mean_t(1:imtf,1,:,n)
#ifdef unfold
       mean_t(imtf-1:imt,1,:,n)= mean_t(1:imt-imtf+2,1,:,n)
#endif
       do k=kk+1,km
         read(io) 
       enddo
      enddo

      close(io)
      end subroutine read_tsuv












      subroutine dvcdf_here(ncid,ivarid,name,iname,unit,iunit,spval)
c-----------------------------------------------------------------------
c define some standard attributes of variable ivarid in NetCDF file ncid 
c-----------------------------------------------------------------------
      use kind_mod
      implicit none
      integer ncid,ivarid,iname,iunit,iret
      character*(*) name, unit
      real (kind=r_knd):: spval
#include "netcdf.inc"
      call ncaptc(ncid,ivarid, 'long_name', NCCHAR,iname , name, iret) 
      call ncaptc(ncid,ivarid, 'units',     NCCHAR,iunit, unit, iret) 
#ifdef real_8
      iret=nf_put_att_double(NCID,iVARID, 'missing_value', NCFLOAT, 
     &                       1,spval)
      iret=nf_put_att_double(NCID,iVARID, '_FillValue', NCFLOAT, 
     &                       1,spval)
#else
      iret=nf_put_att_real(NCID,iVARID, 'missing_value', NCFLOAT, 
     &                       1,spval)
      iret=nf_put_att_real(NCID,iVARID, '_FillValue', NCFLOAT, 
     &                       1,spval)
#endif
c        if (iret.ne.0) print*,nf_strerror(iret)
      end subroutine dvcdf_here


      subroutine barrier
      end

      subroutine halt_stop(s)
      character*(*) s
      print*,s
      stop
      end


