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
125 type (
vol7d) :: clima
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)
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)
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)
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)
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)
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
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
Classi per la gestione delle coordinate temporali.
Oggetto principale per il controllo di qualità
Functions that return a trimmed CHARACTER representation of the input variable.
Compute forward coordinate transformation from geographical system to projected system.
Check data validity based on gross error check.
Operatore di valore assoluto di un intervallo.
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.
Classes for handling georeferenced sparse points in geographical corodinates.
Generic subroutine for checking OPTIONAL parameters.
Classe per la gestione di un volume completo di dati osservati.
Restituiscono il valore dell'oggetto in forma di stringa stampabile.
Restituiscono il valore dell'oggetto nella forma desiderata.
Definisce un oggetto contenente i volumi anagrafica e dati e tutti i descrittori delle loro dimension...
Controllo di qualità climatico.
Oggetto principale per il controllo di qualità
Copy an object, creating a fully new instance.
classe per import ed export di volumi da e in DB-All.e
Utilities and defines for quality control.
Class for expressing a relative time interval.
Controllo di qualità spaziale.
classe per la gestione del logging
Space utilities, derived from NCAR software.
Utilities for CHARACTER variables.
Utilities for managing files.
Emit log message for a category with specific priority.