|
◆ 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 1419 del file grid_transform_class.F90.
1422 CALL proj(in, lout%dim%lon, lout%dim%lat, this%inter_xp, this%inter_yp)
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)
1430 IF (this%trans%trans_type == 'intersearch') THEN
1431 ALLOCATE(this%inter_x(this%innx,this%inny), &
1432 this%inter_y(this%innx,this%inny))
1433 ALLOCATE(this%inter_xp(this%outnx,this%outny), &
1434 this%inter_yp(this%outnx,this%outny))
1437 CALL griddim_gen_coord(in, this%inter_x, this%inter_y)
1439 CALL proj(in, lout%dim%lon, lout%dim%lat, this%inter_xp, this%inter_yp)
1447 ELSE IF (this%trans%trans_type == 'boxinter') THEN
1449 CALL outgrid_setup()
1450 CALL get_val(in, nx=this%innx, ny=this%inny)
1451 CALL get_val(out, nx=this%outnx, ny=this%outny, &
1452 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1455 IF (.NOT.c_e(this%trans%area_info%boxdx)) &
1456 CALL get_val(out, dx=this%trans%area_info%boxdx)
1457 IF (.NOT.c_e(this%trans%area_info%boxdy)) &
1458 CALL get_val(out, dx=this%trans%area_info%boxdy)
1460 this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
1461 this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
1463 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1464 this%inter_index_y(this%innx,this%inny))
1470 CALL this%find_index(out, .true., &
1471 this%outnx, this%outny, xmin, xmax, ymin, ymax, &
1472 lin%dim%lon, lin%dim%lat, .false., &
1473 this%inter_index_x, this%inter_index_y)
1478 ELSE IF (this%trans%trans_type == 'stencilinter') THEN
1480 CALL outgrid_setup()
1482 CALL get_val(in, nx=this%innx, ny=this%inny, &
1483 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1484 CALL get_val(out, nx=this%outnx, ny=this%outny)
1486 ALLOCATE (this%inter_index_x(this%outnx,this%outny), &
1487 this%inter_index_y(this%outnx,this%outny))
1488 CALL copy(out, lout)
1490 CALL this%find_index(in, .true., &
1491 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1492 lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1493 this%inter_index_x, this%inter_index_y)
1496 nr = int(this%trans%area_info%radius)
1499 r2 = this%trans%area_info%radius**2
1500 ALLOCATE(this%stencil(n,n))
1501 this%stencil(:,:) = .true.
1504 IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
1511 xnmax = this%innx - nr
1513 ynmax = this%inny - nr
1514 DO iy = 1, this%outny
1515 DO ix = 1, this%outnx
1516 IF (this%inter_index_x(ix,iy) < xnmin .OR. &
1517 this%inter_index_x(ix,iy) > xnmax .OR. &
1518 this%inter_index_y(ix,iy) < ynmin .OR. &
1519 this%inter_index_y(ix,iy) > ynmax) THEN
1520 this%inter_index_x(ix,iy) = imiss
1521 this%inter_index_y(ix,iy) = imiss
1527 CALL l4f_category_log(this%category, l4f_debug, &
1528 'stencilinter: stencil size '//t2c(n*n))
1529 CALL l4f_category_log(this%category, l4f_debug, &
1530 'stencilinter: stencil points '//t2c(count(this%stencil)))
1536 ELSE IF (this%trans%trans_type == 'maskgen') THEN
1538 IF (this%trans%sub_type == 'poly') THEN
1541 CALL get_val(in, nx=this%innx, ny=this%inny)
1542 this%outnx = this%innx
1543 this%outny = this%inny
1546 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1547 this%inter_index_y(this%innx,this%inny))
1548 this%inter_index_x(:,:) = imiss
1549 this%inter_index_y(:,:) = 1
1558 inside_x: DO i = 1, this%innx
1559 point = georef_coord_new(x=out%dim%lon(i,j), y=out%dim%lat(i,j))
1561 DO n = nprev, this%trans%poly%arraysize
1562 IF (inside(point, this%trans%poly%array(n))) THEN
1563 this%inter_index_x(i,j) = n
1568 DO n = nprev-1, 1, -1
1569 IF (inside(point, this%trans%poly%array(n))) THEN
1570 this%inter_index_x(i,j) = n
1581 ELSE IF (this%trans%sub_type == 'grid') THEN
1584 CALL copy(out, lout)
1587 CALL get_val(in, nx=this%innx, ny=this%inny)
1588 this%outnx = this%innx
1589 this%outny = this%inny
1590 CALL get_val(lout, nx=nx, ny=ny, &
1591 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1594 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1595 this%inter_index_y(this%innx,this%inny))
1601 CALL this%find_index(lout, .true., &
1602 nx, ny, xmin, xmax, ymin, ymax, &
1603 out%dim%lon, out%dim%lat, .false., &
1604 this%inter_index_x, this%inter_index_y)
1606 WHERE(c_e(this%inter_index_x(:,:)))
1607 this%inter_index_x(:,:) = this%inter_index_x(:,:) + &
1608 (this%inter_index_y(:,:)-1)*nx
1616 ELSE IF (this%trans%trans_type == 'polyinter') THEN
1622 CALL get_val(in, nx=this%innx, ny=this%inny)
1623 this%outnx = this%innx
1624 this%outny = this%inny
1627 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1628 this%inter_index_y(this%innx,this%inny))
1629 this%inter_index_x(:,:) = imiss
1630 this%inter_index_y(:,:) = 1
1639 inside_x_2: DO i = 1, this%innx
1640 point = georef_coord_new(x=out%dim%lon(i,j), y=out%dim%lat(i,j))
1642 DO n = nprev, this%trans%poly%arraysize
1643 IF (inside(point, this%trans%poly%array(n))) THEN
1644 this%inter_index_x(i,j) = n
1649 DO n = nprev-1, 1, -1
1650 IF (inside(point, this%trans%poly%array(n))) THEN
1651 this%inter_index_x(i,j) = n
1664 ELSE IF (this%trans%trans_type == 'metamorphosis') THEN
1667 CALL get_val(in, nx=this%innx, ny=this%inny)
1668 this%outnx = this%innx
1669 this%outny = this%inny
1671 IF (this%trans%sub_type == 'maskvalid' .OR. this%trans%sub_type == 'maskinvalid') THEN
1673 IF (.NOT. PRESENT(maskgrid)) THEN
1674 CALL l4f_category_log(this%category,l4f_error, &
1675 'grid_transform_init maskgrid argument missing for metamorphosis:'// &
1676 trim(this%trans%sub_type)// ' transformation')
1681 IF (this%innx /= SIZE(maskgrid,1) .OR. this%inny /= SIZE(maskgrid,2)) THEN
1682 CALL l4f_category_log(this%category,l4f_error, &
1683 'grid_transform_init mask non conformal with input field')
1684 CALL l4f_category_log(this%category,l4f_error, &
1685 'mask: '//t2c( SIZE(maskgrid,1))// 'x'//t2c( SIZE(maskgrid,2))// &
1686 ' input field:'//t2c(this%innx)// 'x'//t2c(this%inny))
1691 ALLOCATE(this%point_mask(this%innx,this%inny))
1693 IF (this%trans%sub_type == 'maskvalid') THEN
1696 IF (.NOT. PRESENT(maskbounds)) THEN
1697 this%point_mask(:,:) = c_e(maskgrid(:,:))
1698 ELSE IF ( SIZE(maskbounds) < 2) THEN
1699 this%point_mask(:,:) = c_e(maskgrid(:,:))
1701 this%point_mask(:,:) = c_e(maskgrid(:,:)) .AND. &
1702 maskgrid(:,:) > maskbounds(1) .AND. &
1703 maskgrid(:,:) <= maskbounds( SIZE(maskbounds))
1706 this%point_mask(:,:) = .NOT.c_e(maskgrid(:,:))
1711 ELSE IF (this%trans%sub_type == 'lemaskinvalid' .OR. &
1712 this%trans%sub_type == 'ltmaskinvalid' .OR. &
1713 this%trans%sub_type == 'gemaskinvalid' .OR. &
1714 this%trans%sub_type == 'gtmaskinvalid') THEN
1717 this%val_mask = maskgrid
1720 ELSE IF (this%trans%sub_type == 'setinvalidto') THEN
1722 IF (.NOT. PRESENT(maskbounds)) THEN
1723 CALL l4f_category_log(this%category,l4f_error, &
1724 'grid_transform_init maskbounds missing for metamorphosis:'// &
1725 trim(this%trans%sub_type)// ' transformation')
1727 ELSE IF ( SIZE(maskbounds) < 1) THEN
1728 CALL l4f_category_log(this%category,l4f_error, &
1729 'grid_transform_init maskbounds empty for metamorphosis:'// &
1730 trim(this%trans%sub_type)// ' transformation')
1733 this%val1 = maskbounds(1)
1735 CALL l4f_category_log(this%category, l4f_debug, &
1736 "grid_transform_init setting invalid data to "//t2c(this%val1))
1742 ELSE IF (this%trans%sub_type == 'settoinvalid') THEN
1744 IF (.NOT. PRESENT(maskbounds)) THEN
1745 CALL l4f_category_log(this%category,l4f_error, &
1746 'grid_transform_init maskbounds missing for metamorphosis:'// &
1747 trim(this%trans%sub_type)// ' transformation')
1750 ELSE IF ( SIZE(maskbounds) < 2) THEN
1751 CALL l4f_category_log(this%category,l4f_error, &
1752 'grid_transform_init maskbounds must have at least 2 elements for metamorphosis:'// &
1753 trim(this%trans%sub_type)// ' transformation')
1757 this%val1 = maskbounds(1)
1758 this%val2 = maskbounds( SIZE(maskbounds))
1760 CALL l4f_category_log(this%category, l4f_debug, &
1761 "grid_transform_init setting to invalid interval ]"//t2c(this%val1)// ','// &
1762 t2c(this%val2)// ']')
1776 SUBROUTINE outgrid_setup()
1779 CALL griddim_setsteps(out)
1781 CALL get_val(in, proj=proj_in, component_flag=cf_i)
1782 CALL get_val(out, proj=proj_out, component_flag=cf_o)
1783 IF (proj_in == proj_out) THEN
1786 CALL set_val(out, component_flag=cf_i)
1791 CALL l4f_category_log(this%category,l4f_warn, &
1792 'trying to interpolate a grid with component flag 1 to a grid on a different projection')
1793 CALL l4f_category_log(this%category,l4f_warn, &
1794 'vector fields will probably be wrong')
1796 CALL set_val(out, component_flag=cf_i)
1800 CALL griddim_set_central_lon(in, griddim_central_lon(out))
1802 END SUBROUTINE outgrid_setup
1804 END SUBROUTINE grid_transform_init
1849 SUBROUTINE grid_transform_grid_vol7d_init(this, trans, in, v7d_out, &
1850 maskgrid, maskbounds, find_index, categoryappend)
1851 TYPE(grid_transform), INTENT(out) :: this
1852 TYPE(transform_def), INTENT(in) :: trans
1853 TYPE(griddim_def), INTENT(in) :: in
1854 TYPE(vol7d), INTENT(inout) :: v7d_out
1855 REAL, INTENT(in), OPTIONAL :: maskgrid(:,:)
1856 REAL, INTENT(in), OPTIONAL :: maskbounds(:)
1857 PROCEDURE(basic_find_index), POINTER, OPTIONAL :: find_index
1858 CHARACTER(len=*), INTENT(in), OPTIONAL :: categoryappend
1860 INTEGER :: ix, iy, n, nm, nr, nprev, nmaskarea, xnmin, xnmax, ynmin, ynmax, &
1862 DOUBLE PRECISION :: xmin, xmax, ymin, ymax, r2, lonref
1863 DOUBLE PRECISION, ALLOCATABLE :: lon1(:), lat1(:), lon(:,:), lat(:,:)
1864 REAL, ALLOCATABLE :: lmaskbounds(:)
1865 TYPE(georef_coord) :: point
1866 TYPE(griddim_def) :: lin
1869 IF ( PRESENT(find_index)) THEN
1870 IF ( ASSOCIATED(find_index)) THEN
1871 this%find_index => find_index
1874 CALL grid_transform_init_common(this, trans, categoryappend)
1876 CALL l4f_category_log(this%category, l4f_debug, "grid_transform vg6d-v7d")
1880 CALL get_val(trans, time_definition=time_definition)
1881 IF (.NOT. c_e(time_definition)) THEN
1885 IF (this%trans%trans_type == 'inter') THEN
1887 IF (this%trans%sub_type == 'near' .OR. this%trans%sub_type == 'bilin' &
1888 .OR. this%trans%sub_type == 'shapiro_near') THEN
1891 CALL get_val(lin, nx=this%innx, ny=this%inny)
1892 this%outnx = SIZE(v7d_out%ana)
1895 ALLOCATE (this%inter_index_x(this%outnx,this%outny),&
1896 this%inter_index_y(this%outnx,this%outny))
1897 ALLOCATE(lon(this%outnx,1),lat(this%outnx,1))
1899 CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
1901 lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1), mask=c_e(lon(:,1))))
1902 CALL griddim_set_central_lon(lin, lonref)
1903 CALL get_val(lin, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1905 IF (this%trans%sub_type == 'bilin') THEN
1906 CALL this%find_index(lin, .false., &
1907 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1908 lon, lat, this%trans%extrap, &
1909 this%inter_index_x, this%inter_index_y)
1911 ALLOCATE(this%inter_x(this%innx,this%inny),this%inter_y(this%innx,this%inny))
1912 ALLOCATE(this%inter_xp(this%outnx,this%outny),this%inter_yp(this%outnx,this%outny))
1914 CALL griddim_gen_coord(lin, this%inter_x, this%inter_y)
1915 CALL proj(lin, lon, lat, this%inter_xp, this%inter_yp)
1918 CALL this%find_index(lin, .true., &
1919 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1920 lon, lat, this%trans%extrap, &
1921 this%inter_index_x, this%inter_index_y)
1923 IF (this%trans%trans_type == 'intersearch') THEN
1924 ALLOCATE(this%inter_x(this%innx,this%inny),this%inter_y(this%innx,this%inny))
1925 ALLOCATE(this%inter_xp(this%outnx,this%outny),this%inter_yp(this%outnx,this%outny))
1927 CALL griddim_gen_coord(lin, this%inter_x, this%inter_y)
1928 CALL proj(lin, lon, lat, this%inter_xp, this%inter_yp)
1939 ELSE IF (this%trans%trans_type == 'polyinter') THEN
1941 CALL get_val(in, nx=this%innx, ny=this%inny)
1943 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1944 this%inter_index_y(this%innx,this%inny))
1945 this%inter_index_x(:,:) = imiss
1946 this%inter_index_y(:,:) = 1
1952 this%outnx = this%trans%poly%arraysize
1954 CALL delete(v7d_out)
1955 CALL init(v7d_out, time_definition=time_definition)
1956 CALL vol7d_alloc(v7d_out, nana=this%outnx)
1959 ALLOCATE(lon(this%outnx,1))
1960 CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1))
1961 lonref = 0.5d0*(maxval(lon(:,1), mask=c_e(lon(:,1))) + minval(lon(:,1), mask=c_e(lon(:,1))))
1962 CALL griddim_set_central_lon(lin, lonref)
1967 CALL poly_to_coordinates(this%trans%poly, v7d_out)
1972 DO iy = 1, this%inny
1973 inside_x: DO ix = 1, this%innx
1974 point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy))
1976 DO n = nprev, this%trans%poly%arraysize
1977 IF (inside(point, this%trans%poly%array(n))) THEN
1978 this%inter_index_x(ix,iy) = n
1983 DO n = nprev-1, 1, -1
1984 IF (inside(point, this%trans%poly%array(n))) THEN
1985 this%inter_index_x(ix,iy) = n
1997 DO n = 1, this%trans%poly%arraysize
1998 CALL l4f_category_log(this%category, l4f_debug, &
1999 'Polygon: '//t2c(n)// ' grid points: '// &
2000 t2c(count(this%inter_index_x(:,:) == n)))
|