SABER
tools_fit.F90
Go to the documentation of this file.
1 !----------------------------------------------------------------------
2 ! Module: tools_fit
3 !> Fit-related tools
4 ! Author: Benjamin Menetrier
5 ! Licensing: this code is distributed under the CeCILL-C license
6 ! Copyright © 2015-... UCAR, CERFACS, METEO-FRANCE and IRIT
7 !----------------------------------------------------------------------
8 module tools_fit
9 
10 use tools_func, only: fit_func
12 use tools_repro, only: inf,sup
13 use type_mpl, only: mpl_type
14 
15 implicit none
16 
17 integer,parameter :: itermax = 10 ! Maximum number of iteration for the threshold definition
18 
19 private
21 
22 contains
23 
24 !----------------------------------------------------------------------
25 ! Subroutine: fast_fit
26 !> Fast fit length-scale estimation based on the value at mid-height
27 !----------------------------------------------------------------------
28 subroutine fast_fit(mpl,n,iz,dist,raw,fit_r)
29 
30 implicit none
31 
32 ! Passed variables
33 type(mpl_type),intent(inout) :: mpl !< MPI data
34 integer,intent(in) :: n !< Vector size
35 integer,intent(in) :: iz !< Zero separation index
36 real(kind_real),intent(in) :: dist(n) !< Distance
37 real(kind_real),intent(in) :: raw(n) !< Raw data
38 real(kind_real),intent(out) :: fit_r !< Fast fit result
39 
40 ! Local variables
41 integer :: di,i,im,ip,iter
42 real(kind_real) :: th,thinv,dthinv,thtest
43 real(kind_real) :: fit_rm,fit_rp,distmin
44 real(kind_real) :: raw_tmp(n)
45 character(len=1024),parameter :: subr = 'fast_fit'
46 
47 if (any(dist<0.0)) call mpl%abort(subr,'negative distance in fast_fit')
48 
49 if (raw(iz)>0.0) then
50  if (n>1) then
51  ! Copy points that are lower than the zero-separation
52  raw_tmp = mpl%msv%valr
53  raw_tmp(iz) = 1.0
54  do i=1,n
55  if (i/=iz) then
56  if (inf(raw(i),raw(iz))) raw_tmp(i) = raw(i)/raw(iz)
57  end if
58  end do
59 
60  if (count(mpl%msv%isnot(raw_tmp))>1) then
61  if (count(raw_tmp>0.0)>1) then
62  ! Curve-dependent threshold
63  th = 0.5*(1.0+minval(raw_tmp,mask=(raw_tmp>0.0)))
64 
65  ! Find inverse threshold with a dichotomy
66  thinv = 0.5
67  dthinv = 0.25
68  do iter=1,itermax
69  thtest = fit_func(mpl,thinv)
70  if (sup(th,thtest)) then
71  thinv = thinv-dthinv
72  else
73  thinv = thinv+dthinv
74  end if
75  dthinv = 0.5*dthinv
76  end do
77 
78  ! Find support radius, lower value
79  fit_rm = mpl%msv%valr
80  ip = iz
81  do di=1,n
82  ! Check whether fit value has been found
83  if (mpl%msv%is(fit_rm)) then
84  ! Index
85  im = iz-di
86 
87  ! Check index validity
88  if (im>=1) then
89  ! Check raw value validity
90  if (raw_tmp(im)>0.0) then
91  ! Check whether threshold has been crossed
92  if (inf(raw_tmp(im),th)) then
93  ! Set fit value
94  fit_rm = dist(im)+(dist(ip)-dist(im))*(th-raw_tmp(im))/(raw_tmp(ip)-raw_tmp(im))
95  else
96  ! Update index
97  ip = im
98  end if
99  end if
100  end if
101  end if
102  end do
103 
104  ! Find support radius, upper value
105  fit_rp = mpl%msv%valr
106  im = iz
107  do di=1,n
108  ! Check whether fit value has been found
109  if (mpl%msv%is(fit_rp)) then
110  ! Index
111  ip = iz+di
112 
113  ! Check index validity
114  if (ip<=n) then
115  ! Check raw value validity
116  if (raw_tmp(ip)>0.0) then
117  ! Check whether threshold has been crossed
118  if (inf(raw_tmp(ip),th)) then
119  ! Set fit value
120  fit_rp = dist(im)+(dist(ip)-dist(im))*(th-raw_tmp(im))/(raw_tmp(ip)-raw_tmp(im))
121  else
122  ! Update index
123  im = ip
124  end if
125  end if
126  end if
127  end if
128  end do
129 
130  ! Gather values
131  if (mpl%msv%isnot(fit_rm).and.mpl%msv%isnot(fit_rp)) then
132  fit_r = 0.5*(fit_rm+fit_rp)
133  elseif (mpl%msv%isnot(fit_rm)) then
134  fit_r = fit_rm
135  elseif (mpl%msv%isnot(fit_rp)) then
136  fit_r = fit_rp
137  end if
138 
139  ! Normalize
140  if (mpl%msv%isnot(fit_r)) fit_r = fit_r/thinv
141 
142  ! Check positivity
143  if (inf(fit_r,0.0_kind_real)) fit_r = mpl%msv%valr
144  else
145  ! All positive-separation points are negative
146  fit_r = 0.0
147  end if
148  else
149  ! All positive-separation points are missing
150  fit_r = mpl%msv%valr
151  end if
152 
153  ! Set minimum distance
154  if (mpl%msv%isnot(fit_r)) then
155  distmin = huge_real
156  if (iz>1) distmin = min(distmin,1.0e-6*abs(dist(iz-1)-dist(iz)))
157  if (iz<n) distmin = min(distmin,1.0e-6*abs(dist(iz+1)-dist(iz)))
158  fit_r = max(fit_r,distmin)
159  end if
160  else
161  ! Only one point, zero radius
162  fit_r = 0.0
163  end if
164 else
165  ! Zero-separation point is negative
166  fit_r = mpl%msv%valr
167 end if
168 
169 end subroutine fast_fit
170 
171 !----------------------------------------------------------------------
172 ! Subroutine: ver_smooth
173 !> Homogeneous smoothing of a vertical profile
174 !----------------------------------------------------------------------
175 subroutine ver_smooth(mpl,n,x,rv,profile)
176 
177 implicit none
178 
179 ! Passed variables
180 type(mpl_type),intent(inout) :: mpl !< MPI data
181 integer,intent(in) :: n !< Vector size
182 real(kind_real),intent(in) :: x(n) !< Coordinate
183 real(kind_real),intent(in) :: rv !< Filtering support radius
184 real(kind_real),intent(inout) :: profile(n) !< Vertical profile
185 
186 ! Local variables
187 integer :: i,j
188 real(kind_real) :: kernel(n,n),distnorm,profile_init(n),norm
189 character(len=1024),parameter :: subr = 'ver_smooth'
190 
191 if (rv<0.0) call mpl%abort(subr,'negative filtering support radius in ver_smooth')
192 
193 if ((rv>0.0).and.mpl%msv%isanynot(profile)) then
194  ! Vertical smoothing kernel
195  kernel = 0.0
196  do i=1,n
197  do j=1,n
198  if (mpl%msv%isnot(profile(j))) then
199  ! Gaspari-Cohn (1999) function
200  distnorm = abs(x(j)-x(i))/rv
201  kernel(i,j) = fit_func(mpl,distnorm)
202  end if
203  end do
204  end do
205 
206  ! Apply kernel
207  profile_init = profile
208  profile = 0.0
209  do i=1,n
210  norm = 0.0
211  do j=1,n
212  profile(i) = profile(i)+kernel(i,j)*profile_init(j)
213  norm = norm+kernel(i,j)
214  end do
215  if (norm>0.0) then
216  profile(i) = profile(i)/norm
217  else
218  profile(i) = mpl%msv%valr
219  end if
220  end do
221 end if
222 
223 end subroutine ver_smooth
224 
225 !----------------------------------------------------------------------
226 ! Subroutine: ver_fill
227 !> Missing values filling of a vertical profile
228 !----------------------------------------------------------------------
229 subroutine ver_fill(mpl,n,x,profile)
230 
231 implicit none
232 
233 ! Passed variables
234 type(mpl_type),intent(inout) :: mpl !< MPI data
235 integer,intent(in) :: n !< Vector size
236 real(kind_real),intent(in) :: x(n) !< Coordinate
237 real(kind_real),intent(inout) :: profile(n) !< Vertical profile
238 
239 ! Local variables
240 integer :: i,j,iinf,isup
241 real(kind_real) :: profile_init(n)
242 character(len=1024),parameter :: subr = 'ver_fill'
243 
244 if (mpl%msv%isanynot(profile)) then
245  ! Initialization
246  profile_init = profile
247  iinf = mpl%msv%vali
248 
249  do i=1,n
250  if (mpl%msv%isnot(profile_init(i))) then
251  ! Valid inferior point
252  iinf = i
253  else
254  ! Look for a superior point
255  isup = mpl%msv%vali
256  j = i+1
257  do while ((j<=n).and.(mpl%msv%is(isup)))
258  if (mpl%msv%isnot(profile_init(j))) isup = j
259  j = j+1
260  end do
261 
262  if (mpl%msv%isnot(iinf).and.mpl%msv%isnot(isup)) then
263  ! Interpolation
264  profile(i) = profile_init(iinf)+(x(i)-x(iinf))*(profile_init(isup)-profile_init(iinf))/(x(isup)-x(iinf))
265  elseif (mpl%msv%isnot(isup)) then
266  ! Extrapolation with nearest superior point
267  profile(i) = profile(isup)
268  elseif (mpl%msv%isnot(iinf)) then
269  ! Extrapolation with nearest inferior point
270  profile(i) = profile(iinf)
271  else
272  call mpl%abort(subr,'ver_fill failed')
273  end if
274  end if
275  end do
276 end if
277 
278 end subroutine ver_fill
279 
280 end module tools_fit
tools_repro::sup
logical function, public sup(x, y)
Superior test for reals.
Definition: tools_repro.F90:92
tools_repro::inf
logical function, public inf(x, y)
Inferior test for reals.
Definition: tools_repro.F90:53
tools_func::fit_func
real(kind_real) function, public fit_func(mpl, distnorm)
Fit_function.
Definition: tools_func.F90:584
tools_fit
Fit-related tools.
Definition: tools_fit.F90:8
tools_fit::itermax
integer, parameter itermax
Definition: tools_fit.F90:17
tools_func
Usual functions.
Definition: tools_func.F90:8
tools_fit::ver_smooth
subroutine, public ver_smooth(mpl, n, x, rv, profile)
Homogeneous smoothing of a vertical profile.
Definition: tools_fit.F90:176
tools_fit::fast_fit
subroutine, public fast_fit(mpl, n, iz, dist, raw, fit_r)
Fast fit length-scale estimation based on the value at mid-height.
Definition: tools_fit.F90:29
tools_repro
Reproducibility functions.
Definition: tools_repro.F90:8
tools_kinds::huge_real
real(kind_real), parameter, public huge_real
Definition: tools_kinds.F90:25
tools_kinds
Kinds definition.
Definition: tools_kinds.F90:8
type_mpl
MPI parameters derived type.
Definition: type_mpl.F90:8
tools_fit::ver_fill
subroutine, public ver_fill(mpl, n, x, profile)
Missing values filling of a vertical profile.
Definition: tools_fit.F90:230
type_mpl::mpl_type
Definition: type_mpl.F90:24
tools_kinds::kind_real
integer, parameter, public kind_real
Definition: tools_kinds.F90:18