Line data Source code
1 : module interp_mod
2 : use shr_kind_mod, only: r8 => shr_kind_r8, r4 => shr_kind_r4
3 : use dimensions_mod, only: nelemd, np, ne
4 : use interpolate_mod, only: interpdata_t
5 : use interpolate_mod, only: interp_lat => lat, interp_lon => lon
6 : use interpolate_mod, only: interp_gweight => gweight
7 : use dyn_grid, only: elem,fvm
8 : use spmd_utils, only: iam
9 : use cam_history_support, only: fillvalue
10 : use hybrid_mod, only: hybrid_t, config_thread_region
11 : use cam_abortutils, only: endrun
12 :
13 : implicit none
14 : private
15 : save
16 :
17 : public :: setup_history_interpolation
18 : public :: set_interp_hfile
19 : public :: write_interpolated
20 :
21 : interface write_interpolated
22 : module procedure write_interpolated_scalar
23 : module procedure write_interpolated_vector
24 : end interface
25 :
26 : ! hybrid is created in setup_history_interpolation
27 : type(hybrid_t) :: hybrid
28 :
29 : ! A type to hold interpdata info for each interpolated history file
30 : type cam_interpolate_t
31 : type(interpdata_t), pointer :: interpdata(:) => NULL()
32 : end type cam_interpolate_t
33 :
34 : type(cam_interpolate_t), pointer :: interpdata_set(:) => NULL() ! all files
35 : type(interpdata_t), pointer :: cam_interpolate(:) => NULL() ! curr. file
36 :
37 : CONTAINS
38 :
39 0 : subroutine setup_history_interpolation(interp_ok, mtapes, interp_output, &
40 0 : interp_info)
41 :
42 : use cam_history_support, only: interp_info_t
43 : use cam_history_support, only: interp_type_native
44 : use cam_history_support, only: interp_type_bilinear
45 : use cam_history_support, only: interp_gridtype_equal_poles
46 : use cam_history_support, only: interp_gridtype_gauss
47 : use cam_history_support, only: interp_gridtype_equal_nopoles
48 : use cam_grid_support, only: horiz_coord_t, horiz_coord_create, iMap
49 : use cam_grid_support, only: cam_grid_register, cam_grid_attribute_register
50 : use cam_grid_support, only: max_hcoordname_len
51 : use interpolate_mod, only: get_interp_lat, get_interp_lon
52 : use interpolate_mod, only: get_interp_parameter, set_interp_parameter
53 : use interpolate_mod, only: get_interp_gweight, setup_latlon_interp
54 : use parallel_mod, only: par
55 :
56 : ! Dummy arguments
57 : logical, intent(inout) :: interp_ok
58 : integer, intent(in) :: mtapes
59 : logical, intent(in) :: interp_output(:)
60 : type(interp_info_t), intent(inout) :: interp_info(:)
61 :
62 : ! Local variables
63 : integer :: i
64 0 : real(r8), pointer :: w(:)
65 0 : integer(iMap), pointer :: grid_map(:,:)
66 : type(horiz_coord_t), pointer :: lat_coord
67 : type(horiz_coord_t), pointer :: lon_coord
68 : character(len=max_hcoordname_len) :: gridname
69 : !---------------------------------------------------------------------------
70 :
71 0 : nullify(grid_map)
72 :
73 : ! For this dycore, interpolated output should be OK
74 0 : interp_ok = (iam < par%nprocs)
75 :
76 0 : if (interp_ok) then
77 0 : hybrid = config_thread_region(par,'serial')
78 :
79 0 : if(any(interp_output)) then
80 0 : allocate(interpdata_set(mtapes))
81 0 : do i = 1, mtapes
82 0 : if (interp_output(i)) then
83 0 : if ( (interp_info(i)%interp_nlon == 0) .or. &
84 : (interp_info(i)%interp_nlat == 0)) then
85 : ! compute interpolation grid based on number of points around equator
86 0 : call set_interp_parameter('auto', (4 * ne * (np-1)))
87 0 : interp_info(i)%interp_nlat = get_interp_parameter('nlat')
88 0 : interp_info(i)%interp_nlon = get_interp_parameter('nlon')
89 : else
90 : call set_interp_parameter('nlat', interp_info(i)%interp_nlat)
91 0 : call set_interp_parameter('nlon', interp_info(i)%interp_nlon)
92 : end if
93 0 : call set_interp_parameter('itype', interp_info(i)%interp_type)
94 0 : call set_interp_parameter('gridtype', interp_info(i)%interp_gridtype)
95 :
96 0 : allocate(interpdata_set(i)%interpdata(nelemd))
97 : ! Reset pointers in the interpolate module so they are not
98 : ! overwritten
99 0 : nullify(interp_lat)
100 0 : nullify(interp_lon)
101 0 : nullify(interp_gweight)
102 : call setup_latlon_interp(elem, interpdata_set(i)%interpdata, par)
103 : ! Create the grid coordinates
104 : lat_coord => horiz_coord_create('lat', '', &
105 0 : interp_info(i)%interp_nlat, 'latitude', 'degrees_north', &
106 0 : 1, interp_info(i)%interp_nlat, get_interp_lat())
107 : lon_coord => horiz_coord_create('lon', '', &
108 0 : interp_info(i)%interp_nlon, 'longitude', 'degrees_east', &
109 0 : 1, interp_info(i)%interp_nlon, get_interp_lon())
110 : ! Create a grid for this history file
111 0 : write(gridname, '(a,i0)') 'interp_out_', i
112 0 : interp_info(i)%grid_id = 200 + i
113 0 : call cam_grid_register(trim(gridname), interp_info(i)%grid_id, &
114 0 : lat_coord, lon_coord, grid_map, unstruct=.false.)
115 0 : interp_info(i)%gridname = trim(gridname)
116 : ! Add grid attributes
117 0 : allocate(w(get_interp_parameter('nlat')))
118 0 : w = get_interp_gweight()
119 0 : select case(interp_info(i)%interp_gridtype)
120 : case(interp_gridtype_equal_poles)
121 : call cam_grid_attribute_register(trim(gridname), &
122 0 : 'interp_outputgridtype', 'equally spaced with poles')
123 : call cam_grid_attribute_register(trim(gridname), 'w', &
124 0 : 'latitude weights', 'lat', w)
125 : case(interp_gridtype_equal_nopoles)
126 : call cam_grid_attribute_register(trim(gridname), &
127 0 : 'interp_outputgridtype', 'equally spaced no poles')
128 : call cam_grid_attribute_register(trim(gridname), 'gw', &
129 0 : 'latitude weights', 'lat', w)
130 : case(interp_gridtype_gauss)
131 : call cam_grid_attribute_register(trim(gridname), &
132 0 : 'interp_outputgridtype', 'Gauss')
133 : call cam_grid_attribute_register(trim(gridname), 'gw', &
134 0 : 'gauss weights', 'lat', w)
135 : case default
136 : call cam_grid_attribute_register(trim(gridname), &
137 : 'interp_outputgridtype', &
138 : 'Unknown interpolation output grid type', &
139 0 : interp_info(i)%interp_gridtype)
140 : end select
141 0 : nullify(w) ! belongs to attribute
142 0 : if(interp_info(i)%interp_type == interp_type_native) then
143 : call cam_grid_attribute_register(trim(gridname), &
144 0 : 'interp_type', 'se basis functions')
145 0 : else if(interp_info(i)%interp_type == interp_type_bilinear) then
146 : call cam_grid_attribute_register(trim(gridname), &
147 0 : 'interp_type', 'bilinear')
148 : else
149 : call cam_grid_attribute_register(trim(gridname), 'interp_type', &
150 0 : 'Unknown interpolation type', interp_info(i)%interp_type)
151 : end if
152 : ! Store the data pointers for reuse later
153 0 : interp_info(i)%interp_lat => interp_lat
154 0 : interp_info(i)%interp_lon => interp_lon
155 0 : interp_info(i)%interp_gweight => interp_gweight
156 : end if
157 : end do
158 : end if
159 : end if
160 :
161 0 : end subroutine setup_history_interpolation
162 :
163 0 : subroutine set_interp_hfile(hfilenum, interp_info)
164 0 : use cam_history_support, only: interp_info_t
165 : use interpolate_mod, only: set_interp_parameter
166 :
167 : ! Dummy arguments
168 : integer, intent(in) :: hfilenum
169 : type(interp_info_t), intent(inout) :: interp_info(:)
170 :
171 0 : if (.not. associated(interpdata_set)) then
172 0 : call endrun('SET_INTERP_HFILE: interpdata_set not allocated')
173 0 : else if ((hfilenum < 1) .or. (hfilenum > size(interpdata_set))) then
174 0 : call endrun('SET_INTERP_HFILE: hfilenum out of range')
175 0 : else if (hfilenum > size(interp_info)) then
176 0 : call endrun('SET_INTERP_HFILE: hfilenum out of range')
177 : else
178 0 : cam_interpolate => interpdata_set(hfilenum)%interpdata
179 0 : interp_lat => interp_info(hfilenum)%interp_lat
180 0 : interp_lon => interp_info(hfilenum)%interp_lon
181 0 : interp_gweight => interp_info(hfilenum)%interp_gweight
182 : call set_interp_parameter('nlat', interp_info(hfilenum)%interp_nlat)
183 : call set_interp_parameter('nlon', interp_info(hfilenum)%interp_nlon)
184 : call set_interp_parameter('itype', interp_info(hfilenum)%interp_type)
185 : call set_interp_parameter('gridtype', interp_info(hfilenum)%interp_gridtype)
186 : end if
187 0 : end subroutine set_interp_hfile
188 :
189 0 : subroutine write_interpolated_scalar(File, varid, fld, numlev, data_type, decomp_type)
190 0 : use pio, only: file_desc_t, var_desc_t
191 : use pio, only: iosystem_desc_t
192 : use pio, only: pio_initdecomp, pio_freedecomp
193 : use pio, only: io_desc_t, pio_write_darray
194 : use pio, only: pio_real
195 : use interpolate_mod, only: interpolate_scalar
196 : use cam_instance, only: atm_id
197 : use spmd_dyn, only: local_dp_map
198 : use ppgrid, only: begchunk
199 : use phys_grid, only: get_dyn_col_p, columns_on_task, get_chunk_info_p
200 : use dimensions_mod, only: fv_nphys, nc, nhc, nhc_phys
201 : use dof_mod, only: PutUniquePoints
202 : use interpolate_mod, only: get_interp_parameter
203 : use shr_pio_mod, only: shr_pio_getiosys
204 : use edge_mod, only: edgevpack, edgevunpack, initedgebuffer, freeedgebuffer
205 : use edgetype_mod, only: EdgeBuffer_t
206 : use bndry_mod, only: bndry_exchange
207 : use parallel_mod, only: par
208 : use thread_mod, only: horz_num_threads
209 : use cam_grid_support, only: cam_grid_id
210 : use hybrid_mod, only: hybrid_t, config_thread_region, get_loop_ranges
211 : use fvm_mod, only: fill_halo_and_extend_panel
212 :
213 : type(file_desc_t), intent(inout) :: File
214 : type(var_desc_t) , intent(inout) :: varid
215 : real(r8), intent(in) :: fld(:,:,:)
216 : integer, intent(in) :: numlev, data_type, decomp_type
217 : !
218 : ! local variables
219 : !
220 : type(io_desc_t) :: iodesc
221 : type(hybrid_t) :: hybrid
222 : type(iosystem_desc_t), pointer :: pio_subsystem
223 0 : type (EdgeBuffer_t) :: edgebuf ! edge buffer
224 :
225 :
226 : integer :: lchnk, i, col_index, icol, ncols, ierr
227 : integer :: nets, nete
228 : integer :: phys_decomp, fvm_decomp, gll_decomp
229 :
230 0 : real(r8), pointer :: dest(:,:,:,:)
231 0 : real(r8), pointer :: fldout(:,:)
232 0 : real(r8), allocatable :: fld_dyn(:,:,:), fld_tmp(:,:,:,:,:)
233 :
234 : integer :: st, en ! Start and end temporaries
235 : integer :: ie, blk_ind(1), ncnt_out, k
236 0 : integer, pointer :: idof(:)
237 : integer :: nlon, nlat, ncol, nsize, nhalo, nhcc
238 : logical :: usefillvalues
239 : character(len=*), parameter :: subname = 'write_interpolated_scalar'
240 :
241 0 : usefillvalues=.false.
242 :
243 0 : phys_decomp = cam_grid_id('physgrid')
244 0 : fvm_decomp = cam_grid_id('FVM')
245 0 : gll_decomp = cam_grid_id('GLL')
246 : !
247 : ! There are 2 main scenarios regarding decomposition:
248 : !
249 : ! decomp_type==phys_decomp: we need to move data from physics decomposition to dynamics decomposition
250 : ! else : data is on dynamics decomposition
251 : !
252 0 : if (decomp_type==phys_decomp) then
253 0 : if(.not. local_dp_map) then
254 0 : call endrun(subname//': weak scaling does not support load balancing')
255 : end if
256 0 : if (fv_nphys > 0) then
257 : !
258 : ! note that even if fv_nphys<4 then SIZE(fld,DIM=1)=PCOLS
259 : !
260 0 : nsize = fv_nphys
261 0 : nhalo = 1!for bilinear only a halo of 1 is needed
262 0 : nhcc = nhc_phys
263 : else
264 0 : nsize = np
265 0 : nhalo = 0!no halo needed (lat-lon point always surrounded by GLL points)
266 0 : nhcc = 0
267 : end if
268 0 : else if (decomp_type == fvm_decomp) then
269 : !
270 : ! CSLAM grid output
271 : !
272 0 : nsize = nc
273 0 : nhalo = 1!for bilinear only a halo of 1 is needed
274 0 : nhcc = nhc
275 0 : else if (decomp_type == gll_decomp) then
276 0 : nsize = np
277 0 : nhalo = 0!no halo needed (lat-lon point always surrounded by GLL points)
278 0 : nhcc = 0
279 : else
280 0 : call endrun(subname//': unknown decomp_type')
281 : end if
282 0 : allocate(fld_tmp(1-nhcc:nsize+nhcc,1-nhcc:nsize+nhcc,numlev,1,nelemd))
283 0 : allocate(dest(1-nhalo:nsize+nhalo,1-nhalo:nsize+nhalo,numlev,nelemd))
284 :
285 0 : nlon=get_interp_parameter('nlon')
286 0 : nlat=get_interp_parameter('nlat')
287 0 : pio_subsystem => shr_pio_getiosys(atm_id)
288 :
289 0 : if(decomp_type == phys_decomp) then
290 :
291 0 : allocate(fld_dyn(nsize*nsize,numlev,nelemd))
292 0 : fld_dyn = -999_R8
293 : !!$omp parallel do num_threads(horz_num_threads) private (col_index, lchnk, icol, ie, blk_ind, k)
294 0 : do col_index = 1, columns_on_task
295 0 : call get_dyn_col_p(col_index, ie, blk_ind)
296 0 : call get_chunk_info_p(col_index, lchnk, icol)
297 0 : do k = 1, numlev
298 0 : fld_dyn(blk_ind(1), k, ie) = fld(icol, k, lchnk-begchunk+1)
299 : end do
300 : end do
301 :
302 0 : if (fv_nphys > 0) then
303 0 : do ie = 1, nelemd
304 0 : fld_tmp(1:nsize,1:nsize,:,1,ie) = RESHAPE(fld_dyn(:,:,ie),(/nsize,nsize,numlev/))
305 : end do
306 : else
307 0 : call initEdgeBuffer(par, edgebuf, elem, numlev,nthreads=1)
308 :
309 0 : do ie = 1, nelemd
310 0 : ncols = elem(ie)%idxp%NumUniquePts
311 0 : call putUniquePoints(elem(ie)%idxP, numlev, fld_dyn(1:ncols,1:numlev,ie), fld_tmp(:,:,1:numlev,1,ie))
312 0 : call edgeVpack(edgebuf, fld_tmp(:,:,1:numlev,1,ie), numlev, 0, ie)
313 : end do
314 0 : if(iam < par%nprocs) then
315 0 : call bndry_exchange(par, edgebuf,location=subname)
316 : end if
317 0 : do ie = 1, nelemd
318 0 : call edgeVunpack(edgebuf, fld_tmp(:,:,1:numlev,1,ie), numlev, 0, ie)
319 : end do
320 0 : call freeEdgeBuffer(edgebuf)
321 : !check if fill values are present:
322 0 : usefillvalues = any(fld_tmp == fillvalue)
323 : end if
324 0 : deallocate(fld_dyn)
325 : else
326 : !
327 : ! not physics decomposition
328 : !
329 0 : do ie = 1, nelemd
330 0 : fld_tmp(1:nsize,1:nsize,1:numlev,1,ie) = RESHAPE(fld(1:nsize*nsize,1:numlev,ie),(/nsize,nsize,numlev/))
331 : end do
332 : !check if fillvalues are present:
333 0 : usefillvalues = any(fld_tmp == fillvalue)
334 : end if
335 : !
336 : ! code for non-GLL grids: need to fill halo and interpolate (if on panel edge/corner) for bilinear interpolation
337 : !
338 0 : if (decomp_type==fvm_decomp.or.(fv_nphys>0.and.decomp_type==phys_decomp)) then
339 : !JMD $OMP PARALLEL NUM_THREADS(horz_num_threads), DEFAULT(SHARED), PRIVATE(hybrid,nets,nete,n)
340 : !JMD hybrid = config_thread_region(par,'horizontal')
341 0 : hybrid = config_thread_region(par,'serial')
342 0 : call get_loop_ranges(hybrid,ibeg=nets,iend=nete)
343 0 : call fill_halo_and_extend_panel(elem(nets:nete),fvm(nets:nete),&
344 0 : fld_tmp(:,:,:,:,nets:nete),hybrid,nets,nete,nsize,nhcc,nhalo,numlev,1,.true.,.true.)
345 :
346 : !check if fill values are present:
347 0 : usefillvalues = any(fld_tmp(:,:,:,:,nets:nete) == fillvalue)
348 : end if
349 : !
350 : ! WARNING - 1:nelemd and nets:nete
351 : !
352 : !!$OMP MASTER !JMD
353 0 : dest(:,:,:,1:nelemd) = fld_tmp(1-nhalo:nsize+nhalo,1-nhalo:nsize+nhalo,:,1,1:nelemd)
354 : !!$OMP END MASTER
355 0 : deallocate(fld_tmp)
356 : !
357 : !***************************************************************************
358 : !
359 : ! now data is on dynamics decomposition
360 : !
361 : !***************************************************************************
362 : !
363 0 : ncnt_out = sum(cam_interpolate(1:nelemd)%n_interp)
364 0 : allocate(fldout(ncnt_out,numlev))
365 0 : allocate(idof(ncnt_out*numlev))
366 0 : fldout = -999_r8
367 0 : idof = 0
368 0 : st = 1
369 :
370 0 : do ie=1,nelemd
371 0 : ncol = cam_interpolate(ie)%n_interp
372 0 : do k=0,numlev-1
373 0 : do i=1,ncol
374 0 : idof(st+i-1+k*ncnt_out)=cam_interpolate(ie)%ilon(i)+nlon*(cam_interpolate(ie)%ilat(i)-1)+nlon*nlat*k
375 : enddo
376 : enddo
377 : ! Now that we have the field on the dyn grid we need to interpolate
378 0 : en = st+cam_interpolate(ie)%n_interp-1
379 0 : if(usefillvalues) then
380 0 : call interpolate_scalar(cam_interpolate(ie),dest(:,:,:,ie), nsize, nhalo, numlev, fldout(st:en,:), fillvalue)
381 : else
382 0 : call interpolate_scalar(cam_interpolate(ie),dest(:,:,:,ie), nsize, nhalo, numlev, fldout(st:en,:))
383 : end if
384 0 : st = en+1
385 : end do
386 :
387 0 : if(numlev==1) then
388 0 : call pio_initdecomp(pio_subsystem, data_type, (/nlon,nlat/), idof, iodesc)
389 : else
390 0 : call pio_initdecomp(pio_subsystem, data_type, (/nlon,nlat,numlev/), idof, iodesc)
391 : end if
392 :
393 0 : if(data_type == pio_real) then
394 0 : call pio_write_darray(File, varid, iodesc, real(fldout, r4), ierr)
395 : else
396 0 : call pio_write_darray(File, varid, iodesc, fldout, ierr)
397 : end if
398 :
399 0 : deallocate(dest)
400 :
401 0 : deallocate(fldout)
402 0 : deallocate(idof)
403 0 : call pio_freedecomp(file,iodesc)
404 :
405 0 : end subroutine write_interpolated_scalar
406 :
407 0 : subroutine write_interpolated_vector(File, varidu, varidv, fldu, fldv, numlev, data_type, decomp_type)
408 0 : use pio, only: file_desc_t, var_desc_t
409 : use pio, only: iosystem_desc_t
410 : use pio, only: pio_initdecomp, pio_freedecomp
411 : use pio, only: io_desc_t, pio_write_darray
412 : use pio, only: pio_real
413 : use cam_instance, only: atm_id
414 : use interpolate_mod, only: interpolate_scalar, vec_latlon_to_contra,get_interp_parameter
415 : use spmd_dyn, only: local_dp_map
416 : use ppgrid, only: begchunk
417 : use phys_grid, only: get_dyn_col_p, columns_on_task, get_chunk_info_p
418 : use dimensions_mod, only: fv_nphys,nc,nhc,nhc_phys
419 : use dof_mod, only: PutUniquePoints
420 : use shr_pio_mod, only: shr_pio_getiosys
421 : use edge_mod, only: edgevpack, edgevunpack, initedgebuffer, freeedgebuffer
422 : use edgetype_mod, only: EdgeBuffer_t
423 : use bndry_mod, only: bndry_exchange
424 : use parallel_mod, only: par
425 : use thread_mod, only: horz_num_threads
426 : use cam_grid_support, only: cam_grid_id
427 : use hybrid_mod, only: hybrid_t,config_thread_region, get_loop_ranges
428 : use fvm_mod, only: fill_halo_and_extend_panel
429 : use control_mod, only: cubed_sphere_map
430 : use cube_mod, only: dmap
431 :
432 : type(file_desc_t), intent(inout) :: File
433 : type(var_desc_t), intent(inout) :: varidu, varidv
434 : real(r8), intent(in) :: fldu(:,:,:), fldv(:,:,:)
435 : integer, intent(in) :: numlev, data_type, decomp_type
436 :
437 : type(hybrid_t) :: hybrid
438 : type(io_desc_t) :: iodesc
439 : type(iosystem_desc_t), pointer :: pio_subsystem
440 0 : type (EdgeBuffer_t) :: edgebuf ! edge buffer
441 :
442 : integer :: lchnk, i, col_index, icol, ncols, ierr
443 : integer :: nets, nete
444 :
445 0 : real(r8), allocatable :: dest(:,:,:,:,:)
446 0 : real(r8), pointer :: fldout(:,:,:)
447 0 : real(r8), allocatable :: fld_dyn(:,:,:,:), fld_tmp(:,:,:,:,:)
448 :
449 : integer :: st, en ! Start and end temporaries
450 : integer :: ie, blk_ind(1), ncnt_out, k
451 0 : integer, pointer :: idof(:)
452 : integer :: nlon, nlat, ncol,nsize,nhalo,nhcc
453 : logical :: usefillvalues
454 : integer :: phys_decomp, fvm_decomp,gll_decomp
455 : real (r8) :: D(2,2) ! derivative of gnomonic mapping
456 : real (r8) :: v1,v2
457 : character(len=*), parameter :: subname = 'write_interpolated_vector'
458 :
459 0 : usefillvalues=.false.
460 :
461 0 : phys_decomp = cam_grid_id('physgrid')
462 0 : fvm_decomp = cam_grid_id('FVM')
463 0 : gll_decomp = cam_grid_id('GLL')
464 : !
465 : ! There are 2 main scenarios regarding decomposition:
466 : !
467 : ! decomp_type==phys_decomp: we need to move data from physics decomposition to dynamics decomposition
468 : ! else : data is on dynamics decomposition
469 : !
470 0 : if (decomp_type==phys_decomp) then
471 0 : if(.not. local_dp_map) then
472 0 : call endrun(subname//': weak scaling does not support load balancing')
473 : end if
474 0 : if (fv_nphys > 0) then
475 : !
476 : ! note that even if fv_nphys<4 then SIZE(fld,DIM=1)=npsq
477 : !
478 0 : nsize = fv_nphys
479 0 : nhalo = 1!for bilinear only a halo of 1 is needed
480 0 : nhcc = nhc_phys
481 : else
482 0 : nsize = np
483 0 : nhalo = 0!no halo needed (lat-lon point always surrounded by GLL points)
484 0 : nhcc = 0
485 : end if
486 0 : else if (decomp_type == fvm_decomp) then
487 : !
488 : ! CSLAM grid output
489 : !
490 0 : nsize = nc
491 0 : nhalo = 1!for bilinear only a halo of 1 is needed
492 0 : nhcc = nhc
493 0 : else if (decomp_type == gll_decomp) then
494 0 : nsize = np
495 0 : nhalo = 0!no halo needed (lat-lon point always surrounded by GLL points)
496 0 : nhcc = 0
497 : else
498 0 : call endrun(subname//': unknown decomp_type')
499 : end if
500 0 : allocate(fld_tmp(1-nhcc:nsize+nhcc,1-nhcc:nsize+nhcc,2,numlev,nelemd))
501 0 : allocate(dest(1-nhalo:nsize+nhalo,1-nhalo:nsize+nhalo,2,numlev,nelemd))
502 :
503 0 : nlon=get_interp_parameter('nlon')
504 0 : nlat=get_interp_parameter('nlat')
505 0 : pio_subsystem => shr_pio_getiosys(atm_id)
506 0 : if(decomp_type == phys_decomp) then
507 0 : allocate(fld_dyn(nsize*nsize,2,numlev,nelemd))
508 0 : fld_dyn = -999_R8
509 : !!$omp parallel do num_threads(horz_num_threads) private (col_index, lchnk, icol, ie, blk_ind, k)
510 0 : do col_index = 1, columns_on_task
511 0 : call get_dyn_col_p(col_index, ie, blk_ind)
512 0 : call get_chunk_info_p(col_index, lchnk, icol)
513 0 : do k = 1, numlev
514 0 : fld_dyn(blk_ind(1), 1, k, ie) = fldu(icol, k, lchnk-begchunk+1)
515 0 : fld_dyn(blk_ind(1), 2, k, ie) = fldv(icol, k, lchnk-begchunk+1)
516 : end do
517 : end do
518 0 : if (fv_nphys > 0) then
519 0 : do ie = 1, nelemd
520 0 : fld_tmp(1:nsize,1:nsize,:,:,ie) = RESHAPE(fld_dyn(:,:,:,ie),(/nsize,nsize,2,numlev/))
521 : end do
522 : else
523 0 : call initEdgeBuffer(par, edgebuf, elem, 2*numlev,nthreads=1)
524 :
525 0 : do ie = 1, nelemd
526 0 : ncols = elem(ie)%idxp%NumUniquePts
527 0 : call putUniquePoints(elem(ie)%idxP, 2, numlev, fld_dyn(1:ncols,:,1:numlev,ie), fld_tmp(:,:,:,1:numlev,ie))
528 0 : call edgeVpack(edgebuf, fld_tmp(:,:,:,:,ie), 2*numlev, 0, ie)
529 : end do
530 0 : if(iam < par%nprocs) then
531 0 : call bndry_exchange(par, edgebuf,location=subname)
532 : end if
533 :
534 0 : do ie = 1, nelemd
535 0 : call edgeVunpack(edgebuf, fld_tmp(:,:,:,:,ie), 2*numlev, 0, ie)
536 : end do
537 0 : call freeEdgeBuffer(edgebuf)
538 : !check if fill values are present:
539 0 : usefillvalues = any(fld_tmp==fillvalue)
540 : end if
541 0 : deallocate(fld_dyn)
542 : else
543 : !
544 : ! not physics decomposition
545 : !
546 : !check if fill values are present:
547 0 : usefillvalues = (any(fldu(1:nsize:1,nsize,:)==fillvalue) .or. any(fldv(1:nsize:1,nsize,:)==fillvalue))
548 0 : do ie = 1, nelemd
549 0 : fld_tmp(1:nsize,1:nsize,1,1:numlev,ie) = RESHAPE(fldu(1:nsize*nsize,1:numlev,ie),(/nsize,nsize,numlev/))
550 0 : fld_tmp(1:nsize,1:nsize,2,1:numlev,ie) = RESHAPE(fldv(1:nsize*nsize,1:numlev,ie),(/nsize,nsize,numlev/))
551 : end do
552 : endif
553 : !
554 : !***************************************************************************
555 : !
556 : ! now data is on dynamics decomposition
557 : !
558 : !***************************************************************************
559 : !
560 0 : if (decomp_type==fvm_decomp.or.(fv_nphys>0.and.decomp_type==phys_decomp)) then
561 : !
562 : !***************************************************************************
563 : !
564 : ! code for non-GLL grids: need to fill halo and interpolate
565 : ! (if on panel edge/corner) for bilinear interpolation
566 : !
567 : !***************************************************************************
568 : !
569 :
570 : !JMD $OMP PARALLEL NUM_THREADS(horz_num_threads), DEFAULT(SHARED), PRIVATE(hybrid,nets,nete,n)
571 : !JMD hybrid = config_thread_region(par,'horizontal')
572 0 : hybrid = config_thread_region(par,'serial')
573 0 : call get_loop_ranges(hybrid,ibeg=nets,iend=nete)
574 0 : call fill_halo_and_extend_panel(elem(nets:nete),fvm(nets:nete),&
575 0 : fld_tmp(:,:,:,:,nets:nete),hybrid,nets,nete,nsize,nhcc,nhalo,2,numlev,.true.,.false.)
576 0 : do ie=nets,nete
577 0 : call vec_latlon_to_contra(elem(ie),nsize,nhcc,numlev,fld_tmp(:,:,:,:,ie),fvm(ie))
578 : end do
579 0 : call fill_halo_and_extend_panel(elem(nets:nete),fvm(nets:nete),&
580 0 : fld_tmp(:,:,:,:,nets:nete),hybrid,nets,nete,nsize,nhcc,nhalo,2,numlev,.false.,.true.)
581 :
582 : !check if fill values are present:
583 0 : usefillvalues = any(fld_tmp(:,:,:,:,nets:nete) == fillvalue)
584 : else
585 0 : do ie=1,nelemd
586 0 : call vec_latlon_to_contra(elem(ie),nsize,nhcc,numlev,fld_tmp(:,:,:,:,ie))
587 : end do
588 : end if
589 : !
590 : ! WARNING - 1:nelemd and nets:nete
591 : !
592 : !!$OMP MASTER !JMD
593 0 : dest(:,:,:,:,1:nelemd) = fld_tmp(1-nhalo:nsize+nhalo,1-nhalo:nsize+nhalo,:,:,1:nelemd)
594 : !!$OMP END MASTER
595 0 : deallocate(fld_tmp)
596 : !
597 : !***************************************************************************
598 : !
599 : ! do mapping from source grid to latlon grid
600 : !
601 : !***************************************************************************
602 : !
603 0 : ncnt_out = sum(cam_interpolate(1:nelemd)%n_interp)
604 0 : allocate(fldout(ncnt_out,numlev,2))
605 0 : allocate(idof(ncnt_out*numlev))
606 :
607 0 : fldout = -999_r8
608 0 : idof = 0
609 0 : st = 1
610 0 : do ie=1,nelemd
611 0 : ncol = cam_interpolate(ie)%n_interp
612 0 : do k=0,numlev-1
613 0 : do i=1,ncol
614 0 : idof(st+i-1+k*ncnt_out)=cam_interpolate(ie)%ilon(i)+nlon*(cam_interpolate(ie)%ilat(i)-1)+nlon*nlat*k
615 : enddo
616 : enddo
617 : ! Now that we have the field on the dyn grid we need to interpolate
618 0 : en = st+cam_interpolate(ie)%n_interp-1
619 0 : if(usefillvalues) then
620 0 : call interpolate_scalar(cam_interpolate(ie),dest(:,:,1,:,ie), nsize, nhalo, numlev, fldout(st:en,:,1), fillvalue)
621 0 : call interpolate_scalar(cam_interpolate(ie),dest(:,:,2,:,ie), nsize, nhalo, numlev, fldout(st:en,:,2), fillvalue)
622 : else
623 0 : call interpolate_scalar(cam_interpolate(ie),dest(:,:,1,:,ie), nsize, nhalo, numlev, fldout(st:en,:,1))
624 0 : call interpolate_scalar(cam_interpolate(ie),dest(:,:,2,:,ie), nsize, nhalo, numlev, fldout(st:en,:,2))
625 : end if
626 : !
627 : ! convert from contravariant components to lat-lon
628 : !
629 0 : do i=1,cam_interpolate(ie)%n_interp
630 : ! convert fld from contra->latlon
631 0 : call dmap(D,cam_interpolate(ie)%interp_xy(i)%x,cam_interpolate(ie)%interp_xy(i)%y,&
632 0 : elem(ie)%corners3D,cubed_sphere_map,elem(ie)%corners,elem(ie)%u2qmap,elem(ie)%facenum)
633 : ! convert fld from contra->latlon
634 0 : do k=1,numlev
635 0 : v1 = fldout(st+i-1,k,1)
636 0 : v2 = fldout(st+i-1,k,2)
637 0 : fldout(st+i-1,k,1)=D(1,1)*v1 + D(1,2)*v2
638 0 : fldout(st+i-1,k,2)=D(2,1)*v1 + D(2,2)*v2
639 : end do
640 : end do
641 0 : st = en+1
642 : end do
643 :
644 0 : if(numlev==1) then
645 0 : call pio_initdecomp(pio_subsystem, data_type, (/nlon,nlat/), idof, iodesc)
646 : else
647 0 : call pio_initdecomp(pio_subsystem, data_type, (/nlon,nlat,numlev/), idof, iodesc)
648 : end if
649 :
650 0 : if(data_type == pio_real) then
651 0 : call pio_write_darray(File, varidu, iodesc, real(fldout(:,:,1), r4), ierr)
652 0 : call pio_write_darray(File, varidv, iodesc, real(fldout(:,:,2), r4), ierr)
653 : else
654 0 : call pio_write_darray(File, varidu, iodesc, fldout(:,:,1), ierr)
655 0 : call pio_write_darray(File, varidv, iodesc, fldout(:,:,2), ierr)
656 : end if
657 :
658 :
659 0 : deallocate(fldout)
660 0 : deallocate(idof)
661 0 : deallocate(dest)
662 0 : call pio_freedecomp(file,iodesc)
663 :
664 0 : end subroutine write_interpolated_vector
665 :
666 0 : end module interp_mod
|
