[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