#include "options.inc"


       module averages_module
c
c=======================================================================
c      perform averages over U/V T and S
c      and write to NetCDF or binary files
c      SPFLAME version:   c.eden
c=======================================================================
c
c      local cpp directives:
c
c      write vertical diffusivities:
#define write_vert_diff
c
c      write GM velocities
#define write_gm_vel
c
c      write CO2 partial pressure at the surface
#define write_pco2
c
c      write oxygen saturation concentration at the surface
#define write_o2
c
c      write tracer stuff in separate netcdf file
c#define splitup_netcdf_file
c
       implicit none
 
       real, allocatable, dimension(:,:,:,:) :: umean,tmean
       real, allocatable, dimension(:,:)     :: psimean
       real, allocatable, dimension(:,:,:)   :: stfmean,smfmean
       integer :: counter 
       logical :: active = .false.
       real, allocatable, dimension(:,:,:) :: ka_mean
       real, allocatable, dimension(:,:,:,:) :: ugm_mean
       real, allocatable, dimension(:,:) :: pco2_mean,o2_mean,dco2_mean
       real :: airco2_mean
       logical :: gm_active = .false.
       logical :: pco2_active = .false.
       logical :: o2_active   = .false.
       end module averages_module


       subroutine init_averages
       use spflame_module
       use averages_module
       implicit none
       if (enable_daily_averages .or.
     &     enable_monthly_averages .or.
     &     enable_seasonal_averages .or.
     &     enable_annual_averages ) active = .true.
       if (active) then
        if (my_pe==0) print*,' initialising diagnostic averages'
        allocate( umean(is_pe:ie_pe,km,js_pe:je_pe,2) )
        allocate( tmean(is_pe:ie_pe,km,js_pe:je_pe,nt) )
        allocate( stfmean(is_pe:ie_pe,js_pe:je_pe,nt) )
        allocate( smfmean(is_pe:ie_pe,js_pe:je_pe,2) )
        allocate( psimean(is_pe:ie_pe,js_pe:je_pe) )
        umean=0.; tmean=0.; stfmean=0.;smfmean=0.; psimean=0.
        counter = 0
#ifdef write_vert_diff
        if (my_pe==0) print*,' (averaging also vertical diffusivity)'
        allocate( ka_mean(is_pe:ie_pe,km,js_pe:je_pe) ); ka_mean=0.
#endif
#ifdef write_gm_vel
        if (enable_diffusion_isopycnic) gm_active = .true.
#endif
#ifdef write_pco2
        pco2_active = .true.
#endif
#ifdef write_o2
        o2_active = .true.
#endif

        if (gm_active) then
         if (my_pe==0) print*,' (averaging also GM bolus velocities)'
         allocate( ugm_mean(is_pe:ie_pe,km,js_pe:je_pe,3) ); ugm_mean=0.
        endif
        if (pco2_active) then
         if (my_pe==0) print*,' (averaging also surface pCO2)'
         allocate( pco2_mean(is_pe:ie_pe,js_pe:je_pe) ); pco2_mean=0.
         allocate( dco2_mean(is_pe:ie_pe,js_pe:je_pe) ); dco2_mean=0.
         airco2_mean=0.0
        endif
        if (o2_active) then
         if (my_pe==0) print*,' (averaging also surface O2 saturation)'
         allocate( o2_mean(is_pe:ie_pe,js_pe:je_pe) ); o2_mean=0.
        endif

        call read_unfinished_averages()
        if (my_pe==0) print*,' done initialising diagnostic averages'
       endif
       end subroutine init_averages


       subroutine averages
       use spflame_module
       use averages_module
       use freesurf_module
       use isopycnic_module
       use passive_tracer_module
       use npzd_module
       implicit none
       if (active) then
        if (.not. eulerback1 ) then
         counter=counter+1
         umean=umean+u(is_pe:ie_pe,1:km,js_pe:je_pe,1:2,tau)
         tmean=tmean+t(is_pe:ie_pe,1:km,js_pe:je_pe,1:nt,tau)
         if (enable_freesurf) then
          psimean=psimean+etat(is_pe:ie_pe,js_pe:je_pe,3)
         else
          psimean=psimean+psi(is_pe:ie_pe,js_pe:je_pe,1)
         endif
         stfmean=stfmean+stf(is_pe:ie_pe,js_pe:je_pe,1:nt)
         smfmean=smfmean+smf(is_pe:ie_pe,js_pe:je_pe,1:2)
#ifdef write_vert_diff
         ka_mean=ka_mean+diff_cbt(is_pe:ie_pe,1:km,js_pe:je_pe)
