<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 @@
        "CC = mpicc" \
        "SFC = gfortran" \
        "SCC = gcc" \
-        "FFLAGS = -O3 -m64 -ffree-line-length-none" \
+        "FFLAGS = -O3 -m64 -ffree-line-length-none -fdefault-real-8" \
        "CFLAGS = -O3 -m64" \
        "LDFLAGS = -O3 -m64" \
        "CORE = $(CORE)" \
        "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) $(EXPAND_LEVELS) -D_MPI -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+g95:
+        ( make all \
+        "FC = mpif90" \
+        "CC = mpicc" \
+        "SFC = g95" \
+        "SCC = gcc" \
+        "FFLAGS = -O3 -ffree-line-length-huge -r8" \
+        "CFLAGS = -O3" \
+        "LDFLAGS = -O3" \
+        "CORE = $(CORE)" \
+        "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) $(EXPAND_LEVELS) -D_MPI -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+g95-serial:
+        ( make all \
+        "FC = g95" \
+        "CC = gcc" \
+        "SFC = g95" \
+        "SCC = gcc" \
+        "FFLAGS = -O3 -ffree-line-length-huge -r8" \
+        "CFLAGS = -O3" \
+        "LDFLAGS = -O3" \
+        "CORE = $(CORE)" \
+        "CPPFLAGS = -DRKIND=8 $(MODEL_FORMULATION) $(EXPAND_LEVELS) -DUNDERSCORE $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+
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
/
&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 @@
&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
/
&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'
/
&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="$(SCC)" )
@@ -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) $< > $*.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 => 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 => state % h % array
scalars => 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), &
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)) &
+ *(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 -> 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)) &
+ ptmp = 0.5*(pressure(k,iCell)+pressure(k+1,iCell))
+ if (znuc(k) >= 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)) &
*sqrt(cos(znuv(k)))* &
((-2.*sin(phi)**6 &
*(cos(phi)**2+1./3.)+10./63.) &
@@ -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), &
+ 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 < 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) >= 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)) &
+ *sqrt(cos(znuv(k)))* &
+ ((-2.*sin(phi)**6 &
+ *(cos(phi)**2+1./3.)+10./63.) &
+ *2.*u0*cos(znuv(k))**1.5 &
+ +(1.6*cos(phi)**3 &
+ *(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)* &
+ (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))* &
+ (0.5*(pressure(k,iCell)+pressure(k+1,iCell))/p0)**cvpm
+
+ if (temperature(k,iCell) > 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), &
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, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
call dmpar_exch_halo_field2dReal(domain % dminfo, block % mesh % cqu % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
call dmpar_exch_halo_field2dReal(domain % dminfo, block % time_levs(2) % state % h % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
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) $< > $*.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 => grid % matrix_reconstruct % array
- normal => grid % normal % array
- rbf_value => grid % rbf_value % array
- coeffs_reconstruct => grid % coeffs_reconstruct % array
-
- !========================================================
- ! temporary variables needed for init procedure
- !========================================================
- xCell => grid % xCell % array
- yCell => grid % yCell % array
- zCell => grid % zCell % array
- cellsOnCell => grid % cellsOnCell % array
- cellsOnEdge => grid % cellsOnEdge % array
- edgesOnCell => grid % edgesOnCell % array
- nEdgesOnCell=> grid % nEdgesOnCell % array
- dcEdge => 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) &
- * 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 => grid % matrix_reconstruct % array
- normal => grid % normal % array
- rbf_value => grid % rbf_value % array
- coeffs_reconstruct => grid % coeffs_reconstruct % array
-
- ! temporary variables
- edgesOnCell => grid % edgesOnCell % array
- nEdgesOnCell=> grid % nEdgesOnCell % array
- nCells = grid % nCells
- nCellsSolve = grid % nCellsSolve
- nEdges = grid % nEdges
- nVertLevels = grid % nVertLevels
- xCell => grid % xCell % array
- yCell => grid % yCell % array
- zCell => grid % zCell % array
-
- u => s % u % array
- uReconstructX => s % uReconstructX % array
- uReconstructY => s % uReconstructY % array
- uReconstructZ => 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) :: &
- p_1 (3), &
- p_2 (3)
-
-!-----------------------------------------------------------------------
-! intent(out)
-!-----------------------------------------------------------------------
- real (kind=RKIND), intent(out) :: &
- 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) $< > $*.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 + &
- 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 => 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 => grid % matrix_reconstruct % array
- normal => grid % normal % array
- rbf_value => grid % rbf_value % array
-
- !========================================================
- ! temporary variables needed for init procedure
- !========================================================
- xCell => grid % xCell % array
- yCell => grid % yCell % array
- zCell => grid % zCell % array
- cellsOnCell => grid % cellsOnCell % array
- cellsOnEdge => grid % cellsOnEdge % array
- edgesOnCell => grid % edgesOnCell % array
- nEdgesOnCell=> grid % nEdgesOnCell % array
- dcEdge => 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 => grid % matrix_reconstruct % array
- normal => grid % normal % array
- rbf_value => grid % rbf_value % array
-
- ! temporary variables
- edgesOnCell => grid % edgesOnCell % array
- nEdgesOnCell=> grid % nEdgesOnCell % array
- nCells = grid % nCells
- nCellsSolve = grid % nCellsSolve
- nEdges = grid % nEdges
- nVertLevels = grid % nVertLevels
- u => s % u % array
- uReconstructX => s % uReconstructX % array
- uReconstructY => s % uReconstructY % array
- uReconstructZ => 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 "netcdf_read_ids.inc"
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) $< > $*.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) $< > $*.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 => grid % coeffs_reconstruct % array
-
- !========================================================
- ! temporary variables needed for init procedure
- !========================================================
- xCell => grid % xCell % array
- yCell => grid % yCell % array
- zCell => grid % zCell % array
- cellsOnCell => grid % cellsOnCell % array
- cellsOnEdge => grid % cellsOnEdge % array
- edgesOnCell => grid % edgesOnCell % array
- nEdgesOnCell=> grid % nEdgesOnCell % array
- dcEdge => 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) &
- + rbfValue * normals(:,i) * invMatrix(i,j)
- enddo
- enddo
- ! polynomial parts
- do j=1,npoints
- coeffs_reconstruct(:,j,iCell) = coeffs_reconstruct(:,j,iCell) &
- + invMatrix(npoints+1,j)*xHatPlane
- coeffs_reconstruct(:,j,iCell) = coeffs_reconstruct(:,j,iCell) &
- + 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 => grid % coeffs_reconstruct % array
-
- ! temporary variables
- edgesOnCell => grid % edgesOnCell % array
- nEdgesOnCell=> grid % nEdgesOnCell % array
- nCellsSolve = grid % nCellsSolve
- u => s % u % array
- uReconstructX => s % uReconstructX % array
- uReconstructY => s % uReconstructY % array
- uReconstructZ => 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) &
- + coeffs_reconstruct(1,i,iCell) * u(:,iEdge)
- uReconstructY(:,iCell) = uReconstructY(:,iCell) &
- + coeffs_reconstruct(2,i,iCell) * u(:,iEdge)
- uReconstructZ(:,iCell) = uReconstructZ(:,iCell) &
- + 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) :: &
- p_1 (3), &
- p_2 (3)
-
-!-----------------------------------------------------------------------
-! intent(out)
-!-----------------------------------------------------------------------
- real (kind=RKIND), intent(out) :: &
- 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, "An Introduction to Computational Physics," 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 => grid % coeffs_reconstruct % array
+
+ !========================================================
+ ! temporary variables needed for init procedure
+ !========================================================
+ xCell => grid % xCell % array
+ yCell => grid % yCell % array
+ zCell => grid % zCell % array
+ cellsOnCell => grid % cellsOnCell % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ edgesOnCell => grid % edgesOnCell % array
+ nEdgesOnCell=> grid % nEdgesOnCell % array
+ dcEdge => 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) &
+ + rbfValue * normals(:,i) * invMatrix(i,j)
+ enddo
+ enddo
+ ! polynomial parts
+ do j=1,npoints
+ coeffs_reconstruct(:,j,iCell) = coeffs_reconstruct(:,j,iCell) &
+ + invMatrix(npoints+1,j)*xHatPlane
+ coeffs_reconstruct(:,j,iCell) = coeffs_reconstruct(:,j,iCell) &
+ + 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 => grid % coeffs_reconstruct % array
+
+ ! temporary variables
+ edgesOnCell => grid % edgesOnCell % array
+ nEdgesOnCell=> grid % nEdgesOnCell % array
+ nCellsSolve = grid % nCellsSolve
+ u => s % u % array
+ uReconstructX => s % uReconstructX % array
+ uReconstructY => s % uReconstructY % array
+ uReconstructZ => 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) &
+ + coeffs_reconstruct(1,i,iCell) * u(:,iEdge)
+ uReconstructY(:,iCell) = uReconstructY(:,iCell) &
+ + coeffs_reconstruct(2,i,iCell) * u(:,iEdge)
+ uReconstructZ(:,iCell) = uReconstructZ(:,iCell) &
+ + 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) :: &
+ p_1 (3), &
+ p_2 (3)
+
+!-----------------------------------------------------------------------
+! intent(out)
+!-----------------------------------------------------------------------
+ real (kind=RKIND), intent(out) :: &
+ 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, "An Introduction to Computational Physics," 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>