SCALE-RM
Data Types | Functions/Subroutines
scale_interp Module Reference

module INTERPOLATION More...

Functions/Subroutines

subroutine, public interp_setup (weight_order, search_limit)
 Setup. More...
 
subroutine, public interp_domain_compatibility (lon_org, lat_org, topc_org, lon_loc, lat_loc, topc_loc, topf_loc, skip_x, skip_y, skip_z)
 
subroutine, public interp_factor1d (KA_ref, KS_ref, KE_ref, KA, KS, KE, hgt_ref, hgt, idx_k, vfact, flag_extrap)
 
subroutine, public interp_factor2d_linear_latlon (IA_ref, JA_ref, IA, JA, lon_ref, lat_ref, lon, lat, idx_i, idx_j, hfact)
 
subroutine, public interp_factor2d_linear_xy (IA_ref, JA_ref, IA, JA, x_ref, y_ref, x, y, idx_i, idx_j, hfact, zonal, pole)
 
subroutine, public interp_factor2d_weight (npoints, IA_ref, JA_ref, IA, JA, lon_ref, lat_ref, lon, lat, idx_i, idx_j, hfact, search_limit, latlon_structure, lon_1d, lat_1d, weight_order)
 
subroutine, public interp_factor3d_linear_latlon (KA_ref, KS_ref, KE_ref, IA_ref, JA_ref, KA, KS, KE, IA, JA, lon_ref, lat_ref, hgt_ref, lon, lat, hgt, idx_i, idx_j, hfact, idx_k, vfact, flag_extrap)
 
subroutine, public interp_factor3d_linear_xy (KA_ref, KS_ref, KE_ref, IA_ref, JA_ref, KA, KS, KE, IA, JA, x_ref, y_ref, hgt_ref, x, y, hgt, idx_i, idx_j, hfact, idx_k, vfact, flag_extrap, zonal, pole)
 
subroutine, public interp_factor3d_weight (npoints, KA_ref, KS_ref, KE_ref, IA_ref, JA_ref, KA, KS, KE, IA, JA, lon_ref, lat_ref, hgt_ref, lon, lat, hgt, idx_i, idx_j, hfact, idx_k, vfact, flag_extrap)
 
subroutine, public interp_interp1d (KA_ref, KS_ref, KE_ref, KA, KS, KE, idx_k, vfact, hgt_ref, hgt, val_ref, val, spline, logwgt)
 
subroutine, public interp_interp2d (npoints, IA_ref, JA_ref, IA, JA, idx_i, idx_j, hfact, val_ref, val, threshold_undef, wsum, val2)
 
subroutine, public interp_interp3d (npoints, KA_ref, KS_ref, KE_ref, IA_ref, JA_ref, KA, KS, KE, IA, JA, idx_i, idx_j, hfact, idx_k, vfact, hgt_ref, hgt, val_ref, val, spline, logwgt, threshold_undef, wsum, val2)
 
subroutine interp_search_horiz_struct (npoints, psizex, psizey, IA_ref, JA_ref, lon_ref, lat_ref, lon_min, lat_min, dlon, dlat, i0, i1, j0, j1, lon, lat, idx_i, idx_j, hfact, search_limit, weight_order)
 
subroutine interp_insert_2d (npoints, lon, lat, lon_ref, lat_ref, i, j, drad, idx_i, idx_j)
 
subroutine interp_div_block (nsize, psize, nidx_max, lon_ref, lat_ref, idx, nidx, lon_min, lon_max, lat_min, lat_max, dlon, dlat)
 
subroutine interp_bilinear_inv (x_ref0, x_ref1, x_ref2, x_ref3, y_ref0, y_ref1, y_ref2, y_ref3, x, y, u, v, error)
 
real(rp) function cross (x0, y0, x1, y1)
 cross product More...
 
subroutine spline_exec (KA_ref, kmax, KA, KS, KE, idx_k, vfact, hgt_ref, hgt, val_ref, idx, idx_r, U, FDZ, val)
 

Detailed Description

module INTERPOLATION

Description
INTERPOLATION module for nesting system
Author
Team SCALE
NAMELIST
  • PARAM_INTERP
    nametypedefault valuecomment
    INTERP_BUFFER_SIZE_FACT real(RP) 2.0_RP
    INTERP_USE_SPLINE_VERT logical .true.
    INTERP_THRESHOLD_UNDEF real(RP) 1.0_RP [0-1]

History Output
No history output

Function/Subroutine Documentation

◆ interp_setup()

subroutine, public scale_interp::interp_setup ( integer, intent(in)  weight_order,
real(rp), intent(in), optional  search_limit 
)

Setup.

Definition at line 88 of file scale_interp.F90.

88  use scale_prc, only: &
89  prc_abort
90  implicit none
91 
92  integer, intent(in) :: weight_order
93  real(RP), intent(in), optional :: search_limit
94 
95  namelist /param_interp/ &
96  interp_buffer_size_fact, &
97  interp_use_spline_vert, &
98  interp_threshold_undef
99 
100  integer :: ierr
101  !---------------------------------------------------------------------------
102 
103  log_newline
104  log_info("INTERP_setup",*) 'Setup'
105 
106  interp_weight_order = weight_order
107 
108  interp_search_limit = large_number
109  if ( present(search_limit) ) then
110  interp_search_limit = search_limit
111  log_info("INTERP_setup",*) 'search limit [m] : ', interp_search_limit
112  endif
113 
114  !--- read namelist
115  rewind(io_fid_conf)
116  read(io_fid_conf,nml=param_interp,iostat=ierr)
117  if( ierr < 0 ) then !--- missing
118  log_info("INTERP_setup",*) 'Not found namelist. Default used.'
119  elseif( ierr > 0 ) then !--- fatal error
120  log_error("INTERP_setup",*) 'Not appropriate names in namelist PARAM_INTERP. Check!'
121  call prc_abort
122  endif
123  log_nml(param_interp)
124 
125  if ( rp == 8 ) then
126  eps_bilinear = 1e-6_rp
127  else
128  eps_bilinear = 1e-2_rp
129  end if
130 
131  return

References scale_io::io_fid_conf, and scale_prc::prc_abort().

Referenced by scale_comm_cartesc_nest::comm_cartesc_nest_setup().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ interp_domain_compatibility()

subroutine, public scale_interp::interp_domain_compatibility ( real(rp), dimension (:,:), intent(in)  lon_org,
real(rp), dimension (:,:), intent(in)  lat_org,
real(rp), dimension(:,:), intent(in)  topc_org,
real(rp), dimension (:,:), intent(in)  lon_loc,
real(rp), dimension (:,:), intent(in)  lat_loc,
real(rp), dimension(:,:), intent(in)  topc_loc,
real(rp), dimension(:,:), intent(in)  topf_loc,
logical, intent(in), optional  skip_x,
logical, intent(in), optional  skip_y,
logical, intent(in), optional  skip_z 
)

Definition at line 146 of file scale_interp.F90.

146  use scale_const, only: &
147  d2r => const_d2r
148  use scale_prc, only: &
149  prc_abort
150  implicit none
151 
152  real(RP), intent(in) :: lon_org (:,:)
153  real(RP), intent(in) :: lat_org (:,:)
154  real(RP), intent(in) :: topc_org(:,:) ! full level
155  real(RP), intent(in) :: lon_loc (:,:)
156  real(RP), intent(in) :: lat_loc (:,:)
157  real(RP), intent(in) :: topc_loc(:,:) ! full level
158  real(RP), intent(in) :: topf_loc(:,:) ! half level (ceil)
159 
160  logical, intent(in), optional :: skip_z
161  logical, intent(in), optional :: skip_x
162  logical, intent(in), optional :: skip_y
163 
164  real(RP) :: max_ref, min_ref
165  real(RP) :: max_loc, min_loc
166 
167  logical :: skip_z_
168  logical :: skip_x_
169  logical :: skip_y_
170 
171  intrinsic size
172  !---------------------------------------------------------------------------
173 
174  skip_z_ = .false.
175  if ( present(skip_z) ) then
176  skip_z_ = skip_z
177  endif
178 
179  skip_x_ = .false.
180  if ( present(skip_x) ) then
181  skip_x_ = skip_x
182  endif
183 
184  skip_y_ = .false.
185  if ( present(skip_y) ) then
186  skip_y_ = skip_y
187  endif
188 
189  if ( .NOT. skip_z_ ) then
190  max_ref = maxval( topc_org(:,:) ) + minval( topf_loc(:,:) - topc_loc(:,:) ) ! not real top boundary (only for check)
191  max_loc = maxval( topc_loc(:,:) )
192 
193  if ( max_ref < max_loc ) then
194  log_error("INTERP_domain_compatibility",*) 'REQUESTED DOMAIN IS TOO MUCH BROAD'
195  log_error_cont(*) '-- VERTICAL direction over the limit'
196  log_error_cont(*) '-- reference max: ', max_ref
197  log_error_cont(*) '-- local max: ', max_loc
198  call prc_abort
199  endif
200  endif
201 
202  if ( .NOT. skip_x_ ) then
203  max_ref = maxval( lon_org(:,:) / d2r )
204  min_ref = minval( lon_org(:,:) / d2r )
205  max_loc = maxval( lon_loc(:,:) / d2r )
206  min_loc = minval( lon_loc(:,:) / d2r )
207 
208  if ( (min_ref+360.0_rp-max_ref) < 360.0_rp / size(lon_org,1) * 2.0_rp ) then
209  ! cyclic OK
210  elseif( max_ref < max_loc .OR. min_ref > min_loc ) then
211  log_error("INTERP_domain_compatibility",*) 'REQUESTED DOMAIN IS TOO MUCH BROAD'
212  log_error_cont(*) '-- LONGITUDINAL direction over the limit'
213  log_error_cont(*) '-- reference max: ', max_ref
214  log_error_cont(*) '-- reference min: ', min_ref
215  log_error_cont(*) '-- local max: ', max_loc
216  log_error_cont(*) '-- local min: ', min_loc
217  call prc_abort
218  endif
219  endif
220 
221  if ( .NOT. skip_y_ ) then
222  max_ref = maxval( lat_org(:,:) / d2r )
223  min_ref = minval( lat_org(:,:) / d2r )
224  max_loc = maxval( lat_loc(:,:) / d2r )
225  min_loc = minval( lat_loc(:,:) / d2r )
226 
227  if ( max_ref < max_loc .OR. min_ref > min_loc ) then
228  log_error("INTERP_domain_compatibility",*) 'REQUESTED DOMAIN IS TOO MUCH BROAD'
229  log_error_cont(*) '-- LATITUDINAL direction over the limit'
230  log_error_cont(*) '-- reference max: ', max_ref
231  log_error_cont(*) '-- reference min: ', min_ref
232  log_error_cont(*) '-- local max: ', max_loc
233  log_error_cont(*) '-- local min: ', min_loc
234  call prc_abort
235  endif
236  endif
237 
238  return

References scale_const::const_d2r, and scale_prc::prc_abort().

Referenced by mod_copytopo::copytopo(), and mod_realinput::realinput_surface().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ interp_factor1d()

subroutine, public scale_interp::interp_factor1d ( integer, intent(in)  KA_ref,
integer, intent(in)  KS_ref,
integer, intent(in)  KE_ref,
integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
real(rp), dimension(ka_ref), intent(in)  hgt_ref,
real(rp), dimension (ka), intent(in)  hgt,
integer, dimension (ka,2), intent(out)  idx_k,
real(rp), dimension (ka), intent(out)  vfact,
logical, intent(in), optional  flag_extrap 
)

