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
2234 DO iy = 1, this%inny
2235 DO ix = 1, this%innx
2236 point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy))
2237 DO n = 1, this%trans%poly%arraysize
2238 IF (
inside(point, this%trans%poly%array(n)))
THEN
2239 this%outnx = this%outnx + 1
2240 this%point_index(ix,iy) = n
2249 IF (this%outnx <= 0)
THEN
2250 CALL l4f_category_log(this%category,l4f_warn, &
2251 "metamorphosis:poly: no points inside polygons")
2254 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2257 DO iy = 1, this%inny
2258 DO ix = 1, this%innx
2259 IF (
c_e(this%point_index(ix,iy)))
THEN
2261 CALL init(v7d_out%ana(n), &
2262 lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
2269 ELSE IF (this%trans%sub_type ==
'mask' )
THEN
2271 IF (.NOT.
PRESENT(maskgrid))
THEN
2272 CALL l4f_category_log(this%category,l4f_error, &
2273 'grid_transform_init maskgrid argument missing for metamorphosis:mask transformation')
2278 IF (this%innx /=
SIZE(maskgrid,1) .OR. this%inny /=
SIZE(maskgrid,2))
THEN
2279 CALL l4f_category_log(this%category,l4f_error, &
2280 'grid_transform_init mask non conformal with input field')
2281 CALL l4f_category_log(this%category,l4f_error, &
2282 'mask: '//t2c(
SIZE(maskgrid,1))//
'x'//t2c(
SIZE(maskgrid,2))// &
2283 ' input field:'//t2c(this%innx)//
'x'//t2c(this%inny))
2289 CALL gen_mask_class()
2297 DO iy = 1, this%inny
2298 DO ix = 1, this%innx
2299 IF (
c_e(maskgrid(ix,iy)))
THEN
2300 IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1))
THEN
2301 DO n = nmaskarea, 1, -1
2302 IF (maskgrid(ix,iy) > lmaskbounds(n))
THEN
2303 this%outnx = this%outnx + 1
2304 this%point_index(ix,iy) = n
2314 IF (this%outnx <= 0)
THEN
2315 CALL l4f_category_log(this%category,l4f_warn, &
2316 "grid_transform_init no points inside mask for metamorphosis:mask transformation")
2320 CALL l4f_category_log(this%category,l4f_info, &
2321 "points in subarea "//t2c(n)//
": "// &
2322 t2c(count(this%point_index(:,:) == n)))
2325 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2328 DO iy = 1, this%inny
2329 DO ix = 1, this%innx
2330 IF (
c_e(this%point_index(ix,iy)))
THEN
2332 CALL init(v7d_out%ana(n),lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
2345SUBROUTINE gen_mask_class()
2346REAL :: startmaskclass, mmin, mmax
2349IF (
PRESENT(maskbounds))
THEN
2350 nmaskarea =
SIZE(maskbounds) - 1
2351 IF (nmaskarea > 0)
THEN
2352 lmaskbounds = maskbounds
2354 ELSE IF (nmaskarea == 0)
THEN
2355 CALL l4f_category_log(this%category,l4f_warn, &
2356 'only one value provided for maskbounds, '//t2c(maskbounds(1)) &
2357 //
', it will be ignored')
2358 CALL l4f_category_log(this%category,l4f_warn, &
2359 'at least 2 values are required for maskbounds')
2363mmin = minval(maskgrid, mask=
c_e(maskgrid))
2364mmax = maxval(maskgrid, mask=
c_e(maskgrid))
2366nmaskarea = int(mmax - mmin + 1.5)
2367startmaskclass = mmin - 0.5
2369ALLOCATE(lmaskbounds(nmaskarea+1))
2370lmaskbounds(:) = (/(startmaskclass+real(i),i=0,nmaskarea)/)
2372CALL l4f_category_log(this%category,l4f_debug, &
2373 'maskinter, '//t2c(nmaskarea)//
' subareas defined, with boundaries:')
2374DO i = 1,
SIZE(lmaskbounds)
2375 CALL l4f_category_log(this%category,l4f_debug, &
2376 'maskinter '//t2c(i)//
' '//t2c(lmaskbounds(i)))
2380END SUBROUTINE gen_mask_class
2382END SUBROUTINE grid_transform_grid_vol7d_init
2403SUBROUTINE grid_transform_vol7d_grid_init(this, trans, v7d_in, out, categoryappend)
2406TYPE(
vol7d),
INTENT(in) :: v7d_in
2408character(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2411DOUBLE PRECISION :: xmin, xmax, ymin, ymax, lonref
2412DOUBLE PRECISION,
ALLOCATABLE :: lon(:,:),lat(:,:)
2416CALL grid_transform_init_common(this, trans, categoryappend)
2418CALL l4f_category_log(this%category, l4f_debug,
"grid_transform v7d-vg6d")
2421IF (this%trans%trans_type ==
'inter')
THEN
2423 IF ( this%trans%sub_type ==
'linear' )
THEN
2425 this%innx=
SIZE(v7d_in%ana)
2427 ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2428 ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2429 CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2431 CALL copy (out, lout)
2433 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
2434 CALL griddim_set_central_lon(lout, lonref)
2436 CALL get_val(lout, nx=nx, ny=ny)
2439 ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2441 CALL get_val(lout, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2442 CALL proj(lout, lon, lat, this%inter_xp, this%inter_yp)
2443 CALL griddim_gen_coord(lout, this%inter_x, this%inter_y)
2452ELSE IF (this%trans%trans_type ==
'boxinter')
THEN
2454 this%innx=
SIZE(v7d_in%ana)
2457 ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2458 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2459 this%inter_index_y(this%innx,this%inny))
2461 CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2463 CALL copy (out, lout)
2465 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
2466 CALL griddim_set_central_lon(lout, lonref)
2468 CALL get_val(lout, nx=this%outnx, ny=this%outny, &
2469 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2472 IF (.NOT.
c_e(this%trans%area_info%boxdx)) &
2473 CALL get_val(out, dx=this%trans%area_info%boxdx)
2474 IF (.NOT.
c_e(this%trans%area_info%boxdy)) &
2475 CALL get_val(out, dx=this%trans%area_info%boxdy)
2477 this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
2478 this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
2481 CALL this%find_index(lout, .true., &
2482 this%outnx, this%outny, xmin, xmax, ymin, ymax, &
2483 lon, lat, .false., &
2484 this%inter_index_x, this%inter_index_y)
2493END SUBROUTINE grid_transform_vol7d_grid_init
2530SUBROUTINE grid_transform_vol7d_vol7d_init(this, trans, v7d_in, v7d_out, &
2531 maskbounds, categoryappend)
2534TYPE(
vol7d),
INTENT(in) :: v7d_in
2535TYPE(
vol7d),
INTENT(inout) :: v7d_out
2536REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
2537CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2540DOUBLE PRECISION,
ALLOCATABLE :: lon(:), lat(:)
2542DOUBLE PRECISION :: lon1, lat1
2546CALL grid_transform_init_common(this, trans, categoryappend)
2548CALL l4f_category_log(this%category, l4f_debug,
"grid_transform v7d-v7d")
2551IF (this%trans%trans_type ==
'inter')
THEN
2553 IF ( this%trans%sub_type ==
'linear' )
THEN
2555 CALL vol7d_alloc_vol(v7d_out)
2556 this%outnx=
SIZE(v7d_out%ana)
2559 this%innx=
SIZE(v7d_in%ana)
2562 ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
3108 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3109 this%trans%sub_type ==
'stddevnm1')
THEN
3111 IF (this%trans%sub_type ==
'stddev')
THEN
3117 navg = this%trans%box_info%npx*this%trans%box_info%npy
3120 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3121 je = j+this%trans%box_info%npy-1
3124 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3125 ie = i+this%trans%box_info%npx-1
3128 reshape(field_in(i:ie,j:je,k),(/navg/)), nm1=nm1)
3133 ELSE IF (this%trans%sub_type ==
'max')
THEN
3136 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3137 je = j+this%trans%box_info%npy-1
3140 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3141 ie = i+this%trans%box_info%npx-1
3143 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3145 field_out(ii,jj,k) = maxval(field_in(i:ie,j:je,k), &
3146 mask=(field_in(i:ie,j:je,k) /= rmiss))
3152 ELSE IF (this%trans%sub_type ==
'min')
THEN
3155 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3156 je = j+this%trans%box_info%npy-1
3159 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3160 ie = i+this%trans%box_info%npx-1
3162 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3164 field_out(ii,jj,k) = minval(field_in(i:ie,j:je,k), &
3165 mask=(field_in(i:ie,j:je,k) /= rmiss))
3171 ELSE IF (this%trans%sub_type ==
'percentile')
THEN
3173 navg = this%trans%box_info%npx*this%trans%box_info%npy
3176 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3177 je = j+this%trans%box_info%npy-1
3180 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3181 ie = i+this%trans%box_info%npx-1
3184 reshape(field_in(i:ie,j:je,k),(/navg/)), &
3185 (/real(this%trans%stat_info%percentile)/))
3190 ELSE IF (this%trans%sub_type ==
'frequency')
THEN
3194 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3195 je = j+this%trans%box_info%npy-1
3198 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3199 ie = i+this%trans%box_info%npx-1
3201 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3203 field_out(ii,jj,k) = sum(interval_info_valid( &
3204 this%trans%interval_info, field_in(i:ie,j:je,k)), &
3205 mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
3213ELSE IF (this%trans%trans_type ==
'inter')
THEN
3215 IF (this%trans%sub_type ==
'near')
THEN
3218 DO j = 1, this%outny
3219 DO i = 1, this%outnx
3221 IF (
c_e(this%inter_index_x(i,j))) field_out(i,j,k) = &
3222 field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3228 ELSE IF (this%trans%sub_type ==
'bilin')
THEN
3231 DO j = 1, this%outny
3232 DO i = 1, this%outnx
3234 IF (
c_e(this%inter_index_x(i,j)))
THEN
3236 z1=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3237 z2=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j),k)
3238 z3=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1,k)
3239 z4=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j)+1,k)
3241 IF (
c_e(z1) .AND.
c_e(z2) .AND.
c_e(z3) .AND.
c_e(z4))
THEN
3243 x1=this%inter_x(this%inter_index_x(i,j),this%inter_index_y(i,j))
3244 y1=this%inter_y(this%inter_index_x(i,j),this%inter_index_y(i,j))
3245 x3=this%inter_x(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3246 y3=this%inter_y(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3248 xp=this%inter_xp(i,j)
3249 yp=this%inter_yp(i,j)
3251 field_out(i,j,k) = hbilin(z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
3259 ELSE IF (this%trans%sub_type ==
'shapiro_near')
THEN
3261 DO j = 1, this%outny
3262 DO i = 1, this%outnx
3264 IF (
c_e(this%inter_index_x(i,j)))
THEN
3266 IF(this%inter_index_x(i,j)-1>0)
THEN
3267 z(1) = field_in(this%inter_index_x(i,j)-1,this%inter_index_y(i,j) ,k)
3271 IF(this%inter_index_x(i,j)+1<=this%outnx)
THEN
3272 z(3)=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j) ,k)
3276 IF(this%inter_index_y(i,j)+1<=this%outny)
THEN
3277 z(2)=field_in(this%inter_index_x(i,j) ,this%inter_index_y(i,j)+1,k)
3281 IF(this%inter_index_y(i,j)-1>0)
THEN
3282 z(4)=field_in(this%inter_index_x(i,j), this%inter_index_y(i,j)-1,k)
3286 field_out(i,j,k) = shapiro(z,field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k))
3295ELSE IF (this%trans%trans_type ==
'intersearch')
THEN
3298 likethis%trans%trans_type =
'inter'
3299 CALL grid_transform_compute(likethis, field_in, field_out, var, coord_3d_in)
3300 CALL getenv(
'LIBSIM_DISABLEOPTSEARCH', env_var)
3301 optsearch = len_trim(env_var) == 0
3304 IF ((.NOT.all(
c_e(field_out(:,:,k)))) .AND. (any(
c_e(field_in(:,:,k)))))
THEN
3305 DO j = 1, this%outny
3306 DO i = 1, this%outnx
3307 IF (.NOT.
c_e(field_out(i,j,k)))
THEN
3311 ix = this%inter_index_x(i,j)
3312 iy = this%inter_index_y(i,j)
3313 DO s = 0, max(this%innx, this%inny)
3315 DO m = iy-s, iy+s, max(2*s, 1)
3316 IF (m < 1 .OR. m > this%inny) cycle
3317 DO l = max(1, ix-s), min(this%innx, ix+s)
3318 disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2
3319 IF (
c_e(field_in(l,m,k)))
THEN
3320 IF (disttmp <
dist)
THEN
3322 field_out(i,j,k) = field_in(l,m,k)
3324 ELSE IF (disttmp ==
dist)
THEN
3325 field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k)
3326 nearcount = nearcount + 1
3329 IF (disttmp <
dist) farenough = .false.
3332 DO m = max(1, iy-s+1), min(this%inny, iy+s-1)
3333 DO l = ix-s, ix+s, 2*s
3334 IF (l < 1 .OR. l > this%innx) cycle
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 IF (s > 0 .AND. farenough)
EXIT
3354 IF (
c_e(field_in(l,m,k)))
THEN
3355 disttmp = (this%inter_xp(i,j) - this%inter_x(l,m))**2 + (this%inter_yp(i,j) - this%inter_y(l,m))**2
3356 IF (disttmp <
dist)
THEN
3358 field_out(i,j,k) = field_in(l,m,k)
3360 ELSE IF (disttmp ==
dist)
THEN
3361 field_out(i,j,k) = field_out(i,j,k) + field_in(l,m,k)
3362 nearcount = nearcount + 1
3369 IF (nearcount > 1) field_out(i,j,k) = field_out(i,j,k)/nearcount
3376ELSE IF (this%trans%trans_type ==
'boxinter' &
3377 .OR. this%trans%trans_type ==
'polyinter' &
3378 .OR. this%trans%trans_type ==
'maskinter')
THEN
3380 IF (this%trans%sub_type ==
'average')
THEN
3382 IF (vartype == var_dir360)
THEN
3384 DO j = 1, this%outny
3385 DO i = 1, this%outnx
3386 field_out(i,j,k) = find_prevailing_direction(field_in(:,:,k), &
3388 mask=this%inter_index_x(:,:) == i .AND. this%inter_index_y(:,:) == j)
3394 ALLOCATE(nval(this%outnx, this%outny))
3395 field_out(:,:,:) = 0.0
3400 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3401 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3402 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3404 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3405 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3409 WHERE (nval(:,:) /= 0)
3410 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3412 field_out(:,:,k) = rmiss
3418 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3419 this%trans%sub_type ==
'stddevnm1')
THEN
3421 IF (this%trans%sub_type ==
'stddev')
THEN
3427 DO j = 1, this%outny
3428 DO i = 1, this%outnx
3431 reshape(field_in(:,:,k), (/
SIZE(field_in(:,:,k))/)), &
3432 mask=reshape((this%inter_index_x == i .AND. &
3433 this%inter_index_y == j), (/
SIZE(field_in(:,:,k))/)), nm1=nm1)
3438 ELSE IF (this%trans%sub_type ==
'max')
THEN
3443 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3444 IF (
c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k)))
THEN
3445 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3446 max(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3449 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3458 ELSE IF (this%trans%sub_type ==
'min')
THEN
3463 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3464 IF (
c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k)))
THEN
3465 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3466 min(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3469 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3477 ELSE IF (this%trans%sub_type ==
'percentile')
THEN
3480 DO j = 1, this%outny
3481 DO i = 1, this%outnx
3484 reshape(field_in(:,:,k), (/
SIZE(field_in(:,:,k))/)), &
3485 (/real(this%trans%stat_info%percentile)/), &
3486 mask=reshape((this%inter_index_x == i .AND. &
3487 this%inter_index_y == j), (/
SIZE(field_in(:,:,k))/)))
3492 ELSE IF (this%trans%sub_type ==
'frequency')
THEN
3494 ALLOCATE(nval(this%outnx, this%outny))
3495 field_out(:,:,:) = 0.0
3500 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3501 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3502 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3503 interval_info_valid(this%trans%interval_info, field_in(i,j,k))
3504 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3505 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3509 WHERE (nval(:,:) /= 0)
3510 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3512 field_out(:,:,k) = rmiss
3519ELSE IF (this%trans%trans_type ==
'stencilinter')
THEN
3520 np =
SIZE(this%stencil,1)/2
3521 ns =
SIZE(this%stencil)
3523 IF (this%trans%sub_type ==
'average')
THEN
3525 IF (vartype == var_dir360)
THEN
3527 DO j = 1, this%outny
3528 DO i = 1, this%outnx
3529 IF (
c_e(this%inter_index_x(i,j)))
THEN
3530 i1 = this%inter_index_x(i,j) - np
3531 i2 = this%inter_index_x(i,j) + np
3532 j1 = this%inter_index_y(i,j) - np
3533 j2 = this%inter_index_y(i,j) + np
3534 field_out(i,j,k) = find_prevailing_direction(field_in(i1:i2,j1:j2,k), &
3536 mask=this%stencil(:,:))
3547 DO j = 1, this%outny
3548 DO i = 1, this%outnx
3549 IF (
c_e(this%inter_index_x(i,j)))
THEN
3550 i1 = this%inter_index_x(i,j) - np
3551 i2 = this%inter_index_x(i,j) + np
3552 j1 = this%inter_index_y(i,j) - np
3553 j2 = this%inter_index_y(i,j) + np
3554 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3556 field_out(i,j,k) = sum(field_in(i1:i2,j1:j2,k), &
3557 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3566 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3567 this%trans%sub_type ==
'stddevnm1')
THEN
3569 IF (this%trans%sub_type ==
'stddev')
THEN
3578 DO j = 1, this%outny
3579 DO i = 1, this%outnx
3580 IF (
c_e(this%inter_index_x(i,j)))
THEN
3581 i1 = this%inter_index_x(i,j) - np
3582 i2 = this%inter_index_x(i,j) + np
3583 j1 = this%inter_index_y(i,j) - np
3584 j2 = this%inter_index_y(i,j) + np
3587 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3588 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3589 this%stencil(:,:), (/ns/)), nm1=nm1)
3596 ELSE IF (this%trans%sub_type ==
'max')
THEN
3601 DO j = 1, this%outny
3602 DO i = 1, this%outnx
3603 IF (
c_e(this%inter_index_x(i,j)))
THEN
3604 i1 = this%inter_index_x(i,j) - np
3605 i2 = this%inter_index_x(i,j) + np
3606 j1 = this%inter_index_y(i,j) - np
3607 j2 = this%inter_index_y(i,j) + np
3608 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3610 field_out(i,j,k) = maxval(field_in(i1:i2,j1:j2,k), &
3611 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3619 ELSE IF (this%trans%sub_type ==
'min')
THEN
3624 DO j = 1, this%outny
3625 DO i = 1, this%outnx
3626 IF (
c_e(this%inter_index_x(i,j)))
THEN
3627 i1 = this%inter_index_x(i,j) - np
3628 i2 = this%inter_index_x(i,j) + np
3629 j1 = this%inter_index_y(i,j) - np
3630 j2 = this%inter_index_y(i,j) + np
3631 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3633 field_out(i,j,k) = minval(field_in(i1:i2,j1:j2,k), &
3634 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3642 ELSE IF (this%trans%sub_type ==
'percentile')
THEN
3647 DO j = 1, this%outny
3648 DO i = 1, this%outnx
3649 IF (
c_e(this%inter_index_x(i,j)))
THEN
3650 i1 = this%inter_index_x(i,j) - np
3651 i2 = this%inter_index_x(i,j) + np
3652 j1 = this%inter_index_y(i,j) - np
3653 j2 = this%inter_index_y(i,j) + np
3656 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3657 (/real(this%trans%stat_info%percentile)/), &
3658 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3659 this%stencil(:,:), (/ns/)))
3666 ELSE IF (this%trans%sub_type ==
'frequency')
THEN
3671 DO j = 1, this%outny
3672 DO i = 1, this%outnx
3673 IF (
c_e(this%inter_index_x(i,j)))
THEN
3674 i1 = this%inter_index_x(i,j) - np
3675 i2 = this%inter_index_x(i,j) + np
3676 j1 = this%inter_index_y(i,j) - np
3677 j2 = this%inter_index_y(i,j) + np
3678 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3680 field_out(i,j,k) = sum(interval_info_valid( &
3681 this%trans%interval_info, field_in(i1:i2,j1:j2,k)), &
3682 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3692ELSE IF (this%trans%trans_type ==
'maskgen')
THEN
3695 WHERE(
c_e(this%inter_index_x(:,:)))
3696 field_out(:,:,k) = real(this%inter_index_x(:,:))
3700ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN
3702 IF (this%trans%sub_type ==
'all')
THEN
3704 field_out(:,:,:) = reshape(field_in(:,:,:), (/this%outnx,this%outny,innz/))
3706 ELSE IF (this%trans%sub_type ==
'coordbb' .OR. this%trans%sub_type ==
'poly' &
3707 .OR. this%trans%sub_type ==
'mask')
THEN
3711 field_out(:,1,k) = pack(field_in(:,:,k),
c_e(this%point_index(:,:)))
3714 ELSE IF (this%trans%sub_type ==
'maskvalid' .OR. &
3715 this%trans%sub_type ==
'maskinvalid')
THEN
3718 WHERE (this%point_mask(:,:))
3719 field_out(:,:,k) = field_in(:,:,k)
3723 ELSE IF (this%trans%sub_type ==
'lemaskinvalid')
THEN
3726 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) > this%val_mask(:,:))
3727 field_out(:,:,k) = field_in(:,:,k)
3729 field_out(:,:,k) = rmiss
3733 ELSE IF (this%trans%sub_type ==
'ltmaskinvalid')
THEN
3736 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) >= this%val_mask(:,:))
3737 field_out(:,:,k) = field_in(:,:,k)
3739 field_out(:,:,k) = rmiss
3743 ELSE IF (this%trans%sub_type ==
'gemaskinvalid')
THEN
3746 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) < this%val_mask(:,:))
3747 field_out(:,:,k) = field_in(:,:,k)
3749 field_out(:,:,k) = rmiss
3753 ELSE IF (this%trans%sub_type ==
'gtmaskinvalid')
THEN
3756 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) <= this%val_mask(:,:))
3757 field_out(:,:,k) = field_in(:,:,k)
3759 field_out(:,:,k) = rmiss
3763 ELSE IF (this%trans%sub_type ==
'setinvalidto')
THEN
3766 WHERE (
c_e(field_in(:,:,k)))
3767 field_out(:,:,k) = field_in(:,:,k)
3769 field_out(:,:,k) = this%val1
3773 ELSE IF (this%trans%sub_type ==
'settoinvalid')
THEN
3774 IF (
c_e(this%val1) .AND.
c_e(this%val2))
THEN
3775 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1 &
3776 .AND. field_in(:,:,:) <= this%val2)
3777 field_out(:,:,:) = rmiss
3779 field_out(:,:,:) = field_in(:,:,:)
3781 ELSE IF (
c_e(this%val1))
THEN
3782 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1)
3783 field_out(:,:,:) = rmiss
3785 field_out(:,:,:) = field_in(:,:,:)
3787 ELSE IF (
c_e(this%val2))
THEN
3788 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) <= this%val2)
3789 field_out(:,:,:) = rmiss
3791 field_out(:,:,:) = field_in(:,:,:)
3796ELSE IF (this%trans%trans_type ==
'vertint')
THEN
3798 IF (this%trans%sub_type ==
'linear')
THEN
3800 alloc_coord_3d_in_act = .false.
3801 IF (
ASSOCIATED(this%inter_index_z))
THEN
3804 IF (
c_e(this%inter_index_z(k)))
THEN
3805 z1 = real(this%inter_zp(k))
3806 z2 = real(1.0d0 - this%inter_zp(k))
3807 SELECT CASE(vartype)
3812 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3813 c_e(field_in(i,j,this%inter_index_z(k)+1)))
THEN
3814 field_out(i,j,k) = &
3815 interp_var_360(field_in(i,j,this%inter_index_z(k)), &
3816 field_in(i,j,this%inter_index_z(k)+1), z1, z2)
3824 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3825 c_e(field_in(i,j,this%inter_index_z(k)+1)) .AND. &
3826 field_in(i,j,this%inter_index_z(k)) > 0. .AND. &
3827 field_in(i,j,this%inter_index_z(k)+1) > 0.)
THEN
3828 field_out(i,j,k) = exp( &
3829 log(field_in(i,j,this%inter_index_z(k)))*z2 + &
3830 log(field_in(i,j,this%inter_index_z(k)+1))*z1)
3838 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3839 c_e(field_in(i,j,this%inter_index_z(k)+1)))
THEN
3840 field_out(i,j,k) = &
3841 field_in(i,j,this%inter_index_z(k))*z2 + &
3842 field_in(i,j,this%inter_index_z(k)+1)*z1
3854 IF (
ALLOCATED(this%coord_3d_in))
THEN
3855 coord_3d_in_act => this%coord_3d_in
3857 CALL l4f_category_log(this%category,l4f_debug, &
3858 "Using external vertical coordinate for vertint")
3861 IF (
PRESENT(coord_3d_in))
THEN
3863 CALL l4f_category_log(this%category,l4f_debug, &
3864 "Using internal vertical coordinate for vertint")
3866 IF (this%dolog)
THEN
3867 ALLOCATE(coord_3d_in_act(
SIZE(coord_3d_in,1), &
3868 SIZE(coord_3d_in,2),
SIZE(coord_3d_in,3)))
3869 alloc_coord_3d_in_act = .true.
3870 WHERE (
c_e(coord_3d_in) .AND. coord_3d_in > 0.0)
3871 coord_3d_in_act = log(coord_3d_in)
3873 coord_3d_in_act = rmiss
3876 coord_3d_in_act => coord_3d_in
3879 CALL l4f_category_log(this%category,l4f_error, &
3880 "Missing internal and external vertical coordinate for vertint")
3886 inused =
SIZE(coord_3d_in_act,3)
3887 IF (inused < 2)
RETURN
3891 IF (.NOT.
c_e(this%vcoord_out(k))) cycle
3895 DO kk = 1, max(inused-kkcache-1, kkcache)
3897 kkdown = kkcache - kk + 1
3899 IF (kkdown >= 1)
THEN
3900 IF (this%vcoord_out(k) >= &
3901 min(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)) .AND. &
3902 this%vcoord_out(k) < &
3903 max(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)))
THEN
3906 kfound = kkcache + this%levshift
3910 IF (kkup < inused)
THEN
3911 IF (this%vcoord_out(k) >= &
3912 min(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)) .AND. &
3913 this%vcoord_out(k) < &
3914 max(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)))
THEN
3917 kfound = kkcache + this%levshift
3924 IF (
c_e(kfound))
THEN
3925 IF (
c_e(field_in(i,j,kfound)) .AND.
c_e(field_in(i,j,kfound+1)))
THEN
3926 z1 = real((this%vcoord_out(k) - coord_3d_in_act(i,j,kfoundin))/ &
3927 (coord_3d_in_act(i,j,kfoundin+1) - coord_3d_in_act(i,j,kfoundin)))
3929 SELECT CASE(vartype)
3932 field_out(i,j,k) = &
3933 interp_var_360(field_in(i,j,kfound), field_in(i,j,kfound+1), z1, z2)
3935 field_out(i,j,k) = exp(log(field_in(i,j,kfound))*z2 + &
3936 log(field_in(i,j,kfound+1))*z1)
3939 field_out(i,j,k) = field_in(i,j,kfound)*z2 + field_in(i,j,kfound+1)*z1
3948 IF (alloc_coord_3d_in_act)
DEALLOCATE(coord_3d_in_act)
3951 ELSE IF (this%trans%sub_type ==
'linearsparse')
THEN
3954 IF (.NOT.
ASSOCIATED(this%vcoord_in) .AND. .NOT.
PRESENT(coord_3d_in))
THEN
3955 CALL l4f_category_log(this%category,l4f_error, &
3956 "linearsparse interpolation, no input vert coord available")
3960 ALLOCATE(coord_in(innz),val_in(innz),mask_in(innz))
3964 IF (
ASSOCIATED(this%vcoord_in))
THEN
3965 mask_in =
c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3966 .AND.
c_e(this%vcoord_in(:))
3968 mask_in =
c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3969 .AND.
c_e(coord_3d_in(i,j,:))
3972 IF (vartype == var_press)
THEN
3973 mask_in(:) = mask_in(:) .AND. &
3974 (field_in(i,j,this%levshift+1:this%levshift+this%levused) > 0.0d0)
3976 inused = count(mask_in)
3977 IF (inused > 1)
THEN
3978 IF (
ASSOCIATED(this%vcoord_in))
THEN
3979 coord_in(1:inused) = pack(this%vcoord_in(:), mask=mask_in)
3981 coord_in(1:inused) = pack(coord_3d_in(i,j,:), mask=mask_in)
3983 IF (vartype == var_press)
THEN
3984 val_in(1:inused) = log(pack( &
3985 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3988 val_in(1:inused) = pack( &
3989 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3996 DO kk = 1, max(inused-kkcache-1, kkcache)
3998 kkdown = kkcache - kk + 1
4000 IF (kkdown >= 1)
THEN
4001 IF (this%vcoord_out(k) >= &
4002 min(coord_in(kkdown), coord_in(kkdown+1)) .AND. &
4003 this%vcoord_out(k) < &
4004 max(coord_in(kkdown), coord_in(kkdown+1)))
THEN
4011 IF (kkup < inused)
THEN
4012 IF (this%vcoord_out(k) >= &
4013 min(coord_in(kkup), coord_in(kkup+1)) .AND. &
4014 this%vcoord_out(k) < &
4015 max(coord_in(kkup), coord_in(kkup+1)))
THEN
4025 IF (
c_e(kfound))
THEN
4026 z1 = real((this%vcoord_out(k) - coord_in(kfound-1))/ &
4027 (coord_in(kfound) - coord_in(kfound-1)))
4029 IF (vartype == var_dir360)
THEN
4030 field_out(i,j,k) = &
4031 interp_var_360(val_in(kfound-1), val_in(kfound), z1, z2)
4032 ELSE IF (vartype == var_press)
THEN
4033 field_out(i,j,k) = exp(val_in(kfound-1)*z2 + val_in(kfound)*z1)
4035 field_out(i,j,k) = val_in(kfound-1)*z2 + val_in(kfound)*z1
4045 DEALLOCATE(coord_in,val_in)
4050ELSE IF (this%trans%trans_type ==
'' .OR. this%trans%trans_type ==
'none')
THEN
4052 field_out(:,:,:) = field_in(:,:,:)
4062FUNCTION interp_var_360(v1, v2, w1, w2)
4063REAL,
INTENT(in) :: v1, v2, w1, w2
4064REAL :: interp_var_360
4068IF (abs(v1 - v2) > 180.)
THEN
4076 interp_var_360 = modulo(lv1*w2 + lv2*w1, 360.)
4078 interp_var_360 = v1*w2 + v2*w1
4081END FUNCTION interp_var_360
4084RECURSIVE FUNCTION find_prevailing_direction(v1, l, h, res, mask)
RESULT(prev)
4085REAL,
INTENT(in) :: v1(:,:)
4086REAL,
INTENT(in) :: l, h, res
4087LOGICAL,
INTENT(in),
OPTIONAL :: mask(:,:)
4095IF ((h - l) <= res)
THEN
4100IF (
PRESENT(mask))
THEN
4101 nl = count(v1 >= l .AND. v1 < m .AND. mask)
4102 nh = count(v1 >= m .AND. v1 < h .AND. mask)
4104 nl = count(v1 >= l .AND. v1 < m)
4105 nh = count(v1 >= m .AND. v1 < h)
4108 prev = find_prevailing_direction(v1, m, h, res)
4109ELSE IF (nl > nh)
THEN
4110 prev = find_prevailing_direction(v1, l, m, res)
4111ELSE IF (nl == 0 .AND. nh == 0)
THEN
4117END FUNCTION find_prevailing_direction
4120END SUBROUTINE grid_transform_compute
4128SUBROUTINE grid_transform_v7d_grid_compute(this, field_in, field_out, var, &
4131REAL,
INTENT(in) :: field_in(:,:)
4132REAL,
INTENT(out):: field_out(:,:,:)
4133TYPE(vol7d_var),
INTENT(in),
OPTIONAL :: var
4134REAL,
INTENT(in),
OPTIONAL,
TARGET :: coord_3d_in(:,:,:)
4136real,
allocatable :: field_in_p(:),x_in_p(:),y_in_p(:)
4137INTEGER :: inn_p, ier, k
4138INTEGER :: innx, inny, innz, outnx, outny, outnz
4141call l4f_category_log(this%category,l4f_debug,
"start v7d_grid_transform_compute")
4144field_out(:,:,:) = rmiss
4146IF (.NOT.this%valid)
THEN
4147 call l4f_category_log(this%category,l4f_error, &
4148 "refusing to perform a non valid transformation")
4153innx =
SIZE(field_in,1); inny = 1; innz =
SIZE(field_in,2)
4154outnx =
SIZE(field_out,1); outny =
SIZE(field_out,2); outnz =
SIZE(field_out,3)
4157IF (this%trans%trans_type ==
'vertint')
THEN
4159 IF (innz /= this%innz)
THEN
4160 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
4161 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
4162 t2c(this%innz)//
" /= "//t2c(innz))
4167 IF (outnz /= this%outnz)
THEN
4168 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
4169 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
4170 t2c(this%outnz)//
" /= "//t2c(outnz))
4175 IF (innx /= outnx .OR. inny /= outny)
THEN
4176 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
4177 CALL l4f_category_log(this%category,l4f_error,
"inconsistent hor. sizes: "//&
4178 t2c(innx)//
","//t2c(inny)//
" /= "//&
4179 t2c(outnx)//
","//t2c(outny))
4186 IF (innx /= this%innx .OR. inny /= this%inny)
THEN
4187 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
4188 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
4189 t2c(this%innx)//
","//t2c(this%inny)//
" /= "//&
4190 t2c(innx)//
","//t2c(inny))
4195 IF (outnx /= this%outnx .OR. outny /= this%outny)
THEN
4196 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
4197 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
4198 t2c(this%outnx)//
","//t2c(this%outny)//
" /= "//&
4199 t2c(outnx)//
","//t2c(outny))
4204 IF (innz /= outnz)
THEN
4205 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
4206 CALL l4f_category_log(this%category,l4f_error,
"inconsistent vert. sizes: "//&
4207 t2c(innz)//
" /= "//t2c(outnz))
4215 call l4f_category_log(this%category,l4f_debug, &
4216 "start grid_transform_v7d_grid_compute "//trim(this%trans%trans_type)//
':'// &
4217 trim(this%trans%sub_type))
4220IF (this%trans%trans_type ==
'inter')
THEN
4222 IF (this%trans%sub_type ==
'linear')
THEN
4224#ifdef HAVE_LIBNGMATH
4226 ALLOCATE(field_in_p(innx*inny), x_in_p(innx*inny), y_in_p(innx*inny))
4228 inn_p = count(
c_e(field_in(:,k)))
4230 CALL l4f_category_log(this%category,l4f_info, &
4231 "Number of sparse data points: "//t2c(inn_p)//
','//t2c(
SIZE(field_in(:,k))))
4235 field_in_p(1:inn_p) = pack(field_in(:,k),
c_e(field_in(:,k)))
4236 x_in_p(1:inn_p) = pack(this%inter_xp(:,1),
c_e(field_in(:,k)))
4237 y_in_p(1:inn_p) = pack(this%inter_yp(:,1),
c_e(field_in(:,k)))
4239 IF (.NOT.this%trans%extrap)
THEN
4240 CALL nnseti(
'ext', 0)
4241 CALL nnsetr(
'nul', rmiss)
4244 CALL natgrids(inn_p, x_in_p, y_in_p, field_in_p, &
4245 this%outnx, this%outny, real(this%inter_x(:,1)), &
4246 REAL(this%inter_y(1,:)), field_out(1,1,k), ier)
4249 CALL l4f_category_log(this%category,l4f_error, &
4250 "Error code from NCAR natgrids: "//t2c(ier))
4256 CALL l4f_category_log(this%category,l4f_info, &
4257 "insufficient data in gridded region to triangulate")
4261 DEALLOCATE(field_in_p, x_in_p, y_in_p)
4264 CALL l4f_category_log(this%category,l4f_error, &
4265 "libsim compiled without NATGRIDD (ngmath ncarg library)")
4272ELSE IF (this%trans%trans_type ==
'boxinter' .OR. &
4273 this%trans%trans_type ==
'polyinter' .OR. &
4274 this%trans%trans_type ==
'vertint' .OR. &
4275 this%trans%trans_type ==
'metamorphosis')
THEN
4278 reshape(field_in, (/
SIZE(field_in,1), 1,
SIZE(field_in,2)/)), field_out, var, &
4283END SUBROUTINE grid_transform_v7d_grid_compute
4298ELEMENTAL FUNCTION hbilin (z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
RESULT(zp)
4299REAL,
INTENT(in) :: z1,z2,z3,z4
4300DOUBLE PRECISION,
INTENT(in):: x1,y1