type_samp.F90

Go to the documentation of this file.

!----------------------------------------------------------------------
! Module: type_samp
! Purpose: sampling derived type
! Author: Benjamin Menetrier
! Licensing: this code is distributed under the CeCILL-C license
! Copyright © 2015-... UCAR, CERFACS, METEO-FRANCE and IRIT
!----------------------------------------------------------------------
module type_samp

use fckit_mpi_module, only: fckit_mpi_sum,fckit_mpi_status
use netcdf
!$ use omp_lib
use tools_const, only: pi,req,deg2rad,rad2deg,msvali,msvalr
use tools_func, only: gc99,sphere_dist,vector_product,vector_triple_product
use tools_kinds, only: kind_real
use tools_missing, only: msi,msr,isnotmsi,isnotmsr,ismsi,ismsr
use tools_nc, only: ncfloat
use tools_qsort, only: qsort
use tools_stripack, only: trans
use type_bpar, only: bpar_type
use type_com, only: com_type
use type_ens, only: ens_type
use type_geom, only: geom_type
use type_io, only: io_type
use type_kdtree, only: kdtree_type
use type_linop, only: linop_type
use type_mesh, only: mesh_type
use type_mpl, only: mpl_type
use type_nam, only: nam_type
use type_rng, only: rng_type

implicit none

integer,parameter :: irmax = 10000                           ! Maximum number of random number draws
real(kind_real),parameter :: lcoast = 1000.0e3_kind_real/req ! Length-scale to increase sampling density along coasts
real(kind_real),parameter :: rcoast = 0.2_kind_real          ! Minimum value to increase sampling density along coasts

! Sampling derived type
type samp_type
   ! Sampling
   real(kind_real),allocatable :: rh_c0(:,:)        ! Sampling radius on subset Sc0
   logical,allocatable :: mask_c0(:,:)              ! Mask on subset Sc0
   logical,allocatable :: mask_hor_c0(:)            ! Union of horizontal masks on subset Sc0, global
   integer,allocatable :: nc0_mask(:)               ! Horizontal mask size on subset Sc0
   integer,allocatable :: c1_to_c0(:)               ! First sampling index
   logical,allocatable :: c1l0_log(:,:)             ! Log for the first sampling index
   integer,allocatable :: c1c3_to_c0(:,:)           ! Second horizontal sampling index
   logical,allocatable :: c1c3l0_log(:,:,:)         ! Log for the second horizontal sampling index
   integer,allocatable :: c2_to_c1(:)               ! Subgrid to diagnostic points
   integer,allocatable :: c2_to_c0(:)               ! Subgrid to grid

   ! Local data
   logical,allocatable ::  vbal_mask(:,:)           ! Vertical balance mask
   logical,allocatable ::  local_mask(:,:)          ! Local mask
   logical,allocatable ::  displ_mask(:,:)          ! Displacement mask
   integer,allocatable :: nn_c2_index(:,:,:)        ! Nearest diagnostic neighbors from diagnostic points
   real(kind_real),allocatable :: nn_c2_dist(:,:,:) ! Nearest diagnostic neighbors distance from diagnostic points
   integer,allocatable :: nn_ldwv_index(:)          ! Nearest diagnostic neighbors for local diagnostics profiles

   ! Sampling mesh
   type(mesh_type) :: mesh                          ! Sampling mesh

   ! Displacement
   real(kind_real),allocatable :: displ_lon(:,:,:)  ! Interpolated displaced longitude
   real(kind_real),allocatable :: displ_lat(:,:,:)  ! Interpolated displaced latitude

   ! Interpolations
   type(linop_type),allocatable :: hfull(:)         ! Horizontal interpolation from Sc2 to Sc0 (global)
   type(linop_type),allocatable :: h(:)             ! Horizontal interpolation from Sc2 to Sc0 (local)
   type(linop_type),allocatable :: d(:,:)           ! Displacement interpolation

   ! MPI distribution
   integer :: nc0c                                  ! Number of points in subset Sc0, halo C
   integer :: nc0d                                  ! Number of points in subset Sc0, halo D
   integer :: nc1a                                  ! Number of points in subset Sc1, halo A
   integer :: nc2a                                  ! Number of points in subset Sc2, halo A
   integer :: nc2b                                  ! Number of points in subset Sc2, halo B
   integer :: nc2f                                  ! Number of points in subset Sc2, halo F
   logical,allocatable :: lcheck_c0a(:)             ! Detection of halo A on subset Sc0
   logical,allocatable :: lcheck_c0c(:)             ! Detection of halo C on subset Sc0
   logical,allocatable :: lcheck_c0d(:)             ! Detection of halo D on subset Sc0
   logical,allocatable :: lcheck_c1a(:)             ! Detection of halo A on subset Sc1
   logical,allocatable :: lcheck_c2a(:)             ! Detection of halo A on subset Sc2
   logical,allocatable :: lcheck_c2b(:)             ! Detection of halo B on subset Sc2
   logical,allocatable :: lcheck_c2f(:)             ! Detection of halo F on subset Sc2
   logical,allocatable :: lcheck_h(:,:)             ! Detection of horizontal interpolation coefficients
   logical,allocatable :: lcheck_d(:,:,:)           ! Detection of displacement interpolation coefficients
   integer,allocatable :: c0c_to_c0(:)              ! Subset Sc0, halo C to global
   integer,allocatable :: c0_to_c0c(:)              ! Subset Sc0, global to halo C
   integer,allocatable :: c0a_to_c0c(:)             ! Subset Sc0, halo A to halo C
   integer,allocatable :: c0d_to_c0(:)              ! Subset Sc0, halo D to global
   integer,allocatable :: c0_to_c0d(:)              ! Subset Sc0, global to halo D
   integer,allocatable :: c0a_to_c0d(:)             ! Subset Sc0, halo A to halo D
   integer,allocatable :: c1a_to_c1(:)              ! Subset Sc1, halo A to global
   integer,allocatable :: c1_to_c1a(:)              ! Subset Sc1, global to halo A
   integer,allocatable :: c2a_to_c2(:)              ! Subset Sc2, halo A to global
   integer,allocatable :: c2_to_c2a(:)              ! Subset Sc2, global to halo A
   integer,allocatable :: c2b_to_c2(:)              ! Subset Sc2, halo B to global
   integer,allocatable :: c2_to_c2b(:)              ! Subset Sc2, global to halo B
   integer,allocatable :: c2a_to_c2b(:)             ! Subset Sc2, halo A to halo B
   integer,allocatable :: c2f_to_c2(:)              ! Subset Sc2, halo B to global
   integer,allocatable :: c2_to_c2f(:)              ! Subset Sc2, global to halo B
   integer,allocatable :: c2a_to_c2f(:)             ! Subset Sc2, halo A to halo B
   integer,allocatable :: c2_to_proc(:)             ! Subset Sc2, global to processor
   integer,allocatable :: proc_to_nc2a(:)           ! Number of points in subset Sc2, halo A, for each processor
   type(com_type) :: com_ac                         ! Communication between halos A and C
   type(com_type) :: com_ab                         ! Communication between halos A and B
   type(com_type) :: com_ad                         ! Communication between halos A and D
   type(com_type) :: com_af                         ! Communication between halos A and F (filtering)
contains
   procedure :: alloc => samp_alloc
   procedure :: dealloc => samp_dealloc
   procedure :: read => samp_read
   procedure :: write => samp_write
   procedure :: setup_sampling => samp_setup_sampling
   procedure :: compute_mask => samp_compute_mask
   procedure :: compute_sampling_zs => samp_compute_sampling_zs
   procedure :: compute_sampling_ps => samp_compute_sampling_ps
   procedure :: compute_sampling_lct => samp_compute_sampling_lct
   procedure :: compute_sampling_mask => samp_compute_sampling_mask
   procedure :: compute_mpi_a => samp_compute_mpi_a
   procedure :: compute_mpi_ab => samp_compute_mpi_ab
   procedure :: compute_mpi_d => samp_compute_mpi_d
   procedure :: compute_mpi_c => samp_compute_mpi_c
   procedure :: compute_mpi_f => samp_compute_mpi_f
   procedure :: diag_filter => samp_diag_filter
   procedure :: diag_fill => samp_diag_fill
end type samp_type

private
public :: samp_type

contains

!----------------------------------------------------------------------
! Subroutine: samp_alloc
! Purpose: sampling allocation
!----------------------------------------------------------------------
subroutine samp_alloc(samp,nam,geom)

implicit none

! Passed variables
class(samp_type),intent(inout) :: samp ! Sampling
type(nam_type),intent(in) :: nam       ! Namelist
type(geom_type),intent(in) :: geom     ! Geometry

! Allocation
allocate(samp%rh_c0(geom%nc0,geom%nl0))
allocate(samp%mask_c0(geom%nc0,geom%nl0))
allocate(samp%mask_hor_c0(geom%nc0))
allocate(samp%nc0_mask(geom%nl0))
allocate(samp%c1_to_c0(nam%nc1))
allocate(samp%c1l0_log(nam%nc1,geom%nl0))
allocate(samp%c1c3_to_c0(nam%nc1,nam%nc3))
allocate(samp%c1c3l0_log(nam%nc1,nam%nc3,geom%nl0))
if (nam%new_vbal.or.nam%new_lct.or.(nam%new_hdiag.and.(nam%var_diag.or.nam%local_diag.or.nam%displ_diag))) then
   allocate(samp%c2_to_c1(nam%nc2))
   allocate(samp%c2_to_c0(nam%nc2))
   allocate(samp%nn_c2_index(nam%nc2,nam%nc2,geom%nl0i))
   allocate(samp%nn_c2_dist(nam%nc2,nam%nc2,geom%nl0i))
   allocate(samp%hfull(geom%nl0i))
end if
if (nam%new_vbal) allocate(samp%vbal_mask(nam%nc1,nam%nc2))
if (nam%local_diag) allocate(samp%local_mask(nam%nc1,nam%nc2))
if (nam%displ_diag) then
   allocate(samp%displ_mask(nam%nc1,nam%nc2))
   allocate(samp%displ_lon(geom%nc0a,geom%nl0,nam%nts))
   allocate(samp%displ_lat(geom%nc0a,geom%nl0,nam%nts))
end if

! Initialization
call msr(samp%rh_c0)
samp%mask_c0 = .false.
samp%mask_hor_c0 = .false.
call msi(samp%nc0_mask)
call msi(samp%c1_to_c0)
samp%c1l0_log = .false.
call msi(samp%c1c3_to_c0)
samp%c1c3l0_log = .false.
if (nam%new_vbal.or.nam%new_lct.or.(nam%new_hdiag.and.(nam%var_diag.or.nam%local_diag.or.nam%displ_diag))) then
   call msi(samp%c2_to_c1)
   call msi(samp%c2_to_c0)
   call msi(samp%nn_c2_index)
   call msr(samp%nn_c2_dist)
end if
if (nam%new_vbal) samp%vbal_mask = .false.
if (nam%local_diag) samp%local_mask = .false.
if (nam%displ_diag) then
   samp%displ_mask = .false.
   call msr(samp%displ_lon)
   call msr(samp%displ_lat)
end if

end subroutine samp_alloc

!----------------------------------------------------------------------
! Subroutine: samp_dealloc
! Purpose: sampling deallocation
!----------------------------------------------------------------------
subroutine samp_dealloc(samp,geom)

implicit none

! Passed variables
class(samp_type),intent(inout) :: samp ! Sampling
type(geom_type),intent(in) :: geom     ! Geometry

! Local variables
integer :: il0

