[Dart-dev] [4293] DART/trunk: Small functional change for all the location_mods, but large amounts of code
nancy at ucar.edu
nancy at ucar.edu
Fri Feb 26 13:07:23 MST 2010
Revision: 4293
Author: nancy
Date: 2010-02-26 13:07:23 -0700 (Fri, 26 Feb 2010)
Log Message:
-----------
Small functional change for all the location_mods, but large amounts of code
reformatting, bringing infrequently used location types up to date
with current code, making versions internally consistent with each
other, and minor bug fixes in little-used types of locations. This is a
combination of work from nancy and Tim. Details follow (and sorry they
got pretty long, but i believe they're accurate).
- added the option to format a location into a human-friendly ascii format,
independent of which location module is being linked with. useful for code
which needs to stay location-independent, such as generic tools.
- removed 'simple_threed_sphere' which was only different from the full
'threed_sphere' in that the vertical could only be one type. added a
'twod' type, which is a periodic (x,y) location in the domain [0-1] in
both dimensions.
- added netcdf write routines which have the same calling interface for
all types of locations. opens the right dimensionality arrays, and
writes into them with an appropriate set of attributes. useful for code
that wants to store a list of locations in a netcdf file but needs to
work with multiple location types (like generic tools).
- renamed the internal location derived type variables in the annulus
to be azm and rad (azimuth and radius), instead of lat and lon. it makes
the code clearer as to what is being used for what. also added a missing
namelist that the comments indicated should be there, for min and max
values of the radius, and added the code to enforce those limits.
- changed the default format for writing ascii locations to be 16 decimal
digits instead of 14. this is intended to avoid problems with roundoff
error. note that this will change, very slightly, the numeric values
if you read/write locations in ascii -- but in my tests they now match
more closely the actual machine values that you'd get if you read/wrote
in binary format.
- made all the code call a consistent routine in the utilities module
to decide if the intent is to write locations in ascii or binary.
the default code is now in a single place instead of being implemented
in every different location module. uses the 'ascii_file_format()'
function from the utilities_mod.
- the interactive_locations subroutine no longer suggests hectopascals
as the units for the vertical if you give a vertical type of 'surface'.
many users create surface observations with a numerical value of meters,
so they can compute differences between observations which should be
located at the surface and where the model thinks the surface is.
users can continue to input pressure if they want. they will need to
input it as pascals instead of hectopascals, and the prompting text
will no longer suggest which units should be used.
- any generic utility can call the 'vert_is_xxx()' functions and
even if the locations module has no vertical component it will return
false, so utility code can be compiled independently of the locations
module it is linked with.
- the order of the subroutines and functions in the code has been
made consistent between each location module, so it's easier to
use diff and xdiff between different types of locations modules.
now for the most part only actual code differences will be found,
instead of comments, white space, and order of code in the source.
- added an overloaded set_location() routine that takes a fortran
array of values, to set a single location in an location-independent
form. the code checks for the proper number of input arguments in
the list, and then calls the standard set_location() to do the actual
work. this allows generic tools to pass through things from a namelist
which are appropriate for different location types without needing
special case code.
- the 2 kind arguments to get_close() are marked as optional, so the
routine can be called with them but they aren't required inside the
location modules - none of the standard code uses the kinds. user code
can intercept this call and is welcome to require and use the kinds,
but the location code in the location_mod does not need them.
- added code to make the 'is_location_in_region()' function do an
accurate computation for annulus, column, and the twod_sphere versions.
- removed obsolete 'alloc_get_close_type()' functions in some of the
little-used location mod types, and added versions of the current
collection of get_close_type functions.
- made the get_close_type private, since no other code outside the
location_mod should be able to access the contents of the derived type.
- added test directories and test cases for each of the location types.
the 'location_test.f90' source is at the top level directory, and the
'testall.csh' script will run the tests for each of the locations types.
Modified Paths:
--------------
DART/trunk/location/annulus/README
DART/trunk/location/annulus/location_mod.f90
DART/trunk/location/column/location_mod.f90
DART/trunk/location/oned/location_mod.f90
DART/trunk/location/threed_sphere/location_mod.f90
DART/trunk/location/threed_sphere/location_mod.html
DART/trunk/location/twod_sphere/location_mod.f90
DART/trunk/models/MITgcm_annulus/work/input.nml
DART/trunk/obs_sequence/obs_sequence_tool.f90
Added Paths:
-----------
DART/trunk/location/README
DART/trunk/location/annulus/test/
DART/trunk/location/annulus/test/input.nml
DART/trunk/location/annulus/test/mkmf_location_test
DART/trunk/location/annulus/test/path_names_location_test
DART/trunk/location/annulus/test/test.in
DART/trunk/location/column/test/
DART/trunk/location/column/test/input.nml
DART/trunk/location/column/test/mkmf_location_test
DART/trunk/location/column/test/path_names_location_test
DART/trunk/location/column/test/test.in
DART/trunk/location/location_test.f90
DART/trunk/location/location_test2.f90
DART/trunk/location/oned/test/
DART/trunk/location/oned/test/input.nml
DART/trunk/location/oned/test/mkmf_location2_test
DART/trunk/location/oned/test/mkmf_location_test
DART/trunk/location/oned/test/path_names_location2_test
DART/trunk/location/oned/test/path_names_location_test
DART/trunk/location/oned/test/test.in
DART/trunk/location/testall.csh
DART/trunk/location/threed_sphere/test/
DART/trunk/location/threed_sphere/test/input.nml
DART/trunk/location/threed_sphere/test/mkmf_location_test
DART/trunk/location/threed_sphere/test/path_names_location_test
DART/trunk/location/threed_sphere/test/test.in
DART/trunk/location/twod/
DART/trunk/location/twod/location_mod.f90
DART/trunk/location/twod/test/
DART/trunk/location/twod/test/input.nml
DART/trunk/location/twod/test/mkmf_location_test
DART/trunk/location/twod/test/path_names_location_test
DART/trunk/location/twod/test/test.in
DART/trunk/location/twod_sphere/test/
DART/trunk/location/twod_sphere/test/input.nml
DART/trunk/location/twod_sphere/test/mkmf_location_test
DART/trunk/location/twod_sphere/test/path_names_location_test
DART/trunk/location/twod_sphere/test/test.in
Removed Paths:
-------------
DART/trunk/location/oned/location_test.f90
DART/trunk/location/simple_threed_sphere/
DART/trunk/location/twod_sphere/location_test.f90
-------------- next part --------------
Added: DART/trunk/location/README
===================================================================
--- DART/trunk/location/README (rev 0)
+++ DART/trunk/location/README 2010-02-26 20:07:23 UTC (rev 4293)
@@ -0,0 +1,52 @@
+# DART software - Copyright © 2004 - 2010 UCAR. This open source software is
+# provided by UCAR, "as is", without charge, subject to all terms of use at
+# http://www.image.ucar.edu/DAReS/DART/DART_download
+#
+# DART $Id: README 4267 2010-02-08 19:13:26Z thoar $
+
+
+The core algorithms of DART work with many different models
+which have a variety of coordinate systems. This directory
+provides code for creating, setting/getting, copying location
+information (coordinates) independently of the actual specific
+coordinate information. It also contains distance routines
+needed by the DART algorithms.
+
+Each of the different location_mod.f90 files provides the same
+set of interfaces and defines a 'module location_mod', so by
+selecting the proper version in your path_names_xxx file you
+can compile your model code with the main DART routines.
+
+threed_sphere:
+
+The most frequently used version for real-world 3d models is
+the 'threed_sphere' module, which uses latitude and longitude
+coordinates, plus a vertical coordinate which can be meters,
+pressure, or model level.
+
+oned:
+
+For small models, the 'oned' version is most frequently used.
+It has a cyclic domain from 0 to 1.
+
+others:
+
+For other uses, there are also:
+ 'column', no x,y but 1d height, pressure, or model level for vertical.
+ 'annulus', a hollow 3d cylinder with azimuth, radius, and depth.
+ 'twod', a periodic 2d domain with x,y coordinates between 0 and 1.
+ 'twod_sphere', a 2d shell with latitude, longitude pairs.
+
+Other schemes can be added, as needed by the models.
+Possible ideas are a non-periodic version of the 1d, 2d
+cartesian versions, or a 3d cartesian - either periodic or
+not. Email 'dart at ucar.edu' if you have a different scheme
+which we might want to support.
+
+testing:
+
+The testall.csh script will cd into each of the locations
+module directories and build and run a simple testcase to
+verify the distance calculations and the read/write options
+are working.
+
Modified: DART/trunk/location/annulus/README
===================================================================
--- DART/trunk/location/annulus/README 2010-02-25 04:41:33 UTC (rev 4292)
+++ DART/trunk/location/annulus/README 2010-02-26 20:07:23 UTC (rev 4293)
@@ -4,6 +4,15 @@
#
# DART $Id$
+Last modified 29 Jan, 2010 by nancy collins
+
+made consistent with other location modules, including the
+revised get_close() routines. this code didn't have any
+min/max limits on the radius, and was still using lon/lat
+in a few subroutines. i put in a stub of a namelist for
+the radius limits, and changed the remaining lon/lats
+to azm and rad.
+
Last modified June 28, 2004 by Jim Hansen.
The location_mod.f90 for the annulus is a modified version
Modified: DART/trunk/location/annulus/location_mod.f90
===================================================================
--- DART/trunk/location/annulus/location_mod.f90 2010-02-25 04:41:33 UTC (rev 4292)
+++ DART/trunk/location/annulus/location_mod.f90 2010-02-26 20:07:23 UTC (rev 4293)
@@ -17,27 +17,27 @@
! direction (longitude-like). The radial distance is latitude-like,
! and the vertical coordinate is zero at the bottom of the annulus.
!
-! Note: in an effort to maintain consistency with the other
-! location_mods, will use lon for azimuthal direction, lat for
-! radial direction, and lev for depth.
-! (Note2: not clear this is a good thing - what consistency do we buy?
-! might want to use real names for the values for clarity.)
use types_mod, only : r8, PI, RAD2DEG, DEG2RAD, MISSING_R8, MISSING_I
-use utilities_mod, only : register_module, error_handler, E_ERR, ascii_file_format
+use utilities_mod, only : register_module, error_handler, E_ERR, ascii_file_format, &
+ nc_check, find_namelist_in_file, check_namelist_read, &
+ do_output, do_nml_file, do_nml_term, nmlfileunit, &
+ open_file, close_file, nc_check, is_longitude_between
use random_seq_mod, only : random_seq_type, init_random_seq, random_uniform
implicit none
private
-public :: location_type, get_dist, get_location, set_location, &
- set_location2, set_location_missing, is_location_in_region, &
- write_location, read_location, interactive_location, &
- vert_is_pressure, vert_is_level, vert_is_height, query_location, &
- LocationDims, LocationName, LocationLName, &
- get_close_obs, alloc_get_close_obs, get_close_type, &
- get_close_maxdist_init, get_close_obs_init, get_close_obs_destroy, &
- operator(==), operator(/=)
+public :: location_type, get_location, set_location, &
+ set_location_missing, is_location_in_region, &
+ write_location, read_location, interactive_location, query_location, &
+ LocationDims, LocationName, LocationLName, get_close_obs, &
+ get_close_maxdist_init, get_close_obs_init, get_close_type, &
+ operator(==), operator(/=), get_dist, get_close_obs_destroy, &
+ nc_write_location_atts, nc_get_location_varids, nc_write_location, &
+ vert_is_height, vert_is_pressure, vert_is_undef, vert_is_level, &
+ vert_is_surface, has_vertical_localization, VERTISSURFACE, &
+ VERTISLEVEL, VERTISHEIGHT
! version controlled file description for error handling, do not edit
character(len=128), parameter :: &
@@ -46,167 +46,129 @@
revdate = "$Date$"
type location_type
- private
+ private
real(r8) :: azm, rad, vloc
integer :: which_vert
- ! which_vert determines if the location is by level or by height/pressure
+ ! which_vert determines if the location is by level or by height
+ ! -1 ==> obs is on surface
! 1 ===> obs is by level
- ! 2 ===> obs is by pressure (Not supported in the annulus)
! 3 ===> obs is by height
end type location_type
+integer, parameter :: VERTISSURFACE = -1 ! surface
+integer, parameter :: VERTISLEVEL = 1 ! by level
+integer, parameter :: VERTISHEIGHT = 3 ! by height
+
type get_close_type
private
- integer :: maxdist
+ integer :: num
+ real(r8) :: maxdist
end type get_close_type
type(random_seq_type) :: ran_seq
logical :: ran_seq_init = .false.
-logical,save :: module_initialized = .false.
+logical, save :: module_initialized = .false.
-character(len=129) :: errstring
-
integer, parameter :: LocationDims = 3
character(len = 129), parameter :: LocationName = "loc_annulus"
character(len = 129), parameter :: LocationLName = &
- "Annulus location: azimuthal angle, radius, and level"
+ "Annulus location: azimuthal angle, radius, and height"
+! really just a placeholder. there was a comment that this code
+! needs a namelist with a min & max limit on the radius, but
+! the code no longer has one.
+real(r8) :: min_radius = 0.0_r8
+real(r8) :: max_radius = 100000.0_r8
+
+! FIXME: if we set mins & maxes, then we have to enforce the values
+! in various places in the code. none of that is there at this point.
+namelist /location_nml/ min_radius, max_radius
+
+character(len = 129) :: errstring
+
interface operator(==); module procedure loc_eq; end interface
interface operator(/=); module procedure loc_ne; end interface
+interface set_location
+ module procedure set_location_single
+ module procedure set_location_array
+end interface set_location
+
contains
+!----------------------------------------------------------------------------
+subroutine initialize_module
+
+! read namelist, set up anything that needs one-time initialization
- subroutine initialize_module
-!----------------------------------------------------------------------------
-! subroutine initialize_module
+integer :: iunit, io
- call register_module(source, revision, revdate)
- module_initialized = .true.
+! only do this code once
+if (module_initialized) return
-end subroutine initialize_module
+call register_module(source, revision, revdate)
+module_initialized = .true.
+! Read the namelist entry
+call find_namelist_in_file("input.nml", "location_nml", iunit)
+read(iunit, nml = location_nml, iostat = io)
+call check_namelist_read(iunit, io, "location_nml")
+! Write the namelist values to the log file
-function get_dist(loc1, loc2)
+if(do_nml_file()) write(nmlfileunit, nml=location_nml)
+if(do_nml_term()) write( * , nml=location_nml)
+
+! copy code from threed sphere module for handing the
+! options in the vertical? e.g. distances?
+
+end subroutine initialize_module
+
!----------------------------------------------------------------------------
-! To comply with current (Guam) configuration, distance depends only on
-! horizontal distance.
-!
+function get_dist(loc1, loc2, kind1, kind2)
+! Compute distance between 2 locations. Right now the distance only
+! depends on the horizontal. A namelist option might need to be added
+! that supports computing a true 3d distance.
+
type(location_type), intent(in) :: loc1, loc2
-real(r8) :: get_dist
+integer, optional, intent(in) :: kind1, kind2
+real(r8) :: get_dist
real(r8) :: x1, y1, x2, y2
if ( .not. module_initialized ) call initialize_module
+! FIXME: this does not take into account any vertical separation
! convert from cylindrical to cartesian coordinates
-x1 = loc1%rad*cos(loc1%azm)
-y1 = loc1%rad*sin(loc1%azm)
-x2 = loc2%rad*cos(loc2%azm)
-y2 = loc2%rad*sin(loc2%azm)
+x1 = loc1%rad * cos(loc1%azm)
+y1 = loc1%rad * sin(loc1%azm)
+x2 = loc2%rad * cos(loc2%azm)
+y2 = loc2%rad * sin(loc2%azm)
! Returns distance in m
get_dist = sqrt((x1 - x2)**2 + (y1 - y2)**2)
end function get_dist
-
-
-function get_location(loc)
!---------------------------------------------------------------------------
-!
-! Given a location type (where the azimuthal angle is in radians), this
-! routine return the azimuthal angle in degrees, the radius, and the level
-type(location_type), intent(in) :: loc
-real(r8), dimension(3) :: get_location
-
- if ( .not. module_initialized ) call initialize_module
-
- get_location(1) = loc%azm * RAD2DEG
- get_location(2) = loc%rad
- get_location(3) = loc%vloc
-
-end function get_location
-
-
-
-function vert_is_pressure(loc)
-!---------------------------------------------------------------------------
-!
-! Given a location, return true if vertical coordinate is pressure, else false.
-!
-! Always returns false, as vertical coordinate is never pressure for the annulus.
-
-logical :: vert_is_pressure
-type(location_type), intent(in) :: loc
-
-if ( .not. module_initialized ) call initialize_module
-
-vert_is_pressure = .false.
-
-end function vert_is_pressure
-
-
-
-function vert_is_height(loc)
-!---------------------------------------------------------------------------
-!
-! Given a location, return true if vertical coordinate is height, else false.
-
-logical :: vert_is_height
-type(location_type), intent(in) :: loc
-
-if ( .not. module_initialized ) call initialize_module
-
-if(loc%which_vert == 3 ) then
- vert_is_height = .true.
-else
- vert_is_height = .false.
-endif
-
-end function vert_is_height
-
-
-
-function vert_is_level(loc)
-!---------------------------------------------------------------------------
-!
-! Given a location, return true if vertical coordinate is pressure, else false.
-
-logical :: vert_is_level
-type(location_type), intent(in) :: loc
-
-if ( .not. module_initialized ) call initialize_module
-
-if(loc%which_vert == 1) then
- vert_is_level = .true.
-else
- vert_is_level = .false.
-endif
-
-end function vert_is_level
-
-
-
function loc_eq(loc1,loc2)
-!---------------------------------------------------------------------------
-!
-! interface operator used to compare two locations.
+
+! Interface operator used to compare two locations.
! Returns true only if all components are 'the same' to within machine
! precision.
type(location_type), intent(in) :: loc1, loc2
-logical :: loc_eq
+logical :: loc_eq
if ( .not. module_initialized ) call initialize_module
loc_eq = .false.
+if ( loc1%which_vert /= loc2%which_vert ) return
if ( abs(loc1%azm - loc2%azm ) > epsilon(loc1%azm ) ) return
if ( abs(loc1%rad - loc2%rad ) > epsilon(loc1%rad ) ) return
if ( abs(loc1%vloc - loc2%vloc) > epsilon(loc1%vloc) ) return
@@ -215,16 +177,15 @@
end function loc_eq
+!---------------------------------------------------------------------------
-
function loc_ne(loc1,loc2)
-!---------------------------------------------------------------------------
-!
-! interface operator used to compare two locations.
+
+! Interface operator used to compare two locations.
! Returns true if locations are not identical to machine precision.
type(location_type), intent(in) :: loc1, loc2
-logical :: loc_ne
+logical :: loc_ne
if ( .not. module_initialized ) call initialize_module
@@ -232,235 +193,255 @@
end function loc_ne
-
-
-function get_location_lon(loc)
!---------------------------------------------------------------------------
-!
-! Given a longitude location, return the azimuthal angle in degrees
-type(location_type), intent(in) :: loc
-real(r8) :: get_location_lon
+function get_location(loc)
+
+! Given a location type (where the azimuthal angle is in radians), this
+! routine return the azimuthal angle in degrees, the radius, and the vert
-if ( .not. module_initialized ) call initialize_module
-
-get_location_lon = loc%azm * RAD2DEG
-
-end function get_location_lon
-
-
-
-function get_location_lat(loc)
-!---------------------------------------------------------------------------
-!
-! Given a latitude location, return the radius
-
type(location_type), intent(in) :: loc
-real(r8) :: get_location_lat
+real(r8), dimension(3) :: get_location
if ( .not. module_initialized ) call initialize_module
-get_location_lat = loc%rad
+get_location(1) = loc%azm * RAD2DEG
+get_location(2) = loc%rad
+get_location(3) = loc%vloc
-end function get_location_lat
+end function get_location
+!----------------------------------------------------------------------------
-
-function set_location(lon, lat, vert_loc, which_vert)
-!----------------------------------------------------------------------------
-!
+function set_location_single(azm, rad, vert_loc, which_vert)
+
! Puts the given azimuthal angle (in degrees), radius, and vertical
! location into a location datatype.
-type (location_type) :: set_location
-real(r8), intent(in) :: lon, lat
+real(r8), intent(in) :: azm, rad
real(r8), intent(in) :: vert_loc
integer, intent(in) :: which_vert
+type (location_type) :: set_location_single
if ( .not. module_initialized ) call initialize_module
-if(lon < 0.0_r8 .or. lon > 360.0_r8) then
- write(errstring,*)'azimuthal angle (',lon,') is not within range [0,360]'
+if(azm < 0.0_r8 .or. azm > 360.0_r8) then
+ write(errstring,*)'azimuthal angle (',azm,') is not within range [0,360]'
call error_handler(E_ERR, 'set_location', errstring, source, revision, revdate)
endif
-set_location%azm = lon * DEG2RAD
-set_location%rad = lat
+set_location_single%azm = azm * DEG2RAD
+set_location_single%rad = rad
-if(which_vert /= -1 .and. which_vert /= 1 .and. which_vert /= 3 ) then
+if(which_vert /= VERTISSURFACE .and. &
+ which_vert /= VERTISLEVEL .and. &
+ which_vert /= VERTISHEIGHT ) then
write(errstring,*)'which_vert (',which_vert,') must be -1, 1 or 3'
call error_handler(E_ERR,'set_location', errstring, source, revision, revdate)
endif
-set_location%which_vert = which_vert
-set_location%vloc = vert_loc
+set_location_single%which_vert = which_vert
+set_location_single%vloc = vert_loc
-end function set_location
+end function set_location_single
+!----------------------------------------------------------------------------
-function set_location2(list)
-!----------------------------------------------------------------------------
-!
-! location semi-independent interface routine
+function set_location_array(list)
+
+! Location semi-independent interface routine
! given 4 float numbers, call the underlying set_location routine
-type (location_type) :: set_location2
real(r8), intent(in) :: list(:)
+type (location_type) :: set_location_array
if ( .not. module_initialized ) call initialize_module
-if (size(list) /= 4) then
+if (size(list) < 4) then
write(errstring,*)'requires 4 input values'
- call error_handler(E_ERR, 'set_location2', errstring, source, revision, revdate)
+ call error_handler(E_ERR, 'set_location', errstring, source, revision, revdate)
endif
-set_location2 = set_location(list(1), list(2), list(3), nint(list(4)))
+set_location_array = set_location_single(list(1), list(2), list(3), nint(list(4)))
-end function set_location2
+end function set_location_array
+!----------------------------------------------------------------------------
-
function set_location_missing()
-!----------------------------------------------------------------------------
-!
+! Initialize a location type to indicate the contents are unset.
+
type (location_type) :: set_location_missing
-set_location_missing%azm = missing_r8
-set_location_missing%rad = missing_r8
-set_location_missing%vloc = missing_r8
-set_location_missing%which_vert = missing_i
+if ( .not. module_initialized ) call initialize_module
+set_location_missing%azm = MISSING_R8
+set_location_missing%rad = MISSING_R8
+set_location_missing%vloc = MISSING_R8
+set_location_missing%which_vert = MISSING_I
+
end function set_location_missing
+!---------------------------------------------------------------------------
-
-function query_location(loc,attr) result(fval)
-!---------------------------------------------------------------------------
-!
+function query_location(loc,attr)
+
! Returns the value of the attribute
!
type(location_type), intent(in) :: loc
character(len=*), optional, intent(in) :: attr
-real(r8) :: fval
+real(r8) :: query_location
-character(len=16) :: attribute
-
if ( .not. module_initialized ) call initialize_module
-attribute = 'which_vert'
+! see the long comment in this routine in the threed_sphere
+! module for warnings about compiler bugs before you change
+! this code.
-if (present(attr)) attribute = attr
-selectcase(adjustl(attribute))
- case ('which_vert','WHICH_VERT')
- fval = loc%which_vert
- case ('lat','LAT')
- fval = loc%rad
- case ('lon','LON')
- fval = loc%azm
- case ('vloc','VLOC')
- fval = loc%vloc
- case default
- fval = loc%which_vert
-end select
+query_location = loc%which_vert
+if (.not. present(attr)) return
+select case(attr)
+ case ('which_vert','WHICH_VERT')
+ query_location = loc%which_vert
+ case ('rad','RAD','radius','RADIUS')
+ query_location = loc%rad
+ case ('azm','AZM','azimuth','AZIMUTH')
+ query_location = loc%azm
+ case ('vloc','VLOC')
+ query_location = loc%vloc
+ case default
+ call error_handler(E_ERR, 'query_location:', &
+ 'Only "azm","rad","vloc","which_vert" are legal attributes to request from location', &
+ source, revision, revdate)
+end select
end function query_location
-subroutine write_location(ifile, loc, fform, charstring)
!----------------------------------------------------------------------------
-!
-! Writes location to the file. Implemented as a subroutine but could
-! rewrite as a function with error control info returned. For initial
-! implementation, file is just an integer file unit number.
-integer, intent(in) :: ifile
-type(location_type), intent(in) :: loc
-character(len=*), intent(in), optional :: fform
-character(len=*), intent(out), optional :: charstring
+subroutine write_location(locfile, loc, fform, charstring)
+
+! Writes a location to the file.
+! additional functionality: if optional argument charstring is specified,
+! it must be long enough to hold the string, and the location information is
+! written into it instead of to a file. fform must be ascii (which is the
+! default if not specified) to use this option.
-logical :: writebuf
-integer :: charlength
+integer, intent(in) :: locfile
+type(location_type), intent(in) :: loc
+character(len = *), intent(in), optional :: fform
+character(len = *), intent(out), optional :: charstring
+integer :: charlength
+logical :: writebuf
+character(len = 128) :: string1
+
+! 10 format(1x,3(f22.14,1x),i4) ! old
+10 format(1X,F21.16,2(1X,G25.16),1X,I2)
+
if ( .not. module_initialized ) call initialize_module
-! incoming char string must be long enough for format
+! writing to a file (normal use) or to a character buffer?
writebuf = present(charstring)
-if (writebuf) then
- charlength = 82
- if (len(charstring) < charlength) then
- write(errstring, *) 'charstring buffer must be at least ', charlength, ' chars long'
- call error_handler(E_ERR, 'write_location', errstring, source, revision, revdate)
+
+! output file; test for ascii or binary, write what's asked, and return
+if (.not. writebuf) then
+ if (ascii_file_format(fform)) then
+ ! Write out pressure or level along with integer tag
+ ! we know azm is between 0, 360, and which_vert is a single digit.
+ write(locfile, '(''loc3a'')' )
+ write(locfile, 10) loc%azm, loc%rad, loc%vloc, loc%which_vert
+ else
+ write(locfile) loc%azm, loc%rad, loc%vloc, loc%which_vert
endif
+ return
endif
+! you only get here if you're writing to a buffer and not
+! to a file, and you can't have binary format set.
+if (.not. ascii_file_format(fform)) then
+ call error_handler(E_ERR, 'write_location', &
+ 'Cannot use string buffer with binary format', &
+ source, revision, revdate)
+endif
-if (ascii_file_format(fform)) then
- if (writebuf) then
- write(charstring,'(3(f22.14,1x),i4)') &
- loc%azm, loc%rad, loc%vloc, loc%which_vert
- else
- ! Write out pressure or level along with integer tag
- write(ifile, '(''loc3d'')' )
- write(ifile, '(1x,3(f22.14,1x),i4)')loc%azm, loc%rad, loc%vloc, loc%which_vert
- endif
-else
- if (writebuf) then
- call error_handler(E_ERR, 'write_location', &
- 'Cannot use string buffer with binary format', &
- source, revision, revdate)
- endif
+! format the location to be more human-friendly; meaning
+! degrees instead of radians, and kilometers for height,
+! hectopascals instead of pascals for pressure, etc.
- write(ifile)loc%azm, loc%rad, loc%vloc, loc%which_vert
+! this must be the sum of the longest of the formats below.
+charlength = 85
+
+if (len(charstring) < charlength) then
+ write(errstring, *) 'charstring buffer must be at least ', charlength, ' chars long'
+ call error_handler(E_ERR, 'write_location', errstring, source, revision, revdate)
endif
-end subroutine write_location
+write(string1, '(A,F12.8,1X,G15.8,A)') 'Azm(deg)/Radius(m): ', &
+ loc%azm*RAD2DEG, loc%rad, ' Depth:'
+! i am attempting to make these line up so if you have a list of mixed
+! vertical units, they all take the same number of columns. thus the extra
+! white space around some of the labels below.
+select case (loc%which_vert)
+ case (VERTISSURFACE)
+ write(charstring, '(A,1X,G15.6,A)') trim(string1), loc%vloc, ' surface (hPa)'
+ case (VERTISLEVEL)
+ write(charstring, '(A,1X,F6.0,A)') trim(string1), loc%vloc, ' level'
+ case (VERTISHEIGHT)
+ write(charstring, '(A,1X,G15.6,A)') trim(string1), loc%vloc / 1000.0_r8, ' km'
+ case default
+ write(errstring, *) 'unrecognized key for vertical type: ', loc%which_vert
+ call error_handler(E_ERR, 'write_location', errstring, source, revision, revdate)
+end select
+end subroutine write_location
-function read_location(ifile, fform)
!----------------------------------------------------------------------------
-!
+
+function read_location(locfile, fform)
+
! Reads location from file that was written by write_location.
! See write_location for additional discussion.
-integer, intent(in) :: ifile
-type(location_type) :: read_location
-character(len=*), intent(in), optional :: fform
+integer, intent(in) :: locfile
+type(location_type) :: read_location
+character(len = *), intent(in), optional :: fform
-character(len=5) :: header
+character(len=5) :: header
if ( .not. module_initialized ) call initialize_module
if (ascii_file_format(fform)) then
- read(ifile, '(a5)' ) header
+ read(locfile, '(a5)' ) header
- if(header /= 'loc3d') then
- write(errstring,*)'Expected location header "loc3d" in input file, got ', header
+ if(header /= 'loc3a') then
+ write(errstring,*)'Expected location header "loc3a" in input file, got ', header
call error_handler(E_ERR, 'read_location', errstring, source, revision, revdate)
endif
! Now read the location data value
- read(ifile, *)read_location%azm, read_location%rad, &
+ read(locfile, *) read_location%azm, read_location%rad, &
+ read_location%vloc, read_location%which_vert
+else
+ read(locfile) read_location%azm, read_location%rad, &
read_location%vloc, read_location%which_vert
-else
- read(ifile)read_location%azm, read_location%rad, &
- read_location%vloc, read_location%which_vert
endif
end function read_location
+!--------------------------------------------------------------------------
-
subroutine interactive_location(location, set_to_default)
-!--------------------------------------------------------------------------
-!
+
! Allows for interactive input of a location. Also gives option of selecting
-! a uniformly distributed random location (what the heck).
+! a uniformly distributed random location.
type(location_type), intent(out) :: location
-logical, intent(in), optional :: set_to_default
+logical, intent(in), optional :: set_to_default
real(r8) :: azm, rad, minazm, maxazm, minrad, maxrad
@@ -472,7 +453,7 @@
location%azm = 0.0_r8
location%rad = 0.0_r8
location%vloc = 0.0_r8
- location%which_vert = 0
+ location%which_vert = 0 ! zero is an invalid vert type
return
endif
endif
@@ -481,17 +462,17 @@
write(*, *)'-1 -> surface, 1 -> model level, 3 -> depth'
100 read(*, *) location%which_vert
-if(location%which_vert == 1 ) then
+if(location%which_vert == VERTISLEVEL ) then
write(*, *) 'Vertical co-ordinate model level'
read(*, *) location%vloc
-else if(location%which_vert == 3 ) then
+else if(location%which_vert == VERTISHEIGHT ) then
write(*, *) 'Vertical co-ordinate depth (in negative m)'
read(*, *) location%vloc
do while (location%vloc > 0)
write(*, *) 'Depth must be negative (zero at top of fluid), please try again'
read(*, *) location%vloc
end do
-else if(location%which_vert == -1 ) then
+else if(location%which_vert == VERTISSURFACE ) then
write(*, *) 'Vertical co-ordinate surface pressure (in hPa)'
read(*, *) location%vloc
location%vloc = 100.0 * location%vloc
@@ -518,38 +499,65 @@
ran_seq_init = .TRUE.
endif
- write(*, *) 'Input minimum azimuthal angle (0 to 360.0)'
- read(*, *) minazm
+ minazm = -1.0
+ do while (minazm < 0.0 .or. minazm > 360.0)
+ write(*, *) 'Input minimum azimuthal angle (0 to 360.0)'
+ read(*, *) minazm
+ if (minazm < 0.0 .or. minazm > 360.0) then
+ write(*, *) 'Angle must be between 0 to 360.0'
+ endif
+ enddo
minazm = minazm * DEG2RAD
- write(*, *) 'Input maximum azimuthal angle(0 to 360.0)'
- read(*, *) maxazm
+ maxazm = -1.0
+ do while (maxazm < 0.0 .or. maxazm > 360.0)
+ write(*, *) 'Input maximum azimuthal angle (0 to 360.0)'
+ read(*, *) maxazm
+ if (maxazm < 0.0 .or. maxazm > 360.0) then
+ write(*, *) 'Angle must be between 0 to 360.0'
+ endif
+ enddo
maxazm = maxazm * DEG2RAD
- ! Longitude is random from minlon to maxlon
+ ! Azimuth is random from minazm to maxazm, handle wrap around 360.0
+ if (minazm > maxazm) maxazm = maxazm + 2.0_r8 * PI
location%azm = random_uniform(ran_seq) * (maxazm-minazm) + minazm
+ if (location%azm > 2.0_r8 * PI) location%azm = location%azm - 2.0_r8 * PI
- write(*, *) 'Input minimum radius '
- read(*, *) minrad
+ minrad = -1.0
+ do while (minrad < min_radius)
+ write(*, *) 'Input minimum radius '
+ read(*, *) minrad
+ if (minrad < min_radius) then
+ write(*, *) 'Radius must be larger or equal to ', min_radius
+ endif
+ enddo
- write(*, *) 'Input maximum radius '
- read(*, *) maxrad
+ maxrad = -1.0
+ do while (maxrad > max_radius .or. maxrad <= minrad)
+ write(*, *) 'Input maximum radius '
+ read(*, *) maxrad
+ if (maxrad > max_radius .or. maxrad <= minrad) then
+ write(*, *) 'Radius must be greater than minrad and less than ', max_radius
+ endif
+ enddo
- ! Latitude must be area weighted to obtain proper random realizations
- location%rad = sqrt(random_uniform(ran_seq) * (maxrad-minrad) + minrad)
+ ! Radius must be area weighted to obtain proper random realizations
+ location%rad = sqrt(random_uniform(ran_seq)) * (maxrad-minrad) + minrad
write(*, *) 'random location is ', location%azm / DEG2RAD, &
location%rad
else
- write(*, *) 'Input radius '
- read(*, *) rad
-
- do while(rad < 0)
- write(*, *) 'Input value is illegal, please try again'
+ rad = -1.0
+ do while (rad < min_radius .or. rad > max_radius)
+ write(*, *) 'Input radius '
read(*, *) rad
- end do
+ if (rad < min_radius .or. rad > max_radius) then
+ write(*, *) 'Radius must be between ', min_radius, ' and ', max_radius
+ endif
+ enddo
location%rad = rad
location%azm = azm*DEG2RAD
@@ -558,95 +566,166 @@
end subroutine interactive_location
+!----------------------------------------------------------------------------
+function nc_write_location_atts( ncFileID, fname, ObsNumDimID ) result (ierr)
+
+! Writes the "location module" -specific attributes to a netCDF file.
-subroutine nc_write_location(ncFileID, LocationVarID, loc, start)
-!----------------------------------------------------------------------------
-!
-! Writes a SINGLE location to the specified netCDF variable and file.
-!
-
use typeSizes
use netcdf
-integer, intent(in) :: ncFileID, LocationVarID
-type(location_type), intent(in) :: loc
-integer, intent(in) :: start
+integer, intent(in) :: ncFileID ! handle to the netcdf file
+character(len=*), intent(in) :: fname ! file name (for printing purposes)
+integer, intent(in) :: ObsNumDimID ! handle to the dimension that grows
+integer :: ierr
+integer :: LocDimID
+integer :: VarID
+
if ( .not. module_initialized ) call initialize_module
-call check(nf90_put_var(ncFileID, LocationVarID, loc%azm, (/ start, 1 /) ))
-call check(nf90_put_var(ncFileID, LocationVarID, loc%rad, (/ start, 2 /) ))
-call check(nf90_put_var(ncFileID, LocationVarID, loc%vloc, (/ start, 3 /) ))
+ierr = -1 ! assume things will fail ...
-contains
-
- ! Internal subroutine - checks error status after each netcdf, prints
- ! text message each time an error code is returned.
- subroutine check(status)
- integer, intent ( in) :: status
+! define the rank/dimension of the location information
+call nc_check(nf90_def_dim(ncid=ncFileID, name='location', len=LocationDims, &
+ dimid = LocDimID), 'nc_write_location_atts', 'def_dim:location '//trim(fname))
- if(status /= nf90_noerr) call error_handler(E_ERR,'nc_write_location', &
- trim(nf90_strerror(status)), source, revision, revdate )
+! Define the observation location variable and attributes
- end subroutine check
+call nc_check(nf90_def_var(ncid=ncFileID, name='location', xtype=nf90_double, &
+ dimids=(/ LocDimID, ObsNumDimID /), varid=VarID), &
+ 'nc_write_location_atts', 'location:def_var')
-end subroutine nc_write_location
+call nc_check(nf90_put_att(ncFileID, VarID, 'description', &
+ 'location coordinates'), 'nc_write_location_atts', 'location:description')
+call nc_check(nf90_put_att(ncFileID, VarID, 'location_type', &
+ trim(LocationName)), 'nc_write_location_atts', 'location:location_type')
+call nc_check(nf90_put_att(ncFileID, VarID, 'long_name', &
+ trim(LocationLName)), 'nc_write_location_atts', 'location:long_name')
+call nc_check(nf90_put_att(ncFileID, VarID, 'storage_order', &
+ 'Azimuth Radius Vertical'), 'nc_write_location_atts', 'location:storage_order')
+call nc_check(nf90_put_att(ncFileID, VarID, 'units', &
+ 'degrees meters which_vert'), 'nc_write_location_atts', 'location:units')
+! Define the ancillary vertical array and attributes
+call nc_check(nf90_def_var(ncid=ncFileID, name='which_vert', xtype=nf90_int, &
+ dimids=(/ ObsNumDimID /), varid=VarID), &
+ 'nc_write_location_atts', 'which_vert:def_var')
+
+call nc_check(nf90_put_att(ncFileID, VarID, 'long_name', 'vertical coordinate system code'), &
+ 'nc_write_location_atts', 'which_vert:long_name')
+call nc_check(nf90_put_att(ncFileID, VarID, 'VERTISSURFACE', VERTISSURFACE), &
+ 'nc_write_location_atts', 'which_vert:VERTISSURFACE')
+call nc_check(nf90_put_att(ncFileID, VarID, 'VERTISLEVEL', VERTISLEVEL), &
+ 'nc_write_location_atts', 'which_vert:VERTISLEVEL')
+call nc_check(nf90_put_att(ncFileID, VarID, 'VERTISHEIGHT', VERTISHEIGHT), &
+ 'nc_write_location_atts', 'which_vert:VERTISHEIGHT')
+
+! If we made it to here without error-ing out ... we're good.
+
+ierr = 0
+
+end function nc_write_location_atts
+
!----------------------------------------------------------------------------
-subroutine get_close_obs_destroy(gc)
+subroutine nc_get_location_varids( ncFileID, fname, LocationVarID, WhichVertVarID )
-type(get_close_type), intent(inout) :: gc
+! Return the LocationVarID and WhichVertVarID variables from a given netCDF file.
+!
+! ncFileId the netcdf file descriptor
+! fname the name of the netcdf file (for error messages only)
+! LocationVarID the integer ID of the 'location' variable in the netCDF file
+! WhichVertVarID the integer ID of the 'which_vert' variable in the netCDF file
+use typeSizes
+use netcdf
-end subroutine get_close_obs_destroy
+integer, intent(in) :: ncFileID ! handle to the netcdf file
+character(len=*), intent(in) :: fname ! file name (for printing purposes)
+integer, intent(out) :: LocationVarID, WhichVertVarID
+if ( .not. module_initialized ) call initialize_module
+
+call nc_check(nf90_inq_varid(ncFileID, 'location', varid=LocationVarID), &
+ 'nc_get_location_varids', 'inq_varid:location '//trim(fname))
+
+call nc_check(nf90_inq_varid(ncFileID, 'which_vert', varid=WhichVertVarID), &
+ 'nc_get_location_varids', 'inq_varid:which_vert '//trim(fname))
+
+end subroutine nc_get_location_varids
+
!----------------------------------------------------------------------------
-subroutine get_close_maxdist_init(gc, maxdist)
+subroutine nc_write_location(ncFileID, LocationVarID, loc, obsindex, WhichVertVarID)
+
+! Writes a SINGLE location to the specified netCDF variable and file.
+! The LocationVarID and WhichVertVarID must be the values returned from
+! the nc_get_location_varids call.
-type(get_close_type), intent(inout) :: gc
-real(r8), intent(in) :: maxdist
+use typeSizes
+use netcdf
-gc%maxdist = maxdist
+integer, intent(in) :: ncFileID, LocationVarID
+type(location_type), intent(in) :: loc
+integer, intent(in) :: obsindex
+integer, intent(in) :: WhichVertVarID
-end subroutine get_close_maxdist_init
+real(r8), dimension(LocationDims) :: locations
+integer, dimension(1) :: intval
+if ( .not. module_initialized ) call initialize_module
+
+locations = get_location( loc ) ! converts from radians to degrees, btw
+
+call nc_check(nf90_put_var(ncFileID, LocationVarId, locations, &
+ start=(/ 1, obsindex /), count=(/ LocationDims, 1 /) ), &
+ 'nc_write_location', 'put_var:location')
+
+intval = loc%which_vert
+call nc_check(nf90_put_var(ncFileID, WhichVertVarID, intval, &
+ start=(/ obsindex /), count=(/ 1 /) ), &
+ 'nc_write_location','put_var:vert' )
+
+end subroutine nc_write_location
+
!----------------------------------------------------------------------------
subroutine get_close_obs_init(gc, num, obs)
-
+
! Initializes part of get_close accelerator that depends on the particular obs
type(get_close_type), intent(inout) :: gc
integer, intent(in) :: num
type(location_type), intent(in) :: obs(num)
+! Set the value of num_obs in the structure
+gc%num = num
+
end subroutine get_close_obs_init
-
!----------------------------------------------------------------------------
-subroutine alloc_get_close_obs(num, obs, cutoff, obs_box)
-! CURRENTLY UNUSED, I BELIEVE.
+subroutine get_close_obs_destroy(gc)
-integer, intent(in) :: num
-type(location_type), intent(in) :: obs(num)
-real(r8), intent(in) :: cutoff
-integer, intent(out) :: obs_box(num)
+type(get_close_type), intent(inout) :: gc
-! There might not need to be code here but if the get_close_obs() call
-! gets too slow, precomputing can be done here.
+end subroutine get_close_obs_destroy
-! set this to satisfy the intent(out) directive.
-obs_box(:) = 0
-return
+!----------------------------------------------------------------------------
-end subroutine alloc_get_close_obs
+subroutine get_close_maxdist_init(gc, maxdist)
+type(get_close_type), intent(inout) :: gc
+real(r8), intent(in) :: maxdist
+! Set the maximum distance in the structure
+gc%maxdist = maxdist
+
+end subroutine get_close_maxdist_init
+
!----------------------------------------------------------------------------
subroutine get_close_obs(gc, base_obs_loc, base_obs_kind, obs, obs_kind, &
@@ -656,20 +735,21 @@
! count, and the actual distances if requested. The kinds are available if
! more sophisticated distance computations are wanted.
-type(get_close_type), intent(in) :: gc
-type(location_type), intent(in) :: base_obs_loc, obs(:)
-integer, intent(in) :: base_obs_kind, obs_kind(:)
-integer, intent(out) :: num_close, close_ind(:)
-real(r8), optional, intent(out) :: dist(:)
+type(get_close_type), intent(in) :: gc
+type(location_type), intent(in) :: base_obs_loc, obs(:)
+integer, intent(in) :: base_obs_kind, obs_kind(:)
+integer, intent(out) :: num_close, close_ind(:)
+real(r8), optional, intent(out) :: dist(:)
integer :: i
real(r8) :: this_dist
-! Return list of obs that are within cutoff and their distances
+! Return list of obs that are within maxdist and their distances
num_close = 0
-do i = 1, size(obs)
- this_dist = get_dist(base_obs_loc, obs(i))
- if(this_dist < gc%maxdist) then
+!do i = 1, size(obs) ! i believe this is right
+do i = 1, gc%num
+ this_dist = get_dist(base_obs_loc, obs(i), base_obs_kind, obs_kind(i))
+ if(this_dist <= gc%maxdist) then
! Add this ob to the list
num_close = num_close + 1
close_ind(num_close) = i
@@ -679,17 +759,15 @@
end subroutine get_close_obs
+!----------------------------------------------------------------------------
-
function is_location_in_region(loc, minl, maxl)
-!----------------------------------------------------------------------------
-!
+
! Returns true if the given location is between the other two.
logical :: is_location_in_region
type(location_type), intent(in) :: loc, minl, maxl
-
if ( .not. module_initialized ) call initialize_module
if ((minl%which_vert /= maxl%which_vert) .or. &
@@ -702,7 +780,8 @@
! set to success only at the bottom after all tests have passed.
is_location_in_region = .false.
-if ((loc%azm < minl%azm) .or. (loc%azm > maxl%azm)) return
+! use the code in the utils module that knows how to wrap longitude/radians.
+if (.not. is_longitude_between(loc%azm, minl%azm, maxl%azm, doradians=.true.)) return
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list