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, missing)
 
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, missing)
 
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 90 of file scale_interp.F90.

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

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

256  !$acc routine vector
257  use scale_prc, only: &
258  prc_abort
259  use scale_const, only: &
260  undef => const_undef, &
261  eps => const_eps
262  implicit none
263  integer, intent(in) :: KA_ref, KS_ref, KE_ref ! number of z-direction (reference)
264  integer, intent(in) :: KA, KS, KE ! number of z-direction (target)
265 
266  real(RP), intent(in) :: hgt_ref(KA_ref) ! height [m] (reference)
267  real(RP), intent(in) :: hgt (KA) ! height [m] (target)
268 
269  integer, intent(out) :: idx_k (KA,2) ! k-index in reference (target)
270  real(RP), intent(out) :: vfact (KA) ! horizontal interp factor (target)
271 
272  logical, intent(in), optional :: flag_extrap ! when true, extrapolation will be executed (just copy)
273 
274  integer :: idx(KA_ref), kmax
275  logical :: flag_extrap_
276 
277  integer :: k, kk
278  !---------------------------------------------------------------------------
279 
280  if ( present(flag_extrap) ) then
281  flag_extrap_ = flag_extrap
282  else
283  flag_extrap_ = .true.
284  end if
285 
286  ! search valid levels
287  kmax = 0
288  do k = ks_ref, ke_ref
289  if ( hgt_ref(k) > undef ) then
290  kmax = kmax + 1
291  idx(kmax) = k
292  end if
293  end do
294 
295  do k = ks, ke
296  idx_k(k,1) = -1
297  idx_k(k,2) = -1
298 
299  if ( hgt(k) < hgt_ref(idx(1)) - eps ) then
300  if ( flag_extrap_ ) then
301  idx_k(k,1) = idx(1)
302  idx_k(k,2) = -1
303  vfact(k) = 1.0_rp
304  else
305  idx_k(k,1) = -1
306  idx_k(k,2) = -1
307  vfact(k) = undef
308  end if
309  elseif( hgt(k) < hgt_ref(idx(1)) ) then
310  idx_k(k,1) = idx(1)
311  idx_k(k,2) = -1
312  vfact(k) = 1.0_rp
313  elseif( hgt(k) > hgt_ref(idx(kmax)) + eps ) then
314  if ( flag_extrap_ ) then
315  idx_k(k,1) = idx(kmax)
316  idx_k(k,2) = -1
317  vfact(k) = 1.0_rp
318  else
319  idx_k(k,1) = -1
320  idx_k(k,2) = -1
321  vfact(k) = undef
322  end if
323  elseif( hgt(k) >= hgt_ref(idx(kmax)) ) then
324  idx_k(k,1) = idx(kmax)
325  idx_k(k,2) = -1
326  vfact(k) = 1.0_rp
327  else
328  !$acc loop seq
329  do kk = 1, kmax-1
330  if ( hgt(k) >= hgt_ref(idx(kk) ) &
331  .AND. hgt(k) < hgt_ref(idx(kk+1)) ) then
332  idx_k(k,1) = idx(kk)
333  idx_k(k,2) = idx(kk+1)
334  vfact(k) = ( hgt_ref(idx(kk+1)) - hgt(k) ) &
335  / ( hgt_ref(idx(kk+1)) - hgt_ref(idx(kk)) )
336 
337  exit
338  endif
339  enddo
340  endif
341 
342  enddo ! k-loop
343 
344  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 357 of file scale_interp.F90.

357  use scale_sort, only: &
358  sort_exec
359  implicit none
360  integer, intent(in) :: IA_ref, JA_ref
361  integer, intent(in) :: IA, JA
362 
363  real(RP), intent(in) :: lon_ref(IA_ref)
364  real(RP), intent(in) :: lat_ref(JA_ref)
365  real(RP), intent(in) :: lon(IA,JA)
366  real(RP), intent(in) :: lat(IA,JA)
367  integer, intent(out) :: idx_i(IA,JA,4)
368  integer, intent(out) :: idx_j(IA,JA,4)
369  real(RP), intent(out) :: hfact(IA,JA,4)
370 
371  real(RP) :: workh(4)
372  integer :: worki(4), workj(4)
373 
374  real(RP) :: f1, f2
375  integer :: i, j, ii, jj
376 
377  call prof_rapstart('INTERP_fact',3)
378 
379  !$omp parallel do schedule(dynamic) collapse(2) &
380  !$omp private(f1,f2,workh,worki,workj)
381  !$acc kernels
382  do j = 1, ja
383  !$acc loop private(workh, worki, workj)
384  do i = 1, ia
385 
386  ! longitude
387  if ( lon(i,j) < lon_ref(1) ) then
388  idx_i(i,j,:) = 1
389  hfact(i,j,:) = 0.0_rp
390  else if ( lon(i,j) > lon_ref(ia_ref) ) then
391  idx_i(i,j,:) = ia_ref
392  hfact(i,j,:) = 0.0_rp
393  else if ( lon(i,j) == lon_ref(ia_ref) ) then
394  idx_i(i,j,:) = ia_ref
395  hfact(i,j,1) = 0.0_rp
396  hfact(i,j,2) = 1.0_rp
397  hfact(i,j,3) = 0.0_rp
398  hfact(i,j,4) = 1.0_rp
399  else
400  do ii = 1, ia_ref
401  if ( lon(i,j) < lon_ref(ii) ) then
402  f1 = ( lon_ref(ii) - lon(i,j) ) / ( lon_ref(ii) - lon_ref(ii-1) )
403  f2 = ( lon(i,j) - lon_ref(ii-1) ) / ( lon_ref(ii) - lon_ref(ii-1) )
404  hfact(i,j,1) = f1
405  hfact(i,j,2) = f2
406  hfact(i,j,3) = f1
407  hfact(i,j,4) = f2
408  idx_i(i,j,1) = ii-1
409  idx_i(i,j,2) = ii
410  idx_i(i,j,3) = ii-1
411  idx_i(i,j,4) = ii
412  exit
413  end if
414  end do
415  end if
416 
417  ! latitude
418  if ( lat(i,j) < lat_ref(1) ) then
419  idx_j(i,j,:) = 1
420  hfact(i,j,:) = 0.0_rp
421  else if ( lat(i,j) > lat_ref(ja_ref) ) then
422  idx_j(i,j,:) = ja_ref
423  hfact(i,j,:) = 0.0_rp
424  else if ( lat(i,j) == lat_ref(ja_ref) ) then
425  idx_j(i,j,:) = ja_ref
426  hfact(i,j,1) = 0.0_rp
427  hfact(i,j,2) = 0.0_rp
428  !hfact(i,j,3) = hfact(i,j,3)
429  !hfact(i,j,4) = hfact(i,j,4)
430  else
431  do jj = 1, ja_ref
432  if ( lat(i,j) < lat_ref(jj) ) then
433  f1 = ( lat_ref(jj) - lat(i,j) ) / ( lat_ref(jj) - lat_ref(jj-1) )
434  f2 = ( lat(i,j) - lat_ref(jj-1) ) / ( lat_ref(jj) - lat_ref(jj-1) )
435  hfact(i,j,1) = hfact(i,j,1) * f1
436  hfact(i,j,2) = hfact(i,j,2) * f1
437  hfact(i,j,3) = hfact(i,j,3) * f2
438  hfact(i,j,4) = hfact(i,j,4) * f2
439  idx_j(i,j,1) = jj-1
440  idx_j(i,j,2) = jj-1
441  idx_j(i,j,3) = jj
442  idx_j(i,j,4) = jj
443  exit
444  end if
445  end do
446  end if
447 
448  do ii = 1, 4
449  workh(ii) = hfact(i,j,ii)
450  worki(ii) = idx_i(i,j,ii)
451  workj(ii) = idx_j(i,j,ii)
452  end do
453  call sort_exec( 4, & ! [IN]
454  workh(:), worki(:), workj(:), & ! [INOUT]
455  reverse = .true. ) ! [IN]
456  do ii = 1, 4
457  hfact(i,j,ii) = workh(ii)
458  idx_i(i,j,ii) = worki(ii)
459  idx_j(i,j,ii) = workj(ii)
460  end do
461 
462  end do
463  end do
464  !$acc end kernels
465 
466  call prof_rapend ('INTERP_fact',3)
467 
468  return

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

Referenced by mod_cnv2d::cnv2d_exec(), 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,
logical, intent(in), optional  missing 
)

Definition at line 483 of file scale_interp.F90.

