libsim Versione 7.1.11

◆ vg6d_rounding()

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

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

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

Parametri
[in]vg6dininput volume
[in]levelalmost equal level list
[in]timerangealmost equal timerange list
[in]mergeif there are data on more then one almost equal levels or timeranges will be merged POINT BY POINT with priority for the fird data found ordered by icreasing var index
[in]nostatprocdo not take in account statistical processing code in timerange and P2

Definizione alla linea 3558 del file volgrid6d_class.F90.

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

Generated with Doxygen.