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 
85 module modqcspa
86 
88 use geo_proj_class
90 use vol7d_class
91 use modqc
92 use modqccli
94 use log4fortran
97 !use array_utilities
98 !use io_units
99 #ifdef HAVE_DBALLE
101 #endif
102 
103 implicit none
104 
105 
106 character (len=255),parameter:: subcategoryspa="QCspa"
107 
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/)
114 
116 type :: 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
129 end type qcspatype
130 
131 
133 interface init
134  module procedure qcspainit
135 end interface
136 
138 interface alloc
139  module procedure qcspaalloc
140 end interface
141 
143 interface delete
144  module procedure qcspadelete
145 end interface
146 
147 PRIVATE
148 PUBLIC spa_nvar, spa_btable, spa_a, spa_b, qcspatype, init, alloc, delete, &
149  qcspatri, quaconspa
150 
151 
152 contains
153 
156 
157 subroutine 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 
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
178 
179 #ifdef HAVE_DBALLE
180 type (vol7d_dballe) :: v7d_dballespa
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
190 #endif
191 
192 integer :: iuni,i
193 TYPE(vol7d_network) :: network
194 character(len=512) :: filepathspa
195 character(len=512) :: a_name
196 
197 call l4f_launcher(a_name,a_name_append=trim(subcategoryspa)//"."//trim(categoryappend))
198 qcspa%category=l4f_category_get(a_name)
199 
200 call delete(qcspa%tri)
201 
202 nullify ( qcspa%data_id_in )
203 nullify ( qcspa%data_id_out )
204 
205 
206 ! riporto il volume dati nel mio oggetto
207 qcspa%v7d => v7d
208 
209 qcspa%operation=optio_c(operation,20)
210 filepathspa=optio_c(spatialpath,512)
211 
212 !check options
213 if (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()
216 end if
217 
218 !inglobe id
219 if (present(data_id_in))then
220  qcspa%data_id_in => data_id_in
221 end if
222 
223 ! load extreme
224 call 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 
234 if (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
273 
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)
297 
298  end if
299 #endif
300 end if
301 
302 
303 return
304 end subroutine qcspainit
305 
306 
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)
312 
313 type(qcspatype),intent(in out) :: qcspa !< Oggetto per il controllo climatico
314 CHARACTER(len=*),INTENT(in),OPTIONAL :: proj_type !< type of projection
315 DOUBLE PRECISION,INTENT(in),OPTIONAL :: lov !< line of view, also known as reference longitude or orientation of the grid (polar projections)
316 INTEGER,INTENT(in),OPTIONAL :: zone !< Earth zone (mainly for UTM), sets lov to the correct zone central meridian
317 DOUBLE PRECISION,INTENT(in),OPTIONAL :: xoff !< offset on x axis (false easting)
318 DOUBLE PRECISION,INTENT(in),OPTIONAL :: yoff !< offset on y axis (false northing)
319 DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_south_pole !< longitude of the southern pole of projection
320 DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_south_pole !< latitude of the southern pole of projection
321 DOUBLE PRECISION,INTENT(in),OPTIONAL :: angle_rotation !< angle of rotation of projection
322 DOUBLE PRECISION,INTENT(in),OPTIONAL :: longitude_stretch_pole !< longitude of the pole of stretching
323 DOUBLE PRECISION,INTENT(in),OPTIONAL :: latitude_stretch_pole !< latitude of the pole of stretching
324 DOUBLE PRECISION,INTENT(in),OPTIONAL :: stretch_factor !< stretching factor
325 DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin1 !< first standard latitude from main pole (Lambert)
326 DOUBLE PRECISION,INTENT(in),OPTIONAL :: latin2 !< second standard latitude from main pole (Lambert)
327 DOUBLE PRECISION,INTENT(in),OPTIONAL :: lad !< latitude at which dx and dy (in m) are specified (Lambert, grib2 only)
328 INTEGER,INTENT(in),OPTIONAL :: projection_center_flag !< flag indicating which pole is represented
329 DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_smaj_axis !< Earth semi-major axis
330 DOUBLE PRECISION,INTENT(in),OPTIONAL :: ellips_flatt !< Earth flattening
331 INTEGER,INTENT(in),OPTIONAL :: ellips_type !< number in the interval [1,nellips] indicating a predefined ellipsoid, alternative to the previous arguments
332 
333 
334 integer :: status
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 ! local type of projection
338 double precision :: lov_l, latin1_l,latin2_l
339 integer :: projection_center_flag_l
340 
341 proj_type_l = optio_c(proj_type,80)
342 
343 lov_l = optio_d(lov)
344 latin1_l = optio_d(latin1)
345 latin2_l = optio_d(latin2)
346 projection_center_flag_l=optio_l(projection_center_flag)
347 
348 .not.if ( 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
354 end if
355 
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)
361 
362 call getval(qcspa%v7d%ana%coord, lon, lat)
363 
364 !print*,"size",size(lon),size(lat)
365 !print*,lat,lon
366 call 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
371 status = triangles_compute(qcspa%co,qcspa%tri)
372 
373 !qcspa%nt,qcspa%ipt,qcspa%nl,qcspa%ipl)
374 
375 if (status /= 0) then
376  call l4f_category_log(qcspa%category,L4F_ERROR,"contng error status="//t2c(status))
377  !call raise_error()
378 end if
379 
380 end subroutine qcspatri
381 
382 
383 !>\brief Allocazioni di memoria
384 subroutine qcspaalloc(qcspa)
385  ! pseudo costruttore con distruttore automatico
386 
387 type(qcspatype),intent(in out) :: qcspa !< Oggetto per il controllo climatico
388 
389 integer :: istatt
390 integer :: sh(5)
391 
392 ! se ti sei dimenticato di deallocare ci penso io
393 call 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 
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)
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
417 end if
418 
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))
423 end if
424 
425 end subroutine qcspaalloc
426 
427 
428 !>\brief Deallocazione della memoria
429 subroutine qcspadealloc(qcspa)
430  ! pseudo distruttore
431 
432 type(qcspatype),intent(in out) :: qcspa !< Oggetto per l controllo climatico
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 
442 if (associated(qcspa%data_id_out)) then
443  deallocate (qcspa%data_id_out)
444  nullify (qcspa%data_id_out)
445 end if
446 call delete(qcspa%tri)
447 if (associated(qcspa%co)) deallocate(qcspa%co)
448 
449 end subroutine qcspadealloc
450 
451 
452 !>\brief Cancellazione
453 
454 
455 subroutine qcspadelete(qcspa)
456  ! decostruttore a mezzo
457 type(qcspatype),intent(in out) :: qcspa !< Oggetto per il controllo climatico
458 
459 call qcspadealloc(qcspa)
460 
461 call delete(qcspa%qccli)
462 
463 qcspa%ndp=imiss
464 
465 !delete logger
466 call l4f_category_delete(qcspa%category)
467 
468 return
469 end subroutine qcspadelete
470 
471 
472 à!>\brief Controllo di Qualit spaziale.
473 èà!!Questo il vero e proprio controllo di qualit spaziale.
474 
475 SUBROUTINE quaconspa (qcspa,timetollerance,noborder,battrinv,battrcli,battrout,&
476  anamask,timemask,levelmask,timerangemask,varmask,networkmask)
477 
478 
479 àtype(qcspatype),intent(in out) :: qcspa !< Oggetto per il controllo di qualit
480 type(timedelta),intent(in) :: timetollerance !< time tollerance to compare nearest stations
481 logical,intent(in),optional :: noborder !< Exclude border from QC
482 character (len=10) ,intent(in),optional :: battrinv !< attributo invalidated in input
483 character (len=10) ,intent(in),optional :: battrcli !< attributo con la confidenza climatologica in input
484 character (len=10) ,intent(in),optional :: battrout !< attributo con la confidenza spaziale in output
485 logical ,intent(in),optional :: anamask(:) !< Filtro sulle anagrafiche
486 logical ,intent(in),optional :: timemask(:) !< Filtro sul tempo
487 logical ,intent(in),optional :: levelmask(:) !< Filtro sui livelli
488 logical ,intent(in),optional :: timerangemask(:) !< filtro sui timerange
489 logical ,intent(in),optional :: varmask(:) !< Filtro sulle variabili
490 logical ,intent(in),optional :: networkmask(:) !< Filtro sui network
491 
492  !REAL(kind=fp_geo) :: lat,lon
493  !local
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))
497 
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
501  !integer, allocatable :: indcanav(:)
502 
503  !TYPE(vol7d_ana) :: ana
504 TYPE(datetime) :: time
505 !YPE(vol7d_level):: level
506 TYPE(vol7d_network):: network
507 type(timedelta) :: deltato,deltat
508 
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
513 !!$CHARACTER(len=vol7d_ana_lenident) :: ident
514 character(len=512) :: filename
515 logical :: exist
516 integer :: ind
517 
518  !call qcspa_validate (qcspa)
519 
520 if (size(qcspa%v7d%ana) < 3 ) then
521  call l4f_category_log(qcspa%category,L4F_WARN,"number of station < 3; do nothing")
522  return
523 end if
524 
525 !localize optional parameter
526 if (present(battrinv))then
527  indbattrinv = index_c(qcspa%v7d%datiattr%b(:)%btable, battrinv)
528 else
529  indbattrinv = index_c(qcspa%v7d%datiattr%b(:)%btable, qcattrvarsbtables(1))
530 end if
531 
532 if (present(battrcli))then
533  indbattrcli = index_c(qcspa%v7d%datiattr%b(:)%btable, battrcli)
534 else
535  indbattrcli = index_c(qcspa%v7d%datiattr%b(:)%btable, qcattrvarsbtables(2))
536 end if
537 
538 if (present(battrout))then
539  indbattrout = index_c(qcspa%v7d%datiattr%b(:)%btable, battrout)
540 else
541  indbattrout = index_c(qcspa%v7d%datiattr%b(:)%btable, qcattrvarsbtables(4))
542 end if
543 
544 
545 ! some checks on input
546 .or..or.!if (indbattrinv <=0 indbattrcli <= 0 indbattrout <= 0 ) then
547 if (indbattrout <= 0 ) then
548 
549  call l4f_category_log(qcspa%category,L4F_ERROR,"error finding attribute index for output")
550  call raise_error()
551 
552 end if
553 
554 !!$if (qcspa%operation == "gradient") then
555 !!$
556 !!$ !check for gradient operation
557 .or.!!$ if ( size(qcspa%v7d%level) > 1 &
558 .or.!!$ size(qcspa%v7d%timerange) > 1 &
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
567 if(present(anamask)) then
568  anamaskl = anamask
569 else
570  anamaskl = .true.
571 endif
572 if(present(timemask)) then
573  timemaskl = timemask
574 else
575  timemaskl = .true.
576 endif
577 if(present(levelmask)) then
578  levelmaskl = levelmask
579 else
580  levelmaskl = .true.
581 endif
582 if(present(timerangemask)) then
583  timerangemaskl = timerangemask
584 else
585  timerangemaskl = .true.
586 endif
587 if(present(varmask)) then
588  varmaskl = varmask
589 else
590  varmaskl = .true.
591 endif
592 if(present(networkmask)) then
593  networkmaskl = networkmask
594 else
595  networkmaskl = .true.
596 endif
597 
598 ! do not touch data that do not pass QC
599 qcspa%v7d%voldatiattrb(:,:,:,:,:,:,indbattrout)=ibmiss
600 
601 !print *,"prima normalize"
602 !print *,qcspa%v7d%voldatir
603 ! normalize data in space and time
604 call vol7d_normalize_data(qcspa%qccli)
605 !print *,"dopo normalize"
606 !print *,qcspa%v7d%voldatir
607 
608 ! triangulate
609 call qcspatri(qcspa)
610 
611 
612 
613 ! compute some index for spatial clima
614 !! compute the conventional generic datetime
615 !!cyclicdt = cyclicdatetime_new(chardate="/////////") !TMMGGhhmm
616 time=cyclicdatetime_to_conventional(cyclicdatetime_new(chardate="/////////")) !TMMGGhhmm
617 !!call init(time, year=1007, month=1, day=1, hour=01, minute=01)
618 
619 
620 if (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)
625 end if
626 
627 do indtime=1,size(qcspa%v7d%time)
628 .not. if (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 .not. if ( 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 .not..or..not..or. if ( anamaskl(indana) levelmaskl(indlevel) &
674 .not..or..not..or..not. timerangemaskl(indtimerange) varmaskl(inddativarr) networkmaskl(indnetwork)) cycle
675 
676  datoqui = qcspa%v7d%voldatir (indana ,indtime ,indlevel ,indtimerange ,inddativarr, indnetwork )
677 .not. if ( 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 .not. if( 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 .or..or..or. if ( indctime <= 0 indclevel <= 0 indctimerange <= 0 indcdativarr <= 0 &
726 .or. indcnetwork <= 0 ) cycle
727  end if
728 
729 .and. if (optio_log(noborder) 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 .EQ. IF(qcspa%tri%IPT(3*IT-2)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 .EQ. IF(qcspa%tri%IPT(3*IT-1)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 .EQ. IF(qcspa%tri%IPT(3*IT)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 .LT.!!$ IF(IVERT(I)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 .NE. IF(IVERT(IV)IVERT(KK))THEN
784  IV=IV+1
785  IVERT(IV)=IVERT(KK)
786  ENDIF
787  END DO
788 .GT. IF (IVITROV)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 .not. if( 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 .not. if ( 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 .not. if ( 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 .or..not..and. if ((deltat < deltato c_e(deltato)) deltat <= timetollerance ) then
846  datola = datila(iindtime)
847  deltato = deltat
848  end if
849  end do
850  end do
851 
852 
853 .NOT. IF(C_E(datola)) cycle
854  ! distanza tra le due stazioni
855  dist = DISTANZA (qcspa%co(INDANA),qcspa%co(IVERT(I)))
856 .EQ. IF (DIST0.)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 .or. IF (ipos == ivb 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 .and. if ( c_e(climaquii) c_e(climaquif )) then
923 
924 .and..or. if ( (gradmin >= climaquii gradmin < climaquif) &
925 .and..or. (indcana == 1 gradmin < climaquii) &
926 .and. (indcana == size(qcspa%clima%ana) 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 
988 !> \example v7d_qcspa.F90
989 !! Sample program for module qcspa
990 
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.
Index method.
Emit log message for a category with specific priority.
Allocazione di memoria.
Definition: modqcspa.F90:336
Cancellazione.
Definition: modqcspa.F90:341
Inizializzazione.
Definition: modqcspa.F90:331
Check for problems return 0 if all check passed print diagnostics with log4f.
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:367
Controllo di qualità climatico.
Definition: modqccli.F90:305
Controllo di qualità spaziale.
Definition: modqcspa.F90:283
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
Oggetto principale per il controllo di qualità
Definition: modqcspa.F90:314

Generated with Doxygen.