#include "options.inc"

c=======================================================================
c       time averaging diagnostics
c       linked in the code in driver only
c       Reference buoyancy profile given by N_0 
c=======================================================================

#ifdef enable_diag_timeave

c=======================================================================
c local preprocessor switches
c=======================================================================

c#define diag_eddy_fluxes
c#define diag_eddy_flux_tendency
c#define diag_potential_vorticity

c#define remove_back_stratification

      module time_averages_module
c=======================================================================
c        Main module for time averages
c=======================================================================
      implicit none
      integer :: nitts
      real, allocatable :: umean(:,:,:),vmean(:,:,:),wmean(:,:,:)
      real, allocatable :: bmean(:,:,:),pmean(:,:,:)
      real, allocatable :: bs(:,:,:),ps(:,:,:),b3(:,:,:)
      real, allocatable :: uvar(:,:,:),vvar(:,:,:),wvar(:,:,:)
      real, allocatable :: bvar(:,:,:), pvar(:,:,:), uv(:,:,:,:)

#ifdef diag_eddy_fluxes
      real, allocatable :: ub(:,:,:,:),vb(:,:,:,:),wb(:,:,:,:)
      real, allocatable :: ubb(:,:,:,:),vbb(:,:,:,:),wbb(:,:,:,:)
      real, allocatable :: ub3(:,:,:,:),vb3(:,:,:,:),wb3(:,:,:,:)
#endif

#ifdef diag_eddy_flux_tendency
      real, allocatable :: ugeo1(:,:,:),ugeo2(:,:,:)
      real, allocatable :: vgeo1(:,:,:),vgeo2(:,:,:)
      real, allocatable :: bpx(:,:,:,:),bpy(:,:,:,:)
      real, allocatable :: bfv(:,:,:,:),bfu(:,:,:,:)
      real, allocatable :: uwN(:,:,:,:),vwN(:,:,:,:)
#endif
      integer :: nr_tr  = 0  ! numbers of tracers to be considered
      real, allocatable :: tmean(:,:,:,:),tvar(:,:,:,:),t3(:,:,:,:)
#ifdef diag_eddy_fluxes
      real, allocatable :: ut(:,:,:,:,:),vt(:,:,:,:,:),wt(:,:,:,:,:)
      real, allocatable :: utt(:,:,:,:,:),vtt(:,:,:,:,:),wtt(:,:,:,:,:)
      real, allocatable :: ut3(:,:,:,:,:),vt3(:,:,:,:,:),wt3(:,:,:,:,:)
#endif
#ifdef diag_potential_vorticity
      real, allocatable :: rv(:,:,:),pv(:,:,:)
      real, allocatable :: qmean(:,:,:),qvar(:,:,:),q3(:,:,:)
      real, allocatable :: uq(:,:,:,:),vq(:,:,:,:),wq(:,:,:,:)
      real, allocatable :: uqq(:,:,:,:),vqq(:,:,:,:),wqq(:,:,:,:)
      real, allocatable :: uq3(:,:,:,:),vq3(:,:,:,:),wq3(:,:,:,:)
#endif
      real,parameter :: epsln     = 1.0e-20
      end module time_averages_module


      subroutine init_time_averages
c=======================================================================
c      Initialize time mean diagnostics
c=======================================================================
      use cpflame_module
      use time_averages_module
#ifdef enable_diag_tracer
      use tracer_module
#endif
      implicit none
      logical :: file_exists = .false.
      integer :: io,j,n
#ifdef enable_diag_tracer
      nr_tr = nt
#endif
      if (my_pe==0) then
       print*,''
       print*,' Initializing time averaging module '
       print*,''
       print*,' found ',nr_tr,' tracer '
       print*,''
      endif
      allocate(bs(imt,jmt,km),ps(imt,jmt,km) );bs=0.;ps=0.
      allocate(umean(imt,jmt,km),vmean(imt,jmt,km),wmean(imt,jmt,km))
      allocate(bmean(imt,jmt,km),pmean(imt,jmt,km))
      allocate(uvar(imt,jmt,km),vvar(imt,jmt,km),wvar(imt,jmt,km))
      allocate(bvar(imt,jmt,km),pvar(imt,jmt,km),uv(imt,jmt,km,4))
      allocate(b3(imt,jmt,km))
#ifdef diag_eddy_fluxes
      allocate(ub(imt,jmt,km,2),vb(imt,jmt,km,2))
      allocate(ubb(imt,jmt,km,2),vbb(imt,jmt,km,2))
      allocate(wb(imt,jmt,km,2),wbb(imt,jmt,km,2))
      allocate(ub3(imt,jmt,km,2),vb3(imt,jmt,km,2))
      allocate(wb3(imt,jmt,km,2))
#endif
#ifdef diag_eddy_flux_tendency
      allocate(bpx(imt,jmt,km,3),bpy(imt,jmt,km,3))
      allocate(bfv(imt,jmt,km,6),bfu(imt,jmt,km,6))
      allocate(ugeo1(imt,jmt,km),ugeo2(imt,jmt,km) );ugeo1=0.;ugeo2=0.
      allocate(vgeo1(imt,jmt,km),vgeo2(imt,jmt,km) );vgeo1=0.;vgeo2=0.
      allocate(uwN(imt,jmt,km,4),vwN(imt,jmt,km,4))
#endif
      allocate(tmean(imt,jmt,km,nr_tr),tvar(imt,jmt,km,nr_tr))
      allocate(t3(imt,jmt,km,nr_tr) )
#ifdef diag_eddy_fluxes
      allocate(ut(imt,jmt,km,2,nr_tr),vt(imt,jmt,km,2,nr_tr))
      allocate(utt(imt,jmt,km,2,nr_tr),vtt(imt,jmt,km,2,nr_tr))
      allocate(wt(imt,jmt,km,2,nr_tr),wtt(imt,jmt,km,2,nr_tr))
      allocate(ut3(imt,jmt,km,2,nr_tr),vt3(imt,jmt,km,2,nr_tr))
      allocate(wt3(imt,jmt,km,2,nr_tr))
#endif
#ifdef diag_potential_vorticity
      allocate(rv(imt,jmt,km),pv(imt,jmt,km))
      allocate(qmean(imt,jmt,km),qvar(imt,jmt,km),q3(imt,jmt,km))
      allocate(uq(imt,jmt,km,2),vq(imt,jmt,km,2))
      allocate(uqq(imt,jmt,km,2),vqq(imt,jmt,km,2))
      allocate(wq(imt,jmt,km,2),wqq(imt,jmt,km,2))
      allocate(uq3(imt,jmt,km,2),vq3(imt,jmt,km,2))
      allocate(wq3(imt,jmt,km,2))
#endif

      inquire( FILE='restart_ave.dta', EXIST=file_exists )
      if (file_exists) then 
       if (my_pe==0) print*,' reading data from restart file '
       call get_free_iounit(io)
       open(io,file='restart_ave.dta',form='unformatted',status='old')
       read(io) nitts 
       if (my_pe==0) then 
        print*,' number of itts = ',nitts
       endif
       do j=1,jmt
        read(io) umean(:,j,:), vmean(:,j,:), wmean(:,j,:)
        read(io) bmean(:,j,:), pmean(:,j,:)
        read(io) uvar(:,j,:), vvar(:,j,:), wvar(:,j,:)
        read(io) bvar(:,j,:), pvar(:,j,:),uv(:,j,:,:),b3(:,j,:)
#ifdef diag_eddy_fluxes
        read(io) ub(:,j,:,:), vb(:,j,:,:), wb(:,j,:,:)
        read(io) ubb(:,j,:,:),vbb(:,j,:,:),wbb(:,j,:,:)
        read(io) ub3(:,j,:,:),vb3(:,j,:,:),wb3(:,j,:,:)
#endif
#ifdef diag_eddy_flux_tendency
        read(io) bpx(:,j,:,:),bpy(:,j,:,:)
        read(io) bfv(:,j,:,:),bfu(:,j,:,:)
        read(io) uwN(:,j,:,:),vwN(:,j,:,:)
#endif
        do n=1,nr_tr
         read(io) tmean(:,j,:,n), tvar(:,j,:,n),t3(:,j,:,n)
#ifdef diag_eddy_fluxes
         read(io) ut(:,j,:,:,n), vt(:,j,:,:,n), wt(:,j,:,:,n)
         read(io) utt(:,j,:,:,n),vtt(:,j,:,:,n),wtt(:,j,:,:,n)
         read(io) ut3(:,j,:,:,n),vt3(:,j,:,:,n),wt3(:,j,:,:,n)
#endif
       enddo
#ifdef diag_potential_vorticity
        read(io) qmean(:,j,:), qvar(:,j,:),q3(:,j,:)
        read(io) uq(:,j,:,:), vq(:,j,:,:), wq(:,j,:,:)
        read(io) uqq(:,j,:,:),vqq(:,j,:,:),wqq(:,j,:,:)
        read(io) uq3(:,j,:,:),vq3(:,j,:,:),wq3(:,j,:,:)
#endif
       enddo
       close(io)
       if (my_pe==0) print*,' done reading restart'
      else
       if (my_pe==0) print*,' cannot read restart file '
       nitts=0
       umean=0;vmean=0;wmean=0;bmean=0;pmean=0
       uvar=0;vvar=0;wvar=0;bvar=0;pvar=0; uv=0.;b3=0
#ifdef diag_eddy_fluxes
       ub=0;vb=0;wb=0; ubb=0;vbb=0;wbb=0
       ub3=0;vb3=0;wb3=0.
#endif
#ifdef diag_eddy_flux_tendency
       bpx=0;bpy=0;bfv=0;bfu=0;uwN=0;vwN=0
#endif
       do n=1,nr_tr
        tmean(:,:,:,n)=0;tvar(:,:,:,n)=0;t3(:,:,:,n)=0
#ifdef diag_eddy_fluxes
        ut(:,:,:,:,n)=0;vt(:,:,:,:,n)=0;wt(:,:,:,:,n)=0; 
        utt(:,:,:,:,n)=0;vtt(:,:,:,:,n)=0;wtt(:,:,:,:,n)=0
        ut3(:,:,:,:,n)=0;vt3(:,:,:,:,n)=0;wt3(:,:,:,:,n)=0
#endif
       enddo
#ifdef diag_potential_vorticity
       qmean=0;qvar=0;q3=0
       uq=0;vq=0;wq=0;wqq=0;uqq=0;vqq=0
       vq3=0;uq3=0;wq3=0
#endif
      endif
      if (my_pe==0) print*,' Done initializing time averaging module'
      end subroutine init_time_averages


      subroutine time_averages_restart
c=======================================================================
c       write restart
c=======================================================================
      use cpflame_module
      use time_averages_module
      implicit none
      integer :: io,n,js,je,j,k

      if (my_pe==0) then
       print*,' writing restart file for time averages'
       call get_free_iounit(io)
       open(io,file='restart_ave.dta',
     &         form='unformatted',status='unknown')
       write(io) nitts
       do j=js_pe,je_pe
        write(io) umean(:,j,:),vmean(:,j,:),wmean(:,j,:)
        write(io) bmean(:,j,:),pmean(:,j,:)
        write(io) uvar(:,j,:),vvar(:,j,:),wvar(:,j,:)
        write(io) bvar(:,j,:),pvar(:,j,:),uv(:,j,:,:),b3(:,j,:)
#ifdef diag_eddy_fluxes
        write(io)  ub(:,j,:,:), vb(:,j,:,:), wb(:,j,:,:)
        write(io) ubb(:,j,:,:),vbb(:,j,:,:),wbb(:,j,:,:)
        write(io) ub3(:,j,:,:),vb3(:,j,:,:),wb3(:,j,:,:)
