libsim  Versione 7.2.6

◆ quaconspa()

subroutine, public modqcspa::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 673 del file modqcspa.F90.

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
964 end do
965 
966 !!$print*,"risultato"
967 !!$print *,qcspa%v7d%voldatiattrb(:,:,:,:,:,:,indbattrout)
968 !!$print*,"fine risultato"
969 
970 return
971 
972 contains
973 
974 elemental double precision function distanza (co1,co2)
975 type(xy), intent(in) :: co1,co2
976 
977 
978 distanza = sqrt((co2%x-co1%x)**2 + (co2%y-co1%y)**2)
979 
980 end function distanza
981 
982 end subroutine quaconspa
983 
984 
985 end module modqcspa
986 
987 
990 
Index method.
Controllo di qualità spaziale.
Definition: modqcspa.F90:283

Generated with Doxygen.