tools_func.F90

Go to the documentation of this file.

# 1 "/Users/miesch/JEDI/code/working_copy/public/fv3-bundle/saber/src/saber/bump/tools_func.fypp"
# 1 "/Users/miesch/JEDI/code/working_copy/public/fv3-bundle/saber/src/saber/bump/../instrumentation.fypp" 1
# 1 "/Users/miesch/JEDI/code/working_copy/public/fv3-bundle/saber/src/saber/bump/../subr_list.fypp" 1
!----------------------------------------------------------------------
! Header: subr_list
!> Subroutines/functions list
! Author: Benjamin Menetrier
! Licensing: this code is distributed under the CeCILL-C license
! Copyright 2015-... UCAR, CERFACS, METEO-FRANCE and IRIT
!----------------------------------------------------------------------
 
# 726 "/Users/miesch/JEDI/code/working_copy/public/fv3-bundle/saber/src/saber/bump/../subr_list.fypp"
# 2 "/Users/miesch/JEDI/code/working_copy/public/fv3-bundle/saber/src/saber/bump/../instrumentation.fypp" 2
!----------------------------------------------------------------------
! Header: instrumentation
!> Instrumentation functions
! Author: Benjamin Menetrier
! Licensing: this code is distributed under the CeCILL-C license
! Copyright 2015-... UCAR, CERFACS, METEO-FRANCE and IRIT
!----------------------------------------------------------------------
 
# 112 "/Users/miesch/JEDI/code/working_copy/public/fv3-bundle/saber/src/saber/bump/../instrumentation.fypp"
# 2 "/Users/miesch/JEDI/code/working_copy/public/fv3-bundle/saber/src/saber/bump/tools_func.fypp" 2
!----------------------------------------------------------------------
! Module: tools_func
!> Usual functions
! Author: Benjamin Menetrier
! Licensing: this code is distributed under the CeCILL-C license
! Copyright 2015-... UCAR, CERFACS, METEO-FRANCE and IRIT
!----------------------------------------------------------------------
module tools_func
 
use atlas_module, only: atlas_geometry
use iso_c_binding, only: c_int16_t,c_int32_t
use tools_asa007, only: cholesky,syminv
use tools_const, only: zero,tenth,hundredth,half,quarter,one,two,three,four,five,eight,ten,thousand,pi,deg2rad,rad2deg
use tools_kinds, only: kind_short,kind_real
use tools_qsort, only: qsort
use tools_repro, only: inf,sup,infeq,small
use type_mpl, only: mpl_type
 
 
implicit none
 
real(kind_real),parameter :: gc2gau = 0.28_kind_real  !< GC99 support radius to Gaussian Daley length-scale (empirical)
real(kind_real),parameter :: gau2gc = 3.57_kind_real  !< Gaussian Daley length-scale to GC99 support radius (empirical)
real(kind_real),parameter :: dmin = 1.0e-12_kind_real !< Minimum tensor diagonal value
real(kind_real),parameter :: condmax = thousand       !< Maximum tensor conditioning number
integer,parameter :: m = 0                            !< Number of implicit iteration for the Matern function (0: Gaussian)
 
interface
   function c_fletcher32(n,var) bind(c,name='fletcher32') result(hash)
   use iso_c_binding, only: c_int16_t,c_int32_t
   integer(c_int32_t) :: n
   integer(c_int16_t) :: var(*)
   integer(c_int32_t) :: hash
   end function c_fletcher32
end interface
interface fletcher32
   module procedure func_fletcher32
end interface
interface lonlatmod
   module procedure func_lonlatmod
end interface
interface gridhash
   module procedure func_gridhash
end interface
interface independent_levels
   module procedure func_independent_levels
end interface
interface lonlathash
   module procedure func_lonlathash
end interface
interface sphere_dist
   module procedure func_sphere_dist
end interface
interface lonlat2xyz
   module procedure func_lonlat2xyz
end interface
interface xyz2lonlat
   module procedure func_xyz2lonlat
end interface
interface vector_product
   module procedure func_vector_product
end interface
interface det
   module procedure func_det
end interface
interface order_cc
   module procedure func_order_cc
end interface
interface add
   module procedure func_add
end interface
interface divide
   module procedure func_divide
end interface
interface fit_diag
   module procedure func_fit_diag
end interface
interface gc99
   module procedure func_gc99
end interface
interface fit_func
   module procedure func_fit_func
end interface
interface fit_lct
   module procedure func_fit_lct
end interface
interface lct_d2h
   module procedure func_lct_d2h
end interface
interface lct_h2r
   module procedure func_lct_h2r
end interface
interface lct_r2d
   module procedure func_lct_r2d
end interface
interface check_cond
   module procedure func_check_cond
end interface
interface matern
   module procedure func_matern
end interface
interface cholesky
   module procedure func_cholesky
end interface
interface syminv
   module procedure func_syminv
end interface
interface histogram
   module procedure func_histogram
end interface
interface cx_to_cxa
   module procedure func_cx_to_cxa
end interface
interface cx_to_proc
   module procedure func_cx_to_proc
end interface
interface cx_to_cxu
   module procedure func_cx_to_cxu
end interface
 
private
public :: gc2gau,gau2gc,dmin,m
public :: fletcher32,lonlatmod,gridhash,independent_levels,lonlathash,sphere_dist,lonlat2xyz,xyz2lonlat,vector_product,det, &
 & order_cc,add,divide,fit_diag,fit_func,fit_lct,lct_d2h,lct_h2r,lct_r2d,check_cond,cholesky,syminv,histogram, &
 & cx_to_cxa,cx_to_proc,cx_to_cxu
 
contains
 
!----------------------------------------------------------------------
! Function: func_fletcher32
!> Fletcher-32 checksum algorithm
!----------------------------------------------------------------------
function func_fletcher32(var) result(value)
 
implicit none
 
! Passed variables
real(kind_real),intent(in) :: var(:) !< Variable
 
! Returned variable
integer :: value
 
! Set name
 
 
! Probe in
 
 
! Call C function
value = c_fletcher32(size(transfer(var,(/0_kind_short/))),transfer(var,(/0_kind_short/)))
 
! Probe out
 
 
end function func_fletcher32
 
!----------------------------------------------------------------------
! Subroutine: func_lonlatmod
!> Set latitude between -pi/2 and pi/2 and longitude between -pi and pi
!----------------------------------------------------------------------
subroutine func_lonlatmod(lon,lat)
 
