Line data Source code
1 : !> \file mo_nc_output.f90
2 : !> \brief \copybrief mo_nc_output
3 : !> \details \copydetails mo_nc_output
4 :
5 : !> \brief Creates NetCDF output for different fluxes and state variables of mHM.
6 : !> \details NetCDF is first initialized and later on variables are put to the NetCDF.
7 : !> \changelog
8 : !! - David Schaefer Aug 2015
9 : !! - major rewrite
10 : !! - Stephan Thober Oct 2015
11 : !! - adapted to mRM
12 : !! - O. Rakovec, R. Kumar Nov 2017
13 : !! - added project description for the netcdf outputs
14 : !! - S. Mueller, Dec 2022
15 : !! - unified module for mHM and mRM
16 : !> \authors Matthias Zink
17 : !> \date Apr 2013
18 : !> \copyright Copyright 2005-\today, the mHM Developers, Luis Samaniego, Sabine Attinger: All rights reserved.
19 : !! mHM is released under the LGPLv3+ license \license_note
20 : !> \ingroup f_common
21 : module mo_nc_output
22 :
23 : use mo_kind, only : i4, dp, sp
24 : use mo_common_types, only: Grid
25 : use mo_common_variables, only : project_details, setup_description, simulation_type, &
26 : Conventions, contact, mHM_details, history, dirOut, iFlag_cordinate_sys
27 : use mo_file, only : version
28 : use mo_common_constants, only : nodata_dp
29 : use mo_netcdf, only : NcDataset, NcDimension, NcVariable
30 : use mo_common_mHM_mRM_variables, only: timeStep
31 :
32 : implicit none
33 :
34 : public :: OutputDataset, OutputVariable, set_attributes, data_dims, data_dtype
35 :
36 : private
37 :
38 : type OutputVariable
39 : type(NcVariable) :: nc !< NcDataset which contains the variable
40 : logical :: avg = .false. !< average data before writing
41 : logical, pointer :: mask(:, :) !< mask to reconstruct data
42 : real(dp), allocatable :: data(:) !< store the data between writes
43 : integer(i4) :: counter = 0 !< count the number of updateVariable calls
44 :
45 : contains
46 : procedure, public :: updateVariable
47 : procedure, public :: writeVariableTimestep
48 :
49 : end type OutputVariable
50 :
51 : ! constructor interface
52 : interface OutputVariable
53 : procedure newOutputVariable
54 : end interface OutputVariable
55 :
56 : type OutputDataset
57 : integer(i4) :: iDomain !< domain id
58 : type(NcDataset) :: nc !< NcDataset to write
59 : type(OutputVariable), allocatable :: vars(:) !< store all created (dynamic) variables
60 : integer(i4) :: counter !< count written time steps
61 : integer(i4) :: previous_time !< previous time steps for bounds
62 : integer(i4) :: time_unit_factor !< possible factor to convert hours to minutes when using center as time reference
63 : integer(i4) :: outputs_frequence !< write out frequence (-3: yearly, -2: monthly, -1: daily, 0: end of run, >0: after n steps)
64 : integer(i4) :: time_reference !< time stamp reference (0: begin, 1: center, 2: end of time interval)
65 : logical :: double_precision !< output precision switch for nc files
66 :
67 : contains
68 : procedure, public :: writeTimestep
69 : procedure, public :: close
70 :
71 : end type OutputDataset
72 :
73 : ! constructor interface
74 : interface OutputDataset
75 : procedure newOutputDataset
76 : end interface OutputDataset
77 :
78 : contains
79 :
80 : !> \brief Output variable dtype for single or double precision.
81 : !> \return "f64" or "f32"
82 38 : character(3) function data_dtype(double_precision)
83 : implicit none
84 : logical, intent(in) :: double_precision !< flag to use double precision
85 38 : if ( double_precision ) then
86 38 : data_dtype = "f64"
87 : else
88 0 : data_dtype = "f32"
89 : end if
90 38 : end function data_dtype
91 :
92 : !> \brief Output variable dimension names.
93 : !> \return (X, Y, T) names tuple
94 76 : function data_dims()
95 : implicit none
96 : character(16), dimension(3) :: data_dims
97 19 : if (iFlag_cordinate_sys == 0) then
98 76 : data_dims = (/"easting ", "northing", "time "/) ! X & Y coordinate system
99 : else
100 0 : data_dims = (/"lon ", "lat ", "time"/) ! lat & lon coordinate system
101 : endif
102 19 : end function data_dims
103 :
104 : !> \brief Initialize OutputVariable
105 : !> \details Modifications:
106 : !! - David Schaefer Nov 2017 - added NcVariable initialization
107 : !! - Robert Schweppe Jun 2018 - refactoring and reformatting
108 : !> \return type(OutputVariable)
109 : !> \authors David Schaefer
110 : !> \date June 2015
111 59 : function newOutputVariable(nc, name, dtype, dims, ncells, mask, deflate_level, avg) result(out)
112 : implicit none
113 :
114 : type(NcDataset), intent(in) :: nc !< NcDataset which contains the variable
115 : character(*), intent(in) :: name !< name of the variable
116 : character(*), intent(in) :: dtype !< data type of the variable
117 : character(16), intent(in), dimension(3) :: dims !< dimensions of the variable (by name)
118 : integer(i4), intent(in) :: ncells !< number of cells in domain
119 : logical, intent(in), target, dimension(:, :) :: mask !< mask of the variable
120 : integer(i4), intent(in) :: deflate_level !< deflate level for compression
121 : logical, intent(in), optional :: avg !< flag to average the data before writing
122 :
123 : type(OutputVariable) :: out
124 :
125 177 : allocate(out%data(ncells))
126 59 : out%nc = nc%setVariable(name, dtype, dims, deflate_level = deflate_level, shuffle = .true.)
127 2484 : out%data = 0
128 59 : out%mask => mask
129 59 : if (present(avg)) out%avg = avg
130 19 : end function newOutputVariable
131 :
132 : !> \brief Update OutputVariable
133 : !> \details Add the array given as actual argument to the derived type's component 'data'
134 : !> \changelog
135 : !! - Robert Schweppe Jun 2018
136 : !! - refactoring and reformatting
137 : !> \return type(OutputVariable)
138 : !> \authors David Schaefer
139 : !> \date June 2015
140 517200 : subroutine updateVariable(self, data)
141 : implicit none
142 :
143 : class(OutputVariable), intent(inout) :: self
144 : real(dp), intent(in), dimension(:) :: data !< data for current time step
145 :
146 23754600 : self%data = self%data + data
147 517200 : self%counter = self%counter + 1
148 :
149 59 : end subroutine updateVariable
150 :
151 : !> \brief Write timestep to file
152 : !> \details Write the content of the derived types's component 'data' to file, average if necessary
153 : !> \changelog
154 : !! - Robert Schweppe Jun 2018
155 : !! - refactoring and reformatting
156 : !> \authors David Schafer
157 : !> \date June 2015
158 1463 : subroutine writeVariableTimestep(self, current_time_step)
159 : implicit none
160 :
161 : class(OutputVariable), intent(inout) :: self
162 :
163 : !> index along the time dimension of the netcdf variable
164 : integer(i4), intent(in) :: current_time_step
165 :
166 1463 : if (self%avg) then
167 45145 : self%data = self%data / real(self%counter, dp)
168 : end if
169 :
170 0 : call self%nc%setData(unpack(self%data, self%mask, nodata_dp), &
171 5852 : (/1, 1, current_time_step/))
172 53028 : self%data = 0
173 1463 : self%counter = 0
174 :
175 517200 : end subroutine writeVariableTimestep
176 :
177 : !> \brief Initialize OutputDataset
178 : !> \details Create and initialize the output file. If new a new output
179 : !! variable needs to be written, this is the first of two
180 : !! procedures to change (second: updateDataset)
181 : !> \changelog
182 : !! - Robert Schweppe Jun 2018
183 : !! - refactoring and reformatting
184 : !! - Sebastian Mueller Jul 2020
185 : !! - added output for river temperature
186 : !> \return type(OutputDataset)
187 : !> \authors Matthias Zink
188 : !> \date Apr 2013
189 19 : function newOutputDataset( iDomain, level, file_name, double_precision, outputs_frequence, time_reference ) result(out)
190 : implicit none
191 :
192 : integer(i4), intent(in) :: iDomain !< domain id
193 : type(Grid), dimension(:), allocatable, target, intent(in) :: level !< level definitions for all domains
194 : character(*), intent(in) :: file_name !< long name of the variable
195 : logical, intent(in) :: double_precision !< mask on desired level
196 : integer(i4), intent(in) :: outputs_frequence !< write out frequence (-3, -2, -1, 0, >0)
197 : integer(i4), intent(in) :: time_reference !< time stamp reference (0: begin, 1: center, 2: end of time interval)
198 :
199 : type(OutputDataset) :: out
200 :
201 19 : out%nc = createOutputFile(iDomain, level, file_name, double_precision, outputs_frequence, time_reference)
202 19 : out%iDomain = iDomain
203 :
204 : ! check if we need minutes instead of hours as time unit and set the time unit factor accordingly
205 : if ( (outputs_frequence > 0) &
206 : .and. (mod(timestep * outputs_frequence, 2) == 1) &
207 19 : .and. (time_reference == 1) &
208 : ) then
209 0 : out%time_unit_factor = 60
210 : else
211 19 : out%time_unit_factor = 1
212 : end if
213 19 : out%outputs_frequence = outputs_frequence
214 19 : out%time_reference = time_reference
215 19 : out%double_precision = double_precision
216 19 : out%previous_time = 0
217 19 : out%counter = 0
218 :
219 1463 : end function newOutputDataset
220 :
221 : !> \brief Write all accumulated data.
222 : !> \details Write all accumulated and potentially averaged data to disk.
223 : !> \changelog
224 : !! - Robert Schweppe Jun 2018
225 : !! - refactoring and reformatting
226 : !> \authors David Schaefer
227 : !> \date June 2015
228 825 : subroutine writeTimestep(self, current_time_step)
229 : implicit none
230 :
231 : class(OutputDataset), intent(inout), target :: self
232 :
233 : !> The model timestep to write
234 : integer(i4), intent(in) :: current_time_step
235 :
236 : integer(i4) :: ii
237 :
238 : type(NcVariable) :: tvar
239 :
240 825 : self%counter = self%counter + 1
241 :
242 : ! add to time variable
243 825 : tvar = self%nc%getVariable("time")
244 1650 : select case( self%time_reference )
245 : case(0)
246 1650 : call tvar%setData(self%previous_time * self%time_unit_factor, (/self%counter/))
247 : case(1)
248 0 : call tvar%setData((self%previous_time + current_time_step) * self%time_unit_factor / 2, (/self%counter/))
249 : case(2)
250 825 : call tvar%setData(current_time_step * self%time_unit_factor, (/self%counter/))
251 : end select
252 : ! add bounds (with current time at end)
253 825 : tvar = self%nc%getVariable("time_bnds")
254 2475 : call tvar%setData(self%previous_time * self%time_unit_factor, (/1, self%counter/))
255 2475 : call tvar%setData(current_time_step * self%time_unit_factor, (/2, self%counter/))
256 825 : self%previous_time = current_time_step
257 :
258 2288 : do ii = 1, size(self%vars)
259 2288 : call self%vars(ii)%writeVariableTimestep(self%counter)
260 : end do
261 :
262 19 : end subroutine writeTimestep
263 :
264 : !> \brief Close the file
265 : !> \details Close the file associated with variable of type(OutputDataset)
266 : !> \changelog
267 : !! - Stephan Thober Oct 2015
268 : !! - adapted to mRM
269 : !! - Robert Schweppe Jun 2018
270 : !! - refactoring and reformatting
271 : !> \authors Rohini Kumar & Stephan Thober
272 : !> \date August 2013
273 19 : subroutine close(self)
274 :
275 825 : use mo_String_utils, only : num2str
276 : use mo_message, only : message
277 :
278 : implicit none
279 :
280 : class(OutputDataset) :: self
281 :
282 19 : call self%nc%close()
283 19 : call message(' OUTPUT: saved netCDF file for domain', trim(num2str(self%iDomain)))
284 19 : call message(' to ', trim(self%nc%fname))
285 :
286 19 : end subroutine close
287 :
288 : !> \brief Create and initialize output file for X & Y coordinate system
289 : !> \details Create output file, write all non-dynamic variables
290 : !! and global attributes for the given domain for X & Y coordinate system
291 : !> \changelog
292 : !! - Stephan Thober Oct 2015
293 : !! - adapted to mRM
294 : !! - Robert Schweppe Jun 2018
295 : !! - refactoring and reformatting
296 : !! - Pallav Shrestha Mar 2020
297 : !! - output file lat and lon are 1d or 2d based on coordinate system
298 : !> \return type(NcDataset)
299 : !> \authors David Schaefer
300 : !> \date June 2015
301 19 : function createOutputFile(iDomain, level, file_name, double_precision, outputs_frequence, time_reference) result(nc)
302 :
303 19 : use mo_common_mhm_mrm_variables, only : evalPer
304 : use mo_grid, only : geoCoordinates, mapCoordinates
305 : use mo_julian, only : dec2date
306 :
307 : implicit none
308 :
309 : ! -> domain id
310 : integer(i4), intent(in) :: iDomain !< selected domain
311 : type(Grid), dimension(:), allocatable, intent(in) :: level !< level definitions for all domains
312 : character(*), intent(in) :: file_name !< long name of the variable
313 : logical, intent(in) :: double_precision !< mask on desired level
314 : integer(i4), intent(in) :: outputs_frequence !< write out frequence (-3, -2, -1, 0, >0)
315 : integer(i4), intent(in) :: time_reference !< time stamp reference (0: begin, 1: center, 2: end of time interval)
316 :
317 : type(NcDataset) :: nc
318 : type(NcDimension), dimension(4) :: dimids1
319 : type(NcVariable) :: var
320 : integer(i4) :: day, month, year
321 : character(1028) :: fname
322 : character(128) :: unit, date, time, datetime
323 19 : real(dp), allocatable, dimension(:) :: easting, northing
324 19 : real(dp), allocatable, dimension(:, :) :: x_bnds, y_bnds
325 19 : real(dp) :: half_step
326 19 : real(dp), allocatable, dimension(:) :: lat1d, lon1d ! 1D lat lon vectors. Used if coordinate system is lat & lon
327 19 : real(dp), allocatable, dimension(:, :) :: lat2d, lon2d ! temporary storage of mHM's 2D latlon array.
328 : ! Used as 2d lat lon arrays if coordinate system is X & Y
329 : character(3) :: dtype
330 :
331 38 : dtype = data_dtype(double_precision)
332 :
333 : ! half cell step to calculate cell bounds from center
334 19 : half_step = level(iDomain)%cellsize / 2.0_dp
335 :
336 19 : fname = trim(dirOut(iDomain)) // trim(file_name)
337 19 : call geoCoordinates(level(iDomain), lat2d, lon2d)
338 :
339 19 : nc = NcDataset(trim(fname), "w")
340 :
341 : ! set the horizonal dimensions
342 19 : if (iFlag_cordinate_sys == 0) then
343 :
344 : ! X & Y coordinate system; 2D lat lon!
345 : !============================================================
346 19 : call mapCoordinates(level(iDomain), northing, easting)
347 57 : allocate(x_bnds(2, size(easting)))
348 57 : allocate(y_bnds(2, size(northing)))
349 143 : x_bnds(1, :) = easting - half_step
350 143 : x_bnds(2, :) = easting + half_step
351 200 : y_bnds(1, :) = northing - half_step
352 200 : y_bnds(2, :) = northing + half_step
353 :
354 : dimids1 = (/&
355 0 : nc%setDimension("easting", size(easting)), &
356 0 : nc%setDimension("northing", size(northing)), &
357 : nc%setDimension("time", 0), &
358 : nc%setDimension("bnds", 2) &
359 114 : /)
360 : ! easting
361 38 : var = nc%setVariable("easting", dtype, (/ dimids1(1) /))
362 : call set_attributes(var, "x-coordinate in the given coordinate system", &
363 19 : unit="m", standard_name="projection_x_coordinate", add_coords=.false., axis="X", bounds="easting_bnds")
364 19 : call var%setData(easting)
365 57 : var = nc%setVariable("easting_bnds", dtype, (/ dimids1(4), dimids1(1) /))
366 19 : call var%setData(x_bnds)
367 : ! northing
368 38 : var = nc%setVariable("northing", dtype, (/ dimids1(2) /))
369 : call set_attributes(var, "y-coordinate in the given coordinate system", &
370 19 : unit="m", standard_name="projection_y_coordinate", add_coords=.false., axis="Y", bounds="northing_bnds")
371 19 : call var%setData(northing)
372 57 : var = nc%setVariable("northing_bnds", dtype, (/ dimids1(4), dimids1(2) /))
373 19 : call var%setData(y_bnds)
374 : ! lon
375 19 : var = nc%setVariable("lon", dtype, dimids1(1 : 2))
376 : call set_attributes(var, "longitude", unit="degrees_east", &
377 19 : double_precision=double_precision, add_coords=.false., standard_name="longitude")
378 19 : call var%setData(lon2d)
379 : ! lat
380 19 : var = nc%setVariable("lat", dtype, dimids1(1 : 2))
381 : call set_attributes(var, "latitude", unit="degrees_north", &
382 19 : double_precision=double_precision, add_coords=.false., standard_name="latitude")
383 38 : call var%setData(lat2d)
384 :
385 : else
386 :
387 : ! lat & lon coordinate system; 1D lat lon!
388 : !============================================================
389 0 : lon1d = lon2d(:, 1) ! first column info is sufficient
390 0 : lat1d = lat2d(1, :) ! first row info is sufficient
391 0 : allocate(x_bnds(2, size(lon1d)))
392 0 : allocate(y_bnds(2, size(lat1d)))
393 : ! cellsize is given in degree in case of lat-lon coordinates
394 0 : x_bnds(1, :) = lon1d - half_step
395 0 : x_bnds(2, :) = lon1d + half_step
396 0 : y_bnds(1, :) = lat1d - half_step
397 0 : y_bnds(2, :) = lat1d + half_step
398 :
399 : dimids1 = (/&
400 0 : nc%setDimension("lon", size(lon1d)), &
401 0 : nc%setDimension("lat", size(lat1d)), &
402 : nc%setDimension("time", 0), &
403 : nc%setDimension("bnds", 2) &
404 0 : /)
405 : ! lon
406 0 : var = nc%setVariable("lon", dtype, (/ dimids1(1) /)) ! sufficient to store lon as vector
407 : call set_attributes(var, "longitude", unit="degrees_east", &
408 0 : add_coords=.false., standard_name="longitude", axis="X", bounds="lon_bnds")
409 0 : call var%setData(lon1d)
410 0 : var = nc%setVariable("lon_bnds", dtype, (/ dimids1(4), dimids1(1) /))
411 0 : call var%setData(x_bnds)
412 : ! lat
413 0 : var = nc%setVariable("lat", dtype, (/ dimids1(2) /)) ! sufficient to store lat as vector
414 : call set_attributes(var, "latitude", unit="degrees_north", &
415 0 : add_coords=.false., standard_name="latitude", axis="Y", bounds="lat_bnds")
416 0 : call var%setData(lat1d)
417 0 : var = nc%setVariable("lat_bnds", dtype, (/ dimids1(4), dimids1(2) /))
418 0 : call var%setData(y_bnds)
419 :
420 : endif
421 :
422 : ! set record dimension
423 : ! time units
424 19 : call dec2date(real(evalPer(iDomain)%julStart, dp), dd = day, mm = month, yy = year)
425 :
426 : ! check if we need minutes instead of hours as time unit
427 : if ( (outputs_frequence > 0) & ! only for output after n steps
428 : .and. (mod(timestep * outputs_frequence, 2) == 1) & ! only for uneven hours
429 19 : .and. (time_reference == 1) & ! only when center of time span is used for time stamp
430 : ) then
431 0 : write(unit, "('minutes since ', i4, '-' ,i2.2, '-', i2.2, 1x, '00:00:00')") year, month, day
432 : else
433 19 : write(unit, "('hours since ', i4, '-' ,i2.2, '-', i2.2, 1x, '00:00:00')") year, month, day
434 : end if
435 :
436 : ! time
437 38 : var = nc%setVariable("time", "i32", (/ dimids1(3) /))
438 19 : call set_attributes(var, "time", unit=unit, add_coords=.false., standard_name="time", axis="T", bounds="time_bnds")
439 57 : var = nc%setVariable("time_bnds", "i32", (/ dimids1(4), dimids1(3) /))
440 :
441 : ! global attributes
442 19 : call date_and_time(date = date, time = time)
443 19 : write(datetime, "(a4,'-',a2,'-',a2,1x,a2,':',a2,':',a2)") date(1 : 4), &
444 38 : date(5 : 6), date(7 : 8), time(1 : 2), time(3 : 4), time(5 : 6)
445 :
446 19 : call nc%setAttribute("project", project_details)
447 19 : call nc%setAttribute("setup_description", setup_description)
448 19 : call nc%setAttribute("simulation_type", simulation_type)
449 19 : call nc%setAttribute("Conventions", Conventions)
450 19 : call nc%setAttribute("contact", contact)
451 19 : call nc%setAttribute("mHM_details", trim(mHM_details) // ", release mHMv" // trim(version))
452 19 : call nc%setAttribute("history", trim(datetime) // ", " // history)
453 19 : call nc%setAttribute("title", "mHMv"//trim(version)//" "//trim(simulation_type)//" outputs")
454 19 : call nc%setAttribute("creation_date", datetime)
455 :
456 19 : end function createOutputFile
457 :
458 : !> \brief Write output variable attributes
459 : !> \details Modifications:
460 : !! - Robert Schweppe Jun 2018 - refactoring and reformatting
461 : !> \authors David Schaefer
462 : !> \date June 2015
463 154 : subroutine set_attributes(var, long_name, unit, double_precision, add_coords, standard_name, axis, bounds)
464 : implicit none
465 :
466 : type(NcVariable), intent(inout) :: var !< NetCDF variable
467 : character(*), intent(in), optional :: long_name !< long name of the variable
468 : character(*), intent(in), optional :: unit !< unit of the variable
469 : logical, intent(in), optional :: double_precision !< precision flag, if missing, no fill-value and missing-value is written
470 : logical, intent(in), optional :: add_coords !< whether to add the "coordiantes" attribute "lat lon", .true. by defult
471 : character(*), intent(in), optional :: standard_name !< standard name of the variable
472 : character(*), intent(in), optional :: axis !< axis attribute
473 : character(*), intent(in), optional :: bounds !< bounds attribute
474 :
475 : logical :: add_coords_
476 :
477 154 : add_coords_ = .true.
478 95 : if (present(add_coords)) add_coords_ = add_coords
479 :
480 154 : if (present(double_precision)) then
481 97 : if (double_precision) then
482 97 : call var%setFillValue(nodata_dp)
483 97 : call var%setAttribute("missing_value", nodata_dp)
484 : ! call var%setAttribute("scale_factor", 1.0_dp) ! not needed if just 1
485 : else
486 0 : call var%setFillValue(real(nodata_dp, kind=sp))
487 0 : call var%setAttribute("missing_value", real(nodata_dp, kind=sp))
488 : ! call var%setAttribute("scale_factor", 1.0_sp) ! not needed if just 1
489 : end if
490 : end if
491 154 : if (present(long_name)) call var%setAttribute("long_name", long_name)
492 154 : if (present(unit)) call var%setAttribute("units", unit)
493 154 : if (add_coords_) call var%setAttribute("coordinates", "lat lon")
494 154 : if (present(standard_name)) call var%setAttribute("standard_name", standard_name)
495 154 : if (present(axis)) call var%setAttribute("axis", axis)
496 154 : if (present(bounds)) call var%setAttribute("bounds", bounds)
497 :
498 19 : end subroutine set_attributes
499 :
500 154 : end module mo_nc_output
|
