libsim  Versione6.3.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
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 USE db_utils
42 
43 implicit none
44 
45 integer :: category,io,ier,i,iun,n,ind
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(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 the form user/password@dsn,' &
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_if_possible)
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  CALL parse_dba_access_info(output_file, dsn, user, password)
299  file=.false.
300  ENDIF
301 
302  IF (output_template == '') output_template = 'generic'
303  ! check whether wipe=file is reasonable
304  CALL init(v7d_dba_out, filename=output_file, format=output_format, &
305  dsn=dsn, user=user, password=password, file=file, write=.true., wipe=file)
306 
307  v7d_dba_out%vol7d = v7dqcspa%clima
308  CALL init(v7dqcspa%clima) ! nullify without deallocating
309  CALL export(v7d_dba_out, template=output_template)
310  CALL delete(v7d_dba_out)
311 #endif
312 
313  end if
314 
315  stop
316 
317 end if
318 
319 !------------------------------------------------------------------------
320 ! read the namelist to define DSN
321 !------------------------------------------------------------------------
322 
323 open(10,file='qc.nml',status='old')
324 read(10,nml=odbc,iostat=io)
325 if ( io == 0 ) read(10,nml=odbcspa,iostat=io)
326 if ( io == 0 ) read(10,nml=switchspa,iostat=io)
327 if ( io == 0 ) read(10,nml=minmax,iostat=io)
328 if ( io == 0 ) read(10,nml=varlist,iostat=io)
329 
330 if (io /= 0 )then
331  call l4f_category_log(category,l4f_error,"Error reading namelist qc.nml")
332  call raise_error("Error reading namelist qc.nml")
333 end if
334 close(10)
335 
336 
337 !------------------------------------------------------------------------
338 ! Define what you want to QC
339 !------------------------------------------------------------------------
340 
341 nvar=count(c_e(var))
342 
343 if (nvar == 0) then
344  call l4f_category_log(category,l4f_error,"0 variables defined")
345  call raise_error()
346 end if
347  ! Definisco le date iniziale e finale
348 CALL init(ti, year=years, month=months, day=days, hour=hours)
349 CALL init(tf, year=yeare, month=monthe, day=daye, hour=houre)
350 print *,"time extreme"
351 call display(ti)
352 call display(tf)
353 
354  ! Define coordinate box
355 CALL init(coordmin,lat=lats,lon=lons)
356 CALL init(coordmax,lat=late,lon=lone)
357 
358 !call getval(coordmin,lon=lon,lat=lat)
359 print*,"lon lat minimum -> ",to_char(coordmin)
360 !call getval(coordmax,lon=lon,lat=lat)
361 print*,"lon lat maximum -> ",to_char(coordmax)
362 
363 !------------------------------------------------------------------------
364 call l4f_category_log(category,l4f_info,"QC on "//t2c(nvar)//" variables")
365 do i=1,nvar
366  call l4f_category_log(category,l4f_info,"QC on "//var(i)//" variable")
367 enddo
368 if (c_e(lons)) call l4f_category_log(category,l4f_info,"QC on "//t2c(lons)//" lon min value")
369 if (c_e(lone)) call l4f_category_log(category,l4f_info,"QC on "//t2c(lone)//" lon max value")
370 if (c_e(lats)) call l4f_category_log(category,l4f_info,"QC on "//t2c(lats)//" lat min value")
371 if (c_e(late)) call l4f_category_log(category,l4f_info,"QC on "//t2c(late)//" lat max value")
372 if (c_e(ti)) call l4f_category_log(category,l4f_info,"QC on "//t2c(ti)//" datetime min value")
373 if (c_e(tf)) call l4f_category_log(category,l4f_info,"QC on "//t2c(tf)//" datetime max value")
374 !------------------------------------------------------------------------
375 
376 time=ti
377 
378 ! aproximate to integer hours (not required reading date from namelist)
379 !time=ti+timedelta_new(minute=30)
380 !CALL getval(time,year, month, day, hour)
381 !call init(time, year, month, day, hour, minute=00, msec=00)
382 
383 !if (time < timei) time=time+timedelta_new(hour=1)
384 !timef=tf
385 !if (time > timef) time=timei
386 
387 #ifdef HAVE_LIBNCARG
388 if (doplot) then
389  call l4f_category_log(category,l4f_info,"start plot")
390  call init(plot,pstype='PS', orient='LANDSCAPE',color='COLOR',file="v7d_qcspa.ps")
391 end if
392 #endif
393 DO WHILE (time <= tf)
394  timei = time - timedelta_new(minute=30)
395  timef = time + timedelta_new(minute=30)
396  timeiqc = time - timedelta_new(minute=15)
397  timefqc = time + timedelta_new(minute=15)
398  call l4f_category_log(category,l4f_info,"elaborate from "//t2c(timeiqc)//" to "//t2c(timefqc))
399 
400  ! Chiamo il costruttore della classe vol7d_dballe per il mio oggetto in import
401  CALL init(v7ddballe,dsn=dsn,user=user,password=password,write=.true.,wipe=.false.,categoryappend="data-"//t2c(time))
402  call l4f_category_log(category,l4f_info,"start data import")
403 
404  CALL import(v7ddballe,var=var(:nvar),varkind=(/("r",i=1,nvar)/),&
405  anavar=(/"B07030"/),anavarkind=(/"r"/),&
406  attr=(/qcattrvarsbtables(1),qcattrvarsbtables(2),qcattrvarsbtables(4)/),attrkind=(/"b","b","b"/)&
407  ,timei=timei,timef=timef,coordmin=coordmin,coordmax=coordmax)
408 
409 ! call display(v7ddballe%vol7d)
410  call l4f_category_log(category,l4f_info,"end data import")
411  call l4f_category_log(category,l4f_info, "input N staz="//t2c(size(v7ddballe%vol7d%ana)))
412 
413  call l4f_category_log(category,l4f_info,"start peeling")
414 
415  !remove data invalidated and gross error only
416  !this change the behaviour of qcsummary_flag to ignore confidence att thresold
417  !qcpar=qcpartype(0_int_b,0_int_b,1_int_b)
418  qcpar%att=bmiss
419  !call vol7d_peeling(v7ddballe%vol7d,v7ddballe%data_id,keep_attr=(/qcattrvarsbtables(4)/),purgeana=.true.)
420  call vol7d_peeling(v7ddballe%vol7d,keep_attr=(/qcattrvarsbtables(4)/),purgeana=.true.)
421 ! call display(v7ddballe%vol7d)
422 
423  call l4f_category_log(category,l4f_info, "filtered N staz="//t2c(size(v7ddballe%vol7d%ana)))
424 
425  call l4f_category_log(category,l4f_info,"initialize QC")
426 
427  ! chiamiamo il "costruttore" per il Q.C.
428 
429 
430  call init(v7dqcspa,v7ddballe%vol7d,var(:nvar),timei=timeiqc,timef=timefqc, coordmin=coordmin, coordmax=coordmax,&
431  !data_id_in=v7ddballe%data_id, &
432  dsne=dsne, usere=usere, passworde=passworde,&
433  dsnspa=dsnspa, userspa=userspa, passwordspa=passwordspa,&
434  height2level=height2level, operation=operation,&
435  categoryappend="space"//t2c(time))
436 
437 ! print *,">>>>>> Clima Spatial Volume <<<<<<"
438 ! call display(v7dqcspa%clima)
439 
440  !print *,">>>>>> Pre Data Volume <<<<<<"
441  !call display(v7dqcspa%v7d)
442 
443  call alloc(v7dqcspa)
444 
445  ! spatial QC
446  !attention: do not exclude the first/last time so we check data two times
447  call l4f_category_log(category,l4f_info,"start spatial QC")
448  !call quaconspa(v7dqcspa,noborder=.true.,timemask= ( v7dqcspa%v7d%time >= timeiqc .and. v7dqcspa%v7d%time < timefqc ))
449  call quaconspa(v7dqcspa,timetollerance=timedelta_new(minute=15),noborder=.true.,&
450  timemask= ( v7dqcspa%v7d%time >= timeiqc .and. v7dqcspa%v7d%time <= timefqc ))
451  call l4f_category_log(category,l4f_info,"end spatial QC")
452 
453 #ifdef HAVE_LIBNCARG
454  if (doplot) then
455  call l4f_category_log(category,l4f_info,"start plot")
456  call plot_triangles(plot,v7dqcspa%co,v7dqcspa%tri,logo="Time: "//t2c(timeiqc)//" to "//t2c(timefqc))
457  call frame()
458  end if
459 #endif
460 
461 
462  ! prepare data_id to be recreated
463  !deallocate(v7ddballe%data_id)
464  !nullify(v7ddballe%data_id)
465 
466  if (v7dqcspa%operation == "run") then
467  call l4f_category_log(category,l4f_info,"start export data")
468  !print *,">>>>>> Post Data Volume <<<<<<"
469  !call display(v7ddballe%vol7d)
470 
471  ! data_id to use is the new one
472  !v7ddballe%data_id => v7dqcspa%data_id_out
473 
474  CALL export(v7ddballe,&
475  filter=dbafilter(datetimemin=dbadatetime(timeiqc),datetimemax=dbadatetime(timefqc)),&
476  attr_only=.true.)
477  !CALL export(v7ddballe)
478  call l4f_category_log(category,l4f_info,"end export data")
479  end if
480 
481  call delete(v7dqcspa)
482  ! data_id was allready deleted
483  !nullify(v7ddballe%data_id)
484  call delete(v7ddballe)
485 
486  time = time + timedelta_new(minute=30)
487 end do
488 
489 #ifdef HAVE_LIBNCARG
490  if (doplot) then
491  call delete(plot)
492  end if
493 #endif
494 
495 !close logger
496 call l4f_category_delete(category)
497 ier=l4f_fini()
498 
499 end program v7d_qcspa

Generated with Doxygen.