libsim  Versione6.3.0
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 
19 module volgrid6d_vapor_class
20 
21 use log4fortran
24 use vol7d_class
28 use vdf4f
29 use err_handling
31 use grid_class
32 
33 implicit none
34 
35 
37 INTERFACE export
38  MODULE PROCEDURE volgrid6d_export_to_vapor
39 END INTERFACE
40 
41 private
42 
43 public export
44 
45 contains
46 
49 subroutine volgrid6d_export_to_vapor (this,normalize,rzscan,filename,filename_auto,reusevdf)
50 
51 TYPE(volgrid6d),INTENT(INOUT) :: this
52 logical,intent(in) :: normalize
53 logical,intent(in),optional :: rzscan
54 character(len=*),intent(in),optional :: filename
55 character(len=*),intent(out),optional :: filename_auto
56 logical,intent(in),optional :: reusevdf
57 
58 character(len=254) :: lfilename
59 integer :: ntime, ntimerange, ntimera, nlevel, nvar
60 logical :: exist
61 integer :: ier, xyzdim(3),ivar,i,j
62 character(len=255),allocatable :: varnames(:),vardescriptions(:),tsdescriptions(:)
63 doubleprecision :: extents(6),vdfmiss=rmiss
64 
65 TYPE(conv_func), pointer :: c_func(:)
66 TYPE(vol7d_var),allocatable :: varbufr(:)
67 type(vol7d_var),pointer :: dballevar(:)
68 CHARACTER(len=255) :: proj_type,mapprojection
69 
70 integer :: zone, irzscan, indele
71 DOUBLE PRECISION :: xoff, yoff, ellips_smaj_axis, ellips_flatt
72 DOUBLE PRECISION :: longitude_south_pole, latitude_south_pole, angle_rotation
73 
74 call l4f_category_log(this%category,l4f_debug,"export to vapor")
75 
76 call vol7d_dballe_import_dballevar(dballevar)
77 
78 if (optio_log(rzscan)) then
79  irzscan=1
80 else
81  irzscan=0
82 end if
83 
84 ntime=imiss
85 ntimerange=imiss
86 nlevel=imiss
87 nvar=imiss
88 ntimera=imiss
89 
90 lfilename="vg6d.vdf"
91 if (present(filename))then
92  if (filename /= "")then
93  lfilename=filename
94  end if
95 end if
96 
97 if (present(filename_auto))filename_auto=lfilename
98 
99 
100 inquire(file=lfilename,exist=exist)
101 
102 if (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
107 else if (exist) then
108  call l4f_category_log(this%category,l4f_error,"file exist; cannot open new file: "//trim(lfilename))
109  CALL raise_error()
110 end if
111 
112 call l4f_category_log(this%category,l4f_info,"writing on file: "//trim(lfilename))
113 
114 if (associated(this%time)) ntime=size(this%time)
115 if (associated(this%timerange)) ntimerange=size(this%timerange)
116 if (associated(this%level)) nlevel=size(this%level)
117 if (associated(this%var)) nvar=size(this%var)
118 
119 if (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 
410 else
411 
412  call l4f_category_log(this%category,l4f_warn,"volume with some dimensions to 0: not exported to vdf")
413 
414 end if
415 
416 contains
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 
429 subroutine fill_underground_missing_values(voldati)
430 real,intent(inout) :: voldati(:,:,:,:,:)
431 
432 integer :: x,y,z,tim,var,zz
433 
434 do 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
456 end do
457 
458 
459 end subroutine fill_underground_missing_values
460 
461 end subroutine volgrid6d_export_to_vapor
462 
463 end module volgrid6d_vapor_class
Functions that return a trimmed CHARACTER representation of the input variable.
Represent level object in a pretty string.
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
Apply the conversion function this to values.
Module for describing geographically referenced regular grids.
Definition: grid_class.F90:262
This module defines objects and methods for managing data volumes on rectangular georeferenced grids...
classe per import ed export di volumi da e in DB-All.e
Method for returning the contents of the object.
Gestione degli errori.
Classe per la gestione di un volume completo di dati osservati.
Definitions of constants and functions for working with missing values.
interface to different architectures (cast some type)
Definition: vdf4f.F03:383
classe per la gestione del logging
Class for managing physical variables in a grib 1/2 fashion.
Utilities for CHARACTER variables.
Emit log message for a category with specific priority.

Generated with Doxygen.