SCALE-RM
scale_interpolation_nest.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
13 !-------------------------------------------------------------------------------
14 module scale_interpolation_nest
15  !-----------------------------------------------------------------------------
16  !
17  !++ used modules
18  !
19  use scale_precision
20  use scale_stdio
21  use scale_prof
22  use scale_debug
23  use scale_grid_index
24  use scale_land_grid_index
25  use scale_index
26  use scale_tracer
27  use scale_const, only: &
28  r_in_m => const_radius
29  !-----------------------------------------------------------------------------
30  implicit none
31  private
32  !-----------------------------------------------------------------------------
33  !
34  !++ Public procedure
35  !
36  public :: intrpnest_setup
37  public :: intrpnest_domain_compatibility
38  public :: intrpnest_interp_fact_latlon
39  public :: intrpnest_interp_fact_llz
40 
41  !-----------------------------------------------------------------------------
42  !
43  !++ Public parameters & variables
44  !
45  !-----------------------------------------------------------------------------
46  !
47  !++ Private procedure
48  !
49  private :: intrpnest_search_nearest_block
50  private :: intrpnest_search_horiz_1points
51  private :: intrpnest_search_horiz_3points
52  private :: intrpnest_search_horiz_4points
53  private :: intrpnest_search_horiz_5points
54  private :: intrpnest_search_horiz_8points
55  private :: intrpnest_search_horiz_12points
56  private :: intrpnest_search_vert_offline
57  private :: intrpnest_search_vert_online
58  private :: intrpnest_interp_2d_1points
59  private :: intrpnest_interp_3d_1points
60  private :: intrpnest_interp_2d_3points
61  private :: intrpnest_interp_3d_3points
62  private :: intrpnest_interp_2d_4points
63  private :: intrpnest_interp_3d_4points
64  private :: intrpnest_interp_2d_5points
65  private :: intrpnest_interp_3d_5points
66  private :: intrpnest_interp_2d_8points
67  private :: intrpnest_interp_3d_8points
68  private :: intrpnest_interp_2d_12points
69  private :: intrpnest_interp_3d_12points
70  private :: intrpnest_haversine
71 
72  !-----------------------------------------------------------------------------
73  abstract interface
74  subroutine intrpnest_intfc_search_h( &
75  hfact, & ! (out)
76  igrd, & ! (out)
77  jgrd, & ! (out)
78  mylat, & ! (in)
79  mylon, & ! (in)
80  inlat, & ! (in)
81  inlon, & ! (in)
82  is, & ! (in)
83  ie, & ! (in)
84  js, & ! (in)
85  je ) ! (in)
86  use scale_precision
87  implicit none
88 
89  real(RP), intent(out) :: hfact(:) ! horizontal interp factor
90  integer, intent(out) :: igrd (:) ! grid points of interp target
91  integer, intent(out) :: jgrd (:) ! grid points of interp target
92  real(RP), intent(in) :: mylat ! latitude data of mine
93  real(RP), intent(in) :: mylon ! longitude data of mine
94  real(RP), intent(in) :: inlat(:,:) ! latitude data of you (input)
95  real(RP), intent(in) :: inlon(:,:) ! longitude data of you (input)
96  integer, intent(in) :: is ! start index for x-direction
97  integer, intent(in) :: ie ! end index for x-direction
98  integer, intent(in) :: js ! start index for y-direction
99  integer, intent(in) :: je ! end index for y-direction
100  end subroutine intrpnest_intfc_search_h
101  end interface
102  procedure(INTRPNEST_intfc_search_h), pointer :: INTRPNEST_search_horiz => null()
103  private :: intrpnest_search_horiz
104 
105  !-----------------------------------------------------------------------------
106  abstract interface
107  subroutine intrpnest_intfc_search_v( &
108  vfact, & ! (out)
109  kgrd, & ! (out)
110  ncopy, & ! (out)
111  igrd, & ! (in)
112  jgrd, & ! (in)
113  myhgt, & ! (in)
114  inhgt, & ! (in)
115  iloc, & ! (in)
116  jloc, & ! (in)
117  ks, & ! (in)
118  ke, & ! (in)
119  inka, & ! (in)
120  lndgrd ) ! (in)
121  use scale_precision
122  implicit none
123 
124  real(RP), intent(inout) :: vfact(:,:,:,:,:)! vertical interp factor
125  integer, intent(inout) :: kgrd (:,:,:,:,:)! grid points of interp target
126  integer, intent(out) :: ncopy(:) ! number of daughter's layers below the parent's lowest layer
127  integer, intent(in) :: igrd(:) ! grid points of interp target
128  integer, intent(in) :: jgrd(:) ! grid points of interp target
129  real(RP), intent(in) :: myhgt(:) ! height data of mine
130  real(RP), intent(in) :: inhgt(:,:,:) ! height data of you (input)
131  integer, intent(in) :: iloc ! locator index for x-direction
132  integer, intent(in) :: jloc ! locator index for y-direction
133  integer, intent(in) :: ks ! start index for z-direction
134  integer, intent(in) :: ke ! end index for z-direction
135  integer, intent(in) :: inKA ! grid number of you (input)
136  logical, intent(in) :: lndgrd ! flag of land grid
137  end subroutine intrpnest_intfc_search_v
138  end interface
139  procedure(INTRPNEST_intfc_search_v), pointer :: INTRPNEST_search_vert => null()
140  private :: intrpnest_search_vert
141 
142  !-----------------------------------------------------------------------------
143  abstract interface
144  subroutine intrpnest_intfc_interp_2d( &
145  intp, & ! (out)
146  ref, & ! (in)
147  hfact, & ! (in)
148  igrd, & ! (in)
149  jgrd, & ! (in)
150  ia, & ! (in)
151  ja ) ! (in)
152  use scale_precision
153  implicit none
154 
155  real(RP), intent(out) :: intp(:,:) ! interpolated data
156  real(RP), intent(in) :: ref (:,:) ! reference data
157  real(RP), intent(in) :: hfact(:,:,:) ! horizontal interp factor
158  integer, intent(in) :: igrd (:,:,:) ! grid points of interp target
159  integer, intent(in) :: jgrd (:,:,:) ! grid points of interp target
160  integer, intent(in) :: ia ! grid number of mine
161  integer, intent(in) :: ja ! grid number of mine
162  end subroutine intrpnest_intfc_interp_2d
163  end interface
164  procedure(intrpnest_intfc_interp_2d), pointer :: intrpnest_interp_2d => null()
165  public :: intrpnest_interp_2d
166 
167  !-----------------------------------------------------------------------------
168  abstract interface
169  subroutine intrpnest_intfc_interp_3d( &
170  intp, & ! (out)
171  ref, & ! (in)
172  hfact, & ! (in)
173  vfact, & ! (in)
174  kgrd, & ! (in)
175  igrd, & ! (in)
176  jgrd, & ! (in)
177  ia, & ! (in)
178  ja, & ! (in)
179  ks, & ! (in)
180  ke, & ! (in)
181  logwegt ) ! (in)
182  use scale_precision
183  implicit none
184 
185  real(RP), intent(out) :: intp(:,:,:) ! interpolated data
186  real(RP), intent(in) :: ref (:,:,:) ! reference data
187  real(RP), intent(in) :: hfact(:,:,:) ! horizontal interp factor
188  real(RP), intent(in) :: vfact(:,:,:,:,:) ! vertical interp factor
189  integer, intent(in) :: kgrd (:,:,:,:,:) ! grid points of interp target
190  integer, intent(in) :: igrd (:,:,:) ! grid points of interp target
191  integer, intent(in) :: jgrd (:,:,:) ! grid points of interp target
192  integer, intent(in) :: ia ! grid number of mine
193  integer, intent(in) :: ja ! grid number of mine
194  integer, intent(in) :: ks ! start grid number of mine
195  integer, intent(in) :: ke ! end grid number of mine
196  logical, intent(in), optional :: logwegt
197  end subroutine intrpnest_intfc_interp_3d
198  end interface
199  procedure(intrpnest_intfc_interp_3d), pointer :: intrpnest_interp_3d => null()
200  public :: intrpnest_interp_3d
201 
202  !-----------------------------------------------------------------------------
203  !
204  !++ Private parameters & variables
205  !
206  real(RP), private, parameter :: large_number_1 = 9.999e+15_rp
207  real(RP), private, parameter :: large_number_2 = 9.888e+15_rp
208  real(RP), private, parameter :: large_number_3 = 9.777e+15_rp
209  real(RP), private, parameter :: large_number_4 = 9.666e+15_rp
210  real(RP), private, parameter :: large_number_5 = 9.555e+15_rp
211  real(RP), private, parameter :: large_number_6 = 9.444e+15_rp
212  real(RP), private, parameter :: large_number_7 = 9.333e+15_rp
213  real(RP), private, parameter :: large_number_8 = 9.222e+15_rp
214  real(RP), private, parameter :: large_number_9 = 9.111e+15_rp
215  real(RP), private, parameter :: large_number_10 = 9.000e+15_rp
216  real(RP), private, parameter :: large_number_11 = 8.999e+15_rp
217  real(RP), private, parameter :: large_number_12 = 8.888e+15_rp
218 
219  integer, private :: weight_order
220 
221  integer, private :: divnum
222  integer, private :: itp_nh
223 
224  !-----------------------------------------------------------------------------
225 contains
226  !-----------------------------------------------------------------------------
228  subroutine intrpnest_setup( &
229  interp_search_divnum, &
230  NEST_INTERP_LEVEL, &
231  NEST_INTERP_WEIGHT_ORDER, &
232  OFFLINE )
233  use scale_process, only: &
234  prc_mpistop
235  implicit none
236 
237  integer, intent(in) :: interp_search_divnum
238  integer, intent(in) :: nest_interp_level
239  integer, intent(in) :: nest_interp_weight_order
240  logical, intent(in) :: offline
241 
242  character(len=7) :: select_type
243  !---------------------------------------------------------------------------
244 
245  if( io_l ) write(io_fid_log,*)
246  if( io_l ) write(io_fid_log,*) '++++++ Module[INTERP NEST] / Categ[ATMOS-RM GRID] / Origin[SCALElib]'
247 
248  divnum = interp_search_divnum
249  weight_order = nest_interp_weight_order
250 
251  select case( nest_interp_level )
252  case( 1 )
253  intrpnest_search_horiz => intrpnest_search_horiz_1points
254  intrpnest_interp_2d => intrpnest_interp_2d_1points
255  intrpnest_interp_3d => intrpnest_interp_3d_1points
256  itp_nh = 1
257 
258  case( 3 )
259  intrpnest_search_horiz => intrpnest_search_horiz_3points
260  intrpnest_interp_2d => intrpnest_interp_2d_3points
261  intrpnest_interp_3d => intrpnest_interp_3d_3points
262  itp_nh = 3
263 
264  case( 4 )
265  intrpnest_search_horiz => intrpnest_search_horiz_4points
266  intrpnest_interp_2d => intrpnest_interp_2d_4points
267  intrpnest_interp_3d => intrpnest_interp_3d_4points
268  itp_nh = 4
269 
270  case( 5 )
271  intrpnest_search_horiz => intrpnest_search_horiz_5points
272  intrpnest_interp_2d => intrpnest_interp_2d_5points
273  intrpnest_interp_3d => intrpnest_interp_3d_5points
274  itp_nh = 5
275 
276  case( 8 )
277  intrpnest_search_horiz => intrpnest_search_horiz_8points
278  intrpnest_interp_2d => intrpnest_interp_2d_8points
279  intrpnest_interp_3d => intrpnest_interp_3d_8points
280  itp_nh = 8
281 
282  case( 12 )
283  intrpnest_search_horiz => intrpnest_search_horiz_12points
284  intrpnest_interp_2d => intrpnest_interp_2d_12points
285  intrpnest_interp_3d => intrpnest_interp_3d_12points
286  itp_nh = 12
287 
288  case default
289  write(*,*) 'xxx invarid NEST_INTERP_LEVEL (', nest_interp_level, &
290  ') [setup: nest/interp]'
291  call prc_mpistop
292  end select
293 
294  if ( offline ) then
295  select_type = "offline"
296  intrpnest_search_vert => intrpnest_search_vert_offline
297  else
298  select_type = "online"
299  intrpnest_search_vert => intrpnest_search_vert_online
300  endif
301 
302  if( io_l ) write(io_fid_log,*)
303  if( io_l ) write(io_fid_log,*) '*** Horizontal interpolation with ', nest_interp_level, " points."
304  if( io_l ) write(io_fid_log,*) '*** Vertical interpolation : ', trim(select_type)
305 
306  return
307  end subroutine intrpnest_setup
308 
309 
310  !-----------------------------------------------------------------------------
311  ! make interpolation factor using Lat-Lon
312  subroutine intrpnest_interp_fact_latlon( &
313  hfact, & ! (out)
314  igrd, & ! (out)
315  jgrd, & ! (out)
316  mylat, & ! (in)
317  mylon, & ! (in)
318  myia, & ! (in)
319  myja, & ! (in)
320  inlat, & ! (in)
321  inlon, & ! (in)
322  inia, & ! (in)
323  inja ) ! (in)
324  use scale_process, only: &
325  prc_mpistop
326  implicit none
327 
328  real(RP), intent(out) :: hfact(:,:,:) ! horizontal interp factor
329  integer, intent(out) :: igrd (:,:,:) ! grid points of interp target
330  integer, intent(out) :: jgrd (:,:,:) ! grid points of interp target
331 
332  real(RP), intent(in) :: mylat(:,:) ! latitude data of mine
333  real(RP), intent(in) :: mylon(:,:) ! longitude data of mine
334  integer, intent(in) :: myia ! grid number of mine
335  integer, intent(in) :: myja ! grid number of mine
336 
337  real(RP), intent(in) :: inlat(:,:) ! latitude data of you (input)
338  real(RP), intent(in) :: inlon(:,:) ! longitude data of you (input)
339  integer, intent(in) :: inia ! grid number of you (input)
340  integer, intent(in) :: inja ! grid number of you (input)
341 
342  integer :: i, j
343  integer :: is, ie
344  integer :: js, je
345  !---------------------------------------------------------------------------
346 
347  hfact(:,:,:) = 0.0_rp
348 
349  do j = 1, myja
350  do i = 1, myia
351  ! nearest block search
352  call intrpnest_search_nearest_block( is, ie, js, je, &
353  mylat(i,j), mylon(i,j), &
354  inlat(:,:), inlon(:,:), &
355  inia, inja )
356 
357  ! main search
358  call intrpnest_search_horiz( hfact(i,j,:), &
359  igrd(i,j,:), &
360  jgrd(i,j,:), &
361  mylat(i,j), &
362  mylon(i,j), &
363  inlat, &
364  inlon, &
365  is, ie, &
366  js, je )
367  enddo
368  enddo
369 
370  return
371  end subroutine intrpnest_interp_fact_latlon
372 
373 
374  !-----------------------------------------------------------------------------
375  ! make interpolation factor using Lat-Lon and Z-Height information
376  subroutine intrpnest_interp_fact_llz( &
377  hfact, & ! (out)
378  vfact, & ! (out)
379  kgrd, & ! (out)
380  igrd, & ! (out)
381  jgrd, & ! (out)
382  ncopy, & ! (out)
383  myhgt, & ! (in)
384  mylat, & ! (in)
385  mylon, & ! (in)
386  myks, & ! (in)
387  myke, & ! (in)
388  myia, & ! (in)
389  myja, & ! (in)
390  inhgt, & ! (in)
391  inlat, & ! (in)
392  inlon, & ! (in)
393  inka, & ! (in)
394  inia, & ! (in)
395  inja, & ! (in)
396  landgrid ) ! (in)
397  use scale_process, only: &
398  prc_mpistop
399  implicit none
400 
401  real(RP), intent(out) :: hfact(:,:,:) ! horizontal interp factor
402  real(RP), intent(out) :: vfact(:,:,:,:,:) ! vertical interp factor
403  integer, intent(out) :: kgrd (:,:,:,:,:) ! grid points of interp target
404  integer, intent(out) :: igrd (:,:,:) ! grid points of interp target
405  integer, intent(out) :: jgrd (:,:,:) ! grid points of interp target
406  integer, intent(out) :: ncopy(:,:,:) ! number of daughter's layers below parent lowest layer
407 
408  real(RP), intent(in) :: myhgt(:,:,:) ! height data of mine
409  real(RP), intent(in) :: mylat(:,:) ! latitude data of mine
410  real(RP), intent(in) :: mylon(:,:) ! longitude data of mine
411  integer, intent(in) :: myks ! start grid number of mine
412  integer, intent(in) :: myke ! end grid number of mine
413  integer, intent(in) :: myia ! grid number of mine
414  integer, intent(in) :: myja ! grid number of mine
415 
416  real(RP), intent(in) :: inhgt(:,:,:) ! height data of you (input)
417  real(RP), intent(in) :: inlat(:,:) ! latitude data of you (input)
418  real(RP), intent(in) :: inlon(:,:) ! longitude data of you (input)
419  integer, intent(in) :: inka ! grid number of you (input)
420  integer, intent(in) :: inia ! grid number of you (input)
421  integer, intent(in) :: inja ! grid number of you (input)
422 
423  logical, intent(in), optional :: landgrid
424 
425  integer :: i, j
426  integer :: is, ie
427  integer :: js, je
428  logical :: lndgrd
429  !---------------------------------------------------------------------------
430 
431  lndgrd = .false.
432  if ( present(landgrid) ) then
433  if ( landgrid ) then
434  lndgrd = .true.
435  endif
436  endif
437 
438  hfact(:,:,:) = 0.0_rp
439  vfact(:,:,:,:,:) = 0.0_rp
440  ncopy(:,:,:) = 0
441 
442  do j = 1, myja
443  do i = 1, myia
444 
445  ! nearest block search
446  call intrpnest_search_nearest_block( is, ie, js, je, &
447  mylat(i,j), mylon(i,j), &
448  inlat(:,:), inlon(:,:), &
449  inia, inja )
450 
451  ! main search
452  call intrpnest_search_horiz( hfact(i,j,:), &
453  igrd(i,j,:), &
454  jgrd(i,j,:), &
455  mylat(i,j), &
456  mylon(i,j), &
457  inlat, &
458  inlon, &
459  is, ie, &
460  js, je )
461 
462 
463  call intrpnest_search_vert( vfact, &
464  kgrd, &
465  ncopy(i,j,:), &
466  igrd(i,j,:), &
467  jgrd(i,j,:), &
468  myhgt(:,i,j), &
469  inhgt, &
470  i, j, &
471  myks, myke, &
472  inka, &
473  lndgrd )
474 
475  enddo
476  enddo
477 
478  return
479  end subroutine intrpnest_interp_fact_llz
480 
481 
482  !-----------------------------------------------------------------------------
483  ! search of nearest region for speed up of interpolation
484  subroutine intrpnest_search_nearest_block( &
485  is, & ! (out)
486  ie, & ! (out)
487  js, & ! (out)
488  je, & ! (out)
489  mylat, & ! (in)
490  mylon, & ! (in)
491  inlat, & ! (in)
492  inlon, & ! (in)
493  inia, & ! (in)
494  inja ) ! (in)
495  implicit none
496 
497  integer, intent(out) :: is ! start index for x-direction
498  integer, intent(out) :: ie ! end index for x-direction
499  integer, intent(out) :: js ! start index for y-direction
500  integer, intent(out) :: je ! end index for y-direction
501 
502  real(RP), intent(in) :: mylat ! latitude data of mine
503  real(RP), intent(in) :: mylon ! longitude data of mine
504  real(RP), intent(in) :: inlat(:,:) ! latitude data of you (input)
505  real(RP), intent(in) :: inlon(:,:) ! longitude data of you (input)
506  integer, intent(in) :: inia ! grid number of you (input)
507  integer, intent(in) :: inja ! grid number of you (input)
508 
509  real(RP) :: distance, dist
510  integer :: ii, jj
511  integer :: iinc, jinc
512  integer :: blk_i, blk_j
513  !---------------------------------------------------------------------------
514 
515  iinc = max( (inia + 1) / divnum, 1 )
516  jinc = max( (inja + 1) / divnum, 1 )
517  dist = large_number_1
518 
519  jj = 1 + (jinc/2)
520  do while (jj <= inja)
521  ii = 1 + (iinc/2)
522  do while (ii <= inia)
523  distance = intrpnest_haversine( mylat, mylon, &
524  inlat(ii,jj), inlon(ii,jj) )
525  if( distance < dist )then
526  dist = distance
527  blk_i = ii
528  blk_j = jj
529  endif
530  ii = ii + iinc
531  enddo
532  jj = jj + jinc
533  enddo
534 
535  ! +- 3 is buffer for 12 points
536  is = blk_i - (iinc/2) - 3
537  if( is < 1 ) is = 1
538  ie = blk_i + (iinc/2) + 3
539  if( ie > inia ) ie = inia
540  js = blk_j - (jinc/2) - 3
541  if( js < 1 ) js = 1
542  je = blk_j + (jinc/2) + 3
543  if( je > inja ) je = inja
544 
545  return
546  end subroutine intrpnest_search_nearest_block
547 
548 
549  !-----------------------------------------------------------------------------
550  ! horizontal search of interpolation points for one-points (nearest-neighbor)
551  subroutine intrpnest_search_horiz_1points( &
552  hfact, & ! (out)
553  igrd, & ! (out)
554  jgrd, & ! (out)
555  mylat, & ! (in)
556  mylon, & ! (in)
557  inlat, & ! (in)
558  inlon, & ! (in)
559  is, & ! (in)
560  ie, & ! (in)
561  js, & ! (in)
562  je ) ! (in)
563  implicit none
564 
565  real(RP), intent(out) :: hfact(:) ! horizontal interp factor
566  integer, intent(out) :: igrd (:) ! grid points of interp target
567  integer, intent(out) :: jgrd (:) ! grid points of interp target
568 
569  real(RP), intent(in) :: mylat ! latitude data of mine
570  real(RP), intent(in) :: mylon ! longitude data of mine
571  real(RP), intent(in) :: inlat(:,:) ! latitude data of you (input)
572  real(RP), intent(in) :: inlon(:,:) ! longitude data of you (input)
573 
574  integer, intent(in) :: is ! start index for x-direction
575  integer, intent(in) :: ie ! end index for x-direction
576  integer, intent(in) :: js ! start index for y-direction
577  integer, intent(in) :: je ! end index for y-direction
578 
579  real(RP) :: distance
580  real(RP) :: dist
581  integer :: ii, jj
582  !---------------------------------------------------------------------------
583 
584  dist = large_number_1
585  igrd(:) = -1
586  jgrd(:) = -1
587 
588  do jj = js, je
589  do ii = is, ie
590  distance = intrpnest_haversine( mylat,mylon,inlat(ii,jj),inlon(ii,jj) )
591  if ( distance <= dist ) then
592  dist = distance; igrd(1) = ii; jgrd(1) = jj
593  endif
594  enddo
595  enddo
596 
597  hfact(1) = 1.0_rp
598 
599  return
600  end subroutine intrpnest_search_horiz_1points
601 
602 
603  !-----------------------------------------------------------------------------
604  ! horizontal search of interpolation points for three-points
605  subroutine intrpnest_search_horiz_3points( &
606  hfact, & ! (out)
607  igrd, & ! (out)
608  jgrd, & ! (out)
609  mylat, & ! (in)
610  mylon, & ! (in)
611  inlat, & ! (in)
612  inlon, & ! (in)
613  is, & ! (in)
614  ie, & ! (in)
615  js, & ! (in)
616  je ) ! (in)
617  use scale_const, only: &
618  eps => const_eps
619  implicit none
620 
621  real(RP), intent(out) :: hfact(:) ! horizontal interp factor
622  integer, intent(out) :: igrd (:) ! grid points of interp target
623  integer, intent(out) :: jgrd (:) ! grid points of interp target
624 
625  real(RP), intent(in) :: mylat ! latitude data of mine
626  real(RP), intent(in) :: mylon ! longitude data of mine
627  real(RP), intent(in) :: inlat(:,:) ! latitude data of you (input)
628  real(RP), intent(in) :: inlon(:,:) ! longitude data of you (input)
629 
630  integer, intent(in) :: is ! start index for x-direction
631  integer, intent(in) :: ie ! end index for x-direction
632  integer, intent(in) :: js ! start index for y-direction
633  integer, intent(in) :: je ! end index for y-direction
634 
635  real(RP) :: distance
636  real(RP) :: denom
637  real(RP) :: dist(3)
638  integer :: ii, jj
639  !---------------------------------------------------------------------------
640 
641  dist(1) = large_number_3
642  dist(2) = large_number_2
643  dist(3) = large_number_1
644  igrd(:) = -1
645  jgrd(:) = -1
646 
647  do jj = js, je
648  do ii = is, ie
649  distance = intrpnest_haversine( mylat,mylon,inlat(ii,jj),inlon(ii,jj) )
650  distance = distance**weight_order
651  if ( distance <= dist(1) ) then
652  dist(3) = dist(2); igrd(3) = igrd(2); jgrd(3) = jgrd(2)
653  dist(2) = dist(1); igrd(2) = igrd(1); jgrd(2) = jgrd(1)
654  dist(1) = distance; igrd(1) = ii; jgrd(1) = jj
655  elseif ( dist(1) < distance .AND. distance <= dist(2) ) then
656  dist(3) = dist(2); igrd(3) = igrd(2); jgrd(3) = jgrd(2)
657  dist(2) = distance; igrd(2) = ii; jgrd(2) = jj
658  elseif ( dist(2) < distance .AND. distance <= dist(3) ) then
659  dist(3) = distance; igrd(3) = ii; jgrd(3) = jj
660  endif
661  enddo
662  enddo
663 
664  if ( abs(dist(1)) < eps ) then
665  hfact(1) = 1.0_rp
666  hfact(2) = 0.0_rp
667  hfact(3) = 0.0_rp
668  else
669  denom = 1.0_rp / ( (1.0_rp/dist(1)) + (1.0_rp/dist(2)) + (1.0_rp/dist(3)) )
670  hfact(1) = ( 1.0_rp/dist(1) ) * denom
671  hfact(2) = ( 1.0_rp/dist(2) ) * denom
672  hfact(3) = ( 1.0_rp/dist(3) ) * denom
673  endif
674 
675  return
676  end subroutine intrpnest_search_horiz_3points
677 
678 
679  !-----------------------------------------------------------------------------
680  ! horizontal search of interpolation points for four-points
681  subroutine intrpnest_search_horiz_4points( &
682  hfact, & ! (out)
683  igrd, & ! (out)
684  jgrd, & ! (out)
685  mylat, & ! (in)
686  mylon, & ! (in)
687  inlat, & ! (in)
688  inlon, & ! (in)
689  is, & ! (in)
690  ie, & ! (in)
691  js, & ! (in)
692  je ) ! (in)
693  use scale_const, only: &
694  eps => const_eps
695  implicit none
696 
697  real(RP), intent(out) :: hfact(:) ! horizontal interp factor
698  integer, intent(out) :: igrd (:) ! grid points of interp target
699  integer, intent(out) :: jgrd (:) ! grid points of interp target
700 
701  real(RP), intent(in) :: mylat ! latitude data of mine
702  real(RP), intent(in) :: mylon ! longitude data of mine
703  real(RP), intent(in) :: inlat(:,:) ! latitude data of you (input)
704  real(RP), intent(in) :: inlon(:,:) ! longitude data of you (input)
705 
706  integer, intent(in) :: is ! start index for x-direction
707  integer, intent(in) :: ie ! end index for x-direction
708  integer, intent(in) :: js ! start index for y-direction
709  integer, intent(in) :: je ! end index for y-direction
710 
711  real(RP) :: distance
712  real(RP) :: denom
713  real(RP) :: dist(4)
714  integer :: ii, jj
715  !---------------------------------------------------------------------------
716 
717  dist(1) = large_number_4
718  dist(2) = large_number_3
719  dist(3) = large_number_2
720  dist(4) = large_number_1
721  igrd(:) = -1
722  jgrd(:) = -1
723 
724  do jj = js, je
725  do ii = is, ie
726  distance = intrpnest_haversine( mylat,mylon,inlat(ii,jj),inlon(ii,jj) )
727  distance = distance**weight_order
728  if ( distance <= dist(1) ) then
729  dist(4) = dist(3); igrd(4) = igrd(3); jgrd(4) = jgrd(3)
730  dist(3) = dist(2); igrd(3) = igrd(2); jgrd(3) = jgrd(2)
731  dist(2) = dist(1); igrd(2) = igrd(1); jgrd(2) = jgrd(1)
732  dist(1) = distance; igrd(1) = ii; jgrd(1) = jj
733  elseif ( dist(1) < distance .AND. distance <= dist(2) ) then
734  dist(4) = dist(3); igrd(4) = igrd(3); jgrd(4) = jgrd(3)
735  dist(3) = dist(2); igrd(3) = igrd(2); jgrd(3) = jgrd(2)
736  dist(2) = distance; igrd(2) = ii; jgrd(2) = jj
737  elseif ( dist(2) < distance .AND. distance <= dist(3) ) then
738  dist(4) = dist(3); igrd(4) = igrd(3); jgrd(4) = jgrd(3)
739  dist(3) = distance; igrd(3) = ii; jgrd(3) = jj
740  elseif ( dist(3) < distance .AND. distance <= dist(4) ) then
741  dist(4) = distance; igrd(4) = ii; jgrd(4) = jj
742  endif
743  enddo
744  enddo
745 
746  if ( abs(dist(1)) < eps ) then
747  hfact(1) = 1.0_rp
748  hfact(2) = 0.0_rp
749  hfact(3) = 0.0_rp
750  hfact(4) = 0.0_rp
751  else
752  denom = 1.0_rp / ( (1.0_rp/dist(1)) + (1.0_rp/dist(2)) &
753  + (1.0_rp/dist(3)) + (1.0_rp/dist(4)) )
754  hfact(1) = ( 1.0_rp/dist(1) ) * denom
755  hfact(2) = ( 1.0_rp/dist(2) ) * denom
756  hfact(3) = ( 1.0_rp/dist(3) ) * denom
757  hfact(4) = ( 1.0_rp/dist(4) ) * denom
758  endif
759 
760  return
761  end subroutine intrpnest_search_horiz_4points
762 
763 
764  !-----------------------------------------------------------------------------
765  ! horizontal search of interpolation points for five-points
766  subroutine intrpnest_search_horiz_5points( &
767  hfact, & ! (out)
768  igrd, & ! (out)
769  jgrd, & ! (out)
770  mylat, & ! (in)
771  mylon, & ! (in)
772  inlat, & ! (in)
773  inlon, & ! (in)
774  is, & ! (in)
775  ie, & ! (in)
776  js, & ! (in)
777  je ) ! (in)
778  use scale_const, only: &
779  eps => const_eps
780  implicit none
781 
782  real(RP), intent(out) :: hfact(:) ! horizontal interp factor
783  integer, intent(out) :: igrd (:) ! grid points of interp target
784  integer, intent(out) :: jgrd (:) ! grid points of interp target
785 
786  real(RP), intent(in) :: mylat ! latitude data of mine
787  real(RP), intent(in) :: mylon ! longitude data of mine
788  real(RP), intent(in) :: inlat(:,:) ! latitude data of you (input)
789  real(RP), intent(in) :: inlon(:,:) ! longitude data of you (input)
790 
791  integer, intent(in) :: is ! start index for x-direction
792  integer, intent(in) :: ie ! end index for x-direction
793  integer, intent(in) :: js ! start index for y-direction
794  integer, intent(in) :: je ! end index for y-direction
795 
796  real(RP) :: distance
797  real(RP) :: denom
798  real(RP) :: dist(5)
799  integer :: ii, jj
800  !---------------------------------------------------------------------------
801 
802  dist(1) = large_number_5
803  dist(2) = large_number_4
804  dist(3) = large_number_3
805  dist(4) = large_number_2
806  dist(5) = large_number_1
807  igrd(:) = -1
808  jgrd(:) = -1
809 
810  do jj = js, je
811  do ii = is, ie
812  distance = intrpnest_haversine( mylat,mylon,inlat(ii,jj),inlon(ii,jj) )
813  distance = distance**weight_order
814  if ( distance <= dist(1) ) then
815  dist(5) = dist(4); igrd(5) = igrd(4); jgrd(5) = jgrd(4)
816  dist(4) = dist(3); igrd(4) = igrd(3); jgrd(4) = jgrd(3)
817  dist(3) = dist(2); igrd(3) = igrd(2); jgrd(3) = jgrd(2)
818  dist(2) = dist(1); igrd(2) = igrd(1); jgrd(2) = jgrd(1)
819  dist(1) = distance; igrd(1) = ii; jgrd(1) = jj
820  elseif ( dist(1) < distance .AND. distance <= dist(2) ) then
821  dist(5) = dist(4); igrd(5) = igrd(4); jgrd(5) = jgrd(4)
822  dist(4) = dist(3); igrd(4) = igrd(3); jgrd(4) = jgrd(3)
823  dist(3) = dist(2); igrd(3) = igrd(2); jgrd(3) = jgrd(2)
824  dist(2) = distance; igrd(2) = ii; jgrd(2) = jj
825  elseif ( dist(2) < distance .AND. distance <= dist(3) ) then
826  dist(5) = dist(4); igrd(5) = igrd(4); jgrd(5) = jgrd(4)
827  dist(4) = dist(3); igrd(4) = igrd(3); jgrd(4) = jgrd(3)
828  dist(3) = distance; igrd(3) = ii; jgrd(3) = jj
829  elseif ( dist(3) < distance .AND. distance <= dist(4) ) then
830  dist(5) = dist(4); igrd(5) = igrd(4); jgrd(5) = jgrd(4)
831  dist(4) = distance; igrd(4) = ii; jgrd(4) = jj
832  elseif ( dist(4) < distance .AND. distance <= dist(5) ) then
833  dist(5) = distance; igrd(5) = ii; jgrd(5) = jj
834  endif
835  enddo
836  enddo
837 
838  if ( abs(dist(1)) < eps ) then
839  hfact(1) = 1.0_rp
840  hfact(2) = 0.0_rp
841  hfact(3) = 0.0_rp
842  hfact(4) = 0.0_rp
843  hfact(5) = 0.0_rp
844  else
845  denom = 1.0_rp / ( (1.0_rp/dist(1)) + (1.0_rp/dist(2)) &
846  + (1.0_rp/dist(3)) + (1.0_rp/dist(4)) + (1.0_rp/dist(5)) )
847  hfact(1) = ( 1.0_rp/dist(1) ) * denom
848  hfact(2) = ( 1.0_rp/dist(2) ) * denom
849  hfact(3) = ( 1.0_rp/dist(3) ) * denom
850  hfact(4) = ( 1.0_rp/dist(4) ) * denom
851  hfact(5) = ( 1.0_rp/dist(5) ) * denom
852  endif
853 
854  return
855  end subroutine intrpnest_search_horiz_5points
856 
857 
858  !-----------------------------------------------------------------------------
859  ! horizontal search of interpolation points for eight-points
860  subroutine intrpnest_search_horiz_8points( &
861  hfact, & ! (out)
862  igrd, & ! (out)
863  jgrd, & ! (out)
864  mylat, & ! (in)
865  mylon, & ! (in)
866  inlat, & ! (in)
867  inlon, & ! (in)
868  is, & ! (in)
869  ie, & ! (in)
870  js, & ! (in)
871  je ) ! (in)
872  use scale_const, only: &
873  eps => const_eps
874  implicit none
875 
876  real(RP), intent(out) :: hfact(:) ! horizontal interp factor
877  integer, intent(out) :: igrd (:) ! grid points of interp target
878  integer, intent(out) :: jgrd (:) ! grid points of interp target
879 
880  real(RP), intent(in) :: mylat ! latitude data of mine
881  real(RP), intent(in) :: mylon ! longitude data of mine
882  real(RP), intent(in) :: inlat(:,:) ! latitude data of you (input)
883  real(RP), intent(in) :: inlon(:,:) ! longitude data of you (input)
884 
885  integer, intent(in) :: is ! start index for x-direction
886  integer, intent(in) :: ie ! end index for x-direction
887  integer, intent(in) :: js ! start index for y-direction
888  integer, intent(in) :: je ! end index for y-direction
889 
890  real(RP) :: distance
891  real(RP) :: denom
892  real(RP) :: dist(8)
893  integer :: ii, jj
894  !---------------------------------------------------------------------------
895 
896  dist(1) = large_number_8
897  dist(2) = large_number_7
898  dist(3) = large_number_6
899  dist(4) = large_number_5
900  dist(5) = large_number_4
901  dist(6) = large_number_3
902  dist(7) = large_number_2
903  dist(8) = large_number_1
904  igrd(:) = -1
905  jgrd(:) = -1
906 
907  do jj = js, je
908  do ii = is, ie
909  distance = intrpnest_haversine( mylat,mylon,inlat(ii,jj),inlon(ii,jj) )
910  distance = distance**weight_order
911  if ( distance <= dist(1) ) then
912  dist(8) = dist(7); igrd(8) = igrd(7); jgrd(8) = jgrd(7)
913  dist(7) = dist(6); igrd(7) = igrd(6); jgrd(7) = jgrd(6)
914  dist(6) = dist(5); igrd(6) = igrd(5); jgrd(6) = jgrd(5)
915  dist(5) = dist(4); igrd(5) = igrd(4); jgrd(5) = jgrd(4)
916  dist(4) = dist(3); igrd(4) = igrd(3); jgrd(4) = jgrd(3)
917  dist(3) = dist(2); igrd(3) = igrd(2); jgrd(3) = jgrd(2)
918  dist(2) = dist(1); igrd(2) = igrd(1); jgrd(2) = jgrd(1)
919  dist(1) = distance; igrd(1) = ii; jgrd(1) = jj
920  elseif ( dist(1) < distance .AND. distance <= dist(2) ) then
921  dist(8) = dist(7); igrd(8) = igrd(7); jgrd(8) = jgrd(7)
922  dist(7) = dist(6); igrd(7) = igrd(6); jgrd(7) = jgrd(6)
923  dist(6) = dist(5); igrd(6) = igrd(5); jgrd(6) = jgrd(5)
924  dist(5) = dist(4); igrd(5) = igrd(4); jgrd(5) = jgrd(4)
925  dist(4) = dist(3); igrd(4) = igrd(3); jgrd(4) = jgrd(3)
926  dist(3) = dist(2); igrd(3) = igrd(2); jgrd(3) = jgrd(2)
927  dist(2) = distance; igrd(2) = ii; jgrd(2) = jj
928  elseif ( dist(2) < distance .AND. distance <= dist(3) ) then
929  dist(8) = dist(7); igrd(8) = igrd(7); jgrd(8) = jgrd(7)
930  dist(7) = dist(6); igrd(7) = igrd(6); jgrd(7) = jgrd(6)
931  dist(6) = dist(5); igrd(6) = igrd(5); jgrd(6) = jgrd(5)
932  dist(5) = dist(4); igrd(5) = igrd(4); jgrd(5) = jgrd(4)
933  dist(4) = dist(3); igrd(4) = igrd(3); jgrd(4) = jgrd(3)
934  dist(3) = distance; igrd(3) = ii; jgrd(3) = jj
935  elseif ( dist(3) < distance .AND. distance <= dist(4) ) then
936  dist(8) = dist(7); igrd(8) = igrd(7); jgrd(8) = jgrd(7)
937  dist(7) = dist(6); igrd(7) = igrd(6); jgrd(7) = jgrd(6)
938  dist(6) = dist(5); igrd(6) = igrd(5); jgrd(6) = jgrd(5)
939  dist(5) = dist(4); igrd(5) = igrd(4); jgrd(5) = jgrd(4)
940  dist(4) = distance; igrd(4) = ii; jgrd(4) = jj
941  elseif ( dist(4) < distance .AND. distance <= dist(5) ) then
942  dist(8) = dist(7); igrd(8) = igrd(7); jgrd(8) = jgrd(7)
943  dist(7) = dist(6); igrd(7) = igrd(6); jgrd(7) = jgrd(6)
944  dist(6) = dist(5); igrd(6) = igrd(5); jgrd(6) = jgrd(5)
945  dist(5) = distance; igrd(5) = ii; jgrd(5) = jj
946  elseif ( dist(5) < distance .AND. distance <= dist(6) ) then
947  dist(8) = dist(7); igrd(8) = igrd(7); jgrd(8) = jgrd(7)
948  dist(7) = dist(6); igrd(7) = igrd(6); jgrd(7) = jgrd(6)
949  dist(6) = distance; igrd(6) = ii; jgrd(6) = jj
950  elseif ( dist(6) < distance .AND. distance <= dist(7) ) then
951  dist(8) = dist(7); igrd(8) = igrd(7); jgrd(8) = jgrd(7)
952  dist(7) = distance; igrd(7) = ii; jgrd(7) = jj
953  elseif ( dist(7) < distance .AND. distance <= dist(8) ) then
954  dist(8) = distance; igrd(8) = ii; jgrd(8) = jj
955  endif
956  enddo
957  enddo
958 
959  if ( abs(dist(1)) < eps ) then
960  hfact(:) = 0.0_rp
961  hfact(1) = 1.0_rp
962  else
963  denom = 1.0_rp / ( (1.0_rp/dist(1)) + (1.0_rp/dist(2)) &
964  + (1.0_rp/dist(3)) + (1.0_rp/dist(4)) &
965  + (1.0_rp/dist(5)) + (1.0_rp/dist(6)) &
966  + (1.0_rp/dist(7)) + (1.0_rp/dist(8)) )
967  hfact(1) = ( 1.0_rp/dist(1) ) * denom
968  hfact(2) = ( 1.0_rp/dist(2) ) * denom
969  hfact(3) = ( 1.0_rp/dist(3) ) * denom
970  hfact(4) = ( 1.0_rp/dist(4) ) * denom
971  hfact(5) = ( 1.0_rp/dist(5) ) * denom
972  hfact(6) = ( 1.0_rp/dist(6) ) * denom
973  hfact(7) = ( 1.0_rp/dist(7) ) * denom
974  hfact(8) = ( 1.0_rp/dist(8) ) * denom
975  endif
976 
977  return
978  end subroutine intrpnest_search_horiz_8points
979 
980 
981  !-----------------------------------------------------------------------------
982  ! horizontal search of interpolation points for twelve-points
983  subroutine intrpnest_search_horiz_12points( &
984  hfact, & ! (out)
985  igrd, & ! (out)
986  jgrd, & ! (out)
987  mylat, & ! (in)
988  mylon, & ! (in)
989  inlat, & ! (in)
990  inlon, & ! (in)
991  is, & ! (in)
992  ie, & ! (in)
993  js, & ! (in)
994  je ) ! (in)
995  use scale_const, only: &
996  eps => const_eps
997  implicit none
998 
999  real(RP), intent(out) :: hfact(:) ! horizontal interp factor
1000  integer, intent(out) :: igrd (:) ! grid points of interp target
1001  integer, intent(out) :: jgrd (:) ! grid points of interp target
1002 
1003  real(RP), intent(in) :: mylat ! latitude data of mine
1004  real(RP), intent(in) :: mylon ! longitude data of mine
1005  real(RP), intent(in) :: inlat(:,:) ! latitude data of you (input)
1006  real(RP), intent(in) :: inlon(:,:) ! longitude data of you (input)
1007 
1008  integer, intent(in) :: is ! start index for x-direction
1009  integer, intent(in) :: ie ! end index for x-direction
1010  integer, intent(in) :: js ! start index for y-direction
1011  integer, intent(in) :: je ! end index for y-direction
1012 
1013  real(RP) :: distance
1014  real(RP) :: denom
1015  real(RP) :: dist(12)
1016  integer :: ii, jj
1017  !---------------------------------------------------------------------------
1018 
1019  dist(1 ) = large_number_12
1020  dist(2 ) = large_number_11
1021  dist(3 ) = large_number_10
1022  dist(4 ) = large_number_9
1023  dist(5 ) = large_number_8
1024  dist(6 ) = large_number_7
1025  dist(7 ) = large_number_6
1026  dist(8 ) = large_number_5
1027  dist(9 ) = large_number_4
1028  dist(10) = large_number_3
1029  dist(11) = large_number_2
1030  dist(12) = large_number_1
1031  igrd(:) = -1
1032  jgrd(:) = -1
1033 
1034  do jj = js, je
1035  do ii = is, ie
1036  distance = intrpnest_haversine( mylat,mylon,inlat(ii,jj),inlon(ii,jj) )
1037  distance = distance**weight_order
1038  if ( distance <= dist(1) ) then
1039  dist(12) = dist(11); igrd(12) = igrd(11); jgrd(12) = jgrd(11)
1040  dist(11) = dist(10); igrd(11) = igrd(10); jgrd(11) = jgrd(10)
1041  dist(10) = dist(9 ); igrd(10) = igrd(9 ); jgrd(10) = jgrd(9 )
1042  dist(9 ) = dist(8 ); igrd(9 ) = igrd(8 ); jgrd(9 ) = jgrd(8 )
1043  dist(8 ) = dist(7 ); igrd(8 ) = igrd(7 ); jgrd(8 ) = jgrd(7 )
1044  dist(7 ) = dist(6 ); igrd(7 ) = igrd(6 ); jgrd(7 ) = jgrd(6 )
1045  dist(6 ) = dist(5 ); igrd(6 ) = igrd(5 ); jgrd(6 ) = jgrd(5 )
1046  dist(5 ) = dist(4 ); igrd(5 ) = igrd(4 ); jgrd(5 ) = jgrd(4 )
1047  dist(4 ) = dist(3 ); igrd(4 ) = igrd(3 ); jgrd(4 ) = jgrd(3 )
1048  dist(3 ) = dist(2 ); igrd(3 ) = igrd(2 ); jgrd(3 ) = jgrd(2 )
1049  dist(2 ) = dist(1 ); igrd(2 ) = igrd(1 ); jgrd(2 ) = jgrd(1 )
1050  dist(1 ) = distance; igrd(1 ) = ii; jgrd(1 ) = jj
1051  elseif ( dist(1) < distance .AND. distance <= dist(2) ) then
1052  dist(12) = dist(11); igrd(12) = igrd(11); jgrd(12) = jgrd(11)
1053  dist(11) = dist(10); igrd(11) = igrd(10); jgrd(11) = jgrd(10)
1054  dist(10) = dist(9 ); igrd(10) = igrd(9 ); jgrd(10) = jgrd(9 )
1055  dist(9 ) = dist(8 ); igrd(9 ) = igrd(8 ); jgrd(9 ) = jgrd(8 )
1056  dist(8 ) = dist(7 ); igrd(8 ) = igrd(7 ); jgrd(8 ) = jgrd(7 )
1057  dist(7 ) = dist(6 ); igrd(7 ) = igrd(6 ); jgrd(7 ) = jgrd(6 )
1058  dist(6 ) = dist(5 ); igrd(6 ) = igrd(5 ); jgrd(6 ) = jgrd(5 )
1059  dist(5 ) = dist(4 ); igrd(5 ) = igrd(4 ); jgrd(5 ) = jgrd(4 )
1060  dist(4 ) = dist(3 ); igrd(4 ) = igrd(3 ); jgrd(4 ) = jgrd(3 )
1061  dist(3 ) = dist(2 ); igrd(3 ) = igrd(2 ); jgrd(3 ) = jgrd(2 )
1062  dist(2 ) = distance; igrd(2 ) = ii; jgrd(2 ) = jj
1063  elseif ( dist(2) < distance .AND. distance <= dist(3) ) then
1064  dist(12) = dist(11); igrd(12) = igrd(11); jgrd(12) = jgrd(11)
1065  dist(11) = dist(10); igrd(11) = igrd(10); jgrd(11) = jgrd(10)
1066  dist(10) = dist(9 ); igrd(10) = igrd(9 ); jgrd(10) = jgrd(9 )
1067  dist(9 ) = dist(8 ); igrd(9 ) = igrd(8 ); jgrd(9 ) = jgrd(8 )
1068  dist(8 ) = dist(7 ); igrd(8 ) = igrd(7 ); jgrd(8 ) = jgrd(7 )
1069  dist(7 ) = dist(6 ); igrd(7 ) = igrd(6 ); jgrd(7 ) = jgrd(6 )
1070  dist(6 ) = dist(5 ); igrd(6 ) = igrd(5 ); jgrd(6 ) = jgrd(5 )
1071  dist(5 ) = dist(4 ); igrd(5 ) = igrd(4 ); jgrd(5 ) = jgrd(4 )
1072  dist(4 ) = dist(3 ); igrd(4 ) = igrd(3 ); jgrd(4 ) = jgrd(3 )
1073  dist(3 ) = distance; igrd(3 ) = ii; jgrd(3 ) = jj
1074  elseif ( dist(3) < distance .AND. distance <= dist(4) ) then
1075  dist(12) = dist(11); igrd(12) = igrd(11); jgrd(12) = jgrd(11)
1076  dist(11) = dist(10); igrd(11) = igrd(10); jgrd(11) = jgrd(10)
1077  dist(10) = dist(9 ); igrd(10) = igrd(9 ); jgrd(10) = jgrd(9 )
1078  dist(9 ) = dist(8 ); igrd(9 ) = igrd(8 ); jgrd(9 ) = jgrd(8 )
1079  dist(8 ) = dist(7 ); igrd(8 ) = igrd(7 ); jgrd(8 ) = jgrd(7 )
1080  dist(7 ) = dist(6 ); igrd(7 ) = igrd(6 ); jgrd(7 ) = jgrd(6 )
1081  dist(6 ) = dist(5 ); igrd(6 ) = igrd(5 ); jgrd(6 ) = jgrd(5 )
1082  dist(5 ) = dist(4 ); igrd(5 ) = igrd(4 ); jgrd(5 ) = jgrd(4 )
1083  dist(4 ) = distance; igrd(4 ) = ii; jgrd(4 ) = jj
1084  elseif ( dist(4) < distance .AND. distance <= dist(5) ) then
1085  dist(12) = dist(11); igrd(12) = igrd(11); jgrd(12) = jgrd(11)
1086  dist(11) = dist(10); igrd(11) = igrd(10); jgrd(11) = jgrd(10)
1087  dist(10) = dist(9 ); igrd(10) = igrd(9 ); jgrd(10) = jgrd(9 )
1088  dist(9 ) = dist(8 ); igrd(9 ) = igrd(8 ); jgrd(9 ) = jgrd(8 )
1089  dist(8 ) = dist(7 ); igrd(8 ) = igrd(7 ); jgrd(8 ) = jgrd(7 )
1090  dist(7 ) = dist(6 ); igrd(7 ) = igrd(6 ); jgrd(7 ) = jgrd(6 )
1091  dist(6 ) = dist(5 ); igrd(6 ) = igrd(5 ); jgrd(6 ) = jgrd(5 )
1092  dist(5 ) = distance; igrd(5 ) = ii; jgrd(5 ) = jj
1093  elseif ( dist(5) < distance .AND. distance <= dist(6) ) then
1094  dist(12) = dist(11); igrd(12) = igrd(11); jgrd(12) = jgrd(11)
1095  dist(11) = dist(10); igrd(11) = igrd(10); jgrd(11) = jgrd(10)
1096  dist(10) = dist(9 ); igrd(10) = igrd(9 ); jgrd(10) = jgrd(9 )
1097  dist(9 ) = dist(8 ); igrd(9 ) = igrd(8 ); jgrd(9 ) = jgrd(8 )
1098  dist(8 ) = dist(7 ); igrd(8 ) = igrd(7 ); jgrd(8 ) = jgrd(7 )
1099  dist(7 ) = dist(6 ); igrd(7 ) = igrd(6 ); jgrd(7 ) = jgrd(6 )
1100  dist(6 ) = distance; igrd(6 ) = ii; jgrd(6 ) = jj
1101  elseif ( dist(6) < distance .AND. distance <= dist(7) ) then
1102  dist(12) = dist(11); igrd(12) = igrd(11); jgrd(12) = jgrd(11)
1103  dist(11) = dist(10); igrd(11) = igrd(10); jgrd(11) = jgrd(10)
1104  dist(10) = dist(9 ); igrd(10) = igrd(9 ); jgrd(10) = jgrd(9 )
1105  dist(9 ) = dist(8 ); igrd(9 ) = igrd(8 ); jgrd(9 ) = jgrd(8 )
1106  dist(8 ) = dist(7 ); igrd(8 ) = igrd(7 ); jgrd(8 ) = jgrd(7 )
1107  dist(7 ) = distance; igrd(7 ) = ii; jgrd(7 ) = jj
1108  elseif ( dist(7) < distance .AND. distance <= dist(8) ) then
1109  dist(12) = dist(11); igrd(12) = igrd(11); jgrd(12) = jgrd(11)
1110  dist(11) = dist(10); igrd(11) = igrd(10); jgrd(11) = jgrd(10)
1111  dist(10) = dist(9 ); igrd(10) = igrd(9 ); jgrd(10) = jgrd(9 )
1112  dist(9 ) = dist(8 ); igrd(9 ) = igrd(8 ); jgrd(9 ) = jgrd(8 )
1113  dist(8 ) = distance; igrd(8 ) = ii; jgrd(8 ) = jj
1114  elseif ( dist(8) < distance .AND. distance <= dist(9) ) then
1115  dist(12) = dist(11); igrd(12) = igrd(11); jgrd(12) = jgrd(11)
1116  dist(11) = dist(10); igrd(11) = igrd(10); jgrd(11) = jgrd(10)
1117  dist(10) = dist(9 ); igrd(10) = igrd(9 ); jgrd(10) = jgrd(9 )
1118  dist(9 ) = distance; igrd(9 ) = ii; jgrd(9 ) = jj
1119  elseif ( dist(9) < distance .AND. distance <= dist(10) ) then
1120  dist(12) = dist(11); igrd(12) = igrd(11); jgrd(12) = jgrd(11)
1121  dist(11) = dist(10); igrd(11) = igrd(10); jgrd(11) = jgrd(10)
1122  dist(10) = distance; igrd(10) = ii; jgrd(10) = jj
1123  elseif ( dist(10) < distance .AND. distance <= dist(11) ) then
1124  dist(12) = dist(11); igrd(12) = igrd(11); jgrd(12) = jgrd(11)
1125  dist(11) = distance; igrd(11) = ii; jgrd(11) = jj
1126  elseif ( dist(11) < distance .AND. distance <= dist(12) ) then
1127  dist(12) = distance; igrd(12) = ii; jgrd(12) = jj
1128  endif
1129  enddo
1130  enddo
1131 
1132  if ( abs(dist(1)) < eps ) then
1133  hfact(:) = 0.0_rp
1134  hfact(1) = 1.0_rp
1135  else
1136  denom = 1.0_rp / ( (1.0_rp/dist(1 )) + (1.0_rp/dist(2 )) &
1137  + (1.0_rp/dist(3 )) + (1.0_rp/dist(4 )) &
1138  + (1.0_rp/dist(5 )) + (1.0_rp/dist(6 )) &
1139  + (1.0_rp/dist(7 )) + (1.0_rp/dist(8 )) &
1140  + (1.0_rp/dist(9 )) + (1.0_rp/dist(10)) &
1141  + (1.0_rp/dist(11)) + (1.0_rp/dist(12)) )
1142  hfact(1 ) = ( 1.0_rp/dist(1 ) ) * denom
1143  hfact(2 ) = ( 1.0_rp/dist(2 ) ) * denom
1144  hfact(3 ) = ( 1.0_rp/dist(3 ) ) * denom
1145  hfact(4 ) = ( 1.0_rp/dist(4 ) ) * denom
1146  hfact(5 ) = ( 1.0_rp/dist(5 ) ) * denom
1147  hfact(6 ) = ( 1.0_rp/dist(6 ) ) * denom
1148  hfact(7 ) = ( 1.0_rp/dist(7 ) ) * denom
1149  hfact(8 ) = ( 1.0_rp/dist(8 ) ) * denom
1150  hfact(9 ) = ( 1.0_rp/dist(9 ) ) * denom
1151  hfact(10) = ( 1.0_rp/dist(10) ) * denom
1152  hfact(11) = ( 1.0_rp/dist(11) ) * denom
1153  hfact(12) = ( 1.0_rp/dist(12) ) * denom
1154  endif
1155 
1156  return
1157  end subroutine intrpnest_search_horiz_12points
1158 
1159 
1160  !-----------------------------------------------------------------------------
1161  ! vertical search of interpolation points for two-points (online)
1162  subroutine intrpnest_search_vert_online( &
1163  vfact, & ! (inout)
1164  kgrd, & ! (inout)
1165  ncopy, & ! (out)
1166  igrd, & ! (in)
1167  jgrd, & ! (in)
1168  myhgt, & ! (in)
1169  inhgt, & ! (in)
1170  iloc, & ! (in)
1171  jloc, & ! (in)
1172  ks, & ! (in)
1173  ke, & ! (in)
1174  inka, & ! (in)
1175  lndgrd ) ! (in)
1176  use scale_process, only: &
1177  prc_mpistop
1178  implicit none
1179 
1180  real(RP), intent(inout) :: vfact(:,:,:,:,:) ! vertical interp factor
1181  integer, intent(inout) :: kgrd (:,:,:,:,:) ! grid points of interp target
1182  integer, intent(out) :: ncopy(:) ! number of daughter's layers below inKS
1183 
1184  integer, intent(in) :: igrd(:) ! grid points of interp target
1185  integer, intent(in) :: jgrd(:) ! grid points of interp target
1186  real(RP), intent(in) :: myhgt(:) ! height data of mine
1187  real(RP), intent(in) :: inhgt(:,:,:) ! height data of you (input)
1188  integer, intent(in) :: iloc ! locator index for x-direction
1189  integer, intent(in) :: jloc ! locator index for y-direction
1190  integer, intent(in) :: ks ! start index for z-direction
1191  integer, intent(in) :: ke ! end index for z-direction
1192  integer, intent(in) :: inka ! grid number of you (input)
1193  logical, intent(in) :: lndgrd ! flag of land grid
1194 
1195  real(RP) :: dist(2)
1196  logical :: dflag ! flag: data found or not
1197  integer :: ii, jj, idx
1198  integer :: k, kk
1199  integer :: inks, inke
1200  logical :: copy
1201  !---------------------------------------------------------------------------
1202 
1203  if ( lndgrd ) then
1204  write(*,*) 'xxx internal error [interporation: nest/interp]'
1205  write(*,*) ' land grid is not araviable in online'
1206  call prc_mpistop
1207  endif
1208 
1209  inks = 1 + khalo
1210  inke = inka - khalo
1211 
1212  do idx = 1, itp_nh
1213  ii = igrd(idx)
1214  jj = jgrd(idx)
1215 
1216  do k = ks, ke
1217  dist(1) = large_number_2
1218  dist(2) = large_number_1
1219  copy = .false.
1220 
1221  kgrd(k,iloc,jloc,idx,:) = -1
1222 
1223  dflag = .false.
1224 
1225  if( myhgt(k) < inhgt(inks,ii,jj) ) then
1226  copy = .true.
1227  ncopy(idx) = ncopy(idx) + 1
1228 
1229  kgrd(k,iloc,jloc,idx,:) = inks
1230 
1231  vfact(k,iloc,jloc,idx,1) = 1.0_rp
1232  vfact(k,iloc,jloc,idx,2) = 0.0_rp
1233 
1234  dflag = .true.
1235  else
1236  do kk = inks, inke
1237  dist(1) = myhgt(k) - inhgt(kk ,ii,jj)
1238  dist(2) = myhgt(k) - inhgt(kk+1,ii,jj)
1239 
1240  if( dist(1) >= 0.0_rp .AND. dist(2) < 0.0_rp ) then
1241  kgrd(k,iloc,jloc,idx,1) = kk
1242  kgrd(k,iloc,jloc,idx,2) = kk+1
1243 
1244  vfact(k,iloc,jloc,idx,1) = abs(dist(2)) / ( abs(dist(1)) + abs(dist(2)) )
1245  vfact(k,iloc,jloc,idx,2) = abs(dist(1)) / ( abs(dist(1)) + abs(dist(2)) )
1246 
1247  dflag = .true.
1248 
1249  exit ! loop end
1250  endif
1251  enddo
1252  endif
1253 
1254  if( .NOT. dflag ) then
1255  write(*,*) 'xxx internal error [INTRPNEST_search_vert_online]'
1256  write(*,*) 'xxx data for interpolation was not found.'
1257  write(*,*) 'xxx iloc=',iloc,' jloc=',jloc,' k=',k,' idx=',idx
1258  call prc_mpistop
1259  endif
1260  enddo
1261 
1262  enddo
1263 
1264  return
1265  end subroutine intrpnest_search_vert_online
1266 
1267  !-----------------------------------------------------------------------------
1268  ! vertical search of interpolation points for two-points (offline)
1269  subroutine intrpnest_search_vert_offline( &
1270  vfact, & ! (inout)
1271  kgrd, & ! (inout)
1272  ncopy, & ! (out)
1273  igrd, & ! (in)
1274  jgrd, & ! (in)
1275  myhgt, & ! (in)
1276  inhgt, & ! (in)
1277  iloc, & ! (in)
1278  jloc, & ! (in)
1279  ks, & ! (in)
1280  ke, & ! (in)
1281  inka, & ! (in)
1282  lndgrd ) ! (in)
1283  use scale_const, only: &
1284  eps => const_eps
1285  use scale_process, only: &
1286  prc_mpistop
1287  implicit none
1288 
1289  real(RP), intent(inout) :: vfact(:,:,:,:,:) ! vertical interp factor
1290  integer, intent(inout) :: kgrd (:,:,:,:,:) ! grid points of interp target
1291  integer, intent(out) :: ncopy(:) ! number of daughter's layers below inKS
1292 
1293  integer, intent(in) :: igrd(:) ! grid points of interp target
1294  integer, intent(in) :: jgrd(:) ! grid points of interp target
1295  real(RP), intent(in) :: myhgt(:) ! height data of mine
1296  real(RP), intent(in) :: inhgt(:,:,:) ! height data of you (input)
1297  integer, intent(in) :: iloc ! locator index for x-direction
1298  integer, intent(in) :: jloc ! locator index for y-direction
1299  integer, intent(in) :: ks ! start index for z-direction
1300  integer, intent(in) :: ke ! end index for z-direction
1301  integer, intent(in) :: inka ! grid number of you (input)
1302  logical, intent(in) :: lndgrd ! flag of land grid
1303 
1304  real(RP) :: denom
1305  real(RP) :: dist(2)
1306  logical :: dflag ! flag: data found or not
1307  integer :: ii, jj
1308  integer :: k, kk, kks, kke
1309  integer :: idx
1310  !---------------------------------------------------------------------------
1311 
1312  if ( lndgrd ) then
1313  kks = 1; kke = lkmax
1314  else
1315  kks = ks; kke = ke
1316  endif
1317 
1318  ncopy = 0 ! dummy
1319  do idx = 1, itp_nh
1320  ii = igrd(idx)
1321  jj = jgrd(idx)
1322 
1323  do k = kks, kke
1324  dist(1) = large_number_2
1325  dist(2) = large_number_1
1326  kgrd(k,iloc,jloc,idx,:) = -1
1327  dflag = .false.
1328 
1329  if( myhgt(k) < inhgt(1,ii,jj) )then ! copy
1330  kgrd(k,iloc,jloc,idx,:) = 1
1331  vfact(k,iloc,jloc,idx,1) = 1.0_rp
1332  vfact(k,iloc,jloc,idx,2) = 0.0_rp
1333  dflag = .true.
1334  else if( abs(inhgt(inka,ii,jj)-myhgt(k))<eps )then
1335  kgrd(k,iloc,jloc,idx,:) = inka
1336  vfact(k,iloc,jloc,idx,1) = 1.0_rp
1337  vfact(k,iloc,jloc,idx,2) = 0.0_rp
1338  dflag = .true.
1339  else if( inhgt(inka,ii,jj) < myhgt(k) )then
1340  if ( lndgrd ) then ! copy
1341  kgrd(k,iloc,jloc,idx,:) = inka
1342  vfact(k,iloc,jloc,idx,1) = 1.0_rp
1343  vfact(k,iloc,jloc,idx,2) = 0.0_rp
1344  dflag = .true.
1345  else
1346  write(*,*) 'xxx internal error [INTRPNEST_search_vert_offline]'
1347  write(*,*) 'xxx data level is beyond parent data'
1348  write(*,*) 'in',ii,jj,inka,inhgt(inka,ii,jj),'my',iloc,jloc,k,myhgt(k)
1349  call prc_mpistop
1350  end if
1351  else
1352 
1353  do kk = 1, inka-1
1354  if( (inhgt(kk,ii,jj)<=myhgt(k)).AND.(myhgt(k)<inhgt(kk+1,ii,jj)) )then
1355  kgrd(k,iloc,jloc,idx,1) = kk
1356  kgrd(k,iloc,jloc,idx,2) = kk+1
1357  dist(1) = abs( myhgt(k) - inhgt(kk,ii,jj) )
1358  dist(2) = abs( myhgt(k) - inhgt(kk+1,ii,jj) )
1359  dflag = .true.
1360  if ( abs(dist(1))<eps )then
1361  vfact(k,iloc,jloc,idx,1) = 1.0_rp
1362  vfact(k,iloc,jloc,idx,2) = 0.0_rp
1363  else
1364  denom = 1.0_rp / ( (1.0_rp/dist(1)) + (1.0_rp/dist(2)) )
1365  vfact(k,iloc,jloc,idx,1) = ( 1.0_rp/dist(1) ) * denom
1366  vfact(k,iloc,jloc,idx,2) = ( 1.0_rp/dist(2) ) * denom
1367  endif
1368  exit
1369  endif
1370  enddo
1371  endif
1372 
1373  if( .NOT. dflag )then
1374  write(*,*) 'xxx internal error [INTRPNEST_search_vert_offline]'
1375  write(*,*) 'xxx data for interpolation was not found.'
1376  write(*,*) 'xxx iloc=',iloc,' jloc=',jloc,' k=',k,' idx=',idx
1377  call prc_mpistop
1378  endif
1379  enddo
1380 
1381  enddo
1382 
1383  return
1384  end subroutine intrpnest_search_vert_offline
1385 
1386  !-----------------------------------------------------------------------------
1387  ! interpolation using one-points for 2D data (nearest-neighbor)
1388  subroutine intrpnest_interp_2d_1points( &
1389  intp, & ! (out)
1390  ref, & ! (in)
1391  hfact, & ! (in)
1392  igrd, & ! (in)
1393  jgrd, & ! (in)
1394  ia, & ! (in)
1395  ja ) ! (in)
1396  implicit none
1397 
1398  real(RP), intent(out) :: intp(:,:) ! interpolated data
1399 
1400  real(RP), intent(in) :: ref (:,:) ! reference data
1401  real(RP), intent(in) :: hfact(:,:,:) ! horizontal interp factor
1402  integer, intent(in) :: igrd (:,:,:) ! grid points of interp target
1403  integer, intent(in) :: jgrd (:,:,:) ! grid points of interp target
1404  integer, intent(in) :: ia ! grid number of mine
1405  integer, intent(in) :: ja ! grid number of mine
1406 
1407  integer :: i, j
1408  !---------------------------------------------------------------------------
1409 
1410 !OCL PREFETCH
1411  do j = 1, ja
1412  do i = 1, ia
1413  intp(i,j) = ref(igrd(i,j,1),jgrd(i,j,1)) * hfact(i,j,1)
1414  end do
1415  end do
1416 
1417  return
1418  end subroutine intrpnest_interp_2d_1points
1419 
1420 
1421  !-----------------------------------------------------------------------------
1422  ! interpolation using one-points for 3D data (nearest-neighbor)
1423  subroutine intrpnest_interp_3d_1points( &
1424  intp, & ! (out)
1425  ref, & ! (in)
1426  hfact, & ! (in)
1427  vfact, & ! (in)
1428  kgrd, & ! (in)
1429  igrd, & ! (in)
1430  jgrd, & ! (in)
1431  ia, & ! (in)
1432  ja, & ! (in)
1433  ks, & ! (in)
1434  ke, & ! (in)
1435  logwegt ) ! (in)
1436  implicit none
1437 
1438  real(RP), intent(out) :: intp(:,:,:) ! interpolated data
1439 
1440  real(RP), intent(in) :: ref (:,:,:) ! reference data
1441  real(RP), intent(in) :: hfact(:,:,:) ! horizontal interp factor
1442  real(RP), intent(in) :: vfact(:,:,:,:,:) ! vertical interp factor
1443  integer, intent(in) :: kgrd (:,:,:,:,:) ! grid points of interp target
1444  integer, intent(in) :: igrd (:,:,:) ! grid points of interp target
1445  integer, intent(in) :: jgrd (:,:,:) ! grid points of interp target
1446  integer, intent(in) :: ia ! grid number of mine
1447  integer, intent(in) :: ja ! grid number of mine
1448  integer, intent(in) :: ks ! start grid number of mine
1449  integer, intent(in) :: ke ! end grid number of mine
1450 
1451  logical, intent(in), optional :: logwegt
1452 
1453  integer :: i, j, k
1454  logical :: logarithmic
1455  !---------------------------------------------------------------------------
1456 
1457  logarithmic = .false.
1458  if ( present(logwegt) ) then
1459  if ( logwegt ) then
1460  logarithmic = .true.
1461  endif
1462  endif
1463 
1464  ! linear interpolation
1465 !OCL PREFETCH
1466  do j = 1, ja
1467  do i = 1, ia
1468  do k = ks, ke
1469  intp(k,i,j) = ref(kgrd(k,i,j,1,1),igrd(i,j,1),jgrd(i,j,1)) &
1470  * hfact(i,j,1) * vfact(k,i,j,1,1) &
1471  + ref(kgrd(k,i,j,1,2),igrd(i,j,1),jgrd(i,j,1)) &
1472  * hfact(i,j,1) * vfact(k,i,j,1,2)
1473  end do
1474  end do
1475  end do
1476 
1477  ! logarithmic weighting (for pres, dens)
1478  if ( logarithmic ) then
1479  do j = 1, ja
1480  do i = 1, ia
1481  do k = ks, ke
1482  intp(k,i,j) = exp( intp(k,i,j) )
1483  end do
1484  end do
1485  end do
1486  endif
1487 
1488  return
1489  end subroutine intrpnest_interp_3d_1points
1490 
1491 
1492  !-----------------------------------------------------------------------------
1493  ! interpolation using three-points for 2D data
1494  subroutine intrpnest_interp_2d_3points( &
1495  intp, & ! (out)
1496  ref, & ! (in)
1497  hfact, & ! (in)
1498  igrd, & ! (in)
1499  jgrd, & ! (in)
1500  ia, & ! (in)
1501  ja ) ! (in)
1502  implicit none
1503 
1504  real(RP), intent(out) :: intp(:,:) ! interpolated data
1505 
1506  real(RP), intent(in) :: ref (:,:) ! reference data
1507  real(RP), intent(in) :: hfact(:,:,:) ! horizontal interp factor
1508  integer, intent(in) :: igrd (:,:,:) ! grid points of interp target
1509  integer, intent(in) :: jgrd (:,:,:) ! grid points of interp target
1510  integer, intent(in) :: ia ! grid number of mine
1511  integer, intent(in) :: ja ! grid number of mine
1512 
1513  integer :: i, j
1514  !---------------------------------------------------------------------------
1515 
1516 !OCL PREFETCH
1517  do j = 1, ja
1518  do i = 1, ia
1519  intp(i,j) = ref(igrd(i,j,1),jgrd(i,j,1)) * hfact(i,j,1) &
1520  + ref(igrd(i,j,2),jgrd(i,j,2)) * hfact(i,j,2) &
1521  + ref(igrd(i,j,3),jgrd(i,j,3)) * hfact(i,j,3)
1522  end do
1523  end do
1524 
1525  return
1526  end subroutine intrpnest_interp_2d_3points
1527 
1528 
1529  !-----------------------------------------------------------------------------
1530  ! interpolation using three-points for 3D data
1531  subroutine intrpnest_interp_3d_3points( &
1532  intp, & ! (out)
1533  ref, & ! (in)
1534  hfact, & ! (in)
1535  vfact, & ! (in)
1536  kgrd, & ! (in)
1537  igrd, & ! (in)
1538  jgrd, & ! (in)
1539  ia, & ! (in)
1540  ja, & ! (in)
1541  ks, & ! (in)
1542  ke, & ! (in)
1543  logwegt ) ! (in)
1544  implicit none
1545 
1546  real(RP), intent(out) :: intp(:,:,:) ! interpolated data
1547 
1548  real(RP), intent(in) :: ref (:,:,:) ! reference data
1549  real(RP), intent(in) :: hfact(:,:,:) ! horizontal interp factor
1550  real(RP), intent(in) :: vfact(:,:,:,:,:) ! vertical interp factor
1551  integer, intent(in) :: kgrd (:,:,:,:,:) ! grid points of interp target
1552  integer, intent(in) :: igrd (:,:,:) ! grid points of interp target
1553  integer, intent(in) :: jgrd (:,:,:) ! grid points of interp target
1554  integer, intent(in) :: ia ! grid number of mine
1555  integer, intent(in) :: ja ! grid number of mine
1556  integer, intent(in) :: ks ! start grid number of mine
1557  integer, intent(in) :: ke ! end grid number of mine
1558 
1559  logical, intent(in), optional :: logwegt
1560 
1561  integer :: i, j, k
1562  logical :: logarithmic
1563  !---------------------------------------------------------------------------
1564 
1565  logarithmic = .false.
1566  if ( present(logwegt) ) then
1567  if ( logwegt ) then
1568  logarithmic = .true.
1569  endif
1570  endif
1571 
1572  ! linear interpolation
1573 !OCL PREFETCH
1574  do j = 1, ja
1575  do i = 1, ia
1576  do k = ks, ke
1577  intp(k,i,j) = ref(kgrd(k,i,j,1,1),igrd(i,j,1),jgrd(i,j,1)) &
1578  * hfact(i,j,1) * vfact(k,i,j,1,1) &
1579  + ref(kgrd(k,i,j,2,1),igrd(i,j,2),jgrd(i,j,2)) &
1580  * hfact(i,j,2) * vfact(k,i,j,2,1) &
1581  + ref(kgrd(k,i,j,3,1),igrd(i,j,3),jgrd(i,j,3)) &
1582  * hfact(i,j,3) * vfact(k,i,j,3,1) &
1583  + ref(kgrd(k,i,j,1,2),igrd(i,j,1),jgrd(i,j,1)) &
1584  * hfact(i,j,1) * vfact(k,i,j,1,2) &
1585  + ref(kgrd(k,i,j,2,2),igrd(i,j,2),jgrd(i,j,2)) &
1586  * hfact(i,j,2) * vfact(k,i,j,2,2) &
1587  + ref(kgrd(k,i,j,3,2),igrd(i,j,3),jgrd(i,j,3)) &
1588  * hfact(i,j,3) * vfact(k,i,j,3,2)
1589  end do
1590  end do
1591  end do
1592 
1593  ! logarithmic weighting (for pres, dens)
1594  if ( logarithmic ) then
1595  do j = 1, ja
1596  do i = 1, ia
1597  do k = ks, ke
1598  intp(k,i,j) = exp( intp(k,i,j) )
1599  end do
1600  end do
1601  end do
1602  endif
1603 
1604  return
1605  end subroutine intrpnest_interp_3d_3points
1606 
1607 
1608  !-----------------------------------------------------------------------------
1609  ! interpolation using four-points for 2D data
1610  subroutine intrpnest_interp_2d_4points( &
1611  intp, & ! (out)
1612  ref, & ! (in)
1613  hfact, & ! (in)
1614  igrd, & ! (in)
1615  jgrd, & ! (in)
1616  ia, & ! (in)
1617  ja ) ! (in)
1618  implicit none
1619 
1620  real(RP), intent(out) :: intp(:,:) ! interpolated data
1621 
1622  real(RP), intent(in) :: ref (:,:) ! reference data
1623  real(RP), intent(in) :: hfact(:,:,:) ! horizontal interp factor
1624  integer, intent(in) :: igrd (:,:,:) ! grid points of interp target
1625  integer, intent(in) :: jgrd (:,:,:) ! grid points of interp target
1626  integer, intent(in) :: ia ! grid number of mine
1627  integer, intent(in) :: ja ! grid number of mine
1628 
1629  integer :: i, j
1630  !---------------------------------------------------------------------------
1631 
1632 !OCL PREFETCH
1633  do j = 1, ja
1634  do i = 1, ia
1635  intp(i,j) = ref(igrd(i,j,1),jgrd(i,j,1)) * hfact(i,j,1) &
1636  + ref(igrd(i,j,2),jgrd(i,j,2)) * hfact(i,j,2) &
1637  + ref(igrd(i,j,3),jgrd(i,j,3)) * hfact(i,j,3) &
1638  + ref(igrd(i,j,4),jgrd(i,j,4)) * hfact(i,j,4)
1639  end do
1640  end do
1641 
1642  return
1643  end subroutine intrpnest_interp_2d_4points
1644 
1645 
1646  !-----------------------------------------------------------------------------
1647  ! interpolation using four-points for 3D data
1648  subroutine intrpnest_interp_3d_4points( &
1649  intp, & ! (out)
1650  ref, & ! (in)
1651  hfact, & ! (in)
1652  vfact, & ! (in)
1653  kgrd, & ! (in)
1654  igrd, & ! (in)
1655  jgrd, & ! (in)
1656  ia, & ! (in)
1657  ja, & ! (in)
1658  ks, & ! (in)
1659  ke, & ! (in)
1660  logwegt ) ! (in)
1661  implicit none
1662 
1663  real(RP), intent(out) :: intp(:,:,:) ! interpolated data
1664 
1665  real(RP), intent(in) :: ref (:,:,:) ! reference data
1666  real(RP), intent(in) :: hfact(:,:,:) ! horizontal interp factor
1667  real(RP), intent(in) :: vfact(:,:,:,:,:) ! vertical interp factor
1668  integer, intent(in) :: kgrd (:,:,:,:,:) ! grid points of interp target
1669  integer, intent(in) :: igrd (:,:,:) ! grid points of interp target
1670  integer, intent(in) :: jgrd (:,:,:) ! grid points of interp target
1671  integer, intent(in) :: ia ! grid number of mine
1672  integer, intent(in) :: ja ! grid number of mine
1673  integer, intent(in) :: ks ! start grid number of mine
1674  integer, intent(in) :: ke ! end grid number of mine
1675 
1676  logical, intent(in), optional :: logwegt
1677 
1678  integer :: i, j, k
1679  logical :: logarithmic
1680  !---------------------------------------------------------------------------
1681 
1682  logarithmic = .false.
1683  if ( present(logwegt) ) then
1684  if ( logwegt ) then
1685  logarithmic = .true.
1686  endif
1687  endif
1688 
1689  ! linear interpolation
1690 !OCL PREFETCH
1691  do j = 1, ja
1692  do i = 1, ia
1693  do k = ks, ke
1694  intp(k,i,j) = ref(kgrd(k,i,j,1,1),igrd(i,j,1),jgrd(i,j,1)) &
1695  * hfact(i,j,1) * vfact(k,i,j,1,1) &
1696  + ref(kgrd(k,i,j,2,1),igrd(i,j,2),jgrd(i,j,2)) &
1697  * hfact(i,j,2) * vfact(k,i,j,2,1) &
1698  + ref(kgrd(k,i,j,3,1),igrd(i,j,3),jgrd(i,j,3)) &
1699  * hfact(i,j,3) * vfact(k,i,j,3,1) &
1700  + ref(kgrd(k,i,j,4,1),igrd(i,j,4),jgrd(i,j,4)) &
1701  * hfact(i,j,4) * vfact(k,i,j,4,1) &
1702  + ref(kgrd(k,i,j,1,2),igrd(i,j,1),jgrd(i,j,1)) &
1703  * hfact(i,j,1) * vfact(k,i,j,1,2) &
1704  + ref(kgrd(k,i,j,2,2),igrd(i,j,2),jgrd(i,j,2)) &
1705  * hfact(i,j,2) * vfact(k,i,j,2,2) &
1706  + ref(kgrd(k,i,j,3,2),igrd(i,j,3),jgrd(i,j,3)) &
1707  * hfact(i,j,3) * vfact(k,i,j,3,2) &
1708  + ref(kgrd(k,i,j,4,2),igrd(i,j,4),jgrd(i,j,4)) &
1709  * hfact(i,j,4) * vfact(k,i,j,4,2)
1710  end do
1711  end do
1712  end do
1713 
1714  ! logarithmic weighting (for pres, dens)
1715  if ( logarithmic ) then
1716  do j = 1, ja
1717  do i = 1, ia
1718  do k = ks, ke
1719  intp(k,i,j) = exp( intp(k,i,j) )
1720  end do
1721  end do
1722  end do
1723  endif
1724 
1725  return
1726  end subroutine intrpnest_interp_3d_4points
1727 
1728 
1729  !-----------------------------------------------------------------------------
1730  ! interpolation using five-points for 2D data
1731  subroutine intrpnest_interp_2d_5points( &
1732  intp, & ! (out)
1733  ref, & ! (in)
1734  hfact, & ! (in)
1735  igrd, & ! (in)
1736  jgrd, & ! (in)
1737  ia, & ! (in)
1738  ja ) ! (in)
1739  implicit none
1740 
1741  real(RP), intent(out) :: intp(:,:) ! interpolated data
1742 
1743  real(RP), intent(in) :: ref (:,:) ! reference data
1744  real(RP), intent(in) :: hfact(:,:,:) ! horizontal interp factor
1745  integer, intent(in) :: igrd (:,:,:) ! grid points of interp target
1746  integer, intent(in) :: jgrd (:,:,:) ! grid points of interp target
1747  integer, intent(in) :: ia ! grid number of mine
1748  integer, intent(in) :: ja ! grid number of mine
1749 
1750  integer :: i, j
1751  !---------------------------------------------------------------------------
1752 
1753 !OCL PREFETCH
1754  do j = 1, ja
1755  do i = 1, ia
1756  intp(i,j) = ref(igrd(i,j,1),jgrd(i,j,1)) * hfact(i,j,1) &
1757  + ref(igrd(i,j,2),jgrd(i,j,2)) * hfact(i,j,2) &
1758  + ref(igrd(i,j,3),jgrd(i,j,3)) * hfact(i,j,3) &
1759  + ref(igrd(i,j,4),jgrd(i,j,4)) * hfact(i,j,4) &
1760  + ref(igrd(i,j,5),jgrd(i,j,5)) * hfact(i,j,5)
1761  end do
1762  end do
1763 
1764  return
1765  end subroutine intrpnest_interp_2d_5points
1766 
1767 
1768  !-----------------------------------------------------------------------------
1769  ! interpolation using five-points for 3D data
1770  subroutine intrpnest_interp_3d_5points( &
1771  intp, & ! (out)
1772  ref, & ! (in)
1773  hfact, & ! (in)
1774  vfact, & ! (in)
1775  kgrd, & ! (in)
1776  igrd, & ! (in)
1777  jgrd, & ! (in)
1778  ia, & ! (in)
1779  ja, & ! (in)
1780  ks, & ! (in)
1781  ke, & ! (in)
1782  logwegt ) ! (in)
1783  implicit none
1784 
1785  real(RP), intent(out) :: intp(:,:,:) ! interpolated data
1786 
1787  real(RP), intent(in) :: ref (:,:,:) ! reference data
1788  real(RP), intent(in) :: hfact(:,:,:) ! horizontal interp factor
1789  real(RP), intent(in) :: vfact(:,:,:,:,:) ! vertical interp factor
1790  integer, intent(in) :: kgrd (:,:,:,:,:) ! grid points of interp target
1791  integer, intent(in) :: igrd (:,:,:) ! grid points of interp target
1792  integer, intent(in) :: jgrd (:,:,:) ! grid points of interp target
1793  integer, intent(in) :: ia ! grid number of mine
1794  integer, intent(in) :: ja ! grid number of mine
1795  integer, intent(in) :: ks ! start grid number of mine
1796  integer, intent(in) :: ke ! end grid number of mine
1797 
1798  logical, intent(in), optional :: logwegt
1799 
1800  integer :: i, j, k
1801  logical :: logarithmic
1802  !---------------------------------------------------------------------------
1803 
1804  logarithmic = .false.
1805  if ( present(logwegt) ) then
1806  if ( logwegt ) then
1807  logarithmic = .true.
1808  endif
1809  endif
1810 
1811  ! linear interpolation
1812 !OCL PREFETCH
1813  do j = 1, ja
1814  do i = 1, ia
1815  do k = ks, ke
1816  intp(k,i,j) = ref(kgrd(k,i,j,1,1),igrd(i,j,1),jgrd(i,j,1)) &
1817  * hfact(i,j,1) * vfact(k,i,j,1,1) &
1818  + ref(kgrd(k,i,j,2,1),igrd(i,j,2),jgrd(i,j,2)) &
1819  * hfact(i,j,2) * vfact(k,i,j,2,1) &
1820  + ref(kgrd(k,i,j,3,1),igrd(i,j,3),jgrd(i,j,3)) &
1821  * hfact(i,j,3) * vfact(k,i,j,3,1) &
1822  + ref(kgrd(k,i,j,4,1),igrd(i,j,4),jgrd(i,j,4)) &
1823  * hfact(i,j,4) * vfact(k,i,j,4,1) &
1824  + ref(kgrd(k,i,j,5,1),igrd(i,j,5),jgrd(i,j,5)) &
1825  * hfact(i,j,5) * vfact(k,i,j,5,1) &
1826  + ref(kgrd(k,i,j,1,2),igrd(i,j,1),jgrd(i,j,1)) &
1827  * hfact(i,j,1) * vfact(k,i,j,1,2) &
1828  + ref(kgrd(k,i,j,2,2),igrd(i,j,2),jgrd(i,j,2)) &
1829  * hfact(i,j,2) * vfact(k,i,j,2,2) &
1830  + ref(kgrd(k,i,j,3,2),igrd(i,j,3),jgrd(i,j,3)) &
1831  * hfact(i,j,3) * vfact(k,i,j,3,2) &
1832  + ref(kgrd(k,i,j,4,2),igrd(i,j,4),jgrd(i,j,4)) &
1833  * hfact(i,j,4) * vfact(k,i,j,4,2) &
1834  + ref(kgrd(k,i,j,5,2),igrd(i,j,5),jgrd(i,j,5)) &
1835  * hfact(i,j,5) * vfact(k,i,j,5,2)
1836  end do
1837  end do
1838  end do
1839 
1840  ! logarithmic weighting (for pres, dens)
1841  if ( logarithmic ) then
1842  do j = 1, ja
1843  do i = 1, ia
1844  do k = ks, ke
1845  intp(k,i,j) = exp( intp(k,i,j) )
1846  end do
1847  end do
1848  end do
1849  endif
1850 
1851  return
1852  end subroutine intrpnest_interp_3d_5points
1853 
1854 
1855  !-----------------------------------------------------------------------------
1856  ! interpolation using eight-points for 2D data
1857  subroutine intrpnest_interp_2d_8points( &
1858  intp, & ! (out)
1859  ref, & ! (in)
1860  hfact, & ! (in)
1861  igrd, & ! (in)
1862  jgrd, & ! (in)
1863  ia, & ! (in)
1864  ja ) ! (in)
1865  implicit none
1866 
1867  real(RP), intent(out) :: intp(:,:) ! interpolated data
1868 
1869  real(RP), intent(in) :: ref (:,:) ! reference data
1870  real(RP), intent(in) :: hfact(:,:,:) ! horizontal interp factor
1871  integer, intent(in) :: igrd (:,:,:) ! grid points of interp target
1872  integer, intent(in) :: jgrd (:,:,:) ! grid points of interp target
1873  integer, intent(in) :: ia ! grid number of mine
1874  integer, intent(in) :: ja ! grid number of mine
1875 
1876  integer :: i, j
1877  !---------------------------------------------------------------------------
1878 
1879 !OCL PREFETCH
1880  do j = 1, ja
1881  do i = 1, ia
1882  intp(i,j) = ref(igrd(i,j,1),jgrd(i,j,1)) * hfact(i,j,1) &
1883  + ref(igrd(i,j,2),jgrd(i,j,2)) * hfact(i,j,2) &
1884  + ref(igrd(i,j,3),jgrd(i,j,3)) * hfact(i,j,3) &
1885  + ref(igrd(i,j,4),jgrd(i,j,4)) * hfact(i,j,4) &
1886  + ref(igrd(i,j,5),jgrd(i,j,5)) * hfact(i,j,5) &
1887  + ref(igrd(i,j,6),jgrd(i,j,6)) * hfact(i,j,6) &
1888  + ref(igrd(i,j,7),jgrd(i,j,7)) * hfact(i,j,7) &
1889  + ref(igrd(i,j,8),jgrd(i,j,8)) * hfact(i,j,8)
1890  end do
1891  end do
1892 
1893  return
1894  end subroutine intrpnest_interp_2d_8points
1895 
1896 
1897  !-----------------------------------------------------------------------------
1898  ! interpolation using eight-points for 3D data
1899  subroutine intrpnest_interp_3d_8points( &
1900  intp, & ! (out)
1901  ref, & ! (in)
1902  hfact, & ! (in)
1903  vfact, & ! (in)
1904  kgrd, & ! (in)
1905  igrd, & ! (in)
1906  jgrd, & ! (in)
1907  ia, & ! (in)
1908  ja, & ! (in)
1909  ks, & ! (in)
1910  ke, & ! (in)
1911  logwegt ) ! (in)
1912  implicit none
1913 
1914  real(RP), intent(out) :: intp(:,:,:) ! interpolated data
1915 
1916  real(RP), intent(in) :: ref (:,:,:) ! reference data
1917  real(RP), intent(in) :: hfact(:,:,:) ! horizontal interp factor
1918  real(RP), intent(in) :: vfact(:,:,:,:,:) ! vertical interp factor
1919  integer, intent(in) :: kgrd (:,:,:,:,:) ! grid points of interp target
1920  integer, intent(in) :: igrd (:,:,:) ! grid points of interp target
1921  integer, intent(in) :: jgrd (:,:,:) ! grid points of interp target
1922  integer, intent(in) :: ia ! grid number of mine
1923  integer, intent(in) :: ja ! grid number of mine
1924  integer, intent(in) :: ks ! start grid number of mine
1925  integer, intent(in) :: ke ! end grid number of mine
1926 
1927  logical, intent(in), optional :: logwegt
1928 
1929  integer :: i, j, k
1930  logical :: logarithmic
1931  !---------------------------------------------------------------------------
1932 
1933  logarithmic = .false.
1934  if ( present(logwegt) ) then
1935  if ( logwegt ) then
1936  logarithmic = .true.
1937  endif
1938  endif
1939 
1940  ! linear interpolation
1941 !OCL PREFETCH
1942  do j = 1, ja
1943  do i = 1, ia
1944  do k = ks, ke
1945  intp(k,i,j) = ref(kgrd(k,i,j,1,1),igrd(i,j,1),jgrd(i,j,1)) &
1946  * hfact(i,j,1) * vfact(k,i,j,1,1) &
1947  + ref(kgrd(k,i,j,2,1),igrd(i,j,2),jgrd(i,j,2)) &
1948  * hfact(i,j,2) * vfact(k,i,j,2,1) &
1949  + ref(kgrd(k,i,j,3,1),igrd(i,j,3),jgrd(i,j,3)) &
1950  * hfact(i,j,3) * vfact(k,i,j,3,1) &
1951  + ref(kgrd(k,i,j,4,1),igrd(i,j,4),jgrd(i,j,4)) &
1952  * hfact(i,j,4) * vfact(k,i,j,4,1) &
1953  + ref(kgrd(k,i,j,5,1),igrd(i,j,5),jgrd(i,j,5)) &
1954  * hfact(i,j,5) * vfact(k,i,j,5,1) &
1955  + ref(kgrd(k,i,j,6,1),igrd(i,j,6),jgrd(i,j,6)) &
1956  * hfact(i,j,6) * vfact(k,i,j,6,1) &
1957  + ref(kgrd(k,i,j,7,1),igrd(i,j,7),jgrd(i,j,7)) &
1958  * hfact(i,j,7) * vfact(k,i,j,7,1) &
1959  + ref(kgrd(k,i,j,8,1),igrd(i,j,8),jgrd(i,j,8)) &
1960  * hfact(i,j,8) * vfact(k,i,j,8,1) &
1961  + ref(kgrd(k,i,j,1,2),igrd(i,j,1),jgrd(i,j,1)) &
1962  * hfact(i,j,1) * vfact(k,i,j,1,2) &
1963  + ref(kgrd(k,i,j,2,2),igrd(i,j,2),jgrd(i,j,2)) &
1964  * hfact(i,j,2) * vfact(k,i,j,2,2) &
1965  + ref(kgrd(k,i,j,3,2),igrd(i,j,3),jgrd(i,j,3)) &
1966  * hfact(i,j,3) * vfact(k,i,j,3,2) &
1967  + ref(kgrd(k,i,j,4,2),igrd(i,j,3),jgrd(i,j,3)) &
1968  * hfact(i,j,4) * vfact(k,i,j,4,2) &
1969  + ref(kgrd(k,i,j,5,2),igrd(i,j,3),jgrd(i,j,3)) &
1970  * hfact(i,j,5) * vfact(k,i,j,5,2) &
1971  + ref(kgrd(k,i,j,6,2),igrd(i,j,3),jgrd(i,j,3)) &
1972  * hfact(i,j,6) * vfact(k,i,j,6,2) &
1973  + ref(kgrd(k,i,j,7,2),igrd(i,j,3),jgrd(i,j,3)) &
1974  * hfact(i,j,7) * vfact(k,i,j,7,2) &
1975  + ref(kgrd(k,i,j,8,2),igrd(i,j,8),jgrd(i,j,8)) &
1976  * hfact(i,j,8) * vfact(k,i,j,8,2)
1977  end do
1978  end do
1979  end do
1980 
1981  ! logarithmic weighting (for pres, dens)
1982  if ( logarithmic ) then
1983  do j = 1, ja
1984  do i = 1, ia
1985  do k = ks, ke
1986  intp(k,i,j) = exp( intp(k,i,j) )
1987  end do
1988  end do
1989  end do
1990  endif
1991 
1992  return
1993  end subroutine intrpnest_interp_3d_8points
1994 
1995 
1996  !-----------------------------------------------------------------------------
1997  ! interpolation using twelve-points for 2D data
1998  subroutine intrpnest_interp_2d_12points( &
1999  intp, & ! (out)
2000  ref, & ! (in)
2001  hfact, & ! (in)
2002  igrd, & ! (in)
2003  jgrd, & ! (in)
2004  ia, & ! (in)
2005  ja ) ! (in)
2006  implicit none
2007 
2008  real(RP), intent(out) :: intp(:,:) ! interpolated data
2009 
2010  real(RP), intent(in) :: ref (:,:) ! reference data
2011  real(RP), intent(in) :: hfact(:,:,:) ! horizontal interp factor
2012  integer, intent(in) :: igrd (:,:,:) ! grid points of interp target
2013  integer, intent(in) :: jgrd (:,:,:) ! grid points of interp target
2014  integer, intent(in) :: ia ! grid number of mine
2015  integer, intent(in) :: ja ! grid number of mine
2016 
2017  integer :: i, j
2018  !---------------------------------------------------------------------------
2019 
2020 !OCL PREFETCH
2021  do j = 1, ja
2022  do i = 1, ia
2023  intp(i,j) = ref(igrd(i,j,1), jgrd(i,j,1)) * hfact(i,j,1) &
2024  + ref(igrd(i,j,2), jgrd(i,j,2)) * hfact(i,j,2) &
2025  + ref(igrd(i,j,3), jgrd(i,j,3)) * hfact(i,j,3) &
2026  + ref(igrd(i,j,4), jgrd(i,j,4)) * hfact(i,j,4) &
2027  + ref(igrd(i,j,5), jgrd(i,j,5)) * hfact(i,j,5) &
2028  + ref(igrd(i,j,6), jgrd(i,j,6)) * hfact(i,j,6) &
2029  + ref(igrd(i,j,7), jgrd(i,j,7)) * hfact(i,j,7) &
2030  + ref(igrd(i,j,8), jgrd(i,j,8)) * hfact(i,j,8) &
2031  + ref(igrd(i,j,9), jgrd(i,j,9)) * hfact(i,j,9) &
2032  + ref(igrd(i,j,10),jgrd(i,j,10)) * hfact(i,j,10) &
2033  + ref(igrd(i,j,11),jgrd(i,j,11)) * hfact(i,j,11) &
2034  + ref(igrd(i,j,12),jgrd(i,j,12)) * hfact(i,j,12)
2035  end do
2036  end do
2037 
2038  return
2039  end subroutine intrpnest_interp_2d_12points
2040 
2041 
2042  !-----------------------------------------------------------------------------
2043  ! interpolation using twelve-points for 3D data
2044  subroutine intrpnest_interp_3d_12points( &
2045  intp, & ! (out)
2046  ref, & ! (in)
2047  hfact, & ! (in)
2048  vfact, & ! (in)
2049  kgrd, & ! (in)
2050  igrd, & ! (in)
2051  jgrd, & ! (in)
2052  ia, & ! (in)
2053  ja, & ! (in)
2054  ks, & ! (in)
2055  ke, & ! (in)
2056  logwegt ) ! (in)
2057  implicit none
2058 
2059  real(RP), intent(out) :: intp(:,:,:) ! interpolated data
2060 
2061  real(RP), intent(in) :: ref (:,:,:) ! reference data
2062  real(RP), intent(in) :: hfact(:,:,:) ! horizontal interp factor
2063  real(RP), intent(in) :: vfact(:,:,:,:,:) ! vertical interp factor
2064  integer, intent(in) :: kgrd (:,:,:,:,:) ! grid points of interp target
2065  integer, intent(in) :: igrd (:,:,:) ! grid points of interp target
2066  integer, intent(in) :: jgrd (:,:,:) ! grid points of interp target
2067  integer, intent(in) :: ia ! grid number of mine
2068  integer, intent(in) :: ja ! grid number of mine
2069  integer, intent(in) :: ks ! start grid number of mine
2070  integer, intent(in) :: ke ! end grid number of mine
2071 
2072  logical, intent(in), optional :: logwegt
2073 
2074  integer :: i, j, k
2075  logical :: logarithmic
2076  !---------------------------------------------------------------------------
2077 
2078  logarithmic = .false.
2079  if ( present(logwegt) ) then
2080  if ( logwegt ) then
2081  logarithmic = .true.
2082  endif
2083  endif
2084 
2085  ! linear interpolation
2086 !OCL PREFETCH
2087  do j = 1, ja
2088  do i = 1, ia
2089  do k = ks, ke
2090  intp(k,i,j) = ref(kgrd(k,i,j,1, 1),igrd(i,j,1 ),jgrd(i,j,1 )) &
2091  * hfact(i,j,1 ) * vfact(k,i,j,1, 1) &
2092  + ref(kgrd(k,i,j,2, 1),igrd(i,j,2 ),jgrd(i,j,2 )) &
2093  * hfact(i,j,2 ) * vfact(k,i,j,2, 1) &
2094  + ref(kgrd(k,i,j,3, 1),igrd(i,j,3 ),jgrd(i,j,3 )) &
2095  * hfact(i,j,3 ) * vfact(k,i,j,3, 1) &
2096  + ref(kgrd(k,i,j,4, 1),igrd(i,j,4 ),jgrd(i,j,4 )) &
2097  * hfact(i,j,4 ) * vfact(k,i,j,4, 1) &
2098  + ref(kgrd(k,i,j,5, 1),igrd(i,j,5 ),jgrd(i,j,5 )) &
2099  * hfact(i,j,5 ) * vfact(k,i,j,5, 1) &
2100  + ref(kgrd(k,i,j,6, 1),igrd(i,j,6 ),jgrd(i,j,6 )) &
2101  * hfact(i,j,6 ) * vfact(k,i,j,6, 1) &
2102  + ref(kgrd(k,i,j,7, 1),igrd(i,j,7 ),jgrd(i,j,7 )) &
2103  * hfact(i,j,7 ) * vfact(k,i,j,7, 1) &
2104  + ref(kgrd(k,i,j,8, 1),igrd(i,j,8 ),jgrd(i,j,8 )) &
2105  * hfact(i,j,8 ) * vfact(k,i,j,8, 1) &
2106  + ref(kgrd(k,i,j,9, 1),igrd(i,j,9 ),jgrd(i,j,9 )) &
2107  * hfact(i,j,9 ) * vfact(k,i,j,9, 1) &
2108  + ref(kgrd(k,i,j,10,1),igrd(i,j,10),jgrd(i,j,10)) &
2109  * hfact(i,j,10) * vfact(k,i,j,10,1) &
2110  + ref(kgrd(k,i,j,11,1),igrd(i,j,11),jgrd(i,j,11)) &
2111  * hfact(i,j,11) * vfact(k,i,j,11,1) &
2112  + ref(kgrd(k,i,j,12,1),igrd(i,j,12),jgrd(i,j,12)) &
2113  * hfact(i,j,12) * vfact(k,i,j,12,1) &
2114  + ref(kgrd(k,i,j,1, 2),igrd(i,j,1 ),jgrd(i,j,1 )) &
2115  * hfact(i,j,1 ) * vfact(k,i,j,1, 2) &
2116  + ref(kgrd(k,i,j,2, 2),igrd(i,j,2 ),jgrd(i,j,2 )) &
2117  * hfact(i,j,2 ) * vfact(k,i,j,2, 2) &
2118  + ref(kgrd(k,i,j,3, 2),igrd(i,j,3 ),jgrd(i,j,3 )) &
2119  * hfact(i,j,3 ) * vfact(k,i,j,3, 2) &
2120  + ref(kgrd(k,i,j,4, 2),igrd(i,j,4 ),jgrd(i,j,4 )) &
2121  * hfact(i,j,4 ) * vfact(k,i,j,4, 2) &
2122  + ref(kgrd(k,i,j,5, 2),igrd(i,j,5 ),jgrd(i,j,5 )) &
2123  * hfact(i,j,5 ) * vfact(k,i,j,5, 2) &
2124  + ref(kgrd(k,i,j,6, 2),igrd(i,j,6 ),jgrd(i,j,6 )) &
2125  * hfact(i,j,6 ) * vfact(k,i,j,6, 2) &
2126  + ref(kgrd(k,i,j,7, 2),igrd(i,j,7 ),jgrd(i,j,7 )) &
2127  * hfact(i,j,7 ) * vfact(k,i,j,7, 2) &
2128  + ref(kgrd(k,i,j,8, 2),igrd(i,j,8 ),jgrd(i,j,8 )) &
2129  * hfact(i,j,8 ) * vfact(k,i,j,8, 2) &
2130  + ref(kgrd(k,i,j,9, 2),igrd(i,j,9 ),jgrd(i,j,9 )) &
2131  * hfact(i,j,9 ) * vfact(k,i,j,9, 2) &
2132  + ref(kgrd(k,i,j,10,2),igrd(i,j,10),jgrd(i,j,10)) &
2133  * hfact(i,j,10) * vfact(k,i,j,10,2) &
2134  + ref(kgrd(k,i,j,11,2),igrd(i,j,11),jgrd(i,j,11)) &
2135  * hfact(i,j,11) * vfact(k,i,j,11,2) &
2136  + ref(kgrd(k,i,j,12,2),igrd(i,j,12),jgrd(i,j,12)) &
2137  * hfact(i,j,12) * vfact(k,i,j,12,2)
2138  end do
2139  end do
2140  end do
2141 
2142  ! logarithmic weighting (for pres, dens)
2143  if ( logarithmic ) then
2144  do j = 1, ja
2145  do i = 1, ia
2146  do k = ks, ke
2147  intp(k,i,j) = exp( intp(k,i,j) )
2148  end do
2149  end do
2150  end do
2151  endif
2152 
2153  return
2154  end subroutine intrpnest_interp_3d_12points
2155 
2156 
2157  !-----------------------------------------------------------------------------
2158  subroutine intrpnest_domain_compatibility( &
2159  lon_org, & ! (in)
2160  lat_org, & ! (in)
2161  lev_org, & ! (in)
2162  lon_loc, & ! (in)
2163  lat_loc, & ! (in)
2164  lev_loc, & ! (in)
2165  skip_x, & ! (in)
2166  skip_y, & ! (in)
2167  skip_z ) ! (in)
2168  use scale_process, only: &
2169  prc_mpistop
2170  use scale_const, only: &
2171  d2r => const_d2r
2172  implicit none
2173  real(RP), intent(in) :: lon_org(:,:)
2174  real(RP), intent(in) :: lat_org(:,:)
2175  real(RP), intent(in) :: lev_org(:,:,:)
2176  real(RP), intent(in) :: lon_loc(:,:)
2177  real(RP), intent(in) :: lat_loc(:,:)
2178  real(RP), intent(in) :: lev_loc(:,:,:)
2179  logical, intent(in), optional :: skip_x
2180  logical, intent(in), optional :: skip_y
2181  logical, intent(in), optional :: skip_z
2182 
2183  real(RP) :: max_ref, min_ref
2184  real(RP) :: max_loc, min_loc
2185 
2186  logical :: do_xdirec
2187  logical :: do_ydirec
2188  logical :: do_zdirec
2189 
2190  intrinsic size
2191  !---------------------------------------------------------------------------
2192 
2193  do_xdirec = .true.
2194  if ( present(skip_x) ) then
2195  if ( skip_x ) then
2196  do_xdirec = .false.
2197  endif
2198  endif
2199 
2200  do_ydirec = .true.
2201  if ( present(skip_y) ) then
2202  if ( skip_y ) then
2203  do_ydirec = .false.
2204  endif
2205  endif
2206 
2207  do_zdirec = .true.
2208  if ( present(skip_z) ) then
2209  if ( skip_z ) then
2210  do_zdirec = .false.
2211  endif
2212  endif
2213 
2214  if ( do_xdirec ) then
2215  max_ref = maxval( lon_org(:,:) / d2r )
2216  min_ref = minval( lon_org(:,:) / d2r )
2217  max_loc = maxval( lon_loc(:,:) / d2r )
2218  min_loc = minval( lon_loc(:,:) / d2r )
2219 
2220  if ( (min_ref+360.0_rp-max_ref) < 360.0_rp / size(lon_org,1) * 2.0_rp) then
2221  ! cyclic OK
2222  else if ( max_ref < max_loc .OR. min_ref > min_loc ) then
2223  write(*,*) 'xxx ERROR: REQUESTED DOMAIN IS TOO MUCH BROAD'
2224  write(*,*) 'xxx -- LONGITUDINAL direction over the limit'
2225  write(*,*) 'xxx -- reference max: ', max_ref
2226  write(*,*) 'xxx -- reference min: ', min_ref
2227  write(*,*) 'xxx -- local max: ', max_loc
2228  write(*,*) 'xxx -- local min: ', min_loc
2229  call prc_mpistop
2230  endif
2231  endif
2232 
2233  if ( do_ydirec ) then
2234  max_ref = maxval( lat_org(:,:) / d2r )
2235  min_ref = minval( lat_org(:,:) / d2r )
2236  max_loc = maxval( lat_loc(:,:) / d2r )
2237  min_loc = minval( lat_loc(:,:) / d2r )
2238 
2239  if ( max_ref < max_loc .OR. min_ref > min_loc ) then
2240  write(*,*) 'xxx ERROR: REQUESTED DOMAIN IS TOO MUCH BROAD'
2241  write(*,*) 'xxx -- LATITUDINAL direction over the limit'
2242  write(*,*) 'xxx -- reference max: ', max_ref
2243  write(*,*) 'xxx -- reference min: ', min_ref
2244  write(*,*) 'xxx -- local max: ', max_loc
2245  write(*,*) 'xxx -- local min: ', min_loc
2246  call prc_mpistop
2247  endif
2248  endif
2249 
2250  if ( do_zdirec ) then
2251  max_ref = maxval( lev_org(:,:,:) )
2252  !max_loc = maxval( lev_loc(KS-1:KE+1,:,:) ) ! HALO + 1
2253  max_loc = maxval( lev_loc(:,:,:) ) ! HALO + 1
2254  !min_ref = minval( lev_org(:,:,:) )
2255  !min_loc = minval( lev_loc(3:KA,:,:) ) ! HALO + 1
2256 
2257  if ( max_ref < max_loc ) then
2258  !if ( max_ref < max_loc .OR. min_ref > min_loc ) then
2259  write(*,*) 'xxx ERROR: REQUESTED DOMAIN IS TOO MUCH BROAD'
2260  write(*,*) 'xxx -- VERTICAL direction over the limit'
2261  write(*,*) 'xxx -- reference max: ', max_ref
2262  !write(*,*) 'xxx -- reference min: ', min_ref
2263  write(*,*) 'xxx -- local max: ', max_loc
2264  !write(*,*) 'xxx -- local min: ', min_loc
2265  call prc_mpistop
2266  endif
2267  endif
2268 
2269  return
2270  end subroutine intrpnest_domain_compatibility
2271 
2272  !-----------------------------------------------------------------------------
2273  ! Haversine Formula (from R.W. Sinnott, "Virtues of the Haversine",
2274  ! Sky and Telescope, vol. 68, no. 2, 1984, p. 159):
2275  function intrpnest_haversine( &
2276  la0, &
2277  lo0, &
2278  la, &
2279  lo ) &
2280  result( d )
2281  implicit none
2282  real(RP), intent(in) :: la0, lo0, la, lo ! la,la0: Lat, lo,lo0: Lon; [rad]
2283  real(RP) :: d, dlon, dlat, work1, work2
2284  !---------------------------------------------------------------------------
2285 
2286  ! output unit : [m]
2287  dlon = lo0 - lo
2288  dlat = la0 - la
2289  work1 = (sin(dlat/2.0_rp))**2.0_rp + &
2290  cos(la0) * cos(la) * (sin(dlon/2.0_rp))**2.0_rp
2291  work2 = 2.0_rp * asin(min( 1.0_rp, sqrt(work1) ))
2292  d = r_in_m * work2
2293 
2294  end function intrpnest_haversine
2295 
2296 end module scale_interpolation_nest
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_debug
module DEBUG
Definition: scale_debug.F90:13
scale_interpolation_nest::intrpnest_domain_compatibility
subroutine, public intrpnest_domain_compatibility(lon_org, lat_org, lev_org, lon_loc, lat_loc, lev_loc, skip_x, skip_y, skip_z)
Definition: scale_interpolation_nest.F90:2168
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_process::prc_mpistop
subroutine, public prc_mpistop
Abort MPI.
Definition: scale_process.F90:267
scale_stdio::io_l
logical, public io_l
output log or not? (this process)
Definition: scale_stdio.F90:61
scale_const::const_radius
real(rp), public const_radius
radius of the planet [m]
Definition: scale_const.F90:46
scale_stdio
module STDIO
Definition: scale_stdio.F90:12
scale_grid_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_grid_index.F90:59
scale_const::const_d2r
real(rp), public const_d2r
degree to radian
Definition: scale_const.F90:35
scale_interpolation_nest::intrpnest_interp_2d
procedure(intrpnest_intfc_interp_2d), pointer, public intrpnest_interp_2d
Definition: scale_interpolation_nest.F90:164
scale_interpolation_nest::intrpnest_interp_fact_latlon
subroutine, public intrpnest_interp_fact_latlon(hfact, igrd, jgrd, mylat, mylon, myIA, myJA, inlat, inlon, inIA, inJA)
Definition: scale_interpolation_nest.F90:324
scale_interpolation_nest::intrpnest_interp_fact_llz
subroutine, public intrpnest_interp_fact_llz(hfact, vfact, kgrd, igrd, jgrd, ncopy, myhgt, mylat, mylon, myKS, myKE, myIA, myJA, inhgt, inlat, inlon, inKA, inIA, inJA, landgrid)
Definition: scale_interpolation_nest.F90:397
scale_grid_index
module grid index
Definition: scale_grid_index.F90:14
scale_tracer
module TRACER
Definition: scale_tracer.F90:15
scale_index
module Index
Definition: scale_index.F90:14
scale_grid_index::ia
integer, public ia
of whole cells: x, local, with HALO
Definition: scale_grid_index.F90:55
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
scale_process
module PROCESS
Definition: scale_process.F90:14
scale_interpolation_nest::intrpnest_interp_3d
procedure(intrpnest_intfc_interp_3d), pointer, public intrpnest_interp_3d
Definition: scale_interpolation_nest.F90:199
scale_const
module CONSTANT
Definition: scale_const.F90:14
scale_grid_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_grid_index.F90:58
scale_grid_index::khalo
integer, parameter, public khalo
of halo cells: z
Definition: scale_grid_index.F90:50
scale_interpolation_nest
module INTERPOLATION (nesting system)
Definition: scale_interpolation_nest.F90:14
scale_prof
module profiler
Definition: scale_prof.F90:10
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:36
scale_precision
module PRECISION
Definition: scale_precision.F90:13
scale_land_grid_index
module land grid index
Definition: scale_land_grid_index.F90:10
scale_stdio::io_fid_log
integer, public io_fid_log
Log file ID.
Definition: scale_stdio.F90:56
scale_grid_index::ja
integer, public ja
of whole cells: y, local, with HALO
Definition: scale_grid_index.F90:56
scale_land_grid_index::lkmax
integer, public lkmax
Definition: scale_land_grid_index.F90:31
scale_interpolation_nest::intrpnest_setup
subroutine, public intrpnest_setup(interp_search_divnum, NEST_INTERP_LEVEL, NEST_INTERP_WEIGHT_ORDER, OFFLINE)
Setup.
Definition: scale_interpolation_nest.F90:233
Generated by   doxygen 1.8.14