display_volgrid6d display_volgrid6dv export_to_gridinfo export_to_gridinfov export_to_gridinfovv import_from_gridinfo import_from_gridinfovv v7d_v7d_transform v7d_volgrid6d_transform vg6d_c2a vg6d_reduce vg6d_rounding vg6d_wind_rot vg6d_wind_unrot vg6dv_rounding volgrid6d_alloc volgrid6d_alloc_vol volgrid6d_delete volgrid6d_export_to_file volgrid6d_import_from_file volgrid6d_init volgrid6d_read_from_file volgrid6d_transform volgrid6d_v7d_transform volgrid6d_v7d_transform_compute volgrid6d_write_on_file volgrid6dv_delete volgrid6dv_transform volgrid6dv_v7d_transform volgrid_get_vol_2d volgrid_get_vol_3d volgrid_set_vol_2d volgrid_set_vol_3d

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] vg6din input volume [out] vg6dout output volume [in] roundlevel new level list to use for rounding [in] roundtimerange new timerange list to use for rounding [in] 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) Definizione alla linea 3663 del file volgrid6d_class.F90. ! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it> ! authors: ! Davide Cesari <dcesari@arpa.emr.it> ! Paolo Patruno <ppatruno@arpa.emr.it> ! This program is free software; you can redistribute it and/or ! modify it under the terms of the GNU General Public License as ! published by the Free Software Foundation; either version 2 of ! the License, or (at your option) any later version. ! This program is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU General Public License for more details. ! You should have received a copy of the GNU General Public License ! along with this program. If not, see <http://www.gnu.org/licenses/>. #include "config.h" !> \defgroup volgrid6d Libsim package, volgrid6d library. !! The libsim volgrid6d library contains classes for managing data on !! regular rectangular grids, tipically the output of a numerical !! weather prediction model, and for their import from a WMO GRIB file !! of from other geophysical data file formats. In order to compile !! and link programs using this library, you have to insert the !! required \c USE statements in the program units involved, specify !! the location of module files when compiling (tipically \c !! -I/usr/lib/gfortran/modules or \c -I/usr/lib64/gfortran/modules or !! \c -I/usr/include) and indicate the library name \c -lsim_volgrid6d !! when linking, assuming that the library has been installed in a !! default location. !> This module defines objects and methods for managing data volumes !! on rectangular georeferenced grids. The data are accomodated in a !! multi-dimensional array with 6 predefined dimensions. Different !! geographic coordinates and projections are supported, mainly !! inspired by grib coding standards. The \a volgrid6d object contains !! information and data on an homogeneous grid definition, while !! different grids are managed as arrays of \a volgrid6d objects. !! Every object contains also an identificator of the grid (\a grid_id !! object), carrying information about the driver used or which has to !! be used for import/export from/to file. With the help of \a !! gridinfo_def class, data can be imported and exported to the !! supported formats, mainly grib1 and grib2 through grib_api and many !! GIS-style formats through gdal. !! !! Simple example program \include example_vg6d_3.f90 !! Example of transformation from volgrid6d to volgrid6d \include example_vg6d_5.f90 !! Example of transformation from volgrid6d to vol7d \include example_vg6d_6.f90 !! Example of transformation from vol7d to volgrid6d \include example_vg6d_7.f90 !! !!\ingroup volgrid6d MODULE volgrid6d_class USE vol7d_level_class USE vol7d_timerange_class USE vol7d_class USE grid_transform_class USE geo_proj_class USE grid_class USE grid_dim_class USE datetime_class USE volgrid6d_var_class USE log4fortran USE array_utilities USE grid_id_class USE gridinfo_class USE optional_values !USE file_utilities #ifdef HAVE_DBALLE USE vol7d_dballe_class #endif IMPLICIT NONE character (len=255),parameter:: subcategory="volgrid6d_class" !> Object describing a rectangular, homogeneous gridded dataset type volgrid6d type(griddim_def) :: griddim !< grid descriptor TYPE(datetime),pointer :: time(:) !< time dimension descriptor TYPE(vol7d_timerange),pointer :: timerange(:) !< timerange (forecast, analysis, statistically processed) dimension descriptor TYPE(vol7d_level),pointer :: level(:) !< vertical level dimension descriptor TYPE(volgrid6d_var),pointer :: var(:) !< physical variable dimension descriptor 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) REAL,POINTER :: voldati(:,:,:,:,:,:) !< array of data, indices are: (x,y,level,time,timerange,var) INTEGER :: time_definition !< time definition; 0=time is reference time ; 1=time is validity time TYPE(l4f_handle) :: category !< log4fortran category end type volgrid6d !> Constructor, it creates a new instance of the object. INTERFACE init MODULE PROCEDURE volgrid6d_init END INTERFACE !> Destructor, it releases every information and memory buffer !! associated with the object. INTERFACE delete MODULE PROCEDURE volgrid6d_delete, volgrid6dv_delete END INTERFACE !> Import an object dirctly from a native file, from a \a gridinfo object !! or from a supported file format through a \a gridinfo object. INTERFACE import MODULE PROCEDURE volgrid6d_read_from_file MODULE PROCEDURE import_from_gridinfo, import_from_gridinfovv, & volgrid6d_import_from_file END INTERFACE !> Export an object dirctly to a native file, to a \a gridinfo object !! or to a supported file format through a \a gridinfo object. INTERFACE export MODULE PROCEDURE volgrid6d_write_on_file MODULE PROCEDURE export_to_gridinfo, export_to_gridinfov, export_to_gridinfovv,& volgrid6d_export_to_file END INTERFACE ! methods for computing transformations through an initialised ! grid_transform object, probably too low level to be interfaced INTERFACE compute MODULE PROCEDURE volgrid6d_transform_compute, volgrid6d_v7d_transform_compute,& v7d_volgrid6d_transform_compute, v7d_v7d_transform_compute END INTERFACE !> Transform between any combination of \a volgrid6d and \a vol7d objects !! by means of a \a transform_def object describing the transformation. INTERFACE transform MODULE PROCEDURE volgrid6d_transform, volgrid6dv_transform,& volgrid6d_v7d_transform, volgrid6dv_v7d_transform, v7d_volgrid6d_transform, & v7d_v7d_transform END INTERFACE INTERFACE wind_rot MODULE PROCEDURE vg6d_wind_rot END INTERFACE INTERFACE wind_unrot MODULE PROCEDURE vg6d_wind_unrot END INTERFACE !> Display on standard output a description of the \a volgrid6d object !! provided. INTERFACE display MODULE PROCEDURE display_volgrid6d,display_volgrid6dv END INTERFACE !> Reduce some dimensions (level and timerage) for semplification (rounding). !! You can use this for simplify and use variables in computation like alchimia !! where fields have to be on the same coordinate !! examples: !! means in time for short periods and istantaneous values !! 2 meter and surface levels !! If there are data on more then one almost equal levels or timeranges, the first var present (at least one point) !! will be taken (order is by icreasing var index). !! You can use predefined values for classic semplification !! almost_equal_levels and almost_equal_timeranges !! The level or timerange in output will be defined by the first element of level and timerange list INTERFACE rounding MODULE PROCEDURE vg6d_rounding, vg6dv_rounding END INTERFACE private PUBLIC volgrid6d,init,delete,export,import,compute,transform, & wind_rot,wind_unrot,vg6d_c2a,display,volgrid6d_alloc,volgrid6d_alloc_vol, & volgrid_get_vol_2d, volgrid_set_vol_2d, volgrid_get_vol_3d, volgrid_set_vol_3d PUBLIC rounding, vg6d_reduce CONTAINS !> Constructor, it creates a new instance of the object. !! The constructor should be explicitly used only in rare cases, !! \a volgrid6d objects are usually created through the \a import !! interface. SUBROUTINE volgrid6d_init(this, griddim, time_definition, categoryappend) TYPE(volgrid6d) :: this !< object to be initialized TYPE(griddim_def),OPTIONAL :: griddim !< grid descriptor INTEGER,INTENT(IN),OPTIONAL :: time_definition !< 0=time is reference time; 1=time is validity time CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category character(len=512) :: a_name if (present(categoryappend))then call l4f_launcher(a_name,a_name_append=trim(subcategory)//"."//trim(categoryappend)) else call l4f_launcher(a_name,a_name_append=trim(subcategory)) endif this%category=l4f_category_get_handle(a_name) #ifdef DEBUG call l4f_category_log(this%category,l4f_debug,"init") #endif call init(this%griddim) if (present(griddim))then call copy (griddim,this%griddim) end if if(present(time_definition)) then this%time_definition = time_definition else this%time_definition = 0 !default to reference time end if nullify (this%time,this%timerange,this%level,this%var) nullify (this%gaid,this%voldati) END SUBROUTINE volgrid6d_init !> Allocate the dimension descriptors of the \a volgrid6d object. !! This method allocates the horizontal grid descriptor and the one !! dimensional arrays of the dimensions !! - time !! - vertical level !! - timerange !! - physical variable !! !! This method should be explicitly used only in rare cases, it is !! usually called implicitly through the \a import interface. SUBROUTINE volgrid6d_alloc(this, dim, ntime, nlevel, ntimerange, nvar, ini) TYPE(volgrid6d),INTENT(inout) :: this !< object whose decriptors should be allocated TYPE(grid_dim),INTENT(in),OPTIONAL :: dim !< horizontal grid size X, Y INTEGER,INTENT(in),OPTIONAL :: ntime !< number of time levels INTEGER,INTENT(in),OPTIONAL :: nlevel !< number of vertical levels INTEGER,INTENT(in),OPTIONAL :: ntimerange !< number of different timeranges INTEGER,INTENT(in),OPTIONAL :: nvar !< number of physical variables LOGICAL,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 INTEGER :: i, stallo LOGICAL :: linit #ifdef DEBUG call l4f_category_log(this%category,l4f_debug,"alloc") #endif IF (PRESENT(ini)) THEN linit = ini ELSE linit = .false. ENDIF if (present(dim)) call copy (dim,this%griddim%dim) IF (PRESENT(ntime)) THEN IF (ntime >= 0) THEN IF (ASSOCIATED(this%time)) DEALLOCATE(this%time) #ifdef DEBUG call l4f_category_log(this%category,l4f_debug,"alloc ntime "//to_char(ntime)) #endif ALLOCATE(this%time(ntime),stat=stallo) if (stallo /=0)then call l4f_category_log(this%category,l4f_fatal,"allocating memory") CALL raise_fatal_error() end if IF (linit) THEN DO i = 1, ntime this%time(i) = datetime_miss ! CALL init(this%time(i)) ! senza argomento inizializza a zero non missing ! baco di datetime? ENDDO ENDIF ENDIF ENDIF IF (PRESENT(nlevel)) THEN IF (nlevel >= 0) THEN IF (ASSOCIATED(this%level)) DEALLOCATE(this%level) #ifdef DEBUG call l4f_category_log(this%category,l4f_debug,"alloc nlevel "//to_char(nlevel)) #endif ALLOCATE(this%level(nlevel),stat=stallo) if (stallo /=0)then call l4f_category_log(this%category,l4f_fatal,"allocating memory") CALL raise_fatal_error() end if IF (linit) THEN DO i = 1, nlevel CALL init(this%level(i)) ENDDO ENDIF ENDIF ENDIF IF (PRESENT(ntimerange)) THEN IF (ntimerange >= 0) THEN IF (ASSOCIATED(this%timerange)) DEALLOCATE(this%timerange) #ifdef DEBUG call l4f_category_log(this%category,l4f_debug,"alloc ntimerange "//to_char(ntimerange)) #endif ALLOCATE(this%timerange(ntimerange),stat=stallo) if (stallo /=0)then call l4f_category_log(this%category,l4f_fatal,"allocating memory") CALL raise_fatal_error() end if IF (linit) THEN DO i = 1, ntimerange CALL init(this%timerange(i)) ENDDO ENDIF ENDIF ENDIF IF (PRESENT(nvar)) THEN IF (nvar >= 0) THEN IF (ASSOCIATED(this%var)) DEALLOCATE(this%var) #ifdef DEBUG call l4f_category_log(this%category,l4f_debug,"alloc nvar "//to_char(nvar)) #endif ALLOCATE(this%var(nvar),stat=stallo) if (stallo /=0)then call l4f_category_log(this%category,l4f_fatal,"allocating memory") CALL raise_fatal_error() end if IF (linit) THEN DO i = 1, nvar CALL init(this%var(i)) ENDDO ENDIF ENDIF ENDIF end SUBROUTINE volgrid6d_alloc !> Allocate the data array of the \a volgrid6d object. !! This method allocates the main 6-dimensional data array !! \a this%voldati and the 4-dimensional \a grid_id array \a this%gaid !! with a shape dictated by the previous call(s) to \a vol7d_alloc(). !! if any descriptor (except horizontal grid) has not been allocated !! yet, it is allocated here with a size of 1. This method should be !! explicitly used only in rare cases, it is usually called implicitly !! through the \a import interface. SUBROUTINE volgrid6d_alloc_vol(this, ini, inivol, decode) TYPE(volgrid6d),INTENT(inout) :: this !< object whose decriptors should be allocated LOGICAL,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 LOGICAL,INTENT(in),OPTIONAL :: inivol !< if provided and \c .FALSE., the allocated volumes will not be initialized to missing values LOGICAL,INTENT(in),OPTIONAL :: decode !< if provided and \c .TRUE., the \a this%voldati volume is allocated, otherwise only \a this%gaid will be allocated INTEGER :: stallo LOGICAL :: linivol #ifdef DEBUG call l4f_category_log(this%category,l4f_debug,"start alloc_vol") #endif IF (PRESENT(inivol)) THEN ! opposite condition, cannot use optio_log linivol = inivol ELSE linivol = .true. ENDIF IF (this%griddim%dim%nx > 0 .AND. this%griddim%dim%ny > 0) THEN ! allocate dimension descriptors to a minimal size for having a ! non-null volume IF (.NOT. ASSOCIATED(this%var)) CALL volgrid6d_alloc(this, nvar=1, ini=ini) IF (.NOT. ASSOCIATED(this%time)) CALL volgrid6d_alloc(this, ntime=1, ini=ini) IF (.NOT. ASSOCIATED(this%level)) CALL volgrid6d_alloc(this, nlevel=1, ini=ini) IF (.NOT. ASSOCIATED(this%timerange)) CALL volgrid6d_alloc(this, ntimerange=1, ini=ini) IF (optio_log(decode)) THEN IF (.NOT.ASSOCIATED(this%voldati)) THEN #ifdef DEBUG CALL l4f_category_log(this%category,l4f_debug,"alloc voldati volume") #endif ALLOCATE(this%voldati(this%griddim%dim%nx,this%griddim%dim%ny,& SIZE(this%level), SIZE(this%time), & SIZE(this%timerange), SIZE(this%var)),stat=stallo) IF (stallo /= 0)THEN CALL l4f_category_log(this%category,l4f_fatal,"allocating memory") CALL raise_fatal_error() ENDIF ! this is not really needed if we can check other flags for a full ! field missing values IF (linivol) this%voldati = rmiss this%voldati = rmiss ENDIF ENDIF IF (.NOT.ASSOCIATED(this%gaid)) THEN #ifdef DEBUG CALL l4f_category_log(this%category,l4f_debug,"alloc gaid volume") #endif ALLOCATE(this%gaid(SIZE(this%level), SIZE(this%time), & SIZE(this%timerange), SIZE(this%var)),stat=stallo) IF (stallo /= 0)THEN CALL l4f_category_log(this%category,l4f_fatal,"allocating memory") CALL raise_fatal_error() ENDIF IF (linivol) THEN !!$ DO i=1 ,SIZE(this%gaid,1) !!$ DO ii=1 ,SIZE(this%gaid,2) !!$ DO iii=1 ,SIZE(this%gaid,3) !!$ DO iiii=1 ,SIZE(this%gaid,4) !!$ this%gaid(i,ii,iii,iiii) = grid_id_new() ! optimize? !!$ ENDDO !!$ ENDDO !!$ ENDDO !!$ ENDDO this%gaid = grid_id_new() ENDIF ENDIF ELSE CALL l4f_category_log(this%category,l4f_fatal,'volgrid6d_alloc_vol: & &trying to allocate a volume with an invalid or unspecified horizontal grid') CALL raise_fatal_error() ENDIF END SUBROUTINE volgrid6d_alloc_vol !> Return a 2-d pointer to a x-y slice of a volume. This method works !! both with volumes having allocated and non-allocated this%voldati !! array, and it returns a pointer to a 2-d slice either from the !! allocated this%voldati array or from the grid_id object on file or !! in memory. In the second case the pointer should be either !! ALLOCATE'd to the expected size or NULLIFY'ed, and if NULLIFY'ed, !! it is allocated within the method, thus it will have to be !! deallocated by the caller when not in use anymore. Since this !! method may be called many times by a program, it is optimized for !! speed and it does not make any check about the matching size of the !! pointer and the array or about the allocation status of \a this, so !! it should be called only when everything has been checked to be in !! good shape. SUBROUTINE volgrid_get_vol_2d(this, ilevel, itime, itimerange, ivar, voldati) TYPE(volgrid6d),INTENT(in) :: this !< object from which the slice has to be retrieved INTEGER,INTENT(in) :: ilevel !< index of vertical level of the slice INTEGER,INTENT(in) :: itime !< index of time level of the slice INTEGER,INTENT(in) :: itimerange !< index of timerange of the slice INTEGER,INTENT(in) :: ivar !< index of physical variable of the slice REAL,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 IF (ASSOCIATED(this%voldati)) THEN voldati => this%voldati(:,:,ilevel,itime,itimerange,ivar) RETURN ELSE IF (.NOT.ASSOCIATED(voldati)) THEN ALLOCATE(voldati(this%griddim%dim%nx,this%griddim%dim%ny)) ENDIF CALL grid_id_decode_data(this%gaid(ilevel,itime,itimerange,ivar), voldati) ENDIF END SUBROUTINE volgrid_get_vol_2d !> Return a 3-d pointer to a x-y-z slice of a volume. This method works !! both with volumes having allocated and non-allocated this%voldati !! array, and it returns a pointer to a 3-d slice either from the !! allocated this%voldati array or from the grid_id object on file or !! in memory. In the second case the pointer should be either !! ALLOCATE'd to the expected size or NULLIFY'ed, and if NULLIFY'ed, !! it is allocated within the method, thus it will have to be !! deallocated by the caller when not in use anymore. Since this !! method may be called many times by a program, it is optimized for !! speed and it does not make any check about the matching size of the !! pointer and the array or about the allocation status of \a this, so !! it should be called only when everything has been checked to be in !! good shape. SUBROUTINE volgrid_get_vol_3d(this, itime, itimerange, ivar, voldati, zlist) TYPE(volgrid6d),INTENT(in) :: this !< object from which the slice has to be retrieved INTEGER,INTENT(in) :: itime !< index of time level of the slice INTEGER,INTENT(in) :: itimerange !< index of timerange of the slice INTEGER,INTENT(in) :: ivar !< index of physical variable of the slice REAL,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 LOGICAL,OPTIONAL,INTENT(out),ALLOCATABLE :: zlist(:) !< list of vertical level present in the data INTEGER :: ilevel LOGICAL :: done IF (ASSOCIATED(this%voldati)) THEN voldati => this%voldati(:,:,:,itime,itimerange,ivar) IF (PRESENT(zlist)) zlist(:) = c_e(this%gaid(:,itime,itimerange,ivar)) RETURN ELSE IF (.NOT.ASSOCIATED(voldati)) THEN ALLOCATE(voldati(this%griddim%dim%nx,this%griddim%dim%ny,SIZE(this%level))) ENDIF IF (PRESENT(zlist)) THEN ALLOCATE(zlist(SIZE(this%level))) ENDIF !$OMP PARALLEL DEFAULT(SHARED) !$OMP MASTER DO ilevel = 1, SIZE(this%level) !$OMP TASK FIRSTPRIVATE(ilevel),PRIVATE(done) CALL grid_id_decode_data(this%gaid(ilevel,itime,itimerange,ivar), & voldati(:,:,ilevel),done) IF (PRESENT(zlist)) zlist(ilevel) = done !$OMP END TASK ENDDO !$OMP END MASTER !$OMP END PARALLEL ENDIF END SUBROUTINE volgrid_get_vol_3d !> Reset a 2-d x-y slice of a volume after the data have been modified. !! This method works both with volumes having allocated and !! non-allocated this%voldati array, and it updates the requested !! slice. In case \a this%voldati is already allocated, this is a !! no-operation while in the other case this method encodes the field !! provided into the grid_id object on file or in memory. Since this !! method may be called many times by a program, it is optimized for !! speed and it does not make any check about the matching size of the !! field and the array or about the allocation status of \a this, so !! it should be called only when everything has been checked to be in !! good shape. SUBROUTINE volgrid_set_vol_2d(this, ilevel, itime, itimerange, ivar, voldati) TYPE(volgrid6d),INTENT(inout) :: this !< object in which slice has to be updated INTEGER,INTENT(in) :: ilevel !< index of vertical level of the slice INTEGER,INTENT(in) :: itime !< index of time level of the slice INTEGER,INTENT(in) :: itimerange !< index of timerange of the slice INTEGER,INTENT(in) :: ivar !< index of physical variable of the slice REAL,INTENT(in) :: voldati(:,:) !< updated values of the slice IF (ASSOCIATED(this%voldati)) THEN RETURN ELSE CALL grid_id_encode_data(this%gaid(ilevel,itime,itimerange,ivar), voldati) ENDIF END SUBROUTINE volgrid_set_vol_2d !> Reset a 3-d x-y-z slice of a volume after the data have been modified. !! This method works both with volumes having allocated and !! non-allocated this%voldati array, and it updates the requested !! slice. In case \a this%voldati is already allocated, this is a !! no-operation while in the other case this method encodes the field !! provided into the grid_id object on file or in memory. Since this !! method may be called many times by a program, it is optimized for !! speed and it does not make any check about the matching size of the !! field and the array or about the allocation status of \a this, so !! it should be called only when everything has been checked to be in !! good shape. SUBROUTINE volgrid_set_vol_3d(this, itime, itimerange, ivar, voldati) TYPE(volgrid6d),INTENT(inout) :: this !< object in which slice has to be updated INTEGER,INTENT(in) :: itime !< index of time level of the slice INTEGER,INTENT(in) :: itimerange !< index of timerange of the slice INTEGER,INTENT(in) :: ivar !< index of physical variable of the slice REAL,INTENT(in) :: voldati(:,:,:) !< updated values of the slice INTEGER :: ilevel IF (ASSOCIATED(this%voldati)) THEN RETURN ELSE !$OMP PARALLEL DEFAULT(SHARED) !$OMP MASTER DO ilevel = 1, SIZE(this%level) !$OMP TASK FIRSTPRIVATE(ilevel) CALL grid_id_encode_data(this%gaid(ilevel,itime,itimerange,ivar), & voldati(:,:,ilevel)) !$OMP END TASK ENDDO !$OMP END MASTER !$OMP END PARALLEL ENDIF END SUBROUTINE volgrid_set_vol_3d !> Destructor, it releases every information and memory buffer !! associated with the object. It should be called also for objects !! crated through the \a import interface. SUBROUTINE volgrid6d_delete(this) TYPE(volgrid6d),INTENT(inout) :: this INTEGER :: i, ii, iii, iiii #ifdef DEBUG call l4f_category_log(this%category,l4f_debug,"delete") #endif if (associated(this%gaid))then DO i=1 ,SIZE(this%gaid,1) DO ii=1 ,SIZE(this%gaid,2) DO iii=1 ,SIZE(this%gaid,3) DO iiii=1 ,SIZE(this%gaid,4) CALL delete(this%gaid(i,ii,iii,iiii)) ENDDO ENDDO ENDDO ENDDO DEALLOCATE(this%gaid) end if call delete(this%griddim) ! call delete(this%time) ! call delete(this%timerange) ! call delete(this%level) ! call delete(this%var) if (associated( this%time )) deallocate(this%time) if (associated( this%timerange )) deallocate(this%timerange) if (associated( this%level )) deallocate(this%level) if (associated( this%var )) deallocate(this%var) if (associated(this%voldati))deallocate(this%voldati) !chiudo il logger call l4f_category_delete(this%category) END SUBROUTINE volgrid6d_delete !>\brief Scrittura su file di un volume Volgrid6d. !! Scrittura su file unformatted di un intero volume Volgrid6d. !! Il volume viene serializzato e scritto su file. !! Il file puo' essere aperto opzionalmente dall'utente. Si possono passare !! opzionalmente unità e nome del file altrimenti assegnati internamente a dei default; se assegnati internamente !! tali parametri saranno in output. !! Se non viene fornito il nome file viene utilizzato un file di default con nome pari al nome del programma in !! esecuzione con postfisso ".vg6d". !! Come parametro opzionale c'è la description che insieme alla data corrente viene inserita nell'header del file. subroutine volgrid6d_write_on_file (this,unit,description,filename,filename_auto) TYPE(volgrid6d),INTENT(IN) :: this !< volume volgrid6d da scrivere integer,optional,intent(inout) :: unit !< unità su cui scrivere; se passata =0 ritorna il valore rielaborato (default =rielaborato internamente con getlun ) character(len=*),intent(in),optional :: filename !< nome del file su cui scrivere character(len=*),intent(out),optional :: filename_auto !< nome del file generato se "filename" è omesso character(len=*),INTENT(IN),optional :: description !< descrizione del volume integer :: lunit character(len=254) :: ldescription,arg,lfilename integer :: ntime, ntimerange, nlevel, nvar integer :: tarray(8) logical :: opened,exist #ifdef DEBUG call l4f_category_log(this%category,l4f_debug,"write on file") #endif ntime=0 ntimerange=0 nlevel=0 nvar=0 !call idate(im,id,iy) call date_and_time(values=tarray) call getarg(0,arg) if (present(description))then ldescription=description else ldescription="Volgrid6d generated by: "//trim(arg) end if if (.not. present(unit))then lunit=getunit() else if (unit==0)then lunit=getunit() unit=lunit else lunit=unit end if end if lfilename=trim(arg)//".vg6d" if (index(arg,'/',back=.true.) > 0) lfilename=lfilename(index(arg,'/',back=.true.)+1 : ) if (present(filename))then if (filename /= "")then lfilename=filename end if end if if (present(filename_auto))filename_auto=lfilename inquire(unit=lunit,opened=opened) if (.not. opened) then inquire(file=lfilename,exist=exist) if (exist) CALL raise_error('file exist; cannot open new file') if (.not.exist) open (unit=lunit,file=lfilename,form="UNFORMATTED") !print *, "opened: ",lfilename end if if (associated(this%time)) ntime=size(this%time) if (associated(this%timerange)) ntimerange=size(this%timerange) if (associated(this%level)) nlevel=size(this%level) if (associated(this%var)) nvar=size(this%var) write(unit=lunit)ldescription write(unit=lunit)tarray call write_unit( this%griddim,lunit) write(unit=lunit) ntime, ntimerange, nlevel, nvar !! prime 4 dimensioni if (associated(this%time)) call write_unit(this%time, lunit) if (associated(this%level)) write(unit=lunit)this%level if (associated(this%timerange)) write(unit=lunit)this%timerange if (associated(this%var)) write(unit=lunit)this%var !! Volumi di valori dati if (associated(this%voldati)) write(unit=lunit)this%voldati if (.not. present(unit)) close(unit=lunit) end subroutine volgrid6d_write_on_file !>\brief Lettura da file di un volume Volgrid6d. !! Lettura da file unformatted di un intero volume Volgrid6d. !! Questa subroutine comprende volgrid6d_alloc e volgrid6d_alloc_vol. !! Il file puo' essere aperto opzionalmente dall'utente. Si possono passare !! opzionalmente unità e nome del file altrimenti assegnati internamente a dei default; se assegnati internamente !! tali parametri saranno in output. subroutine volgrid6d_read_from_file (this,unit,filename,description,tarray,filename_auto) TYPE(volgrid6d),INTENT(OUT) :: this !< Volume volgrid6d da leggere integer,intent(inout),optional :: unit !< unità su cui è stato aperto un file; se =0 rielaborato internamente (default = elaborato internamente con getunit) character(len=*),INTENT(in),optional :: filename !< nome del file eventualmente da aprire (default = nome dell'eseguibile.v7d) character(len=*),intent(out),optional :: filename_auto !< nome del file generato se "filename" è omesso character(len=*),INTENT(out),optional :: description !< descrizione del volume letto integer,intent(out),optional :: tarray(8) !< vettore come definito da "date_and_time" della data di scrittura del volume integer :: ntime, ntimerange, nlevel, nvar character(len=254) :: ldescription,lfilename,arg integer :: ltarray(8),lunit logical :: opened,exist #ifdef DEBUG call l4f_category_log(this%category,l4f_debug,"read from file") #endif call getarg(0,arg) if (.not. present(unit))then lunit=getunit() else if (unit==0)then lunit=getunit() unit=lunit else lunit=unit end if end if lfilename=trim(arg)//".vg6d" if (index(arg,'/',back=.true.) > 0) lfilename=lfilename(index(arg,'/',back=.true.)+1 : ) if (present(filename))then if (filename /= "")then lfilename=filename end if end if if (present(filename_auto))filename_auto=lfilename inquire(unit=lunit,opened=opened) if (.not. opened) then inquire(file=lfilename,exist=exist) IF (.NOT. exist) CALL raise_fatal_error('file '//trim(lfilename)//' does not exist, cannot open') open (unit=lunit,file=lfilename,form="UNFORMATTED") end if read(unit=lunit)ldescription read(unit=lunit)ltarray call l4f_log(l4f_info,"Info: reading volgrid6d from file: "//trim(lfilename)) call l4f_log(l4f_info,"Info: description: "//trim(ldescription)) !call l4f_log("Info: written on ",ltarray) if (present(description))description=ldescription if (present(tarray))tarray=ltarray call read_unit( this%griddim,lunit) read(unit=lunit) ntime, ntimerange, nlevel, nvar call volgrid6d_alloc (this, & ntime=ntime, ntimerange=ntimerange, nlevel=nlevel, nvar=nvar) call volgrid6d_alloc_vol (this) if (associated(this%time)) call read_unit(this%time, lunit) if (associated(this%level)) read(unit=lunit)this%level if (associated(this%timerange)) read(unit=lunit)this%timerange if (associated(this%var)) read(unit=lunit)this%var !! Volumi di valori if (associated(this%voldati)) read(unit=lunit)this%voldati if (.not. present(unit)) close(unit=lunit) end subroutine volgrid6d_read_from_file !> Import a single \a gridinfo object into a \a volgrid6d object. !! This methods imports a single gridded field from a \a gridinfo !! object into a \a volgrid6d object, inserting it into the !! multidimensional structure of volgrid6d. The volgrid6d object must !! have been already initialized and the dimensions specified with !! volgrid6d_alloc(). If the \a force argument is missing or \a !! .FALSE. , the volgrid6d object dimension descriptors (time, !! timerange, vertical level, physical variable) must already have !! space for the corresponding values coming from gridinfo, otherwise !! the object will be rejected; this means that all the volgrid6d !! dimension descriptors should be correctly assigned. If \a force is !! \a .TRUE. , the gridinfo dimension descriptors that do not fit into !! available descriptors in the \a volgrid6d structure, will be !! accomodated in a empty (i.e. equal to missing value) descriptor, if !! available, otherwise the gridinfo will be rejected. The !! descriptor of the grid in the \a volgrid object is assigned to the !! descriptor contained in \a gridinfo if it is missing in \a volgrid, !! otherwise it is checked and the object is rejected if grids do not !! match. SUBROUTINE import_from_gridinfo(this, gridinfo, force, dup_mode, clone, & isanavar) TYPE(volgrid6d),INTENT(inout) :: this !< object in which to import TYPE(gridinfo_def),INTENT(in) :: gridinfo !< gridinfo object to be imported LOGICAL,INTENT(in),OPTIONAL :: force !< if provided and \c .TRUE., the gridinfo is forced into an empty element of \a this, if required and possible INTEGER,INTENT(in),OPTIONAL :: dup_mode !< determines the behavior in case of duplicate metadata: if \a dup_mode is not provided or 0, a duplicate field overwrites, if \a dup_mode is 1, duplicate fields are merged with priority to the last LOGICAL , INTENT(in),OPTIONAL :: clone !< if provided and \c .TRUE. , clone the gaid's from \a gridinfo to \a this LOGICAL,INTENT(IN),OPTIONAL :: isanavar !< if provides and \a .TRUE., the gridinfo object is treated as time-independent and replicated for every time and timerange CHARACTER(len=255) :: type INTEGER :: itime0, itimerange0, itime1, itimerange1, itime, itimerange, & ilevel, ivar, ldup_mode LOGICAL :: dup TYPE(datetime) :: correctedtime TYPE(vol7d_timerange) :: correctedtimerange REAL,ALLOCATABLE :: tmpgrid(:,:) IF (PRESENT(dup_mode)) THEN ldup_mode = dup_mode ELSE ldup_mode = 0 ENDIF call get_val(this%griddim,proj_type=type) #ifdef DEBUG call l4f_category_log(this%category,l4f_debug,"import_from_gridinfo: "//trim(type)) #endif if (.not. c_e(type))then call copy(gridinfo%griddim, this%griddim) ! ho gia` fatto init, altrimenti non potrei fare la get_val(this%griddim) ! per cui meglio non ripetere ! call init(this,gridinfo%griddim,categoryappend) CALL volgrid6d_alloc_vol(this, ini=.true.) ! decode? else if (.not. (this%griddim == gridinfo%griddim ))then CALL l4f_category_log(this%category,l4f_error, & "volgrid and gridinfo grid type or size are different, gridinfo rejected") CALL raise_error() RETURN end if ! Cerco gli indici del campo da inserire, se non trovo metto nel primo missing ilevel = index(this%level, gridinfo%level) IF (ilevel == 0 .AND. optio_log(force)) THEN ilevel = index(this%level, vol7d_level_miss) IF (ilevel /= 0) this%level(ilevel) = gridinfo%level ENDIF IF (ilevel == 0) THEN CALL l4f_category_log(this%category,l4f_error, & "volgrid6d: level not valid for volume, gridinfo rejected") CALL raise_error() RETURN ENDIF IF (optio_log(isanavar)) THEN ! assign to all times and timeranges itime0 = 1 itime1 = SIZE(this%time) itimerange0 = 1 itimerange1 = SIZE(this%timerange) ELSE ! usual case correctedtime = gridinfo%time IF (this%time_definition == 1 .OR. this%time_definition == 2) correctedtime = correctedtime + & timedelta_new(sec=gridinfo%timerange%p1) itime0 = index(this%time, correctedtime) IF (itime0 == 0 .AND. optio_log(force)) THEN itime0 = index(this%time, datetime_miss) IF (itime0 /= 0) this%time(itime0) = correctedtime ENDIF IF (itime0 == 0) THEN CALL l4f_category_log(this%category,l4f_error, & "volgrid6d: time not valid for volume, gridinfo rejected") CALL raise_error() RETURN ENDIF itime1 = itime0 correctedtimerange = gridinfo%timerange IF (this%time_definition == 2) correctedtimerange%p1 = 0 itimerange0 = index(this%timerange, correctedtimerange) IF (itimerange0 == 0 .AND. optio_log(force)) THEN itimerange0 = index(this%timerange, vol7d_timerange_miss) IF (itimerange0 /= 0) this%timerange(itimerange0) = gridinfo%timerange ENDIF IF (itimerange0 == 0) THEN CALL l4f_category_log(this%category,l4f_error, & "volgrid6d: timerange not valid for volume, gridinfo rejected") CALL raise_error() RETURN ENDIF itimerange1 = itimerange0 ENDIF ivar = index(this%var, gridinfo%var) IF (ivar == 0 .AND. optio_log(force)) THEN ivar = index(this%var, volgrid6d_var_miss) IF (ivar /= 0) this%var(ivar) = gridinfo%var ENDIF IF (ivar == 0) THEN CALL l4f_category_log(this%category,l4f_error, & "volgrid6d: var not valid for volume, gridinfo rejected") CALL raise_error() RETURN ENDIF DO itimerange = itimerange0, itimerange1 DO itime = itime0, itime1 IF (ASSOCIATED(this%gaid)) THEN dup = .false. IF (c_e(this%gaid(ilevel,itime,itimerange,ivar))) THEN dup = .true. CALL l4f_category_log(this%category,l4f_warn,"gaid exist: grib duplicated") ! avoid memory leaks IF (optio_log(clone)) CALL delete(this%gaid(ilevel,itime,itimerange,ivar)) ENDIF IF (optio_log(clone)) THEN CALL copy(gridinfo%gaid, this%gaid(ilevel,itime,itimerange,ivar)) #ifdef DEBUG CALL l4f_category_log(this%category,l4f_debug,"cloning to a new gaid") #endif ELSE this%gaid(ilevel,itime,itimerange,ivar) = gridinfo%gaid ENDIF IF (ASSOCIATED(this%voldati))THEN IF (.NOT.dup .OR. ldup_mode == 0) THEN this%voldati(:,:,ilevel,itime,itimerange,ivar) = decode_gridinfo(gridinfo) ELSE IF (ldup_mode == 1) THEN tmpgrid = decode_gridinfo(gridinfo) ! f2003 automatic allocation WHERE(c_e(tmpgrid)) this%voldati(:,:,ilevel,itime,itimerange,ivar) = tmpgrid(:,:) END WHERE ELSE IF (ldup_mode == 2) THEN WHERE(.NOT.c_e(this%voldati(:,:,ilevel,itime,itimerange,ivar))) this%voldati(:,:,ilevel,itime,itimerange,ivar) = decode_gridinfo(gridinfo) END WHERE ENDIF ENDIF ELSE CALL l4f_category_log(this%category,l4f_error, & "gaid not allocated, you probably need to call volgrid6d_alloc_vol first") CALL raise_error() RETURN ENDIF ENDDO ENDDO END SUBROUTINE import_from_gridinfo !> Export a single grid of a \a volgrid6d object to a \a gridinfo_def object. !! A single 2d slice of a \a volgrid6d at a specified location is !! written into a \a gridinfo_def object, including the \a grid_id !! which can be used for the successive export to file. SUBROUTINE export_to_gridinfo(this, gridinfo, itime, itimerange, ilevel, ivar, & gaid_template, clone) TYPE(volgrid6d),INTENT(in) :: this !< volume to be exported TYPE(gridinfo_def),INTENT(inout) :: gridinfo !< output gridinfo_def object INTEGER :: itime !< index within \a this of the element to be exported for the time dimension INTEGER :: itimerange !< index within \a this of the element to be exported for the timerange dimension INTEGER :: ilevel !< index within \a this of the element to be exported for the vertical level dimension INTEGER :: ivar !< index within \a this of the element to be exported for the variable dimension TYPE(grid_id),INTENT(in),OPTIONAL :: gaid_template !< \a grid_id template to be used for output data replacing the one contained in \a this LOGICAL,INTENT(in),OPTIONAL :: clone !< if provided and \c .TRUE., clone the grid_id included in \a this rather than making a shallow copy TYPE(grid_id) :: gaid LOGICAL :: usetemplate REAL,POINTER :: voldati(:,:) TYPE(datetime) :: correctedtime #ifdef DEBUG CALL l4f_category_log(this%category,l4f_debug,"export_to_gridinfo") #endif IF (.NOT.c_e(this%gaid(ilevel,itime,itimerange,ivar))) THEN #ifdef DEBUG CALL l4f_category_log(this%category,l4f_debug,"empty gaid found, skipping export") #endif RETURN ENDIF usetemplate = .false. IF (PRESENT(gaid_template)) THEN CALL copy(gaid_template, gaid) #ifdef DEBUG CALL l4f_category_log(this%category,l4f_debug,"template cloned to a new gaid") #endif usetemplate = c_e(gaid) ENDIF IF (.NOT.usetemplate) THEN IF (optio_log(clone)) THEN CALL copy(this%gaid(ilevel,itime,itimerange,ivar), gaid) #ifdef DEBUG CALL l4f_category_log(this%category,l4f_debug,"original gaid cloned to a new one") #endif ELSE gaid = this%gaid(ilevel,itime,itimerange,ivar) ENDIF ENDIF IF (this%time_definition == 1 .OR. this%time_definition == 2) THEN correctedtime = this%time(itime) - & timedelta_new(sec=this%timerange(itimerange)%p1) ELSE correctedtime = this%time(itime) ENDIF CALL init(gridinfo,gaid, this%griddim, correctedtime, this%timerange(itimerange), & this%level(ilevel), this%var(ivar)) ! reset the gridinfo, bad but necessary at this point for encoding the field CALL export(gridinfo%griddim, gridinfo%gaid) ! encode the field IF (ASSOCIATED(this%voldati)) THEN CALL encode_gridinfo(gridinfo, this%voldati(:,:,ilevel,itime,itimerange,ivar)) ELSE IF (usetemplate) THEN ! field must be forced into template in this case NULLIFY(voldati) CALL volgrid_get_vol_2d(this, ilevel, itime, itimerange, ivar, voldati) CALL encode_gridinfo(gridinfo, voldati) DEALLOCATE(voldati) ENDIF END SUBROUTINE export_to_gridinfo !> Import an array of \a gridinfo objects into an array of \a !! volgrid6d objects. !! This method imports an array of gridded fields from an \a !! arrayof_gridinfo object into a suitable number of \a volgrid6d !! objects. The number of \a volgrid6d allocated is determined by the !! number of different grids encountered in \a arrayof_gridinfo. !! Unlike the import for a single \a gridinfo, here the \a volgrid6d !! object is a non-associated pointer to a 1-d array of uninitialized !! objects, and all the dimension descriptors in each of the objects !! are allocated and assigned within the method according to the data !! contained in \a gridinfov. !! If the \a anavar array argument is provided, all the input messages !! whose variable maps to one of the B-table variables contained in \a !! anavar are treated as time-independent (AKA anagraphic data, !! station data, etc.), thus their time and timerange are ignored and !! they are replicated for every time and timerange present in the !! corresponding data volume. SUBROUTINE import_from_gridinfovv(this, gridinfov, dup_mode, clone, decode, & time_definition, anavar, categoryappend) TYPE(volgrid6d),POINTER :: this(:) !< object in which to import TYPE(arrayof_gridinfo),INTENT(in) :: gridinfov !< array of gridinfo objects to be imported INTEGER,INTENT(in),OPTIONAL :: dup_mode !< determines the behavior in case of duplicate metadata: if \a dup_mode is not provided or 0, a duplicate field overwrites, if \a dup_mode is 1, duplicate fields are merged with priority to the last LOGICAL , INTENT(in),OPTIONAL :: clone !< if provided and \c .TRUE. , clone the gaid's from \a gridinfo to \a this LOGICAL,INTENT(in),OPTIONAL :: decode !< if provided and \a .FALSE. the data volume in the elements of \a this is not allocated and successive work will be performed on grid_id's INTEGER,INTENT(IN),OPTIONAL :: time_definition !< 0=time is reference time; 1=time is validity time CHARACTER(len=*),INTENT(IN),OPTIONAL :: anavar(:) !< list of variables (B-table code equivalent) to be treated as time-independent data CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category INTEGER :: i, j, stallo INTEGER :: ngrid, ntime, ntimerange, nlevel, nvar, ltime_definition TYPE(l4f_handle) :: category CHARACTER(len=512) :: a_name TYPE(datetime),ALLOCATABLE :: correctedtime(:) LOGICAL,ALLOCATABLE :: isanavar(:) TYPE(vol7d_var) :: lvar TYPE(vol7d_timerange),ALLOCATABLE :: correctedtimerange(:) ! category temporanea (altrimenti non possiamo loggare) if (present(categoryappend))then call l4f_launcher(a_name,a_name_append=trim(subcategory)//"."//trim(categoryappend)) else call l4f_launcher(a_name,a_name_append=trim(subcategory)) endif category=l4f_category_get_handle(a_name) #ifdef DEBUG call l4f_category_log(category,l4f_debug,"start import_from_gridinfovv") #endif IF (PRESENT(time_definition)) THEN ltime_definition = max(min(time_definition, 2), 0) ELSE ltime_definition = 0 ENDIF ngrid=count_distinct(gridinfov%array(1:gridinfov%arraysize)%griddim,back=.true.) CALL l4f_category_log(category,l4f_info, t2c(ngrid)// & ' different grid definition(s) found in input data') ALLOCATE(this(ngrid),stat=stallo) IF (stallo /= 0)THEN CALL l4f_category_log(category,l4f_fatal,"allocating memory") CALL raise_fatal_error() ENDIF DO i = 1, ngrid IF (PRESENT(categoryappend))THEN CALL init(this(i), time_definition=ltime_definition, categoryappend=trim(categoryappend)//"-vol"//t2c(i)) ELSE CALL init(this(i), time_definition=ltime_definition, categoryappend="vol"//t2c(i)) ENDIF ENDDO this(:)%griddim=pack_distinct(gridinfov%array(1:gridinfov%arraysize)%griddim, & ngrid, back=.true.) ! mark elements as ana variables (time-independent) ALLOCATE(isanavar(gridinfov%arraysize)) isanavar(:) = .false. IF (PRESENT(anavar)) THEN DO i = 1, gridinfov%arraysize DO j = 1, SIZE(anavar) lvar = convert(gridinfov%array(i)%var) IF (lvar%btable == anavar(j)) THEN isanavar(i) = .true. EXIT ENDIF ENDDO ENDDO CALL l4f_category_log(category,l4f_info,t2c(count(isanavar))//'/'// & t2c(gridinfov%arraysize)//' constant-data messages found in input data') ENDIF IF (ltime_definition == 1 .OR. ltime_definition == 2) THEN ! verification time ALLOCATE(correctedtime(gridinfov%arraysize)) correctedtime(:) = gridinfov%array(1:gridinfov%arraysize)%time DO i = 1, gridinfov%arraysize correctedtime(i) = correctedtime(i) + & timedelta_new(sec=gridinfov%array(i)%timerange%p1) ENDDO ENDIF IF (ltime_definition == 2) THEN ! set all to analysis ALLOCATE(correctedtimerange(gridinfov%arraysize)) correctedtimerange(:) = gridinfov%array(1:gridinfov%arraysize)%timerange correctedtimerange(:)%p1 = 0 ENDIF DO i = 1, ngrid IF (PRESENT(anavar)) THEN j = count((this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) & .AND. .NOT.isanavar(:)) IF (j <= 0) THEN CALL l4f_category_log(category, l4f_fatal, 'grid n.'//t2c(i)// & ' has only constant data, this is not allowed') CALL l4f_category_log(category, l4f_fatal, 'please check anavar argument') CALL raise_fatal_error() ENDIF ENDIF IF (ltime_definition == 1 .OR. ltime_definition == 2) THEN ! verification time ntime = count_distinct(correctedtime, & mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) & .AND. .NOT.isanavar(:), back=.true.) ELSE ntime = count_distinct(gridinfov%array(1:gridinfov%arraysize)%time, & mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) & .AND. .NOT.isanavar(:), back=.true.) ENDIF IF (ltime_definition == 2) THEN ! set all to analysis ntimerange = count_distinct(correctedtimerange, & mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) & .AND. .NOT.isanavar(:), back=.true.) ELSE ntimerange = count_distinct(gridinfov%array(1:gridinfov%arraysize)%timerange, & mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) & .AND. .NOT.isanavar(:), back=.true.) ENDIF nlevel = count_distinct(gridinfov%array(1:gridinfov%arraysize)%level, & mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim), & back=.true.) nvar = count_distinct(gridinfov%array(1:gridinfov%arraysize)%var, & mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim), & back=.true.) #ifdef DEBUG CALL l4f_category_log(this(i)%category,l4f_debug,"alloc volgrid6d index: "//t2c(i)) #endif CALL volgrid6d_alloc(this(i),this(i)%griddim%dim,ntime=ntime, & ntimerange=ntimerange,nlevel=nlevel,nvar=nvar) IF (ltime_definition == 1 .OR. ltime_definition == 2) THEN ! verification time this(i)%time = pack_distinct(correctedtime, ntime, & mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) & .AND. .NOT.isanavar(:), back=.true.) ELSE this(i)%time = pack_distinct(gridinfov%array(1:gridinfov%arraysize)%time, ntime, & mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) & .AND. .NOT.isanavar(:), back=.true.) ENDIF CALL sort(this(i)%time) IF (ltime_definition == 2) THEN ! set all to analysis this(i)%timerange = pack_distinct(correctedtimerange, ntimerange, & mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) & .AND. .NOT.isanavar(:), back=.true.) ELSE this(i)%timerange = pack_distinct(gridinfov%array(1:gridinfov%arraysize)%timerange, & ntimerange, mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) & .AND. .NOT.isanavar(:), back=.true.) ENDIF CALL sort(this(i)%timerange) this(i)%level=pack_distinct(gridinfov%array(1:gridinfov%arraysize)%level, & nlevel,mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim), & back=.true.) CALL sort(this(i)%level) this(i)%var=pack_distinct(gridinfov%array(1:gridinfov%arraysize)%var, nvar, & mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim), & back=.true.) #ifdef DEBUG CALL l4f_category_log(this(i)%category,l4f_debug,"alloc_vol volgrid6d index: "//t2c(i)) #endif CALL volgrid6d_alloc_vol(this(i), decode=decode) ENDDO IF (ltime_definition == 1 .OR. ltime_definition == 2) DEALLOCATE(correctedtime) IF (ltime_definition == 2) DEALLOCATE(correctedtimerange) DO i = 1, gridinfov%arraysize #ifdef DEBUG CALL l4f_category_log(category,l4f_debug,"import from gridinfov index: "//t2c(i)) CALL l4f_category_log(category,l4f_info, & "to volgrid6d index: "//t2c(index(this%griddim, gridinfov%array(i)%griddim))) #endif CALL import(this(index(this%griddim, gridinfov%array(i)%griddim)), & gridinfov%array(i), dup_mode=dup_mode, clone=clone, isanavar=isanavar(i)) ENDDO !chiudo il logger temporaneo CALL l4f_category_delete(category) END SUBROUTINE import_from_gridinfovv !> Export a \a volgrid6d object to an \a arrayof_gridinfo object. !! The multidimensional \a volgrid6d structure is serialized into a !! one-dimensional array of gridinfo_def objects, which is allocated !! to the proper size if not already allocated, or it is extended !! keeping the old data if any. SUBROUTINE export_to_gridinfov(this, gridinfov, gaid_template, clone) TYPE(volgrid6d),INTENT(inout) :: this !< volume to be exported TYPE(arrayof_gridinfo),INTENT(inout) :: gridinfov !< output array of gridinfo_def objects TYPE(grid_id),INTENT(in),OPTIONAL :: gaid_template !< \a grid_id template to be used for output data replacing the one contained in \a this LOGICAL,INTENT(in),OPTIONAL :: clone !< if provided and \c .TRUE., clone the grid_id included in \a this rather than making a shallow copy INTEGER :: i ,itime, itimerange, ilevel, ivar INTEGER :: ntime, ntimerange, nlevel, nvar TYPE(gridinfo_def) :: gridinfol #ifdef DEBUG CALL l4f_category_log(this%category,l4f_debug,"start export_to_gridinfov") #endif ! this is necessary in order not to repeatedly and uselessly copy the ! same array of coordinates for each 2d grid during export, the ! side-effect is that the computed projection in this is lost CALL dealloc(this%griddim%dim) i=0 ntime=size(this%time) ntimerange=size(this%timerange) nlevel=size(this%level) nvar=size(this%var) DO itime=1,ntime DO itimerange=1,ntimerange DO ilevel=1,nlevel DO ivar=1,nvar CALL init(gridinfol) CALL export(this, gridinfol, itime, itimerange, ilevel, ivar, & gaid_template=gaid_template, clone=clone) IF (c_e(gridinfol%gaid)) THEN CALL insert(gridinfov, gridinfol) ELSE CALL delete(gridinfol) ENDIF ENDDO ENDDO ENDDO ENDDO END SUBROUTINE export_to_gridinfov !> Export an array of \a volgrid6d objects to an \a arrayof_gridinfo object. !! The multidimensional \a volgrid6d structures are serialized into a !! one-dimensional array of gridinfo_def objects, which is allocated !! to the proper size if not already allocated, or it is extended !! keeping the old data if any. SUBROUTINE export_to_gridinfovv(this, gridinfov, gaid_template, clone) !, & ! categoryappend) TYPE(volgrid6d),INTENT(inout) :: this(:) !< volume array to be exported TYPE(arrayof_gridinfo),INTENT(inout) :: gridinfov !< output array of gridinfo_def objects TYPE(grid_id),INTENT(in),OPTIONAL :: gaid_template !< \a grid_id template to be used for output data replacing the one contained in \a this LOGICAL,INTENT(in),OPTIONAL :: clone !< if provided and \c .TRUE., clone the grid_id included in \a this rather than making a shallow copy INTEGER :: i DO i = 1, SIZE(this) #ifdef DEBUG CALL l4f_category_log(this(i)%category,l4f_debug, & "export_to_gridinfovv grid index: "//t2c(i)) #endif CALL export(this(i), gridinfov, gaid_template=gaid_template, clone=clone) ENDDO END SUBROUTINE export_to_gridinfovv !> Import the content of a supported file (like grib or gdal-supported !! format) into an array of \a volgrid6d objects. This method imports !! a set of gridded fields from a file object into a suitable number !! of \a volgrid6d objects. The data are imported by creating a !! temporary \a gridinfo object, importing it into the \a volgrid6d !! object cloning the gaid's and then destroying the gridinfo, so it !! works similarly to volgrid6d_class::import_from_gridinfovv() method. !! For a detailed explanation of the \a anavar argument, see the !! documentation of volgrid6d_class::import_from_gridinfovv() method. SUBROUTINE volgrid6d_import_from_file(this, filename, dup_mode, decode, & time_definition, anavar, categoryappend) TYPE(volgrid6d),POINTER :: this(:) !< object in which to import CHARACTER(len=*),INTENT(in) :: filename !< name of file from which to import INTEGER,INTENT(in),OPTIONAL :: dup_mode !< determines the behavior in case of duplicate metadata: if \a dup_mode is not provided or 0, a duplicate field overwrites, if \a dup_mode is 1, duplicate fields are merged with priority to the last LOGICAL,INTENT(in),OPTIONAL :: decode !< if provided and \a .FALSE. the data volume in the elements of \a this is not allocated and successive work will be performed on grid_id's INTEGER,INTENT(IN),OPTIONAL :: time_definition !< 0=time is reference time; 1=time is validity time CHARACTER(len=*),INTENT(IN),OPTIONAL :: anavar(:) !< list of variables (B-table code equivalent) to be treated as time-independent data character(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category TYPE(arrayof_gridinfo) :: gridinfo TYPE(l4f_handle) :: category CHARACTER(len=512) :: a_name NULLIFY(this) IF (PRESENT(categoryappend))THEN CALL l4f_launcher(a_name,a_name_append= & trim(subcategory)//"."//trim(categoryappend)) ELSE CALL l4f_launcher(a_name,a_name_append=trim(subcategory)) ENDIF category=l4f_category_get_handle(a_name) CALL import(gridinfo, filename=filename, categoryappend=categoryappend) IF (gridinfo%arraysize > 0) THEN CALL import(this, gridinfo, dup_mode=dup_mode, clone=.true., decode=decode, & time_definition=time_definition, anavar=anavar, & categoryappend=categoryappend) CALL l4f_category_log(category,l4f_info,"deleting gridinfo") CALL delete(gridinfo) ELSE CALL l4f_category_log(category,l4f_info,"file does not contain gridded data") ENDIF ! close logger CALL l4f_category_delete(category) END SUBROUTINE volgrid6d_import_from_file !> High level method for exporting a volume array to file. !! All the information contained into an array of \a volgrid6d !! objects, i.e. dimension descriptors and data, is exported to a file !! using the proper output driver (typically grib_api for grib !! format). If a template is provided, it will determine the !! characteristic of the output file, otherwise the \a grid_id !! descriptors contained in the volgrid6d object will be used SUBROUTINE volgrid6d_export_to_file(this, filename, gaid_template, categoryappend) TYPE(volgrid6d) :: this(:) !< volume(s) to be exported CHARACTER(len=*),INTENT(in) :: filename !< output file name TYPE(grid_id),INTENT(in),OPTIONAL :: gaid_template !< template for the output file, if provided the grid_id information stored in the volgrid6d objects will be ignored character(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category TYPE(arrayof_gridinfo) :: gridinfo TYPE(l4f_handle) :: category CHARACTER(len=512) :: a_name IF (PRESENT(categoryappend)) THEN CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//"."//trim(categoryappend)) ELSE CALL l4f_launcher(a_name,a_name_append=trim(subcategory)) ENDIF category=l4f_category_get_handle(a_name) #ifdef DEBUG CALL l4f_category_log(category,l4f_debug,"start export to file") #endif CALL l4f_category_log(category,l4f_info,"writing volgrid6d to grib file: "//trim(filename)) !IF (ASSOCIATED(this)) THEN CALL export(this, gridinfo, gaid_template=gaid_template, clone=.true.) IF (gridinfo%arraysize > 0) THEN CALL export(gridinfo, filename) CALL delete(gridinfo) ENDIF !ELSE ! CALL l4f_category_log(category,L4F_INFO,"volume volgrid6d is not associated") !ENDIF ! close logger CALL l4f_category_delete(category) END SUBROUTINE volgrid6d_export_to_file !> Array destructor for \a volgrid6d class. !! Delete an array of \a volgrid6d objects and deallocate the array !! itself. SUBROUTINE volgrid6dv_delete(this) TYPE(volgrid6d),POINTER :: this(:) !< vector of volgrid6d object INTEGER :: i IF (ASSOCIATED(this)) THEN DO i = 1, SIZE(this) #ifdef DEBUG CALL l4f_category_log(this(i)%category,l4f_debug, & "delete volgrid6d vector index: "//trim(to_char(i))) #endif CALL delete(this(i)) ENDDO DEALLOCATE(this) ENDIF END SUBROUTINE volgrid6dv_delete ! Internal method for performing grid to grid computations SUBROUTINE volgrid6d_transform_compute(this, volgrid6d_in, volgrid6d_out, & lev_out, var_coord_vol, clone) TYPE(grid_transform),INTENT(in) :: this ! oggetto di trasformazione per il grigliato type(volgrid6d), INTENT(in) :: volgrid6d_in ! oggetto da trasformare type(volgrid6d), INTENT(inout) :: volgrid6d_out ! oggetto trasformato; deve essere completo (init, alloc, alloc_vol) TYPE(vol7d_level),INTENT(in),OPTIONAL :: lev_out(:) ! vol7d_level object defining target vertical grid, for vertical interpolations INTEGER,INTENT(in),OPTIONAL :: var_coord_vol ! index of variable defining vertical coordinate values in input volume LOGICAL,INTENT(in),OPTIONAL :: clone ! se fornito e \c .TRUE., clona i gaid da volgrid6d_in a volgrid6d_out INTEGER :: ntime, ntimerange, inlevel, onlevel, nvar, & itime, itimerange, ilevel, ivar, levshift, levused, lvar_coord_vol, spos REAL,POINTER :: voldatiin(:,:,:), voldatiout(:,:,:), coord_3d_in(:,:,:) TYPE(vol7d_level) :: output_levtype LOGICAL,ALLOCATABLE :: zlist(:) #ifdef DEBUG call l4f_category_log(volgrid6d_in%category,l4f_debug,"start volgrid6d_transform_compute") #endif ntime=0 ntimerange=0 inlevel=0 onlevel=0 nvar=0 lvar_coord_vol = optio_i(var_coord_vol) if (associated(volgrid6d_in%time))then ntime=size(volgrid6d_in%time) volgrid6d_out%time=volgrid6d_in%time end if if (associated(volgrid6d_in%timerange))then ntimerange=size(volgrid6d_in%timerange) volgrid6d_out%timerange=volgrid6d_in%timerange end if IF (ASSOCIATED(volgrid6d_in%level))THEN inlevel=SIZE(volgrid6d_in%level) ENDIF IF (PRESENT(lev_out)) THEN onlevel=SIZE(lev_out) volgrid6d_out%level=lev_out ELSE IF (ASSOCIATED(volgrid6d_in%level))THEN onlevel=SIZE(volgrid6d_in%level) volgrid6d_out%level=volgrid6d_in%level ENDIF if (associated(volgrid6d_in%var))then nvar=size(volgrid6d_in%var) volgrid6d_out%var=volgrid6d_in%var end if ! allocate once for speed IF (.NOT.ASSOCIATED(volgrid6d_in%voldati)) THEN ALLOCATE(voldatiin(volgrid6d_in%griddim%dim%nx, volgrid6d_in%griddim%dim%ny, & inlevel)) ENDIF IF (.NOT.ASSOCIATED(volgrid6d_out%voldati)) THEN ALLOCATE(voldatiout(volgrid6d_out%griddim%dim%nx, volgrid6d_out%griddim%dim%ny, & onlevel)) ENDIF CALL get_val(this, levshift=levshift, levused=levused) spos = imiss IF (c_e(lvar_coord_vol)) THEN CALL get_val(this%trans, output_levtype=output_levtype) IF (output_levtype%level1 == 103 .OR. output_levtype%level1 == 108) THEN ! here, unlike below, level (1,101) is not accepted because this is ! likely for pressure difference with surface 108, not 103 spos = firsttrue(volgrid6d_in%level(:) == vol7d_level_new(1)) IF (spos == 0) THEN CALL l4f_category_log(volgrid6d_in%category, l4f_error, & 'output level '//t2c(output_levtype%level1)// & ' requested, but height/press of surface not provided in volume') ENDIF IF (.NOT.c_e(levshift) .AND. .NOT.c_e(levused)) THEN CALL l4f_category_log(volgrid6d_in%category, l4f_error, & 'internal inconsistence, levshift and levused undefined when they should be') ENDIF ENDIF ENDIF DO ivar=1,nvar ! IF (c_e(var_coord_vol)) THEN ! IF (ivar == var_coord_vol) CYCLE ! skip coordinate variable in output ! ENDIF DO itimerange=1,ntimerange DO itime=1,ntime ! skip empty columns where possible, improve IF (c_e(levshift) .AND. c_e(levused)) THEN IF (.NOT.any(c_e( & volgrid6d_in%gaid(levshift+1:levshift+levused,itime,itimerange,ivar) & ))) cycle ENDIF DO ilevel = 1, min(inlevel,onlevel) ! if present gaid copy it IF (c_e(volgrid6d_in%gaid(ilevel,itime,itimerange,ivar)) .AND. .NOT. & c_e(volgrid6d_out%gaid(ilevel,itime,itimerange,ivar))) THEN IF (optio_log(clone)) THEN CALL copy(volgrid6d_in%gaid(ilevel,itime,itimerange,ivar),& volgrid6d_out%gaid(ilevel,itime,itimerange,ivar)) #ifdef DEBUG CALL l4f_category_log(volgrid6d_in%category,l4f_debug, & "cloning gaid, level "//t2c(ilevel)) #endif ELSE volgrid6d_out%gaid(ilevel,itime,itimerange,ivar) = & volgrid6d_in%gaid(ilevel,itime,itimerange,ivar) ENDIF ENDIF ENDDO ! if out level are more, we have to clone anyway DO ilevel = min(inlevel,onlevel) + 1, onlevel IF (c_e(volgrid6d_in%gaid(inlevel,itime,itimerange,ivar)) .AND. .NOT. & c_e(volgrid6d_out%gaid(ilevel,itime,itimerange,ivar))) then CALL copy(volgrid6d_in%gaid(inlevel,itime,itimerange,ivar),& volgrid6d_out%gaid(ilevel,itime,itimerange,ivar)) #ifdef DEBUG CALL l4f_category_log(volgrid6d_in%category,l4f_debug, & "forced cloning gaid, level "//t2c(inlevel)//"->"//t2c(ilevel)) #endif ENDIF ENDDO IF (c_e(lvar_coord_vol)) THEN NULLIFY(coord_3d_in) CALL volgrid_get_vol_3d(volgrid6d_in, itime, itimerange, lvar_coord_vol, & coord_3d_in) IF (c_e(spos)) THEN ! compute difference wrt surface coordinate IF (spos == 0) THEN ! error condition, set all to missing and goodnight coord_3d_in(:,:,levshift+1:levshift+levused) = rmiss ELSE DO ilevel = levshift+1, levshift+levused WHERE(c_e(coord_3d_in(:,:,ilevel)) .AND. c_e(coord_3d_in(:,:,spos))) coord_3d_in(:,:,ilevel) = coord_3d_in(:,:,ilevel) - & coord_3d_in(:,:,spos) ELSEWHERE coord_3d_in(:,:,ilevel) = rmiss END WHERE ENDDO ENDIF ENDIF ENDIF CALL volgrid_get_vol_3d(volgrid6d_in, itime, itimerange, ivar, & voldatiin, zlist) IF (ASSOCIATED(volgrid6d_out%voldati)) & ! improve!!!! CALL volgrid_get_vol_3d(volgrid6d_out, itime, itimerange, ivar, & voldatiout) IF (c_e(lvar_coord_vol)) THEN ! this is a vertint case, zlist is useless CALL compute(this, voldatiin, voldatiout, var=convert(volgrid6d_in%var(ivar)), & coord_3d_in=coord_3d_in(:,:,levshift+1:levshift+levused)) ! subset coord_3d_in ELSE CALL compute(this, voldatiin, voldatiout, var=convert(volgrid6d_in%var(ivar)), zlist=zlist) ENDIF CALL volgrid_set_vol_3d(volgrid6d_out, itime, itimerange, ivar, & voldatiout) ENDDO ENDDO ENDDO IF (c_e(lvar_coord_vol)) THEN DEALLOCATE(coord_3d_in) ENDIF IF (.NOT.ASSOCIATED(volgrid6d_in%voldati)) THEN DEALLOCATE(voldatiin) ENDIF IF (.NOT.ASSOCIATED(volgrid6d_out%voldati)) THEN DEALLOCATE(voldatiout) ENDIF END SUBROUTINE volgrid6d_transform_compute !> Performs the specified abstract transformation on the data provided. !! The abstract transformation is specified by \a this parameter; the !! corresponding specifical transformation (\a grid_transform object) !! is created and destroyed internally. The output transformed object !! is created internally and it does not require preliminary !! initialisation. SUBROUTINE volgrid6d_transform(this, griddim, volgrid6d_in, volgrid6d_out, & lev_out, volgrid6d_coord_in, maskgrid, maskbounds, clone, decode, categoryappend) TYPE(transform_def),INTENT(in) :: this !< object specifying the abstract transformation TYPE(griddim_def),INTENT(in),OPTIONAL :: griddim !< griddim specifying the output grid (required by most transformation types) TYPE(volgrid6d),INTENT(inout) :: volgrid6d_in !< object to be transformed, it is not modified, despite the INTENT(inout) TYPE(volgrid6d),INTENT(out) :: volgrid6d_out !< transformed object, it does not require initialisation TYPE(vol7d_level),INTENT(in),OPTIONAL,TARGET :: lev_out(:) !< vol7d_level object defining target vertical grid, for vertical interpolations TYPE(volgrid6d),INTENT(in),OPTIONAL :: volgrid6d_coord_in !< object providing time constant input vertical coordinate for some kind of vertical interpolations REAL,INTENT(in),OPTIONAL :: maskgrid(:,:) !< 2D field to be used for defining subareas according to its values, it must have the same shape as the field to be interpolated (for transformation subtype 'maskfill') REAL,INTENT(in),OPTIONAL :: maskbounds(:) !< array of boundary values for defining a subset of valid points where the values of \a maskgrid are within the first and last value of \a maskbounds (for transformation type 'metamorphosis:maskfill') LOGICAL,INTENT(in),OPTIONAL :: clone !< if provided and \a .TRUE. , clone the \a gaid's from \a volgrid6d_in to \a volgrid6d_out LOGICAL,INTENT(in),OPTIONAL :: decode !< determine whether the data in \a volgrid6d_out should be decoded or remain coded in gaid, if not provided, the decode status is taken from \a volgrid6d_in CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category TYPE(grid_transform) :: grid_trans TYPE(vol7d_level),POINTER :: llev_out(:) TYPE(vol7d_level) :: input_levtype, output_levtype TYPE(vol7d_var) :: vcoord_var INTEGER :: i, k, ntime, ntimerange, nlevel, nvar, var_coord_in, var_coord_vol, & cf_out, nxc, nyc, nxi, nyi, i3, i4, i5, i6, & ulstart, ulend, spos REAL,ALLOCATABLE :: coord_3d_in(:,:,:) TYPE(geo_proj) :: proj_in, proj_out CHARACTER(len=80) :: trans_type LOGICAL :: ldecode LOGICAL,ALLOCATABLE :: mask_in(:) #ifdef DEBUG call l4f_category_log(volgrid6d_in%category, l4f_debug, "start volgrid6d_transform") #endif ntime=0 ntimerange=0 nlevel=0 nvar=0 if (associated(volgrid6d_in%time)) ntime=size(volgrid6d_in%time) if (associated(volgrid6d_in%timerange)) ntimerange=size(volgrid6d_in%timerange) if (associated(volgrid6d_in%level)) nlevel=size(volgrid6d_in%level) if (associated(volgrid6d_in%var)) nvar=size(volgrid6d_in%var) IF (ntime == 0 .OR. ntimerange == 0 .OR. nlevel == 0 .OR. nvar == 0) THEN CALL l4f_category_log(volgrid6d_in%category, l4f_error, & "trying to transform an incomplete volgrid6d object, ntime="//t2c(ntime)// & ' ntimerange='//t2c(ntimerange)//' nlevel='//t2c(nlevel)//' nvar='//t2c(nvar)) CALL init(volgrid6d_out) ! initialize to empty CALL raise_error() RETURN ENDIF CALL get_val(this, trans_type=trans_type) ! store desired output component flag and unrotate if necessary cf_out = imiss IF (PRESENT(griddim) .AND. (trans_type == 'inter' .OR. trans_type == 'boxinter' & .OR. trans_type == 'stencilinter')) THEN ! improve condition!! CALL get_val(volgrid6d_in%griddim, proj=proj_in) CALL get_val(griddim, component_flag=cf_out, proj=proj_out) ! if different projections wind components must be referred to geographical system IF (proj_in /= proj_out) CALL vg6d_wind_unrot(volgrid6d_in) ELSE IF (PRESENT(griddim)) THEN ! just get component_flag, the rest is rubbish CALL get_val(griddim, component_flag=cf_out) ENDIF var_coord_in = imiss var_coord_vol = imiss IF (trans_type == 'vertint') THEN IF (PRESENT(lev_out)) THEN ! if volgrid6d_coord_in provided and allocated, check that it fits IF (PRESENT(volgrid6d_coord_in)) THEN IF (ASSOCIATED(volgrid6d_coord_in%voldati)) THEN ! strictly 1 time and 1 timerange IF (SIZE(volgrid6d_coord_in%voldati,4) /= 1 .OR. & SIZE(volgrid6d_coord_in%voldati,5) /= 1) THEN CALL l4f_category_log(volgrid6d_in%category, l4f_error, & 'volume providing constant input vertical coordinate must have & &only 1 time and 1 timerange') CALL init(volgrid6d_out) ! initialize to empty CALL raise_error() RETURN ENDIF ! search for variable providing vertical coordinate CALL get_val(this, output_levtype=output_levtype) vcoord_var = vol7d_var_new(vol7d_level_to_var(output_levtype)) IF (.NOT.c_e(vcoord_var)) THEN CALL l4f_category_log(volgrid6d_in%category, l4f_error, & 'requested output level type '//t2c(output_levtype%level1)// & ' does not correspond to any known physical variable for & &providing vertical coordinate') CALL init(volgrid6d_out) ! initialize to empty CALL raise_error() RETURN ENDIF DO i = 1, SIZE(volgrid6d_coord_in%var) IF (convert(volgrid6d_coord_in%var(i)) == vcoord_var) THEN var_coord_in = i EXIT ENDIF ENDDO IF (.NOT.c_e(var_coord_in)) THEN CALL l4f_category_log(volgrid6d_in%category, l4f_error, & 'volume providing constant input vertical coordinate contains no & &variables matching output level type '//t2c(output_levtype%level1)) CALL init(volgrid6d_out) ! initialize to empty CALL raise_error() RETURN ENDIF CALL l4f_category_log(volgrid6d_in%category, l4f_info, & 'Coordinate for vertint found in coord volume at position '// & t2c(var_coord_in)) ! check horizontal grid ! this is too strict (component flag and so on) ! IF (volgrid6d_coord_in%griddim /= volgrid6d_in%griddim) THEN CALL get_val(volgrid6d_coord_in%griddim, nx=nxc, ny=nyc) CALL get_val(volgrid6d_in%griddim, nx=nxi, ny=nyi) IF (nxc /= nxi .OR. nyc /= nyi) THEN CALL l4f_category_log(volgrid6d_in%category, l4f_error, & 'volume providing constant input vertical coordinate must have & &the same grid as the input') CALL l4f_category_log(volgrid6d_in%category, l4f_error, & 'vertical coordinate: '//t2c(nxc)//'x'//t2c(nyc)// & ', input volume: '//t2c(nxi)//'x'//t2c(nyi)) CALL init(volgrid6d_out) ! initialize to empty CALL raise_error() RETURN ENDIF ! check vertical coordinate system CALL get_val(this, input_levtype=input_levtype) mask_in = & ! implicit allocation (volgrid6d_coord_in%level(:)%level1 == input_levtype%level1) .AND. & (volgrid6d_coord_in%level(:)%level2 == input_levtype%level2) ulstart = firsttrue(mask_in) ulend = lasttrue(mask_in) IF (ulstart == 0 .OR. ulend == 0) THEN CALL l4f_category_log(volgrid6d_in%category, l4f_error, & 'coordinate file does not contain levels of type '// & t2c(input_levtype%level1)//'/'//t2c(input_levtype%level2)// & ' specified for input data') CALL init(volgrid6d_out) ! initialize to empty CALL raise_error() RETURN ENDIF coord_3d_in = volgrid6d_coord_in%voldati(:,:,ulstart:ulend,1,1,var_coord_in) ! implicit allocation ! special case IF (output_levtype%level1 == 103 .OR. & output_levtype%level1 == 108) THEN ! surface coordinate needed ! allow coding of orography level as in Icon spos = firsttrue(volgrid6d_coord_in%level(:) == vol7d_level_new(level1=1) .OR. & volgrid6d_coord_in%level(:) == vol7d_level_new(level1=1, level2=101)) IF (spos == 0) THEN CALL l4f_category_log(volgrid6d_in%category, l4f_error, & 'output level '//t2c(output_levtype%level1)// & ' requested, but height/press of surface not provided in coordinate file') CALL init(volgrid6d_out) ! initialize to empty CALL raise_error() RETURN ENDIF DO k = 1, SIZE(coord_3d_in,3) WHERE(c_e(coord_3d_in(:,:,k)) .AND. & c_e(volgrid6d_coord_in%voldati(:,:,spos,1,1,var_coord_in))) coord_3d_in(:,:,k) = coord_3d_in(:,:,k) - & volgrid6d_coord_in%voldati(:,:,spos,1,1,var_coord_in) ELSEWHERE coord_3d_in(:,:,k) = rmiss END WHERE ENDDO ENDIF ENDIF ENDIF IF (.NOT.c_e(var_coord_in)) THEN ! search for coordinate within volume ! search for variable providing vertical coordinate CALL get_val(this, output_levtype=output_levtype) vcoord_var = vol7d_var_new(vol7d_level_to_var(output_levtype)) IF (c_e(vcoord_var)) THEN DO i = 1, SIZE(volgrid6d_in%var) IF (convert(volgrid6d_in%var(i)) == vcoord_var) THEN var_coord_vol = i EXIT ENDIF ENDDO IF (c_e(var_coord_vol)) THEN CALL l4f_category_log(volgrid6d_in%category, l4f_info, & 'Coordinate for vertint found in input volume at position '// & t2c(var_coord_vol)) ENDIF ENDIF ENDIF CALL init(volgrid6d_out, griddim=volgrid6d_in%griddim, & time_definition=volgrid6d_in%time_definition, categoryappend=categoryappend) IF (c_e(var_coord_in)) THEN CALL init(grid_trans, this, lev_in=volgrid6d_in%level, lev_out=lev_out, & coord_3d_in=coord_3d_in, categoryappend=categoryappend) ELSE CALL init(grid_trans, this, lev_in=volgrid6d_in%level, lev_out=lev_out, & categoryappend=categoryappend) ENDIF CALL get_val(grid_trans, output_level_auto=llev_out) ! get levels if auto-generated IF (.NOT.ASSOCIATED(llev_out)) llev_out => lev_out nlevel = SIZE(llev_out) ELSE CALL l4f_category_log(volgrid6d_in%category, l4f_error, & 'volgrid6d_transform: vertint requested but lev_out not provided') CALL init(volgrid6d_out) ! initialize to empty CALL raise_error() RETURN ENDIF ELSE CALL init(volgrid6d_out, griddim=griddim, & time_definition=volgrid6d_in%time_definition, categoryappend=categoryappend) CALL init(grid_trans, this, in=volgrid6d_in%griddim, out=volgrid6d_out%griddim, & maskgrid=maskgrid, maskbounds=maskbounds, categoryappend=categoryappend) ENDIF IF (c_e(grid_trans)) THEN ! transformation is valid CALL volgrid6d_alloc(volgrid6d_out, ntime=ntime, nlevel=nlevel, & ntimerange=ntimerange, nvar=nvar) IF (PRESENT(decode)) THEN ! explicitly set decode status ldecode = decode ELSE ! take it from input ldecode = ASSOCIATED(volgrid6d_in%voldati) ENDIF ! force decode if gaid is readonly decode_loop: DO i6 = 1,nvar DO i5 = 1, ntimerange DO i4 = 1, ntime DO i3 = 1, nlevel IF (c_e(volgrid6d_in%gaid(i3,i4,i5,i6))) THEN ldecode = ldecode .OR. grid_id_readonly(volgrid6d_in%gaid(i3,i4,i5,i6)) EXIT decode_loop ENDIF ENDDO ENDDO ENDDO ENDDO decode_loop IF (PRESENT(decode)) THEN IF (ldecode.NEQV.decode) THEN CALL l4f_category_log(volgrid6d_in%category, l4f_warn, & 'volgrid6d_transform: decode status forced to .TRUE. because driver does not allow copy') ENDIF ENDIF CALL volgrid6d_alloc_vol(volgrid6d_out, decode=ldecode) !ensure unproj was called !call griddim_unproj(volgrid6d_out%griddim) IF (trans_type == 'vertint') THEN #ifdef DEBUG CALL l4f_category_log(volgrid6d_in%category, l4f_debug, & "volgrid6d_transform: vertint to "//t2c(nlevel)//" levels") #endif CALL compute(grid_trans, volgrid6d_in, volgrid6d_out, lev_out=llev_out, & var_coord_vol=var_coord_vol, clone=clone) ELSE CALL compute(grid_trans, volgrid6d_in, volgrid6d_out, clone=clone) ENDIF IF (cf_out == 0) THEN ! unrotated components are desired CALL wind_unrot(volgrid6d_out) ! unrotate if necessary ELSE IF (cf_out == 1) THEN ! rotated components are desired CALL wind_rot(volgrid6d_out) ! rotate if necessary ENDIF ELSE ! should log with grid_trans%category, but it is private CALL l4f_category_log(volgrid6d_in%category, l4f_error, & 'volgrid6d_transform: transformation not valid') CALL raise_error() ENDIF CALL delete (grid_trans) END SUBROUTINE volgrid6d_transform !> Performs the specified abstract transformation on the arrays of !! data provided. The abstract transformation is specified by \a this !! parameter; the corresponding specifical transformation (\a !! grid_transform object) is created and destroyed internally. The !! output transformed object is created internally and it does not !! require preliminary initialisation. According to the input data and !! to the transformation type, the output array may have of one or !! more \a volgrid6d elements on different grids. SUBROUTINE volgrid6dv_transform(this, griddim, volgrid6d_in, volgrid6d_out, & lev_out, volgrid6d_coord_in, maskgrid, maskbounds, clone, decode, categoryappend) TYPE(transform_def),INTENT(in) :: this !< object specifying the abstract transformation TYPE(griddim_def),INTENT(in),OPTIONAL :: griddim !< griddim specifying the output grid (required by most transformation types) TYPE(volgrid6d),INTENT(inout) :: volgrid6d_in(:) !< object to be transformed, it is an array of volgrid6d objects, each of which will be transformed, it is not modified, despite the INTENT(inout) TYPE(volgrid6d),POINTER :: volgrid6d_out(:) !< transformed object, it is a non associated pointer to an array of volgrid6d objects which will be allocated by the method TYPE(vol7d_level),INTENT(in),OPTIONAL :: lev_out(:) !< vol7d_level object defining target vertical grid TYPE(volgrid6d),INTENT(in),OPTIONAL :: volgrid6d_coord_in !< object providing time constant input vertical coordinate for some kind of vertical interpolations REAL,INTENT(in),OPTIONAL :: maskgrid(:,:) !< 2D field to be used for defining subareas according to its values, it must have the same shape as the field to be interpolated (for transformation subtype 'maskfill') REAL,INTENT(in),OPTIONAL :: maskbounds(:) !< array of boundary values for defining a subset of valid points where the values of \a maskgrid are within the first and last value of \a maskbounds (for transformation type 'metamorphosis:maskfill') LOGICAL,INTENT(in),OPTIONAL :: clone !< if provided and \a .TRUE. , clone the \a gaid's from \a volgrid6d_in to \a volgrid6d_out LOGICAL,INTENT(in),OPTIONAL :: decode !< if provided and \a .FALSE. the data volume is not allocated, but work is performed on grid_id's CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category INTEGER :: i, stallo allocate(volgrid6d_out(size(volgrid6d_in)),stat=stallo) if (stallo /= 0)then call l4f_log(l4f_fatal,"allocating memory") call raise_fatal_error() end if do i=1,size(volgrid6d_in) call transform(this, griddim, volgrid6d_in(i), volgrid6d_out(i), & lev_out=lev_out, volgrid6d_coord_in=volgrid6d_coord_in, & maskgrid=maskgrid, maskbounds=maskbounds, & clone=clone, decode=decode, categoryappend=categoryappend) end do END SUBROUTINE volgrid6dv_transform ! Internal method for performing grid to sparse point computations SUBROUTINE volgrid6d_v7d_transform_compute(this, volgrid6d_in, vol7d_out, & networkname, noconvert) TYPE(grid_transform),INTENT(in) :: this ! oggetto di trasformazione per grigliato type(volgrid6d), INTENT(in) :: volgrid6d_in ! oggetto da trasformare type(vol7d), INTENT(inout) :: vol7d_out ! oggetto trasformato CHARACTER(len=*),OPTIONAL,INTENT(in) :: networkname ! imposta il network in vol7d_out (default='generic') LOGICAL,OPTIONAL,INTENT(in) :: noconvert !< do not try to match variable and convert values during transform INTEGER :: nntime, nana, ntime, ntimerange, nlevel, nvar, stallo INTEGER :: itime, itimerange, ivar, inetwork REAL,ALLOCATABLE :: voldatir_out(:,:,:) TYPE(conv_func),POINTER :: c_func(:) TYPE(datetime),ALLOCATABLE :: validitytime(:,:) REAL,POINTER :: voldatiin(:,:,:) #ifdef DEBUG call l4f_category_log(volgrid6d_in%category,l4f_debug,"start volgrid6d_v7d_transform_compute") #endif ntime=0 ntimerange=0 nlevel=0 nvar=0 NULLIFY(c_func) if (present(networkname))then call init(vol7d_out%network(1),name=networkname) else call init(vol7d_out%network(1),name='generic') end if if (associated(volgrid6d_in%timerange))then ntimerange=size(volgrid6d_in%timerange) vol7d_out%timerange=volgrid6d_in%timerange end if if (associated(volgrid6d_in%time))then ntime=size(volgrid6d_in%time) if (vol7d_out%time_definition == volgrid6d_in%time_definition) then ! i time sono definiti uguali: assegno vol7d_out%time=volgrid6d_in%time else ! converto reference in validity allocate (validitytime(ntime,ntimerange),stat=stallo) if (stallo /=0)then call l4f_category_log(volgrid6d_in%category,l4f_fatal,"allocating memory") call raise_fatal_error() end if do itime=1,ntime do itimerange=1,ntimerange if (vol7d_out%time_definition > volgrid6d_in%time_definition) then validitytime(itime,itimerange) = & volgrid6d_in%time(itime) + timedelta_new(sec=volgrid6d_in%timerange(itimerange)%p1) else validitytime(itime,itimerange) = & volgrid6d_in%time(itime) - timedelta_new(sec=volgrid6d_in%timerange(itimerange)%p1) end if end do end do nntime = count_distinct(reshape(validitytime,(/ntime*ntimerange/)), back=.true.) vol7d_out%time=pack_distinct(reshape(validitytime,(/ntime*ntimerange/)), nntime,back=.true.) end if end if IF (ASSOCIATED(volgrid6d_in%level)) THEN nlevel = SIZE(volgrid6d_in%level) vol7d_out%level=volgrid6d_in%level ENDIF IF (ASSOCIATED(volgrid6d_in%var)) THEN nvar = SIZE(volgrid6d_in%var) IF (.NOT. optio_log(noconvert)) THEN CALL vargrib2varbufr(volgrid6d_in%var, vol7d_out%dativar%r, c_func) ENDIF ENDIF nana = SIZE(vol7d_out%ana) ! allocate once for speed IF (.NOT.ASSOCIATED(volgrid6d_in%voldati)) THEN ALLOCATE(voldatiin(volgrid6d_in%griddim%dim%nx, volgrid6d_in%griddim%dim%ny, & nlevel)) ENDIF ALLOCATE(voldatir_out(nana,1,nlevel),stat=stallo) IF (stallo /= 0) THEN CALL l4f_category_log(volgrid6d_in%category,l4f_fatal,"allocating memory") CALL raise_fatal_error() ENDIF inetwork=1 do itime=1,ntime do itimerange=1,ntimerange ! do ilevel=1,nlevel do ivar=1,nvar !non è chiaro se questa sezione è utile o no !ossia il compute sotto sembra prevedere voldatir_out solo in out !!$ if (vol7d_out%time_definition == volgrid6d_in%time_definition) then !!$ voldatir_out=reshape(vol7d_out%voldatir(:,itime,ilevel,itimerange,ivar,inetwork),(/nana,1/)) !!$ else !!$ voldatir_out=reshape(vol7d_out%voldatir(:,index(vol7d_out%time,validitytime(itime,ilevel)),ilevel,itimerange,ivar,inetwork),(/nana,1/)) !!$ end if CALL volgrid_get_vol_3d(volgrid6d_in, itime, itimerange, ivar, & voldatiin) CALL compute(this, voldatiin, voldatir_out, vol7d_out%dativar%r(ivar)) if (vol7d_out%time_definition == volgrid6d_in%time_definition) then vol7d_out%voldatir(:,itime,:,itimerange,ivar,inetwork) = & voldatir_out(:,1,:) else vol7d_out%voldatir(:,index(vol7d_out%time,validitytime(itime,itimerange)),:,itimerange,ivar,inetwork)=& reshape(voldatir_out,(/nana,nlevel/)) end if ! 1 indice della dimensione "anagrafica" ! 2 indice della dimensione "tempo" ! 3 indice della dimensione "livello verticale" ! 4 indice della dimensione "intervallo temporale" ! 5 indice della dimensione "variabile" ! 6 indice della dimensione "rete" end do ! end do end do end do deallocate(voldatir_out) IF (.NOT.ASSOCIATED(volgrid6d_in%voldati)) THEN DEALLOCATE(voldatiin) ENDIF if (allocated(validitytime)) deallocate(validitytime) ! Rescale valid data according to variable conversion table IF (ASSOCIATED(c_func)) THEN DO ivar = 1, nvar CALL compute(c_func(ivar), vol7d_out%voldatir(:,:,:,:,ivar,:)) ENDDO DEALLOCATE(c_func) ENDIF end SUBROUTINE volgrid6d_v7d_transform_compute !> Performs the specified abstract transformation on the data provided. !! The abstract transformation is specified by \a this parameter; the !! corresponding specifical transformation (\a grid_transform object) !! is created and destroyed internally. The output transformed object !! is created internally and it does not require preliminary !! initialisation. SUBROUTINE volgrid6d_v7d_transform(this, volgrid6d_in, vol7d_out, v7d, & maskgrid, maskbounds, networkname, noconvert, find_index, categoryappend) TYPE(transform_def),INTENT(in) :: this !< object specifying the abstract transformation TYPE(volgrid6d),INTENT(inout) :: volgrid6d_in !< object to be transformed, it is not modified, despite the INTENT(inout) TYPE(vol7d),INTENT(out) :: vol7d_out !< transformed object, it does not requires initialisation TYPE(vol7d),INTENT(in),OPTIONAL :: v7d !< object containing a list of points over which transformation has to be done (required by some transformation types) REAL,INTENT(in),OPTIONAL :: maskgrid(:,:) !< 2D field to be used for defining subareas according to its values, it must have the same shape as the field to be interpolated (for transformation type 'maskinter') REAL,INTENT(in),OPTIONAL :: maskbounds(:) !< array of boundary values for defining subareas from the values of \a maskgrid, the number of subareas is SIZE(maskbounds) - 1, if not provided a default based on extreme values of \a makgrid is used CHARACTER(len=*),OPTIONAL,INTENT(in) :: networkname !< set the output network name in vol7d_out (default='generic') LOGICAL,OPTIONAL,INTENT(in) :: noconvert !< do not try to match variable and convert values during transform PROCEDURE(basic_find_index),POINTER,OPTIONAL :: find_index CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category type(grid_transform) :: grid_trans INTEGER :: ntime, ntimerange, nlevel, nvar, nana, time_definition, nnetwork, stallo INTEGER :: itime, itimerange, inetwork TYPE(datetime),ALLOCATABLE :: validitytime(:,:) INTEGER,ALLOCATABLE :: point_index(:) TYPE(vol7d) :: v7d_locana #ifdef DEBUG call l4f_category_log(volgrid6d_in%category,l4f_debug,"start volgrid6d_v7d_transform") #endif call vg6d_wind_unrot(volgrid6d_in) ntime=0 ntimerange=0 nlevel=0 nvar=0 nnetwork=1 call get_val(this,time_definition=time_definition) if (.not. c_e(time_definition)) then time_definition=1 ! default to validity time endif IF (PRESENT(v7d)) THEN CALL vol7d_copy(v7d, v7d_locana) ELSE CALL init(v7d_locana) ENDIF if (associated(volgrid6d_in%timerange)) ntimerange=size(volgrid6d_in%timerange) if (associated(volgrid6d_in%time)) then ntime=size(volgrid6d_in%time) if (time_definition /= volgrid6d_in%time_definition) then ! converto reference in validity allocate (validitytime(ntime,ntimerange),stat=stallo) if (stallo /=0)then call l4f_category_log(volgrid6d_in%category,l4f_fatal,"allocating memory") call raise_fatal_error() end if do itime=1,ntime do itimerange=1,ntimerange if (time_definition > volgrid6d_in%time_definition) then validitytime(itime,itimerange) = & volgrid6d_in%time(itime) + timedelta_new(sec=volgrid6d_in%timerange(itimerange)%p1) else validitytime(itime,itimerange) = & volgrid6d_in%time(itime) - timedelta_new(sec=volgrid6d_in%timerange(itimerange)%p1) end if end do end do ntime = count_distinct(reshape(validitytime,(/ntime*ntimerange/)), back=.true.) deallocate (validitytime) end if end if if (associated(volgrid6d_in%level)) nlevel=size(volgrid6d_in%level) if (associated(volgrid6d_in%var)) nvar=size(volgrid6d_in%var) CALL init(grid_trans, this, volgrid6d_in%griddim, v7d_locana, & maskgrid=maskgrid, maskbounds=maskbounds, find_index=find_index, & categoryappend=categoryappend) CALL init (vol7d_out,time_definition=time_definition) IF (c_e(grid_trans)) THEN nana=SIZE(v7d_locana%ana) CALL vol7d_alloc(vol7d_out, nana=nana, ntime=ntime, nlevel=nlevel, & ntimerange=ntimerange, ndativarr=nvar, nnetwork=nnetwork) vol7d_out%ana = v7d_locana%ana CALL get_val(grid_trans, output_point_index=point_index) IF (ALLOCATED(point_index)) THEN ! check that size(point_index) == nana? CALL vol7d_alloc(vol7d_out, nanavari=1) CALL init(vol7d_out%anavar%i(1), 'B01192') ENDIF CALL vol7d_alloc_vol(vol7d_out) IF (ALLOCATED(point_index)) THEN DO inetwork = 1, nnetwork vol7d_out%volanai(:,1,inetwork) = point_index(:) ENDDO ENDIF CALL compute(grid_trans, volgrid6d_in, vol7d_out, networkname, noconvert) ELSE CALL l4f_log(l4f_error, 'vg6d_v7d_transform: transformation not valid') CALL raise_error() ENDIF CALL delete(grid_trans) #ifdef HAVE_DBALLE ! set variables to a conformal state CALL vol7d_dballe_set_var_du(vol7d_out) #endif CALL delete(v7d_locana) END SUBROUTINE volgrid6d_v7d_transform !> Performs the specified abstract transformation on the arrays of !! data provided. The abstract transformation is specified by \a this !! parameter; the corresponding specifical transformation (\a !! grid_transform object) is created and destroyed internally. The !! output transformed object is created internally and it does not !! require preliminary initialisation. The transformation performed on !! each element of the input \a volgrid6d array object is merged into !! a single \a vol7d output object. SUBROUTINE volgrid6dv_v7d_transform(this, volgrid6d_in, vol7d_out, v7d, & maskgrid, maskbounds, networkname, noconvert, find_index, categoryappend) TYPE(transform_def),INTENT(in) :: this !< object specifying the abstract transformation TYPE(volgrid6d),INTENT(inout) :: volgrid6d_in(:) !< object to be transformed, it is an array of volgrid6d objects, each of which will be transformed, it is not modified, despite the INTENT(inout) TYPE(vol7d),INTENT(out) :: vol7d_out !< transformed object, it does not require initialisation TYPE(vol7d),INTENT(in),OPTIONAL :: v7d !< object containing a list of points over which transformation has to be done (required by some transformation types) REAL,INTENT(in),OPTIONAL :: maskgrid(:,:) !< 2D field to be used for defining subareas according to its values, it must have the same shape as the field to be interpolated (for transformation type 'maskinter') REAL,INTENT(in),OPTIONAL :: maskbounds(:) !< array of boundary values for defining subareas from the values of \a maskgrid, the number of subareas is SIZE(maskbounds) - 1, if not provided a default based on extreme values of \a makgrid is used CHARACTER(len=*),OPTIONAL,INTENT(in) :: networkname !< set the output network name in vol7d_out (default='generic') LOGICAL,OPTIONAL,INTENT(in) :: noconvert !< do not try to match variable and convert values during transform PROCEDURE(basic_find_index),POINTER,OPTIONAL :: find_index CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category integer :: i TYPE(vol7d) :: v7dtmp CALL init(v7dtmp) CALL init(vol7d_out) DO i=1,SIZE(volgrid6d_in) CALL transform(this, volgrid6d_in(i), v7dtmp, v7d=v7d, & maskgrid=maskgrid, maskbounds=maskbounds, & networkname=networkname, noconvert=noconvert, find_index=find_index, & categoryappend=categoryappend) CALL vol7d_append(vol7d_out, v7dtmp) ENDDO END SUBROUTINE volgrid6dv_v7d_transform ! Internal method for performing sparse point to grid computations SUBROUTINE v7d_volgrid6d_transform_compute(this, vol7d_in, volgrid6d_out, networkname, gaid_template) TYPE(grid_transform),INTENT(in) :: this ! object specifying the specific transformation type(vol7d), INTENT(in) :: vol7d_in ! object to be transformed type(volgrid6d), INTENT(inout) :: volgrid6d_out ! transformed object CHARACTER(len=*),OPTIONAL,INTENT(in) :: networkname ! select the network to be processed from the \a vol7d input object, default the first network TYPE(grid_id),OPTIONAL,INTENT(in) :: gaid_template ! the template (typically grib_api) to be associated with output data, it also helps in improving variable conversion integer :: nana, ntime, ntimerange, nlevel, nvar INTEGER :: ilevel, itime, itimerange, ivar, inetwork REAL,POINTER :: voldatiout(:,:,:) type(vol7d_network) :: network TYPE(conv_func), pointer :: c_func(:) !TODO category sarebbe da prendere da vol7d #ifdef DEBUG CALL l4f_category_log(volgrid6d_out%category, l4f_debug, & 'start v7d_volgrid6d_transform_compute') #endif ntime=0 ntimerange=0 nlevel=0 nvar=0 IF (PRESENT(networkname)) THEN CALL init(network,name=networkname) inetwork = index(vol7d_in%network,network) IF (inetwork <= 0) THEN CALL l4f_category_log(volgrid6d_out%category,l4f_warn, & 'network '//trim(networkname)//' not found, first network will be transformed') inetwork = 1 ENDIF ELSE inetwork = 1 ENDIF ! no time_definition conversion implemented, output will be the same as input if (associated(vol7d_in%time))then ntime=size(vol7d_in%time) volgrid6d_out%time=vol7d_in%time end if if (associated(vol7d_in%timerange))then ntimerange=size(vol7d_in%timerange) volgrid6d_out%timerange=vol7d_in%timerange end if if (associated(vol7d_in%level))then nlevel=size(vol7d_in%level) volgrid6d_out%level=vol7d_in%level end if if (associated(vol7d_in%dativar%r))then nvar=size(vol7d_in%dativar%r) CALL varbufr2vargrib(vol7d_in%dativar%r, volgrid6d_out%var, c_func, gaid_template) end if nana=SIZE(vol7d_in%voldatir, 1) ! allocate once for speed IF (.NOT.ASSOCIATED(volgrid6d_out%voldati)) THEN ALLOCATE(voldatiout(volgrid6d_out%griddim%dim%nx, volgrid6d_out%griddim%dim%ny, & nlevel)) ENDIF DO ivar=1,nvar DO itimerange=1,ntimerange DO itime=1,ntime ! clone the gaid template where I have data IF (PRESENT(gaid_template)) THEN DO ilevel = 1, nlevel IF (any(c_e(vol7d_in%voldatir(:,itime,ilevel,itimerange,ivar,inetwork)))) THEN CALL copy(gaid_template, volgrid6d_out%gaid(ilevel,itime,itimerange,ivar)) ELSE volgrid6d_out%gaid(ilevel,itime,itimerange,ivar) = grid_id_new() ENDIF ENDDO ENDIF ! get data IF (ASSOCIATED(volgrid6d_out%voldati)) & ! improve!!!! CALL volgrid_get_vol_3d(volgrid6d_out, itime, itimerange, ivar, & voldatiout) ! do the interpolation CALL compute(this, & vol7d_in%voldatir(:,itime,:,itimerange,ivar,inetwork), voldatiout, & vol7d_in%dativar%r(ivar)) ! rescale valid data according to variable conversion table IF (ASSOCIATED(c_func)) THEN CALL compute(c_func(ivar), voldatiout(:,:,:)) ENDIF ! put data CALL volgrid_set_vol_3d(volgrid6d_out, itime, itimerange, ivar, & voldatiout) ENDDO ENDDO ENDDO IF (.NOT.ASSOCIATED(volgrid6d_out%voldati)) THEN DEALLOCATE(voldatiout) ENDIF IF (ASSOCIATED(c_func)) THEN DEALLOCATE(c_func) ENDIF END SUBROUTINE v7d_volgrid6d_transform_compute !> Performs the specified abstract transformation on the data provided. !! The abstract transformation is specified by \a this parameter; the !! corresponding specifical transformation (\a grid_transform object) !! is created and destroyed internally. The output transformed object !! is created internally and it does not require preliminary !! initialisation. SUBROUTINE v7d_volgrid6d_transform(this, griddim, vol7d_in, volgrid6d_out, & networkname, gaid_template, categoryappend) TYPE(transform_def),INTENT(in) :: this !< object specifying the abstract transformation TYPE(griddim_def),INTENT(in),OPTIONAL :: griddim !< griddim specifying the output grid (required by most transformation types) ! TODO ripristinare intent(in) dopo le opportune modifiche in grid_class.F90 TYPE(vol7d),INTENT(inout) :: vol7d_in !< object to be transformed, it is not modified, despite the INTENT(inout) TYPE(volgrid6d),INTENT(out) :: volgrid6d_out !< transformed object, it does not require initialisation CHARACTER(len=*),OPTIONAL,INTENT(in) :: networkname !< select the network to be processed from the \a vol7d input object, default the first network TYPE(grid_id),OPTIONAL,INTENT(in) :: gaid_template !< the template (typically grib_api) to be associated with output data, it also helps in improving variable conversion CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category type(grid_transform) :: grid_trans integer :: ntime, ntimerange, nlevel, nvar !TODO la category sarebbe da prendere da vol7d !call l4f_category_log(vol7d_out%category,L4F_DEBUG,"start volgrid6d_transform") CALL vol7d_alloc_vol(vol7d_in) ! be safe ntime=SIZE(vol7d_in%time) ntimerange=SIZE(vol7d_in%timerange) nlevel=SIZE(vol7d_in%level) nvar=0 if (associated(vol7d_in%dativar%r)) nvar=size(vol7d_in%dativar%r) IF (nvar <= 0) THEN ! use vol7d category once it will be implemented CALL l4f_log(l4f_error, & "trying to transform a vol7d object incomplete or without real variables") CALL init(volgrid6d_out) ! initialize to empty CALL raise_error() RETURN ENDIF CALL init(grid_trans, this, vol7d_in, griddim, categoryappend=categoryappend) CALL init(volgrid6d_out, griddim, time_definition=vol7d_in%time_definition, & categoryappend=categoryappend) IF (c_e(grid_trans)) THEN CALL volgrid6d_alloc(volgrid6d_out, griddim%dim, ntime=ntime, nlevel=nlevel, & ntimerange=ntimerange, nvar=nvar) ! can I avoid decode=.TRUE. here, using gaid_template? CALL volgrid6d_alloc_vol(volgrid6d_out, decode=.true.) CALL compute(grid_trans, vol7d_in, volgrid6d_out, networkname, gaid_template) CALL vg6d_wind_rot(volgrid6d_out) ELSE ! should log with grid_trans%category, but it is private CALL l4f_log(l4f_error, 'v7d_vg6d_transform: transformation not valid') CALL raise_error() ENDIF CALL delete(grid_trans) END SUBROUTINE v7d_volgrid6d_transform ! Internal method for performing sparse point to sparse point computations SUBROUTINE v7d_v7d_transform_compute(this, vol7d_in, vol7d_out, lev_out, & var_coord_vol) TYPE(grid_transform),INTENT(in) :: this ! oggetto di trasformazione per grigliato type(vol7d), INTENT(in) :: vol7d_in ! oggetto da trasformare type(vol7d), INTENT(inout) :: vol7d_out ! oggetto trasformato TYPE(vol7d_level),INTENT(in),OPTIONAL :: lev_out(:) ! vol7d_level object defining target vertical grid, for vertical interpolations INTEGER,INTENT(in),OPTIONAL :: var_coord_vol ! index of variable defining vertical coordinate values in input volume INTEGER :: itime, itimerange, ilevel, ivar, inetwork, & levshift, levused, lvar_coord_vol, spos REAL,ALLOCATABLE :: coord_3d_in(:,:,:) TYPE(vol7d_level) :: output_levtype lvar_coord_vol = optio_i(var_coord_vol) vol7d_out%time(:) = vol7d_in%time(:) vol7d_out%timerange(:) = vol7d_in%timerange(:) IF (PRESENT(lev_out)) THEN vol7d_out%level(:) = lev_out(:) ELSE vol7d_out%level(:) = vol7d_in%level(:) ENDIF vol7d_out%network(:) = vol7d_in%network(:) IF (ASSOCIATED(vol7d_in%dativar%r)) THEN ! work only when real vars are available vol7d_out%dativar%r(:) = vol7d_in%dativar%r(:) CALL get_val(this, levshift=levshift, levused=levused) spos = imiss IF (c_e(lvar_coord_vol)) THEN CALL get_val(this%trans, output_levtype=output_levtype) IF (output_levtype%level1 == 103 .OR. output_levtype%level1 == 108) THEN spos = firsttrue(vol7d_in%level(:) == vol7d_level_new(1)) IF (spos == 0) THEN CALL l4f_log(l4f_error, & 'output level '//t2c(output_levtype%level1)// & ' requested, but height/press of surface not provided in volume') ENDIF IF (.NOT.c_e(levshift) .AND. .NOT.c_e(levused)) THEN CALL l4f_log(l4f_error, & 'internal inconsistence, levshift and levused undefined when they should be') ENDIF ALLOCATE(coord_3d_in(SIZE(vol7d_in%ana),1,SIZE(vol7d_in%level))) ENDIF ENDIF DO inetwork = 1, SIZE(vol7d_in%network) DO ivar = 1, SIZE(vol7d_in%dativar%r) DO itimerange = 1, SIZE(vol7d_in%timerange) DO itime = 1, SIZE(vol7d_in%time) ! dirty trick to make voldatir look like a 2d-array of shape (nana,1) IF (c_e(lvar_coord_vol)) THEN IF (c_e(spos)) THEN ! compute difference wrt surface coordinate IF (spos == 0) THEN ! error condition, set all to missing and goodnight coord_3d_in(:,:,levshift+1:levshift+levused) = rmiss ELSE DO ilevel = levshift+1, levshift+levused WHERE(c_e(vol7d_in%voldatir(:,itime:itime,ilevel,itimerange,lvar_coord_vol,inetwork)) .AND. & c_e(vol7d_in%voldatir(:,itime:itime,spos,itimerange,lvar_coord_vol,inetwork))) coord_3d_in(:,:,ilevel) = vol7d_in%voldatir(:,itime:itime,ilevel,itimerange,lvar_coord_vol,inetwork) - & vol7d_in%voldatir(:,itime:itime,spos,itimerange,lvar_coord_vol,inetwork) ELSEWHERE coord_3d_in(:,:,ilevel) = rmiss END WHERE ENDDO ENDIF CALL compute(this, & vol7d_in%voldatir(:,itime,:,itimerange,ivar,inetwork), & vol7d_out%voldatir(:,itime:itime,:,itimerange,ivar,inetwork), & var=vol7d_in%dativar%r(ivar), & coord_3d_in=coord_3d_in) ELSE CALL compute(this, & vol7d_in%voldatir(:,itime,:,itimerange,ivar,inetwork), & vol7d_out%voldatir(:,itime:itime,:,itimerange,ivar,inetwork), & var=vol7d_in%dativar%r(ivar), & coord_3d_in=vol7d_in%voldatir(:,itime:itime,:,itimerange, & lvar_coord_vol,inetwork)) ENDIF ELSE CALL compute(this, & vol7d_in%voldatir(:,itime,:,itimerange,ivar,inetwork), & vol7d_out%voldatir(:,itime:itime,:,itimerange,ivar,inetwork), & var=vol7d_in%dativar%r(ivar)) ENDIF ENDDO ENDDO ENDDO ENDDO ENDIF END SUBROUTINE v7d_v7d_transform_compute !> Performs the specified abstract transformation on the data provided. !! The abstract transformation is specified by \a this parameter; the !! corresponding specifical transformation (\a grid_transform object) !! is created and destroyed internally. The output transformed object !! is created internally and it does not require preliminary !! initialisation. The success of the transformation can be checked !! with the \a c_e method: c_e(vol7d_out). SUBROUTINE v7d_v7d_transform(this, vol7d_in, vol7d_out, v7d, maskbounds, & lev_out, vol7d_coord_in, categoryappend) TYPE(transform_def),INTENT(in) :: this !< object specifying the abstract transformation TYPE(vol7d),INTENT(inout) :: vol7d_in !< object to be transformed, it is not modified, despite the INTENT(inout) TYPE(vol7d),INTENT(out) :: vol7d_out !< transformed object, it does not require initialisation TYPE(vol7d),INTENT(in),OPTIONAL :: v7d !< object containing a list of points over which transformation has to be done (required by some transformation types) REAL,INTENT(in),OPTIONAL :: maskbounds(:) !< array of boundary values for defining a subset of valid points where the values of \a maskgrid are within the first and last value of \a maskbounds (for transformation type 'metamorphosis:maskfill') TYPE(vol7d_level),INTENT(in),OPTIONAL,TARGET :: lev_out(:) !< vol7d_level object defining target vertical grid, for vertical interpolations TYPE(vol7d),INTENT(in),OPTIONAL :: vol7d_coord_in !< object providing time constant input vertical coordinate for some kind of vertical interpolations CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend !< append this suffix to log4fortran namespace category INTEGER :: nvar, inetwork TYPE(grid_transform) :: grid_trans TYPE(vol7d_level),POINTER :: llev_out(:) TYPE(vol7d_level) :: input_levtype, output_levtype TYPE(vol7d_var) :: vcoord_var REAL,ALLOCATABLE :: coord_3d_in(:,:,:) INTEGER :: var_coord_in, var_coord_vol, i, k, ulstart, ulend, spos INTEGER,ALLOCATABLE :: point_index(:) TYPE(vol7d) :: v7d_locana, vol7d_tmpana CHARACTER(len=80) :: trans_type LOGICAL,ALLOCATABLE :: mask_in(:), point_mask(:) CALL vol7d_alloc_vol(vol7d_in) ! be safe nvar=0 IF (ASSOCIATED(vol7d_in%dativar%r)) nvar=SIZE(vol7d_in%dativar%r) CALL init(v7d_locana) IF (PRESENT(v7d)) v7d_locana = v7d CALL init(vol7d_out, time_definition=vol7d_in%time_definition) CALL get_val(this, trans_type=trans_type) var_coord_vol = imiss IF (trans_type == 'vertint') THEN IF (PRESENT(lev_out)) THEN ! if vol7d_coord_in provided and allocated, check that it fits var_coord_in = -1 IF (PRESENT(vol7d_coord_in)) THEN IF (ASSOCIATED(vol7d_coord_in%voldatir) .AND. & ASSOCIATED(vol7d_coord_in%dativar%r)) THEN ! strictly 1 time, 1 timerange and 1 network IF (SIZE(vol7d_coord_in%voldatir,2) /= 1 .OR. & SIZE(vol7d_coord_in%voldatir,4) /= 1 .OR. & SIZE(vol7d_coord_in%voldatir,6) /= 1) THEN CALL l4f_log(l4f_error, & 'volume providing constant input vertical coordinate must have & &only 1 time, 1 timerange and 1 network') CALL raise_error() RETURN ENDIF ! search for variable providing vertical coordinate CALL get_val(this, output_levtype=output_levtype) vcoord_var = vol7d_var_new(vol7d_level_to_var(output_levtype)) IF (.NOT.c_e(vcoord_var)) THEN CALL l4f_log(l4f_error, & 'requested output level type '//t2c(output_levtype%level1)// & ' does not correspond to any known physical variable for & &providing vertical coordinate') CALL raise_error() RETURN ENDIF var_coord_in = index(vol7d_coord_in%dativar%r, vcoord_var) IF (var_coord_in <= 0) THEN CALL l4f_log(l4f_error, & 'volume providing constant input vertical coordinate contains no & &real variables matching output level type '//t2c(output_levtype%level1)) CALL raise_error() RETURN ENDIF CALL l4f_log(l4f_info, & 'Coordinate for vertint found in coord volume at position '// & t2c(var_coord_in)) ! check vertical coordinate system CALL get_val(this, input_levtype=input_levtype) mask_in = & ! implicit allocation (vol7d_coord_in%level(:)%level1 == input_levtype%level1) .AND. & (vol7d_coord_in%level(:)%level2 == input_levtype%level2) ulstart = firsttrue(mask_in) ulend = lasttrue(mask_in) IF (ulstart == 0 .OR. ulend == 0) THEN CALL l4f_log(l4f_error, & 'coordinate file does not contain levels of type '// & t2c(input_levtype%level1)//'/'//t2c(input_levtype%level2)// & ' specified for input data') CALL raise_error() RETURN ENDIF coord_3d_in = vol7d_coord_in%voldatir(:,1:1,ulstart:ulend,1,var_coord_in,1) ! dirty 1:1, implicit allocation ! special case IF (output_levtype%level1 == 103 & .OR. output_levtype%level1 == 108) THEN ! surface coordinate needed spos = firsttrue(vol7d_coord_in%level(:) == vol7d_level_new(1)) IF (spos == 0) THEN CALL l4f_log(l4f_error, & 'output level '//t2c(output_levtype%level1)// & ' requested, but height/press of surface not provided in coordinate file') CALL raise_error() RETURN ENDIF DO k = 1, SIZE(coord_3d_in,3) WHERE(c_e(coord_3d_in(:,:,k)) .AND. & c_e(vol7d_coord_in%voldatir(:,1:1,spos,1,var_coord_in,1))) coord_3d_in(:,:,k) = coord_3d_in(:,:,k) - & vol7d_coord_in%voldatir(:,1:1,spos,1,var_coord_in,1) ELSEWHERE coord_3d_in(:,:,k) = rmiss END WHERE ENDDO ENDIF ENDIF ENDIF IF (var_coord_in <= 0) THEN ! search for coordinate within volume ! search for variable providing vertical coordinate CALL get_val(this, output_levtype=output_levtype) vcoord_var = vol7d_var_new(vol7d_level_to_var(output_levtype)) IF (c_e(vcoord_var)) THEN DO i = 1, SIZE(vol7d_in%dativar%r) IF (vol7d_in%dativar%r(i) == vcoord_var) THEN var_coord_vol = i EXIT ENDIF ENDDO IF (c_e(var_coord_vol)) THEN CALL l4f_log(l4f_info, & 'Coordinate for vertint found in input volume at position '// & t2c(var_coord_vol)) ENDIF ENDIF ENDIF IF (var_coord_in > 0) THEN CALL init(grid_trans, this, lev_in=vol7d_in%level, lev_out=lev_out, & coord_3d_in=coord_3d_in, categoryappend=categoryappend) ELSE CALL init(grid_trans, this, lev_in=vol7d_in%level, lev_out=lev_out, & categoryappend=categoryappend) ENDIF CALL get_val(grid_trans, output_level_auto=llev_out) ! get levels if auto-generated IF (.NOT.associated(llev_out)) llev_out => lev_out IF (c_e(grid_trans)) then! .AND. nvar > 0) THEN CALL vol7d_alloc(vol7d_out, nana=SIZE(vol7d_in%ana), & ntime=SIZE(vol7d_in%time), ntimerange=SIZE(vol7d_in%timerange), & nlevel=SIZE(llev_out), nnetwork=SIZE(vol7d_in%network), ndativarr=nvar) vol7d_out%ana(:) = vol7d_in%ana(:) CALL vol7d_alloc_vol(vol7d_out) ! no need to check c_e(var_coord_vol) here since the presence of ! this%coord_3d_in (external) has precedence over coord_3d_in internal ! in grid_transform_compute CALL compute(grid_trans, vol7d_in, vol7d_out, llev_out, & var_coord_vol=var_coord_vol) ELSE CALL l4f_log(l4f_error, 'v7d_v7d_transform: transformation not valid') CALL raise_error() ENDIF ELSE CALL l4f_log(l4f_error, & 'v7d_v7d_transform: vertint requested but lev_out not provided') CALL raise_error() ENDIF ELSE CALL init(grid_trans, this, vol7d_in, v7d_locana, maskbounds=maskbounds, & categoryappend=categoryappend) ! if this init is successful, I am sure that v7d_locana%ana is associated IF (c_e(grid_trans)) then! .AND. nvar > 0) THEN CALL vol7d_alloc(vol7d_out, nana=SIZE(v7d_locana%ana), & ntime=SIZE(vol7d_in%time), ntimerange=SIZE(vol7d_in%timerange), & nlevel=SIZE(vol7d_in%level), nnetwork=SIZE(vol7d_in%network), ndativarr=nvar) vol7d_out%ana = v7d_locana%ana CALL get_val(grid_trans, point_mask=point_mask, output_point_index=point_index) IF (ALLOCATED(point_index)) THEN CALL vol7d_alloc(vol7d_out, nanavari=1) CALL init(vol7d_out%anavar%i(1), 'B01192') ENDIF CALL vol7d_alloc_vol(vol7d_out) IF (ALLOCATED(point_index)) THEN DO inetwork = 1, SIZE(vol7d_in%network) vol7d_out%volanai(:,1,inetwork) = point_index(:) ENDDO ENDIF CALL compute(grid_trans, vol7d_in, vol7d_out) IF (ALLOCATED(point_mask)) THEN ! keep full ana IF (SIZE(point_mask) /= SIZE(vol7d_in%ana)) THEN CALL l4f_log(l4f_warn, & 'v7d_v7d_transform: inconsistency in point size: '//t2c(SIZE(point_mask)) & //':'//t2c(SIZE(vol7d_in%ana))) ELSE #ifdef DEBUG CALL l4f_log(l4f_debug, 'v7d_v7d_transform: merging ana from in to out') #endif CALL vol7d_copy(vol7d_in, vol7d_tmpana, & lana=point_mask, lnetwork=(/.true./), & ltime=(/.false./), ltimerange=(/.false./), llevel=(/.false./)) CALL vol7d_append(vol7d_out, vol7d_tmpana) ENDIF ENDIF ELSE CALL l4f_log(l4f_error, 'v7d_v7d_transform: transformation not valid') CALL raise_error() ENDIF ENDIF CALL delete (grid_trans) IF (.NOT. PRESENT(v7d)) CALL delete(v7d_locana) END SUBROUTINE v7d_v7d_transform !> Unrotate the wind components. !! It converts u and v components of vector quantities relative to the !! defined grid in the direction of increasing x and y coordinates to !! u and v components relative to easterly and notherly direction. The !! original fields are overwritten. !! \todo Check and correct wind component flag (to be moved in !! griddim_def?) subroutine vg6d_wind_unrot(this) type(volgrid6d) :: this !< object containing wind to be unrotated integer :: component_flag call get_val(this%griddim,component_flag=component_flag) if (component_flag == 1) then call l4f_category_log(this%category,l4f_info, & "unrotating vector components") call vg6d_wind__un_rot(this,.false.) call set_val(this%griddim,component_flag=0) else call l4f_category_log(this%category,l4f_info, & "no need to unrotate vector components") end if end subroutine vg6d_wind_unrot !> Rotate the wind components. !! It converts u and v components of vector quantities !! relative to easterly and notherly direction to !! defined grid in the direction of increasing x and y coordinates. !! The original fields are overwritten. subroutine vg6d_wind_rot(this) type(volgrid6d) :: this !< object containing wind to be rotated integer :: component_flag call get_val(this%griddim,component_flag=component_flag) if (component_flag == 0) then call l4f_category_log(this%category,l4f_info, & "rotating vector components") call vg6d_wind__un_rot(this,.true.) call set_val(this%griddim,component_flag=1) else call l4f_category_log(this%category,l4f_info, & "no need to rotate vector components") end if end subroutine vg6d_wind_rot ! Generic UnRotate the wind components. SUBROUTINE vg6d_wind__un_rot(this,rot) TYPE(volgrid6d) :: this ! object containing wind to be (un)rotated LOGICAL :: rot ! if .true. rotate else unrotate INTEGER :: i, j, k, l, a11, a12, a21, a22, stallo double precision,pointer :: rot_mat(:,:,:) real,allocatable :: tmp_arr(:,:) REAL,POINTER :: voldatiu(:,:), voldativ(:,:) INTEGER,POINTER :: iu(:), iv(:) IF (.NOT.ASSOCIATED(this%var)) THEN CALL l4f_category_log(this%category, l4f_error, & "trying to unrotate an incomplete volgrid6d object") CALL raise_fatal_error() ! RETURN ENDIF CALL volgrid6d_var_hor_comp_index(this%var, iu, iv) IF (.NOT.ASSOCIATED(iu)) THEN CALL l4f_category_log(this%category,l4f_error, & "unrotation impossible") CALL raise_fatal_error() ! RETURN ENDIF ! Temporary workspace ALLOCATE(tmp_arr(this%griddim%dim%nx, this%griddim%dim%ny),stat=stallo) IF (stallo /= 0) THEN CALL l4f_category_log(this%category, l4f_fatal, "allocating memory") CALL raise_fatal_error() ENDIF ! allocate once for speed IF (.NOT.ASSOCIATED(this%voldati)) THEN ALLOCATE(voldatiu(this%griddim%dim%nx, this%griddim%dim%ny), & voldativ(this%griddim%dim%nx, this%griddim%dim%ny)) ENDIF CALL griddim_unproj(this%griddim) CALL wind_unrot(this%griddim, rot_mat) a11=1 if (rot)then a12=2 a21=3 else a12=3 a21=2 end if a22=4 DO l = 1, SIZE(iu) DO k = 1, SIZE(this%timerange) DO j = 1, SIZE(this%time) DO i = 1, SIZE(this%level) ! get data CALL volgrid_get_vol_2d(this, i, j, k, iu(l), voldatiu) CALL volgrid_get_vol_2d(this, i, j, k, iv(l), voldativ) ! convert units forward ! CALL compute(conv_fwd(iu(l)), voldatiu) ! CALL compute(conv_fwd(iv(l)), voldativ) ! multiply wind components by rotation matrix WHERE(voldatiu /= rmiss .AND. voldativ /= rmiss) tmp_arr(:,:) = real(voldatiu(:,:)*rot_mat(:,:,a11) + & voldativ(:,:)*rot_mat(:,:,a12)) voldativ(:,:) = real(voldatiu(:,:)*rot_mat(:,:,a21) + & voldativ(:,:)*rot_mat(:,:,a22)) voldatiu(:,:) = tmp_arr(:,:) END WHERE ! convert units backward ! CALL uncompute(conv_fwd(iu(l)), voldatiu) ! CALL uncompute(conv_fwd(iv(l)), voldativ) ! put data CALL volgrid_set_vol_2d(this, i, j, k, iu(l), voldatiu) CALL volgrid_set_vol_2d(this, i, j, k, iv(l), voldativ) ENDDO ENDDO ENDDO ENDDO IF (.NOT.ASSOCIATED(this%voldati)) THEN DEALLOCATE(voldatiu, voldativ) ENDIF DEALLOCATE(rot_mat, tmp_arr, iu, iv) END SUBROUTINE vg6d_wind__un_rot !!$ try to understand the problem: !!$ !!$ case: !!$ !!$ 1) we have only one volume: we have to provide the direction of shift !!$ compute H and traslate on it !!$ 2) we have two volumes: !!$ 1) volume U and volume V: compute H and traslate on it !!$ 2) volume U/V and volume H : translate U/V on H !!$ 3) we have tree volumes: translate U and V on H !!$ !!$ strange cases: !!$ 1) do not have U in volume U !!$ 2) do not have V in volume V !!$ 3) we have others variables more than U and V in volumes U e V !!$ !!$ !!$ so the steps are: !!$ 1) find the volumes !!$ 2) define or compute H grid !!$ 3) trasform the volumes in H !!$ N.B. !!$ case 1) for only one vg6d (U or V) is not managed, but !!$ the not pubblic subroutines will work but you have to know what you want to do !> Convert grids type C to type A. !! Use this to interpolate data from grid type C to grid type A !! Grids type are defined by Arakawa. !! !! We need to find these types of area in a vg6d array !! T area of points with temterature etc. !! U area of points with u components of winds !! V area of points with v components of winds !! !! this method works if it finds !! two volumes: !! 1) volume U and volume V: compute H and traslate on it !! 2) volume U/V and volume H : translate U/V on H !! three volumes: translate U and V on H !! !! try to work well on more datasets at once subroutine vg6d_c2a (this) TYPE(volgrid6d),INTENT(inout) :: this(:) !< vettor of volumes volgrid6d to elaborate integer :: ngrid,igrid,jgrid,ugrid,vgrid,tgrid doubleprecision :: xmin, xmax, ymin, ymax doubleprecision :: xmin_t, xmax_t, ymin_t, ymax_t doubleprecision :: step_lon_t,step_lat_t character(len=80) :: type_t,type TYPE(griddim_def):: griddim_t ngrid=size(this) do igrid=1,ngrid call init(griddim_t) call get_val(this(igrid)%griddim,xmin=xmin_t, xmax=xmax_t, ymin=ymin_t, ymax=ymax_t,proj_type=type_t) step_lon_t=(xmax_t-xmin_t)/dble(this(igrid)%griddim%dim%nx-1) step_lat_t=(ymax_t-ymin_t)/dble(this(igrid)%griddim%dim%ny-1) do jgrid=1,ngrid ugrid=imiss vgrid=imiss tgrid=imiss #ifdef DEBUG call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: search U/V/T points:"//to_char(igrid)//to_char(jgrid)) call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: test delta: "//to_char(step_lon_t)//to_char(step_lat_t)) #endif if (this(igrid)%griddim == this(jgrid)%griddim ) cycle if (this(igrid)%griddim%dim%nx == this(jgrid)%griddim%dim%nx .and. & this(igrid)%griddim%dim%ny == this(jgrid)%griddim%dim%ny ) then call get_val(this(jgrid)%griddim,xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,proj_type=type) if (type_t /= type )cycle #ifdef DEBUG call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: test U "//& to_char(xmin)//to_char(xmax)//to_char(ymin)//to_char(ymax)) call l4f_category_log(this(igrid)%category,l4f_debug,"diff coordinate lon"//& to_char(abs(xmin - (xmin_t+step_lon_t/2.d0)))//& to_char(abs(xmax - (xmax_t+step_lon_t/2.d0)))) call l4f_category_log(this(igrid)%category,l4f_debug,"diff coordinate lat"//& to_char(abs(ymin - (ymin_t+step_lat_t/2.d0)))//& to_char(abs(ymax - (ymax_t+step_lat_t/2.d0)))) #endif 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 if ( abs(ymin - ymin_t) < 1.d-3 .and. abs(ymax - ymax_t) < 1.d-3 ) then #ifdef DEBUG call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: found U") #endif ugrid=jgrid tgrid=igrid end if end if #ifdef DEBUG call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: test V "//& to_char(xmin)//to_char(xmax)//to_char(ymin)//to_char(ymax)) #endif 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 if ( abs(xmin - xmin_t) < 1.d-3 .and. abs(xmax - xmax_t) < 1.d-3 ) then #ifdef DEBUG call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: found V") #endif vgrid=jgrid tgrid=igrid end if end if ! test if we have U and V #ifdef DEBUG call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: test U and V"//& to_char(xmin_t)//to_char(xmax_t)//to_char(ymin_t)//to_char(ymax_t)//& to_char(xmin)//to_char(xmax)//to_char(ymin)//to_char(ymax)) call l4f_category_log(this(igrid)%category,l4f_debug,"UV diff coordinate lon"//& to_char(abs(xmin_t - xmin)-step_lon_t/2.d0)//& to_char(abs(xmax_t - xmax)-step_lon_t/2.d0)) call l4f_category_log(this(igrid)%category,l4f_debug,"UV diff coordinate lat"//& to_char(abs(ymin_t - ymin) -step_lat_t/2.d0)//& to_char(abs(ymax_t - ymax)-step_lat_t/2.d0)) #endif 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 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 #ifdef DEBUG call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: found U and V case up and right") #endif vgrid=jgrid ugrid=igrid call init(griddim_t,xmin=xmin, xmax=xmax, ymin=ymin_t, ymax=ymax_t) end if end if end if ! abbiamo almeno due volumi: riportiamo U e/o V su T if (c_e(ugrid)) then #ifdef DEBUG call l4f_category_log(this(igrid)%category,l4f_debug,"C grid U points: interpolate U on T "//to_char(tgrid)//to_char(ugrid)) #endif if (c_e(tgrid))then call vg6d_c2a_grid(this(ugrid),this(tgrid)%griddim,cgrid=1) else call vg6d_c2a_grid(this(ugrid),griddim_t,cgrid=1) end if call vg6d_c2a_mat(this(ugrid),cgrid=1) end if if (c_e(vgrid)) then #ifdef DEBUG call l4f_category_log(this(igrid)%category,l4f_debug,"C grid V points: interpolate V on T "//to_char(tgrid)//to_char(vgrid)) #endif if (c_e(tgrid))then call vg6d_c2a_grid(this(vgrid),this(tgrid)%griddim,cgrid=2) else call vg6d_c2a_grid(this(vgrid),griddim_t,cgrid=2) end if call vg6d_c2a_mat(this(vgrid),cgrid=2) end if end do call delete(griddim_t) end do end subroutine vg6d_c2a ! Convert C grid to A grid subroutine vg6d_c2a_grid(this,griddim_t,cgrid) type(volgrid6d),intent(inout) :: this ! object containing fields to be translated (U or V points) type(griddim_def),intent(in),optional :: griddim_t ! object containing grid of T points integer,intent(in) :: cgrid ! in C grid (Arakawa) we have 0=T,1=U,2=V points doubleprecision :: xmin, xmax, ymin, ymax doubleprecision :: step_lon,step_lat if (present(griddim_t)) then call get_val(griddim_t,xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax) call set_val(this%griddim,xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax) ! improve grid definition precision CALL griddim_setsteps(this%griddim) else select case (cgrid) case(0) #ifdef DEBUG call l4f_category_log(this%category,l4f_debug,"C grid: T points, nothing to do") #endif return case (1) #ifdef DEBUG call l4f_category_log(this%category,l4f_debug,"C grid: U points, we need interpolation") #endif call get_val(this%griddim, xmin=xmin, xmax=xmax) step_lon=(xmax-xmin)/dble(this%griddim%dim%nx-1) xmin=xmin-step_lon/2.d0 xmax=xmax-step_lon/2.d0 call set_val(this%griddim, xmin=xmin, xmax=xmax) ! improve grid definition precision CALL griddim_setsteps(this%griddim) case (2) #ifdef DEBUG call l4f_category_log(this%category,l4f_debug,"C grid: V points, we need interpolation") #endif call get_val(this%griddim, ymin=ymin, ymax=ymax) step_lat=(ymax-ymin)/dble(this%griddim%dim%ny-1) ymin=ymin-step_lat/2.d0 ymax=ymax-step_lat/2.d0 call set_val(this%griddim, ymin=ymin, ymax=ymax) ! improve grid definition precision CALL griddim_setsteps(this%griddim) case default call l4f_category_log(this%category,l4f_fatal,"C grid type not known") call raise_fatal_error () end select end if call griddim_unproj(this%griddim) end subroutine vg6d_c2a_grid ! Convert C grid to A grid subroutine vg6d_c2a_mat(this,cgrid) type(volgrid6d),intent(inout) :: this ! object containing fields to be translated integer,intent(in) :: cgrid ! in C grid (Arakawa) we have 0=T,1=U,2=V points INTEGER :: i, j, k, iv, stallo REAL,ALLOCATABLE :: tmp_arr(:,:) REAL,POINTER :: voldatiuv(:,:) IF (cgrid == 0) RETURN ! nothing to do IF (cgrid /= 1 .AND. cgrid /= 2) THEN ! wrong cgrid CALL l4f_category_log(this%category,l4f_fatal,"C grid type "// & trim(to_char(cgrid))//" not known") CALL raise_error() RETURN ENDIF ! Temporary workspace ALLOCATE(tmp_arr(this%griddim%dim%nx, this%griddim%dim%ny),stat=stallo) if (stallo /=0)then call l4f_log(l4f_fatal,"allocating memory") call raise_fatal_error() end if ! allocate once for speed IF (.NOT.ASSOCIATED(this%voldati)) THEN ALLOCATE(voldatiuv(this%griddim%dim%nx, this%griddim%dim%ny), stat=stallo) IF (stallo /= 0) THEN CALL l4f_log(l4f_fatal,"allocating memory") CALL raise_fatal_error() ENDIF ENDIF IF (cgrid == 1) THEN ! U points to H points DO iv = 1, SIZE(this%var) DO k = 1, SIZE(this%timerange) DO j = 1, SIZE(this%time) DO i = 1, SIZE(this%level) tmp_arr(:,:) = rmiss CALL volgrid_get_vol_2d(this, i, j, k, iv, voldatiuv) ! West boundary extrapolation WHERE(voldatiuv(1,:) /= rmiss .AND. voldatiuv(2,:) /= rmiss) tmp_arr(1,:) = voldatiuv(1,:) - (voldatiuv(2,:) - voldatiuv(1,:)) / 2. END WHERE ! Rest of the field WHERE(voldatiuv(1:this%griddim%dim%nx-1,:) /= rmiss .AND. & voldatiuv(2:this%griddim%dim%nx,:) /= rmiss) tmp_arr(2:this%griddim%dim%nx,:) = & (voldatiuv(1:this%griddim%dim%nx-1,:) + & voldatiuv(2:this%griddim%dim%nx,:)) / 2. END WHERE voldatiuv(:,:) = tmp_arr CALL volgrid_set_vol_2d(this, i, j, k, iv, voldatiuv) ENDDO ENDDO ENDDO ENDDO ELSE IF (cgrid == 2) THEN ! V points to H points DO iv = 1, SIZE(this%var) DO k = 1, SIZE(this%timerange) DO j = 1, SIZE(this%time) DO i = 1, SIZE(this%level) tmp_arr(:,:) = rmiss CALL volgrid_get_vol_2d(this, i, j, k, iv, voldatiuv) ! South boundary extrapolation WHERE(voldatiuv(:,1) /= rmiss .AND. voldatiuv(:,2) /= rmiss) tmp_arr(:,1) = voldatiuv(:,1) - (voldatiuv(:,2) - voldatiuv(:,1)) / 2. END WHERE ! Rest of the field WHERE(voldatiuv(:,1:this%griddim%dim%ny-1) /= rmiss .AND. & voldatiuv(:,2:this%griddim%dim%ny) /= rmiss) tmp_arr(:,2:this%griddim%dim%ny) = & (voldatiuv(:,1:this%griddim%dim%ny-1) + & voldatiuv(:,2:this%griddim%dim%ny)) / 2. END WHERE voldatiuv(:,:) = tmp_arr CALL volgrid_set_vol_2d(this, i, j, k, iv, voldatiuv) ENDDO ENDDO ENDDO ENDDO ENDIF ! tertium non datur IF (.NOT.ASSOCIATED(this%voldati)) THEN DEALLOCATE(voldatiuv) ENDIF DEALLOCATE (tmp_arr) end subroutine vg6d_c2a_mat !> \brief Display object on screen !! !! Show brief content on screen. subroutine display_volgrid6d (this) TYPE(volgrid6d),intent(in) :: this !< object to display integer :: i #ifdef DEBUG call l4f_category_log(this%category,l4f_debug,"ora mostro gridinfo " ) #endif print*,"----------------------- volgrid6d display ---------------------" call display(this%griddim) IF (ASSOCIATED(this%time))then print*,"---- time vector ----" print*,"elements=",size(this%time) do i=1, size(this%time) call display(this%time(i)) end do end IF IF (ASSOCIATED(this%timerange))then print*,"---- timerange vector ----" print*,"elements=",size(this%timerange) do i=1, size(this%timerange) call display(this%timerange(i)) end do end IF IF (ASSOCIATED(this%level))then print*,"---- level vector ----" print*,"elements=",size(this%level) do i=1, size(this%level) call display(this%level(i)) end do end IF IF (ASSOCIATED(this%var))then print*,"---- var vector ----" print*,"elements=",size(this%var) do i=1, size(this%var) call display(this%var(i)) end do end IF IF (ASSOCIATED(this%gaid))then print*,"---- gaid vector (present mask only) ----" print*,"elements=",shape(this%gaid) print*,c_e(reshape(this%gaid,(/SIZE(this%gaid)/))) end IF print*,"--------------------------------------------------------------" end subroutine display_volgrid6d !> \brief Display vector of object on screen !! !! Show brief content on screen. subroutine display_volgrid6dv (this) TYPE(volgrid6d),intent(in) :: this(:) !< vector of object to display integer :: i print*,"----------------------- volgrid6d vector ---------------------" print*,"elements=",size(this) do i=1, size(this) #ifdef DEBUG call l4f_category_log(this(i)%category,l4f_debug,"ora mostro il vettore volgrid6d" ) #endif call display(this(i)) end do print*,"--------------------------------------------------------------" end subroutine display_volgrid6dv !> Reduce some dimensions (level and timerage) for semplification (rounding). !! Vector version of vg6dv_rounding subroutine vg6dv_rounding(vg6din,vg6dout,level,timerange,nostatproc,merge) type(volgrid6d),intent(in) :: vg6din(:) !< input volume type(volgrid6d),intent(out),pointer :: vg6dout(:) !> output volume type(vol7d_level),intent(in),optional :: level(:) !< almost equal level list type(vol7d_timerange),intent(in),optional :: timerange(:) !< almost equal timerange list logical,intent(in),optional :: merge !< if there are data on more then one almost equal levels or timeranges will be merged POINT BY POINT with priority for the fird data found ordered by icreasing var index logical,intent(in),optional :: nostatproc !< do not take in account statistical processing code in timerange and P2 integer :: i allocate(vg6dout(size(vg6din))) do i = 1, size(vg6din) call vg6d_rounding(vg6din(i),vg6dout(i),level,timerange,nostatproc,merge) end do end subroutine vg6dv_rounding !> Reduce some dimensions (level and timerage) for semplification (rounding). !! You can use this for simplify and use variables in computation like alchimia !! where fields have to be on the same coordinate !! examples: !! means in time for short periods and istantaneous values !! 2 meter and surface levels !! If there are data on more then one almost equal levels or timeranges, the first var present (at least one point) !! will be taken (order is by icreasing var index). !! You can use predefined values for classic semplification !! almost_equal_levels and almost_equal_timeranges !! The level or timerange in output will be defined by the first element of level and timerange list subroutine vg6d_rounding(vg6din,vg6dout,level,timerange,nostatproc,merge) type(volgrid6d),intent(in) :: vg6din !< input volume type(volgrid6d),intent(out) :: vg6dout !> output volume type(vol7d_level),intent(in),optional :: level(:) !< almost equal level list type(vol7d_timerange),intent(in),optional :: timerange(:) !< almost equal timerange list logical,intent(in),optional :: merge !< if there are data on more then one almost equal levels or timeranges !! will be merged POINT BY POINT with priority for the fird data found ordered by icreasing var index logical,intent(in),optional :: nostatproc !< do not take in account statistical processing code in timerange and P2 integer :: ilevel,itimerange type(vol7d_level) :: roundlevel(size(vg6din%level)) type(vol7d_timerange) :: roundtimerange(size(vg6din%timerange)) roundlevel=vg6din%level if (present(level))then do ilevel = 1, size(vg6din%level) if ((any(vg6din%level(ilevel) .almosteq. level))) then roundlevel(ilevel)=level(1) end if end do end if roundtimerange=vg6din%timerange if (present(timerange))then do itimerange = 1, size(vg6din%timerange) if ((any(vg6din%timerange(itimerange) .almosteq. timerange))) then roundtimerange(itimerange)=timerange(1) end if end do end if !set istantaneous values everywere !preserve p1 for forecast time if (optio_log(nostatproc)) then roundtimerange(:)%timerange=254 roundtimerange(:)%p2=0 end if call vg6d_reduce(vg6din,vg6dout,roundlevel,roundtimerange,merge) end subroutine vg6d_rounding !> 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 subroutine vg6d_reduce(vg6din,vg6dout,roundlevel,roundtimerange,merge) type(volgrid6d),intent(in) :: vg6din !< input volume type(volgrid6d),intent(out) :: vg6dout !< output volume type(vol7d_level),intent(in) :: roundlevel(:) !< new level list to use for rounding type(vol7d_timerange),intent(in) :: roundtimerange(:) !< new timerange list to use for rounding logical,intent(in),optional :: merge !< if there are data on equal levels or timeranges will be merged POINT BY POINT with priority for the first data found ordered by icreasing var index (require to decode all the data) integer :: nlevel,ntime,ntimerange,nvar,ilevel,itimerange,ivar,indl,indt,itime,nx,ny real,allocatable :: vol2d(:,:) nx=vg6din%griddim%dim%nx ny=vg6din%griddim%dim%ny nlevel=count_distinct(roundlevel,back=.true.) ntime=size(vg6din%time) ntimerange=count_distinct(roundtimerange,back=.true.) nvar=size(vg6din%var) call init(vg6dout, vg6din%griddim, vg6din%time_definition, categoryappend="generated by vg6d_reduce") call volgrid6d_alloc(vg6dout, vg6din%griddim%dim, ntime, nlevel, ntimerange, nvar) if ( ASSOCIATED(vg6din%voldati) .or. optio_log(merge)) then call volgrid6d_alloc_vol(vg6dout,inivol=.true.,decode=.true.) allocate(vol2d(nx,ny)) else call volgrid6d_alloc_vol(vg6dout,inivol=.true.,decode=.false.) end if vg6dout%time=vg6din%time vg6dout%var=vg6din%var vg6dout%timerange=pack_distinct(roundtimerange,ntimerange,back=.true.) vg6dout%level=pack_distinct(roundlevel,nlevel,back=.true.) ! sort modified dimensions CALL sort(vg6dout%timerange) CALL sort(vg6dout%level) do ilevel=1,size(vg6din%level) indl=index(vg6dout%level,roundlevel(ilevel)) do itimerange=1,size(vg6din%timerange) indt=index(vg6dout%timerange,roundtimerange(itimerange)) do ivar=1, nvar do itime=1,ntime if ( ASSOCIATED(vg6din%voldati)) then vol2d=vg6din%voldati(:,:,ilevel,itime,itimerange,ivar) end if if (optio_log(merge)) then if ( .not. ASSOCIATED(vg6din%voldati)) then CALL grid_id_decode_data(vg6din%gaid(ilevel,itime,itimerange,ivar), vol2d) end if !! merge present data point by point where (.not. c_e(vg6dout%voldati(:,:,indl,itime,indt,ivar))) vg6dout%voldati(:,:,indl,itime,indt,ivar)=vol2d end where else if ( ASSOCIATED(vg6din%voldati)) then if (.not. any(c_e(vg6dout%voldati(:,:,indl,itime,indt,ivar))))then vg6dout%voldati(:,:,indl,itime,indt,ivar)=vol2d end if end if if (c_e(vg6din%gaid(ilevel,itime,itimerange,ivar)).and. .not. c_e(vg6dout%gaid(indl,itime,indt,ivar)))then call copy (vg6din%gaid(ilevel,itime,itimerange,ivar), vg6dout%gaid(indl,itime,indt,ivar)) end if end do end do end do end do if ( ASSOCIATED(vg6din%voldati) .or. optio_log(merge)) then deallocate(vol2d) end if end subroutine vg6d_reduce end module volgrid6d_class !>\example example_vg6d_3.f90 !!\brief Programma esempio semplice per gridinfo e volgrid6d. !! !! 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. !>\example example_vg6d_5.f90 !!\brief Programma trasformazione da volgrid6d a volgrid6d !! !! Legge grib da un file e li organizza in un vettore di strutture !! volgrid6d mettendoli a disposizione per eventuali elaborazioni; !! vengono poi riesportati a un file grib !>\example example_vg6d_8.f90 !! \brief Programma scrittura su file vettore di anagrafica !>\example example_vg6d_6.f90 !! \brief Programma trasformazione da volgrid6d a vol7d !>\example example_vg6d_7.f90 !! \brief Programma trasformazione da vol7d a volgrid6d !>\example example_vg6d_9.f90 !! \brief Example to create a grib editionNumber = 2 file from data generated in memory using a grib_api template