libsim Versione 7.2.0
v7d_qcspa.F90

Sample program for module qcspa.

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
22program v7d_qcspa
23
30use modqc
31use modqcspa
32!use vol7d_dballeold_class
37
38#ifdef HAVE_LIBNCARG
39USE ncar_plot_class
40#endif
41
42implicit none
43
44integer :: category,io,ier,i,iun,n,ind
45character(len=512):: a_name,output_format, output_template
46
47 !tipi derivati.
48TYPE(optionparser) :: opt
49TYPE(geo_coord) :: coordmin, coordmax
50TYPE(datetime) :: time,ti, tf, timei, timef, timeiqc, timefqc
51type(qcspatype) :: v7dqcspa
52type(vol7d_dballe) :: v7ddballe
53#ifdef HAVE_LIBNCARG
54type(ncar_plot) :: plot
55#endif
56
57integer, parameter :: maxvar=10
58character(len=6) :: var(maxvar)=cmiss ! variables to elaborate
59#ifdef HAVE_DBALLE
60character(len=512) :: dsn='test1',user='test',password=''
61character(len=512) :: dsne='test',usere='test',passworde=''
62character(len=512) :: dsnspa='test',userspa='test',passwordspa=''
63#endif
64integer :: years=imiss,months=imiss,days=imiss,hours=imiss,yeare=imiss,monthe=imiss,daye=imiss,houre=imiss,nvar=0
65doubleprecision :: lons=dmiss,lats=dmiss,lone=dmiss,late=dmiss
66integer :: year, month, day, hour
67logical :: height2level=.false.,doplot=.false.,version
68CHARACTER(len=512) :: input_file, output_file
69INTEGER :: optind, optstatus, ninput
70CHARACTER(len=20) :: operation
71TYPE(ARRAYOF_REAL):: grad
72real :: val
73integer :: iostat
74REAL, 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.
75REAL, DIMENSION(size(perc_vals)-1) :: ndi
76REAL, DIMENSION(size(perc_vals)) :: nlimbins
77double precision, dimension(2) :: percentile
78TYPE(cyclicdatetime) :: cyclicdt
79type(vol7d_level) :: level,levelo
80type(vol7d_timerange) :: timerange,timerangeo
81type(vol7d_var) :: varia,variao
82CHARACTER(len=vol7d_ana_lenident) :: ident
83integer :: indcana,indctime,indclevel,indctimerange,indcdativarr,indcnetwork,indcattr
84INTEGER,POINTER :: w_s(:), w_e(:)
85TYPE(vol7d_dballe) :: v7d_dba_out
86logical :: file
87integer :: status,grunit
88
89#ifdef HAVE_DBALLE
90namelist /odbc/ dsn,user,password,dsne,usere,passworde
91namelist /odbcspa/ dsnspa,userspa,passwordspa ! namelist to define DSN
92#endif
93namelist /switchspa/ height2level,doplot
94namelist /minmax/ years,months,days,hours,lons,lats,yeare,monthe,daye,houre,lone,late
95namelist /varlist/ var
96
97!init log4fortran
98ier=l4f_init()
99
100! unique name from launcher
101call l4f_launcher(a_name,a_name_force="v7d_qcspa")
102
103! set a_name
104category=l4f_category_get(a_name//".main")
105
106
107! define the option parser
108opt = 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
120CALL 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
127output_template = ''
128CALL 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
138CALL optionparser_add_help(opt, 'h', 'help', help='show an help message and exit')
139CALL optionparser_add(opt, ' ', 'version', version, help='show version and exit')
140
141! parse options and check for errors
142CALL optionparser_parse(opt, optind, optstatus)
143
144IF (optstatus == optionparser_help) THEN
145 CALL exit(0) ! generate a clean manpage
146ELSE IF (optstatus == optionparser_err) THEN
147 CALL l4f_category_log(category,l4f_error,'in command-line parameters')
148 CALL raise_fatal_error()
149ENDIF
150IF (version) THEN
151 WRITE(*,'(A,1X,A)')'v7d_qcspa',version
152 CALL exit(0)
153ENDIF
154
155
156if (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()
161end if
162
163
164! check output format/template
165n = word_split(output_format, w_s, w_e, ':')
166IF (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:) = ' '
169ENDIF
170DEALLOCATE(w_s, w_e)
171
172if (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
315end if
316
317!------------------------------------------------------------------------
318! read the namelist to define DSN
319!------------------------------------------------------------------------
320
321open(10,file='qc.nml',status='old')
322read(10,nml=odbc,iostat=io)
323if ( io == 0 ) read(10,nml=odbcspa,iostat=io)
324if ( io == 0 ) read(10,nml=switchspa,iostat=io)
325if ( io == 0 ) read(10,nml=minmax,iostat=io)
326if ( io == 0 ) read(10,nml=varlist,iostat=io)
327
328if (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")
331end if
332close(10)
333
334
335!------------------------------------------------------------------------
336! Define what you want to QC
337!------------------------------------------------------------------------
338
339nvar=count(c_e(var))
340
341if (nvar == 0) then
342 call l4f_category_log(category,l4f_error,"0 variables defined")
343 call raise_error()
344end if
345 ! Definisco le date iniziale e finale
346CALL init(ti, year=years, month=months, day=days, hour=hours)
347CALL init(tf, year=yeare, month=monthe, day=daye, hour=houre)
348print *,"time extreme"
349call display(ti)
350call display(tf)
351
352 ! Define coordinate box
353CALL init(coordmin,lat=lats,lon=lons)
354CALL init(coordmax,lat=late,lon=lone)
355
356!call getval(coordmin,lon=lon,lat=lat)
357print*,"lon lat minimum -> ",to_char(coordmin)
358!call getval(coordmax,lon=lon,lat=lat)
359print*,"lon lat maximum -> ",to_char(coordmax)
360
361!------------------------------------------------------------------------
362call l4f_category_log(category,l4f_info,"QC on "//t2c(nvar)//" variables")
363do i=1,nvar
364 call l4f_category_log(category,l4f_info,"QC on "//var(i)//" variable")
365enddo
366if (c_e(lons)) call l4f_category_log(category,l4f_info,"QC on "//t2c(lons)//" lon min value")
367if (c_e(lone)) call l4f_category_log(category,l4f_info,"QC on "//t2c(lone)//" lon max value")
368if (c_e(lats)) call l4f_category_log(category,l4f_info,"QC on "//t2c(lats)//" lat min value")
369if (c_e(late)) call l4f_category_log(category,l4f_info,"QC on "//t2c(late)//" lat max value")
370if (c_e(ti)) call l4f_category_log(category,l4f_info,"QC on "//t2c(ti)//" datetime min value")
371if (c_e(tf)) call l4f_category_log(category,l4f_info,"QC on "//t2c(tf)//" datetime max value")
372!------------------------------------------------------------------------
373
374time=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
386if (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")
389end if
390#endif
391DO 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)
485end do
486
487#ifdef HAVE_LIBNCARG
488 if (doplot) then
489 call delete(plot)
490 end if
491#endif
492
493!close logger
494call l4f_category_delete(category)
495ier=l4f_fini()
496
497end program v7d_qcspa
Destructor for finalizing an array object.
Method for inserting elements of the array at a desired position.
Costruttori per le classi datetime e timedelta.
Functions that return a trimmed CHARACTER representation of the input variable.
Restituiscono il valore dell'oggetto in forma di stringa stampabile.
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...
Emit log message for a category with specific priority.
log4fortran destructor
Global log4fortran constructor.
Allocazione di memoria.
Definition: modqcspa.F90:326
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:357
Controllo di qualità spaziale.
Definition: modqcspa.F90:273
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.