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 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 - 200000 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 - 35000 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(200000) character*9 timeob_sfc(200000) character*10 provdr_sfc(200000) integer wmoid_sfc(200000) real lat_sfc(200000) real lon_sfc(200000) real elev_sfc(200000) real ri_sfc(200000) real rj_sfc(200000) C - Variables real p_sfc(200000) real tv_sfc(200000) real q_sfc(200000) real u_sfc(200000) real v_sfc(200000) C - QC information integer p_sfc_qca(200000) integer tv_sfc_qca(200000) integer q_sfc_qca(200000) integer u_sfc_qca(200000) integer v_sfc_qca(200000) integer p_sfc_qcr(200000) integer tv_sfc_qcr(200000) integer q_sfc_qcr(200000) integer u_sfc_qcr(200000) integer v_sfc_qcr(200000) character*1 p_sfc_qcd(200000) character*1 tv_sfc_qcd(200000) character*1 q_sfc_qcd(200000) character*1 u_sfc_qcd(200000) character*1 v_sfc_qcd(200000) 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(60000) character*9 timeob_satw(60000) real lat_satw(60000) real lon_satw(60000) real p_satw(60000) real ri_satw(60000) real rj_satw(60000) C - Variables real u_satw(60000) real v_satw(60000) C - QC information integer u_satw_qca(60000) integer v_satw_qca(60000) integer u_satw_qcr(60000) integer v_satw_qcr(60000) character*1 u_satw_qcd(60000) character*1 v_satw_qcd(60000) 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----------------------------------------------------------------------- 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. CALL MINIT('SFC','FSL',.true.,istatus) if(istatus.ne.success_p)stop CALL MINIT('RAOB','FSL',.true.,istatus) if(istatus.ne.success_p)stop CALL MINIT('NPN','FSL',.true.,istatus) if(istatus.ne.success_p)stop CALL MINIT('ACARS','FSL',.true.,istatus) if(istatus.ne.success_p)stop CALL MINIT('MAP','FSL',.true.,istatus) if(istatus.ne.success_p)stop CALL MINIT('RDMTR','FSL',.true.,istatus) if(istatus.ne.success_p)stop CALL MINIT('SATWND','FSL',.true.,istatus) if(istatus.ne.success_p)stop CALL MINIT('SATSND','FSL',.true.,istatus) if(istatus.ne.success_p)stop C----------------------------------------------------------------------- C OPTIONAL MADIS INITIALIZATION C ---------------------------- C C Set the horizontal domain for the CONUS 212 Lambert conformal C conic projection. CALL MSETDOM(3,40635.25,12.19,-133.459,-95.0,25.0,25.0,dummy, 1 185,129,istatus) if(istatus.ne.success_p)stop C Set QC level so that data are only returned that have passed C all possible QC checks that can be applied to each individual C variable. qctype = 99 CALL MSETQC(qctype,istatus) 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,/) 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) 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 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 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 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 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 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 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 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_satw 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 Add one hour to the time string C ------------------------------- 800 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**---------------------------------------------------------------------