libsim  Versione8.0.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
25 USE missing_values
26 USE simple_stat
27 USE geo_coord_class
28 USE datetime_class
29 USE modqc
30 USE optionparser_class
31 use array_utilities
32 !use vol7d_dballeold_class
33 use vol7d_dballe_class
34 USE vol7d_class
35 use modqctem
36 
37 implicit none
38 
39 integer :: io,ier,i,iun,n,ind
40 TYPE(l4f_handle) :: category
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_handle(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 YRI form,' &
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')
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  file=.false.
306  ENDIF
307 
308  IF (output_template == '') output_template = 'generic'
309  ! check whether wipe=file is reasonable
310  CALL init(v7d_dba_out, filename=output_file, format=output_format, &
311  dsn=output_file, file=file, write=.true., wipe=file)
312 
313  v7d_dba_out%vol7d = v7dqctem%clima
314  CALL init(v7dqctem%clima) ! nullify without deallocating
315  CALL export(v7d_dba_out, template=output_template)
316  CALL delete(v7d_dba_out)
317 #endif
318 
319  end if
320 
321  stop
322 
323 end if
324 
325 !------------------------------------------------------------------------
326 ! read the namelist to define DSN
327 !------------------------------------------------------------------------
328 
329 open(10,file='qc.nml',status='old')
330 read(10,nml=odbc,iostat=io)
331 if ( io == 0 ) read(10,nml=odbctem,iostat=io)
332 if ( io == 0 ) read(10,nml=switchtem,iostat=io)
333 if ( io == 0 ) read(10,nml=minmax,iostat=io)
334 if ( io == 0 ) read(10,nml=varlist,iostat=io)
335 
336 if (io /= 0 )then
337  call l4f_category_log(category,l4f_error,"Error reading namelist qc.nml")
338  call raise_error()
339 end if
340 close(10)
341 
342 
343 !------------------------------------------------------------------------
344 ! Define what you want to QC
345 !------------------------------------------------------------------------
346 
347 nvar=count(c_e(var))
348 
349 if (nvar == 0) then
350  call l4f_category_log(category,l4f_error,"0 variables defined")
351  call raise_error()
352 end if
353  ! Definisco le date iniziale e finale
354 CALL init(timeiqc, year=years, month=months, day=days, hour=hours)
355 CALL init(timefqc, year=yeare, month=monthe, day=daye, hour=houre)
356 !print *,"time extreme"
357 call display(timeiqc)
358 call display(timefqc)
359 
360  ! Define coordinate box
361 CALL init(coordmin,lat=lats,lon=lons)
362 CALL init(coordmax,lat=late,lon=lone)
363 
364 !call getval(coordmin,lon=lon,lat=lat)
365 print*,"lon lat minimum -> ",to_char(coordmin)
366 !call getval(coordmax,lon=lon,lat=lat)
367 print*,"lon lat maximum -> ",to_char(coordmax)
368 
369 !------------------------------------------------------------------------
370 call l4f_category_log(category,l4f_info,"QC on "//t2c(nvar)//" variables")
371 do i=1,nvar
372  call l4f_category_log(category,l4f_info,"QC on "//var(i)//" variable")
373 enddo
374 if (c_e(lons)) call l4f_category_log(category,l4f_info,"QC on "//t2c(lons)//" lon min value")
375 if (c_e(lone)) call l4f_category_log(category,l4f_info,"QC on "//t2c(lone)//" lon max value")
376 if (c_e(lats)) call l4f_category_log(category,l4f_info,"QC on "//t2c(lats)//" lat min value")
377 if (c_e(late)) call l4f_category_log(category,l4f_info,"QC on "//t2c(late)//" lat max value")
378 if (c_e(timeiqc)) call l4f_category_log(category,l4f_info,"QC on "//t2c(timeiqc)//" datetime min value")
379 if (c_e(timefqc)) call l4f_category_log(category,l4f_info,"QC on "//t2c(timefqc)//" datetime max value")
380 !------------------------------------------------------------------------
381 
382 
383 CALL init(v7ddballeana,dsn=dsn,write=.false.,wipe=.false.,categoryappend="QCtarget-anaonly")
384 call l4f_category_log(category,l4f_info,"start ana import")
385 CALL import(v7ddballeana,var=var(:nvar),timei=timeiqc,timef=timefqc,coordmin=coordmin,coordmax=coordmax,anaonly=.true.)
386 call vol7d_copy(v7ddballeana%vol7d,v7dana)
387 call delete(v7ddballeana)
388 
389 !!$ an alternative to copy is:
390 
391 !!$idbhandle=v7ddballeana%idbhandle
392 !!$CALL init(v7ddballe,dsn=dsn,user=user,password=password,write=.true.,wipe=.false.,categoryappend="QCtarget-"//t2c(time),&
393 !!$ idbhandle=idbhandle)
394 !!$call delete(v7ddballeana,preserveidbhandle=.true.)
395 
396 print *,"anagrafica:"
397 call display(v7dana)
398 
399 do iana=1, size(v7dana%ana)
400 
401  call l4f_category_log(category,l4f_info,"elaborate station "//trim(to_char(v7dana%ana(iana))))
402 
403  ! Chiamo il costruttore della classe vol7d_dballe per il mio oggetto in import
404  CALL init(v7ddballe,dsn=dsn,write=.true.,wipe=.false.,categoryappend="QCtarget-"//t2c(time))
405  call l4f_category_log(category,l4f_info,"start data import")
406 
407  CALL import(v7ddballe,var=var(:nvar),varkind=(/("r",i=1,nvar)/),&
408  anavar=(/"B07030"/),anavarkind=(/"r"/),&
409  attr=(/qcattrvarsbtables(1),qcattrvarsbtables(2),qcattrvarsbtables(3)/),attrkind=(/"b","b","b"/)&
410  ,timei=timeiqc,timef=timefqc, ana=v7dana%ana(iana))
411 
412  ! this depend on witch version of dballe is in use: with dballe 6.6 comment this out
413  if (size(v7ddballe%vol7d%time)==0) then
414  CALL l4f_category_log(category,l4f_info,'skip empty volume: you are not using dballe >= 6.6 !')
415  call delete(v7ddballe)
416  cycle
417  end if
418 
419  print *,"input volume"
420  call display(v7ddballe%vol7d)
421  call l4f_category_log(category,l4f_info,"end data import")
422  call l4f_category_log(category,l4f_info, "input N staz="//t2c(size(v7ddballe%vol7d%ana)))
423 
424  call l4f_category_log(category,l4f_info,"start peeling")
425 
426  !remove data invalidated and gross error only
427  !qcpar=qcpartype(0_int_b,0_int_b,0_int_b)
428  qcpar%att=bmiss
429  !call vol7d_peeling(v7ddballe%vol7d,v7ddballe%data_id,keep_attr=(/qcattrvarsbtables(3)/),purgeana=.true.)
430  call vol7d_peeling(v7ddballe%vol7d,keep_attr=(/qcattrvarsbtables(3)/),purgeana=.true.)
431  !call display(v7ddballe%vol7d)
432 
433  call l4f_category_log(category,l4f_info, "filtered N time="//t2c(size(v7ddballe%vol7d%time)))
434 
435  call l4f_category_log(category,l4f_info,"initialize QC")
436 
437  ! chiamiamo il "costruttore" per il Q.C.
438 
439  call init(v7dqctem,v7ddballe%vol7d,var(:nvar),timei=timeiqc,timef=timefqc, &
440  coordmin=v7dana%ana(iana)%coord, coordmax=v7dana%ana(iana)%coord,&
441 ! data_id_in=v7ddballe%data_id, &
442  dsne=dsne, usere=usere, passworde=passworde,&
443  dsntem=dsntem, usertem=usertem, passwordtem=passwordtem,&
444  height2level=height2level, operation=operation,&
445  categoryappend="temporal")
446 
447 ! print *,">>>>>> Clima Temporal Volume <<<<<<"
448 ! call display(v7dqcspa%clima)
449 
450  call alloc(v7dqctem)
451 
452  ! temporal QC
453  call l4f_category_log(category,l4f_info,"start temporal QC")
454  call quacontem(v7dqctem,timemask= ( v7dqctem%v7d%time >= timeiqc .and. v7dqctem%v7d%time <= timefqc ))
455  call l4f_category_log(category,l4f_info,"end temporal QC")
456 
457  ! prepare data_id to be recreated
458  !deallocate(v7ddballe%data_id)
459  !nullify(v7ddballe%data_id)
460 
461  if (v7dqctem%operation == "run") then
462  call l4f_category_log(category,l4f_info,"start export data")
463  print *,"output volume"
464  call display(v7ddballe%vol7d)
465 
466  ! data_id to use is the new one
467  !v7ddballe%data_id => v7dqctem%data_id_out
468  !CALL export(v7ddballe,attr_only=.true.)
469  CALL export(v7ddballe, attr_only=.true.)
470  call l4f_category_log(category,l4f_info,"end export data")
471  end if
472 
473  call delete(v7dqctem)
474  ! data_id was allready deleted
475  !nullify(v7ddballe%data_id)
476  call delete(v7ddballe)
477 
478 end do
479 
480 call delete(v7dana)
481 
482 !close logger
483 call l4f_category_delete(category)
484 ier=l4f_fini()
485 
486 end program v7d_qctem
Generated with Doxygen.