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