Definition at line 252 of file scale_interp.F90.

252  use scale_prc, only: &
253  prc_abort
254  use scale_const, only: &
255  undef => const_undef, &
256  eps => const_eps
257  implicit none
258  integer, intent(in) :: KA_ref, KS_ref, KE_ref ! number of z-direction (reference)
259  integer, intent(in) :: KA, KS, KE ! number of z-direction (target)
260 
261  real(RP), intent(in) :: hgt_ref(KA_ref) ! height [m] (reference)
262  real(RP), intent(in) :: hgt (KA) ! height [m] (target)
263 
264  integer, intent(out) :: idx_k (KA,2) ! k-index in reference (target)
265  real(RP), intent(out) :: vfact (KA) ! horizontal interp factor (target)
266 
267  logical, intent(in), optional :: flag_extrap ! when true, extrapolation will be executed (just copy)
268 
269  logical :: flag_extrap_
270 
271  integer :: k, kk
272  !---------------------------------------------------------------------------
273 
274  if ( present(flag_extrap) ) then
275  flag_extrap_ = flag_extrap
276  else
277  flag_extrap_ = .true.
278  end if
279 
280  do k = ks, ke
281  idx_k(k,1) = -1
282  idx_k(k,2) = -1
283 
284  if ( hgt(k) < hgt_ref(ks_ref) - eps ) then
285  if ( flag_extrap_ ) then
286  idx_k(k,1) = ks_ref
287  idx_k(k,2) = -1
288  vfact(k) = 1.0_rp
289  else
290  idx_k(k,1) = -1
291  idx_k(k,2) = -1
292  vfact(k) = undef
293  end if
294  elseif( hgt(k) < hgt_ref(ks_ref) ) then
295  idx_k(k,1) = ks_ref
296  idx_k(k,2) = -1
297  vfact(k) = 1.0_rp
298  elseif( hgt(k) > hgt_ref(ke_ref) + eps ) then
299  if ( flag_extrap_ ) then
300  idx_k(k,1) = ke_ref
301  idx_k(k,2) = -1
302  vfact(k) = 1.0_rp
303  else
304  idx_k(k,1) = -1
305  idx_k(k,2) = -1
306  vfact(k) = undef
307  end if
308  elseif( hgt(k) >= hgt_ref(ke_ref) ) then
309  idx_k(k,1) = ke_ref
310  idx_k(k,2) = -1
311  vfact(k) = 1.0_rp
312  else
313  do kk = ks_ref, ke_ref-1
314  if ( hgt(k) >= hgt_ref(kk ) &
315  .AND. hgt(k) < hgt_ref(kk+1) ) then
316  idx_k(k,1) = kk
317  idx_k(k,2) = kk + 1
318  vfact(k) = ( hgt_ref(kk+1) - hgt(k) ) &
319  / ( hgt_ref(kk+1) - hgt_ref(kk) )
320 
321  exit
322  endif
323  enddo
324  endif
325 
326  enddo ! k-loop
327 
328  return

References scale_const::const_eps, scale_const::const_undef, and scale_prc::prc_abort().

Referenced by interp_factor3d_linear_latlon(), interp_factor3d_linear_xy(), interp_factor3d_weight(), scale_interp_vert::interp_vert_setcoef(), and scale_interp_vert::interp_vert_setcoef_pres().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ interp_factor2d_linear_latlon()

subroutine, public scale_interp::interp_factor2d_linear_latlon ( integer, intent(in)  IA_ref,
integer, intent(in)  JA_ref,
integer, intent(in)  IA,
integer, intent(in)  JA,
real(rp), dimension(ia_ref), intent(in)  lon_ref,
real(rp), dimension(ja_ref), intent(in)  lat_ref,
real(rp), dimension(ia,ja), intent(in)  lon,
real(rp), dimension(ia,ja), intent(in)  lat,
integer, dimension(ia,ja,4), intent(out)  idx_i,
integer, dimension(ia,ja,4), intent(out)  idx_j,
real(rp), dimension(ia,ja,4), intent(out)  hfact 
)

Definition at line 341 of file scale_interp.F90.

341  use scale_sort, only: &
342  sort_exec
343  implicit none
344  integer, intent(in) :: IA_ref, JA_ref
345  integer, intent(in) :: IA, JA
346 
347  real(RP), intent(in) :: lon_ref(IA_ref)
348  real(RP), intent(in) :: lat_ref(JA_ref)
349  real(RP), intent(in) :: lon(IA,JA)
350  real(RP), intent(in) :: lat(IA,JA)
351  integer, intent(out) :: idx_i(IA,JA,4)
352  integer, intent(out) :: idx_j(IA,JA,4)
353  real(RP), intent(out) :: hfact(IA,JA,4)
354 
355  real(RP) :: f1, f2
356  integer :: i, j, ii, jj
357 
358  call prof_rapstart('INTERP_fact',3)
359 
360  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
361  !$omp private(f1,f2)
362  do j = 1, ja
363  do i = 1, ia
364 
365  ! longitude
366  if ( lon(i,j) < lon_ref(1) ) then
367  idx_i(i,j,:) = 1
368  hfact(i,j,:) = 0.0_rp
369  else if ( lon(i,j) > lon_ref(ia_ref) ) then
370  idx_i(i,j,:) = ia_ref
371  hfact(i,j,:) = 0.0_rp
372  else if ( lon(i,j) == lon_ref(ia_ref) ) then
373  idx_i(i,j,:) = ia_ref
374  hfact(i,j,1) = 0.0_rp
375  hfact(i,j,2) = 1.0_rp
376  hfact(i,j,3) = 0.0_rp
377  hfact(i,j,4) = 1.0_rp
378  else
379  do ii = 1, ia_ref
380  if ( lon(i,j) < lon_ref(ii) ) then
381  f1 = ( lon_ref(ii) - lon(i,j) ) / ( lon_ref(ii) - lon_ref(ii-1) )
382  f2 = ( lon(i,j) - lon_ref(ii-1) ) / ( lon_ref(ii) - lon_ref(ii-1) )
383  hfact(i,j,1) = f1
384  hfact(i,j,2) = f2
385  hfact(i,j,3) = f1
386  hfact(i,j,4) = f2
387  idx_i(i,j,1) = ii-1
388  idx_i(i,j,2) = ii
389  idx_i(i,j,3) = ii-1
390  idx_i(i,j,4) = ii
391  exit
392  end if
393  end do
394  end if
395 
396  ! latitude
397  if ( lat(i,j) < lat_ref(1) ) then
398  idx_j(i,j,:) = 1
399  hfact(i,j,:) = 0.0_rp
400  else if ( lat(i,j) > lat_ref(ja_ref) ) then
401  idx_j(i,j,:) = ja_ref
402  hfact(i,j,:) = 0.0_rp
403  else if ( lat(i,j) == lat_ref(ja_ref) ) then
404  idx_j(i,j,:) = ja_ref
405  hfact(i,j,1) = 0.0_rp
406  hfact(i,j,2) = 0.0_rp
407  !hfact(i,j,3) = hfact(i,j,3)
408  !hfact(i,j,4) = hfact(i,j,4)
409  else
410  do jj = 1, ja_ref
411  if ( lat(i,j) < lat_ref(jj) ) then
412  f1 = ( lat_ref(jj) - lat(i,j) ) / ( lat_ref(jj) - lat_ref(jj-1) )
413  f2 = ( lat(i,j) - lat_ref(jj-1) ) / ( lat_ref(jj) - lat_ref(jj-1) )
414  hfact(i,j,1) = hfact(i,j,1) * f1
415  hfact(i,j,2) = hfact(i,j,2) * f1
416  hfact(i,j,3) = hfact(i,j,3) * f2
417  hfact(i,j,4) = hfact(i,j,4) * f2
418  idx_j(i,j,1) = jj-1
419  idx_j(i,j,2) = jj-1
420  idx_j(i,j,3) = jj
421  idx_j(i,j,4) = jj
422  exit
423  end if
424  end do
425  end if
426 
427  call sort_exec( 4, & ! [IN]
428  hfact(i,j,:), idx_i(i,j,:), idx_j(i,j,:), & ! [INOUT]
429  reverse = .true. ) ! [IN]
430 
431  end do
432  end do
433 
434  call prof_rapend ('INTERP_fact',3)
435 
436  return

References scale_prof::prof_rapend(), and scale_prof::prof_rapstart().

Referenced by mod_cnv2d::cnv2d_exec(), mod_cnvuser::cnvuser_setup(), and interp_factor3d_linear_latlon().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ interp_factor2d_linear_xy()

subroutine, public scale_interp::interp_factor2d_linear_xy ( integer, intent(in)  IA_ref,
integer, intent(in)  JA_ref,
integer, intent(in)  IA,
integer, intent(in)  JA,
real(rp), dimension(ia_ref,ja_ref), intent(in)  x_ref,
real(rp), dimension(ia_ref,ja_ref), intent(in)  y_ref,
real(rp), dimension(ia), intent(in)  x,
real(rp), dimension(ja), intent(in)  y,
integer, dimension(ia,ja,4), intent(out)  idx_i,
integer, dimension(ia,ja,4), intent(out)  idx_j,
real(rp), dimension(ia,ja,4), intent(out)  hfact,
logical, intent(in), optional  zonal,
logical, intent(in), optional  pole 
)

Definition at line 450 of file scale_interp.F90.

