<p><b>ringler@lanl.gov</b> 2010-05-04 14:25:38 -0600 (Tue, 04 May 2010)</p><p><br>
files changed: numerous<br>
<br>
reason: merge trunk to branches through<br>
        svn merge -r 160:240 https://svn-mpas-model.cgd.ucar.edu/trunk/mpas .<br>
<br>
shallow-water test case 5 is bit for bit the same between this branch and the trunk<br>
</p><hr noshade><pre><font color="gray">Modified: branches/lateral_boundary_conditions/Makefile
===================================================================
--- branches/lateral_boundary_conditions/Makefile        2010-05-04 20:03:20 UTC (rev 241)
+++ branches/lateral_boundary_conditions/Makefile        2010-05-04 20:25:38 UTC (rev 242)
@@ -2,7 +2,7 @@
 MODEL_FORMULATION = -DLANL_FORMULATION
 
 #EXPAND_LEVELS = -DEXPAND_LEVELS=26
-#FILE_OFFSET = -DOFFSET64BIT
+FILE_OFFSET = -DOFFSET64BIT
 
 #########################
 # Section for Zoltan TPL
@@ -75,17 +75,40 @@
         &quot;CC = mpicc&quot; \
         &quot;SFC = gfortran&quot; \
         &quot;SCC = gcc&quot; \
-        &quot;FFLAGS = -O3 -m64 -ffree-line-length-none&quot; \
+        &quot;FFLAGS = -O3 -m64 -ffree-line-length-none -fdefault-real-8&quot; \
         &quot;CFLAGS = -O3 -m64&quot; \
         &quot;LDFLAGS = -O3 -m64&quot; \
         &quot;CORE = $(CORE)&quot; \
         &quot;CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) $(EXPAND_LEVELS) -D_MPI -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)&quot; )
 
+g95:
+        ( make all \
+        &quot;FC = mpif90&quot; \
+        &quot;CC = mpicc&quot; \
+        &quot;SFC = g95&quot; \
+        &quot;SCC = gcc&quot; \
+        &quot;FFLAGS = -O3 -ffree-line-length-huge -r8&quot; \
+        &quot;CFLAGS = -O3&quot; \
+        &quot;LDFLAGS = -O3&quot; \
+        &quot;CORE = $(CORE)&quot; \
+        &quot;CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) $(EXPAND_LEVELS) -D_MPI -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)&quot; )
 
+g95-serial:
+        ( make all \
+        &quot;FC = g95&quot; \
+        &quot;CC = gcc&quot; \
+        &quot;SFC = g95&quot; \
+        &quot;SCC = gcc&quot; \
+        &quot;FFLAGS = -O3 -ffree-line-length-huge -r8&quot; \
+        &quot;CFLAGS = -O3&quot; \
+        &quot;LDFLAGS = -O3&quot; \
+        &quot;CORE = $(CORE)&quot; \
+        &quot;CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) $(EXPAND_LEVELS) -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)&quot; )
 
+
 CPPINCLUDES = -I../inc -I$(NETCDF)/include
 FCINCLUDES = -I../inc -I$(NETCDF)/include
-LIBS = -L$(NETCDF)/lib -lnetcdf -L$(BLAS)/lib -L$(LAPACK)/lib -lblas -llapack
+LIBS = -L$(NETCDF)/lib -lnetcdf
 
 RM = rm -f
 CPP = cpp -C -P -traditional

Modified: branches/lateral_boundary_conditions/namelist.input.hyd_atmos
===================================================================
--- branches/lateral_boundary_conditions/namelist.input.hyd_atmos        2010-05-04 20:03:20 UTC (rev 241)
+++ branches/lateral_boundary_conditions/namelist.input.hyd_atmos        2010-05-04 20:25:38 UTC (rev 242)
@@ -13,6 +13,7 @@
    config_v_theta_eddy_visc2 = 0.0
    config_theta_adv_order = 2
    config_scalar_adv_order = 2
+   config_mp_physics = 0
 /
 
 &amp;io

Modified: branches/lateral_boundary_conditions/namelist.input.ocean
===================================================================
--- branches/lateral_boundary_conditions/namelist.input.ocean        2010-05-04 20:03:20 UTC (rev 241)
+++ branches/lateral_boundary_conditions/namelist.input.ocean        2010-05-04 20:25:38 UTC (rev 242)
@@ -1,21 +1,21 @@
 &amp;sw_model
-   config_test_case = 5
+   config_test_case = 0
    config_time_integration = 'RK4'
-   config_dt = 300.0
-   config_ntimesteps = 3000 
-   config_output_interval = 300
-   config_stats_interval = 10 
-   config_visc  = 4.0
+   config_dt = 120.0
+   config_ntimesteps = 30
+   config_output_interval = 30
+   config_stats_interval = 10
+   config_visc  = 100.0
 /
 
 &amp;io
-   config_input_name = grid.nc
-   config_output_name = output.nc
-   config_restart_name = restart.nc
+   config_input_name = 'grid.quad.20km.nc'
+   config_output_name = 'output.quad.20km.nc'
+   config_restart_name = 'restart.nc'
 /
 
 &amp;restart
-   config_restart_interval = 864000
+   config_restart_interval = 3000
    config_do_restart = .false.
    config_restart_time = 1036800.0
 /

Modified: branches/lateral_boundary_conditions/src/Makefile
===================================================================
--- branches/lateral_boundary_conditions/src/Makefile        2010-05-04 20:03:20 UTC (rev 241)
+++ branches/lateral_boundary_conditions/src/Makefile        2010-05-04 20:25:38 UTC (rev 242)
@@ -2,8 +2,8 @@
 
 all: mpas
 