implicit none
 
! Passed variables
real(kind_real),intent(inout) :: lon !< Longitude (radians)
real(kind_real),intent(inout) :: lat !< Latitude (radians)
 
! Set name
 
 
! Probe in
 
 
! Check latitude bounds
if (lat>half*pi) then
   lat = pi-lat
   lon = lon+pi
elseif (lat<-half*pi) then
   lat = -pi-lat
   lon = lon+pi
end if
 
! Check longitude bounds
if (lon>pi) then
   lon = lon-two*pi
elseif (lon<-pi) then
   lon = lon+two*pi
end if
 
! Same zero longitude for poles
if (abs(lat)>(half-1.0e-6_kind_real)*pi) lon = zero
 
! Probe out
 
 
end subroutine func_lonlatmod
 
!----------------------------------------------------------------------
! Subroutine: func_gridhash
!> Compute grid hash profile
!----------------------------------------------------------------------
subroutine func_gridhash(ncx,nlx,lon_cx,lat_cx,mask_cx,grid_hash)
 
implicit none
 
! Passed variables
integer,intent(in) :: ncx                 !< Number of points
integer,intent(in) :: nlx                 !< Number of levels
real(kind_real),intent(in) :: lon_cx(ncx) !< Longitude (radians)
real(kind_real),intent(in) :: lat_cx(ncx) !< Latitude (radians)
logical,intent(in) :: mask_cx(ncx,nlx)    !< Mask
integer,intent(out) :: grid_hash(0:nlx)   !< Grid hash profile
 
! Local variables
integer :: ilx,ncx_eff,icx_eff,icx
real(kind_real),allocatable :: lonlat(:)
 
! Set name
 
 
! Probe in
 
 
do ilx=1,nlx
   ! Count points in the mask
   ncx_eff = count(mask_cx(:,ilx))
 
   if (ncx_eff>0) then
      ! Allocation
      allocate(lonlat(2*ncx_eff))
 
      ! Copy lonlat
      icx_eff = 0
      do icx=1,ncx
         if (mask_cx(icx,ilx)) then
            icx_eff = icx_eff+1
            lonlat(icx_eff) = lon_cx(icx)
            icx_eff = icx_eff+1
            lonlat(icx_eff) = lat_cx(icx)
         end if
      end do
 
      ! Compute hash
      grid_hash(ilx) = fletcher32(lonlat)
 
      ! Release memory
      deallocate(lonlat)
   else
      grid_hash(ilx) = 0
   end if
end do
 
! Final grid hash
grid_hash(0) = fletcher32(real(grid_hash(1:nlx),kind_real))
 
! Probe out
 
 
end subroutine func_gridhash
 
!----------------------------------------------------------------------
! Subroutine: func_independent_levels
!> Compute independent levels
!----------------------------------------------------------------------
subroutine func_independent_levels(mpl,nlx,grid_hash,nlxi,lx_to_lxi,lxi_to_lx,ifmt)
 
implicit none
 
! Passed variables
type(mpl_type),intent(inout) :: mpl   !< MPI data
integer,intent(in) :: nlx             !< Number of levels
integer,intent(in) :: grid_hash(nlx)  !< Grid hash profile
integer,intent(out) :: nlxi           !< Number of independent levels
integer,intent(out) :: lx_to_lxi(nlx) !< Levels to independent levels
integer,intent(out) :: lxi_to_lx(nlx) !< Independent levels to levels
integer,intent(in) :: ifmt            !< Indentation
 
! Local variables
integer :: ilx,ilxi,jlx,jlxi,n
integer :: proc_to_grid_hash(mpl%nproc),grid_hash_glb(nlx)
character(len=1024) :: cfmt
 
! Set name
 
 
! Probe in
 
 
! Compute global hash value
do ilx=1,nlx
   call mpl%f_comm%allgather(grid_hash(ilx),proc_to_grid_hash)
   grid_hash_glb(ilx) = fletcher32(real(proc_to_grid_hash,kind_real))
end do
 
! Count independent levels
nlxi = 1
do ilx=2,nlx
   if (all(grid_hash_glb(1:ilx-1)/=grid_hash_glb(ilx))) nlxi = nlxi+1
end do
 
! Initialization
lxi_to_lx = mpl%msv%vali
 
! Get independent level
ilx = 1
ilxi = 1
lx_to_lxi(ilx) = ilxi
lxi_to_lx(ilxi) = ilx
do ilx=2,nlx
   if (all(grid_hash_glb(1:ilx-1)/=grid_hash_glb(ilx))) then
      ! New independent level
      ilxi = ilxi+1
      lx_to_lxi(ilx) = ilxi
      lxi_to_lx(ilxi) = ilx
   else
      ! Similar level
      do jlx=1,ilx-1
         if (grid_hash_glb(jlx)==grid_hash_glb(ilx)) then
            jlxi = lx_to_lxi(jlx)
            lx_to_lxi(ilx) = jlxi
            exit
         end if
      end do
   end if
end do
 
! Print levels
write(cfmt,'(a,i2.2,a)') '(a',ifmt,',a)'
write(mpl%info,trim(cfmt)) '','Compute independent levels: '
call mpl%flush(.false.)
do ilxi=1,nlxi
   ilx = lxi_to_lx(ilxi)
   n = count(lx_to_lxi==ilxi)
   if (n<10) then
      cfmt = '(i3,a,i1,a)'
   elseif (n<100) then
      cfmt = '(i3,a,i2,a)'
   else
      cfmt = '(i3,a,i3,a)'
   end if
   write(mpl%info,trim(cfmt)) ilx,'[',n,'] '
   call mpl%flush(.false.)
end do
write(mpl%info,'(a)') ''
call mpl%flush
 
! Probe out
 
 
end subroutine func_independent_levels
 
!----------------------------------------------------------------------
! Function: func_lonlathash
!> Define a unique real from a lon/lat pair
!----------------------------------------------------------------------
function func_lonlathash(lon,lat,il) result(value)
 
implicit none
 
! Passed variables
real(kind_real),intent(in) :: lon !< Longitude (radians)
real(kind_real),intent(in) :: lat !< Latitude (radians)
integer,intent(in),optional :: il !< Level
 
