libsim  Versione7.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

Generated with Doxygen.