SCALE-RM
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_factor2d (npoints, IA_ref, JA_ref, lon_ref, lat_ref, IA, JA, lon, lat, idx_i, idx_j, hfact, search_limit, latlon_structure, weight_order)
 
subroutine, public interp_factor3d (npoints, KA_ref, KS_ref, KE_ref, IA_ref, JA_ref, lon_ref, lat_ref, hgt_ref, KA, KS, KE, IA, JA, lon, lat, hgt, idx_i, idx_j, hfact, idx_k, vfact)
 
subroutine, public interp_interp2d (npoints, IA_ref, JA_ref, IA, JA, idx_i, idx_j, hfact, val_ref, val)
 
subroutine, public interp_interp3d (npoints, KA_ref, IA_ref, JA_ref, KA, KS, KE, IA, JA, idx_i, idx_j, hfact, idx_k, vfact, val_ref, val, logwgt)
 
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)
 

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

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

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

Referenced by scale_comm_cartesc_nest::comm_cartesc_nest_setup().

69  use scale_prc, only: &
70  prc_abort
71  implicit none
72 
73  integer, intent(in) :: weight_order
74  real(RP), intent(in), optional :: search_limit
75 
76  namelist /param_interp/ &
77  interp_buffer_size_fact
78 
79  integer :: ierr
80  !---------------------------------------------------------------------------
81 
82  log_newline
83  log_info("INTERP_setup",*) 'Setup'
84 
85  interp_weight_order = weight_order
86 
87  interp_search_limit = large_number
88  if ( present(search_limit) ) then
89  interp_search_limit = search_limit
90  log_info("INTERP_setup",*) 'search limit [m] : ', interp_search_limit
91  endif
92 
93  !--- read namelist
94  rewind(io_fid_conf)
95  read(io_fid_conf,nml=param_interp,iostat=ierr)
96  if( ierr < 0 ) then !--- missing
97  log_info("INTERP_setup",*) 'Not found namelist. Default used.'
98  elseif( ierr > 0 ) then !--- fatal error
99  log_error("INTERP_setup",*) 'Not appropriate names in namelist PARAM_INTERP. Check!'
100  call prc_abort
101  endif
102  log_nml(param_interp)
103 
104  return
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:55
module PROCESS
Definition: scale_prc.F90:11
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:338
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 119 of file scale_interp.F90.

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

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

119  use scale_const, only: &
120  d2r => const_d2r
121  use scale_prc, only: &
122  prc_abort
123  implicit none
124 
125  real(RP), intent(in) :: lon_org (:,:)
126  real(RP), intent(in) :: lat_org (:,:)
127  real(RP), intent(in) :: topc_org(:,:) ! full level
128  real(RP), intent(in) :: lon_loc (:,:)
129  real(RP), intent(in) :: lat_loc (:,:)
130  real(RP), intent(in) :: topc_loc(:,:) ! full level
131  real(RP), intent(in) :: topf_loc(:,:) ! half level (ceil)
132 
133  logical, intent(in), optional :: skip_z
134  logical, intent(in), optional :: skip_x
135  logical, intent(in), optional :: skip_y
136 
137  real(RP) :: max_ref, min_ref
138  real(RP) :: max_loc, min_loc
139 
140  logical :: skip_z_
141  logical :: skip_x_
142  logical :: skip_y_
143 
144  intrinsic size
145  !---------------------------------------------------------------------------
146 
147  skip_z_ = .false.
148  if ( present(skip_z) ) then
149  skip_z_ = skip_z
150  endif
151 
152  skip_x_ = .false.
153  if ( present(skip_x) ) then
154  skip_x_ = skip_x
155  endif
156 
157  skip_y_ = .false.
158  if ( present(skip_y) ) then
159  skip_y_ = skip_y
160  endif
161 
162  if ( .NOT. skip_z_ ) then
163  max_ref = maxval( topc_org(:,:) ) + minval( topf_loc(:,:) - topc_loc(:,:) ) ! not real top boundary (only for check)
164  max_loc = maxval( topc_loc(:,:) )
165 
166  if ( max_ref < max_loc ) then
167  log_error("INTERP_domain_compatibility",*) 'REQUESTED DOMAIN IS TOO MUCH BROAD'
168  log_error_cont(*) '-- VERTICAL direction over the limit'
169  log_error_cont(*) '-- reference max: ', max_ref
170  log_error_cont(*) '-- local max: ', max_loc
171  call prc_abort
172  endif
173  endif
174 
175  if ( .NOT. skip_x_ ) then
176  max_ref = maxval( lon_org(:,:) / d2r )
177  min_ref = minval( lon_org(:,:) / d2r )
178  max_loc = maxval( lon_loc(:,:) / d2r )
179  min_loc = minval( lon_loc(:,:) / d2r )
180 
181  if ( (min_ref+360.0_rp-max_ref) < 360.0_rp / size(lon_org,1) * 2.0_rp ) then
182  ! cyclic OK
183  elseif( max_ref < max_loc .OR. min_ref > min_loc ) then
184  log_error("INTERP_domain_compatibility",*) 'REQUESTED DOMAIN IS TOO MUCH BROAD'
185  log_error_cont(*) '-- LONGITUDINAL direction over the limit'
186  log_error_cont(*) '-- reference max: ', max_ref
187  log_error_cont(*) '-- reference min: ', min_ref
188  log_error_cont(*) '-- local max: ', max_loc
189  log_error_cont(*) '-- local min: ', min_loc
190  call prc_abort
191  endif
192  endif
193 
194  if ( .NOT. skip_y_ ) then
195  max_ref = maxval( lat_org(:,:) / d2r )
196  min_ref = minval( lat_org(:,:) / d2r )
197  max_loc = maxval( lat_loc(:,:) / d2r )
198  min_loc = minval( lat_loc(:,:) / d2r )
199 
200  if ( max_ref < max_loc .OR. min_ref > min_loc ) then
201  log_error("INTERP_domain_compatibility",*) 'REQUESTED DOMAIN IS TOO MUCH BROAD'
202  log_error_cont(*) '-- LATITUDINAL direction over the limit'
203  log_error_cont(*) '-- reference max: ', max_ref
204  log_error_cont(*) '-- reference min: ', min_ref
205  log_error_cont(*) '-- local max: ', max_loc
206  log_error_cont(*) '-- local min: ', min_loc
207  call prc_abort
208  endif
209  endif
210 
211  return
real(rp), public const_d2r
degree to radian
Definition: scale_const.F90:32
module PROCESS
Definition: scale_prc.F90:11
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:338
module CONSTANT
Definition: scale_const.F90:11
Here is the call graph for this function:
Here is the caller graph for this function:

◆ interp_factor2d()

subroutine, public scale_interp::interp_factor2d ( integer, intent(in)  npoints,
integer, intent(in)  IA_ref,
integer, intent(in)  JA_ref,
real(rp), dimension(ia_ref,ja_ref), intent(in)  lon_ref,
real(rp), dimension(ia_ref,ja_ref), intent(in)  lat_ref,
integer, intent(in)  IA,
integer, intent(in)  JA,
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,
integer, intent(in), optional  weight_order 
)

Definition at line 226 of file scale_interp.F90.

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

Referenced by scale_atmos_sfc_ch_rn222::atmos_sfc_ch_rn222_setup(), mod_cnvtopo::cnvtopo(), mod_copytopo::copytopo(), mod_realinput::land_interporation(), and mod_realinput::realinput_surface().

226  use scale_const, only: &
227  undef => const_undef
228  use scale_prc, only: &
229  prc_abort
230  implicit none
231 
232  integer, intent(in) :: npoints ! number of interpolation point for horizontal
233  integer, intent(in) :: ia_ref ! number of x-direction (reference)
234  integer, intent(in) :: ja_ref ! number of y-direction (reference)
235  real(RP), intent(in) :: lon_ref(ia_ref,ja_ref) ! longitude [rad] (reference)
236  real(RP), intent(in) :: lat_ref(ia_ref,ja_ref) ! latitude [rad] (reference)
237  integer, intent(in) :: ia ! number of x-direction (target)
238  integer, intent(in) :: ja ! number of y-direction (target)
239  real(RP), intent(in) :: lon (ia,ja) ! longitude [rad] (target)
240  real(RP), intent(in) :: lat (ia,ja) ! latitude [rad] (target)
241  integer, intent(out) :: idx_i(ia,ja,npoints) ! i-index in reference (target)
242  integer, intent(out) :: idx_j(ia,ja,npoints) ! j-index in reference (target)
243  real(RP), intent(out) :: hfact(ia,ja,npoints) ! horizontal interp factor (target)
244 
245  real(RP), intent(in), optional :: search_limit
246  logical, intent(in), optional :: latlon_structure
247  integer, intent(in), optional :: weight_order
248 
249  logical :: ll_struct_
250 
251  real(RP) :: lon_min, lon_max
252  real(RP) :: lat_min, lat_max
253  real(RP) :: dlon, dlat
254 
255  ! for structure grid
256  integer :: is, ie, js, je
257  integer :: psizex, psizey
258  real(RP) :: lon1d(ia_ref), lat1d(ja_ref)
259  real(RP) :: lon0, lat0
260  integer, allocatable :: i0(:), i1(:), j0(:), j1(:)
261 
262  ! for unstructure grid
263  integer :: nsize, psize, nidx_max
264  integer, allocatable :: idx_blk(:,:,:), nidx(:,:)
265  integer :: idx_ref(npoints)
266 
267 
268  integer :: i, j, ii, jj, n
269  !---------------------------------------------------------------------------
270 
271  call prof_rapstart('INTERP_fact',3)
272 
273  if ( present(latlon_structure) ) then
274  ll_struct_ = latlon_structure
275  else
276  ll_struct_ = .false.
277  end if
278 
279  hfact(:,:,:) = 0.0_rp
280 
281 
282  if ( ll_struct_ ) then
283 
284  lon1d(:) = undef
285  lon_min = undef
286  do i = 1, ia_ref
287  do j = 1, ja_ref
288  if ( lon_ref(i,j) .ne. undef ) then
289  lon1d(i) = lon_ref(i,j)
290  exit
291  end if
292  end do
293  if ( lon_min == undef .and. lon1d(i) .ne. undef ) then
294  lon_min = lon1d(i)
295  is = i
296  endif
297  if ( lon_min .ne. undef ) then
298  if ( lon1d(i) .ne. undef ) then
299  lon_max = lon1d(i)
300  ie = i
301  end if
302  end if
303  end do
304 
305  lat1d(:) = undef
306  lat_min = undef
307  do j = 1, ja_ref
308  do i = 1, ia_ref
309  if ( lat_ref(i,j) .ne. undef ) then
310  lat1d(j) = lat_ref(i,j)
311  exit
312  end if
313  end do
314  if ( lat_min == undef .and. lat1d(j) .ne. undef ) then
315  lat_min = lat1d(j)
316  js = j
317  endif
318  if ( lat_min .ne. undef ) then
319  if ( lat1d(j) .ne. undef ) then
320  lat_max = lat1d(j)
321  je = j
322  end if
323  end if
324  end do
325 
326  ! fill undef
327  dlon = ( lon_max - lon_min ) / ( ie - is )
328  do i = is, ie
329  if ( lon1d(i) == undef ) lon1d(i) = lon_min + dlon * ( i - is )
330  end do
331  dlat = ( lat_max - lat_min ) / ( je - js )
332  do j = js, je
333  if ( lat1d(j) == undef ) lat1d(j) = lat_min + dlat * ( j - js )
334  end do
335 
336 
337 
338  if ( ie-is > 10 ) then
339  psizex = int( 2.0_rp*sqrt(real(ie-is+1,rp)) )
340  else
341  psizex = 1
342  end if
343  if ( je-js > 10 ) then
344  psizey = int( 2.0_rp*sqrt(real(je-js+1,rp)) )
345  else
346  psizey = 1
347  end if
348 
349  allocate( i0(psizex), i1(psizex) )
350  allocate( j0(psizey), j1(psizey) )
351 
352  dlon = ( lon_max - lon_min ) / psizex
353  dlat = ( lat_max - lat_min ) / psizey
354 
355  do ii = 1, psizex
356  lon0 = lon_min + dlon * (ii-1)
357  do i = is, ie
358  if ( lon1d(i) >= lon0 ) then
359  i0(ii) = i
360  exit
361  end if
362  end do
363  end do
364  do ii = 1, psizex-1
365  i1(ii) = i0(ii+1) - 1
366  end do
367  i1(psizex) = ie
368 
369  do jj = 1, psizey
370  lat0 = lat_min + dlat * (jj-1)
371  do j = js, je
372  if ( lat1d(j) >= lat0 ) then
373  j0(jj) = j
374  exit
375  end if
376  end do
377  end do
378  do jj = 1, psizey-1
379  j1(jj) = j0(jj+1) - 1
380  end do
381  j1(psizey) = je
382 
383  !$omp parallel do OMP_SCHEDULE_ collapse(2)
384  do j = 1, ja
385  do i = 1, ia
386  ! main search
387  call interp_search_horiz_struct( npoints, & ! [IN]
388  psizex, psizey, & ! [IN]
389  ia_ref, ja_ref, & ! [IN]
390  lon1d(:), lat1d(:), & ! [IN]
391  lon_min, lat_min, & ! [IN]
392  dlon, dlat, & ! [IN]
393  i0(:), i1(:), j0(:), j1(:), & ! [IN]
394  lon(i,j), lat(i,j), & ! [IN]
395  idx_i(i,j,:), idx_j(i,j,:), & ! [OUT]
396  hfact(i,j,:), & ! [OUT]
397  search_limit = search_limit, & ! [IN]
398  weight_order = weight_order ) ! [IN]
399  enddo
400  enddo
401 
402  deallocate( i0, i1, j0, j1 )
403 
404  else
405 
406  nsize = ia_ref * ja_ref
407  if ( nsize > 100 ) then
408  psize = int( sqrt(2.0_rp*sqrt(real(nsize,rp))) )
409  nidx_max = nsize / psize * interp_buffer_size_fact
410  else
411  psize = 1
412  nidx_max = nsize
413  end if
414 
415  allocate(idx_blk(nidx_max,psize,psize))
416  allocate(nidx( psize,psize))
417 
418  call interp_div_block(nsize, psize, nidx_max, & ! [IN]
419  lon_ref(:,:), lat_ref(:,:), & ! [IN]
420  idx_blk(:,:,:), nidx(:,:), & ! [OUT]
421  lon_min, lon_max, & ! [OUT]
422  lat_min, lat_max, & ! [OUT]
423  dlon, dlat ) ! [OUT]
424 
425  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
426  !$omp private(idx_ref)
427  do j = 1, ja
428  do i = 1, ia
429  ! main search
430  call interp_search_horiz( npoints, & ! [IN]
431  nsize, & ! [IN]
432  lon_ref(:,:), lat_ref(:,:), & ! [IN]
433  lon_min, lon_max, & ! [IN]
434  lat_min, lat_max, & ! [IN]
435  psize, nidx_max, & ! [IN]
436  dlon, dlat, & ! [IN]
437  idx_blk(:,:,:), nidx(:,:), & ! [IN]
438  lon(i,j), lat(i,j), & ! [IN]
439  idx_ref(:), & ! [OUT]
440  hfact(i,j,:), & ! [OUT]
441  search_limit = search_limit, & ! [IN]
442  weight_order = weight_order ) ! [IN]
443  do n = 1, npoints
444  idx_i(i,j,n) = mod(idx_ref(n) - 1, ia_ref) + 1
445  idx_j(i,j,n) = ( idx_ref(n) - 1 ) / ia_ref + 1
446  end do
447  enddo
448  enddo
449 
450  deallocate(idx_blk, nidx)
451 
452  end if
453 
454  call prof_rapend ('INTERP_fact',3)
455 
456  return
integer, public ia
of whole cells: x, local, with HALO
integer, public ja
of whole cells: y, local, with HALO
real(rp), public const_undef
Definition: scale_const.F90:41
integer, public is
start point of inner domain: x, local
integer, public ie
end point of inner domain: x, local
module PROCESS
Definition: scale_prc.F90:11
integer, public je
end point of inner domain: y, local
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:338
module CONSTANT
Definition: scale_const.F90:11
integer, public js
start point of inner domain: y, local
integer, parameter, public rp
Here is the call graph for this function:
Here is the caller graph for this function:

◆ interp_factor3d()

subroutine, public scale_interp::interp_factor3d ( 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,
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,
integer, intent(in)  KA,
integer, intent(in)  KS,
integer, intent(in)  KE,
integer, intent(in)  IA,
integer, intent(in)  JA,
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,2,ia,ja,npoints), intent(out)  vfact 
)

Definition at line 484 of file scale_interp.F90.

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

Referenced by scale_comm_cartesc_nest::comm_cartesc_nest_setup(), mod_realinput::land_interporation(), and mod_realinput::realinput_surface().

484  use scale_prc, only: &
485  prc_abort
486  implicit none
487 
488  integer, intent(in) :: npoints ! number of interpolation point for horizontal
489  integer, intent(in) :: ka_ref, ks_ref, ke_ref ! number of z-direction (reference)
490  integer, intent(in) :: ia_ref ! number of x-direction (reference)
491  integer, intent(in) :: ja_ref ! number of y-direction (reference)
492  real(RP), intent(in) :: lon_ref(ia_ref,ja_ref) ! longitude [rad] (reference)
493  real(RP), intent(in) :: lat_ref(ia_ref,ja_ref) ! latitude [rad] (reference)
494  real(RP), intent(in) :: hgt_ref(ka_ref,ia_ref,ja_ref) ! height [m] (reference)
495  integer, intent(in) :: ka, ks, ke ! number of z-direction (target)
496  integer, intent(in) :: ia ! number of x-direction (target)
497  integer, intent(in) :: ja ! number of y-direction (target)
498  real(RP), intent(in) :: lon (ia,ja) ! longitude [rad] (target)
499  real(RP), intent(in) :: lat (ia,ja) ! latitude [rad] (target)
500  real(RP), intent(in) :: hgt (ka,ia,ja) ! longitude [m] (target)
501  integer, intent(out) :: idx_i(ia,ja,npoints) ! i-index in reference (target)
502  integer, intent(out) :: idx_j(ia,ja,npoints) ! j-index in reference (target)
503  real(RP), intent(out) :: hfact(ia,ja,npoints) ! horizontal interp factor (target)
504  integer, intent(out) :: idx_k(ka,2,ia,ja,npoints) ! i-index in reference (target)
505  real(RP), intent(out) :: vfact(ka,2,ia,ja,npoints) ! horizontal interp factor (target)
506 
507  integer :: nsize, psize, nidx_max
508  integer, allocatable :: idx_blk(:,:,:), nidx(:,:)
509  real(RP) :: lon_min, lon_max
510  real(RP) :: lat_min, lat_max
511  real(RP) :: dlon, dlat
512  integer :: idx_ref(npoints)
513 
514  integer :: i, j, ii, jj, n
515  !---------------------------------------------------------------------------
516 
517  call prof_rapstart('INTERP_fact',3)
518 
519  nsize = ia_ref * ja_ref
520  if ( nsize > 100 ) then
521  psize = int( sqrt(2.0_rp*sqrt(real(nsize,rp))) )
522  nidx_max = nsize / psize * interp_buffer_size_fact
523  else
524  psize = 1
525  nidx_max = nsize
526  end if
527 
528  allocate(idx_blk(nidx_max,psize,psize))
529  allocate(nidx( psize,psize))
530 
531  call interp_div_block(nsize, psize, nidx_max, & ! [IN]
532  lon_ref(:,:), lat_ref(:,:), & ! [IN]
533  idx_blk(:,:,:), nidx(:,:), & ! [OUT]
534  lon_min, lon_max, & ! [OUT]
535  lat_min, lat_max, & ! [OUT]
536  dlon, dlat ) ! [OUT]
537 
538  hfact(:,:,:) = 0.0_rp
539  vfact(:,:,:,:,:) = 0.0_rp
540 
541  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
542  !$omp private(ii,jj,idx_ref)
543  do j = 1, ja
544  do i = 1, ia
545 
546  ! main search
547  call interp_search_horiz( npoints, & ! [IN]
548  nsize, & ! [IN]
549  lon_ref(:,:), lat_ref(:,:), & ! [IN]
550  lon_min, lon_max, & ! [IN]
551  lat_min, lat_max, & ! [IN]
552  psize, nidx_max, & ! [IN]
553  dlon, dlat, & ! [IN]
554  idx_blk(:,:,:), nidx(:,:), & ! [IN]
555  lon(i,j), lat(i,j), & ! [IN]
556  idx_ref(:), hfact(i,j,:) ) ! [OUT]
557 
558  do n = 1, npoints
559  ii = mod(idx_ref(n) - 1, ia_ref) + 1
560  jj = ( idx_ref(n) - 1 ) / ia_ref + 1
561  idx_i(i,j,n) = ii
562  idx_j(i,j,n) = jj
563 
564  call interp_search_vert( ka_ref, ks_ref, ke_ref, & ! [IN]
565  ka, ks, ke, & ! [IN]
566  hgt_ref(:,ii,jj), & ! [IN]
567  hgt(:,i,j), & ! [IN]
568  idx_k(:,:,i,j,n), & ! [OUT]
569  vfact(:,:,i,j,n) ) ! [OUT]
570  enddo
571  enddo
572  enddo
573 
574  deallocate(idx_blk, nidx)
575 
576  call prof_rapend ('INTERP_fact',3)
577 
578  return
integer, public ia
of whole cells: x, local, with HALO
integer, public ja
of whole cells: y, local, with HALO
integer, public ke
end point of inner domain: z, local
module PROCESS
Definition: scale_prc.F90:11
integer, public ks
start point of inner domain: z, local
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:338
integer, public ka
of whole cells: z, local, with HALO
integer, parameter, public rp
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 
)

Definition at line 594 of file scale_interp.F90.

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_cnvtopo::cnvtopo(), mod_copytopo::copytopo(), mod_realinput::land_interporation(), and mod_realinput::realinput_surface().

594  use scale_const, only: &
595  undef => const_undef, &
596  eps => const_eps
597  implicit none
598 
599  integer, intent(in) :: npoints ! number of interpolation point for horizontal
600  integer, intent(in) :: ia_ref ! number of x-direction (reference)
601  integer, intent(in) :: ja_ref ! number of y-direction (reference)
602  integer, intent(in) :: ia ! number of x-direction (target)
603  integer, intent(in) :: ja ! number of y-direction (target)
604  integer, intent(in) :: idx_i (ia,ja,npoints) ! i-index in reference (target)
605  integer, intent(in) :: idx_j (ia,ja,npoints) ! j-index in reference (target)
606  real(RP), intent(in) :: hfact (ia,ja,npoints) ! horizontal interp factor (target)
607  real(RP), intent(in) :: val_ref(ia_ref,ja_ref) ! value (reference)
608  real(RP), intent(out) :: val (ia,ja) ! value (target)
609 
610  integer :: i, j, n
611  !---------------------------------------------------------------------------
612 
613  call prof_rapstart('INTERP_interp',3)
614 
615  !$omp parallel do OMP_SCHEDULE_ collapse(2)
616 !OCL PREFETCH
617  do j = 1, ja
618  do i = 1, ia
619  if ( hfact(i,j,1) < eps .or. val_ref(idx_i(i,j,1),idx_j(i,j,1)) == undef ) then
620  val(i,j) = undef
621  else
622  val(i,j) = hfact(i,j,1) * val_ref(idx_i(i,j,1),idx_j(i,j,1))
623  do n = 2, npoints
624  if ( val_ref(idx_i(i,j,n),idx_j(i,j,n)) == undef ) then
625  val(i,j) = val(i,j) / sum( hfact(i,j,1:n-1) )
626  exit
627  end if
628  val(i,j) = val(i,j) &
629  + hfact(i,j,n) * val_ref(idx_i(i,j,n),idx_j(i,j,n))
630  end do
631  end if
632  enddo
633  enddo
634 
635  call prof_rapend ('INTERP_interp',3)
636 
637  return
integer, public ia
of whole cells: x, local, with HALO
integer, public ja
of whole cells: y, local, with HALO
real(rp), public const_undef
Definition: scale_const.F90:41
module CONSTANT
Definition: scale_const.F90:11
real(rp), public const_eps
small number
Definition: scale_const.F90:33
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)  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,2,ia,ja,npoints), intent(in)  vfact,
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  logwgt 
)

Definition at line 658 of file scale_interp.F90.

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

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

658  implicit none
659 
660  integer, intent(in) :: npoints ! number of interpolation point for horizontal
661  integer, intent(in) :: ka_ref ! number of x-direction (reference)
662  integer, intent(in) :: ia_ref ! number of x-direction (reference)
663  integer, intent(in) :: ja_ref ! number of y-direction (reference)
664  integer, intent(in) :: ka, ks, ke ! number of z-direction (target)
665  integer, intent(in) :: ia ! number of x-direction (target)
666  integer, intent(in) :: ja ! number of y-direction (target)
667  integer, intent(in) :: idx_i (ia,ja,npoints) ! i-index in reference (target)
668  integer, intent(in) :: idx_j (ia,ja,npoints) ! j-index in reference (target)
669  real(RP), intent(in) :: hfact (ia,ja,npoints) ! horizontal interp factor (target)
670  integer, intent(in) :: idx_k (ka,2,ia,ja,npoints) ! k-index in reference (target)
671  real(RP), intent(in) :: vfact (ka,2,ia,ja,npoints) ! vertical interp factor (target)
672 
673  real(RP), intent(in), target :: val_ref(ka_ref,ia_ref,ja_ref) ! value (reference)
674 
675  real(RP), intent(out) :: val (ka,ia,ja) ! value (target)
676 
677  logical, intent(in), optional:: logwgt ! use logarithmic weighted interpolation?
678 
679  logical :: logwgt_
680 
681  real(RP), pointer :: work(:,:,:)
682 
683  integer :: k, i, j, n
684  !---------------------------------------------------------------------------
685 
686  call prof_rapstart('INTERP_interp',3)
687 
688  logwgt_ = .false.
689  if ( present(logwgt) ) then
690  logwgt_ = logwgt
691  endif
692 
693  if ( logwgt_ ) then
694  allocate( work(ka_ref,ia_ref,ja_ref) )
695  !$omp parallel do OMP_SCHEDULE_ collapse(2)
696  do j = 1, ja_ref
697  do i = 1, ia_ref
698  do k = 1, ka_ref
699  work(k,i,j) = log( val_ref(k,i,j) )
700  enddo
701  enddo
702  enddo
703  else
704  work => val_ref
705  endif
706 
707  !$omp parallel do OMP_SCHEDULE_ collapse(2)
708  do j = 1, ja
709  do i = 1, ia
710  do k = ks, ke
711  val(k,i,j) = hfact(i,j,1) * vfact(k,1,i,j,1) * work(idx_k(k,1,i,j,1),idx_i(i,j,1),idx_j(i,j,1)) &
712  + hfact(i,j,1) * vfact(k,2,i,j,1) * work(idx_k(k,2,i,j,1),idx_i(i,j,1),idx_j(i,j,1))
713  enddo
714  enddo
715  enddo
716 
717 !OCL SERIAL
718  do n = 2, npoints
719  !$omp parallel do OMP_SCHEDULE_ collapse(2)
720 !OCL PREFETCH
721  do j = 1, ja
722  do i = 1, ia
723  do k = ks, ke
724  val(k,i,j) = val(k,i,j) &
725  + hfact(i,j,n) * vfact(k,1,i,j,n) * work(idx_k(k,1,i,j,n),idx_i(i,j,n),idx_j(i,j,n)) &
726  + hfact(i,j,n) * vfact(k,2,i,j,n) * work(idx_k(k,2,i,j,n),idx_i(i,j,n),idx_j(i,j,n))
727  enddo
728  enddo
729  enddo
730  enddo
731 
732  if ( logwgt_ ) then
733  deallocate( work )
734  !$omp parallel do OMP_SCHEDULE_ collapse(2)
735  do j = 1, ja
736  do i = 1, ia
737  do k = ks, ke
738  val(k,i,j) = exp( val(k,i,j) )
739  end do
740  end do
741  end do
742  endif
743 
744  call prof_rapend ('INTERP_interp',3)
745 
746  return
integer, public ia
of whole cells: x, local, with HALO
integer, public ja
of whole cells: y, local, with HALO
integer, public ke
end point of inner domain: z, local
integer, public ks
start point of inner domain: z, local
integer, public ka
of whole cells: z, local, with HALO
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 921 of file scale_interp.F90.

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

Referenced by interp_factor2d().

921  use scale_const, only: &
922  eps => const_eps, &
923  radius => const_radius
924  implicit none
925 
926  integer, intent(in) :: npoints ! number of interpolation point for horizontal
927  integer, intent(in) :: psizex ! number of block in x-direction (reference)
928  integer, intent(in) :: psizey ! number of block in y-direction (reference)
929  integer, intent(in) :: ia_ref ! number of grids in x-direction (reference)
930  integer, intent(in) :: ja_ref ! number of grids in y-direction (reference)
931  real(RP), intent(in) :: lon_ref(ia_ref) ! longitude [rad] (reference)
932  real(RP), intent(in) :: lat_ref(ja_ref) ! latitude [rad] (reference)
933  real(RP), intent(in) :: lon_min ! minimum longitude
934  real(RP), intent(in) :: lat_min ! minimum latitude
935  real(RP), intent(in) :: dlon
936  real(RP), intent(in) :: dlat
937  integer, intent(in) :: i0(psizex) ! start index in the block
938  integer, intent(in) :: i1(psizex) ! end index in the block
939  integer, intent(in) :: j0(psizey) ! start index in the block
940  integer, intent(in) :: j1(psizey) ! start index in the block
941  real(RP), intent(in) :: lon ! longitude [rad] (target)
942  real(RP), intent(in) :: lat ! latitude [rad] (target)
943  integer, intent(out) :: idx_i(npoints) ! x-index in reference (target)
944  integer, intent(out) :: idx_j(npoints) ! y-index in reference (target)
945  real(RP), intent(out) :: hfact(npoints) ! horizontal interp factor (target)
946 
947  real(RP), intent(in), optional :: search_limit
948  integer, intent(in), optional :: weight_order
949 
950  real(RP) :: drad(npoints)
951  real(RP) :: sum
952  real(RP) :: search_limit_
953  integer :: weight_order_
954 
955  real(RP) :: lon0, lon1, lat0, lat1
956  real(RP) :: dlon_sl, dlat_sl
957 
958  integer :: i, j, n
959  integer :: ii, jj, ii0, jj0
960  !---------------------------------------------------------------------------
961 
962  if ( present(search_limit) ) then
963  search_limit_ = search_limit
964  else
965  search_limit_ = interp_search_limit
966  end if
967  search_limit_ = search_limit_ / radius ! m to radian
968 
969  if ( present(weight_order) ) then
970  weight_order_ = weight_order
971  else
972  weight_order_ = interp_weight_order
973  end if
974 
975 
976  drad(:) = large_number
977  idx_i(:) = 1
978  idx_j(:) = 1
979 
980  ! find k-nearest points in the nearest block
981  ii0 = max( min( int(( lon - lon_min ) / dlon) + 1, psizex ), 1)
982  jj0 = max( min( int(( lat - lat_min ) / dlat) + 1, psizey ), 1)
983  do j = j0(jj0), j1(jj0)
984  do i = i0(ii0), i1(ii0)
985  call interp_insert_2d( npoints, &
986  lon, lat, &
987  lon_ref(i), lat_ref(j), & ! [IN]
988  i, j, & ! [IN]
989  drad(:), idx_i(:), idx_j(:) ) ! [INOUT]
990  end do
991  end do
992 
993  if ( abs(drad(1)) < eps ) then
994  hfact(:) = 0.0_rp
995  hfact(1) = 1.0_rp
996 
997  return
998  else if ( drad(1) > search_limit_ ) then
999  hfact(:) = 0.0_rp
1000  idx_i(:) = 1 ! dummy
1001  idx_j(:) = 1 ! dummy
1002 
1003  return
1004  end if
1005 
1006  dlon_sl = max(dlon * 0.5_rp, drad(npoints))
1007  dlat_sl = max(dlat * 0.5_rp, drad(npoints))
1008  do jj = 1, psizey
1009  lat0 = lat_min + dlat * (jj-1)
1010  lat1 = lat_min + dlat * jj
1011  if ( lat < lat0 - dlat_sl &
1012  .or. lat >= lat1 + dlat_sl ) cycle
1013  do ii = 1, psizex
1014  if ( ii==ii0 .and. jj==jj0 ) cycle
1015  lon0 = lon_min + dlon * (ii-1)
1016  lon1 = lon_min + dlon * ii
1017  if ( lon < lon0 - dlon_sl &
1018  .or. lon >= lon1 + dlon_sl ) cycle
1019  do j = j0(jj0), j1(jj0)
1020  do i = i0(ii0), i1(ii0)
1021  call interp_insert_2d( npoints, &
1022  lon, lat, &
1023  lon_ref(i), lat_ref(j), & ! [IN]
1024  i, j, & ! [IN]
1025  drad(:), idx_i(:), idx_j(:) ) ! [INOUT]
1026  end do
1027  end do
1028  dlon_sl = max(dlon * 0.5_rp, drad(npoints))
1029  dlat_sl = max(dlat * 0.5_rp, drad(npoints))
1030  end do
1031  end do
1032 
1033 
1034  ! factor = 1 / drad
1035  if ( weight_order_ < 0 ) then
1036  do n = 1, npoints
1037  hfact(n) = 1.0_rp / drad(n)**(-1.0_rp/weight_order_)
1038  enddo
1039  else
1040  do n = 1, npoints
1041  hfact(n) = 1.0_rp / drad(n)**weight_order_
1042  enddo
1043  end if
1044 
1045  ! ignore far point
1046  do n = 1, npoints
1047  if ( drad(n) >= search_limit_ ) then
1048  hfact(n) = 0.0_rp
1049  endif
1050  enddo
1051 
1052  ! normalize factor
1053  sum = 0.0_rp
1054  do n = 1, npoints
1055  sum = sum + hfact(n)
1056  enddo
1057 
1058  if ( sum > 0.0_rp ) then
1059  do n = 1, npoints
1060  hfact(n) = hfact(n) / sum
1061  enddo
1062  endif
1063 
1064  return
real(rp), public const_radius
radius of the planet [m]
Definition: scale_const.F90:44
module CONSTANT
Definition: scale_const.F90:11
real(rp), public const_eps
small number
Definition: scale_const.F90:33
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 1181 of file scale_interp.F90.

