[Dart-dev] DART/branches Revision: 12243

[dart at ucar.edu]

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