#include "options.inc"



       module variances_module
c
c=======================================================================
c      perform averages for some second moments 
c
c      Note that:   sum((T_prime)^2) = sum(T^2) - N bar(T)^2
c      and that N is lost when converting sum(T) to NetCDF
c
c      namelist parameter to be used:
c
c          enable_monthly_variances
c          enable_seasonal_variances
c          enable_daily_variances
c
c      u u_t  +  u u u_x  = - u p_x/rho_0  = - (u p/rho_0)_x + p/rho_0 u_x 
c                = -(u p/rho_0)_x  - p/rho_0 w_z = -(u p/rho_0)_x  - (wp/rho_0)_z + w p_z/rho_0 
c      (u^2/2)_t + u u u_x  + (up/rho_0)_x+(wp/rho_0)_z = w p_z/rho_0 = - g w rho/rho_0 = w b  
c
c
c ------------------------------------------------
c what is needed for the fluxes
c  eddy fluxes :
c sum(ub) = sum( u_bar b_bar + ubar b_prime + u_prime b_bar + u_prime b_prime)
c         = N u_bar b_bar + sum(u_prime b_prime)
c sum(u_prime b_prime)/N = sum(ub)/N - u_bar b_bar
c
c  fluxes of variance:
c sum(ub^2) = u_bar sum(b^2) + sum(u_prime b^2)
c  = u_bar ( sum(b_prime^2) + N b_bar^2) + sum( u_prime (b_bar + b_prime)^2)
c  = N u_bar b_bar^2 + u_bar sum(b_prime^2) 
c  + sum( u_prime (b_bar^2 + b_prime^2 + 2b_bar b_prime ) )
c  = N u_bar b_bar^2 + u_bar sum(b_prime^2) 
c  + sum( u_prime b_prime^2) + 2 b_bar sum( u_prime b_prime  )
c
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
c    =  sum(ub^2)/N - u_bar b_bar^2 - 2 b_bar (sum(ub)/N - u_bar b_bar)
c    =  sum(ub^2)/N + u_bar b_bar^2 - 2 b_bar sum(ub)/N
c  
c  so we need u_bar b_bar sum(ub) sum(ub^2) and sum(b^2)
c  ( the last one only to separate out the triple fluxes)
c
c  we write b_bar, etc averaged on the three edges of the tracer boxes
c ------------------------------------------------
c
c=======================================================================
c
       implicit none
       real, allocatable, dimension(:,:,:,:) :: umean,uvar
       real, allocatable, dimension(:,:,:)   :: uv,T2,T4
       real, allocatable, dimension(:,:,:)   :: bmean
       real, allocatable, dimension(:,:,:)   :: var_tracer
       real, allocatable, dimension(:,:,:,:) :: mean_adv
       real, allocatable, dimension(:,:,:,:) :: mean_eddy_flux
       real, allocatable, dimension(:,:,:,:) :: mean_var_flux
       integer :: counter  = 0
       logical :: active = .false.
       end module variances_module


       subroutine init_variances
c
c      initialize module (should be done in any case)
c
       use spflame_module
       use variances_module
       implicit none
       character (len=80) :: name
       integer :: io,k,n,nr_blk,i,j,m
       real, allocatable ::  buf(:,:)
       logical :: eof , nofile

       if (enable_monthly_variances.or.
     &     enable_seasonal_variances.or.
     &     enable_daily_variances ) active = .true.

       if (.not. active) return ! return if not active

       if (my_pe==0) then
        print*,''
        print*,' setting up averages of second moments'
        if (enable_daily_variances) then
         print*,' based on daily means '
        endif
        if (enable_monthly_variances) then
         print*,' based on monthly means '
        endif
        if (enable_seasonal_variances) then
         print*,' based on seasonal means '
        endif
       endif
c
c      allocate memory and initialize fields and counter
c
       allocate( umean(is_pe-1:ie_pe+1,km,js_pe-1:je_pe+1,2) )
       allocate( uvar(is_pe:ie_pe,km,js_pe:je_pe,2) )
       umean=0.; uvar=0. 
       allocate( uv(is_pe:ie_pe,km,js_pe:je_pe) )
       allocate( T2(is_pe:ie_pe,km,js_pe:je_pe) )
       allocate( T4(is_pe:ie_pe,km,js_pe:je_pe) )
       uv=0.0;T2=0;T4=0.

       allocate(       mean_adv(is_pe:ie_pe,km,js_pe:je_pe,3) )
       allocate(    bmean(is_pe-1:ie_pe+1,km,js_pe-1:je_pe+1) )
       allocate(     var_tracer(is_pe:ie_pe,km,js_pe:je_pe) )
       allocate( mean_eddy_flux(is_pe:ie_pe,km,js_pe:je_pe,3) )
       allocate(  mean_var_flux(is_pe:ie_pe,km,js_pe:je_pe,3) )
       mean_adv=0.0; bmean=0.0; var_tracer=0.0
       mean_eddy_flux=0.0; mean_var_flux = 0.0

       counter = 0
