#include "options.inc"


      module diag_blue_module
      implicit none
      integer :: counter
      real, allocatable, dimension(:,:,:) :: corr
      end module diag_blue_module


      subroutine diag_blue_init() 
      use spflame_module
      use blue_module
      use diag_blue_module
      implicit none
      character (len=80) :: name
      integer :: io,k
      real, allocatable ::  buf(:,:)

      if (my_pe==0) print*,' initializing diagnostics for BLUE'

      if (.not.enable_blue) then
       if (my_pe==0) then
        print*,' '
        print*,' WARNING: '
        print*,' you have to switch on BLUE in order to use this option'
        print*,' '
       endif
       return
      endif

      allocate( corr(is_pe:ie_pe,km,js_pe:je_pe) )
      counter = 0; corr = 0.

      write(name,'("blue_averages_",i3,".dta")') sub_domain
      call replace_space_zero(name)
      if (my_pe==0) then
        print*,''
        print*,' reading unfinished averages for BLUE from file ',
     &          name(1:len_trim(name))
      endif

      call get_free_iounit(io)
      open(io,file=name,form='unformatted',status='old',err=20)
      read(io) counter
      allocate( buf(imt,jmt) )
      do k=1,km
       if (my_pe==0) read(io) buf
       call pe0_send_2D(buf)
       corr(is_pe:ie_pe,k,js_pe:je_pe)=buf(is_pe:ie_pe,js_pe:je_pe)
      enddo
      close(io)
      deallocate(buf)
      if (my_pe==0) print*,' done'

      return

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




      subroutine diag_blue() 
      use spflame_module
      use blue_module
      use diag_blue_module
      implicit none
      character (len=80) :: name
      integer :: io,k,i,j
      integer :: year,month,day
      real, allocatable ::  buf(:,:)

      if (.not.enable_blue) return
c
c     average
c
      counter=counter+1

      if (enable_blue_smooth) then

       do j=js_pe,je_pe
        do k=1,km
         do i=is_pe,ie_pe
          corr(i,k,j)=corr(i,k,j)+blue_rhom(i,k,j)
         enddo
        enddo
       enddo

      else

       do j=js_pe,je_pe
        do k=1,km
         do i=is_pe,ie_pe
          corr(i,k,j)=corr(i,k,j)+
     &        (1.-blue_alpha(i,k,j))*(
     &                blue_rho(i,k,j,blue_p1)*blue_f1+
     &                blue_rho(i,k,j,blue_p2)*blue_f2  - rho(i,k,j) )
         enddo
        enddo
       enddo

      endif
c
c     check for end of month, then write out file and initialize again
c
      if   (end_of_month )  then

#ifdef netcdf_diagnostics
c      not yet implemented
#else
       year  = get_year(get_current_time())
       month = get_month_of_year(get_current_time())
       day   = get_day_of_month(get_current_time())
       write(name,'("blue_averages_",i3,"_y",i4,"m",i2,"d",i2,".dta")') 
     &            sub_domain,year,month,day
       call replace_space_zero(name)

       if (my_pe==0) then
         print*,' writing BLUE averages to file ',name(1:len_trim(name))
         call getunit(io, name,'u s r ieee')
         write(io) imt,jmt,km
         write(io) xt,yt,zt
         write(io) counter
         write(io) kmt_big
       endif
       allocate( buf(imt,jmt) ); buf=0.
       do k=1,km
         buf(is_pe:ie_pe,js_pe:je_pe)=corr(is_pe:ie_pe,k,js_pe:je_pe)
         call pe0_recv_2D(buf)
         if (my_pe==0) write(io) buf
       enddo
       deallocate(buf)
       if (my_pe==0) then
        close(io)
        print*,'done'
       endif
#endif

       counter=0; corr=0.
      endif
c
c     write unfinished averages to be read in the next run
c
      if (last_time_step) then

       write(name,'("blue_averages_",i3,".dta")') sub_domain
       call replace_space_zero(name)
       if (my_pe==0) then
        print*,''
        print*,' writing unfinished averages for BLUE from 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 k=1,km
         buf(is_pe:ie_pe,js_pe:je_pe)=corr(is_pe:ie_pe,k,js_pe:je_pe)
         call pe0_recv_2D(buf)
         if (my_pe==0) write(io) buf
       enddo
       deallocate(buf)

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

      endif

      end subroutine diag_blue 
