125 character (len=255),
parameter:: subcategory=
"QCcli"
127 integer,
parameter :: cli_nlevel=10
129 integer,
parameter :: cli_level1(cli_nlevel) = (/-100,100,250,500,750,1000,1250,1500,1750,2000/)
131 integer,
parameter :: cli_level2(cli_nlevel) = (/100,250,500,750,1000,1250,1500,1750,2000,2250/)
136 type (vol7d),
pointer :: v7d
137 type (vol7d) :: clima
138 type (vol7d) :: extreme
139 integer,
pointer :: data_id_in(:,:,:,:,:)
140 integer,
pointer :: data_id_out(:,:,:,:,:)
141 integer,
pointer :: in_macroa(:)
143 logical :: height2level
150 module procedure qccliinit
155 module procedure qcclialloc
160 module procedure qcclidelete
165 vol7d_normalize_data, quaconcli, cli_level, cli_level_generate, &
166 qc_compute_percentile, qc_compute_normalizeddensityindex
173 subroutine qccliinit(qccli,v7d,var, timei, timef, data_id_in,&
174 macropath, climapath, extremepath, &
176 dsncli,dsnextreme,user,password,&
178 height2level,categoryappend)
180 type(qcclitype),
intent(in out) :: qccli
181 type (vol7d),
intent(in),
target:: v7d
182 character(len=*),
INTENT(in) :: var(:)
184 TYPE(datetime),
INTENT(in),
optional :: timei, timef
185 integer,
intent(in),
optional,
target:: data_id_in(:,:,:,:,:)
186 character(len=*),
intent(in),
optional :: macropath
187 character(len=*),
intent(in),
optional :: climapath
188 character(len=*),
intent(in),
optional :: extremepath
189 logical ,
intent(in),
optional :: height2level
190 character(len=*),
INTENT(in),
OPTIONAL :: categoryappend
196 type (vol7d_dballe) :: v7d_dballecli
197 character(len=*),
intent(in),
optional :: dsncli
198 character(len=*),
intent(in),
optional :: user
199 character(len=*),
intent(in),
optional :: password
200 character(len=512) :: ldsncli
201 character(len=512) :: luser
202 character(len=512) :: lpassword
203 type (vol7d_dballe) :: v7d_dballeextreme
204 character(len=*),
intent(in),
optional :: dsnextreme
205 character(len=512) :: ldsnextreme
206 TYPE(datetime) :: ltimei, ltimef
207 integer :: yeari, yearf, monthi, monthf, dayi, dayf,&
208 houri, minutei, mseci, hourf, minutef, msecf
212 character(len=512) :: filepath
213 character(len=512) :: filepathclima
214 character(len=512) :: filepathextreme
215 character(len=512) :: a_name
216 TYPE(geo_coord) :: lcoordmin,lcoordmax
217 TYPE(vol7d_network) :: network
220 DOUBLE PRECISION :: lon, lat
221 TYPE(georef_coord) :: point
222 TYPE(arrayof_georef_coord_array) :: macroa
225 call l4f_launcher(a_name,a_name_append=trim(subcategory)//
"."//trim(categoryappend))
226 qccli%category=l4f_category_get(a_name)
228 qccli%height2level=optio_log(height2level)
240 nullify ( qccli%in_macroa )
241 nullify ( qccli%data_id_in )
242 nullify ( qccli%data_id_out )
247 if (
present(data_id_in))
then
248 qccli%data_id_in => data_id_in
251 if (qccli%height2level)
then
253 filepath=get_package_filepath(
'share:sup_macroaree.shp', filetype_data)
256 filepath=get_package_filepath(
'share:ens_v8_ll.shp', filetype_data)
259 if (
present(macropath))
then
260 if (
c_e(macropath))
then
265 CALL import(macroa, shpfile=filepath)
267 call optio(climapath,filepathclima)
268 call optio(extremepath,filepathextreme)
272 call init(network,
"qcclima-ndi")
276 if (
present (timei))
then
278 ltimei=ltimei+timedelta_new(minute=30)
281 if (
present (timef))
then
283 ltimef=ltimef+timedelta_new(minute=30)
286 call getval(ltimei, year=yeari, month=monthi, day=dayi, hour=houri, minute=minutei, msec=mseci)
287 call getval(ltimef, year=yearf, month=monthf, day=dayf, hour=hourf, minute=minutef, msec=msecf)
289 if (
c_e(yeari) .and.
c_e(yearf) .and. yeari == yearf .and. monthi == monthf )
then
293 ltimei=cyclicdatetime_to_conventional(cyclicdatetime_new(month=monthi))
294 ltimef=cyclicdatetime_to_conventional(cyclicdatetime_new(month=monthf))
303 call optio(dsncli,ldsncli)
304 call optio(user,luser)
305 call optio(password,lpassword)
307 if ((
c_e(filepathclima) .or.
c_e(filepathextreme)) .and. (
c_e(ldsncli).or.
c_e(luser).or.
c_e(lpassword)))
then
308 call l4f_category_log(qccli%category,l4f_error,
"climapath or extremepath defined together with dba options")
312 if (.not.
c_e(ldsncli))
then
316 if (.not.
c_e(filepathclima))
then
317 filepathclima=get_package_filepath(
'qcclima-ndi.v7d', filetype_data)
320 if (
c_e(filepathclima))
then
322 select case (trim(lowercase(suffixname(filepathclima))))
326 call import(qccli%clima,filename=filepathclima,unit=iuni)
331 call init(v7d_dballecli,file=.true.,filename=filepathclima,categoryappend=trim(a_name)//
".clima")
333 call import(v7d_dballecli,var=var,coordmin=lcoordmin, coordmax=lcoordmax, timei=ltimei, timef=ltimef, &
334 varkind=(/(
"r",i=1,
size(var))/),attr=(/
"*B33209"/),attrkind=(/
"b"/),network=network)
335 call copy(v7d_dballecli%vol7d,qccli%clima)
341 "file type not supported (use .v7d or .bufr suffix only): "//trim(filepathclima))
346 call l4f_category_log(qccli%category,l4f_warn,
"clima volume not iniziatized: clima QC will not be possible")
355 call init(v7d_dballecli,dsn=ldsncli,user=luser,password=lpassword,write=.false.,&
356 file=.false.,categoryappend=trim(a_name)//
".clima")
358 call import(v7d_dballecli,var=var,coordmin=lcoordmin, coordmax=lcoordmax, timei=ltimei, timef=ltimef, &
359 varkind=(/("r
",i=1,size(var))/),attr=(/"*b33209
"/),attrkind=(/"b
"/),network=network)
360 call copy(v7d_dballecli%vol7d,qccli%clima)
361 call delete(v7d_dballecli)
370 .and..and..and.
if ( c_e(yeari) c_e(yearf) yeari == yearf monthi == monthf ) then
372 .and..and..and.
if ( dayi == dayf houri == hourf minutei == minutef mseci == msecf ) then
374 ltimei=cyclicdatetime_to_conventional(cyclicdatetime_new(month=monthi, hour=houri))
375 ltimef=cyclicdatetime_to_conventional(cyclicdatetime_new(month=monthf, hour=hourf))
379 ltimei=cyclicdatetime_to_conventional(cyclicdatetime_new(month=monthi, hour=00))
380 ltimef=cyclicdatetime_to_conventional(cyclicdatetime_new(month=monthf, hour=23))
385 ! if you span years or months or days I read all the climat dataset (should be optimized not so easy)
391 call init(network,"qcclima-perc
")
392 call optio(dsnextreme,ldsnextreme)
394 .not.
if ( c_e(ldsnextreme)) then
398 .not.
if ( c_e(filepathextreme)) then
399 filepathextreme=get_package_filepath('qcclima-extreme.v7d', filetype_data)
402 if (c_e(filepathextreme)) then
404 select case (trim(lowercase(suffixname(filepathextreme))))
408 call import(qccli%extreme,filename=filepathextreme,unit=iuni)
413 call init(v7d_dballeextreme,file=.true.,filename=filepathextreme,categoryappend=trim(a_name)//".climaextreme
")
414 !call import(v7d_dballeextreme)
415 call import(v7d_dballeextreme,var=var,coordmin=lcoordmin, coordmax=lcoordmax, timei=ltimei, timef=ltimef, &
416 varkind=(/("r
",i=1,size(var))/),attr=(/qcattrvarsbtables(2)/),attrkind=(/"b
"/),network=network)
417 call copy(v7d_dballeextreme%vol7d,qccli%extreme)
418 call delete(v7d_dballeextreme)
423 if (c_e(filepathextreme)) then
424 call l4f_category_log(qccli%category,L4F_ERROR,&
425 "file type not supported (user .v7d or .bufr suffix only):
"//trim(filepathextreme))
431 call l4f_category_log(qccli%category,L4F_WARN,"extreme volume not iniziatized: qc or normalize data will not be possible
")
432 ! call raise_fatal_error()
433 call init(qccli%extreme)
440 call l4f_category_log(qccli%category,L4F_DEBUG,"init v7d_dballeextreme
")
441 call init(v7d_dballeextreme,dsn=ldsnextreme,user=luser,password=lpassword,&
442 write=.false.,file=.false.,categoryappend=trim(a_name)//".climaextreme
")
443 call l4f_category_log(qccli%category,L4F_DEBUG,"import v7d_dballeextreme
")
445 call import(v7d_dballeextreme,var=var,coordmin=lcoordmin, coordmax=lcoordmax, timei=ltimei, timef=ltimef, &
446 varkind=(/("r
",i=1,size(var))/),attr=(/qcattrvarsbtables(2)/),attrkind=(/"b
"/),network=network)
447 call copy(v7d_dballeextreme%vol7d,qccli%extreme)
448 call delete(v7d_dballeextreme)
457 call qcclialloc(qccli)
460 ! valuto in quale macroarea sono le stazioni
462 !!$IF (macroa%arraysize <= 0) THEN
463 !!$ CALL l4f_category_log(qccli%category,L4F_ERROR,"maskgen: poly parameter missing or empty
")
464 !!$ CALL raise_fatal_error()
467 if (associated(qccli%in_macroa)) then
468 qccli%in_macroa = imiss
470 DO i = 1, SIZE(qccli%v7d%ana)
471 ! temporary, improve!!!!
472 CALL getval(qccli%v7d%ana(i)%coord,lon=lon,lat=lat)
473 point = georef_coord_new(x=lon, y=lat)
474 DO j = 1, macroa%arraysize
475 IF (inside(point, macroa%array(j))) THEN
476 qccli%in_macroa(i) = j
486 end subroutine qccliinit
489 !>\brief Allocazioni di memoria
490 subroutine qcclialloc(qccli)
491 ! pseudo costruttore con distruttore automatico
493 type(qcclitype),intent(in out) :: qccli !< Oggetto per il controllo climatico
498 ! se ti sei dimenticato di deallocare ci penso io
499 call qcclidealloc(qccli)
502 !!$if (associated (qccli%v7d%dativar%r )) then
503 !!$ nv=size(qccli%v7d%dativar%r)
505 !!$ allocate(qccli%valminr(nv),stat=istat)
506 !!$ istatt=istatt+istat
507 !!$ allocate(qccli%valmaxr(nv),stat=istat)
508 !!$ istatt=istatt+istat
510 !!$ if (istatt /= 0) ier=1
514 if (associated (qccli%v7d%ana )) then
515 allocate (qccli%in_macroa(size(qccli%v7d%ana )),stat=istatt)
516 if (istatt /= 0) then
517 call l4f_category_log(qccli%category,L4F_ERROR,"allocate error
")
518 call raise_error("allocate error
")
522 if (associated(qccli%data_id_in))then
523 sh=shape(qccli%data_id_in)
524 allocate (qccli%data_id_out(sh(1),sh(2),sh(3),sh(4),sh(5)),stat=istatt)
526 call l4f_category_log(qccli%category,L4F_ERROR,"allocate error
")
527 call raise_error("allocate error
")
529 qccli%data_id_out=imiss
535 end subroutine qcclialloc
538 !>\brief Deallocazione della memoria
540 subroutine qcclidealloc(qccli)
543 type(qcclitype),intent(in out) :: qccli !< Oggetto per l controllo climatico
545 !!$if ( associated ( qccli%valminr)) then
546 !!$ deallocate(qccli%valminr)
549 !!$if ( associated ( qccli%valmaxr)) then
550 !!$ deallocate(qccli%valmaxr)
553 if (associated (qccli%in_macroa)) then
554 deallocate (qccli%in_macroa)
557 if (associated(qccli%data_id_out))then
558 deallocate (qccli%data_id_out)
562 end subroutine qcclidealloc
565 !>\brief Cancellazione
568 subroutine qcclidelete(qccli)
569 ! decostruttore a mezzo
570 type(qcclitype),intent(in out) :: qccli !< Oggetto per l controllo climatico
572 call qcclidealloc(qccli)
574 call delete(qccli%clima)
575 call delete(qccli%extreme)
578 call l4f_category_delete(qccli%category)
581 end subroutine qcclidelete
585 !> Modulo 1: Calcolo dei parametri di normalizzazione dei dati
586 °°°
!! I parametri di normalizzazione sono il 25, il 50 e il 75 percentile
588 !!oppure (p15.87,p50.,p84.13)
589 !! Tali parametri verranno calcolati per ogni mese, per ogni ora, per ogni area.
590 !! Modulo 2: Normalizzazione dati
591 Ã
!! Ciascun dato D verr normalizzato come segue:
592 !! DN = (D-p50)*2/(p75-p25)
593 !! oppure DN = (D-p50)*2/(p16-p84)
594 è
!! dove DN il valore normalizzato.
595 !! La scelta dei parametri di normalizzazione dipende dal mese, dall'ora,
597 SUBROUTINE vol7d_normalize_data(qccli)
599 TYPE(qcclitype),INTENT(inout) :: qccli !< volume providing data to be computed, it is modified by the method
600 !character (len=10) ,intent(in),optional :: battrinv !< attributo invalidated in input/output
602 real :: datoqui, perc25, perc50,perc75
603 integer :: indana , indtime ,indlevel ,indtimerange ,inddativarr, indnetwork
604 integer :: indcana, indvar,indctime,indclevel,indctimerange,indcdativarr,indcnetwork
605 !integer :: indbattrinv
606 integer :: iclv(size(qccli%v7d%ana))
608 character(len=1) :: type
609 TYPE(vol7d_var) :: var
610 TYPE(vol7d_ana) :: ana
611 TYPE(datetime) :: time, nintime
612 TYPE(vol7d_level):: level
613 integer :: mese, ora, desc, iarea, k
616 character(len=1) :: mycanc, canc = "#
"
618 CHARACTER(len=vol7d_ana_lenident) :: ident
622 !!$if (associated(qccli%v7d%dativarattr%b))then
623 !!$ if (present(battrinv))then
624 !!$ indbattrinv = index_c(qccli%v7d%dativarattr%b(:)%btable, battrinv)
626 !!$ indbattrinv = index_c(qccli%v7d%dativarattr%b(:)%btable, qcattrvarsbtables(1))
631 .not.
if ( associated(qccli%extreme%voldatir))then
632 call l4f_category_log(qccli%category,L4F_WARN,"extreme data not associated: normalize data not possible
")
633 qccli%v7d%voldatir=rmiss
634 ! call raise_fatal_error()
639 if (qccli%height2level) then
640 call init(var, btable="b07030
") ! height
643 indvar = index(qccli%v7d%anavar, var, type=type)
645 do indana=1,size(qccli%v7d%ana)
648 ! here we take the height fron any network (the first network win)
649 do indnetwork=1,size(qccli%v7d%network)
651 if( indvar > 0 ) then
654 height=realdat(qccli%v7d%volanad(indana,indvar,indnetwork),qccli%v7d%anavar%d(indvar))
656 height=realdat(qccli%v7d%volanar(indana,indvar,indnetwork),qccli%v7d%anavar%r(indvar))
658 height=realdat(qccli%v7d%volanai(indana,indvar,indnetwork),qccli%v7d%anavar%i(indvar))
660 height=realdat(qccli%v7d%volanab(indana,indvar,indnetwork),qccli%v7d%anavar%b(indvar))
662 height=realdat(qccli%v7d%volanac(indana,indvar,indnetwork),qccli%v7d%anavar%c(indvar))
666 if (c_e(height)) exit
669 if (c_e(height)) then
670 .and.
iclv(indana)=firsttrue(cli_level1 <= height height <= cli_level2 )
676 call l4f_category_log(qccli%category,L4F_DEBUG, 'vol7d_normalize_data height has value '//t2c(height,"missing
"))
677 call l4f_category_log(qccli%category,L4F_DEBUG, 'for indana having number '//t2c(indana)//&
678 ' iclv has value '//t2c(iclv(indana),"missing
"))
685 .not.
if ( associated (qccli%in_macroa)) then
686 call l4f_category_log(qccli%category,L4F_WARN,"macroarea data not iniziatized: normalize data not possible
")
687 qccli%v7d%voldatir=rmiss
688 ! call raise_fatal_error()
692 .not.
if ( associated(qccli%extreme%voldatir)) then
693 call l4f_category_log(qccli%category,L4F_WARN,"qccli%extreme%voldatir not iniziatized: normalize data not possible
")
694 qccli%v7d%voldatir=rmiss
695 ! call raise_fatal_error()
699 do indana=1,size(qccli%v7d%ana)
700 iarea= qccli%in_macroa(indana)
702 do indnetwork=1,size(qccli%v7d%network)
703 do indlevel=1,size(qccli%v7d%level)
704 do indtimerange=1,size(qccli%v7d%timerange)
705 do inddativarr=1,size(qccli%v7d%dativar%r)
706 do indtime=1,size(qccli%v7d%time)
708 datoqui = qccli%v7d%voldatir (indana ,indtime ,indlevel ,indtimerange ,inddativarr, indnetwork )
710 .not.
if ( c_e(datoqui)) cycle
712 .not.
if ( c_e(iarea)) then
713 qccli%v7d%voldatir (indana ,indtime ,indlevel ,indtimerange ,&
714 inddativarr, indnetwork ) = rmiss
718 !!$ if (indbattrinv > 0) then
719 !!$ if( invalidated(qccli%v7d%voldatiattrb&
720 !!$ (indana,indtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrinv))) cycle
723 nintime=qccli%v7d%time(indtime)+timedelta_new(minute=30)
724 CALL getval(nintime, month=mese, hour=ora)
726 time=cyclicdatetime_to_conventional(cyclicdatetime_new(month=mese, hour=ora))
727 !call init(time, year=1001, month=mese, day=1, hour=ora, minute=01)
729 level=qccli%v7d%level(indlevel)
733 !indcana = firsttrue(qccli%extreme%ana == ana)
735 indctime = index(qccli%extreme%time , time)
736 indclevel = index(qccli%extreme%level , level)
737 indctimerange = index(qccli%extreme%timerange , qccli%v7d%timerange(indtimerange))
739 ! attenzione attenzione TODO
740 è
! se leggo da bufr il default char e non reale
742 indcdativarr = index(qccli%extreme%dativar%r, qccli%v7d%dativar%r(inddativarr))
744 !!$ print *,"dato
",qccli%v7d%timerange(indtimerange)
745 !!$ print *,"extreme
",qccli%extreme%timerange
746 call l4f_log(L4F_DEBUG,"normalize data
index:
"// to_char(indctime)//to_char(indclevel)//&
747 to_char(indctimerange)//to_char(indcdativarr)//to_char(indcnetwork))
749 .or..or..or..or.
!if (indcana <= 0 indctime <= 0 indclevel <= 0 indctimerange <= 0 indcdativarr <= 0 &
750 .or.
! indcnetwork <= 0 ) cycle
751 .or..or..or.
if (indctime <= 0 indclevel <= 0 indctimerange <= 0 indcdativarr <= 0 &
752 .or.
indcnetwork <= 0 ) then
753 qccli%v7d%voldatir (indana ,indtime ,indlevel ,indtimerange ,&
754 inddativarr, indnetwork ) = rmiss
758 ! find percentile in volume
765 if (qccli%height2level) then
771 .not.
if ( c_e(k)) then
772 qccli%v7d%voldatir (indana ,indtime ,indlevel ,indtimerange ,&
773 inddativarr, indnetwork ) = rmiss
778 write(ident,'("#
",i2.2,2i3.3)')k,iarea,desc ! macro-area e descrittore
780 .and..and..and.
!!$ if (indana == 1996 indtime == 855 indlevel == 1 indtimerange==1 &
781 .and..and.
!!$ inddativarr==1 indnetwork==1 ) then
782 !!$ call l4f_category_log(qccli%category,L4F_WARN,"mmmmmmmmmmmmm:
"//t2c(datoqui)//ident)
785 call init(ana,lat=0.0_fp_geo,lon=0.0_fp_geo,ident=ident)
786 indcana=index(qccli%extreme%ana,ana)
787 call l4f_log(L4F_DEBUG,"normalize data index25:
"//to_char(k)//to_char(indcana)//&
788 to_char(indctime)//to_char(indclevel)//&
789 to_char(indctimerange)//to_char(indcdativarr)//to_char(indcnetwork))
790 if (indcana > 0 )then
791 perc25=qccli%extreme%voldatir(indcana,indctime,indclevel,&
792 indctimerange,indcdativarr,indcnetwork)
796 write(ident,'("#
",i2.2,2i3.3)')k,iarea,desc ! macro-area e descrittore
797 call init(ana,lat=0.d0,lon=0.d0,ident=ident)
798 indcana=index(qccli%extreme%ana,ana)
799 call l4f_log(L4F_DEBUG,"normalize data index50:
"//to_char(k)//to_char(indcana)//&
800 to_char(indctime)//to_char(indclevel)//&
801 to_char(indctimerange)//to_char(indcdativarr)//to_char(indcnetwork))
802 if (indcana > 0 )then
803 perc50=qccli%extreme%voldatir(indcana,indctime,indclevel,&
804 indctimerange,indcdativarr,indcnetwork)
808 write(ident,'("#
",i2.2,2i3.3)')k,iarea,desc ! macro-area e descrittore
809 call init(ana,lat=0.d0,lon=0.d0,ident=ident)
810 indcana=index(qccli%extreme%ana,ana)
811 call l4f_log(L4F_DEBUG,"normalize data index75:
"//to_char(k)//to_char(indcana)//&
812 to_char(indctime)//to_char(indclevel)//&
813 to_char(indctimerange)//to_char(indcdativarr)//to_char(indcnetwork))
814 if (indcana > 0 )then
815 perc75=qccli%extreme%voldatir(indcana,indctime,indclevel,&
816 indctimerange,indcdativarr,indcnetwork)
819 .and..and..and.
if ( c_e(perc25) c_e(perc50) c_e(perc75) &
820 ABS( perc75 - perc25 ) >= SPACING( MAX(ABS(perc75),ABS(perc25))))then
823 call l4f_log(L4F_DEBUG,"normalize data dato in :
"//t2c(datoqui))
824 datoqui = (datoqui - perc50) / (perc75 - perc25) + &
825 base_value(qccli%v7d%dativar%r(inddativarr)%btable)
826 call l4f_log(L4F_DEBUG,"normalize data dato out :
"//t2c(datoqui))
834 qccli%v7d%voldatir (indana ,indtime ,indlevel ,indtimerange ,&
835 inddativarr, indnetwork ) = datoqui
843 read (qccli%v7d%ana(indana)%ident,'(a1,i2.2,2i3.3)') mycanc, clev, iarea, desc
844 if (mycanc == canc) then
847 write (qccli%v7d%ana(indana)%ident,'(a1,i2.2,2i3.3)') canc, clev, iarea, desc
852 end SUBROUTINE vol7d_normalize_data
855 real function base_value(btable)
856 character (len=10) ,intent(in):: btable
858 character (len=10) :: btables(1) =(/"b12101
"/)
859 real :: base_values(1) =(/273.15 /)
862 ind = index_c(btables,btable)
865 base_value = base_values(ind)
867 call l4f_log(L4F_WARN,"modqccli_base_value: variable
"//btable//" do not have base value
")
873 end function base_value
876 Ã
!>\brief Controllo di Qualit climatico.
877 èÃ
!!Questo il vero e proprio controllo di qualit climatico.
878 !!Avendo a disposizione un volume dati climatico
879 !!contenente i percentili suddivisi per area, altezza sul livello
880 !!del mare, per mese dell'anno viene selezionato il percentile e sulla base di questo
881 !!vengono assegnate le opportune confidenze.
882 SUBROUTINE quaconcli (qccli,battrinv,battrout,&
883 anamask,timemask,levelmask,timerangemask,varmask,networkmask)
886 Ã
type(qcclitype),intent(in out) :: qccli !< Oggetto per il controllo di qualit
887 character (len=10) ,intent(in),optional :: battrinv !< attributo invalidated in input/output
888 character (len=10) ,intent(in),optional :: battrout !< attributo con la confidenza climatologica in output
889 logical ,intent(in),optional :: anamask(:) !< Filtro sulle anagrafiche
890 logical ,intent(in),optional :: timemask(:) !< Filtro sul tempo
891 logical ,intent(in),optional :: levelmask(:) !< Filtro sui livelli
892 logical ,intent(in),optional :: timerangemask(:) !< filtro sui timerange
893 logical ,intent(in),optional :: varmask(:) !< Filtro sulle variabili
894 logical ,intent(in),optional :: networkmask(:) !< Filtro sui network
896 CHARACTER(len=vol7d_ana_lenident) :: ident
897 REAL(kind=fp_geo) :: latc,lonc
900 integer :: indbattrinv,indbattrout
901 logical :: anamaskl(size(qccli%v7d%ana)), timemaskl(size(qccli%v7d%time)), levelmaskl(size(qccli%v7d%level)), &
902 timerangemaskl(size(qccli%v7d%timerange)), varmaskl(size(qccli%v7d%dativar%r)), networkmaskl(size(qccli%v7d%network))
904 integer :: indana , indtime ,indlevel ,indtimerange ,inddativarr, indnetwork
905 integer :: indcana, indctime,indclevel,indctimerange,indcdativarr,indcnetwork
906 real :: datoqui,climaquii,climaquif, extremequii,extremequif,perc25,perc50,perc75
907 integer :: iarea,desc,indn,indvar,k
909 !integer, allocatable :: indcanav(:)
912 TYPE(vol7d_network) :: network
913 TYPE(vol7d_ana) :: ana
914 TYPE(datetime) :: time, nintime
915 TYPE(vol7d_level):: level
916 type(vol7d_var) :: anavar
917 integer :: iclv(size(qccli%v7d%ana))
918 character(len=1) :: type
921 !call qccli_validate (qccli)
924 if (associated(qccli%v7d%datiattr%b))then
925 if (present(battrinv))then
926 indbattrinv = index_c(qccli%v7d%datiattr%b(:)%btable, battrinv)
928 indbattrinv = index_c(qccli%v7d%datiattr%b(:)%btable, qcattrvarsbtables(1))
932 if ( indbattrinv <= 0 ) then
934 call l4f_category_log(qccli%category,L4F_ERROR,"error finding attribute
index in/out
"//qcattrvarsbtables(1))
935 call raise_error("error finding attribute
index in/out
"//qcattrvarsbtables(1))
940 if (present(battrout))then
941 indbattrout = index_c(qccli%v7d%datiattr%b(:)%btable, battrout)
943 indbattrout = index_c(qccli%v7d%datiattr%b(:)%btable, qcattrvarsbtables(2))
946 if ( indbattrout <= 0 ) then
948 call l4f_category_log(qccli%category,L4F_ERROR,"error finding attribute
index in/out
"//qcattrvarsbtables(2))
949 call raise_error("error finding attribute
index in/out
"//qcattrvarsbtables(2))
953 if(present(anamask)) then
958 if(present(timemask)) then
963 if(present(levelmask)) then
964 levelmaskl = levelmask
968 if(present(timerangemask)) then
969 timerangemaskl = timerangemask
971 timerangemaskl = .true.
973 if(present(varmask)) then
978 if(present(networkmask)) then
979 networkmaskl = networkmask
981 networkmaskl = .true.
984 qccli%v7d%voldatiattrb(:,:,:,:,:,:,indbattrout)=ibmiss
985 call init(anavar,"b07030
" )
987 indvar = index(qccli%v7d%anavar, anavar, type=type)
989 do indana=1,size(qccli%v7d%ana)
991 ! if (qccli%height2level) then
992 ! iarea= supermacroa(qccli%in_macroa(indana))
994 iarea= qccli%in_macroa(indana)
996 .not.
if ( c_e(iarea)) cycle
999 ! write(ident,'("box
",2i3.3)')iarea,desc ! macro-area e descrittore
1002 !write(ident,'("box-
",2i2.2)')iarea,lperc ! macro-area e percentile
1003 !call init(ana,lat=lat,lon=lon,ident=ident)
1005 !allocate (indcanav(count(match(qccli%clima%ana(:)%ident,ident))))
1006 !indcanav=match(qccli%clima%ana(:)%ident,ident))))
1012 ! use conventional level starting from station height
1013 if (qccli%height2level) then
1017 ! here we take the height fron any network (the first network win)
1018 do indn=1,size(qccli%v7d%network)
1020 .and.
if( indvar > 0 indn > 0 ) then
1023 height=realdat(qccli%v7d%volanad(indana,indvar,indn),qccli%v7d%anavar%d(indvar))
1025 height=realdat(qccli%v7d%volanar(indana,indvar,indn),qccli%v7d%anavar%r(indvar))
1027 height=realdat(qccli%v7d%volanai(indana,indvar,indn),qccli%v7d%anavar%i(indvar))
1029 height=realdat(qccli%v7d%volanab(indana,indvar,indn),qccli%v7d%anavar%b(indvar))
1031 height=realdat(qccli%v7d%volanac(indana,indvar,indn),qccli%v7d%anavar%c(indvar))
1037 if (c_e(height)) exit
1040 if (c_e(height)) then
1041 .and.
iclv(indana)=firsttrue(cli_level1 <= height height <= cli_level2 )
1047 CALL l4f_log(L4F_DEBUG, 'quaconcli height has value '//t2c(height,"missing
"))
1048 CALL l4f_log(L4F_DEBUG, 'for k having number '//t2c(k)//&
1049 ' iclv has value '//t2c(iclv(indana)))
1054 do indnetwork=1,size(qccli%v7d%network)
1055 do indlevel=1,size(qccli%v7d%level)
1056 do indtimerange=1,size(qccli%v7d%timerange)
1057 do inddativarr=1,size(qccli%v7d%dativar%r)
1058 do indtime=1,size(qccli%v7d%time)
1061 call l4f_log(L4F_DEBUG,"index:
"// t2c(indana)//" "//t2c(indnetwork)//" "//t2c(indlevel)//" "//&
1062 t2c(indtimerange)//" "//t2c(inddativarr)//" "//t2c(indtime))
1064 !!$ print *,"elaboro data :
"//t2c(qccli%v7d%time(indtime))
1066 !!$ forall (indnetwork=1:size(qccli%v7d%network), &
1067 !!$ indlevel=1:size(qccli%v7d%level), &
1068 !!$ indtimerange=1:size(qccli%v7d%timerange), &
1069 !!$ inddativarr=1:size(qccli%v7d%dativar%r), &
1070 !!$ indtime=1:size(qccli%v7d%time))
1072 .and..and..and.
if (anamaskl(indana)timemaskl(indtime)levelmaskl(indlevel) &
1073 .and..and..and.
timerangemaskl(indtimerange)varmaskl(inddativarr)networkmaskl(indnetwork)&
1074 c_e(qccli%v7d%voldatir(indana,indtime,indlevel,indtimerange,inddativarr,indnetwork)))then
1076 if (indbattrinv > 0) then
1077 if( invalidated(qccli%v7d%voldatiattrb&
1078 (indana,indtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrinv))) then
1080 ! invalidated flag allready set
1082 call l4f_log (L4F_DEBUG,"qccli: skip station for a preceding
invalidated flag
")
1088 nintime=qccli%v7d%time(indtime)+timedelta_new(minute=30)
1089 CALL getval(nintime, month=mese, hour=ora)
1091 time=cyclicdatetime_to_conventional(cyclicdatetime_new(month=mese, hour=ora))
1092 !call init(time, year=1001, month=mese, day=1, hour=ora, minute=00)
1094 !!$ print *,"data convenzionale per percentili:
",t2c(time)
1096 level=qccli%v7d%level(indlevel)
1098 call init(network,"qcclima-perc
")
1100 indcnetwork = index(qccli%extreme%network , network)
1101 indctime = index(qccli%extreme%time , time)
1102 indclevel = index(qccli%extreme%level , level)
1103 indctimerange = index(qccli%extreme%timerange , qccli%v7d%timerange(indtimerange))
1105 ! attenzione attenzione TODO
1106 è
! se leggo da bufr il default char e non reale
1107 indcdativarr = index(qccli%extreme%dativar%r, qccli%v7d%dativar%r(inddativarr))
1109 !!$ print *,"dato
",qccli%v7d%timerange(indtimerange)
1110 !!$ print *,"clima
",qccli%clima%timerange
1111 !!$ call l4f_log(L4F_INFO,"index:
"// to_char(indctime)//to_char(indclevel)//&
1112 !!$ to_char(indctimerange)//to_char(indcdativarr)//to_char(indcnetwork))
1114 .or..or..or..or.
!if (indcana <= 0 indctime <= 0 indclevel <= 0 indctimerange <= 0 indcdativarr <= 0 &
1115 .or.
! indcnetwork <= 0 ) cycle
1117 !!$ print *,"vector time
"
1118 !!$ do i=1,size(qccli%extreme%time)
1119 !!$ call display(qccli%extreme%time(i))
1122 !!$ call display(time)
1123 !!$ call display(level)
1124 !!$ call display(qccli%v7d%timerange(indtimerange))
1125 !!$ print *,"indici percentili
",indctime,indclevel,indctimerange,indcdativarr,indcnetwork
1126 .or..or..or.
if (indctime <= 0 indclevel <= 0 indctimerange <= 0 indcdativarr <= 0 &
1127 .or.
indcnetwork <= 0 ) cycle
1129 datoqui = qccli%v7d%voldatir (indana ,indtime ,indlevel ,indtimerange ,inddativarr, indnetwork )
1131 if (c_e(datoqui)) then
1133 ! find extreme in volume
1140 !!$ do i=1, size(qccli%extreme%ana)
1142 !!$ call display(qccli%extreme%ana(i))
1145 if (associated(qccli%extreme%voldatir)) then
1147 if (qccli%height2level) then
1154 !!$ do i =1,size(qccli%extreme%ana)
1155 !!$ call display(qccli%extreme%ana(i))
1159 write(ident,'("#
",i2.2,2i3.3)')k,iarea,desc ! macro-area e descrittore
1160 call init(ana,ident=ident,lat=latc,lon=lonc)
1161 indcana=index(qccli%extreme%ana,ana)
1162 if (indcana > 0 )then
1163 perc25=qccli%extreme%voldatir(indcana,indctime,indclevel,indctimerange,indcdativarr,indcnetwork)
1167 write(ident,'("#
",i2.2,2i3.3)')k,iarea,desc ! macro-area e descrittore
1168 call init(ana,ident=ident,lat=latc,lon=lonc)
1169 indcana=index(qccli%extreme%ana,ana)
1170 !!$ call display(ana)
1171 !!$ print *,"indcana 50
",indcana
1172 if (indcana > 0 )then
1173 perc50=qccli%extreme%voldatir(indcana,indctime,indclevel,indctimerange,indcdativarr,indcnetwork)
1177 write(ident,'("#
",i2.2,2i3.3)')k,iarea,desc ! macro-area e descrittore
1178 call init(ana,ident=ident,lat=latc,lon=lonc)
1179 indcana=index(qccli%extreme%ana,ana)
1180 if (indcana > 0 )then
1181 perc75=qccli%extreme%voldatir(indcana,indctime,indclevel,indctimerange,indcdativarr,indcnetwork)
1185 .not..or..not..or..not.
if ( c_e(perc25) c_e(perc50) c_e(perc75)) cycle
1187 !!$ print *, "datoqui:
",datoqui,"clima ->
",perc25,perc50,perc75
1189 !http://it.wikipedia.org/wiki/Funzione_di_ripartizione_della_variabile_casuale_normale
1190 ! 3.65 for 0.01% each side ( 0.02% total )
1191 extremequii=perc50 - (perc75 - perc25) *1.3 * 3.65 ! 1.3 to go to standard deviation and 3.65 to make 3.65 sigma
1192 extremequif=perc50 + (perc75 - perc25) *1.3 * 3.65 ! 1.3 to go to standard deviation and 3.65 to make 3.65 sigma
1195 call l4f_log (L4F_DEBUG,"qccli: gross error
check "//t2c(extremequii)//">
"//t2c(datoqui)//"<
"//t2c(extremequif))
1199 .or.
if ( datoqui <= extremequii extremequif <= datoqui ) then
1200 ! make gross error check
1202 È
!ATTENZIONE TODO : inddativarr UNA GRANDE SEMPLIFICAZIONE NON VERA SE TIPI DI DATO DIVERSI !!!!
1204 call l4f_log (L4F_DEBUG,"qccli: gross error
check flag set to bad
")
1206 qccli%v7d%voldatiattrb(indana,indtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrout)=qcpar%gross_error
1208 if ( associated ( qccli%data_id_in)) then
1210 call l4f_log (L4F_DEBUG,"id:
"//t2c(&
1211 qccli%data_id_in(indana,indtime,indlevel,indtimerange,indnetwork)))
1213 qccli%data_id_out(indana,indtime,indlevel,indtimerange,indnetwork)=&
1214 qccli%data_id_in(indana,indtime,indlevel,indtimerange,indnetwork)
1218 .not.
else if ( vdge(qccli%v7d%voldatiattrb(indana,indtime,indlevel,indtimerange,&
1219 inddativarr,indnetwork,indbattrout))) then
1221 ! gross error check allready done
1223 call l4f_log (L4F_WARN,"qccli: skip station for a preceding gross error
check flagged bad
")
1227 ! normalize wihout call the subroutine to be more fast and do not change data volume
1229 !print *," -------->
",perc25,perc50,perc75,base_value(qccli%v7d%dativar%r(inddativarr)%btable)
1230 datoqui = (datoqui - perc50) / (perc75 - perc25) + base_value(qccli%v7d%dativar%r(inddativarr)%btable)
1231 !print *,"normalizzato=
",datoqui
1233 ! start to compare with clima dataset (NDI)
1234 call init(network,"qcclima-ndi
")
1235 ! reset the level to standard input data for everyone
1236 level=qccli%v7d%level(indlevel)
1237 time=cyclicdatetime_to_conventional(cyclicdatetime_new(month=mese))
1239 indcnetwork = index(qccli%clima%network , network)
1240 indctime = index(qccli%clima%time , time)
1241 indclevel = index(qccli%clima%level , level)
1242 indctimerange = index(qccli%clima%timerange , qccli%v7d%timerange(indtimerange))
1244 ! attenzione attenzione TODO
1245 è
! se leggo da bufr il default char e non reale
1246 indcdativarr = index(qccli%clima%dativar%r, qccli%v7d%dativar%r(inddativarr))
1249 !!$ print *,"indici clima
", indctime, indclevel, indctimerange, indcdativarr, indcnetwork
1250 .or..or..or.
if (indctime <= 0 indclevel <= 0 indctimerange <= 0 indcdativarr <= 0 &
1251 .or.
indcnetwork <= 0 ) cycle
1255 do desc=1,size(qccli%clima%ana)
1260 write(ident,'("#
",i2.2,2i3.3)')0,0,min(desc,size(qccli%clima%ana)-1) *10 ! macro-area e descrittore
1261 call init(ana,ident=ident,lat=0d0,lon=0d0)
1262 indcana=index(qccli%clima%ana,ana)
1263 if (indcana > 0 )then
1264 climaquif=qccli%clima%voldatir(indcana,indctime,indclevel,indctimerange,indcdativarr,indcnetwork)
1268 write(ident,'("#
",i2.2,2i3.3)')0,0,(desc-1)*10 ! macro-area e descrittore
1269 call init(ana,ident=ident,lat=0d0,lon=0d0)
1270 indcana=index(qccli%clima%ana,ana)
1273 if (indcana > 0 )then
1274 climaquii=qccli%clima%voldatir(indcana,indctime,indclevel,indctimerange,indcdativarr,indcnetwork)
1278 !!$ call l4f_log (L4F_INFO,"ident:
"//qccli%clima%ana(indcana)%ident//ident)
1279 .and..and.
!if ( match(qccli%clima%ana(indcana)%ident,ident) c_e(climaquii) c_e(climaquif)) then
1280 .and.
if ( c_e(climaquii) c_e(climaquif )) then
1282 !!$ print *, "son qua
",trim(qccli%clima%ana(indcana)%ident),trim(ident)
1283 .and.
!!$ where (match(qccli%clima%ana(:)%ident,ident) &
1284 !!$ c_e(qccli%clima%voldatir(indcana,indctime,indclevel,indctimerange,indcdativarr,indcnetwork)))
1285 !!$ call l4f_log (L4F_INFO,"macroarea,iarea,mese,altezza,level
"//&
1286 !!$ trim(to_char(qccli%in_macroa(indana)))//" "//trim(to_char(iarea))&
1287 !!$ //" "//trim(to_char(mese))//" "//trim(to_char(altezza))//" "//trim(to_char(level)))
1290 !!$ print*,"ndi=
",climaquii,datoqui,climaquif
1292 .and..or.
if ( (climaquii <= datoqui datoqui < climaquif) &
1293 .and..or.
(desc == 1 datoqui < climaquii) &
1294 .and.
(desc == size(qccli%clima%ana) datoqui >= climaquif) ) then
1296 if (c_e(qccli%clima%voldatiattrb(indcana &
1297 ,indctime,indclevel,indctimerange,indcdativarr,indcnetwork,1))) then
1299 È
!ATTENZIONE TODO : inddativarr UNA GRANDE SEMPLIFICAZIONE NON VERA SE TIPI DI DATO DIVERSI !!!!
1301 !indcana is the left value as descr-1
1302 qccli%v7d%voldatiattrb(indana,indtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrout)=&
1303 max (qccli%clima%voldatiattrb&
1304 (indcana,indctime,indclevel,indctimerange,indcdativarr,indcnetwork,1)&
1305 , 1_int_b) ! 0 reserved for gross error check
1308 call l4f_log (L4F_DEBUG,"data ndi:
"//t2c(datoqui)//"->
"//&
1309 t2c(qccli%clima%voldatiattrb(indcana,indctime,indclevel,indctimerange,indcdativarr,indcnetwork,1))&
1310 //" :
"//t2c(qccli%v7d%time(indtime)))
1311 call l4f_log (L4F_DEBUG,"limits:
"//t2c(indcana)//":
"//qccli%clima%ana(indcana)%ident//&
1312 " :
"//t2c(climaquii)//" -
"//t2c(climaquif)//" :
"//t2c(qccli%clima%time(indctime)))
1313 call l4f_log (L4F_DEBUG,"qccli: clima
check "//t2c(datoqui)//" confidence:
"//&
1314 t2c(qccli%v7d%voldatiattrb(indana,indtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrout))&
1315 //" :
"//t2c(qccli%v7d%time(indtime)))
1319 if ( associated ( qccli%data_id_in)) then
1321 call l4f_log (L4F_DEBUG,"id:
"//t2c(&
1322 qccli%data_id_in(indana,indtime,indlevel,indtimerange,indnetwork)))
1324 qccli%data_id_out(indana,indtime,indlevel,indtimerange,indnetwork)=&
1325 qccli%data_id_in(indana,indtime,indlevel,indtimerange,indnetwork)
1343 !!$print*,"risultato
"
1344 !!$print *,qccli%v7d%voldatiattrb(:,:,:,:,:,:,indbattrout)
1345 !!$print*,"fine risultato
"
1350 end subroutine quaconcli
1353 !>\brief Return a conventional level for climatological definition.
1354 !! Starting from heigth of station in meter return a level where is defined (I hope)
1355 !! the climatological value
1356 subroutine cli_level(heigth,level)
1358 real,intent(in) :: heigth !< heigth of station in meter
1359 TYPE(vol7d_level),intent(out):: level !< level where is defined the climatological value (layer)
1365 if (c_e(heigth)) then
1366 .and.
i=firsttrue(cli_level1 <= heigth heigth <= cli_level2 )
1369 .and.
if (i >= 1 i <= 10 ) then
1370 call init(level, 102,cli_level1(i)*1000,102,cli_level2(i)*1000)
1372 if (c_e(i)) CALL l4f_log(L4F_DEBUG,"cli_level: strange level, heigth:
"//to_char(heigth))
1376 end subroutine cli_level
1379 !> Initialize level according to climate definition at SIMC
1380 subroutine cli_level_generate(level)
1382 TYPE(vol7d_level),intent(out):: level(:) !< level where is defined the climatological value (layer)
1386 if (size(level) /= cli_nlevel ) then
1387 call l4f_log(L4F_ERROR,"cli_level_generate: level dimension /=
"//trim(to_char(cli_nlevel)))
1388 call raise_error("cli_level_generate: level dimension /=
"//trim(to_char(cli_nlevel)))
1392 call init(level(i), 102,cli_level1(i)*1000,102,cli_level2(i)*1000)
1395 end subroutine cli_level_generate
1398 !!$subroutine qccli_validate(qccli)
1399 !!$type(qcclitype),intent(in) :: qccli
1401 !!$!todo da validare
1404 !!$end subroutine qccli_validate
1407 !> Rielabora le macroarea facendole Valentine/Elements thinking
1409 integer function supermacroa(macroa)
1411 integer, intent(in) :: macroa
1412 ! rielabora le macroarea facendole Valentine thinking
1416 .or..or.
if (macroa == 1 macroa == 2 macroa == 4 ) supermacroa=3
1417 .or..or.
if (macroa == 3 macroa == 5 macroa == 6 ) supermacroa=2
1418 .or.
if (macroa == 7 macroa == 8 ) supermacroa=1
1420 !!$ ! rielabora le macroarea facendole Valentine thinking
1422 .or..or.
!!$ if (qccli%in_macroa(indana) == 1 qccli%in_macroa(indana) == 2 qccli%in_macroa(indana) == 4 ) iarea=3
1423 .or..or.
!!$ if (qccli%in_macroa(indana) == 3 qccli%in_macroa(indana) == 5 qccli%in_macroa(indana) == 6 ) iarea=2
1424 .or.
!!$ if (qccli%in_macroa(indana) == 7 qccli%in_macroa(indana) == 8 ) iarea=1
1426 end function supermacroa
1429 SUBROUTINE qc_compute_percentile(this, perc_vals,cyclicdt,presentperc, presentnumb)
1431 TYPE(qcclitype),INTENT(inout) :: this !< volume providing data to be computed and output data. Input data are not modified by the method, apart from performing a \a vol7d_alloc_vol on it
1432 !TYPE(timedelta),INTENT(in) :: step !< length of the step over which the statistical processing is performed
1433 !TYPE(datetime),INTENT(in),OPTIONAL :: start !< start of statistical processing interval
1434 !TYPE(datetime),INTENT(in),OPTIONAL :: stopp !< end of statistical processing interval
1435 real,intent(in) :: perc_vals(:) !< percentile values to use in compute, between 0. and 100.
1436 TYPE(cyclicdatetime),INTENT(in) :: cyclicdt !< cyclic date and time
1437 real, optional :: presentperc !< percentual of data present for compute (default=not used)
1438 integer, optional :: presentnumb !< number of data present for compute (default=not used)
1439 !!$logical, optional :: height2level !< use conventional level starting from station height
1441 integer :: indana,indtime,indvar,indnetwork,indlevel ,indtimerange ,inddativarr,i,j,k,iana,narea
1443 REAL, DIMENSION(:),allocatable :: perc
1444 TYPE(vol7d_var) :: var
1445 character(len=vol7d_ana_lenident) :: ident
1446 character(len=1) :: type
1447 integer :: areav(size(this%v7d%ana)),iclv(size(this%v7d%ana))
1449 logical,allocatable :: mask(:,:,:),maskplus(:,:,:), maskarea(:)
1450 integer,allocatable :: area(:)
1451 real :: lpresentperc
1452 integer :: lpresentnumb
1454 lpresentperc=optio_r(presentperc)
1455 lpresentnumb=optio_i(presentnumb)
1457 allocate (perc(size(perc_vals)))
1459 call delete(this%extreme)
1460 CALL init(this%extreme, time_definition=this%v7d%time_definition)
1462 call init(var, btable="b01192
") ! MeteoDB station ID that here is the number of area
1465 indvar = index(this%v7d%anavar, var, type=type)
1466 indnetwork=min(1,size(this%v7d%network))
1468 .and.
if( indvar > 0 indnetwork > 0 ) then
1471 areav=integerdat(this%v7d%volanad(:,indvar,indnetwork),this%v7d%anavar%d(indvar))
1473 areav=integerdat(this%v7d%volanar(:,indvar,indnetwork),this%v7d%anavar%r(indvar))
1475 areav=integerdat(this%v7d%volanai(:,indvar,indnetwork),this%v7d%anavar%i(indvar))
1477 areav=integerdat(this%v7d%volanab(:,indvar,indnetwork),this%v7d%anavar%b(indvar))
1479 areav=integerdat(this%v7d%volanac(:,indvar,indnetwork),this%v7d%anavar%c(indvar))
1487 allocate(maskarea(size(this%v7d%ana)))
1488 maskarea(:)= areav(:) /= imiss
1489 narea=count_distinct(areav,maskarea)
1490 allocate(area(narea))
1491 area=pack_distinct(areav,narea,maskarea)
1492 deallocate(maskarea)
1493 if (this%height2level) then
1494 call vol7d_alloc(this%extreme,nana=narea*size(perc_vals)*cli_nlevel)
1496 call vol7d_alloc(this%extreme,nana=narea*size(perc_vals))
1499 if (this%height2level) then
1501 call init(var, btable="b07030
") ! height
1504 indvar = index(this%v7d%anavar, var, type=type)
1507 !!$ CALL l4f_log(L4F_DEBUG, 'SIZE this anavar r '//t2c(SIZE(this%v7d%anavar%r)))
1508 !!$ if (ASSOCIATED(this%anavar%r)) then
1509 !!$ CALL l4f_log(L4F_DEBUG, 'SIZE this anavar r '//t2c(SIZE(this%anavar%r)))
1510 !!$ CALL l4f_log(L4F_DEBUG, 'SIZE this anavar r btable '//t2c(this%anavar%r(SIZE(this%anavar%r))%btable))
1512 !!$ CALL l4f_log(L4F_DEBUG, 'SIZE this anavar i '//t2c(SIZE(this%anavar%i)))
1513 !!$ if (ASSOCIATED(this%anavar%i)) then
1514 !!$ CALL l4f_log(L4F_DEBUG, 'SIZE this anavar i '//t2c(SIZE(this%anavar%i)))
1515 !!$ CALL l4f_log(L4F_DEBUG, 'SIZE this anavar i btable '//t2c(this%anavar%i(SIZE(this%anavar%i))%btable))
1517 !!$ CALL l4f_log(L4F_DEBUG, 'SIZE this anavar d '//t2c(SIZE(this%anavar%d)))
1518 !!$ if (ASSOCIATED(this%anavar%d)) then
1519 !!$ CALL l4f_log(L4F_DEBUG, 'SIZE this anavar d '//t2c(SIZE(this%anavar%d)))
1520 !!$ CALL l4f_log(L4F_DEBUG, 'SIZE this anavar d btable '//t2c(this%anavar%d(SIZE(this%anavar%d))%btable))
1522 !!$ CALL l4f_log(L4F_DEBUG, 'SIZE this anavar b '//t2c(SIZE(this%anavar%b)))
1523 !!$ if (ASSOCIATED(this%anavar%b)) then
1524 !!$ CALL l4f_log(L4F_DEBUG, 'SIZE this anavar b '//t2c(SIZE(this%anavar%b)))
1525 !!$ CALL l4f_log(L4F_DEBUG, 'SIZE this anavar b btable '//t2c(this%anavar%d(SIZE(this%anavar%b))%btable))
1527 !!$ CALL l4f_log(L4F_DEBUG, 'SIZE this anavar c '//t2c(SIZE(this%anavar%c)))
1528 !!$ if (ASSOCIATED(this%anavar%c)) then
1529 !!$ CALL l4f_log(L4F_DEBUG, 'SIZE this anavar c '//t2c(SIZE(this%anavar%c)))
1530 !!$ CALL l4f_log(L4F_DEBUG, 'SIZE this anavar c btable '//t2c(this%anavar%c(SIZE(this%anavar%c))%btable))
1532 !!$ CALL l4f_log(L4F_DEBUG, 'indvar has value '//t2c(indvar))
1533 !!$ CALL l4f_log(L4F_DEBUG, 'indnetwork has value '//t2c(indnetwork))
1536 do k=1,size(this%v7d%ana)
1539 ! here we take the height fron any network (the first network win)
1540 do indnetwork=1,size(this%v7d%network)
1542 .and.
if( indvar > 0 indnetwork > 0 ) then
1545 height=realdat(this%v7d%volanad(k,indvar,indnetwork),this%v7d%anavar%d(indvar))
1547 height=realdat(this%v7d%volanar(k,indvar,indnetwork),this%v7d%anavar%r(indvar))
1549 height=realdat(this%v7d%volanai(k,indvar,indnetwork),this%v7d%anavar%i(indvar))
1551 height=realdat(this%v7d%volanab(k,indvar,indnetwork),this%v7d%anavar%b(indvar))
1553 height=realdat(this%v7d%volanac(k,indvar,indnetwork),this%v7d%anavar%c(indvar))
1559 if (c_e(height)) exit
1562 if (c_e(height)) then
1563 .and.
iclv(k)=firsttrue(cli_level1 <= height height <= cli_level2 )
1569 CALL l4f_log(L4F_DEBUG, 'qc_compute_percentile height has value '//t2c(height))
1570 CALL l4f_log(L4F_DEBUG, 'for k having number '//t2c(k)//&
1571 ' iclv has value '//t2c(iclv(k)))
1579 do j=1,size(perc_vals)
1580 if (this%height2level) then
1582 write(ident,'("#
",i2.2,2i3.3)')k,area(i),nint(perc_vals(j))
1583 call init(this%extreme%ana((k-1)*size(perc_vals)*narea + (j-1)*narea + i),ident=ident,lat=0d0,lon=0d0)
1587 write(ident,'("#
",i2.2,2i3.3)')k,area(i),nint(perc_vals(j))
1588 call init(this%extreme%ana((j-1)*narea+i),ident=ident,lat=0d0,lon=0d0)
1593 !!$do i=1,size(that%ana)
1594 !!$ call display(that%ana(i))
1598 CALL l4f_category_log(this%category, L4F_DEBUG, 'nana has value '//t2c(SIZE(this%v7d%ana)))
1599 CALL l4f_category_log(this%category, L4F_DEBUG, 'lpresentperc has value '//t2c(lpresentperc))
1600 CALL l4f_category_log(this%category, L4F_DEBUG, 'lpresentnumb has value '//t2c(lpresentnumb))
1604 call vol7d_alloc(this%extreme,nlevel=size(this%v7d%level), ntimerange=size(this%v7d%timerange), &
1605 ndativarr=size(this%v7d%dativar%r), nnetwork=1,ntime=1)
1607 this%extreme%level=this%v7d%level
1608 this%extreme%timerange=this%v7d%timerange
1609 this%extreme%dativar%r=this%v7d%dativar%r
1610 this%extreme%time(1)=cyclicdatetime_to_conventional(cyclicdt)
1611 call l4f_category_log(this%category, L4F_DEBUG,"vol7d_compute_percentile conventional
datetime "//to_char(this%extreme%time(1)))
1612 call init(this%extreme%network(1),name="qcclima-perc
")
1614 call vol7d_alloc_vol(this%extreme,inivol=.true.)
1616 allocate (mask(size(this%v7d%ana),size(this%v7d%time),size(this%v7d%network)))
1620 do inddativarr=1,size(this%v7d%dativar%r)
1621 do indtimerange=1,size(this%v7d%timerange)
1622 do indlevel=1,size(this%v7d%level) ! all stations, all times, all networks
1625 !this%v7d%voldatir(indana, indtime, indlevel, indtimerange, inddativarr, indnetwork)
1627 !create mask only with valid time
1628 mask = spread(spread((this%v7d%time == cyclicdt ),1,size(this%v7d%ana)),3,size(this%v7d%network))
1631 CALL l4f_category_log(this%category, L4F_DEBUG, 'count has value '//t2c(count &
1632 .and.
(mask c_e(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:)))))
1635 !delete in mask different area
1636 do j=1, size(mask,1)
1637 if (areav(j) /= area(i)) mask(j,:,:) =.false.
1640 if (this%height2level) then
1641 allocate (maskplus(size(this%v7d%ana),size(this%v7d%time),size(this%v7d%network)))
1644 CALL l4f_category_log(this%category, L4F_DEBUG, 'k has value '//t2c(k))
1647 do iana=1,size(mask,1)
1648 if (iclv(iana) /= k) maskplus(iana,:,:) =.false.
1649 if (iclv(iana) == k) maskplus(iana,:,:) = mask(iana,:,:)
1652 call sub_perc(maskplus)
1655 deallocate(maskplus)
1667 deallocate (perc,mask,area)
1671 subroutine sub_perc(mymask)
1673 logical :: mymask(:,:,:)
1675 ! we want more than 30% data present and a number of data bigger than 100 (default)
1677 .and.
( c_e(lpresentperc) ((float(count &
1678 .and.
(mymask c_e(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:)))&
1680 float(count (mymask))) < lpresentperc)) &
1684 .and.
( c_e(lpresentnumb) (count &
1685 .and.
(mymask c_e(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:))) < lpresentnumb)&
1690 perc= stat_percentile (&
1691 pack(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:), &
1696 do j=1,size(perc_vals)
1697 indana=(k-1)*size(perc_vals)*narea + (j-1)*narea + i
1698 this%extreme%voldatir(indana, indtime, indlevel, indtimerange, inddativarr, indnetwork)=&
1703 end subroutine sub_perc
1706 end SUBROUTINE qc_compute_percentile
1709 SUBROUTINE qc_compute_NormalizedDensityIndex(this, perc_vals,cyclicdt,presentperc, presentnumb,data_normalized)
1711 TYPE(qcclitype),INTENT(inout) :: this !< volume providing data to be computed and output data. Input data are not modified by the method, apart from performing a \a vol7d_alloc_vol on it
1712 !TYPE(timedelta),INTENT(in) :: step !< length of the step over which the statistical processing is performed
1713 !TYPE(datetime),INTENT(in),OPTIONAL :: start !< start of statistical processing interval
1714 !TYPE(datetime),INTENT(in),OPTIONAL :: stopp !< end of statistical processing interval
1715 real,intent(in) :: perc_vals(:) !< percentile values to use in compute, between 0. and 100.
1716 TYPE(cyclicdatetime),INTENT(in) :: cyclicdt !< cyclic date and time
1717 real,optional :: presentperc !< rate of data present for compute on expected values (default=not used)
1718 integer,optional :: presentnumb !< number of data present for compute (default= not used)
1719 logical,optional,intent(in) :: data_normalized !< data input are normalized and the output will be not differentiated by area/height.
1721 integer :: indana,indtime,indvar,indnetwork,indlevel ,indtimerange ,inddativarr, indattr
1722 integer :: i,j,k,iana,narea
1724 TYPE(vol7d_var) :: var
1725 character(len=vol7d_ana_lenident) :: ident
1726 character(len=1) :: type
1727 integer :: areav(size(this%v7d%ana)),iclv(size(this%v7d%ana))
1728 logical,allocatable :: mask(:,:,:),maskplus(:,:,:), maskarea(:)
1729 integer,allocatable :: area(:)
1730 REAL, DIMENSION(:),allocatable :: ndi,limbins
1731 real :: lpresentperc
1732 integer :: lpresentnumb
1735 lnorm = optio_log(data_normalized)
1736 lpresentperc=optio_r(presentperc)
1737 lpresentnumb=optio_i(presentnumb)
1739 allocate (ndi(size(perc_vals)-1),limbins(size(perc_vals)))
1740 call delete(this%clima)
1741 CALL init(this%clima, time_definition=this%v7d%time_definition)
1742 call init(var, btable="b01192
") ! MeteoDB station ID that here is the number of area
1744 .NOT.
if ((lnorm)) then
1747 indvar = index(this%v7d%anavar, var, type=type)
1750 !if( ind /= 0 ) then
1753 areav=integerdat(this%v7d%volanad(:,indvar,indnetwork),this%v7d%anavar%d(indvar))
1755 areav=integerdat(this%v7d%volanar(:,indvar,indnetwork),this%v7d%anavar%r(indvar))
1757 areav=integerdat(this%v7d%volanai(:,indvar,indnetwork),this%v7d%anavar%i(indvar))
1759 areav=integerdat(this%v7d%volanab(:,indvar,indnetwork),this%v7d%anavar%b(indvar))
1761 areav=integerdat(this%v7d%volanac(:,indvar,indnetwork),this%v7d%anavar%c(indvar))
1767 allocate(maskarea(size(this%v7d%ana)))
1768 maskarea(:)= areav(:) /= imiss
1769 narea=count_distinct(areav,maskarea)
1770 allocate(area(narea))
1771 area=pack_distinct(areav,narea,maskarea)
1772 deallocate(maskarea)
1774 if (this%height2level) then
1775 call vol7d_alloc(this%clima,nana=narea*(size(perc_vals)-1)*cli_nlevel)
1777 call vol7d_alloc(this%clima,nana=narea*(size(perc_vals)-1))
1784 do j=1,size(perc_vals)-1
1785 if (this%height2level) then
1787 write(ident,'("#
",i2.2,2i3.3)')k,area(i),nint(perc_vals(j))
1788 call init(this%clima%ana((k-1)*(size(perc_vals)-1)*narea + (j-1)*narea + i),ident=ident,lat=0d0,lon=0d0)
1792 write(ident,'("#
",i2.2,2i3.3)')k,area(i),nint(perc_vals(j))
1793 call init(this%clima%ana((j-1)*narea+i),ident=ident,lat=0d0,lon=0d0)
1802 if (this%height2level) then
1804 call init(var, btable="b07030
") ! height
1807 indvar = index(this%v7d%anavar, var, type=type)
1809 do k=1,size(this%v7d%ana)
1812 ! here we take the height fron any network (the first network win)
1813 do indnetwork=1,size(this%v7d%network)
1815 .and.
if( indvar > 0 indnetwork > 0 ) then
1818 height=realdat(this%v7d%volanad(k,indvar,indnetwork),this%v7d%anavar%d(indvar))
1820 height=realdat(this%v7d%volanar(k,indvar,indnetwork),this%v7d%anavar%r(indvar))
1822 height=realdat(this%v7d%volanai(k,indvar,indnetwork),this%v7d%anavar%i(indvar))
1824 height=realdat(this%v7d%volanab(k,indvar,indnetwork),this%v7d%anavar%b(indvar))
1826 height=realdat(this%v7d%volanac(k,indvar,indnetwork),this%v7d%anavar%c(indvar))
1832 if (c_e(height)) exit
1835 if (c_e(height)) then
1836 .and.
iclv(k)=firsttrue(cli_level1 <= height height <= cli_level2 )
1842 CALL l4f_log(L4F_DEBUG, 'qc_compute_NormalizedDensityIndex height has value '//t2c(height,"missing
"))
1843 CALL l4f_log(L4F_DEBUG, 'for k having number '//t2c(k)//&
1844 ' iclv has value '//t2c(iclv(k)))
1855 call vol7d_alloc(this%clima,nana=size(perc_vals)-1)
1856 do j=1,size(perc_vals)-1
1857 write(ident,'("#
",i2.2,2i3.3)')0,0,nint(perc_vals(j))
1858 call init(this%clima%ana(j),ident=ident,lat=0d0,lon=0d0)
1864 !!$do i=1,size(that%ana)
1865 !!$ call display(that%ana(i))
1868 call vol7d_alloc(this%clima,nlevel=size(this%v7d%level), ntimerange=size(this%v7d%timerange), &
1869 ndativarr=size(this%v7d%dativar%r), nnetwork=1,ntime=1,ndativarattrr=size(this%v7d%dativar%r),ndatiattrr=1)
1871 this%clima%level=this%v7d%level
1872 this%clima%timerange=this%v7d%timerange
1873 this%clima%dativar%r=this%v7d%dativar%r
1874 this%clima%dativarattr%r=this%clima%dativar%r
1875 call init(this%clima%datiattr%r(1), btable="*b33209
") ! NDI order number
1876 this%clima%time(1)=cyclicdatetime_to_conventional(cyclicdt)
1878 call l4f_category_log(this%category,L4F_INFO,"vol7d_compute_ndi conventional
datetime "//to_char(this%clima%time(1)))
1879 call init(this%clima%network(1),name="qcclima-ndi
")
1881 call vol7d_alloc_vol(this%clima,inivol=.true.)
1883 allocate (mask(size(this%v7d%ana),size(this%v7d%time),size(this%v7d%network)))
1888 do inddativarr=1,size(this%v7d%dativar%r)
1889 do indtimerange=1,size(this%v7d%timerange)
1890 do indlevel=1,size(this%v7d%level) ! all stations, all times, all networks
1891 .NOT.
if ((lnorm)) then
1894 !this%v7d%voldatir(indana, indtime, indlevel, indtimerange, inddativarr, indnetwork)
1896 !create mask only with valid time
1897 mask = spread(spread((this%v7d%time == cyclicdt ),1,size(this%v7d%ana)),3,size(this%v7d%network))
1898 !delete in mask different area
1899 do j=1, size(mask,1)
1900 if (areav(j) /= area(i)) mask(j,:,:) =.false.
1903 if (this%height2level) then
1904 allocate (maskplus(size(this%v7d%ana),size(this%v7d%time),size(this%v7d%network)))
1907 CALL l4f_category_log(this%category, L4F_DEBUG, 'k has value '//t2c(k))
1910 do iana=1,size(mask,1)
1911 if (iclv(iana) /= k) maskplus(iana,:,:) =.false.
1912 if (iclv(iana) == k) maskplus(iana,:,:) = mask(iana,:,:)
1915 call sub_ndi(mymask=maskplus)
1918 deallocate(maskplus)
1923 call sub_ndi(mymask=mask)
1931 mask = spread(spread((this%v7d%time == cyclicdt ),1,size(this%v7d%ana)),3,size(this%v7d%network))
1935 call sub_ndi(mymask=mask)
1943 deallocate (ndi,limbins,mask)
1944 .NOT.
if ((lnorm)) deallocate(area)
1948 subroutine sub_ndi(mymask)
1950 logical :: mymask(:,:,:)
1952 ! we want more than 30% data present and a number of data bigger than 100 (default)
1954 .and.
( c_e(lpresentperc) ((float(count &
1955 .and.
(mymask c_e(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:)))&
1957 float(count (mymask))) < lpresentperc)) &
1961 .and.
( c_e(lpresentnumb) (count &
1962 .and.
(mymask c_e(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:))) < lpresentnumb)&
1967 !print*,"-------------------------------------------------------------
"
1969 call NormalizedDensityIndex (&
1970 pack(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:), &
1972 perc_vals, ndi, limbins)
1974 do j=1,size(perc_vals)-1
1975 indana=(k-1)*(size(perc_vals)-1)*narea + (j-1)*narea + i
1976 this%clima%voldatir(indana, indtime, indlevel, indtimerange, inddativarr, indnetwork)=&
1978 this%clima%voldatiattrr(indana, indtime, indlevel, indtimerange, inddativarr, indnetwork,indattr)=&
1982 end subroutine sub_ndi
1985 end SUBROUTINE qc_compute_NormalizedDensityIndex
1990 !> \example v7d_qccli.f90
1991 !! Sample program for module qccli
Methods for returning the value of object members.
Import an array of georef_coord_array objects from a file in ESRI/Shapefile format.
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.
This module defines objects describing georeferenced sparse points possibly with topology and project...
classe per la gestione del logging
Utilities and defines for quality control.
Controllo di qualità climatico.
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
Class for expressing an absolute time value.
Oggetto principale per il controllo di qualitÃ