[Dart-dev] DART/branches Revision: 12243
dart at ucar.edu
dart at ucar.edu
Tue Dec 26 11:52:31 MST 2017
thoar at ucar.edu
2017-12-26 11:52:30 -0700 (Tue, 26 Dec 2017)
213
Adding the initial refactoring of Jimmy's interpolation code.
This is largely untested.
When it is tested, the routines should be added to the model_mod
to be incorporated into the model_interpolate() routine.
Added: DART/branches/openggcm/models/openggcm/interp.f90
===================================================================
--- DART/branches/openggcm/models/openggcm/interp.f90 (rev 0)
+++ DART/branches/openggcm/models/openggcm/interp.f90 2017-12-26 18:52:30 UTC (rev 12243)
@@ -0,0 +1,265 @@
+! DART software - Copyright 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$
+
+!> development/test code for interpolation in the CTIM/IPE O+ (and other ions/e-) grid.
+!> Jimmy Raeder, UNH, August 2017
+!>
+!> grid description:
+!> 1. a regular lat/lon grid of footpoints in GEO lat/lon grid
+!> 2. from each footpoint, a field line (dipole only) emerges, which has an
+!> unevenly spaced number of points, however, the spacing of the points along
+!> the field line is the same for all field lines hz(iz).
+!> 3. The dipolar field lines are more regular in SM (no azimuth component), but
+!> that does not help much.
+!> 4. The grid is not completely irregular, but made of 'bricks', i.e., 4 field
+!> lines making a rectangle at the base will make a brick with the 8 points
+!> given at hz(iz) and hz(iz+1). Such grids are common in FEM (Finite Element)
+!> analysis. Such grids are called hexahedron (hex8) element bricks.
+!>
+!> possible simplifications
+!> 1. assume field lines are radial: may work ok over polar cap, but rapidly gets
+!> worse at lower latitudes. Specifically, at mid latitude the horizontal
+!> displacement error becomes roughly equal the height (~600 km) --> not acceptable.
+!> 2. transform everything to SM. Then the azimuthal coordinate is evenly spaced,
+!> leaving us basically only with an irregular grid in latitude and height.
+!> Simplifies the problem somewhat, but we loose generality. In the future
+!> this may change, for example when better coordinates are used, i.e., IGRF/AGGCM
+!> coordinates. Since there is still irregularity we may not save much
+!> computationally --> dismissed.
+!>
+!> Considerations:
+!> 1. We are now having a hex8 grid. Interpolation functions (speficically Lagrange
+!> interpolation) exist for hex8 grids. These functions are defined for
+!> isoparametric coordinates, i.e., the unit brick [-1,1]x[-1,1]x[-1,1].
+!> The mapping from the unit brick into x,y,z space is easy, but not linear, so
+!> the inverse is not easy to come by. This is fine for FEM but not for us,
+!> because we need the inverse mapping (x,y,z)-->(xi,yi,zi), where the latter are
+!> the isoparametric coordinates. --> dismiss
+!> 2. One can instead divide the bricks into tetrahedra (tet4). For tet4, the
+!> isoparametric mapping is linear, so the inverse is easy. The unit tet4 has
+!> the 4 corners (0,0,0) (1,0,0) (0,1,0) and (0,0,1). It is trivial to find out
+!> if a point lies within a tet4: map (x,y,z) --> (xi,yi,zi), then the condition
+!> is xi>0 && yi>0 && zi>0 && xi+yi+zi<1. Interpolation is easily done with
+!> the Lagrange interpolation functions N_j(xi,yi,zi), often called
+!> 'shape functions' in the FEM literature.
+!> 3. Search in tet4 grids is easy and fast, with some pre-computed values
+!> (center, list of neighbors), i.e., descending distance, which scales as
+!> ~(N)**(1/3) where N is the total number of tet4. Also, if the current search
+!> point is close to the previous (as in LOS integration), the search terminates
+!> quickly.
+!> 4. tet4 grids are quite universal, even for completely irregular grids, once the
+!> decomposition is found.
+!> 5. We should do this in (x,y,z) space as opposed to (r,phi,the) spherical
+!> coordinates since DA is all about distance.
+!> 6. The division of bricks into tetrahedra is not unique, see:
+!> https://www.ics.uci.edu/~eppstein/projects/tetra/
+!> Division into 5 tet4 leads to non-conformity: the faces and edges of tet4s
+!> do not line up between bricks. The interpolation is then not continuous
+!> across brick faces. The division into 6 tet4 does not have that problem,
+!> so we use that.
+!>
+!> In this test code I minimize any efforts re/ memory allocation. The grid will
+!> be read in a simple as possible (ascii). The purpose is to demonstrate the
+!> algorithm and test accuracy.
+
+program interp
+
+use types_mod, only : r4, r8
+
+use utilities_mod, only : register_module, initialize_utilities, &
+ find_namelist_in_file, check_namelist_read, &
+ error_handler, E_ERR, E_MSG, nmlfileunit, &
+ do_nml_file, do_nml_term, logfileunit, &
+ finalize_utilities
+
+use openggcm_interp_mod, only : g_oplus_pre, g_oplus_int, nsearch
+
+implicit none
+
+! version controlled file description for error handling, do not edit
+character(len=*), parameter :: source = &
+ "$URL$"
+character(len=*), parameter :: revision = "$Revision$"
+character(len=*), parameter :: revdate = "$Date$"
+
+real(r8), allocatable :: opl_grid_geo_lon(:,:,:)
+real(r8), allocatable :: opl_grid_geo_lat(:,:,:)
+real(r8), allocatable :: opl_grid_geo_hgt(:,:,:)
+real(r8), allocatable :: u(:,:,:)
+
+real(r8) :: a(3,3), b(3,3)
+
+integer :: nphi ! number of longitudes
+integer :: nthe ! number of latitudes
More information about the Dart-dev
mailing list