libsim Versione 7.1.11

◆ qc_compute_percentile()

subroutine, public qc_compute_percentile ( type(qcclitype), intent(inout)  this,
real, dimension(:), intent(in)  perc_vals,
type(cyclicdatetime), intent(in)  cyclicdt,
real, optional  presentperc,
integer, optional  presentnumb 
)
Parametri
[in,out]thisvolume providing data to be computed and output data. Input data are not modified by the method, apart from performing a vol7d_alloc_vol on it
[in]perc_valspercentile values to use in compute, between 0. and 100.
[in]cyclicdtcyclic date and time
presentpercpercentual of data present for compute (default=not used)
presentnumbnumber of data present for compute (default=not used) $logical, optional :: height2level !< use conventional level starting from station height

Definizione alla linea 1614 del file modqccli.F90.

1615
1616allocate (mask(size(this%v7d%ana),size(this%v7d%time),size(this%v7d%network)))
1617
1618indtime=1
1619indnetwork=1
1620do inddativarr=1,size(this%v7d%dativar%r)
1621 do indtimerange=1,size(this%v7d%timerange)
1622 do indlevel=1,size(this%v7d%level) ! all stations, all times, all networks
1623 do i=1,narea
1624
1625 !this%v7d%voldatir(indana, indtime, indlevel, indtimerange, inddativarr, indnetwork)
1626
1627 !create mask only with valid time
1628 mask = spread(spread((this%v7d%time == cyclicdt ),1,size(this%v7d%ana)),3,size(this%v7d%network))
1629
1630#ifdef DEBUG
1631 CALL l4f_category_log(this%category, l4f_debug, 'count has value '//t2c(count &
1632 (mask .and. c_e(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:)))))
1633#endif
1634
1635 !delete in mask different area
1636 do j=1, size(mask,1)
1637 if (areav(j) /= area(i)) mask(j,:,:) =.false.
1638 end do
1639
1640 if (this%height2level) then
1641 allocate (maskplus(size(this%v7d%ana),size(this%v7d%time),size(this%v7d%network)))
1642 do k=1,cli_nlevel
1643#ifdef DEBUG
1644 CALL l4f_category_log(this%category, l4f_debug, 'k has value '//t2c(k))
1645#endif
1646
1647 do iana=1,size(mask,1)
1648 if (iclv(iana) /= k) maskplus(iana,:,:) =.false.
1649 if (iclv(iana) == k) maskplus(iana,:,:) = mask(iana,:,:)
1650 enddo
1651
1652 call sub_perc(maskplus)
1653
1654 enddo
1655 deallocate(maskplus)
1656 else
1657
1658 k=1
1659 call sub_perc(mask)
1660
1661 endif
1662 end do
1663 end do
1664 end do
1665end do
1666
1667deallocate (perc,mask,area)
1668
1669contains
1670
1671subroutine sub_perc(mymask)
1672
1673logical :: mymask(:,:,:)
1674
1675 ! we want more than 30% data present and a number of data bigger than 100 (default)
1676if &
1677 ( c_e(lpresentperc) .and. ((float(count &
1678 (mymask .and. c_e(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:)))&
1679 ) / &
1680 float(count(mymask))) < lpresentperc)) &
1681 return
1682
1683if &
1684 ( c_e(lpresentnumb) .and. (count &
1685 (mymask .and. c_e(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:))) < lpresentnumb)&
1686 ) &
1687 return
1688
1689
1690perc= stat_percentile(&
1691 pack(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:), &
1692 mask=mymask), &
1693 perc_vals)
1694
1695
1696do j=1,size(perc_vals)
1697 indana=(k-1)*size(perc_vals)*narea + (j-1)*narea + i
1698 this%extreme%voldatir(indana, indtime, indlevel, indtimerange, inddativarr, indnetwork)=&
1699 perc(j)
1700enddo
1701
1702
1703end subroutine sub_perc
1704
1705
1706end SUBROUTINE qc_compute_percentile
1707
1708
1709SUBROUTINE qc_compute_normalizeddensityindex(this, perc_vals,cyclicdt,presentperc, presentnumb,data_normalized)
1710
1711TYPE(qcclitype),INTENT(inout) :: this
1712!TYPE(timedelta),INTENT(in) :: step !< length of the step over which the statistical processing is performed
1713!TYPE(datetime),INTENT(in),OPTIONAL :: start !< start of statistical processing interval
1714!TYPE(datetime),INTENT(in),OPTIONAL :: stopp !< end of statistical processing interval
1715real,intent(in) :: perc_vals(:)
1716TYPE(cyclicdatetime),INTENT(in) :: cyclicdt
1717real,optional :: presentperc
1718integer,optional :: presentnumb
1719logical,optional,intent(in) :: data_normalized
1720
1721integer :: indana,indtime,indvar,indnetwork,indlevel ,indtimerange ,inddativarr, indattr
1722integer :: i,j,k,iana,narea
1723real :: height
1724TYPE(vol7d_var) :: var
1725character(len=vol7d_ana_lenident) :: ident
1726character(len=1) :: type
1727integer :: areav(size(this%v7d%ana)),iclv(size(this%v7d%ana))
1728logical,allocatable :: mask(:,:,:),maskplus(:,:,:), maskarea(:)
1729integer,allocatable :: area(:)
1730REAL, DIMENSION(:),allocatable :: ndi,limbins
1731real :: lpresentperc
1732integer :: lpresentnumb
1733logical :: lnorm
1734
1735lnorm = optio_log(data_normalized)
1736lpresentperc=optio_r(presentperc)
1737lpresentnumb=optio_i(presentnumb)
1738
1739allocate (ndi(size(perc_vals)-1),limbins(size(perc_vals)))
1740call delete(this%clima)
1741CALL init(this%clima, time_definition=this%v7d%time_definition)
1742call init(var, btable="B01192") ! MeteoDB station ID that here is the number of area
1743
1744if (.NOT.(lnorm)) then
1745
1746 type=cmiss
1747 indvar = index(this%v7d%anavar, var, type=type)
1748 indnetwork=1
1749
1750!if( ind /= 0 ) then
1751 select case (type)
1752 case("d")
1753 areav=integerdat(this%v7d%volanad(:,indvar,indnetwork),this%v7d%anavar%d(indvar))
1754 case("r")
1755 areav=integerdat(this%v7d%volanar(:,indvar,indnetwork),this%v7d%anavar%r(indvar))
1756 case("i")
1757 areav=integerdat(this%v7d%volanai(:,indvar,indnetwork),this%v7d%anavar%i(indvar))
1758 case("b")
1759 areav=integerdat(this%v7d%volanab(:,indvar,indnetwork),this%v7d%anavar%b(indvar))
1760 case("c")
1761 areav=integerdat(this%v7d%volanac(:,indvar,indnetwork),this%v7d%anavar%c(indvar))
1762 case default
1763 areav=imiss
1764 end select
1765!end if
1766
1767 allocate(maskarea(size(this%v7d%ana)))
1768 maskarea(:)= areav(:) /= imiss
1769 narea=count_distinct(areav,maskarea)
1770 allocate(area(narea))
1771 area=pack_distinct(areav,narea,maskarea)
1772 deallocate(maskarea)
1773
1774 if (this%height2level) then
1775 call vol7d_alloc(this%clima,nana=narea*(size(perc_vals)-1)*cli_nlevel)
1776 else
1777 call vol7d_alloc(this%clima,nana=narea*(size(perc_vals)-1))
1778 endif
1779
1780
1781
1782
1783 do i=1,narea
1784 do j=1,size(perc_vals)-1
1785 if (this%height2level) then
1786 do k=1,cli_nlevel
1787 write(ident,'("#",i2.2,2i3.3)')k,area(i),nint(perc_vals(j))
1788 call init(this%clima%ana((k-1)*(size(perc_vals)-1)*narea + (j-1)*narea + i),ident=ident,lat=0d0,lon=0d0)
1789 enddo
1790 else
1791 k=0
1792 write(ident,'("#",i2.2,2i3.3)')k,area(i),nint(perc_vals(j))
1793 call init(this%clima%ana((j-1)*narea+i),ident=ident,lat=0d0,lon=0d0)
1794 endif
1795 end do
1796 end do
1797
1798
1799
1800
1801
1802 if (this%height2level) then
1803
1804 call init(var, btable="B07030") ! height
1805
1806 type=cmiss
1807 indvar = index(this%v7d%anavar, var, type=type)
1808
1809 do k=1,size(this%v7d%ana)
1810 height=rmiss
1811
1812 ! here we take the height fron any network (the first network win)
1813 do indnetwork=1,size(this%v7d%network)
1814
1815 if( indvar > 0 .and. indnetwork > 0 ) then
1816 select case (type)
1817 case("d")
1818 height=realdat(this%v7d%volanad(k,indvar,indnetwork),this%v7d%anavar%d(indvar))
1819 case("r")
1820 height=realdat(this%v7d%volanar(k,indvar,indnetwork),this%v7d%anavar%r(indvar))
1821 case ("i")
1822 height=realdat(this%v7d%volanai(k,indvar,indnetwork),this%v7d%anavar%i(indvar))
1823 case("b")
1824 height=realdat(this%v7d%volanab(k,indvar,indnetwork),this%v7d%anavar%b(indvar))
1825 case("c")
1826 height=realdat(this%v7d%volanac(k,indvar,indnetwork),this%v7d%anavar%c(indvar))
1827 case default
1828 height=rmiss
1829 end select
1830 end if
1831
1832 if (c_e(height)) exit
1833 end do
1834
1835 if (c_e(height)) then
1836 iclv(k)=firsttrue(cli_level1 <= height .and. height <= cli_level2 )
1837 else
1838 iclv(k)=imiss
1839 endif
1840
1841#ifdef DEBUG
1842 CALL l4f_log(l4f_debug, 'qc_compute_NormalizedDensityIndex height has value '//t2c(height,"Missing"))
1843 CALL l4f_log(l4f_debug, 'for k having number '//t2c(k)//&
1844 ' iclv has value '//t2c(iclv(k)))
1845#endif
1846 end do
1847
1848
1849 endif
1850
1851
1852else
1853
1854 narea=1
1855 call vol7d_alloc(this%clima,nana=size(perc_vals)-1)
1856 do j=1,size(perc_vals)-1
1857 write(ident,'("#",i2.2,2i3.3)')0,0,nint(perc_vals(j))
1858 call init(this%clima%ana(j),ident=ident,lat=0d0,lon=0d0)
1859 enddo
1860
1861endif
1862
1863
1864!!$do i=1,size(that%ana)
1865!!$ call display(that%ana(i))
1866!!$end do
1867
1868call vol7d_alloc(this%clima,nlevel=size(this%v7d%level), ntimerange=size(this%v7d%timerange), &
1869 ndativarr=size(this%v7d%dativar%r), nnetwork=1,ntime=1,ndativarattrr=size(this%v7d%dativar%r),ndatiattrr=1)
1870
1871this%clima%level=this%v7d%level
1872this%clima%timerange=this%v7d%timerange
1873this%clima%dativar%r=this%v7d%dativar%r
1874this%clima%dativarattr%r=this%clima%dativar%r
1875call init(this%clima%datiattr%r(1), btable="*B33209") ! NDI order number
1876this%clima%time(1)=cyclicdatetime_to_conventional(cyclicdt)
1877
1878call l4f_category_log(this%category,l4f_info,"vol7d_compute_ndi conventional datetime "//to_char(this%clima%time(1)))
1879call init(this%clima%network(1),name="qcclima-ndi")
1880
1881call vol7d_alloc_vol(this%clima,inivol=.true.)
1882
1883allocate (mask(size(this%v7d%ana),size(this%v7d%time),size(this%v7d%network)))
1884
1885indtime=1
1886indnetwork=1
1887indattr=1
1888do inddativarr=1,size(this%v7d%dativar%r)
1889 do indtimerange=1,size(this%v7d%timerange)
1890 do indlevel=1,size(this%v7d%level) ! all stations, all times, all networks
Index method.

Generated with Doxygen.