libsim  Versione7.2.3

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

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 
1455 ELSE 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 
1486 ELSE 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 
1544 ELSE 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 
1624 ELSE 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 
1672 ELSE 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 
1778 ENDIF
1779 
1780 CONTAINS
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
1784 SUBROUTINE outgrid_setup()
1785 
1786 ! set increments in new grid in order for all the baraque to work
1787 CALL griddim_setsteps(out)
1788 ! check component flag
1789 CALL get_val(in, proj=proj_in, component_flag=cf_i)
1790 CALL get_val(out, proj=proj_out, component_flag=cf_o)
1791 IF (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)
1795 ELSE
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
1806 ENDIF
1807 ! rotate the input grid by n*360 degrees to bring it closer to the output grid
1808 CALL griddim_set_central_lon(in, griddim_central_lon(out))
1809 
1810 END SUBROUTINE outgrid_setup
1811 
1812 END SUBROUTINE grid_transform_init
1813 
1814 
1857 SUBROUTINE grid_transform_grid_vol7d_init(this, trans, in, v7d_out, &
1858  maskgrid, maskbounds, find_index, categoryappend)
1859 TYPE(grid_transform),INTENT(out) :: this
1860 TYPE(transform_def),INTENT(in) :: trans
1861 TYPE(griddim_def),INTENT(in) :: in
1862 TYPE(vol7d),INTENT(inout) :: v7d_out
1863 REAL,INTENT(in),OPTIONAL :: maskgrid(:,:)
1864 REAL,INTENT(in),OPTIONAL :: maskbounds(:)
1865 PROCEDURE(basic_find_index),POINTER,OPTIONAL :: find_index
1866 CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
1867 
1868 INTEGER :: ix, iy, n, nm, nr, nprev, nmaskarea, xnmin, xnmax, ynmin, ynmax, &
1869  time_definition
1870 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, r2, lonref
1871 DOUBLE PRECISION,ALLOCATABLE :: lon1(:), lat1(:), lon(:,:), lat(:,:)
1872 REAL,ALLOCATABLE :: lmaskbounds(:)
1873 TYPE(georef_coord) :: point
1874 TYPE(griddim_def) :: lin
1875 
1876 
1877 IF (PRESENT(find_index)) THEN ! move in init_common?
1878  IF (ASSOCIATED(find_index)) THEN
1879  this%find_index => find_index
1880  ENDIF
1881 ENDIF
1882 CALL grid_transform_init_common(this, trans, categoryappend)
1883 #ifdef DEBUG
1884 CALL l4f_category_log(this%category, l4f_debug, "grid_transform vg6d-v7d")
1885 #endif
1886 
1887 ! used after in some transformations
1888 CALL get_val(trans, time_definition=time_definition)
1889 IF (.NOT. c_e(time_definition)) THEN
1890  time_definition=1 ! default to validity time
1891 ENDIF
1892 
1893 IF (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 
1947 ELSE 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
2000  ENDDO inside_x
2001  ENDDO
2002 !$OMP END PARALLEL
2003 
2004 #ifdef DEBUG
2005  DO n = 1, this%trans%poly%arraysize
2006  CALL l4f_category_log(this%category, l4f_debug, &
2007  'Polygon: '//t2c(n)//' grid points: '// &
2008  t2c(count(this%inter_index_x(:,:) == n)))
2009  ENDDO

Generated with Doxygen.