C$$$  MAIN PROGRAM DOCUMENTATION BLOCK
C                .      .    .                                       .
C MAIN PROGRAM:  MADIS_INGEST Sample hourly MADIS data assimilation   
C   PRGMMR: M. BARTH         ORG: FSL        DATE: 01-07-31
C
C ABSTRACT: This program provides a sample in how to use the MADIS 
C   library routines as they might be used in an ingest program for an 
C   hourly data assimilation system.  
C
C PROGRAM HISTORY LOG:
C   01-07-30  M. BARTH  Original version -- MADIS Version 2.0
C   02-05-06  M. BARTH  V2.1 -- Added MAP dataset and changed maximum
C                       number of surface station records from
C                       10000 to 25000.
C   03-03-13  M. BARTH  V2.2 - Changed maximum number of ACARS records
C                       from 10000 to 60000.
C   03-05-23  M. BARTH  V2.2 - Changed maximum number of surface records
C                       from 25000 to 50000.
C   03-07-03  M. BARTH  V2.3 - Added radiometer dataset.
C   03-10-02  M. BARTH  V2.4 - Added satellite winds dataset.
C   04-06-21  M. Barth  V2.5 - Changed comments re available datasets.
C   05-04-15  L. BENJAMIN V2.7 - Added satellite soundings dataset.
C                           - Changed maximum number of surface records
C                             from 50000 to 100000.
C                           - Increased maximum number of variable
C                             observation to 9,999,999.  This moves the
C                             count numbers to the right 2 spaces.
C   05-09-16  L. BENJAMIN V2.8 - Changed RAOB levels from 256 to 285 
C   05-10-27  L. BENJAMIN V2.9 - Changed maximum number of sounding
C                                records from 60000 to 80000
C                                Did not add SATRAD to this program,
C                                use SATS example if needed.
C   05-11-29  L. BENJAMIN V3.0 - Changed maximum number of sounding
C                                records from 60000 to 34000, about
C                                the  maxumum number of records in one 
C                                hour.  Changed RDMTR to 600 records,
C                                the maximun number of records in one 
C                                hour.
C   06-01-30  M. BARTH    V3.1 - Changed maximum number of surface
C                                records from 100000 to 200000.
C   06-02-13  L. BENJAMIN      - Added pressure to SATS variables.
C   06-01-30  L. BENJAMIN V3.2 - Changed version number.
C   09-08-14  L. BENJAMIN V3.8 - Initialize MINIT for every time loop
C                                instead of just at the beginning.
C                                Bug caught by MADIS users.  Thanks.  
C                                Changed maximum number of surface
C                                records from 200000 to 300000.
C                                Added wisdom dataset.
C   14-01-02  L. BENJAMIN V4.2 - Change max records from 300000 to
C                                600000. Changed satwnd from 60000 to 
C                                120000.
C
C USAGE:
C   INPUT FILES: (for the entire program)
C
C     MADIS initialization files in MADIS_STATIC directory   
C
C     MADIS observation files under MADIS_DATA directories
C
C   OUTPUT FILES:
C
C     ingest.txt       Text output listing the count of stations and
C                      variables read for each dataset.
C
C     standard output  Error messages
C
C SUBPROGRAMS CALLED:
C     LIBRARY: - MADISLIB:  MINIT, MSETWIN, MSETDOM, MSETQC, MTRNTIM, 
C                           MSFCSTA, MGETIJ, MGETSFC, MRAOBSTA,
C                           MGETRAOB, MNPNSTA, MGETNPN, MACARSSTA,
C                           MGETACARS, MMAPSTA, MGETMAP, MRDMTRSTA,
C                           MGETRDMTR, MSATWNDSTA, MGETSATWND,
C                           MSATSNDSTA, MGETSATSND
C
C REMARKS:
C
C   1) This program provides a sample in how to use the MADIS library
C      routines as they might be used in an ingest program for an 
C      hourly data assimilation system.  For more details on how to use
C      the various calls and options for each individual dataset,
C      refer to the following tutorial programs:
C
C      - sfcdump.f               Surface
C      - raobdump.f              Radiosonde
C      - npndump.f               NOAA Profiler Network 
C      - acarsdump.f             Aircraft data (ACARS, AMDARS, etc.)
C      - acarspdump.f            Aircraft data profiles at airports
C      - hydrodump.f             Hydro     
C      - mapdump.f               Multi-Agency Profilers
C      - rdmtrdump.f             Radiometer
C      - satwnddump.f            Satellite winds
C      - satsnddump.f            Satellite soundings
C      - satraddump.f            Satellite radiance
C      - snowdump.f              Snow dataset
C
C      (The aircraft data profiles & hydro datasets aren't used in this
C       program.)
C
C   2) The MADIS library can calculate a number of variables that aren't
C      actually observations stored in the database for each dataset 
C      (e.g., specific humidity isn't a reported quantity, ACARS 
C      pressures aren't reported).  This program demonstrates this by 
C      getting pressure and height for all datasets, and virtual 
C      temperature and specific humidity for all applicable datasets.
C      Fpr more information on what variables are available, see the
C      *variable_list.txt files in the MADIS doc directory.
C
C   3) See the INSTALL.* files in the MADIS doc directory for 
C      instructions on how to install and build the MADIS library and 
C      programs (INSTALL.unix for Unix/Linux, INSTALL.windows for
C      Windows).
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 77 + EXTENSIONS
C   MACHINE: UNIX
C
C$$$
C-----------------------------------------------------------------------
C     Declarations notes:
C
C     Key constants:
C
C       - 0              successful status return
C       - 999999.0       missing data flag
C       - 600000         maximum number of surface stations
C       - 1000           maximum number of radiosonde stations
C       - 40             maximum number of NPN stations
C       - 60000          maximum number of ACARS reports
C       - 150            maximum number of MAP stations
C       - 600            maximum number of radiometer stations
C       - 12000          maximum number of satellite wind reports

      implicit none

      integer success_p
      real miss_p
      parameter (success_p = 0)
      parameter (miss_p = 999999.0)

      integer istatus,i,j,minbck,minfwd,recwin,qctype,nobs,nhr
      character*9 atime,new_time
      character*13 tstr
      real dummy

