[Dart-dev] [4282] DART/trunk: Updates to the MADIS converters: whether or not the QC values

nancy at ucar.edu nancy at ucar.edu
Tue Feb 16 13:46:34 MST 2010


Revision: 4282
Author:   nancy
Date:     2010-02-16 13:46:34 -0700 (Tue, 16 Feb 2010)
Log Message:
-----------
Updates to the MADIS converters:  whether or not the QC values
are in the files is detected at runtime.  There is a internal
variable which can be set at compile time to force the programs
to ignore QC fields even if present and use all observations, 
but by default the QC will be used if it exists.  

The single surface converter has been split into two versions;
one which reads the standard METAR files and one which reads
the very similar but slightly different Mesonet files.

Add the additional METAR_ALTIMETER and METAR_ moisture
types.  For now, THERE ARE NO MESONET_xxx types -- i still haven't
gotten a clear consensus on whether this was desired or not.  
For now, this version of the MESONET converters will use LAND_SFC_ 
as the output type for mesonet files.  

Modified Paths:
--------------
    DART/trunk/obs_def/obs_def_altimeter_mod.f90
    DART/trunk/obs_def/obs_def_metar_mod.f90
    DART/trunk/observations/MADIS/convert_madis_acars.f90
    DART/trunk/observations/MADIS/convert_madis_marine.f90
    DART/trunk/observations/MADIS/convert_madis_rawin.f90
    DART/trunk/observations/MADIS/data/README
    DART/trunk/observations/MADIS/data/acars.in
    DART/trunk/observations/MADIS/work/input.nml

Added Paths:
-----------
    DART/trunk/observations/MADIS/convert_madis_mesonet.f90
    DART/trunk/observations/MADIS/convert_madis_metar.f90
    DART/trunk/observations/MADIS/work/mkmf_convert_madis_mesonet
    DART/trunk/observations/MADIS/work/mkmf_convert_madis_metar
    DART/trunk/observations/MADIS/work/path_names_convert_madis_mesonet
    DART/trunk/observations/MADIS/work/path_names_convert_madis_metar

Removed Paths:
-------------
    DART/trunk/observations/MADIS/convert_madis_surface.f90
    DART/trunk/observations/MADIS/work/mkmf_convert_madis_surface
    DART/trunk/observations/MADIS/work/path_names_convert_madis_surface

-------------- next part --------------
Modified: DART/trunk/obs_def/obs_def_altimeter_mod.f90
===================================================================
--- DART/trunk/obs_def/obs_def_altimeter_mod.f90	2010-02-16 20:36:48 UTC (rev 4281)
+++ DART/trunk/obs_def/obs_def_altimeter_mod.f90	2010-02-16 20:46:34 UTC (rev 4282)
@@ -6,6 +6,7 @@
 ! RADIOSONDE_SURFACE_ALTIMETER, KIND_SURFACE_PRESSURE
 ! MARINE_SFC_ALTIMETER,         KIND_SURFACE_PRESSURE
 ! LAND_SFC_ALTIMETER,           KIND_SURFACE_PRESSURE
+! METAR_ALTIMETER,              KIND_SURFACE_PRESSURE
 ! END DART PREPROCESS KIND LIST
 
 ! BEGIN DART PREPROCESS USE OF SPECIAL OBS_DEF MODULE
@@ -13,22 +14,26 @@
 ! END DART PREPROCESS USE OF SPECIAL OBS_DEF MODULE
 
 ! BEGIN DART PREPROCESS GET_EXPECTED_OBS_FROM_DEF
-!         case(RADIOSONDE_SURFACE_ALTIMETER, MARINE_SFC_ALTIMETER, LAND_SFC_ALTIMETER)
+!         case(RADIOSONDE_SURFACE_ALTIMETER, MARINE_SFC_ALTIMETER, LAND_SFC_ALTIMETER, &
+               METAR_ALTIMETER)
 !            call get_expected_altimeter(state, location, obs_val, istatus)
 ! END DART PREPROCESS GET_EXPECTED_OBS_FROM_DEF
 
 ! BEGIN DART PREPROCESS READ_OBS_DEF
-!         case(RADIOSONDE_SURFACE_ALTIMETER, MARINE_SFC_ALTIMETER, LAND_SFC_ALTIMETER)
+!         case(RADIOSONDE_SURFACE_ALTIMETER, MARINE_SFC_ALTIMETER, LAND_SFC_ALTIMETER, &
+               METAR_ALTIMETER)
 !            continue
 ! END DART PREPROCESS READ_OBS_DEF
 
 ! BEGIN DART PREPROCESS WRITE_OBS_DEF
-!         case(RADIOSONDE_SURFACE_ALTIMETER, MARINE_SFC_ALTIMETER, LAND_SFC_ALTIMETER)
+!         case(RADIOSONDE_SURFACE_ALTIMETER, MARINE_SFC_ALTIMETER, LAND_SFC_ALTIMETER, &
+               METAR_ALTIMETER)
 !            continue
 ! END DART PREPROCESS WRITE_OBS_DEF
 
 ! BEGIN DART PREPROCESS INTERACTIVE_OBS_DEF
-!         case(RADIOSONDE_SURFACE_ALTIMETER, MARINE_SFC_ALTIMETER, LAND_SFC_ALTIMETER)
+!         case(RADIOSONDE_SURFACE_ALTIMETER, MARINE_SFC_ALTIMETER, LAND_SFC_ALTIMETER, &
+               METAR_ALTIMETER)
 !            continue
 ! END DART PREPROCESS INTERACTIVE_OBS_DEF
 

Modified: DART/trunk/obs_def/obs_def_metar_mod.f90
===================================================================
--- DART/trunk/obs_def/obs_def_metar_mod.f90	2010-02-16 20:36:48 UTC (rev 4281)
+++ DART/trunk/obs_def/obs_def_metar_mod.f90	2010-02-16 20:46:34 UTC (rev 4282)
@@ -7,6 +7,8 @@
 ! METAR_V_10_METER_WIND,           KIND_V_WIND_COMPONENT,       COMMON_CODE
 ! METAR_TEMPERATURE_2_METER,       KIND_TEMPERATURE,            COMMON_CODE
 ! METAR_SPECIFIC_HUMIDITY_2_METER, KIND_SPECIFIC_HUMIDITY,      COMMON_CODE
+! METAR_RELATIVE_HUMIDITY_2_METER, KIND_RELATIVE_HUMIDITY,      COMMON_CODE
+! METAR_DEWPOINT_2_METER,          KIND_DEWPOINT,               COMMON_CODE
 ! METAR_SURFACE_PRESSURE,          KIND_SURFACE_PRESSURE,       COMMON_CODE
 ! METAR_POT_TEMP_2_METER,          KIND_POTENTIAL_TEMPERATURE,  COMMON_CODE
 ! END DART PREPROCESS KIND LIST