c
c      read unfinished averages
c
       if (my_pe==0) then
        write(name,'("variances_",i3,".dta")') sub_domain
        call replace_space_zero(name)
        print*,' reading unfinished variances from file ',
     &          name(1:len_trim(name))
        call get_free_iounit(io);nofile=.false.
        open(io,file=name,form='unformatted',status='old',err=30)
        read(io) counter
        print*,' (counter=',counter,')'
        goto 40
 30     nofile=.true.
 40     continue
       endif

       call bcast_integer(counter,1,0)
       call bcast_logical(nofile,1,0)

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

       allocate( buf(imt,jmt) )

       do n=1,2
        do k=1,km
         if (my_pe==0) read(io) buf
         call pe0_send_2D(buf)
         umean(is_pe:ie_pe,k,js_pe:je_pe,n)=
     &       buf(is_pe:ie_pe,js_pe:je_pe)
        enddo
       enddo

       do n=1,2
        do k=1,km
         if (my_pe==0) read(io) buf
         call pe0_send_2D(buf)
         uvar(is_pe:ie_pe,k,js_pe:je_pe,n)=
     &        buf(is_pe:ie_pe,js_pe:je_pe)
        enddo
       enddo

       do k=1,km
           if (my_pe==0) read(io) buf ! uv
           call pe0_send_2D(buf)
           uv(is_pe:ie_pe,k,js_pe:je_pe)=
     &        buf(is_pe:ie_pe,js_pe:je_pe)
       enddo

       do k=1,km
        do n=1,3  ! sum(u)
           if (my_pe==0) read(io) buf 
           call pe0_send_2D(buf)
           mean_adv(is_pe:ie_pe,k,js_pe:je_pe,n)=
     &        buf(is_pe:ie_pe,js_pe:je_pe)
        enddo
        if (my_pe==0) read(io) buf 
        call pe0_send_2D(buf)
        bmean(is_pe:ie_pe,k,js_pe:je_pe)=buf(is_pe:ie_pe,js_pe:je_pe)

        if (my_pe==0) read(io) buf 
        call pe0_send_2D(buf)
        var_tracer(is_pe:ie_pe,k,js_pe:je_pe)=
     &        buf(is_pe:ie_pe,js_pe:je_pe)
        do m=1,3 ! sum(ub)
            if (my_pe==0) read(io) buf 
            call pe0_send_2D(buf)
            mean_eddy_flux(is_pe:ie_pe,k,js_pe:je_pe,m)=
     &        buf(is_pe:ie_pe,js_pe:je_pe)
        enddo
        do m=1,3 ! sum(ub^2)
            if (my_pe==0) read(io) buf 
            call pe0_send_2D(buf)
            mean_var_flux(is_pe:ie_pe,k,js_pe:je_pe,m)=
     &        buf(is_pe:ie_pe,js_pe:je_pe)
        enddo
       enddo ! k=1,km

       deallocate( buf )

       if (my_pe==0) then
         close(io)
         print*,' done setting up variance averaging'
       endif

       end subroutine init_variances



       subroutine write_unfinished_variances
       use spflame_module
       use variances_module
       implicit none
       character (len=80) :: name
       integer :: i,j,k,n,io,m
       real, allocatable ::  buf(:,:)
c
c      determine file name
c
       write(name,'("variances_",i3,".dta")') sub_domain
       call replace_space_zero(name)
       if (my_pe==0) then
          print*,' writing unfinished variances to file ',
     &             name(1:len_trim(name))
          call get_free_iounit(io)
          open(io,file=name,form='unformatted',status='unknown')
          write(io) counter
       endif

       allocate( buf(imt,jmt) )

       do n=1,2
        do k=1,km
         buf(is_pe:ie_pe,js_pe:je_pe)=umean(is_pe:ie_pe,k,js_pe:je_pe,n)
         call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
        enddo
       enddo

       do n=1,2
        do k=1,km
         buf(is_pe:ie_pe,js_pe:je_pe)=uvar(is_pe:ie_pe,k,js_pe:je_pe,n)
         call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
        enddo
       enddo

       do k=1,km
        buf(is_pe:ie_pe,js_pe:je_pe)=uv(is_pe:ie_pe,k,js_pe:je_pe)
        call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
       enddo

       do k=1,km
          do n=1,3  ! sum(u)
           buf(is_pe:ie_pe,js_pe:je_pe) = 
     &           mean_adv(is_pe:ie_pe,k,js_pe:je_pe,n)
           call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
          enddo
          buf(is_pe:ie_pe,js_pe:je_pe)=bmean(is_pe:ie_pe,k,js_pe:je_pe)
          call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
          buf(is_pe:ie_pe,js_pe:je_pe)=
     &        var_tracer(is_pe:ie_pe,k,js_pe:je_pe)
          call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
          do m=1,3 ! sum(ub)
            buf(is_pe:ie_pe,js_pe:je_pe)=
     &        mean_eddy_flux(is_pe:ie_pe,k,js_pe:je_pe,m)
            call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
          enddo
          do m=1,3 ! sum(ub^2)
            buf(is_pe:ie_pe,js_pe:je_pe)=
     &        mean_var_flux(is_pe:ie_pe,k,js_pe:je_pe,m)
            call pe0_recv_2D(buf); if (my_pe==0) write(io) buf
           enddo
       enddo ! k=1,km

       deallocate(buf)

       if (my_pe==0) close(io)

       end subroutine write_unfinished_variances






       subroutine variances
c
c      perform averages of second moments
c
       use spflame_module
       use variances_module
       implicit none
       integer :: i,j,k,n,m
       real :: sig(is_pe-1:ie_pe+1,km,js_pe-1:je_pe+1)
       real :: fxa

       if (.not. active) return ! return if not active

       if (.not. eulerback1 ) then ! not for first part of euler backward

        counter=counter+1 ! increase counter

        do n=1,2  ! average U and V
         do j=js_pe,je_pe
          do k=1,km
           do i=is_pe,ie_pe
            umean(i,k,j,n)=umean(i,k,j,n)+u(i,k,j,n,tau)
            uvar(i,k,j,n) =uvar(i,k,j,n) +u(i,k,j,n,tau)**2
           enddo
          enddo
         enddo
        enddo
c
c     this is not exact but simpler
c
        do j=js_pe,je_pe
          do k=1,km
           do i=is_pe,ie_pe
            uv(i,k,j)     = uv(i,k,j) + u(i,k,j,1,tau)*u(i,k,j,2,tau)  ! u*v
           enddo
          enddo
        enddo
c
c     sum(u)
c
        do j=js_pe,je_pe
         do k=1,km
           do i=is_pe,ie_pe
            mean_adv(i,k,j,1)=mean_adv(i,k,j,1)+adv_vet(i,k,j)
            mean_adv(i,k,j,2)=mean_adv(i,k,j,2)+adv_vnt(i,k,j)
            mean_adv(i,k,j,3)=mean_adv(i,k,j,3)+adv_vbt(i,k,j)
           enddo
         enddo
        enddo