! Returned variable
real(kind_real) :: value
 
! Local variables
real(kind_real) :: lontmp,lattmp
 
! Set name
 
 
! Probe in
 
 
! Set correct lon/lat
lontmp = lon
lattmp = lat
call lonlatmod(lontmp,lattmp)
 
! Hash value
value = aint((lontmp+pi)*1.0e7_kind_real)+(lattmp+half*pi)*tenth
if (present(il)) value = value+real(il*1e8,kind_real)
 
! Probe out
 
 
end function func_lonlathash
 
!----------------------------------------------------------------------
! Subroutine: func_sphere_dist
!> Compute the great-circle distance between two points
!----------------------------------------------------------------------
subroutine func_sphere_dist(lon_i,lat_i,lon_f,lat_f,dist)
 
implicit none
 
! Passed variable
real(kind_real),intent(in) :: lon_i !< Initial point longitude (radians)
real(kind_real),intent(in) :: lat_i !< Initial point latitude (radians)
real(kind_real),intent(in) :: lon_f !< Final point longitude (radians)
real(kind_real),intent(in) :: lat_f !< Final point latitude (radians)
real(kind_real),intent(out) :: dist !< Great-circle distance
 
! Local variables
type(atlas_geometry) :: ageometry
 
! Set name
 
 
! Probe in
 
 
! Create ATLAS geometry
ageometry = atlas_geometry('UnitSphere')
 
! Compute distance
dist = ageometry%distance(lon_i*rad2deg,lat_i*rad2deg,lon_f*rad2deg,lat_f*rad2deg)
 
! Probe out
 
 
end subroutine func_sphere_dist
 
!----------------------------------------------------------------------
! Subroutine: func_lonlat2xyz
!> Convert longitude/latitude to cartesian coordinates
!----------------------------------------------------------------------
subroutine func_lonlat2xyz(mpl,lon,lat,x,y,z)
 
implicit none
 
! Passed variables
type(mpl_type),intent(inout) :: mpl !< MPI data
real(kind_real),intent(in) :: lon   !< Longitude (radians)
real(kind_real),intent(in) :: lat   !< Latitude (radians)
real(kind_real),intent(out) :: x    !< X coordinate
real(kind_real),intent(out) :: y    !< Y coordinate
real(kind_real),intent(out) :: z    !< Z coordinate
 
! Local variables
type(atlas_geometry) :: ageometry
 
! Set name
 
 
! Probe in
 
 
if (mpl%msv%isnot(lat).and.mpl%msv%isnot(lon)) then
   ! Check longitude/latitude
   if (inf(lon,-pi).and.sup(lon,pi)) call mpl%abort('func_lonlat2xyz','wrong longitude')
   if (inf(lat,-half*pi).and.sup(lat,-half*pi)) call mpl%abort('func_lonlat2xyz','wrong latitude')
 
   ! Create ATLAS geometry
   ageometry = atlas_geometry('UnitSphere')
 
   ! Convert to x/y/z
   call ageometry%lonlat2xyz(lon*rad2deg,lat*rad2deg,x,y,z)
else
   ! Missing values
   x = mpl%msv%valr
   y = mpl%msv%valr
   z = mpl%msv%valr
end if
 
! Probe out
 
 
end subroutine func_lonlat2xyz
 
!----------------------------------------------------------------------
! Subroutine: func_xyz2lonlat
!> Convert longitude/latitude to cartesian coordinates
!----------------------------------------------------------------------
subroutine func_xyz2lonlat(mpl,x,y,z,lon,lat)
 
implicit none
 
! Passed variables
type(mpl_type),intent(in) :: mpl   !< MPI data
real(kind_real),intent(in) :: x    !< X coordinate
real(kind_real),intent(in) :: y    !< Y coordinate
real(kind_real),intent(in) :: z    !< Z coordinate
real(kind_real),intent(out) :: lon !< Longitude (radians)
real(kind_real),intent(out) :: lat !< Latitude (radians)
 
! Local variables
type(atlas_geometry) :: ageometry
 
! Set name
 
 
! Probe in
 
 
if (mpl%msv%isnot(x).and.mpl%msv%isnot(y).and.mpl%msv%isnot(z)) then
   ! Create ATLAS geometry
   ageometry = atlas_geometry('UnitSphere')
 
   ! Convert to lon/lat
   call ageometry%xyz2lonlat(x,y,z,lon,lat)
 
   ! Copy coordinates
   lon = lon*deg2rad
   lat = lat*deg2rad
else
   ! Missing values
   lon = mpl%msv%valr
   lat = mpl%msv%valr
end if
 
! Probe out
 
 
end subroutine func_xyz2lonlat
 
!----------------------------------------------------------------------
! Subroutine: func_vector_product
!> Compute normalized vector product
!----------------------------------------------------------------------
subroutine func_vector_product(v1,v2,vp)
 
implicit none
 
! Passed variables
real(kind_real),intent(in) :: v1(3)  !< First vector
real(kind_real),intent(in) :: v2(3)  !< Second vector
real(kind_real),intent(out) :: vp(3) !< Vector product
 
! Local variable
real(kind_real) :: r
 
! Set name
 
 
! Probe in
 
 
! Vector product
vp(1) = v1(2)*v2(3)-v1(3)*v2(2)
vp(2) = v1(3)*v2(1)-v1(1)*v2(3)
vp(3) = v1(1)*v2(2)-v1(2)*v2(1)
 
! Normalization
r = sqrt(sum(vp**2))
if (r>zero) vp = vp/r
 
! Probe out
 
 
end subroutine func_vector_product
 
!----------------------------------------------------------------------
! Subroutine: func_det
!> Compute determinant (vector triple product)
!----------------------------------------------------------------------
subroutine func_det(v1,v2,v3,p,cflag)
 
implicit none
 
! Passed variables
real(kind_real),intent(in) :: v1(3) !< First vector
real(kind_real),intent(in) :: v2(3) !< Second vector
real(kind_real),intent(in) :: v3(3) !< Third vector
real(kind_real),intent(out) :: p    !< Determinant
logical,intent(out) :: cflag        !< Confidence flag
 
! Local variable
integer :: i
real(kind_real) :: terms(6)
 
