source: flex_extract.git/source/fortran/rwGRIB2.f90 @ a17dbbe

dev
Last change on this file since a17dbbe was a17dbbe, checked in by Anne Philipp <anne.philipp@…>, 10 months ago

fixed bug

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