References scale_const::const_undef.

Referenced by interp_search_horiz_struct().

1181  use scale_const, only: &
1182  undef => const_undef
1183  use scale_sort, only: &
1184  sort_exec
1185 
1186  integer, intent(in) :: npoints
1187  real(RP), intent(in) :: lon, lat
1188  real(RP), intent(in) :: lon_ref, lat_ref
1189  integer, intent(in) :: i, j
1190 
1191  real(RP), intent(inout) :: drad(npoints)
1192  integer, intent(inout) :: idx_i(npoints)
1193  integer, intent(inout) :: idx_j(npoints)
1194 
1195  real(RP) :: dradian
1196 
1197  if ( lon_ref == undef ) return
1198 
1199  dradian = haversine( lon, lat, lon_ref, lat_ref )
1200 
1201  if ( dradian <= drad(npoints) ) then
1202  ! replace last(=longest) value
1203  drad(npoints) = dradian
1204  idx_i(npoints) = i
1205  idx_j(npoints) = j
1206 
1207  ! sort by ascending order
1208  call sort_exec( npoints, & ! [IN]
1209  drad(:), & ! [INOUT]
1210  idx_i(:), & ! [INOUT]
1211  idx_j(:) ) ! [INOUT]
1212  endif
1213 
1214  return
real(rp), public const_undef
Definition: scale_const.F90:41
module CONSTANT
Definition: scale_const.F90:11
module SORT
Definition: scale_sort.F90:11
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 1223 of file scale_interp.F90.

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