#endif
#ifdef diag_eddy_flux_tendency
        write(io) bpx(:,j,:,:),bpy(:,j,:,:)
        write(io) bfv(:,j,:,:),bfu(:,j,:,:)
        write(io) uwN(:,j,:,:),vwN(:,j,:,:)
#endif
        do k=1,nr_tr
         write(io) tmean(:,j,:,k),tvar(:,j,:,k),t3(:,j,:,k)
#ifdef diag_eddy_fluxes
         write(io)  ut(:,j,:,:,k), vt(:,j,:,:,k), wt(:,j,:,:,k)
         write(io) utt(:,j,:,:,k),vtt(:,j,:,:,k),wtt(:,j,:,:,k)
         write(io) ut3(:,j,:,:,k),vt3(:,j,:,:,k),wt3(:,j,:,:,k)
#endif
        enddo
#ifdef diag_potential_vorticity
        write(io)  qmean(:,j,:), qvar(:,j,:), q3(:,j,:)
        write(io)  uq(:,j,:,:), vq(:,j,:,:), wq(:,j,:,:)
        write(io) uqq(:,j,:,:),vqq(:,j,:,:),wqq(:,j,:,:)
        write(io) uq3(:,j,:,:),vq3(:,j,:,:),wq3(:,j,:,:)
#endif
       enddo
      endif

      do n=1,n_pes-1
        if (my_pe==n) then
         call send_integer(js_pe,1,0,0)
         call send_integer(je_pe,1,0,0)
         do j=js_pe,je_pe
          call send_real(umean(:,j,:),imt*km,0,0)
          call send_real(vmean(:,j,:),imt*km,0,0)
          call send_real(wmean(:,j,:),imt*km,0,0)
          call send_real(bmean(:,j,:),imt*km,0,0)
          call send_real(pmean(:,j,:),imt*km,0,0)
          call send_real(uvar(:,j,:),imt*km,0,0)
          call send_real(vvar(:,j,:),imt*km,0,0)
          call send_real(wvar(:,j,:),imt*km,0,0)
          call send_real(bvar(:,j,:),imt*km,0,0)
          call send_real(pvar(:,j,:),imt*km,0,0)
          call send_real(uv(:,j,:,:),imt*km*4,0,0)
          call send_real(b3(:,j,:),imt*km,0,0)
#ifdef diag_eddy_fluxes
          call send_real(ub(:,j,:,:),imt*km*2,0,0)
          call send_real(vb(:,j,:,:),imt*km*2,0,0)
          call send_real(wb(:,j,:,:),imt*km*2,0,0)
          call send_real(ubb(:,j,:,:),imt*km*2,0,0)
          call send_real(vbb(:,j,:,:),imt*km*2,0,0)
          call send_real(wbb(:,j,:,:),imt*km*2,0,0)
          call send_real(ub3(:,j,:,:),imt*km*2,0,0)
          call send_real(vb3(:,j,:,:),imt*km*2,0,0)
          call send_real(wb3(:,j,:,:),imt*km*2,0,0)
#endif
#ifdef diag_eddy_flux_tendency
          call send_real(bpx(:,j,:,:),imt*km*3,0,0)
          call send_real(bpy(:,j,:,:),imt*km*3,0,0)
          call send_real(bfv(:,j,:,:),imt*km*6,0,0)
          call send_real(bfu(:,j,:,:),imt*km*6,0,0)
          call send_real(uwN(:,j,:,:),imt*km*4,0,0)
          call send_real(vwN(:,j,:,:),imt*km*4,0,0)
#endif
          do k=1,nr_tr
           call send_real(tmean(:,j,:,k),imt*km,0,0)
           call send_real(tvar(:,j,:,k),imt*km,0,0)
           call send_real(t3(:,j,:,k),imt*km,0,0)
#ifdef diag_eddy_fluxes
           call send_real(ut(:,j,:,:,k),imt*km*2,0,0)
           call send_real(vt(:,j,:,:,k),imt*km*2,0,0)
           call send_real(wt(:,j,:,:,k),imt*km*2,0,0)
           call send_real(utt(:,j,:,:,k),imt*km*2,0,0)
           call send_real(vtt(:,j,:,:,k),imt*km*2,0,0)
           call send_real(wtt(:,j,:,:,k),imt*km*2,0,0)
           call send_real(ut3(:,j,:,:,k),imt*km*2,0,0)
           call send_real(vt3(:,j,:,:,k),imt*km*2,0,0)
           call send_real(wt3(:,j,:,:,k),imt*km*2,0,0)
#endif
          enddo
#ifdef diag_potential_vorticity
          call send_real(qmean(:,j,:),imt*km,0,0)
          call send_real(qvar(:,j,:),imt*km,0,0)
          call send_real(q3(:,j,:),imt*km,0,0)
          call send_real(uq(:,j,:,:),imt*km*2,0,0)
          call send_real(vq(:,j,:,:),imt*km*2,0,0)
          call send_real(wq(:,j,:,:),imt*km*2,0,0)
          call send_real(uqq(:,j,:,:),imt*km*2,0,0)
          call send_real(vqq(:,j,:,:),imt*km*2,0,0)
          call send_real(wqq(:,j,:,:),imt*km*2,0,0)
          call send_real(uq3(:,j,:,:),imt*km*2,0,0)
          call send_real(vq3(:,j,:,:),imt*km*2,0,0)
          call send_real(wq3(:,j,:,:),imt*km*2,0,0)
#endif
         enddo
        elseif (my_pe==0) then
         call recv_integer(js,1,n,0)
         call recv_integer(je,1,n,0)
         do j=js,je
          call recv_real(umean(:,j,:),imt*km,n,0)
          call recv_real(vmean(:,j,:),imt*km,n,0)
          call recv_real(wmean(:,j,:),imt*km,n,0)
          call recv_real(bmean(:,j,:),imt*km,n,0)
          call recv_real(pmean(:,j,:),imt*km,n,0)
          call recv_real(uvar(:,j,:),imt*km,n,0)
          call recv_real(vvar(:,j,:),imt*km,n,0)
          call recv_real(wvar(:,j,:),imt*km,n,0)
          call recv_real(bvar(:,j,:),imt*km,n,0)
          call recv_real(pvar(:,j,:),imt*km,n,0)
          call recv_real(uv(:,j,:,:),imt*km*4,n,0)
          call recv_real(b3(:,j,:),imt*km,n,0)
#ifdef diag_eddy_fluxes
          call recv_real(ub(:,j,:,:),imt*km*2,n,0)
          call recv_real(vb(:,j,:,:),imt*km*2,n,0)
          call recv_real(wb(:,j,:,:),imt*km*2,n,0)
          call recv_real(ubb(:,j,:,:),imt*km*2,n,0)
          call recv_real(vbb(:,j,:,:),imt*km*2,n,0)
          call recv_real(wbb(:,j,:,:),imt*km*2,n,0)
          call recv_real(ub3(:,j,:,:),imt*km*2,n,0)
          call recv_real(vb3(:,j,:,:),imt*km*2,n,0)
          call recv_real(wb3(:,j,:,:),imt*km*2,n,0)
#endif
#ifdef diag_eddy_flux_tendency
          call recv_real(bpx(:,j,:,:),imt*km*3,n,0)
          call recv_real(bpy(:,j,:,:),imt*km*3,n,0)
          call recv_real(bfv(:,j,:,:),imt*km*6,n,0)
          call recv_real(bfu(:,j,:,:),imt*km*6,n,0)
          call recv_real(uwN(:,j,:,:),imt*km*4,n,0)
          call recv_real(vwN(:,j,:,:),imt*km*4,n,0)
#endif
          do k=1,nr_tr
           call recv_real(tmean(:,j,:,k),imt*km,n,0)
           call recv_real(tvar(:,j,:,k),imt*km,n,0)
           call recv_real(t3(:,j,:,k),imt*km,n,0)
#ifdef diag_eddy_fluxes
           call recv_real(ut(:,j,:,:,k),imt*km*2,n,0)
           call recv_real(vt(:,j,:,:,k),imt*km*2,n,0)
           call recv_real(wt(:,j,:,:,k),imt*km*2,n,0)
           call recv_real(utt(:,j,:,:,k),imt*km*2,n,0)
           call recv_real(vtt(:,j,:,:,k),imt*km*2,n,0)
           call recv_real(wtt(:,j,:,:,k),imt*km*2,n,0)
           call recv_real(ut3(:,j,:,:,k),imt*km*2,n,0)
           call recv_real(vt3(:,j,:,:,k),imt*km*2,n,0)
           call recv_real(wt3(:,j,:,:,k),imt*km*2,n,0)
#endif
          enddo
#ifdef diag_potential_vorticity
          call recv_real(qmean(:,j,:),imt*km,n,0)
          call recv_real(qvar(:,j,:),imt*km,n,0)
          call recv_real(q3(:,j,:),imt*km,n,0)
          call recv_real(uq(:,j,:,:),imt*km*2,n,0)
          call recv_real(vq(:,j,:,:),imt*km*2,n,0)
          call recv_real(wq(:,j,:,:),imt*km*2,n,0)
          call recv_real(uqq(:,j,:,:),imt*km*2,n,0)
          call recv_real(vqq(:,j,:,:),imt*km*2,n,0)
          call recv_real(wqq(:,j,:,:),imt*km*2,n,0)
          call recv_real(uq3(:,j,:,:),imt*km*2,n,0)
          call recv_real(vq3(:,j,:,:),imt*km*2,n,0)
          call recv_real(wq3(:,j,:,:),imt*km*2,n,0)
#endif
          write(io) umean(:,j,:),vmean(:,j,:),wmean(:,j,:)
          write(io) bmean(:,j,:),pmean(:,j,:)
          write(io) uvar(:,j,:),vvar(:,j,:),wvar(:,j,:)
          write(io) bvar(:,j,:),pvar(:,j,:),uv(:,j,:,:),b3(:,j,:)
#ifdef diag_eddy_fluxes
          write(io)  ub(:,j,:,:), vb(:,j,:,:), wb(:,j,:,:)
          write(io) ubb(:,j,:,:),vbb(:,j,:,:),wbb(:,j,:,:)
          write(io) ub3(:,j,:,:),vb3(:,j,:,:),wb3(:,j,:,:)
#endif
#ifdef diag_eddy_flux_tendency
          write(io) bpx(:,j,:,:),bpy(:,j,:,:)
          write(io) bfv(:,j,:,:),bfu(:,j,:,:)
          write(io) uwN(:,j,:,:),vwN(:,j,:,:)
#endif
          do k=1,nr_tr
           write(io) tmean(:,j,:,k),tvar(:,j,:,k),t3(:,j,:,k)
#ifdef diag_eddy_fluxes
           write(io)  ut(:,j,:,:,k), vt(:,j,:,:,k), wt(:,j,:,:,k)
           write(io) utt(:,j,:,:,k),vtt(:,j,:,:,k),wtt(:,j,:,:,k)
           write(io) ut3(:,j,:,:,k),vt3(:,j,:,:,k),wt3(:,j,:,:,k)
#endif
          enddo
#ifdef diag_potential_vorticity
          write(io)  qmean(:,j,:), qvar(:,j,:), q3(:,j,:)
          write(io)  uq(:,j,:,:), vq(:,j,:,:), wq(:,j,:,:)
          write(io) uqq(:,j,:,:),vqq(:,j,:,:),wqq(:,j,:,:)
          write(io) uq3(:,j,:,:),vq3(:,j,:,:),wq3(:,j,:,:)
#endif
         enddo
        endif
        call barrier
      enddo
      if (my_pe==0) then
        close(io)
        print*,' done '
      endif
      end subroutine time_averages_restart


      subroutine time_averages
c=======================================================================
c        Hauptverteiler function fuer time averages
c=======================================================================
      use cpflame_module
      use time_averages_module
      implicit none
      integer :: i,j,k,js,je
      real :: fxa
