[Dart-dev] [4232] DART/trunk/models/wrf: Full set of files to support the Radar experiment, including
nancy at ucar.edu
nancy at ucar.edu
Thu Jan 28 12:31:54 MST 2010
Revision: 4232
Author: nancy
Date: 2010-01-28 12:31:54 -0700 (Thu, 28 Jan 2010)
Log Message:
-----------
Full set of files to support the Radar experiment, including
the additive noise and initial conditions perturbation executables.
Single line modification to the advance_model script - if the
radar noise script is found, call it at the appropriate place.
Modified Paths:
--------------
DART/trunk/models/wrf/shell_scripts/advance_model.csh
DART/trunk/models/wrf/work/input.nml
Added Paths:
-----------
DART/trunk/models/wrf/WRF_DART_utilities/add_pert_where_high_refl.f90
DART/trunk/models/wrf/WRF_DART_utilities/f2kcli.f90
DART/trunk/models/wrf/WRF_DART_utilities/grid_refl_obs.f90
DART/trunk/models/wrf/experiments/
DART/trunk/models/wrf/experiments/Radar/
DART/trunk/models/wrf/experiments/Radar/IC/
DART/trunk/models/wrf/experiments/Radar/IC/init_ideal.ncl
DART/trunk/models/wrf/experiments/Radar/IC/namelist.wps.17May1981.2km
DART/trunk/models/wrf/experiments/Radar/IC/prep_IC_pertsounding.csh
DART/trunk/models/wrf/experiments/Radar/IC/sounding_perturbation/
DART/trunk/models/wrf/experiments/Radar/IC/sounding_perturbation/compile_pert_sounding.csh
DART/trunk/models/wrf/experiments/Radar/IC/sounding_perturbation/input_sounding
DART/trunk/models/wrf/experiments/Radar/IC/sounding_perturbation/pert_sounding.f90
DART/trunk/models/wrf/experiments/Radar/IC/sounding_perturbation/pert_sounding.input
DART/trunk/models/wrf/experiments/Radar/IC/sounding_perturbation/pert_sounding_module.f90
DART/trunk/models/wrf/experiments/Radar/IC/sounding_perturbation/readme.altug
DART/trunk/models/wrf/experiments/Radar/IC/sounding_perturbation/test_sounding.output
DART/trunk/models/wrf/experiments/Radar/input.nml
DART/trunk/models/wrf/experiments/Radar/namelist.input
DART/trunk/models/wrf/experiments/Radar/obs/
DART/trunk/models/wrf/experiments/Radar/obs/README
DART/trunk/models/wrf/experiments/Radar/start_over.csh
DART/trunk/models/wrf/shell_scripts/add_noise.csh
DART/trunk/models/wrf/work/mkmf_add_pert_where_high_refl
DART/trunk/models/wrf/work/mkmf_grid_refl_obs
DART/trunk/models/wrf/work/path_names_add_pert_where_high_refl
DART/trunk/models/wrf/work/path_names_grid_refl_obs
Removed Paths:
-------------
DART/trunk/models/wrf/full_experiment/
-------------- next part --------------
Added: DART/trunk/models/wrf/WRF_DART_utilities/add_pert_where_high_refl.f90
===================================================================
--- DART/trunk/models/wrf/WRF_DART_utilities/add_pert_where_high_refl.f90 (rev 0)
+++ DART/trunk/models/wrf/WRF_DART_utilities/add_pert_where_high_refl.f90 2010-01-28 19:31:54 UTC (rev 4232)
@@ -0,0 +1,674 @@
+PROGRAM add_pert_where_high_refl
+
+! Add 3D random but smooth perturbations to WRF fields in/near where the observations
+! indicate high reflectivity values. The technique is somewhat like that of
+! Caya et al. 2005, Monthly Weather Review, 3081-3094.
+!
+! David Dowell 27 June 2007
+!
+! input parameters from command line:
+! (1) refl_ob_file -- name of text file containing WRF grid indices where observed reflectivity is high
+! (2) wrf_file -- path name of WRF netcdf file
+! (3) lh -- horizontal length scale (m) for perturbations
+! (4) lv -- vertical length scale (m) for perturbations
+! (5) u_sd -- std. dev. of u noise (m/s), before smoothing
+! (6) v_sd -- std. dev. of v noise (m/s), before smoothing
+! (7) w_sd -- std. dev. of w noise (m/s), before smoothing
+! (8) t_sd -- std. dev. of potential temperature noise (K), before smoothing
+! (9) td_sd -- std. dev. of dewpoint noise (K), before smoothing
+! (10) qv_sd -- std. dev. of water vapor mixing ratio noise, before smoothing
+! (input value is in g/kg, value after conversion is in kg/kg)
+!
+! output:
+
+use types_mod, only : r8, gravity, t_kelvin, ps0, gas_constant, gas_constant_v
+use utilities_mod, only : error_handler, E_ERR, E_MSG, initialize_utilities, &
+ timestamp, register_module, logfileunit, file_exist
+use random_nr_mod, only : random_seq_type, init_ran1
+use random_seq_mod, only : random_gaussian
+use netcdf
+use f2kcli
+
+implicit none
+
+! version controlled file description for error handling, do not edit
+character(len=128), parameter :: &
+ source = "$URL$", &
+ revision = "$Revision$", &
+ revdate = "$Date$"
+
+! command-line parameters
+character(len=129) :: refl_ob_file
+character(len=129) :: wrf_file
+real(r8) :: lh
+real(r8) :: lv
+real(r8) :: u_sd
+real(r8) :: v_sd
+real(r8) :: w_sd
+real(r8) :: t_sd
+real(r8) :: td_sd
+real(r8) :: qv_sd
+
+! local variables
+integer :: n_obs ! number of gridded reflectivity observations
+integer, allocatable :: i_ob(:) ! grid indices of observations
+integer, allocatable :: j_ob(:) ! " "
+integer, allocatable :: k_ob(:) ! " "
+real(r8), allocatable :: refl_ob(:) ! observed reflectivity (dBZ)
+
+real(r8), allocatable :: phb(:,:,:) ! base-state geopotential (m^2 s^-2)
+real(r8), allocatable :: ht(:,:,:) ! height MSL of mass grid points (m)
+real(r8), allocatable :: ht_u(:,:,:) ! height MSL of u grid points (m)
+real(r8), allocatable :: ht_v(:,:,:) ! height MSL of v grid points (m)
+real(r8), allocatable :: ht_w(:,:,:) ! height MSL of w grid points (m)
+real(r8), allocatable :: f(:,:,:) ! WRF field
+real(r8), allocatable :: f2(:,:,:) ! Extra WRF field
+real(r8), allocatable :: sd(:,:,:) ! standard deviations of grid-point noise
+
+real(r8), allocatable :: dnw(:) ! d(eta) values between full (w) levels
+real(r8), allocatable :: ph(:,:,:) ! perturbation geopotential (m^2 s^-2)
+real(r8), allocatable :: qv(:,:,:) ! water vapor (kg/kg)
+real(r8), allocatable :: t(:,:,:) ! perturbation potential temperature (K)
+real(r8) :: rho ! density (kg m^-3)
+real(r8), allocatable :: mu(:,:) ! perturbation dry air mass in column (Pa)
+real(r8), allocatable :: mub(:,:) ! base state dry air mass in column (Pa)
+real(r8) :: ph_e
+real(r8) :: qvf1
+real(r8), PARAMETER :: ts0 = 300.0_r8 ! Base potential temperature for all levels.
+real(r8), PARAMETER :: kappa = 2.0_r8/7.0_r8 ! gas_constant / cp
+real(r8), PARAMETER :: rd_over_rv = gas_constant / gas_constant_v
+real(r8), PARAMETER :: cpovcv = 1.4_r8 ! cp / (cp - gas_constant)
+real(r8), allocatable :: p(:,:,:) ! pressure (mb)
+
+
+character(len=8) :: crdate ! needed by F90 DATE_AND_TIME intrinsic
+character(len=10) :: crtime ! needed by F90 DATE_AND_TIME intrinsic
+character(len=5) :: crzone ! needed by F90 DATE_AND_TIME intrinsic
+integer, dimension(8) :: values ! needed by F90 DATE_AND_TIME intrinsic
+real(r8) :: dx, dy ! horizontal grid spacings (m)
+integer :: bt, sn, we ! WRF grid dimensions
+integer :: i, j, k, o
+
+! netcdf stuff
+integer :: var_id, ncid, ierr
+character(len=80) :: varname
+
+! f2kcli stuff
+integer :: status, length
+character(len=120) :: string
+
+! random number generator stuff
+type (random_seq_type) :: rs
+
+
+ ! Get command-line parameters, using the F2KCLI interface. See f2kcli.f90 for details.
+
+if( COMMAND_ARGUMENT_COUNT() .ne. 10 ) then
+ print*, 'INCORRECT # OF ARGUMENTS ON COMMAND LINE: ', COMMAND_ARGUMENT_COUNT()
+ call exit(1)
+else
+
+ call GET_COMMAND_ARGUMENT(1,refl_ob_file,length,status)
+ if( status .ne. 0 ) then
+ print*, 'refl_ob_file NOT RETRIEVED FROM COMMAND LINE: ', status
+ call exit(1)
+ endif
+
+ call GET_COMMAND_ARGUMENT(2,wrf_file,length,status)
+ if( status .ne. 0 ) then
+ print*, 'wrf_file NOT RETRIEVED FROM COMMAND LINE: ', status
+ call exit(1)
+ endif
+
+ call GET_COMMAND_ARGUMENT(3,string,length,status)
+ if( status .ne. 0 ) then
+ print*, 'lh NOT RETRIEVED FROM COMMAND LINE: ', status
+ call exit(1)
+ else
+ read(string,*) lh
+ endif
+
+ call GET_COMMAND_ARGUMENT(4,string,length,status)
+ if( status .ne. 0 ) then
+ print*, 'lv NOT RETRIEVED FROM COMMAND LINE: ', status
+ call exit(1)
+ else
+ read(string,*) lv
+ endif
+
+ call GET_COMMAND_ARGUMENT(5,string,length,status)
+ if( status .ne. 0 ) then
+ print*, 'u_sd NOT RETRIEVED FROM COMMAND LINE: ', status
+ call exit(1)
+ else
+ read(string,*) u_sd
+ endif
+
+ call GET_COMMAND_ARGUMENT(6,string,length,status)
+ if( status .ne. 0 ) then
+ print*, 'v_sd NOT RETRIEVED FROM COMMAND LINE: ', status
+ call exit(1)
+ else
+ read(string,*) v_sd
+ endif
+
+ call GET_COMMAND_ARGUMENT(7,string,length,status)
+ if( status .ne. 0 ) then
+ print*, 'w_sd NOT RETRIEVED FROM COMMAND LINE: ', status
+ call exit(1)
+ else
+ read(string,*) w_sd
+ endif
+
+ call GET_COMMAND_ARGUMENT(8,string,length,status)
+ if( status .ne. 0 ) then
+ print*, 't_sd NOT RETRIEVED FROM COMMAND LINE: ', status
+ call exit(1)
+ else
+ read(string,*) t_sd
+ endif
+
+ call GET_COMMAND_ARGUMENT(9,string,length,status)
+ if( status .ne. 0 ) then
+ print*, 'td_sd NOT RETRIEVED FROM COMMAND LINE: ', status
+ call exit(1)
+ else
+ read(string,*) td_sd
+ endif
+
+ call GET_COMMAND_ARGUMENT(10,string,length,status)
+ if( status .ne. 0 ) then
+ print*, 'qv_sd NOT RETRIEVED FROM COMMAND LINE: ', status
+ call exit(1)
+ else
+ read(string,*) qv_sd
+ endif
+
+endif
+
+qv_sd = 0.001_r8 * qv_sd ! convert g/kg to kg/kg
+
+
+! Read locations where high reflectivity was observed.
+
+! first, count observations
+
+open(unit=11, file=refl_ob_file, status='old')
+n_obs = 0
+998 read(11,*,end=999)
+ n_obs = n_obs + 1
+go to 998
+999 close(11)
+
+! now allocate storage and read the observations
+
+allocate(i_ob(n_obs))
+allocate(j_ob(n_obs))
+allocate(k_ob(n_obs))
+allocate(refl_ob(n_obs))
+
+open(unit=11, file=refl_ob_file, status='old')
+do o=1, n_obs
+ read(11,*) i_ob(o), j_ob(o), k_ob(o), refl_ob(o)
+enddo
+close(11)
+
+
+! Open WRF file and obtain miscellaneous values.
+
+call check ( nf90_open(wrf_file, NF90_WRITE, ncid) )
+
+call check ( nf90_inq_dimid(ncid, 'bottom_top', var_id) )
+call check ( nf90_inquire_dimension(ncid, var_id, varname, bt) )
+
+call check ( nf90_inq_dimid(ncid, 'south_north', var_id) )
+call check ( nf90_inquire_dimension(ncid, var_id, varname, sn) )
+
+call check ( nf90_inq_dimid(ncid, 'west_east', var_id) )
+call check ( nf90_inquire_dimension(ncid, var_id, varname, we) )
+
+call check( nf90_get_att(ncid, nf90_global, 'DX', dx) )
+call check( nf90_get_att(ncid, nf90_global, 'DY', dy) )
+
+
+! Read WRF base-state geopotential height field and compute height (m MSL)
+! of each grid point.
+
+allocate(phb(we,sn,bt+1))
+allocate(ht(we,sn,bt))
+allocate(ht_u(we+1,sn,bt))
+allocate(ht_v(we,sn+1,bt))
+allocate(ht_w(we,sn,bt+1))
+
+call check ( nf90_inq_varid(ncid, 'PHB', var_id))
+call check ( nf90_get_var(ncid, var_id, phb, start = (/ 1, 1, 1, 1/)))
+
+do k=1, bt
+ do j=1, sn
+ do i=1, we
+ ht(i,j,k) = ( phb(i,j,k) + phb(i,j,k+1) ) / (2.0_r8*gravity)
+ enddo
+ enddo
+enddo
+do k=1, bt
+ do j=1, sn
+ do i=2, we
+ ht_u(i,j,k) = ( phb(i-1,j,k) + phb(i-1,j,k+1) + phb(i,j,k) + phb(i,j,k+1) ) / (4.0_r8*gravity)
+ enddo
+ ht_u(1,j,k) = ht_u(2,j,k)
+ ht_u(we+1,j,k) = ht_u(we,j,k)
+ enddo
+enddo
+do k=1, bt
+ do i=1, we
+ do j=2, sn
+ ht_v(i,j,k) = ( phb(i,j-1,k) + phb(i,j-1,k+1) + phb(i,j,k) + phb(i,j,k+1) ) / (4.0_r8*gravity)
+ enddo
+ ht_v(i,1,k) = ht_v(i,2,k)
+ ht_v(i,sn+1,k) = ht_v(i,sn,k)
+ enddo
+enddo
+do k=1, bt+1
+ do j=1, sn
+ do i=1, we
+ ht_w(i,j,k) = phb(i,j,k) / gravity
+ enddo
+ enddo
+enddo
+
+
+! Initialize random number generator with a seed based on the milliseconds
+! portion of the current time.
+
+call date_and_time(crdate,crtime,crzone,values)
+call init_ran1(rs, -int(values(8)))
+
+
+! Add perturbations.
+
+if (u_sd .gt. 0.0_r8) then
+
+ allocate(f(we+1,sn,bt))
+ call check ( nf90_inq_varid(ncid, 'U', var_id))
+ call check ( nf90_get_var(ncid, var_id, f, start = (/ 1, 1, 1, 1/)))
+ allocate(sd(we+1,sn,bt))
+
+ sd(:,:,:) = 0.0_r8
+ do o=1, n_obs
+ sd(i_ob(o), j_ob(o), k_ob(o)) = u_sd
+ sd(i_ob(o)+1, j_ob(o), k_ob(o)) = u_sd
+ enddo
+
+ call add_smooth_perturbations(f, sd, we+1, sn, bt, lh, lv, dx, dy, ht_u)
+
+ call check ( nf90_put_var(ncid, var_id, f, start = (/ 1, 1, 1, 1/)))
+ deallocate(f)
+ deallocate(sd)
+
+end if
+
+
+if (v_sd .gt. 0.0_r8) then
+
+ allocate(f(we,sn+1,bt))
+ call check ( nf90_inq_varid(ncid, 'V', var_id))
+ call check ( nf90_get_var(ncid, var_id, f, start = (/ 1, 1, 1, 1/)))
+ allocate(sd(we,sn+1,bt))
+
+ sd(:,:,:) = 0.0_r8
+ do o=1, n_obs
+ sd(i_ob(o), j_ob(o), k_ob(o)) = v_sd
+ sd(i_ob(o), j_ob(o)+1, k_ob(o)) = v_sd
+ enddo
+
+ call add_smooth_perturbations(f, sd, we, sn+1, bt, lh, lv, dx, dy, ht_v)
+
+ call check ( nf90_put_var(ncid, var_id, f, start = (/ 1, 1, 1, 1/)))
+ deallocate(f)
+ deallocate(sd)
+
+end if
+
+
+! note: Perturbing w is not advised currently because there is no enforcement
+! that w perturbations should be small near the lower and upper boundaries.
+if (w_sd .gt. 0.0_r8) then
+
+ allocate(f(we,sn,bt+1))
+ call check ( nf90_inq_varid(ncid, 'W', var_id))
+ call check ( nf90_get_var(ncid, var_id, f, start = (/ 1, 1, 1, 1/)))
+ allocate(sd(we,sn,bt+1))
+
+ sd(:,:,:) = 0.0_r8
+ do o=1, n_obs
+ sd(i_ob(o), j_ob(o), k_ob(o)) = w_sd
+ sd(i_ob(o), j_ob(o), k_ob(o)+1) = w_sd
+ enddo
+
+ call add_smooth_perturbations(f, sd, we, sn, bt+1, lh, lv, dx, dy, ht_w)
+
+ call check ( nf90_put_var(ncid, var_id, f, start = (/ 1, 1, 1, 1/)))
+ deallocate(f)
+ deallocate(sd)
+
+end if
+
+
+if (t_sd .gt. 0.0_r8) then
+
+ allocate(f(we,sn,bt))
+ call check ( nf90_inq_varid(ncid, 'T', var_id))
+ call check ( nf90_get_var(ncid, var_id, f, start = (/ 1, 1, 1, 1/)))
+ allocate(sd(we,sn,bt))
+
+ sd(:,:,:) = 0.0_r8
+ do o=1, n_obs
+ sd(i_ob(o), j_ob(o), k_ob(o)) = t_sd
+ enddo
+
+ call add_smooth_perturbations(f, sd, we, sn, bt, lh, lv, dx, dy, ht)
+
+ call check ( nf90_put_var(ncid, var_id, f, start = (/ 1, 1, 1, 1/)))
+ deallocate(f)
+ deallocate(sd)
+
+end if
+
+
+! note: Perturbing dewpoint can produce supersaturation.
+if (td_sd .gt. 0.0_r8) then
+
+ allocate(dnw(bt))
+ allocate(ph(we,sn,bt+1))
+ allocate(qv(we,sn,bt))
+ allocate(t(we,sn,bt))
+ allocate(mu(we,sn))
+ allocate(mub(we,sn))
+ allocate(p(we,sn,bt))
+ call check ( nf90_inq_varid(ncid, 'DNW', var_id))
+ call check ( nf90_get_var(ncid, var_id, dnw, start = (/ 1, 1/)))
+ call check ( nf90_inq_varid(ncid, 'PH', var_id))
+ call check ( nf90_get_var(ncid, var_id, ph, start = (/ 1, 1, 1, 1/)))
+ call check ( nf90_inq_varid(ncid, 'QVAPOR', var_id))
+ call check ( nf90_get_var(ncid, var_id, qv, start = (/ 1, 1, 1, 1/)))
+ call check ( nf90_inq_varid(ncid, 'T', var_id))
+ call check ( nf90_get_var(ncid, var_id, t, start = (/ 1, 1, 1, 1/)))
+ call check ( nf90_inq_varid(ncid, 'MU', var_id))
+ call check ( nf90_get_var(ncid, var_id, mu, start = (/ 1, 1, 1/)))
+ call check ( nf90_inq_varid(ncid, 'MUB', var_id))
+ call check ( nf90_get_var(ncid, var_id, mub, start = (/ 1, 1, 1/)))
+ allocate(f(we,sn,bt))
+ allocate(f2(we,sn,bt))
+! f2 is the sensible temperature array
+ allocate(sd(we,sn,bt))
+
+ ! compute pressure
+ ! see model_rho_t and model_pressure_t functions in model_mod.f90
+ do k=1, bt
+ do j=1, sn
+ do i=1, we
+ ph_e = ( (ph(i,j,k+1) + phb(i,j,k+1)) - (ph(i,j,k) + phb(i,j,k)) ) / dnw(k)
+ rho = - ( mub(i,j)+mu(i,j) ) / ph_e
+ qvf1 = 1.0_r8 + qv(i,j,k) / rd_over_rv
+ p(i,j,k) = 0.01_r8 * ps0 * ( (gas_constant*(ts0+t(i,j,k))*qvf1) / (ps0/rho) )**cpovcv
+ f2(i,j,k) = (ts0 + t(i,j,k))*(100.0*p(i,j,k)/ps0)**kappa
+ enddo
+ enddo
+ enddo
+
+ ! compute dewpoint
+ do k=1, bt
+ do j=1, sn
+ do i=1, we
+ call compute_td(f(i,j,k), p(i,j,k), qv(i,j,k))
+ enddo
+ enddo
+ enddo
+
+ ! perturb dewpoint
+ sd(:,:,:) = 0.0_r8
+ do o=1, n_obs
+ sd(i_ob(o), j_ob(o), k_ob(o)) = td_sd
+ enddo
+
+
+ call add_smooth_perturbations(f, sd, we, sn, bt, lh, lv, dx, dy, ht)
+
+! check to make sure the perturbed dewpoint is <= temperature + 4 K
+ do k=1, bt
+ do j=1, sn
+ do i=1, we
+ if ( f(i,j,k) .gt. f2(i,j,k)+4.0 ) then
+! write(*,*) 'supersaturation violation i,j,k ', i, j, k, f(i,j,k), f2(i,j,k)
+ f(i,j,k) = f2(i,j,k)+4.0
+ end if
+ enddo
+ enddo
+ enddo
+
+ ! compute qv
+ do k=1, bt
+ do j=1, sn
+ do i=1, we
+ call compute_qv(qv(i,j,k), p(i,j,k), f(i,j,k))
+ enddo
+ enddo
+ enddo
+
+ call check ( nf90_inq_varid(ncid, 'QVAPOR', var_id))
+ call check ( nf90_put_var(ncid, var_id, qv, start = (/ 1, 1, 1, 1/)))
+ deallocate(dnw)
+ deallocate(ph)
+ deallocate(qv)
+ deallocate(t)
+ deallocate(mu)
+ deallocate(mub)
+ deallocate(p)
+ deallocate(f)
+ deallocate(f2)
+ deallocate(sd)
+
+end if
+
+
+! note: Since negative qv values are set to 0 after the perturbations are added,
+! the following procedure produces a net increase in qv in the domain.
+if (qv_sd .gt. 0.0_r8) then
+
+ allocate(f(we,sn,bt))
+ call check ( nf90_inq_varid(ncid, 'QVAPOR', var_id))
+ call check ( nf90_get_var(ncid, var_id, f, start = (/ 1, 1, 1, 1/)))
+ allocate(sd(we,sn,bt))
+
+ sd(:,:,:) = 0.0_r8
+ do o=1, n_obs
+ sd(i_ob(o), j_ob(o), k_ob(o)) = qv_sd
+ enddo
+
+ call add_smooth_perturbations(f, sd, we, sn, bt, lh, lv, dx, dy, ht)
+
+ f(:,:,:) = max(0.0_r8, f(:,:,:)) ! require qv to be nonnegative
+
+ call check ( nf90_put_var(ncid, var_id, f, start = (/ 1, 1, 1, 1/)))
+ deallocate(f)
+ deallocate(sd)
+
+end if
+
+
+! Close file and deallocate arrays.
+
+ierr = NF90_close(ncid)
+
+deallocate(i_ob)
+deallocate(j_ob)
+deallocate(k_ob)
+deallocate(refl_ob)
+deallocate(phb)
+deallocate(ht)
+deallocate(ht_u)
+deallocate(ht_v)
+deallocate(ht_w)
+
+
+
+contains
+
+
+ ! Internal subroutine - checks error status after each netcdf, prints
+ ! text message each time an error code is returned.
+ subroutine check(istatus)
+ integer, intent ( in) :: istatus
+ if(istatus /= nf90_noerr) call error_handler(E_ERR,'add_pert_where_high_refl', &
+ trim(nf90_strerror(istatus)), source, revision, revdate)
+ end subroutine check
+
+!------------------------------------------------------------------------------------
+
+ ! Compute dewpoint (in Kelvin) from the specified values of pressure and water vapor mixing ratio.
+ ! Author: David Dowell
+ ! Date: July 9, 2007
+
+ subroutine compute_td(td, p, qv)
+ implicit none
+
+!-- returned parameter
+ real(r8), intent(out) :: td ! dewpoint (K)
+
+!-- passed parameters
+ real(r8), intent(in) :: p ! pressure (mb)
+ real(r8), intent(in) :: qv ! water vapor mixing ratio (kg/kg)
+
+!-- local variables
+ real(r8) :: e ! water vapor pressure (mb)
+ real(r8) :: qv_ckd ! checked mixing ratio
+ real(r8), PARAMETER :: e_min = 0.001_r8 ! threshold for minimum vapor pressure (mb),
+ ! to avoid problems near zero in Bolton's equation
+ real(r8), PARAMETER :: qv_min = 1.0e-12 ! threshold for minimum water vapor mixing ratio
+ real(r8), PARAMETER :: qv_max = 0.050 ! threshold for maximum water vapor mixing ratio
+
+! ensure qv is a reasonable number
+! qv_ckd = max (qv, qv_min)
+! qv_ckd = min (qv_ckd, qv_max)
+ qv_ckd = qv
+ e = qv_ckd * p / (0.622_r8 + qv_ckd) ! vapor pressure
+ e = max(e, e_min) ! avoid problems near zero
+ td = t_kelvin + (243.5_r8 / ((17.67_r8 / log(e/6.112_r8)) - 1.0_r8) ) ! Bolton's approximation
+
+ end subroutine compute_td
+
+!------------------------------------------------------------------------------------
+
+ ! Compute water vapor mixing ratio (in kg/kg) from the specified values of pressure and dewpoint.
+ ! Author: David Dowell
+ ! Date: July 9, 2007
+
+ subroutine compute_qv(qv, p, td)
+ implicit none
+
+!-- returned parameter
+ real(r8), intent(out) :: qv ! water vapor mixing ratio (kg/kg)
+
+!-- passed parameters
+ real(r8), intent(in) :: p ! pressure (mb)
+ real(r8), intent(in) :: td ! dewpoint (K)
+
+!-- local variables
+ real(r8) :: tdc ! dewpoint (Celsius)
+ real(r8) :: e ! water vapor pressure (mb)
+
+ tdc = td - t_kelvin
+ e = 6.112_r8 * exp(17.67_r8 * tdc / (tdc+243.5_r8) ) ! Bolton's approximation
+ qv = 0.622_r8 * e / (p-e)
+
+ return
+
+ end subroutine compute_qv
+
+!------------------------------------------------------------------------------------
+
+ ! Add smooth perturbations to an array. The technique is based on
+ ! Caya et al. 2005, Monthly Weather Review, 3081-3094.
+ ! Author: David Dowell
+ ! Date: July 9, 2007
+
+ subroutine add_smooth_perturbations(f, sd, nx, ny, nz, lh, lv, dx, dy, ht)
+ implicit none
+
+!-- passed parameters
+ integer, intent(in) :: nx, ny, nz ! grid dimensions
+ real(r8), intent(in) :: lh, lv ! horizontal and vertical length scales (m)
+ real(r8), intent(in) :: dx, dy ! horizontal grid spacings (m)
+ real(r8), intent(in) :: ht(nx,ny,nz) ! heights MSL of f and sd grid points (m)
+ real(r8), intent(in) :: sd(nx,ny,nz) ! standard deviation of grid-point noise
+
+!-- passed and returned variable
+ real(r8), intent(inout) :: f(nx,ny,nz) ! field to be perturbed
+
+!-- local variables
+
+ real(r8) :: r(nx,ny,nz) ! realization of random, normally distributed noise
+ integer :: i, i0, j, j0, k, k0 ! grid indices
+ integer :: i1, i2, j1, j2, k1, k2 ! more grid indices
+ real(r8) :: rlh, rlv ! reciprocals of lh and lv
+ integer, parameter :: nl = 5 ! number of length scales for computing exponential weight
+
+
+ rlh = 1.0_r8 / lh
+ rlv = 1.0_r8 / lv
+
+! generate random, normally-distributed gridpoint noise
+
+ r(:,:,:) = 0.0_r8
+ do k0=1, nz
+ do j0=1, ny
+ do i0=1, nx
+ if (sd(i0,j0,k0) .ne. 0.0_r8) then
+ r(i0,j0,k0) = random_gaussian(rs, 0.0_r8, sd(i0,j0,k0))
+ endif
+ enddo
+ enddo
+ enddo
+
+! smooth the perturbations with an inverse exponential function
+
+ do k0=1, nz
+ do j0=1, ny
+ do i0=1, nx
+
+ if (r(i0,j0,k0).ne.0.0_r8) then
+
+ i1 = max(1, nint(i0-nl*lh/dx))
+ i2 = min(nx, nint(i0+nl*lh/dx))
+ j1 = max(1, nint(j0-nl*lh/dy))
+ j2 = min(ny, nint(j0+nl*lh/dy))
+ k1 = k0
+ do while ( (k1.gt.1) .and. ( ht(i0,j0,k1) .gt. (ht(i0,j0,k0)-nl*lv) ) )
+ k1 = k1 - 1
+ enddo
+ k2 = k0
+ do while ( (k2.lt.nz) .and. ( ht(i0,j0,k2) .lt. (ht(i0,j0,k0)+nl*lv) ) )
+ k2 = k2 + 1
+ enddo
+
+ do k=k1, k2
+ do j=j1, j2
+ do i=i1, i2
+ f(i,j,k) = f(i,j,k) &
+ + r(i0,j0,k0)*exp( -dx*abs(i0-i)*rlh &
+ -dy*abs(j0-j)*rlh &
+ -abs(ht(i0,j0,k0)-ht(i,j,k))*rlv )
+ enddo
+ enddo
+ enddo
+
+ endif
+
+ enddo
+ enddo
+ enddo
+
+ end subroutine add_smooth_perturbations
+
+
+END PROGRAM add_pert_where_high_refl
Property changes on: DART/trunk/models/wrf/WRF_DART_utilities/add_pert_where_high_refl.f90
___________________________________________________________________
Added: svn:keywords
+ Date Revision Author HeadURL Id
Added: DART/trunk/models/wrf/WRF_DART_utilities/f2kcli.f90
===================================================================
--- DART/trunk/models/wrf/WRF_DART_utilities/f2kcli.f90 (rev 0)
+++ DART/trunk/models/wrf/WRF_DART_utilities/f2kcli.f90 2010-01-28 19:31:54 UTC (rev 4232)
@@ -0,0 +1,187 @@
+! F2KCLI : Fortran 200x Command Line Interface
+! copyright Interactive Software Services Ltd. 2002
+! For conditions of use see manual.txt
+!
+! Platform : Mac OS/X
+! Compiler : Absoft Pro Fortran
+! To compile : f95 -c f2kcli.f90
+! Implementer : Lawson B. Wakefield, I.S.S. Ltd.
+! Date : June 2002
+!
+ MODULE F2KCLI
+!
+ CONTAINS
+!
+ SUBROUTINE GET_COMMAND(COMMAND,LENGTH,STATUS)
+!
+! Description. Returns the entire command by which the program was
+! invoked.
+!
+! Class. Subroutine.
+!
+! Arguments.
+! COMMAND (optional) shall be scalar and of type default character.
+! It is an INTENT(OUT) argument. It is assigned the entire command
+! by which the program was invoked. If the command cannot be
+! determined, COMMAND is assigned all blanks.
+! LENGTH (optional) shall be scalar and of type default integer. It is
+! an INTENT(OUT) argument. It is assigned the significant length
+! of the command by which the program was invoked. The significant
+! length may include trailing blanks if the processor allows commands
+! with significant trailing blanks. This length does not consider any
+! possible truncation or padding in assigning the command to the
+! COMMAND argument; in fact the COMMAND argument need not even be
+! present. If the command length cannot be determined, a length of
+! 0 is assigned.
+! STATUS (optional) shall be scalar and of type default integer. It is
+! an INTENT(OUT) argument. It is assigned the value 0 if the
+! command retrieval is sucessful. It is assigned a processor-dependent
+! non-zero value if the command retrieval fails.
+!
+ CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: COMMAND
+ INTEGER , INTENT(OUT), OPTIONAL :: LENGTH
+ INTEGER , INTENT(OUT), OPTIONAL :: STATUS
+!
+ INTEGER :: IARG,NARG,IPOS
+ INTEGER , SAVE :: LENARG
+ CHARACTER(LEN=2000), SAVE :: ARGSTR
+ LOGICAL , SAVE :: GETCMD = .TRUE.
+!
+! Reconstruct the command line from its constituent parts.
+! This may not be the original command line.
+!
+ IF (GETCMD) THEN
+ NARG = IARGC()
+ IF (NARG > 0) THEN
+ IPOS = 1
+ DO IARG = 1,NARG
+ CALL GETARG(IARG,ARGSTR(IPOS:))
+ LENARG = LEN_TRIM(ARGSTR)
+ IPOS = LENARG + 2
+ IF (IPOS > LEN(ARGSTR)) EXIT
+ END DO
+ ELSE
+ ARGSTR = ' '
+ LENARG = 0
+ ENDIF
+ GETCMD = .FALSE.
+ ENDIF
+ IF (PRESENT(COMMAND)) COMMAND = ARGSTR
+ IF (PRESENT(LENGTH)) LENGTH = LENARG
+ IF (PRESENT(STATUS)) STATUS = 0
+ RETURN
+ END SUBROUTINE GET_COMMAND
+!
+ INTEGER FUNCTION COMMAND_ARGUMENT_COUNT()
+!
+! Description. Returns the number of command arguments.
+!
+! Class. Inquiry function
+!
+! Arguments. None.
+!
+! Result Characteristics. Scalar default integer.
+!
+! Result Value. The result value is equal to the number of command
+! arguments available. If there are no command arguments available
+! or if the processor does not support command arguments, then
+! the result value is 0. If the processor has a concept of a command
+! name, the command name does not count as one of the command
+! arguments.
+!
+ COMMAND_ARGUMENT_COUNT = IARGC()
+ RETURN
+ END FUNCTION COMMAND_ARGUMENT_COUNT
+!
+ SUBROUTINE GET_COMMAND_ARGUMENT(NUMBER,VALUE,LENGTH,STATUS)
+!
+! Description. Returns a command argument.
+!
+! Class. Subroutine.
+!
+! Arguments.
+! NUMBER shall be scalar and of type default integer. It is an
+! INTENT(IN) argument. It specifies the number of the command
+! argument that the other arguments give information about. Useful
+! values of NUMBER are those between 0 and the argument count
+! returned by the COMMAND_ARGUMENT_COUNT intrinsic.
+! Other values are allowed, but will result in error status return
+! (see below). Command argument 0 is defined to be the command
+! name by which the program was invoked if the processor has such
+! a concept. It is allowed to call the GET_COMMAND_ARGUMENT
+! procedure for command argument number 0, even if the processor
+! does not define command names or other command arguments.
+! The remaining command arguments are numbered consecutively from
+! 1 to the argument count in an order determined by the processor.
+! VALUE (optional) shall be scalar and of type default character.
+! It is an INTENT(OUT) argument. It is assigned the value of the
+! command argument specified by NUMBER. If the command argument value
+! cannot be determined, VALUE is assigned all blanks.
+! LENGTH (optional) shall be scalar and of type default integer.
+! It is an INTENT(OUT) argument. It is assigned the significant length
+! of the command argument specified by NUMBER. The significant
+! length may include trailing blanks if the processor allows command
+! arguments with significant trailing blanks. This length does not
+! consider any possible truncation or padding in assigning the
+! command argument value to the VALUE argument; in fact the
+! VALUE argument need not even be present. If the command
+! argument length cannot be determined, a length of 0 is assigned.
+! STATUS (optional) shall be scalar and of type default integer.
+! It is an INTENT(OUT) argument. It is assigned the value 0 if
+! the argument retrieval is sucessful. It is assigned a
+! processor-dependent non-zero value if the argument retrieval fails.
+!
+! NOTE
+! One possible reason for failure is that NUMBER is negative or
+! greater than COMMAND_ARGUMENT_COUNT().
+!
+ INTEGER , INTENT(IN) :: NUMBER
+ CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: VALUE
+ INTEGER , INTENT(OUT), OPTIONAL :: LENGTH
+ INTEGER , INTENT(OUT), OPTIONAL :: STATUS
+!
+! A temporary variable for the rare case case where LENGTH is
+! specified but VALUE is not. An arbitrary maximum argument length
+! of 1000 characters should cover virtually all situations.
+!
+ CHARACTER(LEN=1000) :: TMPVAL
+!
+! Possible error codes:
+! 1 = Argument number is less than minimum
+! 2 = Argument number exceeds maximum
+!
+ IF (NUMBER < 0) THEN
+ IF (PRESENT(VALUE )) VALUE = ' '
+ IF (PRESENT(LENGTH)) LENGTH = 0
+ IF (PRESENT(STATUS)) STATUS = 1
+ RETURN
+ ELSE IF (NUMBER > IARGC()) THEN
+ IF (PRESENT(VALUE )) VALUE = ' '
+ IF (PRESENT(LENGTH)) LENGTH = 0
+ IF (PRESENT(STATUS)) STATUS = 2
+ RETURN
+ END IF
+!
+! Get the argument if VALUE is present
+!
+ IF (PRESENT(VALUE)) CALL GETARG(NUMBER,VALUE)
+!
+! As under Unix, the LENGTH option is probably fairly pointless here,
+! but LEN_TRIM is used to ensure at least some sort of meaningful result.
+!
+ IF (PRESENT(LENGTH)) THEN
+ IF (PRESENT(VALUE)) THEN
+ LENGTH = LEN_TRIM(VALUE)
+ ELSE
+ CALL GETARG(NUMBER,TMPVAL)
+ LENGTH = LEN_TRIM(TMPVAL)
+ END IF
+ END IF
+!
+! Since GETARG does not return a result code, assume success
+!
+ IF (PRESENT(STATUS)) STATUS = 0
+ RETURN
+ END SUBROUTINE GET_COMMAND_ARGUMENT
+!
+ END MODULE F2KCLI
Property changes on: DART/trunk/models/wrf/WRF_DART_utilities/f2kcli.f90
___________________________________________________________________
Added: svn:keywords
+ Date Revision Author HeadURL Id
Added: DART/trunk/models/wrf/WRF_DART_utilities/grid_refl_obs.f90
===================================================================
--- DART/trunk/models/wrf/WRF_DART_utilities/grid_refl_obs.f90 (rev 0)
+++ DART/trunk/models/wrf/WRF_DART_utilities/grid_refl_obs.f90 2010-01-28 19:31:54 UTC (rev 4232)
@@ -0,0 +1,400 @@
+PROGRAM grid_refl_obs
+
+! Extract reflectivity obs from a DART obs_seq.out file and compute the coordinates of
+! these observations on a WRF grid (mass grid points).
+!
+! David Dowell 26 June 2007
+!
+! input parameters from command line:
+! (1) obs_seq_file -- path name of DART obs_seq.out file
+! (2) refl_min -- minimum reflectivity threshold for processing observation
+! (3) days_begin -- start of time range to be processed
+! (4) seconds_begin -- " "
+! (5) days_end -- end of time range to be processed
+! (6) seconds_end -- " "
+! (7) wrf_file -- path name of WRF netcdf file
+!
+! output:
+! (1) text file
+
+use types_mod, only : r8, missing_r8, gravity
+use obs_sequence_mod, only : read_obs_seq, obs_type, obs_sequence_type, get_first_obs, &
+ get_obs_from_key, get_obs_def, get_copy_meta_data, &
+ get_obs_time_range, get_time_range_keys, get_num_obs, &
+ get_next_obs, get_num_times, get_obs_values, init_obs, &
+ get_num_copies, static_init_obs_sequence, &
+ get_qc, destroy_obs_sequence, read_obs_seq_header, &
+ get_last_obs, destroy_obs, get_num_qc, get_qc_meta_data
+use obs_def_mod, only : obs_def_type, get_obs_def_error_variance, get_obs_def_time, &
+ get_obs_def_location, get_obs_kind, get_obs_name
+use obs_kind_mod, only : RADAR_REFLECTIVITY
+use map_utils, only : proj_info, map_init, map_set, latlon_to_ij, &
+ PROJ_LATLON, PROJ_MERC, PROJ_LC, PROJ_PS, &
+ ij_to_latlon, gridwind_to_truewind
+use location_mod, only : location_type, get_location, set_location_missing, &
+ write_location, operator(/=), &
+ vert_is_undef, VERTISUNDEF, &
+ vert_is_surface, VERTISSURFACE, &
+ vert_is_level, VERTISLEVEL, &
+ vert_is_pressure, VERTISPRESSURE, &
+ vert_is_height, VERTISHEIGHT
+use time_manager_mod, only : time_type, set_date, set_time, get_time, print_time, &
+ set_calendar_type, print_date, GREGORIAN, &
+ operator(*), operator(+), operator(-), &
+ operator(>), operator(<), operator(/), &
+ operator(/=), operator(<=)
+use utilities_mod, only : error_handler, E_ERR, E_MSG, initialize_utilities, &
+ timestamp, register_module, logfileunit, file_exist
+use netcdf
+use f2kcli
+
+implicit none
+
+! version controlled file description for error handling, do not edit
+character(len=128), parameter :: &
+ source = "$URL$", &
+ revision = "$Revision$", &
+ revdate = "$Date$"
+
+! command-line parameters
+character(len=129) :: obs_seq_file
+real(r8) :: refl_min
+integer :: seconds_begin
+integer :: days_begin
+integer :: seconds_end
+integer :: days_end
+character(len=129) :: wrf_file
+
+! local variables
+type(obs_sequence_type) :: seq
+type(obs_def_type) :: obs_def
+type(time_type) :: beg_time, end_time
+type(obs_type) :: ob
+type(location_type) :: ob_loc
+
+integer :: num_copies, num_qc, num_obs, max_num_obs, obs_seq_file_id
+character(len=129) :: obs_seq_read_format
+logical :: pre_I_format
+logical :: is_there_one, out_of_range
+integer :: key_bounds(2)
+integer :: num_obs_in_time_period
+integer :: num_refl_obs
+integer :: num_refl_obs_in_domain
+integer, allocatable :: keys(:)
+integer :: obs_kind_ind
+integer :: obs_index
+real(r8) :: ob_value(1)
+
+real(r8), allocatable :: lat(:,:) ! latitude at mass grid points (deg)
+real(r8), allocatable :: lon(:,:) ! longitude at mass grid points (deg)
+real(r8), allocatable :: phb(:,:,:) ! base-state geopotential (m^2 s^-2)
+real(r8), allocatable :: ht(:,:,:) ! height MSL of mass grid points (m)
+type(proj_info) :: proj ! map projection info.
+integer, parameter :: map_sphere = 0, map_lambert = 1, map_polar_stereo = 2, map_mercator = 3
+integer map_proj, proj_code
+real(r8) :: dx
+real(r8) :: stdlon,truelat1,truelat2
+real(r8) :: xyz_loc(3)
+real(r8) :: iloc, jloc, kloc
+
+real(r8), allocatable :: refl_ob(:,:,:) ! gridded reflectivity observations (dBZ)
+
+integer :: bt, sn, we ! WRF grid dimensions
+
+character(len = 150) :: msgstring
+integer :: i, j, k, o
+
+
+! netcdf stuff
+integer :: var_id, ncid, ierr
+character(len=80) :: varname
+
+! f2kcli stuff
+integer :: status, length
+character(len=120) :: string
+
+
+
+! Get command-line parameters, using the F2KCLI interface. See f2kcli.f90 for details.
+
+if( COMMAND_ARGUMENT_COUNT() .ne. 7 ) then
+ print*, 'INCORRECT # OF ARGUMENTS ON COMMAND LINE: ', COMMAND_ARGUMENT_COUNT()
+ call exit(1)
+else
+
+ call GET_COMMAND_ARGUMENT(1,obs_seq_file,length,status)
+ if( status .ne. 0 ) then
+ print*, 'obs_seq_file NOT RETRIEVED FROM COMMAND LINE: ', status
+ call exit(1)
+ endif
+
+ call GET_COMMAND_ARGUMENT(2,string,length,status)
+ if( status .ne. 0 ) then
+ print*, 'refl_min NOT RETRIEVED FROM COMMAND LINE: ', status
+ call exit(1)
+ else
+ read(string,*) refl_min
+ endif
+
+ call GET_COMMAND_ARGUMENT(3,string,length,status)
+ if( status .ne. 0 ) then
+ print*, 'days_begin NOT RETRIEVED FROM COMMAND LINE: ', status
+ call exit(1)
+ else
+ read(string,*) days_begin
+ endif
+
+ call GET_COMMAND_ARGUMENT(4,string,length,status)
+ if( status .ne. 0 ) then
+ print*, 'seconds_begin NOT RETRIEVED FROM COMMAND LINE: ', status
+ call exit(1)
+ else
+ read(string,*) seconds_begin
+ endif
+
+ call GET_COMMAND_ARGUMENT(5,string,length,status)
+ if( status .ne. 0 ) then
+ print*, 'days_end NOT RETRIEVED FROM COMMAND LINE: ', status
+ call exit(1)
+ else
+ read(string,*) days_end
+ endif
+
+ call GET_COMMAND_ARGUMENT(6,string,length,status)
+ if( status .ne. 0 ) then
+ print*, 'seconds_end NOT RETRIEVED FROM COMMAND LINE: ', status
+ call exit(1)
+ else
+ read(string,*) seconds_end
+ endif
+
+ call GET_COMMAND_ARGUMENT(7,wrf_file,length,status)
+ if( status .ne. 0 ) then
+ print*, 'wrf_file NOT RETRIEVED FROM COMMAND LINE: ', status
+ call exit(1)
+ endif
+
+endif
+
+
+! Read observations.
+
+call static_init_obs_sequence() ! Initialize the obs sequence module
+
+obs_seq_file = trim(adjustl(obs_seq_file))
+if ( file_exist(obs_seq_file) ) then
+ write(msgstring,*)'opening ', obs_seq_file
+ call error_handler(E_MSG,'grid_refl_obs',msgstring,source,revision,revdate)
+else
+ write(msgstring,*)obs_seq_file,&
+ ' does not exist. Finishing up.'
+ call error_handler(E_MSG,'grid_refl_obs',msgstring,source,revision,revdate)
+ call exit(1)
+endif
+
+call read_obs_seq_header(obs_seq_file, &
+ num_copies, num_qc, num_obs, max_num_obs, &
+ obs_seq_file_id, obs_seq_read_format, pre_I_format, &
+ close_the_file = .true.)
+
+print*, 'num_copies = ', num_copies
+print*, 'num_qc = ', num_qc
+print*, 'num_obs = ', num_obs
+print*, 'max_num_obs = ', max_num_obs
+print*, 'obs_seq_read_format = ', obs_seq_read_format
+print*, 'pre_I_format = ', pre_I_format
+
+call read_obs_seq(obs_seq_file, 0, 0, 0, seq)
+
+MetaDataLoop : do i=1, get_num_copies(seq)
+ if(index(get_copy_meta_data(seq,i), 'observation') > 0) obs_index = i
+enddo MetaDataLoop
+
+
+! Open WRF file and obtain grid dimensions.
+
+call check ( nf90_open(wrf_file, NF90_NOWRITE, ncid) )
+
+call check ( nf90_inq_dimid(ncid, 'bottom_top', var_id) )
+call check ( nf90_inquire_dimension(ncid, var_id, varname, bt) )
+
+call check ( nf90_inq_dimid(ncid, 'south_north', var_id) )
+call check ( nf90_inquire_dimension(ncid, var_id, varname, sn) )
+
+call check ( nf90_inq_dimid(ncid, 'west_east', var_id) )
+call check ( nf90_inquire_dimension(ncid, var_id, varname, we) )
+
+call check( nf90_get_att(ncid, nf90_global, 'MAP_PROJ', map_proj) )
+call check( nf90_get_att(ncid, nf90_global, 'DX', dx) )
+call check( nf90_get_att(ncid, nf90_global, 'TRUELAT1', truelat1) )
+call check( nf90_get_att(ncid, nf90_global, 'TRUELAT2', truelat2) )
+call check( nf90_get_att(ncid, nf90_global, 'STAND_LON', stdlon) )
+
+! Allocate arrays.
+
+allocate(lat(we,sn))
+allocate(lon(we,sn))
+allocate(phb(we,sn,bt+1))
+allocate(ht(we,sn,bt))
+allocate(refl_ob(we,sn,bt))
+
+refl_ob(:,:,:) = missing_r8
+
+
+! Read WRF grid information.
+
+call check ( nf90_inq_varid(ncid, 'XLAT', var_id))
+call check ( nf90_get_var(ncid, var_id, lat, start = (/ 1, 1, 1/)))
+
+call check ( nf90_inq_varid(ncid, 'XLONG', var_id))
+call check ( nf90_get_var(ncid, var_id, lon, start = (/ 1, 1, 1/)))
+
+call check ( nf90_inq_varid(ncid, 'PHB', var_id))
+call check ( nf90_get_var(ncid, var_id, phb, start = (/ 1, 1, 1, 1/)))
+
+ierr = NF90_close(ncid)
+
+
+! Set up map projection structure.
+
+call map_init(proj)
+if(map_proj == map_sphere) then
+ proj_code = PROJ_LATLON
+elseif(map_proj == map_lambert) then
+ proj_code = PROJ_LC
+elseif(map_proj == map_polar_stereo) then
+ proj_code = PROJ_PS
+elseif(map_proj == map_mercator) then
+ proj_code = PROJ_MERC
+else
+ call error_handler(E_ERR,'grid_refl_obs', &
+ 'Map projection no supported.', source, revision, revdate)
+endif
+!call map_set(proj_code,lat(1,1),lon(1,1), &
+! 1.0_r8,1.0_r8,dx,stdlon,truelat1,truelat2,proj)
+call map_set(proj_code=proj_code, proj=proj, lat1=lat(1,1), lon1=lon(1,1), &
+ knowni=1.0_r8, knownj=1.0_r8, dx=dx, stdlon=stdlon, truelat1=truelat1, truelat2=truelat2)
+
+! Compute height (m MSL) of each grid point.
+
+do k=1, bt
+ do j=1, sn
+ do i=1, we
+ ht(i,j,k) = ( phb(i,j,k) + phb(i,j,k+1) ) / (2.0_r8*gravity)
+ enddo
+ enddo
+enddo
+
+! Make sure longitudes are in the range from 0 to 360.
+
+do j=1, sn
+ do i=1, we
+ do while (lon(i,j) < 0.0_r8)
+ lon(i,j) = lon(i,j) + 360.0_r8
+ end do
+ do while (lon(i,j) > 360.0_r8)
+ lon(i,j) = lon(i,j) - 360.0_r8
+ end do
+ enddo
+enddo
+
+
+! Process observations, assigning to grid points any reflectivity observations that lie within
+! the domain and specified time range.
+
+beg_time = set_time(seconds_begin, days_begin)
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list