370if (
c_e(yeari) .and.
c_e(yearf) .and. yeari == yearf .and. monthi == monthf )
then
372 if ( dayi == dayf .and. houri == hourf .and. minutei == minutef .and. mseci == msecf )
then
374 ltimei=cyclicdatetime_to_conventional(cyclicdatetime_new(month=monthi, hour=houri))
375 ltimef=cyclicdatetime_to_conventional(cyclicdatetime_new(month=monthf, hour=hourf))
379 ltimei=cyclicdatetime_to_conventional(cyclicdatetime_new(month=monthi, hour=00))
380 ltimef=cyclicdatetime_to_conventional(cyclicdatetime_new(month=monthf, hour=23))
391call init(network,
"qcclima-perc")
392call optio(dsnextreme,ldsnextreme)
394if (.not.
c_e(ldsnextreme))
then
398 if (.not.
c_e(filepathextreme))
then
399 filepathextreme=get_package_filepath(
'qcclima-extreme.v7d', filetype_data)
402 if (
c_e(filepathextreme))
then
404 select case (trim(lowercase(suffixname(filepathextreme))))
408 call import(qccli%extreme,filename=filepathextreme,unit=iuni)
413 call init(v7d_dballeextreme,file=.true.,filename=filepathextreme,categoryappend=trim(a_name)//
".climaextreme")
415 call import(v7d_dballeextreme,var=var,coordmin=lcoordmin, coordmax=lcoordmax, timei=ltimei, timef=ltimef, &
416 varkind=(/(
"r",i=1,
size(var))/),attr=(/qcattrvarsbtables(2)/),attrkind=(/
"b"/),network=network)
417 call copy(v7d_dballeextreme%vol7d,qccli%extreme)
418 call delete(v7d_dballeextreme)
423 if (
c_e(filepathextreme))
then
425 "file type not supported (user .v7d or .bufr suffix only): "//trim(filepathextreme))
431 call l4f_category_log(qccli%category,l4f_warn,
"extreme volume not iniziatized: QC or normalize data will not be possible")
433 call init(qccli%extreme)
441 call init(v7d_dballeextreme,dsn=ldsnextreme,user=luser,password=lpassword,&
442 write=.false.,file=.false.,categoryappend=trim(a_name)//
".climaextreme")
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)
457call qcclialloc(qccli)
467if (
associated(qccli%in_macroa))
then
468 qccli%in_macroa = imiss
470 DO i = 1,
SIZE(qccli%v7d%ana)
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
486end subroutine qccliinit
490subroutine qcclialloc(qccli)
499call qcclidealloc(qccli)
514if (
associated (qccli%v7d%ana ))
then
515 allocate (qccli%in_macroa(
size(qccli%v7d%ana )),stat=istatt)
516 if (istatt /= 0)
then
518 call raise_error(
"allocate error")
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)
527 call raise_error(
"allocate error")
529 qccli%data_id_out=imiss
535end subroutine qcclialloc
540subroutine qcclidealloc(qccli)
553if (
associated (qccli%in_macroa))
then
554 deallocate (qccli%in_macroa)
557if (
associated(qccli%data_id_out))
then
558 deallocate (qccli%data_id_out)
562end subroutine qcclidealloc
568subroutine qcclidelete(qccli)
572call qcclidealloc(qccli)
578call l4f_category_delete(qccli%category)
581end subroutine qcclidelete
597SUBROUTINE vol7d_normalize_data(qccli)
602real :: datoqui, perc25, perc50,perc75
603integer :: indana , indtime ,indlevel ,indtimerange ,inddativarr, indnetwork
604integer :: indcana, indvar,indctime,indclevel,indctimerange,indcdativarr,indcnetwork
606integer :: iclv(size(qccli%v7d%ana))
608character(len=1) :: type
609TYPE(vol7d_var) :: var
610TYPE(vol7d_ana) :: ana
612TYPE(vol7d_level):: level
613integer :: mese, ora, desc, iarea, k
616character(len=1) :: mycanc, canc =
"#"
618CHARACTER(len=vol7d_ana_lenident) :: ident
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
639if (qccli%height2level)
then
640 call init(var, btable=
"B07030")
643 indvar =
index(qccli%v7d%anavar, var, type=type)
645 do indana=1,
size(qccli%v7d%ana)
649 do indnetwork=1,
size(qccli%v7d%network)
651 if( indvar > 0 )
then
654 height=
realdat(qccli%v7d%volanad(indana,indvar,indnetwork),qccli%v7d%anavar%d(indvar))
656 height=
realdat(qccli%v7d%volanar(indana,indvar,indnetwork),qccli%v7d%anavar%r(indvar))
658 height=
realdat(qccli%v7d%volanai(indana,indvar,indnetwork),qccli%v7d%anavar%i(indvar))
660 height=
realdat(qccli%v7d%volanab(indana,indvar,indnetwork),qccli%v7d%anavar%b(indvar))
662 height=
realdat(qccli%v7d%volanac(indana,indvar,indnetwork),qccli%v7d%anavar%c(indvar))
666 if (
c_e(height))
exit
669 if (
c_e(height))
then
670 iclv(indana)=firsttrue(cli_level1 <= height .and. height <= cli_level2 )
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)//&
790 if (indcana > 0 )
then
791 perc25=qccli%extreme%voldatir(indcana,indctime,indclevel,&
792 indctimerange,indcdativarr,indcnetwork)
796 write(ident,
'("#",i2.2,2i3.3)')k,iarea,desc
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)//&
802 if (indcana > 0 )
then
803 perc50=qccli%extreme%voldatir(indcana,indctime,indclevel,&
804 indctimerange,indcdativarr,indcnetwork)
808 write(ident,
'("#",i2.2,2i3.3)')k,iarea,desc
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)//&
814 if (indcana > 0 )
then
815 perc75=qccli%extreme%voldatir(indcana,indctime,indclevel,&
816 indctimerange,indcdativarr,indcnetwork)
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
823 call l4f_log(l4f_debug,
"normalize data dato in : "//
t2c(datoqui))
824 datoqui = (datoqui - perc50) / (perc75 - perc25) + &
825 base_value(qccli%v7d%dativar%r(inddativarr)%btable)
826 call l4f_log(l4f_debug,
"normalize data dato out : "//
t2c(datoqui))
834 qccli%v7d%voldatir (indana ,indtime ,indlevel ,indtimerange ,&
835 inddativarr, indnetwork ) = datoqui
843 read (qccli%v7d%ana(indana)%ident,
'(a1,i2.2,2i3.3)') mycanc, clev, iarea, desc
844 if (mycanc == canc)
then
847 write (qccli%v7d%ana(indana)%ident,
'(a1,i2.2,2i3.3)') canc, clev, iarea, desc
852end SUBROUTINE vol7d_normalize_data
855real function base_value(btable)
856character (len=10) ,
intent(in):: btable
858character (len=10) :: btables(1) =(/
"B12101"/)
859real :: base_values(1) =(/273.15 /)
862ind = index_c(btables,btable)
865 base_value = base_values(ind)
867 call l4f_log(l4f_warn,
"modqccli_base_value: variable "//btable//
" do not have base value")
873end function base_value
882SUBROUTINE quaconcli (qccli,battrinv,battrout,&
883 anamask,timemask,levelmask,timerangemask,varmask,networkmask)
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(:)
896CHARACTER(len=vol7d_ana_lenident) :: ident
897REAL(kind=fp_geo) :: latc,lonc
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))
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
912TYPE(vol7d_network) :: network
913TYPE(vol7d_ana) :: ana
915TYPE(vol7d_level):: level
916type(vol7d_var) :: anavar
917integer :: iclv(size(qccli%v7d%ana))
918character(len=1) :: type
924if (
associated(qccli%v7d%datiattr%b))
then
925 if (
present(battrinv))
then
926 indbattrinv = index_c(qccli%v7d%datiattr%b(:)%btable, battrinv)
928 indbattrinv = index_c(qccli%v7d%datiattr%b(:)%btable, qcattrvarsbtables(1))
932if ( indbattrinv <= 0 )
then
934 call l4f_category_log(qccli%category,l4f_error,
"error finding attribute index in/out "//qcattrvarsbtables(1))
935 call raise_error(
"error finding attribute index in/out "//qcattrvarsbtables(1))
940if (
present(battrout))
then
941 indbattrout = index_c(qccli%v7d%datiattr%b(:)%btable, battrout)
943 indbattrout = index_c(qccli%v7d%datiattr%b(:)%btable, qcattrvarsbtables(2))
946if ( indbattrout <= 0 )
then
948 call l4f_category_log(qccli%category,l4f_error,
"error finding attribute index in/out "//qcattrvarsbtables(2))
949 call raise_error(
"error finding attribute index in/out "//qcattrvarsbtables(2))
953if(
present(anamask))
then
958if(
present(timemask))
then
963if(
present(levelmask))
then
964 levelmaskl = levelmask
968if(
present(timerangemask))
then
969 timerangemaskl = timerangemask
971 timerangemaskl = .true.
973if(
present(varmask))
then
978if(
present(networkmask))
then
979 networkmaskl = networkmask
981 networkmaskl = .true.
984qccli%v7d%voldatiattrb(:,:,:,:,:,:,indbattrout)=ibmiss
985call init(anavar,
"B07030" )
987indvar =
index(qccli%v7d%anavar, anavar, type=type)
989do indana=1,
size(qccli%v7d%ana)
994 iarea= qccli%in_macroa(indana)
996 if (.not.
c_e(iarea)) cycle
1013 if (qccli%height2level)
then
1018 do indn=1,
size(qccli%v7d%network)
1020 if( indvar > 0 .and. indn > 0 )
then
1023 height=
realdat(qccli%v7d%volanad(indana,indvar,indn),qccli%v7d%anavar%d(indvar))
1025 height=
realdat(qccli%v7d%volanar(indana,indvar,indn),qccli%v7d%anavar%r(indvar))
1027 height=
realdat(qccli%v7d%volanai(indana,indvar,indn),qccli%v7d%anavar%i(indvar))
1029 height=
realdat(qccli%v7d%volanab(indana,indvar,indn),qccli%v7d%anavar%b(indvar))
1031 height=
realdat(qccli%v7d%volanac(indana,indvar,indn),qccli%v7d%anavar%c(indvar))
1037 if (
c_e(height))
exit
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
1076 if (indbattrinv > 0)
then
1078 (indana,indtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrinv)))
then
1082 call l4f_log (l4f_debug,
"qccli: skip station for a preceding invalidated flag")
1088 nintime=qccli%v7d%time(indtime)+timedelta_new(minute=30)
1089 CALL getval(nintime, month=mese, hour=ora)
1091 time=cyclicdatetime_to_conventional(cyclicdatetime_new(month=mese, hour=ora))
1096 level=qccli%v7d%level(indlevel)
1098 call init(network,
"qcclima-perc")
1100 indcnetwork =
index(qccli%extreme%network , network)
1101 indctime =
index(qccli%extreme%time , time)
1102 indclevel =
index(qccli%extreme%level , level)
1103 indctimerange =
index(qccli%extreme%timerange , qccli%v7d%timerange(indtimerange))
1107 indcdativarr =
index(qccli%extreme%dativar%r, qccli%v7d%dativar%r(inddativarr))
1126 if (indctime <= 0 .or. indclevel <= 0 .or. indctimerange <= 0 .or. indcdativarr <= 0 &
1127 .or. indcnetwork <= 0 ) cycle
1129 datoqui = qccli%v7d%voldatir (indana ,indtime ,indlevel ,indtimerange ,inddativarr, indnetwork )
1131 if (
c_e(datoqui))
then
1145 if (
associated(qccli%extreme%voldatir))
then
1147 if (qccli%height2level)
then
1159 write(ident,
'("#",i2.2,2i3.3)')k,iarea,desc
1160 call init(ana,ident=ident,lat=latc,lon=lonc)
1161 indcana=
index(qccli%extreme%ana,ana)
1162 if (indcana > 0 )
then
1163 perc25=qccli%extreme%voldatir(indcana,indctime,indclevel,indctimerange,indcdativarr,indcnetwork)
1167 write(ident,
'("#",i2.2,2i3.3)')k,iarea,desc
1168 call init(ana,ident=ident,lat=latc,lon=lonc)
1169 indcana=
index(qccli%extreme%ana,ana)
1172 if (indcana > 0 )
then
1173 perc50=qccli%extreme%voldatir(indcana,indctime,indclevel,indctimerange,indcdativarr,indcnetwork)
1177 write(ident,
'("#",i2.2,2i3.3)')k,iarea,desc
1178 call init(ana,ident=ident,lat=latc,lon=lonc)
1179 indcana=
index(qccli%extreme%ana,ana)
1180 if (indcana > 0 )
then
1181 perc75=qccli%extreme%voldatir(indcana,indctime,indclevel,indctimerange,indcdativarr,indcnetwork)
1185 if ( .not.
c_e(perc25) .or. .not.
c_e(perc50) .or. .not.
c_e(perc75)) cycle
1191 extremequii=perc50 - (perc75 - perc25) *1.3 * 3.65
1192 extremequif=perc50 + (perc75 - perc25) *1.3 * 3.65
1195 call l4f_log (l4f_debug,
"qccli: gross error check "//
t2c(extremequii)//
">"//
t2c(datoqui)//
"<"//
t2c(extremequif))
1199 if ( datoqui <= extremequii .or. extremequif <= datoqui )
then
1204 call l4f_log (l4f_debug,
"qccli: gross error check flag set to bad")
1206 qccli%v7d%voldatiattrb(indana,indtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrout)=qcpar%gross_error
1208 if (
associated ( qccli%data_id_in))
then
1210 call l4f_log (l4f_debug,
"id: "//
t2c(&
1211 qccli%data_id_in(indana,indtime,indlevel,indtimerange,indnetwork)))
1213 qccli%data_id_out(indana,indtime,indlevel,indtimerange,indnetwork)=&
1214 qccli%data_id_in(indana,indtime,indlevel,indtimerange,indnetwork)
1218 else if (.not.
vdge(qccli%v7d%voldatiattrb(indana,indtime,indlevel,indtimerange,&
1219 inddativarr,indnetwork,indbattrout)))
then
1223 call l4f_log (l4f_warn,
"qccli: skip station for a preceding gross error check flagged bad")
1230 datoqui = (datoqui - perc50) / (perc75 - perc25) + base_value(qccli%v7d%dativar%r(inddativarr)%btable)
1234 call init(network,
"qcclima-ndi")
1236 level=qccli%v7d%level(indlevel)
1237 time=cyclicdatetime_to_conventional(cyclicdatetime_new(month=mese))
1239 indcnetwork =
index(qccli%clima%network , network)
1240 indctime =
index(qccli%clima%time , time)
1241 indclevel =
index(qccli%clima%level , level)
1242 indctimerange =
index(qccli%clima%timerange , qccli%v7d%timerange(indtimerange))
1246 indcdativarr =
index(qccli%clima%dativar%r, qccli%v7d%dativar%r(inddativarr))
1250 if (indctime <= 0 .or. indclevel <= 0 .or. indctimerange <= 0 .or. indcdativarr <= 0 &
1251 .or. indcnetwork <= 0 ) cycle
1255 do desc=1,
size(qccli%clima%ana)
1260 write(ident,
'("#",i2.2,2i3.3)')0,0,min(desc,
size(qccli%clima%ana)-1) *10
1261 call init(ana,ident=ident,lat=0d0,lon=0d0)
1262 indcana=
index(qccli%clima%ana,ana)
1263 if (indcana > 0 )
then
1264 climaquif=qccli%clima%voldatir(indcana,indctime,indclevel,indctimerange,indcdativarr,indcnetwork)
1268 write(ident,
'("#",i2.2,2i3.3)')0,0,(desc-1)*10
1269 call init(ana,ident=ident,lat=0d0,lon=0d0)
1270 indcana=
index(qccli%clima%ana,ana)
1273 if (indcana > 0 )
then
1274 climaquii=qccli%clima%voldatir(indcana,indctime,indclevel,indctimerange,indcdativarr,indcnetwork)
1280 if (
c_e(climaquii) .and.
c_e(climaquif ))
then
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
1296 if (
c_e(qccli%clima%voldatiattrb(indcana &
1297 ,indctime,indclevel,indctimerange,indcdativarr,indcnetwork,1)))
then
1302 qccli%v7d%voldatiattrb(indana,indtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrout)=&
1303 max(qccli%clima%voldatiattrb&
1304 (indcana,indctime,indclevel,indctimerange,indcdativarr,indcnetwork,1)&
1308 call l4f_log (l4f_debug,
"data ndi: "//
t2c(datoqui)//
"->"//&
1309 t2c(qccli%clima%voldatiattrb(indcana,indctime,indclevel,indctimerange,indcdativarr,indcnetwork,1))&
1310 //
" : "//
t2c(qccli%v7d%time(indtime)))
1311 call l4f_log (l4f_debug,
"limits: "//
t2c(indcana)//
":"//qccli%clima%ana(indcana)%ident//&
1312 " : "//
t2c(climaquii)//
" - "//
t2c(climaquif)//
" : "//
t2c(qccli%clima%time(indctime)))
1313 call l4f_log (l4f_debug,
"qccli: clima check "//
t2c(datoqui)//
" confidence: "//&
1314 t2c(qccli%v7d%voldatiattrb(indana,indtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrout))&
1315 //
" : "//
t2c(qccli%v7d%time(indtime)))
1319 if (
associated ( qccli%data_id_in))
then
1321 call l4f_log (l4f_debug,
"id: "//
t2c(&
1322 qccli%data_id_in(indana,indtime,indlevel,indtimerange,indnetwork)))
1324 qccli%data_id_out(indana,indtime,indlevel,indtimerange,indnetwork)=&
1325 qccli%data_id_in(indana,indtime,indlevel,indtimerange,indnetwork)
1350end subroutine quaconcli
1356subroutine cli_level(heigth,level)
1358real,
intent(in) :: heigth
1359TYPE(vol7d_level),
intent(out):: level
1365if (
c_e(heigth))
then
1366 i=firsttrue(cli_level1 <= heigth .and. heigth <= cli_level2 )
1369if (i >= 1 .and. i <= 10 )
then
1370 call init(level, 102,cli_level1(i)*1000,102,cli_level2(i)*1000)
1372 if (
c_e(i))
CALL l4f_log(l4f_debug,
"cli_level: strange level, heigth: "//
to_char(heigth))
1376end subroutine cli_level
1380subroutine cli_level_generate(level)
1382TYPE(vol7d_level),
intent(out):: level(:)
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)))
1392 call init(level(i), 102,cli_level1(i)*1000,102,cli_level2(i)*1000)
1395end subroutine cli_level_generate
1409integer function supermacroa(macroa)
1411integer,
intent(in) :: macroa
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
1426end function supermacroa
1429SUBROUTINE qc_compute_percentile(this, perc_vals,cyclicdt,presentperc, presentnumb)
1435real,
intent(in) :: perc_vals(:)
1437real,
optional :: presentperc
1438integer,
optional :: presentnumb
1441integer :: indana,indtime,indvar,indnetwork,indlevel ,indtimerange ,inddativarr,i,j,k,iana,narea
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))
1449logical,
allocatable :: mask(:,:,:),maskplus(:,:,:), maskarea(:)
1450integer,
allocatable :: area(:)
1452integer :: lpresentnumb
1454lpresentperc=optio_r(presentperc)
1455lpresentnumb=optio_i(presentnumb)
1457allocate (perc(
size(perc_vals)))
1460CALL init(this%extreme, time_definition=this%v7d%time_definition)
1462call init(var, btable=
"B01192")
1465indvar =
index(this%v7d%anavar, var, type=type)
1466indnetwork=min(1,
size(this%v7d%network))
1468if( indvar > 0 .and. indnetwork > 0 )
then
1471 areav=
integerdat(this%v7d%volanad(:,indvar,indnetwork),this%v7d%anavar%d(indvar))
1473 areav=
integerdat(this%v7d%volanar(:,indvar,indnetwork),this%v7d%anavar%r(indvar))
1475 areav=
integerdat(this%v7d%volanai(:,indvar,indnetwork),this%v7d%anavar%i(indvar))
1477 areav=
integerdat(this%v7d%volanab(:,indvar,indnetwork),this%v7d%anavar%b(indvar))
1479 areav=
integerdat(this%v7d%volanac(:,indvar,indnetwork),this%v7d%anavar%c(indvar))
1487allocate(maskarea(
size(this%v7d%ana)))
1488maskarea(:)= areav(:) /= imiss
1489narea=count_distinct(areav,maskarea)
1490allocate(area(narea))
1491area=pack_distinct(areav,narea,maskarea)
1493if (this%height2level)
then
1494 call vol7d_alloc(this%extreme,nana=narea*
size(perc_vals)*cli_nlevel)
1496 call vol7d_alloc(this%extreme,nana=narea*
size(perc_vals))
1499if (this%height2level)
then
1501 call init(var, btable=
"B07030")
1504 indvar =
index(this%v7d%anavar, var, type=type)
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")
1614call vol7d_alloc_vol(this%extreme,inivol=.true.)
1616allocate (mask(
size(this%v7d%ana),
size(this%v7d%time),
size(this%v7d%network)))
1620do inddativarr=1,
size(this%v7d%dativar%r)
1621 do indtimerange=1,
size(this%v7d%timerange)
1622 do indlevel=1,
size(this%v7d%level)
1628 mask = spread(spread((this%v7d%time == cyclicdt ),1,
size(this%v7d%ana)),3,
size(this%v7d%network))
1632 (mask .and.
c_e(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:)))))
1636 do j=1,
size(mask,1)
1637 if (areav(j) /= area(i)) mask(j,:,:) =.false.
1640 if (this%height2level)
then
1641 allocate (maskplus(
size(this%v7d%ana),
size(this%v7d%time),
size(this%v7d%network)))
1647 do iana=1,
size(mask,1)
1648 if (iclv(iana) /= k) maskplus(iana,:,:) =.false.
1649 if (iclv(iana) == k) maskplus(iana,:,:) = mask(iana,:,:)
1652 call sub_perc(maskplus)
1655 deallocate(maskplus)
1667deallocate (perc,mask,area)
1671subroutine sub_perc(mymask)
1673logical :: mymask(:,:,:)
1677 (
c_e(lpresentperc) .and. ((float(count &
1678 (mymask .and.
c_e(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:)))&
1680 float(count(mymask))) < lpresentperc)) &
1684 (
c_e(lpresentnumb) .and. (count &
1685 (mymask .and.
c_e(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:))) < lpresentnumb)&
1691 pack(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:), &
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)=&
1703end subroutine sub_perc
1706end SUBROUTINE qc_compute_percentile
1709SUBROUTINE qc_compute_normalizeddensityindex(this, perc_vals,cyclicdt,presentperc, presentnumb,data_normalized)
1715real,
intent(in) :: perc_vals(:)
1717real,
optional :: presentperc
1718integer,
optional :: presentnumb
1719logical,
optional,
intent(in) :: data_normalized
1721integer :: indana,indtime,indvar,indnetwork,indlevel ,indtimerange ,inddativarr, indattr
1722integer :: i,j,k,iana,narea
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
1732integer :: lpresentnumb
1735lnorm = optio_log(data_normalized)
1736lpresentperc=optio_r(presentperc)
1737lpresentnumb=optio_i(presentnumb)
1739allocate (ndi(
size(perc_vals)-1),limbins(
size(perc_vals)))
1741CALL init(this%clima, time_definition=this%v7d%time_definition)
1742call init(var, btable=
"B01192")
1744if (.NOT.(lnorm))
then
1747 indvar =
index(this%v7d%anavar, var, type=type)
1753 areav=
integerdat(this%v7d%volanad(:,indvar,indnetwork),this%v7d%anavar%d(indvar))
1755 areav=
integerdat(this%v7d%volanar(:,indvar,indnetwork),this%v7d%anavar%r(indvar))
1757 areav=
integerdat(this%v7d%volanai(:,indvar,indnetwork),this%v7d%anavar%i(indvar))
1759 areav=
integerdat(this%v7d%volanab(:,indvar,indnetwork),this%v7d%anavar%b(indvar))
1761 areav=
integerdat(this%v7d%volanac(:,indvar,indnetwork),this%v7d%anavar%c(indvar))
1767 allocate(maskarea(
size(this%v7d%ana)))
1768 maskarea(:)= areav(:) /= imiss
1769 narea=count_distinct(areav,maskarea)
1770 allocate(area(narea))
1771 area=pack_distinct(areav,narea,maskarea)
1772 deallocate(maskarea)
1774 if (this%height2level)
then
1775 call vol7d_alloc(this%clima,nana=narea*(
size(perc_vals)-1)*cli_nlevel)
1777 call vol7d_alloc(this%clima,nana=narea*(
size(perc_vals)-1))
1784 do j=1,
size(perc_vals)-1
1785 if (this%height2level)
then
1787 write(ident,
'("#",i2.2,2i3.3)')k,area(i),nint(perc_vals(j))
1788 call init(this%clima%ana((k-1)*(
size(perc_vals)-1)*narea + (j-1)*narea + i),ident=ident,lat=0d0,lon=0d0)
1792 write(ident,
'("#",i2.2,2i3.3)')k,area(i),nint(perc_vals(j))
1793 call init(this%clima%ana((j-1)*narea+i),ident=ident,lat=0d0,lon=0d0)
1802 if (this%height2level)
then
1804 call init(var, btable=
"B07030")
1807 indvar =
index(this%v7d%anavar, var, type=type)
1809 do k=1,
size(this%v7d%ana)
1813 do indnetwork=1,
size(this%v7d%network)
1815 if( indvar > 0 .and. indnetwork > 0 )
then
1818 height=
realdat(this%v7d%volanad(k,indvar,indnetwork),this%v7d%anavar%d(indvar))
1820 height=
realdat(this%v7d%volanar(k,indvar,indnetwork),this%v7d%anavar%r(indvar))
1822 height=
realdat(this%v7d%volanai(k,indvar,indnetwork),this%v7d%anavar%i(indvar))
1824 height=
realdat(this%v7d%volanab(k,indvar,indnetwork),this%v7d%anavar%b(indvar))
1826 height=
realdat(this%v7d%volanac(k,indvar,indnetwork),this%v7d%anavar%c(indvar))
1832 if (
c_e(height))
exit
1835 if (
c_e(height))
then
1836 iclv(k)=firsttrue(cli_level1 <= height .and. height <= cli_level2 )
1842 CALL l4f_log(l4f_debug,
'qc_compute_NormalizedDensityIndex height has value '//
t2c(height,
"Missing"))
1843 CALL l4f_log(l4f_debug,
'for k having number '//
t2c(k)//&
1844 ' iclv has value '//
t2c(iclv(k)))
1855 call vol7d_alloc(this%clima,nana=
size(perc_vals)-1)
1856 do j=1,
size(perc_vals)-1
1857 write(ident,
'("#",i2.2,2i3.3)')0,0,nint(perc_vals(j))
1858 call init(this%clima%ana(j),ident=ident,lat=0d0,lon=0d0)
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)
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")
1876this%clima%time(1)=cyclicdatetime_to_conventional(cyclicdt)
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")
1881call vol7d_alloc_vol(this%clima,inivol=.true.)
1883allocate (mask(
size(this%v7d%ana),
size(this%v7d%time),
size(this%v7d%network)))
1888do inddativarr=1,
size(this%v7d%dativar%r)
1889 do indtimerange=1,
size(this%v7d%timerange)
1890 do indlevel=1,
size(this%v7d%level)
1891 if (.NOT.(lnorm))
then
1897 mask = spread(spread((this%v7d%time == cyclicdt ),1,
size(this%v7d%ana)),3,
size(this%v7d%network))
1899 do j=1,
size(mask,1)
1900 if (areav(j) /= area(i)) mask(j,:,:) =.false.
1903 if (this%height2level)
then
1904 allocate (maskplus(
size(this%v7d%ana),
size(this%v7d%time),
size(this%v7d%network)))
1910 do iana=1,
size(mask,1)
1911 if (iclv(iana) /= k) maskplus(iana,:,:) =.false.
1912 if (iclv(iana) == k) maskplus(iana,:,:) = mask(iana,:,:)
1915 call sub_ndi(mymask=maskplus)
1918 deallocate(maskplus)
1923 call sub_ndi(mymask=mask)
1931 mask = spread(spread((this%v7d%time == cyclicdt ),1,
size(this%v7d%ana)),3,
size(this%v7d%network))
1935 call sub_ndi(mymask=mask)
1943deallocate (ndi,limbins,mask)
1944if (.NOT.(lnorm))
deallocate(area)
1948subroutine sub_ndi(mymask)
1950logical :: mymask(:,:,:)
1954 (
c_e(lpresentperc) .and. ((float(count &
1955 (mymask .and.
c_e(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:)))&
1957 float(count(mymask))) < lpresentperc)) &
1961 (
c_e(lpresentnumb) .and. (count &
1962 (mymask .and.
c_e(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:))) < lpresentnumb)&
1969call normalizeddensityindex (&
1970 pack(this%v7d%voldatir(:,:, indlevel, indtimerange, inddativarr,:), &
1972 perc_vals, ndi, limbins)
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)=&
1978 this%clima%voldatiattrr(indana, indtime, indlevel, indtimerange, inddativarr, indnetwork,indattr)=&
1982end subroutine sub_ndi
1985end SUBROUTINE qc_compute_normalizeddensityindex