Dennis,
Thanks for your help. I am close, but not quite there.
I made a couple of changes. I multiplied dp_sigma by
100 to get Pascals since ps was in hPa. I changed
the lat/lons in the "pw = " line to iy/jx/kz to match the data.
Check out the attached plot from ncview. The basic pattern is there
but it is all messy. What might be wrong?
a=addfile("../con9096/Southwest_ATM.1991070100.nc","r")
q=a->qv ; kg per kg
ps=a->ps ; hPa
sigma=a->sigma
A=a->sigma
A=0.0
P0=0.0
dp_sigma = dpres_hybrid_ccm(ps,P0,A,sigma)
dp_sigma=dp_sigma*100.
dp_sigma!0 = ps!0
dp_sigma!1 = "plevel"
dp_sigma!2 = ps!1
dp_sigma!3 = ps!2
pw = prcwater_dp(q(time|:,iy|:,jx|:,kz|:),dp_sigma(time|:,iy|:,jx|:,plevel|:))
pw@long_name = "total column precipitable water"
pw@units = "kg/m2"
copy_VarCoords(ps,pw)
out=addfile("pw.nc","c")
out->pw=pw
Mike
netcdf Southwest_ATM.1991070100 {
dimensions:
iy = 148 ;
jx = 158 ;
time = UNLIMITED ; // (124 currently)
kz = 18 ;
variables:
int rcm_map ;
rcm_map:grid_mapping_name = "lambert_conformal_conic" ;
rcm_map:grid_size_in_meters = 30000.f ;
rcm_map:latitude_of_projection_origin = 30.f ;
rcm_map:longitude_of_projection_origin = -110.f ;
rcm_map:longitude_of_central_meridian = -110.f ;
rcm_map:standard_parallel = 20.f, 40.f ;
float sigma(kz) ;
sigma:standard_name = "atmosphere_sigma_coordinate" ;
sigma:long_name = "Sigma at model layers" ;
sigma:units = "1" ;
sigma:axis = "Z" ;
sigma:positive = "down" ;
sigma:formula_terms = "sigma: sigma ps: ps ptop: ptop" ;
float ptop ;
ptop:standard_name = "air_pressure" ;
ptop:long_name = "Pressure at model top" ;
ptop:units = "hPa" ;
float iy(iy) ;
iy:standard_name = "projection_y_coordinate" ;
iy:long_name = "y-coordinate in Cartesian system" ;
iy:units = "km" ;
float jx(jx) ;
jx:standard_name = "projection_x_coordinate" ;
jx:long_name = "x-coordinate in Cartesian system" ;
jx:units = "km" ;
float xlat(iy, jx) ;
xlat:standard_name = "latitude" ;
xlat:long_name = "Latitude at cross points" ;
xlat:units = "degrees_north" ;
xlat:actual_range = 8.549747f, 49.84976f ;
float xlon(iy, jx) ;
xlon:standard_name = "longitude" ;
xlon:long_name = "Longitude at cross points" ;
xlon:units = "degrees_east" ;
xlon:actual_range = 0.f, 0.f ;
float topo(iy, jx) ;
topo:standard_name = "surface_altitude" ;
topo:long_name = "Domain surface elevation" ;
topo:units = "m" ;
topo:coordinates = "xlat xlon" ;
topo:grid_mapping = "rcm_map" ;
float mask(iy, jx) ;
mask:standard_name = "landmask" ;
mask:long_name = "Domain land/ocean mask" ;
mask:units = "1" ;
mask:coordinates = "xlat xlon" ;
mask:grid_mapping = "rcm_map" ;
double time(time) ;
time:standard_name = "time" ;
time:long_name = "time" ;
time:calendar = "standard" ;
time:units = "hours since 1990-01-01 00:00:00 UTC" ;
float ps(time, iy, jx) ;
ps:standard_name = "surface_air_pressure" ;
ps:long_name = "Surface pressure" ;
ps:units = "hPa" ;
ps:coordinates = "xlat xlon" ;
ps:grid_mapping = "rcm_map" ;
float u(time, kz, iy, jx) ;
u:standard_name = "eastward_wind" ;
u:long_name = "U component (westerly) of wind" ;
u:units = "m s-1" ;
u:coordinates = "xlat xlon" ;
u:grid_mapping = "rcm_map" ;
float v(time, kz, iy, jx) ;
v:standard_name = "northward_wind" ;
v:long_name = "V component (southerly) of wind" ;
v:units = "m s-1" ;
v:coordinates = "xlat xlon" ;
v:grid_mapping = "rcm_map" ;
float omega(time, kz, iy, jx) ;
omega:standard_name = "lagrangian_tendency_of_air_pressure" ;
omega:long_name = "Pressure velocity" ;
omega:units = "hPa s-1" ;
omega:coordinates = "xlat xlon" ;
omega:grid_mapping = "rcm_map" ;
float t(time, kz, iy, jx) ;
t:standard_name = "air_temperature" ;
t:long_name = "Temperature" ;
t:units = "K" ;
t:coordinates = "xlat xlon" ;
t:grid_mapping = "rcm_map" ;
float qv(time, kz, iy, jx) ;
qv:standard_name = "humidity_mixing_ratio" ;
qv:long_name = "Water vapor mixing ratio" ;
qv:units = "kg kg-1" ;
qv:coordinates = "xlat xlon" ;
qv:grid_mapping = "rcm_map" ;
float qc(time, kz, iy, jx) ;
qc:standard_name = "cloud_liquid_water_mixing_ratio" ;
qc:long_name = "Cloud water mixing ratio" ;
qc:units = "kg kg-1" ;
qc:coordinates = "xlat xlon" ;
qc:grid_mapping = "rcm_map" ;
float tpr(time, iy, jx) ;
tpr:standard_name = "precipitation_flux" ;
tpr:long_name = "Total daily precipitation rate" ;
tpr:units = "kg m-2 day-1" ;
tpr:coordinates = "xlat xlon" ;
tpr:grid_mapping = "rcm_map" ;
float tgb(time, iy, jx) ;
tgb:standard_name = "soil_temperature" ;
tgb:long_name = "Lower groud temperature in BATS" ;
tgb:units = "K" ;
tgb:coordinates = "xlat xlon" ;
tgb:grid_mapping = "rcm_map" ;
float swt(time, iy, jx) ;
swt:standard_name = "moisture_content_of_soil_layer" ;
swt:long_name = "Total soil water" ;
swt:units = "kg m-2" ;
swt:coordinates = "xlat xlon" ;
swt:grid_mapping = "rcm_map" ;
swt:_FillValue = -1.e+34f ;
float rno(time, iy, jx) ;
rno:standard_name = "runoff_flux" ;
rno:long_name = "Runoff accumulated infiltration" ;
rno:units = "kg m-2 day-1" ;
rno:coordinates = "xlat xlon" ;
rno:grid_mapping = "rcm_map" ;
rno:_FillValue = -1.e+34f ;
// global attributes:
:title = "ICTP Regional Climatic model V4 ATM output" ;
:institution = "ICTP" ;
:source = "RegCM Model SVN_REV simulation output" ;
:Conventions = "CF-1.4" ;
:history = "2011-02-12 20:02:51 : Created by RegCM model" ;
:references = "http://eforge.escience-lab.org/gf/project/regcm" ;
:experiment = "Southwest" ;
:projection = "LAMCON" ;
:grid_mapping_name = "lambert_conformal_conic" ;
:grid_size_in_meters = 30000.f ;
:latitude_of_projection_origin = 30.f ;
:longitude_of_projection_origin = -110.f ;
:longitude_of_central_meridian = -110.f ;
:standard_parallel = 20.f, 40.f ;
:model_IPCC_scenario = "A1B" ;
:model_boundary_conditions = "Relaxation, linear technique" ;
:model_cumulous_convection_scheme = "Grell" ;
:model_convective_closure_assumption = "Fritsch & Chappell (1980)" ;
:model_boundary_layer_scheme = "Holtslag PBL (Holtslag, 1990)" ;
:model_moist_physics_scheme = "Explicit moisture (SUBEX; Pal et al 2000)" ;
:model_ocean_flux_scheme = "Zeng et al (1998)" ;
:model_pressure_gradient_force_scheme = "Use full fields" ;
:model_use_emission_factor = "No" ;
:model_use_lake_model = "No" ;
:model_chemistry = "Not active" ;
:model_diurnal_cycle_sst = "Not active" ;
:model_seaice_effect = "Not active" ;
:model_seasonal_desert_albedo_effect = "Not active" ;
:model_simulation_initial_start = 1990010100 ;
:model_simulation_start = 1991030100 ;
:model_simulation_expected_end = 1997010100 ;
:model_simulation_is_a_restart = "Yes" ;
:model_timestep_in_seconds = 200. ;
:model_timestep_in_minutes_solar_rad_calc = 30. ;
:model_timestep_in_seconds_bats_calc = 300. ;
:model_timestep_in_hours_radiation_calc = 18. ;
:model_timestep_in_hours_boundary_input = 6 ;
}
On Feb 22, 2011, at 6:24 PM, Dennis Shea wrote:
> One approach
>
> [1]
> http://www.ncl.ucar.edu/Document/Functions/Built-in/dpres_hybrid_ccm.shtml
>
> This uses:
> p(k) = A(k)*PO + B(k)*PS
> The B(k) are the the sigma levels.
> A = sigma
> A = 0.0
> P0 = 0.0
> dp_sigma = dpres_hybrid_ccm(ps,P0,A,sigma)
> dp_sigma!0 = ps!0 ; (say) time
> dp_sigma!1 = "plevel"
> dp_sigma!2 = ps!1 ; lat
> dp_sigma!3 = ps!2 ; lon
> printVarSummary(dp_sigma)
>
>
> [2]
> http://www.ncl.ucar.edu/Document/Functions/Built-in/prcwater_dp.shtml
> pw = prcwater_dp(q(time|:,lat|:,lon|:,sigma|:),dp_sigma(time|:,lat|:,lon|:,plevel|:))
> pw@long_name = "total column precipitable water"
> pw@units = "kg/m2"
>
>
>
>
> On 02/22/2011 02:36 PM, Michael Notaro wrote:
>> How can I compute precipitable water
>> on a lat2d x lon2d grid using sigma
>> coordinates, given water vapor mixing
>> ratio output and surface pressure?
>>
>> Thanks, Michael
>> _______________________________________________
>> ncl-talk mailing list
>> List instructions, subscriber options, unsubscribe:
>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
> --
> ======================================================
> Dennis J. Shea tel: 303-497-1361 |
> P.O. Box 3000 fax: 303-497-1333 |
> Climate Analysis Section |
> Climate & Global Dynamics Div. |
> National Center for Atmospheric Research |
> Boulder, CO 80307 |
> USA email: shea 'at' ucar.edu |
> ======================================================
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
This archive was generated by hypermail 2.1.8 : Wed Feb 23 2011 - 16:47:57 MST