libsim Versione 7.2.1
volgrid6d_vapor_class.F90
1! Copyright (C) 2011 ARPA-SIM <urpsim@smr.arpa.emr.it>
2! authors:
3! Davide Cesari <dcesari@arpa.emr.it>
4! Paolo Patruno <ppatruno@arpa.emr.it>
5
6! This program is free software; you can redistribute it and/or
7! modify it under the terms of the GNU General Public License as
8! published by the Free Software Foundation; either version 2 of
9! the License, or (at your option) any later version.
10
11! This program is distributed in the hope that it will be useful,
12! but WITHOUT ANY WARRANTY; without even the implied warranty of
13! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14! GNU General Public License for more details.
15
16! You should have received a copy of the GNU General Public License
17! along with this program. If not, see <http://www.gnu.org/licenses/>.
18
19module volgrid6d_vapor_class
20
28use vdf4f
31use grid_class
32
33implicit none
34
35
37INTERFACE export
38 MODULE PROCEDURE volgrid6d_export_to_vapor
39END INTERFACE
40
41private
42
43public export
44
45contains
46
49subroutine volgrid6d_export_to_vapor (this,normalize,rzscan,filename,filename_auto,reusevdf)
50
51TYPE(volgrid6d),INTENT(INOUT) :: this
52logical,intent(in) :: normalize
53logical,intent(in),optional :: rzscan
54character(len=*),intent(in),optional :: filename
55character(len=*),intent(out),optional :: filename_auto
56logical,intent(in),optional :: reusevdf
57
58character(len=254) :: lfilename
59integer :: ntime, ntimerange, ntimera, nlevel, nvar
60logical :: exist
61integer :: ier, xyzdim(3),ivar,i,j
62character(len=255),allocatable :: varnames(:),vardescriptions(:),tsdescriptions(:)
63doubleprecision :: extents(6),vdfmiss=rmiss
64
65TYPE(conv_func), pointer :: c_func(:)
66TYPE(vol7d_var),allocatable :: varbufr(:)
67type(vol7d_var),pointer :: dballevar(:)
68CHARACTER(len=255) :: proj_type,mapprojection
69
70integer :: zone, irzscan, indele
71DOUBLE PRECISION :: xoff, yoff, ellips_smaj_axis, ellips_flatt
72DOUBLE PRECISION :: longitude_south_pole, latitude_south_pole, angle_rotation
73
74call l4f_category_log(this%category,l4f_debug,"export to vapor")
75
76call vol7d_dballe_import_dballevar(dballevar)
77
78if (optio_log(rzscan)) then
79 irzscan=1
80else
81 irzscan=0
82end if
83
84ntime=imiss
85ntimerange=imiss
86nlevel=imiss
87nvar=imiss
88ntimera=imiss
89
90lfilename="vg6d.vdf"
91if (present(filename))then
92 if (filename /= "")then
93 lfilename=filename
94 end if
95end if
96
97if (present(filename_auto))filename_auto=lfilename
98
99
100inquire(file=lfilename,exist=exist)
101
102if (optio_log(reusevdf)) then
103 if (.not. exist) then
104 call l4f_category_log(this%category,l4f_error,"file do not exist; cannot reuse file: "//trim(lfilename))
105 CALL raise_error()
106 end if
107else if (exist) then
108 call l4f_category_log(this%category,l4f_error,"file exist; cannot open new file: "//trim(lfilename))
109 CALL raise_error()
110end if
111
112call l4f_category_log(this%category,l4f_info,"writing on file: "//trim(lfilename))
113
114if (associated(this%time)) ntime=size(this%time)
115if (associated(this%timerange)) ntimerange=size(this%timerange)
116if (associated(this%level)) nlevel=size(this%level)
117if (associated(this%var)) nvar=size(this%var)
118
119if (c_e(ntime) .and. c_e(ntimerange) .and. c_e(nlevel) .and. c_e(nvar)) then
120
121 allocate(varnames(nvar),vardescriptions(nvar),varbufr(nvar))
122
123 call get_val (this%griddim, nx=xyzdim(1) , ny=xyzdim(2) )
124
125 !xyzdim(1)=this%griddim%dim%nx
126 !xyzdim(2)=this%griddim%dim%ny
127 xyzdim(3)=nlevel
128
129 if (associated(this%voldati)) then
130
131 if (optio_log(normalize)) then
132 CALL vargrib2varbufr(this%var, varbufr, c_func)
133
134 ! Rescale valid data according to variable conversion table
135 IF (ASSOCIATED(c_func)) THEN
136
137 call l4f_category_log(this%category,l4f_info,"normalize is activated, so the volume data are changed in output")
138 DO ivar = 1, nvar
139 this%voldati(:,:,:,:,:,ivar) = convert(c_func(ivar),this%voldati(:,:,:,:,:,ivar))
140 ENDDO
141 DEALLOCATE(c_func)
142
143 ENDIF
144
145 indele=imiss
146 do ivar=1,nvar
147
148 j=firsttrue(varbufr(ivar)%btable == dballevar(:)%btable)
149
150 if ( j > 0 )then
151 varbufr(ivar)%description = dballevar(j)%description
152 varbufr(ivar)%unit = dballevar(j)%unit
153 varbufr(ivar)%scalefactor = dballevar(j)%scalefactor
154
155 varnames(ivar) = varbufr(ivar)%btable
156 vardescriptions(ivar) = trim(varbufr(ivar)%description)//"_"//trim(varbufr(ivar)%unit)
157 !vardescriptions(ivar) = "V"//wash_char(trim(varbufr(ivar)%description)//"_"//trim(varbufr(ivar)%unit))
158
159 if (varnames(ivar) == "B10007") then
160 indele=ivar
161 end if
162
163 else
164
165 varnames(ivar) = "Vnotnormalized_"//t2c(ivar)
166 vardescriptions(ivar) = "None"
167
168 end if
169
170 end do
171 else
172
173 do ivar=1, nvar
174 varnames(ivar) = "V"//trim(to_char(this%var(ivar)%number))
175 !varnames(ivar) = wash_char("VAR_"//trim(to_char(this%var(ivar)%number)))
176 end do
177 end if
178
179 if (this%time_definition == 1) then
180 ntimera=ntime
181 allocate(tsdescriptions(ntimera))
182
183 do i=1,ntimera
184 tsdescriptions(i)=to_char(this%time(i))
185 end do
186 else
187 ntimera=ntimerange
188 allocate(tsdescriptions(ntimera))
189
190 do i=1,ntimera
191 tsdescriptions(i)=to_char(this%timerange(i))
192 end do
193 end if
194
195 !extents=(/-0.5d0, -14.d0 , 0.d0, 2.563d0, -11.25d0, 1.d0/)
196 call get_val (this%griddim, xmin=extents(1),ymin=extents(2), xmax=extents(4) , ymax=extents(5))
197
198 extents(3)=0.d0
199 extents(6)=10.d0
200
201 call get_val (this%griddim, proj_type=proj_type)
202
203! Alan Norton Wed, 11 May 2011 11:04:36 -0600
204! We do not currently support lat-lon or rotated lat-lon; the only supported
205! projections are lambert, polar stereographic, and mercator.
206! I do expect that we will support lat-lon in our next release.
207
208! VAPOR Release Notes
209! Version 2.1.0
210! New Features
211! Support is provided for rotated lat-lon (Cassini) transform, as used by WRF-ARW version 3.
212
213 select case (proj_type)
214 case ("regular_ll")
215
216 call l4f_category_log(this%category,l4f_info,"VDF: probably vapor do not support this projection ?: "//trim(proj_type))
217 mapprojection = "+proj=latlon +ellps=sphere"
218
219 extents(1)=extents(1)*111177.d0
220 extents(2)=extents(2)*111177.d0
221 extents(3)=extents(3)*100000.d0
222 extents(4)=extents(4)*111177.d0
223 extents(5)=extents(5)*111177.d0
224 extents(6)=extents(6)*100000.d0
225
226 case ("rotated_ll")
227
228 !call l4f_category_log(this%category,L4F_WARN,"VDF: vapor and proj do not support this projection: "//trim(proj_type))
229 !call l4f_category_log(this%category,L4F_WARN,"VDF: you have to considerate it as rotated coordinates")
230
231 call l4f_category_log(this%category,l4f_info,"VDF: vapor probably support this projection: "//trim(proj_type))
232
233 call get_val (this%griddim, longitude_south_pole=longitude_south_pole,&
234 latitude_south_pole=latitude_south_pole, angle_rotation=angle_rotation)
235 !write (mapprojection,*)"-m 57.295779506 +proj=ob_tran +o_proj=latlon +o_lat_p=",32.5"," +o_lon_p=",0"," +lon_0=",10
236 !mapprojection = "+proj=latlon"
237
238 !cassini or rotated lat/lon
239
240 if (angle_rotation /= 0. ) then
241 call l4f_category_log(this%category,l4f_error,"angle_rotation /= 0 not supported in vapor (proj)")
242 call raise_error()
243 end if
244
245 mapprojection = "+proj=ob_tran +o_proj=latlong +o_lat_p="//t2c(-latitude_south_pole)//&
246 "d +o_lon_p=0d +lon_0="//t2c(longitude_south_pole)//"d +ellps=sphere"
247
248 extents(1)=extents(1)*111177.d0
249 extents(2)=extents(2)*111177.d0
250 extents(3)=extents(3)*100000.d0
251 extents(4)=extents(4)*111177.d0
252 extents(5)=extents(5)*111177.d0
253 extents(6)=extents(6)*100000.d0
254
255 case ("UTM")
256
257 call l4f_category_log(this%category,l4f_warn,"VDF: vapor do not support this projection: "//trim(proj_type))
258 call get_val (this%griddim, xmin=extents(1),ymin=extents(2), xmax=extents(4) , ymax=extents(5),&
259 zone=zone, xoff=xoff, yoff=yoff, ellips_smaj_axis=ellips_smaj_axis, ellips_flatt=ellips_flatt)
260
261 extents(3)=0.d0
262 extents(6)=300000.d0
263
264 mapprojection ="+proj=utm +zone="//t2c(zone)
265
266 !mapprojection ="+proj=utm +zone="//t2c(zone)//" +a="//t2c(ellips_smaj_axis)//&
267 ! " +f="//t2c(ellips_flatt)//" +x_0="//t2c(xoff)//" +y_0="//t2c(yoff)
268
269 case default
270
271 call l4f_category_log(this%category,l4f_warn,&
272 "VDF: proj or vdf (vapor) export do not support this projection: "//trim(proj_type))
273 mapprojection = cmiss
274
275 end select
276
277 ier=0
278 if(ier==0) call l4f_category_log(this%category,l4f_info,"VDF: projection parameter "//mapprojection)
279
280 if (reusevdf) then
281 if(ier==0) call l4f_category_log(this%category,l4f_debug,"VDF: call create_metadata_from_file")
282 ier = vdf4f_create_metadata_from_file(lfilename)
283 else
284 if(ier==0) call l4f_category_log(this%category,l4f_debug,"VDF: call create_metadata")
285 ier = vdf4f_create_metadata(xyzdim,vdctype=2)
286
287 if(ier==0) call l4f_category_log(this%category,l4f_debug,"VDF: call set_missing_value")
288 if(ier==0) ier = vdf4f_set_missing_value(vdfmiss)
289
290 if(ier==0) call l4f_category_log(this%category,l4f_debug,"VDF: call set_num_timesteps")
291 if(ier==0) ier = vdf4f_set_num_timesteps(ntimera)
292
293 if(ier==0) call l4f_category_log(this%category,l4f_debug,"VDF: call vdf4f_set_comment")
294 if(ier==0) ier = vdf4f_set_comment("vogrid6d exported")
295 !if(ier==0) ier = vdf4f_set_coord_system_type(coordsystemtype="spherical")
296 if(ier==0) ier = vdf4f_set_coord_system_type(coordsystemtype="cartesian")
297
298 if (c_e(indele)) then
299 if(ier==0) call l4f_category_log(this%category,l4f_info,"VDF: ELEVATION (B10007) found: setting gridtype to layered")
300 extents(6) = maxval(this%voldati(:,:,:,:,:,indele),c_e(this%voldati(:,:,:,:,:,indele)))
301 if(ier==0) ier = vdf4f_set_grid_type(gridtype="layered")
302 else
303 if(ier==0) ier = vdf4f_set_grid_type(gridtype="regular")
304 end if
305
306 if(ier==0) ier = vdf4f_set_grid_extents(extents=extents)
307 if(ier==0) ier = vdf4f_set_map_projection(mapprojection=mapprojection)
308
309 !!!! if(ier==0) int vdf4f_set_grid_permutation(long permutation[3]);
310
311 end if
312
313 if (nlevel > 1) then
314
315 if (c_e(indele)) varnames(indele) = "ELEVATION"
316
317 if(ier==0) call l4f_category_log(this%category,l4f_debug,"VDF: call set_variables_names")
318 if(ier==0) ier = vdf4f_set_variables_names(nvar, varnames)
319
320
321 if(ier==0) call l4f_category_log(this%category,l4f_debug,"VDF: call set_v_comment")
322 if(ier==0) call l4f_category_log(this%category,l4f_debug,"VDF: call set_ts_comment")
323
324 do i=1,ntimera
325 if(ier==0) ier = vdf4f_set_ts_comment(i-1,tsdescriptions(i))
326 do j=1,nvar
327 if(ier==0) ier = vdf4f_set_v_comment(i-1,varnames(j),vardescriptions(j))
328 end do
329 end do
330
331 else
332
333 do ivar=1,nvar
334 varnames(ivar)="XY_"//t2c(varnames(ivar))
335 end do
336
337 if (c_e(indele)) varnames(indele) = "HGT"
338
339 if(ier==0) call l4f_category_log(this%category,l4f_debug,"VDF: call vdf4f_set_variables_2d_xy")
340 if(ier==0) ier = vdf4f_set_variables_2d_xy(nvar, varnames)
341
342 end if
343
344 if(ier==0) call l4f_category_log(this%category,l4f_debug,"VDF: call write_metadata")
345 if(ier==0) ier = vdf4f_write_metadata(lfilename)
346
347 if(ier==0) ier = destroy_metadata_c()
348 ! end metadata section
349
350 ! start data section
351
352 if(ier==0) call l4f_category_log(this%category,l4f_debug,"VDF: call vdf4f_create_writer")
353 if(ier==0) ier = vdf4f_create_writer(lfilename)
354
355 if (this%time_definition == 1) then
356
357 if (ntimerange /= 1) then
358 if(ier==0) call l4f_category_log(this%category,l4f_warn,"VDF: writing only first timerange, there are:"//t2c(ntimerange))
359 end if
360
361 if (.not. c_e(indele)) call fill_underground_missing_values(this%voldati(:,:,:,:,1,:))
362 if(ier==0) call l4f_category_log(this%category,l4f_info,"scan VDF (vapor file) for times")
363
364 if (nlevel > 1) then
365 if(ier==0) call l4f_category_log(this%category,l4f_debug,"VDF: call vdf4f_write")
366 if(ier==0) ier = vdf4f_write(this%voldati(:,:,:,:,1,:), xyzdim, ntime, nvar, varnames, irzscan)
367 else
368 if(ier==0) call l4f_category_log(this%category,l4f_debug,"VDF: call vdf4f_write_2d_xy")
369 if(ier==0) ier = vdf4f_write_2d_xy(this%voldati(:,:,1,:,1,:), xyzdim(:2), ntime, nvar ,varnames)
370
371 end if
372
373 else
374
375 if (ntime /= 1) then
376 if(ier==0) call l4f_category_log(this%category,l4f_warn,"VDF: writing only fisth time, there are:"//t2c(ntime))
377 end if
378
379 if (.not. c_e(indele)) call fill_underground_missing_values(this%voldati(:,:,:,1,:,:))
380 if(ier==0) call l4f_category_log(this%category,l4f_info,"scan VDF (vapor file) for timeranges")
381
382 if (nlevel > 1) then
383 if(ier==0) call l4f_category_log(this%category,l4f_debug,"VDF: call vdf4f_write")
384 if(ier==0) ier = vdf4f_write(this%voldati(:,:,:,1,:,:), xyzdim, ntimerange, nvar, varnames, irzscan)
385 else
386 if(ier==0) call l4f_category_log(this%category,l4f_debug,"VDF: call vdf4f_write_2d_xy")
387 if(ier==0) ier = vdf4f_write_2d_xy(this%voldati(:,:,1,1,:,:), xyzdim(:2), ntimerange, nvar, varnames)
388
389 end if
390
391 end if
392
393 if (ier /= 0) call l4f_category_log(this%category,l4f_error,"export to vdf: "//vdf4f_get_err_msg())
394
395 !todo: check if all are allocated
396 deallocate(varnames,vardescriptions,tsdescriptions,varbufr)
397 if (ier==0) ier = destroy_writer_c()
398
399 if (ier /= 0) then
400 call l4f_category_log(this%category,l4f_error,"exporting to vdf")
401 CALL raise_fatal_error("exporting to vdf")
402 end if
403
404 else
405
406 call l4f_category_log(this%category,l4f_warn,"volume with voldati not associated: not exported to vdf")
407
408 end if
409
410else
411
412 call l4f_category_log(this%category,l4f_warn,"volume with some dimensions to 0: not exported to vdf")
413
414end if
415
416contains
417
418!Vapor does not (yet) have support for missing values. I recommend that
419!you choose a particular floating point value that is different from all
420!normal values, and write that value wherever you find the missing
421!value. You can then make that value map to full transparency when you
422!encounter it in volume rendering.
423!Huge values are OK, but be sure that the values you use are valid
424!floating point numbers. Note also that when these values are displayed
425!at lowered refinement levels, they will be averaged with nearby values.
426!-Alan Norton
427
429subroutine fill_underground_missing_values(voldati)
430real,intent(inout) :: voldati(:,:,:,:,:)
431
432integer :: x,y,z,tim,var,zz
433
434do x=1,size(voldati,1)
435 do y=1,size(voldati,2)
436 do tim=1,size(voldati,4)
437 do var=1,size(voldati,5)
438 do z=1,size(voldati,3)
439
440 if (.not. c_e(voldati(x,y,z,tim,var))) then
441 zz=firsttrue(c_e(voldati(x,y,:,tim,var)))
442 if (zz > 0 ) then
443 voldati(x,y,z,tim,var)=voldati(x,y,firsttrue(c_e(voldati(x,y,:,tim,var))),tim,var)
444 else
445 call l4f_log(l4f_warn,"fill_underground_missing_values: there are only missing values in the full coloumn")
446 exit
447 end if
448 else
449 exit
450 end if
451
452 end do
453 end do
454 end do
455 end do
456end do
457
458
459end subroutine fill_underground_missing_values
460
461end subroutine volgrid6d_export_to_vapor
462
463end module volgrid6d_vapor_class
Set of functions that return a trimmed CHARACTER representation of the input variable.
Set of functions that return a CHARACTER representation of the input variable.
Method for returning the contents of the object.
Definition: grid_class.F90:319
Emit log message for a category with specific priority.
Function to check whether a value is missing or not.
interface to different architectures (cast some type)
Definition: vdf4f.F03:358
Apply the conversion function this to values.
Utilities for CHARACTER variables.
Gestione degli errori.
Module for describing geographically referenced regular grids.
Definition: grid_class.F90:237
classe per la gestione del logging
Definitions of constants and functions for working with missing values.
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
Classe per la gestione di un volume completo di dati osservati.
classe per import ed export di volumi da e in DB-All.e
This module defines objects and methods for managing data volumes on rectangular georeferenced grids.
Class for managing physical variables in a grib 1/2 fashion.
Object describing a rectangular, homogeneous gridded dataset.
Class defining a real conversion function between units.

Generated with Doxygen.