tools_fit.F90

Go to the documentation of this file.

# 1 "/Users/miesch/JEDI/code/working_copy/internal/mpas-bundle/saber/src/saber/bump/tools_fit.fypp"
# 1 "/Users/miesch/JEDI/code/working_copy/internal/mpas-bundle/saber/src/saber/bump/../instrumentation.fypp" 1
# 1 "/Users/miesch/JEDI/code/working_copy/internal/mpas-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
!----------------------------------------------------------------------
 
# 926 "/Users/miesch/JEDI/code/working_copy/internal/mpas-bundle/saber/src/saber/bump/../subr_list.fypp"
# 2 "/Users/miesch/JEDI/code/working_copy/internal/mpas-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/internal/mpas-bundle/saber/src/saber/bump/../instrumentation.fypp"
# 2 "/Users/miesch/JEDI/code/working_copy/internal/mpas-bundle/saber/src/saber/bump/tools_fit.fypp" 2
!----------------------------------------------------------------------
! Module: tools_fit
!> Fit-related tools
! Author: Benjamin Menetrier
! Licensing: this code is distributed under the CeCILL-C license
! Copyright 2015-... UCAR, CERFACS, METEO-FRANCE and IRIT
!----------------------------------------------------------------------
module tools_fit
 
use tools_const, only: zero,quarter,half,one
use tools_func, only: fit_func
use tools_kinds, only: kind_real,huge_real
use tools_repro, only: inf,sup
use type_mpl, only: mpl_type
 
 
implicit none
 
integer,parameter :: itermax = 10 !< Maximum number of iteration for the threshold definition
 
interface fast_fit
   module procedure fit_fast_fit
end interface
interface ver_smooth
   module procedure fit_ver_smooth
end interface
interface ver_fill
   module procedure fit_ver_fill
end interface
 
private
public :: fast_fit,ver_smooth,ver_fill
 
contains
 
!----------------------------------------------------------------------
! Subroutine: fit_fast_fit
!> Fast fit length-scale estimation based on the value at mid-height
!----------------------------------------------------------------------
subroutine fit_fast_fit(mpl,n,iz,dist,raw,fit_r)
 
implicit none
 
! Passed variables
type(mpl_type),intent(inout) :: mpl   !< MPI data
integer,intent(in) :: n               !< Vector size
integer,intent(in) :: iz              !< Zero separation index
real(kind_real),intent(in) :: dist(n) !< Distance
real(kind_real),intent(in) :: raw(n)  !< Raw data
real(kind_real),intent(out) :: fit_r  !< Fast fit result
 
! Local variables
integer :: di,i,im,ip,iter
real(kind_real) :: th,thinv,dthinv,thtest
real(kind_real) :: fit_rm,fit_rp,distmin
real(kind_real) :: raw_tmp(n)
 
! Set name
 
 
! Probe in
 
 
if (any(dist<zero)) call mpl%abort('fit_fast_fit','negative distance')
 
if (raw(iz)>zero) then
   if (n>1) then
      ! Copy points that are lower than the zero-separation
      raw_tmp = mpl%msv%valr
      raw_tmp(iz) = one
      do i=1,n
        if (i/=iz) then
           if (inf(raw(i),raw(iz))) raw_tmp(i) = raw(i)/raw(iz)
        end if
      end do
 
      if (count(mpl%msv%isnot(raw_tmp))>1) then
         if (count(raw_tmp>zero)>1) then
            ! Curve-dependent threshold
            th = half*(one+minval(raw_tmp,mask=(raw_tmp>zero)))
 
            ! Find inverse threshold with a dichotomy
            thinv = half
            dthinv = quarter
            do iter=1,itermax
               thtest = fit_func(mpl,thinv)
               if (sup(th,thtest)) then
                  thinv = thinv-dthinv
               else
                  thinv = thinv+dthinv
               end if
               dthinv = half*dthinv
            end do
 
            ! Find support radius, lower value
            fit_rm = mpl%msv%valr
            ip = iz
            do di=1,n
               ! Check whether fit value has been found
               if (mpl%msv%is(fit_rm)) then
                  ! Index
                  im = iz-di
 
                  ! Check index validity
                  if (im>=1) then
                     ! Check raw value validity
                     if (raw_tmp(im)>zero) then
                        ! Check whether threshold has been crossed
                        if (inf(raw_tmp(im),th)) then
                           ! Set fit value
                           fit_rm = dist(im)+(dist(ip)-dist(im))*(th-raw_tmp(im))/(raw_tmp(ip)-raw_tmp(im))
                        else
                           ! Update index
                           ip = im
                        end if
                     end if
                  end if
               end if
            end do
 
            ! Find support radius, upper value
            fit_rp = mpl%msv%valr
            im = iz
            do di=1,n
               ! Check whether fit value has been found
               if (mpl%msv%is(fit_rp)) then
                  ! Index
                  ip = iz+di
 
                  ! Check index validity
                  if (ip<=n) then
                     ! Check raw value validity
                     if (raw_tmp(ip)>zero) then
                        ! Check whether threshold has been crossed
                        if (inf(raw_tmp(ip),th)) then
                           ! Set fit value
                           fit_rp = dist(im)+(dist(ip)-dist(im))*(th-raw_tmp(im))/(raw_tmp(ip)-raw_tmp(im))
                        else
                           ! Update index
                           im = ip
                        end if
                     end if
                  end if
               end if
            end do
 
            ! Gather values
            if (mpl%msv%isnot(fit_rm).and.mpl%msv%isnot(fit_rp)) then
               fit_r = half*(fit_rm+fit_rp)
            elseif (mpl%msv%isnot(fit_rm)) then
               fit_r = fit_rm
            elseif (mpl%msv%isnot(fit_rp)) then
               fit_r = fit_rp
            end if
 
            ! Normalize
            if (mpl%msv%isnot(fit_r)) fit_r = fit_r/thinv
 
            ! Check positivity
            if (inf(fit_r,zero)) fit_r = mpl%msv%valr
         else
            ! All positive-separation points are negative
            fit_r = zero
         end if
      else
         ! All positive-separation points are missing
         fit_r = mpl%msv%valr
      end if
 
      ! Set minimum distance
      if (mpl%msv%isnot(fit_r)) then
         distmin = huge_real
         if (iz>1) distmin = min(distmin,1.0e-6_kind_real*abs(dist(iz-1)-dist(iz)))
         if (iz<n) distmin = min(distmin,1.0e-6_kind_real*abs(dist(iz+1)-dist(iz)))
         fit_r = max(fit_r,distmin)
      end if
   else
      ! Only one point, zero radius
      fit_r = zero
   end if
