libsim  Versione 7.1.6
v7d_qcspa.F90

Sample program for module qcspa

1 ! Copyright (C) 2010 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 ! Program to quality control data with climatological values
19 
20 #include "config.h"
21 
22 program v7d_qcspa
23 
24 use log4fortran
26 USE simple_stat
29 use dballe_class
30 use modqc
31 use modqcspa
32 !use vol7d_dballeold_class
34 USE vol7d_class
37 
38 #ifdef HAVE_LIBNCARG
39 USE ncar_plot_class
40 #endif
41 
42 implicit none
43 
44 integer :: category,io,ier,i,iun,n,ind
45 character(len=512):: a_name,output_format, output_template
46 
47  !tipi derivati.
48 TYPE(optionparser) :: opt
49 TYPE(geo_coord) :: coordmin, coordmax
50 TYPE(datetime) :: time,ti, tf, timei, timef, timeiqc, timefqc
51 type(qcspatype) :: v7dqcspa
52 type(vol7d_dballe) :: v7ddballe
53 #ifdef HAVE_LIBNCARG
54 type(ncar_plot) :: plot
55 #endif
56 
57 integer, parameter :: maxvar=10
58 character(len=6) :: var(maxvar)=cmiss ! variables to elaborate
59 #ifdef HAVE_DBALLE
60 character(len=512) :: dsn='test1',user='test',password=''
61 character(len=512) :: dsne='test',usere='test',passworde=''
62 character(len=512) :: dsnspa='test',userspa='test',passwordspa=''
63 #endif
64 integer :: years=imiss,months=imiss,days=imiss,hours=imiss,yeare=imiss,monthe=imiss,daye=imiss,houre=imiss,nvar=0
65 doubleprecision :: lons=dmiss,lats=dmiss,lone=dmiss,late=dmiss
66 integer :: year, month, day, hour
67 logical :: height2level=.false.,doplot=.false.,version
68 CHARACTER(len=512) :: input_file, output_file
69 INTEGER :: optind, optstatus, ninput
70 CHARACTER(len=20) :: operation
71 TYPE(ARRAYOF_REAL):: grad
72 real :: val
73 integer :: iostat
74 REAL, DIMENSION(11) :: perc_vals=(/(10.*i,i=0,10)/) !(/0.,10.,20.,30.,40.,50.,60.,70.,80.,90.,100./) !<the percentiles values to be computed, between 0. and 100.
75 REAL, DIMENSION(size(perc_vals)-1) :: ndi
76 REAL, DIMENSION(size(perc_vals)) :: nlimbins
77 double precision, dimension(2) :: percentile
78 TYPE(cyclicdatetime) :: cyclicdt
79 type(vol7d_level) :: level,levelo
80 type(vol7d_timerange) :: timerange,timerangeo
81 type(vol7d_var) :: varia,variao
82 CHARACTER(len=vol7d_ana_lenident) :: ident
83 integer :: indcana,indctime,indclevel,indctimerange,indcdativarr,indcnetwork,indcattr
84 INTEGER,POINTER :: w_s(:), w_e(:)
85 TYPE(vol7d_dballe) :: v7d_dba_out
86 logical :: file
87 integer :: status,grunit
88 
89 #ifdef HAVE_DBALLE
90 namelist /odbc/ dsn,user,password,dsne,usere,passworde
91 namelist /odbcspa/ dsnspa,userspa,passwordspa ! namelist to define DSN
92 #endif
93 namelist /switchspa/ height2level,doplot
94 namelist /minmax/ years,months,days,hours,lons,lats,yeare,monthe,daye,houre,lone,late
95 namelist /varlist/ var
96 
97 !init log4fortran
98 ier=l4f_init()
99 
100 ! unique name from launcher
101 call l4f_launcher(a_name,a_name_force="v7d_qcspa")
102 
103 ! set a_name
104 category=l4f_category_get(a_name//".main")
105 
106 
107 ! define the option parser
108 opt = optionparser_new(description_msg= &
109  'Spatial quality control: compute gradient; compute NDI from gradient; apply quality control', &
110  usage_msg='Usage: v7d_qcspa [options] [inputfile1] [inputfile2...] [outputfile] \n&
111  &If input-format is of file type, inputfile ''-'' indicates stdin,'&
112 #ifdef HAVE_DBALLE
113  //'if output-format is of database type, outputfile specifies &
114  &database access info in URI form,' &
115 #endif
116  //'if empty or ''-'', a suitable default is used. &
117 &')
118 
119 ! options for defining input
120 CALL optionparser_add(opt, ' ', 'operation', operation, "run", help= &
121  'operation to execute: ''gradient'' compute gradient and write on files; '&
122  //'''ndi'' compute NDI from gradient;' &
123  //'''run'' apply quality control ')
124 
125 
126 ! options for defining output
127 output_template = ''
128 CALL optionparser_add(opt, ' ', 'output-format', output_format, 'native', help= &
129  'format of output file for "ndi" operation only, in the form ''name[:template]''; ''native'' for vol7d &
130  &native binary format (no template to be specified)'&
131 #ifdef HAVE_DBALLE
132  //'; ''BUFR'' and ''CREX'' for corresponding formats, with template as an alias like ''synop'', ''metar'', &
133  &''temp'', ''generic'', empty for ''generic'''&
134 #endif
135 )
136 
137 ! help options
138 CALL optionparser_add_help(opt, 'h', 'help', help='show an help message and exit')
139 CALL optionparser_add(opt, ' ', 'version', version, help='show version and exit')
140 
141 ! parse options and check for errors
142 CALL optionparser_parse(opt, optind, optstatus)
143 
144 IF (optstatus == optionparser_help) THEN
145  CALL exit(0) ! generate a clean manpage
146 ELSE IF (optstatus == optionparser_err) THEN
147  CALL l4f_category_log(category,l4f_error,'in command-line parameters')
148  CALL raise_fatal_error()
149 ENDIF
150 IF (version) THEN
151  WRITE(*,'(A,1X,A)')'v7d_qcspa',version
152  CALL exit(0)
153 ENDIF
154 
155 
156 if (operation /= "gradient" .and. operation /= "ndi" .and. operation /= "run") then
157  CALL optionparser_printhelp(opt)
158  CALL l4f_category_log(category, l4f_error, &
159  'argument to --operation is wrong')
160  CALL raise_fatal_error()
161 end if
162 
163 
164 ! check output format/template
165 n = word_split(output_format, w_s, w_e, ':')
166 IF (n >= 2) THEN ! set output template if present
167  output_template = output_format(w_s(2):w_e(2))
168  output_format(w_e(1)+1:) = ' '
169 ENDIF
170 DEALLOCATE(w_s, w_e)
171 
172 if (operation == "ndi") then
173 
174  ! check input/output files
175  i = iargc() - optind
176  IF (i < 0) THEN
177  CALL optionparser_printhelp(opt)
178  CALL l4f_category_log(category,l4f_error,'input file missing')
179  CALL raise_fatal_error()
180  ELSE IF (i < 1) THEN
181  CALL optionparser_printhelp(opt)
182  CALL l4f_category_log(category,l4f_error,'output file missing')
183  CALL raise_fatal_error()
184  ENDIF
185  CALL getarg(iargc(), output_file)
186 
187 
188  !you can define level, timerange, var or get it from file
189  call init(levelo)
190  call init(timerangeo)
191  call init (variao)
192 
193  do ninput = optind, iargc()-1
194  call getarg(ninput, input_file)
195 
196  CALL l4f_category_log(category,l4f_info,"open file: "//t2c(input_file))
197  grunit=getunit()
198  open (grunit,file=input_file,status="old")
199  read (grunit,*,iostat=status) level,timerange,varia
200  if (status /= 0) then
201  CALL l4f_category_log(category,l4f_warn,"error reading: "//t2c(input_file))
202  close(grunit)
203  cycle
204  endif
205 
206  if (c_e(levelo)) then
207 
208  !if (t2c(input_file) /= to_char(levelo)//"_"//to_char(timerangeo)//"_"//variao%btable//".grad")
209  if ( level /= levelo .or. timerange /= timerangeo .or. varia /= variao ) then
210  call l4f_category_log(category,l4f_error,"Error reading grad files: file are incoerent")
211  call raise_error("")
212  end if
213  else
214  levelo = level
215  timerangeo = timerange
216  variao = varia
217  end if
218 
219  do while (.true.)
220  read (grunit,*,iostat=iostat) val
221  if (iostat /= 0) exit
222  if (val /= 0.) call insert(grad,val)
223  end do
224 
225  close(grunit)
226  end do
227 
228  CALL l4f_category_log(category,l4f_info,"compute percentile to remove the tails")
229  percentile = stat_percentile(grad%array(:grad%arraysize),(/10.,90./))
230  !print *,percentile
231 
232  CALL l4f_category_log(category,l4f_info,"compute NDI")
233  call normalizeddensityindex(pack(grad%array(:grad%arraysize),&
234  mask=(grad%array(:grad%arraysize) < percentile(2) )), perc_vals, ndi, nlimbins)
235 
236  call delete(grad)
237  call delete(v7dqcspa%clima)
238  CALL init(v7dqcspa%clima, time_definition=0)
239  call vol7d_alloc(v7dqcspa%clima,nana=size(perc_vals)-1, &
240  nlevel=1, ntimerange=1, &
241  ndativarr=1, nnetwork=1,ntime=1,ndativarattrr=1,ndatiattrr=1)
242 
243  call vol7d_alloc_vol(v7dqcspa%clima,inivol=.true.)
244 
245  call init(v7dqcspa%clima%network(1),name="qcspa-ndi")
246  v7dqcspa%clima%level=level
247  v7dqcspa%clima%timerange=timerange
248  v7dqcspa%clima%dativar%r(1)=varia
249  v7dqcspa%clima%dativarattr%r(1)=varia
250  call init(v7dqcspa%clima%datiattr%r(1), btable="*B33209") ! NDI order number
251  cyclicdt = cyclicdatetime_new(chardate="/////////") !TMMGGhhmm
252  v7dqcspa%clima%time(1)=cyclicdatetime_to_conventional(cyclicdt)
253 
254  indctime=1
255  indclevel=1
256  indcnetwork=1
257  indcattr=1
258  indcdativarr=1
259  indctimerange=1
260 
261  do indcana=1,size(perc_vals)-1
262  write(ident,'("#",i2.2,2i3.3)')0,0,nint(perc_vals(indcana))
263  call init(v7dqcspa%clima%ana(indcana),ident=ident,lat=0d0,lon=0d0)
264  if (c_e(nlimbins(indcana)).and.c_e(ndi(indcana)))then
265  ind=index_c(spa_btable,varia%btable)
266  v7dqcspa%clima%voldatir(indcana, indctime, indclevel, indctimerange, indcdativarr, indcnetwork)=&
267  nlimbins(indcana)*spa_a(ind) + spa_b(ind)
268  v7dqcspa%clima%voldatiattrr(indcana, indctime, indclevel, indctimerange, indcdativarr, indcnetwork,indcattr)=&
269  ndi(indcana)*100.
270  end if
271  end do
272 
273  print *,">>>>>> NDI Volume <<<<<<"
274  call display (v7dqcspa%clima)
275 
276  IF (output_format == 'native') THEN
277  IF (output_file == '-') THEN ! stdout_unit does not work with unformatted
278  CALL l4f_category_log(category, l4f_info, 'trying /dev/stdout as stdout unit.')
279  output_file='/dev/stdout'
280  ENDIF
281  iun = getunit()
282  OPEN(iun, file=output_file, form='UNFORMATTED', access='STREAM')
283  CALL export(v7dqcspa%clima, unit=iun)
284  CLOSE(iun)
285  CALL delete(v7dqcspa%clima)
286 
287 #ifdef HAVE_DBALLE
288  ELSE IF (output_format == 'BUFR' .OR. output_format == 'CREX' .OR. output_format == 'dba') THEN
289  IF (output_format == 'BUFR' .OR. output_format == 'CREX') THEN
290  IF (output_file == '-') THEN
291  CALL l4f_category_log(category, l4f_info, 'trying /dev/stdout as stdout unit.')
292  output_file='/dev/stdout'
293  ENDIF
294  file=.true.
295 
296  ELSE IF (output_format == 'dba') THEN
297  file=.false.
298  ENDIF
299 
300  IF (output_template == '') output_template = 'generic'
301  ! check whether wipe=file is reasonable
302  CALL init(v7d_dba_out, filename=output_file, format=output_format, &
303  dsn=output_file, file=file, write=.true., wipe=file)
304 
305  v7d_dba_out%vol7d = v7dqcspa%clima
306  CALL init(v7dqcspa%clima) ! nullify without deallocating
307  CALL export(v7d_dba_out, template=output_template)
308  CALL delete(v7d_dba_out)
309 #endif
310 
311  end if
312 
313  stop
314 
315 end if
316 
317 !------------------------------------------------------------------------
318 ! read the namelist to define DSN
319 !------------------------------------------------------------------------
320 
321 open(10,file='qc.nml',status='old')
322 read(10,nml=odbc,iostat=io)
323 if ( io == 0 ) read(10,nml=odbcspa,iostat=io)
324 if ( io == 0 ) read(10,nml=switchspa,iostat=io)
325 if ( io == 0 ) read(10,nml=minmax,iostat=io)
326 if ( io == 0 ) read(10,nml=varlist,iostat=io)
327 
328 if (io /= 0 )then
329  call l4f_category_log(category,l4f_error,"Error reading namelist qc.nml")
330  call raise_error("Error reading namelist qc.nml")
331 end if
332 close(10)
333 
334 
335 !------------------------------------------------------------------------
336 ! Define what you want to QC
337 !------------------------------------------------------------------------
338 
339 nvar=count(c_e(var))
340 
341 if (nvar == 0) then
342  call l4f_category_log(category,l4f_error,"0 variables defined")
343  call raise_error()
344 end if
345  ! Definisco le date iniziale e finale
346 CALL init(ti, year=years, month=months, day=days, hour=hours)
347 CALL init(tf, year=yeare, month=monthe, day=daye, hour=houre)
348 print *,"time extreme"
349 call display(ti)
350 call display(tf)
351 
352  ! Define coordinate box
353 CALL init(coordmin,lat=lats,lon=lons)
354 CALL init(coordmax,lat=late,lon=lone)
355 
356 !call getval(coordmin,lon=lon,lat=lat)
357 print*,"lon lat minimum -> ",to_char(coordmin)
358 !call getval(coordmax,lon=lon,lat=lat)
359 print*,"lon lat maximum -> ",to_char(coordmax)
360 
361 !------------------------------------------------------------------------
362 call l4f_category_log(category,l4f_info,"QC on "//t2c(nvar)//" variables")
363 do i=1,nvar
364  call l4f_category_log(category,l4f_info,"QC on "//var(i)//" variable")
365 enddo
366 if (c_e(lons)) call l4f_category_log(category,l4f_info,"QC on "//t2c(lons)//" lon min value")
367 if (c_e(lone)) call l4f_category_log(category,l4f_info,"QC on "//t2c(lone)//" lon max value")
368 if (c_e(lats)) call l4f_category_log(category,l4f_info,"QC on "//t2c(lats)//" lat min value")
369 if (c_e(late)) call l4f_category_log(category,l4f_info,"QC on "//t2c(late)//" lat max value")
370 if (c_e(ti)) call l4f_category_log(category,l4f_info,"QC on "//t2c(ti)//" datetime min value")
371 if (c_e(tf)) call l4f_category_log(category,l4f_info,"QC on "//t2c(tf)//" datetime max value")
372 !------------------------------------------------------------------------
373 
374 time=ti
375 
376 ! aproximate to integer hours (not required reading date from namelist)
377 !time=ti+timedelta_new(minute=30)
378 !CALL getval(time,year, month, day, hour)
379 !call init(time, year, month, day, hour, minute=00, msec=00)
380 
381 !if (time < timei) time=time+timedelta_new(hour=1)
382 !timef=tf
383 !if (time > timef) time=timei
384 
385 #ifdef HAVE_LIBNCARG
386 if (doplot) then
387  call l4f_category_log(category,l4f_info,"start plot")
388  call init(plot,pstype='PS', orient='LANDSCAPE',color='COLOR',file="v7d_qcspa.ps")
389 end if
390 #endif
391 DO WHILE (time <= tf)
392  timei = time - timedelta_new(minute=30)
393  timef = time + timedelta_new(minute=30)
394  timeiqc = time - timedelta_new(minute=15)
395  timefqc = time + timedelta_new(minute=15)
396  call l4f_category_log(category,l4f_info,"elaborate from "//t2c(timeiqc)//" to "//t2c(timefqc))
397 
398  ! Chiamo il costruttore della classe vol7d_dballe per il mio oggetto in import
399  CALL init(v7ddballe,dsn=dsn,write=.true.,wipe=.false.,categoryappend="data-"//t2c(time))
400  call l4f_category_log(category,l4f_info,"start data import")
401 
402  CALL import(v7ddballe,var=var(:nvar),varkind=(/("r",i=1,nvar)/),&
403  anavar=(/"B07030"/),anavarkind=(/"r"/),&
404  attr=(/qcattrvarsbtables(1),qcattrvarsbtables(2),qcattrvarsbtables(4)/),attrkind=(/"b","b","b"/)&
405  ,timei=timei,timef=timef,coordmin=coordmin,coordmax=coordmax)
406 
407 ! call display(v7ddballe%vol7d)
408  call l4f_category_log(category,l4f_info,"end data import")
409  call l4f_category_log(category,l4f_info, "input N staz="//t2c(size(v7ddballe%vol7d%ana)))
410 
411  call l4f_category_log(category,l4f_info,"start peeling")
412 
413  !remove data invalidated and gross error only
414  !this change the behaviour of qcsummary_flag to ignore confidence att thresold
415  !qcpar=qcpartype(0_int_b,0_int_b,1_int_b)
416  qcpar%att=bmiss
417  !call vol7d_peeling(v7ddballe%vol7d,v7ddballe%data_id,keep_attr=(/qcattrvarsbtables(4)/),purgeana=.true.)
418  call vol7d_peeling(v7ddballe%vol7d,keep_attr=(/qcattrvarsbtables(4)/),purgeana=.true.)
419 ! call display(v7ddballe%vol7d)
420 
421  call l4f_category_log(category,l4f_info, "filtered N staz="//t2c(size(v7ddballe%vol7d%ana)))
422 
423  call l4f_category_log(category,l4f_info,"initialize QC")
424 
425  ! chiamiamo il "costruttore" per il Q.C.
426 
427 
428  call init(v7dqcspa,v7ddballe%vol7d,var(:nvar),timei=timeiqc,timef=timefqc, coordmin=coordmin, coordmax=coordmax,&
429  !data_id_in=v7ddballe%data_id, &
430  dsne=dsne, usere=usere, passworde=passworde,&
431  dsnspa=dsnspa, userspa=userspa, passwordspa=passwordspa,&
432  height2level=height2level, operation=operation,&
433  categoryappend="space"//t2c(time))
434 
435 ! print *,">>>>>> Clima Spatial Volume <<<<<<"
436 ! call display(v7dqcspa%clima)
437 
438  !print *,">>>>>> Pre Data Volume <<<<<<"
439  !call display(v7dqcspa%v7d)
440 
441  call alloc(v7dqcspa)
442 
443  ! spatial QC
444  !attention: do not exclude the first/last time so we check data two times
445  call l4f_category_log(category,l4f_info,"start spatial QC")
446  !call quaconspa(v7dqcspa,noborder=.true.,timemask= ( v7dqcspa%v7d%time >= timeiqc .and. v7dqcspa%v7d%time < timefqc ))
447  call quaconspa(v7dqcspa,timetollerance=timedelta_new(minute=15),noborder=.true.,&
448  timemask= ( v7dqcspa%v7d%time >= timeiqc .and. v7dqcspa%v7d%time <= timefqc ))
449  call l4f_category_log(category,l4f_info,"end spatial QC")
450 
451 #ifdef HAVE_LIBNCARG
452  if (doplot) then
453  call l4f_category_log(category,l4f_info,"start plot")
454  call plot_triangles(plot,v7dqcspa%co,v7dqcspa%tri,logo="Time: "//t2c(timeiqc)//" to "//t2c(timefqc))
455  call frame()
456  end if
457 #endif
458 
459 
460  ! prepare data_id to be recreated
461  !deallocate(v7ddballe%data_id)
462  !nullify(v7ddballe%data_id)
463 
464  if (v7dqcspa%operation == "run") then
465  call l4f_category_log(category,l4f_info,"start export data")
466  !print *,">>>>>> Post Data Volume <<<<<<"
467  !call display(v7ddballe%vol7d)
468 
469  ! data_id to use is the new one
470  !v7ddballe%data_id => v7dqcspa%data_id_out
471 
472  CALL export(v7ddballe,&
473  filter=dbafilter(datetimemin=dbadatetime(timeiqc),datetimemax=dbadatetime(timefqc)),&
474  attr_only=.true.)
475  !CALL export(v7ddballe)
476  call l4f_category_log(category,l4f_info,"end export data")
477  end if
478 
479  call delete(v7dqcspa)
480  ! data_id was allready deleted
481  !nullify(v7ddballe%data_id)
482  call delete(v7ddballe)
483 
484  time = time + timedelta_new(minute=30)
485 end do
486 
487 #ifdef HAVE_LIBNCARG
488  if (doplot) then
489  call delete(plot)
490  end if
491 #endif
492 
493 !close logger
494 call l4f_category_delete(category)
495 ier=l4f_fini()
496 
497 end program v7d_qcspa
Functions that return a trimmed CHARACTER representation of the input variable.
Export one or more geo_coordvect objects to a plain text file or to a file in ESRI/Shapefile format.
Import one or more geo_coordvect objects from a plain text file or for a file in ESRI/Shapefile forma...
Constructors for the two classes.
Represent geo_coord object in a pretty string.
Emit log message for a category with specific priority.
log4fortran destructor
Global log4fortran constructor.
Function to check whether a value is missing or not.
Allocazione di memoria.
Definition: modqcspa.F90:332
Add a new option of a specific type.
Compute a set of percentiles for a random variable.
This module defines usefull general purpose function and subroutine.
Classi per la gestione delle coordinate temporali.
class for import and export data from e to DB-All.e.
Classes for handling georeferenced sparse points in geographical corodinates.
classe per la gestione del logging
Definitions of constants and functions for working with missing values.
Utilities and defines for quality control.
Definition: modqc.F90:363
Controllo di qualità spaziale.
Definition: modqcspa.F90:279
Module for parsing command-line optons.
Module for basic statistical computations taking into account missing data.
Definition: simple_stat.f90:25
Classe per la gestione di un volume completo di dati osservati.
classe per import ed export di volumi da e in DB-All.e
filter to apply before ingest data

Generated with Doxygen.