450  use scale_prc, only: &
451  prc_abort
452  use scale_const, only: &
453  undef => const_undef, &
454  eps => const_eps
455  use scale_sort, only: &
456  sort_exec
457  implicit none
458  integer, intent(in) :: IA_ref, JA_ref
459  integer, intent(in) :: IA, JA
460 
461  real(RP), intent(in) :: x_ref(IA_ref,JA_ref)
462  real(RP), intent(in) :: y_ref(IA_ref,JA_ref)
463  real(RP), intent(in) :: x(IA)
464  real(RP), intent(in) :: y(JA)
465 
466  integer, intent(out) :: idx_i(IA,JA,4)
467  integer, intent(out) :: idx_j(IA,JA,4)
468  real(RP), intent(out) :: hfact(IA,JA,4)
469 
470  logical, intent(in), optional :: zonal
471  logical, intent(in), optional :: pole
472 
473  real(RP) :: u, v
474  integer :: inc_i, inc_j
475  logical :: error, err
476  logical :: zonal_, pole_
477 
478  integer :: ii0, jj0
479  integer :: ii, jj
480  integer :: i1, i2, i3, i4
481  integer :: j1, j2, j3, j4
482  integer :: i, j
483  integer :: ite, ite_max
484 
485  call prof_rapstart('INTERP_fact',3)
486 
487  if ( present(zonal) ) then
488  zonal_ = zonal
489  else
490  zonal_ = .false.
491  end if
492  if ( present(pole) ) then
493  pole_ = pole
494  else
495  pole_ = .false.
496  end if
497 
498  ii0 = ia_ref * 0.5_rp
499  jj0 = ja_ref * 0.5_rp
500  error = .false.
501 
502  ite_max = ia_ref * ja_ref / 4
503 
504  !$omp parallel do &
505  !$omp private(inc_i,inc_j,ii,jj,i1,i2,i3,i4,j1,j2,j3,j4,u,v,err) &
506  !$omp firstprivate(ii0,jj0)
507  do j = 1, ja
508  do i = 1, ia
509 
510  if ( i==1 ) then
511  ii = ii0
512  jj = jj0
513  end if
514  do ite = 1, ite_max
515  i1 = ii
516  i2 = ii + 1
517  i3 = ii + 1
518  i4 = ii
519  j1 = jj
520  j2 = jj
521  j3 = jj + 1
522  j4 = jj + 1
523  if ( ii == ia_ref ) then ! zonal_ must be .true.
524  i2 = 1
525  i3 = 1
526  end if
527  if ( jj == 0 ) then ! pole_ must be .true.
528  j1 = 1
529  j2 = 1
530  i2 = ia_ref / 2 + ii
531  if ( i2 > ia_ref ) i2 = i2 - ia_ref
532  i1 = i2 + 1
533  if ( i1 > ia_ref ) i1 = i1 - ia_ref
534  end if
535  if ( jj == ja_ref ) then ! pole_ must be .true.
536  j3 = ja_ref
537  j4 = ja_ref
538  i3 = ia_ref / 2 + ii
539  if ( i3 > ia_ref ) i3 = i3 - ia_ref
540  i4 = i3 + 1
541  if ( i4 > ia_ref ) i4 = i4 - ia_ref
542  end if
543  if ( abs( x_ref(i1,j1) - undef ) < eps .or. &
544  abs( x_ref(i2,j2) - undef ) < eps .or. &
545  abs( x_ref(i3,j3) - undef ) < eps .or. &
546  abs( x_ref(i4,j4) - undef ) < eps ) then
547  if ( ii == ia_ref-1 ) then
548  ii = 1
549  if ( jj == ja_ref-1 ) then
550  jj = 1
551  else
552  jj = jj + 2
553  end if
554  else
555  ii = ii + 2
556  end if
557  if ( ii == ia_ref ) then
558  ii = 1
559  jj = jj + 2
560  end if
561  if ( jj == ja_ref ) then
562  jj = 1
563  ii = ii + 2
564  if ( ii == ia_ref ) ii = 1
565  end if
566  else
567  call interp_check_inside( &
568  x_ref(i1,j1), x_ref(i2,j2), x_ref(i3,j3), x_ref(i4,j4), & ! (in)
569  y_ref(i1,j1), y_ref(i2,j2), y_ref(i3,j3), y_ref(i4,j4), & ! (in)
570  x(i), y(j), & ! (in)
571  inc_i, inc_j, & ! (out)
572  err ) ! (out)
573  if ( err ) error = .true.
574  if ( error ) exit
575  if ( inc_i == 0 .and. inc_j == 0 ) then ! inside the quadrilateral
576  call interp_bilinear_inv( &
577  x_ref(i1,j1), x_ref(i2,j2), x_ref(i3,j3), x_ref(i4,j4), & ! (in)
578  y_ref(i1,j1), y_ref(i2,j2), y_ref(i3,j3), y_ref(i4,j4), & ! (in)
579  x(i), y(j), & ! (in)
580  u, v, & ! (out)
581  err ) ! (out)
582  if ( err ) error = .true.
583  if ( error ) exit
584  idx_i(i,j,1) = i1
585  idx_i(i,j,2) = i2
586  idx_i(i,j,3) = i3
587  idx_i(i,j,4) = i4
588  idx_j(i,j,1) = j1
589  idx_j(i,j,2) = j2
590  idx_j(i,j,3) = j3
591  idx_j(i,j,4) = j4
592  hfact(i,j,1) = ( 1.0_rp - u ) * ( 1.0_rp - v )
593  hfact(i,j,2) = ( u ) * ( 1.0_rp - v )
594  hfact(i,j,3) = ( u ) * ( v )
595  hfact(i,j,4) = ( 1.0_rp - u ) * ( v )
596  exit
597  end if
598  ii = ii + inc_i
599  jj = jj + inc_j
600  if ( zonal_ .or. pole_ ) then
601  if ( ii == 0 ) ii = ia_ref
602  if ( ii == ia_ref+1 ) ii = 1
603  if ( pole_ ) then
604  jj = max( jj, 0 )
605  jj = min( jj, ja_ref )
606  else
607  jj = max( jj, 1 )
608  jj = min( jj, ja_ref-1 )
609  end if
610  else
611  if ( ii == 0 .and. jj == 0 ) then
612  ii = ia_ref - 1
613  jj = ja_ref - 1
614  end if
615  if ( ii == 0 .and. jj == ja_ref ) then
616  ii = ia_ref - 1
617  jj = 1
618  end if
619  if ( ii == ia_ref .and. jj == 0 ) then
620  ii = 1
621  jj = ja_ref - 1
622  end if
623  if ( ii == ia_ref .and. jj == ja_ref ) then
624  ii = 1
625  jj = 1
626  end if
627  ii = max( min( ii, ia_ref-1 ), 1 )
628  jj = max( min( jj, ja_ref-1 ), 1 )
629  end if
630  end if
631  end do
632  if ( ite == ite_max+1 ) then
633  log_error("INTERP_factor2d_linear_xy",*) 'iteration max has been reached', i, j, x(i), y(j)
634  log_error_cont(*) minval(x_ref), maxval(x_ref), minval(y_ref), maxval(y_ref)
635  log_error_cont(*) x_ref(1,1), x_ref(ia_ref,1),x_ref(1,ja_ref)
636  log_error_cont(*) y_ref(1,1), y_ref(ia_ref,1),y_ref(1,ja_ref)
637  idx_i(i,j,:) = 1
638  idx_j(i,j,:) = 1
639  hfact(i,j,:) = 0.0_rp
640  call prc_abort
641  end if
642 
643  if ( i==1 ) then
644  ii0 = ii
645  jj0 = jj
646  end if
647 
648  if ( error ) exit
649 
650  call sort_exec( 4, & ! [IN]
651  hfact(i,j,:), idx_i(i,j,:), idx_j(i,j,:), & ! [INOUT]
652  reverse = .true. ) ! [IN]
653 
654  end do
655  end do
656 
657  if ( error ) call prc_abort
658 
659  call prof_rapend ('INTERP_fact',3)
660 
661  return

References scale_const::const_eps, scale_const::const_undef, interp_bilinear_inv(), scale_prc::prc_abort(), scale_prof::prof_rapend(), and scale_prof::prof_rapstart().

Referenced by mod_cnv2d::cnv2d_exec(), mod_cnvuser::cnvuser_setup(), and interp_factor3d_linear_xy().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ interp_factor2d_weight()

subroutine, public scale_interp::interp_factor2d_weight ( integer, intent(in)  npoints,
integer, intent(in)  IA_ref,
integer, intent(in)  JA_ref,
integer, intent(in)  IA,
integer, intent(in)  JA,
real(rp), dimension(ia_ref,ja_ref), intent(in)  lon_ref,
real(rp), dimension(ia_ref,ja_ref), intent(in)  lat_ref,
real(rp), dimension (ia,ja), intent(in)  lon,
real(rp), dimension (ia,ja), intent(in)  lat,
integer, dimension(ia,ja,npoints), intent(out)  idx_i,
integer, dimension(ia,ja,npoints), intent(out)  idx_j,
real(rp), dimension(ia,ja,npoints), intent(out)  hfact,
real(rp), intent(in), optional  search_limit,
logical, intent(in), optional  latlon_structure,
real(rp), dimension(ia_ref), intent(in), optional  lon_1d,
real(rp), dimension(ja_ref), intent(in), optional  lat_1d,
integer, intent(in), optional  weight_order 
)

Definition at line 677 of file scale_interp.F90.

677  use scale_prc, only: &
678  prc_abort
679  implicit none
680  integer, intent(in) :: npoints ! number of interpolation point for horizontal
681  integer, intent(in) :: IA_ref ! number of x-direction (reference)
682  integer, intent(in) :: JA_ref ! number of y-direction (reference)
683  integer, intent(in) :: IA ! number of x-direction (target)
684  integer, intent(in) :: JA ! number of y-direction (target)
685  real(RP), intent(in) :: lon_ref(IA_ref,JA_ref) ! longitude [rad] (reference)
686  real(RP), intent(in) :: lat_ref(IA_ref,JA_ref) ! latitude [rad] (reference)
687  real(RP), intent(in) :: lon (IA,JA) ! longitude [rad] (target)
688  real(RP), intent(in) :: lat (IA,JA) ! latitude [rad] (target)
689  integer, intent(out) :: idx_i(IA,JA,npoints) ! i-index in reference (target)
690  integer, intent(out) :: idx_j(IA,JA,npoints) ! j-index in reference (target)
691  real(RP), intent(out) :: hfact(IA,JA,npoints) ! horizontal interp factor (target)
692 
693  real(RP), intent(in), optional :: search_limit
694  logical, intent(in), optional :: latlon_structure
695  real(RP), intent(in), optional :: lon_1d(IA_ref), lat_1d(JA_ref)
696  integer, intent(in), optional :: weight_order
697 
698  logical :: ll_struct_
699 
700  real(RP) :: lon_min, lon_max
701  real(RP) :: lat_min, lat_max
702  real(RP) :: dlon, dlat
703 
704  ! for structure grid
705  integer :: psizex, psizey
706  real(RP) :: lon0, lat0
707  integer, allocatable :: i0(:), i1(:), j0(:), j1(:)
708 
709  ! for unstructure grid
710  integer :: nsize, psize, nidx_max
711  integer, allocatable :: idx_blk(:,:,:), nidx(:,:)
712  integer :: idx_ref(npoints)
713 
714 
715  integer :: i, j, ii, jj, n
716  !---------------------------------------------------------------------------
717 
718  call prof_rapstart('INTERP_fact',3)
719 
720  if ( present(latlon_structure) ) then
721  ll_struct_ = latlon_structure
722  else
723  ll_struct_ = .false.
724  end if
725 
726  hfact(:,:,:) = 0.0_rp
727 
728 
729  if ( ll_struct_ ) then
730 
731  if ( ia_ref > 10 ) then
732  psizex = int( 2.0_rp*sqrt(real(ia_ref+1,rp)) )
733  else
734  psizex = 1
735  end if
736  if ( ja_ref > 10 ) then
737  psizey = int( 2.0_rp*sqrt(real(ja_ref+1,rp)) )
738  else
739  psizey = 1
740  end if
741 
742  allocate( i0(psizex), i1(psizex) )
743  allocate( j0(psizey), j1(psizey) )
744 
745  dlon = ( lon_max - lon_min ) / psizex
746  dlat = ( lat_max - lat_min ) / psizey
747 
748  do ii = 1, psizex
749  lon0 = lon_min + dlon * (ii-1)
750  do i = 1, ia_ref
751  if ( lon_1d(i) >= lon0 ) then
752  i0(ii) = i
753  exit
754  end if
755  end do
756  end do
757  do ii = 1, psizex-1
758  i1(ii) = i0(ii+1) - 1
759  end do
760  i1(psizex) = ia_ref
761 
762  do jj = 1, psizey
763  lat0 = lat_min + dlat * (jj-1)
764  do j = 1, ja_ref
765  if ( lat_1d(j) >= lat0 ) then
766  j0(jj) = j
767  exit
768  end if
769  end do
770  end do
771  do jj = 1, psizey-1
772  j1(jj) = j0(jj+1) - 1
773  end do
774  j1(psizey) = ja_ref
775 
776  !$omp parallel do OMP_SCHEDULE_ collapse(2)
777  do j = 1, ja
778  do i = 1, ia
779  ! main search
780  call interp_search_horiz_struct( npoints, & ! [IN]
781  psizex, psizey, & ! [IN]
782  ia_ref, ja_ref, & ! [IN]
783  lon_1d(:), lat_1d(:), & ! [IN]
784  lon_min, lat_min, & ! [IN]
785  dlon, dlat, & ! [IN]
786  i0(:), i1(:), j0(:), j1(:), & ! [IN]
787  lon(i,j), lat(i,j), & ! [IN]
788  idx_i(i,j,:), idx_j(i,j,:), & ! [OUT]
789  hfact(i,j,:), & ! [OUT]
790  search_limit = search_limit, & ! [IN]
791  weight_order = weight_order ) ! [IN]
792  enddo
793  enddo
794 
795  deallocate( i0, i1, j0, j1 )
796 
797  else
798 
799  nsize = ia_ref * ja_ref
800  if ( nsize > 100 ) then
801  psize = int( sqrt(2.0_rp*sqrt(real(nsize,rp))) )
802  nidx_max = nsize / psize * interp_buffer_size_fact
803  else
804  psize = 1
805  nidx_max = nsize
806  end if
807 
808  allocate(idx_blk(nidx_max,psize,psize))
809  allocate(nidx( psize,psize))
810 
811  call interp_div_block(nsize, psize, nidx_max, & ! [IN]
812  lon_ref(:,:), lat_ref(:,:), & ! [IN]
813  idx_blk(:,:,:), nidx(:,:), & ! [OUT]
814  lon_min, lon_max, & ! [OUT]
815  lat_min, lat_max, & ! [OUT]
816  dlon, dlat ) ! [OUT]
817 
818  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
819  !$omp private(idx_ref)
820  do j = 1, ja
821  do i = 1, ia
822  ! main search
823  call interp_search_horiz( npoints, & ! [IN]
824  nsize, & ! [IN]
825  lon_ref(:,:), lat_ref(:,:), & ! [IN]
826  lon_min, lon_max, & ! [IN]
827  lat_min, lat_max, & ! [IN]
828  psize, nidx_max, & ! [IN]
829  dlon, dlat, & ! [IN]
830  idx_blk(:,:,:), nidx(:,:), & ! [IN]
831  lon(i,j), lat(i,j), & ! [IN]
832  idx_ref(:), & ! [OUT]
833  hfact(i,j,:), & ! [OUT]
834  search_limit = search_limit, & ! [IN]
835  weight_order = weight_order ) ! [IN]
836  do n = 1, npoints
837  idx_i(i,j,n) = mod(idx_ref(n) - 1, ia_ref) + 1
838  idx_j(i,j,n) = ( idx_ref(n) - 1 ) / ia_ref + 1
839  end do
840  enddo
841  enddo
842 
843  deallocate(idx_blk, nidx)
844 
845  end if
846 
847  call prof_rapend ('INTERP_fact',3)
848 
849  return

