module hybrid_ensemble_isotropic_global
!$$$   module documentation block
!                .      .    .                                       .
! module:    hybrid_ensemble_isotropic_global
!   prgmmr: parrish          org: np22                date: 2009-10-15
!
! abstract: contains routines for localization of the hybrid ensemble
!            control variable a_en.  this application is for
!            localization with a global model with homogeneous scale in horizontal
!            and vertical scale a function of vertical only.
!
! program history log:
!   2009-10-15  parrish, initial documentation.
!   2010-02-20  parrish, remove modules specmod_special, transform_special, spectral_transforms_special,
!                 as part of streamlining code for dual resolution capability.
!
! subroutines included:
!   sub init_sf_xy               - initialize spectral localization routine sf_xy
!   sub sf_xy                    - spectral localization routine
!
! Variable Definitions:
!   def spectral_filter          - used to guarantee proper form for zonal spectral coefficients
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$ end documentation block

  use kinds, only: r_kind,i_kind
  implicit none

! set default to private
  private
! set subroutines to public
  public :: init_sf_xy
  public :: sf_xy

  real(r_kind),allocatable,dimension(:)::spectral_filter

contains

subroutine init_sf_xy(jcap_in)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    init_sf_xy
!   prgmmr: parrish          org: np23                date: 2009-12-17
!
! abstract: initialize horizontal spectral localization sf_xy
!
! program history log:
!   2009-12-17  parrish
!
!   input argument list:
!     jcap_in - maximum spectral truncation allowed
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!$$$

  use kinds, only: r_kind,i_kind,r_single
  use hybrid_ensemble_parameters,only: s_ens_h,sp_loc,grd_ens
  use general_specmod, only: general_init_spec_vars 
  use constants, only: izero,ione,zero,half,one,two,three,rearth,pi
  use mpimod, only: mype
  implicit none

  integer(i_kind),intent(in   ) :: jcap_in

  integer(i_kind) i,ii,j,l,n,nn,jcap
  real(r_kind),allocatable::g(:),gsave(:),errmax(:)
  real(r_kind) factor
  real(r_kind) rkm(grd_ens%nlat),f(grd_ens%nlat,grd_ens%nlon),f0(grd_ens%nlat,grd_ens%nlon)
  real(r_single) out1(grd_ens%nlon,grd_ens%nlat)
  real(r_single),allocatable::pn0_npole(:)
  real(r_kind) s_ens_h_min
  

!    make sure s_ens_h is within allowable range  ( pi*rearth*.001/jcap_in <= s_ens_h <= 5500 )

  s_ens_h_min=pi*rearth*.001_r_kind/jcap_in
  if(s_ens_h <  s_ens_h_min) then
     if(mype == izero) write(6,*)' s_ens_h = ',s_ens_h,' km--too small, min value = ',s_ens_h_min,' km.'
     if(mype == izero) write(6,*)' s_ens_h reset to min value'
     s_ens_h=s_ens_h_min
  else if(s_ens_h >  5500._r_kind) then
     if(mype == izero) write(6,*)' s_ens_h = ',s_ens_h,' km--too large, max value = 5500 km.'
     if(mype == izero) write(6,*)' s_ens_h reset to max value'
     s_ens_h=5500._r_kind
  end if
    
    
  jcap=nint(1.2_r_kind*pi*rearth*.001_r_kind/s_ens_h)
  jcap=min(jcap,jcap_in)

!  set up spectral variables for jcap

  call general_init_spec_vars(sp_loc,jcap,jcap,grd_ens%nlat,grd_ens%nlon)

!    the following code is used to compute the desired spectrum to get a 
!     gaussian localization of desired length-scale.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!   how spectrum is obtained:
!
!     Correlation matrix:  C = Y*D*Ytrans
!
!       where Y is matrix of spherical harmonics, evaluated on gaussian grid, and D is a diagonal matrix
!
!     To obtain D, exploit fact that for D a function only of total wave-number n, then C is homogeneous
!       and isotropic on the sphere.
!
!     So look at the special case of a test point centered on the north pole.  The correlation function
!       is then only a function of latitude, call it c(phi), where c(pi/2) = 1.
!
!     Now we have C = P*D*Ptrans, where we have reduced the problem to 1 dimension, latitude, and in 
!       spectral space, total wave number n.  P is the zonal component only of Y.
!
!     Next, form the product
!                             C*e1 =P*D*Ptrans*e1,
!
!            where e1 is a vector of all 0, except for 1 at the north pole.
!
!     Then have P*D*Ptrans*e1 = sum(n) p(n,j)*d(n)*p(n,1) = c(j,1)
!
!        where j=1 corresponds to north pole point in this formulation.
!
!     Now if we have available C(j,1), a gaussian of desired length-scale evaluated (note, doesn't have to
!            be gaussian!) on the gaussian grid, then applying the inverse transform subroutine g2s0 to 
!     C yields the product
!
!             Chat(n) = d(n)*p(n,1)
!
!     So finally the desired spectrum is 
!     
!               d(n) = chat(n)/p(n,1)
!
!     To create the full spectral transform of the desired correlation, d(n) is copied to all non-zero 
!      zonal wave numbers on lines of constant total wave number n, multiplied by 1/2 to account for
!      two degrees of freedom for non-zero zonal wave numbers.
!
!     Note that while creating this routine, a bug was discovered in s2g0 routine for evaluation of pole 
!       value.  There was a missing factor of 1/sqrt(2) apparently.
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!                 create reference gaussian centered on north pole with desired length-scale.

