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

Generated with Doxygen.