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