Modified: DART/trunk/observations/MADIS/convert_madis_acars.f90
===================================================================
--- DART/trunk/observations/MADIS/convert_madis_acars.f90	2010-02-16 20:36:48 UTC (rev 4281)
+++ DART/trunk/observations/MADIS/convert_madis_acars.f90	2010-02-16 20:46:34 UTC (rev 4282)
@@ -48,13 +48,13 @@
 character(len=129), parameter :: acars_out_file    = 'obs_seq.acars'
 
 ! the following logical parameters control which water-vapor variables appear in the output file, 
-! whether to use the NCEP error or Lin and Hubbard (2004) moisture error model, and 
-! whether the input data file has QC values.
+! whether to use the NCEP error or Lin and Hubbard (2004) moisture error model, and if the
+! input file has data quality control fields, whether to use or ignore them.
 logical, parameter :: LH_err                    = .false.
 logical, parameter :: include_specific_humidity = .true.
 logical, parameter :: include_relative_humidity = .false.
 logical, parameter :: include_dewpoint          = .false.
-logical, parameter :: input_has_qc              = .true.
+logical, parameter :: use_input_qc              = .true. 
 
 integer, parameter ::   num_copies = 1,   &   ! number of copies in sequence
                         num_qc     = 1        ! number of QC entries
@@ -63,8 +63,9 @@
 character (len=80)  :: name
 character (len=19)  :: datime
 integer :: rcode, ncid, varid, nobs, nvars, n, i, window_sec, dday, dsec, &
-           oday, osec, nused, iyear, imonth, iday, ihour, imin, isec
-logical :: file_exist
+           oday, osec, nused, iyear, imonth, iday, ihour, imin, isec, nfrc
+logical :: file_exist, input_has_qc
+
 real(r8) :: palt_miss, tair_miss, relh_miss, tdew_miss, wdir_miss, wspd_miss, uwnd, &
             vwnd, qobs, qsat, oerr, window_hours, pres, qc, qerr
 
@@ -77,7 +78,7 @@
 type(obs_type)          :: obs
 type(time_type)         :: comp_day0, time_obs, time_anal
 
-print*,'Enter target assimilation time (yyyy-mm-dd_hh:mm:ss), obs window (hours)'
+print*,'Enter target assimilation time (yyyy-mm-dd_hh:mm:ss), obs window (hours): '
 read*,datime,window_hours
 
 call set_calendar_type(GREGORIAN)
@@ -156,8 +157,13 @@
 call check( nf90_inq_varid(ncid, "timeObs", varid) )
 call check( nf90_get_var(ncid, varid, tobs) )
 
+! pick a random QC field and test for it.  if it's there, set
+! the 'has qc' flag to true.  otherwise, set it to false.
+nfrc = nf90_inq_varid(ncid, "altitudeQCR", varid) 
+input_has_qc = (nfrc == NF90_NOERR)
+
 ! read the QC check for each variable
-if (input_has_qc) then
+if (input_has_qc .and. use_input_qc) then
    call check( nf90_inq_varid(ncid, "altitudeQCR", varid) )
    call check( nf90_get_var(ncid, varid, qc_palt) )
    
@@ -176,7 +182,7 @@
    call check( nf90_inq_varid(ncid, "windSpeedQCR", varid) )
    call check( nf90_get_var(ncid, varid, qc_wspd) )
 else
-   ! if input contains no QCs, assume all are ok.
+   ! if input contains no QCs, or user said skip them. assume all are ok.
    qc_palt = 0;  qc_relh = 0
    qc_tair = 0;  qc_tdew = 0
    qc_wdir = 0;  qc_wspd = 0
@@ -322,9 +328,11 @@
 
 if ( get_num_obs(obs_seq) > 0 )  call write_obs_seq(obs_seq, acars_out_file)
 
-stop
-end
+! end of main program
 
+contains
+
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !
 !   check - subroutine that checks the flag from a netCDF function.  If
@@ -398,3 +406,5 @@
 
 return
 end subroutine create_obs_type
+
+end program convert_madis_acars

Modified: DART/trunk/observations/MADIS/convert_madis_marine.f90
===================================================================
--- DART/trunk/observations/MADIS/convert_madis_marine.f90	2010-02-16 20:36:48 UTC (rev 4281)
+++ DART/trunk/observations/MADIS/convert_madis_marine.f90	2010-02-16 20:46:34 UTC (rev 4282)
@@ -56,13 +56,13 @@
 character(len=129), parameter :: marine_out_file    = 'obs_seq.marine_sfc'
 
 ! the following logical parameters control which water-vapor variables appear in the output file,
-! whether to use the NCEP error or Lin and Hubbard (2004) moisture error model, and
-! whether the input data file has QC values.
+! whether to use the NCEP error or Lin and Hubbard (2004) moisture error model, and if the
+! input file has data quality control fields, whether to use or ignore them.
 logical, parameter :: LH_err                    = .false.
 logical, parameter :: include_specific_humidity = .true.
 logical, parameter :: include_relative_humidity = .false.
 logical, parameter :: include_dewpoint          = .false.
-logical, parameter :: input_has_qc              = .true.
+logical, parameter :: use_input_qc              = .true. 
 
 integer, parameter   :: dsecobs = 2700,   &   ! observation window
                         num_copies = 1,   &   ! number of copies in sequence
@@ -74,8 +74,8 @@
 character (len=80)  :: name
 character (len=19)  :: datestr
 integer :: rcode, ncid, varid, nobs, nvars, n, i, dday, dsec, oday, &
-           osec, nused, iyear, imonth, iday, ihour, imin, isec
-logical :: file_exist
+           osec, nused, iyear, imonth, iday, ihour, imin, isec, nfrc
+logical :: file_exist, input_has_qc
 real(r8) :: sfcp_miss, tair_miss, tdew_miss, wdir_miss, wspd_miss, uwnd, &
             vwnd, altim, palt, oerr, qobs, qerr, qsat, qobserr, rh, slp_miss, elev_miss, qc
 
@@ -88,7 +88,7 @@
 type(obs_type)          :: obs
 type(time_type)         :: comp_day0, time_anal, time_obs
 
-print*,'Enter the observation date (yyyy-mm-dd_hh:mm:ss)'
+print*,'Enter the observation date (yyyy-mm-dd_hh:mm:ss): '
 read*,datestr
 
 call set_calendar_type(GREGORIAN)
