Sample program for module qctem.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20#include "config.h"
21
22program v7d_qctem
23
32
36
37implicit none
38
39integer :: category,io,ier,i,iun,n,ind
40character(len=512):: a_name,output_format, output_template
41
42
43TYPE(optionparser) :: opt
44TYPE(geo_coord) :: coordmin, coordmax
45TYPE(datetime) :: time, timeiqc, timefqc
46type(qctemtype) :: v7dqctem
47type(vol7d) :: v7dana
48
49integer, parameter :: maxvar=10
50character(len=6) :: var(maxvar)=cmiss
51#ifdef HAVE_DBALLE
52type(vol7d_dballe) :: v7ddballe,v7ddballeana,v7d_dba_out
53character(len=512) :: dsn='test1',user='test',password=''
54character(len=512) :: dsne='test',usere='test',passworde=''
55character(len=512) :: dsntem='test',usertem='test',passwordtem=''
56#endif
57integer :: years=imiss,months=imiss,days=imiss,hours=imiss,yeare=imiss,monthe=imiss,daye=imiss,houre=imiss,nvar=0
58doubleprecision :: lons=dmiss,lats=dmiss,lone=dmiss,late=dmiss
59logical :: height2level=.false.,version
60CHARACTER(len=512) :: input_file, output_file
61INTEGER :: optind, optstatus, ninput
62CHARACTER(len=20) :: operation
63TYPE(ARRAYOF_REAL):: grad
64real :: val
65integer :: iostat
66REAL, DIMENSION(11) :: perc_vals=(/(10.*i,i=0,10)/)
67REAL, DIMENSION(size(perc_vals)-1) :: ndi
68REAL, DIMENSION(size(perc_vals)) :: nlimbins
69double precision, dimension(2) :: percentile
70TYPE(cyclicdatetime) :: cyclicdt
71type(vol7d_level) :: level,levelo
72type(vol7d_timerange) :: timerange,timerangeo
73type(vol7d_var) :: varia,variao
74CHARACTER(len=vol7d_ana_lenident) :: ident
75integer :: indcana,indctime,indclevel,indctimerange,indcdativarr,indcnetwork,indcattr,iana
76INTEGER,POINTER :: w_s(:), w_e(:)
77logical :: file
78
79#ifdef HAVE_DBALLE
80namelist /odbc/ dsn,user,password,dsne,usere,passworde
81namelist /odbctem/ dsntem,usertem,passwordtem
82#endif
83namelist /switchtem/ height2level
84namelist /minmax/ years,months,days,hours,lons,lats,yeare,monthe,daye,houre,lone,late
85namelist /varlist/ var
86
87
89
90
91call l4f_launcher(a_name,a_name_force="v7d_qctem")
92
93
94category=l4f_category_get(a_name//".main")
95
96
97
98opt = 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,'&
102#ifdef HAVE_DBALLE
103 //'if output-format is of database type, outputfile specifies &
104 &database access info in YRI form,' &
105#endif
106 //'if empty or ''-'', a suitable default is used. &
107&')
108
109
111 'operation to execute: ''gradient'' compute gradient and write on files; '&
112 //'''ndi'' compute NDI from gradient;' &
113 //'''run'' apply quality control ')
114
115
116
117output_template = ''
118CALL 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)'&
121#ifdef HAVE_DBALLE
122 //'; ''BUFR'' and ''CREX'' for corresponding formats, with template as an alias like ''synop'', ''metar'', &
123 &''temp'', ''generic'', empty for ''generic'''&
124#endif
125)
126
127
128CALL optionparser_add_help(opt, 'h', 'help', help='show an help message and exit')
129CALL optionparser_add(opt,
' ',
'version', version, help=
'show version and exit')
130
131
132CALL optionparser_parse(opt, optind, optstatus)
133
134IF (optstatus == optionparser_help) THEN
135 CALL exit(0)
136ELSE IF (optstatus == optionparser_err) THEN
138 CALL raise_fatal_error()
139ENDIF
140IF (version) THEN
141 WRITE(*,'(A,1X,A)')'v7d_qctem',version
142 CALL exit(0)
143ENDIF
144
145
146if (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()
151end if
152
153
154
155n = word_split(output_format, w_s, w_e, ':')
156IF (n >= 2) THEN
157 output_template = output_format(w_s(2):w_e(2))
158 output_format(w_e(1)+1:) = ' '
159ENDIF
160DEALLOCATE(w_s, w_e)
161
162if (operation == "ndi") then
163
164
165 i = iargc() - optind
166 IF (i < 0) THEN
167 CALL optionparser_printhelp(opt)
169 CALL raise_fatal_error()
170 ELSE IF (i < 1) THEN
171 CALL optionparser_printhelp(opt)
173 CALL raise_fatal_error()
174 ENDIF
175 CALL getarg(iargc(), output_file)
176
178 call init(timerangeo)
180
181 do ninput = optind, iargc()-1
182 call getarg(ninput, input_file)
183
185 open (10,file=input_file,status="old")
186 read (10,*) level,timerange,varia
187
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")
191 call raise_error("")
192 end if
193 else
194 levelo = level
195 timerangeo = timerange
196 variao = varia
197 end if
198
199 do while (.true.)
200 read (10,*,iostat=iostat) val
201 if (iostat /= 0) exit
202 if (val /= 0.)
call insert(grad,val)
203 end do
204
205 close(10)
206 end do
207
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)
213
214 call vol7d_alloc_vol(v7dqctem%clima,inivol=.true.)
215
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)
223
224 indctime=1
225 indclevel=1
226 indcattr=1
227 indcdativarr=1
228 indctimerange=1
229
230
231 indcnetwork=1
232 call init(v7dqctem%clima%network(indcnetwork),name=
"qctemsndi")
233
235
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./))
238
239
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)
242
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)=&
251 ndi(indcana)*100.
252 end if
253 end do
254
255 indcnetwork=2
256 call init(v7dqctem%clima%network(indcnetwork),name=
"qctemgndi")
257
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./))
261
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)
264
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)=&
273 ndi(indcana)*100.
274 end if
275 end do
276
279
280
281
282
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'
287 ENDIF
288 iun = getunit()
289 OPEN(iun, file=output_file, form='UNFORMATTED', access='STREAM')
290 CALL export(v7dqctem%clima, unit=iun)
291 CLOSE(iun)
292 CALL delete(v7dqctem%clima)
293
294#ifdef HAVE_DBALLE
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'
300 ENDIF
301 file=.true.
302
303 ELSE IF (output_format == 'dba') THEN
304 file=.false.
305 ENDIF
306
307 IF (output_template == '') output_template = 'generic'
308
309 CALL init(v7d_dba_out, filename=output_file, format=output_format, &
310 dsn=output_file, file=file, write=.true., wipe=file)
311
312 v7d_dba_out%vol7d = v7dqctem%clima
313 CALL init(v7dqctem%clima)
314 CALL export(v7d_dba_out, template=output_template)
316#endif
317
318 end if
319
320 stop
321
322end if
323
324
325
326
327
328open(10,file='qc.nml',status='old')
329read(10,nml=odbc,iostat=io)
330if ( io == 0 ) read(10,nml=odbctem,iostat=io)
331if ( io == 0 ) read(10,nml=switchtem,iostat=io)
332if ( io == 0 ) read(10,nml=minmax,iostat=io)
333if ( io == 0 ) read(10,nml=varlist,iostat=io)
334
335if (io /= 0 )then
337 call raise_error()
338end if
339close(10)
340
341
342
343
344
345
347
348if (nvar == 0) then
350 call raise_error()
351end if
352
353CALL init(timeiqc, year=years, month=months, day=days, hour=hours)
354CALL init(timefqc, year=yeare, month=monthe, day=daye, hour=houre)
355
358
359
360CALL init(coordmin,lat=lats,lon=lons)
361CALL init(coordmax,lat=late,lon=lone)
362
363
364print*,
"lon lat minimum -> ",
to_char(coordmin)
365
366print*,
"lon lat maximum -> ",
to_char(coordmax)
367
368
370do i=1,nvar
372enddo
379
380
381
382CALL init(v7ddballeana,dsn=dsn,write=.false.,wipe=.false.,categoryappend=
"QCtarget-anaonly")
384CALL import(v7ddballeana,var=var(:nvar),timei=timeiqc,timef=timefqc,coordmin=coordmin,coordmax=coordmax,anaonly=.true.)
385call vol7d_copy(v7ddballeana%vol7d,v7dana)
387
388
389
390
391
392
393
394
395print *,"anagrafica:"
397
398do iana=1, size(v7dana%ana)
399
401
402
403 CALL init(v7ddballe,dsn=dsn,write=.true.,wipe=.false.,categoryappend=
"QCtarget-"//
t2c(time))
405
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))
410
411
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 !')
415 cycle
416 end if
417
418 print *,"input volume"
422
424
425
426
427 qcpar%att=bmiss
428
429 call vol7d_peeling(v7ddballe%vol7d,keep_attr=(/qcattrvarsbtables(3)/),purgeana=.true.)
430
431
432 call l4f_category_log(category,l4f_info,
"filtered N time="//
t2c(
size(v7ddballe%vol7d%time)))
433
435
436
437
438 call init(v7dqctem,v7ddballe%vol7d,var(:nvar),timei=timeiqc,timef=timefqc, &
439 coordmin=v7dana%ana(iana)%coord, coordmax=v7dana%ana(iana)%coord,&
440
441 dsne=dsne, usere=usere, passworde=passworde,&
442 dsntem=dsntem, usertem=usertem, passwordtem=passwordtem,&
443 height2level=height2level, operation=operation,&
444 categoryappend="temporal")
445
446
447
448
450
451
453 call quacontem(v7dqctem,timemask= ( v7dqctem%v7d%time >= timeiqc .and. v7dqctem%v7d%time <= timefqc ))
455
456
457
458
459
460 if (v7dqctem%operation == "run") then
462 print *,"output volume"
464
465
466
467
468 CALL export(v7ddballe, attr_only=.true.)
470 end if
471
473
474
476
477end do
478
480
481
482call l4f_category_delete(category)
484
485end program v7d_qctem
Destructor for finalizing an array object.
Method for inserting elements of the array at a desired position.
Operatore di valore assoluto di un intervallo.
Costruttori per le classi datetime e timedelta.
Functions that return a trimmed CHARACTER representation of the input variable.
Restituiscono il valore dell'oggetto in forma di stringa stampabile.
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...
Emit log message for a category with specific priority.
Global log4fortran constructor.
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