libsim Versione 7.2.1
modqccli.F90
1! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
2! authors:
3! Davide Cesari <dcesari@arpa.emr.it>
4! Paolo Patruno <ppatruno@arpa.emr.it>
5
6! This program is free software; you can redistribute it and/or
7! modify it under the terms of the GNU General Public License as
8! published by the Free Software Foundation; either version 2 of
9! the License, or (at your option) any later version.
10
11! This program is distributed in the hope that it will be useful,
12! but WITHOUT ANY WARRANTY; without even the implied warranty of
13! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14! GNU General Public License for more details.
15
16! You should have received a copy of the GNU General Public License
17! along with this program. If not, see <http://www.gnu.org/licenses/>.
18
19#include "config.h"
20
93
100
101!! README
102!! di fatto sono necessari due file v7d o due dsn : uno per il clima e uno per gli estremi del gross error check.
103!! Nel volume extreme ci devono essere due stazioni (differente ident) uno per il valore minimo e uno per il valore massimo.
104!! esiste uno switch hright2level per attivare le altezze convenzionali
105!! bisogna gestire bene le aree e non supporre che ci sia e sia =1
106
107module modqccli
108
110use vol7d_class
111use modqc
113use log4fortran
116use simple_stat
117!use array_utilities
118!use io_units
119#ifdef HAVE_DBALLE
121#endif
122
123implicit none
124
125character (len=255),parameter:: subcategory="QCcli"
126
127integer, parameter :: cli_nlevel=10
129integer, parameter :: cli_level1(cli_nlevel) = (/-100,100,250,500,750,1000,1250,1500,1750,2000/)
131integer, parameter :: cli_level2(cli_nlevel) = (/100,250,500,750,1000,1250,1500,1750,2000,2250/)
132
134type :: qcclitype
135
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(:)
142 integer :: category
143 logical :: height2level
144
145end type qcclitype
146
147
149interface init
150 module procedure qccliinit
151end interface
152
154interface alloc
155 module procedure qcclialloc
156end interface
157
159interface delete
160 module procedure qcclidelete
161end interface
162
163PRIVATE
164PUBLIC cli_nlevel, cli_level1, cli_level2, qcclitype, init, alloc, delete, &
165 vol7d_normalize_data, quaconcli, cli_level, cli_level_generate, &
166 qc_compute_percentile, qc_compute_normalizeddensityindex
167
168contains
169
172
173subroutine qccliinit(qccli,v7d,var, timei, timef, data_id_in,&
174 macropath, climapath, extremepath, &
175#ifdef HAVE_DBALLE
176 dsncli,dsnextreme,user,password,&
177#endif
178 height2level,categoryappend)
179
180type(qcclitype),intent(in out) :: qccli
181type (vol7d),intent(in),target:: v7d
182character(len=*),INTENT(in) :: var(:)
184TYPE(datetime),INTENT(in),optional :: timei, timef
185integer,intent(in),optional,target:: data_id_in(:,:,:,:,:)
186character(len=*),intent(in),optional :: macropath
187character(len=*),intent(in),optional :: climapath
188character(len=*),intent(in),optional :: extremepath
189logical ,intent(in),optional :: height2level
190character(len=*),INTENT(in),OPTIONAL :: categoryappend
194
195#ifdef HAVE_DBALLE
196type (vol7d_dballe) :: v7d_dballecli
197character(len=*),intent(in),optional :: dsncli
198character(len=*),intent(in),optional :: user
199character(len=*),intent(in),optional :: password
200character(len=512) :: ldsncli
201character(len=512) :: luser
202character(len=512) :: lpassword
203type (vol7d_dballe) :: v7d_dballeextreme
204character(len=*),intent(in),optional :: dsnextreme
205character(len=512) :: ldsnextreme
206TYPE(datetime) :: ltimei, ltimef
207integer :: yeari, yearf, monthi, monthf, dayi, dayf,&
208 houri, minutei, mseci, hourf, minutef, msecf
209#endif
210
211integer :: iuni,i,j
212character(len=512) :: filepath
213character(len=512) :: filepathclima
214character(len=512) :: filepathextreme
215character(len=512) :: a_name
216TYPE(geo_coord) :: lcoordmin,lcoordmax
217TYPE(vol7d_network) :: network
218
219! temporary, improve!!!!
220DOUBLE PRECISION :: lon, lat
221TYPE(georef_coord) :: point
222TYPE(arrayof_georef_coord_array) :: macroa
223
224
225call l4f_launcher(a_name,a_name_append=trim(subcategory)//"."//trim(categoryappend))
226qccli%category=l4f_category_get(a_name)
227
228qccli%height2level=optio_log(height2level)
229
230! all coordinates in clima and extreme datasets are zero
231call init(lcoordmin)
232call init(lcoordmax)
233
234!!$! mmm... I am not sure here .... this will be removed
235!!$if (qccli%height2level) then
236!!$ if (present(coordmin)) lcoordmin=coordmin
237!!$ if (present(coordmax)) lcoordmax=coordmax
238!!$end if
239
240nullify ( qccli%in_macroa )
241nullify ( qccli%data_id_in )
242nullify ( qccli%data_id_out )
243
244! riporto il volume dati nel mio oggetto
245qccli%v7d => v7d
246
247if (present(data_id_in))then
248 qccli%data_id_in => data_id_in
249end if
250
251if (qccli%height2level) then
252 !shape file for Emilia Romagna clima !
253 filepath=get_package_filepath('share:sup_macroaree.shp', filetype_data)
254else
255 !shape file for Europa clima !
256 filepath=get_package_filepath('share:ens_v8_ll.shp', filetype_data)
257end if
258
259if (present(macropath))then
260 if (c_e(macropath)) then
261 filepath=macropath
262 end if
263end if
264
265CALL import(macroa, shpfile=filepath)
266
267call optio(climapath,filepathclima)
268call optio(extremepath,filepathextreme)
269
270#ifdef HAVE_DBALLE
271
272call init(network,"qcclima-ndi")
273
274ltimei=datetime_miss
275ltimef=datetime_miss
276if (present (timei)) then
277 ltimei=timei
278 ltimei=ltimei+timedelta_new(minute=30)
279end if
280
281if (present (timef)) then
282 ltimef=timef
283 ltimef=ltimef+timedelta_new(minute=30)
284end if
285
286call getval(ltimei, year=yeari, month=monthi, day=dayi, hour=houri, minute=minutei, msec=mseci)
287call getval(ltimef, year=yearf, month=monthf, day=dayf, hour=hourf, minute=minutef, msec=msecf)
288
289if ( c_e(yeari) .and. c_e(yearf) .and. yeari == yearf .and. monthi == monthf ) then
290
291! call init(ltimei, 1001, monthi, 1, houri, minutei, mseci)
292! call init(ltimef, 1001, monthf, 1, hourf, minutef, msecf)
293 ltimei=cyclicdatetime_to_conventional(cyclicdatetime_new(month=monthi))
294 ltimef=cyclicdatetime_to_conventional(cyclicdatetime_new(month=monthf))
296else
297 ! if you span years or months I read all the climat dataset (should be optimized not so easy)
298 ltimei=datetime_miss
299 ltimef=datetime_miss
300
301end if
302
303call optio(dsncli,ldsncli)
304call optio(user,luser)
305call optio(password,lpassword)
306
307if ((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")
309 call raise_error()
310end if
311
312if (.not. c_e(ldsncli)) then
313
314#endif
315
316 if (.not. c_e(filepathclima)) then
317 filepathclima=get_package_filepath('qcclima-ndi.v7d', filetype_data)
318 end if
320 if (c_e(filepathclima))then
321
322 select case (trim(lowercase(suffixname(filepathclima))))
323
324 case("v7d")
325 iuni=getunit()
326 call import(qccli%clima,filename=filepathclima,unit=iuni)
327 close (unit=iuni)
329#ifdef HAVE_DBALLE
330 case("bufr")
331 call init(v7d_dballecli,file=.true.,filename=filepathclima,categoryappend=trim(a_name)//".clima")
332 !call import(v7d_dballecli)
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)
336 call delete(v7d_dballecli)
337#endif
338
339 case default
340 call l4f_category_log(qccli%category,l4f_error,&
341 "file type not supported (use .v7d or .bufr suffix only): "//trim(filepathclima))
342 call raise_error()
343 end select
344
345 else
346 call l4f_category_log(qccli%category,l4f_warn,"clima volume not iniziatized: clima QC will not be possible")
347 call init(qccli%clima)
348! call raise_fatal_error()
349 end if
350
351#ifdef HAVE_DBALLE
352else
353
354 call l4f_category_log(qccli%category,l4f_debug,"init v7d_dballecli")
355 call init(v7d_dballecli,dsn=ldsncli,user=luser,password=lpassword,write=.false.,&
356 file=.false.,categoryappend=trim(a_name)//".clima")
357 call l4f_category_log(qccli%category,l4f_debug,"import v7d_dballecli")
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)
362
363end if
364#endif
365
366#ifdef HAVE_DBALLE
367call delete(ltimei)
368call delete(ltimef)
369
370if ( c_e(yeari) .and. c_e(yearf) .and. yeari == yearf .and. monthi == monthf ) then
371
372 if ( dayi == dayf .and. houri == hourf .and. minutei == minutef .and. mseci == msecf ) then
373
374 ltimei=cyclicdatetime_to_conventional(cyclicdatetime_new(month=monthi, hour=houri))
375 ltimef=cyclicdatetime_to_conventional(cyclicdatetime_new(month=monthf, hour=hourf))
376
377 else
378
379 ltimei=cyclicdatetime_to_conventional(cyclicdatetime_new(month=monthi, hour=00))
380 ltimef=cyclicdatetime_to_conventional(cyclicdatetime_new(month=monthf, hour=23))
381
382 end if
383
384else
385 ! if you span years or months or days I read all the climat dataset (should be optimized not so easy)
386 ltimei=datetime_miss
387 ltimef=datetime_miss
388
389end if
390
391call init(network,"qcclima-perc")
392call optio(dsnextreme,ldsnextreme)
393
394if (.not. c_e(ldsnextreme)) then
395
396#endif
397
398 if (.not. c_e(filepathextreme)) then
399 filepathextreme=get_package_filepath('qcclima-extreme.v7d', filetype_data)
400 end if
401
402 if (c_e(filepathextreme)) then
403
404 select case (trim(lowercase(suffixname(filepathextreme))))
405
406 case("v7d")
407 iuni=getunit()
408 call import(qccli%extreme,filename=filepathextreme,unit=iuni)
409 close (unit=iuni)
410
411#ifdef HAVE_DBALLE
412 case("bufr")
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)
419#endif
420
421 case default
422
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))
426 call raise_error()
427 end if
428 end select
429
430 else
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)
434 end if
435
436
437#ifdef HAVE_DBALLE
438else
439
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")
444
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)
449
450end if
451
452call delete(ltimei)
453call delete(ltimef)
454#endif
455
456
457call qcclialloc(qccli)
458
459
460! valuto in quale macroarea sono le stazioni
461
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()
465!!$ENDIF
466
467if (associated(qccli%in_macroa)) then
468 qccli%in_macroa = imiss
469
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
477 EXIT
478 ENDIF
479 ENDDO
480 ENDDO
481end if
482
483call delete(macroa)
484
485return
486end subroutine qccliinit
487
488
490subroutine qcclialloc(qccli)
491 ! pseudo costruttore con distruttore automatico
492
493type(qcclitype),intent(in out) :: qccli
494
495integer :: istatt
496integer :: sh(5)
497
498! se ti sei dimenticato di deallocare ci penso io
499call qcclidealloc(qccli)
500
501
502!!$if (associated (qccli%v7d%dativar%r )) then
503!!$ nv=size(qccli%v7d%dativar%r)
504!!$
505!!$ allocate(qccli%valminr(nv),stat=istat)
506!!$ istatt=istatt+istat
507!!$ allocate(qccli%valmaxr(nv),stat=istat)
508!!$ istatt=istatt+istat
509!!$
510!!$ if (istatt /= 0) ier=1
511!!$
512!!$end if
513
514if (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")
519 end if
520end if
521
522if (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)
525 if (istatt /= 0)then
526 call l4f_category_log(qccli%category,l4f_error,"allocate error")
527 call raise_error("allocate error")
528 else
529 qccli%data_id_out=imiss
530 end if
531end if
532
533return
534
535end subroutine qcclialloc
536
537
539
540subroutine qcclidealloc(qccli)
541 ! pseudo distruttore
542
543type(qcclitype),intent(in out) :: qccli
544
545!!$if ( associated ( qccli%valminr)) then
546!!$ deallocate(qccli%valminr)
547!!$end if
548!!$
549!!$if ( associated ( qccli%valmaxr)) then
550!!$ deallocate(qccli%valmaxr)
551!!$end if
552
553if (associated (qccli%in_macroa)) then
554 deallocate (qccli%in_macroa)
555end if
556
557if (associated(qccli%data_id_out))then
558 deallocate (qccli%data_id_out)
559end if
560
561return
562end subroutine qcclidealloc
563
564
566
567
568subroutine qcclidelete(qccli)
569 ! decostruttore a mezzo
570type(qcclitype),intent(in out) :: qccli
571
572call qcclidealloc(qccli)
573
574call delete(qccli%clima)
575call delete(qccli%extreme)
576
577!delete logger
578call l4f_category_delete(qccli%category)
579
580return
581end subroutine qcclidelete
582
583
584
597SUBROUTINE vol7d_normalize_data(qccli)
598
599TYPE(qcclitype),INTENT(inout) :: qccli
600!character (len=10) ,intent(in),optional :: battrinv !< attributo invalidated in input/output
601
602real :: datoqui, perc25, perc50,perc75
603integer :: indana , indtime ,indlevel ,indtimerange ,inddativarr, indnetwork
604integer :: indcana, indvar,indctime,indclevel,indctimerange,indcdativarr,indcnetwork
605!integer :: indbattrinv
606integer :: iclv(size(qccli%v7d%ana))
607real :: height
608character(len=1) :: type
609TYPE(vol7d_var) :: var
610TYPE(vol7d_ana) :: ana
611TYPE(datetime) :: time, nintime
612TYPE(vol7d_level):: level
613integer :: mese, ora, desc, iarea, k
614
615integer :: clev
616character(len=1) :: mycanc, canc = "#"
617
618CHARACTER(len=vol7d_ana_lenident) :: ident
619
620
621!!$indbattrinv=0
622!!$if (associated(qccli%v7d%dativarattr%b))then
623!!$ if (present(battrinv))then
624!!$ indbattrinv = index_c(qccli%v7d%dativarattr%b(:)%btable, battrinv)
625!!$ else
626!!$ indbattrinv = index_c(qccli%v7d%dativarattr%b(:)%btable, qcattrvarsbtables(1))
627!!$ end if
628!!$end if
629
630
631if (.not. 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()
635 return
636end if
637
638
639if (qccli%height2level) then
640 call init(var, btable="B07030") ! height
641
642 type=cmiss
643 indvar = index(qccli%v7d%anavar, var, type=type)
644
645 do indana=1,size(qccli%v7d%ana)
646 height=rmiss
647
648 ! here we take the height fron any network (the first network win)
649 do indnetwork=1,size(qccli%v7d%network)
650
651 if( indvar > 0 ) then
652 select case (type)
653 case("d")
654 height=realdat(qccli%v7d%volanad(indana,indvar,indnetwork),qccli%v7d%anavar%d(indvar))
655 case("r")
656 height=realdat(qccli%v7d%volanar(indana,indvar,indnetwork),qccli%v7d%anavar%r(indvar))
657 case ("i")
658 height=realdat(qccli%v7d%volanai(indana,indvar,indnetwork),qccli%v7d%anavar%i(indvar))
659 case("b")
660 height=realdat(qccli%v7d%volanab(indana,indvar,indnetwork),qccli%v7d%anavar%b(indvar))
661 case("c")
662 height=realdat(qccli%v7d%volanac(indana,indvar,indnetwork),qccli%v7d%anavar%c(indvar))
663 end select
664 end if
665
666 if (c_e(height)) exit
667 end do
668
669 if (c_e(height)) then
670 iclv(indana)=firsttrue(cli_level1 <= height .and. height <= cli_level2 )
671 else
672 iclv(indana)=imiss
673 endif
674
675#ifdef DEBUG
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"))
679#endif
680 end do
681
682endif
683
684
685if (.not. 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()
689 return
690end if
691
692if (.not. 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()
696 return
697end if
698
699do indana=1,size(qccli%v7d%ana)
700 iarea= qccli%in_macroa(indana)
701
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)
707
708 datoqui = qccli%v7d%voldatir (indana ,indtime ,indlevel ,indtimerange ,inddativarr, indnetwork )
709
710 if (.not. c_e(datoqui)) cycle
711
712 if (.not. c_e(iarea)) then
713 qccli%v7d%voldatir (indana ,indtime ,indlevel ,indtimerange ,&
714 inddativarr, indnetwork ) = rmiss
715 cycle
716 end if
717
718!!$ if (indbattrinv > 0) then
719!!$ if( invalidated(qccli%v7d%voldatiattrb&
720!!$ (indana,indtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrinv))) cycle
721!!$ end if
722
723 nintime=qccli%v7d%time(indtime)+timedelta_new(minute=30)
724 CALL getval(nintime, month=mese, hour=ora)
725
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)
730
731 indcnetwork = 1
732
733 !indcana = firsttrue(qccli%extreme%ana == ana)
734
735 indctime = index(qccli%extreme%time , time)
736 indclevel = index(qccli%extreme%level , level)
737 indctimerange = index(qccli%extreme%timerange , qccli%v7d%timerange(indtimerange))
738
739 ! attenzione attenzione TODO
740 ! se leggo da bufr il default è char e non reale
741
742 indcdativarr = index(qccli%extreme%dativar%r, qccli%v7d%dativar%r(inddativarr))
743
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))
748
749 !if (indcana <= 0 .or. indctime <= 0 .or. indclevel <= 0 .or. indctimerange <= 0 .or. indcdativarr <= 0 &
750 ! .or. indcnetwork <= 0 ) cycle
751 if (indctime <= 0 .or. indclevel <= 0 .or. indctimerange <= 0 .or. indcdativarr <= 0 &
752 .or. indcnetwork <= 0 ) then
753 qccli%v7d%voldatir (indana ,indtime ,indlevel ,indtimerange ,&
754 inddativarr, indnetwork ) = rmiss
755 cycle
756 end if
757
758 ! find percentile in volume
759
760 perc25=rmiss
761 perc50=rmiss
762 perc75=rmiss
763
764
765 if (qccli%height2level) then
766 k=iclv(indana)
767 else
768 k=0
769 end if
770
771 if (.not. c_e(k)) then
772 qccli%v7d%voldatir (indana ,indtime ,indlevel ,indtimerange ,&
773 inddativarr, indnetwork ) = rmiss
774 cycle
775 end if
776
777 desc=25
778 write(ident,'("#",i2.2,2i3.3)')k,iarea,desc ! macro-area e descrittore
779
780!!$ if (indana == 1996 .and. indtime == 855 .and. indlevel == 1 .and. indtimerange==1 &
781!!$ .and. inddativarr==1 .and. indnetwork==1 ) then
782!!$ call l4f_category_log(qccli%category,L4F_WARN,"mmmmmmmmmmmmm: "//t2c(datoqui)//ident)
783!!$ end if
784
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)
793 end if
794
795 desc=50
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)
805 end if
806
807 desc=75
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)
817 end if
818
819 if ( c_e(perc25) .and. c_e(perc50) .and. c_e(perc75) .and. &
820 abs( perc75 - perc25 ) >= spacing( max(abs(perc75),abs(perc25))))then
821 ! normalize
822
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))
827
828 else
829
830 datoqui=rmiss
831
832 endif
833
834 qccli%v7d%voldatir (indana ,indtime ,indlevel ,indtimerange ,&
835 inddativarr, indnetwork ) = datoqui
836
837 end do
838 end do
839 end do
840 end do
841 end do
842
843 read (qccli%v7d%ana(indana)%ident,'(a1,i2.2,2i3.3)') mycanc, clev, iarea, desc
844 if (mycanc == canc) then
845 clev=0
846 iarea=0
847 write (qccli%v7d%ana(indana)%ident,'(a1,i2.2,2i3.3)') canc, clev, iarea, desc
848 end if
849
850end do
851
852end SUBROUTINE vol7d_normalize_data
853
854
855real function base_value(btable)
856character (len=10) ,intent(in):: btable
857
858character (len=10) :: btables(1) =(/"B12101"/)
859real :: base_values(1) =(/273.15 /)
860integer :: ind
861
862ind = index_c(btables,btable)
863
864if (ind > 0) then
865 base_value = base_values(ind)
866else
867 call l4f_log(l4f_warn,"modqccli_base_value: variable "//btable//" do not have base value")
868 base_value = 0.
869end if
870
871return
872
873end function base_value
874
875
882SUBROUTINE quaconcli (qccli,battrinv,battrout,&
883 anamask,timemask,levelmask,timerangemask,varmask,networkmask)
884
885
886type(qcclitype),intent(in out) :: qccli
887character (len=10) ,intent(in),optional :: battrinv
888character (len=10) ,intent(in),optional :: battrout
889logical ,intent(in),optional :: anamask(:)
890logical ,intent(in),optional :: timemask(:)
891logical ,intent(in),optional :: levelmask(:)
892logical ,intent(in),optional :: timerangemask(:)
893logical ,intent(in),optional :: varmask(:)
894logical ,intent(in),optional :: networkmask(:)
895
896CHARACTER(len=vol7d_ana_lenident) :: ident
897REAL(kind=fp_geo) :: latc,lonc
898integer :: mese, ora
899 !local
900integer :: indbattrinv,indbattrout
901logical :: 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))
903
904integer :: indana , indtime ,indlevel ,indtimerange ,inddativarr, indnetwork
905integer :: indcana, indctime,indclevel,indctimerange,indcdativarr,indcnetwork
906real :: datoqui,climaquii,climaquif, extremequii,extremequif,perc25,perc50,perc75
907integer :: iarea,desc,indn,indvar,k
908real :: height
909!integer, allocatable :: indcanav(:)
910
911
912TYPE(vol7d_network) :: network
913TYPE(vol7d_ana) :: ana
914TYPE(datetime) :: time, nintime
915TYPE(vol7d_level):: level
916type(vol7d_var) :: anavar
917integer :: iclv(size(qccli%v7d%ana))
918character(len=1) :: type
919
920
921!call qccli_validate (qccli)
922
923indbattrinv=0
924if (associated(qccli%v7d%datiattr%b))then
925 if (present(battrinv))then
926 indbattrinv = index_c(qccli%v7d%datiattr%b(:)%btable, battrinv)
927 else
928 indbattrinv = index_c(qccli%v7d%datiattr%b(:)%btable, qcattrvarsbtables(1))
929 end if
930end if
931
932if ( indbattrinv <= 0 ) then
933
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))
936
937end if
938
939
940if (present(battrout))then
941 indbattrout = index_c(qccli%v7d%datiattr%b(:)%btable, battrout)
942else
943 indbattrout = index_c(qccli%v7d%datiattr%b(:)%btable, qcattrvarsbtables(2))
944end if
945
946if ( indbattrout <= 0 ) then
947
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))
950
951end if
952
953if(present(anamask)) then
954 anamaskl = anamask
955else
956 anamaskl = .true.
957endif
958if(present(timemask)) then
959 timemaskl = timemask
960else
961 timemaskl = .true.
962endif
963if(present(levelmask)) then
964 levelmaskl = levelmask
965else
966 levelmaskl = .true.
967endif
968if(present(timerangemask)) then
969 timerangemaskl = timerangemask
970else
971 timerangemaskl = .true.
972endif
973if(present(varmask)) then
974 varmaskl = varmask
975else
976 varmaskl = .true.
977endif
978if(present(networkmask)) then
979 networkmaskl = networkmask
980else
981 networkmaskl = .true.
982endif
983
984qccli%v7d%voldatiattrb(:,:,:,:,:,:,indbattrout)=ibmiss
985call init(anavar,"B07030" )
986type=cmiss
987indvar = index(qccli%v7d%anavar, anavar, type=type)
988
989do indana=1,size(qccli%v7d%ana)
990
991! if (qccli%height2level) then
992! iarea= supermacroa(qccli%in_macroa(indana))
993! else
994 iarea= qccli%in_macroa(indana)
995
996 if (.not. c_e(iarea)) cycle
997
998! end if
999! write(ident,'("BOX",2i3.3)')iarea,desc ! macro-area e descrittore
1000 !lat=0.0d0
1001 !lon=0.0d0
1002 !write(ident,'("BOX-",2i2.2)')iarea,lperc ! macro-area e percentile
1003 !call init(ana,lat=lat,lon=lon,ident=ident)
1004
1005 !allocate (indcanav(count(match(qccli%clima%ana(:)%ident,ident))))
1006 !indcanav=match(qccli%clima%ana(:)%ident,ident))))
1007
1008 latc=0.d0
1009 lonc=0.d0
1010
1011
1012 ! use conventional level starting from station height
1013 if (qccli%height2level) then
1014
1015 height=rmiss
1016
1017 ! here we take the height fron any network (the first network win)
1018 do indn=1,size(qccli%v7d%network)
1019
1020 if( indvar > 0 .and. indn > 0 ) then
1021 select case (type)
1022 case("d")
1023 height=realdat(qccli%v7d%volanad(indana,indvar,indn),qccli%v7d%anavar%d(indvar))
1024 case("r")
1025 height=realdat(qccli%v7d%volanar(indana,indvar,indn),qccli%v7d%anavar%r(indvar))
1026 case ("i")
1027 height=realdat(qccli%v7d%volanai(indana,indvar,indn),qccli%v7d%anavar%i(indvar))
1028 case("b")
1029 height=realdat(qccli%v7d%volanab(indana,indvar,indn),qccli%v7d%anavar%b(indvar))
1030 case("c")
1031 height=realdat(qccli%v7d%volanac(indana,indvar,indn),qccli%v7d%anavar%c(indvar))
1032 case default
1033 height=rmiss
1034 end select
1035 end if
1036
1037 if (c_e(height)) exit
1038 end do
1039
1040 if (c_e(height)) then
1041 iclv(indana)=firsttrue(cli_level1 <= height .and. height <= cli_level2 )
1042 else
1043 iclv(indana)=imiss
1044 endif
1045
1046#ifdef DEBUG
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)))
1050#endif
1051
1052 end if
1053
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)
1059
1060#ifdef DEBUG
1061 call l4f_log(l4f_debug,"Index:"// t2c(indana)//" "//t2c(indnetwork)//" "//t2c(indlevel)//" "//&
1062 t2c(indtimerange)//" "//t2c(inddativarr)//" "//t2c(indtime))
1063#endif
1064!!$ print *,"elaboro data : "//t2c(qccli%v7d%time(indtime))
1065!!$
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))
1071
1072 if (anamaskl(indana).and.timemaskl(indtime).and.levelmaskl(indlevel).and. &
1073 timerangemaskl(indtimerange).and.varmaskl(inddativarr).and.networkmaskl(indnetwork).and.&
1074 c_e(qccli%v7d%voldatir(indana,indtime,indlevel,indtimerange,inddativarr,indnetwork)))then
1075
1076 if (indbattrinv > 0) then
1077 if( invalidated(qccli%v7d%voldatiattrb&
1078 (indana,indtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrinv))) then
1079
1080 ! invalidated flag allready set
1081#ifdef DEBUG
1082 call l4f_log (l4f_debug,"qccli: skip station for a preceding invalidated flag")
1083#endif
1084 cycle
1085 end if
1086 end if
1087
1088 nintime=qccli%v7d%time(indtime)+timedelta_new(minute=30)
1089 CALL getval(nintime, month=mese, hour=ora)
1090
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)
1093
1094!!$ print *,"data convenzionale per percentili: ",t2c(time)
1095
1096 level=qccli%v7d%level(indlevel)
1097
1098 call init(network,"qcclima-perc")
1099
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))
1104
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))
1108
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))
1113
1114 !if (indcana <= 0 .or. indctime <= 0 .or. indclevel <= 0 .or. indctimerange <= 0 .or. indcdativarr <= 0 &
1115 ! .or. indcnetwork <= 0 ) cycle
1116
1117!!$ print *,"vector time"
1118!!$ do i=1,size(qccli%extreme%time)
1119!!$ call display(qccli%extreme%time(i))
1120!!$ end do
1121!!$ print *,"time"
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 if (indctime <= 0 .or. indclevel <= 0 .or. indctimerange <= 0 .or. indcdativarr <= 0 &
1127 .or. indcnetwork <= 0 ) cycle
1128
1129 datoqui = qccli%v7d%voldatir (indana ,indtime ,indlevel ,indtimerange ,inddativarr, indnetwork )
1130
1131 if (c_e(datoqui)) then
1132
1133 ! find extreme in volume
1134 extremequii=rmiss
1135 extremequif=rmiss
1136 perc25=rmiss
1137 perc50=rmiss
1138 perc75=rmiss
1139
1140!!$ do i=1, size(qccli%extreme%ana)
1141!!$ print *,i
1142!!$ call display(qccli%extreme%ana(i))
1143!!$ end do
1144
1145 if (associated(qccli%extreme%voldatir)) then
1146
1147 if (qccli%height2level) then
1148 k=iclv(indana)
1149 else
1150 k=0
1151 end if
1152
1153
1154!!$ do i =1,size(qccli%extreme%ana)
1155!!$ call display(qccli%extreme%ana(i))
1156!!$ end do
1157
1158 desc=25 ! minimum
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)
1164 end if
1165
1166 desc=50 ! mediana
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)
1174 end if
1175
1176 desc=75 ! maximum
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)
1182 end if
1183 end if
1184
1185 if ( .not. c_e(perc25) .or. .not. c_e(perc50) .or. .not. c_e(perc75)) cycle
1186
1187!!$ print *, "datoqui: ",datoqui,"clima ->",perc25,perc50,perc75
1188
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
1193
1194#ifdef DEBUG
1195 call l4f_log (l4f_debug,"qccli: gross error check "//t2c(extremequii)//">"//t2c(datoqui)//"<"//t2c(extremequif))
1196#endif
1197
1198
1199 if ( datoqui <= extremequii .or. extremequif <= datoqui ) then
1200 ! make gross error check
1201
1202 !ATTENZIONE TODO : inddativarr È UNA GRANDE SEMPLIFICAZIONE NON VERA SE TIPI DI DATO DIVERSI !!!!
1203#ifdef DEBUG
1204 call l4f_log (l4f_debug,"qccli: gross error check flag set to bad")
1205#endif
1206 qccli%v7d%voldatiattrb(indana,indtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrout)=qcpar%gross_error
1207
1208 if ( associated ( qccli%data_id_in)) then
1209#ifdef DEBUG
1210 call l4f_log (l4f_debug,"id: "//t2c(&
1211 qccli%data_id_in(indana,indtime,indlevel,indtimerange,indnetwork)))
1212#endif
1213 qccli%data_id_out(indana,indtime,indlevel,indtimerange,indnetwork)=&
1214 qccli%data_id_in(indana,indtime,indlevel,indtimerange,indnetwork)
1215 end if
1216
1217
1218 else if (.not. vdge(qccli%v7d%voldatiattrb(indana,indtime,indlevel,indtimerange,&
1219 inddativarr,indnetwork,indbattrout))) then
1220
1221 ! gross error check allready done
1222#ifdef DEBUG
1223 call l4f_log (l4f_warn,"qccli: skip station for a preceding gross error check flagged bad")
1224#endif
1225 else
1226
1227 ! normalize wihout call the subroutine to be more fast and do not change data volume
1228
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
1232
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))
1238
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))
1243
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))
1247
1248
1249!!$ print *,"indici clima", indctime, indclevel, indctimerange, indcdativarr, indcnetwork
1250 if (indctime <= 0 .or. indclevel <= 0 .or. indctimerange <= 0 .or. indcdativarr <= 0 &
1251 .or. indcnetwork <= 0 ) cycle
1252
1253 !climat check
1254
1255 do desc=1,size(qccli%clima%ana)
1256
1257 climaquii=rmiss
1258 climaquif=rmiss
1259
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)
1265 end if
1266
1267
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)
1271! call display(ana)
1272! print*,indcana
1273 if (indcana > 0 )then
1274 climaquii=qccli%clima%voldatir(indcana,indctime,indclevel,indctimerange,indcdativarr,indcnetwork)
1275 end if
1276
1277
1278!!$ call l4f_log (L4F_INFO,"ident: "//qccli%clima%ana(indcana)%ident//ident)
1279 !if ( match(qccli%clima%ana(indcana)%ident,ident) .and. c_e(climaquii) .and. c_e(climaquif)) then
1280 if ( c_e(climaquii) .and. c_e(climaquif )) then
1281
1282!!$ print *, "son qua",trim(qccli%clima%ana(indcana)%ident),trim(ident)
1283!!$ where (match(qccli%clima%ana(:)%ident,ident).and. &
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)))
1288!!$
1289
1290!!$ print*,"ndi=",climaquii,datoqui,climaquif
1291
1292 if ( (climaquii <= datoqui.and. datoqui < climaquif) .or. &
1293 (desc == 1 .and. datoqui < climaquii) .or. &
1294 (desc == size(qccli%clima%ana) .and. datoqui >= climaquif) ) then
1295
1296 if (c_e(qccli%clima%voldatiattrb(indcana &
1297 ,indctime,indclevel,indctimerange,indcdativarr,indcnetwork,1))) then
1298
1299 !ATTENZIONE TODO : inddativarr È UNA GRANDE SEMPLIFICAZIONE NON VERA SE TIPI DI DATO DIVERSI !!!!
1300
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
1306
1307#ifdef DEBUG
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)))
1316#endif
1317
1318
1319 if ( associated ( qccli%data_id_in)) then
1320#ifdef DEBUG
1321 call l4f_log (l4f_debug,"id: "//t2c(&
1322 qccli%data_id_in(indana,indtime,indlevel,indtimerange,indnetwork)))
1323#endif
1324 qccli%data_id_out(indana,indtime,indlevel,indtimerange,indnetwork)=&
1325 qccli%data_id_in(indana,indtime,indlevel,indtimerange,indnetwork)
1326 end if
1327 end if
1328 end if
1329!!$ end where
1330 end if
1331 end do
1332 end if
1333 end if
1334 end if
1335 end do
1336 end do
1337 end do
1338 end do
1339 end do
1340!!$ end forall
1341end do
1342
1343!!$print*,"risultato"
1344!!$print *,qccli%v7d%voldatiattrb(:,:,:,:,:,:,indbattrout)
1345!!$print*,"fine risultato"
1346
1347
1348return
1349
1350end subroutine quaconcli
1351
1352
1356subroutine cli_level(heigth,level)
1357
1358real,intent(in) :: heigth
1359TYPE(vol7d_level),intent(out):: level
1360
1361integer :: i
1362
1363i=imiss
1364
1365if (c_e(heigth)) then
1366 i=firsttrue(cli_level1 <= heigth .and. heigth <= cli_level2 )
1367end if
1368
1369if (i >= 1 .and. i <= 10 ) then
1370 call init(level, 102,cli_level1(i)*1000,102,cli_level2(i)*1000)
1371else
1372 if (c_e(i)) CALL l4f_log(l4f_debug,"cli_level: strange level, heigth: "//to_char(heigth))
1373 call init(level)
1374end if
1375
1376end subroutine cli_level
1377
1378
1380subroutine cli_level_generate(level)
1381
1382TYPE(vol7d_level),intent(out):: level(:)
1383
1384integer :: i
1385
1386if (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)))
1389end if
1390
1391do i=1,cli_nlevel
1392 call init(level(i), 102,cli_level1(i)*1000,102,cli_level2(i)*1000)
1393end do
1394
1395end subroutine cli_level_generate
1396
1397
1398!!$subroutine qccli_validate(qccli)
1399!!$type(qcclitype),intent(in) :: qccli
1400!!$
1401!!$!todo da validare
1402!!$
1403!!$return
1404!!$end subroutine qccli_validate
1405
1406
1409integer function supermacroa(macroa)
1410
1411integer, intent(in) :: macroa
1412 ! rielabora le macroarea facendole Valentine thinking
1413
1414supermacroa=imiss
1415
1416if (macroa == 1 .or. macroa == 2 .or. macroa == 4 ) supermacroa=3
1417if (macroa == 3 .or. macroa == 5 .or. macroa == 6 ) supermacroa=2
1418if (macroa == 7 .or. macroa == 8 ) supermacroa=1
1419
1420!!$ ! rielabora le macroarea facendole Valentine thinking
1421!!$
1422!!$ if (qccli%in_macroa(indana) == 1 .or. qccli%in_macroa(indana) == 2 .or. qccli%in_macroa(indana) == 4 ) iarea=3
1423!!$ if (qccli%in_macroa(indana) == 3 .or. qccli%in_macroa(indana) == 5 .or. qccli%in_macroa(indana) == 6 ) iarea=2
1424!!$ if (qccli%in_macroa(indana) == 7 .or. qccli%in_macroa(indana) == 8 ) iarea=1
1425
1426end function supermacroa
1427
1428
1429SUBROUTINE qc_compute_percentile(this, perc_vals,cyclicdt,presentperc, presentnumb)
1430
1431TYPE(qcclitype),INTENT(inout) :: this
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
1435real,intent(in) :: perc_vals(:)
1436TYPE(cyclicdatetime),INTENT(in) :: cyclicdt
1437real, optional :: presentperc
1438integer, optional :: presentnumb
1440
1441integer :: indana,indtime,indvar,indnetwork,indlevel ,indtimerange ,inddativarr,i,j,k,iana,narea
1442
1443REAL, DIMENSION(:),allocatable :: perc
1444TYPE(vol7d_var) :: var
1445character(len=vol7d_ana_lenident) :: ident
1446character(len=1) :: type
1447integer :: areav(size(this%v7d%ana)),iclv(size(this%v7d%ana))
1448real :: height
1449logical,allocatable :: mask(:,:,:),maskplus(:,:,:), maskarea(:)
1450integer,allocatable :: area(:)
1451real :: lpresentperc
1452integer :: lpresentnumb
1453
1454lpresentperc=optio_r(presentperc)
1455lpresentnumb=optio_i(presentnumb)
1456
1457allocate (perc(size(perc_vals)))
1458
1459call delete(this%extreme)
1460CALL init(this%extreme, time_definition=this%v7d%time_definition)
1461
1462call init(var, btable="B01192") ! MeteoDB station ID that here is the number of area
1463
1464type=cmiss
1465indvar = index(this%v7d%anavar, var, type=type)
1466indnetwork=min(1,size(this%v7d%network))
1467
1468if( indvar > 0 .and. indnetwork > 0 ) then
1469 select case (type)
1470 case("d")
1471 areav=integerdat(this%v7d%volanad(:,indvar,indnetwork),this%v7d%anavar%d(indvar))
1472 case("r")
1473 areav=integerdat(this%v7d%volanar(:,indvar,indnetwork),this%v7d%anavar%r(indvar))
1474 case("i")
1475 areav=integerdat(this%v7d%volanai(:,indvar,indnetwork),this%v7d%anavar%i(indvar))
1476 case("b")
1477 areav=integerdat(this%v7d%volanab(:,indvar,indnetwork),this%v7d%anavar%b(indvar))
1478 case("c")
1479 areav=integerdat(this%v7d%volanac(:,indvar,indnetwork),this%v7d%anavar%c(indvar))
1480 case default
1481 areav=imiss
1482 end select
1483else
1484 areav=imiss
1485end if
1486
1487allocate(maskarea(size(this%v7d%ana)))
1488maskarea(:)= areav(:) /= imiss
1489narea=count_distinct(areav,maskarea)
1490allocate(area(narea))
1491area=pack_distinct(areav,narea,maskarea)
1492deallocate(maskarea)
1493if (this%height2level) then
1494 call vol7d_alloc(this%extreme,nana=narea*size(perc_vals)*cli_nlevel)
1495else
1496 call vol7d_alloc(this%extreme,nana=narea*size(perc_vals))
1497endif
1498
1499if (this%height2level) then
1500
1501 call init(var, btable="B07030") ! height
1502
1503 type=cmiss
1504 indvar = index(this%v7d%anavar, var, type=type)
1505
1506!!$#ifdef DEBUG
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))
1511!!$ endif
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))
1516!!$ endif
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))
1521!!$ endif
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))
1526!!$ endif
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))
1531!!$ endif
1532!!$ CALL l4f_log(L4F_DEBUG, 'indvar has value '//t2c(indvar))
1533!!$ CALL l4f_log(L4F_DEBUG, 'indnetwork has value '//t2c(indnetwork))
1534!!$#endif
1536 do k=1,size(this%v7d%ana)
1537 height=rmiss
1538
1539 ! here we take the height fron any network (the first network win)
1540 do indnetwork=1,size(this%v7d%network)
1541
1542 if( indvar > 0 .and. indnetwork > 0 ) then
1543 select case (type)
1544 case("d")
1545 height=realdat(this%v7d%volanad(k,indvar,indnetwork),this%v7d%anavar%d(indvar))
1546 case("r")
1547 height=realdat(this%v7d%volanar(k,indvar,indnetwork),this%v7d%anavar%r(indvar))
1548 case ("i")
1549 height=realdat(this%v7d%volanai(k,indvar,indnetwork),this%v7d%anavar%i(indvar))
1550 case("b")
1551 height=realdat(this%v7d%volanab(k,indvar,indnetwork),this%v7d%anavar%b(indvar))
1552 case("c")
1553 height=realdat(this%v7d%volanac(k,indvar,indnetwork),this%v7d%anavar%c(indvar))
1554 case default
1555 height=rmiss
1556 end select
1557 end if
1558
1559 if (c_e(height)) exit
1560 end do
1561
1562 if (c_e(height)) then
1563 iclv(k)=firsttrue(cli_level1 <= height .and. height <= cli_level2 )
1564 else
1565 iclv(k)=imiss
1566 endif
1567
1568#ifdef DEBUG
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)))
1572#endif
1573 end do
1574
1575
1576endif
1577
1578do i=1,narea
1579 do j=1,size(perc_vals)
1580 if (this%height2level) then
1581 do k=1,cli_nlevel
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)
1584 enddo
1585 else
1586 k=0
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)
1589 endif
1590 end do
1591end do
1592
1593!!$do i=1,size(that%ana)
1594!!$ call display(that%ana(i))
1595!!$end do
1596
1597#ifdef DEBUG
1598CALL l4f_category_log(this%category, l4f_debug, 'nana has value '//t2c(SIZE(this%v7d%ana)))
1599CALL l4f_category_log(this%category, l4f_debug, 'lpresentperc has value '//t2c(lpresentperc))
1600CALL l4f_category_log(this%category, l4f_debug, 'lpresentnumb has value '//t2c(lpresentnumb))
1601#endif
1602
1603
1604call 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)
1606
1607this%extreme%level=this%v7d%level
1608this%extreme%timerange=this%v7d%timerange
1609this%extreme%dativar%r=this%v7d%dativar%r
1610this%extreme%time(1)=cyclicdatetime_to_conventional(cyclicdt)
1611call l4f_category_log(this%category, l4f_debug,"vol7d_compute_percentile conventional datetime "//to_char(this%extreme%time(1)))
1612call init(this%extreme%network(1),name="qcclima-perc")
1613
1614call vol7d_alloc_vol(this%extreme,inivol=.true.)
1615
1616allocate (mask(size(this%v7d%ana),size(this%v7d%time),size(this%v7d%network)))
1617
1618indtime=1
1619indnetwork=1
1620do 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
1623 do i=1,narea
1624
1625 !this%v7d%voldatir(indana, indtime, indlevel, indtimerange, inddativarr, indnetwork)
1626
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))
1629
1630#ifdef DEBUG
1631 CALL l4f_category_log(this%category, l4f_debug, 'count has value '//t2c(count &
1632 (mask .and. c_e(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:)))))
1633#endif
1634
1635 !delete in mask different area
1636 do j=1, size(mask,1)
1637 if (areav(j) /= area(i)) mask(j,:,:) =.false.
1638 end do
1639
1640 if (this%height2level) then
1641 allocate (maskplus(size(this%v7d%ana),size(this%v7d%time),size(this%v7d%network)))
1642 do k=1,cli_nlevel
1643#ifdef DEBUG
1644 CALL l4f_category_log(this%category, l4f_debug, 'k has value '//t2c(k))
1645#endif
1646
1647 do iana=1,size(mask,1)
1648 if (iclv(iana) /= k) maskplus(iana,:,:) =.false.
1649 if (iclv(iana) == k) maskplus(iana,:,:) = mask(iana,:,:)
1650 enddo
1651
1652 call sub_perc(maskplus)
1653
1654 enddo
1655 deallocate(maskplus)
1656 else
1657
1658 k=1
1659 call sub_perc(mask)
1660
1661 endif
1662 end do
1663 end do
1664 end do
1665end do
1666
1667deallocate (perc,mask,area)
1668
1669contains
1670
1671subroutine sub_perc(mymask)
1672
1673logical :: mymask(:,:,:)
1674
1675 ! we want more than 30% data present and a number of data bigger than 100 (default)
1676if &
1677 ( c_e(lpresentperc) .and. ((float(count &
1678 (mymask .and. c_e(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:)))&
1679 ) / &
1680 float(count(mymask))) < lpresentperc)) &
1681 return
1682
1683if &
1684 ( c_e(lpresentnumb) .and. (count &
1685 (mymask .and. c_e(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:))) < lpresentnumb)&
1686 ) &
1687 return
1688
1689
1690perc= stat_percentile(&
1691 pack(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:), &
1692 mask=mymask), &
1693 perc_vals)
1694
1695
1696do 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)=&
1699 perc(j)
1700enddo
1701
1702
1703end subroutine sub_perc
1704
1705
1706end SUBROUTINE qc_compute_percentile
1707
1708
1709SUBROUTINE qc_compute_normalizeddensityindex(this, perc_vals,cyclicdt,presentperc, presentnumb,data_normalized)
1710
1711TYPE(qcclitype),INTENT(inout) :: this
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
1715real,intent(in) :: perc_vals(:)
1716TYPE(cyclicdatetime),INTENT(in) :: cyclicdt
1717real,optional :: presentperc
1718integer,optional :: presentnumb
1719logical,optional,intent(in) :: data_normalized
1720
1721integer :: indana,indtime,indvar,indnetwork,indlevel ,indtimerange ,inddativarr, indattr
1722integer :: i,j,k,iana,narea
1723real :: height
1724TYPE(vol7d_var) :: var
1725character(len=vol7d_ana_lenident) :: ident
1726character(len=1) :: type
1727integer :: areav(size(this%v7d%ana)),iclv(size(this%v7d%ana))
1728logical,allocatable :: mask(:,:,:),maskplus(:,:,:), maskarea(:)
1729integer,allocatable :: area(:)
1730REAL, DIMENSION(:),allocatable :: ndi,limbins
1731real :: lpresentperc
1732integer :: lpresentnumb
1733logical :: lnorm
1734
1735lnorm = optio_log(data_normalized)
1736lpresentperc=optio_r(presentperc)
1737lpresentnumb=optio_i(presentnumb)
1738
1739allocate (ndi(size(perc_vals)-1),limbins(size(perc_vals)))
1740call delete(this%clima)
1741CALL init(this%clima, time_definition=this%v7d%time_definition)
1742call init(var, btable="B01192") ! MeteoDB station ID that here is the number of area
1743
1744if (.NOT.(lnorm)) then
1745
1746 type=cmiss
1747 indvar = index(this%v7d%anavar, var, type=type)
1748 indnetwork=1
1749
1750!if( ind /= 0 ) then
1751 select case (type)
1752 case("d")
1753 areav=integerdat(this%v7d%volanad(:,indvar,indnetwork),this%v7d%anavar%d(indvar))
1754 case("r")
1755 areav=integerdat(this%v7d%volanar(:,indvar,indnetwork),this%v7d%anavar%r(indvar))
1756 case("i")
1757 areav=integerdat(this%v7d%volanai(:,indvar,indnetwork),this%v7d%anavar%i(indvar))
1758 case("b")
1759 areav=integerdat(this%v7d%volanab(:,indvar,indnetwork),this%v7d%anavar%b(indvar))
1760 case("c")
1761 areav=integerdat(this%v7d%volanac(:,indvar,indnetwork),this%v7d%anavar%c(indvar))
1762 case default
1763 areav=imiss
1764 end select
1765!end if
1766
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)
1773
1774 if (this%height2level) then
1775 call vol7d_alloc(this%clima,nana=narea*(size(perc_vals)-1)*cli_nlevel)
1776 else
1777 call vol7d_alloc(this%clima,nana=narea*(size(perc_vals)-1))
1778 endif
1779
1780
1781
1782
1783 do i=1,narea
1784 do j=1,size(perc_vals)-1
1785 if (this%height2level) then
1786 do k=1,cli_nlevel
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)
1789 enddo
1790 else
1791 k=0
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)
1794 endif
1795 end do
1796 end do
1797
1798
1799
1800
1801
1802 if (this%height2level) then
1803
1804 call init(var, btable="B07030") ! height
1805
1806 type=cmiss
1807 indvar = index(this%v7d%anavar, var, type=type)
1808
1809 do k=1,size(this%v7d%ana)
1810 height=rmiss
1811
1812 ! here we take the height fron any network (the first network win)
1813 do indnetwork=1,size(this%v7d%network)
1814
1815 if( indvar > 0 .and. indnetwork > 0 ) then
1816 select case (type)
1817 case("d")
1818 height=realdat(this%v7d%volanad(k,indvar,indnetwork),this%v7d%anavar%d(indvar))
1819 case("r")
1820 height=realdat(this%v7d%volanar(k,indvar,indnetwork),this%v7d%anavar%r(indvar))
1821 case ("i")
1822 height=realdat(this%v7d%volanai(k,indvar,indnetwork),this%v7d%anavar%i(indvar))
1823 case("b")
1824 height=realdat(this%v7d%volanab(k,indvar,indnetwork),this%v7d%anavar%b(indvar))
1825 case("c")
1826 height=realdat(this%v7d%volanac(k,indvar,indnetwork),this%v7d%anavar%c(indvar))
1827 case default
1828 height=rmiss
1829 end select
1830 end if
1831
1832 if (c_e(height)) exit
1833 end do
1834
1835 if (c_e(height)) then
1836 iclv(k)=firsttrue(cli_level1 <= height .and. height <= cli_level2 )
1837 else
1838 iclv(k)=imiss
1839 endif
1840
1841#ifdef DEBUG
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)))
1845#endif
1846 end do
1847
1848
1849 endif
1850
1851
1852else
1853
1854 narea=1
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)
1859 enddo
1860
1861endif
1862
1863
1864!!$do i=1,size(that%ana)
1865!!$ call display(that%ana(i))
1866!!$end do
1867
1868call 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)
1870
1871this%clima%level=this%v7d%level
1872this%clima%timerange=this%v7d%timerange
1873this%clima%dativar%r=this%v7d%dativar%r
1874this%clima%dativarattr%r=this%clima%dativar%r
1875call init(this%clima%datiattr%r(1), btable="*B33209") ! NDI order number
1876this%clima%time(1)=cyclicdatetime_to_conventional(cyclicdt)
1877
1878call l4f_category_log(this%category,l4f_info,"vol7d_compute_ndi conventional datetime "//to_char(this%clima%time(1)))
1879call init(this%clima%network(1),name="qcclima-ndi")
1880
1881call vol7d_alloc_vol(this%clima,inivol=.true.)
1882
1883allocate (mask(size(this%v7d%ana),size(this%v7d%time),size(this%v7d%network)))
1884
1885indtime=1
1886indnetwork=1
1887indattr=1
1888do 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 if (.NOT.(lnorm)) then
1892 do i=1,narea
1893
1894 !this%v7d%voldatir(indana, indtime, indlevel, indtimerange, inddativarr, indnetwork)
1895
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.
1901 end do
1902
1903 if (this%height2level) then
1904 allocate (maskplus(size(this%v7d%ana),size(this%v7d%time),size(this%v7d%network)))
1905 do k=1,cli_nlevel
1906#ifdef DEBUG
1907 CALL l4f_category_log(this%category, l4f_debug, 'k has value '//t2c(k))
1908#endif
1909
1910 do iana=1,size(mask,1)
1911 if (iclv(iana) /= k) maskplus(iana,:,:) =.false.
1912 if (iclv(iana) == k) maskplus(iana,:,:) = mask(iana,:,:)
1913 enddo
1914
1915 call sub_ndi(mymask=maskplus)
1916
1917 enddo
1918 deallocate(maskplus)
1919
1920 else
1921
1922 k=1
1923 call sub_ndi(mymask=mask)
1924
1925 endif
1926
1927 end do
1928
1929 else
1930
1931 mask = spread(spread((this%v7d%time == cyclicdt ),1,size(this%v7d%ana)),3,size(this%v7d%network))
1932
1933 k=1
1934 i=1
1935 call sub_ndi(mymask=mask)
1936
1937 endif
1938
1939 end do
1940 end do
1941end do
1942
1943deallocate (ndi,limbins,mask)
1944if (.NOT.(lnorm)) deallocate(area)
1945
1946contains
1947
1948subroutine sub_ndi(mymask)
1949
1950logical :: mymask(:,:,:)
1951
1952 ! we want more than 30% data present and a number of data bigger than 100 (default)
1953if &
1954 ( c_e(lpresentperc) .and. ((float(count &
1955 (mymask .and. c_e(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:)))&
1956 ) / &
1957 float(count(mymask))) < lpresentperc)) &
1958 return
1959
1960if &
1961 ( c_e(lpresentnumb) .and. (count &
1962 (mymask .and. c_e(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:))) < lpresentnumb)&
1963 ) &
1964 return
1965
1966!print*,"compute"
1967!print*,"-------------------------------------------------------------"
1968
1969call normalizeddensityindex (&
1970 pack(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:), &
1971 mask=mymask), &
1972 perc_vals, ndi, limbins)
1973
1974do 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)=&
1977 limbins(j)
1978 this%clima%voldatiattrr(indana, indtime, indlevel, indtimerange, inddativarr, indnetwork,indattr)=&
1979 ndi(j)*100
1980end do
1981
1982end subroutine sub_ndi
1983
1984
1985end SUBROUTINE qc_compute_normalizeddensityindex
1986
1987end module modqccli
1988
1989
1992
Set of functions that return a trimmed CHARACTER representation of the input variable.
Set of functions that return a CHARACTER representation of the input variable.
Operatore di valore assoluto di un intervallo.
Restituiscono il valore dell'oggetto nella forma desiderata.
Import an array of georef_coord_array objects from a file in ESRI/Shapefile format.
Determine whether a point lies inside a polygon or a rectangle.
Index method.
Emit log message for a category with specific priority.
Test di dato invalidato.
Definition: modqc.F90:411
Check data validity based on gross error check.
Definition: modqc.F90:406
Allocazione di memoria.
Definition: modqccli.F90:342
Cancellazione.
Definition: modqccli.F90:347
Inizializzazione.
Definition: modqccli.F90:337
Compute a set of percentiles for a random variable.
integer data conversion
real data conversion
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.
Definition: modqc.F90:357
Controllo di qualità climatico.
Definition: modqccli.F90:295
Module for basic statistical computations taking into account missing data.
Definition: simple_stat.f90:25
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 a cyclic datetime.
Class for expressing an absolute time value.
Derived type defining a dynamically extensible array of TYPE(georef_coord_array) elements.
Derive type defining a single georeferenced point, either in geodetic or in projected coordinates.
Oggetto principale per il controllo di qualità
Definition: modqccli.F90:322
Definisce un oggetto contenente i volumi anagrafica e dati e tutti i descrittori delle loro dimension...

Generated with Doxygen.