var/da/da_setup_structures/da_setup_obs_structures_madis.inc

   1 subroutine da_setup_obs_structures_madis( ob, iv)
   2
   3 !------------------------------------------------------------------------------
   4 ! PURPOSE: Define, allocate and read of observation structure.
   5 !
   6 ! METHOD:  Define, allocate and read of observation structure.
   7 !
   8 ! HISTORY: 05-21-03  M. Barth         The setup_obs routine to ingest data from
   9 !                                     the FSL MADIS netCDF files.
  10 !          10-20-03  M. Barth         - Changed calls to get SATOB data for the
  11 !                                       changes made in the MADIS API when the
  12 !                                       true SATWND dataset was added.  (This code
  13 !                                       originally accessed a dummy SATW dataset.)
  14 !                                     - Changed max_map from 150 to 2000 to
  15 !                                       accommodate full time window options.
  16 !          04-01-04  M. Barth         Added MADIS COOP surface dataset.
  17 !          08-27-04  M. Barth         - Added sonde_sfc stuff.
  18 !                                     - Initialize new count_obs entries that we
  19 !                                       don't use (radar, profiler, buoy).  Note
  20 !                                       that, as with the BUFR interface, we
  21 !                                       continue to lump profiler in with "sound"
  22 !                                       and buoy with "ships" -- as the downstream
  23 !                                       processing is identical regardless of
  24 !                                       the designation, and we offer MADIS-level
  25 !                                       control of these datasets, this is simpler
  26 !                                       than restructuring this code.
  27 !          02-06-06 M. Barth          Changes for WRFVAR version 2.1.
  28 !                                     - Renamed all "fsl" references to "madis".
  29 !                                       (Exception:  "FSL" is still used as an
  30 !                                        argument to MINIT.)
  31 !                                     - Changed xb to xp.
  32 !                                     - Removed xbx arg, so analysis time is
  33 !                                       now determined by the NAMELIST ANALYSIS_DATE
  34 !                                       variable.
  35 !                                     - Changed satob to geoamv.
  36 !                                     - Changed max_sfc from 50K to 200K.
  37 !                                     - Changed max_raoblev from 256 to 405.
  38 !                                     - Miscellaneous other minor changes to
  39 !                                       conform to current paradigm.
  40 !                                     - Added protection against invalid lat/lons,
  41 !                                       as well as protection against latitude
  42 !                                       at the poles -- some projections can't
  43 !                                       handle this, and the WRFVAR routines
  44 !                                       that do the coordinate conversions aren't
  45 !                                       smart enough to trap this case.  [Are you
  46 !                                       surprised?]
  47 !          08-07-07 M. Barth          - Fixed problem introduced in 2005 version
  48 !                                       of WRFVAR by filling in the ob_num
  49 !                                       and current_ob_time data structures in iv.
  50 !          02-06-06 M. Barth          Changes for WRFVAR version 3:
  51 !                                     - Removed xp argument.
  52 !                                     - Changed DA_ll_to_xy call to da_llxy.
  53 !                                     - Added trace usage.
  54 !                                     - Changed ob_type structure to iv_type.
  55 !                                     - Changes in iv structure:
  56 !                                       - current_ob_time -> time
  57 !                                       - New info structure
  58 !                                     - Fill in global vs. local obs counts
  59 !                                       into iv%info.  This is *presumably* the
  60 !                                       replacement for iv%ob_num.
  61 !                                     - Fill in iv%info with stuff that used
  62 !                                       to go into iv%<dataset>%info.
  63 !                                     - Since DA_Constants and DA_Read_Namelist no
  64 !                                       longer exist, the declaration of, and
  65 !                                       input of, the MADIS namelist variables
  66 !                                       have been moved to this routine.
  67 !                                     - Changed print_detail to print_detail_obs.
  68 !                                     - Added duplication logic (see explanation
  69 !                                       in PROCESS_DATASET).
  70 !          11-28-08 M. Barth          Since I finally got the darn thing to
  71 !                                     compile, I debugged the massive changes
  72 !                                     in version 3.
  73 !
  74 ! PARENT_MODULE: da_setup_structures
  75 !------------------------------------------------------------------------------
  76
  77    use da_tools, only : da_llxy
  78    IMPLICIT NONE
  79
  80 ! ARGUMENTS:
  81
  82    TYPE ( y_type), INTENT(OUT)    :: ob          ! Observation structure.
  83    TYPE (iv_type), INTENT(INOUT)  :: iv          ! O-B structure.
  84
  85 #ifdef MADIS
  86
  87 ! LOCAL:
  88
  89 ! Primarily used by DA routines:
  90
  91    TYPE (multi_level_type)      :: platform
  92
  93    INTEGER                      :: fm
  94
  95    LOGICAL                      :: outside, outside_all
  96
  97 ! Primarily used by MADIS ingest:
  98
  99 ! Constants:
 100
 101                                                          ! Max number of input reports:
 102    INTEGER*4, PARAMETER         :: max_sfc   = 200000, & !  Surface
 103                                    max_raob  =   1000, & !  Radiosonde
 104                                    max_npn   =    160, & !  NOAA Profiler Network (NPN) winds
 105                                    max_acars =  60000, & !  ACARS
 106                                    max_map   =   2000, & !  Multi-Agency profilers (MAP)
 107                                    max_satw  =  60000, & !  Satellite winds
 108
 109                                    max_sta   = max(max_sfc,max_raob,max_npn, &
 110                                                    max_acars,max_map,max_satw), &
 111
 112                                                          ! Max number of input levels:
 113                                    max_raoblev = 405, &  !  Radiosonde
 114                                    max_npnlev_madis = 64, & ! NPN MADIS database
 115                                    max_npnlev_awips = 43, &!NPN AWIPS database
 116                                    max_maplev  = 202, &  !  MAP
 117
 118                                    success_p =      0, & ! Successful status return
 119                                    qcall =         99, & ! Only return data passing all QC
 120                                    max_provider = 100    ! Maximum number of data providers
 121    REAL*4, PARAMETER            :: miss_p    = 999999.0  ! Missing data flag
 122
 123 ! Input arrays:
 124
 125    CHARACTER*10, ALLOCATABLE    :: stanam(:),provdr(:)
 126    CHARACTER*9, ALLOCATABLE     :: timeob(:)
 127    CHARACTER*1, ALLOCATABLE     :: qcd(:)
 128    INTEGER*4, ALLOCATABLE       :: wmoid(:),qca(:),qcr(:),nlevels(:),keep(:)
 129    REAL*4, ALLOCATABLE          :: lat(:),lon(:),elev(:),p(:),q(:),t(:),z(:), &
 130                                    dd(:),ff(:),pw(:),slp(:),levtype(:)
 131
 132 ! Miscellaneous:
 133
 134    INTEGER*4                    :: nsta,nobs,istatus,i,j,ntmp,nlev,ntmp_r, &
 135                                    ntmp_n,ntmp_m,keep_r(max_raob), recwin, &
 136                                    keep_n(max_npn),keep_m(max_map),minbck,minfwd, &
 137                                    nsatw_vartype,satw_var_len(5),pvdriostat,n,iost, &
 138                                    nlocal(num_ob_indexes),ntotal(num_ob_indexes), &
 139                                    nl_raob,nt_raob,nl_npn,nt_npn,nl_map,nt_map, &
 140                                    nlocal_out,ntotal_out
 141    CHARACTER*9                  :: atime
 142    CHARACTER*13                 :: tstr
 143    CHARACTER*10                 :: dpname
 144    CHARACTER*255                :: pvdrfile
 145    INTEGER*4, PARAMETER         :: lun = 19     ! LUN for namelist and data provider files
 146    CHARACTER*5                  :: satw_var(2,5)
 147    CHARACTER*4                  :: pvdriotype
 148    LOGICAL                      :: map_init,npn_init
 149    CHARACTER*20                 :: namelist_file
 150    CHARACTER*29                 :: errmsg(0:8)
 151    INTEGER*4                    :: mfbdebug
 152    INTEGER*4                    :: sfclev(max_map)
 153
 154 !------------------------------------------------------------------------------
 155 ! Namelist options for ingest of MADIS/AWIPS datasets into 3DVAR
 156 !------------------------------------------------------------------------------
 157 ! The user can select either the MADIS or AWIPS database for each dataset that's
 158 ! ingest by the WRFVAR.  Each dataset's time window can be independently set as
 159 ! well, which allows the user to select the number of minutes before and after
 160 ! the nominal time to be used for the window, and how duplicates should be
 161 ! handled.  For the METAR, PROFILER and GEOAMV datasets, the user can also
 162 ! choose to ingest data from all providers of the METAR and PROFILER datasets
 163 ! (e.g., ASOS) as well as mesonets for METAR, and from operational GOES
 164 ! products from both satellites for GEOAMV), or the user can specify which
 165 ! providers should be used.  For more details on how and why to select
 166 ! different settings, see the README.MADIS_WRFVAR file in the
 167 ! var directory.
 168 !
 169 ! The meaning of each variable is the same for all of the datasets, so full
 170 ! comments are only included for the METAR dataset below, and where the variables
 171 ! are specific to only one or two datasets.
 172 !------------------------------------------------------------------------------
 173
 174 ! General
 175
 176    INTEGER*4     :: madis_debug               ! 0 - no debug
 177                                               ! 1 - print all obs
 178                                               ! 2 - print all obs, stop after obs ingest
 179
 180 ! METAR dataset
 181
 182    CHARACTER*6   :: madis_metar_db            ! Database:  'FSL' or 'AWIPS'
 183    INTEGER*4     :: madis_metar_minbck        ! Number of minutes relative to nominal
 184                                               ! time at which to start time window
 185    INTEGER*4     :: madis_metar_minfwd        ! Number of minutes relative to nominal
 186                                               ! time at which to end time window
 187                                               ! (max(minfwd-minbck) is 180)
 188    INTEGER*4     :: madis_metar_recwin        ! 1 - return one record per fixed station,
 189                                               !     closest to nominal time
 190                                               ! 2 - return one record per fixed station,
 191                                               !     closest to start of window
 192                                               ! 3 - return one record per fixed station,
 193                                               !     closest to end of window
 194                                               ! 4 - return all records within *window*
 195    LOGICAL       :: madis_metar_all_providers ! .TRUE./.FALSE. - use all providers
 196    CHARACTER*255 :: madis_metar_provider_list ! Full path to file with list of providers,
 197                                               ! one per line; must be specified if
 198                                               ! all_providers is false
 199
 200 ! SHIPS dataset
 201
 202    CHARACTER*6   :: madis_ships_db            ! Database:  'FSL' or 'AWIPS'
 203    INTEGER*4     :: madis_ships_minbck        ! Time window start
 204    INTEGER*4     :: madis_ships_minfwd        ! Time window end
 205    INTEGER*4     :: madis_ships_recwin        ! How to handle duplicate records
 206
 207 ! GPSPW dataset
 208
 209    CHARACTER*6   :: madis_gpspw_db            ! Database:  'FSL' or 'AWIPS'
 210    INTEGER*4     :: madis_gpspw_minbck        ! Time window start
 211    INTEGER*4     :: madis_gpspw_minfwd        ! Time window end
 212    INTEGER*4     :: madis_gpspw_recwin        ! How to handle duplicate records
 213
 214 ! SOUND dataset
 215
 216    CHARACTER*6   :: madis_sound_db            ! Database:  'FSL' or 'AWIPS'
 217    INTEGER*4     :: madis_sound_minbck        ! Time window start
 218    INTEGER*4     :: madis_sound_minfwd        ! Time window end
 219    INTEGER*4     :: madis_sound_recwin        ! How to handle duplicate records
 220
 221 ! PROFILER dataset
 222
 223    CHARACTER*6   :: madis_profiler_db            ! Database:  'FSL' or 'AWIPS'
 224    INTEGER*4     :: madis_profiler_minbck        ! Time window start
 225    INTEGER*4     :: madis_profiler_minfwd        ! Time window end
 226    INTEGER*4     :: madis_profiler_recwin        ! How to handle duplicate records
 227    LOGICAL       :: madis_profiler_all_providers ! .TRUE./.FALSE. - use all providers
 228    CHARACTER*255 :: madis_profiler_provider_list ! Full path to file with list of providers
 229
 230 ! AIREP dataset
 231
 232    CHARACTER*6   :: madis_airep_db            ! Database:  'FSL' or 'AWIPS'
 233    INTEGER*4     :: madis_airep_minbck        ! Time window start
 234    INTEGER*4     :: madis_airep_minfwd        ! Time window end
 235    INTEGER*4     :: madis_airep_recwin        ! How to handle duplicate records
 236
 237 ! GEOAMV dataset
 238
 239    CHARACTER*6   :: madis_geoamv_db            ! Database:  'FSL' or 'AWIPS'
 240    INTEGER*4     :: madis_geoamv_minbck        ! Time window start
 241    INTEGER*4     :: madis_geoamv_minfwd        ! Time window end
 242    INTEGER*4     :: madis_geoamv_recwin        ! How to handle duplicate records
 243    LOGICAL       :: madis_geoamv_all_providers ! .TRUE./.FALSE. - use default providers
 244    CHARACTER*255 :: madis_geoamv_provider_list ! Full path to file with list of providers
 245    LOGICAL       :: madis_geoamv_wv            ! .TRUE./.FALSE. - use water vapor winds
 246    LOGICAL       :: madis_geoamv_ir            ! .TRUE./.FALSE. - use infrared winds
 247    LOGICAL       :: madis_geoamv_vis           ! .TRUE./.FALSE. - use visible winds
 248    LOGICAL       :: madis_geoamv_s10           ! .TRUE./.FALSE. - use sounder channel 10 winds
 249    LOGICAL       :: madis_geoamv_s11           ! .TRUE./.FALSE. - use sounder channel 11 winds
 250
 251 ! Note that the meaning of "madis_geoamv_all_providers" is different for GEOAMV than
 252 ! for the other datasets.  If the variable is .TRUE., the *default* providers will be
 253 ! selected, i.e., operational products from both GOES satellites.
 254 !------------------------------------------------------------------------------
 255 ! End of MADIS/AWIPS WRFVAR options
 256 !------------------------------------------------------------------------------
 257 ! Namelist contents :
 258
 259    NAMELIST/MADIS_GENERAL/madis_debug
 260    NAMELIST/MADIS_METAR/madis_metar_db,madis_metar_minbck,madis_metar_minfwd,madis_metar_recwin, &
 261                         madis_metar_all_providers,madis_metar_provider_list
 262    NAMELIST/MADIS_SHIPS/madis_ships_db,madis_ships_minbck,madis_ships_minfwd,madis_ships_recwin
 263    NAMELIST/MADIS_GPSPW/madis_gpspw_db,madis_gpspw_minbck,madis_gpspw_minfwd,madis_gpspw_recwin
 264    NAMELIST/MADIS_SOUND/madis_sound_db,madis_sound_minbck,madis_sound_minfwd,madis_sound_recwin
 265    NAMELIST/MADIS_PROFILER/madis_profiler_db,madis_profiler_minbck,madis_profiler_minfwd,&
 266                            madis_profiler_recwin,madis_profiler_all_providers,&
 267                            madis_profiler_provider_list
 268    NAMELIST/MADIS_AIREP/madis_airep_db,madis_airep_minbck,madis_airep_minfwd,madis_airep_recwin
 269    NAMELIST/MADIS_GEOAMV/madis_geoamv_db,madis_geoamv_minbck,madis_geoamv_minfwd,&
 270                          madis_geoamv_recwin,madis_geoamv_all_providers,&
 271                          madis_geoamv_provider_list,madis_geoamv_wv,madis_geoamv_ir,&
 272                          madis_geoamv_vis,madis_geoamv_s10,madis_geoamv_s11
 273
 274    DATA ERRMSG/'Opening namelist.3dvar_madis',&
 275                'Reading madis_general namelist',&
 276                'Reading madis_metar namelist',&
 277                'Reading madis_ships namelist',&
 278                'Reading madis_gpspw namelist',&
 279                'Reading madis_sound namelist',&
 280                'Reading madis_profiler namelist',&
 281                'Reading madis_airep namelist',&
 282                'Reading madis_geoamv namelist'/
 283
 284 ! Note that all of the variables sent into MADIS routines are explicitly
 285 ! declared INTEGER*4 and REAL*4, rather than letting the generic kind be
 286 ! established based on the choices used in the WRF options.  On at least
 287 ! one machine, the WRF-default for REAL is REAL*8, but the MADIS library
 288 ! uses REAL*4.  To stay on the safe side, we'll declare all the int/real
 289 ! MADIS variables to what we want.
 290
 291    if (trace_use) CALL da_trace_entry("da_setup_obs_structures_madis")
 292
 293 !------------------------------------------------------------------------------
 294 ! MADIS-related initialization
 295 !------------------------------------------------------------------------------
 296 ! Note that any variables prefixed with "madis_" are namelist variables read in
 297 ! from namelist.3dvar_madis.  These are used to give the user the ability to
 298 ! set all relevant MADIS options for the various datasets.
 299
 300    write(6,*)' -----------------------------------------------------'
 301    write(6,*)'       Read WRFVAR namelist options for MADIS datasets'
 302    write(6,*)' -----------------------------------------------------'
 303    write(6,*)
 304
 305 !----------------------------------------------------------------------------
 306 ! Set default values.  These represent reasonable choices for an hourly data
 307 ! assimilation cycle, using all datasets, for MADIS database, without
 308 ! debugging info.
 309 !----------------------------------------------------------------------------
 310    madis_debug = 0
 311
 312    madis_metar_db = 'FSL'
 313    madis_metar_minbck = -30
 314    madis_metar_minfwd = 29
 315    madis_metar_recwin = 4
 316    madis_metar_all_providers = .TRUE.
 317
 318    madis_ships_db = 'FSL'
 319    madis_ships_minbck = -30
 320    madis_ships_minfwd = 29
 321    madis_ships_recwin = 4
 322
 323    madis_gpspw_db = 'FSL'
 324    madis_gpspw_minbck = -30
 325    madis_gpspw_minfwd = 29
 326    madis_gpspw_recwin = 4
 327
 328    madis_sound_db = 'FSL'
 329    madis_sound_minbck = -60
 330    madis_sound_minfwd = 60
 331    madis_sound_recwin = 3
 332
 333    madis_profiler_db = 'FSL'
 334    madis_profiler_minbck = -30
 335    madis_profiler_minfwd = 29
 336    madis_profiler_recwin = 4
 337    madis_profiler_all_providers = .TRUE.
 338
 339    madis_airep_db = 'FSL'
 340    madis_airep_minbck = -30
 341    madis_airep_minfwd = 29
 342    madis_airep_recwin = 4
 343
 344    madis_geoamv_db = 'FSL'
 345    madis_geoamv_minbck = -150
 346    madis_geoamv_minfwd = 30
 347    madis_geoamv_recwin = 4
 348    madis_geoamv_all_providers = .TRUE.
 349    madis_geoamv_wv = .TRUE.
 350    madis_geoamv_ir = .TRUE.
 351    madis_geoamv_vis = .FALSE.
 352    madis_geoamv_s10 = .FALSE.
 353    madis_geoamv_s11 = .FALSE.
 354
 355 !----------------------------------------------------------------------------
 356 !  Open namelist file.
 357 !----------------------------------------------------------------------------
 358    namelist_file = 'namelist.3dvar_madis'
 359    write(6,*)' WRFVAR MADIS datasets namelist options used are in: ', &
 360              namelist_file
 361    iost = 0
 362    i = 0
 363
 364    write(6,*) ' Opening WRFVAR MADIS namelist input file...'
 365
 366    open(file=namelist_file,unit=lun,status='old', &
 367         access='sequential',form='formatted',action='read', &
 368         err=8500,iostat=iost)
 369
 370 !----------------------------------------------------------------------------
 371 !  Read namelist and close.
 372 !----------------------------------------------------------------------------
 373    iost = 0
 374    i = i + 1
 375    read(unit=lun,nml=madis_general,err=8500,iostat=iost)
 376    i = i + 1
 377    read(unit=lun,nml=madis_metar,err=8500,iostat=iost)
 378    i = i + 1
 379    read(unit=lun,nml=madis_ships,err=8500,iostat=iost)
 380    i = i + 1
 381    read(unit=lun,nml=madis_gpspw,err=8500,iostat=iost)
 382    i = i + 1
 383    read(unit=lun,nml=madis_sound,err=8500,iostat=iost)
 384    i = i + 1
 385    read(unit=lun,nml=madis_profiler,err=8500,iostat=iost)
 386    i = i + 1
 387    read(unit=lun,nml=madis_airep,err=8500,iostat=iost)
 388    i = i + 1
 389    read(unit=lun,nml=madis_geoamv,err=8500,iostat=iost)
 390    close(lun)
 391
 392    write(6,*) ' WRFVAR MADIS namelist input file successfully read.'
 393
 394 !------------------------------------------------------------------------------
 395 !  The input will be validated when the user's choices are employed in the
 396 !  da_setup_obs_structures_madis routine, and the MADIS library routines that
 397 !  it calls.  Therefore, don't duplicate the validation here.
 398 !------------------------------------------------------------------------------
 399
 400 !----------------------------------------------------------------------------
 401 !  Print namelist.
 402 !----------------------------------------------------------------------------
 403
 404    if (print_detail_obs) then
 405
 406       write(6,*) ' Printing contents of WRFVAR MADIS namelist input file...'
 407
 408       write(6,*)
 409       write(6,*)'madis_debug                  = ',madis_debug
 410       write(6,*)
 411       write(6,*)'madis_metar_db               = ',madis_metar_db
 412       write(6,*)'madis_metar_minbck           = ',madis_metar_minbck
 413       write(6,*)'madis_metar_minfwd           = ',madis_metar_minfwd
 414       write(6,*)'madis_metar_recwin           = ',madis_metar_recwin
 415       write(6,*)'madis_metar_all_providers    = ',madis_metar_all_providers
 416       write(6,*)'madis_metar_provider_list    = ',trim(madis_metar_provider_list)
 417       write(6,*)
 418       write(6,*)'madis_sound_db               = ',madis_sound_db
 419       write(6,*)'madis_sound_minbck           = ',madis_sound_minbck
 420       write(6,*)'madis_sound_minfwd           = ',madis_sound_minfwd
 421       write(6,*)'madis_sound_recwin           = ',madis_sound_recwin
 422       write(6,*)
 423       write(6,*)'madis_profiler_db            = ',madis_profiler_db
 424       write(6,*)'madis_profiler_minbck        = ',madis_profiler_minbck
 425       write(6,*)'madis_profiler_minfwd        = ',madis_profiler_minfwd
 426       write(6,*)'madis_profiler_recwin        = ',madis_profiler_recwin
 427       write(6,*)'madis_profiler_all_providers = ',madis_profiler_all_providers
 428       write(6,*)'madis_profiler_provider_list = ',trim(madis_profiler_provider_list)
 429       write(6,*)
 430       write(6,*)'madis_airep_db               = ',madis_airep_db
 431       write(6,*)'madis_airep_minbck           = ',madis_airep_minbck
 432       write(6,*)'madis_airep_minfwd           = ',madis_airep_minfwd
 433       write(6,*)'madis_airep_recwin           = ',madis_airep_recwin
 434       write(6,*)
 435       write(6,*)'madis_geoamv_db              = ',madis_geoamv_db
 436       write(6,*)'madis_geoamv_minbck          = ',madis_geoamv_minbck
 437       write(6,*)'madis_geoamv_minfwd          = ',madis_geoamv_minfwd
 438       write(6,*)'madis_geoamv_recwin          = ',madis_geoamv_recwin
 439       write(6,*)'madis_geoamv_all_providers   = ',madis_geoamv_all_providers
 440       write(6,*)'madis_geoamv_provider_list   = ',trim(madis_geoamv_provider_list)
 441       write(6,*)'madis_geoamv_wv              = ',madis_geoamv_wv
 442       write(6,*)'madis_geoamv_ir              = ',madis_geoamv_ir
 443       write(6,*)'madis_geoamv_vis             = ',madis_geoamv_vis
 444       write(6,*)'madis_geoamv_s10             = ',madis_geoamv_s10
 445       write(6,*)'madis_geoamv_s11             = ',madis_geoamv_s11
 446
 447    endif
 448
 449    write(6,*)
 450
 451 ! Allocate the local memory needed to read in the MADIS arrays that are
 452 ! dimensioned by station.
 453
 454    write(6,*) ' Allocating memory for MADIS arrays...'
 455
 456    allocate (stanam(max_sta))
 457    allocate (provdr(max_sta))
 458    allocate (timeob(max_sta))
 459    allocate (wmoid(max_sta))
 460    allocate (lat(max_sta))
 461    allocate (lon(max_sta))
 462    allocate (elev(max_sta))
 463    allocate (pw(max_sta))
 464    allocate (slp(max_sta))
 465    allocate (keep(max_sta))
 466    allocate (nlevels(max_sta))
 467
 468 ! Initialize the dataset ob counts.
 469
 470    do i = 1, num_ob_indexes
 471       nlocal(i) = 0
 472       ntotal(i) = 0
 473    enddo
 474
 475    write(6,*) ' Memory allocated, obcounts set to zero'
 476
 477 ! Initialize the MADIS subsystems to be used, using the database
 478 ! selected for each subsystem by the user, and with error messages
 479 ! echoed to standard output.  For subsystems that have optional
 480 ! data provider selection files, read those in as required, and
 481 ! set the subset of data providers to be used.
 482
 483    pvdriostat = 0                  ! Indicates I/O errors on provider files.
 484
 485    if(use_metarobs .or. use_shipsobs .or. use_gpspwobs)then
 486
 487 ! Because these three WRF datasets map into only one MADIS subsystem, and
 488 ! the user has database & time window control over each of them independently,
 489 ! we can only initialize one of them here.  We'll initialize the first one
 490 ! that will be used below when we read in the data.  (This is because something
 491 ! has to be initialized so we can get proper error message handling inside
 492 ! the MADIS routines.)
 493
 494       if(use_metarobs)then
 495
 496          CALL MINIT('SFC',madis_metar_db,.true.,istatus)
 497          if(istatus /= success_p)stop
 498          CALL MSETSFCPVDR('ALL-SFC',.false.,istatus)
 499          if(istatus /= success_p)stop
 500
 501          if(.not. madis_metar_all_providers)then
 502
 503             if(print_detail_obs)then
 504                write(6,*)
 505                write(6,*)'Selecting individual METAR providers from ',&
 506                          trim(madis_metar_provider_list)
 507             endif
 508             if(istatus /= success_p)stop
 509
 510             pvdriotype = 'Open'
 511             pvdrfile = madis_metar_provider_list
 512             open(unit=lun,file=trim(pvdrfile),status='OLD', &
 513                  iostat=pvdriostat,err=8000)
 514
 515             pvdriotype = 'Read'
 516
 517             do i = 1, max_provider
 518                read(lun,1,err=8000,iostat=pvdriostat,end=2)dpname
 519  1             format(a)
 520                if(.not. (dpname(1:1) == '#' .or. len(trim(dpname)) == 0 .or. &
 521                          dpname(1:8) == 'MARITIME'))then
 522                   if(print_detail_obs)write(6,*)'  Selecting:  ',trim(dpname)
 523                   CALL MSETSFCPVDR(trim(dpname),.true.,istatus)
 524                   if(istatus /= success_p)stop
 525                endif
 526             enddo
 527
 528  2          pvdriostat = 0
 529             close(lun)
 530
 531          else
 532
 533 ! Select all of the surface datasets other than maritime.
 534
 535             CALL MSETSFCPVDR('ALL-MTR',.true.,istatus)
 536             if(istatus /= success_p)stop
 537             CALL MSETSFCPVDR('SAO',.true.,istatus)
 538             if(istatus /= success_p)stop
 539             CALL MSETSFCPVDR('ALL-MESO',.true.,istatus)
 540             if(istatus /= success_p)stop
 541             CALL MSETSFCPVDR('COOP',.true.,istatus)
 542             if(istatus /= success_p)stop
 543
 544          endif
 545
 546       else if(use_shipsobs)then
 547
 548          CALL MINIT('SFC',madis_ships_db,.true.,istatus)
 549          if(istatus /= success_p)stop
 550
 551       else if(use_gpspwobs)then
 552
 553          CALL MINIT('SFC',madis_gpspw_db,.true.,istatus)
 554          if(istatus /= success_p)stop
 555
 556       endif
 557
 558    endif
 559
 560    if(use_soundobs)then
 561
 562       CALL MINIT('RAOB',madis_sound_db,.true.,istatus)
 563       if(istatus /= success_p)stop
 564
 565    endif
 566
 567 ! ***** Remember to add use_xxxobs logic for profilers once that's added.
 568
 569    if(use_soundobs)then   ! ***** change to use_profilerobs...
 570
 571       map_init = .false.
 572       npn_init = .false.
 573
 574       if(.not. madis_profiler_all_providers)then
 575
 576          if(print_detail_obs)then
 577             write(6,*)
 578             write(6,*)'Selecting individual PROFILER providers from ',&
 579                          trim(madis_profiler_provider_list)
 580          endif
 581
 582          pvdriotype = 'Open'
 583          pvdrfile = madis_profiler_provider_list
 584          open(unit=lun,file=trim(pvdrfile),status='OLD', &
 585               iostat=pvdriostat,err=8000)
 586
 587          pvdriotype = 'Read'
 588
 589          do i = 1, max_provider
 590
 591 ! Because the WRF profiler dataset maps into two MADIS subsystems (NPN and MAP)
 592 ! we need to make sure that we initialize each subsystem as it's encountered
 593 ! in the user's list.
 594
 595             read(lun,1,err=8000,iostat=pvdriostat,end=3)dpname
 596             if(.not. (dpname(1:1) == '#' .or. len(trim(dpname)) == 0))then
 597
 598                if(print_detail_obs)write(6,*)'  Selecting:  ',trim(dpname)
 599
 600                if(trim(dpname) == 'NOAA')then
 601
 602                   CALL MINIT('NPN',madis_profiler_db,.true.,istatus)
 603                   if(istatus /= success_p)stop
 604                   npn_init = .true.
 605
 606                else if(trim(dpname) == 'ALL-PROF')then
 607
 608                   CALL MINIT('NPN',madis_profiler_db,.true.,istatus)
 609                   if(istatus /= success_p)stop
 610                   npn_init = .true.
 611                   CALL MINIT('MAP','FSL',.true.,istatus)   ! No AWIPS database for MAP.
 612                   if(istatus /= success_p)then
 613
 614 ! We'll be forgiving if the AWIPS profiler database has been selected, and just
 615 ! set things so no MAP data will be read.  If MADIS (FSL) has been selected,
 616 ! however, this is obviously an error on the user's part, and he/she should fix it.
 617
 618                      if(madis_profiler_db(1:3).eq.'FSL')stop
 619                      map_init = .false.
 620                   else
 621                      map_init = .true.
 622                   endif
 623
 624                else if(trim(dpname) == 'ALL-MAP')then
 625
 626                   CALL MINIT('MAP','FSL',.true.,istatus)
 627                   if(istatus /= success_p)then
 628                      if(madis_profiler_db(1:3).eq.'FSL')stop
 629                      map_init = .false.
 630                   else
 631                      map_init = .true.
 632                   endif
 633
 634 ! ***** Get rid of the NO-PROF else once profiler data structures have
 635 ! been added.
 636
 637                else if(trim(dpname) == 'NO-PROF')then
 638
 639                   continue
 640
 641                else
 642
 643 ! An individual MAP provider is being selected.  Make sure we're initialized
 644 ! for MAP, then add in the provider.
 645
 646                   if(.not. map_init)then
 647                      CALL MINIT('MAP','FSL',.true.,istatus)
 648                      if(istatus /= success_p)then
 649                         if(madis_profiler_db(1:3).eq.'FSL')stop
 650                         map_init = .false.
 651                      else
 652                         map_init = .true.
 653                         CALL MSETMAPPVDR('ALL-MAP',.false.,istatus)
 654                         if(istatus /= success_p)stop
 655                      endif
 656                   endif
 657
 658                   if(map_init)then
 659                      CALL MSETMAPPVDR(trim(dpname),.true.,istatus)
 660                      if(istatus /= success_p)stop
 661                   endif
 662
 663                endif
 664
 665             endif
 666
 667          enddo
 668
 669  3       pvdriostat = 0
 670          close(lun)
 671
 672       else
 673
 674 ! "All profilers" consists of both the NPN and the MAP subsystems.  We'll be
 675 ! forgiving if an AWIPS database user naively included MAP (by not specifically
 676 ! setting the providers).  A user of the MADIS database, however, should know
 677 ! better -- we'll stop with a fatal exit.
 678
 679          CALL MINIT('NPN',madis_profiler_db,.true.,istatus)
 680          if(istatus /= success_p)stop
 681          npn_init = .true.
 682          CALL MINIT('MAP','FSL',.true.,istatus) ! No AWIPS database for MAP.
 683          if(istatus /= success_p)then
 684             if(madis_profiler_db(1:3).eq.'FSL')stop
 685             map_init = .false.
 686          else
 687             map_init = .true.
 688          endif
 689
 690       endif
 691
 692    endif
 693
 694    if(use_airepobs)then
 695
 696       CALL MINIT('ACARS',madis_airep_db,.true.,istatus)
 697       if(istatus /= success_p)stop
 698
 699    endif
 700
 701    if(use_geoamvobs)then
 702
 703       CALL MINIT('SATWND',madis_geoamv_db,.true.,istatus)
 704       if(istatus /= success_p)stop
 705
 706       if(.not. madis_geoamv_all_providers)then
 707
 708          if(print_detail_obs)then
 709             write(6,*)
 710             write(6,*)'Selecting individual SATOB providers from ',&
 711                        trim(madis_geoamv_provider_list)
 712          endif
 713          CALL MSETSATWNDPVDR('ALL-GOES-O',.false.,istatus)
 714          if(istatus /= success_p)stop
 715
 716          pvdriotype = 'Open'
 717          pvdrfile = madis_geoamv_provider_list
 718          open(unit=lun,file=trim(pvdrfile),status='OLD', &
 719               iostat=pvdriostat,err=8000)
 720
 721          pvdriotype = 'Read'
 722
 723          do i = 1, max_provider
 724             read(lun,1,err=8000,iostat=pvdriostat,end=4)dpname
 725             if(.not. (dpname(1:1) == '#' .or. len(trim(dpname)) == 0))then
 726                if(print_detail_obs)write(6,*)'  Selecting:  ',trim(dpname)
 727                CALL MSETSATWNDPVDR(trim(dpname),.true.,istatus)
 728                if(istatus /= success_p)stop
 729             endif
 730          enddo
 731
 732  4       pvdriostat = 0
 733          close(lun)
 734
 735       endif
 736
 737 ! Set the number and types of SATW variables to read.
 738
 739       nsatw_vartype = 0
 740
 741       if(madis_geoamv_ir.and.madis_geoamv_wv.and.madis_geoamv_vis.and. &
 742          madis_geoamv_s10.and.madis_geoamv_s11)then
 743
 744          nsatw_vartype = 1
 745          satw_var_len(1) = 2
 746          satw_var(1,1) = 'DD'
 747          satw_var(2,1) = 'FF'
 748
 749       else
 750
 751          if(madis_geoamv_ir)then
 752             nsatw_vartype = nsatw_vartype + 1
 753             satw_var_len(nsatw_vartype) = 4
 754             satw_var(1,nsatw_vartype) = 'DDIR'
 755             satw_var(2,nsatw_vartype) = 'FFIR'
 756          endif
 757
 758          if(madis_geoamv_wv)then
 759             nsatw_vartype = nsatw_vartype + 1
 760             satw_var_len(nsatw_vartype) = 4
 761             satw_var(1,nsatw_vartype) = 'DDWV'
 762             satw_var(2,nsatw_vartype) = 'FFWV'
 763          endif
 764
 765          if(madis_geoamv_vis)then
 766             nsatw_vartype = nsatw_vartype + 1
 767             satw_var_len(nsatw_vartype) = 5
 768             satw_var(1,nsatw_vartype) = 'DDVIS'
 769             satw_var(2,nsatw_vartype) = 'FFVIS'
 770          endif
 771
 772          if(madis_geoamv_s10)then
 773             nsatw_vartype = nsatw_vartype + 1
 774             satw_var_len(nsatw_vartype) = 5
 775             satw_var(1,nsatw_vartype) = 'DDS10'
 776             satw_var(2,nsatw_vartype) = 'FFS10'
 777          endif
 778
 779          if(madis_geoamv_s11)then
 780             nsatw_vartype = nsatw_vartype + 1
 781             satw_var_len(nsatw_vartype) = 5
 782             satw_var(1,nsatw_vartype) = 'DDS11'
 783             satw_var(2,nsatw_vartype) = 'FFS11'
 784          endif
 785
 786       endif
 787
 788    endif
 789
 790 ! Select the QC level to be returned to "max possible QC level for each
 791 ! variable".
 792
 793    CALL MSETQC(qcall,istatus)
 794    if(istatus /= success_p)stop
 795
 796 ! Convert the user's desired time to the format used by MADIS
 797 ! (YYYYMMDD_HHMM).
 798
 799    tstr = analysis_date(1:4)//analysis_date(6:7)//analysis_date(9:10)// &
 800           '_'//analysis_date(12:13)//analysis_date(15:16)
 801
 802 ! ***** This commented-out section is used only by MADIS developers to
 803 !       hardwire the desired time when doing stand-alone testing.
 804 !
 805 !   open(unit=19,file='mfb_time.txt',status='old')
 806 !   read(19,*)tstr
 807 !   close(unit=19)
 808 !
 809 ! ***** End of test code section.
 810
 811 ! Now convert 'YYYYMMDD_HH00' to 'YYJJJHH00' (JJJ = Julian date).
 812
 813    CALL MTRNTIM(tstr,1,atime,istatus)
 814    if(istatus /= success_p)stop
 815
 816 !------------------------------------------------------------------------------
 817 ! MADIS surface datasets --> WRF datasets METAR, GPSPW, SHIPS
 818 !------------------------------------------------------------------------------
 819 ! If we're not processing any of these WRF datasets, skip to the next dataset.
 820
 821    if(.not. use_metarobs .and. .not. use_shipsobs .and. .not. use_gpspwobs)go to 150
 822
 823 ! Allocate the local memory for the MADIS arrays that have a level dimension
 824 ! (only 1 level in this case).
 825
 826    allocate (p(max_sfc))
 827    allocate (q(max_sfc))
 828    allocate (t(max_sfc))
 829    allocate (z(max_sfc))
 830    allocate (dd(max_sfc))
 831    allocate (ff(max_sfc))
 832    allocate (qca(max_sfc))
 833    allocate (qcr(max_sfc))
 834    allocate (qcd(max_sfc))
 835
 836    if(use_metarobs)then
 837
 838 ! Set the time window for this dataset with the user's choices.
 839
 840       CALL MSETWIN(madis_metar_minbck,madis_metar_minfwd,madis_metar_recwin,istatus)
 841       if(istatus /= success_p)stop
 842
 843 ! Load in the stations that have data in this time window.  If we have any
 844 ! problems reading in the stations we'll skip this dataset.  When we're done
 845 ! with all the datasets, if nothing's been read in there's a problem
 846 ! and we'll stop with an error.  If a particular dataset isn't there, but
 847 ! the others aren't, that's not as fundamental of a problem.
 848
 849       CALL MSFCSTA(atime,nsta,stanam,wmoid,lat,lon,elev,timeob,provdr,istatus)
 850       if(istatus /= success_p)go to 50
 851
 852 ! Read in the observations, with QC filtering as specified in MSETQC --
 853 ! this will put MADIS missing flags in all obs that didn't pass the QC.
 854 ! Note that we don't check the status on each variable read.  It's simpler
 855 ! to fill in missing flags into each variable before reading them, then
 856 ! we can sort this all out below.
 857
 858       p(1:nsta) = miss_p
 859       CALL MGETSFC(atime,'P',ntmp,nobs,p,qcd,qca,qcr,istatus)
 860       q(1:nsta) = miss_p
 861       CALL MGETSFC(atime,'Q',ntmp,nobs,q,qcd,qca,qcr,istatus)
 862       t(1:nsta) = miss_p
 863       CALL MGETSFC(atime,'T',ntmp,nobs,t,qcd,qca,qcr,istatus)
 864       dd(1:nsta) = miss_p
 865       CALL MGETSFC(atime,'DD',ntmp,nobs,dd,qcd,qca,qcr,istatus)
 866       ff(1:nsta) = miss_p
 867       CALL MGETSFC(atime,'FF',ntmp,nobs,ff,qcd,qca,qcr,istatus)
 868       slp(1:nsta) = miss_p
 869       CALL MGETSFC(atime,'SLP',ntmp,nobs,slp,qcd,qca,qcr,istatus)
 870
 871 ! Fill in missing flags for the PW variable, copy ELEV to Z
 872 ! and fill in 1 level for each station for use by the PROCESS_DATASET
 873 ! routine.  Also, any stations at sea level that have SLP but not P,
 874 ! put SLP into P -- while a useful idea in general, this particularly
 875 ! makes a difference with maritime obs -- usually at sea level, with
 876 ! a lot of stations measuring SLP but not P.
 877
 878       do i = 1, nsta
 879          pw(i) = miss_p
 880          z(i) = elev(i)
 881          nlevels(i) = 1
 882          if(p(i) == miss_p)then
 883             if(elev(i) == 0.0 .and. slp(i) /= miss_p)p(i) = slp(i)
 884          endif
 885       enddo
 886
 887 ! Determine which ones are inside the user's domain, and then allocate
 888 ! memory for the IV data for this dataset.  Also note that we toss surface obs
 889 ! that don't have P.
 890
 891       ntmp = 0                  ! Count of records inside domain.
 892
 893       do i = 1, nsta
 894
 895 ! Put LAT/LON into PLATFORM and call a DA subroutine to determine if
 896 ! the station is in the domain.  Also check for bad (mesonet) station
 897 ! locations -- we have at least one provider who puts these into his
 898 ! station table when he doesn't have valid coordinates.
 899
 900          if(abs(lat(i)).lt.90.and.abs(lon(i)).le.180)then
 901             platform % info % lat = lat(i)
 902             platform % info % lon = lon(i)
 903
 904             CALL da_llxy(platform%info,platform%loc,outside,outside_all)
 905
 906             if (outside_all .or. p(i) == miss_p)then
 907
 908                 keep(i) = 0         ! 0 = Not in domain or no P.
 909
 910             else
 911
 912                 ntotal(metar) = ntotal(metar) + 1
 913                 if(outside)then
 914
 915                    keep(i) = 0
 916
 917                 else
 918
 919                    keep(i) = 1         ! 1 = Process.
 920                    ntmp = ntmp + 1
 921
 922                    if (global .and. (platform%loc%i < ids.or.platform%loc%i >= ide)) then
 923                       nlocal(metar) = nlocal(metar) + 2
 924                    else
 925                       nlocal(metar) = nlocal(metar) + 1
 926                    endif
 927
 928                 endif
 929
 930              endif
 931
 932          else
 933
 934             keep(i) = 0
 935
 936          endif
 937
 938       enddo
 939
 940       if(nlocal(metar) > 0)then
 941
 942          allocate (iv % metar(1:nlocal(metar)))
 943
 944 ! Call all land surface data METAR, and tell the sub that there's one level.
 945
 946          fm = 15
 947          nlev = 1
 948          iv%info(metar)%max_lev = 1
 949
 950 ! Move the data out of the ingest arrays into the IV structure,
 951 ! including attendant processing that fills in the location within
 952 ! the user's domain, etc.
 953
 954          CALL PROCESS_DATASET(nsta,nlev,fm,keep,lat,lon,elev,timeob, &
 955                               slp,pw,stanam,p,q,t,dd,ff,z,nlocal,ntotal, &
 956                               0,0,nlevels,iv,nlocal_out,ntotal_out,sfclev)
 957
 958 ! Note that if PROCESS_DATASET detects that there's no useful data for
 959 ! a station, it won't put it into IV.  Therefore we'll adjust NLOCAL
 960 ! accordingly (FILL_IV_INFO adjusts NTOTAL in the structure where it's
 961 ! sent back to the data assimilation code).
 962
 963          iv%info(metar)%nlocal = nlocal_out
 964
 965       endif
 966
 967    endif
 968
 969  50 if(use_shipsobs)then
 970
 971 ! Because this WRF dataset is a subset of the MADIS surface subsystem,
 972 ! we have to re-initialize the MADIS subsystem here, with provider selection,
 973 ! QC level selection, etc.
 974
 975       CALL MINIT('SFC',madis_ships_db,.true.,istatus)
 976       if(istatus /= success_p)stop
 977       CALL MSETSFCPVDR('ALL-SFC',.false.,istatus)
 978       if(istatus /= success_p)stop
 979       CALL MSETSFCPVDR('MARITIME',.true.,istatus)
 980       if(istatus /= success_p)stop
 981       CALL MSETQC(qcall,istatus)
 982       if(istatus /= success_p)stop
 983
 984 ! Set the time window for this dataset with the user's choices.
 985
 986       CALL MSETWIN(madis_ships_minbck,madis_ships_minfwd,madis_ships_recwin,istatus)
 987       if(istatus /= success_p)stop
 988
 989 ! Load in the stations that have data in this time window.
 990
 991       CALL MSFCSTA(atime,nsta,stanam,wmoid,lat,lon,elev,timeob,provdr,istatus)
 992       if(istatus /= success_p)go to 75
 993
 994 ! Read in the observations.
 995
 996       p(1:nsta) = miss_p
 997       CALL MGETSFC(atime,'P',ntmp,nobs,p,qcd,qca,qcr,istatus)
 998       q(1:nsta) = miss_p
 999       CALL MGETSFC(atime,'Q',ntmp,nobs,q,qcd,qca,qcr,istatus)
