SCALE-RM
scale_interp.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
10 !-------------------------------------------------------------------------------
11 #include "scalelib.h"
12 module scale_interp
13  !-----------------------------------------------------------------------------
14  !
15  !++ used modules
16  !
17  use scale_precision
18  use scale_io
19  use scale_prof
20  !-----------------------------------------------------------------------------
21  implicit none
22  private
23  !-----------------------------------------------------------------------------
24  !
25  !++ Public procedure
26  !
27  public :: interp_setup
28  public :: interp_domain_compatibility
29  public :: interp_factor1d
30  public :: interp_factor2d
31  public :: interp_factor3d
32  public :: interp_factor2d_linear_latlon
33  public :: interp_factor2d_linear_xy
34  public :: interp_factor2d_weight
35  public :: interp_factor3d_linear_latlon
36  public :: interp_factor3d_linear_xy
37  public :: interp_factor3d_weight
38 
39  public :: interp_interp1d
40  public :: interp_interp2d
41  public :: interp_interp3d
42 
43  interface interp_factor2d
44  procedure INTERP_factor2d_linear_latlon
45  procedure INTERP_factor2d_linear_xy
46  procedure INTERP_factor2d_weight
47  end interface interp_factor2d
48  interface interp_factor3d
49  procedure INTERP_factor3d_linear_latlon
50  procedure INTERP_factor3d_linear_xy
51  procedure INTERP_factor3d_weight
52  end interface interp_factor3d
53 
54  !-----------------------------------------------------------------------------
55  !
56  !++ Public parameters & variables
57  !
58  !-----------------------------------------------------------------------------
59  !
60  !++ Private procedure
61  !
62  private :: interp_search_horiz
63  private :: interp_insert
64 
65  private :: haversine
66 
67  !-----------------------------------------------------------------------------
68  !
69  !++ Private parameters & variables
70  !
71  real(RP), private, parameter :: large_number = 9.999e+15_rp
72 
73  integer, private :: INTERP_weight_order = 2
74  real(RP), private :: INTERP_search_limit
75  real(RP), private :: INTERP_buffer_size_fact = 2.0_rp
76  logical, private :: INTERP_use_spline_vert = .true.
77  real(RP), private :: INTERP_threshold_undef = 1.0_rp ! [0-1]
78 
79  real(RP), private :: EPS_bilinear
80 
81  !$acc declare create(INTERP_weight_order, INTERP_search_limit, INTERP_use_spline_vert, EPS_bilinear)
82 
83  !-----------------------------------------------------------------------------
84 contains
85  !-----------------------------------------------------------------------------
87  subroutine interp_setup( &
88  weight_order, &
89  search_limit )
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
136  end subroutine interp_setup
137 
138  !-----------------------------------------------------------------------------
139  subroutine interp_domain_compatibility( &
140  lon_org, &
141  lat_org, &
142  topc_org, &
143  lon_loc, &
144  lat_loc, &
145  topc_loc, &
146  topf_loc, &
147  skip_x, &
148  skip_y, &
149  skip_z )
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
243  end subroutine interp_domain_compatibility
244 
245  !-----------------------------------------------------------------------------
246  ! vertical search of interpolation points for two-points
247 !OCL SERIAL
248  subroutine interp_factor1d( &
249  KA_ref, KS_ref, KE_ref, &
250  KA, KS, KE, &
251  hgt_ref, &
252  hgt, &
253  idx_k, &
254  vfact, &
255  flag_extrap )
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
345  end subroutine interp_factor1d
346 
347 
348  !-----------------------------------------------------------------------------
349  ! make interpolation factor using bi-linear method on the latlon coordinate
350  ! This can be used only for structured grid
351  subroutine interp_factor2d_linear_latlon( &
352  IA_ref, JA_ref, &
353  IA, JA, &
354  lon_ref, lat_ref, &
355  lon, lat, &
356  idx_i, idx_j, hfact )
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
469  end subroutine interp_factor2d_linear_latlon
470 
471  !-----------------------------------------------------------------------------
472  ! make interpolation factor using bi-linear method on the xy coordinate
473  ! This can be used only for structured grid
474  subroutine interp_factor2d_linear_xy( &
475  IA_ref, JA_ref, &
476  IA, JA, &
477  x_ref, y_ref, &
478  x, y, &
479  idx_i, idx_j, &
480  hfact, &
481  zonal, pole, &
482  missing )
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
733  end subroutine interp_factor2d_linear_xy
734 
735  !-----------------------------------------------------------------------------
736  ! make interpolation factor using Lat-Lon with weight based on the distance
737  subroutine interp_factor2d_weight( &
738  npoints, &
739  IA_ref, JA_ref, &
740  IA, JA, &
741  lon_ref,lat_ref, &
742  lon, lat, &
743  idx_i, idx_j, hfact, &
744  search_limit, &
745  latlon_structure, &
746  lon_1d, lat_1d, &
747  weight_order )
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  !$acc kernels copyin(lon_1d, lat_1d, i0, i1, j0, j1, lon, lat) copyout(idx_i, idx_j, hfact)
852  !$acc loop independent
853  do j = 1, ja
854  !$acc loop private(idx_it,idx_jt,hfactt) independent
855  do i = 1, ia
856  ! main search
857  call interp_search_horiz_struct( npoints, & ! [IN]
858  psizex, psizey, & ! [IN]
859  ia_ref, ja_ref, & ! [IN]
860  lon_1d(:), lat_1d(:), & ! [IN]
861  lon_min, lat_min, & ! [IN]
862  dlon, dlat, & ! [IN]
863  i0(:), i1(:), j0(:), j1(:), & ! [IN]
864  lon(i,j), lat(i,j), & ! [IN]
865  idx_it(:), idx_jt(:), & ! [OUT]
866  hfactt(:), & ! [OUT]
867  search_limit = search_limit, & ! [IN]
868  weight_order = weight_order ) ! [IN]
869  do n = 1, npoints
870  idx_i(i,j,n) = idx_it(n)
871  idx_j(i,j,n) = idx_jt(n)
872  hfact(i,j,n) = hfactt(n)
873  end do
874  enddo
875  enddo
876  !$acc end kernels
877 
878  deallocate( i0, i1, j0, j1 )
879 
880  else
881 
882  nsize = ia_ref * ja_ref
883  if ( nsize > 100 ) then
884  psize = int( sqrt(2.0_rp*sqrt(real(nsize,rp))) )
885  nidx_max = nsize / psize * interp_buffer_size_fact
886  else
887  psize = 1
888  nidx_max = nsize
889  end if
890 
891  allocate(idx_blk(nidx_max,psize,psize))
892  allocate(nidx( psize,psize))
893 
894  !$acc data copyin(lon_ref, lat_ref, lon, lat) copyout(idx_i, idx_j, hfact) create(idx_blk, nidx)
895 
896  call interp_div_block(nsize, psize, nidx_max, & ! [IN]
897  lon_ref(:,:), lat_ref(:,:), & ! [IN]
898  idx_blk(:,:,:), nidx(:,:), & ! [OUT]
899  lon_min, lon_max, & ! [OUT]
900  lat_min, lat_max, & ! [OUT]
901  dlon, dlat ) ! [OUT]
902 
903  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
904  !$omp private(idx_ref)
905  !$acc kernels
906  !$acc loop independent
907  do j = 1, ja
908  !$acc loop private(idx_ref,hfactt) independent
909  do i = 1, ia
910  ! main search
911  call interp_search_horiz( npoints, & ! [IN]
912  nsize, & ! [IN]
913  lon_ref(:,:), lat_ref(:,:), & ! [IN]
914  lon_min, lon_max, & ! [IN]
915  lat_min, lat_max, & ! [IN]
916  psize, nidx_max, & ! [IN]
917  dlon, dlat, & ! [IN]
918  idx_blk(:,:,:), nidx(:,:), & ! [IN]
919  lon(i,j), lat(i,j), & ! [IN]
920  idx_ref(:), & ! [OUT]
921  hfactt(:), & ! [OUT]
922  search_limit = search_limit, & ! [IN]
923  weight_order = weight_order ) ! [IN]
924  do n = 1, npoints
925  idx_i(i,j,n) = mod(idx_ref(n) - 1, ia_ref) + 1
926  idx_j(i,j,n) = ( idx_ref(n) - 1 ) / ia_ref + 1
927  hfact(i,j,n) = hfactt(n)
928  end do
929  enddo
930  enddo
931  !$acc end kernels
932 
933  !$acc end data
934 
935  deallocate(idx_blk, nidx)
936 
937  end if
938 
939  call prof_rapend ('INTERP_fact',3)
940 
941  return
942  end subroutine interp_factor2d_weight
943 
944  !-----------------------------------------------------------------------------
945  ! make interpolation factor using bi-linear method for the horizontal direction on latlon grid
946  ! This can be used only for structured grid
947  subroutine interp_factor3d_linear_latlon( &
948  KA_ref, KS_ref, KE_ref, &
949  IA_ref, JA_ref, &
950  KA, KS, KE, &
951  IA, JA, &
952  lon_ref, lat_ref, &
953  hgt_ref, &
954  lon, lat, &
955  hgt, &
956  idx_i, idx_j, &
957  hfact, &
958  idx_k, &
959  vfact, &
960  flag_extrap )
961  implicit none
962 
963  integer, intent(in) :: ka_ref, ks_ref, ke_ref ! number of z-direction (reference)
964  integer, intent(in) :: ia_ref ! number of x-direction (reference)
965  integer, intent(in) :: ja_ref ! number of y-direction (reference)
966  integer, intent(in) :: ka, ks, ke ! number of z-direction (target)
967  integer, intent(in) :: ia ! number of x-direction (target)
968  integer, intent(in) :: ja ! number of y-direction (target)
969  real(rp), intent(in) :: lon_ref(ia_ref) ! longitude (reference)
970  real(rp), intent(in) :: lat_ref(ja_ref) ! latitude (reference)
971  real(rp), intent(in) :: hgt_ref(ka_ref,ia_ref,ja_ref) ! height [m] (reference)
972  real(rp), intent(in) :: lon(ia,ja) ! longitude (target)
973  real(rp), intent(in) :: lat(ia,ja) ! latitude (target)
974  real(rp), intent(in) :: hgt(ka,ia,ja) ! longitude [m] (target)
975 
976  integer, intent(out) :: idx_i(ia,ja,4) ! i-index in reference (target)
977  integer, intent(out) :: idx_j(ia,ja,4) ! j-index in reference (target)
978  real(rp), intent(out) :: hfact(ia,ja,4) ! horizontal interp factor (target)
979  integer, intent(out) :: idx_k(ka,2,ia,ja,4) ! k-index in reference (target)
980  real(rp), intent(out) :: vfact(ka, ia,ja,4) ! vertical interp factor (target)
981 
982  logical, intent(in), optional :: flag_extrap ! when true, vertical extrapolation will be executed (just copy)
983 
984  integer :: i, j, ii, jj, n
985 
986  !$acc data copyin(lon_ref, lat_ref, hgt_ref, lon, lat, hgt) copyout(idx_i, idx_j, hfact, idx_k, vfact)
987 
988  call interp_factor2d_linear_latlon( &
989  ia_ref, ja_ref, & ! [IN]
990  ia, ja, & ! [IN]
991  lon_ref(:), lat_ref(:), & ! [IN]
992  lon(:,:), lat(:,:), & ! [IN]
993  idx_i(:,:,:), idx_j(:,:,:), hfact(:,:,:) ) ! [OUT]
994 
995  call prof_rapstart('INTERP_fact',3)
996 
997  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
998  !$omp private(ii,jj)
999  !$acc kernels
1000  !$acc loop independent
1001  do j = 1, ja
1002  !$acc loop independent
1003  do i = 1, ia
1004  !$acc loop seq
1005  do n = 1, 4
1006  ii = idx_i(i,j,n)
1007  jj = idx_j(i,j,n)
1008  call interp_factor1d( ka_ref, ks_ref, ke_ref, & ! [IN]
1009  ka, ks, ke, & ! [IN]
1010  hgt_ref(:,ii,jj), & ! [IN]
1011  hgt(:,i,j), & ! [IN]
1012  idx_k(:,:,i,j,n), & ! [OUT]
1013  vfact(:, i,j,n), & ! [OUT]
1014  flag_extrap = flag_extrap ) ! [IN, optional]
1015 
1016  enddo
1017  enddo
1018  enddo
1019  !$acc end kernels
1020 
1021  !$acc end data
1022 
1023  call prof_rapend('INTERP_fact',3)
1024 
1025  return
1026  end subroutine interp_factor3d_linear_latlon
1027 
1028  !-----------------------------------------------------------------------------
1029  ! make interpolation factor using bi-linear method for the horizontal direction on the xy coordinate
1030  ! This can be used only for structured grid
1031  subroutine interp_factor3d_linear_xy( &
1032  KA_ref, KS_ref, KE_ref, &
1033  IA_ref, JA_ref, &
1034  KA, KS, KE, &
1035  IA, JA, &
1036  x_ref, y_ref, &
1037  hgt_ref, &
1038  x, y, &
1039  hgt, &
1040  idx_i, idx_j, &
1041  hfact, &
1042  idx_k, &
1043  vfact, &
1044  flag_extrap, &
1045  zonal, pole, &
1046  missing )
1047  implicit none
1048  integer, intent(in) :: ka_ref, ks_ref, ke_ref ! number of z-direction (reference)
1049  integer, intent(in) :: ia_ref ! number of x-direction (reference)
1050  integer, intent(in) :: ja_ref ! number of y-direction (reference)
1051  integer, intent(in) :: ka, ks, ke ! number of z-direction (target)
1052  integer, intent(in) :: ia ! number of x-direction (target)
1053  integer, intent(in) :: ja ! number of y-direction (target)
1054  real(rp), intent(in) :: x_ref(ia_ref,ja_ref) ! x point (reference)
1055  real(rp), intent(in) :: y_ref(ia_ref,ja_ref) ! y point (reference)
1056  real(rp), intent(in) :: hgt_ref(ka_ref,ia_ref,ja_ref) ! height [m] (reference)
1057  real(rp), intent(in) :: x (ia) ! x point (target)
1058  real(rp), intent(in) :: y (ja) ! y point (target)
1059  real(rp), intent(in) :: hgt(ka,ia,ja) ! longitude [m] (target)
1060 
1061  integer, intent(out) :: idx_i(ia,ja,4) ! i-index in reference
1062  integer, intent(out) :: idx_j(ia,ja,4) ! j-index in reference
1063  real(rp), intent(out) :: hfact(ia,ja,4) ! horizontal interp factor
1064  integer, intent(out) :: idx_k(ka,2,ia,ja,4) ! k-index in reference
1065  real(rp), intent(out) :: vfact(ka, ia,ja,4) ! vertical interp factor
1066 
1067  logical, intent(in), optional :: flag_extrap ! when true, vertical extrapolation will be executed (just copy)
1068  logical, intent(in), optional :: zonal
1069  logical, intent(in), optional :: pole
1070  logical, intent(in), optional :: missing
1071 
1072  integer :: i, j, ii, jj, n
1073 
1074  !$acc data copyin(x_ref, y_ref, hgt_ref, x, y, hgt) copyout(idx_i, idx_j, hfact, idx_k, vfact)
1075 
1076  call interp_factor2d_linear_xy( &
1077  ia_ref, ja_ref, & ! [IN]
1078  ia, ja, & ! [IN]
1079  x_ref(:,:), y_ref(:,:), & ! [IN]
1080  x(:), y(:), & ! [IN]
1081  idx_i(:,:,:), idx_j(:,:,:), hfact(:,:,:), & ! [OUT]
1082  zonal = zonal, pole = pole, & ! [IN]
1083  missing = missing ) ! [IN]
1084 
1085  call prof_rapstart('INTERP_fact',3)
1086 
1087  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
1088  !$omp private(ii,jj)
1089  !$acc kernels
1090  !$acc loop independent
1091  do j = 1, ja
1092  !$acc loop independent
1093  do i = 1, ia
1094  !$acc loop seq
1095  do n = 1, 4
1096  ii = idx_i(i,j,n)
1097  jj = idx_j(i,j,n)
1098  call interp_factor1d( ka_ref, ks_ref, ke_ref, & ! [IN]
1099  ka, ks, ke, & ! [IN]
1100  hgt_ref(:,ii,jj), & ! [IN]
1101  hgt(:,i,j), & ! [IN]
1102  idx_k(:,:,i,j,n), & ! [OUT]
1103  vfact(: ,i,j,n), & ! [OUT]
1104  flag_extrap = flag_extrap ) ! [IN, optional]
1105 
1106  enddo
1107  enddo
1108  enddo
1109  !$acc end kernels
1110 
1111  !$acc end data
1112 
1113  call prof_rapend('INTERP_fact',3)
1114 
1115  return
1116  end subroutine interp_factor3d_linear_xy
1117 
1118  !-----------------------------------------------------------------------------
1119  ! make interpolation factor using Lat-Lon and Z-Height information
1120  subroutine interp_factor3d_weight( &
1121  npoints, &
1122  KA_ref, KS_ref, KE_ref, &
1123  IA_ref, JA_ref, &
1124  KA, KS, KE, &
1125  IA, JA, &
1126  lon_ref, lat_ref, &
1127  hgt_ref, &
1128  lon, lat, &
1129  hgt, &
1130  idx_i, idx_j, &
1131  hfact, &
1132  idx_k, &
1133  vfact, &
1134  flag_extrap )
1135  use scale_prc, only: &
1136  prc_abort
1137  implicit none
1138  integer, intent(in) :: npoints ! number of interpolation point for horizontal
1139  integer, intent(in) :: ka_ref, ks_ref, ke_ref ! number of z-direction (reference)
1140  integer, intent(in) :: ia_ref ! number of x-direction (reference)
1141  integer, intent(in) :: ja_ref ! number of y-direction (reference)
1142  integer, intent(in) :: ka, ks, ke ! number of z-direction (target)
1143  integer, intent(in) :: ia ! number of x-direction (target)
1144  integer, intent(in) :: ja ! number of y-direction (target)
1145 
1146  real(rp), intent(in) :: lon_ref(ia_ref,ja_ref) ! longitude [rad] (reference)
1147  real(rp), intent(in) :: lat_ref(ia_ref,ja_ref) ! latitude [rad] (reference)
1148  real(rp), intent(in) :: hgt_ref(ka_ref,ia_ref,ja_ref) ! height [m] (reference)
1149  real(rp), intent(in) :: lon (ia,ja) ! longitude [rad] (target)
1150  real(rp), intent(in) :: lat (ia,ja) ! latitude [rad] (target)
1151  real(rp), intent(in) :: hgt (ka,ia,ja) ! longitude [m] (target)
1152 
1153  integer, intent(out) :: idx_i(ia,ja,npoints) ! i-index in reference
1154  integer, intent(out) :: idx_j(ia,ja,npoints) ! j-index in reference
1155  real(rp), intent(out) :: hfact(ia,ja,npoints) ! horizontal interp factor
1156  integer, intent(out) :: idx_k(ka,2,ia,ja,npoints) ! k-index in reference
1157  real(rp), intent(out) :: vfact(ka, ia,ja,npoints) ! vertical interp factor
1158 
1159  logical, intent(in), optional :: flag_extrap ! when true, vertical extrapolation will be executed (just copy)
1160 
1161  integer :: nsize, psize, nidx_max
1162  integer, allocatable :: idx_blk(:,:,:), nidx(:,:)
1163  real(rp) :: lon_min, lon_max
1164  real(rp) :: lat_min, lat_max
1165  real(rp) :: dlon, dlat
1166  integer :: idx_ref(npoints)
1167  real(rp) :: hfactt(npoints)
1168 
1169  integer :: i, j, ii, jj, n
1170  !---------------------------------------------------------------------------
1171 
1172  call prof_rapstart('INTERP_fact',3)
1173 
1174  nsize = ia_ref * ja_ref
1175  if ( nsize > 100 ) then
1176  psize = int( sqrt(2.0_rp*sqrt(real(nsize,rp))) )
1177  nidx_max = nsize / psize * interp_buffer_size_fact
1178  else
1179  psize = 1
1180  nidx_max = nsize
1181  end if
1182 
1183  allocate(idx_blk(nidx_max,psize,psize))
1184  allocate(nidx( psize,psize))
1185 
1186  !$acc data copyin(lon_ref, lat_ref, hgt_ref, lon, lat, hgt) &
1187  !$acc copyout(idx_i, idx_j, hfact, idx_k, vfact) &
1188  !$acc create(idx_blk, nidx)
1189 
1190  call interp_div_block(nsize, psize, nidx_max, & ! [IN]
1191  lon_ref(:,:), lat_ref(:,:), & ! [IN]
1192  idx_blk(:,:,:), nidx(:,:), & ! [OUT]
1193  lon_min, lon_max, & ! [OUT]
1194  lat_min, lat_max, & ! [OUT]
1195  dlon, dlat ) ! [OUT]
1196 
1197  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
1198  !$omp private(ii,jj,idx_ref)
1199  !$acc kernels
1200  !$acc loop independent
1201  do j = 1, ja
1202  !$acc loop private(idx_ref,hfactt) independent
1203  do i = 1, ia
1204 
1205  ! main search
1206  call interp_search_horiz( npoints, & ! [IN]
1207  nsize, & ! [IN]
1208  lon_ref(:,:), lat_ref(:,:), & ! [IN]
1209  lon_min, lon_max, & ! [IN]
1210  lat_min, lat_max, & ! [IN]
1211  psize, nidx_max, & ! [IN]
1212  dlon, dlat, & ! [IN]
1213  idx_blk(:,:,:), nidx(:,:), & ! [IN]
1214  lon(i,j), lat(i,j), & ! [IN]
1215  idx_ref(:), hfactt(:) ) ! [OUT]
1216 
1217  !$acc loop seq
1218  do n = 1, npoints
1219  ii = mod(idx_ref(n) - 1, ia_ref) + 1
1220  jj = ( idx_ref(n) - 1 ) / ia_ref + 1
1221  idx_i(i,j,n) = ii
1222  idx_j(i,j,n) = jj
1223  hfact(i,j,n) = hfactt(n)
1224  call interp_factor1d( ka_ref, ks_ref, ke_ref, & ! [IN]
1225  ka, ks, ke, & ! [IN]
1226  hgt_ref(:,ii,jj), & ! [IN]
1227  hgt(:,i,j), & ! [IN]
1228  idx_k(:,:,i,j,n), & ! [OUT]
1229  vfact(: ,i,j,n), & ! [OUT]
1230  flag_extrap = flag_extrap ) ! [IN, optional]
1231 
1232  enddo
1233  enddo
1234  enddo
1235  !$acc end kernels
1236 
1237  !$acc end data
1238 
1239  deallocate(idx_blk, nidx)
1240 
1241  call prof_rapend ('INTERP_fact',3)
1242 
1243  return
1244  end subroutine interp_factor3d_weight
1245 
1246  !-----------------------------------------------------------------------------
1247  ! interpolation for 1D data
1248 !OCL SERIAL
1249  subroutine interp_interp1d( &
1250  KA_ref, KS_ref, KE_ref, &
1251  KA, KS, KE, &
1252 #ifdef _OPENACC
1253  workr, worki, &
1254 #endif
1255  idx_k, &
1256  vfact, &
1257  hgt_ref, &
1258  hgt, &
1259  val_ref, &
1260  val, &
1261  spline, &
1262  logwgt )
1263  !$acc routine vector
1264  use scale_const, only: &
1265  undef => const_undef, &
1266  eps => const_eps
1267  implicit none
1268  integer, intent(in) :: ka_ref, ks_ref, ke_ref ! number of z-direction (reference)
1269  integer, intent(in) :: ka, ks, ke ! number of z-direction (target)
1270 
1271  integer, intent(in) :: idx_k(ka,2) ! k-index in reference
1272  real(rp), intent(in) :: vfact(ka ) ! vertical interp factor
1273  real(rp), intent(in) :: hgt_ref(ka_ref) ! height (reference)
1274  real(rp), intent(in) :: hgt (ka) ! height (target)
1275  real(rp), intent(in), target :: val_ref(ka_ref) ! value (reference)
1276 
1277  real(rp), intent(out) :: val (ka) ! value (target)
1278 
1279  logical, intent(in), optional :: spline
1280  logical, intent(in), optional :: logwgt
1281 
1282  logical :: spline_
1283  logical :: logwgt_
1284 
1285  real(rp), pointer :: work(:)
1286 
1287 #ifdef _OPENACC
1288  real(rp), intent(out), target :: workr(ka_ref,7)
1289  integer, intent(out) :: worki(ka_ref,2)
1290 #define idx_i1(k) worki(k,1)
1291 #define idx_r_i1(k) worki(k,2)
1292 #define FDZ_i1(k) workr(k,1)
1293 #define U_i1(k) workr(k,2)
1294 #else
1295  integer :: idx_i1 (ka_ref)
1296  integer :: idx_r_i1(ka_ref)
1297  real(rp) :: fdz_i1 (ka_ref)
1298  real(rp) :: u_i1 (ka_ref)
1299 #endif
1300 
1301  integer :: kmax
1302  integer :: k
1303  !---------------------------------------------------------------------------
1304 
1305  spline_ = interp_use_spline_vert
1306  if ( present(spline) ) then
1307  spline_ = spline
1308  end if
1309 
1310  logwgt_ = .false.
1311  if ( present(logwgt) ) then
1312  logwgt_ = logwgt
1313  endif
1314 
1315  if ( logwgt_ ) then
1316 #ifdef _OPENACC
1317  work => workr(:,3)
1318 #else
1319  allocate( work(ka_ref) )
1320 #endif
1321  do k = ks_ref, ke_ref
1322  if ( abs( val_ref(k) - undef ) < eps ) then
1323  work(k) = undef
1324  else
1325  work(k) = log( val_ref(k) )
1326  end if
1327  enddo
1328  else
1329  work => val_ref
1330  endif
1331 
1332  call spline_coef( ka_ref, ks_ref, ke_ref, &
1333 #ifdef _OPENACC
1334  workr(:,4:7), &
1335 #endif
1336  hgt_ref(:), work(:), & ! (in)
1337  spline_, & ! (in)
1338  kmax, & ! (out)
1339  idx_i1(:), idx_r_i1(:), & ! (out)
1340  u_i1(:), fdz_i1(:) ) ! (out)
1341 
1342  call spline_exec( ka_ref, kmax, ka, ks, ke, &
1343  idx_k(:,:), vfact(:), & ! (in)
1344  hgt_ref(:), hgt(:), & ! (in)
1345  work(:), & ! (in)
1346  idx_i1(:), idx_r_i1(:), & ! (in)
1347  u_i1(:), fdz_i1(:), & ! (in)
1348  val(:) ) ! (out)
1349 
1350  if ( logwgt_ ) then
1351 #ifndef _OPENACC
1352  deallocate( work )
1353 #endif
1354  do k = ks, ke
1355  if ( abs( val(k) - undef ) > eps ) then
1356  val(k) = exp( val(k) )
1357  endif
1358  end do
1359  endif
1360 
1361  return
1362  end subroutine interp_interp1d
1363 
1364  !-----------------------------------------------------------------------------
1365  ! interpolation for 2D data (nearest-neighbor)
1366  subroutine interp_interp2d( &
1367  npoints, &
1368  IA_ref, JA_ref, &
1369  IA, JA, &
1370  idx_i, idx_j, &
1371  hfact, &
1372  val_ref, &
1373  val, &
1374  threshold_undef, &
1375  wsum, val2 )
1376  use scale_const, only: &
1377  undef => const_undef, &
1378  eps => const_eps
1379  implicit none
1380 
1381  integer, intent(in) :: npoints ! number of interpolation point for horizontal
1382  integer, intent(in) :: ia_ref ! number of x-direction (reference)
1383  integer, intent(in) :: ja_ref ! number of y-direction (reference)
1384  integer, intent(in) :: ia ! number of x-direction (target)
1385  integer, intent(in) :: ja ! number of y-direction (target)
1386  integer, intent(in) :: idx_i (ia,ja,npoints) ! i-index in reference
1387  integer, intent(in) :: idx_j (ia,ja,npoints) ! j-index in reference
1388  real(rp), intent(in) :: hfact (ia,ja,npoints) ! horizontal interp factor
1389  real(rp), intent(in) :: val_ref(ia_ref,ja_ref) ! value (reference)
1390 
1391  real(rp), intent(out) :: val (ia,ja) ! value (target)
1392 
1393  real(rp), intent(in), optional :: threshold_undef
1394  real(rp), intent(out), optional :: wsum(ia,ja)
1395  real(rp), intent(out), optional :: val2(ia,ja)
1396 
1397  real(rp) :: th_undef
1398 
1399  real(rp) :: fact, valn, f, w
1400  real(rp) :: sw
1401  logical :: lval2
1402 
1403  integer :: i, j, n
1404  !---------------------------------------------------------------------------
1405 
1406  call prof_rapstart('INTERP_interp',3)
1407 
1408  if ( present(threshold_undef) ) then
1409  th_undef = threshold_undef
1410  else
1411  th_undef = interp_threshold_undef
1412  end if
1413  th_undef = min( max( th_undef, eps * 2.0_rp ), 1.0_rp - eps * 2.0_rp )
1414 
1415  lval2 = present(wsum) .and. present(val2)
1416 
1417  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
1418  !$omp private(fact,valn,f,w,sw)
1419  !$acc kernels
1420 !OCL PREFETCH
1421  do j = 1, ja
1422  do i = 1, ia
1423  fact = 0.0_rp
1424  valn = 0.0_rp
1425  !$acc loop seq
1426  do n = 1, npoints
1427  f = hfact(i,j,n)
1428  w = val_ref(idx_i(i,j,n),idx_j(i,j,n))
1429  if ( f > eps .and. abs( w - undef ) > eps ) then
1430  fact = fact + f
1431  valn = valn + f * w
1432  end if
1433  end do
1434  sw = 0.5_rp - sign( 0.5_rp, fact - th_undef + eps ) ! 1.0 when fact < threshold
1435  val(i,j) = valn / ( fact + sw ) * ( 1.0_rp - sw ) + undef * sw
1436  if ( lval2 ) then
1437  wsum(i,j) = fact
1438  sw = 0.5_rp - sign( 0.5_rp, fact - eps ) ! 1.0 when fact < 0.0
1439  val2(i,j) = valn / ( fact + sw ) * ( 1.0_rp - sw ) + undef * sw
1440  end if
1441  enddo
1442  enddo
1443  !$acc end kernels
1444 
1445  call prof_rapend ('INTERP_interp',3)
1446 
1447  return
1448  end subroutine interp_interp2d
1449 
1450  !-----------------------------------------------------------------------------
1451  ! interpolation using one-points for 3D data (nearest-neighbor)
1452  subroutine interp_interp3d( &
1453  npoints, &
1454  KA_ref, KS_ref, KE_ref, &
1455  IA_ref, JA_ref, &
1456  KA, KS, KE, &
1457  IA, JA, &
1458  idx_i, idx_j, &
1459  hfact, &
1460  idx_k, &
1461  vfact, &
1462  hgt_ref, &
1463  hgt, &
1464  val_ref, &
1465  val, &
1466  spline, logwgt, &
1467  threshold_undef, &
1468  wsum, val2 )
1469  use scale_const, only: &
1470  undef => const_undef, &
1471  eps => const_eps
1472  implicit none
1473  integer, intent(in) :: npoints ! number of interpolation point for horizontal
1474  integer, intent(in) :: ka_ref, ks_ref, ke_ref ! number of z-direction (reference)
1475  integer, intent(in) :: ia_ref ! number of x-direction (reference)
1476  integer, intent(in) :: ja_ref ! number of y-direction (reference)
1477  integer, intent(in) :: ka, ks, ke ! number of z-direction (target)
1478  integer, intent(in) :: ia ! number of x-direction (target)
1479  integer, intent(in) :: ja ! number of y-direction (target)
1480 
1481  integer, intent(in) :: idx_i (ia,ja,npoints) ! i-index in reference
1482  integer, intent(in) :: idx_j (ia,ja,npoints) ! j-index in reference
1483  real(rp), intent(in) :: hfact (ia,ja,npoints) ! horizontal interp factor
1484  integer, intent(in) :: idx_k (ka,2,ia,ja,npoints) ! k-index in reference
1485  real(rp), intent(in) :: vfact (ka, ia,ja,npoints) ! vertical interp factor
1486  real(rp), intent(in) :: hgt_ref(ka_ref,ia_ref,ja_ref) ! height (reference)
1487  real(rp), intent(in) :: hgt (ka,ia,ja) ! height (target)
1488  real(rp), intent(in), target :: val_ref(ka_ref,ia_ref,ja_ref) ! value (reference)
1489 
1490  real(rp), intent(out) :: val (ka,ia,ja) ! value (target)
1491 
1492  logical, intent(in), optional :: spline
1493  logical, intent(in), optional :: logwgt
1494  real(rp), intent(in), optional :: threshold_undef
1495  real(rp), intent(out), optional :: wsum(ka,ia,ja)
1496  real(rp), intent(out), optional :: val2(ka,ia,ja)
1497 
1498  real(rp) :: th_undef
1499  logical :: spline_
1500  logical :: logwgt_
1501 
1502  real(rp), pointer :: work(:,:,:)
1503 
1504  integer, allocatable :: kmax(:,:)
1505  integer, allocatable :: idx(:,:,:), idx_r(:,:,:)
1506  real(rp), allocatable :: u(:,:,:), fdz(:,:,:)
1507 
1508  real(rp) :: valn
1509  real(rp) :: fact
1510  real(rp) :: f
1511  real(rp) :: sw
1512  real(rp) :: w(ka,npoints)
1513  logical :: lval2
1514 
1515 #ifdef _OPENACC
1516  real(rp) :: workr(ka_ref,4)
1517 #endif
1518 
1519  integer :: imin, imax
1520  integer :: jmin, jmax
1521  integer :: ii, jj
1522  integer :: k, i, j, n
1523  !---------------------------------------------------------------------------
1524 
1525  call prof_rapstart('INTERP_interp',3)
1526 
1527  if ( present(threshold_undef) ) then
1528  th_undef = threshold_undef
1529  else
1530  th_undef = interp_threshold_undef
1531  end if
1532  th_undef = min( max( th_undef, eps * 2.0_rp ), 1.0_rp - eps * 2.0_rp )
1533 
1534 
1535  spline_ = interp_use_spline_vert
1536  if ( present(spline) ) then
1537  spline_ = spline
1538  end if
1539 
1540  logwgt_ = .false.
1541  if ( present(logwgt) ) then
1542  logwgt_ = logwgt
1543  endif
1544 
1545 
1546  lval2 = present(wsum) .and. present(val2)
1547 
1548 
1549  !$acc data copyin(idx_i, idx_j, hfact, idx_k, vfact, hgt_ref, hgt, val_ref) &
1550  !$acc copyout(val)
1551 
1552  !$acc data copyout(wsum) if(present(wsum))
1553  !$acc data copyout(val2) if(present(val2))
1554 
1555  imin = ia_ref
1556  jmin = ja_ref
1557  imax = 1
1558  jmax = 1
1559  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
1560  !$omp reduction(min: imin,jmin) &
1561  !$omp reduction(max: imax,jmax)
1562  !$acc kernels
1563  !$acc loop reduction(min:imin,jmin) reduction(max:imax,jmax)
1564  do n = 1, npoints
1565  !$acc loop reduction(min:imin,jmin) reduction(max:imax,jmax)
1566  do j = 1, ja
1567  !$acc loop reduction(min:imin,jmin) reduction(max:imax,jmax)
1568  do i = 1, ia
1569  imin = min(imin, idx_i(i,j,n))
1570  imax = max(imax, idx_i(i,j,n))
1571  jmin = min(jmin, idx_j(i,j,n))
1572  jmax = max(jmax, idx_j(i,j,n))
1573  end do
1574  end do
1575  end do
1576  !$acc end kernels
1577 
1578  allocate( kmax(imin:imax,jmin:jmax) )
1579  allocate( idx(ka_ref,imin:imax,jmin:jmax), idx_r(ka_ref,imin:imax,jmin:jmax) )
1580  allocate( u(ka_ref,imin:imax,jmin:jmax), fdz(ka_ref,imin:imax,jmin:jmax) )
1581  !$acc enter data create(kmax, idx, idx_r, U, FDZ)
1582 
1583  if ( logwgt_ ) then
1584  allocate( work(ka_ref,imin:imax,jmin:jmax) )
1585  !$acc enter data create(work)
1586  !$omp parallel do OMP_SCHEDULE_ collapse(2)
1587  !$acc kernels
1588  !$acc loop independent
1589  do j = jmin, jmax
1590  !$acc loop independent
1591  do i = imin, imax
1592  !$acc loop independent
1593  do k = ks_ref, ke_ref
1594  if ( abs( val_ref(k,i,j) - undef ) < eps ) then
1595  work(k,i,j) = undef
1596  else
1597  work(k,i,j) = log( val_ref(k,i,j) )
1598  end if
1599  enddo
1600  enddo
1601  enddo
1602  !$acc end kernels
1603  else
1604  work => val_ref
1605  endif
1606 
1607  !$omp parallel do OMP_SCHEDULE_ collapse(2)
1608  !$acc kernels
1609  !$acc loop independent
1610  do j = jmin, jmax
1611  !$acc loop independent private(workr)
1612  do i = imin, imax
1613  call spline_coef( ka_ref, ks_ref, ke_ref, &
1614 #ifdef _OPENACC
1615  workr(:,:), &
1616 #endif
1617  hgt_ref(:,i,j), work(:,i,j), & ! (in)
1618  spline_, & ! (in)
1619  kmax(i,j), & ! (out)
1620  idx(:,i,j), idx_r(:,i,j), & ! (out)
1621  u(:,i,j), fdz(:,i,j) ) ! (out)
1622  end do
1623  end do
1624  !$acc end kernels
1625 
1626  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
1627  !$omp private(valn,fact,w,f,ii,jj,sw)
1628  !$acc kernels
1629  !$acc loop independent
1630  do j = 1, ja
1631  !$acc loop independent private(w)
1632  do i = 1, ia
1633 
1634  !$acc loop seq
1635  do n = 1, npoints
1636  if ( hfact(i,j,n) > eps ) then
1637  ii = idx_i(i,j,n)
1638  jj = idx_j(i,j,n)
1639  call spline_exec( ka_ref, kmax(ii,jj), ka, ks, ke, &
1640  idx_k(:,:,i,j,n), vfact(:, i,j,n), & ! (in)
1641  hgt_ref(:,ii,jj), hgt(:,i,j), & ! (in)
1642  work(:,ii,jj), & ! (in)
1643  idx(:,ii,jj), idx_r(:,ii,jj), & ! (in)
1644  u(:,ii,jj), fdz(:,ii,jj), & ! (in)
1645  w(:,n) ) ! (out)
1646  end if
1647  end do
1648 
1649  do k = ks, ke
1650  fact = 0.0_rp
1651  valn = 0.0_rp
1652  !$acc loop seq
1653  do n = 1, npoints
1654  f = hfact(i,j,n)
1655  if ( f > eps .and. abs( w(k,n) - undef ) > eps ) then
1656  fact = fact + f
1657  valn = valn + f * w(k,n)
1658  endif
1659  enddo
1660  sw = 0.5_rp - sign( 0.5_rp, fact - th_undef ) ! 1.0 when fact < threshold
1661  val(k,i,j) = valn / ( fact + sw ) * ( 1.0_rp - sw ) + undef * sw
1662  if ( lval2 ) then
1663  wsum(k,i,j) = fact
1664  sw = 0.5_rp - sign( 0.5_rp, fact - eps ) ! 1.0 when fact < 0.0
1665  val2(k,i,j) = valn / ( fact + sw ) * ( 1.0_rp - sw ) + undef * sw
1666  end if
1667  enddo
1668 
1669  enddo
1670  enddo
1671  !$acc end kernels
1672 
1673  !$acc exit data delete(kmax, idx, idx_r, U, FDZ)
1674  deallocate( kmax, idx, idx_r )
1675  deallocate( u, fdz )
1676 
1677  if ( logwgt_ ) then
1678  !$acc exit data delete(work)
1679  deallocate( work )
1680  !$omp parallel do OMP_SCHEDULE_ collapse(2)
1681  !$acc kernels
1682  do j = 1, ja
1683  do i = 1, ia
1684  do k = ks, ke
1685  if ( abs( val(k,i,j) - undef ) > eps ) then
1686  val(k,i,j) = exp( val(k,i,j) )
1687  end if
1688  end do
1689  end do
1690  end do
1691  !$acc end kernels
1692  end if
1693 
1694  !$acc end data
1695  !$acc end data
1696  !$acc end data
1697 
1698  call prof_rapend ('INTERP_interp',3)
1699 
1700  return
1701  end subroutine interp_interp3d
1702 
1703 
1704  ! private
1705 
1706 
1707  !-----------------------------------------------------------------------------
1708  ! horizontal search of interpolation points
1709 !OCL SERIAL
1710  subroutine interp_search_horiz( &
1711  npoints, &
1712  nsize, &
1713  lon_ref, lat_ref, &
1714  lon_min, lon_max, &
1715  lat_min, lat_max, &
1716  psize, nidx_max, &
1717  dlon, dlat, &
1718  idx_blk, nidx, &
1719  lon, lat, &
1720  idx_ref, &
1721  hfact, &
1722  search_limit, &
1723  weight_order )
1724  !$acc routine seq
1725  use scale_const, only: &
1726  eps => const_eps, &
1727  radius => const_radius
1728  implicit none
1729 
1730  integer, intent(in) :: npoints ! number of interpolation point for horizontal
1731  integer, intent(in) :: nsize ! number of grids (reference)
1732  real(rp), intent(in) :: lon_ref(nsize) ! longitude [rad] (reference)
1733  real(rp), intent(in) :: lat_ref(nsize) ! latitude [rad] (reference)
1734  real(rp), intent(in) :: lon_min ! minmum longitude [rad] (reference)
1735  real(rp), intent(in) :: lon_max ! maximum longitude [rad] (reference)
1736  real(rp), intent(in) :: lat_min ! minmum latitude [rad] (reference)
1737  real(rp), intent(in) :: lat_max ! maximum latitude [rad] (reference)
1738  integer, intent(in) :: psize ! number of blocks for each dimension
1739  integer, intent(in) :: nidx_max ! maximum number of index in the block
1740  real(rp), intent(in) :: dlon ! block longitude difference
1741  real(rp), intent(in) :: dlat ! block latitude difference
1742  integer, intent(in) :: idx_blk(nidx_max,psize,psize) ! index of the reference in the block
1743  integer, intent(in) :: nidx (psize,psize) ! number of indexes in the block
1744  real(rp), intent(in) :: lon ! longitude [rad] (target)
1745  real(rp), intent(in) :: lat ! latitude [rad] (target)
1746  integer, intent(out) :: idx_ref(npoints) ! index in reference (target)
1747  real(rp), intent(out) :: hfact(npoints) ! horizontal interp factor (target)
1748 
1749  real(rp), intent(in), optional :: search_limit
1750  integer, intent(in), optional :: weight_order
1751 
1752  real(rp) :: drad(npoints)
1753  real(rp) :: sum
1754  real(rp) :: search_limit_
1755  integer :: weight_order_
1756 
1757  real(rp) :: lon0, lon1, lat0, lat1
1758  real(rp) :: dlon_sl, dlat_sl
1759  integer :: i, n
1760  integer :: ii, jj, ii0, jj0
1761  !---------------------------------------------------------------------------
1762 
1763  if ( present(search_limit) ) then
1764  search_limit_ = search_limit
1765  else
1766  search_limit_ = interp_search_limit
1767  end if
1768  search_limit_ = search_limit_ / radius ! m to radian
1769 
1770  if ( present(weight_order) ) then
1771  weight_order_ = weight_order
1772  else
1773  weight_order_ = interp_weight_order
1774  end if
1775 
1776 
1777 
1778  drad(:) = large_number
1779  idx_ref(:) = -1
1780 
1781  ! find k-nearest points in the nearest block
1782  ii0 = max( min( int(( lon - lon_min ) / dlon) + 1, psize ), 1 )
1783  jj0 = max( min( int(( lat - lat_min ) / dlat) + 1, psize ), 1 )
1784  do i = 1, nidx(ii0,jj0)
1785  n = idx_blk(i,ii0,jj0)
1786  call interp_insert( npoints, &
1787  lon, lat, &
1788  lon_ref(n), lat_ref(n), & ! [IN]
1789  n, & ! [IN]
1790  drad(:), idx_ref(:) ) ! [INOUT]
1791  end do
1792 
1793  if ( abs(drad(1)) < eps ) then
1794  hfact(:) = 0.0_rp
1795  hfact(1) = 1.0_rp
1796 
1797  return
1798  else if ( drad(1) > search_limit_ ) then
1799  hfact(:) = 0.0_rp
1800  idx_ref(:) = 1 ! dummy
1801 
1802  return
1803  end if
1804 
1805  dlon_sl = max(dlon * 0.5_rp, drad(npoints))
1806  dlat_sl = max(dlat * 0.5_rp, drad(npoints))
1807  do jj = 1, psize
1808  lat0 = lat_min + dlat * (jj-1)
1809  lat1 = lat_min + dlat * jj
1810  if ( lat < lat0 - dlat_sl &
1811  .or. lat >= lat1 + dlat_sl ) cycle
1812  do ii = 1, psize
1813  if ( ii==ii0 .and. jj==jj0 ) cycle
1814  lon0 = lon_min + dlon * (ii-1)
1815  lon1 = lon_min + dlon * ii
1816  if ( lon < lon0 - dlon_sl &
1817  .or. lon >= lon1 + dlon_sl ) cycle
1818  do i = 1, nidx(ii,jj)
1819  n = idx_blk(i,ii,jj)
1820  call interp_insert( npoints, &
1821  lon, lat, &
1822  lon_ref(n), lat_ref(n), & ! [IN]
1823  n, & ! [IN]
1824  drad(:), idx_ref(:) ) ! [INOUT]
1825  end do
1826  dlon_sl = max(dlon * 0.5_rp, drad(npoints))
1827  dlat_sl = max(dlat * 0.5_rp, drad(npoints))
1828  end do
1829  end do
1830 
1831  ! factor = 1 / dradian
1832  if ( weight_order_ < 0 ) then
1833  do n = 1, npoints
1834  hfact(n) = 1.0_rp / drad(n)**(-1.0_rp/weight_order_)
1835  enddo
1836  else
1837  do n = 1, npoints
1838  hfact(n) = 1.0_rp / drad(n)**weight_order_
1839  enddo
1840  end if
1841 
1842  ! ignore far point
1843  do n = 1, npoints
1844  if ( drad(n) >= search_limit_ ) then
1845  hfact(n) = 0.0_rp
1846  endif
1847  enddo
1848 
1849  ! normalize factor
1850  sum = 0.0_rp
1851  do n = 1, npoints
1852  sum = sum + hfact(n)
1853  enddo
1854 
1855  if ( sum > 0.0_rp ) then
1856  do n = 1, npoints
1857  hfact(n) = hfact(n) / sum
1858  enddo
1859  endif
1860 
1861  return
1862  end subroutine interp_search_horiz
1863 
1864  !-----------------------------------------------------------------------------
1865  ! horizontal search of interpolation points for structure grid
1866 !OCL SERIAL
1867  subroutine interp_search_horiz_struct( &
1868  npoints, &
1869  psizex, psizey, &
1870  IA_ref, JA_ref, &
1871  lon_ref, lat_ref, &
1872  lon_min, lat_min, &
1873  dlon, dlat, &
1874  i0, i1, j0, j1, &
1875  lon, lat, &
1876  idx_i, idx_j, &
1877  hfact, &
1878  search_limit, &
1879  weight_order )
1880  !$acc routine seq
1881  use scale_const, only: &
1882  eps => const_eps, &
1883  radius => const_radius
1884  implicit none
1885 
1886  integer, intent(in) :: npoints ! number of interpolation point for horizontal
1887  integer, intent(in) :: psizex ! number of block in x-direction (reference)
1888  integer, intent(in) :: psizey ! number of block in y-direction (reference)
1889  integer, intent(in) :: IA_ref ! number of grids in x-direction (reference)
1890  integer, intent(in) :: JA_ref ! number of grids in y-direction (reference)
1891  real(RP), intent(in) :: lon_ref(IA_ref) ! longitude [rad] (reference)
1892  real(RP), intent(in) :: lat_ref(JA_ref) ! latitude [rad] (reference)
1893  real(RP), intent(in) :: lon_min ! minimum longitude
1894  real(RP), intent(in) :: lat_min ! minimum latitude
1895  real(RP), intent(in) :: dlon
1896  real(RP), intent(in) :: dlat
1897  integer, intent(in) :: i0(psizex) ! start index in the block
1898  integer, intent(in) :: i1(psizex) ! end index in the block
1899  integer, intent(in) :: j0(psizey) ! start index in the block
1900  integer, intent(in) :: j1(psizey) ! start index in the block
1901  real(RP), intent(in) :: lon ! longitude [rad] (target)
1902  real(RP), intent(in) :: lat ! latitude [rad] (target)
1903  integer, intent(out) :: idx_i(npoints) ! x-index in reference (target)
1904  integer, intent(out) :: idx_j(npoints) ! y-index in reference (target)
1905  real(RP), intent(out) :: hfact(npoints) ! horizontal interp factor (target)
1906 
1907  real(RP), intent(in), optional :: search_limit
1908  integer, intent(in), optional :: weight_order
1909 
1910  real(RP) :: drad(npoints)
1911  real(RP) :: sum
1912  real(RP) :: search_limit_
1913  integer :: weight_order_
1914 
1915  real(RP) :: lon0, lon1, lat0, lat1
1916  real(RP) :: dlon_sl, dlat_sl
1917 
1918  integer :: i, j, n
1919  integer :: ii, jj, ii0, jj0
1920  !---------------------------------------------------------------------------
1921 
1922  if ( present(search_limit) ) then
1923  search_limit_ = search_limit
1924  else
1925  search_limit_ = interp_search_limit
1926  end if
1927  search_limit_ = search_limit_ / radius ! m to radian
1928 
1929  if ( present(weight_order) ) then
1930  weight_order_ = weight_order
1931  else
1932  weight_order_ = interp_weight_order
1933  end if
1934 
1935 
1936  drad(:) = large_number
1937  idx_i(:) = 1
1938  idx_j(:) = 1
1939 
1940  ! find k-nearest points in the nearest block
1941  ii0 = max( min( int(( lon - lon_min ) / dlon) + 1, psizex ), 1)
1942  jj0 = max( min( int(( lat - lat_min ) / dlat) + 1, psizey ), 1)
1943  do j = j0(jj0), j1(jj0)
1944  do i = i0(ii0), i1(ii0)
1945  call interp_insert_2d( npoints, &
1946  lon, lat, &
1947  lon_ref(i), lat_ref(j), & ! [IN]
1948  i, j, & ! [IN]
1949  drad(:), idx_i(:), idx_j(:) ) ! [INOUT]
1950  end do
1951  end do
1952 
1953  if ( abs(drad(1)) < eps ) then
1954  hfact(:) = 0.0_rp
1955  hfact(1) = 1.0_rp
1956 
1957  return
1958  else if ( drad(1) > search_limit_ ) then
1959  hfact(:) = 0.0_rp
1960  idx_i(:) = 1 ! dummy
1961  idx_j(:) = 1 ! dummy
1962 
1963  return
1964  end if
1965 
1966  dlon_sl = max(dlon * 0.5_rp, drad(npoints))
1967  dlat_sl = max(dlat * 0.5_rp, drad(npoints))
1968  do jj = 1, psizey
1969  lat0 = lat_min + dlat * (jj-1)
1970  lat1 = lat_min + dlat * jj
1971  if ( lat < lat0 - dlat_sl &
1972  .or. lat >= lat1 + dlat_sl ) cycle
1973  do ii = 1, psizex
1974  if ( ii==ii0 .and. jj==jj0 ) cycle
1975  lon0 = lon_min + dlon * (ii-1)
1976  lon1 = lon_min + dlon * ii
1977  if ( lon < lon0 - dlon_sl &
1978  .or. lon >= lon1 + dlon_sl ) cycle
1979  do j = j0(jj0), j1(jj0)
1980  do i = i0(ii0), i1(ii0)
1981  call interp_insert_2d( npoints, &
1982  lon, lat, &
1983  lon_ref(i), lat_ref(j), & ! [IN]
1984  i, j, & ! [IN]
1985  drad(:), idx_i(:), idx_j(:) ) ! [INOUT]
1986  end do
1987  end do
1988  dlon_sl = max(dlon * 0.5_rp, drad(npoints))
1989  dlat_sl = max(dlat * 0.5_rp, drad(npoints))
1990  end do
1991  end do
1992 
1993 
1994  ! factor = 1 / drad
1995  if ( weight_order_ < 0 ) then
1996  do n = 1, npoints
1997  hfact(n) = 1.0_rp / drad(n)**(-1.0_rp/weight_order_)
1998  enddo
1999  else
2000  do n = 1, npoints
2001  hfact(n) = 1.0_rp / drad(n)**weight_order_
2002  enddo
2003  end if
2004 
2005  ! ignore far point
2006  do n = 1, npoints
2007  if ( drad(n) >= search_limit_ ) then
2008  hfact(n) = 0.0_rp
2009  endif
2010  enddo
2011 
2012  ! normalize factor
2013  sum = 0.0_rp
2014  do n = 1, npoints
2015  sum = sum + hfact(n)
2016  enddo
2017 
2018  if ( sum > 0.0_rp ) then
2019  do n = 1, npoints
2020  hfact(n) = hfact(n) / sum
2021  enddo
2022  endif
2023 
2024  return
2025  end subroutine interp_search_horiz_struct
2026 
2027 
2028 !OCL SERIAL
2029  subroutine interp_insert( npoints, &
2030  lon, lat, &
2031  lon_ref, lat_ref, &
2032  i, &
2033  drad, idx_i )
2034  !$acc routine seq
2035  use scale_const, only: &
2036  undef => const_undef, &
2037  eps => const_eps
2038  use scale_sort, only: &
2039  sort_exec
2040 
2041  integer, intent(in) :: npoints
2042  real(RP), intent(in) :: lon, lat
2043  real(RP), intent(in) :: lon_ref, lat_ref
2044  integer, intent(in), value :: i
2045 
2046  real(RP), intent(inout) :: drad(npoints)
2047  integer, intent(inout) :: idx_i(npoints)
2048 
2049  real(RP) :: dradian
2050 
2051  if ( abs( lon_ref - undef ) < eps ) return
2052 
2053  dradian = haversine( lon, lat, lon_ref, lat_ref )
2054 
2055  if ( dradian <= drad(npoints) ) then
2056  ! replace last(=longest) value
2057  drad(npoints) = dradian
2058  idx_i(npoints) = i
2059 
2060  ! sort by ascending order
2061  call sort_exec( npoints, & ! [IN]
2062  drad(:), & ! [INOUT]
2063  idx_i(:) ) ! [INOUT]
2064  endif
2065 
2066  return
2067  end subroutine interp_insert
2068 
2069 !OCL SERIAL
2070  subroutine interp_insert_2d( npoints, &
2071  lon, lat, &
2072  lon_ref, lat_ref, &
2073  i, j, &
2074  drad, idx_i, idx_j )
2075  !$acc routine seq
2076  use scale_const, only: &
2077  undef => const_undef, &
2078  eps => const_eps
2079  use scale_sort, only: &
2080  sort_exec
2081 
2082  integer, intent(in) :: npoints
2083  real(RP), intent(in) :: lon, lat
2084  real(RP), intent(in) :: lon_ref, lat_ref
2085  integer, intent(in), value :: i, j
2086 
2087  real(RP), intent(inout) :: drad(npoints)
2088  integer, intent(inout) :: idx_i(npoints)
2089  integer, intent(inout) :: idx_j(npoints)
2090 
2091  real(RP) :: dradian
2092 
2093  if ( abs( lon_ref - undef ) < eps ) return
2094 
2095  dradian = haversine( lon, lat, lon_ref, lat_ref )
2096 
2097  if ( dradian <= drad(npoints) ) then
2098  ! replace last(=longest) value
2099  drad(npoints) = dradian
2100  idx_i(npoints) = i
2101  idx_j(npoints) = j
2102 
2103  ! sort by ascending order
2104  call sort_exec( npoints, & ! [IN]
2105  drad(:), & ! [INOUT]
2106  idx_i(:), & ! [INOUT]
2107  idx_j(:) ) ! [INOUT]
2108  endif
2109 
2110  return
2111  end subroutine interp_insert_2d
2112 
2113 
2114  subroutine interp_div_block(nsize, psize, nidx_max, &
2115  lon_ref, lat_ref, &
2116  idx, nidx, &
2117  lon_min, lon_max, &
2118  lat_min, lat_max, &
2119  dlon, dlat )
2120  use scale_const, only: &
2121  undef => const_undef, &
2122  eps => const_eps
2123  use scale_prc, only: &
2124  prc_abort
2125  integer, intent(in) :: nsize
2126  integer, intent(in) :: psize
2127  integer, intent(in) :: nidx_max
2128  real(RP), intent(in) :: lon_ref(nsize)
2129  real(RP), intent(in) :: lat_ref(nsize)
2130 
2131  integer, intent(out) :: idx (nidx_max,psize,psize)
2132  integer, intent(out) :: nidx(psize,psize)
2133  real(RP), intent(out) :: lon_min, lon_max
2134  real(RP), intent(out) :: lat_min, lat_max
2135  real(RP), intent(out) :: dlon, dlat
2136 
2137  integer :: i, ii, jj, n
2138 
2139  !$acc data copyin(lon_ref, lat_ref) copyout(idx, nidx)
2140 
2141  lon_min = 999.d0
2142  lon_max = -999.d0
2143  lat_min = 999.d0
2144  lat_max = -999.d0
2145  !$omp parallel do reduction(min: lon_min, lat_min) reduction(max: lon_max, lat_max)
2146  !$acc kernels
2147  !$acc loop reduction(min: lon_min, lat_min) reduction(max: lon_max, lat_max)
2148  do i = 1, nsize
2149  if ( abs(lon_ref(i) - undef) > eps ) then
2150  lon_min = min( lon_min, lon_ref(i) )
2151  lon_max = max( lon_max, lon_ref(i) )
2152  end if
2153  if ( abs(lat_ref(i) - undef) > eps ) then
2154  lat_min = min( lat_min, lat_ref(i) )
2155  lat_max = max( lat_max, lat_ref(i) )
2156  end if
2157  end do
2158  !$acc end kernels
2159 
2160  dlon = ( lon_max - lon_min ) / psize
2161  dlat = ( lat_max - lat_min ) / psize
2162 
2163  !$acc kernels
2164  nidx(:,:) = 0
2165  !$acc end kernels
2166  !$acc kernels
2167  !$acc loop independent
2168  do i = 1, nsize
2169  if ( abs( lon_ref(i) - undef ) < eps ) cycle
2170  ii = min(int((lon_ref(i) - lon_min) / dlon) + 1, psize)
2171  jj = min(int((lat_ref(i) - lat_min) / dlat) + 1, psize)
2172  !$acc atomic capture
2173  n = nidx(ii,jj)
2174  nidx(ii,jj) = n + 1
2175  !$acc end atomic
2176  n = n + 1
2177  if ( n <= nidx_max ) idx(n,ii,jj) = i
2178  end do
2179  !$acc end kernels
2180 
2181  !$acc end data
2182 
2183  if ( maxval(nidx) > nidx_max ) then
2184  log_error("INTERP_search_horiz",*) 'Buffer size is not enough'
2185  log_error_cont(*) ' Use larger INTERP_buffer_size_fact: ', interp_buffer_size_fact * maxval(nidx)/nidx_max
2186  call prc_abort
2187  end if
2188 
2189  return
2190  end subroutine interp_div_block
2191 
2192  !-----------------------------------------------------------------------------
2193  subroutine interp_bilinear_inv( &
2194  x_ref0, x_ref1, x_ref2, x_ref3, &
2195  y_ref0, y_ref1, y_ref2, y_ref3, &
2196  x, y, &
2197  u, v, &
2198  error )
2199  !$acc routine seq
2200  implicit none
2201  real(RP), intent(in), value :: x_ref0, x_ref1, x_ref2, x_ref3
2202  real(RP), intent(in), value :: y_ref0, y_ref1, y_ref2, y_ref3
2203  real(RP), intent(in), value :: x, y
2204 
2205  real(RP), intent(out) :: u, v
2206  logical, intent(out) :: error
2207 
2208  real(RP) :: e_x, e_y
2209  real(RP) :: f_x, f_y
2210  real(RP) :: g_x, g_y
2211  real(RP) :: h_x, h_y
2212  real(RP) :: k0, k1, k2
2213  real(RP) :: w
2214  real(RP) :: sig
2215 
2216  e_x = x_ref1 - x_ref0
2217  e_y = y_ref1 - y_ref0
2218 
2219  f_x = x_ref3 - x_ref0
2220  f_y = y_ref3 - y_ref0
2221 
2222  g_x = x_ref0 - x_ref1 + x_ref2 - x_ref3
2223  g_y = y_ref0 - y_ref1 + y_ref2 - y_ref3
2224 
2225  h_x = x - x_ref0
2226  h_y = y - y_ref0
2227 
2228  k0 = cross(h_x, h_y, e_x, e_y)
2229  k1 = cross(e_x, e_y, f_x, f_y) + cross(h_x, h_y, g_x, g_y)
2230  k2 = cross(g_x, g_y, f_x, f_y)
2231 
2232  w = k1**2 - 4.0_rp * k0 * k2
2233  if ( w < 0.0_rp ) then
2234 #ifndef _OPENACC
2235  log_error("INTERP_bilinear_inv",*) 'Outside of the quadrilateral'
2236 #endif
2237  error = .true.
2238  return
2239  end if
2240 
2241  if ( abs(k1) < eps_bilinear ) then
2242 #ifndef _OPENACC
2243  log_error("INTERP_bilinear_inv",*) 'Unexpected error occured', k1
2244  log_error_cont(*) x_ref0, x_ref1, x_ref2, x_ref3
2245  log_error_cont(*) y_ref0, y_ref1, y_ref2, y_ref3
2246  log_error_cont(*) x, y
2247 #endif
2248  error = .true.
2249  return
2250  end if
2251  if ( abs(k2) < eps_bilinear * sqrt( (x_ref2-x_ref0)**2+(y_ref2-y_ref0)**2 ) ) then
2252  v = - k0 / k1
2253  else
2254  sig = sign( 1.0_rp, cross(x_ref1-x_ref0, y_ref1-y_ref0, x_ref3-x_ref0, y_ref3-y_ref0) )
2255  v = ( - k1 + sig * sqrt(w) ) / ( 2.0_rp * k2 )
2256  end if
2257  u = ( h_x - f_x * v ) / ( e_x + g_x * v )
2258 
2259  if ( u < -eps_bilinear .or. u > 1.0_rp+eps_bilinear .or. v < -eps_bilinear .or. v > 1.0_rp+eps_bilinear ) then
2260 #ifndef _OPENACC
2261  log_error("INTERP_bilinear_inv",*) 'Unexpected error occured', u, v
2262  log_error_cont(*) x_ref0, x_ref1, x_ref2, x_ref3
2263  log_error_cont(*) y_ref0, y_ref1, y_ref2, y_ref3
2264  log_error_cont(*) x, y
2265  log_error_cont(*) k0, k1, k2, sig
2266 #endif
2267  error = .true.
2268  return
2269  end if
2270  u = min( max( u, 0.0_rp ), 1.0_rp )
2271  v = min( max( v, 0.0_rp ), 1.0_rp )
2272 
2273  error = .false.
2274 
2275  return
2276  end subroutine interp_bilinear_inv
2277 
2278  !-----------------------------------------------------------------------------
2279  ! check whether the point is inside of the quadrilateral
2280  subroutine interp_check_inside( &
2281  x_ref0, x_ref1, x_ref2, x_ref3, &
2282  y_ref0, y_ref1, y_ref2, y_ref3, &
2283  x, y, &
2284  inc_i, inc_j, &
2285  error )
2286  !$acc routine seq
2287  use scale_const, only: &
2288  eps => const_eps
2289  implicit none
2290  real(RP), intent(in), value :: x_ref0, x_ref1, x_ref2, x_ref3
2291  real(RP), intent(in), value :: y_ref0, y_ref1, y_ref2, y_ref3
2292  real(RP), intent(in), value :: x, y
2293 
2294  integer, intent(out) :: inc_i, inc_j
2295  logical, intent(out) :: error
2296 
2297  real(RP) :: sig
2298  real(RP) :: c1, c2, c3, c4
2299  real(RP) :: th
2300  logical :: fx, fy
2301 
2302  error = .false.
2303 
2304  sig = sign( 1.0_rp, cross(x_ref1-x_ref0, y_ref1-y_ref0, x_ref3-x_ref0, y_ref3-y_ref0) )
2305 
2306  c1 = sig * cross(x_ref1-x_ref0, y_ref1-y_ref0, x-x_ref0, y-y_ref0)
2307  c2 = sig * cross(x_ref2-x_ref1, y_ref2-y_ref1, x-x_ref1, y-y_ref1)
2308  c3 = sig * cross(x_ref3-x_ref2, y_ref3-y_ref2, x-x_ref2, y-y_ref2)
2309  c4 = sig * cross(x_ref0-x_ref3, y_ref0-y_ref3, x-x_ref3, y-y_ref3)
2310 
2311  th = - max( abs(c1), abs(c2), abs(c3), abs(c4) ) * eps * 1.e4_rp
2312 
2313  ! if all the c1 - c4 are positive, the point is inside the quadrilateral
2314  inc_i = 0
2315  inc_j = 0
2316  if ( c1 < th ) inc_j = -1
2317  if ( c2 < th ) inc_i = 1
2318  if ( c3 < th ) inc_j = 1
2319  if ( c4 < th ) inc_i = -1
2320 
2321  fx = c2 < th .and. c4 < th
2322  fy = c1 < th .and. c3 < th
2323  if ( fx .and. fy ) then
2324 #ifndef _OPENACC
2325  log_error("INTERP_check_inside",*) 'Unexpected error occured', c1, c2, c3, c4
2326  log_error_cont(*) x_ref0, x_ref1, x_ref2, x_ref3
2327  log_error_cont(*) y_ref0, y_ref1, y_ref2, y_ref3
2328  log_error_cont(*) x, y
2329 #endif
2330  error = .true.
2331  return
2332  else if ( fx ) then
2333  inc_i = 0
2334  else if ( fy ) then
2335  inc_j = 0
2336  end if
2337 
2338  return
2339  end subroutine interp_check_inside
2340 
2341  !-----------------------------------------------------------------------------
2343  function cross( &
2344  x0, y0, &
2345  x1, y1 )
2346  !$acc routine seq
2347  implicit none
2348  real(rp), intent(in) :: x0, y0
2349  real(rp), intent(in) :: x1, y1
2350  real(rp) :: cross
2351 
2352  cross = x0 * y1 - x1 * y0
2353 
2354  end function cross
2355 
2356  !-----------------------------------------------------------------------------
2357  ! Haversine Formula (from R.W. Sinnott, "Virtues of the Haversine",
2358  ! Sky and Telescope, vol. 68, no. 2, 1984, p. 159):
2359  function haversine( &
2360  lon0, lat0, &
2361  lon1, lat1 ) &
2362  result( d )
2363  !$acc routine seq
2364  implicit none
2365 
2366  real(rp), intent(in) :: lon0 ! [rad]
2367  real(rp), intent(in) :: lat0 ! [rad]
2368  real(rp), intent(in) :: lon1 ! [rad]
2369  real(rp), intent(in) :: lat1 ! [rad]
2370  real(rp) :: d ! [rad]
2371 
2372  real(rp) :: dlonh, dlath
2373  real(rp) :: work1
2374  !---------------------------------------------------------------------------
2375 
2376  dlonh = 0.5_rp * ( lon0 - lon1 )
2377  dlath = 0.5_rp * ( lat0 - lat1 )
2378 
2379  work1 = sin(dlath)**2 + cos(lat0) * cos(lat1) * sin(dlonh)**2
2380 
2381  d = 2.0_rp * asin( min(sqrt(work1),1.0_rp) )
2382 
2383  end function haversine
2384 
2385 
2386 !OCL SERIAL
2387  subroutine spline_coef( &
2388  KA_ref, KS_ref, KE_ref, &
2389 #ifdef _OPENACC
2390  work, &
2391 #endif
2392  hgt_ref, val_ref, &
2393  spline, &
2394  kmax, &
2395  idx, idx_r, &
2396  U, FDZ )
2397  !$acc routine seq
2398  use scale_const, only: &
2399  undef => const_undef, &
2400  eps => const_eps
2401  use scale_matrix, only: &
2402  matrix_solver_tridiagonal_1d_ta
2403  implicit none
2404  integer, intent(in) :: ka_ref, ks_ref, ke_ref
2405 
2406  real(rp), intent(in) :: hgt_ref(ka_ref)
2407  real(rp), intent(in) :: val_ref(ka_ref)
2408  logical, intent(in) :: spline
2409 
2410  integer, intent(out) :: kmax
2411  integer, intent(out) :: idx (ka_ref)
2412  integer, intent(out) :: idx_r(ka_ref)
2413  real(rp), intent(out) :: u (ka_ref)
2414  real(rp), intent(out) :: fdz(ka_ref)
2415 
2416 #ifdef _OPENACC
2417  real(rp), intent(out) :: work(ka_ref,4)
2418 #define MD_sc(k) work(k,1)
2419 #define V_sc(k) work(k,2)
2420 #else
2421  real(rp) :: md_sc(ka_ref)
2422  real(rp) :: v_sc (ka_ref)
2423 #endif
2424  real(rp) :: dz
2425  integer :: k
2426 
2427  if ( spline ) then
2428 
2429  idx(1) = -999
2430  kmax = 1
2431  do k = ks_ref, ke_ref-1
2432  if ( hgt_ref(k) > undef .and. abs( val_ref(k) - undef ) > eps ) then
2433  idx(1) = k
2434  idx_r(k) = 1
2435  exit
2436  end if
2437  end do
2438  if ( idx(1) == -999 ) return ! UNDEF (use linear interpolation)
2439  fdz(1) = 1e10 ! dummy
2440  do k = idx(1)+1, ke_ref
2441  if ( hgt_ref(k) > undef ) then
2442  dz = hgt_ref(k) - hgt_ref(idx(kmax))
2443  if ( abs( val_ref(k) - undef ) > eps .and. dz > eps ) then
2444  do while ( kmax > 1 .and. fdz(kmax) < dz * 0.1_rp )
2445  kmax = kmax - 1 ! marge
2446  end do
2447  kmax = kmax + 1
2448  idx(kmax) = k
2449  if ( idx(kmax-1)+1 <= k-1 ) idx_r(idx(kmax-1)+1:k-1) = kmax-1
2450  idx_r(k) = kmax
2451  fdz(kmax) = hgt_ref(k) - hgt_ref(idx(kmax-1))
2452  end if
2453  end if
2454  end do
2455 
2456  if ( kmax > 3 ) then
2457 
2458  md_sc(2) = 2.0_rp * ( fdz(2) + fdz(3) ) + fdz(2)
2459  do k = 3, kmax-2
2460  md_sc(k) = 2.0_rp * ( fdz(k) + fdz(k+1) )
2461  end do
2462  md_sc(kmax-1) = 2.0_rp * ( fdz(kmax-1) + fdz(kmax) ) + fdz(kmax)
2463 
2464  do k = 2, kmax-1
2465  v_sc(k) = ( val_ref(idx(k+1)) - val_ref(idx(k )) ) / fdz(k+1) &
2466  - ( val_ref(idx(k )) - val_ref(idx(k-1)) ) / fdz(k )
2467  end do
2468 
2469  call matrix_solver_tridiagonal_1d_ta( kmax, 2, kmax-1, &
2470 #ifdef _OPENACC
2471  work(:,3:4), &
2472 #endif
2473  fdz(2:), md_sc(:), fdz(:), & ! (in)
2474  v_sc(:), & ! (in)
2475  u(:) ) ! (out)
2476 ! U(1) = 0.0_RP
2477 ! U(kmax) = 0.0_RP
2478  u(1) = u(2)
2479  u(kmax) = u(kmax-1)
2480 
2481  else
2482 
2483  idx(kmax) = idx(1) ! force linear interpolateion
2484 
2485  end if
2486 
2487  else
2488 
2489  kmax = 1
2490  idx(1) = -999 ! dummy
2491 
2492  end if
2493 
2494  return
2495  end subroutine spline_coef
2496 
2497 !OCL SERIAL
2498  subroutine spline_exec( &
2499  KA_ref, kmax, KA, KS, KE, &
2500  idx_k, vfact, &
2501  hgt_ref, hgt, &
2502  val_ref, &
2503  idx, idx_r, &
2504  U, FDZ, &
2505  val )
2506  !$acc routine vector
2507  use scale_const, only: &
2508  undef => const_undef, &
2509  eps => const_eps
2510  implicit none
2511  integer, intent(in) :: KA_ref, kmax
2512  integer, intent(in) :: KA, KS, KE
2513  integer, intent(in) :: idx_k(KA,2)
2514  real(RP), intent(in) :: vfact(KA)
2515  real(RP), intent(in) :: hgt_ref(KA_ref)
2516  real(RP), intent(in) :: hgt (KA)
2517  real(RP), intent(in) :: val_ref(KA_ref)
2518  integer, intent(in) :: idx (KA_ref)
2519  integer, intent(in) :: idx_r (KA_ref)
2520  real(RP), intent(in) :: U (KA_ref)
2521  real(RP), intent(in) :: FDZ (KA_ref)
2522 
2523  real(RP), intent(out) :: val(KA)
2524 
2525  real(RP) :: c1, c2, c3, d
2526  integer :: k, kk, kk2
2527 
2528  do k = ks, ke
2529  kk = idx_k(k,1)
2530  if ( kk == -1 ) then
2531  val(k) = undef
2532  else if ( idx_k(k,2) == -1 ) then
2533  val(k) = val_ref(kk)
2534  else if ( kk < idx(1) .or. kk >= idx(kmax) ) then ! linear interpolation
2535  if ( abs( val_ref(kk) - undef ) < eps .or. abs( val_ref(kk+1) - undef ) < eps ) then
2536  val(k) = undef
2537  else
2538  val(k) = val_ref(kk) * vfact(k) + val_ref(kk+1) * ( 1.0_rp - vfact(k) )
2539  end if
2540  else
2541  kk2 = idx_r(kk)
2542  kk = idx(kk2)
2543  c3 = ( u(kk2+1) - u(kk2) ) / fdz(kk2+1)
2544  c2 = 3.0_rp * u(kk2)
2545  c1 = ( val_ref(idx(kk2+1)) - val_ref(kk) ) / fdz(kk2+1) - ( 2.0_rp * u(kk2) + u(kk2+1) ) * fdz(kk2+1)
2546  d = hgt(k) - hgt_ref(kk)
2547 
2548  val(k) = ( ( c3 * d + c2 ) * d + c1 ) * d + val_ref(kk)
2549  end if
2550  end do
2551 
2552  return
2553  end subroutine spline_exec
2554 
2555 end module scale_interp
scale_interp::interp_interp2d
subroutine, public interp_interp2d(npoints, IA_ref, JA_ref, IA, JA, idx_i, idx_j, hfact, val_ref, val, threshold_undef, wsum, val2)
Definition: scale_interp.F90:1376
scale_interp::spline_exec
subroutine spline_exec(KA_ref, kmax, KA, KS, KE, idx_k, vfact, hgt_ref, hgt, val_ref, idx, idx_r, U, FDZ, val)
Definition: scale_interp.F90:2506
scale_matrix::matrix_solver_tridiagonal_1d_ta
subroutine, public matrix_solver_tridiagonal_1d_ta( KA, KS, KE, ifdef _OPENACC
solve tridiagonal matrix with Thomas's algorithm
Definition: scale_matrix.F90:133
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:350
scale_interp::interp_factor2d_weight
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)
Definition: scale_interp.F90:748
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_interp::interp_factor3d_linear_xy
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)
Definition: scale_interp.F90:1047
scale_sort
module SORT
Definition: scale_sort.F90:11
scale_interp
module INTERPOLATION
Definition: scale_interp.F90:12
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:35
scale_prof::prof_rapstart
subroutine, public prof_rapstart(rapname_base, level, disable_barrier)
Start raptime.
Definition: scale_prof.F90:174
scale_interp::interp_factor1d
subroutine, public interp_factor1d(KA_ref, KS_ref, KE_ref, KA, KS, KE, hgt_ref, hgt, idx_k, vfact, flag_extrap)
Definition: scale_interp.F90:256
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_io
module STDIO
Definition: scale_io.F90:10
scale_interp::interp_domain_compatibility
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)
Definition: scale_interp.F90:150
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_interp::interp_insert_2d
subroutine interp_insert_2d(npoints, lon, lat, lon_ref, lat_ref, i, j, drad, idx_i, idx_j)
Definition: scale_interp.F90:2075
scale_interp::interp_interp1d
subroutine, public interp_interp1d(KA_ref, KS_ref, KE_ref, KA, KS, KE, idx_k, vfact, hgt_ref, hgt, val_ref, val, spline, logwgt)
Definition: scale_interp.F90:1263
scale_interp::cross
real(rp) function cross(x0, y0, x1, y1)
cross product
Definition: scale_interp.F90:2346
scale_prof
module profiler
Definition: scale_prof.F90:11
scale_interp::interp_factor3d_linear_latlon
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)
Definition: scale_interp.F90:961
scale_interp::interp_interp3d
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)
Definition: scale_interp.F90:1469
scale_interp::interp_factor2d_linear_xy
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)
Definition: scale_interp.F90:483
scale_const::const_radius
real(rp), public const_radius
radius of the planet [m]
Definition: scale_const.F90:47
scale_interp::interp_factor2d_linear_latlon
subroutine, public interp_factor2d_linear_latlon(IA_ref, JA_ref, IA, JA, lon_ref, lat_ref, lon, lat, idx_i, idx_j, hfact)
Definition: scale_interp.F90:357
scale_matrix
module MATRIX
Definition: scale_matrix.F90:17
scale_prof::prof_rapend
subroutine, public prof_rapend(rapname_base, level, disable_barrier)
Save raptime.
Definition: scale_prof.F90:246
scale_interp::interp_setup
subroutine, public interp_setup(weight_order, search_limit)
Setup.
Definition: scale_interp.F90:90
scale_interp::interp_search_horiz_struct
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)
Definition: scale_interp.F90:1880
scale_interp::interp_factor3d_weight
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)
Definition: scale_interp.F90:1135
scale_interp::interp_bilinear_inv
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)
Definition: scale_interp.F90:2199
scale_const::const_d2r
real(rp), public const_d2r
degree to radian
Definition: scale_const.F90:33
scale_interp::interp_div_block
subroutine interp_div_block(nsize, psize, nidx_max, lon_ref, lat_ref, idx, nidx, lon_min, lon_max, lat_min, lat_max, dlon, dlat)
Definition: scale_interp.F90:2120
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:43
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:57
Generated by   doxygen 1.8.17