! Set name
 
 
! Probe in
 
 
! Terms
terms(1) = v1(2)*v2(3)*v3(1)
terms(2) = -v1(3)*v2(2)*v3(1)
terms(3) = v1(3)*v2(1)*v3(2)
terms(4) = -v1(1)*v2(3)*v3(2)
terms(5) = v1(1)*v2(2)*v3(3)
terms(6) = -v1(2)*v2(1)*v3(3)
 
! Sum
p = sum(terms)
 
! Confidence flag
cflag = .true.
do i=1,6
   if ((abs(terms(i))>zero).and.small(p,terms(i))) cflag = .false.
end do
 
! Probe out
 
 
end subroutine func_det
 
!----------------------------------------------------------------------
! Subroutine: func_order_cc
!> Order points in counter-clockwise order with respect to a central point
!----------------------------------------------------------------------
subroutine func_order_cc(mpl,lon,lat,n,x,y,z,order,diff)
 
implicit none
 
! Passed variables
type(mpl_type),intent(inout) :: mpl             !< MPI data
real(kind_real),intent(in) :: lon               !< Longitude of the central point
real(kind_real),intent(in) :: lat               !< Latitude of the central point
integer :: n                                    !< Number of points
real(kind_real),intent(in) :: x(n)              !< List of X-coordinates
real(kind_real),intent(in) :: y(n)              !< List of Y-coordinates
real(kind_real),intent(in) :: z(n)              !< List of Z-coordinates
integer,intent(out) :: order(n)                 !< Counter-clockwise order
real(kind_real),intent(out),optional :: diff(n) !< Angles differences
 
! Local variable
integer :: i
real(kind_real) :: rvec(3),costheta,sintheta,p(3),rvecxv(3),v(3),list(n)
real(kind_real),allocatable :: list_save(:)
 
! Set name
 
 
! Probe in
 
 
! Rotation vector in cartesian coordinates
call lonlat2xyz(mpl,lon-half*pi,zero,rvec(1),rvec(2),rvec(3))
 
! Rotation angle
costheta = cos(half*pi-lat)
sintheta = sin(half*pi-lat)
 
! Compute angle
do i=1,n
   ! Rodrigues' rotation
   p = (/x(i),y(i),z(i)/)
   call vector_product(rvec,p,rvecxv)
   v = p*costheta+rvecxv*sintheta+rvec*sum(rvec*p)*(one-costheta)
 
   ! Angle
   list(i) = atan2(v(2),v(1))
end do
 
if (present(diff)) then
   ! Allocation
   allocate(list_save(n))
 
   ! Copy
   list_save = list
end if
 
! Sort angles in counter-clockwise order
call qsort(n,list,order)
 
if (present(diff)) then
   ! Get angles differences
   diff(order(1)) = list_save(order(1))-list_save(order(n))+two*pi
   do i=2,n
      diff(order(i)) = list_save(order(i))-list_save(order(i-1))
   end do
 
   ! Release memory
   deallocate(list_save)
end if
 
! Probe out
 
 
end subroutine func_order_cc
 
!----------------------------------------------------------------------
! Subroutine: func_add
!> Check if value missing and add if not missing
!----------------------------------------------------------------------
subroutine func_add(mpl,val,cumul,num,wgt)
 
implicit none
 
! Passed variables
type(mpl_type),intent(in) :: mpl           !< MPI data
real(kind_real),intent(in) :: val          !< Value to add
real(kind_real),intent(inout) :: cumul     !< Cumul
real(kind_real),intent(inout) :: num       !< Number of values
real(kind_real),intent(in),optional :: wgt !< Weight
 
! Local variables
real(kind_real) :: lwgt
 
! Set name
 
 
! Probe in
 
 
! Initialize weight
lwgt = one
if (present(wgt)) lwgt = wgt
 
! Add value to cumul
if (mpl%msv%isnot(val)) then
   cumul = cumul+lwgt*val
   num = num+lwgt
end if
 
! Probe out
 
 
end subroutine func_add
 
!----------------------------------------------------------------------
! Subroutine: func_divide
!> Check if value missing and divide if not missing
!----------------------------------------------------------------------
subroutine func_divide(mpl,val,num)
 
implicit none
 
! Passed variables
type(mpl_type),intent(in) :: mpl     !< MPI data
real(kind_real),intent(inout) :: val !< Value to divide
real(kind_real),intent(in) :: num    !< Divider
 
! Set name
 
 
! Probe in
 
 
! Divide cumul by num
if (abs(num)>zero) then
   val = val/num
else
   val = mpl%msv%valr
end if
 
! Probe out
 
 
end subroutine func_divide
 
!----------------------------------------------------------------------
! Subroutine: func_fit_diag
!> Compute diagnostic fit function
!----------------------------------------------------------------------
subroutine func_fit_diag(mpl,nc3,nl0r,nl0,l0rl0_to_l0,disth,distv,coef,rh,rv,fit)
 
implicit none
 
! Passed variables
type(mpl_type),intent(inout) :: mpl              !< MPI data
integer,intent(in) :: nc3                        !< Number of classes
integer,intent(in) :: nl0r                       !< Effective number of levels
integer,intent(in) :: nl0                        !< Number of levels
integer,intent(in) :: l0rl0_to_l0(nl0r,nl0)      !< Effective level to level
real(kind_real),intent(in) :: disth(nc3)         !< Horizontal distance
real(kind_real),intent(in) :: distv(nl0,nl0)     !< Vertical distance
real(kind_real),intent(in) :: coef(nl0)          !< Diagonal coefficient
real(kind_real),intent(in) :: rh(nl0)            !< Horizontal support radius
real(kind_real),intent(in) :: rv(nl0)            !< Vertical support radius
real(kind_real),intent(out) :: fit(nc3,nl0r,nl0) !< Fit
 
! Local variables
integer :: jl0r,jl0,djl0,il0,kl0r,kl0,ic3,jc3,djc3,kc3,ip,jp,np,np_new,nc3max
integer :: plist(nc3*nl0r,2),plist_new(nc3*nl0r,2)
real(kind_real) :: rhsq,rvsq,distvsq,disttest,rhmax
real(kind_real) :: predistnorm(-1:1,nc3,-1:1,nl0),distnorm(nc3,nl0r)
logical :: add_to_front
 
! Set name
 
 
! Probe in
 
 
! Initialization
fit = zero
 
! Find maximum class
rhmax = maxval(rh)
nc3max = 0
do ic3=1,nc3
   if (disth(ic3)>rhmax) exit
   nc3max = ic3
