libsim  Versione 7.2.4

◆ qc_compute_percentile()

subroutine, public modqccli::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 1618 del file modqccli.F90.

1619 indnetwork=1
1620 do 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
1665 end do
1666 
1667 deallocate (perc,mask,area)
1668 
1669 contains
1670 
1671 subroutine sub_perc(mymask)
1672 
1673 logical :: mymask(:,:,:)
1674 
1675  ! we want more than 30% data present and a number of data bigger than 100 (default)
1676 if &
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 
1683 if &
1684  ( c_e(lpresentnumb) .and. (count &
1685  (mymask .and. c_e(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:))) < lpresentnumb)&
1686  ) &
1687  return
1688 
1689 
1690 perc= stat_percentile(&
1691  pack(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:), &
1692  mask=mymask), &
1693  perc_vals)
1694 
1695 
1696 do 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)
1700 enddo
1701 
1702 
1703 end subroutine sub_perc
1704 
1705 
1706 end SUBROUTINE qc_compute_percentile
1707 
1708 
1709 SUBROUTINE qc_compute_normalizeddensityindex(this, perc_vals,cyclicdt,presentperc, presentnumb,data_normalized)
1710 
1711 TYPE(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
1715 real,intent(in) :: perc_vals(:)
1716 TYPE(cyclicdatetime),INTENT(in) :: cyclicdt
1717 real,optional :: presentperc
1718 integer,optional :: presentnumb
1719 logical,optional,intent(in) :: data_normalized
1720 
1721 integer :: indana,indtime,indvar,indnetwork,indlevel ,indtimerange ,inddativarr, indattr
1722 integer :: i,j,k,iana,narea
1723 real :: height
1724 TYPE(vol7d_var) :: var
1725 character(len=vol7d_ana_lenident) :: ident
1726 character(len=1) :: type
1727 integer :: areav(size(this%v7d%ana)),iclv(size(this%v7d%ana))
1728 logical,allocatable :: mask(:,:,:),maskplus(:,:,:), maskarea(:)
1729 integer,allocatable :: area(:)
1730 REAL, DIMENSION(:),allocatable :: ndi,limbins
1731 real :: lpresentperc
1732 integer :: lpresentnumb
1733 logical :: lnorm
1734 
1735 lnorm = optio_log(data_normalized)
1736 lpresentperc=optio_r(presentperc)
1737 lpresentnumb=optio_i(presentnumb)
1738 
1739 allocate (ndi(size(perc_vals)-1),limbins(size(perc_vals)))
1740 call delete(this%clima)
1741 CALL init(this%clima, time_definition=this%v7d%time_definition)
1742 call init(var, btable="B01192") ! MeteoDB station ID that here is the number of area
1743 
1744 if (.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 
1852 else
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 
1861 endif
1862 
1863 
1864 !!$do i=1,size(that%ana)
1865 !!$ call display(that%ana(i))
1866 !!$end do
1867 
1868 call 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 
1871 this%clima%level=this%v7d%level
1872 this%clima%timerange=this%v7d%timerange
1873 this%clima%dativar%r=this%v7d%dativar%r
1874 this%clima%dativarattr%r=this%clima%dativar%r
1875 call init(this%clima%datiattr%r(1), btable="*B33209") ! NDI order number
1876 this%clima%time(1)=cyclicdatetime_to_conventional(cyclicdt)
1877 
1878 call l4f_category_log(this%category,l4f_info,"vol7d_compute_ndi conventional datetime "//to_char(this%clima%time(1)))
1879 call init(this%clima%network(1),name="qcclima-ndi")
1880 
1881 call vol7d_alloc_vol(this%clima,inivol=.true.)
1882 
1883 allocate (mask(size(this%v7d%ana),size(this%v7d%time),size(this%v7d%network)))
1884 
1885 indtime=1
1886 indnetwork=1
1887 indattr=1
1888 do 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
1891  if (.NOT.(lnorm)) then
1892  do i=1,narea
1893 
1894  !this%v7d%voldatir(indana, indtime, indlevel, indtimerange, inddativarr, indnetwork)
Index method.

Generated with Doxygen.