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