program gen_be_cov3d3d_bin3d_contrib
!------------------------------------------------------------------------
!  Purpose: Computes contribution of "psi" Field on the balanced part
!           of "chi" Field for WRFDA cv_options=6
!
!  Auothor: Syed RH Rizvi (MMM/NESL/NCAR)   Date: 02/01/2010
!           & Monika Krysta, (CTBTO, Vienna, Austria)
!
!  Note: Please acknowledge author/institute in work that uses this code.
!------------------------------------------------------------------------

   use da_control, only : stderr, stdout, filename_len
   use da_tools_serial, only : da_get_unit,da_advance_cymdh

   implicit none

   character*10        :: start_date, end_date       ! Starting and ending dates.
   character*10        :: date, new_date             ! Current date (ccyymmddhh).
   character*80        :: variable
   character*80        :: regcoeff_var
   character*10        :: variable1                  ! Variable name
   character*10        :: variable2                  ! Variable name

   real, allocatable   :: regcoeff_psi_chi(:)        ! psi/chi    regression cooefficient.
   real, allocatable   :: regcoeff_psi_t(:,:,:)      ! psi/t      regression cooefficient.
   real, allocatable   :: regcoeff_psi_ps(:,:)       ! psi/ps     regression cooefficient.
   real, allocatable   :: regcoeff_psi_rh(:,:,:)     ! psi/rh     regression cooefficient.
   real, allocatable   :: regcoeff_chi_u_t(:,:,:)    ! chi_u/t    regression coefficient
   real, allocatable   :: regcoeff_chi_u_ps(:,:)     ! chi_u/ps   regression coefficient
   real, allocatable   :: regcoeff_chi_u_rh(:,:,:)   ! chi_u/rh   regression coefficient
   real, allocatable   :: regcoeff_t_u_rh(:,:,:)     ! t_u/rh     regression coefficient
   real, allocatable   :: regcoeff_ps_u_rh(:,:)      ! ps_u/rh    regression coefficient
   
   character(len=filename_len)        :: filename                   ! Input filename.
   character*3         :: ce                         ! Member index -> character.
   integer             :: ni, nj, nk                 ! Grid dimensions.
   integer             :: i, j, k, member            ! Loop counters.
   integer             :: b                          ! Bin marker.
   integer             :: sdate, cdate, edate        ! Starting, current ending dates.
   integer             :: interval                   ! Period between dates (hours).
   integer             :: ne                         ! Number of ensemble members.
   integer             :: bin_type                   ! Type of bin to average over.
   integer             :: num_bins                   ! Number of bins (3D fields).
   integer             :: num_bins2d                 ! Number of bins (3D fields).
   real                :: lat_min, lat_max           ! Used if bin_type = 2 (degrees).
   real                :: binwidth_lat               ! Used if bin_type = 2 (degrees).
   real                :: hgt_min, hgt_max           ! Used if bin_type = 2 (m).
   real                :: binwidth_hgt               ! Used if bin_type = 2 (m).
   real                :: coeffa, coeffb             ! Accumulating mean coefficients.
   logical             :: first_time                 ! True if first file.

   real, allocatable   :: field1(:,:,:)              ! Field 1.
   real, allocatable   :: field2(:,:,:)              ! Field 2 (full)
   real, allocatable   :: field2_balanced(:,:,:)   ! Field 2 (unbalanced)
   real*4, allocatable :: field(:,:) 
   integer, allocatable:: bin(:,:,:)                 ! Bin assigned to each 3D point.
   integer, allocatable:: bin2d(:,:)                 ! Bin assigned to each 2D point.
   integer, allocatable:: bin_pts(:)                 ! Number of points in bin (3D fields).
   real, allocatable   :: covar(:)                   ! Covariance between input fields.
   real, allocatable   :: var(:)                     ! Autocovariance of field.

   real, allocatable   :: regcoeff_field1_field2(:)
   namelist / gen_be_cov3d_nl / start_date, end_date, interval, &
                                ne, variable1, variable2

   integer :: ounit,iunit,namelist_unit

   stderr = 0
   stdout = 6

