Line data Source code
1 : !> \file main.f90
2 : !> \brief \copybrief sm_drought_index
3 : !> \details \copydetails sm_drought_index
4 :
5 : !> \brief SOIL MOISTURE DROUGHT INDEX
6 : !> \details Estimate the soil moisture index using daily values generated by mHM
7 : !! ALGORITHM:
8 : !! 1. Read netcdf files daily/montly
9 : !! 2. Estimate monthly values of mHM-SM(sum over all layers) for each grid cell
10 : !! - write them into netCDF
11 : !! 3. Estimate the empirical density function for each grid cell
12 : !! using a non-parametric approach, e.g. kernel smoother whose
13 : !! bandwith is estimated with an unbiased cross-validation criterium
14 : !! 4. Estimate the mHM-drought index
15 : !! q(s) = \int_a^s f(s)ds
16 : !! where:
17 : !! s = total soil moisture in mm over all soil layers in mHM
18 : !! a = lower limit of soil moisture a a given grid
19 : !! q(s)= quantile corresponding to s
20 : !! 5. Estimate the following indices
21 : !! - start
22 : !! - duration
23 : !! - magnitud
24 : !! - severity
25 : !! - affected area
26 : !! 6. Save results in netCDF format
27 : !!
28 : !> \authors Luis E. Samaniego-Eguiguren, Matthias Zink, Stephan Thober
29 : !> \date 15.02.2011
30 : !! - (LS) main structure
31 : !> \date 20.02.2011
32 : !! - (LS) debuging
33 : !> \date 25.05.2011
34 : !! - (LS) v3. scenarios
35 : !> \date 02.04.2012
36 : !! - (LS) v4. read COSMO SM
37 : !> \date 22.06.2012
38 : !! - (LS) v4. read WRF-NOAH SM
39 : !> \date 07.11.2016
40 : !! - (MZ) modularized version
41 : !> \date 20.03.2017
42 : !! - (LS) daily SMI, SAD flag, restructuring SAD
43 : !> \date 24.07.2018
44 : !! - (ST) bug fix in optimize width
45 : !> \copyright Copyright 2005-\today, the CHS Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
46 : !! SMI is released under the LGPLv3+ license \license_note
47 12 : program SM_Drought_Index
48 :
49 6 : use mo_kind, only : i4, sp, dp
50 : use mo_message, only : message
51 : use mo_smi_write, only : WriteNetCDF, &
52 : nDurations, durList, nClusters, &
53 : writeResultsCluster, WriteResultsBasins, &
54 : WriteSMI, WriteCDF
55 : use mo_smi_read, only : ReadDataMain
56 : use mo_smi, only : optimize_width, calSMI, invSMI
57 : use mo_smi_drought_evaluation, only : droughtIndicator, ClusterEvolution, ClusterStats, calSAD
58 : use mo_constants, only : nodata_sp
59 : use mo_smi_global_variables, only : period
60 : use mo_smi_cli, only : parse_command_line
61 : !$ use omp_lib, ONLY : OMP_GET_NUM_THREADS ! OpenMP routines
62 :
63 : implicit none
64 :
65 : ! variables
66 : logical :: do_cluster ! flag indicating whether cluster should be calculated
67 : logical :: do_sad ! flag indicating whether SAD analysis should be done
68 : logical :: ext_smi ! flag indicating to read external data for clustering
69 : logical :: invert_SMI ! flag for inverting SMI
70 : ! ! calculated or read from file
71 : logical :: read_opt_h ! read kernel width from file
72 : logical :: silverman_h ! flag indicating whether kernel width
73 : ! ! should be optimized
74 : logical :: do_basin ! do_basin flag
75 6 : logical, dimension(:,:), allocatable :: mask
76 :
77 : integer(i4) :: ii
78 6 : type(period) :: per_kde, per_eval, per_smi
79 : integer(i4) :: nCells ! number of effective cells
80 : integer(i4) :: d
81 :
82 : integer(i4) :: nCalendarStepsYear ! Number of calendar time steps per year (month=12, day=365)
83 :
84 : integer(i4) :: thCellClus ! treshold for cluster formation in space ~ 640 km2
85 : integer(i4) :: nCellInter ! number cells for joining clusters in time ~ 6400 km2
86 : integer(i4) :: deltaArea ! number of cells per area interval
87 6 : integer(i4), dimension(:,:), allocatable :: Basin_Id ! IDs for basinwise drought analysis
88 6 : integer(i4), dimension(:,:), allocatable :: cellCoor !
89 :
90 : real(sp) :: SMI_thld ! SMI threshold for clustering
91 : real(sp) :: cellsize ! cell edge lenght of input data
92 6 : real(sp), dimension(:,:), allocatable :: SM_kde ! monthly fields packed for estimation
93 6 : real(sp), dimension(:,:), allocatable :: SM_eval ! monthly fields packed for evaluation
94 6 : real(sp), dimension(:,:), allocatable :: SM_invert ! inverted monthly fields packed
95 6 : real(sp), dimension(:,:), allocatable :: SMI ! soil moisture index at evaluation array
96 6 : real(sp), dimension(:,:,:), allocatable :: dummy_d3_sp
97 6 : integer(i4), dimension(:,:,:), allocatable :: SMIc ! Drought indicator 1 - is under drought
98 : ! ! 0 - no drought
99 6 : real(dp), dimension(:,:), allocatable :: opt_h ! optimized kernel width field
100 6 : real(dp), dimension(:), allocatable :: lats_1d, lons_1d ! latitude and longitude vectors of input
101 6 : real(dp), dimension(:,:), allocatable :: lats_2d, lons_2d ! latitude and longitude vectors of input
102 6 : real(dp), dimension(:), allocatable :: easting ! easting coordinates of input
103 6 : real(dp), dimension(:), allocatable :: northing ! easting coordinates of input
104 :
105 :
106 : ! file handling
107 : character(256) :: outpath ! output path for results
108 :
109 : !parse CLI
110 6 : call parse_command_line()
111 :
112 : !$OMP PARALLEL
113 : !$ ii = OMP_GET_NUM_THREADS()
114 : !$OMP END PARALLEL
115 : !$ print *, 'Run with OpenMP with ', ii, ' threads.'
116 :
117 : call ReadDataMain( SMI, do_cluster, ext_smi, invert_smi, &
118 : read_opt_h, silverman_h, opt_h, lats_1d, lons_1d, lats_2d, lons_2d, easting, northing,&
119 : do_basin, mask, SM_kde, SM_eval, Basin_Id, &
120 : SMI_thld, outpath, cellsize, thCellClus, nCellInter, &
121 6 : do_sad, deltaArea, nCalendarStepsYear, per_kde, per_eval, per_smi )
122 :
123 : ! initialize some variables
124 236 : nCells = count( mask ) ! number of effective cells
125 :
126 6 : call message('FINISHED READING')
127 :
128 : ! optimize kernel width
129 6 : if ( (.NOT. read_opt_h) .AND. (.NOT. ext_smi)) then
130 1 : call optimize_width( opt_h, silverman_h, SM_kde, nCalendarStepsYear, per_kde)
131 1 : call message('optimizing kernel width...ok')
132 : end if
133 :
134 : ! calculate SMI values for SM_eval
135 6 : if (.NOT. ext_smi) then
136 16 : allocate( SMI( size( SM_eval, 1 ), size( SM_eval, 2 ) ) )
137 20579 : SMI(:,:) = nodata_sp
138 4 : call calSMI( opt_h, SM_kde, SM_eval, nCalendarStepsYear, SMI, per_kde, per_eval )
139 4 : call message('calculating SMI... ok')
140 : end if
141 :
142 : ! invert SMI according to given cdf
143 6 : if (invert_smi) then
144 : ! testing with calculated SMI -> SM_invert == SM_kde
145 : call invSMI(SM_kde, opt_h, SMI, nCalendarStepsYear, per_kde, per_smi, SM_invert)
146 : ! write results to file
147 5 : allocate(dummy_D3_sp(size(mask, 1), size(mask, 2), size(SM_invert, 2)))
148 6 : do ii = 1, size(SM_invert, 2)
149 6 : dummy_D3_sp(:, :, ii) = unpack(SM_invert(:, ii), mask, nodata_sp)
150 : end do
151 1 : call WriteNetcdf(outpath, 2, per_smi, lats_1d, lons_1d, lats_2d, lons_2d, easting, northing, SM_invert=dummy_D3_sp)
152 2 : deallocate(dummy_D3_sp)
153 : end if
154 :
155 : ! write output
156 6 : if (.NOT. ext_smi) then
157 4 : call WriteSMI( outpath, SMI, mask, per_eval, lats_1d, lons_1d, lats_2d, lons_2d, easting, northing)
158 4 : call message('write SMI...ok')
159 : end if
160 :
161 6 : if ((.not. read_opt_h) .and. (.not. ext_smi)) then
162 : call WriteCDF( outpath, SM_kde, opt_h, mask, per_kde, nCalendarStepsYear, &
163 1 : lats_1d, lons_1d, lats_2d, lons_2d, easting, northing)
164 1 : call message('write cdf_info file...ok')
165 : end if
166 :
167 : ! calculate drought cluster
168 7 : if ( do_cluster ) then
169 : ! drought indicator
170 1 : call droughtIndicator( SMI, mask, SMI_thld, cellCoor, SMIc )
171 1 : call WriteNetCDF(outpath, 3, per_smi, lats_1d, lons_1d, lats_2d, lons_2d, easting, northing, SMIc=SMIc)
172 :
173 : ! cluster indentification
174 : call ClusterEvolution( SMIc, size( mask, 1), size( mask, 2 ), size(SMI, 2), &
175 1 : nCells, cellCoor, nCellInter, thCellClus)
176 1 : call WriteNetCDF(outpath, 4, per_smi, lats_1d, lons_1d, lats_2d, lons_2d, easting, northing)
177 :
178 : ! statistics
179 1 : call ClusterStats(SMI, mask, size( mask, 1), size( mask, 2 ), size(SMI, 2), nCells, SMI_thld )
180 :
181 : ! write results
182 1 : if (nClusters > 0) call writeResultsCluster(SMIc, outpath, 1, &
183 1 : per_smi%y_start, per_smi%y_end, size(SMI, 2), nCells, deltaArea, cellsize)
184 1 : call message('Cluster evolution ...ok')
185 : end if
186 :
187 : ! SAD analysis
188 6 : if ( do_sad ) then
189 0 : do d = 1, nDurations
190 0 : call calSAD(SMI, mask, d, size( mask, 1), size( mask, 2 ), size(SMI, 2), nCells, deltaArea, cellsize)
191 : ! write SAD for a given duration + percentiles
192 : call writeResultsCluster(SMIc, outpath, 2, per_smi%y_start, per_smi%y_end, size(SMI, 2), &
193 0 : nCells, deltaArea, cellsize, durList(d))
194 : call WriteNetCDF(outpath, 5, per_kde, lats_1d, lons_1d, lats_2d, lons_2d,&
195 0 : easting, northing, duration=durList(d))
196 : end do
197 : end if
198 :
199 : ! make basin averages
200 6 : if ( do_basin ) then
201 : ! write SMI average over major basins
202 0 : call message('calculate Basin Results ...')
203 0 : call WriteResultsBasins( outpath, SMI, mask, per_kde%y_start, per_kde%y_end, size( SM_kde, 2 ), Basin_Id )
204 : end if
205 :
206 : ! print statement for check_cases
207 6 : call message('SMI: finished!')
208 : !
209 6 : end program SM_Drought_Index
|
