106 character (len=255),
parameter:: subcategoryspa=
"QCspa"
108 integer,
parameter :: spa_nvar=1
109 CHARACTER(len=10) :: spa_btable(spa_nvar)=(/
"B12101"/)
111 real,
parameter :: spa_a(spa_nvar) = (/1.e5/)
113 real,
parameter :: spa_b(spa_nvar) = (/273.15/)
117 type (vol7d),
pointer :: v7d => null()
118 integer,
pointer :: data_id_in(:,:,:,:,:) => null()
119 integer,
pointer :: data_id_out(:,:,:,:,:) => null()
122 type(xy),
pointer :: co(:) => null()
123 type (triangles) :: tri
124 type (qcclitype) :: qccli
125 type (vol7d) :: clima
126 character(len=20):: operation
134 module procedure qcspainit
139 module procedure qcspaalloc
144 module procedure qcspadelete
157 subroutine qcspainit(qcspa,v7d,var, timei, timef, coordmin, coordmax, data_id_in,extremepath,spatialpath, &
159 dsne,usere,passworde,&
160 dsnspa,userspa,passwordspa,&
162 height2level,operation,categoryappend)
164 type(qcspatype),
intent(in out) :: qcspa
165 type (vol7d),
intent(in),
target:: v7d
166 character(len=*),
INTENT(in) :: var(:)
169 TYPE(geo_coord),
INTENT(inout),
optional :: coordmin,coordmax
171 TYPE(datetime),
INTENT(in),
optional :: timei, timef
172 integer,
intent(in),
optional,
target:: data_id_in(:,:,:,:,:)
173 character(len=*),
intent(in),
optional :: extremepath
174 character(len=*),
intent(in),
optional :: spatialpath
175 logical ,
intent(in),
optional :: height2level
176 character(len=*),
INTENT(in),
OPTIONAL :: categoryappend
177 character(len=*),
optional :: operation
180 type (vol7d_dballe) :: v7d_dballespa
181 character(len=*),
intent(in),
optional :: dsne
182 character(len=*),
intent(in),
optional :: usere
183 character(len=*),
intent(in),
optional :: passworde
184 character(len=*),
intent(in),
optional :: dsnspa
185 character(len=*),
intent(in),
optional :: userspa
186 character(len=*),
intent(in),
optional :: passwordspa
187 character(len=512) :: ldsnspa
188 character(len=512) :: luserspa
189 character(len=512) :: lpasswordspa
193 TYPE(vol7d_network) :: network
194 character(len=512) :: filepathspa
195 character(len=512) :: a_name
197 call l4f_launcher(a_name,a_name_append=trim(subcategoryspa)//
"."//trim(categoryappend))
198 qcspa%category=l4f_category_get(a_name)
202 nullify ( qcspa%data_id_in )
203 nullify ( qcspa%data_id_out )
209 qcspa%operation=optio_c(operation,20)
210 filepathspa=optio_c(spatialpath,512)
213 if (qcspa%operation /=
"gradient" .and. qcspa%operation /=
"run")
then
214 call l4f_category_log(qcspa%category,l4f_error,
"operation is wrong: "//qcspa%operation)
219 if (
present(data_id_in))
then
220 qcspa%data_id_in => data_id_in
224 call init(qcspa%qccli,v7d,var, timei, timef, data_id_in,&
225 macropath=cmiss, climapath=cmiss, extremepath=extremepath, &
227 dsncli=cmiss,dsnextreme=dsne,user=usere,password=passworde,&
229 height2level=height2level,categoryappend=categoryappend)
234 if (qcspa%operation ==
"run")
then
236 call init(network,
"qcspa-ndi")
239 call optio(dsnspa,ldsnspa)
240 call optio(userspa,luserspa)
241 call optio(passwordspa,lpasswordspa)
243 if (
c_e(filepathspa) .and. (
c_e(ldsnspa).or.
c_e(luserspa).or.
c_e(lpasswordspa)))
then
244 call l4f_category_log(qcspa%category,l4f_error,
"filepath defined together with dba options")
248 if (.not.
c_e(ldsnspa))
then
252 if (.not.
c_e(filepathspa))
then
253 filepathspa=get_package_filepath(
'qcspa-ndi.v7d', filetype_data)
256 if (
c_e(filepathspa))
then
258 select case (trim(lowercase(suffixname(filepathspa))))
262 call import(qcspa%clima,filename=filepathspa,unit=iuni)
267 call init(v7d_dballespa,file=.true.,filename=filepathspa,categoryappend=trim(a_name)//
".clima")
268 call import(v7d_dballespa,var=var, &
269 varkind=(/(
"r",i=1,
size(var))/),attr=(/
"*B33209"/),attrkind=(/
"b"/),network=network)
270 call copy(v7d_dballespa%vol7d,qcspa%clima)
271 call delete(v7d_dballespa)
276 "file type not supported (use .v7d or .bufr suffix only): "//trim(filepathspa))
281 call l4f_category_log(qcspa%category,l4f_warn,
"spatial clima volume not iniziatized: spatial QC will not be possible")
282 call init(qcspa%clima)
283 call raise_fatal_error()
290 call init(v7d_dballespa,dsn=ldsnspa,user=luserspa,password=lpasswordspa,write=.false.,&
291 file=.false.,categoryappend=trim(a_name)//
".spa")
293 call import(v7d_dballespa,var=var, &
294 varkind=(/("r
",i=1,size(var))/),attr=(/"*b33209
"/),attrkind=(/"b
"/),network=network)
295 call copy(v7d_dballespa%vol7d,qcspa%clima)
296 call delete(v7d_dballespa)
304 end subroutine qcspainit
307 subroutine qcspatri(qcspa,proj_type, lov, zone, xoff, yoff, &
308 longitude_south_pole, latitude_south_pole, angle_rotation, &
309 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
310 latin1, latin2, lad, projection_center_flag, &
311 ellips_smaj_axis, ellips_flatt, ellips_type)
313 type(qcspatype),intent(in out) :: qcspa !< Oggetto per il controllo climatico
314 CHARACTER(len=*),INTENT(in),OPTIONAL :: proj_type !< type of projection
315 DOUBLE PRECISION,INTENT(in),OPTIONAL :: lov !< line of view, also known as reference longitude or orientation of the grid (polar projections)
316 INTEGER,INTENT(in),OPTIONAL :: zone !< Earth zone (mainly for UTM), sets lov to the correct zone central meridian
317 DOUBLE PRECISION,INTENT(in),OPTIONAL :: xoff !< offset on x axis (false easting)
318 DOUBLE PRECISION,INTENT(in),OPTIONAL :: yoff !< offset on y axis (false northing)
319 DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_south_pole !< longitude of the southern pole of projection
320 DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_south_pole !< latitude of the southern pole of projection
321 DOUBLE PRECISION,INTENT(in),OPTIONAL :: angle_rotation !< angle of rotation of projection
322 DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_stretch_pole !< longitude of the pole of stretching
323 DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_stretch_pole !< latitude of the pole of stretching
324 DOUBLE PRECISION,INTENT(in),OPTIONAL :: stretch_factor !< stretching factor
325 DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin1 !< first standard latitude from main pole (Lambert)
326 DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin2 !< second standard latitude from main pole (Lambert)
327 DOUBLE PRECISION,INTENT(in),OPTIONAL :: lad !< latitude at which dx and dy (in m) are specified (Lambert, grib2 only)
328 INTEGER,INTENT(in),OPTIONAL :: projection_center_flag !< flag indicating which pole is represented
329 DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_smaj_axis !< Earth semi-major axis
330 DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_flatt !< Earth flattening
331 INTEGER,INTENT(in),OPTIONAL :: ellips_type !< number in the interval [1,nellips] indicating a predefined ellipsoid, alternative to the previous arguments
335 TYPE(geo_proj) :: geoproj
336 REAL(kind=fp_geo) :: lat(size(qcspa%v7d%ana)),lon(size(qcspa%v7d%ana))
337 CHARACTER(len=80) :: proj_type_l ! local type of projection
338 double precision :: lov_l, latin1_l,latin2_l
339 integer :: projection_center_flag_l
341 proj_type_l = optio_c(proj_type,80)
344 latin1_l = optio_d(latin1)
345 latin2_l = optio_d(latin2)
346 projection_center_flag_l=optio_l(projection_center_flag)
348 .not.
if ( c_e(proj_type_l)) then
349 proj_type_l = "lambert
"
353 projection_center_flag_l=1
356 geoproj = geo_proj_new(proj_type_l, lov_l, zone, xoff, yoff, &
357 longitude_south_pole, latitude_south_pole, angle_rotation, &
358 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
359 latin1_l, latin2_l, lad, projection_center_flag_l, &
360 ellips_smaj_axis, ellips_flatt, ellips_type)
362 call getval(qcspa%v7d%ana%coord, lon, lat)
364 !print*,"size
",size(lon),size(lat)
366 call proj(geoproj,lon,lat,qcspa%co%x,qcspa%co%y)
367 !print*,"size x y
",size(qcspa%x),size(qcspa%y)
368 !print*,qcspa%x,qcspa%y
371 status = triangles_compute(qcspa%co,qcspa%tri)
373 !qcspa%nt,qcspa%ipt,qcspa%nl,qcspa%ipl)
375 if (status /= 0) then
376 call l4f_category_log(qcspa%category,L4F_ERROR,"contng error status=
"//t2c(status))
380 end subroutine qcspatri
383 !>\brief Allocazioni di memoria
384 subroutine qcspaalloc(qcspa)
385 ! pseudo costruttore con distruttore automatico
387 type(qcspatype),intent(in out) :: qcspa !< Oggetto per il controllo climatico
392 ! se ti sei dimenticato di deallocare ci penso io
393 call qcspadealloc(qcspa)
396 !!$if (associated (qcspa%v7d%dativar%r )) then
397 !!$ nv=size(qcspa%v7d%dativar%r)
399 !!$ allocate(qcspa%valminr(nv),stat=istat)
400 !!$ istatt=istatt+istat
401 !!$ allocate(qcspa%valmaxr(nv),stat=istat)
402 !!$ istatt=istatt+istat
404 !!$ if (istatt /= 0) ier=1
408 if (associated(qcspa%data_id_in))then
409 sh=shape(qcspa%data_id_in)
410 allocate (qcspa%data_id_out(sh(1),sh(2),sh(3),sh(4),sh(5)),stat=istatt)
412 call l4f_category_log(qcspa%category,L4F_ERROR,"allocate error
")
413 call raise_error("allocate error
")
415 qcspa%data_id_out=imiss
419 if (associated(qcspa%v7d%ana))then
420 qcspa%ndp=size(qcspa%v7d%ana)
421 qcspa%tri = triangles_new(qcspa%ndp)
422 allocate(qcspa%co(qcspa%ndp))
425 end subroutine qcspaalloc
428 !>\brief Deallocazione della memoria
429 subroutine qcspadealloc(qcspa)
432 type(qcspatype),intent(in out) :: qcspa !< Oggetto per l controllo climatico
434 !!$if ( associated ( qcspa%valminr)) then
435 !!$ deallocate(qcspa%valminr)
438 !!$if ( associated ( qcspa%valmaxr)) then
439 !!$ deallocate(qcspa%valmaxr)
442 if (associated(qcspa%data_id_out)) then
443 deallocate (qcspa%data_id_out)
444 nullify (qcspa%data_id_out)
446 call delete(qcspa%tri)
447 if (associated(qcspa%co)) deallocate(qcspa%co)
449 end subroutine qcspadealloc
452 !>\brief Cancellazione
455 subroutine qcspadelete(qcspa)
456 ! decostruttore a mezzo
457 type(qcspatype),intent(in out) :: qcspa !< Oggetto per il controllo climatico
459 call qcspadealloc(qcspa)
461 call delete(qcspa%qccli)
466 call l4f_category_delete(qcspa%category)
469 end subroutine qcspadelete
472 à
!>\brief Controllo di Qualit spaziale.
473 èà
!!Questo il vero e proprio controllo di qualit spaziale.
475 SUBROUTINE quaconspa (qcspa,timetollerance,noborder,battrinv,battrcli,battrout,&
476 anamask,timemask,levelmask,timerangemask,varmask,networkmask)
479 à
type(qcspatype),intent(in out) :: qcspa !< Oggetto per il controllo di qualit
480 type(timedelta),intent(in) :: timetollerance !< time tollerance to compare nearest stations
481 logical,intent(in),optional :: noborder !< Exclude border from QC
482 character (len=10) ,intent(in),optional :: battrinv !< attributo invalidated in input
483 character (len=10) ,intent(in),optional :: battrcli !< attributo con la confidenza climatologica in input
484 character (len=10) ,intent(in),optional :: battrout !< attributo con la confidenza spaziale in output
485 logical ,intent(in),optional :: anamask(:) !< Filtro sulle anagrafiche
486 logical ,intent(in),optional :: timemask(:) !< Filtro sul tempo
487 logical ,intent(in),optional :: levelmask(:) !< Filtro sui livelli
488 logical ,intent(in),optional :: timerangemask(:) !< filtro sui timerange
489 logical ,intent(in),optional :: varmask(:) !< Filtro sulle variabili
490 logical ,intent(in),optional :: networkmask(:) !< Filtro sui network
492 !REAL(kind=fp_geo) :: lat,lon
494 integer :: indbattrinv,indbattrcli,indbattrout
495 logical :: anamaskl(size(qcspa%v7d%ana)), timemaskl(size(qcspa%v7d%time)), levelmaskl(size(qcspa%v7d%level)), &
496 timerangemaskl(size(qcspa%v7d%timerange)), varmaskl(size(qcspa%v7d%dativar%r)), networkmaskl(size(qcspa%v7d%network))
498 integer :: indana , indtime ,indlevel ,indtimerange ,inddativarr, indnetwork,indnet
499 integer :: indcana , indctime ,indclevel ,indctimerange ,indcdativarr, indcnetwork
500 real :: datoqui,datola,datila(size(qcspa%v7d%time)),climaquii, climaquif
501 !integer, allocatable :: indcanav(:)
503 !TYPE(vol7d_ana) :: ana
504 TYPE(datetime) :: time
505 !YPE(vol7d_level):: level
506 TYPE(vol7d_network):: network
507 type(timedelta) :: deltato,deltat
509 integer :: ivert(50),i,ipos,ineg,it,itrov,iv,ivb,kk,iindtime,grunit
510 double precision :: distmin=1000.d0,distscol=100000.d0
511 double precision :: dist,grad,gradmin
512 integer (kind=int_b) :: flag
513 !!$CHARACTER(len=vol7d_ana_lenident) :: ident
514 character(len=512) :: filename
518 !call qcspa_validate (qcspa)
520 if (size(qcspa%v7d%ana) < 3 ) then
521 call l4f_category_log(qcspa%category,L4F_WARN,"number of station < 3; do nothing
")
525 !localize optional parameter
526 if (present(battrinv))then
527 indbattrinv = index_c(qcspa%v7d%datiattr%b(:)%btable, battrinv)
529 indbattrinv = index_c(qcspa%v7d%datiattr%b(:)%btable, qcattrvarsbtables(1))
532 if (present(battrcli))then
533 indbattrcli = index_c(qcspa%v7d%datiattr%b(:)%btable, battrcli)
535 indbattrcli = index_c(qcspa%v7d%datiattr%b(:)%btable, qcattrvarsbtables(2))
538 if (present(battrout))then
539 indbattrout = index_c(qcspa%v7d%datiattr%b(:)%btable, battrout)
541 indbattrout = index_c(qcspa%v7d%datiattr%b(:)%btable, qcattrvarsbtables(4))
545 ! some checks on input
546 .or..or.
!if (indbattrinv <=0 indbattrcli <= 0 indbattrout <= 0 ) then
547 if (indbattrout <= 0 ) then
549 call l4f_category_log(qcspa%category,L4F_ERROR,"error finding attribute
index for output
")
554 !!$if (qcspa%operation == "gradient
") then
556 !!$ !check for gradient operation
557 .or.
!!$ if ( size(qcspa%v7d%level) > 1 &
558 .or.
!!$ size(qcspa%v7d%timerange) > 1 &
559 !!$ size(qcspa%v7d%dativar%r) > 1 ) then
560 !!$ call l4f_category_log(qcspa%category,L4F_ERROR,"gradient operation manage one level/timerange/var only
")
561 !!$ call raise_error()
566 ! set other local variable from optional parameter
567 if(present(anamask)) then
572 if(present(timemask)) then
577 if(present(levelmask)) then
578 levelmaskl = levelmask
582 if(present(timerangemask)) then
583 timerangemaskl = timerangemask
585 timerangemaskl = .true.
587 if(present(varmask)) then
592 if(present(networkmask)) then
593 networkmaskl = networkmask
595 networkmaskl = .true.
598 ! do not touch data that do not pass QC
599 qcspa%v7d%voldatiattrb(:,:,:,:,:,:,indbattrout)=ibmiss
601 !print *,"prima normalize
"
602 !print *,qcspa%v7d%voldatir
603 ! normalize data in space and time
604 call vol7d_normalize_data(qcspa%qccli)
605 !print *,"dopo normalize
"
606 !print *,qcspa%v7d%voldatir
613 ! compute some index for spatial clima
614 !! compute the conventional generic datetime
615 !!cyclicdt = cyclicdatetime_new(chardate="/////////
") !TMMGGhhmm
616 time=cyclicdatetime_to_conventional(cyclicdatetime_new(chardate="/////////
")) !TMMGGhhmm
617 !!call init(time, year=1007, month=1, day=1, hour=01, minute=01)
620 if (qcspa%operation == "run
") then
621 call init(network,"qcspa-ndi
")
622 !!indcana = firsttrue(qcspa%clima%ana == ana)
623 indcnetwork = index(qcspa%clima%network , network)
624 indctime = index(qcspa%clima%time , time)
627 do indtime=1,size(qcspa%v7d%time)
628 .not.
if (timemaskl(indtime)) cycle
629 call l4f_category_log(qcspa%category,L4F_INFO,&
630 "check time:
"//t2c(qcspa%v7d%time(indtime)) )
632 do indlevel=1,size(qcspa%v7d%level)
633 do indtimerange=1,size(qcspa%v7d%timerange)
634 do inddativarr=1,size(qcspa%v7d%dativar%r)
636 ind=index_c(spa_btable,qcspa%v7d%dativar%r(inddativarr)%btable)
638 if (qcspa%operation == "gradient
") then
639 ! open file to write gradient
641 filename=trim(to_char(qcspa%v7d%level(indlevel)))//&
642 "_
"//trim(to_char(qcspa%v7d%timerange(indtimerange)))//&
643 "_
"//trim(qcspa%v7d%dativar%r(inddativarr)%btable)//&
646 call l4f_category_log(qcspa%category,L4F_INFO,"try to open gradient file; filename below
")
647 call l4f_category_log(qcspa%category,L4F_INFO,filename)
649 inquire(file=filename, exist=exist)
652 if (grunit /= -1) then
653 !open (unit=grunit, file=t2c(timei)//"_
"//t2c(timef)//".grad
",STATUS='UNKNOWN', form='FORMATTED')
654 open (grunit, file=filename ,STATUS='UNKNOWN', form='FORMATTED',position='APPEND')
656 ! say we have to write header in file
657 .not.
if ( exist) then
658 call l4f_category_log(qcspa%category,L4F_INFO,"write header in gradient file
")
660 qcspa%v7d%level(indlevel), &
661 qcspa%v7d%timerange(indtimerange), &
662 qcspa%v7d%dativar%r(inddativarr)
667 do indnetwork=1,size(qcspa%v7d%network)
668 do indana=1,size(qcspa%v7d%ana)
670 !!$ call l4f_log(L4F_INFO,"index:
"// t2c(indana)//t2c(indnetwork)//t2c(indlevel)//&
671 !!$ t2c(indtimerange)//t2c(inddativarr)//t2c(indtime))
673 .not..or..not..or.
if ( anamaskl(indana) levelmaskl(indlevel) &
674 .not..or..not..or..not.
timerangemaskl(indtimerange) varmaskl(inddativarr) networkmaskl(indnetwork)) cycle
676 datoqui = qcspa%v7d%voldatir (indana ,indtime ,indlevel ,indtimerange ,inddativarr, indnetwork )
677 .not.
if ( c_e(datoqui)) cycle
680 if (indbattrinv > 0) then
681 if( invalidated(qcspa%v7d%voldatiattrb&
682 (indana,indtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrinv))) then
683 call l4f_category_log(qcspa%category,L4F_WARN,&
684 "it
's better to do a reform on ana to v7d after peeling, before spatial QC")
690 if (indbattrcli > 0) then
691 .not.
if( vdge(qcspa%v7d%voldatiattrb&
692 (indana,indtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrcli))) then
693 call l4f_category_log(qcspa%category,L4F_WARN,&
694 "It's better to do a reform on ana to v7d after peeling, before spatial qc
")
701 !!call init(anavar,"b07031
" )
702 !call init(anavar,"b07030
" )
704 !if (associated (qcspa%v7d%anavar%r)) then
705 ! indanavar = index(qcspa%v7d%anavar%r, anavar)
707 !if (indanavar <= 0 )cycle
708 !altezza= qcspa%v7d%volanar(indana,indanavar,indnetwork)
709 ! call spa_level(altezza,level)
711 if (qcspa%operation == "run
") then
713 indclevel = index(qcspa%clima%level , qcspa%v7d%level(indlevel))
714 indctimerange = index(qcspa%clima%timerange , qcspa%v7d%timerange(indtimerange))
716 ! attenzione attenzione TODO
717 è
! se leggo da bufr il default char e non reale
718 indcdativarr = index(qcspa%clima%dativar%r, qcspa%v7d%dativar%r(inddativarr))
722 call l4f_log(L4F_DEBUG,"index:
"// to_char(indctime)//to_char(indclevel)//&
723 to_char(indctimerange)//to_char(indcdativarr)//to_char(indcnetwork))
725 .or..or..or.
if ( indctime <= 0 indclevel <= 0 indctimerange <= 0 indcdativarr <= 0 &
726 .or.
indcnetwork <= 0 ) cycle
729 .and.
if (optio_log(noborder) any(indana == qcspa%tri%ipl(:3*qcspa%tri%nl:3))) cycle
731 ! ITROV e` il numero di triangoli in cui e` presente il dato
733 ! cicla per tutti i triangoli
735 ! se la stazione considerata e` in prima posizione
736 ! memorizza gli altri due vertici
737 .EQ.
IF(qcspa%tri%IPT(3*IT-2)INDANA)THEN
739 IVERT(2*ITROV)=qcspa%tri%IPT(3*IT)
740 IVERT(2*ITROV-1)=qcspa%tri%IPT(3*IT-1)
743 ! se la stazione considerata e` in seconda posizione
744 ! memorizza gli altri due vertici
745 .EQ.
IF(qcspa%tri%IPT(3*IT-1)INDANA)THEN
747 IVERT(2*ITROV)=qcspa%tri%IPT(3*IT)
748 IVERT(2*ITROV-1)=qcspa%tri%IPT(3*IT-2)
751 ! se la stazione considerata e` in terza posizione
752 ! memorizza gli altri due vertici
753 .EQ.
IF(qcspa%tri%IPT(3*IT)INDANA)THEN
755 IVERT(2*ITROV)=qcspa%tri%IPT(3*IT-1)
756 IVERT(2*ITROV-1)=qcspa%tri%IPT(3*IT-2)
760 ! ITROV ora diviene il numero di vertici nell'intorno
761 ! della stazione trovati
764 ! WRITE(*,*)'NUMERO VERTICI',ITROV
765 ! WRITE(*,*)'VERTICI TROVATI = ',IVERT
767 ! ordina i vettori dei vertici secondo valori decrescenti
769 call sort(ivert(:itrov))
773 .LT.
!!$ IF(IVERT(I)IVERT(KK))THEN
775 !!$ IVERT(KK)=IVERT(I)
780 ! toglie i valori doppi dal vettore dei vertici
783 .NE.
IF(IVERT(IV)IVERT(KK))THEN
788 .GT.
IF (IVITROV)IV=ITROV
790 ! WRITE(*,*)'NUMERO VERTICI puliti',IV
791 ! WRITE(*,*)'VERTICI PULITI = ',IVERT
793 ! inizia il controllo sulla stazione testando i gradienti
794 ! WRITE(*,*)'STAZIONE ',INDANA
798 gradmin=huge(gradmin)
800 !find the nearest data in time
801 datola = qcspa%v7d%voldatir (ivert(i) ,indtime ,indlevel ,indtimerange ,inddativarr, indnetwork )
804 if (indbattrinv > 0) then
805 if( invalidated(qcspa%v7d%voldatiattrb&
806 (ivert(i),indtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrinv))) then
812 if (indbattrcli > 0) then
813 .not.
if( vdge(qcspa%v7d%voldatiattrb&
814 (ivert(i),indtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrcli))) then
819 ! if we do not have data from the same network at the same time
820 ! here we search for the first data found (nearest in time) looping over the network
821 ! we do not have priority for network to take in account
823 deltato=timedelta_miss
824 do indnet=1, size(qcspa%v7d%network)
825 datila = qcspa%v7d%voldatir (ivert(i) ,: ,indlevel ,indtimerange ,inddativarr, indnet )
826 do iindtime=1,size(qcspa%v7d%time)
827 .not.
if ( c_e(datila(iindtime))) cycle
829 if (indbattrinv > 0 ) then
830 if (invalidated(qcspa%v7d%voldatiattrb&
831 (ivert(i),iindtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrinv))) cycle
834 if (indbattrcli > 0 )then
835 .not.
if ( vdge(qcspa%v7d%voldatiattrb&
836 (ivert(i),iindtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrcli))) cycle
839 if (iindtime < indtime) then
840 deltat=qcspa%v7d%time(indtime)-qcspa%v7d%time(iindtime)
841 else if (iindtime >= indtime) then
842 deltat=qcspa%v7d%time(iindtime)-qcspa%v7d%time(indtime)
845 .or..not..and.
if ((deltat < deltato c_e(deltato)) deltat <= timetollerance ) then
846 datola = datila(iindtime)
853 .NOT.
IF(C_E(datola)) cycle
854 ! distanza tra le due stazioni
855 dist = DISTANZA (qcspa%co(INDANA),qcspa%co(IVERT(I)))
857 call l4f_category_log(qcspa%category,L4F_ERROR,"distance from two station == 0.
")
862 call l4f_log (L4F_DEBUG,"distanza:
"//t2c(dist))
865 dist=max(dist,distmin)
867 ! se la distanza supera distscol, stazioni scorrelate - salta -
868 if (dist > distscol) cycle
870 ! valore del gradiente nella direzione delle due stazioni
871 GRAD=(datoqui-datola)/(DIST)
872 IF (GRAD >= 0.d0) Ipos=Ipos+1 ! se il gradiente e` positivo incrementa il contatore di positivi
873 IF (GRAD <= 0.d0) Ineg=Ineg+1 ! se il gradiente e` negativo incrementa il contatore di negativi
875 gradmin=min(gradmin,abs(grad))
880 call l4f_log (L4F_DEBUG,"ivb:
"//t2c(ivb))
883 IF(IVB < 3) cycle ! do nothing if valid gradients < 3
885 .or.
IF (ipos == ivb ineg == ivb)THEN ! se tutti i gradienti sono dello stesso segno
887 gradmin=sign(gradmin,dble(ipos-ineg))
889 if (qcspa%operation == "gradient
") then
890 write(grunit,*)gradmin
893 ! we normalize gradmin or denormalize climaqui after
894 ! gradmin=gradmin*spa_a(ind) + spa_b(ind)
898 call l4f_log (L4F_DEBUG,"gradmin:
"//t2c(gradmin))
903 È
!ATTENZIONE TODO : inddativarr UNA GRANDE SEMPLIFICAZIONE NON VERA SE TIPI DI DATO DIVERSI !!!!
904 if (qcspa%operation == "run
") then
906 do indcana=1,size(qcspa%clima%ana)
907 climaquii=(qcspa%clima%voldatir(indcana &
908 ,indctime,indclevel,indctimerange,indcdativarr,indcnetwork)&
909 - spa_b(ind))/spa_a(ind) ! denormalize
911 ! the last interval have climaquii == climaquif
912 ! the check will be the last extreme to +infinite
913 climaquif=(qcspa%clima%voldatir(min(indcana+1,size(qcspa%clima%ana)) &
914 ,indctime,indclevel,indctimerange,indcdativarr,indcnetwork)&
915 - spa_b(ind))/spa_a(ind) ! denormalize
918 call l4f_log (L4F_DEBUG,"climaquii:
"//t2c(climaquii))
919 call l4f_log (L4F_DEBUG,"climaquif:
"//t2c(climaquif))
922 .and.
if ( c_e(climaquii) c_e(climaquif )) then
924 .and..or.
if ( (gradmin >= climaquii gradmin < climaquif) &
925 .and..or.
(indcana == 1 gradmin < climaquii) &
926 .and.
(indcana == size(qcspa%clima%ana) gradmin >= climaquif) ) then
928 flag=qcspa%clima%voldatiattrb(indcana,indctime,indclevel,indctimerange,indcdativarr,indcnetwork,1)
930 if ( associated ( qcspa%data_id_in)) then
932 call l4f_log (L4F_DEBUG,"id:
"//t2c(&
933 qcspa%data_id_in(indana,indtime,indlevel,indtimerange,indnetwork)))
935 qcspa%data_id_out(indana,indtime,indlevel,indtimerange,indnetwork)=&
936 qcspa%data_id_in(indana,indtime,indlevel,indtimerange,indnetwork)
942 call l4f_log (L4F_INFO,"datoqui:
"//t2c(datoqui))
943 call l4f_log (L4F_INFO,"flag qcspa:
"//t2c(flag))
950 if (qcspa%operation == "run
") then
951 à
!TODO controllare se flag = missing comporta rimozione della precedente flag; risposta: SI quando sar chiusa https://github.com/ARPA-SIMC/dballe/issues/44
952 qcspa%v7d%voldatiattrb( indana, indtime, indlevel, indtimerange, inddativarr, indnetwork, indbattrout)=flag
957 if (qcspa%operation == "gradient
") then
966 !!$print*,"risultato
"
967 !!$print *,qcspa%v7d%voldatiattrb(:,:,:,:,:,:,indbattrout)
968 !!$print*,"fine risultato
"
974 elemental double precision function DISTANZA (co1,co2)
975 type(xy), intent(in) :: co1,co2
978 distanza = sqrt((co2%x-co1%x)**2 + (co2%y-co1%y)**2)
980 end function DISTANZA
982 end subroutine quaconspa
988 !> \example v7d_qcspa.F90
989 !! Sample program for module qcspa
Import one or more geo_coordvect objects from a plain text file or for a file in ESRI/Shapefile forma...
Copy an object, creating a fully new instance.
Emit log message for a category with specific priority.
Check for problems return 0 if all check passed print diagnostics with log4f.
Utilities for CHARACTER variables.
Classi per la gestione delle coordinate temporali.
Utilities for managing files.
Classes for handling georeferenced sparse points in geographical corodinates.
classe per la gestione del logging
Utilities and defines for quality control.
Controllo di qualità climatico.
Controllo di qualità spaziale.
Space utilities, derived from NCAR software.
Classe per la gestione di un volume completo di dati osservati.
classe per import ed export di volumi da e in DB-All.e
Oggetto principale per il controllo di qualità