|
◆ 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 1429 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)))
|