New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
bufrdata.F90 in branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS – NEMO

source: branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/bufrdata.F90 @ 2893

Last change on this file since 2893 was 2893, checked in by djlea, 13 years ago

Adding obs tools to branch

File size: 27.6 KB
Line 
1#ifdef IBM
2@PROCESS NOOPT
3#endif
4MODULE bufrdata
5
6#ifdef HAS_BUFR
7
8   USE toolspar_kind
9   USE obs_utils
10   USE obs_fbm
11
12   IMPLICIT NONE
13
14CONTAINS
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
879END MODULE bufrdata
Note: See TracBrowser for help on using the repository browser.