@@ -175,8 +175,14 @@
 call check( nf90_inq_varid(ncid, "timeObs", varid) )
 call check( nf90_get_var(ncid, varid, tobs) )
 
+! pick a random QC field and test for it.  if it's there, set
+! the 'has qc' flag to true.  otherwise, set it to false.
 ! read the QC check for each variable
-if (input_has_qc) then
+nfrc = nf90_inq_varid(ncid, "stationPressQCR", varid) 
+input_has_qc = (nfrc == nf90_noerr)
+
+! read the QC check for each variable
+if (input_has_qc .and. use_input_qc) then
    call check( nf90_inq_varid(ncid, "stationPressQCR", varid) )
    call check( nf90_get_var(ncid, varid, qc_sfcp) )
    
@@ -195,7 +201,7 @@
    call check( nf90_inq_varid(ncid, "windSpeedQCR", varid) )
    call check( nf90_get_var(ncid, varid, qc_wspd) )
 else
-   ! if input contains no QCs, assume all are ok.
+   ! if input contains no QCs, or user said skip them. assume all are ok.
    qc_sfcp = 0;  qc_slp  = 0
    qc_tair = 0;  qc_tdew = 0
    qc_wdir = 0;  qc_wspd = 0
@@ -399,10 +405,11 @@
 
 if ( get_num_obs(obs_seq) > 0 )  call write_obs_seq(obs_seq, marine_out_file)
 
-stop
-end
+! end of main program
 
+contains
 
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !
 !   check - subroutine that checks the flag from a netCDF function.  If
@@ -480,3 +487,6 @@
 
 return
 end subroutine create_obs_type
+
+
+end program convert_madis_marine