!     compute latitudes in km from north pole
  rkm(grd_ens%nlat)=zero
  rkm(1)=two*asin(one)*rearth*.001_r_kind
  do i=1,(grd_ens%nlat-2_i_kind)/2
     rkm(grd_ens%nlat-i)=(asin(one)-asin(sp_loc%slat(i)))*rearth*.001_r_kind
     rkm(ione+i)=(asin(one)+asin(sp_loc%slat(i)))*rearth*.001_r_kind
  end do
  if(mype == izero) write(6,*)' in init_sf_xy, lat,max(dlat)=', &
           rkm(ione+(grd_ens%nlat-2_i_kind)/2), &
          -rkm(grd_ens%nlat-(grd_ens%nlat-2_i_kind)/2)+rkm(ione+(grd_ens%nlat-2_i_kind)/2),' km'

  do i=1,grd_ens%nlat
     f0(i,1)=exp(-half*(rkm(i)/s_ens_h)**2)
  end do
  
  do j=2,grd_ens%nlon
     do i=1,grd_ens%nlat
        f0(i,j)=f0(i,1)
     end do
  end do

  allocate(g(sp_loc%nc),gsave(sp_loc%nc))
  call general_g2s0(grd_ens,sp_loc,g,f0)

  call general_s2g0(grd_ens,sp_loc,g,f)

!    adjust so value at np = 1
  f=f/f(grd_ens%nlat,1)
  call general_g2s0(grd_ens,sp_loc,g,f)
  call general_s2g0(grd_ens,sp_loc,g,f)
  if(mype == izero) write(6,*)' in init_sf_xy, jcap,s_ens_h,max diff(f0-f)=', &
                                            sp_loc%jcap,s_ens_h,maxval(abs(f0-f))

!            correct spectrum by dividing by pn0_npole
  gsave=g

!    obtain pn0_npole
  allocate(pn0_npole(0:sp_loc%jcap))
  do n=0,sp_loc%jcap
     g=zero
     g(2*n+ione)=one
     call general_s2g0(grd_ens,sp_loc,g,f)
     pn0_npole(n)=f(grd_ens%nlat,1)
  end do

  g=zero
  do n=0,sp_loc%jcap
     g(2*n+ione)=gsave(2*n+ione)/pn0_npole(n)
  end do

!    obtain spectral_filter

  allocate(spectral_filter(sp_loc%nc))
  ii=izero
  do l=0,sp_loc%jcap
     factor=one
     if(l >  izero) factor=half
     do n=l,sp_loc%jcap
        nn=nn+2_i_kind
        ii=ii+ione
        spectral_filter(ii)=factor*g(2*n+ione)
        ii=ii+ione
        if(l == izero) then
           spectral_filter(ii)=zero
        else
           spectral_filter(ii)=factor*g(2*n+ione)
        end if
     end do
  end do

  deallocate(g,gsave,pn0_npole)

  f=zero
  f(grd_ens%nlat/4,grd_ens%nlon/2)=one
  call sf_xy(f,1)
  do j=1,grd_ens%nlon
     do i=1,grd_ens%nlat
        out1(j,i)=f(i,j)
     end do
  end do
  if(mype == izero) call outgrads1(out1,grd_ens%nlon,grd_ens%nlat,'out1')

end subroutine init_sf_xy

subroutine sf_xy(f,nlevs)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    sf_xy       spectral isotropic localization for global domain
!
!   prgrmmr:  parrish        org: np22                              date: 2009-09-30
!
! abstract:  use spherical harmonic transform to do horizontal isotropic localization of hybrid
!             ensemble control variable a_en.
!
! program history log:
!   2009-10-16  parrish  initial documentation
!   2010-03-11  parrish - adjust dimensions for f to allow for nlevs=0
!
!   input argument list:
!     f        - input field to be filtered
!     nlevs    -
!
!   output argument list:
!     f        - filtered output
!
! attributes:
!   language:  f90
!   machine:   ibm RS/6000 SP
!
!$$$ end documentation block

  use kinds, only: r_kind,i_kind
  use hybrid_ensemble_parameters, only: grd_ens,sp_loc
  implicit none

  integer(i_kind),intent(in   ) :: nlevs
  real(r_kind)   ,intent(inout) :: f(grd_ens%nlat*grd_ens%nlon,max(nlevs,1))

  real(r_kind) g(sp_loc%nc)
  integer(i_kind) k

  do k=1,nlevs
     call general_s2g0_ad(grd_ens,sp_loc,g,f(:,k))
     g=g*spectral_filter
     call general_s2g0(grd_ens,sp_loc,g,f(:,k))
  end do

end subroutine sf_xy


end module hybrid_ensemble_isotropic_global