39 integer :: category,io,ier,i,iun,n,ind
40 character(len=512):: a_name,output_format, output_template
44 TYPE(geo_coord
) :: coordmin, coordmax
45 TYPE(datetime
) :: time, timeiqc, timefqc
49 integer,
parameter :: maxvar=10
50 character(len=6) :: var(maxvar)=cmiss
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
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
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)
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)
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 !')
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.)
482 call l4f_category_delete(category)
485 end program v7d_qctem