end do
 
! Precompute normalized local distances
predistnorm = one
do il0=1,nl0
   do djl0=-1,1
      ! Level index
      jl0 = max(1,min(il0+djl0,nl0))
 
      ! Averaged fit parameters
      if (mpl%msv%isnot(rh(il0)).and.mpl%msv%isnot(rh(jl0))) then
         rhsq = half*(rh(il0)**2+rh(jl0)**2)
      else
         rhsq = zero
      end if
      if (mpl%msv%isnot(rv(il0)).and.mpl%msv%isnot(rv(jl0))) then
         rvsq = half*(rv(il0)**2+rv(jl0)**2)
      else
         rvsq = zero
      end if
 
      ! Squared vertical distance
      if (rvsq>zero) then
         distvsq = distv(jl0,il0)**2/rvsq
      elseif (il0/=jl0) then
         distvsq = one
      else
         distvsq = zero
      end if
 
      do ic3=1,nc3max
         do djc3=-1,1
            ! Horizontal index
            jc3 = max(1,min(ic3+djc3,nc3max))
 
            ! Squared normalized distance initialization (squared vertical component)
            predistnorm(djc3,ic3,djl0,il0) = distvsq
 
            ! Add squared horizontal component to squared normalized distance
            if (rhsq>zero) then
               predistnorm(djc3,ic3,djl0,il0) = predistnorm(djc3,ic3,djl0,il0)+(disth(jc3)-disth(ic3))**2/rhsq
            elseif (ic3/=jc3) then
               predistnorm(djc3,ic3,djl0,il0) = predistnorm(djc3,ic3,djl0,il0)+one
            end if
 
            ! Get normalized distance
            predistnorm(djc3,ic3,djl0,il0) = sqrt(predistnorm(djc3,ic3,djl0,il0))
         end do
      end do
   end do
end do
 
! Compte fit
do il0=1,nl0
   ! Initialize the front
   np = 1
   plist = mpl%msv%vali
   plist(1,1) = 1
   do jl0r=1,nl0r
      if (l0rl0_to_l0(jl0r,il0)==il0) plist(1,2) = jl0r
   end do
   distnorm = one
   distnorm(plist(1,1),plist(1,2)) = zero
 
   do while (np>0)
      ! Propagate the front
      np_new = 0
 
      do ip=1,np
         ! Indices of the central point
         jc3 = plist(ip,1)
         jl0r = plist(ip,2)
         jl0 = l0rl0_to_l0(jl0r,il0)
 
         ! Loop over neighbors
         do kc3=max(jc3-1,1),min(jc3+1,nc3max)
            do kl0r=max(jl0r-1,1),min(jl0r+1,nl0r)
               kl0 = l0rl0_to_l0(kl0r,il0)
               disttest = distnorm(jc3,jl0r)+predistnorm(kc3-jc3,jc3,kl0-jl0,jl0)
 
               if (disttest<one) then
                  ! Point is inside the support
                  if (disttest<distnorm(kc3,kl0r)) then
                     ! Update distance
                     distnorm(kc3,kl0r) = disttest
 
                     ! Check if the point should be added to the front (avoid duplicates)
                     add_to_front = .true.
                     do jp=1,np_new
                        if ((plist_new(jp,1)==kc3).and.(plist_new(jp,2)==kl0r)) then
                           add_to_front = .false.
                           exit
                        end if
                     end do
 
                     if (add_to_front) then
                        ! Add point to the front
                        np_new = np_new+1
                        plist_new(np_new,1) = kc3
                        plist_new(np_new,2) = kl0r
                     end if
                  end if
               end if
            end do
         end do
      end do
 
      ! Copy new front
      np = np_new
      plist(1:np,:) = plist_new(1:np,:)
   end do
 
   do jl0r=1,nl0r
      ! Level index
      jl0 = l0rl0_to_l0(jl0r,il0)
 
      do jc3=1,nc3
         ! Fit function
         fit(jc3,jl0r,il0) = fit_func(mpl,distnorm(jc3,jl0r))
      end do
   end do
end do
 
! Diagonal coefficient
do il0=1,nl0
   do jl0r=1,nl0r
      jl0 = l0rl0_to_l0(jl0r,il0)
      fit(:,jl0r,il0) = fit(:,jl0r,il0)*sqrt(coef(il0)*coef(jl0))
   end do
end do
 
! Set to missing values if no value available
do il0=1,nl0
   if (mpl%msv%is(coef(il0)).or.mpl%msv%is(rh(il0)).or.mpl%msv%is(rv(il0))) fit(:,:,il0) = mpl%msv%valr
   do jl0r=1,nl0r
      jl0 = l0rl0_to_l0(jl0r,il0)
      if (mpl%msv%is(coef(jl0)).or.mpl%msv%is(rh(jl0)).or.mpl%msv%is(rv(jl0))) fit(:,jl0r,il0) = mpl%msv%valr
   end do
end do
 
! Probe out
 
 
end subroutine func_fit_diag
 
!----------------------------------------------------------------------
! Function: func_gc99
!> Gaspari and Cohn (1999) function, with the support radius as a parameter
!----------------------------------------------------------------------
function func_gc99(distnorm) result(value)
 
! Passed variables
real(kind_real),intent(in) :: distnorm !< Normalized distance
 
! Returned variable
real(kind_real) :: value
 
! Set name
 
 
! Probe in
 
 
! Gaspari and Cohn (1999) function
if (distnorm<half) then
   value = one-distnorm
   value = one+eight/five*distnorm*value
   value = one-three/four*distnorm*value
   value = one-20.0_kind_real/three*distnorm**2*value
else if (distnorm<one) then
   value = one-distnorm/three
   value = one-eight/five*distnorm*value
   value = one+three/four*distnorm*value
   value = one-two/three*distnorm*value
   value = one-five/two*distnorm*value
   value = one-12.0_kind_real*distnorm*value
   value = -value/(three*distnorm)
else
   value = zero
end if
 
! Probe out
 
 
end function func_gc99
 
!----------------------------------------------------------------------
! Function: func_fit_func
!> Fit_function
!----------------------------------------------------------------------
function func_fit_func(mpl,distnorm) result(value)
 
! Passed variables
type(mpl_type),intent(inout) :: mpl    !< MPI data
real(kind_real),intent(in) :: distnorm !< Normalized distance
 