c
c     calculate model potential density
c
        do j=js_pe-1,je_pe+1
          do k=1,km
           call model_dens(t(is_pe-1,k,j,1,tau), t(is_pe-1,k,j,2,tau),
     &                           sig(is_pe-1,k,j),1,ie_pe-is_pe+3
#ifdef partial_cell
     &                       ,ztp(is_pe-1,k,j)
#endif
     &                       )
          enddo
        enddo
c
c     sum(b), sum(b^2), sum(ub) and sum(ub^2)
c 
        do j=js_pe,je_pe
          do k=1,km
           do i=is_pe,ie_pe
            fxa = sig(i,k,j) 
            bmean(i,k,j)=bmean(i,k,j)+fxa
            var_tracer(i,k,j) =var_tracer(i,k,j)+fxa**2

            fxa = 0.5*(sig(i,k,j) + sig(i+1,k,j))
            mean_eddy_flux(i,k,j,1)=mean_eddy_flux(i,k,j,1)
     &                +fxa*adv_vet(i,k,j)
            mean_var_flux(i,k,j,1)=mean_var_flux(i,k,j,1)
     &                +fxa**2*adv_vet(i,k,j)

            fxa = 0.5*(sig(i,k,j) + sig(i,k,j+1))
            mean_eddy_flux(i,k,j,2)=mean_eddy_flux(i,k,j,2)
     &                +fxa*adv_vnt(i,k,j)
            mean_var_flux(i,k,j,2)=mean_var_flux(i,k,j,2)
     &                +fxa**2*adv_vnt(i,k,j)
           enddo
          enddo

          do k=1,km-1
           do i=is_pe,ie_pe
            fxa = 0.5*(sig(i,k,j) + sig(i,k+1,j))
            mean_eddy_flux(i,k,j,3)=mean_eddy_flux(i,k,j,3)
     &                +fxa*adv_vbt(i,k,j)
            mean_var_flux(i,k,j,3)=mean_var_flux(i,k,j,3)
     &                +fxa**2*adv_vbt(i,k,j)
           enddo
          enddo
        enddo
c
c       write finished averages to a binary file
c
        if ((end_of_month .and.enable_monthly_variances ).or. 
     &      (end_of_day   .and.enable_daily_variances   ).or. 
     &      (end_of_season.and.enable_seasonal_variances)) then
c
c        write netcdf file 
c
         call calculate_variances()
         call write_variances()
c
c        zero out counter and arrays
c
         counter=0
         umean=0.;uvar=0.; uv=0.
         mean_adv=0.0; bmean=0.0; var_tracer=0.0
         mean_eddy_flux=0.0; mean_var_flux = 0.0
        endif !  end_of_month ...
       endif ! eulerback1
       if (last_time_step) call write_unfinished_variances()
       end subroutine variances







       subroutine calculate_variances
       use spflame_module
       use variances_module
       implicit none
       integer :: i,j,k,n
       integer :: is,ie,js,je
       real :: fxa,fxb,fxc
       real :: fx(is_pe-1:ie_pe+1,km,js_pe-1:je_pe+1)

       is=is_pe; ie=ie_pe; js=js_pe; je=je_pe
c
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=js,je
         do k=1,km
          do i=is,ie
           uvar(i,k,j,n)= uvar(i,k,j,n)-umean(i,k,j,n)**2/counter 
           uvar(i,k,j,n)= uvar(i,k,j,n)/counter
           umean(i,k,j,n)= umean(i,k,j,n)/counter
           uv(i,k,j)= uv(i,k,j)/counter-umean(i,k,j,1)*umean(i,k,j,2)
          enddo
         enddo
        enddo
       enddo
       uvar = uvar * 1e-4; uv = uv * 1e-4; umean = umean * 1e-2
c
c  T4 = u_prime^2 (u_bar)_x + v_prime^2 (v_bar)_y + u_prime v_prime ( (u_bar)_y + (v_bar)x )
c
       call border_exchg(umean(:,:,:,1),km,1)
       call border_exchg(umean(:,:,:,2),km,1)
       do k=1,km
        do j=js,je
         do i=is,ie
          fxa=0.5*100*
     &      ((umean(i+1,k,j,1)-umean(i  ,k,j,1))*dxtr(i+1)*csur(j) 
     &      +(umean(i  ,k,j,1)-umean(i-1,k,j,1))*dxtr(i  )*csur(j) )
          fxb=0.5*100*
     &      ((umean(i,k,j+1,2)-umean(i,k,j  ,2))*dytr(j+1) 
     &      +(umean(i,k,j  ,2)-umean(i,k,j-1,2))*dytr(j  ) )
          T4(i,k,j)= uvar(i,k,j,1)*fxa +uvar(i,k,j,2)*fxb 
     &              +uv(i,k,j)*(fxa+fxb)
         enddo
        enddo
       enddo
c
c        Density stuff.  Convert density to buoyancy:   b = - g rho/rho_0
c            b has units of accelaration m/s^2 
c
       fxa = - grav/rho0  /100.
       mean_adv=mean_adv/100;  
       mean_eddy_flux=mean_eddy_flux/100 * fxa; 
       mean_var_flux = mean_var_flux/100 * fxa**2;
       bmean = bmean*fxa 

       mean_adv=mean_adv/counter
       bmean=bmean/counter
       var_tracer=(var_tracer-bmean**2*counter )/counter

       call border_exchg(bmean,km,1)
       do j=js,je
        do k=1,km
         do i=is,ie
          fxa = 0.5*( bmean(i,k,j)+bmean(i+1,k,j) )
          mean_eddy_flux(i,k,j,1) = mean_eddy_flux(i,k,j,1)/counter 
     &         -mean_adv(i,k,j,1)*fxa
          mean_var_flux(i,k,j,1) = mean_var_flux(i,k,j,1)/counter 
     &         -mean_adv(i,k,j,1)*fxa**2-2*fxa*mean_eddy_flux(i,k,j,1)

          fxa = 0.5*( bmean(i,k,j)+bmean(i,k,j+1) )
          mean_eddy_flux(i,k,j,2) = mean_eddy_flux(i,k,j,2)/counter 
     &         -mean_adv(i,k,j,2)*fxa
          mean_var_flux(i,k,j,2) = mean_var_flux(i,k,j,2)/counter 
     &         -mean_adv(i,k,j,2)*fxa**2-2*fxa*mean_eddy_flux(i,k,j,2)
         enddo
        enddo
       enddo

       do j=js,je
        do k=1,km-1
         do i=is,ie
            fxa = 0.5*( bmean(i,k,j)+bmean(i,k+1,j) )
            mean_eddy_flux(i,k,j,3) = mean_eddy_flux(i,k,j,3)/counter 
     &         -mean_adv(i,k,j,3)*fxa
            mean_var_flux(i,k,j,3) = mean_var_flux(i,k,j,3)/counter 
     &         -mean_adv(i,k,j,3)*fxa**2-2*fxa*mean_eddy_flux(i,k,j,3)
         enddo
        enddo
       enddo