References interp_div_block(), interp_search_horiz_struct(), scale_prc::prc_abort(), scale_prof::prof_rapend(), and scale_prof::prof_rapstart().

Referenced by mod_cnv2d::cnv2d_exec(), and mod_cnvuser::cnvuser_setup().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ interp_factor3d_linear_latlon()

subroutine, public scale_interp::interp_factor3d_linear_latlon ( integer, intent(in)  KA_ref,
integer, intent(in)  KS_ref,
integer, intent(in)  KE_ref,
integer, intent(in)  IA_ref,
integer, intent(in)  JA_ref,
integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
integer, intent(in)  IA,
integer, intent(in)  JA,
real(rp), dimension(ia_ref), intent(in)  lon_ref,
real(rp), dimension(ja_ref), intent(in)  lat_ref,
real(rp), dimension(ka_ref,ia_ref,ja_ref), intent(in)  hgt_ref,
real(rp), dimension(ia,ja), intent(in)  lon,
real(rp), dimension(ia,ja), intent(in)  lat,
real(rp), dimension(ka,ia,ja), intent(in)  hgt,
integer, dimension(ia,ja,4), intent(out)  idx_i,
integer, dimension(ia,ja,4), intent(out)  idx_j,
real(rp), dimension(ia,ja,4), intent(out)  hfact,
integer, dimension(ka,2,ia,ja,4), intent(out)  idx_k,
real(rp), dimension(ka, ia,ja,4), intent(out)  vfact,
logical, intent(in), optional  flag_extrap 
)

Definition at line 869 of file scale_interp.F90.

869  implicit none
870 
871  integer, intent(in) :: KA_ref, KS_ref, KE_ref ! number of z-direction (reference)
872  integer, intent(in) :: IA_ref ! number of x-direction (reference)
873  integer, intent(in) :: JA_ref ! number of y-direction (reference)
874  integer, intent(in) :: KA, KS, KE ! number of z-direction (target)
875  integer, intent(in) :: IA ! number of x-direction (target)
876  integer, intent(in) :: JA ! number of y-direction (target)
877  real(RP), intent(in) :: lon_ref(IA_ref) ! longitude (reference)
878  real(RP), intent(in) :: lat_ref(JA_ref) ! latitude (reference)
879  real(RP), intent(in) :: hgt_ref(KA_ref,IA_ref,JA_ref) ! height [m] (reference)
880  real(RP), intent(in) :: lon(IA,JA) ! longitude (target)
881  real(RP), intent(in) :: lat(IA,JA) ! latitude (target)
882  real(RP), intent(in) :: hgt(KA,IA,JA) ! longitude [m] (target)
883 
884  integer, intent(out) :: idx_i(IA,JA,4) ! i-index in reference (target)
885  integer, intent(out) :: idx_j(IA,JA,4) ! j-index in reference (target)
886  real(RP), intent(out) :: hfact(IA,JA,4) ! horizontal interp factor (target)
887  integer, intent(out) :: idx_k(KA,2,IA,JA,4) ! k-index in reference (target)
888  real(RP), intent(out) :: vfact(KA, IA,JA,4) ! vertical interp factor (target)
889 
890  logical, intent(in), optional :: flag_extrap ! when true, vertical extrapolation will be executed (just copy)
891 
892  integer :: i, j, ii, jj, n
893 
894  call interp_factor2d_linear_latlon( &
895  ia_ref, ja_ref, & ! [IN]
896  ia, ja, & ! [IN]
897  lon_ref(:), lat_ref(:), & ! [IN]
898  lon(:,:), lat(:,:), & ! [IN]
899  idx_i(:,:,:), idx_j(:,:,:), hfact(:,:,:) ) ! [OUT]
900 
901  call prof_rapstart('INTERP_fact',3)
902 
903  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
904  !$omp private(ii,jj)
905  do j = 1, ja
906  do i = 1, ia
907  do n = 1, 4
908  ii = idx_i(i,j,n)
909  jj = idx_j(i,j,n)
910  call interp_factor1d( ka_ref, ks_ref, ke_ref, & ! [IN]
911  ka, ks, ke, & ! [IN]
912  hgt_ref(:,ii,jj), & ! [IN]
913  hgt(:,i,j), & ! [IN]
914  idx_k(:,:,i,j,n), & ! [OUT]
915  vfact(:, i,j,n), & ! [OUT]
916  flag_extrap = flag_extrap ) ! [IN, optional]
917 
918  enddo
919  enddo
920  enddo
921 
922  call prof_rapend('INTERP_fact',3)
923 
924  return

References interp_factor1d(), interp_factor2d_linear_latlon(), scale_prof::prof_rapend(), and scale_prof::prof_rapstart().

Here is the call graph for this function:

◆ interp_factor3d_linear_xy()

subroutine, public scale_interp::interp_factor3d_linear_xy ( integer, intent(in)  KA_ref,
integer, intent(in)  KS_ref,
integer, intent(in)  KE_ref,
integer, intent(in)  IA_ref,
integer, intent(in)  JA_ref,
integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
integer, intent(in)  IA,
integer, intent(in)  JA,
real(rp), dimension(ia_ref,ja_ref), intent(in)  x_ref,
real(rp), dimension(ia_ref,ja_ref), intent(in)  y_ref,
real(rp), dimension(ka_ref,ia_ref,ja_ref), intent(in)  hgt_ref,
real(rp), dimension (ia), intent(in)  x,
real(rp), dimension (ja), intent(in)  y,
real(rp), dimension(ka,ia,ja), intent(in)  hgt,
integer, dimension(ia,ja,4), intent(out)  idx_i,
integer, dimension(ia,ja,4), intent(out)  idx_j,
real(rp), dimension(ia,ja,4), intent(out)  hfact,
integer, dimension(ka,2,ia,ja,4), intent(out)  idx_k,
real(rp), dimension(ka, ia,ja,4), intent(out)  vfact,
logical, intent(in), optional  flag_extrap,
logical, intent(in), optional  zonal,
logical, intent(in), optional  pole 
)

Definition at line 945 of file scale_interp.F90.

945  implicit none
946  integer, intent(in) :: KA_ref, KS_ref, KE_ref ! number of z-direction (reference)
947  integer, intent(in) :: IA_ref ! number of x-direction (reference)
948  integer, intent(in) :: JA_ref ! number of y-direction (reference)
949  integer, intent(in) :: KA, KS, KE ! number of z-direction (target)
950  integer, intent(in) :: IA ! number of x-direction (target)
951  integer, intent(in) :: JA ! number of y-direction (target)
952  real(RP), intent(in) :: x_ref(IA_ref,JA_ref) ! x point (reference)
953  real(RP), intent(in) :: y_ref(IA_ref,JA_ref) ! y point (reference)
954  real(RP), intent(in) :: hgt_ref(KA_ref,IA_ref,JA_ref) ! height [m] (reference)
955  real(RP), intent(in) :: x (IA) ! x point (target)
956  real(RP), intent(in) :: y (JA) ! y point (target)
957  real(RP), intent(in) :: hgt(KA,IA,JA) ! longitude [m] (target)
958 
959  integer, intent(out) :: idx_i(IA,JA,4) ! i-index in reference
960  integer, intent(out) :: idx_j(IA,JA,4) ! j-index in reference
961  real(RP), intent(out) :: hfact(IA,JA,4) ! horizontal interp factor
962  integer, intent(out) :: idx_k(KA,2,IA,JA,4) ! k-index in reference
963  real(RP), intent(out) :: vfact(KA, IA,JA,4) ! vertical interp factor
964 
965  logical, intent(in), optional :: flag_extrap ! when true, vertical extrapolation will be executed (just copy)
966  logical, intent(in), optional :: zonal
967  logical, intent(in), optional :: pole
968 
969  integer :: i, j, ii, jj, n
970 
971  call interp_factor2d_linear_xy( &
972  ia_ref, ja_ref, & ! [IN]
973  ia, ja, & ! [IN]
974  x_ref(:,:), y_ref(:,:), & ! [IN]
975  x(:), y(:), & ! [IN]
976  idx_i(:,:,:), idx_j(:,:,:), hfact(:,:,:), & ! [OUT]
977  zonal = zonal, pole = pole ) ! [IN]
978 
979  call prof_rapstart('INTERP_fact',3)
980 
981  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
982  !$omp private(ii,jj)
983  do j = 1, ja
984  do i = 1, ia
985  do n = 1, 4
986  ii = idx_i(i,j,n)
987  jj = idx_j(i,j,n)
988  call interp_factor1d( ka_ref, ks_ref, ke_ref, & ! [IN]
989  ka, ks, ke, & ! [IN]
990  hgt_ref(:,ii,jj), & ! [IN]
991  hgt(:,i,j), & ! [IN]
992  idx_k(:,:,i,j,n), & ! [OUT]
993  vfact(: ,i,j,n), & ! [OUT]
994  flag_extrap = flag_extrap ) ! [IN, optional]
995 
996  enddo
997  enddo
998  enddo
999 
1000  call prof_rapend('INTERP_fact',3)
1001 
1002  return