C     Surface data
C     -------------

C     - Station info

      integer nsta_sfc
      character*10 stanam_sfc(600000)
      character*9  timeob_sfc(600000)
      character*10 provdr_sfc(600000)
      integer       wmoid_sfc(600000)
      real   lat_sfc(600000)
      real   lon_sfc(600000)
      real  elev_sfc(600000)
      real    ri_sfc(600000)
      real    rj_sfc(600000)
      integer source_n(600000)

C     - Variables

      real  p_sfc(600000)
      real tv_sfc(600000)
      real  q_sfc(600000)
      real  u_sfc(600000)
      real  v_sfc(600000)

C     - QC information

      integer  p_sfc_qca(600000)
      integer tv_sfc_qca(600000)
      integer  q_sfc_qca(600000)
      integer  u_sfc_qca(600000)
      integer  v_sfc_qca(600000)

      integer  p_sfc_qcr(600000)
      integer tv_sfc_qcr(600000)
      integer  q_sfc_qcr(600000)
      integer  u_sfc_qcr(600000)
      integer  v_sfc_qcr(600000)

      character*1  p_sfc_qcd(600000)
      character*1 tv_sfc_qcd(600000)
      character*1  q_sfc_qcd(600000)
      character*1  u_sfc_qcd(600000)
      character*1  v_sfc_qcd(600000)

C     Radiosonde data
C     ---------------

C     - Station info

      integer nsta_raob
      character*10 stanam_raob(1000)
      character*9 timeob_raob(1000)
      integer wmoid_raob(1000)
      integer nlevels_raob(1000)
      real lat_raob(1000)
      real lon_raob(1000)
      real elev_raob(1000)
      real ri_raob(1000)
      real rj_raob(1000)

C     - Variables

      real p_raob(285,1000)
      real ht_raob(285,1000)
      real tv_raob(285,1000)
      real q_raob(285,1000)
      real u_raob(285,1000)
      real v_raob(285,1000)

C     - QC information

      integer p_raob_qca(285,1000)
      integer ht_raob_qca(285,1000)
      integer tv_raob_qca(285,1000)
      integer q_raob_qca(285,1000)
      integer u_raob_qca(285,1000)
      integer v_raob_qca(285,1000)

      integer p_raob_qcr(285,1000)
      integer ht_raob_qcr(285,1000)
      integer tv_raob_qcr(285,1000)
      integer q_raob_qcr(285,1000)
      integer u_raob_qcr(285,1000)
      integer v_raob_qcr(285,1000)

      character*1 p_raob_qcd(285,1000)
      character*1 ht_raob_qcd(285,1000)
      character*1 tv_raob_qcd(285,1000)
      character*1 q_raob_qcd(285,1000)
      character*1 u_raob_qcd(285,1000)
      character*1 v_raob_qcd(285,1000)

C     NPN data
C     --------

C     - Station info

      integer nsta_npn
      character*10 stanam_npn(40)
      character*9 timeob_npn(40)
      integer wmoid_npn(40)
      integer nlevels_npn(40)
      real lat_npn(40)
      real lon_npn(40)
      real elev_npn(40)
      real ri_npn(40)
      real rj_npn(40)

C     - Variables

      real u_npn(64,40)
      real v_npn(64,40)
      real ht_npn(64,40)
      real p_npn(64,40)

C     - QC information

      integer u_npn_qca(64,40)
      integer v_npn_qca(64,40)
      integer ht_npn_qca(64,40)
      integer p_npn_qca(64,40)

      integer u_npn_qcr(64,40)
      integer v_npn_qcr(64,40)
      integer ht_npn_qcr(64,40)
      integer p_npn_qcr(64,40)

      character*1 u_npn_qcd(64,40)
      character*1 v_npn_qcd(64,40)
      character*1 ht_npn_qcd(64,40)
      character*1 p_npn_qcd(64,40)

C     ACARS data
C     ----------

C     - Station info

      integer nsta_acars
      character*10 stanam_acars(60000)
      character*9 timeob_acars(60000)
      real ri_acars(60000)
      real rj_acars(60000)

C     - Variables

      real ht_acars(60000)
      real p_acars(60000)
      real lat_acars(60000)
      real lon_acars(60000)
      real tv_acars(60000)
      real q_acars(60000)
      real u_acars(60000)
      real v_acars(60000)

C     - QC information

      integer ht_acars_qca(60000)
      integer p_acars_qca(60000)
      integer lat_acars_qca(60000)
      integer lon_acars_qca(60000)
      integer tv_acars_qca(60000)
      integer q_acars_qca(60000)
      integer u_acars_qca(60000)
      integer v_acars_qca(60000)

      integer ht_acars_qcr(60000)
      integer p_acars_qcr(60000)
      integer lat_acars_qcr(60000)
      integer lon_acars_qcr(60000)
      integer tv_acars_qcr(60000)
      integer q_acars_qcr(60000)
      integer u_acars_qcr(60000)
      integer v_acars_qcr(60000)

      character*1 ht_acars_qcd(60000)
      character*1 p_acars_qcd(60000)
      character*1 lat_acars_qcd(60000)
      character*1 lon_acars_qcd(60000)
      character*1 tv_acars_qcd(60000)
      character*1 q_acars_qcd(60000)
      character*1 u_acars_qcd(60000)
      character*1 v_acars_qcd(60000)

C     MAP data
C     --------

C     - Station info

      integer nsta_map
      character*10 stanam_map(150)
      character*9 timeob_map(150)
      character*10 provdr_map(150)
      integer wmoid_map(150)
      integer nlevels_map(150)
      real lat_map(150)
      real lon_map(150)
      real elev_map(150)
      real ri_map(150)
      real rj_map(150)

C     - Variables

      real u_map(202,150)
      real v_map(202,150)
      real ht_map(202,150)
      real p_map(202,150)
      real tv_map(202,150)

C     - QC information

      integer u_map_qca(202,150)
      integer v_map_qca(202,150)
      integer ht_map_qca(202,150)
      integer p_map_qca(202,150)
      integer tv_map_qca(202,150)

      integer u_map_qcr(202,150)
      integer v_map_qcr(202,150)
      integer ht_map_qcr(202,150)
      integer p_map_qcr(202,150)
      integer tv_map_qcr(202,150)

      character*1 u_map_qcd(202,150)
      character*1 v_map_qcd(202,150)
      character*1 ht_map_qcd(202,150)
      character*1 p_map_qcd(202,150)
      character*1 tv_map_qcd(202,150)

