39 integer :: category,io,ier,i,iun,n,ind
40 character(len=512):: a_name,output_format, output_template
43 TYPE(optionparser) :: opt
44 TYPE(geo_coord) :: coordmin, coordmax
45 TYPE(datetime) :: time, timeiqc, timefqc
46 type(qctemtype) :: v7dqctem
49 integer,
parameter :: maxvar=10
50 character(len=6) :: var(maxvar)=cmiss
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=
''
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
66 REAL,
DIMENSION(11) :: perc_vals=(/(10.*i,i=0,10)/)
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(:)
80 namelist /odbc/ dsn,user,password,dsne,usere,passworde
81 namelist /odbctem/ dsntem,usertem,passwordtem
83 namelist /switchtem/ height2level
84 namelist /minmax/ years,months,days,hours,lons,lats,yeare,monthe,daye,houre,lone,late
85 namelist /varlist/ var
91 call l4f_launcher(a_name,a_name_force=
"v7d_qctem")
94 category=l4f_category_get(a_name//
".main")
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,'&
103 //
'if output-format is of database type, outputfile specifies &
104 &database access info in YRI form,' &
106 //
'if empty or ''-'', a suitable default is used. &
111 'operation to execute: ''gradient'' compute gradient and write on files; '&
112 //
'''ndi'' compute NDI from gradient;' &
113 //
'''run'' apply quality control ')
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)'&
122 //
'; ''BUFR'' and ''CREX'' for corresponding formats, with template as an alias like ''synop'', ''metar'', &
123 &''temp'', ''generic'', empty for ''generic'''&
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')
132 CALL optionparser_parse(opt, optind, optstatus)
134 IF (optstatus == optionparser_help)
THEN
136 ELSE IF (optstatus == optionparser_err)
THEN
138 CALL raise_fatal_error()
141 WRITE(*,
'(A,1X,A)')
'v7d_qctem',version
146 if (operation /=
"gradient" .and. operation /=
"ndi" .and. operation /=
"run")
then
147 CALL optionparser_printhelp(opt)
149 'argument to --operation is wrong')
150 CALL raise_fatal_error()
155 n = word_split(output_format, w_s, w_e,
':')
157 output_template = output_format(w_s(2):w_e(2))
158 output_format(w_e(1)+1:) =
' '
162 if (operation ==
"ndi")
then
167 CALL optionparser_printhelp(opt)
169 CALL raise_fatal_error()
171 CALL optionparser_printhelp(opt)
173 CALL raise_fatal_error()
175 CALL getarg(iargc(), output_file)
178 call init(timerangeo)
181 do ninput = optind, iargc()-1
182 call getarg(ninput, input_file)
185 open (10,file=input_file,status=
"old")
186 read (10,*) level,timerange,varia
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")
195 timerangeo = timerange
200 read (10,*,iostat=iostat) val
201 if (iostat /= 0)
exit
202 if (val /= 0.)
call insert(grad,val)
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)
214 call vol7d_alloc_vol(v7dqctem%clima,inivol=.true.)
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")
221 cyclicdt = cyclicdatetime_new(chardate=
"/////////")
222 v7dqctem%clima%time(1)=cyclicdatetime_to_conventional(cyclicdt)
232 call init(v7dqctem%clima%network(indcnetwork),name=
"qctemsndi")
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./))
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)
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)=&
256 call init(v7dqctem%clima%network(indcnetwork),name=
"qctemgndi")
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./))
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)
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)=&
283 IF (output_format ==
'native')
THEN
284 IF (output_file ==
'-')
THEN
285 CALL l4f_category_log(category, l4f_info,
'trying /dev/stdout as stdout unit.')
286 output_file=
'/dev/stdout'
289 OPEN(iun, file=output_file, form=
'UNFORMATTED', access=
'STREAM')
290 CALL export(v7dqctem%clima, unit=iun)
292 CALL delete(v7dqctem%clima)
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'
303 ELSE IF (output_format ==
'dba')
THEN
307 IF (output_template ==
'') output_template =
'generic'
309 CALL init(v7d_dba_out, filename=output_file, format=output_format, &
310 dsn=output_file, file=file, write=.true., wipe=file)
312 v7d_dba_out%vol7d = v7dqctem%clima
313 CALL init(v7dqctem%clima)
314 CALL export(v7d_dba_out, template=output_template)
315 CALL delete(v7d_dba_out)
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)
353 CALL init(timeiqc, year=years, month=months, day=days, hour=hours)
354 CALL init(timefqc, year=yeare, month=monthe, day=daye, hour=houre)
360 CALL init(coordmin,lat=lats,lon=lons)
361 CALL init(coordmax,lat=late,lon=lone)
364 print*,
"lon lat minimum -> ",
to_char(coordmin)
366 print*,
"lon lat maximum -> ",
to_char(coordmax)
382 CALL init(v7ddballeana,dsn=dsn,write=.false.,wipe=.false.,categoryappend=
"QCtarget-anaonly")
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)
395 print *,
"anagrafica:"
398 do iana=1,
size(v7dana%ana)
403 CALL init(v7ddballe,dsn=dsn,write=.true.,wipe=.false.,categoryappend=
"QCtarget-"//
t2c(time))
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))
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)
418 print *,
"input volume"
429 call vol7d_peeling(v7ddballe%vol7d,keep_attr=(/qcattrvarsbtables(3)/),purgeana=.true.)
432 call l4f_category_log(category,l4f_info,
"filtered N time="//
t2c(
size(v7ddballe%vol7d%time)))
438 call init(v7dqctem,v7ddballe%vol7d,var(:nvar),timei=timeiqc,timef=timefqc, &
439 coordmin=v7dana%ana(iana)%coord, coordmax=v7dana%ana(iana)%coord,&
441 dsne=dsne, usere=usere, passworde=passworde,&
442 dsntem=dsntem, usertem=usertem, passwordtem=passwordtem,&
443 height2level=height2level, operation=operation,&
444 categoryappend=
"temporal")
453 call quacontem(v7dqctem,timemask= ( v7dqctem%v7d%time >= timeiqc .and. v7dqctem%v7d%time <= timefqc ))
460 if (v7dqctem%operation ==
"run")
then
462 print *,
"output volume"
468 CALL export(v7ddballe, attr_only=.true.)
472 call delete(v7dqctem)
475 call delete(v7ddballe)
482 call l4f_category_delete(category)
485 end program v7d_qctem
Operatore di valore assoluto di un intervallo.
Functions that return a trimmed CHARACTER representation of the input variable.
Export one or more geo_coordvect objects to a plain text file or to a file in ESRI/Shapefile format.
Import one or more geo_coordvect objects from a plain text file or for a file in ESRI/Shapefile forma...
Constructors for the two classes.
Represent geo_coord object in a pretty string.
Emit log message for a category with specific priority.
Global log4fortran constructor.
Function to check whether a value is missing or not.
Add a new option of a specific type.
Compute a set of percentiles for a random variable.
This module defines usefull general purpose function and subroutine.
Classi per la gestione delle coordinate temporali.
Classes for handling georeferenced sparse points in geographical corodinates.
classe per la gestione del logging
Definitions of constants and functions for working with missing values.
Utilities and defines for quality control.
Controllo di qualità temporale.
Module for parsing command-line optons.
Module for basic statistical computations taking into account missing data.
Classe per la gestione di un volume completo di dati osservati.
classe per import ed export di volumi da e in DB-All.e