106character (len=255),
parameter:: subcategoryspa=
"QCspa"
108integer,
parameter :: spa_nvar=1
109CHARACTER(len=10) :: spa_btable(spa_nvar)=(/
"B12101"/)
111real,
parameter :: spa_a(spa_nvar) = (/1.e5/)
113real,
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
124 type (qcclitype) :: qccli
125 type (vol7d) :: clima
126 character(len=20):: operation
134 module procedure qcspainit
139 module procedure qcspaalloc
144 module procedure qcspadelete
157subroutine 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)
164type(qcspatype),
intent(in out) :: qcspa
165type (vol7d),
intent(in),
target:: v7d
166character(len=*),
INTENT(in) :: var(:)
169TYPE(geo_coord),
INTENT(inout),
optional :: coordmin,coordmax
171TYPE(datetime),
INTENT(in),
optional :: timei, timef
172integer,
intent(in),
optional,
target:: data_id_in(:,:,:,:,:)
173character(len=*),
intent(in),
optional :: extremepath
174character(len=*),
intent(in),
optional :: spatialpath
175logical ,
intent(in),
optional :: height2level
176character(len=*),
INTENT(in),
OPTIONAL :: categoryappend
177character(len=*),
optional :: operation
180type (vol7d_dballe) :: v7d_dballespa
181character(len=*),
intent(in),
optional :: dsne
182character(len=*),
intent(in),
optional :: usere
183character(len=*),
intent(in),
optional :: passworde
184character(len=*),
intent(in),
optional :: dsnspa
185character(len=*),
intent(in),
optional :: userspa
186character(len=*),
intent(in),
optional :: passwordspa
187character(len=512) :: ldsnspa
188character(len=512) :: luserspa
189character(len=512) :: lpasswordspa
193TYPE(vol7d_network) :: network
194character(len=512) :: filepathspa
195character(len=512) :: a_name
197call l4f_launcher(a_name,a_name_append=trim(subcategoryspa)//
"."//trim(categoryappend))
198qcspa%category=l4f_category_get(a_name)
202nullify ( qcspa%data_id_in )
203nullify ( qcspa%data_id_out )
209qcspa%operation=optio_c(operation,20)
210filepathspa=optio_c(spatialpath,512)
213if (qcspa%operation /=
"gradient" .and. qcspa%operation /=
"run")
then
214 call l4f_category_log(qcspa%category,l4f_error,
"operation is wrong: "//qcspa%operation)
219if (
present(data_id_in))
then
220 qcspa%data_id_in => data_id_in
224call 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)
234if (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)
304end subroutine qcspainit
307subroutine 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)
314CHARACTER(len=*),
INTENT(in),
OPTIONAL :: proj_type
315DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: lov
316INTEGER,
INTENT(in),
OPTIONAL :: zone
317DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: xoff
318DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: yoff
319DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: longitude_south_pole
320DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latitude_south_pole
321DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: angle_rotation
322DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: longitude_stretch_pole
323DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latitude_stretch_pole
324DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: stretch_factor
325DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latin1
326DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: latin2
327DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: lad
328INTEGER,
INTENT(in),
OPTIONAL :: projection_center_flag
329DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: ellips_smaj_axis
330DOUBLE PRECISION,
INTENT(in),
OPTIONAL :: ellips_flatt
331INTEGER,
INTENT(in),
OPTIONAL :: ellips_type
335TYPE(geo_proj) :: geoproj
336REAL(kind=fp_geo) :: lat(
size(qcspa%v7d%ana)),lon(
size(qcspa%v7d%ana))
337CHARACTER(len=80) :: proj_type_l
338double precision :: lov_l, latin1_l,latin2_l
339integer :: projection_center_flag_l
341proj_type_l = optio_c(proj_type,80)
344latin1_l = optio_d(latin1)
345latin2_l = optio_d(latin2)
346projection_center_flag_l=optio_l(projection_center_flag)
348if (.not.
c_e(proj_type_l))
then
349 proj_type_l =
"lambert"
353 projection_center_flag_l=1
356geoproj = 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)
362call getval(qcspa%v7d%ana%coord, lon, lat)
366call proj(geoproj,lon,lat,qcspa%co%x,qcspa%co%y)
371status = triangles_compute(qcspa%co,qcspa%tri)
380end subroutine qcspatri
384subroutine qcspaalloc(qcspa)
393call qcspadealloc(qcspa)
408if (
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
419if (
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))
425end subroutine qcspaalloc
429subroutine qcspadealloc(qcspa)
442if (
associated(qcspa%data_id_out))
then
443 deallocate (qcspa%data_id_out)
444 nullify (qcspa%data_id_out)
447if (
associated(qcspa%co))
deallocate(qcspa%co)
449end subroutine qcspadealloc
455subroutine qcspadelete(qcspa)
459call qcspadealloc(qcspa)
466call l4f_category_delete(qcspa%category)
469end subroutine qcspadelete
475SUBROUTINE quaconspa (qcspa,timetollerance,noborder,battrinv,battrcli,battrout,&
476 anamask,timemask,levelmask,timerangemask,varmask,networkmask)
480type(
timedelta),
intent(in) :: timetollerance
481logical,
intent(in),
optional :: noborder
482character (len=10) ,
intent(in),
optional :: battrinv
483character (len=10) ,
intent(in),
optional :: battrcli
484character (len=10) ,
intent(in),
optional :: battrout
485logical ,
intent(in),
optional :: anamask(:)
486logical ,
intent(in),
optional :: timemask(:)
487logical ,
intent(in),
optional :: levelmask(:)
488logical ,
intent(in),
optional :: timerangemask(:)
489logical ,
intent(in),
optional :: varmask(:)
490logical ,
intent(in),
optional :: networkmask(:)
494integer :: indbattrinv,indbattrcli,indbattrout
495logical :: 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))
498integer :: indana , indtime ,indlevel ,indtimerange ,inddativarr, indnetwork,indnet
499integer :: indcana , indctime ,indclevel ,indctimerange ,indcdativarr, indcnetwork
500real :: datoqui,datola,datila(size(qcspa%v7d%time)),climaquii, climaquif
506TYPE(vol7d_network):: network
509integer :: ivert(50),i,ipos,ineg,it,itrov,iv,ivb,kk,iindtime,grunit
510double precision :: distmin=1000.d0,distscol=100000.d0
511double precision :: dist,grad,gradmin
512integer (kind=int_b) :: flag
514character(len=512) :: filename
520if (
size(qcspa%v7d%ana) < 3 )
then
521 call l4f_category_log(qcspa%category,l4f_warn,
"number of station < 3; do nothing")
526if (
present(battrinv))
then
527 indbattrinv = index_c(qcspa%v7d%datiattr%b(:)%btable, battrinv)
529 indbattrinv = index_c(qcspa%v7d%datiattr%b(:)%btable, qcattrvarsbtables(1))
532if (
present(battrcli))
then
533 indbattrcli = index_c(qcspa%v7d%datiattr%b(:)%btable, battrcli)
535 indbattrcli = index_c(qcspa%v7d%datiattr%b(:)%btable, qcattrvarsbtables(2))
538if (
present(battrout))
then
539 indbattrout = index_c(qcspa%v7d%datiattr%b(:)%btable, battrout)
541 indbattrout = index_c(qcspa%v7d%datiattr%b(:)%btable, qcattrvarsbtables(4))
547if (indbattrout <= 0 )
then
549 call l4f_category_log(qcspa%category,l4f_error,
"error finding attribute index for output")
567if(
present(anamask))
then
572if(
present(timemask))
then
577if(
present(levelmask))
then
578 levelmaskl = levelmask
582if(
present(timerangemask))
then
583 timerangemaskl = timerangemask
585 timerangemaskl = .true.
587if(
present(varmask))
then
592if(
present(networkmask))
then
593 networkmaskl = networkmask
595 networkmaskl = .true.
599qcspa%v7d%voldatiattrb(:,:,:,:,:,:,indbattrout)=ibmiss
604call vol7d_normalize_data(qcspa%qccli)
616time=cyclicdatetime_to_conventional(cyclicdatetime_new(chardate=
"/////////"))
620if (qcspa%operation ==
"run")
then
621 call init(network,
"qcspa-ndi")
623 indcnetwork =
index(qcspa%clima%network , network)
624 indctime =
index(qcspa%clima%time , time)
627do 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
974elemental double precision function distanza (co1,co2)
975type(xy),
intent(in) :: co1,co2
978distanza = sqrt((co2%x-co1%x)**2 + (co2%y-co1%y)**2)
982end subroutine quaconspa
Set of functions that return a trimmed CHARACTER representation of the input variable.
Set of functions that return a CHARACTER representation of the input variable.
Operatore di valore assoluto di un intervallo.
Restituiscono il valore dell'oggetto nella forma desiderata.
Import one or more geo_coordvect objects from a plain text file or for a file in ESRI/Shapefile forma...
Copy an object, creating a fully new instance.
Compute forward coordinate transformation from geographical system to projected system.
Emit log message for a category with specific priority.
Check data validity based on gross error check.
Utilities for CHARACTER variables.
Classi per la gestione delle coordinate temporali.
Utilities for managing files.
Classes for handling georeferenced sparse points in geographical corodinates.
classe per la gestione del logging
Utilities and defines for quality control.
Controllo di qualità climatico.
Controllo di qualità spaziale.
Space utilities, derived from NCAR software.
Classe per la gestione di un volume completo di dati osservati.
classe per import ed export di volumi da e in DB-All.e
Class for expressing an absolute time value.
Class for expressing a relative time interval.
Derived type defining an isolated georeferenced point on Earth in polar geographical coordinates.
Oggetto principale per il controllo di qualità
Oggetto principale per il controllo di qualità
Definisce un oggetto contenente i volumi anagrafica e dati e tutti i descrittori delle loro dimension...