483  use scale_prc, only: &
484  prc_abort
485  use scale_const, only: &
486  undef => const_undef, &
487  eps => const_eps
488  use scale_sort, only: &
489  sort_exec
490  implicit none
491  integer, intent(in) :: IA_ref, JA_ref
492  integer, intent(in) :: IA, JA
493 
494  real(RP), intent(in) :: x_ref(IA_ref,JA_ref)
495  real(RP), intent(in) :: y_ref(IA_ref,JA_ref)
496  real(RP), intent(in) :: x(IA)
497  real(RP), intent(in) :: y(JA)
498 
499  integer, intent(out) :: idx_i(IA,JA,4)
500  integer, intent(out) :: idx_j(IA,JA,4)
501  real(RP), intent(out) :: hfact(IA,JA,4)
502 
503  logical, intent(in), optional :: zonal
504  logical, intent(in), optional :: pole
505  logical, intent(in), optional :: missing
506 
507  real(RP) :: u, v
508  integer :: inc_i, inc_j
509  logical :: error, err
510  logical :: zonal_, pole_
511  logical :: missing_
512 
513  real(RP) :: workh(4)
514  integer :: worki(4), workj(4)
515 
516  integer :: ii0, jj0
517  integer :: ii, jj
518  integer :: i1, i2, i3, i4
519  integer :: j1, j2, j3, j4
520  integer :: i, j, l
521  integer :: ite, ite_max
522 
523  call prof_rapstart('INTERP_fact',3)
524 
525  if ( present(zonal) ) then
526  zonal_ = zonal
527  else
528  zonal_ = .false.
529  end if
530  if ( present(pole) ) then
531  pole_ = pole
532  else
533  pole_ = .false.
534  end if
535  if ( present(missing) ) then
536  missing_ = missing
537  else
538  missing_ = .false.
539  end if
540 
541  ii0 = ( ia_ref + 1 ) / 2
542  jj0 = ( ja_ref + 1 ) / 2
543  error = .false.
544 
545  ite_max = ia_ref + ja_ref
546 
547 
548  !$omp parallel do &
549  !$omp private(inc_i,inc_j,ii,jj,i1,i2,i3,i4,j1,j2,j3,j4,u,v,err,workh,worki,workj) &
550  !$omp firstprivate(ii0,jj0)
551  !$acc kernels
552  !$acc loop reduction(.or.:error)
553  do j = 1, ja
554  !$acc loop private(inc_i,inc_j,ii,jj,i1,i2,i3,i4,j1,j2,j3,j4,u,v,err,workh,worki,workj) reduction(.or.:error)
555  do i = 1, ia
556 
557 #ifndef _OPENACC
558  if ( i==1 ) then
559 #endif
560  ii = ii0
561  jj = jj0
562 #ifndef _OPENACC
563  end if
564 #endif
565  !$acc loop seq
566  do ite = 1, ite_max
567  i1 = ii
568  i2 = ii + 1
569  i3 = ii + 1
570  i4 = ii
571  j1 = jj
572  j2 = jj
573  j3 = jj + 1
574  j4 = jj + 1
575  if ( ii == ia_ref ) then ! zonal_ must be .true.
576  i2 = 1
577  i3 = 1
578  end if
579  if ( jj == 0 ) then ! pole_ must be .true.
580  j1 = 1
581  j2 = 1
582  i2 = ia_ref / 2 + ii
583  if ( i2 > ia_ref ) i2 = i2 - ia_ref
584  i1 = i2 + 1
585  if ( i1 > ia_ref ) i1 = i1 - ia_ref
586  end if
587  if ( jj == ja_ref ) then ! pole_ must be .true.
588  j3 = ja_ref
589  j4 = ja_ref
590  i3 = ia_ref / 2 + ii
591  if ( i3 > ia_ref ) i3 = i3 - ia_ref
592  i4 = i3 + 1
593  if ( i4 > ia_ref ) i4 = i4 - ia_ref
594  end if
595  if ( abs( x_ref(i1,j1) - undef ) < eps .or. &
596  abs( x_ref(i2,j2) - undef ) < eps .or. &
597  abs( x_ref(i3,j3) - undef ) < eps .or. &
598  abs( x_ref(i4,j4) - undef ) < eps ) then
599  if ( ii == ia_ref-1 ) then
600  ii = 1
601  if ( jj == ja_ref-1 ) then
602  jj = 1
603  else
604  jj = jj + 2
605  end if
606  else
607  ii = ii + 2
608  end if
609  if ( ii == ia_ref ) then
610  ii = 1
611  jj = jj + 2
612  end if
613  if ( jj == ja_ref ) then
614  jj = 1
615  ii = ii + 2
616  if ( ii == ia_ref ) ii = 1
617  end if
618  else
619  call interp_check_inside( &
620  x_ref(i1,j1), x_ref(i2,j2), x_ref(i3,j3), x_ref(i4,j4), & ! (in)
621  y_ref(i1,j1), y_ref(i2,j2), y_ref(i3,j3), y_ref(i4,j4), & ! (in)
622  x(i), y(j), & ! (in)
623  inc_i, inc_j, & ! (out)
624  err ) ! (out)
625  if ( err ) error = .true.
626  if ( error ) exit
627  if ( inc_i == 0 .and. inc_j == 0 ) then ! inside the quadrilateral
628  call interp_bilinear_inv( &
629  x_ref(i1,j1), x_ref(i2,j2), x_ref(i3,j3), x_ref(i4,j4), & ! (in)
630  y_ref(i1,j1), y_ref(i2,j2), y_ref(i3,j3), y_ref(i4,j4), & ! (in)
631  x(i), y(j), & ! (in)
632  u, v, & ! (out)
633  err ) ! (out)
634  if ( err ) error = .true.
635  if ( error ) exit
636  idx_i(i,j,1) = i1
637  idx_i(i,j,2) = i2
638  idx_i(i,j,3) = i3
639  idx_i(i,j,4) = i4
640  idx_j(i,j,1) = j1
641  idx_j(i,j,2) = j2
642  idx_j(i,j,3) = j3
643  idx_j(i,j,4) = j4
644  hfact(i,j,1) = ( 1.0_rp - u ) * ( 1.0_rp - v )
645  hfact(i,j,2) = ( u ) * ( 1.0_rp - v )
646  hfact(i,j,3) = ( u ) * ( v )
647  hfact(i,j,4) = ( 1.0_rp - u ) * ( v )
648  exit
649  end if
650  ii = ii + inc_i
651  jj = jj + inc_j
652  if ( zonal_ .or. pole_ ) then
653  if ( ii == 0 ) ii = ia_ref
654  if ( ii == ia_ref+1 ) ii = 1
655  if ( pole_ ) then
656  jj = max( jj, 0 )
657  jj = min( jj, ja_ref )
658  else
659  jj = max( jj, 1 )
660  jj = min( jj, ja_ref-1 )
661  end if
662  else
663  if ( ii == 0 .and. jj == 0 ) then
664  ii = ia_ref - 1
665  jj = ja_ref - 1
666  end if
667  if ( ii == 0 .and. jj == ja_ref ) then
668  ii = ia_ref - 1
669  jj = 1
670  end if
671  if ( ii == ia_ref .and. jj == 0 ) then
672  ii = 1
673  jj = ja_ref - 1
674  end if
675  if ( ii == ia_ref .and. jj == ja_ref ) then
676  ii = 1
677  jj = 1
678  end if
679  ii = max( min( ii, ia_ref-1 ), 1 )
680  jj = max( min( jj, ja_ref-1 ), 1 )
681  end if
682  end if
683  end do
684  if ( ite == ite_max+1 ) then
685  idx_i(i,j,:) = 1
686  idx_j(i,j,:) = 1
687  hfact(i,j,:) = 0.0_rp
688  if ( .not. missing_ ) then
689  error = .true.
690 #ifndef _OPENACC
691  log_error("INTERP_factor2d_linear_xy",*) 'iteration max has been reached'
692  log_error_cont(*) 'The point may be out of region.', i, j, x(i), y(j)
693  log_error_cont(*) minval(x_ref), maxval(x_ref), minval(y_ref), maxval(y_ref)
694  log_error_cont(*) x_ref(1,1), x_ref(ia_ref,1),x_ref(1,ja_ref),x_ref(ia_ref,ja_ref)
695  log_error_cont(*) y_ref(1,1), y_ref(ia_ref,1),y_ref(1,ja_ref),y_ref(ia_ref,ja_ref)
696  call prc_abort
697 #endif
698  end if
699  end if
700 
701 #ifndef _OPENACC
702  if ( i==1 ) then
703  ii0 = ii
704  jj0 = jj
705  end if
706 
707  if ( error ) exit
708 #endif
709 
710  do l = 1, 4
711  workh(l) = hfact(i,j,l)
712  worki(l) = idx_i(i,j,l)
713  workj(l) = idx_j(i,j,l)
714  end do
715  call sort_exec( 4, & ! [IN]
716  workh(:), worki(:), workj(:), & ! [INOUT]
717  reverse = .true. ) ! [IN]
718  do l = 1, 4
719  hfact(i,j,l) = workh(l)
720  idx_i(i,j,l) = worki(l)
721  idx_j(i,j,l) = workj(l)
722  end do
723 
724  end do
725  end do
726  !$acc end kernels
727 
728  if ( error ) call prc_abort
729 
730  call prof_rapend ('INTERP_fact',3)
731 
732  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(), 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 748 of file scale_interp.F90.