! Returned variable
real(kind_real) :: value
 
! Set name
 
 
! Probe in
 
 
! Distance check bound
if (distnorm<zero) call mpl%abort('func_fit_func','negative normalized distance')
 
! Gaspari and Cohn (1999) function
value = gc99(distnorm)
 
! Enforce positivity
value = max(value,zero)
 
! Probe out
 
 
end function func_fit_func
 
!----------------------------------------------------------------------
! Subroutine: func_fit_lct
!> LCT fit
!----------------------------------------------------------------------
subroutine func_fit_lct(mpl,nc,nl0,dxsq,dysq,dxdy,dzsq,dmask,nscales,D,coef,fit)
 
implicit none
 
! Passed variables
type(mpl_type),intent(inout) :: mpl         !< MPI data
integer,intent(in) :: nc                    !< Number of classes
integer,intent(in) :: nl0                   !< Number of levels
real(kind_real),intent(in) :: dxsq(nc,nl0)  !< Zonal separation squared
real(kind_real),intent(in) :: dysq(nc,nl0)  !< Meridian separation squared
real(kind_real),intent(in) :: dxdy(nc,nl0)  !< Zonal x meridian separations product
real(kind_real),intent(in) :: dzsq(nc,nl0)  !< Vertical separation squared
logical,intent(in) :: dmask(nc,nl0)         !< Mask
integer,intent(in) :: nscales               !< Number of LCT scales
real(kind_real),intent(in) :: D(4,nscales)  !< LCT components
real(kind_real),intent(in) :: coef(nscales) !< LCT coefficients
real(kind_real),intent(out) :: fit(nc,nl0)  !< Fit
 
! Local variables
integer :: jl0,jc3,iscales
real(kind_real) :: Dcoef(nscales),D11,D22,D33,D12,H11,H22,H33,H12,rsq
 
! Set name
 
 
! Probe in
 
 
! Initialization
fit = mpl%msv%valr
 
! Coefficients
dcoef = max(dmin,min(coef,one))
dcoef = dcoef/sum(dcoef)
 
do iscales=1,nscales
   ! Ensure positive-definiteness of D
   d11 = max(dmin,d(1,iscales))
   d22 = max(dmin,d(2,iscales))
   if (nl0>1) then
      d33 = max(dmin,d(3,iscales))
   else
      d33 = zero
   end if
   d12 = sqrt(d11*d22)*max(-one+dmin,min(d(4,iscales),one-dmin))
 
   ! Inverse D to get H
   call lct_d2h(mpl,d11,d22,d33,d12,h11,h22,h33,h12)
 
   ! Homogeneous anisotropic approximation
   do jl0=1,nl0
      do jc3=1,nc
         if (dmask(jc3,jl0)) then
            ! Initialization
            if (iscales==1) fit(jc3,jl0) = zero
 
            ! Squared distance
            rsq = h11*dxsq(jc3,jl0)+h22*dysq(jc3,jl0)+h33*dzsq(jc3,jl0)+two*h12*dxdy(jc3,jl0)
 
            if (m==0) then
               ! Gaussian function
               if (rsq<40.0_kind_real) fit(jc3,jl0) = fit(jc3,jl0)+dcoef(iscales)*exp(-half*rsq)
            else
               ! Matern function
               fit(jc3,jl0) = fit(jc3,jl0)+dcoef(iscales)*matern(mpl,m,sqrt(rsq))
            end if
         end if
      end do
   end do
end do
 
! Probe out
 
 
end subroutine func_fit_lct
 
!----------------------------------------------------------------------
! Subroutine: func_lct_d2h
!> From D (Daley tensor) to H (local correlation tensor)
!----------------------------------------------------------------------
subroutine func_lct_d2h(mpl,D11,D22,D33,D12,H11,H22,H33,H12)
 
implicit none
 
! Passed variables
type(mpl_type),intent(inout) :: mpl!< MPI data
real(kind_real),intent(in) :: D11  !< Daley tensor component 11
real(kind_real),intent(in) :: D22  !< Daley tensor component 22
real(kind_real),intent(in) :: D33  !< Daley tensor component 33
real(kind_real),intent(in) :: D12  !< Daley tensor component 12
real(kind_real),intent(out) :: H11 !< Local correlation tensor component 11
real(kind_real),intent(out) :: H22 !< Local correlation tensor component 22
real(kind_real),intent(out) :: H33 !< Local correlation tensor component 33
real(kind_real),intent(out) :: H12 !< Local correlation tensor component 12
 
! Local variables
real(kind_real) :: det
 
! Set name
 
 
! Probe in
 
 
! Compute horizontal determinant
det = d11*d22-d12**2
 
! Inverse D to get H
if (det>zero) then
   h11 = d22/det
   h22 = d11/det
   h12 = -d12/det
else
   call mpl%abort('func_lct_d2h','non-invertible tensor')
end if
if (d33>zero) then
   h33 = one/d33
else
   h33 = zero
end if
 
! Probe out
 
 
end subroutine func_lct_d2h
 
!----------------------------------------------------------------------
! Subroutine: func_lct_h2r
!> From H (local correlation tensor) to support radii
!----------------------------------------------------------------------
subroutine func_lct_h2r(mpl,H11,H22,H33,H12,rh,rv)
 
implicit none
 
! Passed variables
type(mpl_type),intent(inout) :: mpl !< MPI data
real(kind_real),intent(in) :: H11   !< Local correlation tensor component 11
real(kind_real),intent(in) :: H22   !< Local correlation tensor component 22
real(kind_real),intent(in) :: H33   !< Local correlation tensor component 33
real(kind_real),intent(in) :: H12   !< Local correlation tensor component 12
real(kind_real),intent(out) :: rh   !< Horizontal support radius
real(kind_real),intent(out) :: rv   !< Vertical support radius
 
! Local variables
real(kind_real) :: tr,det,diff
 
! Set name
 
 
! Probe in
 
 
! Check diagonal positivity
if ((h11<zero).or.(h22<zero)) call mpl%abort('func_lct_h2r','negative diagonal LCT coefficients')
 
! Compute horizontal trace
tr = h11+h22
 
! Compute horizontal determinant
det = h11*h22-h12**2
 
