mo_EDK.f90 8.27 KB
Newer Older
Stephan Thober's avatar
Stephan Thober committed
1
2
3
4
5
6
7
8
9
module mo_EDK

  implicit none

  private

  public :: EDK

contains
Stephan Thober's avatar
Stephan Thober committed
10
  !****************************************************************************
Matthias Zink's avatar
Matthias Zink committed
11
  !
Stephan Thober's avatar
Stephan Thober committed
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
  !  SUBROUTINE: External Drift Kriging
  !
  !  PURPOSE:  Perform EDK (daily)
  !            Output: output gridsize = cellFactor * gridsize DEM
  !  UPDATES
  !            Created        L. Samaniego            22.03.2006
  !            Last Update       Sa                   14.04.2006
  !            Last Update       Sa                   14.07.2010   DEM ckeck
  !            Last Update       Zi                   05.05.2011   IWD if No stations < 2
  !            Last Update       Zi                   07.02.2012   correct prec data < 0 --> = 0
  !****************************************************************************
  subroutine EDK(k, jStart, jEnd, dCS, MetSta, dS, cell, doOK)
    use mo_kind,    only : i4, dp, sp
    use mainVar,    only : MeteoStation, noDataValue, nSta, thresholdDist
    use kriging,    only : CellCoarser, dtoS
    use runControl, only : correctNeg, distZero

    implicit none

    ! input / output variables
    integer(i4),        intent(in)    :: k         ! cell id
    integer(i4),        intent(in)    :: jStart, jEnd
    real(dp),           intent(in)    :: dCS(:, :) ! distances matrix
    type(MeteoStation), intent(in)    :: MetSta(:) ! MeteoStation input
    type(dtoS),         intent(in)    :: dS(:)     ! distance among stations
    type(CellCoarser),  intent(inout) :: cell(:)   ! cell specification
    logical, optional,  intent(in)    :: doOK   ! switch do ordinary kriging

    ! local variables
    logical                         :: doOK_loc, calc_weights
    integer(i4)                     :: jd        ! julian day
    integer(i4)                     :: i, j, l, ll
    integer(i4)                     :: ii, jj
    integer(i4)                     :: Nk(nSta), Nk_old(nSta), nNmax
    real(dp), allocatable           :: A (:,:), B(:), X(:) , C(:,:)
    real(dp), allocatable           :: lamda(:)
    real(dp)                        :: sumLamda
    !
    ! Check which stations have valid data ( =/ -9, 999 etc. ) store them in Nk(:)
    l      = 0
    ll     = 0
    Nk     = 0
    Nk_old = 0

    ! switch ordinary kriging off if not explicitly given
    doOK_loc = .False.
    if (present(doOK)) doOK_loc = doOK

    ! IF NK changed -> re-estimate weights
    timeloop: do jd = jStart, jEnd

      Nk_old = Nk
      Nk = 0_i4
      l = 0
      ll = 0
      do i = 1, cell(k)%nNS
        j = cell(k)%listNS(i)
        if ( MetSta(j)%z(jd) /= noDataValue ) then

          !***      zero field     ***********************************
          if (MetSta(j)%z(jd) == 0.0_dp ) ll = ll + 1
          !**********************************************************
          l = l + 1
          Nk(l) = j
        end if
      end do
      nNmax = l
      if (all(Nk == Nk_old) .and. allocated(X)) then
        calc_weights = .False.
      else
        calc_weights = .True.
      end if


      !>>>>  no value ! avoid indetermination
      ! avoid 0 value calculations ll == nNmax
      ! avoid calculations where only 1 station is available
      ! avoid numerical instabilities nNmax == 2 (may happen that the solver matrix becomes singular)
      if (.not. ( ll == nNmax .or.  nNmax == 1 .or. nNmax == 2 ) ) then     

        if (calc_weights) then
          ! print *, 'cell:       ', k
          ! print *, '#neighbors: ', nNmax
          ! print *, 'Nk:     ', Nk
          ! print *, 'Nk_old: ', Nk_old
          call get_kriging_weights(X, nNmax, Nk, doOK_loc, dCS(k, :), dS, cell(k), MetSta)
        end if

        ! The BLUE of z is then:
        cell(k)%z(jd) = 0.
        do i=1,nNmax
          ii=Nk(i)
          cell(k)%z(jd) = cell(k)%z(jd) + X(i) * MetSta(ii)%z(jd)
        end do
        !
        ! only one station available /= 0
      else if (ll /= nNmax .AND. nNmax == 1) then
        ii = Nk(1)
        cell(k)%z(jd) = MetSta(ii)%z(jd)
        ! 
        ! for precipitation, distant values are set to zero
        if (distZero .and. dCS(k,ii) .gt. thresholdDist) then
          cell(k)%z(jd) = 0.0_dp
        end if
        !  
        ! if all stations have the value 0
      else if (ll == nNmax ) then 
        cell(k)%z(jd) = 0.0_dp
        ! avoid numerical instabilities --> IWD insverse weighted squared distance
        ! matrix of solver may become instable if only two or three stations are available 
      else if (ll /= nNmax .AND. nNmax == 2) then
        allocate (lamda(nNmax))
        do i=1,nNmax
          ii=Nk(i)
          lamda(i)=1/dCS(k,ii)/dCS(k,ii)
        end do
        sumLamda=sum(lamda)
        lamda=lamda/sum(lamda)
        cell(k)%z = 0.0_dp
        do i=1,nNmax
          ii=Nk(i)
          cell(k)%z(jd) = cell(k)%z(jd) + lamda(i) * MetSta(ii)%z(jd)
        end do
        deallocate (lamda)
      end if
      ! stop 'TESTING'
    end do timeloop
    if (allocated(X)) deallocate (X)
    !
    ! correct negative