c
c        first calculate auxilliary variables
c
      js=max(2,js_pe); je = min(je_pe,jmt-1)
      do k=2,km-1
       do j=js,je 
        do i=2,imt-1
#ifdef  remove_back_stratification
         bs(i,j,k)=b(i,j,k,tau)-b_r(k)
         ps(i,j,k)=p_full(i,j,k,tau)-p_r(k)
#else
         bs(i,j,k)=b(i,j,k,tau)
         ps(i,j,k)=p_full(i,j,k,tau)
#endif
        enddo
       enddo
      enddo
      call border_exchg3D(bs,1); call setcyclic3D(bs)
      call border_exchg3D(ps,1); call setcyclic3D(ps)
#ifdef diag_potential_vorticity
c
c       relative potential vorticity
c
      rv=0.0
      do j=js,je 
       do k=2,km-1
        do i=2,imt-1
          rv(i,j,k)=
     &       (u(i+1,j,k,2,tau)-u(i,j,k,2,tau) )/dx*maskU(i,j,k)
     &      -(u(i,j+1,k,1,tau)-u(i,j,k,1,tau) )/dx*maskV(i,j,k)
        enddo 
       enddo 
      enddo 
      call border_exchg3D(rv,1); call setcyclic3D(rv)
c
c       Ertel potential vorticity  Q N0^2 = b_z (v_x-u_y + f) + u_z b_y - v_z b_x
c
      pv=0.
      do j=js,je 
       do k=2,km-1
        do i=2,imt-1
         fxa= (rv(i,j,k)+rv(i-1,j,k)+rv(i,j-1,k)+rv(i-1,j-1,k))/4.0
         pv(i,j,k)=(fxa+coriolis_t(j)  )*
     &    ((b(i,j,k+1,tau)-b(i,j,k,tau))/dz*maskW(i,j,k)
     &    +(b(i,j,k,tau)-b(i,j,k-1,tau))/dz*maskW(i,j,k-1))/
     &     (maskW(i,j,k)+maskW(i,j,k-1)+epsln) 
         pv(i,j,k)=-N_0**(-2)*pv(i,j,k)
        enddo 
       enddo 
      enddo 
      call border_exchg3D(pv,1); call setcyclic3D(pv)
#endif
c
c       average
c
      nitts=nitts+1      
      call time_averages_means
      call time_averages_eddy_fluxes
      call time_averages_flux_tendency
c
c        write results and restart file
c
      if (last_time_step) then
       call time_averages_restart
       call time_averages_calc
       call time_averages_write
      endif
      end subroutine time_averages


      subroutine time_averages_means
c=======================================================================
c       sum up the mean variables
c=======================================================================
      use cpflame_module
      use time_averages_module
#ifdef enable_diag_tracer
      use tracer_module
#endif
      implicit none
      integer :: i,j,k,js,je
      js=max(2,js_pe); je = min(je_pe,jmt-1)
      do k=2,km-1
       do j=js,je 
        do i=2,imt-1
         umean(i,j,k)=umean(i,j,k)+u(i,j,k,1,tau)
         vmean(i,j,k)=vmean(i,j,k)+u(i,j,k,2,tau)
         wmean(i,j,k)=wmean(i,j,k)+u(i,j,k,3,tau)
         bmean(i,j,k)=bmean(i,j,k)+bs(i,j,k)
         pmean(i,j,k)=pmean(i,j,k)+ps(i,j,k)
         uvar(i,j,k) =uvar(i,j,k) +u(i,j,k,1,tau)**2
         vvar(i,j,k) =vvar(i,j,k) +u(i,j,k,2,tau)**2
         wvar(i,j,k) =wvar(i,j,k) +u(i,j,k,3,tau)**2
         bvar(i,j,k) =bvar(i,j,k) +bs(i,j,k)**2
         b3(i,j,k)    =b3(i,j,k) +bs(i,j,k)**3
         pvar(i,j,k) =pvar(i,j,k) +ps(i,j,k)**2
         uv(i,j,k,1) =uv(i,j,k,1) +u(i  ,j,k,1,tau)*u(i,j  ,k,2,tau)
         uv(i,j,k,2) =uv(i,j,k,2) +u(i-1,j,k,1,tau)*u(i,j  ,k,2,tau)
         uv(i,j,k,3) =uv(i,j,k,3) +u(i  ,j,k,1,tau)*u(i,j-1,k,2,tau)
         uv(i,j,k,4) =uv(i,j,k,4) +u(i-1,j,k,1,tau)*u(i,j-1,k,2,tau)
#ifdef enable_diag_tracer
         tmean(i,j,k,:)=tmean(i,j,k,:)+tr(i,j,k,tau,:)
         tvar(i,j,k,:)=tvar(i,j,k,:)+tr(i,j,k,tau,:)**2
         t3(i,j,k,:)=t3(i,j,k,:)+tr(i,j,k,tau,:)**3
#endif
#ifdef diag_potential_vorticity
         qmean(i,j,k)=qmean(i,j,k)+pv(i,j,k)
         qvar(i,j,k)=qvar(i,j,k)+pv(i,j,k)**2
         q3(i,j,k)=q3(i,j,k)+pv(i,j,k)**3
#endif
        enddo
       enddo
      enddo
      end subroutine time_averages_means


      subroutine time_averages_eddy_fluxes
c=======================================================================
c       sum up the eddy and variance fluxes 
c=======================================================================
      use cpflame_module
      use time_averages_module
#ifdef enable_diag_tracer
      use tracer_module
#endif
      implicit none
      integer :: i,j,k,js,je

#ifdef diag_eddy_fluxes
      js=max(2,js_pe); je = min(je_pe,jmt-1)
      do k=2,km-1
       do j=js,je 
        do i=2,imt-1
c        buoyancy eddy fluxes
          ub(i,j,k,1) =ub(i,j,k,1) +u(i,j,k,1,tau)*bs(i  ,j,k)          
          ub(i,j,k,2) =ub(i,j,k,2) +u(i,j,k,1,tau)*bs(i+1,j,k)          
          vb(i,j,k,1) =vb(i,j,k,1) +u(i,j,k,2,tau)*bs(i,j  ,k)          
          vb(i,j,k,2) =vb(i,j,k,2) +u(i,j,k,2,tau)*bs(i,j+1,k)          
          wb(i,j,k,1) =wb(i,j,k,1) +u(i,j,k,3,tau)*bs(i,j,k  )          
          wb(i,j,k,2) =wb(i,j,k,2) +u(i,j,k,3,tau)*bs(i,j,k+1)          
c        buoyancy variance fluxes
          ubb(i,j,k,1)=ubb(i,j,k,1)+u(i,j,k,1,tau)*bs(i  ,j,k)**2
          ubb(i,j,k,2)=ubb(i,j,k,2)+u(i,j,k,1,tau)*bs(i+1,j,k)**2
          vbb(i,j,k,1)=vbb(i,j,k,1)+u(i,j,k,2,tau)*bs(i,j  ,k)**2
          vbb(i,j,k,2)=vbb(i,j,k,2)+u(i,j,k,2,tau)*bs(i,j+1,k)**2
          wbb(i,j,k,1)=wbb(i,j,k,1)+u(i,j,k,3,tau)*bs(i,j,k  )**2
          wbb(i,j,k,2)=wbb(i,j,k,2)+u(i,j,k,3,tau)*bs(i,j,k+1)**2
c        third order fluxes
          ub3(i,j,k,1)=ub3(i,j,k,1)+u(i,j,k,1,tau)*bs(i  ,j,k)**3
          ub3(i,j,k,2)=ub3(i,j,k,2)+u(i,j,k,1,tau)*bs(i+1,j,k)**3
          vb3(i,j,k,1)=vb3(i,j,k,1)+u(i,j,k,2,tau)*bs(i,j  ,k)**3
          vb3(i,j,k,2)=vb3(i,j,k,2)+u(i,j,k,2,tau)*bs(i,j+1,k)**3
          wb3(i,j,k,1)=wb3(i,j,k,1)+u(i,j,k,3,tau)*bs(i,j,k  )**3
          wb3(i,j,k,2)=wb3(i,j,k,2)+u(i,j,k,3,tau)*bs(i,j,k+1)**3
#ifdef enable_diag_tracer
c        tracer eddy fluxes
          ut(i,j,k,1,:) =ut(i,j,k,1,:) +u(i,j,k,1,tau)*tr(i  ,j,k,tau,:)          
          ut(i,j,k,2,:) =ut(i,j,k,2,:) +u(i,j,k,1,tau)*tr(i+1,j,k,tau,:)          
          vt(i,j,k,1,:) =vt(i,j,k,1,:) +u(i,j,k,2,tau)*tr(i,j  ,k,tau,:)          
          vt(i,j,k,2,:) =vt(i,j,k,2,:) +u(i,j,k,2,tau)*tr(i,j+1,k,tau,:)          
          wt(i,j,k,1,:) =wt(i,j,k,1,:)+u(i,j,k,3,tau)*tr(i,j,k  ,tau,:)          
          wt(i,j,k,2,:) =wt(i,j,k,2,:)+u(i,j,k,3,tau)*tr(i,j,k+1,tau,:)          
c        tracer variance fluxes
          utt(i,j,k,1,:)=utt(i,j,k,1,:)
     &                     +u(i,j,k,1,tau)*tr(i  ,j,k,tau,:)**2
          utt(i,j,k,2,:)=utt(i,j,k,2,:)
     &                     +u(i,j,k,1,tau)*tr(i+1,j,k,tau,:)**2
          vtt(i,j,k,1,:)=vtt(i,j,k,1,:)
     &                      +u(i,j,k,2,tau)*tr(i,j  ,k,tau,:)**2
          vtt(i,j,k,2,:)=vtt(i,j,k,2,:)
     &                      +u(i,j,k,2,tau)*tr(i,j+1,k,tau,:)**2
          wtt(i,j,k,1,:)=wtt(i,j,k,1,:)
     &                      +u(i,j,k,3,tau)*tr(i,j,k  ,tau,:)**2
          wtt(i,j,k,2,:)=wtt(i,j,k,2,:)
     &                       +u(i,j,k,3,tau)*tr(i,j,k+1,tau,:)**2
c        third order fluxes
          ut3(i,j,k,1,:)=ut3(i,j,k,1,:)
     &                     +u(i,j,k,1,tau)*tr(i  ,j,k,tau,:)**3
          ut3(i,j,k,2,:)=ut3(i,j,k,2,:)
     &                     +u(i,j,k,1,tau)*tr(i+1,j,k,tau,:)**3
          vt3(i,j,k,1,:)=vt3(i,j,k,1,:)
     &                      +u(i,j,k,2,tau)*tr(i,j  ,k,tau,:)**3
          vt3(i,j,k,2,:)=vt3(i,j,k,2,:)
     &                      +u(i,j,k,2,tau)*tr(i,j+1,k,tau,:)**3
          wt3(i,j,k,1,:)=wt3(i,j,k,1,:)
     &                      +u(i,j,k,3,tau)*tr(i,j,k  ,tau,:)**3
          wt3(i,j,k,2,:)=wt3(i,j,k,2,:)
     &                       +u(i,j,k,3,tau)*tr(i,j,k+1,tau,:)**3
#endif
#ifdef diag_potential_vorticity
c        pv eddy fluxes
          uq(i,j,k,1) =uq(i,j,k,1) +u(i,j,k,1,tau)*pv(i  ,j,k)          
          uq(i,j,k,2) =uq(i,j,k,2) +u(i,j,k,1,tau)*pv(i+1,j,k)          
          vq(i,j,k,1) =vq(i,j,k,1) +u(i,j,k,2,tau)*pv(i,j  ,k)          
          vq(i,j,k,2) =vq(i,j,k,2) +u(i,j,k,2,tau)*pv(i,j+1,k)          
          wq(i,j,k,1) =wq(i,j,k,1) +u(i,j,k,3,tau)*pv(i,j,k  )          
          wq(i,j,k,2) =wq(i,j,k,2) +u(i,j,k,3,tau)*pv(i,j,k+1)          