!---------------------------------------------------------------------------------------------
   write(6,'(a)')' [1] Initialize namelist variables and other scalars.'
!---------------------------------------------------------------------------------------------

   call da_get_unit(ounit)
   call da_get_unit(iunit)
   call da_get_unit(namelist_unit)


   start_date = '2007070100'
   end_date = '2007070500'
   interval = 24
   ne = 1
   variable1 = 'psi'
   variable2 = 't'

   open(unit=namelist_unit, file='gen_be_cov_contrib_nl.nl', &
        form='formatted', status='old', action='read')
   read(namelist_unit, gen_be_cov3d_nl)
   close(namelist_unit)

   read(start_date(1:10), fmt='(i10)')sdate
   read(end_date(1:10), fmt='(i10)')edate
   write(6,'(4a)')' Computing covariance for fields ', variable1 , ' and ', variable2
   write(6,'(4a)') ' Time period is ', start_date, ' to ', end_date
   write(6,'(a,i8,a)')' Interval between dates = ', interval, 'hours.'
   write(6,'(a,i8)')' Number of ensemble members at each time = ', ne

   date = start_date
   cdate = sdate

   member = 1
   write(ce,'(i3.3)')member

!  Read grid dimensions from first field file
   filename = trim(variable1)//'/'//date(1:10)
   filename = trim(filename)//'.'//trim(variable1)//'.e'//ce
   open (iunit, file = filename, form='unformatted')
   read(iunit)ni, nj, nk
   close(iunit)

!  Allocate
   write(6,'(a,3i8)')'    i, j, k dimensions are ', ni, nj, nk
   allocate( bin(1:ni,1:nj,1:nk) )
   allocate( bin2d(1:ni,1:nj) )
   allocate( field1(1:ni,1:nj,1:nk) )
   allocate( field2(1:ni,1:nj,1:nk) )
   allocate( field2_balanced(1:ni,1:nj,1:nk) )
   allocate( field(1:ni,1:nj) )

!  Read bin info:
   filename = 'bin.data'
   open (iunit, file = filename, form='unformatted')
   read(iunit)bin_type
   read(iunit)lat_min, lat_max, binwidth_lat
   read(iunit)hgt_min, hgt_max, binwidth_hgt
   read(iunit)num_bins, num_bins2d
   read(iunit)bin(1:ni,1:nj,1:nk)
   read(iunit)bin2d(1:ni,1:nj)
   close(iunit)

   allocate  (regcoeff_psi_chi(1:num_bins))
   allocate  (regcoeff_psi_t(1:nk,1:nk,1:num_bins2d))
   allocate  (regcoeff_psi_ps(1:nk,1:num_bins2d))
   allocate  (regcoeff_psi_rh(1:nk,1:nk,1:num_bins2d))
   allocate  (regcoeff_chi_u_t(1:nk,1:nk,1:num_bins2d))
   allocate  (regcoeff_chi_u_ps(1:nk,1:num_bins2d))
   allocate  (regcoeff_chi_u_rh(1:nk,1:nk,1:num_bins2d))
   allocate  (regcoeff_t_u_rh(1:nk,1:nk,1:num_bins2d))
   allocate  (regcoeff_ps_u_rh(1:nk,1:num_bins2d))

   allocate( regcoeff_field1_field2(1:num_bins) )

   allocate( bin_pts(1:num_bins) )
   allocate( covar(1:num_bins) )
   allocate( var(1:num_bins) )
   bin_pts(:) = 0
   covar(:) = 0.0
   var(:) = 0.0

!---------------------------------------------------------------------------------------------
   write(6,'(2a)')' [2] Read fields, and calculate correlations'