C     Radiometer data
C     ---------------

C     - Station info

      integer nsta_rdmtr
      character*10 stanam_rdmtr(600)
      character*9 timeob_rdmtr(600)
      character*10 provdr_rdmtr(600)
      integer wmoid_rdmtr(600)
      integer nlevels_rdmtr(600)
      real lat_rdmtr(600)
      real lon_rdmtr(600)
      real elev_rdmtr(600)
      real ri_rdmtr(600)
      real rj_rdmtr(600)

C     - Variables

      real p_rdmtr(100,600)
      real ht_rdmtr(100,600)
      real tv_rdmtr(100,600)
      real q_rdmtr(100,600)

C     - QC information

      integer p_rdmtr_qca(100,600)
      integer ht_rdmtr_qca(100,600)
      integer tv_rdmtr_qca(100,600)
      integer q_rdmtr_qca(100,600)

      integer p_rdmtr_qcr(100,600)
      integer ht_rdmtr_qcr(100,600)
      integer tv_rdmtr_qcr(100,600)
      integer q_rdmtr_qcr(100,600)

      character*1 p_rdmtr_qcd(100,600)
      character*1 ht_rdmtr_qcd(100,600)
      character*1 tv_rdmtr_qcd(100,600)
      character*1 q_rdmtr_qcd(100,600)

C     Satellite wind data
C     -------------------

C     - Station info

      integer nsta_satw
      character*10 provdr_satw(120000)
      character*9 timeob_satw(120000)
      real lat_satw(120000)
      real lon_satw(120000)
      real p_satw(120000)
      real ri_satw(120000)
      real rj_satw(120000)

C     - Variables

      real u_satw(120000)
      real v_satw(120000)

C     - QC information

      integer u_satw_qca(120000)
      integer v_satw_qca(120000)

      integer u_satw_qcr(120000)
      integer v_satw_qcr(120000)

      character*1 u_satw_qcd(120000)
      character*1 v_satw_qcd(120000)

C     Satellite sounding data
C     -------------------

C     - Station info

      integer nsta_sats
      character*10 provdr_sats(35000)
      character*9 timeob_sats(35000)
      integer nlevels_sats(35000)
      real lat_sats(35000)
      real lon_sats(35000)
      real p_sats(35000)
      real ri_sats(35000)
      real rj_sats(35000)


C     - Variables

      real t_sats(40,35000)
      real mr_sats(40,35000)
      real pr_sats(40,35000)

C     - QC information

      integer t_sats_qca(40,35000)
      integer mr_sats_qca(40,35000)
      integer pr_sats_qca(40,35000)

      integer t_sats_qcr(40,35000)
      integer mr_sats_qcr(40,35000)
      integer pr_sats_qcr(40,35000)

      character*1 t_sats_qcd(40,35000)
      character*1 mr_sats_qcd(40,35000)
      character*1 pr_sats_qcd(40,35000)


C     WISDOM data
C     -------------------

C     - Station info

      integer nsta_wisd
      character*10 stanam_wisd(50000)
      character*10 provdr_wisd(50000)
      real lat_wisd(50000)
      real lon_wisd(50000)
      real elev_wisd(50000)
      character*9 timeob_wisd(50000)
      real ri_wisd(50000)
      real rj_wisd(50000)

C     - Variables

      real u_wisd(50000)
      real v_wisd(50000)

C     - QC information

      integer u_wisd_qca(50000)
      integer v_wisd_qca(50000)

      integer u_wisd_qcr(50000)
      integer v_wisd_qcr(50000)

      character*1 u_wisd_qcd(50000)
      character*1 v_wisd_qcd(50000)



C-----------------------------------------------------------------------
C PROGRAM SETUP
C -------------
C
C     Open an output file for this program.

      open(unit=1,file='ingest.txt',status='UNKNOWN')

C-----------------------------------------------------------------------
C MANDATORY MADIS INITIALIZATION
C -----------------------------

C     Initialize the MADIS subsystems to be used, using the FSL
C     database, and with error messages echoed to standard output.

C-----------------------------------------------------------------------
C OPTIONAL MADIS INITIALIZATION
C ----------------------------
C
C     Set the horizontal domain for the CONUS 212 Lambert conformal
C     conic projection.

C      CALL MSETDOM(3,40635.25,12.19,-133.459,-95.0,25.0,25.0,dummy,
C     1             185,129,istatus)
C      if(istatus.ne.success_p)stop

C     Ask the user what time range he wants to read & inventory.

      write(*,1)
 1    format(' Start time (YYYYMMDD_HHMM):',$)
      read(*,2)tstr
 2    format(a)
      write(*,3)
 3    format(' Number of hours to process:',$)
      read(*,*)nhr

C     Translate the user-specified time string to our internal
C     time string.

      CALL MTRNTIM(tstr,1,atime,istatus)
      if(istatus.ne.success_p)stop


C-----------------------------------------------------------------------
C LOOP THROUGH HOURS 
C ------------------

      do i = 1, nhr

C     Get the surface data for this hour
C     ----------------------------------


         write(1,4)'Inventory of surface data for ',atime
 4       format(/,1x,a,a,/)

         CALL MINIT('SFC','FSL',.true.,istatus)
         if(istatus.ne.success_p)stop

         CALL MSETDOM(3,40635.25,12.19,-133.459,-95.0,25.0,
     1        25.0,dummy,185,129,istatus)
         if(istatus.ne.success_p)stop

         qctype = 99
         CALL MSETQC(qctype,istatus)
         if(istatus.ne.success_p)stop

C     Set the time window to return only one record per station
C     within the specified time window (-30 minutes --> nominal
C     time --> +29 minutes).  If a station has multiple reports
C     within the time window, only the latest one within the
C     window will be returned (RECWIN = 3).  Note that we use +29 
C     minutes and not +30 so that we can't possibly get the same 
C     obs for two different hours in a row -- if +- 30 minutes were 
C     used, then you'd get 11:30 obs when the nominal time is 11:00 
C     and again when it's 12:00.  (Then again, that might be OK for
C     your application.)

         minbck = -30
         minfwd = 29
         recwin = 3
         CALL MSETWIN(minbck,minfwd,recwin,istatus)
         if(istatus.ne.success_p)stop

