|
◆ 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] | 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 1427 del file grid_transform_class.F90.
1430 CALL proj(in, lout%dim%lon, lout%dim%lat, this%inter_xp, this%inter_yp)
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)
1438 IF (this%trans%trans_type == 'intersearch') THEN
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))
1445 CALL griddim_gen_coord(in, this%inter_x, this%inter_y)
1447 CALL proj(in, lout%dim%lon, lout%dim%lat, this%inter_xp, this%inter_yp)
1455 ELSE IF (this%trans%trans_type == 'boxinter') THEN
1457 CALL outgrid_setup()
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)
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)
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
1471 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1472 this%inter_index_y(this%innx,this%inny))
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)
1486 ELSE IF (this%trans%trans_type == 'stencilinter') THEN
1488 CALL outgrid_setup()
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)
1494 ALLOCATE (this%inter_index_x(this%outnx,this%outny), &
1495 this%inter_index_y(this%outnx,this%outny))
1496 CALL copy(out, 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)
1504 nr = int(this%trans%area_info%radius)
1507 r2 = this%trans%area_info%radius**2
1508 ALLOCATE(this%stencil(n,n))
1509 this%stencil(:,:) = .true.
1512 IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
1519 xnmax = this%innx - nr
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
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)))
1544 ELSE IF (this%trans%trans_type == 'maskgen') THEN
1546 IF (this%trans%sub_type == 'poly') THEN
1549 CALL get_val(in, nx=this%innx, ny=this%inny)
1550 this%outnx = this%innx
1551 this%outny = this%inny
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
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))
1569 DO n = nprev, this%trans%poly%arraysize
1570 IF (inside(point, this%trans%poly%array(n))) THEN
1571 this%inter_index_x(i,j) = n
1576 DO n = nprev-1, 1, -1
1577 IF (inside(point, this%trans%poly%array(n))) THEN
1578 this%inter_index_x(i,j) = n
1589 ELSE IF (this%trans%sub_type == 'grid') THEN
1592 CALL copy(out, lout)
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)
1602 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1603 this%inter_index_y(this%innx,this%inny))
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)
1614 WHERE(c_e(this%inter_index_x(:,:)))
1615 this%inter_index_x(:,:) = this%inter_index_x(:,:) + &
1616 (this%inter_index_y(:,:)-1)*nx
1624 ELSE IF (this%trans%trans_type == 'polyinter') THEN
1630 CALL get_val(in, nx=this%innx, ny=this%inny)
1631 this%outnx = this%innx
1632 this%outny = this%inny
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
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))
1650 DO n = nprev, this%trans%poly%arraysize
1651 IF (inside(point, this%trans%poly%array(n))) THEN
1652 this%inter_index_x(i,j) = n
1657 DO n = nprev-1, 1, -1
1658 IF (inside(point, this%trans%poly%array(n))) THEN
1659 this%inter_index_x(i,j) = n
1672 ELSE IF (this%trans%trans_type == 'metamorphosis') THEN
1675 CALL get_val(in, nx=this%innx, ny=this%inny)
1676 this%outnx = this%innx
1677 this%outny = this%inny
1679 IF (this%trans%sub_type == 'maskvalid' .OR. this%trans%sub_type == 'maskinvalid') THEN
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')
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))
1699 ALLOCATE(this%point_mask(this%innx,this%inny))
1701 IF (this%trans%sub_type == 'maskvalid') THEN
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(:,:))
1709 this%point_mask(:,:) = c_e(maskgrid(:,:)) .AND. &
1710 maskgrid(:,:) > maskbounds(1) .AND. &
1711 maskgrid(:,:) <= maskbounds( SIZE(maskbounds))
1714 this%point_mask(:,:) = .NOT.c_e(maskgrid(:,:))
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
1725 this%val_mask = maskgrid
1728 ELSE IF (this%trans%sub_type == 'setinvalidto') THEN
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')
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')
1741 this%val1 = maskbounds(1)
1743 CALL l4f_category_log(this%category, l4f_debug, &
1744 "grid_transform_init setting invalid data to "//t2c(this%val1))
1750 ELSE IF (this%trans%sub_type == 'settoinvalid') THEN
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')
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')
1765 this%val1 = maskbounds(1)
1766 this%val2 = maskbounds( SIZE(maskbounds))
1768 CALL l4f_category_log(this%category, l4f_debug, &
1769 "grid_transform_init setting to invalid interval ]"//t2c(this%val1)// ','// &
1770 t2c(this%val2)// ']')
1784 SUBROUTINE outgrid_setup()
1787 CALL griddim_setsteps(out)
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
1794 CALL set_val(out, component_flag=cf_i)
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')
1804 CALL set_val(out, component_flag=cf_i)
1808 CALL griddim_set_central_lon(in, griddim_central_lon(out))
1810 END SUBROUTINE outgrid_setup
1812 END SUBROUTINE grid_transform_init
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
1868 INTEGER :: ix, iy, n, nm, nr, nprev, nmaskarea, xnmin, xnmax, ynmin, ynmax, &
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
1877 IF ( PRESENT(find_index)) THEN
1878 IF ( ASSOCIATED(find_index)) THEN
1879 this%find_index => find_index
1882 CALL grid_transform_init_common(this, trans, categoryappend)
1884 CALL l4f_category_log(this%category, l4f_debug, "grid_transform vg6d-v7d")
1888 CALL get_val(trans, time_definition=time_definition)
1889 IF (.NOT. c_e(time_definition)) THEN
1893 IF (this%trans%trans_type == 'inter') THEN
1895 IF (this%trans%sub_type == 'near' .OR. this%trans%sub_type == 'bilin' &
1896 .OR. this%trans%sub_type == 'shapiro_near') THEN
1899 CALL get_val(lin, nx=this%innx, ny=this%inny)
1900 this%outnx = SIZE(v7d_out%ana)
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))
1907 CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
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)
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)
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))
1922 CALL griddim_gen_coord(lin, this%inter_x, this%inter_y)
1923 CALL proj(lin, lon, lat, this%inter_xp, this%inter_yp)
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)
1931 IF (this%trans%trans_type == 'intersearch') THEN
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))
1935 CALL griddim_gen_coord(lin, this%inter_x, this%inter_y)
1936 CALL proj(lin, lon, lat, this%inter_xp, this%inter_yp)
1947 ELSE IF (this%trans%trans_type == 'polyinter') THEN
1949 CALL get_val(in, nx=this%innx, ny=this%inny)
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
1960 this%outnx = this%trans%poly%arraysize
1962 CALL delete(v7d_out)
1963 CALL init(v7d_out, time_definition=time_definition)
1964 CALL vol7d_alloc(v7d_out, nana=this%outnx)
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)
1975 CALL poly_to_coordinates(this%trans%poly, v7d_out)
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))
1984 DO n = nprev, this%trans%poly%arraysize
1985 IF (inside(point, this%trans%poly%array(n))) THEN
1986 this%inter_index_x(ix,iy) = n
1991 DO n = nprev-1, 1, -1
1992 IF (inside(point, this%trans%poly%array(n))) THEN
1993 this%inter_index_x(ix,iy) = n
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)))
|