[2893] | 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 |
---|