748  use scale_prc, only: &
749  prc_abort
750  implicit none
751  integer, intent(in) :: npoints ! number of interpolation point for horizontal
752  integer, intent(in) :: IA_ref ! number of x-direction (reference)
753  integer, intent(in) :: JA_ref ! number of y-direction (reference)
754  integer, intent(in) :: IA ! number of x-direction (target)
755  integer, intent(in) :: JA ! number of y-direction (target)
756  real(RP), intent(in) :: lon_ref(IA_ref,JA_ref) ! longitude [rad] (reference)
757  real(RP), intent(in) :: lat_ref(IA_ref,JA_ref) ! latitude [rad] (reference)
758  real(RP), intent(in) :: lon (IA,JA) ! longitude [rad] (target)
759  real(RP), intent(in) :: lat (IA,JA) ! latitude [rad] (target)
760  integer, intent(out) :: idx_i(IA,JA,npoints) ! i-index in reference (target)
761  integer, intent(out) :: idx_j(IA,JA,npoints) ! j-index in reference (target)
762  real(RP), intent(out) :: hfact(IA,JA,npoints) ! horizontal interp factor (target)
763 
764  real(RP), intent(in), optional :: search_limit
765  logical, intent(in), optional :: latlon_structure
766  real(RP), intent(in), optional :: lon_1d(IA_ref), lat_1d(JA_ref)
767  integer, intent(in), optional :: weight_order
768 
769  logical :: ll_struct_
770 
771  real(RP) :: lon_min, lon_max
772  real(RP) :: lat_min, lat_max
773  real(RP) :: dlon, dlat
774 
775  ! for structure grid
776  integer :: psizex, psizey
777  real(RP) :: lon0, lat0
778  integer, allocatable :: i0(:), i1(:), j0(:), j1(:)
779 
780  ! for unstructure grid
781  integer :: nsize, psize, nidx_max
782  integer, allocatable :: idx_blk(:,:,:), nidx(:,:)
783  integer :: idx_ref(npoints)
784 
785  integer :: idx_it(npoints)
786  integer :: idx_jt(npoints)
787  real(RP) :: hfactt(npoints)
788 
789  integer :: i, j, ii, jj, n
790  !---------------------------------------------------------------------------
791 
792  call prof_rapstart('INTERP_fact',3)
793 
794  if ( present(latlon_structure) ) then
795  ll_struct_ = latlon_structure
796  else
797  ll_struct_ = .false.
798  end if
799 
800  hfact(:,:,:) = 0.0_rp
801 
802 
803  if ( ll_struct_ ) then
804 
805  if ( ia_ref > 10 ) then
806  psizex = int( 2.0_rp*sqrt(real(ia_ref+1,rp)) )
807  else
808  psizex = 1
809  end if
810  if ( ja_ref > 10 ) then
811  psizey = int( 2.0_rp*sqrt(real(ja_ref+1,rp)) )
812  else
813  psizey = 1
814  end if
815 
816  allocate( i0(psizex), i1(psizex) )
817  allocate( j0(psizey), j1(psizey) )
818 
819  dlon = ( lon_max - lon_min ) / psizex
820  dlat = ( lat_max - lat_min ) / psizey
821 
822  do ii = 1, psizex
823  lon0 = lon_min + dlon * (ii-1)
824  do i = 1, ia_ref
825  if ( lon_1d(i) >= lon0 ) then
826  i0(ii) = i
827  exit
828  end if
829  end do
830  end do
831  do ii = 1, psizex-1
832  i1(ii) = i0(ii+1) - 1
833  end do
834  i1(psizex) = ia_ref
835 
836  do jj = 1, psizey
837  lat0 = lat_min + dlat * (jj-1)
838  do j = 1, ja_ref
839  if ( lat_1d(j) >= lat0 ) then
840  j0(jj) = j
841  exit
842  end if
843  end do
844  end do
845  do jj = 1, psizey-1
846  j1(jj) = j0(jj+1) - 1
847  end do
848  j1(psizey) = ja_ref
849 
850  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
851  !$omp private(idx_it,idx_jt,hfactt)
852  !$acc kernels copyin(lon_1d, lat_1d, i0, i1, j0, j1, lon, lat) copyout(idx_i, idx_j, hfact)
853  !$acc loop independent
854  do j = 1, ja
855  !$acc loop private(idx_it,idx_jt,hfactt) independent
856  do i = 1, ia
857  ! main search
858  call interp_search_horiz_struct( npoints, & ! [IN]
859  psizex, psizey, & ! [IN]
860  ia_ref, ja_ref, & ! [IN]
861  lon_1d(:), lat_1d(:), & ! [IN]
862  lon_min, lat_min, & ! [IN]
863  dlon, dlat, & ! [IN]
864  i0(:), i1(:), j0(:), j1(:), & ! [IN]
865  lon(i,j), lat(i,j), & ! [IN]
866  idx_it(:), idx_jt(:), & ! [OUT]
867  hfactt(:), & ! [OUT]
868  search_limit = search_limit, & ! [IN]
869  weight_order = weight_order ) ! [IN]
870  do n = 1, npoints
871  idx_i(i,j,n) = idx_it(n)
872  idx_j(i,j,n) = idx_jt(n)
873  hfact(i,j,n) = hfactt(n)
874  end do
875  enddo
876  enddo
877  !$acc end kernels
878 
879  deallocate( i0, i1, j0, j1 )
880 
881  else
882 
883  nsize = ia_ref * ja_ref
884  if ( nsize > 100 ) then
885  psize = int( sqrt(2.0_rp*sqrt(real(nsize,rp))) )
886  nidx_max = nsize / psize * interp_buffer_size_fact
887  else
888  psize = 1
889  nidx_max = nsize
890  end if
891 
892  allocate(idx_blk(nidx_max,psize,psize))
893  allocate(nidx( psize,psize))
894 
895  !$acc data copyin(lon_ref, lat_ref, lon, lat) copyout(idx_i, idx_j, hfact) create(idx_blk, nidx)
896 
897  call interp_div_block(nsize, psize, nidx_max, & ! [IN]
898  lon_ref(:,:), lat_ref(:,:), & ! [IN]
899  idx_blk(:,:,:), nidx(:,:), & ! [OUT]
900  lon_min, lon_max, & ! [OUT]
901  lat_min, lat_max, & ! [OUT]
902  dlon, dlat ) ! [OUT]
903 
904  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
905  !$omp private(idx_ref,hfactt)
906  !$acc kernels
907  !$acc loop independent
908  do j = 1, ja
909  !$acc loop private(idx_ref,hfactt) independent
910  do i = 1, ia
911  ! main search
912  call interp_search_horiz( npoints, & ! [IN]
913  nsize, & ! [IN]
914  lon_ref(:,:), lat_ref(:,:), & ! [IN]
915  lon_min, lon_max, & ! [IN]
916  lat_min, lat_max, & ! [IN]
917  psize, nidx_max, & ! [IN]
918  dlon, dlat, & ! [IN]
919  idx_blk(:,:,:), nidx(:,:), & ! [IN]
920  lon(i,j), lat(i,j), & ! [IN]
921  idx_ref(:), & ! [OUT]
922  hfactt(:), & ! [OUT]
923  search_limit = search_limit, & ! [IN]
924  weight_order = weight_order ) ! [IN]
925  do n = 1, npoints
926  idx_i(i,j,n) = mod(idx_ref(n) - 1, ia_ref) + 1
927  idx_j(i,j,n) = ( idx_ref(n) - 1 ) / ia_ref + 1
928  hfact(i,j,n) = hfactt(n)
929  end do
930  enddo
931  enddo
932  !$acc end kernels
933 
934  !$acc end data
935 
936  deallocate(idx_blk, nidx)
937 
938  end if
939 
940  call prof_rapend ('INTERP_fact',3)
941 
942  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().

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 962 of file scale_interp.F90.

962  implicit none
963 
964  integer, intent(in) :: KA_ref, KS_ref, KE_ref ! number of z-direction (reference)
965  integer, intent(in) :: IA_ref ! number of x-direction (reference)
966  integer, intent(in) :: JA_ref ! number of y-direction (reference)
967  integer, intent(in) :: KA, KS, KE ! number of z-direction (target)
968  integer, intent(in) :: IA ! number of x-direction (target)
969  integer, intent(in) :: JA ! number of y-direction (target)
970  real(RP), intent(in) :: lon_ref(IA_ref) ! longitude (reference)
971  real(RP), intent(in) :: lat_ref(JA_ref) ! latitude (reference)
972  real(RP), intent(in) :: hgt_ref(KA_ref,IA_ref,JA_ref) ! height [m] (reference)
973  real(RP), intent(in) :: lon(IA,JA) ! longitude (target)
974  real(RP), intent(in) :: lat(IA,JA) ! latitude (target)
975  real(RP), intent(in) :: hgt(KA,IA,JA) ! longitude [m] (target)
976 
977  integer, intent(out) :: idx_i(IA,JA,4) ! i-index in reference (target)
978  integer, intent(out) :: idx_j(IA,JA,4) ! j-index in reference (target)
979  real(RP), intent(out) :: hfact(IA,JA,4) ! horizontal interp factor (target)
980  integer, intent(out) :: idx_k(KA,2,IA,JA,4) ! k-index in reference (target)
981  real(RP), intent(out) :: vfact(KA, IA,JA,4) ! vertical interp factor (target)
982 
983  logical, intent(in), optional :: flag_extrap ! when true, vertical extrapolation will be executed (just copy)
984 
985  integer :: i, j, ii, jj, n
986 
987  !$acc data copyin(lon_ref, lat_ref, hgt_ref, lon, lat, hgt) copyout(idx_i, idx_j, hfact, idx_k, vfact)
988 
989  call interp_factor2d_linear_latlon( &
990  ia_ref, ja_ref, & ! [IN]
991  ia, ja, & ! [IN]
992  lon_ref(:), lat_ref(:), & ! [IN]
993  lon(:,:), lat(:,:), & ! [IN]
994  idx_i(:,:,:), idx_j(:,:,:), hfact(:,:,:) ) ! [OUT]
995 
996  call prof_rapstart('INTERP_fact',3)
997 
998  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
999  !$omp private(ii,jj)
1000  !$acc kernels
1001  !$acc loop independent
1002  do j = 1, ja
1003  !$acc loop independent
1004  do i = 1, ia
1005  !$acc loop seq
1006  do n = 1, 4
1007  ii = idx_i(i,j,n)
1008  jj = idx_j(i,j,n)
1009  call interp_factor1d( ka_ref, ks_ref, ke_ref, & ! [IN]
1010  ka, ks, ke, & ! [IN]
1011  hgt_ref(:,ii,jj), & ! [IN]
1012  hgt(:,i,j), & ! [IN]
1013  idx_k(:,:,i,j,n), & ! [OUT]
1014  vfact(:, i,j,n), & ! [OUT]
1015  flag_extrap = flag_extrap ) ! [IN, optional]
1016 
1017  enddo
1018  enddo
1019  enddo
1020  !$acc end kernels
1021 
1022  !$acc end data
1023 
1024  call prof_rapend('INTERP_fact',3)
1025 
1026  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,
logical, intent(in), optional  missing 
)

Definition at line 1048 of file scale_interp.F90.

1048  implicit none
1049  integer, intent(in) :: KA_ref, KS_ref, KE_ref ! number of z-direction (reference)
1050  integer, intent(in) :: IA_ref ! number of x-direction (reference)
1051  integer, intent(in) :: JA_ref ! number of y-direction (reference)
1052  integer, intent(in) :: KA, KS, KE ! number of z-direction (target)
1053  integer, intent(in) :: IA ! number of x-direction (target)
1054  integer, intent(in) :: JA ! number of y-direction (target)
1055  real(RP), intent(in) :: x_ref(IA_ref,JA_ref) ! x point (reference)
1056  real(RP), intent(in) :: y_ref(IA_ref,JA_ref) ! y point (reference)
1057  real(RP), intent(in) :: hgt_ref(KA_ref,IA_ref,JA_ref) ! height [m] (reference)
1058  real(RP), intent(in) :: x (IA) ! x point (target)
1059  real(RP), intent(in) :: y (JA) ! y point (target)
1060  real(RP), intent(in) :: hgt(KA,IA,JA) ! longitude [m] (target)
1061 
1062  integer, intent(out) :: idx_i(IA,JA,4) ! i-index in reference
1063  integer, intent(out) :: idx_j(IA,JA,4) ! j-index in reference
1064  real(RP), intent(out) :: hfact(IA,JA,4) ! horizontal interp factor
1065  integer, intent(out) :: idx_k(KA,2,IA,JA,4) ! k-index in reference
1066  real(RP), intent(out) :: vfact(KA, IA,JA,4) ! vertical interp factor
1067 
1068  logical, intent(in), optional :: flag_extrap ! when true, vertical extrapolation will be executed (just copy)
1069  logical, intent(in), optional :: zonal
1070  logical, intent(in), optional :: pole
1071  logical, intent(in), optional :: missing
1072 
1073  integer :: i, j, ii, jj, n
1074 
1075  !$acc data copyin(x_ref, y_ref, hgt_ref, x, y, hgt) copyout(idx_i, idx_j, hfact, idx_k, vfact)
1076 
1077  call interp_factor2d_linear_xy( &
1078  ia_ref, ja_ref, & ! [IN]
1079  ia, ja, & ! [IN]
1080  x_ref(:,:), y_ref(:,:), & ! [IN]
1081  x(:), y(:), & ! [IN]
1082  idx_i(:,:,:), idx_j(:,:,:), hfact(:,:,:), & ! [OUT]
1083  zonal = zonal, pole = pole, & ! [IN]
1084  missing = missing ) ! [IN]
1085 
1086  call prof_rapstart('INTERP_fact',3)
1087 
1088  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
1089  !$omp private(ii,jj)
1090  !$acc kernels
1091  !$acc loop independent
1092  do j = 1, ja
1093  !$acc loop independent
1094  do i = 1, ia
1095  !$acc loop seq
1096  do n = 1, 4
1097  ii = idx_i(i,j,n)
1098  jj = idx_j(i,j,n)
1099  call interp_factor1d( ka_ref, ks_ref, ke_ref, & ! [IN]
1100  ka, ks, ke, & ! [IN]
1101  hgt_ref(:,ii,jj), & ! [IN]
1102  hgt(:,i,j), & ! [IN]
1103  idx_k(:,:,i,j,n), & ! [OUT]
1104  vfact(: ,i,j,n), & ! [OUT]
1105  flag_extrap = flag_extrap ) ! [IN, optional]
1106 
1107  enddo
1108  enddo
1109  enddo
1110  !$acc end kernels
1111 
1112  !$acc end data
1113 
1114  call prof_rapend('INTERP_fact',3)
1115 
1116  return

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

Referenced by mod_cnvuser::cnvuser().

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 1136 of file scale_interp.F90.

1136  use scale_prc, only: &
1137  prc_abort
1138  implicit none
1139  integer, intent(in) :: npoints ! number of interpolation point for horizontal
1140  integer, intent(in) :: KA_ref, KS_ref, KE_ref ! number of z-direction (reference)
1141  integer, intent(in) :: IA_ref ! number of x-direction (reference)
1142  integer, intent(in) :: JA_ref ! number of y-direction (reference)
1143  integer, intent(in) :: KA, KS, KE ! number of z-direction (target)
1144  integer, intent(in) :: IA ! number of x-direction (target)
1145  integer, intent(in) :: JA ! number of y-direction (target)
1146 
1147  real(RP), intent(in) :: lon_ref(IA_ref,JA_ref) ! longitude [rad] (reference)
1148  real(RP), intent(in) :: lat_ref(IA_ref,JA_ref) ! latitude [rad] (reference)
1149  real(RP), intent(in) :: hgt_ref(KA_ref,IA_ref,JA_ref) ! height [m] (reference)
1150  real(RP), intent(in) :: lon (IA,JA) ! longitude [rad] (target)
1151  real(RP), intent(in) :: lat (IA,JA) ! latitude [rad] (target)
1152  real(RP), intent(in) :: hgt (KA,IA,JA) ! longitude [m] (target)
1153 
1154  integer, intent(out) :: idx_i(IA,JA,npoints) ! i-index in reference
1155  integer, intent(out) :: idx_j(IA,JA,npoints) ! j-index in reference
1156  real(RP), intent(out) :: hfact(IA,JA,npoints) ! horizontal interp factor
1157  integer, intent(out) :: idx_k(KA,2,IA,JA,npoints) ! k-index in reference
1158  real(RP), intent(out) :: vfact(KA, IA,JA,npoints) ! vertical interp factor
1159 
1160  logical, intent(in), optional :: flag_extrap ! when true, vertical extrapolation will be executed (just copy)
1161 
1162  integer :: nsize, psize, nidx_max
1163  integer, allocatable :: idx_blk(:,:,:), nidx(:,:)
1164  real(RP) :: lon_min, lon_max
1165  real(RP) :: lat_min, lat_max
1166  real(RP) :: dlon, dlat
1167  integer :: idx_ref(npoints)
1168  real(RP) :: hfactt(npoints)
1169 
1170  integer :: i, j, ii, jj, n
1171  !---------------------------------------------------------------------------
1172 
1173  call prof_rapstart('INTERP_fact',3)
1174 
1175  nsize = ia_ref * ja_ref
1176  if ( nsize > 100 ) then
1177  psize = int( sqrt(2.0_rp*sqrt(real(nsize,rp))) )
1178  nidx_max = nsize / psize * interp_buffer_size_fact
1179  else
1180  psize = 1
1181  nidx_max = nsize
1182  end if
1183 
1184  allocate(idx_blk(nidx_max,psize,psize))
1185  allocate(nidx( psize,psize))
1186 
1187  !$acc data copyin(lon_ref, lat_ref, hgt_ref, lon, lat, hgt) &
1188  !$acc copyout(idx_i, idx_j, hfact, idx_k, vfact) &
1189  !$acc create(idx_blk, nidx)
1190 
1191  call interp_div_block(nsize, psize, nidx_max, & ! [IN]
1192  lon_ref(:,:), lat_ref(:,:), & ! [IN]
1193  idx_blk(:,:,:), nidx(:,:), & ! [OUT]
1194  lon_min, lon_max, & ! [OUT]
1195  lat_min, lat_max, & ! [OUT]
1196  dlon, dlat ) ! [OUT]
1197 
1198  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
1199  !$omp private(ii,jj,idx_ref)
1200  !$acc kernels
1201  !$acc loop independent
1202  do j = 1, ja
1203  !$acc loop private(idx_ref,hfactt) independent
1204  do i = 1, ia
1205 
1206  ! main search
1207  call interp_search_horiz( npoints, & ! [IN]
1208  nsize, & ! [IN]
1209  lon_ref(:,:), lat_ref(:,:), & ! [IN]
1210  lon_min, lon_max, & ! [IN]
1211  lat_min, lat_max, & ! [IN]
1212  psize, nidx_max, & ! [IN]
1213  dlon, dlat, & ! [IN]
1214  idx_blk(:,:,:), nidx(:,:), & ! [IN]
1215  lon(i,j), lat(i,j), & ! [IN]
1216  idx_ref(:), hfactt(:) ) ! [OUT]
1217 
1218  !$acc loop seq
1219  do n = 1, npoints
1220  ii = mod(idx_ref(n) - 1, ia_ref) + 1
1221  jj = ( idx_ref(n) - 1 ) / ia_ref + 1
1222  idx_i(i,j,n) = ii
1223  idx_j(i,j,n) = jj
1224  hfact(i,j,n) = hfactt(n)
1225  call interp_factor1d( ka_ref, ks_ref, ke_ref, & ! [IN]
1226  ka, ks, ke, & ! [IN]
1227  hgt_ref(:,ii,jj), & ! [IN]
1228  hgt(:,i,j), & ! [IN]
1229  idx_k(:,:,i,j,n), & ! [OUT]
1230  vfact(: ,i,j,n), & ! [OUT]
1231  flag_extrap = flag_extrap ) ! [IN, optional]
1232 
1233  enddo
1234  enddo
1235  enddo
1236  !$acc end kernels
1237 
1238  !$acc end data
1239 
1240  deallocate(idx_blk, nidx)
1241 
1242  call prof_rapend ('INTERP_fact',3)
1243 
1244  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().

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 1264 of file scale_interp.F90.

