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 
)

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

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 
1447 ELSE 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 
1478 ELSE 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 
1536 ELSE 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 
1616 ELSE 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 
1664 ELSE 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 
1770 ENDIF
1771 
1772 CONTAINS
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
1776 SUBROUTINE outgrid_setup()
1777 
1778 ! set increments in new grid in order for all the baraque to work
1779 CALL griddim_setsteps(out)
1780 ! check component flag
1781 CALL get_val(in, proj=proj_in, component_flag=cf_i)
1782 CALL get_val(out, proj=proj_out, component_flag=cf_o)
1783 IF (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)
1787 ELSE
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
1798 ENDIF
1799 ! rotate the input grid by n*360 degrees to bring it closer to the output grid
1800 CALL griddim_set_central_lon(in, griddim_central_lon(out))
1801 
1802 END SUBROUTINE outgrid_setup
1803 
1804 END SUBROUTINE grid_transform_init
1805 
1806 
1849 SUBROUTINE grid_transform_grid_vol7d_init(this, trans, in, v7d_out, &
1850  maskgrid, maskbounds, find_index, categoryappend)
1851 TYPE(grid_transform),INTENT(out) :: this
1852 TYPE(transform_def),INTENT(in) :: trans
1853 TYPE(griddim_def),INTENT(in) :: in
1854 TYPE(vol7d),INTENT(inout) :: v7d_out
1855 REAL,INTENT(in),OPTIONAL :: maskgrid(:,:)
1856 REAL,INTENT(in),OPTIONAL :: maskbounds(:)
1857 PROCEDURE(basic_find_index),POINTER,OPTIONAL :: find_index
1858 CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
1859 
1860 INTEGER :: ix, iy, n, nm, nr, nprev, nmaskarea, xnmin, xnmax, ynmin, ynmax, &
1861  time_definition
1862 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, r2, lonref
1863 DOUBLE PRECISION,ALLOCATABLE :: lon1(:), lat1(:), lon(:,:), lat(:,:)
1864 REAL,ALLOCATABLE :: lmaskbounds(:)
1865 TYPE(georef_coord) :: point
1866 TYPE(griddim_def) :: lin
1867 
1868 
1869 IF (PRESENT(find_index)) THEN ! move in init_common?
1870  IF (ASSOCIATED(find_index)) THEN
1871  this%find_index => find_index
1872  ENDIF
1873 ENDIF
1874 CALL grid_transform_init_common(this, trans, categoryappend)
1875 #ifdef DEBUG
1876 CALL l4f_category_log(this%category, l4f_debug, "grid_transform vg6d-v7d")
1877 #endif
1878 
1879 ! used after in some transformations
1880 CALL get_val(trans, time_definition=time_definition)
1881 IF (.NOT. c_e(time_definition)) THEN
1882  time_definition=1 ! default to validity time
1883 ENDIF
1884 
1885 IF (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 
1939 ELSE 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
1992  ENDDO inside_x
1993  ENDDO
1994 !$OMP END PARALLEL
1995 
1996 #ifdef DEBUG
1997  DO n = 1, this%trans%poly%arraysize
1998  CALL l4f_category_log(this%category, l4f_debug, &
1999  'Polygon: '//t2c(n)//' grid points: '// &
2000  t2c(count(this%inter_index_x(:,:) == n)))
2001  ENDDO

Generated with Doxygen.