c  T2 =   u'b' b_x / b_z  + v'b' b_y /b_z

       fx=0.0;T2=0.0
       do j=js,je
        do k=1,km
         do i=is,ie
          fx(i,k,j)=mean_eddy_flux(i,k,j,1)*
     &          100*(bmean(i+1,k,j)-bmean(i,k,j))*dxur(i)*cstr(j) 
         enddo
        enddo
       enddo
       call border_exchg(fx,km,1)
       do j=js,je
        do k=2,km-1
         do i=is,ie
          T2(i,k,j) = 0.5*(fx(i,k,j)+fx(i-1,k,j))
         enddo
        enddo
       enddo

       fx=0.0
       do j=js,je
        do k=1,km
         do i=is,ie
          fx(i,k,j)=mean_eddy_flux(i,k,j,2)*
     &          100*(bmean(i,k,j+1)-bmean(i,k,j))*dyur(j) 
         enddo
        enddo
       enddo
       call border_exchg(fx,km,1)
       do j=js,je
        do k=2,km-1
         do i=is,ie
          T2(i,k,j) = 0.5*(fx(i,k,j)+fx(i,k,j-1))
         enddo
        enddo
       enddo

       fx=0.0
       do j=js,je
        do k=1,km-1
         do i=is,ie
          fx(i,k,j)=100*(bmean(i,k+1,j)-bmean(i,k,j))*dzwr(k) 
         enddo
        enddo
       enddo
       do j=js,je
        do k=2,km-1
         do i=is,ie
          fxa = 0.5*(fx(i,k,j)+fx(i,k-1,j))
          if (fxa/=0.0) T2(i,k,j) = T2(i,k,j)/ fxa
         enddo
        enddo
       enddo

       end subroutine calculate_variances







       subroutine write_variances
       use spflame_module
       use variances_module
       implicit none
#ifdef netcdf_diagnostics
#include "netcdf.inc"
       character (len=80) :: fname,name,unit
       integer :: year,month,day
       integer :: ncid,iret,i,j,k,dims(4),start(4),count(4)
       real    :: spval=-9.9e12, tt
       integer :: lon_tdim,lon_tid,lon_udim,lon_uid
       integer :: lat_tdim,lat_tid,lat_udim,lat_uid
       integer :: depth_tdim,depth_tid,depth_wdim,depth_wid
       integer :: itimedim,itimeid,uid,vid,trid(nt),npe,n
       integer :: uvarid,umeanid,vvarid,vmeanid,uvid,T4id
       integer :: ubid,vbid,wbid,bmid,T2id
       integer :: ubbid,vbbid,wbbid,bbid
       integer :: is,ie,js,je
#ifdef netcdf_real4
      real (kind=4) :: var(is_pe:ie_pe,js_pe:je_pe)
      real (kind=4) , allocatable :: v2(:)
#else
      real          :: var(is_pe:ie_pe,js_pe:je_pe)
      real          , allocatable :: v2(:)
#endif


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

       year  = get_year(get_current_time())
       month = get_month_of_year(get_current_time())
       day   = get_day_of_month(get_current_time())

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

       if (my_pe == 0) then
        print*,' Writing variances to  NetCDF output file ',
     &         fname(1:len_trim(fname))
        ncid = nccre (fname, NCCLOB, iret)
        iret=nf_set_fill(ncid, NF_NOFILL, iret)
        call store_info_cdf(ncid)
c       dimensions
        lon_tdim  = ncddef(ncid, 'Longitude_t', imt, iret)
        Lon_udim  = ncddef(ncid, 'Longitude_u', imt, iret)
        Lat_tdim  = ncddef(ncid, 'Latitude_t',  jmt, iret)
        Lat_udim  = ncddef(ncid, 'Latitude_u',  jmt, iret)
        depth_wdim = ncddef(ncid, 'depth_w',  km, iret)
        depth_tdim = ncddef(ncid, 'depth_t',  km, iret)
        iTimedim  = ncddef(ncid, 'Time', 1, iret)
c       grid variables
        dims(1)  = Lon_tdim
        Lon_tid  = ncvdef (ncid,'Longitude_t',NCFLOAT,1,dims,iret)
        dims(1)  = Lon_udim
        Lon_uid  = ncvdef (ncid,'Longitude_u',NCFLOAT,1,dims,iret)
        dims(1)  = Lat_tdim
        Lat_tid  = ncvdef (ncid,'Latitude_t', NCFLOAT,1,dims,iret)
        dims(1)  = Lat_udim
        Lat_uid  = ncvdef (ncid,'Latitude_u', NCFLOAT,1,dims,iret)
        dims(1)  = iTimedim
        iTimeid   = ncvdef(ncid,'Time',       NCFLOAT,1,dims,iret)
        dims(1)  = depth_wdim
        depth_wid = ncvdef (ncid,'depth_w', NCFLOAT,1,dims,iret)
        dims(1)  = depth_tdim
        depth_tid = ncvdef (ncid,'depth_t', NCFLOAT,1,dims,iret)
