45 integer :: category,io,ier,i,iun,n,ind
46 character(len=512):: a_name,output_format, output_template
50 TYPE(geo_coord
) :: coordmin, coordmax
51 TYPE(datetime
) :: time,ti, tf, timei, timef, timeiqc, timefqc
55 type(ncar_plot
) :: plot
58 integer,
parameter :: maxvar=10
59 character(len=6) :: var(maxvar)=cmiss
61 character(len=512) :: dsn=
'test1',user=
'test',password=
''
62 character(len=512) :: dsne=
'test',usere=
'test',passworde=
''
63 character(len=512) :: dsnspa=
'test',userspa=
'test',passwordspa=
''
65 integer :: years=imiss,months=imiss,days=imiss,hours=imiss,yeare=imiss,monthe=imiss,daye=imiss,houre=imiss,nvar=0
66 doubleprecision :: lons=dmiss,lats=dmiss,lone=dmiss,late=dmiss
67 integer :: year, month, day, hour
68 logical :: height2level=.false.,doplot=.false.,version
69 CHARACTER(len=512) :: input_file, output_file
70 INTEGER :: optind, optstatus, ninput
71 CHARACTER(len=20) :: operation
75 REAL,
DIMENSION(11) :: perc_vals=(/(10.*i,i=0,10)/)
76 REAL,
DIMENSION(size(perc_vals)-1) :: ndi
77 REAL,
DIMENSION(size(perc_vals)) :: nlimbins
78 double precision,
dimension(2) :: percentile
80 type(vol7d_level) :: level,levelo
81 type(vol7d_timerange
) :: timerange,timerangeo
82 type(vol7d_var
) :: varia,variao
83 CHARACTER(len=vol7d_ana_lenident) :: ident
84 integer :: indcana,indctime,indclevel,indctimerange,indcdativarr,indcnetwork,indcattr
85 INTEGER,
POINTER :: w_s(:), w_e(:)
88 integer :: status,grunit
91 namelist /odbc/ dsn,user,password,dsne,usere,passworde
92 namelist /odbcspa/ dsnspa,userspa,passwordspa
94 namelist /switchspa/ height2level,doplot
95 namelist /minmax/ years,months,days,hours,lons,lats,yeare,monthe,daye,houre,lone,late
96 namelist /varlist/ var
102 call l4f_launcher(a_name,a_name_force=
"v7d_qcspa")
105 category=l4f_category_get(a_name//
".main")
109 opt = optionparser_new(description_msg= &
110 'Spatial quality control: compute gradient; compute NDI from gradient; apply quality control', &
111 usage_msg=
'Usage: v7d_qcspa [options] [inputfile1] [inputfile2...] [outputfile] \n&
112 &If input-format is of file type, inputfile ''-'' indicates stdin,'&
114 //
'if output-format is of database type, outputfile specifies &
115 &database access info in the form user/password@dsn,' &
117 //
'if empty or ''-'', a suitable default is used. &
122 'operation to execute: ''gradient'' compute gradient and write on files; '&
123 //
'''ndi'' compute NDI from gradient;' &
124 //
'''run'' apply quality control ')
129 CALL
optionparser_add(opt,
' ',
'output-format', output_format,
'native', help= &
130 'format of output file for "ndi" operation only, in the form ''name[:template]''; ''native'' for vol7d &
131 &native binary format (no template to be specified)'&
133 //
'; ''BUFR'' and ''CREX'' for corresponding formats, with template as an alias like ''synop'', ''metar'', &
134 &''temp'', ''generic'', empty for ''generic'''&
139 CALL optionparser_add_help(opt,
'h',
'help', help=
'show an help message and exit')
140 CALL
optionparser_add(opt,
' ',
'version', version, help=
'show version and exit')
143 CALL optionparser_parse(opt, optind, optstatus)
145 IF (optstatus == optionparser_help)
THEN
147 ELSE IF (optstatus == optionparser_err)
THEN
149 CALL raise_fatal_error()
152 WRITE(*,
'(A,1X,A)')
'v7d_qcspa',version
157 if (operation /=
"gradient" .and. operation /=
"ndi" .and. operation /=
"run")
then
158 CALL optionparser_printhelp(opt)
160 'argument to --operation is wrong')
161 CALL raise_fatal_error()
166 n = word_split(output_format, w_s, w_e,
':')
168 output_template = output_format(w_s(2):w_e(2))
169 output_format(w_e(1)+1:) =
' '
173 if (operation ==
"ndi")
then
178 CALL optionparser_printhelp(opt)
180 CALL raise_fatal_error()
182 CALL optionparser_printhelp(opt)
184 CALL raise_fatal_error()
186 CALL getarg(iargc(), output_file)
191 call
init(timerangeo)
194 do ninput = optind, iargc()-1
195 call getarg(ninput, input_file)
199 open (grunit,file=input_file,status=
"old")
200 read (grunit,*,iostat=status) level,timerange,varia
201 if (status /= 0)
then
207 if (
c_e(levelo))
then
210 if ( level /= levelo .or. timerange /= timerangeo .or. varia /= variao )
then
211 call
l4f_category_log(category,l4f_error,
"Error reading grad files: file are incoerent")
216 timerangeo = timerange
221 read (grunit,*,iostat=iostat) val
222 if (iostat /= 0)
exit
223 if (val /= 0.) call
insert(grad,val)
229 CALL
l4f_category_log(category,l4f_info,
"compute percentile to remove the tails")
234 call normalizeddensityindex(pack(grad%array(:grad%arraysize),&
235 mask=(grad%array(:grad%arraysize) < percentile(2) )), perc_vals, ndi, nlimbins)
238 call
delete(v7dqcspa%clima)
239 CALL
init(v7dqcspa%clima, time_definition=0)
240 call vol7d_alloc(v7dqcspa%clima,nana=
size(perc_vals)-1, &
241 nlevel=1, ntimerange=1, &
242 ndativarr=1, nnetwork=1,ntime=1,ndativarattrr=1,ndatiattrr=1)
244 call vol7d_alloc_vol(v7dqcspa%clima,inivol=.true.)
246 call
init(v7dqcspa%clima%network(1),name=
"qcspa-ndi")
247 v7dqcspa%clima%level=level
248 v7dqcspa%clima%timerange=timerange
249 v7dqcspa%clima%dativar%r(1)=varia
250 v7dqcspa%clima%dativarattr%r(1)=varia
251 call
init(v7dqcspa%clima%datiattr%r(1), btable=
"*B33209")
252 cyclicdt = cyclicdatetime_new(chardate=
"/////////")
253 v7dqcspa%clima%time(1)=cyclicdatetime_to_conventional(cyclicdt)
262 do indcana=1,
size(perc_vals)-1
263 write(ident,
'("#",i2.2,2i3.3)')0,0,nint(perc_vals(indcana))
264 call
init(v7dqcspa%clima%ana(indcana),ident=ident,lat=0d0,lon=0d0)
265 if (
c_e(nlimbins(indcana)).and.
c_e(ndi(indcana)))
then
266 ind=index_c(spa_btable,varia%btable)
267 v7dqcspa%clima%voldatir(indcana, indctime, indclevel, indctimerange, indcdativarr, indcnetwork)=&
268 nlimbins(indcana)*spa_a(ind) + spa_b(ind)
269 v7dqcspa%clima%voldatiattrr(indcana, indctime, indclevel, indctimerange, indcdativarr, indcnetwork,indcattr)=&
274 print *,
">>>>>> NDI Volume <<<<<<"
277 IF (output_format ==
'native')
THEN
278 IF (output_file ==
'-')
THEN
279 CALL
l4f_category_log(category, l4f_info,
'trying /dev/stdout as stdout unit.')
280 output_file=
'/dev/stdout'
283 OPEN(iun, file=output_file, form=
'UNFORMATTED', access=stream_if_possible)
284 CALL
export(v7dqcspa%clima, unit=iun)
286 CALL
delete(v7dqcspa%clima)
289 ELSE IF (output_format ==
'BUFR' .OR. output_format ==
'CREX' .OR. output_format ==
'dba')
THEN
290 IF (output_format ==
'BUFR' .OR. output_format ==
'CREX')
THEN
291 IF (output_file ==
'-')
THEN
292 CALL
l4f_category_log(category, l4f_info,
'trying /dev/stdout as stdout unit.')
293 output_file=
'/dev/stdout'
297 ELSE IF (output_format ==
'dba')
THEN
298 CALL parse_dba_access_info(output_file, dsn, user, password)
302 IF (output_template ==
'') output_template =
'generic'
304 CALL
init(v7d_dba_out, filename=output_file, format=output_format, &
305 dsn=dsn, user=user, password=password, file=file, write=.true., wipe=file)
307 v7d_dba_out%vol7d = v7dqcspa%clima
308 CALL
init(v7dqcspa%clima)
309 CALL
export(v7d_dba_out, template=output_template)
323 open(10,file=
'qc.nml',status=
'old')
324 read(10,nml=odbc,iostat=io)
325 if ( io == 0 )
read(10,nml=odbcspa,iostat=io)
326 if ( io == 0 )
read(10,nml=switchspa,iostat=io)
327 if ( io == 0 )
read(10,nml=minmax,iostat=io)
328 if ( io == 0 )
read(10,nml=varlist,iostat=io)
332 call raise_error(
"Error reading namelist qc.nml")
348 CALL
init(ti, year=years, month=months, day=days, hour=hours)
349 CALL
init(tf, year=yeare, month=monthe, day=daye, hour=houre)
350 print *,
"time extreme"
355 CALL
init(coordmin,lat=lats,lon=lons)
356 CALL
init(coordmax,lat=late,lon=lone)
359 print*,
"lon lat minimum -> ",
to_char(coordmin)
361 print*,
"lon lat maximum -> ",
to_char(coordmax)
390 call
init(plot,pstype=
'PS', orient=
'LANDSCAPE',color=
'COLOR',file=
"v7d_qcspa.ps")
393 DO WHILE (time <= tf)
394 timei = time - timedelta_new(minute=30)
395 timef = time + timedelta_new(minute=30)
396 timeiqc = time - timedelta_new(minute=15)
397 timefqc = time + timedelta_new(minute=15)
401 CALL
init(v7ddballe,dsn=dsn,user=user,password=password,write=.true.,wipe=.false.,categoryappend=
"data-"//
t2c(time))
404 CALL
import(v7ddballe,var=var(:nvar),varkind=(/(
"r",i=1,nvar)/),&
405 anavar=(/
"B07030"/),anavarkind=(/
"r"/),&
406 attr=(/qcattrvarsbtables(1),qcattrvarsbtables(2),qcattrvarsbtables(4)/),attrkind=(/
"b",
"b",
"b"/)&
407 ,timei=timei,timef=timef,coordmin=coordmin,coordmax=coordmax)
420 call vol7d_peeling(v7ddballe%vol7d,keep_attr=(/qcattrvarsbtables(4)/),purgeana=.true.)
423 call
l4f_category_log(category,l4f_info,
"filtered N staz="//
t2c(
size(v7ddballe%vol7d%ana)))
430 call
init(v7dqcspa,v7ddballe%vol7d,var(:nvar),timei=timeiqc,timef=timefqc, coordmin=coordmin, coordmax=coordmax,&
432 dsne=dsne, usere=usere, passworde=passworde,&
433 dsnspa=dsnspa, userspa=userspa, passwordspa=passwordspa,&
434 height2level=height2level, operation=operation,&
435 categoryappend=
"space"//
t2c(time))
449 call quaconspa(v7dqcspa,timetollerance=timedelta_new(minute=15),noborder=.true.,&
450 timemask= ( v7dqcspa%v7d%time >= timeiqc .and. v7dqcspa%v7d%time <= timefqc ))
456 call plot_triangles(plot,v7dqcspa%co,v7dqcspa%tri,logo=
"Time: "//
t2c(timeiqc)//
" to "//
t2c(timefqc))
466 if (v7dqcspa%operation ==
"run")
then
486 time = time + timedelta_new(minute=30)
496 call l4f_category_delete(category)
499 end program v7d_qcspa