c        pv variance fluxes
          uqq(i,j,k,1)=uqq(i,j,k,1)+u(i,j,k,1,tau)*pv(i  ,j,k)**2
          uqq(i,j,k,2)=uqq(i,j,k,2)+u(i,j,k,1,tau)*pv(i+1,j,k)**2
          vqq(i,j,k,1)=vqq(i,j,k,1)+u(i,j,k,2,tau)*pv(i,j  ,k)**2
          vqq(i,j,k,2)=vqq(i,j,k,2)+u(i,j,k,2,tau)*pv(i,j+1,k)**2
          wqq(i,j,k,1)=wqq(i,j,k,1)+u(i,j,k,3,tau)*pv(i,j,k  )**2
          wqq(i,j,k,2)=wqq(i,j,k,2)+u(i,j,k,3,tau)*pv(i,j,k+1)**2
c        third order fluxes
          uq3(i,j,k,1)=uq3(i,j,k,1)+u(i,j,k,1,tau)*pv(i  ,j,k)**3
          uq3(i,j,k,2)=uq3(i,j,k,2)+u(i,j,k,1,tau)*pv(i+1,j,k)**3
          vq3(i,j,k,1)=vq3(i,j,k,1)+u(i,j,k,2,tau)*pv(i,j  ,k)**3
          vq3(i,j,k,2)=vq3(i,j,k,2)+u(i,j,k,2,tau)*pv(i,j+1,k)**3
          wq3(i,j,k,1)=wq3(i,j,k,1)+u(i,j,k,3,tau)*pv(i,j,k  )**3
          wq3(i,j,k,2)=wq3(i,j,k,2)+u(i,j,k,3,tau)*pv(i,j,k+1)**3
#endif
        enddo
       enddo
      enddo
#endif
      end subroutine time_averages_eddy_fluxes



      subroutine time_averages_flux_tendency
c=======================================================================
c        flux tendencies
c=======================================================================
      use cpflame_module
      use time_averages_module
      implicit none
      integer :: i,j,k,js,je
      js=max(2,js_pe); je = min(je_pe,jmt-1)

#ifdef diag_eddy_flux_tendency
      do k=2,km-1
       do j=js,je 
        do i=2,imt-1
         bpx(i,j,k,1)=bpx(i,j,k,1)+bs(i  ,j,k)*ps(i  ,j,k)
         bpx(i,j,k,2)=bpx(i,j,k,2)+bs(i+1,j,k)*ps(i  ,j,k)
         bpx(i,j,k,3)=bpx(i,j,k,3)+bs(i  ,j,k)*ps(i+1,j,k)

         bpy(i,j,k,1)=bpy(i,j,k,1)+bs(i,j  ,k)*ps(i,j  ,k)
         bpy(i,j,k,2)=bpy(i,j,k,2)+bs(i,j+1,k)*ps(i,j  ,k)
         bpy(i,j,k,3)=bpy(i,j,k,3)+bs(i,j  ,k)*ps(i,j+1,k)

         bfv(i,j,k,1)=bfv(i,j,k,1) + bs(i  ,j,k)*u(i  ,j  ,k,2,tau)
         bfv(i,j,k,2)=bfv(i,j,k,2) + bs(i+1,j,k)*u(i  ,j  ,k,2,tau)
         bfv(i,j,k,3)=bfv(i,j,k,3) + bs(i  ,j,k)*u(i+1,j  ,k,2,tau)
         bfv(i,j,k,4)=bfv(i,j,k,4) + bs(i  ,j,k)*u(i  ,j-1,k,2,tau)
         bfv(i,j,k,5)=bfv(i,j,k,5) + bs(i+1,j,k)*u(i  ,j-1,k,2,tau)
         bfv(i,j,k,6)=bfv(i,j,k,6) + bs(i  ,j,k)*u(i+1,j-1,k,2,tau)

         bfu(i,j,k,1)=bfu(i,j,k,1) + bs(i,j  ,k)*u(i  ,j  ,k,1,tau)
         bfu(i,j,k,2)=bfu(i,j,k,2) + bs(i,j+1,k)*u(i  ,j  ,k,1,tau)
         bfu(i,j,k,3)=bfu(i,j,k,3) + bs(i,j  ,k)*u(i  ,j+1,k,1,tau)
         bfu(i,j,k,4)=bfu(i,j,k,4) + bs(i,j  ,k)*u(i-1,j  ,k,1,tau)
         bfu(i,j,k,5)=bfu(i,j,k,5) + bs(i,j+1,k)*u(i-1,j  ,k,1,tau)
         bfu(i,j,k,6)=bfu(i,j,k,6) + bs(i,j  ,k)*u(i-1,j+1,k,1,tau)

         uwN(i,j,k,1)=uwN(i,j,k,1) + u(i,j,k,1,tau)*u(i  ,j,k  ,3,tau)
         uwN(i,j,k,2)=uwN(i,j,k,2) + u(i,j,k,1,tau)*u(i+1,j,k  ,3,tau)
         uwN(i,j,k,3)=uwN(i,j,k,3) + u(i,j,k,1,tau)*u(i  ,j,k-1,3,tau)
         uwN(i,j,k,4)=uwN(i,j,k,4) + u(i,j,k,1,tau)*u(i+1,j,k-1,3,tau)

         vwN(i,j,k,1)=vwN(i,j,k,1) + u(i,j,k,2,tau)*u(i,j  ,k  ,3,tau)
         vwN(i,j,k,2)=vwN(i,j,k,2) + u(i,j,k,2,tau)*u(i,j+1,k  ,3,tau)
         vwN(i,j,k,3)=vwN(i,j,k,3) + u(i,j,k,2,tau)*u(i,j  ,k-1,3,tau)
         vwN(i,j,k,4)=vwN(i,j,k,4) + u(i,j,k,2,tau)*u(i,j+1,k-1,3,tau)
        enddo
       enddo
      enddo
#endif
      end subroutine time_averages_flux_tendency



      subroutine time_averages_calc
c=======================================================================
c        calculate time mean diagnostics
c=======================================================================
      use cpflame_module
      use time_averages_module
#ifdef enable_diag_tracer
      use tracer_module
#endif
      implicit none
      integer :: i,j,k,js,je,n,m
      js=max(2,js_pe); je = min(je_pe,jmt-1)
c
c        calculate means
c
      n=max(1,nitts)
      umean=umean*maskU/n
      vmean=vmean*maskV/n
      wmean=wmean*maskW/n
      bmean=bmean*maskT/n
      pmean=pmean*maskT/n
      call border_exchg3D(bmean,1); call setcyclic3D(bmean)
      call border_exchg3D(pmean,1); call setcyclic3D(pmean)
      call border_exchg3D(umean,1); call setcyclic3D(umean)
      call border_exchg3D(vmean,1); call setcyclic3D(vmean)
      do k=1,nr_tr
       tmean(:,:,:,k)=tmean(:,:,:,k)*maskT/n
       call border_exchg3D(tmean(:,:,:,k),1)
       call setcyclic3D(tmean(:,:,:,k))
      enddo
#ifdef diag_potential_vorticity
      qmean=qmean*maskT/n
      call border_exchg3D(qmean,1); call setcyclic3D(qmean)
#endif
c
c       Variances
c
      uvar=uvar/n - umean**2
      vvar=vvar/n - vmean**2
      wvar=wvar/n - wmean**2
      bvar=bvar/n - bmean**2
      pvar=pvar/n - pmean**2
!  (b+b')(b+b')(b+b')=(b^2 + 2bb'+b'^2)(b+b')=b^3 + 2b^2b'+b'^2 b + b^2b'+2bb'^2+b'^3=b^3 + 3bb'^2+b'^3
      b3=b3/n - bmean**3-3*bmean*bvar
      do k=1,nr_tr
       tvar(:,:,:,k)=tvar(:,:,:,k)/n - tmean(:,:,:,k)**2
       t3(:,:,:,k)=t3(:,:,:,k)/n - tmean(:,:,:,k)**3
     &                               -3*tmean(:,:,:,k)*tvar(:,:,:,k)
      enddo
#ifdef diag_potential_vorticity
      qvar=qvar/n - qmean**2
      q3=q3/n - qmean**3-3*qmean*qvar
#endif

      do k=2,km-1
       do j=js,je 
        do i=2,imt-1
         uv(i,j,k,1)=uv(i,j,k,1)/n +umean(i  ,j,k)*vmean(i,j  ,k)
         uv(i,j,k,2)=uv(i,j,k,2)/n +umean(i-1,j,k)*vmean(i,j  ,k)
         uv(i,j,k,3)=uv(i,j,k,3)/n +umean(i  ,j,k)*vmean(i,j-1,k)
         uv(i,j,k,4)=uv(i,j,k,4)/n +umean(i-1,j,k)*vmean(i,j-1,k)
         uv(i,j,k,1)=(uv(i,j,k,1)+uv(i,j,k,2)
     &               +uv(i,j,k,3)+uv(i,j,k,4))/4.0
        enddo
       enddo
      enddo

#ifdef diag_eddy_fluxes
c
c       Eddy fluxes
c
      do k=2,km-1
       do j=js,je 
        do i=2,imt-1
c       buoyancy
         ub(i,j,k,1)=ub(i,j,k,1)/n-umean(i,j,k)*bmean(i ,j,k)
         ub(i,j,k,2)=ub(i,j,k,2)/n-umean(i,j,k)*bmean(i+1,j,k)
         vb(i,j,k,1)=vb(i,j,k,1)/n-vmean(i,j,k)*bmean(i,j,k)
         vb(i,j,k,2)=vb(i,j,k,2)/n-vmean(i,j,k)*bmean(i,j+1,k)
         wb(i,j,k,1)=wb(i,j,k,1)/n-wmean(i,j,k)*bmean(i,j,k)
         wb(i,j,k,2)=wb(i,j,k,2)/n-wmean(i,j,k)*bmean(i,j,k+1)
        enddo
       enddo
      enddo

      do m=1,nr_tr
       do k=2,km-1
        do j=js,je 
         do i=2,imt-1
c       tracer
          ut(i,j,k,1,m)=ut(i,j,k,1,m)/n-umean(i,j,k)*tmean(i,j,k,m)
          ut(i,j,k,2,m)=ut(i,j,k,2,m)/n-umean(i,j,k)*tmean(i+1,j,k,m)
          vt(i,j,k,1,m)=vt(i,j,k,1,m)/n-vmean(i,j,k)*tmean(i,j,k,m)
          vt(i,j,k,2,m)=vt(i,j,k,2,m)/n-vmean(i,j,k)*tmean(i,j+1,k,m)
          wt(i,j,k,1,m)=wt(i,j,k,1,m)/n-wmean(i,j,k)*tmean(i,j,k,m)
          wt(i,j,k,2,m)=wt(i,j,k,2,m)/n-wmean(i,j,k)*tmean(i,j,k+1,m)
         enddo
        enddo
       enddo
      enddo
#ifdef diag_potential_vorticity
      do k=2,km-1
       do j=js,je 
        do i=2,imt-1
c       pv
         uq(i,j,k,1)=uq(i,j,k,1)/n-umean(i,j,k)*qmean(i ,j,k)
         uq(i,j,k,2)=uq(i,j,k,2)/n-umean(i,j,k)*qmean(i+1,j,k)
         vq(i,j,k,1)=vq(i,j,k,1)/n-vmean(i,j,k)*qmean(i,j,k)
         vq(i,j,k,2)=vq(i,j,k,2)/n-vmean(i,j,k)*qmean(i,j+1,k)
         wq(i,j,k,1)=wq(i,j,k,1)/n-wmean(i,j,k)*qmean(i,j,k)
         wq(i,j,k,2)=wq(i,j,k,2)/n-wmean(i,j,k)*qmean(i,j,k+1)
        enddo
       enddo
      enddo
