106 character (len=255),
parameter:: subcategoryspa=
"QCspa" 108 integer,
parameter :: spa_nvar=1
109 CHARACTER(len=10) :: spa_btable(spa_nvar)=(/
"B12101"/)
111 real,
parameter :: spa_a(spa_nvar) = (/1.e5/)
113 real,
parameter :: spa_b(spa_nvar) = (/273.15/)
117 type(
vol7d),
pointer :: v7d => null()
118 integer,
pointer :: data_id_in(:,:,:,:,:) => null()
119 integer,
pointer :: data_id_out(:,:,:,:,:) => null()
122 type(xy),
pointer :: co(:) => null()
123 type(triangles) :: tri
126 character(len=20):: operation
134 module procedure qcspainit
139 module procedure qcspaalloc
144 module procedure qcspadelete
157 subroutine qcspainit(qcspa,v7d,var, timei, timef, coordmin, coordmax, data_id_in,extremepath,spatialpath, &
159 dsne,usere,passworde,&
160 dsnspa,userspa,passwordspa,&
162 height2level,operation,categoryappend)
164 type(qcspatype),
intent(in out) :: qcspa
165 type(
vol7d),
intent(in),
target:: v7d
166 character(len=*),
INTENT(in) :: var(:)
169 TYPE(geo_coord),
INTENT(inout),
optional :: coordmin,coordmax
171 TYPE(datetime),
INTENT(in),
optional :: timei, timef
172 integer,
intent(in),
optional,
target:: data_id_in(:,:,:,:,:)
173 character(len=*),
intent(in),
optional :: extremepath
174 character(len=*),
intent(in),
optional :: spatialpath
175 logical ,
intent(in),
optional :: height2level
176 character(len=*),
INTENT(in),
OPTIONAL :: categoryappend
177 character(len=*),
optional :: operation
181 character(len=*),
intent(in),
optional :: dsne
182 character(len=*),
intent(in),
optional :: usere
183 character(len=*),
intent(in),
optional :: passworde
184 character(len=*),
intent(in),
optional :: dsnspa
185 character(len=*),
intent(in),
optional :: userspa
186 character(len=*),
intent(in),
optional :: passwordspa
187 character(len=512) :: ldsnspa
188 character(len=512) :: luserspa
189 character(len=512) :: lpasswordspa
193 TYPE(vol7d_network) :: network
194 character(len=512) :: filepathspa
195 character(len=512) :: a_name
197 call l4f_launcher(a_name,a_name_append=trim(subcategoryspa)//
"."//trim(categoryappend))
198 qcspa%category=l4f_category_get(a_name)
202 nullify ( qcspa%data_id_in )
203 nullify ( qcspa%data_id_out )
209 qcspa%operation=optio_c(operation,20)
210 filepathspa=optio_c(spatialpath,512)
213 if (qcspa%operation /=
"gradient" .and. qcspa%operation /=
"run")
then 214 call l4f_category_log(qcspa%category,l4f_error,
"operation is wrong: "//qcspa%operation)
219 if (
present(data_id_in))
then 220 qcspa%data_id_in => data_id_in
224 call init(qcspa%qccli,v7d,var, timei, timef, data_id_in,&
225 macropath=cmiss, climapath=cmiss, extremepath=extremepath, &
227 dsncli=cmiss,dsnextreme=dsne,user=usere,password=passworde,&
229 height2level=height2level,categoryappend=categoryappend)
234 if (qcspa%operation ==
"run")
then 236 call init(network,
"qcspa-ndi")
239 call optio(dsnspa,ldsnspa)
240 call optio(userspa,luserspa)
241 call optio(passwordspa,lpasswordspa)
243 if (
c_e(filepathspa) .and. (
c_e(ldsnspa).or.
c_e(luserspa).or.
c_e(lpasswordspa)))
then 244 call l4f_category_log(qcspa%category,l4f_error,
"filepath defined together with dba options")
248 if (.not.
c_e(ldsnspa))
then 252 if (.not.
c_e(filepathspa))
then 253 filepathspa=get_package_filepath(
'qcspa-ndi.v7d', filetype_data)
256 if (
c_e(filepathspa))
then 258 select case (trim(lowercase(suffixname(filepathspa))))
262 call import(qcspa%clima,filename=filepathspa,unit=iuni)
267 call init(v7d_dballespa,file=.true.,filename=filepathspa,categoryappend=trim(a_name)//
".clima")
268 call import(v7d_dballespa,var=var, &
269 varkind=(/(
"r",i=1,
size(var))/),attr=(/
"*B33209"/),attrkind=(/
"b"/),network=network)
270 call copy(v7d_dballespa%vol7d,qcspa%clima)
271 call delete(v7d_dballespa)
276 "file type not supported (use .v7d or .bufr suffix only): "//trim(filepathspa))
281 call l4f_category_log(qcspa%category,l4f_warn,
"spatial clima volume not iniziatized: spatial QC will not be possible")
282 call init(qcspa%clima)
283 call raise_fatal_error()
290 call init(v7d_dballespa,dsn=ldsnspa,user=luserspa,password=lpasswordspa,write=.false.,&
291 file=.false.,categoryappend=trim(a_name)//
".spa")
293 call import(v7d_dballespa,var=var, &
294 varkind=(/(
"r",i=1,
size(var))/),attr=(/
"*B33209"/),attrkind=(/
"b"/),network=network)
295 call copy(v7d_dballespa%vol7d,qcspa%clima)
296 call delete(v7d_dballespa)
304 end subroutine qcspainit
307 subroutine qcspatri(qcspa,proj_type, lov, zone, xoff, yoff, &
308 longitude_south_pole, latitude_south_pole, angle_rotation, &
309 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
310 latin1, latin2, lad, projection_center_flag, &
311 ellips_smaj_axis, ellips_flatt, ellips_type)
314 CHARACTER(len=*),
INTENT(in),
OPTIONAL :: proj_type
315 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: lov
316 INTEGER,
INTENT(in),
OPTIONAL :: zone
317 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: xoff
318 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: yoff
319 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: longitude_south_pole
320 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latitude_south_pole
321 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: angle_rotation
322 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: longitude_stretch_pole
323 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latitude_stretch_pole
324 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: stretch_factor
325 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latin1
326 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latin2
327 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: lad
328 INTEGER,
INTENT(in),
OPTIONAL :: projection_center_flag
329 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: ellips_smaj_axis
330 DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: ellips_flatt
331 INTEGER,
INTENT(in),
OPTIONAL :: ellips_type
335 TYPE(geo_proj) :: geoproj
336 REAL(kind=fp_geo) :: lat(size(qcspa%v7d%ana)),lon(size(qcspa%v7d%ana))
337 CHARACTER(len=80) :: proj_type_l
338 double precision :: lov_l, latin1_l,latin2_l
339 integer :: projection_center_flag_l
341 proj_type_l = optio_c(proj_type,80)
344 latin1_l = optio_d(latin1)
345 latin2_l = optio_d(latin2)
346 projection_center_flag_l=optio_l(projection_center_flag)
348 if (.not.
c_e(proj_type_l))
then 349 proj_type_l =
"lambert" 353 projection_center_flag_l=1
356 geoproj = geo_proj_new(proj_type_l, lov_l, zone, xoff, yoff, &
357 longitude_south_pole, latitude_south_pole, angle_rotation, &
358 longitude_stretch_pole, latitude_stretch_pole, stretch_factor, &
359 latin1_l, latin2_l, lad, projection_center_flag_l, &
360 ellips_smaj_axis, ellips_flatt, ellips_type)
362 call getval(qcspa%v7d%ana%coord, lon, lat)
366 call proj(geoproj,lon,lat,qcspa%co%x,qcspa%co%y)
371 status = triangles_compute(qcspa%co,qcspa%tri)
375 if (status /= 0)
then 380 end subroutine qcspatri
384 subroutine qcspaalloc(qcspa)
387 type(qcspatype),
intent(in out) :: qcspa
393 call qcspadealloc(qcspa)
408 if (
associated(qcspa%data_id_in))
then 409 sh=shape(qcspa%data_id_in)
410 allocate (qcspa%data_id_out(sh(1),sh(2),sh(3),sh(4),sh(5)),stat=istatt)
413 call raise_error(
"allocate error")
415 qcspa%data_id_out=imiss
419 if (
associated(qcspa%v7d%ana))
then 420 qcspa%ndp=
size(qcspa%v7d%ana)
421 qcspa%tri = triangles_new(qcspa%ndp)
422 allocate(qcspa%co(qcspa%ndp))
425 end subroutine qcspaalloc
429 subroutine qcspadealloc(qcspa)
432 type(qcspatype),
intent(in out) :: qcspa
442 if (
associated(qcspa%data_id_out))
then 443 deallocate (qcspa%data_id_out)
444 nullify (qcspa%data_id_out)
447 if (
associated(qcspa%co))
deallocate(qcspa%co)
449 end subroutine qcspadealloc
455 subroutine qcspadelete(qcspa)
457 type(qcspatype),
intent(in out) :: qcspa
459 call qcspadealloc(qcspa)
466 call l4f_category_delete(qcspa%category)
469 end subroutine qcspadelete
475 SUBROUTINE quaconspa (qcspa,timetollerance,noborder,battrinv,battrcli,battrout,&
476 anamask,timemask,levelmask,timerangemask,varmask,networkmask)
479 type(qcspatype),
intent(in out) :: qcspa
480 type(timedelta),
intent(in) :: timetollerance
481 logical,
intent(in),
optional :: noborder
482 character (len=10) ,
intent(in),
optional :: battrinv
483 character (len=10) ,
intent(in),
optional :: battrcli
484 character (len=10) ,
intent(in),
optional :: battrout
485 logical ,
intent(in),
optional :: anamask(:)
486 logical ,
intent(in),
optional :: timemask(:)
487 logical ,
intent(in),
optional :: levelmask(:)
488 logical ,
intent(in),
optional :: timerangemask(:)
489 logical ,
intent(in),
optional :: varmask(:)
490 logical ,
intent(in),
optional :: networkmask(:)
494 integer :: indbattrinv,indbattrcli,indbattrout
495 logical :: anamaskl(size(qcspa%v7d%ana)), timemaskl(size(qcspa%v7d%time)), levelmaskl(size(qcspa%v7d%level)), &
496 timerangemaskl(size(qcspa%v7d%timerange)), varmaskl(size(qcspa%v7d%dativar%r)), networkmaskl(size(qcspa%v7d%network))
498 integer :: indana , indtime ,indlevel ,indtimerange ,inddativarr, indnetwork,indnet
499 integer :: indcana , indctime ,indclevel ,indctimerange ,indcdativarr, indcnetwork
500 real :: datoqui,datola,datila(size(qcspa%v7d%time)),climaquii, climaquif
504 TYPE(datetime) :: time
506 TYPE(vol7d_network):: network
507 type(timedelta) :: deltato,deltat
509 integer :: ivert(50),i,ipos,ineg,it,itrov,iv,ivb,kk,iindtime,grunit
510 double precision :: distmin=1000.d0,distscol=100000.d0
511 double precision :: dist,grad,gradmin
512 integer (kind=int_b) :: flag
514 character(len=512) :: filename
520 if (
size(qcspa%v7d%ana) < 3 )
then 521 call l4f_category_log(qcspa%category,l4f_warn,
"number of station < 3; do nothing")
526 if (
present(battrinv))
then 527 indbattrinv = index_c(qcspa%v7d%datiattr%b(:)%btable, battrinv)
529 indbattrinv = index_c(qcspa%v7d%datiattr%b(:)%btable, qcattrvarsbtables(1))
532 if (
present(battrcli))
then 533 indbattrcli = index_c(qcspa%v7d%datiattr%b(:)%btable, battrcli)
535 indbattrcli = index_c(qcspa%v7d%datiattr%b(:)%btable, qcattrvarsbtables(2))
538 if (
present(battrout))
then 539 indbattrout = index_c(qcspa%v7d%datiattr%b(:)%btable, battrout)
541 indbattrout = index_c(qcspa%v7d%datiattr%b(:)%btable, qcattrvarsbtables(4))
547 if (indbattrout <= 0 )
then 549 call l4f_category_log(qcspa%category,l4f_error,
"error finding attribute index for output")
567 if(
present(anamask))
then 572 if(
present(timemask))
then 577 if(
present(levelmask))
then 578 levelmaskl = levelmask
582 if(
present(timerangemask))
then 583 timerangemaskl = timerangemask
585 timerangemaskl = .true.
587 if(
present(varmask))
then 592 if(
present(networkmask))
then 593 networkmaskl = networkmask
595 networkmaskl = .true.
599 qcspa%v7d%voldatiattrb(:,:,:,:,:,:,indbattrout)=ibmiss
604 call vol7d_normalize_data(qcspa%qccli)
616 time=cyclicdatetime_to_conventional(cyclicdatetime_new(chardate=
"/////////"))
620 if (qcspa%operation ==
"run")
then 621 call init(network,
"qcspa-ndi")
623 indcnetwork =
index(qcspa%clima%network , network)
624 indctime =
index(qcspa%clima%time , time)
627 do indtime=1,
size(qcspa%v7d%time)
628 if (.not.timemaskl(indtime)) cycle
630 "Check time:"//
t2c(qcspa%v7d%time(indtime)) )
632 do indlevel=1,
size(qcspa%v7d%level)
633 do indtimerange=1,
size(qcspa%v7d%timerange)
634 do inddativarr=1,
size(qcspa%v7d%dativar%r)
636 ind=index_c(spa_btable,qcspa%v7d%dativar%r(inddativarr)%btable)
638 if (qcspa%operation ==
"gradient")
then 641 filename=trim(
to_char(qcspa%v7d%level(indlevel)))//&
642 "_"//trim(
to_char(qcspa%v7d%timerange(indtimerange)))//&
643 "_"//trim(qcspa%v7d%dativar%r(inddativarr)%btable)//&
646 call l4f_category_log(qcspa%category,l4f_info,
"try to open gradient file; filename below")
649 inquire(file=filename, exist=exist)
652 if (grunit /= -1)
then 654 open (grunit, file=filename ,status=
'UNKNOWN', form=
'FORMATTED',position=
'APPEND')
657 if (.not. exist)
then 658 call l4f_category_log(qcspa%category,l4f_info,
"write header in gradient file")
660 qcspa%v7d%level(indlevel), &
661 qcspa%v7d%timerange(indtimerange), &
662 qcspa%v7d%dativar%r(inddativarr)
667 do indnetwork=1,
size(qcspa%v7d%network)
668 do indana=1,
size(qcspa%v7d%ana)
673 if (.not. anamaskl(indana).or. .not. levelmaskl(indlevel) .or. &
674 .not. timerangemaskl(indtimerange) .or. .not. varmaskl(inddativarr) .or. .not. networkmaskl(indnetwork)) cycle
676 datoqui = qcspa%v7d%voldatir (indana ,indtime ,indlevel ,indtimerange ,inddativarr, indnetwork )
677 if (.not.
c_e(datoqui)) cycle
680 if (indbattrinv > 0)
then 682 (indana,indtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrinv)))
then 684 "It's better to do a reform on ana to v7d after peeling, before spatial QC")
690 if (indbattrcli > 0)
then 691 if( .not.
vdge(qcspa%v7d%voldatiattrb&
692 (indana,indtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrcli)))
then 694 "It's better to do a reform on ana to v7d after peeling, before spatial QC")
711 if (qcspa%operation ==
"run")
then 713 indclevel =
index(qcspa%clima%level , qcspa%v7d%level(indlevel))
714 indctimerange =
index(qcspa%clima%timerange , qcspa%v7d%timerange(indtimerange))
718 indcdativarr =
index(qcspa%clima%dativar%r, qcspa%v7d%dativar%r(inddativarr))
722 call l4f_log(l4f_debug,
"Index:"//
to_char(indctime)//
to_char(indclevel)//&
725 if ( indctime <= 0 .or. indclevel <= 0 .or. indctimerange <= 0 .or. indcdativarr <= 0 &
726 .or. indcnetwork <= 0 ) cycle
729 if (optio_log(noborder) .and. any(indana == qcspa%tri%ipl(:3*qcspa%tri%nl:3))) cycle
737 IF(qcspa%tri%IPT(3*it-2).EQ.indana)
THEN 739 ivert(2*itrov)=qcspa%tri%IPT(3*it)
740 ivert(2*itrov-1)=qcspa%tri%IPT(3*it-1)
745 IF(qcspa%tri%IPT(3*it-1).EQ.indana)
THEN 747 ivert(2*itrov)=qcspa%tri%IPT(3*it)
748 ivert(2*itrov-1)=qcspa%tri%IPT(3*it-2)
753 IF(qcspa%tri%IPT(3*it).EQ.indana)
THEN 755 ivert(2*itrov)=qcspa%tri%IPT(3*it-1)
756 ivert(2*itrov-1)=qcspa%tri%IPT(3*it-2)
769 call sort(ivert(:itrov))
783 IF(ivert(iv).NE.ivert(kk))
THEN 788 IF (iv.GT.itrov)iv=itrov
798 gradmin=huge(gradmin)
801 datola = qcspa%v7d%voldatir (ivert(i) ,indtime ,indlevel ,indtimerange ,inddativarr, indnetwork )
804 if (indbattrinv > 0)
then 806 (ivert(i),indtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrinv)))
then 812 if (indbattrcli > 0)
then 813 if( .not.
vdge(qcspa%v7d%voldatiattrb&
814 (ivert(i),indtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrcli)))
then 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
829 if (indbattrinv > 0 )
then 831 (ivert(i),iindtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrinv))) cycle
834 if (indbattrcli > 0 )
then 835 if (.not.
vdge(qcspa%v7d%voldatiattrb&
836 (ivert(i),iindtime,indlevel,indtimerange,inddativarr,indnetwork,indbattrcli))) cycle
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)
845 if ((deltat < deltato .or. .not.
c_e(deltato)) .and. deltat <= timetollerance )
then 846 datola = datila(iindtime)
853 IF(.NOT.
c_e(datola)) cycle
855 dist = distanza(qcspa%co(indana),qcspa%co(ivert(i)))
857 call l4f_category_log(qcspa%category,l4f_error,
"distance from two station == 0.")
862 call l4f_log (l4f_debug,
"distanza: "//
t2c(dist))
865 dist=max(dist,distmin)
868 if (dist > distscol) cycle
871 grad=(datoqui-datola)/(dist)
872 IF (grad >= 0.d0) ipos=ipos+1
873 IF (grad <= 0.d0) ineg=ineg+1
875 gradmin=min(gradmin,
abs(grad))
880 call l4f_log (l4f_debug,
"ivb: "//
t2c(ivb))
885 IF (ipos == ivb .or. ineg == ivb)
THEN 887 gradmin=sign(gradmin,dble(ipos-ineg))
889 if (qcspa%operation ==
"gradient")
then 890 write(grunit,*)gradmin
898 call l4f_log (l4f_debug,
"gradmin: "//
t2c(gradmin))
904 if (qcspa%operation ==
"run")
then 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)
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)
918 call l4f_log (l4f_debug,
"climaquii: "//
t2c(climaquii))
919 call l4f_log (l4f_debug,
"climaquif: "//
t2c(climaquif))
922 if (
c_e(climaquii) .and.
c_e(climaquif ))
then 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 928 flag=qcspa%clima%voldatiattrb(indcana,indctime,indclevel,indctimerange,indcdativarr,indcnetwork,1)
930 if (
associated ( qcspa%data_id_in))
then 932 call l4f_log (l4f_debug,
"id: "//
t2c(&
933 qcspa%data_id_in(indana,indtime,indlevel,indtimerange,indnetwork)))
935 qcspa%data_id_out(indana,indtime,indlevel,indtimerange,indnetwork)=&
936 qcspa%data_id_in(indana,indtime,indlevel,indtimerange,indnetwork)
942 call l4f_log (l4f_info,
"datoqui: "//
t2c(datoqui))
943 call l4f_log (l4f_info,
"flag qcspa: "//
t2c(flag))
950 if (qcspa%operation ==
"run")
then 952 qcspa%v7d%voldatiattrb( indana, indtime, indlevel, indtimerange, inddativarr, indnetwork, indbattrout)=flag
957 if (qcspa%operation ==
"gradient")
then 974 elemental double precision function distanza (co1,co2)
975 type(xy),
intent(in) :: co1,co2
978 distanza = sqrt((co2%x-co1%x)**2 + (co2%y-co1%y)**2)
980 end function distanza
982 end subroutine quaconspa
Oggetto principale per il controllo di qualità
Compute forward coordinate transformation from geographical system to projected system.
Set of functions that return a trimmed CHARACTER representation of the input variable.
Check data validity based on gross error check.
Operatore di valore assoluto di un intervallo.
Represent geo_coord object in a pretty string.
Import one or more geo_coordvect objects from a plain text file or for a file in ESRI/Shapefile forma...
Oggetto per import ed export da DB-All.e.
Classi per la gestione delle coordinate temporali.
classe per import ed export di volumi da e in DB-All.e
Utilities for managing files.
Classes for handling georeferenced sparse points in geographical corodinates.
Controllo di qualità climatico.
Space utilities, derived from NCAR software.
Classe per la gestione di un volume completo di dati osservati.
Utilities and defines for quality control.
Definisce un oggetto contenente i volumi anagrafica e dati e tutti i descrittori delle loro dimension...
Oggetto principale per il controllo di qualità
Copy an object, creating a fully new instance.
classe per la gestione del logging
Utilities for CHARACTER variables.
Methods for returning the value of object members.
Controllo di qualità spaziale.
Emit log message for a category with specific priority.