libsim Versione 7.2.0

◆ quaconspa()

subroutine, public quaconspa ( type(qcspatype), intent(inout)  qcspa,
type(timedelta), intent(in)  timetollerance,
logical, intent(in), optional  noborder,
character (len=10), intent(in), optional  battrinv,
character (len=10), intent(in), optional  battrcli,
character (len=10), intent(in), optional  battrout,
logical, dimension(:), intent(in), optional  anamask,
logical, dimension(:), intent(in), optional  timemask,
logical, dimension(:), intent(in), optional  levelmask,
logical, dimension(:), intent(in), optional  timerangemask,
logical, dimension(:), intent(in), optional  varmask,
logical, dimension(:), intent(in), optional  networkmask 
)

Controllo di Qualità spaziale.

Questo è il vero e proprio controllo di qualità spaziale.

Parametri
[in,out]qcspaOggetto per il controllo di qualità
[in]timetollerancetime tollerance to compare nearest stations
[in]noborderExclude border from QC
[in]battrinvattributo invalidated in input
[in]battrcliattributo con la confidenza climatologica in input
[in]battroutattributo con la confidenza spaziale in output
[in]anamaskFiltro sulle anagrafiche
[in]timemaskFiltro sul tempo
[in]levelmaskFiltro sui livelli
[in]timerangemaskfiltro sui timerange
[in]varmaskFiltro sulle variabili
[in]networkmaskFiltro sui network

Definizione alla linea 663 del file modqcspa.F90.