c       attributes of the grid
        name = 'Longitude on T grid     '; unit = 'degrees_W       '
        call ncaptc(ncid, Lon_tid, 'long_name', NCCHAR, 24, name, iret) 
        call ncaptc(ncid, Lon_tid, 'units',     NCCHAR, 16, unit, iret) 
        name = 'Longitude on U grid     '; unit = 'degrees_W       '
        call ncaptc(ncid, Lon_uid, 'long_name', NCCHAR, 24, name, iret) 
        call ncaptc(ncid, Lon_uid, 'units',     NCCHAR, 16, unit, iret) 
        name = 'Latitude on T grid      '; unit = 'degrees_N       '
        call ncaptc(ncid, Lat_tid, 'long_name', NCCHAR, 24, name, iret) 
        call ncaptc(ncid, Lat_tid, 'units',     NCCHAR, 16, unit, iret) 
        name = 'Latitude on U grid      '; unit = 'degrees_N       '
        call ncaptc(ncid, Lat_uid, 'long_name', NCCHAR, 24, name, iret) 
        call ncaptc(ncid, Lat_uid, 'units',     NCCHAR, 16, unit, iret) 
        name = 'Depth of T Grid points  '; unit = 'm '
        call ncaptc(ncid, depth_tid, 'long_name',NCCHAR, 24, name, iret) 
        call ncaptc(ncid, depth_tid, 'units',    NCCHAR, 16, unit, iret) 
        name = 'Depth of W Grid points  '; unit = 'm '
        call ncaptc(ncid, depth_wid, 'long_name', NCCHAR,24, name, iret) 
        call ncaptc(ncid, depth_wid, 'units',     NCCHAR,16, unit, iret) 
        name = 'Time                    '; unit = 'days            '
        call ncaptc(ncid, iTimeid, 'long_name', NCCHAR, 24, name, iret) 
        call ncaptc(ncid, iTimeid, 'units',     NCCHAR, 16, unit, iret) 
        call ncaptc(ncid, iTimeid,'time_origin',NCCHAR, 20,
     &  '31-DEC-1899 00:00:00', iret)
c       variables
        dims=(/lon_udim,lat_udim,depth_tdim,itimedim/)
        uvarid= ncvdef (ncid,'uvar' , NCFLOAT,4,dims,iret)
        umeanid=ncvdef (ncid,'umean', NCFLOAT,4,dims,iret)
        vvarid= ncvdef (ncid,'vvar' , NCFLOAT,4,dims,iret)
        vmeanid=ncvdef (ncid,'vmean', NCFLOAT,4,dims,iret)
        uvid  = ncvdef (ncid,'uv'   , NCFLOAT,4,dims,iret)
        T4id  = ncvdef (ncid,'T4'   , NCFLOAT,4,dims,iret)

c       attributes of variables
        name = 'Variance of zonal velocity'; unit = 'm^2/s^2'
        call dvcdf(ncid,uvarid,name,24,unit,16,spval)
        name = 'Mean of zonal velocity'; unit = 'm/s'
        call dvcdf(ncid,umeanid,name,24,unit,16,spval)
        name = 'Variance of meridional velocity'; unit = 'm^2/s^2'
        call dvcdf(ncid,vvarid,name,24,unit,16,spval)
        name = 'Mean of meridional velocity'; unit = 'm/s'
        call dvcdf(ncid,vmeanid,name,24,unit,16,spval)
        name = 'Hor. velocity correlation'; unit = 'm^2/s^2'
        call dvcdf(ncid,uvid,name,24,unit,16,spval)

        name = 'EKE-> MKE'; unit = 'm^2/s^3'
        call dvcdf(ncid,T4id,name,24,unit,16,spval)

c       global attrubutes
        iret= nf_put_att_int(ncid,nf_global,
     &   'Number_of_averaged_time_steps ',nf_int,1,counter)
        call ncendf(ncid, iret)

c       write grid
        allocate( v2(max(imt,jmt,km)) )
        v2(1:imt)=xt
        call ncvpt(ncid, Lon_tid, 1, imt,v2, iret)
        v2(1:imt)=xu
        call ncvpt(ncid, Lon_uid, 1, imt,v2, iret)
        v2(1:jmt)=yt
        call ncvpt(ncid, Lat_tid, 1, jmt,v2, iret)
        v2(1:jmt)=yu
        call ncvpt(ncid, Lat_uid, 1, jmt,v2, iret)
        v2(1:km)=zt/100.
        call ncvpt(ncid, depth_tid, 1, km,v2, iret)
        v2(1:km)=zw/100.
        call ncvpt(ncid, depth_wid, 1, km,v2, iret)
        call read_stamp(current_stamp,tt ); v2(1)=tt
        call ncvpt(ncid, itimeid, 1, 1,v2, iret)
        deallocate(v2)
        call ncclos (ncid, iret)
       endif

c      loop over PEs
       do npe=0,n_pes
        call barrier
        if (my_pe==npe) then
         iret=nf_open(fname,NF_WRITE,ncid)
         iret=nf_set_fill(ncid, NF_NOFILL, iret)
         iret=nf_inq_varid(ncid,'umean',umeanid)
         iret=nf_inq_varid(ncid,'uvar',uvarid)
         iret=nf_inq_varid(ncid,'vmean',vmeanid)
         iret=nf_inq_varid(ncid,'vvar',vvarid)
         iret=nf_inq_varid(ncid,'uv',uvid)
         iret=nf_inq_varid(ncid,'T4',T4id)

         where( umask(is:ie,1:km,js:je) == 0. ) 
          umean(is:ie,:,js:je,1) = spval
          umean(is:ie,:,js:je,2) = spval
          uvar(is:ie,:,js:je,1) = spval
          uvar(is:ie,:,js:je,2) = spval
          uv(is:ie,:,js:je) = spval
          T4(is:ie,:,js:je) = spval
         end where

         do k=1,km
          start=(/is,js,k,1/); count=(/ie-is+1,je-js+1,1,1/)
          var=umean(:,k,:,1)
          iret= nf_put_vara_real (ncid,umeanid,start,count,var )
          var=umean(:,k,:,2)
          iret= nf_put_vara_real (ncid,vmeanid,start,count,var )
          var=uvar(:,k,:,1)
          iret= nf_put_vara_real (ncid,uvarid,start,count,var )
          var=uvar(:,k,:,2)
          iret= nf_put_vara_real (ncid,vvarid,start,count,var )
          var=uv(:,k,:)
          iret= nf_put_vara_real (ncid,uvid,start,count,var )
          var=T4(:,k,:)
          iret= nf_put_vara_real (ncid,t4id,start,count,var )
         enddo
         call ncclos (ncid, iret)
        endif
        call barrier
       enddo


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

       if (my_pe == 0) then
        print*,' Writing variances to  NetCDF output file ',
     &         fname(1:len_trim(fname))
        ncid = nccre (fname, NCCLOB, iret)
        iret=nf_set_fill(ncid, NF_NOFILL, iret)
        call store_info_cdf(ncid)
