Line data Source code
1 : !> \file EDK_driver.f90
2 : !> \brief \copybrief ed_kriging
3 : !> \details \copydetails ed_kriging
4 :
5 : !> \brief External drift kriging - EDK
6 : !> \details Perform EDK (daily)
7 : !! Store results suitable for mHM
8 : !! special version to interpolate values in blocks
9 : !! *** can be parallelized ***
10 : !> \author Luis Samaniego
11 : !> \date 21.03.2006
12 : !> \date 11.06.2010
13 : !! - blocks, whole Germany
14 : !> \date 04.02.2012
15 : !! - changed to general edk version (excluded block seperation)
16 : !> \copyright Copyright 2005-\today, the CHS Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
17 : !! EDK is released under the LGPLv3+ license \license_note
18 2 : program ED_Kriging
19 :
20 2 : use mo_kind , only: i4, dp, sp
21 : use mo_edk_print_message , only: print_start_message, print_end_message
22 : use mo_julian , only: NDAYS, NDYIN, dec2date, julday
23 : use runControl , only: flagEDK, interMth, & ! flag for activate kriging, flag for 'OK' or 'EDK'
24 : correctNeg, & ! pre or temp
25 : distZero, &
26 : flagVario ! flag for activate variogram estimation
27 : use mainVar , only: yStart, yEnd, jStart, jEnd, tBuffer, nSta, DEMNcFlag, & ! interpolation time periods
28 : grid, gridMeteo, & ! grid properties of input and output grid
29 : nCell, MetSta, &
30 : noDataValue
31 : use kriging , only: edk_dist, cell
32 : use mo_edk_setvario , only: setVario, dMatrix
33 : use mo_netcdf , only: NcDataset, NcVariable
34 : use mo_edk_write , only: open_netcdf
35 : use mo_message , only: message
36 : use mo_edk , only: EDK, clean, WriteDataMeteo
37 : use mo_edk_read_data , only: readData
38 : use NetCDFVar , only: invert_y, variable_name
39 : USE mo_timer, ONLY : &
40 : timers_init, timer_start, timer_stop, timer_get ! Timing of processes
41 : use mo_string_utils, ONLY : num2str
42 : use mo_edk_cli, only: parse_command_line
43 : !$ use omp_lib, ONLY : OMP_GET_NUM_THREADS ! OpenMP routines
44 : implicit none
45 :
46 : character(256) :: fname
47 : character(256) :: author_name
48 : character(256) :: vname_data
49 : integer(i4) :: i, j, k, t
50 : integer(i4) :: icell ! loop varaible for cells
51 : integer(i4) :: jDay ! loop variable - current julian day
52 : integer(i4) :: itimer
53 : integer(i4) :: doy ! day of year
54 : integer(i4) :: year, month, day ! current date
55 : integer(i4) :: itime, itemp,sttemp,cnttemp ! loop variable for chunking write
56 : integer(i4) :: jstarttmp, jendtmp ! more loop variables
57 : integer(i4) :: n_threads ! OpenMP number of parallel threads
58 : integer(i4) :: ncell_thread
59 : integer(i4) :: iThread
60 : integer(i4) :: loop_factor
61 2 : real(sp), allocatable :: tmp_array(:, :, :) ! temporal array for output
62 2 : real(sp), allocatable :: tmp_time(:) ! temporal array for time output
63 8 : real(dp) :: param(3) ! variogram parameters
64 : type(NcDataset) :: nc_out
65 : type(NcVariable) :: nc_data, nc_time
66 :
67 2 : call parse_command_line()
68 2 : call print_start_message()
69 :
70 2 : loop_factor = 10 ! factor for setting openMP loop size
71 :
72 2 : n_threads = 1
73 : !$OMP PARALLEL
74 : !$ n_threads = OMP_GET_NUM_THREADS()
75 : !$OMP END PARALLEL
76 : !$ print *, 'Run with OpenMP with ', trim(num2str(n_threads)), ' threads.'
77 :
78 : ! initialize timers
79 2 : call timers_init
80 :
81 : ! start timer for reading
82 2 : itimer = 1
83 2 : call timer_start(itimer)
84 :
85 2 : call message('')
86 2 : call message(' >>> Reading data')
87 2 : call message('')
88 2 : call ReadData
89 2 : call timer_stop(itimer)
90 2 : call message('')
91 2 : call message(' in ', trim(num2str(timer_get(itimer), '(F9.3)')), ' seconds.')
92 2 : call message('')
93 : ! number of cells per thread
94 2 : ncell_thread = ceiling(real(nCell, sp) / real(loop_factor * n_threads, sp))
95 :
96 : ! sanity check
97 2 : if (variable_name == "pre") then
98 2 : if ( .not. correctNeg ) call message( &
99 : '***WARNING: you are probably interpolating precipitation data, but correctNeg is set .False.', &
100 0 : 'It should be set to .True. Check namelist!')
101 2 : if ( .not. distZero ) call message( &
102 : '***WARNING: you are probably interpolating precipitation data, but distZero is set .False. ', &
103 0 : 'It should be set to .True. Check namelist!')
104 : end if
105 :
106 2 : itimer = 2
107 2 : call timer_start(itimer)
108 2 : call message(' >>> Calculating distances and neighborhoods')
109 2 : call message('')
110 : ! call distance matrix
111 2 : call dMatrix
112 2 : call timer_stop(itimer)
113 2 : call message('')
114 2 : call message(' in ', trim(num2str(timer_get(itimer), '(F9.3)')), ' seconds.')
115 2 : call message('')
116 :
117 2 : itimer = 3
118 2 : call timer_start(itimer)
119 2 : call message(' >>> Estimate variogram')
120 2 : call message('')
121 : ! estimate variogram
122 2 : call setVario(param)
123 : ! write variogram
124 2 : if (flagVario) call WriteDataMeteo
125 2 : call timer_stop(itimer)
126 2 : call message('')
127 2 : call message(' in ', trim(num2str(timer_get(itimer), '(F9.3)')), ' seconds.')
128 2 : call message('')
129 :
130 : !write(*,*), "jStart = ",jStart
131 :
132 2 : if (interMth .gt. 0) then
133 2 : itimer = 4
134 2 : call timer_start(itimer)
135 2 : call message(' >>> Perform interpolation')
136 2 : call message('')
137 : ! open netcdf if necessary
138 2 : call open_netcdf(nc_out, nc_data, nc_time)
139 :
140 : ! do iCell = 1, nCell
141 : ! ! initialize cell
142 : ! allocate(cell(iCell)%Nk_old(nSta))
143 : ! cell(iCell)%Nk_old = nint(noDataValue)
144 : ! end do
145 :
146 2 : if (mod((jEnd - jStart + 1),tBuffer) .eq. 0) then ! just use mod
147 0 : iTime = ((jEnd - jStart + 1)/tBuffer)
148 : else
149 2 : iTime = ((jEnd - jStart + 1)/tBuffer) + 1
150 : end if
151 2 : write(*,*) "Total Number of Threads = ",n_threads
152 2 : write(*,*) "Total Number of Time Buffers = ",iTime
153 2 : t = 0
154 76 : bufferloop: do iTemp = 1, iTime
155 74 : write(*,*) " >>> Started buffer #", iTemp
156 74 : jStartTmp = jStart + (iTemp - 1) * tBuffer
157 74 : if (iTemp .lt. iTime) then
158 72 : jEndTmp = jStartTmp + tBuffer - 1
159 : else
160 2 : jEndTmp = jStartTmp + (jEnd-jStartTmp+1)
161 : end if ! use minimum to never exceed jEnd
162 74 : jEndTmp = min(jEndTmp, jEnd)
163 :
164 4514 : do iCell = 1, nCell
165 : ! initialize cell ! deallocate similarly
166 13320 : allocate(cell(iCell)%z(jStartTmp:jEndTmp))
167 48314 : cell(iCell)%z = real(noDataValue, sp)
168 : end do
169 :
170 : !print *, iTemp, iTime
171 :
172 : !$OMP parallel default(shared) &
173 : !$OMP private(iThread, iCell)
174 : !$OMP do SCHEDULE(dynamic)
175 814 : do iThread = 1, loop_factor * n_threads
176 : ! print *, 'thread: ', iThread, " start"
177 :
178 5254 : ncellsloop: do iCell = (iThread - 1) * ncell_thread + 1, min(iThread * ncell_thread, ncell)
179 :
180 : ! check DEM
181 4440 : if (nint(cell(iCell)%h) == grid%nodata_value ) then
182 19296 : cell(iCell)%z = gridMeteo%nodata_value
183 : cycle
184 : end if
185 : ! interploation
186 3404 : call EDK(iCell, jStartTmp, jEndTmp, edk_dist, MetSta, cell(iCell), doOK=(interMth==1))
187 : end do ncellsloop
188 :
189 : ! print *, 'thread: ', iThread, " end"
190 : end do
191 : !$OMP end do
192 : !$OMP end parallel
193 :
194 74 : if (DEMNcFlag == 1) then
195 : ! write output
196 0 : allocate(tmp_array(gridMeteo%nrows, gridMeteo%ncols, jEndTmp - jStartTmp + 1))
197 0 : allocate(tmp_time(jEndTmp - jStartTmp + 1))
198 :
199 0 : k = 0
200 0 : if (invert_y) then
201 0 : do i = 1, gridMeteo%ncols
202 0 : do j = 1, gridMeteo%nrows
203 0 : k = k + 1
204 0 : tmp_array(j,gridMeteo%ncols - i + 1,:) = cell(k)%z
205 : end do
206 : end do
207 : else
208 0 : do i = 1, gridMeteo%ncols
209 0 : do j = 1, gridMeteo%nrows
210 0 : k = k + 1
211 0 : tmp_array(j,i,:) = cell(k)%z
212 : end do
213 : end do
214 : end if
215 :
216 0 : do i = 1, jEndTmp - jStartTmp + 1
217 0 : tmp_time(i) = t
218 0 : t = t + 1
219 : end do
220 :
221 0 : sttemp = nint(tmp_time(1)+1)
222 0 : cnttemp = nint((tmp_time(size(tmp_time)) - sttemp))+2
223 :
224 : !write(*,*),"Final Output ",shape(tmp_array)
225 :
226 0 : call nc_time%setData(values=tmp_time,start=(/sttemp/),cnt=(/cnttemp/))
227 : !call nc_data%setData(values=tmp_array,start=(/1,1,sttemp/),cnt=(/size(tmp_array,1),size(tmp_array,2),cnttemp/))
228 0 : call nc_data%setData(values=tmp_array,start=(/1,1,sttemp/),cnt=(/size(tmp_array,1),size(tmp_array,2),cnttemp/))
229 :
230 : else
231 : ! write output
232 370 : allocate(tmp_array(gridMeteo%ncols, gridMeteo%nrows, jEndTmp - jStartTmp + 1))
233 222 : allocate(tmp_time(jEndTmp - jStartTmp + 1))
234 :
235 74 : k = 0
236 518 : do i = 1, gridMeteo%ncols
237 : ! do j = 1, gridMeteo%nrows
238 : ! k = k + 1
239 : ! tmp_array(i, gridMeteo%nrows - j + 1, :) = cell(k)%z
240 : ! end do
241 : ! end do
242 518 : if (invert_y) then
243 4884 : do j = gridMeteo%nrows, 1, -1
244 4440 : k = k + 1
245 48684 : tmp_array(i, gridMeteo%nrows - j + 1, :) = cell(k)%z
246 : end do
247 : else
248 0 : do j = 1, gridMeteo%nrows
249 0 : k = k + 1
250 0 : tmp_array(i, gridMeteo%nrows - j + 1, :) = cell(k)%z
251 : end do
252 : end if
253 : end do
254 : !t = 0
255 : !do i = 1, jEnd - jStart + 1
256 : ! tmp_time(i) = t
257 : ! t = t + 1
258 : !end do
259 :
260 804 : do i = 1, jEndTmp - jStartTmp + 1
261 730 : tmp_time(i) = t
262 804 : t = t + 1
263 : end do
264 :
265 : !write(*,*),tmp_time
266 74 : sttemp = nint(tmp_time(1)+1)
267 74 : cnttemp = nint((tmp_time(size(tmp_time)) - sttemp))+2
268 :
269 : !write(*,*),"Final Output ",shape(tmp_array)
270 :
271 222 : call nc_time%setData(values=tmp_time,start=(/sttemp/),cnt=(/cnttemp/))
272 : !call nc_data%setData(values=tmp_array,start=(/1,1,sttemp/),cnt=(/size(tmp_array,1),size(tmp_array,2),cnttemp/))
273 518 : call nc_data%setData(values=tmp_array,start=(/1,1,sttemp/),cnt=(/size(tmp_array,1),size(tmp_array,2),cnttemp/))
274 : end if
275 :
276 74 : deallocate(tmp_array, tmp_time)
277 : !deallocate(cell)
278 :
279 4516 : do iCell = 1, nCell
280 : ! initialize cell
281 4514 : deallocate(cell(iCell)%z)
282 : !cell(iCell)%z = noDataValue
283 : end do
284 :
285 : ! close netcdf if necessary
286 : !call nc_out%close() ! outside
287 :
288 : end do bufferloop
289 :
290 : ! close netcdf if necessary
291 2 : call nc_out%close()
292 :
293 2 : call timer_stop(itimer)
294 2 : call message('')
295 2 : call message(' in ', trim(num2str(timer_get(itimer), '(F9.3)')), ' seconds.')
296 2 : call message('')
297 : end if
298 : ! deallocate memory
299 2 : call clean
300 : !
301 2 : call print_end_message()
302 : !
303 2 : end program ED_Kriging
|
