40 integer :: category,io,ier,i,iun,n,ind
41 character(len=512):: a_name,output_format, output_template
45 TYPE(geo_coord
) :: coordmin, coordmax
46 TYPE(datetime
) :: time, timeiqc, timefqc
50 integer,
parameter :: maxvar=10
51 character(len=6) :: var(maxvar)=cmiss
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=
''
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
67 REAL,
DIMENSION(11) :: perc_vals=(/(10.*i,i=0,10)/)
68 REAL,
DIMENSION(size(perc_vals)-1) :: ndi
69 REAL,
DIMENSION(size(perc_vals)) :: nlimbins
70 double precision,
dimension(2) :: percentile
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(:)
81 namelist /odbc/ dsn,user,password,dsne,usere,passworde
82 namelist /odbctem/ dsntem,usertem,passwordtem
84 namelist /switchtem/ height2level
85 namelist /minmax/ years,months,days,hours,lons,lats,yeare,monthe,daye,houre,lone,late
86 namelist /varlist/ var
92 call l4f_launcher(a_name,a_name_force=
"v7d_qctem")
95 category=l4f_category_get(a_name//
".main")
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,'&
104 //
'if output-format is of database type, outputfile specifies &
105 &database access info in the form user/password@dsn,' &
107 //
'if empty or ''-'', a suitable default is used. &
112 'operation to execute: ''gradient'' compute gradient and write on files; '&
113 //
'''ndi'' compute NDI from gradient;' &
114 //
'''run'' apply quality control ')
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)'&
123 //
'; ''BUFR'' and ''CREX'' for corresponding formats, with template as an alias like ''synop'', ''metar'', &
124 &''temp'', ''generic'', empty for ''generic'''&
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')
133 CALL optionparser_parse(opt, optind, optstatus)
135 IF (optstatus == optionparser_help)
THEN
137 ELSE IF (optstatus == optionparser_err)
THEN
139 CALL raise_fatal_error()
142 WRITE(*,
'(A,1X,A)')
'v7d_qctem',version
147 if (operation /=
"gradient" .and. operation /=
"ndi" .and. operation /=
"run")
then
148 CALL optionparser_printhelp(opt)
150 'argument to --operation is wrong')
151 CALL raise_fatal_error()
156 n = word_split(output_format, w_s, w_e,
':')
158 output_template = output_format(w_s(2):w_e(2))
159 output_format(w_e(1)+1:) =
' '
163 if (operation ==
"ndi")
then
168 CALL optionparser_printhelp(opt)
170 CALL raise_fatal_error()
172 CALL optionparser_printhelp(opt)
174 CALL raise_fatal_error()
176 CALL getarg(iargc(), output_file)
179 call
init(timerangeo)
182 do ninput = optind, iargc()-1
183 call getarg(ninput, input_file)
186 open (10,file=input_file,status=
"old")
187 read (10,*) level,timerange,varia
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")
196 timerangeo = timerange
201 read (10,*,iostat=iostat) val
202 if (iostat /= 0)
exit
203 if (val /= 0.) call
insert(grad,val)
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)
215 call vol7d_alloc_vol(v7dqctem%clima,inivol=.true.)
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")
222 cyclicdt = cyclicdatetime_new(chardate=
"/////////")
223 v7dqctem%clima%time(1)=cyclicdatetime_to_conventional(cyclicdt)
233 call
init(v7dqctem%clima%network(indcnetwork),name=
"qctemsndi")
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./))
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)
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)=&
257 call
init(v7dqctem%clima%network(indcnetwork),name=
"qctemgndi")
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./))
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)
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)=&
284 IF (output_format ==
'native')
THEN
285 IF (output_file ==
'-')
THEN
286 CALL
l4f_category_log(category, l4f_info,
'trying /dev/stdout as stdout unit.')
287 output_file=
'/dev/stdout'
290 OPEN(iun, file=output_file, form=
'UNFORMATTED', access=stream_if_possible)
291 CALL
export(v7dqctem%clima, unit=iun)
293 CALL
delete(v7dqctem%clima)
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'
304 ELSE IF (output_format ==
'dba')
THEN
305 CALL parse_dba_access_info(output_file, dsn, user, password)
309 IF (output_template ==
'') output_template =
'generic'
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)
314 v7d_dba_out%vol7d = v7dqctem%clima
315 CALL
init(v7dqctem%clima)
316 CALL
export(v7d_dba_out, template=output_template)
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)
355 CALL
init(timeiqc, year=years, month=months, day=days, hour=hours)
356 CALL
init(timefqc, year=yeare, month=monthe, day=daye, hour=houre)
362 CALL
init(coordmin,lat=lats,lon=lons)
363 CALL
init(coordmax,lat=late,lon=lone)
366 print*,
"lon lat minimum -> ",
to_char(coordmin)
368 print*,
"lon lat maximum -> ",
to_char(coordmax)
384 CALL
init(v7ddballeana,dsn=dsn,user=user,password=password,write=.false.,wipe=.false.,categoryappend=
"QCtarget-anaonly")
386 CALL
import(v7ddballeana,var=var(:nvar),timei=timeiqc,timef=timefqc,coordmin=coordmin,coordmax=coordmax,anaonly=.true.)
387 call vol7d_copy(v7ddballeana%vol7d,v7dana)
397 print *,
"anagrafica:"
400 do iana=1,
size(v7dana%ana)
405 CALL
init(v7ddballe,dsn=dsn,user=user,password=password,write=.true.,wipe=.false.,categoryappend=
"QCtarget-"//
t2c(time))
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))
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 !')
420 print *,
"input volume"
431 call vol7d_peeling(v7ddballe%vol7d,keep_attr=(/qcattrvarsbtables(3)/),purgeana=.true.)
434 call
l4f_category_log(category,l4f_info,
"filtered N time="//
t2c(
size(v7ddballe%vol7d%time)))
440 call
init(v7dqctem,v7ddballe%vol7d,var(:nvar),timei=timeiqc,timef=timefqc, &
441 coordmin=v7dana%ana(iana)%coord, coordmax=v7dana%ana(iana)%coord,&
443 dsne=dsne, usere=usere, passworde=passworde,&
444 dsntem=dsntem, usertem=usertem, passwordtem=passwordtem,&
445 height2level=height2level, operation=operation,&
446 categoryappend=
"temporal")
455 call quacontem(v7dqctem,timemask= ( v7dqctem%v7d%time >= timeiqc .and. v7dqctem%v7d%time <= timefqc ))
462 if (v7dqctem%operation ==
"run")
then
464 print *,
"output volume"
470 CALL
export(v7ddballe, attr_only=.true.)
484 call l4f_category_delete(category)
487 end program v7d_qctem