C     Load in the stations with data for this hour.

         CALL MSFCSTA(atime,nsta_sfc,stanam_sfc,wmoid_sfc,lat_sfc,
     1        lon_sfc,elev_sfc,timeob_sfc,provdr_sfc,istatus,
     1        source_n)

         write(1,5)'number of sfc stations                    ',nsta_sfc
 5       format(1x,a,i7)
         if(istatus.ne.success_p)go to 100

C     Get the floating point grid indices into the CONUS 212 grid.
C     Note that this works for grids laid out this way:

C
C                        NE (I=NX,J=NY)
C
C     SW (I=1,J=1)

         do j = 1, nsta_sfc
            CALL MGETIJ(lat_sfc(j),lon_sfc(j),ri_sfc(j),rj_sfc(j),
     1                  istatus)


         enddo

C     Read the variables in this order:  pressure, virtual temperature, 
C     specific humidity, u, v.  

         CALL MGETSFC(atime,'P',nsta_sfc,nobs,p_sfc,
     1                p_sfc_qcd,p_sfc_qca,p_sfc_qcr,istatus)
         if(istatus.ne.success_p)go to 100
         write(1,5)'number of non-missing pressure            ',nobs

         CALL MGETSFC(atime,'TV',nsta_sfc,nobs,tv_sfc,
     1                tv_sfc_qcd,tv_sfc_qca,tv_sfc_qcr,istatus)
         if(istatus.ne.success_p)go to 100
         write(1,5)'number of non-missing virtual temperature ',nobs

         CALL MGETSFC(atime,'Q',nsta_sfc,nobs,q_sfc,
     1                q_sfc_qcd,q_sfc_qca,q_sfc_qcr,istatus)
         if(istatus.ne.success_p)go to 100
         write(1,5)'number of non-missing specific humidity   ',nobs

         CALL MGETSFC(atime,'U',nsta_sfc,nobs,u_sfc,
     1                u_sfc_qcd,u_sfc_qca,u_sfc_qcr,istatus)
         if(istatus.ne.success_p)go to 100
         write(1,5)'number of non-missing u                   ',nobs

         CALL MGETSFC(atime,'V',nsta_sfc,nobs,v_sfc,
     1                v_sfc_qcd,v_sfc_qca,v_sfc_qcr,istatus)
         if(istatus.ne.success_p)go to 100
         write(1,5)'number of non-missing v                   ',nobs

         

C     Get the radiosonde data for this hour
C     -------------------------------------

 100     write(1,4)'Inventory of radiosonde data for ',atime
         
         CALL MINIT('RAOB','FSL',.true.,istatus)
         if(istatus.ne.success_p)stop

         CALL MSETDOM(3,40635.25,12.19,-133.459,-95.0,25.0,
     1        25.0,dummy,185,129,istatus)
         if(istatus.ne.success_p)stop

         qctype = 99
         CALL MSETQC(qctype,istatus)
         if(istatus.ne.success_p)stop

C     Set the time window to get radiosondes within +- 30 minutes
C     of the nominal time.  Note that the majority of raobs
C     will report at 0000 and 1200.  There will also be a number of
C     stations that report one hour off from those synoptic times 
C     (e.g., 1100 or 1300), and with this hourly cycle example, those 
C     will be returned on off-synoptic hours.  If you want to get *all* 
C     the 0000 and 1200 soundings at those hours, you should probably 
C     set the time window to +- 1 hour.  This isn't done in this
C     example because then you'd have to take into account the case of
C     reading the 0100 raobs at both 0000 and 0100.

         minbck = -30
         minfwd = 30
         recwin = 1
         CALL MSETWIN(minbck,minfwd,recwin,istatus)
         if(istatus.ne.success_p)stop

C     Load in the stations with data for this hour.

         CALL MRAOBSTA(atime,nsta_raob,stanam_raob,wmoid_raob,lat_raob,
     1                 lon_raob,elev_raob,timeob_raob,istatus)

         write(1,5)'number of radiosonde stations             ',
     1             nsta_raob
         if(istatus.ne.success_p)go to 200

C     Get the floating point grid indices into the CONUS 212 grid.

         do j = 1, nsta_raob
            CALL MGETIJ(lat_raob(j),lon_raob(j),ri_raob(j),rj_raob(j),
     1                  istatus)
         enddo

C     Read the variables in this order:  pressure, height, virtual 
C     temperature, specific humidity, u, v.

         CALL MGETRAOB(atime,'P',nsta_raob,nobs,nlevels_raob,p_raob,
     1                 p_raob_qcd,p_raob_qca,p_raob_qcr,istatus)
         if(istatus.ne.success_p)go to 200
         write(1,5)'number of non-missing pressure            ',nobs

         CALL MGETRAOB(atime,'HT',nsta_raob,nobs,nlevels_raob,ht_raob,
     1                 ht_raob_qcd,ht_raob_qca,ht_raob_qcr,istatus)
         if(istatus.ne.success_p)go to 200
         write(1,5)'number of non-missing height              ',nobs

         CALL MGETRAOB(atime,'TV',nsta_raob,nobs,nlevels_raob,tv_raob,
     1                 tv_raob_qcd,tv_raob_qca,tv_raob_qcr,istatus)
         if(istatus.ne.success_p)go to 200
         write(1,5)'number of non-missing virtual temperature ',nobs

         CALL MGETRAOB(atime,'Q',nsta_raob,nobs,nlevels_raob,q_raob,
     1                 q_raob_qcd,q_raob_qca,q_raob_qcr,istatus)
         if(istatus.ne.success_p)go to 200
         write(1,5)'number of non-missing specific humidity   ',nobs

         CALL MGETRAOB(atime,'U',nsta_raob,nobs,nlevels_raob,u_raob,
     1                 u_raob_qcd,u_raob_qca,u_raob_qcr,istatus)
         if(istatus.ne.success_p)go to 200
         write(1,5)'number of non-missing u                   ',nobs

         CALL MGETRAOB(atime,'V',nsta_raob,nobs,nlevels_raob,v_raob,
     1                 v_raob_qcd,v_raob_qca,v_raob_qcr,istatus)
         if(istatus.ne.success_p)go to 200
         write(1,5)'number of non-missing v                   ',nobs


