libsim Versione 7.2.6

◆ vg6d_rounding()

subroutine vg6d_rounding ( type(volgrid6d), intent(in) vg6din,
type(volgrid6d), intent(out) vg6dout,
type(vol7d_level), dimension(:), intent(in), optional level,
type(vol7d_timerange), dimension(:), intent(in), optional timerange,
logical, intent(in), optional nostatproc,
logical, intent(in), optional merge )

Reduce some dimensions (level and timerage) for semplification (rounding).

You can use this for simplify and use variables in computation like alchimia where fields have to be on the same coordinate examples: means in time for short periods and istantaneous values 2 meter and surface levels If there are data on more then one almost equal levels or timeranges, the first var present (at least one point) will be taken (order is by icreasing var index). You can use predefined values for classic semplification almost_equal_levels and almost_equal_timeranges The level or timerange in output will be defined by the first element of level and timerange list

Parametri
[in]vg6dininput volume
[in]levelalmost equal level list
[in]timerangealmost equal timerange list
[in]mergeif 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
[in]nostatprocdo not take in account statistical processing code in timerange and P2

Definizione alla linea 3600 del file volgrid6d_class.F90.

3601! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
3602! authors:
3603! Davide Cesari <dcesari@arpa.emr.it>
3604! Paolo Patruno <ppatruno@arpa.emr.it>
3605
3606! This program is free software; you can redistribute it and/or
3607! modify it under the terms of the GNU General Public License as
3608! published by the Free Software Foundation; either version 2 of
3609! the License, or (at your option) any later version.
3610
3611! This program is distributed in the hope that it will be useful,
3612! but WITHOUT ANY WARRANTY; without even the implied warranty of
3613! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3614! GNU General Public License for more details.
3615
3616! You should have received a copy of the GNU General Public License
3617! along with this program. If not, see <http://www.gnu.org/licenses/>.
3618#include "config.h"
3619!> \defgroup volgrid6d Libsim package, volgrid6d library.
3620!! The libsim volgrid6d library contains classes for managing data on
3621!! regular rectangular grids, tipically the output of a numerical
3622!! weather prediction model, and for their import from a WMO GRIB file
3623!! of from other geophysical data file formats. In order to compile
3624!! and link programs using this library, you have to insert the
3625!! required \c USE statements in the program units involved, specify
3626!! the location of module files when compiling (tipically \c
3627!! -I/usr/lib/gfortran/modules or \c -I/usr/lib64/gfortran/modules or
3628!! \c -I/usr/include) and indicate the library name \c -lsim_volgrid6d
3629!! when linking, assuming that the library has been installed in a
3630!! default location.
3631
3632!> This module defines objects and methods for managing data volumes
3633!! on rectangular georeferenced grids. The data are accomodated in a
3634!! multi-dimensional array with 6 predefined dimensions. Different
3635!! geographic coordinates and projections are supported, mainly
3636!! inspired by grib coding standards. The \a volgrid6d object contains
3637!! information and data on an homogeneous grid definition, while
3638!! different grids are managed as arrays of \a volgrid6d objects.
3639!! Every object contains also an identificator of the grid (\a grid_id
3640!! object), carrying information about the driver used or which has to
3641!! be used for import/export from/to file. With the help of \a
3642!! gridinfo_def class, data can be imported and exported to the
3643!! supported formats, mainly grib1 and grib2 through grib_api and many
3644!! GIS-style formats through gdal.
3645!!
3646!! Simple example program \include example_vg6d_3.f90
3647!! Example of transformation from volgrid6d to volgrid6d \include example_vg6d_5.f90
3648!! Example of transformation from volgrid6d to vol7d \include example_vg6d_6.f90
3649!! Example of transformation from vol7d to volgrid6d \include example_vg6d_7.f90
3650!!
3651!!\ingroup volgrid6d
3652MODULE volgrid6d_class
3655USE vol7d_class
3657USE geo_proj_class
3658USE grid_class
3662USE log4fortran
3664USE grid_id_class
3667!USE file_utilities
3668#ifdef HAVE_DBALLE
3670#endif
3671IMPLICIT NONE
3672
3673character (len=255),parameter:: subcategory="volgrid6d_class"
3674
3675!> Object describing a rectangular, homogeneous gridded dataset
3676type volgrid6d
3677 type(griddim_def) :: griddim !< grid descriptor
3678 TYPE(datetime),pointer :: time(:) !< time dimension descriptor
3679 TYPE(vol7d_timerange),pointer :: timerange(:) !< timerange (forecast, analysis, statistically processed) dimension descriptor
3680 TYPE(vol7d_level),pointer :: level(:) !< vertical level dimension descriptor
3681 TYPE(volgrid6d_var),pointer :: var(:) !< physical variable dimension descriptor
3682 TYPE(grid_id),POINTER :: gaid(:,:,:,:) !< array of grid identifiers, carrying information about the driver for import/export from/to file, indices are: (level,time,timerange,var)
3683 REAL,POINTER :: voldati(:,:,:,:,:,:) !< array of data, indices are: (x,y,level,time,timerange,var)
3684 integer :: time_definition !< time definition; 0=time is reference time ; 1=time is validity time
3685 integer :: category = 0 !< log4fortran category
3686end type volgrid6d
3687
3688!> Constructor, it creates a new instance of the object.
3689INTERFACE init
3690 MODULE PROCEDURE volgrid6d_init
3691END INTERFACE
3692
3693!> Destructor, it releases every information and memory buffer
3694!! associated with the object.
3695INTERFACE delete
3696 MODULE PROCEDURE volgrid6d_delete, volgrid6dv_delete
3697END INTERFACE
3698
3699!> Import an object dirctly from a native file, from a \a gridinfo object
3700!! or from a supported file format through a \a gridinfo object.
3701INTERFACE import
3702 MODULE PROCEDURE volgrid6d_read_from_file
3703 MODULE PROCEDURE import_from_gridinfo, import_from_gridinfovv, &
3704 volgrid6d_import_from_file
3705END INTERFACE
3706
3707!> Export an object dirctly to a native file, to a \a gridinfo object
3708!! or to a supported file format through a \a gridinfo object.
3709INTERFACE export
3710 MODULE PROCEDURE volgrid6d_write_on_file
3711 MODULE PROCEDURE export_to_gridinfo, export_to_gridinfov, export_to_gridinfovv,&
3712 volgrid6d_export_to_file
3713END INTERFACE
3714
3715! methods for computing transformations through an initialised
3716! grid_transform object, probably too low level to be interfaced
3717INTERFACE compute
3718 MODULE PROCEDURE volgrid6d_transform_compute, volgrid6d_v7d_transform_compute,&
3719 v7d_volgrid6d_transform_compute, v7d_v7d_transform_compute
3720END INTERFACE
3721
3722!> Transform between any combination of \a volgrid6d and \a vol7d objects
3723!! by means of a \a transform_def object describing the transformation.
3724INTERFACE transform
3725 MODULE PROCEDURE volgrid6d_transform, volgrid6dv_transform,&
3726 volgrid6d_v7d_transform, volgrid6dv_v7d_transform, v7d_volgrid6d_transform, &
3727 v7d_v7d_transform
3728END INTERFACE
3729
3730INTERFACE wind_rot
3731 MODULE PROCEDURE vg6d_wind_rot
3732END INTERFACE
3733
3734INTERFACE wind_unrot
3735 MODULE PROCEDURE vg6d_wind_unrot
3736END INTERFACE
3737
3738!> Display on standard output a description of the \a volgrid6d object
3739!! provided.
3740INTERFACE display
3741 MODULE PROCEDURE display_volgrid6d,display_volgrid6dv
3742END INTERFACE
3743
3744!> Reduce some dimensions (level and timerage) for semplification (rounding).
3745!! You can use this for simplify and use variables in computation like alchimia
3746!! where fields have to be on the same coordinate
3747!! examples:
3748!! means in time for short periods and istantaneous values
3749!! 2 meter and surface levels
3750!! If there are data on more then one almost equal levels or timeranges, the first var present (at least one point)
3751!! will be taken (order is by icreasing var index).
3752!! You can use predefined values for classic semplification
3753!! almost_equal_levels and almost_equal_timeranges
3754!! The level or timerange in output will be defined by the first element of level and timerange list
3755INTERFACE rounding
3756 MODULE PROCEDURE vg6d_rounding, vg6dv_rounding
3757END INTERFACE
3758
3759private
3760
3761PUBLIC volgrid6d,init,delete,export,import,compute,transform, &
3762 wind_rot,wind_unrot,vg6d_c2a,display,volgrid6d_alloc,volgrid6d_alloc_vol, &
3763 volgrid_get_vol_2d, volgrid_set_vol_2d, volgrid_get_vol_3d, volgrid_set_vol_3d
3764PUBLIC rounding, vg6d_reduce
3765
3766CONTAINS
3767
3768
3769!> Constructor, it creates a new instance of the object.
3770!! The constructor should be explicitly used only in rare cases,
3771!! \a volgrid6d objects are usually created through the \a import
3772!! interface.
3773SUBROUTINE volgrid6d_init(this, griddim, time_definition, categoryappend)
3774TYPE(volgrid6d) :: this !< object to be initialized
3775TYPE(griddim_def),OPTIONAL :: griddim !< grid descriptor
3776INTEGER,INTENT(IN),OPTIONAL :: time_definition !< 0=time is reference time; 1=time is validity time
3777CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
3778
3779character(len=512) :: a_name
3780
3781if (present(categoryappend))then
3782 call l4f_launcher(a_name,a_name_append=trim(subcategory)//"."//trim(categoryappend))
3783else
3784 call l4f_launcher(a_name,a_name_append=trim(subcategory))
3785endif
3786this%category=l4f_category_get(a_name)
3787
3788#ifdef DEBUG
3789call l4f_category_log(this%category,l4f_debug,"init")
3790#endif
3791
3792call init(this%griddim)
3793
3794if (present(griddim))then
3795 call copy (griddim,this%griddim)
3796end if
3797
3798CALL vol7d_var_features_init() ! initialise var features table once
3799
3800if(present(time_definition)) then
3801 this%time_definition = time_definition
3802else
3803 this%time_definition = 0 !default to reference time
3804end if
3805
3806nullify (this%time,this%timerange,this%level,this%var)
3807nullify (this%gaid,this%voldati)
3808
3809END SUBROUTINE volgrid6d_init
3810
3811
3812!> Allocate the dimension descriptors of the \a volgrid6d object.
3813!! This method allocates the horizontal grid descriptor and the one
3814!! dimensional arrays of the dimensions
3815!! - time
3816!! - vertical level
3817!! - timerange
3818!! - physical variable
3819!!
3820!! This method should be explicitly used only in rare cases, it is
3821!! usually called implicitly through the \a import interface.
3822SUBROUTINE volgrid6d_alloc(this, dim, ntime, nlevel, ntimerange, nvar, ini)
3823TYPE(volgrid6d),INTENT(inout) :: this !< object whose decriptors should be allocated
3824TYPE(grid_dim),INTENT(in),OPTIONAL :: dim !< horizontal grid size X, Y
3825INTEGER,INTENT(in),OPTIONAL :: ntime !< number of time levels
3826INTEGER,INTENT(in),OPTIONAL :: nlevel !< number of vertical levels
3827INTEGER,INTENT(in),OPTIONAL :: ntimerange !< number of different timeranges
3828INTEGER,INTENT(in),OPTIONAL :: nvar !< number of physical variables
3829LOGICAL,INTENT(in),OPTIONAL :: ini !< if provided and \c .TRUE., for each allocated dimension descriptor the constructor is called without extra parameters, thus initializing everything as missing value
3830
3831INTEGER :: i, stallo
3832LOGICAL :: linit
3833
3834#ifdef DEBUG
3835call l4f_category_log(this%category,l4f_debug,"alloc")
3836#endif
3837
3838IF (PRESENT(ini)) THEN
3839 linit = ini
3840ELSE
3841 linit = .false.
3842ENDIF
3843
3844
3845if (present(dim)) call copy (dim,this%griddim%dim)
3846
3847
3848IF (PRESENT(ntime)) THEN
3849 IF (ntime >= 0) THEN
3850 IF (ASSOCIATED(this%time)) DEALLOCATE(this%time)
3851#ifdef DEBUG
3852 call l4f_category_log(this%category,l4f_debug,"alloc ntime "//to_char(ntime))
3853#endif
3854 ALLOCATE(this%time(ntime),stat=stallo)
3855 if (stallo /=0)then
3856 call l4f_category_log(this%category,l4f_fatal,"allocating memory")
3857 CALL raise_fatal_error()
3858 end if
3859 IF (linit) THEN
3860 DO i = 1, ntime
3861 this%time(i) = datetime_miss
3862! CALL init(this%time(i)) ! senza argomento inizializza a zero non missing
3863 ! baco di datetime?
3864 ENDDO
3865 ENDIF
3866 ENDIF
3867ENDIF
3868IF (PRESENT(nlevel)) THEN
3869 IF (nlevel >= 0) THEN
3870 IF (ASSOCIATED(this%level)) DEALLOCATE(this%level)
3871#ifdef DEBUG
3872 call l4f_category_log(this%category,l4f_debug,"alloc nlevel "//to_char(nlevel))
3873#endif
3874 ALLOCATE(this%level(nlevel),stat=stallo)
3875 if (stallo /=0)then
3876 call l4f_category_log(this%category,l4f_fatal,"allocating memory")
3877 CALL raise_fatal_error()
3878 end if
3879 IF (linit) THEN
3880 DO i = 1, nlevel
3881 CALL init(this%level(i))
3882 ENDDO
3883 ENDIF
3884 ENDIF
3885ENDIF
3886IF (PRESENT(ntimerange)) THEN
3887 IF (ntimerange >= 0) THEN
3888 IF (ASSOCIATED(this%timerange)) DEALLOCATE(this%timerange)
3889#ifdef DEBUG
3890 call l4f_category_log(this%category,l4f_debug,"alloc ntimerange "//to_char(ntimerange))
3891#endif
3892 ALLOCATE(this%timerange(ntimerange),stat=stallo)
3893 if (stallo /=0)then
3894 call l4f_category_log(this%category,l4f_fatal,"allocating memory")
3895 CALL raise_fatal_error()
3896 end if
3897 IF (linit) THEN
3898 DO i = 1, ntimerange
3899 CALL init(this%timerange(i))
3900 ENDDO
3901 ENDIF
3902 ENDIF
3903ENDIF
3904IF (PRESENT(nvar)) THEN
3905 IF (nvar >= 0) THEN
3906 IF (ASSOCIATED(this%var)) DEALLOCATE(this%var)
3907#ifdef DEBUG
3908 call l4f_category_log(this%category,l4f_debug,"alloc nvar "//to_char(nvar))
3909#endif
3910 ALLOCATE(this%var(nvar),stat=stallo)
3911 if (stallo /=0)then
3912 call l4f_category_log(this%category,l4f_fatal,"allocating memory")
3913 CALL raise_fatal_error()
3914 end if
3915 IF (linit) THEN
3916 DO i = 1, nvar
3917 CALL init(this%var(i))
3918 ENDDO
3919 ENDIF
3920 ENDIF
3921ENDIF
3922
3923end SUBROUTINE volgrid6d_alloc
3924
3925
3926!> Allocate the data array of the \a volgrid6d object.
3927!! This method allocates the main 6-dimensional data array
3928!! \a this%voldati and the 4-dimensional \a grid_id array \a this%gaid
3929!! with a shape dictated by the previous call(s) to \a vol7d_alloc().
3930!! if any descriptor (except horizontal grid) has not been allocated
3931!! yet, it is allocated here with a size of 1. This method should be
3932!! explicitly used only in rare cases, it is usually called implicitly
3933!! through the \a import interface.
3934SUBROUTINE volgrid6d_alloc_vol(this, ini, inivol, decode)
3935TYPE(volgrid6d),INTENT(inout) :: this !< object whose decriptors should be allocated
3936LOGICAL,INTENT(in),OPTIONAL :: ini !< if provided and \c .TRUE., for each dimension descriptor not yet allocated and allocated here the constructor is called without extra parameters, thus initializing the element as missing value
3937LOGICAL,INTENT(in),OPTIONAL :: inivol !< if provided and \c .FALSE., the allocated volumes will not be initialized to missing values
3938LOGICAL,INTENT(in),OPTIONAL :: decode !< if provided and \c .TRUE., the \a this%voldati volume is allocated, otherwise only \a this%gaid will be allocated
3939
3940INTEGER :: stallo
3941LOGICAL :: linivol
3942
3943#ifdef DEBUG
3944call l4f_category_log(this%category,l4f_debug,"start alloc_vol")
3945#endif
3946
3947IF (PRESENT(inivol)) THEN ! opposite condition, cannot use optio_log
3948 linivol = inivol
3949ELSE
3950 linivol = .true.
3951ENDIF
3952
3953IF (this%griddim%dim%nx > 0 .AND. this%griddim%dim%ny > 0) THEN
3954! allocate dimension descriptors to a minimal size for having a
3955! non-null volume
3956 IF (.NOT. ASSOCIATED(this%var)) CALL volgrid6d_alloc(this, nvar=1, ini=ini)
3957 IF (.NOT. ASSOCIATED(this%time)) CALL volgrid6d_alloc(this, ntime=1, ini=ini)
3958 IF (.NOT. ASSOCIATED(this%level)) CALL volgrid6d_alloc(this, nlevel=1, ini=ini)
3959 IF (.NOT. ASSOCIATED(this%timerange)) CALL volgrid6d_alloc(this, ntimerange=1, ini=ini)
3960
3961 IF (optio_log(decode)) THEN
3962 IF (.NOT.ASSOCIATED(this%voldati)) THEN
3963#ifdef DEBUG
3964 CALL l4f_category_log(this%category,l4f_debug,"alloc voldati volume")
3965#endif
3966
3967 ALLOCATE(this%voldati(this%griddim%dim%nx,this%griddim%dim%ny,&
3968 SIZE(this%level), SIZE(this%time), &
3969 SIZE(this%timerange), SIZE(this%var)),stat=stallo)
3970 IF (stallo /= 0)THEN
3971 CALL l4f_category_log(this%category,l4f_fatal,"allocating memory")
3972 CALL raise_fatal_error()
3973 ENDIF
3974
3975! this is not really needed if we can check other flags for a full
3976! field missing values
3977 IF (linivol) this%voldati = rmiss
3978 this%voldati = rmiss
3979 ENDIF
3980 ENDIF
3981
3982 IF (.NOT.ASSOCIATED(this%gaid)) THEN
3983#ifdef DEBUG
3984 CALL l4f_category_log(this%category,l4f_debug,"alloc gaid volume")
3985#endif
3986 ALLOCATE(this%gaid(SIZE(this%level), SIZE(this%time), &
3987 SIZE(this%timerange), SIZE(this%var)),stat=stallo)
3988 IF (stallo /= 0)THEN
3989 CALL l4f_category_log(this%category,l4f_fatal,"allocating memory")
3990 CALL raise_fatal_error()
3991 ENDIF
3992
3993 IF (linivol) THEN
3994!!$ DO i=1 ,SIZE(this%gaid,1)
3995!!$ DO ii=1 ,SIZE(this%gaid,2)
3996!!$ DO iii=1 ,SIZE(this%gaid,3)
3997!!$ DO iiii=1 ,SIZE(this%gaid,4)
3998!!$ this%gaid(i,ii,iii,iiii) = grid_id_new() ! optimize?
3999!!$ ENDDO
4000!!$ ENDDO
4001!!$ ENDDO
4002!!$ ENDDO
4003
4004 this%gaid = grid_id_new()
4005 ENDIF
4006 ENDIF
4007
4008ELSE
4009 CALL l4f_category_log(this%category,l4f_fatal,'volgrid6d_alloc_vol: &
4010 &trying to allocate a volume with an invalid or unspecified horizontal grid')
4011 CALL raise_fatal_error()
4012ENDIF
4013
4014END SUBROUTINE volgrid6d_alloc_vol
4015
4016
4017!> Return a 2-d pointer to a x-y slice of a volume. This method works
4018!! both with volumes having allocated and non-allocated this%voldati
4019!! array, and it returns a pointer to a 2-d slice either from the
4020!! allocated this%voldati array or from the grid_id object on file or
4021!! in memory. In the second case the pointer should be either
4022!! ALLOCATE'd to the expected size or NULLIFY'ed, and if NULLIFY'ed,
4023!! it is allocated within the method, thus it will have to be
4024!! deallocated by the caller when not in use anymore. Since this
4025!! method may be called many times by a program, it is optimized for
4026!! speed and it does not make any check about the matching size of the
4027!! pointer and the array or about the allocation status of \a this, so
4028!! it should be called only when everything has been checked to be in
4029!! good shape.
4030SUBROUTINE volgrid_get_vol_2d(this, ilevel, itime, itimerange, ivar, voldati)
4031TYPE(volgrid6d),INTENT(in) :: this !< object from which the slice has to be retrieved
4032INTEGER,INTENT(in) :: ilevel !< index of vertical level of the slice
4033INTEGER,INTENT(in) :: itime !< index of time level of the slice
4034INTEGER,INTENT(in) :: itimerange !< index of timerange of the slice
4035INTEGER,INTENT(in) :: ivar !< index of physical variable of the slice
4036REAL,POINTER :: voldati(:,:) !< pointer to the data, if \a this%voldati is already allocated, it will just point to the requested slice, otherwise it will be allocated if and only if it is nullified on entry
4037
4038IF (ASSOCIATED(this%voldati)) THEN
4039 voldati => this%voldati(:,:,ilevel,itime,itimerange,ivar)
4040 RETURN
4041ELSE
4042 IF (.NOT.ASSOCIATED(voldati)) THEN
4043 ALLOCATE(voldati(this%griddim%dim%nx,this%griddim%dim%ny))
4044 ENDIF
4045 CALL grid_id_decode_data(this%gaid(ilevel,itime,itimerange,ivar), voldati)
4046ENDIF
4047
4048END SUBROUTINE volgrid_get_vol_2d
4049
4050
4051!> Return a 3-d pointer to a x-y-z slice of a volume. This method works
4052!! both with volumes having allocated and non-allocated this%voldati
4053!! array, and it returns a pointer to a 3-d slice either from the
4054!! allocated this%voldati array or from the grid_id object on file or
4055!! in memory. In the second case the pointer should be either
4056!! ALLOCATE'd to the expected size or NULLIFY'ed, and if NULLIFY'ed,
4057!! it is allocated within the method, thus it will have to be
4058!! deallocated by the caller when not in use anymore. Since this
4059!! method may be called many times by a program, it is optimized for
4060!! speed and it does not make any check about the matching size of the
4061!! pointer and the array or about the allocation status of \a this, so
4062!! it should be called only when everything has been checked to be in
4063!! good shape.
4064SUBROUTINE volgrid_get_vol_3d(this, itime, itimerange, ivar, voldati)
4065TYPE(volgrid6d),INTENT(in) :: this !< object from which the slice has to be retrieved
4066INTEGER,INTENT(in) :: itime !< index of time level of the slice
4067INTEGER,INTENT(in) :: itimerange !< index of timerange of the slice
4068INTEGER,INTENT(in) :: ivar !< index of physical variable of the slice
4069REAL,POINTER :: voldati(:,:,:) !< pointer to the data, if \a this%voldati is already allocated, it will just point to the requested slice, otherwise it will be allocated if and only if it is nullified on entry
4070
4071INTEGER :: ilevel
4072
4073IF (ASSOCIATED(this%voldati)) THEN
4074 voldati => this%voldati(:,:,:,itime,itimerange,ivar)
4075 RETURN
4076ELSE
4077 IF (.NOT.ASSOCIATED(voldati)) THEN
4078 ALLOCATE(voldati(this%griddim%dim%nx,this%griddim%dim%ny,SIZE(this%level)))
4079 ENDIF
4080!$OMP PARALLEL DEFAULT(SHARED)
4081!$OMP MASTER
4082 DO ilevel = 1, SIZE(this%level)
4083!$OMP TASK FIRSTPRIVATE(ilevel)
4084 CALL grid_id_decode_data(this%gaid(ilevel,itime,itimerange,ivar), &
4085 voldati(:,:,ilevel))
4086!$OMP END TASK
4087 ENDDO
4088!$OMP END MASTER
4089!$OMP END PARALLEL
4090ENDIF
4091
4092END SUBROUTINE volgrid_get_vol_3d
4093
4094
4095!> Reset a 2-d x-y slice of a volume after the data have been modified.
4096!! This method works both with volumes having allocated and
4097!! non-allocated this%voldati array, and it updates the requested
4098!! slice. In case \a this%voldati is already allocated, this is a
4099!! no-operation while in the other case this method encodes the field
4100!! provided into the grid_id object on file or in memory. Since this
4101!! method may be called many times by a program, it is optimized for
4102!! speed and it does not make any check about the matching size of the
4103!! field and the array or about the allocation status of \a this, so
4104!! it should be called only when everything has been checked to be in
4105!! good shape.
4106SUBROUTINE volgrid_set_vol_2d(this, ilevel, itime, itimerange, ivar, voldati)
4107TYPE(volgrid6d),INTENT(inout) :: this !< object in which slice has to be updated
4108INTEGER,INTENT(in) :: ilevel !< index of vertical level of the slice
4109INTEGER,INTENT(in) :: itime !< index of time level of the slice
4110INTEGER,INTENT(in) :: itimerange !< index of timerange of the slice
4111INTEGER,INTENT(in) :: ivar !< index of physical variable of the slice
4112REAL,INTENT(in) :: voldati(:,:) !< updated values of the slice
4113
4114IF (ASSOCIATED(this%voldati)) THEN
4115 RETURN
4116ELSE
4117 CALL grid_id_encode_data(this%gaid(ilevel,itime,itimerange,ivar), voldati)
4118ENDIF
4119
4120END SUBROUTINE volgrid_set_vol_2d
4121
4122
4123!> Reset a 3-d x-y-z slice of a volume after the data have been modified.
4124!! This method works both with volumes having allocated and
4125!! non-allocated this%voldati array, and it updates the requested
4126!! slice. In case \a this%voldati is already allocated, this is a
4127!! no-operation while in the other case this method encodes the field
4128!! provided into the grid_id object on file or in memory. Since this
4129!! method may be called many times by a program, it is optimized for
4130!! speed and it does not make any check about the matching size of the
4131!! field and the array or about the allocation status of \a this, so
4132!! it should be called only when everything has been checked to be in
4133!! good shape.
4134SUBROUTINE volgrid_set_vol_3d(this, itime, itimerange, ivar, voldati)
4135TYPE(volgrid6d),INTENT(inout) :: this !< object in which slice has to be updated
4136INTEGER,INTENT(in) :: itime !< index of time level of the slice
4137INTEGER,INTENT(in) :: itimerange !< index of timerange of the slice
4138INTEGER,INTENT(in) :: ivar !< index of physical variable of the slice
4139REAL,INTENT(in) :: voldati(:,:,:) !< updated values of the slice
4140
4141INTEGER :: ilevel
4142
4143IF (ASSOCIATED(this%voldati)) THEN
4144 RETURN
4145ELSE
4146!$OMP PARALLEL DEFAULT(SHARED)
4147!$OMP MASTER
4148 DO ilevel = 1, SIZE(this%level)
4149!$OMP TASK FIRSTPRIVATE(ilevel)
4150 CALL grid_id_encode_data(this%gaid(ilevel,itime,itimerange,ivar), &
4151 voldati(:,:,ilevel))
4152!$OMP END TASK
4153 ENDDO
4154!$OMP END MASTER
4155!$OMP END PARALLEL
4156ENDIF
4157
4158END SUBROUTINE volgrid_set_vol_3d
4159
4160
4161!> Destructor, it releases every information and memory buffer
4162!! associated with the object. It should be called also for objects
4163!! crated through the \a import interface.
4164SUBROUTINE volgrid6d_delete(this)
4165TYPE(volgrid6d),INTENT(inout) :: this
4166
4167INTEGER :: i, ii, iii, iiii
4168
4169#ifdef DEBUG
4170call l4f_category_log(this%category,l4f_debug,"delete")
4171#endif
4172
4173if (associated(this%gaid))then
4174
4175 DO i=1 ,SIZE(this%gaid,1)
4176 DO ii=1 ,SIZE(this%gaid,2)
4177 DO iii=1 ,SIZE(this%gaid,3)
4178 DO iiii=1 ,SIZE(this%gaid,4)
4179 CALL delete(this%gaid(i,ii,iii,iiii))
4180 ENDDO
4181 ENDDO
4182 ENDDO
4183 ENDDO
4184 DEALLOCATE(this%gaid)
4185
4186end if
4187
4188call delete(this%griddim)
4189
4190! call delete(this%time)
4191! call delete(this%timerange)
4192! call delete(this%level)
4193! call delete(this%var)
4194
4195if (associated( this%time )) deallocate(this%time)
4196if (associated( this%timerange )) deallocate(this%timerange)
4197if (associated( this%level )) deallocate(this%level)
4198if (associated( this%var )) deallocate(this%var)
4199
4200if (associated(this%voldati))deallocate(this%voldati)
4201
4202
4203 !chiudo il logger
4204call l4f_category_delete(this%category)
4205
4206END SUBROUTINE volgrid6d_delete
4207
4208
4209!>\brief Scrittura su file di un volume Volgrid6d.
4210!! Scrittura su file unformatted di un intero volume Volgrid6d.
4211!! Il volume viene serializzato e scritto su file.
4212!! Il file puo' essere aperto opzionalmente dall'utente. Si possono passare
4213!! opzionalmente unità e nome del file altrimenti assegnati internamente a dei default; se assegnati internamente
4214!! tali parametri saranno in output.
4215!! Se non viene fornito il nome file viene utilizzato un file di default con nome pari al nome del programma in
4216!! esecuzione con postfisso ".vg6d".
4217!! Come parametro opzionale c'è la description che insieme alla data corrente viene inserita nell'header del file.
4218subroutine volgrid6d_write_on_file (this,unit,description,filename,filename_auto)
4219
4220TYPE(volgrid6d),INTENT(IN) :: this !< volume volgrid6d da scrivere
4221integer,optional,intent(inout) :: unit !< unità su cui scrivere; se passata =0 ritorna il valore rielaborato (default =rielaborato internamente con getlun )
4222character(len=*),intent(in),optional :: filename !< nome del file su cui scrivere
4223character(len=*),intent(out),optional :: filename_auto !< nome del file generato se "filename" è omesso
4224character(len=*),INTENT(IN),optional :: description !< descrizione del volume
4225
4226integer :: lunit
4227character(len=254) :: ldescription,arg,lfilename
4228integer :: ntime, ntimerange, nlevel, nvar
4229integer :: tarray(8)
4230logical :: opened,exist
4231
4232#ifdef DEBUG
4233call l4f_category_log(this%category,l4f_debug,"write on file")
4234#endif
4235
4236ntime=0
4237ntimerange=0
4238nlevel=0
4239nvar=0
4240
4241!call idate(im,id,iy)
4242call date_and_time(values=tarray)
4243call getarg(0,arg)
4244
4245if (present(description))then
4246 ldescription=description
4247else
4248 ldescription="Volgrid6d generated by: "//trim(arg)
4249end if
4250
4251if (.not. present(unit))then
4252 lunit=getunit()
4253else
4254 if (unit==0)then
4255 lunit=getunit()
4256 unit=lunit
4257 else
4258 lunit=unit
4259 end if
4260end if
4261
4262lfilename=trim(arg)//".vg6d"
4263if (index(arg,'/',back=.true.) > 0) lfilename=lfilename(index(arg,'/',back=.true.)+1 : )
4264
4265if (present(filename))then
4266 if (filename /= "")then
4267 lfilename=filename
4268 end if
4269end if
4270
4271if (present(filename_auto))filename_auto=lfilename
4272
4273
4274inquire(unit=lunit,opened=opened)
4275if (.not. opened) then
4276 inquire(file=lfilename,exist=exist)
4277 if (exist) CALL raise_error('file exist; cannot open new file')
4278 if (.not.exist) open (unit=lunit,file=lfilename,form="UNFORMATTED")
4279 !print *, "opened: ",lfilename
4280end if
4281
4282if (associated(this%time)) ntime=size(this%time)
4283if (associated(this%timerange)) ntimerange=size(this%timerange)
4284if (associated(this%level)) nlevel=size(this%level)
4285if (associated(this%var)) nvar=size(this%var)
4286
4287
4288write(unit=lunit)ldescription
4289write(unit=lunit)tarray
4290
4291call write_unit( this%griddim,lunit)
4292write(unit=lunit) ntime, ntimerange, nlevel, nvar
4293
4294!! prime 4 dimensioni
4295if (associated(this%time)) call write_unit(this%time, lunit)
4296if (associated(this%level)) write(unit=lunit)this%level
4297if (associated(this%timerange)) write(unit=lunit)this%timerange
4298if (associated(this%var)) write(unit=lunit)this%var
4299
4300
4301!! Volumi di valori dati
4302
4303if (associated(this%voldati)) write(unit=lunit)this%voldati
4304
4305if (.not. present(unit)) close(unit=lunit)
4306
4307end subroutine volgrid6d_write_on_file
4308
4309
4310!>\brief Lettura da file di un volume Volgrid6d.
4311!! Lettura da file unformatted di un intero volume Volgrid6d.
4312!! Questa subroutine comprende volgrid6d_alloc e volgrid6d_alloc_vol.
4313!! Il file puo' essere aperto opzionalmente dall'utente. Si possono passare
4314!! opzionalmente unità e nome del file altrimenti assegnati internamente a dei default; se assegnati internamente
4315!! tali parametri saranno in output.
4316subroutine volgrid6d_read_from_file (this,unit,filename,description,tarray,filename_auto)
4317
4318TYPE(volgrid6d),INTENT(OUT) :: this !< Volume volgrid6d da leggere
4319integer,intent(inout),optional :: unit !< unità su cui è stato aperto un file; se =0 rielaborato internamente (default = elaborato internamente con getunit)
4320character(len=*),INTENT(in),optional :: filename !< nome del file eventualmente da aprire (default = nome dell'eseguibile.v7d)
4321character(len=*),intent(out),optional :: filename_auto !< nome del file generato se "filename" è omesso
4322character(len=*),INTENT(out),optional :: description !< descrizione del volume letto
4323integer,intent(out),optional :: tarray(8) !< vettore come definito da "date_and_time" della data di scrittura del volume
4324
4325integer :: ntime, ntimerange, nlevel, nvar
4326
4327character(len=254) :: ldescription,lfilename,arg
4328integer :: ltarray(8),lunit
4329logical :: opened,exist
4330
4331#ifdef DEBUG
4332call l4f_category_log(this%category,l4f_debug,"read from file")
4333#endif
4334
4335call getarg(0,arg)
4336
4337if (.not. present(unit))then
4338 lunit=getunit()
4339else
4340 if (unit==0)then
4341 lunit=getunit()
4342 unit=lunit
4343 else
4344 lunit=unit
4345 end if
4346end if
4347
4348lfilename=trim(arg)//".vg6d"
4349if (index(arg,'/',back=.true.) > 0) lfilename=lfilename(index(arg,'/',back=.true.)+1 : )
4350
4351if (present(filename))then
4352 if (filename /= "")then
4353 lfilename=filename
4354 end if
4355end if
4356
4357if (present(filename_auto))filename_auto=lfilename
4358
4359
4360inquire(unit=lunit,opened=opened)
4361if (.not. opened) then
4362 inquire(file=lfilename,exist=exist)
4363 IF (.NOT. exist) CALL raise_fatal_error('file '//trim(lfilename)//' does not exist, cannot open')
4364 open (unit=lunit,file=lfilename,form="UNFORMATTED")
4365end if
4366
4367read(unit=lunit)ldescription
4368read(unit=lunit)ltarray
4369
4370call l4f_log(l4f_info,"Info: reading volgrid6d from file: "//trim(lfilename))
4371call l4f_log(l4f_info,"Info: description: "//trim(ldescription))
4372!call l4f_log("Info: written on ",ltarray)
4373
4374if (present(description))description=ldescription
4375if (present(tarray))tarray=ltarray
4376
4377
4378call read_unit( this%griddim,lunit)
4379read(unit=lunit) ntime, ntimerange, nlevel, nvar
4380
4381
4382call volgrid6d_alloc (this, &
4383 ntime=ntime, ntimerange=ntimerange, nlevel=nlevel, nvar=nvar)
4384
4385call volgrid6d_alloc_vol (this)
4386
4387if (associated(this%time)) call read_unit(this%time, lunit)
4388if (associated(this%level)) read(unit=lunit)this%level
4389if (associated(this%timerange)) read(unit=lunit)this%timerange
4390if (associated(this%var)) read(unit=lunit)this%var
4391
4392
4393!! Volumi di valori
4394
4395if (associated(this%voldati)) read(unit=lunit)this%voldati
4396
4397if (.not. present(unit)) close(unit=lunit)
4398
4399end subroutine volgrid6d_read_from_file
4400
4401
4402!> Import a single \a gridinfo object into a \a volgrid6d object.
4403!! This methods imports a single gridded field from a \a gridinfo
4404!! object into a \a volgrid6d object, inserting it into the
4405!! multidimensional structure of volgrid6d. The volgrid6d object must
4406!! have been already initialized and the dimensions specified with
4407!! volgrid6d_alloc(). If the \a force argument is missing or \a
4408!! .FALSE. , the volgrid6d object dimension descriptors (time,
4409!! timerange, vertical level, physical variable) must already have
4410!! space for the corresponding values coming from gridinfo, otherwise
4411!! the object will be rejected; this means that all the volgrid6d
4412!! dimension descriptors should be correctly assigned. If \a force is
4413!! \a .TRUE. , the gridinfo dimension descriptors that do not fit into
4414!! available descriptors in the \a volgrid6d structure, will be
4415!! accomodated in a empty (i.e. equal to missing value) descriptor, if
4416!! available, otherwise the gridinfo will be rejected. The
4417!! descriptor of the grid in the \a volgrid object is assigned to the
4418!! descriptor contained in \a gridinfo if it is missing in \a volgrid,
4419!! otherwise it is checked and the object is rejected if grids do not
4420!! match.
4421SUBROUTINE import_from_gridinfo(this, gridinfo, force, dup_mode, clone, &
4422 isanavar)
4423TYPE(volgrid6d),INTENT(inout) :: this !< object in which to import
4424TYPE(gridinfo_def),INTENT(in) :: gridinfo !< gridinfo object to be imported
4425LOGICAL,INTENT(in),OPTIONAL :: force !< if provided and \c .TRUE., the gridinfo is forced into an empty element of \a this, if required and possible
4426INTEGER,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
4427LOGICAL , INTENT(in),OPTIONAL :: clone !< if provided and \c .TRUE. , clone the gaid's from \a gridinfo to \a this
4428LOGICAL,INTENT(IN),OPTIONAL :: isanavar !< if provides and \a .TRUE., the gridinfo object is treated as time-independent and replicated for every time and timerange
4429
4430CHARACTER(len=255) :: type
4431INTEGER :: itime0, itimerange0, itime1, itimerange1, itime, itimerange, &
4432 ilevel, ivar, ldup_mode
4433LOGICAL :: dup
4434TYPE(datetime) :: correctedtime
4435TYPE(vol7d_timerange) :: correctedtimerange
4436REAL,ALLOCATABLE :: tmpgrid(:,:)
4437
4438IF (PRESENT(dup_mode)) THEN
4439 ldup_mode = dup_mode
4440ELSE
4441 ldup_mode = 0
4442ENDIF
4443
4444call get_val(this%griddim,proj_type=type)
4445
4446#ifdef DEBUG
4447call l4f_category_log(this%category,l4f_debug,"import_from_gridinfo: "//trim(type))
4448#endif
4449
4450if (.not. c_e(type))then
4451 call copy(gridinfo%griddim, this%griddim)
4452! ho gia` fatto init, altrimenti non potrei fare la get_val(this%griddim)
4453! per cui meglio non ripetere
4454! call init(this,gridinfo%griddim,categoryappend)
4455 CALL volgrid6d_alloc_vol(this, ini=.true.) ! decode?
4456
4457else if (.not. (this%griddim == gridinfo%griddim ))then
4458
4459 CALL l4f_category_log(this%category,l4f_error, &
4460 "volgrid and gridinfo grid type or size are different, gridinfo rejected")
4461 CALL raise_error()
4462 RETURN
4463
4464end if
4465
4466! Cerco gli indici del campo da inserire, se non trovo metto nel primo missing
4467ilevel = index(this%level, gridinfo%level)
4468IF (ilevel == 0 .AND. optio_log(force)) THEN
4469 ilevel = index(this%level, vol7d_level_miss)
4470 IF (ilevel /= 0) this%level(ilevel) = gridinfo%level
4471ENDIF
4472
4473IF (ilevel == 0) THEN
4474 CALL l4f_category_log(this%category,l4f_error, &
4475 "volgrid6d: level not valid for volume, gridinfo rejected")
4476 CALL raise_error()
4477 RETURN
4478ENDIF
4479
4480IF (optio_log(isanavar)) THEN ! assign to all times and timeranges
4481 itime0 = 1
4482 itime1 = SIZE(this%time)
4483 itimerange0 = 1
4484 itimerange1 = SIZE(this%timerange)
4485ELSE ! usual case
4486 correctedtime = gridinfo%time
4487 IF (this%time_definition == 1 .OR. this%time_definition == 2) correctedtime = correctedtime + &
4488 timedelta_new(sec=gridinfo%timerange%p1)
4489 itime0 = index(this%time, correctedtime)
4490 IF (itime0 == 0 .AND. optio_log(force)) THEN
4491 itime0 = index(this%time, datetime_miss)
4492 IF (itime0 /= 0) this%time(itime0) = correctedtime
4493 ENDIF
4494 IF (itime0 == 0) THEN
4495 CALL l4f_category_log(this%category,l4f_error, &
4496 "volgrid6d: time not valid for volume, gridinfo rejected")
4497 CALL raise_error()
4498 RETURN
4499 ENDIF
4500 itime1 = itime0
4501
4502 correctedtimerange = gridinfo%timerange
4503 IF (this%time_definition == 2) correctedtimerange%p1 = 0
4504 itimerange0 = index(this%timerange, correctedtimerange)
4505 IF (itimerange0 == 0 .AND. optio_log(force)) THEN
4506 itimerange0 = index(this%timerange, vol7d_timerange_miss)
4507 IF (itimerange0 /= 0) this%timerange(itimerange0) = gridinfo%timerange
4508 ENDIF
4509 IF (itimerange0 == 0) THEN
4510 CALL l4f_category_log(this%category,l4f_error, &
4511 "volgrid6d: timerange not valid for volume, gridinfo rejected")
4512 CALL raise_error()
4513 RETURN
4514 ENDIF
4515 itimerange1 = itimerange0
4516ENDIF
4517
4518ivar = index(this%var, gridinfo%var)
4519IF (ivar == 0 .AND. optio_log(force)) THEN
4520 ivar = index(this%var, volgrid6d_var_miss)
4521 IF (ivar /= 0) this%var(ivar) = gridinfo%var
4522ENDIF
4523IF (ivar == 0) THEN
4524 CALL l4f_category_log(this%category,l4f_error, &
4525 "volgrid6d: var not valid for volume, gridinfo rejected")
4526 CALL raise_error()
4527 RETURN
4528ENDIF
4529
4530DO itimerange = itimerange0, itimerange1
4531 DO itime = itime0, itime1
4532 IF (ASSOCIATED(this%gaid)) THEN
4533 dup = .false.
4534 IF (c_e(this%gaid(ilevel,itime,itimerange,ivar))) THEN
4535 dup = .true.
4536 CALL l4f_category_log(this%category,l4f_warn,"gaid exist: grib duplicated")
4537! avoid memory leaks
4538 IF (optio_log(clone)) CALL delete(this%gaid(ilevel,itime,itimerange,ivar))
4539 ENDIF
4540
4541 IF (optio_log(clone)) THEN
4542 CALL copy(gridinfo%gaid, this%gaid(ilevel,itime,itimerange,ivar))
4543#ifdef DEBUG
4544 CALL l4f_category_log(this%category,l4f_debug,"cloning to a new gaid")
4545#endif
4546 ELSE
4547 this%gaid(ilevel,itime,itimerange,ivar) = gridinfo%gaid
4548 ENDIF
4549
4550 IF (ASSOCIATED(this%voldati))THEN
4551 IF (.NOT.dup .OR. ldup_mode == 0) THEN
4552 this%voldati(:,:,ilevel,itime,itimerange,ivar) = decode_gridinfo(gridinfo)
4553 ELSE IF (ldup_mode == 1) THEN
4554 tmpgrid = decode_gridinfo(gridinfo) ! f2003 automatic allocation
4555 WHERE(c_e(tmpgrid))
4556 this%voldati(:,:,ilevel,itime,itimerange,ivar) = tmpgrid(:,:)
4557 END WHERE
4558 ELSE IF (ldup_mode == 2) THEN
4559 WHERE(.NOT.c_e(this%voldati(:,:,ilevel,itime,itimerange,ivar)))
4560 this%voldati(:,:,ilevel,itime,itimerange,ivar) = decode_gridinfo(gridinfo)
4561 END WHERE
4562 ENDIF
4563 ENDIF
4564
4565 ELSE
4566 CALL l4f_category_log(this%category,l4f_error, &
4567 "gaid not allocated, you probably need to call volgrid6d_alloc_vol first")
4568 CALL raise_error()
4569 RETURN
4570 ENDIF
4571 ENDDO
4572ENDDO
4573
4574
4575END SUBROUTINE import_from_gridinfo
4576
4577
4578!> Export a single grid of a \a volgrid6d object to a \a gridinfo_def object.
4579!! A single 2d slice of a \a volgrid6d at a specified location is
4580!! written into a \a gridinfo_def object, including the \a grid_id
4581!! which can be used for the successive export to file.
4582SUBROUTINE export_to_gridinfo(this, gridinfo, itime, itimerange, ilevel, ivar, &
4583 gaid_template, clone)
4584TYPE(volgrid6d),INTENT(in) :: this !< volume to be exported
4585TYPE(gridinfo_def),INTENT(inout) :: gridinfo !< output gridinfo_def object
4586INTEGER :: itime !< index within \a this of the element to be exported for the time dimension
4587INTEGER :: itimerange !< index within \a this of the element to be exported for the timerange dimension
4588INTEGER :: ilevel !< index within \a this of the element to be exported for the vertical level dimension
4589INTEGER :: ivar !< index within \a this of the element to be exported for the variable dimension
4590TYPE(grid_id),INTENT(in),OPTIONAL :: gaid_template !< \a grid_id template to be used for output data replacing the one contained in \a this
4591LOGICAL,INTENT(in),OPTIONAL :: clone !< if provided and \c .TRUE., clone the grid_id included in \a this rather than making a shallow copy
4592
4593TYPE(grid_id) :: gaid
4594LOGICAL :: usetemplate
4595REAL,POINTER :: voldati(:,:)
4596TYPE(datetime) :: correctedtime
4597
4598#ifdef DEBUG
4599CALL l4f_category_log(this%category,l4f_debug,"export_to_gridinfo")
4600#endif
4601
4602IF (.NOT.c_e(this%gaid(ilevel,itime,itimerange,ivar))) THEN
4603#ifdef DEBUG
4604 CALL l4f_category_log(this%category,l4f_debug,"empty gaid found, skipping export")
4605#endif
4606 RETURN
4607ENDIF
4608
4609usetemplate = .false.
4610IF (PRESENT(gaid_template)) THEN
4611 CALL copy(gaid_template, gaid)
4612#ifdef DEBUG
4613 CALL l4f_category_log(this%category,l4f_debug,"template cloned to a new gaid")
4614#endif
4615 usetemplate = c_e(gaid)
4616ENDIF
4617
4618IF (.NOT.usetemplate) THEN
4619 IF (optio_log(clone)) THEN
4620 CALL copy(this%gaid(ilevel,itime,itimerange,ivar), gaid)
4621#ifdef DEBUG
4622 CALL l4f_category_log(this%category,l4f_debug,"original gaid cloned to a new one")
4623#endif
4624 ELSE
4625 gaid = this%gaid(ilevel,itime,itimerange,ivar)
4626 ENDIF
4627ENDIF
4628
4629IF (this%time_definition == 1 .OR. this%time_definition == 2) THEN
4630 correctedtime = this%time(itime) - &
4631 timedelta_new(sec=this%timerange(itimerange)%p1)
4632ELSE
4633 correctedtime = this%time(itime)
4634ENDIF
4635
4636CALL init(gridinfo,gaid, this%griddim, correctedtime, this%timerange(itimerange), &
4637 this%level(ilevel), this%var(ivar))
4638
4639! reset the gridinfo, bad but necessary at this point for encoding the field
4640CALL export(gridinfo%griddim, gridinfo%gaid)
4641! encode the field
4642IF (ASSOCIATED(this%voldati)) THEN
4643 CALL encode_gridinfo(gridinfo, this%voldati(:,:,ilevel,itime,itimerange,ivar))
4644ELSE IF (usetemplate) THEN ! field must be forced into template in this case
4645 NULLIFY(voldati)
4646 CALL volgrid_get_vol_2d(this, ilevel, itime, itimerange, ivar, voldati)
4647 CALL encode_gridinfo(gridinfo, voldati)
4648 DEALLOCATE(voldati)
4649ENDIF
4650
4651END SUBROUTINE export_to_gridinfo
4652
4653
4654!> Import an array of \a gridinfo objects into an array of \a
4655!! volgrid6d objects.
4656!! This method imports an array of gridded fields from an \a
4657!! arrayof_gridinfo object into a suitable number of \a volgrid6d
4658!! objects. The number of \a volgrid6d allocated is determined by the
4659!! number of different grids encountered in \a arrayof_gridinfo.
4660!! Unlike the import for a single \a gridinfo, here the \a volgrid6d
4661!! object is a non-associated pointer to a 1-d array of uninitialized
4662!! objects, and all the dimension descriptors in each of the objects
4663!! are allocated and assigned within the method according to the data
4664!! contained in \a gridinfov.
4665!! If the \a anavar array argument is provided, all the input messages
4666!! whose variable maps to one of the B-table variables contained in \a
4667!! anavar are treated as time-independent (AKA anagraphic data,
4668!! station data, etc.), thus their time and timerange are ignored and
4669!! they are replicated for every time and timerange present in the
4670!! corresponding data volume.
4671SUBROUTINE import_from_gridinfovv(this, gridinfov, dup_mode, clone, decode, &
4672 time_definition, anavar, categoryappend)
4673TYPE(volgrid6d),POINTER :: this(:) !< object in which to import
4674TYPE(arrayof_gridinfo),INTENT(in) :: gridinfov !< array of gridinfo objects to be imported
4675INTEGER,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
4676LOGICAL , INTENT(in),OPTIONAL :: clone !< if provided and \c .TRUE. , clone the gaid's from \a gridinfo to \a this
4677LOGICAL,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
4678INTEGER,INTENT(IN),OPTIONAL :: time_definition !< 0=time is reference time; 1=time is validity time
4679CHARACTER(len=*),INTENT(IN),OPTIONAL :: anavar(:) !< list of variables (B-table code equivalent) to be treated as time-independent data
4680CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
4681
4682INTEGER :: i, j, stallo
4683INTEGER :: ngrid, ntime, ntimerange, nlevel, nvar, ltime_definition
4684INTEGER :: category
4685CHARACTER(len=512) :: a_name
4686TYPE(datetime),ALLOCATABLE :: correctedtime(:)
4687LOGICAL,ALLOCATABLE :: isanavar(:)
4688TYPE(vol7d_var) :: lvar
4689TYPE(vol7d_timerange),ALLOCATABLE :: correctedtimerange(:)
4690
4691! category temporanea (altrimenti non possiamo loggare)
4692if (present(categoryappend))then
4693 call l4f_launcher(a_name,a_name_append=trim(subcategory)//"."//trim(categoryappend))
4694else
4695 call l4f_launcher(a_name,a_name_append=trim(subcategory))
4696endif
4697category=l4f_category_get(a_name)
4698
4699#ifdef DEBUG
4700call l4f_category_log(category,l4f_debug,"start import_from_gridinfovv")
4701#endif
4702
4703IF (PRESENT(time_definition)) THEN
4704 ltime_definition = max(min(time_definition, 2), 0)
4705ELSE
4706 ltime_definition = 0
4707ENDIF
4708
4709ngrid=count_distinct(gridinfov%array(1:gridinfov%arraysize)%griddim,back=.true.)
4710CALL l4f_category_log(category,l4f_info, t2c(ngrid)// &
4711 ' different grid definition(s) found in input data')
4712
4713ALLOCATE(this(ngrid),stat=stallo)
4714IF (stallo /= 0)THEN
4715 CALL l4f_category_log(category,l4f_fatal,"allocating memory")
4716 CALL raise_fatal_error()
4717ENDIF
4718DO i = 1, ngrid
4719 IF (PRESENT(categoryappend))THEN
4720 CALL init(this(i), time_definition=ltime_definition, categoryappend=trim(categoryappend)//"-vol"//t2c(i))
4721 ELSE
4722 CALL init(this(i), time_definition=ltime_definition, categoryappend="vol"//t2c(i))
4723 ENDIF
4724ENDDO
4725
4726this(:)%griddim=pack_distinct(gridinfov%array(1:gridinfov%arraysize)%griddim, &
4727 ngrid, back=.true.)
4728
4729! mark elements as ana variables (time-independent)
4730ALLOCATE(isanavar(gridinfov%arraysize))
4731isanavar(:) = .false.
4732IF (PRESENT(anavar)) THEN
4733 DO i = 1, gridinfov%arraysize
4734 DO j = 1, SIZE(anavar)
4735 lvar = convert(gridinfov%array(i)%var)
4736 IF (lvar%btable == anavar(j)) THEN
4737 isanavar(i) = .true.
4738 EXIT
4739 ENDIF
4740 ENDDO
4741 ENDDO
4742 CALL l4f_category_log(category,l4f_info,t2c(count(isanavar))//'/'// &
4743 t2c(gridinfov%arraysize)//' constant-data messages found in input data')
4744ENDIF
4745
4746IF (ltime_definition == 1 .OR. ltime_definition == 2) THEN ! verification time
4747 ALLOCATE(correctedtime(gridinfov%arraysize))
4748 correctedtime(:) = gridinfov%array(1:gridinfov%arraysize)%time
4749 DO i = 1, gridinfov%arraysize
4750 correctedtime(i) = correctedtime(i) + &
4751 timedelta_new(sec=gridinfov%array(i)%timerange%p1)
4752 ENDDO
4753ENDIF
4754IF (ltime_definition == 2) THEN ! set all to analysis
4755 ALLOCATE(correctedtimerange(gridinfov%arraysize))
4756 correctedtimerange(:) = gridinfov%array(1:gridinfov%arraysize)%timerange
4757 correctedtimerange(:)%p1 = 0
4758ENDIF
4759
4760DO i = 1, ngrid
4761 IF (PRESENT(anavar)) THEN
4762 j = count((this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
4763 .AND. .NOT.isanavar(:))
4764 IF (j <= 0) THEN
4765 CALL l4f_category_log(category, l4f_fatal, 'grid n.'//t2c(i)// &
4766 ' has only constant data, this is not allowed')
4767 CALL l4f_category_log(category, l4f_fatal, 'please check anavar argument')
4768 CALL raise_fatal_error()
4769 ENDIF
4770 ENDIF
4771 IF (ltime_definition == 1 .OR. ltime_definition == 2) THEN ! verification time
4772 ntime = count_distinct(correctedtime, &
4773 mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
4774 .AND. .NOT.isanavar(:), back=.true.)
4775 ELSE
4776 ntime = count_distinct(gridinfov%array(1:gridinfov%arraysize)%time, &
4777 mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
4778 .AND. .NOT.isanavar(:), back=.true.)
4779 ENDIF
4780 IF (ltime_definition == 2) THEN ! set all to analysis
4781 ntimerange = count_distinct(correctedtimerange, &
4782 mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
4783 .AND. .NOT.isanavar(:), back=.true.)
4784 ELSE
4785 ntimerange = count_distinct(gridinfov%array(1:gridinfov%arraysize)%timerange, &
4786 mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
4787 .AND. .NOT.isanavar(:), back=.true.)
4788 ENDIF
4789 nlevel = count_distinct(gridinfov%array(1:gridinfov%arraysize)%level, &
4790 mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim), &
4791 back=.true.)
4792 nvar = count_distinct(gridinfov%array(1:gridinfov%arraysize)%var, &
4793 mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim), &
4794 back=.true.)
4795
4796#ifdef DEBUG
4797 CALL l4f_category_log(this(i)%category,l4f_debug,"alloc volgrid6d index: "//t2c(i))
4798#endif
4799
4800 CALL volgrid6d_alloc(this(i),this(i)%griddim%dim,ntime=ntime, &
4801 ntimerange=ntimerange,nlevel=nlevel,nvar=nvar)
4802
4803 IF (ltime_definition == 1 .OR. ltime_definition == 2) THEN ! verification time
4804 this(i)%time = pack_distinct(correctedtime, ntime, &
4805 mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
4806 .AND. .NOT.isanavar(:), back=.true.)
4807 ELSE
4808 this(i)%time = pack_distinct(gridinfov%array(1:gridinfov%arraysize)%time, ntime, &
4809 mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
4810 .AND. .NOT.isanavar(:), back=.true.)
4811 ENDIF
4812 CALL sort(this(i)%time)
4813
4814 IF (ltime_definition == 2) THEN ! set all to analysis
4815 this(i)%timerange = pack_distinct(correctedtimerange, ntimerange, &
4816 mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
4817 .AND. .NOT.isanavar(:), back=.true.)
4818 ELSE
4819 this(i)%timerange = pack_distinct(gridinfov%array(1:gridinfov%arraysize)%timerange, &
4820 ntimerange, mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
4821 .AND. .NOT.isanavar(:), back=.true.)
4822 ENDIF
4823 CALL sort(this(i)%timerange)
4824
4825 this(i)%level=pack_distinct(gridinfov%array(1:gridinfov%arraysize)%level, &
4826 nlevel,mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim), &
4827 back=.true.)
4828 CALL sort(this(i)%level)
4829
4830 this(i)%var=pack_distinct(gridinfov%array(1:gridinfov%arraysize)%var, nvar, &
4831 mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim), &
4832 back=.true.)
4833
4834#ifdef DEBUG
4835 CALL l4f_category_log(this(i)%category,l4f_debug,"alloc_vol volgrid6d index: "//t2c(i))
4836#endif
4837 CALL volgrid6d_alloc_vol(this(i), decode=decode)
4838
4839ENDDO
4840
4841IF (ltime_definition == 1 .OR. ltime_definition == 2) DEALLOCATE(correctedtime)
4842IF (ltime_definition == 2) DEALLOCATE(correctedtimerange)
4843
4844DO i = 1, gridinfov%arraysize
4845
4846#ifdef DEBUG
4847 CALL l4f_category_log(category,l4f_debug,"import from gridinfov index: "//t2c(i))
4848 CALL l4f_category_log(category,l4f_info, &
4849 "to volgrid6d index: "//t2c(index(this%griddim, gridinfov%array(i)%griddim)))
4850#endif
4851
4852 CALL import(this(index(this%griddim, gridinfov%array(i)%griddim)), &
4853 gridinfov%array(i), dup_mode=dup_mode, clone=clone, isanavar=isanavar(i))
4854
4855ENDDO
4856
4857!chiudo il logger temporaneo
4858CALL l4f_category_delete(category)
4859
4860END SUBROUTINE import_from_gridinfovv
4861
4862
4863!> Export a \a volgrid6d object to an \a arrayof_gridinfo object.
4864!! The multidimensional \a volgrid6d structure is serialized into a
4865!! one-dimensional array of gridinfo_def objects, which is allocated
4866!! to the proper size if not already allocated, or it is extended
4867!! keeping the old data if any.
4868SUBROUTINE export_to_gridinfov(this, gridinfov, gaid_template, clone)
4869TYPE(volgrid6d),INTENT(inout) :: this !< volume to be exported
4870TYPE(arrayof_gridinfo),INTENT(inout) :: gridinfov !< output array of gridinfo_def objects
4871TYPE(grid_id),INTENT(in),OPTIONAL :: gaid_template !< \a grid_id template to be used for output data replacing the one contained in \a this
4872LOGICAL,INTENT(in),OPTIONAL :: clone !< if provided and \c .TRUE., clone the grid_id included in \a this rather than making a shallow copy
4873
4874INTEGER :: i ,itime, itimerange, ilevel, ivar
4875INTEGER :: ntime, ntimerange, nlevel, nvar
4876TYPE(gridinfo_def) :: gridinfol
4877
4878#ifdef DEBUG
4879CALL l4f_category_log(this%category,l4f_debug,"start export_to_gridinfov")
4880#endif
4881
4882! this is necessary in order not to repeatedly and uselessly copy the
4883! same array of coordinates for each 2d grid during export, the
4884! side-effect is that the computed projection in this is lost
4885CALL dealloc(this%griddim%dim)
4886
4887i=0
4888ntime=size(this%time)
4889ntimerange=size(this%timerange)
4890nlevel=size(this%level)
4891nvar=size(this%var)
4892
4893DO itime=1,ntime
4894 DO itimerange=1,ntimerange
4895 DO ilevel=1,nlevel
4896 DO ivar=1,nvar
4897
4898 CALL init(gridinfol)
4899 CALL export(this, gridinfol, itime, itimerange, ilevel, ivar, &
4900 gaid_template=gaid_template, clone=clone)
4901 IF (c_e(gridinfol%gaid)) THEN
4902 CALL insert(gridinfov, gridinfol)
4903 ELSE
4904 CALL delete(gridinfol)
4905 ENDIF
4906
4907 ENDDO
4908 ENDDO
4909 ENDDO
4910ENDDO
4911
4912END SUBROUTINE export_to_gridinfov
4913
4914
4915!> Export an array of \a volgrid6d objects to an \a arrayof_gridinfo object.
4916!! The multidimensional \a volgrid6d structures are serialized into a
4917!! one-dimensional array of gridinfo_def objects, which is allocated
4918!! to the proper size if not already allocated, or it is extended
4919!! keeping the old data if any.
4920SUBROUTINE export_to_gridinfovv(this, gridinfov, gaid_template, clone)
4921!, &
4922! categoryappend)
4923TYPE(volgrid6d),INTENT(inout) :: this(:) !< volume array to be exported
4924TYPE(arrayof_gridinfo),INTENT(inout) :: gridinfov !< output array of gridinfo_def objects
4925TYPE(grid_id),INTENT(in),OPTIONAL :: gaid_template !< \a grid_id template to be used for output data replacing the one contained in \a this
4926LOGICAL,INTENT(in),OPTIONAL :: clone !< if provided and \c .TRUE., clone the grid_id included in \a this rather than making a shallow copy
4927
4928INTEGER :: i
4929
4930DO i = 1, SIZE(this)
4931#ifdef DEBUG
4932 CALL l4f_category_log(this(i)%category,l4f_debug, &
4933 "export_to_gridinfovv grid index: "//t2c(i))
4934#endif
4935 CALL export(this(i), gridinfov, gaid_template=gaid_template, clone=clone)
4936ENDDO
4937
4938END SUBROUTINE export_to_gridinfovv
4939
4940
4941!> Import the content of a supported file (like grib or gdal-supported
4942!! format) into an array of \a volgrid6d objects. This method imports
4943!! a set of gridded fields from a file object into a suitable number
4944!! of \a volgrid6d objects. The data are imported by creating a
4945!! temporary \a gridinfo object, importing it into the \a volgrid6d
4946!! object cloning the gaid's and then destroying the gridinfo, so it
4947!! works similarly to volgrid6d_class::import_from_gridinfovv() method.
4948!! For a detailed explanation of the \a anavar argument, see the
4949!! documentation of volgrid6d_class::import_from_gridinfovv() method.
4950SUBROUTINE volgrid6d_import_from_file(this, filename, dup_mode, decode, &
4951 time_definition, anavar, categoryappend)
4952TYPE(volgrid6d),POINTER :: this(:) !< object in which to import
4953CHARACTER(len=*),INTENT(in) :: filename !< name of file from which to import
4954INTEGER,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
4955LOGICAL,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
4956INTEGER,INTENT(IN),OPTIONAL :: time_definition !< 0=time is reference time; 1=time is validity time
4957CHARACTER(len=*),INTENT(IN),OPTIONAL :: anavar(:) !< list of variables (B-table code equivalent) to be treated as time-independent data
4958character(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
4959
4960TYPE(arrayof_gridinfo) :: gridinfo
4961INTEGER :: category
4962CHARACTER(len=512) :: a_name
4963
4964NULLIFY(this)
4965
4966IF (PRESENT(categoryappend))THEN
4967 CALL l4f_launcher(a_name,a_name_append= &
4968 trim(subcategory)//"."//trim(categoryappend))
4969ELSE
4970 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
4971ENDIF
4972category=l4f_category_get(a_name)
4973
4974CALL import(gridinfo, filename=filename, categoryappend=categoryappend)
4975
4976IF (gridinfo%arraysize > 0) THEN
4977
4978 CALL import(this, gridinfo, dup_mode=dup_mode, clone=.true., decode=decode, &
4979 time_definition=time_definition, anavar=anavar, &
4980 categoryappend=categoryappend)
4981
4982 CALL l4f_category_log(category,l4f_info,"deleting gridinfo")
4983 CALL delete(gridinfo)
4984
4985ELSE
4986 CALL l4f_category_log(category,l4f_info,"file does not contain gridded data")
4987ENDIF
4988
4989! close logger
4990CALL l4f_category_delete(category)
4991
4992END SUBROUTINE volgrid6d_import_from_file
4993
4994
4995!> High level method for exporting a volume array to file.
4996!! All the information contained into an array of \a volgrid6d
4997!! objects, i.e. dimension descriptors and data, is exported to a file
4998!! using the proper output driver (typically grib_api for grib
4999!! format). If a template is provided, it will determine the
5000!! characteristic of the output file, otherwise the \a grid_id
5001!! descriptors contained in the volgrid6d object will be used
5002SUBROUTINE volgrid6d_export_to_file(this, filename, gaid_template, categoryappend)
5003TYPE(volgrid6d) :: this(:) !< volume(s) to be exported
5004CHARACTER(len=*),INTENT(in) :: filename !< output file name
5005TYPE(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
5006character(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
5007
5008TYPE(arrayof_gridinfo) :: gridinfo
5009INTEGER :: category
5010CHARACTER(len=512) :: a_name
5011
5012IF (PRESENT(categoryappend)) THEN
5013 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//"."//trim(categoryappend))
5014ELSE
5015 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
5016ENDIF
5017category=l4f_category_get(a_name)
5018
5019#ifdef DEBUG
5020CALL l4f_category_log(category,l4f_debug,"start export to file")
5021#endif
5022
5023CALL l4f_category_log(category,l4f_info,"writing volgrid6d to grib file: "//trim(filename))
5024
5025!IF (ASSOCIATED(this)) THEN
5026 CALL export(this, gridinfo, gaid_template=gaid_template, clone=.true.)
5027 IF (gridinfo%arraysize > 0) THEN
5028 CALL export(gridinfo, filename)
5029 CALL delete(gridinfo)
5030 ENDIF
5031!ELSE
5032! CALL l4f_category_log(category,L4F_INFO,"volume volgrid6d is not associated")
5033!ENDIF
5034
5035! close logger
5036CALL l4f_category_delete(category)
5037
5038END SUBROUTINE volgrid6d_export_to_file
5039
5040
5041!> Array destructor for \a volgrid6d class.
5042!! Delete an array of \a volgrid6d objects and deallocate the array
5043!! itself.
5044SUBROUTINE volgrid6dv_delete(this)
5045TYPE(volgrid6d),POINTER :: this(:) !< vector of volgrid6d object
5046
5047INTEGER :: i
5048
5049IF (ASSOCIATED(this)) THEN
5050 DO i = 1, SIZE(this)
5051#ifdef DEBUG
5052 CALL l4f_category_log(this(i)%category,l4f_debug, &
5053 "delete volgrid6d vector index: "//trim(to_char(i)))
5054#endif
5055 CALL delete(this(i))
5056 ENDDO
5057 DEALLOCATE(this)
5058ENDIF
5059
5060END SUBROUTINE volgrid6dv_delete
5061
5062
5063! Internal method for performing grid to grid computations
5064SUBROUTINE volgrid6d_transform_compute(this, volgrid6d_in, volgrid6d_out, &
5065 lev_out, var_coord_vol, clone)
5066TYPE(grid_transform),INTENT(in) :: this ! oggetto di trasformazione per il grigliato
5067type(volgrid6d), INTENT(in) :: volgrid6d_in ! oggetto da trasformare
5068type(volgrid6d), INTENT(inout) :: volgrid6d_out ! oggetto trasformato; deve essere completo (init, alloc, alloc_vol)
5069TYPE(vol7d_level),INTENT(in),OPTIONAL :: lev_out(:) ! vol7d_level object defining target vertical grid, for vertical interpolations
5070INTEGER,INTENT(in),OPTIONAL :: var_coord_vol ! index of variable defining vertical coordinate values in input volume
5071LOGICAL,INTENT(in),OPTIONAL :: clone ! se fornito e \c .TRUE., clona i gaid da volgrid6d_in a volgrid6d_out
5072
5073INTEGER :: ntime, ntimerange, inlevel, onlevel, nvar, &
5074 itime, itimerange, ilevel, ivar, levshift, levused, lvar_coord_vol, spos
5075REAL,POINTER :: voldatiin(:,:,:), voldatiout(:,:,:), coord_3d_in(:,:,:)
5076TYPE(vol7d_level) :: output_levtype
5077
5078
5079#ifdef DEBUG
5080call l4f_category_log(volgrid6d_in%category,l4f_debug,"start volgrid6d_transform_compute")
5081#endif
5082
5083ntime=0
5084ntimerange=0
5085inlevel=0
5086onlevel=0
5087nvar=0
5088lvar_coord_vol = optio_i(var_coord_vol)
5089
5090if (associated(volgrid6d_in%time))then
5091 ntime=size(volgrid6d_in%time)
5092 volgrid6d_out%time=volgrid6d_in%time
5093end if
5094
5095if (associated(volgrid6d_in%timerange))then
5096 ntimerange=size(volgrid6d_in%timerange)
5097 volgrid6d_out%timerange=volgrid6d_in%timerange
5098end if
5099
5100IF (ASSOCIATED(volgrid6d_in%level))THEN
5101 inlevel=SIZE(volgrid6d_in%level)
5102ENDIF
5103IF (PRESENT(lev_out)) THEN
5104 onlevel=SIZE(lev_out)
5105 volgrid6d_out%level=lev_out
5106ELSE IF (ASSOCIATED(volgrid6d_in%level))THEN
5107 onlevel=SIZE(volgrid6d_in%level)
5108 volgrid6d_out%level=volgrid6d_in%level
5109ENDIF
5110
5111if (associated(volgrid6d_in%var))then
5112 nvar=size(volgrid6d_in%var)
5113 volgrid6d_out%var=volgrid6d_in%var
5114end if
5115! allocate once for speed
5116IF (.NOT.ASSOCIATED(volgrid6d_in%voldati)) THEN
5117 ALLOCATE(voldatiin(volgrid6d_in%griddim%dim%nx, volgrid6d_in%griddim%dim%ny, &
5118 inlevel))
5119ENDIF
5120IF (.NOT.ASSOCIATED(volgrid6d_out%voldati)) THEN
5121 ALLOCATE(voldatiout(volgrid6d_out%griddim%dim%nx, volgrid6d_out%griddim%dim%ny, &
5122 onlevel))
5123ENDIF
5124
5125CALL get_val(this, levshift=levshift, levused=levused)
5126spos = imiss
5127IF (c_e(lvar_coord_vol)) THEN
5128 CALL get_val(this%trans, output_levtype=output_levtype)
5129 IF (output_levtype%level1 == 103 .OR. output_levtype%level1 == 108) THEN
5130 spos = firsttrue(volgrid6d_in%level(:) == vol7d_level_new(1))
5131 IF (spos == 0) THEN
5132 CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
5133 'output level '//t2c(output_levtype%level1)// &
5134 ' requested, but height/press of surface not provided in volume')
5135 ENDIF
5136 IF (.NOT.c_e(levshift) .AND. .NOT.c_e(levused)) THEN
5137 CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
5138 'internal inconsistence, levshift and levused undefined when they should be')
5139 ENDIF
5140 ENDIF
5141ENDIF
5142
5143DO ivar=1,nvar
5144! IF (c_e(var_coord_vol)) THEN
5145! IF (ivar == var_coord_vol) CYCLE ! skip coordinate variable in output
5146! ENDIF
5147 DO itimerange=1,ntimerange
5148 DO itime=1,ntime
5149! skip empty columns where possible, improve
5150 IF (c_e(levshift) .AND. c_e(levused)) THEN
5151 IF (.NOT.any(c_e( &
5152 volgrid6d_in%gaid(levshift+1:levshift+levused,itime,itimerange,ivar) &
5153 ))) cycle
5154 ENDIF
5155 DO ilevel = 1, min(inlevel,onlevel)
5156! if present gaid copy it
5157 IF (c_e(volgrid6d_in%gaid(ilevel,itime,itimerange,ivar)) .AND. .NOT. &
5158 c_e(volgrid6d_out%gaid(ilevel,itime,itimerange,ivar))) THEN
5159
5160 IF (optio_log(clone)) THEN
5161 CALL copy(volgrid6d_in%gaid(ilevel,itime,itimerange,ivar),&
5162 volgrid6d_out%gaid(ilevel,itime,itimerange,ivar))
5163#ifdef DEBUG
5164 CALL l4f_category_log(volgrid6d_in%category,l4f_debug, &
5165 "cloning gaid, level "//t2c(ilevel))
5166#endif
5167 ELSE
5168 volgrid6d_out%gaid(ilevel,itime,itimerange,ivar) = &
5169 volgrid6d_in%gaid(ilevel,itime,itimerange,ivar)
5170 ENDIF
5171 ENDIF
5172 ENDDO
5173! if out level are more, we have to clone anyway
5174 DO ilevel = min(inlevel,onlevel) + 1, onlevel
5175 IF (c_e(volgrid6d_in%gaid(inlevel,itime,itimerange,ivar)) .AND. .NOT. &
5176 c_e(volgrid6d_out%gaid(ilevel,itime,itimerange,ivar))) then
5177
5178 CALL copy(volgrid6d_in%gaid(inlevel,itime,itimerange,ivar),&
5179 volgrid6d_out%gaid(ilevel,itime,itimerange,ivar))
5180#ifdef DEBUG
5181 CALL l4f_category_log(volgrid6d_in%category,l4f_debug, &
5182 "forced cloning gaid, level "//t2c(inlevel)//"->"//t2c(ilevel))
5183#endif
5184 ENDIF
5185 ENDDO
5186
5187 IF (c_e(lvar_coord_vol)) THEN
5188 NULLIFY(coord_3d_in)
5189 CALL volgrid_get_vol_3d(volgrid6d_in, itime, itimerange, lvar_coord_vol, &
5190 coord_3d_in)
5191 IF (c_e(spos)) THEN ! compute difference wrt surface coordinate
5192 IF (spos == 0) THEN ! error condition, set all to missing and goodnight
5193 coord_3d_in(:,:,levshift+1:levshift+levused) = rmiss
5194 ELSE
5195 DO ilevel = levshift+1, levshift+levused
5196 WHERE(c_e(coord_3d_in(:,:,ilevel)) .AND. c_e(coord_3d_in(:,:,spos)))
5197 coord_3d_in(:,:,ilevel) = coord_3d_in(:,:,ilevel) - &
5198 coord_3d_in(:,:,spos)
5199 ELSEWHERE
5200 coord_3d_in(:,:,ilevel) = rmiss
5201 END WHERE
5202 ENDDO
5203 ENDIF
5204 ENDIF
5205 ENDIF
5206 CALL volgrid_get_vol_3d(volgrid6d_in, itime, itimerange, ivar, &
5207 voldatiin)
5208 IF (ASSOCIATED(volgrid6d_out%voldati)) & ! improve!!!!
5209 CALL volgrid_get_vol_3d(volgrid6d_out, itime, itimerange, ivar, &
5210 voldatiout)
5211 IF (c_e(lvar_coord_vol)) THEN
5212 CALL compute(this, voldatiin, voldatiout, convert(volgrid6d_in%var(ivar)), &
5213 coord_3d_in(:,:,levshift+1:levshift+levused)) ! subset coord_3d_in
5214 ELSE
5215 CALL compute(this, voldatiin, voldatiout, convert(volgrid6d_in%var(ivar)))
5216 ENDIF
5217 CALL volgrid_set_vol_3d(volgrid6d_out, itime, itimerange, ivar, &
5218 voldatiout)
5219 ENDDO
5220 ENDDO
5221ENDDO
5222
5223IF (c_e(lvar_coord_vol)) THEN
5224 DEALLOCATE(coord_3d_in)
5225ENDIF
5226IF (.NOT.ASSOCIATED(volgrid6d_in%voldati)) THEN
5227 DEALLOCATE(voldatiin)
5228ENDIF
5229IF (.NOT.ASSOCIATED(volgrid6d_out%voldati)) THEN
5230 DEALLOCATE(voldatiout)
5231ENDIF
5232
5233
5234END SUBROUTINE volgrid6d_transform_compute
5235
5236
5237!> Performs the specified abstract transformation on the data provided.
5238!! The abstract transformation is specified by \a this parameter; the
5239!! corresponding specifical transformation (\a grid_transform object)
5240!! is created and destroyed internally. The output transformed object
5241!! is created internally and it does not require preliminary
5242!! initialisation.
5243SUBROUTINE volgrid6d_transform(this, griddim, volgrid6d_in, volgrid6d_out, &
5244 lev_out, volgrid6d_coord_in, maskgrid, maskbounds, clone, decode, categoryappend)
5245TYPE(transform_def),INTENT(in) :: this !< object specifying the abstract transformation
5246TYPE(griddim_def),INTENT(in),OPTIONAL :: griddim !< griddim specifying the output grid (required by most transformation types)
5247TYPE(volgrid6d),INTENT(inout) :: volgrid6d_in !< object to be transformed, it is not modified, despite the INTENT(inout)
5248TYPE(volgrid6d),INTENT(out) :: volgrid6d_out !< transformed object, it does not require initialisation
5249TYPE(vol7d_level),INTENT(in),OPTIONAL,TARGET :: lev_out(:) !< vol7d_level object defining target vertical grid, for vertical interpolations
5250TYPE(volgrid6d),INTENT(in),OPTIONAL :: volgrid6d_coord_in !< object providing time constant input vertical coordinate for some kind of vertical interpolations
5251REAL,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')
5252REAL,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')
5253LOGICAL,INTENT(in),OPTIONAL :: clone !< if provided and \a .TRUE. , clone the \a gaid's from \a volgrid6d_in to \a volgrid6d_out
5254LOGICAL,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
5255CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
5256
5257TYPE(grid_transform) :: grid_trans
5258TYPE(vol7d_level),POINTER :: llev_out(:)
5259TYPE(vol7d_level) :: input_levtype, output_levtype
5260TYPE(vol7d_var) :: vcoord_var
5261INTEGER :: i, k, ntime, ntimerange, nlevel, nvar, var_coord_in, var_coord_vol, &
5262 cf_out, nxc, nyc, nxi, nyi, i3, i4, i5, i6, &
5263 ulstart, ulend, spos
5264REAL,ALLOCATABLE :: coord_3d_in(:,:,:)
5265TYPE(geo_proj) :: proj_in, proj_out
5266CHARACTER(len=80) :: trans_type
5267LOGICAL :: ldecode
5268LOGICAL,ALLOCATABLE :: mask_in(:)
5269
5270#ifdef DEBUG
5271call l4f_category_log(volgrid6d_in%category, l4f_debug, "start volgrid6d_transform")
5272#endif
5273
5274ntime=0
5275ntimerange=0
5276nlevel=0
5277nvar=0
5278
5279if (associated(volgrid6d_in%time)) ntime=size(volgrid6d_in%time)
5280if (associated(volgrid6d_in%timerange)) ntimerange=size(volgrid6d_in%timerange)
5281if (associated(volgrid6d_in%level)) nlevel=size(volgrid6d_in%level)
5282if (associated(volgrid6d_in%var)) nvar=size(volgrid6d_in%var)
5283
5284IF (ntime == 0 .OR. ntimerange == 0 .OR. nlevel == 0 .OR. nvar == 0) THEN
5285 CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
5286 "trying to transform an incomplete volgrid6d object, ntime="//t2c(ntime)// &
5287 ' ntimerange='//t2c(ntimerange)//' nlevel='//t2c(nlevel)//' nvar='//t2c(nvar))
5288 CALL init(volgrid6d_out) ! initialize to empty
5289 CALL raise_error()
5290 RETURN
5291ENDIF
5292
5293CALL get_val(this, trans_type=trans_type)
5294
5295! store desired output component flag and unrotate if necessary
5296cf_out = imiss
5297IF (PRESENT(griddim) .AND. (trans_type == 'inter' .OR. trans_type == 'boxinter' &
5298 .OR. trans_type == 'stencilinter')) THEN ! improve condition!!
5299 CALL get_val(volgrid6d_in%griddim, proj=proj_in)
5300 CALL get_val(griddim, component_flag=cf_out, proj=proj_out)
5301! if different projections wind components must be referred to geographical system
5302 IF (proj_in /= proj_out) CALL vg6d_wind_unrot(volgrid6d_in)
5303ELSE IF (PRESENT(griddim)) THEN ! just get component_flag, the rest is rubbish
5304 CALL get_val(griddim, component_flag=cf_out)
5305ENDIF
5306
5307
5308var_coord_in = imiss
5309var_coord_vol = imiss
5310IF (trans_type == 'vertint') THEN
5311 IF (PRESENT(lev_out)) THEN
5312
5313! if volgrid6d_coord_in provided and allocated, check that it fits
5314 IF (PRESENT(volgrid6d_coord_in)) THEN
5315 IF (ASSOCIATED(volgrid6d_coord_in%voldati)) THEN
5316
5317! strictly 1 time and 1 timerange
5318 IF (SIZE(volgrid6d_coord_in%voldati,4) /= 1 .OR. &
5319 SIZE(volgrid6d_coord_in%voldati,5) /= 1) THEN
5320 CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
5321 'volume providing constant input vertical coordinate must have &
5322 &only 1 time and 1 timerange')
5323 CALL init(volgrid6d_out) ! initialize to empty
5324 CALL raise_error()
5325 RETURN
5326 ENDIF
5327
5328! search for variable providing vertical coordinate
5329 CALL get_val(this, output_levtype=output_levtype)
5330 vcoord_var = vol7d_var_new(vol7d_level_to_var(output_levtype))
5331 IF (.NOT.c_e(vcoord_var)) THEN
5332 CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
5333 'requested output level type '//t2c(output_levtype%level1)// &
5334 ' does not correspond to any known physical variable for &
5335 &providing vertical coordinate')
5336 CALL init(volgrid6d_out) ! initialize to empty
5337 CALL raise_error()
5338 RETURN
5339 ENDIF
5340
5341 DO i = 1, SIZE(volgrid6d_coord_in%var)
5342 IF (convert(volgrid6d_coord_in%var(i)) == vcoord_var) THEN
5343 var_coord_in = i
5344 EXIT
5345 ENDIF
5346 ENDDO
5347
5348 IF (.NOT.c_e(var_coord_in)) THEN
5349 CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
5350 'volume providing constant input vertical coordinate contains no &
5351 &variables matching output level type '//t2c(output_levtype%level1))
5352 CALL init(volgrid6d_out) ! initialize to empty
5353 CALL raise_error()
5354 RETURN
5355 ENDIF
5356 CALL l4f_category_log(volgrid6d_in%category, l4f_info, &
5357 'Coordinate for vertint found in coord volume at position '// &
5358 t2c(var_coord_in))
5359
5360! check horizontal grid
5361! this is too strict (component flag and so on)
5362! IF (volgrid6d_coord_in%griddim /= volgrid6d_in%griddim) THEN
5363 CALL get_val(volgrid6d_coord_in%griddim, nx=nxc, ny=nyc)
5364 CALL get_val(volgrid6d_in%griddim, nx=nxi, ny=nyi)
5365 IF (nxc /= nxi .OR. nyc /= nyi) THEN
5366 CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
5367 'volume providing constant input vertical coordinate must have &
5368 &the same grid as the input')
5369 CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
5370 'vertical coordinate: '//t2c(nxc)//'x'//t2c(nyc)// &
5371 ', input volume: '//t2c(nxi)//'x'//t2c(nyi))
5372 CALL init(volgrid6d_out) ! initialize to empty
5373 CALL raise_error()
5374 RETURN
5375 ENDIF
5376
5377! check vertical coordinate system
5378 CALL get_val(this, input_levtype=input_levtype)
5379 mask_in = & ! implicit allocation
5380 (volgrid6d_coord_in%level(:)%level1 == input_levtype%level1) .AND. &
5381 (volgrid6d_coord_in%level(:)%level2 == input_levtype%level2)
5382 ulstart = firsttrue(mask_in)
5383 ulend = lasttrue(mask_in)
5384 IF (ulstart == 0 .OR. ulend == 0) THEN
5385 CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
5386 'coordinate file does not contain levels of type '// &
5387 t2c(input_levtype%level1)//'/'//t2c(input_levtype%level2)// &
5388 ' specified for input data')
5389 CALL init(volgrid6d_out) ! initialize to empty
5390 CALL raise_error()
5391 RETURN
5392 ENDIF
5393
5394 coord_3d_in = volgrid6d_coord_in%voldati(:,:,ulstart:ulend,1,1,var_coord_in) ! implicit allocation
5395! special case
5396 IF (output_levtype%level1 == 103 .OR. &
5397 output_levtype%level1 == 108) THEN ! surface coordinate needed
5398 spos = firsttrue(volgrid6d_coord_in%level(:) == vol7d_level_new(1))
5399 IF (spos == 0) THEN
5400 CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
5401 'output level '//t2c(output_levtype%level1)// &
5402 ' requested, but height/press of surface not provided in coordinate file')
5403 CALL init(volgrid6d_out) ! initialize to empty
5404 CALL raise_error()
5405 RETURN
5406 ENDIF
5407 DO k = 1, SIZE(coord_3d_in,3)
5408 WHERE(c_e(coord_3d_in(:,:,k)) .AND. &
5409 c_e(volgrid6d_coord_in%voldati(:,:,spos,1,1,var_coord_in)))
5410 coord_3d_in(:,:,k) = coord_3d_in(:,:,k) - &
5411 volgrid6d_coord_in%voldati(:,:,spos,1,1,var_coord_in)
5412 ELSEWHERE
5413 coord_3d_in(:,:,k) = rmiss
5414 END WHERE
5415 ENDDO
5416 ENDIF
5417
5418 ENDIF
5419 ENDIF
5420
5421 IF (.NOT.c_e(var_coord_in)) THEN ! search for coordinate within volume
5422! search for variable providing vertical coordinate
5423 CALL get_val(this, output_levtype=output_levtype)
5424 vcoord_var = vol7d_var_new(vol7d_level_to_var(output_levtype))
5425 IF (c_e(vcoord_var)) THEN
5426 DO i = 1, SIZE(volgrid6d_in%var)
5427 IF (convert(volgrid6d_in%var(i)) == vcoord_var) THEN
5428 var_coord_vol = i
5429 EXIT
5430 ENDIF
5431 ENDDO
5432
5433 IF (c_e(var_coord_vol)) THEN
5434 CALL l4f_category_log(volgrid6d_in%category, l4f_info, &
5435 'Coordinate for vertint found in input volume at position '// &
5436 t2c(var_coord_vol))
5437 ENDIF
5438
5439 ENDIF
5440 ENDIF
5441
5442 CALL init(volgrid6d_out, griddim=volgrid6d_in%griddim, &
5443 time_definition=volgrid6d_in%time_definition, categoryappend=categoryappend)
5444 IF (c_e(var_coord_in)) THEN
5445 CALL init(grid_trans, this, lev_in=volgrid6d_in%level, lev_out=lev_out, &
5446 coord_3d_in=coord_3d_in, categoryappend=categoryappend)
5447 ELSE
5448 CALL init(grid_trans, this, lev_in=volgrid6d_in%level, lev_out=lev_out, &
5449 categoryappend=categoryappend)
5450 ENDIF
5451
5452 CALL get_val(grid_trans, output_level_auto=llev_out) ! get levels if auto-generated
5453 IF (.NOT.ASSOCIATED(llev_out)) llev_out => lev_out
5454 nlevel = SIZE(llev_out)
5455 ELSE
5456 CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
5457 'volgrid6d_transform: vertint requested but lev_out not provided')
5458 CALL init(volgrid6d_out) ! initialize to empty
5459 CALL raise_error()
5460 RETURN
5461 ENDIF
5462
5463ELSE
5464 CALL init(volgrid6d_out, griddim=griddim, &
5465 time_definition=volgrid6d_in%time_definition, categoryappend=categoryappend)
5466 CALL init(grid_trans, this, in=volgrid6d_in%griddim, out=volgrid6d_out%griddim, &
5467 maskgrid=maskgrid, maskbounds=maskbounds, categoryappend=categoryappend)
5468ENDIF
5469
5470
5471IF (c_e(grid_trans)) THEN ! transformation is valid
5472
5473 CALL volgrid6d_alloc(volgrid6d_out, ntime=ntime, nlevel=nlevel, &
5474 ntimerange=ntimerange, nvar=nvar)
5475
5476 IF (PRESENT(decode)) THEN ! explicitly set decode status
5477 ldecode = decode
5478 ELSE ! take it from input
5479 ldecode = ASSOCIATED(volgrid6d_in%voldati)
5480 ENDIF
5481! force decode if gaid is readonly
5482 decode_loop: DO i6 = 1,nvar
5483 DO i5 = 1, ntimerange
5484 DO i4 = 1, ntime
5485 DO i3 = 1, nlevel
5486 IF (c_e(volgrid6d_in%gaid(i3,i4,i5,i6))) THEN
5487 ldecode = ldecode .OR. grid_id_readonly(volgrid6d_in%gaid(i3,i4,i5,i6))
5488 EXIT decode_loop
5489 ENDIF
5490 ENDDO
5491 ENDDO
5492 ENDDO
5493 ENDDO decode_loop
5494
5495 IF (PRESENT(decode)) THEN
5496 IF (ldecode.NEQV.decode) THEN
5497 CALL l4f_category_log(volgrid6d_in%category, l4f_warn, &
5498 'volgrid6d_transform: decode status forced to .TRUE. because driver does not allow copy')
5499 ENDIF
5500 ENDIF
5501
5502 CALL volgrid6d_alloc_vol(volgrid6d_out, decode=ldecode)
5503
5504!ensure unproj was called
5505!call griddim_unproj(volgrid6d_out%griddim)
5506
5507 IF (trans_type == 'vertint') THEN
5508#ifdef DEBUG
5509 CALL l4f_category_log(volgrid6d_in%category, l4f_debug, &
5510 "volgrid6d_transform: vertint to "//t2c(nlevel)//" levels")
5511#endif
5512 CALL compute(grid_trans, volgrid6d_in, volgrid6d_out, lev_out=llev_out, &
5513 var_coord_vol=var_coord_vol, clone=clone)
5514 ELSE
5515 CALL compute(grid_trans, volgrid6d_in, volgrid6d_out, clone=clone)
5516 ENDIF
5517
5518 IF (cf_out == 0) THEN ! unrotated components are desired
5519 CALL wind_unrot(volgrid6d_out) ! unrotate if necessary
5520 ELSE IF (cf_out == 1) THEN ! rotated components are desired
5521 CALL wind_rot(volgrid6d_out) ! rotate if necessary
5522 ENDIF
5523
5524ELSE
5525! should log with grid_trans%category, but it is private
5526 CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
5527 'volgrid6d_transform: transformation not valid')
5528 CALL raise_error()
5529ENDIF
5530
5531CALL delete (grid_trans)
5532
5533END SUBROUTINE volgrid6d_transform
5534
5535
5536!> Performs the specified abstract transformation on the arrays of
5537!! data provided. The abstract transformation is specified by \a this
5538!! parameter; the corresponding specifical transformation (\a
5539!! grid_transform object) is created and destroyed internally. The
5540!! output transformed object is created internally and it does not
5541!! require preliminary initialisation. According to the input data and
5542!! to the transformation type, the output array may have of one or
5543!! more \a volgrid6d elements on different grids.
5544SUBROUTINE volgrid6dv_transform(this, griddim, volgrid6d_in, volgrid6d_out, &
5545 lev_out, volgrid6d_coord_in, maskgrid, maskbounds, clone, decode, categoryappend)
5546TYPE(transform_def),INTENT(in) :: this !< object specifying the abstract transformation
5547TYPE(griddim_def),INTENT(in),OPTIONAL :: griddim !< griddim specifying the output grid (required by most transformation types)
5548TYPE(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)
5549TYPE(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
5550TYPE(vol7d_level),INTENT(in),OPTIONAL :: lev_out(:) !< vol7d_level object defining target vertical grid
5551TYPE(volgrid6d),INTENT(in),OPTIONAL :: volgrid6d_coord_in !< object providing time constant input vertical coordinate for some kind of vertical interpolations
5552REAL,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')
5553REAL,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')
5554LOGICAL,INTENT(in),OPTIONAL :: clone !< if provided and \a .TRUE. , clone the \a gaid's from \a volgrid6d_in to \a volgrid6d_out
5555LOGICAL,INTENT(in),OPTIONAL :: decode !< if provided and \a .FALSE. the data volume is not allocated, but work is performed on grid_id's
5556CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
5557
5558INTEGER :: i, stallo
5559
5560
5561allocate(volgrid6d_out(size(volgrid6d_in)),stat=stallo)
5562if (stallo /= 0)then
5563 call l4f_log(l4f_fatal,"allocating memory")
5564 call raise_fatal_error()
5565end if
5566
5567do i=1,size(volgrid6d_in)
5568 call transform(this, griddim, volgrid6d_in(i), volgrid6d_out(i), &
5569 lev_out=lev_out, volgrid6d_coord_in=volgrid6d_coord_in, &
5570 maskgrid=maskgrid, maskbounds=maskbounds, &
5571 clone=clone, decode=decode, categoryappend=categoryappend)
5572end do
5573
5574END SUBROUTINE volgrid6dv_transform
5575
5576
5577! Internal method for performing grid to sparse point computations
5578SUBROUTINE volgrid6d_v7d_transform_compute(this, volgrid6d_in, vol7d_out, &
5579 networkname, noconvert)
5580TYPE(grid_transform),INTENT(in) :: this ! oggetto di trasformazione per grigliato
5581type(volgrid6d), INTENT(in) :: volgrid6d_in ! oggetto da trasformare
5582type(vol7d), INTENT(inout) :: vol7d_out ! oggetto trasformato
5583CHARACTER(len=*),OPTIONAL,INTENT(in) :: networkname ! imposta il network in vol7d_out (default='generic')
5584LOGICAL,OPTIONAL,INTENT(in) :: noconvert !< do not try to match variable and convert values during transform
5585
5586INTEGER :: nntime, nana, ntime, ntimerange, nlevel, nvar, stallo
5587INTEGER :: itime, itimerange, ivar, inetwork
5588REAL,ALLOCATABLE :: voldatir_out(:,:,:)
5589TYPE(conv_func),POINTER :: c_func(:)
5590TYPE(datetime),ALLOCATABLE :: validitytime(:,:)
5591REAL,POINTER :: voldatiin(:,:,:)
5592
5593#ifdef DEBUG
5594call l4f_category_log(volgrid6d_in%category,l4f_debug,"start volgrid6d_v7d_transform_compute")
5595#endif
5596
5597ntime=0
5598ntimerange=0
5599nlevel=0
5600nvar=0
5601NULLIFY(c_func)
5602
5603if (present(networkname))then
5604 call init(vol7d_out%network(1),name=networkname)
5605else
5606 call init(vol7d_out%network(1),name='generic')
5607end if
5608
5609if (associated(volgrid6d_in%timerange))then
5610 ntimerange=size(volgrid6d_in%timerange)
5611 vol7d_out%timerange=volgrid6d_in%timerange
5612end if
5613
5614if (associated(volgrid6d_in%time))then
5615 ntime=size(volgrid6d_in%time)
5616
5617 if (vol7d_out%time_definition == volgrid6d_in%time_definition) then
5618
5619 ! i time sono definiti uguali: assegno
5620 vol7d_out%time=volgrid6d_in%time
5621
5622 else
5623 ! converto reference in validity
5624 allocate (validitytime(ntime,ntimerange),stat=stallo)
5625 if (stallo /=0)then
5626 call l4f_category_log(volgrid6d_in%category,l4f_fatal,"allocating memory")
5627 call raise_fatal_error()
5628 end if
5629
5630 do itime=1,ntime
5631 do itimerange=1,ntimerange
5632 if (vol7d_out%time_definition > volgrid6d_in%time_definition) then
5633 validitytime(itime,itimerange) = &
5634 volgrid6d_in%time(itime) + timedelta_new(sec=volgrid6d_in%timerange(itimerange)%p1)
5635 else
5636 validitytime(itime,itimerange) = &
5637 volgrid6d_in%time(itime) - timedelta_new(sec=volgrid6d_in%timerange(itimerange)%p1)
5638 end if
5639 end do
5640 end do
5641
5642 nntime = count_distinct(reshape(validitytime,(/ntime*ntimerange/)), back=.true.)
5643 vol7d_out%time=pack_distinct(reshape(validitytime,(/ntime*ntimerange/)), nntime,back=.true.)
5644
5645 end if
5646end if
5647
5648IF (ASSOCIATED(volgrid6d_in%level)) THEN
5649 nlevel = SIZE(volgrid6d_in%level)
5650 vol7d_out%level=volgrid6d_in%level
5651ENDIF
5652
5653IF (ASSOCIATED(volgrid6d_in%var)) THEN
5654 nvar = SIZE(volgrid6d_in%var)
5655 IF (.NOT. optio_log(noconvert)) THEN
5656 CALL vargrib2varbufr(volgrid6d_in%var, vol7d_out%dativar%r, c_func)
5657 ENDIF
5658ENDIF
5659
5660nana = SIZE(vol7d_out%ana)
5661
5662! allocate once for speed
5663IF (.NOT.ASSOCIATED(volgrid6d_in%voldati)) THEN
5664 ALLOCATE(voldatiin(volgrid6d_in%griddim%dim%nx, volgrid6d_in%griddim%dim%ny, &
5665 nlevel))
5666ENDIF
5667
5668ALLOCATE(voldatir_out(nana,1,nlevel),stat=stallo)
5669IF (stallo /= 0) THEN
5670 CALL l4f_category_log(volgrid6d_in%category,l4f_fatal,"allocating memory")
5671 CALL raise_fatal_error()
5672ENDIF
5673
5674inetwork=1
5675do itime=1,ntime
5676 do itimerange=1,ntimerange
5677! do ilevel=1,nlevel
5678 do ivar=1,nvar
5679
5680 !non è chiaro se questa sezione è utile o no
5681 !ossia il compute sotto sembra prevedere voldatir_out solo in out
5682!!$ if (vol7d_out%time_definition == volgrid6d_in%time_definition) then
5683!!$ voldatir_out=reshape(vol7d_out%voldatir(:,itime,ilevel,itimerange,ivar,inetwork),(/nana,1/))
5684!!$ else
5685!!$ voldatir_out=reshape(vol7d_out%voldatir(:,index(vol7d_out%time,validitytime(itime,ilevel)),ilevel,itimerange,ivar,inetwork),(/nana,1/))
5686!!$ end if
5687
5688 CALL volgrid_get_vol_3d(volgrid6d_in, itime, itimerange, ivar, &
5689 voldatiin)
5690
5691 CALL compute(this, voldatiin, voldatir_out, vol7d_out%dativar%r(ivar))
5692
5693 if (vol7d_out%time_definition == volgrid6d_in%time_definition) then
5694 vol7d_out%voldatir(:,itime,:,itimerange,ivar,inetwork) = &
5695 voldatir_out(:,1,:)
5696 else
5697 vol7d_out%voldatir(:,index(vol7d_out%time,validitytime(itime,itimerange)),:,itimerange,ivar,inetwork)=&
5698 reshape(voldatir_out,(/nana,nlevel/))
5699 end if
5700
5701! 1 indice della dimensione "anagrafica"
5702! 2 indice della dimensione "tempo"
5703! 3 indice della dimensione "livello verticale"
5704! 4 indice della dimensione "intervallo temporale"
5705! 5 indice della dimensione "variabile"
5706! 6 indice della dimensione "rete"
5707
5708 end do
5709! end do
5710 end do
5711end do
5712
5713deallocate(voldatir_out)
5714IF (.NOT.ASSOCIATED(volgrid6d_in%voldati)) THEN
5715 DEALLOCATE(voldatiin)
5716ENDIF
5717if (allocated(validitytime)) deallocate(validitytime)
5718
5719! Rescale valid data according to variable conversion table
5720IF (ASSOCIATED(c_func)) THEN
5721 DO ivar = 1, nvar
5722 CALL compute(c_func(ivar), vol7d_out%voldatir(:,:,:,:,ivar,:))
5723 ENDDO
5724 DEALLOCATE(c_func)
5725ENDIF
5726
5727end SUBROUTINE volgrid6d_v7d_transform_compute
5728
5729
5730!> Performs the specified abstract transformation on the data provided.
5731!! The abstract transformation is specified by \a this parameter; the
5732!! corresponding specifical transformation (\a grid_transform object)
5733!! is created and destroyed internally. The output transformed object
5734!! is created internally and it does not require preliminary
5735!! initialisation.
5736SUBROUTINE volgrid6d_v7d_transform(this, volgrid6d_in, vol7d_out, v7d, &
5737 maskgrid, maskbounds, networkname, noconvert, find_index, categoryappend)
5738TYPE(transform_def),INTENT(in) :: this !< object specifying the abstract transformation
5739TYPE(volgrid6d),INTENT(inout) :: volgrid6d_in !< object to be transformed, it is not modified, despite the INTENT(inout)
5740TYPE(vol7d),INTENT(out) :: vol7d_out !< transformed object, it does not requires initialisation
5741TYPE(vol7d),INTENT(in),OPTIONAL :: v7d !< object containing a list of points over which transformation has to be done (required by some transformation types)
5742REAL,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')
5743REAL,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
5744CHARACTER(len=*),OPTIONAL,INTENT(in) :: networkname !< set the output network name in vol7d_out (default='generic')
5745LOGICAL,OPTIONAL,INTENT(in) :: noconvert !< do not try to match variable and convert values during transform
5746PROCEDURE(basic_find_index),POINTER,OPTIONAL :: find_index
5747CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
5748
5749type(grid_transform) :: grid_trans
5750INTEGER :: ntime, ntimerange, nlevel, nvar, nana, time_definition, nnetwork, stallo
5751INTEGER :: itime, itimerange, inetwork
5752TYPE(datetime),ALLOCATABLE :: validitytime(:,:)
5753INTEGER,ALLOCATABLE :: point_index(:)
5754TYPE(vol7d) :: v7d_locana
5755
5756#ifdef DEBUG
5757call l4f_category_log(volgrid6d_in%category,l4f_debug,"start volgrid6d_v7d_transform")
5758#endif
5759
5760call vg6d_wind_unrot(volgrid6d_in)
5761
5762ntime=0
5763ntimerange=0
5764nlevel=0
5765nvar=0
5766nnetwork=1
5767
5768call get_val(this,time_definition=time_definition)
5769if (.not. c_e(time_definition)) then
5770 time_definition=1 ! default to validity time
5771endif
5772
5773IF (PRESENT(v7d)) THEN
5774 CALL vol7d_copy(v7d, v7d_locana)
5775ELSE
5776 CALL init(v7d_locana)
5777ENDIF
5778
5779if (associated(volgrid6d_in%timerange)) ntimerange=size(volgrid6d_in%timerange)
5780
5781if (associated(volgrid6d_in%time)) then
5782
5783 ntime=size(volgrid6d_in%time)
5784
5785 if (time_definition /= volgrid6d_in%time_definition) then
5786
5787 ! converto reference in validity
5788 allocate (validitytime(ntime,ntimerange),stat=stallo)
5789 if (stallo /=0)then
5790 call l4f_category_log(volgrid6d_in%category,l4f_fatal,"allocating memory")
5791 call raise_fatal_error()
5792 end if
5793
5794 do itime=1,ntime
5795 do itimerange=1,ntimerange
5796 if (time_definition > volgrid6d_in%time_definition) then
5797 validitytime(itime,itimerange) = &
5798 volgrid6d_in%time(itime) + timedelta_new(sec=volgrid6d_in%timerange(itimerange)%p1)
5799 else
5800 validitytime(itime,itimerange) = &
5801 volgrid6d_in%time(itime) - timedelta_new(sec=volgrid6d_in%timerange(itimerange)%p1)
5802 end if
5803 end do
5804 end do
5805
5806 ntime = count_distinct(reshape(validitytime,(/ntime*ntimerange/)), back=.true.)
5807 deallocate (validitytime)
5808
5809 end if
5810end if
5811
5812
5813if (associated(volgrid6d_in%level)) nlevel=size(volgrid6d_in%level)
5814if (associated(volgrid6d_in%var)) nvar=size(volgrid6d_in%var)
5815
5816CALL init(grid_trans, this, volgrid6d_in%griddim, v7d_locana, &
5817 maskgrid=maskgrid, maskbounds=maskbounds, find_index=find_index, &
5818 categoryappend=categoryappend)
5819CALL init (vol7d_out,time_definition=time_definition)
5820
5821IF (c_e(grid_trans)) THEN
5822
5823 nana=SIZE(v7d_locana%ana)
5824 CALL vol7d_alloc(vol7d_out, nana=nana, ntime=ntime, nlevel=nlevel, &
5825 ntimerange=ntimerange, ndativarr=nvar, nnetwork=nnetwork)
5826 vol7d_out%ana = v7d_locana%ana
5827
5828 CALL get_val(grid_trans, output_point_index=point_index)
5829 IF (ALLOCATED(point_index)) THEN
5830! check that size(point_index) == nana?
5831 CALL vol7d_alloc(vol7d_out, nanavari=1)
5832 CALL init(vol7d_out%anavar%i(1), 'B01192')
5833 ENDIF
5834
5835 CALL vol7d_alloc_vol(vol7d_out)
5836
5837 IF (ALLOCATED(point_index)) THEN
5838 DO inetwork = 1, nnetwork
5839 vol7d_out%volanai(:,1,inetwork) = point_index(:)
5840 ENDDO
5841 ENDIF
5842 CALL compute(grid_trans, volgrid6d_in, vol7d_out, networkname, noconvert)
5843ELSE
5844 CALL l4f_log(l4f_error, 'vg6d_v7d_transform: transformation not valid')
5845 CALL raise_error()
5846ENDIF
5847
5848CALL delete(grid_trans)
5849
5850#ifdef HAVE_DBALLE
5851! set variables to a conformal state
5852CALL vol7d_dballe_set_var_du(vol7d_out)
5853#endif
5854
5855CALL delete(v7d_locana)
5856
5857END SUBROUTINE volgrid6d_v7d_transform
5858
5859
5860!> Performs the specified abstract transformation on the arrays of
5861!! data provided. The abstract transformation is specified by \a this
5862!! parameter; the corresponding specifical transformation (\a
5863!! grid_transform object) is created and destroyed internally. The
5864!! output transformed object is created internally and it does not
5865!! require preliminary initialisation. The transformation performed on
5866!! each element of the input \a volgrid6d array object is merged into
5867!! a single \a vol7d output object.
5868SUBROUTINE volgrid6dv_v7d_transform(this, volgrid6d_in, vol7d_out, v7d, &
5869 maskgrid, maskbounds, networkname, noconvert, find_index, categoryappend)
5870TYPE(transform_def),INTENT(in) :: this !< object specifying the abstract transformation
5871TYPE(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)
5872TYPE(vol7d),INTENT(out) :: vol7d_out !< transformed object, it does not require initialisation
5873TYPE(vol7d),INTENT(in),OPTIONAL :: v7d !< object containing a list of points over which transformation has to be done (required by some transformation types)
5874REAL,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')
5875REAL,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
5876CHARACTER(len=*),OPTIONAL,INTENT(in) :: networkname !< set the output network name in vol7d_out (default='generic')
5877LOGICAL,OPTIONAL,INTENT(in) :: noconvert !< do not try to match variable and convert values during transform
5878PROCEDURE(basic_find_index),POINTER,OPTIONAL :: find_index
5879CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
5880
5881integer :: i
5882TYPE(vol7d) :: v7dtmp
5883
5884
5885CALL init(v7dtmp)
5886CALL init(vol7d_out)
5887
5888DO i=1,SIZE(volgrid6d_in)
5889 CALL transform(this, volgrid6d_in(i), v7dtmp, v7d=v7d, &
5890 maskgrid=maskgrid, maskbounds=maskbounds, &
5891 networkname=networkname, noconvert=noconvert, find_index=find_index, &
5892 categoryappend=categoryappend)
5893 CALL vol7d_append(vol7d_out, v7dtmp)
5894ENDDO
5895
5896END SUBROUTINE volgrid6dv_v7d_transform
5897
5898
5899! Internal method for performing sparse point to grid computations
5900SUBROUTINE v7d_volgrid6d_transform_compute(this, vol7d_in, volgrid6d_out, networkname, gaid_template)
5901TYPE(grid_transform),INTENT(in) :: this ! object specifying the specific transformation
5902type(vol7d), INTENT(in) :: vol7d_in ! object to be transformed
5903type(volgrid6d), INTENT(inout) :: volgrid6d_out ! transformed object
5904CHARACTER(len=*),OPTIONAL,INTENT(in) :: networkname ! select the network to be processed from the \a vol7d input object, default the first network
5905TYPE(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
5906
5907integer :: nana, ntime, ntimerange, nlevel, nvar
5908INTEGER :: ilevel, itime, itimerange, ivar, inetwork
5909
5910REAL,POINTER :: voldatiout(:,:,:)
5911type(vol7d_network) :: network
5912TYPE(conv_func), pointer :: c_func(:)
5913!TODO category sarebbe da prendere da vol7d
5914#ifdef DEBUG
5915CALL l4f_category_log(volgrid6d_out%category, l4f_debug, &
5916 'start v7d_volgrid6d_transform_compute')
5917#endif
5918
5919ntime=0
5920ntimerange=0
5921nlevel=0
5922nvar=0
5923
5924IF (PRESENT(networkname)) THEN
5925 CALL init(network,name=networkname)
5926 inetwork = index(vol7d_in%network,network)
5927 IF (inetwork <= 0) THEN
5928 CALL l4f_category_log(volgrid6d_out%category,l4f_warn, &
5929 'network '//trim(networkname)//' not found, first network will be transformed')
5930 inetwork = 1
5931 ENDIF
5932ELSE
5933 inetwork = 1
5934ENDIF
5935
5936! no time_definition conversion implemented, output will be the same as input
5937if (associated(vol7d_in%time))then
5938 ntime=size(vol7d_in%time)
5939 volgrid6d_out%time=vol7d_in%time
5940end if
5941
5942if (associated(vol7d_in%timerange))then
5943 ntimerange=size(vol7d_in%timerange)
5944 volgrid6d_out%timerange=vol7d_in%timerange
5945end if
5946
5947if (associated(vol7d_in%level))then
5948 nlevel=size(vol7d_in%level)
5949 volgrid6d_out%level=vol7d_in%level
5950end if
5951
5952if (associated(vol7d_in%dativar%r))then
5953 nvar=size(vol7d_in%dativar%r)
5954 CALL varbufr2vargrib(vol7d_in%dativar%r, volgrid6d_out%var, c_func, gaid_template)
5955end if
5956
5957nana=SIZE(vol7d_in%voldatir, 1)
5958! allocate once for speed
5959IF (.NOT.ASSOCIATED(volgrid6d_out%voldati)) THEN
5960 ALLOCATE(voldatiout(volgrid6d_out%griddim%dim%nx, volgrid6d_out%griddim%dim%ny, &
5961 nlevel))
5962ENDIF
5963
5964DO ivar=1,nvar
5965 DO itimerange=1,ntimerange
5966 DO itime=1,ntime
5967
5968! clone the gaid template where I have data
5969 IF (PRESENT(gaid_template)) THEN
5970 DO ilevel = 1, nlevel
5971 IF (any(c_e(vol7d_in%voldatir(:,itime,ilevel,itimerange,ivar,inetwork)))) THEN
5972 CALL copy(gaid_template, volgrid6d_out%gaid(ilevel,itime,itimerange,ivar))
5973 ELSE
5974 volgrid6d_out%gaid(ilevel,itime,itimerange,ivar) = grid_id_new()
5975 ENDIF
5976 ENDDO
5977 ENDIF
5978
5979! get data
5980 IF (ASSOCIATED(volgrid6d_out%voldati)) & ! improve!!!!
5981 CALL volgrid_get_vol_3d(volgrid6d_out, itime, itimerange, ivar, &
5982 voldatiout)
5983! do the interpolation
5984 CALL compute(this, &
5985 vol7d_in%voldatir(:,itime,:,itimerange,ivar,inetwork), voldatiout, &
5986 vol7d_in%dativar%r(ivar))
5987! rescale valid data according to variable conversion table
5988 IF (ASSOCIATED(c_func)) THEN
5989 CALL compute(c_func(ivar), voldatiout(:,:,:))
5990 ENDIF
5991! put data
5992 CALL volgrid_set_vol_3d(volgrid6d_out, itime, itimerange, ivar, &
5993 voldatiout)
5994
5995 ENDDO
5996 ENDDO
5997ENDDO
5998
5999IF (.NOT.ASSOCIATED(volgrid6d_out%voldati)) THEN
6000 DEALLOCATE(voldatiout)
6001ENDIF
6002IF (ASSOCIATED(c_func)) THEN
6003 DEALLOCATE(c_func)
6004ENDIF
6005
6006END SUBROUTINE v7d_volgrid6d_transform_compute
6007
6008
6009!> Performs the specified abstract transformation on the data provided.
6010!! The abstract transformation is specified by \a this parameter; the
6011!! corresponding specifical transformation (\a grid_transform object)
6012!! is created and destroyed internally. The output transformed object
6013!! is created internally and it does not require preliminary
6014!! initialisation.
6015SUBROUTINE v7d_volgrid6d_transform(this, griddim, vol7d_in, volgrid6d_out, &
6016 networkname, gaid_template, categoryappend)
6017TYPE(transform_def),INTENT(in) :: this !< object specifying the abstract transformation
6018TYPE(griddim_def),INTENT(in),OPTIONAL :: griddim !< griddim specifying the output grid (required by most transformation types)
6019! TODO ripristinare intent(in) dopo le opportune modifiche in grid_class.F90
6020TYPE(vol7d),INTENT(inout) :: vol7d_in !< object to be transformed, it is not modified, despite the INTENT(inout)
6021TYPE(volgrid6d),INTENT(out) :: volgrid6d_out !< transformed object, it does not require initialisation
6022CHARACTER(len=*),OPTIONAL,INTENT(in) :: networkname !< select the network to be processed from the \a vol7d input object, default the first network
6023TYPE(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
6024CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
6025
6026type(grid_transform) :: grid_trans
6027integer :: ntime, ntimerange, nlevel, nvar
6028
6029
6030!TODO la category sarebbe da prendere da vol7d
6031!call l4f_category_log(vol7d_out%category,L4F_DEBUG,"start volgrid6d_transform")
6032
6033CALL vol7d_alloc_vol(vol7d_in) ! be safe
6034ntime=SIZE(vol7d_in%time)
6035ntimerange=SIZE(vol7d_in%timerange)
6036nlevel=SIZE(vol7d_in%level)
6037nvar=0
6038if (associated(vol7d_in%dativar%r)) nvar=size(vol7d_in%dativar%r)
6039
6040IF (nvar <= 0) THEN ! use vol7d category once it will be implemented
6041 CALL l4f_log(l4f_error, &
6042 "trying to transform a vol7d object incomplete or without real variables")
6043 CALL init(volgrid6d_out) ! initialize to empty
6044 CALL raise_error()
6045 RETURN
6046ENDIF
6047
6048CALL init(grid_trans, this, vol7d_in, griddim, categoryappend=categoryappend)
6049CALL init(volgrid6d_out, griddim, time_definition=vol7d_in%time_definition, &
6050 categoryappend=categoryappend)
6051
6052IF (c_e(grid_trans)) THEN
6053
6054 CALL volgrid6d_alloc(volgrid6d_out, griddim%dim, ntime=ntime, nlevel=nlevel, &
6055 ntimerange=ntimerange, nvar=nvar)
6056! can I avoid decode=.TRUE. here, using gaid_template?
6057 CALL volgrid6d_alloc_vol(volgrid6d_out, decode=.true.)
6058
6059 CALL compute(grid_trans, vol7d_in, volgrid6d_out, networkname, gaid_template)
6060
6061 CALL vg6d_wind_rot(volgrid6d_out)
6062ELSE
6063! should log with grid_trans%category, but it is private
6064 CALL l4f_log(l4f_error, 'v7d_vg6d_transform: transformation not valid')
6065 CALL raise_error()
6066ENDIF
6067
6068CALL delete(grid_trans)
6069
6070END SUBROUTINE v7d_volgrid6d_transform
6071
6072
6073! Internal method for performing sparse point to sparse point computations
6074SUBROUTINE v7d_v7d_transform_compute(this, vol7d_in, vol7d_out, lev_out, &
6075 var_coord_vol)
6076TYPE(grid_transform),INTENT(in) :: this ! oggetto di trasformazione per grigliato
6077type(vol7d), INTENT(in) :: vol7d_in ! oggetto da trasformare
6078type(vol7d), INTENT(inout) :: vol7d_out ! oggetto trasformato
6079TYPE(vol7d_level),INTENT(in),OPTIONAL :: lev_out(:) ! vol7d_level object defining target vertical grid, for vertical interpolations
6080INTEGER,INTENT(in),OPTIONAL :: var_coord_vol ! index of variable defining vertical coordinate values in input volume
6081
6082INTEGER :: itime, itimerange, ilevel, ivar, inetwork, &
6083 levshift, levused, lvar_coord_vol, spos
6084REAL,ALLOCATABLE :: coord_3d_in(:,:,:)
6085TYPE(vol7d_level) :: output_levtype
6086
6087lvar_coord_vol = optio_i(var_coord_vol)
6088vol7d_out%time(:) = vol7d_in%time(:)
6089vol7d_out%timerange(:) = vol7d_in%timerange(:)
6090IF (PRESENT(lev_out)) THEN
6091 vol7d_out%level(:) = lev_out(:)
6092ELSE
6093 vol7d_out%level(:) = vol7d_in%level(:)
6094ENDIF
6095vol7d_out%network(:) = vol7d_in%network(:)
6096IF (ASSOCIATED(vol7d_in%dativar%r)) THEN ! work only when real vars are available
6097 vol7d_out%dativar%r(:) = vol7d_in%dativar%r(:)
6098
6099 CALL get_val(this, levshift=levshift, levused=levused)
6100 spos = imiss
6101 IF (c_e(lvar_coord_vol)) THEN
6102 CALL get_val(this%trans, output_levtype=output_levtype)
6103 IF (output_levtype%level1 == 103 .OR. output_levtype%level1 == 108) THEN
6104 spos = firsttrue(vol7d_in%level(:) == vol7d_level_new(1))
6105 IF (spos == 0) THEN
6106 CALL l4f_log(l4f_error, &
6107 'output level '//t2c(output_levtype%level1)// &
6108 ' requested, but height/press of surface not provided in volume')
6109 ENDIF
6110 IF (.NOT.c_e(levshift) .AND. .NOT.c_e(levused)) THEN
6111 CALL l4f_log(l4f_error, &
6112 'internal inconsistence, levshift and levused undefined when they should be')
6113 ENDIF
6114 ALLOCATE(coord_3d_in(SIZE(vol7d_in%ana),1,SIZE(vol7d_in%level)))
6115 ENDIF
6116
6117 ENDIF
6118
6119 DO inetwork = 1, SIZE(vol7d_in%network)
6120 DO ivar = 1, SIZE(vol7d_in%dativar%r)
6121 DO itimerange = 1, SIZE(vol7d_in%timerange)
6122 DO itime = 1, SIZE(vol7d_in%time)
6123
6124! dirty trick to make voldatir look like a 2d-array of shape (nana,1)
6125 IF (c_e(lvar_coord_vol)) THEN
6126 IF (c_e(spos)) THEN ! compute difference wrt surface coordinate
6127 IF (spos == 0) THEN ! error condition, set all to missing and goodnight
6128 coord_3d_in(:,:,levshift+1:levshift+levused) = rmiss
6129 ELSE
6130 DO ilevel = levshift+1, levshift+levused
6131 WHERE(c_e(vol7d_in%voldatir(:,itime:itime,ilevel,itimerange,lvar_coord_vol,inetwork)) .AND. &
6132 c_e(vol7d_in%voldatir(:,itime:itime,spos,itimerange,lvar_coord_vol,inetwork)))
6133 coord_3d_in(:,:,ilevel) = vol7d_in%voldatir(:,itime:itime,ilevel,itimerange,lvar_coord_vol,inetwork) - &
6134 vol7d_in%voldatir(:,itime:itime,spos,itimerange,lvar_coord_vol,inetwork)
6135 ELSEWHERE
6136 coord_3d_in(:,:,ilevel) = rmiss
6137 END WHERE
6138 ENDDO
6139 ENDIF
6140 CALL compute(this, &
6141 vol7d_in%voldatir(:,itime,:,itimerange,ivar,inetwork), &
6142 vol7d_out%voldatir(:,itime:itime,:,itimerange,ivar,inetwork), &
6143 var=vol7d_in%dativar%r(ivar), &
6144 coord_3d_in=coord_3d_in)
6145 ELSE
6146 CALL compute(this, &
6147 vol7d_in%voldatir(:,itime,:,itimerange,ivar,inetwork), &
6148 vol7d_out%voldatir(:,itime:itime,:,itimerange,ivar,inetwork), &
6149 var=vol7d_in%dativar%r(ivar), &
6150 coord_3d_in=vol7d_in%voldatir(:,itime:itime,:,itimerange, &
6151 lvar_coord_vol,inetwork))
6152 ENDIF
6153 ELSE
6154 CALL compute(this, &
6155 vol7d_in%voldatir(:,itime,:,itimerange,ivar,inetwork), &
6156 vol7d_out%voldatir(:,itime:itime,:,itimerange,ivar,inetwork), &
6157 var=vol7d_in%dativar%r(ivar))
6158 ENDIF
6159 ENDDO
6160 ENDDO
6161 ENDDO
6162 ENDDO
6163
6164ENDIF
6165
6166END SUBROUTINE v7d_v7d_transform_compute
6167
6168
6169!> Performs the specified abstract transformation on the data provided.
6170!! The abstract transformation is specified by \a this parameter; the
6171!! corresponding specifical transformation (\a grid_transform object)
6172!! is created and destroyed internally. The output transformed object
6173!! is created internally and it does not require preliminary
6174!! initialisation. The success of the transformation can be checked
6175!! with the \a c_e method: c_e(vol7d_out).
6176SUBROUTINE v7d_v7d_transform(this, vol7d_in, vol7d_out, v7d, maskbounds, &
6177 lev_out, vol7d_coord_in, categoryappend)
6178TYPE(transform_def),INTENT(in) :: this !< object specifying the abstract transformation
6179TYPE(vol7d),INTENT(inout) :: vol7d_in !< object to be transformed, it is not modified, despite the INTENT(inout)
6180TYPE(vol7d),INTENT(out) :: vol7d_out !< transformed object, it does not require initialisation
6181TYPE(vol7d),INTENT(in),OPTIONAL :: v7d !< object containing a list of points over which transformation has to be done (required by some transformation types)
6182REAL,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')
6183TYPE(vol7d_level),INTENT(in),OPTIONAL,TARGET :: lev_out(:) !< vol7d_level object defining target vertical grid, for vertical interpolations
6184TYPE(vol7d),INTENT(in),OPTIONAL :: vol7d_coord_in !< object providing time constant input vertical coordinate for some kind of vertical interpolations
6185CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category
6186
6187INTEGER :: nvar, inetwork
6188TYPE(grid_transform) :: grid_trans
6189TYPE(vol7d_level),POINTER :: llev_out(:)
6190TYPE(vol7d_level) :: input_levtype, output_levtype
6191TYPE(vol7d_var) :: vcoord_var
6192REAL,ALLOCATABLE :: coord_3d_in(:,:,:)
6193INTEGER :: var_coord_in, var_coord_vol, i, k, ulstart, ulend, spos
6194INTEGER,ALLOCATABLE :: point_index(:)
6195TYPE(vol7d) :: v7d_locana, vol7d_tmpana
6196CHARACTER(len=80) :: trans_type
6197LOGICAL,ALLOCATABLE :: mask_in(:), point_mask(:)
6198
6199CALL vol7d_alloc_vol(vol7d_in) ! be safe
6200nvar=0
6201IF (ASSOCIATED(vol7d_in%dativar%r)) nvar=SIZE(vol7d_in%dativar%r)
6202
6203CALL init(v7d_locana)
6204IF (PRESENT(v7d)) v7d_locana = v7d
6205CALL init(vol7d_out, time_definition=vol7d_in%time_definition)
6206
6207CALL get_val(this, trans_type=trans_type)
6208
6209var_coord_vol = imiss
6210IF (trans_type == 'vertint') THEN
6211
6212 IF (PRESENT(lev_out)) THEN
6213
6214! if vol7d_coord_in provided and allocated, check that it fits
6215 var_coord_in = -1
6216 IF (PRESENT(vol7d_coord_in)) THEN
6217 IF (ASSOCIATED(vol7d_coord_in%voldatir) .AND. &
6218 ASSOCIATED(vol7d_coord_in%dativar%r)) THEN
6219
6220! strictly 1 time, 1 timerange and 1 network
6221 IF (SIZE(vol7d_coord_in%voldatir,2) /= 1 .OR. &
6222 SIZE(vol7d_coord_in%voldatir,4) /= 1 .OR. &
6223 SIZE(vol7d_coord_in%voldatir,6) /= 1) THEN
6224 CALL l4f_log(l4f_error, &
6225 'volume providing constant input vertical coordinate must have &
6226 &only 1 time, 1 timerange and 1 network')
6227 CALL raise_error()
6228 RETURN
6229 ENDIF
6230
6231! search for variable providing vertical coordinate
6232 CALL get_val(this, output_levtype=output_levtype)
6233 vcoord_var = vol7d_var_new(vol7d_level_to_var(output_levtype))
6234 IF (.NOT.c_e(vcoord_var)) THEN
6235 CALL l4f_log(l4f_error, &
6236 'requested output level type '//t2c(output_levtype%level1)// &
6237 ' does not correspond to any known physical variable for &
6238 &providing vertical coordinate')
6239 CALL raise_error()
6240 RETURN
6241 ENDIF
6242
6243 var_coord_in = index(vol7d_coord_in%dativar%r, vcoord_var)
6244
6245 IF (var_coord_in <= 0) THEN
6246 CALL l4f_log(l4f_error, &
6247 'volume providing constant input vertical coordinate contains no &
6248 &real variables matching output level type '//t2c(output_levtype%level1))
6249 CALL raise_error()
6250 RETURN
6251 ENDIF
6252 CALL l4f_log(l4f_info, &
6253 'Coordinate for vertint found in coord volume at position '// &
6254 t2c(var_coord_in))
6255
6256! check vertical coordinate system
6257 CALL get_val(this, input_levtype=input_levtype)
6258 mask_in = & ! implicit allocation
6259 (vol7d_coord_in%level(:)%level1 == input_levtype%level1) .AND. &
6260 (vol7d_coord_in%level(:)%level2 == input_levtype%level2)
6261 ulstart = firsttrue(mask_in)
6262 ulend = lasttrue(mask_in)
6263 IF (ulstart == 0 .OR. ulend == 0) THEN
6264 CALL l4f_log(l4f_error, &
6265 'coordinate file does not contain levels of type '// &
6266 t2c(input_levtype%level1)//'/'//t2c(input_levtype%level2)// &
6267 ' specified for input data')
6268 CALL raise_error()
6269 RETURN
6270 ENDIF
6271
6272 coord_3d_in = vol7d_coord_in%voldatir(:,1:1,ulstart:ulend,1,var_coord_in,1) ! dirty 1:1, implicit allocation
6273! special case
6274 IF (output_levtype%level1 == 103 &
6275 .OR. output_levtype%level1 == 108) THEN ! surface coordinate needed
6276 spos = firsttrue(vol7d_coord_in%level(:) == vol7d_level_new(1))
6277 IF (spos == 0) THEN
6278 CALL l4f_log(l4f_error, &
6279 'output level '//t2c(output_levtype%level1)// &
6280 ' requested, but height/press of surface not provided in coordinate file')
6281 CALL raise_error()
6282 RETURN
6283 ENDIF
6284 DO k = 1, SIZE(coord_3d_in,3)
6285 WHERE(c_e(coord_3d_in(:,:,k)) .AND. &
6286 c_e(vol7d_coord_in%voldatir(:,1:1,spos,1,var_coord_in,1)))
6287 coord_3d_in(:,:,k) = coord_3d_in(:,:,k) - &
6288 vol7d_coord_in%voldatir(:,1:1,spos,1,var_coord_in,1)
6289 ELSEWHERE
6290 coord_3d_in(:,:,k) = rmiss
6291 END WHERE
6292 ENDDO
6293 ENDIF
6294
6295 ENDIF
6296 ENDIF
6297
6298 IF (var_coord_in <= 0) THEN ! search for coordinate within volume
6299! search for variable providing vertical coordinate
6300 CALL get_val(this, output_levtype=output_levtype)
6301 vcoord_var = vol7d_var_new(vol7d_level_to_var(output_levtype))
6302 IF (c_e(vcoord_var)) THEN
6303 DO i = 1, SIZE(vol7d_in%dativar%r)
6304 IF (vol7d_in%dativar%r(i) == vcoord_var) THEN
6305 var_coord_vol = i
6306 EXIT
6307 ENDIF
6308 ENDDO
6309
6310 IF (c_e(var_coord_vol)) THEN
6311 CALL l4f_log(l4f_info, &
6312 'Coordinate for vertint found in input volume at position '// &
6313 t2c(var_coord_vol))
6314 ENDIF
6315
6316 ENDIF
6317 ENDIF
6318
6319 IF (var_coord_in > 0) THEN
6320 CALL init(grid_trans, this, lev_in=vol7d_in%level, lev_out=lev_out, &
6321 coord_3d_in=coord_3d_in, categoryappend=categoryappend)
6322 ELSE
6323 CALL init(grid_trans, this, lev_in=vol7d_in%level, lev_out=lev_out, &
6324 categoryappend=categoryappend)
6325 ENDIF
6326
6327 CALL get_val(grid_trans, output_level_auto=llev_out) ! get levels if auto-generated
6328 IF (.NOT.associated(llev_out)) llev_out => lev_out
6329
6330 IF (c_e(grid_trans)) then! .AND. nvar > 0) THEN
6331
6332 CALL vol7d_alloc(vol7d_out, nana=SIZE(vol7d_in%ana), &
6333 ntime=SIZE(vol7d_in%time), ntimerange=SIZE(vol7d_in%timerange), &
6334 nlevel=SIZE(llev_out), nnetwork=SIZE(vol7d_in%network), ndativarr=nvar)
6335 vol7d_out%ana(:) = vol7d_in%ana(:)
6336
6337 CALL vol7d_alloc_vol(vol7d_out)
6338
6339! no need to check c_e(var_coord_vol) here since the presence of
6340! this%coord_3d_in (external) has precedence over coord_3d_in internal
6341! in grid_transform_compute
6342 CALL compute(grid_trans, vol7d_in, vol7d_out, llev_out, &
6343 var_coord_vol=var_coord_vol)
6344 ELSE
6345 CALL l4f_log(l4f_error, 'v7d_v7d_transform: transformation not valid')
6346 CALL raise_error()
6347 ENDIF
6348 ELSE
6349 CALL l4f_log(l4f_error, &
6350 'v7d_v7d_transform: vertint requested but lev_out not provided')
6351 CALL raise_error()
6352 ENDIF
6353
6354ELSE
6355
6356 CALL init(grid_trans, this, vol7d_in, v7d_locana, maskbounds=maskbounds, &
6357 categoryappend=categoryappend)
6358! if this init is successful, I am sure that v7d_locana%ana is associated
6359
6360 IF (c_e(grid_trans)) then! .AND. nvar > 0) THEN
6361
6362 CALL vol7d_alloc(vol7d_out, nana=SIZE(v7d_locana%ana), &
6363 ntime=SIZE(vol7d_in%time), ntimerange=SIZE(vol7d_in%timerange), &
6364 nlevel=SIZE(vol7d_in%level), nnetwork=SIZE(vol7d_in%network), ndativarr=nvar)
6365 vol7d_out%ana = v7d_locana%ana
6366
6367 CALL get_val(grid_trans, point_mask=point_mask, output_point_index=point_index)
6368
6369 IF (ALLOCATED(point_index)) THEN
6370 CALL vol7d_alloc(vol7d_out, nanavari=1)
6371 CALL init(vol7d_out%anavar%i(1), 'B01192')
6372 ENDIF
6373
6374 CALL vol7d_alloc_vol(vol7d_out)
6375
6376 IF (ALLOCATED(point_index)) THEN
6377 DO inetwork = 1, SIZE(vol7d_in%network)
6378 vol7d_out%volanai(:,1,inetwork) = point_index(:)
6379 ENDDO
6380 ENDIF
6381 CALL compute(grid_trans, vol7d_in, vol7d_out)
6382
6383 IF (ALLOCATED(point_mask)) THEN ! keep full ana
6384 IF (SIZE(point_mask) /= SIZE(vol7d_in%ana)) THEN
6385 CALL l4f_log(l4f_warn, &
6386 'v7d_v7d_transform: inconsistency in point size: '//t2c(SIZE(point_mask)) &
6387 //':'//t2c(SIZE(vol7d_in%ana)))
6388 ELSE
6389#ifdef DEBUG
6390 CALL l4f_log(l4f_debug, 'v7d_v7d_transform: merging ana from in to out')
6391#endif
6392 CALL vol7d_copy(vol7d_in, vol7d_tmpana, &
6393 lana=point_mask, lnetwork=(/.true./), &
6394 ltime=(/.false./), ltimerange=(/.false./), llevel=(/.false./))
6395 CALL vol7d_append(vol7d_out, vol7d_tmpana)
6396 ENDIF
6397 ENDIF
6398
6399 ELSE
6400 CALL l4f_log(l4f_error, 'v7d_v7d_transform: transformation not valid')
6401 CALL raise_error()
6402 ENDIF
6403
6404ENDIF
6405
6406CALL delete (grid_trans)
6407IF (.NOT. PRESENT(v7d)) CALL delete(v7d_locana)
6408
6409END SUBROUTINE v7d_v7d_transform
6410
6411
6412!> Unrotate the wind components.
6413!! It converts u and v components of vector quantities relative to the
6414!! defined grid in the direction of increasing x and y coordinates to
6415!! u and v components relative to easterly and notherly direction. The
6416!! original fields are overwritten.
6417!! \todo Check and correct wind component flag (to be moved in
6418!! griddim_def?)
6419subroutine vg6d_wind_unrot(this)
6420type(volgrid6d) :: this !< object containing wind to be unrotated
6421
6422integer :: component_flag
6423
6424call get_val(this%griddim,component_flag=component_flag)
6425
6426if (component_flag == 1) then
6427 call l4f_category_log(this%category,l4f_info, &
6428 "unrotating vector components")
6429 call vg6d_wind__un_rot(this,.false.)
6430 call set_val(this%griddim,component_flag=0)
6431else
6432 call l4f_category_log(this%category,l4f_info, &
6433 "no need to unrotate vector components")
6434end if
6435
6436end subroutine vg6d_wind_unrot
6437
6438
6439!> Rotate the wind components.
6440!! It converts u and v components of vector quantities
6441!! relative to easterly and notherly direction to
6442!! defined grid in the direction of increasing x and y coordinates.
6443!! The original fields are overwritten.
6444subroutine vg6d_wind_rot(this)
6445type(volgrid6d) :: this !< object containing wind to be rotated
6446
6447integer :: component_flag
6448
6449call get_val(this%griddim,component_flag=component_flag)
6450
6451if (component_flag == 0) then
6452 call l4f_category_log(this%category,l4f_info, &
6453 "rotating vector components")
6454 call vg6d_wind__un_rot(this,.true.)
6455 call set_val(this%griddim,component_flag=1)
6456else
6457 call l4f_category_log(this%category,l4f_info, &
6458 "no need to rotate vector components")
6459end if
6460
6461end subroutine vg6d_wind_rot
6462
6463
6464! Generic UnRotate the wind components.
6465SUBROUTINE vg6d_wind__un_rot(this,rot)
6466TYPE(volgrid6d) :: this ! object containing wind to be (un)rotated
6467LOGICAL :: rot ! if .true. rotate else unrotate
6468
6469INTEGER :: i, j, k, l, a11, a12, a21, a22, stallo
6470double precision,pointer :: rot_mat(:,:,:)
6471real,allocatable :: tmp_arr(:,:)
6472REAL,POINTER :: voldatiu(:,:), voldativ(:,:)
6473INTEGER,POINTER :: iu(:), iv(:)
6474
6475IF (.NOT.ASSOCIATED(this%var)) THEN
6476 CALL l4f_category_log(this%category, l4f_error, &
6477 "trying to unrotate an incomplete volgrid6d object")
6478 CALL raise_fatal_error()
6479! RETURN
6480ENDIF
6481
6482CALL volgrid6d_var_hor_comp_index(this%var, iu, iv)
6483IF (.NOT.ASSOCIATED(iu)) THEN
6484 CALL l4f_category_log(this%category,l4f_error, &
6485 "unrotation impossible")
6486 CALL raise_fatal_error()
6487! RETURN
6488ENDIF
6489
6490! Temporary workspace
6491ALLOCATE(tmp_arr(this%griddim%dim%nx, this%griddim%dim%ny),stat=stallo)
6492IF (stallo /= 0) THEN
6493 CALL l4f_category_log(this%category, l4f_fatal, "allocating memory")
6494 CALL raise_fatal_error()
6495ENDIF
6496! allocate once for speed
6497IF (.NOT.ASSOCIATED(this%voldati)) THEN
6498 ALLOCATE(voldatiu(this%griddim%dim%nx, this%griddim%dim%ny), &
6499 voldativ(this%griddim%dim%nx, this%griddim%dim%ny))
6500ENDIF
6501
6502CALL griddim_unproj(this%griddim)
6503CALL wind_unrot(this%griddim, rot_mat)
6504
6505a11=1
6506if (rot)then
6507 a12=2
6508 a21=3
6509else
6510 a12=3
6511 a21=2
6512end if
6513a22=4
6514
6515DO l = 1, SIZE(iu)
6516 DO k = 1, SIZE(this%timerange)
6517 DO j = 1, SIZE(this%time)
6518 DO i = 1, SIZE(this%level)
6519! get data
6520 CALL volgrid_get_vol_2d(this, i, j, k, iu(l), voldatiu)
6521 CALL volgrid_get_vol_2d(this, i, j, k, iv(l), voldativ)
6522! convert units forward
6523! CALL compute(conv_fwd(iu(l)), voldatiu)
6524! CALL compute(conv_fwd(iv(l)), voldativ)
6525
6526! multiply wind components by rotation matrix
6527 WHERE(voldatiu /= rmiss .AND. voldativ /= rmiss)
6528 tmp_arr(:,:) = real(voldatiu(:,:)*rot_mat(:,:,a11) + &
6529 voldativ(:,:)*rot_mat(:,:,a12))
6530 voldativ(:,:) = real(voldatiu(:,:)*rot_mat(:,:,a21) + &
6531 voldativ(:,:)*rot_mat(:,:,a22))
6532 voldatiu(:,:) = tmp_arr(:,:)
6533 END WHERE
6534! convert units backward
6535! CALL uncompute(conv_fwd(iu(l)), voldatiu)
6536! CALL uncompute(conv_fwd(iv(l)), voldativ)
6537! put data
6538 CALL volgrid_set_vol_2d(this, i, j, k, iu(l), voldatiu)
6539 CALL volgrid_set_vol_2d(this, i, j, k, iv(l), voldativ)
6540 ENDDO
6541 ENDDO
6542 ENDDO
6543ENDDO
6544
6545IF (.NOT.ASSOCIATED(this%voldati)) THEN
6546 DEALLOCATE(voldatiu, voldativ)
6547ENDIF
6548DEALLOCATE(rot_mat, tmp_arr, iu, iv)
6549
6550END SUBROUTINE vg6d_wind__un_rot
6551
6552
6553!!$ try to understand the problem:
6554!!$
6555!!$ case:
6556!!$
6557!!$ 1) we have only one volume: we have to provide the direction of shift
6558!!$ compute H and traslate on it
6559!!$ 2) we have two volumes:
6560!!$ 1) volume U and volume V: compute H and traslate on it
6561!!$ 2) volume U/V and volume H : translate U/V on H
6562!!$ 3) we have tree volumes: translate U and V on H
6563!!$
6564!!$ strange cases:
6565!!$ 1) do not have U in volume U
6566!!$ 2) do not have V in volume V
6567!!$ 3) we have others variables more than U and V in volumes U e V
6568!!$
6569!!$
6570!!$ so the steps are:
6571!!$ 1) find the volumes
6572!!$ 2) define or compute H grid
6573!!$ 3) trasform the volumes in H
6574
6575!!$ N.B.
6576!!$ case 1) for only one vg6d (U or V) is not managed, but
6577!!$ the not pubblic subroutines will work but you have to know what you want to do
6578
6579
6580!> Convert grids type C to type A.
6581!! Use this to interpolate data from grid type C to grid type A
6582!! Grids type are defined by Arakawa.
6583!!
6584!! We need to find these types of area in a vg6d array
6585!! T area of points with temterature etc.
6586!! U area of points with u components of winds
6587!! V area of points with v components of winds
6588!!
6589!! this method works if it finds
6590!! two volumes:
6591!! 1) volume U and volume V: compute H and traslate on it
6592!! 2) volume U/V and volume H : translate U/V on H
6593!! three volumes: translate U and V on H
6594!!
6595!! try to work well on more datasets at once
6596subroutine vg6d_c2a (this)
6597
6598TYPE(volgrid6d),INTENT(inout) :: this(:) !< vettor of volumes volgrid6d to elaborate
6599
6600integer :: ngrid,igrid,jgrid,ugrid,vgrid,tgrid
6601doubleprecision :: xmin, xmax, ymin, ymax
6602doubleprecision :: xmin_t, xmax_t, ymin_t, ymax_t
6603doubleprecision :: step_lon_t,step_lat_t
6604character(len=80) :: type_t,type
6605TYPE(griddim_def):: griddim_t
6606
6607ngrid=size(this)
6608
6609do igrid=1,ngrid
6610
6611 call init(griddim_t)
6612
6613 call get_val(this(igrid)%griddim,xmin=xmin_t, xmax=xmax_t, ymin=ymin_t, ymax=ymax_t,proj_type=type_t)
6614 step_lon_t=(xmax_t-xmin_t)/dble(this(igrid)%griddim%dim%nx-1)
6615 step_lat_t=(ymax_t-ymin_t)/dble(this(igrid)%griddim%dim%ny-1)
6616
6617 do jgrid=1,ngrid
6618
6619 ugrid=imiss
6620 vgrid=imiss
6621 tgrid=imiss
6622
6623#ifdef DEBUG
6624 call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: search U/V/T points:"//to_char(igrid)//to_char(jgrid))
6625 call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: test delta: "//to_char(step_lon_t)//to_char(step_lat_t))
6626#endif
6627
6628 if (this(igrid)%griddim == this(jgrid)%griddim ) cycle
6629
6630 if (this(igrid)%griddim%dim%nx == this(jgrid)%griddim%dim%nx .and. &
6631 this(igrid)%griddim%dim%ny == this(jgrid)%griddim%dim%ny ) then
6632
6633 call get_val(this(jgrid)%griddim,xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,proj_type=type)
6634
6635 if (type_t /= type )cycle
6636
6637#ifdef DEBUG
6638 call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: test U "//&
6639 to_char(xmin)//to_char(xmax)//to_char(ymin)//to_char(ymax))
6640
6641 call l4f_category_log(this(igrid)%category,l4f_debug,"diff coordinate lon"//&
6642 to_char(abs(xmin - (xmin_t+step_lon_t/2.d0)))//&
6643 to_char(abs(xmax - (xmax_t+step_lon_t/2.d0))))
6644 call l4f_category_log(this(igrid)%category,l4f_debug,"diff coordinate lat"//&
6645 to_char(abs(ymin - (ymin_t+step_lat_t/2.d0)))//&
6646 to_char(abs(ymax - (ymax_t+step_lat_t/2.d0))))
6647#endif
6648
6649 if ( abs(xmin - (xmin_t+step_lon_t/2.d0)) < 1.d-3 .and. abs(xmax - (xmax_t+step_lon_t/2.d0)) < 1.d-3 ) then
6650 if ( abs(ymin - ymin_t) < 1.d-3 .and. abs(ymax - ymax_t) < 1.d-3 ) then
6651
6652#ifdef DEBUG
6653 call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: found U")
6654#endif
6655 ugrid=jgrid
6656 tgrid=igrid
6657
6658 end if
6659 end if
6660
6661#ifdef DEBUG
6662 call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: test V "//&
6663 to_char(xmin)//to_char(xmax)//to_char(ymin)//to_char(ymax))
6664#endif
6665
6666 if ( abs(ymin - (ymin_t+step_lat_t/2.d0)) < 1.d-3 .and. abs(ymax - (ymax_t+step_lat_t/2.d0)) < 1.d-3 ) then
6667 if ( abs(xmin - xmin_t) < 1.d-3 .and. abs(xmax - xmax_t) < 1.d-3 ) then
6668
6669#ifdef DEBUG
6670 call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: found V")
6671#endif
6672 vgrid=jgrid
6673 tgrid=igrid
6674
6675 end if
6676 end if
6677
6678
6679 ! test if we have U and V
6680
6681#ifdef DEBUG
6682 call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: test U and V"//&
6683 to_char(xmin_t)//to_char(xmax_t)//to_char(ymin_t)//to_char(ymax_t)//&
6684 to_char(xmin)//to_char(xmax)//to_char(ymin)//to_char(ymax))
6685
6686 call l4f_category_log(this(igrid)%category,l4f_debug,"UV diff coordinate lon"//&
6687 to_char(abs(xmin_t - xmin)-step_lon_t/2.d0)//&
6688 to_char(abs(xmax_t - xmax)-step_lon_t/2.d0))
6689 call l4f_category_log(this(igrid)%category,l4f_debug,"UV diff coordinate lat"//&
6690 to_char(abs(ymin_t - ymin) -step_lat_t/2.d0)//&
6691 to_char(abs(ymax_t - ymax)-step_lat_t/2.d0))
6692#endif
6693
6694 if ( abs(ymin - (ymin_t+step_lat_t/2.d0)) < 2.d-3 .and. abs(ymax - (ymax_t+step_lat_t/2.d0)) < 2.d-3 ) then
6695 if ( abs(xmin_t - (xmin+step_lon_t/2.d0)) < 2.d-3 .and. abs(xmax_t - (xmax+step_lon_t/2.d0)) < 2.d-3 ) then
6696
6697#ifdef DEBUG
6698 call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: found U and V case up and right")
6699#endif
6700
6701 vgrid=jgrid
6702 ugrid=igrid
6703
6704 call init(griddim_t,xmin=xmin, xmax=xmax, ymin=ymin_t, ymax=ymax_t)
6705
6706 end if
6707 end if
6708 end if
6709
6710 ! abbiamo almeno due volumi: riportiamo U e/o V su T
6711 if (c_e(ugrid)) then
6712#ifdef DEBUG
6713 call l4f_category_log(this(igrid)%category,l4f_debug,"C grid U points: interpolate U on T "//to_char(tgrid)//to_char(ugrid))
6714#endif
6715 if (c_e(tgrid))then
6716 call vg6d_c2a_grid(this(ugrid),this(tgrid)%griddim,cgrid=1)
6717 else
6718 call vg6d_c2a_grid(this(ugrid),griddim_t,cgrid=1)
6719 end if
6720 call vg6d_c2a_mat(this(ugrid),cgrid=1)
6721 end if
6722
6723 if (c_e(vgrid)) then
6724#ifdef DEBUG
6725 call l4f_category_log(this(igrid)%category,l4f_debug,"C grid V points: interpolate V on T "//to_char(tgrid)//to_char(vgrid))
6726#endif
6727 if (c_e(tgrid))then
6728 call vg6d_c2a_grid(this(vgrid),this(tgrid)%griddim,cgrid=2)
6729 else
6730 call vg6d_c2a_grid(this(vgrid),griddim_t,cgrid=2)
6731 end if
6732 call vg6d_c2a_mat(this(vgrid),cgrid=2)
6733 end if
6734
6735 end do
6736
6737 call delete(griddim_t)
6738
6739end do
6740
6741
6742end subroutine vg6d_c2a
6743
6744
6745! Convert C grid to A grid
6746subroutine vg6d_c2a_grid(this,griddim_t,cgrid)
6747
6748type(volgrid6d),intent(inout) :: this ! object containing fields to be translated (U or V points)
6749type(griddim_def),intent(in),optional :: griddim_t ! object containing grid of T points
6750integer,intent(in) :: cgrid ! in C grid (Arakawa) we have 0=T,1=U,2=V points
6751
6752doubleprecision :: xmin, xmax, ymin, ymax
6753doubleprecision :: step_lon,step_lat
6754
6755
6756if (present(griddim_t)) then
6757
6758 call get_val(griddim_t,xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
6759 call set_val(this%griddim,xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
6760! improve grid definition precision
6761 CALL griddim_setsteps(this%griddim)
6762
6763else
6764
6765 select case (cgrid)
6766
6767 case(0)
6768
6769#ifdef DEBUG
6770 call l4f_category_log(this%category,l4f_debug,"C grid: T points, nothing to do")
6771#endif
6772 return
6773
6774 case (1)
6775
6776#ifdef DEBUG
6777 call l4f_category_log(this%category,l4f_debug,"C grid: U points, we need interpolation")
6778#endif
6779
6780 call get_val(this%griddim, xmin=xmin, xmax=xmax)
6781 step_lon=(xmax-xmin)/dble(this%griddim%dim%nx-1)
6782 xmin=xmin-step_lon/2.d0
6783 xmax=xmax-step_lon/2.d0
6784 call set_val(this%griddim, xmin=xmin, xmax=xmax)
6785! improve grid definition precision
6786 CALL griddim_setsteps(this%griddim)
6787
6788 case (2)
6789
6790#ifdef DEBUG
6791 call l4f_category_log(this%category,l4f_debug,"C grid: V points, we need interpolation")
6792#endif
6793
6794 call get_val(this%griddim, ymin=ymin, ymax=ymax)
6795 step_lat=(ymax-ymin)/dble(this%griddim%dim%ny-1)
6796 ymin=ymin-step_lat/2.d0
6797 ymax=ymax-step_lat/2.d0
6798 call set_val(this%griddim, ymin=ymin, ymax=ymax)
6799! improve grid definition precision
6800 CALL griddim_setsteps(this%griddim)
6801
6802 case default
6803
6804 call l4f_category_log(this%category,l4f_fatal,"C grid type not known")
6805 call raise_fatal_error ()
6806
6807 end select
6808
6809end if
6810
6811
6812call griddim_unproj(this%griddim)
6813
6814
6815end subroutine vg6d_c2a_grid
6816
6817! Convert C grid to A grid
6818subroutine vg6d_c2a_mat(this,cgrid)
6819
6820type(volgrid6d),intent(inout) :: this ! object containing fields to be translated
6821integer,intent(in) :: cgrid ! in C grid (Arakawa) we have 0=T,1=U,2=V points
6822
6823INTEGER :: i, j, k, iv, stallo
6824REAL,ALLOCATABLE :: tmp_arr(:,:)
6825REAL,POINTER :: voldatiuv(:,:)
6826
6827
6828IF (cgrid == 0) RETURN ! nothing to do
6829IF (cgrid /= 1 .AND. cgrid /= 2) THEN ! wrong cgrid
6830 CALL l4f_category_log(this%category,l4f_fatal,"C grid type "// &
6831 trim(to_char(cgrid))//" not known")
6832 CALL raise_error()
6833 RETURN
6834ENDIF
6835
6836! Temporary workspace
6837ALLOCATE(tmp_arr(this%griddim%dim%nx, this%griddim%dim%ny),stat=stallo)
6838if (stallo /=0)then
6839 call l4f_log(l4f_fatal,"allocating memory")
6840 call raise_fatal_error()
6841end if
6842
6843! allocate once for speed
6844IF (.NOT.ASSOCIATED(this%voldati)) THEN
6845 ALLOCATE(voldatiuv(this%griddim%dim%nx, this%griddim%dim%ny), stat=stallo)
6846 IF (stallo /= 0) THEN
6847 CALL l4f_log(l4f_fatal,"allocating memory")
6848 CALL raise_fatal_error()
6849 ENDIF
6850ENDIF
6851
6852IF (cgrid == 1) THEN ! U points to H points
6853 DO iv = 1, SIZE(this%var)
6854 DO k = 1, SIZE(this%timerange)
6855 DO j = 1, SIZE(this%time)
6856 DO i = 1, SIZE(this%level)
6857 tmp_arr(:,:) = rmiss
6858 CALL volgrid_get_vol_2d(this, i, j, k, iv, voldatiuv)
6859
6860! West boundary extrapolation
6861 WHERE(voldatiuv(1,:) /= rmiss .AND. voldatiuv(2,:) /= rmiss)
6862 tmp_arr(1,:) = voldatiuv(1,:) - (voldatiuv(2,:) - voldatiuv(1,:)) / 2.
6863 END WHERE
6864
6865! Rest of the field
6866 WHERE(voldatiuv(1:this%griddim%dim%nx-1,:) /= rmiss .AND. &
6867 voldatiuv(2:this%griddim%dim%nx,:) /= rmiss)
6868 tmp_arr(2:this%griddim%dim%nx,:) = &
6869 (voldatiuv(1:this%griddim%dim%nx-1,:) + &
6870 voldatiuv(2:this%griddim%dim%nx,:)) / 2.
6871 END WHERE
6872
6873 voldatiuv(:,:) = tmp_arr
6874 CALL volgrid_set_vol_2d(this, i, j, k, iv, voldatiuv)
6875 ENDDO
6876 ENDDO
6877 ENDDO
6878 ENDDO
6879
6880ELSE IF (cgrid == 2) THEN ! V points to H points
6881 DO iv = 1, SIZE(this%var)
6882 DO k = 1, SIZE(this%timerange)
6883 DO j = 1, SIZE(this%time)
6884 DO i = 1, SIZE(this%level)
6885 tmp_arr(:,:) = rmiss
6886 CALL volgrid_get_vol_2d(this, i, j, k, iv, voldatiuv)
6887
6888! South boundary extrapolation
6889 WHERE(voldatiuv(:,1) /= rmiss .AND. voldatiuv(:,2) /= rmiss)
6890 tmp_arr(:,1) = voldatiuv(:,1) - (voldatiuv(:,2) - voldatiuv(:,1)) / 2.
6891 END WHERE
6892
6893! Rest of the field
6894 WHERE(voldatiuv(:,1:this%griddim%dim%ny-1) /= rmiss .AND. &
6895 voldatiuv(:,2:this%griddim%dim%ny) /= rmiss)
6896 tmp_arr(:,2:this%griddim%dim%ny) = &
6897 (voldatiuv(:,1:this%griddim%dim%ny-1) + &
6898 voldatiuv(:,2:this%griddim%dim%ny)) / 2.
6899 END WHERE
6900
6901 voldatiuv(:,:) = tmp_arr
6902 CALL volgrid_set_vol_2d(this, i, j, k, iv, voldatiuv)
6903 ENDDO
6904 ENDDO
6905 ENDDO
6906 ENDDO
6907ENDIF ! tertium non datur
6908
6909IF (.NOT.ASSOCIATED(this%voldati)) THEN
6910 DEALLOCATE(voldatiuv)
6911ENDIF
6912DEALLOCATE (tmp_arr)
6913
6914end subroutine vg6d_c2a_mat
6915
6916
6917!> \brief Display object on screen
6918!!
6919!! Show brief content on screen.
6920subroutine display_volgrid6d (this)
6921
6922TYPE(volgrid6d),intent(in) :: this !< object to display
6923integer :: i
6924
6925#ifdef DEBUG
6926call l4f_category_log(this%category,l4f_debug,"ora mostro gridinfo " )
6927#endif
6928
6929print*,"----------------------- volgrid6d display ---------------------"
6930call display(this%griddim)
6931
6932IF (ASSOCIATED(this%time))then
6933 print*,"---- time vector ----"
6934 print*,"elements=",size(this%time)
6935 do i=1, size(this%time)
6936 call display(this%time(i))
6937 end do
6938end IF
6939
6940IF (ASSOCIATED(this%timerange))then
6941 print*,"---- timerange vector ----"
6942 print*,"elements=",size(this%timerange)
6943 do i=1, size(this%timerange)
6944 call display(this%timerange(i))
6945 end do
6946end IF
6947
6948IF (ASSOCIATED(this%level))then
6949 print*,"---- level vector ----"
6950 print*,"elements=",size(this%level)
6951 do i=1, size(this%level)
6952 call display(this%level(i))
6953 end do
6954end IF
6955
6956IF (ASSOCIATED(this%var))then
6957 print*,"---- var vector ----"
6958 print*,"elements=",size(this%var)
6959 do i=1, size(this%var)
6960 call display(this%var(i))
6961 end do
6962end IF
6963
6964IF (ASSOCIATED(this%gaid))then
6965 print*,"---- gaid vector (present mask only) ----"
6966 print*,"elements=",shape(this%gaid)
6967 print*,c_e(reshape(this%gaid,(/SIZE(this%gaid)/)))
6968end IF
6969
6970print*,"--------------------------------------------------------------"
6971
6972
6973end subroutine display_volgrid6d
6974
6975
6976!> \brief Display vector of object on screen
6977!!
6978!! Show brief content on screen.
6979subroutine display_volgrid6dv (this)
6980
6981TYPE(volgrid6d),intent(in) :: this(:) !< vector of object to display
6982integer :: i
6983
6984print*,"----------------------- volgrid6d vector ---------------------"
6985
6986print*,"elements=",size(this)
6987
6988do i=1, size(this)
6989
6990#ifdef DEBUG
6991 call l4f_category_log(this(i)%category,l4f_debug,"ora mostro il vettore volgrid6d" )
6992#endif
6993
6994 call display(this(i))
6995
6996end do
6997print*,"--------------------------------------------------------------"
6998
6999end subroutine display_volgrid6dv
7000
7001
7002!> Reduce some dimensions (level and timerage) for semplification (rounding).
7003!! Vector version of vg6dv_rounding
7004subroutine vg6dv_rounding(vg6din,vg6dout,level,timerange,nostatproc,merge)
7005type(volgrid6d),intent(in) :: vg6din(:) !< input volume
7006type(volgrid6d),intent(out),pointer :: vg6dout(:) !> output volume
7007type(vol7d_level),intent(in),optional :: level(:) !< almost equal level list
7008type(vol7d_timerange),intent(in),optional :: timerange(:) !< almost equal timerange list
7009logical,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
7010logical,intent(in),optional :: nostatproc !< do not take in account statistical processing code in timerange and P2
7011
7012integer :: i
7013
7014allocate(vg6dout(size(vg6din)))
7015
7016do i = 1, size(vg6din)
7017 call vg6d_rounding(vg6din(i),vg6dout(i),level,timerange,nostatproc,merge)
7018end do
7019
7020end subroutine vg6dv_rounding
7021
7022!> Reduce some dimensions (level and timerage) for semplification (rounding).
7023!! You can use this for simplify and use variables in computation like alchimia
7024!! where fields have to be on the same coordinate
7025!! examples:
7026!! means in time for short periods and istantaneous values
7027!! 2 meter and surface levels
7028!! If there are data on more then one almost equal levels or timeranges, the first var present (at least one point)
7029!! will be taken (order is by icreasing var index).
7030!! You can use predefined values for classic semplification
7031!! almost_equal_levels and almost_equal_timeranges
7032!! The level or timerange in output will be defined by the first element of level and timerange list
7033subroutine vg6d_rounding(vg6din,vg6dout,level,timerange,nostatproc,merge)
7034type(volgrid6d),intent(in) :: vg6din !< input volume
7035type(volgrid6d),intent(out) :: vg6dout !> output volume
7036type(vol7d_level),intent(in),optional :: level(:) !< almost equal level list
7037type(vol7d_timerange),intent(in),optional :: timerange(:) !< almost equal timerange list
7038logical,intent(in),optional :: merge !< if there are data on more then one almost equal levels or timeranges
7039!! will be merged POINT BY POINT with priority for the fird data found ordered by icreasing var index
7040logical,intent(in),optional :: nostatproc !< do not take in account statistical processing code in timerange and P2
7041
7042integer :: ilevel,itimerange
7043type(vol7d_level) :: roundlevel(size(vg6din%level))
7044type(vol7d_timerange) :: roundtimerange(size(vg6din%timerange))
7045
7046roundlevel=vg6din%level
7047
7048if (present(level))then
7049 do ilevel = 1, size(vg6din%level)
7050 if ((any(vg6din%level(ilevel) .almosteq. level))) then
7051 roundlevel(ilevel)=level(1)
7052 end if
7053 end do
7054end if
7055
7056roundtimerange=vg6din%timerange
7057
7058if (present(timerange))then
7059 do itimerange = 1, size(vg6din%timerange)
7060 if ((any(vg6din%timerange(itimerange) .almosteq. timerange))) then
7061 roundtimerange(itimerange)=timerange(1)
7062 end if
7063 end do
7064end if
7065
7066!set istantaneous values everywere
7067!preserve p1 for forecast time
7068if (optio_log(nostatproc)) then
7069 roundtimerange(:)%timerange=254
7070 roundtimerange(:)%p2=0
7071end if
7072
7073
7074call vg6d_reduce(vg6din,vg6dout,roundlevel,roundtimerange,merge)
7075
7076end subroutine vg6d_rounding
7077
7078!> Reduce some dimensions (level and timerage).
7079!! You can pass a volume and specify duplicated levels and timeranges in roundlevel and roundtimerange;
7080!! you get unique levels and timeranges in output.
7081!! If there are data on equal levels or timeranges, the first var present (at least one point)
7082!! will be taken (order is by icreasing var index).
7083!! you can specify merge and if there are data on equal levels or timeranges will be merged POINT BY POINT
7084!! with priority for the first data found ordered by icreasing var index (require to decode all the data)
7085!! Data are decoded only if needed so the output should be with or without voldata allocated
7086subroutine vg6d_reduce(vg6din,vg6dout,roundlevel,roundtimerange,merge)
7087type(volgrid6d),intent(in) :: vg6din !< input volume
7088type(volgrid6d),intent(out) :: vg6dout !< output volume
7089type(vol7d_level),intent(in) :: roundlevel(:) !< new level list to use for rounding
7090type(vol7d_timerange),intent(in) :: roundtimerange(:) !< new timerange list to use for rounding
7091logical,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)
7092
7093integer :: nlevel,ntime,ntimerange,nvar,ilevel,itimerange,ivar,indl,indt,itime,nx,ny
7094real,allocatable :: vol2d(:,:)
7095
7096nx=vg6din%griddim%dim%nx
7097ny=vg6din%griddim%dim%ny
7098nlevel=count_distinct(roundlevel,back=.true.)
7099ntime=size(vg6din%time)
7100ntimerange=count_distinct(roundtimerange,back=.true.)
7101nvar=size(vg6din%var)
7102
7103call init(vg6dout, vg6din%griddim, vg6din%time_definition, categoryappend="generated by vg6d_reduce")
7104call volgrid6d_alloc(vg6dout, vg6din%griddim%dim, ntime, nlevel, ntimerange, nvar)
7105
7106if ( ASSOCIATED(vg6din%voldati) .or. optio_log(merge)) then
7107 call volgrid6d_alloc_vol(vg6dout,inivol=.true.,decode=.true.)
7108 allocate(vol2d(nx,ny))
7109else
7110 call volgrid6d_alloc_vol(vg6dout,inivol=.true.,decode=.false.)
7111end if
7112
7113vg6dout%time=vg6din%time
7114vg6dout%var=vg6din%var
7115vg6dout%timerange=pack_distinct(roundtimerange,ntimerange,back=.true.)
7116vg6dout%level=pack_distinct(roundlevel,nlevel,back=.true.)
7117! sort modified dimensions
7118CALL sort(vg6dout%timerange)
7119CALL sort(vg6dout%level)
7120
7121do ilevel=1,size(vg6din%level)
7122 indl=index(vg6dout%level,roundlevel(ilevel))
7123 do itimerange=1,size(vg6din%timerange)
7124 indt=index(vg6dout%timerange,roundtimerange(itimerange))
7125 do ivar=1, nvar
7126 do itime=1,ntime
7127
7128 if ( ASSOCIATED(vg6din%voldati)) then
7129 vol2d=vg6din%voldati(:,:,ilevel,itime,itimerange,ivar)
7130 end if
7131
7132 if (optio_log(merge)) then
7133
7134 if ( .not. ASSOCIATED(vg6din%voldati)) then
7135 CALL grid_id_decode_data(vg6din%gaid(ilevel,itime,itimerange,ivar), vol2d)
7136 end if
7137
7138 !! merge present data point by point
7139 where (.not. c_e(vg6dout%voldati(:,:,indl,itime,indt,ivar)))
7140
7141 vg6dout%voldati(:,:,indl,itime,indt,ivar)=vol2d
7142
7143 end where
7144 else if ( ASSOCIATED(vg6din%voldati)) then
7145 if (.not. any(c_e(vg6dout%voldati(:,:,indl,itime,indt,ivar))))then
7146 vg6dout%voldati(:,:,indl,itime,indt,ivar)=vol2d
7147 end if
7148 end if
7149
7150 if (c_e(vg6din%gaid(ilevel,itime,itimerange,ivar)).and. .not. c_e(vg6dout%gaid(indl,itime,indt,ivar)))then
7151 call copy (vg6din%gaid(ilevel,itime,itimerange,ivar), vg6dout%gaid(indl,itime,indt,ivar))
7152 end if
7153 end do
7154 end do
7155 end do
7156end do
7157
7158if ( ASSOCIATED(vg6din%voldati) .or. optio_log(merge)) then
7159 deallocate(vol2d)
7160end if
7161
7162end subroutine vg6d_reduce
7163
7164
7165end module volgrid6d_class
7166
7167
7168
7169!>\example example_vg6d_3.f90
7170!!\brief Programma esempio semplice per gridinfo e volgrid6d.
7171!!
7172!! 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.
7173
7174!>\example example_vg6d_5.f90
7175!!\brief Programma trasformazione da volgrid6d a volgrid6d
7176!!
7177!! Legge grib da un file e li organizza in un vettore di strutture
7178!! volgrid6d mettendoli a disposizione per eventuali elaborazioni;
7179!! vengono poi riesportati a un file grib
7180
7181!>\example example_vg6d_8.f90
7182!! \brief Programma scrittura su file vettore di anagrafica
7183
7184!>\example example_vg6d_6.f90
7185!! \brief Programma trasformazione da volgrid6d a vol7d
7186
7187!>\example example_vg6d_7.f90
7188!! \brief Programma trasformazione da vol7d a volgrid6d
7189
7190!>\example example_vg6d_9.f90
7191!! \brief Example to create a grib editionNumber = 2 file from data generated in memory using a grib_api template
Method for inserting elements of the array at a desired position.
Operatore di valore assoluto di un intervallo.
Legge un oggetto datetime/timedelta o un vettore di oggetti datetime/timedelta da un file FORMATTED o...
Functions that return a trimmed CHARACTER representation of the input variable.
Restituiscono il valore dell'oggetto in forma di stringa stampabile.
Scrive un oggetto datetime/timedelta o un vettore di oggetti datetime/timedelta su un file FORMATTED ...
Copy an object, creating a fully new instance.
Method for returning the contents of the object.
Compute forward coordinate transformation from geographical system to projected system.
Method for setting 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.
Convert a level type to a physical variable.
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.
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.
Object describing a rectangular, homogeneous gridded dataset.

Generated with Doxygen.