|
◆ 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] | this | grid_transformation object |
[in] | trans | transformation object |
[in,out] | in | griddim object to transform |
[in,out] | out | griddim object defining target grid (input or output depending on type of transformation) |
[in] | maskgrid | 2D 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] | maskbounds | array 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] | categoryappend | append this suffix to log4fortran namespace category |
Definizione alla linea 1408 del file grid_transform_class.F90.
1410 this%inter_yp(this%outnx,this%outny))
1413 CALL griddim_gen_coord(in, this%inter_x, this%inter_y)
1415 CALL proj(in, lout%dim%lon, lout%dim%lat, this%inter_xp, this%inter_yp
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)
1429 ELSE IF (this%trans%trans_type == 'boxinter') THEN
1431 CALL outgrid_setup()
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)
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)
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
1445 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1446 this%inter_index_y(this%innx,this%inny))
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)
1460 ELSE IF (this%trans%trans_type == 'stencilinter') THEN
1462 CALL outgrid_setup()
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)
1468 ALLOCATE (this%inter_index_x(this%outnx,this%outny), &
1469 this%inter_index_y(this%outnx,this%outny))
1470 CALL copy(out, 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)
1478 nr = int(this%trans%area_info%radius)
1481 r2 = this%trans%area_info%radius**2
1482 ALLOCATE(this%stencil(n,n))
1483 this%stencil(:,:) = .true.
1486 IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
1493 xnmax = this%innx - nr
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
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)))
1518 ELSE IF (this%trans%trans_type == 'maskgen') THEN
1520 IF (this%trans%sub_type == 'poly') THEN
1523 CALL get_val(in, nx=this%innx, ny=this%inny)
1524 this%outnx = this%innx
1525 this%outny = this%inny
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
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))
1543 DO n = nprev, this%trans%poly%arraysize
1544 IF (inside(point, this%trans%poly%array(n))) THEN
1545 this%inter_index_x(i,j) = n
1550 DO n = nprev-1, 1, -1
1551 IF (inside(point, this%trans%poly%array(n))) THEN
1552 this%inter_index_x(i,j) = n
1563 ELSE IF (this%trans%sub_type == 'grid') THEN
1566 CALL copy(out, lout)
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)
1576 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1577 this%inter_index_y(this%innx,this%inny))
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)
1588 WHERE(c_e(this%inter_index_x(:,:)))
1589 this%inter_index_x(:,:) = this%inter_index_x(:,:) + &
1590 (this%inter_index_y(:,:)-1)*nx
1598 ELSE IF (this%trans%trans_type == 'polyinter') THEN
1604 CALL get_val(in, nx=this%innx, ny=this%inny)
1605 this%outnx = this%innx
1606 this%outny = this%inny
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
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))
1624 DO n = nprev, this%trans%poly%arraysize
1625 IF (inside(point, this%trans%poly%array(n))) THEN
1626 this%inter_index_x(i,j) = n
1631 DO n = nprev-1, 1, -1
1632 IF (inside(point, this%trans%poly%array(n))) THEN
1633 this%inter_index_x(i,j) = n
1646 ELSE IF (this%trans%trans_type == 'metamorphosis') THEN
1649 CALL get_val(in, nx=this%innx, ny=this%inny)
1650 this%outnx = this%innx
1651 this%outny = this%inny
1653 IF (this%trans%sub_type == 'maskvalid' .OR. this%trans%sub_type == 'maskinvalid'THEN
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')
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))
1673 ALLOCATE(this%point_mask(this%innx,this%inny))
1675 IF (this%trans%sub_type == 'maskvalid') THEN
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(:,:))
1683 this%point_mask(:,:) = c_e(maskgrid(:,:)) .AND. &
1684 maskgrid(:,:) > maskbounds(1) .AND. &
1685 maskgrid(:,:) <= maskbounds( SIZE(maskbounds))
1688 this%point_mask(:,:) = .NOT.c_e(maskgrid(:,:))
1693 ELSE IF (this%trans%sub_type == 'setinvalidto') THEN
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')
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')
1706 this%val1 = maskbounds(1)
1708 CALL l4f_category_log(this%category, l4f_debug, &
1709 "grid_transform_init setting invalid data to "//t2c(this%val1))
1715 ELSE IF (this%trans%sub_type == 'settoinvalid') THEN
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')
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')
1730 this%val1 = maskbounds(1)
1731 this%val2 = maskbounds( SIZE(maskbounds))
1733 CALL l4f_category_log(this%category, l4f_debug, &
1734 "grid_transform_init setting to invalid interval ]"//t2c(this%val1 ','
1735 t2c(this%val2)// ']')
1749 SUBROUTINE outgrid_setup()
1756 IF (proj_in == proj_out) THEN
1759 CALL set_val(out, component_flag=cf_i)
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 ')
1766 CALL l4f_category_log(this%category,L4F_WARN, &
1767 'vector fields will probably be wrong ')
1769 CALL set_val(out, component_flag=cf_i)
1772 ! rotate the input grid by n*360 degrees to bring it closer to the output grid
1775 END SUBROUTINE outgrid_setup
1777 END SUBROUTINE grid_transform_init
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:
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
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
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
1833 INTEGER :: ix, iy, n, nm, nr, nprev, nmaskarea, xnmin, xnmax, ynmin, ynmax
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
1842 IF ( PRESENT(find_index)) THEN
1843 IF ( ASSOCIATED(find_index)) THEN
1844 this%find_index => find_index
1854 IF (.NOT. c_e(time_definition)) THEN
1858 IF (this%trans%trans_type == 'inter') THEN
1860 IF (this%trans%sub_type == 'near' .OR. this%trans%sub_type == 'bilin'
1861 .OR. this%trans%sub_type == 'shapiro_near') THEN
1864 CALL get_val(lin, nx=this%innx, ny=this%inny)
1865 this%outnx = SIZE(v7d_out%ana)
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))
1872 CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
1874 lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,
1875 CALL griddim_set_central_lon(lin, lonref)
1876 CALL get_val(lin, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
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)
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
1887 CALL griddim_gen_coord(lin, this%inter_x, this%inter_y)
1888 CALL proj(lin, lon, lat, this%inter_xp, this%inter_yp)
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)
1905 ELSE IF (this%trans%trans_type == 'polyinter') THEN
1907 CALL get_val(in, nx=this%innx, ny=this%inny)
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
1918 this%outnx = this%trans%poly%arraysize
1920 CALL delete(v7d_out)
1921 CALL init(v7d_out, time_definition=time_definition)
1922 CALL vol7d_alloc(v7d_out, nana=this%outnx)
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)
1928 CALL griddim_set_central_lon(lin, lonref)
1933 CALL poly_to_coordinates(this%trans%poly, v7d_out)
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
1942 DO n = nprev, this%trans%poly%arraysize
1943 IF (inside(point, this%trans%poly%array(n))) THEN
1944 this%inter_index_x(ix,iy) = n
1949 DO n = nprev-1, 1, -1
1950 IF (inside(point, this%trans%poly%array(n))) THEN
1951 this%inter_index_x(ix,iy) = n
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)))
|