References interp_factor1d(), interp_factor2d_linear_xy(), scale_prof::prof_rapend(), and scale_prof::prof_rapstart().

Referenced by mod_cnvuser::cnvuser_setup().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ interp_factor3d_weight()

subroutine, public scale_interp::interp_factor3d_weight ( integer, intent(in)  npoints,
integer, intent(in)  KA_ref,
integer, intent(in)  KS_ref,
integer, intent(in)  KE_ref,
integer, intent(in)  IA_ref,
integer, intent(in)  JA_ref,
integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
integer, intent(in)  IA,
integer, intent(in)  JA,
real(rp), dimension(ia_ref,ja_ref), intent(in)  lon_ref,
real(rp), dimension(ia_ref,ja_ref), intent(in)  lat_ref,
real(rp), dimension(ka_ref,ia_ref,ja_ref), intent(in)  hgt_ref,
real(rp), dimension (ia,ja), intent(in)  lon,
real(rp), dimension (ia,ja), intent(in)  lat,
real(rp), dimension (ka,ia,ja), intent(in)  hgt,
integer, dimension(ia,ja,npoints), intent(out)  idx_i,
integer, dimension(ia,ja,npoints), intent(out)  idx_j,
real(rp), dimension(ia,ja,npoints), intent(out)  hfact,
integer, dimension(ka,2,ia,ja,npoints), intent(out)  idx_k,
real(rp), dimension(ka, ia,ja,npoints), intent(out)  vfact,
logical, intent(in), optional  flag_extrap 
)

Definition at line 1022 of file scale_interp.F90.

1022  use scale_prc, only: &
1023  prc_abort
1024  implicit none
1025  integer, intent(in) :: npoints ! number of interpolation point for horizontal
1026  integer, intent(in) :: KA_ref, KS_ref, KE_ref ! number of z-direction (reference)
1027  integer, intent(in) :: IA_ref ! number of x-direction (reference)
1028  integer, intent(in) :: JA_ref ! number of y-direction (reference)
1029  integer, intent(in) :: KA, KS, KE ! number of z-direction (target)
1030  integer, intent(in) :: IA ! number of x-direction (target)
1031  integer, intent(in) :: JA ! number of y-direction (target)
1032 
1033  real(RP), intent(in) :: lon_ref(IA_ref,JA_ref) ! longitude [rad] (reference)
1034  real(RP), intent(in) :: lat_ref(IA_ref,JA_ref) ! latitude [rad] (reference)
1035  real(RP), intent(in) :: hgt_ref(KA_ref,IA_ref,JA_ref) ! height [m] (reference)
1036  real(RP), intent(in) :: lon (IA,JA) ! longitude [rad] (target)
1037  real(RP), intent(in) :: lat (IA,JA) ! latitude [rad] (target)
1038  real(RP), intent(in) :: hgt (KA,IA,JA) ! longitude [m] (target)
1039 
1040  integer, intent(out) :: idx_i(IA,JA,npoints) ! i-index in reference
1041  integer, intent(out) :: idx_j(IA,JA,npoints) ! j-index in reference
1042  real(RP), intent(out) :: hfact(IA,JA,npoints) ! horizontal interp factor
1043  integer, intent(out) :: idx_k(KA,2,IA,JA,npoints) ! k-index in reference
1044  real(RP), intent(out) :: vfact(KA, IA,JA,npoints) ! vertical interp factor
1045 
1046  logical, intent(in), optional :: flag_extrap ! when true, vertical extrapolation will be executed (just copy)
1047 
1048  integer :: nsize, psize, nidx_max
1049  integer, allocatable :: idx_blk(:,:,:), nidx(:,:)
1050  real(RP) :: lon_min, lon_max
1051  real(RP) :: lat_min, lat_max
1052  real(RP) :: dlon, dlat
1053  integer :: idx_ref(npoints)
1054 
1055  integer :: i, j, ii, jj, n
1056  !---------------------------------------------------------------------------
1057 
1058  call prof_rapstart('INTERP_fact',3)
1059 
1060  nsize = ia_ref * ja_ref
1061  if ( nsize > 100 ) then
1062  psize = int( sqrt(2.0_rp*sqrt(real(nsize,rp))) )
1063  nidx_max = nsize / psize * interp_buffer_size_fact
1064  else
1065  psize = 1
1066  nidx_max = nsize
1067  end if
1068 
1069  allocate(idx_blk(nidx_max,psize,psize))
1070  allocate(nidx( psize,psize))
1071 
1072  call interp_div_block(nsize, psize, nidx_max, & ! [IN]
1073  lon_ref(:,:), lat_ref(:,:), & ! [IN]
1074  idx_blk(:,:,:), nidx(:,:), & ! [OUT]
1075  lon_min, lon_max, & ! [OUT]
1076  lat_min, lat_max, & ! [OUT]
1077  dlon, dlat ) ! [OUT]
1078 
1079  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
1080  !$omp private(ii,jj,idx_ref)
1081  do j = 1, ja
1082  do i = 1, ia
1083 
1084  ! main search
1085  call interp_search_horiz( npoints, & ! [IN]
1086  nsize, & ! [IN]
1087  lon_ref(:,:), lat_ref(:,:), & ! [IN]
1088  lon_min, lon_max, & ! [IN]
1089  lat_min, lat_max, & ! [IN]
1090  psize, nidx_max, & ! [IN]
1091  dlon, dlat, & ! [IN]
1092  idx_blk(:,:,:), nidx(:,:), & ! [IN]
1093  lon(i,j), lat(i,j), & ! [IN]
1094  idx_ref(:), hfact(i,j,:) ) ! [OUT]
1095 
1096  do n = 1, npoints
1097  ii = mod(idx_ref(n) - 1, ia_ref) + 1
1098  jj = ( idx_ref(n) - 1 ) / ia_ref + 1
1099  idx_i(i,j,n) = ii
1100  idx_j(i,j,n) = jj
1101  call interp_factor1d( ka_ref, ks_ref, ke_ref, & ! [IN]
1102  ka, ks, ke, & ! [IN]
1103  hgt_ref(:,ii,jj), & ! [IN]
1104  hgt(:,i,j), & ! [IN]
1105  idx_k(:,:,i,j,n), & ! [OUT]
1106  vfact(: ,i,j,n), & ! [OUT]
1107  flag_extrap = flag_extrap ) ! [IN, optional]
1108 
1109  enddo
1110  enddo
1111  enddo
1112 
1113  deallocate(idx_blk, nidx)
1114 
1115  call prof_rapend ('INTERP_fact',3)
1116 
1117  return

References interp_div_block(), interp_factor1d(), scale_prc::prc_abort(), scale_prof::prof_rapend(), and scale_prof::prof_rapstart().

Referenced by mod_cnvuser::cnvuser_setup().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ interp_interp1d()

subroutine, public scale_interp::interp_interp1d ( integer, intent(in)  KA_ref,
integer, intent(in)  KS_ref,
integer, intent(in)  KE_ref,
integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
integer, dimension(ka,2), intent(in)  idx_k,
real(rp), dimension(ka ), intent(in)  vfact,
real(rp), dimension(ka_ref), intent(in)  hgt_ref,
real(rp), dimension (ka), intent(in)  hgt,
real(rp), dimension(ka_ref), intent(in), target  val_ref,
real(rp), dimension (ka), intent(out)  val,
logical, intent(in), optional  spline,
logical, intent(in), optional  logwgt 
)

Definition at line 1134 of file scale_interp.F90.

1134  use scale_const, only: &
1135  undef => const_undef, &
1136  eps => const_eps
1137  implicit none
1138  integer, intent(in) :: KA_ref, KS_ref, KE_ref ! number of z-direction (reference)
1139  integer, intent(in) :: KA, KS, KE ! number of z-direction (target)
1140 
1141  integer, intent(in) :: idx_k(KA,2) ! k-index in reference
1142  real(RP), intent(in) :: vfact(KA ) ! vertical interp factor
1143  real(RP), intent(in) :: hgt_ref(KA_ref) ! height (reference)
1144  real(RP), intent(in) :: hgt (KA) ! height (target)
1145  real(RP), intent(in), target :: val_ref(KA_ref) ! value (reference)
1146 
1147  real(RP), intent(out) :: val (KA) ! value (target)
1148 
1149  logical, intent(in), optional :: spline
1150  logical, intent(in), optional :: logwgt
1151 
1152  logical :: spline_
1153  logical :: logwgt_
1154 
1155  real(RP), pointer :: work(:)
1156 
1157  integer :: idx (KA_ref)
1158  integer :: idx_r(KA_ref)
1159  real(RP) :: FDZ (KA_ref)
1160  real(RP) :: U (KA_ref)
1161 
1162  integer :: kmax
1163  integer :: k
1164  !---------------------------------------------------------------------------
1165 
1166  call prof_rapstart('INTERP_interp',3)
1167 
1168  spline_ = interp_use_spline_vert
1169  if ( present(spline) ) then
1170  spline_ = spline
1171  end if
1172 
1173  logwgt_ = .false.
1174  if ( present(logwgt) ) then
1175  logwgt_ = logwgt
1176  endif
1177 
1178  if ( logwgt_ ) then
1179  allocate( work(ka_ref) )
1180  do k = ks_ref, ke_ref
1181  if ( abs( val_ref(k) - undef ) < eps ) then
1182  work(k) = undef
1183  else
1184  work(k) = log( val_ref(k) )
1185  end if
1186  enddo
1187  else
1188  work => val_ref
1189  endif
1190 
1191  call spline_coef( ka_ref, ks_ref, ke_ref, &
1192  hgt_ref(:), work(:), & ! (in)
1193  spline_, & ! (in)
1194  kmax, & ! (out)
1195  idx(:), idx_r(:), & ! (out)
1196  u(:), fdz(:) ) ! (out)
1197 
1198  call spline_exec( ka_ref, kmax, ka, ks, ke, &
1199  idx_k(:,:), vfact(:), & ! (in)
1200  hgt_ref(:), hgt(:), & ! (in)
1201  work(:), & ! (in)
1202  idx(:), idx_r(:), & ! (in)
1203  u(:), fdz(:), & ! (in)
1204  val(:) ) ! (out)
1205 
1206  if ( logwgt_ ) then
1207  deallocate( work )
1208  do k = ks, ke
1209  if ( abs( val(k) - undef ) > eps ) then
1210  val(k) = exp( val(k) )
1211  endif
1212  end do
1213  endif
1214 
1215  call prof_rapend ('INTERP_interp',3)
1216 
1217  return

References scale_const::const_eps, scale_const::const_undef, scale_prof::prof_rapend(), scale_prof::prof_rapstart(), and spline_exec().

Referenced by scale_interp_vert::interp_vert_xi2p(), scale_interp_vert::interp_vert_xi2z(), scale_interp_vert::interp_vert_xih2p(), scale_interp_vert::interp_vert_xih2zh(), scale_interp_vert::interp_vert_z2xi(), and scale_interp_vert::interp_vert_zh2xih().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ interp_interp2d()

subroutine, public scale_interp::interp_interp2d ( integer, intent(in)  npoints,
integer, intent(in)  IA_ref,
integer, intent(in)  JA_ref,
integer, intent(in)  IA,
integer, intent(in)  JA,
integer, dimension (ia,ja,npoints), intent(in)  idx_i,
integer, dimension (ia,ja,npoints), intent(in)  idx_j,
real(rp), dimension (ia,ja,npoints), intent(in)  hfact,
real(rp), dimension(ia_ref,ja_ref), intent(in)  val_ref,
real(rp), dimension (ia,ja), intent(out)  val,
real(rp), intent(in), optional  threshold_undef,
real(rp), dimension(ia,ja), intent(out), optional  wsum,
real(rp), dimension(ia,ja), intent(out), optional  val2 
)

