[Dart-dev] DART/branches Revision: 12637

[dart at ucar.edu]

thoar at ucar.edu
2018-06-06 14:51:20 -0600 (Wed, 06 Jun 2018)
149
Doug sent an email on 5 June 2018:
to me, Jimmy, Nancy 
I believe that the interpolation routine (openggcm_interp_mod.f90) is working properly now.




Modified: DART/branches/openggcm/models/openggcm/interp.f90
===================================================================
--- DART/branches/openggcm/models/openggcm/interp.f90	2018-06-06 16:57:33 UTC (rev 12636)
+++ DART/branches/openggcm/models/openggcm/interp.f90	2018-06-06 20:51:20 UTC (rev 12637)
@@ -75,7 +75,7 @@
                              do_nml_file, do_nml_term, logfileunit,      &
                              finalize_utilities 
 
-use  openggcm_interp_mod, only : g_oplus_pre, g_oplus_int, g_oplus_int_pnt, nsearch
+use  openggcm_interp_mod, only : g_oplus_pre, g_oplus_int, nsearch
 
 implicit none
 

Modified: DART/branches/openggcm/models/openggcm/openggcm_interp_mod.f90
===================================================================
--- DART/branches/openggcm/models/openggcm/openggcm_interp_mod.f90	2018-06-06 16:57:33 UTC (rev 12636)
+++ DART/branches/openggcm/models/openggcm/openggcm_interp_mod.f90	2018-06-06 20:51:20 UTC (rev 12637)
@@ -34,7 +34,7 @@
 ! these are not known a-priory, but could be adjusted later
 integer, parameter :: MTET = 80  !... max num of neighbors for each tet
 integer, parameter :: MPOI = 80  !... max num of tets for each point
-integer, parameter :: MAXDIM = 100  !... max of each dimension
+!integer, parameter :: MAXDIM = 100  !... max of each dimension
 
 !..... x,y,z grid point coord
 real(r8) :: hx(NPOI_MAX),hy(NPOI_MAX),hz(NPOI_MAX) 
@@ -41,9 +41,12 @@
 
 !..... tet mapping matrix (location of center and first corner, 
 real(r8) :: tcen(3,NTET_MAX),tlow(3,NTET_MAX),tmap(3,3,NTET_MAX)
+real(r8) :: tcen_rtp(3,NTET_MAX)
 
 !..... 4 corners + neicount + neighbors + myself
 integer :: itet(MTET+6,NTET_MAX) 
+!..... 4 face neighbors (1:x<0,y&z>0; 2:y<0,x&z>0; 3:z<0,x&y>0; 4:x,y,z>0)
+integer :: ifacetet(4,NTET_MAX)
 !..... for each point, count, and which tet4 it belongs to
 integer :: jtet(MPOI+1,NPOI_MAX) 
 !..... remember the ip,it,iz indices
@@ -51,8 +54,8 @@
 !..... number of tet4 tetrahedra, search steps
 integer :: ntet, nsearch
 
-!..... corner id for grid indices
-integer :: ii(MAXDIM,MAXDIM,MAXDIM)
+!!..... corner id for grid indices
+!integer :: ii(MAXDIM,MAXDIM,MAXDIM)
 
 real(r8), parameter :: PI  = 4.0_r8*atan(1.0_r8)
 real(r8), parameter :: RAD = PI/180.0_r8
@@ -63,7 +66,7 @@
 
 public :: g_oplus_pre, &
           g_oplus_int, &
-          g_oplus_int_pnt, &
+          nsearch, &
           convert_to_cartesian
 
 contains
@@ -87,12 +90,14 @@
 
 real(r8) :: a(3,3), b(3,3), det
 real(r8) :: r, p, t
+real(r8) :: xi, xtmp, yi, ytmp, zi, ztmp
 
-!!!!integer :: ii(np,nt,nz) 
+integer :: ii(np,nt,nz) 
 
 integer :: ip, ip1, it, it1, iz, iz1, k, mytest
 integer :: i1, i2, i3, i4, i5, i6, i7, i8
 integer :: j, jt, kc, kt, l, ll, maxnei, ni, nj
+integer :: mm, nn, npnt, cmatch(3), ctot, jtc, ktc, ktclow
 
 if ( module_initialized ) then
    return ! only need to do this once.
@@ -128,9 +133,12 @@
       hy(k)=it
       hz(k)=iz
    else  ! This is the conversion to cartesian coords
-      hx(k) = r*cos(p)*cos(t)
-      hy(k) = r*sin(p)*cos(t)
-      hz(k) = r*sin(t)
+      !hx(k) = r*cos(p)*cos(t)
+      !hy(k) = r*sin(p)*cos(t)
+      !hz(k) = r*sin(t)
+      hx(k) = r
+      hy(k) = t
+      hz(k) = p
    endif
 
    !... save to map grid values later
@@ -186,9 +194,10 @@
       itet(3,ntet)=i8
       itet(4,ntet)=i7 
       !..... tet4 centers
-      tcen(1,ntet)=0.25*(hx(i1)+hx(i4)+hx(i8)+hx(i7)) 
-      tcen(2,ntet)=0.25*(hy(i1)+hy(i4)+hy(i8)+hy(i7)) 
-      tcen(3,ntet)=0.25*(hz(i1)+hz(i4)+hz(i8)+hz(i7)) 
+      tcen_rtp(1,ntet)=0.25*(hx(i1)+hx(i4)+hx(i8)+hx(i7)) 
+      tcen_rtp(2,ntet)=0.25*(hy(i1)+hy(i4)+hy(i8)+hy(i7))