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

Generated with Doxygen.