#endif
         if (gm_active) then
          ugm_mean(:,:,:,1)=ugm_mean(:,:,:,1)
     &              +adv_vetiso(is_pe:ie_pe,1:km,js_pe:je_pe)
          ugm_mean(:,:,:,2)=ugm_mean(:,:,:,2)
     &              +adv_vntiso(is_pe:ie_pe,1:km,js_pe:je_pe)
          ugm_mean(:,:,:,3)=ugm_mean(:,:,:,3)
     &              +adv_vbtiso(is_pe:ie_pe,1:km,js_pe:je_pe)
         endif
         if (pco2_active .and. enable_co2) then
          pco2_mean(is_pe:ie_pe,js_pe:je_pe)= 
     &       pco2_mean(is_pe:ie_pe,js_pe:je_pe)+
     &       pco2_water(is_pe:ie_pe,js_pe:je_pe)
          dco2_mean(is_pe:ie_pe,js_pe:je_pe)= 
     &       dco2_mean(is_pe:ie_pe,js_pe:je_pe)+
     &       diff_co2(is_pe:ie_pe,js_pe:je_pe)
          airco2_mean = airco2_mean + pco2_air
         endif
         if (o2_active .and. enable_o2) then
          o2_mean(is_pe:ie_pe,js_pe:je_pe)= 
     &       o2_mean(is_pe:ie_pe,js_pe:je_pe)
     &      +o2_sat(is_pe:ie_pe,js_pe:je_pe)
         endif
        endif ! eulerback1

        if ((end_of_day   .and. enable_daily_averages).or. 
     &      (end_of_month .and. enable_monthly_averages).or. 
     &      (end_of_season.and. enable_seasonal_averages).or. 
     &      (end_of_year  .and. enable_annual_averages) ) then
         call write_averages_cdf()
         counter=0 
         umean=0.; tmean=0.; psimean=0.;stfmean=0.;smfmean=0.
#ifdef write_vert_diff
         ka_mean=0.
#endif
         if (gm_active) ugm_mean=0.
         if (pco2_active) then
              pco2_mean=0.0; dco2_mean=0.0; airco2_mean = 0.0
         endif
         if (o2_active) o2_mean=0.0
        endif
        if (last_time_step) call write_unfinished_averages()
       endif
       end subroutine averages


       subroutine write_unfinished_averages()
       use spflame_module
       use averages_module
       implicit none
       character (len=80) :: name,ident
       integer i,j,k,io,n,is,ie,js,je
       real, allocatable ::  buf(:,:)
       is=is_pe; ie=ie_pe; js=js_pe; je=je_pe
       write(name,'("averages_",i3,".dta")') sub_domain
       do i=1,len_trim(name); if (name(i:i) == ' ') name(i:i)='0'; enddo
       if (my_pe==0) then
        print*,''
        print*,' writing unfinished averages of U/T/S to file ',
     &          name(1:len_trim(name))
        call get_free_iounit(io)
        open(io,file=name,form='unformatted',status='unknown')
        write(io) counter
       endif
       allocate( buf(imt,jmt) )
       buf(is:ie,js:je)=psimean(is:ie,js:je)
       call pe0_recv_2D(buf)
       if (my_pe==0) write(io) buf
       do n=1,2
        buf(is:ie,js:je)=smfmean(is:ie,js:je,n)
        call pe0_recv_2D(buf)
        if (my_pe==0) write(io) buf
       enddo
       do n=1,nt
        buf(is:ie,js:je)=stfmean(is:ie,js:je,n)
        call pe0_recv_2D(buf)
        if (my_pe==0) write(io) buf
       enddo
       do k=1,km
        do n=1,2
         buf(is:ie,js:je)=umean(is:ie,k,js:je,n)
         call pe0_recv_2D(buf)
         if (my_pe==0) write(io) buf
        enddo
        do n=1,nt
         buf(is:ie,js:je)=tmean(is:ie,k,js:je,n)
         call pe0_recv_2D(buf)
         if (my_pe==0) write(io) buf
        enddo
       enddo
#ifdef write_vert_diff
       ident = 'diff_cbt'
       if (my_pe==0) write(io) ident
       do k=1,km
         buf(is:ie,js:je)=ka_mean(is:ie,k,js:je)
         call pe0_recv_2D(buf)
         if (my_pe==0) write(io) buf
       enddo
