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  !-----------------------------------------------------------------------------
82 contains
83  !-----------------------------------------------------------------------------
85  subroutine interp_setup( &
86  weight_order, &
87  search_limit )
88  use scale_prc, only: &
89  prc_abort
90  implicit none
91 
92  integer, intent(in) :: weight_order
93  real(rp), intent(in), optional :: search_limit
94 
95  namelist /param_interp/ &
96  interp_buffer_size_fact, &
97  interp_use_spline_vert, &
98  interp_threshold_undef
99 
100  integer :: ierr
101  !---------------------------------------------------------------------------
102 
103  log_newline
104  log_info("INTERP_setup",*) 'Setup'
105 
106  interp_weight_order = weight_order
107 
108  interp_search_limit = large_number
109  if ( present(search_limit) ) then
110  interp_search_limit = search_limit
111  log_info("INTERP_setup",*) 'search limit [m] : ', interp_search_limit
112  endif
113 
114  !--- read namelist
115  rewind(io_fid_conf)
116  read(io_fid_conf,nml=param_interp,iostat=ierr)
117  if( ierr < 0 ) then !--- missing
118  log_info("INTERP_setup",*) 'Not found namelist. Default used.'
119  elseif( ierr > 0 ) then !--- fatal error
120  log_error("INTERP_setup",*) 'Not appropriate names in namelist PARAM_INTERP. Check!'
121  call prc_abort
122  endif
123  log_nml(param_interp)
124 
125  if ( rp == 8 ) then
126  eps_bilinear = 1e-6_rp
127  else
128  eps_bilinear = 1e-2_rp
129  end if
130 
131  return
132  end subroutine interp_setup
133 
134  !-----------------------------------------------------------------------------
135  subroutine interp_domain_compatibility( &
136  lon_org, &
137  lat_org, &
138  topc_org, &
139  lon_loc, &
140  lat_loc, &
141  topc_loc, &
142  topf_loc, &
143  skip_x, &
144  skip_y, &
145  skip_z )
146  use scale_const, only: &
147  d2r => const_d2r
148  use scale_prc, only: &
149  prc_abort
150  implicit none
151 
152  real(rp), intent(in) :: lon_org (:,:)
153  real(rp), intent(in) :: lat_org (:,:)
154  real(rp), intent(in) :: topc_org(:,:) ! full level
155  real(rp), intent(in) :: lon_loc (:,:)
156  real(rp), intent(in) :: lat_loc (:,:)
157  real(rp), intent(in) :: topc_loc(:,:) ! full level
158  real(rp), intent(in) :: topf_loc(:,:) ! half level (ceil)
159 
160  logical, intent(in), optional :: skip_z
161  logical, intent(in), optional :: skip_x
162  logical, intent(in), optional :: skip_y
163 
164  real(rp) :: max_ref, min_ref
165  real(rp) :: max_loc, min_loc
166 
167  logical :: skip_z_
168  logical :: skip_x_
169  logical :: skip_y_
170 
171  intrinsic size
172  !---------------------------------------------------------------------------
173 
174  skip_z_ = .false.
175  if ( present(skip_z) ) then
176  skip_z_ = skip_z
177  endif
178 
179  skip_x_ = .false.
180  if ( present(skip_x) ) then
181  skip_x_ = skip_x
182  endif
183 
184  skip_y_ = .false.
185  if ( present(skip_y) ) then
186  skip_y_ = skip_y
187  endif
188 
189  if ( .NOT. skip_z_ ) then
190  max_ref = maxval( topc_org(:,:) ) + minval( topf_loc(:,:) - topc_loc(:,:) ) ! not real top boundary (only for check)
191  max_loc = maxval( topc_loc(:,:) )
192 
193  if ( max_ref < max_loc ) then
194  log_error("INTERP_domain_compatibility",*) 'REQUESTED DOMAIN IS TOO MUCH BROAD'
195  log_error_cont(*) '-- VERTICAL direction over the limit'
196  log_error_cont(*) '-- reference max: ', max_ref
197  log_error_cont(*) '-- local max: ', max_loc
198  call prc_abort
199  endif
200  endif
201 
202  if ( .NOT. skip_x_ ) then
203  max_ref = maxval( lon_org(:,:) / d2r )
204  min_ref = minval( lon_org(:,:) / d2r )
205  max_loc = maxval( lon_loc(:,:) / d2r )
206  min_loc = minval( lon_loc(:,:) / d2r )
207 
208  if ( (min_ref+360.0_rp-max_ref) < 360.0_rp / size(lon_org,1) * 2.0_rp ) then
209  ! cyclic OK
210  elseif( max_ref < max_loc .OR. min_ref > min_loc ) then
211  log_error("INTERP_domain_compatibility",*) 'REQUESTED DOMAIN IS TOO MUCH BROAD'
212  log_error_cont(*) '-- LONGITUDINAL direction over the limit'
213  log_error_cont(*) '-- reference max: ', max_ref
214  log_error_cont(*) '-- reference min: ', min_ref
215  log_error_cont(*) '-- local max: ', max_loc
216  log_error_cont(*) '-- local min: ', min_loc
217  call prc_abort
218  endif
219  endif
220 
221  if ( .NOT. skip_y_ ) then
222  max_ref = maxval( lat_org(:,:) / d2r )
223  min_ref = minval( lat_org(:,:) / d2r )
224  max_loc = maxval( lat_loc(:,:) / d2r )
225  min_loc = minval( lat_loc(:,:) / d2r )
226 
227  if ( max_ref < max_loc .OR. min_ref > min_loc ) then
228  log_error("INTERP_domain_compatibility",*) 'REQUESTED DOMAIN IS TOO MUCH BROAD'
229  log_error_cont(*) '-- LATITUDINAL direction over the limit'
230  log_error_cont(*) '-- reference max: ', max_ref
231  log_error_cont(*) '-- reference min: ', min_ref
232  log_error_cont(*) '-- local max: ', max_loc
233  log_error_cont(*) '-- local min: ', min_loc
234  call prc_abort
235  endif
236  endif
237 
238  return
239  end subroutine interp_domain_compatibility
240 
241  !-----------------------------------------------------------------------------
242  ! vertical search of interpolation points for two-points
243 !OCL SERIAL
244  subroutine interp_factor1d( &
245  KA_ref, KS_ref, KE_ref, &
246  KA, KS, KE, &
247  hgt_ref, &
248  hgt, &
249  idx_k, &
250  vfact, &
251  flag_extrap )
252  use scale_prc, only: &
253  prc_abort
254  use scale_const, only: &
255  undef => const_undef, &
256  eps => const_eps
257  implicit none
258  integer, intent(in) :: ka_ref, ks_ref, ke_ref ! number of z-direction (reference)
259  integer, intent(in) :: ka, ks, ke ! number of z-direction (target)
260 
261  real(rp), intent(in) :: hgt_ref(ka_ref) ! height [m] (reference)
262  real(rp), intent(in) :: hgt (ka) ! height [m] (target)
263 
264  integer, intent(out) :: idx_k (ka,2) ! k-index in reference (target)
265  real(rp), intent(out) :: vfact (ka) ! horizontal interp factor (target)
266 
267  logical, intent(in), optional :: flag_extrap ! when true, extrapolation will be executed (just copy)
268 
269  logical :: flag_extrap_
270 
271  integer :: k, kk
272  !---------------------------------------------------------------------------
273 
274  if ( present(flag_extrap) ) then
275  flag_extrap_ = flag_extrap
276  else
277  flag_extrap_ = .true.
278  end if
279 
280  do k = ks, ke
281  idx_k(k,1) = -1
282  idx_k(k,2) = -1
283 
284  if ( hgt(k) < hgt_ref(ks_ref) - eps ) then
285  if ( flag_extrap_ ) then
286  idx_k(k,1) = ks_ref
287  idx_k(k,2) = -1
288  vfact(k) = 1.0_rp
289  else
290  idx_k(k,1) = -1
291  idx_k(k,2) = -1
292  vfact(k) = undef
293  end if
294  elseif( hgt(k) < hgt_ref(ks_ref) ) then
295  idx_k(k,1) = ks_ref
296  idx_k(k,2) = -1
297  vfact(k) = 1.0_rp
298  elseif( hgt(k) > hgt_ref(ke_ref) + eps ) then
299  if ( flag_extrap_ ) then
300  idx_k(k,1) = ke_ref
301  idx_k(k,2) = -1
302  vfact(k) = 1.0_rp
303  else
304  idx_k(k,1) = -1
305  idx_k(k,2) = -1
306  vfact(k) = undef
307  end if
308  elseif( hgt(k) >= hgt_ref(ke_ref) ) then
309  idx_k(k,1) = ke_ref
310  idx_k(k,2) = -1
311  vfact(k) = 1.0_rp
312  else
313  do kk = ks_ref, ke_ref-1
314  if ( hgt(k) >= hgt_ref(kk ) &
315  .AND. hgt(k) < hgt_ref(kk+1) ) then
316  idx_k(k,1) = kk
317  idx_k(k,2) = kk + 1
318  vfact(k) = ( hgt_ref(kk+1) - hgt(k) ) &
319  / ( hgt_ref(kk+1) - hgt_ref(kk) )
320 
321  exit
322  endif
323  enddo
324  endif
325 
326  enddo ! k-loop
327 
328  return
329  end subroutine interp_factor1d
330 
331 
332  !-----------------------------------------------------------------------------
333  ! make interpolation factor using bi-linear method on the latlon coordinate
334  ! This can be used only for structured grid
335  subroutine interp_factor2d_linear_latlon( &
336  IA_ref, JA_ref, &
337  IA, JA, &
338  lon_ref, lat_ref, &
339  lon, lat, &
340  idx_i, idx_j, hfact )
341  use scale_sort, only: &
342  sort_exec
343  implicit none
344  integer, intent(in) :: ia_ref, ja_ref
345  integer, intent(in) :: ia, ja
346 
347  real(rp), intent(in) :: lon_ref(ia_ref)
348  real(rp), intent(in) :: lat_ref(ja_ref)
349  real(rp), intent(in) :: lon(ia,ja)
350  real(rp), intent(in) :: lat(ia,ja)
351  integer, intent(out) :: idx_i(ia,ja,4)
352  integer, intent(out) :: idx_j(ia,ja,4)
353  real(rp), intent(out) :: hfact(ia,ja,4)
354 
355  real(rp) :: f1, f2
356  integer :: i, j, ii, jj
357 
358  call prof_rapstart('INTERP_fact',3)
359 
360  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
361  !$omp private(f1,f2)
362  do j = 1, ja
363  do i = 1, ia
364 
365  ! longitude
366  if ( lon(i,j) < lon_ref(1) ) then
367  idx_i(i,j,:) = 1
368  hfact(i,j,:) = 0.0_rp
369  else if ( lon(i,j) > lon_ref(ia_ref) ) then
370  idx_i(i,j,:) = ia_ref
371  hfact(i,j,:) = 0.0_rp
372  else if ( lon(i,j) == lon_ref(ia_ref) ) then
373  idx_i(i,j,:) = ia_ref
374  hfact(i,j,1) = 0.0_rp
375  hfact(i,j,2) = 1.0_rp
376  hfact(i,j,3) = 0.0_rp
377  hfact(i,j,4) = 1.0_rp
378  else
379  do ii = 1, ia_ref
380  if ( lon(i,j) < lon_ref(ii) ) then
381  f1 = ( lon_ref(ii) - lon(i,j) ) / ( lon_ref(ii) - lon_ref(ii-1) )
382  f2 = ( lon(i,j) - lon_ref(ii-1) ) / ( lon_ref(ii) - lon_ref(ii-1) )
383  hfact(i,j,1) = f1
384  hfact(i,j,2) = f2
385  hfact(i,j,3) = f1
386  hfact(i,j,4) = f2
387  idx_i(i,j,1) = ii-1
388  idx_i(i,j,2) = ii
389  idx_i(i,j,3) = ii-1
390  idx_i(i,j,4) = ii
391  exit
392  end if
393  end do
394  end if
395 
396  ! latitude
397  if ( lat(i,j) < lat_ref(1) ) then
398  idx_j(i,j,:) = 1
399  hfact(i,j,:) = 0.0_rp
400  else if ( lat(i,j) > lat_ref(ja_ref) ) then
401  idx_j(i,j,:) = ja_ref
402  hfact(i,j,:) = 0.0_rp
403  else if ( lat(i,j) == lat_ref(ja_ref) ) then
404  idx_j(i,j,:) = ja_ref
405  hfact(i,j,1) = 0.0_rp
406  hfact(i,j,2) = 0.0_rp
407  !hfact(i,j,3) = hfact(i,j,3)
408  !hfact(i,j,4) = hfact(i,j,4)
409  else
410  do jj = 1, ja_ref
411  if ( lat(i,j) < lat_ref(jj) ) then
412  f1 = ( lat_ref(jj) - lat(i,j) ) / ( lat_ref(jj) - lat_ref(jj-1) )
413  f2 = ( lat(i,j) - lat_ref(jj-1) ) / ( lat_ref(jj) - lat_ref(jj-1) )
414  hfact(i,j,1) = hfact(i,j,1) * f1
415  hfact(i,j,2) = hfact(i,j,2) * f1
416  hfact(i,j,3) = hfact(i,j,3) * f2
417  hfact(i,j,4) = hfact(i,j,4) * f2
418  idx_j(i,j,1) = jj-1
419  idx_j(i,j,2) = jj-1
420  idx_j(i,j,3) = jj
421  idx_j(i,j,4) = jj
422  exit
423  end if
424  end do
425  end if
426 
427  call sort_exec( 4, & ! [IN]
428  hfact(i,j,:), idx_i(i,j,:), idx_j(i,j,:), & ! [INOUT]
429  reverse = .true. ) ! [IN]
430 
431  end do
432  end do
433 
434  call prof_rapend ('INTERP_fact',3)
435 
436  return
437  end subroutine interp_factor2d_linear_latlon
438 
439  !-----------------------------------------------------------------------------
440  ! make interpolation factor using bi-linear method on the xy coordinate
441  ! This can be used only for structured grid
442  subroutine interp_factor2d_linear_xy( &
443  IA_ref, JA_ref, &
444  IA, JA, &
445  x_ref, y_ref, &
446  x, y, &
447  idx_i, idx_j, &
448  hfact, &
449  zonal, pole )
450  use scale_prc, only: &
451  prc_abort
452  use scale_const, only: &
453  undef => const_undef, &
454  eps => const_eps
455  use scale_sort, only: &
456  sort_exec
457  implicit none
458  integer, intent(in) :: ia_ref, ja_ref
459  integer, intent(in) :: ia, ja
460 
461  real(rp), intent(in) :: x_ref(ia_ref,ja_ref)
462  real(rp), intent(in) :: y_ref(ia_ref,ja_ref)
463  real(rp), intent(in) :: x(ia)
464  real(rp), intent(in) :: y(ja)
465 
466  integer, intent(out) :: idx_i(ia,ja,4)
467  integer, intent(out) :: idx_j(ia,ja,4)
468  real(rp), intent(out) :: hfact(ia,ja,4)
469 
470  logical, intent(in), optional :: zonal
471  logical, intent(in), optional :: pole
472 
473  real(rp) :: u, v
474  integer :: inc_i, inc_j
475  logical :: error, err
476  logical :: zonal_, pole_
477 
478  integer :: ii0, jj0
479  integer :: ii, jj
480  integer :: i1, i2, i3, i4
481  integer :: j1, j2, j3, j4
482  integer :: i, j
483  integer :: ite, ite_max
484 
485  call prof_rapstart('INTERP_fact',3)
486 
487  if ( present(zonal) ) then
488  zonal_ = zonal
489  else
490  zonal_ = .false.
491  end if
492  if ( present(pole) ) then
493  pole_ = pole
494  else
495  pole_ = .false.
496  end if
497 
498  ii0 = ia_ref * 0.5_rp
499  jj0 = ja_ref * 0.5_rp
500  error = .false.
501 
502  ite_max = ia_ref * ja_ref / 4
503 
504  !$omp parallel do &
505  !$omp private(inc_i,inc_j,ii,jj,i1,i2,i3,i4,j1,j2,j3,j4,u,v,err) &
506  !$omp firstprivate(ii0,jj0)
507  do j = 1, ja
508  do i = 1, ia
509 
510  if ( i==1 ) then
511  ii = ii0
512  jj = jj0
513  end if
514  do ite = 1, ite_max
515  i1 = ii
516  i2 = ii + 1
517  i3 = ii + 1
518  i4 = ii
519  j1 = jj
520  j2 = jj
521  j3 = jj + 1
522  j4 = jj + 1
523  if ( ii == ia_ref ) then ! zonal_ must be .true.
524  i2 = 1
525  i3 = 1
526  end if
527  if ( jj == 0 ) then ! pole_ must be .true.
528  j1 = 1
529  j2 = 1
530  i2 = ia_ref / 2 + ii
531  if ( i2 > ia_ref ) i2 = i2 - ia_ref
532  i1 = i2 + 1
533  if ( i1 > ia_ref ) i1 = i1 - ia_ref
534  end if
535  if ( jj == ja_ref ) then ! pole_ must be .true.
536  j3 = ja_ref
537  j4 = ja_ref
538  i3 = ia_ref / 2 + ii
539  if ( i3 > ia_ref ) i3 = i3 - ia_ref
540  i4 = i3 + 1
541  if ( i4 > ia_ref ) i4 = i4 - ia_ref
542  end if
543  if ( abs( x_ref(i1,j1) - undef ) < eps .or. &
544  abs( x_ref(i2,j2) - undef ) < eps .or. &
545  abs( x_ref(i3,j3) - undef ) < eps .or. &
546  abs( x_ref(i4,j4) - undef ) < eps ) then
547  if ( ii == ia_ref-1 ) then
548  ii = 1
549  if ( jj == ja_ref-1 ) then
550  jj = 1
551  else
552  jj = jj + 2
553  end if
554  else
555  ii = ii + 2
556  end if
557  if ( ii == ia_ref ) then
558  ii = 1
559  jj = jj + 2
560  end if
561  if ( jj == ja_ref ) then
562  jj = 1
563  ii = ii + 2
564  if ( ii == ia_ref ) ii = 1
565  end if
566  else
567  call interp_check_inside( &
568  x_ref(i1,j1), x_ref(i2,j2), x_ref(i3,j3), x_ref(i4,j4), & ! (in)
569  y_ref(i1,j1), y_ref(i2,j2), y_ref(i3,j3), y_ref(i4,j4), & ! (in)
570  x(i), y(j), & ! (in)
571  inc_i, inc_j, & ! (out)
572  err ) ! (out)
573  if ( err ) error = .true.
574  if ( error ) exit
575  if ( inc_i == 0 .and. inc_j == 0 ) then ! inside the quadrilateral
576  call interp_bilinear_inv( &
577  x_ref(i1,j1), x_ref(i2,j2), x_ref(i3,j3), x_ref(i4,j4), & ! (in)
578  y_ref(i1,j1), y_ref(i2,j2), y_ref(i3,j3), y_ref(i4,j4), & ! (in)
579  x(i), y(j), & ! (in)
580  u, v, & ! (out)
581  err ) ! (out)
582  if ( err ) error = .true.
583  if ( error ) exit
584  idx_i(i,j,1) = i1
585  idx_i(i,j,2) = i2
586  idx_i(i,j,3) = i3
587  idx_i(i,j,4) = i4
588  idx_j(i,j,1) = j1
589  idx_j(i,j,2) = j2
590  idx_j(i,j,3) = j3
591  idx_j(i,j,4) = j4
592  hfact(i,j,1) = ( 1.0_rp - u ) * ( 1.0_rp - v )
593  hfact(i,j,2) = ( u ) * ( 1.0_rp - v )
594  hfact(i,j,3) = ( u ) * ( v )
595  hfact(i,j,4) = ( 1.0_rp - u ) * ( v )
596  exit
597  end if
598  ii = ii + inc_i
599  jj = jj + inc_j
600  if ( zonal_ .or. pole_ ) then
601  if ( ii == 0 ) ii = ia_ref
602  if ( ii == ia_ref+1 ) ii = 1
603  if ( pole_ ) then
604  jj = max( jj, 0 )
605  jj = min( jj, ja_ref )
606  else
607  jj = max( jj, 1 )
608  jj = min( jj, ja_ref-1 )
609  end if
610  else
611  if ( ii == 0 .and. jj == 0 ) then
612  ii = ia_ref - 1
613  jj = ja_ref - 1
614  end if
615  if ( ii == 0 .and. jj == ja_ref ) then
616  ii = ia_ref - 1
617  jj = 1
618  end if
619  if ( ii == ia_ref .and. jj == 0 ) then
620  ii = 1
621  jj = ja_ref - 1
622  end if
623  if ( ii == ia_ref .and. jj == ja_ref ) then
624  ii = 1
625  jj = 1
626  end if
627  ii = max( min( ii, ia_ref-1 ), 1 )
628  jj = max( min( jj, ja_ref-1 ), 1 )
629  end if
630  end if
631  end do
632  if ( ite == ite_max+1 ) then
633  log_error("INTERP_factor2d_linear_xy",*) 'iteration max has been reached', i, j, x(i), y(j)
634  log_error_cont(*) minval(x_ref), maxval(x_ref), minval(y_ref), maxval(y_ref)
635  log_error_cont(*) x_ref(1,1), x_ref(ia_ref,1),x_ref(1,ja_ref)
636  log_error_cont(*) y_ref(1,1), y_ref(ia_ref,1),y_ref(1,ja_ref)
637  idx_i(i,j,:) = 1
638  idx_j(i,j,:) = 1
639  hfact(i,j,:) = 0.0_rp
640  call prc_abort
641  end if
642 
643  if ( i==1 ) then
644  ii0 = ii
645  jj0 = jj
646  end if
647 
648  if ( error ) exit
649 
650  call sort_exec( 4, & ! [IN]
651  hfact(i,j,:), idx_i(i,j,:), idx_j(i,j,:), & ! [INOUT]
652  reverse = .true. ) ! [IN]
653 
654  end do
655  end do
656 
657  if ( error ) call prc_abort
658 
659  call prof_rapend ('INTERP_fact',3)
660 
661  return
662  end subroutine interp_factor2d_linear_xy
663 
664  !-----------------------------------------------------------------------------
665  ! make interpolation factor using Lat-Lon with weight based on the distance
666  subroutine interp_factor2d_weight( &
667  npoints, &
668  IA_ref, JA_ref, &
669  IA, JA, &
670  lon_ref,lat_ref, &
671  lon, lat, &
672  idx_i, idx_j, hfact, &
673  search_limit, &
674  latlon_structure, &
675  lon_1d, lat_1d, &
676  weight_order )
677  use scale_prc, only: &
678  prc_abort
679  implicit none
680  integer, intent(in) :: npoints ! number of interpolation point for horizontal
681  integer, intent(in) :: ia_ref ! number of x-direction (reference)
682  integer, intent(in) :: ja_ref ! number of y-direction (reference)
683  integer, intent(in) :: ia ! number of x-direction (target)
684  integer, intent(in) :: ja ! number of y-direction (target)
685  real(rp), intent(in) :: lon_ref(ia_ref,ja_ref) ! longitude [rad] (reference)
686  real(rp), intent(in) :: lat_ref(ia_ref,ja_ref) ! latitude [rad] (reference)
687  real(rp), intent(in) :: lon (ia,ja) ! longitude [rad] (target)
688  real(rp), intent(in) :: lat (ia,ja) ! latitude [rad] (target)
689  integer, intent(out) :: idx_i(ia,ja,npoints) ! i-index in reference (target)
690  integer, intent(out) :: idx_j(ia,ja,npoints) ! j-index in reference (target)
691  real(rp), intent(out) :: hfact(ia,ja,npoints) ! horizontal interp factor (target)
692 
693  real(rp), intent(in), optional :: search_limit
694  logical, intent(in), optional :: latlon_structure
695  real(rp), intent(in), optional :: lon_1d(ia_ref), lat_1d(ja_ref)
696  integer, intent(in), optional :: weight_order
697 
698  logical :: ll_struct_
699 
700  real(rp) :: lon_min, lon_max
701  real(rp) :: lat_min, lat_max
702  real(rp) :: dlon, dlat
703 
704  ! for structure grid
705  integer :: psizex, psizey
706  real(rp) :: lon0, lat0
707  integer, allocatable :: i0(:), i1(:), j0(:), j1(:)
708 
709  ! for unstructure grid
710  integer :: nsize, psize, nidx_max
711  integer, allocatable :: idx_blk(:,:,:), nidx(:,:)
712  integer :: idx_ref(npoints)
713 
714 
715  integer :: i, j, ii, jj, n
716  !---------------------------------------------------------------------------
717 
718  call prof_rapstart('INTERP_fact',3)
719 
720  if ( present(latlon_structure) ) then
721  ll_struct_ = latlon_structure
722  else
723  ll_struct_ = .false.
724  end if
725 
726  hfact(:,:,:) = 0.0_rp
727 
728 
729  if ( ll_struct_ ) then
730 
731  if ( ia_ref > 10 ) then
732  psizex = int( 2.0_rp*sqrt(real(ia_ref+1,rp)) )
733  else
734  psizex = 1
735  end if
736  if ( ja_ref > 10 ) then
737  psizey = int( 2.0_rp*sqrt(real(ja_ref+1,rp)) )
738  else
739  psizey = 1
740  end if
741 
742  allocate( i0(psizex), i1(psizex) )
743  allocate( j0(psizey), j1(psizey) )
744 
745  dlon = ( lon_max - lon_min ) / psizex
746  dlat = ( lat_max - lat_min ) / psizey
747 
748  do ii = 1, psizex
749  lon0 = lon_min + dlon * (ii-1)
750  do i = 1, ia_ref
751  if ( lon_1d(i) >= lon0 ) then
752  i0(ii) = i
753  exit
754  end if
755  end do
756  end do
757  do ii = 1, psizex-1
758  i1(ii) = i0(ii+1) - 1
759  end do
760  i1(psizex) = ia_ref
761 
762  do jj = 1, psizey
763  lat0 = lat_min + dlat * (jj-1)
764  do j = 1, ja_ref
765  if ( lat_1d(j) >= lat0 ) then
766  j0(jj) = j
767  exit
768  end if
769  end do
770  end do
771  do jj = 1, psizey-1
772  j1(jj) = j0(jj+1) - 1
773  end do
774  j1(psizey) = ja_ref
775 
776  !$omp parallel do OMP_SCHEDULE_ collapse(2)
777  do j = 1, ja
778  do i = 1, ia
779  ! main search
780  call interp_search_horiz_struct( npoints, & ! [IN]
781  psizex, psizey, & ! [IN]
782  ia_ref, ja_ref, & ! [IN]
783  lon_1d(:), lat_1d(:), & ! [IN]
784  lon_min, lat_min, & ! [IN]
785  dlon, dlat, & ! [IN]
786  i0(:), i1(:), j0(:), j1(:), & ! [IN]
787  lon(i,j), lat(i,j), & ! [IN]
788  idx_i(i,j,:), idx_j(i,j,:), & ! [OUT]
789  hfact(i,j,:), & ! [OUT]
790  search_limit = search_limit, & ! [IN]
791  weight_order = weight_order ) ! [IN]
792  enddo
793  enddo
794 
795  deallocate( i0, i1, j0, j1 )
796 
797  else
798 
799  nsize = ia_ref * ja_ref
800  if ( nsize > 100 ) then
801  psize = int( sqrt(2.0_rp*sqrt(real(nsize,rp))) )
802  nidx_max = nsize / psize * interp_buffer_size_fact
803  else
804  psize = 1
805  nidx_max = nsize
806  end if
807 
808  allocate(idx_blk(nidx_max,psize,psize))
809  allocate(nidx( psize,psize))
810 
811  call interp_div_block(nsize, psize, nidx_max, & ! [IN]
812  lon_ref(:,:), lat_ref(:,:), & ! [IN]
813  idx_blk(:,:,:), nidx(:,:), & ! [OUT]
814  lon_min, lon_max, & ! [OUT]
815  lat_min, lat_max, & ! [OUT]
816  dlon, dlat ) ! [OUT]
817 
818  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
819  !$omp private(idx_ref)
820  do j = 1, ja
821  do i = 1, ia
822  ! main search
823  call interp_search_horiz( npoints, & ! [IN]
824  nsize, & ! [IN]
825  lon_ref(:,:), lat_ref(:,:), & ! [IN]
826  lon_min, lon_max, & ! [IN]
827  lat_min, lat_max, & ! [IN]
828  psize, nidx_max, & ! [IN]
829  dlon, dlat, & ! [IN]
830  idx_blk(:,:,:), nidx(:,:), & ! [IN]
831  lon(i,j), lat(i,j), & ! [IN]
832  idx_ref(:), & ! [OUT]
833  hfact(i,j,:), & ! [OUT]
834  search_limit = search_limit, & ! [IN]
835  weight_order = weight_order ) ! [IN]
836  do n = 1, npoints
837  idx_i(i,j,n) = mod(idx_ref(n) - 1, ia_ref) + 1
838  idx_j(i,j,n) = ( idx_ref(n) - 1 ) / ia_ref + 1
839  end do
840  enddo
841  enddo
842 
843  deallocate(idx_blk, nidx)
844 
845  end if
846 
847  call prof_rapend ('INTERP_fact',3)
848 
849  return
850  end subroutine interp_factor2d_weight
851 
852  !-----------------------------------------------------------------------------
853  ! make interpolation factor using bi-linear method for the horizontal direction on latlon grid
854  ! This can be used only for structured grid
855  subroutine interp_factor3d_linear_latlon( &
856  KA_ref, KS_ref, KE_ref, &
857  IA_ref, JA_ref, &
858  KA, KS, KE, &
859  IA, JA, &
860  lon_ref, lat_ref, &
861  hgt_ref, &
862  lon, lat, &
863  hgt, &
864  idx_i, idx_j, &
865  hfact, &
866  idx_k, &
867  vfact, &
868  flag_extrap )
869  implicit none
870 
871  integer, intent(in) :: ka_ref, ks_ref, ke_ref ! number of z-direction (reference)
872  integer, intent(in) :: ia_ref ! number of x-direction (reference)
873  integer, intent(in) :: ja_ref ! number of y-direction (reference)
874  integer, intent(in) :: ka, ks, ke ! number of z-direction (target)
875  integer, intent(in) :: ia ! number of x-direction (target)
876  integer, intent(in) :: ja ! number of y-direction (target)
877  real(rp), intent(in) :: lon_ref(ia_ref) ! longitude (reference)
878  real(rp), intent(in) :: lat_ref(ja_ref) ! latitude (reference)
879  real(rp), intent(in) :: hgt_ref(ka_ref,ia_ref,ja_ref) ! height [m] (reference)
880  real(rp), intent(in) :: lon(ia,ja) ! longitude (target)
881  real(rp), intent(in) :: lat(ia,ja) ! latitude (target)
882  real(rp), intent(in) :: hgt(ka,ia,ja) ! longitude [m] (target)
883 
884  integer, intent(out) :: idx_i(ia,ja,4) ! i-index in reference (target)
885  integer, intent(out) :: idx_j(ia,ja,4) ! j-index in reference (target)
886  real(rp), intent(out) :: hfact(ia,ja,4) ! horizontal interp factor (target)
887  integer, intent(out) :: idx_k(ka,2,ia,ja,4) ! k-index in reference (target)
888  real(rp), intent(out) :: vfact(ka, ia,ja,4) ! vertical interp factor (target)
889 
890  logical, intent(in), optional :: flag_extrap ! when true, vertical extrapolation will be executed (just copy)
891 
892  integer :: i, j, ii, jj, n
893 
894  call interp_factor2d_linear_latlon( &
895  ia_ref, ja_ref, & ! [IN]
896  ia, ja, & ! [IN]
897  lon_ref(:), lat_ref(:), & ! [IN]
898  lon(:,:), lat(:,:), & ! [IN]
899  idx_i(:,:,:), idx_j(:,:,:), hfact(:,:,:) ) ! [OUT]
900 
901  call prof_rapstart('INTERP_fact',3)
902 
903  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
904  !$omp private(ii,jj)
905  do j = 1, ja
906  do i = 1, ia
907  do n = 1, 4
908  ii = idx_i(i,j,n)
909  jj = idx_j(i,j,n)
910  call interp_factor1d( ka_ref, ks_ref, ke_ref, & ! [IN]
911  ka, ks, ke, & ! [IN]
912  hgt_ref(:,ii,jj), & ! [IN]
913  hgt(:,i,j), & ! [IN]
914  idx_k(:,:,i,j,n), & ! [OUT]
915  vfact(:, i,j,n), & ! [OUT]
916  flag_extrap = flag_extrap ) ! [IN, optional]
917 
918  enddo
919  enddo
920  enddo
921 
922  call prof_rapend('INTERP_fact',3)
923 
924  return
925  end subroutine interp_factor3d_linear_latlon
926 
927  !-----------------------------------------------------------------------------
928  ! make interpolation factor using bi-linear method for the horizontal direction on the xy coordinate
929  ! This can be used only for structured grid
930  subroutine interp_factor3d_linear_xy( &
931  KA_ref, KS_ref, KE_ref, &
932  IA_ref, JA_ref, &
933  KA, KS, KE, &
934  IA, JA, &
935  x_ref, y_ref, &
936  hgt_ref, &
937  x, y, &
938  hgt, &
939  idx_i, idx_j, &
940  hfact, &
941  idx_k, &
942  vfact, &
943  flag_extrap, &
944  zonal, pole )
945  implicit none
946  integer, intent(in) :: ka_ref, ks_ref, ke_ref ! number of z-direction (reference)
947  integer, intent(in) :: ia_ref ! number of x-direction (reference)
948  integer, intent(in) :: ja_ref ! number of y-direction (reference)
949  integer, intent(in) :: ka, ks, ke ! number of z-direction (target)
950  integer, intent(in) :: ia ! number of x-direction (target)
951  integer, intent(in) :: ja ! number of y-direction (target)
952  real(rp), intent(in) :: x_ref(ia_ref,ja_ref) ! x point (reference)
953  real(rp), intent(in) :: y_ref(ia_ref,ja_ref) ! y point (reference)
954  real(rp), intent(in) :: hgt_ref(ka_ref,ia_ref,ja_ref) ! height [m] (reference)
955  real(rp), intent(in) :: x (ia) ! x point (target)
956  real(rp), intent(in) :: y (ja) ! y point (target)
957  real(rp), intent(in) :: hgt(ka,ia,ja) ! longitude [m] (target)
958 
959  integer, intent(out) :: idx_i(ia,ja,4) ! i-index in reference
960  integer, intent(out) :: idx_j(ia,ja,4) ! j-index in reference
961  real(rp), intent(out) :: hfact(ia,ja,4) ! horizontal interp factor
962  integer, intent(out) :: idx_k(ka,2,ia,ja,4) ! k-index in reference
963  real(rp), intent(out) :: vfact(ka, ia,ja,4) ! vertical interp factor
964 
965  logical, intent(in), optional :: flag_extrap ! when true, vertical extrapolation will be executed (just copy)
966  logical, intent(in), optional :: zonal
967  logical, intent(in), optional :: pole
968 
969  integer :: i, j, ii, jj, n
970 
971  call interp_factor2d_linear_xy( &
972  ia_ref, ja_ref, & ! [IN]
973  ia, ja, & ! [IN]
974  x_ref(:,:), y_ref(:,:), & ! [IN]
975  x(:), y(:), & ! [IN]
976  idx_i(:,:,:), idx_j(:,:,:), hfact(:,:,:), & ! [OUT]
977  zonal = zonal, pole = pole ) ! [IN]
978 
979  call prof_rapstart('INTERP_fact',3)
980 
981  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
982  !$omp private(ii,jj)
983  do j = 1, ja
984  do i = 1, ia
985  do n = 1, 4
986  ii = idx_i(i,j,n)
987  jj = idx_j(i,j,n)
988  call interp_factor1d( ka_ref, ks_ref, ke_ref, & ! [IN]
989  ka, ks, ke, & ! [IN]
990  hgt_ref(:,ii,jj), & ! [IN]
991  hgt(:,i,j), & ! [IN]
992  idx_k(:,:,i,j,n), & ! [OUT]
993  vfact(: ,i,j,n), & ! [OUT]
994  flag_extrap = flag_extrap ) ! [IN, optional]
995 
996  enddo
997  enddo
998  enddo
999 
1000  call prof_rapend('INTERP_fact',3)
1001 
1002  return
1003  end subroutine interp_factor3d_linear_xy
1004 
1005  !-----------------------------------------------------------------------------
1006  ! make interpolation factor using Lat-Lon and Z-Height information
1007  subroutine interp_factor3d_weight( &
1008  npoints, &
1009  KA_ref, KS_ref, KE_ref, &
1010  IA_ref, JA_ref, &
1011  KA, KS, KE, &
1012  IA, JA, &
1013  lon_ref, lat_ref, &
1014  hgt_ref, &
1015  lon, lat, &
1016  hgt, &
1017  idx_i, idx_j, &
1018  hfact, &
1019  idx_k, &
1020  vfact, &
1021  flag_extrap )
1022  use scale_prc, only: &
1023  prc_abort
1024  implicit none
1025  integer, intent(in) :: npoints ! number of interpolation point for horizontal
1026  integer, intent(in) :: ka_ref, ks_ref, ke_ref ! number of z-direction (reference)
1027  integer, intent(in) :: ia_ref ! number of x-direction (reference)
1028  integer, intent(in) :: ja_ref ! number of y-direction (reference)
1029  integer, intent(in) :: ka, ks, ke ! number of z-direction (target)
1030  integer, intent(in) :: ia ! number of x-direction (target)
1031  integer, intent(in) :: ja ! number of y-direction (target)
1032 
1033  real(rp), intent(in) :: lon_ref(ia_ref,ja_ref) ! longitude [rad] (reference)
1034  real(rp), intent(in) :: lat_ref(ia_ref,ja_ref) ! latitude [rad] (reference)
1035  real(rp), intent(in) :: hgt_ref(ka_ref,ia_ref,ja_ref) ! height [m] (reference)
1036  real(rp), intent(in) :: lon (ia,ja) ! longitude [rad] (target)
1037  real(rp), intent(in) :: lat (ia,ja) ! latitude [rad] (target)
1038  real(rp), intent(in) :: hgt (ka,ia,ja) ! longitude [m] (target)
1039 
1040  integer, intent(out) :: idx_i(ia,ja,npoints) ! i-index in reference
1041  integer, intent(out) :: idx_j(ia,ja,npoints) ! j-index in reference
1042  real(rp), intent(out) :: hfact(ia,ja,npoints) ! horizontal interp factor
1043  integer, intent(out) :: idx_k(ka,2,ia,ja,npoints) ! k-index in reference
1044  real(rp), intent(out) :: vfact(ka, ia,ja,npoints) ! vertical interp factor
1045 
1046  logical, intent(in), optional :: flag_extrap ! when true, vertical extrapolation will be executed (just copy)
1047 
1048  integer :: nsize, psize, nidx_max
1049  integer, allocatable :: idx_blk(:,:,:), nidx(:,:)
1050  real(rp) :: lon_min, lon_max
1051  real(rp) :: lat_min, lat_max
1052  real(rp) :: dlon, dlat
1053  integer :: idx_ref(npoints)
1054 
1055  integer :: i, j, ii, jj, n
1056  !---------------------------------------------------------------------------
1057 
1058  call prof_rapstart('INTERP_fact',3)
1059 
1060  nsize = ia_ref * ja_ref
1061  if ( nsize > 100 ) then
1062  psize = int( sqrt(2.0_rp*sqrt(real(nsize,rp))) )
1063  nidx_max = nsize / psize * interp_buffer_size_fact
1064  else
1065  psize = 1
1066  nidx_max = nsize
1067  end if
1068 
1069  allocate(idx_blk(nidx_max,psize,psize))
1070  allocate(nidx( psize,psize))
1071 
1072  call interp_div_block(nsize, psize, nidx_max, & ! [IN]
1073  lon_ref(:,:), lat_ref(:,:), & ! [IN]
1074  idx_blk(:,:,:), nidx(:,:), & ! [OUT]
1075  lon_min, lon_max, & ! [OUT]
1076  lat_min, lat_max, & ! [OUT]
1077  dlon, dlat ) ! [OUT]
1078 
1079  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
1080  !$omp private(ii,jj,idx_ref)
1081  do j = 1, ja
1082  do i = 1, ia
1083 
1084  ! main search
1085  call interp_search_horiz( npoints, & ! [IN]
1086  nsize, & ! [IN]
1087  lon_ref(:,:), lat_ref(:,:), & ! [IN]
1088  lon_min, lon_max, & ! [IN]
1089  lat_min, lat_max, & ! [IN]
1090  psize, nidx_max, & ! [IN]
1091  dlon, dlat, & ! [IN]
1092  idx_blk(:,:,:), nidx(:,:), & ! [IN]
1093  lon(i,j), lat(i,j), & ! [IN]
1094  idx_ref(:), hfact(i,j,:) ) ! [OUT]
1095 
1096  do n = 1, npoints
1097  ii = mod(idx_ref(n) - 1, ia_ref) + 1
1098  jj = ( idx_ref(n) - 1 ) / ia_ref + 1
1099  idx_i(i,j,n) = ii
1100  idx_j(i,j,n) = jj
1101  call interp_factor1d( ka_ref, ks_ref, ke_ref, & ! [IN]
1102  ka, ks, ke, & ! [IN]
1103  hgt_ref(:,ii,jj), & ! [IN]
1104  hgt(:,i,j), & ! [IN]
1105  idx_k(:,:,i,j,n), & ! [OUT]
1106  vfact(: ,i,j,n), & ! [OUT]
1107  flag_extrap = flag_extrap ) ! [IN, optional]
1108 
1109  enddo
1110  enddo
1111  enddo
1112 
1113  deallocate(idx_blk, nidx)
1114 
1115  call prof_rapend ('INTERP_fact',3)
1116 
1117  return
1118  end subroutine interp_factor3d_weight
1119 
1120  !-----------------------------------------------------------------------------
1121  ! interpolation for 1D data
1122 !OCL SERIAL
1123  subroutine interp_interp1d( &
1124  KA_ref, KS_ref, KE_ref, &
1125  KA, KS, KE, &
1126  idx_k, &
1127  vfact, &
1128  hgt_ref, &
1129  hgt, &
1130  val_ref, &
1131  val, &
1132  spline, &
1133  logwgt )
1134  use scale_const, only: &
1135  undef => const_undef, &
1136  eps => const_eps
1137  implicit none
1138  integer, intent(in) :: ka_ref, ks_ref, ke_ref ! number of z-direction (reference)
1139  integer, intent(in) :: ka, ks, ke ! number of z-direction (target)
1140 
1141  integer, intent(in) :: idx_k(ka,2) ! k-index in reference
1142  real(rp), intent(in) :: vfact(ka ) ! vertical interp factor
1143  real(rp), intent(in) :: hgt_ref(ka_ref) ! height (reference)
1144  real(rp), intent(in) :: hgt (ka) ! height (target)
1145  real(rp), intent(in), target :: val_ref(ka_ref) ! value (reference)
1146 
1147  real(rp), intent(out) :: val (ka) ! value (target)
1148 
1149  logical, intent(in), optional :: spline
1150  logical, intent(in), optional :: logwgt
1151 
1152  logical :: spline_
1153  logical :: logwgt_
1154 
1155  real(rp), pointer :: work(:)
1156 
1157  integer :: idx (ka_ref)
1158  integer :: idx_r(ka_ref)
1159  real(rp) :: fdz (ka_ref)
1160  real(rp) :: u (ka_ref)
1161 
1162  integer :: kmax
1163  integer :: k
1164  !---------------------------------------------------------------------------
1165 
1166  call prof_rapstart('INTERP_interp',3)
1167 
1168  spline_ = interp_use_spline_vert
1169  if ( present(spline) ) then
1170  spline_ = spline
1171  end if
1172 
1173  logwgt_ = .false.
1174  if ( present(logwgt) ) then
1175  logwgt_ = logwgt
1176  endif
1177 
1178  if ( logwgt_ ) then
1179  allocate( work(ka_ref) )
1180  do k = ks_ref, ke_ref
1181  if ( abs( val_ref(k) - undef ) < eps ) then
1182  work(k) = undef
1183  else
1184  work(k) = log( val_ref(k) )
1185  end if
1186  enddo
1187  else
1188  work => val_ref
1189  endif
1190 
1191  call spline_coef( ka_ref, ks_ref, ke_ref, &
1192  hgt_ref(:), work(:), & ! (in)
1193  spline_, & ! (in)
1194  kmax, & ! (out)
1195  idx(:), idx_r(:), & ! (out)
1196  u(:), fdz(:) ) ! (out)
1197 
1198  call spline_exec( ka_ref, kmax, ka, ks, ke, &
1199  idx_k(:,:), vfact(:), & ! (in)
1200  hgt_ref(:), hgt(:), & ! (in)
1201  work(:), & ! (in)
1202  idx(:), idx_r(:), & ! (in)
1203  u(:), fdz(:), & ! (in)
1204  val(:) ) ! (out)
1205 
1206  if ( logwgt_ ) then
1207  deallocate( work )
1208  do k = ks, ke
1209  if ( abs( val(k) - undef ) > eps ) then
1210  val(k) = exp( val(k) )
1211  endif
1212  end do
1213  endif
1214 
1215  call prof_rapend ('INTERP_interp',3)
1216 
1217  return
1218  end subroutine interp_interp1d
1219 
1220  !-----------------------------------------------------------------------------
1221  ! interpolation for 2D data (nearest-neighbor)
1222  subroutine interp_interp2d( &
1223  npoints, &
1224  IA_ref, JA_ref, &
1225  IA, JA, &
1226  idx_i, idx_j, &
1227  hfact, &
1228  val_ref, &
1229  val, &
1230  threshold_undef, &
1231  wsum, val2 )
1232  use scale_const, only: &
1233  undef => const_undef, &
1234  eps => const_eps
1235  implicit none
1236 
1237  integer, intent(in) :: npoints ! number of interpolation point for horizontal
1238  integer, intent(in) :: ia_ref ! number of x-direction (reference)
1239  integer, intent(in) :: ja_ref ! number of y-direction (reference)
1240  integer, intent(in) :: ia ! number of x-direction (target)
1241  integer, intent(in) :: ja ! number of y-direction (target)
1242  integer, intent(in) :: idx_i (ia,ja,npoints) ! i-index in reference
1243  integer, intent(in) :: idx_j (ia,ja,npoints) ! j-index in reference
1244  real(rp), intent(in) :: hfact (ia,ja,npoints) ! horizontal interp factor
1245  real(rp), intent(in) :: val_ref(ia_ref,ja_ref) ! value (reference)
1246 
1247  real(rp), intent(out) :: val (ia,ja) ! value (target)
1248 
1249  real(rp), intent(in), optional :: threshold_undef
1250  real(rp), intent(out), optional :: wsum(ia,ja)
1251  real(rp), intent(out), optional :: val2(ia,ja)
1252 
1253  real(rp) :: th_undef
1254 
1255  real(rp) :: fact, valn, f, w
1256  real(rp) :: sw
1257  logical :: lval2
1258 
1259  integer :: i, j, n
1260  !---------------------------------------------------------------------------
1261 
1262  call prof_rapstart('INTERP_interp',3)
1263 
1264  if ( present(threshold_undef) ) then
1265  th_undef = threshold_undef
1266  else
1267  th_undef = interp_threshold_undef
1268  end if
1269  th_undef = min( max( th_undef, eps * 2.0_rp ), 1.0_rp - eps * 2.0_rp )
1270 
1271  lval2 = present(wsum) .and. present(val2)
1272 
1273  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
1274  !$omp private(fact,valn,f,w,sw)
1275 !OCL PREFETCH
1276  do j = 1, ja
1277  do i = 1, ia
1278  fact = 0.0_rp
1279  valn = 0.0_rp
1280  do n = 1, npoints
1281  f = hfact(i,j,n)
1282  w = val_ref(idx_i(i,j,n),idx_j(i,j,n))
1283  if ( f > eps .and. abs( w - undef ) > eps ) then
1284  fact = fact + f
1285  valn = valn + f * w
1286  end if
1287  end do
1288  sw = 0.5_rp - sign( 0.5_rp, fact - th_undef + eps ) ! 1.0 when fact < threshold
1289  val(i,j) = valn / ( fact + sw ) * ( 1.0_rp - sw ) + undef * sw
1290  if ( lval2 ) then
1291  wsum(i,j) = fact
1292  sw = 0.5_rp - sign( 0.5_rp, fact - eps ) ! 1.0 when fact < 0.0
1293  val2(i,j) = valn / ( fact + sw ) * ( 1.0_rp - sw ) + undef * sw
1294  end if
1295  enddo
1296  enddo
1297 
1298  call prof_rapend ('INTERP_interp',3)
1299 
1300  return
1301  end subroutine interp_interp2d
1302 
1303  !-----------------------------------------------------------------------------
1304  ! interpolation using one-points for 3D data (nearest-neighbor)
1305  subroutine interp_interp3d( &
1306  npoints, &
1307  KA_ref, KS_ref, KE_ref, &
1308  IA_ref, JA_ref, &
1309  KA, KS, KE, &
1310  IA, JA, &
1311  idx_i, idx_j, &
1312  hfact, &
1313  idx_k, &
1314  vfact, &
1315  hgt_ref, &
1316  hgt, &
1317  val_ref, &
1318  val, &
1319  spline, logwgt, &
1320  threshold_undef, &
1321  wsum, val2 )
1322  use scale_const, only: &
1323  undef => const_undef, &
1324  eps => const_eps
1325  implicit none
1326  integer, intent(in) :: npoints ! number of interpolation point for horizontal
1327  integer, intent(in) :: ka_ref, ks_ref, ke_ref ! number of z-direction (reference)
1328  integer, intent(in) :: ia_ref ! number of x-direction (reference)
1329  integer, intent(in) :: ja_ref ! number of y-direction (reference)
1330  integer, intent(in) :: ka, ks, ke ! number of z-direction (target)
1331  integer, intent(in) :: ia ! number of x-direction (target)
1332  integer, intent(in) :: ja ! number of y-direction (target)
1333 
1334  integer, intent(in) :: idx_i (ia,ja,npoints) ! i-index in reference
1335  integer, intent(in) :: idx_j (ia,ja,npoints) ! j-index in reference
1336  real(rp), intent(in) :: hfact (ia,ja,npoints) ! horizontal interp factor
1337  integer, intent(in) :: idx_k (ka,2,ia,ja,npoints) ! k-index in reference
1338  real(rp), intent(in) :: vfact (ka, ia,ja,npoints) ! vertical interp factor
1339  real(rp), intent(in) :: hgt_ref(ka_ref,ia_ref,ja_ref) ! height (reference)
1340  real(rp), intent(in) :: hgt (ka,ia,ja) ! height (target)
1341  real(rp), intent(in), target :: val_ref(ka_ref,ia_ref,ja_ref) ! value (reference)
1342 
1343  real(rp), intent(out) :: val (ka,ia,ja) ! value (target)
1344 
1345  logical, intent(in), optional :: spline
1346  logical, intent(in), optional :: logwgt
1347  real(rp), intent(in), optional :: threshold_undef
1348  real(rp), intent(out), optional :: wsum(ka,ia,ja)
1349  real(rp), intent(out), optional :: val2(ka,ia,ja)
1350 
1351  real(rp) :: th_undef
1352  logical :: spline_
1353  logical :: logwgt_
1354 
1355  real(rp), pointer :: work(:,:,:)
1356 
1357  integer, allocatable :: kmax(:,:)
1358  integer, allocatable :: idx(:,:,:), idx_r(:,:,:)
1359  real(rp), allocatable :: u(:,:,:), fdz(:,:,:)
1360 
1361  real(rp) :: valn
1362  real(rp) :: fact
1363  real(rp) :: f
1364  real(rp) :: sw
1365  real(rp) :: w(ka,npoints)
1366  logical :: lval2
1367 
1368  integer :: imin, imax
1369  integer :: jmin, jmax
1370  integer :: ii, jj
1371  integer :: k, i, j, n
1372  !---------------------------------------------------------------------------
1373 
1374  call prof_rapstart('INTERP_interp',3)
1375 
1376  if ( present(threshold_undef) ) then
1377  th_undef = threshold_undef
1378  else
1379  th_undef = interp_threshold_undef
1380  end if
1381  th_undef = min( max( th_undef, eps * 2.0_rp ), 1.0_rp - eps * 2.0_rp )
1382 
1383 
1384  spline_ = interp_use_spline_vert
1385  if ( present(spline) ) then
1386  spline_ = spline
1387  end if
1388 
1389  logwgt_ = .false.
1390  if ( present(logwgt) ) then
1391  logwgt_ = logwgt
1392  endif
1393 
1394 
1395  lval2 = present(wsum) .and. present(val2)
1396 
1397 
1398  imin = ia_ref
1399  jmin = ja_ref
1400  imax = 1
1401  jmax = 1
1402  !$omp parallel do OMP_SCHEDULE_ collapse(3) &
1403  !$omp reduction(min: imin,jmin) &
1404  !$omp reduction(max: imax,jmax)
1405  do n = 1, npoints
1406  do j = 1, ja
1407  do i = 1, ia
1408  imin = min(imin, idx_i(i,j,n))
1409  imax = max(imax, idx_i(i,j,n))
1410  jmin = min(jmin, idx_j(i,j,n))
1411  jmax = max(jmax, idx_j(i,j,n))
1412  end do
1413  end do
1414  end do
1415 
1416  allocate( kmax(imin:imax,jmin:jmax) )
1417  allocate( idx(ka_ref,imin:imax,jmin:jmax), idx_r(ka_ref,imin:imax,jmin:jmax) )
1418  allocate( u(ka_ref,imin:imax,jmin:jmax), fdz(ka_ref,imin:imax,jmin:jmax) )
1419 
1420  if ( logwgt_ ) then
1421  allocate( work(ka_ref,imin:imax,jmin:jmax) )
1422  !$omp parallel do OMP_SCHEDULE_ collapse(2)
1423  do j = jmin, jmax
1424  do i = imin, imax
1425  do k = ks_ref, ke_ref
1426  if ( abs( val_ref(k,i,j) - undef ) < eps ) then
1427  work(k,i,j) = undef
1428  else
1429  work(k,i,j) = log( val_ref(k,i,j) )
1430  end if
1431  enddo
1432  enddo
1433  enddo
1434  else
1435  work => val_ref
1436  endif
1437 
1438  !$omp parallel do OMP_SCHEDULE_ collapse(2)
1439  do j = jmin, jmax
1440  do i = imin, imax
1441  call spline_coef( ka_ref, ks_ref, ke_ref, &
1442  hgt_ref(:,i,j), work(:,i,j), & ! (in)
1443  spline_, & ! (in)
1444  kmax(i,j), & ! (out)
1445  idx(:,i,j), idx_r(:,i,j), & ! (out)
1446  u(:,i,j), fdz(:,i,j) ) ! (out)
1447  end do
1448  end do
1449 
1450  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
1451  !$omp private(valn,fact,w,f,ii,jj,sw)
1452  do j = 1, ja
1453  do i = 1, ia
1454 
1455  do n = 1, npoints
1456  if ( hfact(i,j,n) < eps ) exit
1457  ii = idx_i(i,j,n)
1458  jj = idx_j(i,j,n)
1459  call spline_exec( ka_ref, kmax(ii,jj), ka, ks, ke, &
1460  idx_k(:,:,i,j,n), vfact(:, i,j,n), & ! (in)
1461  hgt_ref(:,ii,jj), hgt(:,i,j), & ! (in)
1462  work(:,ii,jj), & ! (in)
1463  idx(:,ii,jj), idx_r(:,ii,jj), & ! (in)
1464  u(:,ii,jj), fdz(:,ii,jj), & ! (in)
1465  w(:,n) ) ! (out)
1466  end do
1467 
1468  do k = ks, ke
1469  fact = 0.0_rp
1470  valn = 0.0_rp
1471  do n = 1, npoints
1472  f = hfact(i,j,n)
1473  if ( f > eps .and. abs( w(k,n) - undef ) > eps ) then
1474  fact = fact + f
1475  valn = valn + f * w(k,n)
1476  endif
1477  enddo
1478  sw = 0.5_rp - sign( 0.5_rp, fact - th_undef ) ! 1.0 when fact < threshold
1479  val(k,i,j) = valn / ( fact + sw ) * ( 1.0_rp - sw ) + undef * sw
1480  if ( lval2 ) then
1481  wsum(k,i,j) = fact
1482  sw = 0.5_rp - sign( 0.5_rp, fact - eps ) ! 1.0 when fact < 0.0
1483  val2(k,i,j) = valn / ( fact + sw ) * ( 1.0_rp - sw ) + undef * sw
1484  end if
1485  enddo
1486 
1487  enddo
1488  enddo
1489 
1490  deallocate( kmax, idx, idx_r )
1491  deallocate( u, fdz )
1492 
1493  if ( logwgt_ ) then
1494  deallocate( work )
1495  !$omp parallel do OMP_SCHEDULE_ collapse(2)
1496  do j = 1, ja
1497  do i = 1, ia
1498  do k = ks, ke
1499  if ( abs( val(k,i,j) - undef ) > eps ) then
1500  val(k,i,j) = exp( val(k,i,j) )
1501  end if
1502  end do
1503  end do
1504  end do
1505  end if
1506 
1507  call prof_rapend ('INTERP_interp',3)
1508 
1509  return
1510  end subroutine interp_interp3d
1511 
1512 
1513  ! private
1514 
1515 
1516  !-----------------------------------------------------------------------------
1517  ! horizontal search of interpolation points
1518 !OCL SERIAL
1519  subroutine interp_search_horiz( &
1520  npoints, &
1521  nsize, &
1522  lon_ref, lat_ref, &
1523  lon_min, lon_max, &
1524  lat_min, lat_max, &
1525  psize, nidx_max, &
1526  dlon, dlat, &
1527  idx_blk, nidx, &
1528  lon, lat, &
1529  idx_ref, &
1530  hfact, &
1531  search_limit, &
1532  weight_order )
1533  use scale_const, only: &
1534  eps => const_eps, &
1535  radius => const_radius
1536  implicit none
1537 
1538  integer, intent(in) :: npoints ! number of interpolation point for horizontal
1539  integer, intent(in) :: nsize ! number of grids (reference)
1540  real(rp), intent(in) :: lon_ref(nsize) ! longitude [rad] (reference)
1541  real(rp), intent(in) :: lat_ref(nsize) ! latitude [rad] (reference)
1542  real(rp), intent(in) :: lon_min ! minmum longitude [rad] (reference)
1543  real(rp), intent(in) :: lon_max ! maximum longitude [rad] (reference)
1544  real(rp), intent(in) :: lat_min ! minmum latitude [rad] (reference)
1545  real(rp), intent(in) :: lat_max ! maximum latitude [rad] (reference)
1546  integer, intent(in) :: psize ! number of blocks for each dimension
1547  integer, intent(in) :: nidx_max ! maximum number of index in the block
1548  real(rp), intent(in) :: dlon ! block longitude difference
1549  real(rp), intent(in) :: dlat ! block latitude difference
1550  integer, intent(in) :: idx_blk(nidx_max,psize,psize) ! index of the reference in the block
1551  integer, intent(in) :: nidx (psize,psize) ! number of indexes in the block
1552  real(rp), intent(in) :: lon ! longitude [rad] (target)
1553  real(rp), intent(in) :: lat ! latitude [rad] (target)
1554  integer, intent(out) :: idx_ref(npoints) ! index in reference (target)
1555  real(rp), intent(out) :: hfact(npoints) ! horizontal interp factor (target)
1556 
1557  real(rp), intent(in), optional :: search_limit
1558  integer, intent(in), optional :: weight_order
1559 
1560  real(rp) :: drad(npoints)
1561  real(rp) :: sum
1562  real(rp) :: search_limit_
1563  integer :: weight_order_
1564 
1565  real(rp) :: lon0, lon1, lat0, lat1
1566  real(rp) :: dlon_sl, dlat_sl
1567  integer :: i, n
1568  integer :: ii, jj, ii0, jj0
1569  !---------------------------------------------------------------------------
1570 
1571  if ( present(search_limit) ) then
1572  search_limit_ = search_limit
1573  else
1574  search_limit_ = interp_search_limit
1575  end if
1576  search_limit_ = search_limit_ / radius ! m to radian
1577 
1578  if ( present(weight_order) ) then
1579  weight_order_ = weight_order
1580  else
1581  weight_order_ = interp_weight_order
1582  end if
1583 
1584 
1585 
1586  drad(:) = large_number
1587  idx_ref(:) = -1
1588 
1589  ! find k-nearest points in the nearest block
1590  ii0 = max( min( int(( lon - lon_min ) / dlon) + 1, psize ), 1 )
1591  jj0 = max( min( int(( lat - lat_min ) / dlat) + 1, psize ), 1 )
1592  do i = 1, nidx(ii0,jj0)
1593  n = idx_blk(i,ii0,jj0)
1594  call interp_insert( npoints, &
1595  lon, lat, &
1596  lon_ref(n), lat_ref(n), & ! [IN]
1597  n, & ! [IN]
1598  drad(:), idx_ref(:) ) ! [INOUT]
1599  end do
1600 
1601  if ( abs(drad(1)) < eps ) then
1602  hfact(:) = 0.0_rp
1603  hfact(1) = 1.0_rp
1604 
1605  return
1606  else if ( drad(1) > search_limit_ ) then
1607  hfact(:) = 0.0_rp
1608  idx_ref(:) = 1 ! dummy
1609 
1610  return
1611  end if
1612 
1613  dlon_sl = max(dlon * 0.5_rp, drad(npoints))
1614  dlat_sl = max(dlat * 0.5_rp, drad(npoints))
1615  do jj = 1, psize
1616  lat0 = lat_min + dlat * (jj-1)
1617  lat1 = lat_min + dlat * jj
1618  if ( lat < lat0 - dlat_sl &
1619  .or. lat >= lat1 + dlat_sl ) cycle
1620  do ii = 1, psize
1621  if ( ii==ii0 .and. jj==jj0 ) cycle
1622  lon0 = lon_min + dlon * (ii-1)
1623  lon1 = lon_min + dlon * ii
1624  if ( lon < lon0 - dlon_sl &
1625  .or. lon >= lon1 + dlon_sl ) cycle
1626  do i = 1, nidx(ii,jj)
1627  n = idx_blk(i,ii,jj)
1628  call interp_insert( npoints, &
1629  lon, lat, &
1630  lon_ref(n), lat_ref(n), & ! [IN]
1631  n, & ! [IN]
1632  drad(:), idx_ref(:) ) ! [INOUT]
1633  end do
1634  dlon_sl = max(dlon * 0.5_rp, drad(npoints))
1635  dlat_sl = max(dlat * 0.5_rp, drad(npoints))
1636  end do
1637  end do
1638 
1639  ! factor = 1 / dradian
1640  if ( weight_order_ < 0 ) then
1641  do n = 1, npoints
1642  hfact(n) = 1.0_rp / drad(n)**(-1.0_rp/weight_order_)
1643  enddo
1644  else
1645  do n = 1, npoints
1646  hfact(n) = 1.0_rp / drad(n)**weight_order_
1647  enddo
1648  end if
1649 
1650  ! ignore far point
1651  do n = 1, npoints
1652  if ( drad(n) >= search_limit_ ) then
1653  hfact(n) = 0.0_rp
1654  endif
1655  enddo
1656 
1657  ! normalize factor
1658  sum = 0.0_rp
1659  do n = 1, npoints
1660  sum = sum + hfact(n)
1661  enddo
1662 
1663  if ( sum > 0.0_rp ) then
1664  do n = 1, npoints
1665  hfact(n) = hfact(n) / sum
1666  enddo
1667  endif
1668 
1669  return
1670  end subroutine interp_search_horiz
1671 
1672  !-----------------------------------------------------------------------------
1673  ! horizontal search of interpolation points for structure grid
1674 !OCL SERIAL
1675  subroutine interp_search_horiz_struct( &
1676  npoints, &
1677  psizex, psizey, &
1678  IA_ref, JA_ref, &
1679  lon_ref, lat_ref, &
1680  lon_min, lat_min, &
1681  dlon, dlat, &
1682  i0, i1, j0, j1, &
1683  lon, lat, &
1684  idx_i, idx_j, &
1685  hfact, &
1686  search_limit, &
1687  weight_order )
1688  use scale_const, only: &
1689  eps => const_eps, &
1690  radius => const_radius
1691  implicit none
1692 
1693  integer, intent(in) :: npoints ! number of interpolation point for horizontal
1694  integer, intent(in) :: psizex ! number of block in x-direction (reference)
1695  integer, intent(in) :: psizey ! number of block in y-direction (reference)
1696  integer, intent(in) :: IA_ref ! number of grids in x-direction (reference)
1697  integer, intent(in) :: JA_ref ! number of grids in y-direction (reference)
1698  real(RP), intent(in) :: lon_ref(IA_ref) ! longitude [rad] (reference)
1699  real(RP), intent(in) :: lat_ref(JA_ref) ! latitude [rad] (reference)
1700  real(RP), intent(in) :: lon_min ! minimum longitude
1701  real(RP), intent(in) :: lat_min ! minimum latitude
1702  real(RP), intent(in) :: dlon
1703  real(RP), intent(in) :: dlat
1704  integer, intent(in) :: i0(psizex) ! start index in the block
1705  integer, intent(in) :: i1(psizex) ! end index in the block
1706  integer, intent(in) :: j0(psizey) ! start index in the block
1707  integer, intent(in) :: j1(psizey) ! start index in the block
1708  real(RP), intent(in) :: lon ! longitude [rad] (target)
1709  real(RP), intent(in) :: lat ! latitude [rad] (target)
1710  integer, intent(out) :: idx_i(npoints) ! x-index in reference (target)
1711  integer, intent(out) :: idx_j(npoints) ! y-index in reference (target)
1712  real(RP), intent(out) :: hfact(npoints) ! horizontal interp factor (target)
1713 
1714  real(RP), intent(in), optional :: search_limit
1715  integer, intent(in), optional :: weight_order
1716 
1717  real(RP) :: drad(npoints)
1718  real(RP) :: sum
1719  real(RP) :: search_limit_
1720  integer :: weight_order_
1721 
1722  real(RP) :: lon0, lon1, lat0, lat1
1723  real(RP) :: dlon_sl, dlat_sl
1724 
1725  integer :: i, j, n
1726  integer :: ii, jj, ii0, jj0
1727  !---------------------------------------------------------------------------
1728 
1729  if ( present(search_limit) ) then
1730  search_limit_ = search_limit
1731  else
1732  search_limit_ = interp_search_limit
1733  end if
1734  search_limit_ = search_limit_ / radius ! m to radian
1735 
1736  if ( present(weight_order) ) then
1737  weight_order_ = weight_order
1738  else
1739  weight_order_ = interp_weight_order
1740  end if
1741 
1742 
1743  drad(:) = large_number
1744  idx_i(:) = 1
1745  idx_j(:) = 1
1746 
1747  ! find k-nearest points in the nearest block
1748  ii0 = max( min( int(( lon - lon_min ) / dlon) + 1, psizex ), 1)
1749  jj0 = max( min( int(( lat - lat_min ) / dlat) + 1, psizey ), 1)
1750  do j = j0(jj0), j1(jj0)
1751  do i = i0(ii0), i1(ii0)
1752  call interp_insert_2d( npoints, &
1753  lon, lat, &
1754  lon_ref(i), lat_ref(j), & ! [IN]
1755  i, j, & ! [IN]
1756  drad(:), idx_i(:), idx_j(:) ) ! [INOUT]
1757  end do
1758  end do
1759 
1760  if ( abs(drad(1)) < eps ) then
1761  hfact(:) = 0.0_rp
1762  hfact(1) = 1.0_rp
1763 
1764  return
1765  else if ( drad(1) > search_limit_ ) then
1766  hfact(:) = 0.0_rp
1767  idx_i(:) = 1 ! dummy
1768  idx_j(:) = 1 ! dummy
1769 
1770  return
1771  end if
1772 
1773  dlon_sl = max(dlon * 0.5_rp, drad(npoints))
1774  dlat_sl = max(dlat * 0.5_rp, drad(npoints))
1775  do jj = 1, psizey
1776  lat0 = lat_min + dlat * (jj-1)
1777  lat1 = lat_min + dlat * jj
1778  if ( lat < lat0 - dlat_sl &
1779  .or. lat >= lat1 + dlat_sl ) cycle
1780  do ii = 1, psizex
1781  if ( ii==ii0 .and. jj==jj0 ) cycle
1782  lon0 = lon_min + dlon * (ii-1)
1783  lon1 = lon_min + dlon * ii
1784  if ( lon < lon0 - dlon_sl &
1785  .or. lon >= lon1 + dlon_sl ) cycle
1786  do j = j0(jj0), j1(jj0)
1787  do i = i0(ii0), i1(ii0)
1788  call interp_insert_2d( npoints, &
1789  lon, lat, &
1790  lon_ref(i), lat_ref(j), & ! [IN]
1791  i, j, & ! [IN]
1792  drad(:), idx_i(:), idx_j(:) ) ! [INOUT]
1793  end do
1794  end do
1795  dlon_sl = max(dlon * 0.5_rp, drad(npoints))
1796  dlat_sl = max(dlat * 0.5_rp, drad(npoints))
1797  end do
1798  end do
1799 
1800 
1801  ! factor = 1 / drad
1802  if ( weight_order_ < 0 ) then
1803  do n = 1, npoints
1804  hfact(n) = 1.0_rp / drad(n)**(-1.0_rp/weight_order_)
1805  enddo
1806  else
1807  do n = 1, npoints
1808  hfact(n) = 1.0_rp / drad(n)**weight_order_
1809  enddo
1810  end if
1811 
1812  ! ignore far point
1813  do n = 1, npoints
1814  if ( drad(n) >= search_limit_ ) then
1815  hfact(n) = 0.0_rp
1816  endif
1817  enddo
1818 
1819  ! normalize factor
1820  sum = 0.0_rp
1821  do n = 1, npoints
1822  sum = sum + hfact(n)
1823  enddo
1824 
1825  if ( sum > 0.0_rp ) then
1826  do n = 1, npoints
1827  hfact(n) = hfact(n) / sum
1828  enddo
1829  endif
1830 
1831  return
1832  end subroutine interp_search_horiz_struct
1833 
1834 
1835 !OCL SERIAL
1836  subroutine interp_insert( npoints, &
1837  lon, lat, &
1838  lon_ref, lat_ref, &
1839  i, &
1840  drad, idx_i )
1841  use scale_const, only: &
1842  undef => const_undef, &
1843  eps => const_eps
1844  use scale_sort, only: &
1845  sort_exec
1846 
1847  integer, intent(in) :: npoints
1848  real(RP), intent(in) :: lon, lat
1849  real(RP), intent(in) :: lon_ref, lat_ref
1850  integer, intent(in) :: i
1851 
1852  real(RP), intent(inout) :: drad(npoints)
1853  integer, intent(inout) :: idx_i(npoints)
1854 
1855  real(RP) :: dradian
1856 
1857  if ( abs( lon_ref - undef ) < eps ) return
1858 
1859  dradian = haversine( lon, lat, lon_ref, lat_ref )
1860 
1861  if ( dradian <= drad(npoints) ) then
1862  ! replace last(=longest) value
1863  drad(npoints) = dradian
1864  idx_i(npoints) = i
1865 
1866  ! sort by ascending order
1867  call sort_exec( npoints, & ! [IN]
1868  drad(:), & ! [INOUT]
1869  idx_i(:) ) ! [INOUT]
1870  endif
1871 
1872  return
1873  end subroutine interp_insert
1874 
1875 !OCL SERIAL
1876  subroutine interp_insert_2d( npoints, &
1877  lon, lat, &
1878  lon_ref, lat_ref, &
1879  i, j, &
1880  drad, idx_i, idx_j )
1881  use scale_const, only: &
1882  undef => const_undef, &
1883  eps => const_eps
1884  use scale_sort, only: &
1885  sort_exec
1886 
1887  integer, intent(in) :: npoints
1888  real(RP), intent(in) :: lon, lat
1889  real(RP), intent(in) :: lon_ref, lat_ref
1890  integer, intent(in) :: i, j
1891 
1892  real(RP), intent(inout) :: drad(npoints)
1893  integer, intent(inout) :: idx_i(npoints)
1894  integer, intent(inout) :: idx_j(npoints)
1895 
1896  real(RP) :: dradian
1897 
1898  if ( abs( lon_ref - undef ) < eps ) return
1899 
1900  dradian = haversine( lon, lat, lon_ref, lat_ref )
1901 
1902  if ( dradian <= drad(npoints) ) then
1903  ! replace last(=longest) value
1904  drad(npoints) = dradian
1905  idx_i(npoints) = i
1906  idx_j(npoints) = j
1907 
1908  ! sort by ascending order
1909  call sort_exec( npoints, & ! [IN]
1910  drad(:), & ! [INOUT]
1911  idx_i(:), & ! [INOUT]
1912  idx_j(:) ) ! [INOUT]
1913  endif
1914 
1915  return
1916  end subroutine interp_insert_2d
1917 
1918 
1919  subroutine interp_div_block(nsize, psize, nidx_max, &
1920  lon_ref, lat_ref, &
1921  idx, nidx, &
1922  lon_min, lon_max, &
1923  lat_min, lat_max, &
1924  dlon, dlat )
1925  use scale_const, only: &
1926  undef => const_undef, &
1927  eps => const_eps
1928  use scale_prc, only: &
1929  prc_abort
1930  integer, intent(in) :: nsize
1931  integer, intent(in) :: psize
1932  integer, intent(in) :: nidx_max
1933  real(RP), intent(in) :: lon_ref(nsize)
1934  real(RP), intent(in) :: lat_ref(nsize)
1935 
1936  integer, intent(out) :: idx (nidx_max,psize,psize)
1937  integer, intent(out) :: nidx(psize,psize)
1938  real(RP), intent(out) :: lon_min, lon_max
1939  real(RP), intent(out) :: lat_min, lat_max
1940  real(RP), intent(out) :: dlon, dlat
1941 
1942  integer :: i, ii, jj, n
1943 
1944  lon_min = minval(lon_ref(:), mask=abs(lon_ref-undef)>eps)
1945  lon_max = maxval(lon_ref(:), mask=abs(lon_ref-undef)>eps)
1946  lat_min = minval(lat_ref(:), mask=abs(lat_ref-undef)>eps)
1947  lat_max = maxval(lat_ref(:), mask=abs(lat_ref-undef)>eps)
1948 
1949  dlon = ( lon_max - lon_min ) / psize
1950  dlat = ( lat_max - lat_min ) / psize
1951 
1952  nidx(:,:) = 0
1953  do i = 1, nsize
1954  if ( abs( lon_ref(i) - undef ) < eps ) cycle
1955  ii = min(int((lon_ref(i) - lon_min) / dlon) + 1, psize)
1956  jj = min(int((lat_ref(i) - lat_min) / dlat) + 1, psize)
1957  n = nidx(ii,jj) + 1
1958  nidx(ii,jj) = n
1959  if ( n <= nidx_max ) idx(n,ii,jj) = i
1960  end do
1961 
1962  if ( maxval(nidx) > nidx_max ) then
1963  log_error("INTERP_search_horiz",*) 'Buffer size is not enough'
1964  log_error_cont(*) ' Use larger INTERP_buffer_size_fact: ', interp_buffer_size_fact * maxval(nidx)/nidx_max
1965  call prc_abort
1966  end if
1967 
1968  return
1969  end subroutine interp_div_block
1970 
1971  !-----------------------------------------------------------------------------
1972  subroutine interp_bilinear_inv( &
1973  x_ref0, x_ref1, x_ref2, x_ref3, &
1974  y_ref0, y_ref1, y_ref2, y_ref3, &
1975  x, y, &
1976  u, v, &
1977  error )
1978  implicit none
1979  real(RP), intent(in) :: x_ref0, x_ref1, x_ref2, x_ref3
1980  real(RP), intent(in) :: y_ref0, y_ref1, y_ref2, y_ref3
1981  real(RP), intent(in) :: x, y
1982 
1983  real(RP), intent(out) :: u, v
1984  logical, intent(out) :: error
1985 
1986  real(RP) :: e_x, e_y
1987  real(RP) :: f_x, f_y
1988  real(RP) :: g_x, g_y
1989  real(RP) :: h_x, h_y
1990  real(RP) :: k0, k1, k2
1991  real(RP) :: w
1992  real(RP) :: sig
1993 
1994  e_x = x_ref1 - x_ref0
1995  e_y = y_ref1 - y_ref0
1996 
1997  f_x = x_ref3 - x_ref0
1998  f_y = y_ref3 - y_ref0
1999 
2000  g_x = x_ref0 - x_ref1 + x_ref2 - x_ref3
2001  g_y = y_ref0 - y_ref1 + y_ref2 - y_ref3
2002 
2003  h_x = x - x_ref0
2004  h_y = y - y_ref0
2005 
2006  k0 = cross(h_x, h_y, e_x, e_y)
2007  k1 = cross(e_x, e_y, f_x, f_y) + cross(h_x, h_y, g_x, g_y)
2008  k2 = cross(g_x, g_y, f_x, f_y)
2009 
2010  w = k1**2 - 4.0_rp * k0 * k2
2011  if ( w < 0.0_rp ) then
2012  log_error("INTERP_bilinear_inv",*) 'Outside of the quadrilateral'
2013  error = .true.
2014  return
2015  end if
2016 
2017  if ( abs(k1) < eps_bilinear ) then
2018  log_error("INTERP_bilinear_inv",*) 'Unexpected error occured', k1
2019  log_error_cont(*) x_ref0, x_ref1, x_ref2, x_ref3
2020  log_error_cont(*) y_ref0, y_ref1, y_ref2, y_ref3
2021  log_error_cont(*) x, y
2022  error = .true.
2023  return
2024  end if
2025  if ( abs(k2) < eps_bilinear * sqrt( (x_ref2-x_ref0)**2+(y_ref2-y_ref0)**2 ) ) then
2026  v = - k0 / k1
2027  else
2028  sig = sign( 1.0_rp, cross(x_ref1-x_ref0, y_ref1-y_ref0, x_ref3-x_ref0, y_ref3-y_ref0) )
2029  v = ( - k1 + sig * sqrt(w) ) / ( 2.0_rp * k2 )
2030  end if
2031  u = ( h_x - f_x * v ) / ( e_x + g_x * v )
2032 
2033  if ( u < -eps_bilinear .or. u > 1.0_rp+eps_bilinear .or. v < -eps_bilinear .or. v > 1.0_rp+eps_bilinear ) then
2034  log_error("INTERP_bilinear_inv",*) 'Unexpected error occured', u, v
2035  log_error_cont(*) x_ref0, x_ref1, x_ref2, x_ref3
2036  log_error_cont(*) y_ref0, y_ref1, y_ref2, y_ref3
2037  log_error_cont(*) x, y
2038  log_error_cont(*) k0, k1, k2, sig
2039  error = .true.
2040  return
2041  end if
2042  u = min( max( u, 0.0_rp ), 1.0_rp )
2043  v = min( max( v, 0.0_rp ), 1.0_rp )
2044 
2045  error = .false.
2046 
2047  return
2048  end subroutine interp_bilinear_inv
2049 
2050  !-----------------------------------------------------------------------------
2051  ! check whether the point is inside of the quadrilateral
2052  subroutine interp_check_inside( &
2053  x_ref0, x_ref1, x_ref2, x_ref3, &
2054  y_ref0, y_ref1, y_ref2, y_ref3, &
2055  x, y, &
2056  inc_i, inc_j, &
2057  error )
2058  use scale_const, only: &
2059  eps => const_eps
2060  implicit none
2061  real(RP), intent(in) :: x_ref0, x_ref1, x_ref2, x_ref3
2062  real(RP), intent(in) :: y_ref0, y_ref1, y_ref2, y_ref3
2063  real(RP), intent(in) :: x, y
2064 
2065  integer, intent(out) :: inc_i, inc_j
2066  logical, intent(out) :: error
2067 
2068  real(RP) :: sig
2069  real(RP) :: c1, c2, c3, c4
2070  real(RP) :: th
2071  logical :: fx, fy
2072 
2073  error = .false.
2074 
2075  sig = sign( 1.0_rp, cross(x_ref1-x_ref0, y_ref1-y_ref0, x_ref3-x_ref0, y_ref3-y_ref0) )
2076 
2077  c1 = sig * cross(x_ref1-x_ref0, y_ref1-y_ref0, x-x_ref0, y-y_ref0)
2078  c2 = sig * cross(x_ref2-x_ref1, y_ref2-y_ref1, x-x_ref1, y-y_ref1)
2079  c3 = sig * cross(x_ref3-x_ref2, y_ref3-y_ref2, x-x_ref2, y-y_ref2)
2080  c4 = sig * cross(x_ref0-x_ref3, y_ref0-y_ref3, x-x_ref3, y-y_ref3)
2081 
2082  th = - max( abs(c1), abs(c2), abs(c3), abs(c4) ) * eps * 1.e4_rp
2083 
2084  ! if all the c1 - c4 are positive, the point is inside the quadrilateral
2085  inc_i = 0
2086  inc_j = 0
2087  if ( c1 < th ) inc_j = -1
2088  if ( c2 < th ) inc_i = 1
2089  if ( c3 < th ) inc_j = 1
2090  if ( c4 < th ) inc_i = -1
2091 
2092  fx = c2 < th .and. c4 < th
2093  fy = c1 < th .and. c3 < th
2094  if ( fx .and. fy ) then
2095  log_error("INTERP_check_inside",*) 'Unexpected error occured', c1, c2, c3, c4
2096  log_error_cont(*) x_ref0, x_ref1, x_ref2, x_ref3
2097  log_error_cont(*) y_ref0, y_ref1, y_ref2, y_ref3
2098  log_error_cont(*) x, y
2099  error = .true.
2100  return
2101  else if ( fx ) then
2102  inc_i = 0
2103  else if ( fy ) then
2104  inc_j = 0
2105  end if
2106 
2107  return
2108  end subroutine interp_check_inside
2109 
2110  !-----------------------------------------------------------------------------
2112  function cross( &
2113  x0, y0, &
2114  x1, y1 )
2115  implicit none
2116  real(rp), intent(in) :: x0, y0
2117  real(rp), intent(in) :: x1, y1
2118  real(rp) :: cross
2119 
2120  cross = x0 * y1 - x1 * y0
2121 
2122  end function cross
2123 
2124  !-----------------------------------------------------------------------------
2125  ! Haversine Formula (from R.W. Sinnott, "Virtues of the Haversine",
2126  ! Sky and Telescope, vol. 68, no. 2, 1984, p. 159):
2127  function haversine( &
2128  lon0, lat0, &
2129  lon1, lat1 ) &
2130  result( d )
2131  implicit none
2132 
2133  real(rp), intent(in) :: lon0 ! [rad]
2134  real(rp), intent(in) :: lat0 ! [rad]
2135  real(rp), intent(in) :: lon1 ! [rad]
2136  real(rp), intent(in) :: lat1 ! [rad]
2137  real(rp) :: d ! [rad]
2138 
2139  real(rp) :: dlonh, dlath
2140  real(rp) :: work1
2141  !---------------------------------------------------------------------------
2142 
2143  dlonh = 0.5_rp * ( lon0 - lon1 )
2144  dlath = 0.5_rp * ( lat0 - lat1 )
2145 
2146  work1 = sin(dlath)**2 + cos(lat0) * cos(lat1) * sin(dlonh)**2
2147 
2148  d = 2.0_rp * asin( min(sqrt(work1),1.0_rp) )
2149 
2150  end function haversine
2151 
2152 
2153 !OCL SERIAL
2154  subroutine spline_coef( &
2155  KA_ref, KS_ref, KE_ref, &
2156  hgt_ref, val_ref, &
2157  spline, &
2158  kmax, &
2159  idx, idx_r, &
2160  U, FDZ )
2161  use scale_const, only: &
2162  undef => const_undef, &
2163  eps => const_eps
2164  use scale_matrix, only: &
2165  matrix_solver_tridiagonal
2166  implicit none
2167  integer, intent(in) :: ka_ref, ks_ref, ke_ref
2168 
2169  real(rp), intent(in) :: hgt_ref(ka_ref)
2170  real(rp), intent(in) :: val_ref(ka_ref)
2171  logical, intent(in) :: spline
2172 
2173  integer, intent(out) :: kmax
2174  integer, intent(out) :: idx (ka_ref)
2175  integer, intent(out) :: idx_r(ka_ref)
2176  real(rp), intent(out) :: u (ka_ref)
2177  real(rp), intent(out) :: fdz(ka_ref)
2178 
2179  real(rp) :: md(ka_ref)
2180  real(rp) :: v (ka_ref)
2181  real(rp) :: dz
2182  integer :: k
2183 
2184  if ( spline ) then
2185 
2186  idx(1) = -999
2187  kmax = 1
2188  do k = ks_ref, ke_ref-1
2189  if ( abs( val_ref(k) - undef ) > eps ) then
2190  idx(1) = k
2191  idx_r(k) = 1
2192  exit
2193  end if
2194  end do
2195  if ( idx(1) == -999 ) return ! UNDEF (use linear interpolation)
2196  fdz(1) = 1e10 ! dummy
2197  do k = idx(1)+1, ke_ref
2198  dz = hgt_ref(k) - hgt_ref(idx(kmax))
2199  if ( abs( val_ref(k) - undef ) > eps .and. dz > eps ) then
2200  do while ( kmax > 1 .and. fdz(kmax) < dz * 0.1_rp )
2201  kmax = kmax - 1 ! marge
2202  end do
2203  kmax = kmax + 1
2204  idx(kmax) = k
2205  if ( idx(kmax-1)+1 <= k-1 ) idx_r(idx(kmax-1)+1:k-1) = kmax-1
2206  idx_r(k) = kmax
2207  fdz(kmax) = hgt_ref(k) - hgt_ref(idx(kmax-1))
2208  end if
2209  end do
2210 
2211  if ( kmax > 3 ) then
2212 
2213  md(2) = 2.0_rp * ( fdz(2) + fdz(3) ) + fdz(2)
2214  do k = 3, kmax-2
2215  md(k) = 2.0_rp * ( fdz(k) + fdz(k+1) )
2216  end do
2217  md(kmax-1) = 2.0_rp * ( fdz(kmax-1) + fdz(kmax) ) + fdz(kmax)
2218 
2219  do k = 2, kmax-1
2220  v(k) = ( val_ref(idx(k+1)) - val_ref(idx(k )) ) / fdz(k+1) &
2221  - ( val_ref(idx(k )) - val_ref(idx(k-1)) ) / fdz(k )
2222  end do
2223 
2224  call matrix_solver_tridiagonal( kmax, 2, kmax-1, &
2225  fdz(2:), md(:), fdz(:), & ! (in)
2226  v(:), & ! (in)
2227  u(:) ) ! (out)
2228 ! U(1) = 0.0_RP
2229 ! U(kmax) = 0.0_RP
2230  u(1) = u(2)
2231  u(kmax) = u(kmax-1)
2232 
2233  else
2234 
2235  idx(kmax) = idx(1) ! force linear interpolateion
2236 
2237  end if
2238 
2239  else
2240 
2241  kmax = 1
2242  idx(1) = -999 ! dummy
2243 
2244  end if
2245 
2246  return
2247  end subroutine spline_coef
2248 
2249 !OCL SERIAL
2250  subroutine spline_exec( &
2251  KA_ref, kmax, KA, KS, KE, &
2252  idx_k, vfact, &
2253  hgt_ref, hgt, &
2254  val_ref, &
2255  idx, idx_r, &
2256  U, FDZ, &
2257  val )
2258  use scale_const, only: &
2259  undef => const_undef, &
2260  eps => const_eps
2261  implicit none
2262  integer, intent(in) :: KA_ref, kmax
2263  integer, intent(in) :: KA, KS, KE
2264  integer, intent(in) :: idx_k(KA,2)
2265  real(RP), intent(in) :: vfact(KA)
2266  real(RP), intent(in) :: hgt_ref(KA_ref)
2267  real(RP), intent(in) :: hgt (KA)
2268  real(RP), intent(in) :: val_ref(KA_ref)
2269  integer, intent(in) :: idx (KA_ref)
2270  integer, intent(in) :: idx_r (KA_ref)
2271  real(RP), intent(in) :: U (KA_ref)
2272  real(RP), intent(in) :: FDZ (KA_ref)
2273 
2274  real(RP), intent(out) :: val(KA)
2275 
2276  real(RP) :: c1, c2, c3, d
2277  integer :: k, kk, kk2
2278 
2279  do k = ks, ke
2280  kk = idx_k(k,1)
2281  if ( kk == -1 ) then
2282  val(k) = undef
2283  else if ( idx_k(k,2) == -1 ) then
2284  val(k) = val_ref(kk)
2285  else if ( kk < idx(1) .or. kk >= idx(kmax) ) then ! linear interpolation
2286  if ( abs( val_ref(kk) - undef ) < eps .or. abs( val_ref(kk+1) - undef ) < eps ) then
2287  val(k) = undef
2288  else
2289  val(k) = val_ref(kk) * vfact(k) + val_ref(kk+1) * ( 1.0_rp - vfact(k) )
2290  end if
2291  else
2292  kk2 = idx_r(kk)
2293  kk = idx(kk2)
2294  c3 = ( u(kk2+1) - u(kk2) ) / fdz(kk2+1)
2295  c2 = 3.0_rp * u(kk2)
2296  c1 = ( val_ref(idx(kk2+1)) - val_ref(kk) ) / fdz(kk2+1) - ( 2.0_rp * u(kk2) + u(kk2+1) ) * fdz(kk2+1)
2297  d = hgt(k) - hgt_ref(kk)
2298 
2299  val(k) = ( ( c3 * d + c2 ) * d + c1 ) * d + val_ref(kk)
2300  end if
2301  end do
2302 
2303  return
2304  end subroutine spline_exec
2305 
2306 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:1232
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:2258
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:342
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:677
scale_precision
module PRECISION
Definition: scale_precision.F90:14
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:33
scale_prof::prof_rapstart
subroutine, public prof_rapstart(rapname_base, level, disable_barrier)
Start raptime.
Definition: scale_prof.F90:159
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:252
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:146
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:1881
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:1134
scale_interp::cross
real(rp) function cross(x0, y0, x1, y1)
cross product
Definition: scale_interp.F90:2115
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:869
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)
Definition: scale_interp.F90:945
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:1322
scale_const::const_radius
real(rp), public const_radius
radius of the planet [m]
Definition: scale_const.F90:44
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:341
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)
Definition: scale_interp.F90:450
scale_matrix
module MATRIX
Definition: scale_matrix.F90:11
scale_prof::prof_rapend
subroutine, public prof_rapend(rapname_base, level, disable_barrier)
Save raptime.
Definition: scale_prof.F90:217
scale_interp::interp_setup
subroutine, public interp_setup(weight_order, search_limit)
Setup.
Definition: scale_interp.F90:88
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:1688
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:1022
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:1978
scale_const::const_d2r
real(rp), public const_d2r
degree to radian
Definition: scale_const.F90:32
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:1925
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:41
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:56
Generated by   doxygen 1.8.17