1264  !$acc routine vector
1265  use scale_const, only: &
1266  undef => const_undef, &
1267  eps => const_eps
1268  implicit none
1269  integer, intent(in) :: KA_ref, KS_ref, KE_ref ! number of z-direction (reference)
1270  integer, intent(in) :: KA, KS, KE ! number of z-direction (target)
1271 
1272  integer, intent(in) :: idx_k(KA,2) ! k-index in reference
1273  real(RP), intent(in) :: vfact(KA ) ! vertical interp factor
1274  real(RP), intent(in) :: hgt_ref(KA_ref) ! height (reference)
1275  real(RP), intent(in) :: hgt (KA) ! height (target)
1276  real(RP), intent(in), target :: val_ref(KA_ref) ! value (reference)
1277 
1278  real(RP), intent(out) :: val (KA) ! value (target)
1279 
1280  logical, intent(in), optional :: spline
1281  logical, intent(in), optional :: logwgt
1282 
1283  logical :: spline_
1284  logical :: logwgt_
1285 
1286  real(RP), pointer :: work(:)
1287 
1288 #ifdef _OPENACC
1289  real(RP), intent(out), target :: workr(KA_ref,7)
1290  integer, intent(out) :: worki(KA_ref,2)
1291 #define idx_i1(k) worki(k,1)
1292 #define idx_r_i1(k) worki(k,2)
1293 #define FDZ_i1(k) workr(k,1)
1294 #define U_i1(k) workr(k,2)
1295 #else
1296  integer :: idx_i1 (KA_ref)
1297  integer :: idx_r_i1(KA_ref)
1298  real(RP) :: FDZ_i1 (KA_ref)
1299  real(RP) :: U_i1 (KA_ref)
1300 #endif
1301 
1302  integer :: kmax
1303  integer :: k
1304  !---------------------------------------------------------------------------
1305 
1306  spline_ = interp_use_spline_vert
1307  if ( present(spline) ) then
1308  spline_ = spline
1309  end if
1310 
1311  logwgt_ = .false.
1312  if ( present(logwgt) ) then
1313  logwgt_ = logwgt
1314  endif
1315 
1316  if ( logwgt_ ) then
1317 #ifdef _OPENACC
1318  work => workr(:,3)
1319 #else
1320  allocate( work(ka_ref) )
1321 #endif
1322  do k = ks_ref, ke_ref
1323  if ( abs( val_ref(k) - undef ) < eps ) then
1324  work(k) = undef
1325  else
1326  work(k) = log( val_ref(k) )
1327  end if
1328  enddo
1329  else
1330  work => val_ref
1331  endif
1332 
1333  call spline_coef( ka_ref, ks_ref, ke_ref, &
1334 #ifdef _OPENACC
1335  workr(:,4:7), &
1336 #endif
1337  hgt_ref(:), work(:), & ! (in)
1338  spline_, & ! (in)
1339  kmax, & ! (out)
1340  idx_i1(:), idx_r_i1(:), & ! (out)
1341  u_i1(:), fdz_i1(:) ) ! (out)
1342 
1343  call spline_exec( ka_ref, kmax, ka, ks, ke, &
1344  idx_k(:,:), vfact(:), & ! (in)
1345  hgt_ref(:), hgt(:), & ! (in)
1346  work(:), & ! (in)
1347  idx_i1(:), idx_r_i1(:), & ! (in)
1348  u_i1(:), fdz_i1(:), & ! (in)
1349  val(:) ) ! (out)
1350 
1351  if ( logwgt_ ) then
1352 #ifndef _OPENACC
1353  deallocate( work )
1354 #endif
1355  do k = ks, ke
1356  if ( abs( val(k) - undef ) > eps ) then
1357  val(k) = exp( val(k) )
1358  endif
1359  end do
1360  endif
1361 
1362  return

References scale_const::const_eps, scale_const::const_undef, 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 1377 of file scale_interp.F90.

1377  use scale_const, only: &
1378  undef => const_undef, &
1379  eps => const_eps
1380  implicit none
1381 
1382  integer, intent(in) :: npoints ! number of interpolation point for horizontal
1383  integer, intent(in) :: IA_ref ! number of x-direction (reference)
1384  integer, intent(in) :: JA_ref ! number of y-direction (reference)
1385  integer, intent(in) :: IA ! number of x-direction (target)
1386  integer, intent(in) :: JA ! number of y-direction (target)
1387  integer, intent(in) :: idx_i (IA,JA,npoints) ! i-index in reference
1388  integer, intent(in) :: idx_j (IA,JA,npoints) ! j-index in reference
1389  real(RP), intent(in) :: hfact (IA,JA,npoints) ! horizontal interp factor
1390  real(RP), intent(in) :: val_ref(IA_ref,JA_ref) ! value (reference)
1391 
1392  real(RP), intent(out) :: val (IA,JA) ! value (target)
1393 
1394  real(RP), intent(in), optional :: threshold_undef
1395  real(RP), intent(out), optional :: wsum(IA,JA)
1396  real(RP), intent(out), optional :: val2(IA,JA)
1397 
1398  real(RP) :: th_undef
1399 
1400  real(RP) :: fact, valn, f, w
1401  real(RP) :: sw
1402  logical :: lval2
1403 
1404  integer :: i, j, n
1405  !---------------------------------------------------------------------------
1406 
1407  call prof_rapstart('INTERP_interp',3)
1408 
1409  if ( present(threshold_undef) ) then
1410  th_undef = threshold_undef
1411  else
1412  th_undef = interp_threshold_undef
1413  end if
1414  th_undef = min( max( th_undef, eps * 2.0_rp ), 1.0_rp - eps * 2.0_rp )
1415 
1416  lval2 = present(wsum) .and. present(val2)
1417 
1418  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
1419  !$omp private(fact,valn,f,w,sw)
1420  !$acc kernels
1421 !OCL PREFETCH
1422  do j = 1, ja
1423  do i = 1, ia
1424  fact = 0.0_rp
1425  valn = 0.0_rp
1426  !$acc loop seq
1427  do n = 1, npoints
1428  f = hfact(i,j,n)
1429  w = val_ref(idx_i(i,j,n),idx_j(i,j,n))
1430  if ( f > eps .and. abs( w - undef ) > eps ) then
1431  fact = fact + f
1432  valn = valn + f * w
1433  end if
1434  end do
1435  sw = 0.5_rp - sign( 0.5_rp, fact - th_undef + eps ) ! 1.0 when fact < threshold
1436  val(i,j) = valn / ( fact + sw ) * ( 1.0_rp - sw ) + undef * sw
1437  if ( lval2 ) then
1438  wsum(i,j) = fact
1439  sw = 0.5_rp - sign( 0.5_rp, fact - eps ) ! 1.0 when fact < 0.0
1440  val2(i,j) = valn / ( fact + sw ) * ( 1.0_rp - sw ) + undef * sw
1441  end if
1442  enddo
1443  enddo
1444  !$acc end kernels
1445 
1446  call prof_rapend ('INTERP_interp',3)
1447 
1448  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 1470 of file scale_interp.F90.