#endif
       if (gm_active) then
        ident = 'u_gm'
        if (my_pe==0) write(io) ident
        do k=1,km
         buf(is:ie,js:je)=ugm_mean(is:ie,k,js:je,1)
         call pe0_recv_2D(buf)
         if (my_pe==0) write(io) buf
         buf(is:ie,js:je)=ugm_mean(is:ie,k,js:je,2)
         call pe0_recv_2D(buf)
         if (my_pe==0) write(io) buf
         buf(is:ie,js:je)=ugm_mean(is:ie,k,js:je,3)
         call pe0_recv_2D(buf)
         if (my_pe==0) write(io) buf
        enddo
       endif
       if (pco2_active) then
        ident = 'pco2_water'
        if (my_pe==0) write(io) ident
        buf(is:ie,js:je)=pco2_mean(is:ie,js:je)
        call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
        buf(is:ie,js:je)=dco2_mean(is:ie,js:je)
        call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
                                     if (my_pe==0) write(io) airco2_mean
       endif
       if (o2_active) then
        ident = 'o2_sat'
        if (my_pe==0) write(io) ident
        buf(is:ie,js:je)=o2_mean(is:ie,js:je)
        call pe0_recv_2D(buf)
        if (my_pe==0) write(io) buf
       endif
       if (my_pe==0) close(io)
       deallocate( buf )
       end subroutine write_unfinished_averages


       subroutine read_unfinished_averages()
       use spflame_module
       use averages_module
       implicit none
       character (len=80) :: name,ident
       integer i,j,k,io,n,is,ie,js,je
       real, allocatable ::  buf(:,:)
       logical :: ok
       is=is_pe; ie=ie_pe; js=js_pe; je=je_pe
       write(name,'("averages_",i3,".dta")') sub_domain
       do i=1,len_trim(name); if (name(i:i) == ' ') name(i:i)='0'; enddo
       if (my_pe==0) then
        print*,''
        print*,' reading averaged U/T/S from file ',
     &          name(1:len_trim(name))
c        call get_free_iounit(io)
       endif

       call get_free_iounit(io)
       open(io,file=name,form='unformatted',status='old',err=200)
       read(io) counter
       allocate( buf(imt,jmt) )
       if (my_pe==0) read(io) buf
       call pe0_send_2D(buf)
       psimean(is:ie,js:je)=buf(is:ie,js:je)
       do n=1,2
        if (my_pe==0) read(io) buf
        call pe0_send_2D(buf)
        smfmean(is:ie,js:je,n)=buf(is:ie,js:je)
       enddo
       do n=1,nt
        if (my_pe==0) read(io) buf
        call pe0_send_2D(buf)
        stfmean(is:ie,js:je,n)=buf(is:ie,js:je)
       enddo
       do k=1,km
        do n=1,2
         if (my_pe==0) read(io) buf
         call pe0_send_2D(buf)
         umean(is:ie,k,js:je,n)=buf(is:ie,js:je)
        enddo
        do n=1,nt
         if (my_pe==0) read(io) buf
         call pe0_send_2D(buf)
         tmean(is:ie,k,js:je,n)= buf(is:ie,js:je)
        enddo
       enddo


 10    ident=''; ok = .true.
       if (my_pe==0) then 
        read(io,end=20) ident
        goto 30
 20     ok = .false.
 30     continue
       endif
       call bcast_logical(ok,1,0)
       call bcast_char(ident,80,0)

       if (ident(1:8) == 'diff_cbt') then
        do k=1,km
         if (my_pe==0) read(io) buf
#ifdef write_vert_diff
         call pe0_send_2D(buf)
         ka_mean(is:ie,k,js:je)= buf(is:ie,js:je)
#endif
        enddo
       elseif (ident(1:4) == 'u_gm') then
        do k=1,km
         if (my_pe==0) read(io) buf
         if (gm_active) then
          call pe0_send_2D(buf)
          ugm_mean(is:ie,k,js:je,1)= buf(is:ie,js:je)
         endif
         if (my_pe==0) read(io) buf
         if (gm_active) then
          call pe0_send_2D(buf)
          ugm_mean(is:ie,k,js:je,2)= buf(is:ie,js:je)
         endif
         if (my_pe==0) read(io) buf
         if (gm_active) then
          call pe0_send_2D(buf)
          ugm_mean(is:ie,k,js:je,3)= buf(is:ie,js:je)
         endif
        enddo
       elseif (ident(1:10) == 'pco2_water' ) then
        if (my_pe==0) read(io) buf; call pe0_send_2D(buf)
        if (pco2_active) pco2_mean(is:ie,js:je)=buf(is:ie,js:je)
        if (my_pe==0) read(io) buf; call pe0_send_2D(buf)
        if (pco2_active) dco2_mean(is:ie,js:je)=buf(is:ie,js:je)
        if (my_pe==0) read(io) airco2_mean
        call bcast_real(airco2_mean,1,0)
       elseif (ident(1:6) == 'o2_sat' ) then
        if (my_pe==0) read(io) buf; call pe0_send_2D(buf)
        if (o2_active) o2_mean(is:ie,js:je)=buf(is:ie,js:je)
       else
         if (my_pe==0) print*,' cannot read ',ident,' in ',name
       endif
       if (ok) goto 10  ! read next item til end of file

       close(io)
       if (my_pe==0) then
        print*,' done '
        print*,''
       endif
       deallocate( buf )
       return

 200   if (my_pe==0) then
        print*,''
        print*,'---------------------------------------------'
        print*,' WARNING : cannot read file ',name(1:len_trim(name))
        print*,'---------------------------------------------'
        print*,''
       endif

       end subroutine read_unfinished_averages


       subroutine write_averages_cdf()
       use spflame_module
       use averages_module
       implicit none
