1 | #ifdef IBM |
---|
2 | @PROCESS NOOPT |
---|
3 | #endif |
---|
4 | MODULE bufrdata |
---|
5 | |
---|
6 | #ifdef HAS_BUFR |
---|
7 | |
---|
8 | USE toolspar_kind |
---|
9 | USE obs_utils |
---|
10 | USE obs_fbm |
---|
11 | |
---|
12 | IMPLICIT NONE |
---|
13 | |
---|
14 | CONTAINS |
---|
15 | |
---|
16 | SUBROUTINE read_bufrfile( cdfilename, inptsm, inpts1, inpuvm, inpuv1 ) |
---|
17 | !!--------------------------------------------------------------------- |
---|
18 | !! |
---|
19 | !! ** ROUTINE read_bufrfile ** |
---|
20 | !! |
---|
21 | !! ** Purpose : Read from file the profile BUFR observations. |
---|
22 | !! |
---|
23 | !! ** Method : The data file is a BUFR file. |
---|
24 | !! |
---|
25 | !! ** Action : |
---|
26 | !! |
---|
27 | !! History : |
---|
28 | !!---------------------------------------------------------------------- |
---|
29 | !! * Arguments |
---|
30 | CHARACTER(LEN=*) :: cdfilename ! Input file. |
---|
31 | TYPE(obfbdata) :: inptsm ! Output TS data structure (prof). |
---|
32 | TYPE(obfbdata) :: inpts1 ! Output TS data structure (single). |
---|
33 | TYPE(obfbdata) :: inpuvm ! Output UV data structure (prof). |
---|
34 | TYPE(obfbdata) :: inpuv1 ! Output UV data structure (single). |
---|
35 | !! * Local variables |
---|
36 | INTEGER :: iunit, iret |
---|
37 | INTEGER :: ibytesize, isubtype, itype |
---|
38 | INTEGER :: iyea, imon, iday, ihou, imin |
---|
39 | INTEGER :: jk |
---|
40 | REAL :: zlat,zlon |
---|
41 | CHARACTER(len=8) :: csign |
---|
42 | INTEGER, PARAMETER :: imaxlev = 300 |
---|
43 | INTEGER :: ilevts, iobstsm, iobsts1, ilevtsmax |
---|
44 | INTEGER :: ilevuv, iobsuvm, iobsuv1, ilevuvmax |
---|
45 | INTEGER :: iobs |
---|
46 | REAL, DIMENSION(imaxlev) :: zdts, ztem, zsal, zduv, zuv, zvv |
---|
47 | real :: zmdi |
---|
48 | |
---|
49 | ! Initialize data |
---|
50 | |
---|
51 | ibytesize = BIT_SIZE(1)/8 |
---|
52 | zmdi = REAL(fbrmdi) |
---|
53 | |
---|
54 | ! Open file |
---|
55 | CALL pbopen( iunit, cdfilename, 'r', iret ) |
---|
56 | IF ( iret /= 0 ) THEN |
---|
57 | WRITE(*,*)'Error opening file : ', TRIM(cdfilename) |
---|
58 | WRITE(*,*)'Error code : ', iret |
---|
59 | CALL abort |
---|
60 | ENDIF |
---|
61 | |
---|
62 | ! Count number of observations |
---|
63 | |
---|
64 | ilevtsmax = -1 |
---|
65 | ilevuvmax = -1 |
---|
66 | iobstsm = 0 |
---|
67 | iobsts1 = 0 |
---|
68 | iobsuvm = 0 |
---|
69 | iobsuv1 = 0 |
---|
70 | |
---|
71 | DO |
---|
72 | |
---|
73 | CALL decode_bufr( & |
---|
74 | & iunit, ibytesize, & |
---|
75 | & iyea, imon, iday, ihou, imin, & |
---|
76 | & zlat, zlon, & |
---|
77 | & csign, isubtype, & |
---|
78 | & imaxlev, ilevts, ilevuv, & |
---|
79 | & zdts, ztem, zsal, & |
---|
80 | & zduv, zuv, zvv, & |
---|
81 | & zmdi, iret ) |
---|
82 | |
---|
83 | ! No more observations if iret==-1 |
---|
84 | |
---|
85 | IF ( iret == -1 ) EXIT |
---|
86 | |
---|
87 | ! Ignore observations without levels |
---|
88 | |
---|
89 | IF ( ilevts > 0 ) THEN |
---|
90 | |
---|
91 | IF ( ilevts == 1 ) THEN |
---|
92 | iobsts1 = iobsts1 + 1 |
---|
93 | ELSE |
---|
94 | iobstsm = iobstsm + 1 |
---|
95 | IF ( ilevts > ilevtsmax ) ilevtsmax = ilevts |
---|
96 | ENDIF |
---|
97 | |
---|
98 | ENDIF |
---|
99 | |
---|
100 | IF ( ilevuv > 0 ) THEN |
---|
101 | |
---|
102 | IF ( ilevuv == 1 ) THEN |
---|
103 | iobsuv1 = iobsuv1 + 1 |
---|
104 | ELSE |
---|
105 | iobsuvm = iobsuvm + 1 |
---|
106 | IF ( ilevuv > ilevuvmax ) ilevuvmax = ilevuv |
---|
107 | ENDIF |
---|
108 | |
---|
109 | ENDIF |
---|
110 | |
---|
111 | ENDDO |
---|
112 | |
---|
113 | ! Setup data structures |
---|
114 | |
---|
115 | CALL init_obfbdata( inptsm ) |
---|
116 | CALL alloc_obfbdata( inptsm, 2, iobstsm, ilevtsmax, 0, 1, .FALSE. ) |
---|
117 | inptsm%cname(1) = 'POTM' |
---|
118 | inptsm%cname(2) = 'PSAL' |
---|
119 | inptsm%coblong(1) = 'Potential temperature' |
---|
120 | inptsm%coblong(2) = 'Practical salinity' |
---|
121 | inptsm%cobunit(1) = 'Degrees Celsius' |
---|
122 | inptsm%cobunit(2) = 'PSU' |
---|
123 | inptsm%cextname(1) = 'TEMP' |
---|
124 | inptsm%cextlong(1) = 'Insitu temperature' |
---|
125 | inptsm%cextunit(1) = 'Degrees Celsius' |
---|
126 | CALL init_obfbdata( inpts1 ) |
---|
127 | CALL alloc_obfbdata( inpts1, 2, iobsts1, 1, 0, 1, .FALSE. ) |
---|
128 | inpts1%cname(1) = 'POTM' |
---|
129 | inpts1%cname(2) = 'PSAL' |
---|
130 | inpts1%coblong(1) = 'Potential temperature' |
---|
131 | inpts1%coblong(2) = 'Practical salinity' |
---|
132 | inpts1%cobunit(1) = 'Degrees Celsius' |
---|
133 | inpts1%cobunit(2) = 'PSU' |
---|
134 | inpts1%cextname(1) = 'TEMP' |
---|
135 | inpts1%cextlong(1) = 'Insitu temperature' |
---|
136 | inpts1%cextunit(1) = 'Degrees Celsius' |
---|
137 | CALL init_obfbdata( inpuvm ) |
---|
138 | CALL alloc_obfbdata( inpuvm, 2, iobsuvm, ilevuvmax, 0, 0, .FALSE. ) |
---|
139 | inpuvm%cname(1) = 'UVEL' |
---|
140 | inpuvm%cname(2) = 'VVEL' |
---|
141 | inpuvm%coblong(1) = 'Zonal current' |
---|
142 | inpuvm%coblong(2) = 'Meridional current' |
---|
143 | inpuvm%cobunit(1) = 'Meters per second' |
---|
144 | inpuvm%cobunit(2) = 'Meters per second' |
---|
145 | CALL init_obfbdata( inpuv1 ) |
---|
146 | CALL alloc_obfbdata( inpuv1, 2, iobsuv1, 1, 0, 0, .FALSE. ) |
---|
147 | inpuv1%cname(1) = 'UVEL' |
---|
148 | inpuv1%cname(2) = 'VVEL' |
---|
149 | inpuv1%coblong(1) = 'Zonal current' |
---|
150 | inpuv1%coblong(2) = 'Meridional current' |
---|
151 | inpuv1%cobunit(1) = 'Meters per second' |
---|
152 | inpuv1%cobunit(2) = 'Meters per second' |
---|
153 | |
---|
154 | ! Rewind the file |
---|
155 | |
---|
156 | CALL pbseek( iunit, 0, 0, iret ) |
---|
157 | IF ( iret /= 0 ) THEN |
---|
158 | WRITE(*,*)'Error rewinding file : ',TRIM(cdfilename) |
---|
159 | WRITE(*,*)'Error code : ',iret |
---|
160 | CALL abort |
---|
161 | ENDIF |
---|
162 | |
---|
163 | ! Retrieve the observations and put them into the data structure |
---|
164 | |
---|
165 | iobstsm = 0 |
---|
166 | iobsts1 = 0 |
---|
167 | iobsuvm = 0 |
---|
168 | iobsuv1 = 0 |
---|
169 | iret = 0 |
---|
170 | inptsm%cdjuldref = '19500101000000' |
---|
171 | inpts1%cdjuldref = '19500101000000' |
---|
172 | inpuvm%cdjuldref = '19500101000000' |
---|
173 | inpuv1%cdjuldref = '19500101000000' |
---|
174 | iobs = 0 |
---|
175 | DO |
---|
176 | |
---|
177 | CALL decode_bufr( & |
---|
178 | & iunit, ibytesize, & |
---|
179 | & iyea, imon, iday, ihou, imin, & |
---|
180 | & zlat, zlon, & |
---|
181 | & csign, isubtype, & |
---|
182 | & imaxlev, ilevts, ilevuv, & |
---|
183 | & zdts, ztem, zsal, & |
---|
184 | & zduv, zuv, zvv, & |
---|
185 | & zmdi, iret ) |
---|
186 | |
---|
187 | ! No more observations if iret==-1 |
---|
188 | |
---|
189 | IF ( iret == -1 ) EXIT |
---|
190 | iobs = iobs + 1 |
---|
191 | |
---|
192 | ! Convert isubtype to something EN3 like. |
---|
193 | |
---|
194 | IF ( isubtype == 131 ) THEN |
---|
195 | itype = 820 |
---|
196 | ELSEIF ( isubtype == 132 ) THEN |
---|
197 | itype = 401 |
---|
198 | ELSEIF ( isubtype == 133 ) THEN |
---|
199 | IF (LEN_TRIM(CSIGN)>=7) THEN |
---|
200 | itype = 831 |
---|
201 | ELSE |
---|
202 | itype = 741 |
---|
203 | ENDIF |
---|
204 | ELSE |
---|
205 | WRITE(*,*)'Unknown subtype = ',isubtype |
---|
206 | itype = isubtype |
---|
207 | ENDIF |
---|
208 | |
---|
209 | ! TS data |
---|
210 | |
---|
211 | ! Ignore observations without levels |
---|
212 | |
---|
213 | IF ( ilevts > 0 ) THEN |
---|
214 | |
---|
215 | IF ( ilevts==1 ) THEN |
---|
216 | |
---|
217 | iobsts1 = iobsts1 + 1 |
---|
218 | |
---|
219 | ! Position and time |
---|
220 | |
---|
221 | inpts1%kindex(iobsts1) = iobs |
---|
222 | inpts1%pphi(iobsts1) = zlat |
---|
223 | inpts1%plam(iobsts1) = zlon |
---|
224 | CALL greg2jul( 0, imin, ihou, iday, imon, iyea, & |
---|
225 | & inpts1%ptim(iobsts1) ) |
---|
226 | |
---|
227 | ! Call sign and type |
---|
228 | |
---|
229 | inpts1%cdwmo(iobsts1) = csign |
---|
230 | WRITE(inpts1%cdtyp(iobsts1),'(I4.4)') itype |
---|
231 | |
---|
232 | ! Depth, Salinity and Insitu Temperature |
---|
233 | |
---|
234 | DO jk = 1, ilevts |
---|
235 | |
---|
236 | inpts1%pdep(jk,iobsts1) = zdts(jk) |
---|
237 | inpts1%pext(jk,iobsts1,1) = ztem(jk) |
---|
238 | inpts1%pob(jk,iobsts1,2) = zsal(jk) |
---|
239 | |
---|
240 | ! If no insitu temperature set the QC flag to 4 |
---|
241 | |
---|
242 | IF ( ztem(jk) == zmdi ) THEN |
---|
243 | inpts1%ivlqc(jk,iobsts1,1) = 4 |
---|
244 | ELSE |
---|
245 | inpts1%ivlqc(jk,iobsts1,1) = 1 |
---|
246 | ENDIF |
---|
247 | |
---|
248 | ! If no salinity set the QC flag to 4 |
---|
249 | |
---|
250 | IF ( zsal(jk) == zmdi ) THEN |
---|
251 | inpts1%ivlqc(jk,iobsts1,2) = 4 |
---|
252 | ELSE |
---|
253 | inpts1%ivlqc(jk,iobsts1,2) = 1 |
---|
254 | ENDIF |
---|
255 | |
---|
256 | ENDDO |
---|
257 | |
---|
258 | ELSE |
---|
259 | |
---|
260 | iobstsm = iobstsm + 1 |
---|
261 | |
---|
262 | ! Position and time |
---|
263 | |
---|
264 | inptsm%kindex(iobstsm) = iobs |
---|
265 | inptsm%pphi(iobstsm) = zlat |
---|
266 | inptsm%plam(iobstsm) = zlon |
---|
267 | CALL greg2jul( 0, imin, ihou, iday, imon, iyea, & |
---|
268 | & inptsm%ptim(iobstsm) ) |
---|
269 | |
---|
270 | ! Call sign and type |
---|
271 | |
---|
272 | inptsm%cdwmo(iobstsm) = csign |
---|
273 | WRITE(inptsm%cdtyp(iobstsm),'(I4.4)') itype |
---|
274 | |
---|
275 | ! Depth, Salinity and Insitu Temperature |
---|
276 | |
---|
277 | DO jk = 1, ilevts |
---|
278 | |
---|
279 | inptsm%pdep(jk,iobstsm) = zdts(jk) |
---|
280 | inptsm%pext(jk,iobstsm,1) = ztem(jk) |
---|
281 | inptsm%pob(jk,iobstsm,2) = zsal(jk) |
---|
282 | |
---|
283 | ! If no insitu temperature set the QC flag to 4 |
---|
284 | |
---|
285 | IF ( ztem(jk) == zmdi ) THEN |
---|
286 | inptsm%ivlqc(jk,iobstsm,1) = 4 |
---|
287 | ELSE |
---|
288 | inptsm%ivlqc(jk,iobstsm,1) = 1 |
---|
289 | ENDIF |
---|
290 | |
---|
291 | ! If no salinity set the QC flag to 4 |
---|
292 | |
---|
293 | IF ( zsal(jk) == zmdi ) THEN |
---|
294 | inptsm%ivlqc(jk,iobstsm,2) = 4 |
---|
295 | ELSE |
---|
296 | inptsm%ivlqc(jk,iobstsm,2) = 1 |
---|
297 | ENDIF |
---|
298 | |
---|
299 | ENDDO |
---|
300 | |
---|
301 | ENDIF |
---|
302 | |
---|
303 | ENDIF |
---|
304 | |
---|
305 | ! UV data |
---|
306 | |
---|
307 | ! Ignore observations without levels |
---|
308 | |
---|
309 | IF ( ilevuv > 0 ) THEN |
---|
310 | |
---|
311 | ! Skip surface only observations (TBC). |
---|
312 | |
---|
313 | IF ( ilevuv == 1 ) THEN |
---|
314 | |
---|
315 | iobsuv1 = iobsuv1 + 1 |
---|
316 | |
---|
317 | ! Position and time |
---|
318 | |
---|
319 | inpuv1%kindex(iobsuv1) = iobs |
---|
320 | inpuv1%pphi(iobsuv1) = zlat |
---|
321 | inpuv1%plam(iobsuv1) = zlon |
---|
322 | CALL greg2jul( 0, imin, ihou, iday, imon, iyea, & |
---|
323 | & inpuv1%ptim(iobsuv1) ) |
---|
324 | |
---|
325 | ! Call sign and type |
---|
326 | |
---|
327 | inpuv1%cdwmo(iobsuv1) = csign |
---|
328 | WRITE(inpuv1%cdtyp(iobsuv1),'(I4.4)') itype |
---|
329 | |
---|
330 | ! Depth, Salinity and Insitu Temperature |
---|
331 | |
---|
332 | DO jk = 1, ilevuv |
---|
333 | |
---|
334 | inpuv1%pdep(jk,iobsuv1) = zduv(jk) |
---|
335 | inpuv1%pob(jk,iobsuv1,1) = zuv(jk) |
---|
336 | inpuv1%pob(jk,iobsuv1,2) = zvv(jk) |
---|
337 | |
---|
338 | ! If no insitu temperature set the QC flag to 4 |
---|
339 | |
---|
340 | IF ( zuv(jk) == zmdi ) THEN |
---|
341 | inpuv1%ivlqc(jk,iobsuv1,1) = 4 |
---|
342 | ELSE |
---|
343 | inpuv1%ivlqc(jk,iobsuv1,1) = 1 |
---|
344 | ENDIF |
---|
345 | |
---|
346 | ! If no v velocity set the QC flag to 4 |
---|
347 | |
---|
348 | IF ( zvv(jk) == zmdi ) THEN |
---|
349 | inpuv1%ivlqc(jk,iobsuv1,2) = 4 |
---|
350 | ELSE |
---|
351 | inpuv1%ivlqc(jk,iobsuv1,2) = 1 |
---|
352 | ENDIF |
---|
353 | |
---|
354 | ENDDO |
---|
355 | |
---|
356 | ELSE |
---|
357 | |
---|
358 | iobsuvm = iobsuvm + 1 |
---|
359 | |
---|
360 | ! Position and time |
---|
361 | |
---|
362 | inpuvm%kindex(iobsuvm) = iobs |
---|
363 | inpuvm%pphi(iobsuvm) = zlat |
---|
364 | inpuvm%plam(iobsuvm) = zlon |
---|
365 | CALL greg2jul( 0, imin, ihou, iday, imon, iyea, & |
---|
366 | & inpuvm%ptim(iobsuvm) ) |
---|
367 | |
---|
368 | ! Call sign and type |
---|
369 | |
---|
370 | inpuvm%cdwmo(iobsuvm) = csign |
---|
371 | WRITE(inpuvm%cdtyp(iobsuvm),'(I4.4)') itype |
---|
372 | |
---|
373 | ! Depth, Salinity and Insitu Temperature |
---|
374 | |
---|
375 | DO jk = 1, ilevuv |
---|
376 | |
---|
377 | inpuvm%pdep(jk,iobsuvm) = zduv(jk) |
---|
378 | inpuvm%pob(jk,iobsuvm,1) = zuv(jk) |
---|
379 | inpuvm%pob(jk,iobsuvm,2) = zvv(jk) |
---|
380 | |
---|
381 | ! If no insitu temperature set the QC flag to 4 |
---|
382 | |
---|
383 | IF ( zuv(jk) == zmdi ) THEN |
---|
384 | inpuvm%ivlqc(jk,iobsuvm,1) = 4 |
---|
385 | ELSE |
---|
386 | inpuvm%ivlqc(jk,iobsuvm,1) = 1 |
---|
387 | ENDIF |
---|
388 | |
---|
389 | ! If no v velocity set the QC flag to 4 |
---|
390 | |
---|
391 | IF ( zvv(jk) == zmdi ) THEN |
---|
392 | inpuvm%ivlqc(jk,iobsuvm,2) = 4 |
---|
393 | ELSE |
---|
394 | inpuvm%ivlqc(jk,iobsuvm,2) = 1 |
---|
395 | ENDIF |
---|
396 | |
---|
397 | ENDDO |
---|
398 | |
---|
399 | ENDIF |
---|
400 | |
---|
401 | |
---|
402 | ENDIF |
---|
403 | |
---|
404 | |
---|
405 | ENDDO |
---|
406 | |
---|
407 | ! Close the file |
---|
408 | |
---|
409 | CALL pbclose( iunit, iret ) |
---|
410 | IF ( iret /= 0) THEN |
---|
411 | WRITE(*,*)'Error closing file : ',TRIM(cdfilename) |
---|
412 | WRITE(*,*)'Error code : ',iret |
---|
413 | CALL abort |
---|
414 | ENDIF |
---|
415 | |
---|
416 | ! Set the QC flags which can not be read from BUFR to 0 |
---|
417 | |
---|
418 | inptsm%ipqc(:) = 0 |
---|
419 | inptsm%ipqcf(:,:) = 0 |
---|
420 | inptsm%ivqc(:,:) = 0 |
---|
421 | inptsm%ivqcf(:,:,:) = 0 |
---|
422 | inptsm%ivqc(:,:) = 0 |
---|
423 | inptsm%ivqcf(:,:,: ) = 0 |
---|
424 | inptsm%ioqc(:) = 0 |
---|
425 | inptsm%ioqcf(:,:) = 0 |
---|
426 | inptsm%idqc(:,:) = 0 |
---|
427 | inptsm%idqcf(:,:,:) = 0 |
---|
428 | inptsm%ivlqcf(:,:,:,:) = 0 |
---|
429 | inpts1%ipqc(:) = 0 |
---|
430 | inpts1%ipqcf(:,:) = 0 |
---|
431 | inpts1%ivqc(:,:) = 0 |
---|
432 | inpts1%ivqcf(:,:,:) = 0 |
---|
433 | inpts1%ivqc(:,:) = 0 |
---|
434 | inpts1%ivqcf(:,:,: ) = 0 |
---|
435 | inpts1%ioqc(:) = 0 |
---|
436 | inpts1%ioqcf(:,:) = 0 |
---|
437 | inpts1%idqc(:,:) = 0 |
---|
438 | inpts1%idqcf(:,:,:) = 0 |
---|
439 | inpts1%ivlqcf(:,:,:,:) = 0 |
---|
440 | inpuvm%ipqc(:) = 0 |
---|
441 | inpuvm%ipqcf(:,:) = 0 |
---|
442 | inpuvm%ivqc(:,:) = 0 |
---|
443 | inpuvm%ivqcf(:,:,:) = 0 |
---|
444 | inpuvm%ivqc(:,:) = 0 |
---|
445 | inpuvm%ivqcf(:,:,: ) = 0 |
---|
446 | inpuvm%ioqc(:) = 0 |
---|
447 | inpuvm%ioqcf(:,:) = 0 |
---|
448 | inpuvm%idqc(:,:) = 0 |
---|
449 | inpuvm%idqcf(:,:,:) = 0 |
---|
450 | inpuvm%ivlqcf(:,:,:,:) = 0 |
---|
451 | inpuv1%ipqc(:) = 0 |
---|
452 | inpuv1%ipqcf(:,:) = 0 |
---|
453 | inpuv1%ivqc(:,:) = 0 |
---|
454 | inpuv1%ivqcf(:,:,:) = 0 |
---|
455 | inpuv1%ivqc(:,:) = 0 |
---|
456 | inpuv1%ivqcf(:,:,: ) = 0 |
---|
457 | inpuv1%ioqc(:) = 0 |
---|
458 | inpuv1%ioqcf(:,:) = 0 |
---|
459 | inpuv1%idqc(:,:) = 0 |
---|
460 | inpuv1%idqcf(:,:,:) = 0 |
---|
461 | inpuv1%ivlqcf(:,:,:,:) = 0 |
---|
462 | |
---|
463 | END SUBROUTINE read_bufrfile |
---|
464 | |
---|
465 | SUBROUTINE decode_bufr( kunit, kwsize, & |
---|
466 | & kyea, kmon, kday, khou, kmin, & |
---|
467 | & plat,plon, csign, ksty, & |
---|
468 | & kmlevs, klevts, klevve, & |
---|
469 | & pdts, ptem, psal, & |
---|
470 | & pdve, puve, pvve, & |
---|
471 | & pvalnull, ierr) |
---|
472 | !!--------------------------------------------------------------------- |
---|
473 | !! |
---|
474 | !! ** ROUTINE decode_bufr ** |
---|
475 | !! |
---|
476 | !! ** Purpose : Decode a singe BUFR record for ocean data |
---|
477 | !! |
---|
478 | !! ** Method : Call to BUFREX/BUSEL in EMOS lib. |
---|
479 | !! |
---|
480 | !! ** Action : |
---|
481 | !! |
---|
482 | !! History : (??-??) A. Vidard . ODASYS version. |
---|
483 | !! (08-12) K. Mogensen. NEMOVAR version |
---|
484 | !!---------------------------------------------------------------------- |
---|
485 | !! * Arguments |
---|
486 | INTEGER, INTENT(IN) :: kunit ! Bufr file unit number |
---|
487 | INTEGER, INTENT(IN) :: kwsize ! No of bytes in a word |
---|
488 | INTEGER, INTENT(OUT) :: & ! Date |
---|
489 | & kyea, kmon, kday, khou, kmin |
---|
490 | REAL, INTENT(OUT) :: & ! Position |
---|
491 | & plat, plon |
---|
492 | CHARACTER(LEN=*), INTENT(OUT) :: & |
---|
493 | & csign ! Call sign |
---|
494 | INTEGER, INTENT(OUT) :: ksty ! Subtype |
---|
495 | INTEGER, INTENT(IN) :: kmlevs ! Maximum number of levels |
---|
496 | INTEGER, INTENT(OUT) :: klevts ! Actual number of levels |
---|
497 | INTEGER, INTENT(OUT) :: klevve ! Actual number of levels |
---|
498 | REAL, DIMENSION(kmlevs), INTENT(OUT) :: & |
---|
499 | & pdts, & ! Depths for t,s |
---|
500 | & ptem, & ! Insitu temperature |
---|
501 | & psal, & ! Salinity |
---|
502 | & pdve, & ! Depths for currents |
---|
503 | & puve, & ! U velocity |
---|
504 | & pvve ! V velocity |
---|
505 | REAL, INTENT(IN) :: pvalnull ! Missing value |
---|
506 | INTEGER, INTENT(INOUT) :: ierr ! Error code |
---|
507 | !! * Define BUFR parameters |
---|
508 | INTEGER, PARAMETER :: imax_bufr_bytes=80000 ! Max bufr length in bytes |
---|
509 | INTEGER, PARAMETER :: imax_elements=20000 ! Kax no of elements |
---|
510 | INTEGER, PARAMETER :: imax_values=20000 ! Max no of values |
---|
511 | INTEGER, PARAMETER :: imax_sec0=3 ! Length of section 0 |
---|
512 | INTEGER, PARAMETER :: imax_sec1=40 ! Length of section 1 |
---|
513 | INTEGER, PARAMETER :: imax_sec2=64 ! Length of section 2 |
---|
514 | INTEGER, PARAMETER :: imax_sec3=4 ! Length of section 3 |
---|
515 | INTEGER, PARAMETER :: imax_sec4=2 ! Length of section 4 |
---|
516 | INTEGER, PARAMETER :: idim_ksup=9 ! Length of array ksup |
---|
517 | !! * Define local variables |
---|
518 | INTEGER :: ksup(idim_ksup) ! Expanded bufr additional information |
---|
519 | INTEGER :: ksec0(imax_sec0) ! Expanded bufr section 0 |
---|
520 | INTEGER :: ksec1(imax_sec1) ! Expanded bufr section 1 |
---|
521 | INTEGER :: ksec2(imax_sec2) ! Expanded bufr section 2 |
---|
522 | INTEGER :: ksec3(imax_sec3) ! Expanded bufr section 3 |
---|
523 | INTEGER :: ksec4(imax_sec4) ! Expanded bufr section 4 |
---|
524 | INTEGER :: kdtlen,ktdexl ! Bufr lengths |
---|
525 | INTEGER :: & ! BUSEL variables |
---|
526 | & ktdlst(imax_elements), ktdexp(imax_elements) |
---|
527 | INTEGER :: kbuff(imax_bufr_bytes) ! Bufr message array |
---|
528 | REAL :: values(imax_values) ! Bufr message values |
---|
529 | CHARACTER(LEN=64) :: & |
---|
530 | & cnames(imax_elements), & ! Element names |
---|
531 | & cunits(imax_elements), & ! Element units |
---|
532 | & cvals(imax_values) ! Character values |
---|
533 | REAL, DIMENSION((kmlevs)) :: & |
---|
534 | & zspd, & ! Current speed |
---|
535 | & zdir ! Current direction |
---|
536 | INTEGER :: ilenbuf ! Word len. of compressed bufr message |
---|
537 | INTEGER :: kerr ! Error return code |
---|
538 | INTEGER :: ibytes ! Actual no. of bytes in bufr message |
---|
539 | INTEGER :: ji ! Loop counter |
---|
540 | INTEGER :: icur ! Pointer to current data |
---|
541 | ! Define physical constants |
---|
542 | REAL, PARAMETER :: zero_celcius = 273.15 |
---|
543 | |
---|
544 | ! Initialize variables |
---|
545 | |
---|
546 | cvals(1) = ' ' |
---|
547 | |
---|
548 | DO ji = 1, kmlevs |
---|
549 | pdts(ji) = pvalnull |
---|
550 | ptem(ji) = pvalnull |
---|
551 | psal(ji) = pvalnull |
---|
552 | puve(ji) = pvalnull |
---|
553 | pvve(ji) = pvalnull |
---|
554 | END DO |
---|
555 | |
---|
556 | DO ji = 1, imax_elements |
---|
557 | ktdlst(ji) = 0 |
---|
558 | ktdexp(ji) = 0 |
---|
559 | END DO |
---|
560 | values(:) = pvalnull |
---|
561 | |
---|
562 | ! Read bufr message |
---|
563 | |
---|
564 | kerr = 0 |
---|
565 | CALL pbbufr( kunit, kbuff, imax_bufr_bytes, ibytes, kerr) |
---|
566 | IF ( kerr == -1 ) THEN |
---|
567 | ierr=-1 |
---|
568 | RETURN |
---|
569 | ELSE IF ( kerr < -1 ) THEN |
---|
570 | WRITE(*,*)'ERROR in PBBUFR with code ',kerr |
---|
571 | CALL abort |
---|
572 | ENDIF |
---|
573 | |
---|
574 | ! Decode and expand bufr message |
---|
575 | |
---|
576 | ! Calculate length of bufr message in words |
---|
577 | ilenbuf = ibytes/kwsize+1 |
---|
578 | |
---|
579 | ! Expand bufr message |
---|
580 | |
---|
581 | kerr=0 |
---|
582 | CALL bufrex( ilenbuf, kbuff, ksup, ksec0, ksec1, ksec2, ksec3, ksec4, & |
---|
583 | & imax_elements, cnames, cunits, imax_values, values, cvals, & |
---|
584 | & kerr) |
---|
585 | IF ( kerr /= 0 ) THEN |
---|
586 | WRITE(*,*)'ERROR in BUFREX with code ',kerr |
---|
587 | klevts = 0 |
---|
588 | klevve = 0 |
---|
589 | RETURN |
---|
590 | ENDIF |
---|
591 | CALL busel( kdtlen, ktdlst, ktdexl, ktdexp, kerr ) |
---|
592 | IF ( kerr /= 0 ) THEN |
---|
593 | WRITE(*,*)'ERROR in BUSEL with code ',kerr |
---|
594 | klevts = 0 |
---|
595 | klevve = 0 |
---|
596 | RETURN |
---|
597 | ENDIF |
---|
598 | |
---|
599 | ! Set oceanographic variables |
---|
600 | |
---|
601 | ksty=ksec1(7) |
---|
602 | kyea=ksec1(9) |
---|
603 | |
---|
604 | ! Patch for 2 digit only date (only matter if > 2000) |
---|
605 | |
---|
606 | IF ( kyea < 40 ) kyea = 2000 + kyea |
---|
607 | |
---|
608 | SELECT CASE(ksty) |
---|
609 | |
---|
610 | CASE(131) |
---|
611 | |
---|
612 | ! DRIBU |
---|
613 | |
---|
614 | kyea = values(5) |
---|
615 | kmon = values(6) |
---|
616 | kday = values(7) |
---|
617 | khou = values(8) |
---|
618 | kmin = values(9) |
---|
619 | plat = values(10) |
---|
620 | plon = values(11) |
---|
621 | |
---|
622 | ! Check that this dribu contains salinity |
---|
623 | |
---|
624 | IF ( ktdexp(16) == 2033 ) THEN |
---|
625 | |
---|
626 | ! Get number of levels |
---|
627 | IF ( values(17) < 100000. ) THEN |
---|
628 | klevts = INT( values(17) ) |
---|
629 | ELSE |
---|
630 | ! Patch for nasty data.. without data (AV 200307) |
---|
631 | klevts = 0 |
---|
632 | ENDIF |
---|
633 | |
---|
634 | ELSE |
---|
635 | |
---|
636 | klevts = 0 |
---|
637 | |
---|
638 | ENDIF |
---|
639 | |
---|
640 | ! Patch for missing WMO number (AV 200307) |
---|
641 | IF ( values(1) >= 9999999. ) THEN |
---|
642 | csign = 'XXX' |
---|
643 | ELSE |
---|
644 | IF ( INT( values(1) ) .LT. 100000 ) THEN |
---|
645 | WRITE(csign(1:5),'(i5)') INT( values(1) ) |
---|
646 | ELSE |
---|
647 | WRITE(csign, '(i8)') INT( values(1) ) |
---|
648 | ENDIF |
---|
649 | ENDIF |
---|
650 | |
---|
651 | CALL bufr_get_data( 7062, imax_elements, ktdexp, & |
---|
652 | & imax_values, values, klevts, 18, 3, pdts ) |
---|
653 | |
---|
654 | CALL bufr_get_data( 22043, imax_elements, ktdexp, & |
---|
655 | & imax_values, values, klevts, 19, 3, ptem ) |
---|
656 | |
---|
657 | CALL bufr_get_data( 22062, imax_elements, ktdexp, & |
---|
658 | & imax_values, values, klevts, 20, 3, psal ) |
---|
659 | |
---|
660 | IF ( ktdexp(15) == 2031 ) THEN |
---|
661 | |
---|
662 | klevve = INT( values(16) ) |
---|
663 | |
---|
664 | CALL bufr_get_data( 7062, imax_elements, ktdexp, & |
---|
665 | & imax_values, values, klevve, & |
---|
666 | & 17, 3, pdve ) |
---|
667 | |
---|
668 | CALL bufr_get_data( 22004, imax_elements, ktdexp, & |
---|
669 | & imax_values, values, klevve, & |
---|
670 | & 18, 3, zdir ) |
---|
671 | |
---|
672 | CALL bufr_get_data( 22031, imax_elements, ktdexp, & |
---|
673 | & imax_values, values, klevve, & |
---|
674 | & 19, 3, zspd ) |
---|
675 | |
---|
676 | ELSE |
---|
677 | |
---|
678 | klevve = 0 |
---|
679 | |
---|
680 | ENDIF |
---|
681 | |
---|
682 | |
---|
683 | CASE(132) |
---|
684 | |
---|
685 | ! BATHY |
---|
686 | |
---|
687 | kyea = values(2) |
---|
688 | kmon = values(3) |
---|
689 | kday = values(4) |
---|
690 | khou = values(5) |
---|
691 | kmin = values(6) |
---|
692 | plat = values(7) |
---|
693 | plon = values(8) |
---|
694 | |
---|
695 | csign = cvals(1) |
---|
696 | |
---|
697 | klevts = INT( values(10) ) |
---|
698 | |
---|
699 | CALL bufr_get_data( 7062, imax_elements, ktdexp, & |
---|
700 | & imax_values, values, klevts, 11, 2, pdts ) |
---|
701 | |
---|
702 | CALL bufr_get_data( 22042, imax_elements, ktdexp, & |
---|
703 | & imax_values, values, klevts, 12, 2, ptem ) |
---|
704 | |
---|
705 | psal=pvalnull |
---|
706 | |
---|
707 | klevve = 0 |
---|
708 | |
---|
709 | CASE(133) |
---|
710 | |
---|
711 | ! TESAC |
---|
712 | |
---|
713 | kyea = values(2) |
---|
714 | kmon = values(3) |
---|
715 | kday = values(4) |
---|
716 | khou = values(5) |
---|
717 | kmin = values(6) |
---|
718 | plat = values(7) |
---|
719 | plon = values(8) |
---|
720 | |
---|
721 | csign = cvals(1) |
---|
722 | |
---|
723 | klevts = INT( values(11) ) |
---|
724 | |
---|
725 | CALL bufr_get_data( 7062, imax_elements, ktdexp, & |
---|
726 | & imax_values, values, klevts, 12, 3, pdts ) |
---|
727 | |
---|
728 | CALL bufr_get_data( 22043, imax_elements, ktdexp, & |
---|
729 | & imax_values, values, klevts, 13, 3, ptem ) |
---|
730 | |
---|
731 | CALL bufr_get_data( 22062, imax_elements, ktdexp, & |
---|
732 | & imax_values, values, klevts, 14, 3, psal ) |
---|
733 | |
---|
734 | icur = 14 + 3 * (klevts - 1 ) + 1 |
---|
735 | |
---|
736 | IF ( ktdexp(icur) == 2031 ) THEN |
---|
737 | |
---|
738 | klevve = INT( values(icur+1) ) |
---|
739 | |
---|
740 | CALL bufr_get_data( 7062, imax_elements, ktdexp, & |
---|
741 | & imax_values, values, klevve, & |
---|
742 | & icur + 2, 3, pdve ) |
---|
743 | |
---|
744 | CALL bufr_get_data( 22004, imax_elements, ktdexp, & |
---|
745 | & imax_values, values, klevve, & |
---|
746 | & icur + 3, 3, zdir ) |
---|
747 | |
---|
748 | CALL bufr_get_data( 22031, imax_elements, ktdexp, & |
---|
749 | & imax_values, values, klevve, & |
---|
750 | & icur + 4, 3, zspd ) |
---|
751 | |
---|
752 | ELSE |
---|
753 | |
---|
754 | klevve = 0 |
---|
755 | |
---|
756 | ENDIF |
---|
757 | |
---|
758 | CASE DEFAULT |
---|
759 | |
---|
760 | WRITE(*,*)'ERROR: SUBTYPE ',ksty, & |
---|
761 | & ' is not recognised as an Ocean type' |
---|
762 | CALL abort |
---|
763 | |
---|
764 | END SELECT |
---|
765 | |
---|
766 | ! Convert data units |
---|
767 | |
---|
768 | DO ji = 1, klevts |
---|
769 | ! Convert temperature from Kelvin to Celcius |
---|
770 | IF (( ptem(ji) > 250. ) .AND. ( ptem(ji) < 323. )) THEN |
---|
771 | ptem(ji) = ptem(ji) - zero_celcius |
---|
772 | ELSE |
---|
773 | ptem(ji) = pvalnull |
---|
774 | ENDIF |
---|
775 | ! Check for physical salinities |
---|
776 | IF (( psal(ji) < 0.0 ) .OR. ( psal(ji) > 50.0 )) THEN |
---|
777 | psal(ji) = pvalnull |
---|
778 | ENDIF |
---|
779 | END DO |
---|
780 | |
---|
781 | DO ji = 1, klevve |
---|
782 | ! Convert speed and direction to u,v |
---|
783 | IF (( zspd(ji) < 100 ).AND.( zdir(ji) < 361. )) THEN |
---|
784 | CALL spd_to_uv( zdir(ji), zspd(ji), puve(ji), pvve(ji) ) |
---|
785 | ELSE |
---|
786 | puve(ji) = pvalnull |
---|
787 | pvve(ji) = pvalnull |
---|
788 | ENDIF |
---|
789 | END DO |
---|
790 | |
---|
791 | END SUBROUTINE decode_bufr |
---|
792 | |
---|
793 | |
---|
794 | SUBROUTINE bufr_get_data( ktdval, ktdlen, ktdexp, & |
---|
795 | & kmax_values, zvalues, klevs, & |
---|
796 | & istart, istep, zovalues ) |
---|
797 | !!--------------------------------------------------------------------- |
---|
798 | !! |
---|
799 | !! ** ROUTINE bufr_get_data ** |
---|
800 | !! |
---|
801 | !! ** Purpose : Extract data from decode buffer message |
---|
802 | !! |
---|
803 | !! ** Method : |
---|
804 | !! |
---|
805 | !! ** Action : |
---|
806 | !! |
---|
807 | !! History : (??-??) A. Vidard . ODASYS version. |
---|
808 | !! (08-12) K. Mogensen. NEMOVAR version |
---|
809 | !!---------------------------------------------------------------------- |
---|
810 | !! * Arguments |
---|
811 | INTEGER, INTENT(IN) :: ktdval ! Expected descriptor |
---|
812 | INTEGER, INTENT(IN) :: ktdlen ! Size of descriptors |
---|
813 | INTEGER, INTENT(IN), DIMENSION(ktdlen) :: & |
---|
814 | & ktdexp ! Descriptors |
---|
815 | INTEGER, INTENT(IN) :: kmax_values ! Size of input array |
---|
816 | REAL, INTENT(IN) :: zvalues(kmax_values) ! Decoded values |
---|
817 | INTEGER, INTENT(IN) :: klevs ! Number of levels |
---|
818 | INTEGER, INTENT(IN) :: istart ! Start in input array |
---|
819 | INTEGER, INTENT(IN) :: istep ! Step between values |
---|
820 | REAL, INTENT(OUT) :: zovalues(klevs) ! Output values |
---|
821 | !! * Local variables |
---|
822 | INTEGER :: ji, ipos |
---|
823 | |
---|
824 | DO ji = 1, klevs |
---|
825 | ipos = istart + (ji-1) * istep |
---|
826 | IF ( ktdexp(ipos) /= ktdval ) THEN |
---|
827 | WRITE(*,*)'Problem decoding bufr data in bufr_get_data' |
---|
828 | WRITE(*,*)'Expected descriptor : ', ktdval |
---|
829 | WRITE(*,*)'Found descriptor : ', ktdexp(ipos) |
---|
830 | CALL abort |
---|
831 | ENDIF |
---|
832 | zovalues(ji) = zvalues( ipos ) |
---|
833 | END DO |
---|
834 | |
---|
835 | END SUBROUTINE bufr_get_data |
---|
836 | |
---|
837 | SUBROUTINE spd_to_uv( pdir, pspd, pu, pv ) |
---|
838 | !!---------------------------------------------------------------------- |
---|
839 | !! |
---|
840 | !! *** ROUTINE spd_to_uv *** |
---|
841 | !! |
---|
842 | !! ** Purpose : Convert ocean current direction and speed to uv. |
---|
843 | !! |
---|
844 | !! ** Method : According to WMO current directors are coded according |
---|
845 | !! to normal oceanographic convention. |
---|
846 | !! |
---|
847 | !! ** Action : |
---|
848 | !! |
---|
849 | !! ** References: |
---|
850 | !! http://www.wmo.int/pages/prog/www/WMOCodes/OperationalCodes.html |
---|
851 | !! |
---|
852 | !! History |
---|
853 | !!---------------------------------------------------------------------- |
---|
854 | !! * Modules used |
---|
855 | !! * Arguments |
---|
856 | REAL, INTENT(IN) :: & |
---|
857 | & pdir, & ! Current direction |
---|
858 | & pspd ! Current speed |
---|
859 | REAL, INTENT(OUT) :: & |
---|
860 | & pu, & ! Zonal corrent |
---|
861 | & pv ! Meridonal current |
---|
862 | !! * Local variables |
---|
863 | REAL(KIND=dp), PARAMETER :: rad = 3.141592653589793_dp / 180.0_dp |
---|
864 | REAL(KIND=dp) :: zdir |
---|
865 | |
---|
866 | zdir = pdir * rad |
---|
867 | |
---|
868 | pu = pspd * SIN( zdir ) |
---|
869 | pv = pspd * COS( zdir ) |
---|
870 | |
---|
871 | END SUBROUTINE spd_to_uv |
---|
872 | |
---|
873 | #include "ctl_stop.h90" |
---|
874 | |
---|
875 | #include "greg2jul.h90" |
---|
876 | |
---|
877 | #endif |
---|
878 | |
---|
879 | END MODULE bufrdata |
---|