1409 CALL this%find_index(in, .false., &
1410 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1411 lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1412 this%inter_index_x, this%inter_index_y)
1414 ALLOCATE(this%inter_x(this%innx,this%inny), &
1415 this%inter_y(this%innx,this%inny))
1416 ALLOCATE(this%inter_xp(this%outnx,this%outny), &
1417 this%inter_yp(this%outnx,this%outny))
1420 CALL griddim_gen_coord(in, this%inter_x, this%inter_y)
1422 CALL proj(in, lout%dim%lon, lout%dim%lat, this%inter_xp, this%inter_yp)
1425 CALL this%find_index(in, .true., &
1426 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1427 lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1428 this%inter_index_x, this%inter_index_y)
1430 IF (this%trans%trans_type ==
'intersearch')
THEN
1431 ALLOCATE(this%inter_x(this%innx,this%inny), &
1432 this%inter_y(this%innx,this%inny))
1433 ALLOCATE(this%inter_xp(this%outnx,this%outny), &
1434 this%inter_yp(this%outnx,this%outny))
1437 CALL griddim_gen_coord(in, this%inter_x, this%inter_y)
1439 CALL proj(in, lout%dim%lon, lout%dim%lat, this%inter_xp, this%inter_yp)
1447ELSE IF (this%trans%trans_type ==
'boxinter')
THEN
1449 CALL outgrid_setup()
1450 CALL get_val(in, nx=this%innx, ny=this%inny)
1451 CALL get_val(out, nx=this%outnx, ny=this%outny, &
1452 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1455 IF (.NOT.
c_e(this%trans%area_info%boxdx)) &
1456 CALL get_val(out, dx=this%trans%area_info%boxdx)
1457 IF (.NOT.
c_e(this%trans%area_info%boxdy)) &
1458 CALL get_val(out, dx=this%trans%area_info%boxdy)
1460 this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
1461 this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
1463 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1464 this%inter_index_y(this%innx,this%inny))
1470 CALL this%find_index(out, .true., &
1471 this%outnx, this%outny, xmin, xmax, ymin, ymax, &
1472 lin%dim%lon, lin%dim%lat, .false., &
1473 this%inter_index_x, this%inter_index_y)
1478ELSE IF (this%trans%trans_type ==
'stencilinter')
THEN
1480 CALL outgrid_setup()
1482 CALL get_val(in, nx=this%innx, ny=this%inny, &
1483 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1484 CALL get_val(out, nx=this%outnx, ny=this%outny)
1486 ALLOCATE (this%inter_index_x(this%outnx,this%outny), &
1487 this%inter_index_y(this%outnx,this%outny))
1488 CALL copy(out, lout)
1490 CALL this%find_index(in, .true., &
1491 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1492 lout%dim%lon, lout%dim%lat, this%trans%extrap, &
1493 this%inter_index_x, this%inter_index_y)
1496 nr = int(this%trans%area_info%radius)
1499 r2 = this%trans%area_info%radius**2
1500 ALLOCATE(this%stencil(n,n))
1501 this%stencil(:,:) = .true.
1504 IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
1511 xnmax = this%innx - nr
1513 ynmax = this%inny - nr
1514 DO iy = 1, this%outny
1515 DO ix = 1, this%outnx
1516 IF (this%inter_index_x(ix,iy) < xnmin .OR. &
1517 this%inter_index_x(ix,iy) > xnmax .OR. &
1518 this%inter_index_y(ix,iy) < ynmin .OR. &
1519 this%inter_index_y(ix,iy) > ynmax)
THEN
1520 this%inter_index_x(ix,iy) = imiss
1521 this%inter_index_y(ix,iy) = imiss
1527 CALL l4f_category_log(this%category, l4f_debug, &
1528 'stencilinter: stencil size '//t2c(n*n))
1529 CALL l4f_category_log(this%category, l4f_debug, &
1530 'stencilinter: stencil points '//t2c(count(this%stencil)))
1536ELSE IF (this%trans%trans_type ==
'maskgen')
THEN
1538 IF (this%trans%sub_type ==
'poly')
THEN
1541 CALL get_val(in, nx=this%innx, ny=this%inny)
1542 this%outnx = this%innx
1543 this%outny = this%inny
1546 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1547 this%inter_index_y(this%innx,this%inny))
1548 this%inter_index_x(:,:) = imiss
1549 this%inter_index_y(:,:) = 1
1558 inside_x:
DO i = 1, this%innx
1559 point = georef_coord_new(x=out%dim%lon(i,j), y=out%dim%lat(i,j))
1561 DO n = nprev, this%trans%poly%arraysize
1562 IF (
inside(point, this%trans%poly%array(n)))
THEN
1563 this%inter_index_x(i,j) = n
1568 DO n = nprev-1, 1, -1
1569 IF (
inside(point, this%trans%poly%array(n)))
THEN
1570 this%inter_index_x(i,j) = n
1581 ELSE IF (this%trans%sub_type ==
'grid')
THEN
1584 CALL copy(out, lout)
1587 CALL get_val(in, nx=this%innx, ny=this%inny)
1588 this%outnx = this%innx
1589 this%outny = this%inny
1590 CALL get_val(lout, nx=nx, ny=ny, &
1591 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1594 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1595 this%inter_index_y(this%innx,this%inny))
1601 CALL this%find_index(lout, .true., &
1602 nx, ny, xmin, xmax, ymin, ymax, &
1603 out%dim%lon, out%dim%lat, .false., &
1604 this%inter_index_x, this%inter_index_y)
1606 WHERE(
c_e(this%inter_index_x(:,:)))
1607 this%inter_index_x(:,:) = this%inter_index_x(:,:) + &
1608 (this%inter_index_y(:,:)-1)*nx
1616ELSE IF (this%trans%trans_type ==
'polyinter')
THEN
1622 CALL get_val(in, nx=this%innx, ny=this%inny)
1623 this%outnx = this%innx
1624 this%outny = this%inny
1627 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1628 this%inter_index_y(this%innx,this%inny))
1629 this%inter_index_x(:,:) = imiss
1630 this%inter_index_y(:,:) = 1
1639 inside_x_2:
DO i = 1, this%innx
1640 point = georef_coord_new(x=out%dim%lon(i,j), y=out%dim%lat(i,j))
1642 DO n = nprev, this%trans%poly%arraysize
1643 IF (
inside(point, this%trans%poly%array(n)))
THEN
1644 this%inter_index_x(i,j) = n
1649 DO n = nprev-1, 1, -1
1650 IF (
inside(point, this%trans%poly%array(n)))
THEN
1651 this%inter_index_x(i,j) = n
1664ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN
1667 CALL get_val(in, nx=this%innx, ny=this%inny)
1668 this%outnx = this%innx
1669 this%outny = this%inny
1671 IF (this%trans%sub_type ==
'maskvalid' .OR. this%trans%sub_type ==
'maskinvalid')
THEN
1673 IF (.NOT.
PRESENT(maskgrid))
THEN
1674 CALL l4f_category_log(this%category,l4f_error, &
1675 'grid_transform_init maskgrid argument missing for metamorphosis:'// &
1676 trim(this%trans%sub_type)//
' transformation')
1681 IF (this%innx /=
SIZE(maskgrid,1) .OR. this%inny /=
SIZE(maskgrid,2))
THEN
1682 CALL l4f_category_log(this%category,l4f_error, &
1683 'grid_transform_init mask non conformal with input field')
1684 CALL l4f_category_log(this%category,l4f_error, &
1685 'mask: '//t2c(
SIZE(maskgrid,1))//
'x'//t2c(
SIZE(maskgrid,2))// &
1686 ' input field:'//t2c(this%innx)//
'x'//t2c(this%inny))
1691 ALLOCATE(this%point_mask(this%innx,this%inny))
1693 IF (this%trans%sub_type ==
'maskvalid')
THEN
1696 IF (.NOT.
PRESENT(maskbounds))
THEN
1697 this%point_mask(:,:) =
c_e(maskgrid(:,:))
1698 ELSE IF (
SIZE(maskbounds) < 2)
THEN
1699 this%point_mask(:,:) =
c_e(maskgrid(:,:))
1701 this%point_mask(:,:) =
c_e(maskgrid(:,:)) .AND. &
1702 maskgrid(:,:) > maskbounds(1) .AND. &
1703 maskgrid(:,:) <= maskbounds(
SIZE(maskbounds))
1706 this%point_mask(:,:) = .NOT.
c_e(maskgrid(:,:))
1711 ELSE IF (this%trans%sub_type ==
'lemaskinvalid' .OR. &
1712 this%trans%sub_type ==
'ltmaskinvalid' .OR. &
1713 this%trans%sub_type ==
'gemaskinvalid' .OR. &
1714 this%trans%sub_type ==
'gtmaskinvalid')
THEN
1717 this%val_mask = maskgrid
1720 ELSE IF (this%trans%sub_type ==
'setinvalidto')
THEN
1722 IF (.NOT.
PRESENT(maskbounds))
THEN
1723 CALL l4f_category_log(this%category,l4f_error, &
1724 'grid_transform_init maskbounds missing for metamorphosis:'// &
1725 trim(this%trans%sub_type)//
' transformation')
1727 ELSE IF (
SIZE(maskbounds) < 1)
THEN
1728 CALL l4f_category_log(this%category,l4f_error, &
1729 'grid_transform_init maskbounds empty for metamorphosis:'// &
1730 trim(this%trans%sub_type)//
' transformation')
1733 this%val1 = maskbounds(1)
1735 CALL l4f_category_log(this%category, l4f_debug, &
1736 "grid_transform_init setting invalid data to "//t2c(this%val1))
1742 ELSE IF (this%trans%sub_type ==
'settoinvalid')
THEN
1744 IF (.NOT.
PRESENT(maskbounds))
THEN
1745 CALL l4f_category_log(this%category,l4f_error, &
1746 'grid_transform_init maskbounds missing for metamorphosis:'// &
1747 trim(this%trans%sub_type)//
' transformation')
1750 ELSE IF (
SIZE(maskbounds) < 2)
THEN
1751 CALL l4f_category_log(this%category,l4f_error, &
1752 'grid_transform_init maskbounds must have at least 2 elements for metamorphosis:'// &
1753 trim(this%trans%sub_type)//
' transformation')
1757 this%val1 = maskbounds(1)
1758 this%val2 = maskbounds(
SIZE(maskbounds))
1760 CALL l4f_category_log(this%category, l4f_debug, &
1761 "grid_transform_init setting to invalid interval ]"//t2c(this%val1)//
','// &
1762 t2c(this%val2)//
']')
1776SUBROUTINE outgrid_setup()
1779CALL griddim_setsteps(out)
1781CALL get_val(in,
proj=proj_in, component_flag=cf_i)
1782CALL get_val(out,
proj=proj_out, component_flag=cf_o)
1783IF (proj_in == proj_out)
THEN
1786 CALL set_val(out, component_flag=cf_i)
1791 CALL l4f_category_log(this%category,l4f_warn, &
1792 'trying to interpolate a grid with component flag 1 to a grid on a different projection')
1793 CALL l4f_category_log(this%category,l4f_warn, &
1794 'vector fields will probably be wrong')
1796 CALL set_val(out, component_flag=cf_i)
1800CALL griddim_set_central_lon(in, griddim_central_lon(out))
1802END SUBROUTINE outgrid_setup
1804END SUBROUTINE grid_transform_init
1849SUBROUTINE grid_transform_grid_vol7d_init(this, trans, in, v7d_out, &
1850 maskgrid, maskbounds, find_index, categoryappend)
1854TYPE(
vol7d),
INTENT(inout) :: v7d_out
1855REAL,
INTENT(in),
OPTIONAL :: maskgrid(:,:)
1856REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
1857PROCEDURE(basic_find_index),
POINTER,
OPTIONAL :: find_index
1858CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
1860INTEGER :: ix, iy, n, nm, nr, nprev, nmaskarea, xnmin, xnmax, ynmin, ynmax, &
1862DOUBLE PRECISION :: xmin, xmax, ymin, ymax, r2, lonref
1863DOUBLE PRECISION,
ALLOCATABLE :: lon1(:), lat1(:), lon(:,:), lat(:,:)
1864REAL,
ALLOCATABLE :: lmaskbounds(:)
1869IF (
PRESENT(find_index))
THEN
1870 IF (
ASSOCIATED(find_index))
THEN
1871 this%find_index => find_index
1874CALL grid_transform_init_common(this, trans, categoryappend)
1876CALL l4f_category_log(this%category, l4f_debug,
"grid_transform vg6d-v7d")
1880CALL get_val(trans, time_definition=time_definition)
1881IF (.NOT.
c_e(time_definition))
THEN
1885IF (this%trans%trans_type ==
'inter')
THEN
1887 IF (this%trans%sub_type ==
'near' .OR. this%trans%sub_type ==
'bilin' &
1888 .OR. this%trans%sub_type ==
'shapiro_near')
THEN
1891 CALL get_val(lin, nx=this%innx, ny=this%inny)
1892 this%outnx =
SIZE(v7d_out%ana)
1895 ALLOCATE (this%inter_index_x(this%outnx,this%outny),&
1896 this%inter_index_y(this%outnx,this%outny))
1897 ALLOCATE(lon(this%outnx,1),lat(this%outnx,1))
1899 CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
1901 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
1902 CALL griddim_set_central_lon(lin, lonref)
1903 CALL get_val(lin, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
1905 IF (this%trans%sub_type ==
'bilin')
THEN
1906 CALL this%find_index(lin, .false., &
1907 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1908 lon, lat, this%trans%extrap, &
1909 this%inter_index_x, this%inter_index_y)
1911 ALLOCATE(this%inter_x(this%innx,this%inny),this%inter_y(this%innx,this%inny))
1912 ALLOCATE(this%inter_xp(this%outnx,this%outny),this%inter_yp(this%outnx,this%outny))
1914 CALL griddim_gen_coord(lin, this%inter_x, this%inter_y)
1915 CALL proj(lin, lon, lat, this%inter_xp, this%inter_yp)
1918 CALL this%find_index(lin, .true., &
1919 this%innx, this%inny, xmin, xmax, ymin, ymax, &
1920 lon, lat, this%trans%extrap, &
1921 this%inter_index_x, this%inter_index_y)
1923 IF (this%trans%trans_type ==
'intersearch')
THEN
1924 ALLOCATE(this%inter_x(this%innx,this%inny),this%inter_y(this%innx,this%inny))
1925 ALLOCATE(this%inter_xp(this%outnx,this%outny),this%inter_yp(this%outnx,this%outny))
1927 CALL griddim_gen_coord(lin, this%inter_x, this%inter_y)
1928 CALL proj(lin, lon, lat, this%inter_xp, this%inter_yp)
1939ELSE IF (this%trans%trans_type ==
'polyinter')
THEN
1941 CALL get_val(in, nx=this%innx, ny=this%inny)
1943 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
1944 this%inter_index_y(this%innx,this%inny))
1945 this%inter_index_x(:,:) = imiss
1946 this%inter_index_y(:,:) = 1
1952 this%outnx = this%trans%poly%arraysize
1955 CALL init(v7d_out, time_definition=time_definition)
1956 CALL vol7d_alloc(v7d_out, nana=this%outnx)
1959 ALLOCATE(lon(this%outnx,1))
1960 CALL getval(v7d_out%ana(:)%coord,lon=lon(:,1))
1961 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
1962 CALL griddim_set_central_lon(lin, lonref)
1967 CALL poly_to_coordinates(this%trans%poly, v7d_out)
1972 DO iy = 1, this%inny
1973 inside_x:
DO ix = 1, this%innx
1974 point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy))
1976 DO n = nprev, this%trans%poly%arraysize
1977 IF (
inside(point, this%trans%poly%array(n)))
THEN
1978 this%inter_index_x(ix,iy) = n
1983 DO n = nprev-1, 1, -1
1984 IF (
inside(point, this%trans%poly%array(n)))
THEN
1985 this%inter_index_x(ix,iy) = n
2037 this%stencil(:,:) = .true.
2040 IF ((ix-nm)**2+(iy-nm)**2 > r2) this%stencil(ix,iy) = .false.
2047 xnmax = this%innx - nr
2049 ynmax = this%inny - nr
2050 DO iy = 1, this%outny
2051 DO ix = 1, this%outnx
2052 IF (this%inter_index_x(ix,iy) < xnmin .OR. &
2053 this%inter_index_x(ix,iy) > xnmax .OR. &
2054 this%inter_index_y(ix,iy) < ynmin .OR. &
2055 this%inter_index_y(ix,iy) > ynmax)
THEN
2056 this%inter_index_x(ix,iy) = imiss
2057 this%inter_index_y(ix,iy) = imiss
2063 CALL l4f_category_log(this%category, l4f_debug, &
2064 'stencilinter: stencil size '//t2c(n*n))
2065 CALL l4f_category_log(this%category, l4f_debug, &
2066 'stencilinter: stencil points '//t2c(count(this%stencil)))
2074ELSE IF (this%trans%trans_type ==
'maskinter')
THEN
2076 IF (.NOT.
PRESENT(maskgrid))
THEN
2077 CALL l4f_category_log(this%category,l4f_error, &
2078 'grid_transform_init maskgrid argument missing for maskinter transformation')
2079 CALL raise_fatal_error()
2082 CALL get_val(in, nx=this%innx, ny=this%inny)
2083 IF (this%innx /=
SIZE(maskgrid,1) .OR. this%inny /=
SIZE(maskgrid,2))
THEN
2084 CALL l4f_category_log(this%category,l4f_error, &
2085 'grid_transform_init mask non conformal with input field')
2086 CALL l4f_category_log(this%category,l4f_error, &
2087 'mask: '//t2c(
SIZE(maskgrid,1))//
'x'//t2c(
SIZE(maskgrid,2))// &
2088 ' input field:'//t2c(this%innx)//
'x'//t2c(this%inny))
2089 CALL raise_fatal_error()
2092 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2093 this%inter_index_y(this%innx,this%inny))
2094 this%inter_index_x(:,:) = imiss
2095 this%inter_index_y(:,:) = 1
2098 CALL gen_mask_class()
2106 DO iy = 1, this%inny
2107 DO ix = 1, this%innx
2108 IF (
c_e(maskgrid(ix,iy)))
THEN
2109 IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1))
THEN
2110 DO n = nmaskarea, 1, -1
2111 IF (maskgrid(ix,iy) > lmaskbounds(n))
THEN
2112 this%inter_index_x(ix,iy) = n
2122 this%outnx = nmaskarea
2125 CALL init(v7d_out, time_definition=time_definition)
2126 CALL vol7d_alloc(v7d_out, nana=nmaskarea)
2131 CALL init(v7d_out%ana(n), &
2133 mask=(this%inter_index_x(:,:) == n))), &
2135 mask=(this%inter_index_x(:,:) == n))))
2141ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN
2148 CALL get_val(in, nx=this%innx, ny=this%inny)
2150 ALLOCATE(this%point_index(this%innx,this%inny))
2151 this%point_index(:,:) = imiss
2154 CALL init(v7d_out, time_definition=time_definition)
2156 IF (this%trans%sub_type ==
'all' )
THEN
2158 this%outnx = this%innx*this%inny
2160 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2165 CALL init(v7d_out%ana((iy-1)*this%innx+ix), &
2166 lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
2168 this%point_index(ix,iy) = n
2174 ELSE IF (this%trans%sub_type ==
'coordbb' )
THEN
2179 DO iy = 1, this%inny
2180 DO ix = 1, this%innx
2182 IF (lin%dim%lon(ix,iy) > this%trans%rect_coo%ilon .AND. &
2183 lin%dim%lon(ix,iy) < this%trans%rect_coo%flon .AND. &
2184 lin%dim%lat(ix,iy) > this%trans%rect_coo%ilat .AND. &
2185 lin%dim%lat(ix,iy) < this%trans%rect_coo%flat)
THEN
2186 this%outnx = this%outnx + 1
2187 this%point_index(ix,iy) = this%outnx
2192 IF (this%outnx <= 0)
THEN
2193 CALL l4f_category_log(this%category,l4f_warn, &
2194 "metamorphosis:coordbb: no points inside bounding box "//&
2195 trim(
to_char(this%trans%rect_coo%ilon))//
","// &
2196 trim(
to_char(this%trans%rect_coo%flon))//
","// &
2197 trim(
to_char(this%trans%rect_coo%ilat))//
","// &
2198 trim(
to_char(this%trans%rect_coo%flat)))
2201 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2205 DO iy = 1, this%inny
2206 DO ix = 1, this%innx
2207 IF (
c_e(this%point_index(ix,iy)))
THEN
2209 CALL init(v7d_out%ana(n), &
2210 lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
2217 ELSE IF (this%trans%sub_type ==
'poly' )
THEN
2226 DO iy = 1, this%inny
2227 DO ix = 1, this%innx
2228 point = georef_coord_new(x=lin%dim%lon(ix,iy), y=lin%dim%lat(ix,iy))
2229 DO n = 1, this%trans%poly%arraysize
2230 IF (
inside(point, this%trans%poly%array(n)))
THEN
2231 this%outnx = this%outnx + 1
2232 this%point_index(ix,iy) = n
2241 IF (this%outnx <= 0)
THEN
2242 CALL l4f_category_log(this%category,l4f_warn, &
2243 "metamorphosis:poly: no points inside polygons")
2246 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2249 DO iy = 1, this%inny
2250 DO ix = 1, this%innx
2251 IF (
c_e(this%point_index(ix,iy)))
THEN
2253 CALL init(v7d_out%ana(n), &
2254 lon=lin%dim%lon(ix,iy), lat=lin%dim%lat(ix,iy))
2261 ELSE IF (this%trans%sub_type ==
'mask' )
THEN
2263 IF (.NOT.
PRESENT(maskgrid))
THEN
2264 CALL l4f_category_log(this%category,l4f_error, &
2265 'grid_transform_init maskgrid argument missing for metamorphosis:mask transformation')
2270 IF (this%innx /=
SIZE(maskgrid,1) .OR. this%inny /=
SIZE(maskgrid,2))
THEN
2271 CALL l4f_category_log(this%category,l4f_error, &
2272 'grid_transform_init mask non conformal with input field')
2273 CALL l4f_category_log(this%category,l4f_error, &
2274 'mask: '//t2c(
SIZE(maskgrid,1))//
'x'//t2c(
SIZE(maskgrid,2))// &
2275 ' input field:'//t2c(this%innx)//
'x'//t2c(this%inny))
2281 CALL gen_mask_class()
2289 DO iy = 1, this%inny
2290 DO ix = 1, this%innx
2291 IF (
c_e(maskgrid(ix,iy)))
THEN
2292 IF (maskgrid(ix,iy) <= lmaskbounds(nmaskarea+1))
THEN
2293 DO n = nmaskarea, 1, -1
2294 IF (maskgrid(ix,iy) > lmaskbounds(n))
THEN
2295 this%outnx = this%outnx + 1
2296 this%point_index(ix,iy) = n
2306 IF (this%outnx <= 0)
THEN
2307 CALL l4f_category_log(this%category,l4f_warn, &
2308 "grid_transform_init no points inside mask for metamorphosis:mask transformation")
2312 CALL l4f_category_log(this%category,l4f_info, &
2313 "points in subarea "//t2c(n)//
": "// &
2314 t2c(count(this%point_index(:,:) == n)))
2317 CALL vol7d_alloc(v7d_out, nana=this%outnx)
2320 DO iy = 1, this%inny
2321 DO ix = 1, this%innx
2322 IF (
c_e(this%point_index(ix,iy)))
THEN
2324 CALL init(v7d_out%ana(n),lon=lin%dim%lon(ix,iy),lat=lin%dim%lat(ix,iy))
2337SUBROUTINE gen_mask_class()
2338REAL :: startmaskclass, mmin, mmax
2341IF (
PRESENT(maskbounds))
THEN
2342 nmaskarea =
SIZE(maskbounds) - 1
2343 IF (nmaskarea > 0)
THEN
2344 lmaskbounds = maskbounds
2346 ELSE IF (nmaskarea == 0)
THEN
2347 CALL l4f_category_log(this%category,l4f_warn, &
2348 'only one value provided for maskbounds, '//t2c(maskbounds(1)) &
2349 //
', it will be ignored')
2350 CALL l4f_category_log(this%category,l4f_warn, &
2351 'at least 2 values are required for maskbounds')
2355mmin = minval(maskgrid, mask=
c_e(maskgrid))
2356mmax = maxval(maskgrid, mask=
c_e(maskgrid))
2358nmaskarea = int(mmax - mmin + 1.5)
2359startmaskclass = mmin - 0.5
2361ALLOCATE(lmaskbounds(nmaskarea+1))
2362lmaskbounds(:) = (/(startmaskclass+real(i),i=0,nmaskarea)/)
2364CALL l4f_category_log(this%category,l4f_debug, &
2365 'maskinter, '//t2c(nmaskarea)//
' subareas defined, with boundaries:')
2366DO i = 1,
SIZE(lmaskbounds)
2367 CALL l4f_category_log(this%category,l4f_debug, &
2368 'maskinter '//t2c(i)//
' '//t2c(lmaskbounds(i)))
2372END SUBROUTINE gen_mask_class
2374END SUBROUTINE grid_transform_grid_vol7d_init
2395SUBROUTINE grid_transform_vol7d_grid_init(this, trans, v7d_in, out, categoryappend)
2398TYPE(
vol7d),
INTENT(in) :: v7d_in
2400character(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2403DOUBLE PRECISION :: xmin, xmax, ymin, ymax, lonref
2404DOUBLE PRECISION,
ALLOCATABLE :: lon(:,:),lat(:,:)
2408CALL grid_transform_init_common(this, trans, categoryappend)
2410CALL l4f_category_log(this%category, l4f_debug,
"grid_transform v7d-vg6d")
2413IF (this%trans%trans_type ==
'inter')
THEN
2415 IF ( this%trans%sub_type ==
'linear' )
THEN
2417 this%innx=
SIZE(v7d_in%ana)
2419 ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2420 ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
2421 CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2423 CALL copy (out, lout)
2425 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
2426 CALL griddim_set_central_lon(lout, lonref)
2428 CALL get_val(lout, nx=nx, ny=ny)
2431 ALLOCATE(this%inter_x(this%outnx,this%outny),this%inter_y(this%outnx,this%outny))
2433 CALL get_val(lout, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2434 CALL proj(lout, lon, lat, this%inter_xp, this%inter_yp)
2435 CALL griddim_gen_coord(lout, this%inter_x, this%inter_y)
2444ELSE IF (this%trans%trans_type ==
'boxinter')
THEN
2446 this%innx=
SIZE(v7d_in%ana)
2449 ALLOCATE(lon(this%innx,1),lat(this%innx,1))
2450 ALLOCATE(this%inter_index_x(this%innx,this%inny), &
2451 this%inter_index_y(this%innx,this%inny))
2453 CALL getval(v7d_in%ana(:)%coord,lon=lon(:,1),lat=lat(:,1))
2455 CALL copy (out, lout)
2457 lonref = 0.5d0*(maxval(lon(:,1), mask=
c_e(lon(:,1))) + minval(lon(:,1), mask=
c_e(lon(:,1))))
2458 CALL griddim_set_central_lon(lout, lonref)
2460 CALL get_val(lout, nx=this%outnx, ny=this%outny, &
2461 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
2464 IF (.NOT.
c_e(this%trans%area_info%boxdx)) &
2465 CALL get_val(out, dx=this%trans%area_info%boxdx)
2466 IF (.NOT.
c_e(this%trans%area_info%boxdy)) &
2467 CALL get_val(out, dx=this%trans%area_info%boxdy)
2469 this%trans%area_info%boxdx = this%trans%area_info%boxdx*0.5d0
2470 this%trans%area_info%boxdy = this%trans%area_info%boxdy*0.5d0
2473 CALL this%find_index(lout, .true., &
2474 this%outnx, this%outny, xmin, xmax, ymin, ymax, &
2475 lon, lat, .false., &
2476 this%inter_index_x, this%inter_index_y)
2485END SUBROUTINE grid_transform_vol7d_grid_init
2522SUBROUTINE grid_transform_vol7d_vol7d_init(this, trans, v7d_in, v7d_out, &
2523 maskbounds, categoryappend)
2526TYPE(
vol7d),
INTENT(in) :: v7d_in
2527TYPE(
vol7d),
INTENT(inout) :: v7d_out
2528REAL,
INTENT(in),
OPTIONAL :: maskbounds(:)
2529CHARACTER(len=*),
INTENT(in),
OPTIONAL :: categoryappend
2532DOUBLE PRECISION,
ALLOCATABLE :: lon(:), lat(:)
2534DOUBLE PRECISION :: lon1, lat1
2538CALL grid_transform_init_common(this, trans, categoryappend)
2540CALL l4f_category_log(this%category, l4f_debug,
"grid_transform v7d-v7d")
2543IF (this%trans%trans_type ==
'inter')
THEN
2545 IF ( this%trans%sub_type ==
'linear' )
THEN
2547 CALL vol7d_alloc_vol(v7d_out)
2548 this%outnx=
SIZE(v7d_out%ana)
2551 this%innx=
SIZE(v7d_in%ana)
2554 ALLOCATE(this%inter_xp(this%innx,this%inny),this%inter_yp(this%innx,this%inny))
3099 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3100 this%trans%sub_type ==
'stddevnm1')
THEN
3102 IF (this%trans%sub_type ==
'stddev')
THEN
3108 navg = this%trans%box_info%npx*this%trans%box_info%npy
3111 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3112 je = j+this%trans%box_info%npy-1
3115 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3116 ie = i+this%trans%box_info%npx-1
3119 reshape(field_in(i:ie,j:je,k),(/navg/)), nm1=nm1)
3124 ELSE IF (this%trans%sub_type ==
'max')
THEN
3127 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3128 je = j+this%trans%box_info%npy-1
3131 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3132 ie = i+this%trans%box_info%npx-1
3134 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3136 field_out(ii,jj,k) = maxval(field_in(i:ie,j:je,k), &
3137 mask=(field_in(i:ie,j:je,k) /= rmiss))
3143 ELSE IF (this%trans%sub_type ==
'min')
THEN
3146 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3147 je = j+this%trans%box_info%npy-1
3150 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3151 ie = i+this%trans%box_info%npx-1
3153 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3155 field_out(ii,jj,k) = minval(field_in(i:ie,j:je,k), &
3156 mask=(field_in(i:ie,j:je,k) /= rmiss))
3162 ELSE IF (this%trans%sub_type ==
'percentile')
THEN
3164 navg = this%trans%box_info%npx*this%trans%box_info%npy
3167 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3168 je = j+this%trans%box_info%npy-1
3171 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3172 ie = i+this%trans%box_info%npx-1
3175 reshape(field_in(i:ie,j:je,k),(/navg/)), &
3176 (/real(this%trans%stat_info%percentile)/))
3181 ELSE IF (this%trans%sub_type ==
'frequency')
THEN
3185 DO j = 1, this%inny - this%trans%box_info%npy + 1, this%trans%box_info%npy
3186 je = j+this%trans%box_info%npy-1
3189 DO i = 1, this%innx - this%trans%box_info%npx + 1, this%trans%box_info%npx
3190 ie = i+this%trans%box_info%npx-1
3192 navg = count(field_in(i:ie,j:je,k) /= rmiss)
3194 field_out(ii,jj,k) = sum(interval_info_valid( &
3195 this%trans%interval_info, field_in(i:ie,j:je,k)), &
3196 mask=(field_in(i:ie,j:je,k) /= rmiss))/navg
3204ELSE IF (this%trans%trans_type ==
'inter')
THEN
3206 IF (this%trans%sub_type ==
'near')
THEN
3209 DO j = 1, this%outny
3210 DO i = 1, this%outnx
3212 IF (
c_e(this%inter_index_x(i,j))) field_out(i,j,k) = &
3213 field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3219 ELSE IF (this%trans%sub_type ==
'bilin')
THEN
3222 DO j = 1, this%outny
3223 DO i = 1, this%outnx
3225 IF (
c_e(this%inter_index_x(i,j)))
THEN
3227 z1=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k)
3228 z2=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j),k)
3229 z3=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1,k)
3230 z4=field_in(this%inter_index_x(i,j),this%inter_index_y(i,j)+1,k)
3232 IF (
c_e(z1) .AND.
c_e(z2) .AND.
c_e(z3) .AND.
c_e(z4))
THEN
3234 x1=this%inter_x(this%inter_index_x(i,j),this%inter_index_y(i,j))
3235 y1=this%inter_y(this%inter_index_x(i,j),this%inter_index_y(i,j))
3236 x3=this%inter_x(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3237 y3=this%inter_y(this%inter_index_x(i,j)+1,this%inter_index_y(i,j)+1)
3239 xp=this%inter_xp(i,j)
3240 yp=this%inter_yp(i,j)
3242 field_out(i,j,k) = hbilin(z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
3250 ELSE IF (this%trans%sub_type ==
'shapiro_near')
THEN
3252 DO j = 1, this%outny
3253 DO i = 1, this%outnx
3255 IF (
c_e(this%inter_index_x(i,j)))
THEN
3257 IF(this%inter_index_x(i,j)-1>0)
THEN
3258 z(1) = field_in(this%inter_index_x(i,j)-1,this%inter_index_y(i,j) ,k)
3262 IF(this%inter_index_x(i,j)+1<=this%outnx)
THEN
3263 z(3)=field_in(this%inter_index_x(i,j)+1,this%inter_index_y(i,j) ,k)
3267 IF(this%inter_index_y(i,j)+1<=this%outny)
THEN
3268 z(2)=field_in(this%inter_index_x(i,j) ,this%inter_index_y(i,j)+1,k)
3272 IF(this%inter_index_y(i,j)-1>0)
THEN
3273 z(4)=field_in(this%inter_index_x(i,j), this%inter_index_y(i,j)-1,k)
3277 field_out(i,j,k) = shapiro(z,field_in(this%inter_index_x(i,j),this%inter_index_y(i,j),k))
3286ELSE IF (this%trans%trans_type ==
'intersearch')
THEN
3289 likethis%trans%trans_type =
'inter'
3290 CALL grid_transform_compute(likethis, field_in, field_out, var, coord_3d_in)
3293 IF ((.NOT.all(
c_e(field_out(:,:,k)))) .AND. (any(
c_e(field_in(:,:,k)))))
THEN
3294 DO j = 1, this%outny
3295 DO i = 1, this%outnx
3296 IF (.NOT.
c_e(field_out(i,j,k)))
THEN
3300 IF (
c_e(field_in(l,m,k)))
THEN
3301 disttmp = (this%inter_xp(l,m) - this%inter_x(i,j))**2 + (this%inter_yp(l,m) - this%inter_y(i,j))**2
3302 IF (disttmp <
dist)
THEN
3304 field_out(i,j,k) = field_in(l,m,k)
3315ELSE IF (this%trans%trans_type ==
'boxinter' &
3316 .OR. this%trans%trans_type ==
'polyinter' &
3317 .OR. this%trans%trans_type ==
'maskinter')
THEN
3319 IF (this%trans%sub_type ==
'average')
THEN
3321 IF (vartype == var_dir360)
THEN
3323 DO j = 1, this%outny
3324 DO i = 1, this%outnx
3325 field_out(i,j,k) = find_prevailing_direction(field_in(:,:,k), &
3327 mask=this%inter_index_x(:,:) == i .AND. this%inter_index_y(:,:) == j)
3333 ALLOCATE(nval(this%outnx, this%outny))
3334 field_out(:,:,:) = 0.0
3339 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3340 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3341 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3343 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3344 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3348 WHERE (nval(:,:) /= 0)
3349 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3351 field_out(:,:,k) = rmiss
3357 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3358 this%trans%sub_type ==
'stddevnm1')
THEN
3360 IF (this%trans%sub_type ==
'stddev')
THEN
3366 DO j = 1, this%outny
3367 DO i = 1, this%outnx
3370 reshape(field_in(:,:,k), (/
SIZE(field_in(:,:,k))/)), &
3371 mask=reshape((this%inter_index_x == i .AND. &
3372 this%inter_index_y == j), (/
SIZE(field_in(:,:,k))/)), nm1=nm1)
3377 ELSE IF (this%trans%sub_type ==
'max')
THEN
3382 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3383 IF (
c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k)))
THEN
3384 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3385 max(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3388 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3397 ELSE IF (this%trans%sub_type ==
'min')
THEN
3402 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3403 IF (
c_e(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k)))
THEN
3404 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3405 min(field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k), &
3408 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3416 ELSE IF (this%trans%sub_type ==
'percentile')
THEN
3419 DO j = 1, this%outny
3420 DO i = 1, this%outnx
3423 reshape(field_in(:,:,k), (/
SIZE(field_in(:,:,k))/)), &
3424 (/real(this%trans%stat_info%percentile)/), &
3425 mask=reshape((this%inter_index_x == i .AND. &
3426 this%inter_index_y == j), (/
SIZE(field_in(:,:,k))/)))
3431 ELSE IF (this%trans%sub_type ==
'frequency')
THEN
3433 ALLOCATE(nval(this%outnx, this%outny))
3434 field_out(:,:,:) = 0.0
3439 IF (
c_e(this%inter_index_x(i,j)) .AND.
c_e(field_in(i,j,k)))
THEN
3440 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) = &
3441 field_out(this%inter_index_x(i,j),this%inter_index_y(i,j),k) + &
3442 interval_info_valid(this%trans%interval_info, field_in(i,j,k))
3443 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) = &
3444 nval(this%inter_index_x(i,j),this%inter_index_y(i,j)) + 1
3448 WHERE (nval(:,:) /= 0)
3449 field_out(:,:,k) = field_out(:,:,k)/nval(:,:)
3451 field_out(:,:,k) = rmiss
3458ELSE IF (this%trans%trans_type ==
'stencilinter')
THEN
3459 np =
SIZE(this%stencil,1)/2
3460 ns =
SIZE(this%stencil)
3462 IF (this%trans%sub_type ==
'average')
THEN
3464 IF (vartype == var_dir360)
THEN
3466 DO j = 1, this%outny
3467 DO i = 1, this%outnx
3468 IF (
c_e(this%inter_index_x(i,j)))
THEN
3469 i1 = this%inter_index_x(i,j) - np
3470 i2 = this%inter_index_x(i,j) + np
3471 j1 = this%inter_index_y(i,j) - np
3472 j2 = this%inter_index_y(i,j) + np
3473 field_out(i,j,k) = find_prevailing_direction(field_in(i1:i2,j1:j2,k), &
3475 mask=this%stencil(:,:))
3486 DO j = 1, this%outny
3487 DO i = 1, this%outnx
3488 IF (
c_e(this%inter_index_x(i,j)))
THEN
3489 i1 = this%inter_index_x(i,j) - np
3490 i2 = this%inter_index_x(i,j) + np
3491 j1 = this%inter_index_y(i,j) - np
3492 j2 = this%inter_index_y(i,j) + np
3493 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3495 field_out(i,j,k) = sum(field_in(i1:i2,j1:j2,k), &
3496 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3505 ELSE IF (this%trans%sub_type ==
'stddev' .OR. &
3506 this%trans%sub_type ==
'stddevnm1')
THEN
3508 IF (this%trans%sub_type ==
'stddev')
THEN
3517 DO j = 1, this%outny
3518 DO i = 1, this%outnx
3519 IF (
c_e(this%inter_index_x(i,j)))
THEN
3520 i1 = this%inter_index_x(i,j) - np
3521 i2 = this%inter_index_x(i,j) + np
3522 j1 = this%inter_index_y(i,j) - np
3523 j2 = this%inter_index_y(i,j) + np
3526 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3527 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3528 this%stencil(:,:), (/ns/)), nm1=nm1)
3535 ELSE IF (this%trans%sub_type ==
'max')
THEN
3540 DO j = 1, this%outny
3541 DO i = 1, this%outnx
3542 IF (
c_e(this%inter_index_x(i,j)))
THEN
3543 i1 = this%inter_index_x(i,j) - np
3544 i2 = this%inter_index_x(i,j) + np
3545 j1 = this%inter_index_y(i,j) - np
3546 j2 = this%inter_index_y(i,j) + np
3547 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3549 field_out(i,j,k) = maxval(field_in(i1:i2,j1:j2,k), &
3550 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3558 ELSE IF (this%trans%sub_type ==
'min')
THEN
3563 DO j = 1, this%outny
3564 DO i = 1, this%outnx
3565 IF (
c_e(this%inter_index_x(i,j)))
THEN
3566 i1 = this%inter_index_x(i,j) - np
3567 i2 = this%inter_index_x(i,j) + np
3568 j1 = this%inter_index_y(i,j) - np
3569 j2 = this%inter_index_y(i,j) + np
3570 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3572 field_out(i,j,k) = minval(field_in(i1:i2,j1:j2,k), &
3573 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3581 ELSE IF (this%trans%sub_type ==
'percentile')
THEN
3586 DO j = 1, this%outny
3587 DO i = 1, this%outnx
3588 IF (
c_e(this%inter_index_x(i,j)))
THEN
3589 i1 = this%inter_index_x(i,j) - np
3590 i2 = this%inter_index_x(i,j) + np
3591 j1 = this%inter_index_y(i,j) - np
3592 j2 = this%inter_index_y(i,j) + np
3595 reshape(field_in(i1:i2,j1:j2,k), (/ns/)), &
3596 (/real(this%trans%stat_info%percentile)/), &
3597 mask=reshape(field_in(i1:i2,j1:j2,k) /= rmiss .AND. &
3598 this%stencil(:,:), (/ns/)))
3605 ELSE IF (this%trans%sub_type ==
'frequency')
THEN
3610 DO j = 1, this%outny
3611 DO i = 1, this%outnx
3612 IF (
c_e(this%inter_index_x(i,j)))
THEN
3613 i1 = this%inter_index_x(i,j) - np
3614 i2 = this%inter_index_x(i,j) + np
3615 j1 = this%inter_index_y(i,j) - np
3616 j2 = this%inter_index_y(i,j) + np
3617 n = count(field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))
3619 field_out(i,j,k) = sum(interval_info_valid( &
3620 this%trans%interval_info, field_in(i1:i2,j1:j2,k)), &
3621 mask=field_in(i1:i2,j1:j2,k) /= rmiss .AND. this%stencil(:,:))/n
3631ELSE IF (this%trans%trans_type ==
'maskgen')
THEN
3634 WHERE(
c_e(this%inter_index_x(:,:)))
3635 field_out(:,:,k) = real(this%inter_index_x(:,:))
3639ELSE IF (this%trans%trans_type ==
'metamorphosis')
THEN
3641 IF (this%trans%sub_type ==
'all')
THEN
3643 field_out(:,:,:) = reshape(field_in(:,:,:), (/this%outnx,this%outny,innz/))
3645 ELSE IF (this%trans%sub_type ==
'coordbb' .OR. this%trans%sub_type ==
'poly' &
3646 .OR. this%trans%sub_type ==
'mask')
THEN
3650 field_out(:,1,k) = pack(field_in(:,:,k),
c_e(this%point_index(:,:)))
3653 ELSE IF (this%trans%sub_type ==
'maskvalid' .OR. &
3654 this%trans%sub_type ==
'maskinvalid')
THEN
3657 WHERE (this%point_mask(:,:))
3658 field_out(:,:,k) = field_in(:,:,k)
3662 ELSE IF (this%trans%sub_type ==
'lemaskinvalid')
THEN
3665 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) > this%val_mask(:,:))
3666 field_out(:,:,k) = field_in(:,:,k)
3668 field_out(:,:,k) = rmiss
3672 ELSE IF (this%trans%sub_type ==
'ltmaskinvalid')
THEN
3675 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) >= this%val_mask(:,:))
3676 field_out(:,:,k) = field_in(:,:,k)
3678 field_out(:,:,k) = rmiss
3682 ELSE IF (this%trans%sub_type ==
'gemaskinvalid')
THEN
3685 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) < this%val_mask(:,:))
3686 field_out(:,:,k) = field_in(:,:,k)
3688 field_out(:,:,k) = rmiss
3692 ELSE IF (this%trans%sub_type ==
'gtmaskinvalid')
THEN
3695 WHERE (
c_e(field_in(:,:,k)) .AND. field_in(:,:,k) <= this%val_mask(:,:))
3696 field_out(:,:,k) = field_in(:,:,k)
3698 field_out(:,:,k) = rmiss
3702 ELSE IF (this%trans%sub_type ==
'setinvalidto')
THEN
3705 WHERE (
c_e(field_in(:,:,k)))
3706 field_out(:,:,k) = field_in(:,:,k)
3708 field_out(:,:,k) = this%val1
3712 ELSE IF (this%trans%sub_type ==
'settoinvalid')
THEN
3713 IF (
c_e(this%val1) .AND.
c_e(this%val2))
THEN
3714 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1 &
3715 .AND. field_in(:,:,:) <= this%val2)
3716 field_out(:,:,:) = rmiss
3718 field_out(:,:,:) = field_in(:,:,:)
3720 ELSE IF (
c_e(this%val1))
THEN
3721 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) > this%val1)
3722 field_out(:,:,:) = rmiss
3724 field_out(:,:,:) = field_in(:,:,:)
3726 ELSE IF (
c_e(this%val2))
THEN
3727 WHERE (
c_e(field_in(:,:,:)) .AND. field_in(:,:,:) <= this%val2)
3728 field_out(:,:,:) = rmiss
3730 field_out(:,:,:) = field_in(:,:,:)
3735ELSE IF (this%trans%trans_type ==
'vertint')
THEN
3737 IF (this%trans%sub_type ==
'linear')
THEN
3739 alloc_coord_3d_in_act = .false.
3740 IF (
ASSOCIATED(this%inter_index_z))
THEN
3743 IF (
c_e(this%inter_index_z(k)))
THEN
3744 z1 = real(this%inter_zp(k))
3745 z2 = real(1.0d0 - this%inter_zp(k))
3746 SELECT CASE(vartype)
3751 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3752 c_e(field_in(i,j,this%inter_index_z(k)+1)))
THEN
3753 field_out(i,j,k) = &
3754 interp_var_360(field_in(i,j,this%inter_index_z(k)), &
3755 field_in(i,j,this%inter_index_z(k)+1), z1, z2)
3763 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3764 c_e(field_in(i,j,this%inter_index_z(k)+1)) .AND. &
3765 field_in(i,j,this%inter_index_z(k)) > 0. .AND. &
3766 field_in(i,j,this%inter_index_z(k)+1) > 0.)
THEN
3767 field_out(i,j,k) = exp( &
3768 log(field_in(i,j,this%inter_index_z(k)))*z2 + &
3769 log(field_in(i,j,this%inter_index_z(k)+1))*z1)
3777 IF (
c_e(field_in(i,j,this%inter_index_z(k))) .AND. &
3778 c_e(field_in(i,j,this%inter_index_z(k)+1)))
THEN
3779 field_out(i,j,k) = &
3780 field_in(i,j,this%inter_index_z(k))*z2 + &
3781 field_in(i,j,this%inter_index_z(k)+1)*z1
3793 IF (
ALLOCATED(this%coord_3d_in))
THEN
3794 coord_3d_in_act => this%coord_3d_in
3796 CALL l4f_category_log(this%category,l4f_debug, &
3797 "Using external vertical coordinate for vertint")
3800 IF (
PRESENT(coord_3d_in))
THEN
3802 CALL l4f_category_log(this%category,l4f_debug, &
3803 "Using internal vertical coordinate for vertint")
3805 IF (this%dolog)
THEN
3806 ALLOCATE(coord_3d_in_act(
SIZE(coord_3d_in,1), &
3807 SIZE(coord_3d_in,2),
SIZE(coord_3d_in,3)))
3808 alloc_coord_3d_in_act = .true.
3809 WHERE (
c_e(coord_3d_in) .AND. coord_3d_in > 0.0)
3810 coord_3d_in_act = log(coord_3d_in)
3812 coord_3d_in_act = rmiss
3815 coord_3d_in_act => coord_3d_in
3818 CALL l4f_category_log(this%category,l4f_error, &
3819 "Missing internal and external vertical coordinate for vertint")
3825 inused =
SIZE(coord_3d_in_act,3)
3826 IF (inused < 2)
RETURN
3830 IF (.NOT.
c_e(this%vcoord_out(k))) cycle
3834 DO kk = 1, max(inused-kkcache-1, kkcache)
3836 kkdown = kkcache - kk + 1
3838 IF (kkdown >= 1)
THEN
3839 IF (this%vcoord_out(k) >= &
3840 min(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)) .AND. &
3841 this%vcoord_out(k) < &
3842 max(coord_3d_in_act(i,j,kkdown), coord_3d_in_act(i,j,kkdown+1)))
THEN
3845 kfound = kkcache + this%levshift
3849 IF (kkup < inused)
THEN
3850 IF (this%vcoord_out(k) >= &
3851 min(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)) .AND. &
3852 this%vcoord_out(k) < &
3853 max(coord_3d_in_act(i,j,kkup), coord_3d_in_act(i,j,kkup+1)))
THEN
3856 kfound = kkcache + this%levshift
3863 IF (
c_e(kfound))
THEN
3864 IF (
c_e(field_in(i,j,kfound)) .AND.
c_e(field_in(i,j,kfound+1)))
THEN
3865 z1 = real((this%vcoord_out(k) - coord_3d_in_act(i,j,kfoundin))/ &
3866 (coord_3d_in_act(i,j,kfoundin+1) - coord_3d_in_act(i,j,kfoundin)))
3868 SELECT CASE(vartype)
3871 field_out(i,j,k) = &
3872 interp_var_360(field_in(i,j,kfound), field_in(i,j,kfound+1), z1, z2)
3874 field_out(i,j,k) = exp(log(field_in(i,j,kfound))*z2 + &
3875 log(field_in(i,j,kfound+1))*z1)
3878 field_out(i,j,k) = field_in(i,j,kfound)*z2 + field_in(i,j,kfound+1)*z1
3887 IF (alloc_coord_3d_in_act)
DEALLOCATE(coord_3d_in_act)
3890 ELSE IF (this%trans%sub_type ==
'linearsparse')
THEN
3893 IF (.NOT.
ASSOCIATED(this%vcoord_in) .AND. .NOT.
PRESENT(coord_3d_in))
THEN
3894 CALL l4f_category_log(this%category,l4f_error, &
3895 "linearsparse interpolation, no input vert coord available")
3899 ALLOCATE(coord_in(innz),val_in(innz),mask_in(innz))
3903 IF (
ASSOCIATED(this%vcoord_in))
THEN
3904 mask_in =
c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3905 .AND.
c_e(this%vcoord_in(:))
3907 mask_in =
c_e(field_in(i,j,this%levshift+1:this%levshift+this%levused)) &
3908 .AND.
c_e(coord_3d_in(i,j,:))
3911 IF (vartype == var_press)
THEN
3912 mask_in(:) = mask_in(:) .AND. &
3913 (field_in(i,j,this%levshift+1:this%levshift+this%levused) > 0.0d0)
3915 inused = count(mask_in)
3916 IF (inused > 1)
THEN
3917 IF (
ASSOCIATED(this%vcoord_in))
THEN
3918 coord_in(1:inused) = pack(this%vcoord_in(:), mask=mask_in)
3920 coord_in(1:inused) = pack(coord_3d_in(i,j,:), mask=mask_in)
3922 IF (vartype == var_press)
THEN
3923 val_in(1:inused) = log(pack( &
3924 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3927 val_in(1:inused) = pack( &
3928 field_in(i,j,this%levshift+1:this%levshift+this%levused), &
3935 DO kk = 1, max(inused-kkcache-1, kkcache)
3937 kkdown = kkcache - kk + 1
3939 IF (kkdown >= 1)
THEN
3940 IF (this%vcoord_out(k) >= &
3941 min(coord_in(kkdown), coord_in(kkdown+1)) .AND. &
3942 this%vcoord_out(k) < &
3943 max(coord_in(kkdown), coord_in(kkdown+1)))
THEN
3950 IF (kkup < inused)
THEN
3951 IF (this%vcoord_out(k) >= &
3952 min(coord_in(kkup), coord_in(kkup+1)) .AND. &
3953 this%vcoord_out(k) < &
3954 max(coord_in(kkup), coord_in(kkup+1)))
THEN
3964 IF (
c_e(kfound))
THEN
3965 z1 = real((this%vcoord_out(k) - coord_in(kfound-1))/ &
3966 (coord_in(kfound) - coord_in(kfound-1)))
3968 IF (vartype == var_dir360)
THEN
3969 field_out(i,j,k) = &
3970 interp_var_360(val_in(kfound-1), val_in(kfound), z1, z2)
3971 ELSE IF (vartype == var_press)
THEN
3972 field_out(i,j,k) = exp(val_in(kfound-1)*z2 + val_in(kfound)*z1)
3974 field_out(i,j,k) = val_in(kfound-1)*z2 + val_in(kfound)*z1
3984 DEALLOCATE(coord_in,val_in)
3989ELSE IF (this%trans%trans_type ==
'' .OR. this%trans%trans_type ==
'none')
THEN
3991 field_out(:,:,:) = field_in(:,:,:)
4001FUNCTION interp_var_360(v1, v2, w1, w2)
4002REAL,
INTENT(in) :: v1, v2, w1, w2
4003REAL :: interp_var_360
4007IF (abs(v1 - v2) > 180.)
THEN
4015 interp_var_360 = modulo(lv1*w2 + lv2*w1, 360.)
4017 interp_var_360 = v1*w2 + v2*w1
4020END FUNCTION interp_var_360
4023RECURSIVE FUNCTION find_prevailing_direction(v1, l, h, res, mask)
RESULT(prev)
4024REAL,
INTENT(in) :: v1(:,:)
4025REAL,
INTENT(in) :: l, h, res
4026LOGICAL,
INTENT(in),
OPTIONAL :: mask(:,:)
4034IF ((h - l) <= res)
THEN
4039IF (
PRESENT(mask))
THEN
4040 nl = count(v1 >= l .AND. v1 < m .AND. mask)
4041 nh = count(v1 >= m .AND. v1 < h .AND. mask)
4043 nl = count(v1 >= l .AND. v1 < m)
4044 nh = count(v1 >= m .AND. v1 < h)
4047 prev = find_prevailing_direction(v1, m, h, res)
4048ELSE IF (nl > nh)
THEN
4049 prev = find_prevailing_direction(v1, l, m, res)
4050ELSE IF (nl == 0 .AND. nh == 0)
THEN
4056END FUNCTION find_prevailing_direction
4059END SUBROUTINE grid_transform_compute
4067SUBROUTINE grid_transform_v7d_grid_compute(this, field_in, field_out, var, &
4070REAL,
INTENT(in) :: field_in(:,:)
4071REAL,
INTENT(out):: field_out(:,:,:)
4072TYPE(vol7d_var),
INTENT(in),
OPTIONAL :: var
4073REAL,
INTENT(in),
OPTIONAL,
TARGET :: coord_3d_in(:,:,:)
4075real,
allocatable :: field_in_p(:),x_in_p(:),y_in_p(:)
4076INTEGER :: inn_p, ier, k
4077INTEGER :: innx, inny, innz, outnx, outny, outnz
4080call l4f_category_log(this%category,l4f_debug,
"start v7d_grid_transform_compute")
4083field_out(:,:,:) = rmiss
4085IF (.NOT.this%valid)
THEN
4086 call l4f_category_log(this%category,l4f_error, &
4087 "refusing to perform a non valid transformation")
4092innx =
SIZE(field_in,1); inny = 1; innz =
SIZE(field_in,2)
4093outnx =
SIZE(field_out,1); outny =
SIZE(field_out,2); outnz =
SIZE(field_out,3)
4096IF (this%trans%trans_type ==
'vertint')
THEN
4098 IF (innz /= this%innz)
THEN
4099 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
4100 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
4101 t2c(this%innz)//
" /= "//t2c(innz))
4106 IF (outnz /= this%outnz)
THEN
4107 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
4108 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
4109 t2c(this%outnz)//
" /= "//t2c(outnz))
4114 IF (innx /= outnx .OR. inny /= outny)
THEN
4115 CALL l4f_category_log(this%category,l4f_error,
"vertical interpolation")
4116 CALL l4f_category_log(this%category,l4f_error,
"inconsistent hor. sizes: "//&
4117 t2c(innx)//
","//t2c(inny)//
" /= "//&
4118 t2c(outnx)//
","//t2c(outny))
4125 IF (innx /= this%innx .OR. inny /= this%inny)
THEN
4126 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
4127 CALL l4f_category_log(this%category,l4f_error,
"inconsistent input shape: "//&
4128 t2c(this%innx)//
","//t2c(this%inny)//
" /= "//&
4129 t2c(innx)//
","//t2c(inny))
4134 IF (outnx /= this%outnx .OR. outny /= this%outny)
THEN
4135 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
4136 CALL l4f_category_log(this%category,l4f_error,
"inconsistent output shape: "//&
4137 t2c(this%outnx)//
","//t2c(this%outny)//
" /= "//&
4138 t2c(outnx)//
","//t2c(outny))
4143 IF (innz /= outnz)
THEN
4144 CALL l4f_category_log(this%category,l4f_error,
"horizontal interpolation")
4145 CALL l4f_category_log(this%category,l4f_error,
"inconsistent vert. sizes: "//&
4146 t2c(innz)//
" /= "//t2c(outnz))
4154 call l4f_category_log(this%category,l4f_debug, &
4155 "start grid_transform_v7d_grid_compute "//trim(this%trans%trans_type)//
':'// &
4156 trim(this%trans%sub_type))
4159IF (this%trans%trans_type ==
'inter')
THEN
4161 IF (this%trans%sub_type ==
'linear')
THEN
4163#ifdef HAVE_LIBNGMATH
4165 ALLOCATE(field_in_p(innx*inny), x_in_p(innx*inny), y_in_p(innx*inny))
4167 inn_p = count(
c_e(field_in(:,k)))
4169 CALL l4f_category_log(this%category,l4f_info, &
4170 "Number of sparse data points: "//t2c(inn_p)//
','//t2c(
SIZE(field_in(:,k))))
4174 field_in_p(1:inn_p) = pack(field_in(:,k),
c_e(field_in(:,k)))
4175 x_in_p(1:inn_p) = pack(this%inter_xp(:,1),
c_e(field_in(:,k)))
4176 y_in_p(1:inn_p) = pack(this%inter_yp(:,1),
c_e(field_in(:,k)))
4178 IF (.NOT.this%trans%extrap)
THEN
4179 CALL nnseti(
'ext', 0)
4180 CALL nnsetr(
'nul', rmiss)
4183 CALL natgrids(inn_p, x_in_p, y_in_p, field_in_p, &
4184 this%outnx, this%outny, real(this%inter_x(:,1)), &
4185 REAL(this%inter_y(1,:)), field_out(1,1,k), ier)
4188 CALL l4f_category_log(this%category,l4f_error, &
4189 "Error code from NCAR natgrids: "//t2c(ier))
4195 CALL l4f_category_log(this%category,l4f_info, &
4196 "insufficient data in gridded region to triangulate")
4200 DEALLOCATE(field_in_p, x_in_p, y_in_p)
4203 CALL l4f_category_log(this%category,l4f_error, &
4204 "libsim compiled without NATGRIDD (ngmath ncarg library)")
4211ELSE IF (this%trans%trans_type ==
'boxinter' .OR. &
4212 this%trans%trans_type ==
'polyinter' .OR. &
4213 this%trans%trans_type ==
'vertint' .OR. &
4214 this%trans%trans_type ==
'metamorphosis')
THEN
4217 reshape(field_in, (/
SIZE(field_in,1), 1,
SIZE(field_in,2)/)), field_out, var, &
4222END SUBROUTINE grid_transform_v7d_grid_compute
4237ELEMENTAL FUNCTION hbilin (z1,z2,z3,z4,x1,y1,x3,y3,xp,yp)
RESULT(zp)
4238REAL,
INTENT(in) :: z1,z2,z3,z4
4239DOUBLE PRECISION,
INTENT(in):: x1,y1