libsim  Versione 7.1.7

◆ grid_transform_init()

subroutine grid_transform_class::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 1408 del file grid_transform_class.F90.

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

Generated with Doxygen.