1470  use scale_const, only: &
1471  undef => const_undef, &
1472  eps => const_eps
1473  implicit none
1474  integer, intent(in) :: npoints ! number of interpolation point for horizontal
1475  integer, intent(in) :: KA_ref, KS_ref, KE_ref ! number of z-direction (reference)
1476  integer, intent(in) :: IA_ref ! number of x-direction (reference)
1477  integer, intent(in) :: JA_ref ! number of y-direction (reference)
1478  integer, intent(in) :: KA, KS, KE ! number of z-direction (target)
1479  integer, intent(in) :: IA ! number of x-direction (target)
1480  integer, intent(in) :: JA ! number of y-direction (target)
1481 
1482  integer, intent(in) :: idx_i (IA,JA,npoints) ! i-index in reference
1483  integer, intent(in) :: idx_j (IA,JA,npoints) ! j-index in reference
1484  real(RP), intent(in) :: hfact (IA,JA,npoints) ! horizontal interp factor
1485  integer, intent(in) :: idx_k (KA,2,IA,JA,npoints) ! k-index in reference
1486  real(RP), intent(in) :: vfact (KA, IA,JA,npoints) ! vertical interp factor
1487  real(RP), intent(in) :: hgt_ref(KA_ref,IA_ref,JA_ref) ! height (reference)
1488  real(RP), intent(in) :: hgt (KA,IA,JA) ! height (target)
1489  real(RP), intent(in), target :: val_ref(KA_ref,IA_ref,JA_ref) ! value (reference)
1490 
1491  real(RP), intent(out) :: val (KA,IA,JA) ! value (target)
1492 
1493  logical, intent(in), optional :: spline
1494  logical, intent(in), optional :: logwgt
1495  real(RP), intent(in), optional :: threshold_undef
1496  real(RP), intent(out), optional :: wsum(KA,IA,JA)
1497  real(RP), intent(out), optional :: val2(KA,IA,JA)
1498 
1499  real(RP) :: th_undef
1500  logical :: spline_
1501  logical :: logwgt_
1502 
1503  real(RP), pointer :: work(:,:,:)
1504 
1505  integer, allocatable :: kmax(:,:)
1506  integer, allocatable :: idx(:,:,:), idx_r(:,:,:)
1507  real(RP), allocatable :: U(:,:,:), FDZ(:,:,:)
1508 
1509  real(RP) :: valn
1510  real(RP) :: fact
1511  real(RP) :: f
1512  real(RP) :: sw
1513  real(RP) :: w(KA,npoints)
1514  logical :: lval2
1515 
1516 #ifdef _OPENACC
1517  real(RP) :: workr(KA_ref,4)
1518 #endif
1519 
1520  integer :: imin, imax
1521  integer :: jmin, jmax
1522  integer :: ii, jj
1523  integer :: k, i, j, n
1524  !---------------------------------------------------------------------------
1525 
1526  call prof_rapstart('INTERP_interp',3)
1527 
1528  if ( present(threshold_undef) ) then
1529  th_undef = threshold_undef
1530  else
1531  th_undef = interp_threshold_undef
1532  end if
1533  th_undef = min( max( th_undef, eps * 2.0_rp ), 1.0_rp - eps * 2.0_rp )
1534 
1535 
1536  spline_ = interp_use_spline_vert
1537  if ( present(spline) ) then
1538  spline_ = spline
1539  end if
1540 
1541  logwgt_ = .false.
1542  if ( present(logwgt) ) then
1543  logwgt_ = logwgt
1544  endif
1545 
1546 
1547  lval2 = present(wsum) .and. present(val2)
1548 
1549 
1550  !$acc data copyin(idx_i, idx_j, hfact, idx_k, vfact, hgt_ref, hgt, val_ref) &
1551  !$acc copyout(val)
1552 
1553  !$acc data copyout(wsum) if(present(wsum))
1554  !$acc data copyout(val2) if(present(val2))
1555 
1556  imin = ia_ref
1557  jmin = ja_ref
1558  imax = 1
1559  jmax = 1
1560  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
1561  !$omp reduction(min: imin,jmin) &
1562  !$omp reduction(max: imax,jmax)
1563  !$acc kernels
1564  !$acc loop reduction(min:imin,jmin) reduction(max:imax,jmax)
1565  do n = 1, npoints
1566  !$acc loop reduction(min:imin,jmin) reduction(max:imax,jmax)
1567  do j = 1, ja
1568  !$acc loop reduction(min:imin,jmin) reduction(max:imax,jmax)
1569  do i = 1, ia
1570  imin = min(imin, idx_i(i,j,n))
1571  imax = max(imax, idx_i(i,j,n))
1572  jmin = min(jmin, idx_j(i,j,n))
1573  jmax = max(jmax, idx_j(i,j,n))
1574  end do
1575  end do
1576  end do
1577  !$acc end kernels
1578 
1579  allocate( kmax(imin:imax,jmin:jmax) )
1580  allocate( idx(ka_ref,imin:imax,jmin:jmax), idx_r(ka_ref,imin:imax,jmin:jmax) )
1581  allocate( u(ka_ref,imin:imax,jmin:jmax), fdz(ka_ref,imin:imax,jmin:jmax) )
1582  !$acc enter data create(kmax, idx, idx_r, U, FDZ)
1583 
1584  if ( logwgt_ ) then
1585  allocate( work(ka_ref,imin:imax,jmin:jmax) )
1586  !$acc enter data create(work)
1587  !$omp parallel do OMP_SCHEDULE_ collapse(2)
1588  !$acc kernels
1589  !$acc loop independent
1590  do j = jmin, jmax
1591  !$acc loop independent
1592  do i = imin, imax
1593  !$acc loop independent
1594  do k = ks_ref, ke_ref
1595  if ( abs( val_ref(k,i,j) - undef ) < eps ) then
1596  work(k,i,j) = undef
1597  else
1598  work(k,i,j) = log( val_ref(k,i,j) )
1599  end if
1600  enddo
1601  enddo
1602  enddo
1603  !$acc end kernels
1604  else
1605  work => val_ref
1606  endif
1607 
1608  !$omp parallel do OMP_SCHEDULE_ collapse(2)
1609  !$acc kernels
1610  !$acc loop independent
1611  do j = jmin, jmax
1612  !$acc loop independent private(workr)
1613  do i = imin, imax
1614  call spline_coef( ka_ref, ks_ref, ke_ref, &
1615 #ifdef _OPENACC
1616  workr(:,:), &
1617 #endif
1618  hgt_ref(:,i,j), work(:,i,j), & ! (in)
1619  spline_, & ! (in)
1620  kmax(i,j), & ! (out)
1621  idx(:,i,j), idx_r(:,i,j), & ! (out)
1622  u(:,i,j), fdz(:,i,j) ) ! (out)
1623  end do
1624  end do
1625  !$acc end kernels
1626 
1627  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
1628  !$omp private(valn,fact,w,f,ii,jj,sw)
1629  !$acc kernels
1630  !$acc loop independent
1631  do j = 1, ja
1632  !$acc loop independent private(w)
1633  do i = 1, ia
1634 
1635  !$acc loop seq
1636  do n = 1, npoints
1637  if ( hfact(i,j,n) > eps ) then
1638  ii = idx_i(i,j,n)
1639  jj = idx_j(i,j,n)
1640  call spline_exec( ka_ref, kmax(ii,jj), ka, ks, ke, &
1641  idx_k(:,:,i,j,n), vfact(:, i,j,n), & ! (in)
1642  hgt_ref(:,ii,jj), hgt(:,i,j), & ! (in)
1643  work(:,ii,jj), & ! (in)
1644  idx(:,ii,jj), idx_r(:,ii,jj), & ! (in)
1645  u(:,ii,jj), fdz(:,ii,jj), & ! (in)
1646  w(:,n) ) ! (out)
1647  end if
1648  end do
1649 
1650  do k = ks, ke
1651  fact = 0.0_rp
1652  valn = 0.0_rp
1653  !$acc loop seq
1654  do n = 1, npoints
1655  f = hfact(i,j,n)
1656  if ( f > eps .and. abs( w(k,n) - undef ) > eps ) then
1657  fact = fact + f
1658  valn = valn + f * w(k,n)
1659  endif
1660  enddo
1661  sw = 0.5_rp - sign( 0.5_rp, fact - th_undef ) ! 1.0 when fact < threshold
1662  val(k,i,j) = valn / ( fact + sw ) * ( 1.0_rp - sw ) + undef * sw
1663  if ( lval2 ) then
1664  wsum(k,i,j) = fact
1665  sw = 0.5_rp - sign( 0.5_rp, fact - eps ) ! 1.0 when fact < 0.0
1666  val2(k,i,j) = valn / ( fact + sw ) * ( 1.0_rp - sw ) + undef * sw
1667  end if
1668  enddo
1669 
1670  enddo
1671  enddo
1672  !$acc end kernels
1673 
1674  !$acc exit data delete(kmax, idx, idx_r, U, FDZ)
1675  deallocate( kmax, idx, idx_r )
1676  deallocate( u, fdz )
1677 
1678  if ( logwgt_ ) then
1679  !$acc exit data delete(work)
1680  deallocate( work )
1681  !$omp parallel do OMP_SCHEDULE_ collapse(2)
1682  !$acc kernels
1683  do j = 1, ja
1684  do i = 1, ia
1685  do k = ks, ke
1686  if ( abs( val(k,i,j) - undef ) > eps ) then
1687  val(k,i,j) = exp( val(k,i,j) )
1688  end if
1689  end do
1690  end do
1691  end do
1692  !$acc end kernels
1693  end if
1694 
1695  !$acc end data
1696  !$acc end data
1697  !$acc end data
1698 
1699  call prof_rapend ('INTERP_interp',3)
1700 
1701  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 1881 of file scale_interp.F90.