Definition at line 1232 of file scale_interp.F90.

1232  use scale_const, only: &
1233  undef => const_undef, &
1234  eps => const_eps
1235  implicit none
1236 
1237  integer, intent(in) :: npoints ! number of interpolation point for horizontal
1238  integer, intent(in) :: IA_ref ! number of x-direction (reference)
1239  integer, intent(in) :: JA_ref ! number of y-direction (reference)
1240  integer, intent(in) :: IA ! number of x-direction (target)
1241  integer, intent(in) :: JA ! number of y-direction (target)
1242  integer, intent(in) :: idx_i (IA,JA,npoints) ! i-index in reference
1243  integer, intent(in) :: idx_j (IA,JA,npoints) ! j-index in reference
1244  real(RP), intent(in) :: hfact (IA,JA,npoints) ! horizontal interp factor
1245  real(RP), intent(in) :: val_ref(IA_ref,JA_ref) ! value (reference)
1246 
1247  real(RP), intent(out) :: val (IA,JA) ! value (target)
1248 
1249  real(RP), intent(in), optional :: threshold_undef
1250  real(RP), intent(out), optional :: wsum(IA,JA)
1251  real(RP), intent(out), optional :: val2(IA,JA)
1252 
1253  real(RP) :: th_undef
1254 
1255  real(RP) :: fact, valn, f, w
1256  real(RP) :: sw
1257  logical :: lval2
1258 
1259  integer :: i, j, n
1260  !---------------------------------------------------------------------------
1261 
1262  call prof_rapstart('INTERP_interp',3)
1263 
1264  if ( present(threshold_undef) ) then
1265  th_undef = threshold_undef
1266  else
1267  th_undef = interp_threshold_undef
1268  end if
1269  th_undef = min( max( th_undef, eps * 2.0_rp ), 1.0_rp - eps * 2.0_rp )
1270 
1271  lval2 = present(wsum) .and. present(val2)
1272 
1273  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
1274  !$omp private(fact,valn,f,w,sw)
1275 !OCL PREFETCH
1276  do j = 1, ja
1277  do i = 1, ia
1278  fact = 0.0_rp
1279  valn = 0.0_rp
1280  do n = 1, npoints
1281  f = hfact(i,j,n)
1282  w = val_ref(idx_i(i,j,n),idx_j(i,j,n))
1283  if ( f > eps .and. abs( w - undef ) > eps ) then
1284  fact = fact + f
1285  valn = valn + f * w
1286  end if
1287  end do
1288  sw = 0.5_rp - sign( 0.5_rp, fact - th_undef + eps ) ! 1.0 when fact < threshold
1289  val(i,j) = valn / ( fact + sw ) * ( 1.0_rp - sw ) + undef * sw
1290  if ( lval2 ) then
1291  wsum(i,j) = fact
1292  sw = 0.5_rp - sign( 0.5_rp, fact - eps ) ! 1.0 when fact < 0.0
1293  val2(i,j) = valn / ( fact + sw ) * ( 1.0_rp - sw ) + undef * sw
1294  end if
1295  enddo
1296  enddo
1297 
1298  call prof_rapend ('INTERP_interp',3)
1299 
1300  return

References scale_const::const_eps, scale_const::const_undef, scale_prof::prof_rapend(), and scale_prof::prof_rapstart().

Referenced by scale_atmos_sfc_ch_rn222::atmos_sfc_ch_rn222_land_flux(), mod_cnv2d::cnv2d_exec(), mod_copytopo::copytopo(), mod_realinput::land_interporation(), mod_realinput::ocean_interporation(), and mod_realinput::realinput_surface().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ interp_interp3d()

subroutine, public scale_interp::interp_interp3d ( integer, intent(in)  npoints,
integer, intent(in)  KA_ref,
integer, intent(in)  KS_ref,
integer, intent(in)  KE_ref,
integer, intent(in)  IA_ref,
integer, intent(in)  JA_ref,
integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
integer, intent(in)  IA,
integer, intent(in)  JA,
integer, dimension (ia,ja,npoints), intent(in)  idx_i,
integer, dimension (ia,ja,npoints), intent(in)  idx_j,
real(rp), dimension (ia,ja,npoints), intent(in)  hfact,
integer, dimension (ka,2,ia,ja,npoints), intent(in)  idx_k,
real(rp), dimension (ka, ia,ja,npoints), intent(in)  vfact,
real(rp), dimension(ka_ref,ia_ref,ja_ref), intent(in)  hgt_ref,
real(rp), dimension (ka,ia,ja), intent(in)  hgt,
real(rp), dimension(ka_ref,ia_ref,ja_ref), intent(in), target  val_ref,
real(rp), dimension (ka,ia,ja), intent(out)  val,
logical, intent(in), optional  spline,
logical, intent(in), optional  logwgt,
real(rp), intent(in), optional  threshold_undef,
real(rp), dimension(ka,ia,ja), intent(out), optional  wsum,
real(rp), dimension(ka,ia,ja), intent(out), optional  val2 
)

Definition at line 1322 of file scale_interp.F90.

1322  use scale_const, only: &
1323  undef => const_undef, &
1324  eps => const_eps
1325  implicit none
1326  integer, intent(in) :: npoints ! number of interpolation point for horizontal
1327  integer, intent(in) :: KA_ref, KS_ref, KE_ref ! number of z-direction (reference)
1328  integer, intent(in) :: IA_ref ! number of x-direction (reference)
1329  integer, intent(in) :: JA_ref ! number of y-direction (reference)
1330  integer, intent(in) :: KA, KS, KE ! number of z-direction (target)
1331  integer, intent(in) :: IA ! number of x-direction (target)
1332  integer, intent(in) :: JA ! number of y-direction (target)
1333 
1334  integer, intent(in) :: idx_i (IA,JA,npoints) ! i-index in reference
1335  integer, intent(in) :: idx_j (IA,JA,npoints) ! j-index in reference
1336  real(RP), intent(in) :: hfact (IA,JA,npoints) ! horizontal interp factor
1337  integer, intent(in) :: idx_k (KA,2,IA,JA,npoints) ! k-index in reference
1338  real(RP), intent(in) :: vfact (KA, IA,JA,npoints) ! vertical interp factor
1339  real(RP), intent(in) :: hgt_ref(KA_ref,IA_ref,JA_ref) ! height (reference)
1340  real(RP), intent(in) :: hgt (KA,IA,JA) ! height (target)
1341  real(RP), intent(in), target :: val_ref(KA_ref,IA_ref,JA_ref) ! value (reference)
1342 
1343  real(RP), intent(out) :: val (KA,IA,JA) ! value (target)
1344 
1345  logical, intent(in), optional :: spline
1346  logical, intent(in), optional :: logwgt
1347  real(RP), intent(in), optional :: threshold_undef
1348  real(RP), intent(out), optional :: wsum(KA,IA,JA)
1349  real(RP), intent(out), optional :: val2(KA,IA,JA)
1350 
1351  real(RP) :: th_undef
1352  logical :: spline_
1353  logical :: logwgt_
1354 
1355  real(RP), pointer :: work(:,:,:)
1356 
1357  integer, allocatable :: kmax(:,:)
1358  integer, allocatable :: idx(:,:,:), idx_r(:,:,:)
1359  real(RP), allocatable :: U(:,:,:), FDZ(:,:,:)
1360 
1361  real(RP) :: valn
1362  real(RP) :: fact
1363  real(RP) :: f
1364  real(RP) :: sw
1365  real(RP) :: w(KA,npoints)
1366  logical :: lval2
1367 
1368  integer :: imin, imax
1369  integer :: jmin, jmax
1370  integer :: ii, jj
1371  integer :: k, i, j, n
1372  !---------------------------------------------------------------------------
1373 
1374  call prof_rapstart('INTERP_interp',3)
1375 
1376  if ( present(threshold_undef) ) then
1377  th_undef = threshold_undef
1378  else
1379  th_undef = interp_threshold_undef
1380  end if
1381  th_undef = min( max( th_undef, eps * 2.0_rp ), 1.0_rp - eps * 2.0_rp )
1382 
1383 
1384  spline_ = interp_use_spline_vert
1385  if ( present(spline) ) then
1386  spline_ = spline
1387  end if
1388 
1389  logwgt_ = .false.
1390  if ( present(logwgt) ) then
1391  logwgt_ = logwgt
1392  endif
1393 
1394 
1395  lval2 = present(wsum) .and. present(val2)
1396 
1397 
1398  imin = ia_ref
1399  jmin = ja_ref
1400  imax = 1
1401  jmax = 1
1402  !$omp parallel do OMP_SCHEDULE_ collapse(3) &
1403  !$omp reduction(min: imin,jmin) &
1404  !$omp reduction(max: imax,jmax)
1405  do n = 1, npoints
1406  do j = 1, ja
1407  do i = 1, ia
1408  imin = min(imin, idx_i(i,j,n))
1409  imax = max(imax, idx_i(i,j,n))
1410  jmin = min(jmin, idx_j(i,j,n))
1411  jmax = max(jmax, idx_j(i,j,n))
1412  end do
1413  end do
1414  end do
1415 
1416  allocate( kmax(imin:imax,jmin:jmax) )
1417  allocate( idx(ka_ref,imin:imax,jmin:jmax), idx_r(ka_ref,imin:imax,jmin:jmax) )
1418  allocate( u(ka_ref,imin:imax,jmin:jmax), fdz(ka_ref,imin:imax,jmin:jmax) )
1419 
1420  if ( logwgt_ ) then
1421  allocate( work(ka_ref,imin:imax,jmin:jmax) )
1422  !$omp parallel do OMP_SCHEDULE_ collapse(2)
1423  do j = jmin, jmax
1424  do i = imin, imax
1425  do k = ks_ref, ke_ref
1426  if ( abs( val_ref(k,i,j) - undef ) < eps ) then
1427  work(k,i,j) = undef
1428  else
1429  work(k,i,j) = log( val_ref(k,i,j) )
1430  end if
1431  enddo
1432  enddo
1433  enddo
1434  else
1435  work => val_ref
1436  endif
1437 
1438  !$omp parallel do OMP_SCHEDULE_ collapse(2)
1439  do j = jmin, jmax
1440  do i = imin, imax
1441  call spline_coef( ka_ref, ks_ref, ke_ref, &
1442  hgt_ref(:,i,j), work(:,i,j), & ! (in)
1443  spline_, & ! (in)
1444  kmax(i,j), & ! (out)
1445  idx(:,i,j), idx_r(:,i,j), & ! (out)
1446  u(:,i,j), fdz(:,i,j) ) ! (out)
1447  end do
1448  end do
1449 
1450  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
1451  !$omp private(valn,fact,w,f,ii,jj,sw)
1452  do j = 1, ja
1453  do i = 1, ia
1454 
1455  do n = 1, npoints
1456  if ( hfact(i,j,n) < eps ) exit
1457  ii = idx_i(i,j,n)
1458  jj = idx_j(i,j,n)
1459  call spline_exec( ka_ref, kmax(ii,jj), ka, ks, ke, &
1460  idx_k(:,:,i,j,n), vfact(:, i,j,n), & ! (in)
1461  hgt_ref(:,ii,jj), hgt(:,i,j), & ! (in)
1462  work(:,ii,jj), & ! (in)
1463  idx(:,ii,jj), idx_r(:,ii,jj), & ! (in)
1464  u(:,ii,jj), fdz(:,ii,jj), & ! (in)
1465  w(:,n) ) ! (out)
1466  end do
1467 
1468  do k = ks, ke
1469  fact = 0.0_rp
1470  valn = 0.0_rp
1471  do n = 1, npoints
1472  f = hfact(i,j,n)
1473  if ( f > eps .and. abs( w(k,n) - undef ) > eps ) then
1474  fact = fact + f
1475  valn = valn + f * w(k,n)
1476  endif
1477  enddo
1478  sw = 0.5_rp - sign( 0.5_rp, fact - th_undef ) ! 1.0 when fact < threshold
1479  val(k,i,j) = valn / ( fact + sw ) * ( 1.0_rp - sw ) + undef * sw
1480  if ( lval2 ) then
1481  wsum(k,i,j) = fact
1482  sw = 0.5_rp - sign( 0.5_rp, fact - eps ) ! 1.0 when fact < 0.0
1483  val2(k,i,j) = valn / ( fact + sw ) * ( 1.0_rp - sw ) + undef * sw
1484  end if
1485  enddo
1486 
1487  enddo
1488  enddo
1489 
1490  deallocate( kmax, idx, idx_r )
1491  deallocate( u, fdz )
1492 
1493  if ( logwgt_ ) then
1494  deallocate( work )
1495  !$omp parallel do OMP_SCHEDULE_ collapse(2)
1496  do j = 1, ja
1497  do i = 1, ia
1498  do k = ks, ke
1499  if ( abs( val(k,i,j) - undef ) > eps ) then
1500  val(k,i,j) = exp( val(k,i,j) )
1501  end if
1502  end do
1503  end do
1504  end do
1505  end if
1506 
1507  call prof_rapend ('INTERP_interp',3)
1508 
1509  return

