libsim Versione 7.2.1
modqcspa.F90
1! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
2! authors:
3! Davide Cesari <dcesari@arpa.emr.it>
4! Paolo Patruno <ppatruno@arpa.emr.it>
5
6! This program is free software; you can redistribute it and/or
7! modify it under the terms of the GNU General Public License as
8! published by the Free Software Foundation; either version 2 of
9! the License, or (at your option) any later version.
10
11! This program is distributed in the hope that it will be useful,
12! but WITHOUT ANY WARRANTY; without even the implied warranty of
13! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14! GNU General Public License for more details.
15
16! You should have received a copy of the GNU General Public License
17! along with this program. If not, see <http://www.gnu.org/licenses/>.
18
19#include "config.h"
20
77
83
84
85module modqcspa
86
88use geo_proj_class
91use modqc
92use modqccli
97!use array_utilities
98!use io_units
99#ifdef HAVE_DBALLE
101#endif
102
103implicit none
104
105
106character (len=255),parameter:: subcategoryspa="QCspa"
107
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/)
114
116type :: qcspatype
117 type (vol7d),pointer :: v7d => null()
118 integer,pointer :: data_id_in(:,:,:,:,:) => null()
119 integer,pointer :: data_id_out(:,:,:,:,:) => null()
120 integer :: category
121 integer :: ndp
122 type(xy),pointer :: co(:) => null()
123 type (triangles) :: tri
124 type (qcclitype) :: qccli
125 type (vol7d) :: clima
126 character(len=20):: operation
127 !logical :: writeheader !< have to write header in gradient files
128 !integer :: grunit !< unit used internally to write gradient
129end type qcspatype
130
131
133interface init
134 module procedure qcspainit
135end interface
136
138interface alloc
139 module procedure qcspaalloc
140end interface
141
143interface delete
144 module procedure qcspadelete
145end interface
146
147PRIVATE
148PUBLIC spa_nvar, spa_btable, spa_a, spa_b, qcspatype, init, alloc, delete, &
149 qcspatri, quaconspa
150
151
152contains
153
156
157subroutine qcspainit(qcspa,v7d,var, timei, timef, coordmin, coordmax, data_id_in,extremepath,spatialpath, &
158#ifdef HAVE_DBALLE
159 dsne,usere,passworde,&
160 dsnspa,userspa,passwordspa,&
161#endif
162 height2level,operation,categoryappend)
163
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
178
179#ifdef HAVE_DBALLE
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
190#endif
191
192integer :: iuni,i
193TYPE(vol7d_network) :: network
194character(len=512) :: filepathspa
195character(len=512) :: a_name
196
197call l4f_launcher(a_name,a_name_append=trim(subcategoryspa)//"."//trim(categoryappend))
198qcspa%category=l4f_category_get(a_name)
199
200call delete(qcspa%tri)
201
202nullify ( qcspa%data_id_in )
203nullify ( qcspa%data_id_out )
204
205
206! riporto il volume dati nel mio oggetto
207qcspa%v7d => v7d
208
209qcspa%operation=optio_c(operation,20)
210filepathspa=optio_c(spatialpath,512)
211
212!check options
213if (qcspa%operation /= "gradient" .and. qcspa%operation /= "run") then
214 call l4f_category_log(qcspa%category,l4f_error,"operation is wrong: "//qcspa%operation)
215 call raise_error()
216end if
217
218!inglobe id
219if (present(data_id_in))then
220 qcspa%data_id_in => data_id_in
221end if
222
223! load extreme
224call init(qcspa%qccli,v7d,var, timei, timef, data_id_in,&
225 macropath=cmiss, climapath=cmiss, extremepath=extremepath, &
226#ifdef HAVE_DBALLE
227 dsncli=cmiss,dsnextreme=dsne,user=usere,password=passworde,&
228#endif
229 height2level=height2level,categoryappend=categoryappend)
230
231
232! now load spatial clima
233
234if (qcspa%operation == "run") then
235
236 call init(network,"qcspa-ndi")
237
238#ifdef HAVE_DBALLE
239 call optio(dsnspa,ldsnspa)
240 call optio(userspa,luserspa)
241 call optio(passwordspa,lpasswordspa)
242
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")
245 call raise_error()
246 end if
247
248 if (.not. c_e(ldsnspa)) then
249
250#endif
251
252 if (.not. c_e(filepathspa)) then
253 filepathspa=get_package_filepath('qcspa-ndi.v7d', filetype_data)
254 end if
255
256 if (c_e(filepathspa))then
257
258 select case (trim(lowercase(suffixname(filepathspa))))
259
260 case("v7d")
261 iuni=getunit()
262 call import(qcspa%clima,filename=filepathspa,unit=iuni)
263 close (unit=iuni)
264
265#ifdef HAVE_DBALLE
266 case("bufr")
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)
272#endif
274 case default
275 call l4f_category_log(qcspa%category,l4f_error,&
276 "file type not supported (use .v7d or .bufr suffix only): "//trim(filepathspa))
277 call raise_error()
278 end select
279
280 else
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()
284 end if
285
286#ifdef HAVE_DBALLE
287 else
288
289 call l4f_category_log(qcspa%category,l4f_debug,"init v7d_dballespa")
290 call init(v7d_dballespa,dsn=ldsnspa,user=luserspa,password=lpasswordspa,write=.false.,&
291 file=.false.,categoryappend=trim(a_name)//".spa")
292 call l4f_category_log(qcspa%category,l4f_debug,"import v7d_dballespa")
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)
298 end if
299#endif
300end if
302
303return
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)
313type(qcspatype),intent(in out) :: qcspa
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
332
333
334integer :: status
335TYPE(geo_proj) :: geoproj
336REAL(kind=fp_geo) :: lat(size(qcspa%v7d%ana)),lon(size(qcspa%v7d%ana))
337CHARACTER(len=80) :: proj_type_l ! local type of projection
338double precision :: lov_l, latin1_l,latin2_l
339integer :: projection_center_flag_l
340
341proj_type_l = optio_c(proj_type,80)
342
343lov_l = optio_d(lov)
344latin1_l = optio_d(latin1)
345latin2_l = optio_d(latin2)
346projection_center_flag_l=optio_l(projection_center_flag)
347
348if (.not. c_e(proj_type_l)) then
349 proj_type_l = "lambert"
350 lov_l = 10.d0
351 latin1_l = 60.d0
352 latin2_l = 30.d0
353 projection_center_flag_l=1
354end if
355
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)
361
362call getval(qcspa%v7d%ana%coord, lon, lat)
363
364!print*,"size",size(lon),size(lat)
365!print*,lat,lon
366call proj(geoproj,lon,lat,qcspa%co%x,qcspa%co%y)
367!print*,"size x y ",size(qcspa%x),size(qcspa%y)
368!print*,qcspa%x,qcspa%y
369
370!triangulate
371status = triangles_compute(qcspa%co,qcspa%tri)
372
373!qcspa%nt,qcspa%ipt,qcspa%nl,qcspa%ipl)
374
375if (status /= 0) then
376 call l4f_category_log(qcspa%category,l4f_error,"contng error status="//t2c(status))
377 !call raise_error()
378end if
379
380end subroutine qcspatri
381
382
384subroutine qcspaalloc(qcspa)
385 ! pseudo costruttore con distruttore automatico
386
387type(qcspatype),intent(in out) :: qcspa
388
389integer :: istatt
390integer :: sh(5)
391
392! se ti sei dimenticato di deallocare ci penso io
393call qcspadealloc(qcspa)
394
395
396!!$if (associated (qcspa%v7d%dativar%r )) then
397!!$ nv=size(qcspa%v7d%dativar%r)
398!!$
399!!$ allocate(qcspa%valminr(nv),stat=istat)
400!!$ istatt=istatt+istat
401!!$ allocate(qcspa%valmaxr(nv),stat=istat)
402!!$ istatt=istatt+istat
403!!$
404!!$ if (istatt /= 0) ier=1
405!!$
406!!$end if
407
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)
411 if (istatt /= 0)then
412 call l4f_category_log(qcspa%category,l4f_error,"allocate error")
413 call raise_error("allocate error")
414 else
415 qcspa%data_id_out=imiss
416 end if
417end if
418
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))
423end if
424
425end subroutine qcspaalloc
426
427
429subroutine qcspadealloc(qcspa)
430 ! pseudo distruttore
431
432type(qcspatype),intent(in out) :: qcspa
433
434!!$if ( associated ( qcspa%valminr)) then
435!!$ deallocate(qcspa%valminr)
436!!$end if
437!!$
438!!$if ( associated ( qcspa%valmaxr)) then
439!!$ deallocate(qcspa%valmaxr)
440!!$end if
441
442if (associated(qcspa%data_id_out)) then
443 deallocate (qcspa%data_id_out)
444 nullify (qcspa%data_id_out)
445end if
446call delete(qcspa%tri)
447if (associated(qcspa%co)) deallocate(qcspa%co)
448
449end subroutine qcspadealloc
450
451
453
454
455subroutine qcspadelete(qcspa)
456 ! decostruttore a mezzo
457type(qcspatype),intent(in out) :: qcspa
458
459call qcspadealloc(qcspa)
460
461call delete(qcspa%qccli)
462
463qcspa%ndp=imiss
464
465!delete logger
466call l4f_category_delete(qcspa%category)
467
468return
469end subroutine qcspadelete
470
471
474
475SUBROUTINE quaconspa (qcspa,timetollerance,noborder,battrinv,battrcli,battrout,&
476 anamask,timemask,levelmask,timerangemask,varmask,networkmask)
477
478
479type(qcspatype),intent(in out) :: qcspa
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(:)
491
492 !REAL(kind=fp_geo) :: lat,lon
493 !local
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))
497
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
501 !integer, allocatable :: indcanav(:)
502
503 !TYPE(vol7d_ana) :: ana
504TYPE(datetime) :: time
505!YPE(vol7d_level):: level
506TYPE(vol7d_network):: network
507type(timedelta) :: deltato,deltat
508
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
513!!$CHARACTER(len=vol7d_ana_lenident) :: ident
514character(len=512) :: filename
515logical :: exist
516integer :: ind
517
518 !call qcspa_validate (qcspa)
519
520if (size(qcspa%v7d%ana) < 3 ) then
521 call l4f_category_log(qcspa%category,l4f_warn,"number of station < 3; do nothing")
522 return
523end if
524
525!localize optional parameter
526if (present(battrinv))then
527 indbattrinv = index_c(qcspa%v7d%datiattr%b(:)%btable, battrinv)
528else
529 indbattrinv = index_c(qcspa%v7d%datiattr%b(:)%btable, qcattrvarsbtables(1))
530end if
531
532if (present(battrcli))then
533 indbattrcli = index_c(qcspa%v7d%datiattr%b(:)%btable, battrcli)
534else
535 indbattrcli = index_c(qcspa%v7d%datiattr%b(:)%btable, qcattrvarsbtables(2))
536end if
537
538if (present(battrout))then
539 indbattrout = index_c(qcspa%v7d%datiattr%b(:)%btable, battrout)
540else
541 indbattrout = index_c(qcspa%v7d%datiattr%b(:)%btable, qcattrvarsbtables(4))
542end if
543
544
545! some checks on input
546!if (indbattrinv <=0 .or. indbattrcli <= 0 .or. indbattrout <= 0 ) then
547if (indbattrout <= 0 ) then
548
549 call l4f_category_log(qcspa%category,l4f_error,"error finding attribute index for output")
550 call raise_error()
551
552end if
553
554!!$if (qcspa%operation == "gradient") then
555!!$
556!!$ !check for gradient operation
557!!$ if ( size(qcspa%v7d%level) > 1 .or.&
558!!$ size(qcspa%v7d%timerange) > 1 .or.&
559!!$ size(qcspa%v7d%dativar%r) > 1 ) then
560!!$ call l4f_category_log(qcspa%category,L4F_ERROR,"gradient operation manage one level/timerange/var only")
561!!$ call raise_error()
562!!$ end if
563!!$
564!!$end if
565
566! set other local variable from optional parameter
567if(present(anamask)) then
568 anamaskl = anamask
569else
570 anamaskl = .true.
571endif
572if(present(timemask)) then
573 timemaskl = timemask
574else
575 timemaskl = .true.
576endif
577if(present(levelmask)) then
578 levelmaskl = levelmask
579else
580 levelmaskl = .true.
581endif
582if(present(timerangemask)) then
583 timerangemaskl = timerangemask
584else
585 timerangemaskl = .true.
586endif
587if(present(varmask)) then
588 varmaskl = varmask
589else
590 varmaskl = .true.
591endif
592if(present(networkmask)) then
593 networkmaskl = networkmask
594else
595 networkmaskl = .true.
596endif
597
598! do not touch data that do not pass QC
599qcspa%v7d%voldatiattrb(:,:,:,:,:,:,indbattrout)=ibmiss
600
601!print *,"prima normalize"
602!print *,qcspa%v7d%voldatir
603! normalize data in space and time
604call vol7d_normalize_data(qcspa%qccli)
605!print *,"dopo normalize"
606!print *,qcspa%v7d%voldatir
607
608! triangulate
609call qcspatri(qcspa)
610
611
612
613! compute some index for spatial clima
614!! compute the conventional generic datetime
615!!cyclicdt = cyclicdatetime_new(chardate="/////////") !TMMGGhhmm
616time=cyclicdatetime_to_conventional(cyclicdatetime_new(chardate="/////////")) !TMMGGhhmm
617!!call init(time, year=1007, month=1, day=1, hour=01, minute=01)
618
619
620if (qcspa%operation == "run") then
621 call init(network,"qcspa-ndi")
622 !!indcana = firsttrue(qcspa%clima%ana == ana)
623 indcnetwork = index(qcspa%clima%network , network)
624 indctime = index(qcspa%clima%time , time)
625end if
626
627do indtime=1,size(qcspa%v7d%time)
628 if (.not.timemaskl(indtime)) cycle
629 call l4f_category_log(qcspa%category,l4f_info,&
630 "Check time:"//t2c(qcspa%v7d%time(indtime)) )
631
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)
635
636 ind=index_c(spa_btable,qcspa%v7d%dativar%r(inddativarr)%btable)
637
638 if (qcspa%operation == "gradient") then
639 ! open file to write gradient
640
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)//&
644 ".grad"
645
646 call l4f_category_log(qcspa%category,l4f_info,"try to open gradient file; filename below")
647 call l4f_category_log(qcspa%category,l4f_info,filename)
648
649 inquire(file=filename, exist=exist)
650
651 grunit=getunit()
652 if (grunit /= -1) then
653 !open (unit=grunit, file=t2c(timei)//"_"//t2c(timef)//".grad",STATUS='UNKNOWN', form='FORMATTED')
654 open (grunit, file=filename ,status='UNKNOWN', form='FORMATTED',position='APPEND')
655 end if
656 ! say we have to write header in file
657 if (.not. exist) then
658 call l4f_category_log(qcspa%category,l4f_info,"write header in gradient file")
659 write (grunit,*) &
660 qcspa%v7d%level(indlevel), &
661 qcspa%v7d%timerange(indtimerange), &
662 qcspa%v7d%dativar%r(inddativarr)
663 end if
664 end if
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
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.
Index method.
Emit log message for a category with specific priority.
Test di dato invalidato.
Definition modqc.F90:411
Check data validity based on gross error check.
Definition modqc.F90:406
Allocazione di memoria.
Definition modqcspa.F90:326
Cancellazione.
Definition modqcspa.F90:331
Inizializzazione.
Definition modqcspa.F90:321
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.
Definition modqc.F90:357
Controllo di qualità climatico.
Definition modqccli.F90:295
Controllo di qualità spaziale.
Definition modqcspa.F90:273
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à
Definition modqccli.F90:322
Oggetto principale per il controllo di qualità
Definition modqcspa.F90:304
Definisce un oggetto contenente i volumi anagrafica e dati e tutti i descrittori delle loro dimension...

Generated with Doxygen.