1881  !$acc routine seq
1882  use scale_const, only: &
1883  eps => const_eps, &
1884  radius => const_radius
1885  implicit none
1886 
1887  integer, intent(in) :: npoints ! number of interpolation point for horizontal
1888  integer, intent(in) :: psizex ! number of block in x-direction (reference)
1889  integer, intent(in) :: psizey ! number of block in y-direction (reference)
1890  integer, intent(in) :: IA_ref ! number of grids in x-direction (reference)
1891  integer, intent(in) :: JA_ref ! number of grids in y-direction (reference)
1892  real(RP), intent(in) :: lon_ref(IA_ref) ! longitude [rad] (reference)
1893  real(RP), intent(in) :: lat_ref(JA_ref) ! latitude [rad] (reference)
1894  real(RP), intent(in) :: lon_min ! minimum longitude
1895  real(RP), intent(in) :: lat_min ! minimum latitude
1896  real(RP), intent(in) :: dlon
1897  real(RP), intent(in) :: dlat
1898  integer, intent(in) :: i0(psizex) ! start index in the block
1899  integer, intent(in) :: i1(psizex) ! end index in the block
1900  integer, intent(in) :: j0(psizey) ! start index in the block
1901  integer, intent(in) :: j1(psizey) ! start index in the block
1902  real(RP), intent(in) :: lon ! longitude [rad] (target)
1903  real(RP), intent(in) :: lat ! latitude [rad] (target)
1904  integer, intent(out) :: idx_i(npoints) ! x-index in reference (target)
1905  integer, intent(out) :: idx_j(npoints) ! y-index in reference (target)
1906  real(RP), intent(out) :: hfact(npoints) ! horizontal interp factor (target)
1907 
1908  real(RP), intent(in), optional :: search_limit
1909  integer, intent(in), optional :: weight_order
1910 
1911  real(RP) :: drad(npoints)
1912  real(RP) :: sum
1913  real(RP) :: search_limit_
1914  integer :: weight_order_
1915 
1916  real(RP) :: lon0, lon1, lat0, lat1
1917  real(RP) :: dlon_sl, dlat_sl
1918 
1919  integer :: i, j, n
1920  integer :: ii, jj, ii0, jj0
1921  !---------------------------------------------------------------------------
1922 
1923  if ( present(search_limit) ) then
1924  search_limit_ = search_limit
1925  else
1926  search_limit_ = interp_search_limit
1927  end if
1928  search_limit_ = search_limit_ / radius ! m to radian
1929 
1930  if ( present(weight_order) ) then
1931  weight_order_ = weight_order
1932  else
1933  weight_order_ = interp_weight_order
1934  end if
1935 
1936 
1937  drad(:) = large_number
1938  idx_i(:) = 1
1939  idx_j(:) = 1
1940 
1941  ! find k-nearest points in the nearest block
1942  ii0 = max( min( int(( lon - lon_min ) / dlon) + 1, psizex ), 1)
1943  jj0 = max( min( int(( lat - lat_min ) / dlat) + 1, psizey ), 1)
1944  do j = j0(jj0), j1(jj0)
1945  do i = i0(ii0), i1(ii0)
1946  call interp_insert_2d( npoints, &
1947  lon, lat, &
1948  lon_ref(i), lat_ref(j), & ! [IN]
1949  i, j, & ! [IN]
1950  drad(:), idx_i(:), idx_j(:) ) ! [INOUT]
1951  end do
1952  end do
1953 
1954  if ( abs(drad(1)) < eps ) then
1955  hfact(:) = 0.0_rp
1956  hfact(1) = 1.0_rp
1957 
1958  return
1959  else if ( drad(1) > search_limit_ ) then
1960  hfact(:) = 0.0_rp
1961  idx_i(:) = 1 ! dummy
1962  idx_j(:) = 1 ! dummy
1963 
1964  return
1965  end if
1966 
1967  dlon_sl = max(dlon * 0.5_rp, drad(npoints))
1968  dlat_sl = max(dlat * 0.5_rp, drad(npoints))
1969  do jj = 1, psizey
1970  lat0 = lat_min + dlat * (jj-1)
1971  lat1 = lat_min + dlat * jj
1972  if ( lat < lat0 - dlat_sl &
1973  .or. lat >= lat1 + dlat_sl ) cycle
1974  do ii = 1, psizex
1975  if ( ii==ii0 .and. jj==jj0 ) cycle
1976  lon0 = lon_min + dlon * (ii-1)
1977  lon1 = lon_min + dlon * ii
1978  if ( lon < lon0 - dlon_sl &
1979  .or. lon >= lon1 + dlon_sl ) cycle
1980  do j = j0(jj0), j1(jj0)
1981  do i = i0(ii0), i1(ii0)
1982  call interp_insert_2d( npoints, &
1983  lon, lat, &
1984  lon_ref(i), lat_ref(j), & ! [IN]
1985  i, j, & ! [IN]
1986  drad(:), idx_i(:), idx_j(:) ) ! [INOUT]
1987  end do
1988  end do
1989  dlon_sl = max(dlon * 0.5_rp, drad(npoints))
1990  dlat_sl = max(dlat * 0.5_rp, drad(npoints))
1991  end do
1992  end do
1993 
1994 
1995  ! factor = 1 / drad
1996  if ( weight_order_ < 0 ) then
1997  do n = 1, npoints
1998  hfact(n) = 1.0_rp / drad(n)**(-1.0_rp/weight_order_)
1999  enddo
2000  else
2001  do n = 1, npoints
2002  hfact(n) = 1.0_rp / drad(n)**weight_order_
2003  enddo
2004  end if
2005 
2006  ! ignore far point
2007  do n = 1, npoints
2008  if ( drad(n) >= search_limit_ ) then
2009  hfact(n) = 0.0_rp
2010  endif
2011  enddo
2012 
2013  ! normalize factor
2014  sum = 0.0_rp
2015  do n = 1, npoints
2016  sum = sum + hfact(n)
2017  enddo
2018 
2019  if ( sum > 0.0_rp ) then
2020  do n = 1, npoints
2021  hfact(n) = hfact(n) / sum
2022  enddo
2023  endif
2024 
2025  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), value  i,
integer, intent(in), value  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 2076 of file scale_interp.F90.

2076  !$acc routine seq
2077  use scale_const, only: &
2078  undef => const_undef, &
2079  eps => const_eps
2080  use scale_sort, only: &
2081  sort_exec
2082 
2083  integer, intent(in) :: npoints
2084  real(RP), intent(in) :: lon, lat
2085  real(RP), intent(in) :: lon_ref, lat_ref
2086  integer, intent(in), value :: i, j
2087 
2088  real(RP), intent(inout) :: drad(npoints)
2089  integer, intent(inout) :: idx_i(npoints)
2090  integer, intent(inout) :: idx_j(npoints)
2091 
2092  real(RP) :: dradian
2093 
2094  if ( abs( lon_ref - undef ) < eps ) return
2095 
2096  dradian = haversine( lon, lat, lon_ref, lat_ref )
2097 
2098  if ( dradian <= drad(npoints) ) then
2099  ! replace last(=longest) value
2100  drad(npoints) = dradian
2101  idx_i(npoints) = i
2102  idx_j(npoints) = j
2103 
2104  ! sort by ascending order
2105  call sort_exec( npoints, & ! [IN]
2106  drad(:), & ! [INOUT]
2107  idx_i(:), & ! [INOUT]
2108  idx_j(:) ) ! [INOUT]
2109  endif
2110 
2111  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 2121 of file scale_interp.F90.

2121  use scale_const, only: &
2122  undef => const_undef, &
2123  eps => const_eps
2124  use scale_prc, only: &
2125  prc_abort
2126  integer, intent(in) :: nsize
2127  integer, intent(in) :: psize
2128  integer, intent(in) :: nidx_max
2129  real(RP), intent(in) :: lon_ref(nsize)
2130  real(RP), intent(in) :: lat_ref(nsize)
2131 
2132  integer, intent(out) :: idx (nidx_max,psize,psize)
2133  integer, intent(out) :: nidx(psize,psize)
2134  real(RP), intent(out) :: lon_min, lon_max
2135  real(RP), intent(out) :: lat_min, lat_max
2136  real(RP), intent(out) :: dlon, dlat
2137 
2138  integer :: i, ii, jj, n
2139 
2140  !$acc data copyin(lon_ref, lat_ref) copyout(idx, nidx)
2141 
2142  lon_min = 999.d0
2143  lon_max = -999.d0
2144  lat_min = 999.d0
2145  lat_max = -999.d0
2146  !$omp parallel do reduction(min: lon_min, lat_min) reduction(max: lon_max, lat_max)
2147  !$acc kernels
2148  !$acc loop reduction(min: lon_min, lat_min) reduction(max: lon_max, lat_max)
2149  do i = 1, nsize
2150  if ( abs(lon_ref(i) - undef) > eps ) then
2151  lon_min = min( lon_min, lon_ref(i) )
2152  lon_max = max( lon_max, lon_ref(i) )
2153  end if
2154  if ( abs(lat_ref(i) - undef) > eps ) then
2155  lat_min = min( lat_min, lat_ref(i) )
2156  lat_max = max( lat_max, lat_ref(i) )
2157  end if
2158  end do
2159  !$acc end kernels
2160 
2161  dlon = ( lon_max - lon_min ) / psize
2162  dlat = ( lat_max - lat_min ) / psize
2163 
2164  !$acc kernels
2165  nidx(:,:) = 0
2166  !$acc end kernels
2167  !$acc kernels
2168  !$acc loop independent
2169  do i = 1, nsize
2170  if ( abs( lon_ref(i) - undef ) < eps ) cycle
2171  ii = min(int((lon_ref(i) - lon_min) / dlon) + 1, psize)
2172  jj = min(int((lat_ref(i) - lat_min) / dlat) + 1, psize)
2173  !$acc atomic capture
2174  n = nidx(ii,jj)
2175  nidx(ii,jj) = n + 1
2176  !$acc end atomic
2177  n = n + 1
2178  if ( n <= nidx_max ) idx(n,ii,jj) = i
2179  end do
2180  !$acc end kernels
2181 
2182  !$acc end data
2183 
2184  if ( maxval(nidx) > nidx_max ) then
2185  log_error("INTERP_search_horiz",*) 'Buffer size is not enough'
2186  log_error_cont(*) ' Use larger INTERP_buffer_size_fact: ', interp_buffer_size_fact * maxval(nidx)/nidx_max
2187  call prc_abort
2188  end if
2189 
2190  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), value  x_ref0,
real(rp), intent(in), value  x_ref1,
real(rp), intent(in), value  x_ref2,
real(rp), intent(in), value  x_ref3,
real(rp), intent(in), value  y_ref0,
real(rp), intent(in), value  y_ref1,
real(rp), intent(in), value  y_ref2,
real(rp), intent(in), value  y_ref3,
real(rp), intent(in), value  x,
real(rp), intent(in), value  y,
real(rp), intent(out)  u,
real(rp), intent(out)  v,
logical, intent(out)  error 
)

