1417 CALL this%find_index(in, .false., &
1418 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1419 lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1420 this%inter_index_x, this%inter_index_y)
1422 ALLOCATE(this%inter_x(this%innx,this%inny), &
1423 this%inter_y(this%innx,this%inny))
1424 ALLOCATE(this%inter_xp(this%outnx,this%outny), &
1425 this%inter_yp(this%outnx,this%outny))
1428 CALL griddim_gen_coord(in, this%inter_x, this%inter_y)
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)
1455ELSE 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)
1486ELSE 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)))
1544ELSE 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
1624ELSE 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
1672ELSE 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)//
']')
1784SUBROUTINE outgrid_setup()
1787CALL griddim_setsteps(out)
1789CALL get_val(in,
proj=proj_in, component_flag=cf_i)
1790CALL get_val(out,
proj=proj_out, component_flag=cf_o)
1791IF (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)
1808CALL griddim_set_central_lon(in, griddim_central_lon(out))
1810END SUBROUTINE outgrid_setup
1812END SUBROUTINE grid_transform_init
1857SUBROUTINE grid_transform_grid_vol7d_init(this, trans, in, v7d_out, &
1858 maskgrid, maskbounds, find_index, categoryappend)
1862TYPE(
vol7d),
INTENT(inout) :: v7d_out
1863REAL,
INTENT(in),
OPTIONAL :: maskgrid(:,:)
1864REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
1865PROCEDURE(basic_find_index),
POINTER,
OPTIONAL :: find_index
1866CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
1868INTEGER :: ix, iy, n, nm, nr, nprev, nmaskarea, xnmin, xnmax, ynmin, ynmax, &
1870DOUBLE PRECISION :: xmin, xmax, ymin, ymax, r2, lonref
1871DOUBLE PRECISION,
ALLOCATABLE :: lon1(:), lat1(:), lon(:,:), lat(:,:)
1872REAL,
ALLOCATABLE :: lmaskbounds(:)
1877IF (
PRESENT(find_index))
THEN
1878 IF (
ASSOCIATED(find_index))
THEN
1879 this%find_index => find_index
1882CALL grid_transform_init_common(this, trans, categoryappend)
1884CALL l4f_category_log(this%category, l4f_debug,
"grid_transform vg6d-v7d")
1888CALL get_val(trans, time_definition=time_definition)
1889IF (.NOT.
c_e(time_definition))
THEN
1893IF (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)
1947ELSE 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
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
2045 this%stencil(:,:) = .true.
2048 IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
2055 xnmax = this%innx - nr
2057 ynmax = this%inny - nr
2058 DO iy = 1, this%outny
2059 DO ix = 1, this%outnx
2060 IF (this%inter_index_x(ix,iy) < xnmin .OR. &
2061 this%inter_index_x(ix,iy) > xnmax .OR. &
2062 this%inter_index_y(ix,iy) < ynmin .OR. &
2063 this%inter_index_y(ix,iy) > ynmax)
THEN
2064 this%inter_index_x(ix,iy) = imiss
2065 this%inter_index_y(ix,iy) = imiss
2071 CALL l4f_category_log(this%category, l4f_debug, &
2072 'stencilinter: stencil size '//t2c(n*n))
2073 CALL l4f_category_log(this%category, l4f_debug, &
2074 'stencilinter: stencil points '//t2c(count(this%stencil)))
2082ELSE IF (this%trans%trans_type ==
'maskinter')
THEN
2084 IF (.NOT.
PRESENT(maskgrid))
THEN
2085 CALL l4f_category_log(this%category,l4f_error, &
2086 'grid_transform_init maskgrid argument missing for maskinter transformation')
2087 CALL raise_fatal_error()
2090 CALL get_val(in, nx=this%innx, ny=this%inny)
2091 IF (this%innx /=
SIZE(maskgrid,1) .OR. this%inny /=
SIZE(maskgrid,2))
THEN
2092 CALL l4f_category_log(this%category,l4f_error, &
2093 'grid_transform_init mask non conformal with input field')
2094 CALL l4f_category_log(this%category,l4f_error, &
2095 'mask: '//t2c(
SIZE(maskgrid,1))//
'x'//t2c(
SIZE(maskgrid,2))// &
2096 ' input field:'//t2c(this%innx)//
'x'//t2c(this%inny))
2097 CALL raise_fatal_error()
2100 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2101 this%inter_index_y(this%innx,this%inny))
2102 this%inter_index_x(:,:) = imiss
2103 this%inter_index_y(:,:) = 1
2106 CALL gen_mask_class()
2114 DO iy = 1, this%inny
2115 DO ix = 1, this%innx
2116 IF (
c_e(maskgrid(ix,iy)))
THEN
2117 IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1))
THEN
2118 DO n = nmaskarea, 1, -1
2119 IF (maskgrid(ix,iy) > lmaskbounds(n))
THEN
2120 this%inter_index_x(ix,iy) = n
2130 this%outnx = nmaskarea
2133 CALL init(v7d_out, time_definition=time_definition)
2134 CALL vol7d_alloc(v7d_out, nana=nmaskarea)
2139 CALL init(v7d_out%ana(n), &
2141 mask=(this%inter_index_x(:,:) == n))), &
2143 mask=(this%inter_index_x(:,:) == n))))
2149ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN
2156 CALL get_val(in, nx=this%innx, ny=this%inny)
2158 ALLOCATE(this%point_index(this%innx,this%inny))
2159 this%point_index(:,:) = imiss
2162 CALL init(v7d_out, time_definition=time_definition)
2164 IF (this%trans%sub_type ==
'all' )
THEN
2166 this%outnx = this%innx*this%inny
2168 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2173 CALL init(v7d_out%ana((iy-1)*this%innx+ix), &
2174 lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
2176 this%point_index(ix,iy) = n
2182 ELSE IF (this%trans%sub_type ==
'coordbb' )
THEN
2187 DO iy = 1, this%inny
2188 DO ix = 1, this%innx
2190 IF (lin%dim%lon(ix,iy) > this%trans%rect_coo%ilon .AND. &
2191 lin%dim%lon(ix,iy) < this%trans%rect_coo%flon .AND. &
2192 lin%dim%lat(ix,iy) > this%trans%rect_coo%ilat .AND. &
2193 lin%dim%lat(ix,iy) < this%trans%rect_coo%flat)
THEN
2194 this%outnx = this%outnx + 1
2195 this%point_index(ix,iy) = this%outnx
2200 IF (this%outnx <= 0)
THEN
2201 CALL l4f_category_log(this%category,l4f_warn, &
2202 "metamorphosis:coordbb: no points inside bounding box "//&
2203 trim(
to_char(this%trans%rect_coo%ilon))//
","// &
2204 trim(
to_char(this%trans%rect_coo%flon))//
","// &
2205 trim(
to_char(this%trans%rect_coo%ilat))//
","// &
2206 trim(
to_char(this%trans%rect_coo%flat)))
2209 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2213 DO iy = 1, this%inny
2214 DO ix = 1, this%innx
2215 IF (
c_e(this%point_index(ix,iy)))
THEN
2217 CALL init(v7d_out%ana(n), &
2218 lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
2225 ELSE IF (this%trans%sub_type ==
'poly' )
THEN
2238 DO iy = 1, this%inny
2239 DO ix = 1, this%innx
2240 point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy))
2241 DO n = 1, this%trans%poly%arraysize
2242 IF (
inside(point, this%trans%poly%array(n)))
THEN
2246 this%outnx = this%outnx + 1
2248 this%point_index(ix,iy) = n
2257 IF (this%outnx <= 0)
THEN
2258 CALL l4f_category_log(this%category,l4f_warn, &
2259 "metamorphosis:poly: no points inside polygons")
2262 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2265 DO iy = 1, this%inny
2266 DO ix = 1, this%innx
2267 IF (
c_e(this%point_index(ix,iy)))
THEN
2269 CALL init(v7d_out%ana(n), &
2270 lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
2277 ELSE IF (this%trans%sub_type ==
'mask' )
THEN
2279 IF (.NOT.
PRESENT(maskgrid))
THEN
2280 CALL l4f_category_log(this%category,l4f_error, &
2281 'grid_transform_init maskgrid argument missing for metamorphosis:mask transformation')
2286 IF (this%innx /=
SIZE(maskgrid,1) .OR. this%inny /=
SIZE(maskgrid,2))
THEN
2287 CALL l4f_category_log(this%category,l4f_error, &
2288 'grid_transform_init mask non conformal with input field')
2289 CALL l4f_category_log(this%category,l4f_error, &
2290 'mask: '//t2c(
SIZE(maskgrid,1))//
'x'//t2c(
SIZE(maskgrid,2))// &
2291 ' input field:'//t2c(this%innx)//
'x'//t2c(this%inny))
2297 CALL gen_mask_class()
2309 DO iy = 1, this%inny
2310 DO ix = 1, this%innx
2311 IF (
c_e(maskgrid(ix,iy)))
THEN
2312 IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1))
THEN
2313 DO n = nmaskarea, 1, -1
2314 IF (maskgrid(ix,iy) > lmaskbounds(n))
THEN
2318 this%outnx = this%outnx + 1
2320 this%point_index(ix,iy) = n
2331 IF (this%outnx <= 0)
THEN
2332 CALL l4f_category_log(this%category,l4f_warn, &
2333 "grid_transform_init no points inside mask for metamorphosis:mask transformation")
2337 CALL l4f_category_log(this%category,l4f_info, &
2338 "points in subarea "//t2c(n)//
": "// &
2339 t2c(count(this%point_index(:,:) == n)))
2342 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2345 DO iy = 1, this%inny
2346 DO ix = 1, this%innx
2347 IF (
c_e(this%point_index(ix,iy)))
THEN
2349 CALL init(v7d_out%ana(n),lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
2362SUBROUTINE gen_mask_class()
2363REAL :: startmaskclass, mmin, mmax
2366IF (
PRESENT(maskbounds))
THEN
2367 nmaskarea =
SIZE(maskbounds) - 1
2368 IF (nmaskarea > 0)
THEN
2369 lmaskbounds = maskbounds
2371 ELSE IF (nmaskarea == 0)
THEN
2372 CALL l4f_category_log(this%category,l4f_warn, &
2373 'only one value provided for maskbounds, '//t2c(maskbounds(1)) &
2374 //
', it will be ignored')
2375 CALL l4f_category_log(this%category,l4f_warn, &
2376 'at least 2 values are required for maskbounds')
2380mmin = minval(maskgrid, mask=
c_e(maskgrid))
2381mmax = maxval(maskgrid, mask=
c_e(maskgrid))
2383nmaskarea = int(mmax - mmin + 1.5)
2384startmaskclass = mmin - 0.5
2386ALLOCATE(lmaskbounds(nmaskarea+1))
2387lmaskbounds(:) = (/(startmaskclass+real(i),i=0,nmaskarea)/)
2389CALL l4f_category_log(this%category,l4f_debug, &
2390 'maskinter, '//t2c(nmaskarea)//
' subareas defined, with boundaries:')
2391DO i = 1,
SIZE(lmaskbounds)
2392 CALL l4f_category_log(this%category,l4f_debug, &
2393 'maskinter '//t2c(i)//
' '//t2c(lmaskbounds(i)))
2397END SUBROUTINE gen_mask_class
2399END SUBROUTINE grid_transform_grid_vol7d_init
2420SUBROUTINE grid_transform_vol7d_grid_init(this, trans, v7d_in, out, categoryappend)
2423TYPE(
vol7d),
INTENT(in) :: v7d_in
2425character(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2428DOUBLE PRECISION :: xmin, xmax, ymin, ymax, lonref
2429DOUBLE PRECISION,
ALLOCATABLE :: lon(:,:),lat(:,:)
2433CALL grid_transform_init_common(this, trans, categoryappend)
2435CALL l4f_category_log(this%category, l4f_debug,
"grid_transform v7d-vg6d")
2438IF (this%trans%trans_type ==
'inter')
THEN
2440 IF ( this%trans%sub_type ==
'linear' )
THEN
2442 this%innx=
SIZE(v7d_in%ana)
2444 ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2445 ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2446 CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2448 CALL copy (out, lout)
2450 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
2451 CALL griddim_set_central_lon(lout, lonref)
2453 CALL get_val(lout, nx=nx, ny=ny)
2456 ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2458 CALL get_val(lout, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2459 CALL proj(lout, lon, lat, this%inter_xp, this%inter_yp)
2460 CALL griddim_gen_coord(lout, this%inter_x, this%inter_y)
2469ELSE IF (this%trans%trans_type ==
'boxinter')
THEN
2471 this%innx=
SIZE(v7d_in%ana)
2474 ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2475 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2476 this%inter_index_y(this%innx,this%inny))
2478 CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2480 CALL copy (out, lout)
2482 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
2483 CALL griddim_set_central_lon(lout, lonref)
2485 CALL get_val(lout, nx=this%outnx, ny=this%outny, &
2486 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2489 IF (.NOT.
c_e(this%trans%area_info%boxdx)) &
2490 CALL get_val(out, dx=this%trans%area_info%boxdx)
2491 IF (.NOT.
c_e(this%trans%area_info%boxdy)) &
2492 CALL get_val(out, dx=this%trans%area_info%boxdy)
2494 this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
2495 this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
2498 CALL this%find_index(lout, .true., &
2499 this%outnx, this%outny, xmin, xmax, ymin, ymax, &
2500 lon, lat, .false., &
2501 this%inter_index_x, this%inter_index_y)
2510END SUBROUTINE grid_transform_vol7d_grid_init
2547SUBROUTINE grid_transform_vol7d_vol7d_init(this, trans, v7d_in, v7d_out, &
2548 maskbounds, categoryappend)
2551TYPE(
vol7d),
INTENT(in) :: v7d_in
2552TYPE(
vol7d),
INTENT(inout) :: v7d_out
2553REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
2554CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2557DOUBLE PRECISION,
ALLOCATABLE :: lon(:), lat(:)
2559DOUBLE PRECISION :: lon1, lat1
2563CALL grid_transform_init_common(this, trans, categoryappend)
2565CALL l4f_category_log(this%category, l4f_debug,
"grid_transform v7d-v7d")
2568IF (this%trans%trans_type ==
'inter')
THEN
2570 IF ( this%trans%sub_type ==
'linear' )
THEN
2572 CALL vol7d_alloc_vol(v7d_out)
2573 this%outnx=
SIZE(v7d_out%ana)
2576 this%innx=
SIZE(v7d_in%ana)
2579 ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
3125 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3126 this%trans%sub_type ==
'stddevnm1')
THEN
3128 IF (this%trans%sub_type ==
'stddev')
THEN
3134 navg = this%trans%box_info%npx*this%trans%box_info%npy
3137 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3138 je = j+this%trans%box_info%npy-1
3141 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3142 ie = i+this%trans%box_info%npx-1
3145 reshape(field_in(i:ie,j:je,k),(/navg/)), nm1=nm1)
3150 ELSE IF (this%trans%sub_type ==
'max')
THEN
3153 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3154 je = j+this%trans%box_info%npy-1
3157 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3158 ie = i+this%trans%box_info%npx-1
3160 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3162 field_out(ii,jj,k) = maxval(field_in(i:ie,j:je,k), &
3163 mask=(field_in(i:ie,j:je,k) /= rmiss))
3169 ELSE IF (this%trans%sub_type ==
'min')
THEN
3172 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3173 je = j+this%trans%box_info%npy-1
3176 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3177 ie = i+this%trans%box_info%npx-1
3179 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3181 field_out(ii,jj,k) = minval(field_in(i:ie,j:je,k), &
3182 mask=(field_in(i:ie,j:je,k) /= rmiss))
3188 ELSE IF (this%trans%sub_type ==
'percentile')
THEN
3190 navg = this%trans%box_info%npx*this%trans%box_info%npy
3193 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3194 je = j+this%trans%box_info%npy-1
3197 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3198 ie = i+this%trans%box_info%npx-1
3201 reshape(field_in(i:ie,j:je,k),(/navg/)), &
3202 (/real(this%trans%stat_info%percentile)/))
3207 ELSE IF (this%trans%sub_type ==
'frequency')
THEN
3211 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3212 je = j+this%trans%box_info%npy-1
3215 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3216 ie = i+this%trans%box_info%npx-1
3218 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3220 field_out(ii,jj,k) = sum(interval_info_valid( &
3221 this%trans%interval_info, field_in(i:ie,j:je,k)), &
3222 mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
3230ELSE IF (this%trans%trans_type ==
'inter')
THEN
3232 IF (this%trans%sub_type ==
'near')
THEN
3235 DO j = 1, this%outny
3236 DO i = 1, this%outnx
3238 IF (
c_e(this%inter_index_x(i,j))) field_out(i,j,k) = &
3239 field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3245 ELSE IF (this%trans%sub_type ==
'bilin')
THEN
3248 DO j = 1, this%outny
3249 DO i = 1, this%outnx
3251 IF (
c_e(this%inter_index_x(i,j)))
THEN
3253 z1=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3254 z2=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j),k)
3255 z3=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1,k)
3256 z4=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j)+1,k)
3258 IF (
c_e(z1) .AND.
c_e(z2) .AND.
c_e(z3) .AND.
c_e(z4))
THEN
3260 x1=this%inter_x(this%inter_index_x(i,j),this%inter_index_y(i,j))
3261 y1=this%inter_y(this%inter_index_x(i,j),this%inter_index_y(i,j))
3262 x3=this%inter_x(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3263 y3=this%inter_y(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3265 xp=this%inter_xp(i,j)
3266 yp=this%inter_yp(i,j)
3268 field_out(i,j,k) = hbilin(z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
3276 ELSE IF (this%trans%sub_type ==
'shapiro_near')
THEN
3278 DO j = 1, this%outny
3279 DO i = 1, this%outnx
3281 IF (
c_e(this%inter_index_x(i,j)))
THEN
3283 IF(this%inter_index_x(i,j)-1>0)
THEN
3284 z(1) = field_in(this%inter_index_x(i,j)-1,this%inter_index_y(i,j) ,k)
3288 IF(this%inter_index_x(i,j)+1<=this%outnx)
THEN
3289 z(3)=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j) ,k)
3293 IF(this%inter_index_y(i,j)+1<=this%outny)
THEN
3294 z(2)=field_in(this%inter_index_x(i,j) ,this%inter_index_y(i,j)+1,k)
3298 IF(this%inter_index_y(i,j)-1>0)
THEN
3299 z(4)=field_in(this%inter_index_x(i,j), this%inter_index_y(i,j)-1,k)
3303 field_out(i,j,k) = shapiro(z,field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k))
3312ELSE IF (this%trans%trans_type ==
'intersearch')
THEN
3315 likethis%trans%trans_type =
'inter'
3316 CALL grid_transform_compute(likethis, field_in, field_out, var, coord_3d_in)
3317 CALL getenv(
'LIBSIM_DISABLEOPTSEARCH', env_var)
3318 optsearch = len_trim(env_var) == 0
3321 IF ((.NOT.all(
c_e(field_out(:,:,k)))) .AND. (any(
c_e(field_in(:,:,k)))))
THEN
3322 DO j = 1, this%outny
3323 DO i = 1, this%outnx
3324 IF (.NOT.
c_e(field_out(i,j,k)))
THEN
3328 ix = this%inter_index_x(i,j)
3329 iy = this%inter_index_y(i,j)
3330 DO s = 0, max(this%innx, this%inny)
3332 DO m = iy-s, iy+s, max(2*s, 1)
3333 IF (m < 1 .OR. m > this%inny) cycle
3334 DO l = max(1, ix-s), min(this%innx, ix+s)
3335 disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2
3336 IF (
c_e(field_in(l,m,k)))
THEN
3337 IF (disttmp <
dist)
THEN
3339 field_out(i,j,k) = field_in(l,m,k)
3341 ELSE IF (disttmp ==
dist)
THEN
3342 field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k)
3343 nearcount = nearcount + 1
3346 IF (disttmp <
dist) farenough = .false.
3349 DO m = max(1, iy-s+1), min(this%inny, iy+s-1)
3350 DO l = ix-s, ix+s, 2*s
3351 IF (l < 1 .OR. l > this%innx) cycle
3352 disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2
3353 IF (
c_e(field_in(l,m,k)))
THEN
3354 IF (disttmp <
dist)
THEN
3356 field_out(i,j,k) = field_in(l,m,k)
3358 ELSE IF (disttmp ==
dist)
THEN
3359 field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k)
3360 nearcount = nearcount + 1
3363 IF (disttmp <
dist) farenough = .false.
3366 IF (s > 0 .AND. farenough)
EXIT
3371 IF (
c_e(field_in(l,m,k)))
THEN
3372 disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2
3373 IF (disttmp <
dist)
THEN
3375 field_out(i,j,k) = field_in(l,m,k)
3377 ELSE IF (disttmp ==
dist)
THEN
3378 field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k)
3379 nearcount = nearcount + 1
3386 IF (nearcount > 1) field_out(i,j,k) = field_out(i,j,k)/nearcount
3393ELSE IF (this%trans%trans_type ==
'boxinter' &
3394 .OR. this%trans%trans_type ==
'polyinter' &
3395 .OR. this%trans%trans_type ==
'maskinter')
THEN
3397 IF (this%trans%sub_type ==
'average')
THEN
3399 IF (vartype == var_dir360)
THEN
3401 DO j = 1, this%outny
3402 DO i = 1, this%outnx
3403 field_out(i,j,k) = find_prevailing_direction(field_in(:,:,k), &
3405 mask=this%inter_index_x(:,:) == i .AND. this%inter_index_y(:,:) == j)
3411 ALLOCATE(nval(this%outnx, this%outny))
3412 field_out(:,:,:) = 0.0
3417 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3418 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3419 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3421 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3422 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3426 WHERE (nval(:,:) /= 0)
3427 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3429 field_out(:,:,k) = rmiss
3435 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3436 this%trans%sub_type ==
'stddevnm1')
THEN
3438 IF (this%trans%sub_type ==
'stddev')
THEN
3444 DO j = 1, this%outny
3445 DO i = 1, this%outnx
3448 reshape(field_in(:,:,k), (/
SIZE(field_in(:,:,k))/)), &
3449 mask=reshape((this%inter_index_x == i .AND. &
3450 this%inter_index_y == j), (/
SIZE(field_in(:,:,k))/)), nm1=nm1)
3455 ELSE IF (this%trans%sub_type ==
'max')
THEN
3460 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3461 IF (
c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k)))
THEN
3462 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3463 max(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3466 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3475 ELSE IF (this%trans%sub_type ==
'min')
THEN
3480 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3481 IF (
c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k)))
THEN
3482 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3483 min(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3486 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3494 ELSE IF (this%trans%sub_type ==
'percentile')
THEN
3497 DO j = 1, this%outny
3498 DO i = 1, this%outnx
3501 reshape(field_in(:,:,k), (/
SIZE(field_in(:,:,k))/)), &
3502 (/real(this%trans%stat_info%percentile)/), &
3503 mask=reshape((this%inter_index_x == i .AND. &
3504 this%inter_index_y == j), (/
SIZE(field_in(:,:,k))/)))
3509 ELSE IF (this%trans%sub_type ==
'frequency')
THEN
3511 ALLOCATE(nval(this%outnx, this%outny))
3512 field_out(:,:,:) = 0.0
3517 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3518 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3519 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3520 interval_info_valid(this%trans%interval_info, field_in(i,j,k))
3521 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3522 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3526 WHERE (nval(:,:) /= 0)
3527 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3529 field_out(:,:,k) = rmiss
3536ELSE IF (this%trans%trans_type ==
'stencilinter')
THEN
3537 np =
SIZE(this%stencil,1)/2
3538 ns =
SIZE(this%stencil)
3540 IF (this%trans%sub_type ==
'average')
THEN
3542 IF (vartype == var_dir360)
THEN
3544 DO j = 1, this%outny
3545 DO i = 1, this%outnx
3546 IF (
c_e(this%inter_index_x(i,j)))
THEN
3547 i1 = this%inter_index_x(i,j) - np
3548 i2 = this%inter_index_x(i,j) + np
3549 j1 = this%inter_index_y(i,j) - np
3550 j2 = this%inter_index_y(i,j) + np
3551 field_out(i,j,k) = find_prevailing_direction(field_in(i1:i2,j1:j2,k), &
3553 mask=this%stencil(:,:))
3564 DO j = 1, this%outny
3565 DO i = 1, this%outnx
3566 IF (
c_e(this%inter_index_x(i,j)))
THEN
3567 i1 = this%inter_index_x(i,j) - np
3568 i2 = this%inter_index_x(i,j) + np
3569 j1 = this%inter_index_y(i,j) - np
3570 j2 = this%inter_index_y(i,j) + np
3571 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3573 field_out(i,j,k) = sum(field_in(i1:i2,j1:j2,k), &
3574 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3583 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3584 this%trans%sub_type ==
'stddevnm1')
THEN
3586 IF (this%trans%sub_type ==
'stddev')
THEN
3595 DO j = 1, this%outny
3596 DO i = 1, this%outnx
3597 IF (
c_e(this%inter_index_x(i,j)))
THEN
3598 i1 = this%inter_index_x(i,j) - np
3599 i2 = this%inter_index_x(i,j) + np
3600 j1 = this%inter_index_y(i,j) - np
3601 j2 = this%inter_index_y(i,j) + np
3604 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3605 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3606 this%stencil(:,:), (/ns/)), nm1=nm1)
3613 ELSE IF (this%trans%sub_type ==
'max')
THEN
3618 DO j = 1, this%outny
3619 DO i = 1, this%outnx
3620 IF (
c_e(this%inter_index_x(i,j)))
THEN
3621 i1 = this%inter_index_x(i,j) - np
3622 i2 = this%inter_index_x(i,j) + np
3623 j1 = this%inter_index_y(i,j) - np
3624 j2 = this%inter_index_y(i,j) + np
3625 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3627 field_out(i,j,k) = maxval(field_in(i1:i2,j1:j2,k), &
3628 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3636 ELSE IF (this%trans%sub_type ==
'min')
THEN
3641 DO j = 1, this%outny
3642 DO i = 1, this%outnx
3643 IF (
c_e(this%inter_index_x(i,j)))
THEN
3644 i1 = this%inter_index_x(i,j) - np
3645 i2 = this%inter_index_x(i,j) + np
3646 j1 = this%inter_index_y(i,j) - np
3647 j2 = this%inter_index_y(i,j) + np
3648 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3650 field_out(i,j,k) = minval(field_in(i1:i2,j1:j2,k), &
3651 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3659 ELSE IF (this%trans%sub_type ==
'percentile')
THEN
3664 DO j = 1, this%outny
3665 DO i = 1, this%outnx
3666 IF (
c_e(this%inter_index_x(i,j)))
THEN
3667 i1 = this%inter_index_x(i,j) - np
3668 i2 = this%inter_index_x(i,j) + np
3669 j1 = this%inter_index_y(i,j) - np
3670 j2 = this%inter_index_y(i,j) + np
3673 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3674 (/real(this%trans%stat_info%percentile)/), &
3675 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3676 this%stencil(:,:), (/ns/)))
3683 ELSE IF (this%trans%sub_type ==
'frequency')
THEN
3688 DO j = 1, this%outny
3689 DO i = 1, this%outnx
3690 IF (
c_e(this%inter_index_x(i,j)))
THEN
3691 i1 = this%inter_index_x(i,j) - np
3692 i2 = this%inter_index_x(i,j) + np
3693 j1 = this%inter_index_y(i,j) - np
3694 j2 = this%inter_index_y(i,j) + np
3695 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3697 field_out(i,j,k) = sum(interval_info_valid( &
3698 this%trans%interval_info, field_in(i1:i2,j1:j2,k)), &
3699 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3709ELSE IF (this%trans%trans_type ==
'maskgen')
THEN
3712 WHERE(
c_e(this%inter_index_x(:,:)))
3713 field_out(:,:,k) = real(this%inter_index_x(:,:))
3717ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN
3719 IF (this%trans%sub_type ==
'all')
THEN
3721 field_out(:,:,:) = reshape(field_in(:,:,:), (/this%outnx,this%outny,innz/))
3723 ELSE IF (this%trans%sub_type ==
'coordbb' .OR. this%trans%sub_type ==
'poly' &
3724 .OR. this%trans%sub_type ==
'mask')
THEN
3728 field_out(:,1,k) = pack(field_in(:,:,k),
c_e(this%point_index(:,:)))
3731 ELSE IF (this%trans%sub_type ==
'maskvalid' .OR. &
3732 this%trans%sub_type ==
'maskinvalid')
THEN
3735 WHERE (this%point_mask(:,:))
3736 field_out(:,:,k) = field_in(:,:,k)
3740 ELSE IF (this%trans%sub_type ==
'lemaskinvalid')
THEN
3743 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) > this%val_mask(:,:))
3744 field_out(:,:,k) = field_in(:,:,k)
3746 field_out(:,:,k) = rmiss
3750 ELSE IF (this%trans%sub_type ==
'ltmaskinvalid')
THEN
3753 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) >= this%val_mask(:,:))
3754 field_out(:,:,k) = field_in(:,:,k)
3756 field_out(:,:,k) = rmiss
3760 ELSE IF (this%trans%sub_type ==
'gemaskinvalid')
THEN
3763 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) < this%val_mask(:,:))
3764 field_out(:,:,k) = field_in(:,:,k)
3766 field_out(:,:,k) = rmiss
3770 ELSE IF (this%trans%sub_type ==
'gtmaskinvalid')
THEN
3773 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) <= this%val_mask(:,:))
3774 field_out(:,:,k) = field_in(:,:,k)
3776 field_out(:,:,k) = rmiss
3780 ELSE IF (this%trans%sub_type ==
'setinvalidto')
THEN
3783 WHERE (
c_e(field_in(:,:,k)))
3784 field_out(:,:,k) = field_in(:,:,k)
3786 field_out(:,:,k) = this%val1
3790 ELSE IF (this%trans%sub_type ==
'settoinvalid')
THEN
3791 IF (
c_e(this%val1) .AND.
c_e(this%val2))
THEN
3792 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1 &
3793 .AND. field_in(:,:,:) <= this%val2)
3794 field_out(:,:,:) = rmiss
3796 field_out(:,:,:) = field_in(:,:,:)
3798 ELSE IF (
c_e(this%val1))
THEN
3799 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1)
3800 field_out(:,:,:) = rmiss
3802 field_out(:,:,:) = field_in(:,:,:)
3804 ELSE IF (
c_e(this%val2))
THEN
3805 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) <= this%val2)
3806 field_out(:,:,:) = rmiss
3808 field_out(:,:,:) = field_in(:,:,:)
3813ELSE IF (this%trans%trans_type ==
'vertint')
THEN
3815 IF (this%trans%sub_type ==
'linear')
THEN
3817 alloc_coord_3d_in_act = .false.
3818 IF (
ASSOCIATED(this%inter_index_z))
THEN
3821 IF (
c_e(this%inter_index_z(k)))
THEN
3822 z1 = real(this%inter_zp(k))
3823 z2 = real(1.0d0 - this%inter_zp(k))
3824 SELECT CASE(vartype)
3829 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3830 c_e(field_in(i,j,this%inter_index_z(k)+1)))
THEN
3831 field_out(i,j,k) = &
3832 interp_var_360(field_in(i,j,this%inter_index_z(k)), &
3833 field_in(i,j,this%inter_index_z(k)+1), z1, z2)
3841 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3842 c_e(field_in(i,j,this%inter_index_z(k)+1)) .AND. &
3843 field_in(i,j,this%inter_index_z(k)) > 0. .AND. &
3844 field_in(i,j,this%inter_index_z(k)+1) > 0.)
THEN
3845 field_out(i,j,k) = exp( &
3846 log(field_in(i,j,this%inter_index_z(k)))*z2 + &
3847 log(field_in(i,j,this%inter_index_z(k)+1))*z1)
3855 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3856 c_e(field_in(i,j,this%inter_index_z(k)+1)))
THEN
3857 field_out(i,j,k) = &
3858 field_in(i,j,this%inter_index_z(k))*z2 + &
3859 field_in(i,j,this%inter_index_z(k)+1)*z1
3871 IF (
ALLOCATED(this%coord_3d_in))
THEN
3872 coord_3d_in_act => this%coord_3d_in
3874 CALL l4f_category_log(this%category,l4f_debug, &
3875 "Using external vertical coordinate for vertint")
3878 IF (
PRESENT(coord_3d_in))
THEN
3880 CALL l4f_category_log(this%category,l4f_debug, &
3881 "Using internal vertical coordinate for vertint")
3883 IF (this%dolog)
THEN
3884 ALLOCATE(coord_3d_in_act(
SIZE(coord_3d_in,1), &
3885 SIZE(coord_3d_in,2),
SIZE(coord_3d_in,3)))
3886 alloc_coord_3d_in_act = .true.
3887 WHERE (
c_e(coord_3d_in) .AND. coord_3d_in > 0.0)
3888 coord_3d_in_act = log(coord_3d_in)
3890 coord_3d_in_act = rmiss
3893 coord_3d_in_act => coord_3d_in
3896 CALL l4f_category_log(this%category,l4f_error, &
3897 "Missing internal and external vertical coordinate for vertint")
3903 inused =
SIZE(coord_3d_in_act,3)
3904 IF (inused < 2)
RETURN
3908 IF (.NOT.
c_e(this%vcoord_out(k))) cycle
3912 DO kk = 1, max(inused-kkcache-1, kkcache)
3914 kkdown = kkcache - kk + 1
3916 IF (kkdown >= 1)
THEN
3917 IF (this%vcoord_out(k) >= &
3918 min(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)) .AND. &
3919 this%vcoord_out(k) < &
3920 max(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)))
THEN
3923 kfound = kkcache + this%levshift
3927 IF (kkup < inused)
THEN
3928 IF (this%vcoord_out(k) >= &
3929 min(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)) .AND. &
3930 this%vcoord_out(k) < &
3931 max(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)))
THEN
3934 kfound = kkcache + this%levshift
3941 IF (
c_e(kfound))
THEN
3942 IF (
c_e(field_in(i,j,kfound)) .AND.
c_e(field_in(i,j,kfound+1)))
THEN
3943 z1 = real((this%vcoord_out(k) - coord_3d_in_act(i,j,kfoundin))/ &
3944 (coord_3d_in_act(i,j,kfoundin+1) - coord_3d_in_act(i,j,kfoundin)))
3946 SELECT CASE(vartype)
3949 field_out(i,j,k) = &
3950 interp_var_360(field_in(i,j,kfound), field_in(i,j,kfound+1), z1, z2)
3952 field_out(i,j,k) = exp(log(field_in(i,j,kfound))*z2 + &
3953 log(field_in(i,j,kfound+1))*z1)
3956 field_out(i,j,k) = field_in(i,j,kfound)*z2 + field_in(i,j,kfound+1)*z1
3965 IF (alloc_coord_3d_in_act)
DEALLOCATE(coord_3d_in_act)
3968 ELSE IF (this%trans%sub_type ==
'linearsparse')
THEN
3971 IF (.NOT.
ASSOCIATED(this%vcoord_in) .AND. .NOT.
PRESENT(coord_3d_in))
THEN
3972 CALL l4f_category_log(this%category,l4f_error, &
3973 "linearsparse interpolation, no input vert coord available")
3977 ALLOCATE(coord_in(innz),val_in(innz),mask_in(innz))
3981 IF (
ASSOCIATED(this%vcoord_in))
THEN
3982 mask_in =
c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3983 .AND.
c_e(this%vcoord_in(:))
3985 mask_in =
c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3986 .AND.
c_e(coord_3d_in(i,j,:))
3989 IF (vartype == var_press)
THEN
3990 mask_in(:) = mask_in(:) .AND. &
3991 (field_in(i,j,this%levshift+1:this%levshift+this%levused) > 0.0d0)
3993 inused = count(mask_in)
3994 IF (inused > 1)
THEN
3995 IF (
ASSOCIATED(this%vcoord_in))
THEN
3996 coord_in(1:inused) = pack(this%vcoord_in(:), mask=mask_in)
3998 coord_in(1:inused) = pack(coord_3d_in(i,j,:), mask=mask_in)
4000 IF (vartype == var_press)
THEN
4001 val_in(1:inused) = log(pack( &
4002 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
4005 val_in(1:inused) = pack( &
4006 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
4013 DO kk = 1, max(inused-kkcache-1, kkcache)
4015 kkdown = kkcache - kk + 1
4017 IF (kkdown >= 1)
THEN
4018 IF (this%vcoord_out(k) >= &
4019 min(coord_in(kkdown), coord_in(kkdown+1)) .AND. &
4020 this%vcoord_out(k) < &
4021 max(coord_in(kkdown), coord_in(kkdown+1)))
THEN
4028 IF (kkup < inused)
THEN
4029 IF (this%vcoord_out(k) >= &
4030 min(coord_in(kkup), coord_in(kkup+1)) .AND. &
4031 this%vcoord_out(k) < &
4032 max(coord_in(kkup), coord_in(kkup+1)))
THEN
4042 IF (
c_e(kfound))
THEN
4043 z1 = real((this%vcoord_out(k) - coord_in(kfound-1))/ &
4044 (coord_in(kfound) - coord_in(kfound-1)))
4046 IF (vartype == var_dir360)
THEN
4047 field_out(i,j,k) = &
4048 interp_var_360(val_in(kfound-1), val_in(kfound), z1, z2)
4049 ELSE IF (vartype == var_press)
THEN
4050 field_out(i,j,k) = exp(val_in(kfound-1)*z2 + val_in(kfound)*z1)
4052 field_out(i,j,k) = val_in(kfound-1)*z2 + val_in(kfound)*z1
4062 DEALLOCATE(coord_in,val_in)
4067ELSE IF (this%trans%trans_type ==
'' .OR. this%trans%trans_type ==
'none')
THEN
4069 field_out(:,:,:) = field_in(:,:,:)
4079FUNCTION interp_var_360(v1, v2, w1, w2)
4080REAL,
INTENT(in) :: v1, v2, w1, w2
4081REAL :: interp_var_360
4085IF (abs(v1 - v2) > 180.)
THEN
4093 interp_var_360 = modulo(lv1*w2 + lv2*w1, 360.)
4095 interp_var_360 = v1*w2 + v2*w1
4098END FUNCTION interp_var_360
4101RECURSIVE FUNCTION find_prevailing_direction(v1, l, h, res, mask)
RESULT(prev)
4102REAL,
INTENT(in) :: v1(:,:)
4103REAL,
INTENT(in) :: l, h, res
4104LOGICAL,
INTENT(in),
OPTIONAL :: mask(:,:)
4112IF ((h - l) <= res)
THEN
4117IF (
PRESENT(mask))
THEN
4118 nl = count(v1 >= l .AND. v1 < m .AND. mask)
4119 nh = count(v1 >= m .AND. v1 < h .AND. mask)
4121 nl = count(v1 >= l .AND. v1 < m)
4122 nh = count(v1 >= m .AND. v1 < h)
4125 prev = find_prevailing_direction(v1, m, h, res)
4126ELSE IF (nl > nh)
THEN
4127 prev = find_prevailing_direction(v1, l, m, res)
4128ELSE IF (nl == 0 .AND. nh == 0)
THEN
4134END FUNCTION find_prevailing_direction
4137END SUBROUTINE grid_transform_compute
4145SUBROUTINE grid_transform_v7d_grid_compute(this, field_in, field_out, var, &
4148REAL,
INTENT(in) :: field_in(:,:)
4149REAL,
INTENT(out):: field_out(:,:,:)
4150TYPE(vol7d_var),
INTENT(in),
OPTIONAL :: var
4151REAL,
INTENT(in),
OPTIONAL,
TARGET :: coord_3d_in(:,:,:)
4153real,
allocatable :: field_in_p(:),x_in_p(:),y_in_p(:)
4154INTEGER :: inn_p, ier, k
4155INTEGER :: innx, inny, innz, outnx, outny, outnz
4158call l4f_category_log(this%category,l4f_debug,
"start v7d_grid_transform_compute")
4161field_out(:,:,:) = rmiss
4163IF (.NOT.this%valid)
THEN
4164 call l4f_category_log(this%category,l4f_error, &
4165 "refusing to perform a non valid transformation")
4170innx =
SIZE(field_in,1); inny = 1; innz =
SIZE(field_in,2)
4171outnx =
SIZE(field_out,1); outny =
SIZE(field_out,2); outnz =
SIZE(field_out,3)
4174IF (this%trans%trans_type ==
'vertint')
THEN
4176 IF (innz /= this%innz)
THEN
4177 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
4178 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
4179 t2c(this%innz)//
" /= "//t2c(innz))
4184 IF (outnz /= this%outnz)
THEN
4185 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
4186 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
4187 t2c(this%outnz)//
" /= "//t2c(outnz))
4192 IF (innx /= outnx .OR. inny /= outny)
THEN
4193 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
4194 CALL l4f_category_log(this%category,l4f_error,
"inconsistent hor. sizes: "//&
4195 t2c(innx)//
","//t2c(inny)//
" /= "//&
4196 t2c(outnx)//
","//t2c(outny))
4203 IF (innx /= this%innx .OR. inny /= this%inny)
THEN
4204 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
4205 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
4206 t2c(this%innx)//
","//t2c(this%inny)//
" /= "//&
4207 t2c(innx)//
","//t2c(inny))
4212 IF (outnx /= this%outnx .OR. outny /= this%outny)
THEN
4213 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
4214 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
4215 t2c(this%outnx)//
","//t2c(this%outny)//
" /= "//&
4216 t2c(outnx)//
","//t2c(outny))
4221 IF (innz /= outnz)
THEN
4222 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
4223 CALL l4f_category_log(this%category,l4f_error,
"inconsistent vert. sizes: "//&
4224 t2c(innz)//
" /= "//t2c(outnz))
4232 call l4f_category_log(this%category,l4f_debug, &
4233 "start grid_transform_v7d_grid_compute "//trim(this%trans%trans_type)//
':'// &
4234 trim(this%trans%sub_type))
4237IF (this%trans%trans_type ==
'inter')
THEN
4239 IF (this%trans%sub_type ==
'linear')
THEN
4241#ifdef HAVE_LIBNGMATH
4243 ALLOCATE(field_in_p(innx*inny), x_in_p(innx*inny), y_in_p(innx*inny))
4245 inn_p = count(
c_e(field_in(:,k)))
4247 CALL l4f_category_log(this%category,l4f_info, &
4248 "Number of sparse data points: "//t2c(inn_p)//
','//t2c(
SIZE(field_in(:,k))))
4252 field_in_p(1:inn_p) = pack(field_in(:,k),
c_e(field_in(:,k)))
4253 x_in_p(1:inn_p) = pack(this%inter_xp(:,1),
c_e(field_in(:,k)))
4254 y_in_p(1:inn_p) = pack(this%inter_yp(:,1),
c_e(field_in(:,k)))
4256 IF (.NOT.this%trans%extrap)
THEN
4257 CALL nnseti(
'ext', 0)
4258 CALL nnsetr(
'nul', rmiss)
4261 CALL natgrids(inn_p, x_in_p, y_in_p, field_in_p, &
4262 this%outnx, this%outny, real(this%inter_x(:,1)), &
4263 REAL(this%inter_y(1,:)), field_out(1,1,k), ier)
4266 CALL l4f_category_log(this%category,l4f_error, &
4267 "Error code from NCAR natgrids: "//t2c(ier))
4273 CALL l4f_category_log(this%category,l4f_info, &
4274 "insufficient data in gridded region to triangulate")
4278 DEALLOCATE(field_in_p, x_in_p, y_in_p)
4281 CALL l4f_category_log(this%category,l4f_error, &
4282 "libsim compiled without NATGRIDD (ngmath ncarg library)")
4289ELSE IF (this%trans%trans_type ==
'boxinter' .OR. &
4290 this%trans%trans_type ==
'polyinter' .OR. &
4291 this%trans%trans_type ==
'vertint' .OR. &
4292 this%trans%trans_type ==
'metamorphosis')
THEN
4295 reshape(field_in, (/
SIZE(field_in,1), 1,
SIZE(field_in,2)/)), field_out, var, &
4300END SUBROUTINE grid_transform_v7d_grid_compute
4315ELEMENTAL FUNCTION hbilin (z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
RESULT(zp)
4316REAL,
INTENT(in) :: z1,z2,z3,z4
4317DOUBLE PRECISION,
INTENT(in):: x1,y1