LCOV - code coverage report
Current view: top level - src - mo_edk_empvar.f90 (source / functions) Hit Total Coverage
Test: edk coverage Lines: 46 60 76.7 %
Date: 2024-03-11 14:23:05 Functions: 1 1 100.0 %

          Line data    Source code
       1             : !> \file    mo_edk_empvar.f90
       2             : !> \brief   \copybrief mo_edk_empvar
       3             : !> \details \copydetails mo_edk_empvar
       4             : 
       5             : !> \brief   This program calculates the empirical variogram
       6             : !> \details gamma(:,1) : distance
       7             : !!          gamme(:,2) : semivarance
       8             : !!          botm are normalized with gmax(:)
       9             : !> \author  Luis Samaniego
      10             : !> \date    22.11.2001
      11             : !> \date    19.02.2004
      12             : !> \copyright Copyright 2005-\today, the CHS Developers, Sabine Attinger: All rights reserved.
      13             : !! EDK is released under the LGPLv3+ license \license_note
      14             : module mo_edk_empvar
      15             : 
      16             :   implicit none
      17             : 
      18             :   private
      19             : 
      20             :   public :: EmpVar
      21             : 
      22             : contains
      23             : 
      24             :   !> \brief   calculate the empirical variogram
      25         365 :   subroutine EmpVar(jd, flagMax)
      26             :     use mainVar
      27             :     use mo_kind   , only       : i4, dp
      28             :     use kriging
      29             :     use VarFit
      30             :     implicit none
      31             :     integer(i4), intent(in)   :: jd
      32             :     logical, intent(in)       :: flagMax
      33             :     integer(i4)               :: i, j, k, nhk, ni
      34         365 :     real(dp)                  :: gk, delta
      35             : 
      36             :     ! variance h=0 (one pass algorithm)
      37        7300 :     do i = 1,nSta
      38        6935 :       if ( MetSta(i)%z(jd) ==  noDataValue ) cycle
      39        5840 :       nobs  = nobs + 1
      40        5840 :       delta = MetSta(i)%z(jd) - m0
      41        5840 :       m0    = m0 + delta / real(nobs, dp)
      42        7300 :       v0    = v0 + delta * (MetSta(i)%z(jd) - m0)
      43             :       !write (999,'(4i6, 3f15.4)') y,d,i, nobs, MetSta(i)%z(jd), m0, v0
      44             :     end do
      45             :     !
      46             :     !
      47             :     ! calculate the empiric variogram
      48             :     !
      49             :     ! selecting pair - bins
      50         365 :     if (jd == jStart) then
      51           1 :       nbins = ceiling(hMax/dh)
      52           5 :       if (.not. allocated (nH))   allocate (Nh(nbins),      gamma(nbins,2))
      53         151 :       Nh = 0
      54         303 :       gamma = 0.0_dp
      55             :     end if
      56             :     !
      57         365 :     print *, '***WARNING: Removal of outliers in the estimation of the variogram is activated'
      58        6935 :     do i=1,nSta-1
      59       69350 :       do j=i+1,nSta
      60       68985 :         if (edk_dist%getZ2(i,j,jd) /=  noDataValue ) then
      61             :           ! take values up to max distance
      62       43800 :           if (edk_dist%getSS(i,j) > hMax ) cycle
      63             :           ! ! remove outliers for the estimation of the variogram
      64             :           ! if (flagVarTyp == 2 ) then
      65       43800 :             if (  dabs( MetSta(i)%h - MetSta(j)%h ) / edk_dist%getSS(i,j)  > gradE ) then
      66             :               ! write(999,*), 'pair removed', MetSta(i)%id, MetSta(j)%id
      67             :                 cycle
      68             :             end if
      69             :           ! end if
      70             :           !
      71       43800 :           k=max(1, ceiling(edk_dist%getSS(i,j)/dh))
      72       43800 :           Nh(k)=Nh(k)+1
      73       43800 :           gamma(k,2)=gamma(k,2) + edk_dist%getZ2(i,j,jd)
      74             :         end if
      75             :       end do
      76             :     end do
      77             :     !
      78             :     ! only at the end
      79         365 :     if (jd == jEnd) then
      80             :       !
      81             :       ! consolidate bins N(h)>=30
      82             :       !forward
      83         150 :       do k=1,nbins-1
      84         149 :         nhk=Nh(k)
      85         149 :         gk=gamma(k,2)
      86         150 :         if ((nhk>0) .and. (nhk<30)) then
      87           0 :           Nh(k)=-9
      88           0 :           gamma(k,2)=0.0_dp
      89           0 :           if (Nh(k+1)>0) then
      90           0 :             Nh(k+1)=Nh(k+1)+nhk
      91           0 :             gamma(k+1,2)=gamma(k+1,2)+gk
      92             :           else
      93           0 :             Nh(k+1)=nhk
      94           0 :             gamma(k+1,2)=gk
      95             :           end if
      96             :         end if
      97             :       end do
      98             :       !backward
      99         150 :       do k=nbins,2,-1
     100         149 :         nhk=Nh(k)
     101         149 :         gk=gamma(k,2)
     102         150 :         if ((nhk>0) .and. (nhk<30)) then
     103           0 :           Nh(k)=-9
     104           0 :           gamma(k,2)=0_dp
     105           0 :           if (Nh(k-1)>0) then
     106           0 :             Nh(k-1)=Nh(k-1)+nhk
     107           0 :             gamma(k-1,2)=gamma(k-1,2)+gk
     108             :           else
     109           0 :             Nh(k-1)=nhk
     110           0 :             gamma(k-1,2)=gk
     111             :           end if
     112             :         end if
     113             :       end do
     114             :       !
     115           1 :       ni=0
     116             :       ! estimate semi-variance gamma(i,1)=h, gamma(i,2)=g(h)
     117         151 :       do k=1,nbins
     118         151 :         if (Nh(k) > 0) then
     119          29 :           ni=ni+1
     120             :           ! Classsical
     121          29 :           gamma(k,2)=gamma(k,2)/2.0_dp/real(Nh(k), dp)
     122             :           !
     123             :           ! Cressi-Hawkins: adjust bias
     124             :           !gamma(k,2)=0.5_dp/(0.457_dp+0.494_dp/dfloat(Nh(k)))*(gamma(k,2)/dfloat(Nh(k)))**4
     125             :           !
     126          29 :           gamma(k,1)=real(k, dp)*dh-real(ni, dp)*dh/2.0_dp
     127          29 :           ni=0
     128             :         else
     129         121 :           ni=ni+1
     130         121 :           gamma(k,1)=-9.0_dp
     131             :         end if
     132             :       end do
     133             :       !
     134             :       ! scaling
     135           1 :       if (flagMax) gmax=maxval(gamma, dim=1)
     136             : 
     137         151 :       do k=1,nbins
     138         151 :         if (gamma(k,1) > 0.0_dp) then
     139          29 :           gamma(k,1) = gamma(k,1)/gmax(1)
     140          29 :           gamma(k,2) = gamma(k,2)/gmax(2)
     141             :         end if
     142             :       end do
     143             :       !
     144             :       !keep variogram
     145             :     end if
     146             : 
     147         365 :   end subroutine EmpVar
     148             : 
     149             : end module mo_edk_empvar

Generated by: LCOV version 1.16