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  use scale_debug
21  use scale_atmos_grid_cartesc_index
22  use scale_index
23  use scale_tracer
24  !-----------------------------------------------------------------------------
25  implicit none
26  private
27  !-----------------------------------------------------------------------------
28  !
29  !++ Public procedure
30  !
31  public :: interp_setup
32  public :: interp_domain_compatibility
33  public :: interp_factor2d
34  public :: interp_factor3d
35  public :: interp_interp2d
36  public :: interp_interp3d
37 
38  !-----------------------------------------------------------------------------
39  !
40  !++ Public parameters & variables
41  !
42  !-----------------------------------------------------------------------------
43  !
44  !++ Private procedure
45  !
46  private :: interp_search_horiz
47  private :: interp_search_vert
48  private :: interp_insert
49 
50  private :: haversine
51 
52  !-----------------------------------------------------------------------------
53  !
54  !++ Private parameters & variables
55  !
56  real(RP), private, parameter :: large_number = 9.999e+15_rp
57 
58  integer, private :: interp_weight_order = 2
59  real(RP), private :: interp_search_limit
60  real(RP), private :: interp_buffer_size_fact = 2.0_rp
61 
62  !-----------------------------------------------------------------------------
63 contains
64  !-----------------------------------------------------------------------------
66  subroutine interp_setup( &
67  weight_order, &
68  search_limit )
69  use scale_prc, only: &
70  prc_abort
71  implicit none
72 
73  integer, intent(in) :: weight_order
74  real(RP), intent(in), optional :: search_limit
75 
76  namelist /param_interp/ &
77  interp_buffer_size_fact
78 
79  integer :: ierr
80  !---------------------------------------------------------------------------
81 
82  log_newline
83  log_info("INTERP_setup",*) 'Setup'
84 
85  interp_weight_order = weight_order
86 
87  interp_search_limit = large_number
88  if ( present(search_limit) ) then
89  interp_search_limit = search_limit
90  log_info("INTERP_setup",*) 'search limit [m] : ', interp_search_limit
91  endif
92 
93  !--- read namelist
94  rewind(io_fid_conf)
95  read(io_fid_conf,nml=param_interp,iostat=ierr)
96  if( ierr < 0 ) then !--- missing
97  log_info("INTERP_setup",*) 'Not found namelist. Default used.'
98  elseif( ierr > 0 ) then !--- fatal error
99  log_error("INTERP_setup",*) 'Not appropriate names in namelist PARAM_INTERP. Check!'
100  call prc_abort
101  endif
102  log_nml(param_interp)
103 
104  return
105  end subroutine interp_setup
106 
107  !-----------------------------------------------------------------------------
108  subroutine interp_domain_compatibility( &
109  lon_org, &
110  lat_org, &
111  topc_org, &
112  lon_loc, &
113  lat_loc, &
114  topc_loc, &
115  topf_loc, &
116  skip_x, &
117  skip_y, &
118  skip_z )
119  use scale_const, only: &
120  d2r => const_d2r
121  use scale_prc, only: &
122  prc_abort
123  implicit none
124 
125  real(RP), intent(in) :: lon_org (:,:)
126  real(RP), intent(in) :: lat_org (:,:)
127  real(RP), intent(in) :: topc_org(:,:) ! full level
128  real(RP), intent(in) :: lon_loc (:,:)
129  real(RP), intent(in) :: lat_loc (:,:)
130  real(RP), intent(in) :: topc_loc(:,:) ! full level
131  real(RP), intent(in) :: topf_loc(:,:) ! half level (ceil)
132 
133  logical, intent(in), optional :: skip_z
134  logical, intent(in), optional :: skip_x
135  logical, intent(in), optional :: skip_y
136 
137  real(RP) :: max_ref, min_ref
138  real(RP) :: max_loc, min_loc
139 
140  logical :: skip_z_
141  logical :: skip_x_
142  logical :: skip_y_
143 
144  intrinsic size
145  !---------------------------------------------------------------------------
146 
147  skip_z_ = .false.
148  if ( present(skip_z) ) then
149  skip_z_ = skip_z
150  endif
151 
152  skip_x_ = .false.
153  if ( present(skip_x) ) then
154  skip_x_ = skip_x
155  endif
156 
157  skip_y_ = .false.
158  if ( present(skip_y) ) then
159  skip_y_ = skip_y
160  endif
161 
162  if ( .NOT. skip_z_ ) then
163  max_ref = maxval( topc_org(:,:) ) + minval( topf_loc(:,:) - topc_loc(:,:) ) ! not real top boundary (only for check)
164  max_loc = maxval( topc_loc(:,:) )
165 
166  if ( max_ref < max_loc ) then
167  log_error("INTERP_domain_compatibility",*) 'REQUESTED DOMAIN IS TOO MUCH BROAD'
168  log_error_cont(*) '-- VERTICAL direction over the limit'
169  log_error_cont(*) '-- reference max: ', max_ref
170  log_error_cont(*) '-- local max: ', max_loc
171  call prc_abort
172  endif
173  endif
174 
175  if ( .NOT. skip_x_ ) then
176  max_ref = maxval( lon_org(:,:) / d2r )
177  min_ref = minval( lon_org(:,:) / d2r )
178  max_loc = maxval( lon_loc(:,:) / d2r )
179  min_loc = minval( lon_loc(:,:) / d2r )
180 
181  if ( (min_ref+360.0_rp-max_ref) < 360.0_rp / size(lon_org,1) * 2.0_rp ) then
182  ! cyclic OK
183  elseif( max_ref < max_loc .OR. min_ref > min_loc ) then
184  log_error("INTERP_domain_compatibility",*) 'REQUESTED DOMAIN IS TOO MUCH BROAD'
185  log_error_cont(*) '-- LONGITUDINAL direction over the limit'
186  log_error_cont(*) '-- reference max: ', max_ref
187  log_error_cont(*) '-- reference min: ', min_ref
188  log_error_cont(*) '-- local max: ', max_loc
189  log_error_cont(*) '-- local min: ', min_loc
190  call prc_abort
191  endif
192  endif
193 
194  if ( .NOT. skip_y_ ) then
195  max_ref = maxval( lat_org(:,:) / d2r )
196  min_ref = minval( lat_org(:,:) / d2r )
197  max_loc = maxval( lat_loc(:,:) / d2r )
198  min_loc = minval( lat_loc(:,:) / d2r )
199 
200  if ( max_ref < max_loc .OR. min_ref > min_loc ) then
201  log_error("INTERP_domain_compatibility",*) 'REQUESTED DOMAIN IS TOO MUCH BROAD'
202  log_error_cont(*) '-- LATITUDINAL direction over the limit'
203  log_error_cont(*) '-- reference max: ', max_ref
204  log_error_cont(*) '-- reference min: ', min_ref
205  log_error_cont(*) '-- local max: ', max_loc
206  log_error_cont(*) '-- local min: ', min_loc
207  call prc_abort
208  endif
209  endif
210 
211  return
212  end subroutine interp_domain_compatibility
213 
214  !-----------------------------------------------------------------------------
215  ! make interpolation factor using Lat-Lon
216  subroutine interp_factor2d( &
217  npoints, &
218  IA_ref, JA_ref, &
219  lon_ref,lat_ref, &
220  IA, JA, &
221  lon, lat, &
222  idx_i, idx_j, hfact, &
223  search_limit, &
224  latlon_structure, &
225  weight_order )
226  use scale_const, only: &
227  undef => const_undef
228  use scale_prc, only: &
229  prc_abort
230  implicit none
231 
232  integer, intent(in) :: npoints ! number of interpolation point for horizontal
233  integer, intent(in) :: IA_ref ! number of x-direction (reference)
234  integer, intent(in) :: JA_ref ! number of y-direction (reference)
235  real(RP), intent(in) :: lon_ref(ia_ref,ja_ref) ! longitude [rad] (reference)
236  real(RP), intent(in) :: lat_ref(ia_ref,ja_ref) ! latitude [rad] (reference)
237  integer, intent(in) :: IA ! number of x-direction (target)
238  integer, intent(in) :: JA ! number of y-direction (target)
239  real(RP), intent(in) :: lon (ia,ja) ! longitude [rad] (target)
240  real(RP), intent(in) :: lat (ia,ja) ! latitude [rad] (target)
241  integer, intent(out) :: idx_i(ia,ja,npoints) ! i-index in reference (target)
242  integer, intent(out) :: idx_j(ia,ja,npoints) ! j-index in reference (target)
243  real(RP), intent(out) :: hfact(ia,ja,npoints) ! horizontal interp factor (target)
244 
245  real(RP), intent(in), optional :: search_limit
246  logical, intent(in), optional :: latlon_structure
247  integer, intent(in), optional :: weight_order
248 
249  logical :: ll_struct_
250 
251  real(RP) :: lon_min, lon_max
252  real(RP) :: lat_min, lat_max
253  real(RP) :: dlon, dlat
254 
255  ! for structure grid
256  integer :: is, ie, js, je
257  integer :: psizex, psizey
258  real(RP) :: lon1d(ia_ref), lat1d(ja_ref)
259  real(RP) :: lon0, lat0
260  integer, allocatable :: i0(:), i1(:), j0(:), j1(:)
261 
262  ! for unstructure grid
263  integer :: nsize, psize, nidx_max
264  integer, allocatable :: idx_blk(:,:,:), nidx(:,:)
265  integer :: idx_ref(npoints)
266 
267 
268  integer :: i, j, ii, jj, n
269  !---------------------------------------------------------------------------
270 
271  call prof_rapstart('INTERP_fact',3)
272 
273  if ( present(latlon_structure) ) then
274  ll_struct_ = latlon_structure
275  else
276  ll_struct_ = .false.
277  end if
278 
279  hfact(:,:,:) = 0.0_rp
280 
281 
282  if ( ll_struct_ ) then
283 
284  lon1d(:) = undef
285  lon_min = undef
286  do i = 1, ia_ref
287  do j = 1, ja_ref
288  if ( lon_ref(i,j) .ne. undef ) then
289  lon1d(i) = lon_ref(i,j)
290  exit
291  end if
292  end do
293  if ( lon_min == undef .and. lon1d(i) .ne. undef ) then
294  lon_min = lon1d(i)
295  is = i
296  endif
297  if ( lon_min .ne. undef ) then
298  if ( lon1d(i) .ne. undef ) then
299  lon_max = lon1d(i)
300  ie = i
301  end if
302  end if
303  end do
304 
305  lat1d(:) = undef
306  lat_min = undef
307  do j = 1, ja_ref
308  do i = 1, ia_ref
309  if ( lat_ref(i,j) .ne. undef ) then
310  lat1d(j) = lat_ref(i,j)
311  exit
312  end if
313  end do
314  if ( lat_min == undef .and. lat1d(j) .ne. undef ) then
315  lat_min = lat1d(j)
316  js = j
317  endif
318  if ( lat_min .ne. undef ) then
319  if ( lat1d(j) .ne. undef ) then
320  lat_max = lat1d(j)
321  je = j
322  end if
323  end if
324  end do
325 
326  ! fill undef
327  dlon = ( lon_max - lon_min ) / ( ie - is )
328  do i = is, ie
329  if ( lon1d(i) == undef ) lon1d(i) = lon_min + dlon * ( i - is )
330  end do
331  dlat = ( lat_max - lat_min ) / ( je - js )
332  do j = js, je
333  if ( lat1d(j) == undef ) lat1d(j) = lat_min + dlat * ( j - js )
334  end do
335 
336 
337 
338  if ( ie-is > 10 ) then
339  psizex = int( 2.0_rp*sqrt(real(ie-is+1,rp)) )
340  else
341  psizex = 1
342  end if
343  if ( je-js > 10 ) then
344  psizey = int( 2.0_rp*sqrt(real(je-js+1,rp)) )
345  else
346  psizey = 1
347  end if
348 
349  allocate( i0(psizex), i1(psizex) )
350  allocate( j0(psizey), j1(psizey) )
351 
352  dlon = ( lon_max - lon_min ) / psizex
353  dlat = ( lat_max - lat_min ) / psizey
354 
355  do ii = 1, psizex
356  lon0 = lon_min + dlon * (ii-1)
357  do i = is, ie
358  if ( lon1d(i) >= lon0 ) then
359  i0(ii) = i
360  exit
361  end if
362  end do
363  end do
364  do ii = 1, psizex-1
365  i1(ii) = i0(ii+1) - 1
366  end do
367  i1(psizex) = ie
368 
369  do jj = 1, psizey
370  lat0 = lat_min + dlat * (jj-1)
371  do j = js, je
372  if ( lat1d(j) >= lat0 ) then
373  j0(jj) = j
374  exit
375  end if
376  end do
377  end do
378  do jj = 1, psizey-1
379  j1(jj) = j0(jj+1) - 1
380  end do
381  j1(psizey) = je
382 
383  !$omp parallel do OMP_SCHEDULE_ collapse(2)
384  do j = 1, ja
385  do i = 1, ia
386  ! main search
387  call interp_search_horiz_struct( npoints, & ! [IN]
388  psizex, psizey, & ! [IN]
389  ia_ref, ja_ref, & ! [IN]
390  lon1d(:), lat1d(:), & ! [IN]
391  lon_min, lat_min, & ! [IN]
392  dlon, dlat, & ! [IN]
393  i0(:), i1(:), j0(:), j1(:), & ! [IN]
394  lon(i,j), lat(i,j), & ! [IN]
395  idx_i(i,j,:), idx_j(i,j,:), & ! [OUT]
396  hfact(i,j,:), & ! [OUT]
397  search_limit = search_limit, & ! [IN]
398  weight_order = weight_order ) ! [IN]
399  enddo
400  enddo
401 
402  deallocate( i0, i1, j0, j1 )
403 
404  else
405 
406  nsize = ia_ref * ja_ref
407  if ( nsize > 100 ) then
408  psize = int( sqrt(2.0_rp*sqrt(real(nsize,rp))) )
409  nidx_max = nsize / psize * interp_buffer_size_fact
410  else
411  psize = 1
412  nidx_max = nsize
413  end if
414 
415  allocate(idx_blk(nidx_max,psize,psize))
416  allocate(nidx( psize,psize))
417 
418  call interp_div_block(nsize, psize, nidx_max, & ! [IN]
419  lon_ref(:,:), lat_ref(:,:), & ! [IN]
420  idx_blk(:,:,:), nidx(:,:), & ! [OUT]
421  lon_min, lon_max, & ! [OUT]
422  lat_min, lat_max, & ! [OUT]
423  dlon, dlat ) ! [OUT]
424 
425  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
426  !$omp private(idx_ref)
427  do j = 1, ja
428  do i = 1, ia
429  ! main search
430  call interp_search_horiz( npoints, & ! [IN]
431  nsize, & ! [IN]
432  lon_ref(:,:), lat_ref(:,:), & ! [IN]
433  lon_min, lon_max, & ! [IN]
434  lat_min, lat_max, & ! [IN]
435  psize, nidx_max, & ! [IN]
436  dlon, dlat, & ! [IN]
437  idx_blk(:,:,:), nidx(:,:), & ! [IN]
438  lon(i,j), lat(i,j), & ! [IN]
439  idx_ref(:), & ! [OUT]
440  hfact(i,j,:), & ! [OUT]
441  search_limit = search_limit, & ! [IN]
442  weight_order = weight_order ) ! [IN]
443  do n = 1, npoints
444  idx_i(i,j,n) = mod(idx_ref(n) - 1, ia_ref) + 1
445  idx_j(i,j,n) = ( idx_ref(n) - 1 ) / ia_ref + 1
446  end do
447  enddo
448  enddo
449 
450  deallocate(idx_blk, nidx)
451 
452  end if
453 
454  call prof_rapend ('INTERP_fact',3)
455 
456  return
457  end subroutine interp_factor2d
458 
459  !-----------------------------------------------------------------------------
460  ! make interpolation factor using Lat-Lon and Z-Height information
461  subroutine interp_factor3d( &
462  npoints, &
463  KA_ref, &
464  KS_ref, &
465  KE_ref, &
466  IA_ref, &
467  JA_ref, &
468  lon_ref, &
469  lat_ref, &
470  hgt_ref, &
471  KA, &
472  KS, &
473  KE, &
474  IA, &
475  JA, &
476  lon, &
477  lat, &
478  hgt, &
479  idx_i, &
480  idx_j, &
481  hfact, &
482  idx_k, &
483  vfact )
484  use scale_prc, only: &
485  prc_abort
486  implicit none
487 
488  integer, intent(in) :: npoints ! number of interpolation point for horizontal
489  integer, intent(in) :: KA_ref, KS_ref, KE_ref ! number of z-direction (reference)
490  integer, intent(in) :: IA_ref ! number of x-direction (reference)
491  integer, intent(in) :: JA_ref ! number of y-direction (reference)
492  real(RP), intent(in) :: lon_ref(ia_ref,ja_ref) ! longitude [rad] (reference)
493  real(RP), intent(in) :: lat_ref(ia_ref,ja_ref) ! latitude [rad] (reference)
494  real(RP), intent(in) :: hgt_ref(ka_ref,ia_ref,ja_ref) ! height [m] (reference)
495  integer, intent(in) :: KA, KS, KE ! number of z-direction (target)
496  integer, intent(in) :: IA ! number of x-direction (target)
497  integer, intent(in) :: JA ! number of y-direction (target)
498  real(RP), intent(in) :: lon (ia,ja) ! longitude [rad] (target)
499  real(RP), intent(in) :: lat (ia,ja) ! latitude [rad] (target)
500  real(RP), intent(in) :: hgt (ka,ia,ja) ! longitude [m] (target)
501  integer, intent(out) :: idx_i(ia,ja,npoints) ! i-index in reference (target)
502  integer, intent(out) :: idx_j(ia,ja,npoints) ! j-index in reference (target)
503  real(RP), intent(out) :: hfact(ia,ja,npoints) ! horizontal interp factor (target)
504  integer, intent(out) :: idx_k(ka,2,ia,ja,npoints) ! i-index in reference (target)
505  real(RP), intent(out) :: vfact(ka,2,ia,ja,npoints) ! horizontal interp factor (target)
506 
507  integer :: nsize, psize, nidx_max
508  integer, allocatable :: idx_blk(:,:,:), nidx(:,:)
509  real(RP) :: lon_min, lon_max
510  real(RP) :: lat_min, lat_max
511  real(RP) :: dlon, dlat
512  integer :: idx_ref(npoints)
513 
514  integer :: i, j, ii, jj, n
515  !---------------------------------------------------------------------------
516 
517  call prof_rapstart('INTERP_fact',3)
518 
519  nsize = ia_ref * ja_ref
520  if ( nsize > 100 ) then
521  psize = int( sqrt(2.0_rp*sqrt(real(nsize,rp))) )
522  nidx_max = nsize / psize * interp_buffer_size_fact
523  else
524  psize = 1
525  nidx_max = nsize
526  end if
527 
528  allocate(idx_blk(nidx_max,psize,psize))
529  allocate(nidx( psize,psize))
530 
531  call interp_div_block(nsize, psize, nidx_max, & ! [IN]
532  lon_ref(:,:), lat_ref(:,:), & ! [IN]
533  idx_blk(:,:,:), nidx(:,:), & ! [OUT]
534  lon_min, lon_max, & ! [OUT]
535  lat_min, lat_max, & ! [OUT]
536  dlon, dlat ) ! [OUT]
537 
538  hfact(:,:,:) = 0.0_rp
539  vfact(:,:,:,:,:) = 0.0_rp
540 
541  !$omp parallel do OMP_SCHEDULE_ collapse(2) &
542  !$omp private(ii,jj,idx_ref)
543  do j = 1, ja
544  do i = 1, ia
545 
546  ! main search
547  call interp_search_horiz( npoints, & ! [IN]
548  nsize, & ! [IN]
549  lon_ref(:,:), lat_ref(:,:), & ! [IN]
550  lon_min, lon_max, & ! [IN]
551  lat_min, lat_max, & ! [IN]
552  psize, nidx_max, & ! [IN]
553  dlon, dlat, & ! [IN]
554  idx_blk(:,:,:), nidx(:,:), & ! [IN]
555  lon(i,j), lat(i,j), & ! [IN]
556  idx_ref(:), hfact(i,j,:) ) ! [OUT]
557 
558  do n = 1, npoints
559  ii = mod(idx_ref(n) - 1, ia_ref) + 1
560  jj = ( idx_ref(n) - 1 ) / ia_ref + 1
561  idx_i(i,j,n) = ii
562  idx_j(i,j,n) = jj
563 
564  call interp_search_vert( ka_ref, ks_ref, ke_ref, & ! [IN]
565  ka, ks, ke, & ! [IN]
566  hgt_ref(:,ii,jj), & ! [IN]
567  hgt(:,i,j), & ! [IN]
568  idx_k(:,:,i,j,n), & ! [OUT]
569  vfact(:,:,i,j,n) ) ! [OUT]
570  enddo
571  enddo
572  enddo
573 
574  deallocate(idx_blk, nidx)
575 
576  call prof_rapend ('INTERP_fact',3)
577 
578  return
579  end subroutine interp_factor3d
580 
581  !-----------------------------------------------------------------------------
582  ! interpolation using one-points for 2D data (nearest-neighbor)
583  subroutine interp_interp2d( &
584  npoints, &
585  IA_ref, &
586  JA_ref, &
587  IA, &
588  JA, &
589  idx_i, &
590  idx_j, &
591  hfact, &
592  val_ref, &
593  val )
594  use scale_const, only: &
595  undef => const_undef, &
596  eps => const_eps
597  implicit none
598 
599  integer, intent(in) :: npoints ! number of interpolation point for horizontal
600  integer, intent(in) :: IA_ref ! number of x-direction (reference)
601  integer, intent(in) :: JA_ref ! number of y-direction (reference)
602  integer, intent(in) :: IA ! number of x-direction (target)
603  integer, intent(in) :: JA ! number of y-direction (target)
604  integer, intent(in) :: idx_i (ia,ja,npoints) ! i-index in reference (target)
605  integer, intent(in) :: idx_j (ia,ja,npoints) ! j-index in reference (target)
606  real(RP), intent(in) :: hfact (ia,ja,npoints) ! horizontal interp factor (target)
607  real(RP), intent(in) :: val_ref(ia_ref,ja_ref) ! value (reference)
608  real(RP), intent(out) :: val (ia,ja) ! value (target)
609 
610  integer :: i, j, n
611  !---------------------------------------------------------------------------
612 
613  call prof_rapstart('INTERP_interp',3)
614 
615  !$omp parallel do OMP_SCHEDULE_ collapse(2)
616 !OCL PREFETCH
617  do j = 1, ja
618  do i = 1, ia
619  if ( hfact(i,j,1) < eps .or. val_ref(idx_i(i,j,1),idx_j(i,j,1)) == undef ) then
620  val(i,j) = undef
621  else
622  val(i,j) = hfact(i,j,1) * val_ref(idx_i(i,j,1),idx_j(i,j,1))
623  do n = 2, npoints
624  if ( val_ref(idx_i(i,j,n),idx_j(i,j,n)) == undef ) then
625  val(i,j) = val(i,j) / sum( hfact(i,j,1:n-1) )
626  exit
627  end if
628  val(i,j) = val(i,j) &
629  + hfact(i,j,n) * val_ref(idx_i(i,j,n),idx_j(i,j,n))
630  end do
631  end if
632  enddo
633  enddo
634 
635  call prof_rapend ('INTERP_interp',3)
636 
637  return
638  end subroutine interp_interp2d
639 
640  !-----------------------------------------------------------------------------
641  ! interpolation using one-points for 3D data (nearest-neighbor)
642  subroutine interp_interp3d( &
643  npoints, &
644  KA_ref, &
645  IA_ref, &
646  JA_ref, &
647  KA, KS, KE, &
648  IA, &
649  JA, &
650  idx_i, &
651  idx_j, &
652  hfact, &
653  idx_k, &
654  vfact, &
655  val_ref, &
656  val, &
657  logwgt )
658  implicit none
659 
660  integer, intent(in) :: npoints ! number of interpolation point for horizontal
661  integer, intent(in) :: KA_ref ! number of x-direction (reference)
662  integer, intent(in) :: IA_ref ! number of x-direction (reference)
663  integer, intent(in) :: JA_ref ! number of y-direction (reference)
664  integer, intent(in) :: KA, KS, KE ! number of z-direction (target)
665  integer, intent(in) :: IA ! number of x-direction (target)
666  integer, intent(in) :: JA ! number of y-direction (target)
667  integer, intent(in) :: idx_i (ia,ja,npoints) ! i-index in reference (target)
668  integer, intent(in) :: idx_j (ia,ja,npoints) ! j-index in reference (target)
669  real(RP), intent(in) :: hfact (ia,ja,npoints) ! horizontal interp factor (target)
670  integer, intent(in) :: idx_k (ka,2,ia,ja,npoints) ! k-index in reference (target)
671  real(RP), intent(in) :: vfact (ka,2,ia,ja,npoints) ! vertical interp factor (target)
672 
673  real(RP), intent(in), target :: val_ref(ka_ref,ia_ref,ja_ref) ! value (reference)
674 
675  real(RP), intent(out) :: val (ka,ia,ja) ! value (target)
676 
677  logical, intent(in), optional:: logwgt ! use logarithmic weighted interpolation?
678 
679  logical :: logwgt_
680 
681  real(RP), pointer :: work(:,:,:)
682 
683  integer :: k, i, j, n
684  !---------------------------------------------------------------------------
685 
686  call prof_rapstart('INTERP_interp',3)
687 
688  logwgt_ = .false.
689  if ( present(logwgt) ) then
690  logwgt_ = logwgt
691  endif
692 
693  if ( logwgt_ ) then
694  allocate( work(ka_ref,ia_ref,ja_ref) )
695  !$omp parallel do OMP_SCHEDULE_ collapse(2)
696  do j = 1, ja_ref
697  do i = 1, ia_ref
698  do k = 1, ka_ref
699  work(k,i,j) = log( val_ref(k,i,j) )
700  enddo
701  enddo
702  enddo
703  else
704  work => val_ref
705  endif
706 
707  !$omp parallel do OMP_SCHEDULE_ collapse(2)
708  do j = 1, ja
709  do i = 1, ia
710  do k = ks, ke
711  val(k,i,j) = hfact(i,j,1) * vfact(k,1,i,j,1) * work(idx_k(k,1,i,j,1),idx_i(i,j,1),idx_j(i,j,1)) &
712  + hfact(i,j,1) * vfact(k,2,i,j,1) * work(idx_k(k,2,i,j,1),idx_i(i,j,1),idx_j(i,j,1))
713  enddo
714  enddo
715  enddo
716 
717 !OCL SERIAL
718  do n = 2, npoints
719  !$omp parallel do OMP_SCHEDULE_ collapse(2)
720 !OCL PREFETCH
721  do j = 1, ja
722  do i = 1, ia
723  do k = ks, ke
724  val(k,i,j) = val(k,i,j) &
725  + hfact(i,j,n) * vfact(k,1,i,j,n) * work(idx_k(k,1,i,j,n),idx_i(i,j,n),idx_j(i,j,n)) &
726  + hfact(i,j,n) * vfact(k,2,i,j,n) * work(idx_k(k,2,i,j,n),idx_i(i,j,n),idx_j(i,j,n))
727  enddo
728  enddo
729  enddo
730  enddo
731 
732  if ( logwgt_ ) then
733  deallocate( work )
734  !$omp parallel do OMP_SCHEDULE_ collapse(2)
735  do j = 1, ja
736  do i = 1, ia
737  do k = ks, ke
738  val(k,i,j) = exp( val(k,i,j) )
739  end do
740  end do
741  end do
742  endif
743 
744  call prof_rapend ('INTERP_interp',3)
745 
746  return
747  end subroutine interp_interp3d
748 
749  !-----------------------------------------------------------------------------
750  ! horizontal search of interpolation points
751 !OCL SERIAL
752  subroutine interp_search_horiz( &
753  npoints, &
754  nsize, &
755  lon_ref, lat_ref, &
756  lon_min, lon_max, &
757  lat_min, lat_max, &
758  psize, nidx_max, &
759  dlon, dlat, &
760  idx_blk, nidx, &
761  lon, lat, &
762  idx_ref, &
763  hfact, &
764  search_limit, &
765  weight_order )
766  use scale_const, only: &
767  eps => const_eps, &
768  radius => const_radius
769  implicit none
770 
771  integer, intent(in) :: npoints ! number of interpolation point for horizontal
772  integer, intent(in) :: nsize ! number of grids (reference)
773  real(RP), intent(in) :: lon_ref(nsize) ! longitude [rad] (reference)
774  real(RP), intent(in) :: lat_ref(nsize) ! latitude [rad] (reference)
775  real(RP), intent(in) :: lon_min ! minmum longitude [rad] (reference)
776  real(RP), intent(in) :: lon_max ! maximum longitude [rad] (reference)
777  real(RP), intent(in) :: lat_min ! minmum latitude [rad] (reference)
778  real(RP), intent(in) :: lat_max ! maximum latitude [rad] (reference)
779  integer, intent(in) :: psize ! number of blocks for each dimension
780  integer, intent(in) :: nidx_max ! maximum number of index in the block
781  real(RP), intent(in) :: dlon ! block longitude difference
782  real(RP), intent(in) :: dlat ! block latitude difference
783  integer, intent(in) :: idx_blk(nidx_max,psize,psize) ! index of the reference in the block
784  integer, intent(in) :: nidx (psize,psize) ! number of indexes in the block
785  real(RP), intent(in) :: lon ! longitude [rad] (target)
786  real(RP), intent(in) :: lat ! latitude [rad] (target)
787  integer, intent(out) :: idx_ref(npoints) ! index in reference (target)
788  real(RP), intent(out) :: hfact(npoints) ! horizontal interp factor (target)
789 
790  real(RP), intent(in), optional :: search_limit
791  integer, intent(in), optional :: weight_order
792 
793  real(RP) :: drad(npoints)
794  real(RP) :: sum
795  real(RP) :: search_limit_
796  integer :: weight_order_
797 
798  real(RP) :: lon0, lon1, lat0, lat1
799  real(RP) :: dlon_sl, dlat_sl
800  integer :: i, n
801  integer :: ii, jj, ii0, jj0
802  !---------------------------------------------------------------------------
803 
804  if ( present(search_limit) ) then
805  search_limit_ = search_limit
806  else
807  search_limit_ = interp_search_limit
808  end if
809  search_limit_ = search_limit_ / radius ! m to radian
810 
811  if ( present(weight_order) ) then
812  weight_order_ = weight_order
813  else
814  weight_order_ = interp_weight_order
815  end if
816 
817 
818 
819  drad(:) = large_number
820  idx_ref(:) = -1
821 
822  ! find k-nearest points in the nearest block
823  ii0 = max( min( int(( lon - lon_min ) / dlon) + 1, psize ), 1 )
824  jj0 = max( min( int(( lat - lat_min ) / dlat) + 1, psize ), 1 )
825  do i = 1, nidx(ii0,jj0)
826  n = idx_blk(i,ii0,jj0)
827  call interp_insert( npoints, &
828  lon, lat, &
829  lon_ref(n), lat_ref(n), & ! [IN]
830  n, & ! [IN]
831  drad(:), idx_ref(:) ) ! [INOUT]
832  end do
833 
834  if ( abs(drad(1)) < eps ) then
835  hfact(:) = 0.0_rp
836  hfact(1) = 1.0_rp
837 
838  return
839  else if ( drad(1) > search_limit_ ) then
840  hfact(:) = 0.0_rp
841  idx_ref(:) = 1 ! dummy
842 
843  return
844  end if
845 
846  dlon_sl = max(dlon * 0.5_rp, drad(npoints))
847  dlat_sl = max(dlat * 0.5_rp, drad(npoints))
848  do jj = 1, psize
849  lat0 = lat_min + dlat * (jj-1)
850  lat1 = lat_min + dlat * jj
851  if ( lat < lat0 - dlat_sl &
852  .or. lat >= lat1 + dlat_sl ) cycle
853  do ii = 1, psize
854  if ( ii==ii0 .and. jj==jj0 ) cycle
855  lon0 = lon_min + dlon * (ii-1)
856  lon1 = lon_min + dlon * ii
857  if ( lon < lon0 - dlon_sl &
858  .or. lon >= lon1 + dlon_sl ) cycle
859  do i = 1, nidx(ii,jj)
860  n = idx_blk(i,ii,jj)
861  call interp_insert( npoints, &
862  lon, lat, &
863  lon_ref(n), lat_ref(n), & ! [IN]
864  n, & ! [IN]
865  drad(:), idx_ref(:) ) ! [INOUT]
866  end do
867  dlon_sl = max(dlon * 0.5_rp, drad(npoints))
868  dlat_sl = max(dlat * 0.5_rp, drad(npoints))
869  end do
870  end do
871 
872  ! factor = 1 / dradian
873  if ( weight_order_ < 0 ) then
874  do n = 1, npoints
875  hfact(n) = 1.0_rp / drad(n)**(-1.0_rp/weight_order_)
876  enddo
877  else
878  do n = 1, npoints
879  hfact(n) = 1.0_rp / drad(n)**weight_order_
880  enddo
881  end if
882 
883  ! ignore far point
884  do n = 1, npoints
885  if ( drad(n) >= search_limit_ ) then
886  hfact(n) = 0.0_rp
887  endif
888  enddo
889 
890  ! normalize factor
891  sum = 0.0_rp
892  do n = 1, npoints
893  sum = sum + hfact(n)
894  enddo
895 
896  if ( sum > 0.0_rp ) then
897  do n = 1, npoints
898  hfact(n) = hfact(n) / sum
899  enddo
900  endif
901 
902  return
903  end subroutine interp_search_horiz
904 
905  !-----------------------------------------------------------------------------
906  ! horizontal search of interpolation points for structure grid
907 !OCL SERIAL
908  subroutine interp_search_horiz_struct( &
909  npoints, &
910  psizex, psizey, &
911  IA_ref, JA_ref, &
912  lon_ref, lat_ref, &
913  lon_min, lat_min, &
914  dlon, dlat, &
915  i0, i1, j0, j1, &
916  lon, lat, &
917  idx_i, idx_j, &
918  hfact, &
919  search_limit, &
920  weight_order )
921  use scale_const, only: &
922  eps => const_eps, &
923  radius => const_radius
924  implicit none
925 
926  integer, intent(in) :: npoints ! number of interpolation point for horizontal
927  integer, intent(in) :: psizex ! number of block in x-direction (reference)
928  integer, intent(in) :: psizey ! number of block in y-direction (reference)
929  integer, intent(in) :: IA_ref ! number of grids in x-direction (reference)
930  integer, intent(in) :: JA_ref ! number of grids in y-direction (reference)
931  real(RP), intent(in) :: lon_ref(ia_ref) ! longitude [rad] (reference)
932  real(RP), intent(in) :: lat_ref(ja_ref) ! latitude [rad] (reference)
933  real(RP), intent(in) :: lon_min ! minimum longitude
934  real(RP), intent(in) :: lat_min ! minimum latitude
935  real(RP), intent(in) :: dlon
936  real(RP), intent(in) :: dlat
937  integer, intent(in) :: i0(psizex) ! start index in the block
938  integer, intent(in) :: i1(psizex) ! end index in the block
939  integer, intent(in) :: j0(psizey) ! start index in the block
940  integer, intent(in) :: j1(psizey) ! start index in the block
941  real(RP), intent(in) :: lon ! longitude [rad] (target)
942  real(RP), intent(in) :: lat ! latitude [rad] (target)
943  integer, intent(out) :: idx_i(npoints) ! x-index in reference (target)
944  integer, intent(out) :: idx_j(npoints) ! y-index in reference (target)
945  real(RP), intent(out) :: hfact(npoints) ! horizontal interp factor (target)
946 
947  real(RP), intent(in), optional :: search_limit
948  integer, intent(in), optional :: weight_order
949 
950  real(RP) :: drad(npoints)
951  real(RP) :: sum
952  real(RP) :: search_limit_
953  integer :: weight_order_
954 
955  real(RP) :: lon0, lon1, lat0, lat1
956  real(RP) :: dlon_sl, dlat_sl
957 
958  integer :: i, j, n
959  integer :: ii, jj, ii0, jj0
960  !---------------------------------------------------------------------------
961 
962  if ( present(search_limit) ) then
963  search_limit_ = search_limit
964  else
965  search_limit_ = interp_search_limit
966  end if
967  search_limit_ = search_limit_ / radius ! m to radian
968 
969  if ( present(weight_order) ) then
970  weight_order_ = weight_order
971  else
972  weight_order_ = interp_weight_order
973  end if
974 
975 
976  drad(:) = large_number
977  idx_i(:) = 1
978  idx_j(:) = 1
979 
980  ! find k-nearest points in the nearest block
981  ii0 = max( min( int(( lon - lon_min ) / dlon) + 1, psizex ), 1)
982  jj0 = max( min( int(( lat - lat_min ) / dlat) + 1, psizey ), 1)
983  do j = j0(jj0), j1(jj0)
984  do i = i0(ii0), i1(ii0)
985  call interp_insert_2d( npoints, &
986  lon, lat, &
987  lon_ref(i), lat_ref(j), & ! [IN]
988  i, j, & ! [IN]
989  drad(:), idx_i(:), idx_j(:) ) ! [INOUT]
990  end do
991  end do
992 
993  if ( abs(drad(1)) < eps ) then
994  hfact(:) = 0.0_rp
995  hfact(1) = 1.0_rp
996 
997  return
998  else if ( drad(1) > search_limit_ ) then
999  hfact(:) = 0.0_rp
1000  idx_i(:) = 1 ! dummy
1001  idx_j(:) = 1 ! dummy
1002 
1003  return
1004  end if
1005 
1006  dlon_sl = max(dlon * 0.5_rp, drad(npoints))
1007  dlat_sl = max(dlat * 0.5_rp, drad(npoints))
1008  do jj = 1, psizey
1009  lat0 = lat_min + dlat * (jj-1)
1010  lat1 = lat_min + dlat * jj
1011  if ( lat < lat0 - dlat_sl &
1012  .or. lat >= lat1 + dlat_sl ) cycle
1013  do ii = 1, psizex
1014  if ( ii==ii0 .and. jj==jj0 ) cycle
1015  lon0 = lon_min + dlon * (ii-1)
1016  lon1 = lon_min + dlon * ii
1017  if ( lon < lon0 - dlon_sl &
1018  .or. lon >= lon1 + dlon_sl ) cycle
1019  do j = j0(jj0), j1(jj0)
1020  do i = i0(ii0), i1(ii0)
1021  call interp_insert_2d( npoints, &
1022  lon, lat, &
1023  lon_ref(i), lat_ref(j), & ! [IN]
1024  i, j, & ! [IN]
1025  drad(:), idx_i(:), idx_j(:) ) ! [INOUT]
1026  end do
1027  end do
1028  dlon_sl = max(dlon * 0.5_rp, drad(npoints))
1029  dlat_sl = max(dlat * 0.5_rp, drad(npoints))
1030  end do
1031  end do
1032 
1033 
1034  ! factor = 1 / drad
1035  if ( weight_order_ < 0 ) then
1036  do n = 1, npoints
1037  hfact(n) = 1.0_rp / drad(n)**(-1.0_rp/weight_order_)
1038  enddo
1039  else
1040  do n = 1, npoints
1041  hfact(n) = 1.0_rp / drad(n)**weight_order_
1042  enddo
1043  end if
1044 
1045  ! ignore far point
1046  do n = 1, npoints
1047  if ( drad(n) >= search_limit_ ) then
1048  hfact(n) = 0.0_rp
1049  endif
1050  enddo
1051 
1052  ! normalize factor
1053  sum = 0.0_rp
1054  do n = 1, npoints
1055  sum = sum + hfact(n)
1056  enddo
1057 
1058  if ( sum > 0.0_rp ) then
1059  do n = 1, npoints
1060  hfact(n) = hfact(n) / sum
1061  enddo
1062  endif
1063 
1064  return
1065  end subroutine interp_search_horiz_struct
1066 
1067  !-----------------------------------------------------------------------------
1068  ! vertical search of interpolation points for two-points
1069 !OCL SERIAL
1070  subroutine interp_search_vert( &
1071  KA_ref, KS_ref, KE_ref, &
1072  KA, KS, KE, &
1073  hgt_ref, &
1074  hgt, &
1075  idx_k, &
1076  vfact )
1077  use scale_prc, only: &
1078  prc_abort
1079  implicit none
1080 
1081  integer, intent(in) :: KA_ref, KS_ref, KE_ref ! number of z-direction (reference)
1082  integer, intent(in) :: KA, KS, KE ! number of z-direction (target)
1083  real(RP), intent(in) :: hgt_ref(ka_ref) ! height [m] (reference)
1084  real(RP), intent(in) :: hgt (ka) ! height [m] (target)
1085  integer, intent(inout) :: idx_k (ka,2) ! k-index in reference (target)
1086  real(RP), intent(inout) :: vfact (ka,2) ! horizontal interp factor (target)
1087 
1088  real(RP) :: weight
1089  integer :: k, kk
1090  !---------------------------------------------------------------------------
1091 
1092  do k = ks, ke
1093  idx_k(k,1) = -1
1094  idx_k(k,2) = -1
1095 
1096  if ( hgt(k) < hgt_ref(ks_ref) ) then
1097  idx_k(k,1) = ks_ref
1098  idx_k(k,2) = ks_ref
1099  vfact(k,1) = 1.0_rp
1100  vfact(k,2) = 0.0_rp
1101  elseif( hgt(k) >= hgt_ref(ke_ref) ) then
1102  idx_k(k,1) = ke_ref
1103  idx_k(k,2) = ke_ref
1104  vfact(k,1) = 1.0_rp
1105  vfact(k,2) = 0.0_rp
1106  else
1107  do kk = ks_ref, ke_ref-1
1108  if ( hgt(k) >= hgt_ref(kk ) &
1109  .AND. hgt(k) < hgt_ref(kk+1) ) then
1110 
1111  weight = ( hgt(k) - hgt_ref(kk) ) &
1112  / ( hgt_ref(kk+1) - hgt_ref(kk) )
1113 
1114  idx_k(k,1) = kk
1115  idx_k(k,2) = kk + 1
1116  vfact(k,1) = 1.0_rp - weight
1117  vfact(k,2) = weight
1118  exit
1119  endif
1120  enddo
1121  endif
1122 
1123  if ( idx_k(k,1) < 0 ) then
1124  log_error("INTERP_search_vert",*) 'data for interpolation was not found.'
1125  log_error_cont(*) 'k=', k, ', hgt(k)=', hgt(k), ', hgt_ref(:)=', hgt_ref(:)
1126  call prc_abort
1127  endif
1128 
1129  enddo ! k-loop
1130 
1131  return
1132  end subroutine interp_search_vert
1133 
1134  ! private
1135 
1136 !OCL SERIAL
1137  subroutine interp_insert( npoints, &
1138  lon, lat, &
1139  lon_ref, lat_ref, &
1140  i, &
1141  drad, idx_i )
1142  use scale_const, only: &
1143  undef => const_undef
1144  use scale_sort, only: &
1145  sort_exec
1146 
1147  integer, intent(in) :: npoints
1148  real(RP), intent(in) :: lon, lat
1149  real(RP), intent(in) :: lon_ref, lat_ref
1150  integer, intent(in) :: i
1151 
1152  real(RP), intent(inout) :: drad(npoints)
1153  integer, intent(inout) :: idx_i(npoints)
1154 
1155  real(RP) :: dradian
1156 
1157  if ( lon_ref == undef ) return
1158 
1159  dradian = haversine( lon, lat, lon_ref, lat_ref )
1160 
1161  if ( dradian <= drad(npoints) ) then
1162  ! replace last(=longest) value
1163  drad(npoints) = dradian
1164  idx_i(npoints) = i
1165 
1166  ! sort by ascending order
1167  call sort_exec( npoints, & ! [IN]
1168  drad(:), & ! [INOUT]
1169  idx_i(:) ) ! [INOUT]
1170  endif
1171 
1172  return
1173  end subroutine interp_insert
1174 
1175 !OCL SERIAL
1176  subroutine interp_insert_2d( npoints, &
1177  lon, lat, &
1178  lon_ref, lat_ref, &
1179  i, j, &
1180  drad, idx_i, idx_j )
1181  use scale_const, only: &
1182  undef => const_undef
1183  use scale_sort, only: &
1184  sort_exec
1185 
1186  integer, intent(in) :: npoints
1187  real(RP), intent(in) :: lon, lat
1188  real(RP), intent(in) :: lon_ref, lat_ref
1189  integer, intent(in) :: i, j
1190 
1191  real(RP), intent(inout) :: drad(npoints)
1192  integer, intent(inout) :: idx_i(npoints)
1193  integer, intent(inout) :: idx_j(npoints)
1194 
1195  real(RP) :: dradian
1196 
1197  if ( lon_ref == undef ) return
1198 
1199  dradian = haversine( lon, lat, lon_ref, lat_ref )
1200 
1201  if ( dradian <= drad(npoints) ) then
1202  ! replace last(=longest) value
1203  drad(npoints) = dradian
1204  idx_i(npoints) = i
1205  idx_j(npoints) = j
1206 
1207  ! sort by ascending order
1208  call sort_exec( npoints, & ! [IN]
1209  drad(:), & ! [INOUT]
1210  idx_i(:), & ! [INOUT]
1211  idx_j(:) ) ! [INOUT]
1212  endif
1213 
1214  return
1215  end subroutine interp_insert_2d
1216 
1217  subroutine interp_div_block(nsize, psize, nidx_max, &
1218  lon_ref, lat_ref, &
1219  idx, nidx, &
1220  lon_min, lon_max, &
1221  lat_min, lat_max, &
1222  dlon, dlat )
1223  use scale_const, only: &
1224  undef => const_undef
1225  use scale_prc, only: &
1226  prc_abort
1227  integer, intent(in) :: nsize
1228  integer, intent(in) :: psize
1229  integer, intent(in) :: nidx_max
1230  real(RP), intent(in) :: lon_ref(nsize)
1231  real(RP), intent(in) :: lat_ref(nsize)
1232 
1233  integer, intent(out) :: idx (nidx_max,psize,psize)
1234  integer, intent(out) :: nidx(psize,psize)
1235  real(RP), intent(out) :: lon_min, lon_max
1236  real(RP), intent(out) :: lat_min, lat_max
1237  real(RP), intent(out) :: dlon, dlat
1238 
1239  integer :: i, ii, jj, n
1240 
1241  lon_min = minval(lon_ref(:), mask=lon_ref.ne.undef)
1242  lon_max = maxval(lon_ref(:), mask=lon_ref.ne.undef)
1243  lat_min = minval(lat_ref(:), mask=lat_ref.ne.undef)
1244  lat_max = maxval(lat_ref(:), mask=lat_ref.ne.undef)
1245 
1246  dlon = ( lon_max - lon_min ) / psize
1247  dlat = ( lat_max - lat_min ) / psize
1248 
1249  nidx(:,:) = 0
1250  do i = 1, nsize
1251  if ( lon_ref(i) == undef ) cycle
1252  ii = min(int((lon_ref(i) - lon_min) / dlon) + 1, psize)
1253  jj = min(int((lat_ref(i) - lat_min) / dlat) + 1, psize)
1254  n = nidx(ii,jj) + 1
1255  nidx(ii,jj) = n
1256  if ( n <= nidx_max ) idx(n,ii,jj) = i
1257  end do
1258 
1259  if ( maxval(nidx) > nidx_max ) then
1260  log_error("INTERP_search_horiz",*) 'Buffer size is not enough'
1261  log_error_cont(*) ' Use larger INTERP_buffer_size_fact: ', interp_buffer_size_fact * maxval(nidx)/nidx_max
1262  call prc_abort
1263  end if
1264 
1265  return
1266  end subroutine interp_div_block
1267 
1268  !-----------------------------------------------------------------------------
1269  ! Haversine Formula (from R.W. Sinnott, "Virtues of the Haversine",
1270  ! Sky and Telescope, vol. 68, no. 2, 1984, p. 159):
1271  function haversine( &
1272  lon0, lat0, &
1273  lon1, lat1 ) &
1274  result( d )
1275  implicit none
1276 
1277  real(RP), intent(in) :: lon0 ! [rad]
1278  real(RP), intent(in) :: lat0 ! [rad]
1279  real(RP), intent(in) :: lon1 ! [rad]
1280  real(RP), intent(in) :: lat1 ! [rad]
1281  real(RP) :: d ! [rad]
1282 
1283  real(RP) :: dlonh, dlath
1284  real(RP) :: work1
1285  !---------------------------------------------------------------------------
1286 
1287  dlonh = 0.5_rp * ( lon0 - lon1 )
1288  dlath = 0.5_rp * ( lat0 - lat1 )
1289 
1290  work1 = sin(dlath)**2 + cos(lat0) * cos(lat1) * sin(dlonh)**2
1291 
1292  d = 2.0_rp * asin( min(sqrt(work1),1.0_rp) )
1293 
1294  end function haversine
1295 
1296 end module scale_interp
scale_debug
module DEBUG
Definition: scale_debug.F90:11
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:1223
scale_const::const_radius
real(rp), public const_radius
radius of the planet [m]
Definition: scale_const.F90:44
scale_interp
module INTERPOLATION
Definition: scale_interp.F90:12
scale_interp::interp_factor3d
subroutine, public interp_factor3d(npoints, KA_ref, KS_ref, KE_ref, IA_ref, JA_ref, lon_ref, lat_ref, hgt_ref, KA, KS, KE, IA, JA, lon, lat, hgt, idx_i, idx_j, hfact, idx_k, vfact)
Definition: scale_interp.F90:484
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:55
scale_const::const_d2r
real(rp), public const_d2r
degree to radian
Definition: scale_const.F90:32
scale_interp::interp_interp3d
subroutine, public interp_interp3d(npoints, KA_ref, IA_ref, JA_ref, KA, KS, KE, IA, JA, idx_i, idx_j, hfact, idx_k, vfact, val_ref, val, logwgt)
Definition: scale_interp.F90:658
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:41
scale_tracer
module TRACER
Definition: scale_tracer.F90:16
scale_index
module Index
Definition: scale_index.F90:11
scale_atmos_grid_cartesc_index
module atmosphere / grid / cartesC index
Definition: scale_atmos_grid_cartesC_index.F90:12
scale_prc
module PROCESS
Definition: scale_prc.F90:11
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:119
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:1181
scale_interp::interp_interp2d
subroutine, public interp_interp2d(npoints, IA_ref, JA_ref, IA, JA, idx_i, idx_j, hfact, val_ref, val)
Definition: scale_interp.F90:594
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:338
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_interp::interp_setup
subroutine, public interp_setup(weight_order, search_limit)
Setup.
Definition: scale_interp.F90:69
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:921
scale_prof::prof_rapstart
subroutine, public prof_rapstart(rapname_base, level)
Start raptime.
Definition: scale_prof.F90:157
scale_prof
module profiler
Definition: scale_prof.F90:11
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:33
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_interp::interp_factor2d
subroutine, public interp_factor2d(npoints, IA_ref, JA_ref, lon_ref, lat_ref, IA, JA, lon, lat, idx_i, idx_j, hfact, search_limit, latlon_structure, weight_order)
Definition: scale_interp.F90:226
scale_io
module STDIO
Definition: scale_io.F90:10
scale_prof::prof_rapend
subroutine, public prof_rapend(rapname_base, level)
Save raptime.
Definition: scale_prof.F90:210
scale_precision::rp
integer, parameter, public rp
Definition: scale_precision.F90:41
scale_sort
module SORT
Definition: scale_sort.F90:11
Generated by   doxygen 1.8.13