var/obsproc/src/module_write.F90

   1 MODULE module_write
   2
   3 !-----------------------------------------------------------------------------!
   4 ! Output observations for intput into MM5 3D-VAR
   5 !
   6 !-----------------------------------------------------------------------------!
   7 !  HISTORY:
   8 !
   9 !         F. VANDENBERGHE, March 2001
  10 !         01/13/2003 - Updated for Profiler obs.           S. R. H. Rizvi
  11 !
  12 !         02/04/2003 - Updated for Buoy     obs.           S. R. H. Rizvi
  13 !
  14 !         02/11/2003 - Reviewed and modified for Profiler
  15 !                      and Buoy obs.                       Y.-R. Guo
  16 !         11/11/2005 - Added output_prep to write conventional obs
  17 !                      in PREPBUFR if that is chosen       J. F. Drake
  18 !         06/30/2006 -   Updated for AIRS retrievals       Syed  RH  Rizvi
  19 !   References:
  20 !    http://www.emc.ncep.noaa.gov/mmb/data_processing/prepbufr.doc/document.htm
  21 !                                                                      /table_2.htm
  22 !------------------------------------------------------------------------------
  23
  24   USE module_type
  25   USE module_func
  26   USE module_date
  27   USE module_decoded
  28   USE module_map
  29
  30   INCLUDE 'constants.inc'
  31   INCLUDE '../../../inc/version_decl'
  32
  33   REAL (kind=8), parameter             :: bufrlib_missing = 10.D10
  34
  35 CONTAINS
  36
  37 !------------------------------------------------------------------------------!
  38 SUBROUTINE output_ssmi_31 (max_number_of_obs, obs, number_of_obs, index, &
  39                            nssmis, missing_r, time_analysis)
  40
  41 !-------------------------------------------------------------------------------
  42 ! Write observation at MM5 3D-VAR SYSTEM VERSION 2.0
  43 ! Output file is "obs_ascii.ccyymmddhh
  44 !
  45 ! F. VANDENBERGHE, March 2000
  46 !-------------------------------------------------------------------------------
  47
  48
  49   IMPLICIT NONE
  50
  51   INTEGER,                                      INTENT (in) :: max_number_of_obs
  52   TYPE (report), DIMENSION (max_number_of_obs), INTENT (inout) :: obs
  53   INTEGER,                                      INTENT (in) :: number_of_obs
  54   INTEGER,       DIMENSION (max_number_of_obs), INTENT (in) :: index
  55   REAL,                                         INTENT (in) :: missing_r
  56   CHARACTER (LEN =  19),                        INTENT (in) :: time_analysis
  57   INTEGER,                                      INTENT (in) :: nssmis
  58
  59   INTEGER                       :: n, loop_index, fm, i, ssmi_125, ssmi_126
  60   INTEGER                       :: nwrites1, nwrites2
  61   LOGICAL                       :: connected
  62   CHARACTER (LEN = 80)          :: filename1, filename2
  63   CHARACTER (LEN = 120)         :: fmt_info, &
  64                                    fmt_srfc, &
  65                                    fmt_each
  66   REAL                          :: rew_cross, rns_cross
  67 !------------------------------------------------------------------------------
  68
  69   rew_cross=real(nestjx(idd)-1)
  70   rns_cross=real(nestix(idd)-1)
  71
  72       if (nssmis == 0) then
  73         WRITE(0,'(A,/)') "No SSMI observations available."
  74         RETURN
  75       else
  76         ssmi_125 = 0
  77         ssmi_126 = 0
  78 counts: &
  79         DO n = 1, number_of_obs
  80         loop_index = index (n)
  81         IF (obs (loop_index)%info%discard ) THEN
  82           CYCLE  counts
  83         ELSE
  84           READ (obs (loop_index) % info % platform (4:6), '(I3)') fm
  85           if (fm == 125) THEN
  86             ssmi_125 = ssmi_125 + 1
  87           else if (fm == 126) THEN
  88             ssmi_126 = ssmi_126 + 1
  89           endif
  90         ENDIF
  91
  92         ENDDO counts
  93       endif
  94
  95 !------------------------------------------------------------------------------!
  96 !     B. PRINT OBS USING THE NEW STRUCUTURE
  97 !------------------------------------------------------------------------------!
  98
  99       WRITE (0,'(A)')  &
 100 '------------------------------------------------------------------------------'
 101
 102
 103 ! 0.  FORMAT
 104 ! ==========
 105
 106 !     fmt_info = '(a12,1x,a19,1x,a40,1x,i6,3(f12.3,11x),6x,a5)'
 107      fmt_info = '(a12,1x,a19,1x,a40,1x,i6,3(f12.3,11x),6x,a40)'
 108      fmt_srfc = '(7(:,f12.3,i4,f7.2))'
 109      fmt_each = '(3(f12.3,i4,f7.2),11x,3(f12.3,i4,f7.2),11x,1(f12.3,i4,f7.2))'
 110
 111
 112 ! 1. OPEN FILE FOR SSMI RETRIEVAL DATA OUTPUT
 113 ! ==========================================
 114
 115 ! 1.1 Open output file for retrievals in version 3.1 format
 116 !     ----------------------------------------------------
 117
 118       if (ssmi_125 > 0) then
 119 !        filename1 = 'obs_ssmi_retrieval.3dvar'
 120         write(filename1,'("obs_ssmi_retrieval_",a,".",a)') time_analysis, use_for
 121
 122         WRITE (0,'(A,A)') &
 123         "Write 3DVAR SSMI Retrieval obs in files: ", TRIM (filename1)
 124
 125         INQUIRE (UNIT = 125, OPENED = connected )
 126
 127         IF ( .NOT. connected ) THEN
 128           OPEN ( UNIT   = 125,          &
 129                  FILE   = filename1,     &
 130                  FORM   = 'FORMATTED',  &
 131                  ACCESS = 'SEQUENTIAL', &
 132                  STATUS = 'REPLACE')
 133         ENDIF
 134
 135         REWIND ( UNIT = 125)
 136
 137 ! 1.2 FILE HEADER FOR NEW FORMAT
 138 !     --------------------------
 139
 140 ! 1.2.1 TOTAL NUMBER OF STATIONS CONTAINED IN FILE AND MISSING VALUE
 141 !       ------------------------------------------------------------
 142
 143         WRITE  (UNIT = 125, FMT = '(A,I7,A)',ADVANCE='no' ) "SSMI  =",ssmi_125,","
 144         WRITE  (UNIT = 125, FMT = '((A,F8.0),A)') " MISS. = ",missing_r,","
 145
 146 ! 1.2.2 REFERENCE STATE INFO
 147 !       --------------------
 148
 149         WRITE (UNIT = 125, FMT = '(A,F7.2,A,F7.2,4(A,F7.2),A)') &
 150         "PHIC  =", phic,", XLONC =", xlonc,", TRUE1 =", truelat1,&
 151         ", TRUE2 =",truelat2,", XIM11 =",xim11(idd),", XJM11 =",xjm11(idd),","
 152         WRITE (UNIT = 125, FMT = '(2(A,F7.2),4(A,F7.0),A)') &
 153         "base_temp=",  ts0, ", base_lapse=", tlp, &
 154         ", PTOP  =",  ptop,", base_pres=",  ps0, &
 155         ", base_tropo_pres=", pis0, ", base_strat_temp=", tis0, ","
 156
 157 ! 1.2.3 DOMAINS INFO
 158 !       ------------
 159
 160         WRITE (UNIT = 125, FMT = '(5(A,I7),A)' ) &
 161         "IXC   =", ixc,", JXC   =", jxc,", IPROJ =", iproj,&
 162         ", IDD   =", idd,", MAXNES=",maxnes,","
 163         WRITE (UNIT = 125, FMT = '(A,10(:,I7,A))')  &
 164         "NESTIX=",(nestix (i),", ", i = 1, maxnes)
 165         WRITE (UNIT = 125, FMT = '(A,10(:,I7,A))')  &
 166         "NESTJX=",(nestjx (i),", ", i = 1, maxnes)
 167         WRITE (UNIT = 125, FMT = '(A,10(:,I7,A))')  &
 168         "NUMC  =",(numc   (i),", ", i = 1, maxnes)
 169         WRITE (UNIT = 125, FMT = '(A,10(:,F7.2,A))')&
 170         "DIS   =",(dis    (i),", ",i = 1, maxnes)
 171         WRITE (UNIT = 125, FMT = '(A,10(:,I7,A))')  &
 172         "NESTI =",(nesti  (i),", ", i = 1, maxnes)
 173         WRITE (UNIT = 125, FMT = '(A,10(:,I7,A))')  &
 174         "NESTJ =",(nestj  (i),", ", i = 1, maxnes)
 175
 176 ! 1.2.4 VARIABLE NAME AND UNITS
 177 !       -----------------------
 178
 179         WRITE (UNIT = 125, FMT = '(A)' ) &
 180         "INFO  = PLATFORM, DATE, NAME, LEVELS, LATITUDE, LONGITUDE, ELEVATION, ID."
 181         WRITE (UNIT = 125, FMT = '(A)' ) &
 182         "SRFC  = WIND SPEED, PW (DATA,QC,ERROR)."
 183
 184 ! 1.2.5 FORMATS
 185 !       -------
 186
 187        WRITE (UNIT = 125, FMT = '(2A)' ) 'INFO_FMT = ', TRIM (fmt_info)
 188        WRITE (UNIT = 125, FMT = '(2A)' ) 'SRFC_FMT = ', TRIM (fmt_srfc)
 189
 190 ! 1.2.6 END OF HEADER
 191 !       -------------
 192
 193         WRITE (UNIT = 125, FMT = '(A)') &
 194 "#------------------------------------------------------------------------------#"
 195
 196       endif
 197
 198 ! 1.3 Open output file for Tb in version 3.1 format
 199 !     --------------------------------------------
 200
 201       if (ssmi_126 > 0) then
 202 !        filename2 = 'obs_ssmi_Tb.3dvar'
 203         write(filename2,'("obs_ssmi_Tb_",a,".",a)') time_analysis, use_for
 204
 205         WRITE (0,'(A,A)') &
 206         "Write 3DVAR SSMI TB obs in files: ", TRIM (filename2)
 207
 208         INQUIRE (UNIT = 126, OPENED = connected )
 209
 210         IF ( .NOT. connected ) THEN
 211           OPEN ( UNIT   = 126,          &
 212                  FILE   = filename2,     &
 213                  FORM   = 'FORMATTED',  &
 214                  ACCESS = 'SEQUENTIAL', &
 215                  STATUS = 'REPLACE')
 216         ENDIF
 217
 218         REWIND ( UNIT = 126)
 219
 220 ! 1.4 FILE HEADER FOR NEW FORMAT
 221 !     --------------------------
 222
 223 ! 1.4.1 M TOTAL NUMBER OF STATIONS CONTAINED IN FILE AND ISSING VALUE FLAG
 224 !       ------------------------------------------------------------------
 225
 226         WRITE  (UNIT = 126, FMT = '(A,I7,A)',ADVANCE='no' ) "SSMI  =",ssmi_126,","
 227         WRITE  (UNIT = 126, FMT = '((A,F8.0),A)') " MISS. = ",missing_r,","
 228
 229 ! 1.4.2 REFERENCE STATE INFO
 230 !       --------------------
 231
 232         WRITE (UNIT = 126, FMT = '(A,F7.2,A,F7.2,4(A,F7.2),A)') &
 233         "PHIC  =", phic,", XLONC =", xlonc,", TRUE1 =", truelat1,&
 234         ", TRUE2 =",truelat2,", XIM11 =",xim11(idd),", XJM11 =",xjm11(idd),","
 235
 236         WRITE (UNIT = 126, FMT = '(2(A,F7.2),4(A,F7.0),A)') &
 237         "base_temp=",  ts0, ", base_lapse=", tlp, &
 238         ", PTOP  =",  ptop,", base_pres=",  ps0, &
 239         ", base_tropo_pres=", pis0, ", base_strat_temp=", tis0, ","
 240
 241 ! 1.4.3 DOMAINS INFO
 242 !       ------------
 243
 244         WRITE (UNIT = 126, FMT = '(5(A,I7),A)' ) &
 245         "IXC   =", ixc,", JXC   =", jxc,", IPROJ =", iproj,&
 246         ", IDD   =", idd,", MAXNES=",maxnes,","
 247         WRITE (UNIT = 126, FMT = '(A,10(:,I7,A))')  &
 248         "NESTIX=",(nestix (i),", ", i = 1, maxnes)
 249         WRITE (UNIT = 126, FMT = '(A,10(:,I7,A))')  &
 250         "NESTJX=",(nestjx (i),", ", i = 1, maxnes)
 251         WRITE (UNIT = 126, FMT = '(A,10(:,I7,A))')  &
 252         "NUMC  =",(numc   (i),", ", i = 1, maxnes)
 253         WRITE (UNIT = 126, FMT = '(A,10(:,F7.2,A))')&
 254         "DIS   =",(dis    (i),", ",i = 1, maxnes)
 255         WRITE (UNIT = 126, FMT = '(A,10(:,I7,A))')  &
 256         "NESTI =",(nesti  (i),", ", i = 1, maxnes)
 257         WRITE (UNIT = 126, FMT = '(A,10(:,I7,A))')  &
 258         "NESTJ =",(nestj  (i),", ", i = 1, maxnes)
 259
 260 ! 1.4.4 VARIABLE NAME AND UNITS
 261 !       -----------------------
 262
 263         WRITE (UNIT = 126, FMT = '(A)' ) &
 264         "INFO  = PLATFORM, DATE, NAME, LEVELS, LATITUDE, LONGITUDE, ELEVATION, ID."
 265         WRITE (UNIT = 126, FMT = '(A)' ) &
 266         "SRFC  = TB19V, TB19H, TB22V, TB37V, TB37H, TB85V, TB85H (DATA,QC,ERROR)."
 267
 268 ! 1.4.5 FORMATS
 269 !       -------
 270
 271         WRITE (UNIT = 126, FMT = '(2A)' ) 'INFO_FMT = ', TRIM (fmt_info)
 272         WRITE (UNIT = 126, FMT = '(2A)' ) 'SRFC_FMT = ', TRIM (fmt_srfc)
 273
 274 ! 1.4.6 END OF HEADER
 275 !       -------------
 276
 277         WRITE (UNIT = 126, FMT = '(A)') &
 278 "#------------------------------------------------------------------------------#"
 279
 280       endif
 281
 282 ! 2.  WRITE OBSERVATIONS
 283 ! ======================
 284
 285       nwrites1 = 0
 286       nwrites2 = 0
 287
 288
 289 ! 2.1  Loop over stations
 290 !      ------------------
 291
 292 stations: &
 293        DO n = 1, number_of_obs
 294
 295 ! 2.2 Index of obs
 296 !     ------------
 297
 298       loop_index = index (n)
 299
 300 ! 2.3 Check if station is valid
 301 !     -------------------------
 302
 303 stations_valid: &
 304       IF (obs (loop_index)%info%discard ) THEN
 305
 306           CYCLE  stations
 307
 308       ELSE stations_valid
 309
 310       READ (obs (loop_index) % info % platform (4:6), '(I3)') fm
 311
 312       if ((fm /= 125) .AND. (fm /= 126))  CYCLE stations
 313
 314 ! 2.4 Write station info
 315 !     ------------------
 316
 317       IF (fm == 125) THEN
 318
 319       WRITE (UNIT = 125, FMT = TRIM (fmt_info))         &
 320              obs (loop_index) % info     % platform,    &
 321              obs (loop_index) % valid_time % date_mm5,  &
 322              obs (loop_index) % location % name,        &
 323              obs (loop_index) % info % levels,          &
 324              obs (loop_index) % location % latitude,    &
 325              obs (loop_index) % location % longitude,   &
 326              obs (loop_index) % info     % elevation,   &
 327              obs (loop_index) % location % id
 328
 329              nwrites1 = nwrites1 + 1
 330
 331       ELSE IF (fm == 126) THEN
 332
 333       WRITE (UNIT = 126, FMT = TRIM (fmt_info))          &
 334              obs (loop_index) % info     % platform,    &
 335              obs (loop_index) % valid_time % date_mm5,  &
 336              obs (loop_index) % location % name,        &
 337              obs (loop_index) % info % levels,          &
 338              obs (loop_index) % location % latitude,    &
 339              obs (loop_index) % location % longitude,   &
 340              obs (loop_index) % info     % elevation,   &
 341              obs (loop_index) % location % id
 342
 343              nwrites2 = nwrites2 + 1
 344
 345       ENDIF
 346
 347 ! 2.4 Write surface info
 348 !     ------------------
 349
 350       IF (fm == 125) THEN
 351
 352       IF (ASSOCIATED (obs (loop_index) % surface)) THEN
 353
 354          if( domain_check_h .and. &
 355            (obs (loop_index) % location % xjc < 1.0 .or. &
 356             obs (loop_index) % location % xjc > rew_cross .or. &
 357             obs (loop_index) % location % yic < 1.0 .or. &
 358             obs (loop_index) % location % yic > rns_cross) ) then
 359             obs (loop_index) % ground  % pw   % qc = -88
 360          end if
 361
 362          WRITE (UNIT = 125, FMT = TRIM (fmt_srfc))                  &
 363                 obs (loop_index) % surface % meas % speed % data,  &
 364                 obs (loop_index) % surface % meas % speed % qc,    &
 365                 obs (loop_index) % surface % meas % speed % error, &
 366                 obs (loop_index) % ground  % pw   % data, &
 367                 obs (loop_index) % ground  % pw   % qc,   &
 368                 obs (loop_index) % ground  % pw   % error
 369
 370       ELSE
 371
 372          if( domain_check_h .and. &
 373            (obs (loop_index) % location % xjc < 1.0 .or. &
 374             obs (loop_index) % location % xjc > rew_cross .or. &
 375             obs (loop_index) % location % yic < 1.0 .or. &
 376             obs (loop_index) % location % yic > rns_cross) ) then
 377             obs (loop_index) % ground  % pw   % qc = -88
 378          end if
 379
 380          WRITE (UNIT = 125, FMT = TRIM (fmt_srfc))         &
 381                 missing_r, CEILING (missing_r/10000), 2.5, &
 382                 obs (loop_index) % ground % pw  % data,    &
 383                 obs (loop_index) % ground % pw  % qc,      &
 384                 obs (loop_index) % ground % pw  % error
 385
 386       ENDIF
 387
 388                 nwrites1 = nwrites1 + 1
 389
 390       ELSE IF (fm == 126) THEN
 391
 392          if( domain_check_h .and. &
 393            (obs (loop_index) % location % xjc < 1.0 .or. &
 394             obs (loop_index) % location % xjc > rew_cross .or. &
 395             obs (loop_index) % location % yic < 1.0 .or. &
 396             obs (loop_index) % location % yic > rns_cross) ) then
 397             obs (loop_index) % ground % tb19v % qc = -88
 398             obs (loop_index) % ground % tb19h % qc = -88
 399             obs (loop_index) % ground % tb22v % qc = -88
 400             obs (loop_index) % ground % tb37v % qc = -88
 401             obs (loop_index) % ground % tb37h % qc = -88
 402             obs (loop_index) % ground % tb85v % qc = -88
 403             obs (loop_index) % ground % tb85h % qc = -88
 404          end if
 405
 406          WRITE (UNIT = 126, FMT = TRIM (fmt_srfc))         &
 407          obs (loop_index) % ground % tb19v % data,  &
 408          obs (loop_index) % ground % tb19v % qc,    &
 409          obs (loop_index) % ground % tb19v % error, &
 410          obs (loop_index) % ground % tb19h % data,  &
 411          obs (loop_index) % ground % tb19h % qc,    &
 412          obs (loop_index) % ground % tb19h % error, &
 413          obs (loop_index) % ground % tb22v % data,  &
 414          obs (loop_index) % ground % tb22v % qc,    &
 415          obs (loop_index) % ground % tb22v % error, &
 416          obs (loop_index) % ground % tb37v % data,  &
 417          obs (loop_index) % ground % tb37v % qc,    &
 418          obs (loop_index) % ground % tb37v % error, &
 419          obs (loop_index) % ground % tb37h % data,  &
 420          obs (loop_index) % ground % tb37h % qc,    &
 421          obs (loop_index) % ground % tb37h % error, &
 422          obs (loop_index) % ground % tb85v % data,  &
 423          obs (loop_index) % ground % tb85v % qc,    &
 424          obs (loop_index) % ground % tb85v % error, &
 425          obs (loop_index) % ground % tb85h % data,  &
 426          obs (loop_index) % ground % tb85h % qc,    &
 427          obs (loop_index) % ground % tb85h % error
 428
 429          nwrites2 = nwrites2 + 1
 430
 431       ENDIF
 432
 433 ! 3.3 Go to next valid station
 434 !     -----------------
 435
 436       ENDIF stations_valid
 437
 438 ! 3.3 Go to next record
 439 !     -----------------
 440
 441       ENDDO stations
 442
 443
 444 ! 4.   CLOSE OUTPUT FILES
 445 ! ========================
 446
 447       CLOSE (UNIT = 125)
 448       CLOSE (UNIT = 126)
 449
 450
 451 ! 5.  PRINT DIAGNOSTIC
 452 ! =====================
 453
 454       WRITE (0, '(A)') ' '
 455
 456       if (ssmi_125 > 0) WRITE (0, '(A,I7,A,A)') &
 457       'Wrote ',nwrites1,' lines of data in file: ',TRIM (filename1)
 458
 459       if (ssmi_126 > 0) WRITE (0, '(A,I7,A,A)') &
 460       'Wrote ',nwrites2,' lines of data in file: ',TRIM (filename2)
 461
 462       WRITE (0, '(A)') ' '
 463
 464       RETURN
 465
 466 END SUBROUTINE output_ssmi_31
 467
 468 SUBROUTINE output_gts_31 (max_number_of_obs, obs, number_of_obs, windex,&
 469                          nsynops, nshipss, nmetars,                    &
 470                          npilots, nsounds, nsatems,                    &
 471                          nsatobs, naireps, ngpspws, ngpsztd, ngpsref,  &
 472                          ngpseph, nssmt1s, nssmt2s, nssmis,  ntovss,   &
 473                          nothers, namdars, nqscats, nprofls,nbuoyss,   &
 474                          nboguss, nairss, ntamdar, missing_r, time_analysis)
 475
 476 !-------------------------------------------------------------------------------
 477 ! Write observation at MM5 3D-VAR SYSTEM VERSION 2.0
 478 ! Output file is "obs_ascii.ccyymmddhh
 479 !
 480 ! F. VANDENBERGHE, March 2000
 481 !-------------------------------------------------------------------------------
 482
 483
 484   IMPLICIT NONE
 485
 486   INTEGER,                                      INTENT (in) :: max_number_of_obs
 487   TYPE (report), DIMENSION (max_number_of_obs), INTENT (inout) :: obs
 488   INTEGER,                                      INTENT (in) :: number_of_obs
 489   INTEGER,       DIMENSION (max_number_of_obs), INTENT (in) :: windex
 490   REAL,                                         INTENT (in) :: missing_r
 491   CHARACTER (LEN =  19),                        INTENT (in) :: time_analysis
 492   INTEGER,                                      INTENT (in) :: nsynops,nmetars,&
 493                                                                nshipss,nsounds,&
 494                                                                npilots,naireps,&
 495                                                                nsatems,nsatobs,&
 496                                                                ngpspws,nssmt1s,&
 497                                                                nssmt2s, nssmis,&
 498                                                                ntovss, namdars,&
 499                                                                nqscats,nothers,&
 500                                                                nprofls,nbuoyss,&
 501                                                                ngpsztd,ngpsref,&
 502                                                                ngpseph,nboguss,&
 503                                                                nairss,ntamdar
 504
 505   TYPE (measurement ) , POINTER :: current
 506   INTEGER                       :: loop_index
 507   INTEGER                       :: i, ii, n, ntotal, k_levels
 508   INTEGER                       :: nmultis, nsingles, nlevels, nwrites
 509   INTEGER                       :: is_sound, fm
 510   LOGICAL                       :: connected
 511   CHARACTER (LEN = 80)          :: filename
 512   CHARACTER (LEN = 120)         :: fmt_info, &
 513                                    fmt_srfc,  &
 514                                    fmt_each
 515   CHARACTER (LEN=40) :: SID
 516   INTEGER   ::    kx
 517   REAL                          :: val_slp, val_pw, pressure_error
 518   REAL                          :: val_u,  val_v,  val_p, val_t
 519   REAL                          :: val_td, val_rh, val_qv
 520 !------------------------------------------------------------------------------
 521   REAL                          :: rew_cross, rns_cross
 522   LOGICAL                       :: change_qc
 523
 524   rew_cross=real(nestjx(idd)-1)
 525   rns_cross=real(nestix(idd)-1)
 526 !------------------------------------------------------------------------------!
 527
 528       ntotal =   nsynops + nmetars + nshipss + &
 529                  nsounds + npilots + naireps + &
 530                  nsatems + nsatobs + ngpspws + &
 531                  nssmt1s + nssmt2s + ntovss  + &
 532                  namdars + nqscats + nprofls + &
 533                  nbuoyss + nothers + ngpsztd + &
 534                  ngpsref + ngpseph + nboguss + nairss + ntamdar
 535
 536 !------------------------------------------------------------------------------!
 537 !     B. PRINT OBS USING THE NEW STRUCTURE
 538 !------------------------------------------------------------------------------!
 539
 540       WRITE (0,'(A)')  &
 541 '------------------------------------------------------------------------------'
 542
 543
 544 ! 0.  FORMAT
 545 ! ==========
 546
 547 !     fmt_info = '(A12,1X,A19,1X,A40,1X,I6,3(F12.3,11X),6X,A5)'
 548      fmt_info = '(A12,1X,A19,1X,A40,1X,I6,3(F12.3,11X),6X,A40)'
 549      fmt_srfc = '(F12.3,I4,F7.2,F12.3,I4,F7.3)'
 550      fmt_each = '(3(F12.3,I4,F7.2),11X,3(F12.3,I4,F7.2),11X,3(F12.3,I4,F7.2))'
 551
 552
 553 ! 1. OPEN FILE FOR VALID OBSERVATIONS OUTPUT
 554 ! ==========================================
 555
 556 ! 1.1 Name of output file
 557 !     -------------------
 558
 559       write(filename,'("obs_gts_",a,".",a)') time_analysis, use_for
 560
 561 !      filename = 'obs_gts.'//use_for
 562
 563       WRITE (0,'(7A)') 'Write ',use_for,' GTS observations in file ',&
 564                         TRIM (filename),' (WRFDA ',TRIM(release_version),')'
 565
 566
 567 ! 1.2 OPEN FILE AT VERSION 3.1 FORMAT
 568 !     -------------------------------
 569
 570       INQUIRE ( UNIT = 99, OPENED = connected )
 571
 572       IF ( .NOT. connected ) THEN
 573           OPEN ( UNIT   = 99,           &
 574                  FILE   = filename,     &
 575                  FORM   = 'FORMATTED',  &
 576                  ACCESS = 'SEQUENTIAL', &
 577                  STATUS = 'REPLACE')
 578       ENDIF
 579
 580       REWIND ( UNIT = 99)
 581
 582       if (ntotal == 0) then
 583          WRITE(0,'(A,I6,A)') "Ntotal=",ntotal, &
 584                             " No observations other than SSMI is written out."
 585 !        CLOSE (UNIT = 99)
 586 !        RETURN
 587       endif
 588
 589 ! 1.3 FILE HEADER FOR NEW FORMAT
 590 !     --------------------------
 591
 592 ! 1.3.1 TOTAL NUMBER OF STATIONS CONTAINED IN FILE
 593 !       ------------------------------------------
 594
 595       WRITE  (UNIT = 99, FMT = '((A,I7),A)',ADVANCE='no' )    &
 596       "TOTAL =",nsynops + nmetars + nshipss + &
 597                 nsounds + npilots + naireps + &
 598                 nsatems + nsatobs + ngpspws + &
 599                 nssmt1s + nssmt2s + namdars + &
 600                 ntovss  + nqscats + nprofls + &
 601                 nbuoyss + nothers + ngpsztd + &
 602                 ngpsref + ngpseph + nboguss + nairss + ntamdar,", "
 603
 604 ! 1.3.5 MISSING VALUE FLAG
 605 !       ------------------
 606
 607        WRITE  (UNIT = 99, FMT = '((A,F8.0),A)') "MISS. =",missing_r,","
 608
 609 ! 1.3.2 NUMBER OF STATIONS PER TYPE
 610 !       ---------------------------
 611
 612 !      WRITE  (UNIT = 99, FMT = '(3(6(A,I7,A),/,:))' )    &
 613       WRITE  (UNIT = 99, FMT = '(6(A,I7,A))' )    &
 614       "SYNOP =",nsynops,", ", &
 615       "METAR =",nmetars,", ", &
 616       "SHIP  =",nshipss,", ", &
 617       "BUOY  =",nbuoyss,", ", &
 618       "BOGUS =",nboguss,", ", &
 619       "TEMP  =",nsounds,", ", &
 620       "AMDAR =",namdars,", ", &
 621       "AIREP =",naireps,", ", &
 622       "TAMDAR=",ntamdar,", ", &
 623       "PILOT =",npilots,", ", &
 624       "SATEM =",nsatems,", ", &
 625       "SATOB =",nsatobs,", ", &
 626       "GPSPW =",ngpspws,", ", &
 627       "GPSZD =",ngpsztd,", ", &
 628       "GPSRF =",ngpsref,", ", &
 629       "GPSEP =",ngpseph,", ", &
 630       "SSMT1 =",nssmt1s,", ", &
 631       "SSMT2 =",nssmt2s,", ", &
 632 !      "SSMI  =",nssmis, ", ", &
 633       "TOVS  =",ntovss, ", ", &
 634       "QSCAT =",nqscats,", ", &
 635       "PROFL =",nprofls,", ", &
 636       "AIRSR =",nairss ,", ", &
 637       "OTHER =",nothers,", "
 638
 639
 640 ! 1.3.4 REFERENCE STATE INFO
 641 !       --------------------
 642
 643         WRITE (UNIT = 99, FMT = '(A,F7.2,A,F7.2,4(A,F7.2),A)') &
 644         "PHIC  =", phic,", XLONC =", xlonc,", TRUE1 =", truelat1,&
 645       ", TRUE2 =",truelat2,", XIM11 =",xim11(idd),", XJM11 =",xjm11(idd),","
 646
 647         WRITE (UNIT = 99, FMT = '(2(A,F7.2),4(A,F7.0),A)') &
 648         "base_temp=",  ts0, ", base_lapse=", tlp, &
 649         ", PTOP  =",  ptop,", base_pres=",  ps0, &
 650         ", base_tropo_pres=", pis0, ", base_strat_temp=", tis0, ","
 651
 652 ! 1.3.4 DOMAINS INFO
 653 !       ------------
 654
 655         WRITE (UNIT = 99, FMT = '(5(A,I7),A)' ) &
 656         "IXC   =", ixc,", JXC   =", jxc,", IPROJ =", iproj,&
 657         ", IDD   =", idd,", MAXNES=",maxnes,","
 658         WRITE (UNIT = 99, FMT = '(A,10(:,I7,A))')  &
 659        "NESTIX=",(nestix (i),", ", i = 1, maxnes)
 660         WRITE (UNIT = 99, FMT = '(A,10(:,I7,A))')  &
 661        "NESTJX=",(nestjx (i),", ", i = 1, maxnes)
 662         WRITE (UNIT = 99, FMT = '(A,10(:,I7,A))')  &
 663        "NUMC  =",(numc   (i),", ", i = 1, maxnes)
 664         WRITE (UNIT = 99, FMT = '(A,10(:,F7.2,A))')&
 665        "DIS   =",(dis    (i),", ",i = 1, maxnes)
 666         WRITE (UNIT = 99, FMT = '(A,10(:,I7,A))')  &
 667        "NESTI =",(nesti  (i),", ", i = 1, maxnes)
 668         WRITE (UNIT = 99, FMT = '(A,10(:,I7,A))')  &
 669        "NESTJ =",(nestj  (i),", ", i = 1, maxnes)
 670
 671
 672
 673 ! 1.3.6 VARIABLE NAME AND UNITS
 674 !       -----------------------
 675
 676        WRITE (UNIT = 99, FMT = '(A)' ) &
 677       "INFO  = PLATFORM, DATE, NAME, LEVELS, LATITUDE, LONGITUDE, ELEVATION, ID."
 678        WRITE (UNIT = 99, FMT = '(A)' ) &
 679       "SRFC  = SLP, PW (DATA,QC,ERROR)."
 680        WRITE (UNIT = 99, FMT = '(A)' ) &
 681       "EACH  = PRES, SPEED, DIR, HEIGHT, TEMP, DEW PT, HUMID (DATA,QC,ERROR)*LEVELS."
 682
 683 ! 1.3.7 FORMATS
 684 !       -------
 685
 686       WRITE (UNIT = 99, FMT = '(2A)' ) 'INFO_FMT = ', TRIM (fmt_info)
 687       WRITE (UNIT = 99, FMT = '(2A)' ) 'SRFC_FMT = ', TRIM (fmt_srfc)
 688       WRITE (UNIT = 99, FMT = '(2A)' ) 'EACH_FMT = ', TRIM (fmt_each)
 689
 690
 691 ! 1.3.7 END OF HEADER
 692 !       -------------
 693
 694         WRITE (UNIT = 99, FMT = '(A)') &
 695 "#------------------------------------------------------------------------------#"
 696
 697
 698 ! 2.  WRITE OBSERVATIONS
 699 ! ======================
 700
 701       nmultis  = 0
 702       nsingles = 0
 703       nlevels  = 0
 704       nwrites  = 0
 705
 706
 707 ! 2.1  Loop over stations
 708 !      ------------------
 709
 710 stations: &
 711        DO n = 1, number_of_obs
 712
 713 ! 2.2 Index of obs
 714 !     ------------
 715
 716       loop_index = windex (n)
 717
 718 ! 2.3 Check if station is valid
 719 !     -------------------------
 720
 721 stations_valid: &
 722       IF (obs (loop_index)%info%discard ) THEN
 723
 724           CYCLE  stations
 725
 726       ELSE stations_valid
 727
 728       READ (obs (loop_index) % info % platform (4:6), '(I3)') fm
 729
 730       if ((fm == 125) .OR. (fm == 126))  CYCLE stations
 731
 732 ! SATEM reference pressure is assigned to slp:
 733
 734       if (fm == 86) then
 735            obs (loop_index) % ground % slp % data =  &
 736               obs (loop_index) % ground % ref_pres % data
 737            obs (loop_index) % ground % slp % qc =  &
 738               obs (loop_index) % ground % ref_pres % qc
 739            obs (loop_index) % ground % pw % data =  &
 740               obs (loop_index) % ground % cloud_cvr % data
 741 !           obs (loop_index) % ground % pw % qc =  &
 742 !              obs (loop_index) % ground % cloud_cvr % qc
 743       endif
 744
 745 ! 2.4 Write station info
 746 !     ------------------
 747
 748 ! First, fix up the station id
 749 !
 750       SELECT CASE (fm)
 751       CASE ( 86, 87, 88 )
 752         sid = obs (loop_index) % location % id
 753         kx = index(obs (loop_index) % location % id,'GOES')
 754         if (kx .eq. 0) then
 755           kx = index(obs (loop_index) % location % id,'MET')
 756           if (kx .ne. 0) then
 757             sid = 'MET' // sid(kx+3:40)
 758           else
 759             kx = index(obs (loop_index) % location % id,'MODIS')
 760             if (kx .ne. 0) sid = 'MODIS'
 761           endif
 762         else
 763           sid = 'G-' // sid(kx+5:40)
 764         endif
 765       CASE ( 141 )
 766         sid = 'MODIS'
 767       CASE ( 281 )
 768         sid = 'QSCAT'
 769       CASE ( 13, 18, 19 )
 770          kx = index(obs (loop_index) % location % name,' >>>')
 771          sid = obs (loop_index) % location % name
 772          if (kx .ne. 0) then
 773            sid = sid(kx+5:kx+9)
 774          else
 775            if (fm .eq. 13) then
 776              sid = 'SHIP'
 777            else if (fm .eq. 18) then
 778              sid = 'BUOY'
 779            else if (fm .eq. 19) then
 780              sid = 'C-MAN'
 781            else
 782              sid = obs (loop_index) % location % id
 783            endif
 784          endif
 785       CASE ( 15, 16 )
 786         kx = index(obs (loop_index) % location % name,'ICAO ')
 787         sid = obs (loop_index) % location % name
 788         if (kx .ne. 0) then
 789           sid = sid(kx+5:kx+9)
 790         else
 791           sid = obs (loop_index) % location % id
 792         endif
 793       CASE ( 133 )
 794         sid = 'AIRSRET'
 795       CASE DEFAULT
 796         sid = obs (loop_index) % location % id
 797       END SELECT
 798
 799 ! To keep the levels for use later. Because some of the BUOY data have the empty
 800 ! data record from the NCAR archived LITTLE_R file, the number of levels could
 801 ! become 0, which is not cosistent with the usage in WRFVar (YRG, 02/24/2009):
 802
 803       k_levels = obs (loop_index) % info % levels
 804       if ( (fm == 18 .or. fm == 19 ) .and. &
 805            obs (loop_index) % info   % elevation == 0.0 ) then
 806            if ( obs (loop_index) % info % levels == 0 ) &
 807                 obs (loop_index) % info % levels = 1
 808       endif
 809 ! ..............................................................................
 810
 811       WRITE (UNIT = 99, FMT = TRIM (fmt_info))          &
 812              obs (loop_index) % info     % platform,    &
 813              obs (loop_index) % valid_time % date_mm5,  &
 814              obs (loop_index) % location % name,        &
 815              obs (loop_index) % info % levels,          &
 816              obs (loop_index) % location % latitude,    &
 817              obs (loop_index) % location % longitude,   &
 818              obs (loop_index) % info     % elevation,   &
 819              sid
 820 !            obs (loop_index) % location % id
 821
 822 ! 2.4 Write surface info
 823 !     ------------------
 824
 825          change_qc = .false.
 826
 827          if( domain_check_h .and. &
 828            (obs (loop_index) % location % xjc < 1.0 .or. &
 829             obs (loop_index) % location % xjc > rew_cross .or. &
 830             obs (loop_index) % location % yic < 1.0 .or. &
 831             obs (loop_index) % location % yic > rns_cross) ) then
 832             obs (loop_index) % ground % slp % qc = -88
 833             obs (loop_index) % ground % pw  % qc = -88
 834             change_qc = .true.
 835          end if
 836
 837       WRITE (UNIT = 99, FMT = TRIM (fmt_srfc))          &
 838              obs (loop_index) % ground % slp % data,    &
 839              obs (loop_index) % ground % slp % qc,      &
 840              obs (loop_index) % ground % slp % error,   &
 841              obs (loop_index) % ground % pw  % data,    &
 842              obs (loop_index) % ground % pw  % qc,      &
 843              obs (loop_index) % ground % pw  % error
 844
 845              nwrites = nwrites + 1
 846
 847 ! 2.6 Initialise pointer
 848 !     ------------------
 849
 850       current => obs (loop_index) % surface
 851
 852 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 853 ! In the NCAR archie LITTLE_R files (created by Jim Bresch), there is another
 854 ! problem,
 855 !
 856 ! The data record for some of BUOY (FM-18, 19) were filled with all parameters
 857 ! (p, h, t, td,...) missing even though SLP and elevation=0.0m in the header
 858 ! record are available. This caused the BUOY report has only 0-level data,
 859 ! which will be discarded by WRFDA/var/da/da_obs_io/da_read_obs_ascii.inc.
 860 !
 861 ! In tropical storm cases, however, assimilation of the SLP obs from BUOY could
 862 ! be important. This BUOY SLP obs should not be discarded.
 863 !
 864 ! The following modifications attempt to fix these problems.
 865 !
 866 ! In CWB, AFWA, and KMA LITTLE_R files, we did not find the above problems, and
 867 ! the following modifications may not needed.
 868 !
 869 !                                        Yong-Run Guo, 02/19/2009
 870 ! .................................................................................
 871
 872       if ( (fm == 18 .or. fm == 19 ) .and. &
 873            obs (loop_index) % info   % elevation == 0.0 ) then
 874
 875            pressure_error = 100. ! for BUOY
 876
 877            if ( k_levels == 0 ) then
 878
 879              WRITE (UNIT = 99, FMT = TRIM (fmt_each))    &
 880                obs (loop_index) % ground % slp % data,    &
 881                obs (loop_index) % ground % slp % qc,      &
 882                pressure_error, &
 883                missing_r, -88, 1.40,  &  ! speed
 884                missing_r, -88, 5.00,  &  ! direction
 885                obs (loop_index) % info   % elevation, 0, 6.0, &
 886                missing_r, -88,  2.0,  &  ! temperature
 887                missing_r, -88,  2.0,  &  ! dew-point
 888                missing_r, -88, 10.0      ! RH
 889                CYCLE  stations
 890            endif
 891       endif
 892 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 893
 894 ! 2.7 Loop over levels
 895 !     ----------------
 896
 897       is_sound = -1
 898
 899 levels:&
 900       DO WHILE (ASSOCIATED (current))
 901
 902       nlevels  = nlevels  + 1
 903       is_sound = is_sound + 1
 904
 905 ! SATEM thickness is assigned to temperature field:
 906
 907       if (fm == 86) &
 908          current % meas % temperature = current % meas % thickness
 909
 910 ! 2.9 Write height, pressure, temp, mixing ratio, wind and model vertical coord
 911 !     -------------------------------------------------------------------------
 912
 913       if(change_qc) then
 914          current % meas % temperature % qc = -88
 915          current % meas % dew_point   % qc = -88
 916          current % meas % rh          % qc = -88
 917
 918          if(current % meas % height   % qc >= 0) then
 919             current % meas % pressure % qc = -88
 920          end if
 921       end if
 922
 923 ! Because there is no "qc" in pressure and height in data structure in WRFVar,
 924 ! all pressure and height with bad flags(<0), the data value are set to be missing
 925 ! (YRG, 02/11/2009)
 926
 927       if (current % meas % pressure % qc < 0) &
 928           current % meas % pressure % data = missing_r
 929       if (current % meas % height   % qc < 0) &
 930           current % meas % height   % data = missing_r
 931 ! ...............................................................
 932
 933           if (fm == 118 .or. fm == 116) then
 934
 935       WRITE (UNIT = 99, FMT = TRIM (fmt_each))    &
 936              current % meas % pressure    % data, &    ! pressure
 937              current % meas % pressure    % qc,   &
 938              current % meas % pressure    % error,&
 939              current % meas % u       % data, &        ! latitude
 940              current % meas % u       % qc,   &
 941              current % meas % u       % error,&
 942              current % meas % v   % data, &            ! longitude
 943              current % meas % v   % qc,   &
 944              current % meas % v   % error,&
 945              current % meas % height      % data, &    ! height
 946              current % meas % height      % qc,   &
 947              current % meas % height      % error,&
 948              current % meas % temperature % data, &    ! temperature
 949              current % meas % temperature % qc,   &
 950              current % meas % temperature % error,&
 951              current % meas % dew_point   % data, &    ! refractivity
 952              current % meas % dew_point   % qc,   &
 953              current % meas % dew_point   % error,&
 954              current % meas % direction   % data, &    ! azimuth
 955              current % meas % direction   % qc,   &
 956              current % meas % direction   % error,&
 957              current % meas % speed       % data, &    ! impact parameter /1000.
 958              current % meas % speed       % qc,   &
 959              current % meas % speed       % error,&
 960              current % meas % rh          % data, &    ! bending angle*1.e7
 961              current % meas % rh          % qc,   &
 962              current % meas % rh          % error
 963           else
 964
 965       WRITE (UNIT = 99, FMT = TRIM (fmt_each))    &
 966              current % meas % pressure    % data, &
 967              current % meas % pressure    % qc,   &
 968              current % meas % pressure    % error,&
 969              current % meas % speed       % data, &
 970              current % meas % speed       % qc,   &
 971              current % meas % speed       % error,&
 972              current % meas % direction   % data, &
 973              current % meas % direction   % qc,   &
 974              current % meas % direction   % error,&
 975              current % meas % height      % data, &
 976              current % meas % height      % qc,   &
 977              current % meas % height      % error,&
 978              current % meas % temperature % data, &
 979              current % meas % temperature % qc,   &
 980              current % meas % temperature % error,&
 981              current % meas % dew_point   % data, &
 982              current % meas % dew_point   % qc,   &
 983              current % meas % dew_point   % error,&
 984              current % meas % rh          % data, &
 985              current % meas % rh          % qc,   &
 986              current % meas % rh          % error
 987           endif
 988
 989              nwrites = nwrites + 1
 990
 991
 992 ! 3.  GO TO NEXT OBS
 993 !     ==============
 994
 995 ! 3.1 Go to next level
 996 !     -----------------
 997
 998       current => current%next
 999
1000       ENDDO levels
1001
1002 ! 3.2 Count surface and sounding
1003 !     --------------------------
1004
1005       if (is_sound .gt. 0) then
1006           nmultis = nmultis + 1
1007       else
1008           nsingles  = nsingles + 1
1009       endif
1010
1011 ! 3.3 Go to next valid station
1012 !     -----------------
1013
1014       ENDIF stations_valid
1015
1016 ! 3.3 Go to next record
1017 !     -----------------
1018
1019       ENDDO stations
1020
1021
1022 ! 4.   CLOSE OUTPUT FILES
1023 ! ========================
1024
1025       CLOSE (UNIT = 99)
1026
1027
1028 ! 5.  PRINT DIAGNOSTIC
1029 ! =====================
1030
1031       WRITE (0, '(/,A,I8,A,A)') &
1032      'Wrote ',nwrites,' lines of data in file: ',TRIM (filename)
1033       WRITE (0, '(A)') ' '
1034
1035
1036    RETURN
1037
1038 END SUBROUTINE output_gts_31
1039
1040 !----------------------------------------------------------------------
1041 SUBROUTINE output_prep (max_number_of_obs, obs, number_of_obs, windex,&
1042                          prepbufr_table_filename, &
1043                          prepbufr_output_filename, &
1044                          nsynops, nshipss, nmetars,                    &
1045                          npilots, nsounds, nsatems,                    &
1046                          nsatobs, naireps, ngpspws, ngpsztd, ngpsref,  &
1047                          ngpseph, nssmt1s, nssmt2s, nssmis,  ntovss,   &
1048                          nothers, namdars, nqscats, nprofls,nbuoyss,   &
1049                          nboguss, missing_r, time_analysis)
1050
1051 !----------------------------------------------------------------------
1052 ! Write observations in PREPBUFR
1053 ! Jim Bresch, February 2006
1054 ! J. F. Drake, November 2005
1055 !----------------------------------------------------------------------
1056
1057
1058   IMPLICIT NONE
1059
1060
1061   INTEGER,                                      INTENT (in) :: max_number_of_obs
1062   TYPE (report), DIMENSION (max_number_of_obs), INTENT (inout) :: obs
1063   INTEGER,                                      INTENT (in) :: number_of_obs
1064   INTEGER,       DIMENSION (max_number_of_obs), INTENT (in) :: windex
1065   REAL,                                         INTENT (in) :: missing_r
1066   CHARACTER (LEN =  19),                        INTENT (in) :: time_analysis
1067   CHARACTER (LEN =  80),                        INTENT(in)  :: prepbufr_table_filename, &
1068                                                                prepbufr_output_filename
1069   INTEGER,                                      INTENT (in) :: nsynops,nmetars,&
1070                                                                nshipss,nsounds,&
1071                                                                npilots,naireps,&
1072                                                                nsatems,nsatobs,&
1073                                                                ngpspws,nssmt1s,&
1074                                                                nssmt2s, nssmis,&
1075                                                                ntovss, namdars,&
1076                                                                nqscats,nothers,&
1077                                                                nprofls,nbuoyss,&
1078                                                                ngpsztd,ngpsref,&
1079                                                                ngpseph,nboguss
1080
1081   TYPE (measurement ) , POINTER :: current
1082   INTEGER                       :: bfout, bftable, loop_index
1083   INTEGER                       :: i, n, nlv, nmax, ntotal
1084   INTEGER                       :: nmultis, nsingles, nlevels, nwrites
1085   INTEGER                       :: is_sound, fm, idate
1086   INTEGER                       :: year, month, day, hour, minute, second
1087   INTEGER                       :: mxmn, kx
1088   LOGICAL                       :: connected, vld
1089   CHARACTER (LEN = 6)           :: cfm(300)
1090   REAL(kind=8), allocatable, dimension(:,:) :: r8arr
1091   REAL                          :: mixing_ratio, temperat
1092   REAL                          :: press
1093   REAL(kind=8)                  :: wqm, virtual_temperature, &
1094                                    specific_humidity
1095
1096   CHARACTER (LEN=40) :: SID
1097   REAL, parameter  :: PI=3.1415926535
1098
1099 !---------------------------------------------------------------------
1100 !  REAL                          :: rew_cross, rns_cross
1101
1102 #ifdef BUFR
1103
1104 ! rew_cross=real(nestjx(idd)-1)
1105 ! rns_cross=real(nestix(idd)-1)
1106 !---------------------------------------------------------------------
1107
1108 ! Set up missing value for BUFRLIB
1109
1110       bfout = 11
1111       bftable = 10
1112
1113 ! Initialilze cfm
1114 !---------------------------------------------------------------------
1115
1116       cfm = 'UNKNOW'
1117       cfm(12)    = 'ADPSFC'
1118       cfm(13)    = 'SFCSHP'
1119       cfm(14:16) = 'ADPSFC'
1120       cfm(18:19) = 'SFCSHP'
1121       cfm(32:38) = 'ADPUPA'
1122       cfm(42)    = 'AIRCFT'
1123       cfm(86:87) = 'GOESND'
1124       cfm(88)    = 'SATWND'
1125       cfm(96:97) = 'AIRCFT'
1126       cfm(111)   = 'GPSIPW'
1127       cfm(116)   = 'GPSREF'
1128       cfm(125)   = 'SPSSMI'
1129       cfm(132)   = 'PROFLR'
1130       cfm(135)   = 'SYNDAT'
1131       cfm(141)   = 'SATWND'     ! modis used to have fm code 141, but now uses 88
1132       cfm(281)   = 'QKSWND'
1133
1134
1135       ntotal =   nsynops + nmetars + nshipss + &
1136                  nsounds + npilots + naireps + &
1137                  nsatems + nsatobs + ngpspws + &
1138                  nssmis  + &
1139                  nssmt1s + nssmt2s + ntovss  + &
1140                  namdars + nqscats + nprofls + &
1141                  nbuoyss + nothers + ngpsztd + &
1142                  ngpsref + nboguss
1143
1144       write(0,*) 'ntotal = ',ntotal
1145       if (ntotal == 0) then
1146          WRITE(0,'(A,I6,A)') "Ntotal=",ntotal, &
1147                             " No prepbufr observations written out."
1148          RETURN
1149       endif
1150
1151 !     Determine the required size of the real*8 buffer and
1152 !     allocate it.
1153
1154       nmax = max(nsynops + nmetars, nshipss + nbuoyss, &
1155                  nsounds, npilots, naireps + namdars, nprofls)
1156
1157       mxmn = 58
1158       ALLOCATE (r8arr(mxmn, 255))
1159
1160 ! 1. Open unit bftable for the PREPBUFR table
1161 !---------------------------------------------------------------------
1162
1163       INQUIRE (UNIT=bftable, OPENED= connected )
1164       IF (connected) THEN
1165          WRITE(0,'(A, i4, A)') "A file is already connected to unit", bftable, &
1166                         "which is needed for the PREPBUFR tables"
1167          STOP
1168       ENDIF
1169
1170       OPEN (UNIT = bftable, FILE = prepbufr_table_filename)
1171
1172 ! 2. Open unit bfout FOR the PREPBUFR output file
1173 !---------------------------------------------------------------------
1174
1175       OPEN (UNIT = bfout, FILE = prepbufr_output_filename, &
1176             FORM = "UNFORMATTED")
1177
1178       CALL OPENBF (bfout, 'OUT', bftable)
1179
1180       WRITE (0,'(2A)') 'Write 3DVAR PREPBUFR observations in file ',&
1181                         TRIM (prepbufr_output_filename)
1182
1183 ! 3.  WRITE OBSERVATIONS
1184 !---------------------------------------------------------------------
1185
1186       nmultis  = 0
1187       nsingles = 0
1188       nlevels  = 0
1189       nwrites  = 0
1190
1191
1192 ! 3.1  Loop over stations
1193 !      ------------------
1194    print *, 'number of obs=', number_of_obs
1195 stations: &
1196    DO n = 1, number_of_obs
1197
1198 ! 3.2 Index of obs
1199 !     ------------
1200
1201    loop_index = windex (n)
1202
1203 ! 3.3 Check if station is valid
1204 !     -------------------------
1205
1206 stations_valid: &
1207    IF (obs (loop_index)%info%discard ) THEN
1208       CYCLE stations
1209    ELSE stations_valid
1210       READ (obs (loop_index) % info % platform (4:6), '(I3)') fm
1211       IF (fm == 126 .OR. cfm(fm) == 'UNKNOW') &
1212           CYCLE stations
1213
1214 ! SATEM reference pressure is assigned to slp:
1215
1216       IF (fm == 86) then
1217            obs (loop_index) % ground % slp % data =  &
1218               obs (loop_index) % ground % ref_pres % data
1219            obs (loop_index) % ground % slp % qc =  &
1220               obs (loop_index) % ground % ref_pres % qc
1221            obs (loop_index) % ground % pw % data =  &
1222               obs (loop_index) % ground % cloud_cvr % data
1223 !           obs (loop_index) % ground % pw % qc =  &
1224 !              obs (loop_index) % ground % cloud_cvr % qc
1225       endif
1226
1227 SUBSET: &
1228        SELECT CASE (cfm(fm))
1229        CASE ('ADPSFC', 'SFCSHP', 'AIRCFT', 'ADPUPA', 'AIRCAR', 'PROFLR', 'SATWND', &
1230              'SYNDAT', 'GPSIPW', 'GPSREF', 'QKSWND', 'SPSSMI')
1231
1232 ! 3.4 Possibly open a new BUFR message
1233 !     --------------------------------------
1234          CALL split_date_char ( obs (loop_index) % valid_time % date_mm5, &
1235                                 year , month , day , hour , minute , second )
1236
1237          idate = year * 1000000 + month * 10000 + day * 100 + hour
1238
1239 !        write(6,*) 'working on ',cfm(fm),' id = ',obs (loop_index) % location % id
1240
1241          CALL OPENMB (bfout, cfm(fm), idate)
1242
1243 !        Fill r8arr with HEADR Table-B mnemonics
1244 !        SID
1245          SELECT CASE (cfm(fm))
1246          CASE ('SATWND', 'SATEMP','GOESND')
1247            kx = index(obs (loop_index) % location % id,'GOES')
1248            if (kx .eq. 0) kx = index(obs (loop_index) % location % id,'MET')
1249            sid = obs (loop_index) % location % id
1250            sid = sid(kx:40)
1251          CASE ('QKSWND')
1252            sid = 'QUIKSCAT'
1253          CASE ('SPSSMI')
1254            sid = 'SSMI'
1255          CASE ('SFCSHP')
1256            kx = index(obs (loop_index) % location % name,' >>>')
1257            sid = obs (loop_index) % location % name
1258            if (kx .ne. 0) then
1259              sid = sid(kx+5:kx+9)
1260            else
1261              if (fm .eq. 13) then
1262                sid = 'SHIP'
1263              else if (fm .eq. 18) then
1264                sid = 'BUOY'
1265              else if (fm .eq. 19) then
1266                sid = 'C-MAN'
1267              else
1268                sid = obs (loop_index) % location % id
1269              endif
1270            endif
1271          CASE ('ADPSFC')
1272            if (fm .eq. 15 .or. fm .eq. 16) then
1273              kx = index(obs (loop_index) % location % name,'ICAO ')
1274              sid = obs (loop_index) % location % name
1275              if (kx .ne. 0) then
1276                sid = sid(kx+5:kx+9)
1277              else
1278                sid = obs (loop_index) % location % id
1279              endif
1280            else
1281              sid = obs (loop_index) % location % id
1282            endif
1283          CASE DEFAULT
1284            sid = obs (loop_index) % location % id
1285          END SELECT
1286          r8arr(1,1) = stuff (sid)
1287 !        XOB
1288          r8arr(2,1) = obs (loop_index) % location % longitude
1289 !        YOB
1290          IF (r8arr(2,1) < 0.) r8arr(2,1) = 360. + r8arr(2,1)
1291          r8arr(3,1) = obs (loop_index) % location % latitude
1292
1293 !        observation time - cycle time, DHR
1294 !        Assume that time_analysis is the same as cycle time.
1295
1296          r8arr(4,1) = 0
1297
1298 !        ELV
1299 !        r8arr(5,1) = obs (loop_index) % info % elevation
1300          call assignv( obs (loop_index) % info % elevation, r8arr(5,1) )
1301
1302 !        report sequence number SQN
1303          if ( obs (loop_index) % info % seq_num == missing ) then
1304             r8arr(10,1) = bufrlib_missing
1305          else
1306             r8arr(10,1) = obs (loop_index) % info % seq_num
1307          end if
1308
1309 !        reported observation time RPT
1310          r8arr(11,1) = hour + (minute + second / 60.) / 60.
1311
1312 !        indicator whether obs time in dhr was corrected, TCOR
1313          r8arr(12,1) = 0
1314
1315 !        instrument type ITP
1316          r8arr(9,1) = bufrlib_missing
1317
1318 !        input report type T29
1319          r8arr(7,1) = t29(fm, obs (loop_index) % location % id)
1320
1321 !        report subtype TSB
1322 !           no distinction in little_r
1323          r8arr(8,1) = bufrlib_missing
1324
1325          current => obs (loop_index) % surface
1326          IF (ASSOCIATED (current)) then
1327            vld = .true.
1328          else
1329            vld = .false.
1330          endif
1331          press = missing_r
1332          temperat = missing_r
1333          mixing_ratio = missing_r
1334          if (vld) then
1335            press = obs (loop_index) % surface % meas % pressure % data
1336            temperat = obs (loop_index) % surface % meas % temperature % data
1337            mixing_ratio = obs (loop_index) % surface % meas % qv % data
1338          endif
1339 !        write(0,*) 'temperat = ',temperat,' loop_index = ',loop_index
1340
1341          IF ( .not. eps_equal(mixing_ratio, missing_r, 1.) ) THEN       ! Is Humidity available?
1342            specific_humidity = mixing_ratio / (1 + mixing_ratio)
1343            IF ( .not. eps_equal(temperat, missing_r, 1.) ) THEN        ! Is temperature available?
1344              virtual_temperature = temperat * (1 + mixing_ratio/eps) / &    ! TOB is virtual temperature.
1345                                                   (1 + mixing_ratio)
1346              virtual_temperature = virtual_temperature - celkel       ! convert to Celsius
1347              r8arr(6,1) = typ (press, cfm(fm), fm, .true., missing_r)
1348              specific_humidity = 1e6 * specific_humidity              ! convert humidity to mg/kg
1349            else
1350              virtual_temperature = bufrlib_missing
1351              r8arr(6,1) = typ (press, cfm(fm), fm, .false., missing_r)
1352            ENDIF
1353          ELSE                                    ! No mixing_ratio, so store temperature instead of Tv
1354            specific_humidity = bufrlib_missing
1355            if ( .not. eps_equal(temperat, missing_r, 1.) ) then
1356              virtual_temperature = temperat - celkel
1357            else
1358              virtual_temperature = bufrlib_missing
1359            endif
1360            r8arr(6,1) = typ (press, cfm(fm), fm, .false., missing_r)
1361          ENDIF
1362
1363 !        Write out HEADR
1364          CALL UFBSEQ (bfout, r8arr, 12, 1, nlv, 'HEADR')
1365
1366 !        PREPBUFR data category, CAT
1367 !        http://www.emc.ncep.noaa.gov/mmb/data_processing/prepbufr.doc/table_1.htm#c
1368          r8arr(1,1) = 6               ! default value
1369          SELECT CASE (cfm(fm))
1370          CASE ('ADPSFC', 'SFCSHP', 'AIRCAR', 'AIRCFT', 'SATWND', 'QKSWND')
1371             if ( .not. eps_equal(temperat, missing_r, 1.)) then
1372               if ( eps_equal(obs (loop_index) % surface % meas % u % data, missing_r, 1.) ) &
1373                r8arr(1,1) = 0  ! surface temperature only (therefore, a mass report)
1374             endif
1375 !           write(0,*) 'cat is set to ',r8arr(1,1)
1376 !           write(0,*) 'temperat = ',temperat,' missing_r = ',missing_r
1377          END SELECT
1378
1379          CALL UFBINT (bfout, r8arr, mxmn, 1, nlv, 'CAT')
1380
1381          if ( cfm(fm) .ne. 'GPSIPW') then
1382
1383 !        PINFO :  PEVN
1384          call assignv(press, r8arr(1,1))
1385 !              Convert to mb if not missing.
1386          IF ( r8arr(1,1) /= bufrlib_missing ) r8arr(1,1) = r8arr(1,1) * .01    !  POB
1387
1388          r8arr(2,1) = bufrlib_missing
1389          if (vld) r8arr(2,1) = qz (obs (loop_index) % surface % meas % pressure % qc)   !  PQM
1390          r8arr(3,1) = 1                                  !  pressure program code PPC
1391          r8arr(4,1) = 100                                ! pressure reason code PRC
1392          CALL UFBINT (bfout, r8arr, mxmn, 1, nlv, 'POB PQM PPC PRC')
1393
1394 !           PBACKG : POE and PFC
1395 !        Convert error from Pa to mb
1396          r8arr(1,1) = bufrlib_missing
1397          if (vld) r8arr(1,1) = obs (loop_index) % surface % meas % pressure % error * &   ! POE
1398                       0.01
1399          r8arr(2,1) = bufrlib_missing                                            ! PFC
1400 !           PPOSTP :  pressure analyzed value PAN, PCL, PCS
1401          r8arr(3:5,1) = bufrlib_missing
1402          CALL UFBINT (bfout, r8arr, mxmn, 1, nlv, &
1403                    'POE PFC PAN PCL PCS')
1404
1405 !        QINFO  :    QEVN and TDO
1406          r8arr(1,1) = specific_humidity            ! QOB
1407          r8arr(2,1) = bufrlib_missing
1408          if (vld) r8arr(2,1) = qz (obs (loop_index) % surface % meas % qv % qc)     ! QCM
1409          r8arr(3,1) = 1                            ! humidity program code QPC
1410          r8arr(4,1) = 100                          ! humidity reason code QRC
1411          r8arr(5,1) = bufrlib_missing              ! TDO
1412
1413 !           QBACKG : QOE and QFC
1414 !  NCEP uses 2.0, while WRF-Var/obsproc uses 1.0
1415          ! r8arr(6,1) = 2.0
1416          r8arr(6,1) = 1.0
1417          r8arr(7:8,1) = bufrlib_missing
1418
1419 !           QPOSTP
1420          r8arr(9:11,1) = bufrlib_missing     ! specific humidity analyzed value QAN, PCL, PCS
1421          CALL UFBINT (bfout, r8arr, mxmn, 1, nlv, &
1422                       'QOB QQM QPC QRC TDO QOE QFC QAN PCL PCS')
1423
1424 !        TINFO : TEVN and TVO
1425          r8arr(1,1) = virtual_temperature       ! TOB
1426          r8arr(2,1) = bufrlib_missing
1427          if (vld) r8arr(2,1) = qz (obs (loop_index) % surface % meas % temperature % qc)    ! TQM
1428          r8arr(3,1) = 1                         ! temperature program code TPC
1429          r8arr(4,1) = 100                       ! temperature reason code TRC
1430          r8arr(5,1) = bufrlib_missing           ! TVO
1431
1432 !           TBACKG :  TOE and TFC
1433          r8arr(6,1) = bufrlib_missing
1434          if (vld) r8arr(6,1) = obs (loop_index) % surface % meas % temperature % error
1435          r8arr(7,1) = bufrlib_missing
1436
1437 !           QPOSTP
1438          r8arr(9:10,1) = bufrlib_missing        ! temperature analyzed value TAN, PCL, PCS
1439          CALL UFBINT (bfout, r8arr, mxmn, 1, nlv, &
1440                       'TOB TQM TPC TRC TVO TOE TFC TAN PCL PCS')
1441
1442 !        ZINFO :  ZEVN
1443          r8arr(1:2,1) = bufrlib_missing
1444          if (vld) call assignv( obs (loop_index) % surface % meas % height % data, &
1445             r8arr(1,1) )                        ! ZOB
1446          if (vld) then
1447          if ( cfm(fm) .eq. 'SFCSHP' .and. &
1448               eps_equal( obs (loop_index) % surface % meas % height % data, 0., 0.5) ) then
1449            r8arr(2,1) = 0.          ! If a ship's elevation is zero, fix the zqm
1450          else
1451            r8arr(2,1) = qz (obs (loop_index) % surface % meas % height % qc)     ! ZQM
1452          endif
1453          endif
1454          r8arr(3,1) = 1                         ! program code ZPC
1455          r8arr(4,1) = 100                       ! reason code ZRC
1456
1457 !           ZBACKG :  ZOE and ZFC
1458          r8arr(5:6,1) = bufrlib_missing
1459          if (vld) r8arr(5,1) = obs (loop_index) % surface % meas % height % error * 0.1
1460
1461 !           ZPOSTP
1462          r8arr(7:9,1) = bufrlib_missing         ! analyzed value ZAN, PCL, PCS
1463          CALL UFBINT (bfout, r8arr, mxmn, 1, nlv, &
1464                       'ZOB ZQM ZPC ZRC ZOE ZFC ZAN PCL PCS')
1465
1466 !        WINFO : WEVN
1467          if (vld) then
1468          if ( eps_equal( obs (loop_index) % surface % meas % speed % data, missing_r, 1.) .or.  &
1469               eps_equal( obs (loop_index) % surface % meas % direction % data, missing_r, 1.) ) then
1470            r8arr(1,1) = bufrlib_missing   ! UOB
1471            r8arr(5,1) = bufrlib_missing   ! VOB
1472          else
1473            r8arr(1,1) = obs (loop_index) % surface % meas % speed % data * &  ! UOB
1474                               (-sin(pi/180.* obs (loop_index) % surface % meas % direction % data))
1475            r8arr(5,1) = obs (loop_index) % surface % meas % speed % data * &  ! VOB
1476                               (-cos(pi/180.* obs (loop_index) % surface % meas % direction % data))
1477          endif
1478
1479 !        Combination of direction%qc and speed%qc give WQM
1480          r8arr(2,1) = qz (MAX(obs (loop_index) % surface % meas % speed % qc, &    !  WQM
1481                           obs (loop_index) % surface % meas % direction % qc))
1482          wqm = r8arr(2,1)
1483          r8arr(3,1) = 1                         ! wind program code WPC
1484          r8arr(4,1) = 100                       ! wind reason code WRC
1485
1486 !           DFEVN Wind (direction/speed)
1487          call assignv( obs (loop_index) % surface % meas % direction % data, &
1488             r8arr(6,1) )                        ! DDO
1489
1490          call assignv( obs (loop_index) % surface % meas % speed % data, &
1491             r8arr(7,1) )                        ! FFO
1492             r8arr(7,1) = r8arr(7,1) * 1.9425    ! convert to knots
1493
1494          r8arr(8,1) = qz (MAX(obs (loop_index) % surface % meas % direction % qc, &
1495                           obs (loop_index) % surface % meas % speed % qc))        ! DFQ
1496          r8arr(9,1) = 1                         ! wind program code DFP
1497          r8arr(10,1) = 100                      ! wind reason code DFR
1498          endif
1499          if ( cfm(fm) .ne. 'GPSREF' ) then
1500            CALL UFBINT (bfout, r8arr, mxmn, 1, nlv, &
1501                         'UOB WQM WPC WRC VOB DDO FFO DFQ DFP DFR')
1502          endif
1503
1504          SELECT CASE (cfm(fm))
1505          CASE ('ADPSFC', 'ADPUPA', 'AIRCAR', 'AIRCFT', 'SFCSHP', 'SATWND')
1506 !           WBACKG
1507             if (vld) then
1508             r8arr(1,1) = obs (loop_index) % surface % meas % speed % error     ! WOE
1509             r8arr(2:3,1) = bufrlib_missing      ! UFC, VFC
1510             r8arr(4:5,1) = bufrlib_missing      ! WFC_MSQ = UFC_MOD VFC_MOD
1511 !           WPOSTP
1512             r8arr(6:9,1) = bufrlib_missing      ! analyzed values UAN and VAN, PCL, PCS
1513             CALL UFBINT (bfout, r8arr, mxmn, 1, nlv, &
1514                          'WOE UFC VFC UFC_MOD VFC_MOD UAN VAN PCL PCS')
1515             endif
1516          END SELECT
1517
1518          SELECT CASE (cfm(fm))
1519          CASE ('ADPSFC', 'SFCSHP')
1520 !     PMSL
1521 !           write(6,*) 'sid = ',sid,' slp = ',obs (loop_index) % ground % slp % data
1522             if ( eps_equal( obs (loop_index) % ground % slp % data, &
1523                                   missing_r, 1.) ) then
1524               r8arr(1:2,1) = bufrlib_missing      ! PMO PMQ
1525             else
1526               r8arr(1,1) = obs (loop_index) % ground % slp % data * 0.01     ! PMO
1527               r8arr(2,1) = qz (obs (loop_index) % ground % slp % qc)         ! PMQ
1528             endif
1529 !           write(6,*) 'storing PMO = ',r8arr(1,1)
1530             CALL UFBINT (bfout, r8arr, 2, 1, nlv, 'PMO PMQ')
1531
1532             IF (cfm(fm) == 'ADPSFC') THEN
1533                r8arr(1,1) = bufrlib_missing     !     ALTMSQ (ALSE)
1534                CALL UFBINT (bfout, r8arr, 1, 1, nlv, 'ALSE')
1535 !     WSPDSQ
1536                if (vld) then
1537                IF (obs (loop_index) % surface % meas % direction % data ==  &
1538                    missing_r) THEN
1539                   IF ( eps_equal( obs (loop_index) % surface % meas % speed % data, &
1540                                   missing_r, 1.) ) then
1541                      r8arr(1,1) =  bufrlib_missing
1542                   ELSE
1543                      r8arr(1,1) = obs (loop_index) % surface % meas % speed % data  ! SOB
1544                   ENDIF
1545                ENDIF
1546                r8arr(2,1) = wqm                                      ! SQM
1547                CALL UFBINT (bfout, r8arr, 2, 1, nlv, 'SOB SQM')
1548                endif
1549             ENDIF
1550         END SELECT
1551
1552         else
1553
1554         SELECT CASE (cfm(fm))
1555         CASE ('GPSIPW')
1556           r8arr(1,1) = 6.
1557           r8arr(2,1) = obs (loop_index) % ground  % pw   % data
1558           r8arr(3,1) = qz (obs (loop_index) % ground  % pw   % qc)
1559           r8arr(4,1) = 1
1560           r8arr(5,1) = 100
1561           r8arr(6,1) = obs (loop_index) % ground  % pw   % error
1562           CALL UFBINT (bfout, r8arr, 6, 1, nlv, 'CAT PWO PWQ PWP PWR PWE')
1563         END SELECT
1564
1565         endif    ! gpsipw if-test
1566
1567         if (cfm(fm) .eq. 'SPSSMI') then
1568           r8arr(1:11,1) = bufrlib_missing
1569           r8arr(1,1) = 6.
1570           r8arr(2,1) = 10.0*obs (loop_index) % ground  % pw  % data  !cm to mm
1571           r8arr(3,1) = qz (obs (loop_index) % ground  % pw   % qc)
1572           r8arr(4,1) = 1
1573           r8arr(5,1) = 100
1574           r8arr(6,1) = 10.0*obs (loop_index) % ground  % pw  % error !cm to mm
1575           if (vld) then
1576             call assignv( obs (loop_index) % surface % meas % speed % data, &
1577               r8arr(7,1) )                        ! FFO
1578             r8arr(7,1) = r8arr(7,1) * 1.9425    ! convert to knots
1579
1580             r8arr(8,1) = qz (MAX(obs (loop_index) % surface % meas % direction % qc, &
1581                           obs (loop_index) % surface % meas % speed % qc))        ! DFQ
1582            r8arr(9,1) = 1                         ! wind program code DFP
1583            r8arr(10,1) = 100                      ! wind reason code DFR
1584            r8arr(11,1) = obs (loop_index) % surface % meas % speed % error     ! WOE
1585           endif
1586           CALL UFBINT (bfout, r8arr, 11, 1, nlv, 'CAT PWO PWQ PWP PWR PWE FFO DFQ DFP DFR WOE')
1587         endif
1588
1589 UPA: &
1590       SELECT CASE (cfm(fm))
1591       CASE ('ADPUPA', 'PROFLR', 'SYNDAT', 'GPSREF')
1592 !
1593       current => obs (loop_index) % surface
1594 !
1595 !     Split the processing up into 2 do-loops because prebufr has a character string
1596 !     limit of 80 characters in the ufbint call
1597
1598       is_sound = -1
1599
1600       nlevels = 0
1601
1602 levels1:&
1603       DO WHILE (ASSOCIATED (current))
1604
1605       nlevels  = nlevels  + 1
1606       is_sound = is_sound + 1
1607
1608 ! SATEM thickness is assigned to temperature field:
1609
1610 !     if (fm == 86) &
1611 !        current % meas % temperature = current % meas % thickness
1612
1613
1614 !     if(change_qc) then
1615 !        current % meas % temperature % qc = -88
1616 !        current % meas % dew_point   % qc = -88
1617 !        current % meas % rh          % qc = -88
1618 !
1619 !        if(current % meas % height   % qc >= 0) then
1620 !           current % meas % pressure % qc = -88
1621 !        end if
1622 !     end if
1623
1624       r8arr(1,nlevels) = 1           ! assign cat 1 (mandatory level) for all levels
1625       call assignv(current % meas % pressure % data, r8arr(2,nlevels))
1626       if ( r8arr(2,nlevels) /= bufrlib_missing ) r8arr(2,nlevels) = r8arr(2,nlevels) * .01
1627       r8arr(3,nlevels) = qz (current % meas % pressure % qc)
1628       r8arr(4,nlevels) = 1
1629       r8arr(5,nlevels) = 100
1630       r8arr(6,nlevels) = current % meas % pressure % error * 0.01
1631       press = current % meas % pressure % data
1632       temperat = current % meas % temperature % data
1633       mixing_ratio = current % meas % qv % data
1634       if ( .not. eps_equal(mixing_ratio, missing_r, 1.) ) then
1635         specific_humidity = mixing_ratio / (1 + mixing_ratio)
1636         if ( .not. eps_equal(temperat, missing_r, 1.) ) then        !  Is temperature available?
1637           virtual_temperature = temperat * (1 + mixing_ratio/eps) / &     !  TOB is virtual temperature.
1638                                            (1 + mixing_ratio)
1639           virtual_temperature = virtual_temperature - celkel      !  convert to Celsius
1640           r8arr(12,nlevels) = virtual_temperature
1641           specific_humidity = 1e6 * specific_humidity             !  convert humidity to mg/kg
1642         else
1643           virtual_temperature = bufrlib_missing
1644         endif
1645       else                 ! mixing ratio missing
1646         specific_humidity = bufrlib_missing
1647         if ( .not. eps_equal(temperat, missing_r, 1.) ) then
1648           virtual_temperature = temperat - celkel
1649         else
1650           virtual_temperature = bufrlib_missing
1651         endif
1652       endif
1653
1654       r8arr(7,nlevels) = specific_humidity
1655       r8arr(12,nlevels) = virtual_temperature
1656       r8arr(8,nlevels) = qz (current % meas % dew_point % qc)
1657       r8arr(9,nlevels) = 1
1658       r8arr(10,nlevels) = 100
1659 ! convert error to percent divided by 10
1660       r8arr(11,nlevels) = current % meas % rh % error * 0.1
1661       r8arr(13,nlevels) = qz (current % meas % temperature % qc)
1662       r8arr(14,nlevels) = 1
1663       r8arr(15,nlevels) = 100
1664       r8arr(16,nlevels) = current % meas % temperature % error
1665
1666       current => current%next
1667
1668       ENDDO levels1
1669
1670       CALL UFBINT (bfout, r8arr, mxmn, nlevels, nlv, 'CAT POB PQM PPC PRC POE QOB QQM QPC QRC QOE TOB TQM TPC TRC TOE')
1671
1672       current => obs (loop_index) % surface
1673 !
1674       is_sound = -1
1675
1676       nlevels = 0
1677
1678       sid = obs (loop_index) % location % id
1679 !!    if ( index (sid,'47678') .ne. 0 ) then
1680 !     if ( trim(sid) .eq. '47678') then
1681 !       write(6,*) 'sid = ', TRIM(sid)
1682 !     endif
1683
1684 levels2:&
1685       DO WHILE (ASSOCIATED (current))
1686
1687       nlevels  = nlevels  + 1
1688       is_sound = is_sound + 1
1689
1690 !     if ( index (sid,'47678') .ne. 0 ) then
1691 !     if ( trim(sid) .eq. '47678') then
1692 !     if ( trim(sid) .eq. 'ZB678') then
1693 !       write(6,*) 'level = ',nlevels,' z = ',current % meas % height % data
1694 !       write(6,*) '   press  = ',current % meas % pressure % data
1695 !       write(6,*) '   dir = ',current % meas % direction % data,'  speed = ',current % meas % speed % data
1696 !       write(6,*) '   speed qc = ',current % meas % speed % qc,'  dir qc = ',current % meas % direction % qc
1697 !       write(6,*) '   t = ',current % meas % temperature % data
1698 !     endif
1699
1700       call assignv(current % meas % height % data, r8arr(1,nlevels))
1701       if ( r8arr(1,nlevels) /= bufrlib_missing ) r8arr(1,nlevels) = r8arr(1,nlevels)
1702       r8arr(2,nlevels) = qz (current % meas % height % qc)
1703       r8arr(3,nlevels) = 1
1704       r8arr(4,nlevels) = 100
1705       r8arr(5,nlevels) = current % meas % height % error
1706 REF: &
1707       SELECT CASE (cfm(fm))
1708       CASE ('GPSREF')
1709         if ( eps_equal( current % meas % dew_point % data, missing_r, 1.)) then
1710           r8arr(6,nlevels) = bufrlib_missing
1711         else
1712           r8arr(6,nlevels) = current % meas % dew_point % data
1713         endif
1714         r8arr(7,nlevels) = qz (current % meas % dew_point % qc)
1715         r8arr(8,nlevels) = 1
1716         r8arr(9,nlevels) = 100
1717         r8arr(10,nlevels) = current % meas % dew_point % error
1718         r8arr(11,nlevels) = bufrlib_missing
1719       CASE DEFAULT
1720         if ( eps_equal( current % meas % speed % data, missing_r, 1.) .or.  &
1721              eps_equal( current % meas % direction % data, missing_r, 1.) ) then
1722           r8arr(6,nlevels) = bufrlib_missing
1723           r8arr(10,nlevels) = bufrlib_missing
1724         else
1725           r8arr(6,nlevels) = current % meas % speed % data * (-sin(pi/180.* current % meas % direction % data))
1726           r8arr(10,nlevels) = current % meas % speed % data * (-cos(pi/180.* current % meas % direction % data))
1727         endif
1728         r8arr(7,nlevels) = qz (MAX(current % meas % direction % qc, current % meas % speed % qc))
1729         r8arr(8,nlevels) = 1
1730         r8arr(9,nlevels) = 100
1731         r8arr(11,nlevels) = current % meas % u % error
1732
1733       END SELECT REF
1734
1735       current => current%next
1736
1737       ENDDO levels2
1738
1739       SELECT CASE (cfm(fm))
1740       CASE ('GPSREF')
1741         CALL UFBINT (bfout, r8arr, mxmn, nlevels, nlv, 'ZOB ZQM ZPC ZRC ZOE ROB RQM RPC RRC ROE RFC')
1742       CASE DEFAULT
1743         CALL UFBINT (bfout, r8arr, mxmn, nlevels, nlv, 'ZOB ZQM ZPC ZRC ZOE UOB WQM WPC WRC VOB WOE')
1744       END SELECT
1745
1746       END SELECT UPA
1747
1748       SATWIND: &
1749       SELECT CASE (cfm(fm))        ! set the said parameter
1750       CASE ('SATWND')
1751         if (index(obs(loop_index) % location % id, 'GOES') .ne. 0) then
1752           r8arr(1,1) = 252
1753         else if (index(obs(loop_index) % location % id, 'MET') .ne. 0) then
1754           r8arr(1,1) = 54
1755         endif
1756
1757 !       current => obs (loop_index) % surface
1758 !  write(6,*) ' ID = ',obs(loop_index) % location % id,' name = ',obs(loop_index) % location % name
1759 !  write(6,*) 'lat = ',obs(loop_index) % location % latitude,' lon = ',obs(loop_index) % location % longitude
1760 !  write(6,*) 'p = ',current % meas % pressure % data, &
1761 !   ' speed = ',current % meas % speed % data,' dir = ',current % meas % direction % data
1762 !   write(6,*) 'said = ',r8arr(1,1)
1763
1764         CALL UFBINT (bfout, r8arr, 1, 1, nlv, 'SAID')
1765
1766       CASE ('QKSWND')
1767         r8arr(1,1) = 281
1768         CALL UFBINT (bfout, r8arr, 1, 1, nlv, 'SAID')
1769       END SELECT SATWIND
1770
1771
1772       CALL WRITSB (bfout)
1773
1774       CASE DEFAULT
1775 !       PRINT *, 'Erroneous cfm value,  fm = ',fm,' cfm = ',cfm(fm)
1776
1777      END SELECT SUBSET
1778
1779       ENDIF stations_valid
1780
1781 ! 3.3 Go to next record
1782 !     -----------------
1783
1784       ENDDO stations
1785
1786       CALL CLOSBF (bfout)
1787
1788 #else
1789       PRINT *, 'No Prepbufr data output, BUFR compilation is needed to turn it on !!!'
1790 #endif
1791
1792    RETURN
1793
1794 END SUBROUTINE output_prep
1795
1796 !-------------------------------------------------------
1797    FUNCTION stuff (c)
1798    REAL (KIND=8) :: stuff
1799
1800 !  Puts the first 8 characters in c into an 8 byte array
1801
1802 !  Original: James F. Drake, 09/22/05
1803 !  This function may not work on some of the Cray fortran
1804 !  processors, which store character variables in an acceptable
1805 !  but unusual way.
1806
1807    CHARACTER (LEN=40) c
1808    CHARACTER (LEN=1) d(8)
1809    integer (kind=1) :: x(8)
1810    REAL (KIND=8) :: y
1811    INTEGER clen, i
1812    EQUIVALENCE(x, y)
1813
1814    READ (c, '(8a1)') d
1815    x = IACHAR(d)
1816    stuff = y
1817    RETURN
1818
1819 END FUNCTION stuff
1820
1821 !-------------------------------------------------------
1822    FUNCTION t29(fm, sid)
1823    IMPLICIT NONE
1824    REAL (KIND=8) :: t29
1825    INTEGER :: fm
1826    character*40 sid
1827
1828 !     http://www.emc.ncep.noaa.gov/mmb/data_processing/prepbufr.doc/table_6.htm
1829 !     input report type T29
1830
1831       SELECT CASE (fm)
1832       CASE (12, 14)              ! synoptic
1833          t29 = 511
1834       CASE (13)             ! ship
1835          IF (TRIM(sid) /= 'SHIP') THEN
1836 !           ship with name
1837             t29 = 522
1838          ELSE
1839 !           ship without name
1840             t29 = 523
1841          ENDIF
1842       CASE (15:16)              ! METAR
1843             t29 = 512
1844       CASE (18)                 ! buoy
1845 !        What are possible identifiers for drifting buoys?
1846          IF (TRIM(sid) /= 'DRIB') THEN
1847 !           fixed buoy
1848             t29 = 561
1849          ELSE
1850 !           drifting buoy
1851             t29 = 562
1852          ENDIF
1853       CASE (19)               ! C-MAN
1854          t29 = 531
1855       CASE (32, 35)           ! sounding
1856          t29 = 11
1857       CASE (34, 38)           ! land mobil sounding
1858          t29 = 13
1859       CASE (33, 36)           ! ship sounding
1860          t29 = 22
1861       CASE (37)               ! dropsonde
1862          t29 = 22
1863       CASE (42, 96:97)        ! aircraft flight level
1864          t29 = 41
1865       CASE (86:87)            ! satellite
1866          t29 = 61
1867       CASE (88)               ! satellite-derived winds
1868          t29 = 63
1869       CASE (111)              ! assign 583 to gpspwv until we find what NCEP uses
1870          t29 = 583
1871       CASE (116)              ! assign 584 to gpsref
1872          t29 = 584
1873       CASE (125)              ! ssmi pw
1874          t29 = 65
1875       CASE (132)              ! wind profiler
1876          t29 = 71
1877       CASE (141)              ! modis
1878          t29 = 63
1879       CASE (281)              ! quikscat
1880          t29 = 582
1881       CASE DEFAULT
1882          t29 = 500            ! default to 500 (things we haven't added yet)
1883       END SELECT
1884 END FUNCTION t29
1885
1886 !-------------------------------------------------------
1887    FUNCTION typ (press, cfmfm, fm, ltemp, missing_r)
1888    IMPLICIT NONE
1889    REAL (KIND=8) :: typ
1890
1891 REAL :: press, missing_r
1892 INTEGER :: fm
1893 CHARACTER*6 :: cfmfm
1894 LOGICAL :: ltemp   ! temperature available?
1895
1896 !     http://www.emc.ncep.noaa.gov/mmb/data_processing/prepbufr.doc/table_2.htm
1897 !     PREPBUFR report type TYP
1898
1899    if ( cfmfm .eq. 'ADPSFC' .or. cfmfm .eq. 'SFCSHP' ) then
1900    IF ( ltemp ) THEN
1901       IF ( .not. eps_equal(press, missing_r, 1.) ) THEN
1902 !        TYP the PREPBUFR report type.
1903          SELECT CASE(cfmfm)
1904          CASE ('ADPSFC')
1905                typ = 181    ! per direction of Jim Bresch, NCAR
1906          CASE ( 'SFCSHP')
1907                typ = 180
1908          END SELECT
1909       ELSE
1910 !        No.
1911          press = bufrlib_missing
1912 !        Determine the code for TYP the PREPBUFR report type.
1913 !        Altimeter setting always missing.
1914          typ = 183
1915       ENDIF
1916    ELSE
1917 !     temperature, and hence TOB, is missing
1918 !     TYP
1919 !     Station pressure reported?
1920       IF (.not. eps_equal(press, missing_r, 1.) ) THEN
1921 !        Determine the code for TYP the PREPBUFR report type.
1922          SELECT CASE(cfmfm)
1923          CASE ('ADPSFC')
1924             typ = 181
1925          CASE ('SFCSHP')
1926             typ = 183
1927          END SELECT
1928       ELSE
1929 !        No. Determine the code for TYP the PREPBUFR report type.
1930 !        Altimeter setting always missing.
1931          typ = 183
1932          SELECT CASE(cfmfm)
1933          CASE ('ADPSFC', 'SFCSHP')
1934             typ = 183
1935          END SELECT
1936       ENDIF
1937    ENDIF
1938    else
1939      SELECT CASE (fm)
1940      CASE (32, 33)
1941         typ = 221
1942      CASE (34, 38)
1943         typ = 222
1944      CASE (35, 36)
1945         typ = 120
1946      CASE (37)
1947         typ = 232
1948      CASE (42)
1949         typ = 231
1950      CASE (86:87)        ! incomplete (but not necessary)
1951         typ = 164
1952      CASE (88)
1953         typ = 146
1954      CASE (96:97)
1955         typ = 230
1956      CASE (111)
1957         typ = 156       ! use goes entry for gpspw
1958      CASE (116)
1959         typ = 160       ! assign 160 to gpsref
1960      CASE (125)
1961         typ = 152       ! ssmi pw
1962      CASE (132)
1963         typ = 223
1964      CASE (135)         ! tc bogus
1965         typ = 111
1966      CASE (141)         ! use satwnd for modis
1967         typ = 146
1968      CASE (281)
1969         typ = 285
1970      CASE DEFAULT
1971         write(6,*) 'Cannot determine typ for FM = ',fm
1972      END SELECT
1973    endif
1974 END FUNCTION typ
1975
1976 !-------------------------------------------------------
1977 subroutine assignv(x, y)
1978    real (kind=8) :: y
1979    real x
1980
1981    if (eps_equal(x, missing_r, 1.)) then
1982       y = bufrlib_missing
1983    else
1984       y = x
1985    end if
1986 end subroutine assignv
1987
1988 !-------------------------------------------------------
1989    FUNCTION qz (qc)
1990    IMPLICIT NONE
1991    real (kind=8) :: qz
1992    INTEGER :: qc
1993
1994 !  convert qc flag to the ncep prepbufr qc value
1995
1996 if (qc .le. 128 .and. qc .ge. -5 ) then
1997   qz = 0
1998 else
1999   qz = 4    ! >= 4 is reject
2000 endif
2001
2002 END FUNCTION qz
2003
2004 END MODULE module_write