Definition at line 2200 of file scale_interp.F90.

2200  !$acc routine seq
2201  implicit none
2202  real(RP), intent(in), value :: x_ref0, x_ref1, x_ref2, x_ref3
2203  real(RP), intent(in), value :: y_ref0, y_ref1, y_ref2, y_ref3
2204  real(RP), intent(in), value :: x, y
2205 
2206  real(RP), intent(out) :: u, v
2207  logical, intent(out) :: error
2208 
2209  real(RP) :: e_x, e_y
2210  real(RP) :: f_x, f_y
2211  real(RP) :: g_x, g_y
2212  real(RP) :: h_x, h_y
2213  real(RP) :: k0, k1, k2
2214  real(RP) :: w
2215  real(RP) :: sig
2216 
2217  e_x = x_ref1 - x_ref0
2218  e_y = y_ref1 - y_ref0
2219 
2220  f_x = x_ref3 - x_ref0
2221  f_y = y_ref3 - y_ref0
2222 
2223  g_x = x_ref0 - x_ref1 + x_ref2 - x_ref3
2224  g_y = y_ref0 - y_ref1 + y_ref2 - y_ref3
2225 
2226  h_x = x - x_ref0
2227  h_y = y - y_ref0
2228 
2229  k0 = cross(h_x, h_y, e_x, e_y)
2230  k1 = cross(e_x, e_y, f_x, f_y) + cross(h_x, h_y, g_x, g_y)
2231  k2 = cross(g_x, g_y, f_x, f_y)
2232 
2233  w = k1**2 - 4.0_rp * k0 * k2
2234  if ( w < 0.0_rp ) then
2235 #ifndef _OPENACC
2236  log_error("INTERP_bilinear_inv",*) 'Outside of the quadrilateral'
2237 #endif
2238  error = .true.
2239  return
2240  end if
2241 
2242  if ( abs(k1) < eps_bilinear ) then
2243 #ifndef _OPENACC
2244  log_error("INTERP_bilinear_inv",*) 'Unexpected error occured', k1
2245  log_error_cont(*) x_ref0, x_ref1, x_ref2, x_ref3
2246  log_error_cont(*) y_ref0, y_ref1, y_ref2, y_ref3
2247  log_error_cont(*) x, y
2248 #endif
2249  error = .true.
2250  return
2251  end if
2252  if ( abs(k2) < eps_bilinear * sqrt( (x_ref2-x_ref0)**2+(y_ref2-y_ref0)**2 ) ) then
2253  v = - k0 / k1
2254  else
2255  sig = sign( 1.0_rp, cross(x_ref1-x_ref0, y_ref1-y_ref0, x_ref3-x_ref0, y_ref3-y_ref0) )
2256  v = ( - k1 + sig * sqrt(w) ) / ( 2.0_rp * k2 )
2257  end if
2258  u = ( h_x - f_x * v ) / ( e_x + g_x * v )
2259 
2260  if ( u < -eps_bilinear .or. u > 1.0_rp+eps_bilinear .or. v < -eps_bilinear .or. v > 1.0_rp+eps_bilinear ) then
2261 #ifndef _OPENACC
2262  log_error("INTERP_bilinear_inv",*) 'Unexpected error occured', u, v
2263  log_error_cont(*) x_ref0, x_ref1, x_ref2, x_ref3
2264  log_error_cont(*) y_ref0, y_ref1, y_ref2, y_ref3
2265  log_error_cont(*) x, y
2266  log_error_cont(*) k0, k1, k2, sig
2267 #endif
2268  error = .true.
2269  return
2270  end if
2271  u = min( max( u, 0.0_rp ), 1.0_rp )
2272  v = min( max( v, 0.0_rp ), 1.0_rp )
2273 
2274  error = .false.
2275 
2276  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 2347 of file scale_interp.F90.

2347  !$acc routine seq
2348  implicit none
2349  real(RP), intent(in) :: x0, y0
2350  real(RP), intent(in) :: x1, y1
2351  real(RP) :: cross
2352 
2353  cross = x0 * y1 - x1 * y0
2354 

References scale_const::const_eps, scale_const::const_undef, and scale_matrix::matrix_solver_tridiagonal_1d_ta().

Referenced by interp_bilinear_inv().

Here is the call graph for this function:
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 2507 of file scale_interp.F90.

2507  !$acc routine vector
2508  use scale_const, only: &
2509  undef => const_undef, &
2510  eps => const_eps
2511  implicit none
2512  integer, intent(in) :: KA_ref, kmax
2513  integer, intent(in) :: KA, KS, KE
2514  integer, intent(in) :: idx_k(KA,2)
2515  real(RP), intent(in) :: vfact(KA)
2516  real(RP), intent(in) :: hgt_ref(KA_ref)
2517  real(RP), intent(in) :: hgt (KA)
2518  real(RP), intent(in) :: val_ref(KA_ref)
2519  integer, intent(in) :: idx (KA_ref)
2520  integer, intent(in) :: idx_r (KA_ref)
2521  real(RP), intent(in) :: U (KA_ref)
2522  real(RP), intent(in) :: FDZ (KA_ref)
2523 
2524  real(RP), intent(out) :: val(KA)
2525 
2526  real(RP) :: c1, c2, c3, d
2527  integer :: k, kk, kk2
2528 
2529  do k = ks, ke
2530  kk = idx_k(k,1)
2531  if ( kk == -1 ) then
2532  val(k) = undef
2533  else if ( idx_k(k,2) == -1 ) then
2534  val(k) = val_ref(kk)
2535  else if ( kk < idx(1) .or. kk >= idx(kmax) ) then ! linear interpolation
2536  if ( abs( val_ref(kk) - undef ) < eps .or. abs( val_ref(kk+1) - undef ) < eps ) then
2537  val(k) = undef
2538  else
2539  val(k) = val_ref(kk) * vfact(k) + val_ref(kk+1) * ( 1.0_rp - vfact(k) )
2540  end if
2541  else
2542  kk2 = idx_r(kk)
2543  kk = idx(kk2)
2544  c3 = ( u(kk2+1) - u(kk2) ) / fdz(kk2+1)
2545  c2 = 3.0_rp * u(kk2)
2546  c1 = ( val_ref(idx(kk2+1)) - val_ref(kk) ) / fdz(kk2+1) - ( 2.0_rp * u(kk2) + u(kk2+1) ) * fdz(kk2+1)
2547  d = hgt(k) - hgt_ref(kk)
2548 
2549  val(k) = ( ( c3 * d + c2 ) * d + c1 ) * d + val_ref(kk)
2550  end if
2551  end do
2552 
2553  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_atmos_grid_cartesc_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_atmos_grid_cartesC_index.F90:52
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:350
scale_atmos_grid_cartesc_index::ka
integer, public ka
Definition: scale_atmos_grid_cartesC_index.F90:47
scale_sort
module SORT
Definition: scale_sort.F90:11
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:35
scale_atmos_grid_cartesc_index::imax
integer, public imax
Definition: scale_atmos_grid_cartesC_index.F90:37
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_precision::rp
integer, parameter, public rp
Definition: scale_precision.F90:41
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_atmos_grid_cartesc_index::ia
integer, public ia
Definition: scale_atmos_grid_cartesC_index.F90:48
scale_atmos_grid_cartesc_index::kmax
integer, public kmax
Definition: scale_atmos_grid_cartesC_index.F90:36
scale_atmos_grid_cartesc_index::ja
integer, public ja
Definition: scale_atmos_grid_cartesC_index.F90:49
scale_atmos_grid_cartesc_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_atmos_grid_cartesC_index.F90:51
scale_atmos_grid_cartesc_index::jmax
integer, public jmax
Definition: scale_atmos_grid_cartesC_index.F90:38
scale_const::const_radius
real(rp), public const_radius
radius of the planet [m]
Definition: scale_const.F90:47
scale_const::const_d2r
real(rp), public const_d2r
degree to radian
Definition: scale_const.F90:33
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:43
Generated by   doxygen 1.8.17