!--------------------------------------------------------------------------------------------- 
   filename = 'gen_be_stage2.dat'
   open(iunit, file = filename, form='unformatted',status='old')

   do k = 1, 2
   read (iunit)
   write(6,'(a)') 'rizvi read record ',k
   end do

   regcoeff_var='regcoeff_'//trim(variable1)//'_'//trim(variable2)
   do k = 1 , 9
   read (iunit) variable

   select case( trim(adjustl(variable)) )
   case ('regcoeff_psi_chi')
   read (iunit) regcoeff_psi_chi
   regcoeff_field1_field2 = regcoeff_psi_chi
   case ('regcoeff_psi_t')
   read (iunit) regcoeff_psi_t
   case ('regcoeff_psi_ps')
   read (iunit) regcoeff_psi_ps
   case ('regcoeff_psi_rh')
   read (iunit) regcoeff_psi_rh
   case ('regcoeff_chi_u_t')
   read (iunit) regcoeff_chi_u_t
   case ('regcoeff_chi_u_ps')
   read (iunit) regcoeff_chi_u_ps
   case ('regcoeff_chi_u_rh')
   read (iunit) regcoeff_chi_u_rh
   case ('regcoeff_t_u_rh')
   read (iunit) regcoeff_t_u_rh
   case ('regcoeff_ps_u_rh')
   read (iunit) regcoeff_ps_u_rh
   case default;
      write(6,fmt='(a)')'Read problem gen_be_stage2.dat'
      write(6,'(a,a)')' Trying to read regression coefficients:',trim(adjustl(variable))
      stop
   end select
   end do

   close(iunit)
   deallocate (regcoeff_psi_chi)
   deallocate (regcoeff_psi_t)
   deallocate (regcoeff_psi_ps)
   deallocate (regcoeff_psi_rh)
   deallocate (regcoeff_chi_u_t)
   deallocate (regcoeff_chi_u_ps)
   deallocate (regcoeff_chi_u_rh)
   deallocate (regcoeff_t_u_rh)
   deallocate (regcoeff_ps_u_rh)

   do while ( cdate <= edate )

      do member = 1, ne

         write(ce,'(i3.3)')member

!        Read first field:
         filename = trim(variable1)//'/'//date(1:10)
         filename = trim(filename)//'.'//trim(variable1)//'.e'//ce
         open (iunit, file = filename, form='unformatted')
         read(iunit)ni, nj, nk
         read(iunit)field1
         close(iunit)

!        Read second field:
         filename = trim(variable2)//'/'//date(1:10)
         filename = trim(filename)//'.'//trim(variable2)//'.e'//ce
         open (iunit, file = filename, form='unformatted')
         read(iunit)ni, nj, nk
         read(iunit)field2
         close(iunit)

!        Calculate balanced contribution to field2
         do k = 1, nk
            do j = 1, nj
               do i = 1, ni
                  b = bin(i,j,k)
                  field2_balanced(i,j,k) = regcoeff_field1_field2(b) * field1(i,j,k)
               end do
            end do
         end do
         
!        Calculate covariances:
         do k = 1, nk
            do j = 1, nj
               do i = 1, ni
                  b = bin(i,j,k)
                  bin_pts(b) = bin_pts(b) + 1
                  coeffa = 1.0 / real(bin_pts(b))
                  coeffb = real(bin_pts(b)-1) * coeffa
                  covar(b) = coeffb * covar(b) + coeffa * field2_balanced(i,j,k) * field2(i,j,k)
                  var(b) = coeffb * var(b) + coeffa * field2(i,j,k) * field2(i,j,k)
               end do
            end do
         end do

      end do  ! End loop over ensemble members.

!     Calculate next date:
      call da_advance_cymdh( date, interval, new_date )
      date = new_date
      read(date(1:10), fmt='(i10)')cdate

   end do     ! End loop over times.

   filename = trim(variable1)//'.'//trim(variable2)//'.bin'
   open (ounit, file = filename, form = 'unformatted')
   do k = 1, nk
      do j = 1, nj
         do i = 1, ni
            b = bin(i,j,k) ! Take value from center of i dimension.
            if ( var(b) /= 0.0 ) then
               field(i,j) = covar(b) / var(b)
            else
               field(i,j) = 0.
            end if
         end do
      end do
      write(ounit) field
   end do
   close(ounit)


end program gen_be_cov3d3d_bin3d_contrib