References scale_const::const_eps, scale_const::const_radius, scale_const::const_undef, scale_prof::prof_rapend(), scale_prof::prof_rapstart(), and spline_exec().

Referenced by mod_cnvuser::cnvuser(), scale_comm_cartesc_nest::comm_cartesc_nest_intercomm_nestdown_3d(), mod_realinput::land_interporation(), and mod_realinput::realinput_surface().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ interp_search_horiz_struct()

subroutine scale_interp::interp_search_horiz_struct ( integer, intent(in)  npoints,
integer, intent(in)  psizex,
integer, intent(in)  psizey,
integer, intent(in)  IA_ref,
integer, intent(in)  JA_ref,
real(rp), dimension(ia_ref), intent(in)  lon_ref,
real(rp), dimension(ja_ref), intent(in)  lat_ref,
real(rp), intent(in)  lon_min,
real(rp), intent(in)  lat_min,
real(rp), intent(in)  dlon,
real(rp), intent(in)  dlat,
integer, dimension(psizex), intent(in)  i0,
integer, dimension(psizex), intent(in)  i1,
integer, dimension(psizey), intent(in)  j0,
integer, dimension(psizey), intent(in)  j1,
real(rp), intent(in)  lon,
real(rp), intent(in)  lat,
integer, dimension(npoints), intent(out)  idx_i,
integer, dimension(npoints), intent(out)  idx_j,
real(rp), dimension(npoints), intent(out)  hfact,
real(rp), intent(in), optional  search_limit,
integer, intent(in), optional  weight_order 
)

Definition at line 1688 of file scale_interp.F90.

1688  use scale_const, only: &
1689  eps => const_eps, &
1690  radius => const_radius
1691  implicit none
1692 
1693  integer, intent(in) :: npoints ! number of interpolation point for horizontal
1694  integer, intent(in) :: psizex ! number of block in x-direction (reference)
1695  integer, intent(in) :: psizey ! number of block in y-direction (reference)
1696  integer, intent(in) :: IA_ref ! number of grids in x-direction (reference)
1697  integer, intent(in) :: JA_ref ! number of grids in y-direction (reference)
1698  real(RP), intent(in) :: lon_ref(IA_ref) ! longitude [rad] (reference)
1699  real(RP), intent(in) :: lat_ref(JA_ref) ! latitude [rad] (reference)
1700  real(RP), intent(in) :: lon_min ! minimum longitude
1701  real(RP), intent(in) :: lat_min ! minimum latitude
1702  real(RP), intent(in) :: dlon
1703  real(RP), intent(in) :: dlat
1704  integer, intent(in) :: i0(psizex) ! start index in the block
1705  integer, intent(in) :: i1(psizex) ! end index in the block
1706  integer, intent(in) :: j0(psizey) ! start index in the block
1707  integer, intent(in) :: j1(psizey) ! start index in the block
1708  real(RP), intent(in) :: lon ! longitude [rad] (target)
1709  real(RP), intent(in) :: lat ! latitude [rad] (target)
1710  integer, intent(out) :: idx_i(npoints) ! x-index in reference (target)
1711  integer, intent(out) :: idx_j(npoints) ! y-index in reference (target)
1712  real(RP), intent(out) :: hfact(npoints) ! horizontal interp factor (target)
1713 
1714  real(RP), intent(in), optional :: search_limit
1715  integer, intent(in), optional :: weight_order
1716 
1717  real(RP) :: drad(npoints)
1718  real(RP) :: sum
1719  real(RP) :: search_limit_
1720  integer :: weight_order_
1721 
1722  real(RP) :: lon0, lon1, lat0, lat1
1723  real(RP) :: dlon_sl, dlat_sl
1724 
1725  integer :: i, j, n
1726  integer :: ii, jj, ii0, jj0
1727  !---------------------------------------------------------------------------
1728 
1729  if ( present(search_limit) ) then
1730  search_limit_ = search_limit
1731  else
1732  search_limit_ = interp_search_limit
1733  end if
1734  search_limit_ = search_limit_ / radius ! m to radian
1735 
1736  if ( present(weight_order) ) then
1737  weight_order_ = weight_order
1738  else
1739  weight_order_ = interp_weight_order
1740  end if
1741 
1742 
1743  drad(:) = large_number
1744  idx_i(:) = 1
1745  idx_j(:) = 1
1746 
1747  ! find k-nearest points in the nearest block
1748  ii0 = max( min( int(( lon - lon_min ) / dlon) + 1, psizex ), 1)
1749  jj0 = max( min( int(( lat - lat_min ) / dlat) + 1, psizey ), 1)
1750  do j = j0(jj0), j1(jj0)
1751  do i = i0(ii0), i1(ii0)
1752  call interp_insert_2d( npoints, &
1753  lon, lat, &
1754  lon_ref(i), lat_ref(j), & ! [IN]
1755  i, j, & ! [IN]
1756  drad(:), idx_i(:), idx_j(:) ) ! [INOUT]
1757  end do
1758  end do
1759 
1760  if ( abs(drad(1)) < eps ) then
1761  hfact(:) = 0.0_rp
1762  hfact(1) = 1.0_rp
1763 
1764  return
1765  else if ( drad(1) > search_limit_ ) then
1766  hfact(:) = 0.0_rp
1767  idx_i(:) = 1 ! dummy
1768  idx_j(:) = 1 ! dummy
1769 
1770  return
1771  end if
1772 
1773  dlon_sl = max(dlon * 0.5_rp, drad(npoints))
1774  dlat_sl = max(dlat * 0.5_rp, drad(npoints))
1775  do jj = 1, psizey
1776  lat0 = lat_min + dlat * (jj-1)
1777  lat1 = lat_min + dlat * jj
1778  if ( lat < lat0 - dlat_sl &
1779  .or. lat >= lat1 + dlat_sl ) cycle
1780  do ii = 1, psizex
1781  if ( ii==ii0 .and. jj==jj0 ) cycle
1782  lon0 = lon_min + dlon * (ii-1)
1783  lon1 = lon_min + dlon * ii
1784  if ( lon < lon0 - dlon_sl &
1785  .or. lon >= lon1 + dlon_sl ) cycle
1786  do j = j0(jj0), j1(jj0)
1787  do i = i0(ii0), i1(ii0)
1788  call interp_insert_2d( npoints, &
1789  lon, lat, &
1790  lon_ref(i), lat_ref(j), & ! [IN]
1791  i, j, & ! [IN]
1792  drad(:), idx_i(:), idx_j(:) ) ! [INOUT]
1793  end do
1794  end do
1795  dlon_sl = max(dlon * 0.5_rp, drad(npoints))
1796  dlat_sl = max(dlat * 0.5_rp, drad(npoints))
1797  end do
1798  end do
1799 
1800 
1801  ! factor = 1 / drad
1802  if ( weight_order_ < 0 ) then
1803  do n = 1, npoints
1804  hfact(n) = 1.0_rp / drad(n)**(-1.0_rp/weight_order_)
1805  enddo
1806  else
1807  do n = 1, npoints
1808  hfact(n) = 1.0_rp / drad(n)**weight_order_
1809  enddo
1810  end if
1811 
1812  ! ignore far point
1813  do n = 1, npoints
1814  if ( drad(n) >= search_limit_ ) then
1815  hfact(n) = 0.0_rp
1816  endif
1817  enddo
1818 
1819  ! normalize factor
1820  sum = 0.0_rp
1821  do n = 1, npoints
1822  sum = sum + hfact(n)
1823  enddo
1824 
1825  if ( sum > 0.0_rp ) then
1826  do n = 1, npoints
1827  hfact(n) = hfact(n) / sum
1828  enddo
1829  endif
1830 
1831  return

References scale_const::const_eps, scale_const::const_radius, scale_const::const_undef, and interp_insert_2d().

Referenced by interp_factor2d_weight().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ interp_insert_2d()

subroutine scale_interp::interp_insert_2d ( integer, intent(in)  npoints,
real(rp), intent(in)  lon,
real(rp), intent(in)  lat,
real(rp), intent(in)  lon_ref,
real(rp), intent(in)  lat_ref,
integer, intent(in)  i,
integer, intent(in)  j,
real(rp), dimension(npoints), intent(inout)  drad,
integer, dimension(npoints), intent(inout)  idx_i,
integer, dimension(npoints), intent(inout)  idx_j 
)

Definition at line 1881 of file scale_interp.F90.

1881  use scale_const, only: &
1882  undef => const_undef, &
1883  eps => const_eps
1884  use scale_sort, only: &
1885  sort_exec
1886 
1887  integer, intent(in) :: npoints
1888  real(RP), intent(in) :: lon, lat
1889  real(RP), intent(in) :: lon_ref, lat_ref
1890  integer, intent(in) :: i, j
1891 
1892  real(RP), intent(inout) :: drad(npoints)
1893  integer, intent(inout) :: idx_i(npoints)
1894  integer, intent(inout) :: idx_j(npoints)
1895 
1896  real(RP) :: dradian
1897 
1898  if ( abs( lon_ref - undef ) < eps ) return
1899 
1900  dradian = haversine( lon, lat, lon_ref, lat_ref )
1901 
1902  if ( dradian <= drad(npoints) ) then
1903  ! replace last(=longest) value
1904  drad(npoints) = dradian
1905  idx_i(npoints) = i
1906  idx_j(npoints) = j
1907 
1908  ! sort by ascending order
1909  call sort_exec( npoints, & ! [IN]
1910  drad(:), & ! [INOUT]
1911  idx_i(:), & ! [INOUT]
1912  idx_j(:) ) ! [INOUT]
1913  endif
1914 
1915  return

