libsim Versione 7.2.1

◆ grid_transform_init()

subroutine grid_transform_init ( type(grid_transform), intent(out)  this,
type(transform_def), intent(in)  trans,
type(griddim_def), intent(inout)  in,
type(griddim_def), intent(inout)  out,
real, dimension(:,:), intent(in), optional  maskgrid,
real, dimension(:), intent(in), optional  maskbounds,
character(len=*), intent(in), optional  categoryappend 
)
private

Constructor for a grid_transform object, defining a particular grid-to-grid transformation.

It defines an object describing a transformation from one rectangular grid to another; the abstract type of transformation is described in the transformation object trans (type transform_def) which must have been properly initialised. The additional information required here is the description of the input grid in (type griddim_def), the description of the output grid out (type griddim_def as well). The description of the output grid must be initialized for interpolating type transformations ('inter' and 'boxinter'), while it is generated by this constructor and returned in output for 'zoom', 'boxregrid', 'maskgen' and 'polyinter' transformations.

The generated grid_transform object is specific to the input and output grids involved. The function c_e can be used in order to check whether the object has been successfully initialised, if the result is .FALSE., it should not be used further on.

Parametri
[out]thisgrid_transformation object
[in]transtransformation object
[in,out]ingriddim object to transform
[in,out]outgriddim object defining target grid (input or output depending on type of transformation)
[in]maskgrid2D field to be used for defining valid points, it must have the same shape as the field to be interpolated (for transformation type 'metamorphosis:maskvalid')
[in]maskboundsarray of boundary values for defining a subset of valid points where the values of maskgrid are within the first and last value of maskbounds (for transformation type 'metamorphosis:maskvalid/settoinvalid' and others)
[in]categoryappendappend this suffix to log4fortran namespace category

Definizione alla linea 1409 del file grid_transform_class.F90.

