SABER
tools_samp.F90
Go to the documentation of this file.
1 !----------------------------------------------------------------------
2 ! Module: tools_samp
3 !> Sampling functions
4 ! Author: Benjamin Menetrier
5 ! Licensing: this code is distributed under the CeCILL-C license
6 ! Copyright © 2015-... UCAR, CERFACS, METEO-FRANCE and IRIT
7 !----------------------------------------------------------------------
8 module tools_samp
9 
10 use atlas_module
11 use fckit_mpi_module, only: fckit_mpi_sum,fckit_mpi_status
12 use tools_const, only: pi,deg2rad
15 use tools_qsort, only: qsort
16 use tools_repro, only: repro,inf,sup,eq
17 use type_mpl, only: mpl_type
18 use type_rng, only: rng_type
19 use type_tree, only: tree_type
20 
21 implicit none
22 
23 private
25 
26 contains
27 
28 !----------------------------------------------------------------------
29 ! Subroutine: initialize_sampling
30 !> Intialize sampling
31 !----------------------------------------------------------------------
32 subroutine initialize_sampling(mpl,rng,area,n_loc,lon_loc,lat_loc,mask_loc,rh_loc,loc_to_glb,ntry,nrep,ns2_glb,sam2_glb, &
33  & fast,verbosity,n_uni,uni_to_loc,tree_uni)
34 
35 implicit none
36 
37 ! Passed variables
38 type(mpl_type),intent(inout) :: mpl !< MPI data
39 type(rng_type),intent(inout) :: rng !< Random number generator
40 real(kind_real),intent(in) :: area !< Global domain area
41 integer,intent(in) :: n_loc !< Number of points (local)
42 real(kind_real),intent(in) :: lon_loc(n_loc) !< Longitudes (local)
43 real(kind_real),intent(in) :: lat_loc(n_loc) !< Latitudes (local)
44 logical,intent(in) :: mask_loc(n_loc) !< Mask (local)
45 real(kind_real),intent(in) :: rh_loc(n_loc) !< Horizontal support radius (local)
46 integer,intent(in) :: loc_to_glb(n_loc) !< Local to global index
47 integer,intent(in) :: ntry !< Number of tries
48 integer,intent(in) :: nrep !< Number of replacements
49 integer,intent(in) :: ns2_glb !< Number of samplings points (global)
50 integer,intent(out) :: sam2_glb(ns2_glb) !< Horizontal sampling index (global)
51 logical,intent(in),optional :: fast !< Fast sampling flag
52 logical,intent(in),optional :: verbosity !< Verbosity flag
53 integer,intent(in),optional :: n_uni !< Universe size
54 integer,intent(in),optional :: uni_to_loc(:) !< Universe to local index
55 type(tree_type),intent(in),optional :: tree_uni !< Universe KD-tree
56 
57 ! Local variables
58 integer :: n_glb,n_loc_eff,n_glb_eff,i_glb,i_loc,is2_glb,ns1_glb_eff
59 integer,allocatable :: sam1_glb_eff(:),order(:)
60 real(kind_real),allocatable :: hash_glb(:),hash_loc(:)
61 real(kind_real),allocatable :: lon1_glb_eff(:),lat1_glb_eff(:),rh1_glb_eff(:)
62 real(kind_real),allocatable :: list(:)
63 logical :: lfast,lverbosity
64 logical,allocatable :: mask_glb(:)
65 character(len=1024),parameter :: subr = 'initialize_sampling'
66 
67 ! Local flags
68 lfast = .false.
69 lverbosity = .true.
70 if (present(fast)) lfast = fast
71 if (present(verbosity)) lverbosity = verbosity
72 
73 ! Global size
74 call mpl%f_comm%allreduce(n_loc,n_glb,fckit_mpi_sum())
75 
76 ! Number of effective points
77 n_loc_eff = count(mask_loc)
78 call mpl%f_comm%allreduce(n_loc_eff,n_glb_eff,fckit_mpi_sum())
79 
80 ! Check mask size
81 if (n_glb_eff==0) then
82  call mpl%abort(subr,'empty mask in initialize sampling')
83 elseif (n_glb_eff<ns2_glb) then
84  call mpl%abort(subr,'ns2_glb greater than n_glb_eff in initialize_sampling')
85 elseif (n_glb_eff==ns2_glb) then
86  write(mpl%info,'(a)') ' all points are used'
87  if (lverbosity) call mpl%flush
88 
89  ! Allocation
90  allocate(hash_loc(n_loc))
91  if (mpl%main) then
92  allocate(hash_glb(n_glb))
93  allocate(mask_glb(n_glb))
94  allocate(list(ns2_glb))
95  allocate(order(ns2_glb))
96  else
97  allocate(hash_glb(0))
98  allocate(mask_glb(0))
99  end if
100 
101  ! Compute hash
102  do i_loc=1,n_loc
103  hash_loc(i_loc) = lonlathash(lon_loc(i_loc),lat_loc(i_loc))
104  end do
105 
106  ! Communication
107  call mpl%loc_to_glb(n_loc,n_glb,loc_to_glb,hash_loc,hash_glb)
108  call mpl%loc_to_glb(n_loc,n_glb,loc_to_glb,mask_loc,mask_glb)
109 
110  if (mpl%main) then
111  ! Use all valid points
112  is2_glb = 0
113  do i_glb=1,n_glb
114  if (mask_glb(i_glb)) then
115  is2_glb = is2_glb+1
116  sam2_glb(is2_glb) = i_glb
117  list(is2_glb) = hash_glb(i_glb)
118  end if
119  end do
120 
121  ! Define points order
122  call qsort(ns2_glb,list,order)
123 
124  ! Reorder sampling
125  sam2_glb = sam2_glb(order)
126  end if
127 
128  ! Release memory
129  deallocate(hash_loc)
130  deallocate(hash_glb)
131  deallocate(mask_glb)
132  if (mpl%main) then
133  deallocate(list)
134  deallocate(order)
135  end if
136 else
137  ! First subsampling, local
138  if (present(n_uni).and.present(uni_to_loc).and.present(tree_uni)) then
139  call initialize_sampling_local(mpl,area,n_loc,n_loc_eff,n_glb_eff,lon_loc,lat_loc,mask_loc,rh_loc,loc_to_glb,ns2_glb, &
140  & ns1_glb_eff,lon1_glb_eff,lat1_glb_eff,rh1_glb_eff,sam1_glb_eff,lverbosity,n_uni,uni_to_loc,tree_uni)
141  else
142  call initialize_sampling_local(mpl,area,n_loc,n_loc_eff,n_glb_eff,lon_loc,lat_loc,mask_loc,rh_loc,loc_to_glb,ns2_glb, &
143  & ns1_glb_eff,lon1_glb_eff,lat1_glb_eff,rh1_glb_eff,sam1_glb_eff,lverbosity)
144  end if
145 
146  if (mpl%main) then
147  ! Second subsampling, global
148  call initialize_sampling_global(mpl,rng,ns1_glb_eff,lon1_glb_eff,lat1_glb_eff,rh1_glb_eff,sam1_glb_eff, &
149  & ntry,nrep,ns2_glb,sam2_glb,lfast,lverbosity)
150 
151  ! Release memory
152  deallocate(lon1_glb_eff)
153  deallocate(lat1_glb_eff)
154  deallocate(rh1_glb_eff)
155  deallocate(sam1_glb_eff)
156  end if
157 end if
158 
159 ! Broadcast
160 call mpl%f_comm%broadcast(sam2_glb,mpl%rootproc-1)
161 
162 end subroutine initialize_sampling
163 
164 !----------------------------------------------------------------------
165 ! Subroutine: initialize_sampling_local
166 !> Intialize sampling, local, based on ATLAS octahedral grid
167 !----------------------------------------------------------------------
168 subroutine initialize_sampling_local(mpl,area,n_loc,n_loc_eff,n_glb_eff,lon_loc,lat_loc,mask_loc,rh_loc,loc_to_glb,ns2_glb, &
169  & ns1_glb_eff,lon1_glb_eff,lat1_glb_eff,rh1_glb_eff,sam1_glb_eff,verbosity,n_uni,uni_to_loc,tree_uni)
170 
171 implicit none
172 
173 ! Passed variables
174 type(mpl_type),intent(inout) :: mpl !< MPI data
175 real(kind_real),intent(in) :: area !< Global domain area
176 integer,intent(in) :: n_loc !< Number of points (local)
177 integer,intent(in) :: n_loc_eff !< Number of points (local, effective)
178 integer,intent(in) :: n_glb_eff !< Number of points (global, effective)
179 real(kind_real),intent(in) :: lon_loc(n_loc) !< Longitudes (local)
180 real(kind_real),intent(in) :: lat_loc(n_loc) !< Latitudes (local)
181 logical,intent(in) :: mask_loc(n_loc) !< Mask (local)
182 real(kind_real),intent(in) :: rh_loc(n_loc) !< Horizontal support radius (local)
183 integer,intent(in) :: loc_to_glb(n_loc) !< Local to global index
184 integer,intent(in) :: ns2_glb !< Number of samplings points (global)
185 integer,intent(out) :: ns1_glb_eff
186 real(kind_real),allocatable,intent(out) :: lon1_glb_eff(:)
187 real(kind_real),allocatable,intent(out) :: lat1_glb_eff(:)
188 real(kind_real),allocatable,intent(out) :: rh1_glb_eff(:)
189 integer,allocatable,intent(out) :: sam1_glb_eff(:)
190 logical,intent(in),optional :: verbosity !< Verbosity flag
191 integer,intent(in),optional :: n_uni !< Universe size
192 integer,intent(in),optional :: uni_to_loc(:) !< Universe to local index
193 type(tree_type),intent(in),optional :: tree_uni !< Universe KD-tree
194 
195 ! Local variables
196 integer :: i_loc,i_loc_eff,n,ix,iy,iproc,is1_glb,ns1_loc,is1_loc,nfac,ns1_loc_tmp,nn_index(1),ns1_glb,is1_glb_eff,offset
197 integer,allocatable :: s1_loc_to_glb(:),s1_loc_to_glb_tmp(:),sam1_loc(:),sam1_loc_tmp(:),sam1_glb(:)
198 real(kind_real) :: lonlat(2),nn_dist(1)
199 real(kind_real),allocatable :: lon1_loc(:),lat1_loc(:),dist1_loc(:),rh1_loc(:)
200 real(kind_real),allocatable :: lon1_loc_tmp(:),lat1_loc_tmp(:),dist1_loc_tmp(:),rh1_loc_tmp(:)
201 real(kind_real),allocatable :: lon1_glb(:),lat1_glb(:),dist1_glb(:),rh1_glb(:)
202 logical :: lfull_grid,lverbosity,update,retry
203 logical,allocatable :: lmask(:)
204 character(len=6) :: gridid
205 character(len=1024),parameter :: subr = 'initialize_sampling_local'
206 type(fckit_mpi_status) :: status
207 type(atlas_structuredgrid) :: agrid
208 
209 ! Local flags
210 lverbosity = .true.
211 if (present(verbosity)) lverbosity = verbosity
212 lfull_grid = present(n_uni).and.present(uni_to_loc).and.present(tree_uni)
213 
214 ! Initialization
215 nfac = 1
216 retry = .true.
217 
218 do while (retry)
219  ! Update nfac
220  nfac = 2*nfac
221 
222  if (lfull_grid) then
223  ! Number of required points
224  ns1_glb = nfac*int(real(ns2_glb,kind_real)*4.0*pi/area)
225 
226  ! Octahedral grid
227  n = int(-4.5+sqrt(20.25+0.25*real(ns1_glb,kind_real)))+1
228  write(gridid,'(a,i5.5)') 'O',n
229  agrid = atlas_structuredgrid(gridid)
230  ns1_glb = 4*n**2+36*n
231 
232  ! Allocation
233  allocate(lon1_loc(ns1_glb))
234  allocate(lat1_loc(ns1_glb))
235  allocate(rh1_loc(ns1_glb))
236  allocate(sam1_loc(ns1_glb))
237  allocate(dist1_loc(ns1_glb))
238  allocate(s1_loc_to_glb(ns1_glb))
239  allocate(lmask(n_loc))
240 
241  ! Initialization
242  lmask = mask_loc
243  if (lverbosity) call mpl%prog_init(ns1_glb)
244 
245  ! Loop over octahedral grid points
246  is1_loc = 0
247  is1_glb = 0
248  do iy=1,int(agrid%ny(),kind_int)
249  do ix=1,int(agrid%nx(iy),kind_int)
250  ! Global index
251  is1_glb = is1_glb+1
252 
253  ! Get longitude/latitude
254  lonlat = agrid%lonlat(ix,iy)*deg2rad
255 
256  ! Find nearest neighbor in universe
257  call tree_uni%find_nearest_neighbors(lonlat(1),lonlat(2),1,nn_index,nn_dist)
258  i_loc = uni_to_loc(nn_index(1))
259 
260  if (mpl%msv%isnot(i_loc)) then
261  ! Keep valid points
262  if (lmask(i_loc)) then
263  is1_loc = is1_loc+1
264  lon1_loc(is1_loc) = lon_loc(i_loc)
265  lat1_loc(is1_loc) = lat_loc(i_loc)
266  rh1_loc(is1_loc) = rh_loc(i_loc)
267  sam1_loc(is1_loc) = loc_to_glb(i_loc)
268  dist1_loc(is1_loc) = nn_dist(1)
269  s1_loc_to_glb(is1_loc) = is1_glb
270  lmask(i_loc) = .false.
271  end if
272  end if
273 
274  ! Update
275  if (lverbosity) call mpl%prog_print(is1_glb)
276  end do
277  end do
278  if (lverbosity) call mpl%prog_final(.false.)
279 
280  ! Number of local valid points
281  ns1_loc = is1_loc
282 
283  ! Release memory
284  deallocate(lmask)
285  else
286  ! All point are used
287  ns1_glb = n_glb_eff
288  ns1_loc = n_loc_eff
289 
290  ! Allocation
291  allocate(lon1_loc(ns1_loc))
292  allocate(lat1_loc(ns1_loc))
293  allocate(rh1_loc(ns1_loc))
294  allocate(sam1_loc(ns1_loc))
295 
296  ! Copy
297  i_loc_eff = 0
298  do i_loc=1,n_loc
299  if (mask_loc(i_loc)) then
300  i_loc_eff = i_loc_eff+1
301  lon1_loc(i_loc_eff) = lon_loc(i_loc)
302  lat1_loc(i_loc_eff) = lat_loc(i_loc)
303  rh1_loc(i_loc_eff) = rh_loc(i_loc)
304  sam1_loc(i_loc_eff) = loc_to_glb(i_loc)
305  end if
306  end do
307  end if
308 
309  if (mpl%main) then
310  if (lfull_grid) then
311  ! Continue printing
312  write(mpl%info,'(a)') ' => '
313  if (lverbosity) call mpl%flush(.false.)
314  end if
315 
316  ! Allocation
317  allocate(lon1_glb(ns1_glb))
318  allocate(lat1_glb(ns1_glb))
319  allocate(rh1_glb(ns1_glb))
320  allocate(sam1_glb(ns1_glb))
321  if (lfull_grid) allocate(dist1_glb(ns1_glb))
322 
323  ! Initialization
324  lon1_glb = mpl%msv%valr
325  lat1_glb = mpl%msv%valr
326  rh1_glb = mpl%msv%valr
327  sam1_glb = mpl%msv%vali
328  if (lfull_grid) then
329  dist1_glb = mpl%msv%valr
330  else
331  offset = 0
332  end if
333 
334  ! Receive and copy data on rootproc
335  do iproc=1,mpl%nproc
336  if (iproc==mpl%rootproc) then
337  ! Copy dimension
338  ns1_loc_tmp = ns1_loc
339  else
340  ! Receive dimension
341  call mpl%f_comm%receive(ns1_loc_tmp,iproc-1,mpl%tag,status)
342  end if
343 
344  ! Allocation
345  allocate(lon1_loc_tmp(ns1_loc_tmp))
346  allocate(lat1_loc_tmp(ns1_loc_tmp))
347  allocate(rh1_loc_tmp(ns1_loc_tmp))
348  allocate(sam1_loc_tmp(ns1_loc_tmp))
349  if (lfull_grid) then
350  allocate(dist1_loc_tmp(ns1_loc_tmp))
351  allocate(s1_loc_to_glb_tmp(ns1_loc_tmp))
352  end if
353 
354  if (iproc==mpl%rootproc) then
355  ! Copy data
356  lon1_loc_tmp = lon1_loc(1:ns1_loc)
357  lat1_loc_tmp = lat1_loc(1:ns1_loc)
358  rh1_loc_tmp = rh1_loc(1:ns1_loc)
359  sam1_loc_tmp = sam1_loc(1:ns1_loc)
360  if (lfull_grid) then
361  dist1_loc_tmp = dist1_loc(1:ns1_loc)
362  s1_loc_to_glb_tmp = s1_loc_to_glb(1:ns1_loc)
363  end if
364  else
365  ! Receive data
366  call mpl%f_comm%receive(lon1_loc_tmp,iproc-1,mpl%tag+1,status)
367  call mpl%f_comm%receive(lat1_loc_tmp,iproc-1,mpl%tag+2,status)
368  call mpl%f_comm%receive(rh1_loc_tmp,iproc-1,mpl%tag+3,status)
369  call mpl%f_comm%receive(sam1_loc_tmp,iproc-1,mpl%tag+4,status)
370  if (lfull_grid) then
371  call mpl%f_comm%receive(dist1_loc_tmp,iproc-1,mpl%tag+5,status)
372  call mpl%f_comm%receive(s1_loc_to_glb_tmp,iproc-1,mpl%tag+6,status)
373  end if
374  end if
375 
376  ! Update global array
377  if (lfull_grid) then
378  do is1_loc=1,ns1_loc_tmp
379  ! Global index
380  is1_glb = s1_loc_to_glb_tmp(is1_loc)
381 
382  ! Check the global array status
383  update = mpl%msv%is(dist1_glb(is1_glb))
384  if (.not.update) update = (dist1_loc_tmp(is1_loc)<dist1_glb(is1_glb))
385  if (update) then
386  lon1_glb(is1_glb) = lon1_loc_tmp(is1_loc)
387  lat1_glb(is1_glb) = lat1_loc_tmp(is1_loc)
388  rh1_glb(is1_glb) = rh1_loc_tmp(is1_loc)
389  sam1_glb(is1_glb) = sam1_loc_tmp(is1_loc)
390  dist1_glb(is1_glb) = dist1_loc_tmp(is1_loc)
391  end if
392  end do
393  else
394  ! Copy with offset
395  lon1_glb(offset+1:offset+ns1_loc_tmp) = lon1_loc_tmp
396  lat1_glb(offset+1:offset+ns1_loc_tmp) = lat1_loc_tmp
397  rh1_glb(offset+1:offset+ns1_loc_tmp) = rh1_loc_tmp
398  sam1_glb(offset+1:offset+ns1_loc_tmp) = sam1_loc_tmp
399  offset = offset+ns1_loc_tmp
400  end if
401 
402  ! Release memory
403  deallocate(lon1_loc_tmp)
404  deallocate(lat1_loc_tmp)
405  deallocate(rh1_loc_tmp)
406  deallocate(sam1_loc_tmp)
407  if (lfull_grid) then
408  deallocate(dist1_loc_tmp)
409  deallocate(s1_loc_to_glb_tmp)
410  end if
411  end do
412 
413  ! Number of global valid points
414  if (lfull_grid) then
415  ns1_glb_eff = count(mpl%msv%isnot(dist1_glb))
416  else
417  ns1_glb_eff = ns1_glb
418  end if
419  retry = (ns1_glb_eff<ns2_glb)
420  else
421  ! Send data to rootproc
422  call mpl%f_comm%send(ns1_loc,mpl%rootproc-1,mpl%tag)
423  call mpl%f_comm%send(lon1_loc(1:ns1_loc),mpl%rootproc-1,mpl%tag+1)
424  call mpl%f_comm%send(lat1_loc(1:ns1_loc),mpl%rootproc-1,mpl%tag+2)
425  call mpl%f_comm%send(rh1_loc(1:ns1_loc),mpl%rootproc-1,mpl%tag+3)
426  call mpl%f_comm%send(sam1_loc(1:ns1_loc),mpl%rootproc-1,mpl%tag+4)
427  if (lfull_grid) then
428  call mpl%f_comm%send(dist1_loc(1:ns1_loc),mpl%rootproc-1,mpl%tag+5)
429  call mpl%f_comm%send(s1_loc_to_glb(1:ns1_loc),mpl%rootproc-1,mpl%tag+6)
430  end if
431 
432  ! Number of global valid points
433  ns1_glb_eff = 0
434 
435  ! Stop printing
436  write(mpl%info,'(a)') ''
437  if (lverbosity) call mpl%flush
438  end if
439  if (lfull_grid) then
440  call mpl%update_tag(7)
441  else
442  call mpl%update_tag(5)
443  end if
444 
445  ! Broadcast
446  call mpl%f_comm%broadcast(retry,mpl%rootproc-1)
447 
448  ! Check situation
449  if ((.not.lfull_grid).and.retry) call mpl%abort(subr,'retry activated for limited grid')
450  if ((nfac>=8).and.retry) call mpl%abort(subr,'retry activated for nfac>=8')
451 
452  ! Release memory
453  deallocate(lon1_loc)
454  deallocate(lat1_loc)
455  deallocate(rh1_loc)
456  deallocate(sam1_loc)
457  if (lfull_grid) then
458  deallocate(dist1_loc)
459  deallocate(s1_loc_to_glb)
460  end if
461  if (mpl%main.and.retry) then
462  deallocate(lon1_glb)
463  deallocate(lat1_glb)
464  deallocate(rh1_glb)
465  deallocate(sam1_glb)
466  if (lfull_grid) deallocate(dist1_glb)
467  end if
468 end do
469 
470 if (mpl%main) then
471  ! Allocation
472  allocate(lon1_glb_eff(ns1_glb_eff))
473  allocate(lat1_glb_eff(ns1_glb_eff))
474  allocate(rh1_glb_eff(ns1_glb_eff))
475  allocate(sam1_glb_eff(ns1_glb_eff))
476 
477  ! Copy valid data
478  if (lfull_grid) then
479  is1_glb_eff = 0
480  do is1_glb=1,ns1_glb
481  if (mpl%msv%isnot(dist1_glb(is1_glb))) then
482  is1_glb_eff = is1_glb_eff+1
483  lon1_glb_eff(is1_glb_eff) = lon1_glb(is1_glb)
484  lat1_glb_eff(is1_glb_eff) = lat1_glb(is1_glb)
485  rh1_glb_eff(is1_glb_eff) = rh1_glb(is1_glb)
486  sam1_glb_eff(is1_glb_eff) = sam1_glb(is1_glb)
487  end if
488  end do
489  else
490  lon1_glb_eff = lon1_glb
491  lat1_glb_eff = lat1_glb
492  rh1_glb_eff = rh1_glb
493  sam1_glb_eff = sam1_glb
494  end if
495 
496  ! Release memory
497  deallocate(lon1_glb)
498  deallocate(lat1_glb)
499  deallocate(rh1_glb)
500  deallocate(sam1_glb)
501  if (lfull_grid) deallocate(dist1_glb)
502 end if
503 
504 end subroutine initialize_sampling_local
505 
506 
507 !----------------------------------------------------------------------
508 ! Subroutine: initialize_sampling_global
509 !> Intialize sampling, global
510 !----------------------------------------------------------------------
511 subroutine initialize_sampling_global(mpl,rng,ns1_glb_eff,lon1_glb_eff,lat1_glb_eff,rh1_glb_eff,sam1_glb_eff,ntry,nrep, &
512  & ns2_glb,sam2_glb,fast,verbosity)
513 
514 implicit none
515 
516 ! Passed variables
517 type(mpl_type),intent(inout) :: mpl !< MPI data
518 type(rng_type),intent(inout) :: rng !< Random number generator
519 integer,intent(in) :: ns1_glb_eff
520 real(kind_real),intent(in) :: lon1_glb_eff(ns1_glb_eff)
521 real(kind_real),intent(in) :: lat1_glb_eff(ns1_glb_eff)
522 real(kind_real),intent(in) :: rh1_glb_eff(ns1_glb_eff)
523 integer,intent(in) :: sam1_glb_eff(ns1_glb_eff)
524 integer,intent(in) :: ntry !< Number of tries
525 integer,intent(in) :: nrep !< Number of replacements
526 integer,intent(in) :: ns2_glb !< Number of samplings points (global)
527 integer,intent(out) :: sam2_glb(ns2_glb) !< Horizontal sampling index (global)
528 logical,intent(in),optional :: fast !< Fast sampling flag
529 logical,intent(in),optional :: verbosity !< Verbosity flag
530 
531 ! Local variables
532 integer :: is2_glb,js,irep,irmax,itry,ir,irval,irvalmin,irvalmax,is2_glb_min,nrep_eff,nn_index(2),is1_glb_eff,ns1_glb_val
533 integer :: order(ns1_glb_eff),sam1_glb_tmp(ns1_glb_eff)
534 integer,allocatable :: to_valid(:),sam2_glb_tmp(:)
535 real(kind_real) :: d,distmax,distmin,nn_dist(2),cdf_norm,rr
536 real(kind_real) :: list(ns1_glb_eff),lon1_glb_tmp(ns1_glb_eff),lat1_glb_tmp(ns1_glb_eff),rh1_glb_tmp(ns1_glb_eff)
537 real(kind_real),allocatable :: cdf(:),lon_rep(:),lat_rep(:),dist(:)
538 logical :: lfast,lverbosity
539 logical,allocatable :: lmask(:),smask(:),rmask(:)
540 character(len=1024),parameter :: subr = 'initialize_sampling_global'
541 type(tree_type) :: tree
542 
543 ! Local flags
544 lfast = .false.
545 lverbosity = .true.
546 if (present(fast)) lfast = fast
547 if (present(verbosity)) lverbosity = verbosity
548 
549 ! Define points order
550 do is1_glb_eff=1,ns1_glb_eff
551  list(is1_glb_eff) = lonlathash(lon1_glb_eff(is1_glb_eff),lat1_glb_eff(is1_glb_eff))
552 end do
553 call qsort(ns1_glb_eff,list,order)
554 
555 ! Reorder data
556 lon1_glb_tmp = lon1_glb_eff(order)
557 lat1_glb_tmp = lat1_glb_eff(order)
558 rh1_glb_tmp = rh1_glb_eff(order)
559 sam1_glb_tmp = sam1_glb_eff(order)
560 
561 ! Allocation
562 nrep_eff = min(nrep,ns1_glb_eff-ns2_glb)
563 allocate(sam2_glb_tmp(ns2_glb+nrep_eff))
564 allocate(lmask(ns1_glb_eff))
565 allocate(smask(ns1_glb_eff))
566 allocate(to_valid(ns1_glb_eff))
567 
568 ! Initialization
569 sam2_glb_tmp = mpl%msv%vali
570 lmask = .true.
571 smask = .false.
572 to_valid = mpl%msv%vali
573 do is1_glb_eff=1,ns1_glb_eff
574  to_valid(is1_glb_eff) = is1_glb_eff
575 end do
576 ns1_glb_val = ns1_glb_eff
577 if (lverbosity) call mpl%prog_init(ns2_glb+nrep_eff)
578 
579 if (lfast) then
580  ! Define sampling with a cumulative distribution function
581 
582  ! Allocation
583  allocate(cdf(ns1_glb_eff))
584 
585  ! Initialization
586  cdf(1) = 0.0
587  do is1_glb_eff=2,ns1_glb_eff
588  if (lmask(is1_glb_eff)) then
589  cdf(is1_glb_eff) = cdf(is1_glb_eff-1)+1.0/rh1_glb_tmp(is1_glb_eff)**2
590  end if
591  end do
592  cdf_norm = 1.0/cdf(ns1_glb_eff)
593  cdf = cdf*cdf_norm
594 
595  do is2_glb=1,ns2_glb+nrep_eff
596  ! Generate random number
597  call rng%rand_real(0.0_kind_real,1.0_kind_real,rr)
598 
599  ! Dichotomy to find the value
600  irvalmin = 1
601  irvalmax = ns1_glb_val
602  do while (irvalmax-irvalmin>1)
603  irval = (irvalmin+irvalmax)/2
604  if ((sup(cdf(irvalmin),rr).and.sup(cdf(irval),rr)).or.(inf(cdf(irvalmin),rr).and.inf(cdf(irval),rr))) then
605  irvalmin = irval
606  else
607  irvalmax = irval
608  end if
609  end do
610 
611  ! New sampling point
612  ir = to_valid(irval)
613  sam2_glb_tmp(is2_glb) = ir
614 
615  ! Shift valid points array
616  if (irval<ns1_glb_val) then
617  cdf(irval:ns1_glb_val-1) = cdf(irval+1:ns1_glb_val)
618  to_valid(irval:ns1_glb_val-1) = to_valid(irval+1:ns1_glb_val)
619  end if
620  ns1_glb_val = ns1_glb_val-1
621 
622  ! Renormalize cdf
623  cdf_norm = 1.0/cdf(ns1_glb_val)
624  cdf(1:ns1_glb_val) = cdf(1:ns1_glb_val)*cdf_norm
625 
626  ! Update
627  if (lverbosity) call mpl%prog_print(is2_glb)
628  end do
629  if (lverbosity) call mpl%prog_final(.false.)
630 
631  ! Release memory
632  deallocate(cdf)
633 else
634  ! Define sampling with a KD-tree
635  do is2_glb=1,ns2_glb+nrep_eff
636  if (is2_glb>2) then
637  ! Allocation
638  call tree%alloc(mpl,ns1_glb_eff,mask=smask)
639 
640  ! Initialization
641  call tree%init(lon1_glb_tmp,lat1_glb_tmp)
642  end if
643 
644  ! Initialization
645  distmax = 0.0
646  irmax = 0
647  irvalmax = 0
648  itry = 1
649 
650  ! Find a new point
651  do itry=1,ntry
652  ! Generate a random index among valid points
653  call rng%rand_integer(1,ns1_glb_val,irval)
654  ir = to_valid(irval)
655 
656  ! Check point validity
657  if (is2_glb==1) then
658  ! Accept point
659  irvalmax = irval
660  irmax = ir
661  else
662  if (is2_glb==2) then
663  ! Compute distance
664  nn_index(1) = sam2_glb_tmp(1)
665  call sphere_dist(lon1_glb_tmp(ir),lat1_glb_tmp(ir),lon1_glb_tmp(nn_index(1)),lat1_glb_tmp(nn_index(1)), &
666  & nn_dist(1))
667  else
668  ! Find nearest neighbor distance
669  call tree%find_nearest_neighbors(lon1_glb_tmp(ir),lat1_glb_tmp(ir),1,nn_index(1:1),nn_dist(1:1))
670  end if
671  d = nn_dist(1)**2/(rh1_glb_tmp(ir)**2+rh1_glb_tmp(nn_index(1))**2)
672 
673  ! Check distance
674  if (sup(d,distmax)) then
675  distmax = d
676  irvalmax = irval
677  irmax = ir
678  end if
679  end if
680  end do
681 
682  ! Delete tree
683  if (is2_glb>2) call tree%dealloc
684 
685  ! Add point to sampling
686  if (irmax>0) then
687  ! New sampling point
688  sam2_glb_tmp(is2_glb) = irmax
689  lmask(irmax) = .false.
690  smask(irmax) = .true.
691 
692  ! Shift valid points array
693  if (irvalmax<ns1_glb_val) to_valid(irvalmax:ns1_glb_val-1) = to_valid(irvalmax+1:ns1_glb_val)
694  ns1_glb_val = ns1_glb_val-1
695  end if
696 
697  ! Update
698  if (lverbosity) call mpl%prog_print(is2_glb)
699  end do
700  if (lverbosity) call mpl%prog_final(.false.)
701 end if
702 
703 if (nrep_eff>0) then
704  ! Continue printing
705  write(mpl%info,'(a)') ' => '
706  if (lverbosity) call mpl%flush(.false.)
707 
708  ! Allocation
709  allocate(rmask(ns2_glb+nrep_eff))
710  allocate(lon_rep(ns2_glb+nrep_eff))
711  allocate(lat_rep(ns2_glb+nrep_eff))
712  allocate(dist(ns2_glb+nrep_eff))
713 
714  ! Initialization
715  rmask = .true.
716  do is2_glb=1,ns2_glb+nrep_eff
717  lon_rep(is2_glb) = lon1_glb_tmp(sam2_glb_tmp(is2_glb))
718  lat_rep(is2_glb) = lat1_glb_tmp(sam2_glb_tmp(is2_glb))
719  end do
720  dist = mpl%msv%valr
721  if (lverbosity) call mpl%prog_init(nrep_eff)
722 
723  ! Remove closest points
724  do irep=1,nrep_eff
725  ! Allocation
726  call tree%alloc(mpl,ns2_glb+nrep_eff,mask=rmask)
727 
728  ! Initialization
729  call tree%init(lon_rep,lat_rep)
730 
731  ! Get minimum distance
732  do is2_glb=1,ns2_glb+nrep_eff
733  if (rmask(is2_glb)) then
734  ! Find nearest neighbor distance
735  call tree%find_nearest_neighbors(lon1_glb_tmp(sam2_glb_tmp(is2_glb)),lat1_glb_tmp(sam2_glb_tmp(is2_glb)), &
736  & 2,nn_index,nn_dist)
737  if (nn_index(1)==is2_glb) then
738  dist(is2_glb) = nn_dist(2)
739  elseif (nn_index(2)==is2_glb) then
740  dist(is2_glb) = nn_dist(1)
741  else
742  call mpl%abort(subr,'wrong index in replacement')
743  end if
744  dist(is2_glb) = dist(is2_glb)**2/(rh1_glb_tmp(sam2_glb_tmp(nn_index(1)))**2 &
745  & +rh1_glb_tmp(sam2_glb_tmp(nn_index(2)))**2)
746  end if
747  end do
748 
749  ! Delete tree
750  call tree%dealloc
751 
752  ! Remove worst point
753  distmin = huge_real
754  is2_glb_min = mpl%msv%vali
755  do is2_glb=1,ns2_glb+nrep_eff
756  if (rmask(is2_glb)) then
757  if (inf(dist(is2_glb),distmin)) then
758  is2_glb_min = is2_glb
759  distmin = dist(is2_glb)
760  end if
761  end if
762  end do
763  rmask(is2_glb_min) = .false.
764 
765  ! Update
766  if (lverbosity) call mpl%prog_print(irep)
767  end do
768  if (lverbosity) call mpl%prog_final
769 
770  ! Copy sam2_glb
771  js = 0
772  do is2_glb=1,ns2_glb+nrep_eff
773  if (rmask(is2_glb)) then
774  js = js+1
775  sam2_glb(js) = sam2_glb_tmp(is2_glb)
776  end if
777  end do
778 
779  ! Release memory
780  deallocate(rmask)
781  deallocate(lon_rep)
782  deallocate(lat_rep)
783  deallocate(dist)
784 else
785  ! Stop printing
786  write(mpl%info,'(a)') ''
787  if (lverbosity) call mpl%flush
788 
789  ! Copy sam2_glb
790  sam2_glb = sam2_glb_tmp
791 end if
792 
793 ! Apply first sampling step
794 sam2_glb = sam1_glb_tmp(sam2_glb)
795 
796 ! Release memory
797 deallocate(sam2_glb_tmp)
798 deallocate(lmask)
799 deallocate(smask)
800 deallocate(to_valid)
801 
802 end subroutine initialize_sampling_global
803 
804 end module tools_samp
tools_func::sphere_dist
subroutine, public sphere_dist(lon_i, lat_i, lon_f, lat_f, dist)
Compute the great-circle distance between two points.
Definition: tools_func.F90:126
tools_repro::sup
logical function, public sup(x, y)
Superior test for reals.
Definition: tools_repro.F90:92
tools_samp::initialize_sampling_local
subroutine, public initialize_sampling_local(mpl, area, n_loc, n_loc_eff, n_glb_eff, lon_loc, lat_loc, mask_loc, rh_loc, loc_to_glb, ns2_glb, ns1_glb_eff, lon1_glb_eff, lat1_glb_eff, rh1_glb_eff, sam1_glb_eff, verbosity, n_uni, uni_to_loc, tree_uni)
Intialize sampling, local, based on ATLAS octahedral grid.
Definition: tools_samp.F90:170
tools_repro::inf
logical function, public inf(x, y)
Inferior test for reals.
Definition: tools_repro.F90:53
type_rng::rng_type
Definition: type_rng.F90:22
tools_repro::repro
logical, public repro
Definition: tools_repro.F90:16
tools_const
Define usual constants and missing values.
Definition: tools_const.F90:8
tools_func
Usual functions.
Definition: tools_func.F90:8
tools_qsort
Qsort routines.
Definition: tools_qsort.F90:12
type_rng
Random numbers generator derived type.
Definition: type_rng.F90:8
tools_samp::initialize_sampling
subroutine, public initialize_sampling(mpl, rng, area, n_loc, lon_loc, lat_loc, mask_loc, rh_loc, loc_to_glb, ntry, nrep, ns2_glb, sam2_glb, fast, verbosity, n_uni, uni_to_loc, tree_uni)
Intialize sampling.
Definition: tools_samp.F90:34
type_tree::tree_type
Definition: type_tree.F90:22
tools_repro
Reproducibility functions.
Definition: tools_repro.F90:8
tools_samp::initialize_sampling_global
subroutine, public initialize_sampling_global(mpl, rng, ns1_glb_eff, lon1_glb_eff, lat1_glb_eff, rh1_glb_eff, sam1_glb_eff, ntry, nrep, ns2_glb, sam2_glb, fast, verbosity)
Intialize sampling, global.
Definition: tools_samp.F90:513
tools_repro::eq
logical function, public eq(x, y)
Equal test for reals.
Definition: tools_repro.F90:30
tools_samp
Sampling functions.
Definition: tools_samp.F90:8
tools_kinds::huge_real
real(kind_real), parameter, public huge_real
Definition: tools_kinds.F90:25
tools_func::lonlathash
real(kind_real) function, public lonlathash(lon, lat, il)
Define a unique real from a lon/lat pair.
Definition: tools_func.F90:96
tools_qsort::qsort
Definition: tools_qsort.F90:19
tools_kinds
Kinds definition.
Definition: tools_kinds.F90:8
type_tree
Tree derived type.
Definition: type_tree.F90:8
tools_const::deg2rad
real(kind_real), parameter, public deg2rad
Definition: tools_const.F90:15
type_mpl
MPI parameters derived type.
Definition: type_mpl.F90:8
tools_func::lonlatmod
subroutine, public lonlatmod(lon, lat)
Set latitude between -pi/2 and pi/2 and longitude between -pi and pi.
Definition: tools_func.F90:66
tools_const::pi
real(kind_real), parameter, public pi
Definition: tools_const.F90:14
type_mpl::mpl_type
Definition: type_mpl.F90:24
tools_kinds::kind_real
integer, parameter, public kind_real
Definition: tools_kinds.F90:18
tools_kinds::kind_int
integer, parameter, public kind_int
Definition: tools_kinds.F90:16