#endif
c
c       Variance fluxes
c
      do k=2,km-1
       do j=js,je 
        do i=2,imt-1
c       buoyancy
         ubb(i,j,k,1)=ubb(i,j,k,1)/n-umean(i,j,k)*bmean(i,j,k)**2
     &                   -2*bmean(i  ,j,k)*ub(i,j,k,1)
         ubb(i,j,k,2)=ubb(i,j,k,2)/n-umean(i,j,k)*bmean(i+1,j,k)**2
     &                   -2*bmean(i+1,j,k)*ub(i,j,k,2)
         vbb(i,j,k,1)=vbb(i,j,k,1)/n-vmean(i,j,k)*bmean(i,j,k)**2
     &                   -2*bmean(i,j  ,k)*vb(i,j,k,1)
         vbb(i,j,k,2)=vbb(i,j,k,2)/n-vmean(i,j,k)*bmean(i,j+1,k)**2
     &                   -2*bmean(i,j+1,k)*vb(i,j,k,2)
         wbb(i,j,k,1)=wbb(i,j,k,1)/n-wmean(i,j,k)*bmean(i,j,k)**2
     &                   -2*bmean(i,j,k  )*wb(i,j,k,1)
         wbb(i,j,k,2)=wbb(i,j,k,2)/n-wmean(i,j,k)*bmean(i,j,k+1)**2
     &                   -2*bmean(i,j,k+1)*wb(i,j,k,2)
#ifdef diag_potential_vorticity
c       pv
         uqq(i,j,k,1)=uqq(i,j,k,1)/n-umean(i,j,k)*qmean(i,j,k)**2
     &                   -2*qmean(i  ,j,k)*uq(i,j,k,1)
         uqq(i,j,k,2)=uqq(i,j,k,2)/n-umean(i,j,k)*qmean(i+1,j,k)**2
     &                   -2*qmean(i+1,j,k)*uq(i,j,k,2)
         vqq(i,j,k,1)=vqq(i,j,k,1)/n-vmean(i,j,k)*qmean(i,j,k)**2
     &                   -2*qmean(i,j  ,k)*vq(i,j,k,1)
         vqq(i,j,k,2)=vqq(i,j,k,2)/n-vmean(i,j,k)*qmean(i,j+1,k)**2
     &                   -2*qmean(i,j+1,k)*vq(i,j,k,2)
         wqq(i,j,k,1)=wqq(i,j,k,1)/n-wmean(i,j,k)*qmean(i,j,k)**2
     &                   -2*qmean(i,j,k  )*wq(i,j,k,1)
         wqq(i,j,k,2)=wqq(i,j,k,2)/n-wmean(i,j,k)*qmean(i,j,k+1)**2
     &                   -2*qmean(i,j,k+1)*wq(i,j,k,2)
#endif

c       third order fluxes of buoyancy
        ub3(i,j,k,1)=ub3(i,j,k,1)/n-umean(i,j,k)*bmean(i,j,k)**3
     &           -3*bmean(i,j,k)*ubb(i,j,k,1)
     &           -3*bmean(i,j,k)**2*ub(i,j,k,1)
        ub3(i,j,k,2)=ub3(i,j,k,2)/n-umean(i,j,k)*bmean(i+1,j,k)**3
     &           -3*bmean(i+1,j,k)*ubb(i,j,k,2)
     &           -3*bmean(i+1,j,k)**2*ub(i,j,k,2)

        vb3(i,j,k,1)=vb3(i,j,k,1)/n-vmean(i,j,k)*bmean(i,j,k)**3
     &           -3*bmean(i,j,k)*vbb(i,j,k,1)
     &           -3*bmean(i,j,k)**2*vb(i,j,k,1)
        vb3(i,j,k,2)=vb3(i,j,k,2)/n-vmean(i,j,k)*bmean(i,j+1,k)**3
     &           -3*bmean(i,j+1,k)*vbb(i,j,k,2)
     &           -3*bmean(i,j+1,k)**2*vb(i,j,k,2)

        wb3(i,j,k,1)=wb3(i,j,k,1)/n-wmean(i,j,k)*bmean(i,j,k)**3
     &           -3*bmean(i,j,k)*wbb(i,j,k,1)
     &           -3*bmean(i,j,k)**2*wb(i,j,k,1)
        wb3(i,j,k,2)=wb3(i,j,k,2)/n-wmean(i,j,k)*bmean(i,j,k+1)**3
     &           -3*bmean(i,j,k+1)*wbb(i,j,k,2)
     &           -3*bmean(i,j,k+1)**2*wb(i,j,k,2)
#ifdef diag_potential_vorticity
c       third order fluxes of pv
        uq3(i,j,k,1)=uq3(i,j,k,1)/n-umean(i,j,k)*qmean(i,j,k)**3
     &           -3*qmean(i,j,k)*uqq(i,j,k,1)
     &           -3*qmean(i,j,k)**2*uq(i,j,k,1)
        uq3(i,j,k,2)=uq3(i,j,k,2)/n-umean(i,j,k)*qmean(i+1,j,k)**3
     &           -3*qmean(i+1,j,k)*uqq(i,j,k,2)
     &           -3*qmean(i+1,j,k)**2*uq(i,j,k,2)

        vq3(i,j,k,1)=vq3(i,j,k,1)/n-vmean(i,j,k)*qmean(i,j,k)**3
     &           -3*qmean(i,j,k)*vqq(i,j,k,1)
     &           -3*qmean(i,j,k)**2*vq(i,j,k,1)
        vq3(i,j,k,2)=vq3(i,j,k,2)/n-vmean(i,j,k)*qmean(i,j+1,k)**3
     &           -3*qmean(i,j+1,k)*vqq(i,j,k,2)
     &           -3*qmean(i,j+1,k)**2*vq(i,j,k,2)

        wq3(i,j,k,1)=wq3(i,j,k,1)/n-wmean(i,j,k)*qmean(i,j,k)**3
     &           -3*qmean(i,j,k)*wqq(i,j,k,1)
     &           -3*qmean(i,j,k)**2*wq(i,j,k,1)
        wq3(i,j,k,2)=wq3(i,j,k,2)/n-wmean(i,j,k)*qmean(i,j,k+1)**3
     &           -3*qmean(i,j,k+1)*wqq(i,j,k,2)
     &           -3*qmean(i,j,k+1)**2*wq(i,j,k,2)
#endif
        enddo
       enddo
      enddo

c       tracer
      do m=1,nr_tr
      do k=2,km-1
       do j=js,je 
        do i=2,imt-1
         utt(i,j,k,1,m)=utt(i,j,k,1,m)/n
     &                    -umean(i,j,k)*tmean(i,j,k,m)**2
     &                   -2*tmean(i  ,j,k,m)*ut(i,j,k,1,m)
         utt(i,j,k,2,m)=utt(i,j,k,2,m)/n
     &                    -umean(i,j,k)*tmean(i+1,j,k,m)**2
     &                   -2*tmean(i+1,j,k,m)*ut(i,j,k,2,m)
         vtt(i,j,k,1,m)=vtt(i,j,k,1,m)/n
     &                   -vmean(i,j,k)*tmean(i,j,k,m)**2
     &                   -2*tmean(i,j  ,k,m)*vt(i,j,k,1,m)
         vtt(i,j,k,2,m)=vtt(i,j,k,2,m)/n
     &                    -vmean(i,j,k)*tmean(i,j+1,k,m)**2
     &                   -2*tmean(i,j+1,k,m)*vt(i,j,k,2,m)
         wtt(i,j,k,1,m)=wtt(i,j,k,1,m)/n
     &                   -wmean(i,j,k)*tmean(i,j,k,m)**2
     &                   -2*tmean(i,j,k  ,m)*wt(i,j,k,1,m)
         wtt(i,j,k,2,m)=wtt(i,j,k,2,m)/n
     &                    -wmean(i,j,k)*tmean(i,j,k+1,m)**2
     &                   -2*tmean(i,j,k+1,m)*wt(i,j,k,2,m)
c       third order fluxes of tracer
        ut3(i,j,k,1,m)=ut3(i,j,k,1,m)/n-umean(i,j,k)*tmean(i,j,k,m)**3
     &           -3*tmean(i,j,k,m)*utt(i,j,k,1,m)
     &           -3*tmean(i,j,k,m)**2*ut(i,j,k,1,m)
        ut3(i,j,k,2,m)=ut3(i,j,k,2,m)/n-umean(i,j,k)*tmean(i+1,j,k,m)**3
     &           -3*tmean(i+1,j,k,m)*utt(i,j,k,2,m)
     &           -3*tmean(i+1,j,k,m)**2*ut(i,j,k,2,m)

        vt3(i,j,k,1,m)=vt3(i,j,k,1,m)/n-vmean(i,j,k)*tmean(i,j,k,m)**3
     &           -3*tmean(i,j,k,m)*vtt(i,j,k,1,m)
     &           -3*tmean(i,j,k,m)**2*vt(i,j,k,1,m)
        vt3(i,j,k,2,m)=vt3(i,j,k,2,m)/n-vmean(i,j,k)*tmean(i,j+1,k,m)**3
     &           -3*tmean(i,j+1,k,m)*vtt(i,j,k,2,m)
     &           -3*tmean(i,j+1,k,m)**2*vt(i,j,k,2,m)

        wt3(i,j,k,1,m)=wt3(i,j,k,1,m)/n-wmean(i,j,k)*tmean(i,j,k,m)**3
     &           -3*tmean(i,j,k,m)*wtt(i,j,k,1,m)
     &           -3*tmean(i,j,k,m)**2*wt(i,j,k,1,m)
        wt3(i,j,k,2,m)=wt3(i,j,k,2,m)/n-wmean(i,j,k)*tmean(i,j,k+1,m)**3
     &           -3*tmean(i,j,k+1,m)*wtt(i,j,k,2,m)
     &           -3*tmean(i,j,k+1,m)**2*wt(i,j,k,2,m)
        enddo
       enddo
      enddo
      enddo

      do k=2,km-1
       do j=js,je 
        do i=2,imt-1
         ub(i,j,k,1)=0.5*(ub(i,j,k,1)+ub(i,j,k,2))
         vb(i,j,k,1)=0.5*(vb(i,j,k,1)+vb(i,j,k,2))
         wb(i,j,k,1)=0.5*(wb(i,j,k,1)+wb(i,j,k,2))
         ubb(i,j,k,1)=0.5*(ubb(i,j,k,1)+ubb(i,j,k,2))/2.0
         vbb(i,j,k,1)=0.5*(vbb(i,j,k,1)+vbb(i,j,k,2))/2.0
         wbb(i,j,k,1)=0.5*(wbb(i,j,k,1)+wbb(i,j,k,2))/2.0
         ub3(i,j,k,1)=0.5*(ub3(i,j,k,1)+ub3(i,j,k,2))/3.0
         vb3(i,j,k,1)=0.5*(vb3(i,j,k,1)+vb3(i,j,k,2))/3.0
         wb3(i,j,k,1)=0.5*(wb3(i,j,k,1)+wb3(i,j,k,2))/3.0
