#include "options.inc"
c
c
c see also below for a second module blue_mean_module
c
c-----------------------------------------------------------------------
c     This module implements the BLUE approach
c     by J.Sheng and R.J. Greatbatch, JGR, 2000
c     to get a semi diagnostic model in certain
c     regions of the model domain dependend on 
c     the parameter blue_alpha.
c
c     if blue_alpha = 1.  fully prognostic
c     if blue_alpha = 0.5 semi prognostic or diagnostic
c     if blue_alpha = 0.  fully diagnostic
c
c     additional options: 
c
c        use constant (diagnosed) forcing
c        smooth climatological rho with a grid point smoother
c
c     linked in the code in setup, call of blue_init
c                        in clinic, call of add_blue_to_rho
c                        in domain_exchg, send/rec. blue_rho
c
c     four parameters can be set in spflame_module (namelist)
c         enable_blue, 
c         enable_blue_smooth, 
c         enable_blue_const, 
c         enable_blue_tapered, 
c         ismooth_blue
c         ilook_blue
c         enable_blue_simple
c         enable_diag_blue
c
c                                      c eden
c-----------------------------------------------------------------------
c
c
c local defines
c
c#define prognostic_below_200m
c
c

      module blue_module
      implicit none

      real, allocatable :: blue_rho  (:,:,:,:)
      real, allocatable :: blue_alpha(:,:,:)
      real, allocatable :: blue_rhom (:,:,:)

      character (len=80) :: blue_temp_file
      character (len=80) :: blue_salt_file

      integer :: blue_interp_index
      integer :: blue_cache_pointer(2)
      real    :: blue_f1,blue_f2
      integer :: blue_p1,blue_p2

      end module blue_module


      subroutine blue_init
c-----------------------------------------------------------------------
c     initialisation for BLUE
c-----------------------------------------------------------------------
      use spflame_module
      use blue_module
      implicit none
      type( time_type) :: start
      type( time_type), allocatable :: interval(:)
      integer :: io,len,n
      real, allocatable  :: aprec(:)
      character (len=32) :: stamp
c
c     hardwired file names
c 
      blue_temp_file = 'blue_temp.mom'
      if (enable_blue_const) blue_temp_file = 'blue_const.mom'
      blue_salt_file = 'blue_salt.mom'
c
c     informative output about configuration
c
      if (my_pe == 0) then
       print*,  ''
       print*,  ' Initialization of BLUE'
       print*,  ''
       if (mother < 0) then 
        if (enable_blue_simple) then
         print*,' using simple T/S fields for BLUE, as specified'
         print*,' in file setup_template.F'
        else
         if (enable_blue_const) then
          print*,  ' reading momentum corr. from file ',
     &     blue_temp_file(1:len_trim(blue_temp_file))
         else
          print*,  ' reading temp/sal from files ',
     &     blue_temp_file(1:len_trim(blue_temp_file)),' and ',
     &     blue_salt_file(1:len_trim(blue_salt_file))
         endif
        endif
       else
        print*,' receiving blue density from mother domain '
       endif

       if (enable_blue_smooth) then
        print*,' using the smoothed version of BLUE'
        print*,' ismooth_blue =',ismooth_blue
       endif
       if (enable_blue_const) then
        print*,' using the constant forcing '
       endif
       if (enable_blue_tapered) then
        print*,' using boundary tapering '
        print*,' ilook_blue = ',ilook_blue
       endif
       print*,  ''
       if (enable_blue_smooth .and. enable_blue_const) then
        print*,' smoothed version and constant forcing of BLUE'
        print*,' are incompatible '
       endif

      endif
c
c     checks
c
      if (enable_blue_smooth .and. enable_blue_const) then
        call halt_stop(' in blue_init')
      endif
c
c     allocate permanent memory
c
      allocate( blue_rho(is_pe-1:ie_pe+1,km,js_pe-1:je_pe+1,2) )
      allocate( blue_alpha(is_pe-1:ie_pe+1,km,js_pe-1:je_pe+1) )
      blue_rho = 0.; blue_alpha = 1.

      if (enable_blue_smooth) then
       allocate( blue_rhom(is_pe-1:ie_pe+1,km,js_pe-1:je_pe+1) )
       blue_rhom=0.
      endif

      if (mother < 0) then 

       if (.not. enable_blue_simple) then
