15 use atlas_metadata_module,
only: atlas_metadata
16 use fckit_configuration_module,
only : fckit_configuration
17 use fckit_log_module,
only: fckit_log
18 use fckit_mpi_module,
only: fckit_mpi_comm, fckit_mpi_sum,fckit_mpi_min,&
19 fckit_mpi_max,fckit_mpi_status
55 type(atlas_functionspace) :: in_funcspace
56 type(atlas_functionspace) :: out_funcspace
60 integer,
public :: nlev
62 integer,
public :: nout
63 integer,
public :: nout_local
93 #define LISTED_TYPE bump_interpolator
96 #include "saber/util/linkedList_i.f"
106 #include "saber/util/linkedList_c.f"
133 subroutine bint_init(self, config, comm, in_funcspace, out_funcspace, masks)
134 use fckit_configuration_module,
only : fckit_configuration
135 use iso_c_binding,
only : c_char
137 class(fckit_configuration),
intent(in) :: config
138 type(fckit_mpi_comm) ,
intent(in) :: comm
139 class(atlas_functionspace),
intent(in) :: in_funcspace
140 class(atlas_functionspace),
intent(in) :: out_funcspace
144 type(fckit_configuration) :: bump_config
146 real(kind_real) :: msvalr
147 integer,
allocatable :: levels(:)
148 character(len=max_string) :: msg
149 character(len=max_string) :: myname =
"saber::interpolation::bump_interpolation_mod::bint_init "
154 call config%get_or_die(
"bump",bump_config)
157 call self%bump%nam%init(comm%size())
160 self%bump%nam%default_seed = .true.
161 self%bump%nam%prefix =
'bump_interpolator'
162 self%bump%nam%datadir =
'testdata'
165 call self%bump%nam%from_conf(bump_config)
170 If (config%has(
"missingvalue_int")) &
171 call config%get_or_die(
"missingvalue_int",msvali)
172 If (config%has(
"missingvalue_real")) &
173 call config%get_or_die(
"missingvalue_real",msvalr)
176 self%in_funcspace = in_funcspace
177 self%out_funcspace = out_funcspace
183 If (config%has(
"nlevels"))
call config%get_or_die(
"nlevels",self%nlev)
185 self%bump%nam%nl = self%nlev
187 allocate(levels(1:self%nlev))
188 if (config%has(
"levels"))
then
189 call config%get_or_die(
"levels",levels)
195 self%bump%nam%levs(1:self%nlev) = levels
206 if (
present(masks))
then
207 call self%bump%setup(self%bump%mpl%f_comm, in_funcspace, masks, &
208 msvali=msvali, msvalr=msvalr)
210 call self%bump%setup(self%bump%mpl%f_comm, in_funcspace, &
211 msvali=msvali, msvalr=msvalr)
217 call fckit_log%info(
'-------------------------------------------------------------------')
218 call fckit_log%info(
'--- Initialize output grid')
220 self%outgeom%nl0 = self%nlev
223 call self%outgeom%from_atlas(self%bump%mpl, out_funcspace)
224 self%nout_local =
size(self%outgeom%lon_mga)
230 call self%bump%run_drivers
233 call fckit_log%info(
'-------------------------------------------------------------------')
234 call fckit_log%info(
'--- Run bump_interpolation driver')
235 call self%driver(self%bump%mpl,self%bump%rng,self%bump%nam,self%bump%geom)
236 if (self%bump%nam%default_seed)
call self%bump%rng%reseed(self%bump%mpl)
241 if (self%nout < 1)
then
242 write(msg,
'(a)') trim(myname) //
": " //
"Output grid not properly defined"
246 if (self%nlev /= self%bump%geom%nl0)
then
247 write(msg,
'(a)') trim(myname) //
": " // &
248 "number of levels does not match bump geom"
252 if (self%nlev /= self%outgeom%nl0)
then
253 write(msg,
'(a)') trim(myname) //
": " // &
254 "ERROR: Output grid has different number of levels than input grid!"
275 integer :: iouta,iproc,i_s,ic0,ic0u,jc0u,ic0b,ic0a,nouta_eff
276 integer :: nout_eff,nn_index(1),proc_to_nouta(mpl%nproc),proc_to_nouta_eff(mpl%nproc)
277 integer :: c0u_to_c0b(geom%nc0u)
278 integer,
allocatable :: c0b_to_c0(:)
279 real(kind_real) :: nn_dist(1),N_max,C_max
280 logical :: maskouta(self%nout_local),lcheck_nc0b(geom%nc0)
281 character(len=1024),
parameter :: subr =
'bint_driver'
284 if (any(.not.geom%myuniverse))
call mpl%abort(subr,
'universe should be global for interpolation')
287 if (self%nout_local > 0)
then
288 do iouta=1,self%nout_local
289 call geom%mesh_c0u%inside(mpl,self%outgeom%lon_mga(iouta),self%outgeom%lat_mga(iouta),maskouta(iouta))
290 if (.not.maskouta(iouta))
then
292 call geom%tree_c0u%find_nearest_neighbors(self%outgeom%lon_mga(iouta),self%outgeom%lat_mga(iouta), &
293 & 1,nn_index,nn_dist)
294 if (nn_dist(1)<
rth) maskouta(iouta) = .true.
297 nouta_eff = count(maskouta)
303 call mpl%f_comm%allgather(self%nout_local,proc_to_nouta)
304 call mpl%f_comm%allgather(nouta_eff,proc_to_nouta_eff)
305 self%nout = sum(proc_to_nouta)
306 nout_eff = sum(proc_to_nouta_eff)
309 write(mpl%info,
'(a7,a)')
'',
'Number of points in output grid / valid points per MPI task:'
312 write(mpl%info,
'(a10,a,i3,a,i8,a,i8)')
'',
'Task ',iproc,
': ', &
313 proc_to_nouta(iproc),
' / ',proc_to_nouta_eff(iproc)
316 write(mpl%info,
'(a10,a,i8,a,i8)')
'',
'Total : ',self%nout,
' / ',nout_eff
321 write(mpl%info,
'(a7,a)')
'',
'Single level:'
323 call self%h%interp(mpl,rng,nam,geom,0,geom%nc0u,geom%lon_c0u,geom%lat_c0u,&
324 geom%gmask_hor_c0u,self%nout_local,self%outgeom%lon_mga,&
325 self%outgeom%lat_mga,maskouta,10)
328 lcheck_nc0b = .false.
330 ic0u = geom%c0a_to_c0u(ic0a)
331 lcheck_nc0b(ic0u) = .true.
333 do iouta=1,self%nout_local
335 jc0u = self%h%col(i_s)
336 lcheck_nc0b(jc0u) = .true.
339 self%nc0b = count(lcheck_nc0b)
342 allocate(c0b_to_c0(self%nc0b))
345 c0u_to_c0b = mpl%msv%vali
348 if (lcheck_nc0b(ic0u))
then
350 ic0 = geom%c0u_to_c0(ic0u)
351 c0b_to_c0(ic0b) = ic0
352 c0u_to_c0b(ic0u) = ic0b
357 self%h%n_src = self%nc0b
359 self%h%col(i_s) = c0u_to_c0b(self%h%col(i_s))
363 call self%com%setup(mpl,
'com',geom%nc0a,self%nc0b,geom%nc0,geom%c0a_to_c0,c0b_to_c0)
366 if (nout_eff > 0)
then
367 call mpl%f_comm%allreduce(real(self%com%nhalo,kind_real),c_max,fckit_mpi_max())
368 c_max = c_max/(3.0*real(nout_eff,kind_real)/real(mpl%nproc,kind_real))
369 n_max = real(maxval(proc_to_nouta_eff),kind_real)/(real(nout_eff,kind_real)/real(mpl%nproc,kind_real))
372 write(mpl%info,
'(a7,a,f5.1,a)')
'',
'Output grid repartition imbalance: ',100.0*real(maxval(proc_to_nouta_eff) &
373 & -minval(proc_to_nouta_eff),kind_real)/(real(sum(proc_to_nouta_eff),kind_real)/real(mpl%nproc,kind_real)),
' %'
375 write(mpl%info,
'(a7,a,i8,a,i8,a,i8)')
'',
'Number of grid points / halo size / number of received values: ', &
376 & self%com%nred,
' / ',self%com%next,
' / ',self%com%nhalo
378 write(mpl%info,
'(a7,a,f10.2,a,f10.2)')
'',
'Scores (N_max / C_max):',n_max,
' / ',c_max
383 deallocate(c0b_to_c0)
405 type(atlas_field) :: infield, outfield
406 real(kind_real),
allocatable :: infld_mga(:,:), infld_c0a(:,:)
407 real(kind_real),
allocatable :: outfld(:,:)
409 character(len=max_string) :: fieldname
412 allocate(infld_mga(self%bump%geom%nmga,self%nlev))
413 allocate(outfld(self%nout_local,self%nlev))
415 if (.not.self%bump%geom%same_grid)
then
416 allocate(infld_c0a(self%bump%geom%nc0a, self%nlev))
419 do ifield = 1, infields%size()
421 infield = infields%field(ifield)
422 fieldname = infield%name()
427 if (.not. outfields%has_field(fieldname))
then
428 outfield = self%out_funcspace%create_field(name=fieldname, &
429 kind=atlas_real(kind_real),levels=self%nlev)
431 outfield = outfields%field(name=fieldname)
440 if (self%bump%geom%same_grid)
then
441 call self%apply_interp(infld_mga,outfld)
444 infld_c0a = 0.0_kind_real
445 call self%bump%geom%copy_mga_to_c0a(self%bump%mpl,infld_mga,infld_c0a)
448 call self%apply_interp(infld_c0a,outfld)
456 if (.not. outfields%has_field(fieldname))
then
457 call outfields%add(outfield)
462 call outfield%final()
467 if (
allocated(infld_mga))
deallocate(infld_mga)
468 if (
allocated(infld_c0a))
deallocate(infld_c0a)
469 if (
allocated(outfld))
deallocate(outfld)
482 real(kind_real) ,
intent(in) :: infield(self%bump%geom%nc0a,self%nlev)
483 real(kind_real) ,
intent(out) :: outfield(self%nout_local,self%nlev)
486 real(kind_real),
allocatable :: infield_ext(:,:)
488 allocate(infield_ext(self%nc0b,self%nlev))
491 call self%com%ext(self%bump%mpl, self%nlev, infield ,infield_ext)
493 if (self%nout_local > 0)
then
496 do ilev = 1, self%nlev
497 call self%h%apply(self%bump%mpl,infield_ext(:,ilev),outfield(:,ilev))
502 deallocate(infield_ext)
529 type(atlas_field) :: field_ingrid, field_outgrid
530 real(kind_real),
allocatable :: fld_ingrid_mga(:,:), fld_ingrid_c0a(:,:)
531 real(kind_real),
allocatable :: fld_outgrid(:,:)
533 character(len=max_string) :: fieldname
536 allocate(fld_ingrid_mga(self%bump%geom%nmga,self%nlev))
537 allocate(fld_outgrid(self%nout_local,self%nlev))
539 if (.not.self%bump%geom%same_grid)
then
540 allocate(fld_ingrid_c0a(self%bump%geom%nc0a, self%nlev))
543 do ifield = 1, fields_outgrid%size()
545 field_outgrid = fields_outgrid%field(ifield)
546 fieldname = field_outgrid%name()
551 if (.not. fields_ingrid%has_field(fieldname))
then
552 field_ingrid = self%in_funcspace%create_field(name=fieldname, &
553 kind=atlas_real(kind_real),levels=self%nlev)
555 field_ingrid = fields_ingrid%field(name=fieldname)
564 if (self%bump%geom%same_grid)
then
566 call self%apply_interp_ad(fld_outgrid,fld_ingrid_mga)
569 call self%apply_interp_ad(fld_outgrid,fld_ingrid_c0a)
572 call self%bump%geom%copy_c0a_to_mga(self%bump%mpl,fld_ingrid_c0a,fld_ingrid_mga)
581 if (.not. fields_ingrid%has_field(fieldname))
then
582 call fields_ingrid%add(field_ingrid)
586 call field_ingrid%final()
587 call field_outgrid%final()
592 if (
allocated(fld_ingrid_mga))
deallocate(fld_ingrid_mga)
593 if (
allocated(fld_ingrid_c0a))
deallocate(fld_ingrid_c0a)
594 if (
allocated(fld_outgrid))
deallocate(fld_outgrid)
611 real(kind_real),
intent(in) :: fld_outgrid(self%nout_local,self%nlev)
612 real(kind_real),
intent(out) :: fld_ingrid(self%bump%geom%nc0a,self%nlev)
615 real(kind_real) :: fld_ingrid_ext(self%nc0b,self%nlev)
617 if (self%nout_local > 0)
then
620 do ilev = 1, self%nlev
621 call self%h%apply_ad(self%bump%mpl,fld_outgrid(:,ilev),fld_ingrid_ext(:,ilev))
630 call self%com%red(self%bump%mpl,self%nlev,fld_ingrid_ext,fld_ingrid)
641 call self%outgeom%dealloc()
652 call self%bump%dealloc()
654 call self%deallocate_outgrid()
655 call self%h%dealloc()
656 call self%com%dealloc()