C     Get the NPN data for this hour
C     ------------------------------

 200     write(1,4)'Inventory of NPN data for ',atime

         CALL MINIT('NPN','FSL',.true.,istatus)
         if(istatus.ne.success_p)stop

         CALL MSETDOM(3,40635.25,12.19,-133.459,-95.0,25.0,
     1        25.0,dummy,185,129,istatus)
         if(istatus.ne.success_p)stop

         qctype = 99
         CALL MSETQC(qctype,istatus)
         if(istatus.ne.success_p)stop

C     Set the time window to return all profiler stations in
C     the hourly file without any windowing logic.  This default
C     window is good enough for an hourly cycle as all profilers
C     are time-stamped at the top of the hour and there are
C     no duplicates.

         minbck = 0
         minfwd = 0
         recwin = 0
         CALL MSETWIN(minbck,minfwd,recwin,istatus)
         if(istatus.ne.success_p)stop

C     Load in the stations with data for this hour.

         CALL MNPNSTA(atime,nsta_npn,stanam_npn,wmoid_npn,lat_npn,
     1                lon_npn,elev_npn,timeob_npn,istatus)

         write(1,5)'number of NPN stations                    ',
     1             nsta_npn
         if(istatus.ne.success_p)go to 300

C     Get the floating point grid indices into the CONUS 212 grid.

         do j = 1, nsta_npn
            CALL MGETIJ(lat_npn(j),lon_npn(j),ri_npn(j),rj_npn(j),
     1                  istatus)
         enddo

C     Read the variables in this order:  u, v, height, pressure.

         CALL MGETNPN(atime,'U',nsta_npn,nobs,nlevels_npn,
     1                u_npn,u_npn_qcd,u_npn_qca,u_npn_qcr,istatus)
         if(istatus.ne.success_p)go to 300
         write(1,5)'number of non-missing u                   ',nobs

         CALL MGETNPN(atime,'V',nsta_npn,nobs,nlevels_npn,
     1                v_npn,v_npn_qcd,v_npn_qca,v_npn_qcr,istatus)
         if(istatus.ne.success_p)go to 300
         write(1,5)'number of non-missing v                   ',nobs

         CALL MGETNPN(atime,'HT',nsta_npn,nobs,nlevels_npn,
     1                ht_npn,ht_npn_qcd,ht_npn_qca,ht_npn_qcr,istatus)
         if(istatus.ne.success_p)go to 300
         write(1,5)'number of non-missing height              ',nobs

         CALL MGETNPN(atime,'P',nsta_npn,nobs,nlevels_npn,
     1                p_npn,p_npn_qcd,p_npn_qca,p_npn_qcr,istatus)
         if(istatus.ne.success_p)go to 300
         write(1,5)'number of non-missing pressure            ',nobs

C     Get the ACARS data for this hour
C     ------------------------------

 300     write(1,4)'Inventory of ACARS data for ',atime

         CALL MINIT('ACARS','FSL',.true.,istatus)
         if(istatus.ne.success_p)stop

         CALL MSETDOM(3,40635.25,12.19,-133.459,-95.0,25.0,
     1        25.0,dummy,185,129,istatus)
         if(istatus.ne.success_p)stop

         qctype = 99
         CALL MSETQC(qctype,istatus)
         if(istatus.ne.success_p)stop

C     Set the time window to get ACARS within +29 or -30 minutes
C     from the nominal time.  ACARS are generally reported all
C     around the hour, so you should set your time window to
C     whatever you'd really like to get.  
C
C     Also note that the RECWIN choice, used to specify how to handle
C     multiple reports from the same station within the time window,
C     should be set to 4 in order to return all records for each
C     station that's in the time window.  Each ACARS "station" 
C     (aircraft) will typically report several observations at flight 
C     level (spaced 7 minutes apart) or dozens of observations during 
C     ascents/descents (several reports per minute).  

         minbck = -30
         minfwd = 29
         recwin = 4
         CALL MSETWIN(minbck,minfwd,recwin,istatus)
         if(istatus.ne.success_p)stop

C     Load in the stations with data for this hour.

         CALL MACARSSTA(atime,nsta_acars,stanam_acars,timeob_acars,
     1                  istatus)

         write(1,5)'number of ACARS reports                   ',
     1             nsta_acars
         if(istatus.ne.success_p)go to 400

C     Read the latitude and longitude variables and their QC info.  
C     Note that these location coordinates are actually reported 
c     observations and have QC applied to them, unlike with most other 
C     MADIS data sets.

         CALL MGETACARS(atime,'LAT',nsta_acars,nobs,lat_acars,
     1                  lat_acars_qcd,lat_acars_qca,lat_acars_qcr,
     2                  istatus)
         if(istatus.ne.success_p)go to 400
         write(1,5)'number of non-missing lat                 ',nobs


         CALL MGETACARS(atime,'LON',nsta_acars,nobs,lon_acars,
     1                  lon_acars_qcd,lon_acars_qca,lon_acars_qcr,
     2                  istatus)
         if(istatus.ne.success_p)go to 400
         write(1,5)'number of non-missing lon                 ',nobs


C     Get the floating point grid indices into the CONUS 212 grid.

         do j = 1, nsta_acars
            CALL MGETIJ(lat_acars(j),lon_acars(j),ri_acars(j),
     1                  rj_acars(j),istatus)
         enddo

