Line data Source code
1 : !> \file mo_smi_drought_evaluation.f90
2 : !> \brief \copybrief mo_smi_drought_evaluation
3 : !> \details \copydetails mo_smi_drought_evaluation
4 :
5 : !> \brief drought evaluation for SMI
6 : !> \copyright Copyright 2005-\today, the CHS Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
7 : !! SMI is released under the LGPLv3+ license \license_note
8 : module mo_smi_drought_evaluation
9 :
10 : USE mo_kind, only : i4, sp, dp
11 :
12 : IMPLICIT NONE
13 :
14 : PRIVATE
15 :
16 : PUBLIC :: droughtIndicator
17 : PUBLIC :: ClusterEvolution
18 : PUBLIC :: ClusterStats
19 : PUBLIC :: calSAD
20 :
21 : ! ------------------------------------------------------------------
22 :
23 : CONTAINS
24 :
25 : ! ------------------------------------------------------------------
26 :
27 : !> \brief Find drought clusters and statistics
28 : !! 1) truncate SMIp < SMI_th
29 : !> \author Luis Samaniego
30 : !> \date 28.02.2011
31 : !> \date 05.10.2011
32 1 : subroutine droughtIndicator( SMI, mask, SMI_thld, cellCoor , SMIc)
33 :
34 : use mo_constants, only : nodata_sp, nodata_i4
35 : use mo_utils, only : lesserequal, notequal
36 :
37 : implicit none
38 :
39 : ! input variable
40 : real(sp), dimension(:,:), intent(in) :: SMI
41 : logical, dimension(:,:), intent(in) :: mask
42 : real(sp), intent(in) :: SMI_thld
43 : integer(i4), dimension(:,:), allocatable, intent(out) :: cellCoor
44 : integer(i4), dimension(:,:,:), allocatable, intent(out) :: SMIc ! Drought indicator
45 :
46 : ! local variables
47 : real(sp), dimension(size(mask, dim=1),&
48 2 : size(mask, dim=2)) :: dummy_2d_sp
49 : integer(i4) :: i, j, m, k
50 : integer(i4) :: nrows
51 : integer(i4) :: ncols
52 : integer(i4) :: nMonths
53 :
54 : ! initialize
55 1 : nrows = size( mask, 1 )
56 1 : ncols = size( mask, 2 )
57 1 : nMonths = size( SMI, 2 )
58 : !
59 1 : print *, 'Threshold for drought identification: ', SMI_thld
60 5 : allocate ( SMIc( nrows, ncols, nMonths) )
61 :
62 361 : do m = 1, nMonths
63 360 : dummy_2d_sp = unpack(SMI(:,m), mask, nodata_sp)
64 : ! filter for possible error within domain
65 : ! where (SMI(:,:,m) .le. SMI_thld .and. SMI(:,:,m) .ne. nodata_sp )
66 72361 : where ( (lesserequal(dummy_2d_sp, SMI_thld)) .and. (notequal(dummy_2d_sp, nodata_sp)) )
67 : SMIc(:,:,m) = 1
68 : elsewhere ( notequal(dummy_2d_sp, nodata_sp) )
69 : SMIc(:,:,m) = 0
70 : elsewhere
71 : SMIc(:,:,m) = nodata_i4
72 : end where
73 : end do
74 :
75 : !
76 43 : allocate (cellCoor( count(mask),2 ))
77 1 : k = 0
78 6 : do j=1,ncols
79 41 : do i=1,nrows
80 40 : if ( mask(i,j)) then
81 22 : k = k + 1
82 22 : cellCoor(k,1) = i
83 22 : cellCoor(k,2) = j
84 : ! if (k==7645.or.k==2479) print*, i,j,k
85 : end if
86 : end do
87 : end do
88 : !
89 1 : end subroutine droughtIndicator
90 :
91 : !> \brief cluster drought clusters in space and time
92 : !> \date developed in Budapest, 10-11.03.2011
93 1 : subroutine ClusterEvolution( SMIc, nrows, ncols, nMonths, nCells, cellCoor, nCellInter, thCellClus)
94 :
95 : use mo_orderpack, only : sort
96 : use mo_constants, only : nodata_i4
97 : use mo_smi_write, only : idCluster, &
98 : shortCnoList, &
99 : nClusters
100 :
101 :
102 : implicit none
103 :
104 : ! input variables
105 : integer(i4), dimension(:,:,:), intent(inout) :: SMIc ! Drought indicator
106 : integer(i4), intent(in) :: nrows
107 : integer(i4), intent(in) :: ncols
108 : integer(i4), intent(in) :: nMonths
109 : integer(i4), intent(in) :: nCells ! number of effective cells
110 : integer(i4), dimension(:,:), allocatable, intent(in) :: cellCoor
111 : integer(i4), intent(in) :: thCellClus ! treshold for cluster formation in space
112 : integer(i4), intent(in) :: nCellInter ! number cells for joining clusters in time
113 :
114 :
115 : ! local variables
116 : integer(i4) :: i, j, t
117 1 : integer(i4), dimension(:), allocatable :: nC
118 : integer(i4), parameter :: factor = 1000
119 : integer(i4) :: ncInter
120 1 : integer(i4), dimension(:,:), allocatable :: cno
121 1 : integer(i4), dimension(:), allocatable :: cnoList
122 1 : integer(i4), dimension(:), allocatable :: vec
123 : integer(i4) :: maxNc, nTotal, idRep
124 1 : integer(i4), dimension(:,:), allocatable :: idCnew
125 :
126 :
127 5 : if (.not. allocated(idCluster) ) allocate (idCluster(nrows, ncols, nMonths))
128 3 : if (.not. allocated(nC) ) allocate (nC(nMonths))
129 : if ( allocated(cno) ) deallocate (cno, vec, cnoList)
130 4 : if (.not. allocated(idCnew)) allocate (idCnew(nrows, ncols))
131 :
132 : ! ---------------------------------------------
133 : ! CLUSTERING IN SPACE and remove small clusters
134 : ! ---------------------------------------------
135 361 : do t=1,nMonths
136 14760 : if ( count(SMIc(:,:,t) == 1) < thCellClus) then
137 289 : nC(t) = 0
138 11849 : idCluster(:,:,t) = nodata_i4
139 : cycle
140 : end if
141 72 : call findClusters (cellCoor, thCellClus, t, idCluster(:,:,t), nC(t), nrows, ncols, nCells, SMIc)
142 : !print*, 'Finding clusters time :', t, nC(t)
143 : end do
144 1 : print*, 'Clusters in space done ...'
145 : !
146 : ! maximum number of clusters at all time steps
147 361 : maxNc = maxval(nC(:))
148 1 : nTotal = nMonths*maxNc
149 7 : allocate ( cno(nMonths,maxNc), vec(nTotal), cnoList(nTotal) )
150 362 : cno = -9
151 : !
152 : ! set unique cluster numbers and store them
153 : ! NOTE
154 : ! e.g. month*factor + running nr.,
155 : ! 5*1000+1 = 5001 => fisrt cluster in 5th month
156 : !
157 361 : do t=1,nMonths
158 360 : if (nC(t) == 0) cycle
159 139 : do i=1, nC(t)
160 69 : cno(t,i) = t*factor + i
161 3189 : where ( idCluster(:,:,t) == i ) idCluster(:,:,t) = cno(t,i)
162 : end do
163 : end do
164 : ! ---------------------------------
165 : ! FIND CLUSTERS IN TIME
166 : ! determine intersection sets, join
167 : ! ---------------------------------
168 360 : do t=2,nMonths
169 428 : do i=1,nC(t)
170 68 : if (cno(t,i) == -9) cycle
171 477 : do j=1,nC(t-1)
172 50 : if (cno(t-1,j) == -9) cycle
173 2050 : ncInter = count ( idCluster(:,:,t) == cno(t,i) .and. idCluster(:,:,t-1) == cno(t-1,j) )
174 118 : if ( ncInter >= nCellInter ) then
175 : ! renumber all from 1 to t
176 0 : where ( idCluster(:,:,1:t) == cno(t-1,j) ) idCluster(:,:,1:t) = cno(t,i)
177 : ! rename cluster id from cno
178 0 : where ( cno(1:t,:) == cno(t-1,j) ) cno(1:t,:) = cno(t,i)
179 : end if
180 : end do
181 : end do
182 : end do
183 1 : print*, 'Clustering in time done ...'
184 : ! --------------------------
185 : ! COMPILE CLUSTER IDS + list
186 : ! --------------------------
187 361 : vec = -9
188 362 : cnoList = pack (cno, mask = cno > 0, vector=vec)
189 : !
190 : ! 1. sort
191 1 : idRep = -9
192 361 : forall(i=1:nTotal) vec(i) = i
193 : ! call SVIGN(nTotal,cnoList,cnoList)
194 1 : call sort(cnoList)
195 : ! 2. remove repeated values
196 360 : do i=nTotal, 2,-1
197 359 : if (cnoList(i) < 0 ) cycle
198 69 : if (.not. cnoList(i-1) == cnoList(i) ) then
199 69 : if (idRep > 0) idRep = -9
200 : cycle
201 : end if
202 0 : if (idRep == -9) idRep = cnoList(i)
203 291 : cnoList(i) = -9
204 : end do
205 : ! call SVIGN(nTotal,cnoList,cnoList)
206 1 : call sort(cnoList)
207 : !
208 : ! 3. final consolidated list
209 361 : nClusters = count(cnoList > 0)
210 3 : allocate ( shortCnoList(nClusters) )
211 361 : shortCnoList = pack(cnoList, mask = cnoList > 0)
212 : !
213 1 : print*, '# Consolidated Clusters :', nClusters
214 : !
215 1 : deallocate(vec, nC, cnoList, cno )
216 1 : end subroutine ClusterEvolution
217 :
218 : !> \brief SVAT statistics
219 1 : subroutine ClusterStats( SMI, mask, nrows, ncols, nMonths, nCells, SMI_thld )
220 :
221 : use mo_orderpack, only : sort_index
222 : use mo_constants, only: nodata_sp
223 : use mo_smi_write, only: aDA, aDD, TDM, DTMagEvol, DAreaEvol, &
224 : nClusters, idCluster, shortCnoList, &
225 : dASevol, nBasins, nEvents, eIdPerm, eventId
226 :
227 : implicit none
228 :
229 : ! input variables
230 :
231 : real(sp), dimension(:,:), intent(in) :: SMI
232 : logical, dimension(:,:), intent(in) :: mask
233 : integer(i4), intent(in) :: nrows
234 : integer(i4), intent(in) :: ncols
235 : integer(i4), intent(in) :: nMonths
236 : integer(i4), intent(in) :: nCells ! number of effective cell
237 : ! integer(i4), dimension(:,:), intent(in) :: Basin_Id ! IDs for basinwise drought analysiss
238 : real(sp), intent(in) :: SMI_thld ! SMI threshold for clustering
239 :
240 : ! local variables
241 2 : real(sp), dimension(nrows, ncols) :: dummy_2d_sp
242 : integer(i4) :: ic, i
243 : integer(i4) :: t
244 1 : integer(i4) , dimension(:), allocatable :: counterA
245 1 : integer(i4) , dimension(:,:), allocatable :: aDDG
246 1 : real(dp) , dimension(:,:), allocatable :: mSev
247 :
248 : integer(i4) :: eCounter
249 :
250 1 : integer(i4), dimension(:), allocatable :: vec
251 :
252 4 : allocate ( aDDG (nrows, ncols) )
253 4 : allocate ( DAreaEvol(nMonths,nClusters) )
254 3 : allocate ( DTMagEvol(nMonths,nClusters) )
255 3 : allocate ( aDD (nClusters) )
256 2 : allocate ( aDA (nClusters) )
257 3 : allocate ( counterA (nClusters) )
258 2 : allocate ( TDM (nClusters) )
259 4 : allocate ( dASevol (nMonths,2,nBasins+1) ) ! (area, severity) whole Germany => nBasins+1
260 4 : allocate ( mSev (nrows, ncols) ) ! mean severity
261 : !
262 24910 : DAreaEvol = 0.0_dp
263 24910 : DTMagEvol = 0.0_dp
264 70 : counterA = 0
265 70 : do ic = 1, nClusters
266 : !print*, 'Statistics of cluster : 'shortCnoList(ic)
267 2829 : aDDG = 0
268 24909 : do t = 1, nMonths
269 1018440 : where (idCluster(:,:,t) == shortCnoList(ic) )
270 : ! duration
271 : aDDG = aDDG + 1
272 : end where
273 : !
274 : ! drought area evolution
275 1018440 : DAreaEvol(t,ic) = real( count(idCluster(:,:,t) == shortCnoList(ic) ), dp) / real(nCells, dp)
276 24840 : if (DAreaEvol(t,ic) > 0.0_dp) counterA(ic) = counterA(ic) + 1
277 :
278 : ! total magnitude (NEW SM_tr -SMI) !!!
279 24840 : dummy_2d_sp = unpack(SMI(:,t), mask, nodata_sp)
280 1018509 : DTMagEvol(t,ic) = sum( (SMI_thld - dummy_2d_sp), mask = idCluster(:,:,t) == shortCnoList(ic) )
281 : end do
282 :
283 : ! AVERAGE (MEAN) DURATION
284 5658 : aDD(ic) = real( sum(aDDG, mask = aDDG > 0),dp ) / real( count(aDDG > 0), dp)
285 :
286 : ! TOTAL MAGNITUD (sum over space and time of SMI over all
287 : ! cells afected by the event ic )
288 24910 : TDM(ic) = sum( DTMagEvol(:,ic) )
289 : end do
290 : ! AVERAGE (MEAN) DROUGHT AREA per cluster event
291 24911 : aDA = sum(DAreaEvol, DIM=1) / real(counterA, dp)
292 : !
293 : ! (NEW) Time evolution of the area and monthly severity
294 : ! whole domain
295 361 : do t = 1, nMonths
296 14760 : dASevol(t,1,nBasins+1) = real( count( idCluster(:,:,t) > 0 ), dp ) / real( nCells, dp) * 1e2_dp
297 360 : dummy_2d_sp = unpack(SMI(:,t), mask, nodata_sp)
298 14760 : where (idCluster(:,:,t) > 0 )
299 : mSev(:,:) = 1.0_dp - real(dummy_2d_sp, dp)
300 : end where
301 14761 : dASevol(t,2,nBasins+1) = sum( mSev(:,:), mask = idCluster(:,:,t) > 0 ) / real( nCells, dp)
302 : end do
303 : ! ! evolution at basin wise
304 : ! do t = 1, nMonths
305 : ! do i = 1, nBasins
306 : ! dASevol(t,1,i) = real ( count( idCluster(:,:,t) > 0 .and. Basin_Id(:,:) == i ), dp) &
307 : ! / real(count(Basin_Id(:,:) == i), dp) * 1e2_dp
308 : ! where (idCluster(:,:,t) > 0 .and. Basin_Id(:,:) == i )
309 : ! mSev(:,:) = 1.0_dp - SMI(:,:,t)
310 : ! end where
311 : ! dASevol(t,2,i) = sum( mSev(:,:), mask = idCluster(:,:,t) > 0 .and. Basin_Id(:,:) == i ) &
312 : ! / real(count(Basin_Id(:,:) == i), dp)
313 : ! end do
314 : ! end do
315 :
316 :
317 : !!! ORIGINALLY in SAD anaylysis
318 24910 : nEvents = count (DAreaEvol .gt. 0.0_dp )
319 3 : allocate ( eIdPerm (nEvents) )
320 3 : allocate ( eventId (nEvents, 3) ) ! dim1: running Nr.
321 : ! dim2: 1 == cluster Id,
322 : ! 2 == month of ocurrence,
323 : ! 3 == number of cells
324 :
325 : ! keep event identification
326 1 : eCounter = 0
327 70 : do ic = 1, nClusters
328 24910 : do t = 1, nMonths
329 24909 : if (DAreaEvol(t,ic) > 0.0_dp ) then
330 69 : eCounter = eCounter + 1
331 69 : eventId(eCounter,1) = shortCnoList(ic)
332 69 : eventId(eCounter,2) = t
333 2829 : eventId(eCounter,3) = count( idCluster(:,:,t) == shortCnoList(ic) )
334 : end if
335 : end do
336 : end do
337 : ! sort event in order of area magnitude
338 70 : forall(i=1:nEvents) eIdPerm(i) = i
339 3 : allocate (vec(nEvents))
340 : ! call SVIGP (nEvents, eventId(:,3), vec, eIdPerm)
341 71 : eIdPerm = sort_index( eventID(:,3) )
342 1 : deallocate (vec)
343 :
344 1 : print*, 'Cluster statistics were estimated ... '
345 1 : deallocate ( counterA, aDDG, mSev )
346 1 : end subroutine ClusterStats
347 :
348 : !> \brief SAD analysis
349 0 : subroutine calSAD(SMI, mask, iDur, nrows, ncols, nMonths, nCells, deltaArea, cellsize)
350 :
351 : use mo_orderpack, only : sort !, sort_index
352 : use mo_percentile, only : percentile
353 : use mo_constants, only : nodata_dp
354 : use mo_smi_write, only : shortCnoList, &
355 : nInterArea, nEvents, nClusters, idCluster, &
356 : SAD, SADperc, DAreaEvol, severity, nDsteps, &
357 : durList, eventId, eIdPerm, &
358 : nLargerEvents, nQProp, QProp
359 :
360 : implicit none
361 :
362 : ! input variable
363 : real(sp), dimension(:,:), intent(in) :: SMI
364 : logical, dimension(:,:), intent(in) :: mask
365 : integer(i4), intent(in) :: iDur
366 : integer(i4), intent(in) :: nrows
367 : integer(i4), intent(in) :: ncols
368 : integer(i4), intent(in) :: nMonths
369 : integer(i4), intent(in) :: nCells ! number of effective cells
370 : integer(i4), intent(in) :: deltaArea ! number of cells per area interval
371 : real(sp), intent(in) :: cellsize
372 :
373 : ! local variable
374 0 : real(dp), dimension(nrows,ncols, nMonths) :: SMI_unpack
375 : integer(i4) :: t, i, k
376 : ! integer(i4) :: ic
377 : integer(i4) :: iDc
378 : integer(i4) :: d
379 : integer(i4) :: ms, me
380 : integer(i4) :: ke
381 : integer(i4) :: eC
382 : integer(i4) :: eCounter
383 :
384 : integer(i4) :: ncic
385 : integer(i4) :: nIntA
386 0 : real(dp), dimension(:), allocatable :: sevP
387 : !
388 : integer(i4) :: nObs
389 : ! integer(i4), dimension(:), allocatable :: vec
390 : !
391 0 : do i = 1, size(SMI,2)
392 0 : SMI_unpack(:,:,i) = unpack(real(SMI(:,i), dp), mask, nodata_dp)
393 : end do
394 :
395 0 : if (iDur == 1) then
396 : ! set up SAD curves
397 0 : nInterArea = int( ceiling( real(nCells,dp) / real(deltaArea, dp) ), i4 )
398 : !
399 : ! total number of events to be evaluated
400 : ! nEvents = count (DAreaEvol > 0.0_dp )
401 0 : nLargerEvents = nEvents
402 : ! allocate ( eventId (nEvents, 3) ) ! dim1: running Nr.
403 : ! ! dim2: 1 == cluster Id,
404 : ! ! 2 == month of ocurrence,
405 : ! ! 3 == number of cells
406 : !
407 : ! allocate ( eIdPerm (nEvents) )
408 0 : allocate ( SAD (nInterArea, 2, nLargerEvents ) ) ! dim2: 1 == area, 2 == Severity
409 0 : allocate ( SADperc (nInterArea, nQProp ) )
410 : !
411 : ! ! keep event identification
412 : ! eCounter = 0
413 : ! do ic = 1, nClusters
414 : ! do t = 1, nMonths
415 : ! if (DAreaEvol(t,ic) > 0.0_dp ) then
416 : ! eCounter = eCounter + 1
417 : ! eventId(eCounter,1) = shortCnoList(ic)
418 : ! eventId(eCounter,2) = t
419 : ! eventId(eCounter,3) = count( idCluster(:,:,t) == shortCnoList(ic) )
420 : ! end if
421 : ! end do
422 : ! end do
423 : ! ! sort event in order of magnitude
424 : ! forall(i=1:nEvents) eIdPerm(i) = i
425 : ! allocate (vec(nEvents))
426 : ! ! call SVIGP (nEvents, eventId(:,3), vec, eIdPerm)
427 : ! eIdPerm = sort_index( eventID(:,3) )
428 : ! deallocate (vec)
429 : end if
430 : !
431 : ! SAD
432 0 : print*, 'SAD for duration started: ', durList (iDur)
433 0 : SAD = -9.0_dp
434 0 : SADperc = -9.0_dp
435 0 : nDsteps = ceiling (real(nMonths,dp) / real(durList (iDur),dp) )
436 0 : if ( allocated (severity) ) deallocate(severity)
437 0 : allocate ( severity (nrows, ncols, nDsteps ) )
438 0 : severity = nodata_dp
439 : !
440 : ! estimate severities for a given duration over all time steps
441 0 : do d = 1, nDsteps
442 0 : ms = (d-1)*durList(iDur) + 1
443 0 : me = ms + durList(iDur) - 1
444 0 : if (me > nMonths) me = nMonths
445 : !where ( ne(SMI_unpack(:,:,1), nodata_dp) )
446 0 : where ( mask )
447 : severity(:,:,d) = 1.0_dp - sum( SMI_unpack(:,:,ms:me), DIM=3 ) / real( me-ms+1, dp )
448 : end where
449 : end do
450 : !
451 : ! estimate curves for larger events ONLY
452 0 : do eC = 1, nLargerEvents
453 0 : eCounter = eIdPerm( nEvents + 1 - eC )
454 0 : iDc = eventId(eCounter,1)
455 0 : t = eventId(eCounter,2)
456 0 : d = (t-1)/durList (iDur) + 1
457 : ! number of cells of the event in eCounter
458 0 : ncic = eventId(eCounter,3)
459 0 : if (allocated (sevP)) deallocate(sevP)
460 0 : allocate (sevP(ncic))
461 0 : sevP = pack (severity(:,:,d), mask = idCluster(:,:,t) == iDc )
462 : !
463 : ! ----------------------------------------------
464 : ! fast approach
465 : ! rank severities and estimate cummulative areas
466 : ! ----------------------------------------------
467 : ! call DSVRGN (ncic,sevP,sevP)
468 0 : call sort( sevP )
469 0 : nIntA = int( ceiling( real(ncic,dp) / real(deltaArea, dp) ), i4 )
470 0 : do k = 1, nIntA
471 0 : ke = k*deltaArea
472 : if (ke > ncic ) ke = ncic
473 0 : SAD(k, 1, eC ) = real(ke, dp) * real(cellsize,dp)**2.0_dp ! cumulative area
474 0 : SAD(k, 2, eC ) = sum( sevP(ncic+1-ke:ncic) ) / real(ke, dp) ! sorted by incresing value
475 : end do
476 : end do
477 : !
478 : ! estimate percentilles for all intervals
479 0 : do k = 1, nInterArea
480 0 : nObs = count(SAD(k,2,:) > 0.0_dp)
481 0 : if (.not. nObs > 1) cycle
482 0 : if (allocated (sevP)) deallocate(sevP)
483 0 : allocate (sevP(nObs))
484 0 : sevP = pack (SAD(k,2,:), mask = SAD(k,2,:) > 0.0_dp )
485 : ! call DEQTIL (nObs, sevP, nQprop, Qprop, SADperc(k,:), Xlo, Xhi, nMiss)
486 0 : SADperc(k,:) = percentile( sevP, Qprop )
487 : end do
488 : !
489 0 : if ( allocated(sevP) ) deallocate ( sevP )
490 0 : end subroutine calSAD
491 :
492 : !> \brief clustering in space
493 : !> \details assign a unique running number at time t
494 71 : subroutine findClusters (cellCoor, thCellClus, t,iC,nCluster, nrows, ncols, nCells, SMIc)
495 :
496 : use mo_constants, only : nodata_i4
497 :
498 : implicit none
499 :
500 : integer(i4), dimension(:,:), intent(in) :: cellCoor
501 : integer(i4), intent(in) :: thCellClus ! treshold for cluster formation
502 : integer(i4), intent(in) :: nrows
503 : integer(i4), intent(in) :: ncols
504 : integer(i4), intent(in) :: nCells ! number of effective cells
505 : integer(i4), intent(in) :: t ! time step
506 : integer(i4), dimension(:,:,:), intent(inout) :: SMIc ! SMI indicator
507 : integer(i4), intent(out) :: nCluster ! number of clusters larger than a threshold
508 : integer(i4), dimension(nrows,ncols), intent(out) :: iC
509 :
510 : ! local variables
511 : integer(i4) :: i, j, k, klu, klu1
512 : integer(i4) :: iul, idr, jul, jdr
513 : integer(i4) :: krow, kcol
514 71 : integer(i4), dimension(:), allocatable :: cno, vec
515 71 : integer(i4), dimension(:), allocatable :: nCxCluster
516 : integer(i4) :: nClusterR, ncc
517 : !
518 2911 : iC = nodata_i4
519 71 : nCluster = 0
520 1633 : do k = 1, nCells
521 1562 : krow = cellCoor(k,1)
522 1562 : kcol = cellCoor(k,2)
523 1562 : iul = max(krow-1,1)
524 1562 : idr = min(krow+1,nrows)
525 1562 : jul = kcol
526 1562 : jdr = min(kcol+1,ncols)
527 : ! SMIc of k
528 1562 : klu = SMIc(krow, kcol, t)
529 1562 : if (klu /= 1) cycle
530 1244 : if ( (klu == 1) .and. (iC(krow, kcol) == nodata_i4) ) then
531 87 : nCluster=nCluster+1
532 87 : iC(krow, kcol) = nCluster
533 : end if
534 3581 : do j=jul, jdr
535 10244 : do i= iul, idr
536 6416 : if (iC(i,j) == nodata_i4 .and. &
537 2266 : SMIc(i,j,t) == klu ) then
538 1157 : iC(i,j) = iC(krow, kcol)
539 : end if
540 : end do
541 : end do
542 : end do
543 :
544 : ! consolidate clusters
545 71 : if ( nCluster > 0 ) then
546 284 : allocate ( cno(nCluster), vec(nCluster) )
547 474 : cno = (/(i, i=1,nCluster)/)
548 158 : vec = -9
549 : nClusterR = nCluster
550 :
551 1633 : do k = nCells, 1, -1
552 1562 : krow = cellCoor(k,1)
553 1562 : kcol = cellCoor(k,2)
554 1562 : klu = iC(krow, kcol)
555 1562 : iul = max(krow-1,1)
556 1562 : idr = min(krow+1,nrows)
557 1562 : jul = kcol
558 1562 : jdr = min(kcol+1,ncols)
559 4473 : do j=jul, jdr
560 12425 : do i= iul, idr
561 8023 : klu1 = iC(i, j)
562 : if ( klu /= nodata_i4 .and. &
563 8023 : klu1 /= nodata_i4 .and. &
564 2840 : klu /= klu1 ) then
565 7 : cno(klu1) = -9
566 7 : nClusterR = nClusterR - 1
567 287 : where ( iC == klu1 ) iC = klu
568 : end if
569 : end do
570 : end do
571 : end do
572 : !
573 : ! delete small clusters < thesh. area
574 158 : do i=1, nCluster
575 87 : if (cno(i) == -9 ) cycle
576 3280 : ncc = count(iC == cno(i))
577 151 : if ( ncc <= thCellClus ) then
578 451 : where (iC == cno(i)) iC = nodata_i4
579 11 : cno(i) = -9
580 11 : nClusterR = nClusterR -1
581 : end if
582 : end do
583 : !
584 : ! reordering
585 403 : cno = pack (cno, mask = cno > 0, vector=vec)
586 158 : where (cno <= 0) cno = 0
587 : !
588 : !
589 211 : allocate (nCxCluster(nClusterR))
590 140 : do i=1, nClusterR
591 2829 : where (iC == cno(i)) iC = i
592 2900 : nCxCluster(i) = count(iC == i)
593 : end do
594 : !
595 71 : deallocate (cno, nCxCluster)
596 : end if
597 71 : nCluster = nClusterR
598 71 : end subroutine findClusters
599 :
600 : END MODULE mo_smi_drought_evaluation
|