! Compute horizontal support radius
diff = quarter*(h11-h22)**2+h12**2
if ((det>zero).and..not.(diff<zero)) then
   if (half*tr>sqrt(diff)) then
      rh = gau2gc/sqrt(half*tr-sqrt(diff))
   else
      call mpl%abort('func_lct_h2r','non positive-definite LCT (eigenvalue)')
   end if
else
   call mpl%abort('func_lct_h2r','non positive-definite LCT (determinant)')
end if
 
! Compute vertical support radius
if (h33>zero) then
   rv = gau2gc/sqrt(h33)
else
   rv = zero
end if
 
! Probe out
 
 
end subroutine func_lct_h2r
 
!----------------------------------------------------------------------
! Subroutine: func_lct_r2d
!> From support radius to Daley tensor diagonal element
!----------------------------------------------------------------------
subroutine func_lct_r2d(r,D)
 
implicit none
 
! Passed variables
real(kind_real),intent(in) :: r  !< Support radius
real(kind_real),intent(out) :: D !< Daley tensor diagonal element
 
! Set name
 
 
! Probe in
 
 
! Convert from support radius to Daley length-scale and square
d = (gc2gau*r)**2
 
! Probe out
 
 
end subroutine func_lct_r2d
 
!----------------------------------------------------------------------
! Subroutine: func_check_cond
!> Check tensor conditioning
!----------------------------------------------------------------------
subroutine func_check_cond(d1,d2,nod,valid)
 
implicit none
 
! Passed variables
real(kind_real),intent(in) :: d1  !< First diagonal coefficient
real(kind_real),intent(in) :: d2  !< Second diagonal coefficient
real(kind_real),intent(in) :: nod !< Normalized off-diagonal coefficient
logical,intent(out) :: valid      !< Conditioning validity
 
! Local variables
real(kind_real) :: det,tr,diff,ev1,ev2
 
! Set name
 
 
! Probe in
 
 
! Compute trace and determinant
tr = d1+d2
det = d1*d2*(one-nod**2)
diff = quarter*(d1-d2)**2+d1*d2*nod**2
 
if ((det>zero).and..not.(diff<zero)) then
   ! Compute eigenvalues
   ev1 = half*tr+sqrt(diff)
   ev2 = half*tr-sqrt(diff)
 
   if (ev2>zero) then
      ! Check conditioning
      valid = inf(ev1,condmax*ev2)
   else
      ! Lowest negative eigenvalue is negative
      valid = .false.
   end if
else
   ! Non-positive definite tensor
   valid = .false.
end if
 
! Probe out
 
 
end subroutine func_check_cond
 
!----------------------------------------------------------------------
! Function: func_matern
!> Compute the normalized diffusion function from eq. (55) of Mirouze and Weaver (2013), for the 3d case (d = 3)
!----------------------------------------------------------------------
function func_matern(mpl,M,x) result(value)
 
implicit none
 
! Passed variables
type(mpl_type),intent(inout) :: mpl !< MPI data
integer,intent(in) :: m             !< Matern function order
real(kind_real),intent(in) :: x     !< Argument
 
! Returned variable
real(kind_real) :: value
 
! Local variables
integer :: j
real(kind_real) :: xtmp,beta
 
! Set name
 
 
! Probe in
 
 
! Check
if (m<2) call mpl%abort('func_matern','M should be larger than 2')
if (mod(m,2)>0) call mpl%abort('func_matern','M should be even')
 
! Initialization
value = zero
beta = one
xtmp = x*sqrt(real(2*m-5,kind_real))
 
do j=0,m-3
   ! Update sum
   value = value+beta*(xtmp)**(m-2-j)
 
   ! Update beta
   beta = beta*real((j+1+m-2)*(-j+m-2),kind_real)/real(2*(j+1),kind_real)
end do
 
! Last term and normalization
value = value/beta+one
 
! Exponential factor
value = value*exp(-xtmp)
 
! Probe out
 
 
end function func_matern
 
!----------------------------------------------------------------------
! Subroutine: func_cholesky
!> Compute cholesky decomposition
! Author: Original FORTRAN77 version by Michael Healy, modifications by AJ Miller, FORTRAN90 version by John Burkardt.
!----------------------------------------------------------------------
subroutine func_cholesky(mpl,n,a,u)
 
implicit none
 
! Passed variables
type(mpl_type),intent(inout) :: mpl   !< MPI data
integer,intent(in) :: n               !< Matrix rank
real(kind_real),intent(in) :: a(n,n)  !< Matrix
real(kind_real),intent(out) :: u(n,n) !< Matrix square-root
 
! Local variables
integer :: nn,i,j,ij
real(kind_real),allocatable :: apack(:),upack(:)
 
! Set name
 
 
! Probe in
 
 
! Allocation
nn = (n*(n+1))/2
allocate(apack(nn))
allocate(upack(nn))
 
! Pack matrix
ij = 0
do i=1,n
   do j=1,i
      ij = ij+1
      apack(ij) = a(i,j)
   end do
end do
 
! Cholesky decomposition
call cholesky(mpl,n,nn,apack,upack)
 
! Unpack matrix
ij = 0
u = zero
do i=1,n
   do j=1,i
      ij = ij+1
      u(i,j) = upack(ij)
   end do
end do
 
! Release memory
deallocate(apack)
deallocate(upack)
 
! Probe out
 
 
end subroutine func_cholesky
 
!----------------------------------------------------------------------
! Subroutine: func_syminv
!> Compute inverse of a symmetric matrix
! Author: Original FORTRAN77 version by Michael Healy, modifications by AJ Miller, FORTRAN90 version by John Burkardt.
!----------------------------------------------------------------------
subroutine func_syminv(mpl,n,a,c)
 
implicit none
 
! Passed variables
type(mpl_type),intent(inout) :: mpl   !< MPI data
integer,intent(in) :: n               !< Matrix rank
real(kind_real),intent(in) :: a(n,n)  !< Matrix
real(kind_real),intent(out) :: c(n,n) !< Matrix inverse
 
! Local variables
integer :: nn,i,j,ij
real(kind_real),allocatable :: apack(:),cpack(:)
 
! Set name
 
 
! Probe in
 
 
! Allocation
nn = (n*(n+1))/2
allocate(apack(nn))
allocate(cpack(nn))
 
! Pack matrix
ij = 0
do i=1,n
   do j=1,i
      ij = ij+1
      apack(ij) = a(i,j)
   end do
