44 integer :: category,io,ier,i,iun,n,ind
45 character(len=512):: a_name,output_format, output_template
48 TYPE(optionparser) :: opt
49 TYPE(geo_coord) :: coordmin, coordmax
50 TYPE(datetime) :: time,ti, tf, timei, timef, timeiqc, timefqc
51 type(qcspatype) :: v7dqcspa
52 type(vol7d_dballe) :: v7ddballe
54 type(ncar_plot) :: plot
57 integer,
parameter :: maxvar=10
58 character(len=6) :: var(maxvar)=cmiss
60 character(len=512) :: dsn=
'test1',user=
'test',password=
'' 61 character(len=512) :: dsne=
'test',usere=
'test',passworde=
'' 62 character(len=512) :: dsnspa=
'test',userspa=
'test',passwordspa=
'' 64 integer :: years=imiss,months=imiss,days=imiss,hours=imiss,yeare=imiss,monthe=imiss,daye=imiss,houre=imiss,nvar=0
65 doubleprecision :: lons=dmiss,lats=dmiss,lone=dmiss,late=dmiss
66 integer :: year, month, day, hour
67 logical :: height2level=.false.,doplot=.false.,version
68 CHARACTER(len=512) :: input_file, output_file
69 INTEGER :: optind, optstatus, ninput
70 CHARACTER(len=20) :: operation
71 TYPE(ARRAYOF_REAL):: grad
74 REAL,
DIMENSION(11) :: perc_vals=(/(10.*i,i=0,10)/)
75 REAL,
DIMENSION(size(perc_vals)-1) :: ndi
76 REAL,
DIMENSION(size(perc_vals)) :: nlimbins
77 double precision,
dimension(2) :: percentile
78 TYPE(cyclicdatetime) :: cyclicdt
79 type(vol7d_level) :: level,levelo
80 type(vol7d_timerange) :: timerange,timerangeo
81 type(vol7d_var) :: varia,variao
82 CHARACTER(len=vol7d_ana_lenident) :: ident
83 integer :: indcana,indctime,indclevel,indctimerange,indcdativarr,indcnetwork,indcattr
84 INTEGER,
POINTER :: w_s(:), w_e(:)
85 TYPE(vol7d_dballe) :: v7d_dba_out
87 integer :: status,grunit
90 namelist /odbc/ dsn,user,password,dsne,usere,passworde
91 namelist /odbcspa/ dsnspa,userspa,passwordspa
93 namelist /switchspa/ height2level,doplot
94 namelist /minmax/ years,months,days,hours,lons,lats,yeare,monthe,daye,houre,lone,late
95 namelist /varlist/ var
101 call l4f_launcher(a_name,a_name_force=
"v7d_qcspa")
104 category=l4f_category_get(a_name//
".main")
108 opt = optionparser_new(description_msg= &
109 'Spatial quality control: compute gradient; compute NDI from gradient; apply quality control', &
110 usage_msg=
'Usage: v7d_qcspa [options] [inputfile1] [inputfile2...] [outputfile] \n& 111 &If input-format is of file type, inputfile ''-'' indicates stdin,'&
113 //
'if output-format is of database type, outputfile specifies & 114 &database access info in URI form,' &
116 //
'if empty or ''-'', a suitable default is used. & 121 'operation to execute: ''gradient'' compute gradient and write on files; '&
122 //
'''ndi'' compute NDI from gradient;' &
123 //
'''run'' apply quality control ')
128 CALL optionparser_add(opt,
' ',
'output-format', output_format,
'native', help= &
129 'format of output file for "ndi" operation only, in the form ''name[:template]''; ''native'' for vol7d & 130 &native binary format (no template to be specified)'&
132 //
'; ''BUFR'' and ''CREX'' for corresponding formats, with template as an alias like ''synop'', ''metar'', & 133 &''temp'', ''generic'', empty for ''generic'''&
138 CALL optionparser_add_help(opt,
'h',
'help', help=
'show an help message and exit')
139 CALL optionparser_add(opt,
' ',
'version', version, help=
'show version and exit')
142 CALL optionparser_parse(opt, optind, optstatus)
144 IF (optstatus == optionparser_help)
THEN 146 ELSE IF (optstatus == optionparser_err)
THEN 148 CALL raise_fatal_error()
151 WRITE(*,
'(A,1X,A)')
'v7d_qcspa',version
156 if (operation /=
"gradient" .and. operation /=
"ndi" .and. operation /=
"run")
then 157 CALL optionparser_printhelp(opt)
159 'argument to --operation is wrong')
160 CALL raise_fatal_error()
165 n = word_split(output_format, w_s, w_e,
':')
167 output_template = output_format(w_s(2):w_e(2))
168 output_format(w_e(1)+1:) =
' ' 172 if (operation ==
"ndi")
then 177 CALL optionparser_printhelp(opt)
179 CALL raise_fatal_error()
181 CALL optionparser_printhelp(opt)
183 CALL raise_fatal_error()
185 CALL getarg(iargc(), output_file)
190 call init(timerangeo)
193 do ninput = optind, iargc()-1
194 call getarg(ninput, input_file)
198 open (grunit,file=input_file,status=
"old")
199 read (grunit,*,iostat=status) level,timerange,varia
200 if (status /= 0)
then 206 if (
c_e(levelo))
then 209 if ( level /= levelo .or. timerange /= timerangeo .or. varia /= variao )
then 210 call l4f_category_log(category,l4f_error,
"Error reading grad files: file are incoerent")
215 timerangeo = timerange
220 read (grunit,*,iostat=iostat) val
221 if (iostat /= 0)
exit 222 if (val /= 0.)
call insert(grad,val)
228 CALL l4f_category_log(category,l4f_info,
"compute percentile to remove the tails")
233 call normalizeddensityindex(pack(grad%array(:grad%arraysize),&
234 mask=(grad%array(:grad%arraysize) < percentile(2) )), perc_vals, ndi, nlimbins)
237 call delete(v7dqcspa%clima)
238 CALL init(v7dqcspa%clima, time_definition=0)
239 call vol7d_alloc(v7dqcspa%clima,nana=
size(perc_vals)-1, &
240 nlevel=1, ntimerange=1, &
241 ndativarr=1, nnetwork=1,ntime=1,ndativarattrr=1,ndatiattrr=1)
243 call vol7d_alloc_vol(v7dqcspa%clima,inivol=.true.)
245 call init(v7dqcspa%clima%network(1),name=
"qcspa-ndi")
246 v7dqcspa%clima%level=level
247 v7dqcspa%clima%timerange=timerange
248 v7dqcspa%clima%dativar%r(1)=varia
249 v7dqcspa%clima%dativarattr%r(1)=varia
250 call init(v7dqcspa%clima%datiattr%r(1), btable=
"*B33209")
251 cyclicdt = cyclicdatetime_new(chardate=
"/////////")
252 v7dqcspa%clima%time(1)=cyclicdatetime_to_conventional(cyclicdt)
261 do indcana=1,
size(perc_vals)-1
262 write(ident,
'("#",i2.2,2i3.3)')0,0,nint(perc_vals(indcana))
263 call init(v7dqcspa%clima%ana(indcana),ident=ident,lat=0d0,lon=0d0)
264 if (
c_e(nlimbins(indcana)).and.
c_e(ndi(indcana)))
then 265 ind=index_c(spa_btable,varia%btable)
266 v7dqcspa%clima%voldatir(indcana, indctime, indclevel, indctimerange, indcdativarr, indcnetwork)=&
267 nlimbins(indcana)*spa_a(ind) + spa_b(ind)
268 v7dqcspa%clima%voldatiattrr(indcana, indctime, indclevel, indctimerange, indcdativarr, indcnetwork,indcattr)=&
273 print *,
">>>>>> NDI Volume <<<<<<" 276 IF (output_format ==
'native')
THEN 277 IF (output_file ==
'-')
THEN 278 CALL l4f_category_log(category, l4f_info,
'trying /dev/stdout as stdout unit.')
279 output_file=
'/dev/stdout' 282 OPEN(iun, file=output_file, form=
'UNFORMATTED', access=
'STREAM')
283 CALL export(v7dqcspa%clima, unit=iun)
285 CALL delete(v7dqcspa%clima)
288 ELSE IF (output_format ==
'BUFR' .OR. output_format ==
'CREX' .OR. output_format ==
'dba')
THEN 289 IF (output_format ==
'BUFR' .OR. output_format ==
'CREX')
THEN 290 IF (output_file ==
'-')
THEN 291 CALL l4f_category_log(category, l4f_info,
'trying /dev/stdout as stdout unit.')
292 output_file=
'/dev/stdout' 296 ELSE IF (output_format ==
'dba')
THEN 300 IF (output_template ==
'') output_template =
'generic' 302 CALL init(v7d_dba_out, filename=output_file, format=output_format, &
303 dsn=output_file, file=file, write=.true., wipe=file)
305 v7d_dba_out%vol7d = v7dqcspa%clima
306 CALL init(v7dqcspa%clima)
307 CALL export(v7d_dba_out, template=output_template)
308 CALL delete(v7d_dba_out)
321 open(10,file=
'qc.nml',status=
'old')
322 read(10,nml=odbc,iostat=io)
323 if ( io == 0 )
read(10,nml=odbcspa,iostat=io)
324 if ( io == 0 )
read(10,nml=switchspa,iostat=io)
325 if ( io == 0 )
read(10,nml=minmax,iostat=io)
326 if ( io == 0 )
read(10,nml=varlist,iostat=io)
330 call raise_error(
"Error reading namelist qc.nml")
346 CALL init(ti, year=years, month=months, day=days, hour=hours)
347 CALL init(tf, year=yeare, month=monthe, day=daye, hour=houre)
348 print *,
"time extreme" 353 CALL init(coordmin,lat=lats,lon=lons)
354 CALL init(coordmax,lat=late,lon=lone)
357 print*,
"lon lat minimum -> ",
to_char(coordmin)
359 print*,
"lon lat maximum -> ",
to_char(coordmax)
388 call init(plot,pstype=
'PS', orient=
'LANDSCAPE',color=
'COLOR',file=
"v7d_qcspa.ps")
391 DO WHILE (time <= tf)
392 timei = time - timedelta_new(minute=30)
393 timef = time + timedelta_new(minute=30)
394 timeiqc = time - timedelta_new(minute=15)
395 timefqc = time + timedelta_new(minute=15)
399 CALL init(v7ddballe,dsn=dsn,write=.true.,wipe=.false.,categoryappend=
"data-"//
t2c(time))
402 CALL import(v7ddballe,var=var(:nvar),varkind=(/(
"r",i=1,nvar)/),&
403 anavar=(/
"B07030"/),anavarkind=(/
"r"/),&
404 attr=(/qcattrvarsbtables(1),qcattrvarsbtables(2),qcattrvarsbtables(4)/),attrkind=(/
"b",
"b",
"b"/)&
405 ,timei=timei,timef=timef,coordmin=coordmin,coordmax=coordmax)
418 call vol7d_peeling(v7ddballe%vol7d,keep_attr=(/qcattrvarsbtables(4)/),purgeana=.true.)
421 call l4f_category_log(category,l4f_info,
"filtered N staz="//
t2c(
size(v7ddballe%vol7d%ana)))
428 call init(v7dqcspa,v7ddballe%vol7d,var(:nvar),timei=timeiqc,timef=timefqc, coordmin=coordmin, coordmax=coordmax,&
430 dsne=dsne, usere=usere, passworde=passworde,&
431 dsnspa=dsnspa, userspa=userspa, passwordspa=passwordspa,&
432 height2level=height2level, operation=operation,&
433 categoryappend=
"space"//
t2c(time))
447 call quaconspa(v7dqcspa,timetollerance=timedelta_new(minute=15),noborder=.true.,&
448 timemask= ( v7dqcspa%v7d%time >= timeiqc .and. v7dqcspa%v7d%time <= timefqc ))
454 call plot_triangles(plot,v7dqcspa%co,v7dqcspa%tri,logo=
"Time: "//
t2c(timeiqc)//
" to "//
t2c(timefqc))
464 if (v7dqcspa%operation ==
"run")
then 479 call delete(v7dqcspa)
482 call delete(v7ddballe)
484 time = time + timedelta_new(minute=30)
494 call l4f_category_delete(category)
497 end program v7d_qcspa