c
c       initialize time interpolation for blue density
c
        call getunit(io,blue_temp_file,'usr ieee')
        read (io) len
        rewind(io)
        allocate(interval(len), aprec(len))
        read (io)  len,stamp,aprec

        do n=1,len
         call set_time(interval(n),
     &                 int((aprec(n)-int(aprec(n)))*24*60*60),
     &                 int(aprec(n)))
        enddo
        start= get_stamp (stamp)
        call init_forcing_interp(len,.true.,
     &                     start,interval,blue_interp_index)
        deallocate(interval,aprec)
        close(io)

       else
c
c       use simple setup here
c
        call blue_template()

       endif

      else
c
c       child domains receive their value from the mother domain
c
      endif
c
c     set the BLUE  parameter
c     (see also sub_domain_exchg for BLUE and BLUE parameter)
c
      blue_alpha=0.5

#ifdef prognostic_below_200m
c
c     use BLUE only above 200m
c
      blue_alpha(:,13:km,:)=1.0
#endif

      if (enable_blue_tapered) call blue_alpha_tapering()

      if (my_pe == 0) print*,  'done'
      end subroutine blue_init




      subroutine blue_alpha_tapering
      use spflame_module
      use blue_module
      implicit none
      integer :: i,k,j,ii,jj,iii,jjj

c      do j=js_pe-1,je_pe+1
c       do k=1,km
c        do i=is_pe-1,ie_pe+1
c         if (kmt_big(i,j)<k) blue_alpha(i,k,j)=1. ! land
c        enddo
c       enddo
c      enddo

      do j=js_pe-1,je_pe+1
       do k=1,km
        do i=is_pe-1,ie_pe+1
         ii=i
         jj=j
         if (cyclic) then
            if (ii>imt) ii=ii-imt+2
            if (ii<1)   ii=ii+imt-2
         else
            ii = min(imt,ii)
            ii = max(1  ,ii)
         endif
         jj = min(jmt,jj)
         jj = max(1  ,jj)
         if (kmt_big(ii,jj)<k) blue_alpha(i,k,j)=1. ! land
        enddo
       enddo
      enddo


      do jj=-ilook_blue,ilook_blue
       do ii=-ilook_blue,ilook_blue
        do j=js_pe-1,je_pe+1
         do k=1,km
          do i=is_pe-1,ie_pe+1
           iii=ii+i
           jjj=jj+j
           if (cyclic) then
            if (iii>imt) iii=iii-imt+2
            if (iii<1)   iii=iii+imt-2
           else
            iii = min(imt,iii)
            iii = max(1  ,iii)
           endif
           jjj = min(jmt,jjj)
           jjj = max(1  ,jjj)
           ! near land and above land
           if (kmt_big(iii,jjj)<=k) blue_alpha(i,k,j)=1. 
          enddo
         enddo
        enddo
       enddo
      enddo
      end subroutine blue_alpha_tapering



      subroutine read_blue
      use spflame_module
      use blue_module
      implicit none
