libsim Versione 7.2.4

◆ 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 1417 del file grid_transform_class.F90.

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

Generated with Doxygen.