source: flex_extract.git/python/pythontest/TestInstallTar/test_untar/src/rwGRIB2.f90 @ 2fb99de

ctbtodev
Last change on this file since 2fb99de was 2fb99de, checked in by Anne Philipp <anne.philipp@…>, 6 years ago

introduced config with path definitions and changed py files accordingly; Installation works; some tests were added for tarball making; Problems in submission to ecgate

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