[Dart-dev] [6217] DART/branches/development/models: latest revisions of the mpas interface modules.
nancy at ucar.edu
nancy at ucar.edu
Fri May 31 16:30:23 MDT 2013
Revision: 6217
Author: nancy
Date: 2013-05-31 16:30:22 -0600 (Fri, 31 May 2013)
Log Message:
-----------
latest revisions of the mpas interface modules. new code from
soyoung that puts the clamping limits in the namelist, modeled
on code from wrf. several minor changes to both model_mods to
make them as consistent between the two as possible. added an
option to the two atm converters to print out the min/max data
values in each field as they do the conversions. this should
at some point be incorporated into the ocn versions of the
converters.
Modified Paths:
--------------
DART/branches/development/models/mpas_atm/dart_to_model.f90
DART/branches/development/models/mpas_atm/model_mod.f90
DART/branches/development/models/mpas_atm/model_to_dart.f90
DART/branches/development/models/mpas_atm/work/input.nml
DART/branches/development/models/mpas_ocn/model_mod.f90
-------------- next part --------------
Modified: DART/branches/development/models/mpas_atm/dart_to_model.f90
===================================================================
--- DART/branches/development/models/mpas_atm/dart_to_model.f90 2013-05-31 22:26:56 UTC (rev 6216)
+++ DART/branches/development/models/mpas_atm/dart_to_model.f90 2013-05-31 22:30:22 UTC (rev 6217)
@@ -35,7 +35,7 @@
get_time, get_date
use model_mod, only : static_init_model, statevector_to_analysis_file, &
get_model_size, get_model_analysis_filename, &
- write_model_time
+ write_model_time, print_variable_ranges
implicit none
@@ -52,10 +52,12 @@
character(len=256) :: dart_to_model_input_file = 'dart.ic'
logical :: advance_time_present = .false.
character(len=256) :: time_filename = 'mpas_time'
+logical :: print_data_ranges = .true.
namelist /dart_to_model_nml/ dart_to_model_input_file, &
- advance_time_present, &
- time_filename
+ advance_time_present, &
+ time_filename, &
+ print_data_ranges
!----------------------------------------------------------------------
@@ -104,14 +106,22 @@
endif
call close_restart(iunit)
-print *, 'read state vector'
!----------------------------------------------------------------------
+! if requested, print out the data ranges variable by variable
+! (note if we are clamping data values, that happens in the
+! conversion routine and these values are before the clamping happens.)
+!----------------------------------------------------------------------
+if (print_data_ranges) then
+ call print_variable_ranges(statevector)
+endif
+
+
+!----------------------------------------------------------------------
! update the current model state vector
! Convey the amount of time to integrate the model ...
! time_manager_nml: stop_option, stop_count increments
!----------------------------------------------------------------------
-print *, 'calling analysis file to state vector'
call statevector_to_analysis_file(statevector, model_analysis_filename, model_time)
! write time into in text format (YYYY-MM-DD_hh:mm:ss) into a file.
Modified: DART/branches/development/models/mpas_atm/model_mod.f90
===================================================================
--- DART/branches/development/models/mpas_atm/model_mod.f90 2013-05-31 22:26:56 UTC (rev 6216)
+++ DART/branches/development/models/mpas_atm/model_mod.f90 2013-05-31 22:30:22 UTC (rev 6217)
@@ -4,13 +4,7 @@
module model_mod
-! <next few lines under version control, do not edit>
-! $URL$
-! $Id$
-! $Revision$
-! $Date$
-
-! MPAS Atmosphere model interfact to the DART data assimilation system.
+! MPAS Atmosphere 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;
@@ -33,11 +27,11 @@
get_close_maxdist_init, get_close_type, &
set_location, get_location, horiz_dist_only, &
write_location, &
- vert_is_undef, VERTISUNDEF, &
- vert_is_surface, VERTISSURFACE, &
- vert_is_level, VERTISLEVEL, &
- vert_is_pressure, VERTISPRESSURE, &
- vert_is_height, VERTISHEIGHT, &
+ 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
@@ -65,7 +59,7 @@
KIND_U_WIND_COMPONENT, &
KIND_V_WIND_COMPONENT, &
KIND_PRESSURE, &
- KIND_DENSITY, &
+ KIND_DENSITY, &
KIND_VAPOR_MIXING_RATIO, &
KIND_SPECIFIC_HUMIDITY
@@ -82,7 +76,7 @@
! 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
+! 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
@@ -120,7 +114,8 @@
statevector_to_analysis_file, &
get_analysis_time, &
write_model_time, &
- get_grid_dims
+ get_grid_dims, &
+ print_variable_ranges
! version controlled file description for error handling, do not edit
@@ -158,26 +153,21 @@
type(xyz_get_close_type) :: cc_gc
type(xyz_location_type), allocatable :: cell_locs(:)
-! variables which are in the module namelist
+logical :: log_vert_interp = .true. ! if true, interpolate vertical in log space
+
+! 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
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'
-! 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.
@@ -190,7 +180,7 @@
! 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,
+! 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
@@ -216,11 +206,17 @@
! 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
+integer, parameter :: num_bounds_table_columns = 4
character(len=NF90_MAX_NAME) :: mpas_state_variables(max_state_variables * num_state_table_columns ) = ' '
+character(len=NF90_MAX_NAME) :: mpas_state_bounds(num_bounds_table_columns, max_state_variables ) = ' '
character(len=NF90_MAX_NAME) :: variable_table(max_state_variables, num_state_table_columns )
-namelist /mpas_vars_nml/ mpas_state_variables
+namelist /mpas_vars_nml/ mpas_state_variables, mpas_state_bounds
+! FIXME: this shouldn't be a global. the progvar array
+! should be allocated at run time and nfields should be part
+! of a larger derived type that includes nfields and an array
+! of progvartypes.
integer :: nfields
! Everything needed to describe a variable
@@ -232,7 +228,7 @@
character(len=NF90_MAX_NAME) :: units
character(len=NF90_MAX_NAME), dimension(NF90_MAX_VAR_DIMS) :: dimname
integer, dimension(NF90_MAX_VAR_DIMS) :: dimlens
- integer :: xtype ! netCDF variable type (NF90_double, etc.)
+ integer :: xtype ! netCDF variable type (NF90_double, etc.)
integer :: numdims ! number of dims - excluding TIME
integer :: numvertical ! number of vertical levels in variable
integer :: numcells ! number of horizontal locations (cell centers)
@@ -243,13 +239,14 @@
integer :: indexN ! location in dart state vector of last occurrence
integer :: dart_kind
character(len=paramname_length) :: kind_string
- logical :: clamping ! does variable need to be range-restricted before
+ logical :: clamping ! does variable need to be range-restricted before
real(r8) :: range(2) ! being stuffed back into MPAS analysis file.
+ logical :: out_of_range_fail ! is out of range fatal if range-checking?
end type progvartype
type(progvartype), dimension(max_state_variables) :: progvar
-! Grid parameters - the values will be read from an mpas analysis file.
+! Grid parameters - the values will be read from an mpas analysis file.
integer :: nCells = -1 ! Total number of cells making up the grid
integer :: nVertices = -1 ! Unique points in grid that are corners of cells
@@ -264,7 +261,7 @@
! FIXME: we read in a lot of metadata about the grids. if space becomes an
! issue we could consider reading in only the x,y,z arrays for all the items
-! plus the radius, and then compute the lat/lon for locations needed by
+! plus the radius, and then compute the lat/lon for locations needed by
! get_state_meta_data() on demand. most of the computations we need to do
! are actually easier in xyz coords (no issue with poles).
@@ -313,6 +310,19 @@
logical :: has_edge_u = .false. ! true = has original normal u on edges
logical :: has_uvreconstruct = .false. ! true = has reconstructed at centers
+! Do we have any state vector items located on the cell edges?
+! If not, avoid reading in or using the edge arrays to save space.
+! FIXME: this should be set after looking at the fields listed in the
+! namelist which are to be read into the state vector - if any of them
+! are located on the edges then this flag should be changed to .true.
+! however, the way the code is structured these arrays are allocated
+! before the details of field list is examined. since right now the
+! only possible field array that is on the edges is the 'u' edge normal
+! winds, search specifically for that in the state field list and set
+! this based on that. if any other data might be on edges, this routine
+! will need to be updated: is_edgedata_in_state_vector()
+logical :: data_on_edges = .false.
+
! 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.
@@ -361,7 +371,7 @@
! The triangle interpolation keeps a list of how many and which triangles
! overlap each regular lon-lat box. The number is stored in
! array triangle_num. The allocatable array
-! triangle_list lists the uniquen index
+! triangle_list lists the uniquen index
! of each overlapping triangle. The entry in
! triangle_start for a given regular lon-lat box indicates
! where the list of triangles begins in the triangle_list.
@@ -382,7 +392,7 @@
subroutine static_init_model()
! Called to do one time initialization of the model.
-!
+!
! All the grid information comes from the initialization of
! the dart_model_mod module.
@@ -441,61 +451,67 @@
!---------------------------------------------------------------
! 1) get grid dimensions
-! 2) allocate space for the grids
+! 2) allocate space for the grids
! 3) read them from the analysis file
! read_grid_dims() fills in the following module global variables:
! nCells, nVertices, nEdges, maxEdges, nVertLevels, nVertLevelsP1, vertexDegree, nSoilLevels
call read_grid_dims()
-allocate(latCell(nCells), lonCell(nCells))
+allocate(latCell(nCells), lonCell(nCells))
allocate(zGridFace(nVertLevelsP1, nCells))
allocate(zGridCenter(nVertLevels, nCells))
-!allocate(zEdgeFace( nVertLevelsP1, nEdges))
-!allocate(zEdgeCenter(nVertLevels, nEdges))
-allocate(zGridEdge(nVertLevels, nEdges))
-
allocate(cellsOnVertex(vertexDegree, nVertices))
allocate(nEdgesOnCell(nCells))
allocate(edgesOnCell(maxEdges, nCells))
allocate(cellsOnEdge(2, nEdges))
allocate(verticesOnCell(maxEdges, nCells))
allocate(edgeNormalVectors(3, nEdges))
-allocate(latEdge(nEdges), lonEdge(nEdges))
allocate(xVertex(nVertices), yVertex(nVertices), zVertex(nVertices))
-allocate(xEdge(nEdges), yEdge(nEdges), zEdge(nEdges))
+! see if U is in the state vector list. if not, don't read in or
+! use any of the Edge arrays to save space.
+data_on_edges = is_edgedata_in_state_vector(mpas_state_variables)
+
+if(data_on_edges) then
+ allocate(zGridEdge(nVertLevels, nEdges))
+ allocate(xEdge(nEdges), yEdge(nEdges), zEdge(nEdges))
+ allocate(latEdge(nEdges), lonEdge(nEdges))
+endif
+
! this reads in latCell, lonCell, zGridFace, cellsOnVertex
call get_grid()
! read in vert cell face locations and then compute vertical center locations
do kloc=1, nCells
- do iloc=1, nVertLevels
- zGridCenter(iloc,kloc) = (zGridFace(iloc,kloc) + zGridFace(iloc+1,kloc))*0.5_r8
- enddo
+ do iloc=1, nVertLevels
+ zGridCenter(iloc,kloc) = (zGridFace(iloc,kloc) + zGridFace(iloc+1,kloc))*0.5_r8
+ enddo
enddo
-! FIXME: This code is supposed to check whether an edge has 2 neighbours or 1 neighbour and then
-! compute the height accordingly. HOWEVER, the array cellsOnEdge does not change with
-! depth, but it should as an edge may have 2 neighbour cells at the top but not at depth.
-do kloc=1, nEdges
- do iloc=1, nVertLevels
- cel1 = cellsOnEdge(1,kloc)
- cel2 = cellsOnEdge(2,kloc)
- if (cel1>0 .and. cel2>0) then
- zGridEdge(iloc,kloc) = (zGridCenter(iloc,cel1) + zGridCenter(iloc,cel2))*0.5_r8
- else if (cel1>0) then
- zGridEdge(iloc,kloc) = zGridCenter(iloc,cel1)
- else if (cel2>0) then
- zGridEdge(iloc,kloc) = zGridCenter(iloc,cel2)
- else !this is bad...
- write(string1,*)'Edge ',kloc,' at vertlevel ',iloc,' has no neighbouring cells!'
- call error_handler(E_ERR,'static_init_model', string1, source, revision, revdate)
- endif
- enddo
-enddo
-
+if(data_on_edges) then
+ ! FIXME: This code is supposed to check whether an edge has 2 neighbours or 1 neighbour and then
+ ! compute the height accordingly. HOWEVER, the array cellsOnEdge does not change with
+ ! depth, but it should as an edge may have 2 neighbour cells at the top but not at depth.
+ do kloc=1, nEdges
+ do iloc=1, nVertLevels
+ cel1 = cellsOnEdge(1,kloc)
+ cel2 = cellsOnEdge(2,kloc)
+ if (cel1>0 .and. cel2>0) then
+ zGridEdge(iloc,kloc) = (zGridCenter(iloc,cel1) + zGridCenter(iloc,cel2))*0.5_r8
+ else if (cel1>0) then
+ zGridEdge(iloc,kloc) = zGridCenter(iloc,cel1)
+ else if (cel2>0) then
+ zGridEdge(iloc,kloc) = zGridCenter(iloc,cel2)
+ else !this is bad...
+ write(string1,*)'Edge ',kloc,' at vertlevel ',iloc,' has no neighbouring cells!'
+ call error_handler(E_ERR,'static_init_model', string1, source, revision, revdate)
+ endif
+ enddo
+ enddo
+endif
+
!---------------------------------------------------------------
! Compile the list of model variables to use in the creation
! of the DART state vector. Required to determine model_size.
@@ -553,7 +569,7 @@
call nc_check(nf90_inquire_variable(ncid, VarID, xtype=progvar(ivar)%xtype, &
dimids=dimIDs, ndims=numdims), 'static_init_model', 'inquire '//trim(string2))
- ! If the long_name and/or units attributes are set, get them.
+ ! If the long_name and/or units attributes are set, get them.
! They are not REQUIRED to exist but are nice to use if they are present.
if( nf90_inquire_attribute( ncid, VarID, 'long_name') == NF90_NOERR ) then
@@ -602,7 +618,8 @@
enddo DimensionLoop
- call set_variable_clamping(ivar)
+ ! this call sets: clamping, bounds, and out_of_range_fail in the progvar entry
+ call get_variable_bounds(mpas_state_bounds, ivar)
if (progvar(ivar)%numvertical == nVertLevels) then
progvar(ivar)%ZonHalf = .TRUE.
@@ -642,21 +659,21 @@
write( * , *)'static_init_model: grid is a global grid '
else
write(logfileunit, *)'static_init_model: grid has boundaries '
- write( * , *)'static_init_model: grid has boundaries '
+ write( * , *)'static_init_model: grid has boundaries '
endif
if ( all_levels_exist_everywhere ) then
write(logfileunit, *)'static_init_model: all cells have same number of vertical levels '
write( * , *)'static_init_model: all cells have same number of vertical levels '
else
- write(logfileunit, *)'static_init_model: cells have varying number of vertical levels '
+ write(logfileunit, *)'static_init_model: cells have varying number of vertical levels '
write( * , *)'static_init_model: cells have varying number of vertical levels '
endif
endif
! do some sanity checking here:
-! 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
+! 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)
@@ -735,7 +752,7 @@
integer, intent(in) :: index_in
type(location_type), intent(out) :: location
integer, optional, intent(out) :: var_type
-
+
! Local variables
integer :: nxp, nzp, iloc, vloc, nf, n
@@ -778,10 +795,15 @@
iloc = 1 + (myindx-1) / nzp ! cell index
vloc = myindx - (iloc-1)*nzp ! vertical level index
-! the zGrid array contains the location of the cell top and bottom faces, so it has one
+! the zGrid array contains the location of the cell top and bottom faces, so it has one
! more value than the number of cells in each column. for locations of cell centers
! you have to take the midpoint of the top and bottom face of the cell.
if (progvar(nf)%numedges /= MISSING_I) then
+ if (.not. data_on_edges) then
+ call error_handler(E_ERR, 'get_state_meta_data', &
+ 'Internal error: numedges present but data_on_edges false', &
+ source, revision, revdate, text2='variable '//trim(progvar(nf)%varname))
+ endif
if ( progvar(nf)%ZonHalf ) then
height = zGridEdge(vloc,iloc)
else
@@ -799,6 +821,11 @@
endif
if (progvar(nf)%numedges /= MISSING_I) then
+ if (.not. data_on_edges) then
+ call error_handler(E_ERR, 'get_state_meta_data', &
+ 'Internal error: numedges present but data_on_edges false', &
+ source, revision, revdate, text2='variable '//trim(progvar(nf)%varname))
+ endif
if (nzp <= 1) then
location = set_location(lonEdge(iloc),latEdge(iloc), height, VERTISSURFACE)
else
@@ -856,10 +883,10 @@
! ISTATUS = 12: Height vertical coordinate out of model range.
! ISTATUS = 13: Missing value in interpolation.
! ISTATUS = 16: Don't know how to do vertical velocity for now
-! ISTATUS = 17: Unable to compute pressure values
+! ISTATUS = 17: Unable to compute pressure values
! ISTATUS = 18: altitude illegal
! ISTATUS = 19: could not compute u using RBF code
-! ISTATUS = 101: Internal error; reached end of subroutine without
+! ISTATUS = 101: Internal error; reached end of subroutine without
! finding an applicable case.
!
@@ -917,7 +944,7 @@
! at the cell centers (U,V). there are namelist options to control
! which to use if both are in the state vector. we can compute
! specific humidity from the vapor mixing ratio (which we know
-! we have because we require potential temp, mixing ratio, and
+! we have because we require potential temp, mixing ratio, and
! density to be in the state vector in all cases.)
ivar = get_progvar_index_from_kind(obs_kind)
@@ -930,7 +957,7 @@
select case (obs_kind)
case (KIND_TEMPERATURE)
goodkind = .true.
- case (KIND_PRESSURE)
+ case (KIND_PRESSURE)
goodkind = .true.
case (KIND_SPECIFIC_HUMIDITY)
goodkind = .true.
@@ -980,7 +1007,7 @@
call compute_scalar_with_barycentric(x, location, 3, tvars, values, istatus)
if (debug > 4) print *, 'called compute temp, kind, val, istatus: ', obs_kind, interp_val, istatus
if (istatus /= 0) goto 100
-
+
! convert pot_temp, density, vapor mixing ratio into sensible temperature
interp_val = theta_to_tk(values(1), values(2), values(3))
@@ -1002,7 +1029,7 @@
! compute vapor pressure, then: sh = vp / (1.0 + vp)
tvars(1) = get_progvar_index_from_kind(KIND_VAPOR_MIXING_RATIO)
call compute_scalar_with_barycentric(x, location, 1, tvars, values, istatus)
- interp_val = values(1)
+ interp_val = values(1)
if (debug > 4) print *, 'called compute SH, kind, val, istatus: ', obs_kind, interp_val, istatus
if (istatus /= 0) goto 100
@@ -1015,7 +1042,7 @@
endif
if (debug > 4) print *, 'return val is: ', interp_val
-else
+else
! direct interpolation, kind is in the state vector
tvars(1) = ivar
@@ -1064,7 +1091,7 @@
!
! Typical sequence for adding new dimensions,variables,attributes:
! NF90_OPEN ! open existing netCDF dataset
-! NF90_redef ! put into define mode
+! NF90_redef ! put into define mode
! NF90_def_dim ! define additional dimensions (if any)
! NF90_def_var ! define variables: from name, type, and dims
! NF90_put_att ! assign attribute values
@@ -1105,7 +1132,7 @@
integer :: ivar, VarID, mpasFileID
!----------------------------------------------------------------------
-! local variables
+! local variables
!----------------------------------------------------------------------
! we are going to need these to record the creation date in the netCDF file.
@@ -1138,7 +1165,7 @@
write(filename,*) 'ncFileID', ncFileID
!-------------------------------------------------------------------------------
-! make sure ncFileID refers to an open netCDF file,
+! make sure ncFileID refers to an open netCDF file,
! and then put into define mode.
!-------------------------------------------------------------------------------
@@ -1148,7 +1175,7 @@
!-------------------------------------------------------------------------------
! We need the dimension ID for the number of copies/ensemble members, and
-! we might as well check to make sure that Time is the Unlimited dimension.
+! we might as well check to make sure that Time is the Unlimited dimension.
! Our job is create the 'model size' dimension.
!-------------------------------------------------------------------------------
@@ -1171,7 +1198,7 @@
dimid = StateVarDimID),'nc_write_model_atts', 'state def_dim '//trim(filename))
!-------------------------------------------------------------------------------
-! Write Global Attributes
+! Write Global Attributes
!-------------------------------------------------------------------------------
call DATE_AND_TIME(crdate,crtime,crzone,values)
@@ -1202,7 +1229,7 @@
! Create a variable for the state vector
!----------------------------------------------------------------------------
- ! Define the actual (3D) state vector, which gets filled as time goes on ...
+ ! Define the actual (3D) state vector, which gets filled as time goes on ...
call nc_check(nf90_def_var(ncid=ncFileID, name='state', xtype=nf90_real, &
dimids=(/StateVarDimID,MemberDimID,unlimitedDimID/),varid=StateVarID),&
'nc_write_model_atts','state def_var '//trim(filename))
@@ -1320,19 +1347,21 @@
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))
+ if(data_on_edges) then
+ ! 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))
+ ! 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))
+ endif
! Grid relationship information
call nc_check(nf90_def_var(ncFileID,name='nEdgesOnCell',xtype=nf90_int, &
@@ -1368,7 +1397,7 @@
! match shape of the variable to the dimension IDs
- call define_var_dims(ncFileID, ivar, MemberDimID, unlimitedDimID, myndims, mydimids)
+ call define_var_dims(ncFileID, ivar, MemberDimID, unlimitedDimID, myndims, mydimids)
! define the variable and set the attributes
@@ -1390,10 +1419,10 @@
! Finished with dimension/variable definitions, must end 'define' mode to fill.
!----------------------------------------------------------------------------
- call nc_check(nf90_enddef(ncfileID), 'prognostic enddef '//trim(filename))
+ call nc_check(nf90_enddef(ncFileID), 'prognostic enddef '//trim(filename))
!----------------------------------------------------------------------------
- ! Fill the coordinate variables that DART needs and has locally
+ ! Fill the coordinate variables that DART needs and has locally
!----------------------------------------------------------------------------
call nc_check(NF90_inq_varid(ncFileID, 'lonCell', VarID), &
@@ -1406,15 +1435,17 @@
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))
+ if(data_on_edges) then
+ 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, '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))
+ endif
call nc_check(NF90_inq_varid(ncFileID, 'zgrid', VarID), &
'nc_write_model_atts', 'zgrid inq_varid '//trim(filename))
@@ -1439,7 +1470,7 @@
!----------------------------------------------------------------------------
! Fill the coordinate variables needed for plotting only.
! DART has not read these in, so we have to read them from the input file
- ! and parrot them to the DART output file.
+ ! and parrot them to the DART output file.
!----------------------------------------------------------------------------
call nc_check(nf90_open(trim(grid_definition_filename), NF90_NOWRITE, mpasFileID), &
@@ -1492,7 +1523,7 @@
'nc_write_model_atts', 'lonVertex inq_varid '//trim(filename))
call nc_check(nf90_put_var(ncFileID, VarID, data1d ), &
'nc_write_model_atts', 'lonVertex put_var '//trim(filename))
-
+
call nc_check(nf90_inq_varid(mpasFileID, 'latVertex', VarID), &
'nc_write_model_atts', 'latVertex inq_varid ')
call nc_check(nf90_get_var(mpasFileID, VarID, data1d ), &
@@ -1518,7 +1549,7 @@
!------------------------------------------------------------------
-function nc_write_model_vars( ncFileID, state_vec, copyindex, timeindex ) result (ierr)
+function nc_write_model_vars( ncFileID, state_vec, copyindex, timeindex ) result (ierr)
! TJH 29 Aug 2011 -- all errors are fatal, so the
! return code is always '0 == normal', since the fatal errors stop execution.
@@ -1544,7 +1575,7 @@
integer :: ierr ! return value of function
integer, dimension(NF90_MAX_VAR_DIMS) :: dimIDs, mystart, mycount
-character(len=NF90_MAX_NAME) :: varname
+character(len=NF90_MAX_NAME) :: varname
integer :: i, ivar, VarID, ncNdims, dimlen
integer :: TimeDimID, CopyDimID
@@ -1568,7 +1599,7 @@
write(filename,*) 'ncFileID', ncFileID
!-------------------------------------------------------------------------------
-! make sure ncFileID refers to an open netCDF file,
+! make sure ncFileID refers to an open netCDF file,
!-------------------------------------------------------------------------------
call nc_check(nf90_inq_dimid(ncFileID, 'copy', dimid=CopyDimID), &
@@ -1590,7 +1621,7 @@
! We need to process the prognostic variables.
!----------------------------------------------------------------------------
- do ivar = 1,nfields
+ do ivar = 1,nfields
varname = trim(progvar(ivar)%varname)
string2 = trim(filename)//' '//trim(varname)
@@ -1632,7 +1663,7 @@
where(dimIDs == TimeDimID) mystart = timeindex
where(dimIDs == TimeDimID) mycount = 1
- if ( debug > 9 ) then
+ if ((debug > 9) .and. do_output()) 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
@@ -1717,7 +1748,7 @@
function get_model_size()
-! Returns the size of the model as an integer.
+! Returns the size of the model as an integer.
! Required for all applications.
integer :: get_model_size
@@ -1755,17 +1786,13 @@
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 ((debug > 3) .and. do_output()) then
if (do_output()) print *, 'resetting ensemble mean: '
- do ivar = 1, nfields
- call dump_progvar(ivar, x=ens_mean, minmaxonly=.true.)
- enddo
+ call print_variable_ranges(ens_mean)
endif
end subroutine ens_mean_for_model
@@ -1782,6 +1809,20 @@
if (allocated(zGridFace)) deallocate(zGridFace)
if (allocated(zGridCenter)) deallocate(zGridCenter)
if (allocated(cellsOnVertex)) deallocate(cellsOnVertex)
+if (allocated(nEdgesOnCell)) deallocate(nEdgesOnCell)
+if (allocated(edgesOnCell)) deallocate(edgesOnCell)
+if (allocated(cellsOnEdge)) deallocate(cellsOnEdge)
+if (allocated(verticesOnCell)) deallocate(verticesOnCell)
+if (allocated(edgeNormalVectors)) deallocate(edgeNormalVectors)
+if (allocated(xVertex)) deallocate(xVertex)
+if (allocated(yVertex)) deallocate(yVertex)
+if (allocated(zVertex)) deallocate(zVertex)
+if (allocated(zGridEdge)) deallocate(zGridEdge)
+if (allocated(xEdge)) deallocate(xEdge)
+if (allocated(yEdge)) deallocate(yEdge)
+if (allocated(zEdge)) deallocate(zEdge)
+if (allocated(latEdge)) deallocate(latEdge)
+if (allocated(lonEdge)) deallocate(lonEdge)
end subroutine end_model
@@ -1795,7 +1836,7 @@
! A model may choose to provide a NULL INTERFACE by returning
! .false. for the interf_provided argument. This indicates to
! the filter that if it needs to generate perturbed states, it
-! may do so by adding a perturbation to each model state
+! may do so by adding a perturbation to each model state
! variable independently. The interf_provided argument
! should be returned as .true. if the model wants to do its own
! perturbing of states.
@@ -1819,8 +1860,8 @@
! the subsequent code is a pert routine which
! can be used to add minor perturbations which can be spun up.
!
-! if all values in a field are identical (i.e. 0.0) this
-! routine will not change those values since it won't
+! if all values in a field are identical (i.e. 0.0) this
+! routine will not change those values since it won't
! make a new value outside the original min/max of that
! variable in the state vector. to handle this case you can
! remove the min/max limit lines below.
@@ -1885,10 +1926,10 @@
! "base", their indices, and their distances to the "base" ...
! For vertical distance computations, general philosophy is to convert all
-! vertical coordinates to a common coordinate. This coordinate type can be
+! vertical coordinates to a common coordinate. This coordinate type can be
! defined in the namelist with the variable "vert_localization_coord".
! But we first try a single coordinate type as the model level here.
-! FIXME: We need to add more options later.
+! FIXME: We need to add more options later.
! Vertical conversion is carried out by the subroutine vert_convert.
@@ -1977,8 +2018,8 @@
(istatus2 /= 0) ) then
dist(k) = 1.0e9_r8
else
- dist(k) = get_dist(base_obs_loc, local_obs_loc, base_obs_kind, obs_kind(t_ind))
-! if (debug > 4 .and. k < 100 .and. my_task_id() == 0) then
+ dist(k) = get_dist(base_obs_loc, local_obs_loc, base_obs_kind, obs_kind(t_ind))
+! if ((debug > 4) .and. (k < 100) .and. do_output()) then
! print *, 'calling get_dist'
! call write_location(0,base_obs_loc,charstring=string2)
! call error_handler(E_MSG, 'get_close_obs: base_obs_loc',string2,source, revision, revdate)
@@ -1992,7 +2033,7 @@
enddo
endif
-if (debug > 2) then
+if ((debug > 2) .and. do_output()) then
call write_location(0,base_obs_loc,charstring=string2)
print *, 'get_close_obs: nclose, base_obs_loc ', num_close, trim(string2)
endif
@@ -2004,9 +2045,9 @@
subroutine init_time(time)
-! Companion interface to init_conditions. Returns a time that is somehow
+! Companion interface to init_conditions. Returns a time that is somehow
! appropriate for starting up a long integration of the model.
-! At present, this is only used if the namelist parameter
+! At present, this is only used if the namelist parameter
! start_from_restart is set to .false. in the program perfect_model_obs.
type(time_type), intent(out) :: time
@@ -2028,13 +2069,13 @@
! Returns a model state vector, x, that is some sort of appropriate
! initial condition for starting up a long integration of the model.
-! At present, this is only used if the namelist parameter
+! At present, this is only used if the namelist parameter
! start_from_restart is set to .false. in the program perfect_model_obs.
real(r8), intent(out) :: x(:)
if ( .not. module_initialized ) call static_init_model
-
+
! this shuts up the compiler warnings about unused variables
x = 0.0_r8
@@ -2051,13 +2092,13 @@
! Does a single timestep advance of the model. The input value of
! the vector x is the starting condition and x is updated to reflect
! the changed state after a timestep. The time argument is intent
-! in and is used for models that need to know the date/time to
+! in and is used for models that need to know the date/time to
! compute a timestep, for instance for radiation computations.
! This interface is only called IF the namelist parameter
-! async is set to 0 in perfect_model_obs or filter -OR- if the
+! async is set to 0 in perfect_model_obs or filter -OR- if the
! program integrate_model is to be used to advance the model
! state as a separate executable. If none of these options
-! are used (the model will only be advanced as a separate
+! are used (the model will only be advanced as a separate
! model-specific executable), this can be a NULL INTERFACE.
real(r8), intent(inout) :: x(:)
@@ -2211,7 +2252,7 @@
where(dimIDs == TimeDimID) mystart = TimeDimLength ! pick the latest time
where(dimIDs == TimeDimID) mycount = 1 ! only use one time
- if ( debug > 9 ) then
+ if ((debug > 10) .and. do_output()) then
write(*,*)'analysis_file_to_statevector '//trim(varname)//' start = ',mystart(1:ncNdims)
write(*,*)'analysis_file_to_statevector '//trim(varname)//' count = ',mycount(1:ncNdims)
endif
@@ -2359,7 +2400,7 @@
string2 = trim(filename)//' '//trim(varname)
if (( varname == 'uReconstructZonal' .or. &
- varname == 'uReconstructMeridional' ) .and. update_u_from_reconstruct ) then
+ varname == 'uReconstructMeridional' ) .and. update_u_from_reconstruct ) then
if (done_winds) cycle PROGVARLOOP
! this routine updates the edge winds from both the zonal and meridional
@@ -2409,7 +2450,7 @@
where(dimIDs == TimeDimID) mystart = TimeDimLength
where(dimIDs == TimeDimID) mycount = 1 ! only the latest one
- if ( debug > 9 ) then
+ if ((debug > 9) .and. do_output()) then
write(*,*)'statevector_to_analysis_file '//trim(varname)//' start is ',mystart(1:ncNdims)
write(*,*)'statevector_to_analysis_file '//trim(varname)//' count is ',mycount(1:ncNdims)
endif
@@ -2419,17 +2460,16 @@
allocate(data_1d_array(mycount(1)))
call vector_to_prog_var(state_vector, ivar, data_1d_array)
- write(string1, '(A,A32,2F16.7)') 'data min/max ', trim(varname), minval(data_1d_array), maxval(data_1d_array)
+ write(string1, '(A,A32,2(1x,E22.14))') 'min/max ', trim(varname), &
+ minval(data_1d_array), maxval(data_1d_array)
call error_handler(E_MSG, '', string1, source,revision,revdate)
+ ! did the user specify lower and/or upper bounds for this variable?
+ ! if so, follow the instructions to either fail on out-of-range values,
+ ! or set out-of-range values to the given min or max vals
if ( progvar(ivar)%clamping ) then
- where ( data_1d_array < progvar(ivar)%range(1) ) data_1d_array = progvar(ivar)%range(1)
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list