C     Read the variables in this order:  height, pressure, u, v, virtual
C     temperature, specific humidity.

         CALL MGETACARS(atime,'HT',nsta_acars,nobs,ht_acars,
     1                  ht_acars_qcd,ht_acars_qca,ht_acars_qcr,
     2                  istatus)
         if(istatus.ne.success_p)go to 400
         write(1,5)'number of non-missing height              ',nobs

         CALL MGETACARS(atime,'P',nsta_acars,nobs,p_acars,
     1                  p_acars_qcd,p_acars_qca,p_acars_qcr,
     2                  istatus)
         if(istatus.ne.success_p)go to 400
         write(1,5)'number of non-missing pressure            ',nobs

         CALL MGETACARS(atime,'U',nsta_acars,nobs,u_acars,
     1                  u_acars_qcd,u_acars_qca,u_acars_qcr,istatus)
         if(istatus.ne.success_p)go to 400
         write(1,5)'number of non-missing u                   ',nobs

         CALL MGETACARS(atime,'V',nsta_acars,nobs,v_acars,
     1                  v_acars_qcd,v_acars_qca,v_acars_qcr,istatus)
         if(istatus.ne.success_p)go to 400
         write(1,5)'number of non-missing v                   ',nobs

         CALL MGETACARS(atime,'TV',nsta_acars,nobs,tv_acars,
     1                  tv_acars_qcd,tv_acars_qca,tv_acars_qcr,istatus)
         if(istatus.ne.success_p)go to 400
         write(1,5)'number of non-missing virtual temperature ',nobs

         CALL MGETACARS(atime,'Q',nsta_acars,nobs,q_acars,
     1                  q_acars_qcd,q_acars_qca,q_acars_qcr,istatus)
         if(istatus.ne.success_p)go to 400
         write(1,5)'number of non-missing specific humidity   ',nobs

C     Get the MAP data for this hour
C     ------------------------------

 400     write(1,4)'Inventory of MAP data for ',atime

         CALL MINIT('MAP','FSL',.true.,istatus)
         if(istatus.ne.success_p)stop

         CALL MSETDOM(3,40635.25,12.19,-133.459,-95.0,25.0,
     1        25.0,dummy,185,129,istatus)
         if(istatus.ne.success_p)stop

         qctype = 99
         CALL MSETQC(qctype,istatus)
         if(istatus.ne.success_p)stop

C     Set the time window to get MAP data within +29 or -30 minutes
C     from the nominal time.  MAP stations are generally reported all
C     around the hour, so you should set your time window to
C     whatever you'd really like to get.  If a station has multiple 
C     reports within the time window, only the latest one within the
C     window will be returned (RECWIN = 3).  

         minbck = -30
         minfwd = 29
         recwin = 3
         CALL MSETWIN(minbck,minfwd,recwin,istatus)
         if(istatus.ne.success_p)stop

C     Load in the stations with data for this hour.

         CALL MMAPSTA(atime,nsta_map,stanam_map,wmoid_map,lat_map,
     1                lon_map,elev_map,timeob_map,provdr_map,istatus)

         write(1,5)'number of MAP stations                    ',
     1             nsta_map
         if(istatus.ne.success_p)go to 500

C     Get the floating point grid indices into the CONUS 212 grid.

         do j = 1, nsta_map
            CALL MGETIJ(lat_map(j),lon_map(j),ri_map(j),rj_map(j),
     1                  istatus)
         enddo

C     Read the variables in this order:  u, v, height, pressure,
C     virtual temperature.

         CALL MGETMAP(atime,'U',nsta_map,nobs,nlevels_map,
     1                u_map,u_map_qcd,u_map_qca,u_map_qcr,istatus)
         if(istatus.ne.success_p)go to 500
         write(1,5)'number of non-missing u                   ',nobs

         CALL MGETMAP(atime,'V',nsta_map,nobs,nlevels_map,
     1                v_map,v_map_qcd,v_map_qca,v_map_qcr,istatus)
         if(istatus.ne.success_p)go to 500
         write(1,5)'number of non-missing v                   ',nobs

         CALL MGETMAP(atime,'HT',nsta_map,nobs,nlevels_map,
     1                ht_map,ht_map_qcd,ht_map_qca,ht_map_qcr,istatus)
         if(istatus.ne.success_p)go to 500
         write(1,5)'number of non-missing height              ',nobs

         CALL MGETMAP(atime,'P',nsta_map,nobs,nlevels_map,
     1                p_map,p_map_qcd,p_map_qca,p_map_qcr,istatus)
         if(istatus.ne.success_p)go to 500
         write(1,5)'number of non-missing pressure            ',nobs

         CALL MGETMAP(atime,'TV',nsta_map,nobs,nlevels_map,
     1                tv_map,tv_map_qcd,tv_map_qca,tv_map_qcr,istatus)
         if(istatus.ne.success_p)go to 500
         write(1,5)'number of non-missing virtual temperature ',nobs

C     Get the radiometer data for this hour
C     -------------------------------------

 500     write(1,4)'Inventory of radiometer data for ',atime

         CALL MINIT('RDMTR','FSL',.true.,istatus)
         if(istatus.ne.success_p)stop

         CALL MSETDOM(3,40635.25,12.19,-133.459,-95.0,25.0,
     1        25.0,dummy,185,129,istatus)
         if(istatus.ne.success_p)stop

         qctype = 99
         CALL MSETQC(qctype,istatus)
         if(istatus.ne.success_p)stop

C     Set the time window to get radiometer data within +29 or -30 
C     minutes from the nominal time.  Radiometer stations are generally
C     reported all around the hour, so you should set your time window 
C     to whatever you'd really like to get.  If a station has multiple 
C     reports within the time window, only the latest one within the
C     window will be returned (RECWIN = 3).  

         minbck = -29
         minfwd = 30
         recwin = 3
         CALL MSETWIN(minbck,minfwd,recwin,istatus)
         if(istatus.ne.success_p)stop

C     Load in the stations with data for this hour.

         CALL MRDMTRSTA(atime,nsta_rdmtr,stanam_rdmtr,wmoid_rdmtr,
     1                  lat_rdmtr,lon_rdmtr,elev_rdmtr,timeob_rdmtr,
     2                  provdr_rdmtr,istatus)

         write(1,5)'number of radiometer stations             ',
     1             nsta_rdmtr
         if(istatus.ne.success_p)go to 600

C     Get the floating point grid indices into the CONUS 212 grid.

         do j = 1, nsta_rdmtr
            CALL MGETIJ(lat_rdmtr(j),lon_rdmtr(j),ri_rdmtr(j),
     1                  rj_rdmtr(j),istatus)
         enddo

