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