libsim  Versione 7.1.9

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

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 
1436 ELSE 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 
1467 ELSE 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 
1525 ELSE 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 
1605 ELSE 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 
1653 ELSE 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 
1759 ENDIF
1760 
1761 CONTAINS
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
1765 SUBROUTINE outgrid_setup()
1766 
1767 ! set increments in new grid in order for all the baraque to work
1768 CALL griddim_setsteps(out)
1769 ! check component flag
1770 CALL get_val(in, proj=proj_in, component_flag=cf_i)
1771 CALL get_val(out, proj=proj_out, component_flag=cf_o)
1772 IF (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)
1776 ELSE
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
1787 ENDIF
1788 ! rotate the input grid by n*360 degrees to bring it closer to the output grid
1789 CALL griddim_set_central_lon(in, griddim_central_lon(out))
1790 
1791 END SUBROUTINE outgrid_setup
1792 
1793 END 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.
1838 SUBROUTINE grid_transform_grid_vol7d_init(this, trans, in, v7d_out, &
1839  maskgrid, maskbounds, find_index, categoryappend)
1840 TYPE(grid_transform),INTENT(out) :: this
1841 TYPE(transform_def),INTENT(in) :: trans
1842 TYPE(griddim_def),INTENT(in) :: in
1843 TYPE(vol7d),INTENT(inout) :: v7d_out
1844 REAL,INTENT(in),OPTIONAL :: maskgrid(:,:)
1845 REAL,INTENT(in),OPTIONAL :: maskbounds(:)
1846 PROCEDURE(basic_find_index),POINTER,OPTIONAL :: find_index
1847 CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
1848 
1849 INTEGER :: ix, iy, n, nm, nr, nprev, nmaskarea, xnmin, xnmax, ynmin, ynmax, &
1850  time_definition
1851 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, r2, lonref
1852 DOUBLE PRECISION,ALLOCATABLE :: lon1(:), lat1(:), lon(:,:), lat(:,:)
1853 REAL,ALLOCATABLE :: lmaskbounds(:)
1854 TYPE(georef_coord) :: point
1855 TYPE(griddim_def) :: lin
1856 
1857 
1858 IF (PRESENT(find_index)) THEN ! move in init_common?
1859  IF (ASSOCIATED(find_index)) THEN
1860  this%find_index => find_index
1861  ENDIF
1862 ENDIF
1863 CALL grid_transform_init_common(this, trans, categoryappend)
1864 #ifdef DEBUG
1865 CALL l4f_category_log(this%category, L4F_DEBUG, "grid_transform vg6d-v7d")
1866 #endif
1867 
1868 ! used after in some transformations
1869 CALL get_val(trans, time_definition=time_definition)
1870 IF (.NOT. c_e(time_definition)) THEN
1871  time_definition=1 ! default to validity time
1872 ENDIF
1873 
1874 IF (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 
1921 ELSE 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, &
1981  'polygon: '//t2c(n)//' grid points: '// &
1982  t2c(COUNT(this%inter_index_x(:,:) == n)))
1983  ENDDO
1984 #endif
1985 
1986  CALL delete(lin)

Generated with Doxygen.