Referenced by interp_factor2d(), and interp_factor3d().

1223  use scale_const, only: &
1224  undef => const_undef
1225  use scale_prc, only: &
1226  prc_abort
1227  integer, intent(in) :: nsize
1228  integer, intent(in) :: psize
1229  integer, intent(in) :: nidx_max
1230  real(RP), intent(in) :: lon_ref(nsize)
1231  real(RP), intent(in) :: lat_ref(nsize)
1232 
1233  integer, intent(out) :: idx (nidx_max,psize,psize)
1234  integer, intent(out) :: nidx(psize,psize)
1235  real(RP), intent(out) :: lon_min, lon_max
1236  real(RP), intent(out) :: lat_min, lat_max
1237  real(RP), intent(out) :: dlon, dlat
1238 
1239  integer :: i, ii, jj, n
1240 
1241  lon_min = minval(lon_ref(:), mask=lon_ref.ne.undef)
1242  lon_max = maxval(lon_ref(:), mask=lon_ref.ne.undef)
1243  lat_min = minval(lat_ref(:), mask=lat_ref.ne.undef)
1244  lat_max = maxval(lat_ref(:), mask=lat_ref.ne.undef)
1245 
1246  dlon = ( lon_max - lon_min ) / psize
1247  dlat = ( lat_max - lat_min ) / psize
1248 
1249  nidx(:,:) = 0
1250  do i = 1, nsize
1251  if ( lon_ref(i) == undef ) cycle
1252  ii = min(int((lon_ref(i) - lon_min) / dlon) + 1, psize)
1253  jj = min(int((lat_ref(i) - lat_min) / dlat) + 1, psize)
1254  n = nidx(ii,jj) + 1
1255  nidx(ii,jj) = n
1256  if ( n <= nidx_max ) idx(n,ii,jj) = i
1257  end do
1258 
1259  if ( maxval(nidx) > nidx_max ) then
1260  log_error("INTERP_search_horiz",*) 'Buffer size is not enough'
1261  log_error_cont(*) ' Use larger INTERP_buffer_size_fact: ', interp_buffer_size_fact * maxval(nidx)/nidx_max
1262  call prc_abort
1263  end if
1264 
1265  return
real(rp), public const_undef
Definition: scale_const.F90:41
module PROCESS
Definition: scale_prc.F90:11
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:338
module CONSTANT
Definition: scale_const.F90:11
Here is the call graph for this function:
Here is the caller graph for this function: