libsim Versione 7.2.0

◆ vg6d_reduce()

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

Reduce some dimensions (level and timerage).

You can pass a volume and specify duplicated levels and timeranges in roundlevel and roundtimerange; you get unique levels and timeranges in output. If there are data on equal levels or timeranges, the first var present (at least one point) will be taken (order is by icreasing var index). you can specify merge and if there are data on equal levels or timeranges will be merged POINT BY POINT with priority for the first data found ordered by icreasing var index (require to decode all the data) Data are decoded only if needed so the output should be with or without voldata allocated

Parametri
[in]vg6dininput volume
[out]vg6doutoutput volume
[in]roundlevelnew level list to use for rounding
[in]roundtimerangenew timerange list to use for rounding
[in]mergeif there are data on equal levels or timeranges will be merged POINT BY POINT with priority for the first data found ordered by icreasing var index (require to decode all the data)

Definizione alla linea 3605 del file volgrid6d_class.F90.

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