#ifdef diag_potential_vorticity
         uq(i,j,k,1)=0.5*(uq(i,j,k,1)+uq(i,j,k,2))
         vq(i,j,k,1)=0.5*(vq(i,j,k,1)+vq(i,j,k,2))
         wq(i,j,k,1)=0.5*(wq(i,j,k,1)+wq(i,j,k,2))
         uqq(i,j,k,1)=0.5*(uqq(i,j,k,1)+uqq(i,j,k,2))/2.0
         vqq(i,j,k,1)=0.5*(vqq(i,j,k,1)+vqq(i,j,k,2))/2.0
         wqq(i,j,k,1)=0.5*(wqq(i,j,k,1)+wqq(i,j,k,2))/2.0
         uq3(i,j,k,1)=0.5*(uq3(i,j,k,1)+uq3(i,j,k,2))/3.0
         vq3(i,j,k,1)=0.5*(vq3(i,j,k,1)+vq3(i,j,k,2))/3.0
         wq3(i,j,k,1)=0.5*(wq3(i,j,k,1)+wq3(i,j,k,2))/3.0
#endif
        enddo
       enddo
      enddo

      do m=1,nr_tr
      do k=2,km-1
       do j=js,je 
        do i=2,imt-1
         ut(i,j,k,1,m)=0.5*(ut(i,j,k,1,m) +ut(i,j,k,2,m))
         vt(i,j,k,1,m)=0.5*(vt(i,j,k,1,m) +vt(i,j,k,2,m))
         wt(i,j,k,1,m)=0.5*(wt(i,j,k,1,m)+wt(i,j,k,2,m))

         utt(i,j,k,1,m)=0.5*(utt(i,j,k,1,m) +utt(i,j,k,2,m))/2.0
         vtt(i,j,k,1,m)=0.5*(vtt(i,j,k,1,m) +vtt(i,j,k,2,m))/2.0
         wtt(i,j,k,1,m)=0.5*(wtt(i,j,k,1,m)+wtt(i,j,k,2,m))/2.0

         ut3(i,j,k,1,m)=0.5*(ut3(i,j,k,1,m) +ut3(i,j,k,2,m))/2.0
         vt3(i,j,k,1,m)=0.5*(vt3(i,j,k,1,m) +vt3(i,j,k,2,m))/2.0
         wt3(i,j,k,1,m)=0.5*(wt3(i,j,k,1,m)+wt3(i,j,k,2,m))/2.0
        enddo
       enddo
      enddo
      enddo
#endif

#ifdef diag_eddy_flux_tendency
c
c       flux tendencies :  geostrophic terms
c
      do k=2,km-1
       do j=js,je 
        do i=2,imt-1
          bpx(i,j,k,1)=bpx(i,j,k,1)/n-bmean(i  ,j,k)*pmean(i  ,j,k)
          bpx(i,j,k,2)=bpx(i,j,k,2)/n-bmean(i+1,j,k)*pmean(i  ,j,k)
          bpx(i,j,k,3)=bpx(i,j,k,3)/n-bmean(i  ,j,k)*pmean(i+1,j,k)

          bpy(i,j,k,1)=bpy(i,j,k,1)/n-bmean(i,j  ,k)*pmean(i,j  ,k)
          bpy(i,j,k,2)=bpy(i,j,k,2)/n-bmean(i,j+1,k)*pmean(i,j  ,k)
          bpy(i,j,k,3)=bpy(i,j,k,3)/n-bmean(i,j  ,k)*pmean(i,j+1,k)

          bfv(i,j,k,1)=bfv(i,j,k,1)/n-bmean(i  ,j,k)*vmean(i  ,j  ,k)
          bfv(i,j,k,2)=bfv(i,j,k,2)/n-bmean(i+1,j,k)*vmean(i  ,j  ,k)
          bfv(i,j,k,3)=bfv(i,j,k,3)/n-bmean(i  ,j,k)*vmean(i+1,j  ,k)
          bfv(i,j,k,4)=bfv(i,j,k,4)/n-bmean(i  ,j,k)*vmean(i  ,j-1,k)
          bfv(i,j,k,5)=bfv(i,j,k,5)/n-bmean(i+1,j,k)*vmean(i  ,j-1,k)
          bfv(i,j,k,6)=bfv(i,j,k,6)/n-bmean(i  ,j,k)*vmean(i+1,j-1,k)

          bfu(i,j,k,1)=bfu(i,j,k,1)/n-bmean(i,j  ,k)*umean(i  ,j  ,k)
          bfu(i,j,k,2)=bfu(i,j,k,2)/n-bmean(i,j+1,k)*umean(i  ,j  ,k)
          bfu(i,j,k,3)=bfu(i,j,k,3)/n-bmean(i,j  ,k)*umean(i  ,j+1,k)
          bfu(i,j,k,4)=bfu(i,j,k,4)/n-bmean(i,j  ,k)*umean(i-1,j  ,k)
          bfu(i,j,k,5)=bfu(i,j,k,5)/n-bmean(i,j+1,k)*umean(i-1,j  ,k)
          bfu(i,j,k,6)=bfu(i,j,k,6)/n-bmean(i,j  ,k)*umean(i-1,j+1,k)
        enddo
       enddo
      enddo
      call border_exchg3D(bpx(:,:,:,1),1)
      call setcyclic3D(bpx(:,:,:,1))
      call border_exchg3D(bpy(:,:,:,1),1)
      call setcyclic3D(bpy(:,:,:,1))
      call border_exchg3D(bfv(:,:,:,1),1)
      call setcyclic3D(bfv(:,:,:,1))
      call border_exchg3D(bfv(:,:,:,4),1)
      call setcyclic3D(bfv(:,:,:,4))
      call border_exchg3D(bfu(:,:,:,1),1)
      call setcyclic3D(bfu(:,:,:,1))
      call border_exchg3D(bfu(:,:,:,4),1)
      call setcyclic3D(bfu(:,:,:,4))

      ugeo1=0;ugeo2=0; vgeo1=0;vgeo2=0
      do k=2,km-1
       do j=js,je 
        do i=2,imt-1
         ugeo1(i,j,k)=-(bpx(i+1,j,k,1)-bpx(i,j,k,2)
     &                 +bpx(i,j,k,3)-bpx(i,j,k,1))/(2*dx)
         vgeo1(i,j,k)=-(bpy(i,j+1,k,1)-bpy(i,j,k,2)
     &                 +bpy(i,j,k,3)-bpy(i,j,k,1))/(2*dx)
         ugeo2(i,j,k)=(bfv(i,j,k,1)+bfv(i,j,k,2)
     &                +bfv(i,j,k,3)+bfv(i+1,j,k,1) 
     &                +bfv(i,j,k,4)+bfv(i,j,k,5)
     &                +bfv(i,j,k,6)+bfv(i+1,j,k,4) )/8.*coriolis_t(j)
         vgeo2(i,j,k)=-(bfu(i,j,k,1)+bfu(i,j,k,2)
     &                 +bfu(i,j,k,4)+bfu(i,j,k,5) )/8.*coriolis_t(j)
     &                -(bfu(i,j+1,k,1)+bfu(i,j+1,k,4)
     &                 +bfu(i,j,k,6)+bfu(i,j,k,3))/8.*coriolis_t(j+1)
        enddo
       enddo
      enddo
c
c       flux tendencies:   vertical advection
c
      do k=2,km-1
       do j=js,je 
        do i=2,imt-1
          uwN(i,j,k,1)=uwN(i,j,k,1)/n-umean(i,j,k)*wmean(i  ,j,k)
          uwN(i,j,k,2)=uwN(i,j,k,2)/n-umean(i,j,k)*wmean(i+1,j,k)
          uwN(i,j,k,3)=uwN(i,j,k,3)/n-umean(i,j,k)*wmean(i  ,j,k-1)
          uwN(i,j,k,4)=uwN(i,j,k,4)/n-umean(i,j,k)*wmean(i+1,j,k-1)
          vwN(i,j,k,1)=vwN(i,j,k,1)/n-vmean(i,j,k)*wmean(i,j  ,k)
          vwN(i,j,k,2)=vwN(i,j,k,2)/n-vmean(i,j,k)*wmean(i,j+1,k)
          vwN(i,j,k,3)=vwN(i,j,k,3)/n-vmean(i,j,k)*wmean(i,j  ,k-1)
          vwN(i,j,k,4)=vwN(i,j,k,4)/n-vmean(i,j,k)*wmean(i,j+1,k-1)
        enddo
       enddo
      enddo

      do k=2,km-1
       do j=js,je 
        do i=2,imt-1
          uwN(i,j,k,1)=
     &        uwN(i,j,k,1)*(bmean(i,j,k+1)+bmean(i,j,k))/(4*dz)
     &       +uwN(i,j,k,2)*(bmean(i+1,j,k+1)+bmean(i+1,j,k))/(4*dz)
     &       -uwN(i,j,k,3)*(bmean(i,j,k)+bmean(i,j,k-1))/(4*dz)
     &       -uwN(i,j,k,4)*(bmean(i+1,j,k)+bmean(i+1,j,k-1))/(4*dz)
#ifdef  remove_back_stratification
     &       +uwN(i,j,k,1)*(b_r(k+1)+b_r(k))/(4*dz)
     &       +uwN(i,j,k,2)*(b_r(k+1)+b_r(k))/(4*dz)
     &       -uwN(i,j,k,3)*(b_r(k)+b_r(k-1))/(4*dz)
     &       -uwN(i,j,k,4)*(b_r(k)+b_r(k-1))/(4*dz)
#endif
          vwN(i,j,k,1)=
     &        vwN(i,j,k,1)*(bmean(i,j,k+1)+bmean(i,j,k))/(4*dz)
     &       +vwN(i,j,k,2)*(bmean(i,j+1,k+1)+bmean(i,j+1,k))/(4*dz)
     &       -vwN(i,j,k,3)*(bmean(i,j,k)+bmean(i,j,k-1))/(4*dz)
     &       -vwN(i,j,k,4)*(bmean(i,j+1,k)+bmean(i,j+1,k-1))/(4*dz)
#ifdef  remove_back_stratification
     &       +vwN(i,j,k,1)*(b_r(k+1)+b_r(k))/(4*dz)
     &       +vwN(i,j,k,2)*(b_r(k+1)+b_r(k))/(4*dz)
     &       -vwN(i,j,k,3)*(b_r(k)+b_r(k-1))/(4*dz)
     &       -vwN(i,j,k,4)*(b_r(k)+b_r(k-1))/(4*dz)
#endif
        enddo
       enddo
      enddo
#endif
c
c       apply  land mask
c
      where (maskU==0) 
        umean=spval; uvar=spval
      end where
      where (maskV==0) 
        vmean=spval; wvar=spval
      end where
      where (maskW==0) 
        wmean=spval; wvar=spval
      end where
      where (maskT==0) 
       bmean=spval;pmean=spval;b3=spval
       bvar=spval;pvar=spval; uv(:,:,:,1)=spval
      end where
#ifdef diag_eddy_fluxes
      where (maskU==0) 
        ub(:,:,:,1)=spval; ubb(:,:,:,1)=spval;ub3(:,:,:,1)=spval
      end where
      where (maskV==0) 
        vb(:,:,:,1)=spval; vbb(:,:,:,1)=spval;vb3(:,:,:,1)=spval
      end where
      where (maskW==0) 
        wb(:,:,:,1)=spval; wbb(:,:,:,1)=spval;wb3(:,:,:,1)=spval
      end where
#endif
#ifdef diag_potential_vorticity
      where (maskT==0) 
       qmean=spval;q3=spval; qvar=spval;
      end where
      where (maskU==0) 
        uq(:,:,:,1)=spval; uqq(:,:,:,1)=spval;uq3(:,:,:,1)=spval
      end where
      where (maskV==0) 
        vq(:,:,:,1)=spval; vqq(:,:,:,1)=spval;vq3(:,:,:,1)=spval
      end where
      where (maskW==0) 
        wq(:,:,:,1)=spval; wqq(:,:,:,1)=spval;wq3(:,:,:,1)=spval
      end where
#endif

#ifdef diag_eddy_flux_tendency
      where (maskU==0) 
        ugeo1=spval; ugeo2=spval; uwN(:,:,:,1)=spval
      end where
      where (maskV==0) 
        vgeo1=spval; vgeo2=spval; vwN(:,:,:,1)=spval
      end where
#endif
      do k=1,nr_tr
       where (maskT==0) 
        tmean(:,:,:,k)=spval;tvar(:,:,:,k)=spval;t3(:,:,:,k)=spval
       end where
#ifdef diag_eddy_fluxes
       where (maskU==0) 
         ut(:,:,:,1,k)=spval; utt(:,:,:,1,k)=spval;ut3(:,:,:,1,k)=spval
       end where
       where (maskV==0) 
         vt(:,:,:,1,k)=spval; vtt(:,:,:,1,k)=spval;vt3(:,:,:,1,k)=spval
       end where
       where (maskW==0) 
         wt(:,:,:,1,k)=spval; wtt(:,:,:,1,k)=spval;wt3(:,:,:,1,k)=spval
       end where
#endif
      enddo
      end subroutine time_averages_calc



      subroutine time_averages_write
c=======================================================================
c        write time mean diagnostics
c=======================================================================
      use cpflame_module
      use time_averages_module
      implicit none
#include "netcdf.inc"
      integer :: ncid,iret,n,npe, corner(4), edges(4)
      integer :: tdim,ilen,timeid,varid,dims(4)
      integer :: xtdim,xudim,ytdim,yudim,ztdim,zudim
      real :: fxa
      type(time_type) :: time
      character :: name*64, unit*32
c
c        define netcdf file
c
      if (my_pe==0) then 
       print*,' writing time averages to file time_ave.cdf'
       call def_grid_cdf('time_ave.cdf')
       iret=nf_open('time_ave.cdf',NF_WRITE,ncid)
       iret=nf_set_fill(ncid, NF_NOFILL, iret)
       call ncredf(ncid, iret)
       iret=nf_inq_dimid(ncid,'xt',xtdim)
       iret=nf_inq_dimid(ncid,'xu',xudim)
       iret=nf_inq_dimid(ncid,'yt',ytdim)
       iret=nf_inq_dimid(ncid,'yu',yudim)
       iret=nf_inq_dimid(ncid,'zt',ztdim)
       iret=nf_inq_dimid(ncid,'zu',zudim)
       iret=nf_inq_dimid(ncid,'Time',tdim)

       dims = (/xudim, ytdim,ztdim, tdim/)
       varid = ncvdef (ncid,'um', NCFLOAT,4,dims,iret)
       name='Mean zonal velocity'; unit = 'm/s'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xtdim, yudim,ztdim, tdim/)
       varid = ncvdef (ncid,'vm', NCFLOAT,4,dims,iret)
       name='Mean meridional velocity'; unit = 'm/s'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xtdim, ytdim,zudim, tdim/)
       varid = ncvdef (ncid,'wm', NCFLOAT,4,dims,iret)
       name='Mean vertical velocity'; unit = 'm/s'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xtdim, ytdim,ztdim, tdim/)
       varid = ncvdef (ncid,'bm', NCFLOAT,4,dims,iret)
       name='Mean buoyancy'; unit = 'm/s^2'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xtdim, ytdim,ztdim, tdim/)
       varid = ncvdef (ncid,'pm', NCFLOAT,4,dims,iret)
       name='Mean pressure'; unit = 'm^2/s^2'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xudim, ytdim,ztdim, tdim/)
       varid = ncvdef (ncid,'uu', NCFLOAT,4,dims,iret)
       name='Zonal velocity variance'; unit = 'm^2/s^2'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xtdim, yudim,ztdim, tdim/)
       varid = ncvdef (ncid,'vv', NCFLOAT,4,dims,iret)
       name='Meridional velocity variance'; unit = 'm^2/s^2'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xtdim, ytdim,ztdim, tdim/)
       varid = ncvdef (ncid,'uv', NCFLOAT,4,dims,iret)
       name='Cross velocity variance'; unit = 'm^2/s^2'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xtdim, ytdim,zudim, tdim/)
       varid = ncvdef (ncid,'ww', NCFLOAT,4,dims,iret)
       name='Vertical velocity variance'; unit = 'm^2/s^2'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xtdim, ytdim,ztdim, tdim/)
       varid = ncvdef (ncid,'bb', NCFLOAT,4,dims,iret)
       name='Buoyancy variance'; unit = 'm^2/s^4'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xtdim, ytdim,ztdim, tdim/)
       varid = ncvdef (ncid,'b3', NCFLOAT,4,dims,iret)
       name='3. moment buoyancy'; unit = 'm^3/s^6'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xtdim, ytdim,ztdim, tdim/)
       varid = ncvdef (ncid,'pp', NCFLOAT,4,dims,iret)
       name='Pressure variance'; unit = 'm^4/s^4'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

#ifdef diag_eddy_fluxes
       dims = (/xudim, ytdim,ztdim, tdim/)
       varid = ncvdef (ncid,'ub', NCFLOAT,4,dims,iret)
       name='Zonal eddy buoyancy flux'; unit = 'm^2/s^3'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xtdim, yudim,ztdim, tdim/)
       varid = ncvdef (ncid,'vb', NCFLOAT,4,dims,iret)
       name='Meridional eddy buoyancy flux'; unit = 'm^2/s^3'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xtdim, ytdim,zudim, tdim/)
       varid = ncvdef (ncid,'wb', NCFLOAT,4,dims,iret)
       name='Vertical eddy buoyancy flux'; unit = 'm^2/s^3'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xudim, ytdim,ztdim, tdim/)
       varid = ncvdef (ncid,'ubb', NCFLOAT,4,dims,iret)
       name='Zonal buoyancy variance flux'; unit = 'm^3/s^5'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xtdim, yudim,ztdim, tdim/)
       varid = ncvdef (ncid,'vbb', NCFLOAT,4,dims,iret)
       name='Meridional buoyancy variance flux'; unit = 'm^3/s^5'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xtdim, ytdim,zudim, tdim/)
       varid = ncvdef (ncid,'wbb', NCFLOAT,4,dims,iret)
       name='Vertical buoyancy variance flux'; unit = 'm^3/s^5'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xudim, ytdim,ztdim, tdim/)
       varid = ncvdef (ncid,'ub3', NCFLOAT,4,dims,iret)
       name='Zonal buoyancy 3.order flux'; unit = 'm^4/s^7'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xtdim, yudim,ztdim, tdim/)
       varid = ncvdef (ncid,'vb3', NCFLOAT,4,dims,iret)
       name='Meridional buoyancy 3. order flux'; unit = 'm^4/s^7'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xtdim, ytdim,zudim, tdim/)
       varid = ncvdef (ncid,'wb3', NCFLOAT,4,dims,iret)
       name='Vertical buoyancy 3. order flux'; unit = 'm^4/s^7'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

#endif

#ifdef diag_eddy_flux_tendency
       dims = (/xudim, ytdim,ztdim, tdim/)
       varid = ncvdef (ncid,'ugeo1', NCFLOAT,4,dims,iret)
       name='Geostr. flux tendency 1'; unit = 'm^2/s^4'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xudim, ytdim,ztdim, tdim/)
       varid = ncvdef (ncid,'ugeo2', NCFLOAT,4,dims,iret)
       name='Geostr. flux tendency 2'; unit = 'm^2/s^4'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xtdim, yudim,ztdim, tdim/)
       varid = ncvdef (ncid,'vgeo1', NCFLOAT,4,dims,iret)
       name='Geostr. flux tendency 1'; unit = 'm^2/s^4'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xtdim, yudim,ztdim, tdim/)
       varid = ncvdef (ncid,'vgeo2', NCFLOAT,4,dims,iret)
       name='Geostr. flux tendency 2'; unit = 'm^2/s^4'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xudim, ytdim,ztdim, tdim/)
       varid = ncvdef (ncid,'uwN', NCFLOAT,4,dims,iret)
       name='Vertical flux tendency'; unit = 'm^2/s^4'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xtdim, yudim,ztdim, tdim/)
       varid = ncvdef (ncid,'vwN', NCFLOAT,4,dims,iret)
       name='Vertical flux tendency'; unit = 'm^2/s^4'
       call dvcdf(ncid,varid,name,64,unit,32,spval)
#endif

#ifdef  remove_back_stratification
       varid = ncvdef (ncid,'b_r', NCFLOAT,1,ztdim,iret)
       name = 'Reference buoyancy profile'; unit = 'm/s^2'
       call dvcdf(ncid,varid,name,64,unit,32,spval)
#endif
       do n=1,nr_tr
        dims = (/xtdim, ytdim,ztdim, tdim/)
        write(name, '("trm",i2)') n; call replace_space_zero(name)
        varid = ncvdef (ncid,name, NCFLOAT,4,dims,iret)
        name = 'Mean tracer'; unit = ' '
        call dvcdf(ncid,varid,name,64,unit,32,spval)

        write(name, '("trvar",i2)') n; call replace_space_zero(name)
        varid = ncvdef (ncid,name, NCFLOAT,4,dims,iret)
        name = 'Tracer variance'; unit = ' '
        call dvcdf(ncid,varid,name,64,unit,32,spval)

        write(name, '("tr3",i2)') n; call replace_space_zero(name)
        varid = ncvdef (ncid,name, NCFLOAT,4,dims,iret)
        name = '3. moment tracer'; unit = ' '
        call dvcdf(ncid,varid,name,64,unit,32,spval)
#ifdef diag_eddy_fluxes
        dims = (/xudim, ytdim,ztdim, tdim/)
        write(name, '("utr",i2)') n; call replace_space_zero(name)
        varid = ncvdef (ncid,name, NCFLOAT,4,dims,iret)
        name='Zonal eddy tracer flux'; unit = 'm/s '
        call dvcdf(ncid,varid,name,64,unit,32,spval)

        dims = (/xtdim, yudim,ztdim, tdim/)
        write(name, '("vtr",i2)') n; call replace_space_zero(name)
        varid = ncvdef (ncid,name, NCFLOAT,4,dims,iret)
        name='Meridional eddy tracer flux'; unit = 'm/s '
        call dvcdf(ncid,varid,name,64,unit,32,spval)

        dims = (/xtdim, ytdim,zudim, tdim/)
        write(name, '("wtr",i2)') n; call replace_space_zero(name)
        varid = ncvdef (ncid,name, NCFLOAT,4,dims,iret)
        name='Vertical eddy tracer flux'; unit = 'm/s '
        call dvcdf(ncid,varid,name,64,unit,32,spval)

        dims = (/xudim, ytdim,ztdim, tdim/)
        write(name, '("utrvar",i2)') n; call replace_space_zero(name)
        varid = ncvdef (ncid,name, NCFLOAT,4,dims,iret)
        name='Zonal tracer variance flux'; unit = ' m/s'
        call dvcdf(ncid,varid,name,64,unit,32,spval)

        dims = (/xtdim, yudim,ztdim, tdim/)
        write(name, '("vtrvar",i2)') n; call replace_space_zero(name)
        varid = ncvdef (ncid,name, NCFLOAT,4,dims,iret)
        name='Meridional tracer variance flux'; unit = 'm/s'
        call dvcdf(ncid,varid,name,64,unit,32,spval)

        dims = (/xtdim, ytdim,zudim, tdim/)
        write(name, '("wtrvar",i2)') n; call replace_space_zero(name)
        varid = ncvdef (ncid,name, NCFLOAT,4,dims,iret)
        name='Vertical tracer variance flux'; unit = 'm/s'
        call dvcdf(ncid,varid,name,64,unit,32,spval)

        dims = (/xudim, ytdim,ztdim, tdim/)
        write(name, '("ut3_",i2)') n; call replace_space_zero(name)
        varid = ncvdef (ncid,name, NCFLOAT,4,dims,iret)
        name='Zonal tracer  3.order flux'; unit = ' m/s'
        call dvcdf(ncid,varid,name,64,unit,32,spval)

        dims = (/xtdim, yudim,ztdim, tdim/)
        write(name, '("vt3_",i2)') n; call replace_space_zero(name)
        varid = ncvdef (ncid,name, NCFLOAT,4,dims,iret)
        name='Meridional tracer 3.order flux'; unit = 'm/s'
        call dvcdf(ncid,varid,name,64,unit,32,spval)

        dims = (/xtdim, ytdim,zudim, tdim/)
        write(name, '("wt3_",i2)') n; call replace_space_zero(name)
        varid = ncvdef (ncid,name, NCFLOAT,4,dims,iret)
        name='Vertical tracer 3.order flux'; unit = 'm/s'
        call dvcdf(ncid,varid,name,64,unit,32,spval)

