source: flexpart.git/preproc/src/rwGRIB2.f90 @ 5763793

10.4.1_peseiFPv9.3.1FPv9.3.1b_testingFPv9.3.2GFS_025bugfixes+enhancementsdevfp9.3.1-20161214-nc4grib2nc4_repairrelease-10release-10.4.1scaling-bugunivie
Last change on this file since 5763793 was 5763793, checked in by Anne Fouilloux <annefou@…>, 9 years ago

Add pre-processing programs and examples. This is still under development!
  • Property mode set to 100644
File size: 6.4 KB
Line 
1MODULE RWGRIB2
2
3CONTAINS
4
5  SUBROUTINE READLATLON(FILENAME,FELD,MAXL,MAXB,MLEVEL,MPAR)
6
7    USE GRIB_API
8
9    IMPLICIT NONE
10
11    integer                            ::  ifile
12    integer                            ::  iret
13    integer                            ::  n,mk,parid,nm
14    integer                            ::  i,k
15    integer,dimension(:),allocatable   ::  igrib
16    integer                            ::  numberOfPointsAlongAParallel
17    integer                            ::  numberOfPointsAlongAMeridian
18    real, dimension(:), allocatable    ::  values
19    integer                            ::  numberOfValues
20    real,dimension(maxl,maxb,mlevel)   ::  feld 
21    integer:: maxl,maxb,mlevel,mstride,mpar(:),irest,div
22    integer :: l(size(mpar))
23    character*(*):: filename                             
24
25    call grib_open_file(ifile, TRIM(FILENAME),'r')
26
27    ! count the messages in the file
28    call grib_count_in_file(ifile,n)
29    allocate(igrib(n))
30    igrib=-1
31
32    ! Load the messages from the file.
33    DO i=1,n
34       call grib_new_from_file(ifile,igrib(i), iret)
35    END DO
36
37    ! we can close the file
38    call grib_close_file(ifile)
39
40    nm=size(mpar)
41    div=mlevel/nm
42    l=0
43
44    ! Loop on all the messages in memory
45    iloop:  DO i=1,n
46       !      write(*,*) 'processing message number ',i
47       !     get as a integer
48       call grib_get(igrib(i),'numberOfPointsAlongAParallel', &
49            numberOfPointsAlongAParallel)
50
51       !     get as a integer
52       call grib_get(igrib(i),'numberOfPointsAlongAMeridian', &
53            numberOfPointsAlongAMeridian)
54
55       call grib_get(igrib(i),'numberOfVerticalCoordinateValues',mk)
56
57       call grib_get_size(igrib(i),'values',numberOfValues)
58       !      write(*,*) 'numberOfValues=',numberOfValues
59
60       allocate(values(numberOfValues), stat=iret)
61       !     get data values
62       call grib_get(igrib(i),'values',values)
63
64       call grib_get(igrib(i),'paramId',parid)
65
66       kloop:  do k=1,nm
67          if(parid .eq. mpar(k)) then
68             l(k)=l(k)+1
69             feld(:,:,(k-1)*div+l(k))=reshape(values,(/maxl,maxb/))
70             !         print*,(k-1)*div+l(k),parid
71             exit kloop
72          endif
73       enddo kloop
74       if(k .gt. nm .and. parid .ne. mpar(nm)) then
75          write(*,*) k,nm,parid,mpar(nm)
76          write(*,*) 'ERROR readlatlon: parameter ',parid,'is not',mpar
77          stop
78       endif
79
80       !      print*,i
81      deallocate(values)
82    END DO iloop
83    write(*,*) 'readlatlon: ',i-1,' records read'
84
85    DO i=1,n
86       call grib_release(igrib(i))
87    END DO
88
89    deallocate(igrib)
90
91  END SUBROUTINE READLATLON
92
93  SUBROUTINE WRITELATLON(iunit,igrib,FELD,MAXL,MAXB,MLEVEL,&
94       MLEVELIST,MPAR)
95
96    USE GRIB_API
97
98    IMPLICIT NONE
99
100    INTEGER IFIELD,MLEVEL,MNAUF,I,J,K,L,IERR,JOUT
101    INTEGER MPAR,MAXL,MAXB,LEVMIN,LEVMAX
102    INTEGER IUNIT,igrib
103    REAL ZSEC4(MAXL*MAXB)
104    REAL    FELD(MAXL,MAXB,MLEVEL)
105    CHARACTER*(*) MLEVELIST
106    INTEGER ILEVEL(MLEVEL),MLINDEX(MLEVEL+1),LLEN
107
108    ! parse MLEVELIST
109
110    LLEN=len(trim(MLEVELIST))
111    if(index(MLEVELIST,'to') .ne. 0 .or. index(MLEVELIST,'TO') .ne. 0) THEN
112       i=index(MLEVELIST,'/')
113       read(MLEVELIST(1:i-1),*) LEVMIN
114       i=index(MLEVELIST,'/',.true.)
115       read(MLEVELIST(i+1:LLEN),*) LEVMAX
116       l=0
117       do i=LEVMIN,LEVMAX
118          l=l+1
119          ILEVEL(l)=i
120       enddo
121    else
122       l=1
123       MLINDEX(1)=0
124       do i=1,LLEN
125          if(MLEVELIST(i:i) .eq. '/') THEN
126             l=l+1
127             MLINDEX(l)=i
128          endif
129       enddo
130       MLINDEX(l+1)=LLEN+1
131       do i=1,l
132          read(MLEVELIST(MLINDEX(i)+1:MLINDEX(i+1)-1),*) ILEVEL(i)
133       enddo
134    endif
135
136    DO k=1,l
137       call grib_set(igrib,"level",ILEVEL(k))
138       call grib_set(igrib,"paramId",MPAR)
139       if(l .ne. mlevel) then
140          zsec4(1:maxl*maxb)=RESHAPE(FELD(:,:,ILEVEL(k)),(/maxl*maxb/))
141       else
142          zsec4(1:maxl*maxb)=RESHAPE(FELD(:,:,k),(/maxl*maxb/))
143       endif
144       call grib_set(igrib,"values",zsec4)
145
146       call grib_write(igrib,iunit)
147
148    ENDDO
149
150  END SUBROUTINE WRITELATLON
151
152  SUBROUTINE READSPECTRAL(FILENAME,CXMN,MNAUF,MLEVEL,&
153       MAXLEV,MPAR)
154
155    USE GRIB_API
156
157    IMPLICIT NONE
158
159    character*(*),               intent(in)    :: filename     
160
161    integer                            ::  ifile
162    integer                            ::  iret
163    integer                            ::  n,mk,div,nm,k
164    integer                            ::  i,j,parid
165    integer,dimension(:),allocatable   ::  igrib
166    real, dimension(:), allocatable    ::  values
167    integer                            ::  numberOfValues,maxlev 
168    REAL:: CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
169    integer:: maxl,maxb,mlevel,mstride,mpar(:),mnauf,ioffset,ipar,ilev,l(size(mpar))
170
171    call grib_open_file(ifile, TRIM(FILENAME),'r')
172
173    ! count the messages in the file
174    call grib_count_in_file(ifile,n)
175    allocate(igrib(n))
176    igrib=-1
177
178    ! Load the messages from the file.
179    DO i=1,n
180       call grib_new_from_file(ifile,igrib(i), iret)
181    END DO
182
183    ! we can close the file
184    call grib_close_file(ifile)
185
186    l=0
187    ! Loop on all the messages in memory
188    iloop: DO i=1,n
189       ! write(*,*) 'processing message number ',i
190       !     get as a integer
191       call grib_get(igrib(i),'pentagonalResolutionParameterJ', j)
192
193       call grib_get_size(igrib(i),'values',numberOfValues)
194       !   write(*,*) 'numberOfValues=',numberOfValues
195
196       call grib_get(igrib(i),'numberOfVerticalCoordinateValues',mk)
197
198       call grib_get(igrib(i),'level',ilev)
199
200       allocate(values(numberOfValues), stat=iret)
201       !     get data values
202       call grib_get(igrib(i),'values',values)
203
204       call grib_get(igrib(i),'paramId',parid)
205       nm=size(mpar)
206       div=mlevel/nm
207       kloop:  do k=1,nm
208          if(parid .eq. mpar(k)) then
209             l(k)=l(k)+1
210             cxmn(:,(k-1)*div+l(k))=values(1:(MNAUF+1)*(MNAUF+2))
211             !         print*,(k-1)*div+l(k),parid
212             exit kloop
213          endif
214
215       enddo kloop
216       if(k .gt. nm .and. parid .ne. mpar(nm)) then
217          write(*,*) k,nm,parid,mpar(nm)
218          write(*,*) 'ERROR readspectral: parameter ',parid,'is not',mpar
219          stop
220       endif
221
222       !      print*,i
223      deallocate(values)
224
225    END DO iloop
226
227    write(*,*) 'readspectral: ',i-1,' records read'
228
229    DO i=1,n
230       call grib_release(igrib(i))
231    END DO
232
233    deallocate(igrib)
234
235  END SUBROUTINE READSPECTRAL
236
237END MODULE RWGRIB2
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG