[Dart-dev] [3894] DART/trunk/models/wrf/model_mod.f90: Same functionality as before, but added better error checks
nancy at ucar.edu
nancy at ucar.edu
Fri May 29 13:18:53 MDT 2009
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20090529/d6d0b356/attachment.html
-------------- next part --------------
Modified: DART/trunk/models/wrf/model_mod.f90
===================================================================
--- DART/trunk/models/wrf/model_mod.f90 2009-05-28 21:43:58 UTC (rev 3893)
+++ DART/trunk/models/wrf/model_mod.f90 2009-05-29 19:18:52 UTC (rev 3894)
@@ -54,8 +54,8 @@
KIND_ICE_NUMBER_CONCENTRATION, KIND_GEOPOTENTIAL_HEIGHT, &
KIND_POTENTIAL_TEMPERATURE, KIND_SOIL_MOISTURE, &
KIND_VORTEX_LAT, KIND_VORTEX_LON, &
- KIND_VORTEX_PMIN, KIND_VORTEX_WMAX,&
- get_raw_obs_kind_index
+ KIND_VORTEX_PMIN, KIND_VORTEX_WMAX, &
+ get_raw_obs_kind_index, max_obs_kinds, get_raw_obs_kind_name
!nc -- module_map_utils split the declarations of PROJ_* into a separate module called
@@ -188,6 +188,7 @@
character(len = 20) :: wrf_nml_file = 'namelist.input'
logical :: have_wrf_nml_file = .false.
+logical :: in_state_vector(max_obs_kinds) = .false.
!-----------------------------------------------------------------------
@@ -235,6 +236,11 @@
! real(r8), dimension(:,:), pointer :: mapfac_m, mapfac_u, mapfac_v
real(r8), dimension(:,:,:), pointer :: phb
+ ! NEWVAR: Currently you have to add a new type here if you want to use
+ ! NEWVAR: a WRF variable which is not one of these types. This will go
+ ! NEWVAR: away eventually, we hope. Search for NEWVAR for other places
+ ! NEWVAR: the code has to change.
+
! JPH local variables to hold type indices
integer :: type_u, type_v, type_w, type_t, type_qv, type_qr, &
type_qc, type_qg, type_qi, type_qs, type_gz
@@ -324,6 +330,15 @@
endif
!---------------------------
+! set this array so we know exactly which obs kinds are
+! allowed to be interpolated (and can give a reasonably
+! helpful error message if not).
+!---------------------------
+
+call fill_dart_kinds_table(wrf_state_variables, in_state_vector)
+
+
+!---------------------------
! next block to be obsolete
!---------------------------
wrf%dom(:)%n_moist = num_moist_vars
@@ -554,6 +569,11 @@
wrf%dom(id)%domain_size = dart_index - 1 - wrf%dom(id-1)%domain_size
end if
+
+ ! NEWVAR: If you add a new wrf array type which is not yet in this list, currently
+ ! NEWVAR: you will have to add it here, and add a type_xx for it, and also add
+ ! NEWVAR: a in_state_variable case in the select statement. search for NEWVAR.
+
! JPH now that we have the domain ID just go ahead and get type indices once
! NOTE: this is not strictly necessary - can use only stagger info in the future (???)
wrf%dom(id)%type_u = get_type_ind_from_type_string(id,'U')
@@ -585,7 +605,8 @@
wrf%model_size = dart_index - 1
allocate (ens_mean(wrf%model_size))
-if(debug) write(*,*) ' wrf model size is ',wrf%model_size
+write(errstring,*) ' wrf model size is ',wrf%model_size
+call error_handler(E_MSG, 'static_init_model', errstring)
end subroutine static_init_model
@@ -869,6 +890,15 @@
! Otherwise, we need to do interpolation
else
+ ! Is this a valid kind to interpolate? Set up in the static_init_model code,
+ ! based on entries in wrf_state_vector namelist item.
+ if (.not. in_state_vector(obs_kind)) then
+ write(errstring, *) 'cannot interpolate ' // trim(get_raw_obs_kind_name(obs_kind)) &
+ // ' with the current WRF arrays in state vector'
+ call error_handler(E_ERR, 'model_interpolate', errstring, &
+ source, revision, revdate)
+ endif
+
! Unravel location_type information
xyz_loc = get_location(location)
@@ -4877,8 +4907,6 @@
! within filter_assim. In other words, these modifications will only matter within
! filter_assim, but will not propagate backwards to filter.
-implicit none
-
type(get_close_type), intent(in) :: gc
type(location_type), intent(inout) :: base_obs_loc, obs_loc(:)
integer, intent(in) :: base_obs_kind, obs_kind(:)
@@ -5226,8 +5254,6 @@
! and then returns the four cornering gridpoints' 2-element integer indices.
subroutine getCorners(i, j, id, type, ll, ul, lr, ur, rc)
- implicit none
-
integer, intent(in) :: i, j, id, type
integer, dimension(2), intent(out) :: ll, ul, lr, ur
integer, intent(out) :: rc
@@ -5477,8 +5503,6 @@
! ncid: input, file handl
! id: input, domain id
-implicit none
-
integer, intent(in) :: ncid
integer, intent(out) :: bt,bts,sn,sns,we,wes,sls
logical, parameter :: debug = .false.
@@ -5541,8 +5565,6 @@
! ncid: input, file handl
! id: input, domain id
-implicit none
-
integer, intent(in) :: ncid, id
logical, parameter :: debug = .false.
@@ -5554,7 +5576,8 @@
'static_init_model', 'get_att DY')
call nc_check( nf90_get_att(ncid, nf90_global, 'DT', wrf%dom(id)%dt), &
'static_init_model', 'get_att DT')
- if (do_output()) print*,'dt from wrfinput_d0X file is: ', wrf%dom(id)%dt
+ write(errstring, *) 'dt from wrfinput_d0X file is: ', wrf%dom(id)%dt
+ call error_handler(E_MSG, ' ', errstring)
if(debug) write(*,*) ' dx, dy, dt are ',wrf%dom(id)%dx, &
wrf%dom(id)%dy, wrf%dom(id)%dt
@@ -5593,8 +5616,6 @@
! id: input, domain id
-implicit none
-
integer, intent(in) :: id
logical, parameter :: debug = .false.
@@ -5639,8 +5660,6 @@
! ncid: input, file handle
! id: input, domain id
-implicit none
-
integer, intent(in) :: ncid, id
logical, parameter :: debug = .false.
integer :: var_id
@@ -5761,8 +5780,6 @@
! id: input, domain id
-implicit none
-
integer, intent(in) :: id
logical, parameter :: debug = .false.
@@ -5840,8 +5857,6 @@
subroutine fill_default_state_table(default_table)
-implicit none
-
character(len=129), intent(out) :: default_table(num_state_table_columns,max_state_variables)
integer :: row
@@ -5894,10 +5909,111 @@
!--------------------------------------------
!--------------------------------------------
+subroutine fill_dart_kinds_table(wrf_state_variables, in_state_vector)
+
+! for each row in the kinds table, tick them off in an array
+! of all possible kinds. then do some simple error checks for
+! kinds we know have interactions -- like both wind vectors are
+! required to convert directions from projection to lat/lon, etc
+
+character(len=*), intent(in) :: wrf_state_variables(:, :)
+logical, intent(inout) :: in_state_vector(:)
+
+integer :: row, i, nextkind
+
+row = size(wrf_state_variables, 2)
+
+! NEWVAR: see each of part1, part 2, and part 3 below.
+! NEWVAR: if you are adding support for a new kind, see if there are
+! NEWVAR: any interactions that need special case code added.
+
+! part 1: mark off all the kinds that the user specifies, plus the
+! kinds that are related and can be interpolated from the given kind.
+
+do i = 1, row
+
+ ! end of the list?
+ if (wrf_state_variables(2, i) == 'NULL') exit
+
+ nextkind = get_raw_obs_kind_index(trim(wrf_state_variables(2, i)))
+ select case(nextkind)
+
+ ! wrf stores potential temperature (temperature perturbations around a
+ ! threshold) but we can interpolate sensible temperature from it
+ case (KIND_POTENTIAL_TEMPERATURE)
+ in_state_vector(KIND_TEMPERATURE) = .true.
+ in_state_vector(KIND_POTENTIAL_TEMPERATURE) = .true.
+
+ ! we use vapor mixing ratio to compute specific humidity
+ case (KIND_VAPOR_MIXING_RATIO)
+ in_state_vector(KIND_VAPOR_MIXING_RATIO) = .true.
+ in_state_vector(KIND_SPECIFIC_HUMIDITY) = .true.
+
+ case default
+ in_state_vector(nextkind) = .true.
+
+ end select
+enddo
+
+
+! part 2: if you specified one kind but the interpolation is going to
+! require another, make sure the combinations make sense. i.e. if you
+! have U wind, you also have to have V wind, etc.
+
+do i = 1, size(in_state_vector)
+
+ select case(i)
+
+ ! the vortex center computations require wind speeds and phb?
+ case (KIND_VORTEX_LAT, KIND_VORTEX_LON, KIND_VORTEX_PMIN, KIND_VORTEX_WMAX)
+ if ((.not. in_state_vector(KIND_U_WIND_COMPONENT)) .or. &
+ (.not. in_state_vector(KIND_V_WIND_COMPONENT)) .or. &
+ (.not. in_state_vector(KIND_TEMPERATURE)) .or. &
+ (.not. in_state_vector(KIND_VAPOR_MIXING_RATIO)) .or. &
+ (.not. in_state_vector(KIND_PRESSURE))) then
+ write(errstring, *) 'VORTEX kinds require U,V,T,QVAPOR,MU in state vector'
+ call error_handler(E_ERR, 'fill_dart_kinds_table', errstring, &
+ source, revision, revdate)
+ endif
+
+ ! if you have one wind component you have to have both
+ case (KIND_U_WIND_COMPONENT, KIND_V_WIND_COMPONENT)
+ if ((.not. in_state_vector(KIND_U_WIND_COMPONENT)) .or. &
+ (.not. in_state_vector(KIND_V_WIND_COMPONENT))) then
+ write(errstring, *) 'WIND kinds require both U,V in state vector'
+ call error_handler(E_ERR, 'fill_dart_kinds_table', errstring, &
+ source, revision, revdate)
+ endif
+
+ ! by default anything else is fine
+
+ end select
+enddo
+
+
+! part 3: fields you just have to have, always.
+if (.not. in_state_vector(KIND_GEOPOTENTIAL_HEIGHT)) then
+ write(errstring, *) 'PH is always a required field'
+ call error_handler(E_ERR, 'fill_dart_kinds_table', errstring, &
+ source, revision, revdate)
+endif
+
+! FIXME: is this true? or is pressure always required, and surface
+! pressure required only if you have any of the surface obs?
+if ((.not. in_state_vector(KIND_PRESSURE)) .and. &
+ (.not. in_state_vector(KIND_SURFACE_PRESSURE))) then
+ write(errstring, *) 'One of MU or PSFC is a required field'
+ call error_handler(E_ERR, 'fill_dart_kinds_table', errstring, &
+ source, revision, revdate)
+endif
+
+end subroutine fill_dart_kinds_table
+
+!--------------------------------------------
+!--------------------------------------------
+
integer function get_number_of_wrf_variables(id, state_table, var_element_list)
-implicit none
-
integer, intent(in) :: id
character(len=*), intent(in) :: state_table(num_state_table_columns,max_state_variables)
integer, intent(out), optional :: var_element_list(max_state_variables)
@@ -5924,7 +6040,7 @@
enddo ! ivar
if ( debug ) then
- print*,'In functon get_number_of_wrf_variables'
+ print*,'In function get_number_of_wrf_variables'
print*,'Found ',num_vars,' state variables for domain id ',id
endif
@@ -5939,8 +6055,6 @@
subroutine set_variable_bound_defaults(nbounds,lb,ub,instructions)
- implicit none
-
integer, intent(in) :: nbounds
real(r8), dimension(nbounds), intent(out) :: lb, ub
character(len=10), dimension(nbounds), intent(out) :: instructions
@@ -5961,8 +6075,6 @@
! matches WRF variable name in bounds table to input name, and assigns
! the bounds and instructions if they exist
- implicit none
-
character(len=*), intent(in) :: bounds_table(num_bounds_table_columns,max_state_variables)
character(len=*), intent(in) :: wrf_var_name
real(r8), intent(out) :: lb,ub
@@ -6025,8 +6137,6 @@
logical function variable_is_on_domain(domain_id_string, id)
-implicit none
-
integer, intent(in) :: id
character(len=*), intent(in) :: domain_id_string
@@ -6061,8 +6171,6 @@
! ncid: input, file handle
! id: input, domain index
-implicit none
-
integer, intent(in) :: ncid, id
integer, intent(in) :: bt, bts, sn, sns, we, wes
character(len=*), intent(in) :: wrf_var_name
@@ -6124,8 +6232,6 @@
! ncid: input, file handle
! id: input, domain index
-implicit none
-
integer, intent(in) :: ncid, id
character(len=*), intent(in) :: wrf_var_name
character(len=129), intent(out) :: description, units
@@ -6159,8 +6265,6 @@
! simply loop through the state variable table to get the index of the
! type for this domain. returns -1 if not found
- implicit none
-
integer, intent(in) :: id
character(len=*), intent(in) :: wrf_varname
@@ -6195,8 +6299,6 @@
subroutine trans_2Dto1D( a1d, a2d, nx, ny )
-implicit none
-
integer, intent(in) :: nx,ny
real(r8), intent(inout) :: a1d(:)
real(r8), intent(inout) :: a2d(nx,ny)
@@ -6228,8 +6330,6 @@
subroutine trans_3Dto1D( a1d, a3d, nx, ny, nz )
-implicit none
-
integer, intent(in) :: nx,ny,nz
real(r8), intent(inout) :: a1d(:)
real(r8), intent(inout) :: a3d(:,:,:)
@@ -6265,8 +6365,6 @@
subroutine trans_1Dto2D( a1d, a2d, nx, ny )
-implicit none
-
integer, intent(in) :: nx,ny
real(r8), intent(inout) :: a1d(:)
real(r8), intent(inout) :: a2d(nx,ny)
@@ -6298,8 +6396,6 @@
subroutine trans_1Dto3D( a1d, a3d, nx, ny, nz )
-implicit none
-
integer, intent(in) :: nx,ny,nz
real(r8), intent(inout) :: a1d(:)
real(r8), intent(inout) :: a3d(:,:,:)
@@ -6335,7 +6431,6 @@
subroutine get_wrf_date (tstring, year, month, day, hour, minute, second)
-implicit none
!--------------------------------------------------------
! Returns integers taken from tstring
! It is assumed that the tstring char array is as YYYY-MM-DD_hh:mm:ss
@@ -6358,7 +6453,6 @@
subroutine set_wrf_date (tstring, year, month, day, hour, minute, second)
-implicit none
!--------------------------------------------------------
! Returns integers taken from tstring
! It is assumed that the tstring char array is as YYYY-MM-DD_hh:mm:ss
More information about the Dart-dev
mailing list