c       dimensions
        lon_tdim  = ncddef(ncid, 'Longitude_t', imt, iret)
        Lon_udim  = ncddef(ncid, 'Longitude_u', imt, iret)
        Lat_tdim  = ncddef(ncid, 'Latitude_t',  jmt, iret)
        Lat_udim  = ncddef(ncid, 'Latitude_u',  jmt, iret)
        depth_wdim = ncddef(ncid, 'depth_w',  km, iret)
        depth_tdim = ncddef(ncid, 'depth_t',  km, iret)
        iTimedim  = ncddef(ncid, 'Time', 1, iret)
c       grid variables
        dims(1)  = Lon_tdim
        Lon_tid  = ncvdef (ncid,'Longitude_t',NCFLOAT,1,dims,iret)
        dims(1)  = Lon_udim
        Lon_uid  = ncvdef (ncid,'Longitude_u',NCFLOAT,1,dims,iret)
        dims(1)  = Lat_tdim
        Lat_tid  = ncvdef (ncid,'Latitude_t', NCFLOAT,1,dims,iret)
        dims(1)  = Lat_udim
        Lat_uid  = ncvdef (ncid,'Latitude_u', NCFLOAT,1,dims,iret)
        dims(1)  = iTimedim
        iTimeid   = ncvdef(ncid,'Time',       NCFLOAT,1,dims,iret)
        dims(1)  = depth_wdim
        depth_wid = ncvdef (ncid,'depth_w', NCFLOAT,1,dims,iret)
        dims(1)  = depth_tdim
        depth_tid = ncvdef (ncid,'depth_t', NCFLOAT,1,dims,iret)
c       attributes of the grid
        name = 'Longitude on T grid     '; unit = 'degrees_W       '
        call ncaptc(ncid, Lon_tid, 'long_name', NCCHAR, 24, name, iret) 
        call ncaptc(ncid, Lon_tid, 'units',     NCCHAR, 16, unit, iret) 
        name = 'Longitude on U grid     '; unit = 'degrees_W       '
        call ncaptc(ncid, Lon_uid, 'long_name', NCCHAR, 24, name, iret) 
        call ncaptc(ncid, Lon_uid, 'units',     NCCHAR, 16, unit, iret) 
        name = 'Latitude on T grid      '; unit = 'degrees_N       '
        call ncaptc(ncid, Lat_tid, 'long_name', NCCHAR, 24, name, iret) 
        call ncaptc(ncid, Lat_tid, 'units',     NCCHAR, 16, unit, iret) 
        name = 'Latitude on U grid      '; unit = 'degrees_N       '
        call ncaptc(ncid, Lat_uid, 'long_name', NCCHAR, 24, name, iret) 
        call ncaptc(ncid, Lat_uid, 'units',     NCCHAR, 16, unit, iret) 
        name = 'Depth of T Grid points  '; unit = 'm '
        call ncaptc(ncid, depth_tid, 'long_name',NCCHAR, 24, name, iret) 
        call ncaptc(ncid, depth_tid, 'units',    NCCHAR, 16, unit, iret) 
        name = 'Depth of W Grid points  '; unit = 'm '
        call ncaptc(ncid, depth_wid, 'long_name', NCCHAR,24, name, iret) 
        call ncaptc(ncid, depth_wid, 'units',     NCCHAR,16, unit, iret) 
        name = 'Time                    '; unit = 'days            '
        call ncaptc(ncid, iTimeid, 'long_name', NCCHAR, 24, name, iret) 
        call ncaptc(ncid, iTimeid, 'units',     NCCHAR, 16, unit, iret) 
        call ncaptc(ncid, iTimeid,'time_origin',NCCHAR, 20,
     &  '31-DEC-1899 00:00:00', iret)
c       variables
        dims=(/lon_tdim,lat_tdim,depth_tdim,itimedim/)
        T2id  = ncvdef (ncid,'T2'   , NCFLOAT,4,dims,iret)
        bmid  = ncvdef (ncid,'b'   , NCFLOAT,4,dims,iret)
        dims=(/lon_udim,lat_tdim,depth_tdim,itimedim/)
        ubid  = ncvdef (ncid,'ub'   , NCFLOAT,4,dims,iret)
        dims=(/lon_tdim,lat_udim,depth_tdim,itimedim/)
        vbid  = ncvdef (ncid,'vb'   , NCFLOAT,4,dims,iret)
        dims=(/lon_tdim,lat_tdim,depth_wdim,itimedim/)
        wbid  = ncvdef (ncid,'wb'   , NCFLOAT,4,dims,iret)

c       attributes of variables
        name = 'EPE-> MPE'; unit = 'm^2/s^3'
        call dvcdf(ncid,T2id,name,24,unit,16,spval)
        name = 'Zonal eddy buoyancy flux'; unit = 'm^2/s^3'
        call dvcdf(ncid,ubid,name,24,unit,16,spval)
        name = 'Meridional eddy buoyancy flux'; unit = 'm^2/s^3'
        call dvcdf(ncid,vbid,name,24,unit,16,spval)
        name = 'Vertical eddy buoyancy flux'; unit = 'm^2/s^3'
        call dvcdf(ncid,wbid,name,24,unit,16,spval)
        name = 'Mean buoyancy'; unit = 'm/s^2'
        call dvcdf(ncid,bmid,name,24,unit,16,spval)

