libsim Versione 7.2.1
v7d_qctem.F90

Sample program for module qctem.

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
22program v7d_qctem
23
29USE modqc
32!use vol7d_dballeold_class
35use modqctem
36
37implicit none
38
39integer :: category,io,ier,i,iun,n,ind
40character(len=512):: a_name,output_format, output_template
41
42 !tipi derivati.
43TYPE(optionparser) :: opt
44TYPE(geo_coord) :: coordmin, coordmax
45TYPE(datetime) :: time, timeiqc, timefqc
46type(qctemtype) :: v7dqctem
47type(vol7d) :: v7dana
48
49integer, parameter :: maxvar=10
50character(len=6) :: var(maxvar)=cmiss ! variables to elaborate
51#ifdef HAVE_DBALLE
52type(vol7d_dballe) :: v7ddballe,v7ddballeana,v7d_dba_out
53character(len=512) :: dsn='test1',user='test',password=''
54character(len=512) :: dsne='test',usere='test',passworde=''
55character(len=512) :: dsntem='test',usertem='test',passwordtem=''
56#endif
57integer :: years=imiss,months=imiss,days=imiss,hours=imiss,yeare=imiss,monthe=imiss,daye=imiss,houre=imiss,nvar=0
58doubleprecision :: lons=dmiss,lats=dmiss,lone=dmiss,late=dmiss
59logical :: height2level=.false.,version
60CHARACTER(len=512) :: input_file, output_file
61INTEGER :: optind, optstatus, ninput
62CHARACTER(len=20) :: operation
63TYPE(ARRAYOF_REAL):: grad
64real :: val
65integer :: iostat
66REAL, 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.
67REAL, DIMENSION(size(perc_vals)-1) :: ndi
68REAL, DIMENSION(size(perc_vals)) :: nlimbins
69double precision, dimension(2) :: percentile
70TYPE(cyclicdatetime) :: cyclicdt
71type(vol7d_level) :: level,levelo
72type(vol7d_timerange) :: timerange,timerangeo
73type(vol7d_var) :: varia,variao
74CHARACTER(len=vol7d_ana_lenident) :: ident
75integer :: indcana,indctime,indclevel,indctimerange,indcdativarr,indcnetwork,indcattr,iana
76INTEGER,POINTER :: w_s(:), w_e(:)
77logical :: file
78
79#ifdef HAVE_DBALLE
80namelist /odbc/ dsn,user,password,dsne,usere,passworde
81namelist /odbctem/ dsntem,usertem,passwordtem ! namelist to define DSN
82#endif
83namelist /switchtem/ height2level
84namelist /minmax/ years,months,days,hours,lons,lats,yeare,monthe,daye,houre,lone,late
85namelist /varlist/ var
86
87!init log4fortran
88ier=l4f_init()
89
90! unique name from launcher
91call l4f_launcher(a_name,a_name_force="v7d_qctem")
92
93! set a_name
94category=l4f_category_get(a_name//".main")
95
96
97! define the option parser
98opt = 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
110CALL 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
117output_template = ''
118CALL 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
128CALL optionparser_add_help(opt, 'h', 'help', help='show an help message and exit')
129CALL optionparser_add(opt, ' ', 'version', version, help='show version and exit')
130
131! parse options and check for errors
132CALL optionparser_parse(opt, optind, optstatus)
133
134IF (optstatus == optionparser_help) THEN
135 CALL exit(0) ! generate a clean manpage
136ELSE IF (optstatus == optionparser_err) THEN
137 CALL l4f_category_log(category,l4f_error,'in command-line parameters')
138 CALL raise_fatal_error()
139ENDIF
140IF (version) THEN
141 WRITE(*,'(A,1X,A)')'v7d_qctem',version
142 CALL exit(0)
143ENDIF
144
145
146if (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()
151end if
152
153
154! check output format/template
155n = word_split(output_format, w_s, w_e, ':')
156IF (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:) = ' '
159ENDIF
160DEALLOCATE(w_s, w_e)
161
162if (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
322end if
323
324!------------------------------------------------------------------------
325! read the namelist to define DSN
326!------------------------------------------------------------------------
327
328open(10,file='qc.nml',status='old')
329read(10,nml=odbc,iostat=io)
330if ( io == 0 ) read(10,nml=odbctem,iostat=io)
331if ( io == 0 ) read(10,nml=switchtem,iostat=io)
332if ( io == 0 ) read(10,nml=minmax,iostat=io)
333if ( io == 0 ) read(10,nml=varlist,iostat=io)
334
335if (io /= 0 )then
336 call l4f_category_log(category,l4f_error,"Error reading namelist qc.nml")
337 call raise_error()
338end if
339close(10)
340
341
342!------------------------------------------------------------------------
343! Define what you want to QC
344!------------------------------------------------------------------------
345
346nvar=count(c_e(var))
347
348if (nvar == 0) then
349 call l4f_category_log(category,l4f_error,"0 variables defined")
350 call raise_error()
351end if
352 ! Definisco le date iniziale e finale
353CALL init(timeiqc, year=years, month=months, day=days, hour=hours)
354CALL init(timefqc, year=yeare, month=monthe, day=daye, hour=houre)
355!print *,"time extreme"
356call display(timeiqc)
357call display(timefqc)
358
359 ! Define coordinate box
360CALL init(coordmin,lat=lats,lon=lons)
361CALL init(coordmax,lat=late,lon=lone)
362
363!call getval(coordmin,lon=lon,lat=lat)
364print*,"lon lat minimum -> ",to_char(coordmin)
365!call getval(coordmax,lon=lon,lat=lat)
366print*,"lon lat maximum -> ",to_char(coordmax)
367
368!------------------------------------------------------------------------
369call l4f_category_log(category,l4f_info,"QC on "//t2c(nvar)//" variables")
370do i=1,nvar
371 call l4f_category_log(category,l4f_info,"QC on "//var(i)//" variable")
372enddo
373if (c_e(lons)) call l4f_category_log(category,l4f_info,"QC on "//t2c(lons)//" lon min value")
374if (c_e(lone)) call l4f_category_log(category,l4f_info,"QC on "//t2c(lone)//" lon max value")
375if (c_e(lats)) call l4f_category_log(category,l4f_info,"QC on "//t2c(lats)//" lat min value")
376if (c_e(late)) call l4f_category_log(category,l4f_info,"QC on "//t2c(late)//" lat max value")
377if (c_e(timeiqc)) call l4f_category_log(category,l4f_info,"QC on "//t2c(timeiqc)//" datetime min value")
378if (c_e(timefqc)) call l4f_category_log(category,l4f_info,"QC on "//t2c(timefqc)//" datetime max value")
379!------------------------------------------------------------------------
380
381
382CALL init(v7ddballeana,dsn=dsn,write=.false.,wipe=.false.,categoryappend="QCtarget-anaonly")
383call l4f_category_log(category,l4f_info,"start ana import")
384CALL import(v7ddballeana,var=var(:nvar),timei=timeiqc,timef=timefqc,coordmin=coordmin,coordmax=coordmax,anaonly=.true.)
385call vol7d_copy(v7ddballeana%vol7d,v7dana)
386call 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
395print *,"anagrafica:"
396call display(v7dana)
397
398do 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
477end do
478
479call delete(v7dana)
480
481!close logger
482call l4f_category_delete(category)
483ier=l4f_fini()
484
485end program v7d_qctem
Destructor for finalizing an array object.
Method for inserting elements of the array at a desired position.
Operatore di valore assoluto di un intervallo.
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: modqctem.F90:303
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:357
Controllo di qualità temporale.
Definition: modqctem.F90:256
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.