libsim  Versione6.3.0

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

1416  CALL get_val(out, dx=this%trans%area_info%boxdx)
1417  IF (.NOT.c_e(this%trans%area_info%boxdy)) &
1418  CALL get_val(out, dx=this%trans%area_info%boxdy)
1419 ! half size is actually needed
1420  this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
1421  this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
1422 ! unlike before, here index arrays must have the shape of input grid
1423  ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1424  this%inter_index_y(this%innx,this%inny))
1425 
1426 ! compute coordinates of input grid in geo system
1427  CALL unproj(in) ! TODO costringe a dichiarare in INTENT(inout), si puo` evitare?
1428 ! use find_index in the opposite way, here extrap does not make sense
1429  CALL find_index(out, 'near', &
1430  this%outnx, this%outny, xmin, xmax, ymin, ymax, &
1431  in%dim%lon, in%dim%lat, .false., &
1432  this%inter_index_x, this%inter_index_y)
1433 
1434  this%valid = .true. ! warning, no check of subtype
1435 
1436 ELSE IF (this%trans%trans_type == 'stencilinter') THEN
1437 
1438  CALL outgrid_setup() ! common setup for grid-generating methods
1439 ! from inter:near
1440  CALL get_val(in, nx=this%innx, ny=this%inny, &
1441  xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1442  CALL get_val(out, nx=this%outnx, ny=this%outny)
1443 
1444  ALLOCATE (this%inter_index_x(this%outnx,this%outny), &
1445  this%inter_index_y(this%outnx,this%outny))
1446  CALL unproj(out)
1447  CALL find_index(in, 'near', &
1448  this%innx, this%inny, xmin, xmax, ymin, ymax, &
1449  out%dim%lon, out%dim%lat, this%trans%extrap, &
1450  this%inter_index_x, this%inter_index_y)
1451 
1452 ! define the stencil mask
1453  nr = int(this%trans%area_info%radius) ! integer radius
1454  n = nr*2+1 ! n. of points
1455  nm = nr + 1 ! becomes index of center
1456  r2 = this%trans%area_info%radius**2
1457  ALLOCATE(this%stencil(n,n))
1458  this%stencil(:,:) = .true.
1459  DO iy = 1, n
1460  DO ix = 1, n
1461  IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
1462  ENDDO
1463  ENDDO
1464 
1465 ! set to missing the elements of inter_index too close to the domain
1466 ! borders
1467  xnmin = nr + 1
1468  xnmax = this%innx - nr
1469  ynmin = nr + 1
1470  ynmax = this%inny - nr
1471  DO iy = 1, this%outny
1472  DO ix = 1, this%outnx
1473  IF (this%inter_index_x(ix,iy) < xnmin .OR. &
1474  this%inter_index_x(ix,iy) > xnmax .OR. &
1475  this%inter_index_y(ix,iy) < ynmin .OR. &
1476  this%inter_index_y(ix,iy) > ynmax) THEN
1477  this%inter_index_x(ix,iy) = imiss
1478  this%inter_index_y(ix,iy) = imiss
1479  ENDIF
1480  ENDDO
1481  ENDDO
1482 
1483 #ifdef DEBUG
1484  CALL l4f_category_log(this%category, l4f_debug, &
1485  'stencilinter: stencil size '//t2c(n*n))
1486  CALL l4f_category_log(this%category, l4f_debug, &
1487  'stencilinter: stencil points '//t2c(count(this%stencil)))
1488 #endif
1489 
1490  this%valid = .true. ! warning, no check of subtype
1491 
1492 ELSE IF (this%trans%trans_type == 'maskgen' .OR. &
1493  this%trans%trans_type == 'polyinter') THEN
1494 
1495  IF (this%trans%trans_type == 'polyinter') THEN
1496  this%recur = .true. ! grid-to-grid polyinter is done in two steps!
1497  ENDIF
1498 
1499  CALL copy(in,out)
1500  CALL get_val(in, nx=this%innx, ny=this%inny)
1501  this%outnx = this%innx
1502  this%outny = this%inny
1503 
1504 ! unlike before, here index arrays must have the shape of input grid
1505  ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1506  this%inter_index_y(this%innx,this%inny))
1507  this%inter_index_x(:,:) = imiss
1508  this%inter_index_y(:,:) = 1
1509 
1510 ! compute coordinates of input grid in geo system
1511  CALL unproj(in) ! TODO costringe a dichiarare in INTENT(inout), si puo` evitare?
1512 
1513  nprev = 1
1514 !$OMP PARALLEL DEFAULT(SHARED)
1515 !$OMP DO PRIVATE(j, i, n, point) FIRSTPRIVATE(nprev)
1516  DO j = 1, this%inny
1517  inside_x: DO i = 1, this%innx
1518  point = georef_coord_new(x=in%dim%lon(i,j), y=in%dim%lat(i,j))
1519 
1520  DO n = nprev, this%trans%poly%arraysize ! optimize starting from last matched polygon
1521  IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
1522  this%inter_index_x(i,j) = n
1523  nprev = n
1524  cycle inside_x
1525  ENDIF
1526  ENDDO
1527  DO n = nprev-1, 1, -1 ! test the other polygons
1528  IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
1529  this%inter_index_x(i,j) = n
1530  nprev = n
1531  cycle inside_x
1532  ENDIF
1533  ENDDO
1534 
1535 ! CALL delete(point) ! speedup
1536  ENDDO inside_x
1537  ENDDO
1538 !$OMP END PARALLEL
1539 
1540  this%valid = .true. ! warning, no check of subtype
1541 
1542 ELSE IF (this%trans%trans_type == 'metamorphosis') THEN
1543 
1544  CALL copy(in,out)
1545  CALL get_val(in, nx=this%innx, ny=this%inny)
1546  this%outnx = this%innx
1547  this%outny = this%inny
1548 
1549  IF (this%trans%sub_type == 'maskvalid' .OR. this%trans%sub_type == 'maskinvalid') THEN
1550 
1551  IF (.NOT.PRESENT(maskgrid)) THEN
1552  CALL l4f_category_log(this%category,l4f_error, &
1553  'grid_transform_init maskgrid argument missing for metamorphosis:'// &
1554  trim(this%trans%sub_type)//' transformation')
1555  CALL raise_error()
1556  RETURN
1557  ENDIF
1558 
1559  IF (this%innx /= SIZE(maskgrid,1) .OR. this%inny /= SIZE(maskgrid,2)) THEN
1560  CALL l4f_category_log(this%category,l4f_error, &
1561  'grid_transform_init mask non conformal with input field')
1562  CALL l4f_category_log(this%category,l4f_error, &
1563  'mask: '//t2c(SIZE(maskgrid,1))//'x'//t2c(SIZE(maskgrid,2))// &
1564  ' input field:'//t2c(this%innx)//'x'//t2c(this%inny))
1565  CALL raise_error()
1566  RETURN
1567  ENDIF
1568 
1569  ALLOCATE(this%point_mask(this%innx,this%inny))
1570 
1571  IF (this%trans%sub_type == 'maskvalid') THEN
1572 ! behavior depends on the presence/usability of maskbounds,
1573 ! simplified wrt its use in metamorphosis:mask
1574  IF (.NOT.PRESENT(maskbounds)) THEN
1575  this%point_mask(:,:) = c_e(maskgrid(:,:))
1576  ELSE IF (SIZE(maskbounds) < 2) THEN
1577  this%point_mask(:,:) = c_e(maskgrid(:,:))
1578  ELSE
1579  this%point_mask(:,:) = c_e(maskgrid(:,:)) .AND. &
1580  maskgrid(:,:) > maskbounds(1) .AND. &
1581  maskgrid(:,:) <= maskbounds(SIZE(maskbounds))
1582  ENDIF
1583  ELSE ! reverse condition
1584  this%point_mask(:,:) = .NOT.c_e(maskgrid(:,:))
1585  ENDIF
1586 
1587  this%valid = .true.
1588 
1589  ELSE IF (this%trans%sub_type == 'setinvalidto') THEN
1590 
1591  IF (.NOT.PRESENT(maskbounds)) THEN
1592  CALL l4f_category_log(this%category,l4f_error, &
1593  'grid_transform_init maskbounds missing for metamorphosis:'// &
1594  trim(this%trans%sub_type)//' transformation')
1595  RETURN
1596  ELSE IF (SIZE(maskbounds) < 1) THEN
1597  CALL l4f_category_log(this%category,l4f_error, &
1598  'grid_transform_init maskbounds empty for metamorphosis:'// &
1599  trim(this%trans%sub_type)//' transformation')
1600  RETURN
1601  ELSE
1602  this%val1 = maskbounds(1)
1603 #ifdef DEBUG
1604  CALL l4f_category_log(this%category, l4f_debug, &
1605  "grid_transform_init setting invalid data to "//t2c(this%val1))
1606 #endif
1607  ENDIF
1608 
1609  this%valid = .true.
1610 
1611  ELSE IF (this%trans%sub_type == 'settoinvalid') THEN
1612 
1613  IF (.NOT.PRESENT(maskbounds)) THEN
1614  CALL l4f_category_log(this%category,l4f_error, &
1615  'grid_transform_init maskbounds missing for metamorphosis:'// &
1616  trim(this%trans%sub_type)//' transformation')
1617  CALL raise_error()
1618  RETURN
1619  ELSE IF (SIZE(maskbounds) < 2) THEN
1620  CALL l4f_category_log(this%category,l4f_error, &
1621  'grid_transform_init maskbounds must have at least 2 elements for metamorphosis:'// &
1622  trim(this%trans%sub_type)//' transformation')
1623  CALL raise_error()
1624  RETURN
1625  ELSE
1626  this%val1 = maskbounds(1)
1627  this%val2 = maskbounds(SIZE(maskbounds))
1628 #ifdef DEBUG
1629  CALL l4f_category_log(this%category, l4f_debug, &
1630  "grid_transform_init setting to invalid interval ]"//t2c(this%val1)//','// &
1631  t2c(this%val2)//']')
1632 #endif
1633  ENDIF
1634 
1635  this%valid = .true.
1636 
1637  ENDIF
1638 
1639 ENDIF
1640 
1641 CONTAINS
1642 
1643 ! local subroutine to be called by all methods interpolating to a new
1644 ! grid, no parameters passed, used as a macro to avoid repeating code
1645 SUBROUTINE outgrid_setup()
1646 
1647 ! set increments in new grid in order for all the baraque to work
1648 CALL griddim_setsteps(out)
1649 ! check component flag
1650 CALL get_val(in, proj=proj_in, component_flag=cf_i)
1651 CALL get_val(out, proj=proj_out, component_flag=cf_o)
1652 IF (proj_in == proj_out) THEN
1653 ! same projection: set output component flag equal to input regardless
1654 ! of its value
1655  CALL set_val(out, component_flag=cf_i)
1656 ELSE
1657 ! different projection, interpolation possible only with vector data
1658 ! referred to geograpical axes
1659  IF (cf_i == 1) THEN
1660  CALL l4f_category_log(this%category,l4f_warn, &
1661  'trying to interpolate a grid with component flag 1 to a grid on a different projection')
1662  CALL l4f_category_log(this%category,l4f_warn, &
1663  'vector fields will probably be wrong')
1664  ELSE
1665  CALL set_val(out, component_flag=cf_i)
1666  ENDIF
1667 ENDIF
1668 ! rotate the input grid by n*360 degrees to bring it closer to the output grid
1669 CALL griddim_set_central_lon(in, griddim_central_lon(out))
1670 
1671 END SUBROUTINE outgrid_setup
1672 
1673 END SUBROUTINE grid_transform_init
1674 
1675 
1718 SUBROUTINE grid_transform_grid_vol7d_init(this, trans, in, v7d_out, &
1719  maskgrid, maskbounds, categoryappend)
1720 TYPE(grid_transform),INTENT(out) :: this
1721 TYPE(transform_def),INTENT(in) :: trans
1722 TYPE(griddim_def),INTENT(inout) :: in
1723 TYPE(vol7d),INTENT(inout) :: v7d_out
1724 REAL,INTENT(in),OPTIONAL :: maskgrid(:,:)
1725 REAL,INTENT(in),OPTIONAL :: maskbounds(:)
1726 CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
1727 
1728 INTEGER :: ix, iy, n, nm, nr, nprev, nmaskarea, xnmin, xnmax, ynmin, ynmax, &
1729  time_definition
1730 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, r2
1731 DOUBLE PRECISION,ALLOCATABLE :: lon(:), lat(:)
1732 REAL,ALLOCATABLE :: lmaskbounds(:)
1733 TYPE(georef_coord) :: point
1734 
1735 
1736 CALL grid_transform_init_common(this, trans, categoryappend)
1737 #ifdef DEBUG
1738 CALL l4f_category_log(this%category, l4f_debug, "grid_transform vg6d-v7d")
1739 #endif
1740 
1741 ! used after in some transformations
1742 CALL get_val(trans, time_definition=time_definition)
1743 IF (.NOT. c_e(time_definition)) THEN
1744  time_definition=1 ! default to validity time
1745 ENDIF
1746 
1747 IF (this%trans%trans_type == 'inter') THEN
1748 
1749  IF (this%trans%sub_type == 'near' .OR. this%trans%sub_type == 'bilin' &
1750  .OR. this%trans%sub_type == 'shapiro_near') THEN
1751 
1752  CALL get_val(in, nx=this%innx, ny=this%inny)
1753  this%outnx = SIZE(v7d_out%ana)
1754  this%outny = 1
1755 
1756  ALLOCATE (this%inter_index_x(this%outnx,this%outny),&
1757  this%inter_index_y(this%outnx,this%outny))
1758  ALLOCATE(lon(this%outnx),lat(this%outnx))
1759 
1760  CALL get_val(in, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1761  CALL getval(v7d_out%ana(:)%coord,lon=lon,lat=lat)
1762 
1763  CALL find_index(in, this%trans%sub_type,&
1764  this%innx, this%inny, xmin, xmax, ymin, ymax, &
1765  lon, lat, this%trans%extrap, &
1766  this%inter_index_x(:,1), this%inter_index_y(:,1))
1767 
1768  IF ( this%trans%sub_type == 'bilin' ) THEN
1769  ALLOCATE(this%inter_x(this%innx,this%inny),this%inter_y(this%innx,this%inny))
1770  ALLOCATE(this%inter_xp(this%outnx,this%outny),this%inter_yp(this%outnx,this%outny))
1771 
1772  CALL griddim_gen_coord(in, this%inter_x, this%inter_y)
1773  CALL proj(in, reshape(lon,(/SIZE(lon),1/)),reshape(lat,(/SIZE(lat),1/)),&
1774  this%inter_xp,this%inter_yp)
1775  ENDIF
1776 
1777  DEALLOCATE(lon,lat)
1778 
1779  this%valid = .true.
1780 
1781  ENDIF
1782 
1783 ELSE IF (this%trans%trans_type == 'polyinter') THEN
1784 
1785  CALL get_val(in, nx=this%innx, ny=this%inny)
1786 ! unlike before, here index arrays must have the shape of input grid
1787  ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1788  this%inter_index_y(this%innx,this%inny))
1789  this%inter_index_x(:,:) = imiss
1790  this%inter_index_y(:,:) = 1
1791 
1792 ! compute coordinates of input grid in geo system
1793  CALL unproj(in) ! TODO costringe a dichiarare in INTENT(inout), si puo` evitare?
1794 
1795  nprev = 1
1796 !$OMP PARALLEL DEFAULT(SHARED)
1797 !$OMP DO PRIVATE(iy, ix, n, point) FIRSTPRIVATE(nprev)
1798  DO iy = 1, this%inny
1799  inside_x: DO ix = 1, this%innx
1800  point = georef_coord_new(x=in%dim%lon(ix,iy), y=in%dim%lat(ix,iy))
1801 
1802  DO n = nprev, this%trans%poly%arraysize ! optimize starting from last matched polygon
1803  IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
1804  this%inter_index_x(ix,iy) = n
1805  nprev = n
1806  cycle inside_x
1807  ENDIF
1808  ENDDO
1809  DO n = nprev-1, 1, -1 ! test the other polygons
1810  IF (inside(point, this%trans%poly%array(n))) THEN ! stop at the first matching polygon
1811  this%inter_index_x(ix,iy) = n
1812  nprev = n
1813  cycle inside_x
1814  ENDIF
1815  ENDDO
1816 
1817 ! CALL delete(point) ! speedup
1818  ENDDO inside_x
1819  ENDDO
1820 !$OMP END PARALLEL
1821 
1822  this%outnx = this%trans%poly%arraysize
1823  this%outny = 1
1824  CALL delete(v7d_out) ! required to avoid leaks because intent(inout), dangerous
1825  CALL init(v7d_out, time_definition=time_definition)
1826  CALL vol7d_alloc(v7d_out, nana=this%outnx)
1827 
1828 ! setup output point list, equal to average of polygon points
1829 ! warning, in case of concave areas points may coincide!
1830  DO n = 1, this%trans%poly%arraysize
1831  CALL getval(this%trans%poly%array(n), x=lon, y=lat)
1832  CALL init(v7d_out%ana(n), lon=stat_average(lon), lat=stat_average(lat))
1833 ! DEALLOCATE(lon, lat)
1834  ENDDO
1835 
1836 #ifdef DEBUG
1837  DO n = 1, this%trans%poly%arraysize
1838  CALL l4f_category_log(this%category, l4f_debug, &
1839  'Polygon: '//t2c(n)//' grid points: '// &
1840  t2c(count(this%inter_index_x(:,:) == n)))
1841  ENDDO
1842 #endif
1843 
1844  this%valid = .true. ! warning, no check of subtype
1845 
1846 ELSE IF (this%trans%trans_type == 'stencilinter') THEN
1847 
1848 ! from inter:near
1849  CALL get_val(in, nx=this%innx, ny=this%inny)
1850  this%outnx = SIZE(v7d_out%ana)
1851  this%outny = 1
1852 
1853  ALLOCATE (this%inter_index_x(this%outnx,this%outny),&
1854  this%inter_index_y(this%outnx,this%outny))
1855  ALLOCATE(lon(this%outnx),lat(this%outnx))
1856 
1857  CALL get_val(in, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1858  CALL getval(v7d_out%ana(:)%coord,lon=lon,lat=lat)
1859 
1860  CALL find_index(in, 'near',&
1861  this%innx, this%inny, xmin, xmax, ymin, ymax, &
1862  lon, lat, this%trans%extrap, &
1863  this%inter_index_x(:,1), this%inter_index_y(:,1))
1864 
1865 ! define the stencil mask
1866  nr = int(this%trans%area_info%radius) ! integer radius
1867  n = nr*2+1 ! n. of points
1868  nm = nr + 1 ! becomes index of center
1869  r2 = this%trans%area_info%radius**2
1870  ALLOCATE(this%stencil(n,n))
1871  this%stencil(:,:) = .true.
1872  DO iy = 1, n
1873  DO ix = 1, n
1874  IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
1875  ENDDO
1876  ENDDO
1877 
1878 ! set to missing the elements of inter_index too close to the domain
1879 ! borders
1880  xnmin = nr + 1
1881  xnmax = this%innx - nr
1882  ynmin = nr + 1
1883  ynmax = this%inny - nr
1884  DO iy = 1, this%outny
1885  DO ix = 1, this%outnx
Functions that return a trimmed CHARACTER representation of the input variable.
Distruttori per le 2 classi.
Restituiscono il valore dell&#39;oggetto nella forma desiderata.
Costruttori per le classi datetime e timedelta.

Generated with Doxygen.