libsim Versione 7.2.1
volgrid6d_class.F90
1! Copyright (C) 2010 ARPA-SIM <urpsim@smr.arpa.emr.it>
2! authors:
3! Davide Cesari <dcesari@arpa.emr.it>
4! Paolo Patruno <ppatruno@arpa.emr.it>
5
6! This program is free software; you can redistribute it and/or
7! modify it under the terms of the GNU General Public License as
8! published by the Free Software Foundation; either version 2 of
9! the License, or (at your option) any later version.
10
11! This program is distributed in the hope that it will be useful,
12! but WITHOUT ANY WARRANTY; without even the implied warranty of
13! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14! GNU General Public License for more details.
15
16! You should have received a copy of the GNU General Public License
17! along with this program. If not, see <http://www.gnu.org/licenses/>.
18#include "config.h"
19
31
52MODULE volgrid6d_class
57USE geo_proj_class
58USE grid_class
67!USE file_utilities
68#ifdef HAVE_DBALLE
70#endif
71IMPLICIT NONE
72
73character (len=255),parameter:: subcategory="volgrid6d_class"
74
76type volgrid6d
77 type(griddim_def) :: griddim
78 TYPE(datetime),pointer :: time(:)
79 TYPE(vol7d_timerange),pointer :: timerange(:)
80 TYPE(vol7d_level),pointer :: level(:)
81 TYPE(volgrid6d_var),pointer :: var(:)
82 TYPE(grid_id),POINTER :: gaid(:,:,:,:)
83 REAL,POINTER :: voldati(:,:,:,:,:,:)
84 integer :: time_definition
85 integer :: category = 0
86end type volgrid6d
87
89INTERFACE init
90 MODULE PROCEDURE volgrid6d_init
91END INTERFACE
92
95INTERFACE delete
96 MODULE PROCEDURE volgrid6d_delete, volgrid6dv_delete
97END INTERFACE
98
101INTERFACE import
102 MODULE PROCEDURE volgrid6d_read_from_file
103 MODULE PROCEDURE import_from_gridinfo, import_from_gridinfovv, &
104 volgrid6d_import_from_file
105END INTERFACE
106
109INTERFACE export
110 MODULE PROCEDURE volgrid6d_write_on_file
111 MODULE PROCEDURE export_to_gridinfo, export_to_gridinfov, export_to_gridinfovv,&
112 volgrid6d_export_to_file
113END INTERFACE
114
115! methods for computing transformations through an initialised
116! grid_transform object, probably too low level to be interfaced
117INTERFACE compute
118 MODULE PROCEDURE volgrid6d_transform_compute, volgrid6d_v7d_transform_compute,&
119 v7d_volgrid6d_transform_compute, v7d_v7d_transform_compute
120END INTERFACE
121
124INTERFACE transform
125 MODULE PROCEDURE volgrid6d_transform, volgrid6dv_transform,&
126 volgrid6d_v7d_transform, volgrid6dv_v7d_transform, v7d_volgrid6d_transform, &
127 v7d_v7d_transform
128END INTERFACE
129
130INTERFACE wind_rot
131 MODULE PROCEDURE vg6d_wind_rot
132END INTERFACE
133
134INTERFACE wind_unrot
135 MODULE PROCEDURE vg6d_wind_unrot
136END INTERFACE
137
140INTERFACE display
141 MODULE PROCEDURE display_volgrid6d,display_volgrid6dv
142END INTERFACE
143
155INTERFACE rounding
156 MODULE PROCEDURE vg6d_rounding, vg6dv_rounding
157END INTERFACE
158
159private
160
161PUBLIC volgrid6d,init,delete,export,import,compute,transform, &
162 wind_rot,wind_unrot,vg6d_c2a,display,volgrid6d_alloc,volgrid6d_alloc_vol, &
163 volgrid_get_vol_2d, volgrid_set_vol_2d, volgrid_get_vol_3d, volgrid_set_vol_3d
164PUBLIC rounding, vg6d_reduce
165
166CONTAINS
167
168
173SUBROUTINE volgrid6d_init(this, griddim, time_definition, categoryappend)
174TYPE(volgrid6d) :: this
175TYPE(griddim_def),OPTIONAL :: griddim
176INTEGER,INTENT(IN),OPTIONAL :: time_definition
177CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
178
179character(len=512) :: a_name
180
181if (present(categoryappend))then
182 call l4f_launcher(a_name,a_name_append=trim(subcategory)//"."//trim(categoryappend))
183else
184 call l4f_launcher(a_name,a_name_append=trim(subcategory))
185endif
186this%category=l4f_category_get(a_name)
187
188#ifdef DEBUG
189call l4f_category_log(this%category,l4f_debug,"init")
190#endif
191
192call init(this%griddim)
193
194if (present(griddim))then
195 call copy (griddim,this%griddim)
196end if
197
198CALL vol7d_var_features_init() ! initialise var features table once
199
200if(present(time_definition)) then
201 this%time_definition = time_definition
202else
203 this%time_definition = 0 !default to reference time
204end if
205
206nullify (this%time,this%timerange,this%level,this%var)
207nullify (this%gaid,this%voldati)
208
209END SUBROUTINE volgrid6d_init
210
211
222SUBROUTINE volgrid6d_alloc(this, dim, ntime, nlevel, ntimerange, nvar, ini)
223TYPE(volgrid6d),INTENT(inout) :: this
224TYPE(grid_dim),INTENT(in),OPTIONAL :: dim
225INTEGER,INTENT(in),OPTIONAL :: ntime
226INTEGER,INTENT(in),OPTIONAL :: nlevel
227INTEGER,INTENT(in),OPTIONAL :: ntimerange
228INTEGER,INTENT(in),OPTIONAL :: nvar
229LOGICAL,INTENT(in),OPTIONAL :: ini
230
231INTEGER :: i, stallo
232LOGICAL :: linit
233
234#ifdef DEBUG
235call l4f_category_log(this%category,l4f_debug,"alloc")
236#endif
237
238IF (PRESENT(ini)) THEN
239 linit = ini
240ELSE
241 linit = .false.
242ENDIF
243
244
245if (present(dim)) call copy (dim,this%griddim%dim)
246
247
248IF (PRESENT(ntime)) THEN
249 IF (ntime >= 0) THEN
250 IF (ASSOCIATED(this%time)) DEALLOCATE(this%time)
251#ifdef DEBUG
252 call l4f_category_log(this%category,l4f_debug,"alloc ntime "//to_char(ntime))
253#endif
254 ALLOCATE(this%time(ntime),stat=stallo)
255 if (stallo /=0)then
256 call l4f_category_log(this%category,l4f_fatal,"allocating memory")
257 CALL raise_fatal_error()
258 end if
259 IF (linit) THEN
260 DO i = 1, ntime
261 this%time(i) = datetime_miss
262! CALL init(this%time(i)) ! senza argomento inizializza a zero non missing
263 ! baco di datetime?
264 ENDDO
265 ENDIF
266 ENDIF
267ENDIF
268IF (PRESENT(nlevel)) THEN
269 IF (nlevel >= 0) THEN
270 IF (ASSOCIATED(this%level)) DEALLOCATE(this%level)
271#ifdef DEBUG
272 call l4f_category_log(this%category,l4f_debug,"alloc nlevel "//to_char(nlevel))
273#endif
274 ALLOCATE(this%level(nlevel),stat=stallo)
275 if (stallo /=0)then
276 call l4f_category_log(this%category,l4f_fatal,"allocating memory")
277 CALL raise_fatal_error()
278 end if
279 IF (linit) THEN
280 DO i = 1, nlevel
281 CALL init(this%level(i))
282 ENDDO
283 ENDIF
284 ENDIF
285ENDIF
286IF (PRESENT(ntimerange)) THEN
287 IF (ntimerange >= 0) THEN
288 IF (ASSOCIATED(this%timerange)) DEALLOCATE(this%timerange)
289#ifdef DEBUG
290 call l4f_category_log(this%category,l4f_debug,"alloc ntimerange "//to_char(ntimerange))
291#endif
292 ALLOCATE(this%timerange(ntimerange),stat=stallo)
293 if (stallo /=0)then
294 call l4f_category_log(this%category,l4f_fatal,"allocating memory")
295 CALL raise_fatal_error()
296 end if
297 IF (linit) THEN
298 DO i = 1, ntimerange
299 CALL init(this%timerange(i))
300 ENDDO
301 ENDIF
302 ENDIF
303ENDIF
304IF (PRESENT(nvar)) THEN
305 IF (nvar >= 0) THEN
306 IF (ASSOCIATED(this%var)) DEALLOCATE(this%var)
307#ifdef DEBUG
308 call l4f_category_log(this%category,l4f_debug,"alloc nvar "//to_char(nvar))
309#endif
310 ALLOCATE(this%var(nvar),stat=stallo)
311 if (stallo /=0)then
312 call l4f_category_log(this%category,l4f_fatal,"allocating memory")
313 CALL raise_fatal_error()
314 end if
315 IF (linit) THEN
316 DO i = 1, nvar
317 CALL init(this%var(i))
318 ENDDO
319 ENDIF
320 ENDIF
321ENDIF
322
323end SUBROUTINE volgrid6d_alloc
324
325
334SUBROUTINE volgrid6d_alloc_vol(this, ini, inivol, decode)
335TYPE(volgrid6d),INTENT(inout) :: this
336LOGICAL,INTENT(in),OPTIONAL :: ini
337LOGICAL,INTENT(in),OPTIONAL :: inivol
338LOGICAL,INTENT(in),OPTIONAL :: decode
339
340INTEGER :: stallo
341LOGICAL :: linivol
342
343#ifdef DEBUG
344call l4f_category_log(this%category,l4f_debug,"start alloc_vol")
345#endif
346
347IF (PRESENT(inivol)) THEN ! opposite condition, cannot use optio_log
348 linivol = inivol
349ELSE
350 linivol = .true.
351ENDIF
352
353IF (this%griddim%dim%nx > 0 .AND. this%griddim%dim%ny > 0) THEN
354! allocate dimension descriptors to a minimal size for having a
355! non-null volume
356 IF (.NOT. ASSOCIATED(this%var)) CALL volgrid6d_alloc(this, nvar=1, ini=ini)
357 IF (.NOT. ASSOCIATED(this%time)) CALL volgrid6d_alloc(this, ntime=1, ini=ini)
358 IF (.NOT. ASSOCIATED(this%level)) CALL volgrid6d_alloc(this, nlevel=1, ini=ini)
359 IF (.NOT. ASSOCIATED(this%timerange)) CALL volgrid6d_alloc(this, ntimerange=1, ini=ini)
360
361 IF (optio_log(decode)) THEN
362 IF (.NOT.ASSOCIATED(this%voldati)) THEN
363#ifdef DEBUG
364 CALL l4f_category_log(this%category,l4f_debug,"alloc voldati volume")
365#endif
366
367 ALLOCATE(this%voldati(this%griddim%dim%nx,this%griddim%dim%ny,&
368 SIZE(this%level), SIZE(this%time), &
369 SIZE(this%timerange), SIZE(this%var)),stat=stallo)
370 IF (stallo /= 0)THEN
371 CALL l4f_category_log(this%category,l4f_fatal,"allocating memory")
372 CALL raise_fatal_error()
373 ENDIF
374
375! this is not really needed if we can check other flags for a full
376! field missing values
377 IF (linivol) this%voldati = rmiss
378 this%voldati = rmiss
379 ENDIF
380 ENDIF
381
382 IF (.NOT.ASSOCIATED(this%gaid)) THEN
383#ifdef DEBUG
384 CALL l4f_category_log(this%category,l4f_debug,"alloc gaid volume")
385#endif
386 ALLOCATE(this%gaid(SIZE(this%level), SIZE(this%time), &
387 SIZE(this%timerange), SIZE(this%var)),stat=stallo)
388 IF (stallo /= 0)THEN
389 CALL l4f_category_log(this%category,l4f_fatal,"allocating memory")
390 CALL raise_fatal_error()
391 ENDIF
392
393 IF (linivol) THEN
394!!$ DO i=1 ,SIZE(this%gaid,1)
395!!$ DO ii=1 ,SIZE(this%gaid,2)
396!!$ DO iii=1 ,SIZE(this%gaid,3)
397!!$ DO iiii=1 ,SIZE(this%gaid,4)
398!!$ this%gaid(i,ii,iii,iiii) = grid_id_new() ! optimize?
399!!$ ENDDO
400!!$ ENDDO
401!!$ ENDDO
402!!$ ENDDO
403
404 this%gaid = grid_id_new()
405 ENDIF
406 ENDIF
407
408ELSE
409 CALL l4f_category_log(this%category,l4f_fatal,'volgrid6d_alloc_vol: &
410 &trying to allocate a volume with an invalid or unspecified horizontal grid')
411 CALL raise_fatal_error()
412ENDIF
413
414END SUBROUTINE volgrid6d_alloc_vol
415
416
430SUBROUTINE volgrid_get_vol_2d(this, ilevel, itime, itimerange, ivar, voldati)
431TYPE(volgrid6d),INTENT(in) :: this
432INTEGER,INTENT(in) :: ilevel
433INTEGER,INTENT(in) :: itime
434INTEGER,INTENT(in) :: itimerange
435INTEGER,INTENT(in) :: ivar
436REAL,POINTER :: voldati(:,:)
437
438IF (ASSOCIATED(this%voldati)) THEN
439 voldati => this%voldati(:,:,ilevel,itime,itimerange,ivar)
440 RETURN
441ELSE
442 IF (.NOT.ASSOCIATED(voldati)) THEN
443 ALLOCATE(voldati(this%griddim%dim%nx,this%griddim%dim%ny))
444 ENDIF
445 CALL grid_id_decode_data(this%gaid(ilevel,itime,itimerange,ivar), voldati)
446ENDIF
447
448END SUBROUTINE volgrid_get_vol_2d
449
450
464SUBROUTINE volgrid_get_vol_3d(this, itime, itimerange, ivar, voldati)
465TYPE(volgrid6d),INTENT(in) :: this
466INTEGER,INTENT(in) :: itime
467INTEGER,INTENT(in) :: itimerange
468INTEGER,INTENT(in) :: ivar
469REAL,POINTER :: voldati(:,:,:)
470
471INTEGER :: ilevel
472
473IF (ASSOCIATED(this%voldati)) THEN
474 voldati => this%voldati(:,:,:,itime,itimerange,ivar)
475 RETURN
476ELSE
477 IF (.NOT.ASSOCIATED(voldati)) THEN
478 ALLOCATE(voldati(this%griddim%dim%nx,this%griddim%dim%ny,SIZE(this%level)))
479 ENDIF
480 DO ilevel = 1, SIZE(this%level)
481 CALL grid_id_decode_data(this%gaid(ilevel,itime,itimerange,ivar), &
482 voldati(:,:,ilevel))
483 ENDDO
484ENDIF
485
486END SUBROUTINE volgrid_get_vol_3d
487
488
500SUBROUTINE volgrid_set_vol_2d(this, ilevel, itime, itimerange, ivar, voldati)
501TYPE(volgrid6d),INTENT(inout) :: this
502INTEGER,INTENT(in) :: ilevel
503INTEGER,INTENT(in) :: itime
504INTEGER,INTENT(in) :: itimerange
505INTEGER,INTENT(in) :: ivar
506REAL,INTENT(in) :: voldati(:,:)
507
508IF (ASSOCIATED(this%voldati)) THEN
509 RETURN
510ELSE
511 CALL grid_id_encode_data(this%gaid(ilevel,itime,itimerange,ivar), voldati)
512ENDIF
513
514END SUBROUTINE volgrid_set_vol_2d
515
516
528SUBROUTINE volgrid_set_vol_3d(this, itime, itimerange, ivar, voldati)
529TYPE(volgrid6d),INTENT(inout) :: this
530INTEGER,INTENT(in) :: itime
531INTEGER,INTENT(in) :: itimerange
532INTEGER,INTENT(in) :: ivar
533REAL,INTENT(in) :: voldati(:,:,:)
534
535INTEGER :: ilevel
536
537IF (ASSOCIATED(this%voldati)) THEN
538 RETURN
539ELSE
540 DO ilevel = 1, SIZE(this%level)
541 CALL grid_id_encode_data(this%gaid(ilevel,itime,itimerange,ivar), &
542 voldati(:,:,ilevel))
543 ENDDO
544ENDIF
545
546END SUBROUTINE volgrid_set_vol_3d
547
548
552SUBROUTINE volgrid6d_delete(this)
553TYPE(volgrid6d),INTENT(inout) :: this
554
555INTEGER :: i, ii, iii, iiii
556
557#ifdef DEBUG
558call l4f_category_log(this%category,l4f_debug,"delete")
559#endif
560
561if (associated(this%gaid))then
562
563 DO i=1 ,SIZE(this%gaid,1)
564 DO ii=1 ,SIZE(this%gaid,2)
565 DO iii=1 ,SIZE(this%gaid,3)
566 DO iiii=1 ,SIZE(this%gaid,4)
567 CALL delete(this%gaid(i,ii,iii,iiii))
568 ENDDO
569 ENDDO
570 ENDDO
571 ENDDO
572 DEALLOCATE(this%gaid)
573
574end if
575
576call delete(this%griddim)
577
578! call delete(this%time)
579! call delete(this%timerange)
580! call delete(this%level)
581! call delete(this%var)
582
583if (associated( this%time )) deallocate(this%time)
584if (associated( this%timerange )) deallocate(this%timerange)
585if (associated( this%level )) deallocate(this%level)
586if (associated( this%var )) deallocate(this%var)
587
588if (associated(this%voldati))deallocate(this%voldati)
589
590
591 !chiudo il logger
592call l4f_category_delete(this%category)
593
594END SUBROUTINE volgrid6d_delete
595
596
606subroutine volgrid6d_write_on_file (this,unit,description,filename,filename_auto)
607
608TYPE(volgrid6d),INTENT(IN) :: this
609integer,optional,intent(inout) :: unit
610character(len=*),intent(in),optional :: filename
611character(len=*),intent(out),optional :: filename_auto
612character(len=*),INTENT(IN),optional :: description
613
614integer :: lunit
615character(len=254) :: ldescription,arg,lfilename
616integer :: ntime, ntimerange, nlevel, nvar
617integer :: tarray(8)
618logical :: opened,exist
619
620#ifdef DEBUG
621call l4f_category_log(this%category,l4f_debug,"write on file")
622#endif
623
624ntime=0
625ntimerange=0
626nlevel=0
627nvar=0
628
629!call idate(im,id,iy)
630call date_and_time(values=tarray)
631call getarg(0,arg)
632
633if (present(description))then
634 ldescription=description
635else
636 ldescription="Volgrid6d generated by: "//trim(arg)
637end if
638
639if (.not. present(unit))then
640 lunit=getunit()
641else
642 if (unit==0)then
643 lunit=getunit()
644 unit=lunit
645 else
646 lunit=unit
647 end if
648end if
649
650lfilename=trim(arg)//".vg6d"
651if (index(arg,'/',back=.true.) > 0) lfilename=lfilename(index(arg,'/',back=.true.)+1 : )
653if (present(filename))then
654 if (filename /= "")then
655 lfilename=filename
656 end if
657end if
658
659if (present(filename_auto))filename_auto=lfilename
660
661
662inquire(unit=lunit,opened=opened)
663if (.not. opened) then
664 inquire(file=lfilename,exist=exist)
665 if (exist) CALL raise_error('file exist; cannot open new file')
666 if (.not.exist) open (unit=lunit,file=lfilename,form="UNFORMATTED")
667 !print *, "opened: ",lfilename
668end if
669
670if (associated(this%time)) ntime=size(this%time)
671if (associated(this%timerange)) ntimerange=size(this%timerange)
672if (associated(this%level)) nlevel=size(this%level)
673if (associated(this%var)) nvar=size(this%var)
674
675
676write(unit=lunit)ldescription
677write(unit=lunit)tarray
678
679call write_unit( this%griddim,lunit)
680write(unit=lunit) ntime, ntimerange, nlevel, nvar
681
682!! prime 4 dimensioni
683if (associated(this%time)) call write_unit(this%time, lunit)
684if (associated(this%level)) write(unit=lunit)this%level
685if (associated(this%timerange)) write(unit=lunit)this%timerange
686if (associated(this%var)) write(unit=lunit)this%var
687
689!! Volumi di valori dati
690
691if (associated(this%voldati)) write(unit=lunit)this%voldati
692
693if (.not. present(unit)) close(unit=lunit)
694
695end subroutine volgrid6d_write_on_file
696
697
704subroutine volgrid6d_read_from_file (this,unit,filename,description,tarray,filename_auto)
705
706TYPE(volgrid6d),INTENT(OUT) :: this
707integer,intent(inout),optional :: unit
708character(len=*),INTENT(in),optional :: filename
709character(len=*),intent(out),optional :: filename_auto
710character(len=*),INTENT(out),optional :: description
711integer,intent(out),optional :: tarray(8)
712
713integer :: ntime, ntimerange, nlevel, nvar
714
715character(len=254) :: ldescription,lfilename,arg
716integer :: ltarray(8),lunit
717logical :: opened,exist
718
719#ifdef DEBUG
720call l4f_category_log(this%category,l4f_debug,"read from file")
721#endif
722
723call getarg(0,arg)
724
725if (.not. present(unit))then
726 lunit=getunit()
727else
728 if (unit==0)then
729 lunit=getunit()
730 unit=lunit
731 else
732 lunit=unit
733 end if
734end if
735
736lfilename=trim(arg)//".vg6d"
737if (index(arg,'/',back=.true.) > 0) lfilename=lfilename(index(arg,'/',back=.true.)+1 : )
738
739if (present(filename))then
740 if (filename /= "")then
741 lfilename=filename
742 end if
743end if
744
745if (present(filename_auto))filename_auto=lfilename
746
747
748inquire(unit=lunit,opened=opened)
749if (.not. opened) then
750 inquire(file=lfilename,exist=exist)
751 IF (.NOT. exist) CALL raise_fatal_error('file '//trim(lfilename)//' does not exist, cannot open')
752 open (unit=lunit,file=lfilename,form="UNFORMATTED")
753end if
754
755read(unit=lunit)ldescription
756read(unit=lunit)ltarray
757
758call l4f_log(l4f_info,"Info: reading volgrid6d from file: "//trim(lfilename))
759call l4f_log(l4f_info,"Info: description: "//trim(ldescription))
760!call l4f_log("Info: written on ",ltarray)
761
762if (present(description))description=ldescription
763if (present(tarray))tarray=ltarray
764
765
766call read_unit( this%griddim,lunit)
767read(unit=lunit) ntime, ntimerange, nlevel, nvar
768
769
770call volgrid6d_alloc (this, &
771 ntime=ntime, ntimerange=ntimerange, nlevel=nlevel, nvar=nvar)
772
773call volgrid6d_alloc_vol (this)
774
775if (associated(this%time)) call read_unit(this%time, lunit)
776if (associated(this%level)) read(unit=lunit)this%level
777if (associated(this%timerange)) read(unit=lunit)this%timerange
778if (associated(this%var)) read(unit=lunit)this%var
779
780
781!! Volumi di valori
782
783if (associated(this%voldati)) read(unit=lunit)this%voldati
784
785if (.not. present(unit)) close(unit=lunit)
786
787end subroutine volgrid6d_read_from_file
788
789
809SUBROUTINE import_from_gridinfo(this, gridinfo, force, dup_mode, clone, &
810 isanavar)
811TYPE(volgrid6d),INTENT(inout) :: this
812TYPE(gridinfo_def),INTENT(in) :: gridinfo
813LOGICAL,INTENT(in),OPTIONAL :: force
814INTEGER,INTENT(in),OPTIONAL :: dup_mode
815LOGICAL , INTENT(in),OPTIONAL :: clone
816LOGICAL,INTENT(IN),OPTIONAL :: isanavar
817
818CHARACTER(len=255) :: type
819INTEGER :: itime0, itimerange0, itime1, itimerange1, itime, itimerange, &
820 ilevel, ivar, ldup_mode
821LOGICAL :: dup
822TYPE(datetime) :: correctedtime
823REAL,ALLOCATABLE :: tmpgrid(:,:)
824
825IF (PRESENT(dup_mode)) THEN
826 ldup_mode = dup_mode
827ELSE
828 ldup_mode = 0
829ENDIF
830
831call get_val(this%griddim,proj_type=type)
832
833#ifdef DEBUG
834call l4f_category_log(this%category,l4f_debug,"import_from_gridinfo: "//trim(type))
835#endif
836
837if (.not. c_e(type))then
838 call copy(gridinfo%griddim, this%griddim)
839! ho gia` fatto init, altrimenti non potrei fare la get_val(this%griddim)
840! per cui meglio non ripetere
841! call init(this,gridinfo%griddim,categoryappend)
842 CALL volgrid6d_alloc_vol(this, ini=.true.) ! decode?
843
844else if (.not. (this%griddim == gridinfo%griddim ))then
845
846 CALL l4f_category_log(this%category,l4f_error, &
847 "volgrid and gridinfo grid type or size are different, gridinfo rejected")
848 CALL raise_error()
849 RETURN
850
851end if
852
853! Cerco gli indici del campo da inserire, se non trovo metto nel primo missing
854ilevel = index(this%level, gridinfo%level)
855IF (ilevel == 0 .AND. optio_log(force)) THEN
856 ilevel = index(this%level, vol7d_level_miss)
857 IF (ilevel /= 0) this%level(ilevel) = gridinfo%level
858ENDIF
859
860IF (ilevel == 0) THEN
861 CALL l4f_category_log(this%category,l4f_error, &
862 "volgrid6d: level not valid for volume, gridinfo rejected")
863 CALL raise_error()
864 RETURN
865ENDIF
866
867IF (optio_log(isanavar)) THEN ! assign to all times and timeranges
868 itime0 = 1
869 itime1 = SIZE(this%time)
870 itimerange0 = 1
871 itimerange1 = SIZE(this%timerange)
872ELSE ! usual case
873 correctedtime = gridinfo%time
874 IF (this%time_definition == 1) correctedtime = correctedtime + &
875 timedelta_new(sec=gridinfo%timerange%p1)
876 itime0 = index(this%time, correctedtime)
877 IF (itime0 == 0 .AND. optio_log(force)) THEN
878 itime0 = index(this%time, datetime_miss)
879 IF (itime0 /= 0) this%time(itime0) = correctedtime
880 ENDIF
881 IF (itime0 == 0) THEN
882 CALL l4f_category_log(this%category,l4f_error, &
883 "volgrid6d: time not valid for volume, gridinfo rejected")
884 CALL raise_error()
885 RETURN
886 ENDIF
887 itime1 = itime0
888
889 itimerange0 = index(this%timerange,gridinfo%timerange)
890 IF (itimerange0 == 0 .AND. optio_log(force)) THEN
891 itimerange0 = index(this%timerange, vol7d_timerange_miss)
892 IF (itimerange0 /= 0) this%timerange(itimerange0) = gridinfo%timerange
893 ENDIF
894 IF (itimerange0 == 0) THEN
895 CALL l4f_category_log(this%category,l4f_error, &
896 "volgrid6d: timerange not valid for volume, gridinfo rejected")
897 CALL raise_error()
898 RETURN
899 ENDIF
900 itimerange1 = itimerange0
901ENDIF
902
903ivar = index(this%var, gridinfo%var)
904IF (ivar == 0 .AND. optio_log(force)) THEN
905 ivar = index(this%var, volgrid6d_var_miss)
906 IF (ivar /= 0) this%var(ivar) = gridinfo%var
907ENDIF
908IF (ivar == 0) THEN
909 CALL l4f_category_log(this%category,l4f_error, &
910 "volgrid6d: var not valid for volume, gridinfo rejected")
911 CALL raise_error()
912 RETURN
913ENDIF
914
915DO itimerange = itimerange0, itimerange1
916 DO itime = itime0, itime1
917 IF (ASSOCIATED(this%gaid)) THEN
918 dup = .false.
919 IF (c_e(this%gaid(ilevel,itime,itimerange,ivar))) THEN
920 dup = .true.
921 CALL l4f_category_log(this%category,l4f_warn,"gaid exist: grib duplicated")
922! avoid memory leaks
923 IF (optio_log(clone)) CALL delete(this%gaid(ilevel,itime,itimerange,ivar))
924 ENDIF
925
926 IF (optio_log(clone)) THEN
927 CALL copy(gridinfo%gaid, this%gaid(ilevel,itime,itimerange,ivar))
928#ifdef DEBUG
929 CALL l4f_category_log(this%category,l4f_debug,"cloning to a new gaid")
930#endif
931 ELSE
932 this%gaid(ilevel,itime,itimerange,ivar) = gridinfo%gaid
933 ENDIF
934
935 IF (ASSOCIATED(this%voldati))THEN
936 IF (.NOT.dup .OR. ldup_mode == 0) THEN
937 this%voldati(:,:,ilevel,itime,itimerange,ivar) = decode_gridinfo(gridinfo)
938 ELSE IF (ldup_mode == 1) THEN
939 tmpgrid = decode_gridinfo(gridinfo) ! f2003 automatic allocation
940 WHERE(c_e(tmpgrid))
941 this%voldati(:,:,ilevel,itime,itimerange,ivar) = tmpgrid(:,:)
942 END WHERE
943 ELSE IF (ldup_mode == 2) THEN
944 WHERE(.NOT.c_e(this%voldati(:,:,ilevel,itime,itimerange,ivar)))
945 this%voldati(:,:,ilevel,itime,itimerange,ivar) = decode_gridinfo(gridinfo)
946 END WHERE
947 ENDIF
948 ENDIF
949
950 ELSE
951 CALL l4f_category_log(this%category,l4f_error, &
952 "gaid not allocated, you probably need to call volgrid6d_alloc_vol first")
953 CALL raise_error()
954 RETURN
955 ENDIF
956 ENDDO
957ENDDO
958
959
960END SUBROUTINE import_from_gridinfo
961
962
967SUBROUTINE export_to_gridinfo(this, gridinfo, itime, itimerange, ilevel, ivar, &
968 gaid_template, clone)
969TYPE(volgrid6d),INTENT(in) :: this
970TYPE(gridinfo_def),INTENT(inout) :: gridinfo
971INTEGER :: itime
972INTEGER :: itimerange
973INTEGER :: ilevel
974INTEGER :: ivar
975TYPE(grid_id),INTENT(in),OPTIONAL :: gaid_template
976LOGICAL,INTENT(in),OPTIONAL :: clone
977
978TYPE(grid_id) :: gaid
979LOGICAL :: usetemplate
980REAL,POINTER :: voldati(:,:)
981TYPE(datetime) :: correctedtime
982
983#ifdef DEBUG
984CALL l4f_category_log(this%category,l4f_debug,"export_to_gridinfo")
985#endif
986
987IF (.NOT.c_e(this%gaid(ilevel,itime,itimerange,ivar))) THEN
988#ifdef DEBUG
989 CALL l4f_category_log(this%category,l4f_debug,"empty gaid found, skipping export")
990#endif
991 RETURN
992ENDIF
993
994usetemplate = .false.
995IF (PRESENT(gaid_template)) THEN
996 CALL copy(gaid_template, gaid)
997#ifdef DEBUG
998 CALL l4f_category_log(this%category,l4f_debug,"template cloned to a new gaid")
999#endif
1000 usetemplate = c_e(gaid)
1001ENDIF
1002
1003IF (.NOT.usetemplate) THEN
1004 IF (optio_log(clone)) THEN
1005 CALL copy(this%gaid(ilevel,itime,itimerange,ivar), gaid)
1006#ifdef DEBUG
1007 CALL l4f_category_log(this%category,l4f_debug,"original gaid cloned to a new one")
1008#endif
1009 ELSE
1010 gaid = this%gaid(ilevel,itime,itimerange,ivar)
1011 ENDIF
1012ENDIF
1013
1014IF (this%time_definition == 1) THEN
1015 correctedtime = this%time(itime) - &
1016 timedelta_new(sec=this%timerange(itimerange)%p1)
1017ELSE
1018 correctedtime = this%time(itime)
1019ENDIF
1020
1021CALL init(gridinfo,gaid, this%griddim, correctedtime, this%timerange(itimerange), &
1022 this%level(ilevel), this%var(ivar))
1023
1024! reset the gridinfo, bad but necessary at this point for encoding the field
1025CALL export(gridinfo%griddim, gridinfo%gaid)
1026! encode the field
1027IF (ASSOCIATED(this%voldati)) THEN
1028 CALL encode_gridinfo(gridinfo, this%voldati(:,:,ilevel,itime,itimerange,ivar))
1029ELSE IF (usetemplate) THEN ! field must be forced into template in this case
1030 NULLIFY(voldati)
1031 CALL volgrid_get_vol_2d(this, ilevel, itime, itimerange, ivar, voldati)
1032 CALL encode_gridinfo(gridinfo, voldati)
1033 DEALLOCATE(voldati)
1034ENDIF
1035
1036END SUBROUTINE export_to_gridinfo
1037
1038
1056SUBROUTINE import_from_gridinfovv(this, gridinfov, dup_mode, clone, decode, &
1057 time_definition, anavar, categoryappend)
1058TYPE(volgrid6d),POINTER :: this(:)
1059TYPE(arrayof_gridinfo),INTENT(in) :: gridinfov
1060INTEGER,INTENT(in),OPTIONAL :: dup_mode
1061LOGICAL , INTENT(in),OPTIONAL :: clone
1062LOGICAL,INTENT(in),OPTIONAL :: decode
1063INTEGER,INTENT(IN),OPTIONAL :: time_definition
1064CHARACTER(len=*),INTENT(IN),OPTIONAL :: anavar(:)
1065CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
1066
1067INTEGER :: i, j, stallo
1068INTEGER :: ngrid, ntime, ntimerange, nlevel, nvar
1069INTEGER :: category
1070CHARACTER(len=512) :: a_name
1071TYPE(datetime),ALLOCATABLE :: correctedtime(:)
1072LOGICAL,ALLOCATABLE :: isanavar(:)
1073TYPE(vol7d_var) :: lvar
1074
1075! category temporanea (altrimenti non possiamo loggare)
1076if (present(categoryappend))then
1077 call l4f_launcher(a_name,a_name_append=trim(subcategory)//"."//trim(categoryappend))
1078else
1079 call l4f_launcher(a_name,a_name_append=trim(subcategory))
1080endif
1081category=l4f_category_get(a_name)
1082
1083#ifdef DEBUG
1084call l4f_category_log(category,l4f_debug,"start import_from_gridinfovv")
1085#endif
1086
1087ngrid=count_distinct(gridinfov%array(1:gridinfov%arraysize)%griddim,back=.true.)
1088CALL l4f_category_log(category,l4f_info, t2c(ngrid)// &
1089 ' different grid definition(s) found in input data')
1090
1091ALLOCATE(this(ngrid),stat=stallo)
1092IF (stallo /= 0)THEN
1093 CALL l4f_category_log(category,l4f_fatal,"allocating memory")
1094 CALL raise_fatal_error()
1095ENDIF
1096DO i = 1, ngrid
1097 IF (PRESENT(categoryappend))THEN
1098 CALL init(this(i), time_definition=time_definition, categoryappend=trim(categoryappend)//"-vol"//t2c(i))
1099 ELSE
1100 CALL init(this(i), time_definition=time_definition, categoryappend="vol"//t2c(i))
1101 ENDIF
1102ENDDO
1103
1104this(:)%griddim=pack_distinct(gridinfov%array(1:gridinfov%arraysize)%griddim, &
1105 ngrid, back=.true.)
1106
1107! mark elements as ana variables (time-independent)
1108ALLOCATE(isanavar(gridinfov%arraysize))
1109isanavar(:) = .false.
1110IF (PRESENT(anavar)) THEN
1111 DO i = 1, gridinfov%arraysize
1112 DO j = 1, SIZE(anavar)
1113 lvar = convert(gridinfov%array(i)%var)
1114 IF (lvar%btable == anavar(j)) THEN
1115 isanavar(i) = .true.
1116 EXIT
1117 ENDIF
1118 ENDDO
1119 ENDDO
1120 CALL l4f_category_log(category,l4f_info,t2c(count(isanavar))//'/'// &
1121 t2c(gridinfov%arraysize)//' constant-data messages found in input data')
1122ENDIF
1123
1124! create time corrected for time_definition
1125ALLOCATE(correctedtime(gridinfov%arraysize))
1126correctedtime(:) = gridinfov%array(1:gridinfov%arraysize)%time
1127IF (PRESENT(time_definition)) THEN
1128 IF (time_definition == 1) THEN
1129 DO i = 1, gridinfov%arraysize
1130 correctedtime(i) = correctedtime(i) + &
1131 timedelta_new(sec=gridinfov%array(i)%timerange%p1)
1132 ENDDO
1133 ENDIF
1134ENDIF
1135
1136DO i = 1, ngrid
1137 IF (PRESENT(anavar)) THEN
1138 j = count((this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
1139 .AND. .NOT.isanavar(:))
1140 IF (j <= 0) THEN
1141 CALL l4f_category_log(category, l4f_fatal, 'grid n.'//t2c(i)// &
1142 ' has only constant data, this is not allowed')
1143 CALL l4f_category_log(category, l4f_fatal, 'please check anavar argument')
1144 CALL raise_fatal_error()
1145 ENDIF
1146 ENDIF
1147 ntime = count_distinct(correctedtime, &
1148 mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
1149 .AND. .NOT.isanavar(:), back=.true.)
1150 ntimerange = count_distinct(gridinfov%array(1:gridinfov%arraysize)%timerange, &
1151 mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
1152 .AND. .NOT.isanavar(:), back=.true.)
1153 nlevel = count_distinct(gridinfov%array(1:gridinfov%arraysize)%level, &
1154 mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim), &
1155 back=.true.)
1156 nvar = count_distinct(gridinfov%array(1:gridinfov%arraysize)%var, &
1157 mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim), &
1158 back=.true.)
1159
1160#ifdef DEBUG
1161 CALL l4f_category_log(this(i)%category,l4f_debug,"alloc volgrid6d index: "//t2c(i))
1162#endif
1163
1164 CALL volgrid6d_alloc(this(i),this(i)%griddim%dim,ntime=ntime, &
1165 ntimerange=ntimerange,nlevel=nlevel,nvar=nvar)
1166
1167 this(i)%time = pack_distinct(correctedtime, ntime, &
1168 mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
1169 .AND. .NOT.isanavar(:), back=.true.)
1170 CALL sort(this(i)%time)
1171
1172 this(i)%timerange = pack_distinct(gridinfov%array( &
1173 1:gridinfov%arraysize)%timerange, ntimerange, &
1174 mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim) &
1175 .AND. .NOT.isanavar(:), back=.true.)
1176 CALL sort(this(i)%timerange)
1177
1178 this(i)%level=pack_distinct(gridinfov%array(1:gridinfov%arraysize)%level, &
1179 nlevel,mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim), &
1180 back=.true.)
1181 CALL sort(this(i)%level)
1182
1183 this(i)%var=pack_distinct(gridinfov%array(1:gridinfov%arraysize)%var, nvar, &
1184 mask=(this(i)%griddim == gridinfov%array(1:gridinfov%arraysize)%griddim), &
1185 back=.true.)
1186
1187#ifdef DEBUG
1188 CALL l4f_category_log(this(i)%category,l4f_debug,"alloc_vol volgrid6d index: "//t2c(i))
1189#endif
1190 CALL volgrid6d_alloc_vol(this(i), decode=decode)
1191
1192ENDDO
1193
1194DEALLOCATE(correctedtime)
1195
1196DO i = 1, gridinfov%arraysize
1197
1198#ifdef DEBUG
1199 CALL l4f_category_log(category,l4f_debug,"import from gridinfov index: "//t2c(i))
1200 CALL l4f_category_log(category,l4f_info, &
1201 "to volgrid6d index: "//t2c(index(this%griddim, gridinfov%array(i)%griddim)))
1202#endif
1203
1204 CALL import(this(index(this%griddim, gridinfov%array(i)%griddim)), &
1205 gridinfov%array(i), dup_mode=dup_mode, clone=clone, isanavar=isanavar(i))
1206
1207ENDDO
1208
1209!chiudo il logger temporaneo
1210CALL l4f_category_delete(category)
1211
1212END SUBROUTINE import_from_gridinfovv
1213
1214
1220SUBROUTINE export_to_gridinfov(this, gridinfov, gaid_template, clone)
1221TYPE(volgrid6d),INTENT(inout) :: this
1222TYPE(arrayof_gridinfo),INTENT(inout) :: gridinfov
1223TYPE(grid_id),INTENT(in),OPTIONAL :: gaid_template
1224LOGICAL,INTENT(in),OPTIONAL :: clone
1225
1226INTEGER :: i ,itime, itimerange, ilevel, ivar
1227INTEGER :: ntime, ntimerange, nlevel, nvar
1228TYPE(gridinfo_def) :: gridinfol
1229
1230#ifdef DEBUG
1231CALL l4f_category_log(this%category,l4f_debug,"start export_to_gridinfov")
1232#endif
1233
1234! this is necessary in order not to repeatedly and uselessly copy the
1235! same array of coordinates for each 2d grid during export, the
1236! side-effect is that the computed projection in this is lost
1237CALL dealloc(this%griddim%dim)
1238
1239i=0
1240ntime=size(this%time)
1241ntimerange=size(this%timerange)
1242nlevel=size(this%level)
1243nvar=size(this%var)
1245DO itime=1,ntime
1246 DO itimerange=1,ntimerange
1247 DO ilevel=1,nlevel
1248 DO ivar=1,nvar
1249
1250 CALL init(gridinfol)
1251 CALL export(this, gridinfol, itime, itimerange, ilevel, ivar, &
1252 gaid_template=gaid_template, clone=clone)
1253 IF (c_e(gridinfol%gaid)) THEN
1254 CALL insert(gridinfov, gridinfol)
1255 ELSE
1256 CALL delete(gridinfol)
1257 ENDIF
1258
1259 ENDDO
1260 ENDDO
1261 ENDDO
1262ENDDO
1263
1264END SUBROUTINE export_to_gridinfov
1265
1266
1272SUBROUTINE export_to_gridinfovv(this, gridinfov, gaid_template, clone)
1273!, &
1274! categoryappend)
1275TYPE(volgrid6d),INTENT(inout) :: this(:)
1276TYPE(arrayof_gridinfo),INTENT(inout) :: gridinfov
1277TYPE(grid_id),INTENT(in),OPTIONAL :: gaid_template
1278LOGICAL,INTENT(in),OPTIONAL :: clone
1279
1280INTEGER :: i
1281
1282DO i = 1, SIZE(this)
1283#ifdef DEBUG
1284 CALL l4f_category_log(this(i)%category,l4f_debug, &
1285 "export_to_gridinfovv grid index: "//t2c(i))
1286#endif
1287 CALL export(this(i), gridinfov, gaid_template=gaid_template, clone=clone)
1288ENDDO
1289
1290END SUBROUTINE export_to_gridinfovv
1291
1292
1302SUBROUTINE volgrid6d_import_from_file(this, filename, dup_mode, decode, &
1303 time_definition, anavar, categoryappend)
1304TYPE(volgrid6d),POINTER :: this(:)
1305CHARACTER(len=*),INTENT(in) :: filename
1306INTEGER,INTENT(in),OPTIONAL :: dup_mode
1307LOGICAL,INTENT(in),OPTIONAL :: decode
1308INTEGER,INTENT(IN),OPTIONAL :: time_definition
1309CHARACTER(len=*),INTENT(IN),OPTIONAL :: anavar(:)
1310character(len=*),INTENT(in),OPTIONAL :: categoryappend
1311
1312TYPE(arrayof_gridinfo) :: gridinfo
1313INTEGER :: category
1314CHARACTER(len=512) :: a_name
1315
1316NULLIFY(this)
1317
1318IF (PRESENT(categoryappend))THEN
1319 CALL l4f_launcher(a_name,a_name_append= &
1320 trim(subcategory)//"."//trim(categoryappend))
1321ELSE
1322 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
1323ENDIF
1324category=l4f_category_get(a_name)
1325
1326CALL import(gridinfo, filename=filename, categoryappend=categoryappend)
1327
1328IF (gridinfo%arraysize > 0) THEN
1329
1330 CALL import(this, gridinfo, dup_mode=dup_mode, clone=.true., decode=decode, &
1331 time_definition=time_definition, anavar=anavar, &
1332 categoryappend=categoryappend)
1333
1334 CALL l4f_category_log(category,l4f_info,"deleting gridinfo")
1335 CALL delete(gridinfo)
1336
1337ELSE
1338 CALL l4f_category_log(category,l4f_info,"file does not contain gridded data")
1339ENDIF
1340
1341! close logger
1342CALL l4f_category_delete(category)
1343
1344END SUBROUTINE volgrid6d_import_from_file
1345
1346
1354SUBROUTINE volgrid6d_export_to_file(this, filename, gaid_template, categoryappend)
1355TYPE(volgrid6d) :: this(:)
1356CHARACTER(len=*),INTENT(in) :: filename
1357TYPE(grid_id),INTENT(in),OPTIONAL :: gaid_template
1358character(len=*),INTENT(in),OPTIONAL :: categoryappend
1359
1360TYPE(arrayof_gridinfo) :: gridinfo
1361INTEGER :: category
1362CHARACTER(len=512) :: a_name
1363
1364IF (PRESENT(categoryappend)) THEN
1365 CALL l4f_launcher(a_name,a_name_append=trim(subcategory)//"."//trim(categoryappend))
1366ELSE
1367 CALL l4f_launcher(a_name,a_name_append=trim(subcategory))
1368ENDIF
1369category=l4f_category_get(a_name)
1370
1371#ifdef DEBUG
1372CALL l4f_category_log(category,l4f_debug,"start export to file")
1373#endif
1374
1375CALL l4f_category_log(category,l4f_info,"writing volgrid6d to grib file: "//trim(filename))
1376
1377!IF (ASSOCIATED(this)) THEN
1378 CALL export(this, gridinfo, gaid_template=gaid_template, clone=.true.)
1379 IF (gridinfo%arraysize > 0) THEN
1380 CALL export(gridinfo, filename)
1381 CALL delete(gridinfo)
1382 ENDIF
1383!ELSE
1384! CALL l4f_category_log(category,L4F_INFO,"volume volgrid6d is not associated")
1385!ENDIF
1386
1387! close logger
1388CALL l4f_category_delete(category)
1389
1390END SUBROUTINE volgrid6d_export_to_file
1391
1392
1396SUBROUTINE volgrid6dv_delete(this)
1397TYPE(volgrid6d),POINTER :: this(:)
1398
1399INTEGER :: i
1400
1401IF (ASSOCIATED(this)) THEN
1402 DO i = 1, SIZE(this)
1403#ifdef DEBUG
1404 CALL l4f_category_log(this(i)%category,l4f_debug, &
1405 "delete volgrid6d vector index: "//trim(to_char(i)))
1406#endif
1407 CALL delete(this(i))
1408 ENDDO
1409 DEALLOCATE(this)
1410ENDIF
1411
1412END SUBROUTINE volgrid6dv_delete
1413
1414
1415! Internal method for performing grid to grid computations
1416SUBROUTINE volgrid6d_transform_compute(this, volgrid6d_in, volgrid6d_out, &
1417 lev_out, var_coord_vol, clone)
1418TYPE(grid_transform),INTENT(in) :: this ! oggetto di trasformazione per il grigliato
1419type(volgrid6d), INTENT(in) :: volgrid6d_in ! oggetto da trasformare
1420type(volgrid6d), INTENT(inout) :: volgrid6d_out ! oggetto trasformato; deve essere completo (init, alloc, alloc_vol)
1421TYPE(vol7d_level),INTENT(in),OPTIONAL :: lev_out(:) ! vol7d_level object defining target vertical grid, for vertical interpolations
1422INTEGER,INTENT(in),OPTIONAL :: var_coord_vol ! index of variable defining vertical coordinate values in input volume
1423LOGICAL,INTENT(in),OPTIONAL :: clone ! se fornito e \c .TRUE., clona i gaid da volgrid6d_in a volgrid6d_out
1424
1425INTEGER :: ntime, ntimerange, inlevel, onlevel, nvar, &
1426 itime, itimerange, ilevel, ivar, levshift, levused, lvar_coord_vol, spos
1427REAL,POINTER :: voldatiin(:,:,:), voldatiout(:,:,:), coord_3d_in(:,:,:)
1428TYPE(vol7d_level) :: output_levtype
1429
1430
1431#ifdef DEBUG
1432call l4f_category_log(volgrid6d_in%category,l4f_debug,"start volgrid6d_transform_compute")
1433#endif
1434
1435ntime=0
1436ntimerange=0
1437inlevel=0
1438onlevel=0
1439nvar=0
1440lvar_coord_vol = optio_i(var_coord_vol)
1441
1442if (associated(volgrid6d_in%time))then
1443 ntime=size(volgrid6d_in%time)
1444 volgrid6d_out%time=volgrid6d_in%time
1445end if
1446
1447if (associated(volgrid6d_in%timerange))then
1448 ntimerange=size(volgrid6d_in%timerange)
1449 volgrid6d_out%timerange=volgrid6d_in%timerange
1450end if
1451
1452IF (ASSOCIATED(volgrid6d_in%level))THEN
1453 inlevel=SIZE(volgrid6d_in%level)
1454ENDIF
1455IF (PRESENT(lev_out)) THEN
1456 onlevel=SIZE(lev_out)
1457 volgrid6d_out%level=lev_out
1458ELSE IF (ASSOCIATED(volgrid6d_in%level))THEN
1459 onlevel=SIZE(volgrid6d_in%level)
1460 volgrid6d_out%level=volgrid6d_in%level
1461ENDIF
1462
1463if (associated(volgrid6d_in%var))then
1464 nvar=size(volgrid6d_in%var)
1465 volgrid6d_out%var=volgrid6d_in%var
1466end if
1467! allocate once for speed
1468IF (.NOT.ASSOCIATED(volgrid6d_in%voldati)) THEN
1469 ALLOCATE(voldatiin(volgrid6d_in%griddim%dim%nx, volgrid6d_in%griddim%dim%ny, &
1470 inlevel))
1471ENDIF
1472IF (.NOT.ASSOCIATED(volgrid6d_out%voldati)) THEN
1473 ALLOCATE(voldatiout(volgrid6d_out%griddim%dim%nx, volgrid6d_out%griddim%dim%ny, &
1474 onlevel))
1475ENDIF
1476
1477CALL get_val(this, levshift=levshift, levused=levused)
1478spos = imiss
1479IF (c_e(lvar_coord_vol)) THEN
1480 CALL get_val(this%trans, output_levtype=output_levtype)
1481 IF (output_levtype%level1 == 103 .OR. output_levtype%level1 == 108) THEN
1482 spos = firsttrue(volgrid6d_in%level(:) == vol7d_level_new(1))
1483 IF (spos == 0) THEN
1484 CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
1485 'output level '//t2c(output_levtype%level1)// &
1486 ' requested, but height/press of surface not provided in volume')
1487 ENDIF
1488 IF (.NOT.c_e(levshift) .AND. .NOT.c_e(levused)) THEN
1489 CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
1490 'internal inconsistence, levshift and levused undefined when they should be')
1491 ENDIF
1492 ENDIF
1493ENDIF
1494
1495DO ivar=1,nvar
1496! IF (c_e(var_coord_vol)) THEN
1497! IF (ivar == var_coord_vol) CYCLE ! skip coordinate variable in output
1498! ENDIF
1499 DO itimerange=1,ntimerange
1500 DO itime=1,ntime
1501! skip empty columns where possible, improve
1502 IF (c_e(levshift) .AND. c_e(levused)) THEN
1503 IF (.NOT.any(c_e( &
1504 volgrid6d_in%gaid(levshift+1:levshift+levused,itime,itimerange,ivar) &
1505 ))) cycle
1506 ENDIF
1507 DO ilevel = 1, min(inlevel,onlevel)
1508! if present gaid copy it
1509 IF (c_e(volgrid6d_in%gaid(ilevel,itime,itimerange,ivar)) .AND. .NOT. &
1510 c_e(volgrid6d_out%gaid(ilevel,itime,itimerange,ivar))) THEN
1511
1512 IF (optio_log(clone)) THEN
1513 CALL copy(volgrid6d_in%gaid(ilevel,itime,itimerange,ivar),&
1514 volgrid6d_out%gaid(ilevel,itime,itimerange,ivar))
1515#ifdef DEBUG
1516 CALL l4f_category_log(volgrid6d_in%category,l4f_debug, &
1517 "cloning gaid, level "//t2c(ilevel))
1518#endif
1519 ELSE
1520 volgrid6d_out%gaid(ilevel,itime,itimerange,ivar) = &
1521 volgrid6d_in%gaid(ilevel,itime,itimerange,ivar)
1522 ENDIF
1523 ENDIF
1524 ENDDO
1525! if out level are more, we have to clone anyway
1526 DO ilevel = min(inlevel,onlevel) + 1, onlevel
1527 IF (c_e(volgrid6d_in%gaid(inlevel,itime,itimerange,ivar)) .AND. .NOT. &
1528 c_e(volgrid6d_out%gaid(ilevel,itime,itimerange,ivar))) then
1529
1530 CALL copy(volgrid6d_in%gaid(inlevel,itime,itimerange,ivar),&
1531 volgrid6d_out%gaid(ilevel,itime,itimerange,ivar))
1532#ifdef DEBUG
1533 CALL l4f_category_log(volgrid6d_in%category,l4f_debug, &
1534 "forced cloning gaid, level "//t2c(inlevel)//"->"//t2c(ilevel))
1535#endif
1536 ENDIF
1537 ENDDO
1538
1539 IF (c_e(lvar_coord_vol)) THEN
1540 NULLIFY(coord_3d_in)
1541 CALL volgrid_get_vol_3d(volgrid6d_in, itime, itimerange, lvar_coord_vol, &
1542 coord_3d_in)
1543 IF (c_e(spos)) THEN ! compute difference wrt surface coordinate
1544 IF (spos == 0) THEN ! error condition, set all to missing and goodnight
1545 coord_3d_in(:,:,levshift+1:levshift+levused) = rmiss
1546 ELSE
1547 DO ilevel = levshift+1, levshift+levused
1548 WHERE(c_e(coord_3d_in(:,:,ilevel)) .AND. c_e(coord_3d_in(:,:,spos)))
1549 coord_3d_in(:,:,ilevel) = coord_3d_in(:,:,ilevel) - &
1550 coord_3d_in(:,:,spos)
1551 ELSEWHERE
1552 coord_3d_in(:,:,ilevel) = rmiss
1553 END WHERE
1554 ENDDO
1555 ENDIF
1556 ENDIF
1557 ENDIF
1558 CALL volgrid_get_vol_3d(volgrid6d_in, itime, itimerange, ivar, &
1559 voldatiin)
1560 IF (ASSOCIATED(volgrid6d_out%voldati)) & ! improve!!!!
1561 CALL volgrid_get_vol_3d(volgrid6d_out, itime, itimerange, ivar, &
1562 voldatiout)
1563 IF (c_e(lvar_coord_vol)) THEN
1564 CALL compute(this, voldatiin, voldatiout, convert(volgrid6d_in%var(ivar)), &
1565 coord_3d_in(:,:,levshift+1:levshift+levused)) ! subset coord_3d_in
1566 ELSE
1567 CALL compute(this, voldatiin, voldatiout, convert(volgrid6d_in%var(ivar)))
1568 ENDIF
1569 CALL volgrid_set_vol_3d(volgrid6d_out, itime, itimerange, ivar, &
1570 voldatiout)
1571 ENDDO
1572 ENDDO
1573ENDDO
1574
1575IF (c_e(lvar_coord_vol)) THEN
1576 DEALLOCATE(coord_3d_in)
1577ENDIF
1578IF (.NOT.ASSOCIATED(volgrid6d_in%voldati)) THEN
1579 DEALLOCATE(voldatiin)
1580ENDIF
1581IF (.NOT.ASSOCIATED(volgrid6d_out%voldati)) THEN
1582 DEALLOCATE(voldatiout)
1583ENDIF
1585
1586END SUBROUTINE volgrid6d_transform_compute
1587
1588
1595SUBROUTINE volgrid6d_transform(this, griddim, volgrid6d_in, volgrid6d_out, &
1596 lev_out, volgrid6d_coord_in, maskgrid, maskbounds, clone, decode, categoryappend)
1597TYPE(transform_def),INTENT(in) :: this
1598TYPE(griddim_def),INTENT(in),OPTIONAL :: griddim
1599TYPE(volgrid6d),INTENT(inout) :: volgrid6d_in
1600TYPE(volgrid6d),INTENT(out) :: volgrid6d_out
1601TYPE(vol7d_level),INTENT(in),OPTIONAL,TARGET :: lev_out(:)
1602TYPE(volgrid6d),INTENT(in),OPTIONAL :: volgrid6d_coord_in
1603REAL,INTENT(in),OPTIONAL :: maskgrid(:,:)
1604REAL,INTENT(in),OPTIONAL :: maskbounds(:)
1605LOGICAL,INTENT(in),OPTIONAL :: clone
1606LOGICAL,INTENT(in),OPTIONAL :: decode
1607CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
1608
1609TYPE(grid_transform) :: grid_trans
1610TYPE(vol7d_level),POINTER :: llev_out(:)
1611TYPE(vol7d_level) :: input_levtype, output_levtype
1612TYPE(vol7d_var) :: vcoord_var
1613INTEGER :: i, k, ntime, ntimerange, nlevel, nvar, var_coord_in, var_coord_vol, &
1614 cf_out, nxc, nyc, nxi, nyi, i3, i4, i5, i6, &
1615 ulstart, ulend, spos
1616REAL,ALLOCATABLE :: coord_3d_in(:,:,:)
1617TYPE(geo_proj) :: proj_in, proj_out
1618CHARACTER(len=80) :: trans_type
1619LOGICAL :: ldecode
1620LOGICAL,ALLOCATABLE :: mask_in(:)
1621
1622#ifdef DEBUG
1623call l4f_category_log(volgrid6d_in%category, l4f_debug, "start volgrid6d_transform")
1624#endif
1625
1626ntime=0
1627ntimerange=0
1628nlevel=0
1629nvar=0
1630
1631if (associated(volgrid6d_in%time)) ntime=size(volgrid6d_in%time)
1632if (associated(volgrid6d_in%timerange)) ntimerange=size(volgrid6d_in%timerange)
1633if (associated(volgrid6d_in%level)) nlevel=size(volgrid6d_in%level)
1634if (associated(volgrid6d_in%var)) nvar=size(volgrid6d_in%var)
1635
1636IF (ntime == 0 .OR. ntimerange == 0 .OR. nlevel == 0 .OR. nvar == 0) THEN
1637 CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
1638 "trying to transform an incomplete volgrid6d object, ntime="//t2c(ntime)// &
1639 ' ntimerange='//t2c(ntimerange)//' nlevel='//t2c(nlevel)//' nvar='//t2c(nvar))
1640 CALL init(volgrid6d_out) ! initialize to empty
1641 CALL raise_error()
1642 RETURN
1643ENDIF
1644
1645CALL get_val(this, trans_type=trans_type)
1646
1647! store desired output component flag and unrotate if necessary
1648cf_out = imiss
1649IF (PRESENT(griddim) .AND. (trans_type == 'inter' .OR. trans_type == 'boxinter' &
1650 .OR. trans_type == 'stencilinter')) THEN ! improve condition!!
1651 CALL get_val(volgrid6d_in%griddim, proj=proj_in)
1652 CALL get_val(griddim, component_flag=cf_out, proj=proj_out)
1653! if different projections wind components must be referred to geographical system
1654 IF (proj_in /= proj_out) CALL vg6d_wind_unrot(volgrid6d_in)
1655ELSE IF (PRESENT(griddim)) THEN ! just get component_flag, the rest is rubbish
1656 CALL get_val(griddim, component_flag=cf_out)
1657ENDIF
1658
1659
1660var_coord_in = imiss
1661var_coord_vol = imiss
1662IF (trans_type == 'vertint') THEN
1663 IF (PRESENT(lev_out)) THEN
1664
1665! if volgrid6d_coord_in provided and allocated, check that it fits
1666 IF (PRESENT(volgrid6d_coord_in)) THEN
1667 IF (ASSOCIATED(volgrid6d_coord_in%voldati)) THEN
1668
1669! strictly 1 time and 1 timerange
1670 IF (SIZE(volgrid6d_coord_in%voldati,4) /= 1 .OR. &
1671 SIZE(volgrid6d_coord_in%voldati,5) /= 1) THEN
1672 CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
1673 'volume providing constant input vertical coordinate must have &
1674 &only 1 time and 1 timerange')
1675 CALL init(volgrid6d_out) ! initialize to empty
1676 CALL raise_error()
1677 RETURN
1678 ENDIF
1679
1680! search for variable providing vertical coordinate
1681 CALL get_val(this, output_levtype=output_levtype)
1682 vcoord_var = vol7d_var_new(vol7d_level_to_var(output_levtype))
1683 IF (.NOT.c_e(vcoord_var)) THEN
1684 CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
1685 'requested output level type '//t2c(output_levtype%level1)// &
1686 ' does not correspond to any known physical variable for &
1687 &providing vertical coordinate')
1688 CALL init(volgrid6d_out) ! initialize to empty
1689 CALL raise_error()
1690 RETURN
1691 ENDIF
1692
1693 DO i = 1, SIZE(volgrid6d_coord_in%var)
1694 IF (convert(volgrid6d_coord_in%var(i)) == vcoord_var) THEN
1695 var_coord_in = i
1696 EXIT
1697 ENDIF
1698 ENDDO
1699
1700 IF (.NOT.c_e(var_coord_in)) THEN
1701 CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
1702 'volume providing constant input vertical coordinate contains no &
1703 &variables matching output level type '//t2c(output_levtype%level1))
1704 CALL init(volgrid6d_out) ! initialize to empty
1705 CALL raise_error()
1706 RETURN
1707 ENDIF
1708 CALL l4f_category_log(volgrid6d_in%category, l4f_info, &
1709 'Coordinate for vertint found in coord volume at position '// &
1710 t2c(var_coord_in))
1711
1712! check horizontal grid
1713! this is too strict (component flag and so on)
1714! IF (volgrid6d_coord_in%griddim /= volgrid6d_in%griddim) THEN
1715 CALL get_val(volgrid6d_coord_in%griddim, nx=nxc, ny=nyc)
1716 CALL get_val(volgrid6d_in%griddim, nx=nxi, ny=nyi)
1717 IF (nxc /= nxi .OR. nyc /= nyi) THEN
1718 CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
1719 'volume providing constant input vertical coordinate must have &
1720 &the same grid as the input')
1721 CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
1722 'vertical coordinate: '//t2c(nxc)//'x'//t2c(nyc)// &
1723 ', input volume: '//t2c(nxi)//'x'//t2c(nyi))
1724 CALL init(volgrid6d_out) ! initialize to empty
1725 CALL raise_error()
1726 RETURN
1727 ENDIF
1728
1729! check vertical coordinate system
1730 CALL get_val(this, input_levtype=input_levtype)
1731 mask_in = & ! implicit allocation
1732 (volgrid6d_coord_in%level(:)%level1 == input_levtype%level1) .AND. &
1733 (volgrid6d_coord_in%level(:)%level2 == input_levtype%level2)
1734 ulstart = firsttrue(mask_in)
1735 ulend = lasttrue(mask_in)
1736 IF (ulstart == 0 .OR. ulend == 0) THEN
1737 CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
1738 'coordinate file does not contain levels of type '// &
1739 t2c(input_levtype%level1)//'/'//t2c(input_levtype%level2)// &
1740 ' specified for input data')
1741 CALL init(volgrid6d_out) ! initialize to empty
1742 CALL raise_error()
1743 RETURN
1744 ENDIF
1745
1746 coord_3d_in = volgrid6d_coord_in%voldati(:,:,ulstart:ulend,1,1,var_coord_in) ! implicit allocation
1747! special case
1748 IF (output_levtype%level1 == 103 .OR. &
1749 output_levtype%level1 == 108) THEN ! surface coordinate needed
1750 spos = firsttrue(volgrid6d_coord_in%level(:) == vol7d_level_new(1))
1751 IF (spos == 0) THEN
1752 CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
1753 'output level '//t2c(output_levtype%level1)// &
1754 ' requested, but height/press of surface not provided in coordinate file')
1755 CALL init(volgrid6d_out) ! initialize to empty
1756 CALL raise_error()
1757 RETURN
1758 ENDIF
1759 DO k = 1, SIZE(coord_3d_in,3)
1760 WHERE(c_e(coord_3d_in(:,:,k)) .AND. &
1761 c_e(volgrid6d_coord_in%voldati(:,:,spos,1,1,var_coord_in)))
1762 coord_3d_in(:,:,k) = coord_3d_in(:,:,k) - &
1763 volgrid6d_coord_in%voldati(:,:,spos,1,1,var_coord_in)
1764 ELSEWHERE
1765 coord_3d_in(:,:,k) = rmiss
1766 END WHERE
1767 ENDDO
1768 ENDIF
1769
1770 ENDIF
1771 ENDIF
1772
1773 IF (.NOT.c_e(var_coord_in)) THEN ! search for coordinate within volume
1774! search for variable providing vertical coordinate
1775 CALL get_val(this, output_levtype=output_levtype)
1776 vcoord_var = vol7d_var_new(vol7d_level_to_var(output_levtype))
1777 IF (c_e(vcoord_var)) THEN
1778 DO i = 1, SIZE(volgrid6d_in%var)
1779 IF (convert(volgrid6d_in%var(i)) == vcoord_var) THEN
1780 var_coord_vol = i
1781 EXIT
1782 ENDIF
1783 ENDDO
1784
1785 IF (c_e(var_coord_vol)) THEN
1786 CALL l4f_category_log(volgrid6d_in%category, l4f_info, &
1787 'Coordinate for vertint found in input volume at position '// &
1788 t2c(var_coord_vol))
1789 ENDIF
1790
1791 ENDIF
1792 ENDIF
1793
1794 CALL init(volgrid6d_out, griddim=volgrid6d_in%griddim, &
1795 time_definition=volgrid6d_in%time_definition, categoryappend=categoryappend)
1796 IF (c_e(var_coord_in)) THEN
1797 CALL init(grid_trans, this, lev_in=volgrid6d_in%level, lev_out=lev_out, &
1798 coord_3d_in=coord_3d_in, categoryappend=categoryappend)
1799 ELSE
1800 CALL init(grid_trans, this, lev_in=volgrid6d_in%level, lev_out=lev_out, &
1801 categoryappend=categoryappend)
1802 ENDIF
1803
1804 CALL get_val(grid_trans, output_level_auto=llev_out) ! get levels if auto-generated
1805 IF (.NOT.ASSOCIATED(llev_out)) llev_out => lev_out
1806 nlevel = SIZE(llev_out)
1807 ELSE
1808 CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
1809 'volgrid6d_transform: vertint requested but lev_out not provided')
1810 CALL init(volgrid6d_out) ! initialize to empty
1811 CALL raise_error()
1812 RETURN
1813 ENDIF
1814
1815ELSE
1816 CALL init(volgrid6d_out, griddim=griddim, &
1817 time_definition=volgrid6d_in%time_definition, categoryappend=categoryappend)
1818 CALL init(grid_trans, this, in=volgrid6d_in%griddim, out=volgrid6d_out%griddim, &
1819 maskgrid=maskgrid, maskbounds=maskbounds, categoryappend=categoryappend)
1820ENDIF
1821
1822
1823IF (c_e(grid_trans)) THEN ! transformation is valid
1824
1825 CALL volgrid6d_alloc(volgrid6d_out, ntime=ntime, nlevel=nlevel, &
1826 ntimerange=ntimerange, nvar=nvar)
1827
1828 IF (PRESENT(decode)) THEN ! explicitly set decode status
1829 ldecode = decode
1830 ELSE ! take it from input
1831 ldecode = ASSOCIATED(volgrid6d_in%voldati)
1832 ENDIF
1833! force decode if gaid is readonly
1834 decode_loop: DO i6 = 1,nvar
1835 DO i5 = 1, ntimerange
1836 DO i4 = 1, ntime
1837 DO i3 = 1, nlevel
1838 IF (c_e(volgrid6d_in%gaid(i3,i4,i5,i6))) THEN
1839 ldecode = ldecode .OR. grid_id_readonly(volgrid6d_in%gaid(i3,i4,i5,i6))
1840 EXIT decode_loop
1841 ENDIF
1842 ENDDO
1843 ENDDO
1844 ENDDO
1845 ENDDO decode_loop
1846
1847 IF (PRESENT(decode)) THEN
1848 IF (ldecode.NEQV.decode) THEN
1849 CALL l4f_category_log(volgrid6d_in%category, l4f_warn, &
1850 'volgrid6d_transform: decode status forced to .TRUE. because driver does not allow copy')
1851 ENDIF
1852 ENDIF
1853
1854 CALL volgrid6d_alloc_vol(volgrid6d_out, decode=ldecode)
1855
1856!ensure unproj was called
1857!call griddim_unproj(volgrid6d_out%griddim)
1858
1859 IF (trans_type == 'vertint') THEN
1860#ifdef DEBUG
1861 CALL l4f_category_log(volgrid6d_in%category, l4f_debug, &
1862 "volgrid6d_transform: vertint to "//t2c(nlevel)//" levels")
1863#endif
1864 CALL compute(grid_trans, volgrid6d_in, volgrid6d_out, lev_out=llev_out, &
1865 var_coord_vol=var_coord_vol, clone=clone)
1866 ELSE
1867 CALL compute(grid_trans, volgrid6d_in, volgrid6d_out, clone=clone)
1868 ENDIF
1869
1870 IF (cf_out == 0) THEN ! unrotated components are desired
1871 CALL wind_unrot(volgrid6d_out) ! unrotate if necessary
1872 ELSE IF (cf_out == 1) THEN ! rotated components are desired
1873 CALL wind_rot(volgrid6d_out) ! rotate if necessary
1874 ENDIF
1875
1876ELSE
1877! should log with grid_trans%category, but it is private
1878 CALL l4f_category_log(volgrid6d_in%category, l4f_error, &
1879 'volgrid6d_transform: transformation not valid')
1880 CALL raise_error()
1881ENDIF
1882
1883CALL delete (grid_trans)
1884
1885END SUBROUTINE volgrid6d_transform
1886
1887
1896SUBROUTINE volgrid6dv_transform(this, griddim, volgrid6d_in, volgrid6d_out, &
1897 lev_out, volgrid6d_coord_in, maskgrid, maskbounds, clone, decode, categoryappend)
1898TYPE(transform_def),INTENT(in) :: this
1899TYPE(griddim_def),INTENT(in),OPTIONAL :: griddim
1900TYPE(volgrid6d),INTENT(inout) :: volgrid6d_in(:)
1901TYPE(volgrid6d),POINTER :: volgrid6d_out(:)
1902TYPE(vol7d_level),INTENT(in),OPTIONAL :: lev_out(:)
1903TYPE(volgrid6d),INTENT(in),OPTIONAL :: volgrid6d_coord_in
1904REAL,INTENT(in),OPTIONAL :: maskgrid(:,:)
1905REAL,INTENT(in),OPTIONAL :: maskbounds(:)
1906LOGICAL,INTENT(in),OPTIONAL :: clone
1907LOGICAL,INTENT(in),OPTIONAL :: decode
1908CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
1909
1910INTEGER :: i, stallo
1911
1912
1913allocate(volgrid6d_out(size(volgrid6d_in)),stat=stallo)
1914if (stallo /= 0)then
1915 call l4f_log(l4f_fatal,"allocating memory")
1916 call raise_fatal_error()
1917end if
1918
1919do i=1,size(volgrid6d_in)
1920 call transform(this, griddim, volgrid6d_in(i), volgrid6d_out(i), &
1921 lev_out=lev_out, volgrid6d_coord_in=volgrid6d_coord_in, &
1922 maskgrid=maskgrid, maskbounds=maskbounds, &
1923 clone=clone, decode=decode, categoryappend=categoryappend)
1924end do
1925
1926END SUBROUTINE volgrid6dv_transform
1927
1928
1929! Internal method for performing grid to sparse point computations
1930SUBROUTINE volgrid6d_v7d_transform_compute(this, volgrid6d_in, vol7d_out, &
1931 networkname, noconvert)
1932TYPE(grid_transform),INTENT(in) :: this ! oggetto di trasformazione per grigliato
1933type(volgrid6d), INTENT(in) :: volgrid6d_in ! oggetto da trasformare
1934type(vol7d), INTENT(inout) :: vol7d_out ! oggetto trasformato
1935CHARACTER(len=*),OPTIONAL,INTENT(in) :: networkname ! imposta il network in vol7d_out (default='generic')
1936LOGICAL,OPTIONAL,INTENT(in) :: noconvert
1937
1938INTEGER :: nntime, nana, ntime, ntimerange, nlevel, nvar, stallo
1939INTEGER :: itime, itimerange, ivar, inetwork
1940REAL,ALLOCATABLE :: voldatir_out(:,:,:)
1941TYPE(conv_func),POINTER :: c_func(:)
1942TYPE(datetime),ALLOCATABLE :: validitytime(:,:)
1943REAL,POINTER :: voldatiin(:,:,:)
1944
1945#ifdef DEBUG
1946call l4f_category_log(volgrid6d_in%category,l4f_debug,"start volgrid6d_v7d_transform_compute")
1947#endif
1948
1949ntime=0
1950ntimerange=0
1951nlevel=0
1952nvar=0
1953NULLIFY(c_func)
1954
1955if (present(networkname))then
1956 call init(vol7d_out%network(1),name=networkname)
1957else
1958 call init(vol7d_out%network(1),name='generic')
1959end if
1960
1961if (associated(volgrid6d_in%timerange))then
1962 ntimerange=size(volgrid6d_in%timerange)
1963 vol7d_out%timerange=volgrid6d_in%timerange
1964end if
1965
1966if (associated(volgrid6d_in%time))then
1967 ntime=size(volgrid6d_in%time)
1968
1969 if (vol7d_out%time_definition == volgrid6d_in%time_definition) then
1970
1971 ! i time sono definiti uguali: assegno
1972 vol7d_out%time=volgrid6d_in%time
1973
1974 else
1975 ! converto reference in validity
1976 allocate (validitytime(ntime,ntimerange),stat=stallo)
1977 if (stallo /=0)then
1978 call l4f_category_log(volgrid6d_in%category,l4f_fatal,"allocating memory")
1979 call raise_fatal_error()
1980 end if
1981
1982 do itime=1,ntime
1983 do itimerange=1,ntimerange
1984 if (vol7d_out%time_definition > volgrid6d_in%time_definition) then
1985 validitytime(itime,itimerange) = &
1986 volgrid6d_in%time(itime) + timedelta_new(sec=volgrid6d_in%timerange(itimerange)%p1)
1987 else
1988 validitytime(itime,itimerange) = &
1989 volgrid6d_in%time(itime) - timedelta_new(sec=volgrid6d_in%timerange(itimerange)%p1)
1990 end if
1991 end do
1992 end do
1993
1994 nntime = count_distinct(reshape(validitytime,(/ntime*ntimerange/)), back=.true.)
1995 vol7d_out%time=pack_distinct(reshape(validitytime,(/ntime*ntimerange/)), nntime,back=.true.)
1996
1997 end if
1998end if
1999
2000IF (ASSOCIATED(volgrid6d_in%level)) THEN
2001 nlevel = SIZE(volgrid6d_in%level)
2002 vol7d_out%level=volgrid6d_in%level
2003ENDIF
2004
2005IF (ASSOCIATED(volgrid6d_in%var)) THEN
2006 nvar = SIZE(volgrid6d_in%var)
2007 IF (.NOT. optio_log(noconvert)) THEN
2008 CALL vargrib2varbufr(volgrid6d_in%var, vol7d_out%dativar%r, c_func)
2009 ENDIF
2010ENDIF
2011
2012nana = SIZE(vol7d_out%ana)
2013
2014! allocate once for speed
2015IF (.NOT.ASSOCIATED(volgrid6d_in%voldati)) THEN
2016 ALLOCATE(voldatiin(volgrid6d_in%griddim%dim%nx, volgrid6d_in%griddim%dim%ny, &
2017 nlevel))
2018ENDIF
2019
2020ALLOCATE(voldatir_out(nana,1,nlevel),stat=stallo)
2021IF (stallo /= 0) THEN
2022 CALL l4f_category_log(volgrid6d_in%category,l4f_fatal,"allocating memory")
2023 CALL raise_fatal_error()
2024ENDIF
2025
2026inetwork=1
2027do itime=1,ntime
2028 do itimerange=1,ntimerange
2029! do ilevel=1,nlevel
2030 do ivar=1,nvar
2031
2032 !non è chiaro se questa sezione è utile o no
2033 !ossia il compute sotto sembra prevedere voldatir_out solo in out
2034!!$ if (vol7d_out%time_definition == volgrid6d_in%time_definition) then
2035!!$ voldatir_out=reshape(vol7d_out%voldatir(:,itime,ilevel,itimerange,ivar,inetwork),(/nana,1/))
2036!!$ else
2037!!$ voldatir_out=reshape(vol7d_out%voldatir(:,index(vol7d_out%time,validitytime(itime,ilevel)),ilevel,itimerange,ivar,inetwork),(/nana,1/))
2038!!$ end if
2039
2040 CALL volgrid_get_vol_3d(volgrid6d_in, itime, itimerange, ivar, &
2041 voldatiin)
2042
2043 CALL compute(this, voldatiin, voldatir_out, vol7d_out%dativar%r(ivar))
2044
2045 if (vol7d_out%time_definition == volgrid6d_in%time_definition) then
2046 vol7d_out%voldatir(:,itime,:,itimerange,ivar,inetwork) = &
2047 voldatir_out(:,1,:)
2048 else
2049 vol7d_out%voldatir(:,index(vol7d_out%time,validitytime(itime,itimerange)),:,itimerange,ivar,inetwork)=&
2050 reshape(voldatir_out,(/nana,nlevel/))
2051 end if
2052
2053! 1 indice della dimensione "anagrafica"
2054! 2 indice della dimensione "tempo"
2055! 3 indice della dimensione "livello verticale"
2056! 4 indice della dimensione "intervallo temporale"
2057! 5 indice della dimensione "variabile"
2058! 6 indice della dimensione "rete"
2059
2060 end do
2061! end do
2062 end do
2063end do
2064
2065deallocate(voldatir_out)
2066IF (.NOT.ASSOCIATED(volgrid6d_in%voldati)) THEN
2067 DEALLOCATE(voldatiin)
2068ENDIF
2069if (allocated(validitytime)) deallocate(validitytime)
2070
2071! Rescale valid data according to variable conversion table
2072IF (ASSOCIATED(c_func)) THEN
2073 DO ivar = 1, nvar
2074 CALL compute(c_func(ivar), vol7d_out%voldatir(:,:,:,:,ivar,:))
2075 ENDDO
2076 DEALLOCATE(c_func)
2077ENDIF
2078
2079end SUBROUTINE volgrid6d_v7d_transform_compute
2080
2081
2088SUBROUTINE volgrid6d_v7d_transform(this, volgrid6d_in, vol7d_out, v7d, &
2089 maskgrid, maskbounds, networkname, noconvert, find_index, categoryappend)
2090TYPE(transform_def),INTENT(in) :: this
2091TYPE(volgrid6d),INTENT(inout) :: volgrid6d_in
2092TYPE(vol7d),INTENT(out) :: vol7d_out
2093TYPE(vol7d),INTENT(in),OPTIONAL :: v7d
2094REAL,INTENT(in),OPTIONAL :: maskgrid(:,:)
2095REAL,INTENT(in),OPTIONAL :: maskbounds(:)
2096CHARACTER(len=*),OPTIONAL,INTENT(in) :: networkname
2097LOGICAL,OPTIONAL,INTENT(in) :: noconvert
2098PROCEDURE(basic_find_index),POINTER,OPTIONAL :: find_index
2099CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
2100
2101type(grid_transform) :: grid_trans
2102INTEGER :: ntime, ntimerange, nlevel, nvar, nana, time_definition, nnetwork, stallo
2103INTEGER :: itime, itimerange, inetwork
2104TYPE(datetime),ALLOCATABLE :: validitytime(:,:)
2105INTEGER,ALLOCATABLE :: point_index(:)
2106TYPE(vol7d) :: v7d_locana
2107
2108#ifdef DEBUG
2109call l4f_category_log(volgrid6d_in%category,l4f_debug,"start volgrid6d_v7d_transform")
2110#endif
2111
2112call vg6d_wind_unrot(volgrid6d_in)
2113
2114ntime=0
2115ntimerange=0
2116nlevel=0
2117nvar=0
2118nnetwork=1
2119
2120call get_val(this,time_definition=time_definition)
2121if (.not. c_e(time_definition)) then
2122 time_definition=1 ! default to validity time
2123endif
2124
2125IF (PRESENT(v7d)) THEN
2126 CALL vol7d_copy(v7d, v7d_locana)
2127ELSE
2128 CALL init(v7d_locana)
2129ENDIF
2130
2131if (associated(volgrid6d_in%timerange)) ntimerange=size(volgrid6d_in%timerange)
2132
2133if (associated(volgrid6d_in%time)) then
2134
2135 ntime=size(volgrid6d_in%time)
2136
2137 if (time_definition /= volgrid6d_in%time_definition) then
2138
2139 ! converto reference in validity
2140 allocate (validitytime(ntime,ntimerange),stat=stallo)
2141 if (stallo /=0)then
2142 call l4f_category_log(volgrid6d_in%category,l4f_fatal,"allocating memory")
2143 call raise_fatal_error()
2144 end if
2145
2146 do itime=1,ntime
2147 do itimerange=1,ntimerange
2148 if (time_definition > volgrid6d_in%time_definition) then
2149 validitytime(itime,itimerange) = &
2150 volgrid6d_in%time(itime) + timedelta_new(sec=volgrid6d_in%timerange(itimerange)%p1)
2151 else
2152 validitytime(itime,itimerange) = &
2153 volgrid6d_in%time(itime) - timedelta_new(sec=volgrid6d_in%timerange(itimerange)%p1)
2154 end if
2155 end do
2156 end do
2157
2158 ntime = count_distinct(reshape(validitytime,(/ntime*ntimerange/)), back=.true.)
2159 deallocate (validitytime)
2160
2161 end if
2162end if
2163
2164
2165if (associated(volgrid6d_in%level)) nlevel=size(volgrid6d_in%level)
2166if (associated(volgrid6d_in%var)) nvar=size(volgrid6d_in%var)
2167
2168CALL init(grid_trans, this, volgrid6d_in%griddim, v7d_locana, &
2169 maskgrid=maskgrid, maskbounds=maskbounds, find_index=find_index, &
2170 categoryappend=categoryappend)
2171CALL init (vol7d_out,time_definition=time_definition)
2172
2173IF (c_e(grid_trans)) THEN
2174
2175 nana=SIZE(v7d_locana%ana)
2176 CALL vol7d_alloc(vol7d_out, nana=nana, ntime=ntime, nlevel=nlevel, &
2177 ntimerange=ntimerange, ndativarr=nvar, nnetwork=nnetwork)
2178 vol7d_out%ana = v7d_locana%ana
2179
2180 CALL get_val(grid_trans, output_point_index=point_index)
2181 IF (ALLOCATED(point_index)) THEN
2182! check that size(point_index) == nana?
2183 CALL vol7d_alloc(vol7d_out, nanavari=1)
2184 CALL init(vol7d_out%anavar%i(1), 'B01192')
2185 ENDIF
2186
2187 CALL vol7d_alloc_vol(vol7d_out)
2188
2189 IF (ALLOCATED(point_index)) THEN
2190 DO inetwork = 1, nnetwork
2191 vol7d_out%volanai(:,1,inetwork) = point_index(:)
2192 ENDDO
2193 ENDIF
2194 CALL compute(grid_trans, volgrid6d_in, vol7d_out, networkname, noconvert)
2195ELSE
2196 CALL l4f_log(l4f_error, 'vg6d_v7d_transform: transformation not valid')
2197 CALL raise_error()
2198ENDIF
2199
2200CALL delete(grid_trans)
2201
2202#ifdef HAVE_DBALLE
2203! set variables to a conformal state
2204CALL vol7d_dballe_set_var_du(vol7d_out)
2205#endif
2206
2207CALL delete(v7d_locana)
2208
2209END SUBROUTINE volgrid6d_v7d_transform
2210
2211
2220SUBROUTINE volgrid6dv_v7d_transform(this, volgrid6d_in, vol7d_out, v7d, &
2221 maskgrid, maskbounds, networkname, noconvert, find_index, categoryappend)
2222TYPE(transform_def),INTENT(in) :: this
2223TYPE(volgrid6d),INTENT(inout) :: volgrid6d_in(:)
2224TYPE(vol7d),INTENT(out) :: vol7d_out
2225TYPE(vol7d),INTENT(in),OPTIONAL :: v7d
2226REAL,INTENT(in),OPTIONAL :: maskgrid(:,:)
2227REAL,INTENT(in),OPTIONAL :: maskbounds(:)
2228CHARACTER(len=*),OPTIONAL,INTENT(in) :: networkname
2229LOGICAL,OPTIONAL,INTENT(in) :: noconvert
2230PROCEDURE(basic_find_index),POINTER,OPTIONAL :: find_index
2231CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
2232
2233integer :: i
2234TYPE(vol7d) :: v7dtmp
2235
2236
2237CALL init(v7dtmp)
2238CALL init(vol7d_out)
2239
2240DO i=1,SIZE(volgrid6d_in)
2241 CALL transform(this, volgrid6d_in(i), v7dtmp, v7d=v7d, &
2242 maskgrid=maskgrid, maskbounds=maskbounds, &
2243 networkname=networkname, noconvert=noconvert, find_index=find_index, &
2244 categoryappend=categoryappend)
2245 CALL vol7d_append(vol7d_out, v7dtmp)
2246ENDDO
2247
2248END SUBROUTINE volgrid6dv_v7d_transform
2249
2250
2251! Internal method for performing sparse point to grid computations
2252SUBROUTINE v7d_volgrid6d_transform_compute(this, vol7d_in, volgrid6d_out, networkname, gaid_template)
2253TYPE(grid_transform),INTENT(in) :: this ! object specifying the specific transformation
2254type(vol7d), INTENT(in) :: vol7d_in ! object to be transformed
2255type(volgrid6d), INTENT(inout) :: volgrid6d_out ! transformed object
2256CHARACTER(len=*),OPTIONAL,INTENT(in) :: networkname ! select the network to be processed from the \a vol7d input object, default the first network
2257TYPE(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
2258
2259integer :: nana, ntime, ntimerange, nlevel, nvar
2260INTEGER :: ilevel, itime, itimerange, ivar, inetwork
2261
2262REAL,POINTER :: voldatiout(:,:,:)
2263type(vol7d_network) :: network
2264TYPE(conv_func), pointer :: c_func(:)
2265!TODO category sarebbe da prendere da vol7d
2266#ifdef DEBUG
2267CALL l4f_category_log(volgrid6d_out%category, l4f_debug, &
2268 'start v7d_volgrid6d_transform_compute')
2269#endif
2270
2271ntime=0
2272ntimerange=0
2273nlevel=0
2274nvar=0
2275
2276IF (PRESENT(networkname)) THEN
2277 CALL init(network,name=networkname)
2278 inetwork = index(vol7d_in%network,network)
2279 IF (inetwork <= 0) THEN
2280 CALL l4f_category_log(volgrid6d_out%category,l4f_warn, &
2281 'network '//trim(networkname)//' not found, first network will be transformed')
2282 inetwork = 1
2283 ENDIF
2284ELSE
2285 inetwork = 1
2286ENDIF
2287
2288! no time_definition conversion implemented, output will be the same as input
2289if (associated(vol7d_in%time))then
2290 ntime=size(vol7d_in%time)
2291 volgrid6d_out%time=vol7d_in%time
2292end if
2293
2294if (associated(vol7d_in%timerange))then
2295 ntimerange=size(vol7d_in%timerange)
2296 volgrid6d_out%timerange=vol7d_in%timerange
2297end if
2298
2299if (associated(vol7d_in%level))then
2300 nlevel=size(vol7d_in%level)
2301 volgrid6d_out%level=vol7d_in%level
2302end if
2303
2304if (associated(vol7d_in%dativar%r))then
2305 nvar=size(vol7d_in%dativar%r)
2306 CALL varbufr2vargrib(vol7d_in%dativar%r, volgrid6d_out%var, c_func, gaid_template)
2307end if
2308
2309nana=SIZE(vol7d_in%voldatir, 1)
2310! allocate once for speed
2311IF (.NOT.ASSOCIATED(volgrid6d_out%voldati)) THEN
2312 ALLOCATE(voldatiout(volgrid6d_out%griddim%dim%nx, volgrid6d_out%griddim%dim%ny, &
2313 nlevel))
2314ENDIF
2315
2316DO ivar=1,nvar
2317 DO itimerange=1,ntimerange
2318 DO itime=1,ntime
2319
2320! clone the gaid template where I have data
2321 IF (PRESENT(gaid_template)) THEN
2322 DO ilevel = 1, nlevel
2323 IF (any(c_e(vol7d_in%voldatir(:,itime,ilevel,itimerange,ivar,inetwork)))) THEN
2324 CALL copy(gaid_template, volgrid6d_out%gaid(ilevel,itime,itimerange,ivar))
2325 ELSE
2326 volgrid6d_out%gaid(ilevel,itime,itimerange,ivar) = grid_id_new()
2327 ENDIF
2328 ENDDO
2329 ENDIF
2330
2331! get data
2332 IF (ASSOCIATED(volgrid6d_out%voldati)) & ! improve!!!!
2333 CALL volgrid_get_vol_3d(volgrid6d_out, itime, itimerange, ivar, &
2334 voldatiout)
2335! do the interpolation
2336 CALL compute(this, &
2337 vol7d_in%voldatir(:,itime,:,itimerange,ivar,inetwork), voldatiout, &
2338 vol7d_in%dativar%r(ivar))
2339! rescale valid data according to variable conversion table
2340 IF (ASSOCIATED(c_func)) THEN
2341 CALL compute(c_func(ivar), voldatiout(:,:,:))
2342 ENDIF
2343! put data
2344 CALL volgrid_set_vol_3d(volgrid6d_out, itime, itimerange, ivar, &
2345 voldatiout)
2346
2347 ENDDO
2348 ENDDO
2349ENDDO
2350
2351IF (.NOT.ASSOCIATED(volgrid6d_out%voldati)) THEN
2352 DEALLOCATE(voldatiout)
2353ENDIF
2354IF (ASSOCIATED(c_func)) THEN
2355 DEALLOCATE(c_func)
2356ENDIF
2357
2358END SUBROUTINE v7d_volgrid6d_transform_compute
2359
2360
2367SUBROUTINE v7d_volgrid6d_transform(this, griddim, vol7d_in, volgrid6d_out, &
2368 networkname, gaid_template, categoryappend)
2369TYPE(transform_def),INTENT(in) :: this
2370TYPE(griddim_def),INTENT(in),OPTIONAL :: griddim
2371! TODO ripristinare intent(in) dopo le opportune modifiche in grid_class.F90
2372TYPE(vol7d),INTENT(inout) :: vol7d_in
2373TYPE(volgrid6d),INTENT(out) :: volgrid6d_out
2374CHARACTER(len=*),OPTIONAL,INTENT(in) :: networkname
2375TYPE(grid_id),OPTIONAL,INTENT(in) :: gaid_template
2376CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
2377
2378type(grid_transform) :: grid_trans
2379integer :: ntime, ntimerange, nlevel, nvar
2380
2381
2382!TODO la category sarebbe da prendere da vol7d
2383!call l4f_category_log(vol7d_out%category,L4F_DEBUG,"start volgrid6d_transform")
2384
2385CALL vol7d_alloc_vol(vol7d_in) ! be safe
2386ntime=SIZE(vol7d_in%time)
2387ntimerange=SIZE(vol7d_in%timerange)
2388nlevel=SIZE(vol7d_in%level)
2389nvar=0
2390if (associated(vol7d_in%dativar%r)) nvar=size(vol7d_in%dativar%r)
2391
2392IF (nvar <= 0) THEN ! use vol7d category once it will be implemented
2393 CALL l4f_log(l4f_error, &
2394 "trying to transform a vol7d object incomplete or without real variables")
2395 CALL init(volgrid6d_out) ! initialize to empty
2396 CALL raise_error()
2397 RETURN
2398ENDIF
2399
2400CALL init(grid_trans, this, vol7d_in, griddim, categoryappend=categoryappend)
2401CALL init(volgrid6d_out, griddim, time_definition=vol7d_in%time_definition, &
2402 categoryappend=categoryappend)
2403
2404IF (c_e(grid_trans)) THEN
2405
2406 CALL volgrid6d_alloc(volgrid6d_out, griddim%dim, ntime=ntime, nlevel=nlevel, &
2407 ntimerange=ntimerange, nvar=nvar)
2408! can I avoid decode=.TRUE. here, using gaid_template?
2409 CALL volgrid6d_alloc_vol(volgrid6d_out, decode=.true.)
2410
2411 CALL compute(grid_trans, vol7d_in, volgrid6d_out, networkname, gaid_template)
2412
2413 CALL vg6d_wind_rot(volgrid6d_out)
2414ELSE
2415! should log with grid_trans%category, but it is private
2416 CALL l4f_log(l4f_error, 'v7d_vg6d_transform: transformation not valid')
2417 CALL raise_error()
2418ENDIF
2419
2420CALL delete(grid_trans)
2421
2422END SUBROUTINE v7d_volgrid6d_transform
2423
2424
2425! Internal method for performing sparse point to sparse point computations
2426SUBROUTINE v7d_v7d_transform_compute(this, vol7d_in, vol7d_out, lev_out, &
2427 var_coord_vol)
2428TYPE(grid_transform),INTENT(in) :: this ! oggetto di trasformazione per grigliato
2429type(vol7d), INTENT(in) :: vol7d_in ! oggetto da trasformare
2430type(vol7d), INTENT(inout) :: vol7d_out ! oggetto trasformato
2431TYPE(vol7d_level),INTENT(in),OPTIONAL :: lev_out(:) ! vol7d_level object defining target vertical grid, for vertical interpolations
2432INTEGER,INTENT(in),OPTIONAL :: var_coord_vol ! index of variable defining vertical coordinate values in input volume
2433
2434INTEGER :: itime, itimerange, ilevel, ivar, inetwork, &
2435 levshift, levused, lvar_coord_vol, spos
2436REAL,ALLOCATABLE :: coord_3d_in(:,:,:)
2437TYPE(vol7d_level) :: output_levtype
2438
2439lvar_coord_vol = optio_i(var_coord_vol)
2440vol7d_out%time(:) = vol7d_in%time(:)
2441vol7d_out%timerange(:) = vol7d_in%timerange(:)
2442IF (PRESENT(lev_out)) THEN
2443 vol7d_out%level(:) = lev_out(:)
2444ELSE
2445 vol7d_out%level(:) = vol7d_in%level(:)
2446ENDIF
2447vol7d_out%network(:) = vol7d_in%network(:)
2448IF (ASSOCIATED(vol7d_in%dativar%r)) THEN ! work only when real vars are available
2449 vol7d_out%dativar%r(:) = vol7d_in%dativar%r(:)
2450
2451 CALL get_val(this, levshift=levshift, levused=levused)
2452 spos = imiss
2453 IF (c_e(lvar_coord_vol)) THEN
2454 CALL get_val(this%trans, output_levtype=output_levtype)
2455 IF (output_levtype%level1 == 103 .OR. output_levtype%level1 == 108) THEN
2456 spos = firsttrue(vol7d_in%level(:) == vol7d_level_new(1))
2457 IF (spos == 0) THEN
2458 CALL l4f_log(l4f_error, &
2459 'output level '//t2c(output_levtype%level1)// &
2460 ' requested, but height/press of surface not provided in volume')
2461 ENDIF
2462 IF (.NOT.c_e(levshift) .AND. .NOT.c_e(levused)) THEN
2463 CALL l4f_log(l4f_error, &
2464 'internal inconsistence, levshift and levused undefined when they should be')
2465 ENDIF
2466 ALLOCATE(coord_3d_in(SIZE(vol7d_in%ana),1,SIZE(vol7d_in%level)))
2467 ENDIF
2468
2469 ENDIF
2470
2471 DO inetwork = 1, SIZE(vol7d_in%network)
2472 DO ivar = 1, SIZE(vol7d_in%dativar%r)
2473 DO itimerange = 1, SIZE(vol7d_in%timerange)
2474 DO itime = 1, SIZE(vol7d_in%time)
2475
2476! dirty trick to make voldatir look like a 2d-array of shape (nana,1)
2477 IF (c_e(lvar_coord_vol)) THEN
2478 IF (c_e(spos)) THEN ! compute difference wrt surface coordinate
2479 IF (spos == 0) THEN ! error condition, set all to missing and goodnight
2480 coord_3d_in(:,:,levshift+1:levshift+levused) = rmiss
2481 ELSE
2482 DO ilevel = levshift+1, levshift+levused
2483 WHERE(c_e(vol7d_in%voldatir(:,itime:itime,ilevel,itimerange,lvar_coord_vol,inetwork)) .AND. &
2484 c_e(vol7d_in%voldatir(:,itime:itime,spos,itimerange,lvar_coord_vol,inetwork)))
2485 coord_3d_in(:,:,ilevel) = vol7d_in%voldatir(:,itime:itime,ilevel,itimerange,lvar_coord_vol,inetwork) - &
2486 vol7d_in%voldatir(:,itime:itime,spos,itimerange,lvar_coord_vol,inetwork)
2487 ELSEWHERE
2488 coord_3d_in(:,:,ilevel) = rmiss
2489 END WHERE
2490 ENDDO
2491 ENDIF
2492 CALL compute(this, &
2493 vol7d_in%voldatir(:,itime,:,itimerange,ivar,inetwork), &
2494 vol7d_out%voldatir(:,itime:itime,:,itimerange,ivar,inetwork), &
2495 var=vol7d_in%dativar%r(ivar), &
2496 coord_3d_in=coord_3d_in)
2497 ELSE
2498 CALL compute(this, &
2499 vol7d_in%voldatir(:,itime,:,itimerange,ivar,inetwork), &
2500 vol7d_out%voldatir(:,itime:itime,:,itimerange,ivar,inetwork), &
2501 var=vol7d_in%dativar%r(ivar), &
2502 coord_3d_in=vol7d_in%voldatir(:,itime:itime,:,itimerange, &
2503 lvar_coord_vol,inetwork))
2504 ENDIF
2505 ELSE
2506 CALL compute(this, &
2507 vol7d_in%voldatir(:,itime,:,itimerange,ivar,inetwork), &
2508 vol7d_out%voldatir(:,itime:itime,:,itimerange,ivar,inetwork), &
2509 var=vol7d_in%dativar%r(ivar))
2510 ENDIF
2511 ENDDO
2512 ENDDO
2513 ENDDO
2514 ENDDO
2515
2516ENDIF
2517
2518END SUBROUTINE v7d_v7d_transform_compute
2519
2520
2528SUBROUTINE v7d_v7d_transform(this, vol7d_in, vol7d_out, v7d, maskbounds, &
2529 lev_out, vol7d_coord_in, categoryappend)
2530TYPE(transform_def),INTENT(in) :: this
2531TYPE(vol7d),INTENT(inout) :: vol7d_in
2532TYPE(vol7d),INTENT(out) :: vol7d_out
2533TYPE(vol7d),INTENT(in),OPTIONAL :: v7d
2534REAL,INTENT(in),OPTIONAL :: maskbounds(:)
2535TYPE(vol7d_level),INTENT(in),OPTIONAL,TARGET :: lev_out(:)
2536TYPE(vol7d),INTENT(in),OPTIONAL :: vol7d_coord_in
2537CHARACTER(len=*),INTENT(in),OPTIONAL :: categoryappend
2538
2539INTEGER :: nvar, inetwork
2540TYPE(grid_transform) :: grid_trans
2541TYPE(vol7d_level),POINTER :: llev_out(:)
2542TYPE(vol7d_level) :: input_levtype, output_levtype
2543TYPE(vol7d_var) :: vcoord_var
2544REAL,ALLOCATABLE :: coord_3d_in(:,:,:)
2545INTEGER :: var_coord_in, var_coord_vol, i, k, ulstart, ulend, spos
2546INTEGER,ALLOCATABLE :: point_index(:)
2547TYPE(vol7d) :: v7d_locana, vol7d_tmpana
2548CHARACTER(len=80) :: trans_type
2549LOGICAL,ALLOCATABLE :: mask_in(:), point_mask(:)
2550
2551CALL vol7d_alloc_vol(vol7d_in) ! be safe
2552nvar=0
2553IF (ASSOCIATED(vol7d_in%dativar%r)) nvar=SIZE(vol7d_in%dativar%r)
2554
2555CALL init(v7d_locana)
2556IF (PRESENT(v7d)) v7d_locana = v7d
2557CALL init(vol7d_out, time_definition=vol7d_in%time_definition)
2558
2559CALL get_val(this, trans_type=trans_type)
2560
2561var_coord_vol = imiss
2562IF (trans_type == 'vertint') THEN
2563
2564 IF (PRESENT(lev_out)) THEN
2565
2566! if vol7d_coord_in provided and allocated, check that it fits
2567 var_coord_in = -1
2568 IF (PRESENT(vol7d_coord_in)) THEN
2569 IF (ASSOCIATED(vol7d_coord_in%voldatir) .AND. &
2570 ASSOCIATED(vol7d_coord_in%dativar%r)) THEN
2571
2572! strictly 1 time, 1 timerange and 1 network
2573 IF (SIZE(vol7d_coord_in%voldatir,2) /= 1 .OR. &
2574 SIZE(vol7d_coord_in%voldatir,4) /= 1 .OR. &
2575 SIZE(vol7d_coord_in%voldatir,6) /= 1) THEN
2576 CALL l4f_log(l4f_error, &
2577 'volume providing constant input vertical coordinate must have &
2578 &only 1 time, 1 timerange and 1 network')
2579 CALL raise_error()
2580 RETURN
2581 ENDIF
2582
2583! search for variable providing vertical coordinate
2584 CALL get_val(this, output_levtype=output_levtype)
2585 vcoord_var = vol7d_var_new(vol7d_level_to_var(output_levtype))
2586 IF (.NOT.c_e(vcoord_var)) THEN
2587 CALL l4f_log(l4f_error, &
2588 'requested output level type '//t2c(output_levtype%level1)// &
2589 ' does not correspond to any known physical variable for &
2590 &providing vertical coordinate')
2591 CALL raise_error()
2592 RETURN
2593 ENDIF
2594
2595 var_coord_in = index(vol7d_coord_in%dativar%r, vcoord_var)
2596
2597 IF (var_coord_in <= 0) THEN
2598 CALL l4f_log(l4f_error, &
2599 'volume providing constant input vertical coordinate contains no &
2600 &real variables matching output level type '//t2c(output_levtype%level1))
2601 CALL raise_error()
2602 RETURN
2603 ENDIF
2604 CALL l4f_log(l4f_info, &
2605 'Coordinate for vertint found in coord volume at position '// &
2606 t2c(var_coord_in))
2607
2608! check vertical coordinate system
2609 CALL get_val(this, input_levtype=input_levtype)
2610 mask_in = & ! implicit allocation
2611 (vol7d_coord_in%level(:)%level1 == input_levtype%level1) .AND. &
2612 (vol7d_coord_in%level(:)%level2 == input_levtype%level2)
2613 ulstart = firsttrue(mask_in)
2614 ulend = lasttrue(mask_in)
2615 IF (ulstart == 0 .OR. ulend == 0) THEN
2616 CALL l4f_log(l4f_error, &
2617 'coordinate file does not contain levels of type '// &
2618 t2c(input_levtype%level1)//'/'//t2c(input_levtype%level2)// &
2619 ' specified for input data')
2620 CALL raise_error()
2621 RETURN
2622 ENDIF
2623
2624 coord_3d_in = vol7d_coord_in%voldatir(:,1:1,ulstart:ulend,1,var_coord_in,1) ! dirty 1:1, implicit allocation
2625! special case
2626 IF (output_levtype%level1 == 103 &
2627 .OR. output_levtype%level1 == 108) THEN ! surface coordinate needed
2628 spos = firsttrue(vol7d_coord_in%level(:) == vol7d_level_new(1))
2629 IF (spos == 0) THEN
2630 CALL l4f_log(l4f_error, &
2631 'output level '//t2c(output_levtype%level1)// &
2632 ' requested, but height/press of surface not provided in coordinate file')
2633 CALL raise_error()
2634 RETURN
2635 ENDIF
2636 DO k = 1, SIZE(coord_3d_in,3)
2637 WHERE(c_e(coord_3d_in(:,:,k)) .AND. &
2638 c_e(vol7d_coord_in%voldatir(:,1:1,spos,1,var_coord_in,1)))
2639 coord_3d_in(:,:,k) = coord_3d_in(:,:,k) - &
2640 vol7d_coord_in%voldatir(:,1:1,spos,1,var_coord_in,1)
2641 ELSEWHERE
2642 coord_3d_in(:,:,k) = rmiss
2643 END WHERE
2644 ENDDO
2645 ENDIF
2646
2647 ENDIF
2648 ENDIF
2649
2650 IF (var_coord_in <= 0) THEN ! search for coordinate within volume
2651! search for variable providing vertical coordinate
2652 CALL get_val(this, output_levtype=output_levtype)
2653 vcoord_var = vol7d_var_new(vol7d_level_to_var(output_levtype))
2654 IF (c_e(vcoord_var)) THEN
2655 DO i = 1, SIZE(vol7d_in%dativar%r)
2656 IF (vol7d_in%dativar%r(i) == vcoord_var) THEN
2657 var_coord_vol = i
2658 EXIT
2659 ENDIF
2660 ENDDO
2661
2662 IF (c_e(var_coord_vol)) THEN
2663 CALL l4f_log(l4f_info, &
2664 'Coordinate for vertint found in input volume at position '// &
2665 t2c(var_coord_vol))
2666 ENDIF
2667
2668 ENDIF
2669 ENDIF
2670
2671 IF (var_coord_in > 0) THEN
2672 CALL init(grid_trans, this, lev_in=vol7d_in%level, lev_out=lev_out, &
2673 coord_3d_in=coord_3d_in, categoryappend=categoryappend)
2674 ELSE
2675 CALL init(grid_trans, this, lev_in=vol7d_in%level, lev_out=lev_out, &
2676 categoryappend=categoryappend)
2677 ENDIF
2678
2679 CALL get_val(grid_trans, output_level_auto=llev_out) ! get levels if auto-generated
2680 IF (.NOT.associated(llev_out)) llev_out => lev_out
2681
2682 IF (c_e(grid_trans)) then! .AND. nvar > 0) THEN
2683
2684 CALL vol7d_alloc(vol7d_out, nana=SIZE(vol7d_in%ana), &
2685 ntime=SIZE(vol7d_in%time), ntimerange=SIZE(vol7d_in%timerange), &
2686 nlevel=SIZE(llev_out), nnetwork=SIZE(vol7d_in%network), ndativarr=nvar)
2687 vol7d_out%ana(:) = vol7d_in%ana(:)
2688
2689 CALL vol7d_alloc_vol(vol7d_out)
2690
2691! no need to check c_e(var_coord_vol) here since the presence of
2692! this%coord_3d_in (external) has precedence over coord_3d_in internal
2693! in grid_transform_compute
2694 CALL compute(grid_trans, vol7d_in, vol7d_out, llev_out, &
2695 var_coord_vol=var_coord_vol)
2696 ELSE
2697 CALL l4f_log(l4f_error, 'v7d_v7d_transform: transformation not valid')
2698 CALL raise_error()
2699 ENDIF
2700 ELSE
2701 CALL l4f_log(l4f_error, &
2702 'v7d_v7d_transform: vertint requested but lev_out not provided')
2703 CALL raise_error()
2704 ENDIF
2705
2706ELSE
2707
2708 CALL init(grid_trans, this, vol7d_in, v7d_locana, maskbounds=maskbounds, &
2709 categoryappend=categoryappend)
2710! if this init is successful, I am sure that v7d_locana%ana is associated
2711
2712 IF (c_e(grid_trans)) then! .AND. nvar > 0) THEN
2713
2714 CALL vol7d_alloc(vol7d_out, nana=SIZE(v7d_locana%ana), &
2715 ntime=SIZE(vol7d_in%time), ntimerange=SIZE(vol7d_in%timerange), &
2716 nlevel=SIZE(vol7d_in%level), nnetwork=SIZE(vol7d_in%network), ndativarr=nvar)
2717 vol7d_out%ana = v7d_locana%ana
2718
2719 CALL get_val(grid_trans, point_mask=point_mask, output_point_index=point_index)
2720
2721 IF (ALLOCATED(point_index)) THEN
2722 CALL vol7d_alloc(vol7d_out, nanavari=1)
2723 CALL init(vol7d_out%anavar%i(1), 'B01192')
2724 ENDIF
2725
2726 CALL vol7d_alloc_vol(vol7d_out)
2727
2728 IF (ALLOCATED(point_index)) THEN
2729 DO inetwork = 1, SIZE(vol7d_in%network)
2730 vol7d_out%volanai(:,1,inetwork) = point_index(:)
2731 ENDDO
2732 ENDIF
2733 CALL compute(grid_trans, vol7d_in, vol7d_out)
2734
2735 IF (ALLOCATED(point_mask)) THEN ! keep full ana
2736 IF (SIZE(point_mask) /= SIZE(vol7d_in%ana)) THEN
2737 CALL l4f_log(l4f_warn, &
2738 'v7d_v7d_transform: inconsistency in point size: '//t2c(SIZE(point_mask)) &
2739 //':'//t2c(SIZE(vol7d_in%ana)))
2740 ELSE
2741#ifdef DEBUG
2742 CALL l4f_log(l4f_debug, 'v7d_v7d_transform: merging ana from in to out')
2743#endif
2744 CALL vol7d_copy(vol7d_in, vol7d_tmpana, &
2745 lana=point_mask, lnetwork=(/.true./), &
2746 ltime=(/.false./), ltimerange=(/.false./), llevel=(/.false./))
2747 CALL vol7d_append(vol7d_out, vol7d_tmpana)
2748 ENDIF
2749 ENDIF
2750
2751 ELSE
2752 CALL l4f_log(l4f_error, 'v7d_v7d_transform: transformation not valid')
2753 CALL raise_error()
2754 ENDIF
2755
2756ENDIF
2757
2758CALL delete (grid_trans)
2759IF (.NOT. PRESENT(v7d)) CALL delete(v7d_locana)
2760
2761END SUBROUTINE v7d_v7d_transform
2762
2763
2771subroutine vg6d_wind_unrot(this)
2772type(volgrid6d) :: this
2773
2774integer :: component_flag
2775
2776call get_val(this%griddim,component_flag=component_flag)
2777
2778if (component_flag == 1) then
2779 call l4f_category_log(this%category,l4f_info, &
2780 "unrotating vector components")
2781 call vg6d_wind__un_rot(this,.false.)
2782 call set_val(this%griddim,component_flag=0)
2783else
2784 call l4f_category_log(this%category,l4f_info, &
2785 "no need to unrotate vector components")
2786end if
2787
2788end subroutine vg6d_wind_unrot
2789
2790
2796subroutine vg6d_wind_rot(this)
2797type(volgrid6d) :: this
2798
2799integer :: component_flag
2800
2801call get_val(this%griddim,component_flag=component_flag)
2802
2803if (component_flag == 0) then
2804 call l4f_category_log(this%category,l4f_info, &
2805 "rotating vector components")
2806 call vg6d_wind__un_rot(this,.true.)
2807 call set_val(this%griddim,component_flag=1)
2808else
2809 call l4f_category_log(this%category,l4f_info, &
2810 "no need to rotate vector components")
2811end if
2812
2813end subroutine vg6d_wind_rot
2814
2815
2816! Generic UnRotate the wind components.
2817SUBROUTINE vg6d_wind__un_rot(this,rot)
2818TYPE(volgrid6d) :: this ! object containing wind to be (un)rotated
2819LOGICAL :: rot ! if .true. rotate else unrotate
2820
2821INTEGER :: i, j, k, l, a11, a12, a21, a22, stallo
2822double precision,pointer :: rot_mat(:,:,:)
2823real,allocatable :: tmp_arr(:,:)
2824REAL,POINTER :: voldatiu(:,:), voldativ(:,:)
2825INTEGER,POINTER :: iu(:), iv(:)
2826
2827IF (.NOT.ASSOCIATED(this%var)) THEN
2828 CALL l4f_category_log(this%category, l4f_error, &
2829 "trying to unrotate an incomplete volgrid6d object")
2830 CALL raise_fatal_error()
2831! RETURN
2832ENDIF
2833
2834CALL volgrid6d_var_hor_comp_index(this%var, iu, iv)
2835IF (.NOT.ASSOCIATED(iu)) THEN
2836 CALL l4f_category_log(this%category,l4f_error, &
2837 "unrotation impossible")
2838 CALL raise_fatal_error()
2839! RETURN
2840ENDIF
2841
2842! Temporary workspace
2843ALLOCATE(tmp_arr(this%griddim%dim%nx, this%griddim%dim%ny),stat=stallo)
2844IF (stallo /= 0) THEN
2845 CALL l4f_category_log(this%category, l4f_fatal, "allocating memory")
2846 CALL raise_fatal_error()
2847ENDIF
2848! allocate once for speed
2849IF (.NOT.ASSOCIATED(this%voldati)) THEN
2850 ALLOCATE(voldatiu(this%griddim%dim%nx, this%griddim%dim%ny), &
2851 voldativ(this%griddim%dim%nx, this%griddim%dim%ny))
2852ENDIF
2853
2854CALL griddim_unproj(this%griddim)
2855CALL wind_unrot(this%griddim, rot_mat)
2856
2857a11=1
2858if (rot)then
2859 a12=2
2860 a21=3
2861else
2862 a12=3
2863 a21=2
2864end if
2865a22=4
2866
2867DO l = 1, SIZE(iu)
2868 DO k = 1, SIZE(this%timerange)
2869 DO j = 1, SIZE(this%time)
2870 DO i = 1, SIZE(this%level)
2871! get data
2872 CALL volgrid_get_vol_2d(this, i, j, k, iu(l), voldatiu)
2873 CALL volgrid_get_vol_2d(this, i, j, k, iv(l), voldativ)
2874! convert units forward
2875! CALL compute(conv_fwd(iu(l)), voldatiu)
2876! CALL compute(conv_fwd(iv(l)), voldativ)
2877
2878! multiply wind components by rotation matrix
2879 WHERE(voldatiu /= rmiss .AND. voldativ /= rmiss)
2880 tmp_arr(:,:) = real(voldatiu(:,:)*rot_mat(:,:,a11) + &
2881 voldativ(:,:)*rot_mat(:,:,a12))
2882 voldativ(:,:) = real(voldatiu(:,:)*rot_mat(:,:,a21) + &
2883 voldativ(:,:)*rot_mat(:,:,a22))
2884 voldatiu(:,:) = tmp_arr(:,:)
2885 END WHERE
2886! convert units backward
2887! CALL uncompute(conv_fwd(iu(l)), voldatiu)
2888! CALL uncompute(conv_fwd(iv(l)), voldativ)
2889! put data
2890 CALL volgrid_set_vol_2d(this, i, j, k, iu(l), voldatiu)
2891 CALL volgrid_set_vol_2d(this, i, j, k, iv(l), voldativ)
2892 ENDDO
2893 ENDDO
2894 ENDDO
2895ENDDO
2896
2897IF (.NOT.ASSOCIATED(this%voldati)) THEN
2898 DEALLOCATE(voldatiu, voldativ)
2899ENDIF
2900DEALLOCATE(rot_mat, tmp_arr, iu, iv)
2901
2902END SUBROUTINE vg6d_wind__un_rot
2903
2904
2905!!$ try to understand the problem:
2906!!$
2907!!$ case:
2908!!$
2909!!$ 1) we have only one volume: we have to provide the direction of shift
2910!!$ compute H and traslate on it
2911!!$ 2) we have two volumes:
2912!!$ 1) volume U and volume V: compute H and traslate on it
2913!!$ 2) volume U/V and volume H : translate U/V on H
2914!!$ 3) we have tree volumes: translate U and V on H
2915!!$
2916!!$ strange cases:
2917!!$ 1) do not have U in volume U
2918!!$ 2) do not have V in volume V
2919!!$ 3) we have others variables more than U and V in volumes U e V
2920!!$
2921!!$
2922!!$ so the steps are:
2923!!$ 1) find the volumes
2924!!$ 2) define or compute H grid
2925!!$ 3) trasform the volumes in H
2926
2927!!$ N.B.
2928!!$ case 1) for only one vg6d (U or V) is not managed, but
2929!!$ the not pubblic subroutines will work but you have to know what you want to do
2930
2931
2948subroutine vg6d_c2a (this)
2949
2950TYPE(volgrid6d),INTENT(inout) :: this(:)
2951
2952integer :: ngrid,igrid,jgrid,ugrid,vgrid,tgrid
2953doubleprecision :: xmin, xmax, ymin, ymax
2954doubleprecision :: xmin_t, xmax_t, ymin_t, ymax_t
2955doubleprecision :: step_lon_t,step_lat_t
2956character(len=80) :: type_t,type
2957TYPE(griddim_def):: griddim_t
2958
2959ngrid=size(this)
2960
2961do igrid=1,ngrid
2962
2963 call init(griddim_t)
2964
2965 call get_val(this(igrid)%griddim,xmin=xmin_t, xmax=xmax_t, ymin=ymin_t, ymax=ymax_t,proj_type=type_t)
2966 step_lon_t=(xmax_t-xmin_t)/dble(this(igrid)%griddim%dim%nx-1)
2967 step_lat_t=(ymax_t-ymin_t)/dble(this(igrid)%griddim%dim%ny-1)
2968
2969 do jgrid=1,ngrid
2970
2971 ugrid=imiss
2972 vgrid=imiss
2973 tgrid=imiss
2974
2975#ifdef DEBUG
2976 call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: search U/V/T points:"//to_char(igrid)//to_char(jgrid))
2977 call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: test delta: "//to_char(step_lon_t)//to_char(step_lat_t))
2978#endif
2979
2980 if (this(igrid)%griddim == this(jgrid)%griddim ) cycle
2981
2982 if (this(igrid)%griddim%dim%nx == this(jgrid)%griddim%dim%nx .and. &
2983 this(igrid)%griddim%dim%ny == this(jgrid)%griddim%dim%ny ) then
2985 call get_val(this(jgrid)%griddim,xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,proj_type=type)
2986
2987 if (type_t /= type )cycle
2988
2989#ifdef DEBUG
2990 call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: test U "//&
2991 to_char(xmin)//to_char(xmax)//to_char(ymin)//to_char(ymax))
2992
2993 call l4f_category_log(this(igrid)%category,l4f_debug,"diff coordinate lon"//&
2994 to_char(abs(xmin - (xmin_t+step_lon_t/2.d0)))//&
2995 to_char(abs(xmax - (xmax_t+step_lon_t/2.d0))))
2996 call l4f_category_log(this(igrid)%category,l4f_debug,"diff coordinate lat"//&
2997 to_char(abs(ymin - (ymin_t+step_lat_t/2.d0)))//&
2998 to_char(abs(ymax - (ymax_t+step_lat_t/2.d0))))
2999#endif
3000
3001 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
3002 if ( abs(ymin - ymin_t) < 1.d-3 .and. abs(ymax - ymax_t) < 1.d-3 ) then
3003
3004#ifdef DEBUG
3005 call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: found U")
3006#endif
3007 ugrid=jgrid
3008 tgrid=igrid
3009
3010 end if
3011 end if
3012
3013#ifdef DEBUG
3014 call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: test V "//&
3015 to_char(xmin)//to_char(xmax)//to_char(ymin)//to_char(ymax))
3016#endif
3017
3018 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
3019 if ( abs(xmin - xmin_t) < 1.d-3 .and. abs(xmax - xmax_t) < 1.d-3 ) then
3020
3021#ifdef DEBUG
3022 call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: found V")
3023#endif
3024 vgrid=jgrid
3025 tgrid=igrid
3026
3027 end if
3028 end if
3029
3030
3031 ! test if we have U and V
3032
3033#ifdef DEBUG
3034 call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: test U and V"//&
3035 to_char(xmin_t)//to_char(xmax_t)//to_char(ymin_t)//to_char(ymax_t)//&
3036 to_char(xmin)//to_char(xmax)//to_char(ymin)//to_char(ymax))
3037
3038 call l4f_category_log(this(igrid)%category,l4f_debug,"UV diff coordinate lon"//&
3039 to_char(abs(xmin_t - xmin)-step_lon_t/2.d0)//&
3040 to_char(abs(xmax_t - xmax)-step_lon_t/2.d0))
3041 call l4f_category_log(this(igrid)%category,l4f_debug,"UV diff coordinate lat"//&
3042 to_char(abs(ymin_t - ymin) -step_lat_t/2.d0)//&
3043 to_char(abs(ymax_t - ymax)-step_lat_t/2.d0))
3044#endif
3045
3046 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
3047 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
3048
3049#ifdef DEBUG
3050 call l4f_category_log(this(igrid)%category,l4f_debug,"C grid: found U and V case up and right")
3051#endif
3052
3053 vgrid=jgrid
3054 ugrid=igrid
3055
3056 call init(griddim_t,xmin=xmin, xmax=xmax, ymin=ymin_t, ymax=ymax_t)
3057
3058 end if
3059 end if
3060 end if
3061
3062 ! abbiamo almeno due volumi: riportiamo U e/o V su T
3063 if (c_e(ugrid)) then
3064#ifdef DEBUG
3065 call l4f_category_log(this(igrid)%category,l4f_debug,"C grid U points: interpolate U on T "//to_char(tgrid)//to_char(ugrid))
3066#endif
3067 if (c_e(tgrid))then
3068 call vg6d_c2a_grid(this(ugrid),this(tgrid)%griddim,cgrid=1)
3069 else
3070 call vg6d_c2a_grid(this(ugrid),griddim_t,cgrid=1)
3071 end if
3072 call vg6d_c2a_mat(this(ugrid),cgrid=1)
3073 end if
3074
3075 if (c_e(vgrid)) then
3076#ifdef DEBUG
3077 call l4f_category_log(this(igrid)%category,l4f_debug,"C grid V points: interpolate V on T "//to_char(tgrid)//to_char(vgrid))
3078#endif
3079 if (c_e(tgrid))then
3080 call vg6d_c2a_grid(this(vgrid),this(tgrid)%griddim,cgrid=2)
3081 else
3082 call vg6d_c2a_grid(this(vgrid),griddim_t,cgrid=2)
3083 end if
3084 call vg6d_c2a_mat(this(vgrid),cgrid=2)
3085 end if
3086
3087 end do
3088
3089 call delete(griddim_t)
3090
3091end do
3092
3093
3094end subroutine vg6d_c2a
3095
3096
3097! Convert C grid to A grid
3098subroutine vg6d_c2a_grid(this,griddim_t,cgrid)
3099
3100type(volgrid6d),intent(inout) :: this ! object containing fields to be translated (U or V points)
3101type(griddim_def),intent(in),optional :: griddim_t ! object containing grid of T points
3102integer,intent(in) :: cgrid ! in C grid (Arakawa) we have 0=T,1=U,2=V points
3103
3104doubleprecision :: xmin, xmax, ymin, ymax
3105doubleprecision :: step_lon,step_lat
3106
3107
3108if (present(griddim_t)) then
3109
3110 call get_val(griddim_t,xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
3111 call set_val(this%griddim,xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
3112! improve grid definition precision
3113 CALL griddim_setsteps(this%griddim)
3114
3115else
3116
3117 select case (cgrid)
3118
3119 case(0)
3120
3121#ifdef DEBUG
3122 call l4f_category_log(this%category,l4f_debug,"C grid: T points, nothing to do")
3123#endif
3124 return
3125
3126 case (1)
3127
3128#ifdef DEBUG
3129 call l4f_category_log(this%category,l4f_debug,"C grid: U points, we need interpolation")
3130#endif
3131
3132 call get_val(this%griddim, xmin=xmin, xmax=xmax)
3133 step_lon=(xmax-xmin)/dble(this%griddim%dim%nx-1)
3134 xmin=xmin-step_lon/2.d0
3135 xmax=xmax-step_lon/2.d0
3136 call set_val(this%griddim, xmin=xmin, xmax=xmax)
3137! improve grid definition precision
3138 CALL griddim_setsteps(this%griddim)
3139
3140 case (2)
3141
3142#ifdef DEBUG
3143 call l4f_category_log(this%category,l4f_debug,"C grid: V points, we need interpolation")
3144#endif
3145
3146 call get_val(this%griddim, ymin=ymin, ymax=ymax)
3147 step_lat=(ymax-ymin)/dble(this%griddim%dim%ny-1)
3148 ymin=ymin-step_lat/2.d0
3149 ymax=ymax-step_lat/2.d0
3150 call set_val(this%griddim, ymin=ymin, ymax=ymax)
3151! improve grid definition precision
3152 CALL griddim_setsteps(this%griddim)
3153
3154 case default
3155
3156 call l4f_category_log(this%category,l4f_fatal,"C grid type not known")
3157 call raise_fatal_error ()
3158
3159 end select
3160
3161end if
3162
3163
3164call griddim_unproj(this%griddim)
3165
3166
3167end subroutine vg6d_c2a_grid
3168
3169! Convert C grid to A grid
3170subroutine vg6d_c2a_mat(this,cgrid)
3171
3172type(volgrid6d),intent(inout) :: this ! object containing fields to be translated
3173integer,intent(in) :: cgrid ! in C grid (Arakawa) we have 0=T,1=U,2=V points
3174
3175INTEGER :: i, j, k, iv, stallo
3176REAL,ALLOCATABLE :: tmp_arr(:,:)
3177REAL,POINTER :: voldatiuv(:,:)
3178
3179
3180IF (cgrid == 0) RETURN ! nothing to do
3181IF (cgrid /= 1 .AND. cgrid /= 2) THEN ! wrong cgrid
3182 CALL l4f_category_log(this%category,l4f_fatal,"C grid type "// &
3183 trim(to_char(cgrid))//" not known")
3184 CALL raise_error()
3185 RETURN
3186ENDIF
3187
3188! Temporary workspace
3189ALLOCATE(tmp_arr(this%griddim%dim%nx, this%griddim%dim%ny),stat=stallo)
3190if (stallo /=0)then
3191 call l4f_log(l4f_fatal,"allocating memory")
3192 call raise_fatal_error()
3193end if
3194
3195! allocate once for speed
3196IF (.NOT.ASSOCIATED(this%voldati)) THEN
3197 ALLOCATE(voldatiuv(this%griddim%dim%nx, this%griddim%dim%ny), stat=stallo)
3198 IF (stallo /= 0) THEN
3199 CALL l4f_log(l4f_fatal,"allocating memory")
3200 CALL raise_fatal_error()
3201 ENDIF
3202ENDIF
3203
3204IF (cgrid == 1) THEN ! U points to H points
3205 DO iv = 1, SIZE(this%var)
3206 DO k = 1, SIZE(this%timerange)
3207 DO j = 1, SIZE(this%time)
3208 DO i = 1, SIZE(this%level)
3209 tmp_arr(:,:) = rmiss
3210 CALL volgrid_get_vol_2d(this, i, j, k, iv, voldatiuv)
3211
3212! West boundary extrapolation
3213 WHERE(voldatiuv(1,:) /= rmiss .AND. voldatiuv(2,:) /= rmiss)
3214 tmp_arr(1,:) = voldatiuv(1,:) - (voldatiuv(2,:) - voldatiuv(1,:)) / 2.
3215 END WHERE
3216
3217! Rest of the field
3218 WHERE(voldatiuv(1:this%griddim%dim%nx-1,:) /= rmiss .AND. &
3219 voldatiuv(2:this%griddim%dim%nx,:) /= rmiss)
3220 tmp_arr(2:this%griddim%dim%nx,:) = &
3221 (voldatiuv(1:this%griddim%dim%nx-1,:) + &
3222 voldatiuv(2:this%griddim%dim%nx,:)) / 2.
3223 END WHERE
3224
3225 voldatiuv(:,:) = tmp_arr
3226 CALL volgrid_set_vol_2d(this, i, j, k, iv, voldatiuv)
3227 ENDDO
3228 ENDDO
3229 ENDDO
3230 ENDDO
3231
3232ELSE IF (cgrid == 2) THEN ! V points to H points
3233 DO iv = 1, SIZE(this%var)
3234 DO k = 1, SIZE(this%timerange)
3235 DO j = 1, SIZE(this%time)
3236 DO i = 1, SIZE(this%level)
3237 tmp_arr(:,:) = rmiss
3238 CALL volgrid_get_vol_2d(this, i, j, k, iv, voldatiuv)
3239
3240! South boundary extrapolation
3241 WHERE(voldatiuv(:,1) /= rmiss .AND. voldatiuv(:,2) /= rmiss)
3242 tmp_arr(:,1) = voldatiuv(:,1) - (voldatiuv(:,2) - voldatiuv(:,1)) / 2.
3243 END WHERE
3244
3245! Rest of the field
3246 WHERE(voldatiuv(:,1:this%griddim%dim%ny-1) /= rmiss .AND. &
3247 voldatiuv(:,2:this%griddim%dim%ny) /= rmiss)
3248 tmp_arr(:,2:this%griddim%dim%ny) = &
3249 (voldatiuv(:,1:this%griddim%dim%ny-1) + &
3250 voldatiuv(:,2:this%griddim%dim%ny)) / 2.
3251 END WHERE
3252
3253 voldatiuv(:,:) = tmp_arr
3254 CALL volgrid_set_vol_2d(this, i, j, k, iv, voldatiuv)
3255 ENDDO
3256 ENDDO
3257 ENDDO
3258 ENDDO
3259ENDIF ! tertium non datur
3260
3261IF (.NOT.ASSOCIATED(this%voldati)) THEN
3262 DEALLOCATE(voldatiuv)
3263ENDIF
3264DEALLOCATE (tmp_arr)
3265
3266end subroutine vg6d_c2a_mat
3267
3268
3272subroutine display_volgrid6d (this)
3273
3274TYPE(volgrid6d),intent(in) :: this
3275integer :: i
3276
3277#ifdef DEBUG
3278call l4f_category_log(this%category,l4f_debug,"ora mostro gridinfo " )
3279#endif
3280
3281print*,"----------------------- volgrid6d display ---------------------"
3282call display(this%griddim)
3283
3284IF (ASSOCIATED(this%time))then
3285 print*,"---- time vector ----"
3286 print*,"elements=",size(this%time)
3287 do i=1, size(this%time)
3288 call display(this%time(i))
3289 end do
3290end IF
3291
3292IF (ASSOCIATED(this%timerange))then
3293 print*,"---- timerange vector ----"
3294 print*,"elements=",size(this%timerange)
3295 do i=1, size(this%timerange)
3296 call display(this%timerange(i))
3297 end do
3298end IF
3299
3300IF (ASSOCIATED(this%level))then
3301 print*,"---- level vector ----"
3302 print*,"elements=",size(this%level)
3303 do i=1, size(this%level)
3304 call display(this%level(i))
3305 end do
3306end IF
3307
3308IF (ASSOCIATED(this%var))then
3309 print*,"---- var vector ----"
3310 print*,"elements=",size(this%var)
3311 do i=1, size(this%var)
3312 call display(this%var(i))
3313 end do
3314end IF
3315
3316IF (ASSOCIATED(this%gaid))then
3317 print*,"---- gaid vector (present mask only) ----"
3318 print*,"elements=",shape(this%gaid)
3319 print*,c_e(reshape(this%gaid,(/SIZE(this%gaid)/)))
3320end IF
3321
3322print*,"--------------------------------------------------------------"
3323
3324
3325end subroutine display_volgrid6d
3326
3327
3331subroutine display_volgrid6dv (this)
3332
3333TYPE(volgrid6d),intent(in) :: this(:)
3334integer :: i
3335
3336print*,"----------------------- volgrid6d vector ---------------------"
3337
3338print*,"elements=",size(this)
3339
3340do i=1, size(this)
3341
3342#ifdef DEBUG
3343 call l4f_category_log(this(i)%category,l4f_debug,"ora mostro il vettore volgrid6d" )
3344#endif
3345
3346 call display(this(i))
3347
3348end do
3349print*,"--------------------------------------------------------------"
3350
3351end subroutine display_volgrid6dv
3352
3353
3356subroutine vg6dv_rounding(vg6din,vg6dout,level,timerange,nostatproc,merge)
3357type(volgrid6d),intent(in) :: vg6din(:)
3358type(volgrid6d),intent(out),pointer :: vg6dout(:)
3359type(vol7d_level),intent(in),optional :: level(:)
3360type(vol7d_timerange),intent(in),optional :: timerange(:)
3361logical,intent(in),optional :: merge
3362logical,intent(in),optional :: nostatproc
3363
3364integer :: i
3365
3366allocate(vg6dout(size(vg6din)))
3367
3368do i = 1, size(vg6din)
3369 call vg6d_rounding(vg6din(i),vg6dout(i),level,timerange,nostatproc,merge)
3370end do
3371
3372end subroutine vg6dv_rounding
3373
3385subroutine vg6d_rounding(vg6din,vg6dout,level,timerange,nostatproc,merge)
3386type(volgrid6d),intent(in) :: vg6din
3387type(volgrid6d),intent(out) :: vg6dout
3388type(vol7d_level),intent(in),optional :: level(:)
3389type(vol7d_timerange),intent(in),optional :: timerange(:)
3390logical,intent(in),optional :: merge
3392logical,intent(in),optional :: nostatproc
3393
3394integer :: ilevel,itimerange
3395type(vol7d_level) :: roundlevel(size(vg6din%level))
3396type(vol7d_timerange) :: roundtimerange(size(vg6din%timerange))
3397
3398roundlevel=vg6din%level
3399
3400if (present(level))then
3401 do ilevel = 1, size(vg6din%level)
3402 if ((any(vg6din%level(ilevel) .almosteq. level))) then
3403 roundlevel(ilevel)=level(1)
3404 end if
3405 end do
3406end if
3407
3408roundtimerange=vg6din%timerange
3409
3410if (present(timerange))then
3411 do itimerange = 1, size(vg6din%timerange)
3412 if ((any(vg6din%timerange(itimerange) .almosteq. timerange))) then
3413 roundtimerange(itimerange)=timerange(1)
3414 end if
3415 end do
3416end if
3417
3418!set istantaneous values everywere
3419!preserve p1 for forecast time
3420if (optio_log(nostatproc)) then
3421 roundtimerange(:)%timerange=254
3422 roundtimerange(:)%p2=0
3423end if
3424
3425
3426call vg6d_reduce(vg6din,vg6dout,roundlevel,roundtimerange,merge)
3427
3428end subroutine vg6d_rounding
3429
3438subroutine vg6d_reduce(vg6din,vg6dout,roundlevel,roundtimerange,merge)
3439type(volgrid6d),intent(in) :: vg6din
3440type(volgrid6d),intent(out) :: vg6dout
3441type(vol7d_level),intent(in) :: roundlevel(:)
3442type(vol7d_timerange),intent(in) :: roundtimerange(:)
3443logical,intent(in),optional :: merge
3444
3445integer :: nlevel,ntime,ntimerange,nvar,ilevel,itimerange,ivar,indl,indt,itime,nx,ny
3446real,allocatable :: vol2d(:,:)
3447
3448nx=vg6din%griddim%dim%nx
3449ny=vg6din%griddim%dim%ny
3450nlevel=count_distinct(roundlevel,back=.true.)
3451ntime=size(vg6din%time)
3452ntimerange=count_distinct(roundtimerange,back=.true.)
3453nvar=size(vg6din%var)
3454
3455call init(vg6dout, vg6din%griddim, vg6din%time_definition, categoryappend="generated by vg6d_reduce")
3456call volgrid6d_alloc(vg6dout, vg6din%griddim%dim, ntime, nlevel, ntimerange, nvar)
3457
3458if ( ASSOCIATED(vg6din%voldati) .or. optio_log(merge)) then
3459 call volgrid6d_alloc_vol(vg6dout,inivol=.true.,decode=.true.)
3460 allocate(vol2d(nx,ny))
3461else
3462 call volgrid6d_alloc_vol(vg6dout,inivol=.true.,decode=.false.)
3463end if
3464
3465vg6dout%time=vg6din%time
3466vg6dout%var=vg6din%var
3467vg6dout%timerange=pack_distinct(roundtimerange,ntimerange,back=.true.)
3468vg6dout%level=pack_distinct(roundlevel,nlevel,back=.true.)
3469! sort modified dimensions
3470CALL sort(vg6dout%timerange)
3471CALL sort(vg6dout%level)
3472
3473do ilevel=1,size(vg6din%level)
3474 indl=index(vg6dout%level,roundlevel(ilevel))
3475 do itimerange=1,size(vg6din%timerange)
3476 indt=index(vg6dout%timerange,roundtimerange(itimerange))
3477 do ivar=1, nvar
3478 do itime=1,ntime
3479
3480 if ( ASSOCIATED(vg6din%voldati)) then
3481 vol2d=vg6din%voldati(:,:,ilevel,itime,itimerange,ivar)
3482 end if
3483
3484 if (optio_log(merge)) then
3485
3486 if ( .not. ASSOCIATED(vg6din%voldati)) then
3487 CALL grid_id_decode_data(vg6din%gaid(ilevel,itime,itimerange,ivar), vol2d)
3488 end if
3489
3490 !! merge present data point by point
3491 where (.not. c_e(vg6dout%voldati(:,:,indl,itime,indt,ivar)))
3492
3493 vg6dout%voldati(:,:,indl,itime,indt,ivar)=vol2d
3494
3495 end where
3496 else if ( ASSOCIATED(vg6din%voldati)) then
3497 if (.not. any(c_e(vg6dout%voldati(:,:,indl,itime,indt,ivar))))then
3498 vg6dout%voldati(:,:,indl,itime,indt,ivar)=vol2d
3499 end if
3500 end if
3501
3502 if (c_e(vg6din%gaid(ilevel,itime,itimerange,ivar)).and. .not. c_e(vg6dout%gaid(indl,itime,indt,ivar)))then
3503 call copy (vg6din%gaid(ilevel,itime,itimerange,ivar), vg6dout%gaid(indl,itime,indt,ivar))
3504 end if
3505 end do
3506 end do
3507 end do
3508end do
3509
3510if ( ASSOCIATED(vg6din%voldati) .or. optio_log(merge)) then
3511 deallocate(vol2d)
3512end if
3513
3514end subroutine vg6d_reduce
3515
3516
3517end module volgrid6d_class
3518
3519
3520
3525
3532
3535
3538
3541
Method for inserting elements of the array at a desired position.
Operatore di valore assoluto di un intervallo.
Legge un oggetto datetime/timedelta o un vettore di oggetti datetime/timedelta da un file FORMATTED o...
Functions that return a trimmed CHARACTER representation of the input variable.
Restituiscono il valore dell'oggetto in forma di stringa stampabile.
Scrive un oggetto datetime/timedelta o un vettore di oggetti datetime/timedelta su un file FORMATTED ...
Copy an object, creating a fully new instance.
Method for returning the contents of the object.
Compute forward coordinate transformation from geographical system to projected system.
Method for setting the contents of the object.
Clone the object, creating a new independent instance of the object exactly equal to the starting one...
Decode and return the data array from a grid_id object associated to a gridinfo object.
Encode a data array into a grid_id object associated to a gridinfo object.
Index method.
Emit log message for a category with specific priority.
Convert a level type to a physical variable.
Destructor, it releases every information and memory buffer associated with the object.
Display on standard output a description of the volgrid6d object provided.
Export an object dirctly to a native file, to a gridinfo object or to a supported file format through...
Import an object dirctly from a native file, from a gridinfo object or from a supported file format t...
Constructor, it creates a new instance of the object.
Reduce some dimensions (level and timerage) for semplification (rounding).
Transform between any combination of volgrid6d and vol7d objects by means of a transform_def object d...
Apply the conversion function this to values.
This module defines usefull general purpose function and subroutine.
Classi per la gestione delle coordinate temporali.
Module for describing geographically referenced regular grids.
Module for defining the extension and coordinates of a rectangular georeferenced grid.
This module defines an abstract interface to different drivers for access to files containing gridded...
Module for defining transformations between rectangular georeferenced grids and between grids and spa...
Class for managing information about a single gridded georeferenced field, typically imported from an...
classe per la gestione del logging
Module for quickly interpreting the OPTIONAL parameters passed to a subprogram.
Classe per la gestione di un volume completo di dati osservati.
classe per import ed export di volumi da e in DB-All.e
Classe per la gestione dei livelli verticali in osservazioni meteo e affini.
Classe per la gestione degli intervalli temporali di osservazioni meteo e affini.
This module defines objects and methods for managing data volumes on rectangular georeferenced grids.
Class for managing physical variables in a grib 1/2 fashion.
Class for expressing an absolute time value.
This object completely describes a grid on a geographic projection.
Derived type describing the extension of a grid and the geographical coordinates of each point.
Derived type associated to a block/message/record/band of gridded data coming from a file-like object...
This object fully defines a transformation between a couple of particular griddim_def or vol7d object...
This object defines in an abstract way the type of transformation to be applied.
Derived type defining a dynamically extensible array of TYPE(gridinfo_def) elements.
Object describing a single gridded message/band.
Definisce un oggetto contenente i volumi anagrafica e dati e tutti i descrittori delle loro dimension...
Definisce il livello verticale di un'osservazione.
Definisce l'intervallo temporale di un'osservazione meteo.
Object describing a rectangular, homogeneous gridded dataset.
Class defining a real conversion function between units.
Definition of a physical variable in grib coding style.

Generated with Doxygen.