C     Read the variables in this order:  pressure, height, virtual 
C     temperature, specific humidity.

         CALL MGETRDMTR(atime,'P',nsta_rdmtr,nobs,nlevels_rdmtr,
     1                  p_rdmtr,p_rdmtr_qcd,p_rdmtr_qca,p_rdmtr_qcr,
     2                  istatus)
         if(istatus.ne.success_p)go to 600
         write(1,5)'number of non-missing pressure            ',nobs

         CALL MGETRDMTR(atime,'HT',nsta_rdmtr,nobs,nlevels_rdmtr,
     1                  ht_rdmtr,ht_rdmtr_qcd,ht_rdmtr_qca,ht_rdmtr_qcr,
     2                  istatus)
         if(istatus.ne.success_p)go to 600
         write(1,5)'number of non-missing height              ',nobs

         CALL MGETRDMTR(atime,'TV',nsta_rdmtr,nobs,nlevels_rdmtr,
     1                  tv_rdmtr,tv_rdmtr_qcd,tv_rdmtr_qca,tv_rdmtr_qcr,
     2                  istatus)
         if(istatus.ne.success_p)go to 600
         write(1,5)'number of non-missing virtual temperature ',nobs

         CALL MGETRDMTR(atime,'Q',nsta_rdmtr,nobs,nlevels_rdmtr,
     1                  q_rdmtr,q_rdmtr_qcd,q_rdmtr_qca,q_rdmtr_qcr,
     2                  istatus)
         if(istatus.ne.success_p)go to 600
         write(1,5)'number of non-missing specific humidity   ',nobs

C     Get the satellite wind data for this hour
C     -------------------------------------

 600     write(1,4)'Inventory of satellite wind data for ',atime

         CALL MINIT('SATWND','FSL',.true.,istatus)
         if(istatus.ne.success_p)stop

         CALL MSETDOM(3,40635.25,12.19,-133.459,-95.0,25.0,
     1        25.0,dummy,185,129,istatus)
         if(istatus.ne.success_p)stop

         qctype = 99
         CALL MSETQC(qctype,istatus)
         if(istatus.ne.success_p)stop

C     Set the time window to get satellite wind data within +30 or -120 
C     minutes from the nominal time.  The wind products are available on
C     different schedules for the different satellites, and for the 
C     different types of wind products.  The observation times vary, and 
C     the latency ranges from about one and a half to three hours 
C     (latency is defined as the time the data are available minus the 
C     observation time).  RECWIN will be set to return all records in
C     the window.

         minbck = -120
         minfwd = 30
         recwin = 4
         CALL MSETWIN(minbck,minfwd,recwin,istatus)
         if(istatus.ne.success_p)stop

C     Load in the stations with data for this hour.  Note that P is
C     returned in the station call as the vertical coordinate for 
C     satellite winds.

         CALL MSATWNDSTA(atime,nsta_satw,provdr_satw,lat_satw,lon_satw,
     1                   p_satw,timeob_satw,istatus)

         write(1,5)'number of satellite wind reports          ',
     1             nsta_satw
         if(istatus.ne.success_p)go to 700

C     Get the floating point grid indices into the CONUS 212 grid.

         do j = 1, nsta_satw
            CALL MGETIJ(lat_satw(j),lon_satw(j),ri_satw(j),
     1                  rj_satw(j),istatus)
         enddo

C     Read the variables in this order:  u, v.

         CALL MGETSATWND(atime,'U',nsta_satw,nobs,u_satw,u_satw_qcd,
     1                   u_satw_qca,u_satw_qcr,istatus)
         if(istatus.ne.success_p)go to 700
         write(1,5)'number of non-missing u                   ',nobs

         CALL MGETSATWND(atime,'V',nsta_satw,nobs,v_satw,v_satw_qcd,
     1                   v_satw_qca,v_satw_qcr,istatus)
         if(istatus.ne.success_p)go to 700
         write(1,5)'number of non-missing v                   ',nobs

C     Get the satellite soundings data for this hour
C     -------------------------------------

 700     write(1,4)'Inventory of satellite soundings data for ',atime

         CALL MINIT('SATSND','FSL',.true.,istatus)
         if(istatus.ne.success_p)stop

         CALL MSETDOM(3,40635.25,12.19,-133.459,-95.0,25.0,
     1        25.0,dummy,185,129,istatus)
         if(istatus.ne.success_p)stop

         qctype = 99
         CALL MSETQC(qctype,istatus)
         if(istatus.ne.success_p)stop

C     Set the time window to get satellite soundings data within +30 or
C      -120 minutes from the nominal time.  The wind products are 
C     available on different schedules for the different satellites.  
C     The observation times vary, and the latency ranges from about 
C     one and a half to three hours (latency is defined as the time 
C     the data are available minus the observation time).  RECWIN will 
C     be set to return all records in the window.

         minbck = -120
         minfwd = 30
         recwin = 4
         CALL MSETWIN(minbck,minfwd,recwin,istatus)
         if(istatus.ne.success_p)stop

C     Load in the stations with data for this hour.  Note that P is
C     returned in the station call as the vertical coordinate for 
C     satellite soundings.

         CALL MSATSNDSTA(atime,nsta_sats,provdr_sats,lat_sats,lon_sats,
     1                   p_sats,timeob_sats,istatus)

         write(1,5)'number of satellite soundings reports     ',
     1             nsta_sats
         if(istatus.ne.success_p)go to 800

C     Get the floating point grid indices into the CONUS 212 grid.

         do j = 1, nsta_sats
            CALL MGETIJ(lat_sats(j),lon_sats(j),ri_sats(j),
     1                  rj_sats(j),istatus)
         enddo

C     Read the variables:  temperature, mixing ratio, and heights.

         CALL MGETSATSND(atime,'T',nsta_sats,nobs,
     1        nlevels_sats,t_sats,t_sats_qcd,t_sats_qca,t_sats_qcr,
     2        istatus)
         if(istatus.ne.success_p)go to 800
         write(1,5)'number of non-missing t                   ',nobs

         CALL MGETSATSND(atime,'WVMR',nsta_sats,nobs,
     1        nlevels_sats,mr_sats,mr_sats_qcd,mr_sats_qca,mr_sats_qcr,
     2        istatus)
         if(istatus.ne.success_p)go to 800
         write(1,5)'number of non-missing mr                  ',nobs

         CALL MGETSATSND(atime,'P',nsta_sats,nobs,
     1        nlevels_sats,pr_sats,pr_sats_qcd,pr_sats_qca,pr_sats_qcr,
     2        istatus)
         if(istatus.ne.success_p)go to 800
         write(1,5)'number of non-missing p                  ',nobs



