[Dart-dev] [6180] DART/branches/development/observations/MADIS: Add option to read sigwinds on height levels (which are in the MADIS files )
nancy at ucar.edu
nancy at ucar.edu
Wed May 29 10:05:50 MDT 2013
Revision: 6180
Author: nancy
Date: 2013-05-29 10:05:50 -0600 (Wed, 29 May 2013)
Log Message:
-----------
Add option to read sigwinds on height levels (which are in the MADIS files)
instead of the default which is on pressure levels. Also add option to
enable a namelist to set the defaults for sigwinds or no, and which of the
moisture obs to generate. It is off by default so the function is fully
backwards compatible, but if enabled it will no longer prompt at the command
line for input.
Modified Paths:
--------------
DART/branches/development/observations/MADIS/convert_madis_rawin.f90
DART/branches/development/observations/MADIS/work/input.nml
Property Changed:
----------------
DART/branches/development/observations/MADIS/
-------------- next part --------------
Property changes on: DART/branches/development/observations/MADIS
___________________________________________________________________
Added: svn:mergeinfo
+ /DART/trunk/observations/MADIS:4680-5928
Modified: DART/branches/development/observations/MADIS/convert_madis_rawin.f90
===================================================================
--- DART/branches/development/observations/MADIS/convert_madis_rawin.f90 2013-05-29 16:04:07 UTC (rev 6179)
+++ DART/branches/development/observations/MADIS/convert_madis_rawin.f90 2013-05-29 16:05:50 UTC (rev 6180)
@@ -28,11 +28,19 @@
!
! keep original obs times, make source for all converters as similar
! as possbile. nancy collins, ncar/image 26 march 2010
+!
+! add code to use the (new?) pressure levels for vertical coord in the
+! significant wind obs. also make an optional namelist section; if
+! enabled, the selection of mandatory/significant/both and height/pressure
+! can be from a namelist and avoid the prompts and reads from the console.
+! nancy collins, ncar/image, 23 mar 2012
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
use types_mod, only : r8, missing_r8
-use utilities_mod, only : nc_check, initialize_utilities, finalize_utilities
+use utilities_mod, only : nc_check, initialize_utilities, finalize_utilities, &
+ find_namelist_in_file, check_namelist_read, &
+ do_nml_file, do_nml_term, logfileunit, nmlfileunit
use time_manager_mod, only : time_type, set_calendar_type, set_date, &
increment_time, get_time, operator(-), GREGORIAN
use location_mod, only : VERTISSURFACE, VERTISPRESSURE, VERTISHEIGHT
@@ -65,28 +73,20 @@
character(len=19), parameter :: rawin_in_file = 'rawin_input.nc'
character(len=129), parameter :: rawin_out_file = 'obs_seq.rawin'
-! 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 :: num_copies = 1, & ! number of copies in sequence
num_qc = 1 ! number of QC entries
-integer :: oday, osec, nman, nsig, nsound, nused, &
+integer :: oday, osec, nman, nsig, nsound, nused, io, iunit, &
nmaxml, nmaxsw, nmaxst, maxobs, nvars_man, nvars_sigt, k, n, i, ncid
integer, allocatable :: obscnt(:)
-logical :: fexist, sigwnd, sigtmp, first_obs
+logical :: fexist, first_obs
real(r8) :: otime, lat, lon, elev, uwnd, vwnd, qobs, qsat, dptk, oerr, &
pres_miss, wdir_miss, wspd_miss, tair_miss, tdew_miss, prespa, &
- qc, altim, rh, qerr, hght_miss, HGHT_THRESHOLD ! , time_miss
+ qc, altim, rh, qerr, hght_miss
+real(r8), parameter :: HGHT_THRESHOLD = 40000.0_r8
real(r8), allocatable :: latu(:), lonu(:)
real(r8), allocatable :: pres(:), wdir(:), wspd(:), tair(:), tdew(:), hght(:)
@@ -97,11 +97,54 @@
type(time_type) :: comp_day0, time_obs, prev_time
+! 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. more recent
+! MADIS files contain the significant winds in both height and pressure vertical coordinates
+! and the last variable controls which one is used.
+! if you are using the namelist (see 'use_namelist' below) you can change these at runtime.
+! if not, they are compile-time settings.
+logical :: LH_err = .false.
+logical :: include_specific_humidity = .true.
+logical :: include_relative_humidity = .false.
+logical :: include_dewpoint = .false.
+logical :: use_input_qc = .true.
+logical :: wind_use_vert_pressure = .true.
+
+! mandatory levels are always converted. significant levels are controlled
+! by these vars. if you are using the namelist they are set by the
+! namelist in input.nml. if not, you are prompted at runtime for T/F responses.
+logical :: do_significant_level_temps = .true.
+logical :: do_significant_level_winds = .true.
+
+logical :: use_namelist = .false.
+
+! THIS WILL NOT BE READ IN UNLESS use_namelist IS SET TO .true.
+! if true, the code will NOT prompt for input from the console.
+
+namelist /convert_madis_rawin_nml/ &
+ do_significant_level_temps, do_significant_level_winds, &
+ wind_use_vert_pressure, LH_err, include_specific_humidity, &
+ include_relative_humidity, include_dewpoint, use_input_qc
+
+!-------- start of executable code -----------
+
call initialize_utilities('convert_madis_rawin')
-! prompt for optional significant level values
-print*,'Include significant level winds, temperature?: '
-read*, sigwnd, sigtmp
+if (use_namelist) then
+ ! Read the namelist entry
+ call find_namelist_in_file("input.nml", "convert_madis_rawin_nml", iunit)
+ read(iunit, nml = convert_madis_rawin_nml, iostat = io)
+ call check_namelist_read(iunit, io, "convert_madis_rawin_nml")
+
+ ! Record the namelist values used for the run ...
+ if (do_nml_file()) write(nmlfileunit, nml=convert_madis_rawin_nml)
+ if (do_nml_term()) write( * , nml=convert_madis_rawin_nml)
+else
+ ! prompt for optional significant level values
+ print*,'Include significant level winds, temperature?: '
+ read*, do_significant_level_winds, do_significant_level_temps
+endif
! shouldn't need to do this now - the nf90 open will catch it below
inquire(file=trim(rawin_in_file), exist=fexist) ! check for drop file
@@ -132,14 +175,18 @@
if ( obscnt(n) > nmaxml .and. obscnt(n) < 25 ) nmaxml = obscnt(n)
end do
-if ( sigwnd ) then
- call getvar_int(ncid, "numSigW", obscnt)
+if ( do_significant_level_winds ) then
+ if ( wind_use_vert_pressure ) then
+ call getvar_int(ncid, "numSigPresW", obscnt)
+ else
+ call getvar_int(ncid, "numSigW", obscnt)
+ endif
do n = 1, nsound
if ( obscnt(n) > nmaxsw .and. obscnt(n) < 150 ) nmaxsw = obscnt(n)
end do
endif
-if ( sigtmp ) then
+if ( do_significant_level_temps ) then
call getvar_int(ncid, "numSigT", obscnt)
do n = 1, nsound
if ( obscnt(n) > nmaxst .and. obscnt(n) < 150 ) nmaxst = obscnt(n)
@@ -253,7 +300,7 @@
qc_wdir = 0 ; qc_wspd = 0
endif
- if ( pres(1) /= pres_miss .and. qc_pres(1) == 0 ) then
+ if ( pres(1) /= pres_miss .and. qc_pres(1) == 0 .and. elev < 9999.0) then
altim = compute_altimeter(pres(1), elev)
oerr = rawin_pres_error(pres_alt_to_pres(elev) * 0.01_r8)
@@ -392,7 +439,7 @@
! If desired, read the significant-level temperature data, write to obs_seq
call getvar_int_1d_1val(ncid, "numSigT", n, nsig )
- if ( sigtmp .and. nsig <= nmaxst ) then
+ if ( do_significant_level_temps .and. nsig <= nmaxst .and. nsig > 0) then
allocate(pres(nsig)) ; allocate(tair(nsig)) ; allocate(tdew(nsig))
allocate(qc_pres(nsig)) ; allocate(qc_tair(nsig)) ; allocate(qc_tdew(nsig))
@@ -513,12 +560,77 @@
endif
+ ! If desired, read the pressure significant-level wind data, write to obs_seq
+ if ( wind_use_vert_pressure ) then
+ call getvar_int_1d_1val(ncid, "numSigPresW", n, nsig )
+ else
+ nsig = 0
+ endif
- ! If desired, read the significant-level wind data, write to obs_seq
- call getvar_int_1d_1val(ncid, "numSigW", n, nsig )
-
- if ( sigwnd .and. nsig <= nmaxsw ) then
+ if ( do_significant_level_winds .and. nsig > 0 .and. nsig <= nmaxsw) then
+ allocate(pres(nsig)) ; allocate(wdir(nsig)) ; allocate(wspd(nsig))
+ allocate(qc_pres(nsig)) ; allocate(qc_wdir(nsig)) ; allocate(qc_wspd(nsig))
+
+ ! read significant level data
+ call getvar_real_2d(ncid, "prSigW", n, nsig, pres, pres_miss)
+ call getvar_real_2d(ncid, "wdSigPrW", n, nsig, wdir, wdir_miss)
+ call getvar_real_2d(ncid, "wsSigPrW", n, nsig, wspd, wspd_miss)
+
+ if (use_input_qc) then
+ call get_or_fill_QC_2d(ncid, "prSigWQCR", n, nsig, qc_pres)
+ call get_or_fill_QC_2d(ncid, "wdSigPrWQCR", n, nsig, qc_wdir)
+ call get_or_fill_QC_2d(ncid, "wsSigPrWQCR", n, nsig, qc_wspd)
+ else
+ qc_pres = 0
+ qc_wdir = 0
+ qc_wspd = 0
+ endif
+
+ do k = 1, nsig
+
+ prespa = pres(k) * 100.0_r8
+
+ ! add data to the observation sequence here.
+ if ( wdir(k) /= wdir_miss .and. qc_wdir(k) == 0 .and. &
+ wspd(k) /= wspd_miss .and. qc_wspd(k) == 0 .and. &
+ pres(k) /= pres_miss .and. qc_pres(k) == 0 ) then
+
+
+ call wind_dirspd_to_uv(wdir(k), wspd(k), uwnd, vwnd)
+
+ ! expects press in hPa, which is what we already have
+ oerr = rawin_wind_error(pres(k))
+ if ( abs(uwnd) <= 150.0_r8 .and. &
+ abs(vwnd) <= 150.0_r8 .and. oerr /= missing_r8 ) then
+
+ call create_3d_obs(lat, lon, prespa, VERTISPRESSURE, uwnd, &
+ RADIOSONDE_U_WIND_COMPONENT, oerr, oday, osec, qc, obs)
+ call add_obs_to_seq(obs_seq, obs, time_obs, prev_obs, prev_time, first_obs)
+
+ call create_3d_obs(lat, lon, prespa, VERTISPRESSURE, vwnd, &
+ RADIOSONDE_V_WIND_COMPONENT, oerr, oday, osec, qc, obs)
+ call add_obs_to_seq(obs_seq, obs, time_obs, prev_obs, prev_time, first_obs)
+
+ endif
+
+ endif
+
+ end do
+
+ deallocate(pres, wdir, wspd, qc_pres, qc_wdir, qc_wspd)
+
+ endif
+
+ ! If desired, read the height significant-level wind data, write to obs_seq
+ if (.not. wind_use_vert_pressure ) then
+ call getvar_int_1d_1val(ncid, "numSigW", n, nsig )
+ else
+ nsig = 0
+ endif
+
+ if ( do_significant_level_winds .and. nsig > 0 .and. nsig <= nmaxsw) then
+
allocate(hght(nsig)) ; allocate(wdir(nsig)) ; allocate(wspd(nsig))
allocate(qc_hght(nsig)) ; allocate(qc_wdir(nsig)) ; allocate(qc_wspd(nsig))
@@ -528,6 +640,7 @@
call getvar_real_2d(ncid, "wsSigW", n, nsig, wspd, wspd_miss)
if (use_input_qc) then
+ call get_or_fill_QC_2d(ncid, "htSigWQCR", n, nsig, qc_wdir)
call get_or_fill_QC_2d(ncid, "wdSigWQCR", n, nsig, qc_wdir)
call get_or_fill_QC_2d(ncid, "wsSigWQCR", n, nsig, qc_wspd)
else
@@ -540,12 +653,19 @@
! add data to the observation sequence here.
if ( wdir(k) /= wdir_miss .and. qc_wdir(k) == 0 .and. &
wspd(k) /= wspd_miss .and. qc_wspd(k) == 0 .and. &
- hght(k) /= hght_miss ) then
+ hght(k) /= hght_miss .and. qc_hght(k) == 0) then
+ ! make sure elevation is a reasonable value if height comes
+ ! in as zero and use it instead of 0.0
+ if (hght(k) == 0.0 .and. elev < 9999.0) then
+ hght(k) = elev
+ else
+ cycle
+ endif
+
call wind_dirspd_to_uv(wdir(k), wspd(k), uwnd, vwnd)
! the pres_alt_to_pres() fails above 44km, so give any obs above
! this height the same wind error value as 40km.
- HGHT_THRESHOLD = 40000.0_r8
oerr = rawin_wind_error(pres_alt_to_pres(min(hght(k),HGHT_THRESHOLD)) * 0.01_r8)
if ( abs(uwnd) <= 150.0_r8 .and. &
abs(vwnd) <= 150.0_r8 .and. oerr /= missing_r8 ) then
Modified: DART/branches/development/observations/MADIS/work/input.nml
===================================================================
--- DART/branches/development/observations/MADIS/work/input.nml 2013-05-29 16:04:07 UTC (rev 6179)
+++ DART/branches/development/observations/MADIS/work/input.nml 2013-05-29 16:05:50 UTC (rev 6180)
@@ -12,6 +12,17 @@
'../../../obs_def/obs_def_rel_humidity_mod.f90',
/
+&convert_madis_rawin_nml
+ do_significant_level_temps = .true.
+ do_significant_level_winds = .true.
+ wind_use_vert_pressure = .true.
+ LH_err = .false.
+ include_specific_humidity = .true.
+ include_relative_humidity = .false.
+ include_dewpoint = .false.
+ use_input_qc = .true.
+ /
+
&obs_kind_nml
/
More information about the Dart-dev
mailing list