11 use fckit_mpi_module,
only: fckit_mpi_sum,fckit_mpi_status
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)
41 integer,
intent(in) :: n_loc
42 real(
kind_real),
intent(in) :: lon_loc(n_loc)
43 real(
kind_real),
intent(in) :: lat_loc(n_loc)
44 logical,
intent(in) :: mask_loc(n_loc)
45 real(
kind_real),
intent(in) :: rh_loc(n_loc)
46 integer,
intent(in) :: loc_to_glb(n_loc)
47 integer,
intent(in) :: ntry
48 integer,
intent(in) :: nrep
49 integer,
intent(in) :: ns2_glb
50 integer,
intent(out) :: sam2_glb(ns2_glb)
51 logical,
intent(in),
optional :: fast
52 logical,
intent(in),
optional :: verbosity
53 integer,
intent(in),
optional :: n_uni
54 integer,
intent(in),
optional :: uni_to_loc(:)
55 type(
tree_type),
intent(in),
optional :: tree_uni
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(:)
63 logical :: lfast,lverbosity
64 logical,
allocatable :: mask_glb(:)
65 character(len=1024),
parameter :: subr =
'initialize_sampling'
70 if (
present(fast)) lfast = fast
71 if (
present(verbosity)) lverbosity = verbosity
74 call mpl%f_comm%allreduce(n_loc,n_glb,fckit_mpi_sum())
77 n_loc_eff = count(mask_loc)
78 call mpl%f_comm%allreduce(n_loc_eff,n_glb_eff,fckit_mpi_sum())
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
90 allocate(hash_loc(n_loc))
92 allocate(hash_glb(n_glb))
93 allocate(mask_glb(n_glb))
94 allocate(list(ns2_glb))
95 allocate(order(ns2_glb))
103 hash_loc(i_loc) =
lonlathash(lon_loc(i_loc),lat_loc(i_loc))
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)
114 if (mask_glb(i_glb))
then
116 sam2_glb(is2_glb) = i_glb
117 list(is2_glb) = hash_glb(i_glb)
122 call qsort(ns2_glb,list,order)
125 sam2_glb = sam2_glb(order)
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)
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)
149 & ntry,nrep,ns2_glb,sam2_glb,lfast,lverbosity)
152 deallocate(lon1_glb_eff)
153 deallocate(lat1_glb_eff)
154 deallocate(rh1_glb_eff)
155 deallocate(sam1_glb_eff)
160 call mpl%f_comm%broadcast(sam2_glb,mpl%rootproc-1)
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)
176 integer,
intent(in) :: n_loc
177 integer,
intent(in) :: n_loc_eff
178 integer,
intent(in) :: n_glb_eff
179 real(
kind_real),
intent(in) :: lon_loc(n_loc)
180 real(
kind_real),
intent(in) :: lat_loc(n_loc)
181 logical,
intent(in) :: mask_loc(n_loc)
182 real(
kind_real),
intent(in) :: rh_loc(n_loc)
183 integer,
intent(in) :: loc_to_glb(n_loc)
184 integer,
intent(in) :: ns2_glb
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
191 integer,
intent(in),
optional :: n_uni
192 integer,
intent(in),
optional :: uni_to_loc(:)
193 type(
tree_type),
intent(in),
optional :: tree_uni
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(:)
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
211 if (
present(verbosity)) lverbosity = verbosity
212 lfull_grid =
present(n_uni).and.
present(uni_to_loc).and.
present(tree_uni)
224 ns1_glb = nfac*int(real(ns2_glb,
kind_real)*4.0*
pi/area)
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
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))
243 if (lverbosity)
call mpl%prog_init(ns1_glb)
254 lonlat = agrid%lonlat(ix,iy)*
deg2rad
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))
260 if (mpl%msv%isnot(i_loc))
then
262 if (lmask(i_loc))
then
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.
275 if (lverbosity)
call mpl%prog_print(is1_glb)
278 if (lverbosity)
call mpl%prog_final(.false.)
291 allocate(lon1_loc(ns1_loc))
292 allocate(lat1_loc(ns1_loc))
293 allocate(rh1_loc(ns1_loc))
294 allocate(sam1_loc(ns1_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)
312 write(mpl%info,
'(a)')
' => '
313 if (lverbosity)
call mpl%flush(.false.)
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))
324 lon1_glb = mpl%msv%valr
325 lat1_glb = mpl%msv%valr
326 rh1_glb = mpl%msv%valr
327 sam1_glb = mpl%msv%vali
329 dist1_glb = mpl%msv%valr
336 if (iproc==mpl%rootproc)
then
338 ns1_loc_tmp = ns1_loc
341 call mpl%f_comm%receive(ns1_loc_tmp,iproc-1,mpl%tag,status)
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))
350 allocate(dist1_loc_tmp(ns1_loc_tmp))
351 allocate(s1_loc_to_glb_tmp(ns1_loc_tmp))
354 if (iproc==mpl%rootproc)
then
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)
361 dist1_loc_tmp = dist1_loc(1:ns1_loc)
362 s1_loc_to_glb_tmp = s1_loc_to_glb(1:ns1_loc)
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)
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)
378 do is1_loc=1,ns1_loc_tmp
380 is1_glb = s1_loc_to_glb_tmp(is1_loc)
383 update = mpl%msv%is(dist1_glb(is1_glb))
384 if (.not.update) update = (dist1_loc_tmp(is1_loc)<dist1_glb(is1_glb))
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)
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
403 deallocate(lon1_loc_tmp)
404 deallocate(lat1_loc_tmp)
405 deallocate(rh1_loc_tmp)
406 deallocate(sam1_loc_tmp)
408 deallocate(dist1_loc_tmp)
409 deallocate(s1_loc_to_glb_tmp)
415 ns1_glb_eff = count(mpl%msv%isnot(dist1_glb))
417 ns1_glb_eff = ns1_glb
419 retry = (ns1_glb_eff<ns2_glb)
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)
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)
436 write(mpl%info,
'(a)')
''
437 if (lverbosity)
call mpl%flush
440 call mpl%update_tag(7)
442 call mpl%update_tag(5)
446 call mpl%f_comm%broadcast(retry,mpl%rootproc-1)
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')
458 deallocate(dist1_loc)
459 deallocate(s1_loc_to_glb)
461 if (mpl%main.and.retry)
then
466 if (lfull_grid)
deallocate(dist1_glb)
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))
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)
490 lon1_glb_eff = lon1_glb
491 lat1_glb_eff = lat1_glb
492 rh1_glb_eff = rh1_glb
493 sam1_glb_eff = sam1_glb
501 if (lfull_grid)
deallocate(dist1_glb)
512 & ns2_glb,sam2_glb,fast,verbosity)
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
525 integer,
intent(in) :: nrep
526 integer,
intent(in) :: ns2_glb
527 integer,
intent(out) :: sam2_glb(ns2_glb)
528 logical,
intent(in),
optional :: fast
529 logical,
intent(in),
optional :: verbosity
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'
546 if (
present(fast)) lfast = fast
547 if (
present(verbosity)) lverbosity = verbosity
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))
553 call qsort(ns1_glb_eff,list,order)
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)
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))
569 sam2_glb_tmp = mpl%msv%vali
572 to_valid = mpl%msv%vali
573 do is1_glb_eff=1,ns1_glb_eff
574 to_valid(is1_glb_eff) = is1_glb_eff
576 ns1_glb_val = ns1_glb_eff
577 if (lverbosity)
call mpl%prog_init(ns2_glb+nrep_eff)
583 allocate(cdf(ns1_glb_eff))
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
592 cdf_norm = 1.0/cdf(ns1_glb_eff)
595 do is2_glb=1,ns2_glb+nrep_eff
597 call rng%rand_real(0.0_kind_real,1.0_kind_real,rr)
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
613 sam2_glb_tmp(is2_glb) = ir
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)
620 ns1_glb_val = ns1_glb_val-1
623 cdf_norm = 1.0/cdf(ns1_glb_val)
624 cdf(1:ns1_glb_val) = cdf(1:ns1_glb_val)*cdf_norm
627 if (lverbosity)
call mpl%prog_print(is2_glb)
629 if (lverbosity)
call mpl%prog_final(.false.)
635 do is2_glb=1,ns2_glb+nrep_eff
638 call tree%alloc(mpl,ns1_glb_eff,mask=smask)
641 call tree%init(lon1_glb_tmp,lat1_glb_tmp)
653 call rng%rand_integer(1,ns1_glb_val,irval)
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)), &
669 call tree%find_nearest_neighbors(lon1_glb_tmp(ir),lat1_glb_tmp(ir),1,nn_index(1:1),nn_dist(1:1))
671 d = nn_dist(1)**2/(rh1_glb_tmp(ir)**2+rh1_glb_tmp(nn_index(1))**2)
674 if (
sup(d,distmax))
then
683 if (is2_glb>2)
call tree%dealloc
688 sam2_glb_tmp(is2_glb) = irmax
689 lmask(irmax) = .false.
690 smask(irmax) = .true.
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
698 if (lverbosity)
call mpl%prog_print(is2_glb)
700 if (lverbosity)
call mpl%prog_final(.false.)
705 write(mpl%info,
'(a)')
' => '
706 if (lverbosity)
call mpl%flush(.false.)
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))
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))
721 if (lverbosity)
call mpl%prog_init(nrep_eff)
726 call tree%alloc(mpl,ns2_glb+nrep_eff,mask=rmask)
729 call tree%init(lon_rep,lat_rep)
732 do is2_glb=1,ns2_glb+nrep_eff
733 if (rmask(is2_glb))
then
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)
742 call mpl%abort(subr,
'wrong index in replacement')
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)
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)
763 rmask(is2_glb_min) = .false.
766 if (lverbosity)
call mpl%prog_print(irep)
768 if (lverbosity)
call mpl%prog_final
772 do is2_glb=1,ns2_glb+nrep_eff
773 if (rmask(is2_glb))
then
775 sam2_glb(js) = sam2_glb_tmp(is2_glb)
786 write(mpl%info,
'(a)')
''
787 if (lverbosity)
call mpl%flush
790 sam2_glb = sam2_glb_tmp
794 sam2_glb = sam1_glb_tmp(sam2_glb)
797 deallocate(sam2_glb_tmp)