#include "options.inc"


c
c=======================================================================
c  convert binary output of averaged second moments to NetCDF format
c=======================================================================
c
c
c  real variables are 8 or 4 bytes long set by
c#define real_8
c
c#define large_fluxfile
#define define_in_fluxfile_eddyflux
#define define_in_fluxfile_varianceflux

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

      end module kind_mod



      module cv_var_module
      use kind_mod
      integer :: imt,jmt,km,nt
      real (kind=r_knd), allocatable, dimension(:,:,:,:) :: mean_u,var_u
      real (kind=r_knd), allocatable, dimension(:,:) :: mean_ssh,var_ssh
      integer :: counter
      real (kind=r_knd), allocatable, dimension(:)  :: xt,yt,xu,yu,zt,zw
      integer, allocatable, dimension(:,:)     :: kmt,kmu
      logical :: enable_monthly_variances = .false.
      logical :: enable_seasonal_variances = .false.
      logical :: enable_freesurf = .false.
      logical :: energy = .false., energy_ext = .false.
      logical :: var_fluxes = .false.
      real (kind=r_knd), allocatable, dimension(:,:,:) :: uv,rhom
      real (kind=r_knd), allocatable, dimension(:,:,:,:) :: urho,upress
      real (kind=r_knd), allocatable, dimension(:,:,:) :: wrho,pressm
      real (kind=r_knd),allocatable,dimension(:,:,:,:)  ::mean_adv
      real (kind=r_knd),allocatable,dimension(:,:,:,:,:)::mean_tracer
      real (kind=r_knd),allocatable,dimension(:,:,:,:,:)::var_tracer
      real (kind=r_knd),allocatable,dimension(:,:,:,:,:)::mean_eddy_flux
      real (kind=r_knd),allocatable,dimension(:,:,:,:,:)::mean_var_flux
#ifdef large_fluxfile
      real (kind=r_knd),allocatable,dimension(:,:,:,:,:)::theta
      real (kind=r_knd),allocatable,dimension(:,:,:,:,:)::eddy_flux
#endif
      end module cv_var_module


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

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

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

      read(infile,'(a10,i3,a2,i4,a1,i2,a1,i2,a4)') 
     &            c9,domain,c2,year,c1,month,c1,day,c4
      print*,' sub domain = ',domain
      print*,' year/month/day  = ',year,month,day
      call read_variances(infile)

#ifdef notdef
c old version
      write(outfile,'("variances_",i3,"_y",i4,"m",i2,"d",i2,".cdf")') 
     &            domain,year,month,day
      call replace_space_zero(outfile)
      call write_variances(outfile,year,month,day,domain)
#endif

c   new version, split into several files

      write(outfile,'("eke_",i3,"_y",i4,"m",i2,"d",i2,".cdf")') 
     &            domain,year,month,day
      call replace_space_zero(outfile)
      call write_variances_eke(outfile,year,month,day,domain)

      if (energy) then
       write(outfile,'("energy_",i3,"_y",i4,"m",i2,"d",i2,".cdf")') 
     &            domain,year,month,day
       call replace_space_zero(outfile)
       call write_variances_energy(outfile,year,month,day,domain)
      endif

      if (var_fluxes) then
       write(outfile,'("varflux_",i3,"_y",i4,"m",i2,"d",i2,".cdf")') 
     &            domain,year,month,day
       call replace_space_zero(outfile)
       call write_variances_fluxes(outfile,year,month,day,domain)
      endif

      stop

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


      subroutine read_variances(infile)
      use cv_var_module
      implicit none
      character (len=*) :: infile
      character (len=3) :: ssh
      real (kind=r_knd), allocatable :: buf(:,:)
      integer :: i,j,k,io,n,nr,m
      logical :: eof = .false.
      character (len=80) :: name

      call getunit(io,infile,'usr ieee')
      read(io) imt,jmt,km,nt
      print*,' imt=',imt,' jmt=',jmt,' km=',km,' nt=',nt
      allocate( var_u(imt,km,jmt,nt) ); var_u=0.0
      allocate( mean_u(imt,km,jmt,2) ); mean_u=0.0
      allocate( mean_ssh(imt,jmt) );mean_ssh=0.0
      allocate( var_ssh(imt,jmt) ); var_ssh=0.0
      allocate( xt(imt),yt(jmt),xu(imt),yu(jmt),zt(km),zw(0:km))
      zw=0.0
      allocate( kmt(imt,jmt),kmu(imt,jmt) ); kmu=0

      read(io) xt,yt,xu,yu,zt,zw(1:km)
      read(io) counter
      print*,' counter for averages=',counter
      read(io) enable_monthly_variances
      read(io) enable_seasonal_variances
      read(io) kmt
      do j=1,jmt-1
        do i=1,imt-1
           kmu(i,j) = min (kmt(i,j), kmt(i+1,j), 
     &                     kmt(i,j+1), kmt(i+1,j+1))
        enddo
      enddo

      print*,' zt(1:5)=',zt(1:5)
      allocate( buf(imt,jmt) )

      do while (.not. eof) 
       read (io,end=50) nr,name
       goto 60
 50    eof=.true.
 60    continue

       if (.not. eof) then

        if (name(1:5) == 'umean') then
  
         print*,' found mean velocity '
         do n=1,2
          do k=1,km
           print*,'k=',k,'/',km
           read(io) buf
           mean_u(:,k,:,n)=  buf
          enddo
         enddo

        elseif (name(1:4) == 'uvar' ) then

         print*,' found variance of velocity '

         do n=1,2
          do k=1,km
           print*,'k=',k,'/',km
           read(io) buf
           var_u(:,k,:,n)=  buf
          enddo
         enddo

        elseif (name(1:6) == 'energy') then

         print*,' found energy conversion terms '

         energy=.true.
         allocate( uv(imt,km,jmt) ); uv=0.0
         allocate( rhom(imt,km,jmt) ); rhom=0.0
         allocate( urho(imt,km,jmt,2) ); urho=0.0
         do k=1,km
          read(io) buf; uv(:,k,:)=  buf
          read(io) buf; rhom(:,k,:)=  buf
          read(io) buf; urho(:,k,:,1)=  buf
          read(io) buf; urho(:,k,:,2)=  buf
         enddo

        elseif (name(1:6) == 'en_ext') then

         print*,' found extended energy conversion terms '
         energy_ext=.true.
         allocate( wrho(imt,km,jmt) ); wrho=0.0
         allocate( pressm(imt,km,jmt) ); pressm=0.0
         allocate( upress(imt,km,jmt,3) ); upress=0.0
         do k=1,km
          read(io) buf; wrho(:,k,:)=  buf
          read(io) buf; pressm(:,k,:)=  buf
          read(io) buf; upress(:,k,:,1)=  buf
          read(io) buf; upress(:,k,:,2)=  buf
          read(io) buf; upress(:,k,:,3)=  buf
         enddo

        elseif (name(1:10) == 'var_fluxes') then
         print*,' found fluxes of variance '
         var_fluxes = .true.
         allocate(       mean_adv(imt        ,km,jmt        ,3) )
#ifdef large_fluxfile
         allocate(    mean_tracer(imt        ,km,jmt        ,3,nt+1) )
         allocate(     var_tracer(imt        ,km,jmt        ,3,nt+1) )
         allocate( mean_eddy_flux(imt        ,km,jmt        ,3,nt+1) )
         allocate(  mean_var_flux(imt        ,km,jmt        ,3,nt+1) )
#else
         allocate(    mean_tracer(imt,km,jmt,3,nt+1:nt+1) )
         allocate(     var_tracer(imt,km,jmt,3,nt+1:nt+1) )
         allocate( mean_eddy_flux(imt,km,jmt,3,nt+1:nt+1) )
         allocate(  mean_var_flux(imt,km,jmt,3,nt+1:nt+1) )
#endif
         mean_adv=0.0; mean_tracer=0.0; var_tracer=0.0
         mean_eddy_flux=0.0; mean_var_flux = 0.0; 
#ifdef large_fluxfile
         allocate(          theta(imt        ,km,jmt        ,3,nt+1) )
         allocate(      eddy_flux(imt        ,km,jmt        ,3,nt+1) )
         theta=0; eddy_flux=0.
#endif
         do k=1,km
          print*,'k=',k,'/',km
#ifdef large_fluxfile
          do n=1,3  ! sum(u)
           read(io) buf ; mean_adv(:,k,:,n)= buf
          enddo
          do n=1,nt+1 ! sum(b)
           do m=1,3 
            read(io) buf ; mean_tracer(:,k,:,m,n)=buf
           enddo
          enddo
          do n=1,nt+1 ! sum(b^2)
           do m=1,3 
            read(io) buf ; var_tracer(:,k,:,m,n)=buf
           enddo
          enddo
          do n=1,nt+1  ! sum(ub)
           do m=1,3 
            read(io) buf ; mean_eddy_flux(:,k,:,m,n)=buf
           enddo
          enddo
          do n=1,nt+1  ! sum(ub^2)
           do m=1,3 
            read(io) buf ; mean_var_flux(:,k,:,m,n)=buf
           enddo
          enddo
#else
          do n=1,3  ! sum(u)
           read(io) buf ; mean_adv(:,k,:,n)= buf
          enddo
          do m=1,3*nt; read(io) ; enddo
          do m=1,3 
            read(io) buf ; mean_tracer(:,k,:,m,nt+1)=buf
          enddo
          do m=1,3*nt; read(io) ; enddo
          do m=1,3 
            read(io) buf ; var_tracer(:,k,:,m,nt+1)=buf
          enddo
          do m=1,3*nt; read(io) ; enddo
          do m=1,3 
            read(io) buf ; mean_eddy_flux(:,k,:,m,nt+1)=buf
          enddo
          do m=1,3*nt; read(io) ; enddo
          do m=1,3 
            read(io) buf ; mean_var_flux(:,k,:,m,nt+1)=buf
          enddo
#endif
         enddo

        elseif (name(1:6) == 'sshvar' ) then

         print*,' found variance of SSH'
         read(io) var_ssh
         enable_freesurf= .true.

        elseif (name(1:7) == 'sshmean' ) then

         print*,' found mean of SSH'
         read(io) mean_ssh
         enable_freesurf= .true.

        else  

         print*,' WARNING : cannot read variable ',name,' in file'
         print*,' skipping this variable ...'
         do n=1,nr
          read(io)
         enddo

        endif

       endif

      enddo
      close(io)

      end subroutine read_variances










      subroutine write_variances_eke(outfile,year,month,day,domain)
      use cv_var_module
      implicit none
      character (len=80) :: outfile,name,unit,varname,varlname
      integer :: year,month,day,domain
#include "netcdf.inc"
      integer ncid,iret,i,j,k,dims(4),start(4),count(4),m
      real (kind=r_knd) :: spval=-9.9e12, tt
      real (kind=r_knd) :: dx,dy,rhox,rhoy,ux,vy,uy,vx,rhodz
      real, parameter :: pi       = 
     &        3.14159265358979323846264338327950588
      real, parameter :: radius   = 6370.0e3
      real :: degtm   = radius/180*pi
      real (kind=r_knd) :: buf1(imt,jmt),buf2(imt,jmt),buf3(imt,jmt)
      real (kind=r_knd), allocatable :: mean_wt(:,:,:),mean_wu(:,:,:)
      real (kind=r_knd), allocatable :: eke(:,:,:),mke(:,:,:)
      real (kind=r_knd) :: dxt,dxu,dxum,dyt,dyu,dyum,cst,csu,cstp
      real (kind=r_knd) :: due,duw,dun,dus,dzt,dx1,dx2,dz

      integer :: lon_tdim,lon_tid
      integer :: lon_udim,lon_uid
      integer :: lat_tdim,lat_tid
      integer :: lat_udim,lat_uid
      integer :: depth_tdim,depth_tid
      integer :: depth_wdim,depth_wid
      integer :: itimedim,itimeid
      integer :: ekeid,mkeid,sshid,msshid,n,len
      integer :: t2id,t4id,T2xid,T2yid,t1id,t3id
      integer :: ubid(nt+1),vbid(nt+1),wbid(nt+1)
      integer :: ub2id(nt+1),vb2id(nt+1),wb2id(nt+1)
      integer :: ubbid(nt+1),vbbid(nt+1),wbbid(nt+1)
      integer :: tmid(nt+1),tvid(nt+1)
      integer :: 
     &days_in_month(12) = (/31,28,31,30,31,30,31,31,30,31,30,31/)

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

      if (enable_monthly_variances) then
        tt=tt-days_in_month(month)/2.
      elseif (enable_seasonal_variances) then
        tt=tt-(days_in_month(month)    ! month should be always ge 3
     &        +days_in_month(month-1)
     &        +days_in_month(month-2))/2.
      endif

      print*,' writing to file ',outfile(1:len_trim(outfile))

      ncid = nccre (outfile, NCCLOB, iret)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
