[Dart-dev] [5940] DART/branches/development/models/mpas_ocn: fixed up the state vector a bit and added code to read
nancy at ucar.edu
nancy at ucar.edu
Wed Dec 5 16:10:20 MST 2012
Revision: 5940
Author: nancy
Date: 2012-12-05 16:10:19 -0700 (Wed, 05 Dec 2012)
Log Message:
-----------
fixed up the state vector a bit and added code to read
in the zMid() array for now. also committed a sample
obs_seq.in with a single fake argo observation for testing.
Modified Paths:
--------------
DART/branches/development/models/mpas_ocn/model_mod.f90
DART/branches/development/models/mpas_ocn/work/input.nml
Added Paths:
-----------
DART/branches/development/models/mpas_ocn/work/obs_seq.in
-------------- next part --------------
Modified: DART/branches/development/models/mpas_ocn/model_mod.f90
===================================================================
--- DART/branches/development/models/mpas_ocn/model_mod.f90 2012-12-05 21:30:13 UTC (rev 5939)
+++ DART/branches/development/models/mpas_ocn/model_mod.f90 2012-12-05 23:10:19 UTC (rev 5940)
@@ -292,6 +292,8 @@
real(r8), allocatable :: zGridCenter(:,:) ! geometric depth at cell centers (nVertLevels, nCells)
real(r8), allocatable :: zGridEdge(:,:) ! geometric depth at edge centers (nVertLevels, nEdges)
+real(r8), allocatable :: zMid(:,:) ! depths at midpoints - may be able to be computed instead
+ ! of requiring it in the input file (save file space). FIXME
real(r8), allocatable :: hZLevel(:) ! layer thicknesses - maybe - FIXME
!real(r8), allocatable :: zEdgeCenter(:,:) ! geometric height at edges faces (nVertLevels ,nEdges)
@@ -460,7 +462,7 @@
allocate(zGridFace(nVertLevelsP1, nCells))
allocate(zGridCenter(nVertLevels, nCells))
allocate(zGridEdge(nVertLevels, nEdges))
-allocate(hZLevel( nVertLevels))
+allocate(hZLevel(nVertLevels), zMid(nVertLevels, nCells))
!allocate(zEdgeCenter(nVertLevels, nEdges))
allocate(cellsOnVertex(vertexDegree, nVertices))
@@ -473,7 +475,7 @@
allocate(xVertex(nVertices), yVertex(nVertices), zVertex(nVertices))
allocate(xEdge(nEdges), yEdge(nEdges), zEdge(nEdges))
-! this reads in latCell, lonCell, hZLevel, cellsOnVertex
+! this reads in latCell, lonCell, hZLevel, zMid, cellsOnVertex
call get_grid()
! FIXME: This code is supposed to check whether an edge has 2 neighbours or 1 neighbour and then
@@ -903,6 +905,8 @@
! exceptions if the kind isn't directly
! a field in the state vector:
select case (obs_kind)
+ !case (KIND_TEMPERATURE) ! potential temperature is in state vector not in-situ temp
+ ! goodkind = .true.
!case (KIND_PRESSURE)
! goodkind = .true.
end select
@@ -1330,6 +1334,11 @@
call nc_check(nf90_put_var(ncFileID, VarID, hZLevel ), &
'nc_write_model_atts', 'hZLevel put_var '//trim(filename))
+ call nc_check(NF90_inq_varid(ncFileID, 'zMid', VarID), &
+ 'nc_write_model_atts', 'zMid inq_varid '//trim(filename))
+ call nc_check(nf90_put_var(ncFileID, VarID, zMid ), &
+ 'nc_write_model_atts', 'zMid put_var '//trim(filename))
+
call nc_check(NF90_inq_varid(ncFileID, 'nEdgesOnCell', VarID), &
'nc_write_model_atts', 'nEdgesOnCell inq_varid '//trim(filename))
call nc_check(nf90_put_var(ncFileID, VarID, nEdgesOnCell ), &
@@ -2809,6 +2818,11 @@
call nc_check(nf90_get_var( ncid, VarID, hZLevel), &
'get_grid', 'get_var hZLevel '//trim(grid_definition_filename))
+call nc_check(nf90_inq_varid(ncid, 'zMid', VarID), &
+ 'get_grid', 'inq_varid zMid '//trim(grid_definition_filename))
+call nc_check(nf90_get_var( ncid, VarID, zMid), &
+ 'get_grid', 'get_var zMid '//trim(grid_definition_filename))
+
call nc_check(nf90_inq_varid(ncid, 'cellsOnVertex', VarID), &
'get_grid', 'inq_varid cellsOnVertex '//trim(grid_definition_filename))
call nc_check(nf90_get_var( ncid, VarID, cellsOnVertex), &
@@ -2931,7 +2945,7 @@
! FIXME : Use layer thicknesses (perhaps hZLevel) to determine all layer information
zGridFace(:,:) = MISSING_R8
-zGridCenter(:,:) = MISSING_R8
+zGridCenter(:,:) = -1.0_r8 * zMid(:,:) ! all depths are negative in the input file; obs are +
zGridEdge(:,:) = MISSING_R8
! A little sanity check
@@ -4153,7 +4167,7 @@
!------------------------------------------------------------------
-subroutine find_height_bounds(depth, nbounds, bounds, lower, upper, fract, ier)
+subroutine find_depth_bounds(depth, nbounds, bounds, lower, upper, fract, ier)
! Finds position of a given height in an array of height grid points and returns
! the index of the lower and upper bounds and the fractional offset. ier
@@ -4179,10 +4193,21 @@
lower = -1
upper = -1
+if (debug > 7) then
+ print *, 'ready to check height bounds'
+ print *, 'array ranges from 1 to ', nbounds
+ print *, 'depth to check is ' , depth
+ print *, 'array(1) = ', bounds(1)
+ print *, 'array(nbounds) = ', bounds(nbounds)
+endif
+
! Bounds check
if(depth < bounds(1)) ier = 998
if(depth > bounds(nbounds)) ier = 9998
-if (ier /= 0) return
+if (ier /= 0) then
+ if (debug > 7) print *, 'could not find vertical location'
+ return
+endif
! Search two adjacent layers that enclose the given point
do i = 2, nbounds
@@ -4190,6 +4215,7 @@
lower = i - 1
upper = i
fract = (depth - bounds(lower)) / (bounds(upper) - bounds(lower))
+ if (debug > 7) print *, 'found it. lower, upper, fract = ', lower, upper, fract
return
endif
end do
@@ -4198,8 +4224,9 @@
! are tested for at the start of this routine. if you get here
! there is a coding error.
ier = 3
+if (debug > 7) print *, 'internal code inconsistency: could not find vertical location'
-end subroutine find_height_bounds
+end subroutine find_depth_bounds
!------------------------------------------------------------------
@@ -4332,11 +4359,11 @@
! Get the lower and upper bounds and fraction for each column
do i=1, nc
if (oncenters) then
- call find_height_bounds(vert, nVertLevels, zGridCenter(:, ids(i)), &
+ call find_depth_bounds(vert, nVertLevels, zGridCenter(:, ids(i)), &
lower(i), upper(i), fract(i), ier)
else
- call find_height_bounds(vert, nVertLevels, zGridEdge(:, ids(i)), &
+ call find_depth_bounds(vert, nVertLevels, zGridEdge(:, ids(i)), &
lower(i), upper(i), fract(i), ier)
endif
if (ier /= 0) return
Modified: DART/branches/development/models/mpas_ocn/work/input.nml
===================================================================
--- DART/branches/development/models/mpas_ocn/work/input.nml 2012-12-05 21:30:13 UTC (rev 5939)
+++ DART/branches/development/models/mpas_ocn/work/input.nml 2012-12-05 23:10:19 UTC (rev 5940)
@@ -11,7 +11,7 @@
output_interval = 1,
restart_in_file_name = "perfect_ics",
restart_out_file_name = "perfect_restart",
- obs_seq_in_file_name = "obs_seq.1800obs",
+ obs_seq_in_file_name = "obs_seq.in",
obs_seq_out_file_name = "obs_seq.out",
adv_ens_command = "./advance_model.csh",
output_timestamps = .false.,
@@ -104,8 +104,11 @@
nlat = 36,
output_box_info = .false.,
print_box_level = 0,
- /
+ /
+&xyz_location_nml
+ /
+
&cov_cutoff_nml
select_localization = 1
/
@@ -116,10 +119,61 @@
save_reg_diagnostics = .false.,
reg_diagnostics_file = "reg_diagnostics"
/
+
&obs_sequence_nml
write_binary_obs_sequence = .false.
/
+# options for assimilate and evaluate include:
+# (first column is the obs, second is what needs to be in the state vector
+# to compute it.)
+# ARGO_U_CURRENT_COMPONENT, KIND_U_CURRENT_COMPONENT
+# ARGO_V_CURRENT_COMPONENT, KIND_V_CURRENT_COMPONENT
+# ARGO_SALINITY, KIND_SALINITY
+# ARGO_TEMPERATURE, KIND_TEMPERATURE
+# ADCP_U_CURRENT_COMPONENT, KIND_U_CURRENT_COMPONENT
+# ADCP_V_CURRENT_COMPONENT, KIND_V_CURRENT_COMPONENT
+# ADCP_SALINITY, KIND_SALINITY
+# ADCP_TEMPERATURE, KIND_TEMPERATURE
+# FLOAT_SALINITY, KIND_SALINITY
+# FLOAT_TEMPERATURE, KIND_TEMPERATURE
+# DRIFTER_U_CURRENT_COMPONENT, KIND_U_CURRENT_COMPONENT
+# DRIFTER_V_CURRENT_COMPONENT, KIND_V_CURRENT_COMPONENT
+# DRIFTER_SALINITY, KIND_SALINITY
+# DRIFTER_TEMPERATURE, KIND_TEMPERATURE
+# GLIDER_U_CURRENT_COMPONENT, KIND_U_CURRENT_COMPONENT
+# GLIDER_V_CURRENT_COMPONENT, KIND_V_CURRENT_COMPONENT
+# GLIDER_SALINITY, KIND_SALINITY
+# GLIDER_TEMPERATURE, KIND_TEMPERATURE
+# MOORING_U_CURRENT_COMPONENT, KIND_U_CURRENT_COMPONENT
+# MOORING_V_CURRENT_COMPONENT, KIND_V_CURRENT_COMPONENT
+# MOORING_SALINITY, KIND_SALINITY
+# MOORING_TEMPERATURE, KIND_TEMPERATURE
+# MOORING_PRESSURE, KIND_PRESSURE
+# BOTTLE_SALINITY, KIND_SALINITY
+# BOTTLE_TEMPERATURE, KIND_TEMPERATURE
+# CTD_SALINITY, KIND_SALINITY
+# CTD_TEMPERATURE, KIND_TEMPERATURE
+# TCTD_SALINITY, KIND_SALINITY
+# TCTD_TEMPERATURE, KIND_TEMPERATURE
+# STD_SALINITY, KIND_SALINITY
+# STD_TEMPERATURE, KIND_TEMPERATURE
+# XCTD_SALINITY, KIND_SALINITY
+# XCTD_TEMPERATURE, KIND_TEMPERATURE
+# MBT_SALINITY, KIND_SALINITY
+# MBT_TEMPERATURE, KIND_TEMPERATURE
+# XBT_SALINITY, KIND_SALINITY
+# XBT_TEMPERATURE, KIND_TEMPERATURE
+# DBT_SALINITY, KIND_SALINITY
+# DBT_TEMPERATURE, KIND_TEMPERATURE
+# APB_SALINITY, KIND_SALINITY
+# APB_TEMPERATURE, KIND_TEMPERATURE
+# DOPPLER_U_CURRENT_COMPONENT, KIND_U_CURRENT_COMPONENT
+# DOPPLER_V_CURRENT_COMPONENT, KIND_V_CURRENT_COMPONENT
+# DOPPLER_W_CURRENT_COMPONENT, KIND_W_CURRENT_COMPONENT
+# SATELLITE_MICROWAVE_SST, KIND_TEMPERATURE
+# SATELLITE_INFRARED_SST, KIND_TEMPERATURE
+
&obs_kind_nml
assimilate_these_obs_types = 'null',
evaluate_these_obs_types = 'ARGO_TEMPERATURE'
@@ -131,8 +185,8 @@
/
&model_nml
- model_analysis_filename = '../data/mpas_out.nc'
- grid_definition_filename = '../data/mpas_out.nc'
+ model_analysis_filename = '../data/mpas_ocean_in.nc'
+ grid_definition_filename = '../data/mpas_ocean_in.nc'
assimilation_period_days = 0,
assimilation_period_seconds = 3600,
output_state_vector = .true.,
@@ -146,15 +200,23 @@
debug = 2
/
-# NOTE: since h(nVertLevels,nCells,Time) ... it is not really KIND_SEA_SURFACE_HEIGHT
+# NOTE: h(nVertLevels,nCells,Time) is layer thickness, not SSH
+# we don't have a kind defined for thickness yet.
+# examples of other kinds which could be in state vector.
+# the temperature should really be potential temp, but the
+# conversion routines (at the end of the model_mod.f90 file)
+# are commented out here because they're for the atmosphere not ocean.
+# we need to add a case for temp vs potential temp, and conversion code.
+# 'tracer1', 'KIND_TRACER_CONCENTRATION'
+# 'u', 'KIND_EDGE_NORMAL_SPEED',
+# 'temperature', 'KIND_POTENTIAL_TEMPERATURE',
&mpas_vars_nml
- mpas_state_variables = 'temperature', 'KIND_TEMPERATURE',
- 'salinity', 'KIND_SALINITY',
- 'rho', 'KIND_DENSITY',
- 'u', 'KIND_EDGE_NORMAL_SPEED',
- 'h', 'KIND_SEA_SURFACE_HEIGHT'
- 'tracer1', 'KIND_TRACER_CONCENTRATION'
+ mpas_state_variables = 'temperature', 'KIND_TEMPERATURE',
+ 'salinity', 'KIND_SALINITY',
+ 'rho', 'KIND_DENSITY',
+ 'uReconstructMeridional', 'KIND_U_CURRENT_COMPONENT',
+ 'uReconstructZonal', 'KIND_V_CURRENT_COMPONENT',
/
&model_to_dart_nml
Added: DART/branches/development/models/mpas_ocn/work/obs_seq.in
===================================================================
--- DART/branches/development/models/mpas_ocn/work/obs_seq.in (rev 0)
+++ DART/branches/development/models/mpas_ocn/work/obs_seq.in 2012-12-05 23:10:19 UTC (rev 5940)
@@ -0,0 +1,16 @@
+ obs_sequence
+obs_kind_definitions
+ 1
+ 10 ARGO_TEMPERATURE
+ num_copies: 0 num_qc: 0
+ num_obs: 1 max_num_obs: 1
+ first: 1 last: 1
+ OBS 1
+ -1 -1 -1
+obdef
+loc3d
+ 3.839724354387525 0.000000000000000 20.0000000000000 3
+kind
+ 10
+ 0 1
+ 1.00000000000000
Property changes on: DART/branches/development/models/mpas_ocn/work/obs_seq.in
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Added: svn:eol-style
+ native
More information about the Dart-dev
mailing list