[Dart-dev] [3567] DART/trunk/models: Interface between COAMPS and
DART.
nancy at ucar.edu
nancy at ucar.edu
Tue Aug 19 14:53:18 MDT 2008
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20080819/f0c2c7b3/attachment-0001.html
-------------- next part --------------
Added: DART/trunk/models/coamps/.model_mod.f90.swp
===================================================================
(Binary files differ)
Property changes on: DART/trunk/models/coamps/.model_mod.f90.swp
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: DART/trunk/models/coamps/byteswap.c
===================================================================
--- DART/trunk/models/coamps/byteswap.c (rev 0)
+++ DART/trunk/models/coamps/byteswap.c 2008-08-19 20:53:17 UTC (rev 3567)
@@ -0,0 +1,122 @@
+/*******************************
+ * MODULE: coamps_util_mod
+ * AUTHOR: T. R. Whitcomb
+ * Naval Research Laboratory
+ * DART VERSION: Jamaica
+ *
+ * Routines for in-place byteswapping of FORTRAN arrays
+ *******************************/
+
+/******************************
+ * BEGIN DEFINITIONS
+ ******************************/
+
+#define IN
+#define OUT
+#define INOUT
+
+#define INLINE
+
+/******************************
+ * END DEFINITIONS
+ ******************************/
+
+/******************************
+ * BEGIN PROTOTYPES
+ ******************************/
+
+/* Treat as public routine */
+void byteswap_array_( unsigned char*,
+ int*,
+ const int*);
+
+/* Treat as private routines */
+INLINE void byteswap_value( unsigned char*,
+ const int);
+
+INLINE void swap_bytes( unsigned char*,
+ unsigned char*);
+
+/******************************
+ * END PROTOTYPES
+ ******************************/
+
+/******************************
+ * BEGIN ROUTINES
+ ******************************/
+
+/**
+ * c_byteswap_array
+ * ----------------
+ * Reverses the byte ordering of an entire array elementwise in
+ * place. For example:
+ * F90:
+ * real(kind=8), dimension(5) :: test_array
+ * ...
+ * call byteswap_array(test_array, 5, 8)
+ *
+ * PARAMETERS
+ *INOUT array The array to byteswap
+ * IN array_len The length of the array (in units of
+ * element_size)
+ * IN element_size The size (in bytes) of each element of
+ * the array
+ */
+void c_byteswap_array_( INOUT unsigned char *array,
+ IN int *array_len,
+ IN const int *element_size)
+{
+ /* Need to keep track of where we are both in terms of the actual
+ * bytes and in terms of the original elements */
+ int ii;
+ int element_index;
+
+ for (ii=0; ii < *array_len; ii++)
+ {
+ element_index = ii * *element_size;
+ byteswap_value((array + element_index), *element_size);
+ }
+}
+
+/**
+ * byteswap_value
+ * --------------
+ * Reverses the byte ordering of a given value in-place
+ *
+ * PARAMETERS
+ *INOUT value pointer to value to swap bytes in
+ * IN value_size The size in bytes of this value
+ * (assumed even)
+ */
+INLINE void byteswap_value( INOUT unsigned char *value,
+ IN const int value_size)
+{
+ int ii;
+ for (ii = 0; ii < (value_size / 2); ii++)
+ {
+ swap_bytes(value + ii, value + (value_size - ii - 1));
+ }
+}
+
+/**
+ * swap_bytes
+ * ----------
+ * Switches the value of two parameters
+ *
+ * PARAMETERS
+ *INOUT value1 will hold value2 on exit
+ *INOUT value2 will hold value1 on exit
+ */
+INLINE void swap_bytes( INOUT unsigned char *value1,
+ INOUT unsigned char *value2)
+{
+ unsigned char temp_val;
+
+ temp_val = *value1;
+ *value1 = *value2;
+ *value2 = temp_val;
+}
+
+/******************************
+ * END ROUTINES
+ ******************************/
Added: DART/trunk/models/coamps/check_in_grid.f90
===================================================================
--- DART/trunk/models/coamps/check_in_grid.f90 (rev 0)
+++ DART/trunk/models/coamps/check_in_grid.f90 2008-08-19 20:53:17 UTC (rev 3567)
@@ -0,0 +1,42 @@
+! Program wrapper for ll2ij that takes its input from the file
+! domain.dat - this allows us to quickly cal ll2ij many times for a
+! single grid. This is written to print out either a "T" or "F" if
+! the lat/lon point described by the last two lines of domain.dat is
+! in the grid or not - this is for use with actual observations to
+! pare down the list to points that are actually within our domain.
+program check_in_grid
+
+ use coamps_intrinsic_mod, only : ll2ij
+
+ use types_mod, only : r8
+
+ implicit none
+
+ character(len=10), parameter :: FILENAME = 'domain.dat'
+ integer :: fileid
+ real(kind=r8), dimension(16) :: data
+
+ integer :: ii
+ real(kind=r8) :: grid_i, grid_j
+ logical :: in_grid
+
+ ! Pull the numbers out of the file - both domain information and
+ ! also the lat/lon point to convert
+ fileid = 13
+ open(unit=fileid,file=FILENAME);
+ read(fileid,*) data
+ close(fileid)
+
+ call ll2ij(int(data(1)), data(5), data(6), int(data(13)), &
+ int(data(14)), data(2), data(3), data(4), data(7), &
+ data(8), (/ grid_i /), (/ grid_j /), 1, (/data(15)/),&
+ (/ data(16) /))
+
+ ! Check using optimism
+ in_grid = .true.
+ if (grid_i > data(9)) in_grid = .false.
+ if (grid_j > data(10)) in_grid = .false.
+ if (grid_i < 1) in_grid = .false.
+ if (grid_j < 1) in_grid = .false.
+ print *, in_grid
+end program check_in_grid
Property changes on: DART/trunk/models/coamps/check_in_grid.f90
___________________________________________________________________
Name: svn:keywords
+ Date Revision Author HeadURL Id
Added: DART/trunk/models/coamps/coamps_grid_mod.f90
===================================================================
--- DART/trunk/models/coamps/coamps_grid_mod.f90 (rev 0)
+++ DART/trunk/models/coamps/coamps_grid_mod.f90 2008-08-19 20:53:17 UTC (rev 3567)
@@ -0,0 +1,1185 @@
+!------------------------------
+! MODULE: coamps_grid_mod
+! AUTHOR: T. R. Whitcomb
+! Naval Research Laboratory
+! DART VERSION: Jamaica
+!
+! Module containing the definition of a COAMPS grid data structure
+! and routines dealing with COAMPS grid information.
+!------------------------------
+module coamps_grid_mod
+
+ use coamps_intrinsic_mod, only : ij2ll, &
+ ll2ij
+ use coamps_util_mod, only : C_REAL, &
+ check_io_status, &
+ check_alloc_status, &
+ check_dealloc_status, &
+ fix_for_platform
+
+ use location_mod, only : get_location, &
+ location_type, &
+ set_location
+
+ use utilities_mod, only : do_output, &
+ E_ERR, &
+ E_WARN, &
+ error_handler, &
+ get_unit
+ use types_mod, only : r8
+
+ implicit none
+
+ private
+
+ !------------------------------
+ ! BEGIN PUBLIC INTERFACE
+ !------------------------------
+
+ ! COAMPS grid data structure
+ public :: coamps_grid
+ public :: initialize_grid
+
+ ! Location conversion functions
+ public :: gridpt_to_latlon
+ public :: latlon_to_gridpt
+ public :: gridpt_to_location
+ public :: location_to_gridpt
+ public :: grid_ij_to_vector_index
+
+ ! Information about the grid size
+ public :: get_grid_dims
+ public :: get_grid_field_size
+ public :: get_grid_num_levels
+ public :: check_ij_within_grid
+
+ ! Information about horizontal coordinates
+ public :: get_grid_delta_x
+ public :: get_grid_delta_y
+
+ ! Information about vertical coordinates
+ public :: get_grid_msigma
+ public :: get_grid_wsigma
+ public :: get_grid_dsigmaw
+ public :: get_terrain_height_at_points
+
+ ! Domain decomposition tools
+ public :: initialize_decomposition
+ public :: decompose_domain
+ public :: get_grid_iminf
+ public :: get_grid_imaxf
+ public :: get_grid_jminf
+ public :: get_grid_jmaxf
+
+ ! Debugging output
+ public :: dump_grid_info
+
+ !------------------------------
+ ! END PUBLIC INTERFACE
+ !------------------------------
+
+ !------------------------------
+ ! BEGIN EXTERNAL INTERFACES
+ !------------------------------
+ ! [none]
+ !------------------------------
+ ! END EXTERNAL INTERFACES
+ !------------------------------
+
+ !------------------------------
+ ! BEGIN TYPES AND CONSTANTS
+ !------------------------------
+
+ ! Derived data type for the COAMPS grid information - this is all
+ ! contained in the datahd_*_infofld files that the COAMPS analysis
+ ! produces, so we can just read that rather than processing a
+ ! namelist. Most of these are pretty straighforward, but:
+ ! anchor_i - the i-grid point where the grid is anchored to parent
+ ! anchor_j - the j-grid point where the grid is anchored to parent
+ ! iref - the i-grid point where the grid is anchored to Earth
+ ! jref - the j-grid point where the grid is anchored to Earth
+ ! dsigma - the spacing between mass sigma levels
+ ! msigma - actual sigma values at mass levels
+ ! terrain - terrain height
+ ! [ij]minf - the minimum grid point for process N
+ ! [ij]maxf - the maximum grid point for process N
+ ! NOTE: For now, "nests" is currently unused and the pts, anchor,
+ ! and [ij]ref variables are scalars - when multiple grid nests are
+ ! supported, these will be vectors. Currently only a single grid
+ ! is supported to make my life much, much, much easier.
+ type :: coamps_grid
+ real(kind=r8) :: pres_lvls
+ real(kind=r8) :: sigm_lvls
+ real(kind=r8) :: map_proj
+ real(kind=r8) :: stdlat1
+ real(kind=r8) :: stdlat2
+ real(kind=r8) :: stdlon
+ real(kind=r8) :: reflat
+ real(kind=r8) :: reflon
+ real(kind=r8) :: delta_x
+ real(kind=r8) :: delta_y
+ real(kind=r8) :: nests
+ real(kind=r8) :: pts_x
+ real(kind=r8) :: pts_y
+ real(kind=r8) :: anchor_i
+ real(kind=r8) :: anchor_j
+ real(kind=r8) :: iref
+ real(kind=r8) :: jref
+ real(kind=r8), dimension(:), allocatable :: dsigma
+ real(kind=r8), dimension(:), allocatable :: dsigmaw
+ real(kind=r8), dimension(:), allocatable :: msigma
+ real(kind=r8), dimension(:), allocatable :: wsigma
+ real(kind=r8), dimension(:,:), allocatable :: terrain
+
+ ! These fields are used for dealing with domain decomposition
+ integer, dimension(:), allocatable :: iminf
+ integer, dimension(:), allocatable :: imaxf
+ integer, dimension(:), allocatable :: jminf
+ integer, dimension(:), allocatable :: jmaxf
+ end type coamps_grid
+
+ integer, parameter :: SINGLE_POINT = 1
+
+ ! Model vertical coordinate indexes top-down
+ integer, parameter :: TOP_LEVEL = 1
+
+ ! Location of entries in the datahd file
+ integer, parameter :: DATAHD_NUM_P_LEVELS = 1
+ integer, parameter :: DATAHD_NUM_S_LEVELS = 2
+ integer, parameter :: DATAHD_MAP_PROJECTION = 3
+ integer, parameter :: DATAHD_STANDARD_LAT_1 = 4
+ integer, parameter :: DATAHD_STANDARD_LAT_2 = 5
+ integer, parameter :: DATAHD_STANDARD_LON = 6
+ integer, parameter :: DATAHD_REFERENCE_LAT = 7
+ integer, parameter :: DATAHD_REFERENCE_LON = 8
+ integer, parameter :: DATAHD_DELTA_X = 9
+ integer, parameter :: DATAHD_DELTA_Y = 10
+ integer, parameter :: DATAHD_NUM_NESTS = 11
+ integer, parameter :: DATAHD_NUM_X_POINTS = 30
+ integer, parameter :: DATAHD_NUM_Y_POINTS = 31
+ integer, parameter :: DATAHD_ANCHOR_POINT_I = 32
+ integer, parameter :: DATAHD_ANCHOR_POINT_J = 33
+ integer, parameter :: DATAHD_REFERENCE_I = 34
+ integer, parameter :: DATAHD_REFERENCE_J = 35
+ integer, parameter :: DATAHD_DSIGMA_OFFSET = 500
+ integer, parameter :: DATAHD_MSIGMA_OFFSET = 800
+
+ !------------------------------
+ ! END TYPES AND CONSTANTS
+ !------------------------------
+
+ !------------------------------
+ ! BEGIN MODULE VARIABLES
+ !------------------------------
+
+ ! Modified automatically by Subversion
+ character(len=128) :: &
+ source = "$URL$", &
+ revision = "$Revision$", &
+ revdate = "$Date$"
+
+ logical :: module_initialized = .false.
+
+ !------------------------------
+ ! END MODULE VARIABLES
+ !------------------------------
+
+contains
+
+ !------------------------------
+ ! BEGIN PUBLIC ROUTINES
+ !------------------------------
+
+ ! initialize_grid
+ ! ---------------
+ ! Populates a coamps_grid data structure using the COAMPS domain
+ ! data (datahd) file and the COAMPS terrain height (terrht) file
+ ! based on a given date-time group.
+ ! N.B. "grid" parameter is marked as intent(inout) to ensure that
+ ! grid will not be undefined on exit if we don't need to run
+ ! this routine
+ ! PARAMETERS
+ ! IN dtg base date-time group for model run
+ ! INOUT grid coamps_grid structure to initialize
+ subroutine initialize_grid(dtg, grid)
+ character(len=10), intent(in) :: dtg
+ type(coamps_grid), intent(inout) :: grid
+
+
+ character(len=64) :: datahd_filename
+ integer :: datahd_unit
+ character(len=64) :: terrht_filename
+ integer :: terrht_unit
+ integer :: terrht_record_len
+
+ ! Terrain data is stored in a flat file, so it uses the COAMPS
+ ! real number size. Single dimension makes byteswapping easier
+ real(kind=C_REAL), dimension(:), allocatable :: coamps_terrain
+
+ integer :: field_size
+
+ ! Error-checking info
+ character(len=*), parameter :: routine = 'initialize_grid'
+ integer :: io_status
+ integer :: alloc_status
+ integer :: dealloc_status
+
+ if (module_initialized) return
+
+ ! Initialize the domain data section of the coamps_grid structure
+ call generate_datahd_filename(dtg, datahd_filename)
+ datahd_unit = get_unit()
+ open(unit=datahd_unit, file=datahd_filename, status='old', &
+ access='sequential', action='read', form='formatted', &
+ iostat=io_status)
+ call check_io_status(io_status, routine, source, revision, &
+ revdate, 'Opening datahd file')
+ call read_datahd_file(datahd_unit, grid)
+ close(datahd_unit)
+
+ ! Need to set up permanent and temporary storage for terrain data
+ allocate( grid%terrain(int(grid%pts_x), int(grid%pts_y)), &
+ stat = alloc_status )
+ call check_alloc_status(alloc_status, routine, source, &
+ revision, revdate, 'grid%terrain')
+
+ call get_grid_field_size(grid, field_size)
+ allocate (coamps_terrain(field_size), stat = alloc_status)
+ call check_alloc_status(alloc_status, routine, source, revision, &
+ revdate, 'coamps_terrain (temporary)')
+
+ ! Make sure when checking the record length to use what *COAMPS*
+ ! wrote it as, not what we want to store it as
+ call generate_terrht_filename(dtg, grid, terrht_filename)
+ terrht_unit = get_unit()
+ inquire(iolength=terrht_record_len) coamps_terrain
+
+
+ open(unit = terrht_unit, file = terrht_filename, status = 'old',&
+ access = 'direct', action = 'read', form = 'unformatted', &
+ recl = terrht_record_len, iostat = io_status)
+ call check_io_status(io_status, routine, source, revision, &
+ revdate, 'Opening terrain height file')
+ call read_terrht_file(terrht_unit, coamps_terrain)
+ close(terrht_unit)
+
+ ! COAMPS writes big-endian flat files no matter what platform is
+ ! being used, so we may need to change the data we read.
+ call fix_for_platform(coamps_terrain, field_size, C_REAL)
+
+ grid%terrain = reshape(real(coamps_terrain, kind=r8), &
+ (/ int(grid%pts_x), int(grid%pts_y) /))
+
+ deallocate(coamps_terrain, stat=dealloc_status)
+ call check_dealloc_status(dealloc_status, routine, source, revision, &
+ revdate, 'coamps_terrain')
+
+ module_initialized = .true.
+ end subroutine initialize_grid
+
+ ! gridpt_to_latlon
+ ! ----------------
+ ! Given a COAMPS grid structure, convert an (i,j) representation
+ ! of a point on that grid to a (lat, lon) representation.
+ ! Note that this routine will happily accept and convert (i,j)
+ ! points that are beyond the bounds of the domain. For this
+ ! direction of conversion, that's not too big a problem.
+ ! PARAMETERS
+ ! IN grid coamps_grid to base conversion on
+ ! IN ii point's x-direction index
+ ! IN jj point's y-direction index
+ ! OUT lat point's latitude
+ ! OUT lon point's longitude
+ subroutine gridpt_to_latlon(grid, ii, jj, lat, lon)
+ type(coamps_grid), intent(in) :: grid
+ real(kind=r8), intent(in) :: ii
+ real(kind=r8), intent(in) :: jj
+ real(kind=r8), intent(out) :: lat
+ real(kind=r8), intent(out) :: lon
+
+ ! Work-around for not being able to pass in scalars as arrays
+ ! of length 1
+ real(kind=r8), dimension(SINGLE_POINT) :: ii_array, jj_array
+ real(kind=r8), dimension(SINGLE_POINT) :: lat_array, lon_array
+
+ ii_array(SINGLE_POINT) = ii
+ jj_array(SINGLE_POINT) = jj
+
+ call ij2ll(int(grid%map_proj), grid%reflat, grid%reflon, &
+ int(grid%iref), int(grid%jref), grid%stdlat1, &
+ grid%stdlat2, grid%stdlon, grid%delta_x, &
+ grid%delta_y, ii_array, jj_array, &
+ SINGLE_POINT, lat_array, lon_array)
+
+ lat = lat_array(SINGLE_POINT)
+ lon = lon_array(SINGLE_POINT)
+ end subroutine gridpt_to_latlon
+
+ ! latlon_to_gridpt
+ ! ----------------
+ ! Given a COAMPS grid structure, convert a (lat,lon)
+ ! representation of a point on that grid to an (i, j)
+ ! representation.
+ ! Note that this routine will happily accept and convert
+ ! (lat,lon) points that are beyond the bounds of the domain.
+ ! This could happen in several ways, but especially if an obs
+ ! gets included in the sequence that for some reason isn't in
+ ! the domain - in this case, warn the user.
+ ! PARAMETERS
+ ! IN grid coamps_grid to base conversion on
+ ! IN lat point's latitude
+ ! IN lon point's longitude
+ ! OUT ii point's x-direction index
+ ! OUT jj point's y-direction index
+ subroutine latlon_to_gridpt(grid, lat, lon, ii, jj)
+ type(coamps_grid), intent(in) :: grid
+ real(kind=r8), intent(in) :: lat
+ real(kind=r8), intent(in) :: lon
+ real(kind=r8), intent(out) :: ii
+ real(kind=r8), intent(out) :: jj
+
+ ! Work-around for not being able to pass in scalars as arrays
+ ! of length 1
+ real(kind=r8), dimension(SINGLE_POINT) :: ii_array, jj_array
+ real(kind=r8), dimension(SINGLE_POINT) :: lat_array, lon_array
+
+ ! If anything goes wrong in the conversion, pre-define values
+ ! that should set off alarm bells
+ ii = -999.0
+ jj = -999.0
+
+ lat_array(SINGLE_POINT) = lat
+ lon_array(SINGLE_POINT) = lon
+
+ call ll2ij(int(grid%map_proj), grid%reflat, grid%reflon, &
+ int(grid%iref), int(grid%jref), grid%stdlat1, &
+ grid%stdlat2, grid%stdlon, grid%delta_x, &
+ grid%delta_y, ii_array, jj_array, &
+ SINGLE_POINT, lat_array, lon_array)
+
+ ii = ii_array(SINGLE_POINT)
+ jj = jj_array(SINGLE_POINT)
+
+ if ( ii .lt. 1 .or. ii .gt. grid%pts_x .or. &
+ jj .lt. 1 .or. jj .gt. grid%pts_y) then
+ call error_handler(E_WARN, 'latlon_to_gridpt', source, &
+ revision, revdate, 'Point is' // &
+ 'outside domain')
+ end if
+ end subroutine latlon_to_gridpt
+
+ ! gridpt_to_location
+ ! ------------------
+ ! Given a COAMPS grid, an (i,j) horizontal location in that grid,
+ ! and a vertical location and vertical coordinate type, returns
+ ! a DART location structure for that point
+ ! PARAMETERS
+ ! IN grid coamps_grid for lat/lon conversion
+ ! IN ii point's x-direction index
+ ! IN jj point's y-direction index
+ ! IN vert_loc point's vertical direction location
+ ! IN vert_coord point's vertical direction coordinate
+ ! OUT loc point's location as a DART structure
+ subroutine gridpt_to_location(grid, ii, jj, vert_loc, vert_coord,&
+ loc)
+ type(coamps_grid), intent(in) :: grid
+ integer, intent(in) :: ii
+ integer, intent(in) :: jj
+ real(kind=r8), intent(in) :: vert_loc
+ integer, intent(in) :: vert_coord
+ type(location_type), intent(out) :: loc
+
+ real(kind=r8) :: real_ii, real_jj
+ real(kind=r8) :: point_lat, point_lon
+
+ ! Get things squared away in 2D first
+ real_ii = real(ii, kind=r8)
+ real_jj = real(jj, kind=r8)
+ call gridpt_to_latlon(grid, real_ii, real_jj, point_lat, &
+ point_lon)
+
+ ! Do this with keyword arguments because I can't keep the order
+ ! straight
+ loc = set_location( lon = point_lon, &
+ lat = point_lat, &
+ vert_loc = vert_loc, &
+ which_vert = vert_coord )
+ end subroutine gridpt_to_location
+
+ ! location_to_gridpt
+ ! ------------------
+ ! Given a COAMPS grid and a DART location, return the location
+ ! of a point as an (i, j) coordinate and a vertical coordinate.
+ ! Note that DART's get_location function that's used to break
+ ! this down does *not* return any information on what that
+ ! vertical coordinate is, so no information is returned to the
+ ! caller.
+ ! PARAMETERS
+ ! IN grid coamps_grid for lat/loon conversion
+ ! IN loc DART location structure to unpack
+ ! OUT ii point's x-direction index
+ ! OUT jj point's y-direction index
+ ! OUT vert_loc point's vertical location
+ subroutine location_to_gridpt(grid, loc, ii, jj, vert_loc)
+ type(coamps_grid), intent(in) :: grid
+ type(location_type), intent(in) :: loc
+ real(kind=r8), intent(out) :: ii
+ real(kind=r8), intent(out) :: jj
+ real(kind=r8), intent(out) :: vert_loc
+
+ ! The array returned by get_location is (lon, lat, vertical)
+ real(kind=r8), dimension(3) :: loc_array
+
+ integer, parameter :: DART_LOC_LON = 1
+ integer, parameter :: DART_LOC_LAT = 2
+ integer, parameter :: DART_LOC_VERT = 3
+
+ loc_array = get_location(loc)
+
+ ! Deal with the 2D component separately from the vertical
+ ! Again, be explicit since I apparently like to mix these up
+ call latlon_to_gridpt( grid = grid, &
+ lat = loc_array(DART_LOC_LAT), &
+ lon = loc_array(DART_LOC_LON), &
+ ii = ii, &
+ jj = jj )
+ vert_loc = loc_array(DART_LOC_VERT)
+ end subroutine location_to_gridpt
+
+ ! grid_ij_to_vector_index
+ ! -----------------------
+ ! Given a COAMPS grid, converts an horizontal (ii,jj) location
+ ! to an index in the vector representation
+ ! PARAMETERS
+ ! IN grid coamps_grid to use for the conversion
+ ! IN ii zonal index component
+ ! IN jj meridional index component
+ ! OUT index position that the point (ii,jj) has in
+ ! the field re-arranged in vector form
+ subroutine grid_ij_to_vector_index(grid, ii, jj, index)
+ type(coamps_grid), intent(in) :: grid
+ integer, intent(in) :: ii
+ integer, intent(in) :: jj
+ integer, intent(out) :: index
+
+ index = (jj - 1) * grid%pts_x + ii
+ end subroutine grid_ij_to_vector_index
+
+ ! get_grid_dims
+ ! -------------
+ ! Given a COAMPS grid structure, calculate the i (i.e. x) and
+ ! j (i.e. y) dimensions
+ ! PARAMETERS
+ ! IN grid coamps_grid to get limits from
+ ! OUT grid_i number of points on the i dimension
+ ! OUT grid_j number of points on the j dimension
+ subroutine get_grid_dims(grid, grid_i, grid_j)
+ type(coamps_grid), intent(in) :: grid
+ integer, intent(out) :: grid_i
+ integer, intent(out) :: grid_j
+
+ grid_i = int(grid%pts_x)
+ grid_j = int(grid%pts_y)
+ end subroutine get_grid_dims
+
+ ! get_grid_field_size
+ ! -------------------
+ ! Given a coamps grid structure, calculate the field size (the
+ ! total number of points in the i and j directions)
+ ! PARAMETERS
+ ! IN grid coamps_grid to get size from
+ ! OUT field_size number of grid points on a single level
+ subroutine get_grid_field_size(grid, field_size)
+ type(coamps_grid), intent(in) :: grid
+ integer, intent(out) :: field_size
+
+ field_size = int(grid%pts_x) * int(grid%pts_y)
+ end subroutine get_grid_field_size
+
+ ! get_grid_num_levels
+ ! -------------------
+ ! Returns the number of mass sigma levels on a given COAMPS grid
+ ! PARAMETERS
+ ! IN grid coamps_grid to get levels from
+ ! OUT num_sig_lvls number of sigma levels on the grid
+ subroutine get_grid_num_levels(grid, num_sig_lvls)
+ type(coamps_grid), intent(in) :: grid
+ integer, intent(out) :: num_sig_lvls
+
+ num_sig_lvls = grid%sigm_lvls
+ end subroutine get_grid_num_levels
+
+ ! check_ij_within_grid
+ ! --------------------
+ ! Calculates if a given (i,j) index is within a given COAMPS grid
+ ! PARAMETERS
+ ! IN grid coamps_grid to use for checking
+ ! IN ii zonal index
+ ! IN jj meridional index
+ ! OUT is_within_grid True if the (i,j) location listed is
+ ! inside the given grid
+ subroutine check_ij_within_grid(grid, ii, jj, is_within_grid)
+ type(coamps_grid), intent(in) :: grid
+ integer, intent(in) :: ii
+ integer, intent(in) :: jj
+ logical, intent(out) :: is_within_grid
+
+ is_within_grid = .true.
+ if ( (ii < 1) .or. (ii > grid%pts_x) .or. &
+ (jj < 1) .or. (jj > grid%pts_y) ) then
+ is_within_grid = .false.
+ end if
+ end subroutine check_ij_within_grid
+
+ ! get_grid_delta_x
+ ! ----------------
+ ! Returns the x-direction grid spacing (in meters) for the given grid
+ ! PARAMETERS
+ ! IN grid coamps_grid to pull data from
+ ! OUT grid_delta_x x-direction grid spacing
+ subroutine get_grid_delta_x(grid, grid_delta_x)
+ type(coamps_grid), intent(in) :: grid
+ real(kind=r8), intent(out) :: grid_delta_x
+
+ grid_delta_x = grid%delta_x
+ end subroutine get_grid_delta_x
+
+ ! get_grid_delta_y
+ ! ----------------
+ ! Returns the y-direction grid spacing (in meters) for the given grid
+ ! PARAMETERS
+ ! IN grid coamps_grid to pull data from
+ ! OUT grid_delta_y y-direction grid spacing
+ subroutine get_grid_delta_y(grid, grid_delta_y)
+ type(coamps_grid), intent(in) :: grid
+ real(kind=r8), intent(out) :: grid_delta_y
+
+ grid_delta_y = grid%delta_y
+ end subroutine get_grid_delta_y
+
+ ! get_grid_msigma
+ ! ---------------
+ ! Returns the mass sigma levels for a given COAMPS grid
+ ! PARAMETERS
+ ! IN grid coamps_grid to pull msigma out of
+ ! OUT msigma mass sigma levels for the grid
+ subroutine get_grid_msigma(grid, msigma)
+ type(coamps_grid), intent(in) :: grid
+ real(kind=r8), dimension(:), intent(out) :: msigma
+
+ integer :: kka
+
+ ! Don't attempt a copy if the data in the grid won't fit
+ if ( size(grid%msigma) > size(msigma) ) then
+ call error_handler(E_ERR, 'get_grid_msigma', 'msigma &
+ &passed into function is too small', &
+ &source, revision, revdate)
+ end if
+
+ do kka=1,size(grid%msigma)
+ msigma(kka) = grid%msigma(kka)
+ end do
+ end subroutine get_grid_msigma
+
+ ! get_grid_wsigma
+ ! ---------------
+ ! Returns the w sigma levels for a given COAMPS grid
+ ! PARAMETERS
+ ! IN grid coamps_grid to pull wsigma out of
+ ! OUT wsigma w sigma levels for the grid
+ subroutine get_grid_wsigma(grid, wsigma)
+ type(coamps_grid), intent(in) :: grid
+ real(kind=r8), dimension(:), intent(out) :: wsigma
+
+ integer :: kkw
+
+ ! Bail if the data in the grid won't fit
+ if ( size(grid%wsigma) > size(wsigma) ) then
+ call error_handler(E_ERR, 'get_grid_wsigma', 'wsigma ' // &
+ 'passed into function is too small', &
+ source, revision, revdate)
+ end if
+
+ do kkw=1,size(grid%wsigma)
+ wsigma(kkw) = grid%wsigma(kkw)
+ end do
+ end subroutine get_grid_wsigma
+
+ ! get_grid_dsigmaw
+ ! ----------------
+ ! Returns the distance between w sigma levels for a given COAMPS
+ ! grid
+ ! PARAMETERS
+ ! IN grid coamps_grid to pull dsigmaw out of
+ ! OUT dsigmaw distance between w sigma levels
+ subroutine get_grid_dsigmaw(grid, dsigmaw)
+ type(coamps_grid), intent(in) :: grid
+ real(kind=r8), dimension(:), intent(out) :: dsigmaw
+
+ integer :: kkw
+
+ ! Bail if the data in the grid won't fit
+ if ( size(grid%dsigmaw) > size(dsigmaw) ) then
+ call error_handler(E_ERR, 'get_grid_dsigmaw', 'dsigmaw ' // &
+ 'passed into function is too small', &
+ source, revision, revdate)
+ end if
+
+ do kkw=1,size(grid%dsigmaw)
+ dsigmaw(kkw) = grid%dsigmaw(kkw)
+ end do
+ end subroutine get_grid_dsigmaw
+
+ ! get_terrain_height_at_points
+ ! ----------------------------
+ ! Given a COAMPS grid and a vector representation of a set of
+ ! (i,j) points, returns the terrain height for those points
+ ! PARAMETERS
+ ! IN grid coamps_grid to pull terrht from
+ ! IN i_points points's x-direction index
+ ! IN j_points points's y-direction index
+ ! OUT terrain_height points's terrain height
+ subroutine get_terrain_height_at_points(grid, i_points, j_points,&
+ terrain_height)
+ type(coamps_grid), intent(in) :: grid
+ integer, dimension(:), intent(in) :: i_points
+ integer, dimension(:), intent(in) :: j_points
+ real(kind=r8), dimension(:), intent(out) :: terrain_height
+
+ integer :: nn, num_points
+
+ ! Sanity check
+ if (size(i_points) .ne. size(j_points)) then
+ call error_handler(E_ERR, 'get_terrain_height_at_points', &
+ 'Error in getting terrain height - ' // &
+ 'i/j point arrays are different sizes!', &
+ source, revision, revdate)
+ end if
+
+ num_points = size(i_points)
+ do nn=1,num_points
+ terrain_height(nn) = grid%terrain(i_points(nn), j_points(nn))
+ end do
+ end subroutine get_terrain_height_at_points
+
+ ! initialize_decomposition
+ ! ------------------------
+ ! Given a COAMPS grid the total number of non-I/O processors for
+ ! that grid, allocate space for each processor's [ij](min|max)f
+ ! extent
+ ! N.B. The grid is not defined as intent(out) here since that
+ ! could mean that the other components of the coamps_grid
+ ! are undefined after this routine call. This is avoided
+ ! by using intent(inout).
+ ! PARAMETERS
+ ! INOUT grid coamps_grid to allocate space in
+ ! IN totalprocs total number of non-I/O processors
+ subroutine initialize_decomposition(grid, totalprocs)
+ type(coamps_grid), intent(inout) :: grid
+ integer, intent(in) :: totalprocs
+
+ character(len=*), parameter :: routine = 'initialize_decomposition'
+ integer :: alloc_status
+
+ allocate( grid%iminf(totalprocs), stat=alloc_status )
+ call check_alloc_status(alloc_status, routine, source, revision, &
+ revdate, 'grid%iminf')
+ allocate( grid%imaxf(totalprocs), stat=alloc_status )
+ call check_alloc_status(alloc_status, routine, source, revision, &
+ revdate, 'grid%imaxf')
+ allocate( grid%jminf(totalprocs), stat=alloc_status )
+ call check_alloc_status(alloc_status, routine, source, revision, &
+ revdate, 'grid%jminf')
+ allocate( grid%jmaxf(totalprocs), stat=alloc_status )
+ call check_alloc_status(alloc_status, routine, source, revision, &
+ revdate, 'grid%jmaxf')
+ end subroutine initialize_decomposition
+
+ ! decompose_domain
+ ! ----------------
+ ! Given the grid, the number of I/O processors, the number of
+ ! processors in x and y, and the number of boundary points (nbnam
+ ! in the COAMPS namelist), mimics the domain decomposition
+ ! performed by COAMPS routine domdec.F (the math is taken from that
+ ! module and variables renamed for use here). The maximum and
+ ! minimum [ij] points for each processor are given in arrays
+ ! [ij](min|max)f within the coamps_grid structure. This routine
+ ! populates those arrays.
+ ! PARAMETERS
+ ! INOUT grid coamps_grid to decompose
+ ! IN num_io_procs number of I/O processors (0 or 1)
+ ! IN dom_proc_x number of processors in the x direction
+ ! IN dom_proc_y number of processors in the y direction
+ ! IN nbound number of halo boundary points
+ subroutine decompose_domain(grid, num_io_procs, dom_proc_x, &
+ dom_proc_y, nbound)
+ type(coamps_grid), intent(inout) :: grid
+ integer, intent(in) :: num_io_procs
+ integer, intent(in) :: dom_proc_x
+ integer, intent(in) :: dom_proc_y
+ integer, intent(in) :: nbound
+
+ ! Total number of processors we are iterating over, along with
+ ! the current processor and some counters
+ integer :: nproc_x, nproc_y
+ integer :: totalprocs, curproc, x_count, y_count
+
+ ! Number of points per processor (not including the boundary)
+ integer :: xperproc, yperproc
+ integer :: xextra, yextra
+
+ ! The initial boundaries for the domain (not counting the
+ ! halo points)
+ integer, allocatable, dimension(:) :: begin_i
+ integer, allocatable, dimension(:) :: begin_j
+ integer, allocatable, dimension(:) :: end_i
+ integer, allocatable, dimension(:) :: end_j
+
+ ! Error handling
+ character(len=*), parameter :: routine = 'decompose_domain'
+ integer :: alloc_status, dealloc_status
+
+ ! If we are only doing single processor I/O, it doesn't matter
+ ! how many processors there are - we treat it as a single
+ ! decomposition
+ if (num_io_procs .gt. 0) then
+ nproc_x = 1
+ nproc_y = 1
+ else
+ nproc_x = dom_proc_x
+ nproc_y = dom_proc_y
+ endif
+ totalprocs = nproc_x * nproc_y
+
+ ! This sequence of allocations could potentially be offloaded to
+ ! the initialize_decomposition subroutine but as the arrays in
+ ! question are local and temporary, I've kept them here.
+ allocate( begin_i(totalprocs), stat=alloc_status )
+ call check_alloc_status(alloc_status, routine, source, revision, &
+ revdate, 'begin_i')
+ allocate( begin_j(totalprocs), stat=alloc_status )
+ call check_alloc_status(alloc_status, routine, source, revision, &
+ revdate, 'begin_j')
+ allocate( end_i(totalprocs), stat=alloc_status )
+ call check_alloc_status(alloc_status, routine, source, revision, &
+ revdate, 'end_i')
+ allocate( end_j(totalprocs), stat=alloc_status )
+ call check_alloc_status(alloc_status, routine, source, revision, &
+ revdate, 'end_j')
+
+ ! Get the processor-point counts
+ xperproc = (grid%pts_x - 2) / nproc_x
+ yperproc = (grid%pts_y - 2) / nproc_y
+ xextra = (grid%pts_x - 2) - (xperproc * nproc_x)
+ yextra = (grid%pts_y - 2) - (yperproc * nproc_y)
+
+ do y_count = 1,nproc_y
+ do x_count = 1,nproc_x
+
+ ! 2D grid -> 1D number
+ curproc = (y_count - 1) * nproc_x + x_count
+
+ ! X points
+ if (x_count .eq. 1) then
+ begin_i(curproc) = 2
+ else
+ begin_i(curproc) = end_i(curproc-1) + 1
+ endif
+
+ end_i(curproc) = begin_i(curproc) + xperproc - 1
+ if (x_count .gt. nproc_x-xextra) then
+ end_i(curproc) = end_i(curproc) + 1
+ endif
+
+ ! Y Points
+ if (y_count .eq. 1) then
+ begin_j(curproc) = 2
+ else
+ begin_j(curproc) = end_j(curproc-nproc_x) + 1
+ endif
+ end_j(curproc) = begin_j(curproc) + yperproc - 1
+ if (y_count .gt. nproc_y-yextra) then
+ end_j(curproc) = end_j(curproc) + 1
+ endif
+
+ ! Take the boundaries into consideration
+ grid%iminf(curproc) = begin_i(curproc) - nbound
+ grid%imaxf(curproc) = end_i(curproc) + nbound
+ grid%jminf(curproc) = begin_j(curproc) - nbound
+ grid%jmaxf(curproc) = end_j(curproc) + nbound
+ enddo
+ enddo
+
+ deallocate( begin_i, begin_j, end_i, end_j, stat=dealloc_status )
+ call check_dealloc_status(dealloc_status, routine, source, &
+ revision, revdate, '[ij](min|max)f')
+ end subroutine decompose_domain
+
+ ! get_grid_iminf
+ ! --------------
+ ! Accessor method for the index of a subgrid's left extent
+ ! PARAMETERS
+ ! IN grid coamps_grid to get data from
+ ! IN subgrid_num decomposed subgrid number
+ ! OUT iminf i-index in large domain for the
+ ! left edge of this subgrid
+ subroutine get_grid_iminf(grid, subgrid_num, iminf)
+ type(coamps_grid), intent(in) :: grid
+ integer, intent(in) :: subgrid_num
+ integer, intent(out) :: iminf
+
+ iminf = grid%iminf(subgrid_num)
+ end subroutine
+
+ ! get_grid_imaxf
+ ! --------------
+ ! Accessor method for the index of a subgrid's right extent
+ ! PARAMETERS
+ ! IN grid coamps_grid to get data from
+ ! IN subgrid_num decomposed subgrid number
+ ! OUT imaxf i-index in large domain for the
+ ! right edge of this subgrid
+ subroutine get_grid_imaxf(grid, subgrid_num, imaxf)
+ type(coamps_grid), intent(in) :: grid
+ integer, intent(in) :: subgrid_num
+ integer, intent(out) :: imaxf
+
+ imaxf = grid%imaxf(subgrid_num)
+ end subroutine
+
+ ! get_grid_jminf
+ ! --------------
+ ! Accessor method for the index of a subgrid's lower extent
+ ! PARAMETERS
+ ! IN grid coamps_grid to get data from
+ ! IN subgrid_num decomposed subgrid number
+ ! OUT jminf j-index in large domain for the
+ ! lower edge of this subgrid
+ subroutine get_grid_jminf(grid, subgrid_num, jminf)
+ type(coamps_grid), intent(in) :: grid
+ integer, intent(in) :: subgrid_num
+ integer, intent(out) :: jminf
+
+ jminf = grid%jminf(subgrid_num)
+ end subroutine
+
+ ! get_grid_jmaxf
+ ! --------------
+ ! Accessor method for the index of a subgrid's upper extent
+ ! PARAMETERS
+ ! IN grid coamps_grid to get data from
+ ! IN subgrid_num decomposed subgrid number
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list