c       global attrubutes
        iret= nf_put_att_int(ncid,nf_global,
     &   'Number_of_averaged_time_steps ',nf_int,1,counter)
        call ncendf(ncid, iret)

c       write grid
        allocate( v2(max(imt,jmt,km)) )
        v2(1:imt)=xt
        call ncvpt(ncid, Lon_tid, 1, imt,v2, iret)
        v2(1:imt)=xu
        call ncvpt(ncid, Lon_uid, 1, imt,v2, iret)
        v2(1:jmt)=yt
        call ncvpt(ncid, Lat_tid, 1, jmt,v2, iret)
        v2(1:jmt)=yu
        call ncvpt(ncid, Lat_uid, 1, jmt,v2, iret)
        v2(1:km)=zt/100.
        call ncvpt(ncid, depth_tid, 1, km,v2, iret)
        v2(1:km)=zw/100.
        call ncvpt(ncid, depth_wid, 1, km,v2, iret)
        call read_stamp(current_stamp,tt ); v2(1)=tt
        call ncvpt(ncid, itimeid, 1, 1,v2, iret)
        deallocate(v2)
        call ncclos (ncid, iret)
       endif

c      loop over PEs
       do npe=0,n_pes
        call barrier
        if (my_pe==npe) then
         iret=nf_open(fname,NF_WRITE,ncid)
         iret=nf_set_fill(ncid, NF_NOFILL, iret)
         iret=nf_inq_varid(ncid,'T2',T2id)
         iret=nf_inq_varid(ncid,'ub',ubid)
         iret=nf_inq_varid(ncid,'vb',vbid)
         iret=nf_inq_varid(ncid,'wb',wbid)
         iret=nf_inq_varid(ncid,'b', bmid)

         where( tmask(is:ie,1:km,js:je) == 0. ) 
          bmean(is:ie,:,js:je) = spval
         end where
         where( umask(is:ie,1:km,js:je) == 0. ) 
          mean_eddy_flux(is:ie,:,js:je,1)=spval
          mean_eddy_flux(is:ie,:,js:je,2)=spval
          mean_eddy_flux(is:ie,:,js:je,3)=spval
          T2(is:ie,:,js:je)=spval
         end where

         do k=1,km
          start=(/is,js,k,1/); count=(/ie-is+1,je-js+1,1,1/)
          var(is:ie,js:je)=bmean(is:ie,k,js:je)
          iret= nf_put_vara_real(ncid,bmid,start,count,var )
          if (iret.ne.0) print*,nf_strerror(iret)
          var(is:ie,js:je)=mean_eddy_flux(is:ie,k,js:je,1)
          iret= nf_put_vara_real(ncid,ubid,start,count,var )
          var(is:ie,js:je)=mean_eddy_flux(is:ie,k,js:je,2)
          iret= nf_put_vara_real(ncid,vbid,start,count,var )
          var(is:ie,js:je)=mean_eddy_flux(is:ie,k,js:je,3)
          iret= nf_put_vara_real(ncid,wbid,start,count,var )
          var(is:ie,js:je)=T2(is:ie,k,js:je)
          iret= nf_put_vara_real(ncid,t2id,start,count,var )
         enddo
         call ncclos (ncid, iret)
        endif
        call barrier
       enddo



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

       if (my_pe == 0) then
        print*,' Writing variances to  NetCDF output file ',
     &         fname(1:len_trim(fname))
        ncid = nccre (fname, NCCLOB, iret)
        iret=nf_set_fill(ncid, NF_NOFILL, iret)
        call store_info_cdf(ncid)
c       dimensions
        lon_tdim  = ncddef(ncid, 'Longitude_t', imt, iret)
        Lon_udim  = ncddef(ncid, 'Longitude_u', imt, iret)
        Lat_tdim  = ncddef(ncid, 'Latitude_t',  jmt, iret)
        Lat_udim  = ncddef(ncid, 'Latitude_u',  jmt, iret)
        depth_wdim = ncddef(ncid, 'depth_w',  km, iret)
        depth_tdim = ncddef(ncid, 'depth_t',  km, iret)
        iTimedim  = ncddef(ncid, 'Time', 1, iret)
c       grid variables
        dims(1)  = Lon_tdim
        Lon_tid  = ncvdef (ncid,'Longitude_t',NCFLOAT,1,dims,iret)
        dims(1)  = Lon_udim
        Lon_uid  = ncvdef (ncid,'Longitude_u',NCFLOAT,1,dims,iret)
        dims(1)  = Lat_tdim
        Lat_tid  = ncvdef (ncid,'Latitude_t', NCFLOAT,1,dims,iret)
        dims(1)  = Lat_udim
        Lat_uid  = ncvdef (ncid,'Latitude_u', NCFLOAT,1,dims,iret)
        dims(1)  = iTimedim
        iTimeid   = ncvdef(ncid,'Time',       NCFLOAT,1,dims,iret)
        dims(1)  = depth_wdim
        depth_wid = ncvdef (ncid,'depth_w', NCFLOAT,1,dims,iret)
        dims(1)  = depth_tdim
        depth_tid = ncvdef (ncid,'depth_t', NCFLOAT,1,dims,iret)
