libsim  Versione 7.1.6
v7d_qctem.F90

Sample program for module qctem

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_qctem
23 
24 use log4fortran
26 USE simple_stat
29 USE modqc
32 !use vol7d_dballeold_class
34 USE vol7d_class
35 use modqctem
36 
37 implicit none
38 
39 integer :: category,io,ier,i,iun,n,ind
40 character(len=512):: a_name,output_format, output_template
41 
42  !tipi derivati.
43 TYPE(optionparser) :: opt
44 TYPE(geo_coord) :: coordmin, coordmax
45 TYPE(datetime) :: time, timeiqc, timefqc
46 type(qctemtype) :: v7dqctem
47 type(vol7d) :: v7dana
48 
49 integer, parameter :: maxvar=10
50 character(len=6) :: var(maxvar)=cmiss ! variables to elaborate
51 #ifdef HAVE_DBALLE
52 type(vol7d_dballe) :: v7ddballe,v7ddballeana,v7d_dba_out
53 character(len=512) :: dsn='test1',user='test',password=''
54 character(len=512) :: dsne='test',usere='test',passworde=''
55 character(len=512) :: dsntem='test',usertem='test',passwordtem=''
56 #endif
57 integer :: years=imiss,months=imiss,days=imiss,hours=imiss,yeare=imiss,monthe=imiss,daye=imiss,houre=imiss,nvar=0
58 doubleprecision :: lons=dmiss,lats=dmiss,lone=dmiss,late=dmiss
59 logical :: height2level=.false.,version
60 CHARACTER(len=512) :: input_file, output_file
61 INTEGER :: optind, optstatus, ninput
62 CHARACTER(len=20) :: operation
63 TYPE(ARRAYOF_REAL):: grad
64 real :: val
65 integer :: iostat
66 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.
67 REAL, DIMENSION(size(perc_vals)-1) :: ndi
68 REAL, DIMENSION(size(perc_vals)) :: nlimbins
69 double precision, dimension(2) :: percentile
70 TYPE(cyclicdatetime) :: cyclicdt
71 type(vol7d_level) :: level,levelo
72 type(vol7d_timerange) :: timerange,timerangeo
73 type(vol7d_var) :: varia,variao
74 CHARACTER(len=vol7d_ana_lenident) :: ident
75 integer :: indcana,indctime,indclevel,indctimerange,indcdativarr,indcnetwork,indcattr,iana
76 INTEGER,POINTER :: w_s(:), w_e(:)
77 logical :: file
78 
79 #ifdef HAVE_DBALLE
80 namelist /odbc/ dsn,user,password,dsne,usere,passworde
81 namelist /odbctem/ dsntem,usertem,passwordtem ! namelist to define DSN
82 #endif
83 namelist /switchtem/ height2level
84 namelist /minmax/ years,months,days,hours,lons,lats,yeare,monthe,daye,houre,lone,late
85 namelist /varlist/ var
86 
87 !init log4fortran
88 ier=l4f_init()
89 
90 ! unique name from launcher
91 call l4f_launcher(a_name,a_name_force="v7d_qctem")
92 
93 ! set a_name
94 category=l4f_category_get(a_name//".main")
95 
96 
97 ! define the option parser
98 opt = optionparser_new(description_msg= &
99  'Spatial quality control: compute gradient; compute NDI from gradient; apply quality control', &
100  usage_msg='Usage: v7d_qctem [options] [inputfile1] [inputfile2...] [outputfile] \n&
101  &If input-format is of file type, inputfile ''-'' indicates stdin,'&
102 #ifdef HAVE_DBALLE
103  //'if output-format is of database type, outputfile specifies &
104  &database access info in YRI form,' &
105 #endif
106  //'if empty or ''-'', a suitable default is used. &
107 &')
108 
109 ! options for defining input
110 CALL optionparser_add(opt, ' ', 'operation', operation, "run", help= &
111  'operation to execute: ''gradient'' compute gradient and write on files; '&
112  //'''ndi'' compute NDI from gradient;' &
113  //'''run'' apply quality control ')
114 
115 
116 ! options for defining output
117 output_template = ''
118 CALL optionparser_add(opt, ' ', 'output-format', output_format, 'native', help= &
119  'format of output file for "ndi" operation only, in the form ''name[:template]''; ''native'' for vol7d &
120  &native binary format (no template to be specified)'&
121 #ifdef HAVE_DBALLE
122  //'; ''BUFR'' and ''CREX'' for corresponding formats, with template as an alias like ''synop'', ''metar'', &
123  &''temp'', ''generic'', empty for ''generic'''&
124 #endif
125 )
126 
127 ! help options
128 CALL optionparser_add_help(opt, 'h', 'help', help='show an help message and exit')
129 CALL optionparser_add(opt, ' ', 'version', version, help='show version and exit')
130 
131 ! parse options and check for errors
132 CALL optionparser_parse(opt, optind, optstatus)
133 
134 IF (optstatus == optionparser_help) THEN
135  CALL exit(0) ! generate a clean manpage
136 ELSE IF (optstatus == optionparser_err) THEN
137  CALL l4f_category_log(category,l4f_error,'in command-line parameters')
138  CALL raise_fatal_error()
139 ENDIF
140 IF (version) THEN
141  WRITE(*,'(A,1X,A)')'v7d_qctem',version
142  CALL exit(0)
143 ENDIF
144 
145 
146 if (operation /= "gradient" .and. operation /= "ndi" .and. operation /= "run") then
147  CALL optionparser_printhelp(opt)
148  CALL l4f_category_log(category, l4f_error, &
149  'argument to --operation is wrong')
150  CALL raise_fatal_error()
151 end if
152 
153 
154 ! check output format/template
155 n = word_split(output_format, w_s, w_e, ':')
156 IF (n >= 2) THEN ! set output template if present
157  output_template = output_format(w_s(2):w_e(2))
158  output_format(w_e(1)+1:) = ' '
159 ENDIF
160 DEALLOCATE(w_s, w_e)
161 
162 if (operation == "ndi") then
163 
164  ! check input/output files
165  i = iargc() - optind
166  IF (i < 0) THEN
167  CALL optionparser_printhelp(opt)
168  CALL l4f_category_log(category,l4f_error,'input file missing')
169  CALL raise_fatal_error()
170  ELSE IF (i < 1) THEN
171  CALL optionparser_printhelp(opt)
172  CALL l4f_category_log(category,l4f_error,'output file missing')
173  CALL raise_fatal_error()
174  ENDIF
175  CALL getarg(iargc(), output_file)
176 
177  call init(levelo)
178  call init(timerangeo)
179  call init (variao)
180 
181  do ninput = optind, iargc()-1
182  call getarg(ninput, input_file)
183 
184  CALL l4f_category_log(category,l4f_info,"open file:"//t2c(input_file))
185  open (10,file=input_file,status="old")
186  read (10,*) level,timerange,varia
187 
188  if (c_e(levelo)) then
189  if ( level /= levelo .or. timerange /= timerangeo .or. varia /= variao ) then
190  call l4f_category_log(category,l4f_error,"Error reading grad files: file are incoerent")
191  call raise_error("")
192  end if
193  else
194  levelo = level
195  timerangeo = timerange
196  variao = varia
197  end if
198 
199  do while (.true.)
200  read (10,*,iostat=iostat) val
201  if (iostat /= 0) exit
202  if (val /= 0.) call insert(grad,val)
203  end do
204 
205  close(10)
206  end do
207 
208  call delete(v7dqctem%clima)
209  CALL init(v7dqctem%clima, time_definition=0)
210  call vol7d_alloc(v7dqctem%clima,nana=size(perc_vals)-1, &
211  nlevel=1, ntimerange=1, &
212  ndativarr=1, nnetwork=2,ntime=1,ndativarattrr=1,ndatiattrr=1)
213 
214  call vol7d_alloc_vol(v7dqctem%clima,inivol=.true.)
215 
216  v7dqctem%clima%level=level
217  v7dqctem%clima%timerange=timerange
218  v7dqctem%clima%dativar%r(1)=varia
219  v7dqctem%clima%dativarattr%r(1)=varia
220  call init(v7dqctem%clima%datiattr%r(1), btable="*B33209") ! NDI order number
221  cyclicdt = cyclicdatetime_new(chardate="/////////") !TMMGGhhmm
222  v7dqctem%clima%time(1)=cyclicdatetime_to_conventional(cyclicdt)
223 
224  indctime=1
225  indclevel=1
226  indcattr=1
227  indcdativarr=1
228  indctimerange=1
229 
230 
231  indcnetwork=1
232  call init(v7dqctem%clima%network(indcnetwork),name="qctemsndi")
233 
234  CALL l4f_category_log(category,l4f_info,"compute NDI spike")
235 
236  CALL l4f_category_log(category,l4f_info,"compute percentile to remove the tails")
237  percentile = stat_percentile(abs(pack(grad%array(:grad%arraysize),mask=grad%array(:grad%arraysize) >= 0.)),(/10.,90./))
238  !print *,percentile
239 
240  call normalizeddensityindex(abs(pack(grad%array(:grad%arraysize),&
241  mask=(abs(grad%array(:grad%arraysize)) < percentile(2) .and. grad%array(:grad%arraysize) >= 0. ))), perc_vals, ndi, nlimbins)
242 
243  do indcana=1,size(perc_vals)-1
244  write(ident,'("#",i2.2,2i3.3)')0,0,nint(perc_vals(indcana))
245  call init(v7dqctem%clima%ana(indcana),ident=ident,lat=0d0,lon=0d0)
246  if (c_e(nlimbins(indcana)).and.c_e(ndi(indcana)))then
247  ind=index_c(tem_btable,varia%btable)
248  v7dqctem%clima%voldatir(indcana, indctime, indclevel, indctimerange, indcdativarr, indcnetwork)=&
249  nlimbins(indcana)*tem_a(ind) + tem_b(ind)
250  v7dqctem%clima%voldatiattrr(indcana, indctime, indclevel, indctimerange, indcdativarr, indcnetwork,indcattr)=&
251  ndi(indcana)*100.
252  end if
253  end do
254 
255  indcnetwork=2
256  call init(v7dqctem%clima%network(indcnetwork),name="qctemgndi")
257 
258  CALL l4f_category_log(category,l4f_info,"compute NDI gap")
259  CALL l4f_category_log(category,l4f_info,"compute percentile to remove the tails")
260  percentile = stat_percentile(abs(pack(grad%array(:grad%arraysize),mask=grad%array(:grad%arraysize) < 0.)),(/10.,90./))
261 
262  call normalizeddensityindex(abs(pack(grad%array(:grad%arraysize),&
263  mask=(abs(grad%array(:grad%arraysize)) < percentile(2) .and. grad%array(:grad%arraysize) < 0. ))), perc_vals, ndi, nlimbins)
264 
265  do indcana=1,size(perc_vals)-1
266  write(ident,'("#",i2.2,2i3.3)')0,0,nint(perc_vals(indcana))
267  call init(v7dqctem%clima%ana(indcana),ident=ident,lat=0d0,lon=0d0)
268  if (c_e(nlimbins(indcana)).and.c_e(ndi(indcana)))then
269  ind=index_c(tem_btable,varia%btable)
270  v7dqctem%clima%voldatir(indcana, indctime, indclevel, indctimerange, indcdativarr, indcnetwork)=&
271  nlimbins(indcana)*tem_a(ind) + tem_b(ind)
272  v7dqctem%clima%voldatiattrr(indcana, indctime, indclevel, indctimerange, indcdativarr, indcnetwork,indcattr)=&
273  ndi(indcana)*100.
274  end if
275  end do
276 
277  call delete(grad)
278  call display(v7dqctem%clima)
279 
280 ! print *,">>>>>> Clima Temporal Volume <<<<<<"
281 ! call display (v7dqctem%clima)
282 
283  IF (output_format == 'native') THEN
284  IF (output_file == '-') THEN ! stdout_unit does not work with unformatted
285  CALL l4f_category_log(category, l4f_info, 'trying /dev/stdout as stdout unit.')
286  output_file='/dev/stdout'
287  ENDIF
288  iun = getunit()
289  OPEN(iun, file=output_file, form='UNFORMATTED', access='STREAM')
290  CALL export(v7dqctem%clima, unit=iun)
291  CLOSE(iun)
292  CALL delete(v7dqctem%clima)
293 
294 #ifdef HAVE_DBALLE
295  ELSE IF (output_format == 'BUFR' .OR. output_format == 'CREX' .OR. output_format == 'dba') THEN
296  IF (output_format == 'BUFR' .OR. output_format == 'CREX') THEN
297  IF (output_file == '-') THEN
298  CALL l4f_category_log(category, l4f_info, 'trying /dev/stdout as stdout unit.')
299  output_file='/dev/stdout'
300  ENDIF
301  file=.true.
302 
303  ELSE IF (output_format == 'dba') THEN
304  file=.false.
305  ENDIF
306 
307  IF (output_template == '') output_template = 'generic'
308  ! check whether wipe=file is reasonable
309  CALL init(v7d_dba_out, filename=output_file, format=output_format, &
310  dsn=output_file, file=file, write=.true., wipe=file)
311 
312  v7d_dba_out%vol7d = v7dqctem%clima
313  CALL init(v7dqctem%clima) ! nullify without deallocating
314  CALL export(v7d_dba_out, template=output_template)
315  CALL delete(v7d_dba_out)
316 #endif
317 
318  end if
319 
320  stop
321 
322 end if
323 
324 !------------------------------------------------------------------------
325 ! read the namelist to define DSN
326 !------------------------------------------------------------------------
327 
328 open(10,file='qc.nml',status='old')
329 read(10,nml=odbc,iostat=io)
330 if ( io == 0 ) read(10,nml=odbctem,iostat=io)
331 if ( io == 0 ) read(10,nml=switchtem,iostat=io)
332 if ( io == 0 ) read(10,nml=minmax,iostat=io)
333 if ( io == 0 ) read(10,nml=varlist,iostat=io)
334 
335 if (io /= 0 )then
336  call l4f_category_log(category,l4f_error,"Error reading namelist qc.nml")
337  call raise_error()
338 end if
339 close(10)
340 
341 
342 !------------------------------------------------------------------------
343 ! Define what you want to QC
344 !------------------------------------------------------------------------
345 
346 nvar=count(c_e(var))
347 
348 if (nvar == 0) then
349  call l4f_category_log(category,l4f_error,"0 variables defined")
350  call raise_error()
351 end if
352  ! Definisco le date iniziale e finale
353 CALL init(timeiqc, year=years, month=months, day=days, hour=hours)
354 CALL init(timefqc, year=yeare, month=monthe, day=daye, hour=houre)
355 !print *,"time extreme"
356 call display(timeiqc)
357 call display(timefqc)
358 
359  ! Define coordinate box
360 CALL init(coordmin,lat=lats,lon=lons)
361 CALL init(coordmax,lat=late,lon=lone)
362 
363 !call getval(coordmin,lon=lon,lat=lat)
364 print*,"lon lat minimum -> ",to_char(coordmin)
365 !call getval(coordmax,lon=lon,lat=lat)
366 print*,"lon lat maximum -> ",to_char(coordmax)
367 
368 !------------------------------------------------------------------------
369 call l4f_category_log(category,l4f_info,"QC on "//t2c(nvar)//" variables")
370 do i=1,nvar
371  call l4f_category_log(category,l4f_info,"QC on "//var(i)//" variable")
372 enddo
373 if (c_e(lons)) call l4f_category_log(category,l4f_info,"QC on "//t2c(lons)//" lon min value")
374 if (c_e(lone)) call l4f_category_log(category,l4f_info,"QC on "//t2c(lone)//" lon max value")
375 if (c_e(lats)) call l4f_category_log(category,l4f_info,"QC on "//t2c(lats)//" lat min value")
376 if (c_e(late)) call l4f_category_log(category,l4f_info,"QC on "//t2c(late)//" lat max value")
377 if (c_e(timeiqc)) call l4f_category_log(category,l4f_info,"QC on "//t2c(timeiqc)//" datetime min value")
378 if (c_e(timefqc)) call l4f_category_log(category,l4f_info,"QC on "//t2c(timefqc)//" datetime max value")
379 !------------------------------------------------------------------------
380 
381 
382 CALL init(v7ddballeana,dsn=dsn,write=.false.,wipe=.false.,categoryappend="QCtarget-anaonly")
383 call l4f_category_log(category,l4f_info,"start ana import")
384 CALL import(v7ddballeana,var=var(:nvar),timei=timeiqc,timef=timefqc,coordmin=coordmin,coordmax=coordmax,anaonly=.true.)
385 call vol7d_copy(v7ddballeana%vol7d,v7dana)
386 call delete(v7ddballeana)
387 
388 !!$ an alternative to copy is:
389 
390 !!$idbhandle=v7ddballeana%idbhandle
391 !!$CALL init(v7ddballe,dsn=dsn,user=user,password=password,write=.true.,wipe=.false.,categoryappend="QCtarget-"//t2c(time),&
392 !!$ idbhandle=idbhandle)
393 !!$call delete(v7ddballeana,preserveidbhandle=.true.)
394 
395 print *,"anagrafica:"
396 call display(v7dana)
397 
398 do iana=1, size(v7dana%ana)
399 
400  call l4f_category_log(category,l4f_info,"elaborate station "//trim(to_char(v7dana%ana(iana))))
401 
402  ! Chiamo il costruttore della classe vol7d_dballe per il mio oggetto in import
403  CALL init(v7ddballe,dsn=dsn,write=.true.,wipe=.false.,categoryappend="QCtarget-"//t2c(time))
404  call l4f_category_log(category,l4f_info,"start data import")
405 
406  CALL import(v7ddballe,var=var(:nvar),varkind=(/("r",i=1,nvar)/),&
407  anavar=(/"B07030"/),anavarkind=(/"r"/),&
408  attr=(/qcattrvarsbtables(1),qcattrvarsbtables(2),qcattrvarsbtables(3)/),attrkind=(/"b","b","b"/)&
409  ,timei=timeiqc,timef=timefqc, ana=v7dana%ana(iana))
410 
411  ! this depend on witch version of dballe is in use: with dballe 6.6 comment this out
412  if (size(v7ddballe%vol7d%time)==0) then
413  CALL l4f_category_log(category,l4f_info,'skip empty volume: you are not using dballe >= 6.6 !')
414  call delete(v7ddballe)
415  cycle
416  end if
417 
418  print *,"input volume"
419  call display(v7ddballe%vol7d)
420  call l4f_category_log(category,l4f_info,"end data import")
421  call l4f_category_log(category,l4f_info, "input N staz="//t2c(size(v7ddballe%vol7d%ana)))
422 
423  call l4f_category_log(category,l4f_info,"start peeling")
424 
425  !remove data invalidated and gross error only
426  !qcpar=qcpartype(0_int_b,0_int_b,0_int_b)
427  qcpar%att=bmiss
428  !call vol7d_peeling(v7ddballe%vol7d,v7ddballe%data_id,keep_attr=(/qcattrvarsbtables(3)/),purgeana=.true.)
429  call vol7d_peeling(v7ddballe%vol7d,keep_attr=(/qcattrvarsbtables(3)/),purgeana=.true.)
430  !call display(v7ddballe%vol7d)
431 
432  call l4f_category_log(category,l4f_info, "filtered N time="//t2c(size(v7ddballe%vol7d%time)))
433 
434  call l4f_category_log(category,l4f_info,"initialize QC")
435 
436  ! chiamiamo il "costruttore" per il Q.C.
437 
438  call init(v7dqctem,v7ddballe%vol7d,var(:nvar),timei=timeiqc,timef=timefqc, &
439  coordmin=v7dana%ana(iana)%coord, coordmax=v7dana%ana(iana)%coord,&
440 ! data_id_in=v7ddballe%data_id, &
441  dsne=dsne, usere=usere, passworde=passworde,&
442  dsntem=dsntem, usertem=usertem, passwordtem=passwordtem,&
443  height2level=height2level, operation=operation,&
444  categoryappend="temporal")
445 
446 ! print *,">>>>>> Clima Temporal Volume <<<<<<"
447 ! call display(v7dqcspa%clima)
448 
449  call alloc(v7dqctem)
450 
451  ! temporal QC
452  call l4f_category_log(category,l4f_info,"start temporal QC")
453  call quacontem(v7dqctem,timemask= ( v7dqctem%v7d%time >= timeiqc .and. v7dqctem%v7d%time <= timefqc ))
454  call l4f_category_log(category,l4f_info,"end temporal QC")
455 
456  ! prepare data_id to be recreated
457  !deallocate(v7ddballe%data_id)
458  !nullify(v7ddballe%data_id)
459 
460  if (v7dqctem%operation == "run") then
461  call l4f_category_log(category,l4f_info,"start export data")
462  print *,"output volume"
463  call display(v7ddballe%vol7d)
464 
465  ! data_id to use is the new one
466  !v7ddballe%data_id => v7dqctem%data_id_out
467  !CALL export(v7ddballe,attr_only=.true.)
468  CALL export(v7ddballe, attr_only=.true.)
469  call l4f_category_log(category,l4f_info,"end export data")
470  end if
471 
472  call delete(v7dqctem)
473  ! data_id was allready deleted
474  !nullify(v7ddballe%data_id)
475  call delete(v7ddballe)
476 
477 end do
478 
479 call delete(v7dana)
480 
481 !close logger
482 call l4f_category_delete(category)
483 ier=l4f_fini()
484 
485 end program v7d_qctem
Operatore di valore assoluto di un intervallo.
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: modqctem.F90:309
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.
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à temporale.
Definition: modqctem.F90:262
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

Generated with Doxygen.