Added: DART/trunk/observations/MADIS/convert_madis_mesonet.f90
===================================================================
--- DART/trunk/observations/MADIS/convert_madis_mesonet.f90	                        (rev 0)
+++ DART/trunk/observations/MADIS/convert_madis_mesonet.f90	2010-02-16 20:46:34 UTC (rev 4282)
@@ -0,0 +1,432 @@
+! DART software - Copyright \xA9 2004 - 2010 UCAR. This open source software is
+! provided by UCAR, "as is", without charge, subject to all terms of use at
+! http://www.image.ucar.edu/DAReS/DART/DART_download
+
+program convert_madis_mesonet
+
+! <next few lines under version control, do not edit>
+! $URL$
+! $Id$
+! $Revision$
+! $Date$
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+!   convert_madis_mesonet - program that reads a MADIS netCDF land 
+!                           surface observation file and writes a DART
+!                           obs_seq file using the DART library routines.
+!                           This version works on generic MESONET files.
+!
+!     created Dec. 2007 Ryan Torn, NCAR/MMM
+!     modified Dec. 2008 Soyoung Ha and David Dowell, NCAR/MMM
+!     - added dewpoint as an output variable
+!     - added relative humidity as an output variable
+!
+!     modified to include QC_flag check (Soyoung Ha, NCAR/MMM, 08-04-2009)
+!     split from generic surface converter since the METAR and MESONET
+!     files differ slightly in format.  (Glen Romine, NCAR/MMM  Feb 2010)
+!
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+use        types_mod, only : r8, missing_r8
+use time_manager_mod, only : time_type, set_calendar_type, set_date, &
+                             increment_time, get_time, GREGORIAN, operator(-)
+use     location_mod, only : VERTISSURFACE
+use obs_sequence_mod, only : obs_sequence_type, obs_type, read_obs_seq, &
+                             static_init_obs_sequence, init_obs, write_obs_seq, & 
+                             append_obs_to_seq, init_obs_sequence, get_num_obs, & 
+                             set_copy_meta_data, set_qc_meta_data
+use       meteor_mod, only : sat_vapor_pressure, specific_humidity, & 
+                             wind_dirspd_to_uv, invert_altimeter, pres_alt_to_pres, &
+                             temp_and_dewpoint_to_rh
+use      obs_err_mod, only : land_temp_error, land_wind_error, &
+                             land_pres_error, land_rel_hum_error
+use dewpoint_obs_err_mod, only : dewpt_error_from_rh_and_temp, &
+                                 rh_error_from_dewpt_and_temp
+use     obs_kind_mod, only : LAND_SFC_U_WIND_COMPONENT, LAND_SFC_V_WIND_COMPONENT, &
+                             LAND_SFC_TEMPERATURE, LAND_SFC_SPECIFIC_HUMIDITY, & 
+                             LAND_SFC_DEWPOINT, LAND_SFC_RELATIVE_HUMIDITY, &
+                             LAND_SFC_ALTIMETER                
+
+use           netcdf
+
+implicit none
+
+character(len=16),  parameter :: surface_netcdf_file = 'mesonet_input.nc'
+character(len=129), parameter :: surface_out_file    = 'obs_seq.mesonet'
+
+! the following logical parameters control which water-vapor variables appear in the output file,
+! whether to use the NCEP error or Lin and Hubbard (2004) moisture error model, and if the
+! input file has data quality control fields, whether to use or ignore them.
+logical, parameter :: LH_err                    = .false.
+logical, parameter :: include_specific_humidity = .true.
+logical, parameter :: include_relative_humidity = .false.
+logical, parameter :: include_dewpoint          = .false.
+logical, parameter :: use_input_qc              = .true. 
+
+integer, parameter   :: dsecobs    = 420, &   ! observation window
+                        num_copies = 1,   &   ! number of copies in sequence
+                        num_qc     = 1        ! number of QC entries
+
+character (len=129) :: meta_data
+character (len=80)  :: name
+character (len=19)  :: datestr
+character (len=5)   :: rtype
+integer :: rcode, ncid, varid, nobs, nvars, n, i, oday, osec, dday, &
+           dsec, nused, iyear, imonth, iday, ihour, imin, isec, nfrc
+logical :: file_exist, input_has_qc
+real(r8) :: alti_miss, tair_miss, tdew_miss, wdir_miss, wspd_miss, uwnd, &
+            vwnd, palt, qobs, qsat, rh, oerr, pres, qerr, qc
+
+integer,  allocatable :: tobs(:)
+real(r8), allocatable :: lat(:), lon(:), elev(:), alti(:), tair(:), & 
+                         tdew(:), wdir(:), wspd(:), latu(:), lonu(:)
+integer,  allocatable :: qc_alti(:), qc_tair(:), qc_tdew(:), qc_wdir(:), qc_wspd(:)
+
+type(obs_sequence_type) :: obs_seq
+type(obs_type)          :: obs
+type(time_type)         :: comp_day0, time_obs, time_anal
+
+print*,'Enter the analysis time (yyyy-mm-dd_hh:mm:ss):'
+read*,datestr
+
+! put the analysis date into DART format
+call set_calendar_type(GREGORIAN)
+read(datestr(1:4),   fmt='(i4)') iyear
+read(datestr(6:7),   fmt='(i2)') imonth
+read(datestr(9:10),  fmt='(i2)') iday
+read(datestr(12:13), fmt='(i2)') ihour
+read(datestr(15:16), fmt='(i2)') imin
+read(datestr(18:19), fmt='(i2)') isec
+time_anal = set_date(iyear, imonth, iday, ihour, imin, isec)
+call get_time(time_anal, osec, oday)
+comp_day0 = set_date(1970, 1, 1, 0, 0, 0)
+
+rcode = nf90_open(surface_netcdf_file, nf90_nowrite, ncid)
+call check( nf90_inq_dimid(ncid, "recNum", varid) )
+call check( nf90_inquire_dimension(ncid, varid, name, nobs) )
+
+allocate( lat(nobs))  ;  allocate( lon(nobs))
+allocate(latu(nobs))  ;  allocate(lonu(nobs))
+allocate(elev(nobs))  ;  allocate(alti(nobs))
+allocate(tair(nobs))  ;  allocate(tdew(nobs))
+allocate(wdir(nobs))  ;  allocate(wspd(nobs))
+allocate(tobs(nobs))
+
+nvars = 4
+if (include_specific_humidity) nvars = nvars + 1
+if (include_relative_humidity) nvars = nvars + 1
+if (include_dewpoint) nvars = nvars + 1
+
+allocate(qc_alti(nobs))
+allocate(qc_tair(nobs)) ;  allocate(qc_tdew(nobs))
+allocate(qc_wdir(nobs)) ;  allocate(qc_wspd(nobs))
+
+! read the latitude array
+call check( nf90_inq_varid(ncid, "latitude", varid) )
+call check( nf90_get_var(ncid, varid, lat) )
+
+! read the latitude array
+call check( nf90_inq_varid(ncid, "longitude", varid) )
+call check( nf90_get_var(ncid, varid, lon) )
+
+! read the elevation array
+call check( nf90_inq_varid(ncid, "elevation", varid) )
+call check( nf90_get_var(ncid, varid, elev) )
+
+! read the altimeter setting array
+call check( nf90_inq_varid(ncid, "altimeter", varid) )
+call check( nf90_get_var(ncid, varid, alti) )
+call check( nf90_get_att(ncid, varid, '_FillValue', alti_miss) )
+
+! read the air temperature array
+call check( nf90_inq_varid(ncid, "temperature", varid) )
+call check( nf90_get_var(ncid, varid, tair) )
+call check( nf90_get_att(ncid, varid, '_FillValue', tair_miss) )
+
+! read the dew-point temperature array
+call check( nf90_inq_varid(ncid, "dewpoint", varid) )
+call check( nf90_get_var(ncid, varid, tdew) )
+call check( nf90_get_att(ncid, varid, '_FillValue', tdew_miss) )
+
+! read the wind direction array
+call check( nf90_inq_varid(ncid, "windDir", varid) )
+call check( nf90_get_var(ncid, varid, wdir) )
+call check( nf90_get_att(ncid, varid, '_FillValue', wdir_miss) )
+
+! read the wind speed array
+call check( nf90_inq_varid(ncid, "windSpeed", varid) )
+call check( nf90_get_var(ncid, varid, wspd) )
+call check( nf90_get_att(ncid, varid, '_FillValue', wspd_miss) )
+
+! read the observation time array
+call check( nf90_inq_varid(ncid, "observationTime", varid) )
+call check( nf90_get_var(ncid, varid, tobs) )
+
+! pick a random QC field and test for it.  if it's there, set
+! the 'has qc' flag to true.  otherwise, set it to false.
+! read the QC check for each variable
+nfrc = nf90_inq_varid(ncid, "altimeterQCR", varid) 
+input_has_qc = (nfrc == nf90_noerr)
+
+! read the QC check for each variable
+if (input_has_qc .and. use_input_qc) then
+   call check( nf90_inq_varid(ncid, "altimeterQCR", varid) )
+   call check( nf90_get_var(ncid, varid, qc_alti) )
+   
+   call check( nf90_inq_varid(ncid, "temperatureQCR", varid) )
+   call check( nf90_get_var(ncid, varid, qc_tair) )
+   
+   call check( nf90_inq_varid(ncid, "dewpointQCR", varid) )
+   call check( nf90_get_var(ncid, varid, qc_tdew) )
+   
+   call check( nf90_inq_varid(ncid, "windDirQCR", varid) )
+   call check( nf90_get_var(ncid, varid, qc_wdir) )
+   
+   call check( nf90_inq_varid(ncid, "windSpeedQCR", varid) )
+   call check( nf90_get_var(ncid, varid, qc_wspd) )
+else
+   ! if input contains no QCs, or user said skip them. assume all are ok.
+   qc_alti = 0
+   qc_tair = 0 ;  qc_tdew = 0
+   qc_wdir = 0 ;  qc_wspd = 0
+endif
+
+!  either read existing obs_seq or create a new one
+call static_init_obs_sequence()
+call init_obs(obs, num_copies, num_qc)
+inquire(file=surface_out_file, exist=file_exist)
+if ( file_exist ) then
+
+  call read_obs_seq(surface_out_file, 0, 0, nvars*nobs, obs_seq)
+
+else
+
+  call init_obs_sequence(obs_seq, num_copies, num_qc, nvars*nobs)
+  do i = 1, num_copies
+    meta_data = 'MADIS observation'
+    call set_copy_meta_data(obs_seq, i, meta_data)
+  end do
+  do i = 1, num_qc
+    meta_data = 'Data QC'
+    call set_qc_meta_data(obs_seq, i, meta_data)
+  end do
+
+end if
+
+nused = 0
+obsloop: do n = 1, nobs
+
+  ! determine whether observation is close to analysis time
+  time_obs = increment_time(comp_day0, mod(tobs(n),86400), tobs(n) / 86400)
+  call get_time((time_anal - time_obs), dsec, dday)
+  if ( (dsec + dday * 86400) > dsecobs ) cycle obsloop
+  if ( lon(n) < 0.0_r8 )  lon(n) = lon(n) + 360.0_r8
+
+  do i = 1, nused
+    if ( lon(n) == lonu(i) .and. lat(n) == latu(i) ) cycle obsloop
+  end do
+  qc = 1.0_r8
+  palt = pres_alt_to_pres(elev(n)) * 0.01_r8
+
+  ! add altimeter data to text file
+  if ( alti(n) /= alti_miss .and. qc_alti(n) == 0 ) then
+
+    pres = invert_altimeter(alti(n) * 0.01_r8, elev(n))
+    oerr = land_pres_error(palt)
+    if ( alti(n) >= 89000.0_r8 .and. alti(n) <= 110000.0_r8 .and. oerr /= missing_r8 ) then
+
+      call create_obs_type(lat(n), lon(n), elev(n), VERTISSURFACE, alti(n) * 0.01_r8, & 
+                           LAND_SFC_ALTIMETER, oerr, oday, osec, qc, obs)
+      call append_obs_to_seq(obs_seq, obs)
+
+    end if
+
+  end if
+
+  ! add wind component data to text file
+  if ( wdir(n) /= wdir_miss .and. wspd(n) /= wspd_miss .and. qc_wdir(n) == 0 .and. qc_wspd(n) == 0 ) then
+
+    call wind_dirspd_to_uv(wdir(n), wspd(n), uwnd, vwnd)
+    oerr = land_wind_error(palt)
+    if ( abs(uwnd) < 150.0_r8 .and. abs(vwnd) < 150.0_r8 .and. oerr /= missing_r8 ) then
+
+      call create_obs_type(lat(n), lon(n), elev(n), VERTISSURFACE, uwnd, &
+                           LAND_SFC_U_WIND_COMPONENT, oerr, oday, osec, qc, obs)
+      call append_obs_to_seq(obs_seq, obs)
+      call create_obs_type(lat(n), lon(n), elev(n), VERTISSURFACE, vwnd, &
+                           LAND_SFC_V_WIND_COMPONENT, oerr, oday, osec, qc, obs)
+      call append_obs_to_seq(obs_seq, obs)
+
+    end if
+
+  end if
+
+  ! add air temperature data to text file
+  if ( tair(n) /= tair_miss .and. qc_tair(n) == 0 ) then 
+
+    oerr = land_temp_error(palt)
+    if ( tair(n) >= 200.0_r8 .and. tair(n) <= 335.0_r8 .and. oerr /= missing_r8 ) then
+
+      call create_obs_type(lat(n), lon(n), elev(n), VERTISSURFACE, tair(n), &
+                           LAND_SFC_TEMPERATURE, oerr, oday, osec, qc, obs)
+      call append_obs_to_seq(obs_seq, obs)
+
+    end if
+
+  end if
+
+  ! add specific humidity to text file
+  if ( include_specific_humidity .and. tair(n) /= tair_miss .and. tdew(n) /= tdew_miss & 
+       .and. alti(n) /= alti_miss .and. qc_tair(n) == 0 .and. qc_tdew(n) == 0 .and.    &
+       qc_alti(n) == 0 ) then
+
+    qobs = specific_humidity(sat_vapor_pressure(tdew(n)), pres * 100.0_r8)
+    qsat = specific_humidity(sat_vapor_pressure(tair(n)), pres * 100.0_r8)
+    if (LH_err ) then
+      qerr = rh_error_from_dewpt_and_temp(tair(n), tdew(n))
+    else
+      qerr = land_rel_hum_error(pres, tair(n), qobs / qsat)
+    end if
+    oerr = max(qerr * qsat, 0.0001_r8)
+
+    if ( qobs > 0.0_r8 .and. qobs <= 0.07_r8 .and. qerr /= missing_r8 ) then
+
+      call create_obs_type(lat(n), lon(n), elev(n), VERTISSURFACE, qobs, &
+                           LAND_SFC_SPECIFIC_HUMIDITY, oerr, oday, osec, qc, obs)
+      call append_obs_to_seq(obs_seq, obs)
+
+    end if
+
+  end if
+
+  ! add relative humidity data to text file
+  if ( include_relative_humidity .and. tdew(n) /= tdew_miss .and. tair(n) /= tair_miss &
+       .and. qc_tair(n) == 0 .and. qc_tdew(n) == 0 ) then
+
+    rh = temp_and_dewpoint_to_rh(tair(n), tdew(n))
+    if (LH_err ) then
+      oerr = rh_error_from_dewpt_and_temp(tair(n), tdew(n))
+    else
+      oerr = land_rel_hum_error(pres, tair(n), rh)    
+    end if
+
+    if ( rh > 0.0_r8 .and. rh <= 1.5_r8 .and. oerr /= missing_r8 ) then
+
+    call create_obs_type(lat(n), lon(n), elev(n), VERTISSURFACE, rh, &
+                         LAND_SFC_RELATIVE_HUMIDITY, oerr, oday, osec, qc, obs)
+    call append_obs_to_seq(obs_seq, obs)
+
+  end if
+
+  end if
+
+  ! add dew-point temperature data to text file
+  if ( include_dewpoint .and. tdew(n) /= tdew_miss .and. tair(n) /= tair_miss  &
+       .and. qc_tair(n) == 0 .and. qc_tdew(n) == 0 ) then
+
+    rh = temp_and_dewpoint_to_rh(tair(n), tdew(n))
+    oerr = dewpt_error_from_rh_and_temp(tair(n), rh)
+
+    if ( rh > 0.0_r8 .and. rh <= 1.5_r8 .and. oerr /= missing_r8 ) then
+
+    call create_obs_type(lat(n), lon(n), elev(n), VERTISSURFACE, tdew(n), &
+                         LAND_SFC_DEWPOINT, oerr, oday, osec, qc, obs)
+    call append_obs_to_seq(obs_seq, obs)
+
+    end if
+
+!    print*, 'temp (C), rh (%), oerr:  ', tair(n)-273.15_r8, rh*100.0_r8, oerr
+
+  end if
+
+  nused = nused + 1
+  latu(nused) = lat(n)
+  lonu(nused) = lon(n)
+
+end do obsloop
+
+if ( get_num_obs(obs_seq) > 0 )  call write_obs_seq(obs_seq, surface_out_file)
+call check( nf90_close(ncid) )
+
+! end of main program
+
+contains
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+!   check - subroutine that checks the flag from a netCDF function.  If 
+!           there is an error, the message is displayed.
+!
+!    istatus - netCDF output flag
+!
+!     created Dec. 2007 Ryan Torn, NCAR/MMM
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine check( istatus ) 
+
+use netcdf
+
+implicit none
+
+integer, intent (in) :: istatus
+
+if(istatus /= nf90_noerr) print*,'Netcdf error: ',trim(nf90_strerror(istatus))
+
+end subroutine check
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+!   create_obs_type - subroutine that is used to create an observation
+!                     type from observation data.
+!
+!    lat   - latitude of observation
+!    lon   - longitude of observation
+!    pres  - pressure of observation
+!    vcord - vertical coordinate
+!    obsv  - observation value
+!    okind - observation kind
+!    oerr  - observation error
+!    day   - gregorian day
+!    sec   - gregorian second
+!    qc    - quality control value
+!    obs   - observation type
+!
+!     created Oct. 2007 Ryan Torn, NCAR/MMM
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+subroutine create_obs_type(lat, lon, pres, vcord, obsv, okind, oerr, day, sec, qc, obs)
+
+use types_mod,        only : r8
+use obs_sequence_mod, only : obs_type, set_obs_values, set_qc, set_obs_def
+use obs_def_mod,      only : obs_def_type, set_obs_def_time, set_obs_def_kind, &
+                             set_obs_def_error_variance, set_obs_def_location
+use     location_mod, only : location_type, set_location
+use time_manager_mod, only : time_type, set_time
+
+implicit none
+
+integer, intent(in)         :: okind, vcord, day, sec
+real(r8), intent(in)        :: lat, lon, pres, obsv, oerr, qc
+type(obs_type), intent(inout) :: obs
+
+real(r8)              :: obs_val(1), qc_val(1)
+type(obs_def_type)    :: obs_def
+
+call set_obs_def_location(obs_def, set_location(lon, lat, pres, vcord))
+call set_obs_def_kind(obs_def, okind)
+call set_obs_def_time(obs_def, set_time(sec, day))
+call set_obs_def_error_variance(obs_def, oerr * oerr)
+call set_obs_def(obs, obs_def)
+
+obs_val(1) = obsv
+call set_obs_values(obs, obs_val)
+qc_val(1)  = qc
+call set_qc(obs, qc_val)
+
+return
+end subroutine create_obs_type
+
+
+end program convert_madis_mesonet


