wrfv2_fire/var/da/da_obs/da_count_filtered_obs.inc

   1 subroutine da_count_filtered_obs (ptop, map, ds, phic, xlonc, truelat1, truelat2, &
   2    coarse_ix, coarse_jy, start_x, start_y)
   3
   4    !---------------------------------------------------------------------------
   5    ! Purpose: Scans intermediate Filtered Obs file,
   6    !          counts various obs type and writes on filtered_obs_unit
   7    !---------------------------------------------------------------------------
   8
   9    implicit none
  10
  11    real,              intent(in) :: ptop, ds
  12    real,              intent(in) :: phic, xlonc, truelat1, truelat2
  13    integer,           intent(in) :: coarse_ix, coarse_jy
  14    real,              intent(in) :: start_x, start_y
  15    integer,           intent(in) :: map
  16
  17    integer                      :: i, iost, fm
  18    type (multi_level_type)      :: platform
  19    real                         :: height_error
  20    integer                      :: nlocal(num_ob_indexes)
  21    integer                      :: num_others
  22
  23    integer                        :: iunit, files, total_obs
  24    integer                        :: maxnes, numc, idd
  25    character(len=filename_len)    :: filename
  26
  27
  28    if (trace_use) call da_trace_entry("da_count_filtered_obs")
  29
  30    nlocal(:) = 0
  31    num_others = 0
  32
  33    call da_get_unit(iunit)
  34
  35    ! Loop over all data files
  36    do files = 0, num_procs-1
  37       write(unit=filename, fmt='(a,i3.3)') 'filtered_obs.',files
  38       open(unit=iunit, file= trim(filename), form='formatted',iostat=iost)
  39       if (iost /= 0) call da_error(__FILE__,__LINE__, (/"Cannot open "//filename/))
  40
  41       !  loop over records
  42
  43       reports: do
  44          ! read station general info
  45
  46          read (iunit, fmt = fmt_info, iostat = iost) &
  47             platform%info%platform,    &
  48             platform%info%date_char,   &
  49             platform%info%name,        &
  50             platform%info%levels,      &
  51             platform%info%lat,         &
  52             platform%info%lon,         &
  53             platform%info%elv,         &
  54             platform%info%id
  55
  56          if (iost /= 0) then
  57             !write (0,'(/A,I9)') ' end OF OBS unit: ',iunit
  58             !write (0,'(/A,I9)') ' iostat:          ',iost
  59             exit reports
  60          end if
  61
  62          ! Read surface Info
  63
  64          read (iunit, fmt = fmt_srfc)  &
  65             platform%loc%slp%inv, platform%loc%slp%qc, &
  66             platform%loc%slp%error,                    &
  67             platform%loc%pw%inv, platform%loc%pw%qc,   &
  68             platform%loc%pw%error
  69          ! levels < 1 and not GPSPW, go back to reports
  70
  71          if ((platform%info%levels < 1) .AND.            &
  72              (index(platform%info%platform, 'GPSPW') <= 0)) then
  73               cycle reports
  74          end if
  75          read(platform%info%platform(4:6), '(I3)') fm
  76
  77          ! read each level
  78          do i = 1, platform%info%levels
  79             platform%each (i) = each_level_type(missing_r, missing, -1.0, & ! height
  80                field_type(missing_r, missing_data, missing_r, missing_r, missing_r), & ! u
  81                field_type(missing_r, missing_data, missing_r, missing_r, missing_r), & ! v
  82                field_type(missing_r, missing_data, missing_r, missing_r, missing_r), & ! p
  83                field_type(missing_r, missing_data, missing_r, missing_r, missing_r), & ! t
  84                field_type(missing_r, missing_data, missing_r, missing_r, missing_r), & ! q
  85                field_type(missing_r, missing_data, missing_r, missing_r, missing_r), & ! rh
  86                field_type(missing_r, missing_data, missing_r, missing_r, missing_r), & ! td
  87                field_type(missing_r, missing_data, missing_r, missing_r, missing_r))  ! speed
  88
  89             read (unit = iunit, fmt = trim (fmt_each)) &
  90                platform%each(i)%p%inv, platform%each(i)%p%qc, platform%each(i)%p%error, &
  91                platform%each(i)%speed%inv, platform%each(i)%speed%qc,                   &
  92                platform%each(i)%speed%error,                                            &
  93                platform%each(i)%v%inv, platform%each(i)%v%qc, platform%each(i)%v%error, &
  94                platform%each(i)%height,                                                 &
  95                platform%each(i)%height_qc,                                              &
  96                height_error,                                                            &
  97                platform%each(i)%t%inv, platform%each(i)%t%qc, platform%each(i)%t%error, &
  98                platform%each(i)%td%inv, platform%each(i)%td%qc, platform%each(i)%td%error, &
  99                platform%each(i)%rh%inv, platform%each(i)%rh%qc, platform%each(i)%rh%error
 100          end do
 101
 102          if (platform%info%levels < 1) then
 103             if (fm /= 111) then
 104                cycle reports
 105             end if
 106          end if
 107          select case(fm)
 108
 109          case (12) ;
 110
 111             if (.not.use_synopobs) cycle reports
 112             nlocal(synop) = nlocal(synop) + 1
 113
 114          case (13, 17) ;                  ! ships
 115
 116             if (.not.use_shipsobs) cycle reports
 117             nlocal(ships)  = nlocal(ships)  + 1
 118
 119          case (15:16) ;
 120
 121             if (.not.use_metarobs) cycle reports
 122             nlocal(metar) = nlocal(metar) + 1
 123
 124          case (32:34) ;
 125
 126             if (.not.use_pilotobs) cycle reports
 127             nlocal(pilot) = nlocal(pilot) + 1
 128
 129          case (35:38) ;
 130             if (.not.use_soundobs) cycle reports
 131             nlocal(sound) = nlocal(sound) + 1
 132
 133          case (161) ;
 134             if (.not.use_mtgirsobs) cycle reports
 135             nlocal(mtgirs) = nlocal(mtgirs) + 1
 136
 137          case (86) ;
 138
 139             if (.not.use_satemobs) cycle reports
 140
 141             ! Reject cloudy satem obs.
 142
 143             if (platform%loc%pw%inv > 10.0) then
 144                cycle reports
 145             end if
 146
 147             nlocal(satem) = nlocal(satem) + 1
 148
 149          case (88)    ;
 150             ! Geostationary or Polar orbitting Satellite AMVs:
 151
 152             if (index(platform%info%name, 'MODIS') > 0 .or. &
 153                index(platform%info%name, 'modis') > 0)  then
 154                if (.not.use_polaramvobs) cycle reports
 155                nlocal(polaramv) = nlocal(polaramv) + 1
 156             else
 157                if (.not.use_geoamvobs) cycle reports
 158                nlocal(geoamv) = nlocal(geoamv) + 1
 159             end if
 160
 161          case (42,96:97) ;
 162
 163             if (.not.use_airepobs) cycle reports
 164             nlocal(airep) = nlocal(airep) + 1
 165
 166          case (101) ;
 167             if (.not.use_tamdarobs) cycle reports
 168             nlocal(tamdar) = nlocal(tamdar) + 1
 169
 170          case (111) ;
 171
 172             if (.not.use_gpspwobs) cycle reports
 173             nlocal(gpspw) = nlocal(gpspw) + 1
 174
 175          case (116) ;
 176
 177             if (.not.use_gpsrefobs) cycle reports
 178             nlocal(gpsref) = nlocal(gpsref) + 1
 179
 180          case (121) ;
 181             ! SSM/T1 temperatures:
 182
 183             if (.not.use_ssmt1obs) cycle reports
 184             nlocal(ssmt1) = nlocal(ssmt1) + 1
 185
 186          case (122) ;
 187             ! SSM/T2 relative humidities:
 188
 189             if (.not.use_ssmt2obs) cycle reports
 190             nlocal(ssmt2) = nlocal(ssmt2) + 1
 191
 192          case (281)    ;
 193             ! Scatterometer:
 194
 195             if (.not.use_qscatobs) cycle reports
 196             nlocal(qscat)  = nlocal(qscat)  + 1
 197
 198          case (132) ;
 199
 200             if (.not.use_profilerobs) cycle reports
 201             nlocal(profiler) = nlocal(profiler) + 1
 202
 203          case (135) ;
 204
 205             if (.not.use_bogusobs) cycle reports
 206             nlocal(bogus) = nlocal(bogus) + 1
 207
 208          case (18,19) ;             ! bouy
 209
 210             if (.not.use_buoyobs) cycle reports
 211             nlocal(buoy)  = nlocal(buoy)  + 1
 212
 213          case (133) ;      !  AIRS retrievals
 214             if (.not.use_airsretobs) cycle reports
 215             nlocal(airsr) = nlocal(airsr) + 1
 216
 217          case default;
 218             num_others = num_others + 1
 219             write(unit=message(1), fmt='(a)') 'unsaved obs found:'
 220             write(unit=message(2), fmt='(2a)') &
 221                'platform%info%platform=', platform%info%platform
 222             write(unit=message(3), fmt='(a, i3)') &
 223                'platform%info%levels=', platform%info%levels
 224             call da_warning(__FILE__,__LINE__,message(1:3))
 225          end select
 226       end do reports                  !  Loop over reports
 227       close (iunit)
 228    end do               !  Loop over all data files
 229    call da_free_unit (iunit)
 230
 231    ! write counts in the header
 232
 233    total_obs = nlocal(synop) + nlocal(metar) + nlocal(ships) + &
 234       nlocal(buoy) + nlocal(sound) + nlocal(airep) + nlocal(pilot) + &
 235       nlocal(geoamv) + nlocal(polaramv) + nlocal(gpspw) + nlocal(gpsref) + &
 236       nlocal(profiler) + nlocal(qscat) + nlocal(ssmt1) + nlocal(ssmt2) +  &
 237       nlocal(satem)  + nlocal(bogus) +  nlocal(airsr) + nlocal(mtgirs) + nlocal(tamdar) + num_others
 238
 239    write  (unit = filtered_obs_unit, fmt = '(A,I7,A,F8.0,A)')    &
 240       "TOTAL =",total_obs, ", MISS. =",missing_r,","
 241
 242    ! Write other counts
 243
 244    write  (unit = filtered_obs_unit, fmt = '(6(A,I7,A))')    &
 245       "SYNOP =",nlocal(synop),", ", &
 246       "METAR =",nlocal(metar),", ", &
 247       "SHIP  =",nlocal(ships),", ", &
 248       "BUOY  =",nlocal(buoy),", ", &
 249       "TEMP  =",nlocal(sound),", ", &
 250       "AIREP =",nlocal(airep),", ", &
 251       "PILOT =",nlocal(pilot),", ", &
 252       "GeAMV =",nlocal(geoamv),", ", &
 253       "PoAMV =",nlocal(polaramv),", ", &
 254       "GPSPW =",nlocal(gpspw),", ", &
 255       "GPSRF =",nlocal(gpsref),", ", &
 256       "PROFL =",nlocal(profiler),", ", &
 257       "QSCAT =",nlocal(qscat),", ", &
 258       "SSMT1 =",nlocal(ssmt1),", ", &
 259       "SSMT2 =",nlocal(ssmt2),", ", &
 260       "SATEM =",nlocal(satem), ", ", &
 261       "BOGUS =",nlocal(bogus),", ", &
 262       "AIRSR =",nlocal(airsr),", ", &
 263       "MTGIRS=",nlocal(mtgirs),", ", &
 264       "TAMDAR=",nlocal(tamdar),", ", &
 265       "OTHER =",num_others,", "
 266
 267    ! write reference state info
 268
 269    write (unit = filtered_obs_unit, fmt = '(A,F7.2,A,F7.2,4(A,F7.2),A)') &
 270       "PHIC  =", phic,", XLONC =", xlonc,", TRUE1 =", truelat1,&
 271       ", TRUE2 =",truelat2, ", XIM11 =", start_y, ", XJM11 =", start_x, ","
 272    write (unit = filtered_obs_unit, fmt = '(2(A,F7.2),A,F7.0,A,F7.0,A)') &
 273       "TS0   =",  base_temp, ", TLP   =", base_lapse, &
 274       ", PTOP  =",  ptop,", PS0   =",  base_pres,","
 275
 276    ! write domain info
 277
 278    !  hardwire following variables
 279       maxnes = 2; numc = 1
 280    if ( coarse_ix == ide+1 .and. coarse_jy == jde+1 ) then
 281       idd = 2
 282    else
 283       idd = 1
 284    end if
 285    write (unit = filtered_obs_unit, fmt = '(5(A,I7),A)') &
 286       "IXC   =", coarse_jy   ,", JXC   =", coarse_ix   ,", IPROJ =", map,&
 287       ", IDD   =", idd,", MAXNES=",maxnes,","
 288    write (unit = filtered_obs_unit, fmt = '(A,10(:,I7,A))')  &
 289       "NESTIX=", coarse_jy, ", ", jde+1, ", "
 290    write (unit = filtered_obs_unit, fmt = '(A,10(:,I7,A))')  &
 291       "NESTJX=", coarse_ix, ", ", ide+1, ", "
 292    write (unit = filtered_obs_unit, fmt = '(A,10(:,I7,A))')  &
 293       "NUMC  =",(numc ,", ", i = 1, maxnes)
 294    write (unit = filtered_obs_unit, fmt = '(A,10(:,F7.2,A))')&
 295       "DIS   =",(ds/1000.0,", ",i = 1, maxnes)
 296    write (unit = filtered_obs_unit, fmt = '(A,10(:,I7,A))')  &
 297       "NESTI=", 1, ", ", int(start_y), ", "
 298    write (unit = filtered_obs_unit, fmt = '(A,10(:,I7,A))')  &
 299       "NESTJ=", 1, ", ", int(start_x), ", "
 300
 301    ! write variable names and unit
 302
 303    write (unit = filtered_obs_unit, fmt = '(A)') &
 304       "INFO  = PLATFORM, DATE, NAME, LEVELS, LATITUDE, LONGITUDE, ELEVATION, ID."
 305    write (unit = filtered_obs_unit, fmt = '(A)') &
 306       "SRFC  = SLP, PW (DATA,QC,ERROR)."
 307    write (unit = filtered_obs_unit, fmt = '(A)') &
 308       "EACH  = PRES, SPEED, DIR, HEIGHT, TEMP, DEW PT, HUMID (DATA,QC,ERROR)*LEVELS."
 309
 310    ! write format info
 311
 312    write (unit = filtered_obs_unit, fmt = '(2A)') 'INFO_fmt = ', trim (fmt_info)
 313    write (unit = filtered_obs_unit, fmt = '(2A)') 'SRFC_fmt = ', trim (fmt_srfc)
 314    write (unit = filtered_obs_unit, fmt = '(2A)') 'EACH_fmt = ', trim (fmt_each)
 315
 316    ! write end of header record
 317
 318    if (write_mod_filtered_obs) then    !cys_change
 319       write (unit = filtered_obs_unit, fmt = '(A)') &
 320          "MODIFIED FILTERED OBS #--------------------------------------------------------#"
 321    else
 322       write (unit = filtered_obs_unit, fmt = '(A)') &
 323          "FILTERED OBS #-----------------------------------------------------------------#"
 324    end if
 325
 326    if (trace_use) call da_trace_exit("da_count_filtered_obs")
 327
 328 end subroutine da_count_filtered_obs
 329
 330