Line data Source code
1 : module interpolate_mod
2 : use shr_kind_mod, only: r8=>shr_kind_r8
3 : use element_mod, only: element_t
4 : use dimensions_mod, only: np, ne, nelemd, nc, nhe, nhc
5 : use quadrature_mod, only: quadrature_t, legendre, quad_norm
6 : use coordinate_systems_mod, only: spherical_polar_t, cartesian2d_t, &
7 : cartesian3D_t, sphere2cubedsphere, spherical_to_cart, &
8 : cubedsphere2cart, distance, change_coordinates, projectpoint
9 : use physconst, only: PI
10 : use quadrature_mod, only: quadrature_t, gauss, gausslobatto
11 : use parallel_mod, only: syncmp, parallel_t
12 : use cam_abortutils, only: endrun
13 : use spmd_utils, only: MPI_MAX, MPI_SUM, MPI_MIN, mpi_real8, MPI_integer
14 : use cube_mod, only: convert_gbl_index, dmap, ref2sphere
15 : use mesh_mod, only: MeshUseMeshFile
16 : use control_mod, only: cubed_sphere_map
17 : use cam_logfile, only: iulog
18 : use string_utils, only: int2str
19 :
20 : implicit none
21 : private
22 : save
23 :
24 : logical :: debug=.false.
25 :
26 : type, public :: interpolate_t
27 : real (kind=r8), dimension(:,:), pointer :: Imat ! P_k(xj)*wj/gamma(k)
28 : real (kind=r8), dimension(:) , pointer :: rk ! 1/k
29 : real (kind=r8), dimension(:) , pointer :: vtemp ! temp results
30 : real (kind=r8), dimension(:) , pointer :: glp ! GLL pts (nair)
31 : end type interpolate_t
32 :
33 : type, public :: interpdata_t
34 : ! Output Interpolation points. Used to output data on lat-lon (or other grid)
35 : ! with native element interpolation. Each element keeps a list of points from the
36 : ! interpolation grid that are in this element
37 : type (cartesian2D_t),pointer,dimension(:):: interp_xy ! element coordinate
38 : integer, pointer,dimension(:) :: ilat,ilon ! position of interpolation point in lat-lon grid
39 : integer :: n_interp
40 : integer :: nlat
41 : integer :: nlon
42 : logical :: first_entry = .TRUE.
43 : end type interpdata_t
44 :
45 : real (kind=r8), private :: delta = 1.0e-9_r8 ! move tiny bit off center to
46 : ! avoid landing on element edges
47 :
48 :
49 : ! static data for interp_tracers
50 : logical :: interp_tracers_init=.false.
51 : real (kind=r8 ) :: interp_c(np,np)
52 : real (kind=r8 ) :: interp_gll(np)
53 :
54 : public :: interp_init
55 : public :: setup_latlon_interp
56 : public :: interpolate_scalar
57 : public :: interpolate_ce
58 :
59 : public :: interpol_phys_latlon
60 : public :: interpolate_vector
61 : public :: set_interp_parameter
62 : public :: get_interp_parameter
63 : public :: get_interp_gweight
64 : public :: get_interp_lat
65 : public :: get_interp_lon
66 : public :: cube_facepoint_ne
67 : public :: cube_facepoint_unstructured
68 : public :: parametric_coordinates
69 :
70 : public :: interpolate_tracers
71 : public :: interpolate_tracers_init
72 : public :: minmax_tracers
73 : public :: interpolate_2d
74 : public :: interpolate_create
75 : public :: point_inside_quad
76 : public :: vec_latlon_to_contra
77 :
78 :
79 : interface interpolate_scalar
80 : module procedure interpolate_scalar2d
81 : module procedure interpolate_scalar3d
82 : end interface
83 : interface interpolate_vector
84 : module procedure interpolate_vector2d
85 : module procedure interpolate_vector3d
86 : end interface
87 :
88 : type (interpolate_t), target :: interp_p
89 :
90 : ! store the lat-lon grid
91 : ! gridtype = 1 equally spaced, including poles (FV scalars output grid)
92 : ! gridtype = 2 Gauss grid (CAM Eulerian)
93 : ! gridtype = 3 equally spaced, no poles (FV staggered velocity)
94 : ! Seven possible history files, last one is inithist and should be native grid
95 : integer :: nlat,nlon
96 : real (kind=r8), pointer, public :: lat(:) => NULL()
97 : real (kind=r8), pointer, public :: lon(:) => NULL()
98 : real (kind=r8), pointer, public :: gweight(:) => NULL()
99 : integer :: gridtype = 1 !
100 : integer :: itype = 1 ! 0 = native high order
101 : ! 1 = bilinear
102 :
103 : integer :: auto_grid = 0 ! 0 = interpolation grid set by namelist
104 : ! 1 = grid set via mesh resolution
105 :
106 :
107 : ! static data, used by bilin_phys2gll()
108 : ! shared by all threads. only allocate if subroutine will be used
109 : !JMD integer :: nphys_init=0
110 : !JMD integer :: index_l(np),index_r(np)
111 : !JMD real(kind=r8),allocatable :: weights(:,:,:,:,:) ! np,np,2,2,nelemd
112 :
113 : !JMD public :: bilin_phys2gll
114 : !JMD public :: bilin_phys2gll_init
115 : contains
116 :
117 :
118 0 : subroutine set_interp_parameter(parm_name, value)
119 : character*(*), intent(in) :: parm_name
120 : character(len=80) :: msg
121 : integer :: value,power
122 : real (kind=r8) :: value_target
123 :
124 0 : if(parm_name .eq. 'itype') then
125 0 : itype=value
126 0 : else if(parm_name .eq. 'nlon') then
127 0 : nlon=value
128 0 : else if(parm_name .eq. 'nlat') then
129 0 : nlat=value
130 0 : else if(parm_name.eq. 'gridtype') then
131 0 : gridtype=value
132 0 : else if(parm_name.eq. 'auto') then
133 0 : auto_grid=1
134 : ! compute recommended nlat,nlon which has slightly higher
135 : ! resolution than the specifed number of points around equator given in "value"
136 : ! computed recommended lat-lon grid.
137 : ! nlon > peq peq = points around equator cubed sphere grid
138 : ! take nlon power of 2, and at most 1 power of 3
139 0 : if (value.eq.0) then
140 : ! If reading in unstructured mesh, ne = 0
141 : ! This makes it hard to guess how many interpolation points to use
142 : ! So We'll set the default as 720 x 360
143 : ! BUT if you're running with an unstructured mesh, set interp_nlon and interp_nlat
144 0 : nlon = 1536
145 0 : nlat = 768
146 : else
147 0 : value_target=value*1.25_r8
148 0 : power = nint(0.5_r8 + log( value_target)/log(2.0_r8) )
149 0 : power = max(power,7) ! min grid: 64x128
150 0 : if ( 3*2**(power-2) > value_target) then
151 0 : nlon=3*2**(power-2) ! use 1 power of 3
152 : else
153 0 : nlon=2**power
154 : endif
155 : endif
156 0 : nlat=nlon/2
157 0 : if (gridtype==1) nlat=nlat+1
158 : else
159 0 : write(msg,*) 'Did not recognize parameter named ',parm_name,' in interpolate_mod:set_interp_parameter'
160 0 : call endrun(msg)
161 : end if
162 0 : end subroutine set_interp_parameter
163 0 : function get_interp_parameter(parm_name) result(value)
164 : character*(*), intent(in) :: parm_name
165 : integer :: value
166 : character(len=80) :: msg
167 0 : if(parm_name .eq. 'itype') then
168 0 : value=itype
169 0 : else if(parm_name .eq. 'nlon') then
170 0 : value=nlon
171 0 : else if(parm_name .eq. 'nlat') then
172 0 : value=nlat
173 0 : else if(parm_name.eq. 'gridtype') then
174 0 : value=gridtype
175 0 : else if(parm_name.eq. 'auto_grid') then
176 0 : value=auto_grid
177 : else
178 0 : write(msg,*) 'Did not recognize parameter named ',parm_name,' in interpolate_mod:get_interp_parameter'
179 0 : value=-1
180 0 : call endrun(msg)
181 : end if
182 : return
183 : end function get_interp_parameter
184 0 : function get_interp_gweight() result(gw)
185 : real(kind=r8) :: gw(nlat)
186 0 : gw=gweight
187 0 : return
188 0 : end function get_interp_gweight
189 0 : function get_interp_lat() result(thislat)
190 : real(kind=r8) :: thislat(nlat)
191 0 : thislat=lat*180.0_r8/PI
192 0 : return
193 0 : end function get_interp_lat
194 0 : function get_interp_lon() result(thislon)
195 : real(kind=r8) :: thislon(nlon)
196 0 : thislon=lon*180.0_r8/PI
197 0 : return
198 0 : end function get_interp_lon
199 :
200 5216400 : subroutine interpolate_create(gquad,interp)
201 : type (quadrature_t) , intent(in) :: gquad
202 : type (interpolate_t), intent(out) :: interp
203 :
204 :
205 : ! Local variables
206 :
207 : integer k,j
208 : integer npts
209 5216400 : real (kind=r8), dimension(:), allocatable :: gamma
210 5216400 : real (kind=r8), dimension(:), allocatable :: leg
211 :
212 5216400 : npts = size(gquad%points)
213 :
214 20865600 : allocate(interp%Imat(npts,npts))
215 15649200 : allocate(interp%rk(npts))
216 10432800 : allocate(interp%vtemp(npts))
217 10432800 : allocate(interp%glp(npts))
218 10432800 : allocate(gamma(npts))
219 10432800 : allocate(leg(npts))
220 :
221 5216400 : gamma = quad_norm(gquad,npts)
222 :
223 26082000 : do k=1,npts
224 20865600 : interp%rk(k) = 1.0_r8/k
225 26082000 : interp%glp(k) = gquad%points(k) !nair
226 : end do
227 :
228 26082000 : do j=1,npts
229 20865600 : leg=legendre(gquad%points(j),npts-1)
230 109544400 : do k=1,npts
231 104328000 : interp%Imat(j,k)=leg(k)*gquad%weights(j)/gamma(k)
232 : end do
233 : end do
234 :
235 5216400 : deallocate(gamma)
236 5216400 : deallocate(leg)
237 :
238 5216400 : end subroutine interpolate_create
239 :
240 :
241 0 : subroutine interpolate_tracers_init()
242 : use dimensions_mod, only : np, qsize
243 : use quadrature_mod, only : quadrature_t, gausslobatto
244 :
245 :
246 : implicit none
247 :
248 : type (quadrature_t ) :: gll
249 : real (kind=r8 ) :: dp (np)
250 : integer :: i,j
251 :
252 0 : gll=gausslobatto(np)
253 0 : dp = 1
254 0 : do i=1,np
255 0 : do j=1,np
256 0 : if (i /= j) then
257 0 : dp(i) = dp(i) * (gll%points(i) - gll%points(j))
258 : end if
259 : end do
260 : end do
261 0 : do i=1,np
262 0 : do j=1,np
263 0 : interp_c(i,j) = 1/(dp(i)*dp(j))
264 : end do
265 : end do
266 0 : interp_gll(:) = gll%points(:)
267 0 : interp_tracers_init = .true.
268 :
269 0 : deallocate(gll%points)
270 0 : deallocate(gll%weights)
271 :
272 :
273 0 : end subroutine interpolate_tracers_init
274 :
275 :
276 :
277 :
278 0 : subroutine interpolate_tracers(r, tracers, f)
279 : use dimensions_mod, only : np, qsize
280 :
281 :
282 : implicit none
283 : type (cartesian2D_t), intent(in) :: r
284 : real (kind=r8),intent(in) :: tracers(np*np,qsize)
285 : real (kind=r8),intent(out) :: f(qsize)
286 :
287 : real (kind=r8 ) :: x (np)
288 : real (kind=r8 ) :: y (np)
289 : real (kind=r8 ) :: xy (np*np)
290 :
291 : integer :: i,j
292 :
293 :
294 0 : if (.not. interp_tracers_init ) then
295 0 : call endrun('ERROR: interpolate_tracers() was not initialized')
296 : endif
297 :
298 0 : x = 1
299 0 : y = 1
300 0 : do i=1,np
301 0 : do j=1,np
302 0 : if (i /= j) then
303 0 : x(i) = x(i) * (r%x - interp_gll(j))
304 0 : y(i) = y(i) * (r%y - interp_gll(j))
305 : end if
306 : end do
307 : end do
308 :
309 0 : do j=1,np
310 0 : do i=1,np
311 0 : xy(i + (j-1)*np) = x(i)*y(j)*interp_c(i,j)
312 : end do
313 : end do
314 0 : f = MATMUL(xy,tracers)
315 0 : end subroutine interpolate_tracers
316 :
317 :
318 :
319 : subroutine linear_interpolate_2d(x,y,s,v)
320 : use dimensions_mod, only : np, qsize
321 :
322 : real(kind=r8) , intent(in) :: x(np)
323 : real(kind=r8), intent(in) :: y(np,np,qsize)
324 : type (cartesian2D_t), intent(in) :: s
325 : real(kind=r8), intent(inout) :: v(qsize)
326 :
327 : integer :: i,j,q
328 : real (kind=r8) dx, dy(qsize), dydx(qsize)
329 : real (kind=r8) y0(qsize), y1(qsize)
330 : type (cartesian2D_t) :: r
331 :
332 : r = s
333 : if (r%x < -1) r%x = -1
334 : if (r%y < -1) r%y = -1
335 : if ( 1 < r%x) r%x = 1
336 : if ( 1 < r%y) r%y = 1
337 : do i=1,np
338 : if (r%x < x(i)) exit
339 : end do
340 : do j=1,np
341 : if (r%y < x(j)) exit
342 : end do
343 : if (1 < i) i = i-1
344 : if (1 < j) j = j-1
345 : if (np==i) i = i-1
346 : if (np==j) j = j-1
347 :
348 : dx = x(i+1) - x(i)
349 : dy = y(i+1,j,:) - y(i,j,:)
350 : dydx = dy/dx
351 : y0 = y(i,j,:) + (r%x-x(i))*dydx
352 :
353 : dy = y(i+1,j+1,:) - y(i,j+1,:)
354 : dydx = dy/dx
355 : y1 = y(i,j+1,:) + (r%x-x(i))*dydx
356 :
357 : dx = x(j+1) - x(j)
358 : dy = y1 - y0
359 : dydx = dy/dx
360 : v = y0 + (r%y-x(j))*dydx
361 :
362 : end subroutine linear_interpolate_2d
363 :
364 0 : subroutine minmax_tracers(r, tracers, mint, maxt)
365 : use dimensions_mod, only : np, qsize
366 : use quadrature_mod, only : quadrature_t, gausslobatto
367 :
368 :
369 : implicit none
370 :
371 : type (cartesian2D_t), intent(in) :: r
372 : real (kind=r8),intent(in) :: tracers(np,np,qsize)
373 : real (kind=r8),intent(out) :: mint (qsize)
374 : real (kind=r8),intent(out) :: maxt (qsize)
375 :
376 : type (quadrature_t), save :: gll
377 : integer :: i,j
378 : logical , save :: first_time=.true.
379 : real (kind=r8) :: y1 (qsize)
380 : real (kind=r8) :: y2 (qsize)
381 0 : real (kind=r8) :: q_interp (4,qsize)
382 : type (cartesian2D_t) :: s
383 : real (kind=r8) :: delta
384 : integer :: q
385 :
386 0 : do q=1,qsize
387 0 : mint(q) = minval(tracers(:,:,q))
388 0 : maxt(q) = maxval(tracers(:,:,q))
389 : enddo
390 : return
391 :
392 : delta = 1._r8/8._r8
393 :
394 : if (first_time) then
395 : first_time = .false.
396 : gll=gausslobatto(np)
397 : end if
398 :
399 : do i=1,np
400 : if (r%x < gll%points(i)) exit
401 : end do
402 : do j=1,np
403 : if (r%y < gll%points(j)) exit
404 : end do
405 : if (1 < i) i = i-1
406 : if (1 < j) j = j-1
407 : if (np==i) i = i-1
408 : if (np==j) j = j-1
409 :
410 : ! mint(:) = minval(minval(tracers(i:i+1,j:j+1,:),1),1)
411 : ! maxt(:) = maxval(maxval(tracers(i:i+1,j:j+1,:),1),1)
412 :
413 : ! Or check this out:
414 : s = r
415 : s%x = s%x - delta
416 : s%y = s%y - delta
417 : call linear_interpolate_2d(gll%points,tracers,s,q_interp(1,:))
418 : s = r
419 : s%x = s%x + delta
420 : s%y = s%y - delta
421 : call linear_interpolate_2d(gll%points,tracers,s,q_interp(2,:))
422 : s = r
423 : s%x = s%x - delta
424 : s%y = s%y + delta
425 : call linear_interpolate_2d(gll%points,tracers,s,q_interp(3,:))
426 : s = r
427 : s%x = s%x + delta
428 : s%y = s%y + delta
429 : call linear_interpolate_2d(gll%points,tracers,s,q_interp(4,:))
430 :
431 : mint(:) = minval(q_interp(:,:),1)
432 : maxt(:) = maxval(q_interp(:,:),1)
433 : end subroutine minmax_tracers
434 :
435 8732253600 : function interpolate_2d(cart, f, interp, npts, fillvalue) result(fxy)
436 : integer, intent(in) :: npts
437 : type (cartesian2D_t), intent(in) :: cart
438 : real (kind=r8), intent(in) :: f(npts,npts)
439 : type (interpolate_t) :: interp
440 : real (kind=r8) :: fxy ! value of f interpolated to (x,y)
441 : real (kind=r8), intent(in), optional :: fillvalue
442 : ! local variables
443 :
444 : real (kind=r8) :: tmp_1,tmp_2
445 : real (kind=r8) :: fk0,fk1
446 : real (kind=r8) :: pk
447 :
448 : integer :: l,j,k
449 :
450 8732253600 : if(present(fillvalue)) then
451 0 : if (any(f==fillvalue)) then
452 8732253600 : fxy = fillvalue
453 : return
454 : endif
455 : endif
456 :
457 :
458 26196760800 : do l=1,npts,2
459 :
460 : ! Compute Pk(cart%x) for Legendre order 0
461 :
462 17464507200 : pk = 1.0_r8
463 :
464 17464507200 : fk0=0.0_r8
465 17464507200 : fk1=0.0_r8
466 87322536000 : do j=1,npts
467 69858028800 : fk0 = fk0 + interp%Imat(j,1)*f(j,l )
468 87322536000 : fk1 = fk1 + interp%Imat(j,1)*f(j,l+1)
469 : end do
470 17464507200 : interp%vtemp(l ) = pk*fk0
471 17464507200 : interp%vtemp(l+1) = pk*fk1
472 :
473 : ! Compute Pk(cart%x) for Legendre order 1
474 :
475 17464507200 : tmp_2 = pk
476 17464507200 : pk = cart%x
477 :
478 17464507200 : fk0=0.0_r8
479 17464507200 : fk1=0.0_r8
480 87322536000 : do j=1,npts
481 69858028800 : fk0 = fk0 + interp%Imat(j,2)*f(j,l )
482 87322536000 : fk1 = fk1 + interp%Imat(j,2)*f(j,l+1)
483 : end do
484 17464507200 : interp%vtemp(l ) = interp%vtemp(l ) + pk*fk0
485 17464507200 : interp%vtemp(l+1) = interp%vtemp(l+1) + pk*fk1
486 :
487 : ! Compute Pk(cart%x) for Legendre order 2 to npts-1
488 :
489 61125775200 : do k = 2,npts-1
490 :
491 34929014400 : tmp_1 = tmp_2
492 34929014400 : tmp_2 = pk
493 34929014400 : pk = ( (2*k-1)*cart%x*tmp_2 - (k-1)*tmp_1 )*interp%rk(k)
494 :
495 34929014400 : fk0=0.0_r8
496 34929014400 : fk1=0.0_r8
497 >17464*10^7 : do j=1,npts
498 >13971*10^7 : fk0 = fk0 + interp%Imat(j,k+1)*f(j,l )
499 >17464*10^7 : fk1 = fk1 + interp%Imat(j,k+1)*f(j,l+1)
500 : end do
501 34929014400 : interp%vtemp(l ) = interp%vtemp(l ) + pk*fk0
502 52393521600 : interp%vtemp(l+1) = interp%vtemp(l+1) + pk*fk1
503 :
504 : end do
505 :
506 : end do
507 :
508 : ! Compute Pk(cart%y) for Legendre order 0
509 :
510 8732253600 : pk = 1.0_r8
511 :
512 8732253600 : fk0 = 0.0_r8
513 43661268000 : do j=1,npts
514 43661268000 : fk0 = fk0 + interp%Imat(j,1)*interp%vtemp(j)
515 : end do
516 8732253600 : fxy = pk*fk0
517 :
518 : ! Compute Pk(cart%y) for Legendre order 1
519 :
520 8732253600 : tmp_2 = pk
521 8732253600 : pk = cart%y
522 :
523 8732253600 : fk0=0.0_r8
524 43661268000 : do j=1,npts
525 43661268000 : fk0 = fk0 + interp%Imat(j,2)*interp%vtemp(j)
526 : end do
527 8732253600 : fxy = fxy + pk*fk0
528 :
529 : ! Compute Pk(cart%y) for Legendre order 2, npts-1
530 :
531 26196760800 : do k = 2,npts-1
532 17464507200 : tmp_1 = tmp_2
533 17464507200 : tmp_2 = pk
534 17464507200 : pk = ( (2*k-1)*cart%y*tmp_2 - (k-1)*tmp_1 )*interp%rk(k)
535 :
536 17464507200 : fk0 = 0.0_r8
537 87322536000 : do j=1,npts
538 87322536000 : fk0 = fk0 + interp%Imat(j,k+1)*interp%vtemp(j)
539 : end do
540 :
541 26196760800 : fxy = fxy + pk*fk0
542 :
543 : end do
544 :
545 : end function interpolate_2d
546 :
547 : !===============================
548 : !(Nair) Bilinear interpolation for every GLL grid cell
549 : !===============================
550 :
551 0 : function interpol_bilinear(cart, f, xoy, imin, imax, fillvalue) result(fxy)
552 : integer, intent(in) :: imin,imax
553 : type (cartesian2D_t), intent(in) :: cart
554 : real (kind=r8), intent(in) :: f(imin:imax,imin:imax)
555 : real (kind=r8) :: xoy(imin:imax)
556 : real (kind=r8) :: fxy ! value of f interpolated to (x,y)
557 : real (kind=r8), intent(in), optional :: fillvalue
558 : ! local variables
559 :
560 : real (kind=r8) :: p,q,xp,yp ,y4(4)
561 : integer :: l,j,k, ii, jj, na,nb,nm
562 :
563 0 : xp = cart%x
564 0 : yp = cart%y
565 :
566 : ! Search index along "x" (bisection method)
567 :
568 0 : na = imin
569 0 : nb = imax
570 : do
571 0 : if ((nb-na) <= 1) exit
572 0 : nm = (nb + na)/2
573 0 : if (xp > xoy(nm)) then
574 : na = nm
575 : else
576 0 : nb = nm
577 : endif
578 : enddo
579 : ii = na
580 :
581 : ! Search index along "y"
582 :
583 : na = imin
584 : nb = imax
585 : do
586 0 : if ((nb-na) <= 1) exit
587 0 : nm = (nb + na)/2
588 0 : if (yp > xoy(nm)) then
589 : na = nm
590 : else
591 0 : nb = nm
592 : endif
593 : enddo
594 0 : jj = na
595 :
596 : ! GLL cell containing (xp,yp)
597 :
598 0 : y4(1) = f(ii,jj)
599 0 : y4(2) = f(ii+1,jj)
600 0 : y4(3) = f(ii+1,jj+1)
601 0 : y4(4) = f(ii,jj+1)
602 :
603 0 : if(present(fillvalue)) then
604 0 : if (any(y4==fillvalue)) then
605 0 : fxy = fillvalue
606 : return
607 : endif
608 : endif
609 :
610 0 : p = (xp - xoy(ii))/(xoy(ii+1) - xoy(ii))
611 0 : q = (yp - xoy(jj))/(xoy(jj+1) - xoy(jj))
612 :
613 : fxy = (1.0_r8 - p)*(1.0_r8 - q)* y4(1) + p*(1.0_r8 - q) * y4(2) &
614 0 : + p*q* y4(3) + (1.0_r8 - p)*q * y4(4)
615 0 : end function interpol_bilinear
616 :
617 : ! ----------------------------------------------------------------------------------!
618 : !FUNCTION interpol_phys_latlon----------------------------------------CE-for fvm!
619 : ! AUTHOR: CHRISTOPH ERATH, 23. May 2012 !
620 : ! DESCRIPTION: evaluation of the reconstruction for every physics grid cell !
621 : ! !
622 : ! CALLS:
623 : ! INPUT:
624 : !
625 : ! OUTPUT:
626 : !-----------------------------------------------------------------------------------!
627 0 : subroutine interpol_phys_latlon(interpdata,f, fvm, corners, desc, flatlon,lmono)
628 : use fvm_control_volume_mod, only : fvm_struct
629 : ! use fvm_reconstruction_mod, only: reconstruction_gradient, recons_val_cart
630 : use edgetype_mod, only : edgedescriptor_t
631 :
632 : type (interpdata_t), intent(in) :: interpdata
633 : real (kind=r8), intent(inout) :: f(1-nhc:nc+nhc,1-nhc:nc+nhc)
634 : type (fvm_struct), intent(in) :: fvm
635 : type (cartesian2d_t), intent(in) :: corners(:)
636 : type (edgedescriptor_t),intent(in) :: desc
637 : logical, intent(in) :: lmono
638 :
639 : real (kind=r8) :: flatlon(:)
640 : ! local variables
641 : real (kind=r8) :: xp,yp, tmpval
642 : real (kind=r8) :: tmpaxp,tmpaxm, tmpayp, tmpaym
643 : integer :: i, ix, jy, starti,endi,tmpi
644 : real (kind=r8), dimension(1-nhe:nc+nhe,1-nhe:nc+nhe,6) :: recons
645 :
646 : real (kind=r8), dimension(nc+1) :: x, y
647 :
648 : ! call reconstruction_gradient(f, fvm,recons,6,lmono)
649 : ! recons=0.0 ! PCoM
650 :
651 0 : x(1:nc) = fvm%vtx_cart(1,1,1:nc,1 )
652 0 : y(1:nc) = fvm%vtx_cart(1,2,1 ,1:nc)
653 0 : x(nc+1) = fvm%vtx_cart(2,1,nc,1 )
654 0 : y(nc+1) = fvm%vtx_cart(3,2,1 ,nc )
655 :
656 0 : tmpaxp=(corners(1)%x+corners(2)%x)/2
657 0 : tmpaxm=(corners(2)%x-corners(1)%x)/2
658 0 : tmpayp=(corners(1)%y+corners(4)%y)/2
659 0 : tmpaym=(corners(4)%y-corners(1)%y)/2
660 0 : do i=1,interpdata%n_interp
661 : ! caculation phys grid coordinate of xp point, note the interp_xy are on the reference [-1,1]x[-1,1]
662 0 : xp=tan(tmpaxp+interpdata%interp_xy(i)%x*tmpaxm)
663 0 : yp=tan(tmpayp+interpdata%interp_xy(i)%y*tmpaym)
664 :
665 : ! Search index along "x" (bisection method)
666 0 : starti = 1
667 0 : endi = nc+1
668 : do
669 0 : if ((endi-starti) <= 1) exit
670 0 : tmpi = (endi + starti)/2
671 0 : if (xp > x(tmpi)) then
672 : starti = tmpi
673 : else
674 0 : endi = tmpi
675 : endif
676 : enddo
677 : ix = starti
678 :
679 : ! Search index along "y"
680 : starti = 1
681 : endi = nc+1
682 : do
683 0 : if ((endi-starti) <= 1) exit
684 0 : tmpi = (endi + starti)/2
685 0 : if (yp > y(tmpi)) then
686 : starti = tmpi
687 : else
688 0 : endi = tmpi
689 : endif
690 : enddo
691 0 : jy = starti
692 :
693 : ! call recons_val_cart(f(ix,jy), xp,yp, fvm%spherecentroid(ix,jy,:), fvm%recons_metrics(ix,jy,:), &
694 : ! recons(ix,jy,:), tmpval)
695 0 : tmpval=f(ix,jy)
696 0 : flatlon(i)=tmpval
697 : !phl PCoM
698 : ! flatlon(i)=f(ix,jy)
699 : end do
700 0 : end subroutine interpol_phys_latlon
701 :
702 0 : function parametric_coordinates(sphere, corners3D,ref_map_in, corners,u2qmap,facenum) result (ref)
703 : implicit none
704 : type (spherical_polar_t), intent(in) :: sphere
705 : type (cartesian2D_t) :: ref
706 :
707 : type (cartesian3D_t) :: corners3D(4) !x,y,z coords of element corners
708 : integer,optional :: ref_map_in ! default is global variable 'cubed_sphere_map'
709 : ! optional arguments, only needed for ref_map=1 (equi-angle gnomonic projection):
710 : type (cartesian2D_t),optional :: corners(4) ! gnomonic coords of element corners
711 : real (kind=r8),optional :: u2qmap(4,2)
712 : integer,optional :: facenum
713 :
714 :
715 : ! local
716 : integer :: i, MAX_NR_ITER=10
717 : real(kind=r8) :: D(2,2),Dinv(2,2),detD,a,b,resa,resb,dela,delb,costh
718 : real(kind=r8) :: tol_sq = 1.0e-26_r8
719 : type (spherical_polar_t) :: sphere1, sphere_tmp
720 : integer :: ref_map
721 :
722 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
723 : ! newton iteration on: ref=ref - df^-1 (ref2sphere(ref) - sphere)
724 : !
725 : ! Generic version written in terms of HOMME's 'ref2sphere' and 'Dmap' operaters,
726 : ! with no assumption as to the type of map (gnomonic, equi-angular, parametric)
727 : !
728 : ! Note that the coordinate increment from newton iterations is not a direction and thus
729 : ! should not be converted into motion along a great circle arc - this routine
730 : ! correclty applies the increment by just adding it to the coordintes
731 : !
732 : ! f = ref2sphere(xvec) - sphere
733 : ! df = d(ref2sphere)
734 : !
735 : ! D = diag(cos(theta),1) * d(ref2sphere) d(ref2sphere) = diag(1/cos(theta),1)*D
736 : !
737 : ! df = diag(1/cos(theta),1)*D
738 : ! df^-1 = D^-1 * diag(cos(theta),1)
739 : !
740 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
741 0 : if (present(ref_map_in)) then
742 0 : ref_map=ref_map_in
743 : else
744 0 : ref_map=cubed_sphere_map
745 : endif
746 0 : costh=cos(sphere%lat)
747 0 : a=0
748 0 : b=0
749 0 : i=0
750 : do
751 0 : sphere1 = ref2sphere(a,b,corners3D,ref_map,corners,facenum)
752 0 : resa = sphere1%lon - sphere%lon
753 0 : if (resa > pi) resa= resa - 2*pi
754 0 : if (resa < -pi) resa= resa + 2*pi
755 :
756 0 : resb = sphere1%lat - sphere%lat
757 :
758 0 : call Dmap(D,a,b,corners3D,ref_map,corners,u2qmap,facenum)
759 0 : detD = D(1,1)*D(2,2) - D(1,2)*D(2,1)
760 0 : Dinv(1,1) = D(2,2)/detD
761 0 : Dinv(1,2) = -D(1,2)/detD
762 0 : Dinv(2,1) = -D(2,1)/detD
763 0 : Dinv(2,2) = D(1,1)/detD
764 :
765 0 : dela = Dinv(1,1)*costh*resa + Dinv(1,2)*resb
766 0 : delb = Dinv(2,1)*costh*resa + Dinv(2,2)*resb
767 0 : a = a - dela
768 0 : b = b - delb
769 0 : i=i+1
770 0 : if ( (costh*resa)**2 + resb**2 < tol_sq .or. MAX_NR_ITER < i) exit
771 : end do
772 0 : ref%x=a
773 0 : ref%y=b
774 :
775 0 : end function parametric_coordinates
776 :
777 :
778 :
779 :
780 : !
781 : ! find element containing given point, useing HOMME's standard
782 : ! equi-angular gnomonic map.
783 : ! note that with this map, only coordinate lines are great circle arcs
784 : !
785 0 : function point_inside_equiangular(elem, sphere, sphere_xyz) result(inside)
786 : implicit none
787 : type (spherical_polar_t), intent(in) :: sphere
788 : type (cartesian3D_t), intent(in) :: sphere_xyz
789 : type (element_t) , intent(in) :: elem
790 : logical :: inside, inside2
791 : integer :: i,j
792 : type (cartesian2D_t) :: corners(4),sphere_xy,cart
793 : type (cartesian3D_t) :: corners_xyz(4),center,a,b,cross(4)
794 : real (kind=r8) :: yp(4), y, elem_diam,dotprod
795 : real (kind=r8) :: xp(4), x, xc,yc
796 : real (kind=r8) :: tol_inside
797 : real (kind=r8) :: d1,d2
798 :
799 : type (spherical_polar_t) :: sphere_tmp
800 :
801 0 : inside = .false.
802 :
803 :
804 : ! first check if point is near the element:
805 0 : corners_xyz(:) = elem%corners3D(:)
806 : elem_diam = max( distance(corners_xyz(1),corners_xyz(3)), &
807 0 : distance(corners_xyz(2),corners_xyz(4)) )
808 :
809 0 : center%x = sum(corners_xyz(1:4)%x)/4
810 0 : center%y = sum(corners_xyz(1:4)%y)/4
811 0 : center%z = sum(corners_xyz(1:4)%z)/4
812 0 : if ( distance(center,sphere_xyz) > 1.0_r8*elem_diam ) return
813 :
814 0 : tol_inside = 1.0e-10_r8*elem_diam**2
815 : ! the point is close to the element, so project both to cubed sphere
816 : ! and perform contour integral
817 0 : sphere_xy=sphere2cubedsphere(sphere,elem%FaceNum)
818 0 : x = sphere_xy%x
819 0 : y = sphere_xy%y
820 0 : do i=1,4
821 0 : xp(i) = elem%corners(i)%x
822 0 : yp(i) = elem%corners(i)%y
823 : end do
824 :
825 :
826 0 : if (debug) then
827 0 : print *,'point: ',x,y,elem%FaceNum
828 0 : print *,'element:'
829 0 : write(*,'(a,4e16.8,a)') 'x=[',xp(1:4),']'
830 0 : write(*,'(a,4e16.8,a)') 'y=[',yp(1:4),']'
831 :
832 : ! first check if centroid is in this element (sanity check)
833 0 : sphere_tmp=change_coordinates(center)
834 0 : sphere_xy=sphere2cubedsphere(sphere_tmp,elem%FaceNum)
835 0 : xc=sphere_xy%x
836 0 : yc=sphere_xy%y
837 0 : print *,'cross product with centroid: all numbers should be negative'
838 0 : j = 4
839 0 : do i=1,4
840 0 : print *,i,(xc-xp(j))*(yp(i)-yp(j)) - (yc-yp(j))*(xp(i)-xp(j))
841 0 : j = i ! within this loopk j = i-1
842 : end do
843 :
844 0 : print *,'cross product with search point'
845 0 : j = 4
846 0 : do i=1,4
847 0 : print *,i,(x-xp(j))*(yp(i)-yp(j)) - (y-yp(j))*(xp(i)-xp(j))
848 0 : j = i ! within this loopk j = i-1
849 : end do
850 : endif
851 :
852 :
853 0 : j = 4
854 0 : do i=1,4
855 : ! a = x-xp(j), y-yp(j)
856 : ! b = xp(i)-xp(j), yp(i)-yp(j)
857 : ! compute a cross b:
858 0 : if ( -( (x-xp(j))*(yp(i)-yp(j)) - (y-yp(j))*(xp(i)-xp(j))) > tol_inside ) then
859 : return
860 : endif
861 0 : j = i ! within this loopk j = i-1
862 : end do
863 : ! all cross products were negative, must be inside:
864 0 : inside=.true.
865 : end function point_inside_equiangular
866 :
867 :
868 : !
869 : ! find if quad contains given point, with quad edges assumed to be great circle arcs
870 : ! this will work with any map where straight lines are mapped to great circle arcs.
871 : ! (thus it will fail on unstructured grids using the equi-angular gnomonic map)
872 : !
873 0 : function point_inside_quad(corners_xyz, sphere_xyz) result(inside)
874 : implicit none
875 : type (cartesian3D_t), intent(in) :: sphere_xyz
876 : type (cartesian3D_t) , intent(in) :: corners_xyz(4)
877 : logical :: inside, inside2
878 : integer :: i,j,ii
879 : type (cartesian2D_t) :: corners(4),sphere_xy,cart
880 : type (cartesian3D_t) :: center,a,b,cross(4)
881 : real (kind=r8) :: yp(4), y, elem_diam,dotprod
882 : real (kind=r8) :: xp(4), x
883 : real (kind=r8) :: d1,d2, tol_inside = 1.0e-12_r8
884 :
885 : type (spherical_polar_t) :: sphere ! debug
886 :
887 0 : inside = .false.
888 :
889 : ! first check if point is near the corners:
890 : elem_diam = max( distance(corners_xyz(1),corners_xyz(3)), &
891 0 : distance(corners_xyz(2),corners_xyz(4)) )
892 :
893 0 : center%x = sum(corners_xyz(1:4)%x)/4
894 0 : center%y = sum(corners_xyz(1:4)%y)/4
895 0 : center%z = sum(corners_xyz(1:4)%z)/4
896 0 : if ( distance(center,sphere_xyz) > 1.0_r8*elem_diam ) return
897 :
898 : j = 4
899 0 : do i=1,4
900 : ! outward normal to plane containing j->i edge: corner(i) x corner(j)
901 : ! sphere dot (corner(i) x corner(j) ) = negative if inside
902 0 : cross(i)%x = corners_xyz(i)%y*corners_xyz(j)%z - corners_xyz(i)%z*corners_xyz(j)%y
903 0 : cross(i)%y =-(corners_xyz(i)%x*corners_xyz(j)%z - corners_xyz(i)%z*corners_xyz(j)%x)
904 0 : cross(i)%z = corners_xyz(i)%x*corners_xyz(j)%y - corners_xyz(i)%y*corners_xyz(j)%x
905 : dotprod = cross(i)%x*sphere_xyz%x + cross(i)%y*sphere_xyz%y +&
906 0 : cross(i)%z*sphere_xyz%z
907 0 : j = i ! within this loopk j = i-1
908 :
909 : ! dot product is proportional to elem_diam. positive means outside,
910 : ! but allow machine precision tolorence:
911 0 : if (dotprod > tol_inside*elem_diam) return
912 : !if (dotprod > 0) return
913 : end do
914 0 : inside=.true.
915 : return
916 : end function point_inside_quad
917 :
918 : !
919 : ! find element containing given point, with element edges assumed to be great circle arcs
920 : ! this will work with any map where straight lines are mapped to great circle arcs.
921 : ! (thus it will fail on unstructured grids using the equi-angular gnomonic map)
922 : !
923 0 : function point_inside_gc(elem, sphere_xyz) result(inside)
924 : implicit none
925 : type (cartesian3D_t), intent(in) :: sphere_xyz
926 : type (element_t) , intent(in) :: elem
927 : logical :: inside, inside2
928 : integer :: i,j,ii
929 : type (cartesian2D_t) :: corners(4),sphere_xy,cart
930 : type (cartesian3D_t) :: corners_xyz(4),center,a,b,cross(4)
931 : real (kind=r8) :: yp(4), y, elem_diam,dotprod
932 : real (kind=r8) :: xp(4), x
933 : real (kind=r8) :: d1,d2, tol_inside = 1.0e-12_r8
934 :
935 : type (spherical_polar_t) :: sphere ! debug
936 :
937 0 : inside = .false.
938 :
939 : ! first check if point is near the element:
940 0 : corners_xyz(:) = elem%corners3D(:)
941 : elem_diam = max( distance(corners_xyz(1),corners_xyz(3)), &
942 0 : distance(corners_xyz(2),corners_xyz(4)) )
943 :
944 0 : center%x = sum(corners_xyz(1:4)%x)/4
945 0 : center%y = sum(corners_xyz(1:4)%y)/4
946 0 : center%z = sum(corners_xyz(1:4)%z)/4
947 0 : if ( distance(center,sphere_xyz) > 1.0_r8*elem_diam ) return
948 :
949 : j = 4
950 0 : do i=1,4
951 : ! outward normal to plane containing j->i edge: corner(i) x corner(j)
952 : ! sphere dot (corner(i) x corner(j) ) = negative if inside
953 0 : cross(i)%x = corners_xyz(i)%y*corners_xyz(j)%z - corners_xyz(i)%z*corners_xyz(j)%y
954 0 : cross(i)%y =-(corners_xyz(i)%x*corners_xyz(j)%z - corners_xyz(i)%z*corners_xyz(j)%x)
955 0 : cross(i)%z = corners_xyz(i)%x*corners_xyz(j)%y - corners_xyz(i)%y*corners_xyz(j)%x
956 : dotprod = cross(i)%x*sphere_xyz%x + cross(i)%y*sphere_xyz%y +&
957 0 : cross(i)%z*sphere_xyz%z
958 0 : j = i ! within this loopk j = i-1
959 :
960 : !if (dotprod>0 .and. dotprod/elem_diam < 1e-5) print *,dotprod/elem_diam
961 :
962 : ! dot product is proportional to elem_diam. positive means outside,
963 : ! but allow machine precision tolorence:
964 0 : if (dotprod > tol_inside*elem_diam) return
965 : !if (dotprod > 0) return
966 : end do
967 0 : inside=.true.
968 : return
969 : end function point_inside_gc
970 :
971 :
972 : !================================================
973 : ! (Nair) Cube face index and local coordinates
974 : !================================================
975 :
976 0 : subroutine cube_facepoint_ne(sphere, ne, cart, number)
977 : use coordinate_systems_mod, only : cube_face_number_from_sphere, sphere2cubedsphere
978 :
979 : type(spherical_polar_t), intent(in) :: sphere
980 : integer, intent(in) :: ne
981 : type(cartesian2D_t), intent(out) :: cart
982 : integer, intent(out) :: number
983 :
984 : real(kind=r8) :: xp, yp
985 : type(cartesian2D_t) :: cube
986 : integer :: ie, je, face_no
987 : real(kind=r8) :: x1, x2
988 : real(kind=r8) :: dx
989 :
990 0 : face_no = cube_face_number_from_sphere(sphere)
991 0 : cube = sphere2cubedsphere(sphere, face_no)
992 0 : xp = cube%x
993 0 : yp = cube%y
994 :
995 : ! MNL: for uniform grids (on cube face), analytic solution is fine
996 0 : x1 = xp + 0.25_r8*PI
997 0 : x2 = yp + 0.25_r8*PI
998 :
999 0 : dx = (0.5_r8*PI)/ne
1000 0 : ie = INT(ABS(x1)/dx)
1001 0 : je = INT(ABS(x2)/dx)
1002 : ! if we are exactly on an element edge, we can put the point in
1003 : ! either the ie or ie+1 element, EXCEPT if ie==ne.
1004 0 : if ( ABS(x1) < ne*dx ) then
1005 0 : ie = ie + 1
1006 : end if
1007 0 : if ( ABS(x2) < ne*dx ) then
1008 0 : je = je + 1
1009 : end if
1010 0 : if ((ie > ne) .or. (je > ne)) then
1011 0 : write(iulog, *) 'ERROR: ',ie,je,ne
1012 0 : write(iulog, *) 'lat,lon=',sphere%lat,sphere%lon
1013 0 : write(iulog, *) 'face no=',face_no
1014 0 : write(iulog, *) x1,x2,x1/dx,x2/dx
1015 0 : call endrun('interpolate_mod: bad argument')
1016 : endif
1017 :
1018 : ! bug fix MT 1/2009. This was creating a plotting error at
1019 : ! the row of elements in iface=2 at 50 degrees (NE=16 128x256 lat/lon grid)
1020 : ! For point on element edge, we can have ie=2, but x1=dx
1021 : ! but if ie>1, we must execute this statement.
1022 : ! The only time we can skip this statement is if ie=1, but then
1023 : ! the statement has no effect, so lets never skip it:
1024 : ! if (x1 > dx ) then
1025 0 : x1 = x1 - dble(ie-1)*dx
1026 : ! endif
1027 :
1028 0 : x1 = 2.0_r8*(x1/dx)-1.0_r8
1029 :
1030 : ! if (x2 > dx ) then ! removed MT 1/2009, see above
1031 0 : x2 = x2 - dble(je-1)*dx
1032 : ! endif
1033 :
1034 0 : x2 = 2.0_r8*(x2/dx)-1.0_r8
1035 :
1036 : ! coordinates within an element [-1,1]
1037 0 : cart%x = x1
1038 0 : cart%y = x2
1039 0 : number = ie + (je-1)*ne + (face_no-1)*ne*ne
1040 0 : end subroutine cube_facepoint_ne
1041 : !================================================
1042 : ! (Nair) Cube face index and local coordinates
1043 : !================================================
1044 :
1045 :
1046 0 : subroutine cube_facepoint_unstructured(sphere,cart, number, elem)
1047 : use coordinate_systems_mod, only : cube_face_number_from_sphere, &
1048 : sphere2cubedsphere,change_coordinates,cube_face_number_from_cart
1049 : implicit none
1050 :
1051 : type (element_t) , intent(in), target :: elem(:)
1052 : type (spherical_polar_t), intent (in) :: sphere
1053 : type (cartesian2D_t), intent(out) :: cart
1054 : integer , intent(out) :: number
1055 :
1056 : integer :: ii
1057 : Logical :: found
1058 : type (cartesian3D_t) :: sphere_xyz
1059 : type (cartesian2D_t) :: cube
1060 0 : sphere_xyz=spherical_to_cart(sphere)
1061 :
1062 0 : number=-1
1063 : ! print *,'WARNING: using GC map'
1064 0 : do ii = 1,nelemd
1065 : ! for equiangular gnomonic map:
1066 : ! unstructed grid element edges are NOT great circles
1067 0 : if (cubed_sphere_map==0) then
1068 0 : found = point_inside_equiangular(elem(ii), sphere, sphere_xyz)
1069 : else
1070 : ! assume element edges are great circle arcs:
1071 0 : found = point_inside_gc(elem(ii), sphere_xyz)
1072 : endif
1073 :
1074 0 : if (found) then
1075 0 : number = ii
1076 0 : cart = parametric_coordinates(sphere, elem(ii)%corners3D,&
1077 0 : cubed_sphere_map,elem(ii)%corners,elem(ii)%u2qmap,elem(ii)%facenum)
1078 0 : exit
1079 : end if
1080 : end do
1081 0 : end subroutine cube_facepoint_unstructured
1082 :
1083 :
1084 0 : subroutine interp_init()
1085 : type (quadrature_t) :: gp
1086 :
1087 0 : gp = gausslobatto(np)
1088 0 : call interpolate_create(gp,interp_p)
1089 0 : end subroutine interp_init
1090 :
1091 :
1092 0 : subroutine setup_latlon_interp(elem,interpdata,par)
1093 : !
1094 : ! initialize interpolation data structures to interpolate to a lat-lon grid
1095 : !
1096 : !
1097 :
1098 : implicit none
1099 : type (element_t) , intent(in), target :: elem(:)
1100 : type (parallel_t) , intent(in) :: par
1101 : type (interpdata_t) , intent(out) :: interpdata(:)
1102 :
1103 : ! local
1104 : integer i,j,ii,count_total,n_interp,count_max
1105 : integer ngrid, number, elem_num, plat
1106 : integer countx, missing_pts,ierr
1107 : integer :: npts_mult_claims,max_claims
1108 :
1109 0 : real (kind=r8) :: dp,latdeg(nlat+1),clat(nlat+1),w(nlat+1),w_staggered(nlat)
1110 0 : real (kind=r8) :: clat_staggered(nlat),latdeg_st(nlat),err,err2
1111 :
1112 : type (spherical_polar_t) :: sphere
1113 : type (cartesian2D_t) :: cart
1114 : type (cartesian3D_t) :: sphere_xyz,sphere2_xyz
1115 :
1116 : type (quadrature_t) :: gp
1117 :
1118 :
1119 : ! Array to make sure each interp point is on exactly one process
1120 0 : type (cartesian2D_t),allocatable :: cart_vec(:,:)
1121 : integer :: k
1122 0 : integer, allocatable :: global_elem_gid(:,:),local_elem_gid(:,:), local_elem_num(:,:)
1123 :
1124 : ! these arrays often are too large for stack, so lets make sure
1125 : ! they go on the heap:
1126 0 : allocate(local_elem_num(nlat,nlon))
1127 0 : allocate(local_elem_gid(nlat,nlon))
1128 0 : allocate(global_elem_gid(nlat,nlon))
1129 0 : allocate(cart_vec(nlat,nlon))
1130 :
1131 0 : if (par%masterproc) then
1132 0 : write(iulog,'(a,i4,a,i4,a)') 'Initializing ',nlat,' x ',nlon,' lat-lon interpolation grid: '
1133 : endif
1134 :
1135 0 : do ii=1,nelemd
1136 0 : interpdata(ii)%n_interp=0 ! reset counter
1137 : enddo
1138 :
1139 0 : if (associated(lat))then
1140 0 : deallocate(lat)
1141 : nullify(lat)
1142 : endif
1143 0 : if (associated(gweight))then
1144 0 : deallocate(gweight)
1145 : nullify(gweight)
1146 : endif
1147 :
1148 0 : if (associated(lon))then
1149 0 : deallocate(lon)
1150 : nullify(lon)
1151 : endif
1152 :
1153 0 : allocate(lat(nlat))
1154 0 : allocate(gweight(nlat))
1155 0 : allocate(lon(nlon))
1156 0 : call interp_init()
1157 0 : gweight=0
1158 0 : do i=1,nlon
1159 0 : lon(i)=2*pi*(i-1)/nlon
1160 : enddo
1161 0 : if (gridtype==1) then
1162 0 : do j=1,nlat
1163 0 : lat(j) = -pi/2 + pi*(j-1)/(nlat-1)
1164 : end do
1165 : plat=nlat
1166 : endif
1167 0 : if (gridtype==2) then
1168 0 : gp=gauss(nlat)
1169 0 : do j=1,nlat
1170 0 : lat(j) = asin(gp%points(j))
1171 0 : gweight(j) = gp%weights(j)
1172 : end do
1173 : endif
1174 0 : if (gridtype==3) then
1175 0 : do j=1,nlat
1176 0 : lat(j) = -pi/2 + pi*(j-.5_r8)/nlat
1177 : end do
1178 0 : plat=nlat+1
1179 : endif
1180 :
1181 0 : if (gridtype==1 .or. gridtype==3) then
1182 : ! gridtype=1 plat=nlat gweight(1:nlat)=w(1:plat)
1183 : ! gridtype=3 plat=nlat+1 gweight(1:nlat)=w_staggered(1:plat-1)
1184 :
1185 : ! L-R dynamics uses a regular latitude distribution (not gausian).
1186 : ! The algorithm below is a bastardized version of LSM: map.F.
1187 0 : dp = 180.0_r8/(plat-1)
1188 0 : do j = 1, plat
1189 0 : latdeg(j) = -90.0_r8 + (j-1)*dp
1190 0 : clat(j) = latdeg(j)*pi/180.0_r8
1191 : end do
1192 :
1193 : ! Calculate latitudes for the staggered grid
1194 :
1195 0 : do j = 1, plat-1
1196 0 : clat_staggered(j) = (clat(j) + clat(j+1)) / 2
1197 0 : latdeg_st (j) = clat_staggered(j)*180.0_r8/pi
1198 : end do
1199 :
1200 : ! Weights are defined as cos(phi)*(delta-phi)
1201 : ! For a sanity check, the sum of w across all lats should be 2, or 1 across
1202 : ! half of the latitudes.
1203 :
1204 0 : do j = 2, plat-1
1205 0 : w(j) = sin(clat_staggered(j)) - sin(clat_staggered(j-1))
1206 : end do
1207 0 : w(1) = sin(clat_staggered(1)) + 1
1208 0 : w(plat) = w(1)
1209 :
1210 : ! with nlat=2048, this error was 4e-16
1211 0 : if (abs(sum(w(1:plat)) - 2) > 1.0e-8_r8) then
1212 0 : write(iulog,*) 'interpolate_mod: w weights do not sum to 2. sum=',sum(w(1:plat))
1213 0 : call endrun('interpolate_mod: weights do not sum to 2.')
1214 : end if
1215 :
1216 0 : dp = pi / (plat-1)
1217 0 : do j = 1, plat-1
1218 0 : w_staggered(j) = sin(clat(j+1)) - sin(clat(j))
1219 : end do
1220 :
1221 :
1222 0 : if (abs(sum(w_staggered(1:plat-1)) - 2) > 1.0e-8_r8) then
1223 0 : write(iulog,*) 'interpolate_mod: staggered weights do not sum to 2. sum=',sum(w_staggered(1:plat-1))
1224 0 : call endrun('interpolate_mod: weights do not sum to 2.')
1225 : end if
1226 :
1227 0 : if (gridtype==1) then
1228 0 : gweight(1:nlat)=w(1:plat)
1229 : endif
1230 0 : if (gridtype==3) then
1231 0 : gweight(1:nlat)=w_staggered(1:plat-1)
1232 : endif
1233 : endif
1234 :
1235 :
1236 : ! go through once, counting the number of points on each element
1237 0 : sphere%r=1
1238 0 : local_elem_num = -1
1239 0 : local_elem_gid = -1
1240 0 : global_elem_gid = -1
1241 0 : err=0
1242 0 : do j=1,nlat
1243 0 : do i=1,nlon
1244 0 : sphere%lat=lat(j)
1245 0 : sphere%lon=lon(i)
1246 :
1247 0 : number = -1
1248 0 : if ( (cubed_sphere_map /= 0) .or. MeshUseMeshFile) then
1249 0 : call cube_facepoint_unstructured(sphere, cart, number, elem)
1250 0 : if (number /= -1) then
1251 : ! If points are outside element but within tolerance, move to boundary
1252 0 : if (cart%x + 1.0_r8.le.0.0_r8) cart%x = -1.0_r8
1253 0 : if (cart%x - 1.0_r8.ge.0.0_r8) cart%x = 1.0_r8
1254 0 : if (cart%y + 1.0_r8.le.0.0_r8) cart%y = -1.0_r8
1255 0 : if (cart%y - 1.0_r8.ge.0.0_r8) cart%y = 1.0_r8
1256 :
1257 0 : local_elem_num(j,i) = number
1258 0 : local_elem_gid(j,i) = elem(number)%vertex%number
1259 0 : cart_vec(j,i) = cart ! local element coordiante of interpolation point
1260 : endif
1261 : else
1262 0 : call cube_facepoint_ne(sphere, ne, cart, number)
1263 : ! the sphere point belongs to the element number on face = face_no.
1264 : ! do I own this element?
1265 0 : if (number /= -1) then
1266 0 : do ii=1,nelemd
1267 0 : if (number == elem(ii)%vertex%number) then
1268 0 : local_elem_gid(j,i) = number
1269 0 : local_elem_num(j,i) = ii
1270 0 : cart_vec(j,i) = cart ! local element coordinate found above
1271 0 : exit
1272 : endif
1273 : enddo
1274 : endif
1275 : endif
1276 0 : ii=local_elem_num(j,i)
1277 0 : if (ii /= -1) then
1278 : ! compute error: map 'cart' back to sphere and compare with original
1279 : ! interpolation point:
1280 : sphere2_xyz = spherical_to_cart( ref2sphere(cart%x,cart%y, &
1281 0 : elem(ii)%corners3D,cubed_sphere_map,elem(ii)%corners,elem(ii)%facenum ))
1282 0 : sphere_xyz = spherical_to_cart(sphere)
1283 0 : err=max(err,distance(sphere2_xyz,sphere_xyz))
1284 : endif
1285 : enddo
1286 0 : if (par%masterproc) then
1287 0 : if ((MOD(j,64).eq.1).or.(j.eq.nlat)) then
1288 0 : print *,'finished latitude ',j,' of ',nlat
1289 : endif
1290 : endif
1291 : enddo
1292 0 : err2=err
1293 0 : call MPI_Allreduce(err, err2, 1, MPI_real8, MPI_MAX, par%comm, ierr)
1294 0 : if (par%masterproc) then
1295 0 : write(iulog,'(a,e12.4)') 'Max interpolation point search error: ',err2
1296 : endif
1297 :
1298 : ! if multile elements claim a interpolation point, take the one with largest gid:
1299 0 : global_elem_gid = local_elem_gid
1300 0 : call MPI_Allreduce(local_elem_gid, global_elem_gid, nlat*nlon, MPI_integer, MPI_MAX, par%comm,ierr)
1301 :
1302 0 : missing_pts=0
1303 0 : do j=1,nlat
1304 0 : do i=1,nlon
1305 0 : if (global_elem_gid(j,i) == -1 ) then
1306 0 : missing_pts = missing_pts + 1
1307 0 : if (par%masterproc) &
1308 0 : print *,'Error: point not claimed by any element j,i,lat(j),lon(i)=',j,i,lat(j),lon(i)
1309 0 : else if (local_elem_gid(j,i) == global_elem_gid(j,i) ) then
1310 0 : ii = local_elem_num(j,i)
1311 0 : interpdata(ii)%n_interp = interpdata(ii)%n_interp + 1
1312 : endif
1313 : end do
1314 : end do
1315 :
1316 0 : countx=maxval(interpdata(1:nelemd)%n_interp)
1317 0 : count_max = countx
1318 0 : call MPI_Allreduce(countx,count_max,1,MPI_integer,MPI_MAX,par%comm,ierr)
1319 :
1320 0 : if (par%masterproc) then
1321 0 : write(iulog,'(a,i6)') 'Maximum number of interpolation points claimed by an element: ',count_max
1322 : endif
1323 :
1324 : ! allocate storage
1325 0 : do ii=1,nelemd
1326 0 : ngrid = interpdata(ii)%n_interp
1327 0 : if(interpdata(ii)%first_entry)then
1328 0 : NULLIFY(interpdata(ii)%interp_xy)
1329 0 : NULLIFY(interpdata(ii)%ilat)
1330 0 : NULLIFY(interpdata(ii)%ilon)
1331 :
1332 0 : interpdata(ii)%first_entry=.FALSE.
1333 : endif
1334 0 : if(associated(interpdata(ii)%interp_xy))then
1335 0 : if(size(interpdata(ii)%interp_xy)>0)deallocate(interpdata(ii)%interp_xy)
1336 : endif
1337 0 : if(associated(interpdata(ii)%ilat))then
1338 0 : if(size(interpdata(ii)%ilat)>0)deallocate(interpdata(ii)%ilat)
1339 : endif
1340 :
1341 0 : if (associated(interpdata(ii)%ilon))then
1342 0 : if(size(interpdata(ii)%ilon)>0)deallocate(interpdata(ii)%ilon)
1343 : endif
1344 0 : allocate(interpdata(ii)%interp_xy( ngrid ) )
1345 0 : allocate(interpdata(ii)%ilat( ngrid ) )
1346 0 : allocate(interpdata(ii)%ilon( ngrid ) )
1347 0 : interpdata(ii)%n_interp=0 ! reset counter
1348 : enddo
1349 0 : do j=1,nlat
1350 0 : do i=1,nlon
1351 0 : if (local_elem_gid(j,i) == global_elem_gid(j,i) .and. &
1352 0 : local_elem_gid(j,i) /= -1 ) then
1353 0 : ii = local_elem_num(j,i)
1354 0 : ngrid = interpdata(ii)%n_interp + 1
1355 0 : interpdata(ii)%n_interp = ngrid
1356 0 : interpdata(ii)%interp_xy( ngrid ) = cart_vec(j,i)
1357 0 : interpdata(ii)%ilon( ngrid ) = i
1358 0 : interpdata(ii)%ilat( ngrid ) = j
1359 : endif
1360 : enddo
1361 : enddo
1362 :
1363 : ! now lets compute the number of points that were claimed by
1364 : ! more than one element:
1365 0 : do j=1,nlat
1366 0 : do i=1,nlon
1367 0 : if (local_elem_gid(j,i) == -1) then
1368 0 : local_elem_gid(j,i)=0
1369 : else
1370 0 : local_elem_gid(j,i)=1
1371 : endif
1372 : enddo
1373 : enddo
1374 0 : global_elem_gid = local_elem_gid
1375 0 : call MPI_Allreduce(local_elem_gid, global_elem_gid, nlat*nlon, MPI_integer, MPI_SUM, par%comm,ierr)
1376 0 : if (par%masterproc) then
1377 0 : countx=0
1378 0 : do j=1,nlat
1379 0 : do i=1,nlon
1380 0 : if (global_elem_gid(j,i)>1) countx=countx+1
1381 : enddo
1382 : enddo
1383 0 : npts_mult_claims=countx
1384 0 : max_claims=maxval(global_elem_gid)
1385 : endif
1386 :
1387 0 : if (par%masterproc) then
1388 0 : print *,'Number of interpolation points claimed by more than one element: ',npts_mult_claims
1389 0 : print *,'max number of elements which claimed the same interpolation point:',max_claims
1390 : endif
1391 :
1392 0 : deallocate(global_elem_gid)
1393 0 : deallocate(local_elem_num)
1394 0 : deallocate(local_elem_gid)
1395 0 : deallocate(cart_vec)
1396 :
1397 : ! check if every point in interpolation grid was claimed by an element:
1398 0 : if (missing_pts>0) then
1399 0 : count_total = nlat*nlon
1400 0 : if(par%masterproc) then
1401 0 : write(iulog,"(3A,I4,A,I7,a,i5)")"Error:",__FILE__," ",__LINE__," count_total:",count_total," missing:",missing_pts
1402 : end if
1403 0 : call syncmp(par)
1404 0 : call endrun('Error: interpolation points not claimed by any element')
1405 : endif
1406 :
1407 :
1408 0 : end subroutine setup_latlon_interp
1409 :
1410 :
1411 :
1412 : ! interpolate_scalar
1413 : !
1414 : ! Interpolate a scalar field given in an element (fld_cube) to the points in
1415 : ! interpdata%interp_xy(i), i=1 .. interpdata%n_interp.
1416 : !
1417 : ! Note that it is possible the given element contains none of the interpolation points
1418 : ! =======================================
1419 0 : subroutine interpolate_ce(cart,fld_cube,npts,fld, fillvalue)
1420 : type (cartesian2D_t) :: cart
1421 : integer :: npts
1422 : real (kind=r8) :: fld_cube(npts,npts) ! cube field
1423 : real (kind=r8) :: fld ! field at new grid lat,lon coordinates
1424 : real (kind=r8), intent(in), optional :: fillvalue
1425 : ! Local variables
1426 : type (interpolate_t), pointer :: interp ! interpolation structure
1427 :
1428 : integer :: ne
1429 : integer :: i
1430 :
1431 0 : if (npts==np) then
1432 : interp => interp_p
1433 : else
1434 0 : call endrun('Error in interpolate_scalar(): must be called with p or v grid data')
1435 : endif
1436 :
1437 0 : fld=interpolate_2d(cart,fld_cube,interp,npts,fillvalue)
1438 :
1439 0 : end subroutine interpolate_ce
1440 :
1441 :
1442 :
1443 : ! =======================================
1444 : ! interpolate_scalar
1445 : !
1446 : ! Interpolate a scalar field given in an element (fld_cube) to the points in
1447 : ! interpdata%interp_xy(i), i=1 .. interpdata%n_interp.
1448 : !
1449 : ! Note that it is possible the given element contains none of the interpolation points
1450 : ! =======================================
1451 0 : subroutine interpolate_scalar2d(interpdata,fld_cube,nsize,nhalo,fld, fillvalue)
1452 : use dimensions_mod, only: npsq, fv_nphys,nc
1453 : integer, intent(in) :: nsize,nhalo
1454 : real (kind=r8), intent(in) :: fld_cube(1-nhalo:nsize+nhalo,1-nhalo:nsize+nhalo) ! cube field
1455 : real (kind=r8), intent(out):: fld(:) ! field at new grid lat,lon coordinates
1456 : type (interpdata_t), intent(in) :: interpdata
1457 : real (kind=r8), intent(in), optional :: fillvalue
1458 : ! Local variables
1459 : type (interpolate_t), pointer :: interp ! interpolation structure
1460 :
1461 : integer :: i,imin,imax,ne
1462 0 : real (kind=r8):: xoy(1-nhalo:nsize+nhalo),dx
1463 : type (cartesian2D_t) :: cart
1464 :
1465 0 : if (nsize==np.and.nhalo==0) then
1466 : !
1467 : ! GLL grid
1468 : !
1469 0 : interp => interp_p
1470 0 : xoy = interp%glp(:)
1471 0 : imin = 1
1472 0 : imax = np
1473 0 : else if (nhalo>0.and.(nsize==fv_nphys.or.nsize==nc)) then
1474 : !
1475 : ! finite-volume grid
1476 : !
1477 0 : if (itype.ne.1) then
1478 0 : call endrun('itype must be 1 for latlon output from finite-volume (non-GLL) grids')
1479 : end if
1480 0 : imin = 1-nhalo
1481 0 : imax = nsize+nhalo
1482 : !
1483 : ! create normalized coordinates
1484 : !
1485 0 : dx = 2.0_r8/REAL(nsize,KIND=r8)
1486 0 : do i=imin,imax
1487 0 : xoy(i) = -1.0_r8+(i-0.5_r8)*dx
1488 : end do
1489 : else
1490 0 : call endrun('interpolate_scalar2d: resolution not supported')
1491 : endif
1492 :
1493 : ! Choice for Native (high-order) or Bilinear interpolations
1494 0 : if (itype == 0) then
1495 0 : do i=1,interpdata%n_interp
1496 0 : fld(i)=interpolate_2d(interpdata%interp_xy(i),fld_cube,interp,nsize,fillvalue)
1497 : end do
1498 0 : else if (itype == 1) then
1499 0 : do i=1,interpdata%n_interp
1500 0 : fld(i)=interpol_bilinear(interpdata%interp_xy(i),fld_cube,xoy,imin,imax,fillvalue)
1501 : end do
1502 : else
1503 0 : call endrun("interpolate_scalar2d: wrong interpolation type: "//int2str(itype))
1504 : end if
1505 :
1506 0 : end subroutine interpolate_scalar2d
1507 0 : subroutine interpolate_scalar3d(interpdata,fld_cube,nsize,nhalo,nlev,fld, fillvalue)
1508 : use dimensions_mod, only: npsq, fv_nphys,nc
1509 : integer , intent(in) :: nsize, nhalo, nlev
1510 : real (kind=r8), intent(in) :: fld_cube(1-nhalo:nsize+nhalo,1-nhalo:nsize+nhalo,nlev) ! cube field
1511 : real (kind=r8), intent(out) :: fld(:,:) ! field at new grid lat,lon coordinates
1512 : type (interpdata_t), intent(in) :: interpdata
1513 : real (kind=r8), intent(in), optional :: fillvalue
1514 : ! Local variables
1515 : type (interpolate_t), pointer :: interp ! interpolation structure
1516 :
1517 : integer :: ne
1518 :
1519 : integer :: i, k, imin, imax
1520 0 : real (kind=r8) :: xoy(1-nhalo:nsize+nhalo),dx
1521 :
1522 : type (cartesian2D_t) :: cart
1523 :
1524 0 : if (nsize==np.and.nhalo==0) then
1525 : !
1526 : ! GLL grid
1527 : !
1528 0 : interp => interp_p
1529 0 : xoy = interp%glp(:)
1530 0 : imin = 1
1531 0 : imax = np
1532 0 : else if (nhalo>0.and.(nsize==fv_nphys.or.nsize==nc)) then
1533 : !
1534 : ! finite-volume grid
1535 : !
1536 0 : if (itype.ne.1) then
1537 0 : call endrun('itype must be 1 for latlon output from finite-volume (non-GLL) grids')
1538 : end if
1539 0 : imin = 1-nhalo
1540 0 : imax = nsize+nhalo
1541 : !
1542 : ! create normalized coordinates
1543 : !
1544 0 : dx = 2.0_r8/REAL(nsize,KIND=r8)
1545 0 : do i=imin,imax
1546 0 : xoy(i) = -1.0_r8+(i-0.5_r8)*dx
1547 : end do
1548 : else
1549 0 : call endrun('interpolate_scalar3d: resolution not supported')
1550 : endif
1551 :
1552 : ! Choice for Native (high-order) or Bilinear interpolations
1553 0 : if (itype == 0) then
1554 0 : do k=1,nlev
1555 0 : do i=1,interpdata%n_interp
1556 0 : fld(i,k)=interpolate_2d(interpdata%interp_xy(i),fld_cube(:,:,k),interp,nsize,fillvalue)
1557 : end do
1558 : end do
1559 0 : elseif (itype == 1) then
1560 0 : do k=1,nlev
1561 0 : do i=1,interpdata%n_interp
1562 0 : fld(i,k)=interpol_bilinear(interpdata%interp_xy(i),fld_cube(:,:,k),xoy,imin,imax,fillvalue)
1563 : end do
1564 : end do
1565 : else
1566 0 : call endrun("interpolate_scalar3d: wrong interpolation type: "//int2str(itype))
1567 : endif
1568 0 : end subroutine interpolate_scalar3d
1569 :
1570 :
1571 : ! =======================================
1572 : ! interpolate_vector
1573 : !
1574 : ! Interpolate a vector field given in an element (fld_cube)
1575 : ! to the points in interpdata%interp_xy(i), i=1 .. interpdata%n_interp.
1576 : !
1577 : ! input_coords = 0 fld_cube given in lat-lon
1578 : ! input_coords = 1 fld_cube given in contravariant
1579 : !
1580 : ! Note that it is possible the given element contains none of the interpolation points
1581 : ! =======================================
1582 0 : subroutine interpolate_vector2d(interpdata,elem,fld_cube,npts,fld,input_coords, fillvalue)
1583 : implicit none
1584 : integer :: npts
1585 : real (kind=r8) :: fld_cube(npts,npts,2) ! vector field
1586 : real (kind=r8) :: fld(:,:) ! field at new grid lat,lon coordinates
1587 : type (interpdata_t) :: interpdata
1588 : type (element_t), intent(in) :: elem
1589 : real (kind=r8), intent(in), optional :: fillvalue
1590 : integer :: input_coords
1591 :
1592 :
1593 : ! Local variables
1594 0 : real (kind=r8) :: fld_contra(npts,npts,2) ! vector field
1595 : type (interpolate_t), pointer :: interp ! interpolation structure
1596 :
1597 : real (kind=r8) :: v1,v2
1598 : real (kind=r8) :: D(2,2) ! derivative of gnomonic mapping
1599 : real (kind=r8) :: JJ(2,2), tmpD(2,2) ! derivative of gnomonic mapping
1600 :
1601 : integer :: i,j
1602 :
1603 : type (cartesian2D_t) :: cart
1604 :
1605 0 : if(present(fillvalue)) then
1606 0 : if (any(fld_cube==fillvalue)) then
1607 0 : fld = fillvalue
1608 : return
1609 : end if
1610 : end if
1611 :
1612 0 : if (input_coords==0 ) then
1613 : ! convert to contra
1614 0 : do j=1,npts
1615 0 : do i=1,npts
1616 : ! latlon->contra
1617 0 : fld_contra(i,j,1) = elem%Dinv(i,j,1,1)*fld_cube(i,j,1) + elem%Dinv(i,j,1,2)*fld_cube(i,j,2)
1618 0 : fld_contra(i,j,2) = elem%Dinv(i,j,2,1)*fld_cube(i,j,1) + elem%Dinv(i,j,2,2)*fld_cube(i,j,2)
1619 : enddo
1620 : enddo
1621 : else
1622 0 : fld_contra=fld_cube
1623 : endif
1624 :
1625 :
1626 0 : if (npts==np) then
1627 : interp => interp_p
1628 : else
1629 0 : call endrun('interpolate_vector2d: Error in interpolate_vector(): input must be on GLL grid')
1630 : endif
1631 :
1632 :
1633 : ! Choice for Native (high-order) or Bilinear interpolations
1634 :
1635 0 : if (itype == 0) then
1636 0 : do i=1,interpdata%n_interp
1637 0 : fld(i,1)=interpolate_2d(interpdata%interp_xy(i),fld_contra(:,:,1),interp,npts)
1638 0 : fld(i,2)=interpolate_2d(interpdata%interp_xy(i),fld_contra(:,:,2),interp,npts)
1639 : end do
1640 0 : elseif (itype == 1) then
1641 0 : do i=1,interpdata%n_interp
1642 0 : fld(i,1)=interpol_bilinear(interpdata%interp_xy(i),fld_contra(:,:,1),interp%glp(:),1,np)
1643 0 : fld(i,2)=interpol_bilinear(interpdata%interp_xy(i),fld_contra(:,:,2),interp%glp(:),1,np)
1644 : end do
1645 : else
1646 0 : call endrun("interpolate_vector2d: wrong interpolation type: "//int2str(itype))
1647 : endif
1648 0 : do i=1,interpdata%n_interp
1649 : ! convert fld from contra->latlon
1650 0 : call dmap(D,interpdata%interp_xy(i)%x,interpdata%interp_xy(i)%y,&
1651 0 : elem%corners3D,cubed_sphere_map,elem%corners,elem%u2qmap,elem%facenum)
1652 : ! convert fld from contra->latlon
1653 0 : v1 = fld(i,1)
1654 0 : v2 = fld(i,2)
1655 :
1656 0 : fld(i,1)=D(1,1)*v1 + D(1,2)*v2
1657 0 : fld(i,2)=D(2,1)*v1 + D(2,2)*v2
1658 : end do
1659 :
1660 : end subroutine interpolate_vector2d
1661 :
1662 : ! =======================================
1663 : ! interpolate_vector
1664 : !
1665 : ! Interpolate a vector field given in an element (fld_cube)
1666 : ! to the points in interpdata%interp_xy(i), i=1 .. interpdata%n_interp.
1667 : !
1668 : ! input_coords = 0 fld_cube given in lat-lon
1669 : ! input_coords = 1 fld_cube given in contravariant
1670 : !
1671 : ! Note that it is possible the given element contains none of the interpolation points
1672 : ! =======================================
1673 0 : subroutine interpolate_vector3d(interpdata,elem,fld_cube,npts,nlev,fld,input_coords, fillvalue)
1674 : implicit none
1675 : type (interpdata_t),intent(in) :: interpdata
1676 : type (element_t), intent(in) :: elem
1677 : integer, intent(in) :: npts, nlev
1678 : real (kind=r8), intent(in) :: fld_cube(npts,npts,2,nlev) ! vector field
1679 : real (kind=r8), intent(out) :: fld(:,:,:) ! field at new grid lat,lon coordinates
1680 : real (kind=r8), intent(in),optional :: fillvalue
1681 : integer, intent(in) :: input_coords
1682 :
1683 : ! Local variables
1684 0 : real (kind=r8) :: fld_contra(npts,npts,2,nlev) ! vector field
1685 : type (interpolate_t), pointer :: interp ! interpolation structure
1686 :
1687 : real (kind=r8) :: v1,v2
1688 : real (kind=r8) :: D(2,2) ! derivative of gnomonic mapping
1689 : real (kind=r8) :: JJ(2,2), tmpD(2,2) ! derivative of gnomonic mapping
1690 :
1691 :
1692 : integer :: i,j,k
1693 :
1694 : type (cartesian2D_t) :: cart
1695 0 : if(present(fillvalue)) then
1696 0 : if (any(fld_cube==fillvalue)) then
1697 0 : fld = fillvalue
1698 : return
1699 : end if
1700 : end if
1701 0 : if (input_coords==0 ) then
1702 : ! convert to contra
1703 0 : do k=1,nlev
1704 0 : do j=1,npts
1705 0 : do i=1,npts
1706 : ! latlon->contra
1707 0 : fld_contra(i,j,1,k) = elem%Dinv(i,j,1,1)*fld_cube(i,j,1,k) + elem%Dinv(i,j,1,2)*fld_cube(i,j,2,k)
1708 0 : fld_contra(i,j,2,k) = elem%Dinv(i,j,2,1)*fld_cube(i,j,1,k) + elem%Dinv(i,j,2,2)*fld_cube(i,j,2,k)
1709 : enddo
1710 : enddo
1711 : end do
1712 : else
1713 0 : fld_contra=fld_cube
1714 : endif
1715 :
1716 0 : if (npts==np) then
1717 : interp => interp_p
1718 : else
1719 0 : call endrun('interpolate_vector3d: Error in interpolate_vector(): input must be on GLL grid')
1720 : endif
1721 :
1722 :
1723 : ! Choice for Native (high-order) or Bilinear interpolations
1724 :
1725 0 : if (itype == 0) then
1726 0 : do k=1,nlev
1727 0 : do i=1,interpdata%n_interp
1728 0 : fld(i,k,1)=interpolate_2d(interpdata%interp_xy(i),fld_contra(:,:,1,k),interp,npts)
1729 0 : fld(i,k,2)=interpolate_2d(interpdata%interp_xy(i),fld_contra(:,:,2,k),interp,npts)
1730 : end do
1731 : end do
1732 0 : elseif (itype == 1) then
1733 0 : do k=1,nlev
1734 0 : do i=1,interpdata%n_interp
1735 0 : fld(i,k,1)=interpol_bilinear(interpdata%interp_xy(i),fld_contra(:,:,1,k),interp%glp(:),1,np)
1736 0 : fld(i,k,2)=interpol_bilinear(interpdata%interp_xy(i),fld_contra(:,:,2,k),interp%glp(:),1,np)
1737 : end do
1738 : end do
1739 : else
1740 0 : call endrun("interpolate_vector3d: wrong interpolation type: "//int2str(itype))
1741 : endif
1742 :
1743 :
1744 0 : do i=1,interpdata%n_interp
1745 : ! compute D(:,:) at the point elem%interp_cube(i)
1746 0 : call dmap(D,interpdata%interp_xy(i)%x,interpdata%interp_xy(i)%y,&
1747 0 : elem%corners3D,cubed_sphere_map,elem%corners,elem%u2qmap,elem%facenum)
1748 0 : do k=1,nlev
1749 : ! convert fld from contra->latlon
1750 0 : v1 = fld(i,k,1)
1751 0 : v2 = fld(i,k,2)
1752 :
1753 0 : fld(i,k,1)=D(1,1)*v1 + D(1,2)*v2
1754 0 : fld(i,k,2)=D(2,1)*v1 + D(2,2)*v2
1755 : end do
1756 : end do
1757 : end subroutine interpolate_vector3d
1758 :
1759 0 : subroutine vec_latlon_to_contra(elem,nphys,nhalo,nlev,fld,fvm)
1760 : use fvm_control_volume_mod, only: fvm_struct
1761 : use dimensions_mod, only: fv_nphys
1762 : integer , intent(in) :: nphys,nhalo,nlev
1763 : real(kind=r8), intent(inout):: fld(1-nhalo:nphys+nhalo,1-nhalo:nphys+nhalo,2,nlev)
1764 : type (element_t), intent(in) :: elem
1765 : type(fvm_struct), intent(in), optional :: fvm
1766 : !
1767 : ! local variables
1768 : !
1769 : integer :: i,j,k
1770 : real(r8):: v1,v2
1771 :
1772 0 : if (nhalo==0.and.nphys==np) then
1773 0 : do k=1,nlev
1774 0 : do j=1,nphys
1775 0 : do i=1,nphys
1776 : ! latlon->contra
1777 0 : v1 = fld(i,j,1,k)
1778 0 : v2 = fld(i,j,2,k)
1779 0 : fld(i,j,1,k) = elem%Dinv(i,j,1,1)*v1 + elem%Dinv(i,j,1,2)*v2
1780 0 : fld(i,j,2,k) = elem%Dinv(i,j,2,1)*v1 + elem%Dinv(i,j,2,2)*v2
1781 : enddo
1782 : enddo
1783 : end do
1784 0 : else if (nphys==fv_nphys.and.nhalo.le.fv_nphys) then
1785 0 : do k=1,nlev
1786 0 : do j=1-nhalo,nphys+nhalo
1787 0 : do i=1-nhalo,nphys+nhalo
1788 : ! latlon->contra
1789 0 : v1 = fld(i,j,1,k)
1790 0 : v2 = fld(i,j,2,k)
1791 0 : fld(i,j,1,k) = fvm%Dinv_physgrid(i,j,1,1)*v1 + fvm%Dinv_physgrid(i,j,1,2)*v2
1792 0 : fld(i,j,2,k) = fvm%Dinv_physgrid(i,j,2,1)*v1 + fvm%Dinv_physgrid(i,j,2,2)*v2
1793 : enddo
1794 : enddo
1795 : end do
1796 : else
1797 0 : call endrun('ERROR: vec_latlon_to_contra - grid not supported or halo too large')
1798 : end if
1799 0 : end subroutine vec_latlon_to_contra
1800 0 : end module interpolate_mod
|