#ifdef netcdf_diagnostics
c
c      write averaged quantities to a nice NetCDF file
#ifdef splitup_netcdf_file
c     but not the tracer stuff
#endif
c
       character (len=80) :: fname,name,unit
       integer :: year,month,day
#include "netcdf.inc"
       integer :: ncid,iret,i,j,k,dims(4),start(4),count(4)
       real    :: spval=-9.9e12, tt
       integer :: lon_tdim,lon_tid,lon_udim,lon_uid
       integer :: lat_tdim,lat_tid,lat_udim,lat_uid
       integer :: depth_tdim,depth_tid,depth_wdim,depth_wid
       integer :: itimedim,itimeid,uid,vid,trid(nt),npe,n
       integer :: psiid,tauxid, tauyid,strid(nt),diff_cbtid
       integer :: ugm_id,vgm_id,wgm_id,pco2_id,o2_id,dco2_id,airco2_id
       integer :: is,ie,js,je
#ifdef netcdf_real4
      real (kind=4) :: var(is_pe:ie_pe,js_pe:je_pe)
      real (kind=4) , allocatable :: v2(:)
#else
      real          :: var(is_pe:ie_pe,js_pe:je_pe)
      real          , allocatable :: v2(:)
#endif

       is=is_pe; ie=ie_pe; js=js_pe; je=je_pe

       year  = get_year(get_current_time())
       month = get_month_of_year(get_current_time())
       day   = get_day_of_month(get_current_time())
       write(fname,'("averages_",i3,"_y",i4,"m",i2,"d",i2,".cdf")') 
     &            sub_domain,year,month,day
       do i=1,len_trim(fname);if (fname(i:i)==' ')fname(i:i)='0';enddo

       if (my_pe == 0) then
        print*,' Writing averages to  NetCDF output file ',
     &         fname(1:len_trim(fname))
        ncid = nccre (fname, NCCLOB, iret)
        iret=nf_set_fill(ncid, NF_NOFILL, iret)
        call store_info_cdf(ncid)
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)
#ifdef splitup_netcdf_file
        do n=1,2
#else
        do n=1,nt
#endif
         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
#ifdef write_vert_diff
        dims=(/lon_tdim,lat_tdim,depth_wdim,itimedim/)
        diff_cbtid=ncvdef (ncid,'diff_cbt', NCFLOAT,4,dims,iret)
#endif
        if (gm_active) 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
#ifndef splitup_netcdf_file
        if (pco2_active) 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 (o2_active) then
         dims=(/lon_tdim,lat_tdim,itimedim,1/)
         o2_id=ncvdef (ncid,'o2_sat', NCFLOAT,3,dims,iret)
        endif
#endif
c       attributes of variables
        name = 'Zonal velocity'; unit = 'm/s'
        call dvcdf(ncid,uid,name,24,unit,16,spval)
        name = 'Meridional velocity'; unit = 'm/s'
        call dvcdf(ncid,vid,name,24,unit,16,spval)
#ifdef splitup_netcdf_file
        do n=1,2
#else
        do n=1,nt