C     Get the WISDOM data for this hour
C     -------------------------------------

 800     write(1,4)'Inventory of WISDOM data for ',atime

         CALL MINIT('WISDOM','FSL',.true.,istatus)
         if(istatus.ne.success_p)stop

         CALL MSETDOM(3,40635.25,12.19,-133.459,-95.0,25.0,
     1        25.0,dummy,185,129,istatus)
         if(istatus.ne.success_p)stop

         qctype = 99
         CALL MSETQC(qctype,istatus)
         if(istatus.ne.success_p)stop

C     Set the time window to get WISDOM data within +29 or -30 
C     minutes from the nominal time.  WISDOM stations are generally
C     reported all around the hour, so you should set your time window 
C     to whatever you'd really like to get.  If a station has multiple 
C     reports within the time window, only the latest one within the
C     window will be returned (RECWIN = 3).  

         minbck = -29
         minfwd = 30
         recwin = 3
         CALL MSETWIN(minbck,minfwd,recwin,istatus)
         if(istatus.ne.success_p)stop

C     Load in the stations with data for this hour.  

         CALL MWISDOMSTA(atime,nsta_wisd,stanam_wisd,provdr_wisd,
     1        lat_wisd,lon_wisd,elev_wisd,timeob_wisd,istatus)

         write(1,5)'number of WISDOM reports     ',
     1             nsta_wisd
         if(istatus.ne.success_p)go to 900

C     Get the floating point grid indices into the CONUS 212 grid.

         do j = 1, nsta_wisd
            CALL MGETIJ(lat_wisd(j),lon_wisd(j),ri_wisd(j),
     1                  rj_wisd(j),istatus)
         enddo

C     Read the variables in this order:  u, v.

         CALL MGETWISDOM(atime,'U',nsta_wisd,nobs,u_wisd,u_wisd_qcd,
     1                   u_wisd_qca,u_wisd_qcr,istatus)
         if(istatus.ne.success_p)go to 700
         write(1,5)'number of non-missing u                   ',nobs

         CALL MGETWISDOM(atime,'V',nsta_wisd,nobs,v_wisd,v_wisd_qcd,
     1                   v_wisd_qca,v_wisd_qcr,istatus)
         if(istatus.ne.success_p)go to 700
         write(1,5)'number of non-missing v                   ',nobs


C     Add one hour to the time string
C     -------------------------------

 900     CALL MCHGTIM(atime,3600,new_time,istatus)
         if(istatus.ne.success_p)stop
         atime = new_time


      enddo

C-----------------------------------------------------------------------
C OPTIONAL MADIS CLOSE ROUTINE
C ----------------------------
C
C     If you are embedding MADIS library calls in a larger program
C     that does additional I/O using netCDF, you may wish to close
C     all the MADIS netCDF files when you're done with the portion
C     of the program that uses MADIS, simply to avoid any possible
C     issues with using too many netCDF files.  If so, just call
C     this subroutine.  Otherwise, there's no need to call the routine.

      CALL MMADISCLOSE

      stop
      end
C**---------------------------------------------------------------------
C**REQUIRED STANDARD FSL DISCLAIMER (NOVEMBER 2000)
C**---------------------------------------------------------------------
C**  OPEN SOURCE LICENSE/DISCLAIMER,  FORECAST SYSTEMS LABORATORY
C**  NOAA/OAR/FSL, 325 BROADWAY BOULDER, CO 80305
C**
C**  THIS SOFTWARE IS DISTRIBUTED UNDER THE OPEN SOURCE DEFINITION, 
C**  WHICH MAY BE FOUND AT http://www.opensource.org/osd.html.
C**
C**  IN PARTICULAR, REDISTRIBUTION AND USE IN SOURCE AND BINARY FORMS,
C**  WITH OR WITHOUT MODIFICATION, ARE PERMITTED PROVIDED THAT THE 
C**  FOLLOWING CONDITIONS ARE MET:
C**
C**  - REDISTRIBUTIONS OF SOURCE CODE MUST RETAIN THIS NOTICE, THIS LIST
C**    OF CONDITIONS AND THE FOLLOWING DISCLAIMER.
C**
C**  - REDISTRIBUTIONS IN BINARY FORM MUST PROVIDE ACCESS TO THIS 
C**    NOTICE, THIS LIST OF CONDITIONS AND THE FOLLOWING DISCLAIMER, AND
C**    THE UNDERLYING SOURCE CODE.
C**
C**  - ALL MODIFICATIONS TO THIS SOFTWARE MUST BE CLEARLY DOCUMENTED, 
C**    AND ARE SOLELY THE RESPONSIBILITY OF THE AGENT MAKING THE 
C**    MODIFICATIONS.
C**
C**  - IF SIGNIFICANT MODIFICATIONS OR ENHANCEMENTS ARE MADE TO THIS
C**    SOFTWARE, THE FSL SOFTWARE POLICY MANAGER 
C**    (softwaremgr.fsl@noaa.gov) BE NOTIFIED.
C**
C**  THIS SOFTWARE AND ITS DOCUMENTATION ARE IN THE PUBLIC DOMAIN AND 
C**  ARE FURNISHED "AS IS."  THE AUTHORS, THE UNITED STATES GOVERNMENT,
C**  ITS INSTRUMENTALITIES, OFFICERS, EMPLOYEES, AND AGENTS MAKE NO 
C**  WARRANTY, EXPRESS OR IMPLIED, AS TO THE USEFULNESS OF THE SOFTWARE
C**  AND DOCUMENTATION FOR ANY PURPOSE.  THEY ASSUME NO RESPONSIBILITY 
C**  (1) FOR THE USE OF THE SOFTWARE AND DOCUMENTATION; OR (2) TO 
C**  PROVIDE TECHNICAL SUPPORT TO USERS.
C**---------------------------------------------------------------------
