libsim  Versione7.2.1
v7d_qctem.F90

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

Generated with Doxygen.