c       attributes of the grid
        name = 'Longitude on T grid     '; unit = 'degrees_W       '
        call ncaptc(ncid, Lon_tid, 'long_name', NCCHAR, 24, name, iret) 
        call ncaptc(ncid, Lon_tid, 'units',     NCCHAR, 16, unit, iret) 
        name = 'Longitude on U grid     '; unit = 'degrees_W       '
        call ncaptc(ncid, Lon_uid, 'long_name', NCCHAR, 24, name, iret) 
        call ncaptc(ncid, Lon_uid, 'units',     NCCHAR, 16, unit, iret) 
        name = 'Latitude on T grid      '; unit = 'degrees_N       '
        call ncaptc(ncid, Lat_tid, 'long_name', NCCHAR, 24, name, iret) 
        call ncaptc(ncid, Lat_tid, 'units',     NCCHAR, 16, unit, iret) 
        name = 'Latitude on U grid      '; unit = 'degrees_N       '
        call ncaptc(ncid, Lat_uid, 'long_name', NCCHAR, 24, name, iret) 
        call ncaptc(ncid, Lat_uid, 'units',     NCCHAR, 16, unit, iret) 
        name = 'Depth of T Grid points  '; unit = 'm '
        call ncaptc(ncid, depth_tid, 'long_name',NCCHAR, 24, name, iret) 
        call ncaptc(ncid, depth_tid, 'units',    NCCHAR, 16, unit, iret) 
        name = 'Depth of W Grid points  '; unit = 'm '
        call ncaptc(ncid, depth_wid, 'long_name', NCCHAR,24, name, iret) 
        call ncaptc(ncid, depth_wid, 'units',     NCCHAR,16, unit, iret) 
        name = 'Time                    '; unit = 'days            '
        call ncaptc(ncid, iTimeid, 'long_name', NCCHAR, 24, name, iret) 
        call ncaptc(ncid, iTimeid, 'units',     NCCHAR, 16, unit, iret) 
        call ncaptc(ncid, iTimeid,'time_origin',NCCHAR, 20,
     &  '31-DEC-1899 00:00:00', iret)
c       variables
        dims=(/lon_tdim,lat_tdim,depth_tdim,itimedim/)
        bbid  = ncvdef (ncid,'bb'   , NCFLOAT,4,dims,iret)
        dims=(/lon_udim,lat_tdim,depth_tdim,itimedim/)
        ubbid  = ncvdef (ncid,'ubb'   , NCFLOAT,4,dims,iret)
        dims=(/lon_tdim,lat_udim,depth_tdim,itimedim/)
        vbbid  = ncvdef (ncid,'vbb'   , NCFLOAT,4,dims,iret)
        dims=(/lon_tdim,lat_tdim,depth_wdim,itimedim/)
        wbbid  = ncvdef (ncid,'wbb'   , NCFLOAT,4,dims,iret)

c       attributes of variables
        name = 'Zonal variance buoyancy flux'; unit = 'm^3/s^5'
        call dvcdf(ncid,ubbid,name,24,unit,16,spval)
        name = 'Meridional variance buoyancy flux'; unit = 'm^3/s^5'
        call dvcdf(ncid,vbbid,name,24,unit,16,spval)
        name = 'Vertical variance buoyancy flux'; unit = 'm^3/s^5'
        call dvcdf(ncid,wbbid,name,24,unit,16,spval)
        name = 'Buoyancy variance'; unit = 'm^2/s^4'
        call dvcdf(ncid,bbid,name,24,unit,16,spval)

c       global attrubutes
        iret= nf_put_att_int(ncid,nf_global,
     &   'Number_of_averaged_time_steps ',nf_int,1,counter)
        call ncendf(ncid, iret)

c       write grid
        allocate( v2(max(imt,jmt,km)) )
        v2(1:imt)=xt
        call ncvpt(ncid, Lon_tid, 1, imt,v2, iret)
        v2(1:imt)=xu
        call ncvpt(ncid, Lon_uid, 1, imt,v2, iret)
        v2(1:jmt)=yt
        call ncvpt(ncid, Lat_tid, 1, jmt,v2, iret)
        v2(1:jmt)=yu
        call ncvpt(ncid, Lat_uid, 1, jmt,v2, iret)
        v2(1:km)=zt/100.
        call ncvpt(ncid, depth_tid, 1, km,v2, iret)
        v2(1:km)=zw/100.
        call ncvpt(ncid, depth_wid, 1, km,v2, iret)
        call read_stamp(current_stamp,tt ); v2(1)=tt
        call ncvpt(ncid, itimeid, 1, 1,v2, iret)
        deallocate(v2)
        call ncclos (ncid, iret)
       endif

c      loop over PEs
       do npe=0,n_pes
        call barrier
        if (my_pe==npe) then
         iret=nf_open(fname,NF_WRITE,ncid)
         iret=nf_set_fill(ncid, NF_NOFILL, iret)
         iret=nf_inq_varid(ncid,'ubb',ubbid)
         iret=nf_inq_varid(ncid,'vbb',vbbid)
         iret=nf_inq_varid(ncid,'wbb',wbbid)
         iret=nf_inq_varid(ncid,'bb', bbid)

         where( tmask(is:ie,1:km,js:je) == 0. ) 
          var_tracer(is:ie,:,js:je) = spval
         end where
         where( umask(is:ie,1:km,js:je) == 0. ) 
          mean_var_flux(is:ie,:,js:je,1)=spval
          mean_var_flux(is:ie,:,js:je,2)=spval
          mean_var_flux(is:ie,:,js:je,3)=spval
         end where

         do k=1,km
          start=(/is,js,k,1/); count=(/ie-is+1,je-js+1,1,1/)
          var(is:ie,js:je)=var_tracer(is:ie,k,js:je)
          iret= nf_put_vara_real (ncid,bbid,start,count,var )
          var(is:ie,js:je)=mean_var_flux(is:ie,k,js:je,1)
          iret= nf_put_vara_real (ncid,ubbid,start,count,var )
          var(is:ie,js:je)=mean_var_flux(is:ie,k,js:je,2)
          iret= nf_put_vara_real (ncid,vbbid,start,count,var )
          var(is:ie,js:je)=mean_var_flux(is:ie,k,js:je,3)
          iret= nf_put_vara_real (ncid,wbbid,start,count,var )
         enddo
         call ncclos (ncid, iret)
        endif
        call barrier
       enddo

       if (my_pe==0) print*,'done'
#else
       if (my_pe==0) then
         print*,' ERROR : writing variances without netcdf'
         print*,'   is not yet implemented'
       endif
       call halt_stop(' in write_variances')
#endif
       end subroutine write_variances