142
    if (correctNeg) cell(k)%z = merge(0._sp, cell(k)%z, (cell(k)%z .gt. -9999._sp) .and. (cell(k)%z .lt. 0.))
Stephan Thober's avatar
Stephan Thober committed
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
    !
  end subroutine EDK

  subroutine get_kriging_weights(X, nNmax, Nk, doOK_loc, dCS, dS, cell, MetSta)

    use mo_kind,     only : dp, i4
    use mainVar,     only : MeteoStation, noDataValue, nSta, thresholdDist
    use kriging,     only : CellCoarser, dtoS
    use varfit,      only : beta
    use mo_setVario, only : tVar

    implicit none

    real(dp), allocatable, intent(out) :: X(:)
    integer(i4),           intent(in)  :: nNmax
    integer(i4),           intent(in)  :: Nk(nSta)
    logical,               intent(in)  :: doOK_loc
    real(dp),              intent(in)  :: dCS(:)    ! distances matrix
    type(dtoS),            intent(in)  :: dS(:)     ! distance among stations
    type(CellCoarser),     intent(in)  :: cell      ! cell specification
    type(MeteoStation),    intent(in)  :: MetSta(:) ! MeteoStation input

    ! local variables
    integer(i4)              :: i, j, ii, jj, info
    integer(i4), allocatable :: ipvt(:)
    real(dp),    allocatable :: A (:,:), B(:), C(:,:)

Matthias Zink's avatar
Matthias Zink committed
170
    !
Stephan Thober's avatar
Stephan Thober committed
171
172
173
174
175
176
    ! initialize matrices
    if (doOK_loc) then
      allocate (A(nNmax+1,nNmax+1), B(nNmax+1), X(nNmax+1), ipvt(nNmax + 1))
    else
      allocate (A(nNmax+2,nNmax+2), B(nNmax+2), X(nNmax+2), ipvt(nNmax + 2))
    end if
Matthias Zink's avatar
Matthias Zink committed
177
178
179
180
181
    !
    ! assembling the system of equations: OK
    A=0.0_dp

    ! loop over available stations nNmax
182
    do i = 1, nNmax
Matthias Zink's avatar
Matthias Zink committed
183
      ii = Nk(i)
Stephan Thober's avatar
Stephan Thober committed
184

Matthias Zink's avatar
Matthias Zink committed
185
      do j = i + 1, nNmax
Stephan Thober's avatar
Stephan Thober committed
186
187
        jj = Nk(j)

