34 integer,
intent(in) :: n
35 integer,
intent(in) :: iz
41 integer :: di,i,im,ip,iter
45 character(len=1024),
parameter :: subr =
'fast_fit'
47 if (any(dist<0.0))
call mpl%abort(subr,
'negative distance in fast_fit')
52 raw_tmp = mpl%msv%valr
56 if (
inf(raw(i),raw(iz))) raw_tmp(i) = raw(i)/raw(iz)
60 if (count(mpl%msv%isnot(raw_tmp))>1)
then
61 if (count(raw_tmp>0.0)>1)
then
63 th = 0.5*(1.0+minval(raw_tmp,mask=(raw_tmp>0.0)))
70 if (
sup(th,thtest))
then
83 if (mpl%msv%is(fit_rm))
then
90 if (raw_tmp(im)>0.0)
then
92 if (
inf(raw_tmp(im),th))
then
94 fit_rm = dist(im)+(dist(ip)-dist(im))*(th-raw_tmp(im))/(raw_tmp(ip)-raw_tmp(im))
105 fit_rp = mpl%msv%valr
109 if (mpl%msv%is(fit_rp))
then
116 if (raw_tmp(ip)>0.0)
then
118 if (
inf(raw_tmp(ip),th))
then
120 fit_rp = dist(im)+(dist(ip)-dist(im))*(th-raw_tmp(im))/(raw_tmp(ip)-raw_tmp(im))
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
135 elseif (mpl%msv%isnot(fit_rp))
then
140 if (mpl%msv%isnot(fit_r)) fit_r = fit_r/thinv
143 if (
inf(fit_r,0.0_kind_real)) fit_r = mpl%msv%valr
154 if (mpl%msv%isnot(fit_r))
then
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)
181 integer,
intent(in) :: n
184 real(
kind_real),
intent(inout) :: profile(n)
188 real(
kind_real) :: kernel(n,n),distnorm,profile_init(n),norm
189 character(len=1024),
parameter :: subr =
'ver_smooth'
191 if (rv<0.0)
call mpl%abort(subr,
'negative filtering support radius in ver_smooth')
193 if ((rv>0.0).and.mpl%msv%isanynot(profile))
then
198 if (mpl%msv%isnot(profile(j)))
then
200 distnorm = abs(x(j)-x(i))/rv
201 kernel(i,j) =
fit_func(mpl,distnorm)
207 profile_init = profile
212 profile(i) = profile(i)+kernel(i,j)*profile_init(j)
213 norm = norm+kernel(i,j)
216 profile(i) = profile(i)/norm
218 profile(i) = mpl%msv%valr
235 integer,
intent(in) :: n
237 real(
kind_real),
intent(inout) :: profile(n)
240 integer :: i,j,iinf,isup
242 character(len=1024),
parameter :: subr =
'ver_fill'
244 if (mpl%msv%isanynot(profile))
then
246 profile_init = profile
250 if (mpl%msv%isnot(profile_init(i)))
then
257 do while ((j<=n).and.(mpl%msv%is(isup)))
258 if (mpl%msv%isnot(profile_init(j))) isup = j
262 if (mpl%msv%isnot(iinf).and.mpl%msv%isnot(isup))
then
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
267 profile(i) = profile(isup)
268 elseif (mpl%msv%isnot(iinf))
then
270 profile(i) = profile(iinf)
272 call mpl%abort(subr,
'ver_fill failed')