c     dimensions
      lon_tdim  = ncddef(ncid, 'Longitude_t', imt, iret)
      Lon_udim  = ncddef(ncid, 'Longitude_u', imt, iret)
      Lat_tdim  = ncddef(ncid, 'Latitude_t',  jmt, iret)
      Lat_udim  = ncddef(ncid, 'Latitude_u',  jmt, iret)
      depth_wdim = ncddef(ncid, 'depth_w',  km, iret)
      depth_tdim = ncddef(ncid, 'depth_t',  km, iret)
      iTimedim  = ncddef(ncid, 'Time', 1, iret)
c     grid variables
      dims(1)  = Lon_tdim
      Lon_tid  = ncvdef (ncid,'Longitude_t',NCFLOAT,1,dims,iret)
      dims(1)  = Lon_udim
      Lon_uid  = ncvdef (ncid,'Longitude_u',NCFLOAT,1,dims,iret)
      dims(1)  = Lat_tdim
      Lat_tid  = ncvdef (ncid,'Latitude_t', NCFLOAT,1,dims,iret)
      dims(1)  = Lat_udim
      Lat_uid  = ncvdef (ncid,'Latitude_u', NCFLOAT,1,dims,iret)
      dims(1)  = iTimedim
      iTimeid   = ncvdef(ncid,'Time',       NCFLOAT,1,dims,iret)
      dims(1)  = depth_wdim
      depth_wid = ncvdef (ncid,'depth_w', NCFLOAT,1,dims,iret)
      dims(1)  = depth_tdim
      depth_tid = ncvdef (ncid,'depth_t', NCFLOAT,1,dims,iret)
