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
|