-mpas: reg_includes externals frame dycore drver
-        $(FC) $(LDFLAGS) -o $(CORE)_model.exe driver/*.o -L. -ldycore -lframework $(LIBS)
+mpas: reg_includes externals frame ops dycore drver
+        $(FC) $(LDFLAGS) -o $(CORE)_model.exe driver/*.o -L. -ldycore -lops -lframework $(LIBS)
 
 reg_includes: 
         ( cd registry; make CC=&quot;$(SCC)&quot; )
@@ -16,6 +16,10 @@
         ( cd framework; make all ) 
         ln -sf framework/libframework.a libframework.a
 
+ops: 
+        ( cd operators; make all ) 
+        ln -sf operators/libops.a libops.a
+
 dycore: 
         ( cd core_$(CORE); make all ) 
         ln -sf core_$(CORE)/libdycore.a libdycore.a
@@ -24,10 +28,11 @@
         ( cd driver; make all ) 
 
 clean:
-        $(RM) $(CORE)_model.exe libframework.a libdycore.a
+        $(RM) $(CORE)_model.exe libframework.a libops.a libdycore.a
         ( cd registry; make clean )
         ( cd external; make clean )
         ( cd framework; make clean )
+        ( cd operators; make clean )
         ( cd inc; rm -f *.inc )
         if [ -d core_$(CORE) ] ; then \
            ( cd core_$(CORE); make clean ) \

Modified: branches/lateral_boundary_conditions/src/core_hyd_atmos/Makefile
===================================================================
--- branches/lateral_boundary_conditions/src/core_hyd_atmos/Makefile        2010-05-04 20:03:20 UTC (rev 241)
+++ branches/lateral_boundary_conditions/src/core_hyd_atmos/Makefile        2010-05-04 20:25:38 UTC (rev 242)
@@ -24,5 +24,5 @@
 .F.o:
         $(RM) $@ $*.mod
         $(CPP) $(CPPFLAGS) $(CPPINCLUDES) $&lt; &gt; $*.f90
-        $(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework
+        $(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework -I../operators
         $(RM) $*.f90

Modified: branches/lateral_boundary_conditions/src/core_hyd_atmos/Registry
===================================================================
--- branches/lateral_boundary_conditions/src/core_hyd_atmos/Registry        2010-05-04 20:03:20 UTC (rev 241)
+++ branches/lateral_boundary_conditions/src/core_hyd_atmos/Registry        2010-05-04 20:25:38 UTC (rev 242)
@@ -17,6 +17,7 @@
 namelist integer   sw_model config_scalar_adv_order     2
 namelist logical   sw_model config_positive_definite    false
 namelist logical   sw_model config_monotonic            true
+namelist integer   sw_model config_mp_physics           0
 namelist character io       config_input_name           grid.nc
 namelist character io       config_output_name          output.nc
 namelist character io       config_restart_name         restart.nc
@@ -36,6 +37,7 @@
 dim vertexDegree vertexDegree
 dim FIFTEEN 15
 dim TWENTYONE 21
+dim R3 3
 dim nVertLevels nVertLevels
 #dim nTracers nTracers
 dim nVertLevelsP1 nVertLevels+1
@@ -123,6 +125,9 @@
 var real    ke ( nVertLevels nCells Time ) o ke - -
 var real    pv_vertex ( nVertLevels nVertices Time ) o pv_vertex - -
 var real    pv_cell ( nVertLevels nCells Time ) o pv_cell - -
+var real    uReconstructX ( nVertLevels nCells Time ) o uReconstructX - -
+var real    uReconstructY ( nVertLevels nCells Time ) o uReconstructY - -
+var real    uReconstructZ ( nVertLevels nCells Time ) o uReconstructZ - -
 
 # Other diagnostic variables: neither read nor written to any files
 var real    vh ( nVertLevels nEdges Time ) - vh - -
@@ -152,3 +157,6 @@
 var real    deriv_two ( FIFTEEN TWO nEdges ) o deriv_two - -
 var integer advCells ( TWENTYONE nCells ) - advCells - -
 
+# Arrays required for reconstruction of velocity field
+var real    coeffs_reconstruct ( R3 maxEdges nCells ) - coeffs_reconstruct - -
+

Modified: branches/lateral_boundary_conditions/src/core_hyd_atmos/module_test_cases.F
===================================================================
--- branches/lateral_boundary_conditions/src/core_hyd_atmos/module_test_cases.F        2010-05-04 20:03:20 UTC (rev 241)
+++ branches/lateral_boundary_conditions/src/core_hyd_atmos/module_test_cases.F        2010-05-04 20:25:38 UTC (rev 242)
@@ -28,10 +28,11 @@
          write(0,*) ' need hydrostatic test case configuration, error stop '
          stop
 
-      else if ((config_test_case == 1) .or. (config_test_case == 2)) then
+      else if ((config_test_case == 1) .or. (config_test_case == 2) .or. (config_test_case == 3)) then
          write(0,*) ' Jablonowski and Williamson baroclinic wave test case '
          if (config_test_case == 1) write(0,*) ' no initial perturbation '
          if (config_test_case == 2) write(0,*) ' initial perturbation included '
+         if (config_test_case == 3) write(0,*) ' normal-mode perturbation included '
          block_ptr =&gt; domain % blocklist
          do while (associated(block_ptr))
             call hyd_test_case_1(block_ptr % mesh, block_ptr % time_levs(1) % state, config_test_case)
@@ -69,6 +70,8 @@
       real (kind=RKIND), parameter :: u_perturbation = 1., pert_radius = 0.1, latitude_pert = 40., longitude_pert = 20.
       real (kind=RKIND), parameter :: theta_c = pii/4.0
       real (kind=RKIND), parameter :: lambda_c = 3.0*pii/2.0
+      real (kind=RKIND), parameter :: rh_max = 0.4       ! Maximum relative humidity
+      real (kind=RKIND), parameter :: k_x = 9.           ! Normal mode wave number
 
       real (kind=RKIND), dimension(:), pointer :: rdnu, rdnw, fnm, fnp, dbn, dnu, dnw
       real (kind=RKIND), dimension(:), pointer :: surface_pressure
@@ -79,7 +82,12 @@
       real (kind=RKIND) :: u, v, flux, fluxk, lat1, lat2, eta_v, r_pert, u_pert, lat_pert, lon_pert, r
 
       real (kind=RKIND) :: ptop, p0, phi
+      real (kind=RKIND) :: lon_Edge
 
+      real (kind=RKIND), dimension(grid % nVertLevels, grid % nCells) :: rel_hum, temperature
+      real (kind=RKIND) :: ptmp, es, qvs
+      integer :: iter
+
 !      real (kind=RKIND), dimension(grid % nVertLevelsP1 ) :: hyai, hybi, znu, znw, znwc, znwv, hyam, hybm
 !      real (kind=RKIND), dimension(grid % nVertLevelsP1 ) :: znuc, znuv, bn, divh, dpn, teta, phi
       real (kind=RKIND), dimension(grid % nVertLevels + 1 ) :: hyai, hybi, znu, znw, znwc, znwv, hyam, hybm
@@ -141,13 +149,14 @@
       h =&gt; state % h % array
       scalars =&gt; state % scalars % array
 
-      scalars = 0.
+      scalars(:,:,:) = 0.
 
       p0      = 100000.
       bn (1) = 1.
       znw(1) = 1.
       znwc(1) = 1.
-      znwv(1) = (znwc(1)-.252)*pii/2.
+      !znwv(1) = (znwc(1)-.252)*pii/2.
+      znwv(1) = ((znwc(1)-.252)*pii/2.*p0-ptop)/(p0-ptop)
                 
       if (cam26) then
 
@@ -208,8 +217,10 @@
       !
 
       do k=1,nz1
-        znuv(k  ) = (znuc(k  )-.252)*pii/2.
-        znwv(k+1) = (znwc(k+1)-.252)*pii/2.
+        !znuv(k  ) = (znuc(k  )-.252)*pii/2.
+        !znwv(k+1) = (znwc(k+1)-.252)*pii/2.
+        znuv(k  ) = ((znuc(k  )-.252)*pii/2.*p0-ptop)/(p0-ptop)
+        znwv(k+1) = ((znwc(k+1)-.252)*pii/2.*p0-ptop)/(p0-ptop)
         dnw (k) = znw(k+1)-znw(k)
         rdnw(k) = 1./dnw(k)
         dbn (k) = rdnw(k)*(bn(k+1)-bn(k))
@@ -247,6 +258,11 @@
             r_pert = sphere_distance( grid % latEdge % array (iEdge), grid % lonEdge % array (iEdge), &amp;
                                       lat_pert, lon_pert, 1.)/(pert_radius)
             u_pert = u_perturbation*exp(-r_pert**2)*(lat2-lat1)*a/grid % dvEdge % array(iEdge)
+
+         else if (config_test_case == 3) then
+            lon_Edge = grid % lonEdge % array(iEdge)
+            u_pert = u_perturbation*cos(k_x*(lon_Edge - lon_pert)) &amp;
+                         *(0.5*(lat2-lat1) - 0.125*(sin(4.*lat2) - sin(4.*lat1)))*a/grid % dvEdge % array(iEdge)
          else
             u_pert = 0.0
          end if
@@ -274,14 +290,16 @@
                                          )
       end do
 
-      do k=1,nz1
-         if(znuc(k).ge.eta_t)  then
-            teta(k) = t0*znuc(k)**(rgas*dtdz/gravity)
-         else
-            teta(k) = t0*znuc(k)**(rgas*dtdz/gravity) + delta_t*(eta_t-znuc(k))**5
-         end if
-         write(6,*) ' k, reference t ',k,teta(k)
-      end do
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! To get hydrostatic balance with misture -- soln. 2.
+! original scheme by Jablonowski
+!                            T' = -1./R_d *(p/p_0) * d(phi')/d(eta)
+!                               = -1./R_d * eta * d(phi')/d(eta)
+! soln. 2 -&gt; derive temperature profile from hydrostatic balance with moisture
+! 
+!                           T_v = -1/(1+q_v)*(p/R_d)* d(eta)/d(p_d) * d(phi)/d(eta)
+!                           phi'(k) = phi(k+1) + d(eta)* alpha_pert * d(eta)/d(p_d)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                         
       do iCell=1,grid % nCells
 
@@ -299,8 +317,13 @@
          end do
 
          do k=1,nz1
-
-            theta (k,iCell) = teta(k)+.75*znuc(k)*pii*u0/rgas*sin(znuv(k))    &amp;
+            ptmp = 0.5*(pressure(k,iCell)+pressure(k+1,iCell))
+            if (znuc(k) &gt;= eta_t) then
+               teta(k) = t0*(ptmp/p0)**(rgas*dtdz/gravity)
+            else
+               teta(k) = t0*(ptmp/p0)**(rgas*dtdz/gravity) + delta_t*(eta_t-ptmp/p0)**5
+            end if
+            theta (k,iCell) = teta(k)+.75*ptmp/h(k,iCell)*pii*u0/rgas*sin(znuv(k))    &amp;
                               *sqrt(cos(znuv(k)))*                         &amp;
                                 ((-2.*sin(phi)**6                          &amp;
                                      *(cos(phi)**2+1./3.)+10./63.)         &amp;
@@ -333,9 +356,89 @@
          end do
       end do
                 
+      write(6,*) 'ptop_dry = ',ptop,'  zt_dry = ',geopotential(nz,1)/gravity
+
+      write(6,*) ' full sounding for dry'
+      do k=1,nz1
+         write(6,*) k, geopotential(k,1)/gravity, 0.01*pressure(k,1), theta(k,1), &amp;
+                    theta(k,1)*(pressure(k,1)/p0)**(rgas/cp)
+      end do
+
+!
+!     initialization for moisture 
+!
+      if (config_mp_physics /= 0) then
+
+         do iCell=1,grid % nCells
+            do k=1,nz1
+               ptmp = 0.5*(pressure(k,iCell) + pressure(k+1,iCell))
+               if (ptmp &lt; 50000.) then
+                  rel_hum(k,iCell) = 0.0
+               else
+                  rel_hum(k,iCell) = (1.-((p0-ptmp)/50000.)**1.25)
+               end if
+               rel_hum(k,iCell) = min(rh_max,rel_hum(k,iCell))
+            end do
+         end do
+      else
+         rel_hum(:,:) = 0.
+      end if
+
+
+      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
+      ! iteration 
+      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
+
+      do iter=1,30
+         do iCell=1,grid % nCells 
+   
+            phi = grid % latCell % array (iCell)
+            do k=1,nz1
+               ptmp = 0.5*(pressure(k+1,iCell)+pressure(k,iCell))
+   
+               if(znuc(k) &gt;= eta_t)  then
+                  teta(k) = t0*(ptmp/p0)**(rgas*dtdz/gravity)
+               else
+                  teta(k) = t0*(ptmp/p0)**(rgas*dtdz/gravity) + delta_t*(eta_t-ptmp/p0)**5
+               end if
+   
+               temperature (k,iCell) = teta(k)+.75*ptmp/h(k,iCell)*pii*u0/rgas*sin(znuv(k))    &amp;
+                                 *sqrt(cos(znuv(k)))*                         &amp;
+                                   ((-2.*sin(phi)**6                          &amp;
+                                        *(cos(phi)**2+1./3.)+10./63.)         &amp;
+                                        *2.*u0*cos(znuv(k))**1.5              &amp;
+                                   +(1.6*cos(phi)**3                          &amp;
+                                        *(sin(phi)**2+2./3.)-pii/4.)*a*omega)
+   
+               temperature(k,iCell) = temperature(k,iCell)/(1.+0.61*scalars(index_qv,k,iCell))
+   
+               theta (k,iCell) = temperature(k,iCell)*  &amp;
+                                     (0.5*(pressure(k,iCell)+pressure(k+1,iCell))/p0)**(-rgas/cp)
+               alpha (k,iCell) = (rgas/p0)*theta(k,iCell)*(1.+1.61*scalars(index_qv,k,iCell))* &amp;
+                                     (0.5*(pressure(k,iCell)+pressure(k+1,iCell))/p0)**cvpm
+   
+               if (temperature(k,iCell) &gt; 273.15) then
+                   es  = 1000.*0.6112*exp(17.67*(temperature(k,iCell)-273.15)/(temperature(k,iCell)-29.65))
+               else
+                   es  = 1000.*0.6112*exp(21.8745584*(temperature(k,iCell)-273.16)/(temperature(k,iCell)-7.66))
+               end if
+               qvs = (287.04/461.6)*es/(ptmp-es)
+!               qvs =  380.*exp(17.27*(temperature(k,iCell)-273.)/(temperature(k,iCell)-36.))/ptmp
+   
+               scalars(index_qv,k,iCell) = rel_hum(k,iCell)*qvs
+            end do
+   
+            do k=nz1,1,-1
+               pressure(k,iCell) = pressure(k+1,iCell)-dnw(k)*h(k,iCell)*(1.+scalars(index_qv,k,iCell))
+               geopotential(k,iCell) = geopotential(k+1,iCell)+dnw(k)*h(k,iCell)*alpha(k,iCell)
+            end do
+   
+         end do
+      end do
+
       write(6,*) 'ptop = ',ptop,'  zt = ',geopotential(nz,1)/gravity
 
-      write(6,*) ' full sounding '
+      write(6,*) ' full sounding with moisture'
       do k=1,nz1
          write(6,*) k, geopotential(k,1)/gravity, 0.01*pressure(k,1), theta(k,1), &amp;
                     theta(k,1)*(pressure(k,1)/p0)**(rgas/cp)

Modified: branches/lateral_boundary_conditions/src/core_hyd_atmos/module_time_integration.F
===================================================================
--- branches/lateral_boundary_conditions/src/core_hyd_atmos/module_time_integration.F        2010-05-04 20:03:20 UTC (rev 241)
+++ branches/lateral_boundary_conditions/src/core_hyd_atmos/module_time_integration.F        2010-05-04 20:25:38 UTC (rev 242)
@@ -118,8 +118,8 @@
                                             block % mesh % nVertLevels, block % mesh % nCells, &amp;
                                             block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
            call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % cqu % array(:,:), &amp;
-                                            block % mesh % nVertLevels, block % mesh % nCells, &amp;
-                                            block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+                                            block % mesh % nVertLevels, block % mesh % nEdges, &amp;
+                                            block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
            call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % h % array(:,:), &amp;
                                             block % mesh % nVertLevels, block % mesh % nCells, &amp;
                                             block % parinfo % cellsToSend, block % parinfo % cellsToRecv)

Modified: branches/lateral_boundary_conditions/src/core_ocean/Makefile
===================================================================
--- branches/lateral_boundary_conditions/src/core_ocean/Makefile        2010-05-04 20:03:20 UTC (rev 241)
+++ branches/lateral_boundary_conditions/src/core_ocean/Makefile        2010-05-04 20:25:38 UTC (rev 242)
@@ -2,7 +2,6 @@
 
 OBJS = module_test_cases.o \
        module_time_integration.o \
-       module_vector_reconstruction.o \
        module_global_diagnostics.o \
        mpas_interface.o
 
@@ -13,13 +12,11 @@
 
 module_test_cases.o:
 
-module_time_integration.o: module_vector_reconstruction.o
+module_time_integration.o: 
 
-module_vector_reconstruction.o: 
-
 module_global_diagnostics.o: 
 
-mpas_interface.o: module_global_diagnostics.o module_vector_reconstruction.o module_test_cases.o module_time_integration.o
+mpas_interface.o: module_global_diagnostics.o module_test_cases.o module_time_integration.o
 
 clean:
         $(RM) *.o *.mod libdycore.a
@@ -27,5 +24,5 @@
 .F.o:
         $(RM) $@ $*.mod
         $(CPP) $(CPPFLAGS) $(CPPINCLUDES) $&lt; &gt; $*.f90
-        $(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework
+        $(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework -I../operators
         $(RM) $*.f90

Modified: branches/lateral_boundary_conditions/src/core_ocean/Registry
===================================================================
--- branches/lateral_boundary_conditions/src/core_ocean/Registry        2010-05-04 20:03:20 UTC (rev 241)
+++ branches/lateral_boundary_conditions/src/core_ocean/Registry        2010-05-04 20:25:38 UTC (rev 242)
@@ -80,9 +80,6 @@
 var real    rho ( nVertLevels nCells Time ) iro rho - -
 
 # Arrays required for reconstruction of velocity field
-var real    matrix_reconstruct ( maxEdges maxEdges nCells ) - matrix_reconstruct - -
-var real    normal ( R3 maxEdges nCells ) - normal - -
-var real    rbf_value ( maxEdges nCells ) - rbf_value - -
 var real    coeffs_reconstruct ( R3 maxEdges nCells ) - coeffs_reconstruct - -
 
 # Boundary conditions: read from input, saved in restart and written to output

Deleted: branches/lateral_boundary_conditions/src/core_ocean/module_vector_reconstruction.F
===================================================================
--- branches/lateral_boundary_conditions/src/core_ocean/module_vector_reconstruction.F        2010-05-04 20:03:20 UTC (rev 241)
+++ branches/lateral_boundary_conditions/src/core_ocean/module_vector_reconstruction.F        2010-05-04 20:25:38 UTC (rev 242)
@@ -1,354 +0,0 @@
-module vector_reconstruction
-
-   use grid_types
-   use configure
-   use constants
-
-   implicit none
-
-   public :: init_reconstruct, reconstruct
-
-   contains
-
-   subroutine init_reconstruct(grid)
-   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
-   ! Purpose: reconstruct vector field at cell centers based on radial basis functions 
-   !
-   ! Input: grid meta data and vector component data residing at cell edges
-   !
-   ! Output: reconstructed vector field (measured in X,Y,Z) located at cell centers
-   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
-
-      implicit none
-
-      type (grid_meta), intent(inout) :: grid 
-
-      ! temporary arrays needed in the (to be constructed) init procedure
-      integer :: nCells, nEdges, nVertLevels, nCellsSolve
-      integer, dimension(:,:), pointer :: cellsOnCell, edgesOnCell, cellsOnEdge
-      integer, dimension(:), pointer :: nEdgesOnCell
-      integer :: iEdge, iCell, k, cell1, cell2, EdgeMax, j, i, npoints
-      integer :: lwork, info
-      integer, allocatable, dimension(:) :: ipvt
-      real (kind=RKIND), allocatable, dimension(:) :: work
-      real (kind=RKIND), dimension(:), pointer :: dcEdge, xCell, yCell, zCell
-      real (kind=RKIND) :: r, t, v, X1(3), X2(3), alpha
-      real (kind=RKIND), allocatable, dimension(:,:,:) :: xLoc
-
-      real (kind=RKIND), dimension(:,:,:), pointer :: matrix_reconstruct
-      real (kind=RKIND), dimension(:,:,:), pointer :: normal
-      real (kind=RKIND), dimension(:,:), pointer :: rbf_value
-      real (kind=RKIND), dimension(:,:,:), pointer :: coeffs_reconstruct
-      real (kind=RKIND), allocatable, dimension(:,:) :: mwork
-
-      !========================================================
-      ! arrays filled and saved during init procedure
-      !========================================================
-      matrix_reconstruct =&gt; grid % matrix_reconstruct % array
-      normal =&gt; grid % normal % array
-      rbf_value =&gt; grid % rbf_value % array
-      coeffs_reconstruct =&gt; grid % coeffs_reconstruct % array
-
-      !========================================================
-      ! temporary variables needed for init procedure
-      !========================================================
-      xCell       =&gt; grid % xCell % array
-      yCell       =&gt; grid % yCell % array
-      zCell       =&gt; grid % zCell % array
-      cellsOnCell =&gt; grid % cellsOnCell % array
-      cellsOnEdge =&gt; grid % cellsOnEdge % array
-      edgesOnCell =&gt; grid % edgesOnCell % array
-      nEdgesOnCell=&gt; grid % nEdgesOnCell % array
-      dcEdge      =&gt; grid % dcEdge % array
-      nCells      = grid % nCells
-      nCellsSolve = grid % nCellsSolve
-      nEdges      = grid % nEdges
-      nVertLevels = grid % nVertLevels
-
-      ! allocate arrays
-      EdgeMax = maxval(nEdgesOnCell)
-      allocate(xLoc(3,EdgeMax,nCells))
-      allocate(work(EdgeMax*EdgeMax))
-      allocate(ipvt(EdgeMax))
-       
-       ! init arrays
-       matrix_reconstruct = 0.0
-       normal = 0.0
-       rbf_value = 0.0
-
-       ! loop over all cells to be solved on this block
-       do iCell=1,nCellsSolve
-
-         npoints = nEdgesOnCell(iCell)   ! only loop over the number of edges for this cell
-
-         ! fill normal vector and xLoc arrays
-         ! X1 is the location of the reconstruction, X2 are the neighboring centers, xLoc is the edge positions
-         cell1 = iCell
-         X1(1) = xCell(cell1)
-         X1(2) = yCell(cell1)
-         X1(3) = zCell(cell1)
-         do j=1,npoints
-           cell2 = cellsOnCell(j,iCell)
-           iEdge = edgesOnCell(j,iCell)
-           X2(1) = xCell(cell2)
-           X2(2) = yCell(cell2)
-           X2(3) = zCell(cell2)
-           if(cell2.gt.cell1) then
-               normal(:,j,iCell) = X2(:) - X1(:)
-           else
-               normal(:,j,iCell) = X1(:) - X2(:)
-           endif
-           call unit_vector_in_R3(normal(:,j,iCell))
-           xLoc(:,j,iCell) = 0.5*(X2(:) + X1(:))
-         enddo
-
-         alpha = 0.0
-         do i=1,npoints
-           call get_distance(xLoc(:,i,iCell),X1(:),r)
-           alpha = alpha + r
-         enddo
-         alpha = 4*alpha / npoints
-         do j=1,npoints
-           do i=1,npoints
-              call get_distance(xLoc(:,i,iCell),xLoc(:,j,iCell),r)
-              r = r / alpha
-              call evaluate_rbf(r,t)
-              call get_dotproduct(normal(:,i,iCell),normal(:,j,iCell),v)
-              matrix_reconstruct(i,j,iCell) = v*t
-           enddo
-         enddo
- 
-         ! save value of RBF when evaluated between reconstruction location and edge locations
-         do j=1,npoints
-           call get_distance(xLoc(:,j,iCell), X1(:), r)
-           r = r / alpha
-           call evaluate_rbf(r,t)
-           rbf_value(j,iCell) = t
-         enddo
-
-         ! invert matrix_reconstruct using lapack routines
-         ! xsad 10-02-09:
-         !!!!! the lapack stuff isn't working on coyote in parallel !!!!!!!
-         !allocate(mwork(npoints,npoints))
-         !lwork = npoints*npoints
-         !mwork(:,:) = matrix_reconstruct(1:npoints,1:npoints,iCell)
-         !call dgetrf(npoints, npoints, mwork, npoints, ipvt, info)
-         !if(info.ne.0) then
-          !write(6,*) info, 'dgetrf'
-          !stop
-         !endif
-         !call dgetri(npoints, mwork, npoints, ipvt, work, lwork, info)
-         !matrix_reconstruct(1:npoints,1:npoints,iCell) = mwork(:,:)
-         !if(info.ne.0) then
-          !write(6,*) info, 'dgetri'
-          !stop
-         !endif
-          !deallocate(mwork)
-          ! xsad 10-02-09 end
-
-         ! compute the coefficients for reconstructing uX, uY, uZ at cell centers from u_i normal to edges
-         ! uX = sum_j(coeffs(1,j) * u_j) (similarly for Y,Z)
-         ! coeffs(:,j) = sum_i(rbf_values(i) * normal(:,i) * matrix(i,j))
-         do j=1,npoints
-           coeffs_reconstruct(:,j,iCell) = 0.0
-           do i=1,npoints
-              coeffs_reconstruct(:,j,iCell) = coeffs_reconstruct(:,j,iCell) + rbf_value(i,iCell) * normal(:,i,iCell) &amp;
-                * matrix_reconstruct(i,j,iCell)
-           enddo
-         enddo
-
-      enddo   ! iCell
-
-      deallocate(ipvt)
-      deallocate(work)
-      deallocate(xLoc)
-
-      end subroutine init_reconstruct
-
-
-      subroutine reconstruct(s, grid)
-      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-      ! Purpose: reconstruct vector field at cell centers based on radial basis functions
-      !
-      ! Input: grid meta data and vector component data residing at cell edges
-      !
-      ! Output: reconstructed vector field (measured in X,Y,Z) located at cell centers
-      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-      implicit none
-
-      type (grid_state), intent(inout) :: s 
-      type (grid_meta), intent(in) :: grid
-
-      !   temporary arrays needed in the compute procedure
-      integer :: nCells, nEdges, nVertLevels, nCellsSolve
-      integer, dimension(:,:), pointer :: edgesOnCell
-      integer, dimension(:), pointer :: nEdgesOnCell
-      real (kind=RKIND), dimension(:), pointer :: xCell, yCell, zCell
-      integer :: iCell,iEdge, i, j, k, npoints, EdgeMax
-      real (kind=RKIND) :: r, t, p!, t1(3), t2(3), t3(3), east(3)
-      real (kind=RKIND), dimension(:,:), pointer :: u
-      real (kind=RKIND), dimension(:,:), pointer :: uReconstructX, uReconstructY, uReconstructZ
-      real (kind=RKIND), allocatable, dimension(:,:) :: rhs, c
-
-      real (kind=RKIND), dimension(:,:,:), pointer :: matrix_reconstruct
-      real (kind=RKIND), dimension(:,:,:), pointer :: normal
-      real (kind=RKIND), dimension(:,:), pointer :: rbf_value
-      real (kind=RKIND), dimension(:,:,:), pointer :: coeffs_reconstruct
-
-      ! stored arrays used during compute procedure
-      matrix_reconstruct =&gt; grid % matrix_reconstruct % array
-      normal =&gt; grid % normal % array
-      rbf_value =&gt; grid % rbf_value % array
-      coeffs_reconstruct =&gt; grid % coeffs_reconstruct % array
-
-      ! temporary variables
-      edgesOnCell =&gt; grid % edgesOnCell % array
-      nEdgesOnCell=&gt; grid % nEdgesOnCell % array
-      nCells      = grid % nCells
-      nCellsSolve = grid % nCellsSolve
-      nEdges      = grid % nEdges
-      nVertLevels = grid % nVertLevels
-      xCell =&gt; grid % xCell % array
-      yCell =&gt; grid % yCell % array
-      zCell =&gt; grid % zCell % array
-
-      u =&gt; s % u % array
-      uReconstructX =&gt; s % uReconstructX % array
-      uReconstructY =&gt; s % uReconstructY % array
-      uReconstructZ =&gt; s % uReconstructZ % array
-
-      ! allocate space for temporary arrays
-      EdgeMax = maxval(nEdgesOnCell)
-!      allocate(c(nVertLevels,EdgeMax))
-      allocate(rhs(nVertLevels,EdgeMax))
-
-      ! init the intent(out)
-      uReconstructX = 0.0
-      uReconstructY = 0.0
-      uReconstructZ = 0.0
-
-      ! loop over cell centers
-      do iCell=1,nCellsSolve
-
-        npoints = nEdgesOnCell(iCell)
-
-        ! fill the RHS of the matrix equation
-        ! for testing purposes, fill rhs with test function
-        ! test function assumes uniform flow to the east
-        rhs = 0.0
-        do j=1,npoints
-          iEdge = edgesOnCell(j,iCell)
-          do k=1,nVertLevels
-            rhs(k,j) = u(k,iEdge)
-          enddo
-        enddo
-  
-        ! compute c by multiplying matrix_reconstruct * rhs  (Eq 8 in Tex document)
-!        c = 0.0
-!        do i=1,npoints
-!         do j=1,npoints
-!           do k=1,nVertLevels
-!             c(k,i) = c(k,i) + matrix_reconstruct(i,j,iCell)*rhs(k,j)
-!           enddo
-!         enddo
-!        enddo
-  
-        ! reconstruct the velocity field at point X1 (Eq 6 in Tex document)
-!        do i=1,npoints
-!          do k=1,nVertLevels
-!            uReconstructX(k,iCell) = uReconstructX(k,iCell) + c(k,i)*rbf_value(i,iCell)*normal(1,i,iCell)
-!            uReconstructY(k,iCell) = uReconstructY(k,iCell) + c(k,i)*rbf_value(i,iCell)*normal(2,i,iCell)
-!            uReconstructZ(k,iCell) = uReconstructZ(k,iCell) + c(k,i)*rbf_value(i,iCell)*normal(3,i,iCell)
-!          enddo
-!        enddo
-!
-        ! a more efficient reconstruction where rbf_values*matrix_reconstruct has been precomputed
-        ! in coeffs_reconstruct
-        do i=1,npoints
-          uReconstructX(:,iCell) = uReconstructX(:,iCell) + coeffs_reconstruct(1,i,iCell) * rhs(:,i)
-          uReconstructY(:,iCell) = uReconstructY(:,iCell) + coeffs_reconstruct(2,i,iCell) * rhs(:,i)
-          uReconstructZ(:,iCell) = uReconstructZ(:,iCell) + coeffs_reconstruct(3,i,iCell) * rhs(:,i)
-        enddo
-
-        
-      enddo   ! iCell
-
-      ! deallocate
-      deallocate(rhs)
-!      deallocate(c)
- 
-   end subroutine reconstruct
-
-   subroutine get_distance(x1,x2,xout)
-     implicit none
-     real(kind=RKIND), intent(in) :: x1(3), x2(3)
-     real(kind=RKIND), intent(out) :: xout
-     xout = sqrt( (x1(1)-x2(1))**2 + (x1(2)-x2(2))**2 + (x1(3)-x2(3))**2 )
-   end subroutine get_distance
- 
-   subroutine get_dotproduct(x1,x2,xout)
-     implicit none
-     real(kind=RKIND), intent(in) :: x1(3), x2(3)
-     real(kind=RKIND), intent(out) :: xout
-     xout = x1(1)*x2(1) + x1(2)*x2(2) + x1(3)*x2(3)
-   end subroutine get_dotproduct
- 
- 
-   subroutine unit_vector_in_R3(xin)
-     implicit none
-     real (kind=RKIND), intent(inout) :: xin(3)
-     real (kind=RKIND) :: mag
-     mag = sqrt(xin(1)**2+xin(2)**2+xin(3)**2)
-     xin(:) = xin(:) / mag
-   end subroutine unit_vector_in_R3
-
-
-   subroutine evaluate_rbf(r,xout)
-     real(kind=RKIND), intent(in) ::  r
-     real(kind=RKIND), intent(out) :: xout
-
-     ! Gaussian
-       xout = exp(-r**2)
-
-     ! multiquadrics
-     ! xout = 1.0 / sqrt(1.0**2 + r**2)
-
-     ! other
-     ! xout = 1.0 / (1.0**2 + r**2)
-
-   end subroutine evaluate_rbf
-
-!======================================================================
-! BEGINNING OF CROSS_PRODUCT_IN_R3
-!======================================================================
-        subroutine cross_product_in_R3(p_1,p_2,p_out)
-
-!-----------------------------------------------------------------------
-! PURPOSE: compute p_1 cross p_2 and place in p_out
-!-----------------------------------------------------------------------
-
-!-----------------------------------------------------------------------
-! intent(in)
-!-----------------------------------------------------------------------
-        real (kind=RKIND), intent(in) ::                            &amp;
-                        p_1 (3),                                      &amp;
-                        p_2 (3)
-
-!-----------------------------------------------------------------------
-! intent(out)
-!-----------------------------------------------------------------------
-        real (kind=RKIND), intent(out) ::                           &amp;
-                        p_out (3)
-
-        p_out(1) = p_1(2)*p_2(3)-p_1(3)*p_2(2)
-        p_out(2) = p_1(3)*p_2(1)-p_1(1)*p_2(3)
-        p_out(3) = p_1(1)*p_2(2)-p_1(2)*p_2(1)
-
-        end subroutine cross_product_in_R3
-!======================================================================
-! END OF CROSS_PRODUCT_IN_R3
-!======================================================================
-
-
-end module vector_reconstruction

Modified: branches/lateral_boundary_conditions/src/core_sw/Makefile
===================================================================
--- branches/lateral_boundary_conditions/src/core_sw/Makefile        2010-05-04 20:03:20 UTC (rev 241)
+++ branches/lateral_boundary_conditions/src/core_sw/Makefile        2010-05-04 20:25:38 UTC (rev 242)
@@ -1,7 +1,6 @@
 .SUFFIXES: .F .o
 
 OBJS = module_test_cases.o \
-       module_vector_reconstruction.o \
        module_time_integration.o \
        mpas_interface.o
 
@@ -12,17 +11,15 @@
 
 module_test_cases.o: 
 
-module_vector_reconstruction.o:
+module_time_integration.o: 
 
-module_time_integration.o: module_vector_reconstruction.o
+mpas_interface.o: module_test_cases.o module_time_integration.o 
 
-mpas_interface.o: module_test_cases.o module_time_integration.o module_vector_reconstruction.o
-
 clean:
         $(RM) *.o *.mod libdycore.a
 
 .F.o:
         $(RM) $@ $*.mod
         $(CPP) $(CPPFLAGS) $(CPPINCLUDES) $&lt; &gt; $*.f90
-        $(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework
+        $(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework -I../operators
         $(RM) $*.f90

Modified: branches/lateral_boundary_conditions/src/core_sw/Registry
===================================================================
--- branches/lateral_boundary_conditions/src/core_sw/Registry        2010-05-04 20:03:20 UTC (rev 241)
+++ branches/lateral_boundary_conditions/src/core_sw/Registry        2010-05-04 20:25:38 UTC (rev 242)
@@ -78,9 +78,7 @@
 var real    h_s ( nCells ) iro h_s - -
 
 # Arrays required for reconstruction of velocity field
-var real    matrix_reconstruct ( maxEdges maxEdges nCells ) - matrix_reconstruct - -
-var real    normal ( R3 maxEdges nCells ) - normal - -
-var real    rbf_value ( maxEdges nCells ) - rbf_value - -
+var real    coeffs_reconstruct ( R3 maxEdges nCells ) - coeffs_reconstruct - -
 
 # Boundary conditions: read from input, saved in restart and written to output
 var integer boundaryEdge ( nVertLevels nEdges ) iro boundaryEdge - -

Modified: branches/lateral_boundary_conditions/src/core_sw/module_test_cases.F
===================================================================
--- branches/lateral_boundary_conditions/src/core_sw/module_test_cases.F        2010-05-04 20:03:20 UTC (rev 241)
+++ branches/lateral_boundary_conditions/src/core_sw/module_test_cases.F        2010-05-04 20:25:38 UTC (rev 242)
@@ -508,7 +508,7 @@
       real (kind=RKIND), intent(in) :: theta
 
       AA = 0.5 * w * (2.0 * omega + w) * cos(theta)**2.0 + &amp;
-          0.25 * K**2.0 * cos(theta)**(2.0*R) * ((R+1.0)*cos(theta)**2.0 + 2.0*R**2.0 - R - 2.0 - 2.0*(R*cos(theta))**2.0)
+          0.25 * K**2.0 * cos(theta)**(2.0*R) * ((R+1.0)*cos(theta)**2.0 + 2.0*R**2.0 - R - 2.0 - 2.0*R**2.0 * cos(theta)**(-2.0))
 
    end function AA
 

Modified: branches/lateral_boundary_conditions/src/core_sw/module_time_integration.F
===================================================================
--- branches/lateral_boundary_conditions/src/core_sw/module_time_integration.F        2010-05-04 20:03:20 UTC (rev 241)
+++ branches/lateral_boundary_conditions/src/core_sw/module_time_integration.F        2010-05-04 20:25:38 UTC (rev 242)
@@ -1,6 +1,6 @@
 module time_integration
 
-!  use vector_reconstruction
+   use vector_reconstruction
    use grid_types
    use configure
    use constants
@@ -214,7 +214,7 @@
 
          call compute_solve_diagnostics(dt, block % time_levs(2) % state, block % mesh)
 
-!        call reconstruct(block % time_levs(2) % state, block % mesh)
+         call reconstruct(block % time_levs(2) % state, block % mesh)
 
          block =&gt; block % next
       end do

Deleted: branches/lateral_boundary_conditions/src/core_sw/module_vector_reconstruction.F
===================================================================
--- branches/lateral_boundary_conditions/src/core_sw/module_vector_reconstruction.F        2010-05-04 20:03:20 UTC (rev 241)
+++ branches/lateral_boundary_conditions/src/core_sw/module_vector_reconstruction.F        2010-05-04 20:25:38 UTC (rev 242)
@@ -1,286 +0,0 @@
-module vector_reconstruction
-
-   use grid_types
-   use configure
-   use constants
-
-   public :: init_reconstruct, reconstruct
-
-   contains
-
-   subroutine init_reconstruct(grid)
-   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
-   ! Purpose: reconstruct vector field at cell centers based on radial basis functions 
-   !
-   ! Input: grid meta data and vector component data residing at cell edges
-   !
-   ! Output: reconstructed vector field (measured in X,Y,Z) located at cell centers
-   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
-
-      implicit none
-
-      type (grid_meta), intent(inout) :: grid 
-
-      ! temporary arrays needed in the (to be constructed) init procedure
-      integer :: nCells, nEdges, nVertLevels, nCellsSolve
-      integer, dimension(:,:), pointer :: cellsOnCell, edgesOnCell, cellsOnEdge
-      integer, dimension(:), pointer :: nEdgesOnCell
-      integer :: iEdge, iCell, k, cell1, cell2, EdgeMax, j, i, npoints
-      integer :: lwork, info
-      integer, allocatable, dimension(:) :: ipvt
-      real (kind=RKIND), allocatable, dimension(:) :: work
-      real (kind=RKIND), dimension(:), pointer :: dcEdge, xCell, yCell, zCell
-      real (kind=RKIND) :: r, t, v, X1(3), X2(3), alpha
-      real (kind=RKIND), allocatable, dimension(:,:,:) :: xLoc
-
-      real (kind=RKIND), dimension(:,:,:), pointer :: matrix_reconstruct
-      real (kind=RKIND), dimension(:,:,:), pointer :: normal
-      real (kind=RKIND), dimension(:,:), pointer :: rbf_value
-      real (kind=RKIND), allocatable, dimension(:,:) :: mwork
-
-      !========================================================
-      ! arrays filled and saved during init procedure
-      !========================================================
-      matrix_reconstruct =&gt; grid % matrix_reconstruct % array
-      normal =&gt; grid % normal % array
-      rbf_value =&gt; grid % rbf_value % array
-
-      !========================================================
-      ! temporary variables needed for init procedure
-      !========================================================
-      xCell       =&gt; grid % xCell % array
-      yCell       =&gt; grid % yCell % array
-      zCell       =&gt; grid % zCell % array
-      cellsOnCell =&gt; grid % cellsOnCell % array
-      cellsOnEdge =&gt; grid % cellsOnEdge % array
-      edgesOnCell =&gt; grid % edgesOnCell % array
-      nEdgesOnCell=&gt; grid % nEdgesOnCell % array
-      dcEdge      =&gt; grid % dcEdge % array
-      nCells      = grid % nCells
-      nCellsSolve = grid % nCellsSolve
-      nEdges      = grid % nEdges
-      nVertLevels = grid % nVertLevels
-
-      ! allocate arrays
-      EdgeMax = maxval(nEdgesOnCell)
-      allocate(xLoc(3,EdgeMax,nCells))
-      allocate(work(EdgeMax*EdgeMax))
-      allocate(ipvt(EdgeMax))
-       
-       ! init arrays
-       matrix_reconstruct = 0.0
-       normal = 0.0
-       rbf_value = 0.0
-
-       ! loop over all cells to be solved on this block
-       do iCell=1,nCellsSolve
-
-         ! fill normal vector and xLoc arrays
-         ! X1 is the location of the reconstruction, X2 are the neighboring centers, xLoc is the edge positions
-         cell1 = iCell
-         X1(1) = xCell(cell1)
-         X1(2) = yCell(cell1)
-         X1(3) = zCell(cell1)
-         do j=1,nEdgesOnCell(iCell)
-           cell2 = cellsOnCell(j,iCell)
-           iEdge = edgesOnCell(j,iCell)
-           X2(1) = xCell(cell2)
-           X2(2) = yCell(cell2)
-           X2(3) = zCell(cell2)
-           if(cell2.gt.cell1) then
-               normal(:,j,iCell) = X2(:) - X1(:)
-           else
-               normal(:,j,iCell) = X1(:) - X2(:)
-           endif
-           call unit_vector_in_R3(normal(:,j,iCell))
-           xLoc(:,j,iCell) = 0.5*(X2(:) + X1(:))
-         enddo
-
-         npoints = nEdgesOnCell(iCell)   ! only loop over the number of edges for this cell
-         alpha = 0.0
-         do i=1,npoints
-           call get_distance(xLoc(:,i,iCell),X1(:),r)
-           alpha = alpha + r
-         enddo
-         alpha = alpha / npoints
-         do j=1,npoints
-           do i=1,npoints
-              call get_distance(xLoc(:,i,iCell),xLoc(:,j,iCell),r)
-              r = r / alpha
-              call evaluate_rbf(r,t)
-              call get_dotproduct(normal(:,i,iCell),normal(:,j,iCell),v)
-              matrix_reconstruct(i,j,iCell) = v*t
-           enddo
-          enddo
- 
-         ! save value of RBF when evaluated between reconstruction location and edge locations
-         do j=1,npoints
-           call get_distance(xLoc(:,j,iCell), X1(:), r)
-           r = r / alpha
-           call evaluate_rbf(r,t)
-           rbf_value(j,iCell) = t
-         enddo
- 
-         ! invert matrix_reconstruct using lapack routines
-         allocate(mwork(npoints,npoints))
-         lwork = npoints*npoints
-         mwork(:,:) = matrix_reconstruct(1:npoints,1:npoints,iCell)
-         call dgetrf(npoints, npoints, mwork, npoints, ipvt, info)
-         if(info.ne.0) then
-          write(6,*) info, 'dgetrf'
-          stop
-         endif
-         call dgetri(npoints, mwork, npoints, ipvt, work, lwork, info)
-         matrix_reconstruct(1:npoints,1:npoints,iCell) = mwork(:,:)
-         if(info.ne.0) then
-          write(6,*) info, 'dgetri'
-          stop
-         endif
-          deallocate(mwork)
-
-      enddo   ! iCell
-
-      deallocate(ipvt)
-      deallocate(work)
-      deallocate(xLoc)
-
-      end subroutine init_reconstruct
-
-
-      subroutine reconstruct(s, grid)
-      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-      ! Purpose: reconstruct vector field at cell centers based on radial basis functions
-      !
-      ! Input: grid meta data and vector component data residing at cell edges
-      !
-      ! Output: reconstructed vector field (measured in X,Y,Z) located at cell centers
-      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-      implicit none
-
-      type (grid_state), intent(inout) :: s 
-      type (grid_meta), intent(in) :: grid
-
-      !   temporary arrays needed in the compute procedure
-      integer :: nCells, nEdges, nVertLevels, nCellsSolve
-      integer, dimension(:,:), pointer :: edgesOnCell
-      integer, dimension(:), pointer :: nEdgesOnCell
-      integer :: iCell,iEdge, i, j, k, npoints, EdgeMax
-      real (kind=RKIND) :: r
-      real (kind=RKIND), dimension(:,:), pointer :: u
-      real (kind=RKIND), dimension(:,:), pointer :: uReconstructX, uReconstructY, uReconstructZ
-      real (kind=RKIND), allocatable, dimension(:,:) :: rhs, c
-
-      real (kind=RKIND), dimension(:,:,:), pointer :: matrix_reconstruct
-      real (kind=RKIND), dimension(:,:,:), pointer :: normal
-      real (kind=RKIND), dimension(:,:), pointer :: rbf_value
-
-      ! stored arrays used during compute procedure
-      matrix_reconstruct =&gt; grid % matrix_reconstruct % array
-      normal =&gt; grid % normal % array
-      rbf_value =&gt; grid % rbf_value % array
-
-      ! temporary variables
-      edgesOnCell =&gt; grid % edgesOnCell % array
-      nEdgesOnCell=&gt; grid % nEdgesOnCell % array
-      nCells      = grid % nCells
-      nCellsSolve = grid % nCellsSolve
-      nEdges      = grid % nEdges
-      nVertLevels = grid % nVertLevels
-      u =&gt; s % u % array
-      uReconstructX =&gt; s % uReconstructX % array
-      uReconstructY =&gt; s % uReconstructY % array
-      uReconstructZ =&gt; s % uReconstructZ % array
-
-      ! allocate space for temporary arrays
-      EdgeMax = maxval(nEdgesOnCell)
-      allocate(c(nVertLevels,EdgeMax))
-      allocate(rhs(nVertLevels,EdgeMax))
-
-      ! init the intent(out)
-      uReconstructX = 0.0
-      uReconstructY = 0.0
-      uReconstructZ = 0.0
-
-      ! loop over cell centers
-      do iCell=1,nCellsSolve
-        npoints = nEdgesOnCell(iCell)
-
-        ! fill the RHS of the matrix equation
-        ! for testing purposes, fill rhs with test function
-        ! test function assumes uniform flow in direction of normal vector of edge 1
-        rhs = 0.0
-        do j=1,npoints
-          iEdge = edgesOnCell(j,iCell)
-          do k=1,nVertLevels
-            rhs(k,j) = u(k,iEdge)
-          enddo
-        enddo
-  
-        ! compute c by multiplying matrix_reconstruct * rhs  (Eq 8 in Tex document)
-        c = 0.0
-        do i=1,npoints
-         do j=1,npoints
-           do k=1,nVertLevels
-             c(k,i) = c(k,i) + matrix_reconstruct(i,j,iCell)*rhs(k,j)
-           enddo
-         enddo
-        enddo
-  
-        ! reconstruct the velocity field at point X1 (Eq 6 in Tex document)
-        do i=1,npoints
-          do k=1,nVertLevels
-            uReconstructX(k,iCell) = uReconstructX(k,iCell) + c(k,i)*rbf_value(i,iCell)*normal(1,i,iCell)
-            uReconstructY(k,iCell) = uReconstructY(k,iCell) + c(k,i)*rbf_value(i,iCell)*normal(2,i,iCell)
-            uReconstructZ(k,iCell) = uReconstructZ(k,iCell) + c(k,i)*rbf_value(i,iCell)*normal(3,i,iCell)
-          enddo
-        enddo
-
-      enddo   ! iCell
-
-      ! deallocate
-      deallocate(rhs)
-      deallocate(c)
-
-   end subroutine reconstruct
-
-   subroutine get_distance(x1,x2,xout)
-     implicit none
-     real(kind=RKIND), intent(in) :: x1(3), x2(3)
-     real(kind=RKIND), intent(out) :: xout
-     xout = sqrt( (x1(1)-x2(1))**2 + (x1(2)-x2(2))**2 + (x1(3)-x2(3))**2 )
-   end subroutine get_distance
- 
-   subroutine get_dotproduct(x1,x2,xout)
-     implicit none
-     real(kind=RKIND), intent(in) :: x1(3), x2(3)
-     real(kind=RKIND), intent(out) :: xout
-     xout = x1(1)*x2(1) + x1(2)*x2(2) + x1(3)*x2(3)
-   end subroutine get_dotproduct
- 
- 
-   subroutine unit_vector_in_R3(xin)
-     implicit none
-     real (kind=RKIND), intent(inout) :: xin(3)
-     real (kind=RKIND) :: mag
-     mag = sqrt(xin(1)**2+xin(2)**2+xin(3)**2)
-     xin(:) = xin(:) / mag
-   end subroutine unit_vector_in_R3
-
-
-   subroutine evaluate_rbf(xin,xout)
-     real(kind=RKIND), intent(in) ::  xin
-     real(kind=RKIND), intent(out) :: xout
-
-     ! Gaussian
-     ! xout = exp(-r**2)
-
-     ! multiquadrics
-       xout = 1.0 / sqrt(1.0**2 + r**2)
-
-     ! other
-     ! xout = 1.0 / (1.0**2 + r**2)
-
-   end subroutine evaluate_rbf
-
-end module vector_reconstruction

Modified: branches/lateral_boundary_conditions/src/core_sw/mpas_interface.F
===================================================================
--- branches/lateral_boundary_conditions/src/core_sw/mpas_interface.F        2010-05-04 20:03:20 UTC (rev 241)
+++ branches/lateral_boundary_conditions/src/core_sw/mpas_interface.F        2010-05-04 20:25:38 UTC (rev 242)
@@ -2,7 +2,7 @@
 
    use grid_types
    use test_cases
-!  use vector_reconstruction
+   use vector_reconstruction
 
    implicit none
 
@@ -25,7 +25,8 @@
    real (kind=RKIND), intent(in) :: dt
 
    call compute_solve_diagnostics(dt, block % time_levs(1) % state, mesh)
-!  call init_reconstruct(mesh)
+   call init_reconstruct(mesh)
+   call reconstruct(block % time_levs(1) % state, mesh)
 
 end subroutine mpas_init
 

Modified: branches/lateral_boundary_conditions/src/framework/module_dmpar.F
===================================================================
--- branches/lateral_boundary_conditions/src/framework/module_dmpar.F        2010-05-04 20:03:20 UTC (rev 241)
+++ branches/lateral_boundary_conditions/src/framework/module_dmpar.F        2010-05-04 20:25:38 UTC (rev 242)
@@ -110,6 +110,25 @@
    end subroutine dmpar_abort
 
 
+   subroutine dmpar_global_abort(mesg)
+
+      implicit none
+
+      character (len=*), intent(in) :: mesg
+
+#ifdef _MPI
+      integer :: mpi_ierr, mpi_errcode
+
+      write(0,*) trim(mesg)
+      call MPI_Abort(MPI_COMM_WORLD, mpi_errcode, mpi_ierr)
+#endif
+
+      write(0,*) trim(mesg)
+      stop
+
+   end subroutine dmpar_global_abort
+
+
    subroutine dmpar_bcast_int(dminfo, i)
 
       implicit none

Modified: branches/lateral_boundary_conditions/src/framework/module_io_input.F
===================================================================
--- branches/lateral_boundary_conditions/src/framework/module_io_input.F        2010-05-04 20:03:20 UTC (rev 241)
+++ branches/lateral_boundary_conditions/src/framework/module_io_input.F        2010-05-04 20:25:38 UTC (rev 242)
@@ -1001,10 +1001,10 @@
 
       integer, dimension(:), pointer :: super_int1d
       integer, dimension(:,:), pointer :: super_int2d
-      real :: super_real0d
-      real, dimension(:), pointer :: super_real1d
-      real, dimension(:,:), pointer :: super_real2d
-      real, dimension(:,:,:), pointer :: super_real3d
+      real (kind=RKIND) :: super_real0d
+      real (kind=RKIND), dimension(:), pointer :: super_real1d
+      real (kind=RKIND), dimension(:,:), pointer :: super_real2d
+      real (kind=RKIND), dimension(:,:,:), pointer :: super_real3d
 
       integer :: k
 
@@ -1039,6 +1039,17 @@
 #else
       nferr = nf_open(trim(input_obj % filename), NF_SHARE, input_obj % rd_ncid)
 #endif
+
+      if (nferr /= NF_NOERR) then
+         write(0,*) ' '
+         if (config_do_restart) then
+            write(0,*) 'Error opening restart file ''', trim(input_obj % filename), ''''
+         else
+            write(0,*) 'Error opening input file ''', trim(input_obj % filename), ''''
+         end if
+         write(0,*) ' '
+         call dmpar_abort(dminfo)
+      end if
  
 #include &quot;netcdf_read_ids.inc&quot;
 

Modified: branches/lateral_boundary_conditions/src/framework/module_io_output.F
===================================================================
--- branches/lateral_boundary_conditions/src/framework/module_io_output.F        2010-05-04 20:03:20 UTC (rev 241)
+++ branches/lateral_boundary_conditions/src/framework/module_io_output.F        2010-05-04 20:25:38 UTC (rev 242)
@@ -126,9 +126,9 @@
       integer, dimension(:), pointer :: super_int1d
       integer, dimension(:,:), pointer :: super_int2d
       real :: super_real0d
-      real, dimension(:), pointer :: super_real1d
-      real, dimension(:,:), pointer :: super_real2d
-      real, dimension(:,:,:), pointer :: super_real3d
+      real (kind=RKIND), dimension(:), pointer :: super_real1d
+      real (kind=RKIND), dimension(:,:), pointer :: super_real2d
+      real (kind=RKIND), dimension(:,:,:), pointer :: super_real3d
 
       output_obj % time = itime
 

Copied: branches/lateral_boundary_conditions/src/operators (from rev 240, trunk/mpas/src/operators)

Deleted: branches/lateral_boundary_conditions/src/operators/Makefile
===================================================================
--- trunk/mpas/src/operators/Makefile        2010-05-04 15:59:06 UTC (rev 240)
+++ branches/lateral_boundary_conditions/src/operators/Makefile        2010-05-04 20:25:38 UTC (rev 242)
@@ -1,19 +0,0 @@
-.SUFFIXES: .F .o
-
-OBJS = module_vector_reconstruction.o
-
-all: operators
-
-operators: $(OBJS)
-        ar -ru libops.a $(OBJS)
-
-module_vector_reconstruction.o:
-
-clean:
-        $(RM) *.o *.mod libops.a
-
-.F.o:
-        $(RM) $@ $*.mod
-        $(CPP) $(CPPFLAGS) $(CPPINCLUDES) $&lt; &gt; $*.f90
-        $(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework
-        $(RM) $*.f90

Copied: branches/lateral_boundary_conditions/src/operators/Makefile (from rev 240, trunk/mpas/src/operators/Makefile)
===================================================================
--- branches/lateral_boundary_conditions/src/operators/Makefile                                (rev 0)
+++ branches/lateral_boundary_conditions/src/operators/Makefile        2010-05-04 20:25:38 UTC (rev 242)
@@ -0,0 +1,19 @@
+.SUFFIXES: .F .o
+
+OBJS = module_vector_reconstruction.o
+
+all: operators
+
+operators: $(OBJS)
+        ar -ru libops.a $(OBJS)
+
+module_vector_reconstruction.o:
+
+clean:
+        $(RM) *.o *.mod libops.a
+
+.F.o:
+        $(RM) $@ $*.mod
+        $(CPP) $(CPPFLAGS) $(CPPINCLUDES) $&lt; &gt; $*.f90
+        $(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework
+        $(RM) $*.f90

Deleted: branches/lateral_boundary_conditions/src/operators/module_vector_reconstruction.F
===================================================================
--- trunk/mpas/src/operators/module_vector_reconstruction.F        2010-05-04 15:59:06 UTC (rev 240)
+++ branches/lateral_boundary_conditions/src/operators/module_vector_reconstruction.F        2010-05-04 20:25:38 UTC (rev 242)
@@ -1,454 +0,0 @@
-module vector_reconstruction
-
-   use grid_types
-   use configure
-   use constants
-
-   implicit none
-
-   public :: init_reconstruct, reconstruct
-
-   contains
-
-   subroutine init_reconstruct(grid)
-   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
-   ! Purpose: pre-compute coefficients used by the reconstruct() routine
-   !
-   ! Input: grid meta data
-   !
-   ! Output: grid % coeffs_reconstruct - coefficients used to reconstruct 
-   !                                     velocity vectors at cell centers 
-   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
-
-      implicit none
-
-      type (grid_meta), intent(inout) :: grid 
-
-      ! temporary arrays needed in the (to be constructed) init procedure
-      integer :: nCells, nEdges, nVertLevels, nCellsSolve
-      integer, dimension(:,:), pointer :: cellsOnCell, edgesOnCell, cellsOnEdge
-      integer, dimension(:), pointer :: nEdgesOnCell
-      integer :: iEdge, iCell, k, cell1, cell2, EdgeMax, j, i, npoints, matrixSize
-      integer :: lwork, info
-      integer, allocatable, dimension(:) :: pivotingIndices
-      real (kind=RKIND), dimension(:), pointer :: dcEdge, xCell, yCell, zCell
-      real (kind=RKIND) :: r, rbfValue, v, X1(3), X2(3), alpha, rHat(3), normalDotRHat
-      real (kind=RKIND) :: xPlane, yPlane, xNormalPlane, yNormalPlane, xHatPlane(3), yHatPlane(3)
-      real (kind=RKIND), allocatable, dimension(:,:,:) :: xLoc
-
-      real (kind=RKIND), dimension(:,:,:), pointer :: coeffs_reconstruct
-      real (kind=RKIND), allocatable, dimension(:,:) :: mwork
-      real (kind=RKIND), dimension(:,:), pointer :: matrix, invMatrix
-      real (kind=RKIND), dimension(:,:), pointer :: normals
-
-      !========================================================
-      ! arrays filled and saved during init procedure
-      !========================================================
-      coeffs_reconstruct =&gt; grid % coeffs_reconstruct % array
-
-      !========================================================
-      ! temporary variables needed for init procedure
-      !========================================================
-      xCell       =&gt; grid % xCell % array
-      yCell       =&gt; grid % yCell % array
-      zCell       =&gt; grid % zCell % array
-      cellsOnCell =&gt; grid % cellsOnCell % array
-      cellsOnEdge =&gt; grid % cellsOnEdge % array
-      edgesOnCell =&gt; grid % edgesOnCell % array
-      nEdgesOnCell=&gt; grid % nEdgesOnCell % array
-      dcEdge      =&gt; grid % dcEdge % array
-      nCells      = grid % nCells
-      nCellsSolve = grid % nCellsSolve
-      nEdges      = grid % nEdges
-      nVertLevels = grid % nVertLevels
-
-      ! allocate arrays
-      EdgeMax = maxval(nEdgesOnCell)
-      allocate(xLoc(3,EdgeMax,nCells))
-
-      allocate(normals(3,EdgeMax))
-       
-      ! init arrays
-      coeffs_reconstruct = 0.0
-      normals = 0
-
-      ! loop over all cells to be solved on this block
-      do iCell=1,nCellsSolve
-
-         ! fill normal vector and xLoc arrays
-         ! X1 is the location of the reconstruction, X2 are the neighboring centers, 
-         !  xLoc is the edge positions
-         cell1 = iCell
-         X1(1) = xCell(cell1)
-         X1(2) = yCell(cell1)
-         X1(3) = zCell(cell1)
-
-         rHat = X1
-         call unit_vector_in_R3(rHat)
-
-         do j=1,nEdgesOnCell(iCell)
-           iEdge = edgesOnCell(j,iCell)
-           if (iCell == cellsOnEdge(1,iEdge)) then
-               cell2 = cellsOnEdge(2,iEdge)
-               X2(1) = xCell(cell2)
-               X2(2) = yCell(cell2)
-               X2(3) = zCell(cell2)
-               normals(:,j) = X2(:) - X1(:)
-           else
-               cell2 = cellsOnEdge(1,iEdge)
-               X2(1) = xCell(cell2)
-               X2(2) = yCell(cell2)
-               X2(3) = zCell(cell2)
-               normals(:,j) = X1(:) - X2(:)
-           endif
-
-           call unit_vector_in_R3(normals(:,j))
-
-           xLoc(:,j,iCell) = 0.5*(X2(:) + X1(:))
-         enddo
-
-         npoints = nEdgesOnCell(iCell)   ! only loop over the number of edges for this cell
-         matrixSize = npoints+2 ! room for 2 vector components for constant flow
-                                !  basis functions
-         allocate(matrix(matrixSize,matrixSize))
-         matrix = 0.0
-         alpha = 0.0
-         do i=1,npoints
-           call get_distance(xLoc(:,i,iCell),X1(:),r)
-           alpha = alpha + r
-         enddo
-         alpha = alpha / npoints
-         do j=1,npoints
-           do i=1,npoints
-              call get_distance(xLoc(:,i,iCell),xLoc(:,j,iCell),r)
-              r = r / alpha
-              call evaluate_rbf(r,rbfValue)
-              call get_dotproduct(normals(:,i),normals(:,j),v)
-              matrix(i,j) = v*rbfValue
-           enddo
-         enddo
-
-         ! add matrix entries to suppoert constant flow
-         ! xHat and yHat are a local basis in the plane of the horizontal cell
-         ! we arbitrarily choose xHat to point toward the first edge
-         call get_dotproduct(normals(:,1),rHat,normalDotRHat)
-         xHatPlane = normals(:,1) - normalDotRHat*rHat(:)
-         call unit_vector_in_R3(xHatPlane)
-
-         call cross_product_in_R3(rHat, xHatPlane, yHatPlane)
-         call unit_vector_in_R3(yHatPlane) ! just to be sure...
-
-         do j=1,npoints
-           call get_dotproduct(normals(:,j),xHatPlane, xNormalPlane)
-           call get_dotproduct(normals(:,j),yHatPlane, yNormalPlane)
-           matrix(j,npoints+1) = xNormalPlane
-           matrix(j,npoints+2) = yNormalPlane
-           matrix(npoints+1,j) = matrix(j,npoints+1)
-           matrix(npoints+2,j) = matrix(j,npoints+2)
-         end do
-       
- 
-         ! invert matrix
-         allocate(invMatrix(matrixSize,matrixSize))
-         allocate(pivotingIndices(matrixSize))
-         invMatrix = 0.0
-         pivotingIndices = 0
-         call migs(matrix,matrixSize,invMatrix,pivotingIndices)
-
-         ! compute the coefficients for reconstructing uX, uY, uZ at cell centers from 
-         ! u_i normal to edges
-         ! uX = sum_j(coeffs(1,j) * u_j) (similarly for Y,Z)
-         ! coeffs(:,j) = sum_i(rbf_values(i) * normal(:,i) * matrix(i,j))
-         do i=1,npoints
-           ! compute value of RBF when evaluated between reconstruction location and edge locations
-           call get_distance(xLoc(:,i,iCell), X1(:), r)
-           r = r / alpha
-           call evaluate_rbf(r,rbfValue)
-
-           ! project the normals onto tangent plane of the cell
-           ! normal = normal - (normal dot r_hat) r_hat
-           ! rHat, the unit vector pointing from the domain center to the cell center
-           call get_dotproduct(normals(:,i),rHat,normalDotRHat)
-           normals(:,i) = normals(:,i) - normalDotRHat*rHat(:)
-
-           do j=1,npoints
-              coeffs_reconstruct(:,j,iCell) = coeffs_reconstruct(:,j,iCell) &amp;
-                + rbfValue * normals(:,i) * invMatrix(i,j)
-           enddo
-         enddo
-         ! polynomial parts
-         do j=1,npoints
-            coeffs_reconstruct(:,j,iCell) = coeffs_reconstruct(:,j,iCell) &amp;
-              + invMatrix(npoints+1,j)*xHatPlane
-            coeffs_reconstruct(:,j,iCell) = coeffs_reconstruct(:,j,iCell) &amp;
-              + invMatrix(npoints+2,j)*yHatPlane
-         enddo
-
-         deallocate(matrix)
-         deallocate(invMatrix)
-         deallocate(pivotingIndices)
-
-      enddo   ! iCell
-
-      deallocate(xLoc)
-      deallocate(normals)
-
-   end subroutine init_reconstruct
-
-
-   subroutine reconstruct(s, grid)
-   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-   ! Purpose: reconstruct vector field at cell centers based on radial basis functions
-   !
-   ! Input: grid meta data and vector component data residing at cell edges
-   !
-   ! Output: reconstructed vector field (measured in X,Y,Z) located at cell centers
-   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-      implicit none
-
-      type (grid_state), intent(inout) :: s 
-      type (grid_meta), intent(in) :: grid
-
-      !   temporary arrays needed in the compute procedure
-      integer :: nCellsSolve
-      integer, dimension(:,:), pointer :: edgesOnCell
-      integer, dimension(:), pointer :: nEdgesOnCell
-      integer :: iCell,iEdge, i
-      real (kind=RKIND), dimension(:,:), pointer :: u
-      real (kind=RKIND), dimension(:,:), pointer :: uReconstructX, uReconstructY, uReconstructZ
-
-      real (kind=RKIND), dimension(:,:,:), pointer :: coeffs_reconstruct
-
-      ! stored arrays used during compute procedure
-      coeffs_reconstruct =&gt; grid % coeffs_reconstruct % array
-
-      ! temporary variables
-      edgesOnCell =&gt; grid % edgesOnCell % array
-      nEdgesOnCell=&gt; grid % nEdgesOnCell % array
-      nCellsSolve = grid % nCellsSolve
-      u =&gt; s % u % array
-      uReconstructX =&gt; s % uReconstructX % array
-      uReconstructY =&gt; s % uReconstructY % array
-      uReconstructZ =&gt; s % uReconstructZ % array
-
-      ! init the intent(out)
-      uReconstructX = 0.0
-      uReconstructY = 0.0
-      uReconstructZ = 0.0
-
-      ! loop over cell centers
-      do iCell=1,nCellsSolve
-        ! a more efficient reconstruction where rbf_values*matrix_reconstruct has been precomputed
-        ! in coeffs_reconstruct
-        do i=1,nEdgesOnCell(iCell)
-          iEdge = edgesOnCell(i,iCell)
-          uReconstructX(:,iCell) = uReconstructX(:,iCell) &amp;
-            + coeffs_reconstruct(1,i,iCell) * u(:,iEdge)
-          uReconstructY(:,iCell) = uReconstructY(:,iCell) &amp;
-            + coeffs_reconstruct(2,i,iCell) * u(:,iEdge)
-          uReconstructZ(:,iCell) = uReconstructZ(:,iCell) &amp;
-            + coeffs_reconstruct(3,i,iCell) * u(:,iEdge)
-        enddo
-      enddo   ! iCell
-
-   end subroutine reconstruct
-
-   subroutine get_distance(x1,x2,xout)
-     implicit none
-     real(kind=RKIND), intent(in) :: x1(3), x2(3)
-     real(kind=RKIND), intent(out) :: xout
-     xout = sqrt( (x1(1)-x2(1))**2 + (x1(2)-x2(2))**2 + (x1(3)-x2(3))**2 )
-   end subroutine get_distance
- 
-   subroutine get_dotproduct(x1,x2,xout)
-     implicit none
-     real(kind=RKIND), intent(in) :: x1(3), x2(3)
-     real(kind=RKIND), intent(out) :: xout
-     xout = x1(1)*x2(1) + x1(2)*x2(2) + x1(3)*x2(3)
-   end subroutine get_dotproduct
- 
- 
-   subroutine unit_vector_in_R3(xin)
-     implicit none
-     real (kind=RKIND), intent(inout) :: xin(3)
-     real (kind=RKIND) :: mag
-     mag = sqrt(xin(1)**2+xin(2)**2+xin(3)**2)
-     xin(:) = xin(:) / mag
-   end subroutine unit_vector_in_R3
-
-
-   subroutine evaluate_rbf(xin,xout)
-     real(kind=RKIND), intent(in) ::  xin
-     real(kind=RKIND), intent(out) :: xout
-
-     ! Gaussian
-     ! xout = exp(-r**2)
-
-     ! multiquadrics
-       xout = 1.0 / sqrt(1.0**2 + xin**2)
-
-     ! other
-     ! xout = 1.0 / (1.0**2 + r**2)
-
-   end subroutine evaluate_rbf
-
-!======================================================================
-! BEGINNING OF CROSS_PRODUCT_IN_R3
-!======================================================================
-        subroutine cross_product_in_R3(p_1,p_2,p_out)
-
-!-----------------------------------------------------------------------
-! PURPOSE: compute p_1 cross p_2 and place in p_out
-!-----------------------------------------------------------------------
-
-!-----------------------------------------------------------------------
-! intent(in)
-!-----------------------------------------------------------------------
-        real (kind=RKIND), intent(in) ::                            &amp;
-                        p_1 (3),                                      &amp;
-                        p_2 (3)
-
-!-----------------------------------------------------------------------
-! intent(out)
-!-----------------------------------------------------------------------
-        real (kind=RKIND), intent(out) ::                           &amp;
-                        p_out (3)
-
-        p_out(1) = p_1(2)*p_2(3)-p_1(3)*p_2(2)
-        p_out(2) = p_1(3)*p_2(1)-p_1(1)*p_2(3)
-        p_out(3) = p_1(1)*p_2(2)-p_1(2)*p_2(1)
-
-        end subroutine cross_product_in_R3
-!======================================================================
-! END OF CROSS_PRODUCT_IN_R3
-!======================================================================
-
-! Updated 10/24/2001.
-!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!   Program 4.4   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!                                                                       !
-! Please Note:                                                          !
-!                                                                       !
-! (1) This computer program is written by Tao Pang in conjunction with  !
-!     his book, &quot;An Introduction to Computational Physics,&quot; published   !
-!     by Cambridge University Press in 1997.                            !
-!                                                                       !
-! (2) No warranties, express or implied, are made for this program.     !
-!                                                                       !
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!
-SUBROUTINE MIGS (A,N,X,INDX)
-!
-! Subroutine to invert matrix A(N,N) with the inverse stored
-! in X(N,N) in the output.  Copyright (c) Tao Pang 2001.
-!
-  IMPLICIT NONE
-  INTEGER, INTENT (IN) :: N
-  INTEGER :: I,J,K
-  INTEGER, INTENT (OUT), DIMENSION (N) :: INDX
-  REAL (kind=RKIND), INTENT (INOUT), DIMENSION (N,N):: A
-  REAL (kind=RKIND), INTENT (OUT), DIMENSION (N,N):: X
-  REAL (kind=RKIND), DIMENSION (N,N) :: B
-!
-  DO I = 1, N
-    DO J = 1, N
-      B(I,J) = 0.0
-    END DO
-  END DO
-  DO I = 1, N
-    B(I,I) = 1.0
-  END DO
-!
-  CALL ELGS (A,N,INDX)
-!
-  DO I = 1, N-1
-    DO J = I+1, N
-      DO K = 1, N
-        B(INDX(J),K) = B(INDX(J),K)-A(INDX(J),I)*B(INDX(I),K)
-      END DO
-    END DO
-  END DO
-!
-  DO I = 1, N
-    X(N,I) = B(INDX(N),I)/A(INDX(N),N)
-    DO J = N-1, 1, -1
-      X(J,I) = B(INDX(J),I)
-      DO K = J+1, N
-        X(J,I) = X(J,I)-A(INDX(J),K)*X(K,I)
-      END DO
-      X(J,I) =  X(J,I)/A(INDX(J),J)
-    END DO
-  END DO
-END SUBROUTINE MIGS
-
-
-SUBROUTINE ELGS (A,N,INDX)
-!
-! Subroutine to perform the partial-pivoting Gaussian elimination.
-! A(N,N) is the original matrix in the input and transformed matrix
-! plus the pivoting element ratios below the diagonal in the output.
-! INDX(N) records the pivoting order.  Copyright (c) Tao Pang 2001.
-!
-  IMPLICIT NONE
-  INTEGER, INTENT (IN) :: N
-  INTEGER :: I,J,K,ITMP
-  INTEGER, INTENT (OUT), DIMENSION (N) :: INDX
-  REAL (kind=RKIND) :: C1,PI,PI1,PJ
-  REAL (kind=RKIND), INTENT (INOUT), DIMENSION (N,N) :: A
-  REAL (kind=RKIND), DIMENSION (N) :: C
-!
-! Initialize the index
-!
-  DO I = 1, N
-    INDX(I) = I
-  END DO
-!
-! Find the rescaling factors, one from each row
-!
-  DO I = 1, N
-    C1= 0.0
-    DO J = 1, N
-      !C1 = AMAX1(C1,ABS(A(I,J)))
-      C1 = MAX(C1,ABS(A(I,J)))
-    END DO
-    C(I) = C1
-  END DO
-!
-! Search the pivoting (largest) element from each column
-!
-  DO J = 1, N-1
-    PI1 = 0.0
-    DO I = J, N
-      PI = ABS(A(INDX(I),J))/C(INDX(I))
-      IF (PI.GT.PI1) THEN
-        PI1 = PI
-        K   = I
-      ENDIF
-    END DO
-!
-! Interchange the rows via INDX(N) to record pivoting order
-!
-    ITMP    = INDX(J)
-    INDX(J) = INDX(K)
-    INDX(K) = ITMP
-    DO I = J+1, N
-      PJ  = A(INDX(I),J)/A(INDX(J),J)
-!
-! Record pivoting ratios below the diagonal
-!
-      A(INDX(I),J) = PJ
-!
-! Modify other elements accordingly
-!
-      DO K = J+1, N
-        A(INDX(I),K) = A(INDX(I),K)-PJ*A(INDX(J),K)
-      END DO
-    END DO
-  END DO
-!
-END SUBROUTINE ELGS
-
-end module vector_reconstruction

Copied: branches/lateral_boundary_conditions/src/operators/module_vector_reconstruction.F (from rev 240, trunk/mpas/src/operators/module_vector_reconstruction.F)
===================================================================
--- branches/lateral_boundary_conditions/src/operators/module_vector_reconstruction.F                                (rev 0)
+++ branches/lateral_boundary_conditions/src/operators/module_vector_reconstruction.F        2010-05-04 20:25:38 UTC (rev 242)
@@ -0,0 +1,454 @@
+module vector_reconstruction
+
+   use grid_types
+   use configure
+   use constants
+
+   implicit none
+
+   public :: init_reconstruct, reconstruct
+
+   contains
+
+   subroutine init_reconstruct(grid)
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
+   ! Purpose: pre-compute coefficients used by the reconstruct() routine
+   !
+   ! Input: grid meta data
+   !
+   ! Output: grid % coeffs_reconstruct - coefficients used to reconstruct 
+   !                                     velocity vectors at cell centers 
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
+
+      implicit none
+
+      type (grid_meta), intent(inout) :: grid 
+
+      ! temporary arrays needed in the (to be constructed) init procedure
+      integer :: nCells, nEdges, nVertLevels, nCellsSolve
+      integer, dimension(:,:), pointer :: cellsOnCell, edgesOnCell, cellsOnEdge
+      integer, dimension(:), pointer :: nEdgesOnCell
+      integer :: iEdge, iCell, k, cell1, cell2, EdgeMax, j, i, npoints, matrixSize
+      integer :: lwork, info
+      integer, allocatable, dimension(:) :: pivotingIndices
+      real (kind=RKIND), dimension(:), pointer :: dcEdge, xCell, yCell, zCell
+      real (kind=RKIND) :: r, rbfValue, v, X1(3), X2(3), alpha, rHat(3), normalDotRHat
+      real (kind=RKIND) :: xPlane, yPlane, xNormalPlane, yNormalPlane, xHatPlane(3), yHatPlane(3)
+      real (kind=RKIND), allocatable, dimension(:,:,:) :: xLoc
+
+      real (kind=RKIND), dimension(:,:,:), pointer :: coeffs_reconstruct
+      real (kind=RKIND), allocatable, dimension(:,:) :: mwork
+      real (kind=RKIND), dimension(:,:), pointer :: matrix, invMatrix
+      real (kind=RKIND), dimension(:,:), pointer :: normals
+
+      !========================================================
+      ! arrays filled and saved during init procedure
+      !========================================================
+      coeffs_reconstruct =&gt; grid % coeffs_reconstruct % array
+
+      !========================================================
+      ! temporary variables needed for init procedure
+      !========================================================
+      xCell       =&gt; grid % xCell % array
+      yCell       =&gt; grid % yCell % array
+      zCell       =&gt; grid % zCell % array
+      cellsOnCell =&gt; grid % cellsOnCell % array
+      cellsOnEdge =&gt; grid % cellsOnEdge % array
+      edgesOnCell =&gt; grid % edgesOnCell % array
+      nEdgesOnCell=&gt; grid % nEdgesOnCell % array
+      dcEdge      =&gt; grid % dcEdge % array
+      nCells      = grid % nCells
+      nCellsSolve = grid % nCellsSolve
+      nEdges      = grid % nEdges
+      nVertLevels = grid % nVertLevels
+
+      ! allocate arrays
+      EdgeMax = maxval(nEdgesOnCell)
+      allocate(xLoc(3,EdgeMax,nCells))
+
+      allocate(normals(3,EdgeMax))
+       
+      ! init arrays
+      coeffs_reconstruct = 0.0
+      normals = 0
+
+      ! loop over all cells to be solved on this block
+      do iCell=1,nCellsSolve
+
+         ! fill normal vector and xLoc arrays
+         ! X1 is the location of the reconstruction, X2 are the neighboring centers, 
+         !  xLoc is the edge positions
+         cell1 = iCell
+         X1(1) = xCell(cell1)
+         X1(2) = yCell(cell1)
+         X1(3) = zCell(cell1)
+
+         rHat = X1
+         call unit_vector_in_R3(rHat)
+
+         do j=1,nEdgesOnCell(iCell)
+           iEdge = edgesOnCell(j,iCell)
+           if (iCell == cellsOnEdge(1,iEdge)) then
+               cell2 = cellsOnEdge(2,iEdge)
+               X2(1) = xCell(cell2)
+               X2(2) = yCell(cell2)
+               X2(3) = zCell(cell2)
+               normals(:,j) = X2(:) - X1(:)
+           else
+               cell2 = cellsOnEdge(1,iEdge)
+               X2(1) = xCell(cell2)
+               X2(2) = yCell(cell2)
+               X2(3) = zCell(cell2)
+               normals(:,j) = X1(:) - X2(:)
+           endif
+
+           call unit_vector_in_R3(normals(:,j))
+
+           xLoc(:,j,iCell) = 0.5*(X2(:) + X1(:))
+         enddo
+
+         npoints = nEdgesOnCell(iCell)   ! only loop over the number of edges for this cell
+         matrixSize = npoints+2 ! room for 2 vector components for constant flow
+                                !  basis functions
+         allocate(matrix(matrixSize,matrixSize))
+         matrix = 0.0
+         alpha = 0.0
+         do i=1,npoints
+           call get_distance(xLoc(:,i,iCell),X1(:),r)
+           alpha = alpha + r
+         enddo
+         alpha = alpha / npoints
+         do j=1,npoints
+           do i=1,npoints
+              call get_distance(xLoc(:,i,iCell),xLoc(:,j,iCell),r)
+              r = r / alpha
+              call evaluate_rbf(r,rbfValue)
+              call get_dotproduct(normals(:,i),normals(:,j),v)
+              matrix(i,j) = v*rbfValue
+           enddo
+         enddo
+
+         ! add matrix entries to suppoert constant flow
+         ! xHat and yHat are a local basis in the plane of the horizontal cell
+         ! we arbitrarily choose xHat to point toward the first edge
+         call get_dotproduct(normals(:,1),rHat,normalDotRHat)
+         xHatPlane = normals(:,1) - normalDotRHat*rHat(:)
+         call unit_vector_in_R3(xHatPlane)
+
+         call cross_product_in_R3(rHat, xHatPlane, yHatPlane)
+         call unit_vector_in_R3(yHatPlane) ! just to be sure...
+
+         do j=1,npoints
+           call get_dotproduct(normals(:,j),xHatPlane, xNormalPlane)
+           call get_dotproduct(normals(:,j),yHatPlane, yNormalPlane)
+           matrix(j,npoints+1) = xNormalPlane
+           matrix(j,npoints+2) = yNormalPlane
+           matrix(npoints+1,j) = matrix(j,npoints+1)
+           matrix(npoints+2,j) = matrix(j,npoints+2)
+         end do
+       
+ 
+         ! invert matrix
+         allocate(invMatrix(matrixSize,matrixSize))
+         allocate(pivotingIndices(matrixSize))
+         invMatrix = 0.0
+         pivotingIndices = 0
+         call migs(matrix,matrixSize,invMatrix,pivotingIndices)
+
+         ! compute the coefficients for reconstructing uX, uY, uZ at cell centers from 
+         ! u_i normal to edges
+         ! uX = sum_j(coeffs(1,j) * u_j) (similarly for Y,Z)
+         ! coeffs(:,j) = sum_i(rbf_values(i) * normal(:,i) * matrix(i,j))
+         do i=1,npoints
+           ! compute value of RBF when evaluated between reconstruction location and edge locations
+           call get_distance(xLoc(:,i,iCell), X1(:), r)
+           r = r / alpha
+           call evaluate_rbf(r,rbfValue)
+
+           ! project the normals onto tangent plane of the cell
+           ! normal = normal - (normal dot r_hat) r_hat
+           ! rHat, the unit vector pointing from the domain center to the cell center
+           call get_dotproduct(normals(:,i),rHat,normalDotRHat)
+           normals(:,i) = normals(:,i) - normalDotRHat*rHat(:)
+
+           do j=1,npoints
+              coeffs_reconstruct(:,j,iCell) = coeffs_reconstruct(:,j,iCell) &amp;
+                + rbfValue * normals(:,i) * invMatrix(i,j)
+           enddo
+         enddo
+         ! polynomial parts
+         do j=1,npoints
+            coeffs_reconstruct(:,j,iCell) = coeffs_reconstruct(:,j,iCell) &amp;
+              + invMatrix(npoints+1,j)*xHatPlane
+            coeffs_reconstruct(:,j,iCell) = coeffs_reconstruct(:,j,iCell) &amp;
+              + invMatrix(npoints+2,j)*yHatPlane
+         enddo
+
+         deallocate(matrix)
+         deallocate(invMatrix)
+         deallocate(pivotingIndices)
+
+      enddo   ! iCell
+
+      deallocate(xLoc)
+      deallocate(normals)
+
+   end subroutine init_reconstruct
+
+
+   subroutine reconstruct(s, grid)
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+   ! Purpose: reconstruct vector field at cell centers based on radial basis functions
+   !
+   ! Input: grid meta data and vector component data residing at cell edges
+   !
+   ! Output: reconstructed vector field (measured in X,Y,Z) located at cell centers
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+      implicit none
+
+      type (grid_state), intent(inout) :: s 
+      type (grid_meta), intent(in) :: grid
+
+      !   temporary arrays needed in the compute procedure
+      integer :: nCellsSolve
+      integer, dimension(:,:), pointer :: edgesOnCell
+      integer, dimension(:), pointer :: nEdgesOnCell
+      integer :: iCell,iEdge, i
+      real (kind=RKIND), dimension(:,:), pointer :: u
+      real (kind=RKIND), dimension(:,:), pointer :: uReconstructX, uReconstructY, uReconstructZ
+
+      real (kind=RKIND), dimension(:,:,:), pointer :: coeffs_reconstruct
+
+      ! stored arrays used during compute procedure
+      coeffs_reconstruct =&gt; grid % coeffs_reconstruct % array
+
+      ! temporary variables
+      edgesOnCell =&gt; grid % edgesOnCell % array
+      nEdgesOnCell=&gt; grid % nEdgesOnCell % array
+      nCellsSolve = grid % nCellsSolve
+      u =&gt; s % u % array
+      uReconstructX =&gt; s % uReconstructX % array
+      uReconstructY =&gt; s % uReconstructY % array
+      uReconstructZ =&gt; s % uReconstructZ % array
+
+      ! init the intent(out)
+      uReconstructX = 0.0
+      uReconstructY = 0.0
+      uReconstructZ = 0.0
+
+      ! loop over cell centers
+      do iCell=1,nCellsSolve
+        ! a more efficient reconstruction where rbf_values*matrix_reconstruct has been precomputed
+        ! in coeffs_reconstruct
+        do i=1,nEdgesOnCell(iCell)
+          iEdge = edgesOnCell(i,iCell)
+          uReconstructX(:,iCell) = uReconstructX(:,iCell) &amp;
+            + coeffs_reconstruct(1,i,iCell) * u(:,iEdge)
+          uReconstructY(:,iCell) = uReconstructY(:,iCell) &amp;
+            + coeffs_reconstruct(2,i,iCell) * u(:,iEdge)
+          uReconstructZ(:,iCell) = uReconstructZ(:,iCell) &amp;
+            + coeffs_reconstruct(3,i,iCell) * u(:,iEdge)
+        enddo
+      enddo   ! iCell
+
+   end subroutine reconstruct
+
+   subroutine get_distance(x1,x2,xout)
+     implicit none
+     real(kind=RKIND), intent(in) :: x1(3), x2(3)
+     real(kind=RKIND), intent(out) :: xout
+     xout = sqrt( (x1(1)-x2(1))**2 + (x1(2)-x2(2))**2 + (x1(3)-x2(3))**2 )
+   end subroutine get_distance
+ 
+   subroutine get_dotproduct(x1,x2,xout)
+     implicit none
+     real(kind=RKIND), intent(in) :: x1(3), x2(3)
+     real(kind=RKIND), intent(out) :: xout
+     xout = x1(1)*x2(1) + x1(2)*x2(2) + x1(3)*x2(3)
+   end subroutine get_dotproduct
+ 
+ 
+   subroutine unit_vector_in_R3(xin)
+     implicit none
+     real (kind=RKIND), intent(inout) :: xin(3)
+     real (kind=RKIND) :: mag
+     mag = sqrt(xin(1)**2+xin(2)**2+xin(3)**2)
+     xin(:) = xin(:) / mag
+   end subroutine unit_vector_in_R3
+
+
+   subroutine evaluate_rbf(xin,xout)
+     real(kind=RKIND), intent(in) ::  xin
+     real(kind=RKIND), intent(out) :: xout
+
+     ! Gaussian
+     ! xout = exp(-r**2)
+
+     ! multiquadrics
+       xout = 1.0 / sqrt(1.0**2 + xin**2)
+
+     ! other
+     ! xout = 1.0 / (1.0**2 + r**2)
+
+   end subroutine evaluate_rbf
+
+!======================================================================
+! BEGINNING OF CROSS_PRODUCT_IN_R3
+!======================================================================
+        subroutine cross_product_in_R3(p_1,p_2,p_out)
+
+!-----------------------------------------------------------------------
+! PURPOSE: compute p_1 cross p_2 and place in p_out
+!-----------------------------------------------------------------------
+
+!-----------------------------------------------------------------------
+! intent(in)
+!-----------------------------------------------------------------------
+        real (kind=RKIND), intent(in) ::                            &amp;
+                        p_1 (3),                                      &amp;
+                        p_2 (3)
+
+!-----------------------------------------------------------------------
+! intent(out)
+!-----------------------------------------------------------------------
+        real (kind=RKIND), intent(out) ::                           &amp;
+                        p_out (3)
+
+        p_out(1) = p_1(2)*p_2(3)-p_1(3)*p_2(2)
+        p_out(2) = p_1(3)*p_2(1)-p_1(1)*p_2(3)
+        p_out(3) = p_1(1)*p_2(2)-p_1(2)*p_2(1)
+
+        end subroutine cross_product_in_R3
+!======================================================================
+! END OF CROSS_PRODUCT_IN_R3
+!======================================================================
+
+! Updated 10/24/2001.
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!   Program 4.4   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!                                                                       !
+! Please Note:                                                          !
+!                                                                       !
+! (1) This computer program is written by Tao Pang in conjunction with  !
+!     his book, &quot;An Introduction to Computational Physics,&quot; published   !
+!     by Cambridge University Press in 1997.                            !
+!                                                                       !
+! (2) No warranties, express or implied, are made for this program.     !
+!                                                                       !
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+SUBROUTINE MIGS (A,N,X,INDX)
+!
+! Subroutine to invert matrix A(N,N) with the inverse stored
+! in X(N,N) in the output.  Copyright (c) Tao Pang 2001.
+!
+  IMPLICIT NONE
+  INTEGER, INTENT (IN) :: N
+  INTEGER :: I,J,K
+  INTEGER, INTENT (OUT), DIMENSION (N) :: INDX
+  REAL (kind=RKIND), INTENT (INOUT), DIMENSION (N,N):: A
+  REAL (kind=RKIND), INTENT (OUT), DIMENSION (N,N):: X
+  REAL (kind=RKIND), DIMENSION (N,N) :: B
+!
+  DO I = 1, N
+    DO J = 1, N
+      B(I,J) = 0.0
+    END DO
+  END DO
+  DO I = 1, N
+    B(I,I) = 1.0
+  END DO
+!
+  CALL ELGS (A,N,INDX)
+!
+  DO I = 1, N-1
+    DO J = I+1, N
+      DO K = 1, N
+        B(INDX(J),K) = B(INDX(J),K)-A(INDX(J),I)*B(INDX(I),K)
+      END DO
+    END DO
+  END DO
+!
+  DO I = 1, N
+    X(N,I) = B(INDX(N),I)/A(INDX(N),N)
+    DO J = N-1, 1, -1
+      X(J,I) = B(INDX(J),I)
+      DO K = J+1, N
+        X(J,I) = X(J,I)-A(INDX(J),K)*X(K,I)
+      END DO
+      X(J,I) =  X(J,I)/A(INDX(J),J)
+    END DO
+  END DO
+END SUBROUTINE MIGS
+
+
+SUBROUTINE ELGS (A,N,INDX)
+!
+! Subroutine to perform the partial-pivoting Gaussian elimination.
+! A(N,N) is the original matrix in the input and transformed matrix
+! plus the pivoting element ratios below the diagonal in the output.
+! INDX(N) records the pivoting order.  Copyright (c) Tao Pang 2001.
+!
+  IMPLICIT NONE
+  INTEGER, INTENT (IN) :: N
+  INTEGER :: I,J,K,ITMP
+  INTEGER, INTENT (OUT), DIMENSION (N) :: INDX
+  REAL (kind=RKIND) :: C1,PI,PI1,PJ
+  REAL (kind=RKIND), INTENT (INOUT), DIMENSION (N,N) :: A
+  REAL (kind=RKIND), DIMENSION (N) :: C
+!
+! Initialize the index
+!
+  DO I = 1, N
+    INDX(I) = I
+  END DO
+!
+! Find the rescaling factors, one from each row
+!
+  DO I = 1, N
+    C1= 0.0
+    DO J = 1, N
+      !C1 = AMAX1(C1,ABS(A(I,J)))
+      C1 = MAX(C1,ABS(A(I,J)))
+    END DO
+    C(I) = C1
+  END DO
+!
+! Search the pivoting (largest) element from each column
+!
+  DO J = 1, N-1
+    PI1 = 0.0
+    DO I = J, N
+      PI = ABS(A(INDX(I),J))/C(INDX(I))
+      IF (PI.GT.PI1) THEN
+        PI1 = PI
+        K   = I
+      ENDIF
+    END DO
+!
+! Interchange the rows via INDX(N) to record pivoting order
+!
+    ITMP    = INDX(J)
+    INDX(J) = INDX(K)
+    INDX(K) = ITMP
+    DO I = J+1, N
+      PJ  = A(INDX(I),J)/A(INDX(J),J)
+!
+! Record pivoting ratios below the diagonal
+!
+      A(INDX(I),J) = PJ
+!
+! Modify other elements accordingly
+!
+      DO K = J+1, N
+        A(INDX(I),K) = A(INDX(I),K)-PJ*A(INDX(J),K)
+      END DO
+    END DO
+  END DO
+!
+END SUBROUTINE ELGS
+
+end module vector_reconstruction

</font>
</pre>