! Release memory
if (allocated(samp%rh_c0)) deallocate(samp%rh_c0)
if (allocated(samp%mask_c0)) deallocate(samp%mask_c0)
if (allocated(samp%mask_hor_c0)) deallocate(samp%mask_hor_c0)
if (allocated(samp%nc0_mask)) deallocate(samp%nc0_mask)
if (allocated(samp%c1_to_c0)) deallocate(samp%c1_to_c0)
if (allocated(samp%c1l0_log)) deallocate(samp%c1l0_log)
if (allocated(samp%c1c3_to_c0)) deallocate(samp%c1c3_to_c0)
if (allocated(samp%c1c3l0_log)) deallocate(samp%c1c3l0_log)
if (allocated(samp%c2_to_c1)) deallocate(samp%c2_to_c1)
if (allocated(samp%c2_to_c0)) deallocate(samp%c2_to_c0)
if (allocated(samp%c2a_to_c2)) deallocate(samp%c2a_to_c2)
if (allocated(samp%c2_to_c2a)) deallocate(samp%c2_to_c2a)
if (allocated(samp%vbal_mask)) deallocate(samp%vbal_mask)
if (allocated(samp%local_mask)) deallocate(samp%local_mask)
if (allocated(samp%displ_mask)) deallocate(samp%displ_mask)
if (allocated(samp%nn_c2_index)) deallocate(samp%nn_c2_index)
if (allocated(samp%nn_c2_dist)) deallocate(samp%nn_c2_dist)
if (allocated(samp%hfull)) then
   do il0=1,geom%nl0
      call samp%hfull(il0)%dealloc
   end do
   deallocate(samp%hfull)
end if
if (allocated(samp%h)) then
   do il0=1,geom%nl0
      call samp%h(il0)%dealloc
   end do
   deallocate(samp%h)
end if
if (allocated(samp%nn_ldwv_index)) deallocate(samp%nn_ldwv_index)
if (allocated(samp%displ_lon)) deallocate(samp%displ_lon)
if (allocated(samp%displ_lat)) deallocate(samp%displ_lat)

end subroutine samp_dealloc

!----------------------------------------------------------------------
! Subroutine: samp_read
! Purpose: read sampling
!----------------------------------------------------------------------
subroutine samp_read(samp,mpl,nam,geom,bpar,ios)

implicit none

! Passed variables
class(samp_type),intent(inout) :: samp ! Sampling
type(mpl_type),intent(in) :: mpl       ! MPI data
type(nam_type),intent(inout) :: nam    ! Namelist
type(geom_type),intent(in) :: geom     ! Geometry
type(bpar_type),intent(in) :: bpar     ! Block parameters
integer,intent(out) :: ios             ! Status flag

! Local variables
integer :: il0,il0i,ic1,jc3,ic2
integer :: nl0_test,nl0r_test,nc_test,nc1_test,nc2_test,nc2_1_test,nc2_2_test
integer :: info,ncid,nl0_id,nc3_id,nc1_id,nc2_id,nc2_1_id,nc2_2_id
integer :: c1_to_c0_id,c1l0_log_id,c1c3_to_c0_id,c1c3l0_log_id
integer :: c2_to_c1_id,c2_to_c0_id,vbal_mask_id,local_mask_id,displ_mask_id,nn_c2_index_id,nn_c2_dist_id
integer :: c1l0_logint(nam%nc1,geom%nl0),c1c3l0_logint(nam%nc1,nam%nc3,geom%nl0)
integer,allocatable :: vbal_maskint(:,:),local_maskint(:,:),displ_maskint(:,:)
character(len=3) :: il0ichar
character(len=1024) :: subr = 'samp_read'

! Initialization
ios = 0

