libsim  Versione 7.1.9
volgrid6d_class.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 #include "config.h"
19 
31 
52 MODULE volgrid6d_class
55 USE vol7d_class
57 USE geo_proj_class
58 USE grid_class
62 USE log4fortran
64 USE grid_id_class
67 !USE file_utilities
68 #ifdef HAVE_DBALLE
70 #endif
71 IMPLICIT NONE
72 
73 character (len=255),parameter:: subcategory="volgrid6d_class"
74 
76 type volgrid6d
77  type(griddim_def) :: griddim
78  TYPE(datetime),pointer :: time(:)
79  TYPE(vol7d_timerange),pointer :: timerange(:)
80  TYPE(vol7d_level),pointer :: level(:)
81  TYPE(volgrid6d_var),pointer :: var(:)
82  TYPE(grid_id),POINTER :: gaid(:,:,:,:)
83  REAL,POINTER :: voldati(:,:,:,:,:,:)
84  integer :: time_definition
85  integer :: category = 0
86 end type volgrid6d
87 
89 INTERFACE init
90  MODULE PROCEDURE volgrid6d_init
91 END INTERFACE
92 
95 INTERFACE delete
96  MODULE PROCEDURE volgrid6d_delete, volgrid6dv_delete
97 END INTERFACE
98 
101 INTERFACE import
102  MODULE PROCEDURE volgrid6d_read_from_file
103  MODULE PROCEDURE import_from_gridinfo, import_from_gridinfovv, &
104  volgrid6d_import_from_file
105 END INTERFACE
106 
109 INTERFACE export
110  MODULE PROCEDURE volgrid6d_write_on_file
111  MODULE PROCEDURE export_to_gridinfo, export_to_gridinfov, export_to_gridinfovv,&
112  volgrid6d_export_to_file
113 END INTERFACE
114 
115 ! methods for computing transformations through an initialised
116 ! grid_transform object, probably too low level to be interfaced
117 INTERFACE compute
118  MODULE PROCEDURE volgrid6d_transform_compute, volgrid6d_v7d_transform_compute,&
119  v7d_volgrid6d_transform_compute, v7d_v7d_transform_compute
120 END INTERFACE
121 
124 INTERFACE transform
125  MODULE PROCEDURE volgrid6d_transform, volgrid6dv_transform,&
126  volgrid6d_v7d_transform, volgrid6dv_v7d_transform, v7d_volgrid6d_transform, &
127  v7d_v7d_transform
128 END INTERFACE
129 
130 INTERFACE wind_rot
131  MODULE PROCEDURE vg6d_wind_rot
132 END INTERFACE
133 
134 INTERFACE wind_unrot
135  MODULE PROCEDURE vg6d_wind_unrot
136 END INTERFACE
137 
140 INTERFACE display
141  MODULE PROCEDURE display_volgrid6d,display_volgrid6dv
142 END INTERFACE
143 
155 INTERFACE rounding
156  MODULE PROCEDURE vg6d_rounding, vg6dv_rounding
157 END INTERFACE
158 
159 private
160 
161 PUBLIC volgrid6d,init,delete,export,import,compute,transform, &
162  wind_rot,wind_unrot,vg6d_c2a,display,volgrid6d_alloc,volgrid6d_alloc_vol, &
163  volgrid_get_vol_2d, volgrid_set_vol_2d, volgrid_get_vol_3d, volgrid_set_vol_3d
164 PUBLIC rounding, vg6d_reduce
165 
166 CONTAINS
167 
168 
173 SUBROUTINE volgrid6d_init(this, griddim, time_definition, categoryappend)
174 TYPE(volgrid6d) :: this
175 TYPE(griddim_def),OPTIONAL :: griddim
176 INTEGER,INTENT(IN),OPTIONAL :: time_definition
177 CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
178 
179 character(len=512) :: a_name
180 
181 if (present(categoryappend))then
182  call l4f_launcher(a_name,a_name_append=trim(subcategory)//"."//trim(categoryappend))
183 else
184  call l4f_launcher(a_name,a_name_append=trim(subcategory))
185 endif
186 this%category=l4f_category_get(a_name)
187 
188 #ifdef DEBUG
189 call l4f_category_log(this%category,l4f_debug,"init")
190 #endif
191 
192 call init(this%griddim)
193 
194 if (present(griddim))then
195  call copy (griddim,this%griddim)
196 end if
197 
198 CALL vol7d_var_features_init() ! initialise var features table once
199 
200 if(present(time_definition)) then
201  this%time_definition = time_definition
202 else
203  this%time_definition = 0 !default to reference time
204 end if
205 
206 nullify (this%time,this%timerange,this%level,this%var)
207 nullify (this%gaid,this%voldati)
208 
209 END SUBROUTINE volgrid6d_init
210 
211 
222 SUBROUTINE volgrid6d_alloc(this, dim, ntime, nlevel, ntimerange, nvar, ini)
223 TYPE(volgrid6d),INTENT(inout) :: this
224 TYPE(grid_dim),INTENT(in),OPTIONAL :: dim
225 INTEGER,INTENT(in),OPTIONAL :: ntime
226 INTEGER,INTENT(in),OPTIONAL :: nlevel
227 INTEGER,INTENT(in),OPTIONAL :: ntimerange
228 INTEGER,INTENT(in),OPTIONAL :: nvar
229 LOGICAL,INTENT(in),OPTIONAL :: ini
230 
231 INTEGER :: i, stallo
232 LOGICAL :: linit
233 
234 #ifdef DEBUG
235 call l4f_category_log(this%category,l4f_debug,"alloc")
236 #endif
237 
238 IF (PRESENT(ini)) THEN
239  linit = ini
240 ELSE
241  linit = .false.
242 ENDIF
243 
244 
245 if (present(dim)) call copy (dim,this%griddim%dim)
246 
247 
248 IF (PRESENT(ntime)) THEN
249  IF (ntime >= 0) THEN
250  IF (ASSOCIATED(this%time)) DEALLOCATE(this%time)
251 #ifdef DEBUG
252  call l4f_category_log(this%category,l4f_debug,"alloc ntime "//to_char(ntime))
253 #endif
254  ALLOCATE(this%time(ntime),stat=stallo)
255  if (stallo /=0)then
256  call l4f_category_log(this%category,l4f_fatal,"allocating memory")
257  CALL raise_fatal_error()
258  end if
259  IF (linit) THEN
260  DO i = 1, ntime
261  this%time(i) = datetime_miss
262 ! CALL init(this%time(i)) ! senza argomento inizializza a zero non missing
263  ! baco di datetime?
264  ENDDO
265  ENDIF
266  ENDIF
267 ENDIF
268 IF (PRESENT(nlevel)) THEN
269  IF (nlevel >= 0) THEN
270  IF (ASSOCIATED(this%level)) DEALLOCATE(this%level)
271 #ifdef DEBUG
272  call l4f_category_log(this%category,l4f_debug,"alloc nlevel "//to_char(nlevel))
273 #endif
274  ALLOCATE(this%level(nlevel),stat=stallo)
275  if (stallo /=0)then
276  call l4f_category_log(this%category,l4f_fatal,"allocating memory")
277  CALL raise_fatal_error()
278  end if
279  IF (linit) THEN
280  DO i = 1, nlevel
281  CALL init(this%level(i))
282  ENDDO
283  ENDIF
284  ENDIF
285 ENDIF
286 IF (PRESENT(ntimerange)) THEN
287  IF (ntimerange >= 0) THEN
288  IF (ASSOCIATED(this%timerange)) DEALLOCATE(this%timerange)
289 #ifdef DEBUG
290  call l4f_category_log(this%category,l4f_debug,"alloc ntimerange "//to_char(ntimerange))
291 #endif
292  ALLOCATE(this%timerange(ntimerange),stat=stallo)
293  if (stallo /=0)then
294  call l4f_category_log(this%category,l4f_fatal,"allocating memory")
295  CALL raise_fatal_error()
296  end if
297  IF (linit) THEN
298  DO i = 1, ntimerange
299  CALL init(this%timerange(i))
300  ENDDO
301  ENDIF
302  ENDIF
303 ENDIF
304 IF (PRESENT(nvar)) THEN
305  IF (nvar >= 0) THEN
306  IF (ASSOCIATED(this%var)) DEALLOCATE(this%var)
307 #ifdef DEBUG
308  call l4f_category_log(this%category,l4f_debug,"alloc nvar "//to_char(nvar))
309 #endif
310  ALLOCATE(this%var(nvar),stat=stallo)
311  if (stallo /=0)then
312  call l4f_category_log(this%category,l4f_fatal,"allocating memory")
313  CALL raise_fatal_error()
314  end if
315  IF (linit) THEN
316  DO i = 1, nvar
317  CALL init(this%var(i))
318  ENDDO
319  ENDIF
320  ENDIF
321 ENDIF
322 
323 end SUBROUTINE volgrid6d_alloc
324 
325 
334 SUBROUTINE volgrid6d_alloc_vol(this, ini, inivol, decode)
335 TYPE(volgrid6d),INTENT(inout) :: this
336 LOGICAL,INTENT(in),OPTIONAL :: ini
337 LOGICAL,INTENT(in),OPTIONAL :: inivol
338 LOGICAL,INTENT(in),OPTIONAL :: decode
339 
340 INTEGER :: stallo
341 LOGICAL :: linivol
342 
343 #ifdef DEBUG
344 call l4f_category_log(this%category,l4f_debug,"start alloc_vol")
345 #endif
346 
347 IF (PRESENT(inivol)) THEN ! opposite condition, cannot use optio_log
348  linivol = inivol
349 ELSE
350  linivol = .true.
351 ENDIF
352 
353 IF (this%griddim%dim%nx > 0 .AND. this%griddim%dim%ny > 0) THEN
354 ! allocate dimension descriptors to a minimal size for having a
355 ! non-null volume
356  IF (.NOT. ASSOCIATED(this%var)) CALL volgrid6d_alloc(this, nvar=1, ini=ini)
357  IF (.NOT. ASSOCIATED(this%time)) CALL volgrid6d_alloc(this, ntime=1, ini=ini)
358  IF (.NOT. ASSOCIATED(this%level)) CALL volgrid6d_alloc(this, nlevel=1, ini=ini)
359  IF (.NOT. ASSOCIATED(this%timerange)) CALL volgrid6d_alloc(this, ntimerange=1, ini=ini)
360 
361  IF (optio_log(decode)) THEN
362  IF (.NOT.ASSOCIATED(this%voldati)) THEN
363 #ifdef DEBUG
364  CALL l4f_category_log(this%category,l4f_debug,"alloc voldati volume")
365 #endif
366 
367  ALLOCATE(this%voldati(this%griddim%dim%nx,this%griddim%dim%ny,&
368  SIZE(this%level), SIZE(this%time), &
369  SIZE(this%timerange), SIZE(this%var)),stat=stallo)
370  IF (stallo /= 0)THEN
371  CALL l4f_category_log(this%category,l4f_fatal,"allocating memory")
372  CALL raise_fatal_error()
373  ENDIF
374 
375 ! this is not really needed if we can check other flags for a full
376 ! field missing values
377  IF (linivol) this%voldati = rmiss
378  this%voldati = rmiss
379  ENDIF
380  ENDIF
381 
382  IF (.NOT.ASSOCIATED(this%gaid)) THEN
383 #ifdef DEBUG
384  CALL l4f_category_log(this%category,l4f_debug,"alloc gaid volume")
385 #endif
386  ALLOCATE(this%gaid(SIZE(this%level), SIZE(this%time), &
387  SIZE(this%timerange), SIZE(this%var)),stat=stallo)
388  IF (stallo /= 0)THEN
389  CALL l4f_category_log(this%category,l4f_fatal,"allocating memory")
390  CALL raise_fatal_error()
391  ENDIF
392 
393  IF (linivol) THEN
394 !!$ DO i=1 ,SIZE(this%gaid,1)
395 !!$ DO ii=1 ,SIZE(this%gaid,2)
396 !!$ DO iii=1 ,SIZE(this%gaid,3)
397 !!$ DO iiii=1 ,SIZE(this%gaid,4)
398 !!$ this%gaid(i,ii,iii,iiii) = grid_id_new() ! optimize?
399 !!$ ENDDO
400 !!$ ENDDO
401 !!$ ENDDO
402 !!$ ENDDO
403 
404  this%gaid = grid_id_new()
405  ENDIF
406  ENDIF
407 
408 ELSE
409  CALL l4f_category_log(this%category,l4f_fatal,'volgrid6d_alloc_vol: &
410  &trying to allocate a volume with an invalid or unspecified horizontal grid')
411  CALL raise_fatal_error()
412 ENDIF
413 
414 END SUBROUTINE volgrid6d_alloc_vol
415 
416 
430 SUBROUTINE volgrid_get_vol_2d(this, ilevel, itime, itimerange, ivar, voldati)
431 TYPE(volgrid6d),INTENT(in) :: this
432 INTEGER,INTENT(in) :: ilevel
433 INTEGER,INTENT(in) :: itime
434 INTEGER,INTENT(in) :: itimerange
435 INTEGER,INTENT(in) :: ivar
436 REAL,POINTER :: voldati(:,:)
437 
438 IF (ASSOCIATED(this%voldati)) THEN
439  voldati => this%voldati(:,:,ilevel,itime,itimerange,ivar)
440  RETURN
441 ELSE
442  IF (.NOT.ASSOCIATED(voldati)) THEN
443  ALLOCATE(voldati(this%griddim%dim%nx,this%griddim%dim%ny))
444  ENDIF
445  CALL grid_id_decode_data(this%gaid(ilevel,itime,itimerange,ivar), voldati)
446 ENDIF
447 
448 END SUBROUTINE volgrid_get_vol_2d
449 
450 
464 SUBROUTINE volgrid_get_vol_3d(this, itime, itimerange, ivar, voldati)
465 TYPE(volgrid6d),INTENT(in) :: this
466 INTEGER,INTENT(in) :: itime
467 INTEGER,INTENT(in) :: itimerange
468 INTEGER,INTENT(in) :: ivar
469 REAL,POINTER :: voldati(:,:,:)
470 
471 INTEGER :: ilevel
472 
473 IF (ASSOCIATED(this%voldati)) THEN
474  voldati => this%voldati(:,:,:,itime,itimerange,ivar)
475  RETURN
476 ELSE
477  IF (.NOT.ASSOCIATED(voldati)) THEN
478  ALLOCATE(voldati(this%griddim%dim%nx,this%griddim%dim%ny,SIZE(this%level)))
479  ENDIF
480  DO ilevel = 1, SIZE(this%level)
481  CALL grid_id_decode_data(this%gaid(ilevel,itime,itimerange,ivar), &
482  voldati(:,:,ilevel))
483  ENDDO
484 ENDIF
485 
486 END SUBROUTINE volgrid_get_vol_3d
487 
488 
500 SUBROUTINE volgrid_set_vol_2d(this, ilevel, itime, itimerange, ivar, voldati)
501 TYPE(volgrid6d),INTENT(inout) :: this
502 INTEGER,INTENT(in) :: ilevel
503 INTEGER,INTENT(in) :: itime
504 INTEGER,INTENT(in) :: itimerange
505 INTEGER,INTENT(in) :: ivar
506 REAL,INTENT(in) :: voldati(:,:)
507 
508 IF (ASSOCIATED(this%voldati)) THEN
509  RETURN
510 ELSE
511  CALL grid_id_encode_data(this%gaid(ilevel,itime,itimerange,ivar), voldati)
512 ENDIF
513 
514 END SUBROUTINE volgrid_set_vol_2d
515 
516 
528 SUBROUTINE volgrid_set_vol_3d(this, itime, itimerange, ivar, voldati)
529 TYPE(volgrid6d),INTENT(inout) :: this
530 INTEGER,INTENT(in) :: itime
531 INTEGER,INTENT(in) :: itimerange
532 INTEGER,INTENT(in) :: ivar
533 REAL,INTENT(in) :: voldati(:,:,:)
534 
535 INTEGER :: ilevel
536 
537 IF (ASSOCIATED(this%voldati)) THEN
538  RETURN
539 ELSE
540  DO ilevel = 1, SIZE(this%level)
541  CALL grid_id_encode_data(this%gaid(ilevel,itime,itimerange,ivar), &
542  voldati(:,:,ilevel))
543  ENDDO
544 ENDIF
545 
546 END SUBROUTINE volgrid_set_vol_3d
547 
548 
552 SUBROUTINE volgrid6d_delete(this)
553 TYPE(volgrid6d),INTENT(inout) :: this
554 
555 INTEGER :: i, ii, iii, iiii
556 
557 #ifdef DEBUG
558 call l4f_category_log(this%category,l4f_debug,"delete")
559 #endif
560 
561 if (associated(this%gaid))then
562 
563  DO i=1 ,SIZE(this%gaid,1)
564  DO ii=1 ,SIZE(this%gaid,2)
565  DO iii=1 ,SIZE(this%gaid,3)
566  DO iiii=1 ,SIZE(this%gaid,4)
567  CALL delete(this%gaid(i,ii,iii,iiii))
568  ENDDO
569  ENDDO
570  ENDDO
571  ENDDO
572  DEALLOCATE(this%gaid)
573 
574 end if
575 
576 call delete(this%griddim)
577 
578 ! call delete(this%time)
579 ! call delete(this%timerange)
580 ! call delete(this%level)
581 ! call delete(this%var)
582 
583 if (associated( this%time )) deallocate(this%time)
584 if (associated( this%timerange )) deallocate(this%timerange)
585 if (associated( this%level )) deallocate(this%level)
586 if (associated( this%var )) deallocate(this%var)
587 
588 if (associated(this%voldati))deallocate(this%voldati)
589 
590 
591  !chiudo il logger
592 call l4f_category_delete(this%category)
593 
594 END SUBROUTINE volgrid6d_delete
595 
596 
606 subroutine volgrid6d_write_on_file (this,unit,description,filename,filename_auto)
607 
608 TYPE(volgrid6d),INTENT(IN) :: this
609 integer,optional,intent(inout) :: unit
610 character(len=*),intent(in),optional :: filename
611 character(len=*),intent(out),optional :: filename_auto
612 character(len=*),INTENT(IN),optional :: description
613 
614 integer :: lunit
615 character(len=254) :: ldescription,arg,lfilename
616 integer :: ntime, ntimerange, nlevel, nvar
617 integer :: tarray(8)
618 logical :: opened,exist
619 
620 #ifdef DEBUG
621 call l4f_category_log(this%category,l4f_debug,"write on file")
622 #endif
623 
624 ntime=0
625 ntimerange=0
626 nlevel=0
627 nvar=0
628 
629 !call idate(im,id,iy)
630 call date_and_time(values=tarray)
631 call getarg(0,arg)
632 
633 if (present(description))then
634  ldescription=description
635 else
636  ldescription="Volgrid6d generated by: "//trim(arg)
637 end if
638 
639 if (.not. present(unit))then
640  lunit=getunit()
641 else
642  if (unit==0)then
643  lunit=getunit()
644  unit=lunit
645  else
646  lunit=unit
647  end if
648 end if
649 
650 lfilename=trim(arg)//".vg6d"
651 if (index(arg,'/',back=.true.) > 0) lfilename=lfilename(index(arg,'/',back=.true.)+1 : )
652 
653 if (present(filename))then
654  if (filename /= "")then
655  lfilename=filename
656  end if
657 end if
658 
659 if (present(filename_auto))filename_auto=lfilename
660 
661 
662 inquire(unit=lunit,opened=opened)
663 if (.not. opened) then
664  inquire(file=lfilename,exist=exist)
665  if (exist) CALL raise_error('file exist; cannot open new file')
666  if (.not.exist) open (unit=lunit,file=lfilename,form="UNFORMATTED")
667  !print *, "opened: ",lfilename
668 end if
669 
670 if (associated(this%time)) ntime=size(this%time)
671 if (associated(this%timerange)) ntimerange=size(this%timerange)
672 if (associated(this%level)) nlevel=size(this%level)
673 if (associated(this%var)) nvar=size(this%var)
674 
675 
676 write(unit=lunit)ldescription
677 write(unit=lunit)tarray
678 
679 call write_unit( this%griddim,lunit)
680 write(unit=lunit) ntime, ntimerange, nlevel, nvar
681 
682 !! prime 4 dimensioni
683 if (associated(this%time)) call write_unit(this%time, lunit)
684 if (associated(this%level)) write(unit=lunit)this%level
685 if (associated(this%timerange)) write(unit=lunit)this%timerange
686 if (associated(this%var)) write(unit=lunit)this%var
687 
688 
689 !! Volumi di valori dati
690 
691 if (associated(this%voldati)) write(unit=lunit)this%voldati
692 
693 if (.not. present(unit)) close(unit=lunit)
694 
695 end subroutine volgrid6d_write_on_file
696 
697 
704 subroutine volgrid6d_read_from_file (this,unit,filename,description,tarray,filename_auto)
705 
706 TYPE(volgrid6d),INTENT(OUT) :: this
707 integer,intent(inout),optional :: unit
708 character(len=*),INTENT(in),optional :: filename
709 character(len=*),intent(out),optional :: filename_auto
710 character(len=*),INTENT(out),optional :: description
711 integer,intent(out),optional :: tarray(8)
712 
713 integer :: ntime, ntimerange, nlevel, nvar
714 
715 character(len=254) :: ldescription,lfilename,arg
716 integer :: ltarray(8),lunit
717 logical :: opened,exist
718 
719 #ifdef DEBUG
720 call l4f_category_log(this%category,l4f_debug,"read from file")
721 #endif
722 
723 call getarg(0,arg)
724 
725 if (.not. present(unit))then
726  lunit=getunit()
727 else
728  if (unit==0)then
729  lunit=getunit()
730  unit=lunit
731  else
732  lunit=unit
733  end if
734 end if
735 
736 lfilename=trim(arg)//".vg6d"
737 if (index(arg,'/',back=.true.) > 0) lfilename=lfilename(index(arg,'/',back=.true.)+1 : )
738 
739 if (present(filename))then
740  if (filename /= "")then
741  lfilename=filename
742  end if
743 end if
744 
745 if (present(filename_auto))filename_auto=lfilename
746 
747 
748 inquire(unit=lunit,opened=opened)
749 if (.not. opened) then
750  inquire(file=lfilename,exist=exist)
751  IF (.NOT. exist) CALL raise_fatal_error('file '//trim(lfilename)//' does not exist, cannot open')
752  open (unit=lunit,file=lfilename,form="UNFORMATTED")
753 end if
754 
755 read(unit=lunit)ldescription
756 read(unit=lunit)ltarray
757 
758 call l4f_log(l4f_info,"Info: reading volgrid6d from file: "//trim(lfilename))
759 call l4f_log(l4f_info,"Info: description: "//trim(ldescription))
760 !call l4f_log("Info: written on ",ltarray)
761 
762 if (present(description))description=ldescription
763 if (present(tarray))tarray=ltarray
764 
765 
766 call read_unit( this%griddim,lunit)
767 read(unit=lunit) ntime, ntimerange, nlevel, nvar
768 
769 
770 call volgrid6d_alloc (this, &
771  ntime=ntime, ntimerange=ntimerange, nlevel=nlevel, nvar=nvar)
772 
773 call volgrid6d_alloc_vol (this)
774 
775 if (associated(this%time)) call read_unit(this%time, lunit)
776 if (associated(this%level)) read(unit=lunit)this%level
777 if (associated(this%timerange)) read(unit=lunit)this%timerange
778 if (associated(this%var)) read(unit=lunit)this%var
779 
780 
781 !! Volumi di valori
782 
783 if (associated(this%voldati)) read(unit=lunit)this%voldati
784 
785 if (.not. present(unit)) close(unit=lunit)
786 
787 end subroutine volgrid6d_read_from_file
788 
789 
809 SUBROUTINE import_from_gridinfo(this, gridinfo, force, dup_mode, clone, &
810  isanavar)
811 TYPE(volgrid6d),INTENT(inout) :: this
812 TYPE(gridinfo_def),INTENT(in) :: gridinfo
813 LOGICAL,INTENT(in),OPTIONAL :: force
814 INTEGER,INTENT(in),OPTIONAL :: dup_mode
815 LOGICAL , INTENT(in),OPTIONAL :: clone
816 LOGICAL,INTENT(IN),OPTIONAL :: isanavar
817 
818 CHARACTER(len=255) :: type
819 INTEGER :: itime0, itimerange0, itime1, itimerange1, itime, itimerange, &
820  ilevel, ivar, ldup_mode
821 LOGICAL :: dup
822 TYPE(datetime) :: correctedtime
823 REAL,ALLOCATABLE :: tmpgrid(:,:)
824 
825 IF (PRESENT(dup_mode)) THEN
826  ldup_mode = dup_mode
827 ELSE
828  ldup_mode = 0
829 ENDIF
830 
831 call get_val(this%griddim,proj_type=type)
832 
833 #ifdef DEBUG
834 call l4f_category_log(this%category,l4f_debug,"import_from_gridinfo: "//trim(type))
835 #endif
836 
837 if (.not. c_e(type))then
838  call copy(gridinfo%griddim, this%griddim)
839 ! ho gia` fatto init, altrimenti non potrei fare la get_val(this%griddim)
840 ! per cui meglio non ripetere
841 ! call init(this,gridinfo%griddim,categoryappend)
842  CALL volgrid6d_alloc_vol(this, ini=.true.) ! decode?
843 
844 else if (.not. (this%griddim == gridinfo%griddim ))then
845 
846  CALL l4f_category_log(this%category,l4f_error, &
847  "volgrid and gridinfo grid type or size are different, gridinfo rejected")
848  CALL raise_error()
849  RETURN
850 
851 end if
852 
853 ! Cerco gli indici del campo da inserire, se non trovo metto nel primo missing
854 ilevel = index(this%level, gridinfo%level)
855 IF (ilevel == 0 .AND. optio_log(force)) THEN
856  ilevel = index(this%level, vol7d_level_miss)
857  IF (ilevel /= 0) this%level(ilevel) = gridinfo%level
858 ENDIF
859 
860 IF (ilevel == 0) THEN
861  CALL l4f_category_log(this%category,l4f_error, &
862  "volgrid6d: level not valid for volume, gridinfo rejected")
863  CALL raise_error()
864  RETURN
865 ENDIF
866 
867 IF (optio_log(isanavar)) THEN ! assign to all times and timeranges
868  itime0 = 1
869  itime1 = SIZE(this%time)
870  itimerange0 = 1
871  itimerange1 = SIZE(this%timerange)
872 ELSE ! usual case
873  correctedtime = gridinfo%time
874  IF (this%time_definition == 1) correctedtime = correctedtime + &
875  timedelta_new(sec=gridinfo%timerange%p1)
876  itime0 = index(this%time, correctedtime)
877  IF (itime0 == 0 .AND. optio_log(force)) THEN
878  itime0 = index(this%time, datetime_miss)
879  IF (itime0 /= 0) this%time(itime0) = correctedtime
880  ENDIF
881  IF (itime0 == 0) THEN
882  CALL l4f_category_log(this%category,l4f_error, &
883  "volgrid6d: time not valid for volume, gridinfo rejected")
884  CALL raise_error()
885  RETURN
886  ENDIF
887  itime1 = itime0
888 
889  itimerange0 = index(this%timerange,gridinfo%timerange)
890  IF (itimerange0 == 0 .AND. optio_log(force)) THEN
891  itimerange0 = index(this%timerange, vol7d_timerange_miss)
892  IF (itimerange0 /= 0) this%timerange(itimerange0) = gridinfo%timerange
893  ENDIF
894  IF (itimerange0 == 0) THEN
895  CALL l4f_category_log(this%category,l4f_error, &
896  "volgrid6d: timerange not valid for volume, gridinfo rejected")
897  CALL raise_error()
898  RETURN
899  ENDIF
900  itimerange1 = itimerange0
901 ENDIF
902 
903 ivar = index(this%var, gridinfo%var)
904 IF (ivar == 0 .AND. optio_log(force)) THEN
905  ivar = index(this%var, volgrid6d_var_miss)
906  IF (ivar /= 0) this%var(ivar) = gridinfo%var
907 ENDIF
908 IF (ivar == 0) THEN
909  CALL l4f_category_log(this%category,l4f_error, &
910  "volgrid6d: var not valid for volume, gridinfo rejected")
911  CALL raise_error()
912  RETURN
913 ENDIF
914 
915 DO itimerange = itimerange0, itimerange1
916  DO itime = itime0, itime1
917  IF (ASSOCIATED(this%gaid)) THEN
918  dup = .false.
919  IF (c_e(this%gaid(ilevel,itime,itimerange,ivar))) THEN
920  dup = .true.
921  CALL l4f_category_log(this%category,l4f_warn,"gaid exist: grib duplicated")
922 ! avoid memory leaks
923  IF (optio_log(clone)) CALL delete(this%gaid(ilevel,itime,itimerange,ivar))
924  ENDIF
925 
926  IF (optio_log(clone)) THEN
927  CALL copy(gridinfo%gaid, this%gaid(ilevel,itime,itimerange,ivar))
928 #ifdef DEBUG
929  CALL l4f_category_log(this%category,l4f_debug,"cloning to a new gaid")
930 #endif
931  ELSE
932  this%gaid(ilevel,itime,itimerange,ivar) = gridinfo%gaid
933  ENDIF
934 
935  IF (ASSOCIATED(this%voldati))THEN
936  IF (.NOT.dup .OR. ldup_mode == 0) THEN
937  this%voldati(:,:,ilevel,itime,itimerange,ivar) = decode_gridinfo(gridinfo)
938  ELSE IF (ldup_mode == 1) THEN
939  tmpgrid = decode_gridinfo(gridinfo) ! f2003 automatic allocation
940  WHERE(c_e(tmpgrid))
941  this%voldati(:,:,ilevel,itime,itimerange,ivar) = tmpgrid(:,:)
942  END WHERE
943  ELSE IF (ldup_mode == 2) THEN
944  WHERE(.NOT.c_e(this%voldati(:,:,ilevel,itime,itimerange,ivar)))
945  this%voldati(:,:,ilevel,itime,itimerange,ivar) = decode_gridinfo(gridinfo)
946  END WHERE
947  ENDIF
948  ENDIF
949 
950  ELSE
951  CALL l4f_category_log(this%category,l4f_error, &
952  "gaid not allocated, you probably need to call volgrid6d_alloc_vol first")
953  CALL raise_error()
954  RETURN
955  ENDIF
956  ENDDO
957 ENDDO
958 
959 
960 END SUBROUTINE import_from_gridinfo
961 
962 
967 SUBROUTINE export_to_gridinfo(this, gridinfo, itime, itimerange, ilevel, ivar, &
968  gaid_template, clone)
969 TYPE(volgrid6d),INTENT(in) :: this
970 TYPE(gridinfo_def),INTENT(inout) :: gridinfo
971 INTEGER :: itime
972 INTEGER :: itimerange
973 INTEGER :: ilevel
974 INTEGER :: ivar
975 TYPE(grid_id),INTENT(in),OPTIONAL :: gaid_template
976 LOGICAL,INTENT(in),OPTIONAL :: clone
977 
978 TYPE(grid_id) :: gaid
979 LOGICAL :: usetemplate
980 REAL,POINTER :: voldati(:,:)
981 TYPE(datetime) :: correctedtime
982 
983 #ifdef DEBUG
984 CALL l4f_category_log(this%category,l4f_debug,"export_to_gridinfo")
985 #endif
986 
987 IF (.NOT.c_e(this%gaid(ilevel,itime,itimerange,ivar))) THEN
988 #ifdef DEBUG
989  CALL l4f_category_log(this%category,l4f_debug,"empty gaid found, skipping export")
990 #endif
991  RETURN
992 ENDIF
993 
994 usetemplate = .false.
995 IF (PRESENT(gaid_template)) THEN
996  CALL copy(gaid_template, gaid)
997 #ifdef DEBUG
998  CALL l4f_category_log(this%category,l4f_debug,"template cloned to a new gaid")
999 #endif
1000  usetemplate = c_e(gaid)
1001 ENDIF
1002 
1003 IF (.NOT.usetemplate) THEN
1004  IF (optio_log(clone)) THEN
1005  CALL copy(this%gaid(ilevel,itime,itimerange,ivar), gaid)
1006 #ifdef DEBUG
1007  CALL l4f_category_log(this%category,l4f_debug,"original gaid cloned to a new one")
1008 #endif
1009  ELSE
1010  gaid = this%gaid(ilevel,itime,itimerange,ivar)
1011  ENDIF
1012 ENDIF
1013 
1014 IF (this%time_definition == 1) THEN
1015  correctedtime = this%time(itime) - &
1016  timedelta_new(sec=this%timerange(itimerange)%p1)
1017 ELSE
1018  correctedtime = this%time(itime)
1019 ENDIF
1020 
1021 CALL init(gridinfo,gaid, this%griddim, correctedtime, this%timerange(itimerange), &
1022  this%level(ilevel), this%var(ivar))
1023 
1024 ! reset the gridinfo, bad but necessary at this point for encoding the field
1025 CALL export(gridinfo%griddim, gridinfo%gaid)
1026 ! encode the field
1027 IF (ASSOCIATED(this%voldati)) THEN
1028  CALL encode_gridinfo(gridinfo, this%voldati(:,:,ilevel,itime,itimerange,ivar))
1029 ELSE IF (usetemplate) THEN ! field must be forced into template in this case
1030  NULLIFY(voldati)
1031  CALL volgrid_get_vol_2d(this, ilevel, itime, itimerange, ivar, voldati)
1032  CALL encode_gridinfo(gridinfo, voldati)
1033  DEALLOCATE(voldati)
1034 ENDIF
1035 
1036 END SUBROUTINE export_to_gridinfo
1037 
1038 
1056 SUBROUTINE import_from_gridinfovv(this, gridinfov, dup_mode, clone, decode, &
1057  time_definition, anavar, categoryappend)
1058 TYPE(volgrid6d),POINTER :: this(:)
1059 TYPE(arrayof_gridinfo),INTENT(in) :: gridinfov
1060 INTEGER,INTENT(in),OPTIONAL :: dup_mode
1061 LOGICAL , INTENT(in),OPTIONAL :: clone
1062 LOGICAL,INTENT(in),OPTIONAL :: decode
1063 INTEGER,INTENT(IN),OPTIONAL :: time_definition
1064 CHARACTER(len=*),INTENT(IN),OPTIONAL :: anavar(:)
1065 CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
1066 
1067 INTEGER :: i, j, stallo
1068 INTEGER :: ngrid, ntime, ntimerange, nlevel, nvar
1069 INTEGER :: category
1070 CHARACTER(len=512) :: a_name
1071 TYPE(datetime),ALLOCATABLE :: correctedtime(:)
1072 LOGICAL,ALLOCATABLE :: isanavar(:)
1073 TYPE(vol7d_var) :: lvar
1074 
1075 ! category temporanea (altrimenti non possiamo loggare)
1076 if (present(categoryappend))then
1077  call l4f_launcher(a_name,a_name_append=trim(subcategory)//"."//trim(categoryappend))
1078 else
1079  call l4f_launcher(a_name,a_name_append=trim(subcategory))
1080 endif
1081 category=l4f_category_get(a_name)
1082 
1083 #ifdef DEBUG
1084 call l4f_category_log(category,l4f_debug,"start import_from_gridinfovv")
1085 #endif
1086 
1087 ngrid=count_distinct(gridinfov%array(1:gridinfov%arraysize)%griddim,back=.true.)
1088 CALL l4f_category_log(category,l4f_info, t2c(ngrid)// &
1089  ' different grid definition(s) found in input data')
1090 
1091 ALLOCATE(this(ngrid),stat=stallo)
1092 IF (stallo /= 0)THEN
1093  CALL l4f_category_log(category,l4f_fatal,"allocating memory")
1094  CALL raise_fatal_error()
1095 ENDIF
1096 DO i = 1, ngrid
1097  IF (PRESENT(categoryappend))THEN
1098  CALL init(this(i), time_definition=time_definition, categoryappend=trim(categoryappend)//"-vol"//t2c(i))
1099  ELSE
1100  CALL init(this(i), time_definition=time_definition, categoryappend="vol"//t2c(i))
1101  ENDIF
1102 ENDDO
1103 
1104 this(:)%griddim=pack_distinct(gridinfov%array(1:gridinfov%arraysize)%griddim, &
1105  ngrid, back=.true.)
1106 
1107 ! mark elements as ana variables (time-independent)
1108 ALLOCATE(isanavar(gridinfov%arraysize))
1109 isanavar(:) = .false.
1110 IF (PRESENT(anavar)) THEN
1111  DO i = 1, gridinfov%arraysize
1112  DO j = 1, SIZE(anavar)
1113  lvar = convert(gridinfov%array(i)%var)
1114  IF (lvar%btable == anavar(j)) THEN
1115  isanavar(i) = .true.
1116  EXIT
1117  ENDIF
1118  ENDDO
1119  ENDDO
1120  CALL l4f_category_log(category,l4f_info,t2c(count(isanavar))//'/'// &
1121  t2c(gridinfov%arraysize)//' constant-data messages found in input data')
1122 ENDIF
1123 
1124 ! create time corrected for time_definition
1125 ALLOCATE(correctedtime(gridinfov%arraysize))
1126 correctedtime(:) = gridinfov%array(1:gridinfov%arraysize)%time
1127 IF (PRESENT(time_definition)) THEN
1128  IF (time_definition == 1) THEN
1129  DO i = 1, gridinfov%arraysize
1130  correctedtime(i) = correctedtime(i) + &
1131  timedelta_new(sec=gridinfov%array(i)%timerange%p1)
1132  ENDDO
1133  ENDIF
1134 ENDIF
1135 
1136 DO i = 1, ngrid
1137  IF (PRESENT(anavar)) THEN
1138  j = count((this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
1139  .AND. .NOT.isanavar(:))
1140  IF (j <= 0) THEN
1141  CALL l4f_category_log(category, l4f_fatal, 'grid n.'//t2c(i)// &
1142  ' has only constant data, this is not allowed')
1143  CALL l4f_category_log(category, l4f_fatal, 'please check anavar argument')
1144  CALL raise_fatal_error()
1145  ENDIF
1146  ENDIF
1147  ntime = count_distinct(correctedtime, &
1148  mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
1149  .AND. .NOT.isanavar(:), back=.true.)
1150  ntimerange = count_distinct(gridinfov%array(1:gridinfov%arraysize)%timerange, &
1151  mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
1152  .AND. .NOT.isanavar(:), back=.true.)
1153  nlevel = count_distinct(gridinfov%array(1:gridinfov%arraysize)%level, &
1154  mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim), &
1155  back=.true.)
1156  nvar = count_distinct(gridinfov%array(1:gridinfov%arraysize)%var, &
1157  mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim), &
1158  back=.true.)
1159 
1160 #ifdef DEBUG
1161  CALL l4f_category_log(this(i)%category,l4f_debug,"alloc volgrid6d index: "//t2c(i))
1162 #endif
1163 
1164  CALL volgrid6d_alloc(this(i),this(i)%griddim%dim,ntime=ntime, &
1165  ntimerange=ntimerange,nlevel=nlevel,nvar=nvar)
1166 
1167  this(i)%time = pack_distinct(correctedtime, ntime, &
1168  mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
1169  .AND. .NOT.isanavar(:), back=.true.)
1170  CALL sort(this(i)%time)
1171 
1172  this(i)%timerange = pack_distinct(gridinfov%array( &
1173  1:gridinfov%arraysize)%timerange, ntimerange, &
1174  mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
1175  .AND. .NOT.isanavar(:), back=.true.)
1176  CALL sort(this(i)%timerange)
1177 
1178  this(i)%level=pack_distinct(gridinfov%array(1:gridinfov%arraysize)%level, &
1179  nlevel,mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim), &
1180  back=.true.)
1181  CALL sort(this(i)%level)
1182 
1183  this(i)%var=pack_distinct(gridinfov%array(1:gridinfov%arraysize)%var, nvar, &
1184  mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim), &
1185  back=.true.)
1186 
1187 #ifdef DEBUG
1188  CALL l4f_category_log(this(i)%category,l4f_debug,"alloc_vol volgrid6d index: "//t2c(i))
1189 #endif
1190  CALL volgrid6d_alloc_vol(this(i), decode=decode)
1191 
1192 ENDDO
1193 
1194 DEALLOCATE(correctedtime)
1195 
1196 DO i = 1, gridinfov%arraysize
1197 
1198 #ifdef DEBUG
1199  CALL l4f_category_log(category,l4f_debug,import from gridinfov index: "//t2c(i))
1200  CALL l4f_category_log(category,L4F_INFO, &
1201  "to volgrid6d index: "//t2c(index(this%griddim, gridinfov%array(i)%griddim)))
1202 #endif
1203 
1204  CALL import(this(index(this%griddim, gridinfov%array(i)%griddim)), &
1205  gridinfov%array(i), dup_mode=dup_mode, clone=clone, isanavar=isanavar(i))
1206 
1207 ENDDO
1208 
1209 !chiudo il logger temporaneo
1210 CALL l4f_category_delete(category)
1211 
1212 END SUBROUTINE import_from_gridinfovv
1213 
1214 
1215 !> Export a \a volgrid6d object to an \a arrayof_gridinfo object.
1216 !! The multidimensional \a volgrid6d structure is serialized into a
1217 !! one-dimensional array of gridinfo_def objects, which is allocated
1218 !! to the proper size if not already allocated, or it is extended
1219 !! keeping the old data if any.
1220 SUBROUTINE export_to_gridinfov(this, gridinfov, gaid_template, clone)
1221 TYPE(volgrid6d),INTENT(inout) :: this !< volume to be exported
1222 TYPE(arrayof_gridinfo),INTENT(inout) :: gridinfov !< output array of gridinfo_def objects
1223 TYPE(grid_id),INTENT(in),OPTIONAL :: gaid_template !< \a grid_id template to be used for output data replacing the one contained in \a this
1224 LOGICAL,INTENT(in),OPTIONAL :: clone !< if provided and \c .TRUE., clone the grid_id included in \a this rather than making a shallow copy
1225 
1226 INTEGER :: i ,itime, itimerange, ilevel, ivar
1227 INTEGER :: ntime, ntimerange, nlevel, nvar
1228 TYPE(gridinfo_def) :: gridinfol
1229 
1230 #ifdef DEBUG
1231 CALL l4f_category_log(this%category,L4F_DEBUG,"start export_to_gridinfov")
1232 #endif
1233 
1234 ! this is necessary in order not to repeatedly and uselessly copy the
1235 ! same array of coordinates for each 2d grid during export, the
1236 ! side-effect is that the computed projection in this is lost
1237 CALL dealloc(this%griddim%dim)
1238 
1239 i=0
1240 ntime=size(this%time)
1241 ntimerange=size(this%timerange)
1242 nlevel=size(this%level)
1243 nvar=size(this%var)
1244 
1245 DO itime=1,ntime
1246  DO itimerange=1,ntimerange
1247  DO ilevel=1,nlevel
1248  DO ivar=1,nvar
1249 
1250  CALL init(gridinfol)
1251  CALL export(this, gridinfol, itime, itimerange, ilevel, ivar, &
1252  gaid_template=gaid_template, clone=clone)
1253  IF (c_e(gridinfol%gaid)) THEN
1254  CALL insert(gridinfov, gridinfol)
1255  ELSE
1256  CALL delete(gridinfol)
1257  ENDIF
1258 
1259  ENDDO
1260  ENDDO
1261  ENDDO
1262 ENDDO
1263 
1264 END SUBROUTINE export_to_gridinfov
1265 
1266 
1267 !> Export an array of \a volgrid6d objects to an \a arrayof_gridinfo object.
1268 !! The multidimensional \a volgrid6d structures are serialized into a
1269 !! one-dimensional array of gridinfo_def objects, which is allocated
1270 !! to the proper size if not already allocated, or it is extended
1271 !! keeping the old data if any.
1272 SUBROUTINE export_to_gridinfovv(this, gridinfov, gaid_template, clone)
1273 !, &
1274 ! categoryappend)
1275 TYPE(volgrid6d),INTENT(inout) :: this(:) !< volume array to be exported
1276 TYPE(arrayof_gridinfo),INTENT(inout) :: gridinfov !< output array of gridinfo_def objects
1277 TYPE(grid_id),INTENT(in),OPTIONAL :: gaid_template !< \a grid_id template to be used for output data replacing the one contained in \a this
1278 LOGICAL,INTENT(in),OPTIONAL :: clone !< if provided and \c .TRUE., clone the grid_id included in \a this rather than making a shallow copy
1279 
1280 INTEGER :: i
1281 
1282 DO i = 1, SIZE(this)
1283 #ifdef DEBUG
1284  CALL l4f_category_log(this(i)%category,L4F_DEBUG, &
1285  "export_to_gridinfovv grid index: "//t2c(i))
1286 #endif
1287  CALL export(this(i), gridinfov, gaid_template=gaid_template, clone=clone)
1288 ENDDO
1289 
1290 END SUBROUTINE export_to_gridinfovv
1291 
1292 
1293 Import the content of a supported file (like grib or gdal-supported!>
1294 !! format) into an array of \a volgrid6d objects. This method imports
1295 !! a set of gridded fields from a file object into a suitable number
1296 !! of \a volgrid6d objects. The data are imported by creating a
1297 !! temporary \a gridinfo object, importing it into the \a volgrid6d
1298 !! object cloning the gaid's and then destroying the gridinfo, so it
1299 !! works similarly to volgrid6d_class::import_from_gridinfovv() method.
1300 !! For a detailed explanation of the \a anavar argument, see the
1301 !! documentation of volgrid6d_class::import_from_gridinfovv() method.
1302 SUBROUTINE volgrid6d_import_from_file(this, filename, dup_mode, decode, &
1303  time_definition, anavar, categoryappend)
1304 importTYPE(volgrid6d),POINTER :: this(:) !< object in which to
1305 importCHARACTER(len=*),INTENT(in) :: filename !< name of file from which to
1306 INTEGER,INTENT(in),OPTIONAL :: dup_mode !< determines the behavior in case of duplicate metadata: if \a dup_mode is not provided or 0, a duplicate field overwrites, if \a dup_mode is 1, duplicate fields are merged with priority to the last
1307 LOGICAL,INTENT(in),OPTIONAL :: decode !< if provided and \a .FALSE. the data volume in the elements of \a this is not allocated and successive work will be performed on grid_id's
1308 INTEGER,INTENT(IN),OPTIONAL :: time_definition !< 0=time is reference time; 1=time is validity time
1309 CHARACTER(len=*),INTENT(IN),OPTIONAL :: anavar(:) !< list of variables (B-table code equivalent) to be treated as time-independent data
1310 character(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
1311 
1312 TYPE(arrayof_gridinfo) :: gridinfo
1313 INTEGER :: category
1314 CHARACTER(len=512) :: a_name
1315 
1316 NULLIFY(this)
1317 
1318 IF (PRESENT(categoryappend))THEN
1319  CALL l4f_launcher(a_name,a_name_append= &
1320  TRIM(subcategory)//"."//TRIM(categoryappend))
1321 ELSE
1322  CALL l4f_launcher(a_name,a_name_append=TRIM(subcategory))
1323 ENDIF
1324 category=l4f_category_get(a_name)
1325 
1326 CALL import(gridinfo, filename=filename, categoryappend=categoryappend)
1327 
1328 IF (gridinfo%arraysize > 0) THEN
1329 
1330  CALL import(this, gridinfo, dup_mode=dup_mode, clone=.TRUE., decode=decode, &
1331  time_definition=time_definition, anavar=anavar, &
1332  categoryappend=categoryappend)
1333 
1334  CALL l4f_category_log(category,L4F_INFO,"deleting gridinfo")
1335  CALL delete(gridinfo)
1336 
1337 ELSE
1338  CALL l4f_category_log(category,L4F_INFO,"file does not contain gridded data")
1339 ENDIF
1340 
1341 ! close logger
1342 CALL l4f_category_delete(category)
1343 
1344 END SUBROUTINE volgrid6d_import_from_file
1345 
1346 
1347 !> High level method for exporting a volume array to file.
1348 !! All the information contained into an array of \a volgrid6d
1349 !! objects, i.e. dimension descriptors and data, is exported to a file
1350 !! using the proper output driver (typically grib_api for grib
1351 !! format). If a template is provided, it will determine the
1352 !! characteristic of the output file, otherwise the \a grid_id
1353 !! descriptors contained in the volgrid6d object will be used
1354 SUBROUTINE volgrid6d_export_to_file(this, filename, gaid_template, categoryappend)
1355 TYPE(volgrid6d) :: this(:) !< volume(s) to be exported
1356 CHARACTER(len=*),INTENT(in) :: filename !< output file name
1357 TYPE(grid_id),INTENT(in),OPTIONAL :: gaid_template !< template for the output file, if provided the grid_id information stored in the volgrid6d objects will be ignored
1358 character(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
1359 
1360 TYPE(arrayof_gridinfo) :: gridinfo
1361 INTEGER :: category
1362 CHARACTER(len=512) :: a_name
1363 
1364 IF (PRESENT(categoryappend)) THEN
1365  CALL l4f_launcher(a_name,a_name_append=TRIM(subcategory)//"."//TRIM(categoryappend))
1366 ELSE
1367  CALL l4f_launcher(a_name,a_name_append=TRIM(subcategory))
1368 ENDIF
1369 category=l4f_category_get(a_name)
1370 
1371 #ifdef DEBUG
1372 CALL l4f_category_log(category,L4F_DEBUG,"start export to file")
1373 #endif
1374 
1375 CALL l4f_category_log(category,L4F_INFO,"writing volgrid6d to grib file: "//TRIM(filename))
1376 
1377 !IF (ASSOCIATED(this)) THEN
1378  CALL export(this, gridinfo, gaid_template=gaid_template, clone=.TRUE.)
1379  IF (gridinfo%arraysize > 0) THEN
1380  CALL export(gridinfo, filename)
1381  CALL delete(gridinfo)
1382  ENDIF
1383 !ELSE
1384 ! CALL l4f_category_log(category,L4F_INFO,"volume volgrid6d is not associated")
1385 !ENDIF
1386 
1387 ! close logger
1388 CALL l4f_category_delete(category)
1389 
1390 END SUBROUTINE volgrid6d_export_to_file
1391 
1392 
1393 !> Array destructor for \a volgrid6d class.
1394 !! Delete an array of \a volgrid6d objects and deallocate the array
1395 !! itself.
1396 SUBROUTINE volgrid6dv_delete(this)
1397 TYPE(volgrid6d),POINTER :: this(:) !< vector of volgrid6d object
1398 
1399 INTEGER :: i
1400 
1401 IF (ASSOCIATED(this)) THEN
1402  DO i = 1, SIZE(this)
1403 #ifdef DEBUG
1404  CALL l4f_category_log(this(i)%category,L4F_DEBUG, &
1405  "delete volgrid6d vector index: "//TRIM(to_char(i)))
1406 #endif
1407  CALL delete(this(i))
1408  ENDDO
1409  DEALLOCATE(this)
1410 ENDIF
1411 
1412 END SUBROUTINE volgrid6dv_delete
1413 
1415 ! Internal method for performing grid to grid computations
1416 SUBROUTINE volgrid6d_transform_compute(this, volgrid6d_in, volgrid6d_out, &
1417  lev_out, var_coord_vol, clone)
1418 TYPE(grid_transform),INTENT(in) :: this ! oggetto di trasformazione per il grigliato
1419 type(volgrid6d), INTENT(in) :: volgrid6d_in ! oggetto da trasformare
1420 type(volgrid6d), INTENT(inout) :: volgrid6d_out ! oggetto trasformato; deve essere completo (init, alloc, alloc_vol)
1421 TYPE(vol7d_level),INTENT(in),OPTIONAL :: lev_out(:) ! vol7d_level object defining target vertical grid, for vertical interpolations
1422 INTEGER,INTENT(in),OPTIONAL :: var_coord_vol ! index of variable defining vertical coordinate values in input volume
1423 LOGICAL,INTENT(in),OPTIONAL :: clone ! se fornito e \c .TRUE., clona i gaid da volgrid6d_in a volgrid6d_out
1424 
1425 INTEGER :: ntime, ntimerange, inlevel, onlevel, nvar, &
1426  itime, itimerange, ilevel, ivar, levshift, levused, lvar_coord_vol, spos
1427 REAL,POINTER :: voldatiin(:,:,:), voldatiout(:,:,:), coord_3d_in(:,:,:)
1428 TYPE(vol7d_level) :: output_levtype
1429 
1430 
1431 #ifdef DEBUG
1432 call l4f_category_log(volgrid6d_in%category,L4F_DEBUG,"start volgrid6d_transform_compute")
1433 #endif
1434 
1435 ntime=0
1436 ntimerange=0
1437 inlevel=0
1438 onlevel=0
1439 nvar=0
1440 lvar_coord_vol = optio_i(var_coord_vol)
1441 
1442 if (associated(volgrid6d_in%time))then
1443  ntime=size(volgrid6d_in%time)
1444  volgrid6d_out%time=volgrid6d_in%time
1445 end if
1446 
1447 if (associated(volgrid6d_in%timerange))then
1448  ntimerange=size(volgrid6d_in%timerange)
1449  volgrid6d_out%timerange=volgrid6d_in%timerange
1450 end if
1451 
1452 IF (ASSOCIATED(volgrid6d_in%level))THEN
1453  inlevel=SIZE(volgrid6d_in%level)
1454 ENDIF
1455 IF (PRESENT(lev_out)) THEN
1456  onlevel=SIZE(lev_out)
1457  volgrid6d_out%level=lev_out
1458 ELSE IF (ASSOCIATED(volgrid6d_in%level))THEN
1459  onlevel=SIZE(volgrid6d_in%level)
1460  volgrid6d_out%level=volgrid6d_in%level
1461 ENDIF
1462 
1463 if (associated(volgrid6d_in%var))then
1464  nvar=size(volgrid6d_in%var)
1465  volgrid6d_out%var=volgrid6d_in%var
1466 end if
1467 ! allocate once for speed
1468 .NOT.IF (ASSOCIATED(volgrid6d_in%voldati)) THEN
1469  ALLOCATE(voldatiin(volgrid6d_in%griddim%dim%nx, volgrid6d_in%griddim%dim%ny, &
1470  inlevel))
1471 ENDIF
1472 .NOT.IF (ASSOCIATED(volgrid6d_out%voldati)) THEN
1473  ALLOCATE(voldatiout(volgrid6d_out%griddim%dim%nx, volgrid6d_out%griddim%dim%ny, &
1474  onlevel))
1475 ENDIF
1476 
1477 CALL get_val(this, levshift=levshift, levused=levused)
1478 spos = imiss
1479 IF (c_e(lvar_coord_vol)) THEN
1480  CALL get_val(this%trans, output_levtype=output_levtype)
1481 .OR. IF (output_levtype%level1 == 103 output_levtype%level1 == 108) THEN
1482  spos = firsttrue(volgrid6d_in%level(:) == vol7d_level_new(1))
1483  IF (spos == 0) THEN
1484  CALL l4f_category_log(volgrid6d_in%category, L4F_ERROR, &
1485  'output level '//t2c(output_levtype%level1)// &
1486  ' requested, but height/press of surface not provided in volume')
1487  ENDIF
1488 .NOT..AND..NOT. IF (c_e(levshift) c_e(levused)) THEN
1489  CALL l4f_category_log(volgrid6d_in%category, L4F_ERROR, &
1490  'internal inconsistence, levshift and levused undefined when they should be')
1491  ENDIF
1492  ENDIF
1493 ENDIF
1494 
1495 DO ivar=1,nvar
1496 ! IF (c_e(var_coord_vol)) THEN
1497 ! IF (ivar == var_coord_vol) CYCLE ! skip coordinate variable in output
1498 ! ENDIF
1499  DO itimerange=1,ntimerange
1500  DO itime=1,ntime
1501 ! skip empty columns where possible, improve
1502 .AND. IF (c_e(levshift) c_e(levused)) THEN
1503 .NOT. IF (ANY(c_e( &
1504  volgrid6d_in%gaid(levshift+1:levshift+levused,itime,itimerange,ivar) &
1505  ))) CYCLE
1506  ENDIF
1507  DO ilevel = 1, MIN(inlevel,onlevel)
1508 ! if present gaid copy it
1509 .AND..NOT. IF (c_e(volgrid6d_in%gaid(ilevel,itime,itimerange,ivar)) &
1510  c_e(volgrid6d_out%gaid(ilevel,itime,itimerange,ivar))) THEN
1511 
1512  IF (optio_log(clone)) THEN
1513  CALL copy(volgrid6d_in%gaid(ilevel,itime,itimerange,ivar),&
1514  volgrid6d_out%gaid(ilevel,itime,itimerange,ivar))
1515 #ifdef DEBUG
1516  CALL l4f_category_log(volgrid6d_in%category,L4F_DEBUG, &
1517  "cloning gaid, level "//t2c(ilevel))
1518 #endif
1519  ELSE
1520  volgrid6d_out%gaid(ilevel,itime,itimerange,ivar) = &
1521  volgrid6d_in%gaid(ilevel,itime,itimerange,ivar)
1522  ENDIF
1523  ENDIF
1524  ENDDO
1525 ! if out level are more, we have to clone anyway
1526  DO ilevel = MIN(inlevel,onlevel) + 1, onlevel
1527 .AND..NOT. IF (c_e(volgrid6d_in%gaid(inlevel,itime,itimerange,ivar)) &
1528  c_e(volgrid6d_out%gaid(ilevel,itime,itimerange,ivar))) then
1529 
1530  CALL copy(volgrid6d_in%gaid(inlevel,itime,itimerange,ivar),&
1531  volgrid6d_out%gaid(ilevel,itime,itimerange,ivar))
1532 #ifdef DEBUG
1533  CALL l4f_category_log(volgrid6d_in%category,L4F_DEBUG, &
1534  "forced cloning gaid, level "//t2c(inlevel)//"->"//t2c(ilevel))
1535 #endif
1536  ENDIF
1537  ENDDO
1538 
1539  IF (c_e(lvar_coord_vol)) THEN
1540  NULLIFY(coord_3d_in)
1541  CALL volgrid_get_vol_3d(volgrid6d_in, itime, itimerange, lvar_coord_vol, &
1542  coord_3d_in)
1543  IF (c_e(spos)) THEN ! compute difference wrt surface coordinate
1544  IF (spos == 0) THEN ! error condition, set all to missing and goodnight
1545  coord_3d_in(:,:,levshift+1:levshift+levused) = rmiss
1546  ELSE
1547  DO ilevel = levshift+1, levshift+levused
1548 .AND. WHERE(c_e(coord_3d_in(:,:,ilevel)) c_e(coord_3d_in(:,:,spos)))
1549  coord_3d_in(:,:,ilevel) = coord_3d_in(:,:,ilevel) - &
1550  coord_3d_in(:,:,spos)
1551  ELSEWHERE
1552  coord_3d_in(:,:,ilevel) = rmiss
1553  END WHERE
1554  ENDDO
1555  ENDIF
1556  ENDIF
1557  ENDIF
1558  CALL volgrid_get_vol_3d(volgrid6d_in, itime, itimerange, ivar, &
1559  voldatiin)
1560  IF (ASSOCIATED(volgrid6d_out%voldati)) & ! improve!!!!
1561  CALL volgrid_get_vol_3d(volgrid6d_out, itime, itimerange, ivar, &
1562  voldatiout)
1563  IF (c_e(lvar_coord_vol)) THEN
1564  CALL compute(this, voldatiin, voldatiout, convert(volgrid6d_in%var(ivar)), &
1565  coord_3d_in(:,:,levshift+1:levshift+levused)) ! subset coord_3d_in
1566  ELSE
1567  CALL compute(this, voldatiin, voldatiout, convert(volgrid6d_in%var(ivar)))
1568  ENDIF
1569  CALL volgrid_set_vol_3d(volgrid6d_out, itime, itimerange, ivar, &
1570  voldatiout)
1571  ENDDO
1572  ENDDO
1573 ENDDO
1574 
1575 IF (c_e(lvar_coord_vol)) THEN
1576  DEALLOCATE(coord_3d_in)
1577 ENDIF
1578 .NOT.IF (ASSOCIATED(volgrid6d_in%voldati)) THEN
1579  DEALLOCATE(voldatiin)
1580 ENDIF
1581 .NOT.IF (ASSOCIATED(volgrid6d_out%voldati)) THEN
1582  DEALLOCATE(voldatiout)
1583 ENDIF
1584 
1585 
1586 END SUBROUTINE volgrid6d_transform_compute
1587 
1588 
1589 !> Performs the specified abstract transformation on the data provided.
1590 !! The abstract transformation is specified by \a this parameter; the
1591 !! corresponding specifical transformation (\a grid_transform object)
1592 !! is created and destroyed internally. The output transformed object
1593 !! is created internally and it does not require preliminary
1594 !! initialisation.
1595 SUBROUTINE volgrid6d_transform(this, griddim, volgrid6d_in, volgrid6d_out, &
1596  lev_out, volgrid6d_coord_in, maskgrid, maskbounds, clone, decode, categoryappend)
1597 TYPE(transform_def),INTENT(in) :: this !< object specifying the abstract transformation
1598 TYPE(griddim_def),INTENT(in),OPTIONAL :: griddim !< griddim specifying the output grid (required by most transformation types)
1599 TYPE(volgrid6d),INTENT(inout) :: volgrid6d_in !< object to be transformed, it is not modified, despite the INTENT(inout)
1600 TYPE(volgrid6d),INTENT(out) :: volgrid6d_out !< transformed object, it does not require initialisation
1601 TYPE(vol7d_level),INTENT(in),OPTIONAL,TARGET :: lev_out(:) !< vol7d_level object defining target vertical grid, for vertical interpolations
1602 TYPE(volgrid6d),INTENT(in),OPTIONAL :: volgrid6d_coord_in !< object providing time constant input vertical coordinate for some kind of vertical interpolations
1603 REAL,INTENT(in),OPTIONAL :: maskgrid(:,:) !< 2D field to be used for defining subareas according to its values, it must have the same shape as the field to be interpolated (for transformation subtype 'maskfill')
1604 REAL,INTENT(in),OPTIONAL :: maskbounds(:) !< array of boundary values for defining a subset of valid points where the values of \a maskgrid are within the first and last value of \a maskbounds (for transformation type 'metamorphosis:maskfill')
1605 LOGICAL,INTENT(in),OPTIONAL :: clone !< if provided and \a .TRUE. , clone the \a gaid's from \a volgrid6d_in to \a volgrid6d_out
1606 LOGICAL,INTENT(in),OPTIONAL :: decode !< determine whether the data in \a volgrid6d_out should be decoded or remain coded in gaid, if not provided, the decode status is taken from \a volgrid6d_in
1607 CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
1608 
1609 TYPE(grid_transform) :: grid_trans
1610 TYPE(vol7d_level),POINTER :: llev_out(:)
1611 TYPE(vol7d_level) :: input_levtype, output_levtype
1612 TYPE(vol7d_var) :: vcoord_var
1613 INTEGER :: i, k, ntime, ntimerange, nlevel, nvar, var_coord_in, var_coord_vol, &
1614  cf_out, nxc, nyc, nxi, nyi, i3, i4, i5, i6, &
1615  ulstart, ulend, spos
1616 REAL,ALLOCATABLE :: coord_3d_in(:,:,:)
1617 TYPE(geo_proj) :: proj_in, proj_out
1618 CHARACTER(len=80) :: trans_type
1619 LOGICAL :: ldecode
1620 LOGICAL,ALLOCATABLE :: mask_in(:)
1621 
1622 #ifdef DEBUG
1623 call l4f_category_log(volgrid6d_in%category, L4F_DEBUG, "start volgrid6d_transform")
1624 #endif
1625 
1626 ntime=0
1627 ntimerange=0
1628 nlevel=0
1629 nvar=0
1630 
1631 if (associated(volgrid6d_in%time)) ntime=size(volgrid6d_in%time)
1632 if (associated(volgrid6d_in%timerange)) ntimerange=size(volgrid6d_in%timerange)
1633 if (associated(volgrid6d_in%level)) nlevel=size(volgrid6d_in%level)
1634 if (associated(volgrid6d_in%var)) nvar=size(volgrid6d_in%var)
1635 
1636 .OR..OR..OR.IF (ntime == 0 ntimerange == 0 nlevel == 0 nvar == 0) THEN
1637  CALL l4f_category_log(volgrid6d_in%category, L4F_ERROR, &
1638  "trying to transform an incomplete volgrid6d object, ntime="//t2c(ntime)// &
1639  ' ntimerange='//t2c(ntimerange)//' nlevel='//t2c(nlevel)//' nvar='//t2c(nvar))
1640  CALL init(volgrid6d_out) ! initialize to empty
1641  CALL raise_error()
1642  RETURN
1643 ENDIF
1644 
1645 CALL get_val(this, trans_type=trans_type)
1646 
1647 ! store desired output component flag and unrotate if necessary
1648 cf_out = imiss
1649 .AND..OR.IF (PRESENT(griddim) (trans_type == 'inter' trans_type == 'boxinter' &
1650 .OR. trans_type == 'stencilinter')) THEN ! improve condition!!
1651  CALL get_val(volgrid6d_in%griddim, proj=proj_in)
1652  CALL get_val(griddim, component_flag=cf_out, proj=proj_out)
1653 ! if different projections wind components must be referred to geographical system
1654  IF (proj_in /= proj_out) CALL vg6d_wind_unrot(volgrid6d_in)
1655 ELSE IF (PRESENT(griddim)) THEN ! just get component_flag, the rest is rubbish
1656  CALL get_val(griddim, component_flag=cf_out)
1657 ENDIF
1658 
1659 
1660 var_coord_in = imiss
1661 var_coord_vol = imiss
1662 IF (trans_type == 'vertint') THEN
1663  IF (PRESENT(lev_out)) THEN
1664 
1665 ! if volgrid6d_coord_in provided and allocated, check that it fits
1666  IF (PRESENT(volgrid6d_coord_in)) THEN
1667  IF (ASSOCIATED(volgrid6d_coord_in%voldati)) THEN
1668 
1669 ! strictly 1 time and 1 timerange
1670 .OR. IF (SIZE(volgrid6d_coord_in%voldati,4) /= 1 &
1671  SIZE(volgrid6d_coord_in%voldati,5) /= 1) THEN
1672  CALL l4f_category_log(volgrid6d_in%category, L4F_ERROR, &
1673  'volume providing constant input vertical coordinate must have &
1674  &only 1 time and 1 timerange')
1675  CALL init(volgrid6d_out) ! initialize to empty
1676  CALL raise_error()
1677  RETURN
1678  ENDIF
1679 
1680 ! search for variable providing vertical coordinate
1681  CALL get_val(this, output_levtype=output_levtype)
1682  vcoord_var = vol7d_var_new(vol7d_level_to_var(output_levtype))
1683 .NOT. IF (c_e(vcoord_var)) THEN
1684  CALL l4f_category_log(volgrid6d_in%category, L4F_ERROR, &
1685  'requested output level type '//t2c(output_levtype%level1)// &
1686  ' does not correspond to any known physical variable for &
1687  &providing vertical coordinate')
1688  CALL init(volgrid6d_out) ! initialize to empty
1689  CALL raise_error()
1690  RETURN
1691  ENDIF
1692 
1693  DO i = 1, SIZE(volgrid6d_coord_in%var)
1694  IF (convert(volgrid6d_coord_in%var(i)) == vcoord_var) THEN
1695  var_coord_in = i
1696  EXIT
1697  ENDIF
1698  ENDDO
1699 
1700 .NOT. IF (c_e(var_coord_in)) THEN
1701  CALL l4f_category_log(volgrid6d_in%category, L4F_ERROR, &
1702  'volume providing constant input vertical coordinate contains no &
1703  &variables matching output level type '//t2c(output_levtype%level1))
1704  CALL init(volgrid6d_out) ! initialize to empty
1705  CALL raise_error()
1706  RETURN
1707  ENDIF
1708  CALL l4f_category_log(volgrid6d_in%category, L4F_INFO, &
1709  'Coordinate for vertint found in coord volume at position '// &
1710  t2c(var_coord_in))
1711 
1712 ! check horizontal grid
1713 ! this is too strict (component flag and so on)
1714 ! IF (volgrid6d_coord_in%griddim /= volgrid6d_in%griddim) THEN
1715  CALL get_val(volgrid6d_coord_in%griddim, nx=nxc, ny=nyc)
1716  CALL get_val(volgrid6d_in%griddim, nx=nxi, ny=nyi)
1717 .OR. IF (nxc /= nxi nyc /= nyi) THEN
1718  CALL l4f_category_log(volgrid6d_in%category, L4F_ERROR, &
1719  'volume providing constant input vertical coordinate must have &
1720  &the same grid as the input')
1721  CALL l4f_category_log(volgrid6d_in%category, L4F_ERROR, &
1722  'vertical coordinate: '//t2c(nxc)//'x'//t2c(nyc)// &
1723  ', input volume: '//t2c(nxi)//'x'//t2c(nyi))
1724  CALL init(volgrid6d_out) ! initialize to empty
1725  CALL raise_error()
1726  RETURN
1727  ENDIF
1728 
1729 ! check vertical coordinate system
1730  CALL get_val(this, input_levtype=input_levtype)
1731  mask_in = & ! implicit allocation
1732 .AND. (volgrid6d_coord_in%level(:)%level1 == input_levtype%level1) &
1733  (volgrid6d_coord_in%level(:)%level2 == input_levtype%level2)
1734  ulstart = firsttrue(mask_in)
1735  ulend = lasttrue(mask_in)
1736 .OR. IF (ulstart == 0 ulend == 0) THEN
1737  CALL l4f_category_log(volgrid6d_in%category, L4F_ERROR, &
1738  'coordinate file does not contain levels of type '// &
1739  t2c(input_levtype%level1)//'/'//t2c(input_levtype%level2)// &
1740  ' specified for input data')
1741  CALL init(volgrid6d_out) ! initialize to empty
1742  CALL raise_error()
1743  RETURN
1744  ENDIF
1745 
1746  coord_3d_in = volgrid6d_coord_in%voldati(:,:,ulstart:ulend,1,1,var_coord_in) ! implicit allocation
1747 ! special case
1748 .OR. IF (output_levtype%level1 == 103 &
1749  output_levtype%level1 == 108) THEN ! surface coordinate needed
1750  spos = firsttrue(volgrid6d_coord_in%level(:) == vol7d_level_new(1))
1751  IF (spos == 0) THEN
1752  CALL l4f_category_log(volgrid6d_in%category, L4F_ERROR, &
1753  'output level '//t2c(output_levtype%level1)// &
1754  ' requested, but height/press of surface not provided in coordinate file')
1755  CALL init(volgrid6d_out) ! initialize to empty
1756  CALL raise_error()
1757  RETURN
1758  ENDIF
1759  DO k = 1, SIZE(coord_3d_in,3)
1760 .AND. WHERE(c_e(coord_3d_in(:,:,k)) &
1761  c_e(volgrid6d_coord_in%voldati(:,:,spos,1,1,var_coord_in)))
1762  coord_3d_in(:,:,k) = coord_3d_in(:,:,k) - &
1763  volgrid6d_coord_in%voldati(:,:,spos,1,1,var_coord_in)
1764  ELSEWHERE
1765  coord_3d_in(:,:,k) = rmiss
1766  END WHERE
1767  ENDDO
1768  ENDIF
1769 
1770  ENDIF
1771  ENDIF
1772 
1773 .NOT. IF (c_e(var_coord_in)) THEN ! search for coordinate within volume
1774 ! search for variable providing vertical coordinate
1775  CALL get_val(this, output_levtype=output_levtype)
1776  vcoord_var = vol7d_var_new(vol7d_level_to_var(output_levtype))
1777  IF (c_e(vcoord_var)) THEN
1778  DO i = 1, SIZE(volgrid6d_in%var)
1779  IF (convert(volgrid6d_in%var(i)) == vcoord_var) THEN
1780  var_coord_vol = i
1781  EXIT
1782  ENDIF
1783  ENDDO
1784 
1785  IF (c_e(var_coord_vol)) THEN
1786  CALL l4f_category_log(volgrid6d_in%category, L4F_INFO, &
1787  'Coordinate for vertint found in input volume at position '// &
1788  t2c(var_coord_vol))
1789  ENDIF
1790 
1791  ENDIF
1792  ENDIF
1793 
1794  CALL init(volgrid6d_out, griddim=volgrid6d_in%griddim, &
1795  time_definition=volgrid6d_in%time_definition, categoryappend=categoryappend)
1796  IF (c_e(var_coord_in)) THEN
1797  CALL init(grid_trans, this, lev_in=volgrid6d_in%level, lev_out=lev_out, &
1798  coord_3d_in=coord_3d_in, categoryappend=categoryappend)
1799  ELSE
1800  CALL init(grid_trans, this, lev_in=volgrid6d_in%level, lev_out=lev_out, &
1801  categoryappend=categoryappend)
1802  ENDIF
1803 
1804  CALL get_val(grid_trans, output_level_auto=llev_out) ! get levels if auto-generated
1805 .NOT. IF (ASSOCIATED(llev_out)) llev_out => lev_out
1806  nlevel = SIZE(llev_out)
1807  ELSE
1808  CALL l4f_category_log(volgrid6d_in%category, L4F_ERROR, &
1809  'volgrid6d_transform: vertint requested but lev_out not provided')
1810  CALL init(volgrid6d_out) ! initialize to empty
1811  CALL raise_error()
1812  RETURN
1813  ENDIF
1814 
1815 ELSE
1816  CALL init(volgrid6d_out, griddim=griddim, &
1817  time_definition=volgrid6d_in%time_definition, categoryappend=categoryappend)
1818  CALL init(grid_trans, this, in=volgrid6d_in%griddim, out=volgrid6d_out%griddim, &
1819  maskgrid=maskgrid, maskbounds=maskbounds, categoryappend=categoryappend)
1820 ENDIF
1821 
1822 
1823 IF (c_e(grid_trans)) THEN ! transformation is valid
1824 
1825  CALL volgrid6d_alloc(volgrid6d_out, ntime=ntime, nlevel=nlevel, &
1826  ntimerange=ntimerange, nvar=nvar)
1827 
1828  IF (PRESENT(decode)) THEN ! explicitly set decode status
1829  ldecode = decode
1830  ELSE ! take it from input
1831  ldecode = ASSOCIATED(volgrid6d_in%voldati)
1832  ENDIF
1833 ! force decode if gaid is readonly
1834  decode_loop: DO i6 = 1,nvar
1835  DO i5 = 1, ntimerange
1836  DO i4 = 1, ntime
1837  DO i3 = 1, nlevel
1838  IF (c_e(volgrid6d_in%gaid(i3,i4,i5,i6))) THEN
1839 .OR. ldecode = ldecode grid_id_readonly(volgrid6d_in%gaid(i3,i4,i5,i6))
1840  EXIT decode_loop
1841  ENDIF
1842  ENDDO
1843  ENDDO
1844  ENDDO
1845  ENDDO decode_loop
1846 
1847  IF (PRESENT(decode)) THEN
1848 .NEQV. IF (ldecodedecode) THEN
1849  CALL l4f_category_log(volgrid6d_in%category, L4F_WARN, &
1850  'volgrid6d_transform: decode status forced to .TRUE. because driver does not allow copy')
1851  ENDIF
1852  ENDIF
1853 
1854  CALL volgrid6d_alloc_vol(volgrid6d_out, decode=ldecode)
1855 
1856 !ensure unproj was called
1857 !call griddim_unproj(volgrid6d_out%griddim)
1858 
1859  IF (trans_type == 'vertint') THEN
1860 #ifdef DEBUG
1861  CALL l4f_category_log(volgrid6d_in%category, L4F_DEBUG, &
1862  "volgrid6d_transform: vertint to "//t2c(nlevel)//" levels")
1863 #endif
1864  CALL compute(grid_trans, volgrid6d_in, volgrid6d_out, lev_out=llev_out, &
1865  var_coord_vol=var_coord_vol, clone=clone)
1866  ELSE
1867  CALL compute(grid_trans, volgrid6d_in, volgrid6d_out, clone=clone)
1868  ENDIF
1869 
1870  IF (cf_out == 0) THEN ! unrotated components are desired
1871  CALL wind_unrot(volgrid6d_out) ! unrotate if necessary
1872  ELSE IF (cf_out == 1) THEN ! rotated components are desired
1873  CALL wind_rot(volgrid6d_out) ! rotate if necessary
1874  ENDIF
1875 
1876 ELSE
1877 ! should log with grid_trans%category, but it is private
1878  CALL l4f_category_log(volgrid6d_in%category, L4F_ERROR, &
1879  'volgrid6d_transform: transformation not valid')
1880  CALL raise_error()
1881 ENDIF
1882 
1883 CALL delete (grid_trans)
1884 
1885 END SUBROUTINE volgrid6d_transform
1886 
1887 
1888 !> Performs the specified abstract transformation on the arrays of
1889 !! data provided. The abstract transformation is specified by \a this
1890 !! parameter; the corresponding specifical transformation (\a
1891 !! grid_transform object) is created and destroyed internally. The
1892 !! output transformed object is created internally and it does not
1893 !! require preliminary initialisation. According to the input data and
1894 !! to the transformation type, the output array may have of one or
1895 !! more \a volgrid6d elements on different grids.
1896 SUBROUTINE volgrid6dv_transform(this, griddim, volgrid6d_in, volgrid6d_out, &
1897  lev_out, volgrid6d_coord_in, maskgrid, maskbounds, clone, decode, categoryappend)
1898 TYPE(transform_def),INTENT(in) :: this !< object specifying the abstract transformation
1899 TYPE(griddim_def),INTENT(in),OPTIONAL :: griddim !< griddim specifying the output grid (required by most transformation types)
1900 TYPE(volgrid6d),INTENT(inout) :: volgrid6d_in(:) !< object to be transformed, it is an array of volgrid6d objects, each of which will be transformed, it is not modified, despite the INTENT(inout)
1901 TYPE(volgrid6d),POINTER :: volgrid6d_out(:) !< transformed object, it is a non associated pointer to an array of volgrid6d objects which will be allocated by the method
1902 TYPE(vol7d_level),INTENT(in),OPTIONAL :: lev_out(:) !< vol7d_level object defining target vertical grid
1903 TYPE(volgrid6d),INTENT(in),OPTIONAL :: volgrid6d_coord_in !< object providing time constant input vertical coordinate for some kind of vertical interpolations
1904 REAL,INTENT(in),OPTIONAL :: maskgrid(:,:) !< 2D field to be used for defining subareas according to its values, it must have the same shape as the field to be interpolated (for transformation subtype 'maskfill')
1905 REAL,INTENT(in),OPTIONAL :: maskbounds(:) !< array of boundary values for defining a subset of valid points where the values of \a maskgrid are within the first and last value of \a maskbounds (for transformation type 'metamorphosis:maskfill')
1906 LOGICAL,INTENT(in),OPTIONAL :: clone !< if provided and \a .TRUE. , clone the \a gaid's from \a volgrid6d_in to \a volgrid6d_out
1907 LOGICAL,INTENT(in),OPTIONAL :: decode !< if provided and \a .FALSE. the data volume is not allocated, but work is performed on grid_id's
1908 CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
1909 
1910 INTEGER :: i, stallo
1911 
1912 
1913 allocate(volgrid6d_out(size(volgrid6d_in)),stat=stallo)
1914 if (stallo /= 0)then
1915  call l4f_log(L4F_FATAL,"allocating memory")
1916  call raise_fatal_error()
1917 end if
1918 
1919 do i=1,size(volgrid6d_in)
1920  call transform(this, griddim, volgrid6d_in(i), volgrid6d_out(i), &
1921  lev_out=lev_out, volgrid6d_coord_in=volgrid6d_coord_in, &
1922  maskgrid=maskgrid, maskbounds=maskbounds, &
1923  clone=clone, decode=decode, categoryappend=categoryappend)
1924 end do
1925 
1926 END SUBROUTINE volgrid6dv_transform
1927 
1928 
1929 ! Internal method for performing grid to sparse point computations
1930 SUBROUTINE volgrid6d_v7d_transform_compute(this, volgrid6d_in, vol7d_out, &
1931  networkname, noconvert)
1932 TYPE(grid_transform),INTENT(in) :: this ! oggetto di trasformazione per grigliato
1933 type(volgrid6d), INTENT(in) :: volgrid6d_in ! oggetto da trasformare
1934 type(vol7d), INTENT(inout) :: vol7d_out ! oggetto trasformato
1935 CHARACTER(len=*),OPTIONAL,INTENT(in) :: networkname ! imposta il network in vol7d_out (default='generic')
1936 LOGICAL,OPTIONAL,INTENT(in) :: noconvert !< do not try to match variable and convert values during transform
1937 
1938 INTEGER :: nntime, nana, ntime, ntimerange, nlevel, nvar, stallo
1939 INTEGER :: itime, itimerange, ivar, inetwork
1940 REAL,ALLOCATABLE :: voldatir_out(:,:,:)
1941 TYPE(conv_func),POINTER :: c_func(:)
1942 TYPE(datetime),ALLOCATABLE :: validitytime(:,:)
1943 REAL,POINTER :: voldatiin(:,:,:)
1944 
1945 #ifdef DEBUG
1946 call l4f_category_log(volgrid6d_in%category,L4F_DEBUG,"start volgrid6d_v7d_transform_compute")
1947 #endif
1948 
1949 ntime=0
1950 ntimerange=0
1951 nlevel=0
1952 nvar=0
1953 NULLIFY(c_func)
1954 
1955 if (present(networkname))then
1956  call init(vol7d_out%network(1),name=networkname)
1957 else
1958  call init(vol7d_out%network(1),name='generic')
1959 end if
1960 
1961 if (associated(volgrid6d_in%timerange))then
1962  ntimerange=size(volgrid6d_in%timerange)
1963  vol7d_out%timerange=volgrid6d_in%timerange
1964 end if
1965 
1966 if (associated(volgrid6d_in%time))then
1967  ntime=size(volgrid6d_in%time)
1968 
1969  if (vol7d_out%time_definition == volgrid6d_in%time_definition) then
1970 
1971  ! i time sono definiti uguali: assegno
1972  vol7d_out%time=volgrid6d_in%time
1973 
1974  else
1975  ! converto reference in validity
1976  allocate (validitytime(ntime,ntimerange),stat=stallo)
1977  if (stallo /=0)then
1978  call l4f_category_log(volgrid6d_in%category,L4F_FATAL,"allocating memory")
1979  call raise_fatal_error()
1980  end if
1981 
1982  do itime=1,ntime
1983  do itimerange=1,ntimerange
1984  if (vol7d_out%time_definition > volgrid6d_in%time_definition) then
1985  validitytime(itime,itimerange) = &
1986  volgrid6d_in%time(itime) + timedelta_new(sec=volgrid6d_in%timerange(itimerange)%p1)
1987  else
1988  validitytime(itime,itimerange) = &
1989  volgrid6d_in%time(itime) - timedelta_new(sec=volgrid6d_in%timerange(itimerange)%p1)
1990  end if
1991  end do
1992  end do
1993 
1994  nntime = count_distinct(reshape(validitytime,(/ntime*ntimerange/)), back=.TRUE.)
1995  vol7d_out%time=pack_distinct(reshape(validitytime,(/ntime*ntimerange/)), nntime,back=.TRUE.)
1996 
1997  end if
1998 end if
1999 
2000 IF (ASSOCIATED(volgrid6d_in%level)) THEN
2001  nlevel = SIZE(volgrid6d_in%level)
2002  vol7d_out%level=volgrid6d_in%level
2003 ENDIF
2004 
2005 IF (ASSOCIATED(volgrid6d_in%var)) THEN
2006  nvar = SIZE(volgrid6d_in%var)
2007 .NOT. IF ( optio_log(noconvert)) THEN
2008  CALL vargrib2varbufr(volgrid6d_in%var, vol7d_out%dativar%r, c_func)
2009  ENDIF
2010 ENDIF
2011 
2012 nana = SIZE(vol7d_out%ana)
2013 
2014 ! allocate once for speed
2015 .NOT.IF (ASSOCIATED(volgrid6d_in%voldati)) THEN
2016  ALLOCATE(voldatiin(volgrid6d_in%griddim%dim%nx, volgrid6d_in%griddim%dim%ny, &
2017  nlevel))
2018 ENDIF
2019 
2020 ALLOCATE(voldatir_out(nana,1,nlevel),stat=stallo)
2021 IF (stallo /= 0) THEN
2022  CALL l4f_category_log(volgrid6d_in%category,L4F_FATAL,"allocating memory")
2023  CALL raise_fatal_error()
2024 ENDIF
2025 
2026 inetwork=1
2027 do itime=1,ntime
2028  do itimerange=1,ntimerange
2029 ! do ilevel=1,nlevel
2030  do ivar=1,nvar
2031 
2032 Ã¨Ã¨ !non chiaro se questa sezione utile o no
2033  !ossia il compute sotto sembra prevedere voldatir_out solo in out
2034 !!$ if (vol7d_out%time_definition == volgrid6d_in%time_definition) then
2035 !!$ voldatir_out=reshape(vol7d_out%voldatir(:,itime,ilevel,itimerange,ivar,inetwork),(/nana,1/))
2036 !!$ else
2037 !!$ voldatir_out=reshape(vol7d_out%voldatir(:,index(vol7d_out%time,validitytime(itime,ilevel)),ilevel,itimerange,ivar,inetwork),(/nana,1/))
2038 !!$ end if
2039 
2040  CALL volgrid_get_vol_3d(volgrid6d_in, itime, itimerange, ivar, &
2041  voldatiin)
2042 
2043  CALL compute(this, voldatiin, voldatir_out, vol7d_out%dativar%r(ivar))
2044 
2045  if (vol7d_out%time_definition == volgrid6d_in%time_definition) then
2046  vol7d_out%voldatir(:,itime,:,itimerange,ivar,inetwork) = &
2047  voldatir_out(:,1,:)
2048  else
2049  vol7d_out%voldatir(:,index(vol7d_out%time,validitytime(itime,itimerange)),:,itimerange,ivar,inetwork)=&
2050  RESHAPE(voldatir_out,(/nana,nlevel/))
2051  end if
2052 
2053 ! 1 indice della dimensione "anagrafica"
2054 ! 2 indice della dimensione "tempo"
2055 ! 3 indice della dimensione "livello verticale"
2056 ! 4 indice della dimensione "intervallo temporale"
2057 ! 5 indice della dimensione "variabile"
2058 ! 6 indice della dimensione "rete"
2059 
2060  end do
2061 ! end do
2062  end do
2063 end do
2064 
2065 deallocate(voldatir_out)
2066 .NOT.IF (ASSOCIATED(volgrid6d_in%voldati)) THEN
2067  DEALLOCATE(voldatiin)
2068 ENDIF
2069 if (allocated(validitytime)) deallocate(validitytime)
2070 
2071 ! Rescale valid data according to variable conversion table
2072 IF (ASSOCIATED(c_func)) THEN
2073  DO ivar = 1, nvar
2074  CALL compute(c_func(ivar), vol7d_out%voldatir(:,:,:,:,ivar,:))
2075  ENDDO
2076  DEALLOCATE(c_func)
2077 ENDIF
2078 
2079 end SUBROUTINE volgrid6d_v7d_transform_compute
2080 
2081 
2082 !> Performs the specified abstract transformation on the data provided.
2083 !! The abstract transformation is specified by \a this parameter; the
2084 !! corresponding specifical transformation (\a grid_transform object)
2085 !! is created and destroyed internally. The output transformed object
2086 !! is created internally and it does not require preliminary
2087 !! initialisation.
2088 SUBROUTINE volgrid6d_v7d_transform(this, volgrid6d_in, vol7d_out, v7d, &
2089  maskgrid, maskbounds, networkname, noconvert, find_index, categoryappend)
2090 TYPE(transform_def),INTENT(in) :: this !< object specifying the abstract transformation
2091 TYPE(volgrid6d),INTENT(inout) :: volgrid6d_in !< object to be transformed, it is not modified, despite the INTENT(inout)
2092 TYPE(vol7d),INTENT(out) :: vol7d_out !< transformed object, it does not requires initialisation
2093 TYPE(vol7d),INTENT(in),OPTIONAL :: v7d !< object containing a list of points over which transformation has to be done (required by some transformation types)
2094 REAL,INTENT(in),OPTIONAL :: maskgrid(:,:) !< 2D field to be used for defining subareas according to its values, it must have the same shape as the field to be interpolated (for transformation type 'maskinter')
2095 REAL,INTENT(in),OPTIONAL :: maskbounds(:) !< array of boundary values for defining subareas from the values of \a maskgrid, the number of subareas is SIZE(maskbounds) - 1, if not provided a default based on extreme values of \a makgrid is used
2096 CHARACTER(len=*),OPTIONAL,INTENT(in) :: networkname !< set the output network name in vol7d_out (default='generic')
2097 LOGICAL,OPTIONAL,INTENT(in) :: noconvert !< do not try to match variable and convert values during transform
2098 PROCEDURE(basic_find_index),POINTER,OPTIONAL :: find_index
2099 CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
2100 
2101 type(grid_transform) :: grid_trans
2102 INTEGER :: ntime, ntimerange, nlevel, nvar, nana, time_definition, nnetwork, stallo
2103 INTEGER :: itime, itimerange, inetwork
2104 TYPE(datetime),ALLOCATABLE :: validitytime(:,:)
2105 INTEGER,ALLOCATABLE :: point_index(:)
2106 TYPE(vol7d) :: v7d_locana
2107 
2108 #ifdef DEBUG
2109 call l4f_category_log(volgrid6d_in%category,L4F_DEBUG,"start volgrid6d_v7d_transform")
2110 #endif
2111 
2112 call vg6d_wind_unrot(volgrid6d_in)
2113 
2114 ntime=0
2115 ntimerange=0
2116 nlevel=0
2117 nvar=0
2118 nnetwork=1
2119 
2120 call get_val(this,time_definition=time_definition)
2121 .not.if ( c_e(time_definition)) then
2122  time_definition=1 ! default to validity time
2123 endif
2125 IF (PRESENT(v7d)) THEN
2126  CALL vol7d_copy(v7d, v7d_locana)
2127 ELSE
2128  CALL init(v7d_locana)
2129 ENDIF
2130 
2131 if (associated(volgrid6d_in%timerange)) ntimerange=size(volgrid6d_in%timerange)
2132 
2133 if (associated(volgrid6d_in%time)) then
2134 
2135  ntime=size(volgrid6d_in%time)
2136 
2137  if (time_definition /= volgrid6d_in%time_definition) then
2138 
2139  ! converto reference in validity
2140  allocate (validitytime(ntime,ntimerange),stat=stallo)
2141  if (stallo /=0)then
2142  call l4f_category_log(volgrid6d_in%category,L4F_FATAL,"allocating memory")
2143  call raise_fatal_error()
2144  end if
2145 
2146  do itime=1,ntime
2147  do itimerange=1,ntimerange
2148  if (time_definition > volgrid6d_in%time_definition) then
2149  validitytime(itime,itimerange) = &
2150  volgrid6d_in%time(itime) + timedelta_new(sec=volgrid6d_in%timerange(itimerange)%p1)
2151  else
2152  validitytime(itime,itimerange) = &
2153  volgrid6d_in%time(itime) - timedelta_new(sec=volgrid6d_in%timerange(itimerange)%p1)
2154  end if
2155  end do
2156  end do
2157 
2158  ntime = count_distinct(reshape(validitytime,(/ntime*ntimerange/)), back=.TRUE.)
2159  deallocate (validitytime)
2160 
2161  end if
2162 end if
2163 
2164 
2165 if (associated(volgrid6d_in%level)) nlevel=size(volgrid6d_in%level)
2166 if (associated(volgrid6d_in%var)) nvar=size(volgrid6d_in%var)
2167 
2168 CALL init(grid_trans, this, volgrid6d_in%griddim, v7d_locana, &
2169  maskgrid=maskgrid, maskbounds=maskbounds, find_index=find_index, &
2170  categoryappend=categoryappend)
2171 CALL init (vol7d_out,time_definition=time_definition)
2172 
2173 IF (c_e(grid_trans)) THEN
2174 
2175  nana=SIZE(v7d_locana%ana)
2176  CALL vol7d_alloc(vol7d_out, nana=nana, ntime=ntime, nlevel=nlevel, &
2177  ntimerange=ntimerange, ndativarr=nvar, nnetwork=nnetwork)
2178  vol7d_out%ana = v7d_locana%ana
2179 
2180  CALL get_val(grid_trans, output_point_index=point_index)
2181  IF (ALLOCATED(point_index)) THEN
2182 ! check that size(point_index) == nana?
2183  CALL vol7d_alloc(vol7d_out, nanavari=1)
2184  CALL init(vol7d_out%anavar%i(1), 'B01192')
2185  ENDIF
2186 
2187  CALL vol7d_alloc_vol(vol7d_out)
2188 
2189  IF (ALLOCATED(point_index)) THEN
2190  DO inetwork = 1, nnetwork
2191  vol7d_out%volanai(:,1,inetwork) = point_index(:)
2192  ENDDO
2193  ENDIF
2194  CALL compute(grid_trans, volgrid6d_in, vol7d_out, networkname, noconvert)
2195 ELSE
2196  CALL l4f_log(L4F_ERROR, 'vg6d_v7d_transform: transformation not valid')
2197  CALL raise_error()
2198 ENDIF
2199 
2200 CALL delete(grid_trans)
2201 
2202 #ifdef HAVE_DBALLE
2203 ! set variables to a conformal state
2204 CALL vol7d_dballe_set_var_du(vol7d_out)
2205 #endif
2206 
2207 CALL delete(v7d_locana)
2208 
2209 END SUBROUTINE volgrid6d_v7d_transform
2210 
2211 
2212 !> Performs the specified abstract transformation on the arrays of
2213 !! data provided. The abstract transformation is specified by \a this
2214 !! parameter; the corresponding specifical transformation (\a
2215 !! grid_transform object) is created and destroyed internally. The
2216 !! output transformed object is created internally and it does not
2217 !! require preliminary initialisation. The transformation performed on
2218 !! each element of the input \a volgrid6d array object is merged into
2219 !! a single \a vol7d output object.
2220 SUBROUTINE volgrid6dv_v7d_transform(this, volgrid6d_in, vol7d_out, v7d, &
2221  maskgrid, maskbounds, networkname, noconvert, find_index, categoryappend)
2222 TYPE(transform_def),INTENT(in) :: this !< object specifying the abstract transformation
2223 TYPE(volgrid6d),INTENT(inout) :: volgrid6d_in(:) !< object to be transformed, it is an array of volgrid6d objects, each of which will be transformed, it is not modified, despite the INTENT(inout)
2224 TYPE(vol7d),INTENT(out) :: vol7d_out !< transformed object, it does not require initialisation
2225 TYPE(vol7d),INTENT(in),OPTIONAL :: v7d !< object containing a list of points over which transformation has to be done (required by some transformation types)
2226 REAL,INTENT(in),OPTIONAL :: maskgrid(:,:) !< 2D field to be used for defining subareas according to its values, it must have the same shape as the field to be interpolated (for transformation type 'maskinter')
2227 REAL,INTENT(in),OPTIONAL :: maskbounds(:) !< array of boundary values for defining subareas from the values of \a maskgrid, the number of subareas is SIZE(maskbounds) - 1, if not provided a default based on extreme values of \a makgrid is used
2228 CHARACTER(len=*),OPTIONAL,INTENT(in) :: networkname !< set the output network name in vol7d_out (default='generic')
2229 LOGICAL,OPTIONAL,INTENT(in) :: noconvert !< do not try to match variable and convert values during transform
2230 PROCEDURE(basic_find_index),POINTER,OPTIONAL :: find_index
2231 CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
2232 
2233 integer :: i
2234 TYPE(vol7d) :: v7dtmp
2235 
2236 
2237 CALL init(v7dtmp)
2238 CALL init(vol7d_out)
2239 
2240 DO i=1,SIZE(volgrid6d_in)
2241  CALL transform(this, volgrid6d_in(i), v7dtmp, v7d=v7d, &
2242  maskgrid=maskgrid, maskbounds=maskbounds, &
2243  networkname=networkname, noconvert=noconvert, find_index=find_index, &
2244  categoryappend=categoryappend)
2245  CALL vol7d_append(vol7d_out, v7dtmp)
2246 ENDDO
2247 
2248 END SUBROUTINE volgrid6dv_v7d_transform
2249 
2250 
2251 ! Internal method for performing sparse point to grid computations
2252 SUBROUTINE v7d_volgrid6d_transform_compute(this, vol7d_in, volgrid6d_out, networkname, gaid_template)
2253 TYPE(grid_transform),INTENT(in) :: this ! object specifying the specific transformation
2254 type(vol7d), INTENT(in) :: vol7d_in ! object to be transformed
2255 type(volgrid6d), INTENT(inout) :: volgrid6d_out ! transformed object
2256 CHARACTER(len=*),OPTIONAL,INTENT(in) :: networkname ! select the network to be processed from the \a vol7d input object, default the first network
2257 TYPE(grid_id),OPTIONAL,INTENT(in) :: gaid_template ! the template (typically grib_api) to be associated with output data, it also helps in improving variable conversion
2258 
2259 integer :: nana, ntime, ntimerange, nlevel, nvar
2260 INTEGER :: ilevel, itime, itimerange, ivar, inetwork
2261 
2262 REAL,POINTER :: voldatiout(:,:,:)
2263 type(vol7d_network) :: network
2264 TYPE(conv_func), pointer :: c_func(:)
2265 !TODO category sarebbe da prendere da vol7d
2266 #ifdef DEBUG
2267 CALL l4f_category_log(volgrid6d_out%category, L4F_DEBUG, &
2268  'start v7d_volgrid6d_transform_compute')
2269 #endif
2270 
2271 ntime=0
2272 ntimerange=0
2273 nlevel=0
2274 nvar=0
2275 
2276 IF (PRESENT(networkname)) THEN
2277  CALL init(network,name=networkname)
2278  inetwork = INDEX(vol7d_in%network,network)
2279  IF (inetwork <= 0) THEN
2280  CALL l4f_category_log(volgrid6d_out%category,L4F_WARN, &
2281  'network '//TRIM(networkname)//' not found, first network will be transformed')
2282  inetwork = 1
2283  ENDIF
2284 ELSE
2285  inetwork = 1
2286 ENDIF
2287 
2288 ! no time_definition conversion implemented, output will be the same as input
2289 if (associated(vol7d_in%time))then
2290  ntime=size(vol7d_in%time)
2291  volgrid6d_out%time=vol7d_in%time
2292 end if
2293 
2294 if (associated(vol7d_in%timerange))then
2295  ntimerange=size(vol7d_in%timerange)
2296  volgrid6d_out%timerange=vol7d_in%timerange
2297 end if
2298 
2299 if (associated(vol7d_in%level))then
2300  nlevel=size(vol7d_in%level)
2301  volgrid6d_out%level=vol7d_in%level
2302 end if
2303 
2304 if (associated(vol7d_in%dativar%r))then
2305  nvar=size(vol7d_in%dativar%r)
2306  CALL varbufr2vargrib(vol7d_in%dativar%r, volgrid6d_out%var, c_func, gaid_template)
2307 end if
2308 
2309 nana=SIZE(vol7d_in%voldatir, 1)
2310 ! allocate once for speed
2311 .NOT.IF (ASSOCIATED(volgrid6d_out%voldati)) THEN
2312  ALLOCATE(voldatiout(volgrid6d_out%griddim%dim%nx, volgrid6d_out%griddim%dim%ny, &
2313  nlevel))
2314 ENDIF
2315 
2316 DO ivar=1,nvar
2317  DO itimerange=1,ntimerange
2318  DO itime=1,ntime
2319 
2320 ! clone the gaid template where I have data
2321  IF (PRESENT(gaid_template)) THEN
2322  DO ilevel = 1, nlevel
2323  IF (ANY(c_e(vol7d_in%voldatir(:,itime,ilevel,itimerange,ivar,inetwork)))) THEN
2324  CALL copy(gaid_template, volgrid6d_out%gaid(ilevel,itime,itimerange,ivar))
2325  ELSE
2326  volgrid6d_out%gaid(ilevel,itime,itimerange,ivar) = grid_id_new()
2327  ENDIF
2328  ENDDO
2329  ENDIF
2330 
2331 ! get data
2332  IF (ASSOCIATED(volgrid6d_out%voldati)) & ! improve!!!!
2333  CALL volgrid_get_vol_3d(volgrid6d_out, itime, itimerange, ivar, &
2334  voldatiout)
2335 ! do the interpolation
2336  CALL compute(this, &
2337  vol7d_in%voldatir(:,itime,:,itimerange,ivar,inetwork), voldatiout, &
2338  vol7d_in%dativar%r(ivar))
2339 ! rescale valid data according to variable conversion table
2340  IF (ASSOCIATED(c_func)) THEN
2341  CALL compute(c_func(ivar), voldatiout(:,:,:))
2342  ENDIF
2343 ! put data
2344  CALL volgrid_set_vol_3d(volgrid6d_out, itime, itimerange, ivar, &
2345  voldatiout)
2346 
2347  ENDDO
2348  ENDDO
2349 ENDDO
2350 
2351 .NOT.IF (ASSOCIATED(volgrid6d_out%voldati)) THEN
2352  DEALLOCATE(voldatiout)
2353 ENDIF
2354 IF (ASSOCIATED(c_func)) THEN
2355  DEALLOCATE(c_func)
2356 ENDIF
2357 
2358 END SUBROUTINE v7d_volgrid6d_transform_compute
2359 
2360 
2361 !> Performs the specified abstract transformation on the data provided.
2362 !! The abstract transformation is specified by \a this parameter; the
2363 !! corresponding specifical transformation (\a grid_transform object)
2364 !! is created and destroyed internally. The output transformed object
2365 !! is created internally and it does not require preliminary
2366 !! initialisation.
2367 SUBROUTINE v7d_volgrid6d_transform(this, griddim, vol7d_in, volgrid6d_out, &
2368  networkname, gaid_template, categoryappend)
2369 TYPE(transform_def),INTENT(in) :: this !< object specifying the abstract transformation
2370 TYPE(griddim_def),INTENT(in),OPTIONAL :: griddim !< griddim specifying the output grid (required by most transformation types)
2371 ! TODO ripristinare intent(in) dopo le opportune modifiche in grid_class.F90
2372 TYPE(vol7d),INTENT(inout) :: vol7d_in !< object to be transformed, it is not modified, despite the INTENT(inout)
2373 TYPE(volgrid6d),INTENT(out) :: volgrid6d_out !< transformed object, it does not require initialisation
2374 CHARACTER(len=*),OPTIONAL,INTENT(in) :: networkname !< select the network to be processed from the \a vol7d input object, default the first network
2375 TYPE(grid_id),OPTIONAL,INTENT(in) :: gaid_template !< the template (typically grib_api) to be associated with output data, it also helps in improving variable conversion
2376 CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
2377 
2378 type(grid_transform) :: grid_trans
2379 integer :: ntime, ntimerange, nlevel, nvar
2380 
2381 
2382 !TODO la category sarebbe da prendere da vol7d
2383 !call l4f_category_log(vol7d_out%category,L4F_DEBUG,"start volgrid6d_transform")
2384 
2385 CALL vol7d_alloc_vol(vol7d_in) ! be safe
2386 ntime=SIZE(vol7d_in%time)
2387 ntimerange=SIZE(vol7d_in%timerange)
2388 nlevel=SIZE(vol7d_in%level)
2389 nvar=0
2390 if (associated(vol7d_in%dativar%r)) nvar=size(vol7d_in%dativar%r)
2391 
2392 IF (nvar <= 0) THEN ! use vol7d category once it will be implemented
2393  CALL l4f_log(L4F_ERROR, &
2394  "trying to transform a vol7d object incomplete or without real variables")
2395  CALL init(volgrid6d_out) ! initialize to empty
2396  CALL raise_error()
2397  RETURN
2398 ENDIF
2399 
2400 CALL init(grid_trans, this, vol7d_in, griddim, categoryappend=categoryappend)
2401 CALL init(volgrid6d_out, griddim, time_definition=vol7d_in%time_definition, &
2402  categoryappend=categoryappend)
2403 
2404 IF (c_e(grid_trans)) THEN
2405 
2406  CALL volgrid6d_alloc(volgrid6d_out, griddim%dim, ntime=ntime, nlevel=nlevel, &
2407  ntimerange=ntimerange, nvar=nvar)
2408 ! can I avoid decode=.TRUE. here, using gaid_template?
2409  CALL volgrid6d_alloc_vol(volgrid6d_out, decode=.TRUE.)
2410 
2411  CALL compute(grid_trans, vol7d_in, volgrid6d_out, networkname, gaid_template)
2412 
2413  CALL vg6d_wind_rot(volgrid6d_out)
2414 ELSE
2415 ! should log with grid_trans%category, but it is private
2416  CALL l4f_log(L4F_ERROR, 'v7d_vg6d_transform: transformation not valid')
2417  CALL raise_error()
2418 ENDIF
2419 
2420 CALL delete(grid_trans)
2421 
2422 END SUBROUTINE v7d_volgrid6d_transform
2423 
2424 
2425 ! Internal method for performing sparse point to sparse point computations
2426 SUBROUTINE v7d_v7d_transform_compute(this, vol7d_in, vol7d_out, lev_out, &
2427  var_coord_vol)
2428 TYPE(grid_transform),INTENT(in) :: this ! oggetto di trasformazione per grigliato
2429 type(vol7d), INTENT(in) :: vol7d_in ! oggetto da trasformare
2430 type(vol7d), INTENT(inout) :: vol7d_out ! oggetto trasformato
2431 TYPE(vol7d_level),INTENT(in),OPTIONAL :: lev_out(:) ! vol7d_level object defining target vertical grid, for vertical interpolations
2432 INTEGER,INTENT(in),OPTIONAL :: var_coord_vol ! index of variable defining vertical coordinate values in input volume
2433 
2434 INTEGER :: itime, itimerange, ilevel, ivar, inetwork, &
2435  levshift, levused, lvar_coord_vol, spos
2436 REAL,ALLOCATABLE :: coord_3d_in(:,:,:)
2437 TYPE(vol7d_level) :: output_levtype
2438 
2439 lvar_coord_vol = optio_i(var_coord_vol)
2440 vol7d_out%time(:) = vol7d_in%time(:)
2441 vol7d_out%timerange(:) = vol7d_in%timerange(:)
2442 IF (PRESENT(lev_out)) THEN
2443  vol7d_out%level(:) = lev_out(:)
2444 ELSE
2445  vol7d_out%level(:) = vol7d_in%level(:)
2446 ENDIF
2447 vol7d_out%network(:) = vol7d_in%network(:)
2448 IF (ASSOCIATED(vol7d_in%dativar%r)) THEN ! work only when real vars are available
2449  vol7d_out%dativar%r(:) = vol7d_in%dativar%r(:)
2450 
2451  CALL get_val(this, levshift=levshift, levused=levused)
2452  spos = imiss
2453  IF (c_e(lvar_coord_vol)) THEN
2454  CALL get_val(this%trans, output_levtype=output_levtype)
2455 .OR. IF (output_levtype%level1 == 103 output_levtype%level1 == 108) THEN
2456  spos = firsttrue(vol7d_in%level(:) == vol7d_level_new(1))
2457  IF (spos == 0) THEN
2458  CALL l4f_log(L4F_ERROR, &
2459  'output level '//t2c(output_levtype%level1)// &
2460  ' requested, but height/press of surface not provided in volume')
2461  ENDIF
2462 .NOT..AND..NOT. IF (c_e(levshift) c_e(levused)) THEN
2463  CALL l4f_log(L4F_ERROR, &
2464  'internal inconsistence, levshift and levused undefined when they should be')
2465  ENDIF
2466  ALLOCATE(coord_3d_in(SIZE(vol7d_in%ana),1,SIZE(vol7d_in%level)))
2467  ENDIF
2468 
2469  ENDIF
2470 
2471  DO inetwork = 1, SIZE(vol7d_in%network)
2472  DO ivar = 1, SIZE(vol7d_in%dativar%r)
2473  DO itimerange = 1, SIZE(vol7d_in%timerange)
2474  DO itime = 1, SIZE(vol7d_in%time)
2475 
2476 ! dirty trick to make voldatir look like a 2d-array of shape (nana,1)
2477  IF (c_e(lvar_coord_vol)) THEN
2478  IF (c_e(spos)) THEN ! compute difference wrt surface coordinate
2479  IF (spos == 0) THEN ! error condition, set all to missing and goodnight
2480  coord_3d_in(:,:,levshift+1:levshift+levused) = rmiss
2481  ELSE
2482  DO ilevel = levshift+1, levshift+levused
2483 .AND. WHERE(c_e(vol7d_in%voldatir(:,itime:itime,ilevel,itimerange,lvar_coord_vol,inetwork)) &
2484  c_e(vol7d_in%voldatir(:,itime:itime,spos,itimerange,lvar_coord_vol,inetwork)))
2485  coord_3d_in(:,:,ilevel) = vol7d_in%voldatir(:,itime:itime,ilevel,itimerange,lvar_coord_vol,inetwork) - &
2486  vol7d_in%voldatir(:,itime:itime,spos,itimerange,lvar_coord_vol,inetwork)
2487  ELSEWHERE
2488  coord_3d_in(:,:,ilevel) = rmiss
2489  END WHERE
2490  ENDDO
2491  ENDIF
2492  CALL compute(this, &
2493  vol7d_in%voldatir(:,itime,:,itimerange,ivar,inetwork), &
2494  vol7d_out%voldatir(:,itime:itime,:,itimerange,ivar,inetwork), &
2495  var=vol7d_in%dativar%r(ivar), &
2496  coord_3d_in=coord_3d_in)
2497  ELSE
2498  CALL compute(this, &
2499  vol7d_in%voldatir(:,itime,:,itimerange,ivar,inetwork), &
2500  vol7d_out%voldatir(:,itime:itime,:,itimerange,ivar,inetwork), &
2501  var=vol7d_in%dativar%r(ivar), &
2502  coord_3d_in=vol7d_in%voldatir(:,itime:itime,:,itimerange, &
2503  lvar_coord_vol,inetwork))
2504  ENDIF
2505  ELSE
2506  CALL compute(this, &
2507  vol7d_in%voldatir(:,itime,:,itimerange,ivar,inetwork), &
2508  vol7d_out%voldatir(:,itime:itime,:,itimerange,ivar,inetwork), &
2509  var=vol7d_in%dativar%r(ivar))
2510  ENDIF
2511  ENDDO
2512  ENDDO
2513  ENDDO
2514  ENDDO
2515 
2516 ENDIF
2517 
2518 END SUBROUTINE v7d_v7d_transform_compute
2519 
2520 
2521 !> Performs the specified abstract transformation on the data provided.
2522 !! The abstract transformation is specified by \a this parameter; the
2523 !! corresponding specifical transformation (\a grid_transform object)
2524 !! is created and destroyed internally. The output transformed object
2525 !! is created internally and it does not require preliminary
2526 !! initialisation. The success of the transformation can be checked
2527 !! with the \a c_e method: c_e(vol7d_out).
2528 SUBROUTINE v7d_v7d_transform(this, vol7d_in, vol7d_out, v7d, maskbounds, &
2529  lev_out, vol7d_coord_in, categoryappend)
2530 TYPE(transform_def),INTENT(in) :: this !< object specifying the abstract transformation
2531 TYPE(vol7d),INTENT(inout) :: vol7d_in !< object to be transformed, it is not modified, despite the INTENT(inout)
2532 TYPE(vol7d),INTENT(out) :: vol7d_out !< transformed object, it does not require initialisation
2533 TYPE(vol7d),INTENT(in),OPTIONAL :: v7d !< object containing a list of points over which transformation has to be done (required by some transformation types)
2534 REAL,INTENT(in),OPTIONAL :: maskbounds(:) !< array of boundary values for defining a subset of valid points where the values of \a maskgrid are within the first and last value of \a maskbounds (for transformation type 'metamorphosis:maskfill')
2535 TYPE(vol7d_level),INTENT(in),OPTIONAL,TARGET :: lev_out(:) !< vol7d_level object defining target vertical grid, for vertical interpolations
2536 TYPE(vol7d),INTENT(in),OPTIONAL :: vol7d_coord_in !< object providing time constant input vertical coordinate for some kind of vertical interpolations
2537 CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
2538 
2539 INTEGER :: nvar, inetwork
2540 TYPE(grid_transform) :: grid_trans
2541 TYPE(vol7d_level),POINTER :: llev_out(:)
2542 TYPE(vol7d_level) :: input_levtype, output_levtype
2543 TYPE(vol7d_var) :: vcoord_var
2544 REAL,ALLOCATABLE :: coord_3d_in(:,:,:)
2545 INTEGER :: var_coord_in, var_coord_vol, i, k, ulstart, ulend, spos
2546 INTEGER,ALLOCATABLE :: point_index(:)
2547 TYPE(vol7d) :: v7d_locana, vol7d_tmpana
2548 CHARACTER(len=80) :: trans_type
2549 LOGICAL,ALLOCATABLE :: mask_in(:), point_mask(:)
2550 
2551 CALL vol7d_alloc_vol(vol7d_in) ! be safe
2552 nvar=0
2553 IF (ASSOCIATED(vol7d_in%dativar%r)) nvar=SIZE(vol7d_in%dativar%r)
2554 
2555 CALL init(v7d_locana)
2556 IF (PRESENT(v7d)) v7d_locana = v7d
2557 CALL init(vol7d_out, time_definition=vol7d_in%time_definition)
2558 
2559 CALL get_val(this, trans_type=trans_type)
2560 
2561 var_coord_vol = imiss
2562 IF (trans_type == 'vertint') THEN
2563 
2564  IF (PRESENT(lev_out)) THEN
2565 
2566 ! if vol7d_coord_in provided and allocated, check that it fits
2567  var_coord_in = -1
2568  IF (PRESENT(vol7d_coord_in)) THEN
2569 .AND. IF (ASSOCIATED(vol7d_coord_in%voldatir) &
2570  ASSOCIATED(vol7d_coord_in%dativar%r)) THEN
2571 
2572 ! strictly 1 time, 1 timerange and 1 network
2573 .OR. IF (SIZE(vol7d_coord_in%voldatir,2) /= 1 &
2574 .OR. SIZE(vol7d_coord_in%voldatir,4) /= 1 &
2575  SIZE(vol7d_coord_in%voldatir,6) /= 1) THEN
2576  CALL l4f_log(L4F_ERROR, &
2577  'volume providing constant input vertical coordinate must have &
2578  &only 1 time, 1 timerange and 1 network')
2579  CALL raise_error()
2580  RETURN
2581  ENDIF
2582 
2583 ! search for variable providing vertical coordinate
2584  CALL get_val(this, output_levtype=output_levtype)
2585  vcoord_var = vol7d_var_new(vol7d_level_to_var(output_levtype))
2586 .NOT. IF (c_e(vcoord_var)) THEN
2587  CALL l4f_log(L4F_ERROR, &
2588  'requested output level type '//t2c(output_levtype%level1)// &
2589  ' does not correspond to any known physical variable for &
2590  &providing vertical coordinate')
2591  CALL raise_error()
2592  RETURN
2593  ENDIF
2594 
2595  var_coord_in = INDEX(vol7d_coord_in%dativar%r, vcoord_var)
2596 
2597  IF (var_coord_in <= 0) THEN
2598  CALL l4f_log(L4F_ERROR, &
2599  'volume providing constant input vertical coordinate contains no &
2600  &real variables matching output level type '//t2c(output_levtype%level1))
2601  CALL raise_error()
2602  RETURN
2603  ENDIF
2604  CALL l4f_log(L4F_INFO, &
2605  'Coordinate for vertint found in coord volume at position '// &
2606  t2c(var_coord_in))
2607 
2608 ! check vertical coordinate system
2609  CALL get_val(this, input_levtype=input_levtype)
2610  mask_in = & ! implicit allocation
2611 .AND. (vol7d_coord_in%level(:)%level1 == input_levtype%level1) &
2612  (vol7d_coord_in%level(:)%level2 == input_levtype%level2)
2613  ulstart = firsttrue(mask_in)
2614  ulend = lasttrue(mask_in)
2615 .OR. IF (ulstart == 0 ulend == 0) THEN
2616  CALL l4f_log(L4F_ERROR, &
2617  'coordinate file does not contain levels of type '// &
2618  t2c(input_levtype%level1)//'/'//t2c(input_levtype%level2)// &
2619  ' specified for input data')
2620  CALL raise_error()
2621  RETURN
2622  ENDIF
2623 
2624  coord_3d_in = vol7d_coord_in%voldatir(:,1:1,ulstart:ulend,1,var_coord_in,1) ! dirty 1:1, implicit allocation
2625 ! special case
2626  IF (output_levtype%level1 == 103 &
2627 .OR. output_levtype%level1 == 108) THEN ! surface coordinate needed
2628  spos = firsttrue(vol7d_coord_in%level(:) == vol7d_level_new(1))
2629  IF (spos == 0) THEN
2630  CALL l4f_log(L4F_ERROR, &
2631  'output level '//t2c(output_levtype%level1)// &
2632  ' requested, but height/press of surface not provided in coordinate file')
2633  CALL raise_error()
2634  RETURN
2635  ENDIF
2636  DO k = 1, SIZE(coord_3d_in,3)
2637 .AND. WHERE(c_e(coord_3d_in(:,:,k)) &
2638  c_e(vol7d_coord_in%voldatir(:,1:1,spos,1,var_coord_in,1)))
2639  coord_3d_in(:,:,k) = coord_3d_in(:,:,k) - &
2640  vol7d_coord_in%voldatir(:,1:1,spos,1,var_coord_in,1)
2641  ELSEWHERE
2642  coord_3d_in(:,:,k) = rmiss
2643  END WHERE
2644  ENDDO
2645  ENDIF
2646 
2647  ENDIF
2648  ENDIF
2649 
2650  IF (var_coord_in <= 0) THEN ! search for coordinate within volume
2651 ! search for variable providing vertical coordinate
2652  CALL get_val(this, output_levtype=output_levtype)
2653  vcoord_var = vol7d_var_new(vol7d_level_to_var(output_levtype))
2654  IF (c_e(vcoord_var)) THEN
2655  DO i = 1, SIZE(vol7d_in%dativar%r)
2656  IF (vol7d_in%dativar%r(i) == vcoord_var) THEN
2657  var_coord_vol = i
2658  EXIT
2659  ENDIF
2660  ENDDO
2661 
2662  IF (c_e(var_coord_vol)) THEN
2663  CALL l4f_log(L4F_INFO, &
2664  'Coordinate for vertint found in input volume at position '// &
2665  t2c(var_coord_vol))
2666  ENDIF
2667 
2668  ENDIF
2669  ENDIF
2670 
2671  IF (var_coord_in > 0) THEN
2672  CALL init(grid_trans, this, lev_in=vol7d_in%level, lev_out=lev_out, &
2673  coord_3d_in=coord_3d_in, categoryappend=categoryappend)
2674  ELSE
2675  CALL init(grid_trans, this, lev_in=vol7d_in%level, lev_out=lev_out, &
2676  categoryappend=categoryappend)
2677  ENDIF
2678 
2679  CALL get_val(grid_trans, output_level_auto=llev_out) ! get levels if auto-generated
2680 .NOT. IF (associated(llev_out)) llev_out => lev_out
2681 
2682 .AND. IF (c_e(grid_trans)) THEN! nvar > 0) THEN
2683 
2684  CALL vol7d_alloc(vol7d_out, nana=SIZE(vol7d_in%ana), &
2685  ntime=SIZE(vol7d_in%time), ntimerange=SIZE(vol7d_in%timerange), &
2686  nlevel=SIZE(llev_out), nnetwork=SIZE(vol7d_in%network), ndativarr=nvar)
2687  vol7d_out%ana(:) = vol7d_in%ana(:)
2688 
2689  CALL vol7d_alloc_vol(vol7d_out)
2690 
2691 ! no need to check c_e(var_coord_vol) here since the presence of
2692 ! this%coord_3d_in (external) has precedence over coord_3d_in internal
2693 ! in grid_transform_compute
2694  CALL compute(grid_trans, vol7d_in, vol7d_out, llev_out, &
2695  var_coord_vol=var_coord_vol)
2696  ELSE
2697  CALL l4f_log(L4F_ERROR, 'v7d_v7d_transform: transformation not valid')
2698  CALL raise_error()
2699  ENDIF
2700  ELSE
2701  CALL l4f_log(L4F_ERROR, &
2702  'v7d_v7d_transform: vertint requested but lev_out not provided')
2703  CALL raise_error()
2704  ENDIF
2705 
2706 ELSE
2707 
2708  CALL init(grid_trans, this, vol7d_in, v7d_locana, maskbounds=maskbounds, &
2709  categoryappend=categoryappend)
2710 ! if this init is successful, I am sure that v7d_locana%ana is associated
2711 
2712 .AND. IF (c_e(grid_trans)) THEN! nvar > 0) THEN
2713 
2714  CALL vol7d_alloc(vol7d_out, nana=SIZE(v7d_locana%ana), &
2715  ntime=SIZE(vol7d_in%time), ntimerange=SIZE(vol7d_in%timerange), &
2716  nlevel=SIZE(vol7d_in%level), nnetwork=SIZE(vol7d_in%network), ndativarr=nvar)
2717  vol7d_out%ana = v7d_locana%ana
2718 
2719  CALL get_val(grid_trans, point_mask=point_mask, output_point_index=point_index)
2720 
2721  IF (ALLOCATED(point_index)) THEN
2722  CALL vol7d_alloc(vol7d_out, nanavari=1)
2723  CALL init(vol7d_out%anavar%i(1), 'B01192')
2724  ENDIF
2725 
2726  CALL vol7d_alloc_vol(vol7d_out)
2727 
2728  IF (ALLOCATED(point_index)) THEN
2729  DO inetwork = 1, SIZE(vol7d_in%network)
2730  vol7d_out%volanai(:,1,inetwork) = point_index(:)
2731  ENDDO
2732  ENDIF
2733  CALL compute(grid_trans, vol7d_in, vol7d_out)
2734 
2735  IF (ALLOCATED(point_mask)) THEN ! keep full ana
2736  IF (SIZE(point_mask) /= SIZE(vol7d_in%ana)) THEN
2737  CALL l4f_log(L4F_WARN, &
2738  'v7d_v7d_transform: inconsistency in point size: '//t2c(SIZE(point_mask)) &
2739  //':'//t2c(SIZE(vol7d_in%ana)))
2740  ELSE
2741 #ifdef DEBUG
2742  CALL l4f_log(L4F_DEBUG, 'v7d_v7d_transform: merging ana from in to out')
2743 #endif
2744  CALL vol7d_copy(vol7d_in, vol7d_tmpana, &
2745  lana=point_mask, lnetwork=(/.TRUE./), &
2746  ltime=(/.FALSE./), ltimerange=(/.FALSE./), llevel=(/.FALSE./))
2747  CALL vol7d_append(vol7d_out, vol7d_tmpana)
2748  ENDIF
2749  ENDIF
2750 
2751  ELSE
2752  CALL l4f_log(L4F_ERROR, 'v7d_v7d_transform: transformation not valid')
2753  CALL raise_error()
2754  ENDIF
2755 
2756 ENDIF
2757 
2758 CALL delete (grid_trans)
2759 .NOT.IF ( PRESENT(v7d)) CALL delete(v7d_locana)
2760 
2761 END SUBROUTINE v7d_v7d_transform
2762 
2763 
2764 !> Unrotate the wind components.
2765 !! It converts u and v components of vector quantities relative to the
2766 !! defined grid in the direction of increasing x and y coordinates to
2767 !! u and v components relative to easterly and notherly direction. The
2768 !! original fields are overwritten.
2769 !! \todo Check and correct wind component flag (to be moved in
2770 !! griddim_def?)
2771 subroutine vg6d_wind_unrot(this)
2772 type(volgrid6d) :: this !< object containing wind to be unrotated
2773 
2774 integer :: component_flag
2775 
2776 call get_val(this%griddim,component_flag=component_flag)
2777 
2778 if (component_flag == 1) then
2779  call l4f_category_log(this%category,L4F_INFO, &
2780  "unrotating vector components")
2781  call vg6d_wind__un_rot(this,.false.)
2782  call set_val(this%griddim,component_flag=0)
2783 else
2784  call l4f_category_log(this%category,L4F_INFO, &
2785  "no need to unrotate vector components")
2786 end if
2787 
2788 end subroutine vg6d_wind_unrot
2789 
2790 
2791 !> Rotate the wind components.
2792 !! It converts u and v components of vector quantities
2793 !! relative to easterly and notherly direction to
2794 !! defined grid in the direction of increasing x and y coordinates.
2795 !! The original fields are overwritten.
2796 subroutine vg6d_wind_rot(this)
2797 type(volgrid6d) :: this !< object containing wind to be rotated
2798 
2799 integer :: component_flag
2800 
2801 call get_val(this%griddim,component_flag=component_flag)
2802 
2803 if (component_flag == 0) then
2804  call l4f_category_log(this%category,L4F_INFO, &
2805  "rotating vector components")
2806  call vg6d_wind__un_rot(this,.true.)
2807  call set_val(this%griddim,component_flag=1)
2808 else
2809  call l4f_category_log(this%category,L4F_INFO, &
2810  "no need to rotate vector components")
2811 end if
2812 
2813 end subroutine vg6d_wind_rot
2814 
2815 
2816 ! Generic UnRotate the wind components.
2817 SUBROUTINE vg6d_wind__un_rot(this,rot)
2818 TYPE(volgrid6d) :: this ! object containing wind to be (un)rotated
2819 LOGICAL :: rot ! if .true. rotate else unrotate
2820 
2821 INTEGER :: i, j, k, l, a11, a12, a21, a22, stallo
2822 double precision,pointer :: rot_mat(:,:,:)
2823 real,allocatable :: tmp_arr(:,:)
2824 REAL,POINTER :: voldatiu(:,:), voldativ(:,:)
2825 INTEGER,POINTER :: iu(:), iv(:)
2826 
2827 .NOT.IF (ASSOCIATED(this%var)) THEN
2828  CALL l4f_category_log(this%category, L4F_ERROR, &
2829  "trying to unrotate an incomplete volgrid6d object")
2830  CALL raise_fatal_error()
2831 ! RETURN
2832 ENDIF
2833 
2834 CALL volgrid6d_var_hor_comp_index(this%var, iu, iv)
2835 .NOT.IF (ASSOCIATED(iu)) THEN
2836  CALL l4f_category_log(this%category,L4F_ERROR, &
2837  "unrotation impossible")
2838  CALL raise_fatal_error()
2839 ! RETURN
2840 ENDIF
2841 
2842 ! Temporary workspace
2843 ALLOCATE(tmp_arr(this%griddim%dim%nx, this%griddim%dim%ny),stat=stallo)
2844 IF (stallo /= 0) THEN
2845  CALL l4f_category_log(this%category, L4F_FATAL, "allocating memory")
2846  CALL raise_fatal_error()
2847 ENDIF
2848 ! allocate once for speed
2849 .NOT.IF (ASSOCIATED(this%voldati)) THEN
2850  ALLOCATE(voldatiu(this%griddim%dim%nx, this%griddim%dim%ny), &
2851  voldativ(this%griddim%dim%nx, this%griddim%dim%ny))
2852 ENDIF
2853 
2854 CALL griddim_unproj(this%griddim)
2855 CALL wind_unrot(this%griddim, rot_mat)
2856 
2857 a11=1
2858 if (rot)then
2859  a12=2
2860  a21=3
2861 else
2862  a12=3
2863  a21=2
2864 end if
2865 a22=4
2866 
2867 DO l = 1, SIZE(iu)
2868  DO k = 1, SIZE(this%timerange)
2869  DO j = 1, SIZE(this%time)
2870  DO i = 1, SIZE(this%level)
2871 ! get data
2872  CALL volgrid_get_vol_2d(this, i, j, k, iu(l), voldatiu)
2873  CALL volgrid_get_vol_2d(this, i, j, k, iv(l), voldativ)
2874 ! convert units forward
2875 ! CALL compute(conv_fwd(iu(l)), voldatiu)
2876 ! CALL compute(conv_fwd(iv(l)), voldativ)
2877 
2878 ! multiply wind components by rotation matrix
2879 .AND. WHERE(voldatiu /= rmiss voldativ /= rmiss)
2880  tmp_arr(:,:) = real(voldatiu(:,:)*rot_mat(:,:,a11) + &
2881  voldativ(:,:)*rot_mat(:,:,a12))
2882  voldativ(:,:) = real(voldatiu(:,:)*rot_mat(:,:,a21) + &
2883  voldativ(:,:)*rot_mat(:,:,a22))
2884  voldatiu(:,:) = tmp_arr(:,:)
2885  END WHERE
2886 ! convert units backward
2887 ! CALL uncompute(conv_fwd(iu(l)), voldatiu)
2888 ! CALL uncompute(conv_fwd(iv(l)), voldativ)
2889 ! put data
2890  CALL volgrid_set_vol_2d(this, i, j, k, iu(l), voldatiu)
2891  CALL volgrid_set_vol_2d(this, i, j, k, iv(l), voldativ)
2892  ENDDO
2893  ENDDO
2894  ENDDO
2895 ENDDO
2896 
2897 .NOT.IF (ASSOCIATED(this%voldati)) THEN
2898  DEALLOCATE(voldatiu, voldativ)
2899 ENDIF
2900 DEALLOCATE(rot_mat, tmp_arr, iu, iv)
2901 
2902 END SUBROUTINE vg6d_wind__un_rot
2903 
2904 
2905 !!$ try to understand the problem:
2906 !!$
2907 !!$ case:
2908 !!$
2909 !!$ 1) we have only one volume: we have to provide the direction of shift
2910 !!$ compute H and traslate on it
2911 !!$ 2) we have two volumes:
2912 !!$ 1) volume U and volume V: compute H and traslate on it
2913 !!$ 2) volume U/V and volume H : translate U/V on H
2914 !!$ 3) we have tree volumes: translate U and V on H
2915 !!$
2916 !!$ strange cases:
2917 !!$ 1) do not have U in volume U
2918 !!$ 2) do not have V in volume V
2919 !!$ 3) we have others variables more than U and V in volumes U e V
2920 !!$
2921 !!$
2922 !!$ so the steps are:
2923 !!$ 1) find the volumes
2924 !!$ 2) define or compute H grid
2925 !!$ 3) trasform the volumes in H
2926 
2927 !!$ N.B.
2928 !!$ case 1) for only one vg6d (U or V) is not managed, but
2929 !!$ the not pubblic subroutines will work but you have to know what you want to do
2930 
2931 
2932 !> Convert grids type C to type A.
2933 !! Use this to interpolate data from grid type C to grid type A
2934 !! Grids type are defined by Arakawa.
2935 !!
2936 !! We need to find these types of area in a vg6d array
2937 !! T area of points with temterature etc.
2938 !! U area of points with u components of winds
2939 !! V area of points with v components of winds
2940 !!
2941 !! this method works if it finds
2942 !! two volumes:
2943 !! 1) volume U and volume V: compute H and traslate on it
2944 !! 2) volume U/V and volume H : translate U/V on H
2945 !! three volumes: translate U and V on H
2946 !!
2947 !! try to work well on more datasets at once
2948 subroutine vg6d_c2a (this)
2949 
2950 TYPE(volgrid6d),INTENT(inout) :: this(:) !< vettor of volumes volgrid6d to elaborate
2951 
2952 integer :: ngrid,igrid,jgrid,ugrid,vgrid,tgrid
2953 doubleprecision :: xmin, xmax, ymin, ymax
2954 doubleprecision :: xmin_t, xmax_t, ymin_t, ymax_t
2955 doubleprecision :: step_lon_t,step_lat_t
2956 character(len=80) :: type_t,type
2957 TYPE(griddim_def):: griddim_t
2958 
2959 ngrid=size(this)
2960 
2961 do igrid=1,ngrid
2962 
2963  call init(griddim_t)
2964 
2965  call get_val(this(igrid)%griddim,xmin=xmin_t, xmax=xmax_t, ymin=ymin_t, ymax=ymax_t,proj_type=type_t)
2966  step_lon_t=(xmax_t-xmin_t)/dble(this(igrid)%griddim%dim%nx-1)
2967  step_lat_t=(ymax_t-ymin_t)/dble(this(igrid)%griddim%dim%ny-1)
2968 
2969  do jgrid=1,ngrid
2970 
2971  ugrid=imiss
2972  vgrid=imiss
2973  tgrid=imiss
2974 
2975 #ifdef DEBUG
2976  call l4f_category_log(this(igrid)%category,L4F_DEBUG,"c grid: search u/v/t points:"//to_char(igrid)//to_char(jgrid))
2977  call l4f_category_log(this(igrid)%category,L4F_DEBUG,"c grid: test delta: "//to_char(step_lon_t)//to_char(step_lat_t))
2978 #endif
2979 
2980  if (this(igrid)%griddim == this(jgrid)%griddim ) cycle
2981 
2982 .and. if (this(igrid)%griddim%dim%nx == this(jgrid)%griddim%dim%nx &
2983  this(igrid)%griddim%dim%ny == this(jgrid)%griddim%dim%ny ) then
2984 
2985  call get_val(this(jgrid)%griddim,xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,proj_type=type)
2986 
2987  if (type_t /= type )cycle
2988 
2989 #ifdef DEBUG
2990  call l4f_category_log(this(igrid)%category,L4F_DEBUG,"c grid: test u "//&
2991  to_char(xmin)//to_char(xmax)//to_char(ymin)//to_char(ymax))
2992 
2993  call l4f_category_log(this(igrid)%category,L4F_DEBUG,"diff coordinate lon"//&
2994  to_char(abs(xmin - (xmin_t+step_lon_t/2.d0)))//&
2995  to_char(abs(xmax - (xmax_t+step_lon_t/2.d0))))
2996  call l4f_category_log(this(igrid)%category,L4F_DEBUG,"diff coordinate lat"//&
2997  to_char(abs(ymin - (ymin_t+step_lat_t/2.d0)))//&
2998  to_char(abs(ymax - (ymax_t+step_lat_t/2.d0))))
2999 #endif
3000 
3001 .and. if ( abs(xmin - (xmin_t+step_lon_t/2.d0)) < 1.d-3 abs(xmax - (xmax_t+step_lon_t/2.d0)) < 1.d-3 ) then
3002 .and. if ( abs(ymin - ymin_t) < 1.d-3 abs(ymax - ymax_t) < 1.d-3 ) then
3003 
3004 #ifdef DEBUG
3005  call l4f_category_log(this(igrid)%category,L4F_DEBUG,"c grid: found u")
3006 #endif
3007  ugrid=jgrid
3008  tgrid=igrid
3009 
3010  end if
3011  end if
3012 
3013 #ifdef DEBUG
3014  call l4f_category_log(this(igrid)%category,L4F_DEBUG,"c grid: test v "//&
3015  to_char(xmin)//to_char(xmax)//to_char(ymin)//to_char(ymax))
3016 #endif
3017 
3018 .and. if ( abs(ymin - (ymin_t+step_lat_t/2.d0)) < 1.d-3 abs(ymax - (ymax_t+step_lat_t/2.d0)) < 1.d-3 ) then
3019 .and. if ( abs(xmin - xmin_t) < 1.d-3 abs(xmax - xmax_t) < 1.d-3 ) then
3020 
3021 #ifdef DEBUG
3022  call l4f_category_log(this(igrid)%category,L4F_DEBUG,"c grid: found v")
3023 #endif
3024  vgrid=jgrid
3025  tgrid=igrid
3026 
3027  end if
3028  end if
3029 
3030 
3031  ! test if we have U and V
3032 
3033 #ifdef DEBUG
3034  call l4f_category_log(this(igrid)%category,L4F_DEBUG,"c grid: test u and v"//&
3035  to_char(xmin_t)//to_char(xmax_t)//to_char(ymin_t)//to_char(ymax_t)//&
3036  to_char(xmin)//to_char(xmax)//to_char(ymin)//to_char(ymax))
3037 
3038  call l4f_category_log(this(igrid)%category,L4F_DEBUG,"uv diff coordinate lon"//&
3039  to_char(abs(xmin_t - xmin)-step_lon_t/2.d0)//&
3040  to_char(abs(xmax_t - xmax)-step_lon_t/2.d0))
3041  call l4f_category_log(this(igrid)%category,L4F_DEBUG,"uv diff coordinate lat"//&
3042  to_char(abs(ymin_t - ymin) -step_lat_t/2.d0)//&
3043  to_char(abs(ymax_t - ymax)-step_lat_t/2.d0))
3044 #endif
3045 
3046 .and. if ( abs(ymin - (ymin_t+step_lat_t/2.d0)) < 2.d-3 abs(ymax - (ymax_t+step_lat_t/2.d0)) < 2.d-3 ) then
3047 .and. if ( abs(xmin_t - (xmin+step_lon_t/2.d0)) < 2.d-3 abs(xmax_t - (xmax+step_lon_t/2.d0)) < 2.d-3 ) then
3048 
3049 #ifdef DEBUG
3050  call l4f_category_log(this(igrid)%category,L4F_DEBUG,"c grid: found u and v case up and right")
3051 #endif
3052 
3053  vgrid=jgrid
3054  ugrid=igrid
3055 
3056  call init(griddim_t,xmin=xmin, xmax=xmax, ymin=ymin_t, ymax=ymax_t)
3057 
3058  end if
3059  end if
3060  end if
3061 
3062  ! abbiamo almeno due volumi: riportiamo U e/o V su T
3063  if (c_e(ugrid)) then
3064 #ifdef DEBUG
3065  call l4f_category_log(this(igrid)%category,L4F_DEBUG,"c grid u points: interpolate u on t "//to_char(tgrid)//to_char(ugrid))
3066 #endif
3067  if (c_e(tgrid))then
3068  call vg6d_c2a_grid(this(ugrid),this(tgrid)%griddim,cgrid=1)
3069  else
3070  call vg6d_c2a_grid(this(ugrid),griddim_t,cgrid=1)
3071  end if
3072  call vg6d_c2a_mat(this(ugrid),cgrid=1)
3073  end if
3074 
3075  if (c_e(vgrid)) then
3076 #ifdef DEBUG
3077  call l4f_category_log(this(igrid)%category,L4F_DEBUG,"c grid v points: interpolate v on t "//to_char(tgrid)//to_char(vgrid))
3078 #endif
3079  if (c_e(tgrid))then
3080  call vg6d_c2a_grid(this(vgrid),this(tgrid)%griddim,cgrid=2)
3081  else
3082  call vg6d_c2a_grid(this(vgrid),griddim_t,cgrid=2)
3083  end if
3084  call vg6d_c2a_mat(this(vgrid),cgrid=2)
3085  end if
3086 
3087  end do
3088 
3089  call delete(griddim_t)
3090 
3091 end do
3092 
3093 
3094 end subroutine vg6d_c2a
3095 
3096 
3097 ! Convert C grid to A grid
3098 subroutine vg6d_c2a_grid(this,griddim_t,cgrid)
3099 
3100 type(volgrid6d),intent(inout) :: this ! object containing fields to be translated (U or V points)
3101 type(griddim_def),intent(in),optional :: griddim_t ! object containing grid of T points
3102 integer,intent(in) :: cgrid ! in C grid (Arakawa) we have 0=T,1=U,2=V points
3103 
3104 doubleprecision :: xmin, xmax, ymin, ymax
3105 doubleprecision :: step_lon,step_lat
3106 
3107 
3108 if (present(griddim_t)) then
3109 
3110  call get_val(griddim_t,xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
3111  call set_val(this%griddim,xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
3112 ! improve grid definition precision
3113  CALL griddim_setsteps(this%griddim)
3114 
3115 else
3116 
3117  select case (cgrid)
3118 
3119  case(0)
3120 
3121 #ifdef DEBUG
3122  call l4f_category_log(this%category,L4F_DEBUG,"c grid: t points, nothing to do")
3123 #endif
3124  return
3125 
3126  case (1)
3127 
3128 #ifdef DEBUG
3129  call l4f_category_log(this%category,L4F_DEBUG,"c grid: u points, we need interpolation")
3130 #endif
3131 
3132  call get_val(this%griddim, xmin=xmin, xmax=xmax)
3133  step_lon=(xmax-xmin)/dble(this%griddim%dim%nx-1)
3134  xmin=xmin-step_lon/2.d0
3135  xmax=xmax-step_lon/2.d0
3136  call set_val(this%griddim, xmin=xmin, xmax=xmax)
3137 ! improve grid definition precision
3138  CALL griddim_setsteps(this%griddim)
3139 
3140  case (2)
3141 
3142 #ifdef DEBUG
3143  call l4f_category_log(this%category,L4F_DEBUG,"c grid: v points, we need interpolation")
3144 #endif
3145 
3146  call get_val(this%griddim, ymin=ymin, ymax=ymax)
3147  step_lat=(ymax-ymin)/dble(this%griddim%dim%ny-1)
3148  ymin=ymin-step_lat/2.d0
3149  ymax=ymax-step_lat/2.d0
3150  call set_val(this%griddim, ymin=ymin, ymax=ymax)
3151 ! improve grid definition precision
3152  CALL griddim_setsteps(this%griddim)
3153 
3154  case default
3155 
3156  call l4f_category_log(this%category,L4F_FATAL,"c grid type not known")
3157  call raise_fatal_error ()
3158 
3159  end select
3160 
3161 end if
3162 
3163 
3164 call griddim_unproj(this%griddim)
3165 
3166 
3167 end subroutine vg6d_c2a_grid
3168 
3169 ! Convert C grid to A grid
3170 subroutine vg6d_c2a_mat(this,cgrid)
3171 
3172 type(volgrid6d),intent(inout) :: this ! object containing fields to be translated
3173 integer,intent(in) :: cgrid ! in C grid (Arakawa) we have 0=T,1=U,2=V points
3174 
3175 INTEGER :: i, j, k, iv, stallo
3176 REAL,ALLOCATABLE :: tmp_arr(:,:)
3177 REAL,POINTER :: voldatiuv(:,:)
3178 
3179 
3180 IF (cgrid == 0) RETURN ! nothing to do
3181 .AND.IF (cgrid /= 1 cgrid /= 2) THEN ! wrong cgrid
3182  CALL l4f_category_log(this%category,L4F_FATAL,"c grid type "// &
3183  TRIM(to_char(cgrid))//" not known")
3184  CALL raise_error()
3185  RETURN
3186 ENDIF
3187 
3188 ! Temporary workspace
3189 ALLOCATE(tmp_arr(this%griddim%dim%nx, this%griddim%dim%ny),stat=stallo)
3190 if (stallo /=0)then
3191  call l4f_log(L4F_FATAL,"allocating memory")
3192  call raise_fatal_error()
3193 end if
3194 
3195 ! allocate once for speed
3196 .NOT.IF (ASSOCIATED(this%voldati)) THEN
3197  ALLOCATE(voldatiuv(this%griddim%dim%nx, this%griddim%dim%ny), stat=stallo)
3198  IF (stallo /= 0) THEN
3199  CALL l4f_log(L4F_FATAL,"allocating memory")
3200  CALL raise_fatal_error()
3201  ENDIF
3202 ENDIF
3203 
3204 IF (cgrid == 1) THEN ! U points to H points
3205  DO iv = 1, SIZE(this%var)
3206  DO k = 1, SIZE(this%timerange)
3207  DO j = 1, SIZE(this%time)
3208  DO i = 1, SIZE(this%level)
3209  tmp_arr(:,:) = rmiss
3210  CALL volgrid_get_vol_2d(this, i, j, k, iv, voldatiuv)
3211 
3212 ! West boundary extrapolation
3213 .AND. WHERE(voldatiuv(1,:) /= rmiss voldatiuv(2,:) /= rmiss)
3214  tmp_arr(1,:) = voldatiuv(1,:) - (voldatiuv(2,:) - voldatiuv(1,:)) / 2.
3215  END WHERE
3216 
3217 ! Rest of the field
3218 .AND. WHERE(voldatiuv(1:this%griddim%dim%nx-1,:) /= rmiss &
3219  voldatiuv(2:this%griddim%dim%nx,:) /= rmiss)
3220  tmp_arr(2:this%griddim%dim%nx,:) = &
3221  (voldatiuv(1:this%griddim%dim%nx-1,:) + &
3222  voldatiuv(2:this%griddim%dim%nx,:)) / 2.
3223  END WHERE
3224 
3225  voldatiuv(:,:) = tmp_arr
3226  CALL volgrid_set_vol_2d(this, i, j, k, iv, voldatiuv)
3227  ENDDO
3228  ENDDO
3229  ENDDO
3230  ENDDO
3231 
3232 ELSE IF (cgrid == 2) THEN ! V points to H points
3233  DO iv = 1, SIZE(this%var)
3234  DO k = 1, SIZE(this%timerange)
3235  DO j = 1, SIZE(this%time)
3236  DO i = 1, SIZE(this%level)
3237  tmp_arr(:,:) = rmiss
3238  CALL volgrid_get_vol_2d(this, i, j, k, iv, voldatiuv)
3239 
3240 ! South boundary extrapolation
3241 .AND. WHERE(voldatiuv(:,1) /= rmiss voldatiuv(:,2) /= rmiss)
3242  tmp_arr(:,1) = voldatiuv(:,1) - (voldatiuv(:,2) - voldatiuv(:,1)) / 2.
3243  END WHERE
3244 
3245 ! Rest of the field
3246 .AND. WHERE(voldatiuv(:,1:this%griddim%dim%ny-1) /= rmiss &
3247  voldatiuv(:,2:this%griddim%dim%ny) /= rmiss)
3248  tmp_arr(:,2:this%griddim%dim%ny) = &
3249  (voldatiuv(:,1:this%griddim%dim%ny-1) + &
3250  voldatiuv(:,2:this%griddim%dim%ny)) / 2.
3251  END WHERE
3252 
3253  voldatiuv(:,:) = tmp_arr
3254  CALL volgrid_set_vol_2d(this, i, j, k, iv, voldatiuv)
3255  ENDDO
3256  ENDDO
3257  ENDDO
3258  ENDDO
3259 ENDIF ! tertium non datur
3260 
3261 .NOT.IF (ASSOCIATED(this%voldati)) THEN
3262  DEALLOCATE(voldatiuv)
3263 ENDIF
3264 DEALLOCATE (tmp_arr)
3265 
3266 end subroutine vg6d_c2a_mat
3267 
3268 
3269 !> \brief Display object on screen
3270 !!
3271 !! Show brief content on screen.
3272 subroutine display_volgrid6d (this)
3273 
3274 TYPE(volgrid6d),intent(in) :: this !< object to display
3275 integer :: i
3276 
3277 #ifdef DEBUG
3278 call l4f_category_log(this%category,L4F_DEBUG,"ora mostro gridinfo " )
3279 #endif
3280 
3281 print*,"----------------------- volgrid6d display ---------------------"
3282 call display(this%griddim)
3283 
3284 IF (ASSOCIATED(this%time))then
3285  print*,"---- time vector ----"
3286  print*,"elements=",size(this%time)
3287  do i=1, size(this%time)
3288  call display(this%time(i))
3289  end do
3290 end IF
3291 
3292 IF (ASSOCIATED(this%timerange))then
3293  print*,"---- timerange vector ----"
3294  print*,"elements=",size(this%timerange)
3295  do i=1, size(this%timerange)
3296  call display(this%timerange(i))
3297  end do
3298 end IF
3299 
3300 IF (ASSOCIATED(this%level))then
3301  print*,"---- level vector ----"
3302  print*,"elements=",size(this%level)
3303  do i=1, size(this%level)
3304  call display(this%level(i))
3305  end do
3306 end IF
3307 
3308 IF (ASSOCIATED(this%var))then
3309  print*,"---- var vector ----"
3310  print*,"elements=",size(this%var)
3311  do i=1, size(this%var)
3312  call display(this%var(i))
3313  end do
3314 end IF
3315 
3316 IF (ASSOCIATED(this%gaid))then
3317  print*,"---- gaid vector (present mask only) ----"
3318  print*,"elements=",shape(this%gaid)
3319  PRINT*,c_e(RESHAPE(this%gaid,(/SIZE(this%gaid)/)))
3320 end IF
3321 
3322 print*,"--------------------------------------------------------------"
3323 
3324 
3325 end subroutine display_volgrid6d
3326 
3327 
3328 !> \brief Display vector of object on screen
3329 !!
3330 !! Show brief content on screen.
3331 subroutine display_volgrid6dv (this)
3332 
3333 TYPE(volgrid6d),intent(in) :: this(:) !< vector of object to display
3334 integer :: i
3335 
3336 print*,"----------------------- volgrid6d vector ---------------------"
3337 
3338 print*,"elements=",size(this)
3339 
3340 do i=1, size(this)
3341 
3342 #ifdef DEBUG
3343  call l4f_category_log(this(i)%category,L4F_DEBUG,"ora mostro il vettore volgrid6d" )
3344 #endif
3345 
3346  call display(this(i))
3347 
3348 end do
3349 print*,"--------------------------------------------------------------"
3350 
3351 end subroutine display_volgrid6dv
3352 
3353 
3354 !> Reduce some dimensions (level and timerage) for semplification (rounding).
3355 !! Vector version of vg6dv_rounding
3356 subroutine vg6dv_rounding(vg6din,vg6dout,level,timerange,nostatproc,merge)
3357 type(volgrid6d),intent(in) :: vg6din(:) !< input volume
3358 type(volgrid6d),intent(out),pointer :: vg6dout(:) !> output volume
3359 type(vol7d_level),intent(in),optional :: level(:) !< almost equal level list
3360 type(vol7d_timerange),intent(in),optional :: timerange(:) !< almost equal timerange list
3361 logical,intent(in),optional :: merge !< if there are data on more then one almost equal levels or timeranges will be merged POINT BY POINT with priority for the fird data found ordered by icreasing var index
3362 logical,intent(in),optional :: nostatproc !< do not take in account statistical processing code in timerange and P2
3363 
3364 integer :: i
3365 
3366 allocate(vg6dout(size(vg6din)))
3367 
3368 do i = 1, size(vg6din)
3369  call vg6d_rounding(vg6din(i),vg6dout(i),level,timerange,nostatproc,merge)
3370 end do
3371 
3372 end subroutine vg6dv_rounding
3373 
3374 !> Reduce some dimensions (level and timerage) for semplification (rounding).
3375 !! You can use this for simplify and use variables in computation like alchimia
3376 !! where fields have to be on the same coordinate
3377 !! examples:
3378 !! means in time for short periods and istantaneous values
3379 !! 2 meter and surface levels
3380 !! If there are data on more then one almost equal levels or timeranges, the first var present (at least one point)
3381 !! will be taken (order is by icreasing var index).
3382 !! You can use predefined values for classic semplification
3383 !! almost_equal_levels and almost_equal_timeranges
3384 !! The level or timerange in output will be defined by the first element of level and timerange list
3385 subroutine vg6d_rounding(vg6din,vg6dout,level,timerange,nostatproc,merge)
3386 type(volgrid6d),intent(in) :: vg6din !< input volume
3387 type(volgrid6d),intent(out) :: vg6dout !> output volume
3388 type(vol7d_level),intent(in),optional :: level(:) !< almost equal level list
3389 type(vol7d_timerange),intent(in),optional :: timerange(:) !< almost equal timerange list
3390 logical,intent(in),optional :: merge !< if there are data on more then one almost equal levels or timeranges
3391 !! will be merged POINT BY POINT with priority for the fird data found ordered by icreasing var index
3392 logical,intent(in),optional :: nostatproc !< do not take in account statistical processing code in timerange and P2
3393 
3394 integer :: ilevel,itimerange
3395 type(vol7d_level) :: roundlevel(size(vg6din%level))
3396 type(vol7d_timerange) :: roundtimerange(size(vg6din%timerange))
3397 
3398 roundlevel=vg6din%level
3399 
3400 if (present(level))then
3401  do ilevel = 1, size(vg6din%level)
3402  if ((any(vg6din%level(ilevel) .almosteq. level))) then
3403  roundlevel(ilevel)=level(1)
3404  end if
3405  end do
3406 end if
3407 
3408 roundtimerange=vg6din%timerange
3409 
3410 if (present(timerange))then
3411  do itimerange = 1, size(vg6din%timerange)
3412  if ((any(vg6din%timerange(itimerange) .almosteq. timerange))) then
3413  roundtimerange(itimerange)=timerange(1)
3414  end if
3415  end do
3416 end if
3417 
3418 !set istantaneous values everywere
3419 !preserve p1 for forecast time
3420 if (optio_log(nostatproc)) then
3421  roundtimerange(:)%timerange=254
3422  roundtimerange(:)%p2=0
3423 end if
3424 
3425 
3426 call vg6d_reduce(vg6din,vg6dout,roundlevel,roundtimerange,merge)
3427 
3428 end subroutine vg6d_rounding
3429 
3430 !> Reduce some dimensions (level and timerage).
3431 !! You can pass a volume and specify duplicated levels and timeranges in roundlevel and roundtimerange;
3432 !! you get unique levels and timeranges in output.
3433 !! If there are data on equal levels or timeranges, the first var present (at least one point)
3434 !! will be taken (order is by icreasing var index).
3435 !! you can specify merge and if there are data on equal levels or timeranges will be merged POINT BY POINT
3436 !! with priority for the first data found ordered by icreasing var index (require to decode all the data)
3437 !! Data are decoded only if needed so the output should be with or without voldata allocated
3438 subroutine vg6d_reduce(vg6din,vg6dout,roundlevel,roundtimerange,merge)
3439 type(volgrid6d),intent(in) :: vg6din !< input volume
3440 type(volgrid6d),intent(out) :: vg6dout !< output volume
3441 type(vol7d_level),intent(in) :: roundlevel(:) !< new level list to use for rounding
3442 type(vol7d_timerange),intent(in) :: roundtimerange(:) !< new timerange list to use for rounding
3443 logical,intent(in),optional :: merge !< if there are data on equal levels or timeranges will be merged POINT BY POINT with priority for the first data found ordered by icreasing var index (require to decode all the data)
3444 
3445 integer :: nlevel,ntime,ntimerange,nvar,ilevel,itimerange,ivar,indl,indt,itime,nx,ny
3446 real,allocatable :: vol2d(:,:)
3447 
3448 nx=vg6din%griddim%dim%nx
3449 ny=vg6din%griddim%dim%ny
3450 nlevel=count_distinct(roundlevel,back=.true.)
3451 ntime=size(vg6din%time)
3452 ntimerange=count_distinct(roundtimerange,back=.true.)
3453 nvar=size(vg6din%var)
3454 
3455 call init(vg6dout, vg6din%griddim, vg6din%time_definition, categoryappend="generated by vg6d_reduce")
3456 call volgrid6d_alloc(vg6dout, vg6din%griddim%dim, ntime, nlevel, ntimerange, nvar)
3457 
3458 .or.if ( ASSOCIATED(vg6din%voldati) optio_log(merge)) then
3459  call volgrid6d_alloc_vol(vg6dout,inivol=.true.,decode=.true.)
3460  allocate(vol2d(nx,ny))
3461 else
3462  call volgrid6d_alloc_vol(vg6dout,inivol=.true.,decode=.false.)
3463 end if
3464 
3465 vg6dout%time=vg6din%time
3466 vg6dout%var=vg6din%var
3467 vg6dout%timerange=pack_distinct(roundtimerange,ntimerange,back=.true.)
3468 vg6dout%level=pack_distinct(roundlevel,nlevel,back=.true.)
3469 ! sort modified dimensions
3470 CALL sort(vg6dout%timerange)
3471 CALL sort(vg6dout%level)
3472 
3473 do ilevel=1,size(vg6din%level)
3474  indl=index(vg6dout%level,roundlevel(ilevel))
3475  do itimerange=1,size(vg6din%timerange)
3476  indt=index(vg6dout%timerange,roundtimerange(itimerange))
3477  do ivar=1, nvar
3478  do itime=1,ntime
3479 
3480  if ( ASSOCIATED(vg6din%voldati)) then
3481  vol2d=vg6din%voldati(:,:,ilevel,itime,itimerange,ivar)
3482  end if
3483 
3484  if (optio_log(merge)) then
3485 
3486 .not. if ( ASSOCIATED(vg6din%voldati)) then
3487  CALL grid_id_decode_data(vg6din%gaid(ilevel,itime,itimerange,ivar), vol2d)
3488  end if
3489 
3490  !! merge present data point by point
3491 .not. where ( c_e(vg6dout%voldati(:,:,indl,itime,indt,ivar)))
3492 
3493  vg6dout%voldati(:,:,indl,itime,indt,ivar)=vol2d
3494 
3495  end where
3496  else if ( ASSOCIATED(vg6din%voldati)) then
3497 .not. if ( any(c_e(vg6dout%voldati(:,:,indl,itime,indt,ivar))))then
3498  vg6dout%voldati(:,:,indl,itime,indt,ivar)=vol2d
3499  end if
3500  end if
3501 
3502 .and..not. if (c_e(vg6din%gaid(ilevel,itime,itimerange,ivar)) c_e(vg6dout%gaid(indl,itime,indt,ivar)))then
3503  call copy (vg6din%gaid(ilevel,itime,itimerange,ivar), vg6dout%gaid(indl,itime,indt,ivar))
3504  end if
3505  end do
3506  end do
3507  end do
3508 end do
3509 
3510 .or.if ( ASSOCIATED(vg6din%voldati) optio_log(merge)) then
3511  deallocate(vol2d)
3512 end if
3513 
3514 end subroutine vg6d_reduce
3515 
3516 
3517 end module volgrid6d_class
3518 
3519 
3520 
3521 !>\example example_vg6d_3.f90
3522 !!\brief Programma esempio semplice per gridinfo e volgrid6d.
3523 !!
3524 !! Programma che importa da file un vettore di gridinfo poi lo importa in volgrid6d. Da volgrid6d viene di nuovo creato un vettore di gridinfo per poi exportare su file.
3525 
3526 !>\example example_vg6d_5.f90
3527 !!\brief Programma trasformazione da volgrid6d a volgrid6d
3528 !!
3529 !! Legge grib da un file e li organizza in un vettore di strutture
3530 !! volgrid6d mettendoli a disposizione per eventuali elaborazioni;
3531 !! vengono poi riesportati a un file grib
3532 
3533 !>\example example_vg6d_8.f90
3534 !! \brief Programma scrittura su file vettore di anagrafica
3535 
3536 !>\example example_vg6d_6.f90
3537 !! \brief Programma trasformazione da volgrid6d a vol7d
3538 
3539 !>\example example_vg6d_7.f90
3540 !! \brief Programma trasformazione da vol7d a volgrid6d
3541 
3542 !>\example example_vg6d_9.f90
3543 !! \brief Example to create a grib editionNumber = 2 file from data generated in memory using a grib_api template
Functions that return a trimmed CHARACTER representation of the input variable.
Read the object from a formatted or unformatted file.
Write the object on a formatted or unformatted file.
Method for returning the contents of the object.
Clone the object, creating a new independent instance of the object exactly equal to the starting one...
Decode and return the data array from a grid_id object associated to a gridinfo object.
Encode a data array into a grid_id object associated to a gridinfo object.
Index method.
Emit log message for a category with specific priority.
Represent level object in a pretty string.
Destructor, it releases every information and memory buffer associated with the object.
Display on standard output a description of the volgrid6d object provided.
Export an object dirctly to a native file, to a gridinfo object or to a supported file format through...
Import an object dirctly from a native file, from a gridinfo object or from a supported file format t...
Constructor, it creates a new instance of the object.
Reduce some dimensions (level and timerage) for semplification (rounding).
Transform between any combination of volgrid6d and vol7d objects by means of a transform_def object d...
Apply the conversion function this to values.
This module defines usefull general purpose function and subroutine.
Classi per la gestione delle coordinate temporali.
Module for describing geographically referenced regular grids.
Definition: grid_class.F90:243
Module for defining the extension and coordinates of a rectangular georeferenced grid.
This module defines an abstract interface to different drivers for access to files containing gridded...
Module for defining transformations between rectangular georeferenced grids and between grids and spa...
Class for managing information about a single gridded georeferenced field, typically imported from an...
classe per la gestione del logging
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
Classe per la gestione di un volume completo di dati osservati.
classe per import ed export di volumi da e in DB-All.e
Classe per la gestione dei livelli verticali in osservazioni meteo e affini.
Classe per la gestione degli intervalli temporali di osservazioni meteo e affini.
This module defines objects and methods for managing data volumes on rectangular georeferenced grids.
Class for managing physical variables in a grib 1/2 fashion.
Class for expressing an absolute time value.
Derived type associated to a block/message/record/band of gridded data coming from a file-like object...
Derived type defining a dynamically extensible array of TYPE(gridinfo_def) elements.
Object describing a single gridded message/band.
Definisce un oggetto contenente i volumi anagrafica e dati e tutti i descrittori delle loro dimension...
Object describing a rectangular, homogeneous gridded dataset.

Generated with Doxygen.