665
666
667 do indnetwork=1,size(qcspa%v7d%network)
668 do indana=1,size(qcspa%v7d%ana)
669
670!!$ call l4f_log(L4F_INFO,"Index:"// t2c(indana)//t2c(indnetwork)//t2c(indlevel)//&
671!!$ t2c(indtimerange)//t2c(inddativarr)//t2c(indtime))
672
673 if (.not. anamaskl(indana).or. .not. levelmaskl(indlevel) .or. &
674 .not. timerangemaskl(indtimerange) .or. .not. varmaskl(inddativarr) .or. .not. networkmaskl(indnetwork)) cycle
675
676 datoqui = qcspa%v7d%voldatir (indana ,indtime ,indlevel ,indtimerange ,inddativarr, indnetwork )
677 if (.not. c_e(datoqui)) cycle
678
679 ! invalidated
680 if (indbattrinv > 0) then
681 if( invalidated(qcspa%v7d%voldatiattrb&
682 (indana,indtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrinv))) then
683 call l4f_category_log(qcspa%category,l4f_warn,&
684 "It's better to do a reform on ana to v7d after peeling, before spatial QC")
685 cycle
686 end if
687 end if
688
689 ! gross error check
690 if (indbattrcli > 0) then
691 if( .not. vdge(qcspa%v7d%voldatiattrb&
692 (indana,indtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrcli))) then
693 call l4f_category_log(qcspa%category,l4f_warn,&
694 "It's better to do a reform on ana to v7d after peeling, before spatial QC")
695 cycle
696 end if
697 end if
698
699
700
701 !!call init(anavar,"B07031" )
702 !call init(anavar,"B07030" )
703 !indanavar = 0
704 !if (associated (qcspa%v7d%anavar%r)) then
705 ! indanavar = index(qcspa%v7d%anavar%r, anavar)
706 !end if
707 !if (indanavar <= 0 )cycle
708 !altezza= qcspa%v7d%volanar(indana,indanavar,indnetwork)
709 ! call spa_level(altezza,level)
710
711 if (qcspa%operation == "run") then
712
713 indclevel = index(qcspa%clima%level , qcspa%v7d%level(indlevel))
714 indctimerange = index(qcspa%clima%timerange , qcspa%v7d%timerange(indtimerange))
715
716 ! attenzione attenzione TODO
717 ! se leggo da bufr il default è char e non reale
718 indcdativarr = index(qcspa%clima%dativar%r, qcspa%v7d%dativar%r(inddativarr))
719
720
721#ifdef DEBUG
722 call l4f_log(l4f_debug,"Index:"// to_char(indctime)//to_char(indclevel)//&
723 to_char(indctimerange)//to_char(indcdativarr)//to_char(indcnetwork))
724#endif
725 if ( indctime <= 0 .or. indclevel <= 0 .or. indctimerange <= 0 .or. indcdativarr <= 0 &
726 .or. indcnetwork <= 0 ) cycle
727 end if
728
729 if (optio_log(noborder) .and. any(indana == qcspa%tri%ipl(:3*qcspa%tri%nl:3))) cycle
730
731 ! ITROV e` il numero di triangoli in cui e` presente il dato
732 itrov=0
733 ! cicla per tutti i triangoli
734 DO it=1,qcspa%tri%NT
735 ! se la stazione considerata e` in prima posizione
736 ! memorizza gli altri due vertici
737 IF(qcspa%tri%IPT(3*it-2).EQ.indana)THEN
738 itrov=itrov+1
739 ivert(2*itrov)=qcspa%tri%IPT(3*it)
740 ivert(2*itrov-1)=qcspa%tri%IPT(3*it-1)
741 cycle
742 END IF
743 ! se la stazione considerata e` in seconda posizione
744 ! memorizza gli altri due vertici
745 IF(qcspa%tri%IPT(3*it-1).EQ.indana)THEN
746 itrov=itrov+1
747 ivert(2*itrov)=qcspa%tri%IPT(3*it)
748 ivert(2*itrov-1)=qcspa%tri%IPT(3*it-2)
749 cycle
750 END IF
751 ! se la stazione considerata e` in terza posizione
752 ! memorizza gli altri due vertici
753 IF(qcspa%tri%IPT(3*it).EQ.indana)THEN
754 itrov=itrov+1
755 ivert(2*itrov)=qcspa%tri%IPT(3*it-1)
756 ivert(2*itrov-1)=qcspa%tri%IPT(3*it-2)
757 cycle
758 END IF
759 END DO
760 ! ITROV ora diviene il numero di vertici nell'intorno
761 ! della stazione trovati
762 itrov=itrov*2
763
764 ! WRITE(*,*)'NUMERO VERTICI',ITROV
765 ! WRITE(*,*)'VERTICI TROVATI = ',IVERT
766
767 ! ordina i vettori dei vertici secondo valori decrescenti
768
769 call sort(ivert(:itrov))
770
771!!$ DO I=1,ITROV-1
772!!$ DO KK=I+1,ITROV
773!!$ IF(IVERT(I).LT.IVERT(KK))THEN
774!!$ IC=IVERT(KK)
775!!$ IVERT(KK)=IVERT(I)
776!!$ IVERT(I)=IC
777!!$ ENDIF
778!!$ END DO
779!!$ END DO
780 ! toglie i valori doppi dal vettore dei vertici
781 iv=1
782 DO kk=2,itrov
783 IF(ivert(iv).NE.ivert(kk))THEN
784 iv=iv+1
785 ivert(iv)=ivert(kk)
786 ENDIF
787 END DO
788 IF (iv.GT.itrov)iv=itrov
789
790 ! WRITE(*,*)'NUMERO VERTICI puliti',IV
791 ! WRITE(*,*)'VERTICI PULITI = ',IVERT
792
793 ! inizia il controllo sulla stazione testando i gradienti
794 ! WRITE(*,*)'STAZIONE ',INDANA
795 ipos=0
796 ineg=0
797 ivb=0
798 gradmin=huge(gradmin)
799 DO i=1, iv
800 !find the nearest data in time
801 datola = qcspa%v7d%voldatir (ivert(i) ,indtime ,indlevel ,indtimerange ,inddativarr, indnetwork )
802
803 ! invalidated
804 if (indbattrinv > 0) then
805 if( invalidated(qcspa%v7d%voldatiattrb&
806 (ivert(i),indtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrinv))) then
807 datola=rmiss
808 end if
809 end if
810
811 ! gross error check
812 if (indbattrcli > 0) then
813 if( .not. vdge(qcspa%v7d%voldatiattrb&
814 (ivert(i),indtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrcli))) then
815 datola=rmiss
816 end if
817 end if
818 !TODO
819 ! if we do not have data from the same network at the same time
820 ! here we search for the first data found (nearest in time) looping over the network
821 ! we do not have priority for network to take in account
822
823 deltato=timedelta_miss
824 do indnet=1, size(qcspa%v7d%network)
825 datila = qcspa%v7d%voldatir (ivert(i) ,: ,indlevel ,indtimerange ,inddativarr, indnet )
826 do iindtime=1,size(qcspa%v7d%time)
827 if (.not. c_e(datila(iindtime))) cycle
828 ! invalidated
829 if (indbattrinv > 0 ) then
830 if (invalidated(qcspa%v7d%voldatiattrb&
831 (ivert(i),iindtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrinv))) cycle
832 end if
833 ! gross error check
834 if (indbattrcli > 0 )then
835 if (.not. vdge(qcspa%v7d%voldatiattrb&
836 (ivert(i),iindtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrcli))) cycle
837 end if
838
839 if (iindtime < indtime) then
840 deltat=qcspa%v7d%time(indtime)-qcspa%v7d%time(iindtime)
841 else if (iindtime >= indtime) then
842 deltat=qcspa%v7d%time(iindtime)-qcspa%v7d%time(indtime)
843 end if
844
845 if ((deltat < deltato .or. .not. c_e(deltato)) .and. deltat <= timetollerance ) then
846 datola = datila(iindtime)
847 deltato = deltat
848 end if
849 end do
850 end do
851
852
853 IF(.NOT.c_e(datola)) cycle
854 ! distanza tra le due stazioni
855 dist = distanza(qcspa%co(indana),qcspa%co(ivert(i)))
856 IF (dist.EQ.0.)THEN
857 call l4f_category_log(qcspa%category,l4f_error,"distance from two station == 0.")
858 call raise_error()
859 END IF
860
861#ifdef DEBUG
862 call l4f_log (l4f_debug,"distanza: "//t2c(dist))
863#endif
864
865 dist=max(dist,distmin)
866 ! modifica 23/9/1998
867 ! se la distanza supera distscol, stazioni scorrelate - salta -
868 if (dist > distscol) cycle
869 ivb=ivb+1
870 ! valore del gradiente nella direzione delle due stazioni
871 grad=(datoqui-datola)/(dist)
872 IF (grad >= 0.d0) ipos=ipos+1 ! se il gradiente e` positivo incrementa il contatore di positivi
873 IF (grad <= 0.d0) ineg=ineg+1 ! se il gradiente e` negativo incrementa il contatore di negativi
874
875 gradmin=min(gradmin,abs(grad))
876
877 END DO
878
879#ifdef DEBUG
880 call l4f_log (l4f_debug,"ivb: "//t2c(ivb))
881#endif
882
883 IF(ivb < 3) cycle ! do nothing if valid gradients < 3
884
885 IF (ipos == ivb .or. ineg == ivb)THEN ! se tutti i gradienti sono dello stesso segno
886
887 gradmin=sign(gradmin,dble(ipos-ineg))
888
889 if (qcspa%operation == "gradient") then
890 write(grunit,*)gradmin
891 end if
892
893 ! we normalize gradmin or denormalize climaqui after
894 ! gradmin=gradmin*spa_a(ind) + spa_b(ind)
895
896
897#ifdef DEBUG
898 call l4f_log (l4f_debug,"gradmin: "//t2c(gradmin))
899#endif
900
901 flag=bmiss
902
903 !ATTENZIONE TODO : inddativarr È UNA GRANDE SEMPLIFICAZIONE NON VERA SE TIPI DI DATO DIVERSI !!!!
904 if (qcspa%operation == "run") then
905
906 do indcana=1,size(qcspa%clima%ana)
907 climaquii=(qcspa%clima%voldatir(indcana &
908 ,indctime,indclevel,indctimerange,indcdativarr,indcnetwork)&
909 - spa_b(ind))/spa_a(ind) ! denormalize
910
911 ! the last interval have climaquii == climaquif
912 ! the check will be the last extreme to +infinite
913 climaquif=(qcspa%clima%voldatir(min(indcana+1,size(qcspa%clima%ana)) &
914 ,indctime,indclevel,indctimerange,indcdativarr,indcnetwork)&
915 - spa_b(ind))/spa_a(ind) ! denormalize
916
917#ifdef DEBUG
918 call l4f_log (l4f_debug,"climaquii: "//t2c(climaquii))
919 call l4f_log (l4f_debug,"climaquif: "//t2c(climaquif))
920#endif
921
922 if ( c_e(climaquii) .and. c_e(climaquif )) then
923
924 if ( (gradmin >= climaquii .and. gradmin < climaquif) .or. &
925 (indcana == 1 .and. gradmin < climaquii) .or. &
926 (indcana == size(qcspa%clima%ana) .and. gradmin >= climaquif) ) then
927
928 flag=qcspa%clima%voldatiattrb(indcana,indctime,indclevel,indctimerange,indcdativarr,indcnetwork,1)
929
930 if ( associated ( qcspa%data_id_in)) then
931#ifdef DEBUG
932 call l4f_log (l4f_debug,"id: "//t2c(&
933 qcspa%data_id_in(indana,indtime,indlevel,indtimerange,indnetwork)))
934#endif
935 qcspa%data_id_out(indana,indtime,indlevel,indtimerange,indnetwork)=&
936 qcspa%data_id_in(indana,indtime,indlevel,indtimerange,indnetwork)
937 end if
938 end if
939 end if
940 end do
941#ifdef DEBUG
942 call l4f_log (l4f_info,"datoqui: "//t2c(datoqui))
943 call l4f_log (l4f_info,"flag qcspa: "//t2c(flag))
944#endif
945
946 end if
947 else
948 flag=100_int_b
949 end if
950 if (qcspa%operation == "run") then
951 !TODO controllare se flag = missing comporta rimozione della precedente flag; risposta: SI quando sarà chiusa https://github.com/ARPA-SIMC/dballe/issues/44
952 qcspa%v7d%voldatiattrb( indana, indtime, indlevel, indtimerange, inddativarr, indnetwork, indbattrout)=flag
953 end if
954 end do
955 end do
956
957 if (qcspa%operation == "gradient") then
958 close (unit=grunit)
959 end if
960
961 end do
962 end do
963 end do
964end do
965
966!!$print*,"risultato"
967!!$print *,qcspa%v7d%voldatiattrb(:,:,:,:,:,:,indbattrout)
968!!$print*,"fine risultato"
969
970return
971
972contains
973
974elemental double precision function distanza (co1,co2)
975type(xy), intent(in) :: co1,co2
976
977
978distanza = sqrt((co2%x-co1%x)**2 + (co2%y-co1%y)**2)
979
980end function distanza
981
982end subroutine quaconspa
983
984
985end module modqcspa
986
987
990
Index method.
Controllo di qualità spaziale.
Definition: modqcspa.F90:273

Generated with Doxygen.