else
   ! Zero-separation point is negative
   fit_r = mpl%msv%valr
end if
 
! Probe out
 
 
end subroutine fit_fast_fit
 
!----------------------------------------------------------------------
! Subroutine: fit_ver_smooth
!> Homogeneous smoothing of a vertical profile
!----------------------------------------------------------------------
subroutine fit_ver_smooth(mpl,n,x,rv,profile)
 
implicit none
 
! Passed variables
type(mpl_type),intent(inout) :: mpl         !< MPI data
integer,intent(in) :: n                     !< Vector size
real(kind_real),intent(in) :: x(n)          !< Coordinate
real(kind_real),intent(in) :: rv            !< Filtering support radius
real(kind_real),intent(inout) :: profile(n) !< Vertical profile
 
! Local variables
integer :: i,j
real(kind_real) :: kernel(n,n),distnorm,profile_init(n),norm
 
! Set name
 
 
! Probe in
 
 
if (n>1) then
   if (rv<zero) call mpl%abort('fit_ver_smooth','negative filtering support radius')
 
   if ((rv>zero).and.mpl%msv%isanynot(profile)) then
      ! Vertical smoothing kernel
      kernel = zero
      do i=1,n
         do j=1,n
            if (mpl%msv%isnot(profile(i)).and.mpl%msv%isnot(profile(j))) then
               ! Gaspari-Cohn (1999) function
               distnorm = abs(x(j)-x(i))/rv
               kernel(i,j) = fit_func(mpl,distnorm)
            end if
         end do
      end do
 
      ! Apply kernel
      profile_init = profile
      profile = zero
      do i=1,n
         norm = zero
         do j=1,n
            profile(i) = profile(i)+kernel(i,j)*profile_init(j)
            norm = norm+kernel(i,j)
         end do
         if (norm>zero) then
            profile(i) = profile(i)/norm
         else
            profile(i) = mpl%msv%valr
         end if
      end do
   end if
end if
 
! Probe out
 
 
end subroutine fit_ver_smooth
 
!----------------------------------------------------------------------
! Subroutine: fit_ver_fill
!> Missing values filling of a vertical profile
!----------------------------------------------------------------------
subroutine fit_ver_fill(mpl,n,x,profile)
 
implicit none
 
! Passed variables
type(mpl_type),intent(inout) :: mpl         !< MPI data
integer,intent(in) :: n                     !< Vector size
real(kind_real),intent(in) :: x(n)          !< Coordinate
real(kind_real),intent(inout) :: profile(n) !< Vertical profile
 
! Local variables
integer :: i,j,iinf,isup
real(kind_real) :: profile_init(n)
 
! Set name
 
 
! Probe in
 
 
if (mpl%msv%isanynot(profile)) then
   ! Initialization
   profile_init = profile
   iinf = mpl%msv%vali
 
   do i=1,n
      if (mpl%msv%isnot(profile_init(i))) then
         ! Valid inferior point
         iinf = i
      else
         ! Look for a superior point
         isup = mpl%msv%vali
         j = i+1
         do while ((j<=n).and.(mpl%msv%is(isup)))
            if (mpl%msv%isnot(profile_init(j))) isup = j
            j = j+1
         end do
 
         if (mpl%msv%isnot(iinf).and.mpl%msv%isnot(isup)) then
            ! Interpolation
            profile(i) = profile_init(iinf)+(x(i)-x(iinf))*(profile_init(isup)-profile_init(iinf))/(x(isup)-x(iinf))
         elseif (mpl%msv%isnot(isup)) then
            ! Extrapolation with nearest superior point
            profile(i) = profile(isup)
         elseif (mpl%msv%isnot(iinf)) then
            ! Extrapolation with nearest inferior point
            profile(i) = profile(iinf)
         else
            call mpl%abort('fit_ver_fill','too many missing values')
         end if
      end if
   end do
end if
 
! Probe out
 
 
end subroutine fit_ver_fill
 
end module tools_fit