Line data Source code
1 : !> \file mo_smi_write.f90
2 : !> \brief \copybrief mo_smi_write
3 : !> \details \copydetails mo_smi_write
4 :
5 : !> \brief Writing subroutines for SMI
6 : !> \author Luis Samaniego
7 : !> \date 09.02.2011
8 : !> \copyright Copyright 2005-\today, the CHS Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
9 : !! SMI is released under the LGPLv3+ license \license_note
10 : module mo_smi_write
11 :
12 : use mo_kind, only: i4, sp, dp
13 : use mo_smi_global_variables, only: period
14 :
15 : implicit none
16 :
17 : ! ! clusters
18 : integer(i4), dimension(:,:,:), allocatable :: idCluster
19 : integer(i4) :: nInterArea
20 : integer(i4) :: nEvents ! total number of drough events to estimate SAD
21 : integer(i4), dimension(:), allocatable :: shortCnoList ! consolidated cluster no. list
22 : integer(i4), dimension(:,:), allocatable :: eventId ! event identification number, cluster, month of occurence
23 : integer(i4) :: nClusters ! number od clusters
24 : integer(i4), dimension(:), allocatable :: eIdPerm ! permutation of the event Ids with ascending area
25 : !
26 : ! cluster statistics
27 : real(dp), dimension(:,:), allocatable :: DAreaEvol ! drought area evolution (fraction Germany)
28 : real(dp), dimension(:,:), allocatable :: DTMagEvol ! magnitud
29 : real(dp), dimension(:), allocatable :: aDD ! average (mean) drought duration space-time
30 : real(dp), dimension(:), allocatable :: aDA ! average (mean) drought area
31 : real(dp), dimension(:), allocatable :: TDM ! total drought magnitud
32 : real(dp), dimension(:,:,:), allocatable :: SAD ! severity area duration curves
33 : real(dp), dimension(:,:), allocatable :: SADperc ! percentilles of severity area duration curves
34 : integer(i4) :: nDsteps ! number of durations steps
35 : real(dp), dimension(:,:,:), allocatable :: severity ! severity for a given duration
36 : real(dp), dimension(:,:,:), allocatable :: dASevol ! evolution of the drought areas and severity
37 :
38 : integer(i4), parameter :: nDurations = 4 ! number of durations
39 : integer(i4), dimension(nDurations), parameter :: durList = (/3,6,9,12/) !(/3,6,9,12/) ! list of durations to evaluate
40 : integer(i4), parameter :: nQProp = 3 ! number of SAD percetiles for a given duration
41 : real(dp), dimension(nQProp), parameter :: QProp = (/90._dp, 96._dp, 98._dp /)
42 : ! percentiles corresponding
43 : ! to return periods of 10,25,50 years
44 :
45 : integer(i4) :: nLargerEvents = 600 ! n. largest drought events
46 : !
47 : ! Basin summary
48 : integer(i4), parameter :: nBasins = 6
49 :
50 : contains
51 :
52 : !> \brief for writting the SMI to nc file
53 : !> \author Stephan Thober
54 : !> \date 5.7.2014
55 : !> \date 9.8.2019
56 : !! - refactored to new period structure and leap day handling
57 4 : subroutine WriteSMI( outpath, SMI, mask, per, lats_1d, lons_1d, lats_2d, lons_2d, easting, northing)
58 :
59 : use mo_kind, only: i4, sp
60 : use mo_string_utils, only: num2str
61 : use mo_constants, only: nodata_sp, nodata_dp
62 : use mo_netcdf, only: NcDataset, NcVariable, NcDimension
63 :
64 : implicit none
65 :
66 : ! input variables
67 : character(len=*), intent(in) :: outpath ! ouutput path for results
68 :
69 : logical, dimension(:,:), intent(in) :: mask
70 : type(period), intent(in) :: per
71 : real(sp), dimension(:,:), intent(in) :: SMI
72 : real(dp), dimension(:), allocatable, intent(in) :: lats_1d, lons_1d ! latitude and longitude fields of input
73 : real(dp), dimension(:,:), allocatable, intent(in) :: lats_2d, lons_2d ! latitude and longitude fields of input
74 : real(dp), dimension(:), allocatable, intent(in) :: easting ! easting coordinates of input
75 : real(dp), dimension(:), allocatable, intent(in) :: northing ! easting coordinates of input
76 : ! local Variables
77 : type(NcDataset) :: nc_out
78 : type(NcVariable) :: nc_var
79 : type(NcDimension) :: nc_y, nc_x, nc_tim
80 : character(256) :: Fname
81 : integer(i4) :: ii ! month
82 :
83 4 : real(sp), dimension(:,:,:), allocatable :: dummy_D3_sp ! field real unpacked
84 : integer(i4) :: y
85 : integer(i4) :: x
86 :
87 : ! initialize dimension
88 4 : y = size(mask, 2)
89 4 : x = size(mask, 1)
90 4 : print *,'x: ', x
91 4 : print *,'y: ', y
92 : ! fname
93 4 : fName = trim(outpath) //'SMI.nc'
94 4 : nc_out = NcDataset(fName, 'w')
95 4 : if (allocated(lats_1d)) then
96 0 : nc_y = nc_out%setDimension("lat", y)
97 0 : nc_x = nc_out%setDimension("lon", x)
98 4 : else if (allocated(lats_2d)) then
99 3 : nc_y = nc_out%setDimension("northing", y)
100 3 : nc_x = nc_out%setDimension("easting", x)
101 : else
102 1 : nc_y = nc_out%setDimension("y", y)
103 1 : nc_x = nc_out%setDimension("x", x)
104 : end if
105 4 : nc_tim = nc_out%setDimension("time", -1)
106 :
107 : ! unpack estimated SMIp
108 20 : allocate( dummy_d3_sp( x, y, size( SMI, 2) ) )
109 1029 : do ii = 1, size( SMI, 2)
110 1029 : dummy_d3_sp( :, :, ii) = unpack ( SMI(:,ii), mask, nodata_sp )
111 : end do
112 :
113 : ! save SMI (same sequence as in 1)
114 : nc_var = nc_out%setVariable('SMI', "f32", &
115 16 : (/ nc_x, nc_y, nc_tim /))
116 4 : call nc_var%setData(dummy_d3_sp)
117 4 : call nc_var%setAttribute('long_name', 'soil moisture index')
118 4 : call nc_var%setAttribute('missing_value', nodata_sp)
119 4 : call nc_var%setAttribute('units', '-')
120 4 : deallocate( dummy_d3_sp )
121 :
122 : ! add lat and lon
123 4 : if (allocated(lats_1d)) then
124 0 : print *, 'lat1d'
125 0 : nc_var = nc_out%setVariable('lat', "f64", (/ nc_y /))
126 0 : call nc_var%setData(lats_1d)
127 0 : call nc_var%setAttribute('long_name', 'latitude')
128 0 : call nc_var%setAttribute('missing_value', nodata_dp)
129 0 : call nc_var%setAttribute('units', 'degrees_north')
130 4 : else if (allocated(lats_2d)) then
131 3 : print *, 'lats2d'
132 9 : nc_var = nc_out%setVariable('lat', "f64", (/ nc_x, nc_y /))
133 3 : call nc_var%setData(lats_2d)
134 3 : call nc_var%setAttribute('long_name', 'latitude')
135 3 : call nc_var%setAttribute('missing_value', nodata_dp)
136 3 : call nc_var%setAttribute('units', 'degrees_north')
137 :
138 : ! nc_var = nc_out%setVariable('northing', "f64", (/ nc_y/))
139 : ! call nc_var%setData(northing)
140 : ! call nc_var%setAttribute('long_name', 'northing')
141 : ! call nc_var%setAttribute('missing_value', nodata_dp)
142 : ! call nc_var%setAttribute('units', 'meters')
143 : end if
144 :
145 4 : if (allocated(lons_1d)) then
146 0 : nc_var = nc_out%setVariable('lon', "f64", (/ nc_x /))
147 0 : call nc_var%setData(lons_1d)
148 0 : call nc_var%setAttribute('long_name', 'longitude')
149 0 : call nc_var%setAttribute('missing_value', nodata_dp)
150 0 : call nc_var%setAttribute('units', 'degrees_east')
151 4 : else if (allocated(lons_2d)) then
152 9 : nc_var = nc_out%setVariable('lon', "f64", (/ nc_x, nc_y /))
153 3 : call nc_var%setData(lons_2d)
154 3 : call nc_var%setAttribute('long_name', 'longitude')
155 3 : call nc_var%setAttribute('missing_value', nodata_dp)
156 3 : call nc_var%setAttribute('units', 'degrees_east')
157 :
158 : ! nc_var = nc_out%setVariable('easting', "f64", (/ nc_x /))
159 : ! call nc_var%setData(easting)
160 : ! call nc_var%setAttribute('long_name', 'easting')
161 : ! call nc_var%setAttribute('missing_value', nodata_dp)
162 : ! call nc_var%setAttribute('units', 'meters')
163 :
164 :
165 : end if
166 :
167 : ! add time
168 8 : nc_var = nc_out%setVariable('time', "i32", (/ nc_tim /))
169 4 : call nc_var%setData(per%time_points)
170 4 : call nc_var%setAttribute('long_name', 'time')
171 : call nc_var%setAttribute('units', 'days since ' // &
172 : trim(num2str(per%y_start, '(i4)')) // '-' // &
173 : trim(num2str(per%m_start, '(i2.2)')) // '-' // &
174 4 : trim(num2str(per%d_start, '(i2.2)')) // ' 00:00:00' )
175 : ! close file
176 4 : call nc_out%close()
177 :
178 4 : end subroutine WriteSMI
179 :
180 : !> \brief for writting the CDF information (that is the kernel and associated soil moisture values) to nc file
181 : !> \author Stephan Thober
182 : !> \date 8.8.2019
183 1 : subroutine WriteCDF( outpath, SM, hh, mask, per, nCalendarStepsYear, lats_1d, lons_1d, lats_2d, lons_2d, easting, northing)
184 :
185 4 : use mo_kind, only: i4, sp
186 : use mo_message, only: message
187 : use mo_string_utils, only: num2str
188 : use mo_constants, only: nodata_sp, nodata_dp
189 : use mo_netcdf, only: NcDataset, NcVariable, NcDimension
190 :
191 : implicit none
192 :
193 : ! input variables
194 : character(len=*), intent(in) :: outpath ! ouutput path for results
195 :
196 : logical, dimension(:,:), intent(in) :: mask
197 : type(period), intent(in) :: per
198 : real(sp), dimension(:,:), intent(in) :: SM
199 : integer(i4), intent(in) :: nCalendarStepsYear
200 : real(dp), dimension(:), allocatable, intent(in) :: lats_1d, lons_1d ! latitude and longitude fields of input
201 : real(dp), dimension(:,:), allocatable, intent(in) :: lats_2d, lons_2d ! latitude and longitude fields of input
202 : real(dp), dimension(:), allocatable, intent(in) :: easting ! easting coordinates of input
203 : real(dp), dimension(:), allocatable, intent(in) :: northing ! easting coordinates of input
204 : real(dp), dimension(:,:), intent(in) :: hh
205 :
206 : ! local Variables
207 : type(NcDataset) :: nc_out
208 : type(NcVariable) :: nc_var
209 : type(NcDimension) :: nc_y, nc_x, nc_tim, nc_cal
210 : character(256) :: Fname
211 : integer(i4) :: mm ! month
212 :
213 1 : real(sp), dimension(:,:,:), allocatable :: dummy_D3_sp ! field real unpacked
214 1 : real(dp), dimension(:,:,:), allocatable :: dummy_D3_dp ! field real unpacked
215 : integer(i4) :: y
216 : integer(i4) :: x
217 :
218 : ! initialize dimension
219 1 : y = size(mask, 2)
220 1 : x = size(mask, 1)
221 :
222 : ! fname
223 1 : fName = trim(outpath) //'cdf_info.nc'
224 1 : nc_out = NcDataset(fName, 'w')
225 1 : if (allocated(lats_1d)) then
226 0 : nc_y = nc_out%setDimension("lat", y)
227 0 : nc_x = nc_out%setDimension("lon", x)
228 1 : else if (allocated(lats_2d)) then
229 1 : nc_y = nc_out%setDimension("northing", y)
230 1 : nc_x = nc_out%setDimension("easting", x)
231 : else
232 0 : nc_y = nc_out%setDimension("y", y)
233 0 : nc_x = nc_out%setDimension("x", x)
234 : end if
235 1 : nc_tim = nc_out%setDimension("time", -1)
236 :
237 : ! unpack soil moisture estimated
238 5 : allocate( dummy_d3_sp( x, y, size( SM, 2) ) )
239 361 : do mm = 1, size( SM, 2)
240 361 : dummy_d3_sp( :, :, mm) = unpack ( SM(:,mm), mask, nodata_sp )
241 : end do
242 :
243 : ! save SM (same sequence as in 1)
244 : nc_var = nc_out%setVariable('SM', "f32", &
245 4 : (/ nc_x, nc_y, nc_tim /))
246 1 : call nc_var%setData(dummy_d3_sp)
247 1 : call nc_var%setAttribute('long_name', 'soil moisture saturation')
248 1 : call nc_var%setAttribute('missing_value', nodata_sp)
249 1 : call nc_var%setAttribute('units', 'fraction of pore space')
250 1 : deallocate( dummy_d3_sp )
251 :
252 : ! write out kernel width if it has been optimised
253 5 : allocate( dummy_D3_dp( x, y, size( hh, 2 ) ) )
254 13 : do mm = 1, size( hh, 2 )
255 13 : dummy_D3_dp(:, :, mm) = unpack( hh(:,mm), mask, real( nodata_sp, dp ) )
256 : end do
257 1 : nc_cal = nc_out%setDimension("time_steps", size(dummy_d3_dp,3))
258 : nc_var = nc_out%setVariable('kernel_width', "f64", &
259 4 : (/ nc_x, nc_y, nc_cal /))
260 1 : call nc_var%setData(dummy_d3_dp)
261 1 : call nc_var%setAttribute('long_name', 'optimised kernel width')
262 1 : call nc_var%setAttribute('missing_value', nodata_dp)
263 1 : if (nCalendarStepsYear .eq. 12_i4) then
264 1 : call nc_var%setAttribute('units', 'months')
265 0 : else if (nCalendarStepsYear .eq. 365_i4) then
266 0 : call nc_var%setAttribute('units', 'days')
267 : else
268 0 : call message("***ERROR: nCalendarStepsYear has to be 12 or 365 in subroutine WriteCDF")
269 0 : stop 1
270 : end if
271 :
272 1 : deallocate( dummy_D3_dp )
273 :
274 : ! add lat and lon
275 1 : if (allocated(lats_1d)) then
276 0 : nc_var = nc_out%setVariable('lat', "f64", (/ nc_y /))
277 0 : call nc_var%setData(lats_1d)
278 0 : call nc_var%setAttribute('long_name', 'latitude')
279 0 : call nc_var%setAttribute('missing_value', nodata_dp)
280 0 : call nc_var%setAttribute('units', 'degrees_north')
281 1 : else if (allocated(lats_2d)) then
282 3 : nc_var = nc_out%setVariable('lat', "f64", (/ nc_x, nc_y /))
283 1 : call nc_var%setData(lats_2d)
284 1 : call nc_var%setAttribute('long_name', 'latitude')
285 1 : call nc_var%setAttribute('missing_value', nodata_dp)
286 1 : call nc_var%setAttribute('units', 'degrees_north')
287 :
288 : ! nc_var = nc_out%setVariable('northing', "f64", (/ nc_y /))
289 : ! call nc_var%setData(northing)
290 : ! call nc_var%setAttribute('long_name', 'northing')
291 : ! call nc_var%setAttribute('missing_value', nodata_dp)
292 : ! call nc_var%setAttribute('units', 'meters')
293 : end if
294 :
295 1 : if (allocated(lons_1d)) then
296 0 : nc_var = nc_out%setVariable('lon', "f64", (/ nc_x /))
297 0 : call nc_var%setData(lons_1d)
298 0 : call nc_var%setAttribute('long_name', 'longitude')
299 0 : call nc_var%setAttribute('missing_value', nodata_dp)
300 0 : call nc_var%setAttribute('units', 'degrees_east')
301 :
302 1 : else if (allocated(lons_2d)) then
303 3 : nc_var = nc_out%setVariable('lon', "f64", (/ nc_x, nc_y /))
304 1 : call nc_var%setData(lons_2d)
305 1 : call nc_var%setAttribute('long_name', 'longitude')
306 1 : call nc_var%setAttribute('missing_value', nodata_dp)
307 1 : call nc_var%setAttribute('units', 'degrees_east')
308 :
309 : ! nc_var = nc_out%setVariable('easting', "f64", (/ nc_x /))
310 : ! call nc_var%setData(easting)
311 : ! call nc_var%setAttribute('long_name', 'easting')
312 : ! call nc_var%setAttribute('missing_value', nodata_dp)
313 : ! call nc_var%setAttribute('units', 'meters')
314 : end if
315 :
316 : ! add time
317 2 : nc_var = nc_out%setVariable('time', "i32", (/ nc_tim /))
318 1 : call nc_var%setData(per%time_points)
319 1 : call nc_var%setAttribute('long_name', 'time')
320 : call nc_var%setAttribute('units', 'days since ' // &
321 : trim(num2str(per%y_start, '(i4)')) // '-' // &
322 : trim(num2str(per%m_start, '(i2.2)')) // '-' // &
323 1 : trim(num2str(per%d_start, '(i2.2)')) // ' 00:00:00' )
324 : ! close file
325 1 : call nc_out%close()
326 :
327 1 : end subroutine WriteCDF
328 :
329 : !> \brief WRITE netCDF files
330 : !! FORMAT netCDF
331 : !! http://www.unidata.ucar.edu/software/netcdf/
332 : !> \author Luis E. Samaniego-Eguiguren, UFZ
333 : !> \date 16.02.2011
334 3 : subroutine WriteNetCDF(outpath, wFlag, per, lats_1d, lons_1d, lats_2d, lons_2d, easting, northing, &
335 3 : SMIc, SM_invert, duration)
336 : !
337 1 : use mo_kind, only: i4
338 : use mo_string_utils, only: num2str
339 : use mo_constants, only: nodata_i4 , nodata_dp, nodata_sp
340 : use mo_netcdf, only: NcDataset, NcVariable, NcDimension
341 :
342 : implicit none
343 : !
344 : ! input variables
345 : character(len=*), intent(in) :: outpath ! ouutput path for results
346 : type(period), intent(in) :: per
347 : integer(i4), intent(in) :: wFlag
348 : real(dp), dimension(:), allocatable, intent(in) :: lats_1d, lons_1d ! latitude and longitude fields of input
349 : real(dp), dimension(:,:), allocatable, intent(in) :: lats_2d, lons_2d ! latitude and longitude fields of input
350 : real(dp), dimension(:), allocatable, intent(in) :: easting ! easting coordinates of input
351 : real(dp), dimension(:), allocatable, intent(in) :: northing ! easting coordinates of input
352 : integer(i4), dimension(:,:,:), optional, intent(in) :: SMIc ! Drought indicator
353 : real(sp), dimension(:,:,:), optional, intent(in) :: SM_invert
354 : integer(i4), optional, intent(in) :: duration ! optional, duration
355 :
356 : ! local Variables
357 : type(NcDataset) :: nc_out
358 : type(NcVariable) :: nc_var
359 : type(NcDimension) :: nc_y, nc_x, nc_tim
360 : character(256) :: Fname
361 :
362 4 : select case (wFlag)
363 :
364 : case (2)
365 : ! fname
366 1 : fName = trim(outpath)//'SM_invert.nc'
367 1 : nc_out = NcDataset(fname, 'w')
368 :
369 1 : if (allocated(lats_1d)) then
370 0 : nc_y = nc_out%setDimension("lat", size(SM_invert, 2))
371 0 : nc_x = nc_out%setDimension("lon", size(SM_invert, 1))
372 1 : else if (allocated(lats_2d)) then
373 1 : nc_y = nc_out%setDimension("northing", size(SM_invert, 2))
374 1 : nc_x = nc_out%setDimension("easting", size(SM_invert, 1))
375 : else
376 0 : nc_y = nc_out%setDimension("y", size(SM_invert, 2))
377 0 : nc_x = nc_out%setDimension("x", size(SM_invert, 1))
378 : end if
379 1 : nc_tim = nc_out%setDimension("time", -1)
380 :
381 : nc_var = nc_out%setVariable('SM_Lall', "f32", &
382 4 : (/ nc_x, nc_y, nc_tim /))
383 1 : call nc_var%setData(SM_invert)
384 1 : call nc_var%setAttribute('long_name', 'SM according to inverse of SMI')
385 1 : call nc_var%setAttribute('missing_value', nodata_sp)
386 1 : call nc_var%setAttribute('units', 'mm')
387 :
388 :
389 : case (3)
390 : ! fname
391 1 : fName = trim(outpath)//'SMIc.nc'
392 1 : nc_out = NcDataset(fname, 'w')
393 1 : print *, 'y SMIc: ',size(SMIc, 2)
394 1 : print *, 'x SMIc: ',size(SMIc, 1)
395 1 : if (allocated(lats_1d)) then
396 0 : nc_y = nc_out%setDimension("lat", size(SMIc, 2))
397 0 : nc_x = nc_out%setDimension("lon", size(SMIc, 1))
398 1 : else if (allocated(lats_2d)) then
399 1 : nc_y = nc_out%setDimension("northing", size(SMIc, 2))
400 1 : nc_x = nc_out%setDimension("easting", size(SMIc, 1))
401 : else
402 0 : nc_y = nc_out%setDimension("y", size(SMIc, 2))
403 0 : nc_x = nc_out%setDimension("x", size(SMIc, 1))
404 : end if
405 1 : nc_tim = nc_out%setDimension("time", -1)
406 :
407 : nc_var = nc_out%setVariable('mSMIc', "i32", &
408 4 : (/ nc_x, nc_y, nc_tim /))
409 1 : call nc_var%setData(SMIc)
410 1 : call nc_var%setAttribute('long_name', 'monthly SMI indicator SMI < th')
411 1 : call nc_var%setAttribute('missing_value', nodata_i4)
412 1 : call nc_var%setAttribute('units', '-')
413 :
414 : case (4)
415 : ! fname
416 1 : fName = trim(outpath)//'DCluster.nc'
417 1 : nc_out = NcDataset(fname, 'w')
418 :
419 1 : if (allocated(lats_1d)) then
420 0 : nc_y = nc_out%setDimension("lat", size(idCluster, 2))
421 0 : nc_x = nc_out%setDimension("lon", size(idCluster, 1))
422 1 : else if (allocated(lats_2d)) then
423 1 : nc_y = nc_out%setDimension("northing", size(idCluster, 2))
424 1 : nc_x = nc_out%setDimension("easting", size(idCluster, 1))
425 : else
426 0 : nc_y = nc_out%setDimension("y", size(idCluster, 2))
427 0 : nc_x = nc_out%setDimension("x", size(idCluster, 1))
428 : end if
429 1 : nc_tim = nc_out%setDimension("time", -1)
430 :
431 : nc_var = nc_out%setVariable('mDC', "i32", &
432 4 : (/ nc_x, nc_y, nc_tim /))
433 1 : call nc_var%setData(idCluster)
434 1 : call nc_var%setAttribute('long_name', 'consolidated cluster evolution')
435 1 : call nc_var%setAttribute('missing_value', nodata_i4)
436 1 : call nc_var%setAttribute('units', '-')
437 :
438 : case (5)
439 : ! fname
440 0 : write (fName, '(i2.2)') duration
441 0 : fName = trim(outpath)//'sev_'//trim(fName)//'.nc'
442 0 : nc_out = NcDataset(fname, 'w')
443 :
444 0 : if (allocated(lats_1d)) then
445 0 : nc_y = nc_out%setDimension("lat", size(severity, 2))
446 0 : nc_x = nc_out%setDimension("lon", size(severity, 1))
447 0 : else if (allocated(lats_2d)) then
448 0 : nc_y = nc_out%setDimension("northing", size(severity, 2))
449 0 : nc_x = nc_out%setDimension("easting", size(severity, 1))
450 : else
451 0 : nc_y = nc_out%setDimension("y", size(severity, 2))
452 0 : nc_x = nc_out%setDimension("x", size(severity, 1))
453 : end if
454 0 : nc_tim = nc_out%setDimension("time", -1)
455 :
456 : nc_var = nc_out%setVariable('Severity', "f64", &
457 0 : (/ nc_x, nc_y, nc_tim /))
458 0 : call nc_var%setData(idCluster)
459 0 : call nc_var%setAttribute('long_name', 'd-month severity')
460 0 : call nc_var%setAttribute('missing_value', nodata_dp)
461 3 : call nc_var%setAttribute('units', '-')
462 :
463 : end select
464 :
465 : ! add lat and lon
466 3 : if (allocated(lats_1d)) then
467 0 : nc_var = nc_out%setVariable('lat', "f64", (/ nc_y /))
468 0 : call nc_var%setData(lats_1d)
469 0 : call nc_var%setAttribute('long_name', 'latitude')
470 0 : call nc_var%setAttribute('missing_value', nodata_dp)
471 0 : call nc_var%setAttribute('units', 'degrees_north')
472 3 : else if (allocated(lats_2d)) then
473 9 : nc_var = nc_out%setVariable('lat', "f64", (/ nc_x, nc_y /))
474 3 : call nc_var%setData(lats_2d)
475 3 : call nc_var%setAttribute('long_name', 'latitude')
476 3 : call nc_var%setAttribute('missing_value', nodata_dp)
477 3 : call nc_var%setAttribute('units', 'degrees_north')
478 :
479 : ! nc_var = nc_out%setVariable('northing', "f64", (/ nc_y /))
480 : ! call nc_var%setData(northing)
481 : ! call nc_var%setAttribute('long_name', 'northing')
482 : ! call nc_var%setAttribute('missing_value', nodata_dp)
483 : ! call nc_var%setAttribute('units', 'meters')
484 : end if
485 :
486 3 : if (allocated(lons_1d)) then
487 0 : nc_var = nc_out%setVariable('lon', "f64", (/ nc_x /))
488 0 : call nc_var%setData(lons_1d)
489 0 : call nc_var%setAttribute('long_name', 'longitude')
490 0 : call nc_var%setAttribute('missing_value', nodata_dp)
491 0 : call nc_var%setAttribute('units', 'degrees_east')
492 3 : else if (allocated(lons_2d)) then
493 9 : nc_var = nc_out%setVariable('lon', "f64", (/ nc_x, nc_y /))
494 3 : call nc_var%setData(lons_2d)
495 3 : call nc_var%setAttribute('long_name', 'longitude')
496 3 : call nc_var%setAttribute('missing_value', nodata_dp)
497 3 : call nc_var%setAttribute('units', 'degrees_east')
498 :
499 : ! nc_var = nc_out%setVariable('easting', "f64", (/ nc_x /))
500 : ! call nc_var%setData(easting)
501 : ! call nc_var%setAttribute('long_name', 'easting')
502 : ! call nc_var%setAttribute('missing_value', nodata_dp)
503 : ! call nc_var%setAttribute('units', 'meters')
504 : end if
505 :
506 : ! add time
507 : if ((wflag .eq. 2) .or. &
508 3 : (wflag .eq. 3) .or. &
509 : (wflag .eq. 4)) then
510 6 : nc_var = nc_out%setVariable('time', "i32", (/ nc_tim /))
511 3 : call nc_var%setData(per%time_points)
512 3 : call nc_var%setAttribute('long_name', 'time')
513 : call nc_var%setAttribute('units', 'days since ' // &
514 : trim(num2str(per%y_start, '(i4)')) // '-' // &
515 : trim(num2str(per%m_start, '(i2.2)')) // '-' // &
516 3 : trim(num2str(per%d_start, '(i2.2)')) // ' 00:00:00' )
517 : end if
518 : ! close file
519 3 : call nc_out%close()
520 3 : end subroutine WriteNetCDF
521 :
522 : !> \brief WRITE Results of the cluster analysis
523 : !> \author Luis E. Samaniego-Eguiguren, UFZ
524 : !> \date 17.03.2011
525 1 : subroutine WriteResultsCluster(SMIc, outpath, wFlag, yStart, yEnd, nMonths, nCells, &
526 : deltaArea, cellsize, d)
527 :
528 3 : use mo_kind, only : i4, dp
529 : use mo_constants, only : YearMonths
530 :
531 : !
532 : implicit none
533 :
534 : ! input variables
535 : character(len=*), intent(in) :: outpath ! ouutput path for results
536 : integer(i4), dimension(:,:,:),intent(in) :: SMIc ! Drought indicator 1 - is under drought
537 : integer(i4), intent(in) :: wFlag
538 : integer(i4), intent(in) :: yStart
539 : integer(i4), intent(in) :: yEnd
540 : integer(i4), intent(in) :: nMonths
541 : integer(i4), intent(in) :: nCells
542 : integer(i4), intent(in) :: deltaArea
543 : real(sp), intent(in) :: cellsize
544 : integer(i4), optional, intent(in) :: d
545 : real(dp), parameter :: eps = 1.0e-5_dp ! EPSILON(1.0_dp)
546 :
547 :
548 : ! local variables
549 : real(dp) :: pDArea
550 : !
551 : integer(i4) :: i, j, t, k, y, m, mStart, mEnd
552 : character(len=256) :: FMT
553 : character(len=256) :: fName
554 :
555 2 : select case (wFlag)
556 : case (1)
557 : ! main statistics
558 1 : fName = trim(outpath) // 'results_ADM.txt'
559 1 : open (10, file = fName, status='unknown')
560 1 : write (10, 100 ) 'i', 'c_Id', 'mStart', 'mEnd', 'aDD', 'aDA', 'TDM'
561 70 : do i = 1,nClusters
562 : ! find starting and ending months of every cluster
563 69 : mStart = 0
564 69 : mEnd = 0
565 13090 : do t = 1, nMonths
566 13090 : if ( DAreaEvol(t,i) .gt. 0.0_dp) then
567 69 : mEnd = t
568 69 : if (mStart .eq. 0) mStart = t
569 : end if
570 13090 : if ( ( DAreaEvol(t,i) .lt. eps) .and. (mStart .gt. 0) ) exit
571 : end do
572 70 : write (10,110) i, shortCnoList(i), mStart, mEnd, aDD(i), aDA(i), TDM(i)
573 : end do
574 1 : close (10)
575 : !
576 1 : fName = trim(outpath) // 'DArea_evol.txt'
577 1 : open (11, file = fName, status='unknown')
578 1 : write (FMT, 120) '(a5,', nClusters, '(2x,a2,i7.7))'
579 70 : write (11, FMT ) 'm', ('c_', shortCnoList(i), i=1,nClusters)
580 1 : write (FMT, 130) '(i5,', nClusters, 'e11.5)'
581 361 : do t=1, nMonths
582 361 : write (11,FMT) t, (DAreaEvol(t,i), i=1,nClusters)
583 : end do
584 1 : close (11)
585 : !
586 1 : fName = trim(outpath) // 'TDM_evol.txt'
587 1 : open (12, file = fName, status='unknown')
588 1 : write (FMT, 120) '(a5,', nClusters, '(2x,a2,i7.7))'
589 70 : write (12, FMT ) 'm', ('c_', shortCnoList(i), i=1,nClusters)
590 1 : write (FMT, 130) '(i5,', nClusters, 'e11.5)'
591 361 : do t=1, nMonths
592 361 : write (12,FMT) t, (DTMagEvol(t,i), i=1,nClusters)
593 : end do
594 1 : close (12)
595 : !
596 1 : fName = trim(outpath) // 'event_ids.txt'
597 1 : open (13, file = fName, status='unknown')
598 1 : write (13, 140 ) '<event>', 'c_Id', 'month', 'nCells'
599 : ! sorted events in ascending order of areal extend
600 70 : do i=1, nEvents
601 277 : write (13,145) eIdPerm(i), (eventId(eIdPerm(i),j), j=1,3)
602 : end do
603 1 : close (13)
604 : !
605 : ! NEW STATISTICS
606 : ! total drought area evolution (% w.r.t. whole domain)
607 1 : fName = trim(outpath) // 'DArea_evol_total.txt'
608 1 : open (14, file = fName, status='unknown')
609 1 : write (14, 150 ) 'year', 'month', '%AreaEU'
610 1 : t = 0
611 31 : do y =yStart, yEnd
612 391 : do m = 1, int(YearMonths, i4)
613 360 : t = t + 1
614 : pdArea = real( count( SMIc(:,:,t) == 1 ), dp) / &
615 14760 : real( nCells, dp) * 1e2_dp
616 390 : write(14,160) y, m, pdArea
617 : end do
618 : end do
619 1 : close (14)
620 : !
621 : ! Time evolution of the cluster area (less than SMIc) and monthly severity
622 1 : fName = trim(outpath) // 'DcArea_sev_evol.txt'
623 1 : open (15, file = fName, status='unknown')
624 1 : write (15, 170 ) 'year', 'month', '%cAreaEU','SevDE'
625 1 : t = 0
626 31 : do y =yStart, yEnd
627 391 : do m = 1, int(YearMonths, i4)
628 360 : t = t + 1
629 390 : write(15,180) y, m, dASevol(t,1,nBasins+1), dASevol(t,2,nBasins+1)
630 : end do
631 : end do
632 1 : close (15)
633 : !
634 : case(2)
635 : ! SAD curves for each event and duration
636 0 : do k = 1, nLargerEvents
637 0 : write (fName,200) 'SAD_e_', eIdPerm( nEvents + 1 - k ) , '_', d, '.txt'
638 0 : fName = trim(outpath) // trim(fName)
639 0 : open (20, file=fName, status='unknown')
640 0 : write (20,210) 'Area[km2]', 'Severity'
641 0 : write (20,220) (( SAD(i, j, k ), j=1,2 ), i=1, nInterArea)
642 0 : close(20)
643 : end do
644 :
645 : ! SAD percentiles for a given duration
646 0 : write (fName,230) 'SAD_perc_', d, '.txt'
647 0 : fName = trim(outpath) // trim(fName)
648 0 : open (21, file=fName, status='unknown')
649 0 : write (FMT, 240) '(a15,', nQProp, '(9x,a2,f4.2))'
650 0 : write (21,FMT) 'Area[km2]', ('p_', QProp(i), i=1,nQProp)
651 0 : write (FMT, 250) '(f15.0,', nQProp, 'f15.5)'
652 0 : write (21,FMT) ( real(i * deltaArea, dp)*(real(cellsize,dp)**2.0_dp), (SADperc(i,j), j=1,nQProp), i=1, nInterArea)
653 :
654 1 : close (21)
655 : end select
656 :
657 : 100 format (2a15, 2a10, 3a15)
658 : 110 format (2i15, 2i10, 3f15.5)
659 : 120 format (a4,i5,a13)
660 : 130 format (a4,i5,a6)
661 : 140 format (4a12)
662 : 145 format (4i12)
663 : 150 format (2a10,a10)
664 : 160 format (2i10,f10.3)
665 : 170 format (2a10,2a10)
666 : 180 format (2i10,2f10.3)
667 :
668 : 200 format (a6,i5.5,a,i2.2,a4)
669 : 210 format (2a15)
670 : 220 format (f15.0, f15.5)
671 : 230 format (a9,i5.5,a,i2.2,a4)
672 : 240 format (a5,i2,a13)
673 : 250 format (a7,i2,a6)
674 1 : end subroutine WriteResultsCluster
675 :
676 :
677 : !> \brief WRITE Results of the cluster SMI basins
678 : !> \author Luis E. Samaniego-Eguiguren, UFZ
679 : !> \date 24.05.2011
680 0 : subroutine WriteResultsBasins( outpath, SMI, mask, yStart, yEnd, nMonths, Basin_Id )
681 :
682 : use mo_kind, only : i4
683 : use mo_constants, only : nodata_dp, YearMonths
684 :
685 : implicit none
686 :
687 : ! input variables
688 : character(len=*), intent(in) :: outpath ! ouutput path for results
689 : real(sp), dimension(:,:), intent(in) :: SMI
690 : logical, dimension(:,:), intent(in) :: mask
691 : integer(i4), intent(in) :: yStart
692 : integer(i4), intent(in) :: yEnd
693 : integer(i4), intent(in) :: nMonths ! number of simulated months
694 : integer(i4), dimension(:,:), intent(in) :: Basin_Id ! IDs for basinwise drought analysis
695 :
696 : ! local variables
697 : real(dp), dimension(size(mask,1), &
698 0 : size(mask,2), nMonths) :: SMI_unpack
699 : character(len=256) :: fName
700 : integer(i4) :: i, m, y, j
701 0 : real(dp), dimension(:,:), allocatable :: Basin_SMI
702 :
703 : !-------------------------
704 : ! basin wise
705 : !-------------------------
706 0 : allocate( Basin_SMI (nMonths, nBasins+1) )
707 : ! unpack SMI
708 0 : do i = 1, size(SMI,2)
709 0 : SMI_unpack(:,:,i) = unpack(real(SMI(:,i), dp), mask, nodata_dp)
710 : end do
711 :
712 0 : Basin_SMI = nodata_dp
713 : !
714 0 : do m = 1, nMonths
715 0 : do i = 1, nBasins
716 : !
717 0 : Basin_SMI(m,i) = sum( SMI_unpack(:,:,m), Basin_Id(:,:) == i ) / real(count(Basin_Id(:,:) == i), dp)
718 : !
719 : end do
720 : Basin_SMI(m,nBasins+1) = sum(SMI_unpack(:,:,m), mask ) / &
721 0 : real(count(mask), dp)
722 : end do
723 : !
724 : ! Print Results
725 0 : fName = trim(outpath) // 'basin_avg_SMI.txt'
726 0 : open(20, file =fName , status = 'unknown')
727 0 : write(20, 1) 'Month_No', 'year', 'month', ('Basin_', i, i = 1, nBasins), 'Germany_Avg'
728 : !
729 0 : fName = trim(outpath) // 'basin_avg_dArea.txt'
730 0 : open(21, file =fName , status = 'unknown')
731 0 : write(21, 3) 'Month_No', 'year', 'month', ('Basin_', i, i = 1, nBasins)
732 : !
733 0 : fName = trim(outpath) // 'basin_avg_sev.txt'
734 0 : open(22, file =fName , status = 'unknown')
735 0 : write(22, 3) 'Month_No', 'year', 'month', ('Basin_', i, i = 1, nBasins)
736 : !! NEW !!
737 0 : m=0
738 0 : do y = yStart, yEnd
739 0 : do j = 1, int(YearMonths, i4)
740 0 : m = m + 1
741 0 : write(20, 2) m, y, j, (Basin_SMI(m,i), i = 1, nBasins), Basin_SMI(m,nBasins+1)
742 0 : write(21, 4) m, y, j, (dASevol(m,1,i), i = 1, nBasins)
743 0 : write(22, 4) m, y, j, (dASevol(m,2,i), i = 1, nBasins)
744 : end do
745 : end do
746 0 : close(20)
747 0 : close(21)
748 0 : close(22)
749 : !
750 : 1 format(3a8, 2x, 6(a6, i2.2, 2x), a13 )
751 : 2 format(3i8, 2x, 6(f8.4, 2x), f13.4)
752 : 3 format(3a8, 2x, 6(a6, i2.2, 2x) )
753 : 4 format(3i8, 2x, 6(f8.4, 2x) )
754 :
755 0 : end subroutine WriteResultsBasins
756 :
757 : end module mo_smi_write
|