c     attributes of the grid
      name = 'Longitude on T grid     '; unit = 'degrees_W       '
      call ncaptc(ncid, Lon_tid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lon_tid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Longitude on U grid     '; unit = 'degrees_W       '
      call ncaptc(ncid, Lon_uid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lon_uid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Latitude on T grid      '; unit = 'degrees_N       '
      call ncaptc(ncid, Lat_tid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lat_tid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Latitude on U grid      '; unit = 'degrees_N       '
      call ncaptc(ncid, Lat_uid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lat_uid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Depth of T Grid points  '; unit = 'm '
      call ncaptc(ncid, depth_tid, 'long_name',NCCHAR, 24, name, iret) 
      call ncaptc(ncid, depth_tid, 'units',    NCCHAR, 16, unit, iret) 
      name = 'Depth of W Grid points  '; unit = 'm '
      call ncaptc(ncid, depth_wid, 'long_name', NCCHAR,24, name, iret) 
      call ncaptc(ncid, depth_wid, 'units',     NCCHAR,16, unit, iret) 
      name = 'Time                    '; unit = 'days            '
      call ncaptc(ncid, iTimeid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, iTimeid, 'units',     NCCHAR, 16, unit, iret) 
      call ncaptc(ncid, iTimeid,'time_origin',NCCHAR, 20,
     &'31-DEC-1899 00:00:00', iret)
c     variables
      dims=(/lon_udim,lat_udim,depth_tdim,itimedim/)
      ekeid=ncvdef (ncid,'eke', NCFLOAT,4,dims,iret)
      mkeid=ncvdef (ncid,'mke', NCFLOAT,4,dims,iret)
      name = 'Eddy kinetic energy'; unit = 'm^/s^2'
      call dvcdf_here(ncid,ekeid,name,len_trim(name),unit,16,spval)
      name = 'Mean kinetic energy'; unit = 'm^2/s^2'
      call dvcdf_here(ncid,mkeid,name,len_trim(name),unit,16,spval)

      if (enable_freesurf) then
       dims=(/lon_tdim,lat_tdim,itimedim,1/)
       sshid=ncvdef (ncid,'ssh_rms', NCFLOAT,3,dims,iret)
       msshid=ncvdef (ncid,'ssh_mean', NCFLOAT,3,dims,iret)
c      attributes of variables
       name = 'Rms of Sea Surface height'; unit = 'm'
       call dvcdf_here(ncid,sshid,name,len_Trim(name),unit,16,spval)
       name = 'Mean Sea Surface height'; unit = 'm'
       call dvcdf_here(ncid,msshid,name,len_Trim(name),unit,16,spval)
      endif


      iret= nf_put_att_int(ncid,nf_global,
     &   'Number_of_averaged_time_steps ',nf_int,1,counter)
        if (iret.ne.0) print*,nf_strerror(iret)

      write( name,'("SFLAME sub domain # ",i3)') domain
      if (enable_monthly_variances) 
     &   name=name(1:len_trim(name))//' monthly means'
      if (enable_seasonal_variances) 
     &   name=name(1:len_trim(name))//' seasonal means'
      len=len_trim(name)
      iret= nf_put_att_text(ncid,nf_global,'history',len,name)

      call ncendf(ncid, iret)

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

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

c     convert to SI units
      mean_u=mean_u/100.   
      var_u=var_u/100./100. 
      mean_ssh=mean_ssh/100.
      var_ssh=var_ssh/100./100.

c
c    get means and variances:
c     sum(T*) = sum(T^2) - N T_bar^2
c     N T_bar = sum(T_n)
c     1/N [sum(T_n)]^2 = N T_bar^2
c
      do n=1,2
       do j=1,jmt
        do k=1,km
         do i=1,imt
          var_u(i,k,j,n)= var_u(i,k,j,n)-mean_u(i,k,j,n)**2/counter 
          var_u(i,k,j,n)= var_u(i,k,j,n)/counter
          mean_u(i,k,j,n)= mean_u(i,k,j,n)/counter
         enddo
        enddo
       enddo
      enddo

c
c   store energy terms  in var_u and mean_u
c
      allocate(eke(imt,km,jmt) , mke(imt,km,jmt) )
      do j=1,jmt
       do k=1,km
        do i=1,imt
         eke(i,k,j)=(var_u(i,k,j,1)+var_u(i,k,j,2))/2.
         mke(i,k,j)=(mean_u(i,k,j,1)**2+mean_u(i,k,j,2)**2)/2.
        enddo
       enddo
      enddo

c
c   care about external mode
c
      do j=1,jmt
       do i=1,imt
        var_ssh(i,j)= var_ssh(i,j)-mean_ssh(i,j)**2/counter 
        var_ssh(i,j)= sqrt(var_ssh(i,j)/counter)
        mean_ssh(i,j)= mean_ssh(i,j)/counter
       enddo
      enddo

c     apply land mask and write to file
      do k=1,km
         where( kmu <k ) eke(:,k,:) = spval
         where( kmu <k ) mke(:,k,:) = spval
         start=(/1,1,k,1/); count=(/imt,jmt,1,1/)
#ifdef real_8
         iret=nf_put_vara_double(ncid,mkeid,start,count,mke(:,k,:))
         iret=nf_put_vara_double (ncid,ekeid,start,count,eke(:,k,:))
#else
         iret=nf_put_vara_real(ncid,mkeid,start,count,mke(:,k,:))
         iret=nf_put_vara_real(ncid,ekeid,start,count,eke(:,k,:))
#endif
      enddo


      if (enable_freesurf) then
       where( kmt == 0) mean_ssh = spval
       where( kmt == 0) var_ssh = spval
       start=(/1,1,1,1/); count=(/imt,jmt,1,1/)
#ifdef real_8
       iret= nf_put_vara_double(ncid,sshid,start,count,var_ssh)
       iret= nf_put_vara_double(ncid,msshid,start,count,mean_ssh)
#else
       iret= nf_put_vara_real (ncid,sshid,start,count,var_ssh)
       iret= nf_put_vara_real (ncid,msshid,start,count,mean_ssh)
#endif
      endif
      call ncclos (ncid, iret)
      deallocate(eke,mke)
      end subroutine write_variances_eke




      subroutine write_variances_energy(outfile,year,month,day,domain)
      use cv_var_module
      implicit none
      character (len=80) :: outfile,name,unit,varname,varlname
      integer :: year,month,day,domain
#include "netcdf.inc"
      integer ncid,iret,i,j,k,dims(4),start(4),count(4),m
      real (kind=r_knd) :: spval=-9.9e12, tt
      real (kind=r_knd) :: dx,dy,rhox,rhoy,ux,vy,uy,vx,rhodz
      real, parameter :: pi       = 
     &        3.14159265358979323846264338327950588
      real, parameter :: radius   = 6370.0e3
      real :: degtm   = radius/180*pi
      real (kind=r_knd) :: buf1(imt,jmt),buf2(imt,jmt),buf3(imt,jmt)
      real (kind=r_knd), allocatable :: mean_wt(:,:,:),mean_wu(:,:,:)
      real (kind=r_knd) :: dxt,dxu,dxum,dyt,dyu,dyum,cst,csu,cstp
      real (kind=r_knd) :: due,duw,dun,dus,dzt,dx1,dx2,dz

      integer :: lon_tdim,lon_tid
      integer :: lon_udim,lon_uid
      integer :: lat_tdim,lat_tid
      integer :: lat_udim,lat_uid
      integer :: depth_tdim,depth_tid
      integer :: depth_wdim,depth_wid
      integer :: itimedim,itimeid
      integer :: ekeid,mkeid,sshid,msshid,n,len
      integer :: t2id,t4id,T2xid,T2yid,t1id,t3id
      integer :: ubid(nt+1),vbid(nt+1),wbid(nt+1)
      integer :: ub2id(nt+1),vb2id(nt+1),wb2id(nt+1)
      integer :: ubbid(nt+1),vbbid(nt+1),wbbid(nt+1)
      integer :: tmid(nt+1),tvid(nt+1)
      integer :: 
     &days_in_month(12) = (/31,28,31,30,31,30,31,31,30,31,30,31/)

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

      if (enable_monthly_variances) then
        tt=tt-days_in_month(month)/2.
      elseif (enable_seasonal_variances) then
        tt=tt-(days_in_month(month)    ! month should be always ge 3
     &        +days_in_month(month-1)
     &        +days_in_month(month-2))/2.
      endif

      print*,' writing to file ',outfile(1:len_trim(outfile))

      ncid = nccre (outfile, NCCLOB, iret)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
c     dimensions
      lon_tdim  = ncddef(ncid, 'Longitude_t', imt, iret)
      Lon_udim  = ncddef(ncid, 'Longitude_u', imt, iret)
      Lat_tdim  = ncddef(ncid, 'Latitude_t',  jmt, iret)
      Lat_udim  = ncddef(ncid, 'Latitude_u',  jmt, iret)
      depth_wdim = ncddef(ncid, 'depth_w',  km, iret)
      depth_tdim = ncddef(ncid, 'depth_t',  km, iret)
      iTimedim  = ncddef(ncid, 'Time', 1, iret)
c     grid variables
      dims(1)  = Lon_tdim
      Lon_tid  = ncvdef (ncid,'Longitude_t',NCFLOAT,1,dims,iret)
      dims(1)  = Lon_udim
      Lon_uid  = ncvdef (ncid,'Longitude_u',NCFLOAT,1,dims,iret)
      dims(1)  = Lat_tdim
      Lat_tid  = ncvdef (ncid,'Latitude_t', NCFLOAT,1,dims,iret)
      dims(1)  = Lat_udim
      Lat_uid  = ncvdef (ncid,'Latitude_u', NCFLOAT,1,dims,iret)
      dims(1)  = iTimedim
      iTimeid   = ncvdef(ncid,'Time',       NCFLOAT,1,dims,iret)
      dims(1)  = depth_wdim
      depth_wid = ncvdef (ncid,'depth_w', NCFLOAT,1,dims,iret)
      dims(1)  = depth_tdim
      depth_tid = ncvdef (ncid,'depth_t', NCFLOAT,1,dims,iret)
c     attributes of the grid
      name = 'Longitude on T grid     '; unit = 'degrees_W       '
      call ncaptc(ncid, Lon_tid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lon_tid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Longitude on U grid     '; unit = 'degrees_W       '
      call ncaptc(ncid, Lon_uid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lon_uid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Latitude on T grid      '; unit = 'degrees_N       '
      call ncaptc(ncid, Lat_tid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lat_tid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Latitude on U grid      '; unit = 'degrees_N       '
      call ncaptc(ncid, Lat_uid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lat_uid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Depth of T Grid points  '; unit = 'm '
      call ncaptc(ncid, depth_tid, 'long_name',NCCHAR, 24, name, iret) 
      call ncaptc(ncid, depth_tid, 'units',    NCCHAR, 16, unit, iret) 
      name = 'Depth of W Grid points  '; unit = 'm '
      call ncaptc(ncid, depth_wid, 'long_name', NCCHAR,24, name, iret) 
      call ncaptc(ncid, depth_wid, 'units',     NCCHAR,16, unit, iret) 
      name = 'Time                    '; unit = 'days            '
      call ncaptc(ncid, iTimeid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, iTimeid, 'units',     NCCHAR, 16, unit, iret) 
      call ncaptc(ncid, iTimeid,'time_origin',NCCHAR, 20,
     &'31-DEC-1899 00:00:00', iret)
c     variables


       dims=(/lon_tdim,lat_tdim,depth_tdim,itimedim/)
       t2id=ncvdef (ncid,'T2', NCFLOAT,4,dims,iret)
       name = 'Baroclinic conversion rate'; unit = '(kg/m^3)^2/s'
       call dvcdf_here(ncid,T2id,name,len_Trim(name),unit,16,spval)

       dims=(/lon_udim,lat_udim,depth_tdim,itimedim/)
       T4id=ncvdef (ncid,'T4', NCFLOAT,4,dims,iret)
       name = 'Barotropic conversion rate'; unit = '(m/s)^2/s'
       call dvcdf_here(ncid,T4id,name,len_Trim(name),unit,16,spval)

       dims=(/lon_tdim,lat_tdim,depth_tdim,itimedim/)
       t2xid=ncvdef (ncid,'T2x', NCFLOAT,4,dims,iret)
       name = 'Zonal baroclinic conversion rate'; unit = '(kg/m^3)^2/s'
       call dvcdf_here(ncid,T2xid,name,len_Trim(name),unit,16,spval)

      if (energy_ext) then
       dims=(/lon_tdim,lat_tdim,depth_wdim,itimedim/)
       T3id=ncvdef (ncid,'T3', NCFLOAT,4,dims,iret)
       name = 'T3'; unit = ' '
       call dvcdf_here(ncid,T3id,name,len_Trim(name),unit,16,spval)
       dims=(/lon_tdim,lat_tdim,depth_wdim,itimedim/)
       T1id=ncvdef (ncid,'T1', NCFLOAT,4,dims,iret)
       name = 'T1'; unit = ' '
       call dvcdf_here(ncid,T1id,name,len_Trim(name),unit,16,spval)
      endif

      iret= nf_put_att_int(ncid,nf_global,
     &   'Number_of_averaged_time_steps ',nf_int,1,counter)
        if (iret.ne.0) print*,nf_strerror(iret)

      write( name,'("SFLAME sub domain # ",i3)') domain
      if (enable_monthly_variances) 
     &   name=name(1:len_trim(name))//' monthly means'
      if (enable_seasonal_variances) 
     &   name=name(1:len_trim(name))//' seasonal means'
      len=len_trim(name)
      iret= nf_put_att_text(ncid,nf_global,'history',len,name)

      call ncendf(ncid, iret)

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


c     convert to SI units

       urho=urho/100.0 * 1000.0
       uv = uv/100.0/100.0
       rhom = rhom *1000.0

      if (energy_ext) then
        wrho=wrho/100.0*1000.0
        pressm=pressm /100.0 /100.0 *1000.0
        upress=upress/100.0 /100.0/100.0 *1000.0
      endif

c    get means and variances:
c     sum(T*) = sum(T^2) - N T_bar^2
c     N T_bar = sum(T_n)
c     1/N [sum(T_n)]^2 = N T_bar^2
c

       do j=1,jmt
        do k=1,km
         do i=1,imt
           rhom(i,k,j)=rhom(i,k,j)/counter
           urho(i,k,j,1)= urho(i,k,j,1) /counter
     &                  - rhom(i,k,j)*mean_u(i,k,j,1)
           urho(i,k,j,2)= urho(i,k,j,2) /counter
     &                  - rhom(i,k,j)*mean_u(i,k,j,2)
           uv(i,k,j)= uv(i,k,j) /counter
     &                  - mean_u(i,k,j,1)*mean_u(i,k,j,2)
         enddo
        enddo
       enddo

      if (energy_ext) then
c   
c    get T1 and T3
c    first calculate mean vertical velocity
       allocate(mean_wu(imt,km,jmt) ,mean_wt(imt,km,jmt) )
       mean_wu=0.; mean_wt=0.
c
c    vertical velocity on T grid
       do k=1,km
        buf1=0.0;buf2=0.0
        do j=2,jmt-1
          do i=2,imt-1
            dxt = (xu(i)-xu(i-1))*degtm
            dxu = (xt(i+1)-xt(i))*degtm
            dxum = (xt(i)-xt(i-1))*degtm
            dyt = (yu(j)-yu(j-1))*degtm
            dyu = (yt(j+1)-yt(j))*degtm
            dyum = (yt(j)-yt(j-1))*degtm
            csu  = cos( yu(j)/180*pi )
            buf1(i,j) = 0.5*(mean_u(i,k,j,1)*dyu +    ! adv_vet
     &                        mean_u(i,k,j-1,1)*dyum)/dyt
            buf2(i,j) = 0.5*(mean_u(i,k,j,2)*dxu +    ! adv_vnt
     &                        mean_u(i-1,k,j,2)*dxum)*csu/dxt
          enddo
        enddo
        do j=2,jmt
          do i=2,imt
            dxt = (xu(i)-xu(i-1))*degtm
            dyt = (yu(j)-yu(j-1))*degtm
            dzt  = (zw(k)- zw(k-1))
            cst  = cos(yt(j)/180*pi)
            mean_wt(i,k,j) =
     &                   ((buf1(i,j) - buf1(i-1,j))/dxt
     &                   +(buf2(i,j) - buf2(i,j-1))/dyt
     &                   )/cst*dzt
          enddo
        enddo
       end do

       do j=1,jmt
        do k=2,km ! mean_wt(:,0,:)=0.0
          do i=1,imt
            mean_wt(i,k,j) = mean_wt(i,k,j) + mean_wt(i,k-1,j)
          enddo
        enddo
       end do

c    and vertical velocity on the u grid
       do j=2,jmt-1
        do k=1,km
          do i=2,imt-1
            duw = (xu(i)   - xt(i))*degtm
	    due = (xt(i+1) - xu(i))*degtm
	    dun = (yt(j+1) - yu(j))*degtm
            dus = (yu(j)   - yt(j))*degtm
            dxu = (xt(i+1)-xt(i))*degtm
            dyu = (yt(j+1)-yt(j))*degtm
            csu = cos( yu(j)/180*pi )
            cst  = cos( yt(j)/180*pi )
            cstp = cos( yt(j+1)/180*pi )
	    mean_wu(i,k,j) = 1./dyu/csu/dxu*(
     &                    mean_wt(i,k,j)*duw*dus*cst
     &                  + mean_wt(i+1,k,j)*due*dus*cst
     &                  + mean_wt(i,k,j+1)*duw*dun*cstp
     &                  + mean_wt(i+1,k,j+1)*due*dun*cstp )
          enddo
        enddo
       enddo

c    now get means and variances of averaged quantities
c    store T3 in wrho and T1 in pressm

       do j=1,jmt
        do k=1,km
         do i=1,imt
           pressm(i,k,j)=pressm(i,k,j)/counter
           wrho(i,k,j)= wrho(i,k,j) /counter
     &                  - rhom(i,k,j)*mean_wt(i,k,j)
           upress(i,k,j,1)= upress(i,k,j,1) /counter
     &                  - pressm(i,k,j)*mean_u(i,k,j,1)
           upress(i,k,j,2)= upress(i,k,j,2) /counter
     &                  - pressm(i,k,j)*mean_u(i,k,j,2)
           upress(i,k,j,3)= upress(i,k,j,3) /counter
     &                  - pressm(i,k,j)*mean_wu(i,k,j)
           pressm(i,k,j)=rhom(i,k,j)*mean_wt(i,k,j)
         enddo
        enddo
       enddo

       deallocate( mean_wt, mean_wu)
      endif ! energy_ext


c
c    get T2 = u_prime rho_prime (rho_bar)_x and v_prime rho_prime (rho_bar)_y 
c    store T2 in rhom
c
       do j=1,jmt-1
        do k=1,km
         do i=1,imt-1
          dx=(xt(i+1)-xt(i))*cos(yt(j)/180*pi)*degtm
          rhox = (rhom(i+1,k,j) - rhom(i,k,j))/dx
          urho(i,k,j,1)= urho(i,k,j,1) *rhox
          dy=(yt(j+1)-yt(j))*degtm
          rhoy = (rhom(i,k,j+1) - rhom(i,k,j))/dy
          urho(i,k,j,2)= urho(i,k,j,2) *rhoy
         enddo
        enddo
       enddo
       rhom=0.0 ! here we are loosing rho_bar, it becomes T2
       do k=1,km
        buf1=urho(:,k,:,1)
        buf2=urho(:,k,:,2)
        urho(:,k,:,1:2)=0.0
        do j=2,jmt-1
         do i=2,imt-1
           urho(i,k,j,1)=(buf1(i,j)+buf1(i-1,j) )/2.0
           urho(i,k,j,2)=(buf2(i,j)+buf2(i-1,j) )/2.0
           rhom(i,k,j)= urho(i,k,j,1)+urho(i-1,k,j,2) 
         enddo
        enddo
       enddo
c
c    get T4 = 
c    u_prime^2 (u_bar)_x + v_prime^2 (v_bar)_y + 
c    u_prime v_prime ( (u_bar)_y + (v_bar)x )
c     store T4 in uv and first part of T4 in urho(:,:,:,1)
c
       do k=1,km
        buf1=0.0;buf2=0.0;
        buf3=uv(:,k,:);uv(:,k,:)=0.0
        do j=2,jmt-1
         do i=2,imt-1
          dx=(xu(i+1)-xu(i))*cos(yu(j)/180*pi)*degtm
          ux = (mean_u(i+1,k,j,1) - mean_u(i,k,j,1))/dx
          dx=(xu(i)-xu(i-1))*cos(yu(j)/180*pi)*degtm
          ux = (ux + (mean_u(i,k,j,1) - mean_u(i-1,k,j,1))/dx )/2.0

          buf1(i,j)= var_u(i,k,j,1)*ux

          dy=(yu(j+1)-yu(j))*degtm
          vy = (mean_u(i,k,j+1,2) - mean_u(i,k,j,2))/dy
          dy=(yu(j)-yu(j-1))*degtm
          vy = (vy + (mean_u(i,k,j,2) - mean_u(i,k,j-1,2))/dy )/2.0

          buf2(i,j)= var_u(i,k,j,2)*vy

          dx=(xu(i+1)-xu(i))*cos(yu(j)/180*pi)*degtm
          vx = (mean_u(i+1,k,j,2) - mean_u(i,k,j,2))/dx
          dx=(xu(i)-xu(i-1))*cos(yu(j)/180*pi)*degtm
          vx = (vx + (mean_u(i,k,j,2) - mean_u(i-1,k,j,2))/dx )/2.0

          dy=(yu(j+1)-yu(j))*degtm
          uy = (mean_u(i,k,j+1,1) - mean_u(i,k,j,1))/dy
          dy=(yu(j)-yu(j-1))*degtm
          uy = (uy + (mean_u(i,k,j,1) - mean_u(i,k,j-1,1))/dy )/2.0

          uv(i,k,j)=  buf1(i,j) + buf2(i,j) +buf3(i,j)*(vx+uy)

         enddo
        enddo
       enddo



       do k=1,km
         where( kmu <k ) uv(:,k,:) = spval
         where( kmt <k ) rhom(:,k,:) = spval
         where( kmt <k ) urho(:,k,:,1) = spval
         start=(/1,1,k,1/); count=(/imt,jmt,1,1/)
#ifdef real_8
         iret=nf_put_vara_double(ncid,T4id,start,count,uv(:,k,:))
         iret=nf_put_vara_double(ncid,T2id,start,count,rhom(:,k,:))
         iret=nf_put_vara_double(ncid,T2xid,start,count,urho(:,k,:,1))
#else
         iret=nf_put_vara_real(ncid,T4id,start,count,uv(:,k,:))
         iret=nf_put_vara_real(ncid,T2id,start,count,rhom(:,k,:))
         iret=nf_put_vara_real(ncid,T2xid,start,count,urho(:,k,:,1))
#endif
       enddo

      if (energy_ext) then
       do k=1,km
         where( kmt <k ) wrho(:,k,:) = spval
         where( kmt <k ) pressm(:,k,:) = spval
         start=(/1,1,k,1/); count=(/imt,jmt,1,1/)
#ifdef real_8
         iret=nf_put_vara_double(ncid,T1id,start,count,pressm(:,k,:))
         iret=nf_put_vara_double(ncid,T3id,start,count,wrho(:,k,:))
#else
         iret=nf_put_vara_real(ncid,T1id,start,count,pressm(:,k,:))
         iret=nf_put_vara_real(ncid,T3id,start,count,wrho(:,k,:))
#endif
       enddo
      endif


      call ncclos (ncid, iret)
      end subroutine write_variances_energy



      subroutine write_variances_fluxes(outfile,year,month,day,domain)
      use cv_var_module
      implicit none
      character (len=80) :: outfile,name,unit,varname,varlname
      integer :: year,month,day,domain
#include "netcdf.inc"
      integer ncid,iret,i,j,k,dims(4),start(4),count(4),m
      real (kind=r_knd) :: spval=-9.9e12, tt
      real (kind=r_knd) :: dx,dy,rhox,rhoy,ux,vy,uy,vx,rhodz
      real, parameter :: pi       = 
     &        3.14159265358979323846264338327950588
      real, parameter :: radius   = 6370.0e3
      real :: degtm   = radius/180*pi
      real (kind=r_knd) :: buf1(imt,jmt),buf2(imt,jmt),buf3(imt,jmt)
      real (kind=r_knd), allocatable :: mean_wt(:,:,:),mean_wu(:,:,:)
      real (kind=r_knd) :: dxt,dxu,dxum,dyt,dyu,dyum,cst,csu,cstp
      real (kind=r_knd) :: due,duw,dun,dus,dzt,dx1,dx2,dz

      integer :: lon_tdim,lon_tid
      integer :: lon_udim,lon_uid
      integer :: lat_tdim,lat_tid
      integer :: lat_udim,lat_uid
      integer :: depth_tdim,depth_tid
      integer :: depth_wdim,depth_wid
      integer :: itimedim,itimeid
      integer :: ekeid,mkeid,sshid,msshid,n,len
      integer :: t2id,t4id,T2xid,T2yid,t1id,t3id
      integer :: ubid(nt+1),vbid(nt+1),wbid(nt+1)
      integer :: ub2id(nt+1),vb2id(nt+1),wb2id(nt+1)
      integer :: ubbid(nt+1),vbbid(nt+1),wbbid(nt+1)
      integer :: tmid(nt+1),tvid(nt+1)
      integer :: 
     &days_in_month(12) = (/31,28,31,30,31,30,31,31,30,31,30,31/)

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

      if (enable_monthly_variances) then
        tt=tt-days_in_month(month)/2.
      elseif (enable_seasonal_variances) then
        tt=tt-(days_in_month(month)    ! month should be always ge 3
     &        +days_in_month(month-1)
     &        +days_in_month(month-2))/2.
      endif

      print*,' writing to file ',outfile(1:len_trim(outfile))

      ncid = nccre (outfile, NCCLOB, iret)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
c     dimensions
      lon_tdim  = ncddef(ncid, 'Longitude_t', imt, iret)
      Lon_udim  = ncddef(ncid, 'Longitude_u', imt, iret)
      Lat_tdim  = ncddef(ncid, 'Latitude_t',  jmt, iret)
      Lat_udim  = ncddef(ncid, 'Latitude_u',  jmt, iret)
      depth_wdim = ncddef(ncid, 'depth_w',  km, iret)
      depth_tdim = ncddef(ncid, 'depth_t',  km, iret)
      iTimedim  = ncddef(ncid, 'Time', 1, iret)
c     grid variables
      dims(1)  = Lon_tdim
      Lon_tid  = ncvdef (ncid,'Longitude_t',NCFLOAT,1,dims,iret)
      dims(1)  = Lon_udim
      Lon_uid  = ncvdef (ncid,'Longitude_u',NCFLOAT,1,dims,iret)
      dims(1)  = Lat_tdim
      Lat_tid  = ncvdef (ncid,'Latitude_t', NCFLOAT,1,dims,iret)
      dims(1)  = Lat_udim
      Lat_uid  = ncvdef (ncid,'Latitude_u', NCFLOAT,1,dims,iret)
      dims(1)  = iTimedim
      iTimeid   = ncvdef(ncid,'Time',       NCFLOAT,1,dims,iret)
      dims(1)  = depth_wdim
      depth_wid = ncvdef (ncid,'depth_w', NCFLOAT,1,dims,iret)
      dims(1)  = depth_tdim
      depth_tid = ncvdef (ncid,'depth_t', NCFLOAT,1,dims,iret)
c     attributes of the grid
      name = 'Longitude on T grid     '; unit = 'degrees_W       '
      call ncaptc(ncid, Lon_tid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lon_tid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Longitude on U grid     '; unit = 'degrees_W       '
      call ncaptc(ncid, Lon_uid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lon_uid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Latitude on T grid      '; unit = 'degrees_N       '
      call ncaptc(ncid, Lat_tid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lat_tid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Latitude on U grid      '; unit = 'degrees_N       '
      call ncaptc(ncid, Lat_uid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lat_uid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Depth of T Grid points  '; unit = 'm '
      call ncaptc(ncid, depth_tid, 'long_name',NCCHAR, 24, name, iret) 
      call ncaptc(ncid, depth_tid, 'units',    NCCHAR, 16, unit, iret) 
      name = 'Depth of W Grid points  '; unit = 'm '
      call ncaptc(ncid, depth_wid, 'long_name', NCCHAR,24, name, iret) 
      call ncaptc(ncid, depth_wid, 'units',     NCCHAR,16, unit, iret) 
      name = 'Time                    '; unit = 'days            '
      call ncaptc(ncid, iTimeid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, iTimeid, 'units',     NCCHAR, 16, unit, iret) 
      call ncaptc(ncid, iTimeid,'time_origin',NCCHAR, 20,
     &'31-DEC-1899 00:00:00', iret)

#ifdef large_fluxfile
       do n=1,nt+1
#else
       do n=nt+1,nt+1
#endif
        if (n==1) then
         varname='T'        
         varlname = 'temperature'
         unit = 'mK/s'
        elseif (n==2) then
         varname='S'        
         varlname = 'salinity'
         unit = 'm/s'
        elseif (n==nt+1) then
         varname='sig0'        
         varlname = 'pot. density'
         unit = 'kg/m^2/s'
        else
         write(varname,'("tr_",i2)') n
         call replace_space_zero(varname)
         write(varlname,'("tracer ",i2)') n
         unit = 'm/s '
        endif

#ifdef define_in_fluxfile_eddyflux
        dims=(/lon_tdim,lat_tdim,depth_tdim,itimedim/)
        tmid(n)=ncvdef (ncid,'mean_'//varname(1:len_trim(varname))
     &                 , NCFLOAT,4,dims,iret)
        name = 'Mean '//varlname(1:len_trim(varlname)); 
        call dvcdf_here(ncid,tmid(n),name,len_Trim(name),unit,16,spval)

        dims=(/lon_udim,lat_tdim,depth_tdim,itimedim/)
        ubid(n)=ncvdef (ncid,'u'//varname(1:len_trim(varname))
     &                 , NCFLOAT,4,dims,iret)
        name = 'Zonal eddy flux of '//varlname(1:len_trim(varlname)); 
        call dvcdf_here(ncid,ubid(n),name,len_Trim(name),unit,16,spval)

        dims=(/lon_tdim,lat_udim,depth_tdim,itimedim/)
        vbid(n)=ncvdef (ncid,'v'//varname(1:len_trim(varname))
     &                 , NCFLOAT,4,dims,iret)
        name='Meridional eddy flux of '//varlname(1:len_trim(varlname))
        call dvcdf_here(ncid,vbid(n),name,len_Trim(name),unit,16,spval)

        dims=(/lon_tdim,lat_tdim,depth_wdim,itimedim/)
        wbid(n)=ncvdef (ncid,'w'//varname(1:len_trim(varname))
     &                 , NCFLOAT,4,dims,iret)
        name='Vertical eddy flux of '//varlname(1:len_trim(varlname))
        call dvcdf_here(ncid,wbid(n),name,len_Trim(name),unit,16,spval)
#endif

#ifdef define_in_fluxfile_varianceflux
        dims=(/lon_tdim,lat_tdim,depth_tdim,itimedim/)
        tvid(n)=ncvdef (ncid,'var_'//varname(1:len_trim(varname))
     &                 , NCFLOAT,4,dims,iret)
        name = 'Variance of '//varlname(1:len_trim(varlname)); 
        call dvcdf_here(ncid,tvid(n),name,len_Trim(name),unit,16,spval)

        dims=(/lon_udim,lat_tdim,depth_tdim,itimedim/)
        ubbid(n)=ncvdef (ncid,'u'//varname(1:len_trim(varname))//'2'
     &                 , NCFLOAT,4,dims,iret)
        name = 'Zonal variance flux of '//varlname(1:len_trim(varlname)); 
        call dvcdf_here(ncid,ubbid(n),name,len_Trim(name),unit,16,spval)

        dims=(/lon_tdim,lat_udim,depth_tdim,itimedim/)
        vbbid(n)=ncvdef (ncid,'v'//varname(1:len_trim(varname))//'2'
     &                 , NCFLOAT,4,dims,iret)
        name='Meridional variance flux of '//
     &                     varlname(1:len_trim(varlname))
        call dvcdf_here(ncid,vbbid(n),name,len_Trim(name),unit,16,spval)

        dims=(/lon_tdim,lat_tdim,depth_wdim,itimedim/)
        wbbid(n)=ncvdef (ncid,'w'//varname(1:len_trim(varname))//'2'
     &                 , NCFLOAT,4,dims,iret)
        name='Vertical variance flux of '//
     &                    varlname(1:len_trim(varlname))
        call dvcdf_here(ncid,wbbid(n),name,len_Trim(name),unit,16,spval)
#endif

#ifdef large_fluxfile

        ub2id(n)=ncvdef (ncid,'u2'//varname(1:len_trim(varname))
     &                 , NCFLOAT,4,dims,iret)
        name = 'Zonal eddy flux 2 of '//varlname(1:len_trim(varlname)); 
        call dvcdf_here(ncid,ub2id(n),name,len_Trim(name),unit,16,spval)

        vb2id(n)=ncvdef (ncid,'v2'//varname(1:len_trim(varname))
     &                 , NCFLOAT,4,dims,iret)
        name='Meridional eddy flux 2 of '//
     &           varlname(1:len_trim(varlname))
        call dvcdf_here(ncid,vb2id(n),name,len_Trim(name),unit,16,spval)

        wb2id(n)=ncvdef (ncid,'w2'//varname(1:len_trim(varname))
     &                 , NCFLOAT,4,dims,iret)
        name='Vertical eddy flux 2 of '//varlname(1:len_trim(varlname))
        call dvcdf_here(ncid,wb2id(n),name,len_Trim(name),unit,16,spval)
#endif

       enddo

      iret= nf_put_att_int(ncid,nf_global,
     &   'Number_of_averaged_time_steps ',nf_int,1,counter)
        if (iret.ne.0) print*,nf_strerror(iret)

      write( name,'("SFLAME sub domain # ",i3)') domain
      if (enable_monthly_variances) 
     &   name=name(1:len_trim(name))//' monthly means'
      if (enable_seasonal_variances) 
     &   name=name(1:len_trim(name))//' seasonal means'
      len=len_trim(name)
      iret= nf_put_att_text(ncid,nf_global,'history',len,name)

      call ncendf(ncid, iret)

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


         mean_adv=mean_adv/100;  
         mean_eddy_flux=mean_eddy_flux/100; 
         mean_var_flux = mean_var_flux/100;
c
c    get means and variances:
c     sum(T*) = sum(T^2) - N T_bar^2
c     N T_bar = sum(T_n)
c     1/N [sum(T_n)]^2 = N T_bar^2
c

c
c    get the mean quantities and fluxes
c
c sum(u_prime b_prime)/N = sum(ub)/N - u_bar b_bar
c  u_bar sum(b_prime^2)/N + sum( u_prime b_prime^2)/N = 
c       sum(ub^2)/N - u_bar b_bar^2 - 2 b_bar sum(u_prime b_prime)/N
        mean_adv=mean_adv/counter
        mean_tracer=mean_tracer/counter
c        var_tracer=var_tracer/counter
c     sum(T*) = sum(T^2) - N T_bar^2
c     N T_bar = sum(T_n)
c     1/N [sum(T_n)]^2 = N T_bar^2
#ifdef large_fluxfile
        do n=1,nt+1
#else
        do n=nt+1,nt+1
#endif
        do m=1,3
        do j=1,jmt
         do k=1,km
          do i=1,imt
          var_tracer(i,k,j,m,n)= (var_tracer(i,k,j,m,n)-
     &               mean_tracer(i,k,j,m,n)**2*counter )/counter
           mean_eddy_flux(i,k,j,m,n) = mean_eddy_flux(i,k,j,m,n)/counter 
     &           - mean_adv(i,k,j,m)*mean_tracer(i,k,j,m,n)
           mean_var_flux(i,k,j,m,n) = mean_var_flux(i,k,j,m,n)/counter 
     &           - mean_adv(i,k,j,m)*mean_tracer(i,k,j,m,n)**2
     &           - 2*mean_tracer(i,k,j,m,n)*mean_eddy_flux(i,k,j,m,n)
          enddo
         enddo
        enddo
        enddo
        enddo

#ifdef large_fluxfile
c
c     now compute the gauge vector potential for the eddy fluxes
c
        do n=1,nt+1

         do j=1,jmt-1
          do k=1,km-1
           do i=2,imt
c
c     theta_1 = (G_2 b_z - G_3 b_y) /( b_x^2+b_y^2+b_z^2)
c     with G= bar( u phi_prime )
c      is on grid:  T grid in x,  U grid in y, W Grid in z
c
            dx1=(xu(i)-xu(i-1))*cos(yt(j)/180*pi)*degtm
            dx2=(xu(i)-xu(i-1))*cos(yt(j+1)/180*pi)*degtm
            dy=(yt(j+1)-yt(j))*degtm
            dz  = zt(k+1)- zt(k)
            rhodz= 0.25*(
     &        ((mean_tracer(i,k,j,1,n)-mean_tracer(i-1,k,j,1,n))/dx1)**2
     &   +((mean_tracer(i,k,j+1,1,n)-mean_tracer(i-1,k,j+1,1,n))/dx2)**2
     &  +((mean_tracer(i,k+1,j,1,n)-mean_tracer(i-1,k+1,j,1,n))/dx1)**2+
     &((mean_tracer(i,k+1,j+1,1,n)-mean_tracer(i-1,k+1,j+1,1,n))/dx2)**2
     &            )
     &       +((mean_tracer(i,k,j,2,n)-mean_tracer(i,k+1,j,2,n))/dz)**2
     &       +((mean_tracer(i,k,j+1,3,n)-mean_tracer(i,k,j,3,n))/dy)**2
            if (rhodz==0.0) rhodz=1e-10
            theta(i,k,j,1,n)= (
     &       0.5*(mean_var_flux(i,k,j,2,n)+mean_var_flux(i,k+1,j,2,n))
     &       *(mean_tracer(i,k,j,2,n)-mean_tracer(i,k+1,j,2,n))/dz
     &      -0.5*(mean_var_flux(i,k,j,3,n)+mean_var_flux(i,k,j+1,3,n))
     &       *(mean_tracer(i,k,j+1,3,n)-mean_tracer(i,k,j,3,n))/dy )
     &        /rhodz
           enddo
          enddo
         enddo

         do j=2,jmt
          do k=1,km-1
           do i=1,imt-1
c
c     theta_2 = (G_3 b_x - G_1 b_z) /( b_x^2+b_y^2+b_z^2)
c      is on grid:  U grid in x,  T grid in y, W Grid in z
c
            dx=(xt(i+1)-xt(i))*cos(yt(j)/180*pi)*degtm
            dz=zt(k+1)-zt(k)
            dy = (yu(j)-yu(j-1))*degtm
            rhodz=
     &        ((mean_tracer(i+1,k,j,3,n)-mean_tracer(i,k,j,1,n))/dx)**2
     &       +((mean_tracer(i,k,j,1,n)-mean_tracer(i,k+1,j,1,n))/dz)**2
     &       +0.25*(
     &       +((mean_tracer(i,k,j,2,n)-mean_tracer(i,k,j-1,2,n))/dy)**2
     &   +((mean_tracer(i+1,k,j,2,n)-mean_tracer(i+1,k,j-1,2,n))/dy)**2
     & +((mean_tracer(i,k+1,j,2,n)-mean_tracer(i,k+1,j-1,2,n))/dy)**2+
     &((mean_tracer(i+1,k+1,j,2,n)-mean_tracer(i+1,k+1,j-1,2,n))/dy)**2
     &            )
            if (rhodz==0.0) rhodz=1e-10
            theta(i,k,j,2,n)= (
     &       0.5*(mean_var_flux(i,k,j,3,n)+mean_var_flux(i+1,k,j,3,n))
     &       *(mean_tracer(i+1,k,j,3,n)-mean_tracer(i,k,j,3,n))/dx
     &      -0.5*(mean_var_flux(i,k,j,1,n)+mean_var_flux(i,k+1,j,1,n))
     &       *(mean_tracer(i,k,j,1,n)-mean_tracer(i,k+1,j,1,n))/dz )
     &        /rhodz
           enddo
          enddo
         enddo

         do j=1,jmt-1
          do k=2,km
           do i=1,imt-1
c
c     theta_3 = (G_1 b_y - G_2 b_x) /( b_x^2+b_y^2+b_z^2)
c      is on grid:  U grid in x,  U grid in y, T Grid in z
c
            dy = (yt(j+1)-yt(j))*degtm
            dx = (xt(i+1)-xt(i))*cos(yu(j)/180*pi)*degtm
            dz = zw(k)-zw(k-1)
            rhodz=
     &        ((mean_tracer(i+1,k,j,2,n)-mean_tracer(i,k,j,2,n))/dx)**2
     &       +((mean_tracer(i,k,j+1,1,n)-mean_tracer(i,k,j,1,n))/dy)**2
     &       +0.25*(
     &       +((mean_tracer(i,k-1,j,3,n)-mean_tracer(i,k,j,3,n))/dz)**2
     &   +((mean_tracer(i,k-1,j+1,3,n)-mean_tracer(i,k,j+1,3,n))/dz)**2
     &   +((mean_tracer(i+1,k-1,j,3,n)-mean_tracer(i+1,k,j,3,n))/dz)**2
     &+((mean_tracer(i+1,k-1,j+1,3,n)-mean_tracer(i+1,k,j+1,3,n))/dz)**2
     &            )
            if (rhodz==0.0) rhodz=1e-10
            theta(i,k,j,3,n)= (
     &       0.5*(mean_var_flux(i,k,j,1,n)+mean_var_flux(i,k,j+1,1,n))
     &       *(mean_tracer(i,k,j+1,1,n)-mean_tracer(i,k,j,1,n))/dy
     &      -0.5*(mean_var_flux(i+1,k,j,2,n)+mean_var_flux(i,k,j,2,n))
     &       *(mean_tracer(i+1,k,j,2,n)-mean_tracer(i,k,j,2,n))/dx )
     &        /rhodz
           enddo
          enddo
         enddo

c  now remove the vector potential flux from the eddy fluxes

         do j=2,jmt
          do k=2,km
           do i=2,imt
c        F_1 = u_prime b_prime - (theta_2)_z + (theta_3)_y
            dy = (yu(j)-yu(j-1))*degtm
            dz = zw(k)-zw(k-1)
            eddy_flux(i,k,j,1,n)=mean_eddy_flux(i,k,j,1,n)
     &         -(theta(i,k-1,j,2,n)-theta(i,k,j,2,n) )/dz
     &         +(theta(i,k,j,3,n)-theta(i,k,j-1,3,n) )/dy
c        F_2 = v_prime b_prime - (theta_1)_z + (theta_3)_x
            dx = (xu(i)-xu(i-1))*cos(yu(j)/180*pi)*degtm
            dz = zw(k)-zw(k-1)
            eddy_flux(i,k,j,2,n)=mean_eddy_flux(i,k,j,2,n)
     &         -(theta(i,k-1,j,1,n)-theta(i,k,j,1,n) )/dz
     &         +(theta(i,k,j,3,n)-theta(i-1,k,j,3,n) )/dx
c        F_3 = w_prime b_prime - (theta_2)_x + (theta_1)_y
            dy = (yu(j)-yu(j-1))*degtm
            dx = (xu(i)-xu(i-1))*cos(yt(j)/180*pi)*degtm
            eddy_flux(i,k,j,3,n)=mean_eddy_flux(i,k,j,3,n)
     &         -(theta(i,k,j,2,n)-theta(i-1,k,j,2,n) )/dx
     &         +(theta(i,k,j,1,n)-theta(i,k,j-1,1,n) )/dy
           enddo
          enddo
         enddo


        enddo

#endif



       do k=1,km
#ifdef large_fluxfile
         do n=1,nt+1
#else
         do n=nt+1,nt+1
#endif
          do m=1,3
           where( kmt <k ) mean_adv(:,k,:,m) = spval
           where( kmt <k ) mean_tracer(:,k,:,m,n) = spval
           where( kmt <k ) var_tracer(:,k,:,m,n) = spval
           where( kmt <k ) mean_eddy_flux(:,k,:,m,n) = spval
#ifdef large_fluxfile
           where( kmt <k ) eddy_flux(:,k,:,m,n) = spval
#endif
           where( kmt <k ) mean_var_flux(:,k,:,m,n) = spval
          enddo
         enddo

         start=(/1,1,k,1/); count=(/imt,jmt,1,1/)
#ifdef large_fluxfile
         do n=1,nt+1
#else
         do n=nt+1,nt+1
#endif
#ifdef real_8
          iret=nf_put_vara_double(ncid,tmid(n),start,count,
     &           mean_tracer(:,k,:,1,n))
          iret=nf_put_vara_double(ncid,tvid(n),start,count,
     &           var_tracer(:,k,:,1,n))
          iret=nf_put_vara_double(ncid,ubid(n),start,count,
     &           mean_eddy_flux(:,k,:,1,n))
          iret=nf_put_vara_double(ncid,vbid(n),start,count,
     &           mean_eddy_flux(:,k,:,2,n))
          iret=nf_put_vara_double(ncid,ubbid(n),start,count,
     &           mean_var_flux(:,k,:,1,n))
          iret=nf_put_vara_double(ncid,vbbid(n),start,count,
     &           mean_var_flux(:,k,:,2,n))
          iret=nf_put_vara_double(ncid,wbid(n),start,count,
     &           mean_eddy_flux(:,k,:,3,n))
          iret=nf_put_vara_double(ncid,wbbid(n),start,count,
     &           mean_var_flux(:,k,:,3,n))
#ifdef large_fluxfile
          iret=nf_put_vara_double(ncid,ub2id(n),start,count,
     &           eddy_flux(:,k,:,1,n))
          iret=nf_put_vara_double(ncid,vb2id(n),start,count,
     &           eddy_flux(:,k,:,2,n))
          iret=nf_put_vara_double(ncid,wb2id(n),start,count,
     &           eddy_flux(:,k,:,3,n))
#endif
#else
          iret=nf_put_vara_real(ncid,tmid(n),start,count,
     &           mean_tracer(:,k,:,1,n))
          iret=nf_put_vara_real(ncid,tvid(n),start,count,
     &           var_tracer(:,k,:,1,n))
          iret=nf_put_vara_real(ncid,ubid(n),start,count,
     &           mean_eddy_flux(:,k,:,1,n))
          iret=nf_put_vara_real(ncid,vbid(n),start,count,
     &           mean_eddy_flux(:,k,:,2,n))
          iret=nf_put_vara_real(ncid,ubbid(n),start,count,
     &           mean_var_flux(:,k,:,1,n))
          iret=nf_put_vara_real(ncid,vbbid(n),start,count,
     &           mean_var_flux(:,k,:,2,n))
          iret=nf_put_vara_real(ncid,wbid(n),start,count,
     &           mean_eddy_flux(:,k,:,3,n))
          iret=nf_put_vara_real(ncid,wbbid(n),start,count,
     &           mean_var_flux(:,k,:,3,n))
#ifdef large_fluxfile
          iret=nf_put_vara_real(ncid,ub2id(n),start,count,
     &           eddy_flux(:,k,:,1,n))
          iret=nf_put_vara_real(ncid,vb2id(n),start,count,
     &           eddy_flux(:,k,:,2,n))
          iret=nf_put_vara_real(ncid,wb2id(n),start,count,
     &           eddy_flux(:,k,:,3,n))
#endif
#endif
         enddo
       enddo

      call ncclos (ncid, iret)
      end subroutine write_variances_fluxes



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


      subroutine barrier
      end

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

#ifdef notdef
      subroutine write_variances(outfile,year,month,day,domain)
      use cv_var_module
      implicit none
      character (len=80) :: outfile,name,unit,varname,varlname
      integer :: year,month,day,domain
#include "netcdf.inc"
      integer ncid,iret,i,j,k,dims(4),start(4),count(4),m
      real (kind=r_knd) :: spval=-9.9e12, tt
      real (kind=r_knd) :: dx,dy,rhox,rhoy,ux,vy,uy,vx,rhodz
      real, parameter :: pi       = 
     &        3.14159265358979323846264338327950588
      real, parameter :: radius   = 6370.0e3
      real :: degtm   = radius/180*pi
      real (kind=r_knd) :: buf1(imt,jmt),buf2(imt,jmt),buf3(imt,jmt)
      real (kind=r_knd), allocatable :: mean_wt(:,:,:),mean_wu(:,:,:)
      real (kind=r_knd) :: dxt,dxu,dxum,dyt,dyu,dyum,cst,csu,cstp
      real (kind=r_knd) :: due,duw,dun,dus,dzt,dx1,dx2,dz

      integer :: lon_tdim,lon_tid
      integer :: lon_udim,lon_uid
      integer :: lat_tdim,lat_tid
      integer :: lat_udim,lat_uid
      integer :: depth_tdim,depth_tid
      integer :: depth_wdim,depth_wid
      integer :: itimedim,itimeid
      integer :: ekeid,mkeid,sshid,msshid,n,len
      integer :: t2id,t4id,T2xid,T2yid,t1id,t3id
      integer :: ubid(nt+1),vbid(nt+1),wbid(nt+1)
      integer :: ub2id(nt+1),vb2id(nt+1),wb2id(nt+1)
      integer :: ubbid(nt+1),vbbid(nt+1),wbbid(nt+1)
      integer :: tmid(nt+1),tvid(nt+1)
      integer :: 
     &days_in_month(12) = (/31,28,31,30,31,30,31,31,30,31,30,31/)

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

      if (enable_monthly_variances) then
        tt=tt-days_in_month(month)/2.
      elseif (enable_seasonal_variances) then
        tt=tt-(days_in_month(month)    ! month should be always ge 3
     &        +days_in_month(month-1)
     &        +days_in_month(month-2))/2.
      endif

      print*,' writing to file ',outfile(1:len_trim(outfile))

      ncid = nccre (outfile, NCCLOB, iret)
      iret=nf_set_fill(ncid, NF_NOFILL, iret)
c     dimensions
      lon_tdim  = ncddef(ncid, 'Longitude_t', imt, iret)
      Lon_udim  = ncddef(ncid, 'Longitude_u', imt, iret)
      Lat_tdim  = ncddef(ncid, 'Latitude_t',  jmt, iret)
      Lat_udim  = ncddef(ncid, 'Latitude_u',  jmt, iret)
      depth_wdim = ncddef(ncid, 'depth_w',  km, iret)
      depth_tdim = ncddef(ncid, 'depth_t',  km, iret)
      iTimedim  = ncddef(ncid, 'Time', 1, iret)
c     grid variables
      dims(1)  = Lon_tdim
      Lon_tid  = ncvdef (ncid,'Longitude_t',NCFLOAT,1,dims,iret)
      dims(1)  = Lon_udim
      Lon_uid  = ncvdef (ncid,'Longitude_u',NCFLOAT,1,dims,iret)
      dims(1)  = Lat_tdim
      Lat_tid  = ncvdef (ncid,'Latitude_t', NCFLOAT,1,dims,iret)
      dims(1)  = Lat_udim
      Lat_uid  = ncvdef (ncid,'Latitude_u', NCFLOAT,1,dims,iret)
      dims(1)  = iTimedim
      iTimeid   = ncvdef(ncid,'Time',       NCFLOAT,1,dims,iret)
      dims(1)  = depth_wdim
      depth_wid = ncvdef (ncid,'depth_w', NCFLOAT,1,dims,iret)
      dims(1)  = depth_tdim
      depth_tid = ncvdef (ncid,'depth_t', NCFLOAT,1,dims,iret)
c     attributes of the grid
      name = 'Longitude on T grid     '; unit = 'degrees_W       '
      call ncaptc(ncid, Lon_tid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lon_tid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Longitude on U grid     '; unit = 'degrees_W       '
      call ncaptc(ncid, Lon_uid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lon_uid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Latitude on T grid      '; unit = 'degrees_N       '
      call ncaptc(ncid, Lat_tid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lat_tid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Latitude on U grid      '; unit = 'degrees_N       '
      call ncaptc(ncid, Lat_uid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, Lat_uid, 'units',     NCCHAR, 16, unit, iret) 
      name = 'Depth of T Grid points  '; unit = 'm '
      call ncaptc(ncid, depth_tid, 'long_name',NCCHAR, 24, name, iret) 
      call ncaptc(ncid, depth_tid, 'units',    NCCHAR, 16, unit, iret) 
      name = 'Depth of W Grid points  '; unit = 'm '
      call ncaptc(ncid, depth_wid, 'long_name', NCCHAR,24, name, iret) 
      call ncaptc(ncid, depth_wid, 'units',     NCCHAR,16, unit, iret) 
      name = 'Time                    '; unit = 'days            '
      call ncaptc(ncid, iTimeid, 'long_name', NCCHAR, 24, name, iret) 
      call ncaptc(ncid, iTimeid, 'units',     NCCHAR, 16, unit, iret) 
      call ncaptc(ncid, iTimeid,'time_origin',NCCHAR, 20,
     &'31-DEC-1899 00:00:00', iret)
c     variables
      dims=(/lon_udim,lat_udim,depth_tdim,itimedim/)
      ekeid=ncvdef (ncid,'eke', NCFLOAT,4,dims,iret)
      mkeid=ncvdef (ncid,'mke', NCFLOAT,4,dims,iret)
      name = 'Eddy kinetic energy'; unit = 'm^/s^2'
      call dvcdf_here(ncid,ekeid,name,len_trim(name),unit,16,spval)
      name = 'Mean kinetic energy'; unit = 'm^2/s^2'
      call dvcdf_here(ncid,mkeid,name,len_trim(name),unit,16,spval)

      if (enable_freesurf) then
       dims=(/lon_tdim,lat_tdim,itimedim,1/)
       sshid=ncvdef (ncid,'ssh_rms', NCFLOAT,3,dims,iret)
       msshid=ncvdef (ncid,'ssh_mean', NCFLOAT,3,dims,iret)
c      attributes of variables
       name = 'Rms of Sea Surface height'; unit = 'm'
       call dvcdf_here(ncid,sshid,name,len_Trim(name),unit,16,spval)
       name = 'Mean Sea Surface height'; unit = 'm'
       call dvcdf_here(ncid,msshid,name,len_Trim(name),unit,16,spval)
      endif

      if (energy) then
       dims=(/lon_tdim,lat_tdim,depth_tdim,itimedim/)
       t2id=ncvdef (ncid,'T2', NCFLOAT,4,dims,iret)
       name = 'Baroclinic conversion rate'; unit = '(kg/m^3)^2/s'
       call dvcdf_here(ncid,T2id,name,len_Trim(name),unit,16,spval)

       dims=(/lon_udim,lat_udim,depth_tdim,itimedim/)
       T4id=ncvdef (ncid,'T4', NCFLOAT,4,dims,iret)
       name = 'Barotropic conversion rate'; unit = '(m/s)^2/s'
       call dvcdf_here(ncid,T4id,name,len_Trim(name),unit,16,spval)

       dims=(/lon_tdim,lat_tdim,depth_tdim,itimedim/)
       t2xid=ncvdef (ncid,'T2x', NCFLOAT,4,dims,iret)
       name = 'Zonal baroclinic conversion rate'; unit = '(kg/m^3)^2/s'
       call dvcdf_here(ncid,T2xid,name,len_Trim(name),unit,16,spval)
      endif

      if (energy_ext) then
       dims=(/lon_tdim,lat_tdim,depth_wdim,itimedim/)
       T3id=ncvdef (ncid,'T3', NCFLOAT,4,dims,iret)
       name = 'T3'; unit = ' '
       call dvcdf_here(ncid,T3id,name,len_Trim(name),unit,16,spval)
       dims=(/lon_tdim,lat_tdim,depth_wdim,itimedim/)
       T1id=ncvdef (ncid,'T1', NCFLOAT,4,dims,iret)
       name = 'T1'; unit = ' '
       call dvcdf_here(ncid,T1id,name,len_Trim(name),unit,16,spval)
      endif

      if (var_fluxes) then
       do n=1,nt+1
        if (n==1) then
         varname='T'        
         varlname = 'temperature'
         unit = 'mK/s'
        elseif (n==2) then
         varname='S'        
         varlname = 'salinity'
         unit = 'm/s'
        elseif (n==nt+1) then
         varname='sig0'        
         varlname = 'pot. density'
         unit = 'kg/m^2/s'
        else
         write(varname,'("tr_",i2)') n
         call replace_space_zero(varname)
         write(varlname,'("tracer ",i2)') n
         unit = 'm/s '
        endif

        dims=(/lon_tdim,lat_tdim,depth_tdim,itimedim/)
        tmid(n)=ncvdef (ncid,'mean_'//varname(1:len_trim(varname))
     &                 , NCFLOAT,4,dims,iret)
        name = 'Mean '//varlname(1:len_trim(varlname)); 
        call dvcdf_here(ncid,tmid(n),name,len_Trim(name),unit,16,spval)

        dims=(/lon_tdim,lat_tdim,depth_tdim,itimedim/)
        tvid(n)=ncvdef (ncid,'var_'//varname(1:len_trim(varname))
     &                 , NCFLOAT,4,dims,iret)
        name = 'Variance of '//varlname(1:len_trim(varlname)); 
        call dvcdf_here(ncid,tvid(n),name,len_Trim(name),unit,16,spval)

        dims=(/lon_udim,lat_tdim,depth_tdim,itimedim/)
        ubid(n)=ncvdef (ncid,'u'//varname(1:len_trim(varname))
     &                 , NCFLOAT,4,dims,iret)
        name = 'Zonal eddy flux of '//varlname(1:len_trim(varlname)); 
        call dvcdf_here(ncid,ubid(n),name,len_Trim(name),unit,16,spval)

        ubbid(n)=ncvdef (ncid,'u'//varname(1:len_trim(varname))//'2'
     &                 , NCFLOAT,4,dims,iret)
        name = 'Zonal variance flux of '//varlname(1:len_trim(varlname)); 
        call dvcdf_here(ncid,ubbid(n),name,len_Trim(name),unit,16,spval)

        ub2id(n)=ncvdef (ncid,'u2'//varname(1:len_trim(varname))
     &                 , NCFLOAT,4,dims,iret)
        name = 'Zonal eddy flux 2 of '//varlname(1:len_trim(varlname)); 
        call dvcdf_here(ncid,ub2id(n),name,len_Trim(name),unit,16,spval)

        dims=(/lon_tdim,lat_udim,depth_tdim,itimedim/)
        vbid(n)=ncvdef (ncid,'v'//varname(1:len_trim(varname))
     &                 , NCFLOAT,4,dims,iret)
        name='Meridional eddy flux of '//varlname(1:len_trim(varlname))
        call dvcdf_here(ncid,vbid(n),name,len_Trim(name),unit,16,spval)

        vbbid(n)=ncvdef (ncid,'v'//varname(1:len_trim(varname))//'2'
     &                 , NCFLOAT,4,dims,iret)
        name='Meridional variance flux of '//
     &                     varlname(1:len_trim(varlname))
        call dvcdf_here(ncid,vbbid(n),name,len_Trim(name),unit,16,spval)

        vb2id(n)=ncvdef (ncid,'v2'//varname(1:len_trim(varname))
     &                 , NCFLOAT,4,dims,iret)
        name='Meridional eddy flux 2 of '//
     &           varlname(1:len_trim(varlname))
        call dvcdf_here(ncid,vb2id(n),name,len_Trim(name),unit,16,spval)

        dims=(/lon_tdim,lat_tdim,depth_wdim,itimedim/)
        wbid(n)=ncvdef (ncid,'w'//varname(1:len_trim(varname))
     &                 , NCFLOAT,4,dims,iret)
        name='Vertical eddy flux of '//varlname(1:len_trim(varlname))
        call dvcdf_here(ncid,wbid(n),name,len_Trim(name),unit,16,spval)

        wbbid(n)=ncvdef (ncid,'w'//varname(1:len_trim(varname))//'2'
     &                 , NCFLOAT,4,dims,iret)
        name='Vertical variance flux of '//
     &                    varlname(1:len_trim(varlname))
        call dvcdf_here(ncid,wbbid(n),name,len_Trim(name),unit,16,spval)

        wb2id(n)=ncvdef (ncid,'w2'//varname(1:len_trim(varname))
     &                 , NCFLOAT,4,dims,iret)
        name='Vertical eddy flux 2 of '//varlname(1:len_trim(varlname))
        call dvcdf_here(ncid,wb2id(n),name,len_Trim(name),unit,16,spval)

       enddo
      endif

      iret= nf_put_att_int(ncid,nf_global,
     &   'Number_of_averaged_time_steps ',nf_int,1,counter)
        if (iret.ne.0) print*,nf_strerror(iret)

      write( name,'("SFLAME sub domain # ",i3)') domain
      if (enable_monthly_variances) 
     &   name=name(1:len_trim(name))//' monthly means'
      if (enable_seasonal_variances) 
     &   name=name(1:len_trim(name))//' seasonal means'
      len=len_trim(name)
      iret= nf_put_att_text(ncid,nf_global,'history',len,name)

      call ncendf(ncid, iret)

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

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

c     convert to SI units
      mean_u=mean_u/100.   
      var_u=var_u/100./100. 
      mean_ssh=mean_ssh/100.
      var_ssh=var_ssh/100./100.

      if (energy) then
       urho=urho/100.0 * 1000.0
       uv = uv/100.0/100.0
       rhom = rhom *1000.0
      endif

      if (energy_ext) then
        wrho=wrho/100.0*1000.0
        pressm=pressm /100.0 /100.0 *1000.0
        upress=upress/100.0 /100.0/100.0 *1000.0
      endif

      if (var_fluxes) then
         mean_adv=mean_adv/100;  
         mean_eddy_flux=mean_eddy_flux/100; 
         mean_var_flux = mean_var_flux/100;
      endif
c
c    get means and variances:
c     sum(T*) = sum(T^2) - N T_bar^2
c     N T_bar = sum(T_n)
c     1/N [sum(T_n)]^2 = N T_bar^2
c
      do n=1,2
       do j=1,jmt
        do k=1,km
         do i=1,imt
          var_u(i,k,j,n)= var_u(i,k,j,n)-mean_u(i,k,j,n)**2/counter 
          var_u(i,k,j,n)= var_u(i,k,j,n)/counter
          mean_u(i,k,j,n)= mean_u(i,k,j,n)/counter
         enddo
        enddo
       enddo
      enddo

      if (energy) then
       do j=1,jmt
        do k=1,km
         do i=1,imt
           rhom(i,k,j)=rhom(i,k,j)/counter
           urho(i,k,j,1)= urho(i,k,j,1) /counter
     &                  - rhom(i,k,j)*mean_u(i,k,j,1)
           urho(i,k,j,2)= urho(i,k,j,2) /counter
     &                  - rhom(i,k,j)*mean_u(i,k,j,2)
           uv(i,k,j)= uv(i,k,j) /counter
     &                  - mean_u(i,k,j,1)*mean_u(i,k,j,2)
         enddo
        enddo
       enddo
      endif

      if (energy_ext) then
c   
c    get T1 and T3
c    first calculate mean vertical velocity
       allocate(mean_wu(imt,km,jmt) ,mean_wt(imt,km,jmt) )
       mean_wu=0.; mean_wt=0.
c
c    vertical velocity on T grid
       do k=1,km
        buf1=0.0;buf2=0.0
        do j=2,jmt-1
          do i=2,imt-1
            dxt = (xu(i)-xu(i-1))*degtm
            dxu = (xt(i+1)-xt(i))*degtm
            dxum = (xt(i)-xt(i-1))*degtm
            dyt = (yu(j)-yu(j-1))*degtm
            dyu = (yt(j+1)-yt(j))*degtm
            dyum = (yt(j)-yt(j-1))*degtm
            csu  = cos( yu(j)/180*pi )
            buf1(i,j) = 0.5*(mean_u(i,k,j,1)*dyu +    ! adv_vet
     &                        mean_u(i,k,j-1,1)*dyum)/dyt
            buf2(i,j) = 0.5*(mean_u(i,k,j,2)*dxu +    ! adv_vnt
     &                        mean_u(i-1,k,j,2)*dxum)*csu/dxt
          enddo
        enddo
        do j=2,jmt
          do i=2,imt
            dxt = (xu(i)-xu(i-1))*degtm
            dyt = (yu(j)-yu(j-1))*degtm
            dzt  = (zw(k)- zw(k-1))
            cst  = cos(yt(j)/180*pi)
            mean_wt(i,k,j) =
     &                   ((buf1(i,j) - buf1(i-1,j))/dxt
     &                   +(buf2(i,j) - buf2(i,j-1))/dyt
     &                   )/cst*dzt
          enddo
        enddo
       end do

       do j=1,jmt
        do k=2,km ! mean_wt(:,0,:)=0.0
          do i=1,imt
            mean_wt(i,k,j) = mean_wt(i,k,j) + mean_wt(i,k-1,j)
          enddo
        enddo
       end do

c    and vertical velocity on the u grid
       do j=2,jmt-1
        do k=1,km
          do i=2,imt-1
            duw = (xu(i)   - xt(i))*degtm
	    due = (xt(i+1) - xu(i))*degtm
	    dun = (yt(j+1) - yu(j))*degtm
            dus = (yu(j)   - yt(j))*degtm
            dxu = (xt(i+1)-xt(i))*degtm
            dyu = (yt(j+1)-yt(j))*degtm
            csu = cos( yu(j)/180*pi )
            cst  = cos( yt(j)/180*pi )
            cstp = cos( yt(j+1)/180*pi )
	    mean_wu(i,k,j) = 1./dyu/csu/dxu*(
     &                    mean_wt(i,k,j)*duw*dus*cst
     &                  + mean_wt(i+1,k,j)*due*dus*cst
     &                  + mean_wt(i,k,j+1)*duw*dun*cstp
     &                  + mean_wt(i+1,k,j+1)*due*dun*cstp )
          enddo
        enddo
       enddo

c    now get means and variances of averaged quantities
c    store T3 in wrho and T1 in pressm

       do j=1,jmt
        do k=1,km
         do i=1,imt
           pressm(i,k,j)=pressm(i,k,j)/counter
           wrho(i,k,j)= wrho(i,k,j) /counter
     &                  - rhom(i,k,j)*mean_wt(i,k,j)
           upress(i,k,j,1)= upress(i,k,j,1) /counter
     &                  - pressm(i,k,j)*mean_u(i,k,j,1)
           upress(i,k,j,2)= upress(i,k,j,2) /counter
     &                  - pressm(i,k,j)*mean_u(i,k,j,2)
           upress(i,k,j,3)= upress(i,k,j,3) /counter
     &                  - pressm(i,k,j)*mean_wu(i,k,j)
           pressm(i,k,j)=rhom(i,k,j)*mean_wt(i,k,j)
         enddo
        enddo
       enddo

       deallocate( mean_wt, mean_wu)
      endif ! energy_ext


      if (energy) then
c
c    get T2 = u_prime rho_prime (rho_bar)_x and v_prime rho_prime (rho_bar)_y 
c    store T2 in rhom
c
       do j=1,jmt-1
        do k=1,km
         do i=1,imt-1
          dx=(xt(i+1)-xt(i))*cos(yt(j)/180*pi)*degtm
          rhox = (rhom(i+1,k,j) - rhom(i,k,j))/dx
          urho(i,k,j,1)= urho(i,k,j,1) *rhox
          dy=(yt(j+1)-yt(j))*degtm
          rhoy = (rhom(i,k,j+1) - rhom(i,k,j))/dy
          urho(i,k,j,2)= urho(i,k,j,2) *rhoy
         enddo
        enddo
       enddo
       rhom=0.0 ! here we are loosing rho_bar, it becomes T2
       do k=1,km
        buf1=urho(:,k,:,1)
        buf2=urho(:,k,:,2)
        urho(:,k,:,1:2)=0.0
        do j=2,jmt-1
         do i=2,imt-1
           urho(i,k,j,1)=(buf1(i,j)+buf1(i-1,j) )/2.0
           urho(i,k,j,2)=(buf2(i,j)+buf2(i-1,j) )/2.0
           rhom(i,k,j)= urho(i,k,j,1)+urho(i-1,k,j,2) 
         enddo
        enddo
       enddo
c
c    get T4 = 
c    u_prime^2 (u_bar)_x + v_prime^2 (v_bar)_y + 
c    u_prime v_prime ( (u_bar)_y + (v_bar)x )
c     store T4 in uv and first part of T4 in urho(:,:,:,1)
c
       do k=1,km
        buf1=0.0;buf2=0.0;
        buf3=uv(:,k,:);uv(:,k,:)=0.0
        do j=2,jmt-1
         do i=2,imt-1
          dx=(xu(i+1)-xu(i))*cos(yu(j)/180*pi)*degtm
          ux = (mean_u(i+1,k,j,1) - mean_u(i,k,j,1))/dx
          dx=(xu(i)-xu(i-1))*cos(yu(j)/180*pi)*degtm
          ux = (ux + (mean_u(i,k,j,1) - mean_u(i-1,k,j,1))/dx )/2.0

          buf1(i,j)= var_u(i,k,j,1)*ux

          dy=(yu(j+1)-yu(j))*degtm
          vy = (mean_u(i,k,j+1,2) - mean_u(i,k,j,2))/dy
          dy=(yu(j)-yu(j-1))*degtm
          vy = (vy + (mean_u(i,k,j,2) - mean_u(i,k,j-1,2))/dy )/2.0

          buf2(i,j)= var_u(i,k,j,2)*vy

          dx=(xu(i+1)-xu(i))*cos(yu(j)/180*pi)*degtm
          vx = (mean_u(i+1,k,j,2) - mean_u(i,k,j,2))/dx
          dx=(xu(i)-xu(i-1))*cos(yu(j)/180*pi)*degtm
          vx = (vx + (mean_u(i,k,j,2) - mean_u(i-1,k,j,2))/dx )/2.0

          dy=(yu(j+1)-yu(j))*degtm
          uy = (mean_u(i,k,j+1,1) - mean_u(i,k,j,1))/dy
          dy=(yu(j)-yu(j-1))*degtm
          uy = (uy + (mean_u(i,k,j,1) - mean_u(i,k,j-1,1))/dy )/2.0

          uv(i,k,j)=  buf1(i,j) + buf2(i,j) +buf3(i,j)*(vx+uy)

         enddo
        enddo
       enddo

      endif ! energy

      if (var_fluxes) then
c
c    get the mean quantities and fluxes
c
c sum(u_prime b_prime)/N = sum(ub)/N - u_bar b_bar
c  u_bar sum(b_prime^2)/N + sum( u_prime b_prime^2)/N = 
c       sum(ub^2)/N - u_bar b_bar^2 - 2 b_bar sum(u_prime b_prime)/N
        mean_adv=mean_adv/counter
        mean_tracer=mean_tracer/counter
c        var_tracer=var_tracer/counter
c     sum(T*) = sum(T^2) - N T_bar^2
c     N T_bar = sum(T_n)
c     1/N [sum(T_n)]^2 = N T_bar^2
        do n=1,nt+1
        do m=1,3
        do j=1,jmt
         do k=1,km
          do i=1,imt
          var_tracer(i,k,j,m,n)= (var_tracer(i,k,j,m,n)-
     &               mean_tracer(i,k,j,m,n)**2*counter )/counter
           mean_eddy_flux(i,k,j,m,n) = mean_eddy_flux(i,k,j,m,n)/counter 
     &           - mean_adv(i,k,j,m)*mean_tracer(i,k,j,m,n)
           mean_var_flux(i,k,j,m,n) = mean_var_flux(i,k,j,m,n)/counter 
     &           - mean_adv(i,k,j,m)*mean_tracer(i,k,j,m,n)**2
     &           - 2*mean_tracer(i,k,j,m,n)*mean_eddy_flux(i,k,j,m,n)
          enddo
         enddo
        enddo
        enddo
        enddo
c
c     now compute the gauge vector potential for the eddy fluxes
c
        do n=1,nt+1

         do j=1,jmt-1
          do k=1,km-1
           do i=2,imt
c
c     theta_1 = (G_2 b_z - G_3 b_y) /( b_x^2+b_y^2+b_z^2)
c     with G= bar( u phi_prime )
c      is on grid:  T grid in x,  U grid in y, W Grid in z
c
            dx1=(xu(i)-xu(i-1))*cos(yt(j)/180*pi)*degtm
            dx2=(xu(i)-xu(i-1))*cos(yt(j+1)/180*pi)*degtm
            dy=(yt(j+1)-yt(j))*degtm
            dz  = zt(k+1)- zt(k)
            rhodz= 0.25*(
     &        ((mean_tracer(i,k,j,1,n)-mean_tracer(i-1,k,j,1,n))/dx1)**2
     &   +((mean_tracer(i,k,j+1,1,n)-mean_tracer(i-1,k,j+1,1,n))/dx2)**2
     &  +((mean_tracer(i,k+1,j,1,n)-mean_tracer(i-1,k+1,j,1,n))/dx1)**2+
     &((mean_tracer(i,k+1,j+1,1,n)-mean_tracer(i-1,k+1,j+1,1,n))/dx2)**2
     &            )
     &       +((mean_tracer(i,k,j,2,n)-mean_tracer(i,k+1,j,2,n))/dz)**2
     &       +((mean_tracer(i,k,j+1,3,n)-mean_tracer(i,k,j,3,n))/dy)**2
            if (rhodz==0.0) rhodz=1e-10
            theta(i,k,j,1,n)= (
     &       0.5*(mean_var_flux(i,k,j,2,n)+mean_var_flux(i,k+1,j,2,n))
     &       *(mean_tracer(i,k,j,2,n)-mean_tracer(i,k+1,j,2,n))/dz
     &      -0.5*(mean_var_flux(i,k,j,3,n)+mean_var_flux(i,k,j+1,3,n))
     &       *(mean_tracer(i,k,j+1,3,n)-mean_tracer(i,k,j,3,n))/dy )
     &        /rhodz
           enddo
          enddo
         enddo

         do j=2,jmt
          do k=1,km-1
           do i=1,imt-1
c
c     theta_2 = (G_3 b_x - G_1 b_z) /( b_x^2+b_y^2+b_z^2)
c      is on grid:  U grid in x,  T grid in y, W Grid in z
c
            dx=(xt(i+1)-xt(i))*cos(yt(j)/180*pi)*degtm
            dz=zt(k+1)-zt(k)
            dy = (yu(j)-yu(j-1))*degtm
            rhodz=
     &        ((mean_tracer(i+1,k,j,3,n)-mean_tracer(i,k,j,1,n))/dx)**2
     &       +((mean_tracer(i,k,j,1,n)-mean_tracer(i,k+1,j,1,n))/dz)**2
     &       +0.25*(
     &       +((mean_tracer(i,k,j,2,n)-mean_tracer(i,k,j-1,2,n))/dy)**2
     &   +((mean_tracer(i+1,k,j,2,n)-mean_tracer(i+1,k,j-1,2,n))/dy)**2
     & +((mean_tracer(i,k+1,j,2,n)-mean_tracer(i,k+1,j-1,2,n))/dy)**2+
     &((mean_tracer(i+1,k+1,j,2,n)-mean_tracer(i+1,k+1,j-1,2,n))/dy)**2
     &            )
            if (rhodz==0.0) rhodz=1e-10
            theta(i,k,j,2,n)= (
     &       0.5*(mean_var_flux(i,k,j,3,n)+mean_var_flux(i+1,k,j,3,n))
     &       *(mean_tracer(i+1,k,j,3,n)-mean_tracer(i,k,j,3,n))/dx
     &      -0.5*(mean_var_flux(i,k,j,1,n)+mean_var_flux(i,k+1,j,1,n))
     &       *(mean_tracer(i,k,j,1,n)-mean_tracer(i,k+1,j,1,n))/dz )
     &        /rhodz
           enddo
          enddo
         enddo

         do j=1,jmt-1
          do k=2,km
           do i=1,imt-1
c
c     theta_3 = (G_1 b_y - G_2 b_x) /( b_x^2+b_y^2+b_z^2)
c      is on grid:  U grid in x,  U grid in y, T Grid in z
c
            dy = (yt(j+1)-yt(j))*degtm
            dx = (xt(i+1)-xt(i))*cos(yu(j)/180*pi)*degtm
            dz = zw(k)-zw(k-1)
            rhodz=
     &        ((mean_tracer(i+1,k,j,2,n)-mean_tracer(i,k,j,2,n))/dx)**2
     &       +((mean_tracer(i,k,j+1,1,n)-mean_tracer(i,k,j,1,n))/dy)**2
     &       +0.25*(
     &       +((mean_tracer(i,k-1,j,3,n)-mean_tracer(i,k,j,3,n))/dz)**2
     &   +((mean_tracer(i,k-1,j+1,3,n)-mean_tracer(i,k,j+1,3,n))/dz)**2
     &   +((mean_tracer(i+1,k-1,j,3,n)-mean_tracer(i+1,k,j,3,n))/dz)**2
     &+((mean_tracer(i+1,k-1,j+1,3,n)-mean_tracer(i+1,k,j+1,3,n))/dz)**2
     &            )
            if (rhodz==0.0) rhodz=1e-10
            theta(i,k,j,3,n)= (
     &       0.5*(mean_var_flux(i,k,j,1,n)+mean_var_flux(i,k,j+1,1,n))
     &       *(mean_tracer(i,k,j+1,1,n)-mean_tracer(i,k,j,1,n))/dy
     &      -0.5*(mean_var_flux(i+1,k,j,2,n)+mean_var_flux(i,k,j,2,n))
     &       *(mean_tracer(i+1,k,j,2,n)-mean_tracer(i,k,j,2,n))/dx )
     &        /rhodz
           enddo
          enddo
         enddo

c  now remove the vector potential flux from the eddy fluxes

         do j=2,jmt
          do k=2,km
           do i=2,imt
c        F_1 = u_prime b_prime - (theta_2)_z + (theta_3)_y
            dy = (yu(j)-yu(j-1))*degtm
            dz = zw(k)-zw(k-1)
            eddy_flux(i,k,j,1,n)=mean_eddy_flux(i,k,j,1,n)
     &         -(theta(i,k-1,j,2,n)-theta(i,k,j,2,n) )/dz
     &         +(theta(i,k,j,3,n)-theta(i,k,j-1,3,n) )/dy
c        F_2 = v_prime b_prime - (theta_1)_z + (theta_3)_x
            dx = (xu(i)-xu(i-1))*cos(yu(j)/180*pi)*degtm
            dz = zw(k)-zw(k-1)
            eddy_flux(i,k,j,2,n)=mean_eddy_flux(i,k,j,2,n)
     &         -(theta(i,k-1,j,1,n)-theta(i,k,j,1,n) )/dz
     &         +(theta(i,k,j,3,n)-theta(i-1,k,j,3,n) )/dx
c        F_3 = w_prime b_prime - (theta_2)_x + (theta_1)_y
            dy = (yu(j)-yu(j-1))*degtm
            dx = (xu(i)-xu(i-1))*cos(yt(j)/180*pi)*degtm
            eddy_flux(i,k,j,3,n)=mean_eddy_flux(i,k,j,3,n)
     &         -(theta(i,k,j,2,n)-theta(i-1,k,j,2,n) )/dx
     &         +(theta(i,k,j,1,n)-theta(i,k,j-1,1,n) )/dy
           enddo
          enddo
         enddo


        enddo


      endif ! var_fluxes

c
c   store energy terms  in var_u and mean_u
c
      do j=1,jmt
       do k=1,km
        do i=1,imt
         var_u(i,k,j,1)=(var_u(i,k,j,1)+var_u(i,k,j,2))/2.
         mean_u(i,k,j,1)=(mean_u(i,k,j,1)**2+mean_u(i,k,j,2)**2)/2.
        enddo
       enddo
      enddo

c
c   care about external mode
c
      do j=1,jmt
       do i=1,imt
        var_ssh(i,j)= var_ssh(i,j)-mean_ssh(i,j)**2/counter 
        var_ssh(i,j)= sqrt(var_ssh(i,j)/counter)
        mean_ssh(i,j)= mean_ssh(i,j)/counter
       enddo
      enddo

c     apply land mask and write to file
      do k=1,km
         where( kmu <k ) mean_u(:,k,:,1) = spval
         where( kmu <k ) var_u(:,k,:,1) = spval
         start=(/1,1,k,1/); count=(/imt,jmt,1,1/)
#ifdef real_8
         iret=nf_put_vara_double(ncid,mkeid,start,count,mean_u(:,k,:,1))
         iret=nf_put_vara_double (ncid,ekeid,start,count,var_u(:,k,:,1))
#else
         iret=nf_put_vara_real(ncid,mkeid,start,count,mean_u(:,k,:,1))
         iret=nf_put_vara_real(ncid,ekeid,start,count,var_u(:,k,:,1))
#endif
      enddo

      if (energy) then
       do k=1,km
         where( kmu <k ) uv(:,k,:) = spval
         where( kmt <k ) rhom(:,k,:) = spval
         where( kmt <k ) urho(:,k,:,1) = spval
         start=(/1,1,k,1/); count=(/imt,jmt,1,1/)
#ifdef real_8
         iret=nf_put_vara_double(ncid,T4id,start,count,uv(:,k,:))
         iret=nf_put_vara_double(ncid,T2id,start,count,rhom(:,k,:))
         iret=nf_put_vara_double(ncid,T2xid,start,count,urho(:,k,:,1))
#else
         iret=nf_put_vara_real(ncid,T4id,start,count,uv(:,k,:))
         iret=nf_put_vara_real(ncid,T2id,start,count,rhom(:,k,:))
         iret=nf_put_vara_real(ncid,T2xid,start,count,urho(:,k,:,1))
#endif
       enddo
      endif

      if (energy_ext) then
       do k=1,km
         where( kmt <k ) wrho(:,k,:) = spval
         where( kmt <k ) pressm(:,k,:) = spval
         start=(/1,1,k,1/); count=(/imt,jmt,1,1/)
#ifdef real_8
         iret=nf_put_vara_double(ncid,T1id,start,count,pressm(:,k,:))
         iret=nf_put_vara_double(ncid,T3id,start,count,wrho(:,k,:))
#else
         iret=nf_put_vara_real(ncid,T1id,start,count,pressm(:,k,:))
         iret=nf_put_vara_real(ncid,T3id,start,count,wrho(:,k,:))
#endif
       enddo
      endif

      if (var_fluxes) then
       do k=1,km
         do n=1,nt+1
          do m=1,3
           where( kmt <k ) mean_adv(:,k,:,m) = spval
           where( kmt <k ) mean_tracer(:,k,:,m,n) = spval
           where( kmt <k ) var_tracer(:,k,:,m,n) = spval
           where( kmt <k ) mean_eddy_flux(:,k,:,m,n) = spval
           where( kmt <k ) eddy_flux(:,k,:,m,n) = spval
           where( kmt <k ) mean_var_flux(:,k,:,m,n) = spval
          enddo
         enddo
         start=(/1,1,k,1/); count=(/imt,jmt,1,1/)
#ifdef real_8
         do n=1,nt+1
          iret=nf_put_vara_double(ncid,tmid(n),start,count,
     &           mean_tracer(:,k,:,1,n))
          iret=nf_put_vara_double(ncid,tvid(n),start,count,
     &           var_tracer(:,k,:,1,n))
          iret=nf_put_vara_double(ncid,ubid(n),start,count,
     &           mean_eddy_flux(:,k,:,1,n))
          iret=nf_put_vara_double(ncid,vbid(n),start,count,
     &           mean_eddy_flux(:,k,:,2,n))
          iret=nf_put_vara_double(ncid,wbid(n),start,count,
     &           mean_eddy_flux(:,k,:,3,n))
          iret=nf_put_vara_double(ncid,ubbid(n),start,count,
     &           mean_var_flux(:,k,:,1,n))
          iret=nf_put_vara_double(ncid,vbbid(n),start,count,
     &           mean_var_flux(:,k,:,2,n))
          iret=nf_put_vara_double(ncid,wbbid(n),start,count,
     &           mean_var_flux(:,k,:,3,n))
          iret=nf_put_vara_double(ncid,ub2id(n),start,count,
     &           eddy_flux(:,k,:,1,n))
          iret=nf_put_vara_double(ncid,vb2id(n),start,count,
     &           eddy_flux(:,k,:,2,n))
          iret=nf_put_vara_double(ncid,wb2id(n),start,count,
     &           eddy_flux(:,k,:,3,n))
         enddo
#else
         do n=1, nt+1
          iret=nf_put_vara_real(ncid,tmid(n),start,count,
     &           mean_tracer(:,k,:,1,n))
          iret=nf_put_vara_real(ncid,tvid(n),start,count,
     &           var_tracer(:,k,:,1,n))
          iret=nf_put_vara_real(ncid,ubid(n),start,count,
     &           mean_eddy_flux(:,k,:,1,n))
          iret=nf_put_vara_real(ncid,vbid(n),start,count,
     &           mean_eddy_flux(:,k,:,2,n))
          iret=nf_put_vara_real(ncid,wbid(n),start,count,
     &           mean_eddy_flux(:,k,:,3,n))
          iret=nf_put_vara_real(ncid,ubbid(n),start,count,
     &           mean_var_flux(:,k,:,1,n))
          iret=nf_put_vara_real(ncid,vbbid(n),start,count,
     &           mean_var_flux(:,k,:,2,n))
          iret=nf_put_vara_real(ncid,wbbid(n),start,count,
     &           mean_var_flux(:,k,:,3,n))
          iret=nf_put_vara_real(ncid,ub2id(n),start,count,
     &           eddy_flux(:,k,:,1,n))
          iret=nf_put_vara_real(ncid,vb2id(n),start,count,
     &           eddy_flux(:,k,:,2,n))
          iret=nf_put_vara_real(ncid,wb2id(n),start,count,
     &           eddy_flux(:,k,:,3,n))
         enddo
#endif
       enddo
      endif

      if (enable_freesurf) then
       where( kmt == 0) mean_ssh = spval
       where( kmt == 0) var_ssh = spval
       start=(/1,1,1,1/); count=(/imt,jmt,1,1/)
#ifdef real_8
       iret= nf_put_vara_double(ncid,sshid,start,count,var_ssh)
       iret= nf_put_vara_double(ncid,msshid,start,count,mean_ssh)
#else
       iret= nf_put_vara_real (ncid,sshid,start,count,var_ssh)
       iret= nf_put_vara_real (ncid,msshid,start,count,mean_ssh)
#endif
      endif
      call ncclos (ncid, iret)
      end subroutine write_variances

#endif