end do
 
! Matrix inversion
call syminv(mpl,n,nn,apack,cpack)
 
! Unpack matrix
ij = 0
do i=1,n
   do j=1,i
      ij = ij+1
      c(i,j) = cpack(ij)
      c(j,i) = c(i,j)
   end do
end do
 
! Release memory
deallocate(apack)
deallocate(cpack)
 
! Probe out
 
 
end subroutine func_syminv
 
!----------------------------------------------------------------------
! Subroutine: func_histogram
!> Compute bins and histogram from a list of values
!----------------------------------------------------------------------
subroutine func_histogram(mpl,nlist,list,nbins,histmin,histmax,bins,hist)
 
implicit none
 
! Passed variables
type(mpl_type),intent(inout) :: mpl          !< MPI data
integer,intent(in) :: nlist                  !< List size
real(kind_real),intent(in) :: list(nlist)    !< List
integer,intent(in) :: nbins                  !< Number of bins
real(kind_real),intent(in) :: histmin        !< Histogram minimum
real(kind_real),intent(in) :: histmax        !< Histogram maximum
real(kind_real),intent(out) :: bins(nbins+1) !< Bins
real(kind_real),intent(out) :: hist(nbins)   !< Histogram
 
! Local variables
integer :: ibins,ilist
real(kind_real) :: delta
logical :: found
 
! Set name
 
 
! Probe in
 
 
! Check data
if (nbins<=0) call mpl%abort('func_histogram','the number of bins should be positive')
if (histmax>histmin) then
   if (minval(list,mask=mpl%msv%isnot(list))<histmin) call mpl%abort('func_histogram','values below histogram minimum')
   if (maxval(list,mask=mpl%msv%isnot(list))>histmax) call mpl%abort('func_histogram','values over histogram maximum')
 
   ! Compute bins
   delta = (histmax-histmin)/real(nbins,kind_real)
   bins(1) = histmin
   do ibins=2,nbins
      bins(ibins) = histmin+real(ibins-1,kind_real)*delta
   end do
   bins(nbins+1) = histmax
 
   ! Extend first and last bins
   bins(1) = bins(1)-1.0e-6_kind_real*delta
   bins(nbins+1) = bins(nbins+1)+1.0e-6_kind_real*delta
 
   ! Compute histogram
   hist = zero
   do ilist=1,nlist
      if (mpl%msv%isnot(list(ilist))) then
         ibins = 0
         found = .false.
         do while (.not.found)
            ibins = ibins+1
            if (ibins>nbins) call mpl%abort('func_histogram','bin not found')
            if (infeq(bins(ibins),list(ilist)).and.inf(list(ilist),bins(ibins+1))) then
               hist(ibins) = hist(ibins)+one
               found = .true.
            end if
         end do
      end if
   end do
   if (abs(sum(hist)-real(count(mpl%msv%isnot(list)),kind_real))>half) &
 & call mpl%abort('func_histogram','histogram sum is not equal to the number of valid elements')
else
   bins = mpl%msv%valr
   hist = zero
end if
 
! Probe out
 
 
end subroutine func_histogram
 
!----------------------------------------------------------------------
! Function: func_cx_to_cxa
!> Conversion from global to halo A on subset Scx
!----------------------------------------------------------------------
function func_cx_to_cxa(nproc,proc_to_cx_offset,icx) result(icxa)
 
implicit none
 
! Passed variables
integer,intent(in) :: nproc                    !< Number of processors
integer,intent(in) :: proc_to_cx_offset(nproc) !< Processor to offset on subset Scx
integer,intent(in) :: icx                      !< Global index
 
! Returned variable
integer :: icxa
 
! Local variable
integer :: iproc
 
! Set name
 
 
! Probe in
 
 
! Find processor
iproc = cx_to_proc(nproc,proc_to_cx_offset,icx)
 
! Get halo A index
icxa = icx-proc_to_cx_offset(iproc)
 
! Probe out
 
 
end function func_cx_to_cxa
 
!----------------------------------------------------------------------
! Function: func_cx_to_proc
!> Conversion from global to processor on subset Scx
!----------------------------------------------------------------------
function func_cx_to_proc(nproc,proc_to_cx_offset,icx) result(iproc)
 
implicit none
 
! Passed variables
integer,intent(in) :: nproc                    !< Number of processors
integer,intent(in) :: proc_to_cx_offset(nproc) !< Processor to offset on subset Scx
integer,intent(in) :: icx                      !< Global index
 
! Returned variable
integer :: iproc
 
! Set name
 
 
! Probe in
 
 
! Find processor
do iproc=1,nproc-1
   if ((proc_to_cx_offset(iproc)<icx).and.(icx<=proc_to_cx_offset(iproc+1))) then
 
      return
   end if
end do
 
! Probe out
 
 
end function func_cx_to_proc
 
!----------------------------------------------------------------------
! Function: func_cx_to_cxu
!> Conversion from global to universe on subset Scx
!----------------------------------------------------------------------
function func_cx_to_cxu(nproc,proc_to_cx_offset,proc_to_ncxa,myuniverse,icx) result(icxu)
 
implicit none
 
! Passed variables
integer,intent(in) :: nproc                    !< Number of processors
integer,intent(in) :: proc_to_cx_offset(nproc) !< Processor to offset on subset Scx
integer,intent(in) :: proc_to_ncxa(nproc)      !< Processor to halo A size for subset Scx
logical,intent(in) :: myuniverse(nproc)        !< Task universe
integer,intent(in) :: icx                      !< Global index
 
! Returned variable
integer :: icxu
 
! Local variable
integer :: iproc,icxa,offset,jproc
 
! Set name
 
 
! Probe in
 
 
! Find processor
iproc = cx_to_proc(nproc,proc_to_cx_offset,icx)
 
if (myuniverse(iproc)) then
   ! Get halo A index
   icxa = icx-proc_to_cx_offset(iproc)
 
   ! Compute universe offset
   offset = 0
   do jproc=1,iproc-1
      if (myuniverse(jproc)) offset = offset+proc_to_ncxa(jproc)
   end do
 
   ! Get universe index
   icxu = offset+icxa
else
   ! Not in my universe
   icxu = 0
end if
 
! Probe out
 
 
end function func_cx_to_cxu
 
end module tools_func