libsim  Versione7.1.6
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 
107 module modqccli
108 
110 use vol7d_class
111 use modqc
112 use file_utilities
113 use log4fortran
114 use char_utilities
115 use datetime_class
116 use simple_stat
117 !use array_utilities
118 !use io_units
119 #ifdef HAVE_DBALLE
121 #endif
122 
123 implicit none
124 
125 character (len=255),parameter:: subcategory="QCcli"
126 
127 integer, parameter :: cli_nlevel=10
128 
129 integer, parameter :: cli_level1(cli_nlevel) = (/-100,100,250,500,750,1000,1250,1500,1750,2000/)
130 
131 integer, parameter :: cli_level2(cli_nlevel) = (/100,250,500,750,1000,1250,1500,1750,2000,2250/)
132 
134 type :: 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 
145 end type qcclitype
146 
147 
149 interface init
150  module procedure qccliinit
151 end interface
152 
154 interface alloc
155  module procedure qcclialloc
156 end interface
157 
159 interface delete
160  module procedure qcclidelete
161 end interface
162 
163 PRIVATE
164 PUBLIC 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 
168 contains
169 
172 
173 subroutine 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 
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
194 
195 #ifdef HAVE_DBALLE
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
209 #endif
210 
211 integer :: iuni,i,j
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
218 
219 ! temporary, improve!!!!
220 DOUBLE PRECISION :: lon, lat
221 TYPE(georef_coord) :: point
222 TYPE(arrayof_georef_coord_array) :: macroa
223 
224 
225 call l4f_launcher(a_name,a_name_append=trim(subcategory)//"."//trim(categoryappend))
226 qccli%category=l4f_category_get(a_name)
227 
228 qccli%height2level=optio_log(height2level)
229 
230 ! all coordinates in clima and extreme datasets are zero
231 call init(lcoordmin)
232 call 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 
240 nullify( qccli%in_macroa )
241 nullify( qccli%data_id_in )
242 nullify( qccli%data_id_out )
243 
244 ! riporto il volume dati nel mio oggetto
245 qccli%v7d => v7d
246 
247 if (present(data_id_in))then
248  qccli%data_id_in => data_id_in
249 end if
250 
251 if (qccli%height2level) then
252  !shape file for Emilia Romagna clima !
253  filepath=get_package_filepath('share:sup_macroaree.shp', filetype_data)
254 else
255  !shape file for Europa clima !
256  filepath=get_package_filepath('share:ens_v8_ll.shp', filetype_data)
257 end if
258 
259 if (present(macropath))then
260  if (c_e(macropath)) then
261  filepath=macropath
262  end if
263 end if
264 
265 CALL import(macroa, shpfile=filepath)
266 
267 call optio(climapath,filepathclima)
268 call optio(extremepath,filepathextreme)
269 
270 #ifdef HAVE_DBALLE
271 
272 call init(network,"qcclima-ndi")
273 
274 ltimei=datetime_miss
275 ltimef=datetime_miss
276 if (present(timei)) then
277  ltimei=timei
278  ltimei=ltimei+timedelta_new(minute=30)
279 end if
280 
281 if (present(timef)) then
282  ltimef=timef
283  ltimef=ltimef+timedelta_new(minute=30)
284 end if
285 
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)
288 
289 if ( 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))
295 
296 else
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 
301 end if
302 
303 call optio(dsncli,ldsncli)
304 call optio(user,luser)
305 call optio(password,lpassword)
306 
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")
309  call raise_error()
310 end if
311 
312 if (.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
322  select case (trim(lowercase(suffixname(filepathclima))))
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
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
352 else
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 
363 end if
364 #endif
365 
366 #ifdef HAVE_DBALLE
367 call delete(ltimei)
368 call delete(ltimef)
369 
370 if ( 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 
384 else
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 
389 end if
390 
391 call init(network,"qcclima-perc")
392 call optio(dsnextreme,ldsnextreme)
393 
394 if (.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
438 else
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 
450 end if
451 
452 call delete(ltimei)
453 call delete(ltimef)
454 #endif
455 
456 
457 call 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 
467 if (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
481 end if
482 
483 call delete(macroa)
484 
485 return
486 end subroutine qccliinit
487 
488 
490 subroutine qcclialloc(qccli)
491  ! pseudo costruttore con distruttore automatico
492 
493 type(qcclitype),intent(in out) :: qccli
494 
495 integer :: istatt
496 integer :: sh(5)
497 
498 ! se ti sei dimenticato di deallocare ci penso io
499 call 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 
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")
519  end if
520 end if
521 
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)
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
531 end if
532 
533 return
534 
535 end subroutine qcclialloc
536 
537 
539 
540 subroutine qcclidealloc(qccli)
541  ! pseudo distruttore
542 
543 type(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 
553 if (associated (qccli%in_macroa)) then
554  deallocate (qccli%in_macroa)
555 end if
556 
557 if (associated(qccli%data_id_out))then
558  deallocate (qccli%data_id_out)
559 end if
560 
561 return
562 end subroutine qcclidealloc
563 
564 
566 
567 
568 subroutine qcclidelete(qccli)
569  ! decostruttore a mezzo
570 type(qcclitype),intent(in out) :: qccli
571 
572 call qcclidealloc(qccli)
573 
574 call delete(qccli%clima)
575 call delete(qccli%extreme)
576 
577 !delete logger
578 call l4f_category_delete(qccli%category)
579 
580 return
581 end subroutine qcclidelete
582 
583 
584 
597 SUBROUTINE vol7d_normalize_data(qccli)
598 
599 TYPE(qcclitype),INTENT(inout) :: qccli
600 !character (len=10) ,intent(in),optional :: battrinv !< attributo invalidated in input/output
601 
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))
607 real :: height
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
614 
615 integer :: clev
616 character(len=1) :: mycanc, canc = "#"
617 
618 CHARACTER(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 
631 if (.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
636 end if
637 
638 
639 if (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 
682 endif
683 
684 
685 if (.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
690 end if
691 
692 if (.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
697 end if
698 
699 do 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)
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)
728 
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 
850 end do
851 
852 end SUBROUTINE vol7d_normalize_data
853 
854 
855 real function base_value(btable)
856 character (len=10) ,intent(in):: btable
857 
858 character (len=10) :: btables(1) =(/"B12101"/)
859 real :: base_values(1) =(/273.15 /)
860 integer :: ind
861 
862 ind = index_c(btables,btable)
863 
864 if (ind > 0) then
865  base_value = base_values(ind)
866 else
867  call l4f_log(l4f_warn,"modqccli_base_value: variable "//btable//" do not have base value")
868  base_value = 0.
869 end if
870 
871 return
872 
873 end function base_value
874 
875 
882 SUBROUTINE quaconcli (qccli,battrinv,battrout,&
883  anamask,timemask,levelmask,timerangemask,varmask,networkmask)
884 
885 
886 type(qcclitype),intent(in out) :: qccli
887 character (len=10) ,intent(in),optional :: battrinv
888 character (len=10) ,intent(in),optional :: battrout
889 logical ,intent(in),optional :: anamask(:)
890 logical ,intent(in),optional :: timemask(:)
891 logical ,intent(in),optional :: levelmask(:)
892 logical ,intent(in),optional :: timerangemask(:)
893 logical ,intent(in),optional :: varmask(:)
894 logical ,intent(in),optional :: networkmask(:)
895 
896 CHARACTER(len=vol7d_ana_lenident) :: ident
897 REAL(kind=fp_geo) :: latc,lonc
898 integer :: mese, ora
899  !local
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))
903 
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
908 real :: height
909 !integer, allocatable :: indcanav(:)
910 
911 
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
919 
920 
921 !call qccli_validate (qccli)
922 
923 indbattrinv=0
924 if (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
930 end if
931 
932 if ( 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 
937 end if
938 
939 
940 if (present(battrout))then
941  indbattrout = index_c(qccli%v7d%datiattr%b(:)%btable, battrout)
942 else
943  indbattrout = index_c(qccli%v7d%datiattr%b(:)%btable, qcattrvarsbtables(2))
944 end if
945 
946 if ( 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 
951 end if
952 
953 if(present(anamask)) then
954  anamaskl = anamask
955 else
956  anamaskl = .true.
957 endif
958 if(present(timemask)) then
959  timemaskl = timemask
960 else
961  timemaskl = .true.
962 endif
963 if(present(levelmask)) then
964  levelmaskl = levelmask
965 else
966  levelmaskl = .true.
967 endif
968 if(present(timerangemask)) then
969  timerangemaskl = timerangemask
970 else
971  timerangemaskl = .true.
972 endif
973 if(present(varmask)) then
974  varmaskl = varmask
975 else
976  varmaskl = .true.
977 endif
978 if(present(networkmask)) then
979  networkmaskl = networkmask
980 else
981  networkmaskl = .true.
982 endif
983 
984 qccli%v7d%voldatiattrb(:,:,:,:,:,:,indbattrout)=ibmiss
985 call init(anavar,"B07030" )
986 type=cmiss
987 indvar = index(qccli%v7d%anavar, anavar, type=type)
988 
989 do 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
1341 end do
1342 
1343 !!$print*,"risultato"
1344 !!$print *,qccli%v7d%voldatiattrb(:,:,:,:,:,:,indbattrout)
1345 !!$print*,"fine risultato"
1346 
1347 
1348 return
1349 
1350 end subroutine quaconcli
1351 
1352 
1356 subroutine cli_level(heigth,level)
1357 
1358 real,intent(in) :: heigth
1359 type(vol7d_level),intent(out):: level
1360 
1361 integer :: i
1362 
1363 i=imiss
1364 
1365 if (c_e(heigth)) then
1366  i=firsttrue(cli_level1 <= heigth .and. heigth <= cli_level2 )
1367 end if
1368 
1369 if (i >= 1 .and. i <= 10 ) then
1370  call init(level, 102,cli_level1(i)*1000,102,cli_level2(i)*1000)
1371 else
1372  if (c_e(i)) CALL l4f_log(l4f_debug,"cli_level: strange level, heigth: "//to_char(heigth))
1373  call init(level)
1374 end if
1375 
1376 end subroutine cli_level
1377 
1378 
1380 subroutine cli_level_generate(level)
1381 
1382 TYPE(vol7d_level),intent(out):: level(:)
1383 
1384 integer :: i
1385 
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)))
1389 end if
1390 
1391 do i=1,cli_nlevel
1392  call init(level(i), 102,cli_level1(i)*1000,102,cli_level2(i)*1000)
1393 end do
1394 
1395 end 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 
1409 integer function supermacroa(macroa)
1410 
1411 integer, intent(in) :: macroa
1412  ! rielabora le macroarea facendole Valentine thinking
1413 
1414 supermacroa=imiss
1415 
1416 if (macroa == 1 .or. macroa == 2 .or. macroa == 4 ) supermacroa=3
1417 if (macroa == 3 .or. macroa == 5 .or. macroa == 6 ) supermacroa=2
1418 if (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 
1426 end function supermacroa
1427 
1428 
1429 SUBROUTINE qc_compute_percentile(this, perc_vals,cyclicdt,presentperc, presentnumb)
1430 
1431 TYPE(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
1435 real,intent(in) :: perc_vals(:)
1436 type(cyclicdatetime),INTENT(in) :: cyclicdt
1437 real, optional :: presentperc
1438 integer, optional :: presentnumb
1440 
1441 integer :: indana,indtime,indvar,indnetwork,indlevel ,indtimerange ,inddativarr,i,j,k,iana,narea
1442 
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))
1448 real :: height
1449 logical,allocatable :: mask(:,:,:),maskplus(:,:,:), maskarea(:)
1450 integer,allocatable :: area(:)
1451 real :: lpresentperc
1452 integer :: lpresentnumb
1453 
1454 lpresentperc=optio_r(presentperc)
1455 lpresentnumb=optio_i(presentnumb)
1456 
1457 allocate (perc(size(perc_vals)))
1458 
1459 call delete(this%extreme)
1460 CALL init(this%extreme, time_definition=this%v7d%time_definition)
1461 
1462 call init(var, btable="B01192") ! MeteoDB station ID that here is the number of area
1463 
1464 type=cmiss
1465 indvar = index(this%v7d%anavar, var, type=type)
1466 indnetwork=min(1,size(this%v7d%network))
1467 
1468 if( 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
1483 else
1484  areav=imiss
1485 end if
1486 
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)
1495 else
1496  call vol7d_alloc(this%extreme,nana=narea*size(perc_vals))
1497 endif
1498 
1499 if (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
1535 
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 
1576 endif
1577 
1578 do 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
1591 end do
1592 
1593 !!$do i=1,size(that%ana)
1594 !!$ call display(that%ana(i))
1595 !!$end do
1596 
1597 #ifdef DEBUG
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))
1601 #endif
1602 
1603 
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)
1606 
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")
1613 
1614 call vol7d_alloc_vol(this%extreme,inivol=.true.)
1615 
1616 allocate (mask(size(this%v7d%ana),size(this%v7d%time),size(this%v7d%network)))
1617 
1618 indtime=1
1619 indnetwork=1
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
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
1665 end do
1666 
1667 deallocate (perc,mask,area)
1668 
1669 contains
1670 
1671 subroutine sub_perc(mymask)
1672 
1673 logical :: mymask(:,:,:)
1674 
1675  ! we want more than 30% data present and a number of data bigger than 100 (default)
1676 if &
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 
1683 if &
1684  ( c_e(lpresentnumb) .and. (count &
1685  (mymask .and. c_e(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:))) < lpresentnumb)&
1686  ) &
1687  return
1688 
1689 
1690 perc= stat_percentile(&
1691  pack(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:), &
1692  mask=mymask), &
1693  perc_vals)
1694 
1695 
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)=&
1699  perc(j)
1700 enddo
1701 
1702 
1703 end subroutine sub_perc
1704 
1705 
1706 end SUBROUTINE qc_compute_percentile
1707 
1708 
1709 SUBROUTINE qc_compute_normalizeddensityindex(this, perc_vals,cyclicdt,presentperc, presentnumb,data_normalized)
1710 
1711 TYPE(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
1715 real,intent(in) :: perc_vals(:)
1716 type(cyclicdatetime),INTENT(in) :: cyclicdt
1717 real,optional :: presentperc
1718 integer,optional :: presentnumb
1719 logical,optional,intent(in) :: data_normalized
1720 
1721 integer :: indana,indtime,indvar,indnetwork,indlevel ,indtimerange ,inddativarr, indattr
1722 integer :: i,j,k,iana,narea
1723 real :: height
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
1733 logical :: lnorm
1734 
1735 lnorm = optio_log(data_normalized)
1736 lpresentperc=optio_r(presentperc)
1737 lpresentnumb=optio_i(presentnumb)
1738 
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
1743 
1744 if (.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 
1852 else
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 
1861 endif
1862 
1863 
1864 !!$do i=1,size(that%ana)
1865 !!$ call display(that%ana(i))
1866 !!$end do
1867 
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)
1870 
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)
1877 
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")
1880 
1881 call vol7d_alloc_vol(this%clima,inivol=.true.)
1882 
1883 allocate (mask(size(this%v7d%ana),size(this%v7d%time),size(this%v7d%network)))
1884 
1885 indtime=1
1886 indnetwork=1
1887 indattr=1
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  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
1941 end do
1942 
1943 deallocate (ndi,limbins,mask)
1944 if (.NOT.(lnorm)) deallocate(area)
1945 
1946 contains
1947 
1948 subroutine sub_ndi(mymask)
1949 
1950 logical :: mymask(:,:,:)
1951 
1952  ! we want more than 30% data present and a number of data bigger than 100 (default)
1953 if &
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 
1960 if &
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 
1969 call normalizeddensityindex(&
1970  pack(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:), &
1971  mask=mymask), &
1972  perc_vals, ndi, limbins)
1973 
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)=&
1977  limbins(j)
1978  this%clima%voldatiattrr(indana, indtime, indlevel, indtimerange, inddativarr, indnetwork,indattr)=&
1979  ndi(j)*100
1980 end do
1981 
1982 end subroutine sub_ndi
1983 
1984 
1985 end SUBROUTINE qc_compute_normalizeddensityindex
1986 
1987 end module modqccli
1988 
1989 
1992 
Classi per la gestione delle coordinate temporali.
Functions that return a trimmed CHARACTER representation of the input variable.
Check data validity based on gross error check.
Definition: modqc.F90:403
Operatore di valore assoluto di un intervallo.
Compute a set of percentiles for a random variable.
Oggetto per import ed export da DB-All.e.
Derived type defining a dynamically extensible array of TYPE(georef_coord_array) elements.
Inizializzazione.
Definition: modqccli.F90:334
Module for basic statistical computations taking into account missing data.
Definition: simple_stat.f90:25
Determine whether a point lies inside a polygon or a rectangle.
Generic subroutine for checking OPTIONAL parameters.
Import an array of georef_coord_array objects from a file in ESRI/Shapefile format.
Classe per la gestione di un volume completo di dati osservati.
Index method.
Restituiscono il valore dell&#39;oggetto in forma di stringa stampabile.
Restituiscono il valore dell&#39;oggetto nella forma desiderata.
Definisce un oggetto contenente i volumi anagrafica e dati e tutti i descrittori delle loro dimension...
Controllo di qualità climatico.
Definition: modqccli.F90:292
Oggetto principale per il controllo di qualità
Definition: modqccli.F90:319
Test di dato invalidato.
Definition: modqc.F90:408
Allocazione di memoria.
Definition: modqccli.F90:339
Copy an object, creating a fully new instance.
Cancellazione.
Definition: modqccli.F90:344
classe per import ed export di volumi da e in DB-All.e
Utilities and defines for quality control.
Definition: modqc.F90:354
This module defines objects describing georeferenced sparse points possibly with topology and project...
classe per la gestione del logging
Utilities for CHARACTER variables.
Derive type defining a single georeferenced point, either in geodetic or in projected coordinates...
Class for expressing a cyclic datetime.
real data conversion
Utilities for managing files.
integer data conversion
Emit log message for a category with specific priority.

Generated with Doxygen.