[Dart-dev] [5475] DART/branches/development:
nancy at ucar.edu
nancy at ucar.edu
Wed Jan 4 15:48:34 MST 2012
Revision: 5475
Author: nancy
Date: 2012-01-04 15:48:33 -0700 (Wed, 04 Jan 2012)
Log Message:
-----------
add 7 new kinds for the advanced microphysics in wrf.
added code from glen romine for interpolating observation
values using the new microphysics vars. this version also
has the option to interpolate pressure values in a log(p)
scale instead of linear in pressure. the default is to only
use the log scale for the vertical interpolations but there
are options (see the module global vars) for doing log interp
in the horizontal as well.
Modified Paths:
--------------
DART/branches/development/models/wrf/model_mod.f90
DART/branches/development/obs_kind/DEFAULT_obs_kind_mod.F90
-------------- next part --------------
Modified: DART/branches/development/models/wrf/model_mod.f90
===================================================================
--- DART/branches/development/models/wrf/model_mod.f90 2011-12-31 00:03:46 UTC (rev 5474)
+++ DART/branches/development/models/wrf/model_mod.f90 2012-01-04 22:48:33 UTC (rev 5475)
@@ -52,7 +52,8 @@
KIND_SURFACE_PRESSURE, KIND_TEMPERATURE, &
KIND_SPECIFIC_HUMIDITY, KIND_SURFACE_ELEVATION, &
KIND_PRESSURE, KIND_VERTICAL_VELOCITY, &
- KIND_RAINWATER_MIXING_RATIO, KIND_DENSITY, &
+ KIND_DENSITY, KIND_FLASH_RATE_2D, &
+ KIND_RAINWATER_MIXING_RATIO, KIND_HAIL_MIXING_RATIO, &
KIND_GRAUPEL_MIXING_RATIO, KIND_SNOW_MIXING_RATIO, &
KIND_CLOUD_LIQUID_WATER, KIND_CLOUD_ICE, &
KIND_CONDENSATIONAL_HEATING, KIND_VAPOR_MIXING_RATIO, &
@@ -60,8 +61,11 @@
KIND_POTENTIAL_TEMPERATURE, KIND_SOIL_MOISTURE, &
KIND_DROPLET_NUMBER_CONCENTR, KIND_SNOW_NUMBER_CONCENTR, &
KIND_RAIN_NUMBER_CONCENTR, KIND_GRAUPEL_NUMBER_CONCENTR, &
+ KIND_HAIL_NUMBER_CONCENTR, KIND_HAIL_VOLUME, &
+ KIND_GRAUPEL_VOLUME, KIND_DIFFERENTIAL_REFLECTIVITY, &
+ KIND_RADAR_REFLECTIVITY, KIND_POWER_WEIGHTED_FALL_SPEED, &
+ KIND_SPECIFIC_DIFFERENTIAL_PHASE, &
KIND_VORTEX_LAT, KIND_VORTEX_LON, &
- KIND_RADAR_REFLECTIVITY, KIND_POWER_WEIGHTED_FALL_SPEED,&
KIND_VORTEX_PMIN, KIND_VORTEX_WMAX, &
KIND_SKIN_TEMPERATURE, KIND_LANDMASK, &
get_raw_obs_kind_index, get_num_raw_obs_kinds, &
@@ -182,6 +186,11 @@
integer :: vert_localization_coord = VERTISHEIGHT
! Allow observations above the surface but below the lowest sigma level.
logical :: allow_obs_below_vol = .false.
+! Do the interpolation of pressure values only after taking the log (.true.)
+! vs doing a linear interpolation directly in pressure units (.false.)
+logical :: log_vert_interp = .true.
+logical :: log_horz_interpM = .false.
+logical :: log_horz_interpQ = .false.
!nc -- we are adding these to the model.nml until they appear in the NetCDF files
logical :: polar = .false. ! wrap over the poles
logical :: periodic_x = .false. ! wrap in longitude or x
@@ -276,9 +285,12 @@
! JPH local variables to hold type indices
integer :: type_u, type_v, type_w, type_t, type_qv, type_qr, type_hdiab, &
type_qndrp, type_qnsnow, type_qnrain, type_qngraupel, type_qnice, &
- type_qc, type_qg, type_qi, type_qs, type_gz, type_refl, type_fall_spd
+ type_qc, type_qg, type_qi, type_qs, type_gz, type_refl, type_fall_spd, &
+ type_dref, type_spdp, type_qh, type_qnhail, type_qhvol, type_qgvol
+
integer :: type_u10, type_v10, type_t2, type_th2, type_q2, &
- type_ps, type_mu, type_tsk, type_tslb, type_sh2o, type_smois
+ type_ps, type_mu, type_tsk, type_tslb, type_sh2o, &
+ type_smois, type_2dflash
integer :: number_of_wrf_variables
integer, dimension(:,:), pointer :: var_index
@@ -642,13 +654,17 @@
wrf%dom(id)%type_qr = get_type_ind_from_type_string(id,'QRAIN')
wrf%dom(id)%type_qc = get_type_ind_from_type_string(id,'QCLOUD')
wrf%dom(id)%type_qg = get_type_ind_from_type_string(id,'QGRAUP')
+ wrf%dom(id)%type_qh = get_type_ind_from_type_string(id,'QHAIL')
wrf%dom(id)%type_qi = get_type_ind_from_type_string(id,'QICE')
wrf%dom(id)%type_qs = get_type_ind_from_type_string(id,'QSNOW')
+ wrf%dom(id)%type_qgvol = get_type_ind_from_type_string(id,'QVGRAUPEL')
+ wrf%dom(id)%type_qhvol = get_type_ind_from_type_string(id,'QVHAIL')
wrf%dom(id)%type_qnice = get_type_ind_from_type_string(id,'QNICE')
wrf%dom(id)%type_qndrp = get_type_ind_from_type_string(id,'QNDRP')
wrf%dom(id)%type_qnsnow = get_type_ind_from_type_string(id,'QNSNOW')
wrf%dom(id)%type_qnrain = get_type_ind_from_type_string(id,'QNRAIN')
wrf%dom(id)%type_qngraupel = get_type_ind_from_type_string(id,'QNGRAUPEL')
+ wrf%dom(id)%type_qnhail = get_type_ind_from_type_string(id,'QNHAIL')
wrf%dom(id)%type_u10 = get_type_ind_from_type_string(id,'U10')
wrf%dom(id)%type_v10 = get_type_ind_from_type_string(id,'V10')
wrf%dom(id)%type_t2 = get_type_ind_from_type_string(id,'T2')
@@ -657,10 +673,13 @@
wrf%dom(id)%type_ps = get_type_ind_from_type_string(id,'PSFC')
wrf%dom(id)%type_mu = get_type_ind_from_type_string(id,'MU')
wrf%dom(id)%type_tsk = get_type_ind_from_type_string(id,'TSK')
+ wrf%dom(id)%type_2dflash = get_type_ind_from_type_string(id,'FLASH_RATE_2D')
wrf%dom(id)%type_tslb = get_type_ind_from_type_string(id,'TSLB')
wrf%dom(id)%type_smois = get_type_ind_from_type_string(id,'SMOIS')
wrf%dom(id)%type_sh2o = get_type_ind_from_type_string(id,'SH2O')
wrf%dom(id)%type_refl = get_type_ind_from_type_string(id,'REFL_10CM')
+ wrf%dom(id)%type_dref = get_type_ind_from_type_string(id,'DIFF_REFL_10CM')
+ wrf%dom(id)%type_spdp = get_type_ind_from_type_string(id,'SPEC_DIFF_10CM')
wrf%dom(id)%type_fall_spd = get_type_ind_from_type_string(id,'FALL_SPD_Z_WEIGHTED')
wrf%dom(id)%type_hdiab = get_type_ind_from_type_string(id,'H_DIABATIC')
@@ -1216,7 +1235,8 @@
! 1.f Specific Humidity (SH, SH2)
! 1.g Vapor Mixing Ratio (QV, Q2)
! 1.h Rainwater Mixing Ratio (QR)
- ! 1.i Graupel Mixing Ratio (QG)
+ ! 1.i.1 Graupel Mixing Ratio (QG)
+ ! 1.i.2 Hail Mixing Ratio (QH)
! 1.j Snow Mixing Ratio (QS)
! 1.k Ice Mixing Ratio (QI)
! 1.l Cloud Mixing Ratio (QC)
@@ -1224,12 +1244,15 @@
! 1.n Ice Number Concentration (QNICE)
! 1.o Snow Number Concentration (QNSNOW)
! 1.p Rain Number Concentration (QNRAIN)
- ! 1.q Graupel Number Concentration (QNGRAUPEL)
+ ! 1.q.1 Graupel Number Concentration (QNGRAUPEL)
+ ! 1.q.2 Hail Number Concentration (QNHAIL)
! 1.r Previous time step condensational heating (H_DIABATIC)
! 1.s Reflectivity weighted precip fall speed (FALL_SPD_Z_WEIGHTED)
! 1.t Pressure (P)
! 1.u Vortex Center Stuff from Yongsheng
- ! 1.v Radar Reflectivity (REFL_10CM)
+ ! 1.v.1 Radar Reflectivity (REFL_10CM)
+ ! 1.v.2 Differential Reflectivity (DIFF_REFL_10CM)
+ ! 1.v.3 Specific Differential Phase (SPEC_DIFF_10CM)
! 1.w Geopotential Height (GZ)
! 1.x Surface Elevation (HGT)
! 1.y Surface Skin Temperature (TSK)
@@ -1732,6 +1755,9 @@
fld(2) = dym*( dxm*x(ill) + dx*x(ilr) ) + dy*( dxm*x(iul) + dx*x(iur) )
+ ! Don't accept negative water vapor amounts (?)
+ fld = max(0.0_r8, fld)
+
endif
endif
@@ -1758,6 +1784,9 @@
fld(1) = dym*( dxm*x(ill) + dx*x(ilr) ) + dy*( dxm*x(iul) + dx*x(iur) )
+ ! Don't accept negative water vapor amounts (?)
+ fld = max(0.0_r8, fld)
+
endif
endif
endif
@@ -1803,7 +1832,7 @@
!-----------------------------------------------------
- ! 1.i Graupel Mixing Ratio (QG)
+ ! 1.i.1 Graupel Mixing Ratio (QG)
else if( obs_kind == KIND_GRAUPEL_MIXING_RATIO ) then
! Confirm that QG is in the DART state vector
@@ -1834,14 +1863,51 @@
fld(2) = dym*( dxm*x(ill) + dx*x(ilr) ) + dy*( dxm*x(iul) + dx*x(iur) )
- ! Don't accept negative rain amounts (?)
+ ! Don't accept negative graupel amounts (?)
fld = max(0.0_r8, fld)
endif
endif
-
!-----------------------------------------------------
+ ! 1.i.2 Hail Mixing Ratio (QH)
+ else if( obs_kind == KIND_HAIL_MIXING_RATIO ) then
+
+ ! Confirm that QH is in the DART state vector
+ if ( wrf%dom(id)%type_qh >= 0 ) then
+
+ ! Check to make sure retrieved integer gridpoints are in valid range
+ if ( boundsCheck( i, wrf%dom(id)%periodic_x, id, dim=1, type=wrf%dom(id)%type_t ) .and. &
+ boundsCheck( j, wrf%dom(id)%polar, id, dim=2, type=wrf%dom(id)%type_t ) .and. &
+ boundsCheck( k, .false., id, dim=3, type=wrf%dom(id)%type_t ) ) then
+
+ call getCorners(i, j, id, wrf%dom(id)%type_t, ll, ul, lr, ur, rc )
+ if ( rc .ne. 0 ) &
+ print*, 'model_mod.f90 :: model_interpolate :: getCorners QH rc = ', rc
+
+ ! Interpolation for QH field at level k
+ ill = wrf%dom(id)%dart_ind(ll(1), ll(2), k, wrf%dom(id)%type_qh)
+ iul = wrf%dom(id)%dart_ind(ul(1), ul(2), k, wrf%dom(id)%type_qh)
+ ilr = wrf%dom(id)%dart_ind(lr(1), lr(2), k, wrf%dom(id)%type_qh)
+ iur = wrf%dom(id)%dart_ind(ur(1), ur(2), k, wrf%dom(id)%type_qh)
+
+ fld(1) = dym*( dxm*x(ill) + dx*x(ilr) ) + dy*( dxm*x(iul) + dx*x(iur) )
+
+ ! Interpolation for QH field at level k+1
+ ill = wrf%dom(id)%dart_ind(ll(1), ll(2), k+1, wrf%dom(id)%type_qh)
+ iul = wrf%dom(id)%dart_ind(ul(1), ul(2), k+1, wrf%dom(id)%type_qh)
+ ilr = wrf%dom(id)%dart_ind(lr(1), lr(2), k+1, wrf%dom(id)%type_qh)
+ iur = wrf%dom(id)%dart_ind(ur(1), ur(2), k+1, wrf%dom(id)%type_qh)
+
+ fld(2) = dym*( dxm*x(ill) + dx*x(ilr) ) + dy*( dxm*x(iul) + dx*x(iur) )
+
+ ! Don't accept negative hail amounts (?)
+ fld = max(0.0_r8, fld)
+
+ endif
+ endif
+
+ !-----------------------------------------------------
! 1.j Snow Mixing Ratio (QS)
else if( obs_kind == KIND_SNOW_MIXING_RATIO ) then
@@ -1873,7 +1939,7 @@
fld(2) = dym*( dxm*x(ill) + dx*x(ilr) ) + dy*( dxm*x(iul) + dx*x(iur) )
- ! Don't accept negative rain amounts (?)
+ ! Don't accept negative snow amounts (?)
fld = max(0.0_r8, fld)
endif
@@ -1992,7 +2058,7 @@
fld(2) = dym*( dxm*x(ill) + dx*x(ilr) ) + dy*( dxm*x(iul) + dx*x(iur) )
- ! Don't accept negative droplet concentrations (?)
+ ! Don't accept negative droplet number concentrations (?)
fld = max(0.0_r8, fld)
endif
@@ -2031,7 +2097,7 @@
fld(2) = dym*( dxm*x(ill) + dx*x(ilr) ) + dy*( dxm*x(iul) + dx*x(iur) )
- ! Don't accept negative ice concentrations (?)
+ ! Don't accept negative ice number concentrations (?)
fld = max(0.0_r8, fld)
endif
@@ -2070,7 +2136,7 @@
fld(2) = dym*( dxm*x(ill) + dx*x(ilr) ) + dy*( dxm*x(iul) + dx*x(iur) )
- ! Don't accept negative snow concentrations (?)
+ ! Don't accept negative snow number concentrations (?)
fld = max(0.0_r8, fld)
endif
@@ -2109,7 +2175,7 @@
fld(2) = dym*( dxm*x(ill) + dx*x(ilr) ) + dy*( dxm*x(iul) + dx*x(iur) )
- ! Don't accept negative rain concentrations (?)
+ ! Don't accept negative rain number concentrations (?)
fld = max(0.0_r8, fld)
endif
@@ -2117,7 +2183,7 @@
!-----------------------------------------------------
- ! 1.q Graupel Number Concentration (QNGRAUPEL)
+ ! 1.q.1 Graupel Number Concentration (QNGRAUPEL)
else if( obs_kind == KIND_GRAUPEL_NUMBER_CONCENTR ) then
! Confirm that QNGRAUPEL is in the DART state vector
@@ -2148,13 +2214,50 @@
fld(2) = dym*( dxm*x(ill) + dx*x(ilr) ) + dy*( dxm*x(iul) + dx*x(iur) )
- ! Don't accept negative graupel concentrations (?)
+ ! Don't accept negative graupel number concentrations (?)
fld = max(0.0_r8, fld)
endif
endif
-
+ ! 1.q.2 Hail Number Concentration (QNHAIL)
+ else if( obs_kind == KIND_HAIL_NUMBER_CONCENTR ) then
+
+ ! Confirm that QNHAIL is in the DART state vector
+ if ( wrf%dom(id)%type_qnhail >= 0 ) then
+
+ ! Check to make sure retrieved integer gridpoints are in valid range
+ if ( boundsCheck( i, wrf%dom(id)%periodic_x, id, dim=1, type=wrf%dom(id)%type_t ) .and. &
+ boundsCheck( j, wrf%dom(id)%polar, id, dim=2, type=wrf%dom(id)%type_t ) .and. &
+ boundsCheck( k, .false., id, dim=3, type=wrf%dom(id)%type_t ) ) then
+
+ call getCorners(i, j, id, wrf%dom(id)%type_t, ll, ul, lr, ur, rc )
+ if ( rc .ne. 0 ) &
+ print*, 'model_mod.f90 :: model_interpolate :: getCorners QNHAIL rc = ', rc
+
+ ! Interpolation for QNHAIL field at level k
+ ill = wrf%dom(id)%dart_ind(ll(1), ll(2), k, wrf%dom(id)%type_qnhail)
+ iul = wrf%dom(id)%dart_ind(ul(1), ul(2), k, wrf%dom(id)%type_qnhail)
+ ilr = wrf%dom(id)%dart_ind(lr(1), lr(2), k, wrf%dom(id)%type_qnhail)
+ iur = wrf%dom(id)%dart_ind(ur(1), ur(2), k, wrf%dom(id)%type_qnhail)
+
+ fld(1) = dym*( dxm*x(ill) + dx*x(ilr) ) + dy*( dxm*x(iul) + dx*x(iur) )
+
+ ! Interpolation for QNHAIL field at level k+1
+ ill = wrf%dom(id)%dart_ind(ll(1), ll(2), k+1, wrf%dom(id)%type_qnhail)
+ iul = wrf%dom(id)%dart_ind(ul(1), ul(2), k+1, wrf%dom(id)%type_qnhail)
+ ilr = wrf%dom(id)%dart_ind(lr(1), lr(2), k+1, wrf%dom(id)%type_qnhail)
+ iur = wrf%dom(id)%dart_ind(ur(1), ur(2), k+1, wrf%dom(id)%type_qnhail)
+
+ fld(2) = dym*( dxm*x(ill) + dx*x(ilr) ) + dy*( dxm*x(iul) + dx*x(iur) )
+
+ ! Don't accept negative hail number concentrations (?)
+ fld = max(0.0_r8, fld)
+
+ endif
+ endif
+
+
!-----------------------------------------------------
! 1.r Previous time step condensational heating (H_DIABATIC)
else if( obs_kind == KIND_CONDENSATIONAL_HEATING ) then
@@ -2712,7 +2815,7 @@
!*****************************************************************************
!-----------------------------------------------------
- ! 1.v Radar Reflectivity (REFL_10CM)
+ ! 1.v.1 Radar Reflectivity (REFL_10CM)
else if( obs_kind == KIND_RADAR_REFLECTIVITY ) then
! Confirm that REFL is in the DART state vector
@@ -2745,8 +2848,78 @@
endif
endif
-
+
!-----------------------------------------------------
+ ! 1.v.2 Differential Reflectivity (DIFF_REFL_10CM)
+ else if( obs_kind == KIND_DIFFERENTIAL_REFLECTIVITY ) then
+
+ ! Confirm that DREF is in the DART state vector
+ if ( wrf%dom(id)%type_dref >= 0 ) then
+
+ ! Check to make sure retrieved integer gridpoints are in valid range
+ if ( boundsCheck( i, wrf%dom(id)%periodic_x, id, dim=1, type=wrf%dom(id)%type_t ) .and. &
+ boundsCheck( j, wrf%dom(id)%polar, id, dim=2, type=wrf%dom(id)%type_t ) .and. &
+ boundsCheck( k, .false., id, dim=3, type=wrf%dom(id)%type_t ) ) then
+
+ call getCorners(i, j, id, wrf%dom(id)%type_t, ll, ul, lr, ur, rc )
+ if ( rc .ne. 0 ) &
+ print*, 'model_mod.f90 :: model_interpolate :: getCorners DREF rc = ', rc
+
+ ! Interpolation for DREF field at level k
+ ill = wrf%dom(id)%dart_ind(ll(1), ll(2), k, wrf%dom(id)%type_dref)
+ iul = wrf%dom(id)%dart_ind(ul(1), ul(2), k, wrf%dom(id)%type_dref)
+ ilr = wrf%dom(id)%dart_ind(lr(1), lr(2), k, wrf%dom(id)%type_dref)
+ iur = wrf%dom(id)%dart_ind(ur(1), ur(2), k, wrf%dom(id)%type_dref)
+
+ fld(1) = dym*( dxm*x(ill) + dx*x(ilr) ) + dy*( dxm*x(iul) + dx*x(iur) )
+
+ ! Interpolation for DREF field at level k+1
+ ill = wrf%dom(id)%dart_ind(ll(1), ll(2), k+1, wrf%dom(id)%type_dref)
+ iul = wrf%dom(id)%dart_ind(ul(1), ul(2), k+1, wrf%dom(id)%type_dref)
+ ilr = wrf%dom(id)%dart_ind(lr(1), lr(2), k+1, wrf%dom(id)%type_dref)
+ iur = wrf%dom(id)%dart_ind(ur(1), ur(2), k+1, wrf%dom(id)%type_dref)
+
+ fld(2) = dym*( dxm*x(ill) + dx*x(ilr) ) + dy*( dxm*x(iul) + dx*x(iur) )
+
+ endif
+ endif
+
+ !-----------------------------------------------------
+ ! 1.v.3 Specific Differential Phase (SPEC_DIFF_10CM)
+ else if( obs_kind == KIND_SPECIFIC_DIFFERENTIAL_PHASE ) then
+
+ ! Confirm that SPDP is in the DART state vector
+ if ( wrf%dom(id)%type_spdp >= 0 ) then
+
+ ! Check to make sure retrieved integer gridpoints are in valid range
+ if ( boundsCheck( i, wrf%dom(id)%periodic_x, id, dim=1, type=wrf%dom(id)%type_t ) .and. &
+ boundsCheck( j, wrf%dom(id)%polar, id, dim=2, type=wrf%dom(id)%type_t ) .and. &
+ boundsCheck( k, .false., id, dim=3, type=wrf%dom(id)%type_t ) ) then
+
+ call getCorners(i, j, id, wrf%dom(id)%type_t, ll, ul, lr, ur, rc )
+ if ( rc .ne. 0 ) &
+ print*, 'model_mod.f90 :: model_interpolate :: getCorners SPDP rc = ', rc
+
+ ! Interpolation for SPDP field at level k
+ ill = wrf%dom(id)%dart_ind(ll(1), ll(2), k, wrf%dom(id)%type_spdp)
+ iul = wrf%dom(id)%dart_ind(ul(1), ul(2), k, wrf%dom(id)%type_spdp)
+ ilr = wrf%dom(id)%dart_ind(lr(1), lr(2), k, wrf%dom(id)%type_spdp)
+ iur = wrf%dom(id)%dart_ind(ur(1), ur(2), k, wrf%dom(id)%type_spdp)
+
+ fld(1) = dym*( dxm*x(ill) + dx*x(ilr) ) + dy*( dxm*x(iul) + dx*x(iur) )
+
+ ! Interpolation for SPDP field at level k+1
+ ill = wrf%dom(id)%dart_ind(ll(1), ll(2), k+1, wrf%dom(id)%type_spdp)
+ iul = wrf%dom(id)%dart_ind(ul(1), ul(2), k+1, wrf%dom(id)%type_spdp)
+ ilr = wrf%dom(id)%dart_ind(lr(1), lr(2), k+1, wrf%dom(id)%type_spdp)
+ iur = wrf%dom(id)%dart_ind(ur(1), ur(2), k+1, wrf%dom(id)%type_spdp)
+
+ fld(2) = dym*( dxm*x(ill) + dx*x(ilr) ) + dy*( dxm*x(iul) + dx*x(iur) )
+
+ endif
+ endif
+
+ !-----------------------------------------------------
! 1.w Geopotential Height (GZ)
! GZ is on the ZNW grid (bottom_top_stagger), so its bottom-most level is defined to
@@ -4873,7 +5046,11 @@
! sigma value but set lev0 true
if(pres <= mdl_v(0) .and. pres > mdl_v(1)) then
lev0 = .true.
+ if (log_vert_interp) then
+ zk = (log(mdl_v(0)) - log(pres))/(log(mdl_v(0)) - log(mdl_v(1)))
+ else
zk = (mdl_v(0) - pres)/(mdl_v(0) - mdl_v(1))
+ endif
return
endif
@@ -4881,7 +5058,11 @@
! as a real number, including the fraction between the levels.
do k = 1,n3-1
if(pres <= mdl_v(k) .and. pres >= mdl_v(k+1)) then
+ if (log_vert_interp) then
+ zk = real(k) + (log(mdl_v(k)) - log(pres))/(log(mdl_v(k)) - log(mdl_v(k+1)))
+ else
zk = real(k) + (mdl_v(k) - pres)/(mdl_v(k) - mdl_v(k+1))
+ endif
exit
endif
enddo
@@ -4957,7 +5138,8 @@
pres2 = model_pressure_t(lr(1), lr(2), k,id,x)
pres3 = model_pressure_t(ul(1), ul(2), k,id,x)
pres4 = model_pressure_t(ur(1), ur(2), k,id,x)
- v_p(k) = dym*( dxm*pres1 + dx*pres2 ) + dy*( dxm*pres3 + dx*pres4 )
+
+ v_p(k) = interp_4pressure(pres1, pres2, pres3, pres4, dx, dxm, dy, dym)
enddo
if (debug) &
@@ -4966,15 +5148,15 @@
if ( wrf%dom(id)%type_ps >= 0 ) then
ill = wrf%dom(id)%dart_ind(ll(1), ll(2), 1, wrf%dom(id)%type_ps)
+ ilr = wrf%dom(id)%dart_ind(lr(1), lr(2), 1, wrf%dom(id)%type_ps)
iul = wrf%dom(id)%dart_ind(ul(1), ul(2), 1, wrf%dom(id)%type_ps)
- ilr = wrf%dom(id)%dart_ind(lr(1), lr(2), 1, wrf%dom(id)%type_ps)
iur = wrf%dom(id)%dart_ind(ur(1), ur(2), 1, wrf%dom(id)%type_ps)
! I'm not quite sure where this comes from, but I will trust them on it....
if ( x(ill) /= 0.0_r8 .and. x(ilr) /= 0.0_r8 .and. x(iul) /= 0.0_r8 .and. &
x(iur) /= 0.0_r8 ) then
- v_p(0) = dym*( dxm*x(ill) + dx*x(ilr) ) + dy*( dxm*x(iul) + dx*x(iur) )
+ v_p(0) = interp_4pressure(x(ill), x(ilr), x(iul), x(iur), dx, dxm, dy, dym)
else
@@ -4983,8 +5165,8 @@
pres3 = model_pressure_t(ul(1), ul(2), 2,id,x)
pres4 = model_pressure_t(ur(1), ur(2), 2,id,x)
- v_p(0) = (3.0_r8*v_p(1) - &
- dym*( dxm*pres1 + dx*pres2 ) - dy*( dxm*pres3 + dx*pres4 ))/2.0_r8
+ v_p(0) = interp_4pressure(pres1, pres2, pres3, pres4, dx, dxm, dy, dym, &
+ extrapolate=.true., edgep=v_p(1))
endif
@@ -4995,8 +5177,8 @@
pres3 = model_pressure_t(ul(1), ul(2), 2,id,x)
pres4 = model_pressure_t(ur(1), ur(2), 2,id,x)
- v_p(0) = (3.0_r8*v_p(1) - &
- dym*( dxm*pres1 + dx*pres2 ) - dy*( dxm*pres3 + dx*pres4 ))/2.0_r8
+ v_p(0) = interp_4pressure(pres1, pres2, pres3, pres4, dx, dxm, dy, dym, &
+ extrapolate=.true., edgep=v_p(1))
endif
@@ -5034,19 +5216,19 @@
pres1 = model_pressure_t(i,j,k, id,x)
pres2 = model_pressure_t(i,j,k+1,id,x)
- model_pressure = (3.0_r8*pres1 - pres2)/2.0_r8
+ model_pressure = interp_pressure(pres1, pres2, extrapolate=.true.)
elseif( k == wrf%dom(id)%var_size(3,wrf%dom(id)%type_w) ) then
pres1 = model_pressure_t(i,j,k-1,id,x)
pres2 = model_pressure_t(i,j,k-2,id,x)
- model_pressure = (3.0_r8*pres1 - pres2)/2.0_r8
+ model_pressure = interp_pressure(pres1, pres2, extrapolate=.true.)
else
pres1 = model_pressure_t(i,j,k, id,x)
pres2 = model_pressure_t(i,j,k-1,id,x)
- model_pressure = (pres1 + pres2)/2.0_r8
+ model_pressure = interp_pressure(pres1, pres2)
endif
@@ -5062,14 +5244,14 @@
! We are at seam in longitude, take first and last M-grid points
pres1 = model_pressure_t(i-1,j,k,id,x)
pres2 = model_pressure_t(1, j,k,id,x)
- model_pressure = (pres1 + pres2)/2.0_r8
+ model_pressure = interp_pressure(pres1, pres2, vertical=.false.)
else
! If not periodic, then try extrapolating
pres1 = model_pressure_t(i-1,j,k,id,x)
pres2 = model_pressure_t(i-2,j,k,id,x)
- model_pressure = (3.0_r8*pres1 - pres2)/2.0_r8
+ model_pressure = interp_pressure(pres1, pres2, extrapolate=.true., vertical=.false.)
endif
@@ -5081,14 +5263,14 @@
! We are at seam in longitude, take first and last M-grid points
pres1 = model_pressure_t(i, j,k,id,x)
pres2 = model_pressure_t(wrf%dom(id)%we,j,k,id,x)
- model_pressure = (pres1 + pres2)/2.0_r8
+ model_pressure = interp_pressure(pres1, pres2, vertical=.false.)
else
! If not periodic, then try extrapolating
pres1 = model_pressure_t(i, j,k,id,x)
pres2 = model_pressure_t(i+1,j,k,id,x)
- model_pressure = (3.0_r8*pres1 - pres2)/2.0_r8
+ model_pressure = interp_pressure(pres1, pres2, extrapolate=.true., vertical=.false.)
endif
@@ -5096,7 +5278,7 @@
pres1 = model_pressure_t(i, j,k,id,x)
pres2 = model_pressure_t(i-1,j,k,id,x)
- model_pressure = (pres1 + pres2)/2.0_r8
+ model_pressure = interp_pressure(pres1, pres2, vertical=.false.)
endif
@@ -5115,14 +5297,14 @@
pres1 = model_pressure_t(off,j-1,k,id,x)
pres2 = model_pressure_t(i ,j-1,k,id,x)
- model_pressure = (pres1 + pres2)/2.0_r8
+ model_pressure = interp_pressure(pres1, pres2, vertical=.false.)
! If not periodic, then try extrapolating
else
pres1 = model_pressure_t(i,j-1,k,id,x)
pres2 = model_pressure_t(i,j-2,k,id,x)
- model_pressure = (3.0_r8*pres1 - pres2)/2.0_r8
+ model_pressure = interp_pressure(pres1, pres2, extrapolate=.true., vertical=.false.)
endif
@@ -5137,14 +5319,14 @@
pres1 = model_pressure_t(off,j,k,id,x)
pres2 = model_pressure_t(i, j,k,id,x)
- model_pressure = (pres1 + pres2)/2.0_r8
+ model_pressure = interp_pressure(pres1, pres2, vertical=.false.)
! If not periodic, then try extrapolating
else
pres1 = model_pressure_t(i,j, k,id,x)
pres2 = model_pressure_t(i,j+1,k,id,x)
- model_pressure = (3.0_r8*pres1 - pres2)/2.0_r8
+ model_pressure = interp_pressure(pres1, pres2, extrapolate=.true., vertical=.false.)
endif
@@ -5152,7 +5334,7 @@
pres1 = model_pressure_t(i,j, k,id,x)
pres2 = model_pressure_t(i,j-1,k,id,x)
- model_pressure = (pres1 + pres2)/2.0_r8
+ model_pressure = interp_pressure(pres1, pres2, vertical=.false.)
endif
@@ -5202,14 +5384,14 @@
! We are at seam in longitude, take first and last M-grid points
pres1 = model_pressure_s(i-1,j,id,x)
pres2 = model_pressure_s(1, j,id,x)
- model_surface_pressure = (pres1 + pres2)/2.0_r8
+ model_surface_pressure = interp_pressure(pres1, pres2, vertical=.false.)
else
! If not periodic, then try extrapolating
pres1 = model_pressure_s(i-1,j,id,x)
pres2 = model_pressure_s(i-2,j,id,x)
- model_surface_pressure = (3.0_r8*pres1 - pres2)/2.0_r8
+ model_surface_pressure = interp_pressure(pres1, pres2, extrapolate=.true., vertical=.false.)
endif
@@ -5221,14 +5403,14 @@
! We are at seam in longitude, take first and last M-grid points
pres1 = model_pressure_s(i, j,id,x)
pres2 = model_pressure_s(wrf%dom(id)%we,j,id,x)
- model_surface_pressure = (pres1 + pres2)/2.0_r8
+ model_surface_pressure = interp_pressure(pres1, pres2, vertical=.false.)
else
! If not periodic, then try extrapolating
pres1 = model_pressure_s(i, j,id,x)
pres2 = model_pressure_s(i+1,j,id,x)
- model_surface_pressure = (3.0_r8*pres1 - pres2)/2.0_r8
+ model_surface_pressure = interp_pressure(pres1, pres2, extrapolate=.true., vertical=.false.)
endif
@@ -5236,7 +5418,7 @@
pres1 = model_pressure_s(i, j,id,x)
pres2 = model_pressure_s(i-1,j,id,x)
- model_surface_pressure = (pres1 + pres2)/2.0_r8
+ model_surface_pressure = interp_pressure(pres1, pres2, vertical=.false.)
endif
@@ -5255,14 +5437,14 @@
pres1 = model_pressure_s(off,j-1,id,x)
pres2 = model_pressure_s(i ,j-1,id,x)
- model_surface_pressure = (pres1 + pres2)/2.0_r8
+ model_surface_pressure = interp_pressure(pres1, pres2, vertical=.false.)
! If not periodic, then try extrapolating
else
pres1 = model_pressure_s(i,j-1,id,x)
pres2 = model_pressure_s(i,j-2,id,x)
- model_surface_pressure = (3.0_r8*pres1 - pres2)/2.0_r8
+ model_surface_pressure = interp_pressure(pres1, pres2, extrapolate=.true., vertical=.false.)
endif
@@ -5277,14 +5459,14 @@
pres1 = model_pressure_s(off,j,id,x)
pres2 = model_pressure_s(i, j,id,x)
- model_surface_pressure = (pres1 + pres2)/2.0_r8
+ model_surface_pressure = interp_pressure(pres1, pres2, vertical=.false.)
! If not periodic, then try extrapolating
else
pres1 = model_pressure_s(i,j, id,x)
pres2 = model_pressure_s(i,j+1,id,x)
- model_surface_pressure = (3.0_r8*pres1 - pres2)/2.0_r8
+ model_surface_pressure = interp_pressure(pres1, pres2, extrapolate=.true., vertical=.false.)
endif
@@ -5292,7 +5474,7 @@
pres1 = model_pressure_s(i,j, id,x)
pres2 = model_pressure_s(i,j-1,id,x)
- model_surface_pressure = (pres1 + pres2)/2.0_r8
+ model_surface_pressure = interp_pressure(pres1, pres2, vertical=.false.)
endif
@@ -5380,6 +5562,128 @@
!#######################################################
+function interp_pressure(p1, p2, extrapolate, vertical)
+
+! interpolate halfway between pressures 1 and 2 in log units.
+! if extrapolate is true, extrapolate where 1 is the edge and
+! 2 is the inner value, going 1/2 grid cell out.
+
+real(r8), intent(in) :: p1, p2
+logical, intent(in), optional :: extrapolate
+logical, intent(in), optional :: vertical
+real(r8) :: interp_pressure
+
+logical :: do_interp
+logical :: is_vert
+real(r8) :: intermediate
+
+! default is to do interpolation; only extrapolate if the optional
+! arg is specified and if it is true.
+do_interp = .true.
+if (present(extrapolate)) then
+ if (extrapolate) do_interp = .false.
+endif
+
+! if vert is specified and is false, check log_horz_interpM instead
+! of log_vert_interp to decide log vs linear interpolation for the
+! Midpoint value. default is to do vertical interpolation.
+is_vert = .true.
+if (present(vertical)) then
+ is_vert = vertical
+endif
+
+! once we like the results, remove the log_vert_interp test.
+if (do_interp) then
+ if (( is_vert .and. log_vert_interp ) .or. &
+ (.not. is_vert .and. log_horz_interpM)) then
+ interp_pressure = exp((log(p1) + log(p2))/2.0_r8)
+ else
+ interp_pressure = (p1 + p2)/2.0_r8
+ endif
+else
+ if (( is_vert .and. log_vert_interp ) .or. &
+ (.not. is_vert .and. log_horz_interpM)) then
+ intermediate = (3.0_r8*log(p1) - log(p2))/2.0_r8
+ if (intermediate <= 0.0_r8) then
+ interp_pressure = p1
+ else
+ interp_pressure = exp(intermediate)
+ endif
+ else
+ interp_pressure = (3.0_r8*p1 - p2)/2.0_r8
+ endif
+endif
+
+end function interp_pressure
+
+!#######################################################
+
+function interp_4pressure(p1, p2, p3, p4, dx, dxm, dy, dym, extrapolate, edgep)
+
+! given 4 corners of a quad, where the p1, p2, p3 and p4 points are
+! respectively: lower left, lower right, upper left, upper right
+! and dx is the distance in x, dxm is 1.0-dx, dy is distance in y
+! and dym is 1.0-dy, interpolate the pressure while converted to log.
+! if extrapolate is true, extrapolate where edgep is the edge pressure
+! and the 4 points and dx/dy give the location of the inner point.
+
+real(r8), intent(in) :: p1, p2, p3, p4
+real(r8), intent(in) :: dx, dxm, dy, dym
+logical, intent(in), optional :: extrapolate
+real(r8), intent(in), optional :: edgep
+real(r8) :: interp_4pressure
+
+logical :: do_interp
+real(r8) :: intermediate
+real(r8) :: l1, l2, l3, l4
+
+! default is to do interpolation; only extrapolate if the optional
+! arg is specified and if it is true. for extrapolation 'edgep' is
+! required; it is unused for interpolation.
+do_interp = .true.
+if (present(extrapolate)) then
+ if (extrapolate) do_interp = .false.
+endif
+
+if (extrapolate .and. .not. present(edgep)) then
+ call error_handler(E_ERR, 'interp_4pressure:', &
+ 'edgep must be specified for extrapolation. internal error.', &
+ source, revision, revdate)
+endif
+
+if (log_horz_interpQ) then
+ l1 = log(p1)
+ l2 = log(p2)
+ l3 = log(p3)
+ l4 = log(p4)
+endif
+
+! once we like the results, remove the log_horz_interpQ test.
+if (do_interp) then
+ if (log_horz_interpQ) then
+ interp_4pressure = exp(dym*( dxm*l1 + dx*l2 ) + dy*( dxm*l3 + dx*l4 ))
+ else
+ interp_4pressure = dym*( dxm*p1 + dx*p2 ) + dy*( dxm*p3 + dx*p4 )
+ endif
+else
+ if (log_horz_interpQ) then
+ intermediate = (3.0_r8*log(edgep) - &
+ dym*( dxm*l1 + dx*l2 ) - dy*( dxm*l3 + dx*l4 ))/2.0_r8
+ if (intermediate <= 0.0_r8) then
+ interp_4pressure = edgep
+ else
+ interp_4pressure = exp(intermediate)
+ endif
+ else
+ interp_4pressure = (3.0_r8*edgep - &
+ dym*( dxm*p1 + dx*p2 ) - dy*( dxm*p3 + dx*p4 ))/2.0_r8
+ endif
+endif
+
+end function interp_4pressure
+
+!#######################################################
+
function model_rho_t(i,j,k,id,x)
! Calculate the total density on mass point (half (mass) levels, T-point).
Modified: DART/branches/development/obs_kind/DEFAULT_obs_kind_mod.F90
===================================================================
--- DART/branches/development/obs_kind/DEFAULT_obs_kind_mod.F90 2011-12-31 00:03:46 UTC (rev 5474)
+++ DART/branches/development/obs_kind/DEFAULT_obs_kind_mod.F90 2012-01-04 22:48:33 UTC (rev 5475)
@@ -212,9 +212,19 @@
KIND_SMOKE = 98, &
KIND_SEASALT = 99
+! kinds for ZVD (advanced microphysics)
+integer, parameter, public ::&
+ KIND_HAIL_MIXING_RATIO = 100, &
+ KIND_HAIL_NUMBER_CONCENTR = 101, &
+ KIND_GRAUPEL_VOLUME = 102, &
+ KIND_HAIL_VOLUME = 103, &
+ KIND_DIFFERENTIAL_REFLECTIVITY = 104, &
+ KIND_SPECIFIC_DIFFERENTIAL_PHASE = 105, &
+ KIND_FLASH_RATE_2D = 106
+
!! PRIVATE ONLY TO THIS MODULE. see comment below near the max_obs_specific
!! declaration.
-integer, parameter :: max_obs_generic = 99
+integer, parameter :: max_obs_generic = 106
!----------------------------------------------------------------------------
! This list is autogenerated by the 'preprocess' program. To add new
@@ -420,7 +430,6 @@
obs_kind_names(87) = obs_kind_type(KIND_TOTAL_PRECIPITABLE_WATER, 'KIND_TOTAL_PRECIPITABLE_WATER')
obs_kind_names(88) = obs_kind_type(KIND_VERTLEVEL, 'KIND_VERTLEVEL')
obs_kind_names(89) = obs_kind_type(KIND_MICROWAVE_BRIGHT_TEMP, 'KIND_MICROWAVE_BRIGHT_TEMP')
-
obs_kind_names(90) = obs_kind_type(KIND_INTEGRATED_SULFATE, 'KIND_INTEGRATED_SULFATE')
obs_kind_names(91) = obs_kind_type(KIND_INTEGRATED_DUST, 'KIND_INTEGRATED_DUST')
obs_kind_names(92) = obs_kind_type(KIND_INTEGRATED_SMOKE, 'KIND_INTEGRATED_SMOKE')
@@ -431,7 +440,15 @@
obs_kind_names(97) = obs_kind_type(KIND_DUST, 'KIND_DUST')
obs_kind_names(98) = obs_kind_type(KIND_SMOKE, 'KIND_SMOKE')
obs_kind_names(99) = obs_kind_type(KIND_SEASALT, 'KIND_SEASALT')
+obs_kind_names(100) = obs_kind_type(KIND_HAIL_MIXING_RATIO, 'KIND_HAIL_MIXING_RATIO')
+obs_kind_names(101) = obs_kind_type(KIND_HAIL_NUMBER_CONCENTR, 'KIND_HAIL_NUMBER_CONCENTR')
+obs_kind_names(102) = obs_kind_type(KIND_GRAUPEL_VOLUME, 'KIND_GRAUPEL_VOLUME')
+obs_kind_names(103) = obs_kind_type(KIND_HAIL_VOLUME, 'KIND_HAIL_VOLUME')
+obs_kind_names(104) = obs_kind_type(KIND_DIFFERENTIAL_REFLECTIVITY, 'KIND_DIFFERENTIAL_REFLECTIVITY')
+obs_kind_names(105) = obs_kind_type(KIND_SPECIFIC_DIFFERENTIAL_PHASE, 'KIND_SPECIFIC_DIFFERENTIAL_PHASE')
+obs_kind_names(106) = obs_kind_type(KIND_FLASH_RATE_2D, 'KIND_FLASH_RATE_2D')
+
! count here, then output below
num_kind_assimilate = 0
More information about the Dart-dev
mailing list