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