#endif
         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(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(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(ncid,psiid,name,24,unit,16,spval)
        name = 'Zonal wind stress'; unit = 'N/m^2'
        call dvcdf(ncid,tauxid,name,24,unit,16,spval)
        name = 'Meridional wind stress'; unit = 'N/m^2'
        call dvcdf(ncid,tauyid,name,24,unit,16,spval)
#ifdef write_vert_diff
        name = 'Vertical diffusivity'; unit = 'm^2/s'
        call dvcdf(ncid,diff_cbtid,name,24,unit,16,spval)
#endif
        if (gm_active) then
         name = 'Zonal bolus velocity'; unit = 'm/s'
         call dvcdf(ncid,ugm_id,name,24,unit,16,spval)
         name = 'Meridional bolus velocity'; unit = 'm/s'
         call dvcdf(ncid,vgm_id,name,24,unit,16,spval)
         name = 'Vertical bolus velocity'; unit = 'm/s'
         call dvcdf(ncid,wgm_id,name,24,unit,16,spval)
        endif
#ifndef splitup_netcdf_file
        if (pco2_active) then
         name = 'Sea surface CO2 partial pressure'; unit = 'ppmv'
         call dvcdf(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(ncid,dco2_id,name,
     &            len_trim(name),unit,len_trim(unit),spval)
         name = 'Atmospheric CO2 partial pressure'; unit = 'ppmv'
         call dvcdf(ncid,airco2_id,name,
     &            len_trim(name),unit,len_trim(unit),spval)
        endif
        if (o2_active) then
         name = 'Sea surface O2 saturation concentration'; 
         unit = 'mol/m^3'
         call dvcdf(ncid,o2_id,name,
     &            len_trim(name),unit,len_trim(unit),spval)
        endif
#endif
c       global attrubutes
        iret= nf_put_att_int(ncid,nf_global,
     &   'Number_of_averaged_time_steps ',nf_int,1,counter)
        call ncendf(ncid, iret)
c       write grid
        allocate( v2(max(imt,jmt,km)) )
        v2(1:imt)=xt
        call ncvpt(ncid, Lon_tid, 1, imt,v2, iret)
        v2(1:imt)=xu
        call ncvpt(ncid, Lon_uid, 1, imt,v2, iret)
        v2(1:jmt)=yt
        call ncvpt(ncid, Lat_tid, 1, jmt,v2, iret)
        v2(1:jmt)=yu
        call ncvpt(ncid, Lat_uid, 1, jmt,v2, iret)
        v2(1:km)=zt/100.
        call ncvpt(ncid, depth_tid, 1, km,v2, iret)
        v2(1:km)=zw/100.
        call ncvpt(ncid, depth_wid, 1, km,v2, iret)
        call read_stamp(current_stamp,tt ); v2(1)=tt
        call ncvpt(ncid, itimeid, 1, 1,v2, iret)
        deallocate(v2)
        call ncclos (ncid, iret)
       endif
c      loop over PEs
       do npe=0,n_pes
        call barrier
        if (my_pe==npe) then
         iret=nf_open(fname,NF_WRITE,ncid)
         iret=nf_set_fill(ncid, NF_NOFILL, iret)
         iret=nf_inq_varid(ncid,'u',uid)
         iret=nf_inq_varid(ncid,'v',vid)
#ifdef splitup_netcdf_file
         do n=1,2
#else
         do n=1,nt
#endif
          name='temp'
          if (n==2) then
           name='salt'
          elseif (n>2) then
           write(name,'("tracer_",i2)') n
           call replace_space_zero(name)
          endif
          iret=nf_inq_varid(ncid,name,trid(n))
          name = 'heat_fl'
          if (n==2) then
           name='salt_fl'
          elseif (n>2) then
           write(name,'("tr",i2,"_flux")') n
           call replace_space_zero(name)
          endif
          iret=nf_inq_varid(ncid,name,strid(n))
         enddo
         if (enable_freesurf) then
          iret=nf_inq_varid(ncid,'ssh',psiid)
         else
          iret=nf_inq_varid(ncid,'psi',psiid)
         endif
         iret=nf_inq_varid(ncid,'taux',tauxid)
         iret=nf_inq_varid(ncid,'tauy',tauyid)
#ifdef write_vert_diff
         iret=nf_inq_varid(ncid,'diff_cbt',diff_cbtid)
#endif
         if (gm_active) then
          iret=nf_inq_varid(ncid,'u_gm',ugm_id)
          iret=nf_inq_varid(ncid,'v_gm',vgm_id)
          iret=nf_inq_varid(ncid,'w_gm',wgm_id)
         endif
#ifndef splitup_netcdf_file
         if (pco2_active) then
          iret=nf_inq_varid(ncid,'pco2_water',pco2_id)
          iret=nf_inq_varid(ncid,'dco2',dco2_id)
          iret=nf_inq_varid(ncid,'pco2_air',airco2_id)
         endif
         if (o2_active) iret=nf_inq_varid(ncid,'o2_sat',o2_id)
#endif
c        divide with number of averaged time steps 
         umean=umean/counter
         psimean=psimean/counter
         smfmean=smfmean/counter
         tmean=tmean/counter
         stfmean=stfmean/counter
c        convert to SI units
         umean=umean/100.
         tmean(:,:,:,2)=tmean(:,:,:,2)*1000+35.
         stfmean(:,:,1)=stfmean(:,:,1)*41868.
         stfmean(:,:,2)=stfmean(:,:,2)*10.
         smfmean=smfmean/10.
         if (enable_freesurf) then
          psimean=psimean/100.
         else
          psimean=psimean/1.e12
         endif
c        apply land mask
         where( umask(is:ie,1:km,js:je) == 0. ) umean(:,:,:,1) = spval
         where( umask(is:ie,1:km,js:je) == 0. ) umean(:,:,:,2) = spval
         do n=1,nt
          where( tmask(is:ie,1:km,js:je) == 0. ) tmean(:,:,:,n) = spval
          where( tmask(is:ie,1,js:je) == 0. )    stfmean(:,:,n) = spval
         enddo
         where( tmask(is:ie,1,js:je) == 0. ) psimean = spval
         where( umask(is:ie,1,js:je) == 0. ) smfmean(:,:,1) = spval
         where( umask(is:ie,1,js:je) == 0. ) smfmean(:,:,2) = spval
#ifdef write_vert_diff
         ka_mean=ka_mean/counter
         ka_mean=ka_mean/1.0e4
         where( tmask(is:ie,1:km,js:je) == 0. ) ka_mean(:,:,:) = spval
#endif
         if (gm_active) then
          ugm_mean=ugm_mean/counter
          ugm_mean=ugm_mean/100.0
          do j=js,je
            ugm_mean(:,:,j,2)=ugm_mean(:,:,j,2)/csu(j)
          enddo
          where(tmask(is:ie,1:km,js:je)==0.) ugm_mean(:,:,:,1) = spval
          where(tmask(is:ie,1:km,js:je)==0.) ugm_mean(:,:,:,2) = spval
          where(tmask(is:ie,1:km,js:je)==0.) ugm_mean(:,:,:,3) = spval
         endif
         if (pco2_active) then
          pco2_mean=pco2_mean/counter
          dco2_mean=dco2_mean/counter
          airco2_mean=airco2_mean/counter
          where(tmask(is:ie,1,js:je)==0.) pco2_mean(is:ie,js:je) = spval
          where(tmask(is:ie,1,js:je)==0.) dco2_mean(is:ie,js:je) = spval
         endif
         if (o2_active) then 
          o2_mean=o2_mean/counter
          where(tmask(is:ie,1,js:je)==0.) o2_mean(is:ie,js:je) = spval
         endif
         ! write PE domain to file
         do k=1,km
          start=(/is,js,k,1/); count=(/ie-is+1,je-js+1,1,1/)
          var=umean(:,k,:,1)
          iret= nf_put_vara_real (ncid,uid,start,count,var )
          var=umean(:,k,:,2)
          iret= nf_put_vara_real (ncid,vid,start,count,var )
#ifdef splitup_netcdf_file
          do n=1,2
#else
          do n=1,nt
#endif
           var=tmean(:,k,:,n)
           iret= nf_put_vara_real (ncid,trid(n),start,count,var )
          enddo
         enddo
         start=(/is,js,1,1/); count=(/ie-is+1,je-js+1,1,1/)
         var=psimean
         iret= nf_put_vara_real (ncid,psiid,start,count,var)
#ifdef splitup_netcdf_file
         do n=1,2
#else
         do n=1,nt
#endif
          var=stfmean(:,:,n)
          iret= nf_put_vara_real (ncid,strid(n),start,count,var )
         enddo
         var=smfmean(:,:,1)
         iret= nf_put_vara_real (ncid,tauxid,start,count,var )
         var=smfmean(:,:,2)
         iret= nf_put_vara_real (ncid,tauyid,start,count,var )
#ifdef write_vert_diff
         do k=1,km
          start=(/is,js,k,1/); count=(/ie-is+1,je-js+1,1,1/)
          var=ka_mean(:,k,:)
          iret= nf_put_vara_real (ncid,diff_cbtid,start,count,var )
         enddo
#endif
         if (gm_active) then
          do k=1,km
           start=(/is,js,k,1/); count=(/ie-is+1,je-js+1,1,1/)
           var=ugm_mean(:,k,:,1)
           iret= nf_put_vara_real(ncid,ugm_id,start,count,var )
           var=ugm_mean(:,k,:,2)
           iret= nf_put_vara_real(ncid,vgm_id,start,count,var )
           var=ugm_mean(:,k,:,3)
           iret= nf_put_vara_real(ncid,wgm_id,start,count,var )
          enddo
         endif
#ifndef splitup_netcdf_file
         if (pco2_active) then
          var(is:ie,js:je)=pco2_mean(is:ie,js:je)
          start=(/is,js,1,1/); count=(/ie-is+1,je-js+1,1,1/)
          iret= nf_put_vara_real (ncid,pco2_id,start,count,var )
          var(is:ie,js:je)=dco2_mean(is:ie,js:je)
          iret= nf_put_vara_real (ncid,dco2_id,start,count,var )
          var(is,js)=airco2_mean
          iret= nf_put_vara_real(ncid,airco2_id,1,1,var(is,js) )
         endif
         if (o2_active) then
          var=o2_mean(:,:)
          start=(/is,js,1,1/); count=(/ie-is+1,je-js+1,1,1/)
          iret= nf_put_vara_real (ncid,o2_id,start,count,var )
         endif
#endif
         call ncclos (ncid, iret)
        endif
        call barrier
       enddo
#ifdef splitup_netcdf_file
       call write_averages_tracer
#endif
       if (my_pe==0) then
        print*,'done'
       endif
#else
c
c      write averaged quantities to a binary file
c
       integer :: is,ie,js,je
       integer :: year,month,day
       character (len=80) :: fname,ident
       integer :: io,n,i,k,j
       real, allocatable    :: buf(:,:)

       is=is_pe; ie=ie_pe; js=js_pe; je=je_pe

       year  = get_year(get_current_time())
       month = get_month_of_year(get_current_time())
       day   = get_day_of_month(get_current_time())
       write(fname,'("averages_",i3,"_y",i4,"m",i2,"d",i2,".dta")') 
     &            sub_domain,year,month,day
       do i=1,len_trim(fname);if (fname(i:i)==' ')fname(i:i)='0';enddo

       if (my_pe==0) then
         print*,' writing averages to file ',fname(1:len_trim(fname))
         call getunit(io, fname,'u s r ieee')
         write(io) imt,jmt,km,nt
         write(io) xt,yt,xu,yu,zt,zw
         write(io) counter
         write(io) enable_daily_averages
         write(io) enable_monthly_averages
         write(io) enable_seasonal_averages
         write(io) enable_annual_averages
         write(io) kmt_big
       endif
       allocate( buf(imt,jmt) ); buf=0.

       buf(is:ie,js:je)=psimean
       call pe0_recv_2D(buf)
       if (enable_freesurf) then
        if (my_pe==0) write(io) buf,'ssh'
       else
        if (my_pe==0) write(io) buf,'psi'
       endif

       do n=1,2
        buf(is:ie,js:je)=smfmean(:,:,n)
        call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
        do k=1,km
         do j=js,je
          do i=is,ie
            buf(i,j)=umean(i,k,j,n)
          enddo
         enddo
         call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
        enddo
       enddo
       do n=1,nt
        buf(is:ie,js:je)=stfmean(:,:,n)
        call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
        do k=1,km
         do j=js,je
          do i=is,ie
           buf(i,j)=tmean(i,k,j,n)
          enddo
         enddo
         call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
        enddo
      enddo
#ifdef write_vert_diff
      ident = 'diff_cbt'
      if (my_pe==0) write(io) ident
      do k=1,km
         do j=js,je
          do i=is,ie
            buf(i,j)=ka_mean(i,k,j)
          enddo
         enddo
         call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
      enddo
#endif
      if (gm_active) then
       ident = 'u_gm'
       if (my_pe==0) write(io) ident
       do k=1,km
         buf(is:ie,js:je)=ugm_mean(is:ie,k,js:je,1)
         call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
         buf(is:ie,js:je)=ugm_mean(is:ie,k,js:je,2)
         call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
         buf(is:ie,js:je)=ugm_mean(is:ie,k,js:je,3)
         call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
       enddo
      endif
      if (pco2_active) then
       ident = 'pco2_water'
       if (my_pe==0) write(io) ident
       buf(is:ie,js:je)=pco2_mean(is:ie,js:je)
       call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
       buf(is:ie,js:je)=dco2_mean(is:ie,js:je)
       call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
                                    if (my_pe==0) write(io) airco2_mean
      endif
      if (o2_active) then
       ident = 'o2_sat'
       if (my_pe==0) write(io) ident
       buf(is:ie,js:je)=o2_mean(is:ie,js:je)
       call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
      endif

      if (my_pe==0)  close(io)
      deallocate(buf)
      if (my_pe==0) then
        print*,'done'
      endif
#endif
      end subroutine write_averages_cdf





#if defined netcdf_diagnostics && defined splitup_netcdf_file
       subroutine write_averages_tracer
       use spflame_module
       use averages_module
       implicit none
c
c      write averaged quantities to a nice NetCDF file
c      here the tracer stuff
c
       character (len=80) :: fname,name,unit
       integer :: year,month,day
#include "netcdf.inc"
       integer :: ncid,iret,i,j,k,dims(4),start(4),count(4)
       real    :: spval=-9.9e12, tt
       integer :: lon_tdim,lon_tid,lon_udim,lon_uid
       integer :: lat_tdim,lat_tid,lat_udim,lat_uid
       integer :: depth_tdim,depth_tid,depth_wdim,depth_wid
       integer :: itimedim,itimeid,trid(nt),npe,n
       integer :: strid(nt)
       integer :: pco2_id,o2_id,dco2_id,airco2_id
       integer :: is,ie,js,je
#ifdef netcdf_real4
      real (kind=4) :: var(is_pe:ie_pe,js_pe:je_pe)
      real (kind=4) , allocatable :: v2(:)
#else
      real          :: var(is_pe:ie_pe,js_pe:je_pe)
      real          , allocatable :: v2(:)
#endif

       is=is_pe; ie=ie_pe; js=js_pe; je=je_pe

       year  = get_year(get_current_time())
       month = get_month_of_year(get_current_time())
       day   = get_day_of_month(get_current_time())
       write(fname,'("traverages_",i3,"_y",i4,"m",i2,"d",i2,".cdf")') 
     &            sub_domain,year,month,day
       do i=1,len_trim(fname);if (fname(i:i)==' ')fname(i:i)='0';enddo

       if (my_pe == 0) then
        print*,' Writing averages to  NetCDF output file ',
     &         fname(1:len_trim(fname))
        ncid = nccre (fname, NCCLOB, iret)
        iret=nf_set_fill(ncid, NF_NOFILL, iret)
        call store_info_cdf(ncid)
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
        do n=3,nt
         write(name,'("tracer_",i2)') n
         call replace_space_zero(name)
         dims=(/lon_tdim,lat_tdim,depth_tdim,itimedim/)
         trid(n)=ncvdef (ncid,name, NCFLOAT,4,dims,iret)
         write(name,'("tr",i2,"_flux")') n
         call replace_space_zero(name)
         dims=(/lon_tdim,lat_tdim,itimedim,1/)
         strid(n)=ncvdef (ncid,name, NCFLOAT,3,dims,iret)
        enddo
        if (pco2_active) 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 (o2_active) then
         dims=(/lon_tdim,lat_tdim,itimedim,1/)
         o2_id=ncvdef (ncid,'o2_sat', NCFLOAT,3,dims,iret)
        endif
c       attributes of variables
        do n=3,nt
         write(name,'("Passive tracer ",i2)') n
         unit = ' '
         call dvcdf(ncid,trid(n),name,len_trim(name),
     &              unit,len_trim(unit),spval)
         write(name,'("Surface flux passive tracer ",i2)') n
         unit = ' '
         call dvcdf(ncid,strid(n),name,len_trim(name),
     &              unit,len_trim(unit),spval)
        enddo
        if (pco2_active) then
         name = 'Sea surface CO2 partial pressure'; unit = 'ppmv'
         call dvcdf(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(ncid,dco2_id,name,
     &            len_trim(name),unit,len_trim(unit),spval)
         name = 'Atmospheric CO2 partial pressure'; unit = 'ppmv'
         call dvcdf(ncid,airco2_id,name,
     &            len_trim(name),unit,len_trim(unit),spval)
        endif
        if (o2_active) then
         name = 'Sea surface O2 saturation concentration'; 
         unit = 'mol/m^3'
         call dvcdf(ncid,o2_id,name,
     &            len_trim(name),unit,len_trim(unit),spval)
        endif
c       global attrubutes
        iret= nf_put_att_int(ncid,nf_global,
     &   'Number_of_averaged_time_steps ',nf_int,1,counter)
        call ncendf(ncid, iret)
c       write grid
        allocate( v2(max(imt,jmt,km)) )
        v2(1:imt)=xt
        call ncvpt(ncid, Lon_tid, 1, imt,v2, iret)
        v2(1:imt)=xu
        call ncvpt(ncid, Lon_uid, 1, imt,v2, iret)
        v2(1:jmt)=yt
        call ncvpt(ncid, Lat_tid, 1, jmt,v2, iret)
        v2(1:jmt)=yu
        call ncvpt(ncid, Lat_uid, 1, jmt,v2, iret)
        v2(1:km)=zt/100.
        call ncvpt(ncid, depth_tid, 1, km,v2, iret)
        v2(1:km)=zw/100.
        call ncvpt(ncid, depth_wid, 1, km,v2, iret)
        call read_stamp(current_stamp,tt ); v2(1)=tt
        call ncvpt(ncid, itimeid, 1, 1,v2, iret)
        deallocate(v2)
        call ncclos (ncid, iret)
       endif
c      loop over PEs
       do npe=0,n_pes
        call barrier
        if (my_pe==npe) then
         iret=nf_open(fname,NF_WRITE,ncid)
         iret=nf_set_fill(ncid, NF_NOFILL, iret)
         do n=3,nt
          write(name,'("tracer_",i2)') n
          call replace_space_zero(name)
          iret=nf_inq_varid(ncid,name,trid(n))
          write(name,'("tr",i2,"_flux")') n
          call replace_space_zero(name)
          iret=nf_inq_varid(ncid,name,strid(n))
         enddo
         if (pco2_active) then
          iret=nf_inq_varid(ncid,'pco2_water',pco2_id)
          iret=nf_inq_varid(ncid,'dco2',dco2_id)
          iret=nf_inq_varid(ncid,'pco2_air',airco2_id)
         endif
         if (o2_active) iret=nf_inq_varid(ncid,'o2_sat',o2_id)
         ! write PE domain to file
         do k=1,km
          do n=3,nt
           var=tmean(:,k,:,n)
           start=(/is,js,k,1/); count=(/ie-is+1,je-js+1,1,1/)
           iret= nf_put_vara_real (ncid,trid(n),start,count,var )
          enddo
         enddo
         start=(/is,js,1,1/); count=(/ie-is+1,je-js+1,1,1/)
         do n=3,nt
          var=stfmean(:,:,n)
          iret= nf_put_vara_real (ncid,strid(n),start,count,var )
         enddo
         if (pco2_active) then
          var(is:ie,js:je)=pco2_mean(is:ie,js:je)
          start=(/is,js,1,1/); count=(/ie-is+1,je-js+1,1,1/)
          iret= nf_put_vara_real (ncid,pco2_id,start,count,var )
          var(is:ie,js:je)=dco2_mean(is:ie,js:je)
          iret= nf_put_vara_real (ncid,dco2_id,start,count,var )
          var(is,js)=airco2_mean
          iret= nf_put_vara_real(ncid,airco2_id,1,1,var(is,js) )
         endif
         if (o2_active) then
          var=o2_mean(:,:)
          start=(/is,js,1,1/); count=(/ie-is+1,je-js+1,1,1/)
          iret= nf_put_vara_real (ncid,o2_id,start,count,var )
         endif
         call ncclos (ncid, iret)
        endif
        call barrier
       enddo
       end subroutine write_averages_tracer
#endif
