[Dart-dev] [5899] DART/branches/development/models/mpas_ocn: Updated the mpas_ocn to be consistent with the latest mpas_atm.
nancy at ucar.edu
nancy at ucar.edu
Thu Oct 18 08:34:26 MDT 2012
Revision: 5899
Author: thoar
Date: 2012-10-18 08:34:25 -0600 (Thu, 18 Oct 2012)
Log Message:
-----------
Updated the mpas_ocn to be consistent with the latest mpas_atm.
The mpas_atm has seen a lot of work, the mpas_ocn should be
considered a stub at this time.
There are a LOT of references to atmospheric items, I have
commented most of these out and put in fatal errors should these
be encountered.
I reduced the model_interpolate() routine to a stub just to get
things going. Not sure how applicable the atmospheric routine
is to the ocean. The mpas_atm/model_mod:model_interpolate() is
untouched. We can borrow from there as need be.
The code compiles with ifort (IFORT) 11.1 2010080 and throws
a ton of unused variables etc. It is anticipated there will be
a lot of shared code between the mpas_atm and mpas_ocn. When
we consolidate the common components - we can cut out the dead wood.
Modified Paths:
--------------
DART/branches/development/models/mpas_ocn/model_mod.f90
DART/branches/development/models/mpas_ocn/work/input.nml
DART/branches/development/models/mpas_ocn/work/path_names_create_fixed_network_seq
DART/branches/development/models/mpas_ocn/work/path_names_create_obs_sequence
DART/branches/development/models/mpas_ocn/work/path_names_dart_to_model
DART/branches/development/models/mpas_ocn/work/path_names_filter
DART/branches/development/models/mpas_ocn/work/path_names_integrate_model
DART/branches/development/models/mpas_ocn/work/path_names_model_mod_check
DART/branches/development/models/mpas_ocn/work/path_names_model_to_dart
DART/branches/development/models/mpas_ocn/work/path_names_obs_diag
DART/branches/development/models/mpas_ocn/work/path_names_obs_sequence_tool
DART/branches/development/models/mpas_ocn/work/path_names_perfect_model_obs
DART/branches/development/models/mpas_ocn/work/path_names_restart_file_tool
-------------- next part --------------
Modified: DART/branches/development/models/mpas_ocn/model_mod.f90
===================================================================
--- DART/branches/development/models/mpas_ocn/model_mod.f90 2012-10-17 21:46:45 UTC (rev 5898)
+++ DART/branches/development/models/mpas_ocn/model_mod.f90 2012-10-18 14:34:25 UTC (rev 5899)
@@ -10,9 +10,17 @@
! $Revision$
! $Date$
-! This is the interface between the model model and DART.
+! MPAS ocean model interface to the DART data assimilation system.
+! Code in this module is compiled with the DART executables. It isolates
+! all information about the MPAS grids, model variables, and other details.
+! There are a set of 16 subroutine interfaces that are required by DART;
+! these cannot be changed. Additional public routines in this file can
+! be used by converters and utilities and those interfaces can be anything
+! that is useful to other pieces of code.
-! Modules that are absolutely required for use are listed
+
+! Routines in other modules that are used here.
+
use types_mod, only : r4, r8, digits12, SECPERDAY, MISSING_R8, &
rad2deg, deg2rad, PI, MISSING_I
use time_manager_mod, only : time_type, set_time, set_date, get_date, get_time,&
@@ -24,14 +32,21 @@
use location_mod, only : location_type, get_dist, query_location, &
get_close_maxdist_init, get_close_type, &
set_location, get_location, horiz_dist_only, &
- vert_is_undef, VERTISUNDEF, &
- vert_is_surface, VERTISSURFACE, &
- vert_is_level, VERTISLEVEL, &
- vert_is_pressure, VERTISPRESSURE, &
- vert_is_height, VERTISHEIGHT, &
+ write_location, &
+ vert_is_undef, VERTISUNDEF, &
+ vert_is_surface, VERTISSURFACE, &
+ vert_is_level, VERTISLEVEL, &
+ vert_is_pressure, VERTISPRESSURE, &
+ vert_is_height, VERTISHEIGHT, &
+ vert_is_scale_height, VERTISSCALEHEIGHT, &
get_close_obs_init, get_close_obs_destroy, &
loc_get_close_obs => get_close_obs
+use xyz_location_mod, only : xyz_location_type, xyz_get_close_maxdist_init, &
+ xyz_get_close_type, xyz_set_location, xyz_get_location, &
+ xyz_get_close_obs_init, xyz_get_close_obs_destroy, &
+ xyz_find_nearest
+
use utilities_mod, only : register_module, error_handler, &
E_ERR, E_WARN, E_MSG, logfileunit, get_unit, &
nc_check, do_output, to_upper, nmlfileunit, &
@@ -39,26 +54,33 @@
open_file, file_exist, find_textfile_dims, &
file_to_text, close_file, do_nml_file, do_nml_term
-use obs_kind_mod, only : paramname_length, &
- get_raw_obs_kind_index, &
- get_raw_obs_kind_name, &
- KIND_SURFACE_PRESSURE, &
- KIND_VERTICAL_VELOCITY, &
- KIND_POTENTIAL_TEMPERATURE, &
- KIND_EDGE_NORMAL_SPEED, &
- KIND_TEMPERATURE, &
- KIND_U_WIND_COMPONENT, &
- KIND_V_WIND_COMPONENT, &
- KIND_PRESSURE, &
- KIND_DENSITY, &
- KIND_VAPOR_MIXING_RATIO
+use obs_kind_mod, only : paramname_length, &
+ get_raw_obs_kind_index, &
+ get_raw_obs_kind_name, &
+ KIND_TEMPERATURE, &
+ KIND_SALINITY, &
+ KIND_DRY_LAND, &
+ KIND_U_CURRENT_COMPONENT, &
+ KIND_V_CURRENT_COMPONENT, &
+ KIND_SEA_SURFACE_HEIGHT, &
+ KIND_SEA_SURFACE_PRESSURE
use mpi_utilities_mod, only: my_task_id
use random_seq_mod, only: random_seq_type, init_random_seq, random_gaussian
+! netcdf modules
use typesizes
use netcdf
+
+! RBF (radial basis function) modules, donated by LANL. currently deprecated
+! in this version. they did the job but not as well as other techniques and
+! at a much greater execution-time code. they were used to interpolate
+! values at arbitrary locations, not just at cell centers. with too small
+! a set of basis points, the values were discontinuous at cell boundaries;
+! with too many the values were too smoothed out. we went back to
+! barycentric interpolation in triangles formed by the three cell centers
+! that enclosed the given point.
use get_geometry_mod
use get_reconstruct_mod
@@ -103,7 +125,9 @@
revision = '$Revision$', &
revdate = '$Date$'
-character(len=256) :: string1, string2
+! module global storage; maintains values between calls, accessible by
+! any subroutine
+character(len=256) :: string1, string2, string3
logical, save :: module_initialized = .false.
! Real (physical) constants as defined exactly in MPAS.
@@ -114,47 +138,78 @@
real(r8), parameter :: p0 = 100000.0_r8
real(r8), parameter :: rcv = rgas/(cp-rgas)
+! earth radius; needed to convert lat/lon to x,y,z cartesian coords.
! FIXME: one of the example ocean files had a global attr with 6371220.0
! instead of 1229. ??
real(r8), parameter :: radius = 6371229.0 ! meters
-! Storage for a random sequence for perturbing a single initial state
+! roundoff error
+real(r8), parameter :: roundoff = 1.0e-12_r8
+! Storage for a random sequence for perturbing a single initial state
type(random_seq_type) :: random_seq
! Structure for computing distances to cell centers, and assorted arrays
! needed for the get_close code.
-type(get_close_type) :: cc_gc
-type(location_type), allocatable :: cell_locs(:)
-integer, allocatable :: dummy(:), close_ind(:)
+type(xyz_get_close_type) :: cc_gc
+type(xyz_location_type), allocatable :: cell_locs(:)
-! things which can/should be in the model_nml
-
+! variables which are in the module namelist
+integer :: vert_localization_coord = VERTISHEIGHT
integer :: assimilation_period_days = 0
integer :: assimilation_period_seconds = 60
real(r8) :: model_perturbation_amplitude = 0.0001 ! tiny amounts
-logical :: output_state_vector = .true.
+real(r8) :: highest_obs_pressure_mb = 100.0_r8 ! do not assimilate obs higher than this level.
+real(r8) :: cell_size_meters = 100.0_r8 ! size of largest cell in r
+logical :: output_state_vector = .false. ! output prognostic variables (if .false.)
+logical :: logp = .false. ! if true, interpolate pres in log space
integer :: debug = 0 ! turn up for more and more debug messages
+integer :: xyzdebug = 0
character(len=32) :: calendar = 'Gregorian'
character(len=256) :: model_analysis_filename = 'mpas_analysis.nc'
character(len=256) :: grid_definition_filename = 'mpas_analysis.nc'
-logical :: use_new_code = .true.
+! FIXME: these are here to allow regression testing and should be
+! removed when we settle on the final version. for sure the 'old code'
+! for locating a point in a group of triangles should be removed
+! because there are known errors at high latitudes near the poles.
+
+! if .false. use U/V reconstructed winds tri interp at centers for wind forward ops
+! if .true. use edge normal winds (u) with RBF functs for wind forward ops
+logical :: use_u_for_wind = .false.
+
+! if using rbf, options 1,2,3 control how many points go into the rbf.
+! larger numbers use more points from further away
+integer :: use_rbf_option = 2
+
+! if .false. edge normal winds (u) should be in state vector and are written out directly
+! if .true. edge normal winds (u) are updated based on U/V reconstructed winds
+logical :: update_u_from_reconstruct = .true.
+
+! only if update_u_from_reconstruct is true,
+! if .false. use the cell center u,v reconstructed values to update edge winds
+! if .true., read in the original u,v winds, compute the increments after the
+! assimilation, and use only the increments to update the edge winds
+logical :: use_increments_for_u_update = .true.
+
namelist /model_nml/ &
model_analysis_filename, &
grid_definition_filename, &
output_state_vector, &
+ vert_localization_coord, &
assimilation_period_days, &
assimilation_period_seconds, &
model_perturbation_amplitude, &
calendar, &
debug, &
- use_new_code
+ xyzdebug, &
+ use_u_for_wind, &
+ use_rbf_option, &
+ update_u_from_reconstruct, &
+ use_increments_for_u_update, &
+ highest_obs_pressure_mb
-!------------------------------------------------------------------
-! DART state vector are specified in the input.nml:mpas_vars_nml namelist.
-!------------------------------------------------------------------
-
+! DART state vector contents are specified in the input.nml:&mpas_vars_nml namelist.
integer, parameter :: max_state_variables = 80
integer, parameter :: num_state_table_columns = 2
character(len=NF90_MAX_NAME) :: mpas_state_variables(max_state_variables * num_state_table_columns ) = ' '
@@ -199,6 +254,7 @@
integer :: nVertLevels = -1 ! Vertical levels; count of vert cell centers
integer :: nVertLevelsP1 = -1 ! Vert levels plus 1; count of vert cell faces
integer :: vertexDegree = -1 ! Max number of cells/edges that touch any vertex
+integer :: nSoilLevels = -1 ! Number of soil layers
! scalar grid positions
@@ -242,22 +298,23 @@
integer, allocatable :: boundaryVertex(:,:)
integer, allocatable :: maxLevelCell(:)
-real(r8), allocatable :: ens_mean(:) ! may be needed for forward ops
+real(r8), allocatable :: ens_mean(:) ! needed to convert vertical distances consistently
integer :: model_size ! the state vector length
type(time_type) :: model_timestep ! smallest time to adv model
+! useful flags in making decisions when searching for points, etc
logical :: global_grid = .true. ! true = the grid covers the sphere with no holes
logical :: all_levels_exist_everywhere = .true. ! true = cells defined at all levels
-logical :: has_real_u = .false. ! true = has original u on edges
+logical :: has_edge_u = .false. ! true = has original normal u on edges
logical :: has_uvreconstruct = .false. ! true = has reconstructed at centers
-!------------------------------------------------------------------
-! The model analysis manager namelist variables
-!------------------------------------------------------------------
-
+! currently unused; for a regional model it is going to be necessary to know
+! if the grid is continuous around longitudes (wraps in east-west) or not,
+! and if it covers either of the poles.
character(len= 64) :: ew_boundary_type, ns_boundary_type
+! common names that call specific subroutines based on the arg types
INTERFACE vector_to_prog_var
MODULE PROCEDURE vector_to_1d_prog_var
MODULE PROCEDURE vector_to_2d_prog_var
@@ -333,9 +390,10 @@
character(len=paramname_length) :: kind_string
integer :: ncid, VarID, numdims, varsize, dimlen
integer :: iunit, io, ivar, i, index1, indexN, iloc, kloc
-integer :: ss, dd
+integer :: ss, dd, z1, m1
integer :: nDimensions, nVariables, nAttributes, unlimitedDimID, TimeDimID
integer :: cel1, cel2
+logical :: both
if ( module_initialized ) return ! only need to do this once.
@@ -383,7 +441,7 @@
! 3) read them from the analysis file
! read_grid_dims() fills in the following module global variables:
-! nCells, nVertices, nEdges, maxEdges, nVertLevels, nVertLevelsP1, vertexDegree
+! nCells, nVertices, nEdges, maxEdges, nVertLevels, nVertLevelsP1, vertexDegree, nSoilLevels
call read_grid_dims()
allocate(latCell(nCells), lonCell(nCells))
@@ -534,6 +592,8 @@
progvar(ivar)%numedges = dimlen
case ('nVertL') ! nVertLevels, nVertLevelsP1, nVertLevelsP2
progvar(ivar)%numvertical = dimlen
+ case ('nSoilL') ! nSoilLevels
+ progvar(ivar)%numvertical = dimlen
end select
enddo DimensionLoop
@@ -546,7 +606,7 @@
progvar(ivar)%ZonHalf = .FALSE.
endif
- if (varname == 'u') has_real_u = .true.
+ if (varname == 'u') has_edge_u = .true.
if (varname == 'uReconstructZonal' .or. &
varname == 'uReconstructMeridional') has_uvreconstruct = .true.
@@ -555,7 +615,7 @@
progvar(ivar)%indexN = index1 + varsize - 1
index1 = index1 + varsize ! sets up for next variable
- !if ( debug > 0 ) call dump_progvar(ivar)
+ if ( debug > 4 ) call dump_progvar(ivar)
enddo
@@ -567,10 +627,10 @@
if ( debug > 0 .and. do_output()) then
write(logfileunit,*)
write( * ,*)
- write(logfileunit,'(" static_init_model: nCells, nVertices, nVertLevels =",3(1x,i6))') &
- nCells, nVertices, nVertLevels
- write( * ,'(" static_init_model: nCells, nVertices, nVertLevels =",3(1x,i6))') &
- nCells, nVertices, nVertLevels
+ write(logfileunit,'(" static_init_model: nCells, nEdges, nVertices, nVertLevels =",4(1x,i6))') &
+ nCells, nEdges, nVertices, nVertLevels
+ write( * ,'(" static_init_model: nCells, nEdges, nVertices, nVertLevels =",4(1x,i6))') &
+ nCells, nEdges, nVertices, nVertLevels
write(logfileunit, *)'static_init_model: model_size = ', model_size
write( * , *)'static_init_model: model_size = ', model_size
if ( global_grid ) then
@@ -589,11 +649,75 @@
endif
endif
-allocate( ens_mean(model_size) )
+! do some sanity checking here:
-! Initialize the interpolation data structures
-call init_interp()
+! if you have at least one of these wind components in the state vector,
+! you have to have them both. the subroutine will error out if only one
+! is found and not both.
+if (has_uvreconstruct) then
+ call winds_present(z1,m1,both)
+endif
+! if you set the namelist to use the reconstructed cell center winds,
+! they have to be in the state vector.
+if (update_u_from_reconstruct .and. .not. has_uvreconstruct) then
+ write(string1,*) 'update_u_from_reconstruct cannot be True'
+ write(string2,*) 'because state vector does not contain U/V reconstructed winds'
+ call error_handler(E_ERR,'static_init_model',string1,source,revision,revdate, &
+ text2=string2)
+endif
+
+! if you set the namelist to update the edge normal winds based on the
+! updated cell center winds, and you also have the edge normal winds in
+! the state vector, warn that the edge normal winds will be ignored
+! when going back to the mpas_analysis.nc file. not an error, but the
+! updates to the edge normal vectors won't affect the results.
+if (update_u_from_reconstruct .and. has_edge_u) then
+ write(string1,*) 'edge normal winds (u) in MPAS file will be updated based on U/V reconstructed winds'
+ write(string2,*) 'and not from the updated edge normal values in the state vector'
+ write(string3,*) 'because update_u_from_reconstruct is True'
+ call error_handler(E_MSG,'static_init_model',string1,source,revision,revdate, &
+ text2=string2, text3=string3)
+endif
+
+! there are two choices when updating the edge normal winds based on the
+! winds at the cell centers. one is a direct interpolation of the values;
+! the other is to read in the original cell centers, compute the increments
+! changed by the interpolation, and then add or substract only the increments
+! from the original edge normal wind values.
+if (update_u_from_reconstruct) then
+ if (use_increments_for_u_update) then
+ write(string1,*) 'edge normal winds (u) in MPAS file will be updated based on the difference'
+ write(string2,*) 'between the original U/V reconstructed winds and the updated values'
+ write(string3,*) 'because use_increment_for_u_update is True'
+ else
+ write(string1,*) 'edge normal winds (u) in MPAS file will be updated by averaging the'
+ write(string2,*) 'updated U/V reconstructed winds at the corresponding cell centers'
+ write(string3,*) 'because use_increment_for_u_update is False'
+ endif
+ call error_handler(E_MSG,'static_init_model',string1,source,revision,revdate, &
+ text2=string2, text3=string3)
+endif
+
+! basically we cannot do much without having at least these
+! three fields in the state vector. refuse to go further
+! if these are not present:
+!
+! TJH if ((get_progvar_index_from_kind(KIND_POTENTIAL_TEMPERATURE) < 0) .or. &
+! TJH (get_progvar_index_from_kind(KIND_DENSITY) < 0) .or. &
+! TJH (get_progvar_index_from_kind(KIND_VAPOR_MIXING_RATIO) < 0)) then
+! TJH write(string1, *) 'State vector is missing one or more of the following fields:'
+! TJH write(string2, *) 'Potential Temperature (theta), Density (rho), Vapor Mixing Ratio (qv).'
+! TJH write(string3, *) 'Cannot convert between height/pressure nor compute sensible temps.'
+! TJH call error_handler(E_ERR,'static_init_model',string1,source,revision,revdate, &
+! TJH text2=string2, text3=string3)
+! TJH endif
+
+string1 = 'fix block of required variables - detritus from atmosphere'
+call error_handler(E_ERR,'static_init_model',string1,source,revision,revdate)
+
+allocate( ens_mean(model_size) )
+
end subroutine static_init_model
@@ -615,7 +739,9 @@
integer :: nxp, nzp, iloc, vloc, nf, n
integer :: myindx
+integer :: istatus
real(r8) :: height
+type(location_type) :: new_location
if ( .not. module_initialized ) call static_init_model
@@ -671,8 +797,6 @@
endif
endif
-
-
if (progvar(nf)%numedges /= MISSING_I) then
if (nzp <= 1) then
location = set_location(lonEdge(iloc),latEdge(iloc), height, VERTISSURFACE)
@@ -687,8 +811,21 @@
endif
endif
-if (debug > 9) then
+! Let us return the vert location with the requested vertical localization coordinate
+! hoping that the code can run faster when same points are close to many obs
+! since vert_convert does not need to be called repeatedly in get_close_obs any more.
+! FIXME: we should test this to see if it's a win (in computation time). for obs
+! which are not dense relative to the grid, this might be slower than doing the
+! conversions on demand in the localization code (in get_close_obs()).
+if ( .not. horiz_dist_only .and. vert_localization_coord /= VERTISHEIGHT ) then
+ new_location = location
+ call vert_convert(ens_mean, new_location, progvar(nf)%dart_kind, vert_localization_coord, istatus)
+ if(istatus == 0) location = new_location
+endif
+
+if (debug > 12) then
+
write(*,'("INDEX_IN / myindx / IVAR / NX, NZ: ",2(i10,2x),3(i5,2x))') index_in, myindx, nf, nxp, nzp
write(*,'(" ILOC, KLOC: ",2(i5,2x))') iloc, vloc
write(*,'(" LON/LAT/HGT: ",3(f12.3,2x))') lonCell(iloc), latCell(iloc), height
@@ -710,121 +847,6 @@
! interpolated value at that location, and an error code. 0 is success,
! anything positive is an error. (negative reserved for system use)
!
-! for specific error codes, see the local_interpolate() comments
-
-! passed variables
-
-real(r8), intent(in) :: x(:)
-type(location_type), intent(in) :: location
-integer, intent(in) :: obs_type
-real(r8), intent(out) :: interp_val
-integer, intent(out) :: istatus
-
-! local storage
-
-integer :: obs_kinds(3), ivar
-real(r8) :: values(3)
-
-! call the normal interpolate code. if it fails because
-! the kind doesn't exist directly in the state vector, try a few
-! kinds we know how to convert.
-
-interp_val = MISSING_R8
-istatus = 888888 ! must be positive (and integer)
-
-obs_kinds(1) = obs_type
-
-! debug only
-if (.false.) then
- ivar = get_progvar_index_from_kind(obs_kinds(1))
- if ( ivar > 0 .and. debug > 7 ) call dump_progvar(ivar, x)
- ivar = get_progvar_index_from_kind(KIND_POTENTIAL_TEMPERATURE)
- if ( ivar > 0 .and. debug > 7 ) call dump_progvar(ivar, x)
- ivar = get_progvar_index_from_kind(KIND_DENSITY)
- if ( ivar > 0 .and. debug > 7 ) call dump_progvar(ivar, x)
- ivar = get_progvar_index_from_kind(KIND_VAPOR_MIXING_RATIO)
- if ( ivar > 0 .and. debug > 7 ) call dump_progvar(ivar, x)
- ivar = get_progvar_index_from_kind(KIND_EDGE_NORMAL_SPEED)
- if ( ivar > 0 .and. debug > 7 ) call dump_progvar(ivar, x)
-endif
-
-
-call local_interpolate(x, location, 1, obs_kinds, values, istatus)
-if (istatus /= 88) then
- ! this is for debugging - when we're confident the code is
- ! returning consistent values and rc codes, both these tests can
- ! be removed for speed. FIXME.
- if (istatus /= 0 .and. values(1) /= MISSING_R8) then
- write(string1,*) 'interp routine returned a bad status but not a MISSING_R8 value'
- write(string2,*) 'value = ', values(1), ' istatus = ', istatus
- call error_handler(E_ERR,'model_interpolate',string1,source,revision,revdate, &
- text2=string2)
- endif
- if (istatus == 0 .and. values(1) == MISSING_R8) then
- write(string1,*) 'interp routine returned a good status but set value to MISSING_R8'
- call error_handler(E_ERR,'model_interpolate',string1,source,revision,revdate)
- endif
- interp_val = values(1)
-if (debug > 5) print *, 'returning value ', interp_val
- return
-endif
-
-if (obs_type == KIND_TEMPERATURE) then
- obs_kinds(1) = KIND_POTENTIAL_TEMPERATURE
- obs_kinds(2) = KIND_DENSITY
- obs_kinds(3) = KIND_VAPOR_MIXING_RATIO
-
- call local_interpolate(x, location, 3, obs_kinds, values, istatus)
-!print *, '1 local interpolate returns istatus = ', istatus
- if (istatus /= 0) then
- ! this is for debugging - when we're confident the code is
- ! returning consistent values and rc codes, both these tests can
- ! be removed for speed. FIXME.
- if (istatus /= 0 .and. .not. any(values == MISSING_R8)) then
- write(string1,*) 'interp routine returned a bad status but good values'
- call error_handler(E_ERR,'model_interpolate',string1,source,revision,revdate)
- endif
- if (istatus == 0 .and. any(values == MISSING_R8)) then
- write(string1,*) 'interp routine returned a good status but bad values'
- call error_handler(E_ERR,'model_interpolate',string1,source,revision,revdate)
- endif
- interp_val = MISSING_R8
- return
- endif
-
- ! compute sensible temperature as return value
- interp_val = theta_to_tk(values(1), values(2), values(3))
- return
-endif
-
-!print *, '2 local interpolate returns istatus = ', istatus
-
-! add other cases here for kinds we want to handle
-! if (istatus == 88) then
-! could try different interpolating different kinds here.
-! if (obs_type == KIND_xxx) then
-! code goes here
-! endif
-! endif
-
-! shouldn't get here.
-interp_val = MISSING_R8
-istatus = 100
-
-end subroutine model_interpolate
-
-
-!------------------------------------------------------------------
-
-subroutine local_interpolate(x, location, num_kinds, obs_kinds, interp_vals, istatus)
-
-! THIS VERSION IS ONLY CALLED INTERNALLY, so we can mess with the arguments.
-! For a given lat, lon, and height, interpolate the correct state values
-! to that location for the filter from the model state vectors. this version
-! can take multiple kinds at a single location. there is only a single return
-! istatus code - 0 is all interpolations worked, positive means one or more
-! failed.
-!
! ERROR codes:
!
! ISTATUS = 99: general error in case something terrible goes wrong...
@@ -840,295 +862,32 @@
! finding an applicable case.
!
-! Passed variables
+! passed variables
real(r8), intent(in) :: x(:)
type(location_type), intent(in) :: location
-integer, intent(in) :: num_kinds
-integer, intent(in) :: obs_kinds(:)
-real(r8), intent(out) :: interp_vals(:)
+integer, intent(in) :: obs_type
+real(r8), intent(out) :: interp_val
integer, intent(out) :: istatus
-! Local storage
+! local storage
-real(r8) :: loc_array(3), llon, llat, lheight, fract, v_interp
-integer :: i, j, ivar, ier, lower, upper
-integer :: pt_base_offset, density_base_offset, qv_base_offset
+integer :: ivar, obs_kind
+integer :: tvars(3)
+logical :: goodkind
+real(r8) :: values(3), lpres
-real(r8) :: weights(3), lower_interp, upper_interp, ltemp, utemp
-integer :: tri_indices(3), base_offset(num_kinds)
-
if ( .not. module_initialized ) call static_init_model
-! Assume failure. Set return val to missing, then the code can
-! just set istatus to something indicating why it failed, and return.
-! If the interpolation is good, the interp_val will be set to the
-! good value. Make any error codes set here be in the 10s
+interp_val = MISSING_R8
+istatus = 99 ! must be positive (and integer)
-interp_vals(:) = MISSING_R8 ! the DART bad value flag
-istatus = 99 ! unknown error
+write(string1,*) 'routine not written yet for the ocean.'
+call error_handler(E_ERR,'model_interpolate',string1,source,revision,revdate)
-! see if all observation kinds are in the state vector. this sets an
-! error code and returns without a fatal error if any answer is no.
-! exceptions: if the vertical location is specified in pressure
-! we end up computing the sensible temperature as part of converting
-! to height (meters), so if the kind to be computed is sensible temp
-! (which is not in state vector; potential temp is), go ahead and
-! say yes, we know how to compute it. if the vert is any other units
-! then fail and let the calling code call in for the components
-! it needs to compute sensible temp itself. also, if the obs is a
-! wind (or ocean current) component, and the edge normal speeds are
-! in the state vector, we can do it.
+end subroutine model_interpolate
-oktointerp: do i=1, num_kinds
- ivar = get_progvar_index_from_kind(obs_kinds(i))
- if (ivar <= 0) then
- ! exceptions 1 and 2:
- ! FIXME: the new code cannot yet compute sensible temperature in one go. fail for that.
- ! remove the new code test once we add it.
- if ((obs_kinds(i) == KIND_TEMPERATURE) .and. vert_is_pressure(location) .and. &
- .not. use_new_code) cycle oktointerp
- if ((obs_kinds(i) == KIND_U_WIND_COMPONENT) .or. obs_kinds(i) == KIND_V_WIND_COMPONENT) then
- ivar = get_progvar_index_from_kind(KIND_EDGE_NORMAL_SPEED)
- if (ivar > 0) cycle oktointerp
- endif
- istatus = 88 ! this kind not in state vector
- return
- endif
-enddo oktointerp
-
-! Not prepared to do w interpolation at this time
-do i=1, num_kinds
- if(obs_kinds(i) == KIND_VERTICAL_VELOCITY) then
- istatus = 16
- return
- endif
-enddo
-
-! Get the individual locations values
-loc_array = get_location(location)
-llon = loc_array(1)
-llat = loc_array(2)
-lheight = loc_array(3)
-
-if (debug > 5) print *, 'requesting interpolation at ', llon, llat, lheight
-
-
-! if (debug > 2) print *, 'base offset now ', base_offset(1)
-
-! FIXME: this needs to play well with multiple types in the kinds list.
-! for now do only the first one. if 'u' is part of the state vector
-! and if they are asking for U/V components, use the new RBF code.
-! if u isn't in the state vector, default to trying to interpolate
-! in the reconstructed U/V components at the cell centers.
-
-if (use_new_code) then
- kindloop: do i=1, num_kinds
- if ((obs_kinds(i) == KIND_U_WIND_COMPONENT .or. &
- obs_kinds(i) == KIND_V_WIND_COMPONENT) .and. has_real_u) then
- if (obs_kinds(i) == KIND_U_WIND_COMPONENT) then
- call compute_u_with_rbf(x, location, .TRUE., interp_vals(i), istatus)
- else
- call compute_u_with_rbf(x, location, .FALSE., interp_vals(i), istatus)
- endif
-
- else if (obs_kinds(i) /= KIND_TEMPERATURE) then
- ivar = get_progvar_index_from_kind(obs_kinds(i))
- call compute_scalar_with_barycentric(x, location, ivar, interp_vals(i), istatus)
- if (istatus /= 0) return
-
- else if (obs_kinds(i) == KIND_TEMPERATURE) then
- print *, 'need to add case in new code for sensible temperature'
- stop
- ! need to get potential temp, pressure, qv here, but can
- ! use same weights. does this get pushed down into the
- ! scalar code?
- endif
- enddo kindloop
- return
-endif
-
-!! start of original code - should be unneeded if we decide to
-!! use the new code.
-if (debug > 4) print *, 'using original interpolation code'
-
-! Find the start and end offsets for these fields in the state vector x(:)
-do i=1, num_kinds
- if (obs_kinds(i) == KIND_TEMPERATURE) then
- base_offset(i) = 0 ! this won't be used in this case
- else
- call get_index_range(obs_kinds(i), base_offset(i))
- endif
-enddo
-
-! Find the indices of the three cell centers that surround this point in
-! the horizontal along with the barycentric weights.
-call get_cell_indices(llon, llat, tri_indices, weights, ier)
-if (debug > 5) write(*, *) 'tri_inds ', tri_indices
-if (debug > 5) write(*, *) 'weights ', weights
-
-! If istatus is not zero couldn't find a triangle, fail
-if(ier /= 0) then
- istatus = 11
- return
-endif
-
-! FIXME: cannot do both surface pressure and any other 3d field with
-! this code (yet). if num kinds > 1 and any are not surface pressure, this
-! should fail.
-! Surface pressure is a 2D field
-if(obs_kinds(1) == KIND_SURFACE_PRESSURE) then
- lheight = 1.0_r8
- ! the 1 below is max number of vert levels to examine
- call triangle_interp(x, base_offset(1), tri_indices, weights, &
- nint(lheight), 1, interp_vals(1), ier)
- if(ier /= 0) then
- istatus = 13
- else
- istatus = 0
- endif
- return
-endif
-
-! the code below has vert undef returning an error, and vert surface
-! returning level 1 (which i believe is correct - the vertical grid
-! is terrain following). FIXME:
-! vert undef could return anything in the column in theory, right?
-
-if (vert_is_undef(location)) then
- istatus = 12
- return
-endif
-
-if(vert_is_surface(location)) then
- do j=1, num_kinds
- lheight = 1.0_r8 ! first grid level is surface
- call triangle_interp(x, base_offset(j), tri_indices, weights, &
- nint(lheight), nVertLevels, interp_vals(j), ier)
- if(ier /= 0) then
- istatus = 13
- return
- endif
- enddo
- istatus = 0
- return
-endif
-
-! If vertical is on a model level don't need vertical interpolation either
-! (FIXME: since the vertical value is a real, in the wrf model_mod they
-! support non-integer levels and interpolate in the vertical just like
-! any other vertical type. we could easily add that here.)
-if(vert_is_level(location)) then
- ! FIXME: if this is W, the top is nVertLevels+1
- if (lheight > nVertLevels) then
- istatus = 12
- return
- endif
- do j=1, num_kinds
- call triangle_interp(x, base_offset(j), tri_indices, weights, &
- nint(lheight), nVertLevels, interp_vals(j), ier)
- if(ier /= 0) then
- istatus = 13
- return
- endif
- enddo
- istatus = 0
- return
-endif
-
-! Vertical interpolation for pressure coordinates
- if(vert_is_pressure(location) ) then
- ! Need to get base offsets for the potential temperature, density, and water
- ! vapor mixing fields in the state vector
- call get_index_range(KIND_POTENTIAL_TEMPERATURE, pt_base_offset)
- call get_index_range(KIND_DENSITY, density_base_offset)
- call get_index_range(KIND_VAPOR_MIXING_RATIO, qv_base_offset)
-!print *, 'bases: t/rho/v = ', pt_base_offset, density_base_offset, qv_base_offset
- call find_pressure_bounds(x, lheight, tri_indices, weights, nVertLevels, &
- pt_base_offset, density_base_offset, qv_base_offset, lower, upper, fract, &
- ltemp, utemp, ier)
- if(ier /= 0) then
- interp_vals = MISSING_R8
- istatus = 17
- return
- endif
- ! if any of the input kinds are sensible temperature, we had to compute
- ! that value already to convert the pressure to the model levels, so just
- ! interpolate in the vertical and we're done with that one.
- do j=1, num_kinds
- if (obs_kinds(j) == KIND_TEMPERATURE) then
- ! Already have both values, interpolate in the vertical here.
- if (ltemp == MISSING_R8 .or. utemp == MISSING_R8) then
- istatus = 12
- return
- endif
- interp_vals(j) = (1.0_r8 - fract) * ltemp + fract * utemp
- else
- ! Next interpolate the observed quantity to the horizontal point at both levels
- call triangle_interp(x, base_offset(j), tri_indices, weights, lower, nVertLevels, &
- lower_interp, ier)
- if(ier /= 0) then
- istatus = 13
- return
- endif
- call triangle_interp(x, base_offset(j), tri_indices, weights, upper, nVertLevels, &
- upper_interp, ier)
- if(ier /= 0) then
- istatus = 13
- return
- endif
-
- ! Got both values, interpolate and return
- if (lower_interp == MISSING_R8 .or. upper_interp == MISSING_R8) then
- istatus = 12
- return
- endif
- interp_vals(j) = (1.0_r8 - fract) * lower_interp + fract * upper_interp
- endif
- enddo
- istatus = 0
- return
-endif
-
-
-! in this section, unlike the others, we loop 3 times adding successive
-! contributions to the return value. if there's an error the
-! result value has been set to 0 so we have to reset it to the
-! missing value instead of only setting the istatus code.
-if(vert_is_height(location)) then
- ! For height, can do simple vertical search for interpolation for now
- ! Get the lower and upper bounds and fraction for each column
- interp_vals(:) = 0.0_r8
- do i = 1, 3
- call find_height_bounds(lheight, nVertLevels, zGridCenter(:, tri_indices(i)), &
- lower, upper, fract, ier)
- if(ier /= 0) then
- istatus = 12
- interp_vals(:) = MISSING_R8
- return
- endif
- do j=1, num_kinds
- call vert_interp(x, base_offset(j), tri_indices(i), nVertLevels, lower, fract, v_interp, ier)
- if(ier /= 0) then
- istatus = 13
- interp_vals(:) = MISSING_R8
- return
- endif
- interp_vals(j) = interp_vals(j) + weights(i) * v_interp
- enddo
- enddo
- istatus = 0
- return
-endif
-
-! shouldn't get here.
-interp_vals(:) = MISSING_R8
-istatus = 101
-
-end subroutine local_interpolate
-
-
!------------------------------------------------------------------
function nc_write_model_atts( ncFileID ) result (ierr)
@@ -1175,6 +934,7 @@
integer :: nEdgesDimID, maxEdgesDimID
integer :: nVerticesDimID
integer :: VertexDegreeDimID
+integer :: nSoilLevelsDimID
integer :: nVertLevelsDimID
integer :: nVertLevelsP1DimID
@@ -1323,11 +1083,15 @@
len = nVertLevelsP1, dimid = NVertLevelsP1DimID),'nc_write_model_atts', &
'nVertLevelsP1 def_dim '//trim(filename))
+ call nc_check(nf90_def_dim(ncid=ncFileID, name='nSoilLevels', &
+ len = nSoilLevels, dimid = nSoilLevelsDimID),'nc_write_model_atts', &
+ 'nSoilLevels def_dim '//trim(filename))
+
!----------------------------------------------------------------------------
! Create the (empty) Coordinate Variables and the Attributes
!----------------------------------------------------------------------------
- ! Grid Longitudes
+ ! Cell Longitudes
call nc_check(nf90_def_var(ncFileID,name='lonCell', xtype=nf90_double, &
dimids=nCellsDimID, varid=VarID),&
'nc_write_model_atts', 'lonCell def_var '//trim(filename))
@@ -1338,7 +1102,7 @@
call nc_check(nf90_put_att(ncFileID, VarID, 'valid_range', (/ 0.0_r8, 360.0_r8 /)), &
'nc_write_model_atts', 'lonCell valid_range '//trim(filename))
- ! Grid Latitudes
+ ! Cell Latitudes
call nc_check(nf90_def_var(ncFileID,name='latCell', xtype=nf90_double, &
dimids=nCellsDimID, varid=VarID),&
'nc_write_model_atts', 'latCell def_var '//trim(filename))
@@ -1394,6 +1158,20 @@
call nc_check(nf90_put_att(ncFileID, VarID, 'long_name', 'vertex latitudes'), &
'nc_write_model_atts', 'latVertex long_name '//trim(filename))
+ ! Edge Longitudes
+ call nc_check(nf90_def_var(ncFileID,name='lonEdge', xtype=nf90_double, &
+ dimids=nEdgesDimID, varid=VarID),&
+ 'nc_write_model_atts', 'lonEdge def_var '//trim(filename))
+ call nc_check(nf90_put_att(ncFileID, VarID, 'long_name', 'edge longitudes'), &
+ 'nc_write_model_atts', 'lonEdge long_name '//trim(filename))
+
+ ! Edge Latitudes
+ call nc_check(nf90_def_var(ncFileID,name='latEdge', xtype=nf90_double, &
+ dimids=nEdgesDimID, varid=VarID),&
+ 'nc_write_model_atts', 'latEdge def_var '//trim(filename))
+ call nc_check(nf90_put_att(ncFileID, VarID, 'long_name', 'edge latitudes'), &
+ 'nc_write_model_atts', 'latEdge long_name '//trim(filename))
+
! Grid relationship information
call nc_check(nf90_def_var(ncFileID,name='nEdgesOnCell',xtype=nf90_int, &
dimids=nCellsDimID ,varid=VarID), &
@@ -1466,6 +1244,16 @@
call nc_check(nf90_put_var(ncFileID, VarID, latCell ), &
'nc_write_model_atts', 'latCell put_var '//trim(filename))
+ call nc_check(NF90_inq_varid(ncFileID, 'lonEdge', VarID), &
+ 'nc_write_model_atts', 'lonEdge inq_varid '//trim(filename))
+ call nc_check(nf90_put_var(ncFileID, VarID, lonEdge ), &
+ 'nc_write_model_atts', 'lonEdge put_var '//trim(filename))
+
+ call nc_check(NF90_inq_varid(ncFileID, 'latEdge', VarID), &
+ 'nc_write_model_atts', 'latEdge inq_varid '//trim(filename))
+ call nc_check(nf90_put_var(ncFileID, VarID, latEdge ), &
+ 'nc_write_model_atts', 'latEdge put_var '//trim(filename))
+
call nc_check(NF90_inq_varid(ncFileID, 'zgrid', VarID), &
'nc_write_model_atts', 'zgrid inq_varid '//trim(filename))
call nc_check(nf90_put_var(ncFileID, VarID, zGridFace ), &
@@ -1682,7 +1470,7 @@
where(dimIDs == TimeDimID) mystart = timeindex
where(dimIDs == TimeDimID) mycount = 1
- if ( debug > 1 ) then
+ if ( debug > 9 ) then
write(*,*)'nc_write_model_vars '//trim(varname)//' start is ',mystart(1:ncNdims)
write(*,*)'nc_write_model_vars '//trim(varname)//' count is ',mycount(1:ncNdims)
endif
@@ -1805,10 +1593,19 @@
real(r8), intent(in) :: filter_ens_mean(:)
+integer :: ivar
+
if ( .not. module_initialized ) call static_init_model
ens_mean = filter_ens_mean
+if (debug > 3) then
+ if (do_output()) print *, 'resetting ensemble mean: '
+ do ivar = 1, nfields
+ call dump_progvar(ivar, x=ens_mean, minmaxonly=.true.)
+ enddo
+endif
+
end subroutine ens_mean_for_model
@@ -1921,16 +1718,24 @@
subroutine get_close_obs(gc, base_obs_loc, base_obs_kind, &
obs_loc, obs_kind, num_close, close_ind, dist)
!
-! FIXME ... not tested.
-!
! Given a DART location (referred to as "base") and a set of candidate
! locations & kinds (obs, obs_kind), returns the subset close to the
! "base", their indices, and their distances to the "base" ...
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list