#endif
       enddo
#ifdef diag_potential_vorticity
       dims = (/xtdim, ytdim,ztdim, tdim/)
       varid = ncvdef (ncid,'qm', NCFLOAT,4,dims,iret)
       name='Ertel PV'; unit = '1/s'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xtdim, ytdim,ztdim, tdim/)
       varid = ncvdef (ncid,'qq', NCFLOAT,4,dims,iret)
       name='PV variance'; unit = '1/s^2'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xtdim, ytdim,ztdim, tdim/)
       varid = ncvdef (ncid,'q3', NCFLOAT,4,dims,iret)
       name='3. moment pv'; unit = '1/s^3'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xudim, ytdim,ztdim, tdim/)
       varid = ncvdef (ncid,'uq', NCFLOAT,4,dims,iret)
       name='Zonal eddy pv flux'; unit = 'm/s^2'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xtdim, yudim,ztdim, tdim/)
       varid = ncvdef (ncid,'vq', NCFLOAT,4,dims,iret)
       name='Meridional eddy pv flux'; unit = 'm/s^2'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xtdim, ytdim,zudim, tdim/)
       varid = ncvdef (ncid,'wq', NCFLOAT,4,dims,iret)
       name='Vertical eddy pv flux'; unit = 'm/s^2'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xudim, ytdim,ztdim, tdim/)
       varid = ncvdef (ncid,'uqq', NCFLOAT,4,dims,iret)
       name='Zonal pv variance flux'; unit = 'm^2/s^3'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xtdim, yudim,ztdim, tdim/)
       varid = ncvdef (ncid,'vqq', NCFLOAT,4,dims,iret)
       name='Meridional pv variance flux'; unit = 'm^2/s^3'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xtdim, ytdim,zudim, tdim/)
       varid = ncvdef (ncid,'wqq', NCFLOAT,4,dims,iret)
       name='Vertical pv variance flux'; unit = 'm^2/s^3'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xudim, ytdim,ztdim, tdim/)
       varid = ncvdef (ncid,'uq3', NCFLOAT,4,dims,iret)
       name='Zonal pv 3.order flux'; unit = 'm^3/s^4'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xtdim, yudim,ztdim, tdim/)
       varid = ncvdef (ncid,'vq3', NCFLOAT,4,dims,iret)
       name='Meridional pv 3. order flux'; unit = 'm^3/s^4'
       call dvcdf(ncid,varid,name,64,unit,32,spval)

       dims = (/xtdim, ytdim,zudim, tdim/)
       varid = ncvdef (ncid,'wq3', NCFLOAT,4,dims,iret)
       name='Vertical pv 3. order flux'; unit = 'm^3/s^4'
       call dvcdf(ncid,varid,name,64,unit,32,spval)
#endif

       iret= nf_put_att_int(ncid,nf_global,'nitts',nf_int,1,nitts)
       call ncendf(ncid, iret)
#ifdef  remove_back_stratification
       iret=nf_inq_varid(ncid,'b_r',varid)
       iret= nf_put_vara_double(ncid,varid,1,km,b_r)
#endif
       call ncclos (ncid, iret)
      endif ! my_pe ==0
c-----------------------------------------------------------------------
c      write to netcdf file
c-----------------------------------------------------------------------
      do npe=0,n_pes
       if (my_pe==npe) then

         iret=nf_open('time_ave.cdf',NF_WRITE,ncid)
         iret=nf_set_fill(ncid, NF_NOFILL, iret)

         iret=nf_inq_dimid(ncid,'Time',tdim)
         iret=nf_inq_dimlen(ncid,tdim,ilen)
         iret=nf_inq_varid(ncid,'Time',timeid)
         if (my_pe==0) then
          ilen=ilen+1
          time = current_time-initial_time
          fxa = time%days + time%seconds/86400.
          iret= nf_put_vara_double(ncid,timeid,ilen,1,fxa)
         endif

         Corner = (/1,js_pe,1,ilen/); 
         edges  = (/imt,je_pe-js_pe+1,km,1/)
         iret=nf_inq_varid(ncid,'um',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          umean(:,js_pe:je_pe,:))
         iret=nf_inq_varid(ncid,'vm',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          vmean(:,js_pe:je_pe,:))
         iret=nf_inq_varid(ncid,'wm',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          wmean(:,js_pe:je_pe,:))
         iret=nf_inq_varid(ncid,'bm',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          bmean(:,js_pe:je_pe,:))
         iret=nf_inq_varid(ncid,'pm',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          pmean(:,js_pe:je_pe,:))

         iret=nf_inq_varid(ncid,'uu',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          uvar(:,js_pe:je_pe,:))
         iret=nf_inq_varid(ncid,'vv',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          vvar(:,js_pe:je_pe,:))
         iret=nf_inq_varid(ncid,'uv',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          uv(:,js_pe:je_pe,:,1))
         iret=nf_inq_varid(ncid,'ww',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          wvar(:,js_pe:je_pe,:))
         iret=nf_inq_varid(ncid,'bb',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          bvar(:,js_pe:je_pe,:))
         iret=nf_inq_varid(ncid,'b3',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          b3(:,js_pe:je_pe,:))
         iret=nf_inq_varid(ncid,'pp',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          pvar(:,js_pe:je_pe,:))

#ifdef diag_eddy_fluxes
         iret=nf_inq_varid(ncid,'ub',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          ub(:,js_pe:je_pe,:,1))
         iret=nf_inq_varid(ncid,'vb',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          vb(:,js_pe:je_pe,:,1))
         iret=nf_inq_varid(ncid,'wb',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          wb(:,js_pe:je_pe,:,1))

         iret=nf_inq_varid(ncid,'ubb',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          ubb(:,js_pe:je_pe,:,1))
         iret=nf_inq_varid(ncid,'vbb',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          vbb(:,js_pe:je_pe,:,1))
         iret=nf_inq_varid(ncid,'wbb',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          wbb(:,js_pe:je_pe,:,1))

         iret=nf_inq_varid(ncid,'ub3',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          ub3(:,js_pe:je_pe,:,1))
         iret=nf_inq_varid(ncid,'vb3',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          vb3(:,js_pe:je_pe,:,1))
         iret=nf_inq_varid(ncid,'wb3',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          wb3(:,js_pe:je_pe,:,1))
#endif

#ifdef diag_eddy_flux_tendency
         iret=nf_inq_varid(ncid,'ugeo1',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          ugeo1(:,js_pe:je_pe,:))
         iret=nf_inq_varid(ncid,'ugeo2',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          ugeo2(:,js_pe:je_pe,:))
         iret=nf_inq_varid(ncid,'vgeo1',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          vgeo1(:,js_pe:je_pe,:))
         iret=nf_inq_varid(ncid,'vgeo2',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          vgeo2(:,js_pe:je_pe,:))

         iret=nf_inq_varid(ncid,'uwN',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          uwN(:,js_pe:je_pe,:,1))
         iret=nf_inq_varid(ncid,'vwN',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          vwN(:,js_pe:je_pe,:,1))
#endif
        do n=1,nr_tr
         write(name, '("trm",i2)') n; call replace_space_zero(name)
          iret=nf_inq_varid(ncid,name,varid)
          iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          tmean(:,js_pe:je_pe,:,n))

         write(name, '("trvar",i2)') n; call replace_space_zero(name)
          iret=nf_inq_varid(ncid,name,varid)
          iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          tvar(:,js_pe:je_pe,:,n))

         write(name, '("tr3",i2)') n; call replace_space_zero(name)
          iret=nf_inq_varid(ncid,name,varid)
          iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          t3(:,js_pe:je_pe,:,n))
#ifdef diag_eddy_fluxes
          write(name, '("utr",i2)') n; call replace_space_zero(name)
          iret=nf_inq_varid(ncid,name,varid)
          iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          ut(:,js_pe:je_pe,:,1,n))
          write(name, '("vtr",i2)') n; call replace_space_zero(name)
          iret=nf_inq_varid(ncid,name,varid)
          iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          vt(:,js_pe:je_pe,:,1,n))
          write(name, '("wtr",i2)') n; call replace_space_zero(name)
          iret=nf_inq_varid(ncid,name,varid)
          iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          wt(:,js_pe:je_pe,:,1,n))

          write(name, '("utrvar",i2)') n; call replace_space_zero(name)
          iret=nf_inq_varid(ncid,name,varid)
          iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          utt(:,js_pe:je_pe,:,1,n))
          write(name, '("vtrvar",i2)') n; call replace_space_zero(name)
          iret=nf_inq_varid(ncid,name,varid)
          iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          vtt(:,js_pe:je_pe,:,1,n))
          write(name, '("wtrvar",i2)') n; call replace_space_zero(name)
          iret=nf_inq_varid(ncid,name,varid)
          iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          wtt(:,js_pe:je_pe,:,1,n))

          write(name, '("ut3_",i2)') n; call replace_space_zero(name)
          iret=nf_inq_varid(ncid,name,varid)
          iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          ut3(:,js_pe:je_pe,:,1,n))
          write(name, '("vt3_",i2)') n; call replace_space_zero(name)
          iret=nf_inq_varid(ncid,name,varid)
          iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          vt3(:,js_pe:je_pe,:,1,n))
          write(name, '("wt3_",i2)') n; call replace_space_zero(name)
          iret=nf_inq_varid(ncid,name,varid)
          iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          wt3(:,js_pe:je_pe,:,1,n))

#endif
         enddo
#ifdef diag_potential_vorticity
         iret=nf_inq_varid(ncid,'qm',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          qmean(:,js_pe:je_pe,:))
         iret=nf_inq_varid(ncid,'qq',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          qvar(:,js_pe:je_pe,:))
         iret=nf_inq_varid(ncid,'q3',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          q3(:,js_pe:je_pe,:))

         iret=nf_inq_varid(ncid,'uq',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          uq(:,js_pe:je_pe,:,1))
         iret=nf_inq_varid(ncid,'vq',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          vq(:,js_pe:je_pe,:,1))
         iret=nf_inq_varid(ncid,'wq',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          wq(:,js_pe:je_pe,:,1))

         iret=nf_inq_varid(ncid,'uqq',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          uqq(:,js_pe:je_pe,:,1))
         iret=nf_inq_varid(ncid,'vqq',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          vqq(:,js_pe:je_pe,:,1))
         iret=nf_inq_varid(ncid,'wqq',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          wqq(:,js_pe:je_pe,:,1))

         iret=nf_inq_varid(ncid,'uq3',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          uq3(:,js_pe:je_pe,:,1))
         iret=nf_inq_varid(ncid,'vq3',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          vq3(:,js_pe:je_pe,:,1))
         iret=nf_inq_varid(ncid,'wq3',varid)
         iret= nf_put_vara_double(ncid,varid,corner,edges, 
     &          wq3(:,js_pe:je_pe,:,1))
#endif
         call ncclos (ncid, iret)
       endif
       call barrier
      enddo
      end subroutine time_averages_write

#else
      subroutine time_averages_dummy
      end
#endif
