libsim Versione 7.2.6

◆ vg6d_reduce()

subroutine, public vg6d_reduce ( type(volgrid6d), intent(in) vg6din,
type(volgrid6d), intent(out) vg6dout,
type(vol7d_level), dimension(:), intent(in) roundlevel,
type(vol7d_timerange), dimension(:), intent(in) roundtimerange,
logical, intent(in), optional merge )

Reduce some dimensions (level and timerage).

You can pass a volume and specify duplicated levels and timeranges in roundlevel and roundtimerange; you get unique levels and timeranges in output. If there are data on equal levels or timeranges, the first var present (at least one point) will be taken (order is by icreasing var index). you can specify merge and 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) Data are decoded only if needed so the output should be with or without voldata allocated

Parametri
[in]vg6dininput volume
[out]vg6doutoutput volume
[in]roundlevelnew level list to use for rounding
[in]roundtimerangenew timerange list to use for rounding
[in]mergeif 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)

Definizione alla linea 3653 del file volgrid6d_class.F90.

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