20 subroutine horiz_interp_conserve_new_1dx1d_ ( Interp, lon_in, lat_in, lon_out, lat_out, verbose)
21 type(horiz_interp_type),
intent(inout) :: Interp
22 real(FMS_HI_KIND_),
intent(in),
dimension(:) :: lon_in , lat_in
23 real(FMS_HI_KIND_),
intent(in),
dimension(:) :: lon_out, lat_out
24 integer,
intent(in),
optional :: verbose
27 real(FMS_HI_KIND_),
dimension(size(lat_out(:))-1,2) :: sph
28 real(FMS_HI_KIND_),
dimension(size(lon_out(:))-1,2) :: theta
29 real(FMS_HI_KIND_),
dimension(size(lat_in(:))) :: slat_in
30 real(FMS_HI_KIND_),
dimension(size(lon_in(:))-1) :: dlon_in
31 real(FMS_HI_KIND_),
dimension(size(lat_in(:))-1) :: dsph_in
32 real(FMS_HI_KIND_),
dimension(size(lon_out(:))-1) :: dlon_out
33 real(FMS_HI_KIND_),
dimension(size(lat_out(:))-1) :: dsph_out
34 real(FMS_HI_KIND_) :: blon, fac, hpi, tpi, eps
35 integer,
parameter :: num_iters = 4
36 integer :: i, j, m, n, nlon_in, nlat_in, nlon_out, nlat_out, &
37 iverbose, m2, n2, iter
39 character(len=64) :: mesg
40 integer,
parameter :: kindl = fms_hi_kind_
42 if(.not. module_is_initialized)
call mpp_error(fatal, &
43 'HORIZ_INTERP_CONSERVE_NEW_1DX1D_: horiz_interp_conserve_init is not called')
45 if(great_circle_algorithm)
call mpp_error(fatal, &
46 'HORIZ_INTERP_CONSERVE_NEW_1DX1D_: great_circle_algorithm is not implemented, contact developer')
48 iverbose = 0;
if (
present(verbose)) iverbose = verbose
51 root_pe = mpp_root_pe()
53 hpi = 0.5_kindl * real(pi, fms_hi_kind_)
54 tpi = 4.0_kindl * real(hpi, fms_hi_kind_)
56 nlon_in =
size(lon_in(:))-1; nlat_in =
size(lat_in(:))-1
57 nlon_out =
size(lon_out(:))-1; nlat_out =
size(lat_out(:))-1
59 allocate ( interp % HI_KIND_TYPE_ % facj (nlat_out,2), interp % jlat (nlat_out,2), &
60 interp % HI_KIND_TYPE_ % faci (nlon_out,2), interp % ilon (nlon_out,2), &
61 interp % HI_KIND_TYPE_ % area_src (nlon_in, nlat_in), &
62 interp % HI_KIND_TYPE_ % area_dst (nlon_out, nlat_out) )
68 slat_in(j) = sin(lat_in(j))
72 dsph_in(j) = abs(slat_in(j+1)-slat_in(j))
76 dlon_in(i) = abs(lon_in(i+1)-lon_in(i))
81 if (lat_in(1) > lat_in(nlat_in+1)) s2n = .false.
87 dsph_out(n) = abs(sin(lat_out(n+1))-sin(lat_out(n)))
91 theta(m,1) = lon_out(m)
92 theta(m,2) = lon_out(m+1)
93 dlon_out(m) = abs(lon_out(m+1)-lon_out(m))
96 interp%nlon_src = nlon_in; interp%nlat_src = nlat_in
97 interp%nlon_dst = nlon_out; interp%nlat_dst = nlat_out
104 if (lat_out(n) < lat_out(n+1))
then
105 sph(n,1) = sin(lat_out(n))
106 sph(n,2) = sin(lat_out(n+1))
108 sph(n,1) = sin(lat_out(n+1))
109 sph(n,2) = sin(lat_out(n))
120 if ( (s2n .and. (slat_in(j)-sph(n,n2)) <= eps .and. &
121 (sph(n,n2)-slat_in(j+1)) <= eps) .or. &
122 (.not.s2n .and. (slat_in(j+1)-sph(n,n2)) <= eps .and. &
123 (sph(n,n2)-slat_in(j)) <= eps) )
then
124 interp%jlat(n,n2) = j
126 fac = (sph(n,n2)-slat_in(j))/(slat_in(j+1)-slat_in(j))
128 if (n2 == 1) interp%HI_KIND_TYPE_%facj(n,n2) = 1.0_kindl - fac
129 if (n2 == 2) interp%HI_KIND_TYPE_%facj(n,n2) = fac
131 if (n2 == 1) interp%HI_KIND_TYPE_%facj(n,n2) = fac
132 if (n2 == 2) interp%HI_KIND_TYPE_%facj(n,n2) = 1.0_kindl - fac
137 if ( interp%jlat(n,n2) /= 0 )
exit
140 eps = epsilon(sph)*real(10.0_kindl**iter, kindl)
143 if ( interp%jlat(n,n2) == 0 )
then
144 write (mesg,710) n,sph(n,n2)
145 710
format (
': n,sph=',i3,f14.7,40x)
146 call mpp_error(fatal,
'horiz_interp_conserve_mod:no latitude index found'//trim(mesg))
157 if ( blon < lon_in(1) ) blon = blon + tpi
158 if ( blon > lon_in(nlon_in+1) ) blon = blon - tpi
163 if ( (lon_in(i)-blon) <= eps .and. &
164 (blon-lon_in(i+1)) <= eps )
then
165 interp%ilon(m,m2) = i
166 fac = (blon-lon_in(i))/(lon_in(i+1)-lon_in(i))
167 if (m2 == 1) interp%HI_KIND_TYPE_%faci(m,m2) = 1.0_kindl - fac
168 if (m2 == 2) interp%HI_KIND_TYPE_%faci(m,m2) = fac
172 if ( interp%ilon(m,m2) /= 0 )
exit
175 eps = epsilon(blon)*real(10.0_kindl**iter, kindl)
178 if ( interp%ilon(m,m2) == 0 )
then
179 print *,
'lon_out,blon,blon_in,eps=', &
180 theta(m,m2),blon,lon_in(1),lon_in(nlon_in+1),eps
181 call mpp_error(fatal,
'horiz_interp_conserve_mod: no longitude index found')
190 interp%HI_KIND_TYPE_%area_src(i,j) = dlon_in(i) * dsph_in(j)
198 interp%HI_KIND_TYPE_%area_dst(m,n) = dlon_out(m) * dsph_out(n)
205 if (iverbose > 2)
then
206 write (*,801) (i,interp%ilon(i,1),interp%ilon(i,2), &
207 interp%HI_KIND_TYPE_%faci(i,1),interp%HI_KIND_TYPE_%faci(i,2),i=1,nlon_out)
208 write (*,802) (j,interp%jlat(j,1),interp%jlat(j,2), &
209 interp%HI_KIND_TYPE_%facj(j,1),interp%HI_KIND_TYPE_%facj(j,2),j=1,nlat_out)
210 801
format (/,2x,
'i',4x,
'is',5x,
'ie',4x,
'facis',4x,
'facie', &
212 802
format (/,2x,
'j',4x,
'js',5x,
'je',4x,
'facjs',4x,
'facje', &
217 interp% HI_KIND_TYPE_ % is_allocated = .true.
218 interp% interp_method = conserve
220 end subroutine horiz_interp_conserve_new_1dx1d_
224 subroutine horiz_interp_conserve_new_1dx2d_ ( Interp, lon_in, lat_in, lon_out, lat_out, &
225 mask_in, mask_out, verbose)
226 type(horiz_interp_type),
intent(inout) :: Interp
227 real(FMS_HI_KIND_),
intent(in),
dimension(:) :: lon_in , lat_in
228 real(FMS_HI_KIND_),
intent(in),
dimension(:,:) :: lon_out, lat_out
229 real(FMS_HI_KIND_),
intent(in),
optional,
dimension(:,:) :: mask_in
230 real(FMS_HI_KIND_),
intent(inout),
optional,
dimension(:,:) :: mask_out
231 integer,
intent(in),
optional :: verbose
234 integer :: create_xgrid_1DX2D_order1, get_maxxgrid, maxxgrid
235 integer :: create_xgrid_great_circle
236 integer :: nlon_in, nlat_in, nlon_out, nlat_out, nxgrid, i, j
237 real(r8_kind),
dimension(size(lon_in(:))-1, size(lat_in(:))-1) :: mask_src
238 integer,
allocatable,
dimension(:) :: i_src, j_src, i_dst, j_dst
239 real(r8_kind),
allocatable,
dimension(:) :: xgrid_area, clon, clat
240 real(r8_kind),
allocatable,
dimension(:,:) :: dst_area, lon_src, lat_src
241 real(r8_kind),
allocatable,
dimension(:) :: lat_in_flip
242 real(r8_kind),
allocatable,
dimension(:,:) :: mask_src_flip
243 real(r8_kind),
allocatable,
dimension(:) :: lon_in_r8, lat_in_r8
244 real(r8_kind),
allocatable,
dimension(:,:) :: lon_out_r8, lat_out_r8
246 integer :: nincrease, ndecrease
248 integer(kind=1) :: one_byte(8)
249 integer,
parameter :: kindl = fms_hi_kind_
251 if(.not. module_is_initialized)
call mpp_error(fatal, &
252 'HORIZ_INTERP_CONSERVE_NEW_1DX2D_: horiz_interp_conserve_init is not called')
254 if( (
size(lon_out,1) .NE.
size(lat_out,1)) .OR. (
size(lon_out,2) .NE.
size(lat_out,2)) ) &
255 call mpp_error(fatal,
'horiz_interp_conserve_mod: size mismatch between lon_out and lat_out')
256 nlon_in =
size(lon_in(:)) - 1; nlat_in =
size(lat_in(:)) - 1
257 nlon_out =
size(lon_out,1) - 1; nlat_out =
size(lon_out,2) - 1
259 mask_src = 1.0_r8_kind
260 if(
present(mask_in))
then
261 if( (
size(mask_in,1) .NE. nlon_in) .OR. (
size(mask_in,2) .NE. nlat_in))
call mpp_error(fatal, &
262 'horiz_interp_conserve_mod: size mismatch between mask_in and lon_in/lat_in')
263 mask_src = real(mask_in, r8_kind)
266 maxxgrid = get_maxxgrid()
267 allocate(i_src(maxxgrid), j_src(maxxgrid), i_dst(maxxgrid), j_dst(maxxgrid) )
268 allocate( xgrid_area(maxxgrid), dst_area(nlon_out, nlat_out) )
274 if( lat_in(j+1) > lat_in(j) )
then
275 nincrease = nincrease + 1
276 else if ( lat_in(j+1) < lat_in(j) )
then
277 ndecrease = ndecrease + 1
281 if(nincrease == nlat_in)
then
283 else if(ndecrease == nlat_in)
then
286 call mpp_error(fatal,
'horiz_interp_conserve_mod: nlat_in should be equal to nincreaase or ndecrease')
289 allocate(lon_out_r8(
size(lon_out,1),
size(lon_out,2)))
290 allocate(lat_out_r8(
size(lat_out,1),
size(lat_out,2)))
291 lon_out_r8 = real(lon_out, r8_kind)
292 lat_out_r8 = real(lat_out, r8_kind)
294 if( .not. great_circle_algorithm )
then
296 allocate(lat_in_flip(nlat_in+1), mask_src_flip(nlon_in,nlat_in))
298 lat_in_flip(j) = real(lat_in(nlat_in+2-j), r8_kind)
301 mask_src_flip(:,j) = mask_src(:,nlat_in+1-j)
303 allocate(lon_in_r8(
size(lon_in)))
304 lon_in_r8 = real(lon_in, r8_kind)
305 nxgrid = create_xgrid_1dx2d_order1(nlon_in, nlat_in, nlon_out, nlat_out, lon_in_r8, lat_in_flip, &
306 lon_out_r8, lat_out_r8, mask_src_flip, i_src, j_src, i_dst, j_dst, xgrid_area)
307 deallocate(lon_in_r8, lat_in_flip, mask_src_flip)
309 allocate(lon_in_r8(
size(lon_in)))
310 allocate(lat_in_r8(
size(lat_in)))
311 lon_in_r8 = real(lon_in, r8_kind)
312 lat_in_r8 = real(lat_in, r8_kind)
313 nxgrid = create_xgrid_1dx2d_order1(nlon_in, nlat_in, nlon_out, nlat_out, lon_in_r8, lat_in_r8, lon_out_r8, &
314 & lat_out_r8, mask_src, i_src, j_src, i_dst, j_dst, xgrid_area)
315 deallocate(lon_in_r8,lat_in_r8)
318 allocate(lon_src(nlon_in+1,nlat_in+1), lat_src(nlon_in+1,nlat_in+1))
319 allocate(clon(maxxgrid), clat(maxxgrid))
321 allocate(mask_src_flip(nlon_in,nlat_in))
324 lon_src(i,j) = real(lon_in(i), r8_kind)
325 lat_src(i,j) = real(lat_in(nlat_in+2-j), r8_kind)
329 mask_src_flip(:,j) = mask_src(:,nlat_in+1-j)
331 nxgrid = create_xgrid_great_circle(nlon_in, nlat_in, nlon_out, nlat_out, lon_src, lat_src, lon_out_r8, &
332 & lat_out_r8, mask_src_flip, i_src, j_src, i_dst, j_dst, xgrid_area, clon, clat)
333 deallocate(mask_src_flip)
337 lon_src(i,j) = real(lon_in(i), r8_kind)
338 lat_src(i,j) = real(lat_in(j), r8_kind)
341 nxgrid = create_xgrid_great_circle(nlon_in, nlat_in, nlon_out, nlat_out, lon_src, lat_src, lon_out_r8, &
342 & lat_out_r8, mask_src, i_src, j_src, i_dst, j_dst, xgrid_area, clon, clat)
344 deallocate(lon_src, lat_src, clon, clat)
347 deallocate(lon_out_r8, lat_out_r8)
349 allocate(interp%i_src(nxgrid), interp%j_src(nxgrid) )
350 allocate(interp%i_dst(nxgrid), interp%j_dst(nxgrid) )
351 allocate(interp%HI_KIND_TYPE_%area_frac_dst(nxgrid) )
353 interp%nxgrid = nxgrid
354 interp%i_src = i_src(1:nxgrid)+1
355 interp%j_src = j_src(1:nxgrid)+1
356 if(flip_lat) interp%j_src = nlat_in+1-interp%j_src
357 interp%i_dst = i_dst(1:nxgrid)+1
358 interp%j_dst = j_dst(1:nxgrid)+1
361 dst_area = 0.0_r8_kind
363 dst_area(interp%i_dst(i), interp%j_dst(i)) = dst_area(interp%i_dst(i), interp%j_dst(i)) + xgrid_area(i)
367 interp%HI_KIND_TYPE_%area_frac_dst(i) = real(xgrid_area(i)/dst_area(interp%i_dst(i), interp%j_dst(i) ), &
370 interp%nlon_src = nlon_in; interp%nlat_src = nlat_in
371 interp%nlon_dst = nlon_out; interp%nlat_dst = nlat_out
372 if(
present(mask_out))
then
373 if( (
size(mask_out,1) .NE. nlon_out) .OR. (
size(mask_out,2) .NE. nlat_out) )
call mpp_error(fatal, &
374 'horiz_interp_conserve_mod: size mismatch between mask_out and lon_out/lat_out')
377 mask_out(interp%i_dst(i),interp%j_dst(i)) = mask_out(interp%i_dst(i), &
378 & interp%j_dst(i)) + interp%HI_KIND_TYPE_%area_frac_dst(i)
382 deallocate(i_src, j_src, i_dst, j_dst, xgrid_area, dst_area )
384 interp% HI_KIND_TYPE_ % is_allocated = .true.
385 interp% interp_method = conserve
387 end subroutine horiz_interp_conserve_new_1dx2d_
391 subroutine horiz_interp_conserve_new_2dx1d_ ( Interp, lon_in, lat_in, lon_out, lat_out, &
392 mask_in, mask_out, verbose)
393 type(horiz_interp_type),
intent(inout) :: Interp
394 real(FMS_HI_KIND_),
intent(in),
dimension(:,:) :: lon_in , lat_in
395 real(FMS_HI_KIND_),
intent(in),
dimension(:) :: lon_out, lat_out
396 real(FMS_HI_KIND_),
intent(in),
optional,
dimension(:,:) :: mask_in
397 real(FMS_HI_KIND_),
intent(inout),
optional,
dimension(:,:) :: mask_out
398 integer,
intent(in),
optional :: verbose
400 integer :: create_xgrid_2DX1D_order1, get_maxxgrid, maxxgrid
401 integer :: create_xgrid_great_circle
402 integer :: nlon_in, nlat_in, nlon_out, nlat_out, nxgrid, i, j
403 integer,
allocatable,
dimension(:) :: i_src, j_src, i_dst, j_dst
404 real(r8_kind),
allocatable,
dimension(:,:) :: dst_area
405 real(r8_kind),
dimension(size(lon_in,1)-1, size(lon_in,2)-1) :: mask_src
406 real(r8_kind),
allocatable,
dimension(:) :: xgrid_area, clon, clat
407 real(r8_kind),
allocatable,
dimension(:) :: lon_out_r8, lat_out_r8
408 real(r8_kind),
allocatable,
dimension(:,:) :: lon_in_r8, lat_in_r8
409 real(r8_kind),
allocatable,
dimension(:,:) :: lon_dst, lat_dst
410 integer(kind=1) :: one_byte(8)
411 integer,
parameter :: kindl = fms_hi_kind_
413 if(.not. module_is_initialized)
call mpp_error(fatal, &
414 'HORIZ_INTERP_CONSERVE_NEW_2DX1D_: horiz_interp_conserve_init is not called')
416 if( (
size(lon_in,1) .NE.
size(lat_in,1)) .OR. (
size(lon_in,2) .NE.
size(lat_in,2)) ) &
417 call mpp_error(fatal,
'horiz_interp_conserve_mod: size mismatch between lon_in and lat_in')
418 nlon_in =
size(lon_in,1) - 1; nlat_in =
size(lon_in,2) - 1
419 nlon_out =
size(lon_out(:)) - 1; nlat_out =
size(lat_out(:)) - 1
421 mask_src = 1.0_r8_kind
422 if(
present(mask_in))
then
423 if( (
size(mask_in,1) .NE. nlon_in) .OR. (
size(mask_in,2) .NE. nlat_in))
call mpp_error(fatal, &
424 'horiz_interp_conserve_mod: size mismatch between mask_in and lon_in/lat_in')
425 mask_src = real(mask_in, r8_kind)
428 maxxgrid = get_maxxgrid()
429 allocate(i_src(maxxgrid), j_src(maxxgrid), i_dst(maxxgrid), j_dst(maxxgrid) )
430 allocate( xgrid_area(maxxgrid), dst_area(nlon_out, nlat_out) )
432 allocate(lon_in_r8(
size(lon_in,1),
size(lon_in, 2)))
433 allocate(lat_in_r8(
size(lat_in,1),
size(lat_in, 2)))
434 allocate(lon_out_r8(
size(lon_out)))
435 allocate(lat_out_r8(
size(lat_out)))
436 lon_out_r8 = real(lon_out, r8_kind)
437 lat_out_r8 = real(lat_out, r8_kind)
438 lon_in_r8 = real(lon_in, r8_kind)
439 lat_in_r8 = real(lat_in, r8_kind)
441 if( .not. great_circle_algorithm )
then
442 nxgrid = create_xgrid_2dx1d_order1(nlon_in, nlat_in, nlon_out, nlat_out, lon_in_r8, lat_in_r8, &
443 lon_out_r8, lat_out_r8, mask_src, i_src, j_src, i_dst, j_dst, xgrid_area)
445 allocate(lon_dst(nlon_out+1, nlat_out+1), lat_dst(nlon_out+1, nlat_out+1) )
446 allocate(clon(maxxgrid), clat(maxxgrid))
449 lon_dst(i,j) = real(lon_out(i), r8_kind)
450 lat_dst(i,j) = real(lat_out(j), r8_kind)
453 nxgrid = create_xgrid_great_circle(nlon_in, nlat_in, nlon_out, nlat_out, lon_in_r8, lat_in_r8, lon_dst, &
454 & lat_dst, mask_src, i_src, j_src, i_dst, j_dst, xgrid_area, clon, clat)
456 deallocate(lon_out_r8,lat_out_r8, lon_in_r8, lat_in_r8)
457 allocate(interp%i_src(nxgrid), interp%j_src(nxgrid) )
458 allocate(interp%i_dst(nxgrid), interp%j_dst(nxgrid) )
459 allocate(interp%HI_KIND_TYPE_%area_frac_dst(nxgrid) )
461 interp%nxgrid = nxgrid
462 interp%i_src = i_src(1:nxgrid)+1
463 interp%j_src = j_src(1:nxgrid)+1
464 interp%i_dst = i_dst(1:nxgrid)+1
465 interp%j_dst = j_dst(1:nxgrid)+1
468 dst_area = 0.0_r8_kind
470 dst_area(interp%i_dst(i), interp%j_dst(i)) = dst_area(interp%i_dst(i), interp%j_dst(i)) + xgrid_area(i)
474 interp%HI_KIND_TYPE_%area_frac_dst(i) = real(xgrid_area(i)/dst_area(interp%i_dst(i), interp%j_dst(i) ), &
477 interp%nlon_src = nlon_in; interp%nlat_src = nlat_in
478 interp%nlon_dst = nlon_out; interp%nlat_dst = nlat_out
479 if(
present(mask_out))
then
480 if( (
size(mask_out,1) .NE. nlon_out) .OR. (
size(mask_out,2) .NE. nlat_out) )
call mpp_error(fatal, &
481 'horiz_interp_conserve_mod: size mismatch between mask_out and lon_out/lat_out')
484 mask_out(interp%i_dst(i),interp%j_dst(i)) = mask_out(interp%i_dst(i), &
485 & interp%j_dst(i)) + interp%HI_KIND_TYPE_%area_frac_dst(i)
489 deallocate(i_src, j_src, i_dst, j_dst, xgrid_area, dst_area)
491 interp% HI_KIND_TYPE_ % is_allocated = .true.
492 interp% interp_method = conserve
494 end subroutine horiz_interp_conserve_new_2dx1d_
498 subroutine horiz_interp_conserve_new_2dx2d_ ( Interp, lon_in, lat_in, lon_out, lat_out, &
499 mask_in, mask_out, verbose)
500 type(horiz_interp_type),
intent(inout) :: Interp
501 real(FMS_HI_KIND_),
intent(in),
dimension(:,:) :: lon_in , lat_in
502 real(FMS_HI_KIND_),
intent(in),
dimension(:,:) :: lon_out, lat_out
503 real(FMS_HI_KIND_),
intent(in),
optional,
dimension(:,:) :: mask_in
504 real(FMS_HI_KIND_),
intent(inout),
optional,
dimension(:,:) :: mask_out
505 integer,
intent(in),
optional :: verbose
507 integer :: create_xgrid_2DX2D_order1, get_maxxgrid, maxxgrid
508 integer :: create_xgrid_great_circle
509 integer :: nlon_in, nlat_in, nlon_out, nlat_out, nxgrid, i
510 integer,
allocatable,
dimension(:) :: i_src, j_src, i_dst, j_dst
511 real(r8_kind),
dimension(size(lon_in,1)-1, size(lon_in,2)-1) :: mask_src
512 real(r8_kind),
allocatable,
dimension(:) :: xgrid_area, clon, clat
513 real(r8_kind),
allocatable,
dimension(:,:) :: dst_area
514 real(r8_kind),
allocatable,
dimension(:,:) :: lon_in_r8, lat_in_r8
515 real(r8_kind),
allocatable,
dimension(:,:) :: lon_out_r8, lat_out_r8
516 integer(kind=1) :: one_byte(8)
517 integer,
parameter :: kindl = fms_hi_kind_
519 if(.not. module_is_initialized)
call mpp_error(fatal, &
520 'HORIZ_INTERP_CONSERVE_NEW_2DX2D_: horiz_interp_conserve_init is not called')
522 if( (
size(lon_in,1) .NE.
size(lat_in,1)) .OR. (
size(lon_in,2) .NE.
size(lat_in,2)) ) &
523 call mpp_error(fatal,
'horiz_interp_conserve_mod: size mismatch between lon_in and lat_in')
524 if( (
size(lon_out,1) .NE.
size(lat_out,1)) .OR. (
size(lon_out,2) .NE.
size(lat_out,2)) ) &
525 call mpp_error(fatal,
'horiz_interp_conserve_mod: size mismatch between lon_out and lat_out')
526 nlon_in =
size(lon_in,1) - 1; nlat_in =
size(lon_in,2) - 1
527 nlon_out =
size(lon_out,1) - 1; nlat_out =
size(lon_out,2) - 1
529 mask_src = 1.0_r8_kind
530 if(
present(mask_in))
then
531 if( (
size(mask_in,1) .NE. nlon_in) .OR. (
size(mask_in,2) .NE. nlat_in))
call mpp_error(fatal, &
532 'horiz_interp_conserve_mod: size mismatch between mask_in and lon_in/lat_in')
533 mask_src = real(mask_in, r8_kind)
536 maxxgrid = get_maxxgrid()
537 allocate(i_src(maxxgrid), j_src(maxxgrid), i_dst(maxxgrid), j_dst(maxxgrid) )
538 allocate( xgrid_area(maxxgrid), dst_area(nlon_out, nlat_out) )
540 allocate(lon_in_r8(
size(lon_in,1),
size(lon_in,2)))
541 allocate(lat_in_r8(
size(lat_in,1),
size(lat_in,2)))
542 allocate(lon_out_r8(
size(lon_out,1),
size(lon_out,2)))
543 allocate(lat_out_r8(
size(lat_out,1),
size(lat_out,2)))
544 lon_in_r8 = real(lon_in,r8_kind)
545 lat_in_r8 = real(lat_in, r8_kind)
546 lon_out_r8 = real(lon_out, r8_kind)
547 lat_out_r8 = real(lat_out, r8_kind)
549 if( .not. great_circle_algorithm )
then
550 nxgrid = create_xgrid_2dx2d_order1(nlon_in, nlat_in, nlon_out, nlat_out, lon_in_r8, lat_in_r8, lon_out_r8, &
551 & lat_out_r8, mask_src, i_src, j_src, i_dst, j_dst, xgrid_area)
553 allocate(clon(maxxgrid), clat(maxxgrid))
554 nxgrid = create_xgrid_great_circle(nlon_in, nlat_in, nlon_out, nlat_out, lon_in_r8, lat_in_r8, lon_out_r8, &
555 & lat_out_r8, mask_src, i_src, j_src, i_dst, j_dst, xgrid_area, clon, clat)
556 deallocate(clon, clat)
559 deallocate(lon_in_r8, lat_in_r8, lon_out_r8, lat_out_r8)
561 allocate(interp%i_src(nxgrid), interp%j_src(nxgrid) )
562 allocate(interp%i_dst(nxgrid), interp%j_dst(nxgrid) )
563 allocate(interp%HI_KIND_TYPE_%area_frac_dst(nxgrid) )
565 interp%nxgrid = nxgrid
566 interp%i_src = i_src(1:nxgrid)+1
567 interp%j_src = j_src(1:nxgrid)+1
568 interp%i_dst = i_dst(1:nxgrid)+1
569 interp%j_dst = j_dst(1:nxgrid)+1
572 dst_area = 0.0_r8_kind
574 dst_area(interp%i_dst(i), interp%j_dst(i)) = dst_area(interp%i_dst(i), interp%j_dst(i)) + xgrid_area(i)
578 interp%HI_KIND_TYPE_%area_frac_dst(i) = real(xgrid_area(i)/dst_area(interp%i_dst(i), interp%j_dst(i)), &
582 interp%nlon_src = nlon_in; interp%nlat_src = nlat_in
583 interp%nlon_dst = nlon_out; interp%nlat_dst = nlat_out
584 if(
present(mask_out))
then
585 if( (
size(mask_out,1) .NE. nlon_out) .OR. (
size(mask_out,2) .NE. nlat_out) )
call mpp_error(fatal, &
586 'horiz_interp_conserve_mod: size mismatch between mask_out and lon_out/lat_out')
589 mask_out(interp%i_dst(i),interp%j_dst(i)) = mask_out(interp%i_dst(i), &
590 & interp%j_dst(i)) + interp%HI_KIND_TYPE_%area_frac_dst(i)
594 deallocate(i_src, j_src, i_dst, j_dst, xgrid_area, dst_area )
596 interp% HI_KIND_TYPE_ % is_allocated = .true.
597 interp% interp_method = conserve
599 end subroutine horiz_interp_conserve_new_2dx2d_
610 type (horiz_interp_type),
intent(in) :: Interp
611 real(FMS_HI_KIND_),
intent(in),
dimension(:,:) :: data_in
612 real(FMS_HI_KIND_),
intent(out),
dimension(:,:) :: data_out
613 integer,
intent(in),
optional :: verbose
615 real(FMS_HI_KIND_),
intent(in),
dimension(:,:),
optional :: mask_in
620 real(FMS_HI_KIND_),
intent(out),
dimension(:,:),
optional :: mask_out
626 if (
size(data_in,1) /= interp%nlon_src .or.
size(data_in,2) /= interp%nlat_src) &
627 call mpp_error(fatal,
'horiz_interp_conserve_mod: size of input array incorrect')
629 if (
size(data_out,1) /= interp%nlon_dst .or.
size(data_out,2) /= interp%nlat_dst) &
630 call mpp_error(fatal,
'horiz_interp_conserve_mod: size of output array incorrect')
632 select case ( interp%version)
634 call horiz_interp_conserve_version1(interp, data_in, data_out, verbose, mask_in, mask_out)
636 if(
present(mask_in) .OR.
present(mask_out) )
call mpp_error(fatal,
'HORIZ_INTERP_CONSERVE_:'// &
637 &
' for version 2, mask_in and mask_out must be passed in horiz_interp_new, not in horiz_interp')
638 call horiz_interp_conserve_version2(interp, data_in, data_out, verbose)
644 subroutine horiz_interp_conserve_version1_ ( Interp, data_in, data_out, verbose, &
647 type (horiz_interp_type),
intent(in) :: Interp
648 real(FMS_HI_KIND_),
intent(in),
dimension(:,:) :: data_in
649 real(FMS_HI_KIND_),
intent(out),
dimension(:,:) :: data_out
650 integer,
intent(in),
optional :: verbose
651 real(FMS_HI_KIND_),
intent(in),
dimension(:,:),
optional :: mask_in
652 real(FMS_HI_KIND_),
intent(out),
dimension(:,:),
optional :: mask_out
654 integer :: m, n, nlon_in, nlat_in, nlon_out, nlat_out, &
655 miss_in, miss_out, is, ie, js, je, &
657 real(FMS_HI_KIND_) :: dsum, wsum, avg_in, min_in, max_in, &
658 avg_out, min_out, max_out, eps, asum, &
659 dwtsum, wtsum, arsum, fis, fie, fjs, fje
660 integer,
parameter :: kindl = fms_hi_kind_
662 iverbose = 0;
if (
present(verbose)) iverbose = verbose
666 nlon_in = interp%nlon_src; nlat_in = interp%nlat_src
667 nlon_out = interp%nlon_dst; nlat_out = interp%nlat_dst
669 if (
present(mask_in))
then
670 if ( count(mask_in < -.0001_kindl .or. mask_in > 1.0001_kindl) > 0 ) &
671 call mpp_error(fatal,
'horiz_interp_conserve_mod: input mask not between 0,1')
681 if (interp%jlat(n,1) <= interp%jlat(n,2))
then
682 js = interp%jlat(n,1); je = interp%jlat(n,2)
683 fjs = interp%HI_KIND_TYPE_%facj(n,1); fje = interp%HI_KIND_TYPE_%facj(n,2)
685 js = interp%jlat(n,2); je = interp%jlat(n,1)
686 fjs = interp%HI_KIND_TYPE_%facj(n,2); fje = interp%HI_KIND_TYPE_%facj(n,1)
691 is = interp%ilon(m,1); ie = interp%ilon(m,2)
692 fis = interp%HI_KIND_TYPE_%faci(m,1); fie = interp%HI_KIND_TYPE_%faci(m,2)
711 ie = interp%ilon(m,2)
712 fie = interp%HI_KIND_TYPE_%faci(m,2)
716 if (
present(mask_in))
then
717 call data_sum( data_in(is:ie,js:je), interp%HI_KIND_TYPE_%area_src(is:ie,js:je), &
718 fis, fie, fjs,fje, dwtsum, wtsum, arsum, mask_in(is:ie,js:je) )
719 else if(
allocated(interp%HI_KIND_TYPE_%mask_in) )
then
720 call data_sum( data_in(is:ie,js:je), interp%HI_KIND_TYPE_%area_src(is:ie,js:je), &
721 fis, fie, fjs,fje, dwtsum, wtsum, arsum, interp%HI_KIND_TYPE_%mask_in(is:ie,js:je) )
723 call data_sum( data_in(is:ie,js:je), interp%HI_KIND_TYPE_%area_src(is:ie,js:je), &
724 fis, fie, fjs,fje, dwtsum, wtsum, arsum )
728 if (wtsum > eps)
then
729 data_out(m,n) = dwtsum/wtsum
730 if (
present(mask_out)) mask_out(m,n) = wtsum/arsum
732 data_out(m,n) = 0.0_kindl
733 if (
present(mask_out)) mask_out(m,n) = 0.0_kindl
743 if (iverbose > 0)
then
747 call stats(data_in, interp%HI_KIND_TYPE_%area_src, asum, dsum, wsum, min_in, max_in, miss_in, mask_in)
750 if(pe == root_pe)
then
751 if (wsum > 0.0_kindl)
then
754 print *,
'horiz_interp stats: input area equals zero '
757 if (iverbose > 1) print
'(2f16.11)',
'global sum area_in = ', asum, wsum
761 call stats(data_out, interp%HI_KIND_TYPE_%area_dst, asum, dsum, wsum, min_out, max_out, miss_out, mask_out)
763 if(pe == root_pe)
then
764 if (wsum > 0.0_kindl )
then
767 print *,
'horiz_interp stats: output area equals zero '
770 if (iverbose > 1) print
'(2f16.11)',
'global sum area_out = ', asum, wsum
774 if(pe == root_pe)
then
776 write (*,901) min_in ,max_in ,avg_in
777 if (
present(mask_in))
write (*,903) miss_in
778 write (*,902) min_out,max_out,avg_out
779 if (
present(mask_out))
write (*,903) miss_out
782 900
format (/,1x,10(
'-'),
' output from horiz_interp ',10(
'-'))
783 901
format (
' input: min=',f16.9,
' max=',f16.9,
' avg=',f22.15)
784 902
format (
' output: min=',f16.9,
' max=',f16.9,
' avg=',f22.15)
785 903
format (
' number of missing points = ',i6)
790 end subroutine horiz_interp_conserve_version1_
793 subroutine horiz_interp_conserve_version2_ ( Interp, data_in, data_out, verbose )
795 type (horiz_interp_type),
intent(in) :: Interp
796 real(FMS_HI_KIND_),
intent(in),
dimension(:,:) :: data_in
797 real(FMS_HI_KIND_),
intent(out),
dimension(:,:) :: data_out
798 integer,
intent(in),
optional :: verbose
799 integer :: i, i_src, j_src, i_dst, j_dst
800 integer,
parameter :: kindl = fms_hi_kind_
803 do i = 1, interp%nxgrid
804 i_src = interp%i_src(i); j_src = interp%j_src(i)
805 i_dst = interp%i_dst(i); j_dst = interp%j_dst(i)
806 data_out(i_dst, j_dst) = data_out(i_dst, j_dst) + data_in(i_src,j_src)*interp%HI_KIND_TYPE_%area_frac_dst(i)
809 end subroutine horiz_interp_conserve_version2_
814 subroutine stats_ ( dat, area, asum, dsum, wsum, low, high, miss, mask )
815 real(FMS_HI_KIND_),
intent(in) :: dat(:,:), area(:,:)
816 real(FMS_HI_KIND_),
intent(out) :: asum, dsum, wsum, low, high
817 integer,
intent(out) :: miss
818 real(FMS_HI_KIND_),
intent(in),
optional :: mask(:,:)
819 integer,
parameter :: kindl = fms_hi_kind_
821 integer :: pe, root_pe, npes, p, buffer_int(1)
822 real(FMS_HI_KIND_) :: buffer_real(5)
825 root_pe = mpp_root_pe()
830 if (
present(mask))
then
831 asum = sum(area(:,:))
832 dsum = sum(area(:,:)*dat(:,:)*mask(:,:))
833 wsum = sum(area(:,:)*mask(:,:))
834 miss = count(mask(:,:) <= 0.5_kindl)
835 low = minval(dat(:,:),mask=mask(:,:) > 0.5_kindl )
836 high = maxval(dat(:,:),mask=mask(:,:) > 0.5_kindl )
838 asum = sum(area(:,:))
839 dsum = sum(area(:,:)*dat(:,:))
840 wsum = sum(area(:,:))
842 low = minval(dat(:,:))
843 high = maxval(dat(:,:))
849 if(pe == root_pe)
then
852 call mpp_recv(buffer_real(1),glen=5,from_pe=root_pe+p, tag=comm_tag_1)
853 asum = asum + buffer_real(1)
854 dsum = dsum + buffer_real(2)
855 wsum = wsum + buffer_real(3)
856 low = min(low, buffer_real(4))
857 high = max(high, buffer_real(5))
858 call mpp_recv(buffer_int(1),glen=1,from_pe=root_pe+p, tag=comm_tag_2)
859 miss = miss + buffer_int(1)
862 buffer_real(1) = asum
863 buffer_real(2) = dsum
864 buffer_real(3) = wsum
866 buffer_real(5) = high
868 call mpp_send(buffer_real(1),plen=5,to_pe=root_pe, tag=comm_tag_1)
870 call mpp_send(buffer_int(1),plen=1,to_pe=root_pe, tag=comm_tag_2)
880 subroutine data_sum_( grid_data, area, facis, facie, facjs, facje, &
881 dwtsum, wtsum, arsum, mask )
884 real(FMS_HI_KIND_),
intent(in),
dimension(:,:) :: grid_data, area
885 real(FMS_HI_KIND_),
intent(in) :: facis, facie, facjs, facje
886 real(FMS_HI_KIND_),
intent(inout) :: dwtsum, wtsum, arsum
887 real(FMS_HI_KIND_),
intent(in),
optional :: mask(:,:)
895 real(FMS_HI_KIND_),
dimension(size(area,1),size(area,2)) :: wt
896 real(FMS_HI_KIND_) :: asum
900 id=
size(area,1); jd=
size(area,2)
903 wt( 1,:)=wt( 1,:)*facis
904 wt(id,:)=wt(id,:)*facie
905 wt(:, 1)=wt(:, 1)*facjs
906 wt(:,jd)=wt(:,jd)*facje
911 if (
present(mask))
then
913 dwtsum = dwtsum + sum(wt*grid_data)
914 wtsum = wtsum + sum(wt)
916 dwtsum = dwtsum + sum(wt*grid_data)
subroutine horiz_interp_conserve_(Interp, data_in, data_out, verbose, mask_in, mask_out)
Subroutine for performing the horizontal interpolation between two grids.
subroutine stats_(dat, area, asum, dsum, wsum, low, high, miss, mask)
This statistics is for conservative scheme.
subroutine data_sum_(grid_data, area, facis, facie, facjs, facje, dwtsum, wtsum, arsum, mask)
sums up the data and weights for a single output grid box
subroutine mpp_sync_self(pelist, check, request, msg_size, msg_type)
This is to check if current PE's outstanding puts are complete but we can't use shmem_fence because w...
integer function mpp_npes()
Returns processor count for current pelist.
integer function mpp_pe()
Returns processor ID.