libsim Versione 7.2.0

◆ 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 
)

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

Generated with Doxygen.