c-----------------------------------------------------------------------
c     read in BLUE T/S fields from binary files
c     (or momentum corr. for blue_const version
c-----------------------------------------------------------------------
      real    :: f1,f2
      integer :: p1,p2,io1,io2,n,i,k,j
      logical load(2)
      character(len=80) :: text = 'BLUE'
      real, allocatable :: temp(:,:),salt(:,:)
      integer :: is,ie,js,je

      is=max(1,is_pe-1); ie = min(imt,ie_pe+1)
      js=max(1,js_pe-1); je = min(jmt,je_pe+1)

      if (mother < 0) then 
c
c      decide whether to read data or not
c
       call forcing_interp(blue_interp_index,
     &                     blue_p1,blue_p2,blue_f1,blue_f2)

       call update_forcing_cache(blue_interp_index,blue_p1,blue_p2,
     &                           blue_cache_pointer(:),load,text)

       do n=1,2
        if (load(n)) then
c
c        read in data
c
         if (enable_blue_const) then
c
c         read only from first file
c
          if (my_pe==0) then
           call getunit(io1,blue_temp_file,'usr ieee')
           read(io1)  ! len et al
           do k=1,blue_cache_pointer(n)-1
            do j=1,jmt
             read(io1)  ! skip all previous time steps
            enddo
           enddo
          endif ! my_pe==0
          allocate( temp(imt,km) )

          do j=1,jmt
           if (my_pe==0) read (io1) temp  ! these are the data for this time step
           call barrier
           call bcast_real(temp,imt*km,0) 
c
c          transfer to BLUE correction array
c          
           if (j>=js .and. j<=je) then
            do k=1,km
             do i=is,ie
              blue_rho(i,k,j,n)=temp(i,k)
             enddo
            enddo
           endif
          enddo ! do j=1,jmt

          deallocate( temp) 
          if (my_pe==0) close(io1)

         else ! if enable_blue_const
c
c         read T/S from two files
c
          if (my_pe==0) then
           call getunit(io1,blue_temp_file,'usr ieee')
           call getunit(io2,blue_salt_file,'usr ieee')
           read(io1)  ! len et al
           read(io2)  ! len et al
           do k=1,blue_cache_pointer(n)-1
            do j=1,jmt
             read(io1)  ! skip all previous time steps
             read(io2)
            enddo
           enddo
          endif ! my_pe==0

          allocate( temp(imt,km), salt(imt,km) )
          do j=1,jmt
           if (my_pe==0) then
            read (io1)    temp  ! these are the data for this time step
            read (io2)    salt
           endif
           ! assuming that broadcasting is quicker than file IO
           call barrier
           call bcast_real(temp,imt*km,0) 
           call bcast_real(salt,imt*km,0) 
c
c          now calculate BLUE density from T/S fields
c          
           if (j>=js .and. j<=je) then
            do k=1,km
             call model_dens(temp(is,k), salt(is,k),
     &                  blue_rho(is,k,j,n),k,ie-is+1
#ifdef partial_cell
     &                       ,ztp(is,k,j)
#endif
     &                       )
            enddo
           endif
          enddo ! do j=1,jmt

          deallocate( temp,salt) 
          if (my_pe==0) then
           close(io1); close(io2)
          endif

         endif ! if enable_blue_const
        endif ! if load
       enddo ! n

      else ! mother < 0
c
c      childs get their blue_rho in sub_domain_exghc
c
       blue_p1=1; blue_p2=1; blue_f1=1.0; blue_f2=0.0 

      endif
      end subroutine read_blue





      subroutine add_blue_to_rho()
c-----------------------------------------------------------------------
c     add BLUE density in different formulation to model density
c-----------------------------------------------------------------------
      use spflame_module
      use blue_module
      implicit none
      integer :: i,k,j
      real :: a

      if (.not. enable_blue_simple) then
c
c      consider to read in data or not
c
       call read_blue

      endif
 
      if (enable_blue_smooth) then
c
c      smoothed version, calculate smoothed density first
c
       call calc_smooth_blue
       do j=js_pe-1,je_pe+1
        do k=1,km
         do i=is_pe-1,ie_pe+1
          rho(i,k,j) = rho(i,k,j)+blue_rhom(i,k,j) 
         enddo
        enddo
       enddo
      elseif (enable_blue_const) then
c
c      constant forcing, just add to model density
c
       do j=js_pe-1,je_pe+1
        do k=1,km
         do i=is_pe-1,ie_pe+1
          rho(i,k,j) = rho(i,k,j) + 
     &               (blue_rho(i,k,j,blue_p1)*blue_f1+
     &                blue_rho(i,k,j,blue_p2)*blue_f2)
         enddo
        enddo
       enddo
      else
c
c      the normal version
c
       do j=js_pe-1,je_pe+1
        do k=1,km
         do i=is_pe-1,ie_pe+1
          a=blue_alpha(i,k,j)
          rho(i,k,j) = a*rho(i,k,j) + 
     &        (1-a)* (blue_rho(i,k,j,blue_p1)*blue_f1+
     &                blue_rho(i,k,j,blue_p2)*blue_f2)
         enddo
        enddo
       enddo
      endif
      end subroutine add_blue_to_rho



      subroutine calc_smooth_blue
c-----------------------------------------------------------------------
c     calculate a smoothed version of density correction
c-----------------------------------------------------------------------
      use spflame_module
      use blue_module
      implicit none
      real, allocatable :: blue_rhoc(:,:,:)
      integer :: i,j,k,jj,ii
c
c     allocate work space
c
      allocate( blue_rhoc(is_pe-ismooth_blue:ie_pe+ismooth_blue
     &       ,km,js_pe-ismooth_blue:je_pe+ismooth_blue) )
      blue_rhoc=0.
c
c     density correction can be formulated as
c       rho_star = (1-alpha) ( rho_clim - rho_model)
c
      do j=js_pe,je_pe
        do k=1,km
          do i=is_pe,ie_pe
            blue_rhoc(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
c
c     boundary exchanges
c
      call border_exchg(blue_rhoc,km,ismooth_blue)
      call set_cyclic_nx(blue_rhoc,km,ismooth_blue)
c
c     now smooth rho_star
c
      blue_rhom=0.
      do jj=-ismooth_blue,ismooth_blue
       do j=js_pe,je_pe
        do k=1,km
         do ii=-ismooth_blue,ismooth_blue
          do i=is_pe,ie_pe
           blue_rhom(i,k,j)=blue_rhom(i,k,j)+blue_rhoc(i+ii,k,j+jj)
          enddo
         enddo
        enddo
       enddo
      enddo
      blue_rhom=blue_rhom /( 2.*ismooth_blue+1.)**2
      deallocate(blue_rhoc)
c
c     again boundary exchange
c
      call border_exchg(blue_rhom,km,1)
      call set_cyclic(blue_rhom,km,1)
      end subroutine calc_smooth_blue



c
c
c-----------------------------------------------------------------------
c     This module implements the BLUE approach
c     by J.Sheng and R.J. Greatbatch, JGR, 2000
c     to get a semi diagnostic model in certain
c     regions of the model domain dependend on 
c     the parameter blue_alpha.
c
c     if blue_alpha = 1.  fully prognostic
c     if blue_alpha = 0.5 semi prognostic or diagnostic
c     if blue_alpha = 0.  fully diagnostic
c
c     additional options: 
c
c      smooth correction a grid point smoother
c      taper alpha at lateral boundaries
c
c     linked in the code in setup, call of blue_mean_init
c                        in clinic, call of add_blue_mean_to_rho
c                        in domain_exchg, send/rec. blue_rho
c                           (not yet implemented)
c
c     parameters can be set in spflame_module (namelist)
c
c        enable_blue_mean
c        enable_blue_smooth
c        enable_blue_tapered
c        ilook_blue, ismooth_blue
c
c                                      c eden
c-----------------------------------------------------------------------
c

      module blue_mean_module
      implicit none

      real, allocatable :: blue_rho  (:,:,:) ! the clim. density
      real, allocatable :: blue_alpha(:,:,:)   ! the parameter alpha
      real, allocatable :: blue_rhom (:,:,:) ! the mean correction
      real, allocatable :: model_rhom (:,:,:) ! the mean model density

      character (len=80) :: blue_file
      integer :: counter 

      end module blue_mean_module


      subroutine blue_mean_init
c-----------------------------------------------------------------------
c     initialisation for BLUE
c-----------------------------------------------------------------------
      use spflame_module
      use blue_mean_module
      implicit none
      integer :: io,js,je,is,ie,i,j,k
      real, allocatable :: temp(:,:),salt(:,:)
      character (len=80) :: name

      is=max(1,is_pe-1); ie = min(imt,ie_pe+1)
      js=max(1,js_pe-1); je = min(jmt,je_pe+1)
c
c     hardwired file names
c 
      blue_file = 'blue_mean.mom'

      if (my_pe == 0) then
       print*,  ''
       print*,  ' Initialization of BLUE mean'
       print*,  ''
       if (mother < 0) then 
        print*,  ' reading temp/sal from file'
        print*,  blue_file
       else
        print*,' receiving blue density from mother domain '
       endif

       if (enable_blue_smooth) then
        print*,' using the smoothed version of BLUE mean'
        print*,' ismooth_blue =',ismooth_blue
       endif
       if (enable_blue_tapered) then
        print*,' using boundary tapering '
        print*,' ilook_blue = ',ilook_blue
       endif

      endif

      if (mother>0) then
        print*,' sub domain does not work yet'
        call halt_stop(' in blue_mean_init')
      endif
c
c     allocate permanent memory
c
      allocate( blue_rho(is_pe-1:ie_pe+1,km,js_pe-1:je_pe+1) )
      allocate( blue_alpha(is_pe-1:ie_pe+1,km,js_pe-1:je_pe+1) )
      blue_rho = 0.; blue_alpha = 1.

      allocate( blue_rhom(is_pe-1:ie_pe+1,km,js_pe-1:je_pe+1) )
      blue_rhom=0.
      allocate( model_rhom(is_pe-1:ie_pe+1,km,js_pe-1:je_pe+1) )
      model_rhom=0.

      if (mother < 0) then 
c
c      read blue density
c
       call getunit(io,blue_file,'usr ieee')

       allocate( temp(imt,km), salt(imt,km) )

       do j=1,jmt

          if (my_pe==0) read (io)    temp  ,salt

          call barrier
          call bcast_real(temp,imt*km,0) ! assuming that broadcasting is quicker than 
          call bcast_real(salt,imt*km,0) ! file IO
c
c         now calculate BLUE density from T/S fields
c          
          if (j>=js .and. j<=je) then
           do k=1,km
            call model_dens(temp(is,k), salt(is,k),
     &                  blue_rho(is,k,j),k,ie-is+1
#ifdef partial_cell
     &                       ,ztp(is,k,j)
#endif
     &                       )
           enddo
          endif

       enddo
       deallocate( temp,salt) 

       close(io)

      else
c
c       child domains receive their value from the mother domain
c
      endif
c
c     set the BLUE  parameter
c     (see also sub_domain_exchg for BLUE and BLUE parameter)
c
      blue_alpha=0.5


#ifdef prognostic_below_200m
c
c     use BLUE only above 200m
c
c      blue_alpha(:,13:km,:)=1.0
c      blue_alpha=1.0; blue_alpha(:,13:24,:)=0.5
      blue_alpha=1.0; blue_alpha(:,1:24,:)=0.5
#endif

      if (enable_blue_tapered) then
        call blue_mean_alpha_tapering()
      endif
c
c     read written data from previous run (if any)
c
      counter=0; blue_rhom = 0.; model_rhom = 0.

      write(name,'("bluerest_",i3,".dta")') sub_domain
      do i=1,len_trim(name); if (name(i:i) == ' ') name(i:i)='0'; enddo

      if (my_pe==0) then
        print*,''
        print*,' reading averaged model density 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( temp(imt,jmt), salt(imt,jmt) )
      do k=1,km
          if (my_pe==0) read (io) temp,salt
          call pe0_send_2D(temp)
          call pe0_send_2D(salt)
          do j=js,je
            do i=is,ie
              model_rhom(i,k,j) = temp(i,j)
              blue_rhom(i,k,j)  = salt(i,j)
            enddo
          enddo
      enddo
      deallocate( temp, salt) 
      close(io)
c
      call border_exchg(blue_rhom,km,1)
      call set_cyclic(blue_rhom,km,1)
      call border_exchg(model_rhom,km,1)
      call set_cyclic(model_rhom,km,1)
c
      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 blue_mean_init



      subroutine add_blue_mean_to_rho()
c-----------------------------------------------------------------------
c     add BLUE density in different formulation to model density
c-----------------------------------------------------------------------
      use spflame_module
      use blue_mean_module
      implicit none
      integer :: io,js,je,is,ie,i,j,k
      real, allocatable :: temp(:,:),salt(:,:)
      character (len=80) :: name

      is=max(1,is_pe-1); ie = min(imt,ie_pe+1)
      js=max(1,js_pe-1); je = min(jmt,je_pe+1)
c
c     continue to average model density
c
      if (.not. eulerback1 ) then

        counter=counter+1

        do j=js_pe-1,je_pe+1
          do k=1,km
           do i=is_pe-1,ie_pe+1
             model_rhom(i,k,j) = model_rhom(i,k,j) + rho(i,k,j)
           enddo
          enddo
        enddo
 
      endif
c
c     consider to calculate new mean correction
c
      if (end_of_year ) call change_blue_mean()
c
c     add correction to model density
c
      do j=js_pe-1,je_pe+1
       do k=1,km
        do i=is_pe-1,ie_pe+1
          rho(i,k,j) = rho(i,k,j) + blue_rhom(i,k,j)
        enddo
       enddo
      enddo
c
c     write a restart if end of run is reached
c
      if (last_time_step) then

       write(name,'("bluerest_",i3,".dta")') sub_domain
       do i=1,len_trim(name); if (name(i:i) == ' ') name(i:i)='0'; enddo

       if (my_pe==0) then
        print*,''
        print*,' writing averaged model density to file ',
     &          name(1:len_trim(name))
        call get_free_iounit(io)
        open(io,file=name,form='unformatted')
        write(io) counter
       endif

       allocate( temp(imt,jmt), salt(imt,jmt) )
       do k=1,km
         do j=js,je
          do i=is,ie
            temp(i,j)=model_rhom(i,k,j)
            salt(i,j)=blue_rhom(i,k,j)
          enddo
         enddo
         call pe0_recv_2D(temp)
         call pe0_recv_2D(salt)
         if (my_pe==0) write (io) temp,salt
       enddo
       deallocate( temp, salt) 
       if (my_pe == 0) then
        print*,  'done'
        close(io)
       endif

      endif ! last time step

      end subroutine add_blue_mean_to_rho



      subroutine change_blue_mean
c-----------------------------------------------------------------------
c     calculate a smoothed version of density correction
c     at the end of each year
c-----------------------------------------------------------------------
      use spflame_module
      use blue_mean_module
      implicit none
      real, allocatable :: blue_rhoc(:,:,:)
      integer :: i,j,k,jj,ii

      integer :: io,year,month,day
      character (len=80) :: name
      real, allocatable ::  buf(:,:)

      if (enable_blue_smooth) then
c
c      allocate work space for smoothing
c
       allocate( blue_rhoc(is_pe-ismooth_blue:ie_pe+ismooth_blue,
     &           km,js_pe-ismooth_blue:je_pe+ismooth_blue) )
       blue_rhoc=0.
c
c      density correction can be formulated as
c       rho_star = (1-alpha) ( rho_clim - rho_model)
c
       do j=js_pe,je_pe
        do k=1,km
          do i=is_pe,ie_pe
            blue_rhoc(i,k,j) = (1-blue_alpha(i,k,j))*
     &      (  blue_rho(i,k,j) -model_rhom(i,k,j)/counter )
          enddo
        enddo
       enddo
c
c      boundary exchanges
c
       call border_exchg(blue_rhoc,km,ismooth_blue)
       call set_cyclic_nx(blue_rhoc,km,ismooth_blue)
c
c      now smooth rho_star
c
       blue_rhom=0.
       do jj=-ismooth_blue,ismooth_blue
        do j=js_pe,je_pe
         do k=1,km
          do ii=-ismooth_blue,ismooth_blue
           do i=is_pe,ie_pe
            blue_rhom(i,k,j)=blue_rhom(i,k,j)+blue_rhoc(i+ii,k,j+jj)
           enddo
          enddo
         enddo
        enddo
       enddo
       blue_rhom=blue_rhom /( 2.*ismooth_blue+1.)**2
       deallocate(blue_rhoc)
c
c      again boundary exchange
c
       call border_exchg(blue_rhom,km,1)
       call set_cyclic(blue_rhom,km,1)

      else
c
c      no smoothing
c
       do j=js_pe,je_pe
        do k=1,km
          do i=is_pe,ie_pe
            blue_rhom(i,k,j) = (1-blue_alpha(i,k,j))*
     &      (  blue_rho(i,k,j) -model_rhom(i,k,j)/counter )
          enddo
        enddo
       enddo
c
c      boundary exchanges
c
       call border_exchg(blue_rhom,km,1)
       call set_cyclic(blue_rhom,km,1)

      endif
c
c     write the new correction to another binary file
c

#ifdef netcdf_diagnostics
c      not yet implemented
#endif

      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_mean_corr_",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 mean correction 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)=blue_rhom(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
c
c     zerofy counter and average field
c
      counter = 0
      model_rhom = 0.
      end subroutine change_blue_mean




      subroutine blue_mean_alpha_tapering
      use spflame_module
      use blue_mean_module
      implicit none
      integer :: i,k,j,ii,jj,iii,jjj

      do j=js_pe-1,je_pe+1
       do k=1,km
        do i=is_pe-1,ie_pe+1
         ii=i
         jj=j
         if (cyclic) then
            if (ii>imt) ii=ii-imt+2
            if (ii<1)   ii=ii+imt-2
         else
            ii = min(imt,ii)
            ii = max(1  ,ii)
         endif
         jj = min(jmt,jj)
         jj = max(1  ,jj)
         if (kmt_big(ii,jj)<k) blue_alpha(i,k,j)=1. ! land
        enddo
       enddo
      enddo

      do jj=-ilook_blue,ilook_blue
       do ii=-ilook_blue,ilook_blue
        do j=js_pe-1,je_pe+1
         do k=1,km
          do i=is_pe-1,ie_pe+1
           iii=ii+i
           jjj=jj+j
           if (cyclic) then
            if (iii>imt) iii=iii-imt+2
            if (iii<1)   iii=iii+imt-2
           else
            iii = min(imt,iii)
            iii = max(1  ,iii)
           endif
           jjj = min(jmt,jjj)
           jjj = max(1  ,jjj)
           if (kmt_big(iii,jjj)<=k) blue_alpha(i,k,j)=1. ! near land and above land
          enddo
         enddo
        enddo
       enddo
      enddo
      end subroutine blue_mean_alpha_tapering