References scale_const::const_eps, and scale_const::const_undef.

Referenced by interp_search_horiz_struct().

Here is the caller graph for this function:

◆ interp_div_block()

subroutine scale_interp::interp_div_block ( integer, intent(in)  nsize,
integer, intent(in)  psize,
integer, intent(in)  nidx_max,
real(rp), dimension(nsize), intent(in)  lon_ref,
real(rp), dimension(nsize), intent(in)  lat_ref,
integer, dimension (nidx_max,psize,psize), intent(out)  idx,
integer, dimension(psize,psize), intent(out)  nidx,
real(rp), intent(out)  lon_min,
real(rp), intent(out)  lon_max,
real(rp), intent(out)  lat_min,
real(rp), intent(out)  lat_max,
real(rp), intent(out)  dlon,
real(rp), intent(out)  dlat 
)

Definition at line 1925 of file scale_interp.F90.

1925  use scale_const, only: &
1926  undef => const_undef, &
1927  eps => const_eps
1928  use scale_prc, only: &
1929  prc_abort
1930  integer, intent(in) :: nsize
1931  integer, intent(in) :: psize
1932  integer, intent(in) :: nidx_max
1933  real(RP), intent(in) :: lon_ref(nsize)
1934  real(RP), intent(in) :: lat_ref(nsize)
1935 
1936  integer, intent(out) :: idx (nidx_max,psize,psize)
1937  integer, intent(out) :: nidx(psize,psize)
1938  real(RP), intent(out) :: lon_min, lon_max
1939  real(RP), intent(out) :: lat_min, lat_max
1940  real(RP), intent(out) :: dlon, dlat
1941 
1942  integer :: i, ii, jj, n
1943 
1944  lon_min = minval(lon_ref(:), mask=abs(lon_ref-undef)>eps)
1945  lon_max = maxval(lon_ref(:), mask=abs(lon_ref-undef)>eps)
1946  lat_min = minval(lat_ref(:), mask=abs(lat_ref-undef)>eps)
1947  lat_max = maxval(lat_ref(:), mask=abs(lat_ref-undef)>eps)
1948 
1949  dlon = ( lon_max - lon_min ) / psize
1950  dlat = ( lat_max - lat_min ) / psize
1951 
1952  nidx(:,:) = 0
1953  do i = 1, nsize
1954  if ( abs( lon_ref(i) - undef ) < eps ) cycle
1955  ii = min(int((lon_ref(i) - lon_min) / dlon) + 1, psize)
1956  jj = min(int((lat_ref(i) - lat_min) / dlat) + 1, psize)
1957  n = nidx(ii,jj) + 1
1958  nidx(ii,jj) = n
1959  if ( n <= nidx_max ) idx(n,ii,jj) = i
1960  end do
1961 
1962  if ( maxval(nidx) > nidx_max ) then
1963  log_error("INTERP_search_horiz",*) 'Buffer size is not enough'
1964  log_error_cont(*) ' Use larger INTERP_buffer_size_fact: ', interp_buffer_size_fact * maxval(nidx)/nidx_max
1965  call prc_abort
1966  end if
1967 
1968  return

References scale_const::const_eps, scale_const::const_undef, and scale_prc::prc_abort().

Referenced by interp_factor2d_weight(), and interp_factor3d_weight().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ interp_bilinear_inv()

subroutine scale_interp::interp_bilinear_inv ( real(rp), intent(in)  x_ref0,
real(rp), intent(in)  x_ref1,
real(rp), intent(in)  x_ref2,
real(rp), intent(in)  x_ref3,
real(rp), intent(in)  y_ref0,
real(rp), intent(in)  y_ref1,
real(rp), intent(in)  y_ref2,
real(rp), intent(in)  y_ref3,
real(rp), intent(in)  x,
real(rp), intent(in)  y,
real(rp), intent(out)  u,
real(rp), intent(out)  v,
logical, intent(out)  error 
)

Definition at line 1978 of file scale_interp.F90.

1978  implicit none
1979  real(RP), intent(in) :: x_ref0, x_ref1, x_ref2, x_ref3
1980  real(RP), intent(in) :: y_ref0, y_ref1, y_ref2, y_ref3
1981  real(RP), intent(in) :: x, y
1982 
1983  real(RP), intent(out) :: u, v
1984  logical, intent(out) :: error
1985 
1986  real(RP) :: e_x, e_y
1987  real(RP) :: f_x, f_y
1988  real(RP) :: g_x, g_y
1989  real(RP) :: h_x, h_y
1990  real(RP) :: k0, k1, k2
1991  real(RP) :: w
1992  real(RP) :: sig
1993 
1994  e_x = x_ref1 - x_ref0
1995  e_y = y_ref1 - y_ref0
1996 
1997  f_x = x_ref3 - x_ref0
1998  f_y = y_ref3 - y_ref0
1999 
2000  g_x = x_ref0 - x_ref1 + x_ref2 - x_ref3
2001  g_y = y_ref0 - y_ref1 + y_ref2 - y_ref3
2002 
2003  h_x = x - x_ref0
2004  h_y = y - y_ref0
2005 
2006  k0 = cross(h_x, h_y, e_x, e_y)
2007  k1 = cross(e_x, e_y, f_x, f_y) + cross(h_x, h_y, g_x, g_y)
2008  k2 = cross(g_x, g_y, f_x, f_y)
2009 
2010  w = k1**2 - 4.0_rp * k0 * k2
2011  if ( w < 0.0_rp ) then
2012  log_error("INTERP_bilinear_inv",*) 'Outside of the quadrilateral'
2013  error = .true.
2014  return
2015  end if
2016 
2017  if ( abs(k1) < eps_bilinear ) then
2018  log_error("INTERP_bilinear_inv",*) 'Unexpected error occured', k1
2019  log_error_cont(*) x_ref0, x_ref1, x_ref2, x_ref3
2020  log_error_cont(*) y_ref0, y_ref1, y_ref2, y_ref3
2021  log_error_cont(*) x, y
2022  error = .true.
2023  return
2024  end if
2025  if ( abs(k2) < eps_bilinear * sqrt( (x_ref2-x_ref0)**2+(y_ref2-y_ref0)**2 ) ) then
2026  v = - k0 / k1
2027  else
2028  sig = sign( 1.0_rp, cross(x_ref1-x_ref0, y_ref1-y_ref0, x_ref3-x_ref0, y_ref3-y_ref0) )
2029  v = ( - k1 + sig * sqrt(w) ) / ( 2.0_rp * k2 )
2030  end if
2031  u = ( h_x - f_x * v ) / ( e_x + g_x * v )
2032 
2033  if ( u < -eps_bilinear .or. u > 1.0_rp+eps_bilinear .or. v < -eps_bilinear .or. v > 1.0_rp+eps_bilinear ) then
2034  log_error("INTERP_bilinear_inv",*) 'Unexpected error occured', u, v
2035  log_error_cont(*) x_ref0, x_ref1, x_ref2, x_ref3
2036  log_error_cont(*) y_ref0, y_ref1, y_ref2, y_ref3
2037  log_error_cont(*) x, y
2038  log_error_cont(*) k0, k1, k2, sig
2039  error = .true.
2040  return
2041  end if
2042  u = min( max( u, 0.0_rp ), 1.0_rp )
2043  v = min( max( v, 0.0_rp ), 1.0_rp )
2044 
2045  error = .false.
2046 
2047  return

References scale_const::const_eps, and cross().

Referenced by interp_factor2d_linear_xy().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ cross()

real(rp) function scale_interp::cross ( real(rp), intent(in)  x0,
real(rp), intent(in)  y0,
real(rp), intent(in)  x1,
real(rp), intent(in)  y1 
)

cross product

Definition at line 2115 of file scale_interp.F90.

2115  implicit none
2116  real(RP), intent(in) :: x0, y0
2117  real(RP), intent(in) :: x1, y1
2118  real(RP) :: cross
2119 
2120  cross = x0 * y1 - x1 * y0
2121 

References scale_const::const_eps, and scale_const::const_undef.

Referenced by interp_bilinear_inv().

Here is the caller graph for this function:

◆ spline_exec()

subroutine scale_interp::spline_exec ( integer, intent(in)  KA_ref,
integer, intent(in)  kmax,
integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
integer, dimension(ka,2), intent(in)  idx_k,
real(rp), dimension(ka), intent(in)  vfact,
real(rp), dimension(ka_ref), intent(in)  hgt_ref,
real(rp), dimension (ka), intent(in)  hgt,
real(rp), dimension(ka_ref), intent(in)  val_ref,
integer, dimension (ka_ref), intent(in)  idx,
integer, dimension (ka_ref), intent(in)  idx_r,
real(rp), dimension (ka_ref), intent(in)  U,
real(rp), dimension (ka_ref), intent(in)  FDZ,
real(rp), dimension(ka), intent(out)  val 
)

Definition at line 2258 of file scale_interp.F90.

2258  use scale_const, only: &
2259  undef => const_undef, &
2260  eps => const_eps
2261  implicit none
2262  integer, intent(in) :: KA_ref, kmax
2263  integer, intent(in) :: KA, KS, KE
2264  integer, intent(in) :: idx_k(KA,2)
2265  real(RP), intent(in) :: vfact(KA)
2266  real(RP), intent(in) :: hgt_ref(KA_ref)
2267  real(RP), intent(in) :: hgt (KA)
2268  real(RP), intent(in) :: val_ref(KA_ref)
2269  integer, intent(in) :: idx (KA_ref)
2270  integer, intent(in) :: idx_r (KA_ref)
2271  real(RP), intent(in) :: U (KA_ref)
2272  real(RP), intent(in) :: FDZ (KA_ref)
2273 
2274  real(RP), intent(out) :: val(KA)
2275 
2276  real(RP) :: c1, c2, c3, d
2277  integer :: k, kk, kk2
2278 
2279  do k = ks, ke
2280  kk = idx_k(k,1)
2281  if ( kk == -1 ) then
2282  val(k) = undef
2283  else if ( idx_k(k,2) == -1 ) then
2284  val(k) = val_ref(kk)
2285  else if ( kk < idx(1) .or. kk >= idx(kmax) ) then ! linear interpolation
2286  if ( abs( val_ref(kk) - undef ) < eps .or. abs( val_ref(kk+1) - undef ) < eps ) then
2287  val(k) = undef
2288  else
2289  val(k) = val_ref(kk) * vfact(k) + val_ref(kk+1) * ( 1.0_rp - vfact(k) )
2290  end if
2291  else
2292  kk2 = idx_r(kk)
2293  kk = idx(kk2)
2294  c3 = ( u(kk2+1) - u(kk2) ) / fdz(kk2+1)
2295  c2 = 3.0_rp * u(kk2)
2296  c1 = ( val_ref(idx(kk2+1)) - val_ref(kk) ) / fdz(kk2+1) - ( 2.0_rp * u(kk2) + u(kk2+1) ) * fdz(kk2+1)
2297  d = hgt(k) - hgt_ref(kk)
2298 
2299  val(k) = ( ( c3 * d + c2 ) * d + c1 ) * d + val_ref(kk)
2300  end if
2301  end do
2302 
2303  return

References scale_const::const_eps, and scale_const::const_undef.

Referenced by interp_interp1d(), and interp_interp3d().

Here is the caller graph for this function:
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:342
scale_sort
module SORT
Definition: scale_sort.F90:11
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_const
module CONSTANT
Definition: scale_const.F90:11
Generated by   doxygen 1.8.17