Matthias Zink's avatar
Matthias Zink committed
188
189
190
        ! available only the upper triangular matrix
        if (jj > ii) then
          A(i,j)=tVar(dS(ii)%S(jj),beta(1),beta(2),beta(3))
191
          A(j,i)=A(i,j)
Stephan Thober's avatar
Stephan Thober committed
192
          ! print *, 'A: ', A(i,j), dS(ii)%S(jj),beta(1),beta(2),beta(3)
Matthias Zink's avatar
Matthias Zink committed
193
194
        else
          A(i,j)=tVar(dS(jj)%S(ii),beta(1),beta(2),beta(3))
195
          A(j,i)=A(i,j)
Matthias Zink's avatar
Matthias Zink committed
196
197
        end if
      end do
Stephan Thober's avatar
Stephan Thober committed
198
      A(i,i) = tVar(0._dp, beta(1),beta(2),beta(3))
199

Matthias Zink's avatar
Matthias Zink committed
200
      A(i,nNmax+1) = 1.0_dp
201
      A(nNmax+1,i) = 1.0_dp
Stephan Thober's avatar
Stephan Thober committed
202
203
204
205
206
207
      if (.not. doOK_loc) A(i,nNmax+2) = MetSta(ii)%h
      if (.not. doOK_loc) A(nNmax+2,i) = MetSta(ii)%h

      B(i)=tVar(dCS(ii), beta(1),beta(2),beta(3))
      ! print *, 'B: ', B(i), dCS(k,ii), beta(1),beta(2),beta(3)
    end do
208

Matthias Zink's avatar
Matthias Zink committed
209
210
    !
    B(nNmax+1) = 1.0_dp
Stephan Thober's avatar
Stephan Thober committed
211
    if (.not. doOK_loc) B(nNmax+2) = cell%h
Matthias Zink's avatar
Matthias Zink committed
212

213
    ! NOTE: only the upper triangular matrix is needed!
214
    ! call D_LSASF (A, B, X)
215
    ! print *, '<<<<<<<<<<<<<<'
Stephan Thober's avatar
Stephan Thober committed
216
    ! print *, 'rhs = ', B(:10)
217
218
    C = A
    X = B
Stephan Thober's avatar
Stephan Thober committed
219
    ! print *, 'A = ', C(:10, 1)
220
    ! print *, '<<<<<<<<<<<<<<'
221
    call dgesv(size(A, 1), 1, A, size(A, 1), ipvt, B, size(A, 1), info)
222
    ! print *, '<<<<<<<<<<<<<<'
Stephan Thober's avatar
Stephan Thober committed
223
224
    ! print *, 'B = ', B(:10)
    ! print *, 'A = ', A(:10, 1)
225
    ! print *, '<<<<<<<<<<<<<<'
226
    if (maxval(abs(matmul(C, B) - X)) .gt. 1e-10) print *, 'maximum error: ', maxval(abs(matmul(C, B) - X))
227
    X = B
228
229
230
231
232
233
234
    if (info .ne. 0_i4) then
      print *, '***WARNING: calculation of weights failed'
    end if
    if (abs(sum(X(:nNmax)) - 1._dp) .gt. 1.e-4) then
      print *, '***WARNING: sum of weights is not 1, calculation of weights failed'
      print *, 'sum of weights: ', sum(X(:nNmax))
    end if
Stephan Thober's avatar
Stephan Thober committed
235
236
    ! print *, 'easting: ', cell%x
    ! print *, 'northing: ', cell%y
237
    ! print *, 'number of neighbors: ', nNmax
Stephan Thober's avatar
Stephan Thober committed
238
    ! print *, ''
Stephan Thober's avatar
Stephan Thober committed
239
    ! ! print *, 'ipvt: ', ipvt
240
241
    ! print *, 'info: ', info
    ! stop 'testing'
242
    deallocate (A, B, C, ipvt)
Stephan Thober's avatar
Stephan Thober committed
243
  end subroutine get_kriging_weights
Matthias Zink's avatar
Matthias Zink committed
244

Stephan Thober's avatar
Stephan Thober committed
245
end module mo_EDK