1411 lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1412 this%inter_index_x, this%inter_index_y)
1413
1414 ALLOCATE(this%inter_x(this%innx,this%inny), &
1415 this%inter_y(this%innx,this%inny))
1416 ALLOCATE(this%inter_xp(this%outnx,this%outny), &
1417 this%inter_yp(this%outnx,this%outny))
1418
1419! compute coordinates of input grid
1420 CALL griddim_gen_coord(in, this%inter_x, this%inter_y)
1421! compute coordinates of output grid in input system
1422 CALL proj(in, lout%dim%lon, lout%dim%lat, this%inter_xp, this%inter_yp)
1423
1424 ELSE ! near, shapiro_near
1425 CALL this%find_index(in, .true., &
1426 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1427 lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1428 this%inter_index_x, this%inter_index_y)
1429
1430 IF (this%trans%trans_type == 'intersearch') THEN ! replicate code above
1431 ALLOCATE(this%inter_x(this%innx,this%inny), &
1432 this%inter_y(this%innx,this%inny))
1433 ALLOCATE(this%inter_xp(this%outnx,this%outny), &
1434 this%inter_yp(this%outnx,this%outny))
1435
1436! compute coordinates of input grid
1437 CALL griddim_gen_coord(in, this%inter_x, this%inter_y)
1438! compute coordinates of output grid in input system
1439 CALL proj(in, lout%dim%lon, lout%dim%lat, this%inter_xp, this%inter_yp)
1440 ENDIF
1441 ENDIF
1442
1443 CALL delete(lout)
1444 this%valid = .true.
1445 ENDIF
1446
1447ELSE IF (this%trans%trans_type == 'boxinter') THEN
1448
1449 CALL outgrid_setup() ! common setup for grid-generating methods
1450 CALL get_val(in, nx=this%innx, ny=this%inny)
1451 CALL get_val(out, nx=this%outnx, ny=this%outny, &
1452 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1453! TODO now box size is ignored
1454! if box size not provided, use the actual grid step
1455 IF (.NOT.c_e(this%trans%area_info%boxdx)) &
1456 CALL get_val(out, dx=this%trans%area_info%boxdx)
1457 IF (.NOT.c_e(this%trans%area_info%boxdy)) &
1458 CALL get_val(out, dx=this%trans%area_info%boxdy)
1459! half size is actually needed
1460 this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
1461 this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
1462! unlike before, here index arrays must have the shape of input grid
1463 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1464 this%inter_index_y(this%innx,this%inny))
1465
1466! compute coordinates of input grid in geo system
1467 CALL copy(in, lin)
1468 CALL unproj(lin)
1469! use find_index in the opposite way, here extrap does not make sense
1470 CALL this%find_index(out, .true., &
1471 this%outnx, this%outny, xmin, xmax, ymin, ymax, &
1472 lin%dim%lon, lin%dim%lat, .false., &
1473 this%inter_index_x, this%inter_index_y)
1474
1475 CALL delete(lin)
1476 this%valid = .true. ! warning, no check of subtype
1477
1478ELSE IF (this%trans%trans_type == 'stencilinter') THEN
1479
1480 CALL outgrid_setup() ! common setup for grid-generating methods
1481! from inter:near
1482 CALL get_val(in, nx=this%innx, ny=this%inny, &
1483 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1484 CALL get_val(out, nx=this%outnx, ny=this%outny)
1485
1486 ALLOCATE (this%inter_index_x(this%outnx,this%outny), &
1487 this%inter_index_y(this%outnx,this%outny))
1488 CALL copy(out, lout)
1489 CALL unproj(lout)
1490 CALL this%find_index(in, .true., &
1491 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1492 lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1493 this%inter_index_x, this%inter_index_y)
1494
1495! define the stencil mask
1496 nr = int(this%trans%area_info%radius) ! integer radius
1497 n = nr*2+1 ! n. of points
1498 nm = nr + 1 ! becomes index of center
1499 r2 = this%trans%area_info%radius**2
1500 ALLOCATE(this%stencil(n,n))
1501 this%stencil(:,:) = .true.
1502 DO iy = 1, n
1503 DO ix = 1, n
1504 IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
1505 ENDDO
1506 ENDDO
1507
1508! set to missing the elements of inter_index too close to the domain
1509! borders
1510 xnmin = nr + 1
1511 xnmax = this%innx - nr
1512 ynmin = nr + 1
1513 ynmax = this%inny - nr
1514 DO iy = 1, this%outny
1515 DO ix = 1, this%outnx
1516 IF (this%inter_index_x(ix,iy) < xnmin .OR. &
1517 this%inter_index_x(ix,iy) > xnmax .OR. &
1518 this%inter_index_y(ix,iy) < ynmin .OR. &
1519 this%inter_index_y(ix,iy) > ynmax) THEN
1520 this%inter_index_x(ix,iy) = imiss
1521 this%inter_index_y(ix,iy) = imiss
1522 ENDIF
1523 ENDDO
1524 ENDDO
1525
1526#ifdef DEBUG
1527 CALL l4f_category_log(this%category, l4f_debug, &
1528 'stencilinter: stencil size '//t2c(n*n))
1529 CALL l4f_category_log(this%category, l4f_debug, &
1530 'stencilinter: stencil points '//t2c(count(this%stencil)))
1531#endif
1532
1533 CALL delete(lout)
1534 this%valid = .true. ! warning, no check of subtype
1535
1536ELSE IF (this%trans%trans_type == 'maskgen') THEN
1537
1538 IF (this%trans%sub_type == 'poly') THEN
1539
1540 CALL copy(in, out)
1541 CALL get_val(in, nx=this%innx, ny=this%inny)
1542 this%outnx = this%innx
1543 this%outny = this%inny
1544
1545! unlike before, here index arrays must have the shape of input grid
1546 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1547 this%inter_index_y(this%innx,this%inny))
1548 this%inter_index_x(:,:) = imiss
1549 this%inter_index_y(:,:) = 1
1550
1551! compute coordinates of input grid in geo system
1552 CALL unproj(out) ! should be unproj(lin)
1553
1554 nprev = 1
1555!$OMP PARALLEL DEFAULT(SHARED)
1556!$OMP DO PRIVATE(j, i, n, point) FIRSTPRIVATE(nprev)
1557 DO j = 1, this%inny
1558 inside_x: DO i = 1, this%innx
1559 point = georef_coord_new(x=out%dim%lon(i,j), y=out%dim%lat(i,j))
1560
1561 DO n = nprev, this%trans%poly%arraysize ! optimize starting from last matched polygon
1562 IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
1563 this%inter_index_x(i,j) = n
1564 nprev = n
1565 cycle inside_x
1566 ENDIF
1567 ENDDO
1568 DO n = nprev-1, 1, -1 ! test the other polygons
1569 IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
1570 this%inter_index_x(i,j) = n
1571 nprev = n
1572 cycle inside_x
1573 ENDIF
1574 ENDDO
1575
1576! CALL delete(point) ! speedup
1577 ENDDO inside_x
1578 ENDDO
1579!$OMP END PARALLEL
1580
1581 ELSE IF (this%trans%sub_type == 'grid') THEN
1582! here out(put grid) is abused for indicating the box-generating grid
1583! but the real output grid is the input grid
1584 CALL copy(out, lout) ! save out for local use
1585 CALL delete(out) ! needed before copy
1586 CALL copy(in, out)
1587 CALL get_val(in, nx=this%innx, ny=this%inny)
1588 this%outnx = this%innx
1589 this%outny = this%inny
1590 CALL get_val(lout, nx=nx, ny=ny, &
1591 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1592
1593! unlike before, here index arrays must have the shape of input grid
1594 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1595 this%inter_index_y(this%innx,this%inny))
1596
1597! compute coordinates of input/output grid in geo system
1598 CALL unproj(out)
1599
1600! use find_index in the opposite way, here extrap does not make sense
1601 CALL this%find_index(lout, .true., &
1602 nx, ny, xmin, xmax, ymin, ymax, &
1603 out%dim%lon, out%dim%lat, .false., &
1604 this%inter_index_x, this%inter_index_y)
1605! transform indices to 1-d for mask generation
1606 WHERE(c_e(this%inter_index_x(:,:)))
1607 this%inter_index_x(:,:) = this%inter_index_x(:,:) + &
1608 (this%inter_index_y(:,:)-1)*nx
1609 END WHERE
1610
1611 CALL delete(lout)
1612 ENDIF
1613
1614 this%valid = .true.
1615
1616ELSE IF (this%trans%trans_type == 'polyinter') THEN
1617
1618! this is the only difference wrt maskgen:poly
1619 this%recur = .true. ! grid-to-grid polyinter is done in two steps!
1620
1621 CALL copy(in, out)
1622 CALL get_val(in, nx=this%innx, ny=this%inny)
1623 this%outnx = this%innx
1624 this%outny = this%inny
1625
1626! unlike before, here index arrays must have the shape of input grid
1627 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1628 this%inter_index_y(this%innx,this%inny))
1629 this%inter_index_x(:,:) = imiss
1630 this%inter_index_y(:,:) = 1
1631
1632! compute coordinates of input grid in geo system
1633 CALL unproj(out) ! should be unproj(lin)
1634
1635 nprev = 1
1636!$OMP PARALLEL DEFAULT(SHARED)
1637!$OMP DO PRIVATE(j, i, n, point) FIRSTPRIVATE(nprev)
1638 DO j = 1, this%inny
1639 inside_x_2: DO i = 1, this%innx
1640 point = georef_coord_new(x=out%dim%lon(i,j), y=out%dim%lat(i,j))
1641
1642 DO n = nprev, this%trans%poly%arraysize ! optimize starting from last matched polygon
1643 IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
1644 this%inter_index_x(i,j) = n
1645 nprev = n
1646 cycle inside_x_2
1647 ENDIF
1648 ENDDO
1649 DO n = nprev-1, 1, -1 ! test the other polygons
1650 IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
1651 this%inter_index_x(i,j) = n
1652 nprev = n
1653 cycle inside_x_2
1654 ENDIF
1655 ENDDO
1656
1657! CALL delete(point) ! speedup
1658 ENDDO inside_x_2
1659 ENDDO
1660!$OMP END PARALLEL
1661
1662 this%valid = .true. ! warning, no check of subtype
1663
1664ELSE IF (this%trans%trans_type == 'metamorphosis') THEN
1665
1666 CALL copy(in, out)
1667 CALL get_val(in, nx=this%innx, ny=this%inny)
1668 this%outnx = this%innx
1669 this%outny = this%inny
1670
1671 IF (this%trans%sub_type == 'maskvalid' .OR. this%trans%sub_type == 'maskinvalid') THEN
1672
1673 IF (.NOT.PRESENT(maskgrid)) THEN
1674 CALL l4f_category_log(this%category,l4f_error, &
1675 'grid_transform_init maskgrid argument missing for metamorphosis:'// &
1676 trim(this%trans%sub_type)//' transformation')
1677 CALL raise_error()
1678 RETURN
1679 ENDIF
1680
1681 IF (this%innx /= SIZE(maskgrid,1) .OR. this%inny /= SIZE(maskgrid,2)) THEN
1682 CALL l4f_category_log(this%category,l4f_error, &
1683 'grid_transform_init mask non conformal with input field')
1684 CALL l4f_category_log(this%category,l4f_error, &
1685 'mask: '//t2c(SIZE(maskgrid,1))//'x'//t2c(SIZE(maskgrid,2))// &
1686 ' input field:'//t2c(this%innx)//'x'//t2c(this%inny))
1687 CALL raise_error()
1688 RETURN
1689 ENDIF
1690
1691 ALLOCATE(this%point_mask(this%innx,this%inny))
1692
1693 IF (this%trans%sub_type == 'maskvalid') THEN
1694! behavior depends on the presence/usability of maskbounds,
1695! simplified wrt its use in metamorphosis:mask
1696 IF (.NOT.PRESENT(maskbounds)) THEN
1697 this%point_mask(:,:) = c_e(maskgrid(:,:))
1698 ELSE IF (SIZE(maskbounds) < 2) THEN
1699 this%point_mask(:,:) = c_e(maskgrid(:,:))
1700 ELSE
1701 this%point_mask(:,:) = c_e(maskgrid(:,:)) .AND. &
1702 maskgrid(:,:) > maskbounds(1) .AND. &
1703 maskgrid(:,:) <= maskbounds(SIZE(maskbounds))
1704 ENDIF
1705 ELSE ! reverse condition
1706 this%point_mask(:,:) = .NOT.c_e(maskgrid(:,:))
1707 ENDIF
1708
1709 this%valid = .true.
1710
1711 ELSE IF (this%trans%sub_type == 'lemaskinvalid' .OR. &
1712 this%trans%sub_type == 'ltmaskinvalid' .OR. &
1713 this%trans%sub_type == 'gemaskinvalid' .OR. &
1714 this%trans%sub_type == 'gtmaskinvalid') THEN
1715! here i can only store field for computing mask runtime
1716
1717 this%val_mask = maskgrid
1718 this%valid = .true.
1719
1720 ELSE IF (this%trans%sub_type == 'setinvalidto') THEN
1721
1722 IF (.NOT.PRESENT(maskbounds)) THEN
1723 CALL l4f_category_log(this%category,l4f_error, &
1724 'grid_transform_init maskbounds missing for metamorphosis:'// &
1725 trim(this%trans%sub_type)//' transformation')
1726 RETURN
1727 ELSE IF (SIZE(maskbounds) < 1) THEN
1728 CALL l4f_category_log(this%category,l4f_error, &
1729 'grid_transform_init maskbounds empty for metamorphosis:'// &
1730 trim(this%trans%sub_type)//' transformation')
1731 RETURN
1732 ELSE
1733 this%val1 = maskbounds(1)
1734#ifdef DEBUG
1735 CALL l4f_category_log(this%category, l4f_debug, &
1736 "grid_transform_init setting invalid data to "//t2c(this%val1))
1737#endif
1738 ENDIF
1739
1740 this%valid = .true.
1741
1742 ELSE IF (this%trans%sub_type == 'settoinvalid') THEN
1743
1744 IF (.NOT.PRESENT(maskbounds)) THEN
1745 CALL l4f_category_log(this%category,l4f_error, &
1746 'grid_transform_init maskbounds missing for metamorphosis:'// &
1747 trim(this%trans%sub_type)//' transformation')
1748 CALL raise_error()
1749 RETURN
1750 ELSE IF (SIZE(maskbounds) < 2) THEN
1751 CALL l4f_category_log(this%category,l4f_error, &
1752 'grid_transform_init maskbounds must have at least 2 elements for metamorphosis:'// &
1753 trim(this%trans%sub_type)//' transformation')
1754 CALL raise_error()
1755 RETURN
1756 ELSE
1757 this%val1 = maskbounds(1)
1758 this%val2 = maskbounds(SIZE(maskbounds))
1759#ifdef DEBUG
1760 CALL l4f_category_log(this%category, l4f_debug, &
1761 "grid_transform_init setting to invalid interval ]"//t2c(this%val1)//','// &
1762 t2c(this%val2)//']')
1763#endif
1764 ENDIF
1765
1766 this%valid = .true.
1767
1768 ENDIF
1769
1770ENDIF
1771
1772CONTAINS
1773
1774! local subroutine to be called by all methods interpolating to a new
1775! grid, no parameters passed, used as a macro to avoid repeating code
1776SUBROUTINE outgrid_setup()
1777
1778! set increments in new grid in order for all the baraque to work
1779CALL griddim_setsteps(out)
1780! check component flag
1781CALL get_val(in, proj=proj_in, component_flag=cf_i)
1782CALL get_val(out, proj=proj_out, component_flag=cf_o)
1783IF (proj_in == proj_out) THEN
1784! same projection: set output component flag equal to input regardless
1785! of its value
1786 CALL set_val(out, component_flag=cf_i)
1787ELSE
1788! different projection, interpolation possible only with vector data
1789! referred to geograpical axes
1790 IF (cf_i == 1) THEN
1791 CALL l4f_category_log(this%category,l4f_warn, &
1792 'trying to interpolate a grid with component flag 1 to a grid on a different projection')
1793 CALL l4f_category_log(this%category,L4F_WARN, &
1794 'vector fields will probably be wrong')
1795 ELSE
1796 CALL set_val(out, component_flag=cf_i)
1797 ENDIF
1798ENDIF
1799! rotate the input grid by n*360 degrees to bring it closer to the output grid
1800CALL griddim_set_central_lon(in, griddim_central_lon(out))
1801
1802END SUBROUTINE outgrid_setup
1803
1804END SUBROUTINE grid_transform_init
1805
1806
1807!> Constructor for a \a grid_transform object, defining a particular
1808!! grid-to-sparse points transformation.
1809!! It defines an object describing a transformation from a rectangular
1810!! grid to a set of sparse points; the abstract type of transformation
1811!! is described in the transformation object \a trans (type
1812!! transform_def) which must have been properly initialised. The
1813!! additional information required here is the description of the
1814!! input grid \a in (type griddim_def), and, if required by the
1815!! transformation type, the information about the target sparse points
1816!! over which the transformation should take place:
1817!!
1818!! - for 'inter' transformation, this is provided in the form of a
1819!! vol7d object (\a v7d_out argument, input), which must have been
1820!! initialized with the coordinates of desired sparse points
1821!!
1822!! - for 'polyinter' transformation, no target point information has
1823!! to be provided in input (it is calculated on the basis of input
1824!! grid and \a trans object), and the coordinates of the target
1825!! points (polygons' centroids) are returned in output in \a
1826!! v7d_out argument
1827!!
1828!! - for 'maskinter' transformation, this is a two dimensional real
1829!! field (\a maskgrid argument), which, together with the \a
1830!! maskbounds argument (optional with default), divides the input
1831!! grid in a number of subareas according to the values of \a
1832!! maskinter, and, in this case, \a v7d_out is an output argument
1833!! which returns the coordinates of the target points (subareas'
1834!! centroids)
1835!!
1836!! - for 'metamorphosis' transformation, no target point information
1837!! has to be provided in input (it is calculated on the basis of
1838!! input grid and \a trans object), except for 'mask' subtype, for
1839!! which the same information as for 'maskinter' transformation has
1840!! to be provided; in all the cases, as for 'polyinter', the
1841!! information about target points is returned in output in \a
1842!! v7d_out argument.
1843!!
1844!! The generated \a grid_transform object is specific to the grid and
1845!! sparse point list provided or computed. The function \a c_e can be
1846!! used in order to check whether the object has been successfully
1847!! initialised, if the result is \a .FALSE., it should not be used
1848!! further on.
1849SUBROUTINE grid_transform_grid_vol7d_init(this, trans, in, v7d_out, &
1850 maskgrid, maskbounds, find_index, categoryappend)
1851TYPE(grid_transform),INTENT(out) :: this
1852TYPE(transform_def),INTENT(in) :: trans
1853TYPE(griddim_def),INTENT(in) :: in
1854TYPE(vol7d),INTENT(inout) :: v7d_out
1855REAL,INTENT(in),OPTIONAL :: maskgrid(:,:)
1856REAL,INTENT(in),OPTIONAL :: maskbounds(:)
1857PROCEDURE(basic_find_index),POINTER,OPTIONAL :: find_index
1858CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
1859
1860INTEGER :: ix, iy, n, nm, nr, nprev, nmaskarea, xnmin, xnmax, ynmin, ynmax, &
1861 time_definition
1862DOUBLE PRECISION :: xmin, xmax, ymin, ymax, r2, lonref
1863DOUBLE PRECISION,ALLOCATABLE :: lon1(:), lat1(:), lon(:,:), lat(:,:)
1864REAL,ALLOCATABLE :: lmaskbounds(:)
1865TYPE(georef_coord) :: point
1866TYPE(griddim_def) :: lin
1867
1868
1869IF (PRESENT(find_index)) THEN ! move in init_common?
1870 IF (ASSOCIATED(find_index)) THEN
1871 this%find_index => find_index
1872 ENDIF
1873ENDIF
1874CALL grid_transform_init_common(this, trans, categoryappend)
1875#ifdef DEBUG
1876CALL l4f_category_log(this%category, L4F_DEBUG, "grid_transform vg6d-v7d")
1877#endif
1878
1879! used after in some transformations
1880CALL get_val(trans, time_definition=time_definition)
1881IF (.NOT. c_e(time_definition)) THEN
1882 time_definition=1 ! default to validity time
1883ENDIF
1884
1885IF (this%trans%trans_type == 'inter') THEN
1886
1887 IF (this%trans%sub_type == 'near' .OR. this%trans%sub_type == 'bilin' &
1888 .OR. this%trans%sub_type == 'shapiro_near') THEN
1889
1890 CALL copy(in, lin)
1891 CALL get_val(lin, nx=this%innx, ny=this%inny)
1892 this%outnx = SIZE(v7d_out%ana)
1893 this%outny = 1
1894
1895 ALLOCATE (this%inter_index_x(this%outnx,this%outny),&
1896 this%inter_index_y(this%outnx,this%outny))
1897 ALLOCATE(lon(this%outnx,1),lat(this%outnx,1))
1898
1899 CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
1900! equalize in/out coordinates
1901 lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1), mask=c_e(lon(:,1))))
1902 CALL griddim_set_central_lon(lin, lonref)
1903 CALL get_val(lin, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1904
1905 IF (this%trans%sub_type == 'bilin') THEN
1906 CALL this%find_index(lin, .false., &
1907 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1908 lon, lat, this%trans%extrap, &
1909 this%inter_index_x, this%inter_index_y)
1910
1911 ALLOCATE(this%inter_x(this%innx,this%inny),this%inter_y(this%innx,this%inny))
1912 ALLOCATE(this%inter_xp(this%outnx,this%outny),this%inter_yp(this%outnx,this%outny))
1913
1914 CALL griddim_gen_coord(lin, this%inter_x, this%inter_y)
1915 CALL proj(lin, lon, lat, this%inter_xp, this%inter_yp)
1916
1917 ELSE ! near shapiro_near
1918 CALL this%find_index(lin, .true., &
1919 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1920 lon, lat, this%trans%extrap, &
1921 this%inter_index_x, this%inter_index_y)
1922
1923 IF (this%trans%trans_type == 'intersearch') THEN ! replicate code above
1924 ALLOCATE(this%inter_x(this%innx,this%inny),this%inter_y(this%innx,this%inny))
1925 ALLOCATE(this%inter_xp(this%outnx,this%outny),this%inter_yp(this%outnx,this%outny))
1926
1927 CALL griddim_gen_coord(lin, this%inter_x, this%inter_y)
1928 CALL proj(lin, lon, lat, this%inter_xp, this%inter_yp)
1929 ENDIF
1930 ENDIF
1931
1932 DEALLOCATE(lon,lat)
1933 CALL delete(lin)
1934
1935 this%valid = .true.
1936
1937 ENDIF
1938
1939ELSE IF (this%trans%trans_type == 'polyinter') THEN
1940
1941 CALL get_val(in, nx=this%innx, ny=this%inny)
1942! unlike before, here index arrays must have the shape of input grid
1943 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1944 this%inter_index_y(this%innx,this%inny))
1945 this%inter_index_x(:,:) = imiss
1946 this%inter_index_y(:,:) = 1
1947
1948! compute coordinates of input grid in geo system
1949 CALL copy(in, lin)
1950 CALL unproj(lin)
1951
1952 this%outnx = this%trans%poly%arraysize
1953 this%outny = 1
1954 CALL delete(v7d_out) ! required to avoid leaks because intent(inout), dangerous
1955 CALL init(v7d_out, time_definition=time_definition)
1956 CALL vol7d_alloc(v7d_out, nana=this%outnx)
1957
1958! equalize in/out coordinates
1959 ALLOCATE(lon(this%outnx,1))
1960 CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1))
1961 lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1), mask=c_e(lon(:,1))))
1962 CALL griddim_set_central_lon(lin, lonref)
1963 DEALLOCATE(lon)
1964
1965! setup output point list, equal to average of polygon points
1966! warning, in case of concave areas points may coincide!
1967 CALL poly_to_coordinates(this%trans%poly, v7d_out)
1968
1969 nprev = 1
1970!$OMP PARALLEL DEFAULT(SHARED)
1971!$OMP DO PRIVATE(iy, ix, n, point) FIRSTPRIVATE(nprev)
1972 DO iy = 1, this%inny
1973 inside_x: DO ix = 1, this%innx
1974 point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy))
1975
1976 DO n = nprev, this%trans%poly%arraysize ! optimize starting from last matched polygon
1977 IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
1978 this%inter_index_x(ix,iy) = n
1979 nprev = n
1980 cycle inside_x
1981 ENDIF
1982 ENDDO
1983 DO n = nprev-1, 1, -1 ! test the other polygons
1984 IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
1985 this%inter_index_x(ix,iy) = n
1986 nprev = n
1987 cycle inside_x
1988 ENDIF
1989 ENDDO
1990
1991! CALL delete(point) ! speedup

Generated with Doxygen.