Property changes on: DART/trunk/observations/MADIS/convert_madis_mesonet.f90
___________________________________________________________________
Added: svn:keywords
   + Date Revision Author HeadURL Id

Added: DART/trunk/observations/MADIS/convert_madis_metar.f90
===================================================================
--- DART/trunk/observations/MADIS/convert_madis_metar.f90	                        (rev 0)
+++ DART/trunk/observations/MADIS/convert_madis_metar.f90	2010-02-16 20:46:34 UTC (rev 4282)
@@ -0,0 +1,439 @@
+! DART software - Copyright \xA9 2004 - 2010 UCAR. This open source software is
+! provided by UCAR, "as is", without charge, subject to all terms of use at
+! http://www.image.ucar.edu/DAReS/DART/DART_download
+
+program convert_madis_metar
+
+! <next few lines under version control, do not edit>
+! $URL$
+! $Id$
+! $Revision$
+! $Date$
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+!   convert_madis_metar - program that reads a MADIS netCDF land 
+!                         surface observation file and writes a DART
+!                         obs_seq file using the DART library routines.
+!                         This version works on the standard METAR files.
+!
+!     created Dec. 2007 Ryan Torn, NCAR/MMM
+!     modified Dec. 2008 Soyoung Ha and David Dowell, NCAR/MMM
+!     - added dewpoint as an output variable
+!     - added relative humidity as an output variable
+!
+!     modified to include QC_flag check (Soyoung Ha, NCAR/MMM, 08-04-2009)
+!     split from the mesonet version (Glen Romine, NCAR/MMM, Feb 2010)
+!
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+use        types_mod, only : r8, missing_r8
+use time_manager_mod, only : time_type, set_calendar_type, set_date, &
+                             increment_time, get_time, GREGORIAN, operator(-)
+use     location_mod, only : VERTISSURFACE
+use obs_sequence_mod, only : obs_sequence_type, obs_type, read_obs_seq, &
+                             static_init_obs_sequence, init_obs, write_obs_seq, & 
+                             append_obs_to_seq, init_obs_sequence, get_num_obs, & 
+                             set_copy_meta_data, set_qc_meta_data
+use       meteor_mod, only : sat_vapor_pressure, specific_humidity, & 
+                             wind_dirspd_to_uv, invert_altimeter, pres_alt_to_pres, &
+                             temp_and_dewpoint_to_rh
+use      obs_err_mod, only : land_temp_error, land_wind_error, &
+                             land_pres_error, land_rel_hum_error
+use dewpoint_obs_err_mod, only : dewpt_error_from_rh_and_temp, &
+                                 rh_error_from_dewpt_and_temp
+use     obs_kind_mod, only : METAR_U_10_METER_WIND, METAR_V_10_METER_WIND, &
+                             METAR_TEMPERATURE_2_METER, METAR_SPECIFIC_HUMIDITY_2_METER, & 
+                             METAR_DEWPOINT_2_METER, METAR_RELATIVE_HUMIDITY_2_METER, &
+                             METAR_ALTIMETER                
+
+use           netcdf
+
+implicit none
+
+! COMPILE TIME OPTION:  by default this converter only processes hourly
+! observations.  if exclude_special is set to .false., then special obs
+! will be converted as well as the normal obs.  but this is not on by default.
+character(len=16),  parameter :: surface_netcdf_file = 'metar_input.nc'
+character(len=129), parameter :: surface_out_file    = 'obs_seq.metar'
+logical,            parameter :: exclude_special     = .true.
+
+! the following logical parameters control which water-vapor variables appear in the output file,
+! whether to use the NCEP error or Lin and Hubbard (2004) moisture error model, and if the
+! input file has data quality control fields, whether to use or ignore them.
+logical, parameter :: LH_err                    = .false.
+logical, parameter :: include_specific_humidity = .true.
+logical, parameter :: include_relative_humidity = .false.
+logical, parameter :: include_dewpoint          = .false.
+logical, parameter :: use_input_qc              = .true. 
+
+integer, parameter   :: dsecobs    = 420, &   ! observation window
+                        num_copies = 1,   &   ! number of copies in sequence
+                        num_qc     = 1        ! number of QC entries
+
+character (len=129) :: meta_data
+character (len=80)  :: name
+character (len=19)  :: datestr
+character (len=5)   :: rtype
+integer :: rcode, ncid, varid, nobs, nvars, n, i, oday, osec, dday, &
+           dsec, nused, iyear, imonth, iday, ihour, imin, isec, nfrc
+logical :: file_exist, input_has_qc
+real(r8) :: alti_miss, tair_miss, tdew_miss, wdir_miss, wspd_miss, uwnd, &
+            vwnd, palt, qobs, qsat, rh, oerr, pres, qerr, qc
+
+integer,  allocatable :: tobs(:)
+real(r8), allocatable :: lat(:), lon(:), elev(:), alti(:), tair(:), & 
+                         tdew(:), wdir(:), wspd(:), latu(:), lonu(:)
+integer,  allocatable :: qc_alti(:), qc_tair(:), qc_tdew(:), qc_wdir(:), qc_wspd(:)
+
+type(obs_sequence_type) :: obs_seq
+type(obs_type)          :: obs
+type(time_type)         :: comp_day0, time_obs, time_anal
+
+print*,'Enter the analysis time (yyyy-mm-dd_hh:mm:ss):'
+read*,datestr
+
+! put the analysis date into DART format
+call set_calendar_type(GREGORIAN)
+read(datestr(1:4),   fmt='(i4)') iyear
+read(datestr(6:7),   fmt='(i2)') imonth
+read(datestr(9:10),  fmt='(i2)') iday
+read(datestr(12:13), fmt='(i2)') ihour
+read(datestr(15:16), fmt='(i2)') imin
+read(datestr(18:19), fmt='(i2)') isec
+time_anal = set_date(iyear, imonth, iday, ihour, imin, isec)
+call get_time(time_anal, osec, oday)
+comp_day0 = set_date(1970, 1, 1, 0, 0, 0)
+
+rcode = nf90_open(surface_netcdf_file, nf90_nowrite, ncid)
+call check( nf90_inq_dimid(ncid, "recNum", varid) )
+call check( nf90_inquire_dimension(ncid, varid, name, nobs) )
+
+allocate( lat(nobs))  ;  allocate( lon(nobs))
+allocate(latu(nobs))  ;  allocate(lonu(nobs))
+allocate(elev(nobs))  ;  allocate(alti(nobs))
+allocate(tair(nobs))  ;  allocate(tdew(nobs))
+allocate(wdir(nobs))  ;  allocate(wspd(nobs))
+allocate(tobs(nobs))
+
+nvars = 4
+if (include_specific_humidity) nvars = nvars + 1
+if (include_relative_humidity) nvars = nvars + 1
+if (include_dewpoint) nvars = nvars + 1
+
+allocate(qc_alti(nobs))
+allocate(qc_tair(nobs)) ;  allocate(qc_tdew(nobs))
+allocate(qc_wdir(nobs)) ;  allocate(qc_wspd(nobs))
+
+! read the latitude array
+call check( nf90_inq_varid(ncid, "latitude", varid) )
+call check( nf90_get_var(ncid, varid, lat) )
+
+! read the latitude array
+call check( nf90_inq_varid(ncid, "longitude", varid) )
+call check( nf90_get_var(ncid, varid, lon) )
+
+! read the elevation array
+call check( nf90_inq_varid(ncid, "elevation", varid) )
+call check( nf90_get_var(ncid, varid, elev) )
+
+! read the altimeter setting array
+call check( nf90_inq_varid(ncid, "altimeter", varid) )
+call check( nf90_get_var(ncid, varid, alti) )
+call check( nf90_get_att(ncid, varid, '_FillValue', alti_miss) )
+
+! read the air temperature array
+call check( nf90_inq_varid(ncid, "temperature", varid) )
+call check( nf90_get_var(ncid, varid, tair) )
+call check( nf90_get_att(ncid, varid, '_FillValue', tair_miss) )
+
+! read the dew-point temperature array
+call check( nf90_inq_varid(ncid, "dewpoint", varid) )
+call check( nf90_get_var(ncid, varid, tdew) )
+call check( nf90_get_att(ncid, varid, '_FillValue', tdew_miss) )
+
+! read the wind direction array
+call check( nf90_inq_varid(ncid, "windDir", varid) )
+call check( nf90_get_var(ncid, varid, wdir) )
+call check( nf90_get_att(ncid, varid, '_FillValue', wdir_miss) )
+
+! read the wind speed array
+call check( nf90_inq_varid(ncid, "windSpeed", varid) )
+call check( nf90_get_var(ncid, varid, wspd) )
+call check( nf90_get_att(ncid, varid, '_FillValue', wspd_miss) )
+
+! read the observation time array
+call check( nf90_inq_varid(ncid, "timeObs", varid) )
+call check( nf90_get_var(ncid, varid, tobs) )
+
+! pick a random QC field and test for it.  if it's there, set
+! the 'has qc' flag to true.  otherwise, set it to false.
+! read the QC check for each variable
+nfrc = nf90_inq_varid(ncid, "altimeterQCR", varid) 
+input_has_qc = (nfrc == nf90_noerr)
+
+! read the QC check for each variable
+if (input_has_qc .and. use_input_qc) then
+   call check( nf90_inq_varid(ncid, "altimeterQCR", varid) )
+   call check( nf90_get_var(ncid, varid, qc_alti) )
+   
+   call check( nf90_inq_varid(ncid, "temperatureQCR", varid) )
+   call check( nf90_get_var(ncid, varid, qc_tair) )
+   
+   call check( nf90_inq_varid(ncid, "dewpointQCR", varid) )
+   call check( nf90_get_var(ncid, varid, qc_tdew) )
+   
+   call check( nf90_inq_varid(ncid, "windDirQCR", varid) )
+   call check( nf90_get_var(ncid, varid, qc_wdir) )
+   
+   call check( nf90_inq_varid(ncid, "windSpeedQCR", varid) )
+   call check( nf90_get_var(ncid, varid, qc_wspd) )
+else
+   ! if input contains no QCs, or user said skip them. assume all are ok.
+   qc_alti = 0
+   qc_tair = 0 ;  qc_tdew = 0
+   qc_wdir = 0 ;  qc_wspd = 0
+endif
+
+!  either read existing obs_seq or create a new one
+call static_init_obs_sequence()
+call init_obs(obs, num_copies, num_qc)
+inquire(file=surface_out_file, exist=file_exist)
+if ( file_exist ) then
+
+  call read_obs_seq(surface_out_file, 0, 0, nvars*nobs, obs_seq)
+
+else
+
+  call init_obs_sequence(obs_seq, num_copies, num_qc, nvars*nobs)
+  do i = 1, num_copies
+    meta_data = 'MADIS observation'
+    call set_copy_meta_data(obs_seq, i, meta_data)
+  end do
+  do i = 1, num_qc
+    meta_data = 'Data QC'
+    call set_qc_meta_data(obs_seq, i, meta_data)
+  end do
+
+end if
+
+nused = 0
+obsloop: do n = 1, nobs
+
+  ! determine whether observation is close to analysis time
+  time_obs = increment_time(comp_day0, mod(tobs(n),86400), tobs(n) / 86400)
+  call get_time((time_anal - time_obs), dsec, dday)
+  if ( (dsec + dday * 86400) > dsecobs ) cycle obsloop
+  if ( lon(n) < 0.0_r8 )  lon(n) = lon(n) + 360.0_r8
+
+  ! check to make sure this observation has not been used
+  call check( nf90_inq_varid(ncid, "reportType", varid) )
+  call check( nf90_get_var(ncid, varid, rtype, start = (/ 1, n /)) )
+  if ( rtype /= 'METAR' .and. exclude_special )  cycle obsloop
+
+  do i = 1, nused
+    if ( lon(n) == lonu(i) .and. lat(n) == latu(i) ) cycle obsloop
+  end do
+  qc = 1.0_r8
+  palt = pres_alt_to_pres(elev(n)) * 0.01_r8
+
+  ! add altimeter data to text file
+  if ( alti(n) /= alti_miss .and. qc_alti(n) == 0 ) then
+
+    pres = invert_altimeter(alti(n) * 0.01_r8, elev(n))
+    oerr = land_pres_error(palt)
+    if ( alti(n) >= 89000.0_r8 .and. alti(n) <= 110000.0_r8 .and. oerr /= missing_r8 ) then
+
+      call create_obs_type(lat(n), lon(n), elev(n), VERTISSURFACE, alti(n) * 0.01_r8, & 
+                           METAR_ALTIMETER, oerr, oday, osec, qc, obs)
+      call append_obs_to_seq(obs_seq, obs)
+
+    end if
+
+  end if
+
+  ! add wind component data to text file
+  if ( wdir(n) /= wdir_miss .and. wspd(n) /= wspd_miss .and. qc_wdir(n) == 0 .and. qc_wspd(n) == 0 ) then
+
+    call wind_dirspd_to_uv(wdir(n), wspd(n), uwnd, vwnd)
+    oerr = land_wind_error(palt)
+    if ( abs(uwnd) < 150.0_r8 .and. abs(vwnd) < 150.0_r8 .and. oerr /= missing_r8 ) then
+
+      call create_obs_type(lat(n), lon(n), elev(n), VERTISSURFACE, uwnd, &
+                           METAR_U_10_METER_WIND, oerr, oday, osec, qc, obs)
+      call append_obs_to_seq(obs_seq, obs)
+      call create_obs_type(lat(n), lon(n), elev(n), VERTISSURFACE, vwnd, &
+                           METAR_V_10_METER_WIND, oerr, oday, osec, qc, obs)
+      call append_obs_to_seq(obs_seq, obs)
+
+    end if
+
+  end if
+
+  ! add air temperature data to text file
+  if ( tair(n) /= tair_miss .and. qc_tair(n) == 0 ) then 
+
+    oerr = land_temp_error(palt)
+    if ( tair(n) >= 200.0_r8 .and. tair(n) <= 335.0_r8 .and. oerr /= missing_r8 ) then
+
+      call create_obs_type(lat(n), lon(n), elev(n), VERTISSURFACE, tair(n), &
+                           METAR_TEMPERATURE_2_METER, oerr, oday, osec, qc, obs)
+      call append_obs_to_seq(obs_seq, obs)
+
+    end if
+
+  end if
+
+  ! add specific humidity to text file
+  if ( include_specific_humidity .and. tair(n) /= tair_miss .and. tdew(n) /= tdew_miss & 
+       .and. alti(n) /= alti_miss .and. qc_tair(n) == 0 .and. qc_tdew(n) == 0 .and.    &
+       qc_alti(n) == 0 ) then
+
+    qobs = specific_humidity(sat_vapor_pressure(tdew(n)), pres * 100.0_r8)
+    qsat = specific_humidity(sat_vapor_pressure(tair(n)), pres * 100.0_r8)
+    if (LH_err ) then
+      qerr = rh_error_from_dewpt_and_temp(tair(n), tdew(n))
+    else
+      qerr = land_rel_hum_error(pres, tair(n), qobs / qsat)
+    end if
+    oerr = max(qerr * qsat, 0.0001_r8)
+
+    if ( qobs > 0.0_r8 .and. qobs <= 0.07_r8 .and. qerr /= missing_r8 ) then
+
+      call create_obs_type(lat(n), lon(n), elev(n), VERTISSURFACE, qobs, &
+                           METAR_SPECIFIC_HUMIDITY_2_METER, oerr, oday, osec, qc, obs)
+      call append_obs_to_seq(obs_seq, obs)
+
+    end if
+
+  end if
+
+  ! add relative humidity data to text file
+  if ( include_relative_humidity .and. tdew(n) /= tdew_miss .and. tair(n) /= tair_miss &
+       .and. qc_tair(n) == 0 .and. qc_tdew(n) == 0 ) then
+
+    rh = temp_and_dewpoint_to_rh(tair(n), tdew(n))
+    if (LH_err ) then
+      oerr = rh_error_from_dewpt_and_temp(tair(n), tdew(n))
+    else
+      oerr = land_rel_hum_error(pres, tair(n), rh)    
+    end if
+
+    if ( rh > 0.0_r8 .and. rh <= 1.5_r8 .and. oerr /= missing_r8 ) then
+
+    call create_obs_type(lat(n), lon(n), elev(n), VERTISSURFACE, rh, &
+                         METAR_RELATIVE_HUMIDITY_2_METER, oerr, oday, osec, qc, obs)
+    call append_obs_to_seq(obs_seq, obs)
+
+  end if
+
+  end if
+
+  ! add dew-point temperature data to text file
+  if ( include_dewpoint .and. tdew(n) /= tdew_miss .and. tair(n) /= tair_miss  &
+       .and. qc_tair(n) == 0 .and. qc_tdew(n) == 0 ) then
+
+    rh = temp_and_dewpoint_to_rh(tair(n), tdew(n))
+    oerr = dewpt_error_from_rh_and_temp(tair(n), rh)
+
+    if ( rh > 0.0_r8 .and. rh <= 1.5_r8 .and. oerr /= missing_r8 ) then
+
+    call create_obs_type(lat(n), lon(n), elev(n), VERTISSURFACE, tdew(n), &
+                         METAR_DEWPOINT_2_METER, oerr, oday, osec, qc, obs)
+    call append_obs_to_seq(obs_seq, obs)
+
+    end if
+
+!    print*, 'temp (C), rh (%), oerr:  ', tair(n)-273.15_r8, rh*100.0_r8, oerr
+
+  end if
+
+  nused = nused + 1
+  latu(nused) = lat(n)
+  lonu(nused) = lon(n)
+
+end do obsloop
+
+if ( get_num_obs(obs_seq) > 0 )  call write_obs_seq(obs_seq, surface_out_file)
+call check( nf90_close(ncid) )
+
+! end of main program
+
+contains
+

@@ Diff output truncated at 40000 characters. @@


More information about the Dart-dev mailing list