1000       t(1:nsta) = miss_p
1001       CALL MGETSFC(atime,'T',ntmp,nobs,t,qcd,qca,qcr,istatus)
1002       dd(1:nsta) = miss_p
1003       CALL MGETSFC(atime,'DD',ntmp,nobs,dd,qcd,qca,qcr,istatus)
1004       ff(1:nsta) = miss_p
1005       CALL MGETSFC(atime,'FF',ntmp,nobs,ff,qcd,qca,qcr,istatus)
1006       slp(1:nsta) = miss_p
1007       CALL MGETSFC(atime,'SLP',ntmp,nobs,slp,qcd,qca,qcr,istatus)
1008
1009 ! Fill in missing flags for the PW variable, copy ELEV to Z
1010 ! and fill in 1 level for each station for use by the PROCESS_DATASET
1011 ! routine.  Also, any stations at sea level that have SLP but not P,
1012 ! put SLP into P -- while a useful idea in general, this particularly
1013 ! makes a difference with maritime obs -- usually at sea level, with
1014 ! a lot of stations measuring SLP but not P.
1015
1016       do i = 1, nsta
1017          pw(i) = miss_p
1018          z(i) = elev(i)
1019          nlevels(i) = 1
1020          if(p(i) == miss_p)then
1021             if(elev(i) == 0.0 .and. slp(i) /= miss_p)p(i) = slp(i)
1022          endif
1023       enddo
1024
1025 ! Determine which ones are inside the user's domain, and then allocate
1026 ! memory for the IV data for this dataset.  Also note that we toss ships obs
1027 ! that don't have P.
1028
1029       ntmp = 0                  ! Count of records inside domain.
1030
1031       do i = 1, nsta
1032
1033 ! Put LAT/LON into PLATFORM and call a DA subroutine to determine if
1034 ! the station is in the domain.
1035
1036          if(abs(lat(i)).lt.90.and.abs(lon(i)).le.180)then
1037             platform % info % lat = lat(i)
1038             platform % info % lon = lon(i)
1039             CALL da_llxy(platform%info,platform%loc,outside,outside_all)
1040
1041             if (outside_all .or. p(i) == miss_p)then
1042                keep(i) = 0
1043             else
1044                 ntotal(ships) = ntotal(ships) + 1
1045                 if(outside)then
1046                    keep(i) = 0
1047                 else
1048                    keep(i) = 1
1049                    ntmp = ntmp + 1
1050                    if (global .and. (platform%loc%i < ids .or. platform%loc%i >= ide)) then
1051                       nlocal(ships) = nlocal(ships) + 2
1052                    else
1053                       nlocal(ships) = nlocal(ships) + 1
1054                    endif
1055                 endif
1056             endif
1057          else
1058             keep(i) = 0
1059          endif
1060
1061       enddo
1062
1063       if(nlocal(ships) > 0)then
1064
1065          allocate (iv % ships(1:nlocal(ships)))
1066          fm = 13
1067          iv%info(ships)%max_lev = 1
1068          CALL PROCESS_DATASET(nsta,nlev,fm,keep,lat,lon,elev,timeob, &
1069                               slp,pw,stanam,p,q,t,dd,ff,z,nlocal,ntotal, &
1070                               0,0,nlevels,iv,nlocal_out,ntotal_out,sfclev)
1071          iv%info(ships)%nlocal = nlocal_out
1072       endif
1073
1074    endif
1075
1076  75 if(use_gpspwobs)then
1077
1078 ! Reinitialize the MADIS surface subsystem for GPSPW.
1079
1080       CALL MINIT('SFC',madis_gpspw_db,.true.,istatus)
1081       if(istatus /= success_p)stop
1082       CALL MSETSFCPVDR('ALL-SFC',.false.,istatus)
1083       if(istatus /= success_p)stop
1084       CALL MSETSFCPVDR('GPSMET',.true.,istatus)
1085       if(istatus /= success_p)stop
1086       CALL MSETQC(qcall,istatus)
1087       if(istatus /= success_p)stop
1088
1089 ! Set the time window for this dataset with the user's choices.
1090
1091       CALL MSETWIN(madis_gpspw_minbck,madis_gpspw_minfwd,madis_gpspw_recwin,istatus)
1092       if(istatus /= success_p)stop
1093
1094 ! Load in the stations that have data in this time window.
1095
1096       CALL MSFCSTA(atime,nsta,stanam,wmoid,lat,lon,elev,timeob,provdr,istatus)
1097       if(istatus /= success_p)go to 100
1098
1099 ! Read in the observations.  Note that we put missing flags into DD/FF so
1100 ! that the DD/FF->U/V conversion won't be attempted in PROCESS_DATASETS.
1101
1102       dd(1:nsta) = miss_p
1103       ff(1:nsta) = miss_p
1104
1105       pw(1:nsta) = miss_p
1106       CALL MGETSFC(atime,'PWV',ntmp,nobs,pw,qcd,qca,qcr,istatus)
1107       slp(1:nsta) = miss_p
1108       CALL MGETSFC(atime,'SLP',ntmp,nobs,slp,qcd,qca,qcr,istatus)
1109
1110 ! Copy ELEV to Z and fill in 1 level for each station for use by the
1111 ! PROCESS_DATASET routine.
1112
1113       do i = 1, nsta
1114          z(i) = elev(i)
1115          nlevels(i) = 1
1116       enddo
1117
1118 ! Determine which ones are inside the user's domain and have PW, and then
1119 ! allocate memory for the IV data for this dataset.
1120
1121       ntmp = 0
1122
1123       do i = 1, nsta
1124
1125          if(abs(lat(i)).lt.90.and.abs(lon(i)).le.180)then
1126             platform % info % lat = lat(i)
1127             platform % info % lon = lon(i)
1128             CALL da_llxy(platform%info,platform%loc,outside,outside_all)
1129
1130             if (outside_all .or. pw(i) == miss_p)then
1131                keep(i) = 0
1132             else
1133                 ntotal(gpspw) = ntotal(gpspw) + 1
1134                 if(outside)then
1135                    keep(i) = 0
1136                 else
1137                    keep(i) = 1
1138                    pw(i) = pw(i) * 100.0 ! Convert from m to cm.
1139                    ntmp = ntmp + 1
1140                    if (global .and. (platform%loc%i < ids .or. platform%loc%i >= ide)) then
1141                       nlocal(gpspw) = nlocal(gpspw) + 2
1142                    else
1143                       nlocal(gpspw) = nlocal(gpspw) + 1
1144                    endif
1145                 endif
1146             endif
1147          else
1148             keep = 0
1149          endif
1150
1151       enddo
1152
1153       if(nlocal(gpspw) > 0)then
1154
1155          allocate (iv % gpspw(1:nlocal(gpspw)))
1156          fm = 111
1157          iv%info(gpspw)%max_lev = 1
1158          CALL PROCESS_DATASET(nsta,nlev,fm,keep,lat,lon,elev,timeob, &
1159                               slp,pw,stanam,p,q,t,dd,ff,z,nlocal,ntotal, &
1160                               0,0,nlevels,iv,nlocal_out,ntotal_out,sfclev)
1161          iv%info(gpspw)%nlocal = nlocal_out
1162       endif
1163
1164    endif
1165
1166 ! Finish up with the MADIS surface datasets by deallocating the
1167 ! levels-dimensioned arrays.
1168
1169  100  deallocate (p)
1170    deallocate (q)
1171    deallocate (t)
1172    deallocate (z)
1173    deallocate (dd)
1174    deallocate (ff)
1175    deallocate (qca)
1176    deallocate (qcr)
1177    deallocate (qcd)
1178
1179 !------------------------------------------------------------------------------
1180 ! MADIS ACARS dataset --> WRF AIREP dataset
1181 !------------------------------------------------------------------------------
1182 ! If we're not using this WRF dataset, skip to the next one.
1183
1184  150  if(.not. use_airepobs)go to 250
1185
1186 ! Set the time window for this dataset with the user's choices.
1187
1188    CALL MSETWIN(madis_airep_minbck,madis_airep_minfwd,madis_airep_recwin,istatus)
1189    if(istatus /= success_p)stop
1190
1191 ! Load in the stations with data for this time window.
1192
1193    CALL MACARSSTA(atime,nsta,stanam,timeob,istatus)
1194    if(istatus /= success_p)go to 200
1195
1196 ! Allocate the local memory for the MADIS arrays that have a level dimension.
1197
1198    allocate (p(nsta))
1199    allocate (q(nsta))
1200    allocate (t(nsta))
1201    allocate (z(nsta))
1202    allocate (dd(nsta))
1203    allocate (ff(nsta))
1204    allocate (qca(nsta))
1205    allocate (qcr(nsta))
1206    allocate (qcd(nsta))
1207
1208 ! Read in the observations.  Note that unlike the other MADIS datasets,
1209 ! the lat/lon are read in as obs, not station info.
1210
1211    p(1:nsta) = miss_p
1212    CALL MGETACARS(atime,'P',ntmp,nobs,p,qcd,qca,qcr,istatus)
1213    t(1:nsta) = miss_p
1214    CALL MGETACARS(atime,'T',ntmp,nobs,t,qcd,qca,qcr,istatus)
1215    dd(1:nsta) = miss_p
1216    CALL MGETACARS(atime,'DD',ntmp,nobs,dd,qcd,qca,qcr,istatus)
1217    ff(1:nsta) = miss_p
1218    CALL MGETACARS(atime,'FF',ntmp,nobs,ff,qcd,qca,qcr,istatus)
1219    lat(1:nsta) = miss_p
1220    CALL MGETACARS(atime,'LAT',ntmp,nobs,lat,qcd,qca,qcr,istatus)
1221    lon(1:nsta) = miss_p
1222    CALL MGETACARS(atime,'LON',ntmp,nobs,lon,qcd,qca,qcr,istatus)
1223    z(1:nsta) = miss_p
1224    CALL MGETACARS(atime,'HT',ntmp,nobs,z,qcd,qca,qcr,istatus)
1225
1226 ! Fill in missing flags for the PW, SLP and ELEV variables, and fill in
1227 ! 1 level for each station for use by the PROCESS_DATASET routine.
1228
1229    do i = 1, nsta
1230       pw(i) = miss_p
1231       slp(i) = miss_p
1232       elev(i) = miss_p
1233       nlevels(i) = 1
1234    enddo
1235
1236 ! Determine which ones are inside the user's domain, and then allocate
1237 ! memory for the IV data for this dataset.
1238
1239    ntmp = 0                      ! Count of records inside domain.
1240
1241    do i = 1, nsta
1242
1243       if (abs(lat(i)).lt.90.and.abs(lon(i)).le.180) then
1244          platform % info % lat = lat(i)
1245          platform % info % lon = lon(i)
1246          CALL da_llxy(platform%info,platform%loc,outside,outside_all)
1247
1248          if (outside_all) then
1249
1250             keep(i) = 0         ! 0 = Not in domain.
1251
1252          else
1253
1254             ntotal(airep) = ntotal(airep) + 1
1255             if (outside) then
1256                keep(i) = 0
1257             else
1258                keep(i) = 1         ! 1 = Process.
1259                ntmp = ntmp + 1
1260                if (global .and. (platform%loc%i < ids .or. platform%loc%i >= ide)) then
1261                   nlocal(airep) = nlocal(airep) + 2
1262                else
1263                   nlocal(airep) = nlocal(airep) + 1
1264                endif
1265             endif
1266
1267          endif
1268
1269       else
1270
1271          keep(i) = 0
1272
1273       endif
1274
1275    enddo
1276
1277    if(nlocal(airep) > 0)then
1278       allocate (iv % airep(1:nlocal(airep)))
1279    endif
1280
1281    fm = 96
1282    nlev = 1
1283    iv%info(airep)%max_lev = 1
1284
1285 ! Move the data out of the ingest arrays into the IV structure.
1286
1287    CALL PROCESS_DATASET(nsta,nlev,fm,keep,lat,lon,elev,timeob, &
1288                         slp,pw,stanam,p,q,t,dd,ff,z,nlocal,ntotal, &
1289                         0,0,nlevels,iv,nlocal_out,ntotal_out,sfclev)
1290
1291    iv%info(airep)%nlocal = nlocal_out
1292    iv%info(airep)%plocal(1) = nlocal_out
1293    iv%info(airep)%ptotal(1) = ntotal(airep)
1294    iv%info(airep)%ntotal = ntotal(airep)
1295
1296
1297 ! Finish up with the MADIS ACARS datasets by deallocating the
1298 ! levels-dimensioned arrays.
1299
1300  200  deallocate (p)
1301    deallocate (q)
1302    deallocate (t)
1303    deallocate (z)
1304    deallocate (dd)
1305    deallocate (ff)
1306    deallocate (qca)
1307    deallocate (qcr)
1308    deallocate (qcd)
1309
1310 !------------------------------------------------------------------------------
1311 ! MADIS RAOB/NPN/MAP datasets --> WRF SOUND dataset
1312 !------------------------------------------------------------------------------
1313
1314  250  if(.not. use_soundobs)go to 300
1315
1316 ! ***** Split off the RAOB dataset from this section once profiler data
1317 ! structures have been added.
1318
1319 ! ***** Remember to add the use_xxxobs logic for profilers.
1320
1321 ! This MADIS-->WRF dataset transition is more complicated than the other
1322 ! ones we've done above, because we need to map multiple MADIS datasets
1323 ! into one WRF dataset, and we need to know how many reports we're going
1324 ! to have in the domain in order to allocate the IV memory.  Therefore,
1325 ! we'll set the time window for each MADIS dataset, read its station info,
1326 ! and count up the reports inside the domain over all the MADIS datasets.
1327 ! Once that's been done, then the same tasks will be done as for the
1328 ! simpler datasets above.
1329
1330 ! RAOB
1331
1332 ! Set the time window for this dataset with the user's choices.
1333
1334    CALL MSETWIN(madis_sound_minbck,madis_sound_minfwd,madis_sound_recwin,istatus)
1335    if(istatus /= success_p)stop
1336    ntmp_r = 0
1337    CALL MRAOBSTA(atime,nsta,stanam,wmoid,lat,lon,elev,timeob,istatus)
1338    if(istatus == success_p)then
1339       do i = 1, nsta
1340          if(abs(lat(i)).lt.90.and.abs(lon(i)).le.180)then
1341             platform % info % lat = lat(i)
1342             platform % info % lon = lon(i)
1343             CALL da_llxy(platform%info,platform%loc,outside,outside_all)
1344             if (outside_all) then
1345                keep_r(i) = 0    ! 0 = Not in domain.
1346             else
1347                ntotal(sound) = ntotal(sound) + 1
1348                if(outside)then
1349                    keep(i) = 0
1350                else
1351                    keep_r(i) = 1    ! 1 = Process.
1352                    ntmp_r = ntmp_r + 1
1353                    if (global .and. (platform%loc%i < ids .or. platform%loc%i >= ide)) then
1354                       nlocal(sound) = nlocal(sound) + 2
1355                    else
1356                       nlocal(sound) = nlocal(sound) + 1
1357                    endif
1358                endif
1359             endif
1360          else
1361             keep_r(i) = 0
1362          endif
1363
1364       enddo
1365    endif
1366
1367 ! Keep track of nlocal and notal for use by PROCESS_DATASET when it needs to
1368 ! know these values over several calls.
1369
1370    nl_raob = nlocal(sound)
1371    nt_raob = ntotal(sound)
1372
1373 ! NPN
1374
1375    write(unit=message(1),fmt='(A,I10)')'ntotal_raob=',ntotal(sound)
1376    call da_message(message(1:1))
1377    ntmp_n = 0
1378
1379    if(npn_init)then
1380
1381 ! Set the time window for this dataset with the user's choices.
1382
1383       CALL MSETWIN(madis_profiler_minbck,madis_profiler_minfwd,madis_profiler_recwin,istatus)
1384       if(istatus /= success_p)stop
1385       CALL MNPNSTA(atime,nsta,stanam,wmoid,lat,lon,elev,timeob,istatus)
1386       if(istatus == success_p)then
1387          do i = 1, nsta
1388             if(abs(lat(i)).lt.90.and.abs(lon(i)).le.180)then
1389                platform % info % lat = lat(i)
1390                platform % info % lon = lon(i)
1391                CALL da_llxy(platform%info,platform%loc,outside,outside_all)
1392                if (outside_all)then
1393                   keep_n(i) = 0 ! 0 = Not in domain.
1394                else
1395                   ntotal(sound) = ntotal(sound) + 1
1396                   if(outside)then
1397                      keep(i) = 0
1398                   else
1399                      keep_n(i) = 1 ! 1 = Process.
1400                      ntmp_n = ntmp_n + 1
1401                      if (global .and. (platform%loc%i < ids .or. platform%loc%i >= ide)) then
1402                         nlocal(sound) = nlocal(sound) + 2
1403                      else
1404                         nlocal(sound) = nlocal(sound) + 1
1405                      endif
1406                   endif
1407                endif
1408             else
1409                keep_n(i) = 0
1410             endif
1411          enddo
1412       endif
1413
1414    endif
1415
1416 ! Keep track of nlocal and notal for use by PROCESS_DATASET when it needs to
1417 ! know these values over several calls.
1418
1419    nl_npn = nlocal(sound) - nl_raob
1420    nt_npn = ntotal(sound) - nt_raob
1421
1422    write(unit=message(1),fmt='(A,I10)')'ntotal_npn=',ntotal(sound)
1423    call da_message(message(1:1))
1424 ! MAP
1425
1426    ntmp_m = 0
1427
1428    if(map_init)then
1429
1430 ! Set the time window for this dataset with the user's choices.
1431
1432       CALL MSETWIN(madis_profiler_minbck,madis_profiler_minfwd,madis_profiler_recwin,istatus)
1433       if(istatus /= success_p)stop
1434       CALL MMAPSTA(atime,nsta,stanam,wmoid,lat,lon,elev,timeob,provdr, &
1435                    istatus)
1436       if(istatus == success_p)then
1437          do i = 1, nsta
1438             if(abs(lat(i)).lt.90.and.abs(lon(i)).le.180)then
1439                platform % info % lat = lat(i)
1440                platform % info % lon = lon(i)
1441                CALL da_llxy(platform%info,platform%loc,outside,outside_all)
1442                if (outside_all)then
1443                   keep_m(i) = 0 ! 0 = Not in domain.
1444                else
1445                   ntotal(sound) = ntotal(sound) + 1
1446                   if(outside)then
1447                      keep(i) = 0
1448                   else
1449                      keep_m(i) = 1 ! 1 = Process.
1450                      ntmp_m = ntmp_m + 1
1451                      if (global .and. (platform%loc%i < ids .or. platform%loc%i >= ide)) then
1452                         nlocal(sound) = nlocal(sound) + 2
1453                       else
1454                         nlocal(sound) = nlocal(sound) + 1
1455                       endif
1456                   endif
1457                endif
1458             else
1459                keep_m(i) = 0
1460             endif
1461          enddo
1462       endif
1463
1464    endif
1465
1466 ! Keep track of nlocal and notal for use by PROCESS_DATASET when it needs to
1467 ! know these values over several calls.
1468
1469    nl_map = nlocal(sound) - nl_raob - nl_npn
1470    nt_map = ntotal(sound) - nt_raob - nt_npn
1471
1472    write(unit=message(1),fmt='(A,I10)')'ntotal_map=',ntotal(sound)
1473    call da_message(message(1:1))
1474 ! If we didn't find anything, we're done with this WRF dataset.
1475
1476    if (nlocal(sound) > 0) then
1477       allocate (iv % sound(1:nlocal(sound)))
1478       allocate (iv % sonde_sfc(1:nlocal(sound)))
1479    endif
1480
1481    fm = 35
1482    iv%info(sound)%max_lev = max_raoblev
1483    iv%info(sonde_sfc)%max_lev=1
1484
1485 ! Since we have the MAP station info, start with that dataset.
1486
1487    if(ntmp_m /= 0)then
1488
1489 ! Allocate the local memory for the MADIS arrays that have a level dimension.
1490
1491       allocate (p(max_maplev*nsta))
1492       allocate (q(max_maplev*nsta))
1493       allocate (t(max_maplev*nsta))
1494       allocate (z(max_maplev*nsta))
1495       allocate (dd(max_maplev*nsta))
1496       allocate (ff(max_maplev*nsta))
1497       allocate (qca(max_maplev*nsta))
1498       allocate (qcr(max_maplev*nsta))
1499       allocate (qcd(max_maplev*nsta))
1500
1501       nlev = max_maplev
1502
1503 ! Read in the observations.
1504
1505       p(1:nsta*nlev) = miss_p
1506       CALL MGETMAP(atime,'P',ntmp,nobs,nlevels,p,qcd,qca,qcr,istatus)
1507       dd(1:nsta*nlev) = miss_p
1508       CALL MGETMAP(atime,'DD',ntmp,nobs,nlevels,dd,qcd,qca,qcr,istatus)
1509       ff(1:nsta*nlev) = miss_p
1510       CALL MGETMAP(atime,'FF',ntmp,nobs,nlevels,ff,qcd,qca,qcr,istatus)
1511       z(1:nsta*nlev) = miss_p
1512       CALL MGETMAP(atime,'HT',ntmp,nobs,nlevels,z,qcd,qca,qcr,istatus)
1513
1514 ! Fill in missing flags for the variables that aren't in this dataset,
1515 ! process it, deallocate the memory.
1516
1517       slp(1:nsta) = miss_p
1518       pw(1:nsta) = miss_p
1519       q(1:nsta*nlev) = miss_p
1520       t(1:nsta*nlev) = miss_p
1521       sfclev(1:nsta) = 0
1522
1523       CALL PROCESS_DATASET(nsta,nlev,fm,keep_m,lat,lon,elev,timeob, &
1524                            slp,pw,stanam,p,q,t,dd,ff,z,nlocal,ntotal, &
1525                            0,0,nlevels,iv,nlocal_out,ntotal_out,sfclev)
1526       deallocate (p)
1527       deallocate (q)
1528       deallocate (t)
1529       deallocate (z)
1530       deallocate (dd)
1531       deallocate (ff)
1532       deallocate (qca)
1533       deallocate (qcr)
1534       deallocate (qcd)
1535
1536    else
1537
1538       nlocal_out = 0
1539       ntotal_out = 0
1540
1541    endif
1542
1543 ! Now do the NPN dataset.
1544
1545    if(ntmp_n /= 0)then
1546
1547 ! Allocate the local memory for the MADIS arrays that have a level dimension.
1548
1549       if(madis_profiler_db(1:3) == 'FSL')then
1550          nlev = max_npnlev_madis
1551       else
1552          nlev = max_npnlev_awips
1553       endif
1554
1555 ! Set the time window and load in the station info.
1556
1557       CALL MSETWIN(madis_profiler_minbck,madis_profiler_minfwd,madis_profiler_recwin,istatus)
1558       if(istatus /= success_p)stop
1559       CALL MNPNSTA(atime,nsta,stanam,wmoid,lat,lon,elev,timeob,istatus)
1560       if(istatus /= success_p)stop   ! Can't possibly happen, so stop if it does.
1561
1562       allocate (p(nlev*nsta))
1563       allocate (q(nlev*nsta))
1564       allocate (t(nlev*nsta))
1565       allocate (z(nlev*nsta))
1566       allocate (dd(nlev*nsta))
1567       allocate (ff(nlev*nsta))
1568       allocate (qca(nlev*nsta))
1569       allocate (qcr(nlev*nsta))
1570       allocate (qcd(nlev*nsta))
1571
1572 ! Read in the observations.
1573
1574       p(1:nsta*nlev) = miss_p
1575       CALL MGETNPN(atime,'P',ntmp,nobs,nlevels,p,qcd,qca,qcr,istatus)
1576       dd(1:nsta*nlev) = miss_p
1577       CALL MGETNPN(atime,'DD',ntmp,nobs,nlevels,dd,qcd,qca,qcr,istatus)
1578       ff(1:nsta*nlev) = miss_p
1579       CALL MGETNPN(atime,'FF',ntmp,nobs,nlevels,ff,qcd,qca,qcr,istatus)
1580       z(1:nsta*nlev) = miss_p
1581       CALL MGETNPN(atime,'HT',ntmp,nobs,nlevels,z,qcd,qca,qcr,istatus)
1582
1583 ! Fill in missing flags for the variables that aren't in this dataset,
1584 ! process it, deallocate the memory.
1585
1586       slp(1:nsta) = miss_p
1587       pw(1:nsta) = miss_p
1588       q(1:nsta*nlev) = miss_p
1589       t(1:nsta*nlev) = miss_p
1590       sfclev(1:nsta) = 0
1591
1592       CALL PROCESS_DATASET(nsta,nlev,fm,keep_n,lat,lon,elev,timeob, &
1593                            slp,pw,stanam,p,q,t,dd,ff,z,nlocal,ntotal, &
1594                            nlocal_out,ntotal_out,nlevels,iv,nlocal_out, &
1595                            ntotal_out,sfclev)
1596
1597       deallocate (p)
1598       deallocate (q)
1599       deallocate (t)
1600       deallocate (z)
1601       deallocate (dd)
1602       deallocate (ff)
1603       deallocate (qca)
1604       deallocate (qcr)
1605       deallocate (qcd)
1606
1607    endif
1608
1609 ! Do the RAOB dataset.
1610
1611    if(ntmp_r /= 0)then
1612
1613       nlev = max_raoblev
1614
1615 ! Set the time window and load in the station info.
1616
1617       CALL MSETWIN(madis_sound_minbck,madis_sound_minfwd,madis_sound_recwin,istatus)
1618       if(istatus /= success_p)stop
1619       CALL MRAOBSTA(atime,nsta,stanam,wmoid,lat,lon,elev,timeob,istatus)
1620       if(istatus /= success_p)stop   ! Can't possibly happen, so stop if it does.
1621
1622 ! Allocate the local memory for the MADIS arrays that have a level dimension.
1623
1624       allocate (p(max_raoblev*nsta))
1625       allocate (q(max_raoblev*nsta))
1626       allocate (t(max_raoblev*nsta))
1627       allocate (z(max_raoblev*nsta))
1628       allocate (dd(max_raoblev*nsta))
1629       allocate (ff(max_raoblev*nsta))
1630       allocate (qca(max_raoblev*nsta))
1631       allocate (qcr(max_raoblev*nsta))
1632       allocate (qcd(max_raoblev*nsta))
1633       allocate (levtype(max_raoblev*nsta))
1634
1635 ! Read in the observations.
1636
1637       p(1:nsta*nlev) = miss_p
1638       CALL MGETRAOB(atime,'P',ntmp,nobs,nlevels,p,qcd,qca,qcr,istatus)
1639       q(1:nsta*nlev) = miss_p
1640       CALL MGETRAOB(atime,'Q',ntmp,nobs,nlevels,q,qcd,qca,qcr,istatus)
1641       t(1:nsta*nlev) = miss_p
1642       CALL MGETRAOB(atime,'T',ntmp,nobs,nlevels,t,qcd,qca,qcr,istatus)
1643       dd(1:nsta*nlev) = miss_p
1644       CALL MGETRAOB(atime,'DD',ntmp,nobs,nlevels,dd,qcd,qca,qcr,istatus)
1645       ff(1:nsta*nlev) = miss_p
1646       CALL MGETRAOB(atime,'FF',ntmp,nobs,nlevels,ff,qcd,qca,qcr,istatus)
1647       z(1:nsta*nlev) = miss_p
1648       CALL MGETRAOB(atime,'HT',ntmp,nobs,nlevels,z,qcd,qca,qcr,istatus)
1649       levtype(1:nsta*nlev) = miss_p
1650       CALL MGETRAOB(atime,'LEVTYPE',ntmp,nobs,nlevels,levtype,qcd,qca,qcr,istatus)
1651
1652 ! Knock out any stations with zero levels of actual data.  Also find the
1653 ! surface level for each station and fill in its P into SLP.
1654
1655       do i = 1, nsta
1656
1657          slp(i) = miss_p
1658
1659          if(nlevels(i) == 0)then
1660
1661             keep_r(i) = 0
1662
1663          else
1664
1665             CALL FIND_SFC_LEVEL(i,nlev,nsta,nlevels(i),levtype,p,slp(i),sfclev(i))
1666
1667          endif
1668
1669       enddo
1670
1671 ! Fill in missing flags for the variables that aren't in this dataset,
1672 ! process the data, deallocate the memory.
1673
1674  260  pw(1:nsta) = miss_p
1675
1676       CALL PROCESS_DATASET(nsta,nlev,fm,keep_r,lat,lon,elev,timeob, &
1677                            slp,pw,stanam,p,q,t,dd,ff,z,nlocal,ntotal, &
1678                            nlocal_out,ntotal_out,nlevels,iv,nlocal_out, &
1679                            ntotal_out,sfclev)
1680       iv%info(sound)%nlocal = nlocal_out
1681
1682       deallocate (p)
1683       deallocate (q)
1684       deallocate (t)
1685       deallocate (z)
1686       deallocate (dd)
1687       deallocate (ff)
1688       deallocate (qca)
1689       deallocate (qcr)
1690       deallocate (qcd)
1691       deallocate (levtype)
1692
1693    endif
1694
1695    iv%info(sound)%plocal(1) = nlocal_out
1696    iv%info(sound)%ptotal(1) = ntotal(sound)
1697    iv%info(sound)%ntotal = ntotal(sound)
1698    iv%info(sonde_sfc)%plocal(1) = nlocal_out
1699    iv%info(sonde_sfc)%ptotal(1) = ntotal(sound)
1700    iv%info(sonde_sfc)%ntotal = ntotal(sound)
1701    iv%info(sonde_sfc)%nlocal = nlocal_out
1702
1703  300  continue
1704
1705 !***** Put the profiler or raob dataset here (at statement 300) when these
1706 !      get split up.
1707
1708 !------------------------------------------------------------------------------
1709 ! MADIS SATWND dataset --> WRF GeoAMV dataset
1710 !------------------------------------------------------------------------------
1711 ! If we're not using this WRF dataset, skip to the next one.
1712
1713  450  if(.not. use_geoamvobs)go to 550
1714
1715 ! Allocate the local memory for the MADIS arrays that have a level dimension.
1716
1717    allocate (p(max_satw))
1718    allocate (q(max_satw))
1719    allocate (t(max_satw))
1720    allocate (z(max_satw))
1721    allocate (dd(max_satw))
1722    allocate (ff(max_satw))
1723    allocate (qca(max_satw))
1724    allocate (qcr(max_satw))
1725    allocate (qcd(max_satw))
1726
1727 ! Set the time window for this dataset with the user's choices.
1728
1729    CALL MSETWIN(madis_geoamv_minbck,madis_geoamv_minfwd,madis_geoamv_recwin,istatus)
1730    if(istatus /= success_p)stop
1731
1732 ! Load in the stations with data for this hour.  Note that the station call
1733 ! covers all variable types.  Therefore, when we read each variable type's
1734 ! data separately, we'll have to overlay that into the dd/ff arrays so
1735 ! we don't wipe out data read for an earlier variable type.  We'll use the
1736 ! not-applicable pw/slp arrays for temporary storage.
1737
1738    CALL MSATWNDSTA(atime,nsta,provdr,lat,lon,p,timeob,istatus)
1739    if(istatus /= success_p)go to 500
1740
1741 ! Read in the observations, appending each variable type as we go along.
1742
1743    dd(1:nsta) = miss_p
1744    ff(1:nsta) = miss_p
1745
1746    do i = 1, nsatw_vartype
1747
1748       CALL MGETSATWND(atime,satw_var(1,i)(1:satw_var_len(i)),ntmp,nobs, &
1749                       pw,qcd,qca,qcr,istatus)
1750       CALL MGETSATWND(atime,satw_var(2,i)(1:satw_var_len(i)),ntmp,nobs, &
1751                       slp,qcd,qca,qcr,istatus)
1752
1753       if(istatus == success_p)then
1754
1755          do j = 1, ntmp
1756
1757             if(pw(j) /= miss_p)then
1758
1759                dd(j) = pw(j)
1760                ff(j) = slp(j)
1761
1762             endif
1763
1764          enddo
1765
1766       endif
1767
1768    enddo
1769
1770 ! Fill in missing flags for the variables that aren't in this dataset,
1771 ! and fill in 1 level for each station for use by the PROCESS_DATASET routine.
1772
1773    pw(1:nsta) = miss_p
1774    slp(1:nsta) = miss_p
1775    elev(1:nsta) = miss_p
1776    z(1:nsta) = miss_p
1777    nlevels(1:nsta) = 1
1778
1779 ! Determine which ones are inside the user's domain, and then allocate
1780 ! memory for the IV data for this dataset.  Also toss those that don't
1781 ! have any data -- we only have to check one of the wind variables,
1782 ! if it exists, then all variables will exist.
1783
1784    ntmp = 0                      ! Count of records inside domain.
1785
1786    do i = 1, nsta
1787
1788       if(abs(lat(i)).lt.90.and.abs(lon(i)).le.180)then
1789          platform % info % lat = lat(i)
1790          platform % info % lon = lon(i)
1791          CALL da_llxy(platform%info,platform%loc,outside,outside_all)
1792          if (outside_all .or. dd(i) == miss_p)then
1793             keep(i) = 0         ! 0 = Not in domain or no data.
1794          else
1795             ntotal(geoamv) = ntotal(geoamv) + 1
1796             if(outside)then
1797                keep(i) = 0
1798             else
1799                keep(i) = 1         ! 1 = Process.
1800                ntmp = ntmp + 1
1801                if (global .and. (platform%loc%i < ids .or. platform%loc%i >= ide)) then
1802                   nlocal(geoamv) = nlocal(geoamv) + 2
1803                else
1804                   nlocal(geoamv) = nlocal(geoamv) + 1
1805                endif
1806             endif
1807          endif
1808       else
1809          keep(i) = 0
1810       endif
1811
1812    enddo
1813
1814    if(nlocal(geoamv) > 0)then
1815
1816       allocate (iv % geoamv(1:nlocal(geoamv)))
1817
1818       fm = 88
1819       nlev = 1
1820       iv%info(geoamv)%max_lev = 1
1821
1822 ! Move the data out of the ingest arrays into the IV structure.
1823
1824       CALL PROCESS_DATASET(nsta,nlev,fm,keep,lat,lon,elev,timeob, &
1825                            slp,pw,provdr,p,q,t,dd,ff,z,nlocal,ntotal, &
1826                            0,0,nlevels,iv,nlocal_out,ntotal_out,sfclev)
1827       iv%info(geoamv)%nlocal = nlocal_out
1828
1829    endif
1830
1831 ! Finish up with the MADIS SATWND datasets by deallocating the
1832 ! levels-dimensioned arrays.
1833
1834  500  deallocate (p)
1835    deallocate (q)
1836    deallocate (t)
1837    deallocate (z)
1838    deallocate (dd)
1839    deallocate (ff)
1840    deallocate (qca)
1841    deallocate (qcr)
1842    deallocate (qcd)
1843
1844  550  continue              ! Placeholder for the next dataset that's added.
1845
1846    go to 8000
1847
1848 !----------------------------------------------------------------------------
1849 !  Report error reading namelist.
1850 !----------------------------------------------------------------------------
1851 8500  write(0,*)' ERROR:  ', errmsg(i), ' IOSTAT = ', IOST, &
1852                 ' while reading MADIS namelist'
1853    stop
1854
1855 !------------------------------------------------------------------------------
1856 ! Endgame processing
1857 !------------------------------------------------------------------------------
1858
1859  8000 if(pvdriostat /= 0)then
1860
1861 ! I/O error on opening or reading a provider-selection file.
1862
1863       write(0,*)'ERROR:  ',pvdriotype,' failure on ',trim(pvdrfile), &
1864                 ' IOSTAT = ',pvdriostat
1865       stop
1866
1867 ! If we didn't get *anything*, output a suitable message and stop.
1868
1869    else if(iv%info(sound)%nlocal+iv%info(geoamv)%nlocal+iv%info(airep)%nlocal+ &
1870            iv%info(gpspw)%nlocal+iv%info(metar)%nlocal+iv%info(ships)%nlocal &
1871            == 0)then
1872       write(0,*)'ERROR:  No MADIS observations were found.'
1873 !     stop
1874    endif
1875
1876 ! Close the MADIS files and deallocate the memory used for the
1877 ! MADIS arrays that are dimensioned only by station.
1878
1879    CALL MMADISCLOSE
1880    deallocate (stanam)
1881    deallocate (provdr)
1882    deallocate (timeob)
1883    deallocate (wmoid)
1884    deallocate (lat)
1885    deallocate (lon)
1886    deallocate (elev)
1887    deallocate (pw)
1888    deallocate (slp)
1889    deallocate (keep)
1890    deallocate (nlevels)
1891
1892 ! Output the same summary info as the MM5 ingest routine, but just
1893 ! for those datasets we handle.
1894
1895    if ( print_detail_obs ) then
1896
1897       write(6, fmt='(5x,a,i6,a)') &
1898            'Read:  ', iv%info(sound)%nlocal, ' SOUND reports,', &
1899            'Read:  ', iv%info(geoamv)%nlocal, ' Geo AMV reports,', &
1900            'Read:  ', iv%info(airep)%nlocal, ' AIREP reports,', &
1901            'Read:  ', iv%info(gpspw)%nlocal, ' GPSPW/GPSZD reports,', &
1902            'Read:  ', iv%info(metar)%nlocal, ' METAR reports,', &
1903            'Read:  ', iv%info(ships)%nlocal, ' SHIP  reports,'
1904
1905    endif
1906
1907 ! Output a dump of all the obs, when desired, for testing purposes.
1908
1909    if(madis_debug > 0)then
1910
1911       write(6,*)
1912       write(6,*)'METAR:  ',iv%info(metar)%nlocal
1913       write(6,*)
1914
1915       do i = 1, iv%info(metar)%nlocal
1916          write(6,1000)iv%info(metar)%name(i)(1:8), &
1917                       iv%info(metar)%date_char(i)(1:16), &
1918                       iv%info(metar)%platform(i), &
1919                       iv%info(metar)%lat(1,i),iv%info(metar)%lon(1,i), &
1920                       iv%info(metar)%elv(i),iv%info(metar)%pstar(i), &
1921                       iv%info(metar)%levels(i), &
1922                       iv%info(metar)%slp(i)%inv,iv%info(metar)%slp(i)%qc, &
1923                       iv%info(metar)%slp(i)%error,iv%info(metar)%pw(i)%inv, &
1924                       iv%info(metar)%pw(i)%qc,iv%info(metar)%pw(i)%error
1925  1000    format(a,1x,a,1x,a,f6.2,1x,f7.2,1x,f10.2,1x,f8.2,1x,i3,1x, &
1926                 f10.2,1x,i3,1x,f10.2,1x,f12.4,1x,i3,1x,f10.2)
1927          write(6,1001)1,' H/Z/U/V',iv%metar(i)%h,missing_data, &
1928                       iv%info(metar)%zk(1,i), &
1929                       iv%metar(i)%u%inv,iv%metar(i)%u%qc,iv%metar(i)%u%error, &
1930                       iv%metar(i)%v%inv,iv%metar(i)%v%qc,iv%metar(i)%v%error
1931  1001    format(i3,a,3(1x,f11.3,1x,i3,1x,f11.3))
1932          write(6,1001)1,' P/T/Q  ',iv%metar(i)%p%inv,iv%metar(i)%p%qc, &
1933                       iv%metar(i)%p%error,iv%metar(i)%t%inv,iv%metar(i)%t%qc, &
1934                       iv%metar(i)%t%error,iv%metar(i)%q%inv,iv%metar(i)%q%qc, &
1935                       iv%metar(i)%q%error
1936       enddo
1937
1938       write(6,*)
1939       write(6,*)'SHIPS:  ',iv%info(ships)%nlocal
1940       write(6,*)
1941
1942       do i = 1, iv%info(ships)%nlocal
1943          write(6,1000)iv%info(ships)%name(i)(1:8), &
1944                       iv%info(ships)%date_char(i)(1:16), &
1945                       iv%info(ships)%platform(i), &
1946                       iv%info(ships)%lat(1,i),iv%info(ships)%lon(1,i), &
1947                       iv%info(ships)%elv(i),iv%info(ships)%pstar(i), &
1948                       iv%info(ships)%levels(i), &
1949                       iv%info(ships)%slp(i)%inv,iv%info(ships)%slp(i)%qc, &
1950                       iv%info(ships)%slp(i)%error,iv%info(ships)%pw(i)%inv, &
1951                       iv%info(ships)%pw(i)%qc,iv%info(ships)%pw(i)%error
1952          write(6,1001)1,' H/Z/U/V',iv%ships(i)%h,missing_data, &
1953                       iv%info(ships)%zk(1,i), &
1954                       iv%ships(i)%u%inv,iv%ships(i)%u%qc,iv%ships(i)%u%error, &
1955                       iv%ships(i)%v%inv,iv%ships(i)%v%qc,iv%ships(i)%v%error
1956          write(6,1001)1,' P/T/Q  ',iv%ships(i)%p%inv,iv%ships(i)%p%qc, &
1957                       iv%ships(i)%p%error,iv%ships(i)%t%inv,iv%ships(i)%t%qc, &
1958                       iv%ships(i)%t%error,iv%ships(i)%q%inv,iv%ships(i)%q%qc, &
1959                       iv%ships(i)%q%error
1960       enddo
1961
1962       write(6,*)
1963       write(6,*)'GPSPW:  ',iv%info(gpspw)%nlocal
1964       write(6,*)
1965
1966       do i = 1, iv%info(gpspw)%nlocal
1967          write(6,1002)iv%info(gpspw)%name(i)(1:8), &
1968                       iv%info(gpspw)%date_char(i)(1:16), &
1969                       iv%info(gpspw)%platform(i), &
1970                       iv%info(gpspw)%lat(1,i),iv%info(gpspw)%lon(1,i), &
1971                       iv%info(gpspw)%elv(i),iv%info(gpspw)%pstar(i), &
1972                       iv%info(gpspw)%levels(i), &
1973                       iv%info(gpspw)%slp(i)%inv,iv%info(gpspw)%slp(i)%qc, &
1974                       iv%info(gpspw)%slp(i)%error,iv%info(gpspw)%pw(i)%inv, &
1975                       iv%info(gpspw)%pw(i)%qc,iv%info(gpspw)%pw(i)%error, &
1976                       iv%gpspw(i)%tpw%inv,iv%gpspw(i)%tpw%qc, &
1977                       iv%gpspw(i)%tpw%error
1978  1002    format(a,1x,a,1x,a,f6.2,1x,f7.2,1x,f7.2,1x,f8.2,1x,i3,1x, &
1979                 f10.2,1x,i3,1x,f10.2,1x,f12.4,1x,i3,1x,f10.2,1x, &
1980                 f12.4,1x,i3,1x,f10.2)
1981       enddo
1982
1983       write(6,*)
1984       write(6,*)'SOUND:  ',iv%info(sound)%nlocal
1985       write(6,*)
1986
1987       do i = 1, iv%info(sound)%nlocal
1988          write(6,1000)iv%info(sound)%name(i)(1:8), &
1989                       iv%info(sound)%date_char(i)(1:16), &
1990                       iv%info(sound)%platform(i), &
1991                       iv%info(sound)%lat(1,i),iv%info(sound)%lon(1,i), &
1992                       iv%info(sound)%elv(i),iv%info(sound)%pstar(i), &
1993                       iv%info(sound)%levels(i), &
1994                       iv%info(sound)%slp(i)%inv,iv%info(sound)%slp(i)%qc, &
1995                       iv%info(sound)%slp(i)%error,iv%info(sound)%pw(i)%inv, &
1996                       iv%info(sound)%pw(i)%qc,iv%info(sound)%pw(i)%error
1997          do j = 1, iv%info(sound)%levels(i)
1998             write(6,1001)j,' H/Z/U/V',iv%sound(i)%h(j),missing_data, &
1999                          iv%info(sound)%zk(j,i), &
2000                          iv%sound(i)%u(j)%inv,iv%sound(i)%u(j)%qc, &
2001                          iv%sound(i)%u(j)%error,iv%sound(i)%v(j)%inv, &
2002                          iv%sound(i)%v(j)%qc,iv%sound(i)%v(j)%error
2003             write(6,1001)j,' P/T/Q  ',iv%sound(i)%p(j),missing_data,missing_r, &
2004                          iv%sound(i)%t(j)%inv, &
2005                          iv%sound(i)%t(j)%qc,iv%sound(i)%t(j)%error, &
2006                          iv%sound(i)%q(j)%inv,iv%sound(i)%q(j)%qc, &
2007                          iv%sound(i)%q(j)%error
2008          enddo
2009       enddo
2010
2011       write(6,*)
2012       write(6,*)'SONDE_SFC:  ',iv%info(sonde_sfc)%nlocal
2013       write(6,*)
2014
2015       do i = 1, iv%info(sonde_sfc)%nlocal
2016          write(6,1000)iv%info(sonde_sfc)%name(i)(1:8), &
2017                       iv%info(sonde_sfc)%date_char(i)(1:16), &
2018                       iv%info(sonde_sfc)%platform(i), &
2019                       iv%info(sonde_sfc)%lat(1,i),iv%info(sonde_sfc)%lon(1,i), &
2020                       iv%info(sonde_sfc)%elv(i),iv%info(sonde_sfc)%pstar(i), &
2021                       iv%info(sonde_sfc)%levels(i), &
2022                       iv%info(sonde_sfc)%slp(i)%inv,iv%info(sonde_sfc)%slp(i)%qc, &
2023                       iv%info(sonde_sfc)%slp(i)%error,iv%info(sonde_sfc)%pw(i)%inv, &
2024                       iv%info(sonde_sfc)%pw(i)%qc,iv%info(sonde_sfc)%pw(i)%error
2025          write(6,1001)1,' H/Z/U/V',iv%sonde_sfc(i)%h,missing_data, &
2026                       iv%info(sonde_sfc)%zk(1,i), &
2027                       iv%sonde_sfc(i)%u%inv,iv%sonde_sfc(i)%u%qc,iv%sonde_sfc(i)%u%error, &
2028                       iv%sonde_sfc(i)%v%inv,iv%sonde_sfc(i)%v%qc,iv%sonde_sfc(i)%v%error
2029          write(6,1001)1,' P/T/Q  ',iv%sonde_sfc(i)%p%inv,iv%sonde_sfc(i)%p%qc, &
2030                       iv%sonde_sfc(i)%p%error,iv%sonde_sfc(i)%t%inv,iv%sonde_sfc(i)%t%qc, &
2031                       iv%sonde_sfc(i)%t%error,iv%sonde_sfc(i)%q%inv,iv%sonde_sfc(i)%q%qc, &
2032                       iv%sonde_sfc(i)%q%error
2033       enddo
2034
2035       write(6,*)
2036       write(6,*)'AIREP:  ',iv%info(airep)%nlocal
2037       write(6,*)
2038
2039       do i = 1, iv%info(airep)%nlocal
2040          write(6,1000)iv%info(airep)%name(i)(1:8), &
2041                       iv%info(airep)%date_char(i)(1:16), &
2042                       iv%info(airep)%platform(i), &
2043                       iv%info(airep)%lat(1,i),iv%info(airep)%lon(1,i), &
2044                       iv%info(airep)%elv(i),iv%info(airep)%pstar(i), &
2045                       iv%info(airep)%levels(i), &
2046                       iv%info(airep)%slp(i)%inv,iv%info(airep)%slp(i)%qc, &
2047                       iv%info(airep)%slp(i)%error,iv%info(airep)%pw(i)%inv, &
2048                       iv%info(airep)%pw(i)%qc,iv%info(airep)%pw(i)%error
2049          write(6,1001)1,' H/Z/U/V',iv%airep(i)%h(1),missing_data, &
2050                       iv%info(airep)%zk(1,i), &
2051                       iv%airep(i)%u(1)%inv,iv%airep(i)%u(1)%qc, &
2052                       iv%airep(i)%u(1)%error,iv%airep(i)%v(1)%inv, &
2053                       iv%airep(i)%v(1)%qc,iv%airep(i)%v(1)%error
2054          write(6,1001)1,' P/T    ',iv%airep(i)%p(1),missing_data, &
2055                       missing_r,iv%airep(i)%t(1)%inv, &
2056                       iv%airep(i)%t(1)%qc,iv%airep(i)%t(1)%error
2057       enddo
2058
2059       write(6,*)
2060       write(6,*)'GEOAMV:  ',iv%info(geoamv)%nlocal
2061       write(6,*)
2062
2063       do i = 1, iv%info(geoamv)%nlocal
2064          write(6,1000)iv%info(geoamv)%name(i)(1:8), &
2065                       iv%info(geoamv)%date_char(i)(1:16), &
2066                       iv%info(geoamv)%platform(i), &
2067                       iv%info(geoamv)%lat(1,i),iv%info(geoamv)%lon(1,i), &
2068                       iv%info(geoamv)%elv(i),iv%info(geoamv)%pstar(i), &
2069                       iv%info(geoamv)%levels(i), &
2070                       iv%info(geoamv)%slp(i)%inv,iv%info(geoamv)%slp(i)%qc, &
2071                       iv%info(geoamv)%slp(i)%error,iv%info(geoamv)%pw(i)%inv, &
2072                       iv%info(geoamv)%pw(i)%qc,iv%info(geoamv)%pw(i)%error
2073          write(6,1001)1,' P/Z/U/V',iv%geoamv(i)%p(1),missing_data, &
2074                       iv%info(geoamv)%zk(1,i), &
2075                       iv%geoamv(i)%u(1)%inv,iv%geoamv(i)%u(1)%qc, &
2076                       iv%geoamv(i)%u(1)%error,iv%geoamv(i)%v(1)%inv, &
2077                       iv%geoamv(i)%v(1)%qc,iv%geoamv(i)%v(1)%error
2078       enddo
2079
2080    endif         ! if(madis_debug > 0)
2081
2082 ! Fill in the IV time.
2083
2084    iv%time = 1
2085
2086
2087 ! Copy data out of the IV structure into the OB structure.
2088
2089    CALL da_fill_obs_structures(iv,ob)
2090
2091    if(madis_debug > 0)then
2092
2093 ! If we're in debugging mode, indicate that we made it to the bottom of
2094 ! the routine OK.  Also, if requested, just stop after the obs ingest.
2095
2096       write(6,*)
2097       write(6,*)'Successful MADIS ingest.'
2098
2099       if(madis_debug == 2)stop
2100
2101    endif
2102
2103    return
2104
2105 #else
2106    call da_error(__FILE__,__LINE__,(/"Needs to be compiled with a MADIS library"/))
2107 #endif
2108
2109    if (trace_use) CALL da_trace_exit("da_setup_obs_structures_madis")
2110
2111 END SUBROUTINE da_setup_obs_structures_madis
2112
2113 #ifdef MADIS
2114
2115 SUBROUTINE PROCESS_DATASET(nsta,nlev,fm,keep,lat,lon,elev,timeob, &
2116                            slp,pw,stanam,p,q,t,dd,ff,z, &
2117                            nlocal,ntotal,nloc_init,ntot_init,nlevels,iv, &
2118                            nlocal_out,ntotal_out,sfclev)
2119
2120 !------------------------------------------------------------------------------
2121 ! PURPOSE: Move obs out of MADIS ingest arrays into the WRF innovation
2122 !          vector structure, including attendant processing that fills in the
2123 !          location within the user's domain, etc., for one MADIS dataset.
2124 !
2125 ! METHOD:  Move obs & metadata from MADIS ingest arrays to MM5-like
2126 !          PLATFORM structure, manipulate PLATFORM with DA_Obs routines,
2127 !          then move the resulting stuff into IV.
2128 !
2129 ! HISTORY: 05-21-03  M. Barth         Original version.
2130 !          08-27-04  M. Barth         Removed psfc and added sonde_sfc.
2131 !          02-06-06  M. Barth         Removed arguments no longer in use,
2132 !                                     changed xb to xp, removed DA_Obs_Count.
2133 !          05-30-08  M. Barth         Removed xp, changed DA_ll_to_xy to da_llxy,
2134 !                                     used the new iv%info structures,
2135 !                                     changed ob_type structure to iv_type,
2136 !                                     added nlocal and ntotal arrays and
2137 !                                     nloc_init, ntot_init.  Added
2138 !                                     use of FILL_IV_INFO subroutine.
2139 !          11-28-08 M. Barth          Since I finally got the darn thing to
2140 !                                     compile, I debugged the massive changes
2141 !                                     in version 3.
2142 !
2143 !
2144 ! PARENT_MODULE: da_setup_structures
2145 !------------------------------------------------------------------------------
2146
2147    use da_grid_definitions, only : da_ffdduv
2148    use da_tools, only : da_llxy
2149    use da_obs, only : da_obs_proc_station
2150    IMPLICIT NONE
2151
2152 ! ARGUMENTS:
2153
2154    INTEGER*4, INTENT(IN)        :: nsta               ! Number of stations
2155    INTEGER*4, INTENT(IN)        :: nlev               ! Number of levels
2156    INTEGER, INTENT(IN)          :: fm                 ! GTS FM number for this dataset
2157    INTEGER*4, INTENT(IN)        :: keep(nsta)         ! 1 = process this station
2158    REAL*4,    INTENT(IN)        :: lat(nsta)          ! Latitude (deg N)
2159    REAL*4,    INTENT(IN)        :: lon(nsta)          ! Longitude (deg E)
2160    REAL*4,    INTENT(IN)        :: elev(nsta)         ! Elevation (m)
2161    CHARACTER*9, INTENT(IN)      :: timeob(nsta)       ! Time (YYJJJHHMM)
2162    REAL*4,    INTENT(IN)        :: slp(nsta)          ! Sea level pressure (Pa)
2163    REAL*4,    INTENT(IN)        :: pw(nsta)           ! Precipitable water vapor (cm)
2164    CHARACTER*10, INTENT(IN)     :: stanam(nsta)       ! Station ID
2165    REAL*4,    INTENT(IN)        :: p(nlev,nsta)       ! Pressure (Pa)
2166    REAL*4,    INTENT(IN)        :: q(nlev,nsta)       ! Mixing ratio (kg/kg)
2167    REAL*4,    INTENT(IN)        :: t(nlev,nsta)       ! Temperature (K)
2168    REAL*4,    INTENT(INOUT)     :: dd(nlev,nsta)      ! Wind speed (m/s)
2169    REAL*4,    INTENT(INOUT)     :: ff(nlev,nsta)      ! Wind direction (0-360 N)
2170    REAL*4,    INTENT(IN)        :: z(nlev,nsta)       ! Height (m)
2171    INTEGER*4, INTENT(IN)        :: nlocal(*)          ! Obs counts per dataset, in domain
2172    INTEGER*4, INTENT(IN)        :: ntotal(*)          ! Obs counts per dataset, in or out
2173    INTEGER*4, INTENT(IN)        :: nloc_init          ! Initial value for nlocal index
2174    INTEGER*4, INTENT(IN)        :: ntot_init          ! Initial value for ntotal index
2175    INTEGER*4, INTENT(IN)        :: nlevels(nsta)      ! Num levels per station
2176    TYPE (iv_type), INTENT(INOUT) :: iv                ! Innovation vector
2177    INTEGER*4, INTENT(OUT)       :: nlocal_out         ! Obs count in domain, and with useful data
2178    INTEGER*4, INTENT(OUT)       :: ntotal_out         ! Obs count, and with useful data
2179    INTEGER*4, INTENT(IN)        :: sfclev(nsta)       ! Index to raob sfc level per station
2180
2181 ! LOCAL:
2182
2183    TYPE (multi_level_type)      :: platform
2184    INTEGER                      :: i,j,k,n,istatus,obs_index,ntot,ndup,l
2185    REAL*4, PARAMETER            :: miss_p    = 999999.0  ! Missing data flag
2186    CHARACTER*13                 :: tstr
2187    LOGICAL                      :: outside,cycle_report,anyvar,outside_all
2188    REAL                         :: xff,xdd,xlon
2189    INTEGER*4                    :: sfclev_j,m
2190
2191 !------------------------------------------------------------------------------
2192 ! The MM5-like PLATFORM structure is used to move data out of the MADIS
2193 ! ingest arrays into a structure that can then be manipulated by the DA_Obs
2194 ! routines.  Then, the single-station data in PLATFORM is moved to where
2195 ! it belongs in the output innovation vector, IV.
2196 !------------------------------------------------------------------------------
2197 ! Fill in dataset-specific info.
2198
2199    SELECT CASE(fm)
2200       CASE(15)
2201          platform % info % platform = 'FM-015 METAR'
2202          obs_index = metar
2203       CASE(13)
2204          platform % info % platform = 'FM-013 SHIP '
2205          obs_index = ships
2206       CASE(96)
2207          platform % info % platform = 'FM-096 AIREP'
2208          obs_index = airep
2209       CASE(35)
2210          platform % info % platform = 'FM-035 SOUND'
2211          obs_index = sound
2212       CASE(111)
2213          platform % info % platform = 'FM-111 GPSPW'
2214          obs_index = gpspw
2215       CASE(88)
2216          platform % info % platform = 'FM-088 SATOB'
2217          obs_index = geoamv
2218    END SELECT
2219
2220 ! Only allocate stuff the first time we're here for a WRF dataset, since multiple
2221 ! MADIS datasets may comprise a single WRF dataset (e.g., MAP, NPN, raob).
2222
2223    if(nloc_init.eq.0)then
2224
2225      iv%info(obs_index)%nlocal = nlocal(obs_index)
2226      iv%info(obs_index)%ntotal = ntotal(obs_index)
2227      if (iv%info(obs_index)%nlocal > 0) then
2228         allocate (iv%info(obs_index)%name(iv%info(obs_index)%nlocal))
2229         allocate (iv%info(obs_index)%platform(iv%info(obs_index)%nlocal))
2230         allocate (iv%info(obs_index)%id(iv%info(obs_index)%nlocal))
2231         allocate (iv%info(obs_index)%date_char(iv%info(obs_index)%nlocal))
2232         allocate (iv%info(obs_index)%levels(iv%info(obs_index)%nlocal))
2233         allocate (iv%info(obs_index)%lat(iv%info(obs_index)%max_lev,iv%info(obs_index)%nlocal))
2234         allocate (iv%info(obs_index)%lon(iv%info(obs_index)%max_lev,iv%info(obs_index)%nlocal))
2235         allocate (iv%info(obs_index)%elv(iv%info(obs_index)%nlocal))
2236         allocate (iv%info(obs_index)%pstar(iv%info(obs_index)%nlocal))
2237         allocate (iv%info(obs_index)%slp(iv%info(obs_index)%nlocal))
2238         allocate (iv%info(obs_index)%pw(iv%info(obs_index)%nlocal))
2239         allocate (iv%info(obs_index)%x  (kms:kme,iv%info(obs_index)%nlocal))
2240         allocate (iv%info(obs_index)%y  (kms:kme,iv%info(obs_index)%nlocal))
2241         allocate (iv%info(obs_index)%i  (kms:kme,iv%info(obs_index)%nlocal))
2242         allocate (iv%info(obs_index)%j  (kms:kme,iv%info(obs_index)%nlocal))
2243         allocate (iv%info(obs_index)%dx (kms:kme,iv%info(obs_index)%nlocal))
2244         allocate (iv%info(obs_index)%dxm(kms:kme,iv%info(obs_index)%nlocal))
2245         allocate (iv%info(obs_index)%dy (kms:kme,iv%info(obs_index)%nlocal))
2246         allocate (iv%info(obs_index)%dym(kms:kme,iv%info(obs_index)%nlocal))
2247         allocate (iv%info(obs_index)%k  (iv%info(obs_index)%max_lev,iv%info(obs_index)%nlocal))
2248         allocate (iv%info(obs_index)%dz (iv%info(obs_index)%max_lev,iv%info(obs_index)%nlocal))
2249         allocate (iv%info(obs_index)%dzm(iv%info(obs_index)%max_lev,iv%info(obs_index)%nlocal))
2250         allocate (iv%info(obs_index)%zk (iv%info(obs_index)%max_lev,iv%info(obs_index)%nlocal))
2251         allocate (iv%info(obs_index)%proc_domain(iv%info(obs_index)%max_lev,iv%info(obs_index)%nlocal))
2252         allocate (iv%info(obs_index)%obs_global_index(iv%info(obs_index)%nlocal))
2253         iv%info(obs_index)%proc_domain(:,:)  = .false.
2254         iv%info(obs_index)%zk(:,:)           = missing_r
2255
2256         if(fm == 35)then
2257
2258 ! Now allocate the sonde_sfc iv%info arrays.
2259
2260            allocate (iv%info(sonde_sfc)%name(iv%info(obs_index)%nlocal))
2261            allocate (iv%info(sonde_sfc)%platform(iv%info(obs_index)%nlocal))
2262            allocate (iv%info(sonde_sfc)%id(iv%info(obs_index)%nlocal))
2263            allocate (iv%info(sonde_sfc)%date_char(iv%info(obs_index)%nlocal))
2264            allocate (iv%info(sonde_sfc)%levels(iv%info(obs_index)%nlocal))
2265            allocate (iv%info(sonde_sfc)%lat(iv%info(sonde_sfc)%max_lev,iv%info(obs_index)%nlocal))
2266            allocate (iv%info(sonde_sfc)%lon(iv%info(sonde_sfc)%max_lev,iv%info(obs_index)%nlocal))
2267            allocate (iv%info(sonde_sfc)%elv(iv%info(obs_index)%nlocal))
2268            allocate (iv%info(sonde_sfc)%pstar(iv%info(obs_index)%nlocal))
2269            allocate (iv%info(sonde_sfc)%slp(iv%info(obs_index)%nlocal))
2270            allocate (iv%info(sonde_sfc)%pw(iv%info(obs_index)%nlocal))
2271            allocate (iv%info(sonde_sfc)%x  (kms:kme,iv%info(obs_index)%nlocal))
2272            allocate (iv%info(sonde_sfc)%y  (kms:kme,iv%info(obs_index)%nlocal))
2273            allocate (iv%info(sonde_sfc)%i  (kms:kme,iv%info(obs_index)%nlocal))
2274            allocate (iv%info(sonde_sfc)%j  (kms:kme,iv%info(obs_index)%nlocal))
2275            allocate (iv%info(sonde_sfc)%dx (kms:kme,iv%info(obs_index)%nlocal))
2276            allocate (iv%info(sonde_sfc)%dxm(kms:kme,iv%info(obs_index)%nlocal))
2277            allocate (iv%info(sonde_sfc)%dy (kms:kme,iv%info(obs_index)%nlocal))
2278            allocate (iv%info(sonde_sfc)%dym(kms:kme,iv%info(obs_index)%nlocal))
2279            allocate (iv%info(sonde_sfc)%k  (iv%info(sonde_sfc)%max_lev,iv%info(obs_index)%nlocal))
2280            allocate (iv%info(sonde_sfc)%dz (iv%info(sonde_sfc)%max_lev,iv%info(obs_index)%nlocal))
2281            allocate (iv%info(sonde_sfc)%dzm(iv%info(sonde_sfc)%max_lev,iv%info(obs_index)%nlocal))
2282            allocate (iv%info(sonde_sfc)%zk (iv%info(sonde_sfc)%max_lev,iv%info(obs_index)%nlocal))
2283            allocate (iv%info(sonde_sfc)%proc_domain(iv%info(sonde_sfc)%max_lev,iv%info(obs_index)%nlocal))
2284            allocate (iv%info(sonde_sfc)%obs_global_index(iv%info(obs_index)%nlocal))
2285            iv%info(sonde_sfc)%proc_domain(:,:)  = .false.
2286            iv%info(sonde_sfc)%zk(:,:)           = missing_r
2287
2288        endif
2289
2290      endif
2291
2292    endif
2293
2294 ! Process each record that's in the user's domain.
2295 !
2296 ! Note that the PLATFORM structure includes some variables that
2297 ! aren't needed, so we don't deal with them:  RH, TD, SPEED.
2298
2299    n = nloc_init                 ! Index into iv.
2300    ntot = ntot_init              ! Value to use for current ntotal within iv-filling
2301                                  ! loop.  See explanation of duplication logic.
2302
2303    do i = 1, nsta
2304
2305       if(keep(i) == 1)then
2306
2307 ! Fill in the header-type info for this record.
2308
2309          platform % info % lat = lat(i)
2310          platform % info % lon = lon(i)
2311          platform % info % elv = elev(i)
2312          platform % info % name = stanam(i)
2313          platform % info % id = stanam(i)(1:5)
2314          platform % info % levels = nlevels(i)
2315
2316 ! Convert time from YYJJJHHMM to YYYYMMDD_HHMM, then put it into
2317 ! the info % date_char format of YYYY-MM-DD_HH:MM:SS.
2318
2319          CALL MTRNTIM(timeob(i),2,tstr,istatus)
2320          platform % info % date_char = tstr(1:4)//'-'//tstr(5:6)//'-'// &
2321                                        tstr(7:11)//':'//tstr(12:13)//':00'
2322
2323 ! Fill in PW and SLP.
2324
2325 ! ***** Note that I'm not sure if good QC is indicated by qc_good (1) or
2326 !       good_quality (0), or both... I'm guessing qc_good should work
2327 !       since da_read_obs_bufr accepts any value < 4.
2328
2329          if(pw(i) /= miss_p)then
2330             platform % loc % pw % inv = pw(i)
2331             platform % loc % pw % qc = qc_good
2332          else
2333             platform % loc % pw % inv = missing_r
2334             platform % loc % pw % qc = missing_data
2335          endif
2336          platform % loc % pw % error = .2    ! Same reference as in GET_OB_ERRORS.
2337
2338 ! Fill in SLP.
2339
2340          if(slp(i) /= miss_p)then
2341             platform % loc % slp % inv = slp(i)
2342             platform % loc % slp % qc = qc_good
2343          else
2344             platform % loc % slp % inv = missing_r
2345             platform % loc % slp % qc = missing_data
2346          endif
2347          platform % loc % slp % error = 200. ! Same reference as in GET_OB_ERRORS.
2348
2349 ! Fill in the data variables that are a function of level.
2350
2351          j = 1       ! J = level index for what's gonna get used
2352          sfclev_j = 0
2353
2354          do k = 1, nlevels(i)
2355
2356             if(z(k,i) /= miss_p)then
2357                platform % each(j) % height = z(k,i)
2358             else
2359                platform % each(j) % height = missing_r
2360             endif
2361
2362             if(p(k,i) /= miss_p)then
2363                platform % each(j) % p % inv = p(k,i)
2364                platform % each(j) % p % qc = qc_good
2365             else
2366                platform % each(j) % p % inv = missing_r
2367                platform % each(j) % p % qc = missing_data
2368             endif
2369
2370 ! If there isn't any P or Z, this level isn't going to make the cut, so
2371 ! skip the rest of its processing.
2372
2373             if(.not.(p(k,i) == miss_p .and. z(k,i) == miss_p))then
2374
2375                anyvar = .false.
2376
2377 ! If DD and FF exist, convert to U,V with rotation to the user's grid.
2378
2379                if(dd(k,i) /= miss_p .and. ff(k,i) /= miss_p)then
2380                   xff = ff(k,i)
2381                   xdd = dd(k,i)
2382                   xlon = lon(i)
2383                   CALL DA_FFDDUV(xff,xdd,platform%each(j)%u%inv, &
2384                               platform%each(j)%v%inv,xlon,1)
2385                   platform % each(j) % u % qc = qc_good
2386                   platform % each(j) % v % qc = qc_good
2387                   anyvar = .true.
2388                else
2389                   platform % each(j) % u % inv = missing_r
2390                   platform % each(j) % u % qc = missing_data
2391                   platform % each(j) % v % inv = missing_r
2392                   platform % each(j) % v % qc = missing_data
2393                endif
2394
2395                if(t(k,i) /= miss_p)then
2396                   platform % each(j) % t % inv = t(k,i)
2397                   platform % each(j) % t % qc = qc_good
2398                   anyvar = .true.
2399                else
2400                   platform % each(j) % t % inv = missing_r
2401                   platform % each(j) % t % qc = missing_data
2402                endif
2403
2404                if(q(k,i) /= miss_p)then
2405                   platform % each(j) % q % inv = q(k,i)
2406                   platform % each(j) % q % qc = qc_good
2407                   anyvar = .true.
2408                else
2409                   platform % each(j) % q % inv = missing_r
2410                   platform % each(j) % q % qc = missing_data
2411                endif
2412
2413 ! If we don't have any variables for profilers, punt.  (For the
2414 ! other datasets, we can't be here without having at least P or Z,
2415 ! which for the non-profiler datasets, must be actual observations.
2416 ! For profilers, the P/Z variables are only coordinates of the
2417 ! observations, and aren't really observations.)
2418
2419                if(anyvar .or. fm /= 35)then    !***** change 35 when profilers defined
2420
2421 ! Fill in the observation errors for this level.
2422
2423                   CALL GET_OB_ERRORS(platform%each(j))
2424
2425 ! Compute the (x,y,z) and increment the counter of levels that are used.
2426
2427                   CALL da_llxy(platform%info,platform%loc,outside,outside_all)
2428
2429                   if(fm == 35 .and. k == sfclev(i))sfclev_j = j
2430                   j = j + 1
2431
2432                else
2433
2434 ! Subtract this level from the total count for this station.
2435
2436                   platform % info % levels = platform % info % levels - 1
2437
2438                endif  ! If any observed variables exist
2439
2440             else
2441
2442                platform % info % levels = platform % info % levels - 1
2443
2444             endif     ! If P/Z exist
2445
2446          enddo        ! Level loop
2447
2448 ! It's possible that there will be no data at this point.  One example of this
2449 ! is a RASS (temperature) only Multi-Agency Profiler report.
2450
2451          if(platform % info % levels == 0) go to 5
2452
2453 ! Handle surface corrections, and other end-of-station processing.
2454
2455          CALL da_obs_proc_station(platform, fm)
2456
2457 ! Move the data for this station from PLATFORM to IV.
2458
2459          if (global .and. (platform%loc%i < ids .or. platform%loc%i >= ide)) then
2460             ndup = 2
2461          else
2462             ndup = 1
2463          endif
2464
2465 !------------------------------------------------------------------------
2466 ! Duplication explanation:
2467 !
2468 ! - Assume two reports (1 and 2)
2469 !
2470 ! - What this does is that if you have N obs available what goes into
2471 !   iv is 1,1,...,N,N.
2472 !
2473 !   Also note that the ntotal that's put into iv%info%obs_global_index
2474 !   is the "incremented of first pass only" local value, so:
2475 !
2476 !   if ndup = 2 then
2477 !
2478 !   index    orig ob index       iv%info%obs_global_index
2479 !     1           1                  1
2480 !     2           1                  1
2481 !     3           2                  2
2482 !     4           2                  2
2483 !
2484 ! It's a bit more complicated than that, because of the test on
2485 ! platform%loc%i.  This means we're only duplicating obs that are in
2486 ! a particular subset of the domain.
2487 !------------------------------------------------------------------------
2488
2489          do l = 1, ndup
2490
2491           SELECT CASE(fm)
2492
2493             CASE(15)
2494
2495                n = n + 1
2496                if(l.eq.1)ntot = ntot + 1
2497                j = n
2498                iv % metar(j) % h = platform % each(1) % height
2499                iv % info(metar) % zk(1,j) = platform % each(1) % zk
2500                iv % metar(j) % u = platform % each(1) % u
2501                iv % metar(j) % v = platform % each(1) % v
2502                iv % metar(j) % t = platform % each(1) % t
2503                iv % metar(j) % p = platform % each(1) % p
2504                iv % metar(j) % q = platform % each(1) % q
2505
2506             CASE(13)
2507
2508                n = n + 1
2509                if(l.eq.1)ntot = ntot + 1
2510                j = n
2511                iv % ships(j) % h = platform % each(1) % height
2512                iv % info(ships) % zk(1,j) = platform % each(1) % zk
2513                iv % ships(j) % u = platform % each(1) % u
2514                iv % ships(j) % v = platform % each(1) % v
2515                iv % ships(j) % t = platform % each(1) % t
2516                iv % ships(j) % p = platform % each(1) % p
2517                iv % ships(j) % q = platform % each(1) % q
2518
2519             CASE(88)
2520
2521                n = n + 1
2522                if(l.eq.1)ntot = ntot + 1
2523                j = n
2524
2525 ! With the datasets that have levels, we have to allocate the memory.
2526 ! Note that this is needed for GEOAMV and AIREP even though the MADIS data
2527 ! will produce a single-level report for each "station" (as contrasted
2528 ! with multiple-level reports for a single satellite or aircraft).  For example,
2529 ! if you have 10 aircraft reporting 5 levels of info each, this
2530 ! will produce 50 distinct reports.
2531
2532                allocate (iv % geoamv(j) % u(1:1))
2533                allocate (iv % geoamv(j) % v(1:1))
2534                allocate (iv % geoamv(j) % p(1:1))
2535
2536                iv % info(geoamv) % zk(1,j) = platform % each(1) % zk
2537                iv % geoamv(j) % u(1) = platform % each(1) % u
2538                iv % geoamv(j) % v(1) = platform % each(1) % v
2539                iv % geoamv(j) % p(1) = platform % each(1) % p % inv
2540
2541             CASE(96)
2542
2543                n = n + 1
2544                if(l.eq.1)ntot = ntot + 1
2545                j = n
2546
2547                allocate (iv % airep(j) % h(1:1))
2548                allocate (iv % airep(j) % u(1:1))
2549                allocate (iv % airep(j) % v(1:1))
2550                allocate (iv % airep(j) % t(1:1))
2551                allocate (iv % airep(j) % p(1:1))
2552
2553                iv % airep(j) % h(1) = platform % each(1) % height
2554                iv % info(airep) % zk(1,j) = platform % each(1) % zk
2555                iv % airep(j) % u(1) = platform % each(1) % u
2556                iv % airep(j) % v(1) = platform % each(1) % v
2557                iv % airep(j) % t(1) = platform % each(1) % t
2558                iv % airep(j) % p(1) = platform % each(1) % p % inv
2559
2560             CASE(35)
2561
2562 ! If this is a raob sfc level, put it into sonde_sfc and not into sound.
2563
2564                n = n + 1
2565                if(l.eq.1)ntot = ntot + 1
2566                j = n
2567
2568                m = platform%info%levels
2569
2570                if(sfclev_j /= 0)then
2571
2572                   iv%sonde_sfc(j)%h = platform%each(sfclev_j)%height
2573                   iv%sonde_sfc(j)%u = platform%each(sfclev_j)%u
2574                   iv%sonde_sfc(j)%v = platform%each(sfclev_j)%v
2575                   iv%sonde_sfc(j)%t = platform%each(sfclev_j)%t
2576                   iv%sonde_sfc(j)%q = platform%each(sfclev_j)%q
2577                   iv%sonde_sfc(j)%p = platform%each(sfclev_j)%p
2578                   iv%info(sonde_sfc)%zk(1,j) = platform%each(sfclev_j)%zk
2579
2580                   platform%info%levels = platform%info%levels - 1
2581
2582                else
2583
2584                   iv % sonde_sfc(j) % h = missing_r
2585                   iv % sonde_sfc(j) % u % inv = missing_r
2586                   iv % sonde_sfc(j) % u % qc = missing_data
2587                   iv % sonde_sfc(j) % u % error = missing_r
2588                   iv % sonde_sfc(j) % v = iv % sonde_sfc(j) % u
2589                   iv % sonde_sfc(j) % t = iv % sonde_sfc(j) % u
2590                   iv % sonde_sfc(j) % p = iv % sonde_sfc(j) % u
2591                   iv % sonde_sfc(j) % q = iv % sonde_sfc(j) % u
2592
2593                endif
2594
2595                allocate (iv % sound(j) % h        (1:platform % info % levels))
2596                allocate (iv % sound(j) % u        (1:platform % info % levels))
2597                allocate (iv % sound(j) % v        (1:platform % info % levels))
2598                allocate (iv % sound(j) % t        (1:platform % info % levels))
2599                allocate (iv % sound(j) % p        (1:platform % info % levels))
2600                allocate (iv % sound(j) % q        (1:platform % info % levels))
2601
2602                do k = 1, m
2603
2604 ! Note:  I don't have a testcase that involves the duplication logic,
2605 ! and I'm not sure if the "j-1" logic below will work correctly under
2606 ! those circumstances.
2607
2608                   if(sfclev_j.eq.0)then
2609                      iv % sound(j) % h(k) = platform % each(k) % height
2610                      iv % info(sound) % zk(k,j) = platform % each(k) % zk
2611                      iv % sound(j) % u(k) = platform % each(k) % u
2612                      iv % sound(j) % v(k) = platform % each(k) % v
2613                      iv % sound(j) % t(k) = platform % each(k) % t
2614                      iv % sound(j) % p(k) = platform % each(k) % p % inv
2615                      iv % sound(j) % q(k) = platform % each(k) % q
2616                   else if(k.lt.sfclev_j)then
2617                      iv % sound(j) % h(k) = platform % each(k) % height
2618                      iv % info(sound) % zk(k,j) = platform % each(k) % zk
2619                      iv % sound(j) % u(k) = platform % each(k) % u
2620                      iv % sound(j) % v(k) = platform % each(k) % v
2621                      iv % sound(j) % t(k) = platform % each(k) % t
2622                      iv % sound(j) % p(k) = platform % each(k) % p % inv
2623                      iv % sound(j) % q(k) = platform % each(k) % q
2624                   else if(k.ne.sfclev_j)then
2625                      iv % sound(j) % h(k-1) = platform % each(k) % height
2626                      iv % info(sound) % zk(k-1,j) = platform % each(k) % zk
2627                      iv % sound(j) % u(k-1) = platform % each(k) % u
2628                      iv % sound(j) % v(k-1) = platform % each(k) % v
2629                      iv % sound(j) % t(k-1) = platform % each(k) % t
2630                      iv % sound(j) % p(k-1) = platform % each(k) % p % inv
2631                      iv % sound(j) % q(k-1) = platform % each(k) % q
2632                   endif
2633
2634                enddo
2635
2636             CASE(111)
2637
2638                n = n + 1
2639                if(l.eq.1)ntot = ntot + 1
2640                j = n
2641                iv % gpspw(j) % tpw = platform % loc % pw
2642
2643           END SELECT
2644
2645           CALL FILL_IV_INFO(obs_index,n,platform,ntot,platform % info % levels,iv)
2646
2647          enddo      ! duplication loop
2648
2649       endif         ! if(keep(i) == 1)
2650
2651 5  enddo            ! station loop
2652
2653    nlocal_out = n
2654    ntotal_out = ntot
2655
2656 !zx   iv%info(obs_index)%plocal(1) = nlocal_out
2657 !zx   iv%info(obs_index)%ptotal(1) = ntotal(obs_index)
2658 !zx   iv%info(obs_index)%ntotal = ntotal(obs_index)
2659 !zx   if(fm == 35)then
2660 !zx      iv%info(sonde_sfc)%obs_global_index(obs_index) = ntotal_out
2661 !zx      iv%info(sonde_sfc)%plocal(1) = nlocal_out
2662 !zx      iv%info(sonde_sfc)%ptotal(1) = ntotal(sonde_sfc)
2663 !zx      iv%info(sonde_sfc)%ntotal = ntotal(sonde_sfc)
2664 !zx      iv%info(sonde_sfc)%nlocal = nlocal_out
2665 !zx   endif
2666
2667    return
2668
2669 END SUBROUTINE PROCESS_DATASET
2670
2671
2672 SUBROUTINE GET_OB_ERRORS(each)
2673
2674 !------------------------------------------------------------------------------
2675 ! PURPOSE: Fill in the observational errors,
2676 !
2677 ! METHOD:  The method was adapted from the MM5 3DVAR obs preprocessor
2678 !          obs_err_ncep subroutine.  The observational error is defined at
2679 !          1000, 700, 500, 300, 100, 50mb for U, V, T, RH and P.  Values at
2680 !          actual obs position are obtained by interpolation between
2681 !          these levels (logarithmic for U and V) and linear for (T and RH).
2682 !          Errors are extrapolated below 1000mb and above 50mb.
2683 !
2684 !          Also note that the RH error is put into the Q error -- the
2685 !          WRF 3DVAR moisture variable is Q, but the error table is
2686 !          specified in RH.  The da_fill_obs_structures routine (called
2687 !          after this on) will take the RH error in the Q slot, and
2688 !          use that in calculating the actual error to be used for Q.
2689 !
2690 ! REFERENCES:
2691 !
2692 !          Parrish S. and J. Derber 1992: The National Meteorological Centers
2693 !             Spectral Statistical Analysis-Interpolation System.
2694 !             Month. Wea. Rev., 120,1747-1763.
2695 !
2696 ! HISTORY: 07-22-02  M. Barth         Original version.
2697 !
2698 ! PARENT_MODULE: da_setup_structures
2699 !------------------------------------------------------------------------------
2700
2701    IMPLICIT NONE
2702
2703 ! ARGUMENTS:
2704
2705    TYPE (each_level_type), INTENT(INOUT) :: each          ! Data for this level.
2706
2707 ! LOCAL:
2708
2709    real err_k(6),err_u(6),err_v(6),err_t(6),err_rh(6),err_p(6)
2710
2711 !...NCEP errors (U in m/s, V in m/s, T in k, P in Pa)
2712 !...pressure is obtained from surface value linearly decreasing
2713 !...to 10% of its value at 10hPa
2714 !...RH has been divided by 2
2715 !
2716    data err_k/100000.,70000.,50000.,30000.,10000.,5000./
2717    data err_u/1.4,   2.4,   2.8,   3.4,   2.5,  2.7/
2718    data err_v/1.4,   2.4,   2.8,   3.4,   2.5,  2.7/
2719    data err_t/1.8,   1.3,   1.3,   2.0,   3.1,  4.0/
2720    data err_rh/10.0,  10.0,  10.0,  10.0,  10.0,  10.0/
2721    data err_p/100.0,  72.7,  54.5,  36.4,  18.2,  13.6/
2722
2723 !------------------------------------------------------------------------------!
2724 ! At this point, I don't think we can have missing P.  To stay on the safe side,
2725 ! if P is missing --> set all errors to missing.
2726
2727    if(each % p % inv == missing_r)then
2728
2729       each % u % error = missing_r
2730       each % v % error = missing_r
2731       each % p % error = missing_r
2732       each % t % error = missing_r
2733       each % q % error = missing_r
2734
2735    else
2736
2737 ! Vertical interpolation of NCEP observational error.
2738
2739       each % u % error = INTPLIN(each%p%inv,err_k,err_u)
2740       each % v % error = INTPLIN(each%p%inv,err_k,err_v)
2741       each % p % error = INTPLIN(each%p%inv,err_k,err_p)
2742       each % t % error = INTPLOG(each%p%inv,err_k,err_t)
2743       each % q % error = INTPLOG(each%p%inv,err_k,err_rh)
2744
2745     endif
2746
2747     return
2748
2749 END SUBROUTINE GET_OB_ERRORS
2750
2751 FUNCTION INTPLIN (x,xx,yy) RESULT (val)
2752 !------------------------------------------------------------------------------!
2753    implicit none
2754
2755    real, dimension (:) :: xx, yy
2756    real                :: x
2757    real                :: val
2758
2759    integer             :: n,m,jl
2760 !------------------------------------------------------------------------------!
2761    n = size (xx)
2762    m = size (yy)
2763
2764    jl = LOCATE(x,xx)
2765
2766    if (jl .le. 0) then
2767       val = yy (1)
2768    else if (jl .ge. n) then
2769       val = yy (n)
2770    else
2771       val = (xx (jl+1) - x) * yy (jl) + (x - xx (jl)) * yy (jl+1)
2772       val = val / (xx (jl+1) - xx (jl))
2773    endif
2774
2775 END FUNCTION INTPLIN
2776
2777 FUNCTION INTPLOG (x,xx,yy) result (val)
2778 !------------------------------------------------------------------------------!
2779    implicit none
2780
2781    real, dimension (:) :: xx, yy
2782    real                :: x
2783    real                :: val
2784
2785    integer             :: n,m,jl
2786 !------------------------------------------------------------------------------!
2787    n = size (xx)
2788    m = size (yy)
2789
2790    jl = LOCATE(x,xx)
2791
2792    if (jl .le. 0) then
2793       val = yy (1)
2794    else if (jl .ge. n) then
2795       val = yy (n)
2796    else
2797       val = log (xx (jl+1) / x) * yy (jl) + log (x / xx (jl)) * yy (jl+1)
2798       val = val / log (xx (jl+1) / xx (jl))
2799    endif
2800
2801 END FUNCTION INTPLOG
2802
2803 FUNCTION LOCATE (x,xx) result (index)
2804 !------------------------------------------------------------------------------!
2805    implicit none
2806
2807    real, dimension (:) :: xx
2808    real                :: x
2809    integer             :: index
2810
2811    integer             :: n,jl,jm,ju
2812    logical             :: ascnd
2813 !------------------------------------------------------------------------------!
2814    n = size (xx)
2815    ascnd = (xx (n) >= xx (1))   ! true if ascending order, false otherwise
2816    jl = 0                       ! initialize lower limit
2817    ju = n+1                     ! initialize upper limit
2818
2819    do
2820
2821       if (ju-jl <= 1) exit      ! repeat until this condition is satisfied
2822
2823       jm = (ju+jl) / 2.         ! compute a mid point
2824
2825       if (ascnd .eqv. (x >= xx (jm))) then
2826          jl = jm                ! replace mid point by lower limit
2827       else
2828          ju = jm                ! replace mid point by upper limit
2829       endif
2830
2831    enddo
2832
2833    if (x .eq. xx (1)) then      ! set the output, being careful with endpoints
2834       index = 1
2835    else if (x .eq. xx (n)) then
2836       index = n-1
2837    else
2838       index = jl
2839    endif
2840
2841 END FUNCTION LOCATE
2842
2843 SUBROUTINE FIND_SFC_LEVEL(ista,nlev,nsta,nlevels,levtype,p,slp,lev)
2844 !------------------------------------------------------------------------------!
2845    IMPLICIT NONE
2846
2847 ! ARGUMENTS:
2848
2849    INTEGER*4, INTENT(IN)        :: ista               ! Index of current station
2850    INTEGER*4, INTENT(IN)        :: nlev               ! Number of levels dimension
2851    INTEGER*4, INTENT(IN)        :: nsta               ! Number of stations dimension
2852    INTEGER*4, INTENT(IN)        :: nlevels            ! Number of levels for this station
2853    INTEGER*4                    :: lev                ! Index to sfc level
2854    REAL*4,    INTENT(IN)        :: levtype(nlev,nsta) ! Level type (1 is sfc)
2855    REAL*4,    INTENT(IN)        :: p(nlev,nsta)       ! Pressure (Pa)
2856    REAL*4,    INTENT(OUT)       :: slp                ! Surface pressure (Pa)
2857
2858 ! Note that lev is not dimensioned with INTENT(IN) due to some picky compilers...
2859
2860    do lev = 1, nlevels
2861
2862       if(levtype(lev,ista).eq.1)then
2863
2864          slp = p(lev,ista)
2865          return
2866
2867       endif
2868
2869    enddo
2870
2871    return
2872
2873 END SUBROUTINE FIND_SFC_LEVEL
2874
2875 SUBROUTINE FILL_IV_INFO(obs_index,iv_index,platform,ntotal,nlevels,iv)
2876
2877 !------------------------------------------------------------------------------
2878 ! PURPOSE: Fill in the iv%info structures used for each station's reports.
2879 !
2880 ! METHOD:  Take it out of platform and put it into iv%info.
2881 !
2882 ! HISTORY: 06-05-08  M. Barth         Original version.
2883
2884 !
2885 !
2886 ! PARENT_MODULE: da_setup_structures
2887 !------------------------------------------------------------------------------
2888
2889    IMPLICIT NONE
2890
2891 ! ARGUMENTS:
2892
2893    INTEGER*4, INTENT(IN)        :: obs_index       ! iv%info() index
2894    INTEGER*4, INTENT(IN)        :: iv_index        ! iv%info()%<item>() index
2895    TYPE (multi_level_type), INTENT(IN) :: platform ! Structure with the data
2896    INTEGER*4, INTENT(IN)        :: ntotal          ! NTOTAL at this point
2897    INTEGER*4, INTENT(IN)        :: nlevels         ! Num levels for the station
2898    TYPE (iv_type), INTENT(INOUT) :: iv             ! Innovation vector
2899
2900    iv%info(obs_index)%name(iv_index)      = platform%info%name
2901    iv%info(obs_index)%platform(iv_index)  = platform%info%platform
2902    iv%info(obs_index)%id(iv_index)        = platform%info%id
2903    iv%info(obs_index)%date_char(iv_index) = platform%info%date_char
2904
2905    iv%info(obs_index)%levels(iv_index)    = nlevels
2906    iv%info(obs_index)%lat(:,iv_index)     = platform%info%lat
2907    iv%info(obs_index)%lon(:,iv_index)     = platform%info%lon
2908    iv%info(obs_index)%elv(iv_index)       = platform%info%elv
2909    iv%info(obs_index)%pstar(iv_index)     = platform%info%pstar
2910
2911    iv%info(obs_index)%slp(iv_index)           = platform%loc%slp
2912    iv%info(obs_index)%pw(iv_index)            = platform%loc%pw
2913    iv%info(obs_index)%x(:,iv_index)           = platform%loc%x
2914    iv%info(obs_index)%y(:,iv_index)           = platform%loc%y
2915    iv%info(obs_index)%i(:,iv_index)           = platform%loc%i
2916    iv%info(obs_index)%j(:,iv_index)           = platform%loc%j
2917    iv%info(obs_index)%dx(:,iv_index)          = platform%loc%dx
2918    iv%info(obs_index)%dxm(:,iv_index)         = platform%loc%dxm
2919    iv%info(obs_index)%dy(:,iv_index)          = platform%loc%dy
2920    iv%info(obs_index)%dym(:,iv_index)         = platform%loc%dym
2921    iv%info(obs_index)%proc_domain(:,iv_index) = platform%loc%proc_domain
2922
2923    iv%info(obs_index)%obs_global_index(iv_index) = ntotal
2924
2925
2926 ! special case for sonde_sfc, duplicate sound info
2927
2928    if (obs_index == sound) then
2929
2930       iv%info(sonde_sfc)%name(iv_index)      = platform%info%name
2931       iv%info(sonde_sfc)%platform(iv_index)  = platform%info%platform
2932       iv%info(sonde_sfc)%id(iv_index)        = platform%info%id
2933       iv%info(sonde_sfc)%date_char(iv_index) = platform%info%date_char
2934       iv%info(sonde_sfc)%levels(iv_index)    = 1
2935       iv%info(sonde_sfc)%lat(:,iv_index)     = platform%info%lat
2936       iv%info(sonde_sfc)%lon(:,iv_index)     = platform%info%lon
2937       iv%info(sonde_sfc)%elv(iv_index)       = platform%info%elv
2938       iv%info(sonde_sfc)%pstar(iv_index)     = platform%info%pstar
2939
2940       iv%info(sonde_sfc)%slp(iv_index)           = platform%loc%slp
2941       iv%info(sonde_sfc)%pw(iv_index)            = platform%loc%pw
2942       iv%info(sonde_sfc)%x(:,iv_index)           = platform%loc%x
2943       iv%info(sonde_sfc)%y(:,iv_index)           = platform%loc%y
2944       iv%info(sonde_sfc)%i(:,iv_index)           = platform%loc%i
2945       iv%info(sonde_sfc)%j(:,iv_index)           = platform%loc%j
2946       iv%info(sonde_sfc)%dx(:,iv_index)          = platform%loc%dx
2947       iv%info(sonde_sfc)%dxm(:,iv_index)         = platform%loc%dxm
2948       iv%info(sonde_sfc)%dy(:,iv_index)          = platform%loc%dy
2949       iv%info(sonde_sfc)%dym(:,iv_index)         = platform%loc%dym
2950       iv%info(sonde_sfc)%proc_domain(:,iv_index) = platform%loc%proc_domain
2951
2952       iv%info(sonde_sfc)%obs_global_index(iv_index) = ntotal
2953
2954    end if
2955
2956    return
2957
2958 END SUBROUTINE FILL_IV_INFO
2959
2960 #endif