! Open file
info = nf90_open(trim(nam%datadir)//'/'//trim(nam%prefix)//'_sampling.nc',nf90_nowrite,ncid)
if (info/=nf90_noerr) then
   call mpl%warning('cannot find sampling to read, recomputing sampling')
   nam%sam_write = .true.
   ios = 1
   return
end if

! Check dimensions
call mpl%ncerr(subr,nf90_inq_dimid(ncid,'nl0',nl0_id))
call mpl%ncerr(subr,nf90_inquire_dimension(ncid,nl0_id,len=nl0_test))
call mpl%ncerr(subr,nf90_get_att(ncid,nf90_global,'nl0r',nl0r_test))
call mpl%ncerr(subr,nf90_inq_dimid(ncid,'nc3',nc3_id))
call mpl%ncerr(subr,nf90_inquire_dimension(ncid,nc3_id,len=nc_test))
call mpl%ncerr(subr,nf90_inq_dimid(ncid,'nc1',nc1_id))
call mpl%ncerr(subr,nf90_inquire_dimension(ncid,nc1_id,len=nc1_test))
if (nam%new_lct.or.nam%var_diag.or.nam%local_diag.or.nam%displ_diag) then
   info = nf90_inq_dimid(ncid,'nc2',nc2_id)
   if (info==nf90_noerr) then
      call mpl%ncerr(subr,nf90_inquire_dimension(ncid,nc2_id,len=nc2_test))
   else
      call mpl%warning('cannot find nc2 when reading sampling, recomputing sampling')
      nam%sam_write = .true.
      ios = 2
   end if
end if
if ((geom%nl0/=nl0_test).or.(bpar%nl0rmax/=nl0r_test).or.(nam%nc3/=nc_test).or.(nam%nc1/=nc1_test)) then
   call mpl%warning('wrong dimension when reading sampling, recomputing sampling')
   nam%sam_write = .true.
   call mpl%ncerr(subr,nf90_close(ncid))
   ios = 1
   return
end if
if (nam%new_lct.or.nam%var_diag.or.nam%local_diag.or.nam%displ_diag) then
   if (nam%nc2/=nc2_test) then
      call mpl%warning('wrong dimension when reading sampling, recomputing sampling')
      nam%sam_write = .true.
      ios = 2
   end if
end if

write(mpl%info,'(a7,a)') '','Read sampling'
call flush(mpl%info)

! Get arrays ID
call mpl%ncerr(subr,nf90_inq_varid(ncid,'c1_to_c0',c1_to_c0_id))
call mpl%ncerr(subr,nf90_inq_varid(ncid,'c1l0_log',c1l0_log_id))
call mpl%ncerr(subr,nf90_inq_varid(ncid,'c1c3_to_c0',c1c3_to_c0_id))
call mpl%ncerr(subr,nf90_inq_varid(ncid,'c1c3l0_log',c1c3l0_log_id))
if ((ios==0).and.(nam%new_lct.or.nam%var_diag.or.nam%local_diag.or.nam%displ_diag)) then
   call mpl%ncerr(subr,nf90_inq_varid(ncid,'c2_to_c1',c2_to_c1_id))
   call mpl%ncerr(subr,nf90_inq_varid(ncid,'c2_to_c0',c2_to_c0_id))
end if

! Read arrays
call mpl%ncerr(subr,nf90_get_var(ncid,c1_to_c0_id,samp%c1_to_c0))
call mpl%ncerr(subr,nf90_get_var(ncid,c1l0_log_id,c1l0_logint))
do il0=1,geom%nl0
   do ic1=1,nam%nc1
      if (c1l0_logint(ic1,il0)==0) then
         samp%c1l0_log(ic1,il0) = .false.
      else if (c1l0_logint(ic1,il0)==1) then
         samp%c1l0_log(ic1,il0) = .true.
      end if
   end do
end do
call mpl%ncerr(subr,nf90_get_var(ncid,c1c3_to_c0_id,samp%c1c3_to_c0))
call mpl%ncerr(subr,nf90_get_var(ncid,c1c3l0_log_id,c1c3l0_logint))
do il0=1,geom%nl0
   do jc3=1,nam%nc3
      do ic1=1,nam%nc1
         if (c1c3l0_logint(ic1,jc3,il0)==0) then
            samp%c1c3l0_log(ic1,jc3,il0) = .false.
         else if (c1c3l0_logint(ic1,jc3,il0)==1) then
            samp%c1c3l0_log(ic1,jc3,il0) = .true.
         end if
      end do
   end do
end do
if ((ios==0).and.(nam%new_lct.or.nam%var_diag.or.nam%local_diag.or.nam%displ_diag)) then
   call mpl%ncerr(subr,nf90_get_var(ncid,c2_to_c1_id,samp%c2_to_c1))
   call mpl%ncerr(subr,nf90_get_var(ncid,c2_to_c0_id,samp%c2_to_c0))
end if

! Close file
call mpl%ncerr(subr,nf90_close(ncid))

! Read nearest neighbors and interpolation
do il0i=1,geom%nl0i
   if ((ios==0).and.(nam%new_vbal.or.nam%new_lct.or.(nam%new_hdiag.and.(nam%var_diag.or.nam%local_diag.or.nam%displ_diag)))) then
      ! Open file
      write(il0ichar,'(i3.3)') il0i
      info = nf90_open(trim(nam%datadir)//'/'//trim(nam%prefix)//'_sampling_'//il0ichar//'.nc',nf90_nowrite,ncid)
      if (info/=nf90_noerr) then
         call mpl%warning('cannot find nearest neighbors and interpolation data to read, recomputing sampling')
         nam%sam_write = .true.
         ios = 3
         return
      end if

      ! Check dimensions
      call mpl%ncerr(subr,nf90_inq_dimid(ncid,'nc1',nc1_id))
      call mpl%ncerr(subr,nf90_inquire_dimension(ncid,nc1_id,len=nc1_test))
      call mpl%ncerr(subr,nf90_inq_dimid(ncid,'nc2_1',nc2_1_id))
      call mpl%ncerr(subr,nf90_inquire_dimension(ncid,nc2_1_id,len=nc2_1_test))
      call mpl%ncerr(subr,nf90_inq_dimid(ncid,'nc2_2',nc2_2_id))
      call mpl%ncerr(subr,nf90_inquire_dimension(ncid,nc2_2_id,len=nc2_2_test))
      if ((nam%nc1/=nc1_test).or.(nam%nc2/=nc2_1_test).or.(nam%nc2/=nc2_2_test)) then
         call mpl%warning('wrong dimension when reading sampling, recomputing sampling')
         nam%sam_write = .true.
         ios = 2
      end if
      info = nf90_inq_varid(ncid,'nn_c2_index',nn_c2_index_id)
      if (info==nf90_noerr) then
         call mpl%ncerr(subr,nf90_inq_varid(ncid,'nn_c2_dist',nn_c2_dist_id))
         call mpl%ncerr(subr,nf90_get_var(ncid,nn_c2_index_id,samp%nn_c2_index(:,:,il0i)))
         call mpl%ncerr(subr,nf90_get_var(ncid,nn_c2_dist_id,samp%nn_c2_dist(:,:,il0i)))
      else
         call mpl%warning('cannot find nc2 nearest neighbors data to read, recomputing sampling')
         nam%sam_write = .true.
         ios = 4
      end if
      write(samp%hfull(il0i)%prefix,'(a,i3.3)') 'hfull_',il0i
      call samp%hfull(il0i)%read(mpl,ncid)
   end if

   if ((ios==0).and.nam%new_vbal) then
      ! Allocation
      allocate(vbal_maskint(nam%nc1,nam%nc2))

      call mpl%ncerr(subr,nf90_inq_varid(ncid,'vbal_mask',vbal_mask_id))
      call mpl%ncerr(subr,nf90_get_var(ncid,vbal_mask_id,vbal_maskint))
      do ic2=1,nam%nc2
         do ic1=1,nam%nc1
            if (vbal_maskint(ic1,ic2)==1) then
               samp%vbal_mask(ic1,ic2) = .true.
            elseif (vbal_maskint(ic1,ic2)==0) then
               samp%vbal_mask(ic1,ic2) = .false.
            else
               call mpl%abort('wrong vbal_mask')
           end if
         end do
      end do

      ! Release memory
      deallocate(vbal_maskint)
   end if

   if ((ios==0).and.nam%local_diag) then
      ! Allocation
      allocate(local_maskint(nam%nc1,nam%nc2))

      call mpl%ncerr(subr,nf90_inq_varid(ncid,'local_mask',local_mask_id))
      call mpl%ncerr(subr,nf90_get_var(ncid,local_mask_id,local_maskint))
      do ic2=1,nam%nc2
         do ic1=1,nam%nc1
            if (local_maskint(ic1,ic2)==1) then
               samp%local_mask(ic1,ic2) = .true.
            elseif (local_maskint(ic1,ic2)==0) then
               samp%local_mask(ic1,ic2) = .false.
            else
               call mpl%abort('wrong local_mask')
           end if
         end do
      end do

      ! Release memory
      deallocate(local_maskint)
   end if

   if ((ios==0).and.nam%displ_diag) then
      ! Allocation
      allocate(displ_maskint(nam%nc1,nam%nc2))

      call mpl%ncerr(subr,nf90_inq_varid(ncid,'displ_mask',displ_mask_id))
      call mpl%ncerr(subr,nf90_get_var(ncid,displ_mask_id,displ_maskint))
      do ic2=1,nam%nc2
         do ic1=1,nam%nc1
            if (displ_maskint(ic1,ic2)==1) then
               samp%displ_mask(ic1,ic2) = .true.
            elseif (displ_maskint(ic1,ic2)==0) then
               samp%displ_mask(ic1,ic2) = .false.
            else
                call mpl%abort('wrong displ_mask')
            end if
         end do
      end do

      ! Release memory
      deallocate(displ_maskint)
   end if

   if ((ios==0).and.(nam%new_lct.or.nam%var_diag.or.nam%local_diag.or.nam%displ_diag)) then
      ! Close file
      call mpl%ncerr(subr,nf90_close(ncid))
   end if
end do

end subroutine samp_read

!----------------------------------------------------------------------
! Subroutine: samp_write
! Purpose: write sampling
!----------------------------------------------------------------------
subroutine samp_write(samp,mpl,nam,geom,bpar)

implicit none

! Passed variables
class(samp_type),intent(in) :: samp ! Sampling
type(mpl_type),intent(in) :: mpl    ! MPI data
type(nam_type),intent(in) :: nam    ! Namelist
type(geom_type),intent(in) :: geom  ! Geometry
type(bpar_type),intent(in) :: bpar  ! Block parameters

! Local variables
integer :: il0,il0i,ic1,jc3,ic2
integer :: ncid,nl0_id,nc1_id,nc2_id,nc2_1_id,nc2_2_id,nc3_id
integer :: lat_id,lon_id,smax_id,c1_to_c0_id,c1l0_log_id,c1c3_to_c0_id,c1c3l0_log_id
integer :: c2_to_c1_id,c2_to_c0_id,vbal_mask_id,local_mask_id,displ_mask_id,nn_c2_index_id,nn_c2_dist_id
integer :: c1l0_logint(nam%nc1,geom%nl0),c1c3l0_logint(nam%nc1,nam%nc3,geom%nl0)
integer,allocatable :: vbal_maskint(:,:),local_maskint(:,:),displ_maskint(:,:)
real(kind_real) :: lon(nam%nc1,nam%nc3,geom%nl0),lat(nam%nc1,nam%nc3,geom%nl0)
character(len=3) :: il0ichar
character(len=1024) :: subr = 'samp_write'

! Processor verification
if (.not.mpl%main) call mpl%abort('only I/O proc should enter '//trim(subr))

! Create file
write(mpl%info,'(a7,a)') '','Write sampling'
call flush(mpl%info)
call mpl%ncerr(subr,nf90_create(trim(nam%datadir)//'/'//trim(nam%prefix)//'_sampling.nc',or(nf90_clobber,nf90_64bit_offset),ncid))

! Write namelist parameters
call nam%ncwrite(mpl,ncid)

! Define dimensions
call mpl%ncerr(subr,nf90_def_dim(ncid,'nl0',geom%nl0,nl0_id))
call mpl%ncerr(subr,nf90_put_att(ncid,nf90_global,'nl0r',bpar%nl0rmax))
call mpl%ncerr(subr,nf90_def_dim(ncid,'nc3',nam%nc3,nc3_id))
call mpl%ncerr(subr,nf90_def_dim(ncid,'nc1',nam%nc1,nc1_id))
if (nam%new_vbal.or.nam%new_lct.or.(nam%new_hdiag.and.(nam%var_diag.or.nam%local_diag.or.nam%displ_diag))) &
 & call mpl%ncerr(subr,nf90_def_dim(ncid,'nc2',nam%nc2,nc2_id))

! Define variables
call mpl%ncerr(subr,nf90_def_var(ncid,'lat',ncfloat,(/nc1_id,nc3_id,nl0_id/),lat_id))
call mpl%ncerr(subr,nf90_put_att(ncid,lat_id,'_FillValue',msvalr))
call mpl%ncerr(subr,nf90_def_var(ncid,'lon',ncfloat,(/nc1_id,nc3_id,nl0_id/),lon_id))
call mpl%ncerr(subr,nf90_put_att(ncid,lon_id,'_FillValue',msvalr))
call mpl%ncerr(subr,nf90_def_var(ncid,'smax',ncfloat,(/nc3_id,nl0_id/),smax_id))
call mpl%ncerr(subr,nf90_put_att(ncid,smax_id,'_FillValue',msvalr))
call mpl%ncerr(subr,nf90_def_var(ncid,'c1_to_c0',nf90_int,(/nc1_id/),c1_to_c0_id))
call mpl%ncerr(subr,nf90_put_att(ncid,c1_to_c0_id,'_FillValue',msvali))
call mpl%ncerr(subr,nf90_def_var(ncid,'c1l0_log',nf90_int,(/nc1_id,nl0_id/),c1l0_log_id))
call mpl%ncerr(subr,nf90_put_att(ncid,c1l0_log_id,'_FillValue',msvali))
call mpl%ncerr(subr,nf90_def_var(ncid,'c1c3_to_c0',nf90_int,(/nc1_id,nc3_id/),c1c3_to_c0_id))
call mpl%ncerr(subr,nf90_put_att(ncid,c1c3_to_c0_id,'_FillValue',msvali))
call mpl%ncerr(subr,nf90_def_var(ncid,'c1c3l0_log',nf90_int,(/nc1_id,nc3_id,nl0_id/),c1c3l0_log_id))
call mpl%ncerr(subr,nf90_put_att(ncid,c1c3l0_log_id,'_FillValue',msvali))
if (nam%new_lct.or.nam%var_diag.or.nam%local_diag.or.nam%displ_diag) then
   call mpl%ncerr(subr,nf90_def_var(ncid,'c2_to_c1',nf90_int,(/nc2_id/),c2_to_c1_id))
   call mpl%ncerr(subr,nf90_put_att(ncid,c2_to_c1_id,'_FillValue',msvali))
   call mpl%ncerr(subr,nf90_def_var(ncid,'c2_to_c0',nf90_int,(/nc2_id/),c2_to_c0_id))
   call mpl%ncerr(subr,nf90_put_att(ncid,c2_to_c0_id,'_FillValue',msvali))
end if

! End definition mode
call mpl%ncerr(subr,nf90_enddef(ncid))

! Convert data
call msr(lon)
call msr(lat)
do il0=1,geom%nl0
   do jc3=1,nam%nc3
      do ic1=1,nam%nc1
         if (samp%c1c3l0_log(ic1,jc3,il0)) then
            lon(ic1,jc3,il0) = geom%lon(samp%c1c3_to_c0(ic1,jc3))*rad2deg
            lat(ic1,jc3,il0) = geom%lat(samp%c1c3_to_c0(ic1,jc3))*rad2deg
            c1c3l0_logint(ic1,jc3,il0) = 1
         else
            c1c3l0_logint(ic1,jc3,il0) = 0
         end if
      end do
   end do
   do ic1=1,nam%nc1
      if (samp%c1l0_log(ic1,il0)) then
         c1l0_logint(ic1,il0) = 1
      else
         c1l0_logint(ic1,il0) = 0
      end if
   end do
end do

! Write variables
call mpl%ncerr(subr,nf90_put_var(ncid,lon_id,lon))
call mpl%ncerr(subr,nf90_put_var(ncid,lat_id,lat))
call mpl%ncerr(subr,nf90_put_var(ncid,smax_id,real(count(samp%c1c3l0_log,dim=1),kind_real)))
call mpl%ncerr(subr,nf90_put_var(ncid,c1_to_c0_id,samp%c1_to_c0))
call mpl%ncerr(subr,nf90_put_var(ncid,c1l0_log_id,c1l0_logint))
call mpl%ncerr(subr,nf90_put_var(ncid,c1c3_to_c0_id,samp%c1c3_to_c0))
call mpl%ncerr(subr,nf90_put_var(ncid,c1c3l0_log_id,c1c3l0_logint))
if (nam%new_lct.or.nam%var_diag.or.nam%local_diag.or.nam%displ_diag) then
   call mpl%ncerr(subr,nf90_put_var(ncid,c2_to_c1_id,samp%c2_to_c1))
   call mpl%ncerr(subr,nf90_put_var(ncid,c2_to_c0_id,samp%c2_to_c0))
end if

! Close file
call mpl%ncerr(subr,nf90_close(ncid))

! Write nearest neighbors and interpolation
do il0i=1,geom%nl0i
   if (nam%new_vbal.or.nam%new_lct.or.(nam%new_hdiag.and.(nam%var_diag.or.nam%local_diag.or.nam%displ_diag))) then
      ! Create file
      write(il0ichar,'(i3.3)') il0i
      call mpl%ncerr(subr,nf90_create(trim(nam%datadir)//'/'//trim(nam%prefix)//'_sampling_'//il0ichar//'.nc', &
    & or(nf90_clobber,nf90_64bit_offset),ncid))

      ! Write namelist parameters
      call nam%ncwrite(mpl,ncid)

      ! Define dimensions
      call mpl%ncerr(subr,nf90_def_dim(ncid,'nc1',nam%nc1,nc1_id))
      call mpl%ncerr(subr,nf90_def_dim(ncid,'nc2_1',nam%nc2,nc2_1_id))
      call mpl%ncerr(subr,nf90_def_dim(ncid,'nc2_2',nam%nc2,nc2_2_id))

      ! Define variables
      call mpl%ncerr(subr,nf90_def_var(ncid,'nn_c2_index',nf90_int,(/nc2_1_id,nc2_2_id/),nn_c2_index_id))
      call mpl%ncerr(subr,nf90_put_att(ncid,nn_c2_index_id,'_FillValue',msvali))
      call mpl%ncerr(subr,nf90_def_var(ncid,'nn_c2_dist',ncfloat,(/nc2_1_id,nc2_2_id/),nn_c2_dist_id))
      call mpl%ncerr(subr,nf90_put_att(ncid,nn_c2_dist_id,'_FillValue',msvalr))

      ! End definition mode
      call mpl%ncerr(subr,nf90_enddef(ncid))

      ! Write variables
      call mpl%ncerr(subr,nf90_put_var(ncid,nn_c2_index_id,samp%nn_c2_index(:,:,il0i)))
      call mpl%ncerr(subr,nf90_put_var(ncid,nn_c2_dist_id,samp%nn_c2_dist(:,:,il0i)))
      call samp%hfull(il0i)%write(mpl,ncid)
   end if

   if (nam%new_vbal) then
      ! Allocation
      allocate(vbal_maskint(nam%nc1,nam%nc2))

      ! Definition mode
      call mpl%ncerr(subr,nf90_redef(ncid))

      ! Define variables
      call mpl%ncerr(subr,nf90_def_var(ncid,'vbal_mask',nf90_int,(/nc1_id,nc2_1_id/),vbal_mask_id))
      call mpl%ncerr(subr,nf90_put_att(ncid,vbal_mask_id,'_FillValue',msvali))

      ! Convert data
      do ic2=1,nam%nc2
         do ic1=1,nam%nc1
            if (samp%vbal_mask(ic1,ic2)) then
               vbal_maskint(ic1,ic2) = 1
            else
               vbal_maskint(ic1,ic2) = 0
            end if
         end do
      end do

      ! End definition mode
      call mpl%ncerr(subr,nf90_enddef(ncid))

      ! Write variables
      call mpl%ncerr(subr,nf90_put_var(ncid,vbal_mask_id,vbal_maskint))

      ! Release memory
      deallocate(vbal_maskint)
   end if

   if (nam%local_diag) then
      ! Allocation
      allocate(local_maskint(nam%nc1,nam%nc2))

      ! Definition mode
      call mpl%ncerr(subr,nf90_redef(ncid))

      ! Define variables
      call mpl%ncerr(subr,nf90_def_var(ncid,'local_mask',nf90_int,(/nc1_id,nc2_1_id/),local_mask_id))
      call mpl%ncerr(subr,nf90_put_att(ncid,local_mask_id,'_FillValue',msvali))

      ! Convert data
      do ic2=1,nam%nc2
         do ic1=1,nam%nc1
            if (samp%local_mask(ic1,ic2)) then
               local_maskint(ic1,ic2) = 1
            else
               local_maskint(ic1,ic2) = 0
            end if
         end do
      end do

      ! End definition mode
      call mpl%ncerr(subr,nf90_enddef(ncid))

      ! Write variables
      call mpl%ncerr(subr,nf90_put_var(ncid,local_mask_id,local_maskint))

      ! Release memory
      deallocate(local_maskint)
   end if

   if (nam%displ_diag) then
      ! Allocation
      allocate(displ_maskint(nam%nc1,nam%nc2))

      ! Definition mode
      call mpl%ncerr(subr,nf90_redef(ncid))

      ! Define variables
      call mpl%ncerr(subr,nf90_def_var(ncid,'displ_mask',nf90_int,(/nc1_id,nc2_1_id/),displ_mask_id))
      call mpl%ncerr(subr,nf90_put_att(ncid,displ_mask_id,'_FillValue',msvali))

      ! Convert data
      do ic2=1,nam%nc2
         do ic1=1,nam%nc1
            if (samp%displ_mask(ic1,ic2)) then
               displ_maskint(ic1,ic2) = 1
            else
               displ_maskint(ic1,ic2) = 0
            end if
         end do
      end do

      ! End definition mode
      call mpl%ncerr(subr,nf90_enddef(ncid))

      ! Write variables
      call mpl%ncerr(subr,nf90_put_var(ncid,displ_mask_id,displ_maskint))

      ! Release memory
      deallocate(displ_maskint)
   end if

   if (nam%new_lct.or.nam%var_diag.or.nam%local_diag.or.nam%displ_diag) then
      ! Close file
      call mpl%ncerr(subr,nf90_close(ncid))
   end if
end do

end subroutine samp_write

!----------------------------------------------------------------------
! Subroutine: samp_setup_sampling
! Purpose: setup sampling
!----------------------------------------------------------------------
subroutine samp_setup_sampling(samp,mpl,rng,nam,geom,bpar,io,ens)

implicit none

! Passed variables
class(samp_type),intent(inout) :: samp ! Sampling
type(mpl_type),intent(inout) :: mpl    ! MPI data
type(rng_type),intent(inout) :: rng    ! Random number generator
type(nam_type),intent(inout) :: nam    ! Namelist
type(geom_type),intent(in) :: geom     ! Geometry
type(bpar_type),intent(in) :: bpar     ! Block parameters
type(io_type),intent(in) :: io         ! I/O
type(ens_type),intent(in) :: ens       ! Ensemble

! Local variables
integer :: ios,ic0,il0,ic1,ic2,ildw,jc3,il0i,jc1,kc1
integer,allocatable :: vbot(:),vtop(:),nn_c1_index(:)
real(kind_real) :: lon_c1(nam%nc1),lat_c1(nam%nc1)
real(kind_real) :: lon_c2(nam%nc2),lat_c2(nam%nc2)
real(kind_real) :: rh_c0a(geom%nc0a,geom%nl0),nn_dist(1)
real(kind_real),allocatable :: rh_c1(:),nn_c1_dist(:)
logical :: mask_c1(nam%nc1),mask_c2(nam%nc2)
character(len=1024) :: filename
type(kdtree_type) :: kdtree
type(linop_type) :: hbase

! Check subsampling size
if (nam%nc1>maxval(geom%nc0_mask)) then
   call mpl%warning('nc1 is too large for then mask, reset nc1 to the largest possible value')
   nam%nc1 = maxval(geom%nc0_mask)
end if

! Allocation
call samp%alloc(nam,geom)

! Read or compute sampling data
ios = 1
if (nam%sam_read) call samp%read(mpl,nam,geom,bpar,ios)
if (ios==1) then
   ! Compute mask
   call samp%compute_mask(mpl,nam,geom,ens)

   ! Compute zero-separation sampling
   call samp%compute_sampling_zs(mpl,rng,nam,geom)

   if (nam%new_lct) then
      ! Compute LCT sampling
      call samp%compute_sampling_lct(mpl,nam,geom)
   elseif (nam%new_vbal.or.nam%new_hdiag) then
      ! Compute positive separation sampling
      call samp%compute_sampling_ps(mpl,rng,nam,geom)
   end if

   ! Compute sampling mask
   call samp%compute_sampling_mask(mpl,nam,geom)
end if

if (nam%new_vbal.or.nam%new_lct.or.(nam%new_hdiag.and.(nam%var_diag.or.nam%local_diag.or.nam%displ_diag))) then
   if ((ios==1).or.(ios==2)) then
      ! Define subsampling
      write(mpl%info,'(a7,a)',advance='no') '','Define subsampling:'
      call flush(mpl%info)
      allocate(rh_c1(nam%nc1))
      lon_c1 = geom%lon(samp%c1_to_c0)
      lat_c1 = geom%lat(samp%c1_to_c0)
      mask_c1 = any(samp%c1l0_log(:,:),dim=2)
      rh_c1 = 1.0
      call rng%initialize_sampling(mpl,nam%nc1,lon_c1,lat_c1,mask_c1,rh_c1,nam%ntry,nam%nrep,nam%nc2,samp%c2_to_c1)
      samp%c2_to_c0 = samp%c1_to_c0(samp%c2_to_c1)
      deallocate(rh_c1)
   end if

   if ((ios==1).or.(ios==2).or.(ios==3).or.(ios==4)) then
      ! Find nearest neighbors
      write(mpl%info,'(a7,a)') '','Find nearest neighbors'
      call flush(mpl%info)
      do il0=1,geom%nl0
         if ((il0==1).or.(geom%nl0i>1)) then
            write(mpl%info,'(a10,a,i3)') '','Level ',nam%levs(il0)
            call flush(mpl%info)
            lon_c2 = geom%lon(samp%c2_to_c0)
            lat_c2 = geom%lat(samp%c2_to_c0)
            mask_c2 = .true.
            call kdtree%create(mpl,nam%nc2,lon_c2,lat_c2,mask=mask_c2)
            do ic2=1,nam%nc2
               ic1 = samp%c2_to_c1(ic2)
               ic0 = samp%c2_to_c0(ic2)
               call kdtree%find_nearest_neighbors(geom%lon(ic0),geom%lat(ic0),nam%nc2,samp%nn_c2_index(:,ic2,il0), &
             & samp%nn_c2_dist(:,ic2,il0))
            end do
            call kdtree%dealloc
         end if
      end do
   end if

   ! Compute sampling mesh
   write(mpl%info,'(a7,a)') '','Compute sampling mesh'
   call flush(mpl%info)
   lon_c2 = geom%lon(samp%c2_to_c0)
   lat_c2 = geom%lat(samp%c2_to_c0)
   call samp%mesh%create(mpl,rng,nam%nc2,lon_c2,lat_c2)

   ! Compute triangles list
   write(mpl%info,'(a7,a)') '','Compute triangles list '
   call flush(mpl%info)
   call samp%mesh%trlist

   ! Find boundary nodes
   write(mpl%info,'(a7,a)') '','Find boundary nodes'
   call flush(mpl%info)
   call samp%mesh%bnodes

   ! Find boundary arcs
   write(mpl%info,'(a7,a)') '','Find boundary arcs'
   call flush(mpl%info)
   call samp%mesh%barcs

   if ((ios==1).or.(ios==2).or.(ios==3)) then
      if (nam%new_vbal.or.nam%local_diag.or.nam%displ_diag) then
         ! Allocation
         allocate(nn_c1_index(nam%nc1))
         allocate(nn_c1_dist(nam%nc1))

         ! Compute local masks
         write(mpl%info,'(a7,a)') '','Compute local masks'
         call flush(mpl%info)
         lon_c1 = geom%lon(samp%c1_to_c0)
         lat_c1 = geom%lat(samp%c1_to_c0)
         mask_c1 = any(samp%c1l0_log,dim=2)
         call kdtree%create(mpl,nam%nc1,lon_c1,lat_c1,mask=mask_c1)
         do ic2=1,nam%nc2
            ! Inidices
            ic1 = samp%c2_to_c1(ic2)
            ic0 = samp%c2_to_c0(ic2)

            ! Find nearest neighbors
            call kdtree%find_nearest_neighbors(geom%lon(ic0),geom%lat(ic0),nam%nc1,nn_c1_index,nn_c1_dist)
            do jc1=1,nam%nc1
               kc1 = nn_c1_index(jc1)
               if (nam%new_vbal) samp%vbal_mask(kc1,ic2) = (jc1==1) &
             & .or.(nn_c1_dist(jc1)<nam%vbal_rad)
               if (nam%local_diag) samp%local_mask(kc1,ic2) = (jc1==1) &
             & .or.(nn_c1_dist(jc1)<nam%local_rad)
               if (nam%displ_diag) samp%displ_mask(kc1,ic2) = (jc1==1) &
             & .or.(nn_c1_dist(jc1)<min(nam%displ_rad,samp%mesh%bdist(ic2)))
            end do
         end do
         call kdtree%dealloc

         ! Release memory
         deallocate(nn_c1_index)
         deallocate(nn_c1_dist)
      end if

      ! Allocation
      allocate(vbot(nam%nc2))
      allocate(vtop(nam%nc2))

      ! Initialize vbot and vtop
      vbot = 1
      vtop = geom%nl0

      ! Compute grid interpolation
      write(mpl%info,'(a7,a)') '','Compute grid interpolation'
      do il0i=1,geom%nl0i
         ! Compute grid interpolation
         write(samp%hfull(il0i)%prefix,'(a,i3.3)') 'hfull_',il0i
         call samp%hfull(il0i)%interp(mpl,rng,geom,il0i,nam%nc2,samp%c2_to_c0,nam%mask_check,vbot,vtop,nam%diag_interp,hbase)
      end do

      ! Release memory
      deallocate(vbot)
      deallocate(vtop)
   end if
end if

! Write sampling data
if (nam%sam_write) then
   if (mpl%main) call samp%write(mpl,nam,geom,bpar)

   ! Write rh_c0
   if (trim(nam%draw_type)=='random_coast') then
      call mpl%glb_to_loc(geom%nl0,geom%nc0,geom%c0_to_proc,geom%c0_to_c0a,samp%rh_c0,geom%nc0a,rh_c0a)
      filename = trim(nam%prefix)//'_sampling_rh_c0'
      call io%fld_write(mpl,nam,geom,filename,'rh_c0',rh_c0a)
   end if
end if

! Compute nearest neighbors for local diagnostics output
if (nam%local_diag.and.(nam%nldwv>0)) then
   write(mpl%info,'(a7,a)') '','Compute nearest neighbors for local diagnostics output'
   call flush(mpl%info)
   allocate(samp%nn_ldwv_index(nam%nldwv))
   call kdtree%create(mpl,nam%nc2,geom%lon(samp%c2_to_c0),geom%lat(samp%c2_to_c0),mask=samp%c1l0_log(samp%c2_to_c1,1))
   do ildw=1,nam%nldwv
      call kdtree%find_nearest_neighbors(nam%lon_ldwv(ildw),nam%lat_ldwv(ildw),1,samp%nn_ldwv_index(ildw:ildw),nn_dist)
   end do
   call kdtree%dealloc
end if

! Print results
write(mpl%info,'(a7,a)') '','Sampling efficiency (%):'
do il0=1,geom%nl0
   write(mpl%info,'(a10,a,i3,a)',advance='no') '','Level ',nam%levs(il0),' ~> '
   do jc3=1,nam%nc3
      if (count(samp%c1c3l0_log(:,jc3,il0))>=nam%nc1/2) then
         ! Sucessful sampling
         write(mpl%info,'(a,i3,a)',advance='no') trim(mpl%green), &
       & int(100.0*real(count(samp%c1c3l0_log(:,jc3,il0)),kind_real)/real(nam%nc1,kind_real)),trim(mpl%black)
      else
         ! Insufficient sampling
         write(mpl%info,'(a,i3,a)',advance='no') trim(mpl%peach), &
       & int(100.0*real(count(samp%c1c3l0_log(:,jc3,il0)),kind_real)/real(nam%nc1,kind_real)),trim(mpl%black)
      end if
      if (jc3<nam%nc3) write(mpl%info,'(a)',advance='no') '-'
   end do
   write(mpl%info,'(a)') ' '
end do
call flush(mpl%info)

end subroutine samp_setup_sampling

!----------------------------------------------------------------------
! Subroutine: samp_compute_mask
! Purpose: compute mask
!----------------------------------------------------------------------
subroutine samp_compute_mask(samp,mpl,nam,geom,ens)

implicit none

! Passed variables
class(samp_type),intent(inout) :: samp ! Sampling
type(mpl_type),intent(inout) :: mpl    ! MPI data
type(nam_type),intent(in) :: nam       ! Namelist
type(geom_type),intent(in) :: geom     ! Geometry
type(ens_type),intent(in) :: ens       ! Ensemble

! Local variables
integer :: ic0a,il0,iv,its,ie,isub,ie_sub
real(kind_real),allocatable :: var(:,:,:,:)
logical :: valid,mask_c0a(geom%nc0a,geom%nl0)

! Copy geometry mask
mask_c0a = geom%mask_c0a

! Ensemble-based mask
select case (trim(nam%mask_type))
case ('stddev')
   ! Allocation
   allocate(var(geom%nc0a,geom%nl0,nam%nv,nam%nts))

   ! Compute variances
   var = 0.0
   do ie=1,ens%ne
      var = var +ens%fld(:,:,:,:,ie)**2
   end do
   var = var/real(ens%ne-ens%nsub,kind_real)

   ! Check whether standard-deviation is over the threshold
   do its=1,nam%nts
      do iv=1,nam%nv
         mask_c0a = mask_c0a.and.(var(:,:,iv,its)>nam%mask_th**2)
      end do
   end do
case ('all_members')
   ! Check whether all members are over the threshold
   do its=1,nam%nts
      do iv=1,nam%nv
         do il0=1,geom%nl0
            do ic0a=1,geom%nc0a
               valid = .true.
               do isub=1,ens%nsub
                  do ie_sub=1,ens%ne/ens%nsub
                     ie = ie_sub+(isub-1)*ens%ne/ens%nsub
                     valid = valid.and.(ens%mean(ic0a,il0,iv,its,isub)+ens%fld(ic0a,il0,iv,its,ie)>nam%mask_th)
                  end do
               end do
               mask_c0a(ic0a,il0) = mask_c0a(ic0a,il0).and.valid
            end do
         end do
      end do
   end do
end select

! Local to global
call mpl%loc_to_glb(geom%nl0,geom%nc0a,mask_c0a,geom%nc0,geom%c0_to_proc,geom%c0_to_c0a,.true.,samp%mask_c0)

! Other masks
samp%mask_hor_c0 = any(samp%mask_c0,dim=2)
samp%nc0_mask = count(samp%mask_c0,dim=1)

! Print results
select case (trim(nam%mask_type))
case ('stddev','all_members','any_member')
   write(mpl%info,'(a7,a)') '','Number of points excluded from HDIAG sampling:'
   do il0=1,geom%nl0
      write(mpl%info,'(a10,a,i3,a,i8,a,i8)') '','Level',nam%levs(il0),': ',count(geom%mask_c0(:,il0))-count(samp%mask_c0(:,il0)), &
    & ' of ',count(geom%mask_c0(:,il0))
   end do
end select

end subroutine samp_compute_mask

!----------------------------------------------------------------------
! Subroutine: samp_compute_sampling_zs
! Purpose: compute zero-separation sampling
!----------------------------------------------------------------------
subroutine samp_compute_sampling_zs(samp,mpl,rng,nam,geom)

implicit none

! Passed variables
class(samp_type),intent(inout) :: samp ! Sampling
type(mpl_type),intent(inout) :: mpl    ! MPI data
type(rng_type),intent(inout) :: rng    ! Random number generator
type(nam_type),intent(in) :: nam       ! Namelist
type(geom_type),intent(in) :: geom     ! Geometry

! Local variables
integer :: ic0,ic1,il0
integer :: nn_index(1)
real(kind_real) :: nn_dist(1)
type(kdtree_type) :: kdtree

! Compute subset
if (nam%nc1<maxval(samp%nc0_mask)) then
   write(mpl%info,'(a7,a)',advance='no') '','Compute horizontal subset C1: '
   call flush(mpl%info)
   select case (trim(nam%draw_type))
   case ('random_uniform')
      ! Random draw
      do ic0=1,geom%nc0
         if (samp%mask_hor_c0(ic0)) samp%rh_c0(ic0,1) = 1.0
      end do
   case ('random_coast')
      ! More points around coasts
      do ic0=1,geom%nc0
         if (samp%mask_hor_c0(ic0)) samp%rh_c0(ic0,1) = 0.0
      end do
      do il0=1,geom%nl0
         call kdtree%create(mpl,geom%nc0,geom%lon,geom%lat,mask=.not.samp%mask_c0(:,il0))
         do ic0=1,geom%nc0
            if (samp%mask_c0(ic0,il0)) then
               call kdtree%find_nearest_neighbors(geom%lon(ic0),geom%lat(ic0),1,nn_index,nn_dist)
               samp%rh_c0(ic0,1) = samp%rh_c0(ic0,1)+exp(-nn_dist(1)/lcoast)
            else
               samp%rh_c0(ic0,1) = samp%rh_c0(ic0,1)+1.0
            end if
         end do
         call kdtree%dealloc
      end do
      samp%rh_c0(:,1) = rcoast+(1.0-rcoast)*(1.0-samp%rh_c0(:,1)/real(geom%nl0,kind_real))
   end select

   ! Initialize sampling
   call rng%initialize_sampling(mpl,geom%nc0,geom%lon,geom%lat,samp%mask_hor_c0,samp%rh_c0(:,1),nam%ntry,nam%nrep, &
 & nam%nc1,samp%c1_to_c0)
else
   ic1 = 0
   do ic0=1,geom%nc0
      if (samp%mask_hor_c0(ic0)) then
         ic1 = ic1+1
         samp%c1_to_c0(ic1) = ic0
      end if
   end do
end if

end subroutine samp_compute_sampling_zs

!----------------------------------------------------------------------
! Subroutine: samp_compute_sampling_ps
! Purpose: compute positive separation sampling
!----------------------------------------------------------------------
subroutine samp_compute_sampling_ps(samp,mpl,rng,nam,geom)

implicit none

! Passed variables
class(samp_type),intent(inout) :: samp ! Sampling
type(mpl_type),intent(inout) :: mpl    ! MPI data
type(rng_type),intent(inout) :: rng    ! Random number generator
type(nam_type),intent(in) :: nam       ! Namelist
type(geom_type),intent(in) :: geom     ! Geometry

! Local variables
integer :: irmaxloc,jc3,ic1,ir,ic0,jc0,i,nvc0,ivc0,icinf,icsup,ictest
integer,allocatable :: vic0(:)
real(kind_real) :: d
real(kind_real),allocatable :: x(:),y(:),z(:),v1(:),v2(:),va(:),vp(:),t(:)
logical :: found

! First class
samp%c1c3_to_c0(:,1) = samp%c1_to_c0

if (nam%nc3>1) then
   write(mpl%info,'(a7,a)',advance='no') '','Compute positive separation sampling: '
   call flush(mpl%info)

   ! Initialize
   do jc3=1,nam%nc3
      if (jc3/=1) call msi(samp%c1c3_to_c0(:,jc3))
   end do

   ! Define valid nodes vector
   nvc0 = count(samp%mask_hor_c0)
   allocate(vic0(nvc0))
   ivc0 = 0
   do ic0=1,geom%nc0
      if (samp%mask_hor_c0(ic0)) then
         ivc0 = ivc0+1
         vic0(ivc0) = ic0
      end if
   end do

   ! Sample classes of positive separation
   call mpl%prog_init(nam%nc3*nam%nc1)
   ir = 0
   irmaxloc = irmax
   do while ((.not.all(isnotmsi(samp%c1c3_to_c0))).and.(nvc0>1).and.(ir<=irmaxloc))
      ! Try a random point
      if (mpl%main) call rng%rand_integer(1,nvc0,i)
      call mpl%f_comm%broadcast(i,mpl%ioproc-1)
      ir = ir+1
      jc0 = vic0(i)

      !$omp parallel do schedule(static) private(ic1,ic0,d,jc3,icinf,icsup,found,ictest) firstprivate(x,y,z,v1,v2,va,vp,t)
      do ic1=1,nam%nc1
         ! Allocation
         allocate(x(2))
         allocate(y(2))
         allocate(z(2))
         allocate(v1(3))
         allocate(v2(3))
         allocate(va(3))
         allocate(vp(3))
         allocate(t(4))

         ! Check if there is a valid first point
         if (isnotmsi(samp%c1_to_c0(ic1))) then
            ! Compute the distance
            ic0 = samp%c1_to_c0(ic1)
            call sphere_dist(geom%lon(ic0),geom%lat(ic0),geom%lon(jc0),geom%lat(jc0),d)

            ! Find the class (dichotomy method)
            if ((d>0.0).and.(d<(real(nam%nc3,kind_real)-0.5)*nam%dc)) then
               jc3 = 1
               icinf = 1
               icsup = nam%nc3
               found = .false.
               do while (.not.found)
                  ! New value
                  ictest = (icsup+icinf)/2

                  ! Update
                  if (d<(real(ictest,kind_real)-0.5)*nam%dc) icsup = ictest
                  if (d>(real(ictest,kind_real)-0.5)*nam%dc) icinf = ictest

                  ! Exit test
                  if (icsup==icinf+1) then
                     if (abs((real(icinf,kind_real)-0.5)*nam%dc-d)<abs((real(icsup,kind_real)-0.5)*nam%dc-d)) then
                        jc3 = icinf
                     else
                        jc3 = icsup
                     end if
                     found = .true.
                  end if
               end do

               ! Find if this class has not been aready filled
               if ((jc3/=1).and.(ismsi(samp%c1c3_to_c0(ic1,jc3)))) samp%c1c3_to_c0(ic1,jc3) = jc0
            end if
         end if

         ! Release memory
         deallocate(x)
         deallocate(y)
         deallocate(z)
         deallocate(v1)
         deallocate(v2)
         deallocate(va)
         deallocate(vp)
         deallocate(t)
      end do
      !$omp end parallel do

      ! Update valid nodes vector
      vic0(i) = vic0(nvc0)
      nvc0 = nvc0-1

      ! Update
      mpl%done = pack(isnotmsi(samp%c1c3_to_c0),mask=.true.)
      call mpl%prog_print
   end do
   write(mpl%info,'(a)') '100%'
   call flush(mpl%info)

   ! Release memory
   deallocate(vic0)
end if

end subroutine samp_compute_sampling_ps

!----------------------------------------------------------------------
! Subroutine: samp_compute_sampling_lct
! Purpose: compute LCT sampling
!----------------------------------------------------------------------
subroutine samp_compute_sampling_lct(samp,mpl,nam,geom)

implicit none

! Passed variables
class(samp_type),intent(inout) :: samp ! Sampling
type(mpl_type),intent(inout) :: mpl    ! MPI data
type(nam_type),intent(in) :: nam       ! Namelist
type(geom_type),intent(in) :: geom     ! Geometry

! Local variables
integer :: i,il0,ic1,ic0,jc0,ibnd,ic3
integer :: nn_index(nam%nc3)
integer :: iproc,ic1_s(mpl%nproc),ic1_e(mpl%nproc),nc1_loc(mpl%nproc),ic1_loc
integer,allocatable :: sbufi(:),rbufi(:)
real(kind_real) :: nn_dist(nam%nc3)
real(kind_real),allocatable :: x(:),y(:),z(:),v1(:),v2(:),va(:),vp(:),t(:)
logical,allocatable :: sbufl(:),rbufl(:)
type(fckit_mpi_status) :: status

write(mpl%info,'(a7,a)',advance='no') '','Compute LCT sampling: '
call flush(mpl%info)

! MPI splitting
call mpl%split(nam%nc1,ic1_s,ic1_e,nc1_loc)

! Initialization
call mpl%prog_init(nc1_loc(mpl%myproc))

do ic1_loc=1,nc1_loc(mpl%myproc)
   ! MPI offset
   ic1 = ic1_s(mpl%myproc)+ic1_loc-1

   ! Check location validity
   if (isnotmsi(samp%c1_to_c0(ic1))) then
      ! Find neighbors
      call geom%kdtree%find_nearest_neighbors(geom%lon(samp%c1_to_c0(ic1)),geom%lat(samp%c1_to_c0(ic1)), &
    & nam%nc3,nn_index,nn_dist)

      ! Copy neighbor index
      do ic3=1,nam%nc3
         jc0 = nn_index(ic3)
         samp%c1c3_to_c0(ic1,ic3) = nn_index(ic3)
         do il0=1,geom%nl0
            samp%c1c3l0_log(ic1,ic3,il0) = samp%mask_c0(jc0,il0)
         end do
      end do

      if (nam%mask_check) then
         ! Check that great circle to neighbors is not crossing mask boundaries
         do il0=1,geom%nl0
            !$omp parallel do schedule(static) private(ic3,ic0,jc0) firstprivate(x,y,z,v1,v2,va,vp,t)
            do ic3=1,nam%nc3
               ! Allocation
               allocate(x(2))
               allocate(y(2))
               allocate(z(2))
               allocate(v1(3))
               allocate(v2(3))
               allocate(va(3))
               allocate(vp(3))
               allocate(t(4))

               ! Indices
               ic0 = samp%c1_to_c0(ic1)
               jc0 = samp%c1c3_to_c0(ic1,ic3)

               ! Transform to cartesian coordinates
               call trans(2,geom%lat((/ic0,jc0/)),geom%lon((/ic0,jc0/)),x,y,z)

               ! Compute arc orthogonal vector
               v1 = (/x(1),y(1),z(1)/)
               v2 = (/x(2),y(2),z(2)/)
               call vector_product(v1,v2,va)

               ! Check if arc is crossing boundary arcs
               do ibnd=1,geom%nbnd(il0)
                  call vector_product(va,geom%vbnd(:,ibnd,il0),vp)
                  v1 = (/x(1),y(1),z(1)/)
                  call vector_triple_product(v1,va,vp,t(1))
                  v1 = (/x(2),y(2),z(2)/)
                  call vector_triple_product(v1,va,vp,t(2))
                  v1 = (/geom%xbnd(1,ibnd,il0),geom%ybnd(1,ibnd,il0),geom%zbnd(1,ibnd,il0)/)
                  call vector_triple_product(v1,geom%vbnd(:,ibnd,il0),vp,t(3))
                  v1 = (/geom%xbnd(2,ibnd,il0),geom%ybnd(2,ibnd,il0),geom%zbnd(2,ibnd,il0)/)
                  call vector_triple_product(v1,geom%vbnd(:,ibnd,il0),vp,t(4))
                  t(1) = -t(1)
                  t(3) = -t(3)
                  if (all(t>0).or.(all(t<0))) then
                     samp%c1c3l0_log(ic1,ic3,il0) = .false.
                     exit
                  end if
               end do

               ! Memory release
               deallocate(x)
               deallocate(y)
               deallocate(z)
               deallocate(v1)
               deallocate(v2)
               deallocate(va)
               deallocate(vp)
               deallocate(t)
            end do
            !$omp end parallel do
         end do
      end if
   end if

   ! Update
   call mpl%prog_print(ic1_loc)
end do
write(mpl%info,'(a)') '100%'
call flush(mpl%info)

! Communication
write(mpl%info,'(a7,a)') '','Communication'
call flush(mpl%info)
if (mpl%main) then
   do iproc=1,mpl%nproc
      if (iproc/=mpl%ioproc) then
         ! Allocation
         allocate(rbufi(nc1_loc(iproc)*nam%nc3))
         allocate(rbufl(nc1_loc(iproc)*nam%nc3*geom%nl0))

         ! Receive data on ioproc
         call mpl%f_comm%receive(rbufi,iproc-1,mpl%tag,status)
         call mpl%f_comm%receive(rbufl,iproc-1,mpl%tag+1,status)

         ! Format data
         i = 0
         do ic3=1,nam%nc3
            do ic1_loc=1,nc1_loc(iproc)
               i = i+1
               ic1 = ic1_s(iproc)+ic1_loc-1
               samp%c1c3_to_c0(ic1,ic3) = rbufi(i)
            end do
         end do
         i = 0
         do il0=1,geom%nl0
            do ic3=1,nam%nc3
               do ic1_loc=1,nc1_loc(iproc)
                  i = i+1
                  ic1 = ic1_s(iproc)+ic1_loc-1
                  samp%c1c3l0_log(ic1,ic3,il0) = rbufl(i)
               end do
            end do
         end do

         ! Release memory
         deallocate(rbufi)
         deallocate(rbufl)
      end if
   end do
else
   ! Allocation
   allocate(sbufi(nc1_loc(mpl%myproc)*nam%nc3))
   allocate(sbufl(nc1_loc(mpl%myproc)*nam%nc3*geom%nl0))

   ! Prepare buffers
   i = 0
   do ic3=1,nam%nc3
      do ic1_loc=1,nc1_loc(mpl%myproc)
         i = i+1
         ic1 = ic1_s(mpl%myproc)+ic1_loc-1
         sbufi(i) = samp%c1c3_to_c0(ic1,ic3)
      end do
   end do
   i = 0
   do il0=1,geom%nl0
      do ic3=1,nam%nc3
         do ic1_loc=1,nc1_loc(mpl%myproc)
            i = i+1
            ic1 = ic1_s(mpl%myproc)+ic1_loc-1
            sbufl(i) = samp%c1c3l0_log(ic1,ic3,il0)
         end do
      end do
   end do

   ! Send data to ioproc
   call mpl%f_comm%send(sbufi,mpl%ioproc-1,mpl%tag)
   call mpl%f_comm%send(sbufl,mpl%ioproc-1,mpl%tag+1)

   ! Release memory
   deallocate(sbufi)
   deallocate(sbufl)
end if
call mpl%update_tag(2)
call mpl%f_comm%broadcast(samp%c1c3_to_c0,mpl%ioproc-1)
call mpl%f_comm%broadcast(samp%c1c3l0_log,mpl%ioproc-1)

end subroutine samp_compute_sampling_lct

!----------------------------------------------------------------------
! Subroutine: samp_compute_sampling_mask
! Purpose: compute sampling mask
!----------------------------------------------------------------------
subroutine samp_compute_sampling_mask(samp,mpl,nam,geom)

implicit none

! Passed variables
class(samp_type),intent(inout) :: samp ! Sampling
type(mpl_type),intent(inout) :: mpl    ! MPI data
type(nam_type),intent(in) :: nam       ! Namelist
type(geom_type),intent(in) :: geom     ! Geometry

! Local variables
integer :: il0,jc3,ic1,ic0,jc0,i
logical :: valid

! First point
do il0=1,geom%nl0
   samp%c1l0_log(:,il0) = samp%mask_c0(samp%c1_to_c0,il0)
end do

! Second point
if (nam%mask_check) then
   write(mpl%info,'(a7,a)',advance='no') '','Check mask boundaries: '
   call mpl%prog_init(geom%nl0*nam%nc3*nam%nc1)
   i = 0
end if
do il0=1,geom%nl0
   do jc3=1,nam%nc3
      do ic1=1,nam%nc1
         ! Indices
         ic0 = samp%c1c3_to_c0(ic1,jc3)
         jc0 = samp%c1c3_to_c0(ic1,1)

         ! Check point index
         valid = isnotmsi(ic0).and.isnotmsi(jc0)

         if (valid) then
            ! Check mask
            valid = samp%mask_c0(ic0,il0).and.samp%mask_c0(jc0,il0)

            ! Check mask bounds
            if (nam%mask_check) then
               if (valid) call geom%check_arc(il0,geom%lon(ic0),geom%lat(ic0),geom%lon(jc0),geom%lat(jc0),valid)
               i = i+1
               call mpl%prog_print(i)
            end if
         end if
         samp%c1c3l0_log(ic1,jc3,il0) = valid
      end do
   end do
end do
if (nam%mask_check) then
   write(mpl%info,'(a)') '100%'
   call flush(mpl%info)
end if

end subroutine samp_compute_sampling_mask

!----------------------------------------------------------------------
! Subroutine: samp_compute_mpi_a
! Purpose: compute sampling MPI distribution, halo A
!----------------------------------------------------------------------
subroutine samp_compute_mpi_a(samp,mpl,nam,geom)

implicit none

! Passed variables
class(samp_type),intent(inout) :: samp ! Sampling
type(mpl_type),intent(in) :: mpl       ! MPI data
type(nam_type),intent(in) :: nam       ! Namelist
type(geom_type),intent(in) :: geom     ! Geometry

! Local variables
integer :: ic0a,ic0,ic1a,ic1

! Allocation
allocate(samp%lcheck_c0a(geom%nc0))
allocate(samp%lcheck_c1a(nam%nc1))

! Halo definitions

! Halo A
samp%lcheck_c0a = .false.
samp%lcheck_c1a = .false.
do ic0a=1,geom%nc0a
   ic0 = geom%c0a_to_c0(ic0a)
   if (geom%c0_to_proc(ic0)==mpl%myproc) samp%lcheck_c0a(ic0) = .true.
end do
do ic1=1,nam%nc1
   ic0 = samp%c1_to_c0(ic1)
   if (geom%c0_to_proc(ic0)==mpl%myproc) samp%lcheck_c1a(ic1) = .true.
end do

! Halo sizes
samp%nc1a = count(samp%lcheck_c1a)

! Global <-> local conversions for fields

! Halo A
allocate(samp%c1a_to_c1(samp%nc1a))
allocate(samp%c1_to_c1a(nam%nc1))
ic1a = 0
do ic1=1,nam%nc1
   if (samp%lcheck_c1a(ic1)) then
      ic1a = ic1a+1
      samp%c1a_to_c1(ic1a) = ic1
      samp%c1_to_c1a(ic1) = ic1a
   end if
end do

! Print results
write(mpl%info,'(a7,a,i4)') '','Parameters for processor #',mpl%myproc
write(mpl%info,'(a10,a,i8)') '','nc0 =        ',geom%nc0
write(mpl%info,'(a10,a,i8)') '','nc0a =       ',geom%nc0a
write(mpl%info,'(a10,a,i8)') '','nl0 =        ',geom%nl0
write(mpl%info,'(a10,a,i8)') '','nc1 =        ',nam%nc1
write(mpl%info,'(a10,a,i8)') '','nc1a =       ',samp%nc1a
call flush(mpl%info)

end subroutine samp_compute_mpi_a

!----------------------------------------------------------------------
! Subroutine: samp_compute_mpi_ab
! Purpose: compute sampling MPI distribution, halos A-B
!----------------------------------------------------------------------
subroutine samp_compute_mpi_ab(samp,mpl,nam,geom)

implicit none

! Passed variables
class(samp_type),intent(inout) :: samp ! Sampling
type(mpl_type),intent(inout) :: mpl    ! MPI data
type(nam_type),intent(in) :: nam       ! Namelist
type(geom_type),intent(in) :: geom     ! Geometry

! Local variables
integer :: iproc,ic0,ic2a,ic2b,ic2,jc2,i_s,i_s_loc,h_n_s_max,il0i,h_n_s_max_loc
integer,allocatable :: interph_lg(:,:)

! Allocation
h_n_s_max = 0
do il0i=1,geom%nl0i
   h_n_s_max = max(h_n_s_max,samp%hfull(il0i)%n_s)
end do
allocate(samp%c2_to_proc(nam%nc2))
allocate(samp%proc_to_nc2a(mpl%nproc))
allocate(samp%h(geom%nl0i))
allocate(samp%lcheck_c2a(nam%nc2))
allocate(samp%lcheck_c2b(nam%nc2))
allocate(samp%lcheck_h(h_n_s_max,geom%nl0i))

! Halo definitions

! Halo A
samp%lcheck_c2a = .false.
do ic2=1,nam%nc2
   ic0 = samp%c2_to_c0(ic2)
   if (geom%c0_to_proc(ic0)==mpl%myproc) samp%lcheck_c2a(ic2) = .true.
end do

! Halo B
samp%lcheck_h = .false.
samp%lcheck_c2b = .false.
do il0i=1,geom%nl0i
   do i_s=1,samp%hfull(il0i)%n_s
      ic0 = samp%hfull(il0i)%row(i_s)
      iproc = geom%c0_to_proc(ic0)
      if (iproc==mpl%myproc) then
         jc2 = samp%hfull(il0i)%col(i_s)
         samp%lcheck_h(i_s,il0i) = .true.
         samp%lcheck_c2b(jc2) = .true.
      end if
   end do
end do

! Halo sizes
samp%nc2a = count(samp%lcheck_c2a)
do il0i=1,geom%nl0i
   samp%h(il0i)%n_s = count(samp%lcheck_h(:,il0i))
end do
samp%nc2b = count(samp%lcheck_c2b)

! Global <-> local conversions for fields

! Halo A
allocate(samp%c2a_to_c2(samp%nc2a))
allocate(samp%c2_to_c2a(nam%nc2))
ic2a = 0
do ic2=1,nam%nc2
   if (samp%lcheck_c2a(ic2)) then
      ic2a = ic2a+1
      samp%c2a_to_c2(ic2a) = ic2
   end if
end do
call mpl%glb_to_loc_index(samp%nc2a,samp%c2a_to_c2,nam%nc2,samp%c2_to_c2a)

! Halo B
allocate(samp%c2b_to_c2(samp%nc2b))
allocate(samp%c2_to_c2b(nam%nc2))
ic2b = 0
do ic2=1,nam%nc2
   if (samp%lcheck_c2b(ic2)) then
      ic2b = ic2b+1
      samp%c2b_to_c2(ic2b) = ic2
      samp%c2_to_c2b(ic2) = ic2b
   end if
end do

! Inter-halo conversions
allocate(samp%c2a_to_c2b(samp%nc2a))
do ic2a=1,samp%nc2a
   ic2 = samp%c2a_to_c2(ic2a)
   ic2b = samp%c2_to_c2b(ic2)
   samp%c2a_to_c2b(ic2a) = ic2b
end do

! Global <-> local conversions for data
h_n_s_max_loc = 0
do il0i=1,geom%nl0i
   h_n_s_max_loc = max(h_n_s_max_loc,samp%h(il0i)%n_s)
end do
allocate(interph_lg(h_n_s_max_loc,geom%nl0i))
do il0i=1,geom%nl0i
   i_s_loc = 0
   do i_s=1,samp%hfull(il0i)%n_s
      if (samp%lcheck_h(i_s,il0i)) then
         i_s_loc = i_s_loc+1
         interph_lg(i_s_loc,il0i) = i_s
      end if
   end do
end do

! Local data

! Horizontal interpolation
do il0i=1,geom%nl0i
   write(samp%h(il0i)%prefix,'(a,i3.3)') 'h_',il0i
   samp%h(il0i)%n_src = samp%nc2b
   samp%h(il0i)%n_dst = geom%nc0a
   call samp%h(il0i)%alloc
   do i_s_loc=1,samp%h(il0i)%n_s
      i_s = interph_lg(i_s_loc,il0i)
      samp%h(il0i)%row(i_s_loc) = geom%c0_to_c0a(samp%hfull(il0i)%row(i_s))
      samp%h(il0i)%col(i_s_loc) = samp%c2_to_c2b(samp%hfull(il0i)%col(i_s))
      samp%h(il0i)%S(i_s_loc) = samp%hfull(il0i)%S(i_s)
   end do
   call samp%h(il0i)%reorder(mpl)
end do

! MPI splitting
do ic2=1,nam%nc2
   ic0 = samp%c2_to_c0(ic2)
   samp%c2_to_proc(ic2) = geom%c0_to_proc(ic0)
end do
do iproc=1,mpl%nproc
   samp%proc_to_nc2a(iproc) = count(samp%c2_to_proc==iproc)
end do

! Setup communications
call samp%com_AB%setup(mpl,'com_AB',nam%nc2,samp%nc2a,samp%nc2b,samp%c2b_to_c2,samp%c2a_to_c2b,samp%c2_to_proc, &
 & samp%c2_to_c2a)

! Print results
write(mpl%info,'(a7,a,i4)') '','Parameters for processor #',mpl%myproc
write(mpl%info,'(a10,a,i8)') '','nc2 =        ',nam%nc2
write(mpl%info,'(a10,a,i8)') '','nc2a =       ',samp%nc2a
write(mpl%info,'(a10,a,i8)') '','nc2b =       ',samp%nc2b
do il0i=1,geom%nl0i
   write(mpl%info,'(a10,a,i3,a,i8)') '','h(',il0i,')%n_s = ',samp%h(il0i)%n_s
end do
call flush(mpl%info)

end subroutine samp_compute_mpi_ab

!----------------------------------------------------------------------
! Subroutine: samp_compute_mpi_d
! Purpose: compute sampling MPI distribution, halo D
!----------------------------------------------------------------------
subroutine samp_compute_mpi_d(samp,mpl,nam,geom)

implicit none

! Passed variables
class(samp_type),intent(inout) :: samp ! Sampling
type(mpl_type),intent(inout) :: mpl    ! MPI data
type(nam_type),intent(in) :: nam       ! Namelist
type(geom_type),intent(in) :: geom     ! Geometry

! Local variables
integer :: ic0,ic0a,ic0d,ic1,ic2,ic2a,jc0,jc1

! Allocation
allocate(samp%lcheck_c0d(geom%nc0))

! Halo definitions

! Halo D
samp%lcheck_c0d = samp%lcheck_c0a
do ic2a=1,samp%nc2a
   ic2 = samp%c2a_to_c2(ic2a)
   ic1 = samp%c2_to_c1(ic2)
   do jc1=1,nam%nc1
      if (samp%displ_mask(jc1,ic2)) then
         jc0 = samp%c1_to_c0(jc1)
         samp%lcheck_c0d(jc0) = .true.
      end if
   end do
end do
samp%nc0d = count(samp%lcheck_c0d)

! Global <-> local conversions for fields

! Halo D
allocate(samp%c0d_to_c0(samp%nc0d))
allocate(samp%c0_to_c0d(geom%nc0))
call msi(samp%c0_to_c0d)
ic0d = 0
do ic0=1,geom%nc0
   if (samp%lcheck_c0d(ic0)) then
      ic0d = ic0d+1
      samp%c0d_to_c0(ic0d) = ic0
      samp%c0_to_c0d(ic0) = ic0d
   end if
end do

! Inter-halo conversions
allocate(samp%c0a_to_c0d(geom%nc0a))
do ic0a=1,geom%nc0a
   ic0 = geom%c0a_to_c0(ic0a)
   ic0d = samp%c0_to_c0d(ic0)
   samp%c0a_to_c0d(ic0a) = ic0d
end do

! Setup communications
call samp%com_AD%setup(mpl,'com_AD',geom%nc0,geom%nc0a,samp%nc0d,samp%c0d_to_c0,samp%c0a_to_c0d,geom%c0_to_proc,geom%c0_to_c0a)

! Print results
write(mpl%info,'(a7,a,i4)') '','Parameters for processor #',mpl%myproc
write(mpl%info,'(a10,a,i8)') '','nc0d =       ',samp%nc0d
call flush(mpl%info)

end subroutine samp_compute_mpi_d

!----------------------------------------------------------------------
! Subroutine: samp_compute_mpi_c
! Purpose: compute sampling MPI distribution, halo C
!----------------------------------------------------------------------
subroutine samp_compute_mpi_c(samp,mpl,nam,geom)

implicit none

! Passed variables
class(samp_type),intent(inout) :: samp        ! Sampling
type(mpl_type),intent(inout) :: mpl           ! MPI data
type(nam_type),intent(in) :: nam              ! Namelist
type(geom_type),intent(in) :: geom            ! Geometry

! Local variables
integer :: jc3,ic0,ic0a,ic0c,ic1,ic1a,its,il0,d_n_s_max,d_n_s_max_loc,i_s,i_s_loc
integer,allocatable :: interpd_lg(:,:,:)
real(kind_real),allocatable :: displ_lon(:,:),displ_lat(:,:),lon_c1(:),lat_c1(:)
type(linop_type),allocatable :: dfull(:,:)

if (nam%displ_diag) then
   ! Allocation
   allocate(dfull(geom%nl0,nam%nts))
   allocate(displ_lon(geom%nc0,geom%nl0))
   allocate(displ_lat(geom%nc0,geom%nl0))
   allocate(lon_c1(nam%nc1))
   allocate(lat_c1(nam%nc1))

   ! Prepare displacement interpolation
   do its=1,nam%nts
      ! Local to global
      call mpl%loc_to_glb(geom%nl0,geom%nc0a,samp%displ_lon(:,:,its),geom%nc0,geom%c0_to_proc,geom%c0_to_c0a,.true.,displ_lon)
      call mpl%loc_to_glb(geom%nl0,geom%nc0a,samp%displ_lat(:,:,its),geom%nc0,geom%c0_to_proc,geom%c0_to_c0a,.true.,displ_lat)

      do il0=1,geom%nl0
         ! Copy Sc1 points
         do ic1=1,nam%nc1
            ic0 = samp%c1_to_c0(ic1)
            lon_c1(ic1) = displ_lon(ic0,il0)
            lat_c1(ic1) = displ_lat(ic0,il0)
         end do

         ! Compute interpolation
         call dfull(il0,its)%interp(mpl,geom%mesh,geom%kdtree,geom%nc0,geom%mask_c0(:,il0),nam%nc1,lon_c1,lat_c1, &
       & samp%c1l0_log(:,il0),nam%diag_interp)
      end do
   end do
end if

! Allocation
if (nam%displ_diag) then
   d_n_s_max = 0
   do its=1,nam%nts
      do il0=1,geom%nl0
         d_n_s_max = max(d_n_s_max,dfull(il0,its)%n_s)
      end do
   end do
   allocate(samp%lcheck_d(d_n_s_max,geom%nl0,nam%nts))
   allocate(samp%d(geom%nl0,nam%nts))
end if
allocate(samp%lcheck_c0c(geom%nc0))

! Halo C
samp%lcheck_c0c = samp%lcheck_c0a
do jc3=1,nam%nc3
   do ic1a=1,samp%nc1a
      ic1 = samp%c1a_to_c1(ic1a)
      if (any(samp%c1c3l0_log(ic1,jc3,:))) then
         ic0 = samp%c1c3_to_c0(ic1,jc3)
         samp%lcheck_c0c(ic0) = .true.
      end if
   end do
end do
if (nam%displ_diag) then
   samp%lcheck_d = .false.
   do its=1,nam%nts
      do il0=1,geom%nl0
         do i_s=1,dfull(il0,its)%n_s
            ic0 = dfull(il0,its)%col(i_s)
            ic1 = dfull(il0,its)%row(i_s)
            if (samp%lcheck_c1a(ic1)) then
               samp%lcheck_c0c(ic0) = .true.
               samp%lcheck_d(i_s,il0,its) = .true.
            end if
         end do
      end do
   end do
end if
samp%nc0c = count(samp%lcheck_c0c)
if (nam%displ_diag) then
   do its=1,nam%nts
      do il0=1,geom%nl0
         samp%d(il0,its)%n_s = count(samp%lcheck_d(:,il0,its))
      end do
   end do
end if

! Global <-> local conversions for fields

! Halo C
allocate(samp%c0c_to_c0(samp%nc0c))
allocate(samp%c0_to_c0c(geom%nc0))
call msi(samp%c0_to_c0c)
ic0c = 0
do ic0=1,geom%nc0
   if (samp%lcheck_c0c(ic0)) then
      ic0c = ic0c+1
      samp%c0c_to_c0(ic0c) = ic0
      samp%c0_to_c0c(ic0) = ic0c
   end if
end do

! Inter-halo conversions
allocate(samp%c0a_to_c0c(geom%nc0a))
do ic0a=1,geom%nc0a
   ic0 = geom%c0a_to_c0(ic0a)
   ic0c = samp%c0_to_c0c(ic0)
   samp%c0a_to_c0c(ic0a) = ic0c
end do

if (nam%displ_diag) then
   ! Global <-> local conversions for data
   d_n_s_max_loc = 0
   do its=1,nam%nts
      do il0=1,geom%nl0
         d_n_s_max_loc = max(d_n_s_max_loc,samp%d(il0,its)%n_s)
      end do
   end do
   allocate(interpd_lg(d_n_s_max_loc,geom%nl0,nam%nts))
   do its=1,nam%nts
      do il0=1,geom%nl0
         i_s_loc = 0
         do i_s=1,dfull(il0,its)%n_s
            if (samp%lcheck_d(i_s,il0,its)) then
               i_s_loc = i_s_loc+1
               interpd_lg(i_s_loc,il0,its) = i_s
            end if
         end do
      end do
   end do

   ! Local data
   do its=1,nam%nts
      do il0=1,geom%nl0
         write(samp%d(il0,its)%prefix,'(a,i3.3,a,i2.2)') 'd_',il0,'_',its
         samp%d(il0,its)%n_src = samp%nc0c
         samp%d(il0,its)%n_dst = samp%nc1a
         call samp%d(il0,its)%alloc
         do i_s_loc=1,samp%d(il0,its)%n_s
            i_s = interpd_lg(i_s_loc,il0,its)
            samp%d(il0,its)%row(i_s_loc) = samp%c1_to_c1a(dfull(il0,its)%row(i_s))
            samp%d(il0,its)%col(i_s_loc) = samp%c0_to_c0c(dfull(il0,its)%col(i_s))
            samp%d(il0,its)%S(i_s_loc) = dfull(il0,its)%S(i_s)
         end do
         call samp%d(il0,its)%reorder(mpl)
      end do
   end do
end if

! Setup communications
call samp%com_AC%setup(mpl,'com_AC',geom%nc0,geom%nc0a,samp%nc0c,samp%c0c_to_c0,samp%c0a_to_c0c,geom%c0_to_proc,geom%c0_to_c0a)

! Print results
write(mpl%info,'(a7,a,i4)') '','Parameters for processor #',mpl%myproc
write(mpl%info,'(a10,a,i8)') '','nc0c =      ',samp%nc0c
call flush(mpl%info)

end subroutine samp_compute_mpi_c

!----------------------------------------------------------------------
! Subroutine: samp_compute_mpi_f
! Purpose: compute sampling MPI distribution, halo F
!----------------------------------------------------------------------
subroutine samp_compute_mpi_f(samp,mpl,nam,geom)

implicit none

! Passed variables
class(samp_type),intent(inout) :: samp ! Sampling
type(mpl_type),intent(inout) :: mpl    ! MPI data
type(nam_type),intent(in) :: nam       ! Namelist
type(geom_type),intent(in) :: geom     ! Geometry

! Local variables
integer :: ic2,ic1,jc2,ic2a,ic2f,il0,kc2

! Allocation
allocate(samp%lcheck_c2f(nam%nc2))

! Halo definitions

! Halo F
do il0=1,geom%nl0
   samp%lcheck_c2f = samp%lcheck_c2a
   do ic2a=1,samp%nc2a
      ic2 = samp%c2a_to_c2(ic2a)
      ic1 = samp%c2_to_c1(ic2)
      jc2 = 1
      do while (samp%nn_c2_dist(jc2,ic2,min(il0,geom%nl0i))<nam%diag_rhflt)
          kc2 = samp%nn_c2_index(jc2,ic2,min(il0,geom%nl0i))
          samp%lcheck_c2f(kc2) = .true.
          jc2 = jc2+1
          if (jc2>nam%nc2) exit
      end do
   end do
end do
samp%nc2f = count(samp%lcheck_c2f)

! Global <-> local conversions for fields

! Halo F
allocate(samp%c2f_to_c2(samp%nc2f))
allocate(samp%c2_to_c2f(nam%nc2))
call msi(samp%c2_to_c2f)
ic2f = 0
do ic2=1,nam%nc2
   if (samp%lcheck_c2f(ic2)) then
      ic2f = ic2f+1
      samp%c2f_to_c2(ic2f) = ic2
      samp%c2_to_c2f(ic2) = ic2f
   end if
end do

! Inter-halo conversions
allocate(samp%c2a_to_c2f(samp%nc2a))
do ic2a=1,samp%nc2a
   ic2 = samp%c2a_to_c2(ic2a)
   ic2f = samp%c2_to_c2f(ic2)
   samp%c2a_to_c2f(ic2a) = ic2f
end do

! Setup communications
call samp%com_AF%setup(mpl,'com_AF',nam%nc2,samp%nc2a,samp%nc2f,samp%c2f_to_c2,samp%c2a_to_c2f,samp%c2_to_proc, &
 & samp%c2_to_c2a)

! Print results
write(mpl%info,'(a7,a,i4)') '','Parameters for processor #',mpl%myproc
write(mpl%info,'(a10,a,i8)') '','nc2f =       ',samp%nc2f
call flush(mpl%info)

end subroutine samp_compute_mpi_f

!----------------------------------------------------------------------
! Subroutine: samp_diag_filter
! Purpose: filter diagnostics
!----------------------------------------------------------------------
subroutine samp_diag_filter(samp,mpl,nam,geom,il0,filter_type,rflt,diag)

implicit none

! Passed variables
class(samp_type),intent(in) :: samp              ! Sampling
type(mpl_type),intent(inout) :: mpl              ! MPI data
type(nam_type),intent(in) :: nam                 ! Namelist
type(geom_type),intent(in) :: geom               ! Geometry
integer,intent(in) :: il0                        ! Level
character(len=*),intent(in) :: filter_type       ! Filter type
real(kind_real),intent(in) :: rflt               ! Filter support radius
real(kind_real),intent(inout) :: diag(samp%nc2a) ! Filtered diagnostic

! Local variables
integer :: ic2a,ic2,ic1,jc2,nc2eff,kc2,kc2_glb
integer,allocatable :: order(:)
real(kind_real) :: distnorm,norm,wgt
real(kind_real),allocatable :: diag_glb(:),diag_eff(:),diag_eff_dist(:)
logical :: nam_rad

! Check radius
nam_rad = .not.(abs(rflt-nam%diag_rhflt)>0.0)

if (rflt>0.0) then
   if (nam_rad) then
      ! Allocation
      allocate(diag_glb(samp%nc2f))

      ! Communication
      call samp%com_AF%ext(mpl,diag,diag_glb)
   else
      ! Allocation
      allocate(diag_glb(nam%nc2))

      ! Local to global
     call mpl%loc_to_glb(samp%nc2a,diag,nam%nc2,samp%c2_to_proc,samp%c2_to_c2a,.true.,diag_glb)
   end if

   !$omp parallel do schedule(static) private(ic2a,ic2,ic1,nc2eff,jc2,distnorm,norm,wgt) firstprivate(diag_eff,diag_eff_dist,order)
   do ic2a=1,samp%nc2a
      ! Indices
      ic2 = samp%c2a_to_c2(ic2a)
      ic1 = samp%c2_to_c1(ic2)

      ! Allocation
      allocate(diag_eff(nam%nc2))
      allocate(diag_eff_dist(nam%nc2))

      ! Build diag_eff of valid points
      nc2eff = 0
      jc2 = 1
      do while (samp%nn_c2_dist(jc2,ic2,min(il0,geom%nl0i))<rflt)
         ! Check the point validity
         kc2 = samp%nn_c2_index(jc2,ic2,min(il0,geom%nl0i))
         if (nam_rad) then
            kc2_glb = samp%c2_to_c2f(kc2)
         else
            kc2_glb = kc2
         end if
         if (isnotmsr(diag_glb(kc2_glb))) then
            nc2eff = nc2eff+1
            diag_eff(nc2eff) = diag_glb(kc2_glb)
            diag_eff_dist(nc2eff) = samp%nn_c2_dist(jc2,ic2,min(il0,geom%nl0i))
         end if
         jc2 = jc2+1
         if (jc2>nam%nc2) exit
      end do

      ! Apply filter
      if (nc2eff>0) then
         select case (trim(filter_type))
         case ('average')
            ! Compute average
            diag(ic2a) = sum(diag_eff(1:nc2eff))/real(nc2eff,kind_real)
         case ('gc99')
            ! Gaspari-Cohn (1999) kernel
            diag(ic2a) = 0.0
            norm = 0.0
            do jc2=1,nc2eff
               distnorm = diag_eff_dist(jc2)/rflt
               wgt = gc99(mpl,distnorm)
               diag(ic2a) = diag(ic2a)+wgt*diag_eff(jc2)
               norm = norm+wgt
            end do
            if (norm>0.0) diag(ic2a) = diag(ic2a)/norm
         case ('median')
            ! Compute median
            allocate(order(nc2eff))
            call qsort(nc2eff,diag_eff(1:nc2eff),order)
            if (mod(nc2eff,2)==0) then
               diag(ic2a) = 0.5*(diag_eff(nc2eff/2)+diag_eff(nc2eff/2+1))
            else
               diag(ic2a) = diag_eff((nc2eff+1)/2)
            end if
            deallocate(order)
         case default
            ! Wrong filter
            call mpl%abort('wrong filter type')
         end select
      else
         call msr(diag(ic2a))
      end if

      ! Release memory
      deallocate(diag_eff)
      deallocate(diag_eff_dist)
   end do
   !$omp end parallel do
end if

end subroutine samp_diag_filter

!----------------------------------------------------------------------
! Subroutine: samp_diag_fill
! Purpose: fill diagnostics missing values
!----------------------------------------------------------------------
subroutine samp_diag_fill(samp,mpl,nam,geom,il0,diag)

implicit none

! Passed variables
class(samp_type),intent(in) :: samp              ! Sampling
type(mpl_type),intent(inout) :: mpl              ! MPI data
type(nam_type),intent(in) :: nam                 ! Namelist
type(geom_type),intent(in) :: geom               ! Geometry
integer,intent(in) :: il0                        ! Level
real(kind_real),intent(inout) :: diag(samp%nc2a) ! Filtered diagnostic

! Local variables
integer :: nmsr,nmsr_tot,ic2,jc2,kc2
real(kind_real),allocatable :: diag_glb(:)

! Count missing points
if (samp%nc2a>0) then
   nmsr = count(ismsr(diag))
else
   nmsr = 0
end if
call mpl%f_comm%allreduce(nmsr,nmsr_tot,fckit_mpi_sum())

if (nmsr_tot>0) then
   ! Allocation
   allocate(diag_glb(nam%nc2))

   ! Local to global
   call mpl%loc_to_glb(samp%nc2a,diag,nam%nc2,samp%c2_to_proc,samp%c2_to_c2a,.false.,diag_glb)

   if (mpl%main) then
      do ic2=1,nam%nc2
         jc2 = 1
         do while (ismsr(diag_glb(ic2)))
            kc2 = samp%nn_c2_index(jc2,ic2,min(il0,geom%nl0i))
            if (isnotmsr(diag_glb(kc2))) diag_glb(ic2) = diag_glb(kc2)
            jc2 = jc2+1
            if (jc2>nam%nc2) exit
         end do
      end do
   end if

   ! Global to local
   call mpl%glb_to_loc(nam%nc2,samp%c2_to_proc,samp%c2_to_c2a,diag_glb,samp%nc2a,diag)

   